Dispersion¶

DFT-D4 dispersion energies used by place_guest_in_host(method="energy"): cage_isomer_builder.utils.dispersion.

D4PairModel ¶

D4PairModel(host, guests, method='pbe')

Fixed C6 coefficients and r4r2 values for a host and its guest types.

Atoms are indexed in one list: the host first, then each guest type in order. offsets[t] is where guest type t starts.

Parameters:

Name	Type	Description	Default
`host`	`Atoms`		required
`guests`	`sequence of ase.Atoms`		required
`method`	`str or dict`	D4 damping parameters, see :data:`D4_PARAMETERS`.	`"pbe"`

Examples:

>>> import numpy as np
>>> from ase.build import molecule
>>> a, b = molecule("CH4"), molecule("CH4")
>>> b.translate([4.0, 0.0, 0.0])                    # methane dimer, C...C = 4 A
>>> from cage_isomer_builder.utils.dispersion import D4PairModel   # needs dftd4
>>> model = D4PairModel(a, [b])
>>> model.name, model.offsets
('D4(pbe)+UFF repulsion', [0, 5])

Source code in cage_isomer_builder/utils/dispersion.py

def __init__(self, host, guests, method="pbe"):
    _require_dftd4()
    from dftd4.interface import DispersionModel

    self.name = f"D4({method})+UFF repulsion" if isinstance(method, str) else "D4+UFF repulsion"
    self.method = method
    self.params = _damping(method)
    # Place each fragment far from the others so their coordination
    # numbers are those of the isolated molecules.
    fragments = [host] + list(guests)
    combined = Atoms()
    self.offsets = []
    for n, fragment in enumerate(fragments):
        self.offsets.append(len(combined))
        moved = fragment.copy()
        moved.translate(-moved.get_center_of_mass() + [1000.0 * n, 0.0, 0.0])
        combined += moved
    model = DispersionModel(combined.numbers, combined.positions / units.Bohr)
    self.c6 = np.ascontiguousarray(model.get_properties()["c6 coefficients"])
    self.r4r2 = np.array([_r4r2(int(z)) for z in combined.numbers])
    self.uff_x, self.uff_d = _uff_element_params(combined.numbers)

score ¶

score(source_existing, existing_positions, source_guest, guest_positions, min_distance)

Two-body D4 energy plus UFF repulsion of each guest pose with the existing atoms, and the clash penalty.

Parameters:

Name	Type	Description	Default
`source_existing`	`np.ndarray of int, shape (n_e,)`	Index of each existing atom in this model's atom list.	required
`existing_positions`	`(ndarray, shape(n_e, 3))`		required
`source_guest`	`np.ndarray of int, shape (n_g,)`	Index of each guest atom in this model's atom list.	required
`guest_positions`	`(ndarray, shape(n_poses, n_g, 3))`		required
`min_distance`	`(ndarray, shape(n_e, n_g))`	Smallest allowed distance (Angstrom) for each pair.	required

Returns:

Name	Type	Description
`energy`	`(ndarray, shape(n_poses))`	D4 plus repulsion energy (eV).
`penalty`	`(ndarray, shape(n_poses))`	Clash penalty (eV), zero when no pair is too close.
`overlap`	`(ndarray, shape(n_poses))`	How far (Angstrom) the worst pair is inside its minimum distance.

Examples:

>>> import numpy as np
>>> from ase.build import molecule
>>> a, b = molecule("CH4"), molecule("CH4")
>>> b.translate([4.0, 0.0, 0.0])                    # methane dimer, C...C = 4 A
>>> from cage_isomer_builder.utils.dispersion import D4PairModel
>>> model = D4PairModel(a, [b])
>>> host, guest = np.arange(5), model.offsets[1] + np.arange(5)
>>> energy, penalty, overlap = model.score(host, a.positions, guest, b.positions[None],
...                                        min_distance=np.zeros((5, 5)))
>>> bool(energy[0] < 0), float(penalty[0])          # attractive, no clash
(True, 0.0)

Source code in cage_isomer_builder/utils/dispersion.py

