diff --git a/.github/workflows/python-package-conda-cache.yml b/.github/workflows/python-package-conda-cache.yml
index 4c3b3bf..05039de 100644
--- a/.github/workflows/python-package-conda-cache.yml
+++ b/.github/workflows/python-package-conda-cache.yml
@@ -47,7 +47,7 @@ jobs:
       - name: Set cache date
         run: echo "DATE=$(date +'%Y%m%d')" >> $GITHUB_ENV
 
-      - uses: actions/cache@v2
+      - uses: actions/cache@v4
         with:
           path: ${{ matrix.prefix }}
           key: ${{ matrix.label }}-conda-${{ hashFiles('environment.yml') }}-${{ env.DATE }}-${{ env.CACHE_NUMBER }}
diff --git a/qstack/spahm/compute_spahm.py b/qstack/spahm/compute_spahm.py
index 201a0aa..74d6999 100644
--- a/qstack/spahm/compute_spahm.py
+++ b/qstack/spahm/compute_spahm.py
@@ -1,63 +1,135 @@
+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"):
-    """ Compute the guess Hamiltonian.
+
+def get_guess_orbitals(mol, guess, xc="pbe", field=None):
+    """ 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)): 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.
     """
     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)
-    return e,v
+        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 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)): 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]
+    """
+
+    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):
+    """ 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)): applied uniform electric field i.e. $\\vec \\nabla \\phi$ in a.u.
+
+    Returns:
+        numpy 1d array (mol.nao,): eigenvalues
+        numpy 3d ndarray (mol.nao,mol.natm,3): gradient of the eigenvalues in Eh/bohr
+    """
 
-def get_guess_orbitals_grad(mol, guess):
-    e, c = get_guess_orbitals(mol, guess[0])
+    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)
+    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)
 
-def get_guess_dm(mol, guess, xc="pbe", openshell=None):
+    return e, 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.
 
     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.
 
     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:
         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)): applied uniform electric field i.e. $\\vec \\nabla \\phi$ in a.u.
 
     Returns:
         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
 
-def get_spahm_representation_grad(mol, guess_in):
+
+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)): applied uniform electric field i.e. $\\vec \\nabla \\phi$ in a.u.
+
+    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)
-    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/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.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()
diff --git a/tests/test_spahm_grad.py b/tests/test_spahm_grad.py
index 51bb2bb..82c6777 100755
--- a/tests/test_spahm_grad.py
+++ b/tests/test_spahm_grad.py
@@ -39,11 +39,11 @@ 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)
-    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)
 
 
@@ -53,11 +53,11 @@ 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)
-    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,39 @@ 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)
+
+
+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_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)[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)
+
+
+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_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)[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)
 
 
@@ -79,3 +109,5 @@ def spahm_ev(r, mol, guess):
     test_spahm_ev_grad()
     test_spahm_re_grad()
     test_spahm_ev_grad_ecp()
+    test_spahm_ev_grad_field()
+    test_spahm_re_grad_field()