Coverage for /builds/ase/ase/ase/io/espresso.py: 77.53%
868 statements
« prev ^ index » next coverage.py v7.5.3, created at 2025-08-02 00:12 +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 # Extract calculation results
224 # Energy
225 energy = None
226 for energy_index in indexes[_PW_TOTEN]:
227 if image_index < energy_index < next_index:
228 energy = float(
229 pwo_lines[energy_index].split()[-2]) * units['Ry']
231 # Forces
232 forces = None
233 for force_index in indexes[_PW_FORCE]:
234 if image_index < force_index < next_index:
235 # Before QE 5.3 'negative rho' added 2 lines before forces
236 # Use exact lines to stop before 'non-local' forces
237 # in high verbosity
238 if not pwo_lines[force_index + 2].strip():
239 force_index += 4
240 else:
241 force_index += 2
242 # assume contiguous
243 forces = [
244 [float(x) for x in force_line.split()[-3:]] for force_line
245 in pwo_lines[force_index:force_index + len(structure)]]
246 forces = np.array(forces) * units['Ry'] / units['Bohr']
248 # Stress
249 stress = None
250 for stress_index in indexes[_PW_STRESS]:
251 if image_index < stress_index < next_index:
252 sxx, sxy, sxz = pwo_lines[stress_index + 1].split()[:3]
253 _, syy, syz = pwo_lines[stress_index + 2].split()[:3]
254 _, _, szz = pwo_lines[stress_index + 3].split()[:3]
255 stress = np.array([sxx, syy, szz, syz, sxz, sxy], dtype=float)
256 # sign convention is opposite of ase
257 stress *= -1 * units['Ry'] / (units['Bohr'] ** 3)
259 # Magmoms
260 magmoms = None
261 for magmoms_index in indexes[_PW_MAGMOM]:
262 if image_index < magmoms_index < next_index:
263 magmoms = [
264 float(mag_line.split()[-1]) for mag_line
265 in pwo_lines[magmoms_index + 1:
266 magmoms_index + 1 + len(structure)]]
268 # Dipole moment
269 dipole = None
270 if indexes[_PW_DIPOLE]:
271 for dipole_index in indexes[_PW_DIPOLE]:
272 if image_index < dipole_index < next_index:
273 _dipole = float(pwo_lines[dipole_index].split()[-2])
275 for dipole_index in indexes[_PW_DIPOLE_DIRECTION]:
276 if image_index < dipole_index < next_index:
277 _direction = pwo_lines[dipole_index].strip()
278 prefix = 'Computed dipole along edir('
279 _direction = _direction[len(prefix):]
280 _direction = int(_direction[0])
282 dipole = np.eye(3)[_direction - 1] * _dipole * units['Debye']
284 # Fermi level / highest occupied level
285 efermi = None
286 for fermi_index in indexes[_PW_FERMI]:
287 if image_index < fermi_index < next_index:
288 efermi = float(pwo_lines[fermi_index].split()[-2])
290 if efermi is None:
291 for ho_index in indexes[_PW_HIGHEST_OCCUPIED]:
292 if image_index < ho_index < next_index:
293 efermi = float(pwo_lines[ho_index].split()[-1])
295 if efermi is None:
296 for holf_index in indexes[_PW_HIGHEST_OCCUPIED_LOWEST_FREE]:
297 if image_index < holf_index < next_index:
298 efermi = float(pwo_lines[holf_index].split()[-2])
300 # K-points
301 ibzkpts = None
302 weights = None
303 kpoints_warning = "Number of k-points >= 100: " + \
304 "set verbosity='high' to print them."
306 for kpts_index in indexes[_PW_KPTS]:
307 nkpts = int(re.findall(r'\b\d+\b', pwo_lines[kpts_index])[0])
308 kpts_index += 2
310 if pwo_lines[kpts_index].strip() == kpoints_warning:
311 continue
313 # QE prints the k-points in units of 2*pi/alat
314 cell = structure.get_cell()
315 ibzkpts = []
316 weights = []
317 for i in range(nkpts):
318 L = pwo_lines[kpts_index + i].split()
319 weights.append(float(L[-1]))
320 coord = np.array([L[-6], L[-5], L[-4].strip('),')],
321 dtype=float)
322 coord *= 2 * np.pi / cell_alat
323 coord = kpoint_convert(cell, ckpts_kv=coord)
324 ibzkpts.append(coord)
325 ibzkpts = np.array(ibzkpts)
326 weights = np.array(weights)
328 # Bands
329 kpts = None
330 kpoints_warning = "Number of k-points >= 100: " + \
331 "set verbosity='high' to print the bands."
333 for bands_index in indexes[_PW_BANDS] + indexes[_PW_BANDSTRUCTURE]:
334 if image_index < bands_index < next_index:
335 bands_index += 1
336 # skip over the lines with DFT+U occupation matrices
337 if 'enter write_ns' in pwo_lines[bands_index]:
338 while 'exit write_ns' not in pwo_lines[bands_index]:
339 bands_index += 1
340 bands_index += 1
342 if pwo_lines[bands_index].strip() == kpoints_warning:
343 continue
345 assert ibzkpts is not None
346 spin, bands, eigenvalues = 0, [], [[], []]
348 while True:
349 L = pwo_lines[bands_index].replace('-', ' -').split()
350 if len(L) == 0:
351 if len(bands) > 0:
352 eigenvalues[spin].append(bands)
353 bands = []
354 elif L == ['occupation', 'numbers']:
355 # Skip the lines with the occupation numbers
356 bands_index += len(eigenvalues[spin][0]) // 8 + 1
357 elif L[0] == 'k' and L[1].startswith('='):
358 pass
359 elif 'SPIN' in L:
360 if 'DOWN' in L:
361 spin += 1
362 else:
363 try:
364 bands.extend(map(float, L))
365 except ValueError:
366 break
367 bands_index += 1
369 if spin == 1:
370 assert len(eigenvalues[0]) == len(eigenvalues[1])
371 assert len(eigenvalues[0]) == len(ibzkpts), \
372 (np.shape(eigenvalues), len(ibzkpts))
374 kpts = []
375 for s in range(spin + 1):
376 for w, k, e in zip(weights, ibzkpts, eigenvalues[s]):
377 kpt = SinglePointKPoint(w, s, k, eps_n=e)
378 kpts.append(kpt)
380 # Put everything together
381 #
382 # In PW the forces are consistent with the "total energy"; that's why
383 # its value must be assigned to free_energy.
384 # PW doesn't compute the extrapolation of the energy to 0K smearing
385 # the closer thing to this is again the total energy that contains
386 # the correct (i.e. variational) form of the band energy is
387 # Eband = \int e N(e) de for e<Ef , where N(e) is the DOS
388 # This differs by the term (-TS) from the sum of KS eigenvalues:
389 # Eks = \sum wg(n,k) et(n,k)
390 # which is non variational. When a Fermi-Dirac function is used
391 # for a given T, the variational energy is REALLY the free energy F,
392 # and F = E - TS , with E = non variational energy.
393 #
394 calc = SinglePointDFTCalculator(structure, energy=energy,
395 free_energy=energy,
396 forces=forces, stress=stress,
397 magmoms=magmoms, efermi=efermi,
398 ibzkpts=ibzkpts, dipole=dipole)
399 calc.kpts = kpts
400 structure.calc = calc
402 yield structure
405def parse_pwo_start(lines, index=0):
406 """Parse Quantum ESPRESSO calculation info from lines,
407 starting from index. Return a dictionary containing extracted
408 information.
410 - `celldm(1)`: lattice parameters (alat)
411 - `cell`: unit cell in Angstrom
412 - `symbols`: element symbols for the structure
413 - `positions`: cartesian coordinates of atoms in Angstrom
414 - `atoms`: an `ase.Atoms` object constructed from the extracted data
416 Parameters
417 ----------
418 lines : list[str]
419 Contents of PWSCF output file.
420 index : int
421 Line number to begin parsing. Only first calculation will
422 be read.
424 Returns
425 -------
426 info : dict
427 Dictionary of calculation parameters, including `celldm(1)`, `cell`,
428 `symbols`, `positions`, `atoms`.
430 Raises
431 ------
432 KeyError
433 If interdependent values cannot be found (especially celldm(1))
434 an error will be raised as other quantities cannot then be
435 calculated (e.g. cell and positions).
436 """
437 # TODO: extend with extra DFT info?
439 info = {}
441 for idx, line in enumerate(lines[index:], start=index):
442 if 'celldm(1)' in line:
443 # celldm(1) has more digits than alat!!
444 info['celldm(1)'] = float(line.split()[1]) * units['Bohr']
445 info['alat'] = info['celldm(1)']
446 elif 'number of atoms/cell' in line:
447 info['nat'] = int(line.split()[-1])
448 elif 'number of atomic types' in line:
449 info['ntyp'] = int(line.split()[-1])
450 elif 'crystal axes:' in line:
451 info['cell'] = info['celldm(1)'] * np.array([
452 [float(x) for x in lines[idx + 1].split()[3:6]],
453 [float(x) for x in lines[idx + 2].split()[3:6]],
454 [float(x) for x in lines[idx + 3].split()[3:6]]])
455 elif 'positions (alat units)' in line:
456 info['symbols'], info['positions'] = [], []
458 for at_line in lines[idx + 1:idx + 1 + info['nat']]:
459 sym, x, y, z = parse_position_line(at_line)
460 info['symbols'].append(label_to_symbol(sym))
461 info['positions'].append([x * info['celldm(1)'],
462 y * info['celldm(1)'],
463 z * info['celldm(1)']])
464 # This should be the end of interesting info.
465 # Break here to avoid dealing with large lists of kpoints.
466 # Will need to be extended for DFTCalculator info.
467 break
469 # Make atoms for convenience
470 info['atoms'] = Atoms(symbols=info['symbols'],
471 positions=info['positions'],
472 cell=info['cell'], pbc=True)
474 return info
477def parse_position_line(line):
478 """Parse a single line from a pw.x output file.
480 The line must contain information about the atomic symbol and the position,
481 e.g.
483 995 Sb tau( 995) = ( 1.4212023 0.7037863 0.1242640 )
485 Parameters
486 ----------
487 line : str
488 Line to be parsed.
490 Returns
491 -------
492 sym : str
493 Atomic symbol.
494 x : float
495 x-position.
496 y : float
497 y-position.
498 z : float
499 z-position.
500 """
501 pat = re.compile(r'\s*\d+\s*(\S+)\s*tau\(\s*\d+\)\s*='
502 r'\s*\(\s*(\S+)\s+(\S+)\s+(\S+)\s*\)')
503 match = pat.match(line)
504 assert match is not None
505 sym, x, y, z = match.group(1, 2, 3, 4)
506 return sym, float(x), float(y), float(z)
509@reader
510def read_espresso_in(fileobj):
511 """Parse a Quantum ESPRESSO input files, '.in', '.pwi'.
513 ESPRESSO inputs are generally a fortran-namelist format with custom
514 blocks of data. The namelist is parsed as a dict and an atoms object
515 is constructed from the included information.
517 Parameters
518 ----------
519 fileobj : file | str
520 A file-like object that supports line iteration with the contents
521 of the input file, or a filename.
523 Returns
524 -------
525 atoms : Atoms
526 Structure defined in the input file.
528 Raises
529 ------
530 KeyError
531 Raised for missing keys that are required to process the file
532 """
533 # parse namelist section and extract remaining lines
534 data, card_lines = read_fortran_namelist(fileobj)
536 # get the cell if ibrav=0
537 if 'system' not in data:
538 raise KeyError('Required section &SYSTEM not found.')
539 elif 'ibrav' not in data['system']:
540 raise KeyError('ibrav is required in &SYSTEM')
541 elif data['system']['ibrav'] == 0:
542 # celldm(1) is in Bohr, A is in angstrom. celldm(1) will be
543 # used even if A is also specified.
544 if 'celldm(1)' in data['system']:
545 alat = data['system']['celldm(1)'] * units['Bohr']
546 elif 'A' in data['system']:
547 alat = data['system']['A']
548 else:
549 alat = None
550 cell, _ = get_cell_parameters(card_lines, alat=alat)
551 else:
552 raise ValueError(ibrav_error_message)
554 # species_info holds some info for each element
555 species_card = get_atomic_species(
556 card_lines, n_species=data['system']['ntyp'])
557 species_info = {}
558 for ispec, (label, weight, pseudo) in enumerate(species_card):
559 symbol = label_to_symbol(label)
561 # starting_magnetization is in fractions of valence electrons
562 magnet_key = f"starting_magnetization({ispec + 1})"
563 magmom = data["system"].get(magnet_key, 0.0)
564 species_info[symbol] = {"weight": weight, "pseudo": pseudo,
565 "magmom": magmom}
567 positions_card = get_atomic_positions(
568 card_lines, n_atoms=data['system']['nat'], cell=cell, alat=alat)
570 symbols = [label_to_symbol(position[0]) for position in positions_card]
571 positions = [position[1] for position in positions_card]
572 constraint_flags = [position[2] for position in positions_card]
573 magmoms = [species_info[symbol]["magmom"] for symbol in symbols]
575 # TODO: put more info into the atoms object
576 # e.g magmom, forces.
577 atoms = Atoms(symbols=symbols, positions=positions, cell=cell, pbc=True,
578 magmoms=magmoms)
579 atoms.set_constraint(convert_constraint_flags(constraint_flags))
581 return atoms
584def get_atomic_positions(lines, n_atoms, cell=None, alat=None):
585 """Parse atom positions from ATOMIC_POSITIONS card.
587 Parameters
588 ----------
589 lines : list[str]
590 A list of lines containing the ATOMIC_POSITIONS card.
591 n_atoms : int
592 Expected number of atoms. Only this many lines will be parsed.
593 cell : np.array
594 Unit cell of the crystal. Only used with crystal coordinates.
595 alat : float
596 Lattice parameter for atomic coordinates. Only used for alat case.
598 Returns
599 -------
600 positions : list[(str, (float, float, float), (int, int, int))]
601 A list of the ordered atomic positions in the format:
602 label, (x, y, z), (if_x, if_y, if_z)
603 Force multipliers are set to None if not present.
605 Raises
606 ------
607 ValueError
608 Any problems parsing the data result in ValueError
610 """
612 positions = None
613 # no blanks or comment lines, can the consume n_atoms lines for positions
614 trimmed_lines = (line for line in lines if line.strip() and line[0] != '#')
616 for line in trimmed_lines:
617 if line.strip().startswith('ATOMIC_POSITIONS'):
618 if positions is not None:
619 raise ValueError('Multiple ATOMIC_POSITIONS specified')
620 # Priority and behaviour tested with QE 5.3
621 if 'crystal_sg' in line.lower():
622 raise NotImplementedError('CRYSTAL_SG not implemented')
623 elif 'crystal' in line.lower():
624 cell = cell
625 elif 'bohr' in line.lower():
626 cell = np.identity(3) * units['Bohr']
627 elif 'angstrom' in line.lower():
628 cell = np.identity(3)
629 # elif 'alat' in line.lower():
630 # cell = np.identity(3) * alat
631 else:
632 if alat is None:
633 raise ValueError('Set lattice parameter in &SYSTEM for '
634 'alat coordinates')
635 # Always the default, will be DEPRECATED as mandatory
636 # in future
637 cell = np.identity(3) * alat
639 positions = []
640 for _ in range(n_atoms):
641 split_line = next(trimmed_lines).split()
642 # These can be fractions and other expressions
643 position = np.dot((infix_float(split_line[1]),
644 infix_float(split_line[2]),
645 infix_float(split_line[3])), cell)
646 if len(split_line) > 4:
647 force_mult = tuple(int(split_line[i]) for i in (4, 5, 6))
648 else:
649 force_mult = None
651 positions.append((split_line[0], position, force_mult))
653 return positions
656def get_atomic_species(lines, n_species):
657 """Parse atomic species from ATOMIC_SPECIES card.
659 Parameters
660 ----------
661 lines : list[str]
662 A list of lines containing the ATOMIC_POSITIONS card.
663 n_species : int
664 Expected number of atom types. Only this many lines will be parsed.
666 Returns
667 -------
668 species : list[(str, float, str)]
670 Raises
671 ------
672 ValueError
673 Any problems parsing the data result in ValueError
674 """
676 species = None
677 # no blanks or comment lines, can the consume n_atoms lines for positions
678 trimmed_lines = (line.strip() for line in lines
679 if line.strip() and not line.startswith('#'))
681 for line in trimmed_lines:
682 if line.startswith('ATOMIC_SPECIES'):
683 if species is not None:
684 raise ValueError('Multiple ATOMIC_SPECIES specified')
686 species = []
687 for _dummy in range(n_species):
688 label_weight_pseudo = next(trimmed_lines).split()
689 species.append((label_weight_pseudo[0],
690 float(label_weight_pseudo[1]),
691 label_weight_pseudo[2]))
693 return species
696def get_cell_parameters(lines, alat=None):
697 """Parse unit cell from CELL_PARAMETERS card.
699 Parameters
700 ----------
701 lines : list[str]
702 A list with lines containing the CELL_PARAMETERS card.
703 alat : float | None
704 Unit of lattice vectors in Angstrom. Only used if the card is
705 given in units of alat. alat must be None if CELL_PARAMETERS card
706 is in Bohr or Angstrom. For output files, alat will be parsed from
707 the card header and used in preference to this value.
709 Returns
710 -------
711 cell : np.array | None
712 Cell parameters as a 3x3 array in Angstrom. If no cell is found
713 None will be returned instead.
714 cell_alat : float | None
715 If a value for alat is given in the card header, this is also
716 returned, otherwise this will be None.
718 Raises
719 ------
720 ValueError
721 If CELL_PARAMETERS are given in units of bohr or angstrom
722 and alat is not
723 """
725 cell = None
726 cell_alat = None
727 # no blanks or comment lines, can take three lines for cell
728 trimmed_lines = (line for line in lines if line.strip() and line[0] != '#')
730 for line in trimmed_lines:
731 if line.strip().startswith('CELL_PARAMETERS'):
732 if cell is not None:
733 # multiple definitions
734 raise ValueError('CELL_PARAMETERS specified multiple times')
735 # Priority and behaviour tested with QE 5.3
736 if 'bohr' in line.lower():
737 if alat is not None:
738 raise ValueError('Lattice parameters given in '
739 '&SYSTEM celldm/A and CELL_PARAMETERS '
740 'bohr')
741 cell_units = units['Bohr']
742 elif 'angstrom' in line.lower():
743 if alat is not None:
744 raise ValueError('Lattice parameters given in '
745 '&SYSTEM celldm/A and CELL_PARAMETERS '
746 'angstrom')
747 cell_units = 1.0
748 elif 'alat' in line.lower():
749 # Output file has (alat = value) (in Bohrs)
750 if '=' in line:
751 alat = float(line.strip(') \n').split()[-1]) * units['Bohr']
752 cell_alat = alat
753 elif alat is None:
754 raise ValueError('Lattice parameters must be set in '
755 '&SYSTEM for alat units')
756 cell_units = alat
757 elif alat is None:
758 # may be DEPRECATED in future
759 cell_units = units['Bohr']
760 else:
761 # may be DEPRECATED in future
762 cell_units = alat
763 # Grab the parameters; blank lines have been removed
764 cell = [[ffloat(x) for x in next(trimmed_lines).split()[:3]],
765 [ffloat(x) for x in next(trimmed_lines).split()[:3]],
766 [ffloat(x) for x in next(trimmed_lines).split()[:3]]]
767 cell = np.array(cell) * cell_units
769 return cell, cell_alat
772def convert_constraint_flags(constraint_flags):
773 """Convert Quantum ESPRESSO constraint flags to ASE Constraint objects.
775 Parameters
776 ----------
777 constraint_flags : list[tuple[int, int, int]]
778 List of constraint flags (0: fixed, 1: moved) for all the atoms.
779 If the flag is None, there are no constraints on the atom.
781 Returns
782 -------
783 constraints : list[FixAtoms | FixCartesian]
784 List of ASE Constraint objects.
785 """
786 constraints = []
787 for i, constraint in enumerate(constraint_flags):
788 if constraint is None:
789 continue
790 # mask: False (0): moved, True (1): fixed
791 mask = ~np.asarray(constraint, bool)
792 constraints.append(FixCartesian(i, mask))
793 return canonicalize_constraints(constraints)
796def canonicalize_constraints(constraints):
797 """Canonicalize ASE FixCartesian constraints.
799 If the given FixCartesian constraints share the same `mask`, they can be
800 merged into one. Further, if `mask == (True, True, True)`, they can be
801 converted as `FixAtoms`. This method "canonicalizes" FixCartesian objects
802 in such a way.
804 Parameters
805 ----------
806 constraints : List[FixCartesian]
807 List of ASE FixCartesian constraints.
809 Returns
810 -------
811 constrants_canonicalized : List[FixAtoms | FixCartesian]
812 List of ASE Constraint objects.
813 """
814 # https://docs.python.org/3/library/collections.html#defaultdict-examples
815 indices_for_masks = defaultdict(list)
816 for constraint in constraints:
817 key = tuple((constraint.mask).tolist())
818 indices_for_masks[key].extend(constraint.index.tolist())
820 constraints_canonicalized = []
821 for mask, indices in indices_for_masks.items():
822 if mask == (False, False, False): # no directions are fixed
823 continue
824 if mask == (True, True, True): # all three directions are fixed
825 constraints_canonicalized.append(FixAtoms(indices))
826 else:
827 constraints_canonicalized.append(FixCartesian(indices, mask))
829 return constraints_canonicalized
832def str_to_value(string):
833 """Attempt to convert string into int, float (including fortran double),
834 or bool, in that order, otherwise return the string.
835 Valid (case-insensitive) bool values are: '.true.', '.t.', 'true'
836 and 't' (or false equivalents).
838 Parameters
839 ----------
840 string : str
841 Test to parse for a datatype
843 Returns
844 -------
845 value : any
846 Parsed string as the most appropriate datatype of int, float,
847 bool or string.
848 """
850 # Just an integer
851 try:
852 return int(string)
853 except ValueError:
854 pass
855 # Standard float
856 try:
857 return float(string)
858 except ValueError:
859 pass
860 # Fortran double
861 try:
862 return ffloat(string)
863 except ValueError:
864 pass
866 # possible bool, else just the raw string
867 if string.lower() in ('.true.', '.t.', 'true', 't'):
868 return True
869 elif string.lower() in ('.false.', '.f.', 'false', 'f'):
870 return False
871 else:
872 return string.strip("'")
875def read_fortran_namelist(fileobj):
876 """Takes a fortran-namelist formatted file and returns nested
877 dictionaries of sections and key-value data, followed by a list
878 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. All of these are ignored if the value is in 'quotes'.
887 Parameters
888 ----------
889 fileobj : file
890 An open file-like object.
892 Returns
893 -------
894 data : dict[str, dict]
895 Dictionary for each section in the namelist with
896 key = value pairs of data.
897 additional_cards : list[str]
898 Any lines not used to create the data,
899 assumed to belong to 'cards' in the input file.
900 """
902 data = {}
903 card_lines = []
904 in_namelist = False
905 section = 'none' # can't be in a section without changing this
907 for line in fileobj:
908 # leading and trailing whitespace never needed
909 line = line.strip()
910 if line.startswith('&'):
911 # inside a namelist
912 section = line.split()[0][1:].lower() # case insensitive
913 if section in data:
914 # Repeated sections are completely ignored.
915 # (Note that repeated keys overwrite within a section)
916 section = "_ignored"
917 data[section] = {}
918 in_namelist = True
919 if not in_namelist and line:
920 # Stripped line is Truthy, so safe to index first character
921 if line[0] not in ('!', '#'):
922 card_lines.append(line)
923 if in_namelist:
924 # parse k, v from line:
925 key = []
926 value = None
927 in_quotes = False
928 for character in line:
929 if character == ',' and value is not None and not in_quotes:
930 # finished value:
931 data[section][''.join(key).strip()] = str_to_value(
932 ''.join(value).strip())
933 key = []
934 value = None
935 elif character == '=' and value is None and not in_quotes:
936 # start writing value
937 value = []
938 elif character == "'":
939 # only found in value anyway
940 in_quotes = not in_quotes
941 value.append("'")
942 elif character == '!' and not in_quotes:
943 break
944 elif character == '/' and not in_quotes:
945 in_namelist = False
946 break
947 elif value is not None:
948 value.append(character)
949 else:
950 key.append(character)
951 if value is not None:
952 data[section][''.join(key).strip()] = str_to_value(
953 ''.join(value).strip())
955 return Namelist(data), card_lines
958def ffloat(string):
959 """Parse float from fortran compatible float definitions.
961 In fortran exponents can be defined with 'd' or 'q' to symbolise
962 double or quad precision numbers. Double precision numbers are
963 converted to python floats and quad precision values are interpreted
964 as numpy longdouble values (platform specific precision).
966 Parameters
967 ----------
968 string : str
969 A string containing a number in fortran real format
971 Returns
972 -------
973 value : float | np.longdouble
974 Parsed value of the string.
976 Raises
977 ------
978 ValueError
979 Unable to parse a float value.
981 """
983 if 'q' in string.lower():
984 return np.longdouble(string.lower().replace('q', 'e'))
985 else:
986 return float(string.lower().replace('d', 'e'))
989def label_to_symbol(label):
990 """Convert a valid espresso ATOMIC_SPECIES label to a
991 chemical symbol.
993 Parameters
994 ----------
995 label : str
996 chemical symbol X (1 or 2 characters, case-insensitive)
997 or chemical symbol plus a number or a letter, as in
998 "Xn" (e.g. Fe1) or "X_*" or "X-*" (e.g. C1, C_h;
999 max total length cannot exceed 3 characters).
1001 Returns
1002 -------
1003 symbol : str
1004 The best matching species from ase.utils.chemical_symbols
1006 Raises
1007 ------
1008 KeyError
1009 Couldn't find an appropriate species.
1011 Notes
1012 -----
1013 It's impossible to tell whether e.g. He is helium
1014 or hydrogen labelled 'e'.
1015 """
1017 # possibly a two character species
1018 # ase Atoms need proper case of chemical symbols.
1019 if len(label) >= 2:
1020 test_symbol = label[0].upper() + label[1].lower()
1021 if test_symbol in chemical_symbols:
1022 return test_symbol
1023 # finally try with one character
1024 test_symbol = label[0].upper()
1025 if test_symbol in chemical_symbols:
1026 return test_symbol
1027 else:
1028 raise KeyError('Could not parse species from label {}.'
1029 ''.format(label))
1032def infix_float(text):
1033 """Parse simple infix maths into a float for compatibility with
1034 Quantum ESPRESSO ATOMIC_POSITIONS cards. Note: this works with the
1035 example, and most simple expressions, but the capabilities of
1036 the two parsers are not identical. Will also parse a normal float
1037 value properly, but slowly.
1039 >>> infix_float('1/2*3^(-1/2)')
1040 0.28867513459481287
1042 Parameters
1043 ----------
1044 text : str
1045 An arithmetic expression using +, -, *, / and ^, including brackets.
1047 Returns
1048 -------
1049 value : float
1050 Result of the mathematical expression.
1052 """
1054 def middle_brackets(full_text):
1055 """Extract text from innermost brackets."""
1056 start, end = 0, len(full_text)
1057 for (idx, char) in enumerate(full_text):
1058 if char == '(':
1059 start = idx
1060 if char == ')':
1061 end = idx + 1
1062 break
1063 return full_text[start:end]
1065 def eval_no_bracket_expr(full_text):
1066 """Calculate value of a mathematical expression, no brackets."""
1067 exprs = [('+', op.add), ('*', op.mul),
1068 ('/', op.truediv), ('^', op.pow)]
1069 full_text = full_text.lstrip('(').rstrip(')')
1070 try:
1071 return float(full_text)
1072 except ValueError:
1073 for symbol, func in exprs:
1074 if symbol in full_text:
1075 left, right = full_text.split(symbol, 1) # single split
1076 return func(eval_no_bracket_expr(left),
1077 eval_no_bracket_expr(right))
1079 while '(' in text:
1080 middle = middle_brackets(text)
1081 text = text.replace(middle, f'{eval_no_bracket_expr(middle)}')
1083 return float(eval_no_bracket_expr(text))
1086# Number of valence electrons in the pseudopotentials recommended by
1087# http://materialscloud.org/sssp/. These are just used as a fallback for
1088# calculating initial magetization values which are given as a fraction
1089# of valence electrons.
1090SSSP_VALENCE = [
1091 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,
1092 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,
1093 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,
1094 13.0, 14.0, 15.0, 16.0, 17.0, 18.0, 19.0, 12.0, 13.0, 14.0, 15.0, 6.0,
1095 7.0, 18.0, 9.0, 10.0, 11.0, 12.0, 13.0, 14.0, 15.0, 16.0, 17.0, 18.0,
1096 19.0, 20.0, 21.0, 22.0, 23.0, 24.0, 25.0, 36.0, 27.0, 14.0, 15.0, 30.0,
1097 15.0, 32.0, 19.0, 12.0, 13.0, 14.0, 15.0, 16.0, 18.0]
1100def kspacing_to_grid(atoms, spacing, calculated_spacing=None):
1101 """
1102 Calculate the kpoint mesh that is equivalent to the given spacing
1103 in reciprocal space (units Angstrom^-1). The number of kpoints is each
1104 dimension is rounded up (compatible with CASTEP).
1106 Parameters
1107 ----------
1108 atoms: ase.Atoms
1109 A structure that can have get_reciprocal_cell called on it.
1110 spacing: float
1111 Minimum K-Point spacing in $A^{-1}$.
1112 calculated_spacing : list
1113 If a three item list (or similar mutable sequence) is given the
1114 members will be replaced with the actual calculated spacing in
1115 $A^{-1}$.
1117 Returns
1118 -------
1119 kpoint_grid : [int, int, int]
1120 MP grid specification to give the required spacing.
1122 """
1123 # No factor of 2pi in ase, everything in A^-1
1124 # reciprocal dimensions
1125 r_x, r_y, r_z = np.linalg.norm(atoms.cell.reciprocal(), axis=1)
1127 kpoint_grid = [int(r_x / spacing) + 1,
1128 int(r_y / spacing) + 1,
1129 int(r_z / spacing) + 1]
1131 for i, _ in enumerate(kpoint_grid):
1132 if not atoms.pbc[i]:
1133 kpoint_grid[i] = 1
1135 if calculated_spacing is not None:
1136 calculated_spacing[:] = [r_x / kpoint_grid[0],
1137 r_y / kpoint_grid[1],
1138 r_z / kpoint_grid[2]]
1140 return kpoint_grid
1143def format_atom_position(atom, crystal_coordinates, mask='', tidx=None):
1144 """Format one line of atomic positions in
1145 Quantum ESPRESSO ATOMIC_POSITIONS card.
1147 >>> for atom in make_supercell(bulk('Li', 'bcc'), np.ones(3)-np.eye(3)):
1148 >>> format_atom_position(atom, True)
1149 Li 0.0000000000 0.0000000000 0.0000000000
1150 Li 0.5000000000 0.5000000000 0.5000000000
1152 Parameters
1153 ----------
1154 atom : Atom
1155 A structure that has symbol and [position | (a, b, c)].
1156 crystal_coordinates: bool
1157 Whether the atomic positions should be written to the QE input file in
1158 absolute (False, default) or relative (crystal) coordinates (True).
1159 mask, optional : str
1160 String of ndim=3 0 or 1 for constraining atomic positions.
1161 tidx, optional : int
1162 Magnetic type index.
1164 Returns
1165 -------
1166 atom_line : str
1167 Input line for atom position
1168 """
1169 if crystal_coordinates:
1170 coords = [atom.a, atom.b, atom.c]
1171 else:
1172 coords = atom.position
1173 line_fmt = '{atom.symbol}'
1174 inps = dict(atom=atom)
1175 if tidx is not None:
1176 line_fmt += '{tidx}'
1177 inps["tidx"] = tidx
1178 line_fmt += ' {coords[0]:.10f} {coords[1]:.10f} {coords[2]:.10f} '
1179 inps["coords"] = coords
1180 line_fmt += ' ' + mask + '\n'
1181 astr = line_fmt.format(**inps)
1182 return astr
1185@writer
1186def write_espresso_in(fd, atoms, input_data=None, pseudopotentials=None,
1187 kspacing=None, kpts=None, koffset=(0, 0, 0),
1188 crystal_coordinates=False, additional_cards=None,
1189 **kwargs):
1190 """
1191 Create an input file for pw.x.
1193 Use set_initial_magnetic_moments to turn on spin, if nspin is set to 2
1194 with no magnetic moments, they will all be set to 0.0. Magnetic moments
1195 will be converted to the QE units (fraction of valence electrons) using
1196 any pseudopotential files found, or a best guess for the number of
1197 valence electrons.
1199 Units are not converted for any other input data, so use Quantum ESPRESSO
1200 units (Usually Ry or atomic units).
1202 Keys with a dimension (e.g. Hubbard_U(1)) will be incorporated as-is
1203 so the `i` should be made to match the output.
1205 Implemented features:
1207 - Conversion of :class:`ase.constraints.FixAtoms` and
1208 :class:`ase.constraints.FixCartesian`.
1209 - ``starting_magnetization`` derived from the ``magmoms`` and
1210 pseudopotentials (searches default paths for pseudo files.)
1211 - Automatic assignment of options to their correct sections.
1213 Not implemented:
1215 - Non-zero values of ibrav
1216 - Lists of k-points
1217 - Other constraints
1218 - Hubbard parameters
1219 - Validation of the argument types for input
1220 - Validation of required options
1222 Parameters
1223 ----------
1224 fd: file | str
1225 A file to which the input is written.
1226 atoms: Atoms
1227 A single atomistic configuration to write to ``fd``.
1228 input_data: dict
1229 A flat or nested dictionary with input parameters for pw.x
1230 pseudopotentials: dict
1231 A filename for each atomic species, e.g.
1232 {'O': 'O.pbe-rrkjus.UPF', 'H': 'H.pbe-rrkjus.UPF'}.
1233 A dummy name will be used if none are given.
1234 kspacing: float
1235 Generate a grid of k-points with this as the minimum distance,
1236 in A^-1 between them in reciprocal space. If set to None, kpts
1237 will be used instead.
1238 kpts: (int, int, int), dict or np.ndarray
1239 If ``kpts`` is a tuple (or list) of 3 integers, it is interpreted
1240 as the dimensions of a Monkhorst-Pack grid.
1241 If ``kpts`` is set to ``None``, only the Γ-point will be included
1242 and QE will use routines optimized for Γ-point-only calculations.
1243 Compared to Γ-point-only calculations without this optimization
1244 (i.e. with ``kpts=(1, 1, 1)``), the memory and CPU requirements
1245 are typically reduced by half.
1246 If kpts is a dict, it will either be interpreted as a path
1247 in the Brillouin zone (*) if it contains the 'path' keyword,
1248 otherwise it is converted to a Monkhorst-Pack grid (**).
1249 If ``kpts`` is a NumPy array, the raw k-points will be passed to
1250 Quantum Espresso as given in the array (in crystal coordinates).
1251 Must be of shape (n_kpts, 4). The fourth column contains the
1252 k-point weights.
1253 (*) see ase.dft.kpoints.bandpath
1254 (**) see ase.calculators.calculator.kpts2sizeandoffsets
1255 koffset: (int, int, int)
1256 Offset of kpoints in each direction. Must be 0 (no offset) or
1257 1 (half grid offset). Setting to True is equivalent to (1, 1, 1).
1258 crystal_coordinates: bool
1259 Whether the atomic positions should be written to the QE input file in
1260 absolute (False, default) or relative (crystal) coordinates (True).
1262 """
1264 # Convert to a namelist to make working with parameters much easier
1265 # Note that the name ``input_data`` is chosen to prevent clash with
1266 # ``parameters`` in Calculator objects
1267 input_parameters = Namelist(input_data)
1268 input_parameters.to_nested('pw', **kwargs)
1270 # Convert ase constraints to QE constraints
1271 # Nx3 array of force multipliers matches what QE uses
1272 # Do this early so it is available when constructing the atoms card
1273 moved = np.ones((len(atoms), 3), dtype=bool)
1274 for constraint in atoms.constraints:
1275 if isinstance(constraint, FixAtoms):
1276 moved[constraint.index] = False
1277 elif isinstance(constraint, FixCartesian):
1278 moved[constraint.index] = ~constraint.mask
1279 else:
1280 warnings.warn(f'Ignored unknown constraint {constraint}')
1281 masks = []
1282 for atom in atoms:
1283 # only inclued mask if something is fixed
1284 if not all(moved[atom.index]):
1285 mask = ' {:d} {:d} {:d}'.format(*moved[atom.index])
1286 else:
1287 mask = ''
1288 masks.append(mask)
1290 # Species info holds the information on the pseudopotential and
1291 # associated for each element
1292 if pseudopotentials is None:
1293 pseudopotentials = {}
1294 species_info = {}
1295 for species in set(atoms.get_chemical_symbols()):
1296 # Look in all possible locations for the pseudos and try to figure
1297 # out the number of valence electrons
1298 pseudo = pseudopotentials[species]
1299 species_info[species] = {'pseudo': pseudo}
1301 # Convert atoms into species.
1302 # Each different magnetic moment needs to be a separate type even with
1303 # the same pseudopotential (e.g. an up and a down for AFM).
1304 # if any magmom are > 0 or nspin == 2 then use species labels.
1305 # Rememeber: magnetisation uses 1 based indexes
1306 atomic_species = {}
1307 atomic_species_str = []
1308 atomic_positions_str = []
1310 nspin = input_parameters['system'].get('nspin', 1) # 1 is the default
1311 noncolin = input_parameters['system'].get('noncolin', False)
1312 rescale_magmom_fac = kwargs.get('rescale_magmom_fac', 1.0)
1313 if any(atoms.get_initial_magnetic_moments()):
1314 if nspin == 1 and not noncolin:
1315 # Force spin on
1316 input_parameters['system']['nspin'] = 2
1317 nspin = 2
1319 if nspin == 2 or noncolin:
1320 # Magnetic calculation on
1321 for atom, mask, magmom in zip(
1322 atoms, masks, atoms.get_initial_magnetic_moments()):
1323 if (atom.symbol, magmom) not in atomic_species:
1324 # for qe version 7.2 or older magmon must be rescale by
1325 # about a factor 10 to assume sensible values
1326 # since qe-v7.3 magmom values will be provided unscaled
1327 fspin = float(magmom) / rescale_magmom_fac
1328 # Index in the atomic species list
1329 sidx = len(atomic_species) + 1
1330 # Index for that atom type; no index for first one
1331 tidx = sum(atom.symbol == x[0] for x in atomic_species) or ' '
1332 atomic_species[(atom.symbol, magmom)] = (sidx, tidx)
1333 # Add magnetization to the input file
1334 mag_str = f"starting_magnetization({sidx})"
1335 input_parameters['system'][mag_str] = fspin
1336 species_pseudo = species_info[atom.symbol]['pseudo']
1337 atomic_species_str.append(
1338 f"{atom.symbol}{tidx} {atom.mass} {species_pseudo}\n")
1339 # lookup tidx to append to name
1340 sidx, tidx = atomic_species[(atom.symbol, magmom)]
1341 # construct line for atomic positions
1342 atomic_positions_str.append(
1343 format_atom_position(
1344 atom, crystal_coordinates, mask=mask, tidx=tidx)
1345 )
1346 else:
1347 # Do nothing about magnetisation
1348 for atom, mask in zip(atoms, masks):
1349 if atom.symbol not in atomic_species:
1350 atomic_species[atom.symbol] = True # just a placeholder
1351 species_pseudo = species_info[atom.symbol]['pseudo']
1352 atomic_species_str.append(
1353 f"{atom.symbol} {atom.mass} {species_pseudo}\n")
1354 # construct line for atomic positions
1355 atomic_positions_str.append(
1356 format_atom_position(atom, crystal_coordinates, mask=mask)
1357 )
1359 # Add computed parameters
1360 # different magnetisms means different types
1361 input_parameters['system']['ntyp'] = len(atomic_species)
1362 input_parameters['system']['nat'] = len(atoms)
1364 # Use cell as given or fit to a specific ibrav
1365 if 'ibrav' in input_parameters['system']:
1366 ibrav = input_parameters['system']['ibrav']
1367 if ibrav != 0:
1368 raise ValueError(ibrav_error_message)
1369 else:
1370 # Just use standard cell block
1371 input_parameters['system']['ibrav'] = 0
1373 # Construct input file into this
1374 pwi = input_parameters.to_string(list_form=True)
1376 # Pseudopotentials
1377 pwi.append('ATOMIC_SPECIES\n')
1378 pwi.extend(atomic_species_str)
1379 pwi.append('\n')
1381 # KPOINTS - add a MP grid as required
1382 if kspacing is not None:
1383 kgrid = kspacing_to_grid(atoms, kspacing)
1384 elif kpts is not None:
1385 if isinstance(kpts, dict) and 'path' not in kpts:
1386 kgrid, shift = kpts2sizeandoffsets(atoms=atoms, **kpts)
1387 koffset = []
1388 for i, x in enumerate(shift):
1389 assert x == 0 or abs(x * kgrid[i] - 0.5) < 1e-14
1390 koffset.append(0 if x == 0 else 1)
1391 else:
1392 kgrid = kpts
1393 else:
1394 kgrid = "gamma"
1396 # True and False work here and will get converted by ':d' format
1397 if isinstance(koffset, int):
1398 koffset = (koffset, ) * 3
1400 # BandPath object or bandpath-as-dictionary:
1401 if isinstance(kgrid, dict) or hasattr(kgrid, 'kpts'):
1402 pwi.append('K_POINTS crystal_b\n')
1403 assert hasattr(kgrid, 'path') or 'path' in kgrid
1404 kgrid = kpts2ndarray(kgrid, atoms=atoms)
1405 pwi.append(f'{len(kgrid)}\n')
1406 for k in kgrid:
1407 pwi.append(f"{k[0]:.14f} {k[1]:.14f} {k[2]:.14f} 0\n")
1408 pwi.append('\n')
1409 elif isinstance(kgrid, str) and (kgrid == "gamma"):
1410 pwi.append('K_POINTS gamma\n')
1411 pwi.append('\n')
1412 elif isinstance(kgrid, np.ndarray):
1413 if np.shape(kgrid)[1] != 4:
1414 raise ValueError('Only Nx4 kgrids are supported right now.')
1415 pwi.append('K_POINTS crystal\n')
1416 pwi.append(f'{len(kgrid)}\n')
1417 for k in kgrid:
1418 pwi.append(f"{k[0]:.14f} {k[1]:.14f} {k[2]:.14f} {k[3]:.14f}\n")
1419 pwi.append('\n')
1420 else:
1421 pwi.append('K_POINTS automatic\n')
1422 pwi.append(f"{kgrid[0]} {kgrid[1]} {kgrid[2]} "
1423 f" {koffset[0]:d} {koffset[1]:d} {koffset[2]:d}\n")
1424 pwi.append('\n')
1426 # CELL block, if required
1427 if input_parameters['SYSTEM']['ibrav'] == 0:
1428 pwi.append('CELL_PARAMETERS angstrom\n')
1429 pwi.append('{cell[0][0]:.14f} {cell[0][1]:.14f} {cell[0][2]:.14f}\n'
1430 '{cell[1][0]:.14f} {cell[1][1]:.14f} {cell[1][2]:.14f}\n'
1431 '{cell[2][0]:.14f} {cell[2][1]:.14f} {cell[2][2]:.14f}\n'
1432 ''.format(cell=atoms.cell))
1433 pwi.append('\n')
1435 # Positions - already constructed, but must appear after namelist
1436 if crystal_coordinates:
1437 pwi.append('ATOMIC_POSITIONS crystal\n')
1438 else:
1439 pwi.append('ATOMIC_POSITIONS angstrom\n')
1440 pwi.extend(atomic_positions_str)
1441 pwi.append('\n')
1443 # DONE!
1444 fd.write(''.join(pwi))
1446 if additional_cards:
1447 if isinstance(additional_cards, list):
1448 additional_cards = "\n".join(additional_cards)
1449 additional_cards += "\n"
1451 fd.write(additional_cards)
1454def write_espresso_ph(
1455 fd,
1456 input_data=None,
1457 qpts=None,
1458 nat_todo_indices=None,
1459 **kwargs) -> None:
1460 """
1461 Function that write the input file for a ph.x calculation. Normal namelist
1462 cards are passed in the input_data dictionary. Which can be either nested
1463 or flat, ASE style. The q-points are passed in the qpts list. If qplot is
1464 set to True then qpts is expected to be a list of list|tuple of length 4.
1465 Where the first three elements are the coordinates of the q-point in units
1466 of 2pi/alat and the last element is the weight of the q-point. if qplot is
1467 set to False then qpts is expected to be a simple list of length 4 (single
1468 q-point). Finally if ldisp is set to True, the above is discarded and the
1469 q-points are read from the nq1, nq2, nq3 cards in the input_data dictionary.
1471 Additionally, a nat_todo_indices kwargs (list[int]) can be specified in the
1472 kwargs. It will be used if nat_todo is set to True in the input_data
1473 dictionary.
1475 Globally, this function follows the convention set in the ph.x documentation
1476 (https://www.quantum-espresso.org/Doc/INPUT_PH.html)
1478 Parameters
1479 ----------
1480 fd
1481 The file descriptor of the input file.
1483 kwargs
1484 kwargs dictionary possibly containing the following keys:
1486 - input_data: dict
1487 - qpts: list[list[float]] | list[tuple[float]] | list[float]
1488 - nat_todo_indices: list[int]
1490 Returns
1491 -------
1492 None
1493 """
1495 input_data = Namelist(input_data)
1496 input_data.to_nested('ph', **kwargs)
1498 input_ph = input_data["inputph"]
1500 inp_nat_todo = input_ph.get("nat_todo", 0)
1501 qpts = qpts or (0, 0, 0)
1503 pwi = input_data.to_string()
1505 fd.write(pwi)
1507 qplot = input_ph.get("qplot", False)
1508 ldisp = input_ph.get("ldisp", False)
1510 if qplot:
1511 fd.write(f"{len(qpts)}\n")
1512 for qpt in qpts:
1513 fd.write(
1514 f"{qpt[0]:0.8f} {qpt[1]:0.8f} {qpt[2]:0.8f} {qpt[3]:1d}\n"
1515 )
1516 elif not (qplot or ldisp):
1517 fd.write(f"{qpts[0]:0.8f} {qpts[1]:0.8f} {qpts[2]:0.8f}\n")
1518 if inp_nat_todo:
1519 tmp = [str(i) for i in nat_todo_indices]
1520 fd.write(" ".join(tmp))
1521 fd.write("\n")
1524def read_espresso_ph(fileobj):
1525 """
1526 Function that reads the output of a ph.x calculation.
1527 It returns a dictionary where each q-point number is a key and
1528 the value is a dictionary with the following keys if available:
1530 - qpoints: The q-point in cartesian coordinates.
1531 - kpoints: The k-points in cartesian coordinates.
1532 - dieltensor: The dielectric tensor.
1533 - borneffcharge: The effective Born charges.
1534 - borneffcharge_dfpt: The effective Born charges from DFPT.
1535 - polarizability: The polarizability tensor.
1536 - modes: The phonon modes.
1537 - eqpoints: The symmetrically equivalent q-points.
1538 - freqs: The phonon frequencies.
1539 - mode_symmetries: The symmetries of the modes.
1540 - atoms: The atoms object.
1542 Some notes:
1544 - For some reason, the cell is not defined to high level of
1545 precision in ph.x outputs. Be careful when using the atoms object
1546 retrieved from this function.
1547 - This function can be called on incomplete calculations i.e.
1548 if the calculation couldn't diagonalize the dynamical matrix
1549 for some q-points, the results for the other q-points will
1550 still be returned.
1552 Parameters
1553 ----------
1554 fileobj
1555 The file descriptor of the output file.
1557 Returns
1558 -------
1559 dict
1560 The results dictionnary as described above.
1561 """
1562 freg = re.compile(r"-?(?:0|[1-9]\d*)(?:\.\d+)?(?:[eE][+\-]?\d+)?")
1564 QPOINTS = r"(?i)^\s*Calculation\s*of\s*q"
1565 NKPTS = r"(?i)^\s*number\s*of\s*k\s*points\s*"
1566 DIEL = r"(?i)^\s*Dielectric\s*constant\s*in\s*cartesian\s*axis\s*$"
1567 BORN = r"(?i)^\s*Effective\s*charges\s*\(d\s*Force\s*/\s*dE\)"
1568 POLA = r"(?i)^\s*Polarizability\s*(a.u.)\^3"
1569 REPR = r"(?i)^\s*There\s*are\s*\d+\s*irreducible\s*representations\s*$"
1570 EQPOINTS = r"(?i)^\s*Number\s*of\s*q\s*in\s*the\s*star\s*=\s*"
1571 DIAG = r"(?i)^\s*Diagonalizing\s*the\s*dynamical\s*matrix\s*$"
1572 MODE_SYM = r"(?i)^\s*Mode\s*symmetry,\s*"
1573 BORN_DFPT = r"(?i)^\s*Effective\s*charges\s*\(d\s*P\s*/\s*du\)"
1574 POSITIONS = r"(?i)^\s*site\s*n\..*\(alat\s*units\)"
1575 ALAT = r"(?i)^\s*celldm\(1\)="
1576 CELL = (
1577 r"^\s*crystal\s*axes:\s*\(cart.\s*coord.\s*in\s*units\s*of\s*alat\)"
1578 )
1579 ELECTRON_PHONON = r"(?i)^\s*electron-phonon\s*interaction\s*...\s*$"
1581 output = {
1582 QPOINTS: [],
1583 NKPTS: [],
1584 DIEL: [],
1585 BORN: [],
1586 BORN_DFPT: [],
1587 POLA: [],
1588 REPR: [],
1589 EQPOINTS: [],
1590 DIAG: [],
1591 MODE_SYM: [],
1592 POSITIONS: [],
1593 ALAT: [],
1594 CELL: [],
1595 ELECTRON_PHONON: [],
1596 }
1598 names = {
1599 QPOINTS: "qpoints",
1600 NKPTS: "kpoints",
1601 DIEL: "dieltensor",
1602 BORN: "borneffcharge",
1603 BORN_DFPT: "borneffcharge_dfpt",
1604 POLA: "polarizability",
1605 REPR: "representations",
1606 EQPOINTS: "eqpoints",
1607 DIAG: "freqs",
1608 MODE_SYM: "mode_symmetries",
1609 POSITIONS: "positions",
1610 ALAT: "alat",
1611 CELL: "cell",
1612 ELECTRON_PHONON: "ep_data",
1613 }
1615 unique = {
1616 QPOINTS: True,
1617 NKPTS: False,
1618 DIEL: True,
1619 BORN: True,
1620 BORN_DFPT: True,
1621 POLA: True,
1622 REPR: True,
1623 EQPOINTS: True,
1624 DIAG: True,
1625 MODE_SYM: True,
1626 POSITIONS: True,
1627 ALAT: True,
1628 CELL: True,
1629 ELECTRON_PHONON: True,
1630 }
1632 results = {}
1633 fdo_lines = [i for i in fileobj.read().splitlines() if i]
1634 n_lines = len(fdo_lines)
1636 for idx, line in enumerate(fdo_lines):
1637 for key in output:
1638 if bool(re.match(key, line)):
1639 output[key].append(idx)
1641 output = {key: np.array(value) for key, value in output.items()}
1643 def _read_qpoints(idx):
1644 match = re.findall(freg, fdo_lines[idx])
1645 return tuple(round(float(x), 7) for x in match)
1647 def _read_kpoints(idx):
1648 n_kpts = int(re.findall(freg, fdo_lines[idx])[0])
1649 kpts = []
1650 for line in fdo_lines[idx + 2: idx + 2 + n_kpts]:
1651 if bool(re.search(r"^\s*k\(.*wk", line)):
1652 kpts.append([round(float(x), 7)
1653 for x in re.findall(freg, line)[1:]])
1654 return np.array(kpts)
1656 def _read_repr(idx):
1657 n_repr, curr, n = int(re.findall(freg, fdo_lines[idx])[0]), 0, 0
1658 representations = {}
1659 while idx + n < n_lines:
1660 if re.search(r"^\s*Representation.*modes", fdo_lines[idx + n]):
1661 curr = int(re.findall(freg, fdo_lines[idx + n])[0])
1662 representations[curr] = {"done": False, "modes": []}
1663 if re.search(r"Calculated\s*using\s*symmetry", fdo_lines[idx + n]) \
1664 or re.search(r"-\s*Done\s*$", fdo_lines[idx + n]):
1665 representations[curr]["done"] = True
1666 if re.search(r"(?i)^\s*(mode\s*#\s*\d\s*)+", fdo_lines[idx + n]):
1667 representations[curr]["modes"] = _read_modes(idx + n)
1668 if curr == n_repr:
1669 break
1670 n += 1
1671 return representations
1673 def _read_modes(idx):
1674 n = 1
1675 n_modes = len(re.findall(r"mode", fdo_lines[idx]))
1676 modes = []
1677 while not modes or bool(re.match(r"^\s*\(", fdo_lines[idx + n])):
1678 tmp = re.findall(freg, fdo_lines[idx + n])
1679 modes.append([round(float(x), 5) for x in tmp])
1680 n += 1
1681 return np.hsplit(np.array(modes), n_modes)
1683 def _read_eqpoints(idx):
1684 n_star = int(re.findall(freg, fdo_lines[idx])[0])
1685 return np.loadtxt(
1686 fdo_lines[idx + 2: idx + 2 + n_star], usecols=(1, 2, 3)
1687 ).reshape(-1, 3)
1689 def _read_freqs(idx):
1690 n = 0
1691 freqs = []
1692 stop = 0
1693 while not freqs or stop < 2:
1694 if bool(re.search(r"^\s*freq", fdo_lines[idx + n])):
1695 tmp = re.findall(freg, fdo_lines[idx + n])[1]
1696 freqs.append(float(tmp))
1697 if bool(re.search(r"\*{5,}", fdo_lines[idx + n])):
1698 stop += 1
1699 n += 1
1700 return np.array(freqs)
1702 def _read_sym(idx):
1703 n = 1
1704 sym = {}
1705 while bool(re.match(r"^\s*freq", fdo_lines[idx + n])):
1706 r = re.findall("\\d+", fdo_lines[idx + n])
1707 r = tuple(range(int(r[0]), int(r[1]) + 1))
1708 sym[r] = fdo_lines[idx + n].split("-->")[1].strip()
1709 sym[r] = re.sub(r"\s+", " ", sym[r])
1710 n += 1
1711 return sym
1713 def _read_epsil(idx):
1714 epsil = np.zeros((3, 3))
1715 for n in range(1, 4):
1716 tmp = re.findall(freg, fdo_lines[idx + n])
1717 epsil[n - 1] = [round(float(x), 9) for x in tmp]
1718 return epsil
1720 def _read_born(idx):
1721 n = 1
1722 born = []
1723 while idx + n < n_lines:
1724 if re.search(r"^\s*atom\s*\d\s*\S", fdo_lines[idx + n]):
1725 pass
1726 elif re.search(r"^\s*E\*?(x|y|z)\s*\(", fdo_lines[idx + n]):
1727 tmp = re.findall(freg, fdo_lines[idx + n])
1728 born.append([round(float(x), 5) for x in tmp])
1729 else:
1730 break
1731 n += 1
1732 born = np.array(born)
1733 return np.vsplit(born, len(born) // 3)
1735 def _read_born_dfpt(idx):
1736 n = 1
1737 born = []
1738 while idx + n < n_lines:
1739 if re.search(r"^\s*atom\s*\d\s*\S", fdo_lines[idx + n]):
1740 pass
1741 elif re.search(r"^\s*P(x|y|z)\s*\(", fdo_lines[idx + n]):
1742 tmp = re.findall(freg, fdo_lines[idx + n])
1743 born.append([round(float(x), 5) for x in tmp])
1744 else:
1745 break
1746 n += 1
1747 born = np.array(born)
1748 return np.vsplit(born, len(born) // 3)
1750 def _read_pola(idx):
1751 pola = np.zeros((3, 3))
1752 for n in range(1, 4):
1753 tmp = re.findall(freg, fdo_lines[idx + n])[:3]
1754 pola[n - 1] = [round(float(x), 2) for x in tmp]
1755 return pola
1757 def _read_positions(idx):
1758 positions = []
1759 symbols = []
1760 n = 1
1761 while re.findall(r"^\s*\d+", fdo_lines[idx + n]):
1762 symbols.append(fdo_lines[idx + n].split()[1])
1763 positions.append(
1764 [round(float(x), 5)
1765 for x in re.findall(freg, fdo_lines[idx + n])[-3:]]
1766 )
1767 n += 1
1768 atoms = Atoms(positions=positions, symbols=symbols)
1769 atoms.pbc = True
1770 return atoms
1772 def _read_alat(idx):
1773 return round(float(re.findall(freg, fdo_lines[idx])[1]), 5)
1775 def _read_cell(idx):
1776 cell = []
1777 n = 1
1778 while re.findall(r"^\s*a\(\d\)", fdo_lines[idx + n]):
1779 cell.append([round(float(x), 4)
1780 for x in re.findall(freg, fdo_lines[idx + n])[-3:]])
1781 n += 1
1782 return np.array(cell)
1784 def _read_electron_phonon(idx):
1785 results = {}
1787 broad_re = (
1788 r"^\s*Gaussian\s*Broadening:\s+([\d.]+)\s+Ry, ngauss=\s+\d+"
1789 )
1790 dos_re = (
1791 r"^\s*DOS\s*=\s*([\d.]+)\s*states/"
1792 r"spin/Ry/Unit\s*Cell\s*at\s*Ef=\s+([\d.]+)\s+eV"
1793 )
1794 lg_re = (
1795 r"^\s*lambda\(\s+(\d+)\)=\s+([\d.]+)\s+gamma=\s+([\d.]+)\s+GHz"
1796 )
1797 end_re = r"^\s*Number\s*of\s*q\s*in\s*the\s*star\s*=\s+(\d+)$"
1799 lambdas = []
1800 gammas = []
1802 current = None
1804 n = 1
1805 while idx + n < n_lines:
1806 line = fdo_lines[idx + n]
1808 broad_match = re.match(broad_re, line)
1809 dos_match = re.match(dos_re, line)
1810 lg_match = re.match(lg_re, line)
1811 end_match = re.match(end_re, line)
1813 if broad_match:
1814 if lambdas:
1815 results[current]["lambdas"] = lambdas
1816 results[current]["gammas"] = gammas
1817 lambdas = []
1818 gammas = []
1819 current = float(broad_match[1])
1820 results[current] = {}
1821 elif dos_match:
1822 results[current]["dos"] = float(dos_match[1])
1823 results[current]["fermi"] = float(dos_match[2])
1824 elif lg_match:
1825 lambdas.append(float(lg_match[2]))
1826 gammas.append(float(lg_match[3]))
1828 if end_match:
1829 results[current]["lambdas"] = lambdas
1830 results[current]["gammas"] = gammas
1831 break
1833 n += 1
1835 return results
1837 properties = {
1838 NKPTS: _read_kpoints,
1839 DIEL: _read_epsil,
1840 BORN: _read_born,
1841 BORN_DFPT: _read_born_dfpt,
1842 POLA: _read_pola,
1843 REPR: _read_repr,
1844 EQPOINTS: _read_eqpoints,
1845 DIAG: _read_freqs,
1846 MODE_SYM: _read_sym,
1847 POSITIONS: _read_positions,
1848 ALAT: _read_alat,
1849 CELL: _read_cell,
1850 ELECTRON_PHONON: _read_electron_phonon,
1851 }
1853 iblocks = np.append(output[QPOINTS], n_lines)
1855 for qnum, (past, future) in enumerate(zip(iblocks[:-1], iblocks[1:])):
1856 qpoint = _read_qpoints(past)
1857 results[qnum + 1] = curr_result = {"qpoint": qpoint}
1858 for prop in properties:
1859 p = (past < output[prop]) & (output[prop] < future)
1860 selected = output[prop][p]
1861 if len(selected) == 0:
1862 continue
1863 if unique[prop]:
1864 idx = output[prop][p][-1]
1865 curr_result[names[prop]] = properties[prop](idx)
1866 else:
1867 tmp = {k + 1: 0 for k in range(len(selected))}
1868 for k, idx in enumerate(selected):
1869 tmp[k + 1] = properties[prop](idx)
1870 curr_result[names[prop]] = tmp
1871 alat = curr_result.pop("alat", 1.0)
1872 atoms = curr_result.pop("positions", None)
1873 cell = curr_result.pop("cell", np.eye(3))
1874 if atoms:
1875 atoms.positions *= alat * units["Bohr"]
1876 atoms.cell = cell * alat * units["Bohr"]
1877 atoms.wrap()
1878 curr_result["atoms"] = atoms
1880 return results
1883def write_fortran_namelist(
1884 fd,
1885 input_data=None,
1886 binary=None,
1887 additional_cards=None,
1888 **kwargs) -> None:
1889 """
1890 Function which writes input for simple espresso binaries.
1891 List of supported binaries are in the espresso_keys.py file.
1892 Non-exhaustive list (to complete)
1894 Note: "EOF" is appended at the end of the file.
1895 (https://lists.quantum-espresso.org/pipermail/users/2020-November/046269.html)
1897 Additional fields are written 'as is' in the input file. It is expected
1898 to be a string or a list of strings.
1900 Parameters
1901 ----------
1902 fd
1903 The file descriptor of the input file.
1904 input_data: dict
1905 A flat or nested dictionary with input parameters for the binary.x
1906 binary: str
1907 Name of the binary
1908 additional_cards: str | list[str]
1909 Additional fields to be written at the end of the input file, after
1910 the namelist. It is expected to be a string or a list of strings.
1912 Returns
1913 -------
1914 None
1915 """
1916 input_data = Namelist(input_data)
1918 if binary:
1919 input_data.to_nested(binary, **kwargs)
1921 pwi = input_data.to_string()
1923 fd.write(pwi)
1925 if additional_cards:
1926 if isinstance(additional_cards, list):
1927 additional_cards = "\n".join(additional_cards)
1928 additional_cards += "\n"
1930 fd.write(additional_cards)
1932 fd.write("EOF")
1935@deprecated('Please use the ase.io.espresso.Namelist class',
1936 DeprecationWarning)
1937def construct_namelist(parameters=None, keys=None, warn=False, **kwargs):
1938 """
1939 Construct an ordered Namelist containing all the parameters given (as
1940 a dictionary or kwargs). Keys will be inserted into their appropriate
1941 section in the namelist and the dictionary may contain flat and nested
1942 structures. Any kwargs that match input keys will be incorporated into
1943 their correct section. All matches are case-insensitive, and returned
1944 Namelist object is a case-insensitive dict.
1946 If a key is not known to ase, but in a section within `parameters`,
1947 it will be assumed that it was put there on purpose and included
1948 in the output namelist. Anything not in a section will be ignored (set
1949 `warn` to True to see ignored keys).
1951 Keys with a dimension (e.g. Hubbard_U(1)) will be incorporated as-is
1952 so the `i` should be made to match the output.
1954 The priority of the keys is:
1955 kwargs[key] > parameters[key] > parameters[section][key]
1956 Only the highest priority item will be included.
1958 .. deprecated:: 3.23.0
1959 Please use :class:`ase.io.espresso.Namelist` instead.
1961 Parameters
1962 ----------
1963 parameters: dict
1964 Flat or nested set of input parameters.
1965 keys: Namelist | dict
1966 Namelist to use as a template for the output.
1967 warn: bool
1968 Enable warnings for unused keys.
1970 Returns
1971 -------
1972 input_namelist: Namelist
1973 pw.x compatible namelist of input parameters.
1975 """
1977 if keys is None:
1978 keys = deepcopy(pw_keys)
1979 # Convert everything to Namelist early to make case-insensitive
1980 if parameters is None:
1981 parameters = Namelist()
1982 else:
1983 # Maximum one level of nested dict
1984 # Don't modify in place
1985 parameters_namelist = Namelist()
1986 for key, value in parameters.items():
1987 if isinstance(value, dict):
1988 parameters_namelist[key] = Namelist(value)
1989 else:
1990 parameters_namelist[key] = value
1991 parameters = parameters_namelist
1993 # Just a dict
1994 kwargs = Namelist(kwargs)
1996 # Final parameter set
1997 input_namelist = Namelist()
1999 # Collect
2000 for section in keys:
2001 sec_list = Namelist()
2002 for key in keys[section]:
2003 # Check all three separately and pop them all so that
2004 # we can check for missing values later
2005 value = None
2007 if key in parameters.get(section, {}):
2008 value = parameters[section].pop(key)
2009 if key in parameters:
2010 value = parameters.pop(key)
2011 if key in kwargs:
2012 value = kwargs.pop(key)
2014 if value is not None:
2015 sec_list[key] = value
2017 # Check if there is a key(i) version (no extra parsing)
2018 for arg_key in list(parameters.get(section, {})):
2019 if arg_key.split('(')[0].strip().lower() == key.lower():
2020 sec_list[arg_key] = parameters[section].pop(arg_key)
2021 cp_parameters = parameters.copy()
2022 for arg_key in cp_parameters:
2023 if arg_key.split('(')[0].strip().lower() == key.lower():
2024 sec_list[arg_key] = parameters.pop(arg_key)
2025 cp_kwargs = kwargs.copy()
2026 for arg_key in cp_kwargs:
2027 if arg_key.split('(')[0].strip().lower() == key.lower():
2028 sec_list[arg_key] = kwargs.pop(arg_key)
2030 # Add to output
2031 input_namelist[section] = sec_list
2033 unused_keys = list(kwargs)
2034 # pass anything else already in a section
2035 for key, value in parameters.items():
2036 if key in keys and isinstance(value, dict):
2037 input_namelist[key].update(value)
2038 elif isinstance(value, dict):
2039 unused_keys.extend(list(value))
2040 else:
2041 unused_keys.append(key)
2043 if warn and unused_keys:
2044 warnings.warn('Unused keys: {}'.format(', '.join(unused_keys)))
2046 return input_namelist
2049@deprecated('Please use the .to_string() method of Namelist instead.',
2050 DeprecationWarning)
2051def namelist_to_string(input_parameters):
2052 """Format a Namelist object as a string for writing to a file.
2053 Assume sections are ordered (taken care of in namelist construction)
2054 and that repr converts to a QE readable representation (except bools)
2056 .. deprecated:: 3.23.0
2057 Please use the :meth:`ase.io.espresso.Namelist.to_string` method
2058 instead.
2060 Parameters
2061 ----------
2062 input_parameters : Namelist | dict
2063 Expecting a nested dictionary of sections and key-value data.
2065 Returns
2066 -------
2067 pwi : List[str]
2068 Input line for the namelist
2069 """
2070 pwi = []
2071 for section in input_parameters:
2072 pwi.append(f'&{section.upper()}\n')
2073 for key, value in input_parameters[section].items():
2074 if value is True:
2075 pwi.append(f' {key:16} = .true.\n')
2076 elif value is False:
2077 pwi.append(f' {key:16} = .false.\n')
2078 elif isinstance(value, Path):
2079 pwi.append(f' {key:16} = "{value}"\n')
2080 else:
2081 # repr format to get quotes around strings
2082 pwi.append(f' {key:16} = {value!r}\n')
2083 pwi.append('/\n') # terminate section
2084 pwi.append('\n')
2085 return pwi