Mdp parameters for non-bonded interaction using the AMBER force field

User discussions
, ,
1

GROMACS version: 2023.0
GROMACS modification: No

Hey,

I was wondering if there is a clear recommendation for the mdp parameters regarding the non-bonded interactions using the AMBER force field?

The GROAMCS manual recommends following parameters when using the CHARMM force field:

constraints = h-bonds

cutoff-scheme = Verlet

vdwtype = cutoff

vdw-modifier = force-switch

rlist = 1.2

rvdw = 1.2

rvdw-switch = 1.0

coulombtype = PME

rcoulomb = 1.2

DispCorr = no

However, there is no such recommendation in the manual for the AMBER force field.

Doing some research, I came across this post and this post from @alevilla .

In these posts following parameters are recommended for the AMBER force field:

constraints = h-bonds ; bonds involving H are constrained

rcoulomb = 1.0 ; short-range electrostatic cutoff (in nm)

rvdw = 1.0 ; short-range van der Waals cutoff (in nm)

vdw-modifier = Potential-shift-Verlet ; Amber specific

DispCorr = EnerPres ; account for cut-off vdW scheme

coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics

fourierspacing = 0.125 ; grid spacing for FFT

Doing some further research a came across this website that explained how to choose the appropriate cutoff distance. For the AMBER force field, it is 10 Ă… or 1 nm. The cutoff value depends on the cutoff which was used during the parametrization of the force field. As stated in the primary literature ff99SB, ff14SB, ff19SB and lipids21 were all parameterized using a cutoff value of 10 Ă… or 1 nm, respectively.

However, reading the primary literature I couldn’t find a reason why

vdw-modifier= Potential-shift
DispCorr = EnerPres

are used or why specifically these values are chosen.

I also checked the Sorin Lab website for ffAMBER but I couldn’t find an explanation there.

What is the reasoning for using

vdw-modifier = Potential-shift

DispCorr = EnePress

when using the AMBER force field in GROMACS and is there any literature or reference for this available?

Also, if these are the general non-bonded parameters in the mdp file when using the AMBER force field, wouldn’t it be useful to include them in the force field section of the GROMACS manual?

Additinoally, I also use

coulomb-modifier = Potential-shift

when using the AMBER force field. Again in there any literature or reference availbe for this parameter using the AMBER force field?

Cheers
Eddie

1 Like
2

Hey,

In case you are working with ff14SB, I think that one uses 8 Ă… cutoff for non-bonded interactions:
“Nonbonded interactions were calculated directly up to 8 Å.” from the ff14SB publication.

Please correct me if I am missing something.

Regards
Daniel

3

Hey Daniel,

thank you very much for this advice! You are right the cutoff for non-bonded interactions is 8 Ă… in ff14SB. I double cheeked the information for each force field and these are the used cutoffs.

ff99SB 8 Ă…:
„A cutoff of 8-Å was used for nonbonded interactions, and long-range electrostatic interactions were treated with the particle mesh Ewald (PME) method“

ff99SB-ILDN 10 Ă…:
„A cutoff of 10 Å was used for the Lennard-Jones interaction and the short-range electrostatic interactions.“

ff14SB 8 Ă…:
„ Nonbonded interactions were calculated directly up to 8 Å.“

ff19SB 10 Ă…:
„The direct space nonbonded interaction cutoff was 10.0 Å for explicit solvent simulations and 9999.0 (no cutoff) for implicit solvent simulation.“

However, there are separate force fields for lipids in the AMBER force field.

Lipid14 10 Ă…:
„PME was used to treat all electrostatic interactions with a real space cutoff of 10 Å. A long-range analytical dispersion correction was applied to the energy and pressure.“

Lipid21 10 Ă…:
„All restraints were removed and the PME was used to treat all electrostatic interactions with a real space cutoff of 10 Å“

In the current Amber manual (p. 51) it is stated that the behavior of lipids using the AMBER force field is highly sensitive to simulation conditions and a cutoff of 10 Ă… is recommended:

„A word of caution regarding barostat and cut off selection with lipid simulations. It is well known that lipids can be very sensitive to simulation conditions. Previous advice has been to use the Berendsen barostat and a 10A cutoff when simulating lipids with the AMBER Lipid force fields. This was due to bilayer deformation that was regularly seen with simulations run using the Monte Carlo Barostat. Recent work by Gomez et al.[91] investigated this behavior and determined that there are issues when using the MC barostat with hard Lennard-Jones cutoffs. The issue affects all simulations but is most obvious when simulating lipid bilayers. Gomez et al recommend use of an LJ force switch when running simulations with the MC barostat.“

