From 0efa65dff20599051e5d890f3e87c7e704b77f4e Mon Sep 17 00:00:00 2001 From: Michele Ceriotti Date: Thu, 26 Sep 2024 16:17:58 +0200 Subject: [PATCH 1/7] Stub of an o3-averaging example --- examples/o3-averaging/README.rst | 6 + examples/o3-averaging/data/in.lmp | 34 ++ examples/o3-averaging/data/input-noo3.xml | 41 +++ examples/o3-averaging/data/water_32.pdb | 98 ++++++ examples/o3-averaging/data/water_32_data.lmp | 324 +++++++++++++++++++ examples/o3-averaging/environment.yml | 14 + examples/o3-averaging/o3-averaging.py | 90 ++++++ 7 files changed, 607 insertions(+) create mode 100644 examples/o3-averaging/README.rst create mode 100644 examples/o3-averaging/data/in.lmp create mode 100644 examples/o3-averaging/data/input-noo3.xml create mode 100644 examples/o3-averaging/data/water_32.pdb create mode 100644 examples/o3-averaging/data/water_32_data.lmp create mode 100644 examples/o3-averaging/environment.yml create mode 100644 examples/o3-averaging/o3-averaging.py diff --git a/examples/o3-averaging/README.rst b/examples/o3-averaging/README.rst new file mode 100644 index 00000000..649f589f --- /dev/null +++ b/examples/o3-averaging/README.rst @@ -0,0 +1,6 @@ +Rotational averaging of non-equivariant models +============================================== + +This example shows how to assess (and to an extent correct) the +influence of non-equivariant terms in the potential energy computed +by an unconstrained machine-learning model. diff --git a/examples/o3-averaging/data/in.lmp b/examples/o3-averaging/data/in.lmp new file mode 100644 index 00000000..61e74927 --- /dev/null +++ b/examples/o3-averaging/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/o3-averaging/data/input-noo3.xml b/examples/o3-averaging/data/input-noo3.xml new file mode 100644 index 00000000..0a4da305 --- /dev/null +++ b/examples/o3-averaging/data/input-noo3.xml @@ -0,0 +1,41 @@ + + + + [ step, time{picosecond}, conserved, potential, kinetic_cv, + scaledcoords(fd_delta=5e-3) ] + + extras_component_raw(1) + positions + + 2000 + + 32342 + + +
h2o-lammps
1e-4 + + +
h2o-noo3
1e-4 + + + + data/water_32.pdb + 300 + + + + + + + 300 + + + + + 100.0 + + 0.5 + + + + diff --git a/examples/o3-averaging/data/water_32.pdb b/examples/o3-averaging/data/water_32.pdb new file mode 100644 index 00000000..014e40ff --- /dev/null +++ b/examples/o3-averaging/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/o3-averaging/data/water_32_data.lmp b/examples/o3-averaging/data/water_32_data.lmp new file mode 100644 index 00000000..723d35f9 --- /dev/null +++ b/examples/o3-averaging/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/o3-averaging/environment.yml b/examples/o3-averaging/environment.yml new file mode 100644 index 00000000..01315f9a --- /dev/null +++ b/examples/o3-averaging/environment.yml @@ -0,0 +1,14 @@ +channels: + - conda-forge +dependencies: + - python=3.11 + - pip + - gfortran + - make + - git + - lammps + - pip: + - ase==3.22.1 + - chemiscope + - git+https://github.com/lab-cosmo/i-pi.git@feat/equivator#egg=ipi + - matplotlib diff --git a/examples/o3-averaging/o3-averaging.py b/examples/o3-averaging/o3-averaging.py new file mode 100644 index 00000000..c19ad351 --- /dev/null +++ b/examples/o3-averaging/o3-averaging.py @@ -0,0 +1,90 @@ +r""" +Rotational averaging of non-equivariant models +========================== + +:Authors: +Marcel Langer `@sirmarcel `_; +Michele Ceriotti `@ceriottm `_ + +DESCRIPTION +""" + +# %% + +import bz2 +import os +import subprocess +import time +import xml.etree.ElementTree as ET + +import ase +import ase.io +import chemiscope +import ipi +import matplotlib.pyplot as plt +import numpy as np + + +# %% +# Show the problem +# ---------------------------------- + + +# %% +# + +# Open and show the relevant part of the input +xmlroot = ET.parse("data/input-noo3.xml").getroot() +print(" " + ET.tostring(xmlroot.find(".//ffrotations"), encoding="unicode")) + + +# %% +# Installing the Python driver +# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +# +# i-PI comes with a FORTRAN driver, which however has to be installed +# from source. We use a utility function to compile it. Note that this requires +# a functioning build system with `gfortran` and `make`. + +ipi.install_driver() + +# %% +# 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-noo3.out"): + ipi_process = subprocess.Popen(["i-pi", "data/input-noo3.xml"]) + time.sleep(2) # wait for i-PI to start + lmp_process = [subprocess.Popen(["lmp", "-in", "data/in.lmp"]) for i in range(1)] + driver_process = [ + subprocess.Popen( + ["i-pi-driver", "-u", "-a", "h2o-noo3", "-m", "noo3-h2o"], cwd="data/" + ) + for i in range(1) + ] + +# %% +# Skip this cell if you want to run in the background +if ipi_process is not None: + ipi_process.wait() + lmp_process[0].wait() + driver_process[0].wait() + + +# %% +# +# Discuss how molecules get oriented + +traj_data = ipi.read_trajectory( + "water-noo3.pos_0.extxyz", format="ase", dimension="length" +) + +# %% +# then, assemble a visualization +chemiscope.show( + frames=traj_data, + settings=chemiscope.quick_settings(trajectory=True), + mode="structure", +) From bd43c6d1da3fdda1768f1898654b26869ed88831 Mon Sep 17 00:00:00 2001 From: Michele Ceriotti Date: Thu, 26 Sep 2024 16:24:14 +0200 Subject: [PATCH 2/7] Linting and ipi.read -> ase.read --- examples/o3-averaging/o3-averaging.py | 17 ++++++++--------- 1 file changed, 8 insertions(+), 9 deletions(-) diff --git a/examples/o3-averaging/o3-averaging.py b/examples/o3-averaging/o3-averaging.py index c19ad351..d7eb9226 100644 --- a/examples/o3-averaging/o3-averaging.py +++ b/examples/o3-averaging/o3-averaging.py @@ -2,16 +2,15 @@ Rotational averaging of non-equivariant models ========================== -:Authors: -Marcel Langer `@sirmarcel `_; -Michele Ceriotti `@ceriottm `_ +:Authors: +Marcel Langer `@sirmarcel `_; +Michele Ceriotti `@ceriottm `_ DESCRIPTION """ # %% -import bz2 import os import subprocess import time @@ -21,8 +20,10 @@ import ase.io import chemiscope import ipi -import matplotlib.pyplot as plt -import numpy as np + + +# import matplotlib.pyplot as plt +# import numpy as np # %% @@ -77,9 +78,7 @@ # # Discuss how molecules get oriented -traj_data = ipi.read_trajectory( - "water-noo3.pos_0.extxyz", format="ase", dimension="length" -) +traj_data = ase.io.read("water-noo3.pos_0.extxyz", ":") # %% # then, assemble a visualization From 4a18034a70454de257bfa4a84ec9f7e4bc9d7159 Mon Sep 17 00:00:00 2001 From: Michele Ceriotti Date: Sun, 29 Sep 2024 14:55:01 +0200 Subject: [PATCH 3/7] Added the o3-averaging example to relevant sections --- docs/src/software/chemiscope.sec | 1 + docs/src/software/i-pi.sec | 1 + docs/src/software/lammps.sec | 1 + docs/src/topics/ml-models.sec | 1 + docs/src/topics/sampling.sec | 1 + 5 files changed, 5 insertions(+) diff --git a/docs/src/software/chemiscope.sec b/docs/src/software/chemiscope.sec index 8ded4143..811f2b9a 100644 --- a/docs/src/software/chemiscope.sec +++ b/docs/src/software/chemiscope.sec @@ -17,3 +17,4 @@ repository `_. - examples/pi-metad/pi-metad - examples/path-integrals/path-integrals - examples/thermostats/thermostats +- examples/o3-averaging/o3-averaging diff --git a/docs/src/software/i-pi.sec b/docs/src/software/i-pi.sec index 10c69f63..a6ce41ee 100644 --- a/docs/src/software/i-pi.sec +++ b/docs/src/software/i-pi.sec @@ -12,3 +12,4 @@ it on the `ipi-code website `_, the `documentation pages - examples/path-integrals/path-integrals - examples/pi-metad/pi-metad - examples/heat-capacity/heat-capacity +- examples/o3-averaging/o3-averaging diff --git a/docs/src/software/lammps.sec b/docs/src/software/lammps.sec index b937977e..9ba50fed 100644 --- a/docs/src/software/lammps.sec +++ b/docs/src/software/lammps.sec @@ -10,3 +10,4 @@ efficient domain decomposition scheme. Learn more about LAMMPS on its `homepage - examples/thermostats/thermostats - examples/path-integrals/path-integrals - examples/heat-capacity/heat-capacity +- examples/o3-averaging/o3-averaging diff --git a/docs/src/topics/ml-models.sec b/docs/src/topics/ml-models.sec index 1e741685..4e70d36d 100644 --- a/docs/src/topics/ml-models.sec +++ b/docs/src/topics/ml-models.sec @@ -9,3 +9,4 @@ data. - examples/lode-linear/lode-linear - examples/sample-selection/sample-selection - examples/periodic-hamiltonian/periodic-hamiltonian +- examples/o3-averaging/o3-averaging diff --git a/docs/src/topics/sampling.sec b/docs/src/topics/sampling.sec index 8c20d470..e8ec9737 100644 --- a/docs/src/topics/sampling.sec +++ b/docs/src/topics/sampling.sec @@ -10,3 +10,4 @@ set of configurations of an atomistic system. - examples/pi-metad/pi-metad - examples/batch-cp2k/reference-trajectory - examples/heat-capacity/heat-capacity +- examples/o3-averaging/o3-averaging From 8a133e03d4cfc607003c3a2a6e80d526a4f3ca88 Mon Sep 17 00:00:00 2001 From: Michele Ceriotti Date: Thu, 26 Sep 2024 16:17:58 +0200 Subject: [PATCH 4/7] Stub of an o3-averaging example --- examples/o3-averaging/o3-averaging.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/examples/o3-averaging/o3-averaging.py b/examples/o3-averaging/o3-averaging.py index d7eb9226..a3b30c9f 100644 --- a/examples/o3-averaging/o3-averaging.py +++ b/examples/o3-averaging/o3-averaging.py @@ -11,6 +11,7 @@ # %% +import bz2 import os import subprocess import time @@ -24,6 +25,8 @@ # import matplotlib.pyplot as plt # import numpy as np +import matplotlib.pyplot as plt +import numpy as np # %% From b1a6eee2bb48dc9cc851dd5b063701510d62e99c Mon Sep 17 00:00:00 2001 From: Michele Ceriotti Date: Thu, 26 Sep 2024 16:24:14 +0200 Subject: [PATCH 5/7] Linting and ipi.read -> ase.read --- examples/o3-averaging/o3-averaging.py | 3 --- 1 file changed, 3 deletions(-) diff --git a/examples/o3-averaging/o3-averaging.py b/examples/o3-averaging/o3-averaging.py index a3b30c9f..d7eb9226 100644 --- a/examples/o3-averaging/o3-averaging.py +++ b/examples/o3-averaging/o3-averaging.py @@ -11,7 +11,6 @@ # %% -import bz2 import os import subprocess import time @@ -25,8 +24,6 @@ # import matplotlib.pyplot as plt # import numpy as np -import matplotlib.pyplot as plt -import numpy as np # %% From 9e1b48fa5ffa42bc31da86fffd0455750040e328 Mon Sep 17 00:00:00 2001 From: Michele Ceriotti Date: Sat, 26 Oct 2024 16:20:20 +0200 Subject: [PATCH 6/7] Added a visualization of what the rotation averaging does --- examples/o3-averaging/o3-averaging.py | 62 +++++++++++++++++++++++++++ 1 file changed, 62 insertions(+) diff --git a/examples/o3-averaging/o3-averaging.py b/examples/o3-averaging/o3-averaging.py index d7eb9226..eb70293b 100644 --- a/examples/o3-averaging/o3-averaging.py +++ b/examples/o3-averaging/o3-averaging.py @@ -20,11 +20,73 @@ import ase.io import chemiscope import ipi +import numpy as np +from ipi.utils.mathtools import ( + get_rotation_quadrature_lebedev, + get_rotation_quadrature_legendre, +) # import matplotlib.pyplot as plt # import numpy as np +# %% +# show rotations + +# Gets quadrature grid (regular spacing in alpha, +# lebedev grids for beta, gamma) +quad_grid = get_rotation_quadrature_lebedev(3) +quad = { + "rotation_matrices": np.array([q[0] for q in quad_grid]), + "weights": np.array([q[1] for q in quad_grid]), + "angles": np.array([q[2] for q in quad_grid]), +} + +# Display the rotations of an atom with a "force" vector +# The "reference atom" is marked as an O, and its copies as H +vs = [ + np.asarray([[2, 0, 0], [2, 1, 0]]), + np.asarray([[0, 2, 0], [0, 2, 1]]), + np.asarray([[0, 0, 2], [1, 0, 2]]), +] + +frames = [] +for v in vs: + rl = v @ quad["rotation_matrices"] + ats = ase.Atoms("O" + ("H" * (len(rl) - 1)), positions=rl[:, 0]) + ats.arrays["forces"] = rl[:, 1] - rl[:, 0] + ats.arrays["weight"] = quad["weights"] + ats.arrays["alpha"] = quad["angles"][:, 0] + ats.arrays["beta"] = quad["angles"][:, 1] + ats.arrays["gamma"] = quad["angles"][:, 2] + frames.append(ats) + +# Display with chemiscope. Three frames for three +# initial positions. You can also color by grid weight +# and by Euler angles + +chemiscope.show( + frames=frames, + mode="structure", + shapes={"forces": chemiscope.ase_vectors_to_arrows(frames)}, + properties=chemiscope.extract_properties(frames), + settings={ + "structure": [ + dict( + bonds=False, + atoms=True, + shape="forces", + environments={"activated": False}, + color={ + "property": "element", + "transform": "linear", + "palette": "viridis", + }, + ), + ] + }, + environments=chemiscope.all_atomic_environments(frames), +) # %% # Show the problem From 2040804f3a38b0c3ee75a444f2ebe00bca3e789a Mon Sep 17 00:00:00 2001 From: Marcel Langer Date: Fri, 1 Nov 2024 14:15:01 +0100 Subject: [PATCH 7/7] Make some progress --- examples/o3-averaging/data/in-mol.lmp | 35 ++++ examples/o3-averaging/data/init-mol.xyz | 5 + examples/o3-averaging/data/input-mol-grid.xml | 38 ++++ examples/o3-averaging/data/input-mol-noo3.xml | 35 ++++ examples/o3-averaging/data/replay-noo3.xml | 26 +++ examples/o3-averaging/data/water_1_data.lmp | 41 ++++ examples/o3-averaging/o3-averaging.py | 177 +++++++++++++----- 7 files changed, 311 insertions(+), 46 deletions(-) create mode 100644 examples/o3-averaging/data/in-mol.lmp create mode 100644 examples/o3-averaging/data/init-mol.xyz create mode 100644 examples/o3-averaging/data/input-mol-grid.xml create mode 100644 examples/o3-averaging/data/input-mol-noo3.xml create mode 100644 examples/o3-averaging/data/replay-noo3.xml create mode 100644 examples/o3-averaging/data/water_1_data.lmp diff --git a/examples/o3-averaging/data/in-mol.lmp b/examples/o3-averaging/data/in-mol.lmp new file mode 100644 index 00000000..5d620f00 --- /dev/null +++ b/examples/o3-averaging/data/in-mol.lmp @@ -0,0 +1,35 @@ +units electron +atom_style full + +pair_style lj/cut/tip4p/cut 1 2 1 1 0.278072379 17.001 +# high-pppm precision and shift to get meaningful fd estimates +# kspace_style none +# kspace_style pppm/tip4p 1e-5 +pair_modify shift yes +bond_style class2 +angle_style harmonic + + +read_data water_1_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/o3-averaging/data/init-mol.xyz b/examples/o3-averaging/data/init-mol.xyz new file mode 100644 index 00000000..18c674e0 --- /dev/null +++ b/examples/o3-averaging/data/init-mol.xyz @@ -0,0 +1,5 @@ +3 +Lattice="100.0 0.0 0.0 0.0 100.0 0.0 0.0 0.0 100.0" Properties=species:S:1:pos:R:3:momenta:R:3 ipi_comment="Step: 200000 Bead: 0 positions{angstrom} cell{angstrom}" pbc="T T T" +O -2.89770838 6.76492149 25.78733892 0.41094839 -0.09429088 0.29861888 +H -3.22693259 6.92327461 24.89600383 -0.31547195 0.04597911 -0.02131858 +H -3.73951469 6.82619967 26.21403556 -0.09547645 0.04831177 -0.27730030 diff --git a/examples/o3-averaging/data/input-mol-grid.xml b/examples/o3-averaging/data/input-mol-grid.xml new file mode 100644 index 00000000..0d9cd5e5 --- /dev/null +++ b/examples/o3-averaging/data/input-mol-grid.xml @@ -0,0 +1,38 @@ + + + +
h2o-lammps
1e-4 + + +
h2o-noo3
+ 3 + legendre + True + + 10000 + + positions + [ step, time, conserved, temperature{kelvin}, kinetic_md, potential, pressure_md ] + + + 32123 + + + + + + + + init-mol.xyz + + + 300 + + + True + + 0.5 + + + + diff --git a/examples/o3-averaging/data/input-mol-noo3.xml b/examples/o3-averaging/data/input-mol-noo3.xml new file mode 100644 index 00000000..4f452f96 --- /dev/null +++ b/examples/o3-averaging/data/input-mol-noo3.xml @@ -0,0 +1,35 @@ + + + +
h2o-lammps
1e-4 + + +
h2o-noo3
+ + 10000 + + positions + [ step, time, conserved, temperature{kelvin}, kinetic_md, potential, pressure_md ] + + + 32123 + + + + + + + + init-mol.xyz + + + 300 + + + True + + 0.5 + + + + diff --git a/examples/o3-averaging/data/replay-noo3.xml b/examples/o3-averaging/data/replay-noo3.xml new file mode 100644 index 00000000..549c0c71 --- /dev/null +++ b/examples/o3-averaging/data/replay-noo3.xml @@ -0,0 +1,26 @@ + + + + [ step, time{picosecond}, conserved, potential, kinetic_cv, + scaledcoords(fd_delta=5e-3) ] + + + +
h2o-lammps
1e-4 + + +
h2o-noo3
1e-4 + + + + init-mol.xyz + + + + + + + replay.xyz + + + diff --git a/examples/o3-averaging/data/water_1_data.lmp b/examples/o3-averaging/data/water_1_data.lmp new file mode 100644 index 00000000..ccc1220f --- /dev/null +++ b/examples/o3-averaging/data/water_1_data.lmp @@ -0,0 +1,41 @@ +LAMMPS Description + + 3 atoms + 2 bonds + 1 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 + +Bonds + + 1 1 1 2 + 2 1 1 3 + +Angles + + 1 1 2 1 3 diff --git a/examples/o3-averaging/o3-averaging.py b/examples/o3-averaging/o3-averaging.py index eb70293b..511bd2fc 100644 --- a/examples/o3-averaging/o3-averaging.py +++ b/examples/o3-averaging/o3-averaging.py @@ -1,15 +1,36 @@ r""" Rotational averaging of non-equivariant models -========================== +============================================== :Authors: Marcel Langer `@sirmarcel `_; Michele Ceriotti `@ceriottm `_ -DESCRIPTION +Molecular dynamics (MD) simulates the movement of atoms on the potential energy surface +(PES) :math:`U(R)` where :math:`R` stands for all :math:`i` atomic positions. Evaluating +this function makes up the dominating cost in running MD. Machine learning interatomic +potentials (MLPs) are often used to approximate some expensive first-principles PES at +high accuracy and reasonably high speed. + +It is well known that the PES has underlying physical symmetries: It doesn't change under +rotations of the coordinates, overall translations in space, and it doesn't depend on the +order in which we list all the atoms in the system. It has long been believed that it is +essential for MLPs to exactly fulfil these symmetries. However, this places severe +constraints on possible model architectures, and so there has been a lot of interest in +so-called *unconstrained* models that do not respect symmetries exactly, by construction, +but rather learn them approximately from the training data. In the domain of MLPs, such +models have largely focused on dropping rotational invariance. + +This tutorial, based on a `recent paper `_, +explains how to test and control the impact of approximate rotational invariance in +MLPs for MD simulations. """ # %% +# Setting up +# ---------- +# +# First, we import things import os import subprocess @@ -17,7 +38,7 @@ import xml.etree.ElementTree as ET import ase -import ase.io +import ase.io, ase.build import chemiscope import ipi import numpy as np @@ -27,8 +48,9 @@ ) -# import matplotlib.pyplot as plt -# import numpy as np +import matplotlib.pyplot as plt +import numpy as np + # %% # show rotations @@ -89,63 +111,126 @@ ) # %% -# Show the problem -# ---------------------------------- +# Next, we set up the interface with LAMMPS, which has to be compiled. +# We use a utility function to compile it. Note that this requires +# a functioning build system with `gfortran` and `make`. + +ipi.install_driver() # %% +# Non-invariant water potential +# ----------------------------- +# +# Instead of a complex MLP, we'll be using an artificial example in this +# tutorial: We use the classical "tip4p" forcefield for water, and add an +# orientation-dependent additional term to simulate missing rotational invariance. # +# Let's have a look at that potential. We'll take a water, and rotate it +# around 180 degrees in steps of two degrees, and look at the energy: -# Open and show the relevant part of the input -xmlroot = ET.parse("data/input-noo3.xml").getroot() -print(" " + ET.tostring(xmlroot.find(".//ffrotations"), encoding="unicode")) +atoms = ase.build.molecule("H2O") +atoms.set_cell(np.eye(3) * 5) +atoms.positions += 2.5 +frames = [] +for i in range(0, 180, 2): + a = atoms.copy() + a.rotate(i, "z", center="COM") + frames.append(a) +print(len(frames)) +ase.io.write("single/replay.xyz", frames) +# todo: actually do the work, get the results, make plot # %% -# Installing the Python driver -# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -# -# i-PI comes with a FORTRAN driver, which however has to be installed -# from source. We use a utility function to compile it. Note that this requires -# a functioning build system with `gfortran` and `make`. - -ipi.install_driver() +# As we can see, the potential energy changes under rotations. We have a spurious +# angular dependence in our potential! # %% -# 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-noo3.out"): - ipi_process = subprocess.Popen(["i-pi", "data/input-noo3.xml"]) - time.sleep(2) # wait for i-PI to start - lmp_process = [subprocess.Popen(["lmp", "-in", "data/in.lmp"]) for i in range(1)] - driver_process = [ - subprocess.Popen( - ["i-pi-driver", "-u", "-a", "h2o-noo3", "-m", "noo3-h2o"], cwd="data/" - ) - for i in range(1) - ] +# Conservation of angular momentum +# -------------------------------- +# +# The most drastic consequence of a lack of rotational invariance is that angular +# momentum is no longer conserved. We can easily see this by simulating a single +# water molecule, initialised with zero total angular momentum. With our non- +# invariant potential, it starts rotating after just a few timesteps: + +# todo: actually do the things here. # %% -# Skip this cell if you want to run in the background -if ipi_process is not None: - ipi_process.wait() - lmp_process[0].wait() - driver_process[0].wait() +# In the manuscript linked at the beginning, we show that for molecules, this +# leads to a spurious preferred orientation of the molecules, even for NVT +# simulations. Interestingly, we find that for an approximately invariant MLP +# that was trained with data augmentation, there is no impact on observables for +# *liquid* water, which the model was trained on. +# However, this tutorial is focused on something else: What to do to control +# the lack of invariance after the fact. # %% +# Controlling non-invariance with averaging +# ----------------------------------------- # -# Discuss how molecules get oriented +# A not rotationally invariant potential can be made more invariant by *averaging* +# over rotations. Rotations are continuous, so we can't recover exact invariance, +# but we can get close enough for many practical cases, as we show in the paper above. +# +# I-Pi has *builtin* support for different symmetrisation schemes: Two different types +# of grids over rotations (Legendre, Lebedev) and a single randomly chosen rotation. +# They are invoked by using `ffrotations` in place of `ffsocket` to interface with the driver: -traj_data = ase.io.read("water-noo3.pos_0.extxyz", ":") +# Open and show the relevant part of the input +xmlroot = ET.parse("data/input-noo3.xml").getroot() +print(" " + ET.tostring(xmlroot.find(".//ffrotations"), encoding="unicode")) # %% -# then, assemble a visualization -chemiscope.show( - frames=traj_data, - settings=chemiscope.quick_settings(trajectory=True), - mode="structure", -) +# So, let's now run the same dynamics from above *with* this rotational augmentation. + +# todo: run, get results, make animation + +# %% +# +# With this grid, we've managed to successfully control the non-invariance: There is +# only a very slow remaining precession. In a less "strict" setting, for instance NVT, +# this would be more than sufficient to mitigate the impacts of the non-invariance. + +# scratch space below + +# # %% +# # 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-noo3.out"): +# ipi_process = subprocess.Popen(["i-pi", "data/input-noo3.xml"]) +# time.sleep(2) # wait for i-PI to start +# lmp_process = [subprocess.Popen(["lmp", "-in", "data/in.lmp"]) for i in range(1)] +# driver_process = [ +# subprocess.Popen( +# ["i-pi-driver", "-u", "-a", "h2o-noo3", "-m", "noo3-h2o"], cwd="data/" +# ) +# for i in range(1) +# ] + +# # %% +# # Skip this cell if you want to run in the background +# if ipi_process is not None: +# ipi_process.wait() +# lmp_process[0].wait() +# driver_process[0].wait() + + +# # %% +# # +# # Discuss how molecules get oriented + +# traj_data = ase.io.read("water-noo3.pos_0.extxyz", ":") + +# # %% +# # then, assemble a visualization +# chemiscope.show( +# frames=traj_data, +# settings=chemiscope.quick_settings(trajectory=True), +# mode="structure", +# )