From f99ba79e195f62753a54b43ffca17722f1fe0ae8 Mon Sep 17 00:00:00 2001
From: Ksenia <briling@mail.ru>
Date: Tue, 24 Sep 2024 17:51:43 +0200
Subject: [PATCH 1/6] Spahm in field

---
 qstack/spahm/compute_spahm.py | 26 +++++++++++++++++---------
 tests/test_spahm.py           | 10 ++++++++++
 2 files changed, 27 insertions(+), 9 deletions(-)

diff --git a/qstack/spahm/compute_spahm.py b/qstack/spahm/compute_spahm.py
index 201a0aa..80f849a 100644
--- a/qstack/spahm/compute_spahm.py
+++ b/qstack/spahm/compute_spahm.py
@@ -1,32 +1,40 @@
+import numpy as np
 from pyscf import scf, grad
 from qstack.spahm.guesses import solveF, get_guess, get_occ, get_dm, eigenvalue_grad, get_guess_g
 
-def get_guess_orbitals(mol, guess, xc="pbe"):
+def get_guess_orbitals(mol, guess, xc="pbe", field=None):
     """ Compute the guess Hamiltonian.
 
     Args:
         mol (pyscf Mole): pyscf Mole object.
         guess (funct): Method used to compute the guess Hamiltonian. Output of get_guess.
         xc (str): Exchange-correlation functional. Defaults to pbe.
+        field (numpy.array(3)): External electric field i.e. $\\vec \\nabla \\phi$
 
     Returns:
         A 1D numpy array containing the eigenvalues and a 2D numpy array containing the eigenvectors of the guess Hamiltonian.
     """
     if guess == 'huckel':
-        e,v = scf.hf._init_guess_huckel_orbitals(mol)
+        if field is not None:
+            raise NotImplementedError
+        e, v = scf.hf._init_guess_huckel_orbitals(mol)
     else:
         fock = guess(mol, xc)
-        e,v = solveF(mol, fock)
+        if field is not None:
+            with mol.with_common_orig((0,0,0)):
+                ao_dip = mol.intor_symmetric('int1e_r', comp=3)
+            fock += np.einsum('xij,x->ij', ao_dip, field)
+        e, v = solveF(mol, fock)
     return e,v
 
-def get_guess_orbitals_grad(mol, guess):
-    e, c = get_guess_orbitals(mol, guess[0])
+def get_guess_orbitals_grad(mol, guess, field=None):
+    e, c = get_guess_orbitals(mol, guess[0], field=field)
     mf = grad.rhf.Gradients(scf.RHF(mol))
     s1 = mf.get_ovlp(mol)
     h1 = guess[1](mf)
     return eigenvalue_grad(mol, e, c, s1, h1)
 
-def get_guess_dm(mol, guess, xc="pbe", openshell=None):
+def get_guess_dm(mol, guess, xc="pbe", openshell=None, field=None):
     """ Compute the density matrix with the guess Hamiltonian.
 
     Args:
@@ -38,10 +46,10 @@ def get_guess_dm(mol, guess, xc="pbe", openshell=None):
     Returns:
         A numpy ndarray containing the density matrix computed using the guess Hamiltonian.
     """
-    e,v = get_guess_orbitals(mol, guess, xc)
+    e,v = get_guess_orbitals(mol, guess, xc, field=field)
     return get_dm(v, mol.nelec, mol.spin if mol.spin>0 or not openshell is None else None)
 
-def get_spahm_representation(mol, guess_in, xc="pbe"):
+def get_spahm_representation(mol, guess_in, xc="pbe", field=None):
     """ Compute the SPAHM representation.
 
     Args:
@@ -53,7 +61,7 @@ def get_spahm_representation(mol, guess_in, xc="pbe"):
         A numpy ndarray containing the SPAHM representation.
     """
     guess = get_guess(guess_in)
-    e, v  = get_guess_orbitals(mol, guess, xc)
+    e, v  = get_guess_orbitals(mol, guess, xc, field=field)
     e1    = get_occ(e, mol.nelec, mol.spin)
     return e1
 
