Center of mass motion issue

[Rana-1]

GROMACS version: 2020.4
GROMACS modification: Yes/No
Here post your question

Dear users
I ran a 10 ns MD for a DNA molecule.
the .mdp file I used is attached. After the MD is complete, I used the commands below to prepare my .xtc and .gro files I need for trajectory analysis.

gmx trjconv -s “A1_A2_pr.tpr” -f “A1_A2_pr.xtc” -o “A1_A2_pr_PBC.xtc” -pbc mol -center
gmx trjconv -s “A1_A2_pr.tpr” -f “A1_A2_pr.xtc” -o “A1_A2_pr_PBC.gro” -pbc mol -center -dump 1

I noticed that the molecule moves. It translates and rotates. My understanding that rotation is expected but translation is not because the mdout.mdp says that center of mass motion is removed every 100 steps.

md_pr.mdp
title = Simulations on Visual Dynamics

; Run parameters
integrator = md ; leap-frog integrator
nsteps = 5000000 ; 2 * 5000000 = 10000 ps (10 ns)
dt = 0.002 ; 2 fs

; Output control
nstxout = 0 ; save coordinates every 10.0 ps
nstvout = 0 ; save velocities every 10.0 ps
nstenergy = 5000 ; save energies every 10.0 ps
nstlog = 5000 ; update log file every 10.0 ps
nstxout-compressed = 5000 ; save compressed coordinates every 10.0 ps
; nstxout-compressed replaces nstxtcout
compressed-x-grps = System ; replaces xtc-grps

; Bond parameters
continuation = yes ; Restarting after NPT
constraint_algorithm = lincs ; holonomic constraints
constraints = h-bonds ; all bonds (even heavy atom-H bonds) constrained
lincs_iter = 1 ; accuracy of LINCS
lincs_order = 4 ; also related to accuracy

; Nonbonded settings
cutoff-scheme = Verlet ; Buffered neighbor searching
ns_type = grid ; search neighboring grid cells
nstlist = 10 ; 20 fs, largely irrelevant with Verlet
rlist = 1.2
vdwtype = Cut-off
vdw-modifier = Force-switch
rvdw_switch = 1.0
rvdw = 1.2 ; short-range van der Waals cutoff (in nm)
rcoulomb = 1.2 ; short-range electrostatic cutoff (in nm)
DispCorr = no ; account for cut-off vdW scheme

; Electrostatics
coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics
pme_order = 4 ; cubic interpolation
fourierspacing = 0.16 ; grid spacing for FFT

; Temperature coupling is on
tcoupl = V-rescale ; modified Berendsen thermostat
tc-grps = DNA Water_and_ions ; two coupling groups - more accurate
tau_t = 0.1 0.1 ; time constant, in ps
ref_t = 300 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

; Velocity generation
gen_vel = no ; Velocity generation is off

mdout.mdp
; mode for center of mass motion removal
comm-mode = Linear
; number of steps for center of mass motion removal
nstcomm = 100
; group(s) for center of mass motion removal
comm-grps =

[obZehn]

Dear @Rana-1 ,

The two things are different. The COM motion removal should be a calculation of the velocity of the COM of the whole box, and then that velocity is subtracted to all the atoms in the box, that is, you are setting to zero the velocity of the COM of the box. If I remember correctly, this makes sure that the box doesn’t accumulate kinetic energy over time due to numerical integration errors etc, which back in the day in the first MD works gave origin to the so called flying ice cube effect, where the rescaling of velocities to thermostat the systems led to frozen atomic configurations and box huge linear momenta (a flying ice cube!).

The removal of the COM is, very roughly, something to prevent an artifact. However, the ability of the system to diffuse around the box is something else, which is related to the nature of what you simulate. Here you have DNA in water at 300 K and I would say there is no reason to expect it to not diffuse around.