1""" This module defines a FileIOCalculator for DFTB+ 

2 

3http://www.dftbplus.org/ 

4http://www.dftb.org/ 

5 

6Initial development: markus.kaukonen@iki.fi 

7""" 

8 

9import os 

10import numpy as np 

11from ase.calculators.calculator import (FileIOCalculator, kpts2ndarray, 

12 kpts2sizeandoffsets) 

13from ase.units import Hartree, Bohr 

14 

15 

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' 

21 

22 implemented_properties = ['energy', 'forces', 'charges', 

23 'stress', 'dipole'] 

24 discard_results_on_any_change = True 

25 

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: 

34 

35 >>> Hamiltonian = DFTB { 

36 >>> SCC = Yes 

37 >>> SCCTolerance = 1e-8 

38 >>> MaxAngularMomentum = { 

39 >>> H = s 

40 >>> O = p 

41 >>> } 

42 >>> } 

43 

44 This can be generated by the DFTB+ calculator by using the 

45 following settings: 

46 

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') 

53 

54 In addition to keywords specific to DFTB+, also the following keywords 

55 arguments can be used: 

56 

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: 

71 

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) 

81 

82 Additional attribute to be set by the embed() method: 

83 

84 pcpot: PointCharge object 

85 An external point charge potential (for QM/MM calculations) 

86 """ 

87 

88 if slako_dir is None: 

89 slako_dir = os.environ.get('DFTB_PREFIX', './') 

90 if not slako_dir.endswith('/'): 

91 slako_dir += '/' 

92 

93 self.slako_dir = slako_dir 

94 

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') 

109 

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' 

116 

117 FileIOCalculator.__init__(self, restart, ignore_bad_restart_file, 

118 label, atoms, 

119 **kwargs) 

120 

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 

128 

129 # kpoint stuff by ase 

130 self.kpts = kpts 

131 self.kpts_coord = None 

132 

133 if self.kpts is not None: 

134 initkey = 'Hamiltonian_KPointsAndWeights' 

135 mp_mesh = None 

136 offsets = None 

137 

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)) 

161 

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)) 

177 

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 

187 

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 """ 

192 

193 outfile.write('Geometry = GenFormat { \n') 

194 outfile.write(' <<< "geo_end.gen" \n') 

195 outfile.write('} \n') 

196 outfile.write(' \n') 

197 

198 params = self.parameters.copy() 

199 

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]) 

214 

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') 

271 

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 

280 

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') 

298 

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 """ 

303 

304 with open(os.path.join(self.directory, 'results.tag'), 'r') as fd: 

305 self.lines = fd.readlines() 

306 

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 

318 

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 

335 

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 

340 

341 eigenvalues = self.read_eigenvalues() 

342 if eigenvalues is not None: 

343 self.results['eigenvalues'] = eigenvalues 

344 

345 # calculation was carried out with atoms written in write_input 

346 os.remove(os.path.join(self.directory, 'results.tag')) 

347 

348 def read_forces(self): 

349 """Read Forces from dftb output file (results.tag).""" 

350 from ase.units import Hartree, Bohr 

351 

352 # Initialise the indices so their scope 

353 # reaches outside of the for loop 

354 index_force_begin = -1 

355 index_force_end = -1 

356 

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 

366 

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)]) 

371 

372 return np.array(gradients) * Hartree / Bohr 

373 

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() 

380 

381 for line in lines: 

382 if line.strip().startswith('Total energy:'): 

383 energy = float(line.split()[2]) * Hartree 

384 break 

385 

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 

395 

396 lines1 = lines[chargestart:(chargestart + len(self.atoms))] 

397 for line in lines1: 

398 qm_charges.append(float(line.split()[-1])) 

399 

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 

406 

407 return np.array(qm_charges), energy, dipole 

408 

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 

416 

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 

431 

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)) 

441 

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)] 

447 

448 return eigenvalues 

449 

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 

460 

461 fermi_levels = [] 

462 words = self.lines[index_fermi].split() 

463 assert len(words) in [1, 2], 'Expected either 1 or 2 Fermi levels' 

464 

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) 

472 

473 return np.array(fermi_levels) * Hartree 

474 

475 def get_ibz_k_points(self): 

476 return self.kpts_coord.copy() 

477 

478 def get_number_of_spins(self): 

479 return self.nspin 

480 

481 def get_eigenvalues(self, kpt=0, spin=0): 

482 return self.results['eigenvalues'][spin][kpt].copy() 

483 

484 def get_fermi_levels(self): 

485 return self.results['fermi_levels'].copy() 

486 

487 def get_fermi_level(self): 

488 return max(self.get_fermi_levels()) 

489 

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 

495 

496 

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 

505 

506 def set_positions(self, mmpositions): 

507 self.mmpositions = mmpositions 

508 

509 def set_charges(self, mmcharges): 

510 self.mmcharges = mmcharges 

511 

512 def write_mmcharges(self, filename): 

513 """ mok all 

514 write external charges as monopoles for dftb+. 

515 

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)) 

525 

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) 

532 

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() 

538 

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 

552 

553 

554def read_max_angular_momentum(path): 

555 """Read maximum angular momentum from .skf file. 

556 

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 

570 

571 # Sometimes there ar commas, sometimes not: 

572 line = fd.readline().replace(',', ' ') 

573 

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