[gmx-users] Setting attractive part (C6 parameter) zero in non-bonded parameters

[Suman Chakrabarty]

Hello,

Following the Section 5.3.2 of the manual, I have tried to set the C6
term of the non-bonded (Lennard-Jones) interaction to zero as follows:

"When sigma and epsilon need to be supplied (rules 2 and 3), it
would seem it is impossible to have a non-zero C12 combined with a
zero C6 parameter. However, providing a negative sigma will do
exactly that, such that C6 is set to zero and C12 is calculated
normally. This situation represents a special case in reading the
value of sigma, and nothing more."

Unfortunately when I try this the potential energy of the system becomes "nan":
           Step           Time         Lambda
              0        0.00000        0.00000

   Energies (kJ/mol)
        LJ (SR)   Coulomb (SR)   Coul. recip.      Potential    Kinetic En.
            nan   -1.76446e+05   -1.26632e+03            nan    2.30786e+04
   Total Energy  Conserved En.    Temperature Pressure (bar)
            nan            nan    2.25008e+02            nan

With the positive sign for sigma in the topology file, the same
structure gives the following energetics:
   Energies (kJ/mol)
        LJ (SR)   Coulomb (SR)   Coul. recip.      Potential    Kinetic En.
    4.56607e+04   -1.76446e+05   -1.26632e+03   -1.32051e+05    2.30777e+04
   Total Energy  Conserved En.    Temperature Pressure (bar)
   -1.08974e+05   -1.08974e+05    2.24999e+02    3.52545e+04

=========================================================
Here I am trying to look at a single solute atom (I change the sign of
the sigma value for this atom only) in a box of TIP4P water.
Strangely, I encounter this "nan" problem only for comb-rule=3.
Whereas, for comb-rule=2, energetics are reasonable, but it appears
that sigma_ij between the solute and solvent becomes very small (as
evident for rdf of water oxygen atoms around this solute atom). This
would imply that the negative value of sigma_i for the solute is being
retained as it is, so in 0.5*(sigma_i + sigma_j) they cancel each
other to large extent.

It seems there is a bug in the treatment of the negative sigma values
as described in the manual. Please let me know if you think I am
making a mistake somewhere. Thanks.

Best regards,
Suman.