diff --git a/pyscf/neo/efield.py b/pyscf/neo/efield.py
index 0600bcb09f..0b3d3aa4b2 100644
--- a/pyscf/neo/efield.py
+++ b/pyscf/neo/efield.py
@@ -3,10 +3,12 @@
 CNEO with electric field
 '''
 import numpy
+from pyscf import lib
+from pyscf.lib import logger
 from pyscf.neo import cphf, CDFT
 from pyscf.neo.grad import Gradients
 from pyscf.neo.hessian import gen_vind
-from pyscf import lib
+from pyscf.neo.ks import KS
 
 def solve_mo1(mf, mo_energy, mo_coeff, mo_occ, h1ao=None,
               fx=None, atmlst=None, max_memory=4000, verbose=None,
@@ -41,31 +43,60 @@ def solve_mo1(mf, mo_energy, mo_coeff, mo_occ, h1ao=None,
 
     return mo1, e1
 
-def polarizability(mf):
-    'polarizability in CNEO'
-
+def polarizability(polobj, with_cphf=True):
+    log = logger.new_logger(polobj)
+    mf = polobj._scf
+    with_f1 = polobj.with_f1
     mol = mf.mol
-
-    mo1 = solve_mo1(mf, mf.mo_energy, mf.mo_coeff, mf.mo_occ)[0]
+    mo_energy = mf.mo_energy
+    mo_coeff = mf.mo_coeff
+    mo_occ = mf.mo_occ
 
     charges = mol.atom_charges()
     coords  = mol.atom_coords()
     charge_center = numpy.einsum('i,ix->x', charges, coords) / charges.sum()
+    h1 = {}
+    s1 = {}
+    for t, comp in mf.components.items():
+        with comp.mol.with_common_orig(charge_center):
+            int1e_r = comp.mol.intor_symmetric('int1e_r', comp=3)
+        occidx = mo_occ[t] > 0
+        orbo = mo_coeff[t][:, occidx]
+        h1[t] = lib.einsum('xpq,pi,qj->xij', int1e_r, mo_coeff[t], orbo) * comp.charge
+        s1[t] = numpy.zeros_like(h1[t])
+    vind = polobj.gen_vind(mf, mo_coeff, mo_occ)
+    if with_cphf:
+        mo1 = cphf.solve(vind, mo_energy, mo_occ, h1, s1, with_f1=with_f1,
+                         max_cycle=polobj.max_cycle_cphf, tol=polobj.conv_tol,
+                         verbose=log)[0]
+    else:
+        raise NotImplementedError('without cphf is not implemented yet')
+    e2 = 0
+    for t in mf.components.keys():
+        if t == 'e':
+            e2 += 2*numpy.einsum('xpi,ypi->xy', h1['e'], mo1['e'])  #*2 for double occupancy
+        else:
+            ### if with_f1 is False, the contribution from quantum nuclei is included
+            if not with_f1 and t.startswith('n'):
+                e2 += numpy.einsum('xpi,ypi->xy', h1[t], mo1[t])
+    # *-1 from the definition of dipole moment.
+    e2 = (e2 + e2.T) * -1
+
+    if mf.verbose >= logger.INFO:
+        xx, yy, zz = e2.diagonal()
+        log.note('Isotropic polarizability %.12g', (xx+yy+zz)/3)
+        log.note('Polarizability anisotropy %.12g',
+                 (.5 * ((xx-yy)**2 + (yy-zz)**2 + (zz-xx)**2))**.5)
+        log.debug('Static polarizability tensor\n%s', e2)
+    return e2
 
-    with mol.components['e'].with_common_orig(charge_center):
-        int_r = mol.components['e'].intor_symmetric('int1e_r', comp=3)
 
-    mo_coeff = mf.components['e'].mo_coeff
-    mo_occ = mf.components['e'].mo_occ
-    mocc = mo_coeff[:, mo_occ>0]
-    h1 = numpy.einsum('xpq,pi,qj->xij', int_r, mo_coeff, mocc)
+def hyper_polarizability(polobj, with_cphf=True):
+    return NotImplementedError('hyperpolarizability is not implemented yet')
 
-    e2 = numpy.einsum('xpi,ypi->xy', h1, mo1['e'])
 
-    # *-1 from the definition of dipole moment. *2 for double occupancy
-    e2 = (e2 + e2.T) * -2
-
-    return e2
+def polarizability_with_freq(polobj, freq, with_cphf=True):
+    return NotImplementedError('frequency-dependent polarizability is not implemented yet')
 
 
 def dipole_grad(mf):
@@ -201,6 +232,60 @@ def grad_nuc(self, atmlst=None):
 
 SCFwithEfield.Gradients = lib.class_as_method(GradwithEfield)
 
+class NEOwithEfield(KS):
+    '''NEO with electric field'''
+    _keys = {'efield'}
+
+    def __init__(self, mol, *args, **kwargs):
+        super().__init__(mol, *args, **kwargs)
+        self.efield = numpy.array([0, 0, 0])
+        self.mol = mol
+
+    def get_hcore(self, mol=None):
+        if mol is None: mol = self.mol
+        hcore = {}
+        for t, comp in mol.components.items():
+            hcore[t] = self.components[t].get_hcore(mol=comp)
+            comp.set_common_orig([0, 0, 0])  # The gauge origin for dipole integral
+            hcore[t] += (numpy.einsum('x,xij->ij', self.efield,
+                                      comp.intor('int1e_r', comp=3))
+                         * self.components[t].charge)
+
+        return hcore
+
+    def energy_nuc(self):
+        enuc = self.components['e'].energy_nuc()
+
+        nuclear_charges = self.mol.components['e'].atom_charges()
+        nuclear_coords = self.mol.atom_coords()
+
+        E_nuc_field = -numpy.sum([Z * numpy.dot(self.efield, R) for Z, R in zip(nuclear_charges, nuclear_coords)])
+
+        return enuc + E_nuc_field
+
+class Polarizability(lib.StreamObject):
+    def __init__(self, mf):
+        mol = mf.mol
+        self.mol = mol
+        self.verbose = mol.verbose
+        self.stdout = mol.stdout
+        self._scf = mf
+        if isinstance(mf, CDFT):
+            self.with_f1 = True
+        else:
+            self.with_f1 = False
+        self.max_cycle_cphf = 100
+        self.conv_tol = 1e-9
+
+        self._keys = set(self.__dict__.keys())
+
+    def gen_vind(self, mf, mo_coeff, mo_occ):
+        return gen_vind(mf, mo_coeff, mo_occ)
+
+    polarizability = polarizability
+    polarizability_with_freq = polarizability_with_freq
+    hyper_polarizability = hyper_polarizability
+
 if __name__ == '__main__':
     from pyscf import neo
     mol = neo.M(atom='H 0 0 0; F 0 0 0.8', basis='ccpvdz')
diff --git a/pyscf/neo/hf.py b/pyscf/neo/hf.py
index ac20331546..0bb2a5b992 100644
--- a/pyscf/neo/hf.py
+++ b/pyscf/neo/hf.py
@@ -493,12 +493,20 @@ def mulliken_meta(self, mol=None, dm=None, verbose=logger.DEBUG,
     def dip_moment(self, mol=None, dm=None, unit='Debye', origin=None, verbose=logger.NOTE,
                    **kwargs):
         if self.is_nucleus:
-            return None
+            if origin is None:
+                origin = numpy.zeros(3)
+            else:
+                origin = numpy.asarray(origin, dtype=numpy.float64)
+            assert origin.shape == (3,)
+            dip = -self.charge * (lib.einsum('xij,ji->x', self.mol.intor_symmetric('int1e_r', comp=3), dm) - origin)
+            if unit.upper() == 'DEBYE':
+                dip *= nist.AU2DEBYE
         # Temporarily modify charge to supress warning about nonzero charge for a neutral molecule
-        charge = self.mol.charge
-        self.mol.charge = self.mol.super_mol.charge
-        dip = super().dip_moment(mol, dm, unit, origin=origin, verbose=verbose, **kwargs)
-        self.mol.charge = charge
+        else:
+            charge = self.mol.charge
+            self.mol.charge = self.mol.super_mol.charge
+            dip = super().dip_moment(mol, dm, unit, origin=origin, verbose=verbose, **kwargs)
+            self.mol.charge = charge
         return dip
 
     def to_gpu(self):
@@ -1425,7 +1433,26 @@ def canonicalize(self, mo_coeff, mo_occ, fock=None):
 
     def dip_moment(self, mol=None, dm=None, unit='Debye', origin=None,
                    verbose=logger.NOTE, **kwargs):
-        raise TypeError('Use CDFT to calculate total dipole moment.') # CDFT will have this
+        if mol is None: mol = self.mol
+        if dm is None: dm = self.make_rdm1()
+        log = logger.new_logger(mol, verbose)
+
+        # Suppress warning about nonzero charge (if neutral)
+        charge = self.components['e'].mol.charge
+        self.components['e'].mol.charge = self.mol.charge
+        el_dip = self.components['e'].dip_moment(mol.components['e'],
+                                                 dm['e'], unit=unit,
+                                                 origin=origin, verbose=verbose-1)
+        self.components['e'].mol.charge = charge
+
+        # Quantum nuclei
+        nucl_dip = 0
+        for t, comp in self.components.items():
+            if t.startswith('n'):
+                nucl_dip += comp.dip_moment(mol.components[t], dm[t], unit=unit, origin=origin, verbose=verbose-1)
+        mol_dip = nucl_dip + el_dip
+        log.note('Dipole moment(X, Y, Z, Debye): %8.5f, %8.5f, %8.5f', *mol_dip)
+        return mol_dip
 
     def quad_moment(self, mol=None, dm=None, unit='DebyeAngstrom', origin=None,
                     verbose=logger.NOTE, **kwargs):
diff --git a/pyscf/neo/test/test_efield.py b/pyscf/neo/test/test_efield.py
index af9c2ad14e..00ee1bfc30 100644
--- a/pyscf/neo/test/test_efield.py
+++ b/pyscf/neo/test/test_efield.py
@@ -3,7 +3,7 @@
 import unittest
 import numpy
 from pyscf import neo, lib, scf
-from pyscf.neo.efield import polarizability, dipole_grad, SCFwithEfield, GradwithEfield
+from pyscf.neo.efield import Polarizability, dipole_grad, SCFwithEfield, GradwithEfield, NEOwithEfield
 
 
 class KnownValues(unittest.TestCase):
@@ -51,10 +51,11 @@ def test_dipole_grad(self):
 
 
     def test_polarizability(self):
-        mol = neo.M(atom='H 0 0 0; F 0 0 0.9', basis='ccpvdz')
-        mf = neo.CDFT(mol, xc='b3lyp')
-        mf.run()
-        p = polarizability(mf)
+        mol = neo.M(atom='H 0 0 0; F 0 0 0.9', basis='ccpvdz', quantum_nuc=['H'])
+        mfCNEO = neo.CDFT(mol, xc='b3lyp')
+        mfCNEO.run()
+        polobj = Polarizability(mfCNEO)
+        p = polobj.polarizability()
 
         mf1 = SCFwithEfield(mol, xc='b3lyp')
         mf1.efield = numpy.array([0, 0, 0.001])
@@ -66,7 +67,23 @@ def test_polarizability(self):
         mf2.scf()
         dipole2 = mf2.dip_moment(unit='au')
 
+        mfNEO = neo.KS(mol,xc='b3lyp')
+        mfNEO.run()
+        polobj1 = Polarizability(mfNEO)
+        p1 = polobj1.polarizability()
+
+        mf3 = NEOwithEfield(mol, xc='b3lyp')
+        mf3.efield = numpy.array([0, 0, 0.001])
+        mf3.scf()
+        dipole3 = mf3.dip_moment(unit='au')
+
+        mf4 = NEOwithEfield(mol, xc='b3lyp')
+        mf4.efield = numpy.array([0, 0, -0.001])
+        mf4.scf()
+        dipole4 = mf4.dip_moment(unit='au')
+
         self.assertAlmostEqual(p[-1,-1], (dipole1[-1] - dipole2[-1]) / 0.002, 4)
+        self.assertAlmostEqual(p1[-1,-1], (dipole3[-1] - dipole4[-1]) / 0.002, 4)
 
     def test_grad_with_efield(self):
         mol = neo.M(atom='H 0 0 0; F 0 0 0.8', basis='ccpvdz')