vdwcorrection.py 14.6 KB
Newer Older
1
2
3
"""van der Waals correction schemes for DFT"""
import numpy as np
from ase.units import Bohr, Hartree
miwalter's avatar
miwalter committed
4
from ase.calculators.calculator import Calculator
5
from ase.calculators.polarizability import StaticPolarizabilityCalculator
6
7
from scipy.special import erfinv, erfc
from ase.neighborlist import neighbor_list
Michael Walter's avatar
Michael Walter committed
8
from ase.parallel import world, myslice
9
from ase.utils import IOContext
10

11

Jens Jørgen Mortensen's avatar
Jens Jørgen Mortensen committed
12
# dipole polarizabilities and C6 values from
miwalter's avatar
miwalter committed
13
# X. Chu and A. Dalgarno, J. Chem. Phys. 121 (2004) 4083
14
15
16
# atomic units, a_0^3
vdWDB_Chu04jcp = {
    # Element: [alpha, C6]; units [Bohr^3, Hartree * Bohr^6]
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
    'H': [4.5, 6.5],  # [exact, Tkatchenko PRL]
    'He': [1.38, 1.42],
    'Li': [164, 1392],
    'Be': [38, 227],
    'B': [21, 99.5],
    'C': [12, 46.6],
    'N': [7.4, 24.2],
    'O': [5.4, 15.6],
    'F': [3.8, 9.52],
    'Ne': [2.67, 6.20],
    'Na': [163, 1518],
    'Mg': [71, 626],
    'Al': [60, 528],
    'Si': [37, 305],
    'P': [25, 185],
    'S': [19.6, 134],
    'Cl': [15, 94.6],
    'Ar': [11.1, 64.2],
    'Ca': [160, 2163],
    'Sc': [120, 1383],
    'Ti': [98, 1044],
    'V': [84, 832],
    'Cr': [78, 602],
    'Mn': [63, 552],
    'Fe': [56, 482],
    'Co': [50, 408],
    'Ni': [48, 373],
    'Cu': [42, 253],
    'Zn': [40, 284],
    'As': [29, 246],
    'Se': [25, 210],
    'Br': [20, 162],
    'Kr': [16.7, 130],
    'Sr': [199, 3175],
    'Te': [40, 445],
    'I': [35, 385]}

miwalter's avatar
miwalter committed
54
vdWDB_alphaC6 = vdWDB_Chu04jcp
55
56
57
58
59

# dipole polarizabilities and C6 values from
# V. G. Ruiz et al. Phys. Rev. Lett 108 (2012) 146103
# atomic units, a_0^3
vdWDB_Ruiz12prl = {
Jens Jørgen Mortensen's avatar
Jens Jørgen Mortensen committed
60
61
62
63
    'Ag': [50.6, 339],
    'Au': [36.5, 298],
    'Pd': [23.7, 158],
    'Pt': [39.7, 347],
64
65
66
}

vdWDB_alphaC6.update(vdWDB_Ruiz12prl)
67

Jens Jørgen Mortensen's avatar
Jens Jørgen Mortensen committed
68
# C6 values and vdW radii from
69
# S. Grimme, J Comput Chem 27 (2006) 1787-1799
70
71
vdWDB_Grimme06jcc = {
    # Element: [C6, R0]; units [J nm^6 mol^{-1}, Angstrom]
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
    'H': [0.14, 1.001],
    'He': [0.08, 1.012],
    'Li': [1.61, 0.825],
    'Be': [1.61, 1.408],
    'B': [3.13, 1.485],
    'C': [1.75, 1.452],
    'N': [1.23, 1.397],
    'O': [0.70, 1.342],
    'F': [0.75, 1.287],
    'Ne': [0.63, 1.243],
    'Na': [5.71, 1.144],
    'Mg': [5.71, 1.364],
    'Al': [10.79, 1.639],
    'Si': [9.23, 1.716],
    'P': [7.84, 1.705],
    'S': [5.57, 1.683],
    'Cl': [5.07, 1.639],
    'Ar': [4.61, 1.595],
    'K': [10.80, 1.485],
    'Ca': [10.80, 1.474],
    'Sc': [10.80, 1.562],
    'Ti': [10.80, 1.562],
    'V': [10.80, 1.562],
    'Cr': [10.80, 1.562],
    'Mn': [10.80, 1.562],
    'Fe': [10.80, 1.562],
    'Co': [10.80, 1.562],
    'Ni': [10.80, 1.562],
    'Cu': [10.80, 1.562],
    'Zn': [10.80, 1.562],
    'Ga': [16.99, 1.650],
    'Ge': [17.10, 1.727],
    'As': [16.37, 1.760],
    'Se': [12.64, 1.771],
    'Br': [12.47, 1.749],
    'Kr': [12.01, 1.727],
    'Rb': [24.67, 1.628],
    'Sr': [24.67, 1.606],
    'Y-Cd': [24.67, 1.639],
    'In': [37.32, 1.672],
    'Sn': [38.71, 1.804],
    'Sb': [38.44, 1.881],
    'Te': [31.74, 1.892],
    'I': [31.50, 1.892],
    'Xe': [29.99, 1.881]}

