Bosonic PIMD in quadratic time (i-pi#258)

* new implementation of bosonic pimd
* clean-up of the normalmodes code
* documentation and examples
* regtests

---------

Co-authored-by: Barak Hirshberg <[email protected]>
Co-authored-by: Michele Ceriotti <[email protected]>

README.md

@@ -90,11 +90,11 @@ The monitoring can be interrupted with CTRL+C when the run has finished (5000 st
 
 ### Run the automatic test suite
 
-The automatic test suite can be run by calling the i-pi-test script.
+The automatic test suite can be run by calling the i-pi-tests script.
 You need to have the `pytest` package installed
 
 ```
-i-pi-test
+i-pi-tests
 ```
 
 You may also need to install some dependencies, listed in `requirements.txt`.

doc/scripts/help.py

@@ -102,6 +102,7 @@
     "prng": prng.InputRandom(),
     "normal_modes": normalmodes.InputNormalModes(),
     "frequencies": normalmodes.InputNMFrequencies(),
+    "bosons": normalmodes.InputBosons(),
     "initialize": initializer.InputInitializer(),
     "file": initializer.InputInitFile(),
     "positions": initializer.InputInitPositions(),

examples/features/open_paths/one_path/run.sh

@@ -1,6 +1,6 @@
 ipi=i-pi
 driver=i-pi-driver 
-sleep_time=10
+sleep_time=4
 
 ${ipi} input.xml > log.i-pi & 
 echo "# i-PI is running"

examples/features/particle_exchange/README.md

@@ -0,0 +1,10 @@
+Exchange in path integral molecular dynamics. 
+=================================================
+Path integral molecular dynamics including bosonic exchange effects (ordinary path integral molecular dynamics assumes that particles are distinguishable). Exchange manifests in ring polymers that link several particles together.
+Fermionic statistics can also be obtained by reweighting bosonic simulations, when the sign problem is not too large (see Hirshberg et al., doi:10.1063/5.0008720).
+The algorithm is based on the evaluation of an effective spring potential that averages over multiple path connectivities, that scales quadratically with the number of particles and linearly with the number of beads, and is based on Hirshberg et al.'s doi:10.1073/pnas.1913365116 and Feldman and Hirshberg's doi:10.1063/5.0173749.
+
+`basic_trapped_bosons`: trapped bosons and a mixture of bosons and distinguishable particles
+`specify_by_label_and_estimators`: example of specifying bosons and kinetic energy estimators via labels
+`fermions`: obtaining fermionic statistics by reweighting bosonic simulations
+`exchange_probabilities`: properties that indicate the extent that exchange interaction is occurring in the simulation

examples/features/particle_exchange/basic_trapped_bosons/README.md

@@ -0,0 +1,8 @@
+Trapped cold non-interacting bosons
+=================================================
+`trapped_bosons_3`: 3 noninteracting bosons in a 3D anisotropic harmonic trap. The analytical total energy is 0.00058031 Ha.
+`trapped_bosons_2_distinguishable_1`: a mixture of 2 noninteracting bosons and 1 additional non-interacting distinguishable particle, in a 3D anisotropic harmonic trap. The analytical total energy is 0.00062348 Ha.
+
+Both differ from the result for 3 distinguishable particles, which is 0.00065142 Ha.
+
+The total energy is obtained by summing the properties 'potential' for the potential energy, and either -'virial_fq' or 'kinetic_td'.

examples/features/particle_exchange/basic_trapped_bosons/trapped_bosons_2_distinguishable_1/input.xml

@@ -0,0 +1,55 @@
+<simulation threading='False' verbosity='low'>
+   <ffsocket mode='unix' name='driver'>
+       <address>bosons-trapped</address>
+   </ffsocket>
+
+   <total_steps> 100 </total_steps>
+
+   <output prefix="data">
+      <trajectory stride="100" filename="pos" cell_units="angstrom">positions{angstrom}</trajectory>
+      <properties stride="100"> [ step, time{femtosecond}, conserved, temperature{kelvin}, potential, virial_fq, kinetic_td ] </properties>
+   </output>
+
+   <prng>
+      <seed> 18885 </seed>
+   </prng>
+
+   <system>
+
+      <forces> 
+          <force forcefield="driver"></force>
+      </forces>
+
+      <initialize nbeads="32">
+        <positions mode="manual" bead="0"> [78, -21, 58, 17, -84, -93, 52, -56, -13] </positions>
+	      <masses mode="manual"> [1.0, 1.0, 1.0] </masses>
+	      <labels mode="manual"> ['E', 'E', 'E'] </labels>
+        <cell>
+         [   2500, 0, 0, 0, 2500, 0, 0, 0, 2500 ]
+        </cell>
+	<velocities mode='thermal' units='kelvin'> 17.4 </velocities>
+      </initialize>
+
+      <normal_modes propagator='bab'>
+	      <nmts> 10 </nmts>
+	      <bosons> [0, 2] </bosons>
+      </normal_modes>
+
+      <ensemble>
+         <temperature units="kelvin"> 17.4 </temperature>
+      </ensemble>
+
+      <motion mode="dynamics">
+	<fixcom> False </fixcom>
+        <dynamics mode="nvt">
+         <timestep units="femtosecond"> 1 </timestep>
+  	  <thermostat mode='pile_l'>
+		<tau units='femtosecond'>100</tau>
+	  </thermostat>
+
+        </dynamics>
+      </motion>
+
+  </system>
+
+</simulation>

examples/features/particle_exchange/basic_trapped_bosons/trapped_bosons_2_distinguishable_1/run.sh

@@ -0,0 +1,16 @@
+ipi=i-pi
+driver="i-pi-driver -m harm3d -o 1.21647924E-8 -u -h bosons-trapped"
+sleep_time=4
+
+${ipi} input.xml > log.i-pi & 
+echo "# i-PI is running"
+
+echo "# Waiting for ${sleep_time} (s) before executing driver"
+sleep ${sleep_time}
+
+${driver} > /dev/null & 
+echo "# Driver is running"
+
+wait
+
+echo "# Simulation complete"

examples/features/particle_exchange/basic_trapped_bosons/trapped_bosons_3/input.xml

@@ -0,0 +1,55 @@
+<simulation threading='False' verbosity='low'>
+   <ffsocket mode='unix' name='driver'>
+       <address>bosons-trapped</address>
+   </ffsocket>
+
+   <total_steps> 100 </total_steps>
+
+   <output prefix="data">
+      <trajectory stride="100" filename="pos" cell_units="angstrom">positions{angstrom}</trajectory>
+      <properties stride="100"> [ step, time{femtosecond}, conserved, temperature{kelvin}, potential, virial_fq, kinetic_td ] </properties>
+   </output>
+
+   <prng>
+      <seed> 18885 </seed>
+   </prng>
+
+   <system>
+
+      <forces> 
+          <force forcefield="driver"></force>
+      </forces>
+
+      <initialize nbeads="32">
+        <positions mode="manual" bead="0"> [78, -21, 58, 17, -84, -93, 52, -56, -13] </positions>
+	      <masses mode="manual"> [1.0, 1.0, 1.0] </masses>
+	      <labels mode="manual"> ['E', 'E', 'E'] </labels>
+        <cell>
+         [   2500, 0, 0, 0, 2500, 0, 0, 0, 2500 ]
+        </cell>
+	<velocities mode='thermal' units='kelvin'> 17.4 </velocities>
+      </initialize>
+
+      <normal_modes propagator='bab'>
+	      <nmts> 10 </nmts>
+	      <bosons> [0, 1, 2] </bosons>
+      </normal_modes>
+
+      <ensemble>
+         <temperature units="kelvin"> 17.4 </temperature>
+      </ensemble>
+
+      <motion mode="dynamics">
+	<fixcom> False </fixcom>
+        <dynamics mode="nvt">
+         <timestep units="femtosecond"> 1 </timestep>
+  	  <thermostat mode='pile_l'>
+		<tau units='femtosecond'>100</tau>
+	  </thermostat>
+
+        </dynamics>
+      </motion>
+
+  </system>
+
+</simulation>

examples/features/particle_exchange/basic_trapped_bosons/trapped_bosons_3/run.sh

@@ -0,0 +1,16 @@
+ipi=i-pi
+driver="i-pi-driver -m harm3d -o 1.21647924E-8 -u -h bosons-trapped"
+sleep_time=4
+
+${ipi} input.xml > log.i-pi & 
+echo "# i-PI is running"
+
+echo "# Waiting for ${sleep_time} (s) before executing driver"
+sleep ${sleep_time}
+
+${driver} > /dev/null & 
+echo "# Driver is running"
+
+wait
+
+echo "# Simulation complete"

examples/features/particle_exchange/exchange_probabilities/README.md

@@ -0,0 +1,6 @@
+Bosonic exchange probabilities
+=================================================
+An example of recording the extent to which the dynamics are affected by bosonic exchange. Bosonic exchange creates a superposition of multiple ring polymer configurations, some of which combine the rings of several atoms together.
+The property 'exchange_distinct_prob' is the probability of the ordinary, distinguishable configuration, where each atom is in its own separate ring.
+When exchange has a significant effect, 'exchange_distinct_prob' is close to 0; when it has no effect, it is close to 1.
+The property 'exchange_all_prob' is the probability of the other extreme, the configuration that connects all the atoms together into a single ring, scaled by 1/N. In the limit of low temperatures it converges to 1.

examples/features/particle_exchange/exchange_probabilities/input.xml

@@ -0,0 +1,56 @@
+<simulation threading='False' verbosity='low'>
+   <ffsocket mode='unix' name='driver'>
+       <address>bosons-trapped</address>
+   </ffsocket>
+
+   <total_steps> 100 </total_steps>
+
+   <output prefix="data">
+      <trajectory stride="100" filename="pos" cell_units="angstrom">positions{angstrom}</trajectory>
+      <properties stride="100"> [ step, time{femtosecond}, conserved, temperature{kelvin}, potential, virial_fq, kinetic_td,
+                                  exchange_distinct_prob, exchange_all_prob ] </properties>
+   </output>
+
+   <prng>
+      <seed> 18885 </seed>
+   </prng>
+
+   <system>
+
+      <forces> 
+          <force forcefield="driver"></force>
+      </forces>
+
+      <initialize nbeads="32">
+        <positions mode="manual" bead="0"> [78, -21, 58, 17, -84, -93, 52, -56, -13] </positions>
+	      <masses mode="manual"> [1.0, 1.0, 1.0] </masses>
+	      <labels mode="manual"> ['E', 'E', 'E'] </labels>
+        <cell>
+         [   2500, 0, 0, 0, 2500, 0, 0, 0, 2500 ]
+        </cell>
+	<velocities mode='thermal' units='kelvin'> 17.4 </velocities>
+      </initialize>
+
+      <normal_modes propagator='bab'>
+	      <nmts> 10 </nmts>
+	      <bosons> [0, 1, 2] </bosons>
+      </normal_modes>
+
+      <ensemble>
+         <temperature units="kelvin"> 17.4 </temperature>
+      </ensemble>
+
+      <motion mode="dynamics">
+	<fixcom> False </fixcom>
+        <dynamics mode="nvt">
+         <timestep units="femtosecond"> 1 </timestep>
+  	  <thermostat mode='pile_l'>
+		<tau units='femtosecond'>100</tau>
+	  </thermostat>
+
+        </dynamics>
+      </motion>
+
+  </system>
+
+</simulation>

examples/features/particle_exchange/exchange_probabilities/run.sh

@@ -0,0 +1,16 @@
+ipi=i-pi
+driver="i-pi-driver -m harm3d -o 1.21647924E-8 -u -h bosons-trapped"
+sleep_time=4
+
+${ipi} input.xml > log.i-pi & 
+echo "# i-PI is running"
+
+echo "# Waiting for ${sleep_time} (s) before executing driver"
+sleep ${sleep_time}
+
+${driver} > /dev/null & 
+echo "# Driver is running"
+
+wait
+
+echo "# Simulation complete"

examples/features/particle_exchange/fermions/README.md

@@ -0,0 +1,6 @@
+Trapped cold non-interacting fermions
+=================================================
+An example of how to obtain fermionic statistics from a bosonic simulation, here, of 3 non-interacting fermions in an anisotropic harmonic trap.
+The fermions are treated as bosons throughout the simulation. An additional property, 'fermionic_sign', is recorded. Then, in the analysis, the value of the estimator should be multiplied by the fermionic sign in each step, and the total result should be divided the average fermionic sign. See Hirshberg et al.'s doi:10.1063/5.0008720.
+
+When the average sign is close to zero, the results become harder to converge. This is the reason for the choice of a higher temperature compared to other examples of trapped bosons. P=12 beads suffice for convergence in this temperature. The analytical result for the total energy is 0.00091881 Ha, higher than the result for distinguishable particles (0.00078146 Ha) or bosons (0.0006917 Ha) in the same setting.

examples/features/particle_exchange/fermions/input.xml

@@ -0,0 +1,55 @@
+<simulation threading='False' verbosity='low'>
+   <ffsocket mode='unix' name='driver'>
+       <address>bosons-trapped</address>
+   </ffsocket>
+
+   <total_steps> 100 </total_steps>
+
+   <output prefix="data">
+      <trajectory stride="100" filename="pos" cell_units="angstrom">positions{angstrom}</trajectory>
+      <properties stride="100"> [ step, time{femtosecond}, conserved, temperature{kelvin}, potential, virial_fq, kinetic_td, fermionic_sign ] </properties>
+   </output>
+
+   <prng>
+      <seed> 18885 </seed>
+   </prng>
+
+   <system>
+
+      <forces> 
+          <force forcefield="driver"></force>
+      </forces>
+
+      <initialize nbeads="12">
+        <positions mode="manual" bead="0"> [78, -21, 58, 17, -84, -93, 52, -56, -13] </positions>
+	      <masses mode="manual"> [1.0, 1.0, 1.0] </masses>
+	      <labels mode="manual"> ['E', 'E', 'E'] </labels>
+        <cell>
+         [   2500, 0, 0, 0, 2500, 0, 0, 0, 2500 ]
+        </cell>
+	<velocities mode='thermal' units='kelvin'> 23.22 </velocities>
+      </initialize>
+
+      <normal_modes propagator='bab'>
+	      <nmts> 10 </nmts>
+	      <bosons> [0, 1, 2] </bosons>
+      </normal_modes>
+
+      <ensemble>
+         <temperature units="kelvin"> 23.22 </temperature>
+      </ensemble>
+
+      <motion mode="dynamics">
+	<fixcom> False </fixcom>
+        <dynamics mode="nvt">
+         <timestep units="femtosecond"> 1 </timestep>
+  	  <thermostat mode='pile_l'>
+		<tau units='femtosecond'>100</tau>
+	  </thermostat>
+
+        </dynamics>
+      </motion>
+
+  </system>
+
+</simulation>

examples/features/particle_exchange/fermions/run.sh

@@ -0,0 +1,16 @@
+ipi=i-pi
+driver="i-pi-driver -m harm3d -o 1.21647924E-8 -u -h bosons-trapped"
+sleep_time=4
+
+${ipi} input.xml > log.i-pi & 
+echo "# i-PI is running"
+
+echo "# Waiting for ${sleep_time} (s) before executing driver"
+sleep ${sleep_time}
+
+${driver} > /dev/null & 
+echo "# Driver is running"
+
+wait
+
+echo "# Simulation complete"

examples/features/particle_exchange/specify_by_label_and_estimators/README.md

@@ -0,0 +1,3 @@
+Specifying bosons by labels and other estimator concers
+=================================================
+An example of how to specify which atoms participate in bosonic exchange not through their indices but in a list of labels. Here this indicates a mixture of distiguishable particles and bosons. The thermodynamic quantum kinetic energy estimator can be applied to a subset of the atoms, but it must not include some bosons but not others. Here, the value of 'kinetic_td' should always be the sum of 'kinetic_td(B)' (for the bosons) and 'kinetic_td(D)' (for the distinguishable particle). The centroid-virial quantum kinetic energy estimator cannot be applied to bosons; here it is applied to the distiguishable particle alone.