Coverage for /builds/alexhroom/ase/ase/calculators/dftb.py: 75.98%

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

11import numpy as np

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

14 kpts2sizeandoffsets, BadConfiguration)

15from ase.units import Bohr, Hartree

18class Dftb(FileIOCalculator):

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

20 'stress', 'dipole']

21 discard_results_on_any_change = True

23 fileio_rules = FileIOCalculator.ruleset(

24 configspec=dict(skt_path=None),

25 stdout_name='{prefix}.out')

27 def __init__(self, restart=None,

28 ignore_bad_restart_file=FileIOCalculator._deprecated,

29 label='dftb', atoms=None, kpts=None,

30 slako_dir=None,

31 command=None,

32 profile=None,

33 **kwargs):

34 """

35 All keywords for the dftb_in.hsd input file (see the DFTB+ manual)

36 can be set by ASE. Consider the following input file block::

38 Hamiltonian = DFTB {

39 SCC = Yes

40 SCCTolerance = 1e-8

41 MaxAngularMomentum = {

42 H = s

43 O = p

44 }

45 }

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

48 following settings:

50 >>> from ase.calculators.dftb import Dftb

51 >>>

52 >>> calc = Dftb(Hamiltonian_='DFTB', # line is included by default

53 ... Hamiltonian_SCC='Yes',

54 ... Hamiltonian_SCCTolerance=1e-8,

55 ... Hamiltonian_MaxAngularMomentum_='',

56 ... Hamiltonian_MaxAngularMomentum_H='s',

57 ... Hamiltonian_MaxAngularMomentum_O='p')

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

60 arguments can be used:

62 restart: str

63 Prefix for restart file. May contain a directory.

64 Default is None: don't restart.

65 ignore_bad_restart_file: bool

66 Ignore broken or missing restart file. By default, it is an

67 error if the restart file is missing or broken.

68 label: str (default 'dftb')

69 Prefix used for the main output file (<label>.out).

70 atoms: Atoms object (default None)

71 Optional Atoms object to which the calculator will be

72 attached. When restarting, atoms will get its positions and

73 unit-cell updated from file.

74 kpts: (default None)

75 Brillouin zone sampling:

77 * ``(1,1,1)`` or ``None``: Gamma-point only

78 * ``(n1,n2,n3)``: Monkhorst-Pack grid

79 * ``dict``: Interpreted as a path in the Brillouin zone if

80 it contains the 'path_' keyword. Otherwise it is converted

81 into a Monkhorst-Pack grid using

82 ``ase.calculators.calculator.kpts2sizeandoffsets``

83 * ``[(k11,k12,k13),(k21,k22,k23),...]``: Explicit (Nkpts x 3)

84 array of k-points in units of the reciprocal lattice vectors

85 (each with equal weight)

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

89 pcpot: PointCharge object

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

91 """

93 if command is None:

94 if 'DFTB_COMMAND' in self.cfg:

95 command = self.cfg['DFTB_COMMAND'] + ' > PREFIX.out'

97 if command is None and profile is None:

98 try:

99 profile = self.load_argv_profile(self.cfg, 'dftb')

100 except BadConfiguration:

101 pass

102

103 if command is None:

104 command = 'dftb+ > PREFIX.out'

105

106 if slako_dir is None:

107 if profile is not None:

108 slako_dir = profile.configvars.get('skt_path')

109

110 if slako_dir is None:

111 slako_dir = self.cfg.get('DFTB_PREFIX', './')

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

113 slako_dir += '/'

114

115 self.slako_dir = slako_dir

116 if kwargs.get('Hamiltonian_', 'DFTB') == 'DFTB':

117 self.default_parameters = dict(

118 Hamiltonian_='DFTB',

119 Hamiltonian_SlaterKosterFiles_='Type2FileNames',

120 Hamiltonian_SlaterKosterFiles_Prefix=self.slako_dir,

121 Hamiltonian_SlaterKosterFiles_Separator='"-"',

122 Hamiltonian_SlaterKosterFiles_Suffix='".skf"',

123 Hamiltonian_MaxAngularMomentum_='',

124 Options_='',

125 Options_WriteResultsTag='Yes',

126 ParserOptions_='',

127 ParserOptions_ParserVersion=1,

128 ParserOptions_IgnoreUnprocessedNodes='Yes')

