Coverage for gpaw/test/response/test_coulomb.py: 98%

1import pytest 

2from gpaw.grid_descriptor import GridDescriptor 

3from gpaw import PoissonSolver 

4from gpaw.spline import Spline 

5from gpaw.lfc import LocalizedFunctionsCollection as LFC 

6from gpaw.response.qpd import SingleQPWDescriptor 

7import numpy as np 

8from gpaw.response.coulomb_kernels import (get_coulomb_kernel, 

9 get_integrated_kernel) 

10from gpaw.kpt_descriptor import KPointDescriptor 

11from gpaw.poisson_extravacuum import ExtraVacuumPoissonSolver 

12from gpaw.poisson_moment import MomentCorrectionPoissonSolver 

13 

14 

15class ExtraVacuum2DPoisson: 

16 def __init__(self, poisson): 

17 self.poisson = poisson 

18 

19 def set_grid_descriptor(self, gd): 

20 self.N_c = gd.N_c 

21 self.gd = GridDescriptor(gd.N_c * np.array([1, 1, 9]), 

22 gd.cell_cv @ np.diag([1, 1, 9]), 

23 pbc_c=gd.pbc_c) 

24 self.poisson.set_grid_descriptor(self.gd) 

25 

26 def solve(self, v, n): 

27 myv = self.gd.zeros() 

28 myn = self.gd.zeros() 

29 beg, end = self.N_c[2] * 4, self.N_c[2] * 5 - 1 

30 myn[:, :, beg:end] = n 

31 self.poisson.solve(myv, myn) 

32 v[:] = myv[:, :, beg:end] 

33 

34 

35class Grid: 

36 def __init__(self, L, N, q, truncation): 

37 if q == 0: 

38 supercell = 1 

39 else: 

40 supercell = round(1 / q) 

41 assert np.allclose(supercell * q, 1.0) 

42 

43 # First we create a unit cell quantities 

44 pbc_c = [True, True, True] 

45 if truncation is None: 

46 pass 

47 elif truncation == '0D': 

48 pbc_c = [False, False, False] 

49 elif truncation == '2D': 

50 pbc_c = [True, True, False] 

51 else: 

52 raise NotImplementedError 

53 pbc_c = np.array(pbc_c) 

54 self.pbc_c = pbc_c 

55 gd = GridDescriptor((N, N, N), [L, L, L], pbc_c=pbc_c) 

56 self.periodicgd = GridDescriptor((N, N, N), [L, L, L]) 

57 kd = KPointDescriptor([[q, 0, 0]]) 

58 # Create some charge distribution to the unit cell 

59 if truncation == '0D': 

60 cut = L / 9 

61 else: 

62 cut = L / 4 

63 spline = Spline.from_data(l=0, rmax=cut, f_g=np.array([1, 0.5, 0.0])) 

64 c = LFC(gd, [[spline]] * 2, kd=kd, dtype=complex) 

65 

66 if truncation == '0D': 

67 c.set_positions([[0.5 - 0.1, 0.5, 0.5], [0.5 + 0.1, 0.5, 0.5]]) 

68 else: 

69 c.set_positions([[0.11, 0.0, 0.5], [0.5, 0.39, 0.5]]) 

70 n_R = gd.zeros(dtype=complex) 

71 c.add(n_R, {0: np.array([-1], dtype=complex), 

72 1: np.array([1], dtype=complex)}, 0) 

73 

74 # Now we replicate to supercell, to calculate using grid based Poisson 

75 # solvers which cannot handle Bloch-phases at the moment 

76 supergd = GridDescriptor((N * supercell, N, N), 

77 [L * supercell, L, L], 

78 pbc_c=pbc_c) 

79 

80 # Avoid going over 79 characters 

81 EVPS = ExtraVacuumPoissonSolver 

82 MCPS = MomentCorrectionPoissonSolver 

83 PS = PoissonSolver 

84 if truncation is None: 

85 superpoisson = PS(nn=3) 

