Coverage for /builds/ase/ase/ase/calculators/dftd3.py: 87.13%

1# fmt: off

3import os

4import subprocess

5from pathlib import Path

6from warnings import warn

8import numpy as np

10from ase.calculators.calculator import (

11 BaseCalculator,

12 Calculator,

13 FileIOCalculator,

14)

15from ase.io import write

16from ase.io.vasp import write_vasp

17from ase.parallel import world

18from ase.units import Bohr, Hartree

21def dftd3_defaults():

22 default_parameters = {'xc': None, # PBE if no custom damping parameters

23 'grad': True, # calculate forces/stress

24 'abc': False, # ATM 3-body contribution

25 'cutoff': 95 * Bohr, # Cutoff for 2-body calcs

26 'cnthr': 40 * Bohr, # Cutoff for 3-body and CN calcs

27 'old': False, # use old DFT-D2 method instead

28 'damping': 'zero', # Default to zero-damping

29 'tz': False, # 'triple zeta' alt. parameters

30 's6': None, # damping parameters start here

31 'sr6': None,

32 's8': None,

33 'sr8': None,

34 'alpha6': None,

35 'a1': None,

36 'a2': None,

37 'beta': None}

38 return default_parameters

41class DFTD3(BaseCalculator):

42 """Grimme DFT-D3 calculator"""

44 def __init__(self,

45 label='ase_dftd3', # Label for dftd3 output files

46 command=None, # Command for running dftd3

47 dft=None, # DFT calculator

48 comm=world,

49 **kwargs):

51 # Convert from 'func' keyword to 'xc'. Internally, we only store

52 # 'xc', but 'func' is also allowed since it is consistent with the

53 # CLI dftd3 interface.

54 func = kwargs.pop('func', None)

55 if func is not None:

56 if kwargs.get('xc') is not None:

57 raise RuntimeError('Both "func" and "xc" were provided! '

58 'Please provide at most one of these '

59 'two keywords. The preferred keyword '

60 'is "xc"; "func" is allowed for '

61 'consistency with the CLI dftd3 '

62 'interface.')

63 kwargs['xc'] = func

65 # If the user did not supply an XC functional, but did attach a

66 # DFT calculator that has XC set, then we will use that. Note that

67 # DFTD3's spelling convention is different from most, so in general

68 # you will have to explicitly set XC for both the DFT calculator and

69 # for DFTD3 (and DFTD3's will likely be spelled differently...)

70 if dft is not None and kwargs.get('xc') is None:

71 dft_xc = dft.parameters.get('xc')

72 if dft_xc is not None:

73 kwargs['xc'] = dft_xc

75 dftd3 = PureDFTD3(label=label, command=command, comm=comm, **kwargs)

77 # dftd3 only implements energy, forces, and stresses (for periodic

78 # systems). But, if a DFT calculator is attached, and that calculator

79 # implements more properties, we expose those properties.

80 # dftd3 contributions for those properties will be zero.

81 if dft is None:

82 self.implemented_properties = list(dftd3.dftd3_properties)

83 else:

84 self.implemented_properties = list(dft.implemented_properties)

86 # Should our arguments be "parameters" (passed to superclass)

87 # or are they not really "parameters"?

88 #

89 # That's not really well defined. Let's not do anything then.

90 super().__init__()

92 self.dftd3 = dftd3

93 self.dft = dft

95 def todict(self):

96 return {}

98 def calculate(self, atoms, properties, system_changes):

99 common_props = set(self.dftd3.dftd3_properties) & set(properties)

100 dftd3_results = self._get_properties(atoms, common_props, self.dftd3)

101

102 if self.dft is None:

103 results = dftd3_results

104 else:

105 dft_results = self._get_properties(atoms, properties, self.dft)

106 results = dict(dft_results)

107 for name in set(results) & set(dftd3_results):

108 assert np.shape(results[name]) == np.shape(dftd3_results[name])

109 results[name] += dftd3_results[name]

110

