Coverage for ase / calculators / dftb.py: 75.74%
338 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-03-30 08:22 +0000
1# fmt: off
3""" This module defines a FileIOCalculator for DFTB+
5https://www.dftbplus.org/
6https://www.dftb.org/
8Initial development: markus.kaukonen@iki.fi
9"""
11import os
13import numpy as np
15from ase.calculators.calculator import (
16 BadConfiguration,
17 FileIOCalculator,
18 kpts2ndarray,
19 kpts2sizeandoffsets,
20)
21from ase.units import Bohr, Hartree
24class Dftb(FileIOCalculator):
25 implemented_properties = ['energy', 'forces', 'charges',
26 'stress', 'dipole']
27 discard_results_on_any_change = True
29 fileio_rules = FileIOCalculator.ruleset(
30 configspec=dict(skt_path=None),
31 stdout_name='{prefix}.out')
33 def __init__(self, restart=None,
34 ignore_bad_restart_file=FileIOCalculator._deprecated,
35 label='dftb', atoms=None, kpts=None,
36 slako_dir=None,
37 command=None,
38 profile=None,
39 with_tags=False,
40 **kwargs):
41 """
42 All keywords for the dftb_in.hsd input file (see the DFTB+ manual)
43 can be set by ASE. Consider the following input file block::
45 Hamiltonian = DFTB {
46 SCC = Yes
47 SCCTolerance = 1e-8
48 MaxAngularMomentum = {
49 H = s
50 O = p
51 }
52 }
54 This can be generated by the DFTB+ calculator by using the
55 following settings:
57 >>> from ase.calculators.dftb import Dftb
58 >>>
59 >>> calc = Dftb(Hamiltonian_='DFTB', # line is included by default
60 ... Hamiltonian_SCC='Yes',
61 ... Hamiltonian_SCCTolerance=1e-8,
62 ... Hamiltonian_MaxAngularMomentum_='',
63 ... Hamiltonian_MaxAngularMomentum_H='s',
64 ... Hamiltonian_MaxAngularMomentum_O='p')
66 In addition to keywords specific to DFTB+, also the following keywords
67 arguments can be used:
69 restart: str
70 Prefix for restart file. May contain a directory.
71 Default is None: don't restart.
72 ignore_bad_restart_file: bool
73 Ignore broken or missing restart file. By default, it is an
74 error if the restart file is missing or broken.
75 label: str (default 'dftb')
76 Prefix used for the main output file (<label>.out).
77 atoms: Atoms object (default None)
78 Optional Atoms object to which the calculator will be
79 attached. When restarting, atoms will get its positions and
80 unit-cell updated from file.
81 kpts: (default None)
82 Brillouin zone sampling:
84 * ``(1,1,1)`` or ``None``: Gamma-point only
85 * ``(n1,n2,n3)``: Monkhorst-Pack grid
86 * ``dict``: Interpreted as a path in the Brillouin zone if
87 it contains the 'path_' keyword. Otherwise it is converted
88 into a Monkhorst-Pack grid using
89 ``ase.calculators.calculator.kpts2sizeandoffsets``
90 * ``[(k11,k12,k13),(k21,k22,k23),...]``: Explicit (Nkpts x 3)
91 array of k-points in units of the reciprocal lattice vectors
92 (each with equal weight)
93 with_tags: bool (default False)
94 Whether to use atom tags to distinguish between different chemical
95 species. When enabled, atoms with different tags are treated as
96 distinct species, even if they have the same chemical symbol.
97 For example, hydrogen atoms in a nanoparticle (tag 1) and hydrogen
98 atoms in a solvent (tag 2) would be represented as "H1" and "H2",
99 requiring separate Slater-Koster files: "H1-H1.skf", "H1-H2.skf",
100 and "H2-H2.skf".
102 Additional attribute to be set by the embed() method:
104 pcpot: PointCharge object
105 An external point charge potential (for QM/MM calculations)
106 """
108 if command is None:
109 if 'DFTB_COMMAND' in self.cfg:
110 command = self.cfg['DFTB_COMMAND'] + ' > PREFIX.out'
112 if command is None and profile is None:
113 try:
114 profile = self.load_argv_profile(self.cfg, 'dftb')
115 except BadConfiguration:
116 pass
118 if command is None:
119 command = 'dftb+ > PREFIX.out'
121 if slako_dir is None:
122 if profile is not None:
123 slako_dir = profile.configvars.get('skt_path')
125 if slako_dir is None:
126 slako_dir = self.cfg.get('DFTB_PREFIX', './')
127 if not slako_dir.endswith('/'):
128 slako_dir += '/'
130 self.slako_dir = slako_dir
131 if kwargs.get('Hamiltonian_', 'DFTB') == 'DFTB':
132 self.default_parameters = dict(
133 Hamiltonian_='DFTB',
134 Hamiltonian_SlaterKosterFiles_='Type2FileNames',
135 Hamiltonian_SlaterKosterFiles_Prefix=self.slako_dir,
136 Hamiltonian_SlaterKosterFiles_Separator='"-"',
137 Hamiltonian_SlaterKosterFiles_Suffix='".skf"',
138 Hamiltonian_MaxAngularMomentum_='',
139 Options_='',
140 Options_WriteResultsTag='Yes',
141 ParserOptions_='',
142 ParserOptions_ParserVersion=1,
143 ParserOptions_IgnoreUnprocessedNodes='Yes')
144 else:
145 self.default_parameters = dict(
146 Options_='',
147 Options_WriteResultsTag='Yes',
148 ParserOptions_='',
149 ParserOptions_ParserVersion=1,
150 ParserOptions_IgnoreUnprocessedNodes='Yes')
152 self.pcpot = None
153 self.lines = None
154 self.atoms = None
155 self.atoms_input = None
156 self.do_forces = False
157 # Whether to use tags of the atoms for the calculation
158 self.with_tags = with_tags
160 super().__init__(restart, ignore_bad_restart_file,
161 label, atoms, command=command,
162 profile=profile, **kwargs)
164 # Determine number of spin channels
165 try:
166 entry = kwargs['Hamiltonian_SpinPolarisation']
167 spinpol = 'colinear' in entry.lower()
168 except KeyError:
169 spinpol = False
170 self.nspin = 2 if spinpol else 1
172 # kpoint stuff by ase
173 self.kpts = kpts
174 self.kpts_coord = None
176 if self.kpts is not None:
177 initkey = 'Hamiltonian_KPointsAndWeights'
178 mp_mesh = None
179 offsets = None
181 if isinstance(self.kpts, dict):
182 if 'path' in self.kpts:
183 # kpts is path in Brillouin zone
184 self.parameters[initkey + '_'] = 'Klines '
185 self.kpts_coord = kpts2ndarray(self.kpts, atoms=atoms)
186 else:
187 # kpts is (implicit) definition of
188 # Monkhorst-Pack grid
189 self.parameters[initkey + '_'] = 'SupercellFolding '
190 mp_mesh, offsets = kpts2sizeandoffsets(atoms=atoms,
191 **self.kpts)
192 elif np.array(self.kpts).ndim == 1:
193 # kpts is Monkhorst-Pack grid
194 self.parameters[initkey + '_'] = 'SupercellFolding '
195 mp_mesh = self.kpts
196 offsets = [0.] * 3
197 elif np.array(self.kpts).ndim == 2:
198 # kpts is (N x 3) list/array of k-point coordinates
199 # each will be given equal weight
200 self.parameters[initkey + '_'] = ''
201 self.kpts_coord = np.array(self.kpts)
202 else:
203 raise ValueError('Illegal kpts definition:' + str(self.kpts))
205 if mp_mesh is not None:
206 eps = 1e-10
207 for i in range(3):
208 key = initkey + '_empty%03d' % i
209 val = [mp_mesh[i] if j == i else 0 for j in range(3)]
210 self.parameters[key] = ' '.join(map(str, val))
211 offsets[i] *= mp_mesh[i]
212 assert abs(offsets[i]) < eps or abs(offsets[i] - 0.5) < eps
213 # DFTB+ uses a different offset convention, where
214 # the k-point mesh is already Gamma-centered prior
215 # to the addition of any offsets
216 if mp_mesh[i] % 2 == 0:
217 offsets[i] += 0.5
218 key = initkey + '_empty%03d' % 3
219 self.parameters[key] = ' '.join(map(str, offsets))
221 elif self.kpts_coord is not None:
222 for i, c in enumerate(self.kpts_coord):
223 key = initkey + '_empty%09d' % i
224 c_str = ' '.join(map(str, c))
225 if 'Klines' in self.parameters[initkey + '_']:
226 c_str = '1 ' + c_str
227 else:
228 c_str += ' 1.0'
229 self.parameters[key] = c_str
231 def write_dftb_in(self, outfile):
232 """ Write the innput file for the dftb+ calculation.
233 Geometry is taken always from the file 'geo_end.gen'.
234 """
236 outfile.write('Geometry = GenFormat { \n')
237 outfile.write(' <<< "geo_end.gen" \n')
238 outfile.write('} \n')
239 outfile.write(' \n')
241 params = self.parameters.copy()
243 s = 'Hamiltonian_MaxAngularMomentum_'
244 for key in params:
245 if key.startswith(s) and len(key) > len(s):
246 break
247 else:
248 if params.get('Hamiltonian_', 'DFTB') == 'DFTB':
249 # User didn't specify max angular momenta. Get them from
250 # the .skf files:
251 symbols = self.atoms.get_chemical_symbols()
252 if self.with_tags:
253 tags = set(self.atoms.get_tags())
254 symbols = [_s + str(tag) for _s, tag in
255 zip(self.atoms.get_chemical_symbols(), tags)]
256 symbols = set(symbols)
257 for symbol in symbols:
258 path = os.path.join(self.slako_dir,
259 '{0}-{0}.skf'.format(symbol))
260 l = read_max_angular_momentum(path)
261 params[s + symbol] = '"{}"'.format('spdf'[l])
263 if self.do_forces:
264 params['Analysis_'] = ''
265 params['Analysis_CalculateForces'] = 'Yes'
267 # --------MAIN KEYWORDS-------
268 previous_key = 'dummy_'
269 myspace = ' '
270 for key, value in sorted(params.items()):
271 current_depth = key.rstrip('_').count('_')
272 previous_depth = previous_key.rstrip('_').count('_')
273 for my_backslash in reversed(
274 range(previous_depth - current_depth)):
275 outfile.write(3 * (1 + my_backslash) * myspace + '} \n')
276 outfile.write(3 * current_depth * myspace)
277 if key.endswith('_') and len(value) > 0:
278 outfile.write(key.rstrip('_').rsplit('_')[-1] +
279 ' = ' + str(value) + '{ \n')
280 elif (key.endswith('_') and (len(value) == 0)
281 and current_depth == 0): # E.g. 'Options {'
282 outfile.write(key.rstrip('_').rsplit('_')[-1] +
283 ' ' + str(value) + '{ \n')
284 elif (key.endswith('_') and (len(value) == 0)
285 and current_depth > 0): # E.g. 'Hamiltonian_Max... = {'
286 outfile.write(key.rstrip('_').rsplit('_')[-1] +
287 ' = ' + str(value) + '{ \n')
288 elif key.count('_empty') == 1:
289 outfile.write(str(value) + ' \n')
290 elif ((key == 'Hamiltonian_ReadInitialCharges') and
291 (str(value).upper() == 'YES')):
292 f1 = os.path.isfile(self.directory + os.sep + 'charges.dat')
293 f2 = os.path.isfile(self.directory + os.sep + 'charges.bin')
294 if not (f1 or f2):
295 print('charges.dat or .bin not found, switching off guess')
296 value = 'No'
297 outfile.write(key.rsplit('_')[-1] + ' = ' + str(value) + ' \n')
298 else:
299 outfile.write(key.rsplit('_')[-1] + ' = ' + str(value) + ' \n')
300 if self.pcpot is not None and ('DFTB' in str(value)):
301 outfile.write(' ElectricField = { \n')
302 outfile.write(' PointCharges = { \n')
303 outfile.write(
304 ' CoordsAndCharges [Angstrom] = DirectRead { \n')
305 outfile.write(' Records = ' +
306 str(len(self.pcpot.mmcharges)) + ' \n')
307 outfile.write(
308 ' File = "dftb_external_charges.dat" \n')
309 outfile.write(' } \n')
310 outfile.write(' } \n')
311 outfile.write(' } \n')
312 previous_key = key
313 current_depth = key.rstrip('_').count('_')
314 for my_backslash in reversed(range(current_depth)):
315 outfile.write(3 * my_backslash * myspace + '} \n')
317 def check_state(self, atoms):
318 system_changes = FileIOCalculator.check_state(self, atoms)
319 # Ignore unit cell for molecules:
320 if not atoms.pbc.any() and 'cell' in system_changes:
321 system_changes.remove('cell')
322 if self.pcpot and self.pcpot.mmpositions is not None:
323 system_changes.append('positions')
324 return system_changes
326 def write_input(self, atoms, properties=None, system_changes=None):
327 from ase.io.gen import write_gen
328 if properties is not None:
329 if 'forces' in properties or 'stress' in properties:
330 self.do_forces = True
331 FileIOCalculator.write_input(
332 self, atoms, properties, system_changes)
333 with open(os.path.join(self.directory, 'dftb_in.hsd'), 'w') as fd:
334 self.write_dftb_in(fd)
335 write_gen(os.path.join(self.directory, 'geo_end.gen'), atoms,
336 with_tags=self.with_tags)
337 # self.atoms is none until results are read out,
338 # then it is set to the ones at writing input
339 self.atoms_input = atoms
340 self.atoms = None
341 if self.pcpot:
342 self.pcpot.write_mmcharges('dftb_external_charges.dat')
344 def read_results(self):
345 """ all results are read from results.tag file
346 It will be destroyed after it is read to avoid
347 reading it once again after some runtime error """
349 with open(os.path.join(self.directory, 'results.tag')) as fd:
350 self.lines = fd.readlines()
352 self.atoms = self.atoms_input
353 charges, energy, dipole = self.read_charges_energy_dipole()
354 if charges is not None:
355 self.results['charges'] = charges
356 self.results['energy'] = energy
357 if dipole is not None:
358 self.results['dipole'] = dipole
359 if self.do_forces:
360 forces = self.read_forces()
361 self.results['forces'] = forces
362 self.mmpositions = None
364 # stress stuff begins
365 sstring = 'stress'
366 have_stress = False
367 stress = []
368 for iline, line in enumerate(self.lines):
369 if sstring in line:
370 have_stress = True
371 start = iline + 1
372 end = start + 3
373 for i in range(start, end):
374 cell = [float(x) for x in self.lines[i].split()]
375 stress.append(cell)
376 if have_stress:
377 stress = -np.array(stress) * Hartree / Bohr**3
378 self.results['stress'] = stress.flat[[0, 4, 8, 5, 2, 1]]
379 # stress stuff ends
381 # eigenvalues and fermi levels
382 fermi_levels = self.read_fermi_levels()
383 if fermi_levels is not None:
384 self.results['fermi_levels'] = fermi_levels
386 eigenvalues = self.read_eigenvalues()
387 if eigenvalues is not None:
388 self.results['eigenvalues'] = eigenvalues
390 # calculation was carried out with atoms written in write_input
391 os.remove(os.path.join(self.directory, 'results.tag'))
393 def read_forces(self):
394 """Read Forces from dftb output file (results.tag)."""
395 from ase.units import Bohr, Hartree
397 # Initialise the indices so their scope
398 # reaches outside of the for loop
399 index_force_begin = -1
400 index_force_end = -1
402 # Force line indexes
403 for iline, line in enumerate(self.lines):
404 fstring = 'forces '
405 if line.find(fstring) >= 0:
406 index_force_begin = iline + 1
407 line1 = line.replace(':', ',')
408 index_force_end = iline + 1 + \
409 int(line1.split(',')[-1])
410 break
412 gradients = []
413 for j in range(index_force_begin, index_force_end):
414 word = self.lines[j].split()
415 gradients.append([float(word[k]) for k in range(3)])
417 return np.array(gradients) * Hartree / Bohr
419 def read_charges_energy_dipole(self):
420 """Get partial charges on atoms
421 in case we cannot find charges they are set to None
422 """
423 with open(os.path.join(self.directory, 'detailed.out')) as fd:
424 lines = fd.readlines()
426 for line in lines:
427 if line.strip().startswith('Total energy:'):
428 energy = float(line.split()[2]) * Hartree
429 break
431 qm_charges = []
432 for n, line in enumerate(lines):
433 if ('Atom' and 'Charge' in line):
434 chargestart = n + 1
435 break
436 else:
437 # print('Warning: did not find DFTB-charges')
438 # print('This is ok if flag SCC=No')
439 return None, energy, None
441 lines1 = lines[chargestart:(chargestart + len(self.atoms))]
442 for line in lines1:
443 qm_charges.append(float(line.split()[-1]))
445 dipole = None
446 for line in lines:
447 if 'Dipole moment:' in line and 'au' in line:
448 line = line.replace("Dipole moment:", "").replace("au", "")
449 dipole = np.array(line.split(), dtype=float) * Bohr
451 return np.array(qm_charges), energy, dipole
453 def get_charges(self, atoms):
454 """ Get the calculated charges
455 this is inhereted to atoms object """
456 if 'charges' in self.results:
457 return self.results['charges']
458 else:
459 return None
461 def read_eigenvalues(self):
462 """ Read Eigenvalues from dftb output file (results.tag).
463 Unfortunately, the order seems to be scrambled. """
464 # Eigenvalue line indexes
465 index_eig_begin = None
466 for iline, line in enumerate(self.lines):
467 fstring = 'eigenvalues '
468 if line.find(fstring) >= 0:
469 index_eig_begin = iline + 1
470 line1 = line.replace(':', ',')
471 ncol, nband, nkpt, nspin = map(int, line1.split(',')[-4:])
472 break
473 else:
474 return None
476 # Take into account that the last row may lack
477 # columns if nkpt * nspin * nband % ncol != 0
478 nrow = int(np.ceil(nkpt * nspin * nband * 1. / ncol))
479 index_eig_end = index_eig_begin + nrow
480 ncol_last = len(self.lines[index_eig_end - 1].split())
481 # XXX dirty fix
482 self.lines[index_eig_end - 1] = (
483 self.lines[index_eig_end - 1].strip()
484 + ' 0.0 ' * (ncol - ncol_last))
486 eig = np.loadtxt(self.lines[index_eig_begin:index_eig_end]).flatten()
487 eig *= Hartree
488 N = nkpt * nband
489 eigenvalues = [eig[i * N:(i + 1) * N].reshape((nkpt, nband))
490 for i in range(nspin)]
492 return eigenvalues
494 def read_fermi_levels(self):
495 """ Read Fermi level(s) from dftb output file (results.tag). """
496 # Fermi level line indexes
497 for iline, line in enumerate(self.lines):
498 fstring = 'fermi_level '
499 if line.find(fstring) >= 0:
500 index_fermi = iline + 1
501 break
502 else:
503 return None
505 fermi_levels = []
506 words = self.lines[index_fermi].split()
507 assert len(words) in [1, 2], 'Expected either 1 or 2 Fermi levels'
509 for word in words:
510 e = float(word)
511 # In non-spin-polarized calculations with DFTB+ v17.1,
512 # two Fermi levels are given, with the second one being 0,
513 # but we don't want to add that one to the list
514 if abs(e) > 1e-8:
515 fermi_levels.append(e)
517 return np.array(fermi_levels) * Hartree
519 def get_ibz_k_points(self):
520 return self.kpts_coord.copy()
522 def get_number_of_spins(self):
523 return self.nspin
525 def get_eigenvalues(self, kpt=0, spin=0):
526 return self.results['eigenvalues'][spin][kpt].copy()
528 def get_fermi_levels(self):
529 return self.results['fermi_levels'].copy()
531 def get_fermi_level(self):
532 return max(self.get_fermi_levels())
534 def embed(self, mmcharges=None, directory='./'):
535 """Embed atoms in point-charges (mmcharges)
536 """
537 self.pcpot = PointChargePotential(mmcharges, self.directory)
538 return self.pcpot
541class PointChargePotential:
542 def __init__(self, mmcharges, directory='./'):
543 """Point-charge potential for DFTB+.
544 """
545 self.mmcharges = mmcharges
546 self.directory = directory
547 self.mmpositions = None
548 self.mmforces = None
550 def set_positions(self, mmpositions):
551 self.mmpositions = mmpositions
553 def set_charges(self, mmcharges):
554 self.mmcharges = mmcharges
556 def write_mmcharges(self, filename):
557 """ mok all
558 write external charges as monopoles for dftb+.
560 """
561 if self.mmcharges is None:
562 print("DFTB: Warning: not writing exernal charges ")
563 return
564 with open(os.path.join(self.directory, filename), 'w') as charge_file:
565 for [pos, charge] in zip(self.mmpositions, self.mmcharges):
566 [x, y, z] = pos
567 charge_file.write('%12.6f %12.6f %12.6f %12.6f \n'
568 % (x, y, z, charge))
570 def get_forces(self, calc, get_forces=True):
571 """ returns forces on point charges if the flag get_forces=True """
572 if get_forces:
573 return self.read_forces_on_pointcharges()
574 else:
575 return np.zeros_like(self.mmpositions)
577 def read_forces_on_pointcharges(self):
578 """Read Forces from dftb output file (results.tag)."""
579 from ase.units import Bohr, Hartree
580 with open(os.path.join(self.directory, 'detailed.out')) as fd:
581 lines = fd.readlines()
583 external_forces = []
584 for n, line in enumerate(lines):
585 if ('Forces on external charges' in line):
586 chargestart = n + 1
587 break
588 else:
589 raise RuntimeError(
590 'Problem in reading forces on MM external-charges')
591 lines1 = lines[chargestart:(chargestart + len(self.mmcharges))]
592 for line in lines1:
593 external_forces.append(
594 [float(i) for i in line.split()])
595 return np.array(external_forces) * Hartree / Bohr
598def read_max_angular_momentum(path):
599 """Read maximum angular momentum from .skf file.
601 See dftb.org for A detailed description of the Slater-Koster file format.
602 """
603 with open(path) as fd:
604 line = fd.readline()
605 if line[0] == '@':
606 # Extended format
607 fd.readline()
608 l = 3
609 pos = 9
610 else:
611 # Simple format:
612 l = 2
613 pos = 7
615 # Sometimes there ar commas, sometimes not:
616 line = fd.readline().replace(',', ' ')
618 occs = [float(f) for f in line.split()[pos:pos + l + 1]]
619 for f in occs:
620 if f > 0.0:
621 return l
622 l -= 1