86 elif truncation == '0D': 

87 mcps = MCPS(PS(nn=3), moment_corrections=9) 

88 superpoisson = EVPS(gpts=(N * 4, N * 4, N * 4), 

89 poissonsolver_large=mcps, 

90 coarses=1, 

91 poissonsolver_small=PS(nn=3)) 

92 elif truncation == '2D': 

93 superpoisson = ExtraVacuum2DPoisson(PS(nn=3)) 

94 else: 

95 raise NotImplementedError(truncation) 

96 superpoisson.set_grid_descriptor(supergd) 

97 

98 if q != 0: 

99 supern_R = supergd.zeros(dtype=complex) 

100 for i in range(supercell): 

101 pn_R = n_R * np.exp(2j * np.pi * i / supercell) 

102 supern_R[(i * N):((i + 1) * N), :, :] = pn_R 

103 

104 else: 

105 supergd = gd 

106 supern_R = n_R 

107 

108 superv_R = supergd.zeros(dtype=complex) 

109 superpoisson.solve(superv_R.real, supern_R.real) 

110 superpoisson.solve(superv_R.imag, supern_R.imag) 

111 

112 self.Ec = supergd.integrate(supern_R, superv_R) / supercell 

113 

114 self.gd = gd 

115 self.n_R = gd.zero_pad(n_R) 

116 

117 

118@pytest.mark.serial 

119@pytest.mark.response 

120@pytest.mark.parametrize('gridparam', [(32, 1e-5), (64, 1e-7), (96, 1e-8)]) 

121@pytest.mark.parametrize('qtrunc', [ 

122 (0, '2D', 20.7454594963, 506.01293778), 

123 (1 / 3, '2D', 13.174804153, 190.916278817), 

124 (0, None, 20.228908696, 578.42826785), 

125 (1 / 3, None, 13.5467930334, 214.823201910), 

126 (0, '0D', 14.3596834829, 206.74182299)]) 

127def test_coulomb(gridparam, qtrunc): 

128 N, maxdev = gridparam 

129 q, truncation, sqrtV0_dev, V0_dev = qtrunc 

130 L = 10 

131 print() 

132 # Slightly lower tolerance for 0D system, because it uses so localized 

133 # charges to avoid crosstalk 

134 if truncation == '0D': 

135 maxdev *= 1e2 

136 grid = Grid(L, N, q, truncation=truncation) 

137 # Use maximum ecut 

138 ecut0 = 0.5 * np.pi**2 / ((L / N)**2) 

139 qpd = SingleQPWDescriptor.from_q([q, 0, 0], ecut0, grid.periodicgd, 

140 gammacentered=False) 

141 # XXX It is silly that get_coulomb_kernel uses k-point numbers 

142 # per dim to determine non-periodic direction. It should just get 

143 # non_per_dir=2, etc. 

144 N_c = np.array([2, 2, 1]) 

145 v_G = get_coulomb_kernel(qpd, N_c, truncation=truncation, pbc_c=grid.pbc_c) 

146 V0, sqrtV0 = get_integrated_kernel( 

147 qpd=qpd, N_c=N_c, truncation=truncation, pbc_c=grid.pbc_c, N=100) 

148 assert V0 == pytest.approx(V0_dev, maxdev) 

149 assert sqrtV0 == pytest.approx(sqrtV0_dev, maxdev) 

150 

151 n_G = qpd.fft(grid.n_R * grid.periodicgd.plane_wave([-q, 0, 0])) 

152 Ec2 = n_G.conj() @ (v_G * n_G) / N**6 * L**3 

153 dev = np.abs(grid.Ec / Ec2 - 1) 

154 print(f'N={N:4d} q=[{q:.3f},0,0] truncation: {str(truncation):5s}' 

155 f'rel. log10 deviation: {np.log10(dev):7.4f}' 

156 f'coulomb:{Ec2:9.5f}') 

157 assert np.abs(dev) < maxdev