From 8018f9406b88f69787aaec6fe4e8e680dac06175 Mon Sep 17 00:00:00 2001 From: Sebastian Ehlert Date: Sun, 1 Dec 2024 14:38:03 +0100 Subject: [PATCH] Add basic support for GCP in Python API --- python/dftd3/interface.py | 53 +++++++++++++++++++++++++++++++++++++++ python/dftd3/library.py | 19 ++++++++++++++ python/dftd3/pyscf.py | 12 ++++++--- src/dftd3/gcp.f90 | 24 ++++++++++++++++-- src/dftd3/output.f90 | 2 ++ 5 files changed, 104 insertions(+), 6 deletions(-) diff --git a/python/dftd3/interface.py b/python/dftd3/interface.py index 844b22a4..a6e8e250 100644 --- a/python/dftd3/interface.py +++ b/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 ( diff --git a/python/dftd3/library.py b/python/dftd3/library.py index 52cb8a6b..8f780726 100644 --- a/python/dftd3/library.py +++ b/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: diff --git a/python/dftd3/pyscf.py b/python/dftd3/pyscf.py index c8de89ad..505eae04 100644 --- a/python/dftd3/pyscf.py +++ b/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 \ No newline at end of file diff --git a/src/dftd3/gcp.f90 b/src/dftd3/gcp.f90 index 293ba66a..98820c3c 100644 --- a/src/dftd3/gcp.f90 +++ b/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 & diff --git a/src/dftd3/output.f90 b/src/dftd3/output.f90 index 0dc3104c..c36065f4 100644 --- a/src/dftd3/output.f90 +++ b/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