Hide keyboard shortcuts

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

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

99

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

import numpy as np 

from ase.units import kB 

from hiphive.input_output.logging_tools import logger 

 

logger = logger.getChild(__name__) 

 

 

def generate_phonon_rattled_structures(atoms, fc2, n_structures, temperature): 

""" 

Returns list of phonon-rattled configurations. 

 

Configurations are generated by superimposing harmonic phonon 

eigenmodes with random amplitudes and phase factors consistent 

with a certain temperature. 

 

Let :math:`\\boldsymbol{X}_{ai}` be the phonon modes indexed by atom 

:math:`a` and mode :math:`i`, :math:`\\omega_i` the phonon 

frequencies, and let :math:`0 < Q_i \\leq 1` and :math:`0 \\leq 

U_i < 1` be uniformly random numbers. Then 

 

.. math:: 

 

\\boldsymbol{R}_a 

= \\boldsymbol{R}^0_a 

+ \\left<\\frac{k_B T}{m_a} \\right>^{1/2} 

\\sum_i \\frac{1}{\\omega_i} \\boldsymbol{X}_{ai} 

\\sqrt{-2 \\ln Q_i} \\cos(\\pi \\omega_i U_i) 

 

See: West and Estreicher, Physical Review Letters **96**, 115504 (2006) 

 

Parameters 

---------- 

atoms : ase.Atoms 

prototype structure 

fc2 : numpy.ndarray 

second order force constant matrix, with shape `(3N, 3N)` or 

`(N, N, 3, 3)`. The conversion will be done internally if. 

n_structures : int 

number of structures to generate 

temperature : float 

temperature in Kelvin 

 

Returns 

------- 

list(ase.Atoms) 

generated structures 

""" 

structures = [] 

pr = _PhononRattler(atoms.get_masses(), fc2) 

for _ in range(n_structures): 

atoms_tmp = atoms.copy() 

pr(atoms_tmp, temperature) 

structures.append(atoms_tmp) 

return structures 

 

 

def _phonon_rattle(m_a, T, w2_s, e_sai): 

""" Thermal excitation of phonon modes as described by West and 

Estreicher, Physical Review Letters **96**, 115504 (2006). 

 

_s is a mode index 

_i is a Carteesian index 

_a is an atom index 

 

Parameters 

---------- 

m_a : numpy.ndarray 

masses (N,) 

T : float 

temperature in Kelvin 

w2_s : numpy.ndarray 

the squared frequencies from the igenvalue problem (3*N,) 

e_sai : numpy.ndarray 

polarizations (3*N, N, 3) 

 

Returns 

------- 

displacements : numpy.ndarray 

shape (N, 3) 

""" 

n_modes = 3 * len(m_a) 

prefactor_a = np.sqrt(2 * kB * T / m_a) 

 

# The three modes closest to zero are assumed to be zero, i.e. acoustic sum rules are assumed 

argsort = np.argsort(np.abs(w2_s)) 

w2_gamma = w2_s[argsort][:3] 

e_sai = e_sai[argsort][3:] 

w2_s = w2_s[argsort][3:] 

89 ↛ 90line 89 didn't jump to line 90, because the condition on line 89 was never true if np.any(np.abs(w2_gamma) > 1e-6): 

logger.warning('Acoustic sum rules not enforced, squared frequencies: {}'.format(w2_gamma)) 

 

# treat imaginary modes as real 

93 ↛ 94line 93 didn't jump to line 94, because the condition on line 93 was never true if np.any(w2_s < 0): 

logger.warning('Imaginary modes present') 

w_s = np.sqrt(np.abs(w2_s)) 

 

phases_s = np.random.uniform(0, 2 * np.pi, size=n_modes - 3) 

amplitudes_s = np.sqrt(-np.log(1 - np.random.random(n_modes - 3))) 

 

# prefactors are reshaped in order for numpy broadcast to work 

prefactor_a = prefactor_a.reshape(-1, 1) 

u_ai = prefactor_a * np.tensordot(amplitudes_s * np.cos(phases_s) / w_s, e_sai, (0, 0)) 

return u_ai # displacements 

 

 

class _PhononRattler: 

""" 

Class to be able to conveniently save modes and frequencies needed 

for phonon rattle. 

 

Parameters 

---------- 

masses : numpy.ndarray 

masses (N,) 

force_constants : numpy.ndarray 

second order force constant matrix, with shape `(3N, 3N)` or 

`(N, N, 3, 3)`. The conversion will be done internally if. 

""" 

def __init__(self, masses, force_constants): 

n_atoms = len(masses) 

if len(force_constants.shape) == 4: # assume shape = (n_atoms, n_atoms, 3, 3) 

force_constants = force_constants.transpose(0, 2, 1, 3) 

force_constants = force_constants.reshape(3 * n_atoms, 3 * n_atoms) 

# Now the fc should have shape = (n_atoms * 3, n_atoms * 3) 

# Construct the dynamical matrix 

inv_root_masses = (1 / np.sqrt(masses)).repeat(3) 

D = np.outer(inv_root_masses, inv_root_masses) 

D *= force_constants 

# find modes and energies 

w2_s, e_sai = np.linalg.eigh(D) 

# reshape to get atom index and Cartesian index separate 

e_sai = e_sai.T.reshape(-1, n_atoms, 3) 

self.w2_s = w2_s 

self.e_sai = e_sai 

self.masses = masses 

 

def __call__(self, atoms, T): 

""" rattle atoms by adding displacements 

 

Parameters 

---------- 

atoms : ase.Atoms 

Ideal structure to add displacements to. 

T : float 

temperature in Kelvin 

""" 

u_ai = _phonon_rattle(self.masses, T, self.w2_s, self.e_sai) 

atoms.positions += u_ai