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