[gmx-users] Unsteady Density

Fabian Casteblanco fabian.casteblanco at gmail.com
Mon May 2 20:16:46 CEST 2011


Hello,

I have been trying to use CHARMM forcefield to simulate 1000 molecules
of Methanol in a box but I seem to fall short on the density every
time I run.  At first, I had my rlist at 1 but CHARMM recommends
greater than 1.2-1.4 nm so I changed accordingly, then I changed it so
that the short neighbor list would update more freqently and finally I
even added a greater amount of steps, hoping that the pressure and
density would settle down at some value.

Below, I had taken and modified from the Lysozyme tutorial.  One thing
I noticed is the contraint -algorithm and how the constraints are set
to 'all-bonds'.  Is this necessary?  After running NVT (leveled off
near 298K, seemed ok), then NPT, for 300,000 steps, the pressure is
close to the reference 1 bar, ~1.1,1.2 bar, but I noticed the density
is still oscillating, not by much, but it still leaves me questioning
why its not settling around to the literature value of 0.791 g/cm3.  I
posted the last several lines from the analysis *.xvg file.  NPT
script is also below.

Is there still more steps that I have to run for?  The density just
seems to be oscillating too much and although a literature value of
0.791 g/cm3 is within the data, it doesnt seem to be settling around
that value.  I have similar problems for 1-propanol to where my
simulated density is under that of the literature value.  Could it be
possible that maybe I should use different NVT, NPT algorithms
(tcoupl, pcoupl)?

I would really appreciate anyones help. I'm still in my first few
months of learning the program.  Thanks!

  591.800000  38039.125000  297.409546  -206.071747  779.428894
  592.000000  38035.421875  294.951111  -20.557419  777.816956
  592.200000  37982.812500  298.063324  -463.192047  776.430542
  592.400000  37681.421875  295.430878   36.529854  776.460876
  592.600000  37476.148438  297.193787  -28.961288  775.634644
  592.800000  37611.574219  297.059875  -553.369263  774.096924
  593.000000  37682.277344  288.297974  -293.567963  776.055237
  593.200000  37147.824219  295.985901  -221.989441  781.570679
  593.400000  37404.906250  293.938721  399.093018  788.921570
  593.600000  37339.914062  297.464783  415.891235  793.838440
  593.800000  37909.019531  292.707184  574.818726  796.730286
  594.000000  37369.820312  304.095703  -119.569496  796.235962
  594.200000  37293.445312  294.449005  -63.335049  794.578552
  594.400000  37196.917969  295.323761  104.148239  791.275085
  594.600000  37776.199219  298.004333   28.989655  785.244629
  594.800000  37893.097656  302.397308  -268.065765  779.457886
  595.000000  38123.460938  295.343109  -140.914047  775.682678
  595.200000  38002.054688  302.233124  -91.893478  774.813660
  595.400000  37835.652344  297.808319  441.931885  776.142090
  595.600000  38103.964844  296.369019  490.766754  778.422302
  595.800000  37686.902344  303.358093  -19.726471  781.057922
  596.000000  37586.890625  293.751648  -92.387199  783.724854
  596.200000  37408.414062  300.007385  -288.156036  785.990417
  596.400000  37894.898438  295.720520  346.861206  788.753174
  596.600000  37673.378906  288.311432   89.059952  789.179993
  596.800000  36881.062500  295.265564   70.789131  791.552673
  597.000000  37077.789062  296.261444   72.207787  793.970398
  597.200000  37123.230469  296.511475  -232.434616  793.653564
  597.400000  37426.152344  304.281891  -133.300797  791.791077
  597.600000  37514.332031  305.183197  -259.364136  787.592712
  597.800000  37393.667969  301.326111  -150.162338  785.053711
  598.000000  37263.054688  299.811554  326.851837  785.553894
  598.200000  37203.261719  300.688904  225.526001  788.008667
  598.400000  37511.636719  302.506714   66.299194  790.725281
  598.600000  37460.843750  299.512909  198.943665  793.353333
  598.800000  37248.921875  291.575134  -93.283127  793.582703
  599.000000  36870.480469  293.515381  211.069107  792.929443
  599.200000  37398.097656  293.184448  -163.371140  789.150513
  599.400000  37402.238281  297.402008   45.650658  786.303711
  599.600000  37419.816406  298.520508  -147.688004  783.492371
  599.800000  37497.332031  296.529877  -15.991615  782.294250
  600.000000  37691.468750  298.839844  -68.580627  781.974182



title		=Methanol npt equilibration

;Run parameters
integrator	=md		        ;leap-frog algorithm
nsteps	=300000		;2 * 50000 = 100 ps
dt		=0.002		;2fs

;Output control
nstxout		=100		;save coordinates every 0.2 ps
nstvout		=100		;save velocities every 0.2 ps
nstenergy	        =100		;save energies every 0.2 ps
nstlog		=100		;update log file every 0.2 ps

;Bond parameters
continuation	     =yes		;Restarting after NVT
constraint_algorithm=lincs	        ;holonomic constraints
constraints	             =all-bonds	;all bonds (even heavy atom-H
bonds)constraind
lincs_iter	             =1		;accuracy of LINCS
lincs_order              =4		;also related to accuracy

;Neighborhood searching
ns_type		=grid		;search neighboring grid cells
nstlist		=1		;2 fs
rlist		        =1.4		;short-range neighborlist cutoff (in nm)
rcoulomb	        =1.4		;short-range electrostatic cutoff (in nm)
rvdw		        =1.4		;short-range van der Waals cutoff (in nm)

;Electrostatics	
coulombtype	=PME		;Particle Mesh Ewald for long-range electrostat;ics
pme_order         =4		        ;cubic interpolation
fourierspacing	=0.16		;grid spacing for FFT

;Temperature coupling is on
tcoupl	     	=V-rescale	                ;modified Berendsen thermostat
tc-grps		=SYSTEM	                ;one coupling group
tau_t		        =0.1                          ;time constant, in ps
ref_t		        =298                 	;reference temperature, one for
each ;group, in K

;Pressure coupling is on
pcoupl	  	=Parrinello-Rahman	;Pressure coupling on in NPT
pcoupltype	        =isotropic		        ;uniform scaling of box vectors
tau_p		=2.0			        ;time constant, in ps
ref_p		        =1.0			        ;reference pressure, in bar
compressibility	=4.5e-5		        ;isothermal compressibility of h2O, 1/bar

;Periodic boundary conditions
pbc	  	        =xyz		        ; 3-D PBC

;Dispersion correction
DispCorr	        =EnerPres	        ;account for cut-off vdW scheme

;Velocity generation
gen_vel		=no		;Velocity generation is off



--
Best regards,

Fabian F. Casteblanco
Rutgers University --
Chemical Engineering PhD Student
C: +908 917 0723
E:  fabian.casteblanco at gmail.com



More information about the gromacs.org_gmx-users mailing list