Coverage for /builds/kinetik161/ase/ase/calculators/h2morse.py: 100.00%

1from itertools import count 

2 

3import numpy as np 

4 

5from ase import Atoms 

6from ase.calculators.calculator import all_changes 

7from ase.calculators.excitation_list import Excitation, ExcitationList 

8from ase.calculators.morse import MorsePotential 

9from ase.data import atomic_masses 

10from ase.units import Ha, invcm 

11 

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

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

14 

15npa = np.array 

16# data from: 

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

18# X B C C 

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

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

21# electronic transition energy 

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

23 

24# dissociation energy 

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

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

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

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

29De -= Etrans 

30 

31# Morse parameter 

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

33# XXX find scaling factor 

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

35 

36 

37def H2Morse(state=0): 

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

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

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

41 atoms.calc = H2MorseCalculator(state=state) 

42 atoms.get_potential_energy() 

43 return atoms 

44 

45 

46class H2MorseCalculator(MorsePotential): 

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

48 _count = count(0) 

49 

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

51 self.rng = rng 

52 MorsePotential.__init__(self, 

53 restart=restart, 

54 epsilon=De[state], 

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

56 

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

58 system_changes=all_changes): 

59 if atoms is not None: 

60 assert len(atoms) == 2 

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

62 

63 # determine 'wave functions' including 

64 # Berry phase (arbitrary sign) and 

65 # random orientation of wave functions perpendicular 

66 # to the molecular axis 

67 

68 # molecular axis 

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

70 r = np.linalg.norm(vr) 

71 hr = vr / r 

72 # perpendicular axes 

73 vrand = self.rng.random(3) 

74 hx = np.cross(hr, vrand) 

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

76 hy = np.cross(hr, hx) 

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

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

79 # Berry phase 

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

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

82 

83 def read(self, filename): 

84 ms = self 

85 with open(filename) as fd: 

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

87 for i in range(1, 4): 

88 ms.wfs.append( 

89 np.array([float(x) 

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

91 ms.filename = filename 

92 return ms 

93 

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

95 """write calculated state to a file""" 

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

97 fd.write(f'{self.wfs[0]}\n') 

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

99 fd.write('{:g} {:g} {:g}\n'.format(*wf)) 

100 

101 def overlap(self, other): 

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

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

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

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

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

107 return ov 

108 

109 

110class H2MorseExcitedStatesCalculator(): 

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

112 

113 def __init__(self, nstates=3): 

114 """ 

115 Parameters 

116 ---------- 

117 nstates: int 

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

119 """ 

120 assert nstates > 0 and nstates < 4 

121 self.nstates = nstates 

122 

123 def calculate(self, atoms): 

124 """Calculate excitation spectrum 

125 

126 Parameters 

127 ---------- 

128 atoms: Ase atoms object 

129 """ 

130 # central me value and rise, unit Bohr 

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

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

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

134 

135 cgs = atoms.calc 

136 r = atoms.get_distance(0, 1) 

137 E0 = cgs.get_potential_energy(atoms) 

138 

139 exl = H2MorseExcitedStates() 

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

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

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

143 calc = H2MorseCalculator(state=i) 

144 calc.calculate(atoms) 

145 energy += calc.get_potential_energy() 

146 

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

148 muv = mur 

149 

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

151 return exl 

152 

153 

154class H2MorseExcitedStates(ExcitationList): 

155 """First singlet excited states of H2""" 

156 

157 def __init__(self, nstates=3): 

158 """ 

159 Parameters 

160 ---------- 

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

162 Number of excited states to consider, default 3 

163 """ 

164 self.nstates = nstates 

165 super().__init__() 

166 

167 def overlap(self, ov_nn, other): 

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

169 ov_nn[0, 0]) 

170 

171 @classmethod 

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

173 """Read myself from a file""" 

174 exl = cls(nstates) 

175 with open(filename) as fd: 

176 exl.filename = filename 

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

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

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

180 return exl 

181 

182 def write(self, fname): 

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

184 fd.write(f'{len(self)}\n') 

185 for ex in self: 

186 fd.write(ex.outstring()) 

187 

188 

189class H2Excitation(Excitation): 

190 def __eq__(self, other): 

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

192 return self.index == other.index 

193 

194 def __hash__(self): 

195 """Hash similar to __eq__""" 

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

197 self.hash = hash(self.index) 

198 return self.hash 

199 

200 

201class H2MorseExcitedStatesAndCalculator( 

202 H2MorseExcitedStatesCalculator, H2MorseExcitedStates): 

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

204 

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

206 if isinstance(calculator, str): 

207 exlist = H2MorseExcitedStates.read(calculator, nstates) 

208 else: 

209 atoms = calculator.atoms 

210 atoms.calc = calculator 

211 excalc = H2MorseExcitedStatesCalculator(nstates) 

212 exlist = excalc.calculate(atoms) 

213 H2MorseExcitedStates.__init__(self, nstates=nstates) 

214 for ex in exlist: 

215 self.append(ex)