111 # Although DFTD3 may have calculated quantities not provided

112 # by the calculator (e.g. stress), it would be wrong to

113 # return those! Return only what corresponds to the DFT calc.

114 assert set(results) == set(dft_results)

115 self.results = results

116

117 def _get_properties(self, atoms, properties, calc):

118 # We want any and all properties that the calculator

119 # normally produces. So we intend to rob the calc.results

120 # dictionary instead of only getting the requested properties.

121

122 import copy

123 for name in properties:

124 calc.get_property(name, atoms)

125 assert name in calc.results

126

127 # XXX maybe use get_properties() when that makes sense.

128 results = copy.deepcopy(calc.results)

129 assert set(properties) <= set(results)

130 return results

131

132

133class PureDFTD3(FileIOCalculator):

134 """DFTD3 calculator without corresponding DFT contribution.

135

136 This class is an implementation detail."""

137

138 name = 'puredftd3'

139

140 dftd3_properties = {'energy', 'free_energy', 'forces', 'stress'}

141 implemented_properties = list(dftd3_properties)

142 default_parameters = dftd3_defaults()

143 damping_methods = {'zero', 'bj', 'zerom', 'bjm'}

144 _legacy_default_command = 'dftd3'

145

146 def __init__(self,

147 *,

148 label='ase_dftd3', # Label for dftd3 output files

149 command=None, # Command for running dftd3

150 comm=world,

151 **kwargs):

152

153 # FileIOCalculator would default to self.name to get the envvar

154 # which determines the command.

155 # We'll have to overrule that if we want to keep scripts working:

156 command = command or self.cfg.get('ASE_DFTD3_COMMAND')

157

158 super().__init__(label=label,

159 command=command,

160 **kwargs)

161

162 # TARP: This is done because the calculator does not call

163 # FileIOCalculator.calculate, but Calculator.calculate and does not

164 # use the profile defined in FileIOCalculator.__init__

165 self.comm = comm

166

167 def set(self, **kwargs):

168 changed_parameters = {}

169

170 # Check for unknown arguments. Don't raise an error, just let the

171 # user know that we don't understand what they're asking for.

172 unknown_kwargs = set(kwargs) - set(self.default_parameters)

173 if unknown_kwargs:

174 warn('WARNING: Ignoring the following unknown keywords: {}'

175 ''.format(', '.join(unknown_kwargs)))

176

177 changed_parameters.update(FileIOCalculator.set(self, **kwargs))

178

179 # Ensure damping method is valid (zero, bj, zerom, bjm).

180 damping = self.parameters['damping']

181 if damping is not None:

182 damping = damping.lower()

183 if damping not in self.damping_methods:

184 raise ValueError(f'Unknown damping method {damping}!')

185

186 # d2 only is valid with 'zero' damping

187 elif self.parameters['old'] and damping != 'zero':

188 raise ValueError('Only zero-damping can be used with the D2 '

189 'dispersion correction method!')

190

191 # If cnthr (cutoff for three-body and CN calculations) is greater

192 # than cutoff (cutoff for two-body calculations), then set the former

193 # equal to the latter, since that doesn't make any sense.

194 if self.parameters['cnthr'] > self.parameters['cutoff']:

195 warn('WARNING: CN cutoff value of {cnthr} is larger than '

196 'regular cutoff value of {cutoff}! Reducing CN cutoff '

197 'to {cutoff}.'

198 ''.format(cnthr=self.parameters['cnthr'],

199 cutoff=self.parameters['cutoff']))

200 self.parameters['cnthr'] = self.parameters['cutoff']

201

202 # If you only care about the energy, gradient calculations (forces,

203 # stresses) can be bypassed. This will greatly speed up calculations

204 # in dense 3D-periodic systems with three-body corrections. But, we

205 # can no longer say that we implement forces and stresses.

206 # if not self.parameters['grad']:

207 # for val in ['forces', 'stress']:

208 # if val in self.implemented_properties:

209 # self.implemented_properties.remove(val)

210

