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