Coverage for /builds/debichem-team/python-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 (

14 BadConfiguration,

15 FileIOCalculator,

16 kpts2ndarray,

17 kpts2sizeandoffsets,

18)

19from ase.units import Bohr, Hartree

22class Dftb(FileIOCalculator):

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

24 'stress', 'dipole']

25 discard_results_on_any_change = True

27 fileio_rules = FileIOCalculator.ruleset(

28 configspec=dict(skt_path=None),

29 stdout_name='{prefix}.out')

31 def __init__(self, restart=None,

32 ignore_bad_restart_file=FileIOCalculator._deprecated,

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

34 slako_dir=None,

35 command=None,

36 profile=None,

37 **kwargs):

38 """

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

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

42 Hamiltonian = DFTB {

43 SCC = Yes

44 SCCTolerance = 1e-8

45 MaxAngularMomentum = {

46 H = s

47 O = p

48 }

49 }

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

52 following settings:

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

55 >>>

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

57 ... Hamiltonian_SCC='Yes',

58 ... Hamiltonian_SCCTolerance=1e-8,

59 ... Hamiltonian_MaxAngularMomentum_='',

60 ... Hamiltonian_MaxAngularMomentum_H='s',

61 ... Hamiltonian_MaxAngularMomentum_O='p')

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

64 arguments can be used:

66 restart: str

67 Prefix for restart file. May contain a directory.

68 Default is None: don't restart.

69 ignore_bad_restart_file: bool

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

71 error if the restart file is missing or broken.

72 label: str (default 'dftb')

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

74 atoms: Atoms object (default None)

75 Optional Atoms object to which the calculator will be

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

77 unit-cell updated from file.

78 kpts: (default None)

79 Brillouin zone sampling:

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

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

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

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

85 into a Monkhorst-Pack grid using

86 ``ase.calculators.calculator.kpts2sizeandoffsets``

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

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

89 (each with equal weight)

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

93 pcpot: PointCharge object

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

95 """

97 if command is None:

98 if 'DFTB_COMMAND' in self.cfg:

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

100

101 if command is None and profile is None:

102 try:

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

104 except BadConfiguration:

105 pass

106

107 if command is None:

108 command = 'dftb+ > PREFIX.out'

109

110 if slako_dir is None:

111 if profile is not None:

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

113

114 if slako_dir is None:

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

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

117 slako_dir += '/'

118

119 self.slako_dir = slako_dir

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

121 self.default_parameters = dict(

122 Hamiltonian_='DFTB',

123 Hamiltonian_SlaterKosterFiles_='Type2FileNames',

124 Hamiltonian_SlaterKosterFiles_Prefix=self.slako_dir,

125 Hamiltonian_SlaterKosterFiles_Separator='"-"',

126 Hamiltonian_SlaterKosterFiles_Suffix='".skf"',

127 Hamiltonian_MaxAngularMomentum_='',

128 Options_='',

129 Options_WriteResultsTag='Yes',

130 ParserOptions_='',

131 ParserOptions_ParserVersion=1,

132 ParserOptions_IgnoreUnprocessedNodes='Yes')

133 else:

134 self.default_parameters = dict(

135 Options_='',

136 Options_WriteResultsTag='Yes',

137 ParserOptions_='',

138 ParserOptions_ParserVersion=1,

139 ParserOptions_IgnoreUnprocessedNodes='Yes')

140

141 self.pcpot = None

142 self.lines = None

143 self.atoms = None

144 self.atoms_input = None

145 self.do_forces = False

146

147 super().__init__(restart, ignore_bad_restart_file,

148 label, atoms, command=command,

149 profile=profile, **kwargs)

150

151 # Determine number of spin channels

152 try:

153 entry = kwargs['Hamiltonian_SpinPolarisation']

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

155 except KeyError:

156 spinpol = False

157 self.nspin = 2 if spinpol else 1

158

159 # kpoint stuff by ase

160 self.kpts = kpts

161 self.kpts_coord = None

162

163 if self.kpts is not None:

164 initkey = 'Hamiltonian_KPointsAndWeights'

165 mp_mesh = None

