Coverage for gpaw/test/response/test_silicon_chi.py: 100%

1import time 

2import pytest 

3import numpy as np 

4 

5from ase.build import bulk 

6from ase.parallel import parprint 

7from ase.utils.timing import Timer 

8 

9from gpaw import GPAW, PW, FermiDirac 

10from gpaw.test import findpeak 

11from gpaw.mpi import size, world 

12 

13from gpaw.response import ResponseGroundStateAdapter 

14from gpaw.response.df import DielectricFunction, read_response_function 

15from gpaw.response.chiks import ChiKSCalculator 

16from gpaw.response.susceptibility import ChiFactory 

17from gpaw.response.pair_functions import read_pair_function 

18 

19 

20@pytest.mark.dielectricfunction 

21@pytest.mark.kspair 

22@pytest.mark.response 

23def test_response_silicon_chi_RPA(in_tmp_dir): 

24 assert size <= 4**3 

25 

26 # Ground state calculation 

27 

28 t1 = time.time() 

29 

30 a = 5.431 

31 atoms = bulk('Si', 'diamond', a=a) 

32 atoms.center() 

33 calc = GPAW(mode=PW(200), 

34 nbands=8, 

35 kpts=(4, 4, 4), 

36 parallel={'domain': 1}, 

37 occupations=FermiDirac(width=0.05), 

38 xc='LDA') 

39 

40 atoms.calc = calc 

41 atoms.get_potential_energy() 

42 calc.write('Si', 'all') 

43 t2 = time.time() 

44 

45 # Excited state calculation 

46 q = np.array([1 / 4.0, 0, 0]) 

47 w = np.linspace(0, 24, 241) 

48 eta = 0.2 

49 

50 # Using DF 

51 df = DielectricFunction(calc='Si', 

52 frequencies=w, eta=eta, ecut=50, 

53 hilbert=False) 

54 df.get_dynamic_susceptibility(xc='RPA', q_c=q, filename='Si_chi1.csv') 

55 

56 t3 = time.time() 

57 

58 world.barrier() 

59 

60 # Using the ChiFactory 

61 gs = ResponseGroundStateAdapter(calc) 

62 chiks_calc = ChiKSCalculator(gs, ecut=50) 

63 chi_factory = ChiFactory(chiks_calc) 

64 chiks, chi = chi_factory('00', q, w + 1.j * eta) 

65 chi.write_macroscopic_component('Si_chi2.csv') 

66 chi_factory.context.write_timer() 

67 chi_factory.context.set_timer(Timer()) 

68 

69 t4 = time.time() 

70 

71 # Calculate also the ALDA susceptibility manually 

72 hxc_kernel = chi_factory.get_hxc_kernel('ALDA', '00', chiks.qpd) 

73 chi = chi_factory.dyson_solver(chiks, hxc_kernel) 

74 chi.write_macroscopic_component('Si_chi3.csv') 

75 chi_factory.context.write_timer() 

76 

77 t5 = time.time() 

78 

79 world.barrier() 

80 

81 parprint('') 

82 parprint('For ground state calc, it took', (t2 - t1) / 60, 'minutes') 

83 parprint('For excited state calc 1, it took', (t3 - t2) / 60, 'minutes') 

84 parprint('For excited state calc 2, it took', (t4 - t3) / 60, 'minutes') 

85 parprint('For excited state calc 3, it took', (t5 - t4) / 60, 'minutes') 

86 

87 w1_w, _, chi1_w = read_response_function('Si_chi1.csv') 

88 wpeak1, Ipeak1 = findpeak(w1_w, -chi1_w.imag) 

89 w2_w, chi2_w = read_pair_function('Si_chi2.csv') 

90 wpeak2, Ipeak2 = findpeak(w2_w, -chi2_w.imag) 

91 w3_w, chi3_w = read_pair_function('Si_chi3.csv') 

92 wpeak3, Ipeak3 = findpeak(w3_w, -chi3_w.imag) 

93 

94 # The two response codes should hold identical results 

95 assert wpeak1 == pytest.approx(wpeak2, abs=0.02) 

96 assert Ipeak1 == pytest.approx(Ipeak2, abs=1.0) 

97 

98 # Compare to test values 

99 assert wpeak2 == pytest.approx(16.69145, abs=0.02) # RPA 

100 assert wpeak3 == pytest.approx(16.30622, abs=0.02) # ALDA