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

1# fmt: off

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

5http://www.dftbplus.org/

6http://www.dftb.org/

8Initial development: markus.kaukonen@iki.fi

9"""

11import os

13import numpy as np

15from ase.calculators.calculator import (

16 BadConfiguration,

17 FileIOCalculator,

18 kpts2ndarray,

19 kpts2sizeandoffsets,

20)

21from ase.units import Bohr, Hartree

24class Dftb(FileIOCalculator):

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

26 'stress', 'dipole']

27 discard_results_on_any_change = True

29 fileio_rules = FileIOCalculator.ruleset(

30 configspec=dict(skt_path=None),

31 stdout_name='{prefix}.out')

33 def __init__(self, restart=None,

34 ignore_bad_restart_file=FileIOCalculator._deprecated,

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

36 slako_dir=None,

37 command=None,

38 profile=None,

39 **kwargs):

40 """

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

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

44 Hamiltonian = DFTB {

45 SCC = Yes

46 SCCTolerance = 1e-8

47 MaxAngularMomentum = {

48 H = s

49 O = p

50 }

51 }

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

54 following settings:

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

57 >>>

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

59 ... Hamiltonian_SCC='Yes',

60 ... Hamiltonian_SCCTolerance=1e-8,

61 ... Hamiltonian_MaxAngularMomentum_='',

62 ... Hamiltonian_MaxAngularMomentum_H='s',

63 ... Hamiltonian_MaxAngularMomentum_O='p')

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

66 arguments can be used:

68 restart: str

69 Prefix for restart file. May contain a directory.

70 Default is None: don't restart.

71 ignore_bad_restart_file: bool

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

73 error if the restart file is missing or broken.

74 label: str (default 'dftb')

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

76 atoms: Atoms object (default None)

77 Optional Atoms object to which the calculator will be

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

79 unit-cell updated from file.

80 kpts: (default None)

81 Brillouin zone sampling:

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

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

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

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

87 into a Monkhorst-Pack grid using

88 ``ase.calculators.calculator.kpts2sizeandoffsets``

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

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

91 (each with equal weight)

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

95 pcpot: PointCharge object

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

97 """

99 if command is None:

100 if 'DFTB_COMMAND' in self.cfg:

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

102

103 if command is None and profile is None:

104 try:

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

106 except BadConfiguration:

107 pass

108

109 if command is None:

110 command = 'dftb+ > PREFIX.out'

111

112 if slako_dir is None:

113 if profile is not None:

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

115

116 if slako_dir is None:

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

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

119 slako_dir += '/'

120

121 self.slako_dir = slako_dir

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

123 self.default_parameters = dict(

124 Hamiltonian_='DFTB',

125 Hamiltonian_SlaterKosterFiles_='Type2FileNames',

126 Hamiltonian_SlaterKosterFiles_Prefix=self.slako_dir,

127 Hamiltonian_SlaterKosterFiles_Separator='"-"',

128 Hamiltonian_SlaterKosterFiles_Suffix='".skf"',

129 Hamiltonian_MaxAngularMomentum_='',

130 Options_='',

131 Options_WriteResultsTag='Yes',

132 ParserOptions_='',

133 ParserOptions_ParserVersion=1,

134 ParserOptions_IgnoreUnprocessedNodes='Yes')

135 else:

136 self.default_parameters = dict(

137 Options_='',

138 Options_WriteResultsTag='Yes',

139 ParserOptions_='',

140 ParserOptions_ParserVersion=1,

141 ParserOptions_IgnoreUnprocessedNodes='Yes')

142

143 self.pcpot = None

144 self.lines = None

145 self.atoms = None

146 self.atoms_input = None

147 self.do_forces = False

148

149 super().__init__(restart, ignore_bad_restart_file,

150 label, atoms, command=command,

151 profile=profile, **kwargs)

152

153 # Determine number of spin channels

154 try:

155 entry = kwargs['Hamiltonian_SpinPolarisation']

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

157 except KeyError:

158 spinpol = False

159 self.nspin = 2 if spinpol else 1

160

161 # kpoint stuff by ase

162 self.kpts = kpts

163 self.kpts_coord = None

164

165 if self.kpts is not None:

166 initkey = 'Hamiltonian_KPointsAndWeights'

167 mp_mesh = None

168 offsets = None

169

170 if isinstance(self.kpts, dict):

171 if 'path' in self.kpts:

172 # kpts is path in Brillouin zone

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

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

175 else:

176 # kpts is (implicit) definition of

177 # Monkhorst-Pack grid

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

179 mp_mesh, offsets = kpts2sizeandoffsets(atoms=atoms,

180 **self.kpts)

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

182 # kpts is Monkhorst-Pack grid

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

184 mp_mesh = self.kpts

185 offsets = [0.] * 3

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

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

188 # each will be given equal weight

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

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

191 else:

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

193

194 if mp_mesh is not None:

195 eps = 1e-10

196 for i in range(3):

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

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

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

200 offsets[i] *= mp_mesh[i]

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

202 # DFTB+ uses a different offset convention, where

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

204 # to the addition of any offsets

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

206 offsets[i] += 0.5

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

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

209

210 elif self.kpts_coord is not None:

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

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

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

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

215 c_str = '1 ' + c_str

216 else:

217 c_str += ' 1.0'

218 self.parameters[key] = c_str

219

220 def write_dftb_in(self, outfile):

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

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

223 """

