Coverage for gpaw/test/lcao/test_density.py: 100%

1"""Test LCAO density calculation and conversion to grid. 

2Test that the density generated by the following three procedures is the same: 

3 

4 * basis_functions.construct_density as used in a normal calculation 

5 * axpy used on the psit_nG as constructed by lcao_to_grid 

6 * axpy used on the Phit_MG[i] * Phit_MG[j] * rho[j, i], where Phit_MG 

7 are the actual basis functions on the grid, constructed using lcao_to_grid 

8 

9TODO: non-gamma-point test 

10 

11""" 

12import numpy as np 

13import pytest 

14from ase.build import molecule 

15 

16from gpaw import GPAW, ConvergenceError 

17from gpaw.utilities.blas import axpy 

18 

19 

20@pytest.mark.old_gpaw_only 

21def test_lcao_density(): 

22 system = molecule('H2O') 

23 system.center(vacuum=2.5) 

24 

25 calc = GPAW(mode='lcao', 

26 maxiter=1) 

27 

28 system.calc = calc 

29 try: 

30 system.get_potential_energy() 

31 except ConvergenceError: 

32 pass 

33 

34 wfs = calc.wfs 

35 kpt = wfs.kpt_u[0] 

36 nt_G = calc.density.gd.zeros() 

37 bfs = wfs.basis_functions 

38 nao = wfs.setups.nao 

39 f_n = kpt.f_n 

40 rho_MM = np.zeros((nao, nao)) 

41 wfs.calculate_density_matrix(kpt.f_n, kpt.C_nM, rho_MM) 

42 

43 bfs.construct_density(rho_MM, nt_G, -1) 

44 

45 nbands = wfs.bd.nbands 

46 psit_nG = wfs.gd.zeros(nbands) 

47 bfs.lcao_to_grid(kpt.C_nM, psit_nG, -1) 

48 

49 nt2_G = calc.density.gd.zeros() 

50 for f, psit_G in zip(f_n, psit_nG): 

51 axpy(f, psit_G**2, nt2_G) 

52 

53 identity_MM = np.identity(nao) 

54 Phit_MG = calc.wfs.gd.zeros(nao) 

55 bfs.lcao_to_grid(identity_MM, Phit_MG, -1) 

56 

57 nt3_G = calc.density.gd.zeros() 

58 for M1, Phit1_G in enumerate(Phit_MG): 

59 for M2, Phit2_G in enumerate(Phit_MG): 

60 nt3_G += rho_MM[M1, M2] * Phit1_G * Phit2_G 

61 

62 err1_G = nt2_G - nt_G 

63 err2_G = nt3_G - nt_G 

64 

65 maxerr1 = np.abs(err1_G).max() 

66 maxerr2 = np.abs(err2_G).max() 

67 

68 print('err1', maxerr1) 

69 print('err2', maxerr2) 

70 

71 assert max(maxerr1, maxerr2) < 1e-15