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

149

150

151

152

153

154

155

156

157

158

159

160

161

162

163

164

165

166

167

168

169

170

171

172

173

174

175

176

177

178

179

180

181

182

183

184

185

186

187

188

189

190

191

192

193

194

195

196

197

198

199

200

201

202

203

204

205

206

207

208

209

210

211

212

213

214

215

216

217

218

219

220

221

222

223

224

225

226

227

228

229

230

231

232

233

234

235

236

237

238

239

240

241

242

243

244

245

246

247

248

249

250

251

252

253

254

255

""" 

The ``io.shengBTE`` module provides functions for reading and writing data 

files in shengBTE format. 

""" 

 

import numpy as np 

from itertools import product 

from collections import namedtuple 

from ..core.structures import Supercell 

import hiphive 

 

 

def read_shengBTE_fc3(filename, prim, supercell, symprec=1e-5): 

"""Reads third-order force constants file in shengBTE format. 

 

Parameters 

---------- 

filename : str 

input file name 

prim : ase.Atoms 

primitive cell for the force constants 

supercell : ase.Atoms 

supercell onto which to map force constants 

symprec : float 

structural symmetry tolerance 

 

Returns 

------- 

fcs : ForceConstants 

third order force constants for the specified supercell 

""" 

 

raw_sheng = _read_raw_sheng(filename) 

 

sheng = _raw_to_fancy(raw_sheng, prim.cell) 

 

fcs = _sheng_to_fcs(sheng, prim, supercell, symprec) 

 

return fcs 

 

 

def write_shengBTE_fc3(filename, fcs, prim, symprec=1e-5, cutoff=np.inf, 

fc_tol=1e-8): 

"""Writes third-order force constants file in shengBTE format. 

 

Parameters 

----------- 

filename : str 

input file name 

fcs : ForceConstants 

force constants; the supercell associated with this object must be 

based on prim 

prim : ase.Atoms 

primitive configuration (must be equivalent to structure used in the 

shengBTE calculation) 

symprec : float 

structural symmetry tolerance 

cutoff : float 

all atoms in cluster must be within this cutoff 

fc_tol : float 

if the absolute value of the largest entry in a force constant is less 

than fc_tol it will not be written 

""" 

 

sheng = _fcs_to_sheng(fcs, prim, symprec, cutoff, fc_tol) 

 

raw_sheng = _fancy_to_raw(sheng) 

 

_write_raw_sheng(raw_sheng, filename) 

 

 

def _read_raw_sheng(filename): 

""" Read shengBTE fc3 file. 

 

Parameters 

---------- 

filename : str 

input file name 

 

Returns 

------- 

list 

list with shengBTE block entries, where an entry consist of 

[i, j, k, cell_pos2, cell_pos3, fc3] 

""" 

lines_per_fc_block = 32 

 

fc3_data_shengBTE = [] 

with open(filename, 'r') as f: 

lines = f.readlines() 

num_fcs = int(lines[0]) 

lind = 2 

for i in range(1, num_fcs+1): 

# sanity check 

assert int(lines[lind]) == i, (int(lines[lind]), i) 

 

# read cell poses 

cell_pos2 = np.array([float(fld) for fld in lines[lind+1].split()]) 

cell_pos3 = np.array([float(fld) for fld in lines[lind+2].split()]) 

 

# read basis indices 

i, j, k = (int(fld) for fld in lines[lind+3].split()) 

 

# read fc 

fc3_ijk = np.zeros((3, 3, 3)) 

for n in range(27): 

x, y, z = [int(fld) - 1 for fld in lines[lind+4+n].split()[:3]] 

fc3_ijk[x, y, z] = float(lines[lind+4+n].split()[-1]) 

 

# append entry 

entry = [i, j, k, cell_pos2, cell_pos3, fc3_ijk] 

fc3_data_shengBTE.append(entry) 

lind += lines_per_fc_block 

 

return fc3_data_shengBTE 

 

 

_ShengEntry = namedtuple('Entry', ['site_0', 'site_1', 'site_2', 'pos_1', 

'pos_2', 'fc', 'offset_1', 'offset_2']) 

 

 

def _raw_to_fancy(raw_sheng, cell): 

""" 

Converts force constants as read from shengBTE file to the namedtuple 

format defined above (_ShengEntry). 

""" 

sheng = [] 

for raw_entry in raw_sheng: 

p1, p2 = raw_entry[3:5] 

offset_1 = np.linalg.solve(cell.T, p1).round(0).astype(int) 

offset_2 = np.linalg.solve(cell.T, p2).round(0).astype(int) 

entry = _ShengEntry(*(i - 1 for i in raw_entry[:3]), *raw_entry[3:], 

offset_1, offset_2) 

