[gmx-users] System blowing up in final MD run

Justin Lemkul jalemkul at vt.edu
Wed Jul 13 15:01:29 CEST 2016



On 7/13/16 8:56 AM, Deep Bhattacharya wrote:
> Hello Tsjerk,
>
> These are the notes that I got before running the NPT equilibration
>
>
> NOTE 1 [file npt.mdp]:
>   leapfrog does not yet support Nose-Hoover chains, nhchainlength reset to 1
>
> Generated 165 of the 1596 non-bonded parameter combinations
> Excluding 3 bonded neighbours molecule type 'Protein_chain_A'
> Excluding 3 bonded neighbours molecule type 'UNK'
> Excluding 2 bonded neighbours molecule type 'SOL'
> Excluding 1 bonded neighbours molecule type 'NA'
> The center of mass of the position restraint coord's is  3.516  4.410  3.815
>
> NOTE 2 [file topol.top]:
>   The largest charge group contains 11 atoms.
>   Since atoms only see each other when the centers of geometry of the charge
>   groups they belong to are within the cut-off distance, too large charge
>   groups can lead to serious cut-off artifacts.
>   For efficiency and accuracy, charge group should consist of a few atoms.
>   For all-atom force fields use: CH3, CH2, CH, NH2, NH, OH, CO2, CO, etc.
>
>
> NOTE 3 [file topol.top]:
>   The bond in molecule-type Protein_chain_A between atoms 1 N and 2 H1 has
>   an estimated oscillational period of 1.0e-02 ps, which is less than 10
>   times the time step of 1.0e-03 ps.
>   Maybe you forgot to change the constraints mdp option.
>
> Number of degrees of freedom in T-Coupling group Protein_UNK is 5339.83
> Number of degrees of freedom in T-Coupling group Water_and_ions is 89787.17
>
> NOTE 4 [file npt.mdp]:
>   You are using a plain Coulomb cut-off, which might produce artifacts.
>   You might want to consider using PME electrostatics.
>

Here's your biggest problem.  Plain cutoffs are wildly inaccurate and haven't 
generally been used in about two decades.  Use something more accurate like PME.

> I am using GROMOS53a6 forcefield. The system is a protein carbohydrate
> solvent system. While observing the trajectory during the MD run, I found
> that the system was blowing up.
> Please help.
>
> Sincerely,
> Deep Bhattacharya
>
> On Wed, Jul 13, 2016 at 12:37 PM, Tsjerk Wassenaar <tsjerkw at gmail.com>
> wrote:
>
>> Hi Deep,
>>
>> What force field? What system? How many cores? Did you try running on one?
>> What error? Did pdb2gmx or grompp give any warnings?
>>
>> Cheers,
>>
>> Tsjerk
>>
>> On Wed, Jul 13, 2016 at 5:44 AM, Deep Bhattacharya <
>> hypergenetics at gmail.com>
>> wrote:
>>
>>> Hello,
>>> My system is blowing up for the protein carbohydrate solvent complex. I
>>> have attached the mdp file below. The NPT and NVT were in good accordance
>>> to tutorial listed by Justin. 'Please help me solve the issue.
>>> title       = CD44  complex MD simulation
>>> ; Run parameters
>>> integrator  = md        ; leap-frog integrator
>>> nsteps      = 10000000    ; 0.001 * 1000000 = 10000 ps (10 ns)
>>> dt          = 0.001     ; 1 fs
>>> ; Output control
>>> nstxout             = 50000         ; suppress .trr output
>>> nstvout             = 50000         ; suppress .trr output
>>> nstenergy           = 50000    ; save energies every 1000.0 ps
>>> nstlog              = 50000      ; update log file every 1000.0 ps
>>> nstxtcout = 50000
>>> xtc_precision = 1000
>>> energygrps          = Protein LIG SOL
>>> ; Bond parameters
>>> continuation    = yes           ; first dynamics run
>>> constraint_algorithm = lincs    ; holonomic constraints
>>> constraints     = all-bonds     ; all bonds (even heavy atom-H bonds)
>>> constrained
>>> lincs_iter      = 1             ; accuracy of LINCS
>>> lincs_order     = 4             ; also related to accuracy
>>> ; Neighborsearching
>>> ns_type         = grid      ; search neighboring grid cells
>>> nstlist         = 1        ;

