Working on the cell...

drivers/f90/driver.f90

@@ -191,15 +191,15 @@ PROGRAM DRIVER
                   vstyle = 64
                ELSEIF (trim(cmdbuffer) == "qtip4pf-c-json-delta") THEN
                   vstyle = 65
-               ELSEIF (trim(cmdbuffer) == "qtip4pf-noo3") THEN
+               ELSEIF (trim(cmdbuffer) == "noo3-h2o") THEN
                   vstyle = 70
                ELSEIF (trim(cmdbuffer) == "gas") THEN
                   vstyle = 0  ! ideal gas
                ELSEIF (trim(cmdbuffer) == "dummy") THEN
                   vstyle = 99 ! returns non-zero but otherwise meaningless values
                ELSE
                   WRITE(*,*) " Unrecognized potential type ", trim(cmdbuffer)
-                  WRITE(*,*) " Use -m [dummy|gas|lj|sg|harm|harm3d|morse|morsedia|zundel|qtip4pf|pswater|lepsm1|lepsm2|qtip4pf-efield|eckart|ch4hcbe|ljpolymer|MB|doublewell|doublewell_1D|water_dip_pol|harmonic_bath|meanfield_bath|ljmix|qtip4pf-sr|qtip4pf-c-1|qtip4pf-c-2|qtip4pf-c-json|qtip4pf-c-1-delta|qtip4pf-c-2-delta|qtip4pf-c-json-delta|qtip4pf-noo3] "
+                  WRITE(*,*) " Use -m [dummy|gas|lj|sg|harm|harm3d|morse|morsedia|zundel|qtip4pf|pswater|lepsm1|lepsm2|qtip4pf-efield|eckart|ch4hcbe|ljpolymer|MB|doublewell|doublewell_1D|water_dip_pol|harmonic_bath|meanfield_bath|ljmix|qtip4pf-sr|qtip4pf-c-1|qtip4pf-c-2|qtip4pf-c-json|qtip4pf-c-1-delta|qtip4pf-c-2-delta|qtip4pf-c-json-delta|noo3-h2o] "
                   STOP "ENDED"
                ENDIF
             ELSEIF (ccmd == 4) THEN
@@ -793,28 +793,24 @@ PROGRAM DRIVER
                   CALL LJ_getall(rc, 2.5d0, 1d-6, nat, atoms, cell_h, cell_ih, index_list, n_list, pot, forces, virial)
                ENDIF
             ELSEIF (70 == vstyle) THEN 
-               ! qtip4pf potential with rotational noise. adds a orientation-dependent term to 
+               ! potential that can be applied on top of water molecules to add a orientation-dependent term to 
                ! test non-equivariant terms in the potential
                vpars(1) = cell_h(1,1)
                vpars(2) = cell_h(2,2)
                vpars(3) = cell_h(3,3)
-               IF (cell_h(1,2).gt.1d-10 .or. cell_h(1,3).gt.1d-12  .or. cell_h(2,3).gt.1d-12) THEN
-                  WRITE(*,*) " qtip4pf PES only works with orthorhombic cells", cell_h(1,2), cell_h(1,3), cell_h(2,3)
-                  STOP "ENDED"
-               ENDIF
-               CALL qtip4pf(vpars(1:3),atoms,nat,forces,pot,virial)
-
+
+               pot = 0
+               forces = 0
                DO i=1,nat,3
                   dip = 0.5 * (atoms(i+1,:) + atoms(i+2,:)) - atoms(i,:)
-                  vpars(4) = sqrt(dip(1)**2+dip(2)**2+dip(3)**2)
-                  dip = dip/sqrt(dip(1)**2+dip(2)**2+dip(3)**2)
+                  vpars(4) = sqrt(dip(1)**2+dip(2)**2+dip(3)**2)                  
                   vpars(1) = exp(dip(1)/vpars(4)*1.0)
                   vpars(2) = exp(dip(2)/vpars(4)*0.5)
                   vpars(3) = exp(dip(3)/vpars(4)*2.0)
-                  pot = pot + 1e-2*(vpars(1)+vpars(2)+vpars(3))
-                  vpars(1) = 1e-2*vpars(1)*1.0/vpars(4)**3
-                  vpars(2) = 1e-2*vpars(2)*0.5/vpars(4)**3
-                  vpars(3) = 1e-2*vpars(3)*2.0/vpars(4)**3
+                  pot = pot + 2e-2*(vpars(1)+vpars(2)+vpars(3))
+                  vpars(1) = 2e-2*vpars(1)*1.0/vpars(4)**3
+                  vpars(2) = 2e-2*vpars(2)*0.5/vpars(4)**3
+                  vpars(3) = 2e-2*vpars(3)*2.0/vpars(4)**3
                   ! gradients on O
                   vpars(4) = -vpars(1)*(dip(2)**2+dip(3)**2) + vpars(2)*dip(1)*dip(2) + vpars(3)*dip(3)*dip(1) 
                   vpars(5) = -vpars(2)*(dip(1)**2+dip(3)**2) + vpars(1)*dip(1)*dip(2) + vpars(3)*dip(2)*dip(3) 
