-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathstructure2input.py
192 lines (189 loc) · 7.31 KB
/
structure2input.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
import pymatgen
import seekpath
import numpy
import math
from pymatgen.symmetry.analyzer import SpacegroupAnalyzer
from write_openmx import write_openmx
from write_pwx import write_pwx
from write_pp import write_pp
from write_ph import write_ph
from write_wannier import write_wannier
from write_sh import write_sh
# from write_hilapw import write_hilapw
from pymatgen.core.periodic_table import get_el_sp
def structure2input(structure, dk_path, dq_grid, pseudo_kind, host, rel):
if pseudo_kind == "sg15":
if rel:
from sg15_rel import pseudo_dict, ecutwfc_dict, ecutrho_dict, valence_dict, atomwfc_dict
else:
from sg15 import pseudo_dict, ecutwfc_dict, ecutrho_dict, valence_dict, atomwfc_dict
elif pseudo_kind == "pslibrary":
if rel:
from pslibrary_rel import pseudo_dict, ecutwfc_dict, ecutrho_dict, valence_dict, atomwfc_dict
else:
from pslibrary import pseudo_dict, ecutwfc_dict, ecutrho_dict, valence_dict, atomwfc_dict
elif pseudo_kind == "sssp":
from sssp import pseudo_dict, ecutwfc_dict, ecutrho_dict, valence_dict, atomwfc_dict
elif pseudo_kind == "ssspsol":
from ssspsol import pseudo_dict, ecutwfc_dict, ecutrho_dict, valence_dict, atomwfc_dict
elif pseudo_kind == "sssp_us":
from sssp_us import pseudo_dict, ecutwfc_dict, ecutrho_dict, valence_dict, atomwfc_dict
else:
from sssp import pseudo_dict, ecutwfc_dict, ecutrho_dict, valence_dict, atomwfc_dict
print("Unsupported pseudo potential library :", pseudo_kind)
exit(1)
#
# Band path and primitive lattice
#
frac_coord2 = numpy.array(structure.frac_coords)
for ipos in range(len(frac_coord2)):
for iaxis in range(3):
coord3 = frac_coord2[ipos, iaxis] * 6.0
if abs(round(coord3) - coord3) < 0.001:
frac_coord2[ipos, iaxis] = float(round(coord3)) / 6.0
#
skp = seekpath.get_path((structure.lattice.matrix, frac_coord2,
[pymatgen.core.Element(str(spc)).number for spc in structure.species]))
#
# Lattice information
#
avec = skp["primitive_lattice"]
bvec = skp["reciprocal_primitive_lattice"]
atom = [str(get_el_sp(iat)) for iat in skp["primitive_types"]]
typ = sorted(set(atom))
print("Bravais lattice : ", skp["bravais_lattice"])
print("Space group : ", skp['spacegroup_international'])
bz_volume = abs(- bvec[0][2]*bvec[1][1]*bvec[2][0]
+ bvec[0][1]*bvec[1][2]*bvec[2][0]
+ bvec[0][2]*bvec[1][0]*bvec[2][1]
- bvec[0][0]*bvec[1][2]*bvec[2][1]
- bvec[0][1]*bvec[1][0]*bvec[2][2]
+ bvec[0][0]*bvec[1][1]*bvec[2][2])
#
# WFC and Rho cutoff
#
ecutwfc = 0.0
ecutrho = 0.0
for ityp in typ:
if ecutwfc < ecutwfc_dict[str(ityp)]:
ecutwfc = ecutwfc_dict[str(ityp)]
if ecutrho < ecutrho_dict[str(ityp)]:
ecutrho = ecutrho_dict[str(ityp)]
numpw = 4.0*math.pi*math.sqrt(ecutwfc)**3/3.0/bz_volume
print("Estimated number of PW (WFC) :", numpw)
#
# k and q grid
# the number of grid proportional to the Height of b
# b_i * a_i / |a_i| = 2pi / |a_i| (a_i is perpendicular to other b's)
#
nq = numpy.zeros(3, numpy.int_)
for ii in range(3):
norm = numpy.sqrt(numpy.dot(avec[ii][:], avec[ii][:]))
nq[ii] = round(2.0 * numpy.pi / norm / dq_grid)
if nq[ii] == 0:
nq[ii] = 1
print("Coarse grid : ", nq[0], nq[1], nq[2])
#
# Get explicit kpath applicable to bands.x and plotband.x
#
kpath = []
nkpath = []
#
dk_jump = 1.0e10
for ipath in range(len(skp['path'])):
dk = numpy.array(skp['point_coords'][skp['path'][ipath][1]])\
- numpy.array(skp['point_coords'][skp['path'][ipath][0]])
dk = numpy.dot(dk, bvec)
dknorm = math.sqrt(numpy.dot(dk, dk))
nkpath0 = max(2, int(dknorm / dk_path))
nkpath.append(nkpath0)
for ik in range(nkpath0):
xkpath = ik / nkpath0
kpath0 = numpy.array(skp['point_coords'][skp['path'][ipath][0]]) * (1.0 - xkpath) \
+ numpy.array(skp['point_coords'][skp['path'][ipath][1]]) * xkpath
kpath.append(kpath0)
#
# jump case
#
if ipath < len(skp['path']) - 1 and skp['path'][ipath][1] != skp['path'][ipath+1][0]:
kpath0 = numpy.array(skp['point_coords'][skp['path'][ipath][1]])
kpath1 = numpy.array(skp['point_coords'][skp['path'][ipath+1][0]])
dk = kpath1 - kpath0
dk = numpy.dot(dk, bvec)
dk_jump = min(dk_jump, math.sqrt(numpy.dot(dk, dk)))
kpath.append(kpath0)
#
# If the jump is relatively small, shift the last point closer to the end
#
dk = numpy.array(kpath[1]) - numpy.array(kpath[0])
dk = numpy.dot(dk, bvec)
dknorm = math.sqrt(numpy.dot(dk, dk))
if dk_jump < 10.0 * dknorm:
xkpath = dk_jump * 0.09
for i in range(3):
kpath[1][i] = kpath[0][i] + (kpath[1][i] - kpath[0][i])*xkpath/dknorm
#
# Last point
#
kpath.append(numpy.array(skp['point_coords'][skp['path'][len(skp['path']) - 1][1]]))
#
# Band path
#
print("Band path")
for ipath in range(len(skp["path"])):
print("%5d %8s %10.5f %10.5f %10.5f %8s %10.5f %10.5f %10.5f" % (
nkpath[ipath],
skp['path'][ipath][0],
skp['point_coords'][skp['path'][ipath][0]][0],
skp['point_coords'][skp['path'][ipath][0]][1],
skp['point_coords'][skp['path'][ipath][0]][2],
skp['path'][ipath][1],
skp['point_coords'][skp['path'][ipath][1]][0],
skp['point_coords'][skp['path'][ipath][1]][1],
skp['point_coords'][skp['path'][ipath][1]][2]))
#
# Number of electrons
#
nbnd = 0
for iat in atom:
nbnd += valence_dict[iat]
nbnd = nbnd / 2 + len(atom)*20
if rel:
nbnd *= 2
print("Number of Bands : ", nbnd)
#
# Shell scripts
#
structure2 = pymatgen.core.Structure(skp["primitive_lattice"],
skp["primitive_types"], skp["primitive_positions"])
spg_analysis = SpacegroupAnalyzer(structure2)
coarse = spg_analysis.get_ir_reciprocal_mesh(mesh=(nq[0], nq[1], nq[2]), is_shift=(0, 0, 0))
middle = spg_analysis.get_ir_reciprocal_mesh(mesh=(nq[0]*2, nq[1]*2, nq[2]*2), is_shift=(0, 0, 0))
dense = spg_analysis.get_ir_reciprocal_mesh(mesh=(nq[0]*4, nq[1]*4, nq[2]*4), is_shift=(0, 0, 0))
print("Number of irreducible k : ", len(coarse), len(middle), len(dense))
write_sh(nq[0]*nq[1]*nq[2], len(middle), len(dense),
len(kpath), atom, atomwfc_dict, host, numpw*nbnd, rel)
#
# rx.in, scf.in, nscf.in, band.in , nscf_w.in, nscf_r.in
#
write_pwx(skp, ecutwfc, ecutrho, pseudo_dict, nq, nbnd, rel, kpath)
#
# ph.in, elph.in, epmat.in, phdos.in, rpa.in, scdft.in
#
write_ph(nq, ecutrho, nbnd)
#
# bands.in, pp.in, proj.in, pw2wan.in, q2r.in
#
write_pp()
#
# band.gp, pwscf.win, respack.in, disp.in
#
write_wannier(skp, nbnd, nq, atomwfc_dict, kpath)
#
# openmx.in : Input file for openmx
#
write_openmx(skp, nq, rel, nkpath)
#
# HiLAPW input
#
# write_hilapw(skp, nq)