GROMACS version: 2023.3
GROMACS modification: No
Dear all,
I am trying to reproduce using gromacs the procedure described in these two papers: https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.85.3970 and https://pubs.acs.org/doi/full/10.1021/jp027823q .
Basically, I have a protein (2b3p) in vacuum for which I minimized the energy up to Fmax~10**-8 , before performing normal mode analysis using mdrun + nmeig . Until there, nothing special, and I think I handled it correctly.
Then, I want to perturb one and only one normal mode on this zero-enegy, zero-T structure, injecting a precisely controlled amount of kinetic energy on the chosen eigenvector.
First, I had to rescale the eigenvectors to make them orthonormal again (for some reason gmx nmeig outputs rescaled eigenvectors, that are not orthogonals!)
Then I used MDanalysis to inject 4kJ/mol along the said eigenvector; I spare you the detail here (I can send the script if someone wants to check it), it is enough to say that the followig script
Na=6.02214129e23
amu=1.660538921e-27
GMX_UNITS = 1/Na/amu*1e-1 # conversion from kJ/mol to amu.Ang^2/ps^2
u = mda.Universe(TPR, ‘injected.trr’)
vel = u.trajectory.ts.velocities
masses = u.atoms.masses
masses_per_coord = np.repeat(masses, 3)
v_flat = vel.flatten()
KE = 0.5 * np.sum(masses_per_coord * v_flat**2)
print(“KE recomputed from injected.gro:”, KE / GMX_UNITS)
outputs “KE recomputed from injected.gro: 4.000000072775341” (so very close to my target!).
BUT! If I do a mdrun at this point, using
gmx_mpi grompp -f md_inj.mdp -c injected.gro -t injected.trr -p GFPvac.top -o GFPinj.tpr -maxwarn 2
gmx mdrun -deffnm test_inj -s GFPinj.tpr
with the mdp:
integrator = md
nsteps = 10
dt = 0.0005
nstxout = 1
nstvout = 1
nstfout = 1
nstenergy = 1
compressed-x-grps = System
nstlist = 1
cutoff-scheme = Verlet
ns_type = grid
verlet-buffer-tolerance = -1
rlist = 1.4
coulombtype = Reaction-Field-zero
epsilon-rf = 0
rcoulomb = 1.2
rcoulomb-switch = 1.0
vdwtype = Cut-off
vdw-modifier = Potential-switch
rvdw = 1.2
rvdw-switch = 1.0
constraints = none
gen_vel = no
tcoupl = no
pcoupl = no
I get, using gmx energy, at t=0, a kinetic energy of 6.967627kJ/mol instead of 4!
I tested with different injected energies: for E_inj_kJmol = 10 I get 12.962667 , for E_inj_kJmol = 100 I get 102.933517 , for E_inj_kJmol = 1000 I get 1002.840820 … So it seems that gromacs is adding a constant ~2.8 kJ/mol … At the precision I wish to work with, it’s a problem, an I would like to understand what is happening in the background here.
Thanks in advance for your help!