lattice

Lattice energy: what it costs to take a crystal apart.

The energy of the crystal, less the energy of its molecules relaxed on their own, per molecule. Negative for anything that holds together.

result = lattice_energy(crystal, calculator)
print(result)            # kJ/mol, with the decomposition
result.structure         # the relaxed crystal

Both sides are relaxed by default, which is what makes the number comparable with a sublimation enthalpy and with other calculations. Relaxing neither gives the interaction energy of the structure as supplied, which is a different and also useful quantity -- see relax.

The result separates two contributions that are often conflated:

E_lattice = E_interaction + E_strain

E_interaction is what the molecules gain by being packed together, measured with each molecule held in the geometry it has in the crystal. E_strain is what they pay to adopt that geometry rather than their relaxed one, and it is positive by construction. A rigid molecule has none; a flexible one can have tens of kJ/mol, and a lattice energy quoted without saying which of the two it is can be out by that much.

LatticeEnergy dataclass

A lattice energy and the parts it is made of.

Attributes:

Name	Type	Description
energy	float	lattice energy per molecule, eV. Negative for a bound crystal.
interaction	float	the part from packing, with molecules held at their crystal geometry, eV per molecule
strain	float	the part paid to adopt the crystal conformation, eV per molecule, positive by construction
crystal_energy	float	total energy of the relaxed unit cell, eV
molecule_energies	list	relaxed gas-phase energy of each unique molecule, eV
frozen_energies	list	energy of each unique molecule at its crystal geometry, eV
multiplicities	list	how many of each unique molecule the cell holds
z	int	molecules per unit cell
z_prime	int	symmetry-unique molecules per cell
structure	object	the relaxed crystal
molecules	list	the relaxed gas-phase molecules
converged	bool	whether every relaxation reached its tolerance. False means the energies are those of structures still on their way downhill, and the difference between two of those is not a lattice energy.
unconverged	list	which relaxations fell short, as readable names
evaluations	int	calculator evaluations used

Source code in chmpy/opt/lattice.py

@dataclass
class LatticeEnergy:
    """A lattice energy and the parts it is made of.

    Attributes:
        energy: lattice energy per molecule, eV. Negative for a bound crystal.
        interaction: the part from packing, with molecules held at their
            crystal geometry, eV per molecule
        strain: the part paid to adopt the crystal conformation, eV per
            molecule, positive by construction
        crystal_energy: total energy of the relaxed unit cell, eV
        molecule_energies: relaxed gas-phase energy of each unique molecule, eV
        frozen_energies: energy of each unique molecule at its crystal
            geometry, eV
        multiplicities: how many of each unique molecule the cell holds
        z: molecules per unit cell
        z_prime: symmetry-unique molecules per cell
        structure: the relaxed crystal
        molecules: the relaxed gas-phase molecules
        converged: whether every relaxation reached its tolerance. False means
            the energies are those of structures still on their way downhill,
            and the difference between two of those is not a lattice energy.
        unconverged: which relaxations fell short, as readable names
        evaluations: calculator evaluations used
    """

    energy: float
    interaction: float
    strain: float
    crystal_energy: float
    molecule_energies: list
    frozen_energies: list
    multiplicities: list
    z: int
    z_prime: int
    structure: object = None
    molecules: list = field(default_factory=list)
    converged: bool = True
    unconverged: list = field(default_factory=list)
    evaluations: int = 0

    @property
    def kj_per_mol(self) -> float:
        "Lattice energy per molecule in kJ/mol, the usual unit for it"
        return self.energy * EV_TO_KJ_PER_MOL

    @property
    def interaction_kj_per_mol(self) -> float:
        return self.interaction * EV_TO_KJ_PER_MOL

    @property
    def strain_kj_per_mol(self) -> float:
        return self.strain * EV_TO_KJ_PER_MOL

    def __repr__(self) -> str:
        warning = "" if self.converged else " NOT CONVERGED"
        return (
            f"<LatticeEnergy{warning} {self.kj_per_mol:.2f} kJ/mol per molecule, "
            f"Z={self.z} Z'={self.z_prime}, {self.evaluations} evaluations>"
        )

    def __str__(self) -> str:
        return "\n".join(
            [
                repr(self),
                f"  interaction {self.interaction_kj_per_mol:9.2f} kJ/mol",
                f"  strain      {self.strain_kj_per_mol:9.2f} kJ/mol",
                f"  lattice     {self.kj_per_mol:9.2f} kJ/mol",
            ]
        )

kj_per_mol property

Lattice energy per molecule in kJ/mol, the usual unit for it

lattice_energy(crystal, calculator, *, relax_crystal=True, relax_molecules=True, box=None, info=None, fmax=0.01, smax=0.05, steps=300, progress=None, **kwargs)

The lattice energy of a molecular crystal.

Parameters:

Name	Type	Description	Default
crystal		a Crystal whose molecules can be separated	required
calculator		a chmpy.calc.Calculator	required
relax_crystal	bool	relax the crystal before taking its energy. With this off the crystal is used as supplied, which is what you want when comparing a set of structures at fixed geometry.	True
relax_molecules	bool	relax each unique molecule in the gas phase. With this off the result is the interaction energy and the strain is zero.	True
box	float | None	put each isolated molecule in a cubic cell of this size in Angstroms, for models that require a periodic system. The default leaves them genuinely isolated.	None
info		model inputs that are not geometry, e.g. a charge or a spin	None
fmax	float	force convergence, eV/A	0.01
smax	float	stress convergence for the crystal, GPa	0.05
steps	int	iteration cap for each relaxation	300
progress		True to print progress, or a callable given a chmpy.opt.progress.Progress per event. The crystal and each molecule are stages; their relaxation steps are nested inside.	None
**kwargs		passed to relax	{}