129 else:

130 self.default_parameters = dict(

131 Options_='',

132 Options_WriteResultsTag='Yes',

133 ParserOptions_='',

134 ParserOptions_ParserVersion=1,

135 ParserOptions_IgnoreUnprocessedNodes='Yes')

136

137 self.pcpot = None

138 self.lines = None

139 self.atoms = None

140 self.atoms_input = None

141 self.do_forces = False

142

143 super().__init__(restart, ignore_bad_restart_file,

144 label, atoms, command=command,

145 profile=profile, **kwargs)

146

147 # Determine number of spin channels

148 try:

149 entry = kwargs['Hamiltonian_SpinPolarisation']

150 spinpol = 'colinear' in entry.lower()

151 except KeyError:

152 spinpol = False

153 self.nspin = 2 if spinpol else 1

154

155 # kpoint stuff by ase

156 self.kpts = kpts

157 self.kpts_coord = None

158

159 if self.kpts is not None:

160 initkey = 'Hamiltonian_KPointsAndWeights'

161 mp_mesh = None

162 offsets = None

163

164 if isinstance(self.kpts, dict):

165 if 'path' in self.kpts:

166 # kpts is path in Brillouin zone

167 self.parameters[initkey + '_'] = 'Klines '

168 self.kpts_coord = kpts2ndarray(self.kpts, atoms=atoms)

169 else:

170 # kpts is (implicit) definition of

171 # Monkhorst-Pack grid

172 self.parameters[initkey + '_'] = 'SupercellFolding '

173 mp_mesh, offsets = kpts2sizeandoffsets(atoms=atoms,

174 **self.kpts)

175 elif np.array(self.kpts).ndim == 1:

176 # kpts is Monkhorst-Pack grid

177 self.parameters[initkey + '_'] = 'SupercellFolding '

178 mp_mesh = self.kpts

179 offsets = [0.] * 3

180 elif np.array(self.kpts).ndim == 2:

181 # kpts is (N x 3) list/array of k-point coordinates

182 # each will be given equal weight

183 self.parameters[initkey + '_'] = ''

184 self.kpts_coord = np.array(self.kpts)

185 else:

186 raise ValueError('Illegal kpts definition:' + str(self.kpts))

187

188 if mp_mesh is not None:

189 eps = 1e-10

190 for i in range(3):

191 key = initkey + '_empty%03d' % i

192 val = [mp_mesh[i] if j == i else 0 for j in range(3)]

193 self.parameters[key] = ' '.join(map(str, val))

194 offsets[i] *= mp_mesh[i]

195 assert abs(offsets[i]) < eps or abs(offsets[i] - 0.5) < eps

196 # DFTB+ uses a different offset convention, where

197 # the k-point mesh is already Gamma-centered prior

198 # to the addition of any offsets

199 if mp_mesh[i] % 2 == 0:

200 offsets[i] += 0.5

201 key = initkey + '_empty%03d' % 3

202 self.parameters[key] = ' '.join(map(str, offsets))

203

204 elif self.kpts_coord is not None:

205 for i, c in enumerate(self.kpts_coord):

206 key = initkey + '_empty%09d' % i

207 c_str = ' '.join(map(str, c))

208 if 'Klines' in self.parameters[initkey + '_']:

209 c_str = '1 ' + c_str

210 else:

211 c_str += ' 1.0'

212 self.parameters[key] = c_str

213

214 def write_dftb_in(self, outfile):

215 """ Write the innput file for the dftb+ calculation.

216 Geometry is taken always from the file 'geo_end.gen'.

217 """

218

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

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

221 outfile.write('} \n')

222 outfile.write(' \n')

223

224 params = self.parameters.copy()

225

226 s = 'Hamiltonian_MaxAngularMomentum_'

227 for key in params:

228 if key.startswith(s) and len(key) > len(s):

229 break

230 else:

231 if params.get('Hamiltonian_', 'DFTB') == 'DFTB':

232 # User didn't specify max angular mometa. Get them from

233 # the .skf files:

234 symbols = set(self.atoms.get_chemical_symbols())

235 for symbol in symbols:

236 path = os.path.join(self.slako_dir,

237 '{0}-{0}.skf'.format(symbol))