diff --git a/tests/test_spahm.py b/tests/test_spahm.py
index 0dfb0af..0f07445 100755
--- a/tests/test_spahm.py
+++ b/tests/test_spahm.py
@@ -47,7 +47,17 @@ def test_spahm_LB_ecp():
     assert np.allclose(R2, true_R2)
 
 
+def test_spahm_LB_field():
+    path = os.path.dirname(os.path.realpath(__file__))
+    mol = compound.xyz_to_mol(path+'/data/H2O.xyz', 'def2svp')
+    R = compute_spahm.get_spahm_representation(mol, 'lb', field=(0.01,0.01,0.01))
+    true_R = np.array([-18.26790464,  -0.7890498,   -0.32432933,  -0.17412611,  -0.10335613])
+    assert np.allclose(R[0], R[1])
+    assert np.allclose(true_R, R[0])
+
+
 if __name__ == '__main__':
     test_spahm_huckel()
     test_spahm_LB()
     test_spahm_LB_ecp()
+    test_spahm_LB_field()

From d42767082ea4d5081fb718d266111ed2faf60c37 Mon Sep 17 00:00:00 2001
From: Ksenia <briling@mail.ru>
Date: Fri, 27 Sep 2024 14:18:41 +0200
Subject: [PATCH 2/6] Compute Heff eigenvalues derivatives in a field

---
 qstack/spahm/compute_spahm.py | 28 +++++++++++++++++++++++++++-
 tests/test_spahm_grad.py      | 16 ++++++++++++++++
 2 files changed, 43 insertions(+), 1 deletion(-)

diff --git a/qstack/spahm/compute_spahm.py b/qstack/spahm/compute_spahm.py
index 80f849a..5fcbeae 100644
--- a/qstack/spahm/compute_spahm.py
+++ b/qstack/spahm/compute_spahm.py
@@ -27,13 +27,39 @@ def get_guess_orbitals(mol, guess, xc="pbe", field=None):
         e, v = solveF(mol, fock)
     return e,v
 
+
+def ext_field_generator(mol, field):
+    shls_slice = (0, mol.nbas, 0, mol.nbas)
+    with mol.with_common_orig((0,0,0)):
+        int1e_irp = mol.intor('int1e_irp', shls_slice=shls_slice).reshape(3, 3, mol.nao, mol.nao) # ( | rc nabla | )
+    aoslices = mol.aoslice_by_atom()[:,2:]
+    if field is None:
+        field = (0,0,0)
+    def field_deriv(iat):
+        p0, p1 = aoslices[iat]
+        dmu_dr = np.zeros_like(int1e_irp)  # dim(mu)×dim(r)×nao×nao
+        dmu_dr[:,:,p0:p1,:] -= int1e_irp[:,:,:,p0:p1].transpose((0,1,3,2))  # TODO not sure why minus
+        dmu_dr[:,:,:,p0:p1] -= int1e_irp[:,:,:,p0:p1]  # TODO check/fix E definition
+        dhext_dr = np.einsum('x,xypq->ypq', field, dmu_dr)
+        return dhext_dr
+    return field_deriv
+
+
 def get_guess_orbitals_grad(mol, guess, field=None):
     e, c = get_guess_orbitals(mol, guess[0], field=field)
     mf = grad.rhf.Gradients(scf.RHF(mol))
     s1 = mf.get_ovlp(mol)
-    h1 = guess[1](mf)
+    h0 = guess[1](mf)
+
+    if field is None:
+        h1 = h0
+    else:
+        hext = ext_field_generator(mf.mol, field)
+        h1 = lambda iat: h0(iat) + hext(iat)
+
     return eigenvalue_grad(mol, e, c, s1, h1)
 
