Coverage for /builds/kinetik161/ase/ase/calculators/dftb.py: 75.30%

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 def __init__(self, restart=None,

24 ignore_bad_restart_file=FileIOCalculator._deprecated,

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

26 slako_dir=None,

27 command=None,

28 profile=None,

29 **kwargs):

30 """

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

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

34 Hamiltonian = DFTB {

35 SCC = Yes

36 SCCTolerance = 1e-8

37 MaxAngularMomentum = {

38 H = s

39 O = p

40 }

41 }

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

44 following settings:

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

47 >>>

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

49 ... Hamiltonian_SCC='Yes',

50 ... Hamiltonian_SCCTolerance=1e-8,

51 ... Hamiltonian_MaxAngularMomentum_='',

52 ... Hamiltonian_MaxAngularMomentum_H='s',

53 ... Hamiltonian_MaxAngularMomentum_O='p')

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

56 arguments can be used:

58 restart: str

59 Prefix for restart file. May contain a directory.

60 Default is None: don't restart.

61 ignore_bad_restart_file: bool

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

63 error if the restart file is missing or broken.

64 label: str (default 'dftb')

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

66 atoms: Atoms object (default None)

67 Optional Atoms object to which the calculator will be

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

69 unit-cell updated from file.

70 kpts: (default None)

71 Brillouin zone sampling:

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

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

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

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

77 into a Monkhorst-Pack grid using

78 ``ase.calculators.calculator.kpts2sizeandoffsets``

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

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

81 (each with equal weight)

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

85 pcpot: PointCharge object

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

87 """

89 if command is None:

90 if 'DFTB_COMMAND' in self.cfg:

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

92 else:

93 command = 'dftb+ > PREFIX.out'

95 if slako_dir is None:

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

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

98 slako_dir += '/'

100 self.slako_dir = slako_dir

101

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

103 self.default_parameters = dict(

104 Hamiltonian_='DFTB',

105 Hamiltonian_SlaterKosterFiles_='Type2FileNames',

106 Hamiltonian_SlaterKosterFiles_Prefix=self.slako_dir,

107 Hamiltonian_SlaterKosterFiles_Separator='"-"',

108 Hamiltonian_SlaterKosterFiles_Suffix='".skf"',

109 Hamiltonian_MaxAngularMomentum_='',

110 Options_='',

111 Options_WriteResultsTag='Yes')

112 else:

113 self.default_parameters = dict(

114 Options_='',

115 Options_WriteResultsTag='Yes')

116

117 self.pcpot = None

118 self.lines = None

119 self.atoms = None

120 self.atoms_input = None

121 self.do_forces = False

122 self.outfilename = 'dftb.out'

123

124 super().__init__(restart, ignore_bad_restart_file,

125 label, atoms, command=command,

126 profile=profile, **kwargs)

127

128 # Determine number of spin channels

129 try:

130 entry = kwargs['Hamiltonian_SpinPolarisation']

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

132 except KeyError:

133 spinpol = False

134 self.nspin = 2 if spinpol else 1

135

136 # kpoint stuff by ase

137 self.kpts = kpts

138 self.kpts_coord = None

139

140 if self.kpts is not None:

141 initkey = 'Hamiltonian_KPointsAndWeights'

142 mp_mesh = None

143 offsets = None

144

145 if isinstance(self.kpts, dict):

146 if 'path' in self.kpts:

147 # kpts is path in Brillouin zone

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

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

150 else:

151 # kpts is (implicit) definition of

152 # Monkhorst-Pack grid

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

154 mp_mesh, offsets = kpts2sizeandoffsets(atoms=atoms,

155 **self.kpts)

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

157 # kpts is Monkhorst-Pack grid

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

159 mp_mesh = self.kpts

160 offsets = [0.] * 3

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

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

163 # each will be given equal weight

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

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

166 else:

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

168

169 if mp_mesh is not None:

170 eps = 1e-10

171 for i in range(3):

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

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

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

