[gmx-users] compressing the system
Moeed
lecielll at googlemail.com
Thu Jul 29 05:41:05 CEST 2010
Dear Dr.Chaban,
I tried to get the same results as you got for density and total energy
starting from the 8 chains is the box. I used the revised mdp file and since
PME option takes much time I used shift instead for a moment to verify if I
can get the density you obtained. Also I used constraints for bonds. density
I get after 1800ps is 850 kg/m3 while what you got is a slightly different
(after 1000ps density varies around 880!). Also the total energy you have is
bout 2300 whereas I get somethong around 12000!. I dont know why they are so
different..
In the second trial I used "constraints = none" just the "shift"
for electrostatics is not the same as the mdp file you suggested. again
results are not the same for total energy (20000 Kj/mol). Below is the
average values...also simulation crashes at 1550 ps and the last density I
am getting is 834! still less than what you got. dont really know why
results are not the same. The reason simulation crashes could be due to
dt=0.002 while removing constraints?
I have two more inquiries:
1- What I see is that the system is getting compressed to a size of around 3
3 3 nm is the first 20ps abruptly and in the remainder of time system
reaches equilibrium... This is very different approach than what I was
trying doing before. I used berendsen scheme for pressure coupling and box
reduced is size slowly but I had difficulty reaching the size I want and
most of the simulations crashed. Does this mean berenden is not suitable for
compressing system for all systems?
2- You said I dont need NVT at all. But if I want to study the system I get
from NPT giving the density I want, can I not do NVT for a fixed system size
and calcualte energies or other quantities I am after?
Thanks for your attention
Moeed
Energy Average RMSD Fluct. Drift
Tot-Drift
-------------------------------------------------------------------------------
Bond 3885.23 141.365 139.725 0.0478427
74.3668
Angle 5864.23 149.886 149.885 -0.00114897
-1.78597
Ryckaert-Bell. 2243.62 144.068 134.737 -0.113664
-176.679
LJ-14 1167.5 30.8742 30.6402 -0.00845667
-13.1451
Coulomb-14 -553.954 31.3413 30.7865 -0.0130831
-20.3364
LJ (SR) -4898.73 405.434 386.584 -0.272304
-423.27
Coulomb (SR) 1250.22 56.2361 56.1252 0.0078672
12.2288
Potential 8958.11 379.038 344.366 -0.352947
-548.621
Kinetic En. 10834.2 156.085 156.066 -0.00535835
-8.32903
Total Energy 19792.3 414.016 381.523 -0.358305
-556.95
Temperature 300.069 4.32299 4.32248 -0.000148407
-0.230684
Pressure (bar) 29.5896 1263.17 1263.16 0.0112876
17.5454
Box-X 3.67408 1.72247 1.70558 -0.000536305
-0.833634
Box-Y 3.44445 1.61482 1.59898 -0.000502786
-0.781532
Box-Z 2.45126 1.66735 1.64995 -0.000535423
-0.832263
Density (SI) 815.504 72.8491 70.7123 0.0390324
60.6721
pV 54.9014 2111.1 2111.07 -0.0208174
-32.3586
Vir-XX 3606.98 1466.95 1466.92 0.0213665
33.2121
Box-Vel-XX -0.0182685 0.234428 0.232397 6.8633e-05
0.106683
Mu-Y 0.102393 0.966426 0.966077 -5.79406e-05
-0.0900631
Heat Capacity Cv: 12.4757 J/mol K (factor = 0.000207551)
; 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 = 900000 ; 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 = 0
nstfout = 0
nstlog = 1000 ; frequency to write energies to
log file
nstxtcout = 0 ; 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 ;PME rev ; tells gromacs how to model
electrostatics. Shift: 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.9 ;0 ; where to start switching the
Coulomb potential
rvdw-switch = 0.9 ;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.1 ;1.0 ; distance for coulomb cut-off
rvdw = 1.0 ; distance for coulomb cut-off
; Temperature coupling
Tcoupl = v-rescale ;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 = Parrinello-Rahman ;berendsen
Pcoupltype = semiisotropic ;isotropic ; isotropic: means the box
expands and contracts in all directions(x,y,z)in order to maintain the
proper pressure; semiisotropic: isotropic in x & y directions
tau_p = 1 1 ;0.5 ; time constant for
coupling in ps
compressibility = 4.5e-5 4.5e-5 ; compressibility of
solvent used in simulation in 1/bar
ref_p = 30.0 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 = none ;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/20100728/2959980d/attachment.html>
More information about the gromacs.org_gmx-users
mailing list