1import os 

2import subprocess 

3from warnings import warn 

4 

5import numpy as np 

6from ase.calculators.calculator import (Calculator, FileIOCalculator, 

7 PropertyNotImplementedError, 

8 all_changes) 

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 

15class DFTD3(FileIOCalculator): 

16 """Grimme DFT-D3 calculator""" 

17 

18 name = 'DFTD3' 

19 command = 'dftd3' 

20 dftd3_implemented_properties = ['energy', 'forces', 'stress'] 

21 

22 damping_methods = ['zero', 'bj', 'zerom', 'bjm'] 

23 

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

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

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

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

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

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

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

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

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

33 'sr6': None, 

34 's8': None, 

35 'sr8': None, 

36 'alpha6': None, 

37 'a1': None, 

38 'a2': None, 

39 'beta': None} 

40 

41 dftd3_flags = ('grad', 'pbc', 'abc', 'old', 'tz') 

42 

43 def __init__(self, 

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

45 command=None, # Command for running dftd3 

46 dft=None, # DFT calculator 

47 atoms=None, 

48 comm=world, 

49 **kwargs): 

50 

51 self.dft = None 

52 FileIOCalculator.__init__(self, restart=None, 

53 label=label, 

54 atoms=atoms, 

55 command=command, 

56 dft=dft, 

57 **kwargs) 

58 

59 self.comm = comm 

60 

61 def set(self, **kwargs): 

62 changed_parameters = {} 

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

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

65 # CLI dftd3 interface. 

66 if kwargs.get('func'): 

67 if kwargs.get('xc') and kwargs['func'] != kwargs['xc']: 

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

69 'Please provide at most one of these ' 

70 'two keywords. The preferred keyword ' 

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

72 'consistency with the CLI dftd3 ' 

73 'interface.') 

74 if kwargs['func'] != self.parameters['xc']: 

75 changed_parameters['xc'] = kwargs['func'] 

76 self.parameters['xc'] = kwargs['func'] 

77 

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

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

80 # implements more properties, we will expose those properties too. 

81 if 'dft' in kwargs: 

82 dft = kwargs.pop('dft') 

83 if dft is not self.dft: 

84 changed_parameters['dft'] = dft 

85 if dft is None: 

86 self.implemented_properties = self.dftd3_implemented_properties 

87 else: 

88 self.implemented_properties = dft.implemented_properties 

89 self.dft = dft 

90 

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

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

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

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

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

96 if self.parameters['xc'] is None and self.dft is not None: 

97 if self.dft.parameters.get('xc'): 

98 self.parameters['xc'] = self.dft.parameters['xc'] 

99 

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

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

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

103 if unknown_kwargs: 

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

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

106 

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

108 

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

110 if self.parameters['damping'] is not None: 

111 self.parameters['damping'] = self.parameters['damping'].lower() 

112 if self.parameters['damping'] not in self.damping_methods: 

113 raise ValueError('Unknown damping method {}!' 

114 ''.format(self.parameters['damping'])) 

115 

116 # d2 only is valid with 'zero' damping 

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

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

119 'dispersion correction method!') 

120 

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

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

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

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

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

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

127 'to {cutoff}.' 

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

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

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

131 

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

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

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

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

136 if not self.parameters['grad']: 

137 for val in ['forces', 'stress']: 

138 if val in self.implemented_properties: 

139 self.implemented_properties.remove(val) 

140 

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

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

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

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

145 all_damppars = zero_damppars | bj_damppars | zerom_damppars 

146 

147 self.custom_damp = False 

148 damping = self.parameters['damping'] 

149 damppars = set(kwargs) & all_damppars 

150 if damppars: 

151 self.custom_damp = True 

152 if damping == 'zero': 

153 valid_damppars = zero_damppars 

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

155 valid_damppars = bj_damppars 

156 elif damping == 'zerom': 

157 valid_damppars = zerom_damppars 

158 

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

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

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