175 offsets[i] *= mp_mesh[i]

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

177 # DFTB+ uses a different offset convention, where

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

179 # to the addition of any offsets

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

181 offsets[i] += 0.5

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

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

184

185 elif self.kpts_coord is not None:

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

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

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

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

190 c_str = '1 ' + c_str

191 else:

192 c_str += ' 1.0'

193 self.parameters[key] = c_str

194

195 def write_dftb_in(self, outfile):

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

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

198 """

199

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

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

202 outfile.write('} \n')

203 outfile.write(' \n')

204

205 params = self.parameters.copy()

206

207 s = 'Hamiltonian_MaxAngularMomentum_'

208 for key in params:

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

210 break

211 else:

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

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

214 # the .skf files:

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

216 for symbol in symbols:

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

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

219 l = read_max_angular_momentum(path)

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

221

222 # --------MAIN KEYWORDS-------

223 previous_key = 'dummy_'

224 myspace = ' '

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

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

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

228 for my_backsclash in reversed(

229 range(previous_depth - current_depth)):

230 outfile.write(3 * (1 + my_backsclash) * myspace + '} \n')

231 outfile.write(3 * current_depth * myspace)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

249 if not (f1 or f2):

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

251 value = 'No'

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

253 else:

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

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

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

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

258 outfile.write(

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

260 outfile.write(' Records = ' +

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

262 outfile.write(

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

264 outfile.write(' } \n')

265 outfile.write(' } \n')

266 outfile.write(' } \n')

267 previous_key = key

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

269 for my_backsclash in reversed(range(current_depth)):

270 outfile.write(3 * my_backsclash * myspace + '} \n')

271 outfile.write('ParserOptions { \n')

272 outfile.write(' IgnoreUnprocessedNodes = Yes \n')

273 outfile.write('} \n')

274 if self.do_forces:

275 outfile.write('Analysis { \n')

276 outfile.write(' CalculateForces = Yes \n')

277 outfile.write('} \n')

278

279 def check_state(self, atoms):

280 system_changes = FileIOCalculator.check_state(self, atoms)

281 # Ignore unit cell for molecules:

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

283 system_changes.remove('cell')

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

285 system_changes.append('positions')

286 return system_changes

287

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

289 from ase.io import write

290 if properties is not None:

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

292 self.do_forces = True

293 FileIOCalculator.write_input(

294 self, atoms, properties, system_changes)

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

296 self.write_dftb_in(fd)

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

298 parallel=False)

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

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

301 self.atoms_input = atoms

302 self.atoms = None

303 if self.pcpot:

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

305

306 def read_results(self):

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

308 It will be destroyed after it is read to avoid

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

310

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

312 self.lines = fd.readlines()

313

314 self.atoms = self.atoms_input

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

316 if charges is not None:

317 self.results['charges'] = charges

318 self.results['energy'] = energy

319 if dipole is not None:

320 self.results['dipole'] = dipole

321 if self.do_forces:

322 forces = self.read_forces()

323 self.results['forces'] = forces

324 self.mmpositions = None

325

326 # stress stuff begins

327 sstring = 'stress'

328 have_stress = False

329 stress = []

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

331 if sstring in line:

332 have_stress = True

333 start = iline + 1

334 end = start + 3

335 for i in range(start, end):

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

337 stress.append(cell)

338 if have_stress:

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

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

341 # stress stuff ends

342

343 # eigenvalues and fermi levels

344 fermi_levels = self.read_fermi_levels()

345 if fermi_levels is not None:

346 self.results['fermi_levels'] = fermi_levels

347

348 eigenvalues = self.read_eigenvalues()

349 if eigenvalues is not None:

350 self.results['eigenvalues'] = eigenvalues

351

352 # calculation was carried out with atoms written in write_input

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

354

355 def read_forces(self):

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

357 from ase.units import Bohr, Hartree

358

359 # Initialise the indices so their scope

360 # reaches outside of the for loop

361 index_force_begin = -1

362 index_force_end = -1

363

364 # Force line indexes

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

366 fstring = 'forces '

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

368 index_force_begin = iline + 1

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

370 index_force_end = iline + 1 + \

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

372 break

373

374 gradients = []

375 for j in range(index_force_begin, index_force_end):

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

377 gradients.append([float(word[k]) for k in range(0, 3)])

378

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

380

381 def read_charges_energy_dipole(self):

382 """Get partial charges on atoms

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

