# [lammps-users] Fluctuations in NPT ensemble

Dear LAMMPS users,

I am trying to estimate the fluctuations in the enthalpy for the bulk sylicon. I simulated it by using SW potential. The number which I got is basically enormous mainly due to the enormous fluctuations in the pressure. Could you propose any measures to reduce the fluctuations in the npt ensemble.

Dear LAMMPS users,
I am trying to estimate the fluctuations in the enthalpy for the bulk
sylicon. I simulated it by using SW potential. The number which I got is
basically enormous mainly due to the enormous fluctuations in the pressure.
Could you propose any measures to reduce the fluctuations in the npt
ensemble.

the pressure fluctuations are a result of the low compressibility
and finite system size. for any proper analysis, you should use
the average over a large enough chunk of trajectory.

please note that you cannot immediately apply "continuum thinking"
to atomistic simulations. if you have not noticed by now, there is a
whole discipline called "statistical mechanics" devoted to making
this kind of connection.

cheers,
axel.

Well it is rather clear to me that I have to use averages over
sufficient long trajectories but it is not clear to me how the
standard deviation of pressure fluctuations can be equal to thousands
of atm. The second thing which is not clear is why those deviations
are not really decreasing with the increase of the size of the system
( as happens with Temperature fluctuations in NVT ensemble ). I am
aware of those fluctuations but do not understand how they can be so
large and how to make them smaller, because at the present stage I
cannot really calculate what I need by using those fluctuations. I
tried to calculate specific heat capacity through entalpy fluctuations
at the constant pressure but got very large results mainly due to the
large pressure fluctuations.
I am using SW potential;

Initially I relaxed my system by:

fix npt all 300 300 100 xyz 0 0 1000;

then initialized velocities and run one more time

fix npt all 300 300 100 xyz 0 0 1000;

I used reall units the mass of atoms was 28.0855 ( does it matter
for pressure fluctuations ) the number of atoms was 8000.

I also tried different Pdamp but it did not help.

Well it is rather clear to me that I have to use averages over
sufficient long trajectories but it is not clear to me how the
standard deviation of pressure fluctuations can be equal to thousands
of atm. The second thing which is not clear is why those deviations
are not really decreasing with the increase of the size of the system
( as happens with Temperature fluctuations in NVT ensemble ). I am

they should indeed go down. but you also have to consider
the fact that you have a material with very little compressibility.

aware of those fluctuations but do not understand how they can be so
large and how to make them smaller, because at the present stage I
cannot really calculate what I need by using those fluctuations. I
tried to calculate specific heat capacity through entalpy fluctuations
at the constant pressure but got very large results mainly due to the
large pressure fluctuations.

the second thing worth checking out would be whether you have
some intrinsic lattice vibrations that are generated by your starting
conditions and don't die off quickly.

what you could try doing is a three step procedure:
- start with a reasonable lattice constant and a reasonably large system
and run fix nve with fix langevin. i would play with the gamma value a little
bit to that temperature fluctuations are dampened but not too much.
check what is the pressure fluctuation for this system.
- then add press/berendsen to relax the system cell and wait until well
equilibrated and
- then switch to fix npt (i.e. unfix press/berendsen, langevin and nve)
and see if you still get those fluctuations.

if you have lattice phonons generated by your initial conditions, it can
take a very long time until they die off, but with fix langevin you should
be able to effectively dissipate them.

cheers,
axel.

I tried to execute following lines but got errors

fix 2 all press/berendsen xyz 3633.7839 3633.7839 1000.0
fix 3 all nve

I have never used press/berendsen could you tell what is wrong ?

is it implemented in march version of LAMMPS ?

I tried to execute following lines but got errors

_what_ errors. it is impossible to debug something without
knowing any details about the symptoms. you don't go to the
doctor and say "it hurts somewhere", or do you?

fix 2 all press/berendsen xyz 3633.7839 3633.7839 1000.0
fix 3 all nve

try reversing the order. also you should consider adjusting the

I have never used press/berendsen could you tell what is wrong ?

without any information or even better a way to reproduce the problem,
i can only do a wild guess.

axel.

After executing:

velocity all create 600 41211 mom yes rot yes dist gaussian
fix 1 all nve
fix 2 all press/berendsen xyz 3633.7839 3633.7839 1000.0

Lattice spacing in x,y,z = 5.4307 5.4307 5.4307
Created orthogonal box = (0 0 0) to (54.307 54.307 54.307)
2 by 2 by 2 processor grid
Created 8000 atoms
Setting up run ...
Memory usage per processor = 1.97878 Mbytes
Step Atoms Temp TotEng PotEng KinEng Volume Press Enthalpy
0 8000 600 -785693.88 -799999.98 14306.1
160164.93 4221.2018 -775833.82
100 8000 nan nan nan nan
nan nan nan

After executing:

velocity all create 600 41211 mom yes rot yes dist gaussian
fix 1 all nve
fix 2 all press/berendsen xyz 3633.7839 3633.7839 1000.0

Lattice spacing in x,y,z = 5.4307 5.4307 5.4307
Created orthogonal box = (0 0 0) to (54.307 54.307 54.307)
2 by 2 by 2 processor grid
Created 8000 atoms
Setting up run ...
Memory usage per processor = 1.97878 Mbytes
Step Atoms Temp TotEng PotEng KinEng Volume Press Enthalpy
0 8000 600 -785693.88 -799999.98 14306.1
160164.93 4221.2018 -775833.82
100 8000 nan nan nan nan
nan nan nan

