NPT with anisotropic cell fluctuations: conserved quantity and ensemble sampling

Steven_Vandenbrande · March 9, 2018, 8:54am

I am running into problems when simulating in the NPT ensemble with
anisotropic cell fluctuations are allowed. In short, when only
considering isotropic cell variations (using the iso keyword), the
conserved quantity is nearly constant and the volume fluctuations seem
correct. As soon as anisotropic cell variations are allowed (using for
example the aniso or tri keywords), there is a clear drift in the
conserved quantity as well as indications that the ensemble is
incorrectly sampled.

The test system is a metal-organic framework, MOF-5 with an initially
cubic unit cell (I am aware that this is not the most simple test
system, any suggestions for a less complicated case are welcome). The
force field used is UFF4MOF (Lennard-Jones pair interactions, no
electrostatics are considered, standard covalent interaction terms). The
first simulation is done at 298K and an isotropic pressure of 1atm, with
the iso keyword to only allow isotropic fluctuations, with 1ns
equilibration and 4ns for data collection (time step 0.5ns). The plot
md_iso_t0-p0.png shows some quantities during the MD run in blue and the
corresponding running averages in red. Everything seems fine and there
is no systematic drift in the conserved quantity. I also performed a
second simulation at a pressure of 500atm and applied the method of
Shirts (10.1021/ct300688p) to check the ensemble sampling. As shown in
shirts_iso_press.png, the analytical slope corresponds to the calculated
one from the simulation. Up to know everything seems to work perfectly.

Things go haywire when allowing anisotropic cell fluctuations (using for
instance the aniso or tri keywords). The plots of the same simulations
as before (where only the iso keyword is changed, the applied pressure
is still isotropic) show that there is a clear drift in the conserved
quantity (5kcal/mol/ns for aniso, more than 10kcal/mol/ns for tri).
Additionally, by applying the method of Shirts I see a significant
deviation between the analytical and simulated ratio of volume
distributions, which indicates that the wrong ensemble is sampled.

I have tried changing most of the sampling settings (timestep, damping
parameters of thermo/barostat, number of thermostat chains for
particles/barostat particles, using the mtk term or not), but the
conclusions are always the same: for isotropic cell fluctuations
everything seems fine, but things go wrong as soon as anisotropic cell
fluctuations are allowed.

I have looked into the source code, but could not pinpoint the problem.
Some possible problems I encountered:
1) The integration scheme by Tuckerman (10.1088/0305-4470/39/19/S18),
which is said to be followed, does not provide equations for anisotropic
fluctuations. These are given by Yu (10.1016/j.chemphys.2010.02.014). A
difference between both cases is the number of degrees of freedom for
the barostat (1 for iso, 3 for aniso, 6 for tri), but it seems this is
not used in the source code.
2) A separate Nose-Hoover thermostat is used for the barostat particles,
which acts first and last in the integration scheme. I do not understand
fullly why the Liouville operator is written in that way.
3) The derivations by Tuckerman always assume symmetric cell tensors.
The LAMMPS convention is however different, perhaps this leads to some
changes? (Although this would not explain the problems with the aniso
case, where for the test system the cell tensor is diagonal).

Lammps version is from the master branch on 07/03/2018

Does anybody have any suggestions how the mentioned problems could be
solved?

Regards

Steven Vandenbrande

data.system (171 KB)

in.system.aniso (1.98 KB)

in.system.iso (1.98 KB)

in.system.tri (1.98 KB)

Giacomo_Fiorin · March 9, 2018, 2:25pm

Hi Steven, there are some solid-crystalline examples in the LAMMPS distribution that you can adapt, for example adapting the silicon example used in the ELASTIC_T test, which is currently set up for an NVT simulation (Langevin thermostat). Aidan Thompson may be able to comment on whether that is a suitable example, and on the MTK implementation in LAMMPS.

Have you also compared the starting and ending coordinates with each ensemble? I am not familiar with the potential that you’re using, but for sure it is worth checking that out.

Giacomo

athomps · March 9, 2018, 4:25pm

Hi Steven,

Thanks for the careful analysis. There seem to be two distinct problems that you have identified in going from iso to aniso:

Drift in the conserved quantity (large effect)
Distortion of the volume distribution (small effect)

Looking at the results in the paper by Yu et al., they do not see the first problem, but the second problem might still be there, because they did not look at the volume distribution very closely. Moreover, it is well known (but not much talked about) that the non-Hamiltonian equations of motion do not reliably sample strain fluctuations in solids. This seems to be an issue with the system getting trapped in very long-lived orbits.