162 # dftd3 executable & depend on XC functional. 

163 missing_damppars = valid_damppars - damppars 

164 if missing_damppars and missing_damppars != valid_damppars: 

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

166 'parameters for the {} damping method was ' 

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

168 ''.format(damping, 

169 ', '.join(valid_damppars), 

170 ', '.join(damppars))) 

171 

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

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

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

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

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

177 # worry about that here. 

178 if damppars - valid_damppars: 

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

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

181 ''.format(damping, 

182 ', '.join(damppars))) 

183 

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

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

186 # parameters. 

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

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

189 'instead of those parameterized for {}!' 

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

191 

192 if changed_parameters: 

193 self.results.clear() 

194 return changed_parameters 

195 

196 def calculate(self, atoms=None, properties=['energy'], 

197 system_changes=all_changes): 

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

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

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

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

202 

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

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

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

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

207 os.remove(localparfile) 

208 

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

210 # custom damping parameters. 

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

212 command = self._generate_command() 

213 

214 # Finally, call dftd3 and parse results. 

215 # DFTD3 does not run in parallel 

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

217 errorcode = 0 

218 if self.comm.rank == 0: 

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

220 errorcode = subprocess.call(command, 

221 cwd=self.directory, stdout=fd) 

222 

223 errorcode = self.comm.sum(errorcode) 

224 

225 if errorcode: 

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

227 (self.name, errorcode)) 

228 

229 self.read_results() 

230 

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

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

233 system_changes=system_changes) 

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

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

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

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

238 # the user. 

239 pbc = False 

240 if any(atoms.pbc): 

241 if not all(atoms.pbc): 

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

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

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

245 pbc = True 

246 

247 if self.comm.rank == 0: 

248 if pbc: 

249 fname = os.path.join(self.directory, 

250 '{}.POSCAR'.format(self.label)) 

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

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

253 # atomtypes 

254 write_vasp(fname, atoms, sort=True) 

255 else: 

256 fname = os.path.join( 

257 self.directory, '{}.xyz'.format(self.label)) 

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

259 

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

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

262 if self.custom_damp: 

263 damppars = [] 

264 # s6 is always first 

265 damppars.append(str(float(self.parameters['s6']))) 

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

267 if self.parameters['damping'] in ['zero', 'zerom']: 

268 damppars.append(str(float(self.parameters['sr6']))) 

269 elif self.parameters['damping'] in ['bj', 'bjm']: 

270 damppars.append(str(float(self.parameters['a1']))) 

271 # s8 is always third 

272 damppars.append(str(float(self.parameters['s8']))) 

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

274 if self.parameters['damping'] == 'zero': 

275 damppars.append(str(float(self.parameters['sr8']))) 

276 elif self.parameters['damping'] in ['bj', 'bjm']: 

277 damppars.append(str(float(self.parameters['a2']))) 

278 elif self.parameters['damping'] == 'zerom': 

279 damppars.append(str(float(self.parameters['beta']))) 

280 # alpha6 is always fifth 

281 damppars.append(str(int(self.parameters['alpha6']))) 

282 # last is the version number 

283 if self.parameters['old']: 

284 damppars.append('2') 

285 elif self.parameters['damping'] == 'zero': 

286 damppars.append('3') 

287 elif self.parameters['damping'] == 'bj': 

288 damppars.append('4') 

289 elif self.parameters['damping'] == 'zerom': 

290 damppars.append('5') 

291 elif self.parameters['damping'] == 'bjm': 

292 damppars.append('6') 

293 

294 damp_fname = os.path.join(self.directory, '.dftd3par.local') 

295 if self.comm.rank == 0: 

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

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

298 

299 def read_results(self): 

300 # parse the energy 

301 outname = os.path.join(self.directory, self.label + '.out') 

302 energy = 0.0 

303 if self.comm.rank == 0: 

304 with open(outname, 'r') as fd: 

305 for line in fd: 

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

307 if 'functional name unknown' in line: 

