[gmx-users] Unsteady Density

[Justin A. Lemkul]

Fabian Casteblanco wrote:
> 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.
> 

nstlist = 1 is probably overkill.  Such frequent updating should only be 
necessary for EM.

> 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

Contraints should be applied unless you want your timestep to be limited to 
0.5-1 fs.  Being that you're using a 2-fs timestep, yes, constraints are 
necessary.  There is a discussion about all of this in the manual and in the 
LINCS papers.

> 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

How far off from the literature value is your result?

> posted the last several lines from the analysis *.xvg file.  NPT
> script is also below.

The posted data are meaningless without labels to know what we're looking at. 
Raw output from .xvg files is usually not useful.  Averages and fluctuations 
printed from .edr files, however, is more useful, since it is accurate over all 
MD steps, whereas the resulting .xvg files contain only frames corresponding to 
nstenergy intervals.

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

Your settings should correspond to an accurate ensemble, but I suppose you could 
mess around with them.  Berendsen should be avoided.

The one factor you've not really mentioned is your topology.  One of the most 
important factors in determining the accuracy of your simulation is the 
parameters assigned to the actual molecules.  How did you derive the topologies? 
  Are the parameters used (charges, atom types, etc) consistent with any 
literature precedent?

> compressibility	=4.5e-5		        ;isothermal compressibility of h2O, 1/bar

Note that you're using the compressibility for water rather than your actual 
systems.  This could explain some discrepancy.

-Justin

-- 
========================================

Justin A. Lemkul
Ph.D. Candidate
ICTAS Doctoral Scholar
MILES-IGERT Trainee
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin

========================================