118

119
120
# Optimal range parameters sR for different XC functionals
# to be used with the Tkatchenko-Scheffler scheme
121
# Reference: M.A. Caro arXiv:1704.00761 (2017)
Jens Jørgen Mortensen's avatar
Jens Jørgen Mortensen committed
122
123
124
125
126
127
128
sR_opt = {'PBE': 0.940,
          'RPBE': 0.590,
          'revPBE': 0.585,
          'PBEsol': 1.055,
          'BLYP': 0.625,
          'AM05': 0.840,
          'PW91': 0.965}
129
130


131
132
133
134
135
136
137
138
139
def get_logging_file_descriptor(calculator):
    if hasattr(calculator, 'log'):
        fd = calculator.log
        if hasattr(fd, 'write'):
            return fd
        if hasattr(fd, 'fd'):
            return fd.fd
    if hasattr(calculator, 'txt'):
        return calculator.txt
140
141


142
class vdWTkatchenko09prl(Calculator, IOContext):
miwalter's avatar
miwalter committed
143
    """vdW correction after Tkatchenko and Scheffler PRL 102 (2009) 073005."""
144

Jens Jørgen Mortensen's avatar
Jens Jørgen Mortensen committed
145
    def __init__(self,
146
                 hirshfeld=None, vdwradii=None, calculator=None,
147
                 Rmax=10.,  # maximal radius for periodic calculations
Jens Jørgen Mortensen's avatar
Jens Jørgen Mortensen committed
148
                 Ldecay=1.,  # decay length for smoothing in periodic calcs
149
                 vdWDB_alphaC6=vdWDB_alphaC6,
Jens Jørgen Mortensen's avatar
Jens Jørgen Mortensen committed
150
                 txt=None, sR=None):
miwalter's avatar
miwalter committed
151
152
153
154
155
156
157
        """Constructor

        Parameters
        ==========
        hirshfeld: the Hirshfeld partitioning object
        calculator: the calculator to get the PBE energy
        """
158
159
160
161
162
        self.hirshfeld = hirshfeld
        if calculator is None:
            self.calculator = self.hirshfeld.get_calculator()
        else:
            self.calculator = calculator
163

164
165
        if txt is None:
            txt = get_logging_file_descriptor(self.calculator)
166
        if hasattr(self.calculator, 'world'):
Michael Walter's avatar
Michael Walter committed
167
            self.comm = self.calculator.world
168
        else:
Michael Walter's avatar
Michael Walter committed
169
170
            self.comm = world  # the best we know
        self.txt = self.openfile(txt, self.comm)
171

172
        self.vdwradii = vdwradii
miwalter's avatar
miwalter committed
173
        self.vdWDB_alphaC6 = vdWDB_alphaC6
174
        self.Rmax = Rmax
175
        self.Ldecay = Ldecay
miwalter's avatar
miwalter committed
176
        self.atoms = None
177

178
179
        if sR is None:
            try:
Miguel Caro's avatar
Miguel Caro committed
180
                xc_name = self.calculator.get_xc_functional()
181
                self.sR = sR_opt[xc_name]
Miguel Caro's avatar
Miguel Caro committed
182
            except KeyError:
183
184
185
                raise ValueError(
                    'Tkatchenko-Scheffler dispersion correction not ' +
                    'implemented for %s functional' % xc_name)
186
187
        else:
            self.sR = sR
miwalter's avatar
miwalter committed
188
189
        self.d = 20

miwalter's avatar
miwalter committed
190
191
        Calculator.__init__(self)

192
        self.parameters['calculator'] = self.calculator.name