308 message = 'Unknown DFTD3 functional name "{}". ' \ 

309 'Please check the dftd3.f source file ' \ 

310 'for the list of known functionals ' \ 

311 'and their spelling.' \ 

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

313 else: 

314 message = 'dftd3 failed! Please check the {} ' \ 

315 'output file and report any errors ' \ 

316 'to the ASE developers.' \ 

317 ''.format(outname) 

318 raise RuntimeError(message) 

319 

320 if line.startswith(' Edisp'): 

321 # line looks something like this: 

322 # 

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

324 # 

325 parts = line.split() 

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

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

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

329 energy = e_dftd3 

330 break 

331 else: 

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

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

334 

335 self.results['energy'] = self.comm.sum(energy) 

336 self.results['free_energy'] = self.results['energy'] 

337 

338 # FIXME: Calculator.get_potential_energy() simply inspects 

339 # self.results for the free energy rather than calling 

340 # Calculator.get_property('free_energy'). For example, GPAW does 

341 # not actually present free_energy as an implemented property, even 

342 # though it does calculate it. So, we are going to add in the DFT 

343 # free energy to our own results if it is present in the attached 

344 # calculator. TODO: Fix the Calculator interface!!! 

345 if self.dft is not None: 

346 try: 

347 efree = self.dft.get_potential_energy( 

348 force_consistent=True) 

349 self.results['free_energy'] += efree 

350 except PropertyNotImplementedError: 

351 pass 

352 

353 if self.parameters['grad']: 

354 # parse the forces 

355 forces = np.zeros((len(self.atoms), 3)) 

356 forcename = os.path.join(self.directory, 'dftd3_gradient') 

357 if self.comm.rank == 0: 

358 with open(forcename, 'r') as fd: 

359 for i, line in enumerate(fd): 

360 forces[i] = np.array([float(x) for x in line.split()]) 

361 forces *= -Hartree / Bohr 

362 self.comm.broadcast(forces, 0) 

363 if self.atoms.pbc.any(): 

364 ind = np.argsort(self.atoms.get_chemical_symbols()) 

365 forces[ind] = forces.copy() 

366 self.results['forces'] = forces 

367 

368 if any(self.atoms.pbc): 

369 # parse the stress tensor 

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

371 stressname = os.path.join(self.directory, 'dftd3_cellgradient') 

372 if self.comm.rank == 0: 

373 with open(stressname, 'r') as fd: 

374 for i, line in enumerate(fd): 

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

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

377 stress *= Hartree / Bohr / self.atoms.get_volume() 

378 stress = np.dot(stress.T, self.atoms.cell) 

379 self.comm.broadcast(stress, 0) 

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

381 

382 def get_property(self, name, atoms=None, allow_calculation=True): 

383 dft_result = None 

384 if self.dft is not None: 

385 dft_result = self.dft.get_property(name, atoms, allow_calculation) 

386 

387 dftd3_result = FileIOCalculator.get_property(self, name, atoms, 

388 allow_calculation) 

389 

390 if dft_result is None and dftd3_result is None: 

391 return None 

392 elif dft_result is None: 

393 return dftd3_result 

394 elif dftd3_result is None: 

395 return dft_result 

396 else: 

397 return dft_result + dftd3_result 

398 

399 def _generate_command(self): 

400 command = self.command.split() 

401 

402 if any(self.atoms.pbc): 

403 command.append(self.label + '.POSCAR') 

404 else: 

405 command.append(self.label + '.xyz') 

406 

407 if not self.custom_damp: 

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

409 if xc is None: 

410 xc = 'pbe' 

411 command += ['-func', xc.lower()] 

412 

413 for arg in self.dftd3_flags: 

414 if self.parameters.get(arg): 

415 command.append('-' + arg) 

416 

417 if any(self.atoms.pbc): 

418 command.append('-pbc') 

419 

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

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

422 

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

424 command.append('-' + self.parameters['damping']) 

425 

426 return command