Coverage for gpaw/test/test_inducedfield_lrtddft.py: 98%

1import numpy as np 

2import pytest 

3from ase import Atoms 

4 

5from gpaw import GPAW 

6from gpaw.inducedfield.inducedfield_lrtddft import LrTDDFTInducedField 

7from gpaw.lrtddft import LrTDDFT 

8from gpaw.poisson import FDPoissonSolver 

9 

10 

11@pytest.mark.ci 

12@pytest.mark.old_gpaw_only 

13def test_inducedfield_lrtddft(in_tmp_dir): 

14 do_print_values = False # Use this for printing the reference values 

15 poisson_eps = 1e-12 

16 density_eps = 1e-6 

17 

18 # 0) PoissonSolver 

19 poissonsolver = FDPoissonSolver(eps=poisson_eps) 

20 

21 # 1) Ground state calculation with empty states 

22 atoms = Atoms(symbols='Na2', 

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

24 pbc=False) 

25 atoms.center(vacuum=3.0) 

26 

27 calc = GPAW(mode='fd', nbands=20, h=0.6, setups={'Na': '1'}, 

28 poissonsolver=poissonsolver, 

29 convergence={'density': density_eps}) 

30 atoms.calc = calc 

31 atoms.get_potential_energy() 

32 calc.write('na2_gs_casida.gpw', mode='all') 

33 

34 # 2) Casida calculation 

35 calc = GPAW('na2_gs_casida.gpw') 

36 istart = 0 

37 jend = 20 

38 lr = LrTDDFT(calc, xc='LDA', 

39 restrict={'istart': istart, 'jend': jend}) 

40 lr.diagonalize() 

41 lr.write('na2_lr.dat.gz') 

42 

43 # Start from scratch 

44 del lr 

45 del calc 

46 calc = GPAW('na2_gs_casida.gpw') 

47 # calc.initialize_positions() 

48 # calc.set_positions() 

49 lr = LrTDDFT.read('na2_lr.dat.gz') 

50 

51 # 3) Calculate induced field 

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

53 folding = 'Gauss' # Folding function 

54 width = 0.1 # Line width for folding in eV 

55 kickdir = 0 # Kick field direction 0, 1, 2 for x, y, z 

56 ind = LrTDDFTInducedField(paw=calc, lr=lr, frequencies=frequencies, 

57 folding=folding, width=width, kickdir=kickdir) 

58 ind.calculate_induced_field(gridrefinement=2, 

59 from_density='comp', 

60 poissonsolver=poissonsolver, 

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

62 deextend=True) 

63 

64 # Estimate tolerance (worst case error accumulation) 

65 tol = (len(lr) ** 2 * ind.fieldgd.integrate(ind.fieldgd.zeros() + 1.0) * 

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

67 # tol = 0.702253185994 

68 if do_print_values: 

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

70 

71 # Test 

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

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

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

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

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

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

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

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

80 

81 assert val1 == pytest.approx(3175.73307601948, abs=tol) 

82 assert val2 == pytest.approx(1700.88549123764, abs=tol) 

83 assert val3 == pytest.approx(1187.07316974877, abs=tol) 

84 assert val4 == pytest.approx(1187.07316974844, abs=tol) 

85 assert val5 == pytest.approx(10957.7042391289, abs=tol) 

86 assert val6 == pytest.approx(6576.73616797380, abs=tol) 

87 assert val7 == pytest.approx(4589.04913402763, abs=tol) 

88 assert val8 == pytest.approx(4589.04913402650, abs=tol)