Michael Walter's avatar
Michael Walter committed
193
        self.parameters['xc'] = self.calculator.get_xc_functional()
194

195
196
197
198
    @property
    def implemented_properties(self):
        return self.calculator.implemented_properties

199
    def calculation_required(self, atoms, quantities):
200
        if self.calculator.calculation_required(atoms, quantities):
201
            return True
miwalter's avatar
miwalter committed
202
203
204
205
206
        for quantity in quantities:
            if quantity not in self.results:
                return True
        return False

207
    def calculate(self, atoms=None, properties=['energy', 'forces'],
miwalter's avatar
miwalter committed
208
209
210
                  system_changes=[]):
        Calculator.calculate(self, atoms, properties, system_changes)
        self.update(atoms, properties)
211

Michael Walter's avatar
Michael Walter committed
212
213
    def update(self, atoms=None,
               properties=['energy', 'free_energy', 'forces']):
miwalter's avatar
miwalter committed
214
        if not self.calculation_required(atoms, properties):
215
216
            return

217
218
        if atoms is None:
            atoms = self.calculator.get_atoms()
219
220

        properties = list(properties)
Michael Walter's avatar
Michael Walter committed
221
        for name in 'energy', 'free_energy', 'forces':
222
223
224
225
226
            if name not in properties:
                properties.append(name)

        for name in properties:
            self.results[name] = self.calculator.get_property(name, atoms)
227
        self.parameters['uncorrected_energy'] = self.results['energy']
228
        self.atoms = atoms.copy()
229

230
231
232
        if self.vdwradii is not None:
            # external vdW radii
            vdwradii = self.vdwradii
Ask Hjorth Larsen's avatar
Ask Hjorth Larsen committed
233
            assert len(atoms) == len(vdwradii)
234
235
236
        else:
            vdwradii = []
            for atom in atoms:
237
                self.vdwradii.append(vdWDB_Grimme06jcc[atom.symbol][1])
238

239
        if self.hirshfeld is None:
240
            volume_ratios = [1.] * len(atoms)
241
        elif hasattr(self.hirshfeld, '__len__'):  # a list
Ask Hjorth Larsen's avatar
Ask Hjorth Larsen committed
242
            assert len(atoms) == len(self.hirshfeld)
243
            volume_ratios = self.hirshfeld
244
        else:  # should be an object
miwalter's avatar
miwalter committed
245
            self.hirshfeld.initialize()
246
247
            volume_ratios = self.hirshfeld.get_effective_volume_ratios()

248
249
250
251
252
253
254
        # correction for effective C6
        na = len(atoms)
        C6eff_a = np.empty((na))
        alpha_a = np.empty((na))
        R0eff_a = np.empty((na))
        for a, atom in enumerate(atoms):
            # free atom values
miwalter's avatar
miwalter committed
255
            alpha_a[a], C6eff_a[a] = self.vdWDB_alphaC6[atom.symbol]
256
257
            # correction for effective C6
            C6eff_a[a] *= Hartree * volume_ratios[a]**2 * Bohr**6
258
            R0eff_a[a] = vdwradii[a] * volume_ratios[a]**(1 / 3.)
