<div class="gmail_quote"><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex;"><br>
I have been trying to use CHARMM forcefield to simulate 1000 molecules<br>
of Methanol in a box but I seem to fall short on the density every<br>
time I run. At first, I had my rlist at 1 but CHARMM recommends<br>
greater than 1.2-1.4 nm so I changed accordingly, then I changed it so<br>
that the short neighbor list would update more freqently and finally I<br>
even added a greater amount of steps, hoping that the pressure and<br>
density would settle down at some value.<br>
<br>
Below, I had taken and modified from the Lysozyme tutorial. One thing<br>
I noticed is the contraint -algorithm and how the constraints are set<br>
to 'all-bonds'. Is this necessary? After running NVT (leveled off<br>
near 298K, seemed ok), then NPT, for 300,000 steps, the pressure is<br>
close to the reference 1 bar, ~1.1,1.2 bar, but I noticed the density<br>
is still oscillating, not by much, but it still leaves me questioning<br>
why its not settling around to the literature value of 0.791 g/cm3. I<br>
posted the last several lines from the analysis *.xvg file. NPT<br>
script is also below.<br>
<br>
Is there still more steps that I have to run for? The density just<br>
seems to be oscillating too much and although a literature value of<br>
0.791 g/cm3 is within the data, it doesnt seem to be settling around<br>
that value. I have similar problems for 1-propanol to where my<br>
simulated density is under that of the literature value. Could it be<br>
possible that maybe I should use different NVT, NPT algorithms<br>
(tcoupl, pcoupl)?<br>
<br>
I would really appreciate anyones help. I'm still in my first few<br>
months of learning the program. Thanks!<br>
<br>
591.800000 38039.125000 297.409546 -206.071747 779.428894<br>
592.000000 38035.421875 294.951111 -20.557419 777.816956<br>
592.200000 37982.812500 298.063324 -463.192047 776.430542<br>
592.400000 37681.421875 295.430878 36.529854 776.460876<br>
592.600000 37476.148438 297.193787 -28.961288 775.634644<br>
592.800000 37611.574219 297.059875 -553.369263 774.096924<br>
593.000000 37682.277344 288.297974 -293.567963 776.055237<br>
593.200000 37147.824219 295.985901 -221.989441 781.570679<br>
593.400000 37404.906250 293.938721 399.093018 788.921570<br>
593.600000 37339.914062 297.464783 415.891235 793.838440<br>
593.800000 37909.019531 292.707184 574.818726 796.730286<br>
594.000000 37369.820312 304.095703 -119.569496 796.235962<br>
594.200000 37293.445312 294.449005 -63.335049 794.578552<br>
594.400000 37196.917969 295.323761 104.148239 791.275085<br>
594.600000 37776.199219 298.004333 28.989655 785.244629<br>
594.800000 37893.097656 302.397308 -268.065765 779.457886<br>
595.000000 38123.460938 295.343109 -140.914047 775.682678<br>
595.200000 38002.054688 302.233124 -91.893478 774.813660<br>
595.400000 37835.652344 297.808319 441.931885 776.142090<br>
595.600000 38103.964844 296.369019 490.766754 778.422302<br>
595.800000 37686.902344 303.358093 -19.726471 781.057922<br>
596.000000 37586.890625 293.751648 -92.387199 783.724854<br>
596.200000 37408.414062 300.007385 -288.156036 785.990417<br>
596.400000 37894.898438 295.720520 346.861206 788.753174<br>
596.600000 37673.378906 288.311432 89.059952 789.179993<br>
596.800000 36881.062500 295.265564 70.789131 791.552673<br>
597.000000 37077.789062 296.261444 72.207787 793.970398<br>
597.200000 37123.230469 296.511475 -232.434616 793.653564<br>
597.400000 37426.152344 304.281891 -133.300797 791.791077<br>
597.600000 37514.332031 305.183197 -259.364136 787.592712<br>
597.800000 37393.667969 301.326111 -150.162338 785.053711<br>
598.000000 37263.054688 299.811554 326.851837 785.553894<br>
598.200000 37203.261719 300.688904 225.526001 788.008667<br>
598.400000 37511.636719 302.506714 66.299194 790.725281<br>
598.600000 37460.843750 299.512909 198.943665 793.353333<br>
598.800000 37248.921875 291.575134 -93.283127 793.582703<br>
599.000000 36870.480469 293.515381 211.069107 792.929443<br>
599.200000 37398.097656 293.184448 -163.371140 789.150513<br>
599.400000 37402.238281 297.402008 45.650658 786.303711<br>
599.600000 37419.816406 298.520508 -147.688004 783.492371<br>
599.800000 37497.332031 296.529877 -15.991615 782.294250<br>
600.000000 37691.468750 298.839844 -68.580627 781.974182<br>
<br>
<br>
<br>
title =Methanol npt equilibration<br>
<br>
;Run parameters<br>
integrator =md ;leap-frog algorithm<br>
nsteps =300000 ;2 * 50000 = 100 ps<br>
dt =0.002 ;2fs<br>
<br>
;Output control<br>
nstxout =100 ;save coordinates every 0.2 ps<br>
nstvout =100 ;save velocities every 0.2 ps<br>
nstenergy =100 ;save energies every 0.2 ps<br>
nstlog =100 ;update log file every 0.2 ps<br>
<br>
;Bond parameters<br>
continuation =yes ;Restarting after NVT<br>
constraint_algorithm=lincs ;holonomic constraints<br>
constraints =all-bonds ;all bonds (even heavy atom-H<br>
bonds)constraind<br>
lincs_iter =1 ;accuracy of LINCS<br>
lincs_order =4 ;also related to accuracy<br>
<br>
;Neighborhood searching<br>
ns_type =grid ;search neighboring grid cells<br>
nstlist =1 ;2 fs<br>
rlist =1.4 ;short-range neighborlist cutoff (in nm)<br>
rcoulomb =1.4 ;short-range electrostatic cutoff (in nm)<br>
rvdw =1.4 ;short-range van der Waals cutoff (in nm)<br>
<br>
;Electrostatics<br>
coulombtype =PME ;Particle Mesh Ewald for long-range electrostat;ics<br>
pme_order =4 ;cubic interpolation<br>
fourierspacing =0.16 ;grid spacing for FFT<br>
<br>
;Temperature coupling is on<br>
tcoupl =V-rescale ;modified Berendsen thermostat<br>
tc-grps =SYSTEM ;one coupling group<br>
tau_t =0.1 ;time constant, in ps<br>
ref_t =298 ;reference temperature, one for<br>
each ;group, in K<br>
<br>
;Pressure coupling is on<br>
pcoupl =Parrinello-Rahman ;Pressure coupling on in NPT<br>
pcoupltype =isotropic ;uniform scaling of box vectors<br>
tau_p =2.0 ;time constant, in ps<br>
ref_p =1.0 ;reference pressure, in bar<br>
compressibility =4.5e-5 ;isothermal compressibility of h2O, 1/bar<br>
<br>
;Periodic boundary conditions<br>
pbc =xyz ; 3-D PBC<br>
<br>
;Dispersion correction<br>
DispCorr =EnerPres ;account for cut-off vdW scheme<br>
<br>
;Velocity generation<br>
gen_vel =no ;Velocity generation is off<br></blockquote><div><br></div><div><br></div><div><br></div><div>The problem should be with your topology. List your ITP files.</div><div> </div></div>-- <br>Dr. Vitaly V. Chaban, Department of Chemistry<br>
University of Rochester, Rochester, New York 14627-0216<br>