Coverage for gpaw/atom/analyse_setup.py: 2%

1import numpy as np 

2 

3 

4def analyse(generator, show=False): 

5 import matplotlib.pyplot as plt 

6 gen = generator 

7 symbol = gen.symbol 

8 

9 colors = [] 

10 id_j = [] 

11 for j, n in enumerate(gen.vn_j): 

12 if n == -1: 

13 n = '*' 

14 id_j.append(str(n) + 'spdf'[gen.vl_j[j]]) 

15 colors.append('krgbymc'[j]) 

16 

17 r_g = gen.r 

18 rmax = max(gen.rcut_l) 

19 

20 # Construct logarithmic derivatives 

21 if len(gen.logd) > 0: 

22 rlog = gen.rlog 

23 logd = gen.logd 

24 elog = gen.elog 

25 ref = [] 

26 for l, e in zip(gen.vl_j, gen.ve_j): 

27 i = elog.searchsorted(e) 

28 ref.append((elog[i], logd[l][1][i])) 

29 ref = np.array(ref) 

30 

31 dpi = 80 

32 fig = plt.figure(figsize=(8.0, 12.0), dpi=dpi) 

33 fig.subplots_adjust(left=0.1, right=0.99, top=0.97, bottom=0.04, 

34 hspace=0.25, wspace=0.26) 

35 

36 plt.subplot(321) 

37 ymax = 0 

38 s = 1.0 

39 g = r_g.searchsorted(rmax * 2) 

40 for phi_g, phit_g, id, color in zip(gen.vu_j, gen.vs_j, 

41 id_j, colors): 

42 if id[0] != '*': 

43 lw = 2 

44 ymax = max(ymax, phi_g.max(), -phi_g.min()) 

45 else: 

46 lw = 1 

47 s = ymax / max(np.abs(phi_g[:g])) 

48 plt.plot(r_g, s * phi_g, color + '-', label=id, lw=lw) 

49 plt.plot(r_g, s * phit_g, color + '--', label='_nolegend_', lw=lw) 

50 plt.legend(loc='lower right') 

51 plt.axis('tight') 

52 lim = plt.axis(xmin=0, xmax=rmax * 2, ymin=-ymax, ymax=ymax) 

53 plt.plot([rmax, rmax], lim[2:], 'k--', label='_nolegend_') 

54 plt.text(rmax, lim[2], r'$r_c$', ha='left', va='bottom', size=17) 

55 plt.title('Partial Waves') 

56 plt.xlabel('r [Bohr]') 

57 plt.ylabel(r'$r\ \phi,\ r\tilde{\phi}$, [Bohr$^{-1/2}$]') 

58 

59 plt.subplot(322) 

60 ymax = 0 

61 s = 1.0 

62 for pt_g, id, color in zip(gen.vq_j, id_j, colors): 

63 if id[0] != '*': 

64 lw = 2 

65 ymax = max(ymax, pt_g.max(), -pt_g.min()) 

66 else: 

67 lw = 1 

68 s = ymax / max(np.abs(pt_g)) 

69 plt.plot(r_g, s * pt_g, color + '-', label=id, lw=lw) 

70 plt.axis('tight') 

71 lim = plt.axis(xmin=0, xmax=rmax * 1.2, ymin=-ymax, ymax=ymax) 

72 plt.plot([rmax, rmax], lim[2:], 'k--', label='_nolegend_') 

73 plt.text(rmax, lim[2], r'$r_c$', ha='left', va='bottom', size=17) 

74 plt.legend(loc='best') 

75 plt.title('Projectors') 

76 plt.xlabel('r [Bohr]') 

77 plt.ylabel(r'$r\ \tilde{p}$, [Bohr$^{-1/2}$]') 

78 

79 plt.subplot(323) 

80 plt.plot(r_g, gen.nc, colors[0], label=r'$n_c$') 

81 plt.plot(r_g, gen.nct, colors[1], label=r'$\tilde{n}_c$') 

82 plt.axis('tight') 

83 lim = plt.axis(xmin=0, xmax=rmax * 1.2, ymax=max(gen.nct)) 

84 plt.plot([rmax, rmax], lim[2:], 'k--', label='_nolegend_') 

85 plt.text(rmax, lim[2], r'$r_c$', ha='left', va='bottom', size=17) 

86 plt.legend(loc='best') 

87 plt.title('Densities') 

88 plt.xlabel('r [Bohr]') 

89 plt.ylabel(r'density [Bohr$^{-3}$]') 

90 

91 plt.subplot(324) 

92 plt.plot(r_g[1:], gen.vr[1:] / r_g[1:], label=r'$v$') 

93 plt.plot(r_g, gen.vt, label=r'$\tilde{v}$') 

94 plt.plot(r_g, gen.vt - gen.vbar, label=r'$\tilde{v}-\bar{v}$') 

95 plt.axis('tight') 

96 lim = plt.axis(xmin=0, xmax=rmax * 1.2, 

97 ymin=(gen.vt - gen.vbar).min(), ymax=0) 

98 plt.plot([rmax, rmax], lim[2:], 'k--', label='_nolegend_') 

99 plt.text(rmax, lim[2], r'$r_c$', ha='left', va='bottom', size=17) 

100 plt.legend(loc='best') 

101 plt.title('Potentials') 

102 plt.xlabel('r [Bohr]') 

103 plt.ylabel('potential [Hartree]') 

104 

105 plt.subplot(325) 

106 if len(gen.logd) > 0: 

107 plt.plot(ref[:, 0], ref[:, 1], 'ko', label='_nolegend_') 

108 for l, color in enumerate(colors[:4]): 

109 id = 'spdf'[l] 

110 plt.plot(elog, logd[l][0], linestyle='-', color=color, label=id) 

111 plt.plot(elog, logd[l][1], linestyle='--', color=color, 

112 label='_nolegend_') 

113 plt.ylabel('log. deriv. at r=%.2f Bohr' % rlog) 

114 ymin = ref[:, 1].min() 

115 ymax = ref[:, 1].max() 

116 plt.axis(ymin=ymin - (ymax - ymin) * 0.1, 

117 ymax=ymax + (ymax - ymin) * 0.1) 

118 plt.legend(loc='best') 

119 plt.title('Logarithmic Derivatives') 

120 plt.xlabel('Energy [Hartree]') 

121 

122 plt.savefig('%s-setup.png' % symbol, dpi=dpi) 

123 

124 if show: 

125 plt.show()