Coverage for gpaw/test/test_coulomb.py: 91%

1import numpy as np 

2from math import pi 

3from gpaw.coulomb import Coulomb 

4from gpaw.grid_descriptor import GridDescriptor 

5from gpaw.mpi import world 

6from gpaw.utilities.gauss import coordinates 

7import pytest 

8import time 

9 

10 

11def test_coulomb(): 

12 def test_coulomb(N=2**6, a=20): 

13 Nc = (N, N, N) # Number of grid point 

14 gd = GridDescriptor(Nc, (a, a, a), True) # grid-descriptor object 

15 # matrix with the square of the radial coordinate 

16 xyz, r2 = coordinates(gd) 

17 # matrix with the values of the radial coordinate 

18 r = np.sqrt(r2) 

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

20 C = Coulomb(gd) # coulomb calculator 

21 

22 if world.size > 1: 

23 C.load('real') 

24 t0 = time.time() 

25 print('Processor {} of {}: {} Ha in {} sec'.format( 

26 gd.comm.rank + 1, 

27 gd.comm.size, 

28 -0.5 * C.coulomb(nH, method='real'), 

29 time.time() - t0)) 

30 return 

31 else: 

32 C.load('recip_ewald') 

33 C.load('recip_gauss') 

34 C.load('real') 

35 test = {} 

36 t0 = time.time() 

37 test['dual density'] = (-0.5 * C.coulomb(nH, nH.copy()), 

38 time.time() - t0) 

39 for method in ('real', 'recip_gauss', 'recip_ewald'): 

40 t0 = time.time() 

41 test[method] = (-0.5 * C.coulomb(nH, method=method), 

42 time.time() - t0) 

43 return test 

44 

45 analytic = -5 / 16.0 

46 res = test_coulomb(N=48, a=15) 

47 if world.size == 1: 

48 print('Units: Bohr and Hartree') 

49 print('%12s %8s %8s' % ('Method', 'Energy', 'Time')) 

50 print('%12s %2.6f %6s' % ('analytic', analytic, '--')) 

51 for method, et in res.items(): 

52 print('%12s %2.6f %1.7f' % ((method,) + et)) 

53 

54 assert res['real'][0] == pytest.approx(analytic, abs=6e-3) 

55 assert res['recip_gauss'][0] == pytest.approx(analytic, abs=6e-3) 

56 assert res['recip_ewald'][0] == pytest.approx(analytic, abs=2e-2) 

57 assert res['dual density'][0] == pytest.approx(res['recip_gauss'][0], 

58 abs=1e-9) 

59 

60 # mpirun -np 2 python coulomb.py --gpaw-parallel --gpaw-debug