r m x p   toggle line displays

j k   next/prev highlighted chunk

0   (zero) top of page

1   (one) first highlighted chunk

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 ase.calculators.polarizability import StaticPolarizabilityCalculator

6from scipy.special import erfinv, erfc

7from ase.neighborlist import neighbor_list

8from ase.parallel import world, myslice

9from ase.utils import IOContext

12# dipole polarizabilities and C6 values from

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

14# atomic units, a_0^3

15vdWDB_Chu04jcp = {

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

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

18 'He': [1.38, 1.42],

19 'Li': [164, 1392],

20 'Be': [38, 227],

21 'B': [21, 99.5],

22 'C': [12, 46.6],

23 'N': [7.4, 24.2],

24 'O': [5.4, 15.6],

25 'F': [3.8, 9.52],

26 'Ne': [2.67, 6.20],

27 'Na': [163, 1518],

28 'Mg': [71, 626],

29 'Al': [60, 528],

30 'Si': [37, 305],

31 'P': [25, 185],

32 'S': [19.6, 134],

33 'Cl': [15, 94.6],

34 'Ar': [11.1, 64.2],

35 'Ca': [160, 2163],

36 'Sc': [120, 1383],

37 'Ti': [98, 1044],

38 'V': [84, 832],

39 'Cr': [78, 602],

40 'Mn': [63, 552],

41 'Fe': [56, 482],

42 'Co': [50, 408],

43 'Ni': [48, 373],

44 'Cu': [42, 253],

45 'Zn': [40, 284],

46 'As': [29, 246],

47 'Se': [25, 210],

48 'Br': [20, 162],

49 'Kr': [16.7, 130],

50 'Sr': [199, 3175],

51 'Te': [40, 445],

52 'I': [35, 385]}

54vdWDB_alphaC6 = vdWDB_Chu04jcp

56# dipole polarizabilities and C6 values from

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

58# atomic units, a_0^3

59vdWDB_Ruiz12prl = {

60 'Ag': [50.6, 339],

61 'Au': [36.5, 298],

62 'Pd': [23.7, 158],

63 'Pt': [39.7, 347],

64}

66vdWDB_alphaC6.update(vdWDB_Ruiz12prl)

68# C6 values and vdW radii from

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

70vdWDB_Grimme06jcc = {

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

72 'H': [0.14, 1.001],

73 'He': [0.08, 1.012],

74 'Li': [1.61, 0.825],

75 'Be': [1.61, 1.408],

76 'B': [3.13, 1.485],

77 'C': [1.75, 1.452],

78 'N': [1.23, 1.397],

79 'O': [0.70, 1.342],

80 'F': [0.75, 1.287],

81 'Ne': [0.63, 1.243],

82 'Na': [5.71, 1.144],

83 'Mg': [5.71, 1.364],

84 'Al': [10.79, 1.639],

85 'Si': [9.23, 1.716],

86 'P': [7.84, 1.705],

87 'S': [5.57, 1.683],

88 'Cl': [5.07, 1.639],

89 'Ar': [4.61, 1.595],

90 'K': [10.80, 1.485],

91 'Ca': [10.80, 1.474],

92 'Sc': [10.80, 1.562],

93 'Ti': [10.80, 1.562],

94 'V': [10.80, 1.562],

95 'Cr': [10.80, 1.562],

96 'Mn': [10.80, 1.562],

97 'Fe': [10.80, 1.562],

98 'Co': [10.80, 1.562],

99 'Ni': [10.80, 1.562],

100 'Cu': [10.80, 1.562],

101 'Zn': [10.80, 1.562],

102 'Ga': [16.99, 1.650],

103 'Ge': [17.10, 1.727],

104 'As': [16.37, 1.760],

105 'Se': [12.64, 1.771],

106 'Br': [12.47, 1.749],

107 'Kr': [12.01, 1.727],

108 'Rb': [24.67, 1.628],

109 'Sr': [24.67, 1.606],

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

111 'In': [37.32, 1.672],

112 'Sn': [38.71, 1.804],

113 'Sb': [38.44, 1.881],

114 'Te': [31.74, 1.892],

115 'I': [31.50, 1.892],

116 'Xe': [29.99, 1.881]}