sheng.append(entry) 

return sheng 

 

 

def _fancy_to_raw(sheng): 

""" 

Converts force constants namedtuple format defined above (_ShengEntry) to 

format used for writing shengBTE files. 

""" 

raw_sheng = [] 

for entry in sheng: 

raw_entry = list(entry[:6]) 

raw_entry[0] += 1 

raw_entry[1] += 1 

raw_entry[2] += 1 

raw_sheng.append(raw_entry) 

 

return raw_sheng 

 

 

def _write_raw_sheng(raw_sheng, filename): 

""" See corresponding read function. """ 

 

with open(filename, 'w') as f: 

f.write('{}\n\n'.format(len(raw_sheng))) 

 

for index, fc3_row in enumerate(raw_sheng, start=1): 

i, j, k, cell_pos2, cell_pos3, fc3_ijk = fc3_row 

 

f.write('{:5d}\n'.format(index)) 

 

f.write((3*'{:14.10f} '+'\n').format(*cell_pos2)) 

f.write((3*'{:14.10f} '+'\n').format(*cell_pos3)) 

f.write((3*'{:5d}'+'\n').format(i, j, k)) 

 

for x, y, z in product(range(3), repeat=3): 

f.write((3*' {:}').format(x+1, y+1, z+1)) 

f.write(' {:14.10f}\n'.format(fc3_ijk[x, y, z])) 

f.write('\n') 

 

 

def _fcs_to_sheng(fcs, prim, symprec, cutoff, fc_tol): 

""" Produces a list containing fcs in shengBTE format 

 

prim and fcs.supercell must be aligned. 

""" 

 

supercell = Supercell(fcs.supercell, prim, symprec) 

assert all(fcs.supercell.pbc) and all(prim.pbc) 

 

n_atoms = len(supercell) 

 

D = fcs.supercell.get_all_distances(mic=False, vector=True) 

D_mic = fcs.supercell.get_all_distances(mic=True, vector=True) 

M = np.eye(n_atoms, dtype=bool) 

for i in range(n_atoms): 

for j in range(i + 1, n_atoms): 

M[i, j] = (np.allclose(D[i, j], D_mic[i, j], atol=symprec, rtol=0) 

and np.linalg.norm(D[i, j]) < cutoff) 

M[j, i] = M[i, j] 

 

data = {} 

for a0 in supercell: 

for a1 in supercell: 

if not M[a0.index, a1.index]: 

continue 

for a2 in supercell: 

if not (M[a0.index, a2.index] and M[a1.index, a2.index]): 

continue 

 

offset_1 = np.subtract(a1.offset, a0.offset) 

offset_2 = np.subtract(a2.offset, a0.offset) 

 

sites = (a0.site, a1.site, a2.site) 

 

key = sites + tuple(offset_1) + tuple(offset_2) 

 

ijk = (a0.index, a1.index, a2.index) 

 

fc = fcs[ijk] 

 

if key in data: 

assert np.allclose(data[key], fc, atol=fc_tol) 

else: 

data[key] = fc 

 

sheng = [] 

for k, fc in data.items(): 

if np.max(np.abs(fc)) < fc_tol: 

continue 

offset_1 = k[3:6] 

pos_1 = np.dot(offset_1, prim.cell) 

offset_2 = k[6:9] 

pos_2 = np.dot(offset_2, prim.cell) 

entry = _ShengEntry(*k[:3], pos_1, pos_2, fc, offset_1, offset_2) 

sheng.append(entry) 

 

return sheng 

 

 

def _sheng_to_fcs(sheng, prim, supercell, symprec): 

supercell_map = Supercell(supercell, prim, symprec) 

fc_array = np.zeros((len(supercell),) * 3 + (3,) * 3) 

mapped_clusters = np.zeros((len(supercell),) * 3, dtype=bool) 

 

for atom in supercell_map: 

i = atom.index 

for entry in sheng: 

242 ↛ 243line 242 didn't jump to line 243, because the condition on line 242 was never true if atom.site != entry.site_0: 

continue 

j = supercell_map.index(entry.site_1, entry.offset_1 + atom.offset) 

k = supercell_map.index(entry.site_2, entry.offset_2 + atom.offset) 

ijk = (i, j, k) 

247 ↛ 248line 247 didn't jump to line 248, because the condition on line 247 was never true if mapped_clusters[ijk]: 

raise Exception('Ambiguous force constant.' 

' Supercell is too small') 

fc_array[ijk] = entry.fc 

mapped_clusters[ijk] = True 

 

fcs = hiphive.force_constants.ForceConstants.from_arrays( 

supercell, fc3_array=fc_array) 

return fcs