211 # Check to see if we're using custom damping parameters.

212 zero_damppars = {'s6', 'sr6', 's8', 'sr8', 'alpha6'}

213 bj_damppars = {'s6', 'a1', 's8', 'a2', 'alpha6'}

214 zerom_damppars = {'s6', 'sr6', 's8', 'beta', 'alpha6'}

215 all_damppars = zero_damppars | bj_damppars | zerom_damppars

216

217 self.custom_damp = False

218

219 damppars = set(kwargs) & all_damppars

220 if damppars:

221 self.custom_damp = True

222 if damping == 'zero':

223 valid_damppars = zero_damppars

224 elif damping in ['bj', 'bjm']:

225 valid_damppars = bj_damppars

226 elif damping == 'zerom':

227 valid_damppars = zerom_damppars

228

229 # If some but not all damping parameters are provided for the

230 # selected damping method, raise an error. We don't have "default"

231 # values for damping parameters, since those are stored in the

232 # dftd3 executable & depend on XC functional.

233 missing_damppars = valid_damppars - damppars

234 if missing_damppars and missing_damppars != valid_damppars:

235 raise ValueError('An incomplete set of custom damping '

236 'parameters for the {} damping method was '

237 'provided! Expected: {}; got: {}'

238 ''.format(damping,

239 ', '.join(valid_damppars),

240 ', '.join(damppars)))

241

242 # If a user provides damping parameters that are not used in the

243 # selected damping method, let them know that we're ignoring them.

244 # If the user accidentally provided the *wrong* set of parameters,

245 # (e.g., the BJ parameters when they are using zero damping), then

246 # the previous check will raise an error, so we don't need to

247 # worry about that here.

248 if damppars - valid_damppars:

249 warn('WARNING: The following damping parameters are not '

250 'valid for the {} damping method and will be ignored: {}'

251 ''.format(damping,

252 ', '.join(damppars)))

253

254 # The default XC functional is PBE, but this is only set if the user

255 # did not provide their own value for xc or any custom damping

256 # parameters.

257 if self.parameters['xc'] and self.custom_damp:

258 warn('WARNING: Custom damping parameters will be used '

259 'instead of those parameterized for {}!'

260 ''.format(self.parameters['xc']))

261

262 if changed_parameters:

263 self.results.clear()

264 return changed_parameters

265

266 def calculate(self, atoms, properties, system_changes):

267 # We don't call FileIOCalculator.calculate here, because that method

268 # calls subprocess.call(..., shell=True), which we don't want to do.

269 # So, we reproduce some content from that method here.

270 Calculator.calculate(self, atoms, properties, system_changes)

271

272 # If a parameter file exists in the working directory, delete it

273 # first. If we need that file, we'll recreate it later.

274 localparfile = os.path.join(self.directory, '.dftd3par.local')

275 if self.comm.rank == 0 and os.path.isfile(localparfile):

276 os.remove(localparfile)

277

278 # Write XYZ or POSCAR file and .dftd3par.local file if we are using

279 # custom damping parameters.

280 self.write_input(self.atoms, properties, system_changes)

281 # command = self._generate_command()

282

283 inputs = DFTD3Inputs(command=self.command, prefix=self.label,

284 atoms=self.atoms, parameters=self.parameters)

285 command = inputs.get_argv(custom_damp=self.custom_damp)

286

287 # Finally, call dftd3 and parse results.

288 # DFTD3 does not run in parallel

289 # so we only need it to run on 1 core

290 errorcode = 0

291 if self.comm.rank == 0:

292 with open(self.label + '.out', 'w') as fd:

293 errorcode = subprocess.call(command,

294 cwd=self.directory, stdout=fd)

295

296 errorcode = self.comm.sum_scalar(errorcode)

297

298 if errorcode:

299 raise RuntimeError('%s returned an error: %d' %

300 (self.name, errorcode))

301

302 self.read_results()

303

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

305 FileIOCalculator.write_input(self, atoms, properties=properties,

306 system_changes=system_changes)

