Coverage for gpaw/gauss.py: 79%

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 R.any(): 

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)))