Coverage for gpaw/test/lrtddft/test_pes.py: 100%

1import pytest 

2from ase import Atom, Atoms 

3from ase.parallel import parprint 

4 

5from gpaw import GPAW, mpi 

6from gpaw.lrtddft import LrTDDFT 

7from gpaw.poisson import FDPoissonSolver 

8from gpaw.pes.dos import DOSPES 

9from gpaw.pes.tddft import TDDFTPES 

10 

11 

12@pytest.mark.lrtddft 

13def test_lrtddft_pes(in_tmp_dir): 

14 txt = None 

15 R = 0.7 # approx. experimental bond length 

16 a = 3.0 

17 c = 3.0 

18 H2 = Atoms([Atom('H', (a / 2, a / 2, (c - R) / 2)), 

19 Atom('H', (a / 2, a / 2, (c + R) / 2))], 

20 cell=(a, a, c)) 

21 

22 H2_plus = Atoms([Atom('H', (a / 2, a / 2, (c - R) / 2)), 

23 Atom('H', (a / 2, a / 2, (c + R) / 2))], 

24 cell=(a, a, c)) 

25 

26 xc = 'LDA' 

27 calc = GPAW(mode='fd', gpts=(12, 12, 12), xc=xc, nbands=1, 

28 poissonsolver=FDPoissonSolver(), 

29 parallel={'domain': mpi.world.size}, 

30 spinpol=True, txt=txt) 

31 H2.calc = calc 

32 e_H2 = H2.get_potential_energy() 

33 

34 calc_plus = GPAW(mode='fd', gpts=(12, 12, 12), xc=xc, nbands=2, 

35 poissonsolver=FDPoissonSolver(), 

36 parallel={'domain': mpi.world.size}, 

37 spinpol=True, txt=txt) 

38 calc_plus = calc_plus.new(charge=+1) 

39 H2_plus.calc = calc_plus 

40 e_H2_plus = H2_plus.get_potential_energy() 

41 

42 out = 'dospes.dat' 

43 pes = DOSPES(calc, calc_plus, shift=True) 

44 pes.save_folded_pes(filename=out, folding=None) 

45 parprint('DOS:') 

46 pes.save_folded_pes(filename=None, folding=None) 

47 

48 # check for correct shift 

49 VDE = calc_plus.get_potential_energy() - calc.get_potential_energy() 

50 BE_HOMO = 1.e23 

51 be_n, f_n = pes.get_energies_and_weights() 

52 for be, f in zip(be_n, f_n): 

53 if f > 0.1 and be < BE_HOMO: 

54 BE_HOMO = be 

55 assert BE_HOMO == pytest.approx(VDE, abs=0.0) 

56 

57 lr = LrTDDFT(calc_plus, xc=xc) 

58 

59 out = 'lrpes.dat' 

60 pes = TDDFTPES(calc, lr) 

61 pes.save_folded_pes(filename=out, folding='Gauss') 

62 parprint('Linear response:') 

63 pes.save_folded_pes(filename=None, folding=None) 

64 

65 energy_tolerance = 0.001 

66 assert e_H2 == pytest.approx(-3.90059, abs=energy_tolerance) 

67 assert e_H2_plus == pytest.approx(10.5659703, abs=energy_tolerance) 

68 

69 # io 

70 out = 'lrpes.dat.gz' 

71 lr.write(out) 

72 lr = LrTDDFT.read(out) 

73 lr.calculator = calc_plus 

74 

75 pes = TDDFTPES(calc, lr) 

76 parprint('Linear response:') 

77 pes.save_folded_pes(filename=None, folding=None)