[gmx-users] Strange energy values with GMX 5.1, PME and Verlet

[Mark Abraham]

Hi,

On Thu, Dec 28, 2017 at 6:40 AM Diana Lousa <dlousa at itqb.unl.pt> wrote:

> Dear all,
>
> For testing purposes, I performed a short simulation (10 ps) of a CL- and a
> Na+ ion in vacuum with gmx 5.1.4 and compared the Group and Verlet cutoff
> schemes. I ran gmx energy to monitor the interaction energy between the
> ions and obtained the following results:
>
> 1) Verlet
> Energy                      Average   Err.Est.       RMSD  Tot-Drift
>
> -------------------------------------------------------------------------------
> Pressure                   0.128465      0.044    15.8209   -0.27251  (bar)
> Coul-SR:NA-NA              -244.832          0       -nan -4.68375e-14
> (kJ/mol)
> LJ-SR:NA-NA                       0          0          0          0
> (kJ/mol)
> Coul-SR:NA-CL              -153.193        7.4    32.7214   -43.0706
> (kJ/mol)
> LJ-SR:NA-CL                 42.5252        1.7    20.9531    9.65428
> (kJ/mol)
> Coul-SR:CL-CL              -244.832          0       -nan -4.68375e-14
> (kJ/mol)
> LJ-SR:CL-CL                       0          0          0          0
> (kJ/mol)
>
> 2) Group
> Energy                      Average   Err.Est.       RMSD  Tot-Drift
>
> -------------------------------------------------------------------------------
> Pressure                  0.0159313      0.039    17.1949   0.303503  (bar)
> Coul-SR:NA-NA                     0          0          0          0
> (kJ/mol)
> LJ-SR:NA-NA                       0          0          0          0
> (kJ/mol)
> Coul-SR:NA-CL              -152.361          8    33.2265   -46.2306
> (kJ/mol)
> LJ-SR:NA-CL                 42.2222        1.9    22.5877    11.0263
> (kJ/mol)
> Coul-SR:CL-CL                     0          0          0          0
> (kJ/mol)
> LJ-SR:CL-CL                       0          0          0          0
> (kJ/mol)
>
> As you can see, the Coul-SR:NA-NA and Coul-SR:CL-CL have the value of
> -244.832 when using Verlet, which is very strange.


PME for Verlet and group is implemented differently. The various kinds of
atom pairs whose non-bonded interactions must be excluded can be computed
in different ways. The approach in the Verlet scheme is more efficient, but
does not reproduce that of the group scheme, But anyway since such
decomposition is physically meaningless, it doesn't matter.


> This does not happen
> when using Group cutoff or when using Verlet in GMX 4.5 (the value is 0 in
> these cases).
>

4.5 did not have such a scheme.


> GMX 2016 gives similar results to GMX 5.1.4:
>
> 3) Verlet GMX 2016:
> Energy                      Average   Err.Est.       RMSD  Tot-Drift
>
> -------------------------------------------------------------------------------
> Pressure                 -0.0106255       0.22    19.2119   0.756146  (bar)
> Coul-SR:NA-NA              -244.832          0       -nan -4.68375e-14
> (kJ/mol)
> LJ-SR:NA-NA                       0          0          0          0
> (kJ/mol)
> Coul-SR:NA-CL              -150.921        7.9     35.159   -45.7216
> (kJ/mol)
> LJ-SR:NA-CL                 42.0629          2    25.1668    11.2282
> (kJ/mol)
> Coul-SR:CL-CL              -244.832          0       -nan -4.68375e-14
> (kJ/mol)
> LJ-SR:CL-CL                       0          0          0          0
> (kJ/mol)
>
> I am not using GPUs, since I know you cannot monitor energy groups when
> using GPUs.
> Does anyone have any idea why this happens and what consequences does it
> have?
>

None.

Mark


>
> Here's the mdp:
>
> ****************************
> define              =
> integrator          =  md
> tinit               =  0.0
> dt                  =  0.002
> nsteps              =  5000
> nstcomm             =  1
> nstxtcout           =  500
> xtc-precision       =  1000
> nstxout             =  0
> nstvout             =  0
> nstfout             =  0
> nstlog              =  10
> nstenergy           =  1
>
> nstlist             =  20
> ns_type             =  grid
>
> ;PME
> rlist               = 1.0
> coulombtype         = PME
> rcoulomb            = 1.0
> fourierspacing      = 0.16
> pme_order           = 4
> vdwtype             = cut-off
> rvdw                = 1.0
> cutoff-scheme         = Verlet ; or Group
>
> ; Energy monitoring
> energygrps          =  NA CL
>
> Tcoupl              =  v-rescale
> tc-grps             =  Ion
> tau_t               =  0.10
> ref_t               =  300.0
>
> ; Isotropic pressure coupling is now of
> ; Pcoupl              =  berendsen
> Pcoupl              =  No
>
>
> ; lincs_
> constraint_algorithm=lincs
> lincs_order         =  8
> constraints         = none
> lincs-warnangle     =  90
>
>
> ​Thanks in advance
> Best regards​
>
> --
> Diana Lousa
> Post-doc
> Protein Modeling Laboratory
> ITQB/UNL
> Oeiras, Portugal
> --
> Gromacs Users mailing list
>
> * Please search the archive at
> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
> posting!
>
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
> * For (un)subscribe requests visit
> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
> send a mail to gmx-users-request at gromacs.org.