Based on that I would use a cutoff of 10 Ă… as soon as there is a membrane system.

Still I couldn’t find the information that recommends the use of

vdw-modifier = Potential-shift 
DispCorr = EnePress

Additionally, I tried CHARMM-GUI to generate the inputs and these are the parameter in the mdp file:

cutoff-scheme           = Verlet
nstlist                 = 20
rlist                   = 0.9
vdwtype                 = Cut-off
vdw-modifier            = None
DispCorr                = EnerPres
rvdw                    = 0.9
coulombtype             = PME
rcoulomb                = 0.9

I would be very interested in the most suited vdw-modifier for the AMBER force field.

Best
Eddie

1 Like
4

Hi Eddie,

maybe this helps you (from the AMBER22 manual):

In the section: 21.7.2. Particle Mesh Ewald

The Particle Mesh Ewald (PME) method is always on … Please take care in changing any values from their defaults."

Then further down in the section describing the options of PME:

vdwmeth: Determines the method used for van der Waals interactions beyond those included in the direct sum. A value of 0 includes no correction; the default value of 1 uses a continuum model correction for energy and pressure.

We can probably assume this is equivalent to gmx DispCorr.

Further in: 21.7. Potential function parameters

fswitch: When off, fswitch<=0 , uses a truncation cutoff. When on fswitch>0, sets a force switching region where the force cutoff smoothly approaches 0 between the region of the fswitch value to the cut value. Force values below the fswitch value follow the standard Lennard-Jones force. Default is -1. This option is not supported for use with GB (i.e., only igb=0 and ntb>0), nor is it compatible with the 12-6-4 Lennard-Jones model (lj1264=1). Due to performance regressions (about 20%) with running with the force switching on, it is recommended that simulations run with fswitch off unless using a force field that requires or recommends using the force switch.

Therefore, the only gmx vdw-modifier that fits the current default behavior of AMBER:

vdw-modifier = Potential-shift
DispCorr = EnerPress

These are more general recommendations and not really force field specific. However, I guess we are looking for a reference one can cite that explicitly states this. There are references scattered in this section of the AMBER manual and one might find something there.

And just to make sure we are on the same page, the ff14sb publication also implies (again, please correct me if I am wrong):

rcoulomb = 0.8 ; short-range electrostatic cutoff (in nm)

What makes all of this confusing to me is that there are publications out there (at least one) that use GROMACS in combination with the ff14sb forcefield without stating the specific options used.

Still, until ff14sb is officially supported with recommended parameters we probably should no be working with it if we want to publish. As far as I know, the GROMACS devs have stated that this depends on a validation that should be performed by the users once, which is understandable given that GROMACS is free and open source (which is great!). However, that this has not happened yet is astonishing. Are there any test systems/protocols in place that might constitute a proof? Interestingly, some work in this direction on a modified version for RNA/DNA simulations of ff14sb is done here.

I welcome any corrections if something I said is inaccurate.

Regards,
Daniel

5

Hi Daniel,

thank you very much for your informative answers!

Regarding the DispCorr:

I guess the DispCorr is the GROMACS equivalent of the AMBER’s vdwmeth. As pointed out in the AMBER22 manual vdwmeth is set to 1 by default using a continuum model correction for energy and pressure.

So far, I could only find in the ff14SB publication that a continuum model correction for energy and pressure was used:

A continuum model correction for energy and pressure was applied to long-range van der Waals interactions. The production timesteps were 2 fs for NVT and 1 fs for NVE.

I think this would also be an excellent source to cite.

Additionally, I found section 3.5.1 of the AMBER22 manual:

The rest of parameters are set to current Amber defaults; note that these include accounting for the van der Waals interactions beyond the cut-off via a continuum model (vdwmeth=1).

I assume that the defaults apply to both the AMBER MD engine and the AMBER force field.

Regarding the vdw-modifier:

I think AMBER’s fswitch corresponds more to GROMACS’ vdw-modifier options None or Force-switch. But correct me if I’m wrong.

In the current GROMACS manual in the section Non-bonded Interactions it says:

In the Verlet cut-off scheme, the potential is shifted by a constant such that it is zero at the cut-off distance.

