Coverage for ase / io / espresso.py: 77.82%
879 statements
« prev ^ index » next coverage.py v7.14.0, created at 2026-05-21 15:52 +0000
1# fmt: off
3"""Reads Quantum ESPRESSO files.
5Read multiple structures and results from pw.x output files. Read
6structures from pw.x input files.
8Built for PWSCF v.5.3.0 but should work with earlier and later versions.
9Can deal with most major functionality, with the notable exception of ibrav,
10for which we only support ibrav == 0 and force CELL_PARAMETERS to be provided
11explicitly.
13Units are converted using CODATA 2006, as used internally by Quantum
14ESPRESSO.
15"""
17import operator as op
18import re
19import warnings
20from collections import defaultdict
21from copy import deepcopy
22from pathlib import Path
24import numpy as np
26from ase.atoms import Atoms
27from ase.calculators.calculator import kpts2ndarray, kpts2sizeandoffsets
28from ase.calculators.singlepoint import (
29 SinglePointDFTCalculator,
30 SinglePointKPoint,
31)
32from ase.constraints import FixAtoms, FixCartesian
33from ase.data import chemical_symbols
34from ase.dft.kpoints import kpoint_convert
35from ase.io.espresso_namelist.keys import pw_keys
36from ase.io.espresso_namelist.namelist import Namelist
37from ase.units import create_units
38from ase.utils import deprecated, reader, writer
40# Quantum ESPRESSO uses CODATA 2006 internally
41units = create_units('2006')
43# Section identifiers
44_PW_START = 'Program PWSCF'
45_PW_END = 'End of self-consistent calculation'
46_PW_CELL = 'CELL_PARAMETERS'
47_PW_POS = 'ATOMIC_POSITIONS'
48_PW_MAGMOM = 'Magnetic moment per site'
49_PW_FORCE = 'Forces acting on atoms'
50_PW_TOTEN = '! total energy'
51_PW_STRESS = 'total stress'
52_PW_FERMI = 'the Fermi energy is'
53_PW_HIGHEST_OCCUPIED = 'highest occupied level'
54_PW_HIGHEST_OCCUPIED_LOWEST_FREE = 'highest occupied, lowest unoccupied level'
55_PW_KPTS = 'number of k points='
56_PW_BANDS = _PW_END
57_PW_BANDSTRUCTURE = 'End of band structure calculation'
58_PW_DIPOLE = "Debye"
59_PW_DIPOLE_DIRECTION = "Computed dipole along edir"
61# ibrav error message
62ibrav_error_message = (
63 'ASE does not support ibrav != 0. Note that with ibrav '
64 '== 0, Quantum ESPRESSO will still detect the symmetries '
65 'of your system because the CELL_PARAMETERS are defined '
66 'to a high level of precision.')
69@reader
70def read_espresso_out(fileobj, index=slice(None), results_required=True):
71 """Reads Quantum ESPRESSO output files.
73 The atomistic configurations as well as results (energy, force, stress,
74 magnetic moments) of the calculation are read for all configurations
75 within the output file.
77 Will probably raise errors for broken or incomplete files.
79 Parameters
80 ----------
81 fileobj : file|str
82 A file like object or filename
83 index : slice
84 The index of configurations to extract.
85 results_required : bool
86 If True, atomistic configurations that do not have any
87 associated results will not be included. This prevents double
88 printed configurations and incomplete calculations from being
89 returned as the final configuration with no results data.
91 Yields
92 ------
93 structure : Atoms
94 The next structure from the index slice. The Atoms has a
95 SinglePointCalculator attached with any results parsed from
96 the file.
99 """
100 # work with a copy in memory for faster random access
101 pwo_lines = fileobj.readlines()
103 # TODO: index -1 special case?
104 # Index all the interesting points
105 indexes = {
106 _PW_START: [],
107 _PW_END: [],
108 _PW_CELL: [],
109 _PW_POS: [],
110 _PW_MAGMOM: [],
111 _PW_FORCE: [],
112 _PW_TOTEN: [],
113 _PW_STRESS: [],
114 _PW_FERMI: [],
115 _PW_HIGHEST_OCCUPIED: [],
116 _PW_HIGHEST_OCCUPIED_LOWEST_FREE: [],
117 _PW_KPTS: [],
118 _PW_BANDS: [],
119 _PW_BANDSTRUCTURE: [],
120 _PW_DIPOLE: [],
121 _PW_DIPOLE_DIRECTION: [],
122 }
124 for idx, line in enumerate(pwo_lines):
125 for identifier in indexes:
126 if identifier in line:
127 indexes[identifier].append(idx)
129 # Configurations are either at the start, or defined in ATOMIC_POSITIONS
130 # in a subsequent step. Can deal with concatenated output files.
131 all_config_indexes = sorted(indexes[_PW_START] +
132 indexes[_PW_POS])
134 # Slice only requested indexes
135 # setting results_required argument stops configuration-only
136 # structures from being returned. This ensures the [-1] structure
137 # is one that has results. Two cases:
138 # - SCF of last configuration is not converged, job terminated
139 # abnormally.
140 # - 'relax' and 'vc-relax' re-prints the final configuration but
141 # only 'vc-relax' recalculates.
142 if results_required:
143 results_indexes = sorted(indexes[_PW_TOTEN] + indexes[_PW_FORCE] +
144 indexes[_PW_STRESS] + indexes[_PW_MAGMOM] +
145 indexes[_PW_BANDS] +
146 indexes[_PW_BANDSTRUCTURE])
148 # Prune to only configurations with results data before the next
149 # configuration
150 results_config_indexes = []
151 for config_index, config_index_next in zip(
152 all_config_indexes,
153 all_config_indexes[1:] + [len(pwo_lines)]):
154 if any(config_index < results_index < config_index_next
155 for results_index in results_indexes):
156 results_config_indexes.append(config_index)
158 # slice from the subset
159 image_indexes = results_config_indexes[index]
160 else:
161 image_indexes = all_config_indexes[index]
163 # Extract initialisation information each time PWSCF starts
164 # to add to subsequent configurations. Use None so slices know
165 # when to fill in the blanks.
166 pwscf_start_info = {idx: None for idx in indexes[_PW_START]}
168 for image_index in image_indexes:
169 # Find the nearest calculation start to parse info. Needed in,
170 # for example, relaxation where cell is only printed at the
171 # start.
172 if image_index in indexes[_PW_START]:
173 prev_start_index = image_index
174 else:
175 # The greatest start index before this structure
176 prev_start_index = [idx for idx in indexes[_PW_START]
177 if idx < image_index][-1]
179 # add structure to reference if not there
180 if pwscf_start_info[prev_start_index] is None:
181 pwscf_start_info[prev_start_index] = parse_pwo_start(
182 pwo_lines, prev_start_index)
184 # Get the bounds for information for this structure. Any associated
185 # values will be between the image_index and the following one,
186 # EXCEPT for cell, which will be 4 lines before if it exists.
187 for next_index in all_config_indexes:
188 if next_index > image_index:
189 break
190 else:
191 # right to the end of the file
192 next_index = len(pwo_lines)
194 # Get the structure
195 # Use this for any missing data
196 prev_structure = pwscf_start_info[prev_start_index]['atoms']
197 cell_alat = pwscf_start_info[prev_start_index]['alat']
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, _ = 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 def find_thing(get_thing, indices, **kwargs):
224 for index in indices:
225 if image_index < index < next_index:
226 return get_thing(pwo_lines, index, **kwargs)
227 return None
229 natoms = len(structure)
231 energy = find_thing(_get_energy, indexes[_PW_TOTEN])
232 forces = find_thing(_get_forces, indexes[_PW_FORCE], natoms=natoms)
233 stress = find_thing(_get_stress, indexes[_PW_STRESS])
234 magmoms = find_thing(_get_magmoms, indexes[_PW_MAGMOM], natoms=natoms)
236 # Dipole moment
237 dipole = None
238 if indexes[_PW_DIPOLE]:
239 for dipole_index in indexes[_PW_DIPOLE]:
240 if image_index < dipole_index < next_index:
241 _dipole = float(pwo_lines[dipole_index].split()[-2])
243 for dipole_index in indexes[_PW_DIPOLE_DIRECTION]:
244 if image_index < dipole_index < next_index:
245 _direction = pwo_lines[dipole_index].strip()
246 prefix = 'Computed dipole along edir('
247 _direction = _direction[len(prefix):]
248 _direction = int(_direction[0])
250 dipole = np.eye(3)[_direction - 1] * _dipole * units['Debye']
252 # Fermi level / highest occupied level
253 efermi = None
254 for fermi_index in indexes[_PW_FERMI]:
255 if image_index < fermi_index < next_index:
256 efermi = float(pwo_lines[fermi_index].split()[-2])
258 if efermi is None:
259 for ho_index in indexes[_PW_HIGHEST_OCCUPIED]:
260 if image_index < ho_index < next_index:
261 efermi = float(pwo_lines[ho_index].split()[-1])
263 if efermi is None:
264 for holf_index in indexes[_PW_HIGHEST_OCCUPIED_LOWEST_FREE]:
265 if image_index < holf_index < next_index:
266 efermi = float(pwo_lines[holf_index].split()[-2])
268 # K-points
269 ibzkpts = None
270 weights = None
271 kpoints_warning = "Number of k-points >= 100: " + \
272 "set verbosity='high' to print them."
274 for kpts_index in indexes[_PW_KPTS]:
275 nkpts = int(re.findall(r'\b\d+\b', pwo_lines[kpts_index])[0])
276 kpts_index += 2
278 if pwo_lines[kpts_index].strip() == kpoints_warning:
279 continue
281 # QE prints the k-points in units of 2*pi/alat
282 cell = structure.get_cell()
283 ibzkpts = []
284 weights = []
285 for i in range(nkpts):
286 L = pwo_lines[kpts_index + i].split()
287 weights.append(float(L[-1]))
288 coord = np.array([L[-6], L[-5], L[-4].strip('),')],
289 dtype=float)
290 coord *= 2 * np.pi / cell_alat
291 coord = kpoint_convert(cell, ckpts_kv=coord)
292 ibzkpts.append(coord)
293 ibzkpts = np.array(ibzkpts)
294 weights = np.array(weights)
296 # Bands
297 kpts = None
298 kpoints_warning = "Number of k-points >= 100: " + \
299 "set verbosity='high' to print the bands."
301 for bands_index in indexes[_PW_BANDS] + indexes[_PW_BANDSTRUCTURE]:
302 if image_index < bands_index < next_index:
303 bands_index += 1
304 # skip over the lines with DFT+U occupation matrices
305 if 'enter write_ns' in pwo_lines[bands_index]:
306 while 'exit write_ns' not in pwo_lines[bands_index]:
307 bands_index += 1
308 bands_index += 1
310 if pwo_lines[bands_index].strip() == kpoints_warning:
311 continue
313 assert ibzkpts is not None
314 spin, bands, eigenvalues = 0, [], [[], []]
316 while True:
317 L = pwo_lines[bands_index].replace('-', ' -').split()
318 if len(L) == 0:
319 if len(bands) > 0:
320 eigenvalues[spin].append(bands)
321 bands = []
322 elif L == ['occupation', 'numbers']:
323 # Skip the lines with the occupation numbers
324 bands_index += len(eigenvalues[spin][0]) // 8 + 1
325 elif L[0] == 'k' and L[1].startswith('='):
326 pass
327 elif 'SPIN' in L:
328 if 'DOWN' in L:
329 spin += 1
330 else:
331 try:
332 bands.extend(map(float, L))
333 except ValueError:
334 break
335 bands_index += 1
337 if spin == 1:
338 assert len(eigenvalues[0]) == len(eigenvalues[1])
339 assert len(eigenvalues[0]) == len(ibzkpts), \
340 (np.shape(eigenvalues), len(ibzkpts))
342 kpts = []
343 for s in range(spin + 1):
344 for w, k, e in zip(weights, ibzkpts, eigenvalues[s]):
345 kpt = SinglePointKPoint(w, s, k, eps_n=e)
346 kpts.append(kpt)
348 # Put everything together
349 #
350 # In PW the forces are consistent with the "total energy"; that's why
351 # its value must be assigned to free_energy.
352 # PW doesn't compute the extrapolation of the energy to 0K smearing
353 # the closer thing to this is again the total energy that contains
354 # the correct (i.e. variational) form of the band energy is
355 # Eband = \int e N(e) de for e<Ef , where N(e) is the DOS
356 # This differs by the term (-TS) from the sum of KS eigenvalues:
357 # Eks = \sum wg(n,k) et(n,k)
358 # which is non variational. When a Fermi-Dirac function is used
359 # for a given T, the variational energy is REALLY the free energy F,
360 # and F = E - TS , with E = non variational energy.
361 #
362 calc = SinglePointDFTCalculator(structure, energy=energy,
363 free_energy=energy,
364 forces=forces, stress=stress,
365 magmoms=magmoms, efermi=efermi,
366 ibzkpts=ibzkpts, dipole=dipole)
367 calc.kpts = kpts
368 structure.calc = calc
370 yield structure
373def parse_pwo_start(lines, index=0):
374 """Parse Quantum ESPRESSO calculation info from lines,
375 starting from index. Return a dictionary containing extracted
376 information.
378 - `celldm(1)`: lattice parameters (alat)
379 - `cell`: unit cell in Angstrom
380 - `symbols`: element symbols for the structure
381 - `positions`: cartesian coordinates of atoms in Angstrom
382 - `atoms`: an `ase.Atoms` object constructed from the extracted data
384 Parameters
385 ----------
386 lines : list[str]
387 Contents of PWSCF output file.
388 index : int
389 Line number to begin parsing. Only first calculation will
390 be read.
392 Returns
393 -------
394 info : dict
395 Dictionary of calculation parameters, including `celldm(1)`, `cell`,
396 `symbols`, `positions`, `atoms`.
398 Raises
399 ------
400 KeyError
401 If interdependent values cannot be found (especially celldm(1))
402 an error will be raised as other quantities cannot then be
403 calculated (e.g. cell and positions).
404 """
405 # TODO: extend with extra DFT info?
407 info = {}
409 for idx, line in enumerate(lines[index:], start=index):
410 if 'celldm(1)' in line:
411 # celldm(1) has more digits than alat!!
412 info['celldm(1)'] = float(line.split()[1]) * units['Bohr']
413 info['alat'] = info['celldm(1)']
414 elif 'number of atoms/cell' in line:
415 info['nat'] = int(line.split()[-1])
416 elif 'number of atomic types' in line:
417 info['ntyp'] = int(line.split()[-1])
418 elif 'crystal axes:' in line:
419 info['cell'] = info['celldm(1)'] * np.array([
420 [float(x) for x in lines[idx + 1].split()[3:6]],
421 [float(x) for x in lines[idx + 2].split()[3:6]],
422 [float(x) for x in lines[idx + 3].split()[3:6]]])
423 elif 'positions (alat units)' in line:
424 info['symbols'], info['positions'] = [], []
426 for at_line in lines[idx + 1:idx + 1 + info['nat']]:
427 sym, x, y, z = parse_position_line(at_line)
428 info['symbols'].append(label_to_symbol(sym))
429 info['positions'].append([x * info['celldm(1)'],
430 y * info['celldm(1)'],
431 z * info['celldm(1)']])
432 # This should be the end of interesting info.
433 # Break here to avoid dealing with large lists of kpoints.
434 # Will need to be extended for DFTCalculator info.
435 break
437 # Make atoms for convenience
438 info['atoms'] = Atoms(symbols=info['symbols'],
439 positions=info['positions'],
440 cell=info['cell'], pbc=True)
442 return info
445def parse_position_line(line):
446 """Parse a single line from a pw.x output file.
448 The line must contain information about the atomic symbol and the position,
449 e.g.
451 995 Sb tau( 995) = ( 1.4212023 0.7037863 0.1242640 )
453 Parameters
454 ----------
455 line : str
456 Line to be parsed.
458 Returns
459 -------
460 sym : str
461 Atomic symbol.
462 x : float
463 x-position.
464 y : float
465 y-position.
466 z : float
467 z-position.
468 """
469 pat = re.compile(r'\s*\d+\s*(\S+)\s*tau\(\s*\d+\)\s*='
470 r'\s*\(\s*(\S+)\s+(\S+)\s+(\S+)\s*\)')
471 match = pat.match(line)
472 assert match is not None
473 sym, x, y, z = match.group(1, 2, 3, 4)
474 return sym, float(x), float(y), float(z)
477@reader
478def read_espresso_in(fileobj):
479 """Parse a Quantum ESPRESSO input files, '.in', '.pwi'.
481 ESPRESSO inputs are generally a fortran-namelist format with custom
482 blocks of data. The namelist is parsed as a dict and an atoms object
483 is constructed from the included information.
485 Parameters
486 ----------
487 fileobj : file | str
488 A file-like object that supports line iteration with the contents
489 of the input file, or a filename.
491 Returns
492 -------
493 atoms : Atoms
494 Structure defined in the input file.
496 Raises
497 ------
498 KeyError
499 Raised for missing keys that are required to process the file
500 """
501 # parse namelist section and extract remaining lines
502 data, card_lines = read_fortran_namelist(fileobj)
504 # get the cell if ibrav=0
505 if 'system' not in data:
506 raise KeyError('Required section &SYSTEM not found.')
507 elif 'ibrav' not in data['system']:
508 raise KeyError('ibrav is required in &SYSTEM')
509 elif data['system']['ibrav'] == 0:
510 # celldm(1) is in Bohr, A is in angstrom. celldm(1) will be
511 # used even if A is also specified.
512 if 'celldm(1)' in data['system']:
513 alat = data['system']['celldm(1)'] * units['Bohr']
514 elif 'A' in data['system']:
515 alat = data['system']['A']
516 else:
517 alat = None
518 cell, _ = get_cell_parameters(card_lines, alat=alat)
519 else:
520 raise ValueError(ibrav_error_message)
522 # species_info holds some info for each element
523 species_card = get_atomic_species(
524 card_lines, n_species=data['system']['ntyp'])
525 species_info = {}
526 for ispec, (label, weight, pseudo) in enumerate(species_card):
527 symbol = label_to_symbol(label)
529 # starting_magnetization is in fractions of valence electrons
530 magnet_key = f"starting_magnetization({ispec + 1})"
531 magmom = data["system"].get(magnet_key, 0.0)
532 species_info[symbol] = {"weight": weight, "pseudo": pseudo,
533 "magmom": magmom}
535 positions_card = get_atomic_positions(
536 card_lines, n_atoms=data['system']['nat'], cell=cell, alat=alat)
538 symbols = [label_to_symbol(position[0]) for position in positions_card]
539 positions = [position[1] for position in positions_card]
540 constraint_flags = [position[2] for position in positions_card]
541 magmoms = [species_info[symbol]["magmom"] for symbol in symbols]
543 # TODO: put more info into the atoms object
544 # e.g magmom, forces.
545 atoms = Atoms(symbols=symbols, positions=positions, cell=cell, pbc=True,
546 magmoms=magmoms)
547 atoms.set_constraint(convert_constraint_flags(constraint_flags))
549 return atoms
552def get_atomic_positions(lines, n_atoms, cell=None, alat=None):
553 """Parse atom positions from ATOMIC_POSITIONS card.
555 Parameters
556 ----------
557 lines : list[str]
558 A list of lines containing the ATOMIC_POSITIONS card.
559 n_atoms : int
560 Expected number of atoms. Only this many lines will be parsed.
561 cell : np.array
562 Unit cell of the crystal. Only used with crystal coordinates.
563 alat : float
564 Lattice parameter for atomic coordinates. Only used for alat case.
566 Returns
567 -------
568 positions : list[(str, (float, float, float), (int, int, int))]
569 A list of the ordered atomic positions in the format:
570 label, (x, y, z), (if_x, if_y, if_z)
571 Force multipliers are set to None if not present.
573 Raises
574 ------
575 ValueError
576 Any problems parsing the data result in ValueError
578 """
580 positions = None
581 # no blanks or comment lines, can the consume n_atoms lines for positions
582 trimmed_lines = (line for line in lines if line.strip() and line[0] != '#')
584 for line in trimmed_lines:
585 if line.strip().startswith('ATOMIC_POSITIONS'):
586 if positions is not None:
587 raise ValueError('Multiple ATOMIC_POSITIONS specified')
588 # Priority and behaviour tested with QE 5.3
589 if 'crystal_sg' in line.lower():
590 raise NotImplementedError('CRYSTAL_SG not implemented')
591 elif 'crystal' in line.lower():
592 cell = cell
593 elif 'bohr' in line.lower():
594 cell = np.identity(3) * units['Bohr']
595 elif 'angstrom' in line.lower():
596 cell = np.identity(3)
597 # elif 'alat' in line.lower():
598 # cell = np.identity(3) * alat
599 else:
600 if alat is None:
601 raise ValueError('Set lattice parameter in &SYSTEM for '
602 'alat coordinates')
603 # Always the default, will be DEPRECATED as mandatory
604 # in future
605 cell = np.identity(3) * alat
607 positions = []
608 for _ in range(n_atoms):
609 split_line = next(trimmed_lines).split()
610 # These can be fractions and other expressions
611 position = np.dot((infix_float(split_line[1]),
612 infix_float(split_line[2]),
613 infix_float(split_line[3])), cell)
614 if len(split_line) > 4:
615 force_mult = tuple(int(split_line[i]) for i in (4, 5, 6))
616 else:
617 force_mult = None
619 positions.append((split_line[0], position, force_mult))
621 return positions
624def get_atomic_species(lines, n_species):
625 """Parse atomic species from ATOMIC_SPECIES card.
627 Parameters
628 ----------
629 lines : list[str]
630 A list of lines containing the ATOMIC_POSITIONS card.
631 n_species : int
632 Expected number of atom types. Only this many lines will be parsed.
634 Returns
635 -------
636 species : list[(str, float, str)]
638 Raises
639 ------
640 ValueError
641 Any problems parsing the data result in ValueError
642 """
644 species = None
645 # no blanks or comment lines, can the consume n_atoms lines for positions
646 trimmed_lines = (line.strip() for line in lines
647 if line.strip() and not line.startswith('#'))
649 for line in trimmed_lines:
650 if line.startswith('ATOMIC_SPECIES'):
651 if species is not None:
652 raise ValueError('Multiple ATOMIC_SPECIES specified')
654 species = []
655 for _dummy in range(n_species):
656 label_weight_pseudo = next(trimmed_lines).split()
657 species.append((label_weight_pseudo[0],
658 float(label_weight_pseudo[1]),
659 label_weight_pseudo[2]))
661 return species
664def get_cell_parameters(lines, alat=None):
665 """Parse unit cell from CELL_PARAMETERS card.
667 Parameters
668 ----------
669 lines : list[str]
670 A list with lines containing the CELL_PARAMETERS card.
671 alat : float | None
672 Unit of lattice vectors in Angstrom. Only used if the card is
673 given in units of alat. alat must be None if CELL_PARAMETERS card
674 is in Bohr or Angstrom. For output files, alat will be parsed from
675 the card header and used in preference to this value.
677 Returns
678 -------
679 cell : np.array | None
680 Cell parameters as a 3x3 array in Angstrom. If no cell is found
681 None will be returned instead.
682 cell_alat : float | None
683 If a value for alat is given in the card header, this is also
684 returned, otherwise this will be None.
686 Raises
687 ------
688 ValueError
689 If CELL_PARAMETERS are given in units of bohr or angstrom
690 and alat is not
691 """
693 cell = None
694 cell_alat = None
695 # no blanks or comment lines, can take three lines for cell
696 trimmed_lines = (line for line in lines if line.strip() and line[0] != '#')
698 for line in trimmed_lines:
699 if line.strip().startswith('CELL_PARAMETERS'):
700 if cell is not None:
701 # multiple definitions
702 raise ValueError('CELL_PARAMETERS specified multiple times')
703 # Priority and behaviour tested with QE 5.3
704 if 'bohr' in line.lower():
705 if alat is not None:
706 raise ValueError('Lattice parameters given in '
707 '&SYSTEM celldm/A and CELL_PARAMETERS '
708 'bohr')
709 cell_units = units['Bohr']
710 elif 'angstrom' in line.lower():
711 if alat is not None:
712 raise ValueError('Lattice parameters given in '
713 '&SYSTEM celldm/A and CELL_PARAMETERS '
714 'angstrom')
715 cell_units = 1.0
716 elif 'alat' in line.lower():
717 # Output file has (alat = value) (in Bohrs)
718 if '=' in line:
719 alat = float(line.strip(') \n').split()[-1]) * units['Bohr']
720 cell_alat = alat
721 elif alat is None:
722 raise ValueError('Lattice parameters must be set in '
723 '&SYSTEM for alat units')
724 cell_units = alat
725 elif alat is None:
726 # may be DEPRECATED in future
727 cell_units = units['Bohr']
728 else:
729 # may be DEPRECATED in future
730 cell_units = alat
731 # Grab the parameters; blank lines have been removed
732 cell = [[ffloat(x) for x in next(trimmed_lines).split()[:3]],
733 [ffloat(x) for x in next(trimmed_lines).split()[:3]],
734 [ffloat(x) for x in next(trimmed_lines).split()[:3]]]
735 cell = np.array(cell) * cell_units
737 return cell, cell_alat
740def convert_constraint_flags(constraint_flags):
741 """Convert Quantum ESPRESSO constraint flags to ASE Constraint objects.
743 Parameters
744 ----------
745 constraint_flags : list[tuple[int, int, int]]
746 List of constraint flags (0: fixed, 1: moved) for all the atoms.
747 If the flag is None, there are no constraints on the atom.
749 Returns
750 -------
751 constraints : list[FixAtoms | FixCartesian]
752 List of ASE Constraint objects.
753 """
754 constraints = []
755 for i, constraint in enumerate(constraint_flags):
756 if constraint is None:
757 continue
758 # mask: False (0): moved, True (1): fixed
759 mask = ~np.asarray(constraint, bool)
760 constraints.append(FixCartesian(i, mask))
761 return canonicalize_constraints(constraints)
764def canonicalize_constraints(constraints):
765 """Canonicalize ASE FixCartesian constraints.
767 If the given FixCartesian constraints share the same `mask`, they can be
768 merged into one. Further, if `mask == (True, True, True)`, they can be
769 converted as `FixAtoms`. This method "canonicalizes" FixCartesian objects
770 in such a way.
772 Parameters
773 ----------
774 constraints : List[FixCartesian]
775 List of ASE FixCartesian constraints.
777 Returns
778 -------
779 constrants_canonicalized : List[FixAtoms | FixCartesian]
780 List of ASE Constraint objects.
781 """
782 # https://docs.python.org/3/library/collections.html#defaultdict-examples
783 indices_for_masks = defaultdict(list)
784 for constraint in constraints:
785 key = tuple((constraint.mask).tolist())
786 indices_for_masks[key].extend(constraint.index.tolist())
788 constraints_canonicalized = []
789 for mask, indices in indices_for_masks.items():
790 if mask == (False, False, False): # no directions are fixed
791 continue
792 if mask == (True, True, True): # all three directions are fixed
793 constraints_canonicalized.append(FixAtoms(indices))
794 else:
795 constraints_canonicalized.append(FixCartesian(indices, mask))
797 return constraints_canonicalized
800def str_to_value(string):
801 """Attempt to convert string into int, float (including fortran double),
802 or bool, in that order, otherwise return the string.
803 Valid (case-insensitive) bool values are: '.true.', '.t.', 'true'
804 and 't' (or false equivalents).
806 Parameters
807 ----------
808 string : str
809 Test to parse for a datatype
811 Returns
812 -------
813 value : any
814 Parsed string as the most appropriate datatype of int, float,
815 bool or string.
816 """
818 # Just an integer
819 try:
820 return int(string)
821 except ValueError:
822 pass
823 # Standard float
824 try:
825 return float(string)
826 except ValueError:
827 pass
828 # Fortran double
829 try:
830 return ffloat(string)
831 except ValueError:
832 pass
834 # possible bool, else just the raw string
835 if string.lower() in ('.true.', '.t.', 'true', 't'):
836 return True
837 elif string.lower() in ('.false.', '.f.', 'false', 'f'):
838 return False
839 else:
840 return string.strip("'")
843def read_fortran_namelist(fileobj):
844 """Takes a fortran-namelist formatted file and returns nested
845 dictionaries of sections and key-value data, followed by a list
846 of lines of text that do not fit the specifications.
847 Behaviour is taken from Quantum ESPRESSO 5.3. Parses fairly
848 convoluted files the same way that QE should, but may not get
849 all the MANDATORY rules and edge cases for very non-standard files
850 Ignores anything after '!' in a namelist, split pairs on ','
851 to include multiple key=values on a line, read values on section
852 start and end lines, section terminating character, '/', can appear
853 anywhere on a line. All of these are ignored if the value is in 'quotes'.
855 Parameters
856 ----------
857 fileobj : file
858 An open file-like object.
860 Returns
861 -------
862 data : dict[str, dict]
863 Dictionary for each section in the namelist with
864 key = value pairs of data.
865 additional_cards : list[str]
866 Any lines not used to create the data,
867 assumed to belong to 'cards' in the input file.
868 """
870 data = {}
871 card_lines = []
872 in_namelist = False
873 section = 'none' # can't be in a section without changing this
875 for line in fileobj:
876 # leading and trailing whitespace never needed
877 line = line.strip()
878 if line.startswith('&'):
879 # inside a namelist
880 section = line.split()[0][1:].lower() # case insensitive
881 if section in data:
882 # Repeated sections are completely ignored.
883 # (Note that repeated keys overwrite within a section)
884 section = "_ignored"
885 data[section] = {}
886 in_namelist = True
887 if not in_namelist and line:
888 # Stripped line is Truthy, so safe to index first character
889 if line[0] not in ('!', '#'):
890 card_lines.append(line)
891 if in_namelist:
892 # parse k, v from line:
893 key = []
894 value = None
895 in_quotes = False
896 for character in line:
897 if character == ',' and value is not None and not in_quotes:
898 # finished value:
899 data[section][''.join(key).strip()] = str_to_value(
900 ''.join(value).strip())
901 key = []
902 value = None
903 elif character == '=' and value is None and not in_quotes:
904 # start writing value
905 value = []
906 elif character == "'":
907 # only found in value anyway
908 in_quotes = not in_quotes
909 value.append("'")
910 elif character == '!' and not in_quotes:
911 break
912 elif character == '/' and not in_quotes:
913 in_namelist = False
914 break
915 elif value is not None:
916 value.append(character)
917 else:
918 key.append(character)
919 if value is not None:
920 data[section][''.join(key).strip()] = str_to_value(
921 ''.join(value).strip())
923 return Namelist(data), card_lines
926def ffloat(string):
927 """Parse float from fortran compatible float definitions.
929 In fortran exponents can be defined with 'd' or 'q' to symbolise
930 double or quad precision numbers. Double precision numbers are
931 converted to python floats and quad precision values are interpreted
932 as numpy longdouble values (platform specific precision).
934 Parameters
935 ----------
936 string : str
937 A string containing a number in fortran real format
939 Returns
940 -------
941 value : float | np.longdouble
942 Parsed value of the string.
944 Raises
945 ------
946 ValueError
947 Unable to parse a float value.
949 """
951 if 'q' in string.lower():
952 return np.longdouble(string.lower().replace('q', 'e'))
953 else:
954 return float(string.lower().replace('d', 'e'))
957def label_to_symbol(label):
958 """Convert a valid espresso ATOMIC_SPECIES label to a
959 chemical symbol.
961 Parameters
962 ----------
963 label : str
964 chemical symbol X (1 or 2 characters, case-insensitive)
965 or chemical symbol plus a number or a letter, as in
966 "Xn" (e.g. Fe1) or "X_*" or "X-*" (e.g. C1, C_h;
967 max total length cannot exceed 3 characters).
969 Returns
970 -------
971 symbol : str
972 The best matching species from ase.utils.chemical_symbols
974 Raises
975 ------
976 KeyError
977 Couldn't find an appropriate species.
979 Notes
980 -----
981 It's impossible to tell whether e.g. He is helium
982 or hydrogen labelled 'e'.
983 """
985 # possibly a two character species
986 # ase Atoms need proper case of chemical symbols.
987 if len(label) >= 2:
988 test_symbol = label[0].upper() + label[1].lower()
989 if test_symbol in chemical_symbols:
990 return test_symbol
991 # finally try with one character
992 test_symbol = label[0].upper()
993 if test_symbol in chemical_symbols:
994 return test_symbol
995 else:
996 raise KeyError('Could not parse species from label {}.'
997 ''.format(label))
1000def infix_float(text):
1001 """Parse simple infix maths into a float for compatibility with
1002 Quantum ESPRESSO ATOMIC_POSITIONS cards. Note: this works with the
1003 example, and most simple expressions, but the capabilities of
1004 the two parsers are not identical. Will also parse a normal float
1005 value properly, but slowly.
1007 >>> infix_float('1/2*3^(-1/2)')
1008 0.28867513459481287
1010 Parameters
1011 ----------
1012 text : str
1013 An arithmetic expression using +, -, *, / and ^, including brackets.
1015 Returns
1016 -------
1017 value : float
1018 Result of the mathematical expression.
1020 """
1022 def middle_brackets(full_text):
1023 """Extract text from innermost brackets."""
1024 start, end = 0, len(full_text)
1025 for (idx, char) in enumerate(full_text):
1026 if char == '(':
1027 start = idx
1028 if char == ')':
1029 end = idx + 1
1030 break
1031 return full_text[start:end]
1033 def eval_no_bracket_expr(full_text):
1034 """Calculate value of a mathematical expression, no brackets."""
1035 exprs = [('+', op.add), ('*', op.mul),
1036 ('/', op.truediv), ('^', op.pow)]
1037 full_text = full_text.lstrip('(').rstrip(')')
1038 try:
1039 return float(full_text)
1040 except ValueError:
1041 for symbol, func in exprs:
1042 if symbol in full_text:
1043 left, right = full_text.split(symbol, 1) # single split
1044 return func(eval_no_bracket_expr(left),
1045 eval_no_bracket_expr(right))
1047 while '(' in text:
1048 middle = middle_brackets(text)
1049 text = text.replace(middle, f'{eval_no_bracket_expr(middle)}')
1051 return float(eval_no_bracket_expr(text))
1054# Number of valence electrons in the pseudopotentials recommended by
1055# http://materialscloud.org/sssp/. These are just used as a fallback for
1056# calculating initial magetization values which are given as a fraction
1057# of valence electrons.
1058SSSP_VALENCE = [
1059 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,
1060 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,
1061 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,
1062 13.0, 14.0, 15.0, 16.0, 17.0, 18.0, 19.0, 12.0, 13.0, 14.0, 15.0, 6.0,
1063 7.0, 18.0, 9.0, 10.0, 11.0, 12.0, 13.0, 14.0, 15.0, 16.0, 17.0, 18.0,
1064 19.0, 20.0, 21.0, 22.0, 23.0, 24.0, 25.0, 36.0, 27.0, 14.0, 15.0, 30.0,
1065 15.0, 32.0, 19.0, 12.0, 13.0, 14.0, 15.0, 16.0, 18.0]
1068def kspacing_to_grid(atoms, spacing, calculated_spacing=None):
1069 """
1070 Calculate the kpoint mesh that is equivalent to the given spacing
1071 in reciprocal space (units Angstrom^-1). The number of kpoints is each
1072 dimension is rounded up (compatible with CASTEP).
1074 Parameters
1075 ----------
1076 atoms: ase.Atoms
1077 A structure that can have get_reciprocal_cell called on it.
1078 spacing: float
1079 Minimum K-Point spacing in $A^{-1}$.
1080 calculated_spacing : list
1081 If a three item list (or similar mutable sequence) is given the
1082 members will be replaced with the actual calculated spacing in
1083 $A^{-1}$.
1085 Returns
1086 -------
1087 kpoint_grid : [int, int, int]
1088 MP grid specification to give the required spacing.
1090 """
1091 # No factor of 2pi in ase, everything in A^-1
1092 # reciprocal dimensions
1093 r_x, r_y, r_z = np.linalg.norm(atoms.cell.reciprocal(), axis=1)
1095 kpoint_grid = [int(r_x / spacing) + 1,
1096 int(r_y / spacing) + 1,
1097 int(r_z / spacing) + 1]
1099 for i, _ in enumerate(kpoint_grid):
1100 if not atoms.pbc[i]:
1101 kpoint_grid[i] = 1
1103 if calculated_spacing is not None:
1104 calculated_spacing[:] = [r_x / kpoint_grid[0],
1105 r_y / kpoint_grid[1],
1106 r_z / kpoint_grid[2]]
1108 return kpoint_grid
1111def format_atom_position(atom, crystal_coordinates, mask='', tidx=None):
1112 """Format one line of atomic positions in
1113 Quantum ESPRESSO ATOMIC_POSITIONS card.
1115 >>> for atom in make_supercell(bulk('Li', 'bcc'), np.ones(3)-np.eye(3)):
1116 >>> format_atom_position(atom, True)
1117 Li 0.0000000000 0.0000000000 0.0000000000
1118 Li 0.5000000000 0.5000000000 0.5000000000
1120 Parameters
1121 ----------
1122 atom : Atom
1123 A structure that has symbol and [position | (a, b, c)].
1124 crystal_coordinates: bool
1125 Whether the atomic positions should be written to the QE input file in
1126 absolute (False, default) or relative (crystal) coordinates (True).
1127 mask, optional : str
1128 String of ndim=3 0 or 1 for constraining atomic positions.
1129 tidx, optional : int
1130 Magnetic type index.
1132 Returns
1133 -------
1134 atom_line : str
1135 Input line for atom position
1136 """
1137 if crystal_coordinates:
1138 coords = [atom.a, atom.b, atom.c]
1139 else:
1140 coords = atom.position
1141 line_fmt = '{atom.symbol}'
1142 inps = dict(atom=atom)
1143 if tidx is not None:
1144 line_fmt += '{tidx}'
1145 inps["tidx"] = tidx
1146 line_fmt += ' {coords[0]:.10f} {coords[1]:.10f} {coords[2]:.10f} '
1147 inps["coords"] = coords
1148 line_fmt += ' ' + mask + '\n'
1149 astr = line_fmt.format(**inps)
1150 return astr
1153@writer
1154def write_espresso_in(fd, atoms, input_data=None, pseudopotentials=None,
1155 kspacing=None, kpts=None, koffset=(0, 0, 0),
1156 crystal_coordinates=False, additional_cards=None,
1157 **kwargs):
1158 """
1159 Create an input file for pw.x.
1161 Use set_initial_magnetic_moments to turn on spin, if nspin is set to 2
1162 with no magnetic moments, they will all be set to 0.0. Magnetic moments
1163 will be converted to the QE units (fraction of valence electrons) using
1164 any pseudopotential files found, or a best guess for the number of
1165 valence electrons.
1167 Units are not converted for any other input data, so use Quantum ESPRESSO
1168 units (Usually Ry or atomic units).
1170 Keys with a dimension (e.g. Hubbard_U(1)) will be incorporated as-is
1171 so the `i` should be made to match the output.
1173 Implemented features:
1175 - Conversion of :class:`ase.constraints.FixAtoms` and
1176 :class:`ase.constraints.FixCartesian`.
1177 - ``starting_magnetization`` derived from the ``magmoms`` and
1178 pseudopotentials (searches default paths for pseudo files.)
1179 - Automatic assignment of options to their correct sections.
1181 Not implemented:
1183 - Non-zero values of ibrav
1184 - Lists of k-points
1185 - Other constraints
1186 - Hubbard parameters
1187 - Validation of the argument types for input
1188 - Validation of required options
1190 Parameters
1191 ----------
1192 fd: file | str
1193 A file to which the input is written.
1194 atoms: Atoms
1195 A single atomistic configuration to write to ``fd``.
1196 input_data: dict
1197 A flat or nested dictionary with input parameters for pw.x
1198 pseudopotentials: dict
1199 A filename for each atomic species, e.g.
1200 {'O': 'O.pbe-rrkjus.UPF', 'H': 'H.pbe-rrkjus.UPF'}.
1201 A dummy name will be used if none are given.
1202 kspacing: float
1203 Generate a grid of k-points with this as the minimum distance,
1204 in A^-1 between them in reciprocal space. If set to None, kpts
1205 will be used instead.
1206 kpts: (int, int, int), dict or np.ndarray
1207 If ``kpts`` is a tuple (or list) of 3 integers, it is interpreted
1208 as the dimensions of a Monkhorst-Pack grid.
1209 If ``kpts`` is set to ``None``, only the Γ-point will be included
1210 and QE will use routines optimized for Γ-point-only calculations.
1211 Compared to Γ-point-only calculations without this optimization
1212 (i.e. with ``kpts=(1, 1, 1)``), the memory and CPU requirements
1213 are typically reduced by half.
1214 If kpts is a dict, it will either be interpreted as a path
1215 in the Brillouin zone (*) if it contains the 'path' keyword,
1216 otherwise it is converted to a Monkhorst-Pack grid (**).
1217 If ``kpts`` is a NumPy array, the raw k-points will be passed to
1218 Quantum Espresso as given in the array (in crystal coordinates).
1219 Must be of shape (n_kpts, 4). The fourth column contains the
1220 k-point weights.
1221 (*) see ase.dft.kpoints.bandpath
1222 (**) see ase.calculators.calculator.kpts2sizeandoffsets
1223 koffset: (int, int, int)
1224 Offset of kpoints in each direction. Must be 0 (no offset) or
1225 1 (half grid offset). Setting to True is equivalent to (1, 1, 1).
1226 crystal_coordinates: bool
1227 Whether the atomic positions should be written to the QE input file in
1228 absolute (False, default) or relative (crystal) coordinates (True).
1230 """
1232 # Convert to a namelist to make working with parameters much easier
1233 # Note that the name ``input_data`` is chosen to prevent clash with
1234 # ``parameters`` in Calculator objects
1235 input_parameters = Namelist(input_data)
1236 input_parameters.to_nested('pw', **kwargs)
1238 # Convert ase constraints to QE constraints
1239 # Nx3 array of force multipliers matches what QE uses
1240 # Do this early so it is available when constructing the atoms card
1241 moved = np.ones((len(atoms), 3), dtype=bool)
1242 for constraint in atoms.constraints:
1243 if isinstance(constraint, FixAtoms):
1244 moved[constraint.index] = False
1245 elif isinstance(constraint, FixCartesian):
1246 moved[constraint.index] = ~constraint.mask
1247 else:
1248 warnings.warn(f'Ignored unknown constraint {constraint}')
1249 masks = []
1250 for atom in atoms:
1251 # only inclued mask if something is fixed
1252 if not all(moved[atom.index]):
1253 mask = ' {:d} {:d} {:d}'.format(*moved[atom.index])
1254 else:
1255 mask = ''
1256 masks.append(mask)
1258 # Species info holds the information on the pseudopotential and
1259 # associated for each element
1260 if pseudopotentials is None:
1261 pseudopotentials = {}
1262 species_info = {}
1263 for species in set(atoms.get_chemical_symbols()):
1264 # Look in all possible locations for the pseudos and try to figure
1265 # out the number of valence electrons
1266 pseudo = pseudopotentials[species]
1267 species_info[species] = {'pseudo': pseudo}
1269 # Convert atoms into species.
1270 # Each different magnetic moment needs to be a separate type even with
1271 # the same pseudopotential (e.g. an up and a down for AFM).
1272 # if any magmom are > 0 or nspin == 2 then use species labels.
1273 # Rememeber: magnetisation uses 1 based indexes
1274 atomic_species = {}
1275 atomic_species_str = []
1276 atomic_positions_str = []
1278 nspin = input_parameters['system'].get('nspin', 1) # 1 is the default
1279 noncolin = input_parameters['system'].get('noncolin', False)
1280 rescale_magmom_fac = kwargs.get('rescale_magmom_fac', 1.0)
1281 if any(atoms.get_initial_magnetic_moments()):
1282 if nspin == 1 and not noncolin:
1283 # Force spin on
1284 input_parameters['system']['nspin'] = 2
1285 nspin = 2
1287 if nspin == 2 or noncolin:
1288 # Magnetic calculation on
1289 for atom, mask, magmom in zip(
1290 atoms, masks, atoms.get_initial_magnetic_moments()):
1291 if (atom.symbol, magmom) not in atomic_species:
1292 # for qe version 7.2 or older magmon must be rescale by
1293 # about a factor 10 to assume sensible values
1294 # since qe-v7.3 magmom values will be provided unscaled
1295 fspin = float(magmom) / rescale_magmom_fac
1296 # Index in the atomic species list
1297 sidx = len(atomic_species) + 1
1298 # Index for that atom type; no index for first one
1299 tidx = sum(atom.symbol == x[0] for x in atomic_species) or ' '
1300 atomic_species[(atom.symbol, magmom)] = (sidx, tidx)
1301 # Add magnetization to the input file
1302 mag_str = f"starting_magnetization({sidx})"
1303 input_parameters['system'][mag_str] = fspin
1304 species_pseudo = species_info[atom.symbol]['pseudo']
1305 atomic_species_str.append(
1306 f"{atom.symbol}{tidx} {atom.mass} {species_pseudo}\n")
1307 # lookup tidx to append to name
1308 sidx, tidx = atomic_species[(atom.symbol, magmom)]
1309 # construct line for atomic positions
1310 atomic_positions_str.append(
1311 format_atom_position(
1312 atom, crystal_coordinates, mask=mask, tidx=tidx)
1313 )
1314 else:
1315 # Do nothing about magnetisation
1316 for atom, mask in zip(atoms, masks):
1317 if atom.symbol not in atomic_species:
1318 atomic_species[atom.symbol] = True # just a placeholder
1319 species_pseudo = species_info[atom.symbol]['pseudo']
1320 atomic_species_str.append(
1321 f"{atom.symbol} {atom.mass} {species_pseudo}\n")
1322 # construct line for atomic positions
1323 atomic_positions_str.append(
1324 format_atom_position(atom, crystal_coordinates, mask=mask)
1325 )
1327 # Add computed parameters
1328 # different magnetisms means different types
1329 input_parameters['system']['ntyp'] = len(atomic_species)
1330 input_parameters['system']['nat'] = len(atoms)
1332 # Use cell as given or fit to a specific ibrav
1333 if 'ibrav' in input_parameters['system']:
1334 ibrav = input_parameters['system']['ibrav']
1335 if ibrav != 0:
1336 raise ValueError(ibrav_error_message)
1337 else:
1338 # Just use standard cell block
1339 input_parameters['system']['ibrav'] = 0
1341 # Construct input file into this
1342 pwi = input_parameters.to_string(list_form=True)
1344 # Pseudopotentials
1345 pwi.append('ATOMIC_SPECIES\n')
1346 pwi.extend(atomic_species_str)
1347 pwi.append('\n')
1349 # KPOINTS - add a MP grid as required
1350 if kspacing is not None:
1351 kgrid = kspacing_to_grid(atoms, kspacing)
1352 elif kpts is not None:
1353 if isinstance(kpts, dict) and 'path' not in kpts:
1354 kgrid, shift = kpts2sizeandoffsets(atoms=atoms, **kpts)
1355 koffset = []
1356 for i, x in enumerate(shift):
1357 assert x == 0 or abs(x * kgrid[i] - 0.5) < 1e-14
1358 koffset.append(0 if x == 0 else 1)
1359 else:
1360 kgrid = kpts
1361 else:
1362 kgrid = "gamma"
1364 # True and False work here and will get converted by ':d' format
1365 if isinstance(koffset, int):
1366 koffset = (koffset, ) * 3
1368 # BandPath object or bandpath-as-dictionary:
1369 if isinstance(kgrid, dict) or hasattr(kgrid, 'kpts'):
1370 pwi.append('K_POINTS crystal_b\n')
1371 assert hasattr(kgrid, 'path') or 'path' in kgrid
1372 kgrid = kpts2ndarray(kgrid, atoms=atoms)
1373 pwi.append(f'{len(kgrid)}\n')
1374 for k in kgrid:
1375 pwi.append(f"{k[0]:.14f} {k[1]:.14f} {k[2]:.14f} 0\n")
1376 pwi.append('\n')
1377 elif isinstance(kgrid, str) and (kgrid == "gamma"):
1378 pwi.append('K_POINTS gamma\n')
1379 pwi.append('\n')
1380 elif isinstance(kgrid, np.ndarray):
1381 if np.shape(kgrid)[1] != 4:
1382 raise ValueError('Only Nx4 kgrids are supported right now.')
1383 pwi.append('K_POINTS crystal\n')
1384 pwi.append(f'{len(kgrid)}\n')
1385 for k in kgrid:
1386 pwi.append(f"{k[0]:.14f} {k[1]:.14f} {k[2]:.14f} {k[3]:.14f}\n")
1387 pwi.append('\n')
1388 else:
1389 pwi.append('K_POINTS automatic\n')
1390 pwi.append(f"{kgrid[0]} {kgrid[1]} {kgrid[2]} "
1391 f" {koffset[0]:d} {koffset[1]:d} {koffset[2]:d}\n")
1392 pwi.append('\n')
1394 # CELL block, if required
1395 if input_parameters['SYSTEM']['ibrav'] == 0:
1396 pwi.append('CELL_PARAMETERS angstrom\n')
1397 pwi.append('{cell[0][0]:.14f} {cell[0][1]:.14f} {cell[0][2]:.14f}\n'
1398 '{cell[1][0]:.14f} {cell[1][1]:.14f} {cell[1][2]:.14f}\n'
1399 '{cell[2][0]:.14f} {cell[2][1]:.14f} {cell[2][2]:.14f}\n'
1400 ''.format(cell=atoms.cell))
1401 pwi.append('\n')
1403 # Positions - already constructed, but must appear after namelist
1404 if crystal_coordinates:
1405 pwi.append('ATOMIC_POSITIONS crystal\n')
1406 else:
1407 pwi.append('ATOMIC_POSITIONS angstrom\n')
1408 pwi.extend(atomic_positions_str)
1409 pwi.append('\n')
1411 # DONE!
1412 fd.write(''.join(pwi))
1414 if additional_cards:
1415 if isinstance(additional_cards, list):
1416 additional_cards = "\n".join(additional_cards)
1417 additional_cards += "\n"
1419 fd.write(additional_cards)
1422@writer
1423def write_espresso_ph(
1424 fd,
1425 input_data=None,
1426 qpts=None,
1427 nat_todo_indices=None,
1428 **kwargs) -> None:
1429 """
1430 Function that write the input file for a ph.x calculation. Normal namelist
1431 cards are passed in the input_data dictionary. Which can be either nested
1432 or flat, ASE style. The q-points are passed in the qpts list. If qplot is
1433 set to True then qpts is expected to be a list of list|tuple of length 4.
1434 Where the first three elements are the coordinates of the q-point in units
1435 of 2pi/alat and the last element is the weight of the q-point. if qplot is
1436 set to False then qpts is expected to be a simple list of length 4 (single
1437 q-point). Finally if ldisp is set to True, the above is discarded and the
1438 q-points are read from the nq1, nq2, nq3 cards in the input_data dictionary.
1440 Additionally, a nat_todo_indices kwargs (list[int]) can be specified in the
1441 kwargs. It will be used if nat_todo is set to True in the input_data
1442 dictionary.
1444 Globally, this function follows the convention set in the ph.x documentation
1445 (https://www.quantum-espresso.org/Doc/INPUT_PH.html)
1447 Parameters
1448 ----------
1449 fd
1450 The file descriptor of the input file.
1452 kwargs
1453 kwargs dictionary possibly containing the following keys:
1455 - input_data: dict
1456 - qpts: list[list[float]] | list[tuple[float]] | list[float]
1457 - nat_todo_indices: list[int]
1459 Returns
1460 -------
1461 None
1462 """
1464 input_data = Namelist(input_data)
1465 input_data.to_nested('ph', **kwargs)
1467 input_ph = input_data["inputph"]
1469 inp_nat_todo = input_ph.get("nat_todo", 0)
1470 qpts = qpts or (0, 0, 0)
1472 pwi = input_data.to_string()
1474 fd.write(pwi)
1476 qplot = input_ph.get("qplot", False)
1477 ldisp = input_ph.get("ldisp", False)
1479 if qplot:
1480 fd.write(f"{len(qpts)}\n")
1481 for qpt in qpts:
1482 fd.write(
1483 f"{qpt[0]:0.8f} {qpt[1]:0.8f} {qpt[2]:0.8f} {qpt[3]:1d}\n"
1484 )
1485 elif not (qplot or ldisp):
1486 fd.write(f"{qpts[0]:0.8f} {qpts[1]:0.8f} {qpts[2]:0.8f}\n")
1487 if inp_nat_todo:
1488 tmp = [str(i) for i in nat_todo_indices]
1489 fd.write(" ".join(tmp))
1490 fd.write("\n")
1493class _PHHelper:
1494 freg = re.compile(r"-?(?:0|[1-9]\d*)(?:\.\d+)?(?:[eE][+\-]?\d+)?")
1496 def __init__(self, fdo_lines):
1497 self.fdo_lines = fdo_lines
1498 self.n_lines = len(fdo_lines)
1500 def _read_qpoints(self, idx):
1501 match = re.findall(self.freg, self.fdo_lines[idx])
1502 return tuple(round(float(x), 7) for x in match)
1504 def _read_kpoints(self, idx):
1505 n_kpts = int(re.findall(self.freg, self.fdo_lines[idx])[0])
1506 kpts = []
1507 for line in self.fdo_lines[idx + 2: idx + 2 + n_kpts]:
1508 if bool(re.search(r"^\s*k\(.*wk", line)):
1509 kpts.append([round(float(x), 7)
1510 for x in re.findall(self.freg, line)[1:]])
1511 return np.array(kpts)
1513 def _read_repr(self, idx):
1514 n_repr, curr, n = int(re.findall(self.freg,
1515 self.fdo_lines[idx])[0]), 0, 0
1516 representations = {}
1517 while idx + n < self.n_lines:
1518 if re.search(r"^\s*Representation.*modes", self.fdo_lines[idx + n]):
1519 curr = int(re.findall(self.freg, self.fdo_lines[idx + n])[0])
1520 representations[curr] = {"done": False, "modes": []}
1521 if re.search(r"Calculated\s*using\s*symmetry",
1522 self.fdo_lines[idx + n]) \
1523 or re.search(r"-\s*Done\s*$", self.fdo_lines[idx + n]):
1524 representations[curr]["done"] = True
1525 if re.search(r"(?i)^\s*(mode\s*#\s*\d\s*)+",
1526 self.fdo_lines[idx + n]):
1527 representations[curr]["modes"] = self._read_modes(idx + n)
1528 if curr == n_repr:
1529 break
1530 n += 1
1531 return representations
1533 def _read_modes(self, idx):
1534 n = 1
1535 n_modes = len(re.findall(r"mode", self.fdo_lines[idx]))
1536 modes = []
1537 while not modes or bool(re.match(r"^\s*\(", self.fdo_lines[idx + n])):
1538 tmp = re.findall(self.freg, self.fdo_lines[idx + n])
1539 modes.append([round(float(x), 5) for x in tmp])
1540 n += 1
1541 return np.hsplit(np.array(modes), n_modes)
1543 def _read_eqpoints(self, idx):
1544 n_star = int(re.findall(self.freg, self.fdo_lines[idx])[0])
1545 return np.loadtxt(
1546 self.fdo_lines[idx + 2: idx + 2 + n_star], usecols=(1, 2, 3)
1547 ).reshape(-1, 3)
1549 def _read_freqs(self, idx):
1550 n = 0
1551 freqs = []
1552 stop = 0
1553 while not freqs or stop < 2:
1554 if bool(re.search(r"^\s*freq", self.fdo_lines[idx + n])):
1555 tmp = re.findall(self.freg, self.fdo_lines[idx + n])[1]
1556 freqs.append(float(tmp))
1557 if bool(re.search(r"\*{5,}", self.fdo_lines[idx + n])):
1558 stop += 1
1559 n += 1
1560 return np.array(freqs)
1562 def _read_sym(self, idx):
1563 n = 1
1564 sym = {}
1565 while bool(re.match(r"^\s*freq", self.fdo_lines[idx + n])):
1566 r = re.findall("\\d+", self.fdo_lines[idx + n])
1567 r = tuple(range(int(r[0]), int(r[1]) + 1))
1568 sym[r] = self.fdo_lines[idx + n].split("-->")[1].strip()
1569 sym[r] = re.sub(r"\s+", " ", sym[r])
1570 n += 1
1571 return sym
1573 def _read_epsil(self, idx):
1574 epsil = np.zeros((3, 3))
1575 for n in range(1, 4):
1576 tmp = re.findall(self.freg, self.fdo_lines[idx + n])
1577 epsil[n - 1] = [round(float(x), 9) for x in tmp]
1578 return epsil
1580 def _read_born(self, idx):
1581 n = 1
1582 born = []
1583 while idx + n < self.n_lines:
1584 if re.search(r"^\s*atom\s*\d\s*\S", self.fdo_lines[idx + n]):
1585 pass
1586 elif re.search(r"^\s*E\*?(x|y|z)\s*\(", self.fdo_lines[idx + n]):
1587 tmp = re.findall(self.freg, self.fdo_lines[idx + n])
1588 born.append([round(float(x), 5) for x in tmp])
1589 else:
1590 break
1591 n += 1
1592 born = np.array(born)
1593 return np.vsplit(born, len(born) // 3)
1595 def _read_born_dfpt(self, idx):
1596 n = 1
1597 born = []
1598 while idx + n < self.n_lines:
1599 if re.search(r"^\s*atom\s*\d\s*\S", self.fdo_lines[idx + n]):
1600 pass
1601 elif re.search(r"^\s*P(x|y|z)\s*\(", self.fdo_lines[idx + n]):
1602 tmp = re.findall(self.freg, self.fdo_lines[idx + n])
1603 born.append([round(float(x), 5) for x in tmp])
1604 else:
1605 break
1606 n += 1
1607 born = np.array(born)
1608 return np.vsplit(born, len(born) // 3)
1610 def _read_pola(self, idx):
1611 pola = np.zeros((3, 3))
1612 for n in range(1, 4):
1613 tmp = re.findall(self.freg, self.fdo_lines[idx + n])[:3]
1614 pola[n - 1] = [round(float(x), 2) for x in tmp]
1615 return pola
1617 def _read_positions(self, idx):
1618 positions = []
1619 symbols = []
1620 n = 1
1621 while re.findall(r"^\s*\d+", self.fdo_lines[idx + n]):
1622 symbols.append(self.fdo_lines[idx + n].split()[1])
1623 positions.append(
1624 [round(float(x), 5)
1625 for x in re.findall(self.freg, self.fdo_lines[idx + n])[-3:]]
1626 )
1627 n += 1
1628 atoms = Atoms(positions=positions, symbols=symbols)
1629 atoms.pbc = True
1630 return atoms
1632 def _read_alat(self, idx):
1633 return round(float(re.findall(self.freg, self.fdo_lines[idx])[1]), 5)
1635 def _read_cell(self, idx):
1636 cell = []
1637 n = 1
1638 while re.findall(r"^\s*a\(\d\)", self.fdo_lines[idx + n]):
1639 cell.append(
1640 [round(float(x), 4)
1641 for x in re.findall(self.freg, self.fdo_lines[idx + n])[-3:]])
1642 n += 1
1643 return np.array(cell)
1645 def _read_electron_phonon(self, idx):
1646 results = {}
1648 broad_re = (
1649 r"^\s*Gaussian\s*Broadening:\s+([\d.]+)\s+Ry, ngauss=\s+\d+"
1650 )
1651 dos_re = (
1652 r"^\s*DOS\s*=\s*([\d.]+)\s*states/"
1653 r"spin/Ry/Unit\s*Cell\s*at\s*Ef=\s+([\d.]+)\s+eV"
1654 )
1655 lg_re = (
1656 r"^\s*lambda\(\s+(\d+)\)=\s+([\d.]+)\s+gamma=\s+([\d.]+)\s+GHz"
1657 )
1658 end_re = r"^\s*Number\s*of\s*q\s*in\s*the\s*star\s*=\s+(\d+)$"
1660 lambdas = []
1661 gammas = []
1663 current = None
1665 n = 1
1666 while idx + n < self.n_lines:
1667 line = self.fdo_lines[idx + n]
1669 broad_match = re.match(broad_re, line)
1670 dos_match = re.match(dos_re, line)
1671 lg_match = re.match(lg_re, line)
1672 end_match = re.match(end_re, line)
1674 if broad_match:
1675 if lambdas:
1676 results[current]["lambdas"] = lambdas
1677 results[current]["gammas"] = gammas
1678 lambdas = []
1679 gammas = []
1680 current = float(broad_match[1])
1681 results[current] = {}
1682 elif dos_match:
1683 results[current]["dos"] = float(dos_match[1])
1684 results[current]["fermi"] = float(dos_match[2])
1685 elif lg_match:
1686 lambdas.append(float(lg_match[2]))
1687 gammas.append(float(lg_match[3]))
1689 if end_match:
1690 results[current]["lambdas"] = lambdas
1691 results[current]["gammas"] = gammas
1692 break
1694 n += 1
1696 return results
1699@reader
1700def read_espresso_ph(fileobj):
1701 """
1702 Function that reads the output of a ph.x calculation.
1703 It returns a dictionary where each q-point number is a key and
1704 the value is a dictionary with the following keys if available:
1706 - qpoints: The q-point in cartesian coordinates.
1707 - kpoints: The k-points in cartesian coordinates.
1708 - dieltensor: The dielectric tensor.
1709 - borneffcharge: The effective Born charges.
1710 - borneffcharge_dfpt: The effective Born charges from DFPT.
1711 - polarizability: The polarizability tensor.
1712 - modes: The phonon modes.
1713 - eqpoints: The symmetrically equivalent q-points.
1714 - freqs: The phonon frequencies.
1715 - mode_symmetries: The symmetries of the modes.
1716 - atoms: The atoms object.
1718 Some notes:
1720 - For some reason, the cell is not defined to high level of
1721 precision in ph.x outputs. Be careful when using the atoms object
1722 retrieved from this function.
1723 - This function can be called on incomplete calculations i.e.
1724 if the calculation couldn't diagonalize the dynamical matrix
1725 for some q-points, the results for the other q-points will
1726 still be returned.
1728 Parameters
1729 ----------
1730 fileobj
1731 The file descriptor of the output file.
1733 Returns
1734 -------
1735 dict
1736 The results dictionnary as described above.
1737 """
1739 QPOINTS = r"(?i)^\s*Calculation\s*of\s*q"
1740 NKPTS = r"(?i)^\s*number\s*of\s*k\s*points\s*"
1741 DIEL = r"(?i)^\s*Dielectric\s*constant\s*in\s*cartesian\s*axis\s*$"
1742 BORN = r"(?i)^\s*Effective\s*charges\s*\(d\s*Force\s*/\s*dE\)"
1743 POLA = r"(?i)^\s*Polarizability\s*(a.u.)\^3"
1744 REPR = r"(?i)^\s*There\s*are\s*\d+\s*irreducible\s*representations\s*$"
1745 EQPOINTS = r"(?i)^\s*Number\s*of\s*q\s*in\s*the\s*star\s*=\s*"
1746 DIAG = r"(?i)^\s*Diagonalizing\s*the\s*dynamical\s*matrix\s*$"
1747 MODE_SYM = r"(?i)^\s*Mode\s*symmetry,\s*"
1748 BORN_DFPT = r"(?i)^\s*Effective\s*charges\s*\(d\s*P\s*/\s*du\)"
1749 POSITIONS = r"(?i)^\s*site\s*n\..*\(alat\s*units\)"
1750 ALAT = r"(?i)^\s*celldm\(1\)="
1751 CELL = (
1752 r"^\s*crystal\s*axes:\s*\(cart.\s*coord.\s*in\s*units\s*of\s*alat\)"
1753 )
1754 ELECTRON_PHONON = r"(?i)^\s*electron-phonon\s*interaction\s*...\s*$"
1756 output = {
1757 QPOINTS: [],
1758 NKPTS: [],
1759 DIEL: [],
1760 BORN: [],
1761 BORN_DFPT: [],
1762 POLA: [],
1763 REPR: [],
1764 EQPOINTS: [],
1765 DIAG: [],
1766 MODE_SYM: [],
1767 POSITIONS: [],
1768 ALAT: [],
1769 CELL: [],
1770 ELECTRON_PHONON: [],
1771 }
1773 names = {
1774 QPOINTS: "qpoints",
1775 NKPTS: "kpoints",
1776 DIEL: "dieltensor",
1777 BORN: "borneffcharge",
1778 BORN_DFPT: "borneffcharge_dfpt",
1779 POLA: "polarizability",
1780 REPR: "representations",
1781 EQPOINTS: "eqpoints",
1782 DIAG: "freqs",
1783 MODE_SYM: "mode_symmetries",
1784 POSITIONS: "positions",
1785 ALAT: "alat",
1786 CELL: "cell",
1787 ELECTRON_PHONON: "ep_data",
1788 }
1790 unique = {
1791 QPOINTS: True,
1792 NKPTS: False,
1793 DIEL: True,
1794 BORN: True,
1795 BORN_DFPT: True,
1796 POLA: True,
1797 REPR: True,
1798 EQPOINTS: True,
1799 DIAG: True,
1800 MODE_SYM: True,
1801 POSITIONS: True,
1802 ALAT: True,
1803 CELL: True,
1804 ELECTRON_PHONON: True,
1805 }
1807 results = {}
1808 fdo_lines = [i for i in fileobj.read().splitlines() if i]
1809 n_lines = len(fdo_lines)
1811 for idx, line in enumerate(fdo_lines):
1812 for key in output:
1813 if bool(re.match(key, line)):
1814 output[key].append(idx)
1816 output = {key: np.array(value) for key, value in output.items()}
1818 helper = _PHHelper(fdo_lines)
1820 properties = {
1821 NKPTS: helper._read_kpoints,
1822 DIEL: helper._read_epsil,
1823 BORN: helper._read_born,
1824 BORN_DFPT: helper._read_born_dfpt,
1825 POLA: helper._read_pola,
1826 REPR: helper._read_repr,
1827 EQPOINTS: helper._read_eqpoints,
1828 DIAG: helper._read_freqs,
1829 MODE_SYM: helper._read_sym,
1830 POSITIONS: helper._read_positions,
1831 ALAT: helper._read_alat,
1832 CELL: helper._read_cell,
1833 ELECTRON_PHONON: helper._read_electron_phonon,
1834 }
1836 iblocks = np.append(output[QPOINTS], n_lines)
1838 for qnum, (past, future) in enumerate(zip(iblocks[:-1], iblocks[1:])):
1839 qpoint = helper._read_qpoints(past)
1840 results[qnum + 1] = curr_result = {"qpoint": qpoint}
1841 for prop in properties:
1842 p = (past < output[prop]) & (output[prop] < future)
1843 selected = output[prop][p]
1844 if len(selected) == 0:
1845 continue
1846 if unique[prop]:
1847 idx = output[prop][p][-1]
1848 curr_result[names[prop]] = properties[prop](idx)
1849 else:
1850 tmp = {k + 1: 0 for k in range(len(selected))}
1851 for k, idx in enumerate(selected):
1852 tmp[k + 1] = properties[prop](idx)
1853 curr_result[names[prop]] = tmp
1854 alat = curr_result.pop("alat", 1.0)
1855 atoms = curr_result.pop("positions", None)
1856 cell = curr_result.pop("cell", np.eye(3))
1857 if atoms:
1858 atoms.positions *= alat * units["Bohr"]
1859 atoms.cell = cell * alat * units["Bohr"]
1860 atoms.wrap()
1861 curr_result["atoms"] = atoms
1863 return results
1866@writer
1867def write_fortran_namelist(
1868 fd,
1869 input_data=None,
1870 binary=None,
1871 additional_cards=None,
1872 **kwargs) -> None:
1873 """
1874 Function which writes input for simple espresso binaries.
1875 List of supported binaries are in the espresso_keys.py file.
1876 Non-exhaustive list (to complete)
1878 Note: "EOF" is appended at the end of the file.
1879 (https://lists.quantum-espresso.org/pipermail/users/2020-November/046269.html)
1881 Additional fields are written 'as is' in the input file. It is expected
1882 to be a string or a list of strings.
1884 Parameters
1885 ----------
1886 fd
1887 The file descriptor of the input file.
1888 input_data: dict
1889 A flat or nested dictionary with input parameters for the binary.x
1890 binary: str
1891 Name of the binary
1892 additional_cards: str | list[str]
1893 Additional fields to be written at the end of the input file, after
1894 the namelist. It is expected to be a string or a list of strings.
1896 Returns
1897 -------
1898 None
1899 """
1900 input_data = Namelist(input_data)
1902 if binary:
1903 input_data.to_nested(binary, **kwargs)
1905 pwi = input_data.to_string()
1907 fd.write(pwi)
1909 if additional_cards:
1910 if isinstance(additional_cards, list):
1911 additional_cards = "\n".join(additional_cards)
1912 additional_cards += "\n"
1914 fd.write(additional_cards)
1916 fd.write("EOF")
1919@deprecated('Please use the ase.io.espresso.Namelist class',
1920 DeprecationWarning)
1921def construct_namelist(parameters=None, keys=None, warn=False, **kwargs):
1922 """
1923 Construct an ordered Namelist containing all the parameters given (as
1924 a dictionary or kwargs). Keys will be inserted into their appropriate
1925 section in the namelist and the dictionary may contain flat and nested
1926 structures. Any kwargs that match input keys will be incorporated into
1927 their correct section. All matches are case-insensitive, and returned
1928 Namelist object is a case-insensitive dict.
1930 If a key is not known to ase, but in a section within `parameters`,
1931 it will be assumed that it was put there on purpose and included
1932 in the output namelist. Anything not in a section will be ignored (set
1933 `warn` to True to see ignored keys).
1935 Keys with a dimension (e.g. Hubbard_U(1)) will be incorporated as-is
1936 so the `i` should be made to match the output.
1938 The priority of the keys is:
1939 kwargs[key] > parameters[key] > parameters[section][key]
1940 Only the highest priority item will be included.
1942 .. deprecated:: 3.23.0
1943 Please use :class:`ase.io.espresso.Namelist` instead.
1945 Parameters
1946 ----------
1947 parameters: dict
1948 Flat or nested set of input parameters.
1949 keys: Namelist | dict
1950 Namelist to use as a template for the output.
1951 warn: bool
1952 Enable warnings for unused keys.
1954 Returns
1955 -------
1956 input_namelist: Namelist
1957 pw.x compatible namelist of input parameters.
1959 """
1961 if keys is None:
1962 keys = deepcopy(pw_keys)
1963 # Convert everything to Namelist early to make case-insensitive
1964 if parameters is None:
1965 parameters = Namelist()
1966 else:
1967 # Maximum one level of nested dict
1968 # Don't modify in place
1969 parameters_namelist = Namelist()
1970 for key, value in parameters.items():
1971 if isinstance(value, dict):
1972 parameters_namelist[key] = Namelist(value)
1973 else:
1974 parameters_namelist[key] = value
1975 parameters = parameters_namelist
1977 # Just a dict
1978 kwargs = Namelist(kwargs)
1980 # Final parameter set
1981 input_namelist = Namelist()
1983 # Collect
1984 for section in keys:
1985 sec_list = Namelist()
1986 for key in keys[section]:
1987 # Check all three separately and pop them all so that
1988 # we can check for missing values later
1989 value = None
1991 if key in parameters.get(section, {}):
1992 value = parameters[section].pop(key)
1993 if key in parameters:
1994 value = parameters.pop(key)
1995 if key in kwargs:
1996 value = kwargs.pop(key)
1998 if value is not None:
1999 sec_list[key] = value
2001 # Check if there is a key(i) version (no extra parsing)
2002 for arg_key in list(parameters.get(section, {})):
2003 if arg_key.split('(')[0].strip().lower() == key.lower():
2004 sec_list[arg_key] = parameters[section].pop(arg_key)
2005 cp_parameters = parameters.copy()
2006 for arg_key in cp_parameters:
2007 if arg_key.split('(')[0].strip().lower() == key.lower():
2008 sec_list[arg_key] = parameters.pop(arg_key)
2009 cp_kwargs = kwargs.copy()
2010 for arg_key in cp_kwargs:
2011 if arg_key.split('(')[0].strip().lower() == key.lower():
2012 sec_list[arg_key] = kwargs.pop(arg_key)
2014 # Add to output
2015 input_namelist[section] = sec_list
2017 unused_keys = list(kwargs)
2018 # pass anything else already in a section
2019 for key, value in parameters.items():
2020 if key in keys and isinstance(value, dict):
2021 input_namelist[key].update(value)
2022 elif isinstance(value, dict):
2023 unused_keys.extend(list(value))
2024 else:
2025 unused_keys.append(key)
2027 if warn and unused_keys:
2028 warnings.warn('Unused keys: {}'.format(', '.join(unused_keys)))
2030 return input_namelist
2033@deprecated('Please use the .to_string() method of Namelist instead.',
2034 DeprecationWarning)
2035def namelist_to_string(input_parameters):
2036 """Format a Namelist object as a string for writing to a file.
2037 Assume sections are ordered (taken care of in namelist construction)
2038 and that repr converts to a QE readable representation (except bools)
2040 .. deprecated:: 3.23.0
2041 Please use the :meth:`ase.io.espresso.Namelist.to_string` method
2042 instead.
2044 Parameters
2045 ----------
2046 input_parameters : Namelist | dict
2047 Expecting a nested dictionary of sections and key-value data.
2049 Returns
2050 -------
2051 pwi : List[str]
2052 Input line for the namelist
2053 """
2054 pwi = []
2055 for section in input_parameters:
2056 pwi.append(f'&{section.upper()}\n')
2057 for key, value in input_parameters[section].items():
2058 if value is True:
2059 pwi.append(f' {key:16} = .true.\n')
2060 elif value is False:
2061 pwi.append(f' {key:16} = .false.\n')
2062 elif isinstance(value, Path):
2063 pwi.append(f' {key:16} = "{value}"\n')
2064 else:
2065 # repr format to get quotes around strings
2066 pwi.append(f' {key:16} = {value!r}\n')
2067 pwi.append('/\n') # terminate section
2068 pwi.append('\n')
2069 return pwi
2072def _get_energy(pwo_lines, energy_index):
2073 return float(pwo_lines[energy_index].split()[-2]) * units['Ry']
2076def _get_forces(pwo_lines, force_index, natoms):
2077 # Before QE 5.3 'negative rho' added 2 lines before forces
2078 # Use exact lines to stop before 'non-local' forces
2079 # in high verbosity
2080 if not pwo_lines[force_index + 2].strip():
2081 force_index += 4
2082 else:
2083 force_index += 2
2084 # assume contiguous
2085 forces = [
2086 [float(x) for x in force_line.split()[-3:]] for force_line
2087 in pwo_lines[force_index:force_index + natoms]]
2088 return np.array(forces) * units['Ry'] / units['Bohr']
2091def _get_stress(pwo_lines, stress_index):
2092 sxx, sxy, sxz = pwo_lines[stress_index + 1].split()[:3]
2093 _, syy, syz = pwo_lines[stress_index + 2].split()[:3]
2094 _, _, szz = pwo_lines[stress_index + 3].split()[:3]
2095 stress = np.array([sxx, syy, szz, syz, sxz, sxy], dtype=float)
2096 # sign convention is opposite of ase
2097 stress *= -1 * units['Ry'] / (units['Bohr'] ** 3)
2098 return stress
2101def _get_magmoms(pwo_lines, magmoms_index, natoms):
2102 return [
2103 float(mag_line.split()[-1]) for mag_line
2104 in pwo_lines[magmoms_index + 1:
2105 magmoms_index + 1 + natoms]]