1import os 

2import subprocess 

3from warnings import warn 

4from pathlib import Path 

5 

6import numpy as np 

7from ase.calculators.calculator import (BaseCalculator, FileIOCalculator, 

8 Calculator) 

9from ase.io import write 

10from ase.io.vasp import write_vasp 

11from ase.parallel import world 

12from ase.units import Bohr, Hartree 

13 

14 

15def dftd3_defaults(): 

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

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

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

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

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

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

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

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

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

25 'sr6': None, 

26 's8': None, 

27 'sr8': None, 

28 'alpha6': None, 

29 'a1': None, 

30 'a2': None, 

31 'beta': None} 

32 return default_parameters 

33 

34 

35class DFTD3(BaseCalculator): 

36 """Grimme DFT-D3 calculator""" 

37 

38 def __init__(self, 

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

40 command=None, # Command for running dftd3 

41 dft=None, # DFT calculator 

42 comm=world, 

43 **kwargs): 

44 

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

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

47 # CLI dftd3 interface. 

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

49 if func is not None: 

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

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

52 'Please provide at most one of these ' 

53 'two keywords. The preferred keyword ' 

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

55 'consistency with the CLI dftd3 ' 

56 'interface.') 

57 kwargs['xc'] = func 

58 

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

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

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

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

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

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

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

66 if dft_xc is not None: 

67 kwargs['xc'] = dft_xc 

68 

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

70 

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

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

73 # implements more properties, we expose those properties. 

74 # dftd3 contributions for those properties will be zero. 

75 if dft is None: 

76 self.implemented_properties = list(dftd3.dftd3_properties) 

77 else: 

78 self.implemented_properties = list(dft.implemented_properties) 

79 

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

81 # or are they not really "parameters"? 

82 # 

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

84 super().__init__() 

85 

86 self.dftd3 = dftd3 

87 self.dft = dft 

88 

89 def todict(self): 

90 return {} 

91 

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

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

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

95 

96 if self.dft is None: 

97 results = dftd3_results 

98 else: 

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

100 results = dict(dft_results) 

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

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

103 results[name] += dftd3_results[name] 

104 

105 # Although DFTD3 may have calculated quantities not provided 

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

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

108 assert set(results) == set(dft_results) 

109 self.results = results 

110 

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

112 # We want any and all properties that the calculator 

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

114 # dictionary instead of only getting the requested properties. 

115 

116 import copy 

117 for name in properties: 

118 calc.get_property(name, atoms) 

119 assert name in calc.results 

120 

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

122 results = copy.deepcopy(calc.results) 

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

124 return results 

125 

126 

127class PureDFTD3(FileIOCalculator): 

128 """DFTD3 calculator without corresponding DFT contribution. 

129 

130 This class is an implementation detail.""" 

131 

132 name = 'puredftd3' 

133 command = 'dftd3' 

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 

140 def __init__(self, 

141 *, 

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

143 command=None, # Command for running dftd3 

144 comm=world, 

145 **kwargs): 

146 

147 super().__init__(label=label, 

148 command=command, 

149 **kwargs) 

150 

151 self.comm = comm 

152 

153 def set(self, **kwargs): 

154 changed_parameters = {} 

155 

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

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

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

159 if unknown_kwargs: 

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

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

162 

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

164 

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

166 damping = self.parameters['damping'] 

167 if damping is not None: 

168 damping = damping.lower() 

169 if damping not in self.damping_methods: 

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

171 

172 # d2 only is valid with 'zero' damping 

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

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

175 'dispersion correction method!') 

176 

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

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

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

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

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

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

183 'to {cutoff}.' 

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

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

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

187 

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

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

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

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

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

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

194 # if val in self.implemented_properties: 

195 # self.implemented_properties.remove(val) 

196 

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

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

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

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

201 all_damppars = zero_damppars | bj_damppars | zerom_damppars 

202 

203 self.custom_damp = False 

204 

205 damppars = set(kwargs) & all_damppars 

206 if damppars: 

207 self.custom_damp = True 

208 if damping == 'zero': 

209 valid_damppars = zero_damppars 

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

211 valid_damppars = bj_damppars 

212 elif damping == 'zerom': 

213 valid_damppars = zerom_damppars 

214 

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

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

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

218 # dftd3 executable & depend on XC functional. 

219 missing_damppars = valid_damppars - damppars 

220 if missing_damppars and missing_damppars != valid_damppars: 

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

222 'parameters for the {} damping method was ' 

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

224 ''.format(damping, 

225 ', '.join(valid_damppars), 

226 ', '.join(damppars))) 

227 

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

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

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

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

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

233 # worry about that here. 

234 if damppars - valid_damppars: 

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

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

237 ''.format(damping, 

238 ', '.join(damppars))) 

239 

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

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

242 # parameters. 

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

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

245 'instead of those parameterized for {}!' 

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

247 

248 if changed_parameters: 

249 self.results.clear() 

250 return changed_parameters 

251 

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

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

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

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

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

257 

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

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

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

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

262 os.remove(localparfile) 

263 

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

265 # custom damping parameters. 

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

267 # command = self._generate_command() 

268 

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

270 atoms=self.atoms, parameters=self.parameters) 

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

272 

273 # Finally, call dftd3 and parse results. 

274 # DFTD3 does not run in parallel 

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

276 errorcode = 0 

277 if self.comm.rank == 0: 

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

279 errorcode = subprocess.call(command, 

280 cwd=self.directory, stdout=fd) 

281 