you didn't follow my suggestion to first equilibrate
with fix langevin. i did suggest this for a reason and
even explained it. you are free to do whatever you
want, but don't ask for help if you don't intend to follow
what was recommended to do.

axel.

With this I tried to understand why press/berendsen did not work for
Axel really thanks for your input, I am just a newbie in the field and
there is almost no literature about how to do the simulation. There
are a lot of sophisticated math texts about simulation but for a
newbie it is very difficult to read those texts.

I understand that this is the relaxation time but again it is not
clear how to choose it ( as a function of time step ) what is even
more striking I tried to track the functional dependence of
fluctuations from this parameter and unfortunately did not get
monotonic function.

Again I am asking this question but could you please suggest some
texts where it is clearly explained how this langevin thermostat,
berendssen barostat, Nose - hoover type barostats and thermostats are
implemented in velocity verlet algorithm. I understand your suggestion
to read Physical Review or JCP but sometimes such journals are not

I am asking those questions because I really want to understand what
is happening and not only use LAMMPS as a black box. Sometimes there
are those drag parameters ( NPT ), which one can use to achieve the
result but then I really do not understand whether it is really
physical.

I would really appreciate if you could provide sources of
implementation of verlet with different barostats and thermostats and
on the intuition behind those barostats and thermostats.

It seems that your approach worked quite well for me. Could you please
explain a bit more details behind such methodology or provide some
sources so I could you understand this phenomena much better.

If you have problems with NPT or berendsen, I would run straight NVE
first, or NVT (or langevin) next. If your system is stable in those modes,
with reasonable pressures, then turning on NPT will likely work. If not,
then NPT will likely have problems.

Steve

There
are a lot of sophisticated math texts about simulation but for a
newbie it is very difficult to read those texts.

to read Physical Review or JCP but sometimes such journals are not

I am asking those questions because I really want to understand what
is happening and not only use LAMMPS as a black box.

Mezakuilis,

v/r,

Dan

> > I am asking those questions because I really want to understand
what
> > is happening and not only use LAMMPS as a black box.

Mezakuilis,

well, you probably know the saying: no pain, no gain.

but enough of the cheap shots. there are a number of

fist off, you have to distinguish between learning the
method and learning how this method is implemented in
a particular piece of software.

a mailing list is _not_ the place to learn to learn the
first thing. there are quite a few things to learn and
those are covered by lectures and text books on statistical
mechanics and thermodynamics and molecular simulations.

if you have not taken these classes or are willing to
learn these things in self-study, you _will_ have a lot
of pain. there are many subtle issues that are trivial
once you understand the principles, but confusing if you
the method, you are practically putting the cart before
the horse and the result is rarely satisfactory.

to me personally, it is quite troublesome to see this
happening more and more over the last years as it
basically means that either advisors don't do their
job properly, or that students don't know how to learn
by themselves anymore. which of the two applies to
your cases, i don't know and it doesn't really matter.

with that all in mind, i can only encourage you to not
but rather spend some time to study the principles of
stat mech and molecular simulations from text books.

one of the best exercises that i went through (and for which
you can find a number of tutorials on the web) is to write
a toy MD code for yourself (for a simple one atom type
lennard-jones potential) and implement whatever you want
to understand on top of that. you have to keep in mind
that MD is pretty much like experiments and for those some
fundamental practical skills and practicing on simple
cases is a good preparation.

one thing that you will have to accept is that there is
not going to be _the_ textbook or article that will contain
exactly and only what you want to learn in the way that
you can learn it best. you will have to pick up some pieces
here and there and also have to learn how to look for
what you need to learn. the latter is as much a useful
skill as understanding the simulation techniques as well.

i hope that will help you to deal better with the pain
and frustration that you are feeling.

cheers,
axel.

Hoover thermostats and Barostats in velocity verlet. Is iterative
velocity - verlet used or Liuville operator ?

if you want to know this for certain, there is no
alternative to reading the source code. if you
cannot deal with that, then you are out of luck.

as i wrote before, you most likely won't find a
customized translation of the theory or implementation
so that it is easy for you to understand.

using langevin and equilibration at fixed volume helps,
i can only paraphrase the explanation, that i motivated
when making the original suggestion. your initial conditions
are very likely not commensurate with the potential that
you are using. as a consequence you would get strong
lattice vibrations. those don't die off quickly in periodic
boundary conditions. due to its dissipative nature, the
langevin thermostat helps to dampen and randomize
these vibrations and as a consequence the pressure
will fluctuate less. once your system is equilibrated
at NVT, the pressure and thus the volume changes
in the variable cell procedure will be less. due to the
very low compressibility of your material, it is very easy
to get large pressure fluctuation and resulting volume
changes. if you don't pre-equilibrate at NVT, they
are even larger. you also seem to add a large pressure
to the system to begin with. that just makes matter worse.

in general, it is always a good idea to smoothly guide
a system to equilibrium and not just set something up,
let it go and assume that it will go into the right direction
by itself. it will do, but it would take a long, looonng time.
it is better to be more gentle.

cheers,
axel.

LAMMPS is velocity Verlet. The NH work is in src/fix_nh.cpp.
You're welcome to look at the source code.

Steve