diff --git a/docs/src/conf.py b/docs/src/conf.py index cd199aba..39ba8519 100644 --- 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 index 360bf713..766bfe7c 100644 --- a/docs/src/software/chemiscope.sec +++ b/docs/src/software/chemiscope.sec @@ -16,3 +16,4 @@ repository `_. - 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 index 210b860a..0940f475 100644 --- a/docs/src/software/i-pi.sec +++ b/docs/src/software/i-pi.sec @@ -8,7 +8,7 @@ it on the `ipi-code website `_, the `documentation pages `_ or `the github repository `_. - +- 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 index 3da23643..3286f139 100644 --- 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 `_. +- 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 index 3cdda29e..ead8b039 100644 --- 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 new file mode 100644 index 00000000..58a69280 --- /dev/null +++ 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 `_ and +`LAMMPS `_ as practical codes, and it +analizes and visualize the output using ``matplotlib`` and +`chemiscope `_. diff --git a/examples/thermostats/data/gle.lmp b/examples/thermostats/data/gle.lmp new file mode 100644 index 00000000..7a72c52b --- /dev/null +++ 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 new file mode 100644 index 00000000..2b49b207 --- /dev/null +++ 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 new file mode 100644 index 00000000..f6fcb9b3 --- /dev/null +++ b/examples/thermostats/data/input_cvv_sample.xml @@ -0,0 +1,45 @@ + + + + [ step, time{picosecond}, conserved{electronvolt}, + temperature{kelvin}, kinetic_md{electronvolt}, potential{electronvolt}, + temperature(H){kelvin}, temperature(O){kelvin} ] + + + 20000 + + 32342 + + +
h2o-lammps
1e-4 + + + + data/water_32.pdb + 300 + + + lmpserial + + + 300 + + + + 1.0 + + + [ 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 + ] + + + + + + diff --git a/examples/thermostats/data/input_cvv_traj.xml b/examples/thermostats/data/input_cvv_traj.xml new file mode 100644 index 00000000..d9e30463 --- /dev/null +++ b/examples/thermostats/data/input_cvv_traj.xml @@ -0,0 +1,48 @@ + + + + [ step, time{picosecond}, conserved{electronvolt}, + temperature{kelvin}, kinetic_md{electronvolt}, potential{electronvolt}, + temperature(H){kelvin}, temperature(O){kelvin} ] + positions + velocities + + 2000 + + 32342 + + +
h2o-lammps
1e-4 + + + ['ITRAJ'] + [1] + [2] + [3] + [4] + [5] + [6] + [7] + [8 ] + + + diff --git a/examples/thermostats/data/input_gle.xml b/examples/thermostats/data/input_gle.xml new file mode 100644 index 00000000..497d54ce --- /dev/null +++ b/examples/thermostats/data/input_gle.xml @@ -0,0 +1,46 @@ + + + + [ step, time{picosecond}, conserved{electronvolt}, + temperature{kelvin}, kinetic_md{electronvolt}, potential{electronvolt}, + temperature(H){kelvin}, temperature(O){kelvin} ] + positions + velocities + + 2000 + + 32342 + + +
h2o-lammps
1e-4 + + + + data/water_32.pdb + 300 + + + lmpserial + + + 300 + + + + 1.0 + + + [ 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 + ] + + + + + + diff --git a/examples/thermostats/data/input_higamma.xml b/examples/thermostats/data/input_higamma.xml new file mode 100644 index 00000000..190f089a --- /dev/null +++ b/examples/thermostats/data/input_higamma.xml @@ -0,0 +1,37 @@ + + + + [ step, time{picosecond}, conserved{electronvolt}, + temperature{kelvin}, kinetic_md{electronvolt}, potential{electronvolt}, + temperature(H){kelvin}, temperature(O){kelvin} ] + positions + velocities + + 2000 + + 32342 + + +
h2o-lammps
1e-4 + + + + data/water_32.pdb + 300 + + + lmpserial + + + 300 + + + + 1.0 + + 10 + + + + + diff --git a/examples/thermostats/data/input_nve.xml b/examples/thermostats/data/input_nve.xml new file mode 100644 index 00000000..2a667964 --- /dev/null +++ b/examples/thermostats/data/input_nve.xml @@ -0,0 +1,34 @@ + + + + [ step, time{picosecond}, conserved{electronvolt}, + temperature{kelvin}, kinetic_md{electronvolt}, potential{electronvolt}, + temperature(H){kelvin}, temperature(O){kelvin} ] + positions + velocities + + 2000 + + 32342 + + +
h2o-lammps
1e-4 + + + + data/water_32.pdb + 300 + + + lmpserial + + + 300 + + + + 1.0 + + + + diff --git a/examples/thermostats/data/input_svr.xml b/examples/thermostats/data/input_svr.xml new file mode 100644 index 00000000..e4457fb5 --- /dev/null +++ b/examples/thermostats/data/input_svr.xml @@ -0,0 +1,37 @@ + + + + [ step, time{picosecond}, conserved{electronvolt}, + temperature{kelvin}, kinetic_md{electronvolt}, potential{electronvolt}, + temperature(H){kelvin}, temperature(O){kelvin} ] + positions + velocities + + 2000 + + 32342 + + +
h2o-lammps
1e-4 + + + + data/water_32.pdb + 300 + + + lmpserial + + + 300 + + + + 1.0 + + 10 + + + + + diff --git a/examples/thermostats/data/run_traj.sh b/examples/thermostats/data/run_traj.sh new file mode 100644 index 00000000..445bca42 --- /dev/null +++ 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 new file mode 100644 index 00000000..60dac0ef --- /dev/null +++ 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 diff --git a/examples/thermostats/data/traj-all_facf.data b/examples/thermostats/data/traj-all_facf.data new file mode 100644 index 00000000..05f908b1 --- /dev/null +++ b/examples/thermostats/data/traj-all_facf.data @@ -0,0 +1,1000 @@ +0.0000e+00 2.3051e-05 3.7476e-06 +5.0661e-05 2.7021e-05 3.5291e-06 +1.0132e-04 3.6464e-05 3.4979e-06 +1.5198e-04 4.5866e-05 3.5446e-06 +2.0264e-04 5.0656e-05 3.2834e-06 +2.5331e-04 4.9747e-05 3.6485e-06 +3.0397e-04 4.4888e-05 4.5096e-06 +3.5463e-04 3.8662e-05 4.7062e-06 +4.0529e-04 3.3333e-05 4.1907e-06 +4.5595e-04 3.0269e-05 3.6770e-06 +5.0661e-04 2.9355e-05 3.6804e-06 +5.5727e-04 2.9358e-05 4.1701e-06 +6.0793e-04 2.9342e-05 4.6282e-06 +6.5859e-04 2.9408e-05 4.5458e-06 +7.0925e-04 3.0026e-05 3.9879e-06 +7.5992e-04 3.1230e-05 3.5481e-06 +8.1058e-04 3.2642e-05 3.4626e-06 +8.6124e-04 3.3972e-05 3.4608e-06 +9.1190e-04 3.5249e-05 3.4353e-06 +9.6256e-04 3.6614e-05 3.4057e-06 +1.0132e-03 3.7940e-05 3.3900e-06 +1.0639e-03 3.8790e-05 3.2214e-06 +1.1145e-03 3.8782e-05 2.6005e-06 +1.1652e-03 3.7859e-05 2.1349e-06 +1.2159e-03 3.6203e-05 2.5377e-06 +1.2665e-03 3.4201e-05 3.0619e-06 +1.3172e-03 3.2552e-05 3.3960e-06 +1.3679e-03 3.2017e-05 3.7150e-06 +1.4185e-03 3.2847e-05 4.4079e-06 +1.4692e-03 3.4674e-05 5.9052e-06 +1.5198e-03 3.7078e-05 7.6977e-06 +1.5705e-03 4.0038e-05 8.9979e-06 +1.6212e-03 4.3817e-05 9.6154e-06 +1.6718e-03 4.8592e-05 9.7467e-06 +1.7225e-03 5.4141e-05 9.3565e-06 +1.7731e-03 5.9836e-05 8.5929e-06 +1.8238e-03 6.5222e-05 8.4988e-06 +1.8745e-03 7.0658e-05 9.7542e-06 +1.9251e-03 7.6873e-05 1.1362e-05 +1.9758e-03 8.3859e-05 1.1724e-05 +2.0264e-03 9.0876e-05 1.0283e-05 +2.0771e-03 9.7616e-05 8.2095e-06 +2.1278e-03 1.0463e-04 7.6559e-06 +2.1784e-03 1.1231e-04 9.2875e-06 +2.2291e-03 1.1974e-04 1.1721e-05 +2.2797e-03 1.2465e-04 1.3154e-05 +2.3304e-03 1.2513e-04 1.1728e-05 +2.3811e-03 1.2180e-04 7.9989e-06 +2.4317e-03 1.1779e-04 5.9961e-06 +2.4824e-03 1.1580e-04 7.5271e-06 +2.5330e-03 1.1594e-04 8.8668e-06 +2.5837e-03 1.1658e-04 9.2704e-06 +2.6344e-03 1.1642e-04 8.7243e-06 +2.6850e-03 1.1514e-04 6.4799e-06 +2.7357e-03 1.1300e-04 5.4857e-06 +2.7863e-03 1.1052e-04 7.1882e-06 +2.8370e-03 1.0817e-04 7.1246e-06 +2.8877e-03 1.0601e-04 5.6219e-06 +2.9383e-03 1.0357e-04 5.9077e-06 +2.9890e-03 1.0012e-04 7.7540e-06 +3.0397e-03 9.5356e-05 8.4763e-06 +3.0903e-03 8.9881e-05 8.4957e-06 +3.1410e-03 8.4845e-05 9.1169e-06 +3.1916e-03 8.1077e-05 8.8058e-06 +3.2423e-03 7.8510e-05 7.4460e-06 +3.2930e-03 7.6499e-05 7.3026e-06 +3.3436e-03 7.4522e-05 8.3409e-06 +3.3943e-03 7.2433e-05 9.2325e-06 +3.4450e-03 7.0360e-05 9.3820e-06 +3.4956e-03 6.8528e-05 8.6551e-06 +3.5463e-03 6.6803e-05 8.0439e-06 +3.5969e-03 6.4489e-05 8.3172e-06 +3.6476e-03 6.1236e-05 8.4782e-06 +3.6983e-03 5.7785e-05 7.9908e-06 +3.7489e-03 5.4921e-05 7.4818e-06 +3.7996e-03 5.2333e-05 7.4519e-06 +3.8502e-03 4.9525e-05 7.3419e-06 +3.9009e-03 4.6937e-05 6.4210e-06 +3.9516e-03 4.5092e-05 5.1268e-06 +4.0022e-03 4.3586e-05 4.6526e-06 +4.0529e-03 4.1775e-05 5.0913e-06 +4.1035e-03 3.9480e-05 5.9640e-06 +4.1542e-03 3.6525e-05 6.5822e-06 +4.2049e-03 3.2697e-05 6.2151e-06 +4.2555e-03 2.8466e-05 5.2396e-06 +4.3062e-03 2.4768e-05 4.5093e-06 +4.3568e-03 2.1929e-05 3.8279e-06 +4.4075e-03 1.9566e-05 3.0954e-06 +4.4582e-03 1.7413e-05 2.5890e-06 +4.5088e-03 1.5584e-05 2.3539e-06 +4.5595e-03 1.4173e-05 2.2033e-06 +4.6102e-03 1.3070e-05 1.9763e-06 +4.6608e-03 1.2104e-05 1.6519e-06 +4.7115e-03 1.1166e-05 1.3025e-06 +4.7621e-03 1.0283e-05 1.0579e-06 +4.8128e-03 9.5635e-06 9.7946e-07 +4.8635e-03 9.0516e-06 9.0923e-07 +4.9141e-03 8.6664e-06 8.3999e-07 +4.9648e-03 8.3079e-06 8.1296e-07 +5.0154e-03 7.9440e-06 8.0726e-07 +5.0661e-03 7.6152e-06 8.2460e-07 +5.1168e-03 7.3959e-06 8.4667e-07 +5.1674e-03 7.3188e-06 8.2944e-07 +5.2181e-03 7.3214e-06 7.7089e-07 +5.2687e-03 7.2941e-06 7.3119e-07 +5.3194e-03 7.1827e-06 7.4021e-07 +5.3701e-03 7.0245e-06 7.6849e-07 +5.4207e-03 6.8803e-06 8.1023e-07 +5.4714e-03 6.7576e-06 8.4096e-07 +5.5220e-03 6.6225e-06 8.3245e-07 +5.5727e-03 6.4667e-06 7.9557e-07 +5.6234e-03 6.3178e-06 7.3789e-07 +5.6740e-03 6.1930e-06 6.5525e-07 +5.7247e-03 6.0783e-06 5.8981e-07 +5.7754e-03 5.9574e-06 5.8345e-07 +5.8260e-03 5.8315e-06 5.8274e-07 +5.8767e-03 5.7097e-06 5.4620e-07 +5.9273e-03 5.6042e-06 4.9871e-07 +5.9780e-03 5.5390e-06 4.9525e-07 +6.0287e-03 5.5308e-06 6.0693e-07 +6.0793e-03 5.5491e-06 7.3430e-07 +6.1300e-03 5.5299e-06 7.6070e-07 +6.1806e-03 5.4518e-06 6.9647e-07 +6.2313e-03 5.3679e-06 6.2021e-07 +6.2820e-03 5.3460e-06 5.9492e-07 +6.3326e-03 5.4130e-06 6.4836e-07 +6.3833e-03 5.5708e-06 7.5009e-07 +6.4339e-03 5.8090e-06 8.2377e-07 +6.4846e-03 6.0809e-06 8.7562e-07 +6.5353e-03 6.3337e-06 9.4956e-07 +6.5859e-03 6.5777e-06 9.8707e-07 +6.6366e-03 6.8506e-06 9.5639e-07 +6.6872e-03 7.1356e-06 9.6201e-07 +6.7379e-03 7.4136e-06 1.0285e-06 +6.7886e-03 7.7503e-06 1.0944e-06 +6.8392e-03 8.2472e-06 1.1743e-06 +6.8899e-03 8.9670e-06 1.2400e-06 +6.9406e-03 9.9545e-06 1.3938e-06 +6.9912e-03 1.1273e-05 1.6516e-06 +7.0419e-03 1.3036e-05 1.8867e-06 +7.0925e-03 1.5431e-05 1.9018e-06 +7.1432e-03 1.8717e-05 1.6586e-06 +7.1939e-03 2.3360e-05 1.5028e-06 +7.2445e-03 3.0325e-05 2.1126e-06 +7.2952e-03 4.1000e-05 4.2263e-06 +7.3458e-03 5.6563e-05 7.4587e-06 +7.3965e-03 7.7158e-05 9.8889e-06 +7.4472e-03 1.0095e-04 1.1086e-05 +7.4978e-03 1.2348e-04 1.4046e-05 +7.5485e-03 1.3884e-04 1.7532e-05 +7.5991e-03 1.4291e-04 1.8959e-05 +7.6498e-03 1.3563e-04 1.9262e-05 +7.7005e-03 1.1995e-04 1.9370e-05 +7.7511e-03 9.9352e-05 1.8063e-05 +7.8018e-03 7.6985e-05 1.4969e-05 +7.8525e-03 5.5917e-05 1.1553e-05 +7.9031e-03 3.8726e-05 8.6078e-06 +7.9538e-03 2.6460e-05 5.6072e-06 +8.0044e-03 1.8470e-05 2.9282e-06 +8.0551e-03 1.3360e-05 1.5099e-06 +8.1058e-03 9.9902e-06 1.0497e-06 +8.1564e-03 7.7273e-06 7.5944e-07 +8.2071e-03 6.2117e-06 5.2931e-07 +8.2577e-03 5.2064e-06 4.5573e-07 +8.3084e-03 4.5795e-06 5.6332e-07 +8.3591e-03 4.2233e-06 7.3450e-07 +8.4097e-03 3.9818e-06 8.4358e-07 +8.4604e-03 3.7232e-06 8.2997e-07 +8.5110e-03 3.4248e-06 6.9780e-07 +8.5617e-03 3.1325e-06 5.2913e-07 +8.6124e-03 2.8852e-06 4.3412e-07 +8.6630e-03 2.6964e-06 4.0455e-07 +8.7137e-03 2.5594e-06 3.7778e-07 +8.7644e-03 2.4513e-06 3.6420e-07 +8.8150e-03 2.3545e-06 3.6152e-07 +8.8657e-03 2.2719e-06 3.4211e-07 +8.9163e-03 2.2109e-06 3.1316e-07 +8.9670e-03 2.1609e-06 2.8267e-07 +9.0177e-03 2.1013e-06 2.3965e-07 +9.0683e-03 2.0268e-06 1.9001e-07 +9.1190e-03 1.9477e-06 1.7532e-07 +9.1696e-03 1.8695e-06 1.9812e-07 +9.2203e-03 1.7903e-06 2.1098e-07 +9.2710e-03 1.7169e-06 2.0625e-07 +9.3216e-03 1.6645e-06 1.9836e-07 +9.3723e-03 1.6386e-06 1.9405e-07 +9.4229e-03 1.6289e-06 1.9112e-07 +9.4736e-03 1.6185e-06 1.8639e-07 +9.5243e-03 1.5958e-06 1.8686e-07 +9.5749e-03 1.5615e-06 1.9313e-07 +9.6256e-03 1.5284e-06 1.8892e-07 +9.6763e-03 1.5067e-06 1.8347e-07 +9.7269e-03 1.4944e-06 1.9766e-07 +9.7776e-03 1.4822e-06 2.2244e-07 +9.8282e-03 1.4638e-06 2.3441e-07 +9.8789e-03 1.4383e-06 2.2592e-07 +9.9296e-03 1.4087e-06 2.0875e-07 +9.9802e-03 1.3798e-06 1.9847e-07 +1.0031e-02 1.3523e-06 1.9658e-07 +1.0082e-02 1.3192e-06 1.9005e-07 +1.0132e-02 1.2702e-06 1.7070e-07 +1.0183e-02 1.1998e-06 1.4491e-07 +1.0233e-02 1.1153e-06 1.2930e-07 +1.0284e-02 1.0351e-06 1.3154e-07 +1.0335e-02 9.7561e-07 1.3815e-07 +1.0386e-02 9.3625e-07 1.3821e-07 +1.0436e-02 9.0372e-07 1.2801e-07 +1.0487e-02 8.6788e-07 1.1320e-07 +1.0537e-02 8.2726e-07 1.1330e-07 +1.0588e-02 7.8295e-07 1.2945e-07 +1.0639e-02 7.3609e-07 1.4575e-07 +1.0689e-02 6.9044e-07 1.4869e-07 +1.0740e-02 6.5051e-07 1.3736e-07 +1.0791e-02 6.1797e-07 1.2163e-07 +1.0842e-02 5.9281e-07 1.0909e-07 +1.0892e-02 5.7468e-07 1.0113e-07 +1.0943e-02 5.6032e-07 9.5103e-08 +1.0993e-02 5.4446e-07 8.6445e-08 +1.1044e-02 5.2591e-07 7.5516e-08 +1.1095e-02 5.0781e-07 6.8237e-08 +1.1145e-02 4.9107e-07 6.7794e-08 +1.1196e-02 4.7332e-07 7.1711e-08 +1.1247e-02 4.5451e-07 7.4689e-08 +1.1297e-02 4.3785e-07 7.4606e-08 +1.1348e-02 4.2469e-07 7.3655e-08 +1.1399e-02 4.1337e-07 7.3453e-08 +1.1449e-02 4.0336e-07 7.2831e-08 +1.1500e-02 3.9600e-07 7.1766e-08 +1.1551e-02 3.9105e-07 7.2144e-08 +1.1601e-02 3.8682e-07 7.2354e-08 +1.1652e-02 3.8427e-07 7.0905e-08 +1.1703e-02 3.8625e-07 6.9668e-08 +1.1753e-02 3.9189e-07 7.0784e-08 +1.1804e-02 3.9628e-07 7.4614e-08 +1.1855e-02 3.9788e-07 7.9396e-08 +1.1905e-02 4.0170e-07 8.4152e-08 +1.1956e-02 4.1262e-07 9.1593e-08 +1.2007e-02 4.2861e-07 1.0197e-07 +1.2057e-02 4.4368e-07 1.1098e-07 +1.2108e-02 4.5605e-07 1.1733e-07 +1.2159e-02 4.7057e-07 1.2379e-07 +1.2209e-02 4.9232e-07 1.3305e-07 +1.2260e-02 5.2010e-07 1.4391e-07 +1.2311e-02 5.4763e-07 1.5186e-07 +1.2361e-02 5.7021e-07 1.5345e-07 +1.2412e-02 5.8804e-07 1.4961e-07 +1.2463e-02 6.0444e-07 1.4488e-07 +1.2513e-02 6.2312e-07 1.4179e-07 +1.2564e-02 6.4665e-07 1.3895e-07 +1.2615e-02 6.7543e-07 1.3349e-07 +1.2665e-02 7.0730e-07 1.2552e-07 +1.2716e-02 7.3908e-07 1.2229e-07 +1.2767e-02 7.6896e-07 1.3103e-07 +1.2817e-02 7.9759e-07 1.4773e-07 +1.2868e-02 8.2719e-07 1.6540e-07 +1.2919e-02 8.5966e-07 1.8090e-07 +1.2969e-02 8.9558e-07 1.9232e-07 +1.3020e-02 9.3468e-07 2.0261e-07 +1.3071e-02 9.7710e-07 2.1711e-07 +1.3121e-02 1.0220e-06 2.3501e-07 +1.3172e-02 1.0645e-06 2.5026e-07 +1.3222e-02 1.0987e-06 2.5906e-07 +1.3273e-02 1.1260e-06 2.6210e-07 +1.3324e-02 1.1540e-06 2.6019e-07 +1.3374e-02 1.1870e-06 2.5406e-07 +1.3425e-02 1.2268e-06 2.4917e-07 +1.3476e-02 1.2803e-06 2.5423e-07 +1.3527e-02 1.3515e-06 2.6958e-07 +1.3577e-02 1.4286e-06 2.8348e-07 +1.3628e-02 1.4954e-06 2.8929e-07 +1.3678e-02 1.5526e-06 2.9509e-07 +1.3729e-02 1.6142e-06 3.1098e-07 +1.3780e-02 1.6858e-06 3.3655e-07 +1.3831e-02 1.7579e-06 3.5693e-07 +1.3881e-02 1.8191e-06 3.5856e-07 +1.3932e-02 1.8697e-06 3.4747e-07 +1.3982e-02 1.9202e-06 3.3969e-07 +1.4033e-02 1.9754e-06 3.4454e-07 +1.4084e-02 2.0225e-06 3.5906e-07 +1.4134e-02 2.0484e-06 3.7591e-07 +1.4185e-02 2.0652e-06 3.9318e-07 +1.4236e-02 2.1027e-06 4.1386e-07 +1.4286e-02 2.1724e-06 4.4170e-07 +1.4337e-02 2.2589e-06 4.7456e-07 +1.4388e-02 2.3448e-06 5.0154e-07 +1.4438e-02 2.4261e-06 5.1393e-07 +1.4489e-02 2.5052e-06 5.1893e-07 +1.4540e-02 2.5910e-06 5.3104e-07 +1.4590e-02 2.7035e-06 5.4978e-07 +1.4641e-02 2.8594e-06 5.5959e-07 +1.4692e-02 3.0531e-06 5.5559e-07 +1.4742e-02 3.2639e-06 5.5228e-07 +1.4793e-02 3.4862e-06 5.5602e-07 +1.4844e-02 3.7471e-06 5.7708e-07 +1.4894e-02 4.0818e-06 6.5658e-07 +1.4945e-02 4.4966e-06 7.7501e-07 +1.4996e-02 4.9670e-06 8.4208e-07 +1.5046e-02 5.4623e-06 8.6037e-07 +1.5097e-02 5.9487e-06 9.0608e-07 +1.5148e-02 6.3955e-06 1.0041e-06 +1.5198e-02 6.7975e-06 1.1554e-06 +1.5249e-02 7.1528e-06 1.3926e-06 +1.5300e-02 7.4252e-06 1.6761e-06 +1.5350e-02 7.5993e-06 1.9006e-06 +1.5401e-02 7.7493e-06 2.0448e-06 +1.5452e-02 7.9822e-06 2.1739e-06 +1.5502e-02 8.3680e-06 2.3625e-06 +1.5553e-02 8.9751e-06 2.6904e-06 +1.5604e-02 9.8668e-06 3.2094e-06 +1.5654e-02 1.1019e-05 3.7917e-06 +1.5705e-02 1.2394e-05 4.2168e-06 +1.5756e-02 1.4185e-05 4.5039e-06 +1.5806e-02 1.6834e-05 5.0815e-06 +1.5857e-02 2.0697e-05 6.5664e-06 +1.5908e-02 2.5777e-05 9.1576e-06 +1.5958e-02 3.1896e-05 1.2184e-05 +1.6009e-02 3.8901e-05 1.4983e-05 +1.6060e-02 4.6521e-05 1.7545e-05 +1.6110e-02 5.4407e-05 1.9626e-05 +1.6161e-02 6.2629e-05 2.0706e-05 +1.6212e-02 7.1559e-05 2.1220e-05 +1.6262e-02 8.1057e-05 2.2043e-05 +1.6313e-02 9.0580e-05 2.3696e-05 +1.6363e-02 1.0015e-04 2.5867e-05 +1.6414e-02 1.1011e-04 2.7340e-05 +1.6465e-02 1.1998e-04 2.7077e-05 +1.6515e-02 1.2853e-04 2.4728e-05 +1.6566e-02 1.3458e-04 2.1439e-05 +1.6617e-02 1.3708e-04 1.8761e-05 +1.6667e-02 1.3570e-04 1.7668e-05 +1.6718e-02 1.3166e-04 1.7744e-05 +1.6769e-02 1.2669e-04 1.7296e-05 +1.6820e-02 1.2097e-04 1.5430e-05 +1.6870e-02 1.1351e-04 1.2460e-05 +1.6921e-02 1.0429e-04 9.8761e-06 +1.6971e-02 9.4406e-05 9.0052e-06 +1.7022e-02 8.4584e-05 9.2360e-06 +1.7073e-02 7.4995e-05 8.4495e-06 +1.7123e-02 6.5998e-05 6.3262e-06 +1.7174e-02 5.8160e-05 5.0162e-06 +1.7225e-02 5.1820e-05 5.0566e-06 +1.7275e-02 4.7058e-05 5.0931e-06 +1.7326e-02 4.3708e-05 4.9178e-06 +1.7377e-02 4.1248e-05 4.9520e-06 +1.7427e-02 3.9070e-05 5.2966e-06 +1.7478e-02 3.6826e-05 5.6220e-06 +1.7529e-02 3.4257e-05 5.7692e-06 +1.7579e-02 3.1001e-05 5.4544e-06 +1.7630e-02 2.6978e-05 4.6567e-06 +1.7681e-02 2.2697e-05 3.9411e-06 +1.7731e-02 1.8784e-05 3.6316e-06 +1.7782e-02 1.5414e-05 3.3379e-06 +1.7833e-02 1.2444e-05 2.8140e-06 +1.7883e-02 9.8118e-06 2.2141e-06 +1.7934e-02 7.5822e-06 1.6983e-06 +1.7985e-02 5.8035e-06 1.3211e-06 +1.8035e-02 4.4748e-06 1.0937e-06 +1.8086e-02 3.5453e-06 9.4220e-07 +1.8137e-02 2.9108e-06 7.8540e-07 +1.8187e-02 2.4690e-06 6.3001e-07 +1.8238e-02 2.1659e-06 5.0777e-07 +1.8289e-02 1.9699e-06 4.2328e-07 +1.8339e-02 1.8475e-06 3.7263e-07 +1.8390e-02 1.7736e-06 3.5126e-07 +1.8441e-02 1.7414e-06 3.5407e-07 +1.8491e-02 1.7508e-06 3.7166e-07 +1.8542e-02 1.7913e-06 3.8716e-07 +1.8593e-02 1.8410e-06 3.9971e-07 +1.8643e-02 1.8868e-06 4.1137e-07 +1.8694e-02 1.9292e-06 4.1079e-07 +1.8745e-02 1.9647e-06 3.9854e-07 +1.8795e-02 1.9833e-06 3.7823e-07 +1.8846e-02 1.9874e-06 3.4910e-07 +1.8897e-02 1.9943e-06 3.1990e-07 +1.8947e-02 2.0139e-06 3.0364e-07 +1.8998e-02 2.0389e-06 3.0310e-07 +1.9048e-02 2.0551e-06 3.0779e-07 +1.9099e-02 2.0499e-06 2.9920e-07 +1.9150e-02 2.0134e-06 2.6714e-07 +1.9200e-02 1.9462e-06 2.2775e-07 +1.9251e-02 1.8642e-06 2.0805e-07 +1.9302e-02 1.7831e-06 2.1575e-07 +1.9353e-02 1.7040e-06 2.2642e-07 +1.9403e-02 1.6198e-06 2.2150e-07 +1.9454e-02 1.5311e-06 2.1084e-07 +1.9505e-02 1.4460e-06 2.0208e-07 +1.9555e-02 1.3658e-06 1.8621e-07 +1.9606e-02 1.2807e-06 1.6548e-07 +1.9657e-02 1.1824e-06 1.4784e-07 +1.9707e-02 1.0782e-06 1.2697e-07 +1.9758e-02 9.8424e-07 1.1055e-07 +1.9808e-02 9.0523e-07 1.0438e-07 +1.9859e-02 8.3279e-07 9.6613e-08 +1.9910e-02 7.6018e-07 8.6772e-08 +1.9960e-02 6.8957e-07 8.0287e-08 +2.0011e-02 6.2518e-07 8.0050e-08 +2.0062e-02 5.6789e-07 8.0660e-08 +2.0112e-02 5.1764e-07 7.5962e-08 +2.0163e-02 4.7497e-07 6.5444e-08 +2.0214e-02 4.3848e-07 5.2453e-08 +2.0264e-02 4.0456e-07 4.2760e-08 +2.0315e-02 3.7176e-07 3.9818e-08 +2.0366e-02 3.4222e-07 4.0146e-08 +2.0416e-02 3.1734e-07 4.1453e-08 +2.0467e-02 2.9482e-07 4.1012e-08 +2.0518e-02 2.7140e-07 3.7152e-08 +2.0568e-02 2.4664e-07 3.2705e-08 +2.0619e-02 2.2282e-07 2.9760e-08 +2.0670e-02 2.0211e-07 2.7865e-08 +2.0720e-02 1.8466e-07 2.5105e-08 +2.0771e-02 1.6933e-07 2.0789e-08 +2.0822e-02 1.5521e-07 1.7210e-08 +2.0872e-02 1.4206e-07 1.6278e-08 +2.0923e-02 1.2998e-07 1.6842e-08 +2.0974e-02 1.1916e-07 1.7255e-08 +2.1024e-02 1.0993e-07 1.7612e-08 +2.1075e-02 1.0235e-07 1.8047e-08 +2.1126e-02 9.5888e-08 1.7985e-08 +2.1176e-02 8.9948e-08 1.6861e-08 +2.1227e-02 8.4568e-08 1.4970e-08 +2.1278e-02 7.9989e-08 1.3090e-08 +2.1328e-02 7.5947e-08 1.1557e-08 +2.1379e-02 7.2057e-08 1.0491e-08 +2.1430e-02 6.8391e-08 1.0033e-08 +2.1480e-02 6.5181e-08 9.9310e-09 +2.1531e-02 6.2390e-08 9.7864e-09 +2.1582e-02 5.9827e-08 9.4018e-09 +2.1632e-02 5.7431e-08 8.8306e-09 +2.1683e-02 5.5291e-08 8.2540e-09 +2.1734e-02 5.3438e-08 7.7638e-09 +2.1784e-02 5.1798e-08 7.3067e-09 +2.1835e-02 5.0383e-08 6.8381e-09 +2.1885e-02 4.9251e-08 6.3802e-09 +2.1936e-02 4.8256e-08 6.0528e-09 +2.1987e-02 4.7145e-08 5.9063e-09 +2.2038e-02 4.5872e-08 5.9082e-09 +2.2088e-02 4.4684e-08 6.0791e-09 +2.2139e-02 4.3869e-08 6.3509e-09 +2.2190e-02 4.3386e-08 6.4774e-09 +2.2240e-02 4.2840e-08 6.2821e-09 +2.2291e-02 4.1953e-08 5.9744e-09 +2.2342e-02 4.0856e-08 5.7523e-09 +2.2392e-02 3.9842e-08 5.4469e-09 +2.2443e-02 3.9085e-08 5.0654e-09 +2.2493e-02 3.8605e-08 4.9690e-09 +2.2544e-02 3.8335e-08 5.2505e-09 +2.2595e-02 3.8167e-08 5.6028e-09 +2.2645e-02 3.7929e-08 5.7737e-09 +2.2696e-02 3.7427e-08 5.6979e-09 +2.2747e-02 3.6681e-08 5.3528e-09 +2.2797e-02 3.5967e-08 4.8847e-09 +2.2848e-02 3.5496e-08 4.6080e-09 +2.2899e-02 3.5195e-08 4.6382e-09 +2.2949e-02 3.4862e-08 4.8154e-09 +2.3000e-02 3.4435e-08 4.8857e-09 +2.3051e-02 3.4032e-08 4.8031e-09 +2.3101e-02 3.3760e-08 4.8221e-09 +2.3152e-02 3.3608e-08 5.0989e-09 +2.3203e-02 3.3555e-08 5.5295e-09 +2.3253e-02 3.3569e-08 5.9201e-09 +2.3304e-02 3.3532e-08 6.1854e-09 +2.3355e-02 3.3417e-08 6.3115e-09 +2.3405e-02 3.3395e-08 6.2923e-09 +2.3456e-02 3.3573e-08 6.2065e-09 +2.3507e-02 3.3807e-08 6.1747e-09 +2.3557e-02 3.3954e-08 6.2466e-09 +2.3608e-02 3.4163e-08 6.3762e-09 +2.3659e-02 3.4755e-08 6.5420e-09 +2.3709e-02 3.5822e-08 6.7884e-09 +2.3760e-02 3.7129e-08 7.0870e-09 +2.3811e-02 3.8477e-08 7.3475e-09 +2.3861e-02 3.9936e-08 7.5887e-09 +2.3912e-02 4.1590e-08 7.8522e-09 +2.3963e-02 4.3281e-08 8.0602e-09 +2.4013e-02 4.4750e-08 8.1439e-09 +2.4064e-02 4.5844e-08 8.1190e-09 +2.4115e-02 4.6584e-08 8.0945e-09 +2.4165e-02 4.7167e-08 8.2049e-09 +2.4216e-02 4.7806e-08 8.2170e-09 +2.4267e-02 4.8437e-08 7.7108e-09 +2.4317e-02 4.8706e-08 6.9444e-09 +2.4368e-02 4.8391e-08 6.6576e-09 +2.4419e-02 4.7587e-08 6.7933e-09 +2.4469e-02 4.6391e-08 6.7367e-09 +2.4520e-02 4.4776e-08 6.2176e-09 +2.4571e-02 4.2814e-08 5.1608e-09 +2.4621e-02 4.0667e-08 3.8300e-09 +2.4672e-02 3.8382e-08 2.9417e-09 +2.4723e-02 3.6001e-08 2.6792e-09 +2.4773e-02 3.3697e-08 2.6990e-09 +2.4824e-02 3.1618e-08 2.7601e-09 +2.4875e-02 2.9756e-08 2.8196e-09 +2.4925e-02 2.8066e-08 2.9160e-09 +2.4976e-02 2.6566e-08 2.8810e-09 +2.5027e-02 2.5261e-08 2.6302e-09 +2.5077e-02 2.4048e-08 2.4308e-09 +2.5128e-02 2.2817e-08 2.5035e-09 +2.5178e-02 2.1612e-08 2.6644e-09 +2.5229e-02 2.0563e-08 2.6770e-09 +2.5280e-02 1.9673e-08 2.5395e-09 +2.5330e-02 1.8840e-08 2.3964e-09 +2.5381e-02 1.8050e-08 2.2902e-09 +2.5432e-02 1.7382e-08 2.1684e-09 +2.5483e-02 1.6885e-08 2.0224e-09 +2.5533e-02 1.6502e-08 1.9088e-09 +2.5584e-02 1.6123e-08 1.8437e-09 +2.5635e-02 1.5715e-08 1.7751e-09 +2.5685e-02 1.5365e-08 1.6957e-09 +2.5736e-02 1.5154e-08 1.6474e-09 +2.5786e-02 1.5055e-08 1.6518e-09 +2.5837e-02 1.4987e-08 1.6958e-09 +2.5888e-02 1.4906e-08 1.7207e-09 +2.5938e-02 1.4801e-08 1.6687e-09 +2.5989e-02 1.4665e-08 1.5457e-09 +2.6040e-02 1.4502e-08 1.4505e-09 +2.6090e-02 1.4354e-08 1.4739e-09 +2.6141e-02 1.4268e-08 1.5823e-09 +2.6192e-02 1.4232e-08 1.6745e-09 +2.6242e-02 1.4190e-08 1.6824e-09 +2.6293e-02 1.4121e-08 1.6168e-09 +2.6344e-02 1.4037e-08 1.5302e-09 +2.6394e-02 1.3944e-08 1.4640e-09 +2.6445e-02 1.3820e-08 1.4151e-09 +2.6496e-02 1.3649e-08 1.3550e-09 +2.6546e-02 1.3478e-08 1.2958e-09 +2.6597e-02 1.3385e-08 1.2681e-09 +2.6648e-02 1.3391e-08 1.2932e-09 +2.6698e-02 1.3446e-08 1.3863e-09 +2.6749e-02 1.3490e-08 1.5013e-09 +2.6800e-02 1.3483e-08 1.5549e-09 +2.6850e-02 1.3396e-08 1.5403e-09 +2.6901e-02 1.3216e-08 1.5068e-09 +2.6952e-02 1.2975e-08 1.4653e-09 +2.7002e-02 1.2749e-08 1.4235e-09 +2.7053e-02 1.2593e-08 1.4030e-09 +2.7104e-02 1.2487e-08 1.3951e-09 +2.7154e-02 1.2385e-08 1.3501e-09 +2.7205e-02 1.2279e-08 1.2426e-09 +2.7256e-02 1.2205e-08 1.1284e-09 +2.7306e-02 1.2182e-08 1.0886e-09 +2.7357e-02 1.2163e-08 1.1511e-09 +2.7408e-02 1.2087e-08 1.2618e-09 +2.7458e-02 1.1969e-08 1.3401e-09 +2.7509e-02 1.1875e-08 1.3621e-09 +2.7560e-02 1.1810e-08 1.3603e-09 +2.7610e-02 1.1713e-08 1.3606e-09 +2.7661e-02 1.1564e-08 1.3632e-09 +2.7712e-02 1.1411e-08 1.3604e-09 +2.7762e-02 1.1291e-08 1.3565e-09 +2.7813e-02 1.1203e-08 1.3580e-09 +2.7863e-02 1.1137e-08 1.3487e-09 +2.7914e-02 1.1098e-08 1.3083e-09 +2.7965e-02 1.1075e-08 1.2536e-09 +2.8015e-02 1.1036e-08 1.2156e-09 +2.8066e-02 1.0974e-08 1.2018e-09 +2.8117e-02 1.0910e-08 1.1999e-09 +2.8168e-02 1.0838e-08 1.1916e-09 +2.8218e-02 1.0742e-08 1.1909e-09 +2.8269e-02 1.0642e-08 1.2213e-09 +2.8320e-02 1.0580e-08 1.2597e-09 +2.8370e-02 1.0569e-08 1.2735e-09 +2.8421e-02 1.0584e-08 1.2777e-09 +2.8472e-02 1.0595e-08 1.3097e-09 +2.8522e-02 1.0590e-08 1.3614e-09 +2.8573e-02 1.0568e-08 1.3875e-09 +2.8623e-02 1.0534e-08 1.3753e-09 +2.8674e-02 1.0507e-08 1.3494e-09 +2.8725e-02 1.0503e-08 1.3285e-09 +2.8775e-02 1.0521e-08 1.3210e-09 +2.8826e-02 1.0556e-08 1.3402e-09 +2.8877e-02 1.0609e-08 1.3942e-09 +2.8927e-02 1.0659e-08 1.4717e-09 +2.8978e-02 1.0679e-08 1.5432e-09 +2.9029e-02 1.0666e-08 1.5838e-09 +2.9079e-02 1.0643e-08 1.5939e-09 +2.9130e-02 1.0652e-08 1.5997e-09 +2.9181e-02 1.0722e-08 1.6405e-09 +2.9231e-02 1.0827e-08 1.7248e-09 +2.9282e-02 1.0904e-08 1.7892e-09 +2.9333e-02 1.0931e-08 1.8056e-09 +2.9383e-02 1.0970e-08 1.8470e-09 +2.9434e-02 1.1084e-08 1.9530e-09 +2.9485e-02 1.1253e-08 2.0423e-09 +2.9535e-02 1.1410e-08 2.0525e-09 +2.9586e-02 1.1521e-08 2.0396e-09 +2.9637e-02 1.1592e-08 2.0526e-09 +2.9687e-02 1.1662e-08 2.1234e-09 +2.9738e-02 1.1808e-08 2.2723e-09 +2.9789e-02 1.2067e-08 2.4785e-09 +2.9839e-02 1.2368e-08 2.7155e-09 +2.9890e-02 1.2616e-08 2.9646e-09 +2.9941e-02 1.2809e-08 3.2032e-09 +2.9991e-02 1.3004e-08 3.4097e-09 +3.0042e-02 1.3246e-08 3.6142e-09 +3.0093e-02 1.3571e-08 3.8779e-09 +3.0143e-02 1.4018e-08 4.2213e-09 +3.0194e-02 1.4575e-08 4.6233e-09 +3.0245e-02 1.5148e-08 5.0563e-09 +3.0295e-02 1.5624e-08 5.4654e-09 +3.0346e-02 1.6033e-08 5.7990e-09 +3.0397e-02 1.6567e-08 6.1141e-09 +3.0447e-02 1.7351e-08 6.4955e-09 +3.0498e-02 1.8253e-08 6.8670e-09 +3.0549e-02 1.9039e-08 7.0632e-09 +3.0599e-02 1.9666e-08 7.0902e-09 +3.0650e-02 2.0291e-08 7.1530e-09 +3.0701e-02 2.0981e-08 7.3855e-09 +3.0751e-02 2.1608e-08 7.7325e-09 +3.0802e-02 2.2070e-08 8.1293e-09 +3.0853e-02 2.2457e-08 8.5985e-09 +3.0903e-02 2.2909e-08 9.1390e-09 +3.0954e-02 2.3499e-08 9.7208e-09 +3.1005e-02 2.4307e-08 1.0404e-08 +3.1055e-02 2.5428e-08 1.1263e-08 +3.1106e-02 2.6847e-08 1.2311e-08 +3.1157e-02 2.8477e-08 1.3663e-08 +3.1207e-02 3.0334e-08 1.5403e-08 +3.1258e-02 3.2519e-08 1.7278e-08 +3.1309e-02 3.5044e-08 1.8844e-08 +3.1359e-02 3.7882e-08 1.9959e-08 +3.1410e-02 4.1177e-08 2.1120e-08 +3.1461e-02 4.5255e-08 2.3164e-08 +3.1511e-02 5.0174e-08 2.6450e-08 +3.1562e-02 5.5279e-08 3.0156e-08 +3.1613e-02 5.9636e-08 3.2642e-08 +3.1663e-02 6.3214e-08 3.3041e-08 +3.1714e-02 6.7097e-08 3.2408e-08 +3.1764e-02 7.2109e-08 3.2782e-08 +3.1815e-02 7.7850e-08 3.4651e-08 +3.1866e-02 8.3419e-08 3.6282e-08 +3.1916e-02 8.8294e-08 3.6610e-08 +3.1967e-02 9.2397e-08 3.6915e-08 +3.2018e-02 9.6265e-08 3.8973e-08 +3.2068e-02 1.0113e-07 4.3698e-08 +3.2119e-02 1.0778e-07 5.0636e-08 +3.2170e-02 1.1531e-07 5.7427e-08 +3.2220e-02 1.2178e-07 6.0859e-08 +3.2271e-02 1.2592e-07 5.9757e-08 +3.2322e-02 1.2813e-07 5.6538e-08 +3.2372e-02 1.3030e-07 5.4506e-08 +3.2423e-02 1.3455e-07 5.4905e-08 +3.2474e-02 1.4103e-07 5.6741e-08 +3.2524e-02 1.4759e-07 5.8315e-08 +3.2575e-02 1.5266e-07 5.8888e-08 +3.2626e-02 1.5754e-07 5.8783e-08 +3.2676e-02 1.6387e-07 5.8165e-08 +3.2727e-02 1.7036e-07 5.6697e-08 +3.2778e-02 1.7478e-07 5.4385e-08 +3.2828e-02 1.7803e-07 5.2074e-08 +3.2879e-02 1.8297e-07 5.1604e-08 +3.2930e-02 1.8943e-07 5.3888e-08 +3.2980e-02 1.9433e-07 5.6293e-08 +3.3031e-02 1.9680e-07 5.7017e-08 +3.3082e-02 1.9925e-07 5.7604e-08 +3.3132e-02 2.0244e-07 5.9612e-08 +3.3183e-02 2.0318e-07 6.1567e-08 +3.3234e-02 1.9887e-07 6.0003e-08 +3.3284e-02 1.9123e-07 5.4324e-08 +3.3335e-02 1.8369e-07 4.7145e-08 +3.3386e-02 1.7737e-07 4.0850e-08 +3.3436e-02 1.7176e-07 3.6445e-08 +3.3487e-02 1.6698e-07 3.5104e-08 +3.3538e-02 1.6352e-07 3.6770e-08 +3.3588e-02 1.6046e-07 3.8698e-08 +3.3639e-02 1.5551e-07 3.9570e-08 +3.3690e-02 1.4733e-07 3.8876e-08 +3.3740e-02 1.3759e-07 3.6510e-08 +3.3791e-02 1.2988e-07 3.4728e-08 +3.3841e-02 1.2621e-07 3.5351e-08 +3.3892e-02 1.2539e-07 3.6536e-08 +3.3943e-02 1.2498e-07 3.5034e-08 +3.3994e-02 1.2316e-07 3.0955e-08 +3.4044e-02 1.1868e-07 2.6568e-08 +3.4095e-02 1.1153e-07 2.3409e-08 +3.4146e-02 1.0397e-07 2.3072e-08 +3.4196e-02 9.8728e-08 2.4815e-08 +3.4247e-02 9.6056e-08 2.5266e-08 +3.4298e-02 9.4009e-08 2.2516e-08 +3.4348e-02 9.1207e-08 1.8512e-08 +3.4399e-02 8.7633e-08 1.6091e-08 +3.4450e-02 8.3577e-08 1.5015e-08 +3.4500e-02 7.9353e-08 1.4506e-08 +3.4551e-02 7.5627e-08 1.4657e-08 +3.4602e-02 7.3021e-08 1.5848e-08 +3.4652e-02 7.1411e-08 1.9018e-08 +3.4703e-02 6.9778e-08 2.2192e-08 +3.4753e-02 6.6890e-08 2.1783e-08 +3.4804e-02 6.2462e-08 1.7926e-08 +3.4855e-02 5.7556e-08 1.4109e-08 +3.4905e-02 5.3413e-08 1.2707e-08 +3.4956e-02 5.0097e-08 1.2907e-08 +3.5007e-02 4.6709e-08 1.2927e-08 +3.5057e-02 4.2616e-08 1.1797e-08 +3.5108e-02 3.8085e-08 9.8300e-09 +3.5159e-02 3.3934e-08 8.1475e-09 +3.5209e-02 3.0846e-08 7.6148e-09 +3.5260e-02 2.8901e-08 7.9585e-09 +3.5311e-02 2.7598e-08 8.2572e-09 +3.5361e-02 2.6280e-08 8.0131e-09 +3.5412e-02 2.4554e-08 7.1419e-09 +3.5463e-02 2.2453e-08 5.8743e-09 +3.5513e-02 2.0319e-08 4.6199e-09 +3.5564e-02 1.8505e-08 3.7725e-09 +3.5615e-02 1.7135e-08 3.4137e-09 +3.5665e-02 1.6118e-08 3.2226e-09 +3.5716e-02 1.5277e-08 2.9606e-09 +3.5767e-02 1.4448e-08 2.6663e-09 +3.5817e-02 1.3598e-08 2.4536e-09 +3.5868e-02 1.2863e-08 2.3375e-09 +3.5919e-02 1.2395e-08 2.3058e-09 +3.5969e-02 1.2167e-08 2.3460e-09 +3.6020e-02 1.2001e-08 2.3696e-09 +3.6071e-02 1.1756e-08 2.2751e-09 +3.6121e-02 1.1437e-08 2.0828e-09 +3.6172e-02 1.1130e-08 1.8991e-09 +3.6223e-02 1.0884e-08 1.7292e-09 +3.6273e-02 1.0699e-08 1.5323e-09 +3.6324e-02 1.0568e-08 1.3986e-09 +3.6375e-02 1.0461e-08 1.4247e-09 +3.6425e-02 1.0321e-08 1.5306e-09 +3.6476e-02 1.0108e-08 1.5780e-09 +3.6527e-02 9.8191e-09 1.5050e-09 +3.6577e-02 9.4553e-09 1.3281e-09 +3.6628e-02 9.0276e-09 1.1409e-09 +3.6679e-02 8.5822e-09 1.0283e-09 +3.6729e-02 8.1773e-09 9.7364e-10 +3.6780e-02 7.8309e-09 9.3100e-10 +3.6831e-02 7.5164e-09 8.6989e-10 +3.6881e-02 7.2070e-09 7.7432e-10 +3.6932e-02 6.9076e-09 6.5802e-10 +3.6983e-02 6.6422e-09 5.8036e-10 +3.7033e-02 6.4198e-09 5.7880e-10 +3.7084e-02 6.2213e-09 6.1259e-10 +3.7135e-02 6.0230e-09 6.2850e-10 +3.7185e-02 5.8254e-09 6.1034e-10 +3.7236e-02 5.6552e-09 5.7611e-10 +3.7287e-02 5.5389e-09 5.3938e-10 +3.7337e-02 5.4756e-09 5.0254e-10 +3.7388e-02 5.4387e-09 4.7824e-10 +3.7438e-02 5.4026e-09 4.7746e-10 +3.7489e-02 5.3555e-09 4.8703e-10 +3.7540e-02 5.2871e-09 4.8809e-10 +3.7590e-02 5.1890e-09 4.8087e-10 +3.7641e-02 5.0759e-09 4.6980e-10 +3.7692e-02 4.9794e-09 4.5171e-10 +3.7742e-02 4.9157e-09 4.2592e-10 +3.7793e-02 4.8771e-09 3.9797e-10 +3.7844e-02 4.8490e-09 3.8108e-10 +3.7894e-02 4.8204e-09 3.8079e-10 +3.7945e-02 4.7861e-09 3.8235e-10 +3.7996e-02 4.7436e-09 3.6920e-10 +3.8046e-02 4.6893e-09 3.4389e-10 +3.8097e-02 4.6239e-09 3.2228e-10 +3.8148e-02 4.5612e-09 3.1467e-10 +3.8198e-02 4.5191e-09 3.2196e-10 +3.8249e-02 4.5010e-09 3.3779e-10 +3.8300e-02 4.4944e-09 3.5023e-10 +3.8350e-02 4.4869e-09 3.5107e-10 +3.8401e-02 4.4740e-09 3.3994e-10 +3.8452e-02 4.4534e-09 3.2375e-10 +3.8502e-02 4.4266e-09 3.1344e-10 +3.8553e-02 4.4030e-09 3.1523e-10 +3.8604e-02 4.3910e-09 3.2648e-10 +3.8654e-02 4.3858e-09 3.3869e-10 +3.8705e-02 4.3770e-09 3.4430e-10 +3.8756e-02 4.3630e-09 3.4292e-10 +3.8806e-02 4.3465e-09 3.3795e-10 +3.8857e-02 4.3254e-09 3.3088e-10 +3.8908e-02 4.2979e-09 3.2174e-10 +3.8958e-02 4.2700e-09 3.1043e-10 +3.9009e-02 4.2520e-09 2.9915e-10 +3.9060e-02 4.2457e-09 2.9187e-10 +3.9110e-02 4.2439e-09 2.8983e-10 +3.9161e-02 4.2425e-09 2.9231e-10 +3.9212e-02 4.2428e-09 2.9535e-10 +3.9262e-02 4.2398e-09 2.9325e-10 +3.9313e-02 4.2221e-09 2.8760e-10 +3.9364e-02 4.1881e-09 2.8491e-10 +3.9414e-02 4.1521e-09 2.8658e-10 +3.9465e-02 4.1300e-09 2.8755e-10 +3.9516e-02 4.1242e-09 2.8801e-10 +3.9566e-02 4.1258e-09 2.9667e-10 +3.9617e-02 4.1264e-09 3.1473e-10 +3.9668e-02 4.1195e-09 3.2915e-10 +3.9718e-02 4.0990e-09 3.2592e-10 +3.9769e-02 4.0678e-09 3.0852e-10 +3.9820e-02 4.0395e-09 2.9416e-10 +3.9870e-02 4.0248e-09 2.9174e-10 +3.9921e-02 4.0219e-09 2.9540e-10 +3.9972e-02 4.0237e-09 2.9843e-10 +4.0022e-02 4.0240e-09 3.0006e-10 +4.0073e-02 4.0200e-09 2.9930e-10 +4.0123e-02 4.0108e-09 2.9446e-10 +4.0174e-02 3.9989e-09 2.8829e-10 +4.0225e-02 3.9898e-09 2.8680e-10 +4.0275e-02 3.9862e-09 2.9144e-10 +4.0326e-02 3.9833e-09 2.9948e-10 +4.0377e-02 3.9758e-09 3.0943e-10 +4.0427e-02 3.9637e-09 3.2122e-10 +4.0478e-02 3.9494e-09 3.2994e-10 +4.0529e-02 3.9352e-09 3.2792e-10 +4.0579e-02 3.9260e-09 3.1668e-10 +4.0630e-02 3.9238e-09 3.0675e-10 +4.0681e-02 3.9233e-09 3.0570e-10 +4.0731e-02 3.9171e-09 3.1031e-10 +4.0782e-02 3.9020e-09 3.1039e-10 +4.0833e-02 3.8799e-09 3.0159e-10 +4.0883e-02 3.8556e-09 2.9029e-10 +4.0934e-02 3.8344e-09 2.8522e-10 +4.0985e-02 3.8200e-09 2.8711e-10 +4.1035e-02 3.8110e-09 2.8926e-10 +4.1086e-02 3.8026e-09 2.8616e-10 +4.1137e-02 3.7931e-09 2.7976e-10 +4.1187e-02 3.7844e-09 2.7404e-10 +4.1238e-02 3.7762e-09 2.6756e-10 +4.1289e-02 3.7646e-09 2.6021e-10 +4.1339e-02 3.7475e-09 2.5727e-10 +4.1390e-02 3.7277e-09 2.6017e-10 +4.1441e-02 3.7112e-09 2.6539e-10 +4.1491e-02 3.7026e-09 2.6979e-10 +4.1542e-02 3.7006e-09 2.7206e-10 +4.1593e-02 3.6981e-09 2.7206e-10 +4.1643e-02 3.6903e-09 2.6815e-10 +4.1694e-02 3.6796e-09 2.5949e-10 +4.1745e-02 3.6712e-09 2.5006e-10 +4.1795e-02 3.6640e-09 2.4491e-10 +4.1846e-02 3.6513e-09 2.4524e-10 +4.1897e-02 3.6319e-09 2.4916e-10 +4.1947e-02 3.6130e-09 2.5323e-10 +4.1998e-02 3.6002e-09 2.5480e-10 +4.2049e-02 3.5933e-09 2.5514e-10 +4.2099e-02 3.5891e-09 2.5622e-10 +4.2150e-02 3.5838e-09 2.5775e-10 +4.2201e-02 3.5743e-09 2.5895e-10 +4.2251e-02 3.5605e-09 2.5906e-10 +4.2302e-02 3.5457e-09 2.5751e-10 +4.2353e-02 3.5336e-09 2.5427e-10 +4.2403e-02 3.5245e-09 2.4882e-10 +4.2454e-02 3.5163e-09 2.4106e-10 +4.2505e-02 3.5076e-09 2.3427e-10 +4.2555e-02 3.4973e-09 2.3341e-10 +4.2606e-02 3.4838e-09 2.3817e-10 +4.2657e-02 3.4683e-09 2.4234e-10 +4.2707e-02 3.4555e-09 2.4202e-10 +4.2758e-02 3.4476e-09 2.3813e-10 +4.2808e-02 3.4431e-09 2.3376e-10 +4.2859e-02 3.4403e-09 2.3168e-10 +4.2910e-02 3.4365e-09 2.3227e-10 +4.2960e-02 3.4280e-09 2.3381e-10 +4.3011e-02 3.4148e-09 2.3473e-10 +4.3062e-02 3.4028e-09 2.3477e-10 +4.3112e-02 3.3980e-09 2.3430e-10 +4.3163e-02 3.3985e-09 2.3328e-10 +4.3214e-02 3.3965e-09 2.3175e-10 +4.3264e-02 3.3883e-09 2.3036e-10 +4.3315e-02 3.3774e-09 2.2946e-10 +4.3366e-02 3.3677e-09 2.2854e-10 +4.3416e-02 3.3603e-09 2.2772e-10 +4.3467e-02 3.3545e-09 2.2775e-10 +4.3518e-02 3.3485e-09 2.2863e-10 +4.3569e-02 3.3405e-09 2.3072e-10 +4.3619e-02 3.3310e-09 2.3394e-10 +4.3670e-02 3.3219e-09 2.3589e-10 +4.3720e-02 3.3143e-09 2.3442e-10 +4.3771e-02 3.3082e-09 2.3081e-10 +4.3822e-02 3.3020e-09 2.2814e-10 +4.3872e-02 3.2948e-09 2.2716e-10 +4.3923e-02 3.2868e-09 2.2599e-10 +4.3974e-02 3.2784e-09 2.2374e-10 +4.4024e-02 3.2704e-09 2.2195e-10 +4.4075e-02 3.2642e-09 2.2163e-10 +4.4126e-02 3.2593e-09 2.2137e-10 +4.4176e-02 3.2536e-09 2.1967e-10 +4.4227e-02 3.2465e-09 2.1721e-10 +4.4278e-02 3.2384e-09 2.1553e-10 +4.4328e-02 3.2296e-09 2.1481e-10 +4.4379e-02 3.2211e-09 2.1468e-10 +4.4430e-02 3.2146e-09 2.1586e-10 +4.4480e-02 3.2101e-09 2.1867e-10 +4.4531e-02 3.2050e-09 2.2126e-10 +4.4582e-02 3.1980e-09 2.2147e-10 +4.4632e-02 3.1911e-09 2.1988e-10 +4.4683e-02 3.1863e-09 2.1896e-10 +4.4734e-02 3.1820e-09 2.1947e-10 +4.4784e-02 3.1765e-09 2.1981e-10 +4.4835e-02 3.1702e-09 2.1730e-10 +4.4886e-02 3.1633e-09 2.1085e-10 +4.4936e-02 3.1566e-09 2.0392e-10 +4.4987e-02 3.1520e-09 2.0215e-10 +4.5038e-02 3.1497e-09 2.0709e-10 +4.5088e-02 3.1470e-09 2.1456e-10 +4.5139e-02 3.1415e-09 2.1981e-10 +4.5190e-02 3.1331e-09 2.2202e-10 +4.5240e-02 3.1237e-09 2.2244e-10 +4.5291e-02 3.1155e-09 2.2120e-10 +4.5342e-02 3.1094e-09 2.1719e-10 +4.5392e-02 3.1053e-09 2.1139e-10 +4.5443e-02 3.1032e-09 2.0686e-10 +4.5494e-02 3.1028e-09 2.0437e-10 +4.5544e-02 3.1026e-09 2.0244e-10 +4.5595e-02 3.1007e-09 2.0076e-10 +4.5646e-02 3.0972e-09 2.0160e-10 +4.5696e-02 3.0944e-09 2.0665e-10 +4.5747e-02 3.0935e-09 2.1307e-10 +4.5797e-02 3.0922e-09 2.1539e-10 +4.5848e-02 3.0864e-09 2.1196e-10 +4.5899e-02 3.0761e-09 2.0684e-10 +4.5949e-02 3.0661e-09 2.0573e-10 +4.6000e-02 3.0627e-09 2.1156e-10 +4.6051e-02 3.0660e-09 2.2274e-10 +4.6101e-02 3.0699e-09 2.3348e-10 +4.6152e-02 3.0699e-09 2.3801e-10 +4.6203e-02 3.0666e-09 2.3583e-10 +4.6254e-02 3.0613e-09 2.3039e-10 +4.6304e-02 3.0550e-09 2.2505e-10 +4.6355e-02 3.0524e-09 2.2190e-10 +4.6406e-02 3.0568e-09 2.2122e-10 +4.6456e-02 3.0641e-09 2.1983e-10 +4.6507e-02 3.0672e-09 2.1496e-10 +4.6558e-02 3.0644e-09 2.0955e-10 +4.6608e-02 3.0591e-09 2.0637e-10 +4.6659e-02 3.0542e-09 2.0685e-10 +4.6709e-02 3.0505e-09 2.1283e-10 +4.6760e-02 3.0466e-09 2.2408e-10 +4.6811e-02 3.0403e-09 2.3444e-10 +4.6861e-02 3.0311e-09 2.3887e-10 +4.6912e-02 3.0229e-09 2.3867e-10 +4.6963e-02 3.0186e-09 2.3787e-10 +4.7013e-02 3.0176e-09 2.3918e-10 +4.7064e-02 3.0187e-09 2.4244e-10 +4.7115e-02 3.0220e-09 2.4882e-10 +4.7165e-02 3.0256e-09 2.5732e-10 +4.7216e-02 3.0282e-09 2.6565e-10 +4.7267e-02 3.0319e-09 2.7418e-10 +4.7317e-02 3.0422e-09 2.8485e-10 +4.7368e-02 3.0615e-09 2.9691e-10 +4.7419e-02 3.0840e-09 3.0592e-10 +4.7469e-02 3.0991e-09 3.0695e-10 +4.7520e-02 3.1025e-09 3.0121e-10 +4.7571e-02 3.1001e-09 2.9781e-10 +4.7621e-02 3.1015e-09 3.0394e-10 +4.7672e-02 3.1140e-09 3.1649e-10 +4.7723e-02 3.1363e-09 3.2937e-10 +4.7773e-02 3.1548e-09 3.4041e-10 +4.7824e-02 3.1546e-09 3.4676e-10 +4.7875e-02 3.1367e-09 3.4682e-10 +4.7925e-02 3.1190e-09 3.4037e-10 +4.7976e-02 3.1172e-09 3.2906e-10 +4.8027e-02 3.1323e-09 3.1567e-10 +4.8077e-02 3.1536e-09 3.0249e-10 +4.8128e-02 3.1680e-09 2.9581e-10 +4.8179e-02 3.1668e-09 2.9908e-10 +4.8229e-02 3.1535e-09 3.1154e-10 +4.8280e-02 3.1441e-09 3.3428e-10 +4.8331e-02 3.1537e-09 3.6266e-10 +4.8381e-02 3.1818e-09 3.8088e-10 +4.8432e-02 3.2156e-09 3.8299e-10 +4.8483e-02 3.2430e-09 3.7764e-10 +4.8533e-02 3.2577e-09 3.7077e-10 +4.8584e-02 3.2587e-09 3.6607e-10 +4.8635e-02 3.2488e-09 3.6937e-10 +4.8685e-02 3.2329e-09 3.7890e-10 +4.8736e-02 3.2159e-09 3.8079e-10 +4.8786e-02 3.2012e-09 3.6826e-10 +4.8837e-02 3.1880e-09 3.5452e-10 +4.8888e-02 3.1758e-09 3.5573e-10 +4.8939e-02 3.1675e-09 3.6734e-10 +4.8989e-02 3.1631e-09 3.6076e-10 +4.9040e-02 3.1555e-09 3.2724e-10 +4.9091e-02 3.1413e-09 2.9474e-10 +4.9141e-02 3.1292e-09 2.8679e-10 +4.9192e-02 3.1276e-09 3.0548e-10 +4.9243e-02 3.1306e-09 3.2941e-10 +4.9293e-02 3.1258e-09 3.3821e-10 +4.9344e-02 3.1125e-09 3.3586e-10 +4.9395e-02 3.1013e-09 3.3160e-10 +4.9445e-02 3.1017e-09 3.3532e-10 +4.9496e-02 3.1171e-09 3.5592e-10 +4.9547e-02 3.1442e-09 3.8322e-10 +4.9597e-02 3.1697e-09 3.8795e-10 +4.9648e-02 3.1774e-09 3.5353e-10 +4.9698e-02 3.1652e-09 3.0229e-10 +4.9749e-02 3.1458e-09 2.6921e-10 +4.9800e-02 3.1279e-09 2.6114e-10 +4.9850e-02 3.1065e-09 2.6349e-10 +4.9901e-02 3.0755e-09 2.6692e-10 +4.9952e-02 3.0392e-09 2.7335e-10 +5.0002e-02 3.0077e-09 2.8960e-10 +5.0053e-02 2.9886e-09 3.1572e-10 +5.0104e-02 2.9836e-09 3.3542e-10 +5.0154e-02 2.9871e-09 3.2696e-10 +5.0205e-02 2.9869e-09 2.9247e-10 +5.0256e-02 2.9735e-09 2.6090e-10 +5.0306e-02 2.9479e-09 2.4489e-10 +5.0357e-02 2.9182e-09 2.3522e-10 +5.0408e-02 2.8903e-09 2.2500e-10 +5.0458e-02 2.8667e-09 2.1541e-10 +5.0509e-02 2.8509e-09 2.1445e-10 +5.0560e-02 2.8450e-09 2.2186e-10 +5.0610e-02 2.8458e-09 2.3093e-10 diff --git a/examples/thermostats/data/water_32.pdb b/examples/thermostats/data/water_32.pdb new file mode 100644 index 00000000..014e40ff --- /dev/null +++ b/examples/thermostats/data/water_32.pdb @@ -0,0 +1,98 @@ +CRYST1 10.260 10.260 10.260 90.00 90.00 90.00 P 1 +ATOM 1 O 1 1 3.756 4.710 9.494 0.00 0.00 O +ATOM 2 H 1 1 4.604 4.272 9.671 0.00 0.00 H +ATOM 3 H 1 1 3.998 5.320 8.788 0.00 0.00 H +ATOM 4 O 1 1 9.933 8.841 0.366 0.00 0.00 O +ATOM 5 H 1 1 10.132 8.196 1.120 0.00 0.00 H +ATOM 6 H 1 1 9.368 8.449 -0.316 0.00 0.00 H +ATOM 7 O 1 1 0.321 1.492 5.796 0.00 0.00 O +ATOM 8 H 1 1 -0.287 1.993 5.241 0.00 0.00 H +ATOM 9 H 1 1 0.791 2.061 6.364 0.00 0.00 H +ATOM 10 O 1 1 8.035 9.735 4.307 0.00 0.00 O +ATOM 11 H 1 1 7.203 9.789 4.847 0.00 0.00 H +ATOM 12 H 1 1 8.636 9.307 4.920 0.00 0.00 H +ATOM 13 O 1 1 5.663 9.082 0.660 0.00 0.00 O +ATOM 14 H 1 1 6.378 9.721 0.814 0.00 0.00 H +ATOM 15 H 1 1 5.213 8.991 1.552 0.00 0.00 H +ATOM 16 O 1 1 8.130 0.215 1.201 0.00 0.00 O +ATOM 17 H 1 1 8.196 -0.065 2.118 0.00 0.00 H +ATOM 18 H 1 1 8.938 -0.161 0.818 0.00 0.00 H +ATOM 19 O 1 1 8.177 4.165 0.716 0.00 0.00 O +ATOM 20 H 1 1 7.895 5.066 0.840 0.00 0.00 H +ATOM 21 H 1 1 7.722 3.562 1.341 0.00 0.00 H +ATOM 22 O 1 1 6.341 3.256 9.678 0.00 0.00 O +ATOM 23 H 1 1 7.133 3.423 10.189 0.00 0.00 H +ATOM 24 H 1 1 6.507 2.398 9.350 0.00 0.00 H +ATOM 25 O 1 1 0.136 7.798 2.738 0.00 0.00 O +ATOM 26 H 1 1 -0.006 8.300 3.584 0.00 0.00 H +ATOM 27 H 1 1 0.314 6.905 2.966 0.00 0.00 H +ATOM 28 O 1 1 5.027 2.563 6.169 0.00 0.00 O +ATOM 29 H 1 1 5.538 3.336 6.262 0.00 0.00 H +ATOM 30 H 1 1 5.313 2.031 6.922 0.00 0.00 H +ATOM 31 O 1 1 7.164 2.542 2.418 0.00 0.00 O +ATOM 32 H 1 1 6.248 2.467 2.667 0.00 0.00 H +ATOM 33 H 1 1 7.384 1.681 2.057 0.00 0.00 H +ATOM 34 O 1 1 3.336 9.051 3.265 0.00 0.00 O +ATOM 35 H 1 1 2.818 8.341 2.787 0.00 0.00 H +ATOM 36 H 1 1 2.733 9.794 3.398 0.00 0.00 H +ATOM 37 O 1 1 1.476 1.420 0.819 0.00 0.00 O +ATOM 38 H 1 1 1.060 0.568 0.631 0.00 0.00 H +ATOM 39 H 1 1 1.582 1.509 1.770 0.00 0.00 H +ATOM 40 O 1 1 1.222 4.946 3.218 0.00 0.00 O +ATOM 41 H 1 1 2.111 4.703 3.522 0.00 0.00 H +ATOM 42 H 1 1 1.342 4.956 2.245 0.00 0.00 H +ATOM 43 O 1 1 6.790 6.491 4.488 0.00 0.00 O +ATOM 44 H 1 1 7.083 5.755 5.074 0.00 0.00 H +ATOM 45 H 1 1 6.747 7.293 4.975 0.00 0.00 H +ATOM 46 O 1 1 9.330 3.465 4.430 0.00 0.00 O +ATOM 47 H 1 1 9.974 3.956 3.918 0.00 0.00 H +ATOM 48 H 1 1 8.583 3.143 3.966 0.00 0.00 H +ATOM 49 O 1 1 7.484 4.543 6.379 0.00 0.00 O +ATOM 50 H 1 1 7.604 4.450 7.328 0.00 0.00 H +ATOM 51 H 1 1 8.241 4.098 5.987 0.00 0.00 H +ATOM 52 O 1 1 0.448 5.701 8.219 0.00 0.00 O +ATOM 53 H 1 1 0.573 4.840 7.871 0.00 0.00 H +ATOM 54 H 1 1 0.720 5.651 9.123 0.00 0.00 H +ATOM 55 O 1 1 0.736 4.082 0.545 0.00 0.00 O +ATOM 56 H 1 1 1.032 3.147 0.501 0.00 0.00 H +ATOM 57 H 1 1 -0.250 3.978 0.476 0.00 0.00 H +ATOM 58 O 1 1 4.229 2.582 3.562 0.00 0.00 O +ATOM 59 H 1 1 4.650 2.364 4.424 0.00 0.00 H +ATOM 60 H 1 1 4.306 3.526 3.496 0.00 0.00 H +ATOM 61 O 1 1 3.749 5.318 3.591 0.00 0.00 O +ATOM 62 H 1 1 3.577 5.600 2.690 0.00 0.00 H +ATOM 63 H 1 1 4.519 5.812 3.831 0.00 0.00 H +ATOM 64 O 1 1 0.242 8.540 5.195 0.00 0.00 O +ATOM 65 H 1 1 0.243 9.452 5.454 0.00 0.00 H +ATOM 66 H 1 1 0.725 8.102 5.864 0.00 0.00 H +ATOM 67 O 1 1 6.065 0.243 8.171 0.00 0.00 O +ATOM 68 H 1 1 6.685 -0.360 8.559 0.00 0.00 H +ATOM 69 H 1 1 5.235 -0.020 8.586 0.00 0.00 H +ATOM 70 O 1 1 7.362 8.029 9.049 0.00 0.00 O +ATOM 71 H 1 1 6.719 7.831 8.312 0.00 0.00 H +ATOM 72 H 1 1 6.813 7.839 9.757 0.00 0.00 H +ATOM 73 O 1 1 1.971 1.171 3.631 0.00 0.00 O +ATOM 74 H 1 1 1.520 1.399 4.474 0.00 0.00 H +ATOM 75 H 1 1 2.821 1.577 3.533 0.00 0.00 H +ATOM 76 O 1 1 7.674 6.605 2.134 0.00 0.00 O +ATOM 77 H 1 1 7.143 6.490 2.918 0.00 0.00 H +ATOM 78 H 1 1 8.563 6.681 2.409 0.00 0.00 H +ATOM 79 O 1 1 2.769 7.384 7.275 0.00 0.00 O +ATOM 80 H 1 1 2.353 8.088 7.748 0.00 0.00 H +ATOM 81 H 1 1 2.249 6.554 7.408 0.00 0.00 H +ATOM 82 O 1 1 2.253 2.791 7.349 0.00 0.00 O +ATOM 83 H 1 1 2.718 3.230 8.076 0.00 0.00 H +ATOM 84 H 1 1 2.955 2.671 6.707 0.00 0.00 H +ATOM 85 O 1 1 3.652 9.953 9.290 0.00 0.00 O +ATOM 86 H 1 1 4.114 9.509 9.983 0.00 0.00 H +ATOM 87 H 1 1 3.199 10.639 9.784 0.00 0.00 H +ATOM 88 O 1 1 5.263 6.708 7.589 0.00 0.00 O +ATOM 89 H 1 1 4.318 7.082 7.707 0.00 0.00 H +ATOM 90 H 1 1 5.455 7.057 6.733 0.00 0.00 H +ATOM 91 O 1 1 3.507 6.721 1.020 0.00 0.00 O +ATOM 92 H 1 1 4.097 7.303 0.607 0.00 0.00 H +ATOM 93 H 1 1 3.586 5.984 0.419 0.00 0.00 H +ATOM 94 O 1 1 5.455 9.196 5.888 0.00 0.00 O +ATOM 95 H 1 1 5.473 9.602 6.758 0.00 0.00 H +ATOM 96 H 1 1 4.549 9.238 5.484 0.00 0.00 H +END diff --git a/examples/thermostats/data/water_32_data.lmp b/examples/thermostats/data/water_32_data.lmp new file mode 100644 index 00000000..385c7d7c --- /dev/null +++ b/examples/thermostats/data/water_32_data.lmp @@ -0,0 +1,227 @@ +LAMMPS Description + + 96 atoms + 64 bonds + 32 angles + + 2 atom types + 1 bond types + 1 angle types + + -9.694295 9.694295 xlo xhi + -9.694295 9.694295 ylo yhi + -9.694295 9.694295 zlo zhi + +Masses + + 1 15.9994 + 2 1.0080 + +Bond Coeffs + + 1 1.78 0.2708585 -0.327738785 0.231328959 + +Angle Coeffs + + 1 0.0700 107.400000 + +Atoms + + 1 1 1 -1.1128 7.098 8.901 17.941 + 2 1 2 0.5564 8.700 8.073 18.276 + 3 1 2 0.5564 7.555 10.053 16.607 + 4 2 1 -1.1128 18.771 16.707 0.692 + 5 2 2 0.5564 19.147 15.488 2.116 + 6 2 2 0.5564 17.703 15.966 -0.597 + 7 3 1 -1.1128 0.607 2.819 10.953 + 8 3 2 0.5564 -0.542 3.766 9.904 + 9 3 2 0.5564 1.495 3.895 12.026 + 10 4 1 -1.1128 15.184 18.396 8.139 + 11 4 2 0.5564 13.612 18.499 9.160 + 12 4 2 0.5564 16.320 17.588 9.297 + 13 5 1 -1.1128 10.702 17.162 1.247 + 14 5 2 0.5564 12.053 18.370 1.538 + 15 5 2 0.5564 9.851 16.991 2.933 + 16 6 1 -1.1128 15.363 0.406 2.270 + 17 6 2 0.5564 15.488 -0.123 4.002 + 18 6 2 0.5564 16.890 -0.304 1.546 + 19 7 1 -1.1128 15.452 7.871 1.353 + 20 7 2 0.5564 14.919 9.573 1.587 + 21 7 2 0.5564 14.592 6.731 2.534 + 22 8 1 -1.1128 11.983 6.153 18.289 + 23 8 2 0.5564 13.479 6.469 19.254 + 24 8 2 0.5564 12.296 4.532 17.669 + 25 9 1 -1.1128 0.257 14.736 5.174 + 26 9 2 0.5564 -0.011 15.685 6.773 + 27 9 2 0.5564 0.593 13.049 5.605 + 28 10 1 -1.1128 9.500 4.843 11.658 + 29 10 2 0.5564 10.465 6.304 11.833 + 30 10 2 0.5564 10.040 3.838 13.081 + 31 11 1 -1.1128 13.538 4.804 4.569 + 32 11 2 0.5564 11.807 4.662 5.040 + 33 11 2 0.5564 13.954 3.177 3.887 + 34 12 1 -1.1128 6.304 17.104 6.170 + 35 12 2 0.5564 5.325 15.762 5.267 + 36 12 2 0.5564 5.165 18.508 6.421 + 37 13 1 -1.1128 2.789 2.683 1.548 + 38 13 2 0.5564 2.003 1.073 1.192 + 39 13 2 0.5564 2.990 2.852 3.345 + 40 14 1 -1.1128 2.309 9.347 6.081 + 41 14 2 0.5564 3.989 8.887 6.656 + 42 14 2 0.5564 2.536 9.365 4.242 + 43 15 1 -1.1128 12.831 12.266 8.481 + 44 15 2 0.5564 13.385 10.875 9.588 + 45 15 2 0.5564 12.750 13.782 9.401 + 46 16 1 -1.1128 17.631 6.548 8.371 + 47 16 2 0.5564 18.848 7.476 7.404 + 48 16 2 0.5564 16.220 5.939 7.495 + 49 17 1 -1.1128 14.143 8.585 12.055 + 50 17 2 0.5564 14.369 8.409 13.848 + 51 17 2 0.5564 15.573 7.744 11.314 + 52 18 1 -1.1128 0.847 10.773 15.532 + 53 18 2 0.5564 1.083 9.146 14.874 + 54 18 2 0.5564 1.361 10.679 17.240 + 55 19 1 -1.1128 1.391 7.714 1.030 + 56 19 2 0.5564 1.950 5.947 0.947 + 57 19 2 0.5564 -0.472 7.517 0.900 + 58 20 1 -1.1128 7.992 4.879 6.731 + 59 20 2 0.5564 8.787 4.467 8.360 + 60 20 2 0.5564 8.137 6.663 6.606 + 61 21 1 -1.1128 7.085 10.050 6.786 + 62 21 2 0.5564 6.760 10.582 5.083 + 63 21 2 0.5564 8.540 10.983 7.240 + 64 22 1 -1.1128 0.457 16.138 9.817 + 65 22 2 0.5564 0.459 17.862 10.307 + 66 22 2 0.5564 1.370 15.311 11.081 + 67 23 1 -1.1128 11.461 0.459 15.441 + 68 23 2 0.5564 12.633 -0.680 16.174 + 69 23 2 0.5564 9.893 -0.038 16.225 + 70 24 1 -1.1128 13.912 15.173 17.100 + 71 24 2 0.5564 12.697 14.798 15.707 + 72 24 2 0.5564 12.875 14.814 18.438 + 73 25 1 -1.1128 3.725 2.213 6.862 + 74 25 2 0.5564 2.872 2.644 8.455 + 75 25 2 0.5564 5.331 2.980 6.676 + 76 26 1 -1.1128 14.502 12.482 4.033 + 77 26 2 0.5564 13.498 12.264 5.514 + 78 26 2 0.5564 16.182 12.625 4.552 + 79 27 1 -1.1128 5.233 13.954 13.748 + 80 27 2 0.5564 4.447 15.284 14.642 + 81 27 2 0.5564 4.250 12.385 13.999 + 82 28 1 -1.1128 4.258 5.274 13.888 + 83 28 2 0.5564 5.136 6.104 15.261 + 84 28 2 0.5564 5.584 5.047 12.674 + 85 29 1 -1.1128 6.901 18.808 17.556 + 86 29 2 0.5564 7.774 17.969 18.865 + 87 29 2 0.5564 6.045 20.105 18.489 + 88 30 1 -1.1128 9.946 12.676 14.341 + 89 30 2 0.5564 8.160 13.383 14.564 + 90 30 2 0.5564 10.308 13.336 12.724 + 91 31 1 -1.1128 6.627 12.701 1.928 + 92 31 2 0.5564 7.742 13.801 1.147 + 93 31 2 0.5564 6.777 11.308 0.792 + 94 32 1 -1.1128 10.308 17.378 11.127 + 95 32 2 0.5564 10.342 18.145 12.771 + 96 32 2 0.5564 8.596 17.457 10.363 + +Bonds + + 1 1 1 2 + 2 1 1 3 + 3 1 4 5 + 4 1 4 6 + 5 1 7 8 + 6 1 7 9 + 7 1 10 11 + 8 1 10 12 + 9 1 13 14 + 10 1 13 15 + 11 1 16 17 + 12 1 16 18 + 13 1 19 20 + 14 1 19 21 + 15 1 22 23 + 16 1 22 24 + 17 1 25 26 + 18 1 25 27 + 19 1 28 29 + 20 1 28 30 + 21 1 31 32 + 22 1 31 33 + 23 1 34 35 + 24 1 34 36 + 25 1 37 38 + 26 1 37 39 + 27 1 40 41 + 28 1 40 42 + 29 1 43 44 + 30 1 43 45 + 31 1 46 47 + 32 1 46 48 + 33 1 49 50 + 34 1 49 51 + 35 1 52 53 + 36 1 52 54 + 37 1 55 56 + 38 1 55 57 + 39 1 58 59 + 40 1 58 60 + 41 1 61 62 + 42 1 61 63 + 43 1 64 65 + 44 1 64 66 + 45 1 67 68 + 46 1 67 69 + 47 1 70 71 + 48 1 70 72 + 49 1 73 74 + 50 1 73 75 + 51 1 76 77 + 52 1 76 78 + 53 1 79 80 + 54 1 79 81 + 55 1 82 83 + 56 1 82 84 + 57 1 85 86 + 58 1 85 87 + 59 1 88 89 + 60 1 88 90 + 61 1 91 92 + 62 1 91 93 + 63 1 94 95 + 64 1 94 96 + +Angles + + 1 1 2 1 3 + 2 1 5 4 6 + 3 1 8 7 9 + 4 1 11 10 12 + 5 1 14 13 15 + 6 1 17 16 18 + 7 1 20 19 21 + 8 1 23 22 24 + 9 1 26 25 27 + 10 1 29 28 30 + 11 1 32 31 33 + 12 1 35 34 36 + 13 1 38 37 39 + 14 1 41 40 42 + 15 1 44 43 45 + 16 1 47 46 48 + 17 1 50 49 51 + 18 1 53 52 54 + 19 1 56 55 57 + 20 1 59 58 60 + 21 1 62 61 63 + 22 1 65 64 66 + 23 1 68 67 69 + 24 1 71 70 72 + 25 1 74 73 75 + 26 1 77 76 78 + 27 1 80 79 81 + 28 1 83 82 84 + 29 1 86 85 87 + 30 1 89 88 90 + 31 1 92 91 93 + 32 1 95 94 96 diff --git a/examples/thermostats/environment.yml b/examples/thermostats/environment.yml new file mode 100644 index 00000000..6d18827e --- /dev/null +++ b/examples/thermostats/environment.yml @@ -0,0 +1,12 @@ +channels: + - conda-forge +dependencies: + - python=3.11 + - pip + - lammps + - pip: + - ase==3.22.1 + - chemiscope>=0.7 + - git+https://github.com/i-pi/i-pi.git@main + - matplotlib + - skmatter diff --git a/examples/thermostats/thermostats.py b/examples/thermostats/thermostats.py new file mode 100644 index 00000000..2a282aa2 --- /dev/null +++ b/examples/thermostats/thermostats.py @@ -0,0 +1,902 @@ +""" +Constant-temperature MD and thermostats +======================================= + +:Authors: Michele Ceriotti `@ceriottm `_ + +This recipe gives a practical introduction to finite-temperature +molecular dynamics simulations, and provides a guide to choose the +most appropriate thermostat for the simulation at hand. + +As for other examples in the cookbook, a small simulation of liquid +water is used as an archetypal example. Molecular dynamics, sampling, +and constant-temperature simulations are discussed in much detail in +the book "Understanding Molecular Simulations" by Daan Frenkel and Berend Smit. +This +`seminal paper by H.C.Andersen `_ +provides a good historical introduction to the problem of +thermostatting, and this +`PhD thesis +`_ +provides a more detailed background to several of the techniques +discussed in this recipe. +""" + +import os + +# %% +import subprocess +import time +import xml.etree.ElementTree as ET + +import chemiscope +import ipi +import matplotlib as mpl +import matplotlib.pyplot as plt +import numpy as np +from ipi.utils.tools.acf_xyz import compute_acf_xyz +from ipi.utils.tools.gle import get_gle_matrices, gle_frequency_kernel, isra_deconvolute + + +# %% +# Constant-temperature sampling of (thermo)dynamics +# ------------------------------------------------- +# +# Even though Hamilton's equations in classical mechanics conserve the total +# energy of the group of atoms in a simulation, experimental boundary conditions +# usually involve exchange of heat with the surroundings, especially when considering +# the relatively small supercells that are often used in simulations. +# +# The goal of a constant-temperature MD simulation is to compute efficiently thermal +# averages of the form :math:`\langle A(q,p)\rangle_\beta`, where the average +# of the observable :math:`A(q,p)` is +# evaluated over the Boltzmann distribution at inverse temperature +# :math:`\beta=1/k_\mathrm{B}T`, +# :math:`P(q,p)=Q^{-1} \exp(-\beta(p^2/2m + V(q)))` +# In all these scenarios, optimizing the simulation involves reducing as much as +# possible the *autocorrelation time* of the observable. +# +# Constant-temperature sampling is also important when one wants to compute +# *dynamical* properties. In principle these would require +# constant-energy trajectories, as any thermostatting procedure modifies +# the dynamics of the system. However, the initial conditions +# should usually be determined from constant-temperature conditions, +# averaging over multiple constant-energy trajectories. +# As we shall see, this protocol can often be simplified greatly, by choosing +# thermostats that don't interfere with the natural microscopic dynamics. + +# %% +# Running simulations +# ~~~~~~~~~~~~~~~~~~~ +# +# We use `i-PI `_ together with a ``LAMMPS`` driver to run +# all the simulations in this recipe. The two codes need to be run separately, +# and communicate atomic positions, energy and forces through a socket interface. +# +# The LAMMPS input defines the parameters of the +# `q-TIP4P/f water model `_, +# while the XML-formatted input of i-PI describes the setup of the +# MD simulation. +# +# We begin running a constant-energy calculation, that +# we will use to illustrate the metrics that can be applied to +# assess the performance of a thermostatting scheme. If it is the +# first time you see an ``i-PI`` input, you may want to look at +# the input file side-by-sidewith the +# `input reference `_. + +# Open and read the XML file +with open("data/input_nve.xml", "r") as file: + xml_content = file.read() +print(xml_content) + +# %% +# The part of the input that describes the molecular dynamics integrator +# is the ``motion`` class. For this run, it specifies an *NVE* ensemble, and +# a ``timestep`` of 1 fs for the integrator. + +xmlroot = ET.parse("data/input_nve.xml").getroot() +print(" " + ET.tostring(xmlroot.find(".//motion"), encoding="unicode")) + +# %% +# Note that this -- and other runs in this example -- are too short to +# provide quantitative results, and you may want to increase the +# ```` parameter so that the simulation runs for at least +# a few tens of ps. The time step of 1 fs is also at the limit of what +# is acceptable for running simulations of water. 0.5 fs would be a +# safer, stabler value. + + +# %% +# To launch i-PI and LAMMPS from the command line you can just +# execute the following commands +# +# .. code-block:: bash +# +# i-pi data/input_nve.xml > log & +# sleep 2 +# lmp -in data/in.lmp & +# +# To launch the external processes from a Python script +# proceed as follows: + +ipi_process = None +if not os.path.exists("simulation_nve.out"): + ipi_process = subprocess.Popen(["i-pi", "data/input_nve.xml"]) + time.sleep(4) # wait for i-PI to start + lmp_process = [subprocess.Popen(["lmp", "-in", "data/in.lmp"]) for i in range(1)] + +# %% +# If you run this in a notebook, you can go ahead and start loading +# output files *before* i-PI and LAMMPS have finished running, by +# skipping this cell + +if ipi_process is not None: + ipi_process.wait() + lmp_process[0].wait() + + +# %% +# Analyzing the simulation +# ~~~~~~~~~~~~~~~~~~~~~~~~ +# +# After the simulation is finished, we can look at the outputs. +# The outputs include the trajectory of positions, the velocities +# and a number of energetic observables + +output_data, output_desc = ipi.read_output("simulation_nve.out") +traj_data = ipi.read_trajectory("simulation_nve.pos_0.xyz") + +# %% +# The trajectory shows mostly local vibrations on this short time scale, +# but if you re-run with a longer ```` settings you should be +# able to observe diffusing molecules in the liquid. + +chemiscope.show( + traj_data, + mode="structure", + settings=chemiscope.quick_settings( + trajectory=True, structure_settings={"unitCell": True} + ), +) + +# %% +# Potential and kinetic energy fluctuate, but the total energy is +# (almost) constant, the small fluctuations being due to integration +# errors, that are quite large with the long time step used for this +# example. If you run with smaller ```` values, you should +# see that the energy conservation condition is fulfilled with higher +# accuracy. + +fig, ax = plt.subplots(1, 1, figsize=(4, 3), constrained_layout=True) +ax.plot( + output_data["time"], + output_data["potential"] - output_data["potential"][0], + "b-", + label="Potential, $V$", +) +ax.plot( + output_data["time"], + output_data["kinetic_md"], + "r-", + label="Kinetic, $K$", +) +ax.plot( + output_data["time"], + output_data["conserved"] - output_data["conserved"][0], + "k-", + label="Conserved, $H$", +) +ax.set_xlabel(r"$t$ / ps") +ax.set_ylabel(r"energy / eV") +ax.legend() +plt.show() + +# %% +# In a classical MD simulation, based on the momentum :math:`\mathbf{p}` +# of each atom, it is possible to evaluate its *kinetic temperature +# estimator* :math:`T=\langle \mathbf{p}^2/m \rangle /3k_B` the average is to +# be intended over a converged trajectory. Keep in mind that +# +# 1. The *instantaneous* value of this estimator is meaningless +# 2. It is only well-defined in a constant-temperature simulation, so here +# it only gives a sense of whether atomic momenta are close to what one +# would expect at 300 K. +# +# With these caveats in mind, we can observe that the simulation has higher +# velocities than expected at 300 K, and that there is no equipartition, the +# O atoms having on average a higher energy than the H atoms. + +fig, ax = plt.subplots(1, 1, figsize=(4, 3), constrained_layout=True) +ax.plot( + output_data["time"], + output_data["temperature(O)"], + "r-", + label="O atoms", +) +ax.plot( + output_data["time"], + output_data["temperature(H)"], + "c-", + label="H atoms", +) +ax.plot(output_data["time"], output_data["temperature"], "k-", label="All atoms") +ax.set_xlabel(r"$t$ / ps") +ax.set_ylabel(r"$\tilde{T}$ / K") +plt.show() + +# %% +# In order to investigate the dynamics more carefully, we +# compute the velocity-velocity autocorrelation function +# :math:`c_{vv}(t)=\sum_i \langle \mathbf{v}_i(t) \cdot \mathbf{v}_i(0) \rangle`. +# We use a utility function that reads the outputs of ``i-PI`` +# and computes both the autocorrelation function and its Fourier +# transform. +# :math:`c_{vv}(t)` contains information on the time scale and amplitude +# of molecular motion, and is closely related to the vibrational density +# of states and to spectroscopic observables such as IR and Raman spectra. + +acf_nve = compute_acf_xyz( + "simulation_nve.vel_0.xyz", + maximum_lag=600, + length_zeropadding=2000, + spectral_windowing="cosine-blackman", + timestep=1, + time_units="femtosecond", + skip=100, +) + +fig, ax = plt.subplots(1, 1, figsize=(4, 3), constrained_layout=True) +ax.plot( + acf_nve[0][:1200] * 2.4188843e-05, # atomic time to ps + acf_nve[1][:1200] * 1e5, + "r-", +) +ax.set_xlabel(r"$t$ / ps$") +ax.set_ylabel(r"$c_{vv}$ / arb. units") +ax.legend() +plt.show() + +# %% +# The power spectrum (that can be computed as the Fourier transform of +# :math:`c_{vv}`) reveals the frequencies of stretching, bending and libration +# modes of water; the :math:`\omega\rightarrow 0` limit is proportional +# to the diffusion coefficient. +# We also load the results from a reference calculation (average of 8 +# trajectories initiated from NVT-equilibrated samples, shown as the +# confidence interval). You can see how to run these reference calculations +# from the script ``data/run_traj.sh``. +# The differences are due to the short trajectory, and to the fact that the +# NVE trajectory is not equilibrated at 300 K. + +ha2cm1 = 219474.63 + +# Loads reference trajectory +acf_ref = np.loadtxt("data/traj-all_facf.data") + +fig, ax = plt.subplots(1, 1, figsize=(4, 3), constrained_layout=True) + +ax.fill_between( + acf_ref[:1200, 0] * ha2cm1, + (acf_ref[:1200, 1] - acf_ref[:1200, 2]) * 1e5, + (acf_ref[:1200, 1] + acf_ref[:1200, 2]) * 1e5, + color="gray", + label="reference", +) + +ax.loglog(acf_nve[3][:1200] * ha2cm1, acf_nve[4][:1200] * 1e5, "r-", label="NVE") +ax.set_xlabel(r"$\omega$ / cm$^{-1}$") +ax.set_ylabel(r"$\hat{c}_{vv}$ / arb. units") +ax.legend() +plt.show() + +# %% +# Langevin thermostatting +# ----------------------- +# +# In order to perform a simulations that samples configurations +# consistent with a Boltzmann distribution :math:`e^{-V(x)/k_B T}` +# one needs to modify the equations of motion. There are many +# different approaches to do this, some of which lead to deterministic +# dynamics; the two more widely used deterministic thermostats +# are the +# `Berendsen thermostat `_ +# which does not sample the Boltzmann distribution exactly and +# should never be used given the many more rigorous alternatives, +# and the Nosé-Hoover thermostat, that requires a +# `"chain" implementation `_ +# to be ergodic, which amounts essentially to a complicated way +# to generate poor-quality pseudo-random numbers. +# +# Given the limitations of deterministic thermostats, in this +# recipe we focus on stochastic thermostats, that model the +# coupling to the chaotic dynamics of an external bath through +# explicit random numbers. Langevin dynamics amounts to adding +# to Hamilton's equations of motion, for each degree of freedom, +# a term of the form +# +# .. math:: +# +# \dot{p} = -\gamma p + \sqrt{2\gamma m k_B T} \, \xi(t) +# +# where :math:` +# \gamma` is a friction coefficient, and :math:`\xi` +# uncorrelated random numbers that mimic collisions with the bath +# particles. The friction can be seen as the inverse of a +# characteristic *coupling time scale* +# :math:`\tau=1/\gamma` that describes how strongly the bath +# interacts with the system. + +# %% +# Setting up a thermostat in ``i-PI`` +# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +# +# In order to set up a thermostat in ``i-PI``, one simply needs +# to adjust the ```` block, to perform ``nvt`` dynamics +# and include an appropriate ```` section. +# Here we use a very-strongly coupled Langevin thermostat, +# with :math:`\tau=10~fs`. + +xmlroot = ET.parse("data/input_higamma.xml").getroot() +print(" " + ET.tostring(xmlroot.find(".//dynamics"), encoding="unicode")) + +# %% +# ``i-PI`` and ``LAMMPS`` are launched as above ... + +ipi_process = None +if not os.path.exists("simulation_higamma.out"): + ipi_process = subprocess.Popen(["i-pi", "data/input_higamma.xml"]) + time.sleep(4) # wait for i-PI to start + lmp_process = [subprocess.Popen(["lmp", "-in", "data/in.lmp"]) for i in range(1)] + +# %% +# ... and you should probably wait until they're done, +# it'll take less than a minute. + +if ipi_process is not None: + ipi_process.wait() + lmp_process[0].wait() + +# %% +# Analysis of the trajectory +# ~~~~~~~~~~~~~~~~~~~~~~~~~~ +# +# The temperature converges very quickly to the target value +# (fluctuations are to be expected, given that as discussed above +# the temperature estimator is just the instantaneous kinetic energy, +# that is not constant). There is also equipartition between O and H. + +output_data, output_desc = ipi.read_output("simulation_higamma.out") +traj_data = ipi.read_trajectory("simulation_higamma.pos_0.xyz") + +fig, ax = plt.subplots(1, 1, figsize=(4, 3), constrained_layout=True) +ax.plot( + output_data["time"], + output_data["temperature(O)"], + "r-", + label="O atoms", +) +ax.plot( + output_data["time"], + output_data["temperature(H)"], + "c-", + label="H atoms", +) +ax.plot(output_data["time"], output_data["temperature"], "k-", label="All atoms") +ax.set_xlabel(r"$t$ / ps") +ax.set_ylabel(r"$\tilde{T}$ / K") +ax.legend() +plt.show() + + +# %% +# The velocity-velocity correlation function shows how much +# this thermostat affects the system dynamics. The high-frequency peaks, +# corresponding to stretches and bending, are +# greatly broadened, and the :math:`\omega\rightarrow 0` +# limit of :math:`\hat{c}_{vv}`, corresponding to the +# diffusion coefficient, is reduced by almost a factor of 5. +# This last observation highlights that a too-aggressive +# thermostat is not only disrupting the dynamics: +# it also slows down diffusion through phase space, +# making the dynamics less efficient at sampling slow, +# collective motions. We shall see further down various +# methods to counteract this effect, but in general one should +# use a weaker coupling, that improves the sampling of configuration +# space even though it slows down the convergence of the +# kinetic energy. If you want a thermostat that equilibrates +# aggressively the temperature while disturbing less the diffusive +# modes, you may try the *fast-forward Langevin* thermostat +# `(Hijazi et al., JCP (2018)) `_ +# that can be activated with the option ``mode="ffl"``. + +# compute the v-v acf +acf_higamma = compute_acf_xyz( + "simulation_higamma.vel_0.xyz", + maximum_lag=600, + length_zeropadding=2000, + spectral_windowing="cosine-blackman", + timestep=1, + time_units="femtosecond", + skip=100, +) + +# and plot +fig, ax = plt.subplots(1, 1, figsize=(4, 3), constrained_layout=True) +ax.fill_between( + acf_ref[:1200, 0] * ha2cm1, + (acf_ref[:1200, 1] - acf_ref[:1200, 2]) * 1e5, + (acf_ref[:1200, 1] + acf_ref[:1200, 2]) * 1e5, + color="gray", + label="reference", +) +ax.loglog( + acf_higamma[3][:1200] * ha2cm1, + acf_higamma[4][:1200] * 1e5, + "b-", + label=r"Langevin, $\tau=10$fs", +) +ax.set_xlabel(r"$\omega$ / cm$^{-1}$") +ax.set_ylabel(r"$\hat{c}_{vv}$ / arb. units") +ax.legend() +plt.show() + +# %% +# Global thermostats: stochastic velocity rescaling +# ------------------------------------------------- +# +# An alternative approach to sample the canonical Boltzmann +# distribution while introducing fewer disturbances to the system +# dynamics is to use a *global* thermostat, i.e. a scheme that +# targets the *total* kinetic energy of the system, rather than that +# of individual degrees of freedom. +# We recommend the "stochastic velocity rescaling" thermostat +# `(Bussi, Donadio, Parrinello, JCP (2007)) `_ +# that acts by rescaling the total momentum vector, adding a +# suitably distributed random noise term. + + +# %% +# Setting up a thermostat in ``i-PI`` +# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +# +# Stochastic velocity rescaling is implemented in ``i-PI`` +# can be selected by setting ``mode="svr"``, and has a +# ``tau`` parameter that corresponds to the time scale of the +# coupling. + +xmlroot = ET.parse("data/input_svr.xml").getroot() +print(" " + ET.tostring(xmlroot.find(".//thermostat"), encoding="unicode")) + +# %% +# We run a simulation with the usual set up ... + +ipi_process = None +if not os.path.exists("simulation_svr.out"): + ipi_process = subprocess.Popen(["i-pi", "data/input_svr.xml"]) + time.sleep(4) # wait for i-PI to start + lmp_process = [subprocess.Popen(["lmp", "-in", "data/in.lmp"]) for i in range(1)] + +# %% +# ... and wait for it to finish. + +if ipi_process is not None: + ipi_process.wait() + lmp_process[0].wait() + +# %% +# Analysis of the trajectory +# ~~~~~~~~~~~~~~~~~~~~~~~~~~ +# +# The kinetic temperature of the trajectory equilibrates very +# rapidly to the target value. However, it takes a bit longer +# (approximately 0.5 ps) to reach equipartition between O and H +# atoms. This is an important shortcoming of global thermostats: +# since they only target the total kinetic energy, they must rely +# on internal energy redistribution to reach equilibrium between +# different degrees of freedom. +# Liquid water is a very ergodic system, in which all degrees of +# freedom are strongly coupled, so this is not a major issue. However +# care must be taken when modeling a quasi-harmonic crystal (e.g. +# diamond, a metal, or an inorganic crystal), or a molecular system +# in which the coupling between molecules is weaker (e.g. methane, +# or another apolar compound). + +output_data, output_desc = ipi.read_output("simulation_svr.out") +traj_data = ipi.read_trajectory("simulation_svr.pos_0.xyz") + +fig, ax = plt.subplots(1, 1, figsize=(4, 3), constrained_layout=True) +ax.plot( + output_data["time"], + output_data["temperature(O)"], + "r-", + label="O atoms", +) +ax.plot( + output_data["time"], + output_data["temperature(H)"], + "c-", + label="H atoms", +) +ax.plot(output_data["time"], output_data["temperature"], "k-", label="All atoms") +ax.set_xlabel(r"$t$ / ps") +ax.set_ylabel(r"$\tilde{T}$ / K") +ax.legend() +plt.show() + + +# %% +# The velocity-velocity autocorrelation function is +# essentially indistinguishable from the reference, computed +# with an ensemble of NVE trajectories starting from canonical +# samples. In fact, the small discrepancies are mostly due to +# incomplete convergence of the averages in the short trajectory. +# +# This highlights the advantages of a global thermostat, that +# does not disrupt the natural diffusion in configuration space, +# and can often be used to compute dynamical, time-dependent +# observables out of a single trajectory -- which is far more +# practical than performing a collection of NVE trajectories. + +acf_svr = compute_acf_xyz( + "simulation_svr.vel_0.xyz", + maximum_lag=600, + length_zeropadding=2000, + spectral_windowing="cosine-blackman", + timestep=1, + time_units="femtosecond", + skip=100, +) + +fig, ax = plt.subplots(1, 1, figsize=(4, 3), constrained_layout=True) +ax.fill_between( + acf_ref[:1200, 0] * ha2cm1, + (acf_ref[:1200, 1] - acf_ref[:1200, 2]) * 1e5, + (acf_ref[:1200, 1] + acf_ref[:1200, 2]) * 1e5, + color="gray", + label="reference", +) +ax.loglog( + acf_svr[3][:1200] * ha2cm1, + acf_svr[4][:1200] * 1e5, + "b-", + label=r"SVR, $\tau=10$fs", +) +ax.set_xlabel(r"$\omega$ / cm$^{-1}$") +ax.set_ylabel(r"$\hat{c}_{vv}$ / arb. units") +ax.legend() +plt.show() + +# %% +# Generalized Langevin Equation thermostat +# ---------------------------------------- +# +# The issue with a Langevin thermostat is that, for a given coupling +# time :math:`\tau`, only molecular motions with a comparable time scale +# are sampled efficiently: faster modes are *underdamped*, and slower modes +# are *overdamped*, cf. the slowing down of diffusive behavior. +# +# A possible solution to this problem is using a +# *Generalized Langevin Equation (GLE)* thermostat. A GLE +# thermostat uses a matrix generalization of the Langevin term, +# in which the physical momentum is supplemented by a few fictitious +# momenta :math:`\mathbf{s}`, i.e. +# +# .. math:: +# +# (\dot{p},\dot\mathbf{s}) = -\mathbf{A}_p (p,\mathbf{s})+ +# \mathbf{B}_p (\xi,\boldsymbol{\xi}) +# +# Here :math:`\mathbf{A}_p` is the *drift matrix* and :math:`\mathbf{B}_p` +# is a diffusion matrix which, for canonical sampling, is determined by the target +# temperature and the drift matrix through a fluctuation-dissipation relation. +# The key idea is that :math:`\mathbf{A}_p` provides a lot of flexibility in defining +# the behavior of the GLE, that can be tuned to achieve near-optimal sampling +# for every degree of freedom (effectively acting as if the coupling constant was +# tuned separately for slow and fast molecular motions). +# The general idea and the practical implementation are discussed in +# `(Ceriotti et al. JCTC (2010)) `_ +# which also discusses other applications of the same principle, including +# performing simulations with a non-equilibrium *quantum thermostat* that +# mimics the quantum the quantum mechanical behavior of light nuclei. + +# %% +# Setting up a thermostat in ``i-PI`` +# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +# +# A GLE thermostat can be activated using ``mode="gle"``. +# The drift matrix used here has been generated from the +# `GLE4MD website `_, using parameters that +# aim for the most efficient sampling possible with the short +# simulation time (2 ps). The `online generator +# `_ +# can be tuned to provide the best possible sampling for the system +# of interest, the most important parameter being the slowest time scale +# that one is interested in sampling (typically a fraction of the total +# simulation time). The range of frequencies that is optimized can then +# be tuned so as to reach, roughly, the maximum frequency present in the +# system. + +xmlroot = ET.parse("data/input_gle.xml").getroot() +print(" " + ET.tostring(xmlroot.find(".//thermostat"), encoding="unicode")) + +# %% +# We launch ``i-PI`` as usual ... + +ipi_process = None +if not os.path.exists("simulation_gle.out"): + ipi_process = subprocess.Popen(["i-pi", "data/input_gle.xml"]) + time.sleep(4) # wait for i-PI to start + lmp_process = [subprocess.Popen(["lmp", "-in", "data/in.lmp"]) for i in range(1)] + +# %% +# ... and wait for simulations to finish. + +if ipi_process is not None: + ipi_process.wait() + lmp_process[0].wait() + +# %% +# Analysis of the trajectory +# ~~~~~~~~~~~~~~~~~~~~~~~~~~ +# +# The kinetic temperature equilibrates quickly to the target value. +# Since the GLE is a local thermostat, targeting each degree of freedom +# separately, equipartition is also reached quickly. Sampling is less +# fast than with an aggressive Langevin thermostat, because the GLE targets +# each vibrational frequency separately, to minimize the impact on diffusion. +output_data, output_desc = ipi.read_output("simulation_gle.out") +traj_data = ipi.read_trajectory("simulation_gle.pos_0.xyz") + +fig, ax = plt.subplots(1, 1, figsize=(4, 3), constrained_layout=True) +ax.plot( + output_data["time"], + output_data["temperature(O)"], + "r-", + label="O atoms", +) +ax.plot( + output_data["time"], + output_data["temperature(H)"], + "c-", + label="H atoms", +) +ax.plot(output_data["time"], output_data["temperature"], "k-", label="All atoms") +ax.set_xlabel(r"$t$ / ps") +ax.set_ylabel(r"$\tilde{T}$ / K") +ax.legend() +plt.show() + + +# %% +# :math:`\hat{c}_{vv}` reflects the adaptive behavior of the GLE. +# The fast modes are damped aggressively, leading to a large +# broadening of the high frequency peaks, but librations and diffusive +# modes are much less dampened than in the high-coupling Langevin case. +# An optimal-coupling GLE is a safe choice to sample any system, from +# molecular liquids to harmonic crystals, although a stochastic velocity +# rescaling is preferable if one is interested in preserving the natural +# dynamics. + +acf_gle = compute_acf_xyz( + "simulation_gle.vel_0.xyz", + maximum_lag=600, + length_zeropadding=2000, + spectral_windowing="cosine-blackman", + timestep=1, + time_units="femtosecond", + skip=100, +) + +fig, ax = plt.subplots(1, 1, figsize=(4, 3), constrained_layout=True) +ax.fill_between( + acf_ref[:1200, 0] * ha2cm1, + (acf_ref[:1200, 1] - acf_ref[:1200, 2]) * 1e5, + (acf_ref[:1200, 1] + acf_ref[:1200, 2]) * 1e5, + color="gray", + label="reference", +) +ax.loglog( + acf_gle[3][:1200] * ha2cm1, + acf_gle[4][:1200] * 1e5, + "b-", + label=r"GLE", +) +ax.set_xlabel(r"$\omega$ / cm$^{-1}$") +ax.set_ylabel(r"$\hat{c}_{vv}$ / arb. units") +ax.legend() +plt.show() + +# %% +# R-L purification +# ~~~~~~~~~~~~~~~~ +# +# What if you also want to extract dynamical information from a GLE +# (or Langevin) trajectory? It is actually possible to post-process the +# power spectrum, performing a deconvolution based on the amount of +# disturbance introduced by the GLE, that can be predicted analytically +# in the harmonic limit. +# The idea, discussed in `(Rossi et al., JCP (2018)) +# `_ +# is that if :math:`\hat{y}(\omega)` is the "natural" NVE power +# spectrum, and :math:`k_{\mathrm{GLE}}(\omega_0, \omega)` is the power +# spectrum predicted for a harmonic oscillator of frequency :math:`\omega_0`, +# then the spectrum from the GLE dynamics will be approximately +# +# .. math:: +# +# \hat{y}_{\mathrm{GLE}}(\omega) = \int \mathrm{d}\omega' +# k_{\mathrm{GLE}}(\omega', \omega) \hat{y}(\omega') +# +# The kernel can be computed analytically for all frequencies that +# are relevant for the power spectrum, based on the GLE parameters +# extracted from the input of ``i-PI``. + +n_omega = 1200 +Ap, Cp, Dp = get_gle_matrices("data/input_gle.xml") +gle_kernel = gle_frequency_kernel(acf_gle[3][:n_omega], Ap, Dp) + + +lomega = acf_gle[3][:n_omega] * ha2cm1 +fig, ax = plt.subplots(1, 1, figsize=(4, 3), constrained_layout=True) +levels = np.logspace(np.log10(gle_kernel.min()), np.log10(gle_kernel.max()), num=50) +contour = ax.contourf(lomega, lomega, gle_kernel, norm=mpl.colors.LogNorm()) +ax.set_xscale("log") +ax.set_yscale("log") +ax.set_xlabel(r"$\omega_0$ / cm$^{-1}$") +ax.set_ylabel(r"$\omega$ / cm$^{-1}$") +ax.set_xlim(10, 5000) +ax.set_ylim(10, 5000) +cbar = fig.colorbar(contour, ticks=[1e1, 1e3, 1e5, 1e7]) + +# %% +# The deconvolution is based on the Iterative Image Space +# Reconstruction Algorithm, which preserves the positive-definiteness +# of the spectrum + +isra_acf, history, errors, laplace = isra_deconvolute( + acf_gle[3][:n_omega], acf_gle[4][:n_omega], gle_kernel, 64, 4 +) + +# %% +# Even though the ISRA algorithm is less prone to enhancing noise than +# other deconvolution algorithms, successive iterations sharpen the spectrum +# but introduce higher and higher levles of noise, particularly on the +# low-frequency end of the spectrum so one has to choose when to stop. + +fig, ax = plt.subplots(1, 1, figsize=(4, 3), constrained_layout=True) +ax.loglog( + acf_gle[3][:1200] * ha2cm1, + acf_gle[4][:1200] * 1e5, + "b-", + label=r"GLE", +) +ax.loglog( + acf_gle[3][:1200] * ha2cm1, + history[0] * 1e5, + ":", + color="#4000D0", + label=r"iter[1]", +) +ax.loglog( + acf_gle[3][:1200] * ha2cm1, + history[2] * 1e5, + ":", + color="#A000A0", + label=r"iter[9]", +) +ax.loglog( + acf_gle[3][:1200] * ha2cm1, + history[4] * 1e5, + ":", + color="#D00040", + label=r"iter[17]", +) +ax.loglog( + acf_gle[3][:1200] * ha2cm1, + history[12] * 1e5, + ":", + color="#FF0000", + label=r"iter[49]", +) + +ax.set_xlabel(r"$\omega$ / cm$^{-1}$") +ax.set_ylabel(r"$\hat{c}_{vv}$ / arb. units") +ax.legend() +plt.show() + +# %% +# Especially in the high-frequency region, the deconvolution +# algorithm succees in recovering the underlying NVE dynamics, +# which can be useful whenever one wants to optimize statistical +# efficiency while still being able to estimate dynamical +# properties. + +fig, ax = plt.subplots(1, 1, figsize=(4, 3), constrained_layout=True) +ax.fill_between( + acf_ref[:1200, 0] * ha2cm1, + (acf_ref[:1200, 1] - acf_ref[:1200, 2]) * 1e5, + (acf_ref[:1200, 1] + acf_ref[:1200, 2]) * 1e5, + color="gray", + label="reference", +) +ax.loglog( + acf_gle[3][:1200] * ha2cm1, + acf_gle[4][:1200] * 1e5, + "b-", + label=r"GLE", +) +ax.loglog( + acf_gle[3][:1200] * ha2cm1, + history[2] * 1e5, + "r-", + label=r"GLE$\rightarrow$ NVE (iter[5])", +) +ax.set_xlabel(r"$\omega$ / cm$^{-1}$") +ax.set_ylabel(r"$\hat{c}_{vv}$ / arb. units") +ax.legend() +plt.show() + + +# %% +# Running with LAMMPS +# ~~~~~~~~~~~~~~~~~~~ +# +# GLE thermostats (as well as conventional Langevin, and +# stochastic velocity rescaling) are also implemented natively +# in ``LAMMPS``. +# +# An example of ``LAMMPS`` input containing a GLE thermostat can +# be found in ``data/gle.lmp``. See also the +# `documentation of the fix gle command +# `_ +# +# .. code-block:: text +# +# fix 1 all gle 6 300 300 31415 data/smart.A +# +# The drift matrix can be obtained from the same website, simply +# asking to output the matrix in raw format, choosing units consistent +# with the ``LAMMPS`` settings, e.g. for this `optimal sampling setup +# `_ +# + +# %% +# We can run ``LAMMPS`` from the command line +# +# .. code-block:: bash +# +# lmp -in data/gle.lmp & +# +# or from Python + +lmp_process = None +if not os.path.exists("lammps_out.dat"): + lmp_process = subprocess.Popen(["lmp", "-in", "data/gle.lmp"]) + +# %% +# ... and wait +# + +if lmp_process is not None: + lmp_process.wait() + +# %% +# The simulation is much faster (for such a small system and +# cheap potential the overhead of ``i-PI``'s client-server mechanism +# is substantial) and leads to similar results for the kinetic temperature +# + +traj_data = np.loadtxt("lammps_out.dat") + +fig, ax = plt.subplots(1, 1, figsize=(4, 3), constrained_layout=True) +ax.plot( + traj_data[:, 0] * 1e-3, + traj_data[:, 1], + "k-", + label="All atoms", +) +ax.set_xlabel(r"$t$ / ps") +ax.set_ylabel(r"$\tilde{T}$ / K") +ax.legend() +plt.show()