-
Notifications
You must be signed in to change notification settings - Fork 835
/
shengbte.py
270 lines (224 loc) · 8.83 KB
/
shengbte.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
"""This module implements reading and writing of ShengBTE CONTROL files."""
from __future__ import annotations
import warnings
from typing import Any
import numpy as np
from monty.dev import requires
from monty.json import MSONable
from pymatgen.core.structure import Structure
from pymatgen.io.vasp import Kpoints
try:
import f90nml
except ImportError:
f90nml = None
__author__ = "Rees Chang, Alex Ganose"
__copyright__ = "Copyright 2019, The Materials Project"
__version__ = "0.1"
__email__ = "rc564@cornell.edu, aganose@lbl.gov"
__date__ = "June 27, 2019"
class Control(MSONable, dict):
"""
Class for reading, updating, and writing ShengBTE CONTROL files.
See https://bitbucket.org/sousaw/shengbte/src/master/ for more
detailed description and default values of CONTROL arguments.
"""
required_params = data_keys = (
"nelements",
"natoms",
"ngrid",
"lattvec",
"types",
"elements",
"positions",
"scell",
)
allocations_keys = ("nelements", "natoms", "ngrid", "norientations")
crystal_keys = (
"lfactor",
"lattvec",
"types",
"elements",
"positions",
"masses",
"gfactors",
"epsilon",
"born",
"scell",
"orientations",
)
params_keys = (
"t",
"t_min",
"t_max",
"t_step",
"omega_max",
"scalebroad",
"rmin",
"rmax",
"dr",
"maxiter",
"nticks",
"eps",
)
flags_keys = (
"nonanalytic",
"convergence",
"isotopes",
"autoisotopes",
"nanowires",
"onlyharmonic",
"espresso",
)
def __init__(self, ngrid: list[int] | None = None, temperature: float | dict[str, float] = 300, **kwargs):
"""
Args:
ngrid: Reciprocal space grid density as a list of 3 ints.
temperature: The temperature to calculate the lattice thermal
conductivity for. Can be given as a single float, or a dictionary
with the keys "min", "max", "step".
**kwargs: Other ShengBTE parameters. Several parameters are required
for ShengBTE to run - we have listed these parameters below:
- nelements (int): number of different elements in the compound
- natoms (int): number of atoms in the unit cell
- lattvec (size 3x3 array): real-space lattice vectors, in units of lfactor
- lfactor (float): unit of measurement for lattice vectors (nm).
I.e., set to 0.1 if lattvec given in Angstrom.
- types (size natom list): a vector of natom integers, ranging
from 1 to nelements, assigning an element to each atom in the system
- elements (size natom list): a vector of element names
- positions (size natomx3 array): atomic positions in lattice
coordinates
- scell (size 3 list): supercell sizes along each crystal axis
used for the 2nd-order force constant calculation
"""
super().__init__()
if ngrid is None:
ngrid = [25, 25, 25]
self["ngrid"] = ngrid
if isinstance(temperature, (int, float)):
self["t"] = temperature
elif isinstance(temperature, dict):
self["t_min"] = temperature["min"]
self["t_max"] = temperature["max"]
self["t_step"] = temperature["step"]
else:
raise ValueError("Unsupported temperature type, must be float or dict")
self.update(kwargs)
@classmethod
@requires(
f90nml,
"ShengBTE Control object requires f90nml to be installed. Please get it at https://pypi.org/project/f90nml.",
)
def from_file(cls, filepath: str):
"""
Read a CONTROL namelist file and output a 'Control' object.
Args:
filepath: Path of the CONTROL file.
Returns:
'Control' object with parameters instantiated.
"""
nml = f90nml.read(filepath)
sdict = nml.todict()
all_dict: dict[str, Any] = {}
all_dict.update(sdict["allocations"])
all_dict.update(sdict["crystal"])
all_dict.update(sdict["parameters"])
all_dict.update(sdict["flags"])
all_dict.pop("_start_index") # remove unnecessary cruft
return cls.from_dict(all_dict)
@classmethod
def from_dict(cls, control_dict: dict):
"""
Write a CONTROL file from a Python dictionary. Description and default
parameters can be found at
https://bitbucket.org/sousaw/shengbte/src/master/.
Note some parameters are mandatory. Optional parameters default here to
None and will not be written to file.
Args:
control_dict: A Python dictionary of ShengBTE input parameters.
"""
return cls(**control_dict)
@requires(
f90nml,
"ShengBTE Control object requires f90nml to be installed. Please get it at https://pypi.org/project/f90nml.",
)
def to_file(self, filename: str = "CONTROL"):
"""
Writes ShengBTE CONTROL file from 'Control' object.
Args:
filename: A file name.
"""
for param in self.required_params:
if param not in self.as_dict():
warnings.warn(f"Required parameter {param!r} not specified!")
alloc_dict = _get_subdict(self, self.allocations_keys)
alloc_nml = f90nml.Namelist({"allocations": alloc_dict})
control_str = str(alloc_nml) + "\n"
crystal_dict = _get_subdict(self, self.crystal_keys)
crystal_nml = f90nml.Namelist({"crystal": crystal_dict})
control_str += str(crystal_nml) + "\n"
params_dict = _get_subdict(self, self.params_keys)
params_nml = f90nml.Namelist({"parameters": params_dict})
control_str += str(params_nml) + "\n"
flags_dict = _get_subdict(self, self.flags_keys)
flags_nml = f90nml.Namelist({"flags": flags_dict})
control_str += str(flags_nml) + "\n"
with open(filename, mode="w") as file:
file.write(control_str)
@classmethod
def from_structure(cls, structure: Structure, reciprocal_density: int | None = 50000, **kwargs):
"""
Get a ShengBTE control object from a structure.
Args:
structure: A structure object.
reciprocal_density: If not None, the q-point grid ("ngrid") will be
set using this density.
kwargs: Additional options to be passed to the Control constructor.
See the docstring of the __init__ method for more details
Returns:
A ShengBTE control object.
"""
elements = list(map(str, structure.elements))
unique_nums = np.unique(structure.atomic_numbers)
types_dict = dict(zip(unique_nums, range(len(unique_nums))))
types = [types_dict[i] + 1 for i in structure.atomic_numbers]
control_dict = {
"nelements": structure.n_elems,
"natoms": len(structure),
"norientations": 0,
"lfactor": 0.1,
"lattvec": structure.lattice.matrix.tolist(),
"elements": elements,
"types": types,
"positions": structure.frac_coords.tolist(),
}
if reciprocal_density:
kpoints = Kpoints.automatic_density(structure, reciprocal_density)
control_dict["ngrid"] = kpoints.kpts[0]
control_dict.update(**kwargs)
return Control(**control_dict)
def get_structure(self) -> Structure:
"""
Get a pymatgen Structure from a ShengBTE control object.
The control object must have the "lattvec", "types", "elements", and
"positions" settings otherwise an error will be thrown.
Returns:
The structure.
"""
required = ["lattvec", "types", "elements", "positions"]
if not all(r in self for r in required):
raise ValueError("All of ['lattvec', 'types', 'elements', 'positions'] must be in control object")
unique_elements = self["elements"]
n_unique_elements = len(unique_elements)
element_map = dict(zip(range(1, n_unique_elements + 1), unique_elements))
species = [element_map[i] for i in self["types"]]
cell = np.array(self["lattvec"])
if "lfactor" in self:
cell *= self["lfactor"] * 10 # to nm then to Angstrom
return Structure(cell, species, self["positions"])
def as_dict(self):
"""Returns: MSONable dict."""
return dict(self)
def _get_subdict(master_dict, subkeys):
"""Helper method to get a set of keys from a larger dictionary."""
return {k: master_dict[k] for k in subkeys if k in master_dict and master_dict[k] is not None}