Coverage for /builds/kinetik161/ase/ase/io/espresso.py: 83.08%
597 statements
« prev ^ index » next coverage.py v7.2.7, created at 2023-12-10 11:04 +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 OrderedDict
19from copy import deepcopy
21import numpy as np
22from ase.atoms import Atoms
23from ase.calculators.calculator import kpts2ndarray, kpts2sizeandoffsets
24from ase.calculators.singlepoint import (SinglePointDFTCalculator,
25 SinglePointKPoint)
26from ase.constraints import FixAtoms, FixCartesian
27from ase.data import chemical_symbols
28from ase.dft.kpoints import kpoint_convert
29from ase.units import create_units
30from ase.utils import iofunction
32# Quantum ESPRESSO uses CODATA 2006 internally
33units = create_units('2006')
35# Section identifiers
36_PW_START = 'Program PWSCF'
37_PW_END = 'End of self-consistent calculation'
38_PW_CELL = 'CELL_PARAMETERS'
39_PW_POS = 'ATOMIC_POSITIONS'
40_PW_MAGMOM = 'Magnetic moment per site'
41_PW_FORCE = 'Forces acting on atoms'
42_PW_TOTEN = '! total energy'
43_PW_STRESS = 'total stress'
44_PW_FERMI = 'the Fermi energy is'
45_PW_HIGHEST_OCCUPIED = 'highest occupied level'
46_PW_HIGHEST_OCCUPIED_LOWEST_FREE = 'highest occupied, lowest unoccupied level'
47_PW_KPTS = 'number of k points='
48_PW_BANDS = _PW_END
49_PW_BANDSTRUCTURE = 'End of band structure calculation'
50_PW_DIPOLE = "Debye"
51_PW_DIPOLE_DIRECTION = "Computed dipole along edir"
53# ibrav error message
54ibrav_error_message = (
55 'ASE does not support ibrav != 0. Note that with ibrav '
56 '== 0, Quantum ESPRESSO will still detect the symmetries '
57 'of your system because the CELL_PARAMETERS are defined '
58 'to a high level of precision.')
61class Namelist(OrderedDict):
62 """Case insensitive dict that emulates Fortran Namelists."""
64 def __contains__(self, key):
65 return super().__contains__(key.lower())
67 def __delitem__(self, key):
68 return super().__delitem__(key.lower())
70 def __getitem__(self, key):
71 return super().__getitem__(key.lower())
73 def __setitem__(self, key, value):
74 super().__setitem__(key.lower(), value)
76 def get(self, key, default=None):
77 return super().get(key.lower(), default)
80@iofunction('r')
81def read_espresso_out(fileobj, index=slice(None), results_required=True):
82 """Reads Quantum ESPRESSO output files.
84 The atomistic configurations as well as results (energy, force, stress,
85 magnetic moments) of the calculation are read for all configurations
86 within the output file.
88 Will probably raise errors for broken or incomplete files.
90 Parameters
91 ----------
92 fileobj : file|str
93 A file like object or filename
94 index : slice
95 The index of configurations to extract.
96 results_required : bool
97 If True, atomistic configurations that do not have any
98 associated results will not be included. This prevents double
99 printed configurations and incomplete calculations from being
100 returned as the final configuration with no results data.
102 Yields
103 ------
104 structure : Atoms
105 The next structure from the index slice. The Atoms has a
106 SinglePointCalculator attached with any results parsed from
107 the file.
110 """
111 # work with a copy in memory for faster random access
112 pwo_lines = fileobj.readlines()
114 # TODO: index -1 special case?
115 # Index all the interesting points
116 indexes = {
117 _PW_START: [],
118 _PW_END: [],
119 _PW_CELL: [],
120 _PW_POS: [],
121 _PW_MAGMOM: [],
122 _PW_FORCE: [],
123 _PW_TOTEN: [],
124 _PW_STRESS: [],
125 _PW_FERMI: [],
126 _PW_HIGHEST_OCCUPIED: [],
127 _PW_HIGHEST_OCCUPIED_LOWEST_FREE: [],
128 _PW_KPTS: [],
129 _PW_BANDS: [],
130 _PW_BANDSTRUCTURE: [],
131 _PW_DIPOLE: [],
132 _PW_DIPOLE_DIRECTION: [],
133 }
135 for idx, line in enumerate(pwo_lines):
136 for identifier in indexes:
137 if identifier in line:
138 indexes[identifier].append(idx)
140 # Configurations are either at the start, or defined in ATOMIC_POSITIONS
141 # in a subsequent step. Can deal with concatenated output files.
142 all_config_indexes = sorted(indexes[_PW_START] +
143 indexes[_PW_POS])
145 # Slice only requested indexes
146 # setting results_required argument stops configuration-only
147 # structures from being returned. This ensures the [-1] structure
148 # is one that has results. Two cases:
149 # - SCF of last configuration is not converged, job terminated
150 # abnormally.
151 # - 'relax' and 'vc-relax' re-prints the final configuration but
152 # only 'vc-relax' recalculates.
153 if results_required:
154 results_indexes = sorted(indexes[_PW_TOTEN] + indexes[_PW_FORCE] +
155 indexes[_PW_STRESS] + indexes[_PW_MAGMOM] +
156 indexes[_PW_BANDS] +
157 indexes[_PW_BANDSTRUCTURE])
159 # Prune to only configurations with results data before the next
160 # configuration
161 results_config_indexes = []
162 for config_index, config_index_next in zip(
163 all_config_indexes,
164 all_config_indexes[1:] + [len(pwo_lines)]):
165 if any(config_index < results_index < config_index_next
166 for results_index in results_indexes):
167 results_config_indexes.append(config_index)
169 # slice from the subset
170 image_indexes = results_config_indexes[index]
171 else:
172 image_indexes = all_config_indexes[index]
174 # Extract initialisation information each time PWSCF starts
175 # to add to subsequent configurations. Use None so slices know
176 # when to fill in the blanks.
177 pwscf_start_info = {idx: None for idx in indexes[_PW_START]}
179 for image_index in image_indexes:
180 # Find the nearest calculation start to parse info. Needed in,
181 # for example, relaxation where cell is only printed at the
182 # start.
183 if image_index in indexes[_PW_START]:
184 prev_start_index = image_index
185 else:
186 # The greatest start index before this structure
187 prev_start_index = [idx for idx in indexes[_PW_START]
188 if idx < image_index][-1]
190 # add structure to reference if not there
191 if pwscf_start_info[prev_start_index] is None:
192 pwscf_start_info[prev_start_index] = parse_pwo_start(
193 pwo_lines, prev_start_index)
195 # Get the bounds for information for this structure. Any associated
196 # values will be between the image_index and the following one,
197 # EXCEPT for cell, which will be 4 lines before if it exists.
198 for next_index in all_config_indexes:
199 if next_index > image_index:
200 break
201 else:
202 # right to the end of the file
203 next_index = len(pwo_lines)
205 # Get the structure
206 # Use this for any missing data
207 prev_structure = pwscf_start_info[prev_start_index]['atoms']
208 if image_index in indexes[_PW_START]:
209 structure = prev_structure.copy() # parsed from start info
210 else:
211 if _PW_CELL in pwo_lines[image_index - 5]:
212 # CELL_PARAMETERS would be just before positions if present
213 cell, cell_alat = get_cell_parameters(
214 pwo_lines[image_index - 5:image_index])
215 else:
216 cell = prev_structure.cell
217 cell_alat = pwscf_start_info[prev_start_index]['alat']
219 # give at least enough lines to parse the positions
220 # should be same format as input card
221 n_atoms = len(prev_structure)
222 positions_card = get_atomic_positions(
223 pwo_lines[image_index:image_index + n_atoms + 1],
224 n_atoms=n_atoms, cell=cell, alat=cell_alat)
226 # convert to Atoms object
227 symbols = [label_to_symbol(position[0]) for position in
228 positions_card]
229 positions = [position[1] for position in positions_card]
230 structure = Atoms(symbols=symbols, positions=positions, cell=cell,
231 pbc=True)
233 # Extract calculation results
234 # Energy
235 energy = None
236 for energy_index in indexes[_PW_TOTEN]:
237 if image_index < energy_index < next_index:
238 energy = float(
239 pwo_lines[energy_index].split()[-2]) * units['Ry']
241 # Forces
242 forces = None
243 for force_index in indexes[_PW_FORCE]:
244 if image_index < force_index < next_index:
245 # Before QE 5.3 'negative rho' added 2 lines before forces
246 # Use exact lines to stop before 'non-local' forces
247 # in high verbosity
248 if not pwo_lines[force_index + 2].strip():
249 force_index += 4
250 else:
251 force_index += 2
252 # assume contiguous
253 forces = [
254 [float(x) for x in force_line.split()[-3:]] for force_line
255 in pwo_lines[force_index:force_index + len(structure)]]
256 forces = np.array(forces) * units['Ry'] / units['Bohr']
258 # Stress
259 stress = None
260 for stress_index in indexes[_PW_STRESS]:
261 if image_index < stress_index < next_index:
262 sxx, sxy, sxz = pwo_lines[stress_index + 1].split()[:3]
263 _, syy, syz = pwo_lines[stress_index + 2].split()[:3]
264 _, _, szz = pwo_lines[stress_index + 3].split()[:3]
265 stress = np.array([sxx, syy, szz, syz, sxz, sxy], dtype=float)
266 # sign convention is opposite of ase
267 stress *= -1 * units['Ry'] / (units['Bohr'] ** 3)
269 # Magmoms
270 magmoms = None
271 for magmoms_index in indexes[_PW_MAGMOM]:
272 if image_index < magmoms_index < next_index:
273 magmoms = [
274 float(mag_line.split()[-1]) for mag_line
275 in pwo_lines[magmoms_index + 1:
276 magmoms_index + 1 + len(structure)]]
278 # Dipole moment
279 dipole = None
280 if indexes[_PW_DIPOLE]:
281 for dipole_index in indexes[_PW_DIPOLE]:
282 if image_index < dipole_index < next_index:
283 _dipole = float(pwo_lines[dipole_index].split()[-2])
285 for dipole_index in indexes[_PW_DIPOLE_DIRECTION]:
286 if image_index < dipole_index < next_index:
287 _direction = pwo_lines[dipole_index].strip()
288 prefix = 'Computed dipole along edir('
289 _direction = _direction[len(prefix):]
290 _direction = int(_direction[0])
292 dipole = np.eye(3)[_direction - 1] * _dipole * units['Debye']
294 # Fermi level / highest occupied level
295 efermi = None
296 for fermi_index in indexes[_PW_FERMI]:
297 if image_index < fermi_index < next_index:
298 efermi = float(pwo_lines[fermi_index].split()[-2])
300 if efermi is None:
301 for ho_index in indexes[_PW_HIGHEST_OCCUPIED]:
302 if image_index < ho_index < next_index:
303 efermi = float(pwo_lines[ho_index].split()[-1])
305 if efermi is None:
306 for holf_index in indexes[_PW_HIGHEST_OCCUPIED_LOWEST_FREE]:
307 if image_index < holf_index < next_index:
308 efermi = float(pwo_lines[holf_index].split()[-2])
310 # K-points
311 ibzkpts = None
312 weights = None
313 kpoints_warning = "Number of k-points >= 100: " + \
314 "set verbosity='high' to print them."
316 for kpts_index in indexes[_PW_KPTS]:
317 nkpts = int(re.findall(r'\b\d+\b', pwo_lines[kpts_index])[0])
318 kpts_index += 2
320 if pwo_lines[kpts_index].strip() == kpoints_warning:
321 continue
323 # QE prints the k-points in units of 2*pi/alat
324 # with alat defined as the length of the first
325 # cell vector
326 cell = structure.get_cell()
327 alat = np.linalg.norm(cell[0])
328 ibzkpts = []
329 weights = []
330 for i in range(nkpts):
331 L = pwo_lines[kpts_index + i].split()
332 weights.append(float(L[-1]))
333 coord = np.array([L[-6], L[-5], L[-4].strip('),')],
334 dtype=float)
335 coord *= 2 * np.pi / alat
336 coord = kpoint_convert(cell, ckpts_kv=coord)
337 ibzkpts.append(coord)
338 ibzkpts = np.array(ibzkpts)
339 weights = np.array(weights)
341 # Bands
342 kpts = None
343 kpoints_warning = "Number of k-points >= 100: " + \
344 "set verbosity='high' to print the bands."
346 for bands_index in indexes[_PW_BANDS] + indexes[_PW_BANDSTRUCTURE]:
347 if image_index < bands_index < next_index:
348 bands_index += 1
349 # skip over the lines with DFT+U occupation matrices
350 if 'enter write_ns' in pwo_lines[bands_index]:
351 while 'exit write_ns' not in pwo_lines[bands_index]:
352 bands_index += 1
353 bands_index += 1
355 if pwo_lines[bands_index].strip() == kpoints_warning:
356 continue
358 assert ibzkpts is not None
359 spin, bands, eigenvalues = 0, [], [[], []]
361 while True:
362 L = pwo_lines[bands_index].replace('-', ' -').split()
363 if len(L) == 0:
364 if len(bands) > 0:
365 eigenvalues[spin].append(bands)
366 bands = []
367 elif L == ['occupation', 'numbers']:
368 # Skip the lines with the occupation numbers
369 bands_index += len(eigenvalues[spin][0]) // 8 + 1
370 elif L[0] == 'k' and L[1].startswith('='):
371 pass
372 elif 'SPIN' in L:
373 if 'DOWN' in L:
374 spin += 1
375 else:
376 try:
377 bands.extend(map(float, L))
378 except ValueError:
379 break
380 bands_index += 1
382 if spin == 1:
383 assert len(eigenvalues[0]) == len(eigenvalues[1])
384 assert len(eigenvalues[0]) == len(ibzkpts), \
385 (np.shape(eigenvalues), len(ibzkpts))
387 kpts = []
388 for s in range(spin + 1):
389 for w, k, e in zip(weights, ibzkpts, eigenvalues[s]):
390 kpt = SinglePointKPoint(w, s, k, eps_n=e)
391 kpts.append(kpt)
393 # Put everything together
394 #
395 # In PW the forces are consistent with the "total energy"; that's why
396 # its value must be assigned to free_energy.
397 # PW doesn't compute the extrapolation of the energy to 0K smearing
398 # the closer thing to this is again the total energy that contains
399 # the correct (i.e. variational) form of the band energy is
400 # Eband = \int e N(e) de for e<Ef , where N(e) is the DOS
401 # This differs by the term (-TS) from the sum of KS eigenvalues:
402 # Eks = \sum wg(n,k) et(n,k)
403 # which is non variational. When a Fermi-Dirac function is used
404 # for a given T, the variational energy is REALLY the free energy F,
405 # and F = E - TS , with E = non variational energy.
406 #
407 calc = SinglePointDFTCalculator(structure, energy=energy,
408 free_energy=energy,
409 forces=forces, stress=stress,
410 magmoms=magmoms, efermi=efermi,
411 ibzkpts=ibzkpts, dipole=dipole)
412 calc.kpts = kpts
413 structure.calc = calc
415 yield structure
418def parse_pwo_start(lines, index=0):
419 """Parse Quantum ESPRESSO calculation info from lines,
420 starting from index. Return a dictionary containing extracted
421 information.
423 - `celldm(1)`: lattice parameters (alat)
424 - `cell`: unit cell in Angstrom
425 - `symbols`: element symbols for the structure
426 - `positions`: cartesian coordinates of atoms in Angstrom
427 - `atoms`: an `ase.Atoms` object constructed from the extracted data
429 Parameters
430 ----------
431 lines : list[str]
432 Contents of PWSCF output file.
433 index : int
434 Line number to begin parsing. Only first calculation will
435 be read.
437 Returns
438 -------
439 info : dict
440 Dictionary of calculation parameters, including `celldm(1)`, `cell`,
441 `symbols`, `positions`, `atoms`.
443 Raises
444 ------
445 KeyError
446 If interdependent values cannot be found (especially celldm(1))
447 an error will be raised as other quantities cannot then be
448 calculated (e.g. cell and positions).
449 """
450 # TODO: extend with extra DFT info?
452 info = {}
454 for idx, line in enumerate(lines[index:], start=index):
455 if 'celldm(1)' in line:
456 # celldm(1) has more digits than alat!!
457 info['celldm(1)'] = float(line.split()[1]) * units['Bohr']
458 info['alat'] = info['celldm(1)']
459 elif 'number of atoms/cell' in line:
460 info['nat'] = int(line.split()[-1])
461 elif 'number of atomic types' in line:
462 info['ntyp'] = int(line.split()[-1])
463 elif 'crystal axes:' in line:
464 info['cell'] = info['celldm(1)'] * np.array([
465 [float(x) for x in lines[idx + 1].split()[3:6]],
466 [float(x) for x in lines[idx + 2].split()[3:6]],
467 [float(x) for x in lines[idx + 3].split()[3:6]]])
468 elif 'positions (alat units)' in line:
469 info['symbols'], info['positions'] = [], []
471 for at_line in lines[idx + 1:idx + 1 + info['nat']]:
472 sym, x, y, z = parse_position_line(at_line)
473 info['symbols'].append(label_to_symbol(sym))
474 info['positions'].append([x * info['celldm(1)'],
475 y * info['celldm(1)'],
476 z * info['celldm(1)']])
477 # This should be the end of interesting info.
478 # Break here to avoid dealing with large lists of kpoints.
479 # Will need to be extended for DFTCalculator info.
480 break
482 # Make atoms for convenience
483 info['atoms'] = Atoms(symbols=info['symbols'],
484 positions=info['positions'],
485 cell=info['cell'], pbc=True)
487 return info
490def parse_position_line(line):
491 """Parse a single line from a pw.x output file.
493 The line must contain information about the atomic symbol and the position,
494 e.g.
496 995 Sb tau( 995) = ( 1.4212023 0.7037863 0.1242640 )
498 Parameters
499 ----------
500 line : str
501 Line to be parsed.
503 Returns
504 -------
505 sym : str
506 Atomic symbol.
507 x : float
508 x-position.
509 y : float
510 y-position.
511 z : float
512 z-position.
513 """
514 pat = re.compile(r'\s*\d+\s*(\S+)\s*tau\(\s*\d+\)\s*='
515 r'\s*\(\s*(\S+)\s+(\S+)\s+(\S+)\s*\)')
516 match = pat.match(line)
517 assert match is not None
518 sym, x, y, z = match.group(1, 2, 3, 4)
519 return sym, float(x), float(y), float(z)
522@iofunction('r')
523def read_espresso_in(fileobj):
524 """Parse a Quantum ESPRESSO input files, '.in', '.pwi'.
526 ESPRESSO inputs are generally a fortran-namelist format with custom
527 blocks of data. The namelist is parsed as a dict and an atoms object
528 is constructed from the included information.
530 Parameters
531 ----------
532 fileobj : file | str
533 A file-like object that supports line iteration with the contents
534 of the input file, or a filename.
536 Returns
537 -------
538 atoms : Atoms
539 Structure defined in the input file.
541 Raises
542 ------
543 KeyError
544 Raised for missing keys that are required to process the file
545 """
546 # parse namelist section and extract remaining lines
547 data, card_lines = read_fortran_namelist(fileobj)
549 # get the cell if ibrav=0
550 if 'system' not in data:
551 raise KeyError('Required section &SYSTEM not found.')
552 elif 'ibrav' not in data['system']:
553 raise KeyError('ibrav is required in &SYSTEM')
554 elif data['system']['ibrav'] == 0:
555 # celldm(1) is in Bohr, A is in angstrom. celldm(1) will be
556 # used even if A is also specified.
557 if 'celldm(1)' in data['system']:
558 alat = data['system']['celldm(1)'] * units['Bohr']
559 elif 'A' in data['system']:
560 alat = data['system']['A']
561 else:
562 alat = None
563 cell, _ = get_cell_parameters(card_lines, alat=alat)
564 else:
565 raise ValueError(ibrav_error_message)
567 # species_info holds some info for each element
568 species_card = get_atomic_species(
569 card_lines, n_species=data['system']['ntyp'])
570 species_info = {}
571 for ispec, (label, weight, pseudo) in enumerate(species_card):
572 symbol = label_to_symbol(label)
574 # starting_magnetization is in fractions of valence electrons
575 magnet_key = f"starting_magnetization({ispec + 1})"
576 magmom = data["system"].get(magnet_key, 0.0)
577 species_info[symbol] = {"weight": weight, "pseudo": pseudo,
578 "magmom": magmom}
580 positions_card = get_atomic_positions(
581 card_lines, n_atoms=data['system']['nat'], cell=cell, alat=alat)
583 symbols = [label_to_symbol(position[0]) for position in positions_card]
584 positions = [position[1] for position in positions_card]
585 magmoms = [species_info[symbol]["magmom"] for symbol in symbols]
587 # TODO: put more info into the atoms object
588 # e.g magmom, forces.
589 atoms = Atoms(symbols=symbols, positions=positions, cell=cell, pbc=True,
590 magmoms=magmoms)
592 return atoms
595def get_atomic_positions(lines, n_atoms, cell=None, alat=None):
596 """Parse atom positions from ATOMIC_POSITIONS card.
598 Parameters
599 ----------
600 lines : list[str]
601 A list of lines containing the ATOMIC_POSITIONS card.
602 n_atoms : int
603 Expected number of atoms. Only this many lines will be parsed.
604 cell : np.array
605 Unit cell of the crystal. Only used with crystal coordinates.
606 alat : float
607 Lattice parameter for atomic coordinates. Only used for alat case.
609 Returns
610 -------
611 positions : list[(str, (float, float, float), (float, float, float))]
612 A list of the ordered atomic positions in the format:
613 label, (x, y, z), (if_x, if_y, if_z)
614 Force multipliers are set to None if not present.
616 Raises
617 ------
618 ValueError
619 Any problems parsing the data result in ValueError
621 """
623 positions = None
624 # no blanks or comment lines, can the consume n_atoms lines for positions
625 trimmed_lines = (line for line in lines
626 if line.strip() and not line[0] == '#')
628 for line in trimmed_lines:
629 if line.strip().startswith('ATOMIC_POSITIONS'):
630 if positions is not None:
631 raise ValueError('Multiple ATOMIC_POSITIONS specified')
632 # Priority and behaviour tested with QE 5.3
633 if 'crystal_sg' in line.lower():
634 raise NotImplementedError('CRYSTAL_SG not implemented')
635 elif 'crystal' in line.lower():
636 cell = cell
637 elif 'bohr' in line.lower():
638 cell = np.identity(3) * units['Bohr']
639 elif 'angstrom' in line.lower():
640 cell = np.identity(3)
641 # elif 'alat' in line.lower():
642 # cell = np.identity(3) * alat
643 else:
644 if alat is None:
645 raise ValueError('Set lattice parameter in &SYSTEM for '
646 'alat coordinates')
647 # Always the default, will be DEPRECATED as mandatory
648 # in future
649 cell = np.identity(3) * alat
651 positions = []
652 for _dummy in range(n_atoms):
653 split_line = next(trimmed_lines).split()
654 # These can be fractions and other expressions
655 position = np.dot((infix_float(split_line[1]),
656 infix_float(split_line[2]),
657 infix_float(split_line[3])), cell)
658 if len(split_line) > 4:
659 force_mult = (float(split_line[4]),
660 float(split_line[5]),
661 float(split_line[6]))
662 else:
663 force_mult = None
665 positions.append((split_line[0], position, force_mult))
667 return positions
670def get_atomic_species(lines, n_species):
671 """Parse atomic species from ATOMIC_SPECIES card.
673 Parameters
674 ----------
675 lines : list[str]
676 A list of lines containing the ATOMIC_POSITIONS card.
677 n_species : int
678 Expected number of atom types. Only this many lines will be parsed.
680 Returns
681 -------
682 species : list[(str, float, str)]
684 Raises
685 ------
686 ValueError
687 Any problems parsing the data result in ValueError
688 """
690 species = None
691 # no blanks or comment lines, can the consume n_atoms lines for positions
692 trimmed_lines = (line.strip() for line in lines
693 if line.strip() and not line.startswith('#'))
695 for line in trimmed_lines:
696 if line.startswith('ATOMIC_SPECIES'):
697 if species is not None:
698 raise ValueError('Multiple ATOMIC_SPECIES specified')
700 species = []
701 for _dummy in range(n_species):
702 label_weight_pseudo = next(trimmed_lines).split()
703 species.append((label_weight_pseudo[0],
704 float(label_weight_pseudo[1]),
705 label_weight_pseudo[2]))
707 return species
710def get_cell_parameters(lines, alat=None):
711 """Parse unit cell from CELL_PARAMETERS card.
713 Parameters
714 ----------
715 lines : list[str]
716 A list with lines containing the CELL_PARAMETERS card.
717 alat : float | None
718 Unit of lattice vectors in Angstrom. Only used if the card is
719 given in units of alat. alat must be None if CELL_PARAMETERS card
720 is in Bohr or Angstrom. For output files, alat will be parsed from
721 the card header and used in preference to this value.
723 Returns
724 -------
725 cell : np.array | None
726 Cell parameters as a 3x3 array in Angstrom. If no cell is found
727 None will be returned instead.
728 cell_alat : float | None
729 If a value for alat is given in the card header, this is also
730 returned, otherwise this will be None.
732 Raises
733 ------
734 ValueError
735 If CELL_PARAMETERS are given in units of bohr or angstrom
736 and alat is not
737 """
739 cell = None
740 cell_alat = None
741 # no blanks or comment lines, can take three lines for cell
742 trimmed_lines = (line for line in lines
743 if line.strip() and not line[0] == '#')
745 for line in trimmed_lines:
746 if line.strip().startswith('CELL_PARAMETERS'):
747 if cell is not None:
748 # multiple definitions
749 raise ValueError('CELL_PARAMETERS specified multiple times')
750 # Priority and behaviour tested with QE 5.3
751 if 'bohr' in line.lower():
752 if alat is not None:
753 raise ValueError('Lattice parameters given in '
754 '&SYSTEM celldm/A and CELL_PARAMETERS '
755 'bohr')
756 cell_units = units['Bohr']
757 elif 'angstrom' in line.lower():
758 if alat is not None:
759 raise ValueError('Lattice parameters given in '
760 '&SYSTEM celldm/A and CELL_PARAMETERS '
761 'angstrom')
762 cell_units = 1.0
763 elif 'alat' in line.lower():
764 # Output file has (alat = value) (in Bohrs)
765 if '=' in line:
766 alat = float(line.strip(') \n').split()[-1]) * units['Bohr']
767 cell_alat = alat
768 elif alat is None:
769 raise ValueError('Lattice parameters must be set in '
770 '&SYSTEM for alat units')
771 cell_units = alat
772 elif alat is None:
773 # may be DEPRECATED in future
774 cell_units = units['Bohr']
775 else:
776 # may be DEPRECATED in future
777 cell_units = alat
778 # Grab the parameters; blank lines have been removed
779 cell = [[ffloat(x) for x in next(trimmed_lines).split()[:3]],
780 [ffloat(x) for x in next(trimmed_lines).split()[:3]],
781 [ffloat(x) for x in next(trimmed_lines).split()[:3]]]
782 cell = np.array(cell) * cell_units
784 return cell, cell_alat
787def str_to_value(string):
788 """Attempt to convert string into int, float (including fortran double),
789 or bool, in that order, otherwise return the string.
790 Valid (case-insensitive) bool values are: '.true.', '.t.', 'true'
791 and 't' (or false equivalents).
793 Parameters
794 ----------
795 string : str
796 Test to parse for a datatype
798 Returns
799 -------
800 value : any
801 Parsed string as the most appropriate datatype of int, float,
802 bool or string.
804 """
806 # Just an integer
807 try:
808 return int(string)
809 except ValueError:
810 pass
811 # Standard float
812 try:
813 return float(string)
814 except ValueError:
815 pass
816 # Fortran double
817 try:
818 return ffloat(string)
819 except ValueError:
820 pass
822 # possible bool, else just the raw string
823 if string.lower() in ('.true.', '.t.', 'true', 't'):
824 return True
825 elif string.lower() in ('.false.', '.f.', 'false', 'f'):
826 return False
827 else:
828 return string.strip("'")
831def read_fortran_namelist(fileobj):
832 """Takes a fortran-namelist formatted file and returns nested
833 dictionaries of sections and key-value data, followed by a list
834 of lines of text that do not fit the specifications.
836 Behaviour is taken from Quantum ESPRESSO 5.3. Parses fairly
837 convoluted files the same way that QE should, but may not get
838 all the MANDATORY rules and edge cases for very non-standard files:
839 Ignores anything after '!' in a namelist, split pairs on ','
840 to include multiple key=values on a line, read values on section
841 start and end lines, section terminating character, '/', can appear
842 anywhere on a line.
843 All of these are ignored if the value is in 'quotes'.
845 Parameters
846 ----------
847 fileobj : file
848 An open file-like object.
850 Returns
851 -------
852 data : dict of dict
853 Dictionary for each section in the namelist with key = value
854 pairs of data.
855 card_lines : list of str
856 Any lines not used to create the data, assumed to belong to 'cards'
857 in the input file.
859 """
860 # Espresso requires the correct order
861 data = Namelist()
862 card_lines = []
863 in_namelist = False
864 section = 'none' # can't be in a section without changing this
866 for line in fileobj:
867 # leading and trailing whitespace never needed
868 line = line.strip()
869 if line.startswith('&'):
870 # inside a namelist
871 section = line.split()[0][1:].lower() # case insensitive
872 if section in data:
873 # Repeated sections are completely ignored.
874 # (Note that repeated keys overwrite within a section)
875 section = "_ignored"
876 data[section] = Namelist()
877 in_namelist = True
878 if not in_namelist and line:
879 # Stripped line is Truthy, so safe to index first character
880 if line[0] not in ('!', '#'):
881 card_lines.append(line)
882 if in_namelist:
883 # parse k, v from line:
884 key = []
885 value = None
886 in_quotes = False
887 for character in line:
888 if character == ',' and value is not None and not in_quotes:
889 # finished value:
890 data[section][''.join(key).strip()] = str_to_value(
891 ''.join(value).strip())
892 key = []
893 value = None
894 elif character == '=' and value is None and not in_quotes:
895 # start writing value
896 value = []
897 elif character == "'":
898 # only found in value anyway
899 in_quotes = not in_quotes
900 value.append("'")
901 elif character == '!' and not in_quotes:
902 break
903 elif character == '/' and not in_quotes:
904 in_namelist = False
905 break
906 elif value is not None:
907 value.append(character)
908 else:
909 key.append(character)
910 if value is not None:
911 data[section][''.join(key).strip()] = str_to_value(
912 ''.join(value).strip())
914 return data, card_lines
917def ffloat(string):
918 """Parse float from fortran compatible float definitions.
920 In fortran exponents can be defined with 'd' or 'q' to symbolise
921 double or quad precision numbers. Double precision numbers are
922 converted to python floats and quad precision values are interpreted
923 as numpy longdouble values (platform specific precision).
925 Parameters
926 ----------
927 string : str
928 A string containing a number in fortran real format
930 Returns
931 -------
932 value : float | np.longdouble
933 Parsed value of the string.
935 Raises
936 ------
937 ValueError
938 Unable to parse a float value.
940 """
942 if 'q' in string.lower():
943 return np.longdouble(string.lower().replace('q', 'e'))
944 else:
945 return float(string.lower().replace('d', 'e'))
948def label_to_symbol(label):
949 """Convert a valid espresso ATOMIC_SPECIES label to a
950 chemical symbol.
952 Parameters
953 ----------
954 label : str
955 chemical symbol X (1 or 2 characters, case-insensitive)
956 or chemical symbol plus a number or a letter, as in
957 "Xn" (e.g. Fe1) or "X_*" or "X-*" (e.g. C1, C_h;
958 max total length cannot exceed 3 characters).
960 Returns
961 -------
962 symbol : str
963 The best matching species from ase.utils.chemical_symbols
965 Raises
966 ------
967 KeyError
968 Couldn't find an appropriate species.
970 Notes
971 -----
972 It's impossible to tell whether e.g. He is helium
973 or hydrogen labelled 'e'.
974 """
976 # possibly a two character species
977 # ase Atoms need proper case of chemical symbols.
978 if len(label) >= 2:
979 test_symbol = label[0].upper() + label[1].lower()
980 if test_symbol in chemical_symbols:
981 return test_symbol
982 # finally try with one character
983 test_symbol = label[0].upper()
984 if test_symbol in chemical_symbols:
985 return test_symbol
986 else:
987 raise KeyError('Could not parse species from label {}.'
988 ''.format(label))
991def infix_float(text):
992 """Parse simple infix maths into a float for compatibility with
993 Quantum ESPRESSO ATOMIC_POSITIONS cards. Note: this works with the
994 example, and most simple expressions, but the capabilities of
995 the two parsers are not identical. Will also parse a normal float
996 value properly, but slowly.
998 >>> infix_float('1/2*3^(-1/2)')
999 0.28867513459481287
1001 Parameters
1002 ----------
1003 text : str
1004 An arithmetic expression using +, -, *, / and ^, including brackets.
1006 Returns
1007 -------
1008 value : float
1009 Result of the mathematical expression.
1011 """
1013 def middle_brackets(full_text):
1014 """Extract text from innermost brackets."""
1015 start, end = 0, len(full_text)
1016 for (idx, char) in enumerate(full_text):
1017 if char == '(':
1018 start = idx
1019 if char == ')':
1020 end = idx + 1
1021 break
1022 return full_text[start:end]
1024 def eval_no_bracket_expr(full_text):
1025 """Calculate value of a mathematical expression, no brackets."""
1026 exprs = [('+', op.add), ('*', op.mul),
1027 ('/', op.truediv), ('^', op.pow)]
1028 full_text = full_text.lstrip('(').rstrip(')')
1029 try:
1030 return float(full_text)
1031 except ValueError:
1032 for symbol, func in exprs:
1033 if symbol in full_text:
1034 left, right = full_text.split(symbol, 1) # single split
1035 return func(eval_no_bracket_expr(left),
1036 eval_no_bracket_expr(right))
1038 while '(' in text:
1039 middle = middle_brackets(text)
1040 text = text.replace(middle, f'{eval_no_bracket_expr(middle)}')
1042 return float(eval_no_bracket_expr(text))
1044###
1045# Input file writing
1046###
1049# Ordered and case insensitive
1050KEYS = Namelist((
1051 ('CONTROL', [
1052 'calculation', 'title', 'verbosity', 'restart_mode', 'wf_collect',
1053 'nstep', 'iprint', 'tstress', 'tprnfor', 'dt', 'outdir', 'wfcdir',
1054 'prefix', 'lkpoint_dir', 'max_seconds', 'etot_conv_thr',
1055 'forc_conv_thr', 'disk_io', 'pseudo_dir', 'tefield', 'dipfield',
1056 'lelfield', 'nberrycyc', 'lorbm', 'lberry', 'gdir', 'nppstr',
1057 'lfcpopt', 'monopole']),
1058 ('SYSTEM', [
1059 'ibrav', 'nat', 'ntyp', 'nbnd', 'tot_charge', 'tot_magnetization',
1060 'starting_magnetization', 'ecutwfc', 'ecutrho', 'ecutfock', 'nr1',
1061 'nr2', 'nr3', 'nr1s', 'nr2s', 'nr3s', 'nosym', 'nosym_evc', 'noinv',
1062 'no_t_rev', 'force_symmorphic', 'use_all_frac', 'occupations',
1063 'one_atom_occupations', 'starting_spin_angle', 'degauss', 'smearing',
1064 'nspin', 'noncolin', 'ecfixed', 'qcutz', 'q2sigma', 'input_dft',
1065 'exx_fraction', 'screening_parameter', 'exxdiv_treatment',
1066 'x_gamma_extrapolation', 'ecutvcut', 'nqx1', 'nqx2', 'nqx3',
1067 'lda_plus_u', 'lda_plus_u_kind', 'Hubbard_U', 'Hubbard_J0',
1068 'Hubbard_alpha', 'Hubbard_beta', 'Hubbard_J',
1069 'starting_ns_eigenvalue', 'U_projection_type', 'edir',
1070 'emaxpos', 'eopreg', 'eamp', 'angle1', 'angle2',
1071 'constrained_magnetization', 'fixed_magnetization', 'lambda',
1072 'report', 'lspinorb', 'assume_isolated', 'esm_bc', 'esm_w',
1073 'esm_efield', 'esm_nfit', 'fcp_mu', 'vdw_corr', 'london',
1074 'london_s6', 'london_c6', 'london_rvdw', 'london_rcut',
1075 'ts_vdw_econv_thr', 'ts_vdw_isolated', 'xdm', 'xdm_a1', 'xdm_a2',
1076 'space_group', 'uniqueb', 'origin_choice', 'rhombohedral', 'zmon',
1077 'realxz', 'block', 'block_1', 'block_2', 'block_height']),
1078 ('ELECTRONS', [
1079 'electron_maxstep', 'scf_must_converge', 'conv_thr', 'adaptive_thr',
1080 'conv_thr_init', 'conv_thr_multi', 'mixing_mode', 'mixing_beta',
1081 'mixing_ndim', 'mixing_fixed_ns', 'diagonalization', 'ortho_para',
1082 'diago_thr_init', 'diago_cg_maxiter', 'diago_david_ndim',
1083 'diago_full_acc', 'efield', 'efield_cart', 'efield_phase',
1084 'startingpot', 'startingwfc', 'tqr']),
1085 ('IONS', [
1086 'ion_dynamics', 'ion_positions', 'pot_extrapolation',
1087 'wfc_extrapolation', 'remove_rigid_rot', 'ion_temperature', 'tempw',
1088 'tolp', 'delta_t', 'nraise', 'refold_pos', 'upscale', 'bfgs_ndim',
1089 'trust_radius_max', 'trust_radius_min', 'trust_radius_ini', 'w_1',
1090 'w_2']),
1091 ('CELL', [
1092 'cell_dynamics', 'press', 'wmass', 'cell_factor', 'press_conv_thr',
1093 'cell_dofree'])))
1096# Number of valence electrons in the pseudopotentials recommended by
1097# http://materialscloud.org/sssp/. These are just used as a fallback for
1098# calculating initial magetization values which are given as a fraction
1099# of valence electrons.
1100SSSP_VALENCE = [
1101 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,
1102 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,
1103 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,
1104 13.0, 14.0, 15.0, 16.0, 17.0, 18.0, 19.0, 12.0, 13.0, 14.0, 15.0, 6.0,
1105 7.0, 18.0, 9.0, 10.0, 11.0, 12.0, 13.0, 14.0, 15.0, 16.0, 17.0, 18.0,
1106 19.0, 20.0, 21.0, 22.0, 23.0, 24.0, 25.0, 36.0, 27.0, 14.0, 15.0, 30.0,
1107 15.0, 32.0, 19.0, 12.0, 13.0, 14.0, 15.0, 16.0, 18.0]
1110def construct_namelist(parameters=None, keys=None, warn=False, **kwargs):
1111 """
1112 Construct an ordered Namelist containing all the parameters given (as
1113 a dictionary or kwargs). Keys will be inserted into their appropriate
1114 section in the namelist and the dictionary may contain flat and nested
1115 structures. Any kwargs that match input keys will be incorporated into
1116 their correct section. All matches are case-insensitive, and returned
1117 Namelist object is a case-insensitive dict.
1119 If a key is not known to ase, but in a section within `parameters`,
1120 it will be assumed that it was put there on purpose and included
1121 in the output namelist. Anything not in a section will be ignored (set
1122 `warn` to True to see ignored keys).
1124 Keys with a dimension (e.g. Hubbard_U(1)) will be incorporated as-is
1125 so the `i` should be made to match the output.
1127 The priority of the keys is:
1128 kwargs[key] > parameters[key] > parameters[section][key]
1129 Only the highest priority item will be included.
1131 Parameters
1132 ----------
1133 parameters: dict
1134 Flat or nested set of input parameters.
1135 keys: Namelist | dict
1136 Namelist to use as a template for the output.
1137 warn: bool
1138 Enable warnings for unused keys.
1140 Returns
1141 -------
1142 input_namelist: Namelist
1143 pw.x compatible namelist of input parameters.
1145 """
1147 if keys is None:
1148 keys = deepcopy(KEYS)
1149 # Convert everything to Namelist early to make case-insensitive
1150 if parameters is None:
1151 parameters = Namelist()
1152 else:
1153 # Maximum one level of nested dict
1154 # Don't modify in place
1155 parameters_namelist = Namelist()
1156 for key, value in parameters.items():
1157 if isinstance(value, dict):
1158 parameters_namelist[key] = Namelist(value)
1159 else:
1160 parameters_namelist[key] = value
1161 parameters = parameters_namelist
1163 # Just a dict
1164 kwargs = Namelist(kwargs)
1166 # Final parameter set
1167 input_namelist = Namelist()
1169 # Collect
1170 for section in keys:
1171 sec_list = Namelist()
1172 for key in keys[section]:
1173 # Check all three separately and pop them all so that
1174 # we can check for missing values later
1175 value = None
1177 if key in parameters.get(section, {}):
1178 value = parameters[section].pop(key)
1179 if key in parameters:
1180 value = parameters.pop(key)
1181 if key in kwargs:
1182 value = kwargs.pop(key)
1184 if value is not None:
1185 sec_list[key] = value
1187 # Check if there is a key(i) version (no extra parsing)
1188 for arg_key in list(parameters.get(section, {})):
1189 if arg_key.split('(')[0].strip().lower() == key.lower():
1190 sec_list[arg_key] = parameters[section].pop(arg_key)
1191 cp_parameters = parameters.copy()
1192 for arg_key in cp_parameters:
1193 if arg_key.split('(')[0].strip().lower() == key.lower():
1194 sec_list[arg_key] = parameters.pop(arg_key)
1195 cp_kwargs = kwargs.copy()
1196 for arg_key in cp_kwargs:
1197 if arg_key.split('(')[0].strip().lower() == key.lower():
1198 sec_list[arg_key] = kwargs.pop(arg_key)
1200 # Add to output
1201 input_namelist[section] = sec_list
1203 unused_keys = list(kwargs)
1204 # pass anything else already in a section
1205 for key, value in parameters.items():
1206 if key in keys and isinstance(value, dict):
1207 input_namelist[key].update(value)
1208 elif isinstance(value, dict):
1209 unused_keys.extend(list(value))
1210 else:
1211 unused_keys.append(key)
1213 if warn and unused_keys:
1214 warnings.warn('Unused keys: {}'.format(', '.join(unused_keys)))
1216 return input_namelist
1219def kspacing_to_grid(atoms, spacing, calculated_spacing=None):
1220 """
1221 Calculate the kpoint mesh that is equivalent to the given spacing
1222 in reciprocal space (units Angstrom^-1). The number of kpoints is each
1223 dimension is rounded up (compatible with CASTEP).
1225 Parameters
1226 ----------
1227 atoms: ase.Atoms
1228 A structure that can have get_reciprocal_cell called on it.
1229 spacing: float
1230 Minimum K-Point spacing in $A^{-1}$.
1231 calculated_spacing : list
1232 If a three item list (or similar mutable sequence) is given the
1233 members will be replaced with the actual calculated spacing in
1234 $A^{-1}$.
1236 Returns
1237 -------
1238 kpoint_grid : [int, int, int]
1239 MP grid specification to give the required spacing.
1241 """
1242 # No factor of 2pi in ase, everything in A^-1
1243 # reciprocal dimensions
1244 r_x, r_y, r_z = np.linalg.norm(atoms.cell.reciprocal(), axis=1)
1246 kpoint_grid = [int(r_x / spacing) + 1,
1247 int(r_y / spacing) + 1,
1248 int(r_z / spacing) + 1]
1250 for i, _ in enumerate(kpoint_grid):
1251 if not atoms.pbc[i]:
1252 kpoint_grid[i] = 1
1254 if calculated_spacing is not None:
1255 calculated_spacing[:] = [r_x / kpoint_grid[0],
1256 r_y / kpoint_grid[1],
1257 r_z / kpoint_grid[2]]
1259 return kpoint_grid
1262def format_atom_position(atom, crystal_coordinates, mask='', tidx=None):
1263 """Format one line of atomic positions in
1264 Quantum ESPRESSO ATOMIC_POSITIONS card.
1266 >>> for atom in make_supercell(bulk('Li', 'bcc'), np.ones(3)-np.eye(3)):
1267 >>> format_atom_position(atom, True)
1268 Li 0.0000000000 0.0000000000 0.0000000000
1269 Li 0.5000000000 0.5000000000 0.5000000000
1271 Parameters
1272 ----------
1273 atom : Atom
1274 A structure that has symbol and [position | (a, b, c)].
1275 crystal_coordinates: bool
1276 Whether the atomic positions should be written to the QE input file in
1277 absolute (False, default) or relative (crystal) coordinates (True).
1278 mask, optional : str
1279 String of ndim=3 0 or 1 for constraining atomic positions.
1280 tidx, optional : int
1281 Magnetic type index.
1283 Returns
1284 -------
1285 atom_line : str
1286 Input line for atom position
1287 """
1288 if crystal_coordinates:
1289 coords = [atom.a, atom.b, atom.c]
1290 else:
1291 coords = atom.position
1292 line_fmt = '{atom.symbol}'
1293 inps = dict(atom=atom)
1294 if tidx is not None:
1295 line_fmt += '{tidx}'
1296 inps["tidx"] = tidx
1297 line_fmt += ' {coords[0]:.10f} {coords[1]:.10f} {coords[2]:.10f} '
1298 inps["coords"] = coords
1299 line_fmt += ' ' + mask + '\n'
1300 astr = line_fmt.format(**inps)
1301 return astr
1304def namelist_to_string(input_parameters):
1305 """Format a Namelist object as a string for writing to a file.
1306 Assume sections are ordered (taken care of in namelist construction)
1307 and that repr converts to a QE readable representation (except bools)
1309 Parameters
1310 ----------
1311 input_parameters : Namelist | dict
1312 Expecting a nested dictionary of sections and key-value data.
1314 Returns
1315 -------
1316 pwi : List[str]
1317 Input line for the namelist
1318 """
1319 pwi = []
1320 for section in input_parameters:
1321 pwi.append(f'&{section.upper()}\n')
1322 for key, value in input_parameters[section].items():
1323 if value is True:
1324 pwi.append(f' {key:16} = .true.\n')
1325 elif value is False:
1326 pwi.append(f' {key:16} = .false.\n')
1327 else:
1328 # repr format to get quotes around strings
1329 pwi.append(f' {key:16} = {value!r}\n')
1330 pwi.append('/\n') # terminate section
1331 pwi.append('\n')
1332 return pwi
1335def write_espresso_in(fd, atoms, input_data=None, pseudopotentials=None,
1336 kspacing=None, kpts=None, koffset=(0, 0, 0),
1337 crystal_coordinates=False, **kwargs):
1338 """
1339 Create an input file for pw.x.
1341 Use set_initial_magnetic_moments to turn on spin, if ispin is set to 2
1342 with no magnetic moments, they will all be set to 0.0. Magnetic moments
1343 will be converted to the QE units (fraction of valence electrons) using
1344 any pseudopotential files found, or a best guess for the number of
1345 valence electrons.
1347 Units are not converted for any other input data, so use Quantum ESPRESSO
1348 units (Usually Ry or atomic units).
1350 Keys with a dimension (e.g. Hubbard_U(1)) will be incorporated as-is
1351 so the `i` should be made to match the output.
1353 Implemented features:
1355 - Conversion of :class:`ase.constraints.FixAtoms` and
1356 :class:`ase.constraints.FixCartesian`.
1357 - `starting_magnetization` derived from the `mgmoms` and pseudopotentials
1358 (searches default paths for pseudo files.)
1359 - Automatic assignment of options to their correct sections.
1361 Not implemented:
1363 - Non-zero values of ibrav
1364 - Lists of k-points
1365 - Other constraints
1366 - Hubbard parameters
1367 - Validation of the argument types for input
1368 - Validation of required options
1370 Parameters
1371 ----------
1372 fd: file
1373 A file like object to write the input file to.
1374 atoms: Atoms
1375 A single atomistic configuration to write to `fd`.
1376 input_data: dict
1377 A flat or nested dictionary with input parameters for pw.x
1378 pseudopotentials: dict
1379 A filename for each atomic species, e.g.
1380 {'O': 'O.pbe-rrkjus.UPF', 'H': 'H.pbe-rrkjus.UPF'}.
1381 A dummy name will be used if none are given.
1382 kspacing: float
1383 Generate a grid of k-points with this as the minimum distance,
1384 in A^-1 between them in reciprocal space. If set to None, kpts
1385 will be used instead.
1386 kpts: (int, int, int) or dict
1387 If kpts is a tuple (or list) of 3 integers, it is interpreted
1388 as the dimensions of a Monkhorst-Pack grid.
1389 If ``kpts`` is set to ``None``, only the Γ-point will be included
1390 and QE will use routines optimized for Γ-point-only calculations.
1391 Compared to Γ-point-only calculations without this optimization
1392 (i.e. with ``kpts=(1, 1, 1)``), the memory and CPU requirements
1393 are typically reduced by half.
1394 If kpts is a dict, it will either be interpreted as a path
1395 in the Brillouin zone (*) if it contains the 'path' keyword,
1396 otherwise it is converted to a Monkhorst-Pack grid (**).
1397 (*) see ase.dft.kpoints.bandpath
1398 (**) see ase.calculators.calculator.kpts2sizeandoffsets
1399 koffset: (int, int, int)
1400 Offset of kpoints in each direction. Must be 0 (no offset) or
1401 1 (half grid offset). Setting to True is equivalent to (1, 1, 1).
1402 crystal_coordinates: bool
1403 Whether the atomic positions should be written to the QE input file in
1404 absolute (False, default) or relative (crystal) coordinates (True).
1406 """
1408 # Convert to a namelist to make working with parameters much easier
1409 # Note that the name ``input_data`` is chosen to prevent clash with
1410 # ``parameters`` in Calculator objects
1411 input_parameters = construct_namelist(input_data, **kwargs)
1413 # Convert ase constraints to QE constraints
1414 # Nx3 array of force multipliers matches what QE uses
1415 # Do this early so it is available when constructing the atoms card
1416 constraint_mask = np.ones((len(atoms), 3), dtype='int')
1417 for constraint in atoms.constraints:
1418 if isinstance(constraint, FixAtoms):
1419 constraint_mask[constraint.index] = 0
1420 elif isinstance(constraint, FixCartesian):
1421 constraint_mask[constraint.a] = constraint.mask
1422 else:
1423 warnings.warn(f'Ignored unknown constraint {constraint}')
1424 masks = []
1425 for atom in atoms:
1426 # only inclued mask if something is fixed
1427 if not all(constraint_mask[atom.index]):
1428 mask = ' {mask[0]} {mask[1]} {mask[2]}'.format(
1429 mask=constraint_mask[atom.index])
1430 else:
1431 mask = ''
1432 masks.append(mask)
1434 # Species info holds the information on the pseudopotential and
1435 # associated for each element
1436 if pseudopotentials is None:
1437 pseudopotentials = {}
1438 species_info = {}
1439 for species in set(atoms.get_chemical_symbols()):
1440 # Look in all possible locations for the pseudos and try to figure
1441 # out the number of valence electrons
1442 pseudo = pseudopotentials[species]
1443 species_info[species] = {'pseudo': pseudo}
1445 # Convert atoms into species.
1446 # Each different magnetic moment needs to be a separate type even with
1447 # the same pseudopotential (e.g. an up and a down for AFM).
1448 # if any magmom are > 0 or nspin == 2 then use species labels.
1449 # Rememeber: magnetisation uses 1 based indexes
1450 atomic_species = OrderedDict()
1451 atomic_species_str = []
1452 atomic_positions_str = []
1454 nspin = input_parameters['system'].get('nspin', 1) # 1 is the default
1455 noncolin = input_parameters['system'].get('noncolin', False)
1456 rescale_magmom_fac = kwargs.get('rescale_magmom_fac', 1.0)
1457 if any(atoms.get_initial_magnetic_moments()):
1458 if nspin == 1 and not noncolin:
1459 # Force spin on
1460 input_parameters['system']['nspin'] = 2
1461 nspin = 2
1463 if nspin == 2 or noncolin:
1464 # Magnetic calculation on
1465 for atom, mask, magmom in zip(
1466 atoms, masks, atoms.get_initial_magnetic_moments()):
1467 if (atom.symbol, magmom) not in atomic_species:
1468 # for qe version 7.2 or older magmon must be rescale by
1469 # about a factor 10 to assume sensible values
1470 # since qe-v7.3 magmom values will be provided unscaled
1471 fspin = float(magmom) / rescale_magmom_fac
1472 # Index in the atomic species list
1473 sidx = len(atomic_species) + 1
1474 # Index for that atom type; no index for first one
1475 tidx = sum(atom.symbol == x[0] for x in atomic_species) or ' '
1476 atomic_species[(atom.symbol, magmom)] = (sidx, tidx)
1477 # Add magnetization to the input file
1478 mag_str = f"starting_magnetization({sidx})"
1479 input_parameters['system'][mag_str] = fspin
1480 species_pseudo = species_info[atom.symbol]['pseudo']
1481 atomic_species_str.append(
1482 f"{atom.symbol}{tidx} {atom.mass} {species_pseudo}\n")
1483 # lookup tidx to append to name
1484 sidx, tidx = atomic_species[(atom.symbol, magmom)]
1485 # construct line for atomic positions
1486 atomic_positions_str.append(
1487 format_atom_position(
1488 atom, crystal_coordinates, mask=mask, tidx=tidx)
1489 )
1490 else:
1491 # Do nothing about magnetisation
1492 for atom, mask in zip(atoms, masks):
1493 if atom.symbol not in atomic_species:
1494 atomic_species[atom.symbol] = True # just a placeholder
1495 species_pseudo = species_info[atom.symbol]['pseudo']
1496 atomic_species_str.append(
1497 f"{atom.symbol} {atom.mass} {species_pseudo}\n")
1498 # construct line for atomic positions
1499 atomic_positions_str.append(
1500 format_atom_position(atom, crystal_coordinates, mask=mask)
1501 )
1503 # Add computed parameters
1504 # different magnetisms means different types
1505 input_parameters['system']['ntyp'] = len(atomic_species)
1506 input_parameters['system']['nat'] = len(atoms)
1508 # Use cell as given or fit to a specific ibrav
1509 if 'ibrav' in input_parameters['system']:
1510 ibrav = input_parameters['system']['ibrav']
1511 if ibrav != 0:
1512 raise ValueError(ibrav_error_message)
1513 else:
1514 # Just use standard cell block
1515 input_parameters['system']['ibrav'] = 0
1517 # Construct input file into this
1518 pwi = namelist_to_string(input_parameters)
1520 # Pseudopotentials
1521 pwi.append('ATOMIC_SPECIES\n')
1522 pwi.extend(atomic_species_str)
1523 pwi.append('\n')
1525 # KPOINTS - add a MP grid as required
1526 if kspacing is not None:
1527 kgrid = kspacing_to_grid(atoms, kspacing)
1528 elif kpts is not None:
1529 if isinstance(kpts, dict) and 'path' not in kpts:
1530 kgrid, shift = kpts2sizeandoffsets(atoms=atoms, **kpts)
1531 koffset = []
1532 for i, x in enumerate(shift):
1533 assert x == 0 or abs(x * kgrid[i] - 0.5) < 1e-14
1534 koffset.append(0 if x == 0 else 1)
1535 else:
1536 kgrid = kpts
1537 else:
1538 kgrid = "gamma"
1540 # True and False work here and will get converted by ':d' format
1541 if isinstance(koffset, int):
1542 koffset = (koffset, ) * 3
1544 # BandPath object or bandpath-as-dictionary:
1545 if isinstance(kgrid, dict) or hasattr(kgrid, 'kpts'):
1546 pwi.append('K_POINTS crystal_b\n')
1547 assert hasattr(kgrid, 'path') or 'path' in kgrid
1548 kgrid = kpts2ndarray(kgrid, atoms=atoms)
1549 pwi.append(f'{len(kgrid)}\n')
1550 for k in kgrid:
1551 pwi.append(f"{k[0]:.14f} {k[1]:.14f} {k[2]:.14f} 0\n")
1552 pwi.append('\n')
1553 elif isinstance(kgrid, str) and (kgrid == "gamma"):
1554 pwi.append('K_POINTS gamma\n')
1555 pwi.append('\n')
1556 else:
1557 pwi.append('K_POINTS automatic\n')
1558 pwi.append(f"{kgrid[0]} {kgrid[1]} {kgrid[2]} "
1559 f" {koffset[0]:d} {koffset[1]:d} {koffset[2]:d}\n")
1560 pwi.append('\n')
1562 # CELL block, if required
1563 if input_parameters['SYSTEM']['ibrav'] == 0:
1564 pwi.append('CELL_PARAMETERS angstrom\n')
1565 pwi.append('{cell[0][0]:.14f} {cell[0][1]:.14f} {cell[0][2]:.14f}\n'
1566 '{cell[1][0]:.14f} {cell[1][1]:.14f} {cell[1][2]:.14f}\n'
1567 '{cell[2][0]:.14f} {cell[2][1]:.14f} {cell[2][2]:.14f}\n'
1568 ''.format(cell=atoms.cell))
1569 pwi.append('\n')
1571 # Positions - already constructed, but must appear after namelist
1572 if crystal_coordinates:
1573 pwi.append('ATOMIC_POSITIONS crystal\n')
1574 else:
1575 pwi.append('ATOMIC_POSITIONS angstrom\n')
1576 pwi.extend(atomic_positions_str)
1577 pwi.append('\n')
1579 # DONE!
1580 fd.write(''.join(pwi))