Integral driver

src/dxtb/basis/indexhelper.py

@@ -129,6 +129,8 @@ def __init__(
         **_,
     ):
         super().__init__(device, dtype)
+        self.clear_cache()
+
         self.unique_angular = unique_angular
         self.angular = angular
         self.atom_to_unique = atom_to_unique

src/dxtb/basis/type.py

@@ -72,14 +72,14 @@ def __init__(
         self,
         numbers: Tensor,
         par: Param,
-        angular: Tensor,
+        ihelp: IndexHelper,
         device: torch.device | None = None,
         dtype: torch.dtype | None = None,
     ) -> None:
         super().__init__(device, dtype)
         self.numbers = numbers
         self.meta = par.meta
-        self.angular = angular
+        self.ihelp = ihelp
 
         self.ngauss = get_elem_param(
             numbers,
@@ -112,9 +112,12 @@ def create_cgtos(self) -> tuple[list[Tensor], list[Tensor]]:
 
         # maybe this can be batched too, but this loop is rather small
         # so it is probably not worth it
-        for i in range(self.angular.size(0)):
+        for i in range(self.ihelp.unique_angular.size(0)):
             alpha, coeff = slater.to_gauss(
-                self.ngauss[i], self.pqn[i], self.angular[i], self.slater[i]
+                self.ngauss[i],
+                self.pqn[i],
+                self.ihelp.unique_angular[i],
+                self.slater[i],
             )
 
             # FIXME:
@@ -138,7 +141,6 @@ def create_cgtos(self) -> tuple[list[Tensor], list[Tensor]]:
 
     def unique_shell_pairs(
         self,
-        ihelp: IndexHelper,
         mask: Tensor | None = None,
         uplo: Literal["n", "N", "u", "U", "l", "L"] = "l",
     ) -> tuple[Tensor, Tensor]:
@@ -147,8 +149,6 @@ def unique_shell_pairs(
 
         Parameters
         ----------
-        ihelp : IndexHelper
-            Helper class for indexing.
         mask : Tensor
             Mask for
         uplo : Literal["n", "N", "u", "U", "l", "L"]
@@ -168,20 +168,20 @@ def unique_shell_pairs(
         """
         # NOTE: Zeros correspond to padding only if batched. They have
         # meaning for single runs and cannot be ignored in this case.
-        sh2ush = ihelp.spread_shell_to_orbital(ihelp.shells_to_ushell)
+        sh2ush = self.ihelp.spread_shell_to_orbital(self.ihelp.shells_to_ushell)
 
         # FIXME: Maybe a bitwise operation is easier to understand? For now,
         # we convert unique shell indices to prime numbers to obtain unique
         # products upon multiplication (fundamental theorem of arithmetic).
         orbs = primes.to(self.device)[sh2ush]
         orbs = orbs.unsqueeze(-2) * orbs.unsqueeze(-1)
-        sh2orb = ihelp.spread_shell_to_orbital(ihelp.orbitals_per_shell)
+        sh2orb = self.ihelp.spread_shell_to_orbital(self.ihelp.orbitals_per_shell)
 
         # extra offset along only one dimension to distinguish (n, m) and
         # (m, n) of the same orbital block (e.g. 1x3 sp and 3x1 ps block)
         offset = 10000
 
-        if ihelp.batched:
+        if self.ihelp.batched:
             orbs = torch.where(
                 real_pairs(sh2ush),
                 orbs + sh2orb.unsqueeze(-1) * offset,
@@ -192,7 +192,7 @@ def unique_shell_pairs(
 
         # catch systems with one single orbital (e.g. Helium)
         if orbs.size(-1) == 1:
-            if ihelp.batched:
+            if self.ihelp.batched:
                 raise NotImplementedError()
             return torch.tensor([[0]]), torch.tensor(1)
 
@@ -223,7 +223,6 @@ def unique_shell_pairs(
 
     def to_bse(
         self,
-        ihelp: IndexHelper,
         qcformat: Literal["gaussian94", "nwchem"] = "nwchem",
         save: bool = False,
         overwrite: bool = False,
@@ -273,10 +272,10 @@ def to_bse(
             if qcformat == "gaussian94":
                 txt += f"{PSE[number]}\n"
 
-            shells = ihelp.shells_per_atom[i]
+            shells = self.ihelp.shells_per_atom[i]
             for _ in range(shells):
                 alpha, coeff = slater.to_gauss(
-                    self.ngauss[s], self.pqn[s], self.angular[s], self.slater[s]
+                    self.ngauss[s], self.pqn[s], self.ihelp.angular[s], self.slater[s]
                 )
                 if self.valence[s].item() is False:
                     alpha, coeff = orthogonalize(
@@ -286,7 +285,7 @@ def to_bse(
                 alphas.append(alpha)
                 coeffs.append(coeff)
 
-                l = angular2label[self.angular.tolist()[s]]
+                l = angular2label[self.ihelp.angular.tolist()[s]]
                 if qcformat == "gaussian94":
                     txt += f"{l}    {len(alpha)}    1.00\n"
                 elif qcformat == "nwchem":
@@ -336,7 +335,7 @@ def to_bse(
 
         return fulltxt
 
-    def create_dqc(self, positions: Tensor, ihelp: IndexHelper) -> list[AtomCGTOBasis]:
+    def create_dqc(self, positions: Tensor) -> list[AtomCGTOBasis]:
         if self.numbers.ndim > 1:
             raise NotImplementedError("Batch mode not implemented.")
 
@@ -350,11 +349,11 @@ def create_dqc(self, positions: Tensor, ihelp: IndexHelper) -> list[AtomCGTOBasi
         s = 0
         for i, number in enumerate(self.numbers):
             bases: list[CGTOBasis] = []
-            shells = ihelp.shells_per_atom[i]
+            shells = self.ihelp.shells_per_atom[i]
 
             for _ in range(shells):
                 alpha, coeff = slater.to_gauss(
-                    self.ngauss[s], self.pqn[s], ihelp.angular[s], self.slater[s]
+                    self.ngauss[s], self.pqn[s], self.ihelp.angular[s], self.slater[s]
                 )
 
                 # orthogonalize the H2s against the H1s
@@ -367,7 +366,7 @@ def create_dqc(self, positions: Tensor, ihelp: IndexHelper) -> list[AtomCGTOBasi
                 coeffs.append(coeff)
 
                 cgto = CGTOBasis(
-                    angmom=ihelp.angular.tolist()[s],  # int!
+                    angmom=self.ihelp.angular.tolist()[s],  # int!
                     alphas=alpha,
                     coeffs=coeff,
                     normalized=True,

src/dxtb/constants/defaults.py

@@ -40,6 +40,9 @@
 INTCUTOFF = 50.0
 """Real-space cutoff (in Angstrom) for integral evaluation. (50.0)"""
 
+INTDRIVER = "libcint"
+"""Integral driver."""
+
 # SCF settings
 
 GUESS = "eeq"

src/dxtb/integral/__init__.py

@@ -5,4 +5,4 @@
 from .mmd import overlap_gto
 from .mmd.explicit import mmd_explicit
 from .mmd.recursion import mmd_recursion
-from .overlap import Overlap
+from .overlap import *

src/dxtb/integral/libcint/intor.py

@@ -13,7 +13,11 @@
 
 import numpy as np
 import torch
-from dxtblibs import CGTO, CINT
+
+try:
+    from dxtblibs import CGTO, CINT
+except ImportError:
+    pass
 
 from ..._types import Callable, Tensor
 from .namemanager import IntorNameManager

src/dxtb/integral/libcint/namemanager.py

@@ -37,6 +37,9 @@ class IntorNameManager:
             "r_origj": (3,),
             "rr_origj": (9,),
             "rrr_origj": (27,),
+            "j": (3,),
+            "jj": (9,),
+            "jjj": (27,),
         },
     )
     op_comp = defaultdict(

src/dxtb/integral/libcint/symmetry.py

@@ -5,6 +5,7 @@
 In tight-binding, we do not require anything special here. Only S1 symmetry is
 of intereset.
 """
+from __future__ import annotations
 
 from abc import abstractmethod
 

src/dxtb/integral/libcint/utils.py

@@ -4,7 +4,6 @@
 
 This module contains helpers required for calculating integrals with libcint.
 """
-
 from __future__ import annotations
 
 import ctypes

src/dxtb/integral/libcint/wrapper.py

@@ -12,7 +12,11 @@
 
 import numpy as np
 import torch
-from dxtblibs import CINT
+
+try:
+    from dxtblibs import CINT  # type: ignore
+except ImportError:
+    pass
 
 from ..._types import Iterator, Tensor
 from ...basis import AtomCGTOBasis, CGTOBasis, IndexHelper