Coverage for /builds/ase/ase/ase/io/espresso.py : 74.29%
Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
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 os
16import operator as op
17import re
18import warnings
19from collections import OrderedDict
20from os import path
22import numpy as np
24from ase.atoms import Atoms
25from ase.cell import Cell
26from ase.calculators.singlepoint import (SinglePointDFTCalculator,
27 SinglePointKPoint)
28from ase.calculators.calculator import kpts2ndarray, kpts2sizeandoffsets
29from ase.dft.kpoints import kpoint_convert
30from ase.constraints import FixAtoms, FixCartesian
31from ase.data import chemical_symbols, atomic_numbers
32from ase.units import create_units
33from ase.utils import iofunction
36# Quantum ESPRESSO uses CODATA 2006 internally
37units = create_units('2006')
39# Section identifiers
40_PW_START = 'Program PWSCF'
41_PW_END = 'End of self-consistent calculation'
42_PW_CELL = 'CELL_PARAMETERS'
43_PW_POS = 'ATOMIC_POSITIONS'
44_PW_MAGMOM = 'Magnetic moment per site'
45_PW_FORCE = 'Forces acting on atoms'
46_PW_TOTEN = '! total energy'
47_PW_STRESS = 'total stress'
48_PW_FERMI = 'the Fermi energy is'
49_PW_HIGHEST_OCCUPIED = 'highest occupied level'
50_PW_HIGHEST_OCCUPIED_LOWEST_FREE = 'highest occupied, lowest unoccupied level'
51_PW_KPTS = 'number of k points='
52_PW_BANDS = _PW_END
53_PW_BANDSTRUCTURE = 'End of band structure calculation'
55# ibrav error message
56ibrav_error_message = (
57 'ASE does not support ibrav != 0. Note that with ibrav '
58 '== 0, Quantum ESPRESSO will still detect the symmetries '
59 'of your system because the CELL_PARAMETERS are defined '
60 'to a high level of precision.')
63class Namelist(OrderedDict):
64 """Case insensitive dict that emulates Fortran Namelists."""
66 def __contains__(self, key):
67 return super(Namelist, self).__contains__(key.lower())
69 def __delitem__(self, key):
70 return super(Namelist, self).__delitem__(key.lower())
72 def __getitem__(self, key):
73 return super(Namelist, self).__getitem__(key.lower())
75 def __setitem__(self, key, value):
76 super(Namelist, self).__setitem__(key.lower(), value)
78 def get(self, key, default=None):
79 return super(Namelist, self).get(key.lower(), default)
82@iofunction('rU')
83def read_espresso_out(fileobj, index=-1, results_required=True):
84 """Reads Quantum ESPRESSO output files.
86 The atomistic configurations as well as results (energy, force, stress,
87 magnetic moments) of the calculation are read for all configurations
88 within the output file.
90 Will probably raise errors for broken or incomplete files.
92 Parameters
93 ----------
94 fileobj : file|str
95 A file like object or filename
96 index : slice
97 The index of configurations to extract.
98 results_required : bool
99 If True, atomistic configurations that do not have any
100 associated results will not be included. This prevents double
101 printed configurations and incomplete calculations from being
102 returned as the final configuration with no results data.
104 Yields
105 ------
106 structure : Atoms
107 The next structure from the index slice. The Atoms has a
108 SinglePointCalculator attached with any results parsed from
109 the file.
112 """
113 # work with a copy in memory for faster random access
114 pwo_lines = fileobj.readlines()
116 # TODO: index -1 special case?
117 # Index all the interesting points
118 indexes = {
119 _PW_START: [],
120 _PW_END: [],
121 _PW_CELL: [],
122 _PW_POS: [],
123 _PW_MAGMOM: [],
124 _PW_FORCE: [],
125 _PW_TOTEN: [],
126 _PW_STRESS: [],
127 _PW_FERMI: [],
128 _PW_HIGHEST_OCCUPIED: [],
129 _PW_HIGHEST_OCCUPIED_LOWEST_FREE: [],
130 _PW_KPTS: [],
131 _PW_BANDS: [],
132 _PW_BANDSTRUCTURE: [],
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 = dict((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 # Fermi level / highest occupied level
279 efermi = None
280 for fermi_index in indexes[_PW_FERMI]:
281 if image_index < fermi_index < next_index:
282 efermi = float(pwo_lines[fermi_index].split()[-2])
284 if efermi is None:
285 for ho_index in indexes[_PW_HIGHEST_OCCUPIED]:
286 if image_index < ho_index < next_index:
287 efermi = float(pwo_lines[ho_index].split()[-1])
289 if efermi is None:
290 for holf_index in indexes[_PW_HIGHEST_OCCUPIED_LOWEST_FREE]:
291 if image_index < holf_index < next_index:
292 efermi = float(pwo_lines[holf_index].split()[-2])
294 # K-points
295 ibzkpts = None
296 weights = None
297 kpoints_warning = "Number of k-points >= 100: " + \
298 "set verbosity='high' to print them."
300 for kpts_index in indexes[_PW_KPTS]:
301 nkpts = int(pwo_lines[kpts_index].split()[4])
302 kpts_index += 2
304 if pwo_lines[kpts_index].strip() == kpoints_warning:
305 continue
307 # QE prints the k-points in units of 2*pi/alat
308 # with alat defined as the length of the first
309 # cell vector
310 cell = structure.get_cell()
311 alat = np.linalg.norm(cell[0])
312 ibzkpts = []
313 weights = []
314 for i in range(nkpts):
315 L = pwo_lines[kpts_index + i].split()
316 weights.append(float(L[-1]))
317 coord = np.array([L[-6], L[-5], L[-4].strip('),')],
318 dtype=float)
319 coord *= 2 * np.pi / alat
320 coord = kpoint_convert(cell, ckpts_kv=coord)
321 ibzkpts.append(coord)
322 ibzkpts = np.array(ibzkpts)
323 weights = np.array(weights)
325 # Bands
326 kpts = None
327 kpoints_warning = "Number of k-points >= 100: " + \
328 "set verbosity='high' to print the bands."
330 for bands_index in indexes[_PW_BANDS] + indexes[_PW_BANDSTRUCTURE]:
331 if image_index < bands_index < next_index:
332 bands_index += 1
333 # skip over the lines with DFT+U occupation matrices
334 if 'enter write_ns' in pwo_lines[bands_index]:
335 while 'exit write_ns' not in pwo_lines[bands_index]:
336 bands_index += 1
337 bands_index += 1
339 if pwo_lines[bands_index].strip() == kpoints_warning:
340 continue
342 assert ibzkpts is not None
343 spin, bands, eigenvalues = 0, [], [[], []]
345 while True:
346 L = pwo_lines[bands_index].replace('-', ' -').split()
347 if len(L) == 0:
348 if len(bands) > 0:
349 eigenvalues[spin].append(bands)
350 bands = []
351 elif L == ['occupation', 'numbers']:
352 # Skip the lines with the occupation numbers
353 bands_index += len(eigenvalues[spin][0]) // 8 + 1
354 elif L[0] == 'k' and L[1].startswith('='):
355 pass
356 elif 'SPIN' in L:
357 if 'DOWN' in L:
358 spin += 1
359 else:
360 try:
361 bands.extend(map(float, L))
362 except ValueError:
363 break
364 bands_index += 1
366 if spin == 1:
367 assert len(eigenvalues[0]) == len(eigenvalues[1])
368 assert len(eigenvalues[0]) == len(ibzkpts), \
369 (np.shape(eigenvalues), len(ibzkpts))
371 kpts = []
372 for s in range(spin + 1):
373 for w, k, e in zip(weights, ibzkpts, eigenvalues[s]):
374 kpt = SinglePointKPoint(w, s, k, eps_n=e)
375 kpts.append(kpt)
377 # Put everything together
378 #
379 # I have added free_energy. Can and should we distinguish
380 # energy and free_energy? --askhl
381 calc = SinglePointDFTCalculator(structure, energy=energy,
382 free_energy=energy,
383 forces=forces, stress=stress,
384 magmoms=magmoms, efermi=efermi,
385 ibzkpts=ibzkpts)
386 calc.kpts = kpts
387 structure.calc = calc
389 yield structure
392def parse_pwo_start(lines, index=0):
393 """Parse Quantum ESPRESSO calculation info from lines,
394 starting from index. Return a dictionary containing extracted
395 information.
397 - `celldm(1)`: lattice parameters (alat)
398 - `cell`: unit cell in Angstrom
399 - `symbols`: element symbols for the structure
400 - `positions`: cartesian coordinates of atoms in Angstrom
401 - `atoms`: an `ase.Atoms` object constructed from the extracted data
403 Parameters
404 ----------
405 lines : list[str]
406 Contents of PWSCF output file.
407 index : int
408 Line number to begin parsing. Only first calculation will
409 be read.
411 Returns
412 -------
413 info : dict
414 Dictionary of calculation parameters, including `celldm(1)`, `cell`,
415 `symbols`, `positions`, `atoms`.
417 Raises
418 ------
419 KeyError
420 If interdependent values cannot be found (especially celldm(1))
421 an error will be raised as other quantities cannot then be
422 calculated (e.g. cell and positions).
423 """
424 # TODO: extend with extra DFT info?
426 info = {}
428 for idx, line in enumerate(lines[index:], start=index):
429 if 'celldm(1)' in line:
430 # celldm(1) has more digits than alat!!
431 info['celldm(1)'] = float(line.split()[1]) * units['Bohr']
432 info['alat'] = info['celldm(1)']
433 elif 'number of atoms/cell' in line:
434 info['nat'] = int(line.split()[-1])
435 elif 'number of atomic types' in line:
436 info['ntyp'] = int(line.split()[-1])
437 elif 'crystal axes:' in line:
438 info['cell'] = info['celldm(1)'] * np.array([
439 [float(x) for x in lines[idx + 1].split()[3:6]],
440 [float(x) for x in lines[idx + 2].split()[3:6]],
441 [float(x) for x in lines[idx + 3].split()[3:6]]])
442 elif 'positions (alat units)' in line:
443 info['symbols'], info['positions'] = [], []
445 for at_line in lines[idx + 1:idx + 1 + info['nat']]:
446 sym, x, y, z = parse_position_line(at_line)
447 info['symbols'].append(label_to_symbol(sym))
448 info['positions'].append([x * info['celldm(1)'],
449 y * info['celldm(1)'],
450 z * info['celldm(1)']])
451 # This should be the end of interesting info.
452 # Break here to avoid dealing with large lists of kpoints.
453 # Will need to be extended for DFTCalculator info.
454 break
456 # Make atoms for convenience
457 info['atoms'] = Atoms(symbols=info['symbols'],
458 positions=info['positions'],
459 cell=info['cell'], pbc=True)
461 return info
464def parse_position_line(line):
465 """Parse a single line from a pw.x output file.
467 The line must contain information about the atomic symbol and the position,
468 e.g.
470 995 Sb tau( 995) = ( 1.4212023 0.7037863 0.1242640 )
472 Parameters
473 ----------
474 line : str
475 Line to be parsed.
477 Returns
478 -------
479 sym : str
480 Atomic symbol.
481 x : float
482 x-position.
483 y : float
484 y-position.
485 z : float
486 z-position.
487 """
488 pat = re.compile(r'\s*\d+\s*(\S+)\s*tau\(\s*\d+\)\s*='
489 r'\s*\(\s*(\S+)\s+(\S+)\s+(\S+)\s*\)')
490 match = pat.match(line)
491 assert match is not None
492 sym, x, y, z = match.group(1, 2, 3, 4)
493 return sym, float(x), float(y), float(z)
496@iofunction('rU')
497def read_espresso_in(fileobj):
498 """Parse a Quantum ESPRESSO input files, '.in', '.pwi'.
500 ESPRESSO inputs are generally a fortran-namelist format with custom
501 blocks of data. The namelist is parsed as a dict and an atoms object
502 is constructed from the included information.
504 Parameters
505 ----------
506 fileobj : file | str
507 A file-like object that supports line iteration with the contents
508 of the input file, or a filename.
510 Returns
511 -------
512 atoms : Atoms
513 Structure defined in the input file.
515 Raises
516 ------
517 KeyError
518 Raised for missing keys that are required to process the file
519 """
520 # parse namelist section and extract remaining lines
521 data, card_lines = read_fortran_namelist(fileobj)
523 # get the cell if ibrav=0
524 if 'system' not in data:
525 raise KeyError('Required section &SYSTEM not found.')
526 elif 'ibrav' not in data['system']:
527 raise KeyError('ibrav is required in &SYSTEM')
528 elif data['system']['ibrav'] == 0:
529 # celldm(1) is in Bohr, A is in angstrom. celldm(1) will be
530 # used even if A is also specified.
531 if 'celldm(1)' in data['system']:
532 alat = data['system']['celldm(1)'] * units['Bohr']
533 elif 'A' in data['system']:
534 alat = data['system']['A']
535 else:
536 alat = None
537 cell, _ = get_cell_parameters(card_lines, alat=alat)
538 else:
539 raise ValueError(ibrav_error_message)
541 # species_info holds some info for each element
542 species_card = get_atomic_species(
543 card_lines, n_species=data['system']['ntyp'])
544 species_info = {}
545 for ispec, (label, weight, pseudo) in enumerate(species_card):
546 symbol = label_to_symbol(label)
547 valence = get_valence_electrons(symbol, data, pseudo)
549 # starting_magnetization is in fractions of valence electrons
550 magnet_key = "starting_magnetization({0})".format(ispec + 1)
551 magmom = valence * data["system"].get(magnet_key, 0.0)
552 species_info[symbol] = {"weight": weight, "pseudo": pseudo,
553 "valence": valence, "magmom": magmom}
555 positions_card = get_atomic_positions(
556 card_lines, n_atoms=data['system']['nat'], cell=cell, alat=alat)
558 symbols = [label_to_symbol(position[0]) for position in positions_card]
559 positions = [position[1] for position in positions_card]
560 magmoms = [species_info[symbol]["magmom"] for symbol in symbols]
562 # TODO: put more info into the atoms object
563 # e.g magmom, forces.
564 atoms = Atoms(symbols=symbols, positions=positions, cell=cell, pbc=True,
565 magmoms=magmoms)
567 return atoms
570def ibrav_to_cell(system):
571 """
572 Convert a value of ibrav to a cell. Any unspecified lattice dimension
573 is set to 0.0, but will not necessarily raise an error. Also return the
574 lattice parameter.
576 Parameters
577 ----------
578 system : dict
579 The &SYSTEM section of the input file, containing the 'ibrav' setting,
580 and either celldm(1)..(6) or a, b, c, cosAB, cosAC, cosBC.
582 Returns
583 -------
584 cell : Cell
585 The cell as an ASE Cell object
587 Raises
588 ------
589 KeyError
590 Raise an error if any required keys are missing.
591 NotImplementedError
592 Only a limited number of ibrav settings can be parsed. An error
593 is raised if the ibrav interpretation is not implemented.
594 """
595 if 'celldm(1)' in system and 'a' in system:
596 raise KeyError('do not specify both celldm and a,b,c!')
597 elif 'celldm(1)' in system:
598 # celldm(x) in bohr
599 alat = system['celldm(1)'] * units['Bohr']
600 b_over_a = system.get('celldm(2)', 0.0)
601 c_over_a = system.get('celldm(3)', 0.0)
602 cosab = system.get('celldm(4)', 0.0)
603 cosac = system.get('celldm(5)', 0.0)
604 cosbc = 0.0
605 if system['ibrav'] == 14:
606 cosbc = system.get('celldm(4)', 0.0)
607 cosac = system.get('celldm(5)', 0.0)
608 cosab = system.get('celldm(6)', 0.0)
609 elif 'a' in system:
610 # a, b, c, cosAB, cosAC, cosBC in Angstrom
611 raise NotImplementedError(
612 'params_to_cell() does not yet support A/B/C/cosAB/cosAC/cosBC')
613 else:
614 raise KeyError("Missing celldm(1)")
616 if system['ibrav'] == 1:
617 cell = np.identity(3) * alat
618 elif system['ibrav'] == 2:
619 cell = np.array([[-1.0, 0.0, 1.0],
620 [0.0, 1.0, 1.0],
621 [-1.0, 1.0, 0.0]]) * (alat / 2)
622 elif system['ibrav'] == 3:
623 cell = np.array([[1.0, 1.0, 1.0],
624 [-1.0, 1.0, 1.0],
625 [-1.0, -1.0, 1.0]]) * (alat / 2)
626 elif system['ibrav'] == -3:
627 cell = np.array([[-1.0, 1.0, 1.0],
628 [1.0, -1.0, 1.0],
629 [1.0, 1.0, -1.0]]) * (alat / 2)
630 elif system['ibrav'] == 4:
631 cell = np.array([[1.0, 0.0, 0.0],
632 [-0.5, 0.5 * 3**0.5, 0.0],
633 [0.0, 0.0, c_over_a]]) * alat
634 elif system['ibrav'] == 5:
635 tx = ((1.0 - cosab) / 2.0)**0.5
636 ty = ((1.0 - cosab) / 6.0)**0.5
637 tz = ((1 + 2 * cosab) / 3.0)**0.5
638 cell = np.array([[tx, -ty, tz],
639 [0, 2 * ty, tz],
640 [-tx, -ty, tz]]) * alat
641 elif system['ibrav'] == -5:
642 ty = ((1.0 - cosab) / 6.0)**0.5
643 tz = ((1 + 2 * cosab) / 3.0)**0.5
644 a_prime = alat / 3**0.5
645 u = tz - 2 * 2**0.5 * ty
646 v = tz + 2**0.5 * ty
647 cell = np.array([[u, v, v],
648 [v, u, v],
649 [v, v, u]]) * a_prime
650 elif system['ibrav'] == 6:
651 cell = np.array([[1.0, 0.0, 0.0],
652 [0.0, 1.0, 0.0],
653 [0.0, 0.0, c_over_a]]) * alat
654 elif system['ibrav'] == 7:
655 cell = np.array([[1.0, -1.0, c_over_a],
656 [1.0, 1.0, c_over_a],
657 [-1.0, -1.0, c_over_a]]) * (alat / 2)
658 elif system['ibrav'] == 8:
659 cell = np.array([[1.0, 0.0, 0.0],
660 [0.0, b_over_a, 0.0],
661 [0.0, 0.0, c_over_a]]) * alat
662 elif system['ibrav'] == 9:
663 cell = np.array([[1.0 / 2.0, b_over_a / 2.0, 0.0],
664 [-1.0 / 2.0, b_over_a / 2.0, 0.0],
665 [0.0, 0.0, c_over_a]]) * alat
666 elif system['ibrav'] == -9:
667 cell = np.array([[1.0 / 2.0, -b_over_a / 2.0, 0.0],
668 [1.0 / 2.0, b_over_a / 2.0, 0.0],
669 [0.0, 0.0, c_over_a]]) * alat
670 elif system['ibrav'] == 10:
671 cell = np.array([[1.0 / 2.0, 0.0, c_over_a / 2.0],
672 [1.0 / 2.0, b_over_a / 2.0, 0.0],
673 [0.0, b_over_a / 2.0, c_over_a / 2.0]]) * alat
674 elif system['ibrav'] == 11:
675 cell = np.array([[1.0 / 2.0, b_over_a / 2.0, c_over_a / 2.0],
676 [-1.0 / 2.0, b_over_a / 2.0, c_over_a / 2.0],
677 [-1.0 / 2.0, -b_over_a / 2.0, c_over_a / 2.0]]) * alat
678 elif system['ibrav'] == 12:
679 sinab = (1.0 - cosab**2)**0.5
680 cell = np.array([[1.0, 0.0, 0.0],
681 [b_over_a * cosab, b_over_a * sinab, 0.0],
682 [0.0, 0.0, c_over_a]]) * alat
683 elif system['ibrav'] == -12:
684 sinac = (1.0 - cosac**2)**0.5
685 cell = np.array([[1.0, 0.0, 0.0],
686 [0.0, b_over_a, 0.0],
687 [c_over_a * cosac, 0.0, c_over_a * sinac]]) * alat
688 elif system['ibrav'] == 13:
689 sinab = (1.0 - cosab**2)**0.5
690 cell = np.array([[1.0 / 2.0, 0.0, -c_over_a / 2.0],
691 [b_over_a * cosab, b_over_a * sinab, 0.0],
692 [1.0 / 2.0, 0.0, c_over_a / 2.0]]) * alat
693 elif system['ibrav'] == 14:
694 sinab = (1.0 - cosab**2)**0.5
695 v3 = [c_over_a * cosac,
696 c_over_a * (cosbc - cosac * cosab) / sinab,
697 c_over_a * ((1 + 2 * cosbc * cosac * cosab
698 - cosbc**2 - cosac**2 - cosab**2)**0.5) / sinab]
699 cell = np.array([[1.0, 0.0, 0.0],
700 [b_over_a * cosab, b_over_a * sinab, 0.0],
701 v3]) * alat
702 else:
703 raise NotImplementedError('ibrav = {0} is not implemented'
704 ''.format(system['ibrav']))
706 return Cell(cell)
709def get_pseudo_dirs(data):
710 """Guess a list of possible locations for pseudopotential files.
712 Parameters
713 ----------
714 data : Namelist
715 Namelist representing the quantum espresso input parameters
717 Returns
718 -------
719 pseudo_dirs : list[str]
720 A list of directories where pseudopotential files could be located.
721 """
722 pseudo_dirs = []
723 if 'pseudo_dir' in data['control']:
724 pseudo_dirs.append(data['control']['pseudo_dir'])
725 if 'ESPRESSO_PSEUDO' in os.environ:
726 pseudo_dirs.append(os.environ['ESPRESSO_PSEUDO'])
727 pseudo_dirs.append(path.expanduser('~/espresso/pseudo/'))
728 return pseudo_dirs
731def get_valence_electrons(symbol, data, pseudo=None):
732 """The number of valence electrons for a atomic symbol.
734 Parameters
735 ----------
736 symbol : str
737 Chemical symbol
739 data : Namelist
740 Namelist representing the quantum espresso input parameters
742 pseudo : str, optional
743 File defining the pseudopotential to be used. If missing a fallback
744 to the number of valence electrons recommended at
745 http://materialscloud.org/sssp/ is employed.
746 """
747 if pseudo is None:
748 pseudo = '{}_dummy.UPF'.format(symbol)
749 for pseudo_dir in get_pseudo_dirs(data):
750 if path.exists(path.join(pseudo_dir, pseudo)):
751 valence = grep_valence(path.join(pseudo_dir, pseudo))
752 break
753 else: # not found in a file
754 valence = SSSP_VALENCE[atomic_numbers[symbol]]
755 return valence
758def get_atomic_positions(lines, n_atoms, cell=None, alat=None):
759 """Parse atom positions from ATOMIC_POSITIONS card.
761 Parameters
762 ----------
763 lines : list[str]
764 A list of lines containing the ATOMIC_POSITIONS card.
765 n_atoms : int
766 Expected number of atoms. Only this many lines will be parsed.
767 cell : np.array
768 Unit cell of the crystal. Only used with crystal coordinates.
769 alat : float
770 Lattice parameter for atomic coordinates. Only used for alat case.
772 Returns
773 -------
774 positions : list[(str, (float, float, float), (float, float, float))]
775 A list of the ordered atomic positions in the format:
776 label, (x, y, z), (if_x, if_y, if_z)
777 Force multipliers are set to None if not present.
779 Raises
780 ------
781 ValueError
782 Any problems parsing the data result in ValueError
784 """
786 positions = None
787 # no blanks or comment lines, can the consume n_atoms lines for positions
788 trimmed_lines = (line for line in lines
789 if line.strip() and not line[0] == '#')
791 for line in trimmed_lines:
792 if line.strip().startswith('ATOMIC_POSITIONS'):
793 if positions is not None:
794 raise ValueError('Multiple ATOMIC_POSITIONS specified')
795 # Priority and behaviour tested with QE 5.3
796 if 'crystal_sg' in line.lower():
797 raise NotImplementedError('CRYSTAL_SG not implemented')
798 elif 'crystal' in line.lower():
799 cell = cell
800 elif 'bohr' in line.lower():
801 cell = np.identity(3) * units['Bohr']
802 elif 'angstrom' in line.lower():
803 cell = np.identity(3)
804 # elif 'alat' in line.lower():
805 # cell = np.identity(3) * alat
806 else:
807 if alat is None:
808 raise ValueError('Set lattice parameter in &SYSTEM for '
809 'alat coordinates')
810 # Always the default, will be DEPRECATED as mandatory
811 # in future
812 cell = np.identity(3) * alat
814 positions = []
815 for _dummy in range(n_atoms):
816 split_line = next(trimmed_lines).split()
817 # These can be fractions and other expressions
818 position = np.dot((infix_float(split_line[1]),
819 infix_float(split_line[2]),
820 infix_float(split_line[3])), cell)
821 if len(split_line) > 4:
822 force_mult = (float(split_line[4]),
823 float(split_line[5]),
824 float(split_line[6]))
825 else:
826 force_mult = None
828 positions.append((split_line[0], position, force_mult))
830 return positions
833def get_atomic_species(lines, n_species):
834 """Parse atomic species from ATOMIC_SPECIES card.
836 Parameters
837 ----------
838 lines : list[str]
839 A list of lines containing the ATOMIC_POSITIONS card.
840 n_species : int
841 Expected number of atom types. Only this many lines will be parsed.
843 Returns
844 -------
845 species : list[(str, float, str)]
847 Raises
848 ------
849 ValueError
850 Any problems parsing the data result in ValueError
851 """
853 species = None
854 # no blanks or comment lines, can the consume n_atoms lines for positions
855 trimmed_lines = (line.strip() for line in lines
856 if line.strip() and not line.startswith('#'))
858 for line in trimmed_lines:
859 if line.startswith('ATOMIC_SPECIES'):
860 if species is not None:
861 raise ValueError('Multiple ATOMIC_SPECIES specified')
863 species = []
864 for _dummy in range(n_species):
865 label_weight_pseudo = next(trimmed_lines).split()
866 species.append((label_weight_pseudo[0],
867 float(label_weight_pseudo[1]),
868 label_weight_pseudo[2]))
870 return species
873def get_cell_parameters(lines, alat=None):
874 """Parse unit cell from CELL_PARAMETERS card.
876 Parameters
877 ----------
878 lines : list[str]
879 A list with lines containing the CELL_PARAMETERS card.
880 alat : float | None
881 Unit of lattice vectors in Angstrom. Only used if the card is
882 given in units of alat. alat must be None if CELL_PARAMETERS card
883 is in Bohr or Angstrom. For output files, alat will be parsed from
884 the card header and used in preference to this value.
886 Returns
887 -------
888 cell : np.array | None
889 Cell parameters as a 3x3 array in Angstrom. If no cell is found
890 None will be returned instead.
891 cell_alat : float | None
892 If a value for alat is given in the card header, this is also
893 returned, otherwise this will be None.
895 Raises
896 ------
897 ValueError
898 If CELL_PARAMETERS are given in units of bohr or angstrom
899 and alat is not
900 """
902 cell = None
903 cell_alat = None
904 # no blanks or comment lines, can take three lines for cell
905 trimmed_lines = (line for line in lines
906 if line.strip() and not line[0] == '#')
908 for line in trimmed_lines:
909 if line.strip().startswith('CELL_PARAMETERS'):
910 if cell is not None:
911 # multiple definitions
912 raise ValueError('CELL_PARAMETERS specified multiple times')
913 # Priority and behaviour tested with QE 5.3
914 if 'bohr' in line.lower():
915 if alat is not None:
916 raise ValueError('Lattice parameters given in '
917 '&SYSTEM celldm/A and CELL_PARAMETERS '
918 'bohr')
919 cell_units = units['Bohr']
920 elif 'angstrom' in line.lower():
921 if alat is not None:
922 raise ValueError('Lattice parameters given in '
923 '&SYSTEM celldm/A and CELL_PARAMETERS '
924 'angstrom')
925 cell_units = 1.0
926 elif 'alat' in line.lower():
927 # Output file has (alat = value) (in Bohrs)
928 if '=' in line:
929 alat = float(line.strip(') \n').split()[-1]) * units['Bohr']
930 cell_alat = alat
931 elif alat is None:
932 raise ValueError('Lattice parameters must be set in '
933 '&SYSTEM for alat units')
934 cell_units = alat
935 elif alat is None:
936 # may be DEPRECATED in future
937 cell_units = units['Bohr']
938 else:
939 # may be DEPRECATED in future
940 cell_units = alat
941 # Grab the parameters; blank lines have been removed
942 cell = [[ffloat(x) for x in next(trimmed_lines).split()[:3]],
943 [ffloat(x) for x in next(trimmed_lines).split()[:3]],
944 [ffloat(x) for x in next(trimmed_lines).split()[:3]]]
945 cell = np.array(cell) * cell_units
947 return cell, cell_alat
950def str_to_value(string):
951 """Attempt to convert string into int, float (including fortran double),
952 or bool, in that order, otherwise return the string.
953 Valid (case-insensitive) bool values are: '.true.', '.t.', 'true'
954 and 't' (or false equivalents).
956 Parameters
957 ----------
958 string : str
959 Test to parse for a datatype
961 Returns
962 -------
963 value : any
964 Parsed string as the most appropriate datatype of int, float,
965 bool or string.
967 """
969 # Just an integer
970 try:
971 return int(string)
972 except ValueError:
973 pass
974 # Standard float
975 try:
976 return float(string)
977 except ValueError:
978 pass
979 # Fortran double
980 try:
981 return ffloat(string)
982 except ValueError:
983 pass
985 # possible bool, else just the raw string
986 if string.lower() in ('.true.', '.t.', 'true', 't'):
987 return True
988 elif string.lower() in ('.false.', '.f.', 'false', 'f'):
989 return False
990 else:
991 return string.strip("'")
994def read_fortran_namelist(fileobj):
995 """Takes a fortran-namelist formatted file and returns nested
996 dictionaries of sections and key-value data, followed by a list
997 of lines of text that do not fit the specifications.
999 Behaviour is taken from Quantum ESPRESSO 5.3. Parses fairly
1000 convoluted files the same way that QE should, but may not get
1001 all the MANDATORY rules and edge cases for very non-standard files:
1002 Ignores anything after '!' in a namelist, split pairs on ','
1003 to include multiple key=values on a line, read values on section
1004 start and end lines, section terminating character, '/', can appear
1005 anywhere on a line.
1006 All of these are ignored if the value is in 'quotes'.
1008 Parameters
1009 ----------
1010 fileobj : file
1011 An open file-like object.
1013 Returns
1014 -------
1015 data : dict of dict
1016 Dictionary for each section in the namelist with key = value
1017 pairs of data.
1018 card_lines : list of str
1019 Any lines not used to create the data, assumed to belong to 'cards'
1020 in the input file.
1022 """
1023 # Espresso requires the correct order
1024 data = Namelist()
1025 card_lines = []
1026 in_namelist = False
1027 section = 'none' # can't be in a section without changing this
1029 for line in fileobj:
1030 # leading and trailing whitespace never needed
1031 line = line.strip()
1032 if line.startswith('&'):
1033 # inside a namelist
1034 section = line.split()[0][1:].lower() # case insensitive
1035 if section in data:
1036 # Repeated sections are completely ignored.
1037 # (Note that repeated keys overwrite within a section)
1038 section = "_ignored"
1039 data[section] = Namelist()
1040 in_namelist = True
1041 if not in_namelist and line:
1042 # Stripped line is Truthy, so safe to index first character
1043 if line[0] not in ('!', '#'):
1044 card_lines.append(line)
1045 if in_namelist:
1046 # parse k, v from line:
1047 key = []
1048 value = None
1049 in_quotes = False
1050 for character in line:
1051 if character == ',' and value is not None and not in_quotes:
1052 # finished value:
1053 data[section][''.join(key).strip()] = str_to_value(
1054 ''.join(value).strip())
1055 key = []
1056 value = None
1057 elif character == '=' and value is None and not in_quotes:
1058 # start writing value
1059 value = []
1060 elif character == "'":
1061 # only found in value anyway
1062 in_quotes = not in_quotes
1063 value.append("'")
1064 elif character == '!' and not in_quotes:
1065 break
1066 elif character == '/' and not in_quotes:
1067 in_namelist = False
1068 break
1069 elif value is not None:
1070 value.append(character)
1071 else:
1072 key.append(character)
1073 if value is not None:
1074 data[section][''.join(key).strip()] = str_to_value(
1075 ''.join(value).strip())
1077 return data, card_lines
1080def ffloat(string):
1081 """Parse float from fortran compatible float definitions.
1083 In fortran exponents can be defined with 'd' or 'q' to symbolise
1084 double or quad precision numbers. Double precision numbers are
1085 converted to python floats and quad precision values are interpreted
1086 as numpy longdouble values (platform specific precision).
1088 Parameters
1089 ----------
1090 string : str
1091 A string containing a number in fortran real format
1093 Returns
1094 -------
1095 value : float | np.longdouble
1096 Parsed value of the string.
1098 Raises
1099 ------
1100 ValueError
1101 Unable to parse a float value.
1103 """
1105 if 'q' in string.lower():
1106 return np.longdouble(string.lower().replace('q', 'e'))
1107 else:
1108 return float(string.lower().replace('d', 'e'))
1111def label_to_symbol(label):
1112 """Convert a valid espresso ATOMIC_SPECIES label to a
1113 chemical symbol.
1115 Parameters
1116 ----------
1117 label : str
1118 chemical symbol X (1 or 2 characters, case-insensitive)
1119 or chemical symbol plus a number or a letter, as in
1120 "Xn" (e.g. Fe1) or "X_*" or "X-*" (e.g. C1, C_h;
1121 max total length cannot exceed 3 characters).
1123 Returns
1124 -------
1125 symbol : str
1126 The best matching species from ase.utils.chemical_symbols
1128 Raises
1129 ------
1130 KeyError
1131 Couldn't find an appropriate species.
1133 Notes
1134 -----
1135 It's impossible to tell whether e.g. He is helium
1136 or hydrogen labelled 'e'.
1137 """
1139 # possibly a two character species
1140 # ase Atoms need proper case of chemical symbols.
1141 if len(label) >= 2:
1142 test_symbol = label[0].upper() + label[1].lower()
1143 if test_symbol in chemical_symbols:
1144 return test_symbol
1145 # finally try with one character
1146 test_symbol = label[0].upper()
1147 if test_symbol in chemical_symbols:
1148 return test_symbol
1149 else:
1150 raise KeyError('Could not parse species from label {0}.'
1151 ''.format(label))
1154def infix_float(text):
1155 """Parse simple infix maths into a float for compatibility with
1156 Quantum ESPRESSO ATOMIC_POSITIONS cards. Note: this works with the
1157 example, and most simple expressions, but the capabilities of
1158 the two parsers are not identical. Will also parse a normal float
1159 value properly, but slowly.
1161 >>> infix_float('1/2*3^(-1/2)')
1162 0.28867513459481287
1164 Parameters
1165 ----------
1166 text : str
1167 An arithmetic expression using +, -, *, / and ^, including brackets.
1169 Returns
1170 -------
1171 value : float
1172 Result of the mathematical expression.
1174 """
1176 def middle_brackets(full_text):
1177 """Extract text from innermost brackets."""
1178 start, end = 0, len(full_text)
1179 for (idx, char) in enumerate(full_text):
1180 if char == '(':
1181 start = idx
1182 if char == ')':
1183 end = idx + 1
1184 break
1185 return full_text[start:end]
1187 def eval_no_bracket_expr(full_text):
1188 """Calculate value of a mathematical expression, no brackets."""
1189 exprs = [('+', op.add), ('*', op.mul),
1190 ('/', op.truediv), ('^', op.pow)]
1191 full_text = full_text.lstrip('(').rstrip(')')
1192 try:
1193 return float(full_text)
1194 except ValueError:
1195 for symbol, func in exprs:
1196 if symbol in full_text:
1197 left, right = full_text.split(symbol, 1) # single split
1198 return func(eval_no_bracket_expr(left),
1199 eval_no_bracket_expr(right))
1201 while '(' in text:
1202 middle = middle_brackets(text)
1203 text = text.replace(middle, '{}'.format(eval_no_bracket_expr(middle)))
1205 return float(eval_no_bracket_expr(text))
1207###
1208# Input file writing
1209###
1212# Ordered and case insensitive
1213KEYS = Namelist((
1214 ('CONTROL', [
1215 'calculation', 'title', 'verbosity', 'restart_mode', 'wf_collect',
1216 'nstep', 'iprint', 'tstress', 'tprnfor', 'dt', 'outdir', 'wfcdir',
1217 'prefix', 'lkpoint_dir', 'max_seconds', 'etot_conv_thr',
1218 'forc_conv_thr', 'disk_io', 'pseudo_dir', 'tefield', 'dipfield',
1219 'lelfield', 'nberrycyc', 'lorbm', 'lberry', 'gdir', 'nppstr',
1220 'lfcpopt', 'monopole']),
1221 ('SYSTEM', [
1222 'ibrav', 'nat', 'ntyp', 'nbnd', 'tot_charge', 'tot_magnetization',
1223 'starting_magnetization', 'ecutwfc', 'ecutrho', 'ecutfock', 'nr1',
1224 'nr2', 'nr3', 'nr1s', 'nr2s', 'nr3s', 'nosym', 'nosym_evc', 'noinv',
1225 'no_t_rev', 'force_symmorphic', 'use_all_frac', 'occupations',
1226 'one_atom_occupations', 'starting_spin_angle', 'degauss', 'smearing',
1227 'nspin', 'noncolin', 'ecfixed', 'qcutz', 'q2sigma', 'input_dft',
1228 'exx_fraction', 'screening_parameter', 'exxdiv_treatment',
1229 'x_gamma_extrapolation', 'ecutvcut', 'nqx1', 'nqx2', 'nqx3',
1230 'lda_plus_u', 'lda_plus_u_kind', 'Hubbard_U', 'Hubbard_J0',
1231 'Hubbard_alpha', 'Hubbard_beta', 'Hubbard_J',
1232 'starting_ns_eigenvalue', 'U_projection_type', 'edir',
1233 'emaxpos', 'eopreg', 'eamp', 'angle1', 'angle2',
1234 'constrained_magnetization', 'fixed_magnetization', 'lambda',
1235 'report', 'lspinorb', 'assume_isolated', 'esm_bc', 'esm_w',
1236 'esm_efield', 'esm_nfit', 'fcp_mu', 'vdw_corr', 'london',
1237 'london_s6', 'london_c6', 'london_rvdw', 'london_rcut',
1238 'ts_vdw_econv_thr', 'ts_vdw_isolated', 'xdm', 'xdm_a1', 'xdm_a2',
1239 'space_group', 'uniqueb', 'origin_choice', 'rhombohedral', 'zmon',
1240 'realxz', 'block', 'block_1', 'block_2', 'block_height']),
1241 ('ELECTRONS', [
1242 'electron_maxstep', 'scf_must_converge', 'conv_thr', 'adaptive_thr',
1243 'conv_thr_init', 'conv_thr_multi', 'mixing_mode', 'mixing_beta',
1244 'mixing_ndim', 'mixing_fixed_ns', 'diagonalization', 'ortho_para',
1245 'diago_thr_init', 'diago_cg_maxiter', 'diago_david_ndim',
1246 'diago_full_acc', 'efield', 'efield_cart', 'efield_phase',
1247 'startingpot', 'startingwfc', 'tqr']),
1248 ('IONS', [
1249 'ion_dynamics', 'ion_positions', 'pot_extrapolation',
1250 'wfc_extrapolation', 'remove_rigid_rot', 'ion_temperature', 'tempw',
1251 'tolp', 'delta_t', 'nraise', 'refold_pos', 'upscale', 'bfgs_ndim',
1252 'trust_radius_max', 'trust_radius_min', 'trust_radius_ini', 'w_1',
1253 'w_2']),
1254 ('CELL', [
1255 'cell_dynamics', 'press', 'wmass', 'cell_factor', 'press_conv_thr',
1256 'cell_dofree'])))
1259# Number of valence electrons in the pseudopotentials recommended by
1260# http://materialscloud.org/sssp/. These are just used as a fallback for
1261# calculating initial magetization values which are given as a fraction
1262# of valence electrons.
1263SSSP_VALENCE = [
1264 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,
1265 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,
1266 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,
1267 13.0, 14.0, 15.0, 16.0, 17.0, 18.0, 19.0, 12.0, 13.0, 14.0, 15.0, 6.0,
1268 7.0, 18.0, 9.0, 10.0, 11.0, 12.0, 13.0, 14.0, 15.0, 16.0, 17.0, 18.0,
1269 19.0, 20.0, 21.0, 22.0, 23.0, 24.0, 25.0, 36.0, 27.0, 14.0, 15.0, 30.0,
1270 15.0, 32.0, 19.0, 12.0, 13.0, 14.0, 15.0, 16.0, 18.0]
1273def construct_namelist(parameters=None, warn=False, **kwargs):
1274 """
1275 Construct an ordered Namelist containing all the parameters given (as
1276 a dictionary or kwargs). Keys will be inserted into their appropriate
1277 section in the namelist and the dictionary may contain flat and nested
1278 structures. Any kwargs that match input keys will be incorporated into
1279 their correct section. All matches are case-insensitive, and returned
1280 Namelist object is a case-insensitive dict.
1282 If a key is not known to ase, but in a section within `parameters`,
1283 it will be assumed that it was put there on purpose and included
1284 in the output namelist. Anything not in a section will be ignored (set
1285 `warn` to True to see ignored keys).
1287 Keys with a dimension (e.g. Hubbard_U(1)) will be incorporated as-is
1288 so the `i` should be made to match the output.
1290 The priority of the keys is:
1291 kwargs[key] > parameters[key] > parameters[section][key]
1292 Only the highest priority item will be included.
1294 Parameters
1295 ----------
1296 parameters: dict
1297 Flat or nested set of input parameters.
1298 warn: bool
1299 Enable warnings for unused keys.
1301 Returns
1302 -------
1303 input_namelist: Namelist
1304 pw.x compatible namelist of input parameters.
1306 """
1307 # Convert everything to Namelist early to make case-insensitive
1308 if parameters is None:
1309 parameters = Namelist()
1310 else:
1311 # Maximum one level of nested dict
1312 # Don't modify in place
1313 parameters_namelist = Namelist()
1314 for key, value in parameters.items():
1315 if isinstance(value, dict):
1316 parameters_namelist[key] = Namelist(value)
1317 else:
1318 parameters_namelist[key] = value
1319 parameters = parameters_namelist
1321 # Just a dict
1322 kwargs = Namelist(kwargs)
1324 # Final parameter set
1325 input_namelist = Namelist()
1327 # Collect
1328 for section in KEYS:
1329 sec_list = Namelist()
1330 for key in KEYS[section]:
1331 # Check all three separately and pop them all so that
1332 # we can check for missing values later
1333 if key in parameters.get(section, {}):
1334 sec_list[key] = parameters[section].pop(key)
1335 if key in parameters:
1336 sec_list[key] = parameters.pop(key)
1337 if key in kwargs:
1338 sec_list[key] = kwargs.pop(key)
1340 # Check if there is a key(i) version (no extra parsing)
1341 for arg_key in list(parameters.get(section, {})):
1342 if arg_key.split('(')[0].strip().lower() == key.lower():
1343 sec_list[arg_key] = parameters[section].pop(arg_key)
1344 cp_parameters = parameters.copy()
1345 for arg_key in cp_parameters:
1346 if arg_key.split('(')[0].strip().lower() == key.lower():
1347 sec_list[arg_key] = parameters.pop(arg_key)
1348 cp_kwargs = kwargs.copy()
1349 for arg_key in cp_kwargs:
1350 if arg_key.split('(')[0].strip().lower() == key.lower():
1351 sec_list[arg_key] = kwargs.pop(arg_key)
1353 # Add to output
1354 input_namelist[section] = sec_list
1356 unused_keys = list(kwargs)
1357 # pass anything else already in a section
1358 for key, value in parameters.items():
1359 if key in KEYS and isinstance(value, dict):
1360 input_namelist[key].update(value)
1361 elif isinstance(value, dict):
1362 unused_keys.extend(list(value))
1363 else:
1364 unused_keys.append(key)
1366 if warn and unused_keys:
1367 warnings.warn('Unused keys: {}'.format(', '.join(unused_keys)))
1369 return input_namelist
1372def grep_valence(pseudopotential):
1373 """
1374 Given a UPF pseudopotential file, find the number of valence atoms.
1376 Parameters
1377 ----------
1378 pseudopotential: str
1379 Filename of the pseudopotential.
1381 Returns
1382 -------
1383 valence: float
1384 Valence as reported in the pseudopotential.
1386 Raises
1387 ------
1388 ValueError
1389 If valence cannot be found in the pseudopotential.
1390 """
1392 # Example lines
1393 # Sr.pbe-spn-rrkjus_psl.1.0.0.UPF: z_valence="1.000000000000000E+001"
1394 # C.pbe-n-kjpaw_psl.1.0.0.UPF (new ld1.x):
1395 # ...PBC" z_valence="4.000000000000e0" total_p...
1396 # C_ONCV_PBE-1.0.upf: z_valence=" 4.00"
1397 # Ta_pbe_v1.uspp.F.UPF: 13.00000000000 Z valence
1399 with open(pseudopotential) as psfile:
1400 for line in psfile:
1401 if 'z valence' in line.lower():
1402 return float(line.split()[0])
1403 elif 'z_valence' in line.lower():
1404 if line.split()[0] == '<PP_HEADER':
1405 line = list(filter(lambda x: 'z_valence' in x,
1406 line.split(' ')))[0]
1407 return float(line.split('=')[-1].strip().strip('"'))
1408 else:
1409 raise ValueError('Valence missing in {}'.format(pseudopotential))
1412def kspacing_to_grid(atoms, spacing, calculated_spacing=None):
1413 """
1414 Calculate the kpoint mesh that is equivalent to the given spacing
1415 in reciprocal space (units Angstrom^-1). The number of kpoints is each
1416 dimension is rounded up (compatible with CASTEP).
1418 Parameters
1419 ----------
1420 atoms: ase.Atoms
1421 A structure that can have get_reciprocal_cell called on it.
1422 spacing: float
1423 Minimum K-Point spacing in $A^{-1}$.
1424 calculated_spacing : list
1425 If a three item list (or similar mutable sequence) is given the
1426 members will be replaced with the actual calculated spacing in
1427 $A^{-1}$.
1429 Returns
1430 -------
1431 kpoint_grid : [int, int, int]
1432 MP grid specification to give the required spacing.
1434 """
1435 # No factor of 2pi in ase, everything in A^-1
1436 # reciprocal dimensions
1437 r_x, r_y, r_z = np.linalg.norm(atoms.cell.reciprocal(), axis=1)
1439 kpoint_grid = [int(r_x / spacing) + 1,
1440 int(r_y / spacing) + 1,
1441 int(r_z / spacing) + 1]
1443 for i, _ in enumerate(kpoint_grid):
1444 if not atoms.pbc[i]:
1445 kpoint_grid[i] = 1
1447 if calculated_spacing is not None:
1448 calculated_spacing[:] = [r_x / kpoint_grid[0],
1449 r_y / kpoint_grid[1],
1450 r_z / kpoint_grid[2]]
1452 return kpoint_grid
1455def format_atom_position(atom, crystal_coordinates, mask='', tidx=None):
1456 """Format one line of atomic positions in
1457 Quantum ESPRESSO ATOMIC_POSITIONS card.
1459 >>> for atom in make_supercell(bulk('Li', 'bcc'), np.ones(3)-np.eye(3)):
1460 >>> format_atom_position(atom, True)
1461 Li 0.0000000000 0.0000000000 0.0000000000
1462 Li 0.5000000000 0.5000000000 0.5000000000
1464 Parameters
1465 ----------
1466 atom : Atom
1467 A structure that has symbol and [position | (a, b, c)].
1468 crystal_coordinates: bool
1469 Whether the atomic positions should be written to the QE input file in
1470 absolute (False, default) or relative (crystal) coordinates (True).
1471 mask, optional : str
1472 String of ndim=3 0 or 1 for constraining atomic positions.
1473 tidx, optional : int
1474 Magnetic type index.
1476 Returns
1477 -------
1478 atom_line : str
1479 Input line for atom position
1480 """
1481 if crystal_coordinates:
1482 coords = [atom.a, atom.b, atom.c]
1483 else:
1484 coords = atom.position
1485 line_fmt = '{atom.symbol}'
1486 inps = dict(atom=atom)
1487 if tidx is not None:
1488 line_fmt += '{tidx}'
1489 inps["tidx"] = tidx
1490 line_fmt += ' {coords[0]:.10f} {coords[1]:.10f} {coords[2]:.10f} '
1491 inps["coords"] = coords
1492 line_fmt += ' ' + mask + '\n'
1493 astr = line_fmt.format(**inps)
1494 return astr
1497def write_espresso_in(fd, atoms, input_data=None, pseudopotentials=None,
1498 kspacing=None, kpts=None, koffset=(0, 0, 0),
1499 crystal_coordinates=False, **kwargs):
1500 """
1501 Create an input file for pw.x.
1503 Use set_initial_magnetic_moments to turn on spin, if ispin is set to 2
1504 with no magnetic moments, they will all be set to 0.0. Magnetic moments
1505 will be converted to the QE units (fraction of valence electrons) using
1506 any pseudopotential files found, or a best guess for the number of
1507 valence electrons.
1509 Units are not converted for any other input data, so use Quantum ESPRESSO
1510 units (Usually Ry or atomic units).
1512 Keys with a dimension (e.g. Hubbard_U(1)) will be incorporated as-is
1513 so the `i` should be made to match the output.
1515 Implemented features:
1517 - Conversion of :class:`ase.constraints.FixAtoms` and
1518 :class:`ase.constraints.FixCartesian`.
1519 - `starting_magnetization` derived from the `mgmoms` and pseudopotentials
1520 (searches default paths for pseudo files.)
1521 - Automatic assignment of options to their correct sections.
1523 Not implemented:
1525 - Non-zero values of ibrav
1526 - Lists of k-points
1527 - Other constraints
1528 - Hubbard parameters
1529 - Validation of the argument types for input
1530 - Validation of required options
1531 - Noncollinear magnetism
1533 Parameters
1534 ----------
1535 fd: file
1536 A file like object to write the input file to.
1537 atoms: Atoms
1538 A single atomistic configuration to write to `fd`.
1539 input_data: dict
1540 A flat or nested dictionary with input parameters for pw.x
1541 pseudopotentials: dict
1542 A filename for each atomic species, e.g.
1543 {'O': 'O.pbe-rrkjus.UPF', 'H': 'H.pbe-rrkjus.UPF'}.
1544 A dummy name will be used if none are given.
1545 kspacing: float
1546 Generate a grid of k-points with this as the minimum distance,
1547 in A^-1 between them in reciprocal space. If set to None, kpts
1548 will be used instead.
1549 kpts: (int, int, int) or dict
1550 If kpts is a tuple (or list) of 3 integers, it is interpreted
1551 as the dimensions of a Monkhorst-Pack grid.
1552 If ``kpts`` is set to ``None``, only the Γ-point will be included
1553 and QE will use routines optimized for Γ-point-only calculations.
1554 Compared to Γ-point-only calculations without this optimization
1555 (i.e. with ``kpts=(1, 1, 1)``), the memory and CPU requirements
1556 are typically reduced by half.
1557 If kpts is a dict, it will either be interpreted as a path
1558 in the Brillouin zone (*) if it contains the 'path' keyword,
1559 otherwise it is converted to a Monkhorst-Pack grid (**).
1560 (*) see ase.dft.kpoints.bandpath
1561 (**) see ase.calculators.calculator.kpts2sizeandoffsets
1562 koffset: (int, int, int)
1563 Offset of kpoints in each direction. Must be 0 (no offset) or
1564 1 (half grid offset). Setting to True is equivalent to (1, 1, 1).
1565 crystal_coordinates: bool
1566 Whether the atomic positions should be written to the QE input file in
1567 absolute (False, default) or relative (crystal) coordinates (True).
1569 """
1571 # Convert to a namelist to make working with parameters much easier
1572 # Note that the name ``input_data`` is chosen to prevent clash with
1573 # ``parameters`` in Calculator objects
1574 input_parameters = construct_namelist(input_data, **kwargs)
1576 # Convert ase constraints to QE constraints
1577 # Nx3 array of force multipliers matches what QE uses
1578 # Do this early so it is available when constructing the atoms card
1579 constraint_mask = np.ones((len(atoms), 3), dtype='int')
1580 for constraint in atoms.constraints:
1581 if isinstance(constraint, FixAtoms):
1582 constraint_mask[constraint.index] = 0
1583 elif isinstance(constraint, FixCartesian):
1584 constraint_mask[constraint.a] = constraint.mask
1585 else:
1586 warnings.warn('Ignored unknown constraint {}'.format(constraint))
1587 masks = []
1588 for atom in atoms:
1589 # only inclued mask if something is fixed
1590 if not all(constraint_mask[atom.index]):
1591 mask = ' {mask[0]} {mask[1]} {mask[2]}'.format(
1592 mask=constraint_mask[atom.index])
1593 else:
1594 mask = ''
1595 masks.append(mask)
1597 # Species info holds the information on the pseudopotential and
1598 # associated for each element
1599 if pseudopotentials is None:
1600 pseudopotentials = {}
1601 species_info = {}
1602 for species in set(atoms.get_chemical_symbols()):
1603 # Look in all possible locations for the pseudos and try to figure
1604 # out the number of valence electrons
1605 pseudo = pseudopotentials.get(species, None)
1606 valence = get_valence_electrons(species, input_parameters, pseudo)
1607 species_info[species] = {'pseudo': pseudo, 'valence': valence}
1609 # Convert atoms into species.
1610 # Each different magnetic moment needs to be a separate type even with
1611 # the same pseudopotential (e.g. an up and a down for AFM).
1612 # if any magmom are > 0 or nspin == 2 then use species labels.
1613 # Rememeber: magnetisation uses 1 based indexes
1614 atomic_species = OrderedDict()
1615 atomic_species_str = []
1616 atomic_positions_str = []
1618 nspin = input_parameters['system'].get('nspin', 1) # 1 is the default
1619 if any(atoms.get_initial_magnetic_moments()):
1620 if nspin == 1:
1621 # Force spin on
1622 input_parameters['system']['nspin'] = 2
1623 nspin = 2
1625 if nspin == 2:
1626 # Spin on
1627 for atom, mask, magmom in zip(
1628 atoms, masks, atoms.get_initial_magnetic_moments()):
1629 if (atom.symbol, magmom) not in atomic_species:
1630 # spin as fraction of valence
1631 fspin = float(magmom) / species_info[atom.symbol]['valence']
1632 # Index in the atomic species list
1633 sidx = len(atomic_species) + 1
1634 # Index for that atom type; no index for first one
1635 tidx = sum(atom.symbol == x[0] for x in atomic_species) or ' '
1636 atomic_species[(atom.symbol, magmom)] = (sidx, tidx)
1637 # Add magnetization to the input file
1638 mag_str = 'starting_magnetization({0})'.format(sidx)
1639 input_parameters['system'][mag_str] = fspin
1640 atomic_species_str.append(
1641 '{species}{tidx} {mass} {pseudo}\n'.format(
1642 species=atom.symbol, tidx=tidx, mass=atom.mass,
1643 pseudo=species_info[atom.symbol]['pseudo']))
1644 # lookup tidx to append to name
1645 sidx, tidx = atomic_species[(atom.symbol, magmom)]
1646 # construct line for atomic positions
1647 atomic_positions_str.append(
1648 format_atom_position(
1649 atom, crystal_coordinates, mask=mask, tidx=tidx)
1650 )
1651 else:
1652 # Do nothing about magnetisation
1653 for atom, mask in zip(atoms, masks):
1654 if atom.symbol not in atomic_species:
1655 atomic_species[atom.symbol] = True # just a placeholder
1656 atomic_species_str.append(
1657 '{species} {mass} {pseudo}\n'.format(
1658 species=atom.symbol, mass=atom.mass,
1659 pseudo=species_info[atom.symbol]['pseudo']))
1660 # construct line for atomic positions
1661 atomic_positions_str.append(
1662 format_atom_position(atom, crystal_coordinates, mask=mask)
1663 )
1665 # Add computed parameters
1666 # different magnetisms means different types
1667 input_parameters['system']['ntyp'] = len(atomic_species)
1668 input_parameters['system']['nat'] = len(atoms)
1670 # Use cell as given or fit to a specific ibrav
1671 if 'ibrav' in input_parameters['system']:
1672 ibrav = input_parameters['system']['ibrav']
1673 if ibrav != 0:
1674 raise ValueError(ibrav_error_message)
1675 else:
1676 # Just use standard cell block
1677 input_parameters['system']['ibrav'] = 0
1679 # Construct input file into this
1680 pwi = []
1682 # Assume sections are ordered (taken care of in namelist construction)
1683 # and that repr converts to a QE readable representation (except bools)
1684 for section in input_parameters:
1685 pwi.append('&{0}\n'.format(section.upper()))
1686 for key, value in input_parameters[section].items():
1687 if value is True:
1688 pwi.append(' {0:16} = .true.\n'.format(key))
1689 elif value is False:
1690 pwi.append(' {0:16} = .false.\n'.format(key))
1691 else:
1692 # repr format to get quotes around strings
1693 pwi.append(' {0:16} = {1!r:}\n'.format(key, value))
1694 pwi.append('/\n') # terminate section
1695 pwi.append('\n')
1697 # Pseudopotentials
1698 pwi.append('ATOMIC_SPECIES\n')
1699 pwi.extend(atomic_species_str)
1700 pwi.append('\n')
1702 # KPOINTS - add a MP grid as required
1703 if kspacing is not None:
1704 kgrid = kspacing_to_grid(atoms, kspacing)
1705 elif kpts is not None:
1706 if isinstance(kpts, dict) and 'path' not in kpts:
1707 kgrid, shift = kpts2sizeandoffsets(atoms=atoms, **kpts)
1708 koffset = []
1709 for i, x in enumerate(shift):
1710 assert x == 0 or abs(x * kgrid[i] - 0.5) < 1e-14
1711 koffset.append(0 if x == 0 else 1)
1712 else:
1713 kgrid = kpts
1714 else:
1715 kgrid = "gamma"
1717 # True and False work here and will get converted by ':d' format
1718 if isinstance(koffset, int):
1719 koffset = (koffset, ) * 3
1721 # BandPath object or bandpath-as-dictionary:
1722 if isinstance(kgrid, dict) or hasattr(kgrid, 'kpts'):
1723 pwi.append('K_POINTS crystal_b\n')
1724 assert hasattr(kgrid, 'path') or 'path' in kgrid
1725 kgrid = kpts2ndarray(kgrid, atoms=atoms)
1726 pwi.append('%s\n' % len(kgrid))
1727 for k in kgrid:
1728 pwi.append('{k[0]:.14f} {k[1]:.14f} {k[2]:.14f} 0\n'.format(k=k))
1729 pwi.append('\n')
1730 elif isinstance(kgrid, str) and (kgrid == "gamma"):
1731 pwi.append('K_POINTS gamma\n')
1732 pwi.append('\n')
1733 else:
1734 pwi.append('K_POINTS automatic\n')
1735 pwi.append('{0[0]} {0[1]} {0[2]} {1[0]:d} {1[1]:d} {1[2]:d}\n'
1736 ''.format(kgrid, koffset))
1737 pwi.append('\n')
1739 # CELL block, if required
1740 if input_parameters['SYSTEM']['ibrav'] == 0:
1741 pwi.append('CELL_PARAMETERS angstrom\n')
1742 pwi.append('{cell[0][0]:.14f} {cell[0][1]:.14f} {cell[0][2]:.14f}\n'
1743 '{cell[1][0]:.14f} {cell[1][1]:.14f} {cell[1][2]:.14f}\n'
1744 '{cell[2][0]:.14f} {cell[2][1]:.14f} {cell[2][2]:.14f}\n'
1745 ''.format(cell=atoms.cell))
1746 pwi.append('\n')
1748 # Positions - already constructed, but must appear after namelist
1749 if crystal_coordinates:
1750 pwi.append('ATOMIC_POSITIONS crystal\n')
1751 else:
1752 pwi.append('ATOMIC_POSITIONS angstrom\n')
1753 pwi.extend(atomic_positions_str)
1754 pwi.append('\n')
1756 # DONE!
1757 fd.write(''.join(pwi))