Coverage for gpaw/test/response/test_nicl2_pair_potential.py: 95%

1import pytest 

2 

3import numpy as np 

4 

5from gpaw import GPAW 

6from gpaw.mpi import world 

7from gpaw.response import ResponseContext, ResponseGroundStateAdapter 

8from gpaw.response.pw_parallelization import block_partition 

9from gpaw.response.matrix_elements import TransversePairPotentialCalculator 

10 

11from gpaw.test.response.test_parallel_kptpair_extraction import \ 

12 initialize_extractor, initialize_transitions, initialize_integral 

13 

14 

15@pytest.mark.response 

16@pytest.mark.kspair 

17@pytest.mark.old_gpaw_only 

18def test_nicl2_pair_potential(gpw_files): 

19 """Test that the transverse pair potential vanishes in vacuum.""" 

20 

21 # ---------- Inputs ---------- # 

22 

23 wfs = 'nicl2_pw_evac' 

24 q_qc = [[0., 0., 0.], 

25 [1. / 3., 1. / 3., 0.]] 

26 rshewmin = 1e-8 

27 nblocks = world.size // 2 if world.size % 2 == 0 else 1 

28 

29 # ---------- Script ---------- # 

30 

31 context = ResponseContext() 

32 calc = GPAW(gpw_files[wfs], parallel=dict(domain=1)) 

33 gs = ResponseGroundStateAdapter(calc) 

34 

35 # Set up extractor and transitions 

36 tcomm, kcomm = block_partition(context.comm, nblocks) 

37 extractor = initialize_extractor(gs, context, tcomm, kcomm) 

38 transitions = initialize_transitions(gs, '+-', nbands=20) 

39 

40 # Set up calculator 

41 pair_potential_calc = TransversePairPotentialCalculator(gs, context, 

42 rshewmin=rshewmin) 

43 

44 # Loop over k-point pairs 

45 for q_c in q_qc: 

46 integral = initialize_integral(extractor, context, q_c) 

47 for kptpair, weight in integral.weighted_kpoint_pairs(transitions): 

48 if weight is None: 

49 assert kptpair is None 

50 continue 

51 _, _, ut1_mytR, ut2_mytR = \ 

52 pair_potential_calc.extract_pseudo_waves(kptpair) 

53 wt_mytR = \ 

54 pair_potential_calc._evaluate_pseudo_matrix_element( 

55 ut1_mytR, ut2_mytR) 

56 

57 # Average out the in-plane degrees of freedom 

58 wt_mytz = np.average(wt_mytR, axis=(1, 2)) 

59 

60 # Gather all the transitions 

61 wt_tz = kptpair.tblocks.all_gather(wt_mytz) 

62 

63 # import matplotlib.pyplot as plt 

64 # for wt_z in wt_tz: 

65 # plt.subplot(1, 2, 1) 

66 # plt.plot(np.arange(len(wt_z)), wt_z.real) 

67 # plt.subplot(1, 2, 2) 

68 # plt.plot(np.arange(len(wt_z)), wt_z.imag) 

69 # plt.show() 

70 

71 # Find the maximum absolute value of the pair potential and make 

72 # sure that the pair potential is much smaller than that close to 

73 # the cell boundary for all the transitions 

74 abswt_tz = np.abs(wt_tz) 

75 max_wt = np.max(abswt_tz) 

76 assert np.all(abswt_tz[:, :10] < max_wt / 25.) 

77 assert np.all(abswt_tz[:, -10:] < max_wt / 25.)