238 l = read_max_angular_momentum(path)

239 params[s + symbol] = '"{}"'.format('spdf'[l])

240

241 if self.do_forces:

242 params['Analysis_'] = ''

243 params['Analysis_CalculateForces'] = 'Yes'

244

245 # --------MAIN KEYWORDS-------

246 previous_key = 'dummy_'

247 myspace = ' '

248 for key, value in sorted(params.items()):

249 current_depth = key.rstrip('_').count('_')

250 previous_depth = previous_key.rstrip('_').count('_')

251 for my_backslash in reversed(

252 range(previous_depth - current_depth)):

253 outfile.write(3 * (1 + my_backslash) * myspace + '} \n')

254 outfile.write(3 * current_depth * myspace)

255 if key.endswith('_') and len(value) > 0:

256 outfile.write(key.rstrip('_').rsplit('_')[-1] +

257 ' = ' + str(value) + '{ \n')

258 elif (key.endswith('_') and (len(value) == 0)

259 and current_depth == 0): # E.g. 'Options {'

260 outfile.write(key.rstrip('_').rsplit('_')[-1] +

261 ' ' + str(value) + '{ \n')

262 elif (key.endswith('_') and (len(value) == 0)

263 and current_depth > 0): # E.g. 'Hamiltonian_Max... = {'

264 outfile.write(key.rstrip('_').rsplit('_')[-1] +

265 ' = ' + str(value) + '{ \n')

266 elif key.count('_empty') == 1:

267 outfile.write(str(value) + ' \n')

268 elif ((key == 'Hamiltonian_ReadInitialCharges') and

269 (str(value).upper() == 'YES')):

270 f1 = os.path.isfile(self.directory + os.sep + 'charges.dat')

271 f2 = os.path.isfile(self.directory + os.sep + 'charges.bin')

272 if not (f1 or f2):

273 print('charges.dat or .bin not found, switching off guess')

274 value = 'No'

275 outfile.write(key.rsplit('_')[-1] + ' = ' + str(value) + ' \n')

276 else:

277 outfile.write(key.rsplit('_')[-1] + ' = ' + str(value) + ' \n')

278 if self.pcpot is not None and ('DFTB' in str(value)):

279 outfile.write(' ElectricField = { \n')

280 outfile.write(' PointCharges = { \n')

281 outfile.write(

282 ' CoordsAndCharges [Angstrom] = DirectRead { \n')

283 outfile.write(' Records = ' +

284 str(len(self.pcpot.mmcharges)) + ' \n')

285 outfile.write(

286 ' File = "dftb_external_charges.dat" \n')

287 outfile.write(' } \n')

288 outfile.write(' } \n')

289 outfile.write(' } \n')

290 previous_key = key

291 current_depth = key.rstrip('_').count('_')

292 for my_backslash in reversed(range(current_depth)):

293 outfile.write(3 * my_backslash * myspace + '} \n')

294

295 def check_state(self, atoms):

296 system_changes = FileIOCalculator.check_state(self, atoms)

297 # Ignore unit cell for molecules:

298 if not atoms.pbc.any() and 'cell' in system_changes:

299 system_changes.remove('cell')

300 if self.pcpot and self.pcpot.mmpositions is not None:

301 system_changes.append('positions')

302 return system_changes

303

304 def write_input(self, atoms, properties=None, system_changes=None):

305 from ase.io import write

306 if properties is not None:

307 if 'forces' in properties or 'stress' in properties:

308 self.do_forces = True

309 FileIOCalculator.write_input(

310 self, atoms, properties, system_changes)

311 with open(os.path.join(self.directory, 'dftb_in.hsd'), 'w') as fd:

312 self.write_dftb_in(fd)

313 write(os.path.join(self.directory, 'geo_end.gen'), atoms,

314 parallel=False)

315 # self.atoms is none until results are read out,

316 # then it is set to the ones at writing input

317 self.atoms_input = atoms

318 self.atoms = None

319 if self.pcpot:

320 self.pcpot.write_mmcharges('dftb_external_charges.dat')

321

322 def read_results(self):

323 """ all results are read from results.tag file

324 It will be destroyed after it is read to avoid

325 reading it once again after some runtime error """

