LAMMPS WWW Site - LAMMPS Documentation - LAMMPS Mailing List Archives
Re: [lammps-users] NPT with anisotropic cell fluctuations: conserved quantity and ensemble sampling
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: [lammps-users] NPT with anisotropic cell fluctuations: conserved quantity and ensemble sampling


From: Steven Vandenbrande <Steven.Vandenbrande@...6377...>
Date: Mon, 26 Mar 2018 14:48:33 +0000

On 09-03-18 17:25, Thompson, Aidan wrote:

Hi Steven,

 

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

 

1. Drift in the conserved quantity (large effect)

2. 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.  

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 some 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).

 

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. 

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.

 

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.

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

 

Aidan

 

 

--

      Aidan P. Thompson

      01444 Multiscale Science

      Sandia National Laboratories

      PO Box 5800, MS 1322      Phone: 505-844-9702

      Albuquerque, NM 87185     Fax  : 505-845-7442

      E-mail:athomps@...3... Cell : 505-218-1011

 

From: Giacomo Fiorin <giacomo.fiorin@...24...>
Date: Friday, March 9, 2018 at 7:25 AM
To: Steven Vandenbrande <Steven.Vandenbrande@...6381...>
Cc: Aidan Thompson <aidan.thompson@...24...>, lammps/lammps LAMMPS Users List <lammps-users@lists.sourceforge.net>
Subject: [EXTERNAL] Re: [lammps-users] NPT with anisotropic cell fluctuations: conserved quantity and ensemble sampling

 

 

On Fri, Mar 9, 2018 at 3:54 AM, Steven Vandenbrande <Steven.Vandenbrande@...6381...> wrote:

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

--
Steven Vandenbrande
Center for Molecular Modeling
Ghent University
Technologiepark 903,
B9052 Zwijnaarde
Belgium
Tel:
+32 9 264 66 37
E-mail: Steven.Vandenbrande@...6377...
http://molmod.UGent.be/


------------------------------------------------------------------------------
Check out the vibrant tech community on one of the world's most
engaging tech sites, Slashdot.org!
http://sdm.link/slashdot
_______________________________________________
lammps-users mailing list
lammps-users@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/lammps-users




--

Giacomo Fiorin

Associate Professor of Research, Temple University, Philadelphia, PA

Contractor, National Institutes of Health, Bethesda, MD


#############################
##### General Settings ######
#############################
log             log.system append
units           real
atom_style      full
boundary        p p p


#############################
##### PES Settings ##########
#############################

lattice         hcp 4.0
region          box block 0 5 0 5 0 5
create_box      1 box
create_atoms    1 box
mass            1 16.0

pair_style      lj/cut 12.500

dielectric      1.0
pair_modify     tail yes mix arithmetic
pair_coeff      1 1 0.294072658 3.73

#############################
##### Sampling Settings #####
#############################
variable        dt          equal 2.50
variable        pdamp       equal 1000*${dt}
variable        tdamp       equal 100*${dt}
variable        temperature equal 140.0
variable        pressure    equal 200.0
timestep        ${dt}
fix             1 all npt temp ${temperature} ${temperature} ${tdamp} aniso ${pressure} ${pressure} ${pdamp} mtk yes
fix_modify      1 energy yes # Add barostat/thermostat corrections so we get the conserved quantity


#############################
#### Equilibration ##########
#############################
thermo          0
run             400000

#############################
##### Output Settings ####### 
#############################
thermo_style 	custom step time etotal ke pe emol epair temp pxx pyy pzz pxy pxz pyz cella cellb cellc cellalpha cellbeta cellgamma vol press
variable ra_etotal equal etotal
variable ra_temp equal temp
variable ra_press equal press
variable ra_volume equal vol
fix runningavgs all ave/time 1 1000 1000 v_ra_etotal v_ra_temp v_ra_press v_ra_volume ave running file running_averages.out format %20.10f

#############################
#### Run ####################
#############################

thermo          100
run             1600000
unfix           1
#############################
##### General Settings ######
#############################
log             log.system append
units           real
atom_style      full
boundary        p p p


#############################
##### PES Settings ##########
#############################

lattice         hcp 4.0
region          box block 0 5 0 5 0 5
create_box      1 box
create_atoms    1 box
mass            1 16.0

pair_style      lj/cut 12.500

dielectric      1.0
pair_modify     tail yes mix arithmetic
pair_coeff      1 1 0.294072658 3.73