def score(self, source_existing, existing_positions, source_guest,
          guest_positions, min_distance):
    """
    Two-body D4 energy plus UFF repulsion of each guest pose with the
    existing atoms, and the clash penalty.

    Parameters
    ----------
    source_existing : np.ndarray of int, shape (n_e,)
        Index of each existing atom in this model's atom list.
    existing_positions : np.ndarray, shape (n_e, 3)
    source_guest : np.ndarray of int, shape (n_g,)
        Index of each guest atom in this model's atom list.
    guest_positions : np.ndarray, shape (n_poses, n_g, 3)
    min_distance : np.ndarray, shape (n_e, n_g)
        Smallest allowed distance (Angstrom) for each pair.

    Returns
    -------
    energy : np.ndarray, shape (n_poses,)
        D4 plus repulsion energy (eV).
    penalty : np.ndarray, shape (n_poses,)
        Clash penalty (eV), zero when no pair is too close.
    overlap : np.ndarray, shape (n_poses,)
        How far (Angstrom) the worst pair is inside its minimum distance.

    Examples
    --------
    >>> import numpy as np
    >>> from ase.build import molecule
    >>> a, b = molecule("CH4"), molecule("CH4")
    >>> b.translate([4.0, 0.0, 0.0])                    # methane dimer, C...C = 4 A
    >>> from cage_isomer_builder.utils.dispersion import D4PairModel
    >>> model = D4PairModel(a, [b])
    >>> host, guest = np.arange(5), model.offsets[1] + np.arange(5)
    >>> energy, penalty, overlap = model.score(host, a.positions, guest, b.positions[None],
    ...                                        min_distance=np.zeros((5, 5)))
    >>> bool(energy[0] < 0), float(penalty[0])          # attractive, no clash
    (True, 0.0)
    """
    p = self.params
    x = np.sqrt(self.uff_x[source_existing][:, None] * self.uff_x[source_guest][None, :])
    d = np.sqrt(self.uff_d[source_existing][:, None] * self.uff_d[source_guest][None, :])
    energy, penalty, overlap = _score_poses(
        np.ascontiguousarray(existing_positions / units.Bohr),
        np.ascontiguousarray(guest_positions / units.Bohr),
        np.ascontiguousarray(self.c6[np.ix_(source_existing, source_guest)]),
        self.r4r2[source_existing],
        self.r4r2[source_guest],
        np.ascontiguousarray(min_distance / units.Bohr),
        np.ascontiguousarray(x / units.Bohr),
        np.ascontiguousarray(d / units.Hartree),
        p["s6"], p["s8"], p["a1"], p["a2"], _CLASH_PENALTY,
    )
    return energy * units.Hartree, penalty * units.Hartree, overlap * units.Bohr

interaction_energy ¶

interaction_energy(fragments)

Full D4 interaction energy (see d4_interaction_energy) plus the UFF repulsion between the fragments, in eV.

Examples:

>>> import numpy as np
>>> from ase.build import molecule
>>> a, b = molecule("CH4"), molecule("CH4")
>>> b.translate([4.0, 0.0, 0.0])                    # methane dimer, C...C = 4 A
>>> from cage_isomer_builder.utils.dispersion import D4PairModel
>>> energy = D4PairModel(a, [b]).interaction_energy([a, b])
>>> bool(-0.05 < energy < 0)                         # eV: a weak dispersion contact
True

Source code in cage_isomer_builder/utils/dispersion.py

def interaction_energy(self, fragments):
    """
    Full D4 interaction energy (see d4_interaction_energy) plus the UFF
    repulsion between the fragments, in eV.

    Examples
    --------
    >>> import numpy as np
    >>> from ase.build import molecule
    >>> a, b = molecule("CH4"), molecule("CH4")
    >>> b.translate([4.0, 0.0, 0.0])                    # methane dimer, C...C = 4 A
    >>> from cage_isomer_builder.utils.dispersion import D4PairModel
    >>> energy = D4PairModel(a, [b]).interaction_energy([a, b])
    >>> bool(-0.05 < energy < 0)                         # eV: a weak dispersion contact
    True
    """
    return (d4_interaction_energy(fragments, self.method)
            + uff_interaction_energy(fragments, repulsion_only=True))

