Coverage for /builds/ase/ase/ase/calculators/vdwcorrection.py : 87.29%
Hot-keys on this page
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,
146 hirshfeld=None, vdwradii=None, calculator=None,
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)
172 self.vdwradii = vdwradii
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:
231 # external vdW radii
232 vdwradii = self.vdwradii
233 assert len(atoms) == len(vdwradii)
234 else:
235 vdwradii = []
236 for atom in atoms:
237 self.vdwradii.append(vdWDB_Grimme06jcc[atom.symbol][1])
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
271 # use the cutoff radius instead
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():
277 # Effective cutoff radius
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)