119# Optimal range parameters sR for different XC functionals

120# to be used with the Tkatchenko-Scheffler scheme

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

122sR_opt = {'PBE': 0.940,

123 'RPBE': 0.590,

124 'revPBE': 0.585,

125 'PBEsol': 1.055,

126 'BLYP': 0.625,

127 'AM05': 0.840,

128 'PW91': 0.965}

131def get_logging_file_descriptor(calculator):

132 if hasattr(calculator, 'log'):

133 fd = calculator.log

134 if hasattr(fd, 'write'):

135 return fd

136 if hasattr(fd, 'fd'):

137 return fd.fd

138 if hasattr(calculator, 'txt'):

139 return calculator.txt

142class vdWTkatchenko09prl(Calculator, IOContext):

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

145 def __init__(self,

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

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

149 vdWDB_alphaC6=vdWDB_alphaC6,

150 txt=None, sR=None):

151 """Constructor

153 Parameters

154 ==========

155 hirshfeld: the Hirshfeld partitioning object

156 calculator: the calculator to get the PBE energy

157 """

158 self.hirshfeld = hirshfeld

159 if calculator is None:

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

161 else:

162 self.calculator = calculator

164 if txt is None:

165 txt = get_logging_file_descriptor(self.calculator)

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

167 self.comm = self.calculator.world

168 else:

169 self.comm = world # the best we know

170 self.txt = self.openfile(txt, self.comm)

173 self.vdWDB_alphaC6 = vdWDB_alphaC6

174 self.Rmax = Rmax

175 self.Ldecay = Ldecay

176 self.atoms = None

178 if sR is None:

179 try:

180 xc_name = self.calculator.get_xc_functional()

181 self.sR = sR_opt[xc_name]

182 except KeyError:

183 raise ValueError(

184 'Tkatchenko-Scheffler dispersion correction not ' +

185 'implemented for %s functional' % xc_name)

186 else:

187 self.sR = sR

188 self.d = 20

190 Calculator.__init__(self)

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

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

195 @property

196 def implemented_properties(self):

197 return self.calculator.implemented_properties

199 def calculation_required(self, atoms, quantities):

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

201 return True

202 for quantity in quantities:

203 if quantity not in self.results:

204 return True

205 return False

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

208 system_changes=[]):

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

210 self.update(atoms, properties)

212 def update(self, atoms=None,

213 properties=['energy', 'free_energy', 'forces']):

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

215 return

217 if atoms is None:

218 atoms = self.calculator.get_atoms()

220 properties = list(properties)

221 for name in 'energy', 'free_energy', 'forces':

222 if name not in properties:

223 properties.append(name)

225 for name in properties:

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

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

228 self.atoms = atoms.copy()

230 if self.vdwradii is not None:

234 else:

236 for atom in atoms:

239 if self.hirshfeld is None:

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

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

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

243 volume_ratios = self.hirshfeld

244 else: # should be an object

245 self.hirshfeld.initialize()

246 volume_ratios = self.hirshfeld.get_effective_volume_ratios()

248 # correction for effective C6

249 na = len(atoms)

250 C6eff_a = np.empty((na))

251 alpha_a = np.empty((na))

252 R0eff_a = np.empty((na))

253 for a, atom in enumerate(atoms):

254 # free atom values

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

256 # correction for effective C6

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

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

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

260 for a in range(na):

261 for b in range(a, na):

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

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

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

265 C6eff_aa[b, a] = C6eff_aa[a, b]

267 # New implementation by Miguel Caro

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

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

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

272 pbc_c = atoms.get_pbc()

