Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Add transform= keyword to cbasis integrals (#178) #179

Merged
merged 4 commits into from
May 2, 2024
Merged

Add transform= keyword to cbasis integrals (#178) #179

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
51eb3d5
Add transform= keyword to cbasis integrals
msricher Apr 30, 2024
7d93540
Expose transform= keyword
msricher May 1, 2024
d51fec6
[pre-commit.ci] auto fixes from pre-commit.com hooks
pre-commit-ci[bot] May 1, 2024
ecdd01c
Add tests for CBasis transform argument
leila-pujal May 2, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
138 changes: 105 additions & 33 deletions gbasis/integrals/libcint.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -670,7 +670,7 @@ def make_int1e(
)

# Make instance-bound integral method
def int1e(notation="physicist", origin=None, inv_origin=None):
def int1e(notation="physicist", transform=None, origin=None, inv_origin=None):
# Handle ``notation`` argument
if notation not in ("physicist", "chemist"):
raise ValueError("``notation`` must be one of 'physicist' or 'chemist'")
Expand Down Expand Up @@ -755,11 +755,17 @@ def int1e(notation="physicist", origin=None, inv_origin=None):
# Apply permutation
out = out[self._permutations, :][:, self._permutations]

# Return normalized integrals
# Normalize integrals
if self.coord_type == "cartesian":
return np.einsum(norm_einsum, self._ovlp_minhalf, self._ovlp_minhalf, out)
else:
return out
out = np.einsum(norm_einsum, self._ovlp_minhalf, self._ovlp_minhalf, out)

# Apply transformation
if transform is not None:
out = np.tensordot(transform, out, (1, 0))
out = np.tensordot(transform, out, (1, 1))
out = np.swapaxes(out, 0, 1)

return out

# Return instance-bound integral method
return int1e
Expand Down Expand Up @@ -823,7 +829,7 @@ def make_int2e(
)

# Make instance-bound integral method
def int2e(notation="physicist", origin=None, inv_origin=None):
def int2e(notation="physicist", transform=None, origin=None, inv_origin=None):
# Handle ``notation`` argument
if notation == "physicist":
physicist = True
Expand Down Expand Up @@ -952,101 +958,133 @@ def int2e(notation="physicist", origin=None, inv_origin=None):
if constant is not None:
out *= constant

# Transpose integrals in `out` array to proper notation
if physicist:
out = out.transpose(0, 2, 1, 3)

# Apply permutation
out = out[self._permutations]
out = out[:, self._permutations]
out = out[:, :, self._permutations]
out = out[:, :, :, self._permutations]

# Return normalized integrals
# Normalize integrals
if self.coord_type == "cartesian":
return np.einsum(
out = np.einsum(
norm_einsum,
self._ovlp_minhalf,
self._ovlp_minhalf,
self._ovlp_minhalf,
self._ovlp_minhalf,
out,
)
else:
return out

# Transpose integrals in `out` array to proper notation
if physicist:
out = out.transpose(0, 2, 1, 3)

# Apply transformation
if transform is not None:
out = np.tensordot(transform, out, (1, 0))
out = np.tensordot(transform, out, (1, 1))
out = np.tensordot(transform, out, (1, 2))
out = np.tensordot(transform, out, (1, 3))
out = np.swapaxes(np.swapaxes(out, 0, 3), 1, 2)

return out

# Return instance-bound integral method
return int2e

def overlap_integral(self, notation="physicist"):
def overlap_integral(self, notation="physicist", transform=None):
r"""
Compute the overlap integrals.

Parameters
----------
notation : ("physicist" | "chemist"), default="physicist"
Axis order convention.
transform : np.ndarray(K, K_cont)
Transformation matrix from the basis set in the given coordinate system (e.g. AO) to linear
combinations of contractions (e.g. MO).
Transformation is applied to the left, i.e. the sum is over the index 1 of `transform`
and index 0 of the array for contractions.
Default is no transformation.

Returns
-------
out : np.ndarray(Nbasis, Nbasis, dtype=float)
Integral array.

"""
return self._ovlp(notation=notation)
return self._ovlp(notation=notation, transform=transform)

def kinetic_energy_integral(self, notation="physicist"):
def kinetic_energy_integral(self, notation="physicist", transform=None):
r"""
Compute the kinetic energy integrals.

Parameters
----------
notation : ("physicist" | "chemist"), default="physicist"
Axis order convention.
transform : np.ndarray(K, K_cont)
Transformation matrix from the basis set in the given coordinate system (e.g. AO) to linear
combinations of contractions (e.g. MO).
Transformation is applied to the left, i.e. the sum is over the index 1 of `transform`
and index 0 of the array for contractions.
Default is no transformation.

Returns
-------
out : np.ndarray(Nbasis, Nbasis, dtype=float)
Integral array.

"""
return self._kin(notation=notation)
return self._kin(notation=notation, transform=transform)

def nuclear_attraction_integral(self, notation="physicist"):
def nuclear_attraction_integral(self, notation="physicist", transform=None):
r"""
Compute the nuclear attraction integrals.

Parameters
----------
notation : ("physicist" | "chemist"), default="physicist"
Axis order convention.
transform : np.ndarray(K, K_cont)
Transformation matrix from the basis set in the given coordinate system (e.g. AO) to linear
combinations of contractions (e.g. MO).
Transformation is applied to the left, i.e. the sum is over the index 1 of `transform`
and index 0 of the array for contractions.
Default is no transformation.

Returns
-------
out : np.ndarray(Nbasis, Nbasis, dtype=float)
Integral array.

"""
return self._nuc(notation=notation)
return self._nuc(notation=notation, transform=transform)

def electron_repulsion_integral(self, notation="physicist"):
def electron_repulsion_integral(self, notation="physicist", transform=None):
r"""
Compute the electron repulsion integrals.

Parameters
----------
notation : ("physicist" | "chemist"), default="physicist"
Axis order convention.
transform : np.ndarray(K, K_cont)
Transformation matrix from the basis set in the given coordinate system (e.g. AO) to linear
combinations of contractions (e.g. MO).
Transformation is applied to the left, i.e. the sum is over the index 1 of `transform`
and index 0 of the array for contractions.
Default is no transformation.

Returns
-------
out : np.ndarray(Nbasis, Nbasis, Nbasis, Nbasis, dtype=float)
Integral array.

"""
return self._eri(notation=notation)
return self._eri(notation=notation, transform=transform)

def r_inv_integral(self, origin=None, notation="physicist"):
def r_inv_integral(self, origin=None, notation="physicist", transform=None):
r"""
Compute the :math:`1/\left|\mathbf{r} - \mathbf{R}_\text{inv}\right|` integrals.

Expand All @@ -1056,16 +1094,22 @@ def r_inv_integral(self, origin=None, notation="physicist"):
Origin about which to evaluate integrals.
notation : ("physicist" | "chemist"), default="physicist"
Axis order convention.
transform : np.ndarray(K, K_cont)
Transformation matrix from the basis set in the given coordinate system (e.g. AO) to linear
combinations of contractions (e.g. MO).
Transformation is applied to the left, i.e. the sum is over the index 1 of `transform`
and index 0 of the array for contractions.
Default is no transformation.

Returns
-------
out : np.ndarray(Nbasis, Nbasis, dtype=float)
Integral array.

"""
return self._rinv(inv_origin=origin, notation=notation)
return self._rinv(inv_origin=origin, notation=notation, transform=transform)

def momentum_integral(self, origin=None, notation="physicist"):
def momentum_integral(self, origin=None, notation="physicist", transform=None):
r"""
Compute the momentum integrals.

Expand All @@ -1075,16 +1119,22 @@ def momentum_integral(self, origin=None, notation="physicist"):
Origin about which to evaluate integrals.
notation : ("physicist" | "chemist"), default="physicist"
Axis order convention.
transform : np.ndarray(K, K_cont)
Transformation matrix from the basis set in the given coordinate system (e.g. AO) to linear
combinations of contractions (e.g. MO).
Transformation is applied to the left, i.e. the sum is over the index 1 of `transform`
and index 0 of the array for contractions.
Default is no transformation.

Returns
-------
out : np.ndarray(Nbasis, Nbasis, 3, dtype=complex)
Integral array.

"""
return self._mom(origin=origin, notation=notation)
return self._mom(origin=origin, notation=notation, transform=transform)

def angular_momentum_integral(self, origin=None, notation="physicist"):
def angular_momentum_integral(self, origin=None, notation="physicist", transform=None):
r"""
Compute the angular momentum integrals.

Expand All @@ -1094,6 +1144,12 @@ def angular_momentum_integral(self, origin=None, notation="physicist"):
Origin about which to evaluate integrals.
notation : ("physicist" | "chemist"), default="physicist"
Axis order convention.
transform : np.ndarray(K, K_cont)
Transformation matrix from the basis set in the given coordinate system (e.g. AO) to linear
combinations of contractions (e.g. MO).
Transformation is applied to the left, i.e. the sum is over the index 1 of `transform`
and index 0 of the array for contractions.
Default is no transformation.

Returns
-------
Expand All @@ -1102,9 +1158,11 @@ def angular_momentum_integral(self, origin=None, notation="physicist"):

"""
raise NotImplementedError("Angular momentum integral doesn't work; see Issue #149")
# return self._amom(origin=origin, notation=notation)
# return self._amom(origin=origin, notation=notation, transform=transform)

def point_charge_integral(self, point_coords, point_charges, notation="physicist"):
def point_charge_integral(
self, point_coords, point_charges, notation="physicist", transform=None
):
r"""
Compute the point charge integrals.

Expand All @@ -1116,6 +1174,12 @@ def point_charge_integral(self, point_coords, point_charges, notation="physicist
Charges of point charges.
notation : ("physicist" | "chemist"), default="physicist"
Axis order convention.
transform : np.ndarray(K, K_cont)
Transformation matrix from the basis set in the given coordinate system (e.g. AO) to linear
combinations of contractions (e.g. MO).
Transformation is applied to the left, i.e. the sum is over the index 1 of `transform`
and index 0 of the array for contractions.
Default is no transformation.

Returns
-------
Expand All @@ -1127,13 +1191,13 @@ def point_charge_integral(self, point_coords, point_charges, notation="physicist
out = np.zeros((self.nbfn, self.nbfn, len(point_charges)), dtype=c_double, order="F")
# Compute 1/|r - r_{inv}| for each charge
for icharge, (coord, charge) in enumerate(zip(point_coords, point_charges)):
val = self._rinv(inv_origin=coord, notation=notation)
val = self._rinv(inv_origin=coord, notation=notation, transform=transform)
val *= -charge
out[:, :, icharge] = val
# Return integrals in `out` array
return out

def moment_integral(self, orders, origin=None, notation="physicist"):
def moment_integral(self, orders, origin=None, notation="physicist", transform=None):
r"""
Compute the moment integrals.

Expand All @@ -1145,6 +1209,12 @@ def moment_integral(self, orders, origin=None, notation="physicist"):
Origin about which to evaluate integrals.
notation : ("physicist" | "chemist"), default="physicist"
Axis order convention.
transform : np.ndarray(K, K_cont)
Transformation matrix from the basis set in the given coordinate system (e.g. AO) to linear
combinations of contractions (e.g. MO).
Transformation is applied to the left, i.e. the sum is over the index 1 of `transform`
and index 0 of the array for contractions.
Default is no transformation.

Returns
-------
Expand All @@ -1164,9 +1234,11 @@ def moment_integral(self, orders, origin=None, notation="physicist"):
try:
for i, order in enumerate(orders):
if sum(order) == 0:
out[:, :, i] = self._ovlp(notation=notation)
out[:, :, i] = self._ovlp(notation=notation, transform=transform)
else:
out[:, :, i] = self._moments[tuple(order)](origin=origin, notation=notation)
out[:, :, i] = self._moments[tuple(order)](
origin=origin, notation=notation, transform=transform
)
except KeyError:
raise ValueError(
"Invalid order; can use up to order 4 for any XYZ component,"
Expand Down
Loading
Loading