need practical advice

LAMMPS LAMMPS Mailing List Mirror
1

Hello,

I am a new to classical MD and LAMMPS. I want to calculate some properties of my system in NVT ensemble at 300K. The strategy I am using is described as follows (in serial order):

(i) minimize: the structure to really small forces.
(ii) NVE: to equilibrate the minimized structure.
(iii) NVT: to 300K using Nose-Hoover thermostat.
(iv) NPT: since the volume was always fixed in the above calculations, I thought, I will perform NPT at 300K and 0 bar to get the correct volume, and atomic coordinates.
(v) NVE: to equilibrate the output system from NPT.
(vi) Finally, NVT at 300K to perform measurements I need.

Is the above wrong/too much ? What should be the correct strategy ?

In the above I face some problems. The output pressure is actually 0 at the end of the step (ii) and (iii). But the pressure fluctuates to several tens of thousand bars in NPT stage (step iv). As a result I get wrong atomic coordinates and structure. I read the forum and found that Nose Hoover barostat can fluctuate a lot. So, I have been playing with drag parameter in Nose Hoover. Upto certain steps, quantities are fine, but a sharp change in volume and pressure happens at some step in step (iv), ie. NPT. Then pressure starts fluctuating quite a bit and never stabilizes and the structure also goes bad.

Finally, I decided to change from NPT to:

fix 1 all nvt temp 300.0 300.0 100.0
fix 2 all press/berendsen iso 0.0 0.0 100000.0

I figured, the simulation is very sensitive to Pdamp above. When Pdamp = 1000, several quantities such as volume, pressure, etc. become “nan” very quickly, so, changed Pdamp = 10000. It ran well for some time and then the program stopped. Currently, I am running with Pdamp = 100000. It is still running.

My question is:
(i) If my strategy described above in steps (i) to (vi) is ok ? Do I need to relax volume after NVE and NVT for calculation of actual quantities or since the pressure is zero after NVE and NVT, it is alright to go ahead with other calculations without any need to perform NPT ?
(ii) How to choose the parameters of NPT (or press/berendsen) correctly for a given a system ?

Thanks,

George

2

Hi George,

The steps look fine. Huge pressure fluctuation is actually normal. The smaller the system, the larger the fluctuation, and fluid system has a higher fluctuation than solid. What properties are you intending to calculate? Quantities extraction by averaging should be performed in NPT, I’m assuming that you are comparing with experimental data. If after a reasonable equilibration, and your system still goes wild, probably your potential that you are using does not captured the system well. Also, experiments are preformed in 1 atm or 1 bar, not 0 bar. Hope this helps.

wf

3

Hi George,

The steps look fine. Huge pressure fluctuation is actually normal. The
smaller the system, the larger the fluctuation, and fluid system has a
higher fluctuation than solid. What properties are you intending to
calculate? Quantities extraction by averaging should be performed in NPT,
I'm assuming that you are comparing with experimental data. If after a
reasonable equilibration, and your system still goes wild, probably your
potential that you are using does not captured the system well. Also,
experiments are preformed in 1 atm or 1 bar, not 0 bar. Hope this helps.

considering that condensed systems are not very compressible.
the difference between 0 bar and 1 kbar is negligible. with the huge
fluctuations, you'll almost never get a good enough sampling to have
this one converged that well. :wink:

axel.

4

Thanks for your posts. How is Pdamp implemented within LAMMPS ? Is there a paper where it is described ? Just wonder if there is a way to systematically evaluate Pdamp for a given system for a good initial guess?

Thanks,

George

5

I too was a little bit confused by the difference between "damping"
and "dragging". I don't know if it helps, but from a casual look at
the fix_hp.cpp file, my guess is is looks like the Pdamp and Tdamp
parameters are used to calculate to the Nose-Hoover artificial
mass(es). (It's more convenient for users to just specify the
timescale of Nose-Hoover oscillations than to specify all the
Nose-Hoover masses.)

The "drag" parameter looks quite different. Using the drag parameter
will speed up equilibration but will generate data whose fluctuations
are not consistent with the NVT, NPT, or NPH ensemble.

If you want more specific details, let us know.
If I'm wrong, please someone please correct me!

Andrew

"Note that use of the drag keyword will interfere with energy
conservation and will also change the distribution of positions and
velocities so that they do not correspond to the nominal NVT, NPT, or
NPH ensembles."

6

Hi George,

If you are referring Pdamp to the pressure damping param in NPT, the document says it well.

The Tdamp parameter is specified in time units and determines how rapidly the temperature is relaxed. For example, a value of 10.0 means to relax the temperature in a timespan of (roughly) 10 time units.

and similar to Pdamp…

IMPORTANT NOTE: A Nose-Hoover barostat will not work well for arbitrary values of Pdamp. If Pdamp is too small, the pressure and volume can fluctuate wildly; if it is too large, the pressure will take a very long time to equilibrate. A good choice for many models is a Pdamp of around 1000 timesteps. Note that this is NOT the same as 1000 time units for most units settings.

This is the reason why your system collapsed when you used a small pdamp. Choosing pdamp can be rigorous. In my case where timestep = 1fs, I used a pdamp value of 1000.
Hope this helps.

wf

7

Hello !

Thanks for your emails. I had some opportunity to investigate the matter. In one literature I found the effective relaxation time in Nose Hoover (NVT) for temperature as