The bottom line is, we might be able to eliminate problem 1 by tweaking something in the code, using the Yu paper for inspiration. However, this will probably not affect problem 2.

If you could follow Giacomo’s suggestion, and construct a simpler demo of Problems 1 and 2, perhaps using the ELASTIC_T test, or else even LJ, I can take a closer look.

Aidan

Steven_Vandenbrande · March 10, 2018, 12:04pm

I have already performed simulations of a very similar system using our in-house MD code. Both for the isotropic and anisotropic case, it satisfies the test proposed by Shirts, providing evidence that the volume distribution is correctly sampled. FYI: this MD code also uses the MTK barostat, but a different expansion of the Liouville operator. There are also some other conceptual differences (for instance a symmetric cell tensor is used). I am not sure the two problems are independent. If there is something wrong with the implementation of the equations of motion, it is possible that this affects both the conserved quantity and the ensemble sampling. Again, I am not sure, but we should not yet rule out the possibility that both problems are connected. I was able to reproduce the same problems for a simulation of solid Lennard-Jones in the hcp form. Input files and resulting plots are attached. This gives an indication that the problem is not related to the potential energy surface, but to the sampling. Thanks for your help. Steven

in.system.aniso.t0-p0 (1.81 KB)

in.system.aniso.t0-p2 (1.81 KB)

in.system.iso.t0-p0 (1.81 KB)

in.system.iso.t0-p2 (1.81 KB)

Steven_Vandenbrande · March 26, 2018, 2:48pm

I have already performed simulations of a very similar system using our in-house MD code. Both for the isotropic and anisotropic case, it satisfies the test proposed by Shirts, providing evidence that the volume distribution is correctly sampled. FYI: this MD code also uses the MTK barostat, but a different expansion of the Liouville operator. There are also some other conceptual differences (for instance a symmetric cell tensor is used). I am not sure the two problems are independent. If there is something wrong with the implementation of the equations of motion, it is possible that this affects both the conserved quantity and the ensemble sampling. Again, I am not sure, but we should not yet rule out the possibility that both problems are connected. I was able to reproduce the same problems for a simulation of solid Lennard-Jones in the hcp form. Input files and resulting plots are attached. This gives an indication that the problem is not related to the potential energy surface, but to the sampling. Thanks for your help. Steven

in.system.aniso.t0-p0 (1.81 KB)

in.system.aniso.t0-p2 (1.81 KB)

in.system.iso.t0-p0 (1.81 KB)

in.system.iso.t0-p2 (1.81 KB)

Steven_Vandenbrande · April 4, 2018, 3:30pm

Dear LAMMPS users,

To further investigate the problems with the sampling of anisotropic cells, I implemented another version of the MTTK barostat with NHC temperature control in LAMMPS.

The results I report here are for fcc Lennard-Jones argon (at 60K and 1bar), inspired by the following paper: https://doi.org/10.1080/08927022.2017.1313418
In this paper, elastic constants of this system were computed by looking at cell fluctuations in NPT simulations using LAMMPS. It was concluded that this method is not reliable, as the results deviate from other approaches (such as MC) and strongly depend on the barostat relaxation time. This is shown in Figure 2 of the paper and my reproduction is attached as C_vs_taup_orig.png (symbols: NPT MD, full lines NPT MC).

When I use my new implementation of the NPT ensemble, things look much better (see C_vs_taup_newbarostat.png). Except for small values of the barostat relaxation time, the elastic constants are constant and quite close to the MC results.

This confirms that there is something wrong with the sampling of the NPT ensemble in LAMMPS. Additionally, these results shows that there can be important consequences when studying fluctuations.
The fact that things are fixed with another implementation, indicate that the problems are not related to the PES or due to issues with the system getting trapped in very long-lived orbits, as suggested below.
Are there developers who could help resolve these issues? Other users who experienced similar problems?

