# [lammps-users] which ensemble should I use

Hi, everyone,
I want to relax my system and find the stable positions at the certain temperature.
Both NVE AND NVT and fix the temperature.
But, there are some diffcults about fix a certain temperature by using NVT.
AND, my question is which ensemble is proper for my problem?
Thanks,
–Nwei

A typical use of MD is to equilibrate with NVT, then turn it
off to run NVE and see if the temperature stays constant. If so,
then you're at equilibrium.

Steve

Thinking out loud:

Typically in an MD simulation, one would input atomic positions, assign velocities and make a choice of interatomic potentials. These may be inconsistent with the normal pressure of the system (especially when you are simulating solids). Therefore I would think that a “fix NPT at the required pressure” would be an appropriate choice for finding the stable positions at a certain temperature. Then as Steve says run NVE to see if temperature stays constant.

If one does a “fix NVT” directly on a system of atoms, the system may change volume when one does the “fix NVE” (and this will cause a resultant change in temperature too).

Finally, if one wants to find the stable positions at a given temperature, perhaps a “simulated annealing” method might be more appropriate. This is because that has a better chance than MD to take you to a global minima. In MD, if you DO NOT start with a high enough temperature and slowly cool your system to the desired temperature, you will probably end up at a local minima of potential energy which may not be the most stable position for your system at the desired temperature.

Best’
Manoj

2009/9/14 Steve Plimpton <[email protected]…33…24…>

Thinking out loud:

Typically in an MD simulation, one would input atomic positions,
assign velocities and make a choice of interatomic potentials. These
may be inconsistent with the normal pressure of the system (especially
when you are simulating solids). Therefore I would think that a "fix
NPT at the required pressure" would be an appropriate choice for
finding the stable positions at a certain temperature. Then as Steve
says run NVE to see if temperature stays constant.

manoj,

there is no "one rule that always works(TM)" in MD. what is the best
thing to do, always depends on the circumstances and the kind of system
that you are studying.

If one does a "fix NVT" directly on a system of atoms, the system may
change volume when one does the "fix NVE" (and this will cause a
resultant change in temperature too).

sorry, but how can a system change its volume, when you don't allow it
to change, as is the case in both NVE and NVT?

Finally, if one wants to find the stable positions at a given
temperature, perhaps a "simulated annealing" method might be more
appropriate. This is because that has a better chance than MD to take
you to a global minima. In MD, if you DO NOT start with a high enough

i would like to point out that a "global minimum" can only exists
for 0K. for a system in equilibrium at a temperature > 0K there
will (in a classical system) always be multiple equivalent states.
in fact, the whole point of (many) MD studies is that you sample
those states and average over them, and the problem is always
whether you get a representative sample of them (ergodicity?).

temperature and slowly cool your system to the desired temperature,
you will probably end up at a local minima of potential energy which
may not be the most stable position for your system at the desired
temperature.

if you'd follow this protocol with any weakly bound species, e.g.
biomolecules, you would destroy the rather delicate structure and
get a meaningless simulation. the art of equilibrating biomolecules
lies often in how to get the system moving without losing the structural
information of the starting structure.

the opposite case exists, when you start from a crystal like structure
and want to study the liquid phase. there you need to destroy the
"memory" of the starting structure.

so i would recommend to always be careful with general statements
(like this one ),

cheers,
axel.

Axel,

2009/9/15 Axel Kohlmeyer <[email protected]>

Thinking out loud:

Typically in an MD simulation, one would input atomic positions,
assign velocities and make a choice of interatomic potentials. These
may be inconsistent with the normal pressure of the system (especially
when you are simulating solids). Therefore I would think that a “fix
NPT at the required pressure” would be an appropriate choice for
finding the stable positions at a certain temperature. Then as Steve
says run NVE to see if temperature stays constant.

manoj,

there is no “one rule that always works™” in MD. what is the best
thing to do, always depends on the circumstances and the kind of system
that you are studying.

An example where NPT must be done for equilibration may be helpful in appreciating this (like the example you gave about bio-molecules in the later part of your mail).

If one does a “fix NVT” directly on a system of atoms, the system may
change volume when one does the “fix NVE” (and this will cause a
resultant change in temperature too).

sorry, but how can a system change its volume, when you don’t allow it
to change, as is the case in both NVE and NVT?

Regarding volume change when one does NVT and then NVE, AFAIK there is no active control in MD for keeping volume constant unlike the case with pressure (barostatting) and temperature (thermostatting). Therefore if you do not have ¨periodic boundary conditions¨ in an NVE MD simulation, how will one keep the volume constant?

