diff --git a/python/dftd3/pyscf.py b/python/dftd3/pyscf.py
index d0b6be89..22d0cc4c 100644
--- a/python/dftd3/pyscf.py
+++ b/python/dftd3/pyscf.py
@@ -21,12 +21,13 @@
 """
 
 try:
-    from pyscf import lib, gto
+    from pyscf import gto, lib, mcscf, scf
+    from pyscf.grad import rhf as rhf_grad
 except ModuleNotFoundError:
     raise ModuleNotFoundError("This submodule requires pyscf installed")
 
 import numpy as np
-from typing import Optional, Tuple
+from typing import Dict, Optional, Tuple
 
 from .interface import (
     DispersionModel,
@@ -210,7 +211,7 @@ class _DFTD3Grad:
     pass
 
 
-def energy(mf: scf.hf.Scf, **kwargs) -> scf.hf.Scf:
+def energy(mf: scf.hf.SCF, **kwargs) -> scf.hf.SCF:
     """
     Apply DFT-D3 corrections to SCF or MCSCF methods by returning an
     instance of a new class built from the original instances class.
@@ -248,10 +249,7 @@ def energy(mf: scf.hf.Scf, **kwargs) -> scf.hf.Scf:
     -110.93260361702605
     """
 
-    from pyscf.scf import hf
-    from pyscf.mcscf import casci
-
-    if not isinstance(mf, (hf.SCF, casci.CASCI)):
+    if not isinstance(mf, (scf.hf.SCF, mcscf.casci.CASCI)):
         raise TypeError("mf must be an instance of SCF or CASCI")
 
     with_dftd3 = DFTD3Dispersion(
@@ -342,7 +340,6 @@ def grad(scf_grad, **kwargs):
     5 H    -0.0154527822     0.0229409425    -0.0215141991
     ----------------------------------------------
     """
-    from pyscf.grad import rhf as rhf_grad
 
     if not isinstance(scf_grad, rhf_grad.Gradients):
         raise TypeError("scf_grad must be an instance of Gradients")
diff --git a/python/dftd3/test_pyscf.py b/python/dftd3/test_pyscf.py
index 91a6fc77..973650ac 100644
--- a/python/dftd3/test_pyscf.py
+++ b/python/dftd3/test_pyscf.py
@@ -56,6 +56,49 @@ def test_energy_r2scan_d3():
     assert d3.kernel()[0] == approx(-0.00578401192369041, abs=1.0e-7)
 
 
+@pytest.mark.skipif(pyscf is None, reason="requires pyscf")
+@pytest.mark.parametrize("atm", [True, False])
+def test_energy_bp_d3zero(atm):
+    thr = 1e-9
+
+    mol = gto.M(
+        atom=[
+            ("C", [-1.42754169820131, -1.50508961850828, -1.93430551124333]),
+            ("C", [+1.19860572924150, -1.66299114873979, -2.03189643761298]),
+            ("C", [+2.65876001301880, +0.37736955363609, -1.23426391650599]),
+            ("C", [+1.50963368042358, +2.57230374419743, -0.34128058818180]),
+            ("C", [-1.12092277855371, +2.71045691257517, -0.25246348639234]),
+            ("I", [-2.60071517756218, +0.67879949508239, -1.04550707592673]),
+            ("H", [-2.86169588073340, +5.99660765711210, +1.08394899986031]),
+            ("H", [+2.09930989272956, -3.36144811062374, -2.72237695164263]),
+            ("H", [+2.64405246349916, +4.15317840474646, +0.27856972788526]),
+            ("H", [+4.69864865613751, +0.26922271535391, -1.30274048619151]),
+            ("H", [-4.63786461351839, +0.79856258572808, -0.96906659938432]),
+            ("H", [-2.57447518692275, -3.08132039046931, -2.54875517521577]),
+            ("S", [-5.88211879210329, 11.88491819358157, +2.31866455902233]),
+            ("H", [-8.18022701418703, 10.95619984550779, +1.83940856333092]),
+            ("C", [-5.08172874482867, 12.66714386256482, -0.92419491629867]),
+            ("H", [-3.18311711399702, 13.44626574330220, -0.86977613647871]),
+            ("H", [-5.07177399637298, 10.99164969235585, -2.10739192258756]),
+            ("H", [-6.35955320518616, 14.08073002965080, -1.68204314084441]),
+        ],
+        unit="bohr",
+    )
+
+    d3 = disp.DFTD3Dispersion(
+        mol,
+        param={
+            "s8": 1.683,
+            "rs6": 1.139,
+            "s9": 1.0 if atm else 0.0,
+        },
+        version="d3zero",
+    )
+    ref = -0.01410721853585842 if atm else -0.014100267345314462
+
+    assert approx(d3.kernel()[0], abs=thr) == ref
+
+
 @pytest.mark.skipif(pyscf is None, reason="requires pyscf")
 @pytest.mark.parametrize("xc", ["b3lyp", "b3lypg", "b3lyp5", "b3lyp3"])
 def test_energy_b3lyp_d3(xc: str):