<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
<html>
<head>
<meta content="text/html; charset=ISO-8859-1"
http-equiv="Content-Type">
<title></title>
</head>
<body bgcolor="#ffffff" text="#000000">
On 14/07/2010 1:09 PM, Moeed wrote:
<blockquote
cite="mid:AANLkTinRRO4ycVXzzszRu3GdtauQQ2SK8jgJt7DDBu1P@mail.gmail.com"
type="cite">Hello,<br>
<br>
Thank you for suggestions. Yes simulation crashes. I am not able see
the whole run (1000ps). For all pressures I tired the final simulation
box is 2.2*2.2*2.2 nm but simulation crashes before 1000ps. The longest
simulation is for 120 bar (727ps). and I think that is why I am not
getting the _after_md.gro file after mdrun in none of the simulations.
I wanted to take that convoluted structure and replicate the box.. </blockquote>
<br>
I don't recall why you have a large box that you have to compress, but
the way to do that is to be gentler with it. What you're doing is
integrating equations of motion with forces that are much larger than
those normally experienced, and that is not a numerically stable
situation if you keep the timestep the same size, since the atoms will
travel so far that the approximations in the integration scheme are no
longer reasonable. Think about how (badly) a similar model of a simple
pendulum would work if the timestep was comparable with the period of
the motion, and someone was bashing the pendulum each time step...<br>
<br>
A few bars (maybe 10?) of overpressure will compress the box, and do so
in such a way that the system does not get too far from equilbrium.
That maximizes your chance of reaching your desired box size with a
system that is not greatly unhappy. Even if such an unhappy system was
not exploding, you'll still have to equilibrate such a system for
longer than one that had a more gentle compression regime. Write the
coordinates fairly often, and then even if you overshoot the
compression, you can get a structure of about the right size to start
again to equilibrate at about the right size with 1 bar P-coupling.<br>
<br>
There's no value in doing the compression fast if fast doesn't work!<br>
<br>
Also, generating velocities and then immediately trying to compress the
box before anything has equilibrated is asking for trouble. If your
system density is very far from equilibrium, you might need to use a
small timestep to get anything to work<br>
<br>
Mark<br>
<br>
<blockquote
cite="mid:AANLkTinRRO4ycVXzzszRu3GdtauQQ2SK8jgJt7DDBu1P@mail.gmail.com"
type="cite">
<meta http-equiv="CONTENT-TYPE"
content="text/html; charset=ISO-8859-1">
<title></title>
<meta name="GENERATOR" content="OpenOffice.org 3.1 (Unix)">
<style type="text/css">
<!--
@page { margin: 2cm }
PRE.western { font-family: "Nimbus Roman No9 L", "Arial Unicode MS", serif }
PRE.cjk { font-family: "DejaVu Sans Mono", "Arial Unicode MS", monospace }
PRE.ctl { font-family: "DejaVu Sans Mono", "Arial Unicode MS", monospace }
P { margin-bottom: 0.21cm }
-->
</style>
<pre class="western" style="text-align: justify;">mdrun -s PE60-P120_md.tpr -o PE60-P120_md.tpr -c PE60-P120_after_md -v >& output.mdrun_md</pre>
As you suggested I do the runs in stepwise manner. I need still a
little smaller box size than 2*2*2 but I noticed sth during these
simulations. When I view the trajectory in ngmx on the menu bar I can
select type of box. When I select triclinic or octahedron the final
confromation of the chain is more convoluted and fits completely in the
box while in case of rectangular box except for pressure of 120 bar for
all other pressures I tried most parts of the chain remain outside of
the small box. molecules is more extended...I am wondeirng if positions
of atoms are calculated and stored how changing the box type in
trajectory viewer can affect the structure?!<br>
***********************************************************************************************mdp
file I used:<br>
title = PE-Hexane
<br>
;define = -DPOSRES ; tells gromacs to perform position
restrained dynamics/include posre.itp into topology used for position
restraint<br>
pbc = xyz ; use priodic BCs in all directions<br>
<br>
; Run control<br>
integrator = md ; type of dynamics algorithm. Here md
uses a leap-frog algorithm for integrating Newtons's eq of motion<br>
dt = 0.002 ; in ps !<br>
nsteps = 500000 ; length of simulation=
nsteps*dt <br>
nstcomm = 1 ; frequency for center of mass motion
removal <br>
<br>
; Output control<br>
nstenergy = 100 ; frequency to write energies to
energy file. i.e., energies and other statistical data are stored every
10 steps<br>
nstxout = 100 ; frequency to write
coordinates/velocity/force to output trajectory file. how often
snapshots are collected= nstxout*dt<br>
nstvout = 100<br>
nstfout = 0<br>
nstlog = 100 ; frequency to write energies to log
file<br>
nstxtcout = 10 ; frequency to write coordinates to xtc
trajectory<br>
<br>
; Neighbor searching<br>
nstlist = 10 ; frequency to update neighbor list.
Neighborlist will be updated at least every 10 steps. Manual p80<br>
ns_type = grid ; make a grid in the box and only
check atoms in neighboring grid cells when constructing a new neighbor
list every nstlist steps<br>
<br>
; Electrostatics/VdW<br>
coulombtype = Shift ; tells gromacs how to model
electrostatics. Coulomb/LJ potential is decreased over the whole range
and forces decay smoothly to zero between <br>
vdw-type = Shift ; rcoulomb-switch/rvw-switch &
rcoulomb/rvdw<br>
rcoulomb-switch = 0 ; where to start switching the Coulomb
potential <br>
rvdw-switch = 0 ; where to start switching the LJ
potential <br>
<br>
; Cut-offs<br>
rlist = 1.1 ; in nm. Cut-off distance for
short-range neighbor list<br>
rcoulomb = 1.0 ; distance for coulomb cut-off <br>
rvdw = 1.0 ; distance for coulomb cut-off <br>
<br>
; Temperature coupling <br>
Tcoupl = berendsen<br>
tc-grps = System ;HEX ; groups to couple to
thermostat; Berendsen temperature coupling is on in these groups
<br>
tau_t = 0.1 ;0.1 ; time constant for T
coupling <br>
ref_t = 300 ;300 ; reference T for
coupling. When you alter the T, don't forget to change the gen_temp for
velocity generation<br>
<br>
; Pressure coupling<br>
Pcoupl = berendsen ; Pressure coupling is not on<br>
Pcoupltype = isotropic ; means the box expands and
contracts in all directions (x,y,z) in order to maintain the proper
pressure<br>
tau_p = 0.5 ; time constant for coupling in ps
<br>
compressibility = 4.5e-5 ; compressibility of solvent used
in simulation in 1/bar<br>
ref_p = 160.0 ; reference P for coupling in bar<br>
<br>
; Velocity generation Generate velocites is on at 300 K.
Manual p155<br>
gen_vel = yes ; generate velocites according to
Maxwell distribution at T: gen_temp with random gen seed gen_seed<br>
gen_temp = 300.0 ; T for Maxwell distribution <br>
gen_seed = 173529 ; used to initialize random
generator for random velocities<br>
<br>
; Bonds<br>
constraints = all-bonds ; sets the LINCS constraint for all
bonds<br>
constraint-algorithm = lincs<br>
*******************************************************************************************<br>
<br>
<br>
Thanks<br>
<br>
<br>
<br>
<br clear="all">
</blockquote>
<br>
</body>
</html>