+
 def get_guess_dm(mol, guess, xc="pbe", openshell=None, field=None):
     """ Compute the density matrix with the guess Hamiltonian.
 
diff --git a/tests/test_spahm_grad.py b/tests/test_spahm_grad.py
index 51bb2bb..78dd7b5 100755
--- a/tests/test_spahm_grad.py
+++ b/tests/test_spahm_grad.py
@@ -75,7 +75,23 @@ def spahm_ev(r, mol, guess):
         assert(np.linalg.norm(g1-g2)<1e-6)
 
 
+def test_spahm_ev_grad_field():
+    def spahm_ev(r, mol, guess):
+        mymol = build_mol(mol, r)
+        e, c = spahm.compute_spahm.get_guess_orbitals(mymol, guess[0], field=field)
+        return e
+    path  = os.path.dirname(os.path.realpath(__file__))
+    mol   = compound.xyz_to_mol(path+'/data/H2O_dist.xyz', 'def2svp', charge=0, spin=0)
+    field = np.array((0.01, 0.01, 0.01))
+    guess = spahm.guesses.get_guess_g('lb')
+    agrad = spahm.compute_spahm.get_guess_orbitals_grad(mol, guess, field=field)
+    ngrad = grad_num(spahm_ev, mol, guess)
+    for g1, g2 in zip(ngrad.T, agrad.reshape(-1, mol.natm*3)):
+        assert(np.linalg.norm(g1-g2)<1e-6)
+
+
 if __name__ == '__main__':
     test_spahm_ev_grad()
     test_spahm_re_grad()
     test_spahm_ev_grad_ecp()
+    test_spahm_ev_grad_field()

From d1e77a6be9d643bdaf6ad955608fcec597653d51 Mon Sep 17 00:00:00 2001
From: Ksenia <briling@mail.ru>
Date: Fri, 27 Sep 2024 14:22:21 +0200
Subject: [PATCH 3/6] Compute SPAHM representation derivatives in a field

---
 qstack/spahm/compute_spahm.py |  4 ++--
 tests/test_spahm_grad.py      | 16 ++++++++++++++++
 2 files changed, 18 insertions(+), 2 deletions(-)

diff --git a/qstack/spahm/compute_spahm.py b/qstack/spahm/compute_spahm.py
index 5fcbeae..adcd37c 100644
--- a/qstack/spahm/compute_spahm.py
+++ b/qstack/spahm/compute_spahm.py
@@ -91,7 +91,7 @@ def get_spahm_representation(mol, guess_in, xc="pbe", field=None):
     e1    = get_occ(e, mol.nelec, mol.spin)
     return e1
 
-def get_spahm_representation_grad(mol, guess_in):
+def get_spahm_representation_grad(mol, guess_in, field=None):
     guess = get_guess_g(guess_in)
-    agrad = get_guess_orbitals_grad(mol, guess)
+    agrad = get_guess_orbitals_grad(mol, guess, field=field)
     return get_occ(agrad, mol.nelec, mol.spin)
diff --git a/tests/test_spahm_grad.py b/tests/test_spahm_grad.py
index 78dd7b5..8014697 100755
--- a/tests/test_spahm_grad.py
+++ b/tests/test_spahm_grad.py
@@ -90,8 +90,24 @@ def spahm_ev(r, mol, guess):
         assert(np.linalg.norm(g1-g2)<1e-6)
 
 
+def test_spahm_re_grad_field():
+    def spahm_re(r, mol, guess_in):
+        mymol = build_mol(mol, r)
+        e = spahm.compute_spahm.get_spahm_representation(mymol, guess_in, field=field)
+        return e
+    path  = os.path.dirname(os.path.realpath(__file__))
+    mol   = compound.xyz_to_mol(path+'/data/H2O_dist.xyz', 'def2svp', charge=1, spin=1)
+    field = np.array((0.01, 0.01, 0.01))
+    guess = 'lb'
+    agrad = spahm.compute_spahm.get_spahm_representation_grad(mol, guess, field=field)
+    ngrad = grad_num(spahm_re, mol, guess)
+    for g1, g2 in zip(ngrad.reshape(mol.natm*3, -1).T, agrad.reshape(-1, mol.natm*3)):
+        assert(np.linalg.norm(g1-g2)<1e-6)
+
+
 if __name__ == '__main__':
     test_spahm_ev_grad()
     test_spahm_re_grad()
     test_spahm_ev_grad_ecp()
     test_spahm_ev_grad_field()
+    test_spahm_re_grad_field()

From 678c0e7fb502f49b66167463f4141b036ef4edde Mon Sep 17 00:00:00 2001
From: Ksenia <briling@mail.ru>
Date: Fri, 27 Sep 2024 14:23:54 +0200
Subject: [PATCH 4/6] Change H2O_dist to H2O_dist_rot for gradient tests

More nontrivial components
---
 tests/data/H2O_dist_rot.xyz | 5 +++++
 tests/test_spahm_grad.py    | 8 ++++----
 2 files changed, 9 insertions(+), 4 deletions(-)
 create mode 100644 tests/data/H2O_dist_rot.xyz

diff --git a/tests/data/H2O_dist_rot.xyz b/tests/data/H2O_dist_rot.xyz
new file mode 100644
index 0000000..53ddde3
--- /dev/null
+++ b/tests/data/H2O_dist_rot.xyz
@@ -0,0 +1,5 @@
+3
+
+ O      0.187252   -0.055745   -0.112535
+ H     -0.637721   -0.524458   -0.381093
+ H      0.450469    0.580203    0.493628
diff --git a/tests/test_spahm_grad.py b/tests/test_spahm_grad.py
index 8014697..3cba716 100755
--- a/tests/test_spahm_grad.py
+++ b/tests/test_spahm_grad.py
@@ -39,7 +39,7 @@ def spahm_ev(r, mol, guess):
         e, c = spahm.compute_spahm.get_guess_orbitals(mymol, guess[0])
         return e
     path  = os.path.dirname(os.path.realpath(__file__))
-    mol   = compound.xyz_to_mol(path+'/data/H2O_dist.xyz', 'def2svp', charge=0, spin=0)
+    mol   = compound.xyz_to_mol(path+'/data/H2O_dist_rot.xyz', 'def2svp', charge=0, spin=0)
     guess = spahm.guesses.get_guess_g('lb')
     agrad = spahm.compute_spahm.get_guess_orbitals_grad(mol, guess)
     ngrad = grad_num(spahm_ev, mol, guess)
@@ -53,7 +53,7 @@ def spahm_re(r, mol, guess_in):
         e = spahm.compute_spahm.get_spahm_representation(mymol, guess_in)
         return e
     path  = os.path.dirname(os.path.realpath(__file__))
-    mol   = compound.xyz_to_mol(path+'/data/H2O_dist.xyz', 'def2svp', charge=1, spin=1)
+    mol   = compound.xyz_to_mol(path+'/data/H2O_dist_rot.xyz', 'def2svp', charge=1, spin=1)
     guess = 'lb'
     agrad = spahm.compute_spahm.get_spahm_representation_grad(mol, guess)
     ngrad = grad_num(spahm_re, mol, guess)
@@ -81,7 +81,7 @@ def spahm_ev(r, mol, guess):
         e, c = spahm.compute_spahm.get_guess_orbitals(mymol, guess[0], field=field)
         return e
     path  = os.path.dirname(os.path.realpath(__file__))
-    mol   = compound.xyz_to_mol(path+'/data/H2O_dist.xyz', 'def2svp', charge=0, spin=0)
+    mol   = compound.xyz_to_mol(path+'/data/H2O_dist_rot.xyz', 'def2svp', charge=0, spin=0)
     field = np.array((0.01, 0.01, 0.01))
     guess = spahm.guesses.get_guess_g('lb')
     agrad = spahm.compute_spahm.get_guess_orbitals_grad(mol, guess, field=field)
@@ -96,7 +96,7 @@ def spahm_re(r, mol, guess_in):
         e = spahm.compute_spahm.get_spahm_representation(mymol, guess_in, field=field)
         return e
     path  = os.path.dirname(os.path.realpath(__file__))
-    mol   = compound.xyz_to_mol(path+'/data/H2O_dist.xyz', 'def2svp', charge=1, spin=1)
+    mol   = compound.xyz_to_mol(path+'/data/H2O_dist_rot.xyz', 'def2svp', charge=1, spin=1)
     field = np.array((0.01, 0.01, 0.01))
     guess = 'lb'
     agrad = spahm.compute_spahm.get_spahm_representation_grad(mol, guess, field=field)

From 49babe1d722c29e9d3fdddcdfd894b520a3e94dc Mon Sep 17 00:00:00 2001
From: Ksenia <briling@mail.ru>
Date: Fri, 27 Sep 2024 14:50:41 +0200
Subject: [PATCH 5/6] Return repr along with the grad; update docs

---
 qstack/spahm/compute_spahm.py | 52 ++++++++++++++++++++++++++++++-----
 tests/test_spahm_grad.py      | 30 ++++++++++----------
 2 files changed, 60 insertions(+), 22 deletions(-)

diff --git a/qstack/spahm/compute_spahm.py b/qstack/spahm/compute_spahm.py
index adcd37c..4aa3224 100644
--- a/qstack/spahm/compute_spahm.py
+++ b/qstack/spahm/compute_spahm.py
@@ -2,12 +2,13 @@
 from pyscf import scf, grad
 from qstack.spahm.guesses import solveF, get_guess, get_occ, get_dm, eigenvalue_grad, get_guess_g
 
+
 def get_guess_orbitals(mol, guess, xc="pbe", field=None):
-    """ Compute the guess Hamiltonian.
+    """ Compute the guess Hamiltonian orbitals
 
     Args:
         mol (pyscf Mole): pyscf Mole object.
-        guess (funct): Method used to compute the guess Hamiltonian. Output of get_guess.
+        guess (func): Method used to compute the guess Hamiltonian. Output of get_guess.
         xc (str): Exchange-correlation functional. Defaults to pbe.
         field (numpy.array(3)): External electric field i.e. $\\vec \\nabla \\phi$
 
@@ -25,10 +26,20 @@ def get_guess_orbitals(mol, guess, xc="pbe", field=None):
                 ao_dip = mol.intor_symmetric('int1e_r', comp=3)
             fock += np.einsum('xij,x->ij', ao_dip, field)
         e, v = solveF(mol, fock)
-    return e,v
+    return e, v
 
 
 def ext_field_generator(mol, field):
+    """ Generator for Hext (i.e. applied uniform electiric field interaction) gradient
+
+    Args:
+        mol (pyscf Mole): pyscf Mole object.
+        field (numpy.array(3)): External electric field i.e. $\\vec \\nabla \\phi$
+
+    Returns:
+        func(int: iat): returns the derivative of Hext wrt the coordinates of atom iat, i.e. dHext/dr[iat]
+    """
+
     shls_slice = (0, mol.nbas, 0, mol.nbas)
     with mol.with_common_orig((0,0,0)):
         int1e_irp = mol.intor('int1e_irp', shls_slice=shls_slice).reshape(3, 3, mol.nao, mol.nao) # ( | rc nabla | )
@@ -46,6 +57,18 @@ def field_deriv(iat):
 
 
 def get_guess_orbitals_grad(mol, guess, field=None):
+    """ Compute the guess Hamiltonian eigenvalues and their derivatives
+
+    Args:
+        mol (pyscf Mole): pyscf Mole object.
+        guess (func): Tuple of methods used to compute the guess Hamiltonian and its eigenvalue derivatives. Output of get_guess_g
+        field (numpy.array(3)): External electric field i.e. $\\vec \\nabla \\phi$
+
+    Returns:
+        numpy 1d array (mol.nao,): eigenvalues
+        numpy 3d ndarray (mol.nao,mol.natm,3): gradient of the eigenvalues in Eh/bohr
+    """
+
     e, c = get_guess_orbitals(mol, guess[0], field=field)
     mf = grad.rhf.Gradients(scf.RHF(mol))
     s1 = mf.get_ovlp(mol)
@@ -57,7 +80,7 @@ def get_guess_orbitals_grad(mol, guess, field=None):
         hext = ext_field_generator(mf.mol, field)
         h1 = lambda iat: h0(iat) + hext(iat)
 
-    return eigenvalue_grad(mol, e, c, s1, h1)
+    return e, eigenvalue_grad(mol, e, c, s1, h1)
 
 
 def get_guess_dm(mol, guess, xc="pbe", openshell=None, field=None):
@@ -65,7 +88,7 @@ def get_guess_dm(mol, guess, xc="pbe", openshell=None, field=None):
 
     Args:
         mol (pyscf Mole): pyscf Mole object.
-        guess (funct): Method used to compute the guess Hamiltonian. Output of get_guess.
+        guess (func): Method used to compute the guess Hamiltonian. Output of get_guess.
         xc (str): Exchange-correlation functional. Defaults to pbe
         openshell (bool): . Defaults to None.
 
@@ -75,6 +98,7 @@ def get_guess_dm(mol, guess, xc="pbe", openshell=None, field=None):
     e,v = get_guess_orbitals(mol, guess, xc, field=field)
     return get_dm(v, mol.nelec, mol.spin if mol.spin>0 or not openshell is None else None)
 
+
 def get_spahm_representation(mol, guess_in, xc="pbe", field=None):
     """ Compute the SPAHM representation.
 