Finally, if one wants to find the stable positions at a given
temperature, perhaps a “simulated annealing” method might be more
appropriate. This is because that has a better chance than MD to take
you to a global minima. In MD, if you DO NOT start with a high enough

i would like to point out that a “global minimum” can only exists
for 0K. for a system in equilibrium at a temperature > 0K there
will (in a classical system) always be multiple equivalent states.
in fact, the whole point of (many) MD studies is that you sample
those states and average over them, and the problem is always
whether you get a representative sample of them (ergodicity?).

For a group of atoms, say at 300 K, a particular molecular structure may be the most stable state. I would call this structure the global minima which is then computationally best accessed with either a genetic algorithm or a simulated annealing proceedure (http://prola.aps.org/abstract/PRL/v75/i2/p288_1). Of course, in an MD simulation, various degenerate states of this structure can be accessed. From a point of view of accessible states, yes, only one global minimum exists at 0 K. However if one is interested in the minimum energy structure of a system of atoms at some temperature (T), then I would define the minimum energy structure as the global minimum for that problem.

temperature and slowly cool your system to the desired temperature,
you will probably end up at a local minima of potential energy which
may not be the most stable position for your system at the desired
temperature.

if you’d follow this protocol with any weakly bound species, e.g.
biomolecules, you would destroy the rather delicate structure and
get a meaningless simulation. the art of equilibrating biomolecules
lies often in how to get the system moving without losing the structural
information of the starting structure.

the opposite case exists, when you start from a crystal like structure
and want to study the liquid phase. there you need to destroy the
“memory” of the starting structure.

so i would recommend to always be careful with general statements
(like this one ),

Thank you. I was not thinking of bio-molecules when I made the statement.

manoj,

there is no "one rule that always works(TM)" in MD. what is
the best
thing to do, always depends on the circumstances and the kind
of system
that you are studying.
An example where NPT must be done for equilibration may be helpful in
appreciating this (like the example you gave about bio-molecules in
the later part of your mail).

actually, there are many cases where NPT _must_ be done for production
even. some people in our group work on self-assembly and particularly
with smaller size systems, one has to allow for changes in the size and
especially the shape of the unit cell, to form, e.g. a bilayer, as
you have a fixed number of bilayer forming molecules and then there
is an equilibrium surface area per molecule and that has to be found.

running a system at NVE after equilibration is usually meant to be
a test for getting close enough to the equilibrium to start production
and to test whether the simulation is running stable with the given
choice of time step and potentials. in the long run you will always
accumulate kinetic energy from numerical noise and due to using cutoffs
and so on and _that_ has to be taken care of by the thermostat. with
NVE it is difficult to average over the whole trajectory since your E
is not perfectly constant.

> If one does a "fix NVT" directly on a system of atoms, the
system may
> change volume when one does the "fix NVE" (and this will
cause a
> resultant change in temperature too).

sorry, but how can a system change its volume, when you don't
allow it
to change, as is the case in both NVE and NVT?

Regarding volume change when one does NVT and then NVE, AFAIK there is
no active control in MD for keeping volume constant unlike the case
with pressure (barostatting) and temperature (thermostatting).

please have a look at the implementation of MD programs or

if you don't _actively_ change the volume in a simulation
it will stay the same. the NPT (or NPH) integrators take
care of that by coupling the volume change to the computed
pressure (or stress).

Therefore if you do not have ¨periodic boundary conditions¨ in an
NVE MD simulation, how will one keep the volume constant?

if you don't have periodic boundary conditions, you do _not_
have a constant volume. you could say you have an infinite volume.
please don't confuse using the NVT/NVE integrator with having
a simulation in the NVT/NVE ensemble. this reasoning only works
one way, not both.

For a group of atoms, say at 300 K, a particular molecular structure
may be the most stable state. I would call this structure the global

your system is _moving_ all the time, there cannot be _the_ structure,
it always has to be an ensemble. you can average over this ensemble

minima which is then computationally best accessed with either a
genetic algorithm or a simulated annealing proceedure
(http://prola.aps.org/abstract/PRL/v75/i2/p288_1). Of course, in an MD

optimization. it refers to simulated annealing, but also in the
sense of reaching the 0K structure. btw: there are better performing
linear scaling geometry optimization algorithms around now. e.g. L-BFGS.

simulation, various degenerate states of this structure can be
accessed. From a point of view of accessible states, yes, only one
global minimum exists at 0 K. However if one is interested in the
minimum energy structure of a system of atoms at some temperature (T),
then I would define the minimum energy structure as the global minimum
for that problem.

well, you should at the very least call this a minimum of _free_ energy,
as using the plain energy term - at least for me - in the context that
you were using it, would imply talking about the potential energy, which
is at 0K. in any case, if you have a very deep valley in the free energy
surface, you may get something that is close to what a regular
minimization would produce, but you cannot neglect that at finite
temperature your system state is defined by integrating over _all_
accessible configurations. some configurations are more likely than
others, but if you have a shallow potential surface, then you cannot
make statements like you were alluding to.

assume a simple lennard-jones liquid. at 0K you have a minimum energy
structure, at finite temperature, all you can do to describe the system
structure is to look at distribution functions and only because you make
the assumption of looking only at pairwise interactions, you can
describe the system "structure" by a radial distribution function.
yet, there are many, many configurations that contribute to it.

hope this helps to clear up matters.

cheers,
axel.

Axel,

2009/9/16 Axel Kohlmeyer <[email protected]…24…>

manoj,

there is no “one rule that always works™” in MD. what is
the best
thing to do, always depends on the circumstances and the kind
of system
that you are studying.
An example where NPT must be done for equilibration may be helpful in
appreciating this (like the example you gave about bio-molecules in
the later part of your mail).

actually, there are many cases where NPT must be done for production
even. some people in our group work on self-assembly and particularly
with smaller size systems, one has to allow for changes in the size and
especially the shape of the unit cell, to form, e.g. a bilayer, as
you have a fixed number of bilayer forming molecules and then there
is an equilibrium surface area per molecule and that has to be found.

running a system at NVE after equilibration is usually meant to be
a test for getting close enough to the equilibrium to start production
and to test whether the simulation is running stable with the given
choice of time step and potentials. in the long run you will always
accumulate kinetic energy from numerical noise and due to using cutoffs
and so on and that has to be taken care of by the thermostat. with
NVE it is difficult to average over the whole trajectory since your E
is not perfectly constant.

OOps that was a costly mistake. I meant an example where ¨NVT must be done … ¨. Sorry for wasting your time on this. Please note that I was the one who felt that NPT is necessary to relax a system as opposed to NVT.

Peferably it should be an example where periodic boundary conditions (PBC) exist.

If one does a “fix NVT” directly on a system of atoms, the
system may
change volume when one does the “fix NVE” (and this will
cause a
resultant change in temperature too).

sorry, but how can a system change its volume, when you don’t
allow it
to change, as is the case in both NVE and NVT?

Regarding volume change when one does NVT and then NVE, AFAIK there is
no active control in MD for keeping volume constant unlike the case
with pressure (barostatting) and temperature (thermostatting).

please have a look at the implementation of MD programs or

Please note that in the above statement, I meant that there is no active control in MD for keeping Volume constant when one does NVT. Of course when one does NPT one has to change the volume of the system (as you point out below) to maintain the pressure.

This lack of active volume control during NVT is the main reason why I feel that when one starts off a MD simulation with some user-assigned particle positions one must do an NPT relaxation. If one does not do it, the particle positions most probably are not in a relaxed state with respect to the chosen interatomic potential.

if you don’t actively change the volume in a simulation
it will stay the same. the NPT (or NPH) integrators take
care of that by coupling the volume change to the computed
pressure (or stress).

Therefore if you do not have ¨periodic boundary conditions¨ in an
NVE MD simulation, how will one keep the volume constant?

if you don’t have periodic boundary conditions, you do not
have a constant volume. you could say you have an infinite volume.
please don’t confuse using the NVT/NVE integrator with having
a simulation in the NVT/NVE ensemble. this reasoning only works
one way, not both.

I would not say that I have infinite volume if there are no PBCs, at least not when I am simulating a solid / molecule. Based on the potential cutoffs used, one can define a volume.

This is also the reason why I feel LAMMPS should allow NPT relaxation even without PBCs for solid simulations. The atoms on the surface will not evaporate away as long as there is an attractive part to the potential of sufficient depth. An estimate can be made as follows:
if ‘Em’ is the depth of the attractive part of the potential, if ‘Natoms’ is the number of atoms on the surface and if \omega_o is the typical vibration frequency of the surface atoms, then for an MD simulation at a temperature ‘T’, the average number of atoms evaporating from the surface per second is:

Natoms * \omega_o exp^{\frac{-Em}{k_B T}}

Assuming 1000 atoms on the surface, with \omega_o = 1.0e12 (s^{-1}), at a temperature of 300 K, there would be 1.59e-11 desorptions per nano-second simulation from the surface for a Em = 1.0 eV. For the same parameters with an Em of 0.5 eV there would be 0.004 desorptions from the surface in a nanosecond simulation. This would however quickly escalate to 20000 desorptions from the surface in a nanosecond simulation at an Em of 0.1. Therefore depending on Em one can judiciously use non-periodic boundary conditions with NPT.

For a group of atoms, say at 300 K, a particular molecular structure
may be the most stable state. I would call this structure the global

your system is moving all the time, there cannot be the structure,
it always has to be an ensemble. you can average over this ensemble

minima which is then computationally best accessed with either a
genetic algorithm or a simulated annealing proceedure
(http://prola.aps.org/abstract/PRL/v75/i2/p288_1). Of course, in an MD

optimization. it refers to simulated annealing, but also in the
sense of reaching the 0K structure. btw: there are better performing
linear scaling geometry optimization algorithms around now. e.g. L-BFGS.

simulation, various degenerate states of this structure can be
accessed. From a point of view of accessible states, yes, only one
global minimum exists at 0 K. However if one is interested in the
minimum energy structure of a system of atoms at some temperature (T),
then I would define the minimum energy structure as the global minimum
for that problem.

well, you should at the very least call this a minimum of free energy,
as using the plain energy term - at least for me - in the context that
you were using it, would imply talking about the potential energy, which
is at 0K. in any case, if you have a very deep valley in the free energy
surface, you may get something that is close to what a regular
minimization would produce, but you cannot neglect that at finite
temperature your system state is defined by integrating over all
accessible configurations. some configurations are more likely than
others, but if you have a shallow potential surface, then you cannot
make statements like you were alluding to.

assume a simple lennard-jones liquid. at 0K you have a minimum energy
structure, at finite temperature, all you can do to describe the system
structure is to look at distribution functions and only because you make
the assumption of looking only at pairwise interactions, you can
describe the system “structure” by a radial distribution function.
yet, there are many, many configurations that contribute to it.

hope this helps to clear up matters.

cheers,
axel.

Dr. Axel Kohlmeyer [email protected]
Research Associate Professor
Institute for Computational Molecular Science
College of Science and Technology

Thank you!
With Best Regards,
Manoj

This is also the reason why I feel LAMMPS should allow NPT relaxation even without PBCs for solid simulations.

This point has come up several times in the past. I've never understood it.
Physically, if there is no periodic boundary, then there is nothing acting
on the atoms to provide an external pressure. Moreover, it doesn't make
sense to compress the atomic lattice when there is a free surface.
Atoms will just
expand back out it they want to.

Steve

2009/9/16 Steve Plimpton <[email protected]>

This is also the reason why I feel LAMMPS should allow NPT relaxation even without PBCs for solid simulations.

This point has come up several times in the past. I’ve never understood it.
Physically, if there is no periodic boundary, then there is nothing acting
on the atoms to provide an external pressure. Moreover, it doesn’t make
sense to compress the atomic lattice when there is a free surface.
Atoms will just
expand back out it they want to.

NPT relaxation is usually carried out at zero Pressure. For a solid/molecule if there is a free surface (no PBCs), and there exists an attractive part in the interatomic potential, the attractive part prevents the surface atoms from expanding out. Because they have a temperature, the atoms do try to escape the attractive potential (evaporate from the surface). The formula

Natoms * \omega_o exp^{\frac{-Em}{k_B T}} (terms explained earlier in this thread).

yields the number of surface atoms escaping from the attractive part of the interatomic potential per second. If the number of atoms escaping from the surface during the course of the MD simulation is << 1, then one can safely assume that the atoms do not expand out.

Steve

Manoj

dear manoj,

[...]

done .. ¨. Sorry for wasting your time on this. Please note that I
was the one who felt that NPT is necessary to relax a system as
opposed to NVT.

whether you need NVT or NPT depends on your system and your potentials
and what you want to achieve. say you want to simulate a liquid at a
given density, you can compute the required volume up front, so NVT
would work just fine. quite often using the NVE integrator and the
berendsen fixes for temperature or pressure can be a much better
solution than using the NPT integrator. also sometimes using langevin
with a large damping constant can be very helpful to establish
equipartitioning of kinetic energy throughout the system. since
nose-hoover thermostats tend to be somewhat sensitive to feedback
oscillations, it is usually better to apply them (i.e. NPT/NVT)
only after the system is getting close to equilibrium and you want
to start production. in short. it doesn't matter _how_ you get to
equilibrium as long as you get there, but once you start production,
you have to be in a defined ensemble throughout the simulation, or
else you cannot extract meaningful thermodynamic properties from it
(disregarding non-equilibrium or perturbation methods).

[...]

Please note that in the above statement, I meant that there is no
active control in MD for keeping Volume constant when one does NVT. Of
course when one does NPT one has to change the volume of the system
(as you point out below) to maintain the pressure.

why do you keep saying this? the volume is controlled,
it is well known, as it is _constant_, it is an input parameter.

This lack of active volume control during NVT is the main reason why I
feel that when one starts off a MD simulation with some user-assigned
particle positions one must do an NPT relaxation. If one does not do
it, the particle positions most probably are not in a relaxed state
with respect to the chosen interatomic potential.

as written above, this entirely depends on the system
and the state you want to study. for some it is required,
for others not.

[...]

I would not say that I have infinite volume if there are no PBCs, at
least not when I am simulating a solid / molecule. Based on the
potential cutoffs used, one can define a volume.

a) using cutoff to define the volume of an individual
molecule or crystal is plain nonsense as this has usually
no relation to the (relative) particle size in a system.
other - effectively arbitrary - parameter, _but_ ...

if you don't use PBC or don't have walls, or a mixture
of that your _system_ has an infinite volume. this is
the volume that matters to the _ensemble_.

This is also the reason why I feel LAMMPS should allow NPT relaxation
even without PBCs for solid simulations. The atoms on the surface will

this is also complete nonsense. if you have no boundaries, the system
will relax any which way it likes, and if you have no other particles
in the system (gas or liquid) that would interact with this solid and
thus affect the equilibrium structure of said solid, then it will
always behave the same way. but in the case of a two phase system, you
need - again - walls or PBC to have a finite volume that a pressure in
the system can build at all.

not evaporate away as long as there is an attractive part to the
potential of sufficient depth. An estimate can be made as follows:
if 'Em' is the depth of the attractive part of the potential, if
'Natoms' is the number of atoms on the surface and if \omega_o is the
typical vibration frequency of the surface atoms, then for an MD
simulation at a temperature 'T', the average number of atoms
evaporating from the surface per second is:

Natoms * \omega_o exp^{\frac{-Em}{k_B T}}

Assuming 1000 atoms on the surface, with \omega_o = 1.0e12 (s^{-1}),
at a temperature of 300 K, there would be 1.59e-11 desorptions per
nano-second simulation from the surface for a Em = 1.0 eV. For the
same parameters with an Em of 0.5 eV there would be 0.004 desorptions
from the surface in a nanosecond simulation. This would however
quickly escalate to 20000 desorptions from the surface in a nanosecond
simulation at an Em of 0.1. Therefore depending on Em one can
judiciously use non-periodic boundary conditions with NPT.

no one cannot. the particles that desorb have to be able to come back
(via PBC or a wall) or else you have an infinite volume again and
no need for NPT. i suggest you spend some quality time with your
favorite stat mech book and look up how ensembles are defined and
what they mean in the practical application of MD simulations.
once you understand how ensembles relate to their implementation
in MD codes, it should become very obvious that NPT for a
non-boundary system is meaningless.

cheers,
axel.

Dear Axel,

Firstly, thank you for the example where NVT works just fine.

Please bear with me as I briefly state what I wish to say.
Carrying out an MD simulation usually consists of two main parts:

1. You set up the system … i.e. assign particle positions, velocities, etc and relax the system.
2. Choose an appropriate ensemble, and sample the required physical quantities of interest.

An NPT is usually required in part (1) above. If you are simulating any system with a surface (say surface bombardment at room temperature and pressure), you would want to do an NPT of the system so that the positions chosen by you are consistent with the interatomic potential and would expand / contract such that the pressure of the system oscillates around the ambient pressure. Because you have a surface at one plane using PBC is ruled out.

Part (2) depends on what one wants to measure and is not relevant to our discussion.

LAMMPS does not allow one to carry out Part (1) in a straightforward way. Other codes (HCParCas, Maintainer: K. Nordlund, Univ. Of Helsinki) allow NPT with non-periodic BCs and I do not think that is wrong -> Note that in Part (1) we are only relaxing the system and not carrying out production runs of an appropriate ensemble to glean physical quantities of interest. The atoms on the surface do not evaporate away during the NPT simulation, for reasons based on statistical mechanics, presented as a formula in my earlier mail.

That said, LAMMPS is my favourite tool and I am highly motivated to convert my codes into C++ with a similar flexible structure and great documentation.

With Best Regards,
Manoj

2009/9/16 Axel Kohlmeyer <[email protected]>

Dear Axel,

dear manoj,

Please bear with me as I briefly state what I wish to say.
Carrying out an MD simulation usually consists of two main parts:

an _equilibrium_ MD simulation, right?

1) You set up the system .. i.e. assign particle positions,
velocities, etc and relax the system.

ok.

2) Choose an appropriate ensemble, and sample the required physical
quantities of interest.

