1"""van der Waals correction schemes for DFT""" 

2import numpy as np 

3from ase.units import Bohr, Hartree 

4from ase.calculators.calculator import Calculator 

5from scipy.special import erfinv, erfc 

6from ase.neighborlist import neighbor_list 

7from ase.parallel import world 

8from ase.utils import IOContext 

9 

10 

11# dipole polarizabilities and C6 values from 

12# X. Chu and A. Dalgarno, J. Chem. Phys. 121 (2004) 4083 

13# atomic units, a_0^3 

14vdWDB_Chu04jcp = { 

15 # Element: [alpha, C6]; units [Bohr^3, Hartree * Bohr^6] 

16 'H': [4.5, 6.5], # [exact, Tkatchenko PRL] 

17 'He': [1.38, 1.42], 

18 'Li': [164, 1392], 

19 'Be': [38, 227], 

20 'B': [21, 99.5], 

21 'C': [12, 46.6], 

22 'N': [7.4, 24.2], 

23 'O': [5.4, 15.6], 

24 'F': [3.8, 9.52], 

25 'Ne': [2.67, 6.20], 

26 'Na': [163, 1518], 

27 'Mg': [71, 626], 

28 'Al': [60, 528], 

29 'Si': [37, 305], 

30 'P': [25, 185], 

31 'S': [19.6, 134], 

32 'Cl': [15, 94.6], 

33 'Ar': [11.1, 64.2], 

34 'Ca': [160, 2163], 

35 'Sc': [120, 1383], 

36 'Ti': [98, 1044], 

37 'V': [84, 832], 

38 'Cr': [78, 602], 

39 'Mn': [63, 552], 

40 'Fe': [56, 482], 

41 'Co': [50, 408], 

42 'Ni': [48, 373], 

43 'Cu': [42, 253], 

44 'Zn': [40, 284], 

45 'As': [29, 246], 

46 'Se': [25, 210], 

47 'Br': [20, 162], 

48 'Kr': [16.7, 130], 

49 'Sr': [199, 3175], 

50 'Te': [40, 445], 

51 'I': [35, 385]} 

52 

53vdWDB_alphaC6 = vdWDB_Chu04jcp 

54 

55# dipole polarizabilities and C6 values from 

56# V. G. Ruiz et al. Phys. Rev. Lett 108 (2012) 146103 

57# atomic units, a_0^3 

58vdWDB_Ruiz12prl = { 

59 'Ag': [50.6, 339], 

60 'Au': [36.5, 298], 

61 'Pd': [23.7, 158], 

62 'Pt': [39.7, 347], 

63} 

64 

65vdWDB_alphaC6.update(vdWDB_Ruiz12prl) 

66 

67# C6 values and vdW radii from 

68# S. Grimme, J Comput Chem 27 (2006) 1787-1799 

69vdWDB_Grimme06jcc = { 

70 # Element: [C6, R0]; units [J nm^6 mol^{-1}, Angstrom] 

71 'H': [0.14, 1.001], 

72 'He': [0.08, 1.012], 

73 'Li': [1.61, 0.825], 

74 'Be': [1.61, 1.408], 

75 'B': [3.13, 1.485], 

76 'C': [1.75, 1.452], 

77 'N': [1.23, 1.397], 

78 'O': [0.70, 1.342], 

79 'F': [0.75, 1.287], 

80 'Ne': [0.63, 1.243], 

81 'Na': [5.71, 1.144], 

82 'Mg': [5.71, 1.364], 

83 'Al': [10.79, 1.639], 

84 'Si': [9.23, 1.716], 

85 'P': [7.84, 1.705], 

86 'S': [5.57, 1.683], 

87 'Cl': [5.07, 1.639], 

88 'Ar': [4.61, 1.595], 

89 'K': [10.80, 1.485], 

90 'Ca': [10.80, 1.474], 

91 'Sc': [10.80, 1.562], 

92 'Ti': [10.80, 1.562], 

93 'V': [10.80, 1.562], 

94 'Cr': [10.80, 1.562], 

95 'Mn': [10.80, 1.562], 

96 'Fe': [10.80, 1.562], 

97 'Co': [10.80, 1.562], 

98 'Ni': [10.80, 1.562], 

99 'Cu': [10.80, 1.562], 

100 'Zn': [10.80, 1.562], 

101 'Ga': [16.99, 1.650], 

102 'Ge': [17.10, 1.727], 

103 'As': [16.37, 1.760], 

104 'Se': [12.64, 1.771], 

105 'Br': [12.47, 1.749], 

106 'Kr': [12.01, 1.727], 

107 'Rb': [24.67, 1.628], 

108 'Sr': [24.67, 1.606], 

109 'Y-Cd': [24.67, 1.639], 

110 'In': [37.32, 1.672], 

111 'Sn': [38.71, 1.804], 

112 'Sb': [38.44, 1.881], 

113 'Te': [31.74, 1.892], 

114 'I': [31.50, 1.892], 

115 'Xe': [29.99, 1.881]} 

