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