[gmx-users] mdrun WARING and crash
L.Liu at utwente.nl
L.Liu at utwente.nl
Mon Mar 11 16:59:13 CET 2013
Hallo Justin,
Thank you for your comments.
Taking your suggestions, I set nrexcl=1, and comment the [pairs] section, because there is no special case of non-bonded interactions to declare, then we try to see what happens.
Now we minimize it by steep, then by cg, both of the processor are very quick done, because after around 7000 steps, the energy can not go further down any more.
Then we finish the mdrun, the energies output are like:
Step Time Lambda
10000 10.00000 0.00000
Energies (kJ/mol)
LJ (SR) Coulomb (SR) Potential Kinetic En. Total Energy
1.49468e+05 0.00000e+00 1.49468e+05 8.73890e+03 1.58207e+05
Temperature Pressure (bar)
4.38208e+02 9.93028e+04
Although this time no running error, we find the outputs are extremely wired, for example through
VMD conf.gro traj.xtc
we find 0 frame a homogeneous box, starting from the first step, the box becomes a lattice, which is far away from our expected the Polymer melt system should be.
The force parameters are taken from literature PRL 85(5), 1128(2000), I am still very worried about the format of my input files. Could you please give me, a very beginner a help.
Thanks a lot.
Kind regards,
Li
________________________________________
From: gmx-users-bounces at gromacs.org [gmx-users-bounces at gromacs.org] on behalf of Justin Lemkul [jalemkul at vt.edu]
Sent: Thursday, February 28, 2013 3:02 PM
To: Discussion list for GROMACS users
Subject: Re: [gmx-users] mdrun WARING and crash
On 2/28/13 6:59 AM, L.Liu at utwente.nl wrote:
> Hallo Justin,
>
> Thank you for you help. I have read the previous discussions on this topic, which is very helpful.
> The link is: http://gromacs.5086.n6.nabble.com/What-is-the-purpose-of-the-pairs-section-td5003598.html
> Well, there are still something I want to make sure, which might be the reason of mdrun crash of my system.
>
> ###################Introduction of system##############################
> Linear Polyethylene melt: each chain contains 16 beads, each bead coarse grained 3 monomers. Number of chain in the box is 64.
>
> ####################Force Field##########################
> ffbonded.itp
> [ bondtypes ]
> ; FENE, Ro = 1.5 sigma and kb = 30 epsilon/sigma^2
> ; i j func b0 (nm) kb (kJ/mol nm^2)
> CH2 CH2 7 0.795 393.
>
> ffnonbonded.itp
> [ atomtypes ]
> ; epsilon / kB = 443K
> ;name at.num mass (au) charge ptype sigma (nm) epsilon (kJ/mol)
> CH2 1 42.30000 0.000 A 0.5300 3.68133
>
> [ nonbond_params ]
> ; i j func sigma epsilon
> CH2 CH2 1 0.530 3.68133
>
> [ pairtypes ]
> ; i j func sigma epsilon
> CH2 CH2 1 0.53 3.68133
>
> ############################topology######################
> [ defaults ]
> ; nbfunc comb-rule gen-pairs fudgeLJ fudgeQQ
> 1 2 no 1.0 1.0
>
> ; The force field files to be included
> #include "../forcefield/forcefield.itp"
>
> [ moleculetype ]
> ; name nrexcl
> PE 0
> [atoms]
> ; nr type resnr residu atom cgnr charge
> 1 CH2 1 PE C 1 0.0
> 2 CH2 1 PE C 2 0.0
> 3 CH2 1 PE C 3 0.0
> 4 CH2 1 PE C 4 0.0
> ......
> 15 CH2 1 PE C 15 0.0
> 16 CH2 1 PE C 16 0.0
>
> [ bonds ]
> ; ai aj funct c0 c1
> 1 2 7 0.795 393.
> 2 3 7 0.795 393.
> 3 4 7 0.795 393.
> ......
> 14 15 7 0.795 393.
> 15 16 7 0.795 393.
>
> [ pairs ]
> ; ai aj funct c0 c1
> 1 2 1 0.53 3.68133
> 1 3 1 0.53 3.68133
> 1 4 1 0.53 3.68133
> 2 3 1 0.53 3.68133
> 2 4 1 0.53 3.68133
> 2 5 1 0.53 3.68133
> ......
> 14 15 1 0.53 3.68133
> 14 16 1 0.53 3.68133
> 15 16 1 0.53 3.68133
>
> ##############################mdp###########################
>
> ;directories to include in your topology. Format
> include = -I/home/otterw/Install/gromacs-4.6/src/gmxlib
> integrator = md
> ;****************************************
> ; RUN CONTROL *
> ;****************************************
> tinit = 0 ; start time to run
> dt = 1e-3 ; 0.014
> nsteps = 10000000
> init-step = 0 ; start time step, for step i, t = tinit + dt*(init-step + i)
> ;COM restriction: Linear, Angular or None
> comm-mode = Linear ; Linear: remove COM translation
> ; Angular: remove COM translation and rotation
> ; None: no restriction on COM
> nstcomm = 1000
>
> ;****************************************
> ; OUTPUT CONTROL *
> ;****************************************
> nstxout = 100 ;frequency to write coordinates to output trajectory file
> nstvout = 1000 ;frequency to write velocities to output trajectory
> nstfout = 1000 ;frequency to write forces to output trajectory
> nstlog = 1000 ;frequency to write energies to log file
> nstcalcenergy = 1000 ;frequency for calculating the energies
> nstenergy = 1000 ;frequency to write energies to energy file
> nstxtcout = 100 ;frequency to write coordinates to xtc trajectory
> xtc-precision = 1.0 ;precision to write to xtc trajectory
> xtc-grps = PE ;grps: group(s) to write, use default: the whole system
> energygrps = PE ;group(s) to write to energy file
>
> ;****************************************
> ; NEIGHBOR SEARCHING *
> ;****************************************
> cutoff-scheme = Verlet ;group or verlet(generate a pair list with buffering)
> verlet-buffer-drift = -1 ;sets the target energy drift per particle caused by the Verlet buffer
> ; set to -1, rlist will be used
> rlist = 5.95e-1 ;Cut-off distance for the short-range neighbor list
> nstlist = 10 ;>0 Frequency to update the neighbor list
> ;=0 The neighbor list is only constructed once and never updated
> ;-1 Automated update frequency, only supported with cutoff-scheme=group
> ns-type = grid ;grid, check atoms in neighboring grid cells, when every nstlist steps
> ;simple, check every atom in the box, when every nstlist steps.
> pbc = xyz ;xyz,Use periodic boundary conditions in all directions
> ;no
> ;xy
> periodic-molecules = yes ;no, molecules are finite
> ;yes, molecules that couple to themselves through PBC
>
> ;****************************************
> ; Electrostatics *
> ;****************************************
> coulombtype = cut-off ; Twin range cut-offs with neighborlist cut-off rlist and Coulomb cut-off rcoulomb, where rcoulomb≥rlist.
> ; Shift; Switch; User... check VwD options below.
> coulomb-modifier = Potential-shift-Verlet
> rcoulomb = 5.95e-1 ; distance for the Coulomb cut-off
>
> ;****************************************
> ; VdW *
> ;****************************************
> vdwtype = cut-off ;Twin range cut-offs with neighbor list cut-off rlist and VdW cut-off rvdw, where rvdw ≥ rlist
> ;Shift:The LJ (not Buckingham) potential is decreased over the whole range and the forces decay smoothly to zero between rvdw-switch and rvdw, rlist>rvdw(0.1-0.3nm)
> ;Switch:The LJ (not Buckingham) potential is normal out to rvdw-switch, after which it is switched off to reach zero at rvdw,rlist>rvdw(0.1-0.3nm)
> ; User
> vdw-modifier = potential-shift-Verlet ;Potential-shift:Shift the Van der Waals potential by a constant such that it is zero at the cut-off
> ;None:Use an unmodified Van der Waals potential
> ;Potential-shift-Verlet:Selects Potential-shift with the Verlet cutoff-scheme
> rvdw = 5.95e-1 ;distance for the LJ or Buckingham cut-off
> Dispcorr = no ;don't apply any correction
> ;EnerPres:apply long range dispersion corrections for Energy and Pressure
> ;Ener:apply long range dispersion corrections for Energy only
>
> ;****************************************
> ; Temperature coupling *
> ;****************************************
> tcoupl = berendsen ; berendsen: ref-t [K],tau-t [ps],tc-grps
> ; nose-hoover
> ; no
> nsttcouple = -1 ; frequency for coupling the temperature
> ; default -1, nsttcouple=nstlist, unless nstlist≤0, then 10 is used.
> tc-grps = PE ; groups to couple separately to temperature bath
> tau-t = 2.0 ; time constant for coupling, -1 means no temperature coupling
> ref-t = 443 ; reference temperature for coupling
>
> ;****************************************
> ; Bonds *
> ;****************************************
> ;constraints = none ; all-bonds
> ;constraint-algorithm ; LINCS;SHAKE
>
> ;****************************************
> ; NEMD *
> ;****************************************
> ;acc-grps
> ;accelerate
> ;freezegrps
> ;freezedim
> ;cos-acceleration
> ;deform
>
> ###############Questions and Suspicions ##########################################
>
> 1. command grompp, the screen writes the following information. I am very confused by
> "Generated 0 of the 1 non-bonded parameter combinations", what does this mean? where are the numbers 0 and 1 from?
>
You have one atomtype (CH2), which means there is only type of nonbonded
interaction (CH2-CH2). Since you have specified parameters in [pairtypes], then
grompp does not need to generate the interaction. If there was no specification
in [pairtypes] then grompp would apply normal combination rules and fudge
factors in calculating the parameters for the interaction. All it's telling you
is that it found 1 interaction, and didn't have to calculate any of them.
> "NOTE 1 [file grompp.mdp]:
> The Berendsen thermostat does not generate the correct kinetic energy
> distribution. You might want to consider using the V-rescale thermostat.
>
> Generated 0 of the 1 non-bonded parameter combinations
> Excluding 0 bonded neighbours molecule type 'PE'
> Removing all charge groups because cutoff-scheme=Verlet
> Analysing residue names:
> There are: 64 Other residues
> Analysing residues not classified as Protein/DNA/RNA/Water and splitting into groups...
> Number of degrees of freedom in T-Coupling group PE is 3069.00
> Determining Verlet buffer for an energy drift of 0.005 kJ/mol/ps at 443 K
> Calculated rlist for 1x1 atom pair-list as 0.595 nm, buffer size 0.000 nm
> Set rlist, assuming 4x4 atom pair-list, to 0.595 nm, buffer size 0.000 nm"
>
> 2. command mdrun, it runs a while then complains with
>
> "WARNING: Listed nonbonded interaction between particles 140 and 143
> at distance 3f which is larger than the table limit 3f nm.
>
> This is likely either a 1,4 interaction, or a listed interaction inside
> a smaller molecule you are decoupling during a free energy calculation.
> Since interactions at distances beyond the table cannot be computed,
> they are skipped until they are inside the table limit again. You will
> only see this message once, even if it occurs for several interactions.
>
> IMPORTANT: This should not happen in a stable simulation, so there is
> probably something wrong with your system. Only change the table-extension
> distance in the mdp file if you are really sure that is the reason."
>
> At last it crashes because of very large bonds, and from md.log it can be seen that the LJ(SR) goes extremely large, so does the temperature.
> Following your suggestions on your previous emails, I checked the Blowing up reasons.
> I smaller the time step and make nstlist=1 to update the neighbourlist every step, run energy minimization and fully done EQ run.
> The mdrun does survive after all these, but with some strange output files:
> dd_dump_err_0_n0.pdb
> dd_dump_err_0_n1.pdb
> .......
> dd_dump_err_0_n7.pdb
>
> and if I check the RDF, by command 'g_rdf -f traj.xtc -n index.ndx', I find discrete peaks, moreover, one peak is even negative. These are of course not acceptable.
> I am very afraid that something wrong in my coarse grained system. Could you Please help me to check my input files for simulation? To myself there is a suspicious of [paris] listed on topology file,
> that according to what I understood from earlier discussions, if set nrexcl=0, then only 1-2, 1-3, 1-4 nonbonded interactions need to be listed, more beyond are calculated automatically.
>
Pair interactions and nonbonded interactions are different. In Gromacs, a pair
is actually considered a bonded interaction, while everything else beyond nrexcl
bonds is assigned to the normal nonbonded interaction list. The explanation in
the above thread from Bogdan is about as good of an explanation as you will
find, so I won't repeat it.
What's odd about what you're doing is you are calculating both bonded and
nonbonded forces for particles connected via bonds. This is not how most normal
MM force fields work, even CG ones. I'm assuming that LJ(SR) is a very large,
positive number, indicating repulsion among your particles. That, in turn, is
causing the bonds to become excessively strained and your system blows up. For
instance, the MARTINI CG force field excludes one bonded neighbor. I suppose
you can parameterize anything to work if you balance the forces, but it's not a
trivial task.
-Justin
--
========================================
Justin A. Lemkul, Ph.D.
Research Scientist
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin
========================================
--
gmx-users mailing list gmx-users at gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
* Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
* Please don't post (un)subscribe requests to the list. Use the
www interface or send it to gmx-users-request at gromacs.org.
* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
More information about the gromacs.org_gmx-users
mailing list