@@ -1097,7 +1093,7 @@ PROGRAM DRIVER
       SUBROUTINE helpmessage
          ! Help banner
 
-         WRITE(*,*) " SYNTAX: driver.x [-u] -a address [-p port] -m [dummy|gas|lj|sg|harm|harm3d|morse|morsedia|zundel|qtip4pf|pswater|lepsm1|lepsm2|qtip4p-efield|eckart|ch4hcbe|ljpolymer|MB|doublewell|doublewell_1D|water_dip_pol|harmonic_bath|meanfield_bath|ljmix|qtip4pf-sr|qtip4pf-c-1|qtip4pf-c-2|qtip4pf-c-json|qtip4pf-c-1-delta|qtip4pf-c-2-delta|qtip4pf-c-json-delta|qtip4pf-noo3]"
+         WRITE(*,*) " SYNTAX: driver.x [-u] -a address [-p port] -m [dummy|gas|lj|sg|harm|harm3d|morse|morsedia|zundel|qtip4pf|pswater|lepsm1|lepsm2|qtip4p-efield|eckart|ch4hcbe|ljpolymer|MB|doublewell|doublewell_1D|water_dip_pol|harmonic_bath|meanfield_bath|ljmix|qtip4pf-sr|qtip4pf-c-1|qtip4pf-c-2|qtip4pf-c-json|qtip4pf-c-1-delta|qtip4pf-c-2-delta|qtip4pf-c-json-delta|noo3-h2o]"
          WRITE(*,*) "         -o 'comma_separated_parameters' [-S sockets_prefix] [-v] "
          WRITE(*,*) ""
          WRITE(*,*) " For LJ potential use -o sigma,epsilon,cutoff "

examples/features/o3_averaging/random_liquid_water/input.xml

@@ -1,15 +1,19 @@
 <simulation verbosity='high'>
   <output prefix='simulation'>
-    <properties stride='1' filename='out'>  [ step, time{picosecond}, conserved{electronvolt}, temperature{kelvin}, kinetic_cv{electronvolt}, potential{electronvolt}, pressure_cv{megapascal} ] </properties>
+      <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='pos' stride='20' cell_units='angstrom'> positions{angstrom} </trajectory>
     <checkpoint stride='200'/>
   </output>
-  <total_steps>100</total_steps>
+  <total_steps>1000</total_steps>
   <prng>
     <seed>31415</seed>
   </prng>
-  <ffrotations name='driver-qtip4pf' mode='unix' pbc='false'>
-      <address>h2o-md.1</address>
+  <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>
+      <random> True </random>
       <grid> 2 </grid>
   </ffrotations>
   <system>
@@ -18,7 +22,8 @@
       <velocities mode='thermal' units='kelvin'> 300 </velocities>
     </initialize>
     <forces>
-      <force forcefield='driver-qtip4pf'> </force>
+      <force forcefield='driver-qtip4pf'/>
+      <force forcefield='driver-noo3' weight="1.0"/>
     </forces>
     <motion mode='dynamics'>
       <dynamics mode='nvt'>

examples/features/o3_averaging/random_liquid_water/run.sh

@@ -1,5 +1,6 @@
 ipi=i-pi
-driver="i-pi-driver -m qtip4pf -u -a h2o-md.1"
+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 & 
@@ -8,7 +9,8 @@ echo "# i-PI is running"
 echo "# Waiting for ${sleep_time} (s) before executing driver"
 sleep ${sleep_time}
 
-${driver} > /dev/null & 
+${driver1} > /dev/null & 
+${driver2} > /dev/null &
 echo "# Driver is running"
 
 wait

examples/temp/lammps/h2o-piglet.2/data.lmp

@@ -11,6 +11,7 @@ LAMMPS Description
          0 35.233 xlo xhi
          0 35.233 ylo yhi
          0 35.233 zlo zhi
+         1e-4 1e-4 1e-4 xy xz yz
 
 Masses
 

ipi/engine/atoms.py

@@ -148,7 +148,7 @@ def __init__(self, natoms, _prebind=None):
             name="kstress", func=self.get_kstress, dependencies=[self._p, self._m]
         )
 
