Coverage for gpaw/test/pw/test_expert_diag.py: 98%

1import pytest 

2import numpy as np 

3from ase.build import bulk 

4 

5from gpaw import GPAW, PW 

6from gpaw.mpi import world 

7 

8# This test is asserting that the expert diagonalization 

9# routine gives the same result as the non-expert version 

10# in terms of eigenvalues and wavefunctions 

11 

12 

13def test_pw_expert_diag(in_tmp_dir, scalapack): 

14 wfs_e = [] 

15 for i, nbands in enumerate([None, 48]): 

16 si = bulk('Si') 

17 name = f'si_{i}' 

18 si.center() 

19 calc = GPAW(mode=PW(120), kpts=(1, 1, 2), 

20 eigensolver='rmm-diis', 

21 parallel={'domain': 1}, 

22 symmetry='off', txt=name + '.txt') 

23 si.calc = calc 

24 si.get_potential_energy() 

25 calc.diagonalize_full_hamiltonian(nbands=nbands) 

26 string = name + '.gpw' 

27 calc.write(string, 'all') 

28 wfs_e.append(calc.wfs) 

29 

30 # Test against values from reference revision 

31 epsn_n = np.array([-0.21072488, 0.05651796, 

32 0.17629109, 0.17629122, 

33 0.26459912]) 

34 wfsold_G = np.array([7.52087176, 54.64580691, 

35 -3.813492786, -0.14781141, 

36 -0.64662074]) 

37 

38 kpt_u = wfs_e[0].kpt_u 

39 for kpt in kpt_u: 

40 if kpt.k == 0: 

41 psit = kpt.psit_nG[1, 0:5].copy() 

42 if wfsold_G[0] * psit[0] < 0: 

43 psit *= -1. 

44 if world.rank == 0: 

45 assert np.allclose(epsn_n, kpt.eps_n[0:5], atol=1e-4), \ 

46 'Eigenvalues have changed' 

47 assert np.allclose(wfsold_G, psit, atol=5e-3), \ 

48 'Wavefunctions have changed' 

49 

50 wfstmp, wfs = wfs_e 

51 for kpt, kpttmp in zip(wfs.kpt_u, wfstmp.kpt_u): 

52 if wfs.bd.comm.rank == 0: 

53 for m, (psi_G, eps) in enumerate(zip(kpt.psit_nG, kpt.eps_n)): 

54 # Have to do like this if bands are degenerate 

55 booleanarray = np.abs(kpttmp.eps_n - eps) < 1e-10 

56 inds = np.argwhere(booleanarray) 

57 count = len(inds) 

58 assert count > 0, 'Difference between eigenvalues!' 

59 

60 psitmp_nG = kpttmp.psit_nG[inds][:, 0, :] 

61 fidelity = 0 

62 for psitmp_G in psitmp_nG: 

63 fidelity += ( 

64 np.abs(np.dot(psitmp_G.conj(), psi_G))**2 / 

65 np.dot(psitmp_G.conj(), psitmp_G) / 

66 np.dot(psi_G.conj(), psi_G)) 

67 

68 assert fidelity == pytest.approx(1, abs=1e-10) 

69 

70 

71if __name__ == '__main__': 

72 test_pw_expert_diag(1, 1)