Coverage for ase / io / espresso.py: 77.59%
870 statements
« prev ^ index » next coverage.py v7.13.3, created at 2026-02-04 10:20 +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)
1454@writer
1455def write_espresso_ph(
1456 fd,
1457 input_data=None,
1458 qpts=None,
1459 nat_todo_indices=None,
1460 **kwargs) -> None:
1461 """
1462 Function that write the input file for a ph.x calculation. Normal namelist
1463 cards are passed in the input_data dictionary. Which can be either nested
1464 or flat, ASE style. The q-points are passed in the qpts list. If qplot is
1465 set to True then qpts is expected to be a list of list|tuple of length 4.
1466 Where the first three elements are the coordinates of the q-point in units
1467 of 2pi/alat and the last element is the weight of the q-point. if qplot is
1468 set to False then qpts is expected to be a simple list of length 4 (single
1469 q-point). Finally if ldisp is set to True, the above is discarded and the
1470 q-points are read from the nq1, nq2, nq3 cards in the input_data dictionary.
1472 Additionally, a nat_todo_indices kwargs (list[int]) can be specified in the
1473 kwargs. It will be used if nat_todo is set to True in the input_data
1474 dictionary.
1476 Globally, this function follows the convention set in the ph.x documentation
1477 (https://www.quantum-espresso.org/Doc/INPUT_PH.html)
1479 Parameters
1480 ----------
1481 fd
1482 The file descriptor of the input file.
1484 kwargs
1485 kwargs dictionary possibly containing the following keys:
1487 - input_data: dict
1488 - qpts: list[list[float]] | list[tuple[float]] | list[float]
1489 - nat_todo_indices: list[int]
1491 Returns
1492 -------
1493 None
1494 """
1496 input_data = Namelist(input_data)
1497 input_data.to_nested('ph', **kwargs)
1499 input_ph = input_data["inputph"]
1501 inp_nat_todo = input_ph.get("nat_todo", 0)
1502 qpts = qpts or (0, 0, 0)
1504 pwi = input_data.to_string()
1506 fd.write(pwi)
1508 qplot = input_ph.get("qplot", False)
1509 ldisp = input_ph.get("ldisp", False)
1511 if qplot:
1512 fd.write(f"{len(qpts)}\n")
1513 for qpt in qpts:
1514 fd.write(
1515 f"{qpt[0]:0.8f} {qpt[1]:0.8f} {qpt[2]:0.8f} {qpt[3]:1d}\n"
1516 )
1517 elif not (qplot or ldisp):
1518 fd.write(f"{qpts[0]:0.8f} {qpts[1]:0.8f} {qpts[2]:0.8f}\n")
1519 if inp_nat_todo:
1520 tmp = [str(i) for i in nat_todo_indices]
1521 fd.write(" ".join(tmp))
1522 fd.write("\n")
1525@reader
1526def read_espresso_ph(fileobj):
1527 """
1528 Function that reads the output of a ph.x calculation.
1529 It returns a dictionary where each q-point number is a key and
1530 the value is a dictionary with the following keys if available:
1532 - qpoints: The q-point in cartesian coordinates.
1533 - kpoints: The k-points in cartesian coordinates.
1534 - dieltensor: The dielectric tensor.
1535 - borneffcharge: The effective Born charges.
1536 - borneffcharge_dfpt: The effective Born charges from DFPT.
1537 - polarizability: The polarizability tensor.
1538 - modes: The phonon modes.
1539 - eqpoints: The symmetrically equivalent q-points.
1540 - freqs: The phonon frequencies.
1541 - mode_symmetries: The symmetries of the modes.
1542 - atoms: The atoms object.
1544 Some notes:
1546 - For some reason, the cell is not defined to high level of
1547 precision in ph.x outputs. Be careful when using the atoms object
1548 retrieved from this function.
1549 - This function can be called on incomplete calculations i.e.
1550 if the calculation couldn't diagonalize the dynamical matrix
1551 for some q-points, the results for the other q-points will
1552 still be returned.
1554 Parameters
1555 ----------
1556 fileobj
1557 The file descriptor of the output file.
1559 Returns
1560 -------
1561 dict
1562 The results dictionnary as described above.
1563 """
1564 freg = re.compile(r"-?(?:0|[1-9]\d*)(?:\.\d+)?(?:[eE][+\-]?\d+)?")
1566 QPOINTS = r"(?i)^\s*Calculation\s*of\s*q"
1567 NKPTS = r"(?i)^\s*number\s*of\s*k\s*points\s*"
1568 DIEL = r"(?i)^\s*Dielectric\s*constant\s*in\s*cartesian\s*axis\s*$"
1569 BORN = r"(?i)^\s*Effective\s*charges\s*\(d\s*Force\s*/\s*dE\)"
1570 POLA = r"(?i)^\s*Polarizability\s*(a.u.)\^3"
1571 REPR = r"(?i)^\s*There\s*are\s*\d+\s*irreducible\s*representations\s*$"
1572 EQPOINTS = r"(?i)^\s*Number\s*of\s*q\s*in\s*the\s*star\s*=\s*"
1573 DIAG = r"(?i)^\s*Diagonalizing\s*the\s*dynamical\s*matrix\s*$"
1574 MODE_SYM = r"(?i)^\s*Mode\s*symmetry,\s*"
1575 BORN_DFPT = r"(?i)^\s*Effective\s*charges\s*\(d\s*P\s*/\s*du\)"
1576 POSITIONS = r"(?i)^\s*site\s*n\..*\(alat\s*units\)"
1577 ALAT = r"(?i)^\s*celldm\(1\)="
1578 CELL = (
1579 r"^\s*crystal\s*axes:\s*\(cart.\s*coord.\s*in\s*units\s*of\s*alat\)"
1580 )
1581 ELECTRON_PHONON = r"(?i)^\s*electron-phonon\s*interaction\s*...\s*$"
1583 output = {
1584 QPOINTS: [],
1585 NKPTS: [],
1586 DIEL: [],
1587 BORN: [],
1588 BORN_DFPT: [],
1589 POLA: [],
1590 REPR: [],
1591 EQPOINTS: [],
1592 DIAG: [],
1593 MODE_SYM: [],
1594 POSITIONS: [],
1595 ALAT: [],
1596 CELL: [],
1597 ELECTRON_PHONON: [],
1598 }
1600 names = {
1601 QPOINTS: "qpoints",
1602 NKPTS: "kpoints",
1603 DIEL: "dieltensor",
1604 BORN: "borneffcharge",
1605 BORN_DFPT: "borneffcharge_dfpt",
1606 POLA: "polarizability",
1607 REPR: "representations",
1608 EQPOINTS: "eqpoints",
1609 DIAG: "freqs",
1610 MODE_SYM: "mode_symmetries",
1611 POSITIONS: "positions",
1612 ALAT: "alat",
1613 CELL: "cell",
1614 ELECTRON_PHONON: "ep_data",
1615 }
1617 unique = {
1618 QPOINTS: True,
1619 NKPTS: False,
1620 DIEL: True,
1621 BORN: True,
1622 BORN_DFPT: True,
1623 POLA: True,
1624 REPR: True,
1625 EQPOINTS: True,
1626 DIAG: True,
1627 MODE_SYM: True,
1628 POSITIONS: True,
1629 ALAT: True,
1630 CELL: True,
1631 ELECTRON_PHONON: True,
1632 }
1634 results = {}
1635 fdo_lines = [i for i in fileobj.read().splitlines() if i]
1636 n_lines = len(fdo_lines)
1638 for idx, line in enumerate(fdo_lines):
1639 for key in output:
1640 if bool(re.match(key, line)):
1641 output[key].append(idx)
1643 output = {key: np.array(value) for key, value in output.items()}
1645 def _read_qpoints(idx):
1646 match = re.findall(freg, fdo_lines[idx])
1647 return tuple(round(float(x), 7) for x in match)
1649 def _read_kpoints(idx):
1650 n_kpts = int(re.findall(freg, fdo_lines[idx])[0])
1651 kpts = []
1652 for line in fdo_lines[idx + 2: idx + 2 + n_kpts]:
1653 if bool(re.search(r"^\s*k\(.*wk", line)):
1654 kpts.append([round(float(x), 7)
1655 for x in re.findall(freg, line)[1:]])
1656 return np.array(kpts)
1658 def _read_repr(idx):
1659 n_repr, curr, n = int(re.findall(freg, fdo_lines[idx])[0]), 0, 0
1660 representations = {}
1661 while idx + n < n_lines:
1662 if re.search(r"^\s*Representation.*modes", fdo_lines[idx + n]):
1663 curr = int(re.findall(freg, fdo_lines[idx + n])[0])
1664 representations[curr] = {"done": False, "modes": []}
1665 if re.search(r"Calculated\s*using\s*symmetry", fdo_lines[idx + n]) \
1666 or re.search(r"-\s*Done\s*$", fdo_lines[idx + n]):
1667 representations[curr]["done"] = True
1668 if re.search(r"(?i)^\s*(mode\s*#\s*\d\s*)+", fdo_lines[idx + n]):
1669 representations[curr]["modes"] = _read_modes(idx + n)
1670 if curr == n_repr:
1671 break
1672 n += 1
1673 return representations
1675 def _read_modes(idx):
1676 n = 1
1677 n_modes = len(re.findall(r"mode", fdo_lines[idx]))
1678 modes = []
1679 while not modes or bool(re.match(r"^\s*\(", fdo_lines[idx + n])):
1680 tmp = re.findall(freg, fdo_lines[idx + n])
1681 modes.append([round(float(x), 5) for x in tmp])
1682 n += 1
1683 return np.hsplit(np.array(modes), n_modes)
1685 def _read_eqpoints(idx):
1686 n_star = int(re.findall(freg, fdo_lines[idx])[0])
1687 return np.loadtxt(
1688 fdo_lines[idx + 2: idx + 2 + n_star], usecols=(1, 2, 3)
1689 ).reshape(-1, 3)
1691 def _read_freqs(idx):
1692 n = 0
1693 freqs = []
1694 stop = 0
1695 while not freqs or stop < 2:
1696 if bool(re.search(r"^\s*freq", fdo_lines[idx + n])):
1697 tmp = re.findall(freg, fdo_lines[idx + n])[1]
1698 freqs.append(float(tmp))
1699 if bool(re.search(r"\*{5,}", fdo_lines[idx + n])):
1700 stop += 1
1701 n += 1
1702 return np.array(freqs)
1704 def _read_sym(idx):
1705 n = 1
1706 sym = {}
1707 while bool(re.match(r"^\s*freq", fdo_lines[idx + n])):
1708 r = re.findall("\\d+", fdo_lines[idx + n])
1709 r = tuple(range(int(r[0]), int(r[1]) + 1))
1710 sym[r] = fdo_lines[idx + n].split("-->")[1].strip()
1711 sym[r] = re.sub(r"\s+", " ", sym[r])
1712 n += 1
1713 return sym
1715 def _read_epsil(idx):
1716 epsil = np.zeros((3, 3))
1717 for n in range(1, 4):
1718 tmp = re.findall(freg, fdo_lines[idx + n])
1719 epsil[n - 1] = [round(float(x), 9) for x in tmp]
1720 return epsil
1722 def _read_born(idx):
1723 n = 1
1724 born = []
1725 while idx + n < n_lines:
1726 if re.search(r"^\s*atom\s*\d\s*\S", fdo_lines[idx + n]):
1727 pass
1728 elif re.search(r"^\s*E\*?(x|y|z)\s*\(", fdo_lines[idx + n]):
1729 tmp = re.findall(freg, fdo_lines[idx + n])
1730 born.append([round(float(x), 5) for x in tmp])
1731 else:
1732 break
1733 n += 1
1734 born = np.array(born)
1735 return np.vsplit(born, len(born) // 3)
1737 def _read_born_dfpt(idx):
1738 n = 1
1739 born = []
1740 while idx + n < n_lines:
1741 if re.search(r"^\s*atom\s*\d\s*\S", fdo_lines[idx + n]):
1742 pass
1743 elif re.search(r"^\s*P(x|y|z)\s*\(", fdo_lines[idx + n]):
1744 tmp = re.findall(freg, fdo_lines[idx + n])
1745 born.append([round(float(x), 5) for x in tmp])
1746 else:
1747 break
1748 n += 1
1749 born = np.array(born)
1750 return np.vsplit(born, len(born) // 3)
1752 def _read_pola(idx):
1753 pola = np.zeros((3, 3))
1754 for n in range(1, 4):
1755 tmp = re.findall(freg, fdo_lines[idx + n])[:3]
1756 pola[n - 1] = [round(float(x), 2) for x in tmp]
1757 return pola
1759 def _read_positions(idx):
1760 positions = []
1761 symbols = []
1762 n = 1
1763 while re.findall(r"^\s*\d+", fdo_lines[idx + n]):
1764 symbols.append(fdo_lines[idx + n].split()[1])
1765 positions.append(
1766 [round(float(x), 5)
1767 for x in re.findall(freg, fdo_lines[idx + n])[-3:]]
1768 )
1769 n += 1
1770 atoms = Atoms(positions=positions, symbols=symbols)
1771 atoms.pbc = True
1772 return atoms
1774 def _read_alat(idx):
1775 return round(float(re.findall(freg, fdo_lines[idx])[1]), 5)
1777 def _read_cell(idx):
1778 cell = []
1779 n = 1
1780 while re.findall(r"^\s*a\(\d\)", fdo_lines[idx + n]):
1781 cell.append([round(float(x), 4)
1782 for x in re.findall(freg, fdo_lines[idx + n])[-3:]])
1783 n += 1
1784 return np.array(cell)
1786 def _read_electron_phonon(idx):
1787 results = {}
1789 broad_re = (
1790 r"^\s*Gaussian\s*Broadening:\s+([\d.]+)\s+Ry, ngauss=\s+\d+"
1791 )
1792 dos_re = (
1793 r"^\s*DOS\s*=\s*([\d.]+)\s*states/"
1794 r"spin/Ry/Unit\s*Cell\s*at\s*Ef=\s+([\d.]+)\s+eV"
1795 )
1796 lg_re = (
1797 r"^\s*lambda\(\s+(\d+)\)=\s+([\d.]+)\s+gamma=\s+([\d.]+)\s+GHz"
1798 )
1799 end_re = r"^\s*Number\s*of\s*q\s*in\s*the\s*star\s*=\s+(\d+)$"
1801 lambdas = []
1802 gammas = []
1804 current = None
1806 n = 1
1807 while idx + n < n_lines:
1808 line = fdo_lines[idx + n]
1810 broad_match = re.match(broad_re, line)
1811 dos_match = re.match(dos_re, line)
1812 lg_match = re.match(lg_re, line)
1813 end_match = re.match(end_re, line)
1815 if broad_match:
1816 if lambdas:
1817 results[current]["lambdas"] = lambdas
1818 results[current]["gammas"] = gammas
1819 lambdas = []
1820 gammas = []
1821 current = float(broad_match[1])
1822 results[current] = {}
1823 elif dos_match:
1824 results[current]["dos"] = float(dos_match[1])
1825 results[current]["fermi"] = float(dos_match[2])
1826 elif lg_match:
1827 lambdas.append(float(lg_match[2]))
1828 gammas.append(float(lg_match[3]))
1830 if end_match:
1831 results[current]["lambdas"] = lambdas
1832 results[current]["gammas"] = gammas
1833 break
1835 n += 1
1837 return results
1839 properties = {
1840 NKPTS: _read_kpoints,
1841 DIEL: _read_epsil,
1842 BORN: _read_born,
1843 BORN_DFPT: _read_born_dfpt,
1844 POLA: _read_pola,
1845 REPR: _read_repr,
1846 EQPOINTS: _read_eqpoints,
1847 DIAG: _read_freqs,
1848 MODE_SYM: _read_sym,
1849 POSITIONS: _read_positions,
1850 ALAT: _read_alat,
1851 CELL: _read_cell,
1852 ELECTRON_PHONON: _read_electron_phonon,
1853 }
1855 iblocks = np.append(output[QPOINTS], n_lines)
1857 for qnum, (past, future) in enumerate(zip(iblocks[:-1], iblocks[1:])):
1858 qpoint = _read_qpoints(past)
1859 results[qnum + 1] = curr_result = {"qpoint": qpoint}
1860 for prop in properties:
1861 p = (past < output[prop]) & (output[prop] < future)
1862 selected = output[prop][p]
1863 if len(selected) == 0:
1864 continue
1865 if unique[prop]:
1866 idx = output[prop][p][-1]
1867 curr_result[names[prop]] = properties[prop](idx)
1868 else:
1869 tmp = {k + 1: 0 for k in range(len(selected))}
1870 for k, idx in enumerate(selected):
1871 tmp[k + 1] = properties[prop](idx)
1872 curr_result[names[prop]] = tmp
1873 alat = curr_result.pop("alat", 1.0)
1874 atoms = curr_result.pop("positions", None)
1875 cell = curr_result.pop("cell", np.eye(3))
1876 if atoms:
1877 atoms.positions *= alat * units["Bohr"]
1878 atoms.cell = cell * alat * units["Bohr"]
1879 atoms.wrap()
1880 curr_result["atoms"] = atoms
1882 return results
1885def write_fortran_namelist(
1886 fd,
1887 input_data=None,
1888 binary=None,
1889 additional_cards=None,
1890 **kwargs) -> None:
1891 """
1892 Function which writes input for simple espresso binaries.
1893 List of supported binaries are in the espresso_keys.py file.
1894 Non-exhaustive list (to complete)
1896 Note: "EOF" is appended at the end of the file.
1897 (https://lists.quantum-espresso.org/pipermail/users/2020-November/046269.html)
1899 Additional fields are written 'as is' in the input file. It is expected
1900 to be a string or a list of strings.
1902 Parameters
1903 ----------
1904 fd
1905 The file descriptor of the input file.
1906 input_data: dict
1907 A flat or nested dictionary with input parameters for the binary.x
1908 binary: str
1909 Name of the binary
1910 additional_cards: str | list[str]
1911 Additional fields to be written at the end of the input file, after
1912 the namelist. It is expected to be a string or a list of strings.
1914 Returns
1915 -------
1916 None
1917 """
1918 input_data = Namelist(input_data)
1920 if binary:
1921 input_data.to_nested(binary, **kwargs)
1923 pwi = input_data.to_string()
1925 fd.write(pwi)
1927 if additional_cards:
1928 if isinstance(additional_cards, list):
1929 additional_cards = "\n".join(additional_cards)
1930 additional_cards += "\n"
1932 fd.write(additional_cards)
1934 fd.write("EOF")
1937@deprecated('Please use the ase.io.espresso.Namelist class',
1938 DeprecationWarning)
1939def construct_namelist(parameters=None, keys=None, warn=False, **kwargs):
1940 """
1941 Construct an ordered Namelist containing all the parameters given (as
1942 a dictionary or kwargs). Keys will be inserted into their appropriate
1943 section in the namelist and the dictionary may contain flat and nested
1944 structures. Any kwargs that match input keys will be incorporated into
1945 their correct section. All matches are case-insensitive, and returned
1946 Namelist object is a case-insensitive dict.
1948 If a key is not known to ase, but in a section within `parameters`,
1949 it will be assumed that it was put there on purpose and included
1950 in the output namelist. Anything not in a section will be ignored (set
1951 `warn` to True to see ignored keys).
1953 Keys with a dimension (e.g. Hubbard_U(1)) will be incorporated as-is
1954 so the `i` should be made to match the output.
1956 The priority of the keys is:
1957 kwargs[key] > parameters[key] > parameters[section][key]
1958 Only the highest priority item will be included.
1960 .. deprecated:: 3.23.0
1961 Please use :class:`ase.io.espresso.Namelist` instead.
1963 Parameters
1964 ----------
1965 parameters: dict
1966 Flat or nested set of input parameters.
1967 keys: Namelist | dict
1968 Namelist to use as a template for the output.
1969 warn: bool
1970 Enable warnings for unused keys.
1972 Returns
1973 -------
1974 input_namelist: Namelist
1975 pw.x compatible namelist of input parameters.
1977 """
1979 if keys is None:
1980 keys = deepcopy(pw_keys)
1981 # Convert everything to Namelist early to make case-insensitive
1982 if parameters is None:
1983 parameters = Namelist()
1984 else:
1985 # Maximum one level of nested dict
1986 # Don't modify in place
1987 parameters_namelist = Namelist()
1988 for key, value in parameters.items():
1989 if isinstance(value, dict):
1990 parameters_namelist[key] = Namelist(value)
1991 else:
1992 parameters_namelist[key] = value
1993 parameters = parameters_namelist
1995 # Just a dict
1996 kwargs = Namelist(kwargs)
1998 # Final parameter set
1999 input_namelist = Namelist()
2001 # Collect
2002 for section in keys:
2003 sec_list = Namelist()
2004 for key in keys[section]:
2005 # Check all three separately and pop them all so that
2006 # we can check for missing values later
2007 value = None
2009 if key in parameters.get(section, {}):
2010 value = parameters[section].pop(key)
2011 if key in parameters:
2012 value = parameters.pop(key)
2013 if key in kwargs:
2014 value = kwargs.pop(key)
2016 if value is not None:
2017 sec_list[key] = value
2019 # Check if there is a key(i) version (no extra parsing)
2020 for arg_key in list(parameters.get(section, {})):
2021 if arg_key.split('(')[0].strip().lower() == key.lower():
2022 sec_list[arg_key] = parameters[section].pop(arg_key)
2023 cp_parameters = parameters.copy()
2024 for arg_key in cp_parameters:
2025 if arg_key.split('(')[0].strip().lower() == key.lower():
2026 sec_list[arg_key] = parameters.pop(arg_key)
2027 cp_kwargs = kwargs.copy()
2028 for arg_key in cp_kwargs:
2029 if arg_key.split('(')[0].strip().lower() == key.lower():
2030 sec_list[arg_key] = kwargs.pop(arg_key)
2032 # Add to output
2033 input_namelist[section] = sec_list
2035 unused_keys = list(kwargs)
2036 # pass anything else already in a section
2037 for key, value in parameters.items():
2038 if key in keys and isinstance(value, dict):
2039 input_namelist[key].update(value)
2040 elif isinstance(value, dict):
2041 unused_keys.extend(list(value))
2042 else:
2043 unused_keys.append(key)
2045 if warn and unused_keys:
2046 warnings.warn('Unused keys: {}'.format(', '.join(unused_keys)))
2048 return input_namelist
2051@deprecated('Please use the .to_string() method of Namelist instead.',
2052 DeprecationWarning)
2053def namelist_to_string(input_parameters):
2054 """Format a Namelist object as a string for writing to a file.
2055 Assume sections are ordered (taken care of in namelist construction)
2056 and that repr converts to a QE readable representation (except bools)
2058 .. deprecated:: 3.23.0
2059 Please use the :meth:`ase.io.espresso.Namelist.to_string` method
2060 instead.
2062 Parameters
2063 ----------
2064 input_parameters : Namelist | dict
2065 Expecting a nested dictionary of sections and key-value data.
2067 Returns
2068 -------
2069 pwi : List[str]
2070 Input line for the namelist
2071 """
2072 pwi = []
2073 for section in input_parameters:
2074 pwi.append(f'&{section.upper()}\n')
2075 for key, value in input_parameters[section].items():
2076 if value is True:
2077 pwi.append(f' {key:16} = .true.\n')
2078 elif value is False:
2079 pwi.append(f' {key:16} = .false.\n')
2080 elif isinstance(value, Path):
2081 pwi.append(f' {key:16} = "{value}"\n')
2082 else:
2083 # repr format to get quotes around strings
2084 pwi.append(f' {key:16} = {value!r}\n')
2085 pwi.append('/\n') # terminate section
2086 pwi.append('\n')
2087 return pwi