Coverage for /builds/hweiske/ase/ase/calculators/dftb.py: 75.00%

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)

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 stdout_name='{prefix}.out')

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 command=None,

31 profile=None,

32 **kwargs):

33 """

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

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

37 Hamiltonian = DFTB {

38 SCC = Yes

39 SCCTolerance = 1e-8

40 MaxAngularMomentum = {

41 H = s

42 O = p

43 }

44 }

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

47 following settings:

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

50 >>>

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

52 ... Hamiltonian_SCC='Yes',

53 ... Hamiltonian_SCCTolerance=1e-8,

54 ... Hamiltonian_MaxAngularMomentum_='',

55 ... Hamiltonian_MaxAngularMomentum_H='s',

56 ... Hamiltonian_MaxAngularMomentum_O='p')

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

59 arguments can be used:

61 restart: str

62 Prefix for restart file. May contain a directory.

63 Default is None: don't restart.

64 ignore_bad_restart_file: bool

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

66 error if the restart file is missing or broken.

67 label: str (default 'dftb')

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

69 atoms: Atoms object (default None)

70 Optional Atoms object to which the calculator will be

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

72 unit-cell updated from file.

73 kpts: (default None)

74 Brillouin zone sampling:

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

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

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

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

80 into a Monkhorst-Pack grid using

81 ``ase.calculators.calculator.kpts2sizeandoffsets``

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

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

84 (each with equal weight)

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

88 pcpot: PointCharge object

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

90 """

92 if command is None:

93 if 'DFTB_COMMAND' in self.cfg:

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

95 else:

96 command = 'dftb+ > PREFIX.out'

98 if slako_dir is None:

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

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

101 slako_dir += '/'

102

103 self.slako_dir = slako_dir

104

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

106 self.default_parameters = dict(

107 Hamiltonian_='DFTB',

108 Hamiltonian_SlaterKosterFiles_='Type2FileNames',

109 Hamiltonian_SlaterKosterFiles_Prefix=self.slako_dir,

110 Hamiltonian_SlaterKosterFiles_Separator='"-"',

111 Hamiltonian_SlaterKosterFiles_Suffix='".skf"',

112 Hamiltonian_MaxAngularMomentum_='',

113 Options_='',

114 Options_WriteResultsTag='Yes',

115 ParserOptions_='',

116 ParserOptions_ParserVersion=1,

117 ParserOptions_IgnoreUnprocessedNodes='Yes')

118 else:

119 self.default_parameters = dict(

120 Options_='',

121 Options_WriteResultsTag='Yes',

122 ParserOptions_='',

123 ParserOptions_ParserVersion=1,

124 ParserOptions_IgnoreUnprocessedNodes='Yes')

125

126 self.pcpot = None

127 self.lines = None

128 self.atoms = None

129 self.atoms_input = None

130 self.do_forces = False

131

132 super().__init__(restart, ignore_bad_restart_file,

133 label, atoms, command=command,

134 profile=profile, **kwargs)

135

136 # Determine number of spin channels

137 try:

138 entry = kwargs['Hamiltonian_SpinPolarisation']

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

140 except KeyError:

141 spinpol = False

142 self.nspin = 2 if spinpol else 1

143

144 # kpoint stuff by ase

145 self.kpts = kpts

146 self.kpts_coord = None

147

148 if self.kpts is not None:

149 initkey = 'Hamiltonian_KPointsAndWeights'

150 mp_mesh = None

151 offsets = None

152

153 if isinstance(self.kpts, dict):

154 if 'path' in self.kpts:

155 # kpts is path in Brillouin zone

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

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

158 else:

159 # kpts is (implicit) definition of

160 # Monkhorst-Pack grid

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

162 mp_mesh, offsets = kpts2sizeandoffsets(atoms=atoms,

163 **self.kpts)

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

165 # kpts is Monkhorst-Pack grid

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

167 mp_mesh = self.kpts

168 offsets = [0.] * 3

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

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

171 # each will be given equal weight

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

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

174 else:

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

176

177 if mp_mesh is not None:

178 eps = 1e-10

179 for i in range(3):

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

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

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

183 offsets[i] *= mp_mesh[i]

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

185 # DFTB+ uses a different offset convention, where

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

187 # to the addition of any offsets

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

189 offsets[i] += 0.5

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

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

192

193 elif self.kpts_coord is not None:

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

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

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

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

198 c_str = '1 ' + c_str

199 else:

200 c_str += ' 1.0'

201 self.parameters[key] = c_str

202

203 def write_dftb_in(self, outfile):

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

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

