Coverage for /builds/kinetik161/ase/ase/calculators/vdwcorrection.py: 87.29%
181 statements
« prev ^ index » next coverage.py v7.2.7, created at 2023-12-10 11:04 +0000
1"""van der Waals correction schemes for DFT"""
2import numpy as np
3from scipy.special import erfc, erfinv
5from ase.calculators.calculator import Calculator
6from ase.calculators.polarizability import StaticPolarizabilityCalculator
7from ase.neighborlist import neighbor_list
8from ase.parallel import myslice, world
9from ase.units import Bohr, Hartree
10from 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 f'implemented for {xc_name} functional')
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_scalar(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)