116 

117 

118# Optimal range parameters sR for different XC functionals 

119# to be used with the Tkatchenko-Scheffler scheme 

120# Reference: M.A. Caro arXiv:1704.00761 (2017) 

121sR_opt = {'PBE': 0.940, 

122 'RPBE': 0.590, 

123 'revPBE': 0.585, 

124 'PBEsol': 1.055, 

125 'BLYP': 0.625, 

126 'AM05': 0.840, 

127 'PW91': 0.965} 

128 

129 

130def get_logging_file_descriptor(calculator): 

131 if hasattr(calculator, 'log'): 

132 fd = calculator.log 

133 if hasattr(fd, 'write'): 

134 return fd 

135 if hasattr(fd, 'fd'): 

136 return fd.fd 

137 if hasattr(calculator, 'txt'): 

138 return calculator.txt 

139 

140 

141class vdWTkatchenko09prl(Calculator, IOContext): 

142 """vdW correction after Tkatchenko and Scheffler PRL 102 (2009) 073005.""" 

143 def __init__(self, 

144 hirshfeld=None, vdwradii=None, calculator=None, 

145 Rmax=10., # maximal radius for periodic calculations 

146 Ldecay=1., # decay length for smoothing in periodic calcs 

147 vdWDB_alphaC6=vdWDB_alphaC6, 

148 txt=None, sR=None): 

149 """Constructor 

150 

151 Parameters 

152 ========== 

153 hirshfeld: the Hirshfeld partitioning object 

154 calculator: the calculator to get the PBE energy 

155 """ 

156 self.hirshfeld = hirshfeld 

157 if calculator is None: 

158 self.calculator = self.hirshfeld.get_calculator() 

159 else: 

160 self.calculator = calculator 

161 

162 if txt is None: 

163 txt = get_logging_file_descriptor(self.calculator) 

164 if hasattr(self.calculator, 'world'): 

165 myworld = self.calculator.world 

166 else: 

167 myworld = world # the best we know 

168 self.txt = self.openfile(txt, myworld) 

169 

170 self.vdwradii = vdwradii 

171 self.vdWDB_alphaC6 = vdWDB_alphaC6 

172 self.Rmax = Rmax 

173 self.Ldecay = Ldecay 

174 self.atoms = None 

175 

176 if sR is None: 

177 try: 

178 xc_name = self.calculator.get_xc_functional() 

179 self.sR = sR_opt[xc_name] 

180 except KeyError: 

181 raise ValueError( 

182 'Tkatchenko-Scheffler dispersion correction not ' + 

183 'implemented for %s functional' % xc_name) 

184 else: 

185 self.sR = sR 

186 self.d = 20 

187 

188 Calculator.__init__(self) 

189 

190 self.parameters['calculator'] = self.calculator.name 

191 self.parameters['xc'] = self.calculator.get_xc_functional() 

192 

193 @property 

194 def implemented_properties(self): 

195 return self.calculator.implemented_properties 

196 

197 def calculation_required(self, atoms, quantities): 

198 if self.calculator.calculation_required(atoms, quantities): 

199 return True 

200 for quantity in quantities: 

