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
6from gpaw.test import equal
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('  %s\'th moment = %2.6f' % (mL, m))
37            equal(m, gL == mL, 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('  %s\'th moment = %2.6f' % (L, m))
44        equal(m, (L == 0) / sqrt(4 * pi), 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    equal(m, 0., 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=%s, processor %s of %s: %s' % (
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
80