Coverage for gpaw/test/lcaotddft/demo_tddft.py: 19%

1import numpy as np 

2 

3from ase import Atoms 

4 

5from gpaw import LCAO 

6from gpaw.calculator import GPAW as old_GPAW 

7from gpaw.lcaotddft import LCAOTDDFT 

8from gpaw.lcaotddft.dipolemomentwriter import DipoleMomentWriter 

9from gpaw.new.ase_interface import GPAW as new_GPAW 

10from gpaw.new.rttddft import RTTDDFT 

11from gpaw.tddft.units import as_to_au, autime_to_asetime 

12 

13 

14def main(): 

15 atoms = Atoms('H2', positions=[(0, 0, 0), (1, 0, 0)]) 

16 atoms.center(vacuum=5) 

17 atoms.pbc = False 

18 

19 kick_v = [1e-5, 0, 0] 

20 

21 run_old_gs = True 

22 run_new_gs = True 

23 run_old_td = True 

24 run_new_td = True 

25 

26 def assert_equal(a, b): 

27 from gpaw.core.matrix import Matrix 

28 from gpaw.core.atom_arrays import AtomArrays 

29 

30 def extract(o): 

31 if isinstance(o, Matrix): 

32 return o.data 

33 elif isinstance(o, AtomArrays): 

34 return o.data 

35 else: 

36 return o 

37 

38 a = extract(a) 

39 b = extract(b) 

40 

41 assert np.allclose(a, b), f'{str(a)} != {str(b)}' 

42 

43 if run_old_gs: 

44 old_calc = old_GPAW(mode=LCAO(), basis='sz(dzp)', xc='LDA', 

45 symmetry={'point_group': False}, 

46 txt='old.out', convergence={'density': 1e-12}) 

47 atoms.calc = old_calc 

48 atoms.get_potential_energy() 

49 old_calc.write('old_gs.gpw', mode='all') 

50 

51 if run_new_gs: 

52 new_calc = new_GPAW(mode='lcao', basis='sz(dzp)', xc='LDA', 

53 txt='new.out', force_complex_dtype=True, 

54 convergence={'density': 1e-12}) 

55 atoms.calc = new_calc 

56 atoms.get_potential_energy() 

57 new_calc.write('new_gs.gpw', mode='all') 

58 

59 new_restart_calc = new_GPAW('new_gs.gpw') 

60 

61 assert_equal( 

62 new_calc.dft.ibzwfs.wfs_qs[0][0].P_ain, 

63 new_restart_calc.dft.ibzwfs.wfs_qs[0][0].P_ain) 

64 

65 assert_equal( 

66 new_calc.dft.ibzwfs.wfs_qs[0][0].C_nM, 

67 new_restart_calc.dft.ibzwfs.wfs_qs[0][0].C_nM) 

68 

69 if run_old_td: 

70 old_tddft = LCAOTDDFT('old_gs.gpw', propagator='ecn', txt='/dev/null') 

71 DipoleMomentWriter(old_tddft, 'old_dm.out') 

72 old_tddft.absorption_kick(kick_v) 

73 old_tddft.propagate(10, 10) 

74 # old_C_nM = old_tddft.wfs.kpt_u[0].C_nM 

75 # old_f_n = old_tddft.get_occupation_numbers() 

76 # old_rho_MM = old_C_nM.T.conj() @ (old_f_n[:, None] * old_C_nM) 

77 # print('rho_MM', old_rho_MM) 

78 

79 if run_new_td: 

80 # new_tddft = RTTDDFT.from_dft_calculation(new_calc) 

81 new_tddft = RTTDDFT.from_dft_file('new_gs.gpw') 

82 

83 new_tddft.absorption_kick(kick_v) 

84 dt = 10 * as_to_au * autime_to_asetime 

85 with open('new_dm.out', 'w') as fp: 

86 for result in new_tddft.ipropagate(dt, 10): 

87 dm = result.dipolemoment 

88 fp.write('%20.8lf %20.8le %22.12le %22.12le %22.12le\n' % 

89 (result.time, 0, dm[0], dm[1], dm[2])) 

90 print(result) 

91 # wfs = new_tddft.state.ibzwfs.wfs_qs[0][0] 

92 # new_rho_MM = wfs.calculate_density_matrix() 

93 # print('rho_MM', new_rho_MM) 

94 

95 

96if __name__ == '__main__': 

97 main()