201 if quantity not in self.results: 

202 return True 

203 return False 

204 

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

206 system_changes=[]): 

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

208 self.update(atoms, properties) 

209 

210 def update(self, atoms=None, properties=['energy', 'forces']): 

211 if not self.calculation_required(atoms, properties): 

212 return 

213 

214 if atoms is None: 

215 atoms = self.calculator.get_atoms() 

216 

217 properties = list(properties) 

218 for name in 'energy', 'forces': 

219 if name not in properties: 

220 properties.append(name) 

221 

222 for name in properties: 

223 self.results[name] = self.calculator.get_property(name, atoms) 

224 self.parameters['uncorrected_energy'] = self.results['energy'] 

225 self.atoms = atoms.copy() 

226 

227 if self.vdwradii is not None: 

228 # external vdW radii 

229 vdwradii = self.vdwradii 

230 assert(len(atoms) == len(vdwradii)) 

231 else: 

232 vdwradii = [] 

233 for atom in atoms: 

234 self.vdwradii.append(vdWDB_Grimme06jcc[atom.symbol][1]) 

235 

236 if self.hirshfeld is None: 

237 volume_ratios = [1.] * len(atoms) 

238 elif hasattr(self.hirshfeld, '__len__'): # a list 

239 assert(len(atoms) == len(self.hirshfeld)) 

240 volume_ratios = self.hirshfeld 

241 else: # should be an object 

242 self.hirshfeld.initialize() 

243 volume_ratios = self.hirshfeld.get_effective_volume_ratios() 

244 

245 # correction for effective C6 

246 na = len(atoms) 

247 C6eff_a = np.empty((na)) 

248 alpha_a = np.empty((na)) 

249 R0eff_a = np.empty((na)) 

250 for a, atom in enumerate(atoms): 

251 # free atom values 

252 alpha_a[a], C6eff_a[a] = self.vdWDB_alphaC6[atom.symbol] 

253 # correction for effective C6 

254 C6eff_a[a] *= Hartree * volume_ratios[a]**2 * Bohr**6 

255 R0eff_a[a] = vdwradii[a] * volume_ratios[a]**(1 / 3.) 

256 C6eff_aa = np.empty((na, na)) 

257 for a in range(na): 

258 for b in range(a, na): 

259 C6eff_aa[a, b] = (2 * C6eff_a[a] * C6eff_a[b] / 

260 (alpha_a[b] / alpha_a[a] * C6eff_a[a] + 

261 alpha_a[a] / alpha_a[b] * C6eff_a[b])) 

262 C6eff_aa[b, a] = C6eff_aa[a, b] 

263 

264 # New implementation by Miguel Caro 

265 # (complaints etc to mcaroba@gmail.com) 

266 # If all 3 PBC are False, we do the summation over the atom 

267 # pairs in the simulation box. If any of them is True, we 

268 # use the cutoff radius instead 

269 pbc_c = atoms.get_pbc() 

270 EvdW = 0.0 

271 forces = 0. * self.results['forces'] 

272 # PBC: we build a neighbor list according to the Reff criterion 

273 if pbc_c.any(): 

274 # Effective cutoff radius 

275 tol = 1.e-5 

276 Reff = self.Rmax + self.Ldecay * erfinv(1. - 2. * tol) 

277 # Build list of neighbors 

278 n_list = neighbor_list(quantities="ijdDS", 

279 a=atoms, 

280 cutoff=Reff, 

281 self_interaction=False) 

282 atom_list = [[] for _ in range(0, len(atoms))] 

283 d_list = [[] for _ in range(0, len(atoms))] 

284 v_list = [[] for _ in range(0, len(atoms))] 

285 # r_list = [[] for _ in range(0, len(atoms))] 

286 # Look for neighbor pairs 

287 for k in range(0, len(n_list[0])): 

288 i = n_list[0][k] 

289 j = n_list[1][k] 

290 dist = n_list[2][k] 

291 vect = n_list[3][k] # vect is the distance rj - ri 

292 # repl = n_list[4][k] 

293 if j >= i: 

