Skip to content

Commit

Permalink
Merge pull request #372 from tovrstra/convert-module
Browse files Browse the repository at this point in the history 
Move convert functions to new convert module
  • Loading branch information
@tovrstra
tovrstra authored Jul 11, 2024
2 parents d27e68c + 0e0ac03 commit cdb30cc
Show file tree
Hide file tree
Showing 19 changed files with 558 additions and 461 deletions.
195 changes: 1 addition & 194 deletions iodata/basis.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -29,22 +29,11 @@
from typing import Union

import attrs
import numpy as np
from numpy.typing import NDArray

from .attrutils import convert_array_to, validate_shape

__all__ = (
"angmom_sti",
"angmom_its",
"Shell",
"MolecularBasis",
"convert_convention_shell",
"convert_conventions",
"iter_cart_alphabet",
"HORTON2_CONVENTIONS",
"CCA_CONVENTIONS",
)
__all__ = ("angmom_sti", "angmom_its", "Shell", "MolecularBasis")

ANGMOM_CHARS = "spdfghiklmnoqrtuvwxyzabce"

Expand Down Expand Up @@ -260,185 +249,3 @@ def get_segmented(self):
Shell(shell.icenter, [angmom], [kind], shell.exponents, coeffs.reshape(-1, 1))
)
return attrs.evolve(self, shells=shells)


def convert_convention_shell(
conv1: list[str], conv2: list[str], reverse=False
) -> tuple[NDArray[int], NDArray[int]]:
"""Return a permutation vector and sign changes to convert from 1 to 2.
The transformation from convention 1 to convention 2 can be done applying
the results of this function as follows:
.. code-block:: python
vector2 = vector1[permutation]*signs
When using the option ``reverse=True``, one can use the results to convert
in the opposite sense:
.. code-block:: python
vector1 = vector2[permutation]*signs
Parameters
----------
conv1, conv2
Two lists, with the same strings (in different order), where each string
may be prefixed with a '-'.
reverse:
When, true the conversion from 2 to 1 is returned.
Returns
-------
permutation
An integer array that permutes basis function from 1 to 2.
signs
Sign changes when going from 1 to 2, must be applied after permutation
"""
if len(conv1) != len(conv2):
raise TypeError("conv1 and conv2 must contain the same number of elements.")
# Get signs from both
signs1 = [1 - 2 * el1.startswith("-") for el1 in conv1]
signs2 = [1 - 2 * el2.startswith("-") for el2 in conv2]
# Strip signs from both
conv1 = [el1.lstrip("-") for el1 in conv1]
conv2 = [el2.lstrip("-") for el2 in conv2]
if len(conv1) != len(set(conv1)):
raise TypeError("Argument conv1 contains duplicates.")
if len(conv2) != len(set(conv2)):
raise TypeError("Argument conv2 contains duplicates.")
if set(conv1) != set(conv2):
raise TypeError(
"Without the minus signs, conv1 and conv2 must contain "
f"the same elements. Got {conv1} and {conv2}."
)
# Get the permutation
if reverse:
permutation = [conv2.index(el1) for el1 in conv1]
signs = [signs2[i] * sign1 for i, sign1 in zip(permutation, signs1)]
else:
permutation = [conv1.index(el2) for el2 in conv2]
signs = [signs1[i] * sign2 for i, sign2 in zip(permutation, signs2)]
return permutation, signs


def convert_conventions(
molbasis: MolecularBasis, new_conventions: dict[str, list[str]], reverse=False
) -> tuple[NDArray[int], NDArray[int]]:
"""Return a permutation vector and sign changes to convert from 1 to 2.
The transformation from molbasis.convention to the new convention can be done
applying the results of this function as follows:
.. code-block:: python
vector2 = vector1[permutation]*signs
When using the option ``reverse=True``, one can use the results to convert
in the opposite sense:
.. code-block:: python
vector1 = vector2[permutation]*signs
Parameters
----------
molbasis
The description of a molecular basis set.
new_conventions
The new conventions for ordering and signs, to which data for the
orbital basis needs to be converted.
reverse:
When, true the conversion from 2 to 1 is returned.
Returns
-------
permutation
An integer array that permutes basis function from 1 to 2.
signs
Sign changes when going from 1 to 2, must be applied after permutation
"""
permutation = []
signs = []
for shell in molbasis.shells:
for angmom, kind in zip(shell.angmoms, shell.kinds):
key = (angmom, kind)
conv1 = molbasis.conventions[key]
conv2 = new_conventions[key]
shell_permutation, shell_signs = convert_convention_shell(conv1, conv2, reverse)
offset = len(permutation)
permutation.extend(i + offset for i in shell_permutation)
signs.extend(shell_signs)
return np.array(permutation), np.array(signs)


def iter_cart_alphabet(n: int) -> NDArray[int]:
"""Loop over powers of Cartesian basis functions in alphabetical order.
See https://theochem.github.io/horton/2.1.1/tech_ref_gaussian_basis.html
for details.
Parameters
----------
n
The angular momentum, i.e. sum of Cartesian powers in this case.
"""
for nx in range(n, -1, -1):
for ny in range(n - nx, -1, -1):
nz = n - nx - ny
yield np.array((nx, ny, nz), dtype=int)


def get_default_conventions() -> tuple[dict, dict]:
"""Produce conventions dictionaries compatible with HORTON2 and CCA.
Do not change this! Both conventions are also used by several file formats
from other QC codes.
Common Component Architecture (CCA) conventions are defined in appendix B of
the following article:
Kenny, J. P.; Janssen, C. L.; Valeev, E. F.; Windus, T. L. Components for
Integral Evaluation in Quantum Chemistry: Components for Integral Evaluation
in Quantum Chemistry. J. Comput. Chem. 2008, 29 (4), 562-577.
https://doi.org/10.1002/jcc.20815.
The ordering of the spherical harmonics within one shell is rather vague
in appendix B and a more precise description is given on the LibInt Wiki:
https://github.com/evaleev/libint/wiki/using-modern-CPlusPlus-API
Returns
-------
horton2_conventions
A conventions dictionary for HORTON2, of which parts are used by various
file formats.
cca_conventions
A conventions dictionary compatible with the Common Component
Architecture (CCA).
"""
horton2 = {(0, "c"): ["1"]}
cca = horton2.copy()
for angmom in range(1, 25):
conv_cart = ["x" * nx + "y" * ny + "z" * nz for nx, ny, nz in iter_cart_alphabet(angmom)]
key = (angmom, "c")
horton2[key] = conv_cart
cca[key] = conv_cart
if angmom > 1:
conv_pure = ["c0"]
for absm in range(1, angmom + 1):
conv_pure.append(f"c{absm}")
conv_pure.append(f"s{absm}")
key = (angmom, "p")
horton2[key] = conv_pure
cca[key] = conv_pure[:1:-2] + conv_pure[:1] + conv_pure[1::2]
return horton2, cca


HORTON2_CONVENTIONS, CCA_CONVENTIONS = get_default_conventions()
Loading

0 comments on commit cdb30cc

Please sign in to comment.