Coverage for /builds/kinetik161/ase/ase/calculators/dftd3.py: 86.91%

1import os

2import subprocess

3from pathlib import Path

4from warnings import warn

6import numpy as np

8from ase.calculators.calculator import (BaseCalculator, Calculator,

9 FileIOCalculator)

10from ase.io import write

11from ase.io.vasp import write_vasp

12from ase.parallel import world

13from ase.units import Bohr, Hartree

16def dftd3_defaults():

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

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

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

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

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

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

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

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

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

26 'sr6': None,

27 's8': None,

28 'sr8': None,

29 'alpha6': None,

30 'a1': None,

31 'a2': None,

32 'beta': None}

33 return default_parameters

36class DFTD3(BaseCalculator):

37 """Grimme DFT-D3 calculator"""

39 def __init__(self,

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

41 command=None, # Command for running dftd3

42 dft=None, # DFT calculator

43 comm=world,

44 **kwargs):

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

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

48 # CLI dftd3 interface.

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

50 if func is not None:

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

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

53 'Please provide at most one of these '

54 'two keywords. The preferred keyword '

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

56 'consistency with the CLI dftd3 '

57 'interface.')

58 kwargs['xc'] = func

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

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

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

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

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

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

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

67 if dft_xc is not None:

68 kwargs['xc'] = dft_xc

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

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

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

74 # implements more properties, we expose those properties.

75 # dftd3 contributions for those properties will be zero.

76 if dft is None:

77 self.implemented_properties = list(dftd3.dftd3_properties)

78 else:

79 self.implemented_properties = list(dft.implemented_properties)

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

82 # or are they not really "parameters"?

83 #

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

85 super().__init__()

87 self.dftd3 = dftd3

88 self.dft = dft

90 def todict(self):

91 return {}

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

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

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

97 if self.dft is None:

98 results = dftd3_results

99 else:

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

101 results = dict(dft_results)

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

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

104 results[name] += dftd3_results[name]

105

106 # Although DFTD3 may have calculated quantities not provided

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

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

109 assert set(results) == set(dft_results)

110 self.results = results

111

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

113 # We want any and all properties that the calculator

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

115 # dictionary instead of only getting the requested properties.

116

117 import copy

118 for name in properties:

119 calc.get_property(name, atoms)

120 assert name in calc.results

121

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

123 results = copy.deepcopy(calc.results)

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

125 return results

126

127

128class PureDFTD3(FileIOCalculator):

129 """DFTD3 calculator without corresponding DFT contribution.

130

131 This class is an implementation detail."""

132

133 name = 'puredftd3'

134

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

136 implemented_properties = list(dftd3_properties)

137 default_parameters = dftd3_defaults()

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

139 _legacy_default_command = 'dftd3'

140

141 def __init__(self,

142 *,

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

144 command=None, # Command for running dftd3

145 comm=world,

146 **kwargs):

147

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

149 # which determines the command.

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

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

152

153 super().__init__(label=label,

154 command=command,

155 **kwargs)

156

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

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

159 # use the profile defined in FileIOCalculator.__init__

160 self.comm = comm

161

162 def set(self, **kwargs):

163 changed_parameters = {}

164

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

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

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

168 if unknown_kwargs:

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

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

171

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

173

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

175 damping = self.parameters['damping']

176 if damping is not None:

177 damping = damping.lower()

178 if damping not in self.damping_methods:

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

180

181 # d2 only is valid with 'zero' damping

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

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

184 'dispersion correction method!')

185

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

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

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

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

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

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

192 'to {cutoff}.'

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

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

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

196

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

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

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

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

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

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

203 # if val in self.implemented_properties:

204 # self.implemented_properties.remove(val)

205

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

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

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

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

210 all_damppars = zero_damppars | bj_damppars | zerom_damppars

211

212 self.custom_damp = False

213

214 damppars = set(kwargs) & all_damppars

215 if damppars:

216 self.custom_damp = True

217 if damping == 'zero':

218 valid_damppars = zero_damppars

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

220 valid_damppars = bj_damppars

221 elif damping == 'zerom':

222 valid_damppars = zerom_damppars

223

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

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

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

227 # dftd3 executable & depend on XC functional.

228 missing_damppars = valid_damppars - damppars

229 if missing_damppars and missing_damppars != valid_damppars:

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

231 'parameters for the {} damping method was '

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

233 ''.format(damping,

234 ', '.join(valid_damppars),

235 ', '.join(damppars)))

236

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

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

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

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

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

242 # worry about that here.

243 if damppars - valid_damppars:

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

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

246 ''.format(damping,

247 ', '.join(damppars)))

248

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

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

251 # parameters.

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

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

254 'instead of those parameterized for {}!'

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

256

257 if changed_parameters:

258 self.results.clear()

259 return changed_parameters

260

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

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

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

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

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

266

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

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

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

270 if world.rank == 0 and os.path.isfile(localparfile):

271 os.remove(localparfile)

272

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

274 # custom damping parameters.

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

276 # command = self._generate_command()

277

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

279 atoms=self.atoms, parameters=self.parameters)

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

281

282 # Finally, call dftd3 and parse results.

283 # DFTD3 does not run in parallel

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

285 errorcode = 0

286 if self.comm.rank == 0:

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