#############################
##### Sampling Settings #####
#############################
variable        dt          equal 2.50
variable        pdamp       equal 1000*${dt}
variable        tdamp       equal 100*${dt}
variable        temperature equal 140.0
variable        pressure    equal 500.0
timestep        ${dt}
fix             1 all npt temp ${temperature} ${temperature} ${tdamp} aniso ${pressure} ${pressure} ${pdamp} mtk yes
fix_modify      1 energy yes # Add barostat/thermostat corrections so we get the conserved quantity


#############################
#### Equilibration ##########
#############################
thermo          0
run             400000

#############################
##### Output Settings ####### 
#############################
thermo_style 	custom step time etotal ke pe emol epair temp pxx pyy pzz pxy pxz pyz cella cellb cellc cellalpha cellbeta cellgamma vol press
variable ra_etotal equal etotal
variable ra_temp equal temp
variable ra_press equal press
variable ra_volume equal vol
fix runningavgs all ave/time 1 1000 1000 v_ra_etotal v_ra_temp v_ra_press v_ra_volume ave running file running_averages.out format %20.10f

#############################
#### Run ####################
#############################

thermo          100
run             1600000
unfix           1
#############################
##### General Settings ######
#############################
log             log.system append
units           real
atom_style      full
boundary        p p p


#############################
##### PES Settings ##########
#############################

lattice         hcp 4.0
region          box block 0 5 0 5 0 5
create_box      1 box
create_atoms    1 box
mass            1 16.0

pair_style      lj/cut 12.500

dielectric      1.0
pair_modify     tail yes mix arithmetic
pair_coeff      1 1 0.294072658 3.73

#############################
##### Sampling Settings #####
#############################
variable        dt          equal 2.50
variable        pdamp       equal 1000*${dt}
variable        tdamp       equal 100*${dt}
variable        temperature equal 140.0
variable        pressure    equal 200.0
timestep        ${dt}
fix             1 all npt temp ${temperature} ${temperature} ${tdamp} iso ${pressure} ${pressure} ${pdamp} mtk yes
fix_modify      1 energy yes # Add barostat/thermostat corrections so we get the conserved quantity


#############################
#### Equilibration ##########
#############################
thermo          0
run             400000

#############################
##### Output Settings ####### 
#############################
thermo_style 	custom step time etotal ke pe emol epair temp pxx pyy pzz pxy pxz pyz cella cellb cellc cellalpha cellbeta cellgamma vol press
variable ra_etotal equal etotal
variable ra_temp equal temp
variable ra_press equal press
variable ra_volume equal vol
fix runningavgs all ave/time 1 1000 1000 v_ra_etotal v_ra_temp v_ra_press v_ra_volume ave running file running_averages.out format %20.10f

#############################
#### Run ####################
#############################

thermo          100
run             1600000
unfix           1
#############################
##### General Settings ######
#############################
log             log.system append
units           real
atom_style      full
boundary        p p p


#############################
##### PES Settings ##########
#############################

lattice         hcp 4.0
region          box block 0 5 0 5 0 5
create_box      1 box
create_atoms    1 box
mass            1 16.0

pair_style      lj/cut 12.500

dielectric      1.0
pair_modify     tail yes mix arithmetic
pair_coeff      1 1 0.294072658 3.73

#############################
##### Sampling Settings #####
#############################
variable        dt          equal 2.50
variable        pdamp       equal 1000*${dt}
variable        tdamp       equal 100*${dt}
variable        temperature equal 140.0
variable        pressure    equal 500.0
timestep        ${dt}
fix             1 all npt temp ${temperature} ${temperature} ${tdamp} iso ${pressure} ${pressure} ${pdamp} mtk yes
fix_modify      1 energy yes # Add barostat/thermostat corrections so we get the conserved quantity


#############################
#### Equilibration ##########
#############################
thermo          0
run             400000

#############################
##### Output Settings ####### 
#############################
thermo_style 	custom step time etotal ke pe emol epair temp pxx pyy pzz pxy pxz pyz cella cellb cellc cellalpha cellbeta cellgamma vol press
variable ra_etotal equal etotal
variable ra_temp equal temp
variable ra_press equal press
variable ra_volume equal vol
fix runningavgs all ave/time 1 1000 1000 v_ra_etotal v_ra_temp v_ra_press v_ra_volume ave running file running_averages.out format %20.10f

#############################
#### Run ####################
#############################

thermo          100
run             1600000
unfix           1

Attachment: md_hcp_lj_aniso_t0-p0.png
Description: md_hcp_lj_aniso_t0-p0.png

Attachment: md_hcp_lj_iso_t0-p0.png
Description: md_hcp_lj_iso_t0-p0.png

Attachment: shirts_hcp_lj_aniso_press.png
Description: shirts_hcp_lj_aniso_press.png

Attachment: shirts_hcp_lj_iso_press.png
Description: shirts_hcp_lj_iso_press.png