Vector gradient

qcelemental/models/align.py

@@ -126,6 +126,26 @@ def align_hessian(self, hess) -> NDArray:
         alhess = blockwise_contract(alhess)
         return alhess
 
+    def align_vector_gradient(self, mu_derivatives):
+        """Align the nuclear gradients of vector components (e.g. dipole derivatives)."""
+        # Input data is assumed to be organized into outermost x, y, z vector components.
+        # Organize derivatives for each atom into 3x3 and transform it.
+        mu_x, mu_y, mu_z = mu_derivatives
+        nat = mu_x.shape[0] // 3
+        al_mu = np.zeros( (3,3*nat) )
+
+        Datom = np.zeros( (3,3) ) # atom whose nuclear derivatives are taken
+        for at in range(nat):
+            Datom.fill(0)
+            Datom[0,:] = mu_x[3*self.atommap[at]:3*self.atommap[at]+3]
+            Datom[1,:] = mu_y[3*self.atommap[at]:3*self.atommap[at]+3]
+            Datom[2,:] = mu_z[3*self.atommap[at]:3*self.atommap[at]+3]
+            Datom[:] = np.dot( self.rotation.T , np.dot(Datom, self.rotation) )
+            al_mu[0, 3*at : 3*at+3] = Datom[0,:]
+            al_mu[1, 3*at : 3*at+3] = Datom[1,:]
+            al_mu[2, 3*at : 3*at+3] = Datom[2,:]
+        return al_mu
+
     def align_system(self, geom, mass, elem, elez, uniq, *, reverse: bool = False):
         """For AlignmentRecipe `ar`, apply its translation, rotation, and atom map."""
 

qcelemental/molutil/__init__.py

@@ -1 +1,2 @@
 from .align import B787, kabsch_align, compute_scramble
+

qcelemental/molutil/test_align.py

@@ -442,3 +442,78 @@ def test_hessian_align():
 
     p2chess = mill.align_hessian(p4_hooh_hess_native)
     assert compare_values(c4_hooh_hess, p2chess, atol=1.e-4)
+
+@using_networkx
+def test_vector_gradient_align():
+    # HOOH TS (optimized to be very nearly planar)
+    p4_hooh_xyz = """
+    units au
+    H  1.8327917647 -1.5752960165 -0.0000055594
+    O  1.3171390326  0.1388012713  0.0000003503
+    O -1.3171390326 -0.1388012713  0.0000003503
+    H -1.8327917647  1.5752960165 -0.0000055594
+    """
+
+    # from C4 GRD file, produced by p4_hooh_xyz in ZMAT
+    c4_hooh_xyz = """
+    units au
+    H  1.8319897007 -1.5762287045  -0.0000055594
+    H -1.8319897007  1.5762287045  -0.0000055594
+    O  1.3172095119  0.1381308288   0.0000003503
+    O -1.3172095119 -0.1381308288   0.0000003503
+    """
+
+    # From C4 DIPDER file, analytic
+    c4_hooh_dipder_x = np.array([
+        0.2780463810,      -0.0627423838,      -0.0000001663,
+        0.2780463810,      -0.0627423838,       0.0000001663,
+       -0.2780463810,       0.0627423838,       0.0000007872,
+       -0.2780463810,       0.0627423838,      -0.0000007872])
+    c4_hooh_dipder_y = np.array([
+       -0.0452364698,       0.2701572972,      -0.0000004246,
+       -0.0452364698,       0.2701572972,       0.0000004246,
+        0.0452364698,      -0.2701572972,       0.0000007936,
+        0.0452364698,      -0.2701572972,      -0.0000007936])
+    c4_hooh_dipder_z = np.array([
+       -0.0000001575,      -0.0000004725,       0.4019549601,
+        0.0000001575,       0.0000004725,       0.4019549601,
+       -0.0000000523,       0.0000008401,      -0.4019549601,
+        0.0000000523,      -0.0000008401,      -0.4019549601])
+
+    # Generated from fixing orientation/com in psi4.  Then using
+    # a 5-point finite differences of nuclear gradients computed
+    # with an applied electric field to produce a file17.dat.
+    p4_hooh_dipder_x  = np.array([
+        0.2781013514,      -0.0627383175,      -0.0000001660,
+       -0.2781013514,       0.0627383175,       0.0000007867,
+       -0.2781013514,       0.0627383175,      -0.0000007867,
+        0.2781013514,      -0.0627383175,       0.0000001660])
+    p4_hooh_dipder_y  = np.array([
+       -0.0452324587,       0.2701024305,      -0.0000004247,
+        0.0452324587,      -0.2701024305,       0.0000007939,
+        0.0452324587,      -0.2701024305,      -0.0000007939,
+       -0.0452324587,       0.2701024305,       0.0000004247])
+    p4_hooh_dipder_z  = np.array([
+       -0.0000001572,      -0.0000004726,       0.4019549470,
+       -0.0000000527,       0.0000008401,      -0.4019549470,
+        0.0000000527,      -0.0000008401,      -0.4019549470,
+        0.0000001572,       0.0000004726,       0.4019549470])
+    p4_hooh_dipder = np.concatenate((p4_hooh_dipder_x, p4_hooh_dipder_y,
+                                     p4_hooh_dipder_z)).reshape(3,-1)
+
+    p4mol = qcel.models.Molecule.from_data(p4_hooh_xyz)
+    c4mol = qcel.models.Molecule.from_data(c4_hooh_xyz)
+    aqmol, data = p4mol.align(c4mol, atoms_map=False, mols_align=True, verbose=4)
+    mill = data['mill']
+
+    assert compare([0, 3, 1, 2], mill.atommap)
+
+    p2cgeom = mill.align_coordinates(p4mol.geometry)
+    assert compare_values(c4mol.geometry, p2cgeom, atol=1.e-6)
+
+    p2c_dipder_x, p2c_dipder_y, p2c_dipder_z = mill.align_vector_gradient(p4_hooh_dipder)
+
+    assert compare_values(c4_hooh_dipder_x, p2c_dipder_x, atol=1.e-5)
+    assert compare_values(c4_hooh_dipder_y, p2c_dipder_y, atol=1.e-5)
+    assert compare_values(c4_hooh_dipder_z, p2c_dipder_z, atol=1.e-5)
+