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