diff --git a/.github/workflows/ci_linux/python_deps.sh b/.github/workflows/ci_linux/python_deps.sh
index 03722ee456..3f52c7cb8a 100755
--- a/.github/workflows/ci_linux/python_deps.sh
+++ b/.github/workflows/ci_linux/python_deps.sh
@@ -7,6 +7,8 @@ version=$(python -c 'import sys; version=sys.version_info[:2]; print("{0}.{1}".f
 if [ $version != '3.12' ]; then
     pip install geometric
     pip install spglib
+    pip install dftd3
+    pip install dftd4
 fi
 
 #cppe
diff --git a/pyscf/dft/rks.py b/pyscf/dft/rks.py
index af9312b094..70f5fdea37 100644
--- a/pyscf/dft/rks.py
+++ b/pyscf/dft/rks.py
@@ -243,6 +243,7 @@ def energy_elec(ks, dm=None, h1e=None, vhf=None):
     ecoul = vhf.ecoul.real
     exc = vhf.exc.real
     e2 = ecoul + exc
+
     ks.scf_summary['e1'] = e1
     ks.scf_summary['coul'] = ecoul
     ks.scf_summary['exc'] = exc
diff --git a/pyscf/dft/test/test_h2o.py b/pyscf/dft/test/test_h2o.py
index 4e71c9bae5..7a147c9bae 100644
--- a/pyscf/dft/test/test_h2o.py
+++ b/pyscf/dft/test/test_h2o.py
@@ -19,6 +19,18 @@
 from pyscf import lib
 from pyscf import dft
 
+
+import sys
+try:
+    import dftd3
+except ImportError:
+    pass
+
+try:
+    import dftd4
+except ImportError:
+    pass
+
 def setUpModule():
     global h2o, h2osym, h2o_cation, h2osym_cation
     h2o = gto.Mole()
@@ -505,6 +517,14 @@ def test_camb3lyp_rsh_omega(self):
         mf2.kernel()
         self.assertAlmostEqual(mf1.e_tot, -76.36649222362115, 9)
 
+    @unittest.skipIf('dftd3' not in sys.modules, "requires the dftd3 library")
+    def test_dispersion(self):
+        mf = dft.RKS(h2o)
+        mf.xc = 'B3LYP'
+        mf.disp = 'd3bj'
+        mf.run(xc='B3LYP')
+        self.assertAlmostEqual(mf.e_tot, -76.38945547396322, 9)
+
     def test_reset(self):
         mf = dft.RKS(h2o).newton()
         mf.reset(h2osym)
diff --git a/pyscf/grad/__init__.py b/pyscf/grad/__init__.py
index dc36fa1697..4bd6a695c1 100644
--- a/pyscf/grad/__init__.py
+++ b/pyscf/grad/__init__.py
@@ -32,6 +32,7 @@
 from . import dhf
 from . import uhf
 from . import rohf
+from . import dispersion
 from .rhf import Gradients as RHF
 from .dhf import Gradients as DHF
 from .uhf import Gradients as UHF
diff --git a/pyscf/grad/dispersion.py b/pyscf/grad/dispersion.py
new file mode 100644
index 0000000000..b4865db0a8
--- /dev/null
+++ b/pyscf/grad/dispersion.py
@@ -0,0 +1,59 @@
+#!/usr/bin/env python
+# Copyright 2014-2023 The PySCF Developers. All Rights Reserved.
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+#
+# Author: Xiaojie Wu <wxj6000@gmail.com>
+#
+
+'''
+gradient of dispersion correction for HF and DFT
+'''
+
+import numpy
+from pyscf.scf.hf import KohnShamDFT
+
+def get_dispersion(mf_grad, disp_version=None):
+    '''gradient of dispersion correction for RHF/RKS'''
+    if disp_version is None:
+        disp_version = mf_grad.base.disp
+    mol = mf_grad.base.mol
+    disp_version = mf_grad.base.disp
+    if disp_version is None:
+        return numpy.zeros([mol.natm,3])
+
+    if isinstance(mf_grad.base, KohnShamDFT):
+        method = mf_grad.base.xc
+    else:
+        method = 'hf'
+
+    if disp_version[:2].upper() == 'D3':
+        # raised error in SCF module, assuming dftd3 installed
+        import dftd3.pyscf as disp
+        d3 = disp.DFTD3Dispersion(mol, xc=method, version=disp_version)
+        _, g_d3 = d3.kernel()
+        return g_d3
+    elif disp_version[:2].upper() == 'D4':
+        from pyscf.data.elements import charge
+        atoms = numpy.array([ charge(a[0]) for a in mol._atom])
+        coords = mol.atom_coords()
+        from dftd4.interface import DampingParam, DispersionModel
+        model = DispersionModel(atoms, coords)
+        res = model.get_dispersion(DampingParam(method=method), grad=True)
+        return res.get("gradient")
+    else:
+        raise RuntimeError(f'dispersion correction: {disp_version} is not supported.')
+
+# Inject to Gradient
+from pyscf import grad
+grad.rhf.GradientsBase.get_dispersion = get_dispersion
diff --git a/pyscf/grad/rhf.py b/pyscf/grad/rhf.py
index 543890f5dd..e5eaf7f7f8 100644
--- a/pyscf/grad/rhf.py
+++ b/pyscf/grad/rhf.py
@@ -416,6 +416,8 @@ def kernel(self, mo_energy=None, mo_coeff=None, mo_occ=None, atmlst=None):
         self.de = de + self.grad_nuc(atmlst=atmlst)
         if self.mol.symmetry:
             self.de = self.symmetrize(self.de, atmlst)
+        if self.base.disp is not None:
+            self.de += self.get_dispersion()
         logger.timer(self, 'SCF gradients', *cput0)
         self._finalize()
         return self.de
diff --git a/pyscf/grad/test/test_rhf.py b/pyscf/grad/test/test_rhf.py
index b70679ab3a..263be449b5 100644
--- a/pyscf/grad/test/test_rhf.py
+++ b/pyscf/grad/test/test_rhf.py
@@ -18,6 +18,17 @@
 from pyscf import gto, scf, lib
 from pyscf import grad
 
+import sys
+try:
+    import dftd3
+except ImportError:
+    pass
+
+try:
+    import dftd4
+except ImportError:
+    pass
+
 def setUpModule():
     global mol
     mol = gto.Mole()
@@ -71,6 +82,30 @@ def test_df_rhf_grad(self):
         e2 = mfs('O  0.  0.  0.001; H  0.  -0.757  0.587; H  0.  0.757   0.587')
         self.assertAlmostEqual(g[0,2], (e2-e1)/0.002*lib.param.BOHR, 5)
 
+    @unittest.skipIf('dftd3' not in sys.modules, "requires the dftd3 library")
+    def test_rhf_d3_grad(self):
+        mf = scf.RHF(mol)
+        mf.disp = 'd3bj'
+        g_scan = mf.nuc_grad_method().as_scanner()
+        g = g_scan(mol)[1]
+
+        mf_scan = mf.as_scanner()
+        e1 = mf_scan('O  0.  0. -0.001; H  0.  -0.757  0.587; H  0.  0.757   0.587')
+        e2 = mf_scan('O  0.  0.  0.001; H  0.  -0.757  0.587; H  0.  0.757   0.587')
+        self.assertAlmostEqual((e2-e1)/0.002*lib.param.BOHR, g[0,2], 5)
+
+    @unittest.skipIf('dftd4' not in sys.modules, "requires the dftd4 library")
+    def test_rhf_d4_grad(self):
+        mf = scf.RHF(mol)
+        mf.disp = 'd4'
+        g_scan = mf.nuc_grad_method().as_scanner()
+        g = g_scan(mol)[1]
+
+        mf_scan = mf.as_scanner()
+        e1 = mf_scan('O  0.  0. -0.001; H  0.  -0.757  0.587; H  0.  0.757   0.587')
+        e2 = mf_scan('O  0.  0.  0.001; H  0.  -0.757  0.587; H  0.  0.757   0.587')
+        self.assertAlmostEqual((e2-e1)/0.002*lib.param.BOHR, g[0,2], 5)
+
     def test_x2c_rhf_grad(self):
         h2o = gto.Mole()
         h2o.verbose = 0
diff --git a/pyscf/hessian/__init__.py b/pyscf/hessian/__init__.py
index 25c10776b9..11f2c7ebe0 100644
--- a/pyscf/hessian/__init__.py
+++ b/pyscf/hessian/__init__.py
@@ -18,6 +18,7 @@
 
 from pyscf.hessian import rhf
 from pyscf.hessian import uhf
+from pyscf.hessian import dispersion
 from pyscf.hessian.rhf import Hessian as RHF
 from pyscf.hessian.uhf import Hessian as UHF
 from pyscf.hessian.rhf import hess_nuc
diff --git a/pyscf/hessian/dispersion.py b/pyscf/hessian/dispersion.py
new file mode 100644
index 0000000000..8751da7e1c
--- /dev/null
+++ b/pyscf/hessian/dispersion.py
@@ -0,0 +1,94 @@
+#!/usr/bin/env python
+# Copyright 2014-2023 The PySCF Developers. All Rights Reserved.
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+#
+# Author: Xiaojie Wu <wxj6000@gmail.com>
+#
+
+'''
+Hessian of dispersion correction for HF and DFT
+'''
+
+
+import numpy
+from pyscf.scf.hf import KohnShamDFT
+
+def get_dispersion(hessobj, disp_version=None):
+    if disp_version is None:
+        disp_version = hessobj.base.disp
+    mol = hessobj.base.mol
+    natm = mol.natm
+    mf = hessobj.base
+    h_disp = numpy.zeros([natm,natm,3,3])
+    if disp_version is None:
+        return h_disp
+    if isinstance(hessobj.base, KohnShamDFT):
+        method = hessobj.base.xc
+    else:
+        method = 'hf'
+
+    if mf.disp[:2].upper() == 'D3':
+        import dftd3.pyscf as disp
+        coords = hessobj.mol.atom_coords()
+        mol = mol.copy()
+        eps = 1e-5
+        for i in range(natm):
+            for j in range(3):
+                coords[i,j] += eps
+                mol.set_geom_(coords, unit='Bohr')
+                d3 = disp.DFTD3Dispersion(mol, xc=method, version=mf.disp)
+                _, g1 = d3.kernel()
+
+                coords[i,j] -= 2.0*eps
+                mol.set_geom_(coords, unit='Bohr')
+                d3 = disp.DFTD3Dispersion(mol, xc=method, version=mf.disp)
+                _, g2 = d3.kernel()
+
+                coords[i,j] += eps
+                h_disp[i,:,j,:] = (g1 - g2)/(2.0*eps)
+            return h_disp
+
+    elif mf.disp[:2].upper() == 'D4':
+        from pyscf.data.elements import charge
+        atoms = numpy.array([ charge(a[0]) for a in mol._atom])
+        coords = mol.atom_coords()
+        natm = mol.natm
+        from dftd4.interface import DampingParam, DispersionModel
+        params = DampingParam(method=method)
+        mol = mol.copy()
+        eps = 1e-5
+        for i in range(natm):
+            for j in range(3):
+                coords[i,j] += eps
+                mol.set_geom_(coords, unit='Bohr')
+                model = DispersionModel(atoms, coords)
+                res = model.get_dispersion(params, grad=True)
+                g1 = res.get("gradient")
+
+                coords[i,j] -= 2.0*eps
+                mol.set_geom_(coords, unit='Bohr')
+                model = DispersionModel(atoms, coords)
+                res = model.get_dispersion(params, grad=True)
+                g2 = res.get("gradient")
+
+                coords[i,j] += eps
+                h_disp[i,:,j,:] = (g1 - g2)/(2.0*eps)
+
+        return h_disp
+    else:
+        raise RuntimeError(f'dispersion correction: {disp_version} is not supported.')
+
+# Inject to SCF class
+from pyscf import hessian
+hessian.rhf.HessianBase.get_dispersion = get_dispersion
diff --git a/pyscf/hessian/rhf.py b/pyscf/hessian/rhf.py
index 9b1d69fc00..90e7db492d 100644
--- a/pyscf/hessian/rhf.py
+++ b/pyscf/hessian/rhf.py
@@ -27,7 +27,7 @@
 from pyscf import lib
 from pyscf import gto
 from pyscf.lib import logger
-from pyscf.scf import _vhf
+from pyscf.scf import _vhf, hf
 from pyscf.scf import cphf
 
 # import _response_functions to load gen_response methods in SCF class
@@ -583,6 +583,8 @@ def kernel(self, mo_energy=None, mo_coeff=None, mo_occ=None, atmlst=None):
 
         de = self.hess_elec(mo_energy, mo_coeff, mo_occ, atmlst=atmlst)
         self.de = de + self.hess_nuc(self.mol, atmlst=atmlst)
+        if self.base.disp is not None:
+            self.de += self.get_dispersion()
         return self.de
     hess = kernel
 
diff --git a/pyscf/hessian/test/test_rks.py b/pyscf/hessian/test/test_rks.py
index 46ceb795d5..ab7389f015 100644
--- a/pyscf/hessian/test/test_rks.py
+++ b/pyscf/hessian/test/test_rks.py
@@ -18,6 +18,17 @@
 from pyscf import gto, dft, lib
 from pyscf import grad, hessian
 
+import sys
+try:
+    import dftd3
+except ImportError:
+    pass
+
+try:
+    import dftd4
+except ImportError:
+    pass
+
 def setUpModule():
     global mol, h4
     mol = gto.Mole()
@@ -116,6 +127,40 @@ def test_finite_diff_b3lyp_hess(self):
         #FIXME: errors seems too big
         self.assertAlmostEqual(abs(hess[0,:,2] - (e1-e2)/2e-4*lib.param.BOHR).max(), 0, 3)
 
+    @unittest.skipIf('dftd3' not in sys.modules, "requires the dftd3 library")
+    def test_finite_diff_b3lyp_d3_hess(self):
+        mf = dft.RKS(mol)
+        mf.conv_tol = 1e-14
+        mf.xc = 'b3lyp'
+        mf.disp = 'd3bj'
+        e0 = mf.kernel()
+        hess = mf.Hessian().kernel()
+        self.assertAlmostEqual(lib.fp(hess), -0.7586078053657133, 6)
+
+        g_scanner = mf.nuc_grad_method().as_scanner()
+        pmol = mol.copy()
+        e1 = g_scanner(pmol.set_geom_('O  0. 0. 0.0001; 1  0. -0.757 0.587; 1  0. 0.757 0.587'))[1]
+        e2 = g_scanner(pmol.set_geom_('O  0. 0. -.0001; 1  0. -0.757 0.587; 1  0. 0.757 0.587'))[1]
+        #FIXME: errors seems too big
+        self.assertAlmostEqual(abs(hess[0,:,2] - (e1-e2)/2e-4*lib.param.BOHR).max(), 0, 3)
+
+    @unittest.skipIf('dftd4' not in sys.modules, "requires the dftd4 library")
+    def test_finite_diff_b3lyp_d4_hess(self):
+        mf = dft.RKS(mol)
+        mf.conv_tol = 1e-14
+        mf.xc = 'b3lyp'
+        mf.disp = 'd4'
+        e0 = mf.kernel()
+        hess = mf.Hessian().kernel()
+        self.assertAlmostEqual(lib.fp(hess), -0.7588415571313422, 6)
+
+        g_scanner = mf.nuc_grad_method().as_scanner()
+        pmol = mol.copy()
+        e1 = g_scanner(pmol.set_geom_('O  0. 0. 0.0001; 1  0. -0.757 0.587; 1  0. 0.757 0.587'))[1]
+        e2 = g_scanner(pmol.set_geom_('O  0. 0. -.0001; 1  0. -0.757 0.587; 1  0. 0.757 0.587'))[1]
+        #FIXME: errors seems too big
+        self.assertAlmostEqual(abs(hess[0,:,2] - (e1-e2)/2e-4*lib.param.BOHR).max(), 0, 3)
+
     def test_finite_diff_wb97x_hess(self):
         mf = dft.RKS(mol)
         mf.conv_tol = 1e-14
diff --git a/pyscf/scf/__init__.py b/pyscf/scf/__init__.py
index 32ffa220f0..77d895d3d9 100644
--- a/pyscf/scf/__init__.py
+++ b/pyscf/scf/__init__.py
@@ -92,7 +92,7 @@
         SCF converged or not
     e_tot : float
         Total HF energy (electronic energy plus nuclear repulsion)
-    mo_energy : 
+    mo_energy :
         Orbital energies
     mo_occ
         Orbital occupancy
@@ -115,6 +115,7 @@
 from pyscf.scf import chkfile
 from pyscf.scf import addons
 from pyscf.scf import diis
+from pyscf.scf import dispersion
 from pyscf.scf.diis import DIIS, CDIIS, EDIIS, ADIIS
 from pyscf.scf.uhf import spin_square
 from pyscf.scf.hf import get_init_guess
diff --git a/pyscf/scf/dispersion.py b/pyscf/scf/dispersion.py
new file mode 100644
index 0000000000..607ab5c9c7
--- /dev/null
+++ b/pyscf/scf/dispersion.py
@@ -0,0 +1,80 @@
+#!/usr/bin/env python
+# Copyright 2014-2023 The PySCF Developers. All Rights Reserved.
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+#
+# Author: Xiaojie Wu <wxj6000@gmail.com>
+#
+
+'''
+dispersion correction for HF and DFT
+'''
+
+
+import numpy
+from pyscf.scf.hf import KohnShamDFT
+
+def get_dispersion(mf, disp_version=None):
+    if disp_version is None:
+        disp_version = mf.disp
+    mol = mf.mol
+    if disp_version is None:
+        return 0.0
+    if isinstance(mf, KohnShamDFT):
+        method = mf.xc
+    else:
+        method = 'hf'
+
+    # for dftd3
+    if disp_version[:2].upper() == 'D3':
+        try:
+            import dftd3.pyscf as disp
+        except ImportError:
+            raise ImportError("\n \
+cannot find dftd3 in the current environment.\n \
+please install dftd3 via \n \
+**************************************\n\
+        pip3 install dftd3 \n \
+**************************************")
+
+        d3 = disp.DFTD3Dispersion(mol, xc=method, version=disp_version)
+        e_d3, _ = d3.kernel()
+        mf.scf_summary['dispersion'] = e_d3
+        return e_d3
+
+    # for dftd4
+    elif disp_version[:2].upper() == 'D4':
+        from pyscf.data.elements import charge
+        atoms = numpy.array([ charge(a[0]) for a in mol._atom])
+        coords = mol.atom_coords()
+        try:
+            from dftd4.interface import DampingParam, DispersionModel
+        except ImportError:
+            raise ImportError("\n \
+cannot find dftd4 in the current environment. \n \
+please install dftd4 via \n \
+***************************************\n \
+        pip3 install dftd4 \n \
+***************************************")
+
+        model = DispersionModel(atoms, coords)
+        res = model.get_dispersion(DampingParam(method=method), grad=False)
+        e_d4 = res.get("energy")
+        mf.scf_summary['dispersion'] = e_d4
+        return e_d4
+    else:
+        raise RuntimeError(f'dipersion correction: {disp_version} is not supported.')
+
+# Inject to SCF class
+from pyscf import scf
+scf.hf.SCF.get_dispersion = get_dispersion
diff --git a/pyscf/scf/hf.py b/pyscf/scf/hf.py
index 434af8d571..8536176c47 100644
--- a/pyscf/scf/hf.py
+++ b/pyscf/scf/hf.py
@@ -231,6 +231,11 @@ def kernel(mf, conv_tol=1e-10, conv_tol_grad=None,
         if dump_chk:
             mf.dump_chk(locals())
 
+    if mf.disp is not None:
+        e_disp = mf.get_dispersion()
+        mf.scf_summary['dispersion'] = e_disp
+        e_tot += e_disp
+
     logger.timer(mf, 'scf_cycle', *cput0)
     # A post-processing hook before return
     mf.post_kernel(locals())
@@ -1478,7 +1483,7 @@ class SCF(lib.StreamObject):
         'diis_file', 'diis_space_rollback', 'damp', 'level_shift',
         'direct_scf', 'direct_scf_tol', 'conv_check', 'callback',
         'mol', 'chkfile', 'mo_energy', 'mo_coeff', 'mo_occ',
-        'e_tot', 'converged', 'scf_summary', 'opt',
+        'e_tot', 'converged', 'scf_summary', 'opt', 'disp',
     }
 
     def __init__(self, mol):
@@ -1490,6 +1495,7 @@ def __init__(self, mol):
         self.verbose = mol.verbose
         self.max_memory = mol.max_memory
         self.stdout = mol.stdout
+        self.disp = None
 
         # If chkfile is muted, SCF intermediates will not be dumped anywhere.
         if MUTE_CHKFILE:
diff --git a/pyscf/scf/test/test_rhf.py b/pyscf/scf/test/test_rhf.py
index 048192187d..200a495e6e 100644
--- a/pyscf/scf/test/test_rhf.py
+++ b/pyscf/scf/test/test_rhf.py
@@ -24,6 +24,17 @@
 from pyscf import scf
 from pyscf.scf import atom_hf
 
+import sys
+try:
+    import dftd3
+except ImportError:
+    pass
+
+try:
+    import dftd4
+except ImportError:
+    pass
+
 def setUpModule():
     global mol, mf, n2sym, n2mf, re_ecp1, re_ecp2
     mol = gto.M(
@@ -344,6 +355,24 @@ def test_analyze(self):
     def test_scf(self):
         self.assertAlmostEqual(mf.e_tot, -76.026765673119627, 9)
 
+    @unittest.skipIf('dftd3' not in sys.modules, "requires the dftd3 library")
+    def test_scf_d3(self):
+        mf = scf.RHF(mol)
+        mf.disp = 'd3bj'
+        mf.conv_tol = 1e-10
+        mf.chkfile = None
+        e_tot = mf.kernel()
+        self.assertAlmostEqual(e_tot, -76.03127458778653, 9)
+
+    @unittest.skipIf('dftd4' not in sys.modules, "requires the dftd4 library")
+    def test_scf_d4(self):
+        mf = scf.RHF(mol)
+        mf.disp = 'd4'
+        mf.conv_tol = 1e-10
+        mf.chkfile = None
+        e_tot = mf.kernel()
+        self.assertAlmostEqual(e_tot, -76.0277467547733, 9)
+
     def test_scf_negative_spin(self):
         mol = gto.M(atom = '''
         O     0    0        0
diff --git a/pyscf/solvent/grad/pcm.py b/pyscf/solvent/grad/pcm.py
index 6c92278de4..2d666a39bf 100644
--- a/pyscf/solvent/grad/pcm.py
+++ b/pyscf/solvent/grad/pcm.py
@@ -174,7 +174,7 @@ def grad_kernel(pcmobj, dm):
     dvj = numpy.zeros([nao,3])
     dq = numpy.zeros([ngrids,3])
     for p0, p1 in lib.prange(0, ngrids, blksize):
-        fakemol = gto.fakemol_for_charges(grid_coords[p0:p1], expnt=exponents**2)
+        fakemol = gto.fakemol_for_charges(grid_coords[p0:p1], expnt=exponents[p0:p1]**2)
         # charge response
         v_nj_ip1 = df.incore.aux_e2(mol, fakemol, intor='int3c2e_ip1', aosym='s1', comp=3)
         vj = numpy.einsum('xijn,n->xij', v_nj_ip1, q_sym)
diff --git a/pyscf/solvent/pcm.py b/pyscf/solvent/pcm.py
index 349b6b8c5c..fc6412d292 100644
--- a/pyscf/solvent/pcm.py
+++ b/pyscf/solvent/pcm.py
@@ -374,7 +374,7 @@ def _get_v(self, surface, dms):
         cintopt = gto.moleintor.make_cintopt(mol._atm, mol._bas, mol._env, int3c2e)
         v_grids_e = numpy.empty(ngrids)
         for p0, p1 in lib.prange(0, ngrids, blksize):
-            fakemol = gto.fakemol_for_charges(grid_coords[p0:p1], expnt=exponents**2)
+            fakemol = gto.fakemol_for_charges(grid_coords[p0:p1], expnt=exponents[p0:p1]**2)
             v_nj = df.incore.aux_e2(mol, fakemol, intor=int3c2e, aosym='s1', cintopt=cintopt)
             v_grids_e[p0:p1] = numpy.einsum('ijL,ij->L',v_nj, dms[0])