[gmx-users] mdrun WARING and crash
L.Liu at utwente.nl
L.Liu at utwente.nl
Thu Feb 28 12:59:05 CET 2013
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?
"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.
Would my understanding stupid wrong? Please ask me for more details If I didn't post my problem clearly.
Thanks a billion.
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: Wednesday, February 27, 2013 4:12 PM
To: Discussion list for GROMACS users
Subject: Re: [gmx-users] mdrun WARING and crash
On 2/27/13 9:43 AM, L.Liu at utwente.nl wrote:
> Dear Justin,
>
> Thank you for your reply. I have some other questions,
>
> 1) Is it allowed to set nrexcl=0 in topology file in order to make every coarse-grained particles interact with each.
>
Are there bonded interactions? If so, such an exclusion setup doesn't make much
sense.
> 2) If it is allowed, then in [pairs] section, do I need to list all the pairs?
> For example, a linear chain with 10 beads:
> [ pairs ]
> 1 2
> 1 3
> 1 4
> 1 5
> ......
> 2 3
> 2 4
> ......
> 9 10
>
> or it can recognize the 1-5 1-6 or above automatically.
>
Well, with nrexcl=0, then all nonbonded interactions will be calculated, but
pairs are usually used to modify the default interactions. Please see the very
recent, extensive discussions in the list archive on this topic.
-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