 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. _

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. _

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 # [-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

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

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