294 atom_list[i].append(j) 

295 d_list[i].append(dist) 

296 v_list[i].append(vect) 

297 # r_list[i].append( repl ) 

298 # Not PBC: we loop over all atoms in the unit cell only 

299 else: 

300 atom_list = [] 

301 d_list = [] 

302 v_list = [] 

303 # r_list = [] 

304 # Do this to avoid double counting 

305 for i in range(0, len(atoms)): 

306 atom_list.append(range(i + 1, len(atoms))) 

307 d_list.append([atoms.get_distance(i, j) 

308 for j in range(i + 1, len(atoms))]) 

309 v_list.append([atoms.get_distance(i, j, vector=True) 

310 for j in range(i + 1, len(atoms))]) 

311 # r_list.append( [[0,0,0] for j in range(i+1, len(atoms))]) 

312 # No PBC means we are in the same cell 

313 # Here goes the calculation, valid with and without 

314 # PBC because we loop over 

315 # independent pairwise *interactions* 

316 for i in range(len(atoms)): 

317 # for j, r, vect, repl in zip(atom_list[i], d_list[i], 

318 # v_list[i], r_list[i]): 

319 for j, r, vect in zip(atom_list[i], d_list[i], v_list[i]): 

320 r6 = r**6 

321 Edamp, Fdamp = self.damping(r, 

322 R0eff_a[i], 

323 R0eff_a[j], 

324 d=self.d, 

325 sR=self.sR) 

326 if pbc_c.any(): 

327 smooth = 0.5 * erfc((r - self.Rmax) / self.Ldecay) 

328 smooth_der = -1. / np.sqrt(np.pi) / self.Ldecay * np.exp( 

329 -((r - self.Rmax) / self.Ldecay)**2) 

330 else: 

331 smooth = 1. 

332 smooth_der = 0. 

333 # Here we compute the contribution to the energy 

334 # Self interactions (only possible in PBC) are double counted. 

335 # We correct it here 

336 if i == j: 

337 EvdW -= (Edamp * C6eff_aa[i, j] / r6) / 2. * smooth 

338 else: 

339 EvdW -= (Edamp * C6eff_aa[i, j] / r6) * smooth 

340 # Here we compute the contribution to the forces 

341 # We neglect the C6eff contribution to the forces 

342 # (which can actually be larger 

343 # than the other contributions) 

344 # Self interactions do not contribute to the forces 

345 if i != j: 

346 # Force on i due to j 

347 force_ij = -( 

348 (Fdamp - 6 * Edamp / r) * C6eff_aa[i, j] / r6 * smooth 

349 + (Edamp * C6eff_aa[i, j] / r6) * smooth_der) * vect / r 

350 # Forces go both ways for every interaction 

351 forces[i] += force_ij 

352 forces[j] -= force_ij 

353 self.results['energy'] += EvdW 

354 self.results['forces'] += forces 

355 

356 if self.txt: 

357 print(('\n' + self.__class__.__name__), file=self.txt) 

358 print('vdW correction: %g' % (EvdW), file=self.txt) 

359 print('Energy: %g' % self.results['energy'], 

360 file=self.txt) 

361 print('\nForces in eV/Ang:', file=self.txt) 

362 symbols = self.atoms.get_chemical_symbols() 

363 for ia, symbol in enumerate(symbols): 

364 print('%3d %-2s %10.5f %10.5f %10.5f' % 

365 ((ia, symbol) + tuple(self.results['forces'][ia])), 

366 file=self.txt) 

367 self.txt.flush() 

368 

369 def damping(self, RAB, R0A, R0B, 

370 d=20, # steepness of the step function for PBE 

371 sR=0.94): 

372 """Damping factor. 

373 

374 Standard values for d and sR as given in 

375 Tkatchenko and Scheffler PRL 102 (2009) 073005.""" 

376 scale = 1.0 / (sR * (R0A + R0B)) 

377 x = RAB * scale 

378 chi = np.exp(-d * (x - 1.0)) 

379 return 1.0 / (1.0 + chi), d * scale * chi / (1.0 + chi)**2