MD minimization

User discussions
1

GROMACS version: 2018.1
GROMACS modification: Yes/No
I have a system with more than 100,000 atoms. After minimization using the following parameters:

; minim.mdp - used as input into grompp to generate em.tpr
; Parameters describing what to do, when to stop and what to save
integrator = steep ;CG ; Algorithm (steep = steepest descent minimization)
;define = -DPOSRES
emtol = 100.0 ; Stop minimization when the maximum force < 1000.0 kJ/mol/nm
emstep = 0.01 ; Energy step size
nsteps = 2000 ; Maximum number of (minimization) steps to perform

; Parameters describing how to find the neighbors of each atom and how to calculate the interactions
nstlist = 10 ; Frequency to update the neighbor list and long range forces
ns_type = grid ; Method to determine neighbor list (simple, grid)
rlist = 1.2 ; Cut-off for making neighbor list (short range forces)
cutoff-scheme = verlet
coulombtype = PME ; Treatment of long range electrostatic interactions
rcoulomb = 1.2 ; Short-range electrostatic cut-off
rvdw = 1.2 ; Short-range Van der Waals cut-off
pbc = xyz ; Periodic Boundary Conditions
periodic-molecules = yes

I got this result (See below). Kindly advise on how to get a system that satisfies -ve PE and Fmax < tolerance (emtol)

Steepest Descents:
Tolerance (Fmax) = 1.00000e+02
Number of steps = 2000

Energy minimization reached the maximum number of steps before the forces
reached the requested precision Fmax < 100.

writing lowest energy coordinates.

Back Off! I just backed up em.gro to ./#em.gro.1#

Steepest Descents did not converge to Fmax < 100 in 2001 steps.
Potential Energy = -1.7481004e+07
Maximum force = 9.7902646e+03 on atom 11248
Norm of force = 2.4687517e+01

GROMACS reminds you: “Everything Must Go” (Red Hot Chili Peppers)

2

Hi,

This means that your energy minimization is completed. Check the potential energy of your system using energy command.

3

I have checked it. The potential energy is -ve. My issue is that the minimization satisfied the conditon " Epot should be negative" but not "Fmax < emtol (100 KJ/mol). On Justin’s mdtutorial page, he wrote this “It is possible to arrive at a reasonable Epot with Fmax > emtol. If this happens, your system may not be stable enough for simulation. Evaluate why it may be happening, and perhaps change your minimization parameters (integrator, emstep, etc).”. I have modified the emstep from 1000 to 2000. The Fmax change from 4.7062959e+03 to 9.7902646e+03 KJ/mol. At emstep of 1000, the PE is -1.7104652e+07. I am not sure about how to change integrator or go from here.

4

Try several rounds of minimization with different algorithms (e.g. try conjugate gradient next) and perhaps use a smaller value of emstep to step more carefully through the potential energy surface.

5

Thanks, Dr Justin. I am currently working on using CG.

While reading the GMX page, i saw this In GROMACS conjugate gradient can not be used with constraints, including the SETTLE algorithm for water 47, as this has not been implemented. If water is present it must be of a flexible model, which can be specified in the mdp file by define = -DFLEXIBLE . My system has water (SPCE) and I am working a position restraint file attached to the topology file. Does it mean the porse.itp won’t work and I need to use define = -DFLEXIBLE in the mdp file because I have water.

Also, I will be happy to have a sample mdp file where CG was implemented.

6

The .mdp file you posted above has position restraints disabled, but in any case, restraints and constraints are separate. As for the CG .mdp file, it’s simply:

integrator = cg
define = -DFLEXIBLE

and everything else is the same.

You will then want to likely perform another round of steepest descent EM with rigid water in case the use of flexible water leads to distorted structures, which can happen from time to time.

7

I tried your suggestion by modifiying the mdp files as follows:

minim.mdp - used as input into grompp to generate em.tpr
; Parameters describing what to do, when to stop and what to save
integrator = CG ;steep CG ; Algorithm (steep = steepest descent minimization)
define = -DFLEXIBLE
emtol = 100.0 ; Stop minimization when the maximum force < 1000.0 kJ/mol/nm
emstep = 0.01 ; Energy step size
nstcgsteep = 1000
nsteps = 1000 ; Maximum number of (minimization) steps to perform

