Coverage for /builds/kinetik161/ase/ase/calculators/siesta/siesta.py: 62.80%
656 statements
« prev ^ index » next coverage.py v7.2.7, created at 2023-12-10 11:04 +0000
1"""
2This module defines the ASE interface to SIESTA.
4Written by Mads Engelund (see www.espeem.com)
6Home of the SIESTA package:
7http://www.uam.es/departamentos/ciencias/fismateriac/siesta
92017.04 - Pedro Brandimarte: changes for python 2-3 compatible
11"""
13import os
14import re
15import shutil
16import tempfile
17from typing import Any, Dict, List
18import warnings
19from os.path import isfile, islink, join
21import numpy as np
23from ase import Atoms
24from ase.calculators.calculator import (FileIOCalculator, Parameters,
25 ReadError, all_changes)
26from ase.calculators.siesta.import_functions import (get_valence_charge,
27 read_rho,
28 read_vca_synth_block)
29from ase.calculators.siesta.parameters import (PAOBasisBlock, Species,
30 format_fdf)
31from ase.data import atomic_numbers
32from ase.io.siesta import read_siesta_xv
33from ase.units import Bohr, Ry, eV
34from ase.utils import deprecated
36meV = 0.001 * eV
39def parse_siesta_version(output: bytes) -> str:
40 match = re.search(rb'Siesta Version\s*:\s*(\S+)', output)
42 if match is None:
43 raise RuntimeError('Could not get Siesta version info from output '
44 '{!r}'.format(output))
46 string = match.group(1).decode('ascii')
47 return string
50def get_siesta_version(executable: str) -> str:
51 """ Return SIESTA version number.
53 Run the command, for instance 'siesta' and
54 then parse the output in order find the
55 version number.
56 """
57 # XXX We need a test of this kind of function. But Siesta().command
58 # is not enough to tell us how to run Siesta, because it could contain
59 # all sorts of mpirun and other weird parts.
61 temp_dirname = tempfile.mkdtemp(prefix='siesta-version-check-')
62 try:
63 from subprocess import PIPE, Popen
64 proc = Popen([executable],
65 stdin=PIPE,
66 stdout=PIPE,
67 stderr=PIPE,
68 cwd=temp_dirname)
69 output, _ = proc.communicate()
70 # We are not providing any input, so Siesta will give us a failure
71 # saying that it has no Chemical_species_label and exit status 1
72 # (as of siesta-4.1-b4)
73 finally:
74 shutil.rmtree(temp_dirname)
76 return parse_siesta_version(output)
79def bandpath2bandpoints(path):
80 lines = []
81 add = lines.append
83 add('BandLinesScale ReciprocalLatticeVectors\n')
84 add('%block BandPoints\n')
85 for kpt in path.kpts:
86 add(' {:18.15f} {:18.15f} {:18.15f}\n'.format(*kpt))
87 add('%endblock BandPoints')
88 return ''.join(lines)
91def read_bands_file(fd):
92 efermi = float(next(fd))
93 next(fd) # Appears to be max/min energy. Not important for us
94 header = next(fd) # Array shape: nbands, nspins, nkpoints
95 nbands, nspins, nkpts = np.array(header.split()).astype(int)
97 # three fields for kpt coords, then all the energies
98 ntokens = nbands * nspins + 3
100 # Read energies for each kpoint:
101 data = []
102 for i in range(nkpts):
103 line = next(fd)
104 tokens = line.split()
105 while len(tokens) < ntokens:
106 # Multirow table. Keep adding lines until the table ends,
107 # which should happen exactly when we have all the energies
108 # for this kpoint.
109 line = next(fd)
110 tokens += line.split()
111 assert len(tokens) == ntokens
112 values = np.array(tokens).astype(float)
113 data.append(values)
115 data = np.array(data)
116 assert len(data) == nkpts
117 kpts = data[:, :3]
118 energies = data[:, 3:]
119 energies = energies.reshape(nkpts, nspins, nbands)
120 assert energies.shape == (nkpts, nspins, nbands)
121 return kpts, energies, efermi
124def resolve_band_structure(path, kpts, energies, efermi):
125 """Convert input BandPath along with Siesta outputs into BS object."""
126 # Right now this function doesn't do much.
127 #
128 # Not sure how the output kpoints in the siesta.bands file are derived.
129 # They appear to be related to the lattice parameter.
130 #
131 # We should verify that they are consistent with our input path,
132 # but since their meaning is unclear, we can't quite do so.
133 #
134 # Also we should perhaps verify the cell. If we had the cell, we
135 # could construct the bandpath from scratch (i.e., pure outputs).
136 from ase.spectrum.band_structure import BandStructure
137 ksn2e = energies
138 skn2e = np.swapaxes(ksn2e, 0, 1)
139 bs = BandStructure(path, skn2e, reference=efermi)
140 return bs
143def is_along_cartesian(norm_dir: np.ndarray) -> bool:
144 """Return whether `norm_dir` is along a Cartesian coordidate."""
145 directions = [
146 [+1, 0, 0],
147 [-1, 0, 0],
148 [0, +1, 0],
149 [0, -1, 0],
150 [0, 0, +1],
151 [0, 0, -1],
152 ]
153 for direction in directions:
154 if np.allclose(norm_dir, direction, rtol=0.0, atol=1e-6):
155 return True
156 return False
159class SiestaParameters(Parameters):
160 """Parameters class for the calculator.
161 Documented in BaseSiesta.__init__
163 """
165 def __init__(
166 self,
167 label='siesta',
168 mesh_cutoff=200 * Ry,
169 energy_shift=100 * meV,
170 kpts=None,
171 xc='LDA',
172 basis_set='DZP',
173 spin='non-polarized',
174 species=(),
175 pseudo_qualifier=None,
176 pseudo_path=None,
177 symlink_pseudos=None,
178 atoms=None,
179 restart=None,
180 fdf_arguments=None,
181 atomic_coord_format='xyz',
182 bandpath=None):
183 kwargs = locals()
184 kwargs.pop('self')
185 Parameters.__init__(self, **kwargs)
188_DEPRECATED_SIESTA_COMMAND_MESSAGE = (
189 'Please use ``$ASE_SIESTA_COMMAND`` and not '
190 '``$SIESTA_COMMAND``, which will be ignored '
191 'in the future. The new command format will not '
192 'work with the "<%s > %s" specification. Use '
193 'instead e.g. "ASE_SIESTA_COMMAND=siesta'
194 ' < PREFIX.fdf > PREFIX.out", where PREFIX will '
195 'automatically be replaced by calculator label.'
196)
199def _nonpolarized_alias(_: List, kwargs: Dict[str, Any]) -> bool:
200 if kwargs.get("spin", None) == "UNPOLARIZED":
201 kwargs["spin"] = "non-polarized"
202 return True
203 return False
206class Siesta(FileIOCalculator):
207 """Calculator interface to the SIESTA code.
208 """
209 # Siesta manual does not document many of the basis names.
210 # basis_specs.f has a ton of aliases for each.
211 # Let's just list one of each type then.
212 #
213 # Maybe we should be less picky about these keyword names.
214 allowed_basis_names = ['SZ', 'SZP',
215 'DZ', 'DZP', 'DZP2',
216 'TZ', 'TZP', 'TZP2', 'TZP3']
217 allowed_spins = ['non-polarized', 'collinear',
218 'non-collinear', 'spin-orbit']
219 allowed_xc = {
220 'LDA': ['PZ', 'CA', 'PW92'],
221 'GGA': ['PW91', 'PBE', 'revPBE', 'RPBE',
222 'WC', 'AM05', 'PBEsol', 'PBEJsJrLO',
223 'PBEGcGxLO', 'PBEGcGxHEG', 'BLYP'],
224 'VDW': ['DRSLL', 'LMKLL', 'KBM', 'C09', 'BH', 'VV']}
226 name = 'siesta'
227 _legacy_default_command = 'siesta < PREFIX.fdf > PREFIX.out'
228 implemented_properties = [
229 'energy',
230 'free_energy',
231 'forces',
232 'stress',
233 'dipole',
234 'eigenvalues',
235 'density',
236 'fermi_energy']
238 # Dictionary of valid input vaiables.
239 default_parameters = SiestaParameters()
241 # XXX Not a ASE standard mechanism (yet). We need to communicate to
242 # ase.spectrum.band_structure.calculate_band_structure() that we expect
243 # it to use the bandpath keyword.
244 accepts_bandpath_keyword = True
246 def __init__(self, command=None, **kwargs):
247 f"""ASE interface to the SIESTA code.
249 Parameters:
250 - label : The basename of all files created during
251 calculation.
252 - mesh_cutoff : Energy in eV.
253 The mesh cutoff energy for determining number of
254 grid points in the matrix-element calculation.
255 - energy_shift : Energy in eV
256 The confining energy of the basis set generation.
257 - kpts : Tuple of 3 integers, the k-points in different
258 directions.
259 - xc : The exchange-correlation potential. Can be set to
260 any allowed value for either the Siesta
261 XC.funtional or XC.authors keyword. Default "LDA"
262 - basis_set : "SZ"|"SZP"|"DZ"|"DZP"|"TZP", strings which specify
263 the type of functions basis set.
264 - spin : "non-polarized"|"collinear"|
265 "non-collinear|spin-orbit".
266 The level of spin description to be used.
267 - species : None|list of Species objects. The species objects
268 can be used to to specify the basis set,
269 pseudopotential and whether the species is ghost.
270 The tag on the atoms object and the element is used
271 together to identify the species.
272 - pseudo_path : None|path. This path is where
273 pseudopotentials are taken from.
274 If None is given, then then the path given
275 in $SIESTA_PP_PATH will be used.
276 - pseudo_qualifier: None|string. This string will be added to the
277 pseudopotential path that will be retrieved.
278 For hydrogen with qualifier "abc" the
279 pseudopotential "H.abc.psf" will be retrieved.
280 - symlink_pseudos: None|bool
281 If true, symlink pseudopotentials
282 into the calculation directory, else copy them.
283 Defaults to true on Unix and false on Windows.
284 - atoms : The Atoms object.
285 - restart : str. Prefix for restart file.
286 May contain a directory.
287 Default is None, don't restart.
288 - fdf_arguments: Explicitly given fdf arguments. Dictonary using
289 Siesta keywords as given in the manual. List values
290 are written as fdf blocks with each element on a
291 separate line, while tuples will write each element
292 in a single line. ASE units are assumed in the
293 input.
294 - atomic_coord_format: "xyz"|"zmatrix", strings to switch between
295 the default way of entering the system's geometry
296 (via the block AtomicCoordinatesAndAtomicSpecies)
297 and a recent method via the block Zmatrix. The
298 block Zmatrix allows to specify basic geometry
299 constrains such as realized through the ASE classes
300 FixAtom, FixedLine and FixedPlane.
301 .. deprecated:: 3.18.2
302 {_DEPRECATED_SIESTA_COMMAND_MESSAGE}
303 """
305 # Put in the default arguments.
306 parameters = self.default_parameters.__class__(**kwargs)
308 # Call the base class.
309 FileIOCalculator.__init__(
310 self,
311 command=command,
312 **parameters)
314 commandvar = self.cfg.get("SIESTA_COMMAND")
315 runfile = self.prefix + ".fdf"
316 outfile = self.prefix + ".out"
317 if commandvar is not None:
318 warnings.warn(_DEPRECATED_SIESTA_COMMAND_MESSAGE)
319 try:
320 self.command = commandvar % (runfile, outfile)
321 except TypeError as err:
322 msg = (
323 "The 'SIESTA_COMMAND' environment must be a format string "
324 "with two string arguments.\n"
325 "Example : 'siesta < %s > %s'.\n"
326 f"Got '{commandvar}'"
327 )
328 raise ValueError(msg) from err
330 def __getitem__(self, key):
331 """Convenience method to retrieve a parameter as
332 calculator[key] rather than calculator.parameters[key]
334 Parameters:
335 -key : str, the name of the parameters to get.
336 """
337 return self.parameters[key]
339 def species(self, atoms):
340 """Find all relevant species depending on the atoms object and
341 species input.
343 Parameters :
344 - atoms : An Atoms object.
345 """
346 # For each element use default species from the species input, or set
347 # up a default species from the general default parameters.
348 symbols = np.array(atoms.get_chemical_symbols())
349 tags = atoms.get_tags()
350 species = list(self['species'])
351 default_species = [
352 s for s in species
353 if (s['tag'] is None) and s['symbol'] in symbols]
354 default_symbols = [s['symbol'] for s in default_species]
355 for symbol in symbols:
356 if symbol not in default_symbols:
357 spec = Species(symbol=symbol,
358 basis_set=self['basis_set'],
359 tag=None)
360 default_species.append(spec)
361 default_symbols.append(symbol)
362 assert len(default_species) == len(np.unique(symbols))
364 # Set default species as the first species.
365 species_numbers = np.zeros(len(atoms), int)
366 i = 1
367 for spec in default_species:
368 mask = symbols == spec['symbol']
369 species_numbers[mask] = i
370 i += 1
372 # Set up the non-default species.
373 non_default_species = [s for s in species if s['tag'] is not None]
374 for spec in non_default_species:
375 mask1 = (tags == spec['tag'])
376 mask2 = (symbols == spec['symbol'])
377 mask = np.logical_and(mask1, mask2)
378 if sum(mask) > 0:
379 species_numbers[mask] = i
380 i += 1
381 all_species = default_species + non_default_species
383 return all_species, species_numbers
385 @deprecated(
386 "The keyword 'UNPOLARIZED' has been deprecated,"
387 "and replaced by 'non-polarized'",
388 category=FutureWarning,
389 callback=_nonpolarized_alias,
390 )
391 def set(self, **kwargs):
392 """Set all parameters.
394 Parameters:
395 -kwargs : Dictionary containing the keywords defined in
396 SiestaParameters.
398 .. deprecated:: 3.18.2
399 The keyword 'UNPOLARIZED' has been deprecated and replaced by
400 'non-polarized'
401 """
403 # XXX Inserted these next few lines because set() would otherwise
404 # discard all previously set keywords to their defaults! --askhl
405 current = self.parameters.copy()
406 current.update(kwargs)
407 kwargs = current
409 # Find not allowed keys.
410 default_keys = list(self.__class__.default_parameters)
411 offending_keys = set(kwargs) - set(default_keys)
412 if len(offending_keys) > 0:
413 mess = "'set' does not take the keywords: %s "
414 raise ValueError(mess % list(offending_keys))
416 # Use the default parameters.
417 parameters = self.__class__.default_parameters.copy()
418 parameters.update(kwargs)
419 kwargs = parameters
421 # Check energy inputs.
422 for arg in ['mesh_cutoff', 'energy_shift']:
423 value = kwargs.get(arg)
424 if value is None:
425 continue
426 if not (isinstance(value, (float, int)) and value > 0):
427 mess = "'{}' must be a positive number(in eV), \
428 got '{}'".format(arg, value)
429 raise ValueError(mess)
431 # Check the basis set input.
432 if 'basis_set' in kwargs:
433 basis_set = kwargs['basis_set']
434 allowed = self.allowed_basis_names
435 if not (isinstance(basis_set, PAOBasisBlock) or
436 basis_set in allowed):
437 mess = f"Basis must be either {allowed}, got {basis_set}"
438 raise ValueError(mess)
440 # Check the spin input.
441 if 'spin' in kwargs:
442 spin = kwargs['spin']
443 if spin is not None and (spin.lower() not in self.allowed_spins):
444 mess = f"Spin must be {self.allowed_spins}, got '{spin}'"
445 raise ValueError(mess)
447 # Check the functional input.
448 xc = kwargs.get('xc', 'LDA')
449 if isinstance(xc, (tuple, list)) and len(xc) == 2:
450 functional, authors = xc
451 if functional.lower() not in [k.lower() for k in self.allowed_xc]:
452 mess = f"Unrecognized functional keyword: '{functional}'"
453 raise ValueError(mess)
455 lsauthorslower = [a.lower() for a in self.allowed_xc[functional]]
456 if authors.lower() not in lsauthorslower:
457 mess = "Unrecognized authors keyword for %s: '%s'"
458 raise ValueError(mess % (functional, authors))
460 elif xc in self.allowed_xc:
461 functional = xc
462 authors = self.allowed_xc[xc][0]
463 else:
464 found = False
465 for key, value in self.allowed_xc.items():
466 if xc in value:
467 found = True
468 functional = key
469 authors = xc
470 break
472 if not found:
473 raise ValueError(f"Unrecognized 'xc' keyword: '{xc}'")
474 kwargs['xc'] = (functional, authors)
476 # Check fdf_arguments.
477 if kwargs['fdf_arguments'] is None:
478 kwargs['fdf_arguments'] = {}
480 if not isinstance(kwargs['fdf_arguments'], dict):
481 raise TypeError("fdf_arguments must be a dictionary.")
483 # Call baseclass.
484 FileIOCalculator.set(self, **kwargs)
486 def set_fdf_arguments(self, fdf_arguments):
487 """ Set the fdf_arguments after the initialization of the
488 calculator.
489 """
490 self.validate_fdf_arguments(fdf_arguments)
491 FileIOCalculator.set(self, fdf_arguments=fdf_arguments)
493 def validate_fdf_arguments(self, fdf_arguments):
494 """ Raises error if the fdf_argument input is not a
495 dictionary of allowed keys.
496 """
497 # None is valid
498 if fdf_arguments is None:
499 return
501 # Type checking.
502 if not isinstance(fdf_arguments, dict):
503 raise TypeError("fdf_arguments must be a dictionary.")
505 def calculate(self,
506 atoms=None,
507 properties=['energy'],
508 system_changes=all_changes):
509 """Capture the RuntimeError from FileIOCalculator.calculate
510 and add a little debug information from the Siesta output.
512 See base FileIocalculator for documentation.
513 """
515 FileIOCalculator.calculate(
516 self,
517 atoms=atoms,
518 properties=properties,
519 system_changes=system_changes)
521 # The below snippet would run if calculate() failed but I have
522 # disabled it for now since it looks to be just for debugging.
523 # --askhl
524 """
525 # Here a test to check if the potential are in the right place!!!
526 except RuntimeError as e:
527 try:
528 fname = os.path.join(self.directory, self.label+'.out')
529 with open(fname, 'r') as fd:
530 lines = fd.readlines()
531 debug_lines = 10
532 print('##### %d last lines of the Siesta output' % debug_lines)
533 for line in lines[-20:]:
534 print(line.strip())
535 print('##### end of siesta output')
536 raise e
537 except:
538 raise e
539 """
541 def write_input(self, atoms, properties=None, system_changes=None):
542 """Write input (fdf)-file.
543 See calculator.py for further details.
545 Parameters:
546 - atoms : The Atoms object to write.
547 - properties : The properties which should be calculated.
548 - system_changes : List of properties changed since last run.
549 """
550 # Call base calculator.
551 FileIOCalculator.write_input(
552 self,
553 atoms=atoms,
554 properties=properties,
555 system_changes=system_changes)
557 if system_changes is None and properties is None:
558 return
560 filename = self.getpath(ext='fdf')
562 # On any changes, remove all analysis files.
563 if system_changes is not None:
564 self.remove_analysis()
566 # Start writing the file.
567 with open(filename, 'w') as fd:
568 # Write system name and label.
569 fd.write(format_fdf('SystemName', self.prefix))
570 fd.write(format_fdf('SystemLabel', self.prefix))
571 fd.write("\n")
573 # Write explicitly given options first to
574 # allow the user to override anything.
575 fdf_arguments = self['fdf_arguments']
576 keys = sorted(fdf_arguments.keys())
577 for key in keys:
578 fd.write(format_fdf(key, fdf_arguments[key]))
580 # Force siesta to return error on no convergence.
581 # as default consistent with ASE expectations.
582 if 'SCFMustConverge' not in fdf_arguments.keys():
583 fd.write(format_fdf('SCFMustConverge', True))
584 fd.write("\n")
586 # Write spin level.
587 fd.write(format_fdf('Spin ', self['spin']))
588 # Spin backwards compatibility.
589 if self['spin'] == 'collinear':
590 fd.write(
591 format_fdf(
592 'SpinPolarized',
593 (True,
594 "# Backwards compatibility.")))
595 elif self['spin'] == 'non-collinear':
596 fd.write(
597 format_fdf(
598 'NonCollinearSpin',
599 (True,
600 "# Backwards compatibility.")))
602 # Write functional.
603 functional, authors = self['xc']
604 fd.write(format_fdf('XC.functional', functional))
605 fd.write(format_fdf('XC.authors', authors))
606 fd.write("\n")
608 # Write mesh cutoff and energy shift.
609 fd.write(format_fdf('MeshCutoff',
610 (self['mesh_cutoff'], 'eV')))
611 fd.write(format_fdf('PAO.EnergyShift',
612 (self['energy_shift'], 'eV')))
613 fd.write("\n")
615 # Write the minimal arg
616 self._write_species(fd, atoms)
617 self._write_structure(fd, atoms)
619 # Use the saved density matrix if only 'cell' and 'positions'
620 # have changed.
621 if (system_changes is None or
622 ('numbers' not in system_changes and
623 'initial_magmoms' not in system_changes and
624 'initial_charges' not in system_changes)):
625 fd.write(format_fdf('DM.UseSaveDM', True))
627 # Save density.
628 if 'density' in properties:
629 fd.write(format_fdf('SaveRho', True))
631 self._write_kpts(fd)
633 if self['bandpath'] is not None:
634 lines = bandpath2bandpoints(self['bandpath'])
635 fd.write(lines)
636 fd.write('\n')
638 def read(self, filename):
639 """Read structural parameters from file .XV file
640 Read other results from other files
641 filename : siesta.XV
642 """
644 fname = self.getpath(filename)
645 if not os.path.exists(fname):
646 raise ReadError(f"The restart file '{fname}' does not exist")
647 with open(fname) as fd:
648 self.atoms = read_siesta_xv(fd)
649 self.read_results()
651 def getpath(self, fname=None, ext=None):
652 """ Returns the directory/fname string """
653 if fname is None:
654 fname = self.prefix
655 if ext is not None:
656 fname = f'{fname}.{ext}'
657 return os.path.join(self.directory, fname)
659 def remove_analysis(self):
660 """ Remove all analysis files"""
661 filename = self.getpath(ext='RHO')
662 if os.path.exists(filename):
663 os.remove(filename)
665 def _write_structure(self, fd, atoms):
666 """Translate the Atoms object to fdf-format.
668 Parameters
669 ----------
670 fd : IO
671 An open file object.
672 atoms: Atoms
673 An atoms object.
674 """
675 cell = atoms.cell
676 fd.write('\n')
678 if cell.rank in [1, 2]:
679 raise ValueError('Expected 3D unit cell or no unit cell. You may '
680 'wish to add vacuum along some directions.')
682 # Write lattice vectors
683 if np.any(cell):
684 fd.write(format_fdf('LatticeConstant', '1.0 Ang'))
685 fd.write('%block LatticeVectors\n')
686 for i in range(3):
687 for j in range(3):
688 s = (' %.15f' % cell[i, j]).rjust(16) + ' '
689 fd.write(s)
690 fd.write('\n')
691 fd.write('%endblock LatticeVectors\n')
692 fd.write('\n')
694 self._write_atomic_coordinates(fd, atoms)
696 # Write magnetic moments.
697 magmoms = atoms.get_initial_magnetic_moments()
699 # The DM.InitSpin block must be written to initialize to
700 # no spin. SIESTA default is FM initialization, if the
701 # block is not written, but we must conform to the
702 # atoms object.
703 if magmoms is not None:
704 if len(magmoms) == 0:
705 fd.write('#Empty block forces ASE initialization.\n')
707 fd.write('%block DM.InitSpin\n')
708 if len(magmoms) != 0 and isinstance(magmoms[0], np.ndarray):
709 for n, M in enumerate(magmoms):
710 if M[0] != 0:
711 fd.write(
712 ' %d %.14f %.14f %.14f \n' %
713 (n + 1, M[0], M[1], M[2]))
714 elif len(magmoms) != 0 and isinstance(magmoms[0], float):
715 for n, M in enumerate(magmoms):
716 if M != 0:
717 fd.write(' %d %.14f \n' % (n + 1, M))
718 fd.write('%endblock DM.InitSpin\n')
719 fd.write('\n')
721 def _write_atomic_coordinates(self, fd, atoms: Atoms):
722 """Write atomic coordinates.
724 Parameters
725 ----------
726 fd : IO
727 An open file object.
728 atoms : Atoms
729 An atoms object.
730 """
731 af = self.parameters["atomic_coord_format"].lower()
732 if af == 'xyz':
733 self._write_atomic_coordinates_xyz(fd, atoms)
734 elif af == 'zmatrix':
735 self._write_atomic_coordinates_zmatrix(fd, atoms)
736 else:
737 raise RuntimeError(f'Unknown atomic_coord_format: {af}')
739 def _write_atomic_coordinates_xyz(self, fd, atoms: Atoms):
740 """Write atomic coordinates.
742 Parameters
743 ----------
744 fd : IO
745 An open file object.
746 atoms : Atoms
747 An atoms object.
748 """
749 species, species_numbers = self.species(atoms)
750 fd.write('\n')
751 fd.write('AtomicCoordinatesFormat Ang\n')
752 fd.write('%block AtomicCoordinatesAndAtomicSpecies\n')
753 for atom, number in zip(atoms, species_numbers):
754 xyz = atom.position
755 line = (' %.9f' % xyz[0]).rjust(16) + ' '
756 line += (' %.9f' % xyz[1]).rjust(16) + ' '
757 line += (' %.9f' % xyz[2]).rjust(16) + ' '
758 line += str(number) + '\n'
759 fd.write(line)
760 fd.write('%endblock AtomicCoordinatesAndAtomicSpecies\n')
761 fd.write('\n')
763 origin = tuple(-atoms.get_celldisp().flatten())
764 if any(origin):
765 fd.write('%block AtomicCoordinatesOrigin\n')
766 fd.write(' %.4f %.4f %.4f\n' % origin)
767 fd.write('%endblock AtomicCoordinatesOrigin\n')
768 fd.write('\n')
770 def _write_atomic_coordinates_zmatrix(self, fd, atoms: Atoms):
771 """Write atomic coordinates in Z-matrix format.
773 Parameters
774 ----------
775 fd : IO
776 An open file object.
777 atoms : Atoms
778 An atoms object.
779 """
780 species, species_numbers = self.species(atoms)
781 fd.write('\n')
782 fd.write('ZM.UnitsLength Ang\n')
783 fd.write('%block Zmatrix\n')
784 fd.write(' cartesian\n')
785 fstr = "{:5d}" + "{:20.10f}" * 3 + "{:3d}" * 3 + "{:7d} {:s}\n"
786 a2constr = self.make_xyz_constraints(atoms)
787 a2p, a2s = atoms.get_positions(), atoms.get_chemical_symbols()
788 for ia, (sp, xyz, ccc, sym) in enumerate(zip(species_numbers,
789 a2p,
790 a2constr,
791 a2s)):
792 fd.write(fstr.format(
793 sp, xyz[0], xyz[1], xyz[2], ccc[0],
794 ccc[1], ccc[2], ia + 1, sym))
795 fd.write('%endblock Zmatrix\n')
797 origin = tuple(-atoms.get_celldisp().flatten())
798 if any(origin):
799 fd.write('%block AtomicCoordinatesOrigin\n')
800 fd.write(' %.4f %.4f %.4f\n' % origin)
801 fd.write('%endblock AtomicCoordinatesOrigin\n')
802 fd.write('\n')
804 def make_xyz_constraints(self, atoms: Atoms):
805 """ Create coordinate-resolved list of constraints [natoms, 0:3]
806 The elements of the list must be integers 0 or 1
807 1 -- means that the coordinate will be updated during relaxation
808 0 -- mains that the coordinate will be fixed during relaxation
809 """
810 import sys
812 from ase.constraints import (FixAtoms, FixCartesian, FixedLine,
813 FixedPlane)
815 a2c = np.ones((len(atoms), 3), dtype=int) # (0: fixed, 1: updated)
816 for c in atoms.constraints:
817 if isinstance(c, FixAtoms):
818 a2c[c.get_indices()] = 0
819 elif isinstance(c, FixedLine):
820 norm_dir = c.dir / np.linalg.norm(c.dir)
821 if not is_along_cartesian(norm_dir):
822 raise RuntimeError(
823 'norm_dir: {} -- must be one of the Cartesian axes...'
824 .format(norm_dir))
825 a2c[c.get_indices()] = norm_dir.round().astype(int)
826 elif isinstance(c, FixedPlane):
827 norm_dir = c.dir / np.linalg.norm(c.dir)
828 if not is_along_cartesian(norm_dir):
829 raise RuntimeError(
830 'norm_dir: {} -- must be one of the Cartesian axes...'
831 .format(norm_dir))
832 a2c[c.get_indices()] = abs(1 - norm_dir.round().astype(int))
833 elif isinstance(c, FixCartesian):
834 a2c[c.get_indices()] = c.mask.astype(int)
835 else:
836 warnings.warn('Constraint {} is ignored at {}'
837 .format(str(c), sys._getframe().f_code))
838 return a2c
840 def _write_kpts(self, fd):
841 """Write kpts.
843 Parameters:
844 - f : Open filename.
845 """
846 if self["kpts"] is None:
847 return
848 kpts = np.array(self['kpts'])
849 fd.write('\n')
850 fd.write('#KPoint grid\n')
851 fd.write('%block kgrid_Monkhorst_Pack\n')
853 for i in range(3):
854 s = ''
855 if i < len(kpts):
856 number = kpts[i]
857 displace = 0.0
858 else:
859 number = 1
860 displace = 0
861 for j in range(3):
862 if j == i:
863 write_this = number
864 else:
865 write_this = 0
866 s += ' %d ' % write_this
867 s += '%1.1f\n' % displace
868 fd.write(s)
869 fd.write('%endblock kgrid_Monkhorst_Pack\n')
870 fd.write('\n')
872 def _write_species(self, fd, atoms):
873 """Write input related the different species.
875 Parameters:
876 - f: An open file object.
877 - atoms: An atoms object.
878 """
879 species, species_numbers = self.species(atoms)
881 if self['pseudo_path'] is not None:
882 pseudo_path = self['pseudo_path']
883 elif 'SIESTA_PP_PATH' in self.cfg:
884 pseudo_path = self.cfg['SIESTA_PP_PATH']
885 else:
886 mess = "Please set the environment variable 'SIESTA_PP_PATH'"
887 raise Exception(mess)
889 fd.write(format_fdf('NumberOfSpecies', len(species)))
890 fd.write(format_fdf('NumberOfAtoms', len(atoms)))
892 pao_basis = []
893 chemical_labels = []
894 basis_sizes = []
895 synth_blocks = []
896 for species_number, spec in enumerate(species):
897 species_number += 1
898 symbol = spec['symbol']
899 atomic_number = atomic_numbers[symbol]
901 if spec['pseudopotential'] is None:
902 if self.pseudo_qualifier() == '':
903 label = symbol
904 pseudopotential = label + '.psf'
905 else:
906 label = '.'.join([symbol, self.pseudo_qualifier()])
907 pseudopotential = label + '.psf'
908 else:
909 pseudopotential = spec['pseudopotential']
910 label = os.path.basename(pseudopotential)
911 label = '.'.join(label.split('.')[:-1])
913 if not os.path.isabs(pseudopotential):
914 pseudopotential = join(pseudo_path, pseudopotential)
916 if not os.path.exists(pseudopotential):
917 mess = f"Pseudopotential '{pseudopotential}' not found"
918 raise RuntimeError(mess)
920 name = os.path.basename(pseudopotential)
921 name = name.split('.')
922 name.insert(-1, str(species_number))
923 if spec['ghost']:
924 name.insert(-1, 'ghost')
925 atomic_number = -atomic_number
927 name = '.'.join(name)
928 pseudo_targetpath = self.getpath(name)
930 if join(os.getcwd(), name) != pseudopotential:
931 if islink(pseudo_targetpath) or isfile(pseudo_targetpath):
932 os.remove(pseudo_targetpath)
933 symlink_pseudos = self['symlink_pseudos']
935 if symlink_pseudos is None:
936 symlink_pseudos = not os.name == 'nt'
938 if symlink_pseudos:
939 os.symlink(pseudopotential, pseudo_targetpath)
940 else:
941 shutil.copy(pseudopotential, pseudo_targetpath)
943 if spec['excess_charge'] is not None:
944 atomic_number += 200
945 n_atoms = sum(np.array(species_numbers) == species_number)
947 paec = float(spec['excess_charge']) / n_atoms
948 vc = get_valence_charge(pseudopotential)
949 fraction = float(vc + paec) / vc
950 pseudo_head = name[:-4]
951 fractional_command = self.cfg['SIESTA_UTIL_FRACTIONAL']
952 cmd = '{} {} {:.7f}'.format(fractional_command,
953 pseudo_head,
954 fraction)
955 os.system(cmd)
957 pseudo_head += '-Fraction-%.5f' % fraction
958 synth_pseudo = pseudo_head + '.psf'
959 synth_block_filename = pseudo_head + '.synth'
960 os.remove(name)
961 shutil.copyfile(synth_pseudo, name)
962 synth_block = read_vca_synth_block(
963 synth_block_filename,
964 species_number=species_number)
965 synth_blocks.append(synth_block)
967 if len(synth_blocks) > 0:
968 fd.write(format_fdf('SyntheticAtoms', list(synth_blocks)))
970 label = '.'.join(np.array(name.split('.'))[:-1])
971 string = ' %d %d %s' % (species_number, atomic_number, label)
972 chemical_labels.append(string)
973 if isinstance(spec['basis_set'], PAOBasisBlock):
974 pao_basis.append(spec['basis_set'].script(label))
975 else:
976 basis_sizes.append((" " + label, spec['basis_set']))
977 fd.write(format_fdf('ChemicalSpecieslabel', chemical_labels))
978 fd.write('\n')
979 fd.write(format_fdf('PAO.Basis', pao_basis))
980 fd.write(format_fdf('PAO.BasisSizes', basis_sizes))
981 fd.write('\n')
983 def pseudo_qualifier(self):
984 """Get the extra string used in the middle of the pseudopotential.
985 The retrieved pseudopotential for a specific element will be
986 'H.xxx.psf' for the element 'H' with qualifier 'xxx'. If qualifier
987 is set to None then the qualifier is set to functional name.
988 """
989 if self['pseudo_qualifier'] is None:
990 return self['xc'][0].lower()
991 else:
992 return self['pseudo_qualifier']
994 def read_results(self):
995 """Read the results.
996 """
997 self.read_number_of_grid_points()
998 self.read_energy()
999 self.read_forces_stress()
1000 self.read_eigenvalues()
1001 self.read_kpoints()
1002 self.read_dipole()
1003 self.read_pseudo_density()
1004 self.read_hsx()
1005 self.read_dim()
1006 if self.results['hsx'] is not None:
1007 self.read_pld(self.results['hsx'].norbitals,
1008 len(self.atoms))
1009 self.atoms.cell = self.results['pld'].cell * Bohr
1010 else:
1011 self.results['pld'] = None
1013 self.read_wfsx()
1014 self.read_ion(self.atoms)
1016 self.read_bands()
1018 def read_bands(self):
1019 bandpath = self['bandpath']
1020 if bandpath is None:
1021 return
1023 if len(bandpath.kpts) < 1:
1024 return
1026 fname = self.getpath(ext='bands')
1027 with open(fname) as fd:
1028 kpts, energies, efermi = read_bands_file(fd)
1029 bs = resolve_band_structure(bandpath, kpts, energies, efermi)
1030 self.results['bandstructure'] = bs
1032 def band_structure(self):
1033 return self.results['bandstructure']
1035 def read_ion(self, atoms):
1036 """
1037 Read the ion.xml file of each specie
1038 """
1039 from ase.calculators.siesta.import_ion_xml import get_ion
1041 species, species_numbers = self.species(atoms)
1043 self.results['ion'] = {}
1044 for species_number, spec in enumerate(species):
1045 species_number += 1
1047 symbol = spec['symbol']
1048 atomic_number = atomic_numbers[symbol]
1050 if spec['pseudopotential'] is None:
1051 if self.pseudo_qualifier() == '':
1052 label = symbol
1053 else:
1054 label = '.'.join([symbol, self.pseudo_qualifier()])
1055 pseudopotential = self.getpath(label, 'psf')
1056 else:
1057 pseudopotential = spec['pseudopotential']
1058 label = os.path.basename(pseudopotential)
1059 label = '.'.join(label.split('.')[:-1])
1061 name = os.path.basename(pseudopotential)
1062 name = name.split('.')
1063 name.insert(-1, str(species_number))
1064 if spec['ghost']:
1065 name.insert(-1, 'ghost')
1066 atomic_number = -atomic_number
1067 name = '.'.join(name)
1069 label = '.'.join(np.array(name.split('.'))[:-1])
1071 if label not in self.results['ion']:
1072 fname = self.getpath(label, 'ion.xml')
1073 if os.path.isfile(fname):
1074 self.results['ion'][label] = get_ion(fname)
1076 def read_hsx(self):
1077 """
1078 Read the siesta HSX file.
1079 return a namedtuple with the following arguments:
1080 'norbitals', 'norbitals_sc', 'nspin', 'nonzero',
1081 'is_gamma', 'sc_orb2uc_orb', 'row2nnzero', 'sparse_ind2column',
1082 'H_sparse', 'S_sparse', 'aB2RaB_sparse', 'total_elec_charge', 'temp'
1083 """
1084 from ase.calculators.siesta.import_functions import readHSX
1086 filename = self.getpath(ext='HSX')
1087 if isfile(filename):
1088 self.results['hsx'] = readHSX(filename)
1089 else:
1090 self.results['hsx'] = None
1092 def read_dim(self):
1093 """
1094 Read the siesta DIM file
1095 Retrun a namedtuple with the following arguments:
1096 'natoms_sc', 'norbitals_sc', 'norbitals', 'nspin',
1097 'nnonzero', 'natoms_interacting'
1098 """
1099 from ase.calculators.siesta.import_functions import readDIM
1101 filename = self.getpath(ext='DIM')
1102 if isfile(filename):
1103 self.results['dim'] = readDIM(filename)
1104 else:
1105 self.results['dim'] = None
1107 def read_pld(self, norb, natms):
1108 """
1109 Read the siesta PLD file
1110 Return a namedtuple with the following arguments:
1111 'max_rcut', 'orb2ao', 'orb2uorb', 'orb2occ', 'atm2sp',
1112 'atm2shift', 'coord_sc', 'cell', 'nunit_cells'
1113 """
1114 from ase.calculators.siesta.import_functions import readPLD
1116 filename = self.getpath(ext='PLD')
1117 if isfile(filename):
1118 self.results['pld'] = readPLD(filename, norb, natms)
1119 else:
1120 self.results['pld'] = None
1122 def read_wfsx(self):
1123 """
1124 Read the siesta WFSX file
1125 Return a namedtuple with the following arguments:
1126 """
1127 from ase.calculators.siesta.import_functions import readWFSX
1129 fname_woext = os.path.join(self.directory, self.prefix)
1131 if isfile(fname_woext + '.WFSX'):
1132 filename = fname_woext + '.WFSX'
1133 self.results['wfsx'] = readWFSX(filename)
1134 elif isfile(fname_woext + '.fullBZ.WFSX'):
1135 filename = fname_woext + '.fullBZ.WFSX'
1136 readWFSX(filename)
1137 self.results['wfsx'] = readWFSX(filename)
1138 else:
1139 self.results['wfsx'] = None
1141 def read_pseudo_density(self):
1142 """Read the density if it is there."""
1143 filename = self.getpath(ext='RHO')
1144 if isfile(filename):
1145 self.results['density'] = read_rho(filename)
1147 def read_number_of_grid_points(self):
1148 """Read number of grid points from SIESTA's text-output file. """
1150 fname = self.getpath(ext='out')
1151 with open(fname) as fd:
1152 for line in fd:
1153 line = line.strip().lower()
1154 if line.startswith('initmesh: mesh ='):
1155 n_points = [int(word) for word in line.split()[3:8:2]]
1156 self.results['n_grid_point'] = n_points
1157 break
1158 else:
1159 raise RuntimeError
1161 def read_energy(self):
1162 """Read energy from SIESTA's text-output file.
1163 """
1164 fname = self.getpath(ext='out')
1165 with open(fname) as fd:
1166 text = fd.read().lower()
1168 assert 'final energy' in text
1169 lines = iter(text.split('\n'))
1171 # Get the energy and free energy the last time it appears
1172 for line in lines:
1173 has_energy = line.startswith('siesta: etot =')
1174 if has_energy:
1175 self.results['energy'] = float(line.split()[-1])
1176 line = next(lines)
1177 self.results['free_energy'] = float(line.split()[-1])
1179 if ('energy' not in self.results or
1180 'free_energy' not in self.results):
1181 raise RuntimeError
1183 def read_forces_stress(self):
1184 """Read the forces and stress from the FORCE_STRESS file.
1185 """
1186 fname = self.getpath('FORCE_STRESS')
1187 with open(fname) as fd:
1188 lines = fd.readlines()
1190 stress_lines = lines[1:4]
1191 stress = np.empty((3, 3))
1192 for i in range(3):
1193 line = stress_lines[i].strip().split(' ')
1194 line = [s for s in line if len(s) > 0]
1195 stress[i] = [float(s) for s in line]
1197 self.results['stress'] = np.array(
1198 [stress[0, 0], stress[1, 1], stress[2, 2],
1199 stress[1, 2], stress[0, 2], stress[0, 1]])
1201 self.results['stress'] *= Ry / Bohr**3
1203 start = 5
1204 self.results['forces'] = np.zeros((len(lines) - start, 3), float)
1205 for i in range(start, len(lines)):
1206 line = [s for s in lines[i].strip().split(' ') if len(s) > 0]
1207 self.results['forces'][i - start] = [float(s) for s in line[2:5]]
1209 self.results['forces'] *= Ry / Bohr
1211 def read_eigenvalues(self):
1212 """ A robust procedure using the suggestion by Federico Marchesin """
1214 file_name = self.getpath(ext='EIG')
1215 try:
1216 with open(file_name) as fd:
1217 self.results['fermi_energy'] = float(fd.readline())
1218 n, num_hamilton_dim, nkp = map(int, fd.readline().split())
1219 _ee = np.split(
1220 np.array(fd.read().split()).astype(float), nkp)
1221 except OSError:
1222 return 1
1224 n_spin = 1 if num_hamilton_dim > 2 else num_hamilton_dim
1225 ksn2e = np.delete(_ee, 0, 1).reshape([nkp, n_spin, n])
1227 eig_array = np.empty((n_spin, nkp, n))
1228 eig_array[:] = np.inf
1230 for k, sn2e in enumerate(ksn2e):
1231 for s, n2e in enumerate(sn2e):
1232 eig_array[s, k, :] = n2e
1234 assert np.isfinite(eig_array).all()
1236 self.results['eigenvalues'] = eig_array
1237 return 0
1239 def read_kpoints(self):
1240 """ Reader of the .KP files """
1242 fname = self.getpath(ext='KP')
1243 try:
1244 with open(fname) as fd:
1245 nkp = int(next(fd))
1246 kpoints = np.empty((nkp, 3))
1247 kweights = np.empty(nkp)
1249 for i in range(nkp):
1250 line = next(fd)
1251 tokens = line.split()
1252 numbers = np.array(tokens[1:]).astype(float)
1253 kpoints[i] = numbers[:3]
1254 kweights[i] = numbers[3]
1256 except (OSError):
1257 return 1
1259 self.results['kpoints'] = kpoints
1260 self.results['kweights'] = kweights
1262 return 0
1264 def read_dipole(self):
1265 """Read dipole moment. """
1266 dipole = np.zeros([1, 3])
1267 with open(self.getpath(ext='out')) as fd:
1268 for line in fd:
1269 if line.rfind('Electric dipole (Debye)') > -1:
1270 dipole = np.array([float(f) for f in line.split()[5:8]])
1271 # debye to e*Ang
1272 self.results['dipole'] = dipole * 0.2081943482534
1274 def get_fermi_level(self):
1275 return self.results['fermi_energy']
1277 def get_k_point_weights(self):
1278 return self.results['kweights']
1280 def get_ibz_k_points(self):
1281 return self.results['kpoints']
1284class Siesta3_2(Siesta):
1285 @deprecated(
1286 "The Siesta3_2 calculator class will no longer be supported. "
1287 "Use 'ase.calculators.siesta.Siesta instead. "
1288 "If using the ASE interface with SIESTA 3.2 you must explicitly "
1289 "include the keywords 'SpinPolarized', 'NonCollinearSpin' and "
1290 "'SpinOrbit' if needed.",
1291 FutureWarning
1292 )
1293 def __init__(self, *args, **kwargs):
1294 """
1295 .. deprecated: 3.18.2
1296 The Siesta3_2 calculator class will no longer be supported.
1297 Use :class:`~ase.calculators.siesta.Siesta` instead.
1298 If using the ASE interface with SIESTA 3.2 you must explicitly
1299 include the keywords 'SpinPolarized', 'NonCollinearSpin' and
1300 'SpinOrbit' if needed.
1301 """
1302 Siesta.__init__(self, *args, **kwargs)