1from typing import Tuple
2import numpy as np
3
4from ase.units import Bohr, Ha
5from ase.data import covalent_radii
6from ase.neighborlist import NeighborList
7
8
9class LippincottStuttman:
10    # atomic polarizability values from:
11    #   Lippincott and Stutman J. Phys. Chem. 68 (1964) 2926-2940
12    #   DOI: 10.1021/j100792a033
13    # see also:
14    #   Marinov and Zotov Phys. Rev. B 55 (1997) 2938-2944
15    #   DOI: 10.1103/PhysRevB.55.2938
16    # unit: Angstrom^3
17    atomic_polarizability = {
18        'B': 1.358,
19        'C': 0.978,
20        'N': 0.743,
21        'O': 0.592,
22        'Al': 3.918,
23        'Si': 2.988,
24    }
25    # reduced electronegativity Table I
26    reduced_eletronegativity = {
27        'B': 0.538,
28        'C': 0.846,
29        'N': 0.927,
30        'O': 1.0,
31        'Al': 0.533,
32        'Si': 0.583,
33    }
34
35    def __call__(self, el1: str, el2: str,
36                 length: float) -> Tuple[float, float]:
37        """Bond polarizability
38
39        Parameters
40        ----------
41        el1: element string
42        el2: element string
43        length: float
44
45        Returns
46        -------
47        alphal: float
48          Parallel component
49        alphap: float
50          Perpendicular component
51        """
52        alpha1 = self.atomic_polarizability[el1]
53        alpha2 = self.atomic_polarizability[el2]
54        ren1 = self.reduced_eletronegativity[el1]
55        ren2 = self.reduced_eletronegativity[el2]
56
57        sigma = 1.
58        if el1 != el2:
59            sigma = np.exp(- (ren1 - ren2)**2 / 4)
60
61        # parallel component
62        alphal = sigma * length**4 / (4**4 * alpha1 * alpha2)**(1. / 6)
63        # XXX consider fractional covalency ?
64
65        # prependicular component
66        alphap = ((ren1**2 * alpha1 + ren2**2 * alpha2)
67                  / (ren1**2 + ren2**2))
68        # XXX consider fractional covalency ?
69
70        return alphal, alphap
71
72
73class Linearized:
74    def __init__(self):
75        self._data = {
76            # L. Wirtz, M. Lazzeri, F. Mauri, A. Rubio,
77            # Phys. Rev. B 2005, 71, 241402.
78            #      R0     al    al'   ap    ap'
79            'CC': (1.53, 1.69, 7.43, 0.71, 0.37),
80            'BN': (1.56, 1.58, 4.22, 0.42, 0.90),
81        }
82
83    def __call__(self, el1: str, el2: str,
84                 length: float) -> Tuple[float, float]:
85        """Bond polarizability
86
87        Parameters
88        ----------
89        el1: element string
90        el2: element string
91        length: float
92
93        Returns
94        -------
95        alphal: float
96          Parallel component
97        alphap: float
98          Perpendicular component
99        """
100        if el1 > el2:
101            bond = el2 + el1
102        else:
103            bond = el1 + el2
104        assert bond in self._data
105        length0, al, ald, ap, apd = self._data[bond]
106
107        return al + ald * (length - length0), ap + apd * (length - length0)
108
109
110class BondPolarizability:
111    def __init__(self, model=LippincottStuttman()):
112        self.model = model
113
114    def __call__(self, *args, **kwargs):
115        """Shorthand for calculate"""
116        return self.calculate(*args, **kwargs)
117
118    def calculate(self, atoms, radiicut=1.5):
119        """Sum up the bond polarizability from all bonds
120
121        Parameters
122        ----------
123        atoms: Atoms object
124        radiicut: float
125          Bonds are counted up to
126          radiicut * (sum of covalent radii of the pairs)
127          Default: 1.5
128
129        Returns
130        -------
131        polarizability tensor with unit (e^2 Angstrom^2 / eV).
132        Multiply with Bohr * Ha to get (Angstrom^3)
133        """
134        radii = np.array([covalent_radii[z]
135                          for z in atoms.numbers])
136        nl = NeighborList(radii * 1.5, skin=0,
137                          self_interaction=False)
138        nl.update(atoms)
139        pos_ac = atoms.get_positions()
140
141        alpha = 0
142        for ia, atom in enumerate(atoms):
143            indices, offsets = nl.get_neighbors(ia)
144            pos_ac = atoms.get_positions() - atoms.get_positions()[ia]
145
146            for ib, offset in zip(indices, offsets):
147                weight = 1
148                if offset.any():  # this comes from a periodic image
149                    weight = 0.5  # count half the bond only
150
151                dist_c = pos_ac[ib] + np.dot(offset, atoms.get_cell())
152                dist = np.linalg.norm(dist_c)
153                al, ap = self.model(atom.symbol, atoms[ib].symbol, dist)
154
155                eye3 = np.eye(3) / 3
156                alpha += weight * (al + 2 * ap) * eye3
157                alpha += weight * (al - ap) * (
158                    np.outer(dist_c, dist_c) / dist**2 - eye3)
159        return alpha / Bohr / Ha
160