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

104 

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' 

111 

112 FileIOCalculator.__init__(self, restart, ignore_bad_restart_file, 

113 label, atoms, 

114 **kwargs) 

115 

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 

123 

124 # kpoint stuff by ase 

125 self.kpts = kpts 

126 self.kpts_coord = None 

127 

128 if self.kpts is not None: 

129 initkey = 'Hamiltonian_KPointsAndWeights' 

130 mp_mesh = None 

131 offsets = None 

132 

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

156 

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

172 

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 

182 

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

187 

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

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

190 outfile.write('} \n') 

191 outfile.write(' \n') 

192 

193 params = self.parameters.copy() 

194 

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

208 

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

265 

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 

274 

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

292 

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

297 

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

299 self.lines = fd.readlines() 

300 

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 

312 

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 

329 

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 

334 

335 eigenvalues = self.read_eigenvalues() 

336 if eigenvalues is not None: 

337 self.results['eigenvalues'] = eigenvalues 

338 

339 # calculation was carried out with atoms written in write_input 

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

341 

342 def read_forces(self): 

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

344 from ase.units import Hartree, Bohr 

345 

346 # Initialise the indices so their scope 

347 # reaches outside of the for loop 

348 index_force_begin = -1 

349 index_force_end = -1 

350 

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 

360 

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

365 

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

367 

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

374 

375 for line in lines: 

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

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

378 break 

379 

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 

389 

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

391 for line in lines1: 

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

393 

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 

400 

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

402 

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 

410 

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 

425 

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

435 

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

441 

442 return eigenvalues 

443 

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 

454 

455 fermi_levels = [] 

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

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

458 

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) 

466 

467 return np.array(fermi_levels) * Hartree 

468 

469 def get_ibz_k_points(self): 

470 return self.kpts_coord.copy() 

471 

472 def get_number_of_spins(self): 

473 return self.nspin 

474 

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

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

477 

478 def get_fermi_levels(self): 

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

480 

481 def get_fermi_level(self): 

482 return max(self.get_fermi_levels()) 

483 

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 

489 

490 

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 

499 

500 def set_positions(self, mmpositions): 

501 self.mmpositions = mmpositions 

502 

503 def set_charges(self, mmcharges): 

504 self.mmcharges = mmcharges 

505 

506 def write_mmcharges(self, filename): 

507 """ mok all 

508 write external charges as monopoles for dftb+. 

509 

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

519 

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) 

526 

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

532 

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 

546 

547 

548def read_max_angular_momentum(path): 

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

550 

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 

564 

565 # Sometimes there ar commas, sometimes not: 

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

567 

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