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

User discussions
1

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
2

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.

3

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.

4

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

5

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.

6

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.

7

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).

8

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.

9

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?

10

@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.

11

Hi,

I am facing the same problem. Iwas wondering if you could solve it?

```
integrator = md
dt = 0.004
nsteps = 250000
nstxout-compressed = 50000
nstxout = 0
nstvout = 0
nstfout = 0
nstcalcenergy = 100
nstenergy = 1000
nstlog = 1000

mass-repartition-factor = 3 ;

; cutoff-scheme = Verlet
nstlist = 20
rlist = 1.2
vdwtype = Cut-off
vdw-modifier = Force-switch
rvdw_switch = 1.0
rvdw = 1.2
coulombtype = PME
rcoulomb = 1.2
;
tcoupl = v-rescale
tc_grps = SOLU_MEMB SOLV
tau_t = 1.0 1.0
ref_t = 310 310
;
pcoupl = C-rescale
pcoupltype = semiisotropic
tau_p = 5.0
compressibility = 4.5e-5 4.5e-5
ref_p = 1.0 1.0
;
constraints = h-bonds
constraint_algorithm = LINCS
continuation = yes
;
nstcomm = 100
comm_mode = linear
comm_grps = System

```

I cant get ridd off this note

```

ommand line:
gmx grompp -f ../step7_production.mdp -o step7_1.tpr -c ../equilibration/step6.6_equilibration.gro -p ../topol.top -n ../index.ndx

Generating 1-4 interactions: fudge = 1

NOTE 1 [file topol.top, line 36]:
The bond in molecule-type POPC between atoms 31 C21 and 32 O22 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.

Number of degrees of freedom in T-Coupling group SOLU_MEMB is 406549.19
Number of degrees of freedom in T-Coupling group SOLV is 1091502.88

There was 1 NOTE

GROMACS reminds you: “Heavier-than-air flying machines are impossible.” (Lord Kelvin, President of Royal Society, 1895.)

````

12

This doesn’t look like a suitable system for a 4 fs time step with mass repartitioning. Are the examples in the literature of people doing this?

13

For GPCRs, it is a quite common protocol:

https://www.nature.com/articles/s41592-020-0884-y

https://www.nature.com/articles/s41467-025-57034-y

https://www.nature.com/articles/s41467-024-54103-6#Sec9

14

Ah, I see now that at least one of those publications uses AceMD which by default turns on constraints for all bonds with a time step of more than 2 fs. That would also work in GROMACS, but we advise against constraining all bonds, as the force field has been parametrized with constraints on bonds involving hydrogens only. Constraining all bonds alters the dihedral distributions and transition rates.

You could use multiple time stepping with PME every second step for a modest gain in performance.

15

Thanks for your reply @hess . Would you have any reference about this

```
Constraining all bonds after the dihedral distributions and transition rates.
```

I wanted to use a tmiestep of 4fs because I'm trying to reproduce something that I saw with ACEMD on GROMACs
16

I don’t actually know of any references. This is common knowledge among several of my colleagues.

If you want to reproduce results from AceMD, you need to turn on constraints for all bonds.

17

Thanks!!