259
260
261
262
263
        C6eff_aa = np.empty((na, na))
        for a in range(na):
            for b in range(a, na):
                C6eff_aa[a, b] = (2 * C6eff_a[a] * C6eff_a[b] /
                                  (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]
266

Jens Jørgen Mortensen's avatar
Jens Jørgen Mortensen committed
267
268
        # New implementation by Miguel Caro
        # (complaints etc to mcaroba@gmail.com)
Miguel Caro's avatar
Miguel Caro committed
269
270
271
272
273
274
275
276
277
278
        # If all 3 PBC are False, we do the summation over the atom
        # pairs in the simulation box. If any of them is True, we
        # use the cutoff radius instead
        pbc_c = atoms.get_pbc()
        EvdW = 0.0
        forces = 0. * self.results['forces']
        # PBC: we build a neighbor list according to the Reff criterion
        if pbc_c.any():
            # Effective cutoff radius
            tol = 1.e-5
Jens Jørgen Mortensen's avatar
Jens Jørgen Mortensen committed
279
            Reff = self.Rmax + self.Ldecay * erfinv(1. - 2. * tol)
Miguel Caro's avatar
Miguel Caro committed
280
            # Build list of neighbors
Jens Jørgen Mortensen's avatar
Jens Jørgen Mortensen committed
281
282
283
            n_list = neighbor_list(quantities="ijdDS",
                                   a=atoms,
                                   cutoff=Reff,
Miguel Caro's avatar
Miguel Caro committed
284
285
286
287
                                   self_interaction=False)
            atom_list = [[] for _ in range(0, len(atoms))]
            d_list = [[] for _ in range(0, len(atoms))]
            v_list = [[] for _ in range(0, len(atoms))]
Jens Jørgen Mortensen's avatar
Jens Jørgen Mortensen committed
288
            # r_list = [[] for _ in range(0, len(atoms))]
Miguel Caro's avatar
Miguel Caro committed
289
290
291
292
293
            # Look for neighbor pairs
            for k in range(0, len(n_list[0])):
                i = n_list[0][k]
                j = n_list[1][k]
                dist = n_list[2][k]
Jens Jørgen Mortensen's avatar
Jens Jørgen Mortensen committed
294
295
                vect = n_list[3][k]  # vect is the distance rj - ri
                # repl = n_list[4][k]
Miguel Caro's avatar
Miguel Caro committed
296
                if j >= i:
Jens Jørgen Mortensen's avatar
Jens Jørgen Mortensen committed
297
298
299
300
                    atom_list[i].append(j)
                    d_list[i].append(dist)
                    v_list[i].append(vect)
                    # r_list[i].append( repl )
Miguel Caro's avatar
Miguel Caro committed
301
302
303
304
305
        # Not PBC: we loop over all atoms in the unit cell only
        else:
            atom_list = []
            d_list = []
            v_list = []
Jens Jørgen Mortensen's avatar
Jens Jørgen Mortensen committed
306
            # r_list = []
Miguel Caro's avatar
Miguel Caro committed
307
308
            # Do this to avoid double counting
            for i in range(0, len(atoms)):
Jens Jørgen Mortensen's avatar
Jens Jørgen Mortensen committed
309
310
311
312
313
314
315
                atom_list.append(range(i + 1, len(atoms)))
                d_list.append([atoms.get_distance(i, j)
                               for j in range(i + 1, len(atoms))])
                v_list.append([atoms.get_distance(i, j, vector=True)
                               for j in range(i + 1, len(atoms))])
                # r_list.append( [[0,0,0] for j in range(i+1, len(atoms))])
                # No PBC means we are in the same cell
316

Jens Jørgen Mortensen's avatar
Jens Jørgen Mortensen committed
317
318
        # Here goes the calculation, valid with and without
        # PBC because we loop over
Miguel Caro's avatar
Miguel Caro committed
319
        # independent pairwise *interactions*
Michael Walter's avatar
Michael Walter committed
320
321
        ms = myslice(len(atoms), self.comm)
        for i in range(len(atoms))[ms]:
Jens Jørgen Mortensen's avatar
Jens Jørgen Mortensen committed
322
323
            # for j, r, vect, repl in zip(atom_list[i], d_list[i],
            #                             v_list[i], r_list[i]):
Miguel Caro's avatar
Miguel Caro committed
324
325
326
327
328
329
330
331
332
            for j, r, vect in zip(atom_list[i], d_list[i], v_list[i]):
                r6 = r**6
                Edamp, Fdamp = self.damping(r,
                                            R0eff_a[i],
                                            R0eff_a[j],
                                            d=self.d,
                                            sR=self.sR)
                if pbc_c.any():
                    smooth = 0.5 * erfc((r - self.Rmax) / self.Ldecay)
333
                    smooth_der = -1. / np.sqrt(np.pi) / self.Ldecay * np.exp(
Jens Jørgen Mortensen's avatar
Jens Jørgen Mortensen committed
334
                        -((r - self.Rmax) / self.Ldecay)**2)
Miguel Caro's avatar
Miguel Caro committed
335
336
337
338
                else:
                    smooth = 1.
                    smooth_der = 0.
                # Here we compute the contribution to the energy
Jens Jørgen Mortensen's avatar
Jens Jørgen Mortensen committed
339
340
                # Self interactions (only possible in PBC) are double counted.
                # We correct it here
Miguel Caro's avatar
Miguel Caro committed
341
342
343
344
345
                if i == j:
                    EvdW -= (Edamp * C6eff_aa[i, j] / r6) / 2. * smooth
                else:
                    EvdW -= (Edamp * C6eff_aa[i, j] / r6) * smooth
                # Here we compute the contribution to the forces
Jens Jørgen Mortensen's avatar
Jens Jørgen Mortensen committed
346
347
                # We neglect the C6eff contribution to the forces
                # (which can actually be larger
Miguel Caro's avatar
Miguel Caro committed
348
349
350
351
352
                # than the other contributions)
                # Self interactions do not contribute to the forces
                if i != j:
                    # Force on i due to j
                    force_ij = -(
Jens Jørgen Mortensen's avatar
Jens Jørgen Mortensen committed
353
354
                        (Fdamp - 6 * Edamp / r) * C6eff_aa[i, j] / r6 * smooth
                        + (Edamp * C6eff_aa[i, j] / r6) * smooth_der) * vect / r
Miguel Caro's avatar
Miguel Caro committed
355
356
357
                    # Forces go both ways for every interaction
                    forces[i] += force_ij
                    forces[j] -= force_ij
Michael Walter's avatar
Michael Walter committed
358
359
        EvdW = self.comm.sum(EvdW)
        self.comm.sum(forces)
360

Miguel Caro's avatar
Miguel Caro committed
361
        self.results['energy'] += EvdW
Michael Walter's avatar
Michael Walter committed
362
        self.results['free_energy'] += EvdW
Miguel Caro's avatar
Miguel Caro committed
363
        self.results['forces'] += forces
364

miwalter's avatar
miwalter committed
365
        if self.txt:
366
            print(('\n' + self.__class__.__name__), file=self.txt)
Michael Walter's avatar
Michael Walter committed
367
368
            print(f'vdW correction: {EvdW}', file=self.txt)
            print(f'Energy:         {self.results["energy"]}',
369
370
                  file=self.txt)
            print('\nForces in eV/Ang:', file=self.txt)
miwalter's avatar
miwalter committed
371
372
            symbols = self.atoms.get_chemical_symbols()
            for ia, symbol in enumerate(symbols):
373
374
375
                print('%3d %-2s %10.5f %10.5f %10.5f' %
                      ((ia, symbol) + tuple(self.results['forces'][ia])),
                      file=self.txt)
Jens Jørgen Mortensen's avatar
Jens Jørgen Mortensen committed
376
            self.txt.flush()
377

378
    def damping(self, RAB, R0A, R0B,
379
380
                d=20,   # steepness of the step function for PBE
                sR=0.94):
381
382
        """Damping factor.

Jens Jørgen Mortensen's avatar
Jens Jørgen Mortensen committed
383
        Standard values for d and sR as given in
384
        Tkatchenko and Scheffler PRL 102 (2009) 073005."""
miwalter's avatar
miwalter committed
385
386
387
388
        scale = 1.0 / (sR * (R0A + R0B))
        x = RAB * scale
        chi = np.exp(-d * (x - 1.0))
        return 1.0 / (1.0 + chi), d * scale * chi / (1.0 + chi)**2
Michael Walter's avatar
Michael Walter committed
389
390


391
392
def calculate_ts09_polarizability(atoms):
    """Calculate polarizability tensor
393

394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
    atoms: Atoms object
    The atoms object must have a vdWTkatchenko90prl calculator attached.

    Returns
    -------
      polarizability tensor:
      Unit (e^2 Angstrom^2 / eV).
      Multiply with Bohr * Ha to get (Angstrom^3)
    """
    calc = atoms.calc
    assert isinstance(calc, vdWTkatchenko09prl)
    atoms.get_potential_energy()

    volume_ratios = calc.hirshfeld.get_effective_volume_ratios()

    na = len(atoms)
    alpha_a = np.empty((na))
    alpha_eff_a = np.empty((na))
    for a, atom in enumerate(atoms):
        # free atom values
        alpha_a[a], _ = calc.vdWDB_alphaC6[atom.symbol]
        # effective polarizability assuming linear combination
        # of atomic polarizability from ts09
        alpha_eff_a[a] = volume_ratios[a] * alpha_a[a]

    alpha = np.sum(alpha_eff_a) * Bohr**2 / Hartree
    return np.diag([alpha] * 3)


class TS09Polarizability(StaticPolarizabilityCalculator):
    """Class interface as expected by Displacement"""
425

426
    def __call__(self, atoms):
427
        return calculate_ts09_polarizability(atoms)