[gmx-users] energy minimization - infinite force on atom

[Mark Abraham]

On 8/04/2010 4:23 PM, Evelyne Deplazes wrote:
> Hi gmx-users
>
> I am somehow famliar with gromacs but only used it with coarse grained
> force fields (Martini) so far, not with all atom ones. I am now trying
> to run a simulation on a membrane protein system (protein embedded in a
> POPC bilayer solvated with water) using the gromos force field (Gromacs
> 4.0.2  ffG53a). Using gromacs tutorials I managed to create the required
> input files (.gro, .top and .itp). I at at the stage of running a energy
> minimization using position restraints on the protein and the peptide
> molecules. I can grompp the input files
>
> grompp -f min_posre.mdp -p topol.top -c system.gro -o topol.tpr -n index.ndx
>
> with the following input file:
>
> ------------------------------
> -------------------------------------------------------------------------------------
> define          = -DPOSRES -DPOSRES_LIPID
> integrator      = steep         ; Algorithm (steep = steepest descent
> minimization)
> emtol           = 100.0         ; Stop minimization when the maximum
> force < 1000.0 kJ/mol/nm
> emstep          = 0.01          ; Energy step size
> nsteps          = 50000         ; Maximum number of (minimization) steps
> to perform
>
> ; Parameters describing how to find the neighbors of each atom and how
> to calculate the interactions
> nstlist         = 1             ; 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)
> 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 (yes/no)
> -------------------------------------------------------------------------------------------------------------------
>
> and the following topology file
>
> -------------------------------------------------------------------------------------------------------------------
> ; Include forcefield parameters
> #include "ffG53a6_lipid.itp"
>
> ; Include chain topologies
> #include "protein.itp"
>
> ; Include Position restraint file
> #ifdef POSRES
> #include "posre_protein.itp"
> #endif
>
> ; Include POPC topology
> #include "popc.itp"
>
> #ifdef POSRES_LIPID
> #include "lipid_posre.itp"
> #endif
>
> ; Include water topology
> #include "spc.itp"
>
> ; Include generic topology for ions
> #include "ions.itp"
>
> [ system ]
> ; Name
> mscl protein with POPC bilayer
>
> [ molecules ]
> ; Compound        #mols
> Protein             1
> POPC                242
> SOL                 32145
> CL-                 91
> NA+                 91
> -------------------------------------------------------------------------------------------------------------------
>
> When I run the simulation I get the following output after 16 steps of
> the minimisation
> -------------------------------------------------------------------------------------------------------------------
> Stepsize too small, or no change in energy.
> Converged to machine precision,
> but not to the requested precision Fmax < 100
>
> Double precision normally gives you higher accuracy.
>
> Steepest Descents converged to machine precision in 16 steps,
> but did not reach the requested Fmax < 100.
> Potential Energy  =  8.0175742e+17
> Maximum force     =            inf on atom 4443
> Norm of force     =            inf
> -------------------------------------------------------------------------------------------------------------------
>
> It seems there is an infinite force on atom 4443 but I don't really
> understand what that means. The confout.gro structure looks ok ie the
> system has not "exploded". Should the protein not be restrained? How do
> I know my position restraints actually work? I tried increasing the
> force constant of the restraints or removing the restraints but still
> get the same output. I also tried adding -DFLEXIBLE since I am using the
> spc model, but still the same output.
>
> can anyone tell me what that message actually means and what is wrong
> with my system?

A useful diagnostic will be to remove the attempts to use position 
restraints in the "define = " line. If you still get an error like this, 
then your problem (massive positive energy, infinite forces) is almost 
certainly a severe atomic overlap, possibly near atom 4443. If so, get 
out a visualization program and look for water overlapping lipid, or 
something. Position restraints won't help resolve gross errors like that...

Mark