""" This module defines a FileIOCalculator for DFTB+
http://www.dftbplus.org/
http://www.dftb.org/
Initial development: markus.kaukonen@iki.fi
"""
import os
import numpy as np
from ase.calculators.calculator import (FileIOCalculator, kpts2ndarray,
kpts2sizeandoffsets)
from ase.units import Hartree, Bohr
[docs]class Dftb(FileIOCalculator):
if 'DFTB_COMMAND' in os.environ:
command = os.environ['DFTB_COMMAND'] + ' > PREFIX.out'
else:
command = 'dftb+ > PREFIX.out'
implemented_properties = ['energy', 'forces', 'charges', 'stress']
def __init__(self, restart=None, ignore_bad_restart_file=False,
label='dftb', atoms=None, kpts=None,
slako_dir=None,
**kwargs):
"""
All keywords for the dftb_in.hsd input file (see the DFTB+ manual)
can be set by ASE. Consider the following input file block:
>>> Hamiltonian = DFTB {
>>> SCC = Yes
>>> SCCTolerance = 1e-8
>>> MaxAngularMomentum = {
>>> H = s
>>> O = p
>>> }
>>> }
This can be generated by the DFTB+ calculator by using the
following settings:
>>> calc = Dftb(Hamiltonian_='DFTB', # line is included by default
>>> Hamiltonian_SCC='Yes',
>>> Hamiltonian_SCCTolerance=1e-8,
>>> Hamiltonian_MaxAngularMomentum_='',
>>> Hamiltonian_MaxAngularMomentum_H='s',
>>> Hamiltonian_MaxAngularMomentum_O='p')
In addition to keywords specific to DFTB+, also the following keywords
arguments can be used:
restart: str
Prefix for restart file. May contain a directory.
Default is None: don't restart.
ignore_bad_restart_file: bool
Ignore broken or missing restart file. By default, it is an
error if the restart file is missing or broken.
label: str (default 'dftb')
Prefix used for the main output file (<label>.out).
atoms: Atoms object (default None)
Optional Atoms object to which the calculator will be
attached. When restarting, atoms will get its positions and
unit-cell updated from file.
kpts: (default None)
Brillouin zone sampling:
* ``(1,1,1)`` or ``None``: Gamma-point only
* ``(n1,n2,n3)``: Monkhorst-Pack grid
* ``dict``: Interpreted as a path in the Brillouin zone if
it contains the 'path_' keyword. Otherwise it is converted
into a Monkhorst-Pack grid using
``ase.calculators.calculator.kpts2sizeandoffsets``
* ``[(k11,k12,k13),(k21,k22,k23),...]``: Explicit (Nkpts x 3)
array of k-points in units of the reciprocal lattice vectors
(each with equal weight)
Additional attribute to be set by the embed() method:
pcpot: PointCharge object
An external point charge potential (for QM/MM calculations)
"""
if slako_dir is None:
slako_dir = os.environ.get('DFTB_PREFIX', './')
if not slako_dir.endswith('/'):
slako_dir += '/'
self.slako_dir = slako_dir
self.default_parameters = dict(
Hamiltonian_='DFTB',
Hamiltonian_SlaterKosterFiles_='Type2FileNames',
Hamiltonian_SlaterKosterFiles_Prefix=self.slako_dir,
Hamiltonian_SlaterKosterFiles_Separator='"-"',
Hamiltonian_SlaterKosterFiles_Suffix='".skf"',
Hamiltonian_MaxAngularMomentum_='',
Options_='',
Options_WriteResultsTag='Yes')
self.pcpot = None
self.lines = None
self.atoms = None
self.atoms_input = None
self.do_forces = False
self.outfilename = 'dftb.out'
FileIOCalculator.__init__(self, restart, ignore_bad_restart_file,
label, atoms,
**kwargs)
# Determine number of spin channels
try:
entry = kwargs['Hamiltonian_SpinPolarisation']
spinpol = 'colinear' in entry.lower()
except KeyError:
spinpol = False
self.nspin = 2 if spinpol else 1
# kpoint stuff by ase
self.kpts = kpts
self.kpts_coord = None
if self.kpts is not None:
initkey = 'Hamiltonian_KPointsAndWeights'
mp_mesh = None
offsets = None
if isinstance(self.kpts, dict):
if 'path' in self.kpts:
# kpts is path in Brillouin zone
self.parameters[initkey + '_'] = 'Klines '
self.kpts_coord = kpts2ndarray(self.kpts, atoms=atoms)
else:
# kpts is (implicit) definition of
# Monkhorst-Pack grid
self.parameters[initkey + '_'] = 'SupercellFolding '
mp_mesh, offsets = kpts2sizeandoffsets(atoms=atoms,
**self.kpts)
elif np.array(self.kpts).ndim == 1:
# kpts is Monkhorst-Pack grid
self.parameters[initkey + '_'] = 'SupercellFolding '
mp_mesh = self.kpts
offsets = [0.] * 3
elif np.array(self.kpts).ndim == 2:
# kpts is (N x 3) list/array of k-point coordinates
# each will be given equal weight
self.parameters[initkey + '_'] = ''
self.kpts_coord = np.array(self.kpts)
else:
raise ValueError('Illegal kpts definition:' + str(self.kpts))
if mp_mesh is not None:
eps = 1e-10
for i in range(3):
key = initkey + '_empty%03d' % i
val = [mp_mesh[i] if j == i else 0 for j in range(3)]
self.parameters[key] = ' '.join(map(str, val))
offsets[i] *= mp_mesh[i]
assert abs(offsets[i]) < eps or abs(offsets[i] - 0.5) < eps
# DFTB+ uses a different offset convention, where
# the k-point mesh is already Gamma-centered prior
# to the addition of any offsets
if mp_mesh[i] % 2 == 0:
offsets[i] += 0.5
key = initkey + '_empty%03d' % 3
self.parameters[key] = ' '.join(map(str, offsets))
elif self.kpts_coord is not None:
for i, c in enumerate(self.kpts_coord):
key = initkey + '_empty%09d' % i
c_str = ' '.join(map(str, c))
if 'Klines' in self.parameters[initkey + '_']:
c_str = '1 ' + c_str
else:
c_str += ' 1.0'
self.parameters[key] = c_str
def write_dftb_in(self, filename):
""" Write the innput file for the dftb+ calculation.
Geometry is taken always from the file 'geo_end.gen'.
"""
outfile = open(filename, 'w')
outfile.write('Geometry = GenFormat { \n')
outfile.write(' <<< "geo_end.gen" \n')
outfile.write('} \n')
outfile.write(' \n')
params = self.parameters.copy()
s = 'Hamiltonian_MaxAngularMomentum_'
for key in params:
if key.startswith(s) and len(key) > len(s):
break
else:
# User didn't specify max angular mometa. Get them from
# the .skf files:
symbols = set(self.atoms.get_chemical_symbols())
for symbol in symbols:
path = os.path.join(self.slako_dir,
'{0}-{0}.skf'.format(symbol))
l = read_max_angular_momentum(path)
params[s + symbol] = '"{}"'.format('spdf'[l])
# --------MAIN KEYWORDS-------
previous_key = 'dummy_'
myspace = ' '
for key, value in sorted(params.items()):
current_depth = key.rstrip('_').count('_')
previous_depth = previous_key.rstrip('_').count('_')
for my_backsclash in reversed(
range(previous_depth - current_depth)):
outfile.write(3 * (1 + my_backsclash) * myspace + '} \n')
outfile.write(3 * current_depth * myspace)
if key.endswith('_') and len(value) > 0:
outfile.write(key.rstrip('_').rsplit('_')[-1] +
' = ' + str(value) + '{ \n')
elif (key.endswith('_') and (len(value) == 0)
and current_depth == 0): # E.g. 'Options {'
outfile.write(key.rstrip('_').rsplit('_')[-1] +
' ' + str(value) + '{ \n')
elif (key.endswith('_') and (len(value) == 0)
and current_depth > 0): # E.g. 'Hamiltonian_Max... = {'
outfile.write(key.rstrip('_').rsplit('_')[-1] +
' = ' + str(value) + '{ \n')
elif key.count('_empty') == 1:
outfile.write(str(value) + ' \n')
elif ((key == 'Hamiltonian_ReadInitialCharges') and
(str(value).upper() == 'YES')):
f1 = os.path.isfile(self.directory + os.sep + 'charges.dat')
f2 = os.path.isfile(self.directory + os.sep + 'charges.bin')
if not (f1 or f2):
print('charges.dat or .bin not found, switching off guess')
value = 'No'
outfile.write(key.rsplit('_')[-1] + ' = ' + str(value) + ' \n')
else:
outfile.write(key.rsplit('_')[-1] + ' = ' + str(value) + ' \n')
if self.pcpot is not None and ('DFTB' in str(value)):
outfile.write(' ElectricField = { \n')
outfile.write(' PointCharges = { \n')
outfile.write(
' CoordsAndCharges [Angstrom] = DirectRead { \n')
outfile.write(' Records = ' +
str(len(self.pcpot.mmcharges)) + ' \n')
outfile.write(
' File = "dftb_external_charges.dat" \n')
outfile.write(' } \n')
outfile.write(' } \n')
outfile.write(' } \n')
previous_key = key
current_depth = key.rstrip('_').count('_')
for my_backsclash in reversed(range(current_depth)):
outfile.write(3 * my_backsclash * myspace + '} \n')
outfile.write('ParserOptions { \n')
outfile.write(' IgnoreUnprocessedNodes = Yes \n')
outfile.write('} \n')
if self.do_forces:
outfile.write('Analysis { \n')
outfile.write(' CalculateForces = Yes \n')
outfile.write('} \n')
outfile.close()
def set(self, **kwargs):
changed_parameters = FileIOCalculator.set(self, **kwargs)
if changed_parameters:
self.reset()
return changed_parameters
def check_state(self, atoms):
system_changes = FileIOCalculator.check_state(self, atoms)
# Ignore unit cell for molecules:
if not atoms.pbc.any() and 'cell' in system_changes:
system_changes.remove('cell')
if self.pcpot and self.pcpot.mmpositions is not None:
system_changes.append('positions')
return system_changes
def write_input(self, atoms, properties=None, system_changes=None):
from ase.io import write
if properties is not None:
if 'forces' in properties or 'stress' in properties:
self.do_forces = True
FileIOCalculator.write_input(
self, atoms, properties, system_changes)
self.write_dftb_in(os.path.join(self.directory, 'dftb_in.hsd'))
write(os.path.join(self.directory, 'geo_end.gen'), atoms)
# self.atoms is none until results are read out,
# then it is set to the ones at writing input
self.atoms_input = atoms
self.atoms = None
if self.pcpot:
self.pcpot.write_mmcharges('dftb_external_charges.dat')
def read_results(self):
""" all results are read from results.tag file
It will be destroyed after it is read to avoid
reading it once again after some runtime error """
myfile = open(os.path.join(self.directory, 'results.tag'), 'r')
self.lines = myfile.readlines()
myfile.close()
self.atoms = self.atoms_input
charges, energy = self.read_charges_and_energy()
if charges is not None:
self.results['charges'] = charges
self.results['energy'] = energy
if self.do_forces:
forces = self.read_forces()
self.results['forces'] = forces
self.mmpositions = None
# stress stuff begins
sstring = 'stress'
have_stress = False
stress = list()
for iline, line in enumerate(self.lines):
if sstring in line:
have_stress = True
start = iline + 1
end = start + 3
for i in range(start, end):
cell = [float(x) for x in self.lines[i].split()]
stress.append(cell)
if have_stress:
stress = -np.array(stress) * Hartree / Bohr**3
self.results['stress'] = stress.flat[[0, 4, 8, 5, 2, 1]]
# stress stuff ends
# eigenvalues and fermi levels
fermi_levels = self.read_fermi_levels()
if fermi_levels is not None:
self.results['fermi_levels'] = fermi_levels
eigenvalues = self.read_eigenvalues()
if eigenvalues is not None:
self.results['eigenvalues'] = eigenvalues
# calculation was carried out with atoms written in write_input
os.remove(os.path.join(self.directory, 'results.tag'))
def read_forces(self):
"""Read Forces from dftb output file (results.tag)."""
from ase.units import Hartree, Bohr
# Initialise the indices so their scope
# reaches outside of the for loop
index_force_begin = -1
index_force_end = -1
# Force line indexes
for iline, line in enumerate(self.lines):
fstring = 'forces '
if line.find(fstring) >= 0:
index_force_begin = iline + 1
line1 = line.replace(':', ',')
index_force_end = iline + 1 + \
int(line1.split(',')[-1])
break
gradients = []
for j in range(index_force_begin, index_force_end):
word = self.lines[j].split()
gradients.append([float(word[k]) for k in range(0, 3)])
return np.array(gradients) * Hartree / Bohr
def read_charges_and_energy(self):
"""Get partial charges on atoms
in case we cannot find charges they are set to None
"""
infile = open(os.path.join(self.directory, 'detailed.out'), 'r')
lines = infile.readlines()
infile.close()
for line in lines:
if line.strip().startswith('Total energy:'):
energy = float(line.split()[2]) * Hartree
break
qm_charges = []
for n, line in enumerate(lines):
if ('Atom' and 'Charge' in line):
chargestart = n + 1
break
else:
# print('Warning: did not find DFTB-charges')
# print('This is ok if flag SCC=No')
return None, energy
lines1 = lines[chargestart:(chargestart + len(self.atoms))]
for line in lines1:
qm_charges.append(float(line.split()[-1]))
return np.array(qm_charges), energy
def get_charges(self, atoms):
""" Get the calculated charges
this is inhereted to atoms object """
if 'charges' in self.results:
return self.results['charges']
else:
return None
def read_eigenvalues(self):
""" Read Eigenvalues from dftb output file (results.tag).
Unfortunately, the order seems to be scrambled. """
# Eigenvalue line indexes
index_eig_begin = None
for iline, line in enumerate(self.lines):
fstring = 'eigenvalues '
if line.find(fstring) >= 0:
index_eig_begin = iline + 1
line1 = line.replace(':', ',')
ncol, nband, nkpt, nspin = map(int, line1.split(',')[-4:])
break
else:
return None
# Take into account that the last row may lack
# columns if nkpt * nspin * nband % ncol != 0
nrow = int(np.ceil(nkpt * nspin * nband * 1. / ncol))
index_eig_end = index_eig_begin + nrow
ncol_last = len(self.lines[index_eig_end - 1].split())
self.lines[index_eig_end - 1] += ' 0.0 ' * (ncol - ncol_last)
eig = np.loadtxt(self.lines[index_eig_begin:index_eig_end]).flatten()
eig *= Hartree
N = nkpt * nband
eigenvalues = [eig[i * N:(i + 1) * N].reshape((nkpt, nband))
for i in range(nspin)]
return eigenvalues
def read_fermi_levels(self):
""" Read Fermi level(s) from dftb output file (results.tag). """
# Fermi level line indexes
for iline, line in enumerate(self.lines):
fstring = 'fermi_level '
if line.find(fstring) >= 0:
index_fermi = iline + 1
break
else:
return None
fermi_levels = []
words = self.lines[index_fermi].split()
assert len(words) == 2
for word in words:
e = float(word)
if abs(e) > 1e-8:
# Without spin polarization, one of the Fermi
# levels is equal to 0.000000000000000E+000
fermi_levels.append(e)
return np.array(fermi_levels) * Hartree
def get_ibz_k_points(self):
return self.kpts_coord.copy()
def get_number_of_spins(self):
return self.nspin
def get_eigenvalues(self, kpt=0, spin=0):
return self.results['eigenvalues'][spin][kpt].copy()
def get_fermi_levels(self):
return self.results['fermi_levels'].copy()
def get_fermi_level(self):
return max(self.get_fermi_levels())
def embed(self, mmcharges=None, directory='./'):
"""Embed atoms in point-charges (mmcharges)
"""
self.pcpot = PointChargePotential(mmcharges, self.directory)
return self.pcpot
class PointChargePotential:
def __init__(self, mmcharges, directory='./'):
"""Point-charge potential for DFTB+.
"""
self.mmcharges = mmcharges
self.directory = directory
self.mmpositions = None
self.mmforces = None
def set_positions(self, mmpositions):
self.mmpositions = mmpositions
def set_charges(self, mmcharges):
self.mmcharges = mmcharges
def write_mmcharges(self, filename='dftb_external_charges.dat'):
""" mok all
write external charges as monopoles for dftb+.
"""
if self.mmcharges is None:
print("DFTB: Warning: not writing exernal charges ")
return
charge_file = open(os.path.join(self.directory, filename), 'w')
for [pos, charge] in zip(self.mmpositions, self.mmcharges):
[x, y, z] = pos
charge_file.write('%12.6f %12.6f %12.6f %12.6f \n'
% (x, y, z, charge))
charge_file.close()
def get_forces(self, calc, get_forces=True):
""" returns forces on point charges if the flag get_forces=True """
if get_forces:
return self.read_forces_on_pointcharges()
else:
return np.zeros_like(self.mmpositions)
def read_forces_on_pointcharges(self):
"""Read Forces from dftb output file (results.tag)."""
from ase.units import Hartree, Bohr
infile = open(os.path.join(self.directory, 'detailed.out'), 'r')
lines = infile.readlines()
infile.close()
external_forces = []
for n, line in enumerate(lines):
if ('Forces on external charges' in line):
chargestart = n + 1
break
else:
raise RuntimeError(
'Problem in reading forces on MM external-charges')
lines1 = lines[chargestart:(chargestart + len(self.mmcharges))]
for line in lines1:
external_forces.append(
[float(i) for i in line.split()])
return np.array(external_forces) * Hartree / Bohr
def read_max_angular_momentum(path):
"""Read maximum angular momentum from .skf file.
See dftb.org for A detailed description of the Slater-Koster file format.
"""
with open(path, 'r') as fd:
line = fd.readline()
if line[0] == '@':
# Extended format
fd.readline()
l = 3
pos = 9
else:
# Simple format:
l = 2
pos = 7
# Sometimes there ar commas, sometimes not:
line = fd.readline().replace(',', ' ')
occs = [float(f) for f in line.split()[pos:pos + l + 1]]
for f in occs:
if f > 0.0:
return l
l -= 1