206 """

207

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

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

210 outfile.write('} \n')

211 outfile.write(' \n')

212

213 params = self.parameters.copy()

214

215 s = 'Hamiltonian_MaxAngularMomentum_'

216 for key in params:

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

218 break

219 else:

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

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

222 # the .skf files:

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

224 for symbol in symbols:

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

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

227 l = read_max_angular_momentum(path)

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

229

230 if self.do_forces:

231 params['Analysis_'] = ''

232 params['Analysis_CalculateForces'] = 'Yes'

233

234 # --------MAIN KEYWORDS-------

235 previous_key = 'dummy_'

236 myspace = ' '

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

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

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

240 for my_backslash in reversed(

241 range(previous_depth - current_depth)):

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

243 outfile.write(3 * current_depth * myspace)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

261 if not (f1 or f2):

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

263 value = 'No'

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

265 else:

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

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

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

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

270 outfile.write(

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

272 outfile.write(' Records = ' +

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

274 outfile.write(

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

276 outfile.write(' } \n')

277 outfile.write(' } \n')

278 outfile.write(' } \n')

279 previous_key = key

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

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

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

283

284 def check_state(self, atoms):

285 system_changes = FileIOCalculator.check_state(self, atoms)

286 # Ignore unit cell for molecules:

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

288 system_changes.remove('cell')

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

290 system_changes.append('positions')

291 return system_changes

292

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

294 from ase.io import write

295 if properties is not None:

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

297 self.do_forces = True

298 FileIOCalculator.write_input(

299 self, atoms, properties, system_changes)

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

301 self.write_dftb_in(fd)

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

303 parallel=False)

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

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

306 self.atoms_input = atoms

307 self.atoms = None

308 if self.pcpot:

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

310

311 def read_results(self):

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

313 It will be destroyed after it is read to avoid

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

315

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

317 self.lines = fd.readlines()

318

319 self.atoms = self.atoms_input

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

321 if charges is not None:

322 self.results['charges'] = charges

323 self.results['energy'] = energy

324 if dipole is not None:

325 self.results['dipole'] = dipole

326 if self.do_forces:

327 forces = self.read_forces()

328 self.results['forces'] = forces

329 self.mmpositions = None

330

331 # stress stuff begins

332 sstring = 'stress'

333 have_stress = False

334 stress = []

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

336 if sstring in line:

337 have_stress = True

338 start = iline + 1

339 end = start + 3

340 for i in range(start, end):

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

342 stress.append(cell)

343 if have_stress:

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

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

346 # stress stuff ends

347

348 # eigenvalues and fermi levels

349 fermi_levels = self.read_fermi_levels()

350 if fermi_levels is not None:

351 self.results['fermi_levels'] = fermi_levels

352

353 eigenvalues = self.read_eigenvalues()

354 if eigenvalues is not None:

355 self.results['eigenvalues'] = eigenvalues

356

357 # calculation was carried out with atoms written in write_input

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

359

360 def read_forces(self):

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

362 from ase.units import Bohr, Hartree

363

364 # Initialise the indices so their scope

365 # reaches outside of the for loop

366 index_force_begin = -1

367 index_force_end = -1

368

369 # Force line indexes

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

371 fstring = 'forces '

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

373 index_force_begin = iline + 1

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

375 index_force_end = iline + 1 + \

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

377 break

378

379 gradients = []

380 for j in range(index_force_begin, index_force_end):

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

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

383

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

385

386 def read_charges_energy_dipole(self):

387 """Get partial charges on atoms

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

389 """

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

391 lines = fd.readlines()

392

393 for line in lines:

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

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

396 break

397

398 qm_charges = []

399 for n, line in enumerate(lines):

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

401 chargestart = n + 1

402 break

403 else:

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

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

406 return None, energy, None

407

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

409 for line in lines1:

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

411

412 dipole = None

413 for line in lines:

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

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

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

417

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

419

420 def get_charges(self, atoms):

421 """ Get the calculated charges

422 this is inhereted to atoms object """

423 if 'charges' in self.results:

424 return self.results['charges']

425 else:

426 return None

427

428 def read_eigenvalues(self):

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

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

431 # Eigenvalue line indexes

432 index_eig_begin = None

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

434 fstring = 'eigenvalues '

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

436 index_eig_begin = iline + 1

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

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

439 break

440 else:

441 return None

442

443 # Take into account that the last row may lack

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

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

446 index_eig_end = index_eig_begin + nrow

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

448 # XXX dirty fix

449 self.lines[index_eig_end - 1] = (

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

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

452

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

454 eig *= Hartree

455 N = nkpt * nband

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

457 for i in range(nspin)]

458

459 return eigenvalues

460

461 def read_fermi_levels(self):

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

463 # Fermi level line indexes

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

465 fstring = 'fermi_level '

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

467 index_fermi = iline + 1

468 break

469 else:

470 return None

471

472 fermi_levels = []

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

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

475

476 for word in words:

477 e = float(word)

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

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

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

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

482 fermi_levels.append(e)

483

484 return np.array(fermi_levels) * Hartree

485

486 def get_ibz_k_points(self):

487 return self.kpts_coord.copy()

488

489 def get_number_of_spins(self):

490 return self.nspin

491

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

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

494

495 def get_fermi_levels(self):

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

497

498 def get_fermi_level(self):

499 return max(self.get_fermi_levels())

500

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

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

503 """

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

505 return self.pcpot

506

507

508class PointChargePotential:

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

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

511 """

512 self.mmcharges = mmcharges

513 self.directory = directory

514 self.mmpositions = None

515 self.mmforces = None

516

517 def set_positions(self, mmpositions):

518 self.mmpositions = mmpositions

519

520 def set_charges(self, mmcharges):

521 self.mmcharges = mmcharges

522

523 def write_mmcharges(self, filename):

524 """ mok all

525 write external charges as monopoles for dftb+.

526

527 """

528 if self.mmcharges is None:

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

530 return

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

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

533 [x, y, z] = pos

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

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

536

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

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

539 if get_forces:

540 return self.read_forces_on_pointcharges()

541 else:

542 return np.zeros_like(self.mmpositions)

543

544 def read_forces_on_pointcharges(self):

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

546 from ase.units import Bohr, Hartree

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

548 lines = fd.readlines()

549

550 external_forces = []

551 for n, line in enumerate(lines):

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

553 chargestart = n + 1

554 break

555 else:

556 raise RuntimeError(

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

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

559 for line in lines1:

560 external_forces.append(

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

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

563

564

565def read_max_angular_momentum(path):

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

567

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

569 """

570 with open(path) as fd:

571 line = fd.readline()

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

573 # Extended format

574 fd.readline()

575 l = 3

576 pos = 9

577 else:

578 # Simple format:

579 l = 2

580 pos = 7

581

582 # Sometimes there ar commas, sometimes not:

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

584

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

586 for f in occs:

587 if f > 0.0:

588 return l

589 l -= 1