Added an example, and options, to use lebedev grids in o3-averaging

examples/clients/xtb/README.md

@@ -15,9 +15,5 @@ the example can be run as usual
 
 ```bash
 i-pi input.xml &> log.ipi &
-<<<<<<< HEAD
 i-pi-py_driver -m xtb -o template=ch4.xyz,method=GFN2-xTB -u -a xtb &> log.xtb
-=======
-i-pi-py_driver -m xtb -o '{"method": "GFN2-xTB", "numbers": [6,1,1,1,1], "periodic": false}' -u -a xtb &> log.xtb &
->>>>>>> d72ab1f8 (Recovered weird changes wrt main)
 ```

examples/clients/xtb/run.sh

@@ -3,9 +3,5 @@ echo "Running i-PI"
 i-pi input.xml &> log.ipi &
 sleep 1
 echo "Running driver"
-<<<<<<< HEAD
 i-pi-py_driver -m xtb -o template=ch4.xyz,method=GFN2-xTB -u -a xtb &> log.xtb 
-=======
-i-pi-py_driver -m xtb -o '{"method": "GFN2-xTB", "numbers": [6,1,1,1,1], "periodic": false}' -u -a xtb &> log.xtb 
->>>>>>> d72ab1f8 (Recovered weird changes wrt main)
 

examples/features/o3_averaging/grid_liquid_water/input.xml

@@ -14,7 +14,8 @@
   </ffsocket>
   <ffrotations name='driver-noo3' mode='unix' pbc='false'>
       <address>h2o-noo3</address>
-      <grid> 2 </grid>
+      <grid_order> 2 </grid_order> 
+      <grid_mode> legendre </grid_mode>
       <inversion> True </inversion>
   </ffrotations>
   <system>

examples/features/o3_averaging/lebedev_liquid_water/input.xml

@@ -0,0 +1,41 @@
+<simulation verbosity='medium'>
+  <output prefix='simulation'>
+    <properties stride='1' filename='out'>  [ step, time{picosecond}, conserved{electronvolt}, temperature{kelvin}, kinetic_md{electronvolt}, potential{electronvolt}, pressure_md{megapascal}, stress_md{megapascal} ] </properties>
+    <trajectory filename='noo3_pots' stride='2' cell_units='angstrom' extra_type="o3grid_pots"> extras_component_raw(1) </trajectory>  
+    <trajectory filename='pos' stride='20' cell_units='angstrom'> positions{angstrom} </trajectory>
+    <checkpoint stride='200'/>
+  </output>
+  <total_steps>1000</total_steps>
+  <prng>
+    <seed>31415</seed>
+  </prng>
+  <ffsocket name='driver-qtip4pf' mode='unix' pbc='false'>
+      <address>h2o-base</address>
+  </ffsocket>
+  <ffrotations name='driver-noo3' mode='unix' pbc='false'>
+      <address>h2o-noo3</address>
+      <grid_order> 3 </grid_order> <grid_mode> lebedev </grid_mode>
+      <inversion> True </inversion>
+  </ffrotations>
+  <system>
+    <initialize nbeads='1'>
+      <file mode='xyz'> water_216.xyz </file>
+      <velocities mode='thermal' units='kelvin'> 300 </velocities>
+    </initialize>
+    <forces>
+      <force forcefield='driver-qtip4pf'/>
+      <force forcefield='driver-noo3' weight="1.0"/>
+    </forces>
+    <motion mode='dynamics'>
+      <dynamics mode='nvt'>
+        <timestep units='femtosecond'> 0.5 </timestep>
+        <thermostat mode='langevin'>
+          <tau units='femtosecond'> 100 </tau>
+        </thermostat>
+      </dynamics>
+    </motion>
+    <ensemble>
+      <temperature units='kelvin'> 300 </temperature>
+    </ensemble>
+  </system>
+</simulation>

examples/features/o3_averaging/lebedev_liquid_water/run.sh