@@ -82,6 +106,7 @@ def get_spahm_representation(mol, guess_in, xc="pbe", field=None):
         mol (pyscf Mole): pyscf Mole object.
         guess_in (str): Method used to obtain the guess Hamiltoninan.
         xc (str): Exchange-correlation functional. Defaults to pbe.
+        field (numpy.array(3)): External electric field i.e. $\\vec \\nabla \\phi$
 
     Returns:
         A numpy ndarray containing the SPAHM representation.
@@ -91,7 +116,20 @@ def get_spahm_representation(mol, guess_in, xc="pbe", field=None):
     e1    = get_occ(e, mol.nelec, mol.spin)
     return e1
 
+
 def get_spahm_representation_grad(mol, guess_in, field=None):
+    """ Compute the SPAHM representation and its gradient
+
+    Args:
+        mol (pyscf Mole): pyscf Mole object.
+        guess_in (str): Method used to obtain the guess Hamiltoninan.
+        xc (str): Exchange-correlation functional. Defaults to pbe.
+        field (numpy.array(3)): External electric field i.e. $\\vec \\nabla \\phi$
+
+    Returns:
+        numpy 1d array (occ,): the SPAHM representation (Eh).
+        numpy 3d array (occ,mol.natm,3): gradient of the representation (Eh/bohr)
+    """
     guess = get_guess_g(guess_in)
