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

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 if kwargs.get('Hamiltonian_', 'DFTB') == 'DFTB':
96 self.default_parameters = dict(
97 Hamiltonian_='DFTB',
98 Hamiltonian_SlaterKosterFiles_='Type2FileNames',
99 Hamiltonian_SlaterKosterFiles_Prefix=self.slako_dir,
100 Hamiltonian_SlaterKosterFiles_Separator='"-"',
101 Hamiltonian_SlaterKosterFiles_Suffix='".skf"',
102 Hamiltonian_MaxAngularMomentum_='',
103 Options_='',
104 Options_WriteResultsTag='Yes')
105 else:
106 self.default_parameters = dict(
107 Options_='',
108 Options_WriteResultsTag='Yes')
110 self.pcpot = None
111 self.lines = None
112 self.atoms = None
113 self.atoms_input = None
114 self.do_forces = False
115 self.outfilename = 'dftb.out'
117 FileIOCalculator.__init__(self, restart, ignore_bad_restart_file,
118 label, atoms,
119 **kwargs)
121 # Determine number of spin channels
122 try:
123 entry = kwargs['Hamiltonian_SpinPolarisation']
124 spinpol = 'colinear' in entry.lower()
125 except KeyError:
126 spinpol = False
127 self.nspin = 2 if spinpol else 1
129 # kpoint stuff by ase
130 self.kpts = kpts
131 self.kpts_coord = None
133 if self.kpts is not None:
134 initkey = 'Hamiltonian_KPointsAndWeights'
135 mp_mesh = None
136 offsets = None
138 if isinstance(self.kpts, dict):
139 if 'path' in self.kpts:
140 # kpts is path in Brillouin zone
141 self.parameters[initkey + '_'] = 'Klines '
142 self.kpts_coord = kpts2ndarray(self.kpts, atoms=atoms)
143 else:
144 # kpts is (implicit) definition of
145 # Monkhorst-Pack grid
146 self.parameters[initkey + '_'] = 'SupercellFolding '
147 mp_mesh, offsets = kpts2sizeandoffsets(atoms=atoms,
148 **self.kpts)
149 elif np.array(self.kpts).ndim == 1:
150 # kpts is Monkhorst-Pack grid
151 self.parameters[initkey + '_'] = 'SupercellFolding '
152 mp_mesh = self.kpts
153 offsets = [0.] * 3
154 elif np.array(self.kpts).ndim == 2:
155 # kpts is (N x 3) list/array of k-point coordinates
156 # each will be given equal weight
157 self.parameters[initkey + '_'] = ''
158 self.kpts_coord = np.array(self.kpts)
159 else:
160 raise ValueError('Illegal kpts definition:' + str(self.kpts))
162 if mp_mesh is not None:
163 eps = 1e-10
164 for i in range(3):
165 key = initkey + '_empty%03d' % i
166 val = [mp_mesh[i] if j == i else 0 for j in range(3)]
167 self.parameters[key] = ' '.join(map(str, val))
168 offsets[i] *= mp_mesh[i]
169 assert abs(offsets[i]) < eps or abs(offsets[i] - 0.5) < eps
170 # DFTB+ uses a different offset convention, where
171 # the k-point mesh is already Gamma-centered prior
172 # to the addition of any offsets
173 if mp_mesh[i] % 2 == 0:
174 offsets[i] += 0.5
175 key = initkey + '_empty%03d' % 3
176 self.parameters[key] = ' '.join(map(str, offsets))
178 elif self.kpts_coord is not None:
179 for i, c in enumerate(self.kpts_coord):
180 key = initkey + '_empty%09d' % i
181 c_str = ' '.join(map(str, c))
182 if 'Klines' in self.parameters[initkey + '_']:
183 c_str = '1 ' + c_str
184 else:
185 c_str += ' 1.0'
186 self.parameters[key] = c_str
188 def write_dftb_in(self, outfile):
189 """ Write the innput file for the dftb+ calculation.
190 Geometry is taken always from the file 'geo_end.gen'.
191 """
193 outfile.write('Geometry = GenFormat { \n')
194 outfile.write(' <<< "geo_end.gen" \n')
195 outfile.write('} \n')
196 outfile.write(' \n')
198 params = self.parameters.copy()
200 s = 'Hamiltonian_MaxAngularMomentum_'
201 for key in params:
202 if key.startswith(s) and len(key) > len(s):
203 break
204 else:
205 if params.get('Hamiltonian_', 'DFTB') == 'DFTB':
206 # User didn't specify max angular mometa. Get them from
207 # the .skf files:
208 symbols = set(self.atoms.get_chemical_symbols())
209 for symbol in symbols:
210 path = os.path.join(self.slako_dir,
211 '{0}-{0}.skf'.format(symbol))
212 l = read_max_angular_momentum(path)
213 params[s + symbol] = '"{}"'.format('spdf'[l])
215 # --------MAIN KEYWORDS-------
216 previous_key = 'dummy_'
217 myspace = ' '
218 for key, value in sorted(params.items()):
219 current_depth = key.rstrip('_').count('_')
220 previous_depth = previous_key.rstrip('_').count('_')
221 for my_backsclash in reversed(
222 range(previous_depth - current_depth)):
223 outfile.write(3 * (1 + my_backsclash) * myspace + '} \n')
224 outfile.write(3 * current_depth * myspace)
225 if key.endswith('_') and len(value) > 0:
226 outfile.write(key.rstrip('_').rsplit('_')[-1] +
227 ' = ' + str(value) + '{ \n')
228 elif (key.endswith('_') and (len(value) == 0)
229 and current_depth == 0): # E.g. 'Options {'
230 outfile.write(key.rstrip('_').rsplit('_')[-1] +
231 ' ' + str(value) + '{ \n')
232 elif (key.endswith('_') and (len(value) == 0)
233 and current_depth > 0): # E.g. 'Hamiltonian_Max... = {'
234 outfile.write(key.rstrip('_').rsplit('_')[-1] +
235 ' = ' + str(value) + '{ \n')
236 elif key.count('_empty') == 1:
237 outfile.write(str(value) + ' \n')
238 elif ((key == 'Hamiltonian_ReadInitialCharges') and
239 (str(value).upper() == 'YES')):
240 f1 = os.path.isfile(self.directory + os.sep + 'charges.dat')
241 f2 = os.path.isfile(self.directory + os.sep + 'charges.bin')
242 if not (f1 or f2):
243 print('charges.dat or .bin not found, switching off guess')
244 value = 'No'
245 outfile.write(key.rsplit('_')[-1] + ' = ' + str(value) + ' \n')
246 else:
247 outfile.write(key.rsplit('_')[-1] + ' = ' + str(value) + ' \n')
248 if self.pcpot is not None and ('DFTB' in str(value)):
249 outfile.write(' ElectricField = { \n')
250 outfile.write(' PointCharges = { \n')
251 outfile.write(
252 ' CoordsAndCharges [Angstrom] = DirectRead { \n')
253 outfile.write(' Records = ' +
254 str(len(self.pcpot.mmcharges)) + ' \n')
255 outfile.write(
256 ' File = "dftb_external_charges.dat" \n')
257 outfile.write(' } \n')
258 outfile.write(' } \n')
259 outfile.write(' } \n')
260 previous_key = key
261 current_depth = key.rstrip('_').count('_')
262 for my_backsclash in reversed(range(current_depth)):
263 outfile.write(3 * my_backsclash * myspace + '} \n')
264 outfile.write('ParserOptions { \n')
265 outfile.write(' IgnoreUnprocessedNodes = Yes \n')
266 outfile.write('} \n')
267 if self.do_forces:
268 outfile.write('Analysis { \n')
269 outfile.write(' CalculateForces = Yes \n')
270 outfile.write('} \n')
272 def check_state(self, atoms):
273 system_changes = FileIOCalculator.check_state(self, atoms)
274 # Ignore unit cell for molecules:
275 if not atoms.pbc.any() and 'cell' in system_changes:
276 system_changes.remove('cell')
277 if self.pcpot and self.pcpot.mmpositions is not None:
278 system_changes.append('positions')
279 return system_changes
281 def write_input(self, atoms, properties=None, system_changes=None):
282 from ase.io import write
283 if properties is not None:
284 if 'forces' in properties or 'stress' in properties:
285 self.do_forces = True
286 FileIOCalculator.write_input(
287 self, atoms, properties, system_changes)
288 with open(os.path.join(self.directory, 'dftb_in.hsd'), 'w') as fd:
289 self.write_dftb_in(fd)
290 write(os.path.join(self.directory, 'geo_end.gen'), atoms,
291 parallel=False)
292 # self.atoms is none until results are read out,
293 # then it is set to the ones at writing input
294 self.atoms_input = atoms
295 self.atoms = None
296 if self.pcpot:
297 self.pcpot.write_mmcharges('dftb_external_charges.dat')
299 def read_results(self):
300 """ all results are read from results.tag file
301 It will be destroyed after it is read to avoid
302 reading it once again after some runtime error """
304 with open(os.path.join(self.directory, 'results.tag'), 'r') as fd:
305 self.lines = fd.readlines()
307 self.atoms = self.atoms_input
308 charges, energy, dipole = self.read_charges_energy_dipole()
309 if charges is not None:
310 self.results['charges'] = charges
311 self.results['energy'] = energy
312 if dipole is not None:
313 self.results['dipole'] = dipole
314 if self.do_forces:
315 forces = self.read_forces()
316 self.results['forces'] = forces
317 self.mmpositions = None
319 # stress stuff begins
320 sstring = 'stress'
321 have_stress = False
322 stress = list()
323 for iline, line in enumerate(self.lines):
324 if sstring in line:
325 have_stress = True
326 start = iline + 1
327 end = start + 3
328 for i in range(start, end):
329 cell = [float(x) for x in self.lines[i].split()]
330 stress.append(cell)
331 if have_stress:
332 stress = -np.array(stress) * Hartree / Bohr**3
333 self.results['stress'] = stress.flat[[0, 4, 8, 5, 2, 1]]
334 # stress stuff ends
336 # eigenvalues and fermi levels
337 fermi_levels = self.read_fermi_levels()
338 if fermi_levels is not None:
339 self.results['fermi_levels'] = fermi_levels
341 eigenvalues = self.read_eigenvalues()
342 if eigenvalues is not None:
343 self.results['eigenvalues'] = eigenvalues
345 # calculation was carried out with atoms written in write_input
346 os.remove(os.path.join(self.directory, 'results.tag'))
348 def read_forces(self):
349 """Read Forces from dftb output file (results.tag)."""
350 from ase.units import Hartree, Bohr
352 # Initialise the indices so their scope
353 # reaches outside of the for loop
354 index_force_begin = -1
355 index_force_end = -1
357 # Force line indexes
358 for iline, line in enumerate(self.lines):
359 fstring = 'forces '
360 if line.find(fstring) >= 0:
361 index_force_begin = iline + 1
362 line1 = line.replace(':', ',')
363 index_force_end = iline + 1 + \
364 int(line1.split(',')[-1])
365 break
367 gradients = []
368 for j in range(index_force_begin, index_force_end):
369 word = self.lines[j].split()
370 gradients.append([float(word[k]) for k in range(0, 3)])
372 return np.array(gradients) * Hartree / Bohr
374 def read_charges_energy_dipole(self):
375 """Get partial charges on atoms
376 in case we cannot find charges they are set to None
377 """
378 with open(os.path.join(self.directory, 'detailed.out'), 'r') as fd:
379 lines = fd.readlines()
381 for line in lines:
382 if line.strip().startswith('Total energy:'):
383 energy = float(line.split()[2]) * Hartree
384 break
386 qm_charges = []
387 for n, line in enumerate(lines):
388 if ('Atom' and 'Charge' in line):
389 chargestart = n + 1
390 break
391 else:
392 # print('Warning: did not find DFTB-charges')
393 # print('This is ok if flag SCC=No')
394 return None, energy, None
396 lines1 = lines[chargestart:(chargestart + len(self.atoms))]
397 for line in lines1:
398 qm_charges.append(float(line.split()[-1]))
400 dipole = None
401 for line in lines:
402 if 'Dipole moment:' in line and 'au' in line:
403 words = line.split()
404 dipole = np.array(
405 [float(w) for w in words[-4:-1]]) * Bohr
407 return np.array(qm_charges), energy, dipole
409 def get_charges(self, atoms):
410 """ Get the calculated charges
411 this is inhereted to atoms object """
412 if 'charges' in self.results:
413 return self.results['charges']
414 else:
415 return None
417 def read_eigenvalues(self):
418 """ Read Eigenvalues from dftb output file (results.tag).
419 Unfortunately, the order seems to be scrambled. """
420 # Eigenvalue line indexes
421 index_eig_begin = None
422 for iline, line in enumerate(self.lines):
423 fstring = 'eigenvalues '
424 if line.find(fstring) >= 0:
425 index_eig_begin = iline + 1
426 line1 = line.replace(':', ',')
427 ncol, nband, nkpt, nspin = map(int, line1.split(',')[-4:])
428 break
429 else:
430 return None
432 # Take into account that the last row may lack
433 # columns if nkpt * nspin * nband % ncol != 0
434 nrow = int(np.ceil(nkpt * nspin * nband * 1. / ncol))
435 index_eig_end = index_eig_begin + nrow
436 ncol_last = len(self.lines[index_eig_end - 1].split())
437 # XXX dirty fix
438 self.lines[index_eig_end - 1] = (
439 self.lines[index_eig_end - 1].strip()
440 + ' 0.0 ' * (ncol - ncol_last))
442 eig = np.loadtxt(self.lines[index_eig_begin:index_eig_end]).flatten()
443 eig *= Hartree
444 N = nkpt * nband
445 eigenvalues = [eig[i * N:(i + 1) * N].reshape((nkpt, nband))
446 for i in range(nspin)]
448 return eigenvalues
450 def read_fermi_levels(self):
451 """ Read Fermi level(s) from dftb output file (results.tag). """
452 # Fermi level line indexes
453 for iline, line in enumerate(self.lines):
454 fstring = 'fermi_level '
455 if line.find(fstring) >= 0:
456 index_fermi = iline + 1
457 break
458 else:
459 return None
461 fermi_levels = []
462 words = self.lines[index_fermi].split()
463 assert len(words) in [1, 2], 'Expected either 1 or 2 Fermi levels'
465 for word in words:
466 e = float(word)
467 # In non-spin-polarized calculations with DFTB+ v17.1,
468 # two Fermi levels are given, with the second one being 0,
469 # but we don't want to add that one to the list
470 if abs(e) > 1e-8:
471 fermi_levels.append(e)
473 return np.array(fermi_levels) * Hartree
475 def get_ibz_k_points(self):
476 return self.kpts_coord.copy()
478 def get_number_of_spins(self):
479 return self.nspin
481 def get_eigenvalues(self, kpt=0, spin=0):
482 return self.results['eigenvalues'][spin][kpt].copy()
484 def get_fermi_levels(self):
485 return self.results['fermi_levels'].copy()
487 def get_fermi_level(self):
488 return max(self.get_fermi_levels())
490 def embed(self, mmcharges=None, directory='./'):
491 """Embed atoms in point-charges (mmcharges)
492 """
493 self.pcpot = PointChargePotential(mmcharges, self.directory)
494 return self.pcpot
497class PointChargePotential:
498 def __init__(self, mmcharges, directory='./'):
499 """Point-charge potential for DFTB+.
500 """
501 self.mmcharges = mmcharges
502 self.directory = directory
503 self.mmpositions = None
504 self.mmforces = None
506 def set_positions(self, mmpositions):
507 self.mmpositions = mmpositions
509 def set_charges(self, mmcharges):
510 self.mmcharges = mmcharges
512 def write_mmcharges(self, filename):
513 """ mok all
514 write external charges as monopoles for dftb+.
516 """
517 if self.mmcharges is None:
518 print("DFTB: Warning: not writing exernal charges ")
519 return
520 with open(os.path.join(self.directory, filename), 'w') as charge_file:
521 for [pos, charge] in zip(self.mmpositions, self.mmcharges):
522 [x, y, z] = pos
523 charge_file.write('%12.6f %12.6f %12.6f %12.6f \n'
524 % (x, y, z, charge))
526 def get_forces(self, calc, get_forces=True):
527 """ returns forces on point charges if the flag get_forces=True """
528 if get_forces:
529 return self.read_forces_on_pointcharges()
530 else:
531 return np.zeros_like(self.mmpositions)
533 def read_forces_on_pointcharges(self):
534 """Read Forces from dftb output file (results.tag)."""
535 from ase.units import Hartree, Bohr
536 with open(os.path.join(self.directory, 'detailed.out'), 'r') as fd:
537 lines = fd.readlines()
539 external_forces = []
540 for n, line in enumerate(lines):
541 if ('Forces on external charges' in line):
542 chargestart = n + 1
543 break
544 else:
545 raise RuntimeError(
546 'Problem in reading forces on MM external-charges')
547 lines1 = lines[chargestart:(chargestart + len(self.mmcharges))]
548 for line in lines1:
549 external_forces.append(
550 [float(i) for i in line.split()])
551 return np.array(external_forces) * Hartree / Bohr
554def read_max_angular_momentum(path):
555 """Read maximum angular momentum from .skf file.
557 See dftb.org for A detailed description of the Slater-Koster file format.
558 """
559 with open(path, 'r') as fd:
560 line = fd.readline()
561 if line[0] == '@':
562 # Extended format
563 fd.readline()
564 l = 3
565 pos = 9
566 else:
567 # Simple format:
568 l = 2
569 pos = 7
571 # Sometimes there ar commas, sometimes not:
572 line = fd.readline().replace(',', ' ')
574 occs = [float(f) for f in line.split()[pos:pos + l + 1]]
575 for f in occs:
576 if f > 0.0:
577 return l
578 l -= 1