Advice for keeping pressure fluctuations low

LAMMPS LAMMPS Mailing List Mirror
1

I am running a simulation of silicon atoms, first in NPT ensemble and the pressure fluctuations are very low (plus or minus 20 bars) in my 14mil atom simulation, which is great. But when I switch from NPT ensemble to NVE, I get fairly large pressure fluctuations (for instance right now the pressure is at ~800 bars). I believe it will swing back down the other way as well. Are there any measures that can be taken to reduce this? I understand that with the NVE ensemble, I am fixing the volume and energy and thus pressure fluctuates, but I would not except such a huge pressure wave in the system under equilibration. Is this sort of thing normal to see, or is there a procedure for reducing these fluctuations?

Jackson

2

I am running a simulation of silicon atoms, first in NPT ensemble and the
pressure fluctuations are very low (plus or minus 20 bars) in my 14mil atom
simulation, which is great. But when I switch from NPT ensemble to NVE, I
get fairly large pressure fluctuations (for instance right now the pressure
is at ~800 bars). I believe it will swing back down the other way as well.

that doesn't sound like a fluctuation. perhaps, your volume for NVE
was by chance significantly away from the equilibrium volume.

Are there any measures that can be taken to reduce this? I understand that

please keep in mind, that different materials have a different
compressibility leading to either larger or lesser changes in pressure
due to a given change in volume.
when you switch from NPT to NVE, you don't want to switch to the
*instantaneous* box dimensions, but rather sufficiently averaged ones.
there are other factors, that have an impact as well.

with the NVE ensemble, I am fixing the volume and energy and thus pressure
fluctuates, but I would not except such a huge pressure wave in the system
under equilibration. Is this sort of thing normal to see, or is there a
procedure for reducing these fluctuations?

in a properly tuned and equilibrated NPT simulation, the switch to NVE
should not incur significantly larger pressure fluctuations. the time
constant for adjusting the volume with fix npt should be
*significantly* larger than typical pressure fluctuations, or else the
system will couple to strongly to the fictitious DOFs of the
thermostats, and unphysical things may happen. if you use a "drag"
factor to artificially suppress too strong reactions to pressure
fluctuations, then the switch to NVE will simply demonstrate, how much
you have impacted the simulation before.

it is difficult to make more specific statements without knowing more
details. please keep in mind, that you get out what you put in, since
computer software doesn't know anything about "reasonable". it simply
blindly does exactly what you tell it, even if that is not what you
mean.

axel.

3

Axel,

What would be the best way to switch to NVE without using the instantaneous box dimensions? Should I write a new data file and just manually change the x y and z dimensions of the NPT to be that which is averaged?

Or is there a way within the script itself to change the box dimensions based on the average volume. I did notice that the volume of my simulation was slowly decreasing throughout NPT equilibration and probably would have wanted to continue to decrease, but I started the NVE ensemble.

I think the first course of action that I should do would be to equilibrate longer.

I wish I could say I understand this sentence:

" the time
constant for adjusting the volume with fix npt should be
significantly larger than typical pressure fluctuations, or else the
system will couple to strongly to the fictitious DOFs of the
thermostats, and unphysical things may happen. "

But I cannot.

In my fix npt command I specify the pressure and temperature and use damping factors 100x and 1000x the timestep for temperature and pressure, respectively, as the LAMMPS manual suggests. Which of these is the time constant for adjusting the volume?

Thanks for your assistance.

Jackson

4

Axel,

What would be the best way to switch to NVE without using the instantaneous
box dimensions? Should I write a new data file and just manually change the
x y and z dimensions of the NPT to be that which is averaged?

that is one way to do it and not the worst. as all automatic
procedures depend on heuristics and guesses and thus may not make the
best choice unless you are running a system, that is well understood
and well behaved. you can compute (time window) averaged box
dimensions with fix ave/time and then use the result in a change_box
command.

Or is there a way within the script itself to change the box dimensions
based on the average volume. I did notice that the volume of my simulation
was slowly decreasing throughout NPT equilibration and probably would have
wanted to continue to decrease, but I started the NVE ensemble.

I think the first course of action that I should do would be to equilibrate
longer.

yes. outside of the volume, you should also monitor the potential
energy to identify, when you have reached equilibrium.

I wish I could say I understand this sentence:

" the time
constant for adjusting the volume with fix npt should be
*significantly* larger than typical pressure fluctuations, or else the
system will couple to strongly to the fictitious DOFs of the
thermostats, and unphysical things may happen. "

But I cannot.

In my fix npt command I specify the pressure and temperature and use damping
factors 100x and 1000x the timestep for temperature and pressure,
respectively, as the LAMMPS manual suggests. Which of these is the time
constant for adjusting the volume?

please read the fix npt docs again. carefully. it should be obvious.
if not, you should probably look into some additional papers or
textbooks that explain this.

you should also pay attention to the change of volume during
equilibration. if the volume first increases and later decreases, you
probably should not start with fix npt, but with fix nvt.
there are many more subtle issues to pay attention to, especially if
your system is not very compressible and has few or no defects. there
should be multiple useful discussions on this topic in the mailing
archives.

