1# Copyright (C) 2003  CAMP
2# Please see the accompanying LICENSE file for further information.
3
4import numpy as np
5from scipy.special import erf
6
7
8def I(R, a, b, alpha, beta):  # noqa
9    """Calculate integral and derivatives wrt. positions of Gaussian product.
10
11    ::
12
13                /                              2             2
14               |  a0+b0  a1+b1  a2+b2  -alpha r   -beta (r-R)
15       value = | x      y      z      e          e             dr,
16               |
17              /
18
19    Returns the tuple (value, d[value]/dx, d[value]/dy, d[value]/dz).
20    """
21    result = np.zeros(4)
22    R = np.array(R)
23    result[0] = I1(R, a, b, alpha, beta)
24    a = np.array(a)
25    for i in range(3):
26        a[i] += 1
27        result[1 + i] = 2 * alpha * I1(R, tuple(a), b, alpha, beta)
28        a[i] -= 2
29        if a[i] >= 0:
30            result[1 + i] -= (a[i] + 1) * I1(R, tuple(a), b, alpha, beta)
31        a[i] += 1
32    return result
33
34
35def I1(R, ap1, b, alpha, beta, m=0):
36    if ap1 == (0, 0, 0):
37        if b != (0, 0, 0):
38            return I1(-R, b, ap1, beta, alpha, m)
39        else:
40            f = 2 * np.sqrt(np.pi**5 / (alpha + beta)) / (alpha * beta)
41            if np.sometrue(R):
42                T = alpha * beta / (alpha + beta) * np.dot(R, R)
43                f1 = f * erf(T**0.5) * (np.pi / T)**0.5
44                if m == 0:
45                    return 0.5 * f1
46                f2 = f * np.exp(-T) / T**m
47                if m == 1:
48                    return 0.25 * f1 / T - 0.5 * f2
49                if m == 2:
50                    return 0.375 * f1 / T**2 - 0.5 * f2 * (1.5 + T)
51                if m == 3:
52                    return 0.9375 * f1 / T**3 - 0.25 * f2 * (7.5 +
53                                                             T * (5 + 2 * T))
54                if m == 4:
55                    return 3.28125 * f1 / T**4 - 0.125 * f2 * \
56                        (52.5 + T * (35 + 2 * T * (7 + 2 * T)))
57                if m == 5:
58                    return 14.7656 * f1 / T**5 - 0.03125 * f2 * \
59                        (945 + T * (630 + T * (252 + T * (72 + T * 16))))
60                if m == 6:
61                    return 81.2109 * f1 / T**6 - 0.015625 * f2 * \
62                        (10395 + T *
63                         (6930 + T *
64                          (2772 + T * (792 + T * (176 + T * 32)))))
65                else:
66                    raise NotImplementedError
67
68            return f / (1 + 2 * m)
69    for i in range(3):
70        if ap1[i] > 0:
71            break
72    a = ap1[:i] + (ap1[i] - 1,) + ap1[i + 1:]
73    result = beta / (alpha + beta) * R[i] * I1(R, a, b, alpha, beta, m + 1)
74    if a[i] > 0:
75        am1 = a[:i] + (a[i] - 1,) + a[i + 1:]
76        result += a[i] / (2 * alpha) * (I1(R, am1, b, alpha, beta, m) -
77                                        beta / (alpha + beta) *
78                                        I1(R, am1, b, alpha, beta, m + 1))
79    if b[i] > 0:
80        bm1 = b[:i] + (b[i] - 1,) + b[i + 1:]
81        result += b[i] / (2 * (alpha + beta)) * I1(R,
82                                                   a, bm1, alpha, beta, m + 1)
83    return result
84
85
86def test_derivatives(R, a, b, alpha, beta, i):
87    R = np.array(R)
88    a = np.array(a)
89    a[i] += 1
90    dIdRi = 2 * alpha * I1(R, tuple(a), b, alpha, beta)
91    a[i] -= 2
92    if a[i] >= 0:
93        dIdRi -= (a[i] + 1) * I1(R, tuple(a), b, alpha, beta)
94    a[i] += 1
95    dr = 0.001
96    R[i] += 0.5 * dr
97    dIdRi2 = I1(R, tuple(a), b, alpha, beta)
98    R[i] -= dr
99    dIdRi2 -= I1(R, tuple(a), b, alpha, beta)
100    dIdRi2 /= -dr
101    R[i] += 0.5 * dr
102    return dIdRi, dIdRi2
103
104
105class Gauss:
106    """Normalised Gauss distribution
107
108    from gpaw.gauss import Gauss
109
110    width=0.4
111    gs = Gauss(width)
112
113    for i in range(4):
114        print('Gauss(i)=', gs.get(i))
115    """
116    def __init__(self, width=0.08):
117        self.dtype = float
118        self.set_width(width)
119
120    def get(self, x, x0=0):
121        return self.norm * np.exp(-((x - x0) * self.wm1)**2)
122
123    def set_width(self, width=0.08):
124        self.norm = 1. / width / np.sqrt(2 * np.pi)
125        self.wm1 = np.sqrt(.5) / width
126        self._fwhm = np.sqrt(8 * np.log(2)) * width
127
128    @property
129    def fwhm(self):
130        return self._fwhm
131
132    @fwhm.setter
133    def fwhm(self, value):
134        self.set_width(value / np.sqrt(8 * np.log(2)))
135