Further down in this section it states:

This can either be done as a switching function that changes the shape of the potential and force over a small range, or by shifting the entire potential by a constant factor such that it becomes zero at the cutoff. The advantage of the shifted interaction modification is that it does not influence the force at all, and since only forces enter the equations of motion it will not influence the dynamics of the system. The drawback is that the total change in the potential is larger. Presently GROMACS only supports this shifted modification, and it is even applied by default (but possible to turn off). Note that we also shift the direct-space component of the PME interaction

I would assume that explains:

vdw-modifier=Potential-shift

Then in the Long Range Electoristatics sectiontion of the GROMACS manual it states for the PME method:

With the Verlet cut-off scheme, the PME direct space potential is shifted by a constant such that the potential is zero at the cut-off. This shift is small and since the net system charge is close to zero, the total shift is very small, unlike in the case of the Lennard-Jones potential where all shifts add up. We apply the shift anyhow, such that the potential is the exact integral of the force.

I find this very interesting. This is almost identical to coulomb-modifier=Potential-shift in the mdp options section of the GORMACS manual:

Shift the Coulomb potential by a constant such that it is zero at the cut-off. This makes the potential the integral of the force. Note that this does not affect the forces or the sampling.“

If I understand this correctly this means that once cutoff-scheme=Verlet is used vdw-modifiert=Potential-shift and coulomb-modifier=Potential-shift are used by default.

So both settings vdw-modifier=Porential-shift and coulomb-modifier=Potential-shift are rather determined by GROMACS as MD engine than by the AMBER force field.

Regarding ff14SB:

I agree that in the ff14SB publication

rvdw= 0.8
rcoulomb= 0.8

are implied. Yet ff19SB and the lipid21 use:

rvdw= 1.0
rcoulomb = 1.0

Best
Eddie

6

@alevilla @milosz.wieczor
I am interested in conducting DNA simulations using the bsc0 and bsc1 force fields. However, I have encountered a challenge related to the nonbonded parameters, and I find the information somewhat unclear.
In articles 1 and 2, a cutoff of 10 angstroms is mentioned.I also refer to another article 3 (which serves as a reference for amberparm99) stating that the cutoff is 9 angstroms. Given that bsc0 is developed based on the amberparm99 force field, I am unsure which cutoff value I should choose for my simulations.
Could you kindly provide clarification on this matter? Additionally, can I use these parameters (expect rcoulomb = 1.0 and rvdw = 1 that should be different) for both bsc0 and bsc1 force fields?

  • constraints = h-bonds
  • rcoulomb = 1.0
  • rvdw = 1.0
  • Potential-shift-Verlet (vdw-modifier = Potential-shift-Verlet)
  • EnerPres (DispCorr = EnerPres)
  • coulombtype = PME
  • fourierspacing = 0.125

I greatly appreciate your time and assistance in resolving these queries.
Thank you and best regards,

Ref:

  1. Galindo-Murillo, R., Robertson, J. C., Zgarbova, M., Sponer, J., Otyepka, M., Jurecka, P., & Cheatham III, T. E. (2016). Assessing the current state of amber force field modifications for DNA. Journal of chemical theory and computation , 12 (8), 4114-4127.
  2. Love, O., Galindo-Murillo, R., Zgarbová, M., Šponer, J., Jurečka, P., & Cheatham III, T. E. (2023). Assessing the Current State of Amber Force Field Modifications for DNA─ 2023 Edition. Journal of Chemical Theory and Computation .
  3. Wang, J., Cieplak, P., & Kollman, P. A. (2000). How well does a restrained electrostatic potential (RESP) model perform in calculating conformational energies of organic and biological molecules?. Journal of computational chemistry , 21 (12), 1049-1074
7

Hi Fatemeh,

I personally used vdw cut-offs of 1.0 and 1.2 nm for nucleic acid simulations with bsc0/bsc1 corrections, and am not aware of any specific artifacts that would be introduced by any value between 0.8 and 1.2 nm. It’s indeed a very prominent topic in lipid FFs as the small differences multiply in the dense environment of a bilayer to affect phase transitions etc, but DNA/RNA are generally robust to these details. (Or at least considered to be, unless someone proves otherwise.)

I also asked around in our group (where bsc0/1 were developed) and there seems to be no known issues with that.

I never used DispCorr with these FFs though, and I don’t think they should be applied.

8

thank you so much