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