Returns:

Type	Description
LatticeEnergy	LatticeEnergy

Source code in chmpy/opt/lattice.py

def lattice_energy(
    crystal,
    calculator,
    *,
    relax_crystal: bool = True,
    relax_molecules: bool = True,
    box: float | None = None,
    info=None,
    fmax: float = 0.01,
    smax: float = 0.05,
    steps: int = 300,
    progress=None,
    **kwargs,
) -> LatticeEnergy:
    """The lattice energy of a molecular crystal.

    Args:
        crystal: a `Crystal` whose molecules can be separated
        calculator: a `chmpy.calc.Calculator`
        relax_crystal: relax the crystal before taking its energy. With this
            off the crystal is used as supplied, which is what you want when
            comparing a set of structures at fixed geometry.
        relax_molecules: relax each unique molecule in the gas phase. With this
            off the result is the interaction energy and the strain is zero.
        box: put each isolated molecule in a cubic cell of this size in
            Angstroms, for models that require a periodic system. The default
            leaves them genuinely isolated.
        info: model inputs that are not geometry, e.g. a charge or a spin
        fmax: force convergence, eV/A
        smax: stress convergence for the crystal, GPa
        steps: iteration cap for each relaxation
        progress: True to print progress, or a callable given a
            `chmpy.opt.progress.Progress` per event. The crystal and each
            molecule are stages; their relaxation steps are nested inside.
        **kwargs: passed to `relax`

    Returns:
        LatticeEnergy
    """
    started = calculator.stats.calls

    unique = crystal.symmetry_unique_molecules()
    in_cell = crystal.unit_cell_molecules()
    _check_separation(crystal, unique, in_cell)

    counts = np.bincount(
        [m.properties["asym_mol_idx"] for m in in_cell], minlength=len(unique)
    )
    z = len(in_cell)

    report = reporter(progress)
    stages = (["crystal"] if relax_crystal else []) + (
        [f"molecule {index}" for index in range(len(unique))] if relax_molecules else []
    )

    def announce(stage, message, done=False):
        report(
            "lattice energy",
            stage,
            message,
            index=stages.index(stage),
            total=len(stages),
            done=done,
        )

    unconverged = []
    if relax_crystal:
        announce("crystal", f"relaxing the crystal ({crystal}, Z={z})")
        outcome = relax(
            crystal,
            calculator,
            info=info,
            fmax=fmax,
            smax=smax,
            steps=steps,
            progress=report.nested("lattice energy", "crystal"),
            **kwargs,
        )
        if not outcome.converged:
            unconverged.append("crystal")
            LOG.warning(
                "the crystal relaxation did not converge (fmax %.4g eV/A, "
                "smax %.4g GPa); the lattice energy is that of an unrelaxed "
                "structure",
                outcome.measures.get("fmax", float("nan")),
                outcome.measures.get("smax", float("nan")),
            )
        crystal = outcome.structure
        crystal_energy = outcome.energy
        # the relaxation can change the bonding, and a cell that collapsed is
        # not caught by anything downstream: the arithmetic still works, and
        # the molecules cut out of the wreckage give a confident nonsense
        # number rather than an error
        _check_still_molecular(crystal, z)
        unique = crystal.symmetry_unique_molecules()
        announce("crystal", f"crystal: {outcome!r}", done=True)
    else:
        from chmpy.calc.system import System

        crystal_energy = calculator.energy(System.from_crystal(crystal, **(info or {})))

    frozen, relaxed, molecules = [], [], []
    for index, molecule in enumerate(unique):
        frozen.append(_molecule_energy(molecule, calculator, box, info))
        if relax_molecules:
            stage = f"molecule {index}"
            announce(stage, f"relaxing {stage} ({molecule.molecular_formula})")
            outcome = relax(
                molecule,
                calculator,
                info=info,
                fmax=fmax,
                steps=steps,
                progress=report.nested("lattice energy", stage),
                **kwargs,
            )
            if not outcome.converged:
                unconverged.append(f"molecule {index}")
                LOG.warning(
                    "molecule %d did not relax to fmax %g eV/A; its strain "
                    "energy is a lower bound",
                    index,
                    fmax,
                )
            relaxed.append(outcome.energy)
            molecules.append(outcome.structure)
            announce(stage, f"{stage}: {outcome!r}", done=True)
        else:
            relaxed.append(frozen[-1])
            molecules.append(molecule)

    total_frozen = float(np.dot(counts, frozen))
    total_relaxed = float(np.dot(counts, relaxed))
    interaction = (crystal_energy - total_frozen) / z
    strain = (total_frozen - total_relaxed) / z

    return LatticeEnergy(
        energy=(crystal_energy - total_relaxed) / z,
        interaction=interaction,
        strain=strain,
        crystal_energy=crystal_energy,
        molecule_energies=relaxed,
        frozen_energies=frozen,
        multiplicities=[int(c) for c in counts],
        z=z,
        z_prime=len(unique),
        structure=crystal,
        molecules=molecules,
        converged=not unconverged,
        unconverged=unconverged,
        evaluations=calculator.stats.calls - started,
    )