Coverage for /builds/debichem-team/python-ase/ase/calculators/aims.py : 16.76%
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"""This module defines an ASE interface to FHI-aims.
3Felix Hanke hanke@liverpool.ac.uk
4Jonas Bjork j.bjork@liverpool.ac.uk
5Simon P. Rittmeyer simon.rittmeyer@tum.de
6"""
8import os
9import warnings
10import time
11from typing import Optional
12import re
14import numpy as np
16from ase.units import Hartree
17from ase.io.aims import write_aims, read_aims
18from ase.data import atomic_numbers
19from ase.calculators.calculator import FileIOCalculator, Parameters, kpts2mp, \
20 ReadError, PropertyNotImplementedError
23def get_aims_version(string):
24 match = re.search(r'\s*FHI-aims version\s*:\s*(\S+)', string, re.M)
25 return match.group(1)
28float_keys = [
29 'charge',
30 'charge_mix_param',
31 'default_initial_moment',
32 'fixed_spin_moment',
33 'hartree_convergence_parameter',
34 'harmonic_length_scale',
35 'ini_linear_mix_param',
36 'ini_spin_mix_parma',
37 'initial_moment',
38 'MD_MB_init',
39 'MD_time_step',
40 'prec_mix_param',
41 'set_vacuum_level',
42 'spin_mix_param',
43]
45exp_keys = [
46 'basis_threshold',
47 'occupation_thr',
48 'sc_accuracy_eev',
49 'sc_accuracy_etot',
50 'sc_accuracy_forces',
51 'sc_accuracy_rho',
52 'sc_accuracy_stress',
53]
55string_keys = [
56 'communication_type',
57 'density_update_method',
58 'KS_method',
59 'mixer',
60 'output_level',
61 'packed_matrix_format',
62 'relax_unit_cell',
63 'restart',
64 'restart_read_only',
65 'restart_write_only',
66 'spin',
67 'total_energy_method',
68 'qpe_calc',
69 'xc',
70 'species_dir',
71 'run_command',
72 'plus_u',
73]
75int_keys = [
76 'empty_states',
77 'ini_linear_mixing',
78 'max_relaxation_steps',
79 'max_zeroin',
80 'multiplicity',
81 'n_max_pulay',
82 'sc_iter_limit',
83 'walltime',
84]
86bool_keys = [
87 'collect_eigenvectors',
88 'compute_forces',
89 'compute_kinetic',
90 'compute_numerical_stress',
91 'compute_analytical_stress',
92 'compute_heat_flux',
93 'distributed_spline_storage',
94 'evaluate_work_function',
95 'final_forces_cleaned',
96 'hessian_to_restart_geometry',
97 'load_balancing',
98 'MD_clean_rotations',
99 'MD_restart',
100 'override_illconditioning',
101 'override_relativity',
102 'restart_relaxations',
103 'squeeze_memory',
104 'symmetry_reduced_k_grid',
105 'use_density_matrix',
106 'use_dipole_correction',
107 'use_local_index',
108 'use_logsbt',
109 'vdw_correction_hirshfeld',
110]
112list_keys = [
113 'init_hess',
114 'k_grid',
115 'k_offset',
116 'MD_run',
117 'MD_schedule',
118 'MD_segment',
119 'mixer_threshold',
120 'occupation_type',
121 'output',
122 'cube',
123 'preconditioner',
124 'relativistic',
125 'relax_geometry',
126]
129class Aims(FileIOCalculator):
130 # was "command" before the refactoring to dynamical commands
131 __command_default = 'aims.version.serial.x > aims.out'
132 __outfilename_default = 'aims.out'
134 implemented_properties = ['energy', 'forces', 'stress', 'stresses',
135 'dipole', 'magmom']
137 def __init__(self, restart=None,
138 ignore_bad_restart_file=FileIOCalculator._deprecated,
139 label=os.curdir, atoms=None, cubes=None, radmul=None,
140 tier=None, aims_command=None,
141 outfilename=None, **kwargs):
142 """Construct the FHI-aims calculator.
144 The keyword arguments (kwargs) can be one of the ASE standard
145 keywords: 'xc', 'kpts' and 'smearing' or any of FHI-aims'
146 native keywords.
148 .. note:: The behavior of command/run_command has been refactored ase X.X.X
149 It is now possible to independently specify the command to call
150 FHI-aims and the outputfile into which stdout is directed. In
151 general, we replaced
153 <run_command> = <aims_command> + " > " + <outfilename
155 That is,what used to be, e.g.,
157 >>> calc = Aims(run_command = "mpiexec -np 4 aims.x > aims.out")
159 can now be achieved with the two arguments
161 >>> calc = Aims(aims_command = "mpiexec -np 4 aims.x"
162 >>> outfilename = "aims.out")
164 Backward compatibility, however, is provided. Also, the command
165 actually used to run FHI-aims is dynamically updated (i.e., the
166 "command" member variable). That is, e.g.,
168 >>> calc = Aims()
169 >>> print(calc.command)
170 aims.version.serial.x > aims.out
171 >>> calc.outfilename = "systemX.out"
172 >>> print(calc.command)
173 aims.version.serial.x > systemX.out
174 >>> calc.aims_command = "mpiexec -np 4 aims.version.scalapack.mpi.x"
175 >>> print(calc.command)
176 mpiexec -np 4 aims.version.scalapack.mpi > systemX.out
179 Arguments:
181 cubes: AimsCube object
182 Cube file specification.
184 radmul: int
185 Set radial multiplier for the basis set of all atomic species.
187 tier: int or array of ints
188 Set basis set tier for all atomic species.
190 aims_command : str
191 The full command as executed to run FHI-aims *without* the
192 redirection to stdout. For instance "mpiexec -np 4 aims.x". Note
193 that this is not the same as "command" or "run_command".
194 .. note:: Added in ase X.X.X
196 outfilename : str
197 The file (incl. path) to which stdout is redirected. Defaults to
198 "aims.out"
199 .. note:: Added in ase X.X.X
201 run_command : str, optional (default=None)
202 Same as "command", see FileIOCalculator documentation.
203 .. note:: Deprecated in ase X.X.X
205 outfilename : str, optional (default=aims.out)
206 File into which the stdout of the FHI aims run is piped into. Note
207 that this will be only of any effect, if the <run_command> does not
208 yet contain a '>' directive.
209 plus_u : dict
210 For DFT+U. Adds a +U term to one specific shell of the species.
212 kwargs : dict
213 Any of the base class arguments.
215 """
216 # yes, we pop the key and run it through our legacy filters
217 command = kwargs.pop('command', None)
219 # Check for the "run_command" (deprecated keyword)
220 # Consistently, the "command" argument should be used as suggested by the FileIO base class.
221 # For legacy reasons, however, we here also accept "run_command"
222 run_command = kwargs.pop('run_command', None)
223 if run_command:
224 # this warning is debatable... in my eyes it is more consistent to
225 # use 'command'
226 warnings.warn('Argument "run_command" is deprecated and will be replaced with "command". Alternatively, use "aims_command" and "outfile". See documentation for more details.')
227 if command:
228 warnings.warn('Caution! Argument "command" overwrites "run_command.')
229 else:
230 command = run_command
232 # this is the fallback to the default value for empty init
233 if np.all([i is None for i in (command, aims_command, outfilename)]):
234 # we go for the FileIOCalculator default way (env variable) with the former default as fallback
235 command = os.environ.get('ASE_AIMS_COMMAND', Aims.__command_default)
237 # filter the command and set the member variables "aims_command" and "outfilename"
238 self.__init_command(command=command,
239 aims_command=aims_command,
240 outfilename=outfilename)
242 FileIOCalculator.__init__(self, restart, ignore_bad_restart_file,
243 label, atoms,
244 # well, this is not nice, but cannot work around it...
245 command=self.command,
246 **kwargs)
248 self.cubes = cubes
249 self.radmul = radmul
250 self.tier = tier
252 # handling the filtering for dynamical commands with properties,
253 @property # type: ignore
254 def command(self) -> Optional[str]: # type: ignore
255 return self.__command
257 @command.setter
258 def command(self, x):
259 self.__update_command(command=x)
261 @property
262 def aims_command(self):
263 return self.__aims_command
265 @aims_command.setter
266 def aims_command(self, x):
267 self.__update_command(aims_command=x)
269 @property
270 def outfilename(self):
271 return self.__outfilename
273 @outfilename.setter
274 def outfilename(self, x):
275 self.__update_command(outfilename=x)
277 def __init_command(self, command=None, aims_command=None,
278 outfilename=None):
279 """
280 Create the private variables for which properties are defines and set
281 them accordingly.
282 """
283 # new class variables due to dynamical command handling
284 self.__aims_command = None
285 self.__outfilename = None
286 self.__command: Optional[str] = None
288 # filter the command and set the member variables "aims_command" and "outfilename"
289 self.__update_command(command=command,
290 aims_command=aims_command,
291 outfilename=outfilename)
293 # legacy handling of the (run_)command behavior a.k.a. a universal setter routine
294 def __update_command(self, command=None, aims_command=None,
295 outfilename=None):
296 """
297 Abstracted generic setter routine for a dynamic behavior of "command".
299 The command that is actually called on the command line and enters the
300 base class, is <command> = <aims_command> > <outfilename>.
302 This new scheme has been introduced in order to conveniently change the
303 outfile name from the outside while automatically updating the
304 <command> member variable.
306 Obiously, changing <command> conflicts with changing <aims_command>
307 and/or <outfilename>, which thus raises a <ValueError>. This should,
308 however, not happen if this routine is not used outside the property
309 definitions.
311 Parameters
312 ----------
313 command : str
314 The full command as executed to run FHI-aims. This includes
315 any potential mpiexec call, as well as the redirection of stdout.
316 For instance "mpiexec -np 4 aims.x > aims.out".
318 aims_command : str
319 The full command as executed to run FHI-aims *without* the
320 redirection to stdout. For instance "mpiexec -np 4 aims.x"
322 outfilename : str
323 The file (incl. path) to which stdout is redirected.
324 """
325 # disentangle the command if given
326 if command:
327 if aims_command:
328 raise ValueError('Cannot specify "command" and "aims_command" simultaneously.')
329 if outfilename:
330 raise ValueError('Cannot specify "command" and "outfilename" simultaneously.')
332 # check if the redirection of stdout is included
333 command_spl = command.split('>')
334 if len(command_spl) > 1:
335 self.__aims_command = command_spl[0].strip()
336 self.__outfilename = command_spl[-1].strip()
337 else:
338 # this should not happen if used correctly
339 # but just to ensure legacy behavior of how "run_command" was handled
340 self.__aims_command = command.strip()
341 self.__outfilename = Aims.__outfilename_default
342 else:
343 if aims_command is not None:
344 self.__aims_command = aims_command
345 elif outfilename is None:
346 # nothing to do here, empty call with 3x None
347 return
348 if outfilename is not None:
349 self.__outfilename = outfilename
350 else:
351 # default to 'aims.out'
352 if not self.outfilename:
353 self.__outfilename = Aims.__outfilename_default
355 self.__command = '{0:s} > {1:s}'.format(self.aims_command,
356 self.outfilename)
358 def set_atoms(self, atoms):
359 self.atoms = atoms
361 def set_label(self, label, update_outfilename=False):
362 msg = "Aims.set_label is not supported anymore, please use `directory`"
363 raise RuntimeError(msg)
365 @property
366 def out(self):
367 return os.path.join(self.label, self.outfilename)
369 def check_state(self, atoms):
370 system_changes = FileIOCalculator.check_state(self, atoms)
371 # Ignore unit cell for molecules:
372 if not atoms.pbc.any() and 'cell' in system_changes:
373 system_changes.remove('cell')
374 return system_changes
376 def set(self, **kwargs):
377 xc = kwargs.get('xc')
378 if xc:
379 kwargs['xc'] = {'LDA': 'pw-lda', 'PBE': 'pbe'}.get(xc, xc)
381 changed_parameters = FileIOCalculator.set(self, **kwargs)
383 if changed_parameters:
384 self.reset()
385 return changed_parameters
387 def write_input(self, atoms, properties=None, system_changes=None,
388 ghosts=None, geo_constrain=None, scaled=None, velocities=None):
389 FileIOCalculator.write_input(self, atoms, properties, system_changes)
391 if geo_constrain is None:
392 geo_constrain = "relax_geometry" in self.parameters
394 if scaled is None:
395 scaled = np.all(atoms.get_pbc())
396 if velocities is None:
397 velocities = atoms.has('momenta')
399 have_lattice_vectors = atoms.pbc.any()
400 have_k_grid = ('k_grid' in self.parameters or
401 'kpts' in self.parameters)
402 if have_lattice_vectors and not have_k_grid:
403 raise RuntimeError('Found lattice vectors but no k-grid!')
404 if not have_lattice_vectors and have_k_grid:
405 raise RuntimeError('Found k-grid but no lattice vectors!')
406 write_aims(
407 os.path.join(self.directory, 'geometry.in'),
408 atoms,
409 scaled,
410 geo_constrain,
411 velocities=velocities,
412 ghosts=ghosts
413 )
414 self.write_control(atoms, os.path.join(self.directory, 'control.in'))
415 self.write_species(atoms, os.path.join(self.directory, 'control.in'))
416 self.parameters.write(os.path.join(self.directory, 'parameters.ase'))
418 def prepare_input_files(self):
419 """
420 Wrapper function to prepare input filesi, e.g., to a run on a remote
421 machine
422 """
423 if self.atoms is None:
424 raise ValueError('No atoms object attached')
425 self.write_input(self.atoms)
427 def write_control(self, atoms, filename, debug=False):
428 lim = '#' + '='*79
429 output = open(filename, 'w')
430 output.write(lim + '\n')
431 for line in ['FHI-aims file: ' + filename,
432 'Created using the Atomic Simulation Environment (ASE)',
433 time.asctime(),
434 ]:
435 output.write('# ' + line + '\n')
436 if debug:
437 output.write('# \n# List of parameters used to initialize the calculator:',)
438 for p, v in self.parameters.items():
439 s = '# {} : {}\n'.format(p, v)
440 output.write(s)
441 output.write(lim + '\n')
443 assert not ('kpts' in self.parameters and 'k_grid' in self.parameters)
444 assert not ('smearing' in self.parameters and
445 'occupation_type' in self.parameters)
447 for key, value in self.parameters.items():
448 if key == 'kpts':
449 mp = kpts2mp(atoms, self.parameters.kpts)
450 output.write('%-35s%d %d %d\n' % (('k_grid',) + tuple(mp)))
451 dk = 0.5 - 0.5 / np.array(mp)
452 output.write('%-35s%f %f %f\n' % (('k_offset',) + tuple(dk)))
453 elif key == 'species_dir' or key == 'run_command':
454 continue
455 elif key == 'plus_u':
456 continue
457 elif key == 'smearing':
458 name = self.parameters.smearing[0].lower()
459 if name == 'fermi-dirac':
460 name = 'fermi'
461 width = self.parameters.smearing[1]
462 output.write('%-35s%s %f' % ('occupation_type', name, width))
463 if name == 'methfessel-paxton':
464 order = self.parameters.smearing[2]
465 output.write(' %d' % order)
466 output.write('\n' % order)
467 elif key == 'output':
468 for output_type in value:
469 output.write('%-35s%s\n' % (key, output_type))
470 elif key == 'vdw_correction_hirshfeld' and value:
471 output.write('%-35s\n' % key)
472 elif key in bool_keys:
473 output.write('%-35s.%s.\n' % (key, repr(bool(value)).lower()))
474 elif isinstance(value, (tuple, list)):
475 output.write('%-35s%s\n' %
476 (key, ' '.join(str(x) for x in value)))
477 elif isinstance(value, str):
478 output.write('%-35s%s\n' % (key, value))
479 else:
480 output.write('%-35s%r\n' % (key, value))
481 if self.cubes:
482 self.cubes.write(output)
483 output.write(lim + '\n\n')
484 output.close()
486 def read(self, label=None):
487 if label is None:
488 label = self.label
489 FileIOCalculator.read(self, label)
490 geometry = os.path.join(self.directory, 'geometry.in')
491 control = os.path.join(self.directory, 'control.in')
493 for filename in [geometry, control, self.out]:
494 if not os.path.isfile(filename):
495 raise ReadError
497 self.atoms, symmetry_block = read_aims(geometry, True)
498 self.parameters = Parameters.read(os.path.join(self.directory,
499 'parameters.ase'))
500 if symmetry_block:
501 self.parameters["symmetry_block"] = symmetry_block
502 self.read_results()
504 def read_results(self):
505 converged = self.read_convergence()
506 if not converged:
507 os.system('tail -20 ' + self.out)
508 raise RuntimeError('FHI-aims did not converge!\n' +
509 'The last lines of output are printed above ' +
510 'and should give an indication why.')
511 self.read_energy()
512 if ('compute_forces' in self.parameters or
513 'sc_accuracy_forces' in self.parameters):
514 self.read_forces()
516 if ('sc_accuracy_stress' in self.parameters or
517 ('compute_numerical_stress' in self.parameters
518 and self.parameters['compute_numerical_stress']) or
519 ('compute_analytical_stress' in self.parameters
520 and self.parameters['compute_analytical_stress']) or
521 ('compute_heat_flux' in self.parameters
522 and self.parameters['compute_heat_flux'])):
523 self.read_stress()
525 if ('compute_heat_flux' in self.parameters
526 and self.parameters['compute_heat_flux']):
527 self.read_stresses()
529 if ('dipole' in self.parameters.get('output', []) and
530 not self.atoms.pbc.any()):
531 self.read_dipole()
533 def write_species(self, atoms, filename):
534 self.ctrlname = filename
535 species_path = self.parameters.get('species_dir')
536 if species_path is None:
537 species_path = os.environ.get('AIMS_SPECIES_DIR')
538 if species_path is None:
539 raise RuntimeError(
540 'Missing species directory! Use species_dir ' +
541 'parameter or set $AIMS_SPECIES_DIR environment variable.')
542 control = open(filename, 'a')
543 symbols = atoms.get_chemical_symbols()
544 symbols2 = []
545 for n, symbol in enumerate(symbols):
546 if symbol not in symbols2:
547 symbols2.append(symbol)
548 if self.tier is not None:
549 if isinstance(self.tier, int):
550 self.tierlist = np.ones(len(symbols2), 'int') * self.tier
551 elif isinstance(self.tier, list):
552 assert len(self.tier) == len(symbols2)
553 self.tierlist = self.tier
555 for i, symbol in enumerate(symbols2):
556 fd = os.path.join(species_path, '%02i_%s_default' %
557 (atomic_numbers[symbol], symbol))
558 reached_tiers = False
559 for line in open(fd, 'r'):
560 if self.tier is not None:
561 if 'First tier' in line:
562 reached_tiers = True
563 self.targettier = self.tierlist[i]
564 self.foundtarget = False
565 self.do_uncomment = True
566 if reached_tiers:
567 line = self.format_tiers(line)
568 control.write(line)
569 if self.tier is not None and not self.foundtarget:
570 raise RuntimeError(
571 "Basis tier %i not found for element %s" %
572 (self.targettier, symbol))
573 if self.parameters.get('plus_u') is not None:
574 if symbol in self.parameters.plus_u.keys():
575 control.write('plus_u %s \n' %
576 self.parameters.plus_u[symbol])
577 control.close()
579 if self.radmul is not None:
580 self.set_radial_multiplier()
582 def format_tiers(self, line):
583 if 'meV' in line:
584 assert line[0] == '#'
585 if 'tier' in line and 'Further' not in line:
586 tier = line.split(" tier")[0]
587 tier = tier.split('"')[-1]
588 current_tier = self.translate_tier(tier)
589 if current_tier == self.targettier:
590 self.foundtarget = True
591 elif current_tier > self.targettier:
592 self.do_uncomment = False
593 else:
594 self.do_uncomment = False
595 return line
596 elif self.do_uncomment and line[0] == '#':
597 return line[1:]
598 elif not self.do_uncomment and line[0] != '#':
599 return '#' + line
600 else:
601 return line
603 def translate_tier(self, tier):
604 if tier.lower() == 'first':
605 return 1
606 elif tier.lower() == 'second':
607 return 2
608 elif tier.lower() == 'third':
609 return 3
610 elif tier.lower() == 'fourth':
611 return 4
612 else:
613 return -1
615 def set_radial_multiplier(self):
616 assert isinstance(self.radmul, int)
617 newctrl = self.ctrlname + '.new'
618 fin = open(self.ctrlname, 'r')
619 fout = open(newctrl, 'w')
620 newline = " radial_multiplier %i\n" % self.radmul
621 for line in fin:
622 if ' radial_multiplier' in line:
623 fout.write(newline)
624 else:
625 fout.write(line)
626 fin.close()
627 fout.close()
628 os.rename(newctrl, self.ctrlname)
630 def get_dipole_moment(self, atoms):
631 if ('dipole' not in self.parameters.get('output', []) or
632 atoms.pbc.any()):
633 raise PropertyNotImplementedError
634 return FileIOCalculator.get_dipole_moment(self, atoms)
636 def get_stress(self, atoms):
637 if ('compute_numerical_stress' not in self.parameters and
638 'compute_analytical_stress' not in self.parameters):
639 raise PropertyNotImplementedError
640 return FileIOCalculator.get_stress(self, atoms)
642 def get_forces(self, atoms):
643 if ('compute_forces' not in self.parameters and
644 'sc_accuracy_forces' not in self.parameters):
645 raise PropertyNotImplementedError
646 return FileIOCalculator.get_forces(self, atoms)
648 def read_dipole(self):
649 "Method that reads the electric dipole moment from the output file."
650 for line in open(self.out, 'r'):
651 if line.rfind('Total dipole moment [eAng]') > -1:
652 dipolemoment = np.array([float(f)
653 for f in line.split()[6:9]])
654 self.results['dipole'] = dipolemoment
656 def read_energy(self):
657 for line in open(self.out, 'r'):
658 if line.rfind('Total energy corrected') > -1:
659 E0 = float(line.split()[5])
660 elif line.rfind('Total energy uncorrected') > -1:
661 F = float(line.split()[5])
662 self.results['free_energy'] = F
663 self.results['energy'] = E0
665 def read_forces(self):
666 """Method that reads forces from the output file.
668 If 'all' is switched on, the forces for all ionic steps
669 in the output file will be returned, in other case only the
670 forces for the last ionic configuration are returned."""
671 lines = open(self.out, 'r').readlines()
672 forces = np.zeros([len(self.atoms), 3])
673 for n, line in enumerate(lines):
674 if line.rfind('Total atomic forces') > -1:
675 for iatom in range(len(self.atoms)):
676 data = lines[n + iatom + 1].split()
677 for iforce in range(3):
678 forces[iatom, iforce] = float(data[2 + iforce])
679 self.results['forces'] = forces
681 def read_stress(self):
682 lines = open(self.out, 'r').readlines()
683 stress = None
684 for n, line in enumerate(lines):
685 if (line.rfind('| Analytical stress tensor') > -1 or
686 line.rfind('Numerical stress tensor') > -1):
687 stress = []
688 for i in [n + 5, n + 6, n + 7]:
689 data = lines[i].split()
690 stress += [float(data[2]), float(data[3]), float(data[4])]
691 # rearrange in 6-component form and return
692 self.results['stress'] = np.array([stress[0], stress[4], stress[8],
693 stress[5], stress[2], stress[1]])
695 def read_stresses(self):
696 """ Read stress per atom """
697 with open(self.out) as fd:
698 next(l for l in fd if
699 'Per atom stress (eV) used for heat flux calculation' in l)
700 # scroll to boundary
701 next(l for l in fd if '-------------' in l)
703 stresses = []
704 for l in [next(fd) for _ in range(len(self.atoms))]:
705 # Read stresses and rearrange from
706 # (xx, yy, zz, xy, xz, yz) to (xx, yy, zz, yz, xz, xy)
707 xx, yy, zz, xy, xz, yz = [float(d) for d in l.split()[2:8]]
708 stresses.append([xx, yy, zz, yz, xz, xy])
710 self.results['stresses'] = np.array(stresses)
712 def get_stresses(self, voigt=False):
713 """ Return stress per atom
715 Returns an array of the six independent components of the
716 symmetric stress tensor per atom, in the traditional Voigt order
717 (xx, yy, zz, yz, xz, xy) or as a 3x3 matrix. Default is 3x3 matrix.
718 """
720 voigt_stresses = self.results['stresses']
722 if voigt:
723 return voigt_stresses
724 else:
725 stresses = np.zeros((len(self.atoms), 3, 3))
726 for ii, stress in enumerate(voigt_stresses):
727 xx, yy, zz, yz, xz, xy = stress
728 stresses[ii] = np.array([(xx, xy, xz),
729 (xy, yy, yz),
730 (xz, yz, zz)])
731 return stresses
733 def read_convergence(self):
734 converged = False
735 lines = open(self.out, 'r').readlines()
736 for n, line in enumerate(lines):
737 if line.rfind('Have a nice day') > -1:
738 converged = True
739 return converged
741 def get_number_of_iterations(self):
742 return self.read_number_of_iterations()
744 def read_number_of_iterations(self):
745 niter = None
746 lines = open(self.out, 'r').readlines()
747 for n, line in enumerate(lines):
748 if line.rfind('| Number of self-consistency cycles') > -1:
749 niter = int(line.split(':')[-1].strip())
750 return niter
752 def get_electronic_temperature(self):
753 return self.read_electronic_temperature()
755 def read_electronic_temperature(self):
756 width = None
757 lines = open(self.out, 'r').readlines()
758 for n, line in enumerate(lines):
759 if line.rfind('Occupation type:') > -1:
760 width = float(line.split('=')[-1].strip().split()[0])
761 return width
763 def get_number_of_electrons(self):
764 return self.read_number_of_electrons()
766 def read_number_of_electrons(self):
767 nelect = None
768 lines = open(self.out, 'r').readlines()
769 for n, line in enumerate(lines):
770 if line.rfind('The structure contains') > -1:
771 nelect = float(line.split()[-2].strip())
772 return nelect
774 def get_number_of_bands(self):
775 return self.read_number_of_bands()
777 def read_number_of_bands(self):
778 nband = None
779 lines = open(self.out, 'r').readlines()
780 for n, line in enumerate(lines):
781 if line.rfind('Number of Kohn-Sham states') > -1:
782 nband = int(line.split(':')[-1].strip())
783 return nband
785 def get_k_point_weights(self):
786 return self.read_kpts(mode='k_point_weights')
788 def get_bz_k_points(self):
789 raise NotImplementedError
791 def get_ibz_k_points(self):
792 return self.read_kpts(mode='ibz_k_points')
794 def get_spin_polarized(self):
795 return self.read_number_of_spins()
797 def get_number_of_spins(self):
798 return 1 + self.get_spin_polarized()
800 def get_magnetic_moment(self, atoms=None):
801 return self.read_magnetic_moment()
803 def read_number_of_spins(self):
804 spinpol = None
805 lines = open(self.out, 'r').readlines()
806 for n, line in enumerate(lines):
807 if line.rfind('| Number of spin channels') > -1:
808 spinpol = int(line.split(':')[-1].strip()) - 1
809 return spinpol
811 def read_magnetic_moment(self):
812 magmom = None
813 if not self.get_spin_polarized():
814 magmom = 0.0
815 else: # only for spinpolarized system Magnetisation is printed
816 for line in open(self.out, 'r').readlines():
817 if line.find('N_up - N_down') != -1: # last one
818 magmom = float(line.split(':')[-1].strip())
819 return magmom
821 def get_fermi_level(self):
822 return self.read_fermi()
824 def get_eigenvalues(self, kpt=0, spin=0):
825 return self.read_eigenvalues(kpt, spin, 'eigenvalues')
827 def get_occupations(self, kpt=0, spin=0):
828 return self.read_eigenvalues(kpt, spin, 'occupations')
830 def read_fermi(self):
831 E_f = None
832 lines = open(self.out, 'r').readlines()
833 for n, line in enumerate(lines):
834 if line.rfind('| Chemical potential (Fermi level) in eV') > -1:
835 E_f = float(line.split(':')[-1].strip())
836 return E_f
838 def read_kpts(self, mode='ibz_k_points'):
839 """ Returns list of kpts weights or kpts coordinates. """
840 values = []
841 assert mode in ['ibz_k_points', 'k_point_weights']
842 lines = open(self.out, 'r').readlines()
843 kpts = None
844 kptsstart = None
845 for n, line in enumerate(lines):
846 if line.rfind('| Number of k-points') > -1:
847 kpts = int(line.split(':')[-1].strip())
848 for n, line in enumerate(lines):
849 if line.rfind('K-points in task') > -1:
850 kptsstart = n # last occurrence of (
851 assert kpts is not None
852 assert kptsstart is not None
853 text = lines[kptsstart + 1:]
854 values = []
855 for line in text[:kpts]:
856 if mode == 'ibz_k_points':
857 b = [float(c.strip()) for c in line.split()[4:7]]
858 else:
859 b = float(line.split()[-1])
860 values.append(b)
861 if len(values) == 0:
862 values = None
863 return np.array(values)
865 def read_eigenvalues(self, kpt=0, spin=0, mode='eigenvalues'):
866 """ Returns list of last eigenvalues, occupations
867 for given kpt and spin. """
868 values = []
869 assert mode in ['eigenvalues', 'occupations']
870 lines = open(self.out, 'r').readlines()
871 # number of kpts
872 kpts = None
873 for n, line in enumerate(lines):
874 if line.rfind('| Number of k-points') > -1:
875 kpts = int(line.split(':')[-1].strip())
876 break
877 assert kpts is not None
878 assert kpt + 1 <= kpts
879 # find last (eigenvalues)
880 eigvalstart = None
881 for n, line in enumerate(lines):
882 # eigenvalues come after Preliminary charge convergence reached
883 if line.rfind('Preliminary charge convergence reached') > -1:
884 eigvalstart = n
885 break
886 assert eigvalstart is not None
887 lines = lines[eigvalstart:]
888 for n, line in enumerate(lines):
889 if line.rfind('Writing Kohn-Sham eigenvalues') > -1:
890 eigvalstart = n
891 break
892 assert eigvalstart is not None
893 text = lines[eigvalstart + 1:] # remove first 1 line
894 # find the requested k-point
895 nbands = self.read_number_of_bands()
896 sppol = self.get_spin_polarized()
897 beg = ((nbands + 4 + int(sppol) * 1) * kpt * (sppol + 1) +
898 3 + sppol * 2 + kpt * sppol)
899 if self.get_spin_polarized():
900 if spin == 0:
901 beg = beg
902 end = beg + nbands
903 else:
904 beg = beg + nbands + 5
905 end = beg + nbands
906 else:
907 end = beg + nbands
908 values = []
909 for line in text[beg:end]:
910 # aims prints stars for large values ...
911 line = line.replace('**************', ' 10000')
912 line = line.replace('***************', ' 10000')
913 line = line.replace('****************', ' 10000')
914 b = [float(c.strip()) for c in line.split()[1:]]
915 values.append(b)
916 if mode == 'eigenvalues':
917 values = [Hartree * v[1] for v in values]
918 else:
919 values = [v[0] for v in values]
920 if len(values) == 0:
921 values = None
922 return np.array(values)
925class AimsCube:
926 "Object to ensure the output of cube files, can be attached to Aims object"
927 def __init__(self, origin=(0, 0, 0),
928 edges=[(0.1, 0.0, 0.0), (0.0, 0.1, 0.0), (0.0, 0.0, 0.1)],
929 points=(50, 50, 50), plots=None):
930 """parameters:
932 origin, edges, points:
933 Same as in the FHI-aims output
934 plots:
935 what to print, same names as in FHI-aims """
937 self.name = 'AimsCube'
938 self.origin = origin
939 self.edges = edges
940 self.points = points
941 self.plots = plots
943 def ncubes(self):
944 """returns the number of cube files to output """
945 if self.plots:
946 number = len(self.plots)
947 else:
948 number = 0
949 return number
951 def set(self, **kwargs):
952 """ set any of the parameters ... """
953 # NOT IMPLEMENTED AT THE MOMENT!
955 def move_to_base_name(self, basename):
956 """ when output tracking is on or the base namem is not standard,
957 this routine will rename add the base to the cube file output for
958 easier tracking """
959 for plot in self.plots:
960 found = False
961 cube = plot.split()
962 if (cube[0] == 'total_density' or
963 cube[0] == 'spin_density' or
964 cube[0] == 'delta_density'):
965 found = True
966 old_name = cube[0] + '.cube'
967 new_name = basename + '.' + old_name
968 if cube[0] == 'eigenstate' or cube[0] == 'eigenstate_density':
969 found = True
970 state = int(cube[1])
971 s_state = cube[1]
972 for i in [10, 100, 1000, 10000]:
973 if state < i:
974 s_state = '0' + s_state
975 old_name = cube[0] + '_' + s_state + '_spin_1.cube'
976 new_name = basename + '.' + old_name
977 if found:
978 os.system('mv ' + old_name + ' ' + new_name)
980 def add_plot(self, name):
981 """ in case you forgot one ... """
982 self.plots += [name]
984 def write(self, file):
985 """ write the necessary output to the already opened control.in """
986 file.write('output cube ' + self.plots[0] + '\n')
987 file.write(' cube origin ')
988 for ival in self.origin:
989 file.write(str(ival) + ' ')
990 file.write('\n')
991 for i in range(3):
992 file.write(' cube edge ' + str(self.points[i]) + ' ')
993 for ival in self.edges[i]:
994 file.write(str(ival) + ' ')
995 file.write('\n')
996 if self.ncubes() > 1:
997 for i in range(self.ncubes() - 1):
998 file.write('output cube ' + self.plots[i + 1] + '\n')