I'm entering this discussion late so I may be discussing issues that
have already been addressed.
As I understand it, Federico, you began by describing a model for data
in which two factors have a fixed set of levels and one factor has an
extensible, or "random", set of levels and you wanted to fit a model
that you described as
y ~ effect1 * effect2 * effect3
The problem is that this specification is not complete. An
interaction of factors with fixed levels and a factor with random
levels can mean, in the lmer specification,
lmer(y ~ effect1 * effect2 + (1| effect3) + (1|effect1:effect2:effect3), ...)
or
lmer(y ~ effect1 * effect2 + (effect1*effect2 | effect3), ...)
or other variations. When you specify a random effect or an random
interaction term you must, either explicitly or implicitly, specify
the form of the variance-covariance matrix associated with those
random effects.
The "advantage" that other software may provide for you is that it
chooses the model for you but that, of course, means that you only
have the one choice.
If you can describe how many variance components you think should be
estimated in your model and what they would represent then I think it
will be easier to describe how to fit the model.