Coverage for /builds/debichem-team/python-ase/ase/calculators/harmonic.py : 91.89%
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 ase import units
3from ase.calculators.calculator import Calculator, all_changes
6class SpringCalculator(Calculator):
7 """
8 Spring calculator corresponding to independent oscillators with a fixed
9 spring constant.
12 Energy for an atom is given as
14 E = k / 2 * (r - r_0)**2
16 where k is the spring constant and, r_0 the ideal positions.
19 Parameters
20 ----------
21 ideal_positions : array
22 array of the ideal crystal positions
23 k : float
24 spring constant in eV/Angstrom
25 """
26 implemented_properties = ['forces', 'energy']
28 def __init__(self, ideal_positions, k):
29 Calculator.__init__(self)
30 self.ideal_positions = ideal_positions.copy()
31 self.k = k
33 def calculate(self, atoms=None, properties=['energy'],
34 system_changes=all_changes):
35 Calculator.calculate(self, atoms, properties, system_changes)
36 energy, forces = self.compute_energy_and_forces(atoms)
37 self.results['energy'], self.results['forces'] = energy, forces
39 def compute_energy_and_forces(self, atoms):
40 disps = atoms.positions - self.ideal_positions
41 forces = - self.k * disps
42 energy = sum(self.k / 2.0 * np.linalg.norm(disps, axis=1)**2)
43 return energy, forces
45 def get_free_energy(self, T, method='classical'):
46 """Get analytic vibrational free energy for the spring system.
48 Parameters
49 ----------
50 T : float
51 temperature (K)
52 method : str
53 method for free energy computation; 'classical' or 'QM'.
54 """
55 F = 0.0
56 masses, counts = np.unique(self.atoms.get_masses(), return_counts=True)
57 for m, c in zip(masses, counts):
58 F += c * SpringCalculator.compute_Einstein_solid_free_energy(self.k, m, T, method)
59 return F
61 @staticmethod
62 def compute_Einstein_solid_free_energy(k, m, T, method='classical'):
63 """ Get free energy (per atom) for an Einstein crystal.
65 Free energy of a Einstein solid given by classical (1) or QM (2)
66 1. F_E = 3NkbT log( hw/kbT )
67 2. F_E = 3NkbT log( 1-exp(hw/kbT) ) + zeropoint
69 Parameters
70 -----------
71 k : float
72 spring constant (eV/A^2)
73 m : float
74 mass (grams/mole or AMU)
75 T : float
76 temperature (K)
77 method : str
78 method for free energy computation, classical or QM.
80 Returns
81 --------
82 float
83 free energy of the Einstein crystal (eV/atom)
84 """
85 assert method in ['classical', 'QM']
87 hbar = units._hbar * units.J # eV/s
88 m = m / units.kg # mass kg
89 k = k * units.m**2 / units.J # spring constant J/m2
90 omega = np.sqrt(k / m) # angular frequency 1/s
92 if method == 'classical':
93 F_einstein = 3 * units.kB * T * np.log(hbar * omega / (units.kB * T))
94 elif method == 'QM':
95 log_factor = np.log(1.0 - np.exp(-hbar * omega / (units.kB * T)))
96 F_einstein = 3 * units.kB * T * log_factor + 1.5 * hbar * omega
98 return F_einstein