GROMACS version: 2020.4
GROMACS modification: No
I am interested in studying interactions between a model lipid bilayer and a compound. In order to speed up the simulation I am planning to use Hydrogen Mass Repartitioning for the lipids. Should I also perform HMR for the interacting compound? Can a mixed simulation with HMR for lipids and no HMR for the compound be performed?
Regards,
RPS
The whole system needs to have hydrogen masses repartitioned.
This mass redistribution is performed to decrease the frequency / speed of movement of the fastest moving atoms/bonds in the system, which is the lightest atoms (hydrogen). So if you only redistribute the mass of one set of molecules and leave another set, then the latter will still have the same set of fast/high frequency movement and you will not be able to use a longer time step i.e. you will achieve nothing.
Thank you Dr. Lemkul and Dr_DBW for the clarification.
Regards,
Raman
Dear Dr. Lemkul, Dr_DBW and Gromacs community,
I tried building the HMR membrane using CHARMM-GUI with CHARMM36 FF. The hydrogen and carbon masses of the lipid and protein were increased/decreased as expected (hydrogen mass increased to 3.024 and carbon mass decreased depending on nummber of hydrogen attached to it) but the hydrogen mass remained at 1.008 for the TIP3P water. As pointed above, shouldn’t the mass of hydrogen and oxygen be also changed in water?
Regards,
RPS
If CHARMM-GUI does not produce the expected output, I suggest you contact their developers directly. But the mass of H and O in water is controlled via the use of define = -DHEAVY_H in the .mdp file in the CHARMM36 port. I do not know if the CHARMM-GUI files provide this option, but this is the way it’s normally handled in the CHARMM force field for GROMACS. Other species have to have their masses explicitly reset; water is controlled at the input file level rather than modifying the topology.
Hello everyone,
So, to run HMR in GROMACS there are two things we need to change i.e. define = -DHEAVY_H and dt = 0.004? Because this is giving me error in the production run. Is there any tutorial or manual for HMR?
Thanks for your help,
Regards,
Gaurav
Can you provide more detail? Which force field are you using? Some only require HEAVY_H when using constraints, others need it for the water topology to function properly. Help us to help you.
Hello Dr. Lemkul,
So, I am performing protein-DNA MD simulation.
Force field: DNA_AAMBER99SB_ILDN_STAR_AA_MUTATION_FORCEFIELD for both protein and DNA and TIP3P for water
Then I have included define = -DHEAVY_H and dt = 0.004 in the mdp file.
After 1ns of production run I am getting the following error.
Fatal error:
Step 10100: The total potential energy is nan, which is not finite. The LJ and
electrostatic contributions to the energy are 0 and 0, respectively. A
non-finite potential energy can be caused by overlapping interactions in
bonded interactions or very large or Nan coordinate values. Usually this is
caused by a badly- or non-equilibrated initial configuration, incorrect
interactions or parameters in the topology.
I have also energy minimized and equilibrated the system again but getting the same error.
Thank you,
Gaurav
Is the simulation stable with normal H masses and dt?
Yes, it was stable
I have tried HMR with CHARMM-GUI and it works fine. In the CHARMM-GUI topology file DNA restrains are mentioned as ===>
#ifdef POSRES
[ position_restraints ]
** 2 1 POSRES_FC_BB POSRES_FC_BB POSRES_FC_BB**
** 3 1 POSRES_FC_BB POSRES_FC_BB POSRES_FC_BB**
** 6 1 POSRES_FC_BB POSRES_FC_BB POSRES_FC_BB**
** 8 1 POSRES_FC_SC POSRES_FC_SC POSRES_FC_SC**
** 9 1 POSRES_FC_SC POSRES_FC_SC POSRES_FC_SC**
However, my topology files the restrains are like
[ position_restraints ]
7 ; atom type fx fy fz
8 1 1 1000 1000 1000
9 3 1 1000 1000 1000
10 6 1 1000 1000 1000
11 8 1 1000 1000 1000
12 9 1 1000 1000 1000
Could this be the reason for error? I assume I am making mistake which creating the topology files
Thanks,
Gaurav
All CHARMM-GUI is doing is allowing one to set a variable for the value of the restraint constant; they modify this in different stages of equilibration.
Please provide the .mdp file(s) for successful and crashing simulations.
HMR .mdp ==>
integrator = md
dt = 0.004
nsteps = 25000000
nstxtcout = 25000
nstvout = 25000
nstfout = 25000
nstcalcenergy = 100
nstenergy = 1000
nstlog = 1000
;
cutoff-scheme = Verlet
nstlist = 20
vdwtype = Cut-off
vdw-modifier = Force-switch
rvdw_switch = 1.0
rvdw = 1.2
rlist = 1.2
rcoulomb = 1.2
coulombtype = PME
;
tcoupl = V-rescale
tc_grps = Protein_NucleicAcid Water_and_ions
tau_t = 1.0 1.0
ref_t = 303.15 303.15
;
pcoupl = Parrinello-Rahman
pcoupltype = isotropic
tau_p = 5.0
compressibility = 4.5e-5
ref_p = 1.0
;
constraints = h-bonds
constraint_algorithm = LINCS
continuation = yes
;
nstcomm = 100
comm_mode = linear
comm_grps = Protein_NucleicAcid Water_and_ions
Suceessful Without HMR .mdp
title = Protein-ligand complex MD simulation
; Run parameters
integrator = md ; leap-frog integrator
nsteps = 25000000
dt = 0.002 ; 2 fs
; Output control
nstenergy = 5000 ; save energies every 10.0 ps
nstlog = 5000 ; update log file every 10.0 ps
nstxout-compressed = 5000 ; save coordinates every 10.0 ps
; Bond parameters
continuation = yes ; continuing from NPT
constraint_algorithm = lincs ; holonomic constraints
constraints = h-bonds ; bonds to H are constrained
lincs_iter = 1 ; accuracy of LINCS
lincs_order = 4 ; also related to accuracy
; Neighbor searching and vdW
cutoff-scheme = Verlet
ns_type = grid ; search neighboring grid cells
nstlist = 20 ; largely irrelevant with Verlet
rlist = 1.2
vdwtype = cutoff
vdw-modifier = force-switch
rvdw-switch = 1.0
rvdw = 1.2 ; short-range van der Waals cutoff (in nm)
; Electrostatics
coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics
rcoulomb = 1.2
pme_order = 4 ; cubic interpolation
fourierspacing = 0.16 ; grid spacing for FFT
; Temperature coupling
tcoupl = V-rescale ; modified Berendsen thermostat
tc-grps = Protein_NucleicAcid 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
pcoupl = Parrinello-Rahman ; pressure coupling is on for 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 is not used for proteins with the C36 additive FF
DispCorr = no
; Velocity generation
gen_vel = no ; continuing from NPT equilibration
It is not appropriate to couple groups separately for COM motion removal for a system like this. You may also have to look at coupling constants for e.g. the barostat, because you have tau_p = 5.0 in the failing run but 2.0 in the successful run. You should really try using the exact same .mdp file from the successful run in the HMR system (obviously with altered dt) to see if it runs, and if it doesn’t it is easier to determine what changes might be needed. In principle, there really shouldn’t be anything special about these settings for a reasonably well-equilibrated system.
It appears all your settings are for the CHARMM force field, but earlier you said you were using AMBER. If you are actually using AMBER, then the nonbonded settings are all incorrect, but that’s a fundamental problem in the inputs, not related to HMR at all.
In 2022.6 this feature seems to be buggy too.
pdb2gmx does create modified topology
grommp shows warning about oscillational period for leucine
mdrun fails with segfault with dt = 0.004 at first unrestrained simulation steps, does not produce errors with dt = 0.003
In the current main branch and the upcoming 2024 release, hydrogen mass repartitioning can now be done by grompp using a single mdp option using a standard topology.