How does hydrogen mass repartitioning with an increased time step of 4 fs affect other mdp options?

GROMACS version: 2024.2
GROMACS modification: No

Hey everyone,

I ran some test simulations using the new mdp option “mass-repartition-factor,” and it worked fine with my other default mdp options. I have a question regarding the impact of increasing the time step from 2 fs to 4 fs on the following MDP parameters:

  • nstcalcenergy
  • nstlist
  • nstcomm
  • tau-p
  • tau-t

For nstcalcenergy, I expect it might be necessary to decrease the default number of steps from 100 to 50 due to the increased time step. Similarly, I am considering adjustments for nstlist and nstcomm.

However, I believe tau-p and tau-t can remain at 5 ps and 1 ps, respectively, as they are not affected by the increased time step.

Is it necessary to adjust these parameters when using hydrogen mass repartitioning, or can I leave them at their default values? Were there any investigations during the development phase that addressed this?

For reference, here are my default MDP settings for simulations of membrane proteins with the AMBER force field:

; Run parameters
integrator             = md                               
nsteps                 = 50000000                         
dt                     = 0.002                            
init_step = 0

; Output control
nstxout                = 0                                
nstvout                = 0                                
nstcalcenergy          = 100                              
nstenergy              = 250000                           
nstxout-compressed     = 250000                           
compressed-x-precision = 1000                             
compressed-x-grps      = System                           
nstlog                 = 250000                           

; Bond parameters
constraint-algorithm   = lincs                            
constraints            = h-bonds                          
lincs-iter             = 1                                
lincs-order            = 4                                

; Neighborsearching
cutoff-scheme          = Verlet                           
nstlist                = 20                               
rlist                  = 1.0                              

;van der Waals Interactions
vdwtype                = Cut-off                          
vdw-modifier           = Potential-shift                  
rvdw-switch            = 0                                
rvdw                   = 1.0                              
DispCorr               = EnerPres                         

; Electrostatics
coulombtype            = PME                              
pme-order              = 4                                
fourierspacing         = 0.125                            
coulomb-modifier       = Potential-shift                  
rcoulomb-switch        = 0                                
rcoulomb               = 1.0                              

; Temperature coupling is on
tcoupl                 = V-rescale                        
tc-grps                = SOLU        MEMB           SOLV  
tau-t                  = 1.0         1.0            1.0   
ref-t                  = 310         310            310   
nsttcouple             = -1                               

; Pressure coupling is on
pcoupl                 = C-rescale                        
pcoupltype             = semiisotropic                    
tau-p                  = 5.0                              
ref-p                  = 1.0 1.0                          
compressibility        = 4.5E-5 4.5E-5                    
nstpcouple             = -1                               

;Position restraints
;refcoord-scaling       = no                             

;Center of Mass Motion
comm-mode             = linear                            
nstcomm               = 100                               
comm-grps             = SOLU_MEMB SOLV                    

; Periodic boundary conditions
pbc                    = xyz                              

; Continuation of velocity generation 
gen-vel                = no                               
continuation           = yes                              
gen-temp               = 310                              

nstlist is optimized by mdrun, so not relevant.
nstcomm is uncritical, I would leave it at 100 to avoid performance overhead.
nstcalcenergy only affects the statistics of the energy file output, no other results. It affects performance a little bit.

Your tau values look fine. nsttcouple will be lowered when doubling the time-step, but I wouldn’t suggest to use tau-t=2 ps, even if it’s probably fine in most cases.

Just as a follow-up to your parameters. Why are you running three separate temperature coupling groups? Is SOLU very large? I think you should run with two groups, SOLU_MEMB SOLV, or even just use one group.

Thank you Berk and Magnus!

@hess : Perfect, then I will leave the values as they are!

@MagnusL : This is a mdp file or a protein (GPCR) in membrane system. The system is around 100k atoms big. I generated the system using CHARMM-GUI (AMBER force field / GROMACS input files). By default, CHARMM-GUI creates 3 temperature coupling groups (SOLU, SOLV, MEMB). When I modified the mdp file I just went with it. Would there be any benefit using just 2 groups (SOLV, SOLU_MEMB)?

Just as a follow up question - when I grompp’d a mdp file with

mass-repartition-factor = 3

to generate a trp file I received following note:

NOTE 1 [file topol.top, line 28]:
  The bond in molecule-type MOL between atoms 30 O4 and 39 C10 has an
  estimated oscillational period of 2.1e-02 ps, which is less than 10 times
  the time step of 4.0e-03 ps.
  Maybe you forgot to change the constraints mdp option.

During the simulations I receive following warning:

Step 8151, time 32.604 (ps)  LINCS WARNING
relative constraint deviation after LINCS:
rms 0.000118, max 0.005614 (between atoms 6511 and 6513)
bonds that rotated more than 30 degrees:
 atom 1 atom 2  angle  previous, current, constraint length
   6511   6514   30.7    0.1091   0.1084      0.1090

I did some research and noticed that other users experience similar problems:

Is there any solution to it? I use the AMBER force field and from my understanding I can only constrain the bonds to hydrogen atoms (constraints = h-bonds). So, using (constraints = all-bonds) would not be an option since the force field was designed for that.

Regards
Eddie

Using a factor 3 will only allow you to increase the time step by about a factor of sqrt(3). I think that AMBER uses a factor of 5.

In general, temperature coupling groups should be large - or you can get very large fluctuations of temperature in the group. If parts of the system would not “communicate temperature/velocities” it might be good to couple them separately. In most cases it’s fine using one temperature coupling group, but possibly separate comm-grps, e.g., solvent, non-solvent.

To address your LINCS warnings (and possibly crashes), try dt = 0.003. It will not give you the same speed-up, but it should be more stable. You could possibly try mass-repartition-factor = 4 as well. But that will at least not fix the note from gmx grompp, it might even turn into a warning.

Thanks again for the fast reply!

@hess: I tired mass-repartition-factor = 4 but I got this error from gmx grompp :

ERROR 1 [file topol.top, line 28]:
  Light atoms are bound to at least one atom that has a too low mass for
  repartioning

In the original publication for AMBER they used a factor of 3. In one of your earlier publications a factor of 4 was used.

@MagnusL : dt = 0.003 did the trick. I still received the note from gmx grompp but the LINCS warnings disappeared and the simulations did not crash (so far).

Good to hear that it works. It’s unfortunate that you cannot always increase the timestep with a factor 2. I think a mass repartition factor 3 is a good balance, and I’ve often used that successfully with dt=0.003, even when there is a note from grompp regarding the oscillational period. FYI, it will be a warning if you use a far too high time step, so it often works well when you get a note, but not always, as you saw in your case.

Sorry, I’m not on top of the common mass repartitioning settings.

But what time step is then typically used with Amber when using a factor of 3?

@hess : Usually a time step of 4 fs is used.

Also setting lincs-warnangle = 90 helps to prevent simulations crashing due to too many LINCS warnings in the log file. That way mass-repartition-factor = 3 and dt = 4 worked. However, I need to evaluate how stable these simulations are long term.