-    def copy(self):
+    def clone(self):
         """Creates a new Atoms object.
 
         Returns:

ipi/engine/forcefields.py

@@ -16,6 +16,7 @@
 
 import numpy as np
 
+from ipi.utils.prng import Random
 from ipi.utils.softexit import softexit
 from ipi.utils.messages import info, verbosity, warning
 from ipi.interfaces.sockets import InterfaceSocket
@@ -1402,6 +1403,8 @@ def __init__(
             interface,
         )
 
+        # TODO: initialize and save (probably better to use different PRNG than the one from the simulation)
+        self.prng = Random()
         self.random = random
         self.inversion = inversion
         self.grid = grid
@@ -1411,16 +1414,25 @@ def __init__(
     def queue(self, atoms, cell, reqid=-1):
         # launches requests for all of the rotations FF objects
         ffh = []
-        for r in self._rotations:
-            print("submitting rotation ", r)
-            rot_atoms = atoms.copy()
-            rot_cell = cell.copy()  # todo update to clone
+        rots = []
+        if self.random:
+            R_random = random_rotation(self.prng, improper=True)
+
+        for R, w in self._rotations:
+            print("submitting rotation ", R, R_random)
+
+            R = R@R_random
+            rot_atoms = atoms.clone()
+            rot_cell = cell.clone()
+            rot_atoms.q[:] = (dstrip(rot_atoms.q).reshape(-1,3)@R.T).flatten()
+            rot_cell.h[:] = dstrip(rot_cell.h)@R.T
+            rots.append((R, w))
             ffh.append(super(FFRotations, self).queue(rot_atoms, rot_cell, reqid))
 
         # creates the request with the help of the base class,
         # making sure it already contains a handle to the list of FF
         # requests
-        req = ForceField.queue(self, atoms, cell, reqid, template=dict(ff_handles=ffh))
+        req = ForceField.queue(self, atoms, cell, reqid, template=dict(ff_handles=ffh, rots=rots))
         req["status"] = "Running"
         req["t_dispatched"] = time.time()
         return req
@@ -1436,7 +1448,6 @@ def check_finish(self, r):
 
     def gather(self, r):
         """Collects results from all sub-requests, and assemble the committee of models."""
-        print("GATHERING ")
         r["result"] = [
             0.0,
             np.zeros(len(r["pos"]), float),
@@ -1445,28 +1456,32 @@ def gather(self, r):
         ]
 
         # list of pointers to the forcefield requests. shallow copy so we can remove stuff
-        com_handles = r["ff_handles"].copy()
+        rot_handles = r["ff_handles"].copy()
+        rots = r["rots"].copy()
 
         # Gathers the forcefield energetics and extras
         pots = []
         frcs = []
         virs = []
         xtrs = []
-
-        for ff_r in com_handles:
-            pots.append(ff_r["result"][0])
-            frcs.append(ff_r["result"][1])
-            virs.append(ff_r["result"][2])
+        quad_w = []
+        for ff_r, (R, w) in zip(rot_handles, rots):
+            pots.append(ff_r["result"][0])            
+            frcs.append((ff_r["result"][1].reshape(-1,3)@R).flatten())
+            print("force component", frcs[-1][:3])
+            virs.append((R.T@ff_r["result"][2]@R))
             xtrs.append(ff_r["result"][3])
+            quad_w.append(w)
 
+        quad_w = np.array(quad_w)
         pots = np.array(pots)
         frcs = np.array(frcs).reshape(len(frcs), -1)
         virs = np.array(virs).reshape(-1, 3, 3)
 
         # Computes the mean energetics
-        mean_pot = np.mean(pots, axis=0)
-        mean_frc = np.mean(frcs, axis=0)
-        mean_vir = np.mean(virs, axis=0)
+        mean_pot = np.mean(pots*quad_w, axis=0)/quad_w.mean()
+        mean_frc = np.mean(frcs*quad_w[:,np.newaxis], axis=0)/quad_w.mean()
+        mean_vir = np.mean(virs*quad_w[:,np.newaxis,np.newaxis], axis=0)/quad_w.mean()
         # Sets the output of the committee model.
         r["result"][0] = mean_pot
         r["result"][1] = mean_frc
@@ -1479,10 +1494,8 @@ def gather(self, r):
             for x in xtrs:
                 r["result"][3][k].append(x[k])
 
-        print("RELEASING")
         for ff_r in r["ff_handles"]:
             self.release(ff_r)
-        print("RESEASIG SELF ")
         r["status"] = "Done"
         self.release(r)
 
@@ -1491,7 +1504,6 @@ def poll(self):
 
         super().poll()
 
-        print(len(self.requests))
         for r in self.requests:
             if "ff_handles" in r and r["status"] != "Done" and self.check_finish(r):
                 r["t_finished"] = time.time()

ipi/interfaces/sockets.py

@@ -842,7 +842,6 @@ def pool_distribute(self):
         if len(self.prlist) == 0 or len(freec) > len(self.prlist):
             self.prlist = [r for r in self.requests if r["status"] == "Queued"]
 
-        print(len(self.prlist), len(self.requests))
         if self.match_mode == "auto":
             match_seq = ["match", "none", "free", "any"]
         elif self.match_mode == "any":

ipi/utils/mathtools.py

@@ -598,6 +598,9 @@ def roots_legendre(L):
 
 
 def get_rotation_quadrature(L):
+    if L==1:
+        # returns the identity (for some reason this algo below generates a different rotation)
+        return [(np.eye(3), 2.0)]
     quads = []
     for theta_index in range(0, 2 * L - 1):
         for w_index in range(0, 2 * L - 1):