@@ -0,0 +1,18 @@
+ipi=i-pi
+driver1="i-pi-driver -m noo3-h2o -u -a h2o-noo3"
+driver2="i-pi-driver -m qtip4pf -u -a h2o-base"
+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}
+
+${driver1} > /dev/null & 
+${driver2} > /dev/null &
+echo "# Driver is running"
+
+wait
+
+echo "# Simulation complete"

examples/features/o3_averaging/lebedev_liquid_water/test_settings.dat

@@ -0,0 +1,4 @@
+#Settings file for automatic test
+driver_model qtip4pf
+driver_model2 noo3-h2o
+nsteps 10

examples/features/o3_averaging/lebedev_liquid_water/water_216.xyz

@@ -0,0 +1 @@
+../../../init_files/water_216.xyz

ipi/engine/forcefields.py

@@ -26,7 +26,7 @@
 from ipi.utils.units import unit_to_internal
 from ipi.utils.distance import vector_separation
 from ipi.pes import __drivers__
-from ipi.utils.mathtools import get_rotation_quadrature, random_rotation
+from ipi.utils.mathtools import get_rotation_quadrature_legendre, get_rotation_quadrature_lebedev, random_rotation
 
 plumed = None
 
@@ -1462,7 +1462,8 @@ def __init__(
         interface=None,
         random=False,
         inversion=False,
-        grid=1,
+        grid_order=1,
+        grid_mode="lebedev"
     ):
         super(FFRotations, self).__init__(
             latency,
@@ -1479,9 +1480,19 @@ def __init__(
         self.prng = Random()
         self.random = random
         self.inversion = inversion
-        self.grid = grid
+        self.grid_order = grid_order
+        self.grid_mode = grid_mode
 
-        self._rotations = get_rotation_quadrature(self.grid)
+        if self.grid_mode == "lebedev":
+            self._rotations = get_rotation_quadrature_lebedev(self.grid_order)
+        elif self.grid_mode == "legendre":
+            self._rotations = get_rotation_quadrature_legendre(self.grid_order)
+        else:
+            raise ValueError(f"Invalid quadrature {self.grid_mode}")
+        info(f"""
+# Generating {self.grid_mode} rotation quadrature of order {self.grid_order}.
+# Grid contains {len(self._rotations)} proper rotations.
+""", verbosity.low)
 
     def queue(self, atoms, cell, reqid=-1):
 
@@ -1493,7 +1504,7 @@ def queue(self, atoms, cell, reqid=-1):
         else:
             R_random = np.eye(3)
 
-        for R, w in self._rotations:
+        for R, w, _ in self._rotations:
             R = R @ R_random
 
             rot_atoms = atoms.clone()

ipi/inputs/forcefields.py

@@ -473,7 +473,7 @@ class InputFFDebye(InputForceField):
                 "default": input_default(factory=np.zeros, args=(0,)),
                 "help": "Specifies the Hessian of the harmonic potential. "
                 "Default units are atomic. Units can be specified only by xml attribute. "
-                r"Implemented options are: 'atomic_unit', 'ev/ang\^2'",
+                r"Implemented options are: 'atomic_unit', 'ev/ang^2'",
                 "dimension": "hessian",
             },
         ),
@@ -943,18 +943,33 @@ class InputFFRotations(InputFFSocket):
         },
     )
 
