Coverage for gpaw/test/test_gauss_func.py: 100%

1from math import pi, sqrt 

2import numpy as np 

3from gpaw.utilities.tools import coordinates 

4from gpaw.utilities.gauss import Gaussian 

5from gpaw.grid_descriptor import GridDescriptor 

6import pytest 

7from gpaw.poisson import PoissonSolver 

8 

9 

10def test_gauss_func(): 

11 def norm(a): 

12 return np.sqrt(np.sum(a.ravel()**2)) / len(a.ravel()) 

13 

14 # Initialize classes 

15 a = 20 # Size of cell 

16 N = 48 # Number of grid points 

17 Nc = (N, N, N) # Number of grid points along each axis 

18 gd = GridDescriptor(Nc, (a, a, a), 0) # Grid-descriptor object 

19 solver = PoissonSolver(nn=3) # Numerical poisson solver 

20 solver.set_grid_descriptor(gd) 

21 solve = solver.solve 

22 xyz, r2 = coordinates(gd) # Matrix with square of the radial coordinate 

23 print(r2.shape) 

24 r = np.sqrt(r2) # Matrix with the values of the radial coordinate 

25 nH = np.exp(-2 * r) / pi # Density of the hydrogen atom 

26 gauss = Gaussian(gd) # An instance of Gaussian 

27 

28 # /------------------------------------------------\ 

29 # | Check if Gaussian densities are made correctly | 

30 # \------------------------------------------------/ 

31 for gL in range(2, 9): 

32 g = gauss.get_gauss(gL) # a gaussian of gL'th order 

33 print('\nGaussian of order', gL) 

34 for mL in range(9): 

35 m = gauss.get_moment(g, mL) # the mL'th moment of g 

36 print(f' {mL}\'th moment = {m:2.6f}') 

37 assert m == pytest.approx(float(gL == mL), abs=1e-4) 

38 

39 # Check the moments of the constructed 1s density 

40 print('\nDensity of Hydrogen atom') 

41 for L in range(4): 

42 m = gauss.get_moment(nH, L) 

43 print(f' {L}\'th moment = {m:2.6f}') 

44 assert m == pytest.approx((L == 0) / sqrt(4 * pi), abs=1.5e-3) 

45 

46 # Check that it is removed correctly 

47 gauss.remove_moment(nH, 0) 

48 m = gauss.get_moment(nH, 0) 

49 print('\nZero\'th moment of compensated Hydrogen density =', m) 

50 assert m == pytest.approx(0., abs=1e-7) 

51 

52 # /-------------------------------------------------\ 

53 # | Check if Gaussian potentials are made correctly | 

54 # \-------------------------------------------------/ 

55 

56 # Array for storing the potential 

57 pot = gd.zeros(dtype=float, global_array=False) 

58 for L in range(7): # Angular index of gaussian 

59 # Get analytic functions 

60 ng = gauss.get_gauss(L) 

61 vg = gauss.get_gauss_pot(L) 

62 

63 # Solve potential numerically 

64 solve(pot, ng, charge=None, zero_initial_phi=True) 

65 

66 # Determine residual 

67 residual = norm(pot - vg) 

68 residual = gd.integrate((pot - vg)**2)**0.5 

69 

70 # print result 

71 print('L={}, processor {} of {}: {}'.format( 

72 L, 

73 gd.comm.rank + 1, 

74 gd.comm.size, 

75 residual)) 

76 

77 assert residual < 0.6 

78 

79 # mpirun -np 2 python gauss_func.py --gpaw-parallel --gpaw-debug