[gmx-users] How to avoid the infinite potential energy in simulations of dimers?

Justin Lemkul jalemkul at vt.edu
Sat Feb 22 19:49:12 CET 2020



On 2/21/20 1:01 PM, Qing Liu wrote:
> Dear Justin Lemkul,
>
>   Thanks for your reply. I just downloaded the Cryo-EM structure of the dimer from RCSB, added missing residues using Modeller and built systems by the following commands:
>
>      
>
>       gmx pdb2gmx -f dimer.pdb -o dimer.gro -p dimer.top -water tip3p -ff amber99sb-ildn -ignh
>        gmx editconf -f dimer.gro -o vac-min-pbc.gro -bt cubic -d 1.5  -c
>        gmx solvate -cp vac-min-pbc.gro -cs spc216.gro -p dimer.top -o vac-min-pbc-solv.gro
>        gmx grompp -v -f sol-min.mdp -c vac-min-pbc-solv.gro -p dimer.top -o vac-min-pbc-solv.tpr
>        gmx genion -s vac-min-pbc-solv.tpr -o vac-min-pbc-solv-salt.gro -conc 0.15 -neutral -pname NA -nname CL  -p dimer.top
>        gmx grompp -f sol-min.mdp -c vac-min-pbc-solv-salt.gro -p dimer.top -o vac-min-pbc-solv-salt-min.tpr
>        gmx  mdrun -v -deffnm vac-min-pbc-solv-salt-min
>        gmx grompp -f pr-md.mdp -c vac-min-pbc-solv-salt-min.gro -p dimer.top -o pr-md.tpr  -r vac-min-pbc-solv-salt-min.gro
>        gmx mdrun -nb gpu -v -deffnm pr-md
>
>
> The contents of sol-min.mdp file are:
>
>
>        ; Define can be used to control processes
>        define          = -DFLEXIBLE
>        ; Parameters describing what to do, when to stop and what to save
>         integrator      = steep         ; Algorithm (steep = steepest descent minimization)
>         emtol           = 1.0           ; Stop minimization when the maximum force < 1.0 kJ/mol
>         nsteps          = 5000          ; Maximum number of (minimization) steps to perform
>         nstenergy       = 1             ; Write energies to disk every nstenergy steps
>         energygrps      = System        ; Which energy group(s) to write to disk
>        ; Parameters describing how to find the neighbors of each atom and how to calculate the interactions
>       ns_type         = grid          ; Method to determine neighbor list (simple, grid)
>       coulombtype     = cut-off       ; Treatment of long range electrostatic interactions
>       rcoulomb        = 1.0           ; long range electrostatic cut-off
>       rvdw            = 1.0           ; long range Van der Waals cut-off
>       constraints     = none          ; Bond types to replace by constraints
>       pbc             = xyz           ; Periodic Boundary Conditions (yes/no)
>       sc-coul         = yes
>
>
>
> The log file of energy minimization shows:
>        Steepest Descents converged to machine precision in 2108 steps,
>        but did not reach the requested Fmax < 1.
>        Potential Energy  = -2.0181694e+07
>        Maximum force     =  3.4778253e+02 on atom 114839
>         Norm of force     =  3.9966155e+01
>
>
> The contents of pr-md.mdp are:
>       title           = PR MD
>      define          =-DPOSRES
>      ;run control
>      integrator      =md
>      tinit           =0
>     dt              =0.002
>      nsteps          =500000
>      comm_mode       =Linear
>      nstcomm         =10
>     comm_grps       =System
>     ;output control
>     nstxout         =0
>    nstvout         =0
>    nstlog          =50000
>     nstcalcenergy   =1
>     nstenergy       =50000
>     nstxtcout       =50000
>    xtc_grps        =System
>    energygrps      =System
>   ;neighbor searching
>    nstlist         =10
>    ns_type         =grid
>    pbc             =xyz
>    rlist           =1.4
>   ;electrostatics
>    coulombtype     =PME
>    rcoulomb        =1.4
>    ;vdw
>    vdwtype         =Cut-off
>    rvdw            =1.4
>    dispCorr        =EnerPres
>    ;Ewald
>    fourierspacing  =0.1
>    pme_order       =4
>    ewald_rtol      =1e-5
>   ;temperature coupling
>    tcoupl          =v-rescale
>    tc_grps         =System
>    tau_t           =0.1
>    ref_t           =200
>    ;velocity generation
>    gen-vel         =no
>    ;bonds
>    constraints     =all-bonds
>    constraint_algorithm    =SHAKE
>    shake-tol       =0.0001
>    morse                   =no
>    continuation            =yes
>    sc-coul         =yes
>
>
>
> The pr-md.log shows:
>
>     Fatal error:
>     Step 2460: The total potential energy is nan, which is not finite. The LJ and
>     electrostatic contributions to the energy are 3.25168e+06 and -1.99067e+07,
>     respectively. A non-finite potential energy can be caused by overlapping
>     interactions in bonded interactions or very large or Nan coordinate values.
>     Usually this is caused by a badly- or non-equilibrated initial configuration,
>     incorrect interactions or parameters in the topology.
>
>

