Hide keyboard shortcuts

Hot-keys on this page

r m x p   toggle line displays

j k   next/prev highlighted chunk

0   (zero) top of page

1   (one) first highlighted chunk

1from itertools import count 

2import numpy as np 

3 

4from ase import Atoms 

5from ase.units import invcm, Ha 

6from ase.data import atomic_masses 

7from ase.calculators.calculator import all_changes 

8from ase.calculators.morse import MorsePotential 

9from ase.calculators.excitation_list import Excitation, ExcitationList 

10 

11"""The H2 molecule represented by Morse-Potentials for 

12gound and first 3 excited singlet states B + C(doubly degenerate)""" 

13 

14npa = np.array 

15# data from: 

16# https://webbook.nist.gov/cgi/cbook.cgi?ID=C1333740&Mask=1000#Diatomic 

17# X B C C 

18Re = npa([0.74144, 1.2928, 1.0327, 1.0327]) # eq. bond length 

19ome = npa([4401.21, 1358.09, 2443.77, 2443.77]) # vibrational frequency 

20# electronic transition energy 

21Etrans = npa([0, 91700.0, 100089.9, 100089.9]) * invcm 

22 

23# dissociation energy 

24# GS: https://aip.scitation.org/doi/10.1063/1.3120443 

25De = np.ones(4) * 36118.069 * invcm 

26# B, C separated energy E(1s) - E(2p) 

27De[1:] += Ha / 2 - Ha / 8 

28De -= Etrans 

29 

30# Morse parameter 

31m = atomic_masses[1] * 0.5 # reduced mass 

32# XXX find scaling factor 

33rho0 = Re * ome * invcm * np.sqrt(m / 2 / De) * 4401.21 / 284.55677429605862 

34 

35 

36def H2Morse(state=0): 

37 """Return H2 as a Morse-Potential with calculator attached.""" 

38 atoms = Atoms('H2', positions=np.zeros((2, 3))) 

39 atoms[1].position[2] = Re[state] 

40 atoms.calc = H2MorseCalculator(state=state) 

41 atoms.get_potential_energy() 

42 return atoms 

43 

44 

45class H2MorseCalculator(MorsePotential): 

46 """H2 ground or excited state as Morse potential""" 

47 _count = count(0) 

48 

49 def __init__(self, restart=None, state=0, rng=np.random): 

50 self.rng = rng 

51 MorsePotential.__init__(self, 

52 restart=restart, 

53 epsilon=De[state], 

54 r0=Re[state], rho0=rho0[state]) 

55 

56 def calculate(self, atoms=None, properties=['energy'], 

57 system_changes=all_changes): 

58 if atoms is not None: 

59 assert len(atoms) == 2 

60 MorsePotential.calculate(self, atoms, properties, system_changes) 

61 

62 # determine 'wave functions' including 

63 # Berry phase (arbitrary sign) and 

64 # random orientation of wave functions perpendicular 

65 # to the molecular axis 

66 

67 # molecular axis 

68 vr = atoms[1].position - atoms[0].position 

69 r = np.linalg.norm(vr) 

70 hr = vr / r 

71 # perpendicular axes 

72 vrand = self.rng.rand(3) 

73 hx = np.cross(hr, vrand) 

74 hx /= np.linalg.norm(hx) 

75 hy = np.cross(hr, hx) 

76 hy /= np.linalg.norm(hy) 

77 wfs = [1, hr, hx, hy] 

78 # Berry phase 

79 berry = (-1)**self.rng.randint(0, 2, 4) 

80 self.wfs = [wf * b for wf, b in zip(wfs, berry)] 

81 

82 def read(self, filename): 

83 ms = self 

84 with open(filename) as fd: 

85 ms.wfs = [int(fd.readline().split()[0])] 

86 for i in range(1, 4): 

87 ms.wfs.append( 

88 np.array([float(x) 

89 for x in fd.readline().split()[:4]])) 

90 ms.filename = filename 

91 return ms 

92 

93 def write(self, filename, option=None): 

94 """write calculated state to a file""" 

95 with open(filename, 'w') as fd: 

96 fd.write('{}\n'.format(self.wfs[0])) 

97 for wf in self.wfs[1:]: 

98 fd.write('{0:g} {1:g} {2:g}\n'.format(*wf)) 

