Simulating Protein in 8M Urea & Water: Forces Won't Converge

User discussions
1

GROMACS version: 2020.3
GROMACS modification: No

Hi all,

I have been attempting to simulate a GFP construct (with a zwitterionic TYG chromophore) in a mixed solvent box with ~8M urea and water in the Amber99 force field, using the TIP3P water model. When I run energy minimization using the steepest descent algorithm, the forces will not converge below 1000 kJ/mol/nm, though the system converges to machine precision. I am hoping someone could advise me on how to fix this such that the forces will converge below 1000 kJ/mol/nm (I am an undergraduate - not an expert).

The specifics of my protocol is as follows:

  • I use pdb2gmx to generate a standard topology file, position restraint file and coordinate file
  • I use editconf to place the GFP construct in the center of a periodic box such that there is a distance of 1.5 nm between the construct and the box edges
  • I use insert-molecules to place 2700 urea molecules in the box
  • I solvate with 19000 molecules of TIP3P using the -maxsol option of solvate
  • I neutralize the charge of the system with 6 NA ions using genion

Initially, I hypothesized that the forces may be too large as a result of the pressure in the box being too high. As such, I ran short, successive NPT equilibrations (using isotropic Berendsen pressure coupling, PME electrostatics, electrostatic cutoff of 1.2 nm, VDW cutoff of 1.2 nm) to allow the box to expand and contract. This was only successful if I used a time step of 0.4 fs - anything larger produced too many LINCS errors. I ultimately equilibrated the system for 1 ns in total, and was able to reach a desirable potential energy value as well as a desirable average pressure value.

That being said, I remain bothered by the fact that I cannot get the forces to converge below 1000 kJ/mol/nm, and I would like to use a time step that is larger than 0.4 fs if at all possible.

Hence, if anyone has advice for me, I am all ears, and I am eager to learn.

Best regards,
Elise

Elise White
Rensselaer Polytechnic Institute ’22
Undergraduate Researcher, Bystroff Lab
Co-Op Research Assistant, RNA Institute

2

Hi,

When did you energy minimize? Usually, one perform energy minimization to remove “bad” contact before and then run the MD simulations to get the system to the correct temperature first and then pressure.

The time step depends on the force field you are using. Usually is 2 fs for biomolecular force field. Note the setting for bond contraints depends also on the force field setting (and they go together with time step). You can check what is the setting required for the force field you are using.

Best regards
Alessandra

3

Hi Alessandra,

I initially attempted to run this simulation following the standard protocol, which you described - i.e., I ran energy minimization, NVT equilibration, NPT equilibration and then production. However, since the forces remained large after energy minimization, I was unable to successfully run NVT equilibration or NPT equilibration using a 2 fs time step without generating several LINCS warnings.

To troubleshoot this, I decided to run NPT equilibration first to allow the box to expand and contract, thinking that stabilizing (and lowering) the pressure acting on the system would also lower the forces. However, as I said above, I found that this was not conducive to solving my problem of not being able to run equilibrations or production using a time step greater than 0.4 fs.

I understand that 2 fs is standard for a biomolecular force field (such as Amber99, which I am using). I have ran many simulations using this time step and force field on other systems - including this one in water - and have not encountered a situation where I am forced to use a much smaller time step.

I am simply looking to speed up the simulation process, to understand why these LINCS warnings are being generated, and to understand why the forces in my system remain so large after energy minimization.

4

Hi Elise,

Happy New Year! Post the mdp file for your NVT equilibration. There might be a bad setting.

-AB

5

Hi Anthony! It’s good to hear from you!

Here are the mdp files for my NVT and NPT equilibration runs. Thank you for being willing to look at them. Since new users can’t upload files, I have copied and pasted their contents below:

  1. NVT.mdp

title = Amber99 LII-GFP NVT Equilibration
define = -DPOSRES ; position restrain the protein
; Run parameters
integrator = md ; leap-frog integrator
nsteps = 2500000 ; 0.4 * 2500000 = 1000 ps = 1 ns
dt = 0.0004 ; 0.4 fs

; Output control
nstxout = 500 ; save coordinates every 1.0 ps
nstvout = 500 ; save velocities every 1.0 ps
nstenergy = 500 ; save energies every 1.0 ps
nstlog = 500 ; update log file every 1.0 ps

; Bond parameters
continuation = no ; first dynamics run
constraint_algorithm = lincs ; holonomic constraints
constraints = h-bonds ; bonds involving H are constrained
lincs_iter = 1 ; accuracy of LINCS
lincs_order = 4 ; also related to accuracy

