Source code for ase.calculators.harmonic

import numpy as np
from numpy.linalg import eigh, norm, pinv
from scipy.linalg import lstsq  # performs better than numpy.linalg.lstsq

from ase import units
from ase.calculators.calculator import Calculator, BaseCalculator, all_changes
from ase.calculators.calculator import CalculatorSetupError, CalculationFailed


[docs]class HarmonicCalculator(BaseCalculator): """Class for calculations with a Hessian-based harmonic force field. See :class:`HarmonicForceField` and the literature. [1]_ """ implemented_properties = ['energy', 'forces'] def __init__(self, harmonicforcefield): """ Parameters ---------- harmonicforcefield: :class:`HarmonicForceField` Class for calculations with a Hessian-based harmonic force field. """ super().__init__() # parameters have been passed to the force field self.harmonicforcefield = harmonicforcefield def calculate(self, atoms, properties, system_changes): self.atoms = atoms.copy() # for caching of results energy, forces_x = self.harmonicforcefield.get_energy_forces(atoms) self.results['energy'] = energy self.results['forces'] = forces_x
[docs]class HarmonicForceField: def __init__(self, ref_atoms, hessian_x, ref_energy=0.0, get_q_from_x=None, get_jacobian=None, cartesian=True, variable_orientation=False, hessian_limit=0.0, constrained_q=None, rcond=1e-7, zero_thresh=0.0): """ Class that represents a Hessian-based harmonic force field. Energy and forces of this force field are based on the Cartesian Hessian for a local reference configuration, i.e. if desired, on the Hessian matrix transformed to a user-defined coordinate system. The required Hessian has to be passed as an argument, e.g. predetermined numerically via central finite differences in Cartesian coordinates. Note that a potential being harmonic in Cartesian coordinates **x** is not necessarily equivalently harmonic in another coordinate system **q**, e.g. when the transformation between the coordinate systems is non-linear. By default, the force field is evaluated in Cartesian coordinates in which energy and forces are not rotationally and translationally invariant. Systems with variable orientation, require rotationally and translationally invariant calculations for which a set of appropriate coordinates has to be defined. This can be a set of (redundant) internal coordinates (bonds, angles, dihedrals, coordination numbers, ...) or any other user-defined coordinate system. Together with the :class:`HarmonicCalculator` this :class:`HarmonicForceField` can be used to compute Anharmonic Corrections to the Harmonic Approximation. [1]_ Parameters ---------- ref_atoms: :class:`~ase.Atoms` object Reference structure for which energy (``ref_energy``) and Hessian matrix in Cartesian coordinates (``hessian_x``) are provided. hessian_x: numpy array Cartesian Hessian matrix for the reference structure ``ref_atoms``. If a user-defined coordinate system is provided via ``get_q_from_x`` and ``get_jacobian``, the Cartesian Hessian matrix is transformed to the user-defined coordinate system and back to Cartesian coordinates, thereby eliminating rotational and translational traits from the Hessian. The Hessian matrix obtained after this double-transformation is then used as the reference Hessian matrix to evaluate energy and forces for ``cartesian = True``. For ``cartesian = False`` the reference Hessian matrix transformed to the user-defined coordinates is used to compute energy and forces. ref_energy: float Energy of the reference structure ``ref_atoms``, typically in `eV`. get_q_from_x: python function, default: None (Cartesian coordinates) Function that returns a vector of user-defined coordinates **q** for a given :class:`~ase.Atoms` object 'atoms'. The signature should be: :obj:`get_q_from_x(atoms)`. get_jacobian: python function, default: None (Cartesian coordinates) Function that returns the geometric Jacobian matrix of the user-defined coordinates **q** w.r.t. Cartesian coordinates **x** defined as `dq/dx` (Wilson B-matrix) for a given :class:`~ase.Atoms` object 'atoms'. The signature should be: :obj:`get_jacobian(atoms)`. cartesian: bool Set to True to evaluate energy and forces based on the reference Hessian (system harmonic in Cartesian coordinates). Set to False to evaluate energy and forces based on the reference Hessian transformed to user-defined coordinates (system harmonic in user-defined coordinates). hessian_limit: float Reconstruct the reference Hessian matrix with a lower limit for the eigenvalues, typically in `eV/A^2`. Eigenvalues in the interval [``zero_thresh``, ``hessian_limit``] are set to ``hessian_limit`` while the eigenvectors are left untouched. variable_orientation: bool Set to True if the investigated :class:`~ase.Atoms` has got rotational degrees of freedom such that the orientation with respect to ``ref_atoms`` might be different (typically for molecules). Set to False to speed up the calculation when ``cartesian = True``. constrained_q: list A list of indices 'i' of constrained coordinates `q_i` to be projected out from the Hessian matrix (e.g. remove forces along imaginary mode of a transition state). rcond: float Cutoff for singular value decomposition in the computation of the Moore-Penrose pseudo-inverse during transformation of the Hessian matrix. Equivalent to the rcond parameter in scipy.linalg.lstsq. zero_thresh: float Reconstruct the reference Hessian matrix with absolute eigenvalues below this threshold set to zero. """ self.check_input([get_q_from_x, get_jacobian], variable_orientation, cartesian) self.parameters = {'ref_atoms': ref_atoms, 'ref_energy': ref_energy, 'hessian_x': hessian_x, 'hessian_limit': hessian_limit, 'get_q_from_x': get_q_from_x, 'get_jacobian': get_jacobian, 'cartesian': cartesian, 'variable_orientation': variable_orientation, 'constrained_q': constrained_q, 'rcond': rcond, 'zero_thresh': zero_thresh} # set up user-defined coordinate system or Cartesian coordinates self.get_q_from_x = (self.parameters['get_q_from_x'] or (lambda atoms: atoms.get_positions())) self.get_jacobian = (self.parameters['get_jacobian'] or (lambda atoms: np.diagflat(np.ones(3 * len(atoms))))) # reference Cartesian coords. x0; reference user-defined coords. q0 self.x0 = self.parameters['ref_atoms'].get_positions().ravel() self.q0 = self.get_q_from_x(self.parameters['ref_atoms']).ravel() self.setup_reference_hessians() # self._hessian_x and self._hessian_q # store number of zero eigenvalues of G-matrix for redundancy check jac0 = self.get_jacobian(self.parameters['ref_atoms']) Gmat = jac0.T @ jac0 self.Gmat_eigvals, _ = eigh(Gmat) # stored for inspection purposes self.zero_eigvals = len(np.flatnonzero(np.abs(self.Gmat_eigvals) < self.parameters['zero_thresh'])) @staticmethod def check_input(coord_functions, variable_orientation, cartesian): if None in coord_functions: if not all([func is None for func in coord_functions]): msg = ('A user-defined coordinate system requires both ' '`get_q_from_x` and `get_jacobian`.') raise CalculatorSetupError(msg) if variable_orientation: msg = ('The use of `variable_orientation` requires a ' 'user-defined, translationally and rotationally ' 'invariant coordinate system.') raise CalculatorSetupError(msg) if not cartesian: msg = ('A user-defined coordinate system is required for ' 'calculations with cartesian=False.') raise CalculatorSetupError(msg) def setup_reference_hessians(self): """Prepare projector to project out constrained user-defined coordinates **q** from Hessian. Then do transformation to user-defined coordinates and back. Relevant literature: * Peng, C. et al. J. Comput. Chem. 1996, 17 (1), 49-56. * Baker, J. et al. J. Chem. Phys. 1996, 105 (1), 192–212.""" jac0 = self.get_jacobian(self.parameters['ref_atoms']) # Jacobian (dq/dx) jac0 = self.constrain_jac(jac0) # for reference Cartesian coordinates ijac0 = self.get_ijac(jac0, self.parameters['rcond']) self.transform2reference_hessians(jac0, ijac0) # perform projection def constrain_jac(self, jac): """Procedure by Peng, Ayala, Schlegel and Frisch adjusted for redundant coordinates. Peng, C. et al. J. Comput. Chem. 1996, 17 (1), 49–56. """ proj = jac @ jac.T # build non-redundant projector constrained_q = self.parameters['constrained_q'] or [] Cmat = np.zeros(proj.shape) # build projector for constraints Cmat[constrained_q, constrained_q] = 1.0 proj = proj - proj @ Cmat @ pinv(Cmat @ proj @ Cmat) @ Cmat @ proj jac = pinv(jac) @ proj # come back to redundant projector return jac.T def transform2reference_hessians(self, jac0, ijac0): """Transform Cartesian Hessian matrix to user-defined coordinates and back to Cartesian coordinates. For suitable coordinate systems (e.g. internals) this removes rotational and translational degrees of freedom. Furthermore, apply the lower limit to the force constants and reconstruct Hessian matrix.""" hessian_x = self.parameters['hessian_x'] hessian_x = 0.5 * (hessian_x + hessian_x.T) # guarantee symmetry hessian_q = ijac0.T @ hessian_x @ ijac0 # forward transformation hessian_x = jac0.T @ hessian_q @ jac0 # backward transformation hessian_x = 0.5 * (hessian_x + hessian_x.T) # guarantee symmetry w, v = eigh(hessian_x) # rot. and trans. degrees of freedom are removed w[np.abs(w) < self.parameters['zero_thresh']] = 0.0 # noise-cancelling w[(0.0 < w) & # substitute small eigenvalues by lower limit (w < self.parameters['hessian_limit'])] = self.parameters['hessian_limit'] # reconstruct Hessian from new eigenvalues and preserved eigenvectors hessian_x = v @ np.diagflat(w) @ v.T # v.T == inv(v) due to symmetry self._hessian_x = 0.5 * (hessian_x + hessian_x.T) # guarantee symmetry self._hessian_q = ijac0.T @ self._hessian_x @ ijac0 @staticmethod def get_ijac(jac, rcond): # jac is the Wilson B-matrix """Compute Moore-Penrose pseudo-inverse of Wilson B-matrix.""" jac_T = jac.T # btw. direct Jacobian inversion is slow, hence form Gmat Gmat = jac_T @ jac # avoid: numpy.linalg.pinv(Gmat, rcond) @ jac_T ijac = lstsq(Gmat, jac_T, rcond, lapack_driver='gelsy') return ijac[0] # [-1] would be eigenvalues of Gmat def get_energy_forces(self, atoms): """Return a tuple with energy and forces in Cartesian coordinates for a given :class:`~ase.Atoms` object.""" q = self.get_q_from_x(atoms).ravel() if self.parameters['cartesian']: x = atoms.get_positions().ravel() x0 = self.x0 hessian_x = self._hessian_x if self.parameters['variable_orientation']: # determine x0 for present orientation x0 = self.back_transform(x, q, self.q0, atoms.copy()) ref_atoms = atoms.copy() ref_atoms.set_positions(x0.reshape(int(len(x0) / 3), 3)) x0 = ref_atoms.get_positions().ravel() # determine jac0 for present orientation jac0 = self.get_jacobian(ref_atoms) self.check_redundancy(jac0) # check for coordinate failure # determine hessian_x for present orientation hessian_x = jac0.T @ self._hessian_q @ jac0 xdiff = x - x0 forces_x = -hessian_x @ xdiff energy = -0.5 * (forces_x * xdiff).sum() else: jac = self.get_jacobian(atoms) self.check_redundancy(jac) # check for coordinate failure qdiff = q - self.q0 forces_q = -self._hessian_q @ qdiff forces_x = forces_q @ jac energy = -0.5 * (forces_q * qdiff).sum() energy += self.parameters['ref_energy'] forces_x = forces_x.reshape(int(forces_x.size / 3), 3) return energy, forces_x def back_transform(self, x, q, q0, atoms_copy): """Find the right orientation in Cartesian reference coordinates.""" xk = 1 * x qk = 1 * q dq = qk - q0 err = abs(dq).max() count = 0 atoms_copy.set_constraint() # helpful for back-transformation while err > 1e-7: # back-transformation tolerance for convergence count += 1 if count > 99: # maximum number of iterations during back-transf. msg = ('Back-transformation from user-defined to Cartesian ' 'coordinates failed.') raise CalculationFailed(msg) jac = self.get_jacobian(atoms_copy) ijac = self.get_ijac(jac, self.parameters['rcond']) dx = ijac @ dq xk = xk - dx atoms_copy.set_positions(xk.reshape(int(len(xk) / 3), 3)) qk = self.get_q_from_x(atoms_copy).ravel() dq = qk - q0 err = abs(dq).max() return xk def check_redundancy(self, jac): """Compare number of zero eigenvalues of G-matrix to initial number.""" Gmat = jac.T @ jac self.Gmat_eigvals, _ = eigh(Gmat) zero_eigvals = len(np.flatnonzero(np.abs(self.Gmat_eigvals) < self.parameters['zero_thresh'])) if zero_eigvals != self.zero_eigvals: raise CalculationFailed('Suspected coordinate failure: ' f'G-matrix has got {zero_eigvals} ' 'zero eigenvalues, but had ' f'{self.zero_eigvals} during setup') @property def hessian_x(self): return self._hessian_x @property def hessian_q(self): return self._hessian_q
class SpringCalculator(Calculator): """ Spring calculator corresponding to independent oscillators with a fixed spring constant. Energy for an atom is given as E = k / 2 * (r - r_0)**2 where k is the spring constant and, r_0 the ideal positions. Parameters ---------- ideal_positions : array array of the ideal crystal positions k : float spring constant in eV/Angstrom """ implemented_properties = ['forces', 'energy', 'free_energy'] def __init__(self, ideal_positions, k): Calculator.__init__(self) self.ideal_positions = ideal_positions.copy() self.k = k def calculate(self, atoms=None, properties=['energy'], system_changes=all_changes): Calculator.calculate(self, atoms, properties, system_changes) energy, forces = self.compute_energy_and_forces(atoms) self.results['energy'], self.results['forces'] = energy, forces def compute_energy_and_forces(self, atoms): disps = atoms.positions - self.ideal_positions forces = - self.k * disps energy = sum(self.k / 2.0 * norm(disps, axis=1)**2) return energy, forces def get_free_energy(self, T, method='classical'): """Get analytic vibrational free energy for the spring system. Parameters ---------- T : float temperature (K) method : str method for free energy computation; 'classical' or 'QM'. """ F = 0.0 masses, counts = np.unique(self.atoms.get_masses(), return_counts=True) for m, c in zip(masses, counts): F += c * SpringCalculator.compute_Einstein_solid_free_energy(self.k, m, T, method) return F @staticmethod def compute_Einstein_solid_free_energy(k, m, T, method='classical'): """ Get free energy (per atom) for an Einstein crystal. Free energy of a Einstein solid given by classical (1) or QM (2) 1. F_E = 3NkbT log( hw/kbT ) 2. F_E = 3NkbT log( 1-exp(hw/kbT) ) + zeropoint Parameters ----------- k : float spring constant (eV/A^2) m : float mass (grams/mole or AMU) T : float temperature (K) method : str method for free energy computation, classical or QM. Returns -------- float free energy of the Einstein crystal (eV/atom) """ assert method in ['classical', 'QM'] hbar = units._hbar * units.J # eV/s m = m / units.kg # mass kg k = k * units.m**2 / units.J # spring constant J/m2 omega = np.sqrt(k / m) # angular frequency 1/s if method == 'classical': F_einstein = 3 * units.kB * T * np.log(hbar * omega / (units.kB * T)) elif method == 'QM': log_factor = np.log(1.0 - np.exp(-hbar * omega / (units.kB * T))) F_einstein = 3 * units.kB * T * log_factor + 1.5 * hbar * omega return F_einstein