[gmx-users] question regarding sc-power in gromacs 3.3

[David Mobley]

Dear all,

Concerning the new soft-core free energy parameter sc-power in GROMACS
3.3: What exactly is it supposed to do? I'm concerned that perhaps it
works as one would expect for sc-power=1 or 2 but not for other
values, but perhaps I'm mistaken. In particular, consider the
following code snippet from gmxlib/nonbonded/nb_free_energy.c:

alfA  = alpha*(lam_power==2 ? lambda*lambda : lambda);
alfB  = alpha*(lam_power==2 ? L1*L1 : L1);
dalfA = alpha*lam_power/6.0*(lam_power==2 ? lambda : 1);
dalfB = alpha*lam_power/6.0*(lam_power==2 ? L1 : 1);

lam_power is the same as sc_power (as read from the input mdp file)
and is set in the function call. There is no other reference to
lam_power or sc_power in the function.

This alfA, B, ends up controlling the denominator of the potential --
that is, r=(alfA*sigma^6+r^6)^(1/6), I think. For GROMACS 3.2, alfA is
always alpha*lambda^2. Here, I assumed that sc_power is supposed to
control the exponent (i.e. alfA=alpha*lambda^sc_power). But based on
the code snippet above, it looks like the only ALLOWED values of
sc_power (lam_power) are 1 or 2, as there is no definition of alfA or
alfB which incorporates other powers.

Is this handled elsewhere in the code, and is it handled correctly?
I've tried values of sc_power that are not 1 or 2, and I get results
different from those with sc_power=1 and 2, but it's not clear from
this that these are handled correctly. It's clear from dalfA that
changing lam_power will affect the computed free energies (there IS a
lam_power depndence there), but why aren't these something like:

alfA=alpha*pow(lambda,lam_power);
dalfA=alpha*lam_power/6*pow(lambda,lam_power);
or some such?

Thanks,
David