UFFPairModel ¶

UFFPairModel(host, guests)

UFF Lennard-Jones pair energies for a host and its guest types, with the same interface as :class:D4PairModel. Needs no optional packages.

Atoms are indexed in one list: the host first, then each guest type in order. offsets[t] is where guest type t starts.

Examples:

>>> import numpy as np
>>> from ase.build import molecule
>>> a, b = molecule("CH4"), molecule("CH4")
>>> b.translate([4.0, 0.0, 0.0])                    # methane dimer, C...C = 4 A
>>> from cage_isomer_builder.utils.dispersion import UFFPairModel
>>> model = UFFPairModel(a, [b])
>>> model.name, [int(o) for o in model.offsets]
('UFF', [0, 5])

Source code in cage_isomer_builder/utils/dispersion.py

def __init__(self, host, guests):
    numbers = np.concatenate([host.numbers] + [g.numbers for g in guests])
    self.offsets = list(np.cumsum([0, len(host)] + [len(g) for g in guests])[:-1])
    self.x, self.d = _uff_element_params(numbers)

score ¶

score(source_existing, existing_positions, source_guest, guest_positions, min_distance)

Lennard-Jones energy of each guest pose with the existing atoms, plus the clash penalty. Same arguments and returns as :meth:D4PairModel.score, in eV and Angstrom.

Examples:

>>> import numpy as np
>>> from ase.build import molecule
>>> a, b = molecule("CH4"), molecule("CH4")
>>> b.translate([4.0, 0.0, 0.0])                    # methane dimer, C...C = 4 A
>>> from cage_isomer_builder.utils.dispersion import UFFPairModel, uff_interaction_energy
>>> model = UFFPairModel(a, [b])
>>> host, guest = np.arange(5), model.offsets[1] + np.arange(5)
>>> energy, penalty, overlap = model.score(host, a.positions, guest, b.positions[None],
...                                        min_distance=np.zeros((5, 5)))
>>> bool(np.isclose(energy[0], uff_interaction_energy([a, b]))), float(penalty[0])
(True, 0.0)

Source code in cage_isomer_builder/utils/dispersion.py

def score(self, source_existing, existing_positions, source_guest,
          guest_positions, min_distance):
    """
    Lennard-Jones energy of each guest pose with the existing atoms, plus
    the clash penalty. Same arguments and returns as
    :meth:`D4PairModel.score`, in eV and Angstrom.

    Examples
    --------
    >>> import numpy as np
    >>> from ase.build import molecule
    >>> a, b = molecule("CH4"), molecule("CH4")
    >>> b.translate([4.0, 0.0, 0.0])                    # methane dimer, C...C = 4 A
    >>> from cage_isomer_builder.utils.dispersion import UFFPairModel, uff_interaction_energy
    >>> model = UFFPairModel(a, [b])
    >>> host, guest = np.arange(5), model.offsets[1] + np.arange(5)
    >>> energy, penalty, overlap = model.score(host, a.positions, guest, b.positions[None],
    ...                                        min_distance=np.zeros((5, 5)))
    >>> bool(np.isclose(energy[0], uff_interaction_energy([a, b]))), float(penalty[0])
    (True, 0.0)
    """
    x = np.sqrt(self.x[source_existing][:, None] * self.x[source_guest][None, :])
    d = np.sqrt(self.d[source_existing][:, None] * self.d[source_guest][None, :])
    k_clash = _CLASH_PENALTY * units.Hartree / units.Bohr ** 2
    return _score_poses_lj(
        np.ascontiguousarray(existing_positions, dtype=float),
        np.ascontiguousarray(guest_positions, dtype=float),
        np.ascontiguousarray(x), np.ascontiguousarray(d),
        np.ascontiguousarray(min_distance, dtype=float), k_clash,
    )

interaction_energy ¶

interaction_energy(fragments)

Lennard-Jones energy (eV) between all pairs of fragments.

Examples:

