Coverage for gpaw/test/lcaotddft/test_lcaotddft_vs_lrtddft2_rpa.py: 100%

35 statements  

« prev     ^ index     » next       coverage.py v7.7.1, created at 2025-07-20 00:19 +0000

1import numpy as np 

2import pytest 

3from ase.build import molecule 

4from gpaw import GPAW 

5from gpaw.lcaotddft import LCAOTDDFT 

6from gpaw.lcaotddft.dipolemomentwriter import DipoleMomentWriter 

7from gpaw.lrtddft2 import LrTDDFT2 

8from gpaw.mpi import world 

9from gpaw.tddft.spectrum import photoabsorption_spectrum 

10 

11 

12@pytest.mark.rttddft 

13def test_lcaotddft_lcaotddft_vs_lrtddft2_rpa(in_tmp_dir): 

14 atoms = molecule('Na2') 

15 atoms.center(vacuum=4.0) 

16 

17 # Ground-state calculation 

18 calc = GPAW(nbands=2, h=0.4, setups=dict(Na='1'), 

19 basis='sz(dzp)', mode='lcao', xc='oldLDA', 

20 convergence={'density': 1e-8}, 

21 symmetry={'point_group': False}, 

22 txt='gs.out') 

23 atoms.calc = calc 

24 atoms.get_potential_energy() 

25 calc.write('gs.gpw', mode='all') 

26 

27 # Time-propagation calculation 

28 td_calc = LCAOTDDFT('gs.gpw', fxc='RPA', txt='td.out') 

29 DipoleMomentWriter(td_calc, 'dm.dat') 

30 td_calc.absorption_kick([0, 0, 1e-5]) 

31 td_calc.propagate(30, 150) 

32 photoabsorption_spectrum('dm.dat', 'spec.dat', 

33 e_max=10, width=0.5, delta_e=0.1) 

34 

35 # LrTDDFT2 calculation 

36 calc = GPAW('gs.gpw', txt='lr.out') 

37 # It doesn't matter which fxc is here 

38 lr = LrTDDFT2('lr2', calc, fxc='LDA') 

39 lr.K_matrix.fxc_pre = 0.0 # Ignore fxc part 

40 lr.calculate() 

41 lr.get_spectrum('lr_spec.dat', 0, 10.1, 0.1, width=0.5) 

42 world.barrier() 

43 

44 # Scale spectra due to slightly different definitions 

45 spec1_ej = np.loadtxt('spec.dat') 

46 data1_e = spec1_ej[:, 3] 

47 data1_e[1:] /= spec1_ej[1:, 0] 

48 spec2_ej = np.loadtxt('lr_spec.dat') 

49 E = lr.lr_transitions.get_transitions()[0][0] 

50 data2_e = spec2_ej[:, 5] / E 

51 

52 # One can decrease the tolerance by decreasing the time step 

53 # and other parameters 

54 assert data1_e == pytest.approx(data2_e, abs=0.01)