"""
This module defines the ASE interface to SIESTA.
Written by Mads Engelund (see www.espeem.com)
Home of the SIESTA package:
http://www.uam.es/departamentos/ciencias/fismateriac/siesta
2017.04 - Pedro Brandimarte: changes for python 2-3 compatible
"""
import os
import re
import tempfile
import warnings
import shutil
from os.path import join, isfile, islink
import numpy as np
from ase.units import Ry, eV, Bohr
from ase.data import atomic_numbers
from ase.io.siesta import read_siesta_xv
from ase.calculators.siesta.import_functions import read_rho
from ase.calculators.siesta.import_functions import \
get_valence_charge, read_vca_synth_block
from ase.calculators.calculator import FileIOCalculator, ReadError
from ase.calculators.calculator import Parameters, all_changes
from ase.calculators.siesta.parameters import PAOBasisBlock, Species
from ase.calculators.siesta.parameters import format_fdf
meV = 0.001 * eV
def parse_siesta_version(output: bytes) -> str:
match = re.search(rb'Siesta Version\s*:\s*(\S+)', output)
if match is None:
raise RuntimeError('Could not get Siesta version info from output '
'{!r}'.format(output))
string = match.group(1).decode('ascii')
return string
def get_siesta_version(executable: str) -> str:
""" Return SIESTA version number.
Run the command, for instance 'siesta' and
then parse the output in order find the
version number.
"""
# XXX We need a test of this kind of function. But Siesta().command
# is not enough to tell us how to run Siesta, because it could contain
# all sorts of mpirun and other weird parts.
temp_dirname = tempfile.mkdtemp(prefix='siesta-version-check-')
try:
from subprocess import Popen, PIPE
proc = Popen([executable],
stdin=PIPE,
stdout=PIPE,
stderr=PIPE,
cwd=temp_dirname)
output, _ = proc.communicate()
# We are not providing any input, so Siesta will give us a failure
# saying that it has no Chemical_species_label and exit status 1
# (as of siesta-4.1-b4)
finally:
shutil.rmtree(temp_dirname)
return parse_siesta_version(output)
def bandpath2bandpoints(path):
lines = []
add = lines.append
add('BandLinesScale ReciprocalLatticeVectors\n')
add('%block BandPoints\n')
for kpt in path.kpts:
add(' {:18.15f} {:18.15f} {:18.15f}\n'.format(*kpt))
add('%endblock BandPoints')
return ''.join(lines)
def read_bands_file(fd):
efermi = float(next(fd))
next(fd) # Appears to be max/min energy. Not important for us
header = next(fd) # Array shape: nbands, nspins, nkpoints
nbands, nspins, nkpts = np.array(header.split()).astype(int)
# three fields for kpt coords, then all the energies
ntokens = nbands * nspins + 3
# Read energies for each kpoint:
data = []
for i in range(nkpts):
line = next(fd)
tokens = line.split()
while len(tokens) < ntokens:
# Multirow table. Keep adding lines until the table ends,
# which should happen exactly when we have all the energies
# for this kpoint.
line = next(fd)
tokens += line.split()
assert len(tokens) == ntokens
values = np.array(tokens).astype(float)
data.append(values)
data = np.array(data)
assert len(data) == nkpts
kpts = data[:, :3]
energies = data[:, 3:]
energies = energies.reshape(nkpts, nspins, nbands)
assert energies.shape == (nkpts, nspins, nbands)
return kpts, energies, efermi
def resolve_band_structure(path, kpts, energies, efermi):
"""Convert input BandPath along with Siesta outputs into BS object."""
# Right now this function doesn't do much.
#
# Not sure how the output kpoints in the siesta.bands file are derived.
# They appear to be related to the lattice parameter.
#
# We should verify that they are consistent with our input path,
# but since their meaning is unclear, we can't quite do so.
#
# Also we should perhaps verify the cell. If we had the cell, we
# could construct the bandpath from scratch (i.e., pure outputs).
from ase.spectrum.band_structure import BandStructure
ksn2e = energies
skn2e = np.swapaxes(ksn2e, 0, 1)
bs = BandStructure(path, skn2e, reference=efermi)
return bs
class SiestaParameters(Parameters):
"""Parameters class for the calculator.
Documented in BaseSiesta.__init__
"""
def __init__(
self,
label='siesta',
mesh_cutoff=200 * Ry,
energy_shift=100 * meV,
kpts=None,
xc='LDA',
basis_set='DZP',
spin='non-polarized',
species=tuple(),
pseudo_qualifier=None,
pseudo_path=None,
symlink_pseudos=None,
atoms=None,
restart=None,
fdf_arguments=None,
atomic_coord_format='xyz',
bandpath=None):
kwargs = locals()
kwargs.pop('self')
Parameters.__init__(self, **kwargs)
[docs]class Siesta(FileIOCalculator):
"""Calculator interface to the SIESTA code.
"""
# Siesta manual does not document many of the basis names.
# basis_specs.f has a ton of aliases for each.
# Let's just list one of each type then.
#
# Maybe we should be less picky about these keyword names.
allowed_basis_names = ['SZ', 'SZP',
'DZ', 'DZP', 'DZP2',
'TZ', 'TZP', 'TZP2', 'TZP3']
allowed_spins = ['non-polarized', 'collinear',
'non-collinear', 'spin-orbit']
allowed_xc = {
'LDA': ['PZ', 'CA', 'PW92'],
'GGA': ['PW91', 'PBE', 'revPBE', 'RPBE',
'WC', 'AM05', 'PBEsol', 'PBEJsJrLO',
'PBEGcGxLO', 'PBEGcGxHEG', 'BLYP'],
'VDW': ['DRSLL', 'LMKLL', 'KBM', 'C09', 'BH', 'VV']}
name = 'siesta'
command = 'siesta < PREFIX.fdf > PREFIX.out'
implemented_properties = [
'energy',
'forces',
'stress',
'dipole',
'eigenvalues',
'density',
'fermi_energy']
# Dictionary of valid input vaiables.
default_parameters = SiestaParameters()
# XXX Not a ASE standard mechanism (yet). We need to communicate to
# ase.spectrum.band_structure.calculate_band_structure() that we expect
# it to use the bandpath keyword.
accepts_bandpath_keyword = True
def __init__(self, command=None, **kwargs):
"""ASE interface to the SIESTA code.
Parameters:
- label : The basename of all files created during
calculation.
- mesh_cutoff : Energy in eV.
The mesh cutoff energy for determining number of
grid points in the matrix-element calculation.
- energy_shift : Energy in eV
The confining energy of the basis set generation.
- kpts : Tuple of 3 integers, the k-points in different
directions.
- xc : The exchange-correlation potential. Can be set to
any allowed value for either the Siesta
XC.funtional or XC.authors keyword. Default "LDA"
- basis_set : "SZ"|"SZP"|"DZ"|"DZP"|"TZP", strings which specify
the type of functions basis set.
- spin : "non-polarized"|"collinear"|
"non-collinear|spin-orbit".
The level of spin description to be used.
- species : None|list of Species objects. The species objects
can be used to to specify the basis set,
pseudopotential and whether the species is ghost.
The tag on the atoms object and the element is used
together to identify the species.
- pseudo_path : None|path. This path is where
pseudopotentials are taken from.
If None is given, then then the path given
in $SIESTA_PP_PATH will be used.
- pseudo_qualifier: None|string. This string will be added to the
pseudopotential path that will be retrieved.
For hydrogen with qualifier "abc" the
pseudopotential "H.abc.psf" will be retrieved.
- symlink_pseudos: None|bool
If true, symlink pseudopotentials
into the calculation directory, else copy them.
Defaults to true on Unix and false on Windows.
- atoms : The Atoms object.
- restart : str. Prefix for restart file.
May contain a directory.
Default is None, don't restart.
- fdf_arguments: Explicitly given fdf arguments. Dictonary using
Siesta keywords as given in the manual. List values
are written as fdf blocks with each element on a
separate line, while tuples will write each element
in a single line. ASE units are assumed in the
input.
- atomic_coord_format: "xyz"|"zmatrix", strings to switch between
the default way of entering the system's geometry
(via the block AtomicCoordinatesAndAtomicSpecies)
and a recent method via the block Zmatrix. The
block Zmatrix allows to specify basic geometry
constrains such as realized through the ASE classes
FixAtom, FixedLine and FixedPlane.
"""
# Put in the default arguments.
parameters = self.default_parameters.__class__(**kwargs)
# Call the base class.
FileIOCalculator.__init__(
self,
command=command,
**parameters)
# For compatibility with old variable name:
commandvar = os.environ.get('SIESTA_COMMAND')
if commandvar is not None:
warnings.warn('Please use $ASE_SIESTA_COMMAND and not '
'$SIESTA_COMMAND, which will be ignored '
'in the future. The new command format will not '
'work with the "<%s > %s" specification. Use '
'instead e.g. "ASE_SIESTA_COMMAND=siesta'
' < PREFIX.fdf > PREFIX.out", where PREFIX will '
'automatically be replaced by calculator label',
np.VisibleDeprecationWarning)
runfile = self.prefix + '.fdf'
outfile = self.prefix + '.out'
try:
self.command = commandvar % (runfile, outfile)
except TypeError:
raise ValueError(
"The 'SIESTA_COMMAND' environment must " +
"be a format string" +
" with two string arguments.\n" +
"Example : 'siesta < %s > %s'.\n" +
"Got '%s'" % commandvar)
def __getitem__(self, key):
"""Convenience method to retrieve a parameter as
calculator[key] rather than calculator.parameters[key]
Parameters:
-key : str, the name of the parameters to get.
"""
return self.parameters[key]
def species(self, atoms):
"""Find all relevant species depending on the atoms object and
species input.
Parameters :
- atoms : An Atoms object.
"""
# For each element use default species from the species input, or set
# up a default species from the general default parameters.
symbols = np.array(atoms.get_chemical_symbols())
tags = atoms.get_tags()
species = list(self['species'])
default_species = [
s for s in species
if (s['tag'] is None) and s['symbol'] in symbols]
default_symbols = [s['symbol'] for s in default_species]
for symbol in symbols:
if symbol not in default_symbols:
spec = Species(symbol=symbol,
basis_set=self['basis_set'],
tag=None)
default_species.append(spec)
default_symbols.append(symbol)
assert len(default_species) == len(np.unique(symbols))
# Set default species as the first species.
species_numbers = np.zeros(len(atoms), int)
i = 1
for spec in default_species:
mask = symbols == spec['symbol']
species_numbers[mask] = i
i += 1
# Set up the non-default species.
non_default_species = [s for s in species if not s['tag'] is None]
for spec in non_default_species:
mask1 = (tags == spec['tag'])
mask2 = (symbols == spec['symbol'])
mask = np.logical_and(mask1, mask2)
if sum(mask) > 0:
species_numbers[mask] = i
i += 1
all_species = default_species + non_default_species
return all_species, species_numbers
def set(self, **kwargs):
"""Set all parameters.
Parameters:
-kwargs : Dictionary containing the keywords defined in
SiestaParameters.
"""
# XXX Inserted these next few lines because set() would otherwise
# discard all previously set keywords to their defaults! --askhl
current = self.parameters.copy()
current.update(kwargs)
kwargs = current
# Find not allowed keys.
default_keys = list(self.__class__.default_parameters)
offending_keys = set(kwargs) - set(default_keys)
if len(offending_keys) > 0:
mess = "'set' does not take the keywords: %s "
raise ValueError(mess % list(offending_keys))
# Use the default parameters.
parameters = self.__class__.default_parameters.copy()
parameters.update(kwargs)
kwargs = parameters
# Check energy inputs.
for arg in ['mesh_cutoff', 'energy_shift']:
value = kwargs.get(arg)
if value is None:
continue
if not (isinstance(value, (float, int)) and value > 0):
mess = "'%s' must be a positive number(in eV), \
got '%s'" % (arg, value)
raise ValueError(mess)
# Check the basis set input.
if 'basis_set' in kwargs:
basis_set = kwargs['basis_set']
allowed = self.allowed_basis_names
if not (isinstance(basis_set, PAOBasisBlock) or
basis_set in allowed):
mess = "Basis must be either %s, got %s" % (allowed, basis_set)
raise ValueError(mess)
# Check the spin input.
if 'spin' in kwargs:
if kwargs['spin'] == 'UNPOLARIZED':
warnings.warn("The keyword 'UNPOLARIZED' is deprecated,"
"and replaced by 'non-polarized'",
np.VisibleDeprecationWarning)
kwargs['spin'] = 'non-polarized'
spin = kwargs['spin']
if spin is not None and (spin.lower() not in self.allowed_spins):
mess = "Spin must be %s, got '%s'" % (self.allowed_spins, spin)
raise ValueError(mess)
# Check the functional input.
xc = kwargs.get('xc', 'LDA')
if isinstance(xc, (tuple, list)) and len(xc) == 2:
functional, authors = xc
if functional.lower() not in [k.lower() for k in self.allowed_xc]:
mess = "Unrecognized functional keyword: '%s'" % functional
raise ValueError(mess)
lsauthorslower = [a.lower() for a in self.allowed_xc[functional]]
if authors.lower() not in lsauthorslower:
mess = "Unrecognized authors keyword for %s: '%s'"
raise ValueError(mess % (functional, authors))
elif xc in self.allowed_xc:
functional = xc
authors = self.allowed_xc[xc][0]
else:
found = False
for key, value in self.allowed_xc.items():
if xc in value:
found = True
functional = key
authors = xc
break
if not found:
raise ValueError("Unrecognized 'xc' keyword: '%s'" % xc)
kwargs['xc'] = (functional, authors)
# Check fdf_arguments.
if kwargs['fdf_arguments'] is None:
kwargs['fdf_arguments'] = {}
if not isinstance(kwargs['fdf_arguments'], dict):
raise TypeError("fdf_arguments must be a dictionary.")
# Call baseclass.
FileIOCalculator.set(self, **kwargs)
def set_fdf_arguments(self, fdf_arguments):
""" Set the fdf_arguments after the initialization of the
calculator.
"""
self.validate_fdf_arguments(fdf_arguments)
FileIOCalculator.set(self, fdf_arguments=fdf_arguments)
def validate_fdf_arguments(self, fdf_arguments):
""" Raises error if the fdf_argument input is not a
dictionary of allowed keys.
"""
# None is valid
if fdf_arguments is None:
return
# Type checking.
if not isinstance(fdf_arguments, dict):
raise TypeError("fdf_arguments must be a dictionary.")
def calculate(self,
atoms=None,
properties=['energy'],
system_changes=all_changes):
"""Capture the RuntimeError from FileIOCalculator.calculate
and add a little debug information from the Siesta output.
See base FileIocalculator for documentation.
"""
FileIOCalculator.calculate(
self,
atoms=atoms,
properties=properties,
system_changes=system_changes)
# The below snippet would run if calculate() failed but I have
# disabled it for now since it looks to be just for debugging.
# --askhl
"""
# Here a test to check if the potential are in the right place!!!
except RuntimeError as e:
try:
fname = os.path.join(self.directory, self.label+'.out')
with open(fname, 'r') as fd:
lines = fd.readlines()
debug_lines = 10
print('##### %d last lines of the Siesta output' % debug_lines)
for line in lines[-20:]:
print(line.strip())
print('##### end of siesta output')
raise e
except:
raise e
"""
def write_input(self, atoms, properties=None, system_changes=None):
"""Write input (fdf)-file.
See calculator.py for further details.
Parameters:
- atoms : The Atoms object to write.
- properties : The properties which should be calculated.
- system_changes : List of properties changed since last run.
"""
# Call base calculator.
FileIOCalculator.write_input(
self,
atoms=atoms,
properties=properties,
system_changes=system_changes)
if system_changes is None and properties is None:
return
filename = self.getpath(ext='fdf')
# On any changes, remove all analysis files.
if system_changes is not None:
self.remove_analysis()
# Start writing the file.
with open(filename, 'w') as fd:
# Write system name and label.
fd.write(format_fdf('SystemName', self.prefix))
fd.write(format_fdf('SystemLabel', self.prefix))
fd.write("\n")
# Write explicitly given options first to
# allow the user to override anything.
fdf_arguments = self['fdf_arguments']
keys = sorted(fdf_arguments.keys())
for key in keys:
fd.write(format_fdf(key, fdf_arguments[key]))
# Force siesta to return error on no convergence.
# as default consistent with ASE expectations.
if 'SCFMustConverge' not in fdf_arguments.keys():
fd.write(format_fdf('SCFMustConverge', True))
fd.write("\n")
# Write spin level.
fd.write(format_fdf('Spin ', self['spin']))
# Spin backwards compatibility.
if self['spin'] == 'collinear':
fd.write(format_fdf('SpinPolarized', (True, "# Backwards compatibility.")))
elif self['spin'] == 'non-collinear':
fd.write(format_fdf('NonCollinear', (True, "# Backwards compatibility.")))
# Write functional.
functional, authors = self['xc']
fd.write(format_fdf('XC.functional', functional))
fd.write(format_fdf('XC.authors', authors))
fd.write("\n")
# Write mesh cutoff and energy shift.
fd.write(format_fdf('MeshCutoff',
(self['mesh_cutoff'], 'eV')))
fd.write(format_fdf('PAO.EnergyShift',
(self['energy_shift'], 'eV')))
fd.write("\n")
# Write the minimal arg
self._write_species(fd, atoms)
self._write_structure(fd, atoms)
# Use the saved density matrix if only 'cell' and 'positions'
# have changed.
if (system_changes is None or
('numbers' not in system_changes and
'initial_magmoms' not in system_changes and
'initial_charges' not in system_changes)):
fd.write(format_fdf('DM.UseSaveDM', True))
# Save density.
if 'density' in properties:
fd.write(format_fdf('SaveRho', True))
self._write_kpts(fd)
if self['bandpath'] is not None:
lines = bandpath2bandpoints(self['bandpath'])
fd.write(lines)
fd.write('\n')
def read(self, filename):
"""Read structural parameters from file .XV file
Read other results from other files
filename : siesta.XV
"""
fname = self.getpath(filename)
if not os.path.exists(fname):
raise ReadError("The restart file '%s' does not exist" % fname)
with open(fname) as fd:
self.atoms = read_siesta_xv(fd)
self.read_results()
def getpath(self, fname=None, ext=None):
""" Returns the directory/fname string """
if fname is None:
fname = self.prefix
if ext is not None:
fname = '{}.{}'.format(fname, ext)
return os.path.join(self.directory, fname)
def remove_analysis(self):
""" Remove all analysis files"""
filename = self.getpath(ext='RHO')
if os.path.exists(filename):
os.remove(filename)
def _write_structure(self, fd, atoms):
"""Translate the Atoms object to fdf-format.
Parameters:
- f: An open file object.
- atoms: An atoms object.
"""
cell = atoms.cell
fd.write('\n')
if cell.rank in [1, 2]:
raise ValueError('Expected 3D unit cell or no unit cell. You may '
'wish to add vacuum along some directions.')
# Write lattice vectors
if np.any(cell):
fd.write(format_fdf('LatticeConstant', '1.0 Ang'))
fd.write('%block LatticeVectors\n')
for i in range(3):
for j in range(3):
s = (' %.15f' % cell[i, j]).rjust(16) + ' '
fd.write(s)
fd.write('\n')
fd.write('%endblock LatticeVectors\n')
fd.write('\n')
self._write_atomic_coordinates(fd, atoms)
# Write magnetic moments.
magmoms = atoms.get_initial_magnetic_moments()
# The DM.InitSpin block must be written to initialize to
# no spin. SIESTA default is FM initialization, if the
# block is not written, but we must conform to the
# atoms object.
if magmoms is not None:
if len(magmoms) == 0:
fd.write('#Empty block forces ASE initialization.\n')
fd.write('%block DM.InitSpin\n')
if len(magmoms) != 0 and isinstance(magmoms[0], np.ndarray):
for n, M in enumerate(magmoms):
if M[0] != 0:
fd.write(' %d %.14f %.14f %.14f \n' % (n + 1, M[0], M[1], M[2]))
elif len(magmoms) != 0 and isinstance(magmoms[0], float):
for n, M in enumerate(magmoms):
if M != 0:
fd.write(' %d %.14f \n' % (n + 1, M))
fd.write('%endblock DM.InitSpin\n')
fd.write('\n')
def _write_atomic_coordinates(self, fd, atoms):
"""Write atomic coordinates.
Parameters:
- f: An open file object.
- atoms: An atoms object.
"""
af = self.parameters.atomic_coord_format.lower()
if af == 'xyz':
self._write_atomic_coordinates_xyz(fd, atoms)
elif af == 'zmatrix':
self._write_atomic_coordinates_zmatrix(fd, atoms)
else:
raise RuntimeError('Unknown atomic_coord_format: {}'.format(af))
def _write_atomic_coordinates_xyz(self, fd, atoms):
"""Write atomic coordinates.
Parameters:
- f: An open file object.
- atoms: An atoms object.
"""
species, species_numbers = self.species(atoms)
fd.write('\n')
fd.write('AtomicCoordinatesFormat Ang\n')
fd.write('%block AtomicCoordinatesAndAtomicSpecies\n')
for atom, number in zip(atoms, species_numbers):
xyz = atom.position
line = (' %.9f' % xyz[0]).rjust(16) + ' '
line += (' %.9f' % xyz[1]).rjust(16) + ' '
line += (' %.9f' % xyz[2]).rjust(16) + ' '
line += str(number) + '\n'
fd.write(line)
fd.write('%endblock AtomicCoordinatesAndAtomicSpecies\n')
fd.write('\n')
origin = tuple(-atoms.get_celldisp().flatten())
if any(origin):
fd.write('%block AtomicCoordinatesOrigin\n')
fd.write(' %.4f %.4f %.4f\n' % origin)
fd.write('%endblock AtomicCoordinatesOrigin\n')
fd.write('\n')
def _write_atomic_coordinates_zmatrix(self, fd, atoms):
"""Write atomic coordinates in Z-matrix format.
Parameters:
- f: An open file object.
- atoms: An atoms object.
"""
species, species_numbers = self.species(atoms)
fd.write('\n')
fd.write('ZM.UnitsLength Ang\n')
fd.write('%block Zmatrix\n')
fd.write(' cartesian\n')
fstr = "{:5d}" + "{:20.10f}" * 3 + "{:3d}" * 3 + "{:7d} {:s}\n"
a2constr = self.make_xyz_constraints(atoms)
a2p, a2s = atoms.get_positions(), atoms.get_chemical_symbols()
for ia, (sp, xyz, ccc, sym) in enumerate(zip(species_numbers,
a2p,
a2constr,
a2s)):
fd.write(fstr.format(
sp, xyz[0], xyz[1], xyz[2], ccc[0],
ccc[1], ccc[2], ia + 1, sym))
fd.write('%endblock Zmatrix\n')
origin = tuple(-atoms.get_celldisp().flatten())
if any(origin):
fd.write('%block AtomicCoordinatesOrigin\n')
fd.write(' %.4f %.4f %.4f\n' % origin)
fd.write('%endblock AtomicCoordinatesOrigin\n')
fd.write('\n')
def make_xyz_constraints(self, atoms):
""" Create coordinate-resolved list of constraints [natoms, 0:3]
The elements of the list must be integers 0 or 1
1 -- means that the coordinate will be updated during relaxation
0 -- mains that the coordinate will be fixed during relaxation
"""
from ase.constraints import FixAtoms, FixedLine, FixedPlane
import sys
import warnings
a = atoms
a2c = np.ones((len(a), 3), dtype=int)
for c in a.constraints:
if isinstance(c, FixAtoms):
a2c[c.get_indices()] = 0
elif isinstance(c, FixedLine):
norm_dir = c.dir / np.linalg.norm(c.dir)
if (max(norm_dir) - 1.0) > 1e-6:
raise RuntimeError(
'norm_dir: {} -- must be one of the Cartesian axes...'
.format(norm_dir))
a2c[c.a] = norm_dir.round().astype(int)
elif isinstance(c, FixedPlane):
norm_dir = c.dir / np.linalg.norm(c.dir)
if (max(norm_dir) - 1.0) > 1e-6:
raise RuntimeError(
'norm_dir: {} -- must be one of the Cartesian axes...'
.format(norm_dir))
a2c[c.a] = abs(1 - norm_dir.round().astype(int))
else:
warnings.warn('Constraint {} is ignored at {}'
.format(str(c), sys._getframe().f_code))
return a2c
def _write_kpts(self, fd):
"""Write kpts.
Parameters:
- f : Open filename.
"""
if self["kpts"] is None:
return
kpts = np.array(self['kpts'])
fd.write('\n')
fd.write('#KPoint grid\n')
fd.write('%block kgrid_Monkhorst_Pack\n')
for i in range(3):
s = ''
if i < len(kpts):
number = kpts[i]
displace = 0.0
else:
number = 1
displace = 0
for j in range(3):
if j == i:
write_this = number
else:
write_this = 0
s += ' %d ' % write_this
s += '%1.1f\n' % displace
fd.write(s)
fd.write('%endblock kgrid_Monkhorst_Pack\n')
fd.write('\n')
def _write_species(self, fd, atoms):
"""Write input related the different species.
Parameters:
- f: An open file object.
- atoms: An atoms object.
"""
species, species_numbers = self.species(atoms)
if self['pseudo_path'] is not None:
pseudo_path = self['pseudo_path']
elif 'SIESTA_PP_PATH' in os.environ:
pseudo_path = os.environ['SIESTA_PP_PATH']
else:
mess = "Please set the environment variable 'SIESTA_PP_PATH'"
raise Exception(mess)
fd.write(format_fdf('NumberOfSpecies', len(species)))
fd.write(format_fdf('NumberOfAtoms', len(atoms)))
pao_basis = []
chemical_labels = []
basis_sizes = []
synth_blocks = []
for species_number, spec in enumerate(species):
species_number += 1
symbol = spec['symbol']
atomic_number = atomic_numbers[symbol]
if spec['pseudopotential'] is None:
if self.pseudo_qualifier() == '':
label = symbol
pseudopotential = label + '.psf'
else:
label = '.'.join([symbol, self.pseudo_qualifier()])
pseudopotential = label + '.psf'
else:
pseudopotential = spec['pseudopotential']
label = os.path.basename(pseudopotential)
label = '.'.join(label.split('.')[:-1])
if not os.path.isabs(pseudopotential):
pseudopotential = join(pseudo_path, pseudopotential)
if not os.path.exists(pseudopotential):
mess = "Pseudopotential '%s' not found" % pseudopotential
raise RuntimeError(mess)
name = os.path.basename(pseudopotential)
name = name.split('.')
name.insert(-1, str(species_number))
if spec['ghost']:
name.insert(-1, 'ghost')
atomic_number = -atomic_number
name = '.'.join(name)
pseudo_targetpath = self.getpath(name)
if join(os.getcwd(), name) != pseudopotential:
if islink(pseudo_targetpath) or isfile(pseudo_targetpath):
os.remove(pseudo_targetpath)
symlink_pseudos = self['symlink_pseudos']
if symlink_pseudos is None:
symlink_pseudos = not os.name == 'nt'
if symlink_pseudos:
os.symlink(pseudopotential, pseudo_targetpath)
else:
shutil.copy(pseudopotential, pseudo_targetpath)
if not spec['excess_charge'] is None:
atomic_number += 200
n_atoms = sum(np.array(species_numbers) == species_number)
paec = float(spec['excess_charge']) / n_atoms
vc = get_valence_charge(pseudopotential)
fraction = float(vc + paec) / vc
pseudo_head = name[:-4]
fractional_command = os.environ['SIESTA_UTIL_FRACTIONAL']
cmd = '%s %s %.7f' % (fractional_command,
pseudo_head,
fraction)
os.system(cmd)
pseudo_head += '-Fraction-%.5f' % fraction
synth_pseudo = pseudo_head + '.psf'
synth_block_filename = pseudo_head + '.synth'
os.remove(name)
shutil.copyfile(synth_pseudo, name)
synth_block = read_vca_synth_block(
synth_block_filename,
species_number=species_number)
synth_blocks.append(synth_block)
if len(synth_blocks) > 0:
fd.write(format_fdf('SyntheticAtoms', list(synth_blocks)))
label = '.'.join(np.array(name.split('.'))[:-1])
string = ' %d %d %s' % (species_number, atomic_number, label)
chemical_labels.append(string)
if isinstance(spec['basis_set'], PAOBasisBlock):
pao_basis.append(spec['basis_set'].script(label))
else:
basis_sizes.append((" " + label, spec['basis_set']))
fd.write((format_fdf('ChemicalSpecieslabel', chemical_labels)))
fd.write('\n')
fd.write((format_fdf('PAO.Basis', pao_basis)))
fd.write((format_fdf('PAO.BasisSizes', basis_sizes)))
fd.write('\n')
def pseudo_qualifier(self):
"""Get the extra string used in the middle of the pseudopotential.
The retrieved pseudopotential for a specific element will be
'H.xxx.psf' for the element 'H' with qualifier 'xxx'. If qualifier
is set to None then the qualifier is set to functional name.
"""
if self['pseudo_qualifier'] is None:
return self['xc'][0].lower()
else:
return self['pseudo_qualifier']
def read_results(self):
"""Read the results.
"""
self.read_number_of_grid_points()
self.read_energy()
self.read_forces_stress()
self.read_eigenvalues()
self.read_kpoints()
self.read_dipole()
self.read_pseudo_density()
self.read_hsx()
self.read_dim()
if self.results['hsx'] is not None:
self.read_pld(self.results['hsx'].norbitals,
len(self.atoms))
self.atoms.cell = self.results['pld'].cell * Bohr
else:
self.results['pld'] = None
self.read_wfsx()
self.read_ion(self.atoms)
self.read_bands()
def read_bands(self):
bandpath = self['bandpath']
if bandpath is None:
return
if len(bandpath.kpts) < 1:
return
fname = self.getpath(ext='bands')
with open(fname) as fd:
kpts, energies, efermi = read_bands_file(fd)
bs = resolve_band_structure(bandpath, kpts, energies, efermi)
self.results['bandstructure'] = bs
def band_structure(self):
return self.results['bandstructure']
def read_ion(self, atoms):
"""
Read the ion.xml file of each specie
"""
from ase.calculators.siesta.import_ion_xml import get_ion
species, species_numbers = self.species(atoms)
self.results['ion'] = {}
for species_number, spec in enumerate(species):
species_number += 1
symbol = spec['symbol']
atomic_number = atomic_numbers[symbol]
if spec['pseudopotential'] is None:
if self.pseudo_qualifier() == '':
label = symbol
else:
label = '.'.join([symbol, self.pseudo_qualifier()])
pseudopotential = self.getpath(label, 'psf')
else:
pseudopotential = spec['pseudopotential']
label = os.path.basename(pseudopotential)
label = '.'.join(label.split('.')[:-1])
name = os.path.basename(pseudopotential)
name = name.split('.')
name.insert(-1, str(species_number))
if spec['ghost']:
name.insert(-1, 'ghost')
atomic_number = -atomic_number
name = '.'.join(name)
label = '.'.join(np.array(name.split('.'))[:-1])
if label not in self.results['ion']:
fname = self.getpath(label, 'ion.xml')
if os.path.isfile(fname):
self.results['ion'][label] = get_ion(fname)
def read_hsx(self):
"""
Read the siesta HSX file.
return a namedtuple with the following arguments:
'norbitals', 'norbitals_sc', 'nspin', 'nonzero',
'is_gamma', 'sc_orb2uc_orb', 'row2nnzero', 'sparse_ind2column',
'H_sparse', 'S_sparse', 'aB2RaB_sparse', 'total_elec_charge', 'temp'
"""
from ase.calculators.siesta.import_functions import readHSX
filename = self.getpath(ext='HSX')
if isfile(filename):
self.results['hsx'] = readHSX(filename)
else:
self.results['hsx'] = None
def read_dim(self):
"""
Read the siesta DIM file
Retrun a namedtuple with the following arguments:
'natoms_sc', 'norbitals_sc', 'norbitals', 'nspin',
'nnonzero', 'natoms_interacting'
"""
from ase.calculators.siesta.import_functions import readDIM
filename = self.getpath(ext='DIM')
if isfile(filename):
self.results['dim'] = readDIM(filename)
else:
self.results['dim'] = None
def read_pld(self, norb, natms):
"""
Read the siesta PLD file
Return a namedtuple with the following arguments:
'max_rcut', 'orb2ao', 'orb2uorb', 'orb2occ', 'atm2sp',
'atm2shift', 'coord_sc', 'cell', 'nunit_cells'
"""
from ase.calculators.siesta.import_functions import readPLD
filename = self.getpath(ext='PLD')
if isfile(filename):
self.results['pld'] = readPLD(filename, norb, natms)
else:
self.results['pld'] = None
def read_wfsx(self):
"""
Read the siesta WFSX file
Return a namedtuple with the following arguments:
"""
from ase.calculators.siesta.import_functions import readWFSX
fname_woext = os.path.join(self.directory, self.prefix)
if isfile(fname_woext + '.WFSX'):
filename = fname_woext + '.WFSX'
self.results['wfsx'] = readWFSX(filename)
elif isfile(fname_woext + '.fullBZ.WFSX'):
filename = fname_woext + '.fullBZ.WFSX'
readWFSX(filename)
self.results['wfsx'] = readWFSX(filename)
else:
self.results['wfsx'] = None
def read_pseudo_density(self):
"""Read the density if it is there."""
filename = self.getpath(ext='RHO')
if isfile(filename):
self.results['density'] = read_rho(filename)
def read_number_of_grid_points(self):
"""Read number of grid points from SIESTA's text-output file. """
fname = self.getpath(ext='out')
with open(fname, 'r') as fd:
for line in fd:
line = line.strip().lower()
if line.startswith('initmesh: mesh ='):
n_points = [int(word) for word in line.split()[3:8:2]]
self.results['n_grid_point'] = n_points
break
else:
raise RuntimeError
def read_energy(self):
"""Read energy from SIESTA's text-output file.
"""
fname = self.getpath(ext='out')
with open(fname, 'r') as fd:
text = fd.read().lower()
assert 'final energy' in text
lines = iter(text.split('\n'))
# Get the energy and free energy the last time it appears
for line in lines:
has_energy = line.startswith('siesta: etot =')
if has_energy:
self.results['energy'] = float(line.split()[-1])
line = next(lines)
self.results['free_energy'] = float(line.split()[-1])
if ('energy' not in self.results or
'free_energy' not in self.results):
raise RuntimeError
def read_forces_stress(self):
"""Read the forces and stress from the FORCE_STRESS file.
"""
fname = self.getpath('FORCE_STRESS')
with open(fname, 'r') as fd:
lines = fd.readlines()
stress_lines = lines[1:4]
stress = np.empty((3, 3))
for i in range(3):
line = stress_lines[i].strip().split(' ')
line = [s for s in line if len(s) > 0]
stress[i] = [float(s) for s in line]
self.results['stress'] = np.array(
[stress[0, 0], stress[1, 1], stress[2, 2],
stress[1, 2], stress[0, 2], stress[0, 1]])
self.results['stress'] *= Ry / Bohr**3
start = 5
self.results['forces'] = np.zeros((len(lines) - start, 3), float)
for i in range(start, len(lines)):
line = [s for s in lines[i].strip().split(' ') if len(s) > 0]
self.results['forces'][i - start] = [float(s) for s in line[2:5]]
self.results['forces'] *= Ry / Bohr
def read_eigenvalues(self):
""" A robust procedure using the suggestion by Federico Marchesin """
fname = self.getpath(ext='EIG')
try:
with open(fname, "r") as fd:
self.results['fermi_energy'] = float(fd.readline())
n, nspin, nkp = map(int, fd.readline().split())
_ee = np.split(
np.array(fd.read().split()).astype(float), nkp)
except (IOError):
return 1
ksn2e = np.delete(_ee, 0, 1).reshape([nkp, nspin, n])
eigarray = np.empty((nspin, nkp, n))
eigarray[:] = np.inf
for k, sn2e in enumerate(ksn2e):
for s, n2e in enumerate(sn2e):
eigarray[s, k, :] = n2e
assert np.isfinite(eigarray).all()
self.results['eigenvalues'] = eigarray
return 0
def read_kpoints(self):
""" Reader of the .KP files """
fname = self.getpath(ext='KP')
try:
with open(fname, "r") as fd:
nkp = int(next(fd))
kpoints = np.empty((nkp, 3))
kweights = np.empty(nkp)
for i in range(nkp):
line = next(fd)
tokens = line.split()
numbers = np.array(tokens[1:]).astype(float)
kpoints[i] = numbers[:3]
kweights[i] = numbers[3]
except (IOError):
return 1
self.results['kpoints'] = kpoints
self.results['kweights'] = kweights
return 0
def read_dipole(self):
"""Read dipole moment. """
dipole = np.zeros([1, 3])
with open(self.getpath(ext='out'), 'r') as fd:
for line in fd:
if line.rfind('Electric dipole (Debye)') > -1:
dipole = np.array([float(f) for f in line.split()[5:8]])
# debye to e*Ang
self.results['dipole'] = dipole * 0.2081943482534
def get_fermi_level(self):
return self.results['fermi_energy']
def get_k_point_weights(self):
return self.results['kweights']
def get_ibz_k_points(self):
return self.results['kpoints']
class Siesta3_2(Siesta):
def __init__(self, *args, **kwargs):
warnings.warn(
"The Siesta3_2 calculator class will no longer be supported. "
"Use 'ase.calculators.siesta.Siesta in stead. "
"If using the ASE interface with SIESTA 3.2 you must explicitly "
"include the keywords 'SpinPolarized', 'NonCollinearSpin' and "
"'SpinOrbit' if needed.",
np.VisibleDeprecationWarning)
Siesta.__init__(self, *args, **kwargs)