GROMACS version: 5.1.4 and 18.6
GROMACS modification: No
Hi,
Does anyone know why short-range Coulomb energies for cation-cation and anion-anion interactions in an aqueous solution are negative?
My system is simply an OPLS_aa alkane with water and NaCl and when I used energy groups I got this result for the crossed ion interactions:
Energy Average Err.Est. RMSD Tot-Drift
Coul-SR:NA-NA -6705.58 17 73.3361 -59.4862 (kJ/mol)
LJ-SR:NA-NA -0.0627519 0.0064 0.0266583 0.021759 (kJ/mol)
Coul-14:NA-NA 0 0 0 0 (kJ/mol)
LJ-14:NA-NA 0 0 0 0 (kJ/mol)
Coul-SR:NA-CL -2413.54 98 504.244 480.765 (kJ/mol)
LJ-SR:NA-CL 272.935 11 66.7708 -54.1533 (kJ/mol)
Coul-14:NA-CL 0 0 0 0 (kJ/mol)
LJ-14:NA-CL 0 0 0 0 (kJ/mol)
Coul-SR:CL-CL -6756.83 10 49.5624 -33.212 (kJ/mol)
LJ-SR:CL-CL -1.01777 0.34 3.58415 -0.957724 (kJ/mol)
Coul-14:CL-CL 0 0 0 0 (kJ/mol)
LJ-14:CL-CL 0 0 0 0 (kJ/mol)
I got a similar result with AMBER99sb.
I tried PME, Ewald sum, and a simple cut-off and got qualitatively similar results, although a cut-off gives larger values (less negative) for the like ion interactions.
I cannot see how the real space Ewald sum or the direct Coulomb formula can give a negative energy for particles of like charge. Perhaps there is some simple explanation for this?
Any hints would be highly appreciated.
Thanks in advance.
Nuno
****** mdp file *****
title = Solute in Water NPT Equilibration
;define = -DPOSRES ; position restrain the solute
; Run parameters
integrator = md ; leap-frog
nsteps = 5000000
dt = 0.002 ; 2 fs
tinit = 0
energygrps = MOL SOL NA CL
; Output control
nstxout = 0
nstvout = 0
nstenergy = 100
nstlog = 1000
; Output frequency and precision for .xtc file
nstxout-compressed = 1000 ; save coordinates in .xtc file every 5 time-steps (10 fs)
compressed-x-precision = 10000
; Bond parameters
continuation = yes ; continue from NVT
constraint_algorithm = lincs
;constraints = all-bonds ; all bonds (even heavy atom-H bonds) constrained
constraints = h-bonds ;
lincs_iter = 2
lincs_order = 4
; Neighborsearching
cutoff-scheme = Verlet
ns_type = grid ; search neighboring grid cels
;rlist = 1.5 ; short-range neighborlist cutoff (in nm)
rcoulomb = 1.0 ; short-range electrostatic cutoff (in nm)
rvdw = 1.0 ; short-range van der Waals cutoff (in nm)
verlet-buffer-tolerance = 1.0e-06
; Electrostatics
coulombtype = PME
pme_order = 4
fourierspacing = 0.16
; Temperature coupling is on
tcoupl = V-rescale
tc-grps = System
tau_t = 0.1
ref_t = 298
; Pressure coupling is on
pcoupl = Parrinello-Rahman
pcoupltype = isotropic ; uniform scaling of x-y-z box vectors
tau_p = 2.0
ref_p = 1
compressibility = 4.5e-5 ; isothermal compressibility, bar^-1
refcoord_scaling = com
; Periodic boundary conditions
pbc = xyz
; Dispersion correction
DispCorr = EnerPres
; Velocity generation
gen_vel = no
;gen_temp = 298
;gen_seed = -1
; COM motion removal
; These options remove COM motion of the system
nstcomm = 10
comm-mode = Linear
comm-grps = System
**** end mdp file ****
Last part of the log file
<====== ############### ==>
<==== A V E R A G E S ====>
<== ############### ======>
Statistics over 5000001 steps using 50001 frames
Energies (kJ/mol)
Bond Angle Ryckaert-Bell. LJ-14 Coulomb-14
1.45222e+01 7.61212e+01 2.25249e+01 1.28000e+01 -5.62811e+00
LJ (SR) Disper. corr. Coulomb (SR) Coul. recip. Potential
1.22024e+04 -5.84756e+02 -9.29796e+04 4.41796e+02 -8.07998e+04
Kinetic En. Total Energy Temperature Pres. DC (bar) Pressure (bar)
9.56039e+03 -7.12394e+04 2.98005e+02 -2.51559e+02 7.50236e-01
Constr. rmsd
0.00000e+00
Box-X Box-Y Box-Z
3.38067e+00 3.38067e+00 3.38067e+00
Total Virial (kJ/mol)
3.18669e+03 1.03961e+00 6.27159e+00
1.04365e+00 3.18816e+03 1.19012e+00
6.27510e+00 1.19317e+00 3.18664e+03
Pressure (bar)
1.46362e+00 -9.36082e-01 -5.47850e+00
-9.39604e-01 -4.51542e-01 -8.74366e-01
-5.48149e+00 -8.76964e-01 1.23863e+00
Epot (kJ/mol) Coul-SR LJ-SR Coul-14 LJ-14
MOL-MOL 1.84332e+01 -6.48157e+00 -5.62811e+00 1.28000e+01
MOL-SOL 4.09928e+00 -9.82426e+01 0.00000e+00 0.00000e+00
MOL-NA 7.50592e-01 -6.79239e-02 0.00000e+00 0.00000e+00
MOL-CL -1.03322e+00 -1.06480e+00 0.00000e+00 0.00000e+00
SOL-SOL -6.30204e+04 9.97591e+03 0.00000e+00 0.00000e+00
SOL-NA -6.41298e+03 1.06244e+03 0.00000e+00 0.00000e+00
SOL-CL -7.69247e+03 9.98018e+02 0.00000e+00 0.00000e+00
NA-NA -6.70558e+03 -6.27519e-02 0.00000e+00 0.00000e+00
NA-CL -2.41354e+03 2.72935e+02 0.00000e+00 0.00000e+00
CL-CL -6.75683e+03 -1.01777e+00 0.00000e+00 0.00000e+00
M E G A - F L O P S A C C O U N T I N G
NB=Group-cutoff nonbonded kernels NxN=N-by-N cluster Verlet kernels
RF=Reaction-Field VdW=Van der Waals QSTab=quadratic-spline table
W3=SPC/TIP3p W4=TIP4p (single or pairs)
V&F=Potential and force V=Potential only F=Force only