Coverage for /builds/ase/ase/ase/calculators/dftb.py : 74.07%

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""" This module defines a FileIOCalculator for DFTB+
3http://www.dftbplus.org/
4http://www.dftb.org/
6Initial development: markus.kaukonen@iki.fi
7"""
9import os
10import numpy as np
11from ase.calculators.calculator import (FileIOCalculator, kpts2ndarray,
12 kpts2sizeandoffsets)
13from ase.units import Hartree, Bohr
16class Dftb(FileIOCalculator):
17 if 'DFTB_COMMAND' in os.environ:
18 command = os.environ['DFTB_COMMAND'] + ' > PREFIX.out'
19 else:
20 command = 'dftb+ > PREFIX.out'
22 implemented_properties = ['energy', 'forces', 'charges',
23 'stress', 'dipole']
24 discard_results_on_any_change = True
26 def __init__(self, restart=None,
27 ignore_bad_restart_file=FileIOCalculator._deprecated,
28 label='dftb', atoms=None, kpts=None,
29 slako_dir=None,
30 **kwargs):
31 """
32 All keywords for the dftb_in.hsd input file (see the DFTB+ manual)
33 can be set by ASE. Consider the following input file block:
35 >>> Hamiltonian = DFTB {
36 >>> SCC = Yes
37 >>> SCCTolerance = 1e-8
38 >>> MaxAngularMomentum = {
39 >>> H = s
40 >>> O = p
41 >>> }
42 >>> }
44 This can be generated by the DFTB+ calculator by using the
45 following settings:
47 >>> calc = Dftb(Hamiltonian_='DFTB', # line is included by default
48 >>> Hamiltonian_SCC='Yes',
49 >>> Hamiltonian_SCCTolerance=1e-8,
50 >>> Hamiltonian_MaxAngularMomentum_='',
51 >>> Hamiltonian_MaxAngularMomentum_H='s',
52 >>> Hamiltonian_MaxAngularMomentum_O='p')
54 In addition to keywords specific to DFTB+, also the following keywords
55 arguments can be used:
57 restart: str
58 Prefix for restart file. May contain a directory.
59 Default is None: don't restart.
60 ignore_bad_restart_file: bool
61 Ignore broken or missing restart file. By default, it is an
62 error if the restart file is missing or broken.
63 label: str (default 'dftb')
64 Prefix used for the main output file (<label>.out).
65 atoms: Atoms object (default None)
66 Optional Atoms object to which the calculator will be
67 attached. When restarting, atoms will get its positions and
68 unit-cell updated from file.
69 kpts: (default None)
70 Brillouin zone sampling:
72 * ``(1,1,1)`` or ``None``: Gamma-point only
73 * ``(n1,n2,n3)``: Monkhorst-Pack grid
74 * ``dict``: Interpreted as a path in the Brillouin zone if
75 it contains the 'path_' keyword. Otherwise it is converted
76 into a Monkhorst-Pack grid using
77 ``ase.calculators.calculator.kpts2sizeandoffsets``
78 * ``[(k11,k12,k13),(k21,k22,k23),...]``: Explicit (Nkpts x 3)
79 array of k-points in units of the reciprocal lattice vectors
80 (each with equal weight)
82 Additional attribute to be set by the embed() method:
84 pcpot: PointCharge object
85 An external point charge potential (for QM/MM calculations)
86 """
88 if slako_dir is None:
89 slako_dir = os.environ.get('DFTB_PREFIX', './')
90 if not slako_dir.endswith('/'):
91 slako_dir += '/'
93 self.slako_dir = slako_dir
95 self.default_parameters = dict(
96 Hamiltonian_='DFTB',
97 Hamiltonian_SlaterKosterFiles_='Type2FileNames',
98 Hamiltonian_SlaterKosterFiles_Prefix=self.slako_dir,
99 Hamiltonian_SlaterKosterFiles_Separator='"-"',
100 Hamiltonian_SlaterKosterFiles_Suffix='".skf"',
101 Hamiltonian_MaxAngularMomentum_='',
102 Options_='',
103 Options_WriteResultsTag='Yes')
105 self.pcpot = None
106 self.lines = None
107 self.atoms = None
108 self.atoms_input = None
109 self.do_forces = False
110 self.outfilename = 'dftb.out'
112 FileIOCalculator.__init__(self, restart, ignore_bad_restart_file,
113 label, atoms,
114 **kwargs)
116 # Determine number of spin channels
117 try:
118 entry = kwargs['Hamiltonian_SpinPolarisation']
119 spinpol = 'colinear' in entry.lower()
120 except KeyError:
121 spinpol = False
122 self.nspin = 2 if spinpol else 1
124 # kpoint stuff by ase
125 self.kpts = kpts
126 self.kpts_coord = None
128 if self.kpts is not None:
129 initkey = 'Hamiltonian_KPointsAndWeights'
130 mp_mesh = None
131 offsets = None
133 if isinstance(self.kpts, dict):
134 if 'path' in self.kpts:
135 # kpts is path in Brillouin zone
136 self.parameters[initkey + '_'] = 'Klines '
137 self.kpts_coord = kpts2ndarray(self.kpts, atoms=atoms)
138 else:
139 # kpts is (implicit) definition of
140 # Monkhorst-Pack grid
141 self.parameters[initkey + '_'] = 'SupercellFolding '
142 mp_mesh, offsets = kpts2sizeandoffsets(atoms=atoms,
143 **self.kpts)
144 elif np.array(self.kpts).ndim == 1:
145 # kpts is Monkhorst-Pack grid
146 self.parameters[initkey + '_'] = 'SupercellFolding '
147 mp_mesh = self.kpts
148 offsets = [0.] * 3
149 elif np.array(self.kpts).ndim == 2:
150 # kpts is (N x 3) list/array of k-point coordinates
151 # each will be given equal weight
152 self.parameters[initkey + '_'] = ''
153 self.kpts_coord = np.array(self.kpts)
154 else:
155 raise ValueError('Illegal kpts definition:' + str(self.kpts))
157 if mp_mesh is not None:
158 eps = 1e-10
159 for i in range(3):
160 key = initkey + '_empty%03d' % i
161 val = [mp_mesh[i] if j == i else 0 for j in range(3)]
162 self.parameters[key] = ' '.join(map(str, val))
163 offsets[i] *= mp_mesh[i]
164 assert abs(offsets[i]) < eps or abs(offsets[i] - 0.5) < eps
165 # DFTB+ uses a different offset convention, where
166 # the k-point mesh is already Gamma-centered prior
167 # to the addition of any offsets
168 if mp_mesh[i] % 2 == 0:
169 offsets[i] += 0.5
170 key = initkey + '_empty%03d' % 3
171 self.parameters[key] = ' '.join(map(str, offsets))
173 elif self.kpts_coord is not None:
174 for i, c in enumerate(self.kpts_coord):
175 key = initkey + '_empty%09d' % i
176 c_str = ' '.join(map(str, c))
177 if 'Klines' in self.parameters[initkey + '_']:
178 c_str = '1 ' + c_str
179 else:
180 c_str += ' 1.0'
181 self.parameters[key] = c_str
183 def write_dftb_in(self, outfile):
184 """ Write the innput file for the dftb+ calculation.
185 Geometry is taken always from the file 'geo_end.gen'.
186 """
188 outfile.write('Geometry = GenFormat { \n')
189 outfile.write(' <<< "geo_end.gen" \n')
190 outfile.write('} \n')
191 outfile.write(' \n')
193 params = self.parameters.copy()
195 s = 'Hamiltonian_MaxAngularMomentum_'
196 for key in params:
197 if key.startswith(s) and len(key) > len(s):
198 break
199 else:
200 # User didn't specify max angular mometa. Get them from
201 # the .skf files:
202 symbols = set(self.atoms.get_chemical_symbols())
203 for symbol in symbols:
204 path = os.path.join(self.slako_dir,
205 '{0}-{0}.skf'.format(symbol))
206 l = read_max_angular_momentum(path)
207 params[s + symbol] = '"{}"'.format('spdf'[l])
209 # --------MAIN KEYWORDS-------
210 previous_key = 'dummy_'
211 myspace = ' '
212 for key, value in sorted(params.items()):
213 current_depth = key.rstrip('_').count('_')
214 previous_depth = previous_key.rstrip('_').count('_')
215 for my_backsclash in reversed(
216 range(previous_depth - current_depth)):
217 outfile.write(3 * (1 + my_backsclash) * myspace + '} \n')
218 outfile.write(3 * current_depth * myspace)
219 if key.endswith('_') and len(value) > 0:
220 outfile.write(key.rstrip('_').rsplit('_')[-1] +
221 ' = ' + str(value) + '{ \n')
222 elif (key.endswith('_') and (len(value) == 0)
223 and current_depth == 0): # E.g. 'Options {'
224 outfile.write(key.rstrip('_').rsplit('_')[-1] +
225 ' ' + str(value) + '{ \n')
226 elif (key.endswith('_') and (len(value) == 0)
227 and current_depth > 0): # E.g. 'Hamiltonian_Max... = {'
228 outfile.write(key.rstrip('_').rsplit('_')[-1] +
229 ' = ' + str(value) + '{ \n')
230 elif key.count('_empty') == 1:
231 outfile.write(str(value) + ' \n')
232 elif ((key == 'Hamiltonian_ReadInitialCharges') and
233 (str(value).upper() == 'YES')):
234 f1 = os.path.isfile(self.directory + os.sep + 'charges.dat')
235 f2 = os.path.isfile(self.directory + os.sep + 'charges.bin')
236 if not (f1 or f2):
237 print('charges.dat or .bin not found, switching off guess')
238 value = 'No'
239 outfile.write(key.rsplit('_')[-1] + ' = ' + str(value) + ' \n')
240 else:
241 outfile.write(key.rsplit('_')[-1] + ' = ' + str(value) + ' \n')
242 if self.pcpot is not None and ('DFTB' in str(value)):
243 outfile.write(' ElectricField = { \n')
244 outfile.write(' PointCharges = { \n')
245 outfile.write(
246 ' CoordsAndCharges [Angstrom] = DirectRead { \n')
247 outfile.write(' Records = ' +
248 str(len(self.pcpot.mmcharges)) + ' \n')
249 outfile.write(
250 ' File = "dftb_external_charges.dat" \n')
251 outfile.write(' } \n')
252 outfile.write(' } \n')
253 outfile.write(' } \n')
254 previous_key = key
255 current_depth = key.rstrip('_').count('_')
256 for my_backsclash in reversed(range(current_depth)):
257 outfile.write(3 * my_backsclash * myspace + '} \n')
258 outfile.write('ParserOptions { \n')
259 outfile.write(' IgnoreUnprocessedNodes = Yes \n')
260 outfile.write('} \n')
261 if self.do_forces:
262 outfile.write('Analysis { \n')
263 outfile.write(' CalculateForces = Yes \n')
264 outfile.write('} \n')
266 def check_state(self, atoms):
267 system_changes = FileIOCalculator.check_state(self, atoms)
268 # Ignore unit cell for molecules:
269 if not atoms.pbc.any() and 'cell' in system_changes:
270 system_changes.remove('cell')
271 if self.pcpot and self.pcpot.mmpositions is not None:
272 system_changes.append('positions')
273 return system_changes
275 def write_input(self, atoms, properties=None, system_changes=None):
276 from ase.io import write
277 if properties is not None:
278 if 'forces' in properties or 'stress' in properties:
279 self.do_forces = True
280 FileIOCalculator.write_input(
281 self, atoms, properties, system_changes)
282 with open(os.path.join(self.directory, 'dftb_in.hsd'), 'w') as fd:
283 self.write_dftb_in(fd)
284 write(os.path.join(self.directory, 'geo_end.gen'), atoms,
285 parallel=False)
286 # self.atoms is none until results are read out,
287 # then it is set to the ones at writing input
288 self.atoms_input = atoms
289 self.atoms = None
290 if self.pcpot:
291 self.pcpot.write_mmcharges('dftb_external_charges.dat')
293 def read_results(self):
294 """ all results are read from results.tag file
295 It will be destroyed after it is read to avoid
296 reading it once again after some runtime error """
298 with open(os.path.join(self.directory, 'results.tag'), 'r') as fd:
299 self.lines = fd.readlines()
301 self.atoms = self.atoms_input
302 charges, energy, dipole = self.read_charges_energy_dipole()
303 if charges is not None:
304 self.results['charges'] = charges
305 self.results['energy'] = energy
306 if dipole is not None:
307 self.results['dipole'] = dipole
308 if self.do_forces:
309 forces = self.read_forces()
310 self.results['forces'] = forces
311 self.mmpositions = None
313 # stress stuff begins
314 sstring = 'stress'
315 have_stress = False
316 stress = list()
317 for iline, line in enumerate(self.lines):
318 if sstring in line:
319 have_stress = True
320 start = iline + 1
321 end = start + 3
322 for i in range(start, end):
323 cell = [float(x) for x in self.lines[i].split()]
324 stress.append(cell)
325 if have_stress:
326 stress = -np.array(stress) * Hartree / Bohr**3
327 self.results['stress'] = stress.flat[[0, 4, 8, 5, 2, 1]]
328 # stress stuff ends
330 # eigenvalues and fermi levels
331 fermi_levels = self.read_fermi_levels()
332 if fermi_levels is not None:
333 self.results['fermi_levels'] = fermi_levels
335 eigenvalues = self.read_eigenvalues()
336 if eigenvalues is not None:
337 self.results['eigenvalues'] = eigenvalues
339 # calculation was carried out with atoms written in write_input
340 os.remove(os.path.join(self.directory, 'results.tag'))
342 def read_forces(self):
343 """Read Forces from dftb output file (results.tag)."""
344 from ase.units import Hartree, Bohr
346 # Initialise the indices so their scope
347 # reaches outside of the for loop
348 index_force_begin = -1
349 index_force_end = -1
351 # Force line indexes
352 for iline, line in enumerate(self.lines):
353 fstring = 'forces '
354 if line.find(fstring) >= 0:
355 index_force_begin = iline + 1
356 line1 = line.replace(':', ',')
357 index_force_end = iline + 1 + \
358 int(line1.split(',')[-1])
359 break
361 gradients = []
362 for j in range(index_force_begin, index_force_end):
363 word = self.lines[j].split()
364 gradients.append([float(word[k]) for k in range(0, 3)])
366 return np.array(gradients) * Hartree / Bohr
368 def read_charges_energy_dipole(self):
369 """Get partial charges on atoms
370 in case we cannot find charges they are set to None
371 """
372 with open(os.path.join(self.directory, 'detailed.out'), 'r') as fd:
373 lines = fd.readlines()
375 for line in lines:
376 if line.strip().startswith('Total energy:'):
377 energy = float(line.split()[2]) * Hartree
378 break
380 qm_charges = []
381 for n, line in enumerate(lines):
382 if ('Atom' and 'Charge' in line):
383 chargestart = n + 1
384 break
385 else:
386 # print('Warning: did not find DFTB-charges')
387 # print('This is ok if flag SCC=No')
388 return None, energy, None
390 lines1 = lines[chargestart:(chargestart + len(self.atoms))]
391 for line in lines1:
392 qm_charges.append(float(line.split()[-1]))
394 dipole = None
395 for line in lines:
396 if 'Dipole moment:' in line and 'au' in line:
397 words = line.split()
398 dipole = np.array(
399 [float(w) for w in words[-4:-1]]) * Bohr
401 return np.array(qm_charges), energy, dipole
403 def get_charges(self, atoms):
404 """ Get the calculated charges
405 this is inhereted to atoms object """
406 if 'charges' in self.results:
407 return self.results['charges']
408 else:
409 return None
411 def read_eigenvalues(self):
412 """ Read Eigenvalues from dftb output file (results.tag).
413 Unfortunately, the order seems to be scrambled. """
414 # Eigenvalue line indexes
415 index_eig_begin = None
416 for iline, line in enumerate(self.lines):
417 fstring = 'eigenvalues '
418 if line.find(fstring) >= 0:
419 index_eig_begin = iline + 1
420 line1 = line.replace(':', ',')
421 ncol, nband, nkpt, nspin = map(int, line1.split(',')[-4:])
422 break
423 else:
424 return None
426 # Take into account that the last row may lack
427 # columns if nkpt * nspin * nband % ncol != 0
428 nrow = int(np.ceil(nkpt * nspin * nband * 1. / ncol))
429 index_eig_end = index_eig_begin + nrow
430 ncol_last = len(self.lines[index_eig_end - 1].split())
431 # XXX dirty fix
432 self.lines[index_eig_end - 1] = (
433 self.lines[index_eig_end - 1].strip()
434 + ' 0.0 ' * (ncol - ncol_last))
436 eig = np.loadtxt(self.lines[index_eig_begin:index_eig_end]).flatten()
437 eig *= Hartree
438 N = nkpt * nband
439 eigenvalues = [eig[i * N:(i + 1) * N].reshape((nkpt, nband))
440 for i in range(nspin)]
442 return eigenvalues
444 def read_fermi_levels(self):
445 """ Read Fermi level(s) from dftb output file (results.tag). """
446 # Fermi level line indexes
447 for iline, line in enumerate(self.lines):
448 fstring = 'fermi_level '
449 if line.find(fstring) >= 0:
450 index_fermi = iline + 1
451 break
452 else:
453 return None
455 fermi_levels = []
456 words = self.lines[index_fermi].split()
457 assert len(words) in [1, 2], 'Expected either 1 or 2 Fermi levels'
459 for word in words:
460 e = float(word)
461 # In non-spin-polarized calculations with DFTB+ v17.1,
462 # two Fermi levels are given, with the second one being 0,
463 # but we don't want to add that one to the list
464 if abs(e) > 1e-8:
465 fermi_levels.append(e)
467 return np.array(fermi_levels) * Hartree
469 def get_ibz_k_points(self):
470 return self.kpts_coord.copy()
472 def get_number_of_spins(self):
473 return self.nspin
475 def get_eigenvalues(self, kpt=0, spin=0):
476 return self.results['eigenvalues'][spin][kpt].copy()
478 def get_fermi_levels(self):
479 return self.results['fermi_levels'].copy()
481 def get_fermi_level(self):
482 return max(self.get_fermi_levels())
484 def embed(self, mmcharges=None, directory='./'):
485 """Embed atoms in point-charges (mmcharges)
486 """
487 self.pcpot = PointChargePotential(mmcharges, self.directory)
488 return self.pcpot
491class PointChargePotential:
492 def __init__(self, mmcharges, directory='./'):
493 """Point-charge potential for DFTB+.
494 """
495 self.mmcharges = mmcharges
496 self.directory = directory
497 self.mmpositions = None
498 self.mmforces = None
500 def set_positions(self, mmpositions):
501 self.mmpositions = mmpositions
503 def set_charges(self, mmcharges):
504 self.mmcharges = mmcharges
506 def write_mmcharges(self, filename):
507 """ mok all
508 write external charges as monopoles for dftb+.
510 """
511 if self.mmcharges is None:
512 print("DFTB: Warning: not writing exernal charges ")
513 return
514 with open(os.path.join(self.directory, filename), 'w') as charge_file:
515 for [pos, charge] in zip(self.mmpositions, self.mmcharges):
516 [x, y, z] = pos
517 charge_file.write('%12.6f %12.6f %12.6f %12.6f \n'
518 % (x, y, z, charge))
520 def get_forces(self, calc, get_forces=True):
521 """ returns forces on point charges if the flag get_forces=True """
522 if get_forces:
523 return self.read_forces_on_pointcharges()
524 else:
525 return np.zeros_like(self.mmpositions)
527 def read_forces_on_pointcharges(self):
528 """Read Forces from dftb output file (results.tag)."""
529 from ase.units import Hartree, Bohr
530 with open(os.path.join(self.directory, 'detailed.out'), 'r') as fd:
531 lines = fd.readlines()
533 external_forces = []
534 for n, line in enumerate(lines):
535 if ('Forces on external charges' in line):
536 chargestart = n + 1
537 break
538 else:
539 raise RuntimeError(
540 'Problem in reading forces on MM external-charges')
541 lines1 = lines[chargestart:(chargestart + len(self.mmcharges))]
542 for line in lines1:
543 external_forces.append(
544 [float(i) for i in line.split()])
545 return np.array(external_forces) * Hartree / Bohr
548def read_max_angular_momentum(path):
549 """Read maximum angular momentum from .skf file.
551 See dftb.org for A detailed description of the Slater-Koster file format.
552 """
553 with open(path, 'r') as fd:
554 line = fd.readline()
555 if line[0] == '@':
556 # Extended format
557 fd.readline()
558 l = 3
559 pos = 9
560 else:
561 # Simple format:
562 l = 2
563 pos = 7
565 # Sometimes there ar commas, sometimes not:
566 line = fd.readline().replace(',', ' ')
568 occs = [float(f) for f in line.split()[pos:pos + l + 1]]
569 for f in occs:
570 if f > 0.0:
571 return l
572 l -= 1