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