-
Notifications
You must be signed in to change notification settings - Fork 2
/
Copy pathcastep2shengbte-python3.py
320 lines (282 loc) · 13.4 KB
/
castep2shengbte-python3.py
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
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
#
# CASTEP interface to ShengBTE-v1.1
# Written by Genadi Naydenov <genadi.naydenov@york.ac.uk> University of York (2016)
#
# USAGE: ./castep2shengbte.py <seedname> nx ny nz
#
# Input files: <seedname>.castep
#
# The <seedname>.castep file must to be generated by a Phonon calculation.
# Additionally, 'PHONON_WRITE_FORCE_CONSTANTS: true' needs to be set
# in the <seedname>.param file for the phonon run.
#
# Output files: CONTROL, FORCE_CONSTANTS_2ND
#
# The dielectric tensor, epsilon_0, is not calculated during in CASTEP.
# If needed, please use optados to obtain it and check the ShengBTE manual
# for information on how to add it to the CONTROL file.
#
# Please use thirdorder_castep.py to calculate FORCE_CONSTANTS_3RD
# and then proceed to a ShengBTE calculation.
import os.path
import glob
import os
import sys
import math
import numpy as np
def main(argv = None):
print('=========================================================')
print('|| CASTEP 2 ShengBTE Interface ||')
print('|| version 1.3 ||')
print('|| 01 May 2020 ||')
print('||-----------------------------------------------------||')
if argv is None:
argv = sys.argv
if len(argv) < 2:
# Avoid ugly errors
print('|| Usage: castep2shengbte.py <seedname> nx ny nz ||')
print('||-----------------------------------------------------||')
print('|| UNSUCCESSFUL! READ Usage above ||')
print('...')
sys.exit()
# Define <seedname>
prefix = argv[1]
# Help menu, it shows the message and stops the process
help = ['h', '-h','--h', 'help', '-help', '--help']
for i in help:
if i in argv:
print('|| Usage: castep2shengbte.py <seedname> nx ny nz ||')
print('=========================================================')
sys.exit()
"""
Return all the relevant information contained in a .castep file.
"""
castep_file = prefix + '.castep'
castep_file = open(castep_file, 'r')
castep_data = castep_file.readlines()
castep_file.close()
ifc_present = False
for line in castep_data:
if 'Total number of ions in cell' in line:
ions = int(line.split()[7])
elif 'Total number of species in cell' in line:
species = int(line.split()[7])
# elif 'MP grid size for SCF calculation is' in line:
# nx = int(line.split()[7])
# ny = int(line.split()[8])
# nz = int(line.split()[9])
# supercell_size = nx*ny*nz
elif 'Force Constant matrix' in line:
ifc_present = True
if ifc_present == False:
raise ValueError("Force Coonstant matrix not found in .castep."
"\nPlease add 'phonon_write_force_constants: true' to the .param file \nand perform a continuation run.")
nx = int(argv[2])
ny = int(argv[3])
nz = int(argv[4])
supercell_size = nx*ny*nz
# Get unit cell from <seedname>.castep.
lat_vectors = []
for index, line in enumerate(castep_data):
if 'Real Lattice(A)' in line:
start = index + 1
for j in range (3):
lat_vectors.append(
[float(castep_data[start].split()[0]),
float(castep_data[start].split()[1]),
float(castep_data[start].split()[2])])
start += 1
break # avoid double counting
# Species names
species_names = []
for index, line in enumerate(castep_data):
if ' Mass of species in AMU' in line:
start1 = index + 1
for j1 in range (int(species)):
species_names.append(str(castep_data[start1].split()[0]))
start1 += 1
break # avoid double counting
# Species coordinates
positions = []
for index, line in enumerate(castep_data):
if 'Cell Contents' in line:
for i in range(0, ions):
positions.append([str(castep_data[index+10+i].split()[1]),
float(castep_data[index+10+i].split()[3]),
float(castep_data[index+10+i].split()[4]),
float(castep_data[index+10+i].split()[5])])
break # avoid double counting
# Create a list for types which is used in the CONTROL file
types = []
for number_ions in range(len(positions)):
for position, item in enumerate(species_names):
if positions[number_ions][0] in item:
types.append(position+1)
# Born Effective Charges
born_charges = []
for index, line in enumerate(castep_data):
if 'Born Effective Charges' in line:
for i in range(0, ions):
for j in range(3):
if len(castep_data[index+2+i*3+j].split()) >= 4:
born_charges.append([float(castep_data[index+2+i*3+j].split()[2]),
float(castep_data[index+2+i*3+j].split()[3]),
float(castep_data[index+2+i*3+j].split()[4])])
elif len(castep_data[index+2+i*3+j].split())<=3:
born_charges.append([float(castep_data[index+2+i*3+j].split()[0]),
float(castep_data[index+2+i*3+j].split()[1]),
float(castep_data[index+2+i*3+j].split()[2])])
break # avoid double counting
# 2nd order IFC data
ifc_data=dict()
for i in range(supercell_size):
supercell = "Supercell" + "%4s" %(i+1)
ifc_data[supercell]=[]
for n in range(ions*ions):
ifc_data[supercell].append(np.empty((3,3)))
# Define a number which gives the number of lines we have
# per atom per space index in the force constants matrix
l_castep = int(math.ceil(float(ions)/float(2)))
for index, line in enumerate(castep_data):
# Extract information for all supercells
for z in range(supercell_size):
# supercell_3nx3n_matrix used to avoid castep format
# for systems with more than 2 atoms in the unit cell
supercell_3nx3n_matrix = []
supercell = "Supercell" + "%4s" %(z+1)
if supercell in line:
# Go through all lines for a supercell
for n in range(ions):
for i in range(3):
# CASTEP data is stored in 6 columns with ceil(atoms/2) rows.
# line_data stores this data into a single row.
line_data = []
for f in range(l_castep):
new_line = index + (n*3)*l_castep + i*l_castep + f + 1
line_split = castep_data[new_line].split()
if len(line_split)==8:
line_data+=(float(line_split[j]) for j in range(2, 8))
elif len(line_split)<8 and len(line_split)>3:
line_data+=(float(line_split[j]) for j in range(0, 6))
elif len(line_split)<4:
line_data+=(float(line_split[j]) for j in range(0, 3))
supercell_3nx3n_matrix += [line_data]
for m in range(ions):
ifc_data[supercell][ions*n+m][i]=[float(supercell_3nx3n_matrix[3*n+i][j])
for j in range(m*3, 3*m+3)]
"""
All of the necessary information is loaded into the dictionary ifc_data.
It would be used from now on to write out all permutations.
"""
# Required format in the FORCE_CONSTANTS_2ND file:
# line 1: nx*xy*nz*ions
# line 2: atom1 atom2
# lines 3,4,5: IFC matrix for atom 1 and 2
#
# Permutaions use the following pattern: Pick atom 1 and write out its interations with all
# other atoms. For example, in a 5x5x5 supercell with 2 atoms, e.g. NaCl in the unit cell,
# the interaction between atoms 1 and 2# is the interaction between Na in cell (1,1,1) and
# Na in (2,1,1). 1-125 corresponds to Na(1,1,1) - Na(5,5,5), respectively. Therefore, the
# interaction between Na and Cl in (1,1,1) is given by the pair 1-126. The supercell is
# constructured by increasing x first, then y and then z.
ion_reset=0
indeces_array = []
for every_atom in range(0,ions*nx*ny*nz):
# pick 1 ion, go through all cells and then
# reset the number to 1 for the next ion
ion_reset+=1
if ion_reset > nx*ny*nz:
ion_reset=1
divnx, dummy1 = divmod(every_atom, nx)
divnxny, dummy2 = divmod(every_atom, nx*ny)
for every_ion in range(0,ions):
for z in range(0,nz):
counter_z = z + divnxny + every_ion*nz
dummy3, modnz = divmod(counter_z, nz)
for y in range(0,ny):
counter_y = y + divnx + z*ny + every_ion*nz*ny
dummy4, modny = divmod(counter_y, ny)
for x in range(0,nx):
counter_x = (x+ion_reset) + y*nx + z*ny*nx + every_ion*nz*ny*nx
dummy5, modnx = divmod(counter_x-1,nx)
indeces_array.append([int(every_atom+1),
int((modnx+1+nx*modny+modnz*ny*nx + every_ion*nz*ny*nx))])
#keylist = ifc_data.keys()
#keylist.sort()
keylist = sorted(ifc_data)
ifc_for_sheng=dict()
i=0
# go through all columns
for x in range(0,ions):
# do the process for all copies of atom (x)
# [just keep the same column and print its data again]
for lp in range(0,(nx*ny*nz)):
# extract data from a column (given by x) for the interaction between
# atom N with a given atom (defined by x above)
for y in range(0,ions):
# go through all supercells
for key in keylist:
i+=1
# here [i-1][0] and [i-1][1] define the direction of the forces
ion_couple = '%6s' %indeces_array[i-1][0] + '%6s' %indeces_array[i-1][1]
ifc_for_sheng[ion_couple]= [ifc_data[key][ions*x+y]]
# Create the 2nd order file
f = open('FORCE_CONSTANTS_2ND', 'w')
f_struct = '%6s' %str(nx*ny*nz*ions) + '\n'
for key1 in sorted(ifc_for_sheng):
f_struct += str(key1) + '\n'
for m in range(3):
for l in range (3):
f_struct += '%22s' %str("{0:.15f}".format(ifc_for_sheng[key1][0][m][l]))
f_struct += '\n'
f.write(f_struct)
f.close()
# Create the control file
f = open('CONTROL', 'w')
f_control = '&allocations' +'\n'
f_control += '\tnelements=' + '%s,' %str(species) + '\n'
f_control += '\tnatoms=' + '%s,' %str(ions) + '\n'
f_control += '\tngrid(:)=' + '%s ' %str(nx) + '%s ' %str(ny) +'%s' %str(nz) +'\n'
f_control += '&end' +'\n'
f_control += '&crystal' +'\n'
f_control += '\tlattvec(:,1)=' + '%s ' %str(lat_vectors[0][0]/10) + '%s ' %str(lat_vectors[0][1]/10)+ '%s,' %str(lat_vectors[0][2]/10) + '\n'
f_control += '\tlattvec(:,2)=' + '%s ' %str(lat_vectors[1][0]/10) + '%s ' %str(lat_vectors[1][1]/10)+ '%s,' %str(lat_vectors[1][2]/10) + '\n'
f_control += '\tlattvec(:,3)=' + '%s ' %str(lat_vectors[2][0]/10) + '%s ' %str(lat_vectors[2][1]/10)+ '%s,' %str(lat_vectors[2][2]/10) + '\n'
f_control += '\telements='
for el in range(len(species_names)):
f_control +='"%s" ' %str(species_names[el])
f_control += '\n'
f_control += '\ttypes='
for el_num in range(len(types)):
f_control += ' %s' %str(types[el_num])
f_control += ',\n'
for atoms in range(len(positions)):
f_control += '\tpositions(:,%s)=' %str(atoms+1)
for i in range(1,4):
f_control += ' %s' %str("{0:.8f}".format(positions[atoms][i]))
f_control += ',\n'
if len(born_charges)!=0:
for atoms in range(len(positions)):
for i in range(1,4):
f_control += '\tborn(:,%s,' %str(i) + '%s)=' %str(atoms+1)
for j in range(3):
f_control += ' %s' %str("{0:.5f}".format(born_charges[atoms*3 + (i-1)][j]))
f_control += ',\n'
f_control += '\tscell(:)=' + '%s ' %str(nx) + '%s ' %str(ny) +'%s' %str(nz) +'\n'
f_control += '&end' +'\n'
f_control += '¶meters' +'\n'
f_control += '\tT=300.'+'\n'
f_control += '\tscalebroad=1.0'+'\n'
f_control += '&end' +'\n'
f_control += '&flags' +'\n'
f_control += '\tnonanalytic=.TRUE.'+'\n'
f_control += '&end' +'\n'
f.write(f_control)
f.close()
print('|| DONE ||')
print('=========================================================')
if __name__ == "__main__":
import sys
sys.exit(main())