Add basic support for GCP in Python API

python/dftd3/interface.py

@@ -459,6 +459,59 @@ def get_pairwise_dispersion(self, param: DampingParam) -> dict:
         }
 
 
+class GeometricCounterpoise(Structure):
+    """
+    .. Counterpoise correction parameters
+
+    Contains the required information to evaluate the counterpoise correction
+    for a given geometry. The counterpoise correction is a method to correct
+    the interaction energy in a supermolecular calculation for the basis set
+    superposition error (BSSE).
+    """
+
+    _gcp = library.ffi.NULL
+
+    def __init__(
+        self,
+        numbers: np.ndarray,
+        positions: np.ndarray,
+        lattice: Optional[np.ndarray] = None,
+        periodic: Optional[np.ndarray] = None,
+        method: Optional[str] = None,
+        basis: Optional[str] = None,
+    ):
+        Structure.__init__(self, numbers, positions, lattice, periodic)
+
+        _method = library.ffi.new("char[]", method.encode()) if method else library.ffi.NULL
+        _basis = library.ffi.new("char[]", basis.encode()) if basis else library.ffi.NULL
+        self._gcp = library.load_gcp_param(self._mol, _method, _basis)
+
+    def set_realspace_cutoff(self, bas: float, srb: float):
+        """Set realspace cutoff for evaluation of interactions"""
+
+        library.set_gcp_realspace_cutoff(self._disp, bas, srb)
+
+    def get_counterpoise(self, grad: bool) -> dict:
+        """Evaluate the counterpoise corrected interaction energy"""
+
+        _energy = np.array(0.0)
+        if grad:
+            _gradient = np.zeros((len(self), 3))
+            _sigma = np.zeros((3, 3))
+        else:
+            _gradient = None
+            _sigma = None
+
+        library.get_counterpoise(self._mol, self._gcp, _cast("double*", _energy), _cast("double*", _gradient), _cast("double*", _sigma))
+
+        results = dict(energy=_energy)
+        if _gradient is not None:
+            results.update(gradient=_gradient)
+        if _sigma is not None:
+            results.update(virial=_sigma)
+        return results
+
+
 def _cast(ctype, array):
     """Cast a numpy array to a FFI pointer"""
     return (

python/dftd3/library.py

@@ -191,6 +191,25 @@ def load_optimizedpower_damping(method, atm):
 get_pairwise_dispersion = error_check(lib.dftd3_get_pairwise_dispersion)
 
 
+def _delete_gcp(gcp):
+    """Delete a geometric counter-poise object"""
+    ptr = ffi.new("dftd3_gcp *")
+    ptr[0] = gcp
+    lib.dftd3_delete_gcp(ptr)
+
+
+def load_gcp_param(mol, method: str, basis: str):
+    """Load GCP parameters from internal storage"""
+    return ffi.gc(
+        error_check(lib.dftd3_load_gcp_param)(mol, method.encode(), basis.encode()),
+        _delete_gcp,
+    )
+
+
+set_gcp_realspace_cutoff = error_check(lib.dftd3_set_gcp_realspace_cutoff)
+get_counterpoise = error_check(lib.dftd3_get_counterpoise)
+
+
 def _ref(ctype, value):
     """Create a reference to a value"""
     if value is None:

python/dftd3/pyscf.py

@@ -267,7 +267,7 @@ class _DFTD3Grad:
     pass
 
 
-def energy(mf: scf.hf.SCF, **kwargs) -> scf.hf.SCF:
+def d3_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.
@@ -301,7 +301,7 @@ def energy(mf: scf.hf.SCF, **kwargs) -> scf.hf.SCF:
     ...          H   2.15862174  -0.13639605   0.80956529
     ...          '''
     ... )
-    >>> mf = disp.energy(scf.RHF(mol)).run()
+    >>> mf = disp.d3_energy(scf.RHF(mol)).run()
     converged SCF energy = -110.932603617026
     >>> mf.kernel()
     -110.93260361702605
@@ -355,7 +355,7 @@ def nuc_grad_method(self):
     return DFTD3(mf, with_dftd3)
 
 
-def grad(scf_grad: GradientsBase, **kwargs):
+def d3_grad(scf_grad: GradientsBase, **kwargs):
     """
     Apply DFT-D3 corrections to SCF or MCSCF nuclear gradients methods
     by returning an instance of a new class built from the original class.
@@ -387,7 +387,7 @@ def grad(scf_grad: GradientsBase, **kwargs):
     ...          H   1.57598558  -0.38252146   0.75856129
     ...          '''
     ... )
-    >>> grad = disp.energy(scf.RHF(mol)).run().nuc_grad_method()
+    >>> grad = disp.d3_energy(scf.RHF(mol)).run().nuc_grad_method()
     converged SCF energy = -149.947191000075
     >>> g = grad.kernel()
     --------------- DFTD3 gradients ---------------
