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