On 2/28/13 6:59 AM, [email protected] 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 [email protected]
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 [email protected].
* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
