Fixing ammonia molecule and adding symmetric orthogonalization

meta.yaml

@@ -1,6 +1,6 @@
 package:
   name: "pyqint"
-  version: "0.14.0.1"
+  version: "0.14.1.0"
 
 source:
   path: .

pyqint/_version.py

@@ -1,2 +1,2 @@
-__version__ = "0.14.0.1"
+__version__ = "0.14.1.0"
 

pyqint/hf.py

@@ -14,7 +14,7 @@ class HF:
     """
     def rhf(self, mol, basis, calc_forces=False, itermax=100,
             use_diis=True, verbose=False, tolerance=1e-7,
-            orbc_init=None):
+            orbc_init=None, ortho='canonical'):
         """
         Performs a Hartree-Fock type calculation
 
@@ -51,7 +51,14 @@ def rhf(self, mol, basis, calc_forces=False, itermax=100,
         s, U = np.linalg.eigh(S)
 
         # construct transformation matrix X
-        X = U.dot(np.diag(1.0/np.sqrt(s)))
+        if ortho == 'canonical': # perform canonical orthogonalization
+            X = U @ np.diag(1.0/np.sqrt(s))
+        elif ortho == 'symmetric':
+            X = U @ np.diag(1.0/np.sqrt(s)) @ U.transpose()
+        else:
+            raise Exception("Invalid orthogonalization option selected: ", ortho)
+
+
 
         # create empty P matrix as initial guess
         if orbc_init is None:
@@ -83,9 +90,9 @@ def rhf(self, mol, basis, calc_forces=False, itermax=100,
                     continue
 
                 F = self.extrapolate_fock_from_diis_coefficients(fmats_diis, diis_coeff)
-                Fprime = X.transpose().dot(F).dot(X)
+                Fprime = X.transpose() @ F @ X
                 e, Cprime = np.linalg.eigh(Fprime)
-                C = X.dot(Cprime)
+                C = X @ Cprime
                 P = np.einsum('ik,jk,k->ij', C, C, occ)
 
             # calculate G
@@ -100,13 +107,13 @@ def rhf(self, mol, basis, calc_forces=False, itermax=100,
             F = T + V + G
 
             # transform Fock matrix
-            Fprime = X.transpose().dot(F).dot(X)
+            Fprime = X.transpose() @ F @ X
 
             # diagonalize F
             orbe, Cprime = np.linalg.eigh(Fprime)
 
             # back-transform
-            C = X.dot(Cprime)
+            C = X @ Cprime
 
             # calculate energy E
             energy = 0.0

pyqint/molecules/nh3.xyz

@@ -1,7 +1,6 @@
-5
+4
 
-  H        0.7034100        0.7034100        0.7034100
-  H       -0.7034100       -0.7034100        0.7034100
-  H       -0.7034100        0.7034100       -0.7034100
-  H        0.7034100       -0.7034100       -0.7034100
-  N        0.0000000        0.0000000        0.0000000
+N  0.000   0.000   0.000
+H  0.000  -0.937  -0.3816
+H  0.812   0.469  -0.3816
+H -0.812   0.469  -0.3816

tests/test_hf.py

@@ -53,6 +53,36 @@ def test_hartree_fock_ch4(self):
 
         en = -39.35007843284954
         np.testing.assert_almost_equal(results['energies'][-1], en, 4)
+
+    def test_hartree_fock_ch4_symmetric(self):
+        """
+        Test Hartree-Fock calculation on CH4 using an STO-3g basis set and
+        symmetric orthogonalization
+        """
+        mol = Molecule()
+        dist = 1.78/2
+        mol.add_atom('C', 0.0, 0.0, 0.0, unit='angstrom')
+        mol.add_atom('H', dist, dist, dist, unit='angstrom')
+        mol.add_atom('H', -dist, -dist, dist, unit='angstrom')
+        mol.add_atom('H', -dist, dist, -dist, unit='angstrom')
+        mol.add_atom('H', dist, -dist, -dist, unit='angstrom')
+
+        results = HF().rhf(mol, 'sto3g', ortho='symmetric')
+
+        # check that orbital energies are correctly approximated
+        ans = np.array([-11.0707,
+                         -0.7392,
+                         -0.3752,
+                         -0.3752,
+                         -0.3752,
+                          0.2865,
+                          0.4092,
+                          0.4092,
+                          0.4092])
+        np.testing.assert_almost_equal(results['orbe'], ans, 4)
+
+        en = -39.35007843284954
+        np.testing.assert_almost_equal(results['energies'][-1], en, 4)
 
     def test_hartree_fock_restart(self):
         """