[gmx-users] NPT for compressing the system
Moeed
lecielll at googlemail.com
Wed Jul 14 20:05:01 CEST 2010
Hello,
Thanks for your reply. As for the simulation box, I did not change anything
and I am still saying that I am having a cubic box of size 30nm. When I view
the molecule in ngmx I see the cubic box.. no problem with that. What I was
referring to was the option on menu bar where one can select type of box to
view the trajecotry and using say triclinic box I see chain fits better
(almost completely)in the box while in the case of cubic (which is my real
box) I see molecules is not fitted completely inside the box.. I was just
asking why is this happening. as you said "none of that probably means
anything" . I think also it is unnecessary for now so please dont get
confused.
As you requested I am including ALL the commands I used from the very
beginning starting from replicating a 4C unit.
editconf -f Eth4C.pdb -o Eth4C-d.gro -princ -d 0.058 -bt cubic
editconf -f Eth4C-d.gro -o Eth4C-d-index.gro -n index.ndx -princ
genconf -f Eth4C-d-index.gro -o Eth4C-d-index-rep30.gro -nbox 30 1 1
now I have PE chain of 60 units.
using proper rtp and hdb files:
pdb2gmx -f
Eth4C-d-index-rep30-res.gro -o Eth4C-d-index-rep30-res-pdb2.gro -p
Eth4C-d-index-rep30-res-pdb2.top -ff oplsaa >& output.pdb2gmx
Then I EM this last gro file but since dimensions are 15 * 0.5 0.5 nm
I am getting the error message :
The cut-off length is
longer than half the shortest box vector or longer than the smallest
box diagonal element. Increase the box size or decrease rlist.
then increased the size >
editconf -f Eth4C-d-index-rep30-res-pdb2.gro -o
Eth4C-d-index-rep30-res-pdb2-boxsize30.gro -box 30 30 30 -angles 90 90
90
Then I EM the molecule with new box size 30nm and call this energy
minimized structure PE60.gro. Now I want to replicate this chain.
density demands 1 molecule in 5.76 nm3. As you said I forgot the
density for a moment and took a large box.
PE60 is already in a large enough box (30nm). this is all I did before
doing any NPT.
As for compressing the system to get more realistic structure and
approaching the density I am after, I did several NPT runs. All show
system is not equilbrated and they all crash and box does not become
smaller than 2.2 nm.
As you proposed I tried doing runs in stepwise manner. Below you see my trails:
First trial :P 120 bar for 100 ps then P80 bar for 500 ps then P50 for
500 ps. this last run crashes at 2.2 nm size
As Mark suggested I tried started with lower pressure:
Second trail: P 50 for 500ps > P40 bar for 500 ps > P35 for 1000ps >
P30 for 1000 ps . Below is the results of this last run. before I had
negative average Pressure but below the pressure is 30bar and real
values fluctuate between -500 to 500 bar.
T thermostat is working fine but LJ Columb terms not reaching
equilibrium. Also box size is 2.7 nm which is larger than former
trials.Molecule in not inside the box for almost the entire
simulation.Only a small portion in inside and since it is long can not
fit in.
Do i need to start with lower pressures and go to higher values.
Because the size is not changing that much with 30 bar?
Energy Average RMSD Fluct. Drift Tot-Drift
-------------------------------------------------------------------------------
Angle 736.582 30.0543 30.0543 -9.65991e-05
-0.0965993
Ryckaert-Bell. 147.603 13.5226 13.5185 -0.00115216 -1.15216
LJ-14 153.926 6.02125 6.02079 0.000258266 0.258266
Coulomb-14 -106.141 1.90986 1.90613 0.00041316 0.413161
LJ (SR) -354.892 15.5047 15.4964 -0.00176674 -1.76675
Coulomb (SR) 217.529 2.63621 2.6302 -0.000615852
-0.615853
Potential 794.607 37.4862 37.4764 -0.00295993 -2.95994
Kinetic En. 899.88 30.6454 30.6399 0.00201524 2.01525
Total Energy 1694.49 18.7071 18.7051 -0.000944686
-0.944688
Temperature 299.806 10.2099 10.2081 0.000671404 0.671405
Pressure (bar) 30.0416 269.98 269.948 0.0145652 14.5652
Cons. rmsd () 5.11735e-06 3.49122e-07 3.49102e-07
1.29073e-11 1.29073e-08
Box-X 2.79128 0.00907586 0.00906645 1.43079e-06 0.00143079
Box-Y 2.79128 0.00907586 0.00906645 1.43079e-06 0.00143079
Box-Z 2.79128 0.00907586 0.00906645 1.43079e-06 0.00143079
Density (SI) 128.686 1.2601 1.25879 -0.000199256
-0.199256
pV 39.3202 353.469 353.428 0.0188014 18.8014
Vir-XX 327.606 298.203 298.201 -0.00356238 -3.56239
Mu-X -0.028926 0.304356 0.30434 -1.08836e-05
-0.0108837
Lamb-System 1.00002 0.000368363 0.000368248
-3.18082e-08 -3.18083e-05
Heat Capacity Cv: 12.4935 J/mol K (factor = 0.00115974)
***********************************************************************************
itle = PE-Hexane
;define = -DPOSRES ; tells gromacs to perform position restrained
dynamics/include posre.itp into topology used for position restraint
pbc = xyz ; use priodic BCs in all directions
; Run control
integrator = md ; type of dynamics algorithm. Here md uses
a leap-frog algorithm for integrating Newtons's eq of motion
dt = 0.002 ; in ps !
nsteps = 500000 ; length of simulation= nsteps*dt
nstcomm = 1 ; frequency for center of mass motion removal
; Output control
nstenergy = 100 ; frequency to write energies to energy
file. i.e., energies and other statistical data are stored every 10
steps
nstxout = 100 ; frequency to write
coordinates/velocity/force to output trajectory file. how often
snapshots are collected= nstxout*dt
nstvout = 100
nstfout = 0
nstlog = 100 ; frequency to write energies to log file
nstxtcout = 10 ; frequency to write coordinates to xtc trajectory
; Neighbor searching
nstlist = 10 ; frequency to update neighbor list.
Neighborlist will be updated at least every 10 steps. Manual p80
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
; Electrostatics/VdW
coulombtype = Shift ; tells gromacs how to model
electrostatics. Coulomb/LJ potential is decreased over the whole range
and forces decay smoothly to zero between
vdw-type = Shift ; rcoulomb-switch/rvw-switch & rcoulomb/rvdw
rcoulomb-switch = 0 ; where to start switching the Coulomb potential
rvdw-switch = 0 ; where to start switching the LJ potential
; Cut-offs
rlist = 1.1 ; in nm. Cut-off distance for short-range
neighbor list
rcoulomb = 1.0 ; distance for coulomb cut-off
rvdw = 1.0 ; distance for coulomb cut-off
; Temperature coupling
Tcoupl = berendsen
tc-grps = System ;HEX ; groups to couple to
thermostat; Berendsen temperature coupling is on in these groups
tau_t = 0.1 ;0.1 ; time constant for T coupling
ref_t = 300 ;300 ; reference T for coupling. When
you alter the T, don't forget to change the gen_temp for velocity
generation
; Pressure coupling
Pcoupl = berendsen ; Pressure coupling is not on
Pcoupltype = isotropic ; means the box expands and contracts in
all directions (x,y,z) in order to maintain the proper pressure
tau_p = 0.5 ; time constant for coupling in ps
compressibility = 4.5e-5 ; compressibility of solvent used in
simulation in 1/bar
ref_p = 30.0 ; reference P for coupling in bar
; Velocity generation Generate velocites is on at 300 K. Manual p155
gen_vel = yes ; generate velocites according to Maxwell
distribution at T: gen_temp with random gen seed gen_seed
gen_temp = 300.0 ; T for Maxwell distribution
gen_seed = 173529 ; used to initialize random generator
for random velocities
; Bonds
constraints = all-bonds ; sets the LINCS constraint for all bonds
constraint-algorithm = lincs
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-users/attachments/20100714/78e28022/attachment.html>
More information about the gromacs.org_gmx-users
mailing list