Coverage for /builds/ase/ase/ase/calculators/vdwcorrection.py: 87.29%
181 statements
« prev ^ index » next coverage.py v7.5.3, created at 2025-08-02 00:12 +0000
1# fmt: off
3"""van der Waals correction schemes for DFT"""
4import numpy as np
5from scipy.special import erfc, erfinv
7from ase.calculators.calculator import Calculator
8from ase.calculators.polarizability import StaticPolarizabilityCalculator
9from ase.neighborlist import neighbor_list
10from ase.parallel import myslice, world
11from ase.units import Bohr, Hartree
12from ase.utils import IOContext
14# dipole polarizabilities and C6 values from
15# X. Chu and A. Dalgarno, J. Chem. Phys. 121 (2004) 4083
16# atomic units, a_0^3
17vdWDB_Chu04jcp = {
18 # Element: [alpha, C6]; units [Bohr^3, Hartree * Bohr^6]
19 'H': [4.5, 6.5], # [exact, Tkatchenko PRL]
20 'He': [1.38, 1.42],
21 'Li': [164, 1392],
22 'Be': [38, 227],
23 'B': [21, 99.5],
24 'C': [12, 46.6],
25 'N': [7.4, 24.2],
26 'O': [5.4, 15.6],
27 'F': [3.8, 9.52],
28 'Ne': [2.67, 6.20],
29 'Na': [163, 1518],
30 'Mg': [71, 626],
31 'Al': [60, 528],
32 'Si': [37, 305],
33 'P': [25, 185],
34 'S': [19.6, 134],
35 'Cl': [15, 94.6],
36 'Ar': [11.1, 64.2],
37 'Ca': [160, 2163],
38 'Sc': [120, 1383],
39 'Ti': [98, 1044],
40 'V': [84, 832],
41 'Cr': [78, 602],
42 'Mn': [63, 552],
43 'Fe': [56, 482],
44 'Co': [50, 408],
45 'Ni': [48, 373],
46 'Cu': [42, 253],
47 'Zn': [40, 284],
48 'As': [29, 246],
49 'Se': [25, 210],
50 'Br': [20, 162],
51 'Kr': [16.7, 130],
52 'Sr': [199, 3175],
53 'Te': [40, 445],
54 'I': [35, 385]}
56vdWDB_alphaC6 = vdWDB_Chu04jcp
58# dipole polarizabilities and C6 values from
59# V. G. Ruiz et al. Phys. Rev. Lett 108 (2012) 146103
60# atomic units, a_0^3
61vdWDB_Ruiz12prl = {
62 'Ag': [50.6, 339],
63 'Au': [36.5, 298],
64 'Pd': [23.7, 158],
65 'Pt': [39.7, 347],
66}
68vdWDB_alphaC6.update(vdWDB_Ruiz12prl)
70# C6 values and vdW radii from
71# S. Grimme, J Comput Chem 27 (2006) 1787-1799
72vdWDB_Grimme06jcc = {
73 # Element: [C6, R0]; units [J nm^6 mol^{-1}, Angstrom]
74 'H': [0.14, 1.001],
75 'He': [0.08, 1.012],
76 'Li': [1.61, 0.825],
77 'Be': [1.61, 1.408],
78 'B': [3.13, 1.485],
79 'C': [1.75, 1.452],
80 'N': [1.23, 1.397],
81 'O': [0.70, 1.342],
82 'F': [0.75, 1.287],
83 'Ne': [0.63, 1.243],
84 'Na': [5.71, 1.144],
85 'Mg': [5.71, 1.364],
86 'Al': [10.79, 1.639],
87 'Si': [9.23, 1.716],
88 'P': [7.84, 1.705],
89 'S': [5.57, 1.683],
90 'Cl': [5.07, 1.639],
91 'Ar': [4.61, 1.595],
92 'K': [10.80, 1.485],
93 'Ca': [10.80, 1.474],
94 'Sc': [10.80, 1.562],
95 'Ti': [10.80, 1.562],
96 'V': [10.80, 1.562],
97 'Cr': [10.80, 1.562],
98 'Mn': [10.80, 1.562],
99 'Fe': [10.80, 1.562],
100 'Co': [10.80, 1.562],
101 'Ni': [10.80, 1.562],
102 'Cu': [10.80, 1.562],
103 'Zn': [10.80, 1.562],
104 'Ga': [16.99, 1.650],
105 'Ge': [17.10, 1.727],
106 'As': [16.37, 1.760],
107 'Se': [12.64, 1.771],
108 'Br': [12.47, 1.749],
109 'Kr': [12.01, 1.727],
110 'Rb': [24.67, 1.628],
111 'Sr': [24.67, 1.606],
112 'Y-Cd': [24.67, 1.639],
113 'In': [37.32, 1.672],
114 'Sn': [38.71, 1.804],
115 'Sb': [38.44, 1.881],
116 'Te': [31.74, 1.892],
117 'I': [31.50, 1.892],
118 'Xe': [29.99, 1.881]}
121# Optimal range parameters sR for different XC functionals
122# to be used with the Tkatchenko-Scheffler scheme
123# Reference: M.A. Caro arXiv:1704.00761 (2017)
124sR_opt = {'PBE': 0.940,
125 'RPBE': 0.590,
126 'revPBE': 0.585,
127 'PBEsol': 1.055,
128 'BLYP': 0.625,
129 'AM05': 0.840,
130 'PW91': 0.965}
133def get_logging_file_descriptor(calculator):
134 if hasattr(calculator, 'log'):
135 fd = calculator.log
136 if hasattr(fd, 'write'):
137 return fd
138 if hasattr(fd, 'fd'):
139 return fd.fd
140 if hasattr(calculator, 'txt'):
141 return calculator.txt
144class vdWTkatchenko09prl(Calculator, IOContext):
145 """vdW correction after Tkatchenko and Scheffler PRL 102 (2009) 073005."""
147 def __init__(self,
148 hirshfeld=None, vdwradii=None, calculator=None,
149 Rmax=10., # maximal radius for periodic calculations
150 Ldecay=1., # decay length for smoothing in periodic calcs
151 vdWDB_alphaC6=vdWDB_alphaC6,
152 txt=None, sR=None, comm=world):
153 """Constructor
155 Parameters
156 ==========
157 hirshfeld: the Hirshfeld partitioning object
158 calculator: the calculator to get the PBE energy
159 """
160 self.hirshfeld = hirshfeld
161 if calculator is None:
162 self.calculator = self.hirshfeld.get_calculator()
163 else:
164 self.calculator = calculator
166 if txt is None:
167 txt = get_logging_file_descriptor(self.calculator)
168 if hasattr(self.calculator, 'world'):
169 self.comm = self.calculator.world
170 else:
171 self.comm = comm # the best we know
172 self.txt = self.openfile(file=txt, comm=self.comm)
174 self.vdwradii = vdwradii
175 self.vdWDB_alphaC6 = vdWDB_alphaC6
176 self.Rmax = Rmax
177 self.Ldecay = Ldecay
178 self.atoms = None
180 if sR is None:
181 try:
182 xc_name = self.calculator.get_xc_functional()
183 self.sR = sR_opt[xc_name]
184 except KeyError:
185 raise ValueError(
186 'Tkatchenko-Scheffler dispersion correction not ' +
187 f'implemented for {xc_name} functional')
188 else:
189 self.sR = sR
190 self.d = 20
192 Calculator.__init__(self)
194 self.parameters['calculator'] = self.calculator.name
195 self.parameters['xc'] = self.calculator.get_xc_functional()
197 @property
198 def implemented_properties(self):
199 return self.calculator.implemented_properties
201 def calculation_required(self, atoms, quantities):
202 if self.calculator.calculation_required(atoms, quantities):
203 return True
204 for quantity in quantities:
205 if quantity not in self.results:
206 return True
207 return False
209 def calculate(self, atoms=None, properties=['energy', 'forces'],
210 system_changes=[]):
211 Calculator.calculate(self, atoms, properties, system_changes)
212 self.update(atoms, properties)
214 def update(self, atoms=None,
215 properties=['energy', 'free_energy', 'forces']):
216 if not self.calculation_required(atoms, properties):
217 return
219 if atoms is None:
220 atoms = self.calculator.get_atoms()
222 properties = list(properties)
223 for name in 'energy', 'free_energy', 'forces':
224 if name not in properties:
225 properties.append(name)
227 for name in properties:
228 self.results[name] = self.calculator.get_property(name, atoms)
229 self.parameters['uncorrected_energy'] = self.results['energy']
230 self.atoms = atoms.copy()
232 if self.vdwradii is not None:
233 # external vdW radii
234 vdwradii = self.vdwradii
235 assert len(atoms) == len(vdwradii)
236 else:
237 vdwradii = []
238 for atom in atoms:
239 self.vdwradii.append(vdWDB_Grimme06jcc[atom.symbol][1])
241 if self.hirshfeld is None:
242 volume_ratios = [1.] * len(atoms)
243 elif hasattr(self.hirshfeld, '__len__'): # a list
244 assert len(atoms) == len(self.hirshfeld)
245 volume_ratios = self.hirshfeld
246 else: # should be an object
247 self.hirshfeld.initialize()
248 volume_ratios = self.hirshfeld.get_effective_volume_ratios()
250 # correction for effective C6
251 na = len(atoms)
252 C6eff_a = np.empty(na)
253 alpha_a = np.empty(na)
254 R0eff_a = np.empty(na)
255 for a, atom in enumerate(atoms):
256 # free atom values
257 alpha_a[a], C6eff_a[a] = self.vdWDB_alphaC6[atom.symbol]
258 # correction for effective C6
259 C6eff_a[a] *= Hartree * volume_ratios[a]**2 * Bohr**6
260 R0eff_a[a] = vdwradii[a] * volume_ratios[a]**(1 / 3.)
261 C6eff_aa = np.empty((na, na))
262 for a in range(na):
263 for b in range(a, na):
264 C6eff_aa[a, b] = (2 * C6eff_a[a] * C6eff_a[b] /
265 (alpha_a[b] / alpha_a[a] * C6eff_a[a] +
266 alpha_a[a] / alpha_a[b] * C6eff_a[b]))
267 C6eff_aa[b, a] = C6eff_aa[a, b]
269 # New implementation by Miguel Caro
270 # (complaints etc to mcaroba@gmail.com)
271 # If all 3 PBC are False, we do the summation over the atom
272 # pairs in the simulation box. If any of them is True, we
273 # use the cutoff radius instead
274 pbc_c = atoms.get_pbc()
275 EvdW = 0.0
276 forces = 0. * self.results['forces']
277 # PBC: we build a neighbor list according to the Reff criterion
278 if pbc_c.any():
279 # Effective cutoff radius
280 tol = 1.e-5
281 Reff = self.Rmax + self.Ldecay * erfinv(1. - 2. * tol)
282 # Build list of neighbors
283 n_list = neighbor_list(quantities="ijdDS",
284 a=atoms,
285 cutoff=Reff,
286 self_interaction=False)
287 atom_list = [[] for _ in range(len(atoms))]
288 d_list = [[] for _ in range(len(atoms))]
289 v_list = [[] for _ in range(len(atoms))]
290 # r_list = [[] for _ in range(0, len(atoms))]
291 # Look for neighbor pairs
292 for k in range(len(n_list[0])):
293 i = n_list[0][k]
294 j = n_list[1][k]
295 dist = n_list[2][k]
296 vect = n_list[3][k] # vect is the distance rj - ri
297 # repl = n_list[4][k]
298 if j >= i:
299 atom_list[i].append(j)
300 d_list[i].append(dist)
301 v_list[i].append(vect)
302 # r_list[i].append( repl )
303 # Not PBC: we loop over all atoms in the unit cell only
304 else:
305 atom_list = []
306 d_list = []
307 v_list = []
308 # r_list = []
309 # Do this to avoid double counting
310 for i in range(len(atoms)):
311 atom_list.append(range(i + 1, len(atoms)))
312 d_list.append([atoms.get_distance(i, j)
313 for j in range(i + 1, len(atoms))])
314 v_list.append([atoms.get_distance(i, j, vector=True)
315 for j in range(i + 1, len(atoms))])
316 # r_list.append( [[0,0,0] for j in range(i+1, len(atoms))])
317 # No PBC means we are in the same cell
319 # Here goes the calculation, valid with and without
320 # PBC because we loop over
321 # independent pairwise *interactions*
322 ms = myslice(len(atoms), self.comm)
323 for i in range(len(atoms))[ms]:
324 # for j, r, vect, repl in zip(atom_list[i], d_list[i],
325 # v_list[i], r_list[i]):
326 for j, r, vect in zip(atom_list[i], d_list[i], v_list[i]):
327 r6 = r**6
328 Edamp, Fdamp = self.damping(r,
329 R0eff_a[i],
330 R0eff_a[j],
331 d=self.d,
332 sR=self.sR)
333 if pbc_c.any():
334 smooth = 0.5 * erfc((r - self.Rmax) / self.Ldecay)
335 smooth_der = -1. / np.sqrt(np.pi) / self.Ldecay * np.exp(
336 -((r - self.Rmax) / self.Ldecay)**2)
337 else:
338 smooth = 1.
339 smooth_der = 0.
340 # Here we compute the contribution to the energy
341 # Self interactions (only possible in PBC) are double counted.
342 # We correct it here
343 if i == j:
344 EvdW -= (Edamp * C6eff_aa[i, j] / r6) / 2. * smooth
345 else:
346 EvdW -= (Edamp * C6eff_aa[i, j] / r6) * smooth
347 # Here we compute the contribution to the forces
348 # We neglect the C6eff contribution to the forces
349 # (which can actually be larger
350 # than the other contributions)
351 # Self interactions do not contribute to the forces
352 if i != j:
353 # Force on i due to j
354 force_ij = -(
355 (Fdamp - 6 * Edamp / r) * C6eff_aa[i, j] / r6 * smooth
356 + (Edamp * C6eff_aa[i, j] / r6) * smooth_der) * vect / r
357 # Forces go both ways for every interaction
358 forces[i] += force_ij
359 forces[j] -= force_ij
360 EvdW = self.comm.sum_scalar(EvdW)
361 self.comm.sum(forces)
363 self.results['energy'] += EvdW
364 self.results['free_energy'] += EvdW
365 self.results['forces'] += forces
367 if self.txt:
368 print(('\n' + self.__class__.__name__), file=self.txt)
369 print(f'vdW correction: {EvdW}', file=self.txt)
370 print(f'Energy: {self.results["energy"]}',
371 file=self.txt)
372 print('\nForces in eV/Ang:', file=self.txt)
373 symbols = self.atoms.get_chemical_symbols()
374 for ia, symbol in enumerate(symbols):
375 print('%3d %-2s %10.5f %10.5f %10.5f' %
376 ((ia, symbol) + tuple(self.results['forces'][ia])),
377 file=self.txt)
378 self.txt.flush()
380 def damping(self, RAB, R0A, R0B,
381 d=20, # steepness of the step function for PBE
382 sR=0.94):
383 """Damping factor.
385 Standard values for d and sR as given in
386 Tkatchenko and Scheffler PRL 102 (2009) 073005."""
387 scale = 1.0 / (sR * (R0A + R0B))
388 x = RAB * scale
389 chi = np.exp(-d * (x - 1.0))
390 return 1.0 / (1.0 + chi), d * scale * chi / (1.0 + chi)**2
393def calculate_ts09_polarizability(atoms):
394 """Calculate polarizability tensor
396 atoms: Atoms object
397 The atoms object must have a vdWTkatchenko90prl calculator attached.
399 Returns
400 -------
401 polarizability tensor:
402 Unit (e^2 Angstrom^2 / eV).
403 Multiply with Bohr * Ha to get (Angstrom^3)
404 """
405 calc = atoms.calc
406 assert isinstance(calc, vdWTkatchenko09prl)
407 atoms.get_potential_energy()
409 volume_ratios = calc.hirshfeld.get_effective_volume_ratios()
411 na = len(atoms)
412 alpha_a = np.empty(na)
413 alpha_eff_a = np.empty(na)
414 for a, atom in enumerate(atoms):
415 # free atom values
416 alpha_a[a], _ = calc.vdWDB_alphaC6[atom.symbol]
417 # effective polarizability assuming linear combination
418 # of atomic polarizability from ts09
419 alpha_eff_a[a] = volume_ratios[a] * alpha_a[a]
421 alpha = np.sum(alpha_eff_a) * Bohr**2 / Hartree
422 return np.diag([alpha] * 3)
425class TS09Polarizability(StaticPolarizabilityCalculator):
426 """Class interface as expected by Displacement"""
428 def __call__(self, atoms):
429 return calculate_ts09_polarizability(atoms)