nstlist should not be set to 1.  That is only for energy minimization and you 
will be losing a huge amount of performance if you do this during MD.

>>> rcoulomb        = 0.81
>>> rvdw            = 0.81

These values are incorrect.  rvdw should be 1.4 for GROMOS force fields, and 
rcoulomb should be 0.8 (used in the original parametrized, with RF for 
electrostatics) or 0.9 (I've seen this more recently).  With PME, rcoulomb 
becomes a bit more flexible.

>>> coulombtype     = Cut-off       ; Particle Mesh Ewald for long-
>>> vdwtype        = cut-off
>>> rlistlong      = 1.4
>>> epsilon-rf     = 61
>>> rlist          = 0.8

As with rvdw, rlist should be 1.4.

-Justin

>>> ; Temperature coupling is on
>>> tcoupl = nose-hoover            ; modified Berendsen thermostat
>>> tc-grps       =Protein_LIG Water_and_ions   ; two coupling groups - more
>>> accurate
>>> tau_t = 2.0  2.0        ; time constant, in ps
>>> ref_t = 300  300        ; reference temperature, one for each group, in K
>>> ; Pressure coupling is on
>>> pcoupl        = Parrinello-Rahman            ; Pressure coupling on in
>> NPT
>>> pcoupltype        = isotropic            ; uniform scaling of box vectors
>>> tau_p        = 6.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       ; 3-D PBC
>>> ; Dispersion correction
>>> DispCorr    = EnerPres  ; account for cut-off vdW scheme
>>> ; Velocity generation
>>> gen_vel     = no        ; assign velocities from Maxwell distribution
>>>
>>>
>>> Steepest Descents converged to Fmax < 1000 in 573 steps
>>> Potential Energy  = -7.58370311664555e+05
>>> Maximum force     =  9.53896348096555e+02 on atom 1163
>>> Norm of force     =  2.44608995242254e+01
>>>
>>>
>>> *Deep S Bhattacharya*
>>>
>>> *Graduate Research Assistant*
>>>
>>> Mohs Biomedical Imaging & Nanotechnology Group
>>>
>>> Pharmaceutical Sciences
>>>
>>> *University of Nebraska Medical Center*
>>>
>>> 4018 Eppley Science Hall  |  Omaha, NE 68198-6805
>>>
>>> office 402.559.4349  | cell 402.906.1640
>>>
>>> deep.bhattacharya at unmc.edu.deep.bhattacharya199@gmail.com
>>> --
>>> Gromacs Users mailing list
>>>
>>> * Please search the archive at
>>> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
>>> posting!
>>>
>>> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>>>
>>> * For (un)subscribe requests visit
>>> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
>>> send a mail to gmx-users-request at gromacs.org.
>>>
>>
>>
>>
>> --
>> Tsjerk A. Wassenaar, Ph.D.
>> --
>> Gromacs Users mailing list
>>
>> * Please search the archive at
>> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
>> posting!
>>
>> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>>
>> * For (un)subscribe requests visit
>> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
>> send a mail to gmx-users-request at gromacs.org.
>>
>
>
>

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

Justin A. Lemkul, Ph.D.
Ruth L. Kirschstein NRSA Postdoctoral Fellow

Department of Pharmaceutical Sciences
School of Pharmacy
Health Sciences Facility II, Room 629
University of Maryland, Baltimore
20 Penn St.
Baltimore, MD 21201

jalemkul at outerbanks.umaryland.edu | (410) 706-7441
http://mackerell.umaryland.edu/~jalemkul

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


More information about the gromacs.org_gmx-users mailing list