>>> import numpy as np
>>> from ase.build import molecule
>>> a, b = molecule("CH4"), molecule("CH4")
>>> b.translate([4.0, 0.0, 0.0])                    # methane dimer, C...C = 4 A
>>> from cage_isomer_builder.utils.dispersion import UFFPairModel, uff_interaction_energy
>>> UFFPairModel(a, [b]).interaction_energy([a, b]) == uff_interaction_energy([a, b])
True

Source code in cage_isomer_builder/utils/dispersion.py

def interaction_energy(self, fragments):
    """Lennard-Jones energy (eV) between all pairs of fragments.

    Examples
    --------
    >>> import numpy as np
    >>> from ase.build import molecule
    >>> a, b = molecule("CH4"), molecule("CH4")
    >>> b.translate([4.0, 0.0, 0.0])                    # methane dimer, C...C = 4 A
    >>> from cage_isomer_builder.utils.dispersion import UFFPairModel, uff_interaction_energy
    >>> UFFPairModel(a, [b]).interaction_energy([a, b]) == uff_interaction_energy([a, b])
    True
    """
    return uff_interaction_energy(fragments)

dftd4_available ¶

dftd4_available()

True if the optional dftd4 package can be imported.

Examples:

>>> from cage_isomer_builder.utils.dispersion import dftd4_available
>>> isinstance(dftd4_available(), bool)
True

Source code in cage_isomer_builder/utils/dispersion.py

def dftd4_available():
    """True if the optional dftd4 package can be imported.

    Examples
    --------
    >>> from cage_isomer_builder.utils.dispersion import dftd4_available
    >>> isinstance(dftd4_available(), bool)
    True
    """
    return importlib.util.find_spec("dftd4") is not None

d4_interaction_energy ¶

d4_interaction_energy(fragments, method='pbe')

Full D4 interaction energy (eV), including the three-body term: E(all fragments together) - sum of E(each fragment alone).

Parameters:

Name	Type	Description	Default
`fragments`	`sequence of ase.Atoms`	For example the host followed by each placed guest.	required
`method`	`str or dict`		`"pbe"`

Examples:

>>> import numpy as np
>>> from ase.build import molecule
>>> a, b = molecule("CH4"), molecule("CH4")
>>> b.translate([4.0, 0.0, 0.0])                    # methane dimer, C...C = 4 A
>>> from cage_isomer_builder.utils.dispersion import d4_interaction_energy
>>> far = b.copy()
>>> far.translate([50.0, 0.0, 0.0])
>>> bool(d4_interaction_energy([a, b]) < 0), round(float(d4_interaction_energy([a, far])), 6)
(True, 0.0)

Source code in cage_isomer_builder/utils/dispersion.py

def d4_interaction_energy(fragments, method="pbe"):
    """
    Full D4 interaction energy (eV), including the three-body term:
    E(all fragments together) - sum of E(each fragment alone).

    Parameters
    ----------
    fragments : sequence of ase.Atoms
        For example the host followed by each placed guest.
    method : str or dict, default "pbe"

    Examples
    --------
    >>> import numpy as np
    >>> from ase.build import molecule
    >>> a, b = molecule("CH4"), molecule("CH4")
    >>> b.translate([4.0, 0.0, 0.0])                    # methane dimer, C...C = 4 A
    >>> from cage_isomer_builder.utils.dispersion import d4_interaction_energy
    >>> far = b.copy()
    >>> far.translate([50.0, 0.0, 0.0])
    >>> bool(d4_interaction_energy([a, b]) < 0), round(float(d4_interaction_energy([a, far])), 6)
    (True, 0.0)
    """
    _require_dftd4()
    from dftd4.interface import DampingParam, DispersionModel

    param = DampingParam(**_damping(method))

    def energy(atoms):
        model = DispersionModel(atoms.numbers, atoms.positions / units.Bohr)
        return model.get_dispersion(param, grad=False)["energy"]

    combined = Atoms()
    for fragment in fragments:
        combined += fragment
    total = energy(combined) - sum(energy(f) for f in fragments)
    return total * units.Hartree

uff_interaction_energy ¶

uff_interaction_energy(fragments, repulsion_only=False)

UFF Lennard-Jones interaction energy (eV): the sum over every atom pair that belongs to two different fragments. With repulsion_only, only the repulsive D*(x/r)^12 part.