@@ -422,3 +422,7 @@ def grad_nuc(self, mol=None, atmlst=None):
     mfgrad = DFTD3Grad.__new__(DFTD3Grad)
     mfgrad.__dict__.update(scf_grad.__dict__)
     return mfgrad
+
+
+energy = d3_energy
+grad = d3_grad

src/dftd3/gcp.f90

@@ -173,6 +173,12 @@ subroutine gcp_energy(mol, trans, cutoff, iz, emiss, slater, xv, rvdw, escal, al
    real(wp) :: dampval, grd_dmp
    real(wp) :: dE
 
+   !$omp parallel do default(none) &
+   !$omp reduction(+:energies) &
+   !$omp shared(mol, iz, xv, emiss, rvdw, trans, cutoff, alpha, beta, &
+   !$omp&       dmp_scal, dmp_exp, escal, slater, damp) &
+   !$omp private(izp, jat, jzp, xvi, xvj, emi, emj, r0, jtr, vec, r1, sij, &
+   !$omp&        expv, bsse, argv, dE, dampval, grd_dmp, rscal, rscalexp)
    do iat = 1, mol%nat
       izp = mol%id(iat)
       xvi = merge(1.0_wp / sqrt(xv(izp)), 0.0_wp, xv(izp) >= 0.5_wp)
@@ -262,6 +268,13 @@ subroutine gcp_deriv(mol, trans, cutoff, iz, emiss, slater, xv, rvdw, escal, alp
    real(wp) :: dampval, grd_dmp
    real(wp) :: dE, dG(3), dS(3, 3)
 
+   !$omp parallel do default(none) &
+   !$omp reduction(+:energies, gradient, sigma) &
+   !$omp shared(mol, iz, xv, emiss, rvdw, trans, cutoff, alpha, beta, &
+   !$omp&       dmp_scal, dmp_exp, escal, slater, damp) &
+   !$omp private(izp, jat, jzp, xvi, xvj, emi, emj, r0, jtr, vec, r1, sij, gij, emij, &
+   !$omp&        expv, bsse, argv, dE, dampval, grd_dmp, dG, dS, expd, argd, ovlpd, gs, &
+   !$omp&        rscal, rscalexp, rscalexpm1)
    do iat = 1, mol%nat
       izp = mol%id(iat)
       xvi = merge(1.0_wp / sqrt(xv(izp)), 0.0_wp, xv(izp) >= 0.5_wp)
@@ -358,6 +371,10 @@ subroutine srb_energy(mol, trans, cutoff, iz, r0ab, rscal, qscal, rexp, zexp, en
    real(wp) :: r0, vec(3), dE
    integer iat, jat, jtr, izp, jzp
 
+   !$omp parallel do default(none) &
+   !$omp reduction(+:energies) &
+   !$omp shared(mol, trans, cutoff, iz, r0ab, rexp, zexp, rscal, qscal) &
+   !$omp private(izp, jat, jzp, r0, vec, r1, fi, fj, ff, expt, dE)
    do iat = 1, mol%nat
       izp = mol%id(iat)
       fi = real(iz(izp), wp)
@@ -416,6 +433,10 @@ subroutine srb_deriv(mol, trans, cutoff, iz, r0ab, rscal, qscal, rexp, zexp, ene
    real(wp) :: r0, vec(3), dE, dG(3), dS(3, 3)
    integer iat, jat, jtr, izp, jzp
 
+   !$omp parallel do default(none) &
+   !$omp reduction(+:energies, gradient, sigma) &
+   !$omp shared(mol, trans, cutoff, iz, r0ab, rexp, zexp, rscal, qscal) &
+   !$omp private(izp, jat, jzp, r0, vec, r1, fi, fj, ff, expt, rf, dE, dG, dS)
    do iat = 1, mol%nat
       izp = mol%id(iat)
       fi = real(iz(izp), wp)
@@ -468,8 +489,7 @@ subroutine ssovl(r, iat, jat, iz, xza, xzb, ovl)
    logical debug
    real(wp) za, zb, R, ovl, ax, bx, norm, R05
    integer na, nb
-   real(wp) Bxx0, Bxx1, Bxx2, xx, Bxx4, Bxx6
-   real(wp) Bxx3, Bxx5
+   real(wp) xx
    data shell/                 &
    !          h, he
             1, 1               &

src/dftd3/output.f90

@@ -470,6 +470,7 @@ subroutine turbomole_gradlatt(mol, fname, energy, sigma, stat)
    icycle = 1
    i = 0
    escf = 0.0_wp
+   line_number = 0
 
    inquire(file=fname,exist=exist)
    if (exist) then
@@ -549,6 +550,7 @@ subroutine turbomole_gradient(mol, fname, energy, gradient, stat)
    icycle = 1
    i = 0
    escf = 0.0_wp
+   line_number = 0
 
    inquire(file=fname,exist=exist)
    if (exist) then