[gmx-users] twin range and switched cutoffs energy differences

[David Osguthorpe]

Hi,

I am having a problem understanding the GROMACS implementation of cutoffs
and switching functions - at least Im getting results I cant make sense of.

I have been running dynamics of a protein-water system rhombic dodecahedron
periodic boundary system with the twin range cutoff scheme.
I have been attempting to minimize coordinate snapshots from this simulation
and to get to a reasonable maximum force was using the switching function version
of cutoffs (in addition to flexible waters and no constraints).

What Im seeing is a large difference in the coulomb energies which I do not
understand where it is coming from.
To pin this down I am running the following test cases from a minimized system and
computing the energy using two inputs.
First, twin range cutoffs (I got the rlist = 0.01 from some gmxtest runs)

rlist               =  0.01
coulombtype         =  cut-off
rcoulomb            =  1.2
vdwtype             =  cut-off
rvdw                =  0.9

Second, switched cutoffs

coulombtype         =  switch
rcoulomb_switch     =  1.199
rcoulomb            =  1.201
vdwtype             =  switch
rvdw_switch         =  0.899
rvdw                =  0.901

Twin cutoffs log energy:

   Energies (kJ/mol)
           Bond          Angle    Proper Dih. Ryckaert-Bell.          LJ-14
    7.03530e+04    4.33556e+04    4.98766e+02    4.81356e+03    5.41441e+03
     Coulomb-14        LJ (SR)        LJ (LR)   Coulomb (SR)   Coulomb (LR)
    2.27345e+04    0.00000e+00    2.44898e+04    0.00000e+00   -3.16338e+05
      Potential    Kinetic En.   Total Energy    Temperature Pressure (bar)
   -1.44678e+05    0.00000e+00    0.00000e+00    0.00000e+00    0.00000e+00

Switched cutoff log energy:
      
   Energies (kJ/mol) 
           Bond          Angle    Proper Dih. Ryckaert-Bell.          LJ-14
    7.03530e+04    4.33556e+04    4.98766e+02    4.81356e+03    5.41441e+03
     Coulomb-14        LJ (SR)        LJ (LR)   Coulomb (SR)   Coulomb (LR)
    2.27345e+04    0.00000e+00    2.46563e+04    0.00000e+00   -1.06533e+06
      Potential    Kinetic En.   Total Energy    Temperature Pressure (bar)
   -8.93502e+05    0.00000e+00    0.00000e+00    0.00000e+00    0.00000e+00


In the initial computation of the energies there is a large difference between
the Coulomb (LR) component - Ive tried using switching ranges from -/+ 0.15 nm but nothing
seems to reduce the difference in the coulomb LR from a near order of magnitude
- yet the above switch functions should be close to the cutoff function so I dont
expect such a large difference.

Thanks

David