224

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

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

227 outfile.write('} \n')

228 outfile.write(' \n')

229

230 params = self.parameters.copy()

231

232 s = 'Hamiltonian_MaxAngularMomentum_'

233 for key in params:

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

235 break

236 else:

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

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

239 # the .skf files:

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

241 for symbol in symbols:

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

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

244 l = read_max_angular_momentum(path)

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

246

247 if self.do_forces:

248 params['Analysis_'] = ''

249 params['Analysis_CalculateForces'] = 'Yes'

250

251 # --------MAIN KEYWORDS-------

252 previous_key = 'dummy_'

253 myspace = ' '

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

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

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

257 for my_backslash in reversed(

258 range(previous_depth - current_depth)):

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

260 outfile.write(3 * current_depth * myspace)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

278 if not (f1 or f2):

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

280 value = 'No'

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

282 else:

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

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

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

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

287 outfile.write(

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

289 outfile.write(' Records = ' +

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

291 outfile.write(

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

293 outfile.write(' } \n')

294 outfile.write(' } \n')

295 outfile.write(' } \n')

296 previous_key = key

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

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

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

300

301 def check_state(self, atoms):

302 system_changes = FileIOCalculator.check_state(self, atoms)

303 # Ignore unit cell for molecules:

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

305 system_changes.remove('cell')

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

307 system_changes.append('positions')

308 return system_changes

309

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

311 from ase.io import write

312 if properties is not None:

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

314 self.do_forces = True

315 FileIOCalculator.write_input(

316 self, atoms, properties, system_changes)

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

318 self.write_dftb_in(fd)

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

320 parallel=False)

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

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

323 self.atoms_input = atoms

324 self.atoms = None

325 if self.pcpot:

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

327

328 def read_results(self):

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

330 It will be destroyed after it is read to avoid

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

332

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

334 self.lines = fd.readlines()

335

336 self.atoms = self.atoms_input

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

338 if charges is not None:

339 self.results['charges'] = charges

340 self.results['energy'] = energy

341 if dipole is not None:

342 self.results['dipole'] = dipole

343 if self.do_forces:

344 forces = self.read_forces()

345 self.results['forces'] = forces

346 self.mmpositions = None

347

348 # stress stuff begins

349 sstring = 'stress'

350 have_stress = False

351 stress = []

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

353 if sstring in line:

354 have_stress = True

355 start = iline + 1

356 end = start + 3

357 for i in range(start, end):

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

359 stress.append(cell)

360 if have_stress:

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

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

363 # stress stuff ends

364

365 # eigenvalues and fermi levels

366 fermi_levels = self.read_fermi_levels()

367 if fermi_levels is not None:

368 self.results['fermi_levels'] = fermi_levels

369

370 eigenvalues = self.read_eigenvalues()

371 if eigenvalues is not None:

372 self.results['eigenvalues'] = eigenvalues

373

374 # calculation was carried out with atoms written in write_input

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

376

377 def read_forces(self):

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

379 from ase.units import Bohr, Hartree

380

381 # Initialise the indices so their scope

382 # reaches outside of the for loop

383 index_force_begin = -1

384 index_force_end = -1

385

386 # Force line indexes

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

388 fstring = 'forces '

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

390 index_force_begin = iline + 1

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

392 index_force_end = iline + 1 + \

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

394 break

395

396 gradients = []

397 for j in range(index_force_begin, index_force_end):

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

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

400

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

402

403 def read_charges_energy_dipole(self):

404 """Get partial charges on atoms

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

406 """

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

408 lines = fd.readlines()

409

410 for line in lines:

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

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

413 break

414

415 qm_charges = []

416 for n, line in enumerate(lines):

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

418 chargestart = n + 1

419 break

420 else:

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

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

423 return None, energy, None

424

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

426 for line in lines1:

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

428

429 dipole = None

430 for line in lines:

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

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

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

434

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

436

437 def get_charges(self, atoms):

438 """ Get the calculated charges

439 this is inhereted to atoms object """

440 if 'charges' in self.results:

441 return self.results['charges']

442 else:

443 return None

444

445 def read_eigenvalues(self):

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

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

448 # Eigenvalue line indexes

449 index_eig_begin = None

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

451 fstring = 'eigenvalues '

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

453 index_eig_begin = iline + 1

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

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

456 break

457 else:

458 return None

459

460 # Take into account that the last row may lack

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

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

463 index_eig_end = index_eig_begin + nrow

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

465 # XXX dirty fix

466 self.lines[index_eig_end - 1] = (

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

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

469

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

471 eig *= Hartree

472 N = nkpt * nband

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

474 for i in range(nspin)]

475

476 return eigenvalues

477

478 def read_fermi_levels(self):

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

480 # Fermi level line indexes

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

482 fstring = 'fermi_level '

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

484 index_fermi = iline + 1

485 break

486 else:

487 return None

488

489 fermi_levels = []

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

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

492

493 for word in words:

494 e = float(word)

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

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

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

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

499 fermi_levels.append(e)

500

501 return np.array(fermi_levels) * Hartree

502

503 def get_ibz_k_points(self):

504 return self.kpts_coord.copy()

505

506 def get_number_of_spins(self):

507 return self.nspin

508

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

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

511

512 def get_fermi_levels(self):

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

514

515 def get_fermi_level(self):

516 return max(self.get_fermi_levels())

517

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

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

520 """

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

522 return self.pcpot

523

524

525class PointChargePotential:

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

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

528 """

529 self.mmcharges = mmcharges

530 self.directory = directory

531 self.mmpositions = None

532 self.mmforces = None

533

534 def set_positions(self, mmpositions):

535 self.mmpositions = mmpositions

536

537 def set_charges(self, mmcharges):

538 self.mmcharges = mmcharges

539

540 def write_mmcharges(self, filename):

541 """ mok all

542 write external charges as monopoles for dftb+.

543

544 """

545 if self.mmcharges is None:

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

547 return

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

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

550 [x, y, z] = pos

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

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

553

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

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

556 if get_forces:

557 return self.read_forces_on_pointcharges()

558 else:

559 return np.zeros_like(self.mmpositions)

560

561 def read_forces_on_pointcharges(self):

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

563 from ase.units import Bohr, Hartree

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

565 lines = fd.readlines()

566

567 external_forces = []

568 for n, line in enumerate(lines):

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

570 chargestart = n + 1

571 break

572 else:

573 raise RuntimeError(

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

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

576 for line in lines1:

577 external_forces.append(

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

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

580

581

582def read_max_angular_momentum(path):

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

584

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

586 """

587 with open(path) as fd:

588 line = fd.readline()

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

590 # Extended format

591 fd.readline()

592 l = 3

593 pos = 9

594 else:

595 # Simple format:

596 l = 2

597 pos = 7

598

599 # Sometimes there ar commas, sometimes not:

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

601

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

603 for f in occs:

604 if f > 0.0:

605 return l

606 l -= 1