166 offsets = None

167

168 if isinstance(self.kpts, dict):

169 if 'path' in self.kpts:

170 # kpts is path in Brillouin zone

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

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

173 else:

174 # kpts is (implicit) definition of

175 # Monkhorst-Pack grid

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

177 mp_mesh, offsets = kpts2sizeandoffsets(atoms=atoms,

178 **self.kpts)

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

180 # kpts is Monkhorst-Pack grid

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

182 mp_mesh = self.kpts

183 offsets = [0.] * 3

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

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

186 # each will be given equal weight

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

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

189 else:

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

191

192 if mp_mesh is not None:

193 eps = 1e-10

194 for i in range(3):

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

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

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

198 offsets[i] *= mp_mesh[i]

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

200 # DFTB+ uses a different offset convention, where

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

202 # to the addition of any offsets

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

204 offsets[i] += 0.5

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

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

207

208 elif self.kpts_coord is not None:

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

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

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

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

213 c_str = '1 ' + c_str

214 else:

215 c_str += ' 1.0'

216 self.parameters[key] = c_str

217

218 def write_dftb_in(self, outfile):

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

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

221 """

222

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

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

225 outfile.write('} \n')

226 outfile.write(' \n')

227

228 params = self.parameters.copy()

229

230 s = 'Hamiltonian_MaxAngularMomentum_'

231 for key in params:

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

233 break

234 else:

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

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

237 # the .skf files:

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

239 for symbol in symbols:

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

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

242 l = read_max_angular_momentum(path)

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

244

245 if self.do_forces:

246 params['Analysis_'] = ''

247 params['Analysis_CalculateForces'] = 'Yes'

248

249 # --------MAIN KEYWORDS-------

250 previous_key = 'dummy_'

251 myspace = ' '

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

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

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

255 for my_backslash in reversed(

256 range(previous_depth - current_depth)):

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

258 outfile.write(3 * current_depth * myspace)

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

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. 'Options {'

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

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

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

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

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

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

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

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

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

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

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

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

276 if not (f1 or f2):

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

278 value = 'No'

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

280 else:

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

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

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

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

285 outfile.write(

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

287 outfile.write(' Records = ' +

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

289 outfile.write(

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

291 outfile.write(' } \n')

292 outfile.write(' } \n')

293 outfile.write(' } \n')

294 previous_key = key

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

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

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

298

299 def check_state(self, atoms):

300 system_changes = FileIOCalculator.check_state(self, atoms)

301 # Ignore unit cell for molecules:

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

303 system_changes.remove('cell')

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

305 system_changes.append('positions')

306 return system_changes

307

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

309 from ase.io import write

310 if properties is not None:

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

312 self.do_forces = True

313 FileIOCalculator.write_input(

314 self, atoms, properties, system_changes)

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

316 self.write_dftb_in(fd)

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

318 parallel=False)

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

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

321 self.atoms_input = atoms

322 self.atoms = None

323 if self.pcpot:

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

325

326 def read_results(self):

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

328 It will be destroyed after it is read to avoid

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

330

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

332 self.lines = fd.readlines()

333

334 self.atoms = self.atoms_input

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

336 if charges is not None:

337 self.results['charges'] = charges

338 self.results['energy'] = energy

339 if dipole is not None:

340 self.results['dipole'] = dipole

341 if self.do_forces:

342 forces = self.read_forces()

343 self.results['forces'] = forces

344 self.mmpositions = None

345

346 # stress stuff begins

347 sstring = 'stress'

348 have_stress = False

349 stress = []

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

351 if sstring in line:

352 have_stress = True

353 start = iline + 1

354 end = start + 3

355 for i in range(start, end):

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

357 stress.append(cell)

358 if have_stress:

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

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

361 # stress stuff ends

362

363 # eigenvalues and fermi levels

364 fermi_levels = self.read_fermi_levels()

365 if fermi_levels is not None:

366 self.results['fermi_levels'] = fermi_levels

367

368 eigenvalues = self.read_eigenvalues()

369 if eigenvalues is not None:

370 self.results['eigenvalues'] = eigenvalues

371

372 # calculation was carried out with atoms written in write_input

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

374

375 def read_forces(self):

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

377 from ase.units import Bohr, Hartree

378

379 # Initialise the indices so their scope

380 # reaches outside of the for loop

381 index_force_begin = -1

382 index_force_end = -1

383

384 # Force line indexes

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

386 fstring = 'forces '

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

388 index_force_begin = iline + 1

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

390 index_force_end = iline + 1 + \

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

392 break

393

394 gradients = []

395 for j in range(index_force_begin, index_force_end):

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

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

398

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

400

401 def read_charges_energy_dipole(self):

402 """Get partial charges on atoms

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