307 # dftd3 can either do fully 3D periodic or non-periodic calculations.

308 # It cannot do calculations that are only periodic in 1 or 2

309 # dimensions. If the atoms object is periodic in only 1 or 2

310 # dimensions, then treat it as a fully 3D periodic system, but warn

311 # the user.

312

313 if self.custom_damp:

314 damppars = _get_damppars(self.parameters)

315 else:

316 damppars = None

317

318 if self.comm.rank == 0:

319 self._actually_write_input(

320 directory=Path(self.directory), atoms=atoms,

321 properties=properties, prefix=self.label,

322 damppars=damppars, pbc=any(atoms.pbc))

323

324 def _actually_write_input(self, directory, prefix, atoms, properties,

325 damppars, pbc: bool):

326 if pbc:

327 fname = directory / f'{prefix}.POSCAR'

328 # We sort the atoms so that the atomtypes list becomes as

329 # short as possible. The dftd3 program can only handle 10

330 # atomtypes

331 write_vasp(fname, atoms, sort=True)

332 else:

333 fname = directory / f'{prefix}.xyz'

334 write(fname, atoms, format='xyz', parallel=False)

335

336 # Generate custom damping parameters file. This is kind of ugly, but

337 # I don't know of a better way of doing this.

338 if damppars is not None:

339 damp_fname = directory / '.dftd3par.local'

340 with open(damp_fname, 'w') as fd:

341 fd.write(' '.join(damppars))

342

343 def _outname(self):

344 return Path(self.directory) / f'{self.label}.out'

345

346 def _read_and_broadcast_results(self):

347 from ase.parallel import broadcast

348 if self.comm.rank == 0:

349 output = DFTD3Output(directory=self.directory,

350 stdout_path=self._outname())

351 dct = output.read(atoms=self.atoms,

352 read_forces=bool(self.parameters['grad']))

353 else:

354 dct = None

355

356 dct = broadcast(dct, root=0, comm=self.comm)

357 return dct

358

359 def read_results(self):

360 results = self._read_and_broadcast_results()

361 self.results = results

362

363

364class DFTD3Inputs:

365 dftd3_flags = {'grad', 'pbc', 'abc', 'old', 'tz'}

366

367 def __init__(self, command, prefix, atoms, parameters):

368 self.command = command

369 self.prefix = prefix

370 self.atoms = atoms

371 self.parameters = parameters

372

373 @property

374 def pbc(self):

375 return any(self.atoms.pbc)

376

377 @property

378 def inputformat(self):

379 if self.pbc:

380 return 'POSCAR'

381 else:

382 return 'xyz'

383

384 def get_argv(self, custom_damp):

385 argv = self.command.split()

386

387 argv.append(f'{self.prefix}.{self.inputformat}')

388

389 if not custom_damp:

390 xc = self.parameters.get('xc')

391 if xc is None:

392 xc = 'pbe'

393 argv += ['-func', xc.lower()]

394

395 for arg in self.dftd3_flags:

396 if self.parameters.get(arg):

397 argv.append('-' + arg)

398

399 if self.pbc:

400 argv.append('-pbc')

401

402 argv += ['-cnthr', str(self.parameters['cnthr'] / Bohr)]

403 argv += ['-cutoff', str(self.parameters['cutoff'] / Bohr)]

404

405 if not self.parameters['old']:

406 argv.append('-' + self.parameters['damping'])

407

408 return argv

409

410

411class DFTD3Output:

412 def __init__(self, directory, stdout_path):

413 self.directory = Path(directory)

414 self.stdout_path = Path(stdout_path)

415

416 def read(self, *, atoms, read_forces):

417 results = {}

418

419 energy = self.read_energy()

420 results['energy'] = energy

421 results['free_energy'] = energy

422

423 if read_forces:

424 results['forces'] = self.read_forces(atoms)

425

426 if any(atoms.pbc):

427 results['stress'] = self.read_stress(atoms.cell)

428

429 return results

430

