[gmx-users] non-zero Coul-SR interaction between metal ion in the system with only one metal (GMX 5.0.7)

[qiaobf]

Dear All,

I just got one problem with the energy printed by the g_energy. Could 
anyone help me with it? Thanks!

In the test system, there is ONLY ONE trivalent metal (Me) and 2240 
modified TIP3P waters. After one MD step, I print the energies using 
g_energy. The results are:

** GROMACS 5.0.7 ************

*Coul-SR:Me-Me              -521.008 *     -- 0          0  (kJ/mol)
LJ-SR:Me-Me                       0                 -- 0          0  
(kJ/mol)
Coul-SR:Me-SOL             -78.5025         --    30.1404 -60.2807  (kJ/mol)
LJ-SR:Me-SOL                195.268            --    104.218 -208.437  
(kJ/mol)
Coul-SR:SOL-SOL             1958.38         --     452.82 -905.64  (kJ/mol)
LJ-SR:SOL-SOL               91953.5           --    1966.23 -3932.45  
(kJ/mol)

*********************************

The weird thing is the Coul-SR: Me-Me is not zero, instead it is negative!

Once I change the atomic charge of the trivalent metal from +3 to 0, 
this Coul-SR: Me-Me becomes zero.

I also rerun the simulation using GROAMCS 4.5.5. The Coul-SR is zero 
(see below). But there exist some difference between the other energies 
between 5.0.7 and 4.5.5.

** GROMACS 4.5.5 ******************

Coul-SR:Me-Me                     0            -- 0          0  (kJ/mol)
LJ-SR:Me-Me                         0            -- 0          0  (kJ/mol)
Coul-SR:Me-SOL             -57.4723      --          0          0 (kJ/mol)
LJ-SR:Me-SOL                94.4865         --          0 0  (kJ/mol)
Coul-SR:SOL-SOL             1771.36      --          0          0 (kJ/mol)
LJ-SR:SOL-SOL               89843.4        --          0 0  (kJ/mol)

**************************


The topology and md.mdp files are listed below:

***** system.top *****************

#include "charmm36-nov2016.ff/forcefield.itp"

[ atomtypes ]
;type at.num  mass     charge ptype   sigma, nm  epsilon, kJ/mol
  Me      1     1.000   0.0     A      3.385415e-01   2.092000e-01

;------------------------------ metal-----------------------
[ moleculetype ]
; molname       nrexcl
  Metal            1

[ atoms ]
; id    at type       res.nr  residu.name     at.name  cg.nr  charge
1       Me            1       Me              Me       1       3.0

#include "charmm36-nov2016.ff/tip3p.itp"

[ system ]
; Name
metal ion in water

[ molecules ]
; Compound        #mols
  Metal             1
  SOL             2240
*****************

****** md.mdp******

; Run control
integrator               = md
tinit                    = 0
dt                       = 0.001
nsteps                   = 1    ; 10ps
nstcomm                  = 10

; Output control
nstenergy                = 50
energygrps               = ME SOL


; Neighborsearching and short-range nonbonded interactions
nstlist                  = 10
ns_type                  = grid
pbc                      = xyz
rlist                    = 1.2
; Electrostatics
coulombtype              = cutoff
rcoulomb                 = 1.2
; van der Waals
vdwtype                  = cutoff
rvdw-switch              = 1.0
rvdw                     = 1.2
; Apply long range dispersion corrections for Energy and Pressure
DispCorr                  = EnerPres
; Spacing for the PME/PPPM FFT grid
fourierspacing           = 0.12
; EWALD/PME/PPPM parameters
pme_order                = 4
ewald_rtol               = 1e-06
epsilon_surface          = 0

; Temperature coupling
; tcoupl is implicitly handled by the sd integrator
tc_grps                  = system
tau_t                    = 1.0
ref_t                    = 298.15

; Pressure coupling is off for NVT
Pcoupl                   = No

***************