404 """

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

406 lines = fd.readlines()

407

408 for line in lines:

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

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

411 break

412

413 qm_charges = []

414 for n, line in enumerate(lines):

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

416 chargestart = n + 1

417 break

418 else:

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

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

421 return None, energy, None

422

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

424 for line in lines1:

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

426

427 dipole = None

428 for line in lines:

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

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

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

432

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

434

435 def get_charges(self, atoms):

436 """ Get the calculated charges

437 this is inhereted to atoms object """

438 if 'charges' in self.results:

439 return self.results['charges']

440 else:

441 return None

442

443 def read_eigenvalues(self):

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

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

446 # Eigenvalue line indexes

447 index_eig_begin = None

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

449 fstring = 'eigenvalues '

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

451 index_eig_begin = iline + 1

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

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

454 break

455 else:

456 return None

457

458 # Take into account that the last row may lack

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

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

461 index_eig_end = index_eig_begin + nrow

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

463 # XXX dirty fix

464 self.lines[index_eig_end - 1] = (

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

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

467

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

469 eig *= Hartree

470 N = nkpt * nband

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

472 for i in range(nspin)]

473

474 return eigenvalues

475

476 def read_fermi_levels(self):

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

478 # Fermi level line indexes

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

480 fstring = 'fermi_level '

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

482 index_fermi = iline + 1

483 break

484 else:

485 return None

486

487 fermi_levels = []

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

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

490

491 for word in words:

492 e = float(word)

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

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

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

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

497 fermi_levels.append(e)

498

499 return np.array(fermi_levels) * Hartree

500

501 def get_ibz_k_points(self):

502 return self.kpts_coord.copy()

503

504 def get_number_of_spins(self):

505 return self.nspin

506

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

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

509

510 def get_fermi_levels(self):

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

512

513 def get_fermi_level(self):

514 return max(self.get_fermi_levels())

515

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

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

518 """

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

520 return self.pcpot

521

522

523class PointChargePotential:

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

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

526 """

527 self.mmcharges = mmcharges

528 self.directory = directory

529 self.mmpositions = None

530 self.mmforces = None

531

532 def set_positions(self, mmpositions):

533 self.mmpositions = mmpositions

534

535 def set_charges(self, mmcharges):

536 self.mmcharges = mmcharges

537

538 def write_mmcharges(self, filename):

539 """ mok all

540 write external charges as monopoles for dftb+.

541

542 """

543 if self.mmcharges is None:

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

545 return

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

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

548 [x, y, z] = pos

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

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

551

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

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

554 if get_forces:

555 return self.read_forces_on_pointcharges()

556 else:

557 return np.zeros_like(self.mmpositions)

558

559 def read_forces_on_pointcharges(self):

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

561 from ase.units import Bohr, Hartree

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

563 lines = fd.readlines()

564

565 external_forces = []

566 for n, line in enumerate(lines):

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

568 chargestart = n + 1

569 break

570 else:

571 raise RuntimeError(

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

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

574 for line in lines1:

575 external_forces.append(

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

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

578

579

580def read_max_angular_momentum(path):

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

582

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

584 """

585 with open(path) as fd:

586 line = fd.readline()

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

588 # Extended format

589 fd.readline()

590 l = 3

591 pos = 9

592 else:

593 # Simple format:

594 l = 2

595 pos = 7

596

597 # Sometimes there ar commas, sometimes not:

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

599

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

601 for f in occs:

602 if f > 0.0:

603 return l

604 l -= 1