99 

100 def overlap(self, other): 

101 ov = np.zeros((4, 4)) 

102 ov[0, 0] = self.wfs[0] * other.wfs[0] 

103 wfs = np.array(self.wfs[1:]) 

104 owfs = np.array(other.wfs[1:]) 

105 ov[1:, 1:] = np.dot(wfs, owfs.T) 

106 return ov 

107 

108 

109class H2MorseExcitedStatesCalculator(): 

110 """First singlet excited states of H2 from Morse potentials""" 

111 def __init__(self, nstates=3): 

112 """ 

113 Parameters 

114 ---------- 

115 nstates: int 

116 Numer of states to calculate 0 < nstates < 4, default 3 

117 """ 

118 assert nstates > 0 and nstates < 4 

119 self.nstates = nstates 

120 

121 def calculate(self, atoms): 

122 """Calculate excitation spectrum 

123 

124 Parameters 

125 ---------- 

126 atoms: Ase atoms object 

127 """ 

128 # central me value and rise, unit Bohr 

129 # from DOI: 10.1021/acs.jctc.9b00584 

130 mc = [0, 0.8, 0.7, 0.7] 

131 mr = [0, 1.0, 0.5, 0.5] 

132 

133 cgs = atoms.calc 

134 r = atoms.get_distance(0, 1) 

135 E0 = cgs.get_potential_energy(atoms) 

136 

137 exl = H2MorseExcitedStates() 

138 for i in range(1, self.nstates + 1): 

139 hvec = cgs.wfs[0] * cgs.wfs[i] 

140 energy = Ha * (0.5 - 1. / 8) - E0 

141 calc = H2MorseCalculator(state=i) 

142 calc.calculate(atoms) 

143 energy += calc.get_potential_energy() 

144 

145 mur = hvec * (mc[i] + (r - Re[0]) * mr[i]) 

146 muv = mur 

147 

148 exl.append(H2Excitation(energy, i, mur, muv)) 

149 return exl 

150 

151 

152class H2MorseExcitedStates(ExcitationList): 

153 """First singlet excited states of H2""" 

154 def __init__(self, nstates=3): 

155 """ 

156 Parameters 

157 ---------- 

158 nstates: int, 1 <= nstates <= 3 

159 Number of excited states to consider, default 3 

160 """ 

161 self.nstates = nstates 

162 super().__init__() 

163 

164 def overlap(self, ov_nn, other): 

165 return (ov_nn[1:len(self) + 1, 1:len(self) + 1] * 

166 ov_nn[0, 0]) 

167 

168 @classmethod 

169 def read(cls, filename, nstates=3): 

170 """Read myself from a file""" 

171 exl = cls(nstates) 

172 with open(filename, 'r') as fd: 

173 exl.filename = filename 

174 n = int(fd.readline().split()[0]) 

175 for i in range(min(n, exl.nstates)): 

176 exl.append(H2Excitation.fromstring(fd.readline())) 

177 return exl 

178 

179 def write(self, fname): 

180 with open(fname, 'w') as fd: 

181 fd.write('{0}\n'.format(len(self))) 

182 for ex in self: 

183 fd.write(ex.outstring()) 

184 

185 

186class H2Excitation(Excitation): 

187 def __eq__(self, other): 

188 """Considered to be equal when their indices are equal.""" 

189 return self.index == other.index 

190 

191 def __hash__(self): 

192 """Hash similar to __eq__""" 

193 if not hasattr(self, 'hash'): 

194 self.hash = hash(self.index) 

195 return self.hash 

196 

197 

198class H2MorseExcitedStatesAndCalculator( 

199 H2MorseExcitedStatesCalculator, H2MorseExcitedStates): 

200 """Traditional joined object for backward compatibility only""" 

201 def __init__(self, calculator, nstates=3): 

202 if isinstance(calculator, str): 

203 exlist = H2MorseExcitedStates.read(calculator, nstates) 

204 else: 

205 atoms = calculator.atoms 

206 atoms.calc = calculator 

207 excalc = H2MorseExcitedStatesCalculator(nstates) 

208 exlist = excalc.calculate(atoms) 

209 H2MorseExcitedStates.__init__(self, nstates=nstates) 

210 for ex in exlist: 

211 self.append(ex)