431 def read_forces(self, atoms):

432 forcename = self.directory / 'dftd3_gradient'

433 with open(forcename) as fd:

434 forces = self.parse_forces(fd)

435

436 assert len(forces) == len(atoms)

437

438 forces *= -Hartree / Bohr

439 # XXXX ordering!

440 if any(atoms.pbc):

441 # This seems to be due to vasp file sorting.

442 # If that sorting rule changes, we will get garbled

443 # forces!

444 ind = np.argsort(atoms.symbols)

445 forces[ind] = forces.copy()

446 return forces

447

448 def read_stress(self, cell):

449 volume = cell.volume

450 assert volume > 0

451

452 stress = self.read_cellgradient()

453 stress *= Hartree / Bohr / volume

454 stress = stress.T @ cell

455 return stress.flat[[0, 4, 8, 5, 2, 1]]

456

457 def read_cellgradient(self):

458 with (self.directory / 'dftd3_cellgradient').open() as fd:

459 return self.parse_cellgradient(fd)

460

461 def read_energy(self) -> float:

462 with self.stdout_path.open() as fd:

463 return self.parse_energy(fd, self.stdout_path)

464

465 def parse_energy(self, fd, outname):

466 for line in fd:

467 if line.startswith(' program stopped'):

468 if 'functional name unknown' in line:

469 message = ('Unknown DFTD3 functional name. '

470 'Please check the dftd3.f source file '

471 'for the list of known functionals '

472 'and their spelling.')

473 else:

474 message = ('dftd3 failed! Please check the {} '

475 'output file and report any errors '

476 'to the ASE developers.'

477 ''.format(outname))

478 raise RuntimeError(message)

479

480 if line.startswith(' Edisp'):

481 # line looks something like this:

482 #

483 # Edisp /kcal,au,ev: xxx xxx xxx

484 #

485 parts = line.split()

486 assert parts[1][0] == '/'

487 index = 2 + parts[1][1:-1].split(',').index('au')

488 e_dftd3 = float(parts[index]) * Hartree

489 return e_dftd3

490

491 raise RuntimeError('Could not parse energy from dftd3 '

492 'output, see file {}'.format(outname))

493

494 def parse_forces(self, fd):

495 forces = []

496 for i, line in enumerate(fd):

497 forces.append(line.split())

498 return np.array(forces, dtype=float)

499

500 def parse_cellgradient(self, fd):

501 stress = np.zeros((3, 3))

502 for i, line in enumerate(fd):

503 for j, x in enumerate(line.split()):

504 stress[i, j] = float(x)

505 # Check if all stress elements are present?

506 # Check if file is longer?

507 return stress

508

509

510def _get_damppars(par):

511 damping = par['damping']

512

513 damppars = []

514

515 # s6 is always first

516 damppars.append(str(float(par['s6'])))

517

518 # sr6 is the second value for zero{,m} damping, a1 for bj{,m}

519 if damping in ['zero', 'zerom']:

520 damppars.append(str(float(par['sr6'])))

521 elif damping in ['bj', 'bjm']:

522 damppars.append(str(float(par['a1'])))

523

524 # s8 is always third

525 damppars.append(str(float(par['s8'])))

526

527 # sr8 is fourth for zero, a2 for bj{,m}, beta for zerom

528 if damping == 'zero':

529 damppars.append(str(float(par['sr8'])))

530 elif damping in ['bj', 'bjm']:

531 damppars.append(str(float(par['a2'])))

532 elif damping == 'zerom':

533 damppars.append(str(float(par['beta'])))

534 # alpha6 is always fifth

535 damppars.append(str(int(par['alpha6'])))

536

537 # last is the version number

538 if par['old']:

539 damppars.append('2')

540 elif damping == 'zero':

541 damppars.append('3')

542 elif damping == 'bj':

543 damppars.append('4')

544 elif damping == 'zerom':

545 damppars.append('5')

546 elif damping == 'bjm':

547 damppars.append('6')

548 return damppars