Asking my students to use the margins command after mlogit for computing marginal effects (using dydx), I've found that including variables into the mlogit model by using the ## operator would give you different results from not using the operator.

For example:
mlogit y x1 x2 c.x3##c.x3 /* model 1*/
and
mlogit y x1 x2 x3 x3_2 /* model 2 */
where x3_2 is generated by: gen x3_2 = x3*x3
After model 1 or model 2 run,
margins, dydx(x1) at(x2=(0 1)) atmeans
would give you different results depending on which model is run first.

After model 2, the sample mean values of x3 and x3_2 are explicitly used while after model 1, only x3 mean is specified while x3_2 is not set at its mean (presumably individual x3_2 values are used before averaging the final result?).
It appears that if x3 is entered as at() such as "... at(x3=(20(1)80)), then the squared computation is passed on correctly but not when it's involved in "atmeans"... It is hard to imagine if that's a reasonable way to compute.
Anyone knows any reasons for this that we failed to see? Thank you.