Thermostatting tutorial (#89)

Example of running MD with different thermostats, discussing sampling efficiency and impact on dynamics. Based on both i-PI and LAMMPS. --------- Co-authored-by: Davide Tisi <[email protected]>
lab-cosmo · Oct 14, 2024 · 6d4d1a9 · 6d4d1a9
1 parent 791e47c
commit 6d4d1a9
Show file tree

Hide file tree

Showing 21 changed files with 2,591 additions and 2 deletions.
diff --git a/docs/src/conf.py b/docs/src/conf.py
@@ -70,7 +70,7 @@ def post_build_tweaks(app, exception):
 
     # adds custom urls to the sitemap
     custom_urls = [
-        {"loc": "https://atomistic-cookbook.org", "priority": 2.0},
+        {"loc": "https://atomistic-cookbook.org", "priority": 1.0},
     ]
 
     sitemap_file = os.path.join(build_folder, "sitemap.xml")

diff --git a/docs/src/software/chemiscope.sec b/docs/src/software/chemiscope.sec
@@ -16,3 +16,4 @@ repository <https://github.com/lab-cosmo/chemiscope>`_.
 - examples/gaas-map/gaas-map
 - examples/pi-metad/pi-metad
 - examples/path-integrals/path-integrals
+- examples/thermostatting/thermostatting
diff --git a/docs/src/software/i-pi.sec b/docs/src/software/i-pi.sec
@@ -8,7 +8,7 @@ it on the `ipi-code website <http://ipi-code.org>`_, the `documentation pages
 <http://ipi-code.org/i-pi>`_ or `the github repository
 <https://github.com/i-pi/i-pi>`_.
 
-
+- examples/thermostatting/thermostatting
 - examples/path-integrals/path-integrals
 - examples/pi-metad/pi-metad
 - examples/heat-capacity/heat-capacity
diff --git a/docs/src/software/lammps.sec b/docs/src/software/lammps.sec
@@ -7,5 +7,6 @@ advanced molecular simulations techniques, and is highly parallelized using an
 efficient domain decomposition scheme. Learn more about LAMMPS on its `homepage
 <https://www.lammps.org>`_.
 
+- examples/thermostatting/thermostatting
 - examples/path-integrals/path-integrals
 - examples/heat-capacity/heat-capacity
diff --git a/docs/src/topics/sampling.sec b/docs/src/topics/sampling.sec
@@ -5,6 +5,7 @@ This section contains recipes that compute thermodynamic averages by sampling,
 evaluates dynamical properties, or otherwise computes the properties of a
 set of configurations of an atomistic system.
 
+- examples/thermostatting/thermostatting
 - examples/path-integrals/path-integrals
 - examples/pi-metad/pi-metad
 - examples/batch-cp2k/reference-trajectory

diff --git a/examples/thermostats/README.rst b/examples/thermostats/README.rst
@@ -0,0 +1,9 @@
+Thermostatting molecular dynamics simulations, a primer
+======================================================
+
+This example provides some practical guidelines to set up and 
+run simulations using different types of (primarily stochastic)
+thermostats. It discusses both `i-PI <http://ipi-code.org>`_ and
+`LAMMPS <http://lammps.org>`_ as practical codes, and it 
+analizes and visualize the output using ``matplotlib`` and 
+`chemiscope <https://chemiscope.org>`_.
diff --git a/examples/thermostats/data/gle.lmp b/examples/thermostats/data/gle.lmp
@@ -0,0 +1,29 @@
+units		electron
+atom_style	full
+
+pair_style      lj/cut/tip4p/long 1 2 1 1 0.278072379 17.007
+bond_style      class2 
+angle_style     harmonic
+kspace_style	pppm/tip4p 0.0001
+
+read_data	data/water_32_data.lmp
+pair_coeff  * * 0 0
+pair_coeff  1  1  0.000295147 5.96946
+
+neighbor	2.0 bin
+
+velocity all create 300.0 2345187
+
+dump		1 all xyz 100 lammps_pos.xyz
+thermo     100
+thermo_style custom step temp press pe etotal
+variable T equal temp
+variable P equal press
+variable PE equal pe
+variable ETOTAL equal etotal
+fix thermo_out all ave/time 10 1 10 v_T v_PE v_ETOTAL v_P file lammps_out.dat
+
+timestep	1.0
+fix 1 all gle 6 300 300 31415 data/smart.A
+run		10000
+
diff --git a/examples/thermostats/data/in.lmp b/examples/thermostats/data/in.lmp
@@ -0,0 +1,30 @@
+units		electron
+atom_style	full
+
+pair_style      lj/cut/tip4p/long 1 2 1 1 0.278072379 17.007
+bond_style      class2 
+angle_style     harmonic
+kspace_style	pppm/tip4p 0.0001
+
+read_data	data/water_32_data.lmp
+pair_coeff  * * 0 0
+pair_coeff  1  1  0.000295147 5.96946
+
+neighbor	2.0 bin
+
+timestep	0.00025
+
+#velocity all create 298.0 2345187
+
+#thermo_style	multi
+#thermo		1
+
+#fix		1 all nvt temp 298.0 298.0 30.0 tchain 1
+#fix 1 all nve
+fix 1 all ipi h2o-lammps 32342 unix
+
+
+#dump		1 all xyz 25 dump.xyz
+
+run		100000000
+
diff --git a/examples/thermostats/data/input_cvv_sample.xml b/examples/thermostats/data/input_cvv_sample.xml
@@ -0,0 +1,45 @@
+<simulation verbosity='medium' safe_stride='100'>
+  <output prefix='sample_nvt'>
+      <properties stride='1' filename='out'>  
+      [ step, time{picosecond}, conserved{electronvolt}, 
+      temperature{kelvin}, kinetic_md{electronvolt}, potential{electronvolt}, 
+              temperature(H){kelvin}, temperature(O){kelvin} ] </properties>
+      <checkpoint filename='chk' stride='2000' overwrite='False'/>
+  </output>
+  <total_steps> 20000 </total_steps>
+  <prng>
+    <seed> 32342 </seed>
+  </prng>
+  <ffsocket name='lmpserial' mode='unix' pbc='false'>
+      <address>h2o-lammps</address> <latency> 1e-4 </latency>
+  </ffsocket>
+  <system>
+    <initialize nbeads='1'>
+      <file mode='pdb' units='angstrom'> data/water_32.pdb </file>
+      <velocities mode='thermal' units='kelvin'> 300 </velocities>
+    </initialize>
+    <forces>
+      <force forcefield='lmpserial'> lmpserial </force>
+    </forces>
+    <ensemble>
+      <temperature units='kelvin'>300</temperature>
+    </ensemble>
+    <motion mode='dynamics'>
+      <dynamics mode='nvt'>
+        <timestep units='femtosecond'> 1.0 </timestep>
+  <thermostat mode='gle'>
+    <A shape='(7,7)'>
+      [   8.191023526179e-4,    8.328506066524e-3,    1.657771834013e-3,    9.736989925341e-4,    2.841803794895e-4,   -3.176846864198e-5,   -2.967010478210e-4, 
+  -8.389856546341e-4,    2.405526974742e-2,   -1.507872374848e-2,    2.589784240185e-3,    1.516783633362e-3,   -5.958833418565e-4,    4.198422349789e-4, 
+   7.798710586406e-4,    1.507872374848e-2,    8.569039501219e-3,    6.001000899602e-3,    1.062029383877e-3,    1.093939147968e-3,   -2.661575532976e-3, 
+  -9.676783161546e-4,   -2.589784240185e-3,   -6.001000899602e-3,    2.680459336535e-5,   -5.214694469742e-5,    4.231304910751e-4,   -2.104894919743e-5, 
+  -2.841997149166e-4,   -1.516783633362e-3,   -1.062029383877e-3,    5.214694469742e-5,    1.433903506353e-9,   -4.241574212449e-5,    7.910178912362e-5, 
+   3.333208286893e-5,    5.958833418565e-4,   -1.093939147968e-3,   -4.231304910751e-4,    4.241574212449e-5,    2.385554468441e-8,   -3.139255482869e-5, 
+   2.967533789056e-4,   -4.198422349789e-4,    2.661575532976e-3,    2.104894919743e-5,   -7.910178912362e-5,    3.139255482869e-5,   2.432567259684e-11
+     ]
+    </A>
+  </thermostat>        
+      </dynamics>
+    </motion>
+  </system>
+</simulation>
diff --git a/examples/thermostats/data/input_cvv_traj.xml b/examples/thermostats/data/input_cvv_traj.xml
@@ -0,0 +1,48 @@
+<simulation verbosity='medium' safe_stride='100' threading="True">
+  <output prefix='traj'>
+      <properties stride='1' filename='out'>  
+      [ step, time{picosecond}, conserved{electronvolt}, 
+      temperature{kelvin}, kinetic_md{electronvolt}, potential{electronvolt}, 
+      temperature(H){kelvin}, temperature(O){kelvin} ] </properties>
+    <trajectory filename='pos' stride='10'> positions </trajectory>
+    <trajectory filename='vel' stride='1'> velocities </trajectory>
+  </output>
+  <total_steps> 2000 </total_steps>
+  <prng>
+    <seed> 32342 </seed>
+  </prng>
+  <ffsocket name='lmpserial' mode='unix' pbc='false'>
+      <address>h2o-lammps</address> <latency> 1e-4 </latency>
+  </ffsocket>
+  <system_template>
+  <labels> ['ITRAJ'] </labels>
+  <instance> [1] </instance>
+  <instance> [2] </instance>
+  <instance> [3] </instance>
+  <instance> [4] </instance>
+  <instance> [5] </instance>
+  <instance> [6] </instance>
+  <instance> [7] </instance>
+  <instance> [8 ] </instance>
+  <template>
+  <system prefix='traj-ITRAJ'>
+    <initialize nbeads='1'>
+      <file mode='chk'> sample_nvt.chk_ITRAJ </file>
+    </initialize>
+    <forces>
+      <force forcefield='lmpserial'> lmpserial </force>
+    </forces>
+    <ensemble>
+      <temperature units='kelvin'>300</temperature>
+    </ensemble>
+    <motion mode='dynamics'>
+      <dynamics mode='nve'>
+        <!-- this is still too long to be accurate but we use it
+        to be consistent with the setup in the example. -->
+        <timestep units='femtosecond'> 1.0 </timestep>
+      </dynamics>
+    </motion>
+  </system>
+  </template>
+  </system_template>
+</simulation>
diff --git a/examples/thermostats/data/input_gle.xml b/examples/thermostats/data/input_gle.xml
@@ -0,0 +1,46 @@
+<simulation verbosity='medium' safe_stride='100'>
+  <output prefix='simulation_gle'>
+      <properties stride='1' filename='out'>  
+      [ step, time{picosecond}, conserved{electronvolt}, 
+      temperature{kelvin}, kinetic_md{electronvolt}, potential{electronvolt}, 
+      temperature(H){kelvin}, temperature(O){kelvin} ] </properties>
+    <trajectory filename='pos' stride='10'> positions </trajectory>
+    <trajectory filename='vel' stride='1'> velocities </trajectory>
+  </output>
+  <total_steps> 2000 </total_steps>
+  <prng>
+    <seed> 32342 </seed>
+  </prng>
+  <ffsocket name='lmpserial' mode='unix' pbc='false'>
+      <address>h2o-lammps</address> <latency> 1e-4 </latency>
+  </ffsocket>
+  <system>
+    <initialize nbeads='1'>
+      <file mode='pdb' units='angstrom'> data/water_32.pdb </file>
+      <velocities mode='thermal' units='kelvin'> 300 </velocities>
+    </initialize>
+    <forces>
+      <force forcefield='lmpserial'> lmpserial </force>
+    </forces>
+    <ensemble>
+      <temperature units='kelvin'>300</temperature>
+    </ensemble>
+    <motion mode='dynamics'>
+      <dynamics mode='nvt'>
+        <timestep units='femtosecond'> 1.0 </timestep>
+  <thermostat mode='gle'>
+    <A shape='(7,7)'>
+      [   8.191023526179e-4,    8.328506066524e-3,    1.657771834013e-3,    9.736989925341e-4,    2.841803794895e-4,   -3.176846864198e-5,   -2.967010478210e-4, 
+  -8.389856546341e-4,    2.405526974742e-2,   -1.507872374848e-2,    2.589784240185e-3,    1.516783633362e-3,   -5.958833418565e-4,    4.198422349789e-4, 
+   7.798710586406e-4,    1.507872374848e-2,    8.569039501219e-3,    6.001000899602e-3,    1.062029383877e-3,    1.093939147968e-3,   -2.661575532976e-3, 
+  -9.676783161546e-4,   -2.589784240185e-3,   -6.001000899602e-3,    2.680459336535e-5,   -5.214694469742e-5,    4.231304910751e-4,   -2.104894919743e-5, 
+  -2.841997149166e-4,   -1.516783633362e-3,   -1.062029383877e-3,    5.214694469742e-5,    1.433903506353e-9,   -4.241574212449e-5,    7.910178912362e-5, 
+   3.333208286893e-5,    5.958833418565e-4,   -1.093939147968e-3,   -4.231304910751e-4,    4.241574212449e-5,    2.385554468441e-8,   -3.139255482869e-5, 
+   2.967533789056e-4,   -4.198422349789e-4,    2.661575532976e-3,    2.104894919743e-5,   -7.910178912362e-5,    3.139255482869e-5,   2.432567259684e-11
+     ]
+    </A>
+  </thermostat>        
+      </dynamics>
+    </motion>
+  </system>
+</simulation>
diff --git a/examples/thermostats/data/input_higamma.xml b/examples/thermostats/data/input_higamma.xml
@@ -0,0 +1,37 @@
+<simulation verbosity='medium' safe_stride='100'>
+  <output prefix='simulation_higamma'>
+      <properties stride='1' filename='out'>  
+      [ step, time{picosecond}, conserved{electronvolt}, 
+      temperature{kelvin}, kinetic_md{electronvolt}, potential{electronvolt}, 
+      temperature(H){kelvin}, temperature(O){kelvin} ] </properties>
+    <trajectory filename='pos' stride='10'> positions </trajectory>
+    <trajectory filename='vel' stride='1'> velocities </trajectory>
+  </output>
+  <total_steps> 2000 </total_steps>
+  <prng>
+    <seed> 32342 </seed>
+  </prng>
+  <ffsocket name='lmpserial' mode='unix' pbc='false'>
+      <address>h2o-lammps</address> <latency> 1e-4 </latency>
+  </ffsocket>
+  <system>
+    <initialize nbeads='1'>
+      <file mode='pdb' units='angstrom'> data/water_32.pdb </file>
+      <velocities mode='thermal' units='kelvin'> 300 </velocities>
+    </initialize>
+    <forces>
+      <force forcefield='lmpserial'> lmpserial </force>
+    </forces>
+    <ensemble>
+      <temperature units='kelvin'>300</temperature>
+    </ensemble>
+    <motion mode='dynamics'>
+      <dynamics mode='nvt'>
+        <timestep units='femtosecond'> 1.0 </timestep>
+        <thermostat mode='langevin'>
+          <tau units='femtosecond'> 10 </tau>
+        </thermostat>
+      </dynamics>
+    </motion>
+  </system>
+</simulation>
diff --git a/examples/thermostats/data/input_nve.xml b/examples/thermostats/data/input_nve.xml
@@ -0,0 +1,34 @@
+<simulation verbosity='medium' safe_stride='100'>
+  <output prefix='simulation_nve'>
+      <properties stride='1' filename='out'>  
+      [ step, time{picosecond}, conserved{electronvolt}, 
+      temperature{kelvin}, kinetic_md{electronvolt}, potential{electronvolt}, 
+      temperature(H){kelvin}, temperature(O){kelvin} ] </properties>
+    <trajectory filename='pos' stride='10'> positions </trajectory>
+    <trajectory filename='vel' stride='1'> velocities </trajectory>
+  </output>
+  <total_steps> 2000 </total_steps>
+  <prng>
+    <seed> 32342 </seed>
+  </prng>
+  <ffsocket name='lmpserial' mode='unix' pbc='false'>
+      <address>h2o-lammps</address> <latency> 1e-4 </latency>
+  </ffsocket>
+  <system>
+    <initialize nbeads='1'>
+      <file mode='pdb' units='angstrom'> data/water_32.pdb </file>
+      <velocities mode='thermal' units='kelvin'> 300 </velocities>
+    </initialize>
+    <forces>
+      <force forcefield='lmpserial'> lmpserial </force>
+    </forces>
+    <ensemble>
+      <temperature units='kelvin'>300</temperature>
+    </ensemble>
+    <motion mode='dynamics'>
+      <dynamics mode='nve'>
+        <timestep units='femtosecond'> 1.0 </timestep>
+      </dynamics>
+    </motion>
+  </system>
+</simulation>
diff --git a/examples/thermostats/data/input_svr.xml b/examples/thermostats/data/input_svr.xml
@@ -0,0 +1,37 @@
+<simulation verbosity='medium' safe_stride='100'>
+  <output prefix='simulation_svr'>
+      <properties stride='1' filename='out'>  
+      [ step, time{picosecond}, conserved{electronvolt}, 
+      temperature{kelvin}, kinetic_md{electronvolt}, potential{electronvolt}, 
+      temperature(H){kelvin}, temperature(O){kelvin} ] </properties>
+    <trajectory filename='pos' stride='10'> positions </trajectory>
+    <trajectory filename='vel' stride='1'> velocities </trajectory>
+  </output>
+  <total_steps> 2000 </total_steps>
+  <prng>
+    <seed> 32342 </seed>
+  </prng>
+  <ffsocket name='lmpserial' mode='unix' pbc='false'>
+      <address>h2o-lammps</address> <latency> 1e-4 </latency>
+  </ffsocket>
+  <system>
+    <initialize nbeads='1'>
+      <file mode='pdb' units='angstrom'> data/water_32.pdb </file>
+      <velocities mode='thermal' units='kelvin'> 300 </velocities>
+    </initialize>
+    <forces>
+      <force forcefield='lmpserial'> lmpserial </force>
+    </forces>
+    <ensemble>
+      <temperature units='kelvin'>300</temperature>
+    </ensemble>
+    <motion mode='dynamics'>
+      <dynamics mode='nvt'>
+        <timestep units='femtosecond'> 1.0 </timestep>
+        <thermostat mode='svr'>
+          <tau units='femtosecond'> 10 </tau>
+        </thermostat>
+      </dynamics>
+    </motion>
+  </system>
+</simulation>
diff --git a/examples/thermostats/data/run_traj.sh b/examples/thermostats/data/run_traj.sh
@@ -0,0 +1,25 @@
+#!/bin/bash
+
+# This is a simple script to run some reference
+# NVE trajectories for the autocorrelation function,
+# starting from samples collected from a NVT trajectory
+
+# First, we launch the sampling trajectory. This 
+# outputs checkpoint files every 2 ps.
+
+i-pi data/input_cvv_sample.xml & sleep 4
+lmp < data/in.lmp
+
+# Then, we run an i-PI input that launches multiple
+# NVE trajectories
+
+i-pi data/input_cvv_traj.xml & sleep 4
+for i in {1..8}; do lmp < data/in.lmp & done
+wait
+
+# Finally, run the i-PI post-processing on the 
+# concatenated velocity trajectories. The block size
+# is chosen so each trajectory is a separate block.
+
+cat traj-*_traj.vel_0.xyz &> traj-all.xyz
+i-pi-getacf -ifile traj-all.xyz -mlag 1000 -bsize 2001 -ftpad 2000 -ftwin cosine-blackman -dt "1 femtosecond" -oprefix traj-all
diff --git a/examples/thermostats/data/smart.A b/examples/thermostats/data/smart.A
@@ -0,0 +1,7 @@
+   3.386281653149e-2    3.443118824048e-1    6.853456504773e-2    4.025405400887e-2    1.174840729213e-2   -1.313352136850e-3   -1.226602892172e-2 
+  -3.468481955231e-2    9.944778982366e-1   -6.233751547554e-1    1.070652383078e-1    6.270591914470e-2   -2.463463597066e-2    1.735685476891e-2 
+   3.224094094292e-2    6.233751547554e-1    3.542558650374e-1    2.480896213019e-1    4.390575373437e-2    4.522494721916e-2   -1.100331889779e-1 
+  -4.000515097620e-2   -1.070652383078e-1   -2.480896213019e-1    1.108138713594e-3   -2.155826332720e-3    1.749279579329e-2   -8.701924766486e-4 
+  -1.174920664525e-2   -6.270591914470e-2   -4.390575373437e-2    2.155826332720e-3    5.927954083432e-8   -1.753525049730e-3    3.270176631583e-3 
+   1.377994097069e-3    2.463463597066e-2   -4.522494721916e-2   -1.749279579329e-2    1.753525049730e-3    9.862209897519e-7   -1.297811343382e-3 
+   1.226819236065e-2   -1.735685476891e-2    1.100331889779e-1    8.701924766486e-4   -3.270176631583e-3    1.297811343382e-3    1.005656723509e-9