-    agrad = get_guess_orbitals_grad(mol, guess, field=field)
-    return get_occ(agrad, mol.nelec, mol.spin)
+    e, agrad = get_guess_orbitals_grad(mol, guess, field=field)
+    return get_occ(e, mol.nelec, mol.spin), get_occ(agrad, mol.nelec, mol.spin)
diff --git a/tests/test_spahm_grad.py b/tests/test_spahm_grad.py
index 3cba716..82c6777 100755
--- a/tests/test_spahm_grad.py
+++ b/tests/test_spahm_grad.py
@@ -41,9 +41,9 @@ def spahm_ev(r, mol, guess):
     path  = os.path.dirname(os.path.realpath(__file__))
     mol   = compound.xyz_to_mol(path+'/data/H2O_dist_rot.xyz', 'def2svp', charge=0, spin=0)
     guess = spahm.guesses.get_guess_g('lb')
-    agrad = spahm.compute_spahm.get_guess_orbitals_grad(mol, guess)
-    ngrad = grad_num(spahm_ev, mol, guess)
-    for g1, g2 in zip(ngrad.T, agrad.reshape(-1, 9)):
+    agrad = spahm.compute_spahm.get_guess_orbitals_grad(mol, guess)[1].reshape(-1, mol.natm*3)
+    ngrad = grad_num(spahm_ev, mol, guess).T
+    for g1, g2 in zip(ngrad, agrad):
         assert(np.linalg.norm(g1-g2)<1e-6)
 
 
