Coverage for /builds/ase/ase/ase/calculators/harmonic.py: 97.56%
164 statements
« prev ^ index » next coverage.py v7.5.3, created at 2025-08-02 00:12 +0000
1# fmt: off
3import numpy as np
4from numpy.linalg import eigh, norm, pinv
5from scipy.linalg import lstsq # performs better than numpy.linalg.lstsq
7from ase import units
8from ase.calculators.calculator import (
9 BaseCalculator,
10 CalculationFailed,
11 Calculator,
12 CalculatorSetupError,
13 all_changes,
14)
17class HarmonicCalculator(BaseCalculator):
18 """Class for calculations with a Hessian-based harmonic force field.
20 See :class:`HarmonicForceField` and the literature. [1]_
21 """
23 implemented_properties = ['energy', 'forces']
25 def __init__(self, harmonicforcefield):
26 """
27 Parameters
28 ----------
29 harmonicforcefield: :class:`HarmonicForceField`
30 Class for calculations with a Hessian-based harmonic force field.
31 """
32 super().__init__() # parameters have been passed to the force field
33 self.harmonicforcefield = harmonicforcefield
35 def calculate(self, atoms, properties, system_changes):
36 self.atoms = atoms.copy() # for caching of results
37 energy, forces_x = self.harmonicforcefield.get_energy_forces(atoms)
38 self.results['energy'] = energy
39 self.results['forces'] = forces_x
42class HarmonicForceField:
43 def __init__(self, ref_atoms, hessian_x, ref_energy=0.0, get_q_from_x=None,
44 get_jacobian=None, cartesian=True, variable_orientation=False,
45 hessian_limit=0.0, constrained_q=None, rcond=1e-7,
46 zero_thresh=0.0):
47 """Class that represents a Hessian-based harmonic force field.
49 Energy and forces of this force field are based on the
50 Cartesian Hessian for a local reference configuration, i.e. if
51 desired, on the Hessian matrix transformed to a user-defined
52 coordinate system. The required Hessian has to be passed as
53 an argument, e.g. predetermined numerically via central finite
54 differences in Cartesian coordinates. Note that a potential
55 being harmonic in Cartesian coordinates **x** is not
56 necessarily equivalently harmonic in another coordinate system
57 **q**, e.g. when the transformation between the coordinate
58 systems is non-linear. By default, the force field is
59 evaluated in Cartesian coordinates in which energy and forces
60 are not rotationally and translationally invariant. Systems
61 with variable orientation, require rotationally and
62 translationally invariant calculations for which a set of
63 appropriate coordinates has to be defined. This can be a set
64 of (redundant) internal coordinates (bonds, angles, dihedrals,
65 coordination numbers, ...) or any other user-defined
66 coordinate system.
68 Together with the :class:`HarmonicCalculator` this
69 :class:`HarmonicForceField` can be used to compute
70 Anharmonic Corrections to the Harmonic Approximation. [1]_
72 Parameters
73 ----------
74 ref_atoms: :class:`~ase.Atoms` object
75 Reference structure for which energy (``ref_energy``) and Hessian
76 matrix in Cartesian coordinates (``hessian_x``) are provided.
78 hessian_x: numpy array
79 Cartesian Hessian matrix for the reference structure ``ref_atoms``.
80 If a user-defined coordinate system is provided via
81 ``get_q_from_x`` and ``get_jacobian``, the Cartesian Hessian matrix
82 is transformed to the user-defined coordinate system and back to
83 Cartesian coordinates, thereby eliminating rotational and
84 translational traits from the Hessian. The Hessian matrix
85 obtained after this double-transformation is then used as
86 the reference Hessian matrix to evaluate energy and forces for
87 ``cartesian = True``. For ``cartesian = False`` the reference
88 Hessian matrix transformed to the user-defined coordinates is used
89 to compute energy and forces.
91 ref_energy: float
92 Energy of the reference structure ``ref_atoms``, typically in `eV`.
94 get_q_from_x: python function, default: None (Cartesian coordinates)
95 Function that returns a vector of user-defined coordinates **q** for
96 a given :class:`~ase.Atoms` object 'atoms'. The signature should be:
97 :obj:`get_q_from_x(atoms)`.
99 get_jacobian: python function, default: None (Cartesian coordinates)
100 Function that returns the geometric Jacobian matrix of the
101 user-defined coordinates **q** w.r.t. Cartesian coordinates **x**
102 defined as `dq/dx` (Wilson B-matrix) for a given :class:`~ase.Atoms`
103 object 'atoms'. The signature should be: :obj:`get_jacobian(atoms)`.
105 cartesian: bool
106 Set to True to evaluate energy and forces based on the reference
107 Hessian (system harmonic in Cartesian coordinates).
108 Set to False to evaluate energy and forces based on the reference
109 Hessian transformed to user-defined coordinates (system harmonic in
110 user-defined coordinates).
112 hessian_limit: float
113 Reconstruct the reference Hessian matrix with a lower limit for the
114 eigenvalues, typically in `eV/A^2`. Eigenvalues in the interval
115 [``zero_thresh``, ``hessian_limit``] are set to ``hessian_limit``
116 while the eigenvectors are left untouched.
118 variable_orientation: bool
119 Set to True if the investigated :class:`~ase.Atoms` has got
120 rotational degrees of freedom such that the orientation with respect
121 to ``ref_atoms`` might be different (typically for molecules).
122 Set to False to speed up the calculation when ``cartesian = True``.
124 constrained_q: list
125 A list of indices 'i' of constrained coordinates `q_i` to be
126 projected out from the Hessian matrix
127 (e.g. remove forces along imaginary mode of a transition state).
129 rcond: float
130 Cutoff for singular value decomposition in the computation of the
131 Moore-Penrose pseudo-inverse during transformation of the Hessian
132 matrix. Equivalent to the rcond parameter in scipy.linalg.lstsq.
134 zero_thresh: float
135 Reconstruct the reference Hessian matrix with absolute eigenvalues
136 below this threshold set to zero.
138 """
139 self.check_input([get_q_from_x, get_jacobian],
140 variable_orientation, cartesian)
142 self.parameters = {'ref_atoms': ref_atoms,
143 'ref_energy': ref_energy,
144 'hessian_x': hessian_x,
145 'hessian_limit': hessian_limit,
146 'get_q_from_x': get_q_from_x,
147 'get_jacobian': get_jacobian,
148 'cartesian': cartesian,
149 'variable_orientation': variable_orientation,
150 'constrained_q': constrained_q,
151 'rcond': rcond,
152 'zero_thresh': zero_thresh}
154 # set up user-defined coordinate system or Cartesian coordinates
155 self.get_q_from_x = (self.parameters['get_q_from_x'] or
156 (lambda atoms: atoms.get_positions()))
157 self.get_jacobian = (self.parameters['get_jacobian'] or
158 (lambda atoms: np.diagflat(
159 np.ones(3 * len(atoms)))))
161 # reference Cartesian coords. x0; reference user-defined coords. q0
162 self.x0 = self.parameters['ref_atoms'].get_positions().ravel()
163 self.q0 = self.get_q_from_x(self.parameters['ref_atoms']).ravel()
164 self.setup_reference_hessians() # self._hessian_x and self._hessian_q
166 # store number of zero eigenvalues of G-matrix for redundancy check
167 jac0 = self.get_jacobian(self.parameters['ref_atoms'])
168 Gmat = jac0.T @ jac0
169 self.Gmat_eigvals, _ = eigh(Gmat) # stored for inspection purposes
170 self.zero_eigvals = len(np.flatnonzero(np.abs(self.Gmat_eigvals) <
171 self.parameters['zero_thresh']))
173 @staticmethod
174 def check_input(coord_functions, variable_orientation, cartesian):
175 if None in coord_functions:
176 if not all(func is None for func in coord_functions):
177 msg = ('A user-defined coordinate system requires both '
178 '`get_q_from_x` and `get_jacobian`.')
179 raise CalculatorSetupError(msg)
180 if variable_orientation:
181 msg = ('The use of `variable_orientation` requires a '
182 'user-defined, translationally and rotationally '
183 'invariant coordinate system.')
184 raise CalculatorSetupError(msg)
185 if not cartesian:
186 msg = ('A user-defined coordinate system is required for '
187 'calculations with cartesian=False.')
188 raise CalculatorSetupError(msg)
190 def setup_reference_hessians(self):
191 """Prepare projector to project out constrained user-defined coordinates
192 **q** from Hessian. Then do transformation to user-defined coordinates
193 and back. Relevant literature:
194 * Peng, C. et al. J. Comput. Chem. 1996, 17 (1), 49-56.
195 * Baker, J. et al. J. Chem. Phys. 1996, 105 (1), 192–212."""
196 jac0 = self.get_jacobian(
197 self.parameters['ref_atoms']) # Jacobian (dq/dx)
198 jac0 = self.constrain_jac(jac0) # for reference Cartesian coordinates
199 ijac0 = self.get_ijac(jac0, self.parameters['rcond'])
200 self.transform2reference_hessians(jac0, ijac0) # perform projection
202 def constrain_jac(self, jac):
203 """Procedure by Peng, Ayala, Schlegel and Frisch adjusted for redundant
204 coordinates.
205 Peng, C. et al. J. Comput. Chem. 1996, 17 (1), 49–56.
206 """
207 proj = jac @ jac.T # build non-redundant projector
208 constrained_q = self.parameters['constrained_q'] or []
209 Cmat = np.zeros(proj.shape) # build projector for constraints
210 Cmat[constrained_q, constrained_q] = 1.0
211 proj = proj - proj @ Cmat @ pinv(Cmat @ proj @ Cmat) @ Cmat @ proj
212 jac = pinv(jac) @ proj # come back to redundant projector
213 return jac.T
215 def transform2reference_hessians(self, jac0, ijac0):
216 """Transform Cartesian Hessian matrix to user-defined coordinates
217 and back to Cartesian coordinates. For suitable coordinate systems
218 (e.g. internals) this removes rotational and translational degrees of
219 freedom. Furthermore, apply the lower limit to the force constants
220 and reconstruct Hessian matrix."""
221 hessian_x = self.parameters['hessian_x']
222 hessian_x = 0.5 * (hessian_x + hessian_x.T) # guarantee symmetry
223 hessian_q = ijac0.T @ hessian_x @ ijac0 # forward transformation
224 hessian_x = jac0.T @ hessian_q @ jac0 # backward transformation
225 hessian_x = 0.5 * (hessian_x + hessian_x.T) # guarantee symmetry
226 w, v = eigh(hessian_x) # rot. and trans. degrees of freedom are removed
227 w[np.abs(w) < self.parameters['zero_thresh']] = 0.0 # noise-cancelling
228 w[(w > 0.0) & (w < self.parameters['hessian_limit'])] = self.parameters[
229 'hessian_limit'
230 ]
231 # reconstruct Hessian from new eigenvalues and preserved eigenvectors
232 hessian_x = v @ np.diagflat(w) @ v.T # v.T == inv(v) due to symmetry
233 self._hessian_x = 0.5 * (hessian_x + hessian_x.T) # guarantee symmetry
234 self._hessian_q = ijac0.T @ self._hessian_x @ ijac0
236 @staticmethod
237 def get_ijac(jac, rcond): # jac is the Wilson B-matrix
238 """Compute Moore-Penrose pseudo-inverse of Wilson B-matrix."""
239 jac_T = jac.T # btw. direct Jacobian inversion is slow, hence form Gmat
240 Gmat = jac_T @ jac # avoid: numpy.linalg.pinv(Gmat, rcond) @ jac_T
241 ijac = lstsq(Gmat, jac_T, rcond, lapack_driver='gelsy')
242 return ijac[0] # [-1] would be eigenvalues of Gmat
244 def get_energy_forces(self, atoms):
245 """Return a tuple with energy and forces in Cartesian coordinates for
246 a given :class:`~ase.Atoms` object."""
247 q = self.get_q_from_x(atoms).ravel()
249 if self.parameters['cartesian']:
250 x = atoms.get_positions().ravel()
251 x0 = self.x0
252 hessian_x = self._hessian_x
254 if self.parameters['variable_orientation']:
255 # determine x0 for present orientation
256 x0 = self.back_transform(x, q, self.q0, atoms.copy())
257 ref_atoms = atoms.copy()
258 ref_atoms.set_positions(x0.reshape(int(len(x0) / 3), 3))
259 x0 = ref_atoms.get_positions().ravel()
260 # determine jac0 for present orientation
261 jac0 = self.get_jacobian(ref_atoms)
262 self.check_redundancy(jac0) # check for coordinate failure
263 # determine hessian_x for present orientation
264 hessian_x = jac0.T @ self._hessian_q @ jac0
266 xdiff = x - x0
267 forces_x = -hessian_x @ xdiff
268 energy = -0.5 * (forces_x * xdiff).sum()
270 else:
271 jac = self.get_jacobian(atoms)
272 self.check_redundancy(jac) # check for coordinate failure
273 qdiff = q - self.q0
274 forces_q = -self._hessian_q @ qdiff
275 forces_x = forces_q @ jac
276 energy = -0.5 * (forces_q * qdiff).sum()
278 energy += self.parameters['ref_energy']
279 forces_x = forces_x.reshape(int(forces_x.size / 3), 3)
280 return energy, forces_x
282 def back_transform(self, x, q, q0, atoms_copy):
283 """Find the right orientation in Cartesian reference coordinates."""
284 xk = 1 * x
285 qk = 1 * q
286 dq = qk - q0
287 err = abs(dq).max()
288 count = 0
289 atoms_copy.set_constraint() # helpful for back-transformation
290 while err > 1e-7: # back-transformation tolerance for convergence
291 count += 1
292 if count > 99: # maximum number of iterations during back-transf.
293 msg = ('Back-transformation from user-defined to Cartesian '
294 'coordinates failed.')
295 raise CalculationFailed(msg)
296 jac = self.get_jacobian(atoms_copy)
297 ijac = self.get_ijac(jac, self.parameters['rcond'])
298 dx = ijac @ dq
299 xk = xk - dx
300 atoms_copy.set_positions(xk.reshape(int(len(xk) / 3), 3))
301 qk = self.get_q_from_x(atoms_copy).ravel()
302 dq = qk - q0
303 err = abs(dq).max()
304 return xk
306 def check_redundancy(self, jac):
307 """Compare number of zero eigenvalues of G-matrix to initial number."""
308 Gmat = jac.T @ jac
309 self.Gmat_eigvals, _ = eigh(Gmat)
310 zero_eigvals = len(np.flatnonzero(np.abs(self.Gmat_eigvals) <
311 self.parameters['zero_thresh']))
312 if zero_eigvals != self.zero_eigvals:
313 raise CalculationFailed('Suspected coordinate failure: '
314 f'G-matrix has got {zero_eigvals} '
315 'zero eigenvalues, but had '
316 f'{self.zero_eigvals} during setup')
318 @property
319 def hessian_x(self):
320 return self._hessian_x
322 @property
323 def hessian_q(self):
324 return self._hessian_q
327class SpringCalculator(Calculator):
328 """
329 Spring calculator corresponding to independent oscillators with a fixed
330 spring constant.
333 Energy for an atom is given as
335 E = k / 2 * (r - r_0)**2
337 where k is the spring constant and, r_0 the ideal positions.
340 Parameters
341 ----------
342 ideal_positions : array
343 array of the ideal crystal positions
344 k : float
345 spring constant in eV/Angstrom
346 """
347 implemented_properties = ['forces', 'energy', 'free_energy']
349 def __init__(self, ideal_positions, k):
350 Calculator.__init__(self)
351 self.ideal_positions = ideal_positions.copy()
352 self.k = k
354 def calculate(self, atoms=None, properties=['energy'],
355 system_changes=all_changes):
356 Calculator.calculate(self, atoms, properties, system_changes)
357 energy, forces = self.compute_energy_and_forces(atoms)
358 self.results['energy'], self.results['forces'] = energy, forces
360 def compute_energy_and_forces(self, atoms):
361 disps = atoms.positions - self.ideal_positions
362 forces = - self.k * disps
363 energy = sum(self.k / 2.0 * norm(disps, axis=1)**2)
364 return energy, forces
366 def get_free_energy(self, T, method='classical'):
367 """Get analytic vibrational free energy for the spring system.
369 Parameters
370 ----------
371 T : float
372 temperature (K)
373 method : str
374 method for free energy computation; 'classical' or 'QM'.
375 """
376 F = 0.0
377 masses, counts = np.unique(self.atoms.get_masses(), return_counts=True)
378 for m, c in zip(masses, counts):
379 F += c * \
380 SpringCalculator.compute_Einstein_solid_free_energy(
381 self.k, m, T, method)
382 return F
384 @staticmethod
385 def compute_Einstein_solid_free_energy(k, m, T, method='classical'):
386 """ Get free energy (per atom) for an Einstein crystal.
388 Free energy of a Einstein solid given by classical (1) or QM (2)
389 1. F_E = 3NkbT log( hw/kbT )
390 2. F_E = 3NkbT log( 1-exp(hw/kbT) ) + zeropoint
392 Parameters
393 -----------
394 k : float
395 spring constant (eV/A^2)
396 m : float
397 mass (grams/mole or AMU)
398 T : float
399 temperature (K)
400 method : str
401 method for free energy computation, classical or QM.
403 Returns
404 --------
405 float
406 free energy of the Einstein crystal (eV/atom)
407 """
408 assert method in ['classical', 'QM']
410 hbar = units._hbar * units.J # eV/s
411 m = m / units.kg # mass kg
412 k = k * units.m**2 / units.J # spring constant J/m2
413 omega = np.sqrt(k / m) # angular frequency 1/s
415 if method == 'classical':
416 F_einstein = 3 * units.kB * T * \
417 np.log(hbar * omega / (units.kB * T))
418 elif method == 'QM':
419 log_factor = np.log(1.0 - np.exp(-hbar * omega / (units.kB * T)))
420 F_einstein = 3 * units.kB * T * log_factor + 1.5 * hbar * omega
422 return F_einstein