\tau^2 = Q/(N_df * k_b * T_0)

where Q = coupling strength and N_df = number of degrees of freedom, T_0 = temperature of heat bath (the desired final temperature).

For my system with 768 atoms at 300K, N_dfk_b * T_0 ~ 19.2 eV and I have chosen a timestep of 0.1 fs. My system works well for Tdamp = 100, fs^2. So, Q = 19.2(100)^2 = 192000 eV.

Wonder, if the above seem ok. I do not have such a formula for barometer.

A few remarks:

(i) I have 760 atoms in a box (orthogonal) with no vacuum (with periodic boundary condition)
(ii) I tried large damping (upto 100,000) for Berendsen and NPT. For one case, I had nice result for Berendsen but never for NPT.
(ii) The dilemma I have is that I guess, I need to relax volume when temperature is not zero ? Am I correct ? Its clear that from my NPT and Berendsen calculations that volume does change ?
(iv) If NPT does not work very well, then how do experts in the field take care of volume ? Is there another method to do that ?
(v) Or, I am ok if I don’t relax volume and simply perform calculations within NVT as that is eventually my aim.

Thanks,

George

8

Hi George,

I’m not sure about the formula of barastat, maybe Dr. Kohlmeyer or other experts may want to comment.
However, I have a few suggestions which may be useful.
You said you used the pdamp up to 100,000 and occasionally got nice results. I think this sounds unreasonable because you must always end up in good result if the parameters you used are correct. Perhaps you did not performed proper sampling?
Another is whether you should relax your system in NPT before performing your ultimate NVT. My suggestion is if you are certain that your system is in correct arrangement, then I think it is ok to straight away proceed to NVT. In the case of measuring thermal conductivity, in which the calculation must be performed in NVT, I did not relax my system in NPT prior to NVT, I arranged the atoms according to the experimental structure and straight away use NVT.
Hope this helps.

Christopher

9

Hello !

Thanks for your emails. I had some opportunity to investigate the matter. In
one literature I found the effective relaxation time in Nose Hoover (NVT)
for temperature as

\tau^2 = Q/(N_df * k_b * T_0)

where Q = coupling strength and N_df = number of degrees of freedom, T_0 =
temperature of heat bath (the desired final temperature).

For my system with 768 atoms at 300K, N_df*k_b * T_0 ~ 19.2 eV and I have
chosen a timestep of 0.1 fs. My system works well for Tdamp = 100, fs^2.
So, Q = 19.2*(100)^2 = 192000 eV.

Wonder, if the above seem ok. I do not have such a formula for barometer.

A few remarks:

some comments

(i) I have 760 atoms in a box (orthogonal) with no vacuum (with periodic
boundary condition)

this is a very small system and hence you will have very large
pressure fluctuations. those will "upset" the nose-hoover algorithm.
you can limit this effect using the "drag" flag, but that is not a
real solution for a production calculation, since you won't be
sampling a well defined statistical mechanical ensemble anymore.
a better approach would be to enlarge the system using the
replicate command, e.g. 'replicate 2 2 2' or 'replicate 3 3 3'
which will multiply the number of atoms by 8 or 27, respectively.

(ii) I tried large damping (upto 100,000) for Berendsen and NPT. For one
case, I had nice result for Berendsen but never for NPT.

that is by construction. the berendsen algorithm is simple
but robust. however with such a large damping factor, it
will take forever to converge. better to enlarge your system.

(ii) The dilemma I have is that I guess, I need to relax volume when
temperature is not zero ? Am I correct ? Its clear that from my NPT and
Berendsen calculations that volume does change ?

whether you need to "optimize" the volume or not, depends on
what you want to get out of the simulation. if you know exactly
what density you want to simulate, you can just rescale the
box to a suitable value and then relax.

(iv) If NPT does not work very well, then how do experts in the field take
care of volume ? Is there another method to do that ?

since classical MD is "cheap" nobody runs such tiny systems
except for debugging or other special purposes. with a system
or a few 10,000 atoms or 100,000 atoms you'll see that problems
resulting from pressure fluctuations are much less. also, you
can use a more robust scheme like berendsen for as long as you
are far from equilibrium and then switch to npt with nose-hoover,
or first use the damp flag and then turn it off. if CPU time is a
problem, you can start equilibration with a smaller system and
then replicate to full size when you are close to equilibrium.

(v) Or, I am ok if I don't relax volume and simply perform calculations
within NVT as that is eventually my aim.

impossible to say in general terms. you should not do things
because somebody says so, but because there is a physical
or technical reason to do so. what you typically need in MD
is a configuration in equilibrium. how you get there, is not
very relevant. different approaches are differently efficient.

axel.

10

From the fix nvt doc page:

IMPORTANT NOTE: A Nose-Hoover thermostat will not work well for
arbitrary values of Tdamp. If Tdamp is too small, the temperature can
fluctuate wildly; if it is too large, the temperature will take a very
long time to equilibrate. A good choice for many models is a Tdamp of
around 100 timesteps. Note that this is NOT the same as 100 time units
for most units settings.

You don't say what units you are using. If it's metal then time units
are psecs. And a Tdamp of 100 psecs with a timestep of 0.1 fsec
would mean relaxation on the order of 100,000 timesteps, which is
far too slow. You should be able to use a Tdamp = 0.1 for metal
units.

Steve