@@ -55,9 +55,9 @@ def spahm_re(r, mol, guess_in):
     path  = os.path.dirname(os.path.realpath(__file__))
     mol   = compound.xyz_to_mol(path+'/data/H2O_dist_rot.xyz', 'def2svp', charge=1, spin=1)
     guess = 'lb'
-    agrad = spahm.compute_spahm.get_spahm_representation_grad(mol, guess)
-    ngrad = grad_num(spahm_re, mol, guess)
-    for g1, g2 in zip(ngrad.reshape(9, -1).T, agrad.reshape(-1, 9)):
+    agrad = spahm.compute_spahm.get_spahm_representation_grad(mol, guess)[1].reshape(-1, mol.natm*3)
+    ngrad = grad_num(spahm_re, mol, guess).reshape(mol.natm*3, -1).T
+    for g1, g2 in zip(ngrad, agrad):
         assert(np.linalg.norm(g1-g2)<1e-6)
 
 
@@ -69,9 +69,9 @@ def spahm_ev(r, mol, guess):
     path  = os.path.dirname(os.path.realpath(__file__))
     mol   = compound.xyz_to_mol(path+'/data/H2Te.xyz', 'minao', charge=0, spin=0, ecp='def2-svp')
     guess = spahm.guesses.get_guess_g('lb')