-    fields["grid"] = (
+    fields["grid_order"] = (
         InputValue,
         {
             "dtype": int,
             "default": 1,
             "help": """Sums over a grid of rotations of the given order.
-            Note that the number of rotations increases rapidly with the order, e.g.
-            '1' leads to a single rotation, '2' to 18, '3' to 75 rotations.
+            Note that the number of rotations increases rapidly with the order:
+            e.g. for a legendre grid
+            '1' leads to a single rotation, '2' to 18, '3' to 75 rotations, while 
+            for a lebedev grid '3' contains 18 rotations, '5' 70 rotations and so on.
             """,
         },
     )
 
+    fields["grid_mode"] = (
+        InputValue,
+        {
+            "dtype": str,
+            "options": ["legendre", "lebedev"],
+            "default": "legendre",
+            "help": """Defines the type of integration grid. 
+            Lebedev grids are usually more efficient in integrating. 
+            """,
+        },
+    )
+
+
     fields["inversion"] = (
         InputValue,
         {
@@ -967,11 +982,12 @@ class InputFFRotations(InputFFSocket):
     )
 
     def store(self, ff):
-        """Store all the sub-forcefields"""
+        """Store the base and rotation quadrature parameters"""
 
         super(InputFFRotations, self).store(ff)
         self.inversion.store(ff.inversion)
-        self.grid.store(ff.grid)
+        self.grid_order.store(ff.grid_order)
+        self.grid_mode.store(ff.grid_mode)
         self.random.store(ff.random)
 
     def fetch(self):
@@ -994,7 +1010,8 @@ def fetch(self):
                 match_mode=self.matching.fetch(),
                 exit_on_disconnect=self.exit_on_disconnect.fetch(),
             ),
-            grid=self.grid.fetch(),
+            grid_order=self.grid_order.fetch(),
+            grid_mode=self.grid_mode.fetch(),
             random=self.random.fetch(),
             inversion=self.inversion.fetch(),
         )

ipi/pes/xtb.py

@@ -1,38 +1,15 @@
-<<<<<<< HEAD:ipi/pes/xtb.py
 import numpy as np
 
 from ipi.utils.units import unit_to_internal, unit_to_user
 from ipi.utils.io import read_file_name
 from .dummy import Dummy_driver
 
 tb = None
-=======
-"""
-Interface with tblite to provide GFN1-xTB and GFN2-xTB calculators.
-
-Example:: 
-
-    python driver.py -m xtb -u -o '{"method": "GFN2-xTB", "numbers": [8,1,1], "periodic": false}'
-"""
-
-import numpy as np
-import json
-
-from ipi.utils.units import unit_to_internal, unit_to_user
-
-try:
-    import tblite.interface as tb
-except ImportError as e:
-    raise ModuleNotFoundError(
-        "Could not find tblite for xtb driver. Please install tblite-python with mamba"
-    ) from e
->>>>>>> d72ab1f8 (Recovered weird changes wrt main):drivers/py/pes/xtb.py
 
 __DRIVER_NAME__ = "xtb"
 __DRIVER_CLASS__ = "TBLiteDriver"
 
 
-<<<<<<< HEAD:ipi/pes/xtb.py
 class TBLiteDriver(Dummy_driver):
     """
     Interface with tblite to provide GFN1-xTB and GFN2-xTB calculators.
@@ -76,27 +53,6 @@ def __init__(
             np.repeat(self.periodic, 3) if self.periodic else None,
         )
         self.calc.set("verbosity", self.verbosity)
-=======
-class TBLiteDriver(object):
-    """Base class providing the structure of a PES for the python driver."""
-
-    def __init__(
-        self,
-        args="",
-        verbose=False,
-    ):
-        """Initialized dummy drivers"""
-        config = json.loads(args)
-        try:
-            self.method = config["method"]
-            self.numbers = np.asarray(config["numbers"])
-        except KeyError as e:
-            raise ValueError(f"Required key {str(e)} not found.")
-        self.charge = config.get("charge")
-        self.uhf = config.get("uhf")
-        self.periodic = config.get("periodic")
-        self.verbosity = 1 if verbose else 0
->>>>>>> d72ab1f8 (Recovered weird changes wrt main):drivers/py/pes/xtb.py
 
     def __call__(self, cell, pos):
         """
@@ -109,26 +65,10 @@ def __call__(self, cell, pos):
         pos = unit_to_user("length", "atomic_unit", pos)
         cell = unit_to_user("length", "atomic_unit", cell.T)
 
