r m x p   toggle line displays

j k   next/prev highlighted chunk

0   (zero) top of page

1   (one) first highlighted chunk

1import numpy as np

3from ase.units import Hartree, Bohr

6class Excitation:

7 """Base class for a single excitation"""

9 def __init__(self, energy, index, mur, muv=None, magn=None):

10 """

11 Parameters

12 ----------

13 energy: float

14 Energy realtive to the ground state

15 index: int

16 Excited state index

17 mur: list of three floats or complex numbers

18 Length form dipole matrix element

19 muv: list of three floats or complex numbers or None

20 Velocity form dipole matrix element, default None

21 magn: list of three floats or complex numbers or None

22 Magnetic matrix element, default None

23 """

24 self.energy = energy

25 self.index = index

26 self.mur = mur

27 self.muv = muv

28 self.magn = magn

29 self.fij = 1.

31 def outstring(self):

32 """Format yourself as a string"""

33 string = '{0:g} {1} '.format(self.energy, self.index)

35 def format_me(me):

36 string = ''

37 if me.dtype == float:

38 for m in me:

39 string += ' {0:g}'.format(m)

40 else:

41 for m in me:

42 string += ' {0.real:g}{0.imag:+g}j'.format(m)

43 return string

45 string += ' ' + format_me(self.mur)

46 if self.muv is not None:

47 string += ' ' + format_me(self.muv)

48 if self.magn is not None:

49 string += ' ' + format_me(self.magn)

50 string += '\n'

52 return string

54 @classmethod

55 def fromstring(cls, string):

56 """Initialize yourself from a string"""

57 l = string.split()

58 energy = float(l.pop(0))

59 index = int(l.pop(0))

60 mur = np.array([float(l.pop(0)) for i in range(3)])

61 try:

62 muv = np.array([float(l.pop(0)) for i in range(3)])

63 except IndexError:

64 muv = None

65 try:

66 magn = np.array([float(l.pop(0)) for i in range(3)])

67 except IndexError:

68 magn = None

70 return cls(energy, index, mur, muv, magn)

72 def get_dipole_me(self, form='r'):

73 """Return the excitations dipole matrix element

74 including the occupation factor sqrt(fij)"""

75 if form == 'r': # length form

76 me = - self.mur

77 elif form == 'v': # velocity form

78 me = - self.muv

79 else:

80 raise RuntimeError('Unknown form >' + form + '<')

81 return np.sqrt(self.fij) * me

83 def get_dipole_tensor(self, form='r'):

84 """Return the oscillator strength tensor"""

85 me = self.get_dipole_me(form)

86 return 2 * np.outer(me, me.conj()) * self.energy / Hartree

88 def get_oscillator_strength(self, form='r'):

89 """Return the excitations dipole oscillator strength."""

90 me2_c = self.get_dipole_tensor(form).diagonal().real

91 return np.array([np.sum(me2_c) / 3.] + me2_c.tolist())

94class ExcitationList(list):

95 """Base class for excitions from the ground state"""

97 def __init__(self):

98 # initialise empty list

99 super().__init__()

101 # set default energy scale to get eV

102 self.energy_to_eV_scale = 1.

105def polarizability(exlist, omega, form='v',

106 tensor=False, index=0):

107 """Evaluate the photon energy dependent polarizability

108 from the sum over states

110 Parameters

111 ----------

112 exlist: ExcitationList

113 omega:

114 Photon energy (eV)

115 form: {'v', 'r'}

116 Form of the dipole matrix element, default 'v'

117 index: {0, 1, 2, 3}

118 0: averaged, 1,2,3:alpha_xx, alpha_yy, alpha_zz, default 0

119 tensor: boolean

120 if True returns alpha_ij, i,j=x,y,z

121 index is ignored, default False

123 Returns

124 -------

125 alpha:

126 Unit (e^2 Angstrom^2 / eV).

127 Multiply with Bohr * Ha to get (Angstrom^3)

128 shape = (omega.shape,) if tensor == False

129 shape = (omega.shape, 3, 3) else

130 """

131 omega = np.asarray(omega)

132 om2 = 1. * omega**2

133 esc = exlist.energy_to_eV_scale

135 if tensor:

136 if not np.isscalar(om2):

137 om2 = om2[:, None, None]

138 alpha = np.zeros(omega.shape + (3, 3),

139 dtype=om2.dtype)

140 for ex in exlist:

141 alpha += ex.get_dipole_tensor(form=form) / (

142 (ex.energy * esc)**2 - om2)

143 else:

144 alpha = np.zeros_like(om2)

145 for ex in exlist:

146 alpha += ex.get_oscillator_strength(form=form)[index] / (

147 (ex.energy * esc)**2 - om2)

149 return alpha * Bohr**2 * Hartree