Coverage for /builds/debichem-team/python-ase/ase/calculators/siesta/siesta.py : 62.11%
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 """
145 def __init__(
146 self,
147 label='siesta',
148 mesh_cutoff=200 * Ry,
149 energy_shift=100 * meV,
150 kpts=None,
151 xc='LDA',
152 basis_set='DZP',
153 spin='non-polarized',
154 species=tuple(),
155 pseudo_qualifier=None,
156 pseudo_path=None,
157 symlink_pseudos=None,
158 atoms=None,
159 restart=None,
160 fdf_arguments=None,
161 atomic_coord_format='xyz',
162 bandpath=None):
163 kwargs = locals()
164 kwargs.pop('self')
165 Parameters.__init__(self, **kwargs)
168class Siesta(FileIOCalculator):
169 """Calculator interface to the SIESTA code.
170 """
171 # Siesta manual does not document many of the basis names.
172 # basis_specs.f has a ton of aliases for each.
173 # Let's just list one of each type then.
174 #
175 # Maybe we should be less picky about these keyword names.
176 allowed_basis_names = ['SZ', 'SZP',
177 'DZ', 'DZP', 'DZP2',
178 'TZ', 'TZP', 'TZP2', 'TZP3']
179 allowed_spins = ['non-polarized', 'collinear',
180 'non-collinear', 'spin-orbit']
181 allowed_xc = {
182 'LDA': ['PZ', 'CA', 'PW92'],
183 'GGA': ['PW91', 'PBE', 'revPBE', 'RPBE',
184 'WC', 'AM05', 'PBEsol', 'PBEJsJrLO',
185 'PBEGcGxLO', 'PBEGcGxHEG', 'BLYP'],
186 'VDW': ['DRSLL', 'LMKLL', 'KBM', 'C09', 'BH', 'VV']}
188 name = 'siesta'
189 command = 'siesta < PREFIX.fdf > PREFIX.out'
190 implemented_properties = [
191 'energy',
192 'forces',
193 'stress',
194 'dipole',
195 'eigenvalues',
196 'density',
197 'fermi_energy']
199 # Dictionary of valid input vaiables.
200 default_parameters = SiestaParameters()
202 # XXX Not a ASE standard mechanism (yet). We need to communicate to
203 # ase.spectrum.band_structure.calculate_band_structure() that we expect
204 # it to use the bandpath keyword.
205 accepts_bandpath_keyword = True
207 def __init__(self, command=None, **kwargs):
208 """ASE interface to the SIESTA code.
210 Parameters:
211 - label : The basename of all files created during
212 calculation.
213 - mesh_cutoff : Energy in eV.
214 The mesh cutoff energy for determining number of
215 grid points in the matrix-element calculation.
216 - energy_shift : Energy in eV
217 The confining energy of the basis set generation.
218 - kpts : Tuple of 3 integers, the k-points in different
219 directions.
220 - xc : The exchange-correlation potential. Can be set to
221 any allowed value for either the Siesta
222 XC.funtional or XC.authors keyword. Default "LDA"
223 - basis_set : "SZ"|"SZP"|"DZ"|"DZP"|"TZP", strings which specify
224 the type of functions basis set.
225 - spin : "non-polarized"|"collinear"|
226 "non-collinear|spin-orbit".
227 The level of spin description to be used.
228 - species : None|list of Species objects. The species objects
229 can be used to to specify the basis set,
230 pseudopotential and whether the species is ghost.
231 The tag on the atoms object and the element is used
232 together to identify the species.
233 - pseudo_path : None|path. This path is where
234 pseudopotentials are taken from.
235 If None is given, then then the path given
236 in $SIESTA_PP_PATH will be used.
237 - pseudo_qualifier: None|string. This string will be added to the
238 pseudopotential path that will be retrieved.
239 For hydrogen with qualifier "abc" the
240 pseudopotential "H.abc.psf" will be retrieved.
241 - symlink_pseudos: None|bool
242 If true, symlink pseudopotentials
243 into the calculation directory, else copy them.
244 Defaults to true on Unix and false on Windows.
245 - atoms : The Atoms object.
246 - restart : str. Prefix for restart file.
247 May contain a directory.
248 Default is None, don't restart.
249 - fdf_arguments: Explicitly given fdf arguments. Dictonary using
250 Siesta keywords as given in the manual. List values
251 are written as fdf blocks with each element on a
252 separate line, while tuples will write each element
253 in a single line. ASE units are assumed in the
254 input.
255 - atomic_coord_format: "xyz"|"zmatrix", strings to switch between
256 the default way of entering the system's geometry
257 (via the block AtomicCoordinatesAndAtomicSpecies)
258 and a recent method via the block Zmatrix. The
259 block Zmatrix allows to specify basic geometry
260 constrains such as realized through the ASE classes
261 FixAtom, FixedLine and FixedPlane.
262 """
264 # Put in the default arguments.
265 parameters = self.default_parameters.__class__(**kwargs)
267 # Call the base class.
268 FileIOCalculator.__init__(
269 self,
270 command=command,
271 **parameters)
273 # For compatibility with old variable name:
274 commandvar = os.environ.get('SIESTA_COMMAND')
275 if commandvar is not None:
276 warnings.warn('Please use $ASE_SIESTA_COMMAND and not '
277 '$SIESTA_COMMAND, which will be ignored '
278 'in the future. The new command format will not '
279 'work with the "<%s > %s" specification. Use '
280 'instead e.g. "ASE_SIESTA_COMMAND=siesta'
281 ' < PREFIX.fdf > PREFIX.out", where PREFIX will '
282 'automatically be replaced by calculator label',
283 np.VisibleDeprecationWarning)
284 runfile = self.prefix + '.fdf'
285 outfile = self.prefix + '.out'
286 try:
287 self.command = commandvar % (runfile, outfile)
288 except TypeError:
289 raise ValueError(
290 "The 'SIESTA_COMMAND' environment must " +
291 "be a format string" +
292 " with two string arguments.\n" +
293 "Example : 'siesta < %s > %s'.\n" +
294 "Got '%s'" % commandvar)
296 def __getitem__(self, key):
297 """Convenience method to retrieve a parameter as
298 calculator[key] rather than calculator.parameters[key]
300 Parameters:
301 -key : str, the name of the parameters to get.
302 """
303 return self.parameters[key]
305 def species(self, atoms):
306 """Find all relevant species depending on the atoms object and
307 species input.
309 Parameters :
310 - atoms : An Atoms object.
311 """
312 # For each element use default species from the species input, or set
313 # up a default species from the general default parameters.
314 symbols = np.array(atoms.get_chemical_symbols())
315 tags = atoms.get_tags()
316 species = list(self['species'])
317 default_species = [
318 s for s in species
319 if (s['tag'] is None) and s['symbol'] in symbols]
320 default_symbols = [s['symbol'] for s in default_species]
321 for symbol in symbols:
322 if symbol not in default_symbols:
323 spec = Species(symbol=symbol,
324 basis_set=self['basis_set'],
325 tag=None)
326 default_species.append(spec)
327 default_symbols.append(symbol)
328 assert len(default_species) == len(np.unique(symbols))
330 # Set default species as the first species.
331 species_numbers = np.zeros(len(atoms), int)
332 i = 1
333 for spec in default_species:
334 mask = symbols == spec['symbol']
335 species_numbers[mask] = i
336 i += 1
338 # Set up the non-default species.
339 non_default_species = [s for s in species if not s['tag'] is None]
340 for spec in non_default_species:
341 mask1 = (tags == spec['tag'])
342 mask2 = (symbols == spec['symbol'])
343 mask = np.logical_and(mask1, mask2)
344 if sum(mask) > 0:
345 species_numbers[mask] = i
346 i += 1
347 all_species = default_species + non_default_species
349 return all_species, species_numbers
351 def set(self, **kwargs):
352 """Set all parameters.
354 Parameters:
355 -kwargs : Dictionary containing the keywords defined in
356 SiestaParameters.
357 """
359 # XXX Inserted these next few lines because set() would otherwise
360 # discard all previously set keywords to their defaults! --askhl
361 current = self.parameters.copy()
362 current.update(kwargs)
363 kwargs = current
365 # Find not allowed keys.
366 default_keys = list(self.__class__.default_parameters)
367 offending_keys = set(kwargs) - set(default_keys)
368 if len(offending_keys) > 0:
369 mess = "'set' does not take the keywords: %s "
370 raise ValueError(mess % list(offending_keys))
372 # Use the default parameters.
373 parameters = self.__class__.default_parameters.copy()
374 parameters.update(kwargs)
375 kwargs = parameters
377 # Check energy inputs.
378 for arg in ['mesh_cutoff', 'energy_shift']:
379 value = kwargs.get(arg)
380 if value is None:
381 continue
382 if not (isinstance(value, (float, int)) and value > 0):
383 mess = "'%s' must be a positive number(in eV), \
384 got '%s'" % (arg, value)
385 raise ValueError(mess)
387 # Check the basis set input.
388 if 'basis_set' in kwargs:
389 basis_set = kwargs['basis_set']
390 allowed = self.allowed_basis_names
391 if not (isinstance(basis_set, PAOBasisBlock) or
392 basis_set in allowed):
393 mess = "Basis must be either %s, got %s" % (allowed, basis_set)
394 raise ValueError(mess)
396 # Check the spin input.
397 if 'spin' in kwargs:
398 if kwargs['spin'] == 'UNPOLARIZED':
399 warnings.warn("The keyword 'UNPOLARIZED' is deprecated,"
400 "and replaced by 'non-polarized'",
401 np.VisibleDeprecationWarning)
402 kwargs['spin'] = 'non-polarized'
404 spin = kwargs['spin']
405 if spin is not None and (spin.lower() not in self.allowed_spins):
406 mess = "Spin must be %s, got '%s'" % (self.allowed_spins, spin)
407 raise ValueError(mess)
409 # Check the functional input.
410 xc = kwargs.get('xc', 'LDA')
411 if isinstance(xc, (tuple, list)) and len(xc) == 2:
412 functional, authors = xc
413 if functional.lower() not in [k.lower() for k in self.allowed_xc]:
414 mess = "Unrecognized functional keyword: '%s'" % functional
415 raise ValueError(mess)
417 lsauthorslower = [a.lower() for a in self.allowed_xc[functional]]
418 if authors.lower() not in lsauthorslower:
419 mess = "Unrecognized authors keyword for %s: '%s'"
420 raise ValueError(mess % (functional, authors))
422 elif xc in self.allowed_xc:
423 functional = xc
424 authors = self.allowed_xc[xc][0]
425 else:
426 found = False
427 for key, value in self.allowed_xc.items():
428 if xc in value:
429 found = True
430 functional = key
431 authors = xc
432 break
434 if not found:
435 raise ValueError("Unrecognized 'xc' keyword: '%s'" % xc)
436 kwargs['xc'] = (functional, authors)
438 # Check fdf_arguments.
439 if kwargs['fdf_arguments'] is None:
440 kwargs['fdf_arguments'] = {}
442 if not isinstance(kwargs['fdf_arguments'], dict):
443 raise TypeError("fdf_arguments must be a dictionary.")
445 # Call baseclass.
446 FileIOCalculator.set(self, **kwargs)
448 def set_fdf_arguments(self, fdf_arguments):
449 """ Set the fdf_arguments after the initialization of the
450 calculator.
451 """
452 self.validate_fdf_arguments(fdf_arguments)
453 FileIOCalculator.set(self, fdf_arguments=fdf_arguments)
455 def validate_fdf_arguments(self, fdf_arguments):
456 """ Raises error if the fdf_argument input is not a
457 dictionary of allowed keys.
458 """
459 # None is valid
460 if fdf_arguments is None:
461 return
463 # Type checking.
464 if not isinstance(fdf_arguments, dict):
465 raise TypeError("fdf_arguments must be a dictionary.")
467 def calculate(self,
468 atoms=None,
469 properties=['energy'],
470 system_changes=all_changes):
471 """Capture the RuntimeError from FileIOCalculator.calculate
472 and add a little debug information from the Siesta output.
474 See base FileIocalculator for documentation.
475 """
477 FileIOCalculator.calculate(
478 self,
479 atoms=atoms,
480 properties=properties,
481 system_changes=system_changes)
483 # The below snippet would run if calculate() failed but I have
484 # disabled it for now since it looks to be just for debugging.
485 # --askhl
486 """
487 # Here a test to check if the potential are in the right place!!!
488 except RuntimeError as e:
489 try:
490 fname = os.path.join(self.directory, self.label+'.out')
491 with open(fname, 'r') as fd:
492 lines = fd.readlines()
493 debug_lines = 10
494 print('##### %d last lines of the Siesta output' % debug_lines)
495 for line in lines[-20:]:
496 print(line.strip())
497 print('##### end of siesta output')
498 raise e
499 except:
500 raise e
501 """
503 def write_input(self, atoms, properties=None, system_changes=None):
504 """Write input (fdf)-file.
505 See calculator.py for further details.
507 Parameters:
508 - atoms : The Atoms object to write.
509 - properties : The properties which should be calculated.
510 - system_changes : List of properties changed since last run.
511 """
512 # Call base calculator.
513 FileIOCalculator.write_input(
514 self,
515 atoms=atoms,
516 properties=properties,
517 system_changes=system_changes)
519 if system_changes is None and properties is None:
520 return
522 filename = self.getpath(ext='fdf')
524 # On any changes, remove all analysis files.
525 if system_changes is not None:
526 self.remove_analysis()
528 # Start writing the file.
529 with open(filename, 'w') as fd:
530 # Write system name and label.
531 fd.write(format_fdf('SystemName', self.prefix))
532 fd.write(format_fdf('SystemLabel', self.prefix))
533 fd.write("\n")
535 # Write explicitly given options first to
536 # allow the user to override anything.
537 fdf_arguments = self['fdf_arguments']
538 keys = sorted(fdf_arguments.keys())
539 for key in keys:
540 fd.write(format_fdf(key, fdf_arguments[key]))
542 # Force siesta to return error on no convergence.
543 # as default consistent with ASE expectations.
544 if 'SCFMustConverge' not in fdf_arguments.keys():
545 fd.write(format_fdf('SCFMustConverge', True))
546 fd.write("\n")
548 # Write spin level.
549 fd.write(format_fdf('Spin ', self['spin']))
550 # Spin backwards compatibility.
551 if self['spin'] == 'collinear':
552 fd.write(format_fdf('SpinPolarized', (True, "# Backwards compatibility.")))
553 elif self['spin'] == 'non-collinear':
554 fd.write(format_fdf('NonCollinear', (True, "# Backwards compatibility.")))
556 # Write functional.
557 functional, authors = self['xc']
558 fd.write(format_fdf('XC.functional', functional))
559 fd.write(format_fdf('XC.authors', authors))
560 fd.write("\n")
562 # Write mesh cutoff and energy shift.
563 fd.write(format_fdf('MeshCutoff',
564 (self['mesh_cutoff'], 'eV')))
565 fd.write(format_fdf('PAO.EnergyShift',
566 (self['energy_shift'], 'eV')))
567 fd.write("\n")
569 # Write the minimal arg
570 self._write_species(fd, atoms)
571 self._write_structure(fd, atoms)
573 # Use the saved density matrix if only 'cell' and 'positions'
574 # have changed.
575 if (system_changes is None or
576 ('numbers' not in system_changes and
577 'initial_magmoms' not in system_changes and
578 'initial_charges' not in system_changes)):
579 fd.write(format_fdf('DM.UseSaveDM', True))
581 # Save density.
582 if 'density' in properties:
583 fd.write(format_fdf('SaveRho', True))
585 self._write_kpts(fd)
587 if self['bandpath'] is not None:
588 lines = bandpath2bandpoints(self['bandpath'])
589 fd.write(lines)
590 fd.write('\n')
592 def read(self, filename):
593 """Read structural parameters from file .XV file
594 Read other results from other files
595 filename : siesta.XV
596 """
598 fname = self.getpath(filename)
599 if not os.path.exists(fname):
600 raise ReadError("The restart file '%s' does not exist" % fname)
601 with open(fname) as fd:
602 self.atoms = read_siesta_xv(fd)
603 self.read_results()
605 def getpath(self, fname=None, ext=None):
606 """ Returns the directory/fname string """
607 if fname is None:
608 fname = self.prefix
609 if ext is not None:
610 fname = '{}.{}'.format(fname, ext)
611 return os.path.join(self.directory, fname)
613 def remove_analysis(self):
614 """ Remove all analysis files"""
615 filename = self.getpath(ext='RHO')
616 if os.path.exists(filename):
617 os.remove(filename)
619 def _write_structure(self, fd, atoms):
620 """Translate the Atoms object to fdf-format.
622 Parameters:
623 - f: An open file object.
624 - atoms: An atoms object.
625 """
626 cell = atoms.cell
627 fd.write('\n')
629 if cell.rank in [1, 2]:
630 raise ValueError('Expected 3D unit cell or no unit cell. You may '
631 'wish to add vacuum along some directions.')
633 # Write lattice vectors
634 if np.any(cell):
635 fd.write(format_fdf('LatticeConstant', '1.0 Ang'))
636 fd.write('%block LatticeVectors\n')
637 for i in range(3):
638 for j in range(3):
639 s = (' %.15f' % cell[i, j]).rjust(16) + ' '
640 fd.write(s)
641 fd.write('\n')
642 fd.write('%endblock LatticeVectors\n')
643 fd.write('\n')
645 self._write_atomic_coordinates(fd, atoms)
647 # Write magnetic moments.
648 magmoms = atoms.get_initial_magnetic_moments()
650 # The DM.InitSpin block must be written to initialize to
651 # no spin. SIESTA default is FM initialization, if the
652 # block is not written, but we must conform to the
653 # atoms object.
654 if magmoms is not None:
655 if len(magmoms) == 0:
656 fd.write('#Empty block forces ASE initialization.\n')
658 fd.write('%block DM.InitSpin\n')
659 if len(magmoms) != 0 and isinstance(magmoms[0], np.ndarray):
660 for n, M in enumerate(magmoms):
661 if M[0] != 0:
662 fd.write(' %d %.14f %.14f %.14f \n' % (n + 1, M[0], M[1], M[2]))
663 elif len(magmoms) != 0 and isinstance(magmoms[0], float):
664 for n, M in enumerate(magmoms):
665 if M != 0:
666 fd.write(' %d %.14f \n' % (n + 1, M))
667 fd.write('%endblock DM.InitSpin\n')
668 fd.write('\n')
670 def _write_atomic_coordinates(self, fd, atoms):
671 """Write atomic coordinates.
673 Parameters:
674 - f: An open file object.
675 - atoms: An atoms object.
676 """
677 af = self.parameters.atomic_coord_format.lower()
678 if af == 'xyz':
679 self._write_atomic_coordinates_xyz(fd, atoms)
680 elif af == 'zmatrix':
681 self._write_atomic_coordinates_zmatrix(fd, atoms)
682 else:
683 raise RuntimeError('Unknown atomic_coord_format: {}'.format(af))
685 def _write_atomic_coordinates_xyz(self, fd, atoms):
686 """Write atomic coordinates.
688 Parameters:
689 - f: An open file object.
690 - atoms: An atoms object.
691 """
692 species, species_numbers = self.species(atoms)
693 fd.write('\n')
694 fd.write('AtomicCoordinatesFormat Ang\n')
695 fd.write('%block AtomicCoordinatesAndAtomicSpecies\n')
696 for atom, number in zip(atoms, species_numbers):
697 xyz = atom.position
698 line = (' %.9f' % xyz[0]).rjust(16) + ' '
699 line += (' %.9f' % xyz[1]).rjust(16) + ' '
700 line += (' %.9f' % xyz[2]).rjust(16) + ' '
701 line += str(number) + '\n'
702 fd.write(line)
703 fd.write('%endblock AtomicCoordinatesAndAtomicSpecies\n')
704 fd.write('\n')
706 origin = tuple(-atoms.get_celldisp().flatten())
707 if any(origin):
708 fd.write('%block AtomicCoordinatesOrigin\n')
709 fd.write(' %.4f %.4f %.4f\n' % origin)
710 fd.write('%endblock AtomicCoordinatesOrigin\n')
711 fd.write('\n')
713 def _write_atomic_coordinates_zmatrix(self, fd, atoms):
714 """Write atomic coordinates in Z-matrix format.
716 Parameters:
717 - f: An open file object.
718 - atoms: An atoms object.
719 """
720 species, species_numbers = self.species(atoms)
721 fd.write('\n')
722 fd.write('ZM.UnitsLength Ang\n')
723 fd.write('%block Zmatrix\n')
724 fd.write(' cartesian\n')
725 fstr = "{:5d}" + "{:20.10f}" * 3 + "{:3d}" * 3 + "{:7d} {:s}\n"
726 a2constr = self.make_xyz_constraints(atoms)
727 a2p, a2s = atoms.get_positions(), atoms.get_chemical_symbols()
728 for ia, (sp, xyz, ccc, sym) in enumerate(zip(species_numbers,
729 a2p,
730 a2constr,
731 a2s)):
732 fd.write(fstr.format(
733 sp, xyz[0], xyz[1], xyz[2], ccc[0],
734 ccc[1], ccc[2], ia + 1, sym))
735 fd.write('%endblock Zmatrix\n')
737 origin = tuple(-atoms.get_celldisp().flatten())
738 if any(origin):
739 fd.write('%block AtomicCoordinatesOrigin\n')
740 fd.write(' %.4f %.4f %.4f\n' % origin)
741 fd.write('%endblock AtomicCoordinatesOrigin\n')
742 fd.write('\n')
744 def make_xyz_constraints(self, atoms):
745 """ Create coordinate-resolved list of constraints [natoms, 0:3]
746 The elements of the list must be integers 0 or 1
747 1 -- means that the coordinate will be updated during relaxation
748 0 -- mains that the coordinate will be fixed during relaxation
749 """
750 from ase.constraints import FixAtoms, FixedLine, FixedPlane
751 import sys
752 import warnings
754 a = atoms
755 a2c = np.ones((len(a), 3), dtype=int)
756 for c in a.constraints:
757 if isinstance(c, FixAtoms):
758 a2c[c.get_indices()] = 0
759 elif isinstance(c, FixedLine):
760 norm_dir = c.dir / np.linalg.norm(c.dir)
761 if (max(norm_dir) - 1.0) > 1e-6:
762 raise RuntimeError(
763 'norm_dir: {} -- must be one of the Cartesian axes...'
764 .format(norm_dir))
765 a2c[c.a] = norm_dir.round().astype(int)
766 elif isinstance(c, FixedPlane):
767 norm_dir = c.dir / np.linalg.norm(c.dir)
768 if (max(norm_dir) - 1.0) > 1e-6:
769 raise RuntimeError(
770 'norm_dir: {} -- must be one of the Cartesian axes...'
771 .format(norm_dir))
772 a2c[c.a] = abs(1 - norm_dir.round().astype(int))
773 else:
774 warnings.warn('Constraint {} is ignored at {}'
775 .format(str(c), sys._getframe().f_code))
776 return a2c
778 def _write_kpts(self, fd):
779 """Write kpts.
781 Parameters:
782 - f : Open filename.
783 """
784 if self["kpts"] is None:
785 return
786 kpts = np.array(self['kpts'])
787 fd.write('\n')
788 fd.write('#KPoint grid\n')
789 fd.write('%block kgrid_Monkhorst_Pack\n')
791 for i in range(3):
792 s = ''
793 if i < len(kpts):
794 number = kpts[i]
795 displace = 0.0
796 else:
797 number = 1
798 displace = 0
799 for j in range(3):
800 if j == i:
801 write_this = number
802 else:
803 write_this = 0
804 s += ' %d ' % write_this
805 s += '%1.1f\n' % displace
806 fd.write(s)
807 fd.write('%endblock kgrid_Monkhorst_Pack\n')
808 fd.write('\n')
810 def _write_species(self, fd, atoms):
811 """Write input related the different species.
813 Parameters:
814 - f: An open file object.
815 - atoms: An atoms object.
816 """
817 species, species_numbers = self.species(atoms)
819 if self['pseudo_path'] is not None:
820 pseudo_path = self['pseudo_path']
821 elif 'SIESTA_PP_PATH' in os.environ:
822 pseudo_path = os.environ['SIESTA_PP_PATH']
823 else:
824 mess = "Please set the environment variable 'SIESTA_PP_PATH'"
825 raise Exception(mess)
827 fd.write(format_fdf('NumberOfSpecies', len(species)))
828 fd.write(format_fdf('NumberOfAtoms', len(atoms)))
830 pao_basis = []
831 chemical_labels = []
832 basis_sizes = []
833 synth_blocks = []
834 for species_number, spec in enumerate(species):
835 species_number += 1
836 symbol = spec['symbol']
837 atomic_number = atomic_numbers[symbol]
839 if spec['pseudopotential'] is None:
840 if self.pseudo_qualifier() == '':
841 label = symbol
842 pseudopotential = label + '.psf'
843 else:
844 label = '.'.join([symbol, self.pseudo_qualifier()])
845 pseudopotential = label + '.psf'
846 else:
847 pseudopotential = spec['pseudopotential']
848 label = os.path.basename(pseudopotential)
849 label = '.'.join(label.split('.')[:-1])
851 if not os.path.isabs(pseudopotential):
852 pseudopotential = join(pseudo_path, pseudopotential)
854 if not os.path.exists(pseudopotential):
855 mess = "Pseudopotential '%s' not found" % pseudopotential
856 raise RuntimeError(mess)
858 name = os.path.basename(pseudopotential)
859 name = name.split('.')
860 name.insert(-1, str(species_number))
861 if spec['ghost']:
862 name.insert(-1, 'ghost')
863 atomic_number = -atomic_number
865 name = '.'.join(name)
866 pseudo_targetpath = self.getpath(name)
868 if join(os.getcwd(), name) != pseudopotential:
869 if islink(pseudo_targetpath) or isfile(pseudo_targetpath):
870 os.remove(pseudo_targetpath)
871 symlink_pseudos = self['symlink_pseudos']
873 if symlink_pseudos is None:
874 symlink_pseudos = not os.name == 'nt'
876 if symlink_pseudos:
877 os.symlink(pseudopotential, pseudo_targetpath)
878 else:
879 shutil.copy(pseudopotential, pseudo_targetpath)
881 if not spec['excess_charge'] is None:
882 atomic_number += 200
883 n_atoms = sum(np.array(species_numbers) == species_number)
885 paec = float(spec['excess_charge']) / n_atoms
886 vc = get_valence_charge(pseudopotential)
887 fraction = float(vc + paec) / vc
888 pseudo_head = name[:-4]
889 fractional_command = os.environ['SIESTA_UTIL_FRACTIONAL']
890 cmd = '%s %s %.7f' % (fractional_command,
891 pseudo_head,
892 fraction)
893 os.system(cmd)
895 pseudo_head += '-Fraction-%.5f' % fraction
896 synth_pseudo = pseudo_head + '.psf'
897 synth_block_filename = pseudo_head + '.synth'
898 os.remove(name)
899 shutil.copyfile(synth_pseudo, name)
900 synth_block = read_vca_synth_block(
901 synth_block_filename,
902 species_number=species_number)
903 synth_blocks.append(synth_block)
905 if len(synth_blocks) > 0:
906 fd.write(format_fdf('SyntheticAtoms', list(synth_blocks)))
908 label = '.'.join(np.array(name.split('.'))[:-1])
909 string = ' %d %d %s' % (species_number, atomic_number, label)
910 chemical_labels.append(string)
911 if isinstance(spec['basis_set'], PAOBasisBlock):
912 pao_basis.append(spec['basis_set'].script(label))
913 else:
914 basis_sizes.append((" " + label, spec['basis_set']))
915 fd.write((format_fdf('ChemicalSpecieslabel', chemical_labels)))
916 fd.write('\n')
917 fd.write((format_fdf('PAO.Basis', pao_basis)))
918 fd.write((format_fdf('PAO.BasisSizes', basis_sizes)))
919 fd.write('\n')
921 def pseudo_qualifier(self):
922 """Get the extra string used in the middle of the pseudopotential.
923 The retrieved pseudopotential for a specific element will be
924 'H.xxx.psf' for the element 'H' with qualifier 'xxx'. If qualifier
925 is set to None then the qualifier is set to functional name.
926 """
927 if self['pseudo_qualifier'] is None:
928 return self['xc'][0].lower()
929 else:
930 return self['pseudo_qualifier']
932 def read_results(self):
933 """Read the results.
934 """
935 self.read_number_of_grid_points()
936 self.read_energy()
937 self.read_forces_stress()
938 self.read_eigenvalues()
939 self.read_kpoints()
940 self.read_dipole()
941 self.read_pseudo_density()
942 self.read_hsx()
943 self.read_dim()
944 if self.results['hsx'] is not None:
945 self.read_pld(self.results['hsx'].norbitals,
946 len(self.atoms))
947 self.atoms.cell = self.results['pld'].cell * Bohr
948 else:
949 self.results['pld'] = None
951 self.read_wfsx()
952 self.read_ion(self.atoms)
954 self.read_bands()
956 def read_bands(self):
957 bandpath = self['bandpath']
958 if bandpath is None:
959 return
961 if len(bandpath.kpts) < 1:
962 return
964 fname = self.getpath(ext='bands')
965 with open(fname) as fd:
966 kpts, energies, efermi = read_bands_file(fd)
967 bs = resolve_band_structure(bandpath, kpts, energies, efermi)
968 self.results['bandstructure'] = bs
970 def band_structure(self):
971 return self.results['bandstructure']
973 def read_ion(self, atoms):
974 """
975 Read the ion.xml file of each specie
976 """
977 from ase.calculators.siesta.import_ion_xml import get_ion
979 species, species_numbers = self.species(atoms)
981 self.results['ion'] = {}
982 for species_number, spec in enumerate(species):
983 species_number += 1
985 symbol = spec['symbol']
986 atomic_number = atomic_numbers[symbol]
988 if spec['pseudopotential'] is None:
989 if self.pseudo_qualifier() == '':
990 label = symbol
991 else:
992 label = '.'.join([symbol, self.pseudo_qualifier()])
993 pseudopotential = self.getpath(label, 'psf')
994 else:
995 pseudopotential = spec['pseudopotential']
996 label = os.path.basename(pseudopotential)
997 label = '.'.join(label.split('.')[:-1])
999 name = os.path.basename(pseudopotential)
1000 name = name.split('.')
1001 name.insert(-1, str(species_number))
1002 if spec['ghost']:
1003 name.insert(-1, 'ghost')
1004 atomic_number = -atomic_number
1005 name = '.'.join(name)
1007 label = '.'.join(np.array(name.split('.'))[:-1])
1009 if label not in self.results['ion']:
1010 fname = self.getpath(label, 'ion.xml')
1011 if os.path.isfile(fname):
1012 self.results['ion'][label] = get_ion(fname)
1014 def read_hsx(self):
1015 """
1016 Read the siesta HSX file.
1017 return a namedtuple with the following arguments:
1018 'norbitals', 'norbitals_sc', 'nspin', 'nonzero',
1019 'is_gamma', 'sc_orb2uc_orb', 'row2nnzero', 'sparse_ind2column',
1020 'H_sparse', 'S_sparse', 'aB2RaB_sparse', 'total_elec_charge', 'temp'
1021 """
1022 from ase.calculators.siesta.import_functions import readHSX
1024 filename = self.getpath(ext='HSX')
1025 if isfile(filename):
1026 self.results['hsx'] = readHSX(filename)
1027 else:
1028 self.results['hsx'] = None
1030 def read_dim(self):
1031 """
1032 Read the siesta DIM file
1033 Retrun a namedtuple with the following arguments:
1034 'natoms_sc', 'norbitals_sc', 'norbitals', 'nspin',
1035 'nnonzero', 'natoms_interacting'
1036 """
1037 from ase.calculators.siesta.import_functions import readDIM
1039 filename = self.getpath(ext='DIM')
1040 if isfile(filename):
1041 self.results['dim'] = readDIM(filename)
1042 else:
1043 self.results['dim'] = None
1045 def read_pld(self, norb, natms):
1046 """
1047 Read the siesta PLD file
1048 Return a namedtuple with the following arguments:
1049 'max_rcut', 'orb2ao', 'orb2uorb', 'orb2occ', 'atm2sp',
1050 'atm2shift', 'coord_sc', 'cell', 'nunit_cells'
1051 """
1052 from ase.calculators.siesta.import_functions import readPLD
1054 filename = self.getpath(ext='PLD')
1055 if isfile(filename):
1056 self.results['pld'] = readPLD(filename, norb, natms)
1057 else:
1058 self.results['pld'] = None
1060 def read_wfsx(self):
1061 """
1062 Read the siesta WFSX file
1063 Return a namedtuple with the following arguments:
1064 """
1065 from ase.calculators.siesta.import_functions import readWFSX
1067 fname_woext = os.path.join(self.directory, self.prefix)
1069 if isfile(fname_woext + '.WFSX'):
1070 filename = fname_woext + '.WFSX'
1071 self.results['wfsx'] = readWFSX(filename)
1072 elif isfile(fname_woext + '.fullBZ.WFSX'):
1073 filename = fname_woext + '.fullBZ.WFSX'
1074 readWFSX(filename)
1075 self.results['wfsx'] = readWFSX(filename)
1076 else:
1077 self.results['wfsx'] = None
1079 def read_pseudo_density(self):
1080 """Read the density if it is there."""
1081 filename = self.getpath(ext='RHO')
1082 if isfile(filename):
1083 self.results['density'] = read_rho(filename)
1085 def read_number_of_grid_points(self):
1086 """Read number of grid points from SIESTA's text-output file. """
1088 fname = self.getpath(ext='out')
1089 with open(fname, 'r') as fd:
1090 for line in fd:
1091 line = line.strip().lower()
1092 if line.startswith('initmesh: mesh ='):
1093 n_points = [int(word) for word in line.split()[3:8:2]]
1094 self.results['n_grid_point'] = n_points
1095 break
1096 else:
1097 raise RuntimeError
1099 def read_energy(self):
1100 """Read energy from SIESTA's text-output file.
1101 """
1102 fname = self.getpath(ext='out')
1103 with open(fname, 'r') as fd:
1104 text = fd.read().lower()
1106 assert 'final energy' in text
1107 lines = iter(text.split('\n'))
1109 # Get the energy and free energy the last time it appears
1110 for line in lines:
1111 has_energy = line.startswith('siesta: etot =')
1112 if has_energy:
1113 self.results['energy'] = float(line.split()[-1])
1114 line = next(lines)
1115 self.results['free_energy'] = float(line.split()[-1])
1117 if ('energy' not in self.results or
1118 'free_energy' not in self.results):
1119 raise RuntimeError
1121 def read_forces_stress(self):
1122 """Read the forces and stress from the FORCE_STRESS file.
1123 """
1124 fname = self.getpath('FORCE_STRESS')
1125 with open(fname, 'r') as fd:
1126 lines = fd.readlines()
1128 stress_lines = lines[1:4]
1129 stress = np.empty((3, 3))
1130 for i in range(3):
1131 line = stress_lines[i].strip().split(' ')
1132 line = [s for s in line if len(s) > 0]
1133 stress[i] = [float(s) for s in line]
1135 self.results['stress'] = np.array(
1136 [stress[0, 0], stress[1, 1], stress[2, 2],
1137 stress[1, 2], stress[0, 2], stress[0, 1]])
1139 self.results['stress'] *= Ry / Bohr**3
1141 start = 5
1142 self.results['forces'] = np.zeros((len(lines) - start, 3), float)
1143 for i in range(start, len(lines)):
1144 line = [s for s in lines[i].strip().split(' ') if len(s) > 0]
1145 self.results['forces'][i - start] = [float(s) for s in line[2:5]]
1147 self.results['forces'] *= Ry / Bohr
1149 def read_eigenvalues(self):
1150 """ A robust procedure using the suggestion by Federico Marchesin """
1152 fname = self.getpath(ext='EIG')
1153 try:
1154 with open(fname, "r") as fd:
1155 self.results['fermi_energy'] = float(fd.readline())
1156 n, nspin, nkp = map(int, fd.readline().split())
1157 _ee = np.split(
1158 np.array(fd.read().split()).astype(float), nkp)
1159 except (IOError):
1160 return 1
1162 ksn2e = np.delete(_ee, 0, 1).reshape([nkp, nspin, n])
1164 eigarray = np.empty((nspin, nkp, n))
1165 eigarray[:] = np.inf
1167 for k, sn2e in enumerate(ksn2e):
1168 for s, n2e in enumerate(sn2e):
1169 eigarray[s, k, :] = n2e
1171 assert np.isfinite(eigarray).all()
1173 self.results['eigenvalues'] = eigarray
1174 return 0
1176 def read_kpoints(self):
1177 """ Reader of the .KP files """
1179 fname = self.getpath(ext='KP')
1180 try:
1181 with open(fname, "r") as fd:
1182 nkp = int(next(fd))
1183 kpoints = np.empty((nkp, 3))
1184 kweights = np.empty(nkp)
1186 for i in range(nkp):
1187 line = next(fd)
1188 tokens = line.split()
1189 numbers = np.array(tokens[1:]).astype(float)
1190 kpoints[i] = numbers[:3]
1191 kweights[i] = numbers[3]
1193 except (IOError):
1194 return 1
1196 self.results['kpoints'] = kpoints
1197 self.results['kweights'] = kweights
1199 return 0
1201 def read_dipole(self):
1202 """Read dipole moment. """
1203 dipole = np.zeros([1, 3])
1204 with open(self.getpath(ext='out'), 'r') as fd:
1205 for line in fd:
1206 if line.rfind('Electric dipole (Debye)') > -1:
1207 dipole = np.array([float(f) for f in line.split()[5:8]])
1208 # debye to e*Ang
1209 self.results['dipole'] = dipole * 0.2081943482534
1211 def get_fermi_level(self):
1212 return self.results['fermi_energy']
1214 def get_k_point_weights(self):
1215 return self.results['kweights']
1217 def get_ibz_k_points(self):
1218 return self.results['kpoints']
1221class Siesta3_2(Siesta):
1222 def __init__(self, *args, **kwargs):
1223 warnings.warn(
1224 "The Siesta3_2 calculator class will no longer be supported. "
1225 "Use 'ase.calculators.siesta.Siesta in stead. "
1226 "If using the ASE interface with SIESTA 3.2 you must explicitly "
1227 "include the keywords 'SpinPolarized', 'NonCollinearSpin' and "
1228 "'SpinOrbit' if needed.",
1229 np.VisibleDeprecationWarning)
1230 Siesta.__init__(self, *args, **kwargs)