Coverage for gpaw/test/fdtd/test_ed_inducedfield.py: 98%

1import numpy as np 

2import pytest 

3from gpaw import GPAW 

4from gpaw.inducedfield.inducedfield_base import BaseInducedField 

5from gpaw.inducedfield.inducedfield_fdtd import ( 

6 FDTDInducedField, calculate_hybrid_induced_field) 

7from gpaw.inducedfield.inducedfield_tddft import TDDFTInducedField 

8from gpaw.tddft import TDDFT, DipoleMomentWriter 

9 

10 

11@pytest.mark.old_gpaw_only 

12def test_fdtd_ed_inducedfield(in_tmp_dir, gpw_files): 

13 do_print_values = 0 # Use this for printing the reference values 

14 

15 # Accuracy 

16 energy_eps = 0.0005 

17 poisson_eps = 1e-12 

18 density_eps = 1e-6 

19 

20 # load gpw file 

21 gs_calc = GPAW(gpw_files['na2_isolated']) 

22 energy = gs_calc.get_potential_energy() 

23 

24 # Test ground state 

25 assert energy == pytest.approx( 

26 -0.631881, 

27 abs=energy_eps * gs_calc.get_number_of_electrons()) 

28 

29 # Test floating point arithmetic errors 

30 assert gs_calc.hamiltonian.poisson.shift_indices_1 == pytest.approx( 

31 [4, 4, 10], abs=0) 

32 assert gs_calc.hamiltonian.poisson.shift_indices_2 == pytest.approx( 

33 [8, 8, 16], abs=0) 

34 

35 # Initialize TDDFT and FDTD 

36 kick = [0.0, 0.0, 1.0e-3] 

37 time_step = 10.0 

38 iterations = 10 

39 

40 td_calc = TDDFT(gpw_files['na2_isolated']) 

41 DipoleMomentWriter(td_calc, 'dm.dat') 

42 

43 # Attach InducedFields to the calculation 

44 frequencies = [2.05, 4.0] 

45 width = 0.15 

46 cl_ind = FDTDInducedField(paw=td_calc, 

47 frequencies=frequencies, 

48 width=width) 

49 qm_ind = TDDFTInducedField(paw=td_calc, 

50 frequencies=frequencies, 

51 width=width) 

52 

53 # Propagate TDDFT and FDTD 

54 td_calc.absorption_kick(kick_strength=kick) 

55 td_calc.propagate(time_step, iterations // 2) 

56 td_calc.write('td.gpw', 'all') 

57 cl_ind.write('cl.ind') 

58 qm_ind.write('qm.ind') 

59 

60 # Restart 

61 td_calc = TDDFT('td.gpw') 

62 DipoleMomentWriter(td_calc, 'dm.dat') 

63 cl_ind = FDTDInducedField(filename='cl.ind', 

64 paw=td_calc) 

65 qm_ind = TDDFTInducedField(filename='qm.ind', 

66 paw=td_calc) 

67 td_calc.propagate(time_step, iterations // 2) 

68 td_calc.write('td.gpw', 'all') 

69 cl_ind.write('cl.ind') 

70 qm_ind.write('qm.ind') 

71 

72 # Test 

73 ref_cl_dipole_moment = [5.25374117e-14, 5.75811267e-14, 3.08349334e-02] 

74 ref_qm_dipole_moment = [1.78620337e-11, -1.57782578e-11, 5.21368300e-01] 

75 

76 tol = 1e-4 

77 assert td_calc.hamiltonian.poisson.get_classical_dipole_moment() == ( 

78 pytest.approx(ref_cl_dipole_moment, abs=tol)) 

79 assert td_calc.hamiltonian.poisson.get_quantum_dipole_moment() == ( 

80 pytest.approx(ref_qm_dipole_moment, abs=tol)) 

81 

82 # Calculate induced fields 

83 td_calc = TDDFT('td.gpw') 

84 

85 # Classical subsystem 

86 cl_ind = FDTDInducedField(filename='cl.ind', paw=td_calc) 

87 cl_ind.calculate_induced_field(gridrefinement=2) 

88 cl_ind.write('cl_field.ind', mode='all') 

89 

90 # Quantum subsystem 

91 qm_ind = TDDFTInducedField(filename='qm.ind', paw=td_calc) 

92 qm_ind.calculate_induced_field(gridrefinement=2) 

93 qm_ind.write('qm_field.ind', mode='all') 

94 

95 # Total system, interpolate/extrapolate to a grid with spacing h 

96 tot_ind = calculate_hybrid_induced_field(cl_ind, qm_ind, h=0.4) 

97 tot_ind.write('tot_field.ind', mode='all') 

98 

99 # Test induced fields 

100 ref_values = [72404.467117024149, 

101 0.520770766296, 

102 0.520770766299, 

103 0.830247064075, 

104 72416.234345610734, 

105 0.517294132489, 

106 0.517294132492, 

107 0.824704513888, 

108 2451.767847927681, 

109 0.088037476748, 

110 0.088037476316, 

111 0.123334033914, 

112 2454.462292798476, 

113 0.087537484422, 

114 0.087537483971, 

115 0.122592730690, 

116 76582.089818637178, 

117 0.589941751987, 

118 0.589941751804, 

119 0.869526245360, 

120 76592.175846021099, 

121 0.586223386358, 

122 0.586223386102, 

123 0.864478308364] 

124 

125 for fname in ['cl_field.ind', 'qm_field.ind', 'tot_field.ind']: 

126 ind = BaseInducedField(filename=fname, 

127 readmode='field') 

128 # Estimate tolerance (worst case error accumulation) 

129 tol = (iterations * ind.fieldgd.integrate(ind.fieldgd.zeros() + 1.0) * 

130 max(density_eps, np.sqrt(poisson_eps))) 

131 if do_print_values: 

132 print('tol = %.12f' % tol) 

133 for w in range(len(frequencies)): 

134 val = ind.fieldgd.integrate(ind.Ffe_wg[w]) 

135 assert val == pytest.approx(ref_values.pop(0), abs=tol) 

136 for v in range(3): 

137 val = ind.fieldgd.integrate(np.abs(ind.Fef_wvg[w][v])) 

138 assert val == pytest.approx(ref_values.pop(0), abs=tol)