288 errorcode = subprocess.call(command,

289 cwd=self.directory, stdout=fd)

290

291 errorcode = self.comm.sum_scalar(errorcode)

292

293 if errorcode:

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

295 (self.name, errorcode))

296

297 self.read_results()

298

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

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

301 system_changes=system_changes)

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

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

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

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

306 # the user.

307

308 if self.custom_damp:

309 damppars = _get_damppars(self.parameters)

310 else:

311 damppars = None

312

313 pbc = any(atoms.pbc)

314 if pbc and not all(atoms.pbc):

315 warn('WARNING! dftd3 can only calculate the dispersion energy '

316 'of non-periodic or 3D-periodic systems. We will treat '

317 'this system as 3D-periodic!')

318

319 if self.comm.rank == 0:

320 self._actually_write_input(

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

322 properties=properties, prefix=self.label,

323 damppars=damppars, pbc=pbc)

324

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

326 damppars, pbc):

327 if pbc:

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

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

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

331 # atomtypes

332 write_vasp(fname, atoms, sort=True)

333 else:

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

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

336

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

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

339 if damppars is not None:

340 damp_fname = directory / '.dftd3par.local'

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

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

343

344 def _outname(self):

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

346

347 def _read_and_broadcast_results(self):

348 from ase.parallel import broadcast

349 if self.comm.rank == 0:

350 output = DFTD3Output(directory=self.directory,

351 stdout_path=self._outname())

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

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

354 else:

355 dct = None

356

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

358 return dct

359

360 def read_results(self):

361 results = self._read_and_broadcast_results()

362 self.results = results

363

364

365class DFTD3Inputs:

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

367

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

369 self.command = command

370 self.prefix = prefix

371 self.atoms = atoms

372 self.parameters = parameters

373

374 @property

375 def pbc(self):

376 return any(self.atoms.pbc)

377

378 @property

379 def inputformat(self):

380 if self.pbc:

381 return 'POSCAR'

382 else:

383 return 'xyz'

384

385 def get_argv(self, custom_damp):

386 argv = self.command.split()

387

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

389

390 if not custom_damp:

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

392 if xc is None:

393 xc = 'pbe'

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

395

396 for arg in self.dftd3_flags:

397 if self.parameters.get(arg):

398 argv.append('-' + arg)

399

400 if self.pbc:

401 argv.append('-pbc')

402

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

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

405

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

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

408

409 return argv

410

411

412class DFTD3Output:

413 def __init__(self, directory, stdout_path):

414 self.directory = Path(directory)

415 self.stdout_path = Path(stdout_path)

416

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

418 results = {}

419

420 energy = self.read_energy()

421 results['energy'] = energy

422 results['free_energy'] = energy

423

424 if read_forces:

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

426

427 if any(atoms.pbc):

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

429

430 return results

431

432 def read_forces(self, atoms):

433 forcename = self.directory / 'dftd3_gradient'

434 with open(forcename) as fd:

435 forces = self.parse_forces(fd)

436

437 assert len(forces) == len(atoms)

438

439 forces *= -Hartree / Bohr

440 # XXXX ordering!

441 if any(atoms.pbc):

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

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

444 # forces!

445 ind = np.argsort(atoms.symbols)

446 forces[ind] = forces.copy()

447 return forces

448

449 def read_stress(self, cell):

450 volume = cell.volume

451 assert volume > 0

452

453 stress = self.read_cellgradient()

454 stress *= Hartree / Bohr / volume

455 stress = stress.T @ cell

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

457

458 def read_cellgradient(self):

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

460 return self.parse_cellgradient(fd)

461

462 def read_energy(self) -> float:

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

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

465

466 def parse_energy(self, fd, outname):

467 for line in fd:

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

469 if 'functional name unknown' in line:

470 message = ('Unknown DFTD3 functional name. '

471 'Please check the dftd3.f source file '

472 'for the list of known functionals '

473 'and their spelling.')

474 else:

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

476 'output file and report any errors '

477 'to the ASE developers.'

478 ''.format(outname))

479 raise RuntimeError(message)

480

481 if line.startswith(' Edisp'):

482 # line looks something like this:

483 #

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

485 #

486 parts = line.split()

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

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

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

490 return e_dftd3

491

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

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

494

495 def parse_forces(self, fd):

496 forces = []

497 for i, line in enumerate(fd):

498 forces.append(line.split())

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

500

501 def parse_cellgradient(self, fd):

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

503 for i, line in enumerate(fd):

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

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

506 # Check if all stress elements are present?

507 # Check if file is longer?

508 return stress

509

510

511def _get_damppars(par):

512 damping = par['damping']

513

514 damppars = []

515

516 # s6 is always first

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

518

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

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

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

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

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

524

525 # s8 is always third

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

527

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

529 if damping == 'zero':

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

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

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

533 elif damping == 'zerom':

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

535 # alpha6 is always fifth

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

537

538 # last is the version number

539 if par['old']:

540 damppars.append('2')

541 elif damping == 'zero':

542 damppars.append('3')

543 elif damping == 'bj':

544 damppars.append('4')

545 elif damping == 'zerom':

546 damppars.append('5')

547 elif damping == 'bjm':

548 damppars.append('6')

549 return damppars