diff --git a/docs/src/sampling.rst b/docs/src/sampling.rst index b09b7544..dbb14fe1 100644 --- a/docs/src/sampling.rst +++ b/docs/src/sampling.rst @@ -12,3 +12,4 @@ set of configurations of an atomistic system. examples/path-integrals/path-integrals examples/pi-metad/pi-metad examples/batch-cp2k/reference-trajectory + examples/heat-capacity/heat-capacity diff --git a/examples/heat-capacity/README.rst b/examples/heat-capacity/README.rst new file mode 100644 index 00000000..2e12bbf7 --- /dev/null +++ b/examples/heat-capacity/README.rst @@ -0,0 +1,7 @@ +Heat capacity of water from path integral molecular dynamics +============================================================ + +This example shows how to calculate the constant-volume heat +capacity from a path integral molecular dynamics simulation +using ``i-PI``. `LAMMPS `_ is used +as the driver for the evaluation of the q-TIP4P/f water model. diff --git a/examples/heat-capacity/data/in.lmp b/examples/heat-capacity/data/in.lmp new file mode 100644 index 00000000..61e74927 --- /dev/null +++ b/examples/heat-capacity/data/in.lmp @@ -0,0 +1,34 @@ +units electron +atom_style full + +pair_style lj/cut/tip4p/long 1 2 1 1 0.278072379 17.001 +# high-pppm precision and shift to get meaningful fd estimates +kspace_style pppm/tip4p 1e-5 +pair_modify shift yes +bond_style class2 +angle_style harmonic + + +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/heat-capacity/data/input.xml b/examples/heat-capacity/data/input.xml new file mode 100644 index 00000000..ae6b8209 --- /dev/null +++ b/examples/heat-capacity/data/input.xml @@ -0,0 +1,35 @@ + + + + [ step, time{picosecond}, conserved, potential, kinetic_cv, + scaledcoords(fd_delta=5e-3) ] + + + 2000 + + 32342 + + +
h2o-lammps
1e-4 + + + + data/water_32.pdb + 298 + + + lmpserial + + + 298 + + + + + 5.0 + + 0.5 + + + + diff --git a/examples/heat-capacity/data/water_32.pdb b/examples/heat-capacity/data/water_32.pdb new file mode 100644 index 00000000..014e40ff --- /dev/null +++ b/examples/heat-capacity/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/heat-capacity/data/water_32_data.lmp b/examples/heat-capacity/data/water_32_data.lmp new file mode 100644 index 00000000..723d35f9 --- /dev/null +++ b/examples/heat-capacity/data/water_32_data.lmp @@ -0,0 +1,324 @@ +LAMMPS Description + + 96 atoms + 64 bonds + 32 angles + + 2 atom types + 1 bond types + 1 angle types + + -19 19.38859 xlo xhi + -10 19.38859 ylo yhi + -19 19.38859 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/heat-capacity/environment.yml b/examples/heat-capacity/environment.yml new file mode 100644 index 00000000..8c91ad59 --- /dev/null +++ b/examples/heat-capacity/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@heat-capacity-clean + - matplotlib + - skmatter diff --git a/examples/heat-capacity/heat-capacity.py b/examples/heat-capacity/heat-capacity.py new file mode 100644 index 00000000..e0174ea2 --- /dev/null +++ b/examples/heat-capacity/heat-capacity.py @@ -0,0 +1,296 @@ +""" +Quantum heat capacity of water +============================== + +:Authors: +Filippo Bigi `@frostedoyster `_; +Michele Ceriotti `@ceriottm `_ + +This example shows how to estimate the heat capacity of liquid water +from a path integral molecular dynamics simulation. The dynamics are +run with `i-PI `_, and +`LAMMPS `_ is used +as the driver to simulate the `q-TIP4P/f water +model `_. +""" + +import os +import subprocess +import time +import xml.etree.ElementTree as ET + +import ipi +import matplotlib.pyplot as plt +import numpy as np + + +# %% +# A non-trivial energy estimator +# ------------------------------ +# +# As introduced in the +# `path-integrals example +# `_, +# path-integral estimators +# for observables that depend on momenta are generally not trivial to compute. +# +# In this example, we will focus on the constant-volume heat capacity, +# :math:`c_V`, which is one such +# observable, and we will calculate it for liquid water at room temperature. +# Because of the presence of high-frequency vibrations, many of the nuclear +# degrees of freedom are trapped in the vibrational ground state, which reduces +# substantially the heat capacity from the value that would be obtained +# in the classical limit. See `this review +# `_ +# for an overview of the impact of quantum nuclei on the properties of water. +# From a computational perspective, this means that it is necessary to use +# specialized simulations and estimators to evaluate the correct value. + + +# %% +# Running the PIMD calculation +# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +# +# This follows the same steps as the ``path-integrals`` example. One important +# difference is that we will request the ``scaledcoords`` output to the relevant +# section of the ``i-PI`` input XML file, which +# contains estimators that can be used to calculate the total energy and +# heat capacity as following +# `Yamamoto, J. Chem. Phys. (2005) `_. +# +# The input file is shown below. It should be noted that ``scaledcoords`` +# is given a finite differences displacement as a parameter. This is necessary +# as the estimators require higher order derivatives of the potential energy, +# which are calculated using finite differences. This also means that +# evaluating the estimator adds substantial overhead (so it is wise to only +# compute it every few simulation steps, to eliminate correlations between +# snapshots) and that one should be careful to use well-converged simulation +# parameters to avoid discontinuities and noise (for instance, we increase +# the accuracy of the particle-mesh electrostatic calculation, and use a +# shifted Lennard-Jones potential to avoid a discontinuity at the cutoff). + +# Open and show the relevant part of the input +xmlroot = ET.parse("data/input.xml").getroot() +print(" " + ET.tostring(xmlroot.find(".//properties"), encoding="unicode")) + +# %% +# We launch the i-PI and LAMMPS processes, exactly as in the +# ``path-integrals`` example. + +# don't rerun if the outputs already exist +ipi_process = None +if not os.path.exists("water-cv.out"): + ipi_process = subprocess.Popen(["i-pi", "data/input.xml"]) + time.sleep(2) # wait for i-PI to start + lmp_process = [subprocess.Popen(["lmp", "-in", "data/in.lmp"]) for i in range(2)] + +# %% +# Skip this cell if you want to run in the background +if ipi_process is not None: + ipi_process.wait() + lmp_process[0].wait() + lmp_process[1].wait() + +# %% +# Analyzing the results +# ~~~~~~~~~~~~~~~~~~~~~ +# Let's plot the potential and conserved energy as a function of time, +# just to check that the simulation ran sensibly. + +output_data, output_desc = ipi.read_output("water-cv.out") +fix, 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["conserved"] - output_data["conserved"][0], + "r-", + label="Conserved, $H$", +) +ax.set_xlabel(r"$t$ / ps") +ax.set_ylabel(r"energy / a.u.") +ax.legend() +plt.show() + +# %% +# As described in the `i-PI documentation +# `_, +# the two quantities returned by the ``scaledcoords`` output are ``eps_v`` +# and ``eps_v'``, defined in the aforementioned +# `paper `_. +# +# These estimators (:math:`\epsilon_v` and :math:`\epsilon_v'`) are derived in the +# "scaled coordinates" formalism, which is a useful trick to avoid the +# growth of the error in the instantaneous values of the estimators with +# the number of beads used in the path integral simulation. +# +# The same paper contains the formulas to calculate the total energy and +# heat capacity from these estimators: +# +# .. math:: +# E = \langle \epsilon_v \rangle \quad +# C_V = k_B \beta^2 \left( \langle \epsilon_v^2 \rangle - \langle +# \epsilon_v \rangle^2 - \langle \epsilon_v' \rangle \right) +# +# First the energy, whose estimator will be compared to the total energy +# calculated as the sum of the potential and kinetic energy estimators. +# Since the kinetic energy is itself calculated from a scaled-coordinates +# estimator (the "centroid virial" estimator), the two total energies are +# the same. + +eps_v = output_data["scaledcoords(fd_delta=5e-3)"][:, 0] +eps_v_prime = output_data["scaledcoords(fd_delta=5e-3)"][:, 1] + +energy_estimator = eps_v # first formula + +fix, ax = plt.subplots(1, 1, figsize=(4, 3), constrained_layout=True) +ax.plot( + output_data["time"], + energy_estimator - energy_estimator[0], + "b-", + label="scaled coordinates estimator", +) +ax.plot( + output_data["time"][:], + ( + output_data["potential"] + - output_data["potential"][0] + + output_data["kinetic_cv"] + - output_data["kinetic_cv"][0] + ), + "r.", + label="potential + virial kinetic", +) +ax.set_xlabel(r"$t$ / ps") +ax.set_ylabel(r"total energy / a.u.") +ax.legend() +plt.show() + + +# %% +# And, finally, the heat capacity. Note that normally the simulation +# requires a few ps for equilibration. Here we discard a few dozen steps +# to eliminate the initial jump, which is due to the relaxation of the +# ring polymers starting from a single atomic configuration. + +# i-PI scaledcoords outputs are in atomic units (see docs) +kB = 3.16681e-6 # Boltzmann constant in atomic units +T = 298.0 # temperature in K, as defined in the input file +beta = 1.0 / (kB * T) + +skip = 20 +heat_capacity = ( # second formula + kB + * (beta**2) + * ( + np.mean(eps_v[skip:] ** 2) + - np.mean(eps_v[skip:]) ** 2 + - np.mean(eps_v_prime[skip:]) + ) +) +heat_capacity_per_molecule = heat_capacity / 32 # 32 molecules in the simulation +print(f"Heat capacity (per water molecule): {(heat_capacity_per_molecule/kB):.2f} kB") + +# %% +# You may recognize that the first part of the estimator is reminiscent +# of the classical estimator for the heat capacity as the fluctuations of the +# (quantum) total energy, which in this case however requires a correction given +# by the mean of the second part of the scaled-coordinates estimator. + +# %% +# Estimating the statistical error +# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +# +# Especially with such an underconverged simulation, it is important to +# estimate the statistical error in the heat capacity. +# +# Generally, errors on measurements are computed +# as "standard errors", i.e. the standard deviation of a series of data points +# divided by the square root of the number of data points. In our case, +# however, this is made more complicated by the correlation between +# close steps in the molecular dynamics trajectory, which would lead to an +# overestimation of the number of independent samples. To fix this, we can +# calculate the autocorrelation time of the estimators whose errors we +# want to estimate, and apply a correction factor to the number of samples. + + +def autocorrelate(x): + n = len(x) + xo = x - x.mean() # remove mean + acov = (np.correlate(xo, xo, "full"))[n - 1 :] + return acov[: len(acov) // 2] + + +def autocorrelation_time(x): + acov = autocorrelate(x) + return 1.0 + np.sum(acov) / acov[0] + + +# %% +# Using these helper functions, we can now calculate the error on the various +# parts of the heat capacity estimator. Note also the autocorrelation times, that +# are just a little larger than one, indicating that the stride used to print out +# the estimators is appropriate (as there is little correlation between the samples). + +# Autocorrelation times (i.e. number of steps needed to have independent samples) +autocorr_time_error_delta_eps_v = autocorrelation_time( + (eps_v[skip:] - eps_v[skip:].mean()) ** 2 +) +autocorr_time_error_eps_v_prime = autocorrelation_time(eps_v_prime[skip:]) + +print( + f""" +Autocorrelation times (in number of samples): +(eps-)^2: {autocorr_time_error_delta_eps_v:.2f} +eps': {autocorr_time_error_eps_v_prime:.2f} +""" +) + +# Effective number of samples +effective_samples_delta_eps_v = len(eps_v[skip:]) / autocorr_time_error_delta_eps_v +effective_samples_eps_v_prime = len(eps_v[skip:]) / autocorr_time_error_eps_v_prime + +# Standard errors using the effective number of samples +error_delta_eps_v = np.std((eps_v[skip:] - eps_v[skip:].mean()) ** 2) / np.sqrt( + effective_samples_delta_eps_v +) +error_eps_v_prime = np.std(eps_v_prime[skip:]) / np.sqrt(effective_samples_eps_v_prime) + +# Error on the heat capacity (assuming quadrature sum) +error_heat_capacity = ( + kB * (beta**2) * np.sqrt(error_delta_eps_v**2 + error_eps_v_prime**2) +) + +error_heat_capacity_per_molecule = ( + error_heat_capacity / 32 +) # 32 molecules in the simulation + +print( + "Error on the heat capacity (per water molecule): " + f"{(error_heat_capacity_per_molecule/kB):.2f} kB" +) + + +# %% +# The obtained heat capacity is consistent with the values from the literature +# (see e.g. `Ceriotti et al., J. Chem. Phys. (2011) +# `_ +# where the convergence of the heat capacity with number of beads is shown for the same +# water model used in this example). +# However, the error is quite large, which is expected given the short simulation time. +# To reduce the error, one would need to run a longer simulation. Other important error +# sources, which are not accounted for in the error estimate, are the finite size of the +# system and number of beads. Both of these are too small in this example to give +# reliable results. +# +# In a realistic simulation, up to a few 100s of picoseconds might be needed to reduce +# the sampling error to a small value (1-2% of the heat capacity). For water at room +# temperature, you will need 32 beads at the very least (8 were used in this example). +# It is more difficult to give a general rule for the system size: (quantum) energy +# fluctuations are usually localized, but to guarantee accurate sampling of the +# liquid structure, a few hundred water molecules would be a reasonable guess +# (32 were used in this example). diff --git a/examples/pi-metad/pi-metad.py b/examples/pi-metad/pi-metad.py index 1dc87f51..52a17382 100644 --- a/examples/pi-metad/pi-metad.py +++ b/examples/pi-metad/pi-metad.py @@ -467,7 +467,7 @@ def moving_average(arr, window_size): # ----------------------------------- # # You can see `this recipe -# `_ +# `_ # for a brief introduction to path integral simulations with `i-PI`. # From a practical perspective, very little needs to change with respect # to the classical case.