; Parameters describing how to find the neighbors of each atom and how to calculate the interactions
nstlist = 10 ; Frequency to update the neighbor list and long range forces
ns_type = grid ; Method to determine neighbor list (simple, grid)
rlist = 1.2 ; Cut-off for making neighbor list (short range forces)
cutoff-scheme = verlet
coulombtype = PME ; Treatment of long range electrostatic interactions
rcoulomb = 1.2 ; Short-range electrostatic cut-off
rvdw = 1.2 ; Short-range Van der Waals cut-off
pbc = xyz ; Periodic Boundary Conditions
periodic-molecules = yes

I got the following:

Using 1 MPI process
Using 56 OpenMP threads

Back Off! I just backed up em.trr to ./#em.trr.2#

Back Off! I just backed up em.edr to ./#em.edr.2#

Polak-Ribiere Conjugate Gradients:
Tolerance (Fmax) = 1.00000e+02
Number of steps = 1000
F-max = 1.56774e+06 on atom 3482
F-Norm = 1.99633e+04

Energy minimization has stopped, but the forces have not converged to the
requested precision Fmax < 100 (which may not be possible for your system).
It stopped because the algorithm tried to make a new step whose size was too
small, or there was no change in the energy since last step. Either way, we
regard the minimization as converged to within the available machine
precision, given your starting configuration and EM parameters.

Double precision normally gives you higher accuracy, but this is often not
needed for preparing to run molecular dynamics.

writing lowest energy coordinates.

Back Off! I just backed up em.gro to ./#em.gro.2#

Polak-Ribiere Conjugate Gradients converged to machine precision in 44 steps,
but did not reach the requested Fmax < 100.
Potential Energy = 2.8119632e+08
Maximum force = 3.3050301e+08 on atom 226084
Norm of force = 5.3002691e+05

After this minimization, i used the gro files (em.gro) to create an em.tpr file, i got the following result.

Back Off! I just backed up em.trr to ./#em.trr.3#

Back Off! I just backed up em.edr to ./#em.edr.3#

Polak-Ribiere Conjugate Gradients:
Tolerance (Fmax) = 1.00000e+02
Number of steps = 1000
F-max = 3.16473e+08 on atom 226084
F-Norm = 5.07357e+05

Energy minimization has stopped, but the forces have not converged to the
requested precision Fmax < 100 (which may not be possible for your system).
It stopped because the algorithm tried to make a new step whose size was too
small, or there was no change in the energy since last step. Either way, we
regard the minimization as converged to within the available machine
precision, given your starting configuration and EM parameters.

Double precision normally gives you higher accuracy, but this is often not
needed for preparing to run molecular dynamics.

writing lowest energy coordinates.

Back Off! I just backed up em.gro to ./#em.gro.3#

Polak-Ribiere Conjugate Gradients converged to machine precision in 34 steps,
but did not reach the requested Fmax < 100.
Potential Energy = 1.6851430e+08
Maximum force = 1.5776958e+08 on atom 246673
Norm of force = 2.5440825e+05

If I plan to carry out a long time(about 200ns MD) simulation, can I make do with a minimization output that gives a -ve PE but a Fmax > emtol. Also, what’s the basis for selecting the emtol. Can one have an emtol > 1,000?

8

You should investigate what is going on in the vicinity of the atoms with huge forces.

There are no hard and fast rules here for standard MD simulations, but you need to make sure you’re in a plausible starting configuration. The value of emtol is the gradient; if it is set to zero (impossible to achieve in practice) then you are exactly at the energy minimum. Values of 100 - 1000 are common and somewhat permissive in terms of their position on the potential energy surface, but are reasonable for getting a stable simulation started.

1 Like
9

Thank you very much, Dr Justin.

10

Hi Dr Jalemkul

I am very new to GROMACS and MD simulations. I am trying to run both SD and CG energy minimization for my system. How do I do that? Commands will be greatly appreciated! Thank you