273 EvdW = 0.0

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

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

276 if pbc_c.any():

278 tol = 1.e-5

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

280 # Build list of neighbors

281 n_list = neighbor_list(quantities="ijdDS",

282 a=atoms,

283 cutoff=Reff,

284 self_interaction=False)

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

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

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

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

289 # Look for neighbor pairs

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

291 i = n_list[0][k]

292 j = n_list[1][k]

293 dist = n_list[2][k]

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

295 # repl = n_list[4][k]

296 if j >= i:

297 atom_list[i].append(j)

298 d_list[i].append(dist)

299 v_list[i].append(vect)

300 # r_list[i].append( repl )

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

302 else:

303 atom_list = []

304 d_list = []

305 v_list = []

306 # r_list = []

307 # Do this to avoid double counting

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

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

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

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

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

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

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

315 # No PBC means we are in the same cell

317 # Here goes the calculation, valid with and without

318 # PBC because we loop over

319 # independent pairwise *interactions*

320 ms = myslice(len(atoms), self.comm)

321 for i in range(len(atoms))[ms]:

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

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

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

325 r6 = r**6

326 Edamp, Fdamp = self.damping(r,

327 R0eff_a[i],

328 R0eff_a[j],

329 d=self.d,

330 sR=self.sR)

331 if pbc_c.any():

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

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

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

335 else:

336 smooth = 1.

337 smooth_der = 0.

338 # Here we compute the contribution to the energy

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

340 # We correct it here

341 if i == j:

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

343 else:

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

345 # Here we compute the contribution to the forces

346 # We neglect the C6eff contribution to the forces

347 # (which can actually be larger

348 # than the other contributions)

349 # Self interactions do not contribute to the forces

350 if i != j:

351 # Force on i due to j

352 force_ij = -(

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

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

355 # Forces go both ways for every interaction

356 forces[i] += force_ij

357 forces[j] -= force_ij

358 EvdW = self.comm.sum(EvdW)

359 self.comm.sum(forces)

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

362 self.results['free_energy'] += EvdW

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

365 if self.txt:

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

367 print(f'vdW correction: {EvdW}', file=self.txt)

368 print(f'Energy: {self.results["energy"]}',

369 file=self.txt)

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

371 symbols = self.atoms.get_chemical_symbols()

372 for ia, symbol in enumerate(symbols):

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

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

375 file=self.txt)

376 self.txt.flush()

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

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

380 sR=0.94):

381 """Damping factor.

383 Standard values for d and sR as given in

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

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

386 x = RAB * scale

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

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

391def calculate_ts09_polarizability(atoms):

392 """Calculate polarizability tensor

394 atoms: Atoms object

395 The atoms object must have a vdWTkatchenko90prl calculator attached.

397 Returns

398 -------

399 polarizability tensor:

400 Unit (e^2 Angstrom^2 / eV).

401 Multiply with Bohr * Ha to get (Angstrom^3)

402 """

403 calc = atoms.calc

404 assert isinstance(calc, vdWTkatchenko09prl)

405 atoms.get_potential_energy()

407 volume_ratios = calc.hirshfeld.get_effective_volume_ratios()

409 na = len(atoms)

410 alpha_a = np.empty((na))

411 alpha_eff_a = np.empty((na))

412 for a, atom in enumerate(atoms):

413 # free atom values

414 alpha_a[a], _ = calc.vdWDB_alphaC6[atom.symbol]

415 # effective polarizability assuming linear combination

416 # of atomic polarizability from ts09

417 alpha_eff_a[a] = volume_ratios[a] * alpha_a[a]

419 alpha = np.sum(alpha_eff_a) * Bohr**2 / Hartree

420 return np.diag([alpha] * 3)

423class TS09Polarizability(StaticPolarizabilityCalculator):

424 """Class interface as expected by Displacement"""

426 def __call__(self, atoms):

427 return calculate_ts09_polarizability(atoms)