From e1592ba182e180b27f0d9151de2b6ec06a8b60ba Mon Sep 17 00:00:00 2001
From: ifilot <ivo@ivofilot.nl>
Date: Fri, 1 Dec 2023 20:39:32 +0100
Subject: [PATCH] Fixing ammonia molecule and adding symmetric
 orthogonalization

---
 meta.yaml                |  2 +-
 pyqint/_version.py       |  2 +-
 pyqint/hf.py             | 19 +++++++++++++------
 pyqint/molecules/nh3.xyz | 11 +++++------
 tests/test_hf.py         | 30 ++++++++++++++++++++++++++++++
 5 files changed, 50 insertions(+), 14 deletions(-)

diff --git a/meta.yaml b/meta.yaml
index 9ba977a..5152033 100644
--- a/meta.yaml
+++ b/meta.yaml
@@ -1,6 +1,6 @@
 package:
   name: "pyqint"
-  version: "0.14.0.1"
+  version: "0.14.1.0"
 
 source:
   path: .
diff --git a/pyqint/_version.py b/pyqint/_version.py
index 2b1964a..d417133 100644
--- a/pyqint/_version.py
+++ b/pyqint/_version.py
@@ -1,2 +1,2 @@
-__version__ = "0.14.0.1"
+__version__ = "0.14.1.0"
 
diff --git a/pyqint/hf.py b/pyqint/hf.py
index 9ac6c8b..f9a3249 100644
--- a/pyqint/hf.py
+++ b/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
diff --git a/pyqint/molecules/nh3.xyz b/pyqint/molecules/nh3.xyz
index 9f35c83..dd0ef52 100644
--- a/pyqint/molecules/nh3.xyz
+++ b/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
\ No newline at end of file
diff --git a/tests/test_hf.py b/tests/test_hf.py
index cc38650..90ed714 100644
--- a/tests/test_hf.py
+++ b/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):
         """