Examples:

>>> import numpy as np
>>> from ase.build import molecule
>>> a, b = molecule("CH4"), molecule("CH4")
>>> b.translate([4.0, 0.0, 0.0])                    # methane dimer, C...C = 4 A
>>> from cage_isomer_builder.utils.dispersion import uff_interaction_energy
>>> bool(uff_interaction_energy([a, b]) < 0 < uff_interaction_energy([a, b], repulsion_only=True))
True

Source code in cage_isomer_builder/utils/dispersion.py

def uff_interaction_energy(fragments, repulsion_only=False):
    """
    UFF Lennard-Jones interaction energy (eV): the sum over every atom pair
    that belongs to two different fragments. With ``repulsion_only``, only
    the repulsive D*(x/r)^12 part.

    Examples
    --------
    >>> import numpy as np
    >>> from ase.build import molecule
    >>> a, b = molecule("CH4"), molecule("CH4")
    >>> b.translate([4.0, 0.0, 0.0])                    # methane dimer, C...C = 4 A
    >>> from cage_isomer_builder.utils.dispersion import uff_interaction_energy
    >>> bool(uff_interaction_energy([a, b]) < 0 < uff_interaction_energy([a, b], repulsion_only=True))
    True
    """
    total = 0.0
    params = [_uff_element_params(f.numbers) for f in fragments]
    for i in range(len(fragments)):
        for j in range(i + 1, len(fragments)):
            r = np.linalg.norm(
                fragments[i].positions[:, None, :] - fragments[j].positions[None, :, :], axis=2
            )
            x = np.sqrt(params[i][0][:, None] * params[j][0][None, :])
            d = np.sqrt(params[i][1][:, None] * params[j][1][None, :])
            attraction = 0.0 if repulsion_only else 2 * (x / r) ** 6
            total += float(np.sum(d * ((x / r) ** 12 - attraction)))
    return total

pair_model ¶

pair_model(host, guests, dispersion='auto')

The pair model for dispersion:

"auto": D4 with PBE parameters (plus UFF repulsion) if dftd4 is installed, otherwise UFF Lennard-Jones (with a warning).
"uff": UFF Lennard-Jones.
a D4 method name (e.g. "pbe") or a dict of D4 damping parameters.

Examples:

>>> import numpy as np
>>> from ase.build import molecule
>>> a, b = molecule("CH4"), molecule("CH4")
>>> b.translate([4.0, 0.0, 0.0])                    # methane dimer, C...C = 4 A
>>> from cage_isomer_builder.utils.dispersion import pair_model
>>> pair_model(a, [b], dispersion="uff").name
'UFF'

Source code in cage_isomer_builder/utils/dispersion.py

def pair_model(host, guests, dispersion="auto"):
    """
    The pair model for ``dispersion``:

    * ``"auto"``: D4 with PBE parameters (plus UFF repulsion) if dftd4 is
      installed, otherwise UFF Lennard-Jones (with a warning).
    * ``"uff"``: UFF Lennard-Jones.
    * a D4 method name (e.g. ``"pbe"``) or a dict of D4 damping parameters.

    Examples
    --------
    >>> import numpy as np
    >>> from ase.build import molecule
    >>> a, b = molecule("CH4"), molecule("CH4")
    >>> b.translate([4.0, 0.0, 0.0])                    # methane dimer, C...C = 4 A
    >>> from cage_isomer_builder.utils.dispersion import pair_model
    >>> pair_model(a, [b], dispersion="uff").name
    'UFF'
    """
    if isinstance(dispersion, str) and dispersion.lower() == "auto":
        if dftd4_available():
            return D4PairModel(host, guests, "pbe")
        warnings.warn(
            "dftd4 is not installed, so guest sites are scored with the UFF "
            "Lennard-Jones term instead of D4. Install it with "
            "'pip install cage_isomer_builder[d4]' (no Windows wheels), or pass "
            "dispersion='uff' to silence this warning."
        )
        return UFFPairModel(host, guests)
    if isinstance(dispersion, str) and dispersion.lower() == "uff":
        return UFFPairModel(host, guests)
    return D4PairModel(host, guests, dispersion)