-    agrad = spahm.compute_spahm.get_guess_orbitals_grad(mol, guess)
-    ngrad = grad_num(spahm_ev, mol, guess)
-    for g1, g2 in zip(ngrad.T, agrad.reshape(-1, 9)):
+    agrad = spahm.compute_spahm.get_guess_orbitals_grad(mol, guess)[1].reshape(-1, mol.natm*3)
+    ngrad = grad_num(spahm_ev, mol, guess).T
+    for g1, g2 in zip(ngrad, agrad):
         assert(np.linalg.norm(g1-g2)<1e-6)
 
 
@@ -84,9 +84,9 @@ def spahm_ev(r, mol, guess):
     mol   = compound.xyz_to_mol(path+'/data/H2O_dist_rot.xyz', 'def2svp', charge=0, spin=0)
     field = np.array((0.01, 0.01, 0.01))
     guess = spahm.guesses.get_guess_g('lb')
-    agrad = spahm.compute_spahm.get_guess_orbitals_grad(mol, guess, field=field)
-    ngrad = grad_num(spahm_ev, mol, guess)
-    for g1, g2 in zip(ngrad.T, agrad.reshape(-1, mol.natm*3)):
+    agrad = spahm.compute_spahm.get_guess_orbitals_grad(mol, guess, field=field)[1].reshape(-1, mol.natm*3)
+    ngrad = grad_num(spahm_ev, mol, guess).T
+    for g1, g2 in zip(ngrad, agrad):
         assert(np.linalg.norm(g1-g2)<1e-6)
 
 
@@ -99,9 +99,9 @@ def spahm_re(r, mol, guess_in):
     mol   = compound.xyz_to_mol(path+'/data/H2O_dist_rot.xyz', 'def2svp', charge=1, spin=1)
     field = np.array((0.01, 0.01, 0.01))
     guess = 'lb'
