Coverage for gpaw/test/test_inducedfield_td.py: 100%

1import pytest 

2import numpy as np 

3from ase import Atoms 

4 

5from gpaw import GPAW 

6from gpaw.tddft import TDDFT, DipoleMomentWriter 

7from gpaw.inducedfield.inducedfield_base import BaseInducedField 

8from gpaw.inducedfield.inducedfield_tddft import TDDFTInducedField 

9from gpaw.poisson import PoissonSolver 

10 

11 

12@pytest.mark.ci 

13def test_inducedfield_td(in_tmp_dir): 

14 poisson_eps = 1e-12 

15 density_eps = 1e-6 

16 

17 # PoissonSolver 

18 poissonsolver = PoissonSolver('fd', eps=poisson_eps) 

19 

20 # Na2 cluster 

21 atoms = Atoms(symbols='Na2', 

22 positions=[(0, 0, 0), (3.0, 0, 0)], 

23 pbc=False) 

24 atoms.center(vacuum=3.0) 

25 

26 # Standard ground state calculation 

27 calc = GPAW(mode='fd', nbands=2, 

28 h=0.6, 

29 setups={'Na': '1'}, 

30 poissonsolver=poissonsolver, 

31 symmetry={'point_group': False}, 

32 convergence={'density': density_eps}) 

33 atoms.calc = calc 

34 _ = atoms.get_potential_energy() 

35 calc.write('na2_gs.gpw', mode='all') 

36 

37 # Standard time-propagation initialization 

38 time_step = 10.0 

39 iterations = 20 

40 kick_strength = [1.0e-3, 1.0e-3, 0.0] 

41 td_calc = TDDFT('na2_gs.gpw') 

42 DipoleMomentWriter(td_calc, 'na2_td_dm.dat') 

43 

44 # Create and attach InducedField object 

45 frequencies = [1.0, 2.08] # Frequencies of interest in eV 

46 folding = 'Gauss' # Folding function 

47 width = 0.1 # Line width for folding in eV 

48 ind = TDDFTInducedField(paw=td_calc, 

49 frequencies=frequencies, 

50 folding=folding, 

51 width=width, 

52 restart_file='na2_td.ind') 

53 

54 # Propagate as usual 

55 td_calc.absorption_kick(kick_strength=kick_strength) 

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

57 

58 # Save TDDFT and InducedField objects 

59 td_calc.write('na2_td.gpw', mode='all') 

60 ind.write('na2_td.ind') 

61 ind.paw = None 

62 

63 # Restart and continue 

64 td_calc = TDDFT('na2_td.gpw') 

65 DipoleMomentWriter(td_calc, 'na2_td_dm.dat') 

66 

67 # Load and attach InducedField object 

68 ind = TDDFTInducedField(filename='na2_td.ind', 

69 paw=td_calc, 

70 restart_file='na2_td.ind') 

71 

72 # Continue propagation as usual 

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

74 

75 # Calculate induced electric field 

76 ind.calculate_induced_field(gridrefinement=2, 

77 from_density='comp', 

78 poissonsolver=poissonsolver, 

79 extend_N_cd=3 * np.ones((3, 2), int), 

80 deextend=True) 

81 

82 # Save 

83 ind.write('na2_td_field.ind', 'all') 

84 ind.paw = None 

85 td_calc = None 

86 ind = None 

87 

88 # Read data (test also field data I/O) 

89 ind = BaseInducedField(filename='na2_td_field.ind', 

90 readmode='field') 

91 

92 # Estimate tolerance (worst case error accumulation) 

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

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

95 # tol = 0.038905993684 

96 # print('tol = %.12f' % tol) 

97 

98 # Test 

99 val1 = ind.fieldgd.integrate(ind.Ffe_wg[0]) 

100 val2 = ind.fieldgd.integrate(np.abs(ind.Fef_wvg[0][0])) 

101 val3 = ind.fieldgd.integrate(np.abs(ind.Fef_wvg[0][1])) 

102 val4 = ind.fieldgd.integrate(np.abs(ind.Fef_wvg[0][2])) 

103 val5 = ind.fieldgd.integrate(ind.Ffe_wg[1]) 

104 val6 = ind.fieldgd.integrate(np.abs(ind.Fef_wvg[1][0])) 

105 val7 = ind.fieldgd.integrate(np.abs(ind.Fef_wvg[1][1])) 

106 val8 = ind.fieldgd.integrate(np.abs(ind.Fef_wvg[1][2])) 

107 

108 assert val1 == pytest.approx(1926.232999117403, abs=tol) 

109 assert val2 == pytest.approx(0.427606450419, abs=tol) 

110 assert val3 == pytest.approx(0.565823985683, abs=tol) 

111 assert val4 == pytest.approx(0.372493489423, abs=tol) 

112 assert val5 == pytest.approx(1945.618902611449, abs=tol) 

113 assert val6 == pytest.approx(0.423899965987, abs=tol) 

114 assert val7 == pytest.approx(0.560882533828, abs=tol) 

115 assert val8 == pytest.approx(0.369203021329, abs=tol)