Coverage for /builds/debichem-team/python-ase/ase/io/espresso.py : 65.88%
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"""Reads Quantum ESPRESSO files.
3Read multiple structures and results from pw.x output files. Read
4structures from pw.x input files.
6Built for PWSCF v.5.3.0 but should work with earlier and later versions.
7Can deal with most major functionality, but might fail with ibrav =/= 0
8or crystal_sg positions.
10Units are converted using CODATA 2006, as used internally by Quantum
11ESPRESSO.
12"""
14import os
15import operator as op
16import re
17import warnings
18from collections import OrderedDict
19from os import path
21import numpy as np
23from ase.atoms import Atoms
24from ase.calculators.singlepoint import (SinglePointDFTCalculator,
25 SinglePointKPoint)
26from ase.calculators.calculator import kpts2ndarray, kpts2sizeandoffsets
27from ase.dft.kpoints import kpoint_convert
28from ase.constraints import FixAtoms, FixCartesian
29from ase.data import chemical_symbols, atomic_numbers
30from ase.units import create_units
31from ase.utils import iofunction
34# Quantum ESPRESSO uses CODATA 2006 internally
35units = create_units('2006')
37# Section identifiers
38_PW_START = 'Program PWSCF'
39_PW_END = 'End of self-consistent calculation'
40_PW_CELL = 'CELL_PARAMETERS'
41_PW_POS = 'ATOMIC_POSITIONS'
42_PW_MAGMOM = 'Magnetic moment per site'
43_PW_FORCE = 'Forces acting on atoms'
44_PW_TOTEN = '! total energy'
45_PW_STRESS = 'total stress'
46_PW_FERMI = 'the Fermi energy is'
47_PW_HIGHEST_OCCUPIED = 'highest occupied level'
48_PW_HIGHEST_OCCUPIED_LOWEST_FREE = 'highest occupied, lowest unoccupied level'
49_PW_KPTS = 'number of k points='
50_PW_BANDS = _PW_END
51_PW_BANDSTRUCTURE = 'End of band structure calculation'
54class Namelist(OrderedDict):
55 """Case insensitive dict that emulates Fortran Namelists."""
56 def __contains__(self, key):
57 return super(Namelist, self).__contains__(key.lower())
59 def __delitem__(self, key):
60 return super(Namelist, self).__delitem__(key.lower())
62 def __getitem__(self, key):
63 return super(Namelist, self).__getitem__(key.lower())
65 def __setitem__(self, key, value):
66 super(Namelist, self).__setitem__(key.lower(), value)
68 def get(self, key, default=None):
69 return super(Namelist, self).get(key.lower(), default)
72@iofunction('rU')
73def read_espresso_out(fileobj, index=-1, results_required=True):
74 """Reads Quantum ESPRESSO output files.
76 The atomistic configurations as well as results (energy, force, stress,
77 magnetic moments) of the calculation are read for all configurations
78 within the output file.
80 Will probably raise errors for broken or incomplete files.
82 Parameters
83 ----------
84 fileobj : file|str
85 A file like object or filename
86 index : slice
87 The index of configurations to extract.
88 results_required : bool
89 If True, atomistic configurations that do not have any
90 associated results will not be included. This prevents double
91 printed configurations and incomplete calculations from being
92 returned as the final configuration with no results data.
94 Yields
95 ------
96 structure : Atoms
97 The next structure from the index slice. The Atoms has a
98 SinglePointCalculator attached with any results parsed from
99 the file.
102 """
103 # work with a copy in memory for faster random access
104 pwo_lines = fileobj.readlines()
106 # TODO: index -1 special case?
107 # Index all the interesting points
108 indexes = {
109 _PW_START: [],
110 _PW_END: [],
111 _PW_CELL: [],
112 _PW_POS: [],
113 _PW_MAGMOM: [],
114 _PW_FORCE: [],
115 _PW_TOTEN: [],
116 _PW_STRESS: [],
117 _PW_FERMI: [],
118 _PW_HIGHEST_OCCUPIED: [],
119 _PW_HIGHEST_OCCUPIED_LOWEST_FREE: [],
120 _PW_KPTS: [],
121 _PW_BANDS: [],
122 _PW_BANDSTRUCTURE: [],
123 }
125 for idx, line in enumerate(pwo_lines):
126 for identifier in indexes:
127 if identifier in line:
128 indexes[identifier].append(idx)
130 # Configurations are either at the start, or defined in ATOMIC_POSITIONS
131 # in a subsequent step. Can deal with concatenated output files.
132 all_config_indexes = sorted(indexes[_PW_START] +
133 indexes[_PW_POS])
135 # Slice only requested indexes
136 # setting results_required argument stops configuration-only
137 # structures from being returned. This ensures the [-1] structure
138 # is one that has results. Two cases:
139 # - SCF of last configuration is not converged, job terminated
140 # abnormally.
141 # - 'relax' and 'vc-relax' re-prints the final configuration but
142 # only 'vc-relax' recalculates.
143 if results_required:
144 results_indexes = sorted(indexes[_PW_TOTEN] + indexes[_PW_FORCE] +
145 indexes[_PW_STRESS] + indexes[_PW_MAGMOM] +
146 indexes[_PW_BANDS] +
147 indexes[_PW_BANDSTRUCTURE])
149 # Prune to only configurations with results data before the next
150 # configuration
151 results_config_indexes = []
152 for config_index, config_index_next in zip(
153 all_config_indexes,
154 all_config_indexes[1:] + [len(pwo_lines)]):
155 if any([config_index < results_index < config_index_next
156 for results_index in results_indexes]):
157 results_config_indexes.append(config_index)
159 # slice from the subset
160 image_indexes = results_config_indexes[index]
161 else:
162 image_indexes = all_config_indexes[index]
164 # Extract initialisation information each time PWSCF starts
165 # to add to subsequent configurations. Use None so slices know
166 # when to fill in the blanks.
167 pwscf_start_info = dict((idx, None) for idx in indexes[_PW_START])
169 for image_index in image_indexes:
170 # Find the nearest calculation start to parse info. Needed in,
171 # for example, relaxation where cell is only printed at the
172 # start.
173 if image_index in indexes[_PW_START]:
174 prev_start_index = image_index
175 else:
176 # The greatest start index before this structure
177 prev_start_index = [idx for idx in indexes[_PW_START]
178 if idx < image_index][-1]
180 # add structure to reference if not there
181 if pwscf_start_info[prev_start_index] is None:
182 pwscf_start_info[prev_start_index] = parse_pwo_start(
183 pwo_lines, prev_start_index)
185 # Get the bounds for information for this structure. Any associated
186 # values will be between the image_index and the following one,
187 # EXCEPT for cell, which will be 4 lines before if it exists.
188 for next_index in all_config_indexes:
189 if next_index > image_index:
190 break
191 else:
192 # right to the end of the file
193 next_index = len(pwo_lines)
195 # Get the structure
196 # Use this for any missing data
197 prev_structure = pwscf_start_info[prev_start_index]['atoms']
198 if image_index in indexes[_PW_START]:
199 structure = prev_structure.copy() # parsed from start info
200 else:
201 if _PW_CELL in pwo_lines[image_index - 5]:
202 # CELL_PARAMETERS would be just before positions if present
203 cell, cell_alat = get_cell_parameters(
204 pwo_lines[image_index - 5:image_index])
205 else:
206 cell = prev_structure.cell
207 cell_alat = pwscf_start_info[prev_start_index]['alat']
209 # give at least enough lines to parse the positions
210 # should be same format as input card
211 n_atoms = len(prev_structure)
212 positions_card = get_atomic_positions(
213 pwo_lines[image_index:image_index + n_atoms + 1],
214 n_atoms=n_atoms, cell=cell, alat=cell_alat)
216 # convert to Atoms object
217 symbols = [label_to_symbol(position[0]) for position in
218 positions_card]
219 positions = [position[1] for position in positions_card]
220 structure = Atoms(symbols=symbols, positions=positions, cell=cell,
221 pbc=True)
223 # Extract calculation results
224 # Energy
225 energy = None
226 for energy_index in indexes[_PW_TOTEN]:
227 if image_index < energy_index < next_index:
228 energy = float(
229 pwo_lines[energy_index].split()[-2]) * units['Ry']
231 # Forces
232 forces = None
233 for force_index in indexes[_PW_FORCE]:
234 if image_index < force_index < next_index:
235 # Before QE 5.3 'negative rho' added 2 lines before forces
236 # Use exact lines to stop before 'non-local' forces
237 # in high verbosity
238 if not pwo_lines[force_index + 2].strip():
239 force_index += 4
240 else:
241 force_index += 2
242 # assume contiguous
243 forces = [
244 [float(x) for x in force_line.split()[-3:]] for force_line
245 in pwo_lines[force_index:force_index + len(structure)]]
246 forces = np.array(forces) * units['Ry'] / units['Bohr']
248 # Stress
249 stress = None
250 for stress_index in indexes[_PW_STRESS]:
251 if image_index < stress_index < next_index:
252 sxx, sxy, sxz = pwo_lines[stress_index + 1].split()[:3]
253 _, syy, syz = pwo_lines[stress_index + 2].split()[:3]
254 _, _, szz = pwo_lines[stress_index + 3].split()[:3]
255 stress = np.array([sxx, syy, szz, syz, sxz, sxy], dtype=float)
256 # sign convention is opposite of ase
257 stress *= -1 * units['Ry'] / (units['Bohr'] ** 3)
259 # Magmoms
260 magmoms = None
261 for magmoms_index in indexes[_PW_MAGMOM]:
262 if image_index < magmoms_index < next_index:
263 magmoms = [
264 float(mag_line.split()[5]) for mag_line
265 in pwo_lines[magmoms_index + 1:
266 magmoms_index + 1 + len(structure)]]
268 # Fermi level / highest occupied level
269 efermi = None
270 for fermi_index in indexes[_PW_FERMI]:
271 if image_index < fermi_index < next_index:
272 efermi = float(pwo_lines[fermi_index].split()[-2])
274 if efermi is None:
275 for ho_index in indexes[_PW_HIGHEST_OCCUPIED]:
276 if image_index < ho_index < next_index:
277 efermi = float(pwo_lines[ho_index].split()[-1])
279 if efermi is None:
280 for holf_index in indexes[_PW_HIGHEST_OCCUPIED_LOWEST_FREE]:
281 if image_index < holf_index < next_index:
282 efermi = float(pwo_lines[holf_index].split()[-2])
284 # K-points
285 ibzkpts = None
286 weights = None
287 kpoints_warning = "Number of k-points >= 100: " + \
288 "set verbosity='high' to print them."
290 for kpts_index in indexes[_PW_KPTS]:
291 nkpts = int(pwo_lines[kpts_index].split()[4])
292 kpts_index += 2
294 if pwo_lines[kpts_index].strip() == kpoints_warning:
295 continue
297 # QE prints the k-points in units of 2*pi/alat
298 # with alat defined as the length of the first
299 # cell vector
300 cell = structure.get_cell()
301 alat = np.linalg.norm(cell[0])
302 ibzkpts = []
303 weights = []
304 for i in range(nkpts):
305 L = pwo_lines[kpts_index + i].split()
306 weights.append(float(L[-1]))
307 coord = np.array([L[-6], L[-5], L[-4].strip('),')],
308 dtype=float)
309 coord *= 2 * np.pi / alat
310 coord = kpoint_convert(cell, ckpts_kv=coord)
311 ibzkpts.append(coord)
312 ibzkpts = np.array(ibzkpts)
313 weights = np.array(weights)
315 # Bands
316 kpts = None
317 kpoints_warning = "Number of k-points >= 100: " + \
318 "set verbosity='high' to print the bands."
320 for bands_index in indexes[_PW_BANDS] + indexes[_PW_BANDSTRUCTURE]:
321 if image_index < bands_index < next_index:
322 bands_index += 2
324 if pwo_lines[bands_index].strip() == kpoints_warning:
325 continue
327 assert ibzkpts is not None
328 spin, bands, eigenvalues = 0, [], [[], []]
330 while True:
331 L = pwo_lines[bands_index].replace('-', ' -').split()
332 if len(L) == 0:
333 if len(bands) > 0:
334 eigenvalues[spin].append(bands)
335 bands = []
336 elif L == ['occupation', 'numbers']:
337 # Skip the lines with the occupation numbers
338 bands_index += len(eigenvalues[spin][0]) // 8 + 1
339 elif L[0] == 'k' and L[1].startswith('='):
340 pass
341 elif 'SPIN' in L:
342 if 'DOWN' in L:
343 spin += 1
344 else:
345 try:
346 bands.extend(map(float, L))
347 except ValueError:
348 break
349 bands_index += 1
351 if spin == 1:
352 assert len(eigenvalues[0]) == len(eigenvalues[1])
353 assert len(eigenvalues[0]) == len(ibzkpts), \
354 (np.shape(eigenvalues), len(ibzkpts))
356 kpts = []
357 for s in range(spin + 1):
358 for w, k, e in zip(weights, ibzkpts, eigenvalues[s]):
359 kpt = SinglePointKPoint(w, s, k, eps_n=e)
360 kpts.append(kpt)
362 # Put everything together
363 #
364 # I have added free_energy. Can and should we distinguish
365 # energy and free_energy? --askhl
366 calc = SinglePointDFTCalculator(structure, energy=energy,
367 free_energy=energy,
368 forces=forces, stress=stress,
369 magmoms=magmoms, efermi=efermi,
370 ibzkpts=ibzkpts)
371 calc.kpts = kpts
372 structure.calc = calc
374 yield structure
377def parse_pwo_start(lines, index=0):
378 """Parse Quantum ESPRESSO calculation info from lines,
379 starting from index. Return a dictionary containing extracted
380 information.
382 - `celldm(1)`: lattice parameters (alat)
383 - `cell`: unit cell in Angstrom
384 - `symbols`: element symbols for the structure
385 - `positions`: cartesian coordinates of atoms in Angstrom
386 - `atoms`: an `ase.Atoms` object constructed from the extracted data
388 Parameters
389 ----------
390 lines : list[str]
391 Contents of PWSCF output file.
392 index : int
393 Line number to begin parsing. Only first calculation will
394 be read.
396 Returns
397 -------
398 info : dict
399 Dictionary of calculation parameters, including `celldm(1)`, `cell`,
400 `symbols`, `positions`, `atoms`.
402 Raises
403 ------
404 KeyError
405 If interdependent values cannot be found (especially celldm(1))
406 an error will be raised as other quantities cannot then be
407 calculated (e.g. cell and positions).
408 """
409 # TODO: extend with extra DFT info?
411 info = {}
413 for idx, line in enumerate(lines[index:], start=index):
414 if 'celldm(1)' in line:
415 # celldm(1) has more digits than alat!!
416 info['celldm(1)'] = float(line.split()[1]) * units['Bohr']
417 info['alat'] = info['celldm(1)']
418 elif 'number of atoms/cell' in line:
419 info['nat'] = int(line.split()[-1])
420 elif 'number of atomic types' in line:
421 info['ntyp'] = int(line.split()[-1])
422 elif 'crystal axes:' in line:
423 info['cell'] = info['celldm(1)'] * np.array([
424 [float(x) for x in lines[idx + 1].split()[3:6]],
425 [float(x) for x in lines[idx + 2].split()[3:6]],
426 [float(x) for x in lines[idx + 3].split()[3:6]]])
427 elif 'positions (alat units)' in line:
428 info['symbols'], info['positions'] = [], []
430 for at_line in lines[idx + 1:idx + 1 + info['nat']]:
431 sym, x, y, z = parse_position_line(at_line)
432 info['symbols'].append(label_to_symbol(sym))
433 info['positions'].append([x * info['celldm(1)'],
434 y * info['celldm(1)'],
435 z * info['celldm(1)']])
436 # This should be the end of interesting info.
437 # Break here to avoid dealing with large lists of kpoints.
438 # Will need to be extended for DFTCalculator info.
439 break
441 # Make atoms for convenience
442 info['atoms'] = Atoms(symbols=info['symbols'],
443 positions=info['positions'],
444 cell=info['cell'], pbc=True)
446 return info
449def parse_position_line(line):
450 """Parse a single line from a pw.x output file.
452 The line must contain information about the atomic symbol and the position,
453 e.g.
455 995 Sb tau( 995) = ( 1.4212023 0.7037863 0.1242640 )
457 Parameters
458 ----------
459 line : str
460 Line to be parsed.
462 Returns
463 -------
464 sym : str
465 Atomic symbol.
466 x : float
467 x-position.
468 y : float
469 y-position.
470 z : float
471 z-position.
472 """
473 pat = re.compile(r'\s*\d+\s*(\S+)\s*tau\(\s*\d+\)\s*='
474 r'\s*\(\s*(\S+)\s+(\S+)\s+(\S+)\s*\)')
475 match = pat.match(line)
476 assert match is not None
477 sym, x, y, z = match.group(1, 2, 3, 4)
478 return sym, float(x), float(y), float(z)
481@iofunction('rU')
482def read_espresso_in(fileobj):
483 """Parse a Quantum ESPRESSO input files, '.in', '.pwi'.
485 ESPRESSO inputs are generally a fortran-namelist format with custom
486 blocks of data. The namelist is parsed as a dict and an atoms object
487 is constructed from the included information.
489 Parameters
490 ----------
491 fileobj : file | str
492 A file-like object that supports line iteration with the contents
493 of the input file, or a filename.
495 Returns
496 -------
497 atoms : Atoms
498 Structure defined in the input file.
500 Raises
501 ------
502 KeyError
503 Raised for missing keys that are required to process the file
504 """
505 # parse namelist section and extract remaining lines
506 data, card_lines = read_fortran_namelist(fileobj)
508 # get the cell if ibrav=0
509 if 'system' not in data:
510 raise KeyError('Required section &SYSTEM not found.')
511 elif 'ibrav' not in data['system']:
512 raise KeyError('ibrav is required in &SYSTEM')
513 elif data['system']['ibrav'] == 0:
514 # celldm(1) is in Bohr, A is in angstrom. celldm(1) will be
515 # used even if A is also specified.
516 if 'celldm(1)' in data['system']:
517 alat = data['system']['celldm(1)'] * units['Bohr']
518 elif 'A' in data['system']:
519 alat = data['system']['A']
520 else:
521 alat = None
522 cell, cell_alat = get_cell_parameters(card_lines, alat=alat)
523 else:
524 alat, cell = ibrav_to_cell(data['system'])
526 # species_info holds some info for each element
527 species_card = get_atomic_species(card_lines, n_species=data['system']['ntyp'])
528 species_info = {}
529 for ispec, (label, weight, pseudo) in enumerate(species_card):
530 symbol = label_to_symbol(label)
531 valence = get_valence_electrons(symbol, data, pseudo)
533 # starting_magnetization is in fractions of valence electrons
534 magnet_key = "starting_magnetization({0})".format(ispec + 1)
535 magmom = valence * data["system"].get(magnet_key, 0.0)
536 species_info[symbol] = {"weight": weight, "pseudo": pseudo,
537 "valence": valence, "magmom": magmom}
539 positions_card = get_atomic_positions(
540 card_lines, n_atoms=data['system']['nat'], cell=cell, alat=alat)
542 symbols = [label_to_symbol(position[0]) for position in positions_card]
543 positions = [position[1] for position in positions_card]
544 magmoms = [species_info[symbol]["magmom"] for symbol in symbols]
546 # TODO: put more info into the atoms object
547 # e.g magmom, forces.
548 atoms = Atoms(symbols=symbols, positions=positions, cell=cell, pbc=True,
549 magmoms=magmoms)
551 return atoms
554def ibrav_to_cell(system):
555 """
556 Convert a value of ibrav to a cell. Any unspecified lattice dimension
557 is set to 0.0, but will not necessarily raise an error. Also return the
558 lattice parameter.
560 Parameters
561 ----------
562 system : dict
563 The &SYSTEM section of the input file, containing the 'ibrav' setting,
564 and either celldm(1)..(6) or a, b, c, cosAB, cosAC, cosBC.
566 Returns
567 -------
568 alat, cell : float, np.array
569 Cell parameter in Angstrom, and
570 The 3x3 array representation of the cell.
572 Raises
573 ------
574 KeyError
575 Raise an error if any required keys are missing.
576 NotImplementedError
577 Only a limited number of ibrav settings can be parsed. An error
578 is raised if the ibrav interpretation is not implemented.
579 """
580 if 'celldm(1)' in system and 'a' in system:
581 raise KeyError('do not specify both celldm and a,b,c!')
582 elif 'celldm(1)' in system:
583 # celldm(x) in bohr
584 alat = system['celldm(1)'] * units['Bohr']
585 b_over_a = system.get('celldm(2)', 0.0)
586 c_over_a = system.get('celldm(3)', 0.0)
587 cosab = system.get('celldm(4)', 0.0)
588 cosac = system.get('celldm(5)', 0.0)
589 cosbc = 0.0
590 if system['ibrav'] == 14:
591 cosbc = system.get('celldm(4)', 0.0)
592 cosac = system.get('celldm(5)', 0.0)
593 cosab = system.get('celldm(6)', 0.0)
594 elif 'a' in system:
595 # a, b, c, cosAB, cosAC, cosBC in Angstrom
596 alat = system['a']
597 b_over_a = system.get('b', 0.0) / alat
598 c_over_a = system.get('c', 0.0) / alat
599 cosab = system.get('cosab', 0.0)
600 cosac = system.get('cosac', 0.0)
601 cosbc = system.get('cosbc', 0.0)
602 else:
603 raise KeyError("Missing celldm(1) or a cell parameter.")
605 if system['ibrav'] == 1:
606 cell = np.identity(3) * alat
607 elif system['ibrav'] == 2:
608 cell = np.array([[-1.0, 0.0, 1.0],
609 [0.0, 1.0, 1.0],
610 [-1.0, 1.0, 0.0]]) * (alat / 2)
611 elif system['ibrav'] == 3:
612 cell = np.array([[1.0, 1.0, 1.0],
613 [-1.0, 1.0, 1.0],
614 [-1.0, -1.0, 1.0]]) * (alat / 2)
615 elif system['ibrav'] == -3:
616 cell = np.array([[-1.0, 1.0, 1.0],
617 [1.0, -1.0, 1.0],
618 [1.0, 1.0, -1.0]]) * (alat / 2)
619 elif system['ibrav'] == 4:
620 cell = np.array([[1.0, 0.0, 0.0],
621 [-0.5, 0.5 * 3**0.5, 0.0],
622 [0.0, 0.0, c_over_a]]) * alat
623 elif system['ibrav'] == 5:
624 tx = ((1.0 - cosab) / 2.0)**0.5
625 ty = ((1.0 - cosab) / 6.0)**0.5
626 tz = ((1 + 2 * cosab) / 3.0)**0.5
627 cell = np.array([[tx, -ty, tz],
628 [0, 2 * ty, tz],
629 [-tx, -ty, tz]]) * alat
630 elif system['ibrav'] == -5:
631 ty = ((1.0 - cosab) / 6.0)**0.5
632 tz = ((1 + 2 * cosab) / 3.0)**0.5
633 a_prime = alat / 3**0.5
634 u = tz - 2 * 2**0.5 * ty
635 v = tz + 2**0.5 * ty
636 cell = np.array([[u, v, v],
637 [v, u, v],
638 [v, v, u]]) * a_prime
639 elif system['ibrav'] == 6:
640 cell = np.array([[1.0, 0.0, 0.0],
641 [0.0, 1.0, 0.0],
642 [0.0, 0.0, c_over_a]]) * alat
643 elif system['ibrav'] == 7:
644 cell = np.array([[1.0, -1.0, c_over_a],
645 [1.0, 1.0, c_over_a],
646 [-1.0, -1.0, c_over_a]]) * (alat / 2)
647 elif system['ibrav'] == 8:
648 cell = np.array([[1.0, 0.0, 0.0],
649 [0.0, b_over_a, 0.0],
650 [0.0, 0.0, c_over_a]]) * alat
651 elif system['ibrav'] == 9:
652 cell = np.array([[1.0 / 2.0, b_over_a / 2.0, 0.0],
653 [-1.0 / 2.0, b_over_a / 2.0, 0.0],
654 [0.0, 0.0, c_over_a]]) * alat
655 elif system['ibrav'] == -9:
656 cell = np.array([[1.0 / 2.0, -b_over_a / 2.0, 0.0],
657 [1.0 / 2.0, b_over_a / 2.0, 0.0],
658 [0.0, 0.0, c_over_a]]) * alat
659 elif system['ibrav'] == 10:
660 cell = np.array([[1.0 / 2.0, 0.0, c_over_a / 2.0],
661 [1.0 / 2.0, b_over_a / 2.0, 0.0],
662 [0.0, b_over_a / 2.0, c_over_a / 2.0]]) * alat
663 elif system['ibrav'] == 11:
664 cell = np.array([[1.0 / 2.0, b_over_a / 2.0, c_over_a / 2.0],
665 [-1.0 / 2.0, b_over_a / 2.0, c_over_a / 2.0],
666 [-1.0, 2.0, -b_over_a / 2.0, c_over_a / 2.0]]) * alat
667 elif system['ibrav'] == 12:
668 sinab = (1.0 - cosab**2)**0.5
669 cell = np.array([[1.0, 0.0, 0.0],
670 [b_over_a * cosab, b_over_a * sinab, 0.0],
671 [0.0, 0.0, c_over_a]]) * alat
672 elif system['ibrav'] == -12:
673 sinac = (1.0 - cosac**2)**0.5
674 cell = np.array([[1.0, 0.0, 0.0],
675 [0.0, b_over_a, 0.0],
676 [c_over_a * cosac, 0.0, c_over_a * sinac]]) * alat
677 elif system['ibrav'] == 13:
678 sinab = (1.0 - cosab**2)**0.5
679 cell = np.array([[1.0 / 2.0, 0.0, -c_over_a / 2.0],
680 [b_over_a * cosab, b_over_a * sinab, 0.0],
681 [1.0 / 2.0, 0.0, c_over_a / 2.0]]) * alat
682 elif system['ibrav'] == 14:
683 sinab = (1.0 - cosab**2)**0.5
684 v3 = [c_over_a * cosac,
685 c_over_a * (cosbc - cosac * cosab) / sinab,
686 c_over_a * ((1 + 2 * cosbc * cosac * cosab
687 - cosbc**2 - cosac**2 - cosab**2)**0.5) / sinab]
688 cell = np.array([[1.0, 0.0, 0.0],
689 [b_over_a * cosab, b_over_a * sinab, 0.0],
690 v3]) * alat
691 else:
692 raise NotImplementedError('ibrav = {0} is not implemented'
693 ''.format(system['ibrav']))
695 return alat, cell
698def get_pseudo_dirs(data):
699 """Guess a list of possible locations for pseudopotential files.
701 Parameters
702 ----------
703 data : Namelist
704 Namelist representing the quantum espresso input parameters
706 Returns
707 -------
708 pseudo_dirs : list[str]
709 A list of directories where pseudopotential files could be located.
710 """
711 pseudo_dirs = []
712 if 'pseudo_dir' in data['control']:
713 pseudo_dirs.append(data['control']['pseudo_dir'])
714 if 'ESPRESSO_PSEUDO' in os.environ:
715 pseudo_dirs.append(os.environ['ESPRESSO_PSEUDO'])
716 pseudo_dirs.append(path.expanduser('~/espresso/pseudo/'))
717 return pseudo_dirs
720def get_valence_electrons(symbol, data, pseudo=None):
721 """The number of valence electrons for a atomic symbol.
723 Parameters
724 ----------
725 symbol : str
726 Chemical symbol
728 data : Namelist
729 Namelist representing the quantum espresso input parameters
731 pseudo : str, optional
732 File defining the pseudopotential to be used. If missing a fallback
733 to the number of valence electrons recommended at
734 http://materialscloud.org/sssp/ is employed.
735 """
736 if pseudo is None:
737 pseudo = '{}_dummy.UPF'.format(symbol)
738 for pseudo_dir in get_pseudo_dirs(data):
739 if path.exists(path.join(pseudo_dir, pseudo)):
740 valence = grep_valence(path.join(pseudo_dir, pseudo))
741 break
742 else: # not found in a file
743 valence = SSSP_VALENCE[atomic_numbers[symbol]]
744 return valence
747def get_atomic_positions(lines, n_atoms, cell=None, alat=None):
748 """Parse atom positions from ATOMIC_POSITIONS card.
750 Parameters
751 ----------
752 lines : list[str]
753 A list of lines containing the ATOMIC_POSITIONS card.
754 n_atoms : int
755 Expected number of atoms. Only this many lines will be parsed.
756 cell : np.array
757 Unit cell of the crystal. Only used with crystal coordinates.
758 alat : float
759 Lattice parameter for atomic coordinates. Only used for alat case.
761 Returns
762 -------
763 positions : list[(str, (float, float, float), (float, float, float))]
764 A list of the ordered atomic positions in the format:
765 label, (x, y, z), (if_x, if_y, if_z)
766 Force multipliers are set to None if not present.
768 Raises
769 ------
770 ValueError
771 Any problems parsing the data result in ValueError
773 """
775 positions = None
776 # no blanks or comment lines, can the consume n_atoms lines for positions
777 trimmed_lines = (line for line in lines
778 if line.strip() and not line[0] == '#')
780 for line in trimmed_lines:
781 if line.strip().startswith('ATOMIC_POSITIONS'):
782 if positions is not None:
783 raise ValueError('Multiple ATOMIC_POSITIONS specified')
784 # Priority and behaviour tested with QE 5.3
785 if 'crystal_sg' in line.lower():
786 raise NotImplementedError('CRYSTAL_SG not implemented')
787 elif 'crystal' in line.lower():
788 cell = cell
789 elif 'bohr' in line.lower():
790 cell = np.identity(3) * units['Bohr']
791 elif 'angstrom' in line.lower():
792 cell = np.identity(3)
793 # elif 'alat' in line.lower():
794 # cell = np.identity(3) * alat
795 else:
796 if alat is None:
797 raise ValueError('Set lattice parameter in &SYSTEM for '
798 'alat coordinates')
799 # Always the default, will be DEPRECATED as mandatory
800 # in future
801 cell = np.identity(3) * alat
803 positions = []
804 for _dummy in range(n_atoms):
805 split_line = next(trimmed_lines).split()
806 # These can be fractions and other expressions
807 position = np.dot((infix_float(split_line[1]),
808 infix_float(split_line[2]),
809 infix_float(split_line[3])), cell)
810 if len(split_line) > 4:
811 force_mult = (float(split_line[4]),
812 float(split_line[5]),
813 float(split_line[6]))
814 else:
815 force_mult = None
817 positions.append((split_line[0], position, force_mult))
819 return positions
822def get_atomic_species(lines, n_species):
823 """Parse atomic species from ATOMIC_SPECIES card.
825 Parameters
826 ----------
827 lines : list[str]
828 A list of lines containing the ATOMIC_POSITIONS card.
829 n_species : int
830 Expected number of atom types. Only this many lines will be parsed.
832 Returns
833 -------
834 species : list[(str, float, str)]
836 Raises
837 ------
838 ValueError
839 Any problems parsing the data result in ValueError
840 """
842 species = None
843 # no blanks or comment lines, can the consume n_atoms lines for positions
844 trimmed_lines = (line.strip() for line in lines
845 if line.strip() and not line.startswith('#'))
847 for line in trimmed_lines:
848 if line.startswith('ATOMIC_SPECIES'):
849 if species is not None:
850 raise ValueError('Multiple ATOMIC_SPECIES specified')
852 species = []
853 for _dummy in range(n_species):
854 label_weight_pseudo = next(trimmed_lines).split()
855 species.append((label_weight_pseudo[0],
856 float(label_weight_pseudo[1]),
857 label_weight_pseudo[2]))
859 return species
862def get_cell_parameters(lines, alat=None):
863 """Parse unit cell from CELL_PARAMETERS card.
865 Parameters
866 ----------
867 lines : list[str]
868 A list with lines containing the CELL_PARAMETERS card.
869 alat : float | None
870 Unit of lattice vectors in Angstrom. Only used if the card is
871 given in units of alat. alat must be None if CELL_PARAMETERS card
872 is in Bohr or Angstrom. For output files, alat will be parsed from
873 the card header and used in preference to this value.
875 Returns
876 -------
877 cell : np.array | None
878 Cell parameters as a 3x3 array in Angstrom. If no cell is found
879 None will be returned instead.
880 cell_alat : float | None
881 If a value for alat is given in the card header, this is also
882 returned, otherwise this will be None.
884 Raises
885 ------
886 ValueError
887 If CELL_PARAMETERS are given in units of bohr or angstrom
888 and alat is not
889 """
891 cell = None
892 cell_alat = None
893 # no blanks or comment lines, can take three lines for cell
894 trimmed_lines = (line for line in lines
895 if line.strip() and not line[0] == '#')
897 for line in trimmed_lines:
898 if line.strip().startswith('CELL_PARAMETERS'):
899 if cell is not None:
900 # multiple definitions
901 raise ValueError('CELL_PARAMETERS specified multiple times')
902 # Priority and behaviour tested with QE 5.3
903 if 'bohr' in line.lower():
904 if alat is not None:
905 raise ValueError('Lattice parameters given in '
906 '&SYSTEM celldm/A and CELL_PARAMETERS '
907 'bohr')
908 cell_units = units['Bohr']
909 elif 'angstrom' in line.lower():
910 if alat is not None:
911 raise ValueError('Lattice parameters given in '
912 '&SYSTEM celldm/A and CELL_PARAMETERS '
913 'angstrom')
914 cell_units = 1.0
915 elif 'alat' in line.lower():
916 # Output file has (alat = value) (in Bohrs)
917 if '=' in line:
918 alat = float(line.strip(') \n').split()[-1]) * units['Bohr']
919 cell_alat = alat
920 elif alat is None:
921 raise ValueError('Lattice parameters must be set in '
922 '&SYSTEM for alat units')
923 cell_units = alat
924 elif alat is None:
925 # may be DEPRECATED in future
926 cell_units = units['Bohr']
927 else:
928 # may be DEPRECATED in future
929 cell_units = alat
930 # Grab the parameters; blank lines have been removed
931 cell = [[ffloat(x) for x in next(trimmed_lines).split()[:3]],
932 [ffloat(x) for x in next(trimmed_lines).split()[:3]],
933 [ffloat(x) for x in next(trimmed_lines).split()[:3]]]
934 cell = np.array(cell) * cell_units
936 return cell, cell_alat
939def str_to_value(string):
940 """Attempt to convert string into int, float (including fortran double),
941 or bool, in that order, otherwise return the string.
942 Valid (case-insensitive) bool values are: '.true.', '.t.', 'true'
943 and 't' (or false equivalents).
945 Parameters
946 ----------
947 string : str
948 Test to parse for a datatype
950 Returns
951 -------
952 value : any
953 Parsed string as the most appropriate datatype of int, float,
954 bool or string.
956 """
958 # Just an integer
959 try:
960 return int(string)
961 except ValueError:
962 pass
963 # Standard float
964 try:
965 return float(string)
966 except ValueError:
967 pass
968 # Fortran double
969 try:
970 return ffloat(string)
971 except ValueError:
972 pass
974 # possible bool, else just the raw string
975 if string.lower() in ('.true.', '.t.', 'true', 't'):
976 return True
977 elif string.lower() in ('.false.', '.f.', 'false', 'f'):
978 return False
979 else:
980 return string.strip("'")
983def read_fortran_namelist(fileobj):
984 """Takes a fortran-namelist formatted file and returns nested
985 dictionaries of sections and key-value data, followed by a list
986 of lines of text that do not fit the specifications.
988 Behaviour is taken from Quantum ESPRESSO 5.3. Parses fairly
989 convoluted files the same way that QE should, but may not get
990 all the MANDATORY rules and edge cases for very non-standard files:
991 Ignores anything after '!' in a namelist, split pairs on ','
992 to include multiple key=values on a line, read values on section
993 start and end lines, section terminating character, '/', can appear
994 anywhere on a line.
995 All of these are ignored if the value is in 'quotes'.
997 Parameters
998 ----------
999 fileobj : file
1000 An open file-like object.
1002 Returns
1003 -------
1004 data : dict of dict
1005 Dictionary for each section in the namelist with key = value
1006 pairs of data.
1007 card_lines : list of str
1008 Any lines not used to create the data, assumed to belong to 'cards'
1009 in the input file.
1011 """
1012 # Espresso requires the correct order
1013 data = Namelist()
1014 card_lines = []
1015 in_namelist = False
1016 section = 'none' # can't be in a section without changing this
1018 for line in fileobj:
1019 # leading and trailing whitespace never needed
1020 line = line.strip()
1021 if line.startswith('&'):
1022 # inside a namelist
1023 section = line.split()[0][1:].lower() # case insensitive
1024 if section in data:
1025 # Repeated sections are completely ignored.
1026 # (Note that repeated keys overwrite within a section)
1027 section = "_ignored"
1028 data[section] = Namelist()
1029 in_namelist = True
1030 if not in_namelist and line:
1031 # Stripped line is Truthy, so safe to index first character
1032 if line[0] not in ('!', '#'):
1033 card_lines.append(line)
1034 if in_namelist:
1035 # parse k, v from line:
1036 key = []
1037 value = None
1038 in_quotes = False
1039 for character in line:
1040 if character == ',' and value is not None and not in_quotes:
1041 # finished value:
1042 data[section][''.join(key).strip()] = str_to_value(
1043 ''.join(value).strip())
1044 key = []
1045 value = None
1046 elif character == '=' and value is None and not in_quotes:
1047 # start writing value
1048 value = []
1049 elif character == "'":
1050 # only found in value anyway
1051 in_quotes = not in_quotes
1052 value.append("'")
1053 elif character == '!' and not in_quotes:
1054 break
1055 elif character == '/' and not in_quotes:
1056 in_namelist = False
1057 break
1058 elif value is not None:
1059 value.append(character)
1060 else:
1061 key.append(character)
1062 if value is not None:
1063 data[section][''.join(key).strip()] = str_to_value(
1064 ''.join(value).strip())
1066 return data, card_lines
1069def ffloat(string):
1070 """Parse float from fortran compatible float definitions.
1072 In fortran exponents can be defined with 'd' or 'q' to symbolise
1073 double or quad precision numbers. Double precision numbers are
1074 converted to python floats and quad precision values are interpreted
1075 as numpy longdouble values (platform specific precision).
1077 Parameters
1078 ----------
1079 string : str
1080 A string containing a number in fortran real format
1082 Returns
1083 -------
1084 value : float | np.longdouble
1085 Parsed value of the string.
1087 Raises
1088 ------
1089 ValueError
1090 Unable to parse a float value.
1092 """
1094 if 'q' in string.lower():
1095 return np.longdouble(string.lower().replace('q', 'e'))
1096 else:
1097 return float(string.lower().replace('d', 'e'))
1100def label_to_symbol(label):
1101 """Convert a valid espresso ATOMIC_SPECIES label to a
1102 chemical symbol.
1104 Parameters
1105 ----------
1106 label : str
1107 chemical symbol X (1 or 2 characters, case-insensitive)
1108 or chemical symbol plus a number or a letter, as in
1109 "Xn" (e.g. Fe1) or "X_*" or "X-*" (e.g. C1, C_h;
1110 max total length cannot exceed 3 characters).
1112 Returns
1113 -------
1114 symbol : str
1115 The best matching species from ase.utils.chemical_symbols
1117 Raises
1118 ------
1119 KeyError
1120 Couldn't find an appropriate species.
1122 Notes
1123 -----
1124 It's impossible to tell whether e.g. He is helium
1125 or hydrogen labelled 'e'.
1126 """
1128 # possibly a two character species
1129 # ase Atoms need proper case of chemical symbols.
1130 if len(label) >= 2:
1131 test_symbol = label[0].upper() + label[1].lower()
1132 if test_symbol in chemical_symbols:
1133 return test_symbol
1134 # finally try with one character
1135 test_symbol = label[0].upper()
1136 if test_symbol in chemical_symbols:
1137 return test_symbol
1138 else:
1139 raise KeyError('Could not parse species from label {0}.'
1140 ''.format(label))
1143def infix_float(text):
1144 """Parse simple infix maths into a float for compatibility with
1145 Quantum ESPRESSO ATOMIC_POSITIONS cards. Note: this works with the
1146 example, and most simple expressions, but the capabilities of
1147 the two parsers are not identical. Will also parse a normal float
1148 value properly, but slowly.
1150 >>> infix_float('1/2*3^(-1/2)')
1151 0.28867513459481287
1153 Parameters
1154 ----------
1155 text : str
1156 An arithmetic expression using +, -, *, / and ^, including brackets.
1158 Returns
1159 -------
1160 value : float
1161 Result of the mathematical expression.
1163 """
1165 def middle_brackets(full_text):
1166 """Extract text from innermost brackets."""
1167 start, end = 0, len(full_text)
1168 for (idx, char) in enumerate(full_text):
1169 if char == '(':
1170 start = idx
1171 if char == ')':
1172 end = idx + 1
1173 break
1174 return full_text[start:end]
1176 def eval_no_bracket_expr(full_text):
1177 """Calculate value of a mathematical expression, no brackets."""
1178 exprs = [('+', op.add), ('*', op.mul),
1179 ('/', op.truediv), ('^', op.pow)]
1180 full_text = full_text.lstrip('(').rstrip(')')
1181 try:
1182 return float(full_text)
1183 except ValueError:
1184 for symbol, func in exprs:
1185 if symbol in full_text:
1186 left, right = full_text.split(symbol, 1) # single split
1187 return func(eval_no_bracket_expr(left),
1188 eval_no_bracket_expr(right))
1190 while '(' in text:
1191 middle = middle_brackets(text)
1192 text = text.replace(middle, '{}'.format(eval_no_bracket_expr(middle)))
1194 return float(eval_no_bracket_expr(text))
1196###
1197# Input file writing
1198###
1201# Ordered and case insensitive
1202KEYS = Namelist((
1203 ('CONTROL', [
1204 'calculation', 'title', 'verbosity', 'restart_mode', 'wf_collect',
1205 'nstep', 'iprint', 'tstress', 'tprnfor', 'dt', 'outdir', 'wfcdir',
1206 'prefix', 'lkpoint_dir', 'max_seconds', 'etot_conv_thr',
1207 'forc_conv_thr', 'disk_io', 'pseudo_dir', 'tefield', 'dipfield',
1208 'lelfield', 'nberrycyc', 'lorbm', 'lberry', 'gdir', 'nppstr',
1209 'lfcpopt', 'monopole']),
1210 ('SYSTEM', [
1211 'ibrav', 'celldm', 'A', 'B', 'C', 'cosAB', 'cosAC', 'cosBC', 'nat',
1212 'ntyp', 'nbnd', 'tot_charge', 'tot_magnetization',
1213 'starting_magnetization', 'ecutwfc', 'ecutrho', 'ecutfock', 'nr1',
1214 'nr2', 'nr3', 'nr1s', 'nr2s', 'nr3s', 'nosym', 'nosym_evc', 'noinv',
1215 'no_t_rev', 'force_symmorphic', 'use_all_frac', 'occupations',
1216 'one_atom_occupations', 'starting_spin_angle', 'degauss', 'smearing',
1217 'nspin', 'noncolin', 'ecfixed', 'qcutz', 'q2sigma', 'input_dft',
1218 'exx_fraction', 'screening_parameter', 'exxdiv_treatment',
1219 'x_gamma_extrapolation', 'ecutvcut', 'nqx1', 'nqx2', 'nqx3',
1220 'lda_plus_u', 'lda_plus_u_kind', 'Hubbard_U', 'Hubbard_J0',
1221 'Hubbard_alpha', 'Hubbard_beta', 'Hubbard_J',
1222 'starting_ns_eigenvalue', 'U_projection_type', 'edir',
1223 'emaxpos', 'eopreg', 'eamp', 'angle1', 'angle2',
1224 'constrained_magnetization', 'fixed_magnetization', 'lambda',
1225 'report', 'lspinorb', 'assume_isolated', 'esm_bc', 'esm_w',
1226 'esm_efield', 'esm_nfit', 'fcp_mu', 'vdw_corr', 'london',
1227 'london_s6', 'london_c6', 'london_rvdw', 'london_rcut',
1228 'ts_vdw_econv_thr', 'ts_vdw_isolated', 'xdm', 'xdm_a1', 'xdm_a2',
1229 'space_group', 'uniqueb', 'origin_choice', 'rhombohedral', 'zmon',
1230 'realxz', 'block', 'block_1', 'block_2', 'block_height']),
1231 ('ELECTRONS', [
1232 'electron_maxstep', 'scf_must_converge', 'conv_thr', 'adaptive_thr',
1233 'conv_thr_init', 'conv_thr_multi', 'mixing_mode', 'mixing_beta',
1234 'mixing_ndim', 'mixing_fixed_ns', 'diagonalization', 'ortho_para',
1235 'diago_thr_init', 'diago_cg_maxiter', 'diago_david_ndim',
1236 'diago_full_acc', 'efield', 'efield_cart', 'efield_phase',
1237 'startingpot', 'startingwfc', 'tqr']),
1238 ('IONS', [
1239 'ion_dynamics', 'ion_positions', 'pot_extrapolation',
1240 'wfc_extrapolation', 'remove_rigid_rot', 'ion_temperature', 'tempw',
1241 'tolp', 'delta_t', 'nraise', 'refold_pos', 'upscale', 'bfgs_ndim',
1242 'trust_radius_max', 'trust_radius_min', 'trust_radius_ini', 'w_1',
1243 'w_2']),
1244 ('CELL', [
1245 'cell_dynamics', 'press', 'wmass', 'cell_factor', 'press_conv_thr',
1246 'cell_dofree'])))
1249# Number of valence electrons in the pseudopotentials recommended by
1250# http://materialscloud.org/sssp/. These are just used as a fallback for
1251# calculating initial magetization values which are given as a fraction
1252# of valence electrons.
1253SSSP_VALENCE = [
1254 0, 1.0, 2.0, 3.0, 4.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 3.0, 4.0,
1255 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0, 13.0, 14.0, 15.0, 16.0, 17.0,
1256 18.0, 19.0, 20.0, 13.0, 14.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0,
1257 13.0, 14.0, 15.0, 16.0, 17.0, 18.0, 19.0, 12.0, 13.0, 14.0, 15.0, 6.0,
1258 7.0, 18.0, 9.0, 10.0, 11.0, 12.0, 13.0, 14.0, 15.0, 16.0, 17.0, 18.0,
1259 19.0, 20.0, 21.0, 22.0, 23.0, 24.0, 25.0, 36.0, 27.0, 14.0, 15.0, 30.0,
1260 15.0, 32.0, 19.0, 12.0, 13.0, 14.0, 15.0, 16.0, 18.0]
1263def construct_namelist(parameters=None, warn=False, **kwargs):
1264 """
1265 Construct an ordered Namelist containing all the parameters given (as
1266 a dictionary or kwargs). Keys will be inserted into their appropriate
1267 section in the namelist and the dictionary may contain flat and nested
1268 structures. Any kwargs that match input keys will be incorporated into
1269 their correct section. All matches are case-insensitive, and returned
1270 Namelist object is a case-insensitive dict.
1272 If a key is not known to ase, but in a section within `parameters`,
1273 it will be assumed that it was put there on purpose and included
1274 in the output namelist. Anything not in a section will be ignored (set
1275 `warn` to True to see ignored keys).
1277 Keys with a dimension (e.g. Hubbard_U(1)) will be incorporated as-is
1278 so the `i` should be made to match the output.
1280 The priority of the keys is:
1281 kwargs[key] > parameters[key] > parameters[section][key]
1282 Only the highest priority item will be included.
1284 Parameters
1285 ----------
1286 parameters: dict
1287 Flat or nested set of input parameters.
1288 warn: bool
1289 Enable warnings for unused keys.
1291 Returns
1292 -------
1293 input_namelist: Namelist
1294 pw.x compatible namelist of input parameters.
1296 """
1297 # Convert everything to Namelist early to make case-insensitive
1298 if parameters is None:
1299 parameters = Namelist()
1300 else:
1301 # Maximum one level of nested dict
1302 # Don't modify in place
1303 parameters_namelist = Namelist()
1304 for key, value in parameters.items():
1305 if isinstance(value, dict):
1306 parameters_namelist[key] = Namelist(value)
1307 else:
1308 parameters_namelist[key] = value
1309 parameters = parameters_namelist
1311 # Just a dict
1312 kwargs = Namelist(kwargs)
1314 # Final parameter set
1315 input_namelist = Namelist()
1317 # Collect
1318 for section in KEYS:
1319 sec_list = Namelist()
1320 for key in KEYS[section]:
1321 # Check all three separately and pop them all so that
1322 # we can check for missing values later
1323 if key in parameters.get(section, {}):
1324 sec_list[key] = parameters[section].pop(key)
1325 if key in parameters:
1326 sec_list[key] = parameters.pop(key)
1327 if key in kwargs:
1328 sec_list[key] = kwargs.pop(key)
1330 # Check if there is a key(i) version (no extra parsing)
1331 for arg_key in parameters.get(section, {}):
1332 if arg_key.split('(')[0].strip().lower() == key.lower():
1333 sec_list[arg_key] = parameters[section].pop(arg_key)
1334 cp_parameters = parameters.copy()
1335 for arg_key in cp_parameters:
1336 if arg_key.split('(')[0].strip().lower() == key.lower():
1337 sec_list[arg_key] = parameters.pop(arg_key)
1338 cp_kwargs = kwargs.copy()
1339 for arg_key in cp_kwargs:
1340 if arg_key.split('(')[0].strip().lower() == key.lower():
1341 sec_list[arg_key] = kwargs.pop(arg_key)
1343 # Add to output
1344 input_namelist[section] = sec_list
1346 unused_keys = list(kwargs)
1347 # pass anything else already in a section
1348 for key, value in parameters.items():
1349 if key in KEYS and isinstance(value, dict):
1350 input_namelist[key].update(value)
1351 elif isinstance(value, dict):
1352 unused_keys.extend(list(value))
1353 else:
1354 unused_keys.append(key)
1356 if warn and unused_keys:
1357 warnings.warn('Unused keys: {}'.format(', '.join(unused_keys)))
1359 return input_namelist
1362def grep_valence(pseudopotential):
1363 """
1364 Given a UPF pseudopotential file, find the number of valence atoms.
1366 Parameters
1367 ----------
1368 pseudopotential: str
1369 Filename of the pseudopotential.
1371 Returns
1372 -------
1373 valence: float
1374 Valence as reported in the pseudopotential.
1376 Raises
1377 ------
1378 ValueError
1379 If valence cannot be found in the pseudopotential.
1380 """
1382 # Example lines
1383 # Sr.pbe-spn-rrkjus_psl.1.0.0.UPF: z_valence="1.000000000000000E+001"
1384 # C.pbe-n-kjpaw_psl.1.0.0.UPF (new ld1.x):
1385 # ...PBC" z_valence="4.000000000000e0" total_p...
1386 # C_ONCV_PBE-1.0.upf: z_valence=" 4.00"
1387 # Ta_pbe_v1.uspp.F.UPF: 13.00000000000 Z valence
1389 with open(pseudopotential) as psfile:
1390 for line in psfile:
1391 if 'z valence' in line.lower():
1392 return float(line.split()[0])
1393 elif 'z_valence' in line.lower():
1394 if line.split()[0] == '<PP_HEADER':
1395 line = list(filter(lambda x: 'z_valence' in x,
1396 line.split(' ')))[0]
1397 return float(line.split('=')[-1].strip().strip('"'))
1398 else:
1399 raise ValueError('Valence missing in {}'.format(pseudopotential))
1402def cell_to_ibrav(cell, ibrav):
1403 """
1404 Calculate the appropriate `celldm(..)` parameters for the given ibrav
1405 using the given cell. The units for `celldm(..)` are Bohr.
1407 Does minimal checking of the cell shape, so it is possible to create
1408 a nonsense structure if the ibrav is inapproprite for the cell. These
1409 are derived to be symmetric with the routine for constructing the cell
1410 from ibrav parameters so directions of some vectors may be unexpected.
1412 Parameters
1413 ----------
1414 cell : np.array
1415 A 3x3 representation of a unit cell
1416 ibrav : int
1417 Bravais-lattice index according to the pw.x designations.
1419 Returns
1420 -------
1421 parameters : dict
1422 A dictionary with all the necessary `celldm(..)` keys assigned
1423 necessary values (in units of Bohr). Also includes `ibrav` so it
1424 can be passed back to `ibrav_to_cell`.
1426 Raises
1427 ------
1428 NotImplementedError
1429 Only a limited number of ibrav settings can be parsed. An error
1430 is raised if the ibrav interpretation is not implemented.
1431 """
1432 parameters = {'ibrav': ibrav}
1434 if ibrav == 1:
1435 parameters['celldm(1)'] = cell[0][0] / units['Bohr']
1436 elif ibrav in [2, 3, -3]:
1437 parameters['celldm(1)'] = cell[0][2] * 2 / units['Bohr']
1438 elif ibrav in [4, 6]:
1439 parameters['celldm(1)'] = cell[0][0] / units['Bohr']
1440 parameters['celldm(3)'] = cell[2][2] / cell[0][0]
1441 elif ibrav in [5, -5]:
1442 # Manually derive
1443 a = np.linalg.norm(cell[0])
1444 cosab = np.dot(cell[0], cell[1]) / (a ** 2)
1445 parameters['celldm(1)'] = a / units['Bohr']
1446 parameters['celldm(4)'] = cosab
1447 elif ibrav == 7:
1448 parameters['celldm(1)'] = cell[0][0] * 2 / units['Bohr']
1449 parameters['celldm(3)'] = cell[2][2] / cell[0][0]
1450 elif ibrav == 8:
1451 parameters['celldm(1)'] = cell[0][0] / units['Bohr']
1452 parameters['celldm(2)'] = cell[1][1] / cell[0][0]
1453 parameters['celldm(3)'] = cell[2][2] / cell[0][0]
1454 elif ibrav in [9, -9]:
1455 parameters['celldm(1)'] = cell[0][0] * 2 / units['Bohr']
1456 parameters['celldm(2)'] = cell[1][1] / cell[0][0]
1457 parameters['celldm(3)'] = cell[2][2] * 2 / cell[0][0]
1458 elif ibrav in [10, 11]:
1459 parameters['celldm(1)'] = cell[0][0] * 2 / units['Bohr']
1460 parameters['celldm(2)'] = cell[1][1] / cell[0][0]
1461 parameters['celldm(3)'] = cell[2][2] / cell[0][0]
1462 elif ibrav == 12:
1463 # cos^2 + sin^2
1464 b = (cell[1][0]**2 + cell[1][1]**2)**0.5
1465 parameters['celldm(1)'] = cell[0][0] / units['Bohr']
1466 parameters['celldm(2)'] = b / cell[0][0]
1467 parameters['celldm(3)'] = cell[2][2] / cell[0][0]
1468 parameters['celldm(4)'] = cell[1][0] / b
1469 elif ibrav == -12:
1470 # cos^2 + sin^2
1471 c = (cell[2][0]**2 + cell[2][2]**2)**0.5
1472 parameters['celldm(1)'] = cell[0][0] / units['Bohr']
1473 parameters['celldm(2)'] = cell[1][1] / cell[0][0]
1474 parameters['celldm(3)'] = c / cell[0][0]
1475 parameters['celldm(4)'] = cell[2][0] / c
1476 elif ibrav == 13:
1477 b = (cell[1][0]**2 + cell[1][1]**2)**0.5
1478 parameters['celldm(1)'] = cell[0][0] * 2 / units['Bohr']
1479 parameters['celldm(2)'] = b / (cell[0][0] * 2)
1480 parameters['celldm(3)'] = cell[2][2] / cell[0][0]
1481 parameters['celldm(4)'] = cell[1][0] / b
1482 elif ibrav == 14:
1483 # Manually derive
1484 a, b, c = np.linalg.norm(cell, axis=1)
1485 cosbc = np.dot(cell[1], cell[2]) / (b * c)
1486 cosac = np.dot(cell[0], cell[2]) / (a * c)
1487 cosab = np.dot(cell[0], cell[1]) / (a * b)
1488 parameters['celldm(1)'] = a / units['Bohr']
1489 parameters['celldm(2)'] = b / a
1490 parameters['celldm(3)'] = c / a
1491 parameters['celldm(4)'] = cosbc
1492 parameters['celldm(5)'] = cosac
1493 parameters['celldm(6)'] = cosab
1494 else:
1495 raise NotImplementedError('ibrav = {0} is not implemented'
1496 ''.format(ibrav))
1498 return parameters
1501def kspacing_to_grid(atoms, spacing, calculated_spacing=None):
1502 """
1503 Calculate the kpoint mesh that is equivalent to the given spacing
1504 in reciprocal space (units Angstrom^-1). The number of kpoints is each
1505 dimension is rounded up (compatible with CASTEP).
1507 Parameters
1508 ----------
1509 atoms: ase.Atoms
1510 A structure that can have get_reciprocal_cell called on it.
1511 spacing: float
1512 Minimum K-Point spacing in $A^{-1}$.
1513 calculated_spacing : list
1514 If a three item list (or similar mutable sequence) is given the
1515 members will be replaced with the actual calculated spacing in
1516 $A^{-1}$.
1518 Returns
1519 -------
1520 kpoint_grid : [int, int, int]
1521 MP grid specification to give the required spacing.
1523 """
1524 # No factor of 2pi in ase, everything in A^-1
1525 # reciprocal dimensions
1526 r_x, r_y, r_z = np.linalg.norm(atoms.cell.reciprocal(), axis=1)
1528 kpoint_grid = [int(r_x / spacing) + 1,
1529 int(r_y / spacing) + 1,
1530 int(r_z / spacing) + 1]
1532 for i, _ in enumerate(kpoint_grid):
1533 if not atoms.pbc[i]:
1534 kpoint_grid[i] = 1
1536 if calculated_spacing is not None:
1537 calculated_spacing[:] = [r_x / kpoint_grid[0],
1538 r_y / kpoint_grid[1],
1539 r_z / kpoint_grid[2]]
1541 return kpoint_grid
1544def write_espresso_in(fd, atoms, input_data=None, pseudopotentials=None,
1545 kspacing=None, kpts=None, koffset=(0, 0, 0),
1546 crystal_coordinates=False, **kwargs):
1547 """
1548 Create an input file for pw.x.
1550 Use set_initial_magnetic_moments to turn on spin, if ispin is set to 2
1551 with no magnetic moments, they will all be set to 0.0. Magnetic moments
1552 will be converted to the QE units (fraction of valence electrons) using
1553 any pseudopotential files found, or a best guess for the number of
1554 valence electrons.
1556 Units are not converted for any other input data, so use Quantum ESPRESSO
1557 units (Usually Ry or atomic units).
1559 Keys with a dimension (e.g. Hubbard_U(1)) will be incorporated as-is
1560 so the `i` should be made to match the output.
1562 Implemented features:
1564 - Conversion of :class:`ase.constraints.FixAtoms` and
1565 :class:`ase.constraints.FixCartesian`.
1566 - `starting_magnetization` derived from the `mgmoms` and pseudopotentials
1567 (searches default paths for pseudo files.)
1568 - Automatic assignment of options to their correct sections.
1569 - Interpretation of ibrav (cell must exactly match the vectors defined
1570 in the QE docs).
1572 Not implemented:
1574 - Lists of k-points
1575 - Other constraints
1576 - Hubbard parameters
1577 - Validation of the argument types for input
1578 - Validation of required options
1579 - Reorientation for ibrav settings
1580 - Noncollinear magnetism
1582 Parameters
1583 ----------
1584 fd: file
1585 A file like object to write the input file to.
1586 atoms: Atoms
1587 A single atomistic configuration to write to `fd`.
1588 input_data: dict
1589 A flat or nested dictionary with input parameters for pw.x
1590 pseudopotentials: dict
1591 A filename for each atomic species, e.g.
1592 {'O': 'O.pbe-rrkjus.UPF', 'H': 'H.pbe-rrkjus.UPF'}.
1593 A dummy name will be used if none are given.
1594 kspacing: float
1595 Generate a grid of k-points with this as the minimum distance,
1596 in A^-1 between them in reciprocal space. If set to None, kpts
1597 will be used instead.
1598 kpts: (int, int, int) or dict
1599 If kpts is a tuple (or list) of 3 integers, it is interpreted
1600 as the dimensions of a Monkhorst-Pack grid.
1601 If ``kpts`` is set to ``None``, only the Γ-point will be included
1602 and QE will use routines optimized for Γ-point-only calculations.
1603 Compared to Γ-point-only calculations without this optimization
1604 (i.e. with ``kpts=(1, 1, 1)``), the memory and CPU requirements
1605 are typically reduced by half.
1606 If kpts is a dict, it will either be interpreted as a path
1607 in the Brillouin zone (*) if it contains the 'path' keyword,
1608 otherwise it is converted to a Monkhorst-Pack grid (**).
1609 (*) see ase.dft.kpoints.bandpath
1610 (**) see ase.calculators.calculator.kpts2sizeandoffsets
1611 koffset: (int, int, int)
1612 Offset of kpoints in each direction. Must be 0 (no offset) or
1613 1 (half grid offset). Setting to True is equivalent to (1, 1, 1).
1614 crystal_coordinates: bool
1615 Whether the atomic positions should be written to the QE input file in
1616 absolute (False, default) or relative (crystal) coordinates (True).
1618 """
1620 # Convert to a namelist to make working with parameters much easier
1621 # Note that the name ``input_data`` is chosen to prevent clash with
1622 # ``parameters`` in Calculator objects
1623 input_parameters = construct_namelist(input_data, **kwargs)
1625 # Convert ase constraints to QE constraints
1626 # Nx3 array of force multipliers matches what QE uses
1627 # Do this early so it is available when constructing the atoms card
1628 constraint_mask = np.ones((len(atoms), 3), dtype='int')
1629 for constraint in atoms.constraints:
1630 if isinstance(constraint, FixAtoms):
1631 constraint_mask[constraint.index] = 0
1632 elif isinstance(constraint, FixCartesian):
1633 constraint_mask[constraint.a] = constraint.mask
1634 else:
1635 warnings.warn('Ignored unknown constraint {}'.format(constraint))
1637 # Species info holds the information on the pseudopotential and
1638 # associated for each element
1639 if pseudopotentials is None:
1640 pseudopotentials = {}
1641 species_info = {}
1642 for species in set(atoms.get_chemical_symbols()):
1643 # Look in all possible locations for the pseudos and try to figure
1644 # out the number of valence electrons
1645 pseudo = pseudopotentials.get(species, None)
1646 valence = get_valence_electrons(species, input_parameters, pseudo)
1647 species_info[species] = {'pseudo': pseudo, 'valence': valence}
1649 # Convert atoms into species.
1650 # Each different magnetic moment needs to be a separate type even with
1651 # the same pseudopotential (e.g. an up and a down for AFM).
1652 # if any magmom are > 0 or nspin == 2 then use species labels.
1653 # Rememeber: magnetisation uses 1 based indexes
1654 atomic_species = OrderedDict()
1655 atomic_species_str = []
1656 atomic_positions_str = []
1658 nspin = input_parameters['system'].get('nspin', 1) # 1 is the default
1659 if any(atoms.get_initial_magnetic_moments()):
1660 if nspin == 1:
1661 # Force spin on
1662 input_parameters['system']['nspin'] = 2
1663 nspin = 2
1665 if nspin == 2:
1666 # Spin on
1667 for atom, magmom in zip(atoms, atoms.get_initial_magnetic_moments()):
1668 if (atom.symbol, magmom) not in atomic_species:
1669 # spin as fraction of valence
1670 fspin = float(magmom) / species_info[atom.symbol]['valence']
1671 # Index in the atomic species list
1672 sidx = len(atomic_species) + 1
1673 # Index for that atom type; no index for first one
1674 tidx = sum(atom.symbol == x[0] for x in atomic_species) or ' '
1675 atomic_species[(atom.symbol, magmom)] = (sidx, tidx)
1676 # Add magnetization to the input file
1677 mag_str = 'starting_magnetization({0})'.format(sidx)
1678 input_parameters['system'][mag_str] = fspin
1679 atomic_species_str.append(
1680 '{species}{tidx} {mass} {pseudo}\n'.format(
1681 species=atom.symbol, tidx=tidx, mass=atom.mass,
1682 pseudo=species_info[atom.symbol]['pseudo']))
1683 # lookup tidx to append to name
1684 sidx, tidx = atomic_species[(atom.symbol, magmom)]
1686 # only inclued mask if something is fixed
1687 if not all(constraint_mask[atom.index]):
1688 mask = ' {mask[0]} {mask[1]} {mask[2]}'.format(
1689 mask=constraint_mask[atom.index])
1690 else:
1691 mask = ''
1693 # construct line for atomic positions
1694 atomic_positions_str.append(
1695 '{atom.symbol}{tidx} '
1696 '{atom.x:.10f} {atom.y:.10f} {atom.z:.10f}'
1697 '{mask}\n'.format(atom=atom, tidx=tidx, mask=mask))
1699 else:
1700 # Do nothing about magnetisation
1701 for atom in atoms:
1702 if atom.symbol not in atomic_species:
1703 atomic_species[atom.symbol] = True # just a placeholder
1704 atomic_species_str.append(
1705 '{species} {mass} {pseudo}\n'.format(
1706 species=atom.symbol, mass=atom.mass,
1707 pseudo=species_info[atom.symbol]['pseudo']))
1709 # only inclued mask if something is fixed
1710 if not all(constraint_mask[atom.index]):
1711 mask = ' {mask[0]} {mask[1]} {mask[2]}'.format(
1712 mask=constraint_mask[atom.index])
1713 else:
1714 mask = ''
1716 if crystal_coordinates:
1717 coords = [atom.a, atom.b, atom.c]
1718 else:
1719 coords = atom.position
1720 atomic_positions_str.append(
1721 '{atom.symbol} '
1722 '{coords[0]:.10f} {coords[1]:.10f} {coords[2]:.10f} '
1723 '{mask}\n'.format(atom=atom, coords=coords, mask=mask))
1725 # Add computed parameters
1726 # different magnetisms means different types
1727 input_parameters['system']['ntyp'] = len(atomic_species)
1728 input_parameters['system']['nat'] = len(atoms)
1730 # Use cell as given or fit to a specific ibrav
1731 if 'ibrav' in input_parameters['system']:
1732 ibrav = input_parameters['system']['ibrav']
1733 if ibrav != 0:
1734 celldm = cell_to_ibrav(atoms.cell, ibrav)
1735 regen_cell = ibrav_to_cell(celldm)[1]
1736 if not np.allclose(atoms.cell, regen_cell):
1737 warnings.warn('Input cell does not match requested ibrav'
1738 '{} != {}'.format(regen_cell, atoms.cell))
1739 input_parameters['system'].update(celldm)
1740 else:
1741 # Just use standard cell block
1742 input_parameters['system']['ibrav'] = 0
1744 # Construct input file into this
1745 pwi = []
1747 # Assume sections are ordered (taken care of in namelist construction)
1748 # and that repr converts to a QE readable representation (except bools)
1749 for section in input_parameters:
1750 pwi.append('&{0}\n'.format(section.upper()))
1751 for key, value in input_parameters[section].items():
1752 if value is True:
1753 pwi.append(' {0:16} = .true.\n'.format(key))
1754 elif value is False:
1755 pwi.append(' {0:16} = .false.\n'.format(key))
1756 else:
1757 # repr format to get quotes around strings
1758 pwi.append(' {0:16} = {1!r:}\n'.format(key, value))
1759 pwi.append('/\n') # terminate section
1760 pwi.append('\n')
1762 # Pseudopotentials
1763 pwi.append('ATOMIC_SPECIES\n')
1764 pwi.extend(atomic_species_str)
1765 pwi.append('\n')
1767 # KPOINTS - add a MP grid as required
1768 if kspacing is not None:
1769 kgrid = kspacing_to_grid(atoms, kspacing)
1770 elif kpts is not None:
1771 if isinstance(kpts, dict) and 'path' not in kpts:
1772 kgrid, shift = kpts2sizeandoffsets(atoms=atoms, **kpts)
1773 koffset = []
1774 for i, x in enumerate(shift):
1775 assert x == 0 or abs(x * kgrid[i] - 0.5) < 1e-14
1776 koffset.append(0 if x == 0 else 1)
1777 else:
1778 kgrid = kpts
1779 else:
1780 kgrid = "gamma"
1782 # True and False work here and will get converted by ':d' format
1783 if isinstance(koffset, int):
1784 koffset = (koffset, ) * 3
1786 # BandPath object or bandpath-as-dictionary:
1787 if isinstance(kgrid, dict) or hasattr(kgrid, 'kpts'):
1788 pwi.append('K_POINTS crystal_b\n')
1789 assert hasattr(kgrid, 'path') or 'path' in kgrid
1790 kgrid = kpts2ndarray(kgrid, atoms=atoms)
1791 pwi.append('%s\n' % len(kgrid))
1792 for k in kgrid:
1793 pwi.append('{k[0]:.14f} {k[1]:.14f} {k[2]:.14f} 0\n'.format(k=k))
1794 pwi.append('\n')
1795 elif isinstance(kgrid, str) and (kgrid == "gamma"):
1796 pwi.append('K_POINTS gamma\n')
1797 pwi.append('\n')
1798 else:
1799 pwi.append('K_POINTS automatic\n')
1800 pwi.append('{0[0]} {0[1]} {0[2]} {1[0]:d} {1[1]:d} {1[2]:d}\n'
1801 ''.format(kgrid, koffset))
1802 pwi.append('\n')
1804 # CELL block, if required
1805 if input_parameters['SYSTEM']['ibrav'] == 0:
1806 pwi.append('CELL_PARAMETERS angstrom\n')
1807 pwi.append('{cell[0][0]:.14f} {cell[0][1]:.14f} {cell[0][2]:.14f}\n'
1808 '{cell[1][0]:.14f} {cell[1][1]:.14f} {cell[1][2]:.14f}\n'
1809 '{cell[2][0]:.14f} {cell[2][1]:.14f} {cell[2][2]:.14f}\n'
1810 ''.format(cell=atoms.cell))
1811 pwi.append('\n')
1813 # Positions - already constructed, but must appear after namelist
1814 if crystal_coordinates:
1815 pwi.append('ATOMIC_POSITIONS crystal\n')
1816 else:
1817 pwi.append('ATOMIC_POSITIONS angstrom\n')
1818 pwi.extend(atomic_positions_str)
1819 pwi.append('\n')
1821 # DONE!
1822 fd.write(''.join(pwi))