-<<<<<<< HEAD:ipi/pes/xtb.py
         self.calc.update(np.asarray(pos), np.asarray(cell) if self.periodic else None)
 
         # Do the actual calculation
         results = self.calc.singlepoint()
-=======
-        calc = tb.Calculator(
-            self.method,
-            self.numbers,
-            np.asarray(pos),
-            self.charge,
-            self.uhf,  # unpaired electrons
-            np.asarray(cell) if self.periodic else None,
-            np.repeat(self.periodic, 3) if self.periodic else None,
-        )
-        calc.set("verbosity", self.verbosity)
-
-        # Do the actual calculation
-        results = calc.singlepoint()
->>>>>>> d72ab1f8 (Recovered weird changes wrt main):drivers/py/pes/xtb.py
         pot = results.get("energy")
         grad = results.get("gradient")
         vir = results.get("virial")

ipi/utils/mathtools.py

@@ -6,6 +6,7 @@
 
 
 import math
+from importlib import resources
 
 import numpy as np
 
@@ -25,6 +26,8 @@
     "logsumlog",
     "sinch",
     "gaussian_inv",
+    "get_rotation_quadrature_legendre",
+    "get_rotation_quadrature_lebedev",
 ]
 
 
@@ -309,7 +312,6 @@ def det_ut3x3(h):
 
 MINSERIES = 1e-8
 
-
 def exp_ut3x3(h):
     """Computes the matrix exponential for a 3x3 upper-triangular matrix.
 
@@ -597,7 +599,7 @@ def roots_legendre(L):
     return roots, weights
 
 
-def get_rotation_quadrature(L):
+def get_rotation_quadrature_legendre(L):
     if L == 1:
         # returns the identity (for some reason this algo below generates a different rotation)
         return [(np.eye(3), 2.0)]
@@ -611,10 +613,24 @@ def get_rotation_quadrature(L):
             for v, weight in zip(all_v, weights_now):
                 angles = [theta, v, w]
                 rotation_matrix = euler_zxz_to_matrix(*angles)
-                quads.append((rotation_matrix, weight))
+                quads.append((rotation_matrix, weight, angles))
 
     return quads
 
+def get_rotation_quadrature_lebedev(L):    
+    with resources.path("ipi.utils", "lebedev_grids.npy") as file_path:
+        lebedev = np.load(file_path, allow_pickle=True).item()
+    if not L in lebedev:
+        raise ValueError(f"There is no Lebedev grid of order {L} available")
+    leb_quad = lebedev[L]
+    quads = []
+    for theta_index in range(0, L):
+        theta = 2 * np.pi * theta_index / L
+        for w, v, weight in leb_quad:
+            angles = [theta, v, w] # [w, v, theta]
+            rotation_matrix = euler_zxz_to_matrix(*angles)
+            quads.append((rotation_matrix, weight, angles))
+    return quads   
 
 def random_rotation(prng, improper=True):
     """Generates a (uniform) random rotation matrix"""

setup.cfg

@@ -21,10 +21,10 @@ classifiers =
     Operating System :: POSIX :: Linux
     License :: OSI Approved :: GNU General Public License v3 (GPLv3)
     Programming Language :: Python :: 3
-    Programming Language :: Python :: 3.6
     Programming Language :: Python :: 3.7
     Programming Language :: Python :: 3.8
     Programming Language :: Python :: 3.9
+    Programming Language :: Python :: 3.10
     Topic :: Scientific/Engineering :: Chemistry
     Topic :: Scientific/Engineering :: Physics
 

setup.py

@@ -7,5 +7,8 @@
 
 setup(
     packages=[*find_packages(exclude=["ipi_tests*"])],
+    package_data={
+        "ipi.utils": ["lebedev_grids.npy"],  # Include the npy file here
+    },
     scripts=scripts,
 )