384 """

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

386 lines = fd.readlines()

387

388 for line in lines:

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

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

391 break

392

393 qm_charges = []

394 for n, line in enumerate(lines):

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

396 chargestart = n + 1

397 break

398 else:

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

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

401 return None, energy, None

402

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

404 for line in lines1:

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

406

407 dipole = None

408 for line in lines:

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

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

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

412

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

414

415 def get_charges(self, atoms):

416 """ Get the calculated charges

417 this is inhereted to atoms object """

418 if 'charges' in self.results:

419 return self.results['charges']

420 else:

421 return None

422

423 def read_eigenvalues(self):

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

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

426 # Eigenvalue line indexes

427 index_eig_begin = None

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

429 fstring = 'eigenvalues '

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

431 index_eig_begin = iline + 1

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

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

434 break

435 else:

436 return None

437

438 # Take into account that the last row may lack

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

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

441 index_eig_end = index_eig_begin + nrow

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

443 # XXX dirty fix

444 self.lines[index_eig_end - 1] = (

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

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

447

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

449 eig *= Hartree

450 N = nkpt * nband

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

452 for i in range(nspin)]

453

454 return eigenvalues

455

456 def read_fermi_levels(self):

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

458 # Fermi level line indexes

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

460 fstring = 'fermi_level '

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

462 index_fermi = iline + 1

463 break

464 else:

465 return None

466

467 fermi_levels = []

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

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

470

471 for word in words:

472 e = float(word)

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

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

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

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

477 fermi_levels.append(e)

478

479 return np.array(fermi_levels) * Hartree

480

481 def get_ibz_k_points(self):

482 return self.kpts_coord.copy()

483

484 def get_number_of_spins(self):

485 return self.nspin

486

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

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

489

490 def get_fermi_levels(self):

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

492

493 def get_fermi_level(self):

494 return max(self.get_fermi_levels())

495

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

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

498 """

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

500 return self.pcpot

501

502

503class PointChargePotential:

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

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

506 """

507 self.mmcharges = mmcharges

508 self.directory = directory

509 self.mmpositions = None

510 self.mmforces = None

511

512 def set_positions(self, mmpositions):

513 self.mmpositions = mmpositions

514

515 def set_charges(self, mmcharges):

516 self.mmcharges = mmcharges

517

518 def write_mmcharges(self, filename):

519 """ mok all

520 write external charges as monopoles for dftb+.

521

522 """

523 if self.mmcharges is None:

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

525 return

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

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

528 [x, y, z] = pos

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

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

531

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

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

534 if get_forces:

535 return self.read_forces_on_pointcharges()

536 else:

537 return np.zeros_like(self.mmpositions)

538

539 def read_forces_on_pointcharges(self):

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

541 from ase.units import Bohr, Hartree

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

543 lines = fd.readlines()

544

545 external_forces = []

546 for n, line in enumerate(lines):

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

548 chargestart = n + 1

549 break

550 else:

551 raise RuntimeError(

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

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

554 for line in lines1:

555 external_forces.append(

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

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

558

559

560def read_max_angular_momentum(path):

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

562

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

564 """

565 with open(path) as fd:

566 line = fd.readline()

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

568 # Extended format

569 fd.readline()

570 l = 3

571 pos = 9

572 else:

573 # Simple format:

574 l = 2

575 pos = 7

576

577 # Sometimes there ar commas, sometimes not:

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

579

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

581 for f in occs:

582 if f > 0.0:

583 return l

584 l -= 1