GROMACS version: 2021
GROMACS modification: Yes
Hello,
I want to calculate a Coulomb force contribution from my simulations for every atom in simulation to use it later in visualization. What is the easiest way to extract the force components?
I started to do it myself by applying Coulomb law in coordinate processing and then comparing the electrostatic forces with total forces which I calculated from an xvg dump file. As a result, I got some values that I decided to double-check. After all, simplifying the system step-by-step, I ended up with a simple synthetic example system with just one H and one Cl ions with turned off LJ interactions. I hoped that now I would get a 100% match so I can make sure that my calculations are correct. But it appears to be that they still mismatch quite significantly:
H ion: total force (by gromacs): 9.4910180592178843 kJ/mol/nm vs coulomb (my utility) 11.3263893059850318 kJ/mol/nm
Cl ion: total force (by gromacs): 9.4894741077024918 kJ/mol/nm vs coulomb (my utility) 11.3263893059850318 kJ/mol/nm
What is strange to me is that:
- Why are the force values for H ion and Cl ion are still different (9.491 vs 9.489 kJ/mol/nm)?
- Why is there such a discrepancy between my and gromacs values (11.326 vs 9.491 kJ/mol/nm for NPT and ~11 vs 13 kJ/mol/nm for NVT)? I understand that I apply a slightly different approach (PME vs Coulomb law), but why is there such a difference? For more complex examples the difference can be up to 50%. Are there different forces there except electrostatics? Does the thermostat/barostat/anything else modify the forces?
Here is the code on how I calculate the interactions by Coulomb law (C++):
double k_e_si = 8.9875517923e9; // Coulomb constant; N * m^2 / C^2
double N_a = 6.0221409e23;
double e_to_coulomb = 1.602176565e-19; // e to C
// force in SI -> gromacs: N / mol = J / (mol * m) = 10^-3 kJ / (mol * 10^9 nm) = 10^-12 kJ / (mol * nm)
const double N_per_m_to_kJ_per_nm = 1e-12;
for (int i = 0; i < in->N_atoms; ++i) {
f_coul[i][0] = 0;
f_coul[i][1] = 0;
f_coul[i][2] = 0;
for (int j = 0; j < in->N_atoms; ++j) {
if (i != j) {
// position delta vector components between particles i and j
double delta_x = pbc(coords[i][0] - coords[j][0], in->box_x); // Angstrom
double delta_y = pbc(coords[i][1] - coords[j][1], in->box_y); // Angstrom
double delta_z = pbc(coords[i][2] - coords[j][2], in->box_z); // Angstrom
// square distance between particles i and j in Angstrom and meters
double delta_r_sq_A2 = delta_x * delta_x + delta_y * delta_y + delta_z * delta_z; // Angstrom^2
double delta_r_sq_m2 = 1e-20 * delta_r_sq_A2; // m^2
// force absolute value between particles i and j
double force_ij_abs = N_per_m_to_kJ_per_nm * N_a * e_to_coulomb * e_to_coulomb
* k_e_si * charges[i] * charges[j] / delta_r_sq_m2;
// force vector components from given delta r vector and the force modulus
double delta_r = sqrt(delta_r_sq_A2);
double fx = force_ij_abs * delta_x / delta_r;
double fy = force_ij_abs * delta_y / delta_r;
double fz = force_ij_abs * delta_z / delta_r;
// add to the total vector
f_coul[i][0] += fx;
f_coul[i][1] += fy;
f_coul[i][2] += fz;
}
}
}
The param file (a slightly modified file from a gromacs tutorial on lysozyme in water simulation; I also tried NVT with a similar result):
title = OPLS Lysozyme NPT equilibration
; Run parameters
integrator = md ; leap-frog integrator
nsteps = 5000 ; 2 * 500000000 = 1000000 ps = 1000 ns = 1us
dt = 0.001 ; 2 fs
; Output control
nstxout = 1000
nstvout = 1000
nstfout = 1000
nstenergy = 1000
nstlog = 1000
nstxout-compressed = 1000 ; save compressed coordinates
compressed-x-grps = System ; save the whole system
; Bond parameters
continuation = yes ; Restarting after NPT
constraint_algorithm = lincs ; holonomic constraints
constraints = h-bonds ; bonds involving H are constrained
lincs_iter = 1 ; accuracy of LINCS
lincs_order = 4 ; also related to accuracy
; Neighborsearching
cutoff-scheme = Verlet ; Buffered neighbor searching
ns_type = grid ; search neighboring grid cells
nstlist = 10 ; 20 fs, largely irrelevant with Verlet scheme
rcoulomb = 4.99 ; short-range electrostatic cutoff (in nm)
rvdw = 4.99 ; short-range van der Waals cutoff (in nm)
; Electrostatics
coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics
pme_order = 4 ; cubic interpolation
fourierspacing = 0.16 ; grid spacing for FFT
;coulombtype = Cut-off
; Temperature coupling is on
tcoupl = V-rescale ; modified Berendsen thermostat
tc-grps = Water ; two coupling groups - more accurate
tau_t = 0.1 ; time constant, in ps
ref_t = 300 ; reference temperature, one for each group, in K
; Pressure coupling is on
pcoupl = Parrinello-Rahman ; Pressure coupling on in NPT
pcoupltype = isotropic ; uniform scaling of box vectors
tau_p = 2.0 ; time constant, in ps
ref_p = 1.0 ; reference pressure, in bar
compressibility = 4.5e-5 ; isothermal compressibility of water, bar^-1
; Periodic boundary conditions
pbc = xyz ; 3-D PBC
; Dispersion correction
DispCorr = EnerPres ; account for cut-off vdW scheme
; Velocity generation
gen_vel = no ; Velocity generation is off
Here is how I get the xvg file (analyze.sh file):
OUTPREFIX="traj_comp"
echo 0 | gmx trjconv -s "${OUTPREFIX}.gro" -f "${OUTPREFIX}.xtc" -o "${OUTPREFIX}_movie.gro" -dt 1
echo 0 | gmx traj -f "${OUTPREFIX}.trr" -s "${OUTPREFIX}.tpr" -of "${OUTPREFIX}_force_dump.xvg" -dt 1
(then I call my utility)
The rest of the files, including the xyz files (can be opened in text editor and Ovito), are also attached.
topol.top (749 Bytes)
npt.mdp (2.5 KB)
The whole archive with xyz files