-    agrad = spahm.compute_spahm.get_spahm_representation_grad(mol, guess, field=field)
-    ngrad = grad_num(spahm_re, mol, guess)
-    for g1, g2 in zip(ngrad.reshape(mol.natm*3, -1).T, agrad.reshape(-1, mol.natm*3)):
+    agrad = spahm.compute_spahm.get_spahm_representation_grad(mol, guess, field=field)[1].reshape(-1, mol.natm*3)
+    ngrad = grad_num(spahm_re, mol, guess).reshape(mol.natm*3, -1).T
+    for g1, g2 in zip(ngrad, agrad):
         assert(np.linalg.norm(g1-g2)<1e-6)
 
 

From ed5cc6c093b12f3131616ac050e0f77293b0e701 Mon Sep 17 00:00:00 2001
From: Ksenia <briling@mail.ru>
Date: Fri, 27 Sep 2024 15:04:55 +0200
Subject: [PATCH 6/6] Clarify that field is in atomic units

---
 qstack/spahm/compute_spahm.py | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/qstack/spahm/compute_spahm.py b/qstack/spahm/compute_spahm.py
index 4aa3224..74d6999 100644
--- a/qstack/spahm/compute_spahm.py
+++ b/qstack/spahm/compute_spahm.py
@@ -10,7 +10,7 @@ def get_guess_orbitals(mol, guess, xc="pbe", field=None):
         mol (pyscf Mole): pyscf Mole object.
         guess (func): Method used to compute the guess Hamiltonian. Output of get_guess.
         xc (str): Exchange-correlation functional. Defaults to pbe.
-        field (numpy.array(3)): External electric field i.e. $\\vec \\nabla \\phi$
+        field (numpy.array(3)): applied uniform electric field i.e. $\\vec \\nabla \\phi$ in a.u.
 
     Returns:
         A 1D numpy array containing the eigenvalues and a 2D numpy array containing the eigenvectors of the guess Hamiltonian.
@@ -34,7 +34,7 @@ def ext_field_generator(mol, field):
 
     Args:
         mol (pyscf Mole): pyscf Mole object.
-        field (numpy.array(3)): External electric field i.e. $\\vec \\nabla \\phi$
+        field (numpy.array(3)): applied uniform electric field i.e. $\\vec \\nabla \\phi$ in a.u.
 
     Returns:
         func(int: iat): returns the derivative of Hext wrt the coordinates of atom iat, i.e. dHext/dr[iat]
@@ -62,7 +62,7 @@ def get_guess_orbitals_grad(mol, guess, field=None):
     Args:
         mol (pyscf Mole): pyscf Mole object.
         guess (func): Tuple of methods used to compute the guess Hamiltonian and its eigenvalue derivatives. Output of get_guess_g
-        field (numpy.array(3)): External electric field i.e. $\\vec \\nabla \\phi$
+        field (numpy.array(3)): applied uniform electric field i.e. $\\vec \\nabla \\phi$ in a.u.
 
     Returns:
         numpy 1d array (mol.nao,): eigenvalues
@@ -106,7 +106,7 @@ def get_spahm_representation(mol, guess_in, xc="pbe", field=None):
         mol (pyscf Mole): pyscf Mole object.
         guess_in (str): Method used to obtain the guess Hamiltoninan.
         xc (str): Exchange-correlation functional. Defaults to pbe.
-        field (numpy.array(3)): External electric field i.e. $\\vec \\nabla \\phi$
+        field (numpy.array(3)): applied uniform electric field i.e. $\\vec \\nabla \\phi$ in a.u.
 
     Returns:
         A numpy ndarray containing the SPAHM representation.
@@ -124,7 +124,7 @@ def get_spahm_representation_grad(mol, guess_in, field=None):
         mol (pyscf Mole): pyscf Mole object.
         guess_in (str): Method used to obtain the guess Hamiltoninan.
         xc (str): Exchange-correlation functional. Defaults to pbe.
-        field (numpy.array(3)): External electric field i.e. $\\vec \\nabla \\phi$
+        field (numpy.array(3)): applied uniform electric field i.e. $\\vec \\nabla \\phi$ in a.u.
 
     Returns:
         numpy 1d array (occ,): the SPAHM representation (Eh).