326

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

328 self.lines = fd.readlines()

329

330 self.atoms = self.atoms_input

331 charges, energy, dipole = self.read_charges_energy_dipole()

332 if charges is not None:

333 self.results['charges'] = charges

334 self.results['energy'] = energy

335 if dipole is not None:

336 self.results['dipole'] = dipole

337 if self.do_forces:

338 forces = self.read_forces()

339 self.results['forces'] = forces

340 self.mmpositions = None

341

342 # stress stuff begins

343 sstring = 'stress'

344 have_stress = False

345 stress = []

346 for iline, line in enumerate(self.lines):

347 if sstring in line:

348 have_stress = True

349 start = iline + 1

350 end = start + 3

351 for i in range(start, end):

352 cell = [float(x) for x in self.lines[i].split()]

353 stress.append(cell)

354 if have_stress:

355 stress = -np.array(stress) * Hartree / Bohr**3

356 self.results['stress'] = stress.flat[[0, 4, 8, 5, 2, 1]]

357 # stress stuff ends

358

359 # eigenvalues and fermi levels

360 fermi_levels = self.read_fermi_levels()

361 if fermi_levels is not None:

362 self.results['fermi_levels'] = fermi_levels

363

364 eigenvalues = self.read_eigenvalues()

365 if eigenvalues is not None:

366 self.results['eigenvalues'] = eigenvalues

367

368 # calculation was carried out with atoms written in write_input

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

370

371 def read_forces(self):

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

373 from ase.units import Bohr, Hartree

374

375 # Initialise the indices so their scope

376 # reaches outside of the for loop

377 index_force_begin = -1

378 index_force_end = -1

379

380 # Force line indexes

381 for iline, line in enumerate(self.lines):

382 fstring = 'forces '

383 if line.find(fstring) >= 0:

384 index_force_begin = iline + 1

385 line1 = line.replace(':', ',')

386 index_force_end = iline + 1 + \

387 int(line1.split(',')[-1])

388 break

389

390 gradients = []

391 for j in range(index_force_begin, index_force_end):

392 word = self.lines[j].split()

393 gradients.append([float(word[k]) for k in range(3)])

394

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

396

397 def read_charges_energy_dipole(self):

398 """Get partial charges on atoms

399 in case we cannot find charges they are set to None

400 """

401 with open(os.path.join(self.directory, 'detailed.out')) as fd:

402 lines = fd.readlines()

403

404 for line in lines:

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

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

407 break

408

409 qm_charges = []

410 for n, line in enumerate(lines):

411 if ('Atom' and 'Charge' in line):

412 chargestart = n + 1

413 break

414 else:

415 # print('Warning: did not find DFTB-charges')

416 # print('This is ok if flag SCC=No')

417 return None, energy, None

418

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

420 for line in lines1:

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

422

423 dipole = None

424 for line in lines:

425 if 'Dipole moment:' in line and 'au' in line:

426 line = line.replace("Dipole moment:", "").replace("au", "")

427 dipole = np.array(line.split(), dtype=float) * Bohr

428

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

430

431 def get_charges(self, atoms):

432 """ Get the calculated charges

433 this is inhereted to atoms object """

434 if 'charges' in self.results:

435 return self.results['charges']

436 else:

437 return None

438

439 def read_eigenvalues(self):

440 """ Read Eigenvalues from dftb output file (results.tag).

441 Unfortunately, the order seems to be scrambled. """

442 # Eigenvalue line indexes

443 index_eig_begin = None

444 for iline, line in enumerate(self.lines):

445 fstring = 'eigenvalues '

446 if line.find(fstring) >= 0:

447 index_eig_begin = iline + 1

448 line1 = line.replace(':', ',')

449 ncol, nband, nkpt, nspin = map(int, line1.split(',')[-4:])

450 break

451 else:

452 return None

453

454 # Take into account that the last row may lack

455 # columns if nkpt * nspin * nband % ncol != 0

456 nrow = int(np.ceil(nkpt * nspin * nband * 1. / ncol))

457 index_eig_end = index_eig_begin + nrow

458 ncol_last = len(self.lines[index_eig_end - 1].split())

459 # XXX dirty fix