; Nonbonded settings
cutoff-scheme = Verlet ; Buffered neighbor searching
ns_type = grid ; search neighboring grid cells
nstlist = 10 ; 20 fs, largely irrelevant with Verlet
rcoulomb = 1.2 ; short-range electrostatic cutoff (in nm)
rvdw = 1.2 ; short-range van der Waals cutoff (in nm)
DispCorr = EnerPres ; account for cut-off vdW scheme

; Electrostatics
coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics
pme_order = 4 ; cubic interpolation
fourierspacing = 0.10 ; grid spacing for FFT

; Temperature coupling
tcoupl = V-rescale ; modified Berendsen thermostat
tc-grps = Protein Non-Protein ; 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 is off
pcoupl = no ; no pressure coupling in NVT ensemble

; Periodic boundary conditions
pbc = xyz ; 3D PBC

; Velocity generation
gen_vel = yes ; assign velocities from Maxwell distribution
gen_temp = 300 ; temperature for Maxwell distribution
gen_seed = -1 ; generate a random seed

  1. NPT.mdp

title = Amber99 LII-GFP NPT Equilibration
define = -DPOSRES ; position restrain the protein
; Run parameters
integrator = md ; leap-frog integrator
nsteps = 2500000 ; 0.4 * 2500000 = 1000 ps = 1 ns
dt = 0.00004 ; 0.4 fs

; Output control
nstxout = 500 ; save coordinates every 1.0 ps
nstvout = 500 ; save velocities every 1.0 ps
nstenergy = 500 ; save energies every 1.0 ps
nstlog = 500 ; update log file every 1.0 ps

; Bond parameters
continuation = no ; Restarting after NVT
constraint_algorithm = lincs ; holonomic constraints
constraints = h-bonds ; bonds involving H are constrained
lincs_iter = 1 ; accuracy of LINCS
lincs_order = 4 ; also related to accuracy

; Nonbonded settings
cutoff-scheme = Verlet ; Buffered neighbor searching
ns_type = grid ; search neighboring grid cells
nstlist = 10 ; 20 fs, largely irrelevant with Verlet scheme
rcoulomb = 1.2 ; short-range electrostatic cutoff (in nm)
rvdw = 1.2 ; short-range van der Waals cutoff (in nm)
DispCorr = EnerPres ; account for cut-off vdW scheme

; Electrostatics
coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics
pme_order = 4 ; cubic interpolation
fourierspacing = 0.10 ; grid spacing for FFT

; Temperature coupling
tcoupl = V-rescale ; modified Berendsen thermostat
tc-grps = Protein Non-Protein ; 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 = Berendsen ; Pressure coupling algorithm
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
refcoord_scaling = com
; Periodic boundary conditions
pbc = xyz ; 3D PBC
; Velocity generation
gen_vel = no ; Velocity generation is off

Elise

6

It’s good to hear from you too!

I don’t see anything obviously wrong with your parameters. If you haven’t checked out this page yet, I would try some of these suggestions.

http://www.gromacs.org/Documentation_of_outdated_versions/Terminology/Blowing_Up

Where I would start looking:

Make sure your urea topology is correct. I have never used urea in a sim, so forgive my ignorance. Is that standard in the AMBER force-field or did you generate it yourself?

Try equilibrating the protein in water without urea. If that is successful I would double check your urea topology.

Confirm the urea is not the issue by equilibrating some urea in water without any protein.

Sometimes genion places ions at an unfortunate place and regenerating your starting box can fix it.

Sorry I couldn’t be more help.
Good Luck,
AB

7

Well that’s certainly reassuring to hear!

As far as I am aware, my urea topology is correct. It is standard in the Amber99 force field and comes with the latest version of GROMACS.

I have successfully equilibrated and simulated the protein in water (without urea). I have also successfully equilibrated a mixed solvent box of urea and water without the protein.

I am wondering if perhaps I need to specify an ion concentration when I run genion, as opposed to simply adding the correct number of counterions needed to neutralize the charge on the system.

I am also curious as to whether I should change my LINCS constraints from “h-bonds” to “all-bonds” - it seems many people suggest using the latter option instead.

You were very helpful - I really appreciate you taking a look at this.

Elise

8

Constraints are a function of the force field being used; it’s not really a freely tunable parameters. For AMBER, this should be h-bonds.

9

Ah, that makes sense. I didn’t know that. Thank you so much for your help!

10

It would be nice to add this clarification to the mdp options page (ideally for all packaged ff families). Especially since older tutorials and mailing list answers recommended “all-bonds”.

11

We have not put such information in the manual because it may change with different force field variants, and the user can always find this by reading the primary literature for the chosen force field. That’s a requisite first step before using any force field - know how it is implemented!

Unfortunately, yes, there are incorrect answers on the mailing list and in this forum. Such is the nature of crowd-sourcing help.