[gmx-users] NPT for compressing the system

Moeed lecielll at googlemail.com
Wed Jul 14 20:05:01 CEST 2010


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

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

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
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
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
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
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
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
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

;		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