282 errorcode = self.comm.sum(errorcode) 

283 

284 if errorcode: 

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

286 (self.name, errorcode)) 

287 

288 self.read_results() 

289 

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

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

292 system_changes=system_changes) 

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

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

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

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

297 # the user. 

298 

299 if self.custom_damp: 

300 damppars = _get_damppars(self.parameters) 

301 else: 

302 damppars = None 

303 

304 pbc = any(atoms.pbc) 

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

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

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

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

309 

310 if self.comm.rank == 0: 

311 self._actually_write_input( 

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

313 properties=properties, prefix=self.label, 

314 damppars=damppars, pbc=pbc) 

315 

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

317 damppars, pbc): 

318 if pbc: 

319 fname = directory / '{}.POSCAR'.format(prefix) 

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

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

322 # atomtypes 

323 write_vasp(fname, atoms, sort=True) 

324 else: 

325 fname = directory / '{}.xyz'.format(prefix) 

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

327 

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

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

330 if damppars is not None: 

331 damp_fname = directory / '.dftd3par.local' 

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

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

334 

335 def _outname(self): 

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

337 

338 def _read_and_broadcast_results(self): 

339 from ase.parallel import broadcast 

340 if self.comm.rank == 0: 

341 output = DFTD3Output(directory=self.directory, 

342 stdout_path=self._outname()) 

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

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

345 else: 

346 dct = None 

347 

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

349 return dct 

350 

351 def read_results(self): 

352 results = self._read_and_broadcast_results() 

353 self.results = results 

354 

355 

356class DFTD3Inputs: 

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

358 

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

360 self.command = command 

361 self.prefix = prefix 

362 self.atoms = atoms 

363 self.parameters = parameters 

364 

365 @property 

366 def pbc(self): 

367 return any(self.atoms.pbc) 

368 

369 @property 

370 def inputformat(self): 

371 if self.pbc: 

372 return 'POSCAR' 

373 else: 

374 return 'xyz' 

375 

376 def get_argv(self, custom_damp): 

377 argv = self.command.split() 

378 

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

380 

381 if not custom_damp: 

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

383 if xc is None: 

384 xc = 'pbe' 

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

386 

387 for arg in self.dftd3_flags: 

388 if self.parameters.get(arg): 

389 argv.append('-' + arg) 

390 

391 if self.pbc: 

392 argv.append('-pbc') 

393 

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

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

396 

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

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

399 

400 return argv 

401 

402 

403class DFTD3Output: 

404 def __init__(self, directory, stdout_path): 

405 self.directory = Path(directory) 

406 self.stdout_path = Path(stdout_path) 

407 

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

409 results = {} 

410 

411 energy = self.read_energy() 

412 results['energy'] = energy 

413 results['free_energy'] = energy 

414 

415 if read_forces: 

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

417 

418 if any(atoms.pbc): 

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

420 

421 return results 

422 

423 def read_forces(self, atoms): 

424 forcename = self.directory / 'dftd3_gradient' 

425 with open(forcename) as fd: 

426 forces = self.parse_forces(fd) 

427 

428 assert len(forces) == len(atoms) 

429 

430 forces *= -Hartree / Bohr 

431 # XXXX ordering! 

432 if any(atoms.pbc): 

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

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

435 # forces! 

436 ind = np.argsort(atoms.symbols) 

437 forces[ind] = forces.copy() 

438 return forces 

439 

440 def read_stress(self, cell): 

441 volume = cell.volume 

442 assert volume > 0 

443 

444 stress = self.read_cellgradient() 

445 stress *= Hartree / Bohr / volume 

446 stress = stress.T @ cell 

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

448 

449 def read_cellgradient(self): 

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

451 return self.parse_cellgradient(fd) 

452 

453 def read_energy(self) -> float: 

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

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

456 

457 def parse_energy(self, fd, outname): 

458 for line in fd: 

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

460 if 'functional name unknown' in line: 

461 message = ('Unknown DFTD3 functional name. ' 

462 'Please check the dftd3.f source file ' 

463 'for the list of known functionals ' 

464 'and their spelling.') 

465 else: 

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

467 'output file and report any errors ' 

468 'to the ASE developers.' 

469 ''.format(outname)) 

470 raise RuntimeError(message) 

471 

472 if line.startswith(' Edisp'): 

473 # line looks something like this: 

474 # 

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

476 # 

477 parts = line.split() 

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

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

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

481 return e_dftd3 

482 

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

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

485 

486 def parse_forces(self, fd): 

487 forces = [] 

488 for i, line in enumerate(fd): 

489 forces.append(line.split()) 

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

491 

492 def parse_cellgradient(self, fd): 

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

494 for i, line in enumerate(fd): 

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

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

497 # Check if all stress elements are present? 

498 # Check if file is longer? 

499 return stress 

500 

501 

502def _get_damppars(par): 

503 damping = par['damping'] 

504 

505 damppars = [] 

506 

507 # s6 is always first 

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

509 

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

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

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

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

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

515 

516 # s8 is always third 

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

518 

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

520 if damping == 'zero': 

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

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

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

524 elif damping == 'zerom': 

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

526 # alpha6 is always fifth 

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

528 

529 # last is the version number 

530 if par['old']: 

531 damppars.append('2') 

532 elif damping == 'zero': 

533 damppars.append('3') 

534 elif damping == 'bj': 

535 damppars.append('4') 

536 elif damping == 'zerom': 

537 damppars.append('5') 

538 elif damping == 'bjm': 

539 damppars.append('6') 

540 return damppars