Coverage for gpaw/lrtddft/convergence.py: 4%

1import os 

2 

3from ase.units import Hartree 

4 

5from gpaw.lrtddft.spectrum import spectrum, rotatory_spectrum 

6 

7 

8def check_convergence(lr, # LrTDDFT object 

9 dirname='conv', # directory name to store the files 

10 title=None, # title for gnuplot 

11 dn=None, # steps to vary istart/jend 

12 dE=None, # steps to vary energy range 

13 emin=None, 

14 emax=None, 

15 folding='Gauss', 

16 istart=None, 

17 jend=None, 

18 width=0.03): # Gauss/Lorentz width 

19 """Study the convergence of a LrTDDFT calcualtion by varying istart/jend. 

20 A gnuplot file will be created with the name 'dirname'/conv.gpl.""" 

21 

22 if istart is None: 

23 istart0 = lr.kss.restrict['istart'] 

24 else: 

25 if istart < lr.kss.restrict['istart']: 

26 raise RuntimeError 

27 istart0 = istart 

28 if jend is None: 

29 jend0 = lr.kss.restrict['jend'] 

30 else: 

31 if jend > lr.kss.restrict['jend']: 

32 raise RuntimeError 

33 jend0 = jend 

34 

35 # create subdirectory for the files 

36 if not os.path.isdir(dirname): 

37 if not os.path.exists(dirname): 

38 os.makedirs(dirname) 

39 else: 

40 raise RuntimeError('Can\'t create directory ' + dirname) 

41 

42 def fname(filename): 

43 return dirname + '/' + filename 

44 

45 with open(fname('conv.gpl'), 'w') as fgpl: 

46 print('set xlabel "omega [eV]"', file=fgpl) 

47 print('set ylabel "Folded osc. strength [1/eV]"', file=fgpl) 

48 if not emin: 

49 emin_gpl = '*' 

50 else: 

51 emin_gpl = str(emin) 

52 if not emax: 

53 emax_gpl = '*' 

54 else: 

55 emax_gpl = str(emax) 

56 print('set xrange [' + emin_gpl + ':' + emax_gpl + ']', file=fgpl) 

57 if title: 

58 print('set title "' + str(title) + '"', file=fgpl) 

59 

60 # kss 

61 spectrum(lr.kss, fname('kss.dat'), width=width) 

62 spectrum(lr.kss, fname('ksssticks.dat'), folding=None) 

63 

64 # full 

65 lr.diagonalize(restrict=dict(istart=istart0, jend=jend0)) 

66 spectrum(lr, fname('full.dat'), width=width) 

67 spectrum(lr, fname('fullsticks.dat'), folding=None) 

68 print('plot "' + fname('full.dat') + '" t "full" w l lt 1, \\', 

69 file=fgpl) 

70 print(' "' + fname('fullsticks.dat') + 

71 '" u 1:($2*20) t "" w impulses lt 1, \\', file=fgpl) 

72 print(' "' + fname('kss.dat') + '" t "Kohn-Sham" w l lt 2', 

73 file=fgpl) 

74 print('pause -10', file=fgpl) 

75 

76 if dn is None: 

77 dn = -istart0 + jend0 

78 dn = int(dn / 10.) 

79 

80 if dE is None: 

81 # vary istart 

82 print('plot "' + fname('full.dat') + '" t "istart=' + 

83 str(istart0) + '" w l lt 1, \\', file=fgpl) 

84 for i in range(1, 4): 

85 istart = istart0 + i * dn 

86 lr.diagonalize(restrict=dict(istart=istart, jend=jend0)) 

87 fn = fname('istart' + str(istart) + '.dat') 

88 spectrum(lr, fn, width=width) 

89 print(' "' + fn + '" t "istart=' + str(istart) + 

90 '" w l lt', i + 1, end=' ', file=fgpl) 

91 if i < 3: 

92 print(', \\', end=' ', file=fgpl) 

93 print(file=fgpl) 

94 print('pause -10', file=fgpl) 

95 

96 # vary jend 

97 print('plot "' + fname('full.dat') + '" t "jend=' + 

98 str(jend0) + '" w l lt 1, \\', file=fgpl) 

99 for i in range(1, 4): 

100 jend = jend0 - i * dn 

101 lr.diagonalize(restrict=dict(jend=jend, istart=istart0)) 

102 fn = fname('jend' + str(jend) + '.dat') 

103 spectrum(lr, fn, width=width) 

104 print(' "' + fn + '" t "jend=' + str(jend) + 

105 '" w l lt', i + 1, end=' ', file=fgpl) 

106 if i < 3: 

107 print(', \\', end=' ', file=fgpl) 

108 print(file=fgpl) 

109 print('pause -10', file=fgpl) 

110 else: 

111 # vary the enrgy range 

112 max_kss_energy = 0. 

113 for kss in lr.Om.kss: 

114 max_kss_energy = max(max_kss_energy, 

115 kss.get_energy() * Hartree) 

116 print('plot "' + fname('full.dat') + '" t "full"' + 

117 ' w l lt 1, \\', file=fgpl) 

118 for i in range(1, 4): 

119 max_energy = max_kss_energy - i * dE 

120 lr.diagonalize(restrict=dict(energy_range=max_energy)) 

121 fn = fname('max_energy' + str(max_energy) + '.dat') 

122 spectrum(lr, fn, width=width) 

123 print(' "' + fn + '" t "dE=-' + str( 

124 i * dE) + '" w l lt', end=' ', file=fgpl) 

125 print(i + 1, end=' ', file=fgpl) 

126 if i < 3: 

127 print(', \\', end=' ', file=fgpl) 

128 print(file=fgpl) 

129 print('pause -10', file=fgpl) 

130 

131 # plot different directions 

132 print('plot "' + fname('full.dat') + 

133 '" u 1:3 t "x" w l lt 1, \\', file=fgpl) 

134 print(' "' + fname('full.dat') + 

135 '" u 1:4 t "y" w l lt 2, \\', file=fgpl) 

136 print(' "' + fname('full.dat') + '" u 1:5 t "z" w l lt 3', 

137 file=fgpl) 

138 print('pause -10', file=fgpl) 

139 

140 # plot rotary strength 

141 if lr[0].magn is not None: 

142 print('set ylabel "Folded rot. strength [cgs/eV]"', file=fgpl) 

143 rotatory_spectrum(lr, fname('rotatory.dat'), width=width) 

144 rotatory_spectrum(lr, fname('rotatory_sticks.dat'), folding=None) 

145 print('plot "' + fname('rotatory.dat') + 

146 '" t "rotatory" w l lt 1, \\', file=fgpl) 

147 print(' "' + fname('rotatory_sticks.dat') + 

148 '" u 1:($2*20) t "" w impulses lt 1', file=fgpl) 

149 print('pause -10', file=fgpl)