ok.

An NPT is usually required in part (1) above. If you are simulating

i object to the term _required_. ...and also we should be specific
here again. are you talking about "using fix NPT in LAMMPS" or
"running in NPT ensemble"? strictly speaking, you are not running
in NPT ensemble when you are still equilibrating, ... and - as i was
pointing out before - when relaxing a system in LAMMPS 'fix npt' can
be quite often less efficient to other methods, e.g. using fix nve
in combination with fix temp/berendsen or langevin and fix
press/berendsen as those fixes are much better suited to remove or
add larger quantities for energy or volume and are less sensitive
to feedback fluctuations.

any system with a surface (say surface bombardment at room
temperature and pressure), you would want to do an NPT of the system
so that the positions chosen by you are consistent with the
interatomic potential and would expand / contract such that the
pressure of the system oscillates around the ambient pressure. Because
you have a surface at one plane using PBC is ruled out.

this is wrong again. if you have a surface, you still (can) have
PBC in two dimensions and for _those_ two dimensions you can
run 'fix npt'. the free surface is free and if there is nothing
colliding with it, it is not experiencing any pressure.

LAMMPS does not allow one to carry out Part (1) in a straightforward
way. Other codes (HCParCas, Maintainer: K. Nordlund, Univ. Of
Helsinki) allow NPT with non-periodic BCs and I do not think that is
wrong -> Note that in Part (1) we are only relaxing the system and not

well, at the very least it doesn't make sense to call it NPT.
you can call it relaxation with a simulated pressure, but even
then i would be very doubtful of its foundation in the statistical
mechanics framework. to me it simply doesn't make sense to consider
this NPT. also, if you bombard a surface with particles, you
are _not_ running a non-equilibrium MD unless you have a way to do
this as a grand canonical ensemble (which i would be curious to
see how this could be set up), so your ensemble is something
else altogether.

carrying out production runs of an appropriate ensemble to glean
physical quantities of interest. The atoms on the surface do not
evaporate away during the NPT simulation, for reasons based on
statistical mechanics, presented as a formula in my earlier mail.

this is irrelevant. you system is the _total_ system. if particles
stay or evaporate from the surface, doesn't make the tiniest difference.
if you remove or add particles to a system you do not have a fixed
'N' to begin with. again, the issue here seems to be that you
consistently confuse the "fix" commands in LAMMPS or other codes
with the statistical mechanical ensembles, and have a tendency to
make oversimplified and generalized statements, that don't hold water.
please keep in mind that there are many beginners subscribed on this
list and they are easily confused by sloppy wording, so it is very
important to try and be precise and careful with general statements.
the latter is the reason why i am so insisting and hardheaded on
this subject.

thanks for you understanding,
axel.

I commend Axel and Manoj on keeping this lively discussion civil, for the most part. I agree with Axel, in that Manoj’s request for an NPT fix to relax the pressure in non-periodic systems is not strictly necessary. Atoms in a non-periodic system are free to move father apart or closer together, until an equilibrium pressure of zero is approached. Adding in some kind of effective barostat in the non-periodic direction might make it a little easier to get to that zero-pressure equilibrium state, but there are lots of other solutions already available in LAMMPS.

That said, at the risk of pushing Axel over the edge, I will mention that there are circumstances imposing a non-zero pressure in a non-periodic system would be useful. Experimentally, a finite cluster of condensed material can be pressurized using an inert gas. The gas can be thought of as a stochastic momentum bath that exerts a positive pressure on the free surfaces of the cluster. Here is a link to an MD paper where this situation was simulated explicitly.

http://rafiki.tau.ac.il/~rabani/papers/paper61.pdf

They needed up to 350,000 gas particles to pressurize nanocrystals composed of only 300-4500 atoms. I imagine that non-periodic shrink-wrap boundary conditions, combined with “fix npt aniso”, would give the same result alot more easily.

Aidan