# Simulating Protein in 8M Urea & Water: Forces Won't Converge **URL:** \ **Category:** User discussions\ **Created:** [January 6, 2021, 1:27am UTC](https://gromacs.bioexcel.eu/t/simulating-protein-in-8m-urea-water-forces-wont-converge/1389 "2021-01-06T01:27:00Z")\ **Posts on this page:** 11\ **Page:** 1
**Author:** ![elisenwhite](https://avatars.discourse-cdn.com/v4/letter/e/b19c9b/32.png) [@elisenwhite](https://gromacs.bioexcel.eu/u/elisenwhite)\ **Post date:** [January 6, 2021, 1:27am UTC](https://gromacs.bioexcel.eu/t/simulating-protein-in-8m-urea-water-forces-wont-converge/1389/1 "2021-01-06T01:27:00Z")
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 ---
**Author:** ![alevilla](https://dub1.discourse-cdn.com/flex017/user_avatar/gromacs.bioexcel.eu/alevilla/32/439_2.png) [@alevilla](https://gromacs.bioexcel.eu/u/alevilla)\ **Post date:** [January 7, 2021, 11:17am UTC](https://gromacs.bioexcel.eu/t/simulating-protein-in-8m-urea-water-forces-wont-converge/1389/2 "2021-01-07T11:17:11Z")
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 ---
**Author:** ![elisenwhite](https://avatars.discourse-cdn.com/v4/letter/e/b19c9b/32.png) [@elisenwhite](https://gromacs.bioexcel.eu/u/elisenwhite)\ **Post date:** [January 7, 2021, 6:43pm UTC](https://gromacs.bioexcel.eu/t/simulating-protein-in-8m-urea-water-forces-wont-converge/1389/3 "2021-01-07T18:43:09Z")
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. ---
**Author:** ![anthonycbishop](https://avatars.discourse-cdn.com/v4/letter/a/8e8cbc/32.png) [@anthonycbishop](https://gromacs.bioexcel.eu/u/anthonycbishop)\ **Post date:** [January 7, 2021, 10:36pm UTC](https://gromacs.bioexcel.eu/t/simulating-protein-in-8m-urea-water-forces-wont-converge/1389/4 "2021-01-07T22:36:57Z")
Hi Elise, Happy New Year! Post the mdp file for your NVT equilibration. There might be a bad setting. -AB ---
**Author:** ![elisenwhite](https://avatars.discourse-cdn.com/v4/letter/e/b19c9b/32.png) [@elisenwhite](https://gromacs.bioexcel.eu/u/elisenwhite)\ **Post date:** [January 7, 2021, 11:56pm UTC](https://gromacs.bioexcel.eu/t/simulating-protein-in-8m-urea-water-forces-wont-converge/1389/5 "2021-01-07T23:56:04Z")
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 ---
**Author:** ![anthonycbishop](https://avatars.discourse-cdn.com/v4/letter/a/8e8cbc/32.png) [@anthonycbishop](https://gromacs.bioexcel.eu/u/anthonycbishop)\ **Post date:** [January 8, 2021, 5:19pm UTC](https://gromacs.bioexcel.eu/t/simulating-protein-in-8m-urea-water-forces-wont-converge/1389/6 "2021-01-08T17:19:23Z")
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](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 ---
**Author:** ![elisenwhite](https://avatars.discourse-cdn.com/v4/letter/e/b19c9b/32.png) [@elisenwhite](https://gromacs.bioexcel.eu/u/elisenwhite)\ **Post date:** [January 8, 2021, 8:25pm UTC](https://gromacs.bioexcel.eu/t/simulating-protein-in-8m-urea-water-forces-wont-converge/1389/7 "2021-01-08T20:25:35Z")
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 ---
**Author:** ![jalemkul](https://dub1.discourse-cdn.com/flex017/user_avatar/gromacs.bioexcel.eu/jalemkul/32/18_2.png) [@jalemkul](https://gromacs.bioexcel.eu/u/jalemkul)\ **Post date:** [January 8, 2021, 9:34pm UTC](https://gromacs.bioexcel.eu/t/simulating-protein-in-8m-urea-water-forces-wont-converge/1389/8 "2021-01-08T21:34:12Z")
> [@elisenwhite](#): > > 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. 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`. ---
**Author:** ![elisenwhite](https://avatars.discourse-cdn.com/v4/letter/e/b19c9b/32.png) [@elisenwhite](https://gromacs.bioexcel.eu/u/elisenwhite)\ **Post date:** [January 11, 2021, 8:35pm UTC](https://gromacs.bioexcel.eu/t/simulating-protein-in-8m-urea-water-forces-wont-converge/1389/9 "2021-01-11T20:35:45Z")
Ah, that makes sense. I didn’t know that. Thank you so much for your help! ---
**Author:** ![shelved](https://avatars.discourse-cdn.com/v4/letter/s/e68b1a/32.png) [@shelved](https://gromacs.bioexcel.eu/u/shelved)\ **Post date:** [January 13, 2021, 5:17am UTC](https://gromacs.bioexcel.eu/t/simulating-protein-in-8m-urea-water-forces-wont-converge/1389/10 "2021-01-13T05:17:24Z")
> [@jalemkul](#): > > For AMBER, this should be `h-bonds`. 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”. ---
**Author:** ![jalemkul](https://dub1.discourse-cdn.com/flex017/user_avatar/gromacs.bioexcel.eu/jalemkul/32/18_2.png) [@jalemkul](https://gromacs.bioexcel.eu/u/jalemkul)\ **Post date:** [January 13, 2021, 12:02pm UTC](https://gromacs.bioexcel.eu/t/simulating-protein-in-8m-urea-water-forces-wont-converge/1389/11 "2021-01-13T12:02:48Z")
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.