460 self.lines[index_eig_end - 1] = (

461 self.lines[index_eig_end - 1].strip()

462 + ' 0.0 ' * (ncol - ncol_last))

463

464 eig = np.loadtxt(self.lines[index_eig_begin:index_eig_end]).flatten()

465 eig *= Hartree

466 N = nkpt * nband

467 eigenvalues = [eig[i * N:(i + 1) * N].reshape((nkpt, nband))

468 for i in range(nspin)]

469

470 return eigenvalues

471

472 def read_fermi_levels(self):

473 """ Read Fermi level(s) from dftb output file (results.tag). """

474 # Fermi level line indexes

475 for iline, line in enumerate(self.lines):

476 fstring = 'fermi_level '

477 if line.find(fstring) >= 0:

478 index_fermi = iline + 1

479 break

480 else:

481 return None

482

483 fermi_levels = []

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

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

486

487 for word in words:

488 e = float(word)

489 # In non-spin-polarized calculations with DFTB+ v17.1,

490 # two Fermi levels are given, with the second one being 0,

491 # but we don't want to add that one to the list

492 if abs(e) > 1e-8:

493 fermi_levels.append(e)

494

495 return np.array(fermi_levels) * Hartree

496

497 def get_ibz_k_points(self):

498 return self.kpts_coord.copy()

499

500 def get_number_of_spins(self):

501 return self.nspin

502

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

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

505

506 def get_fermi_levels(self):

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

508

509 def get_fermi_level(self):

510 return max(self.get_fermi_levels())

511

512 def embed(self, mmcharges=None, directory='./'):

513 """Embed atoms in point-charges (mmcharges)

514 """

515 self.pcpot = PointChargePotential(mmcharges, self.directory)

516 return self.pcpot

517

518

519class PointChargePotential:

520 def __init__(self, mmcharges, directory='./'):

521 """Point-charge potential for DFTB+.

522 """

523 self.mmcharges = mmcharges

524 self.directory = directory

525 self.mmpositions = None

526 self.mmforces = None

527

528 def set_positions(self, mmpositions):

529 self.mmpositions = mmpositions

530

531 def set_charges(self, mmcharges):

532 self.mmcharges = mmcharges

533

534 def write_mmcharges(self, filename):

535 """ mok all

536 write external charges as monopoles for dftb+.

537

538 """

539 if self.mmcharges is None:

540 print("DFTB: Warning: not writing exernal charges ")

541 return

542 with open(os.path.join(self.directory, filename), 'w') as charge_file:

543 for [pos, charge] in zip(self.mmpositions, self.mmcharges):

544 [x, y, z] = pos

545 charge_file.write('%12.6f %12.6f %12.6f %12.6f \n'

546 % (x, y, z, charge))

547

548 def get_forces(self, calc, get_forces=True):

549 """ returns forces on point charges if the flag get_forces=True """

550 if get_forces:

551 return self.read_forces_on_pointcharges()

552 else:

553 return np.zeros_like(self.mmpositions)

554

555 def read_forces_on_pointcharges(self):

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

557 from ase.units import Bohr, Hartree

558 with open(os.path.join(self.directory, 'detailed.out')) as fd:

559 lines = fd.readlines()

560

561 external_forces = []

562 for n, line in enumerate(lines):

563 if ('Forces on external charges' in line):

564 chargestart = n + 1

565 break

566 else:

567 raise RuntimeError(

568 'Problem in reading forces on MM external-charges')

569 lines1 = lines[chargestart:(chargestart + len(self.mmcharges))]

570 for line in lines1:

571 external_forces.append(

572 [float(i) for i in line.split()])

573 return np.array(external_forces) * Hartree / Bohr

574

575

576def read_max_angular_momentum(path):

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

578

579 See dftb.org for A detailed description of the Slater-Koster file format.

580 """

581 with open(path) as fd:

582 line = fd.readline()

583 if line[0] == '@':

584 # Extended format

585 fd.readline()

586 l = 3

587 pos = 9

588 else:

589 # Simple format:

590 l = 2

591 pos = 7

592

593 # Sometimes there ar commas, sometimes not:

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

595

596 occs = [float(f) for f in line.split()[pos:pos + l + 1]]

597 for f in occs:

598 if f > 0.0:

599 return l

600 l -= 1