Some technical notes:
    the new version of the barostat is the same as used in YAFF (https://github.com/molmod/yaff). It uses a different expansion of the Liouville operator than the current LAMMPS implementation. Unfortunately, this expansion requires evaluation of the forces 3 times per time step, which is a significant additional computional cost.
    the LAMMPS version with the new barostat can be found here: https://github.com/stevenvdb/lammps (``newbarostat'' branch). It is not quite ready for a pull-request yet...
    example input files for tau_p = 1ps are attached

Regards

Steven Vandenbrande

data.system (80.1 KB)

in_newbarostat.system (2.08 KB)

in_original.system (2.16 KB)

athomps · April 5, 2018, 8:49pm

Hi Steven,

This is very interesting. The problem with Nose-Hoover barostat and elastic constants is not unique to LAMMPS and has been rediscoved in the literature every few years, going back to Sprik, Impey, and Klein, PRB 29 4368 (1984). This is the first time I have heard of any integration scheme which does not exhibit the problem, so that is very intriguing. It would be interesting to know if there is a way to rearrange the Trotter factorization of the Liouville operator so that the force need only be evaluated once per timestep, without giving up the correct sampling of strain.

Aidan

Steven_Vandenbrande · April 7, 2018, 7:34pm

In that paper they conclude that extracting elastic constants from NPT MD requires long sampling. They however only used 50 000 steps, which is easy to surpass with modern computational power. I must admit that my knowledge about relevant other papers is very limited, as this is not my area of research. I did another test for elastic constants of fcc LJ, this time using the HOOMD-blue code. This uses a two step algorithm (so the force is only evaluated once per timestep) and seems to provide the same results as LAMMPS with the new barostat, while the current LAMMPS barostat provides results that depend on the barostat relaxation time and deviate from HOOMD-blue. See attachment for input files and resulting plots. The system is very similar to , with two slight modifications to run it in HOOMD-blue: 1) no tail corrections 2) 7x7x7 instead of 5x5x5 unit cells.

in_newbaro.system (2.02 KB)

in_original.system (2.02 KB)

md_hoomdblue.py (1.17 KB)

Steven_Vandenbrande · April 10, 2018, 7:16am

Dear LAMMPS users,

I repeated the simulations mentioned below (extracting elastic constants from NPT MD for fcc LJ), this time using DL-POLY-4.
The results are similar to those obtained using HOOMD-blue and LAMMPS with my version of the barostat (see elastic_constants.png). The only outlier is now the current LAMMPS code.

It seems to me that there is a bug in the LAMMPS MTK barostat. It is probably something subtle, because equilibrium properties seem fine but fluctuations are incorrectly sampled.
This is shown in volume_dists.png: the mean of the volume distribution is correct but the width is incorrect (and depends on the barostat relaxation time).

If anyone with a good knowledge of that part of the LAMMPS code would be willing to debug this, I would be very thankful.

Regards,

Steven

DL-POLY-4 input files are attached. Note that I turned off long-range corrections in the source code, for consistency with earlier calculations.

CONFIG (130 KB)

CONTROL (389 Bytes)

FIELD (163 Bytes)

David_Nicholson · April 10, 2018, 7:27pm

Hello Steven,

I ran your script “in_original.system” and I noticed that there is a significant center of mass velocity for the system (around 122 m/s). I think this is likely to be significantly affecting your results, since “compute temp/com” registers a temperature substantially lower than the set value (see attached plots).

I find that I can fix this issue with a “velocity create” before the start of the run.

Best,
-David

in_original_vel.system (2.35 KB)

Steven_Vandenbrande1 · April 10, 2018, 7:51pm

Hi David,

The file you mentioned comes directly from the supporting information of https://doi.org/10.1080/08927022.2017.1313418
I used it to replicate as closely as possible those reported results and already noticed that this was an issue.

In all following simulations I issued “velocity all zero linear” (see input files in this post: ) after creating the velocities. Thanks for thinking along, I am not yet excluding the possibility that there is something wrong with my input files.

David_Nicholson · April 10, 2018, 8:59pm

Hi Steven,

My mistake; I didn’t see that you already caught it. Clearly this does not resolve the fluctuation issues, but is it possible that it was the source of the drift in the conserved quantity?

-David

Steven_Vandenbrande · April 11, 2018, 9:52am

There is still a drift, but for large barostat relaxation times it becomes very small.
The attached plot shows the running average of the conserved quantity during the production run of the simulation, offset to 0 at t=1ns.

Steven

Steven_Vandenbrande · December 20, 2018, 8:40am

It seems that these issues are fixed by this pull request: https://github.com/lammps/lammps/pull/1259
The attached plot shows that with this LAMMPS version the elastic constants are nearly independent of barostat relaxation time and now do correspond with other codes. The previously observed energy drift is gone as well.

Thanks for solving this!

Steven