You're changing cutoffs, treatment of constraints, etc. at various 
points in your protocol. Don't do that. Be consistent. You probably have 
a water molecule that got stretched too far via -DFLEXIBLE and now it 
can't be constrained properly.

-Justin

>
>
>
>
>
> --
>
>
> Best wishes,
>
> ------------------------------------------------------------
>
> Qing Liu
>
> Fudan Univ.
>
> Mobile: +86—13358129621
>
> E-mail: 15110700076 at fudan.edu.cn
>
>
>
>
>
>
>> Message: 5
>> Date: Fri, 21 Feb 2020 09:21:34 -0500
>> From: Justin Lemkul <jalemkul at vt.edu>
>> To: gmx-users at gromacs.org
>> Subject: Re: [gmx-users] How to avoid the infinite potential energy in
>> 	simulations of dimers?
>> Message-ID: <b1dbc8bb-8640-5834-dd59-9a69a98e591a at vt.edu>
>> Content-Type: text/plain; charset=gbk; format=flowed
>>
>>
>>
>> On 2/21/20 5:11 AM, Qing Liu wrote:
>>> Dear Gromacs users,
>>>         I run some simulations of dimers, a protein binding another protein. When these simulations step into equilibriumstages,  they will crash with the following:
>>>
>>>
>>>
>>>
>>>
>>>                     Step 2029: The total potential energy is nan, which is not finite. The LJ and
>>>                     electrostatic contributions to the energy are 746119 and -5.71229e+06,
>>>                     respectively. A non-finite potential energy can be caused by overlapping
>>>                    interactions in bonded interactions or very large or Nan coordinate values.
>>>                   Usually this is caused by a badly- or non-equilibrated initial configuration,
>>>                    incorrect interactions or parameters in the topology.
>>>          
>>>
>>>         I try to turn on soft-core potential but it's not working. However, the simulations of monomers with same mdp files are normal.  How should I do to solve the problem?
>>>
>> Your initial configuration is unreasonable as you have a massive LJ
>> repulsion. How did you construct the coordinates of the system? Likely
>> energy minimization will have reported unreasonable forces, as well.
>>
>> -Justin
>>
>> -- 
>> ==================================================
>>
>> Justin A. Lemkul, Ph.D.
>> Assistant Professor
>> Office: 301 Fralin Hall
>> Lab: 303 Engel Hall
>>
>> Virginia Tech Department of Biochemistry
>> 340 West Campus Dr.
>> Blacksburg, VA 24061
>>
>> jalemkul at vt.edu | (540) 231-3129
>> http://www.thelemkullab.com
>>
>> ==================================================

-- 
==================================================

Justin A. Lemkul, Ph.D.
Assistant Professor
Office: 301 Fralin Hall
Lab: 303 Engel Hall

Virginia Tech Department of Biochemistry
340 West Campus Dr.
Blacksburg, VA 24061

jalemkul at vt.edu | (540) 231-3129
http://www.thelemkullab.com

==================================================



More information about the gromacs.org_gmx-users mailing list