axel.

5

OK sounds good!

I will check through the archives and take all of your suggestions into account.
I appreciate all of your help!

Jackson

6

Why do you want small pressure fluctuations? Fluctuations are not intrinsically bad. Also, fluctuations in NPT should be of similar magnitude to those in NVE.

7

I understand that pressure fluctuations are generally normal. I just thought this was more of a concern, since the pressure fluctuations were much larger in the NVE ensemble, compared to that with NPT.

Currently, I am equilibrating the 14mil atom sample with NVT ensemble for 50ps, as per the suggestion of Dr. Kohlmeyer. Once that finishes, I am going to equilibrate in NPT ensemble for 50ps and then analyze the potential energy convergence and volume convergence. If they are not yet converged, I will equilibrate for another 50ps in NPT. I will repeat this until I see adequate convergence, and only then try my hand at switching to the NVE ensemble.

After inspecting my previous NPT ensemble, I could easily see that the volume and potential energy were not yet converged, which would easily create large fluctuations under an instantaneous switch to NVE. Maybe I would run smaller systems, but currently my main computer to test this on is busy running the 14mil atom case, and I don’t want to disturb it.

I will let you guys know the results once they come in. Thanks again for all your help.

Jackson

8

I know that my KE + PE should remain constant throughout the NPT simulation, However, it “fluctuates” around being constant. Let me explain.

The PE is -6.8e7 and KE is 5.75e5 at time step one

If I look at timestep 2, the additions do not equal that at time step one, but the different is only a few 100 eV. When the difference is only a few 100 eV, and the total PE and KE is order of magnitudes higher, is this a reasonable error?

Thanks,

Jackson

9

I know that my KE + PE should remain constant throughout the NPT
simulation, However, it "fluctuates" around being constant. Let me explain.

not quite.

"KE + PE = const" only applies to the NVE ensemble. but since you
solve newton's equation of motion by discretizing them, even then the
numerical solution of the discretized differential equations will have
some fluctuation around the average (which depends on the length of
the time step). in fact, one the reasons why the verlet family of
algorithms is popular in MD (and less so elsewhere) is that its long
time error behavior is better than for higher order integrators (which
in turn have better accuracy on the short time scale).

for NVT, the conserved quantity is KE + PE + E_thermostat (i.e. the
net amount of energy exchanged between your system and the heat bath,
which can be accessed as global scalar from fix nvt)

for NPT, you have in addition to include the work against the pressure
(plus the net exchange of the barostat).

The PE is -6.8e7 and KE is 5.75e5 at time step one

If I look at timestep 2, the additions do not equal that at time step one,
but the different is only a few 100 eV. When the difference is only a few
100 eV, and the total PE and KE is order of magnitudes higher, is this a
reasonable error?

please have a look at a good MD text book, where various sources of
errors and how to assess them should be explained.

axel.

10

The change on one time step should be small, but it will never be zero. In the end, you have to figure out how big you can make the time step without changing the outcome of the simulation significantly. As the scientist, the onus is on you to figure out the meaning of significantly

11

Axel, I have been carrying out the procedure you suggested. First I equilibrated in NVT for 50ps. Then I switched to NPT and have been equilibrating for 40 ps so far. I attached the volume change in the unit cell as well as PE as a function of time. In my experience, I am a little surprised that both PE and Volume are still fluctuating after 40ps for such a large system (14 million atoms). Would you consider this to be reasonable? It might be that I need to equilibrate for a full ns, but I had previously considered that to be overkill for a 14 mil atom system.

Jackson

PE.PNG

volume.PNG

12

Axel, I have been carrying out the procedure you suggested. First I
equilibrated in NVT for 50ps. Then I switched to NPT and have been
equilibrating for 40 ps so far. I attached the volume change in the unit
cell as well as PE as a function of time. In my experience, I am a little
surprised that both PE and Volume are still fluctuating after 40ps for such
a large system (14 million atoms). Would you consider this to be reasonable?
It might be that I need to equilibrate for a full ns, but I had previously
considered that to be overkill for a 14 mil atom system.

why? the larger the system, the longer the time it takes to
equilibrate it, since a larger system can accommodate lower frequency
fluctuations.

what i prefer to do for large bulk systems, is to equilibrate a
smaller system first (so it is fast) and then use the replicate
command to enlarge the system until the finite size effects are gone.
if needed, this can be repeated.

a few 10s of picoseconds is not a lot of time...

axel.

13

If I replicate a smaller structure, would I not need to equilibrate the larger sample for at least some period of time? I guess 10s of ps is not that long, but it feels long when the NVT + NPT simulation takes 96 hours and does not seem yet close to being equilibrated.

Jackson

14

If I replicate a smaller structure, would I not need to equilibrate the
larger sample for at least some period of time? I guess 10s of ps is not

yes, of course, but you will start from a position much closer to
equilibrium for a fraction of the computational effort.

axel.