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

1# fmt: off 

2 

3from itertools import count 

4 

5import numpy as np 

6 

7from ase import Atoms 

8from ase.calculators.calculator import all_changes 

9from ase.calculators.excitation_list import Excitation, ExcitationList 

10from ase.calculators.morse import MorsePotential 

11from ase.data import atomic_masses 

12from ase.units import Ha, invcm 

13 

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

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

16 

17npa = np.array 

18# data from: 

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

20# X B C C 

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

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

23# electronic transition energy 

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

25 

26# dissociation energy 

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

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

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

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

31De -= Etrans 

32 

33# Morse parameter 

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

35# XXX find scaling factor 

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

37 

38 

39def H2Morse(state=0): 

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

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

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

43 atoms.calc = H2MorseCalculator(state=state) 

44 atoms.get_potential_energy() 

45 return atoms 

46 

47 

48class H2MorseCalculator(MorsePotential): 

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

50 _count = count(0) 

51 

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

53 self.rng = rng 

54 MorsePotential.__init__(self, 

55 restart=restart, 

56 epsilon=De[state], 

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

58 

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

60 system_changes=all_changes): 

61 if atoms is not None: 

62 assert len(atoms) == 2 

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

64 

65 # determine 'wave functions' including 

66 # Berry phase (arbitrary sign) and 

67 # random orientation of wave functions perpendicular 

68 # to the molecular axis 

69 

70 # molecular axis 

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

72 r = np.linalg.norm(vr) 

73 hr = vr / r 

74 # perpendicular axes 

75 vrand = self.rng.random(3) 

76 hx = np.cross(hr, vrand) 

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

78 hy = np.cross(hr, hx) 

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

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

81 # Berry phase 

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

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

84 

85 def read(self, filename): 

86 ms = self 

87 with open(filename) as fd: 

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

89 for _ in range(1, 4): 

90 ms.wfs.append( 

91 np.array([float(x) 

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

93 ms.filename = filename 

94 return ms 

95 

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

97 """write calculated state to a file""" 

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

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

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

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

102 

103 def overlap(self, other): 

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

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

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

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

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

109 return ov 

110 

111 

112class H2MorseExcitedStatesCalculator(): 

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

114 

115 def __init__(self, nstates=3): 

116 """ 

117 Parameters 

118 ---------- 

119 nstates: int 

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

121 """ 

122 assert nstates > 0 and nstates < 4 

123 self.nstates = nstates 

124 

125 def calculate(self, atoms): 

126 """Calculate excitation spectrum 

127 

128 Parameters 

129 ---------- 

130 atoms: Ase atoms object 

131 """ 

132 # central me value and rise, unit Bohr 

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

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

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

136 

137 cgs = atoms.calc 

138 r = atoms.get_distance(0, 1) 

139 E0 = cgs.get_potential_energy(atoms) 

140 

141 exl = H2MorseExcitedStates() 

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

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

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

145 calc = H2MorseCalculator(state=i) 

146 calc.calculate(atoms) 

147 energy += calc.get_potential_energy() 

148 

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

150 muv = mur 

151 

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

153 return exl 

154 

155 

156class H2MorseExcitedStates(ExcitationList): 

157 """First singlet excited states of H2""" 

158 

159 def __init__(self, nstates=3): 

160 """ 

161 Parameters 

162 ---------- 

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

164 Number of excited states to consider, default 3 

165 """ 

166 self.nstates = nstates 

167 super().__init__() 

168 

169 def overlap(self, ov_nn, other): 

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

171 ov_nn[0, 0]) 

172 

173 @classmethod 

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

175 """Read myself from a file""" 

176 exl = cls(nstates) 

177 with open(filename) as fd: 

178 exl.filename = filename 

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

180 for _ in range(min(n, exl.nstates)): 

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

182 return exl 

183 

184 def write(self, fname): 

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

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

187 for ex in self: 

188 fd.write(ex.outstring()) 

189 

190 

191class H2Excitation(Excitation): 

192 def __eq__(self, other): 

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

194 return self.index == other.index 

195 

196 def __hash__(self): 

197 """Hash similar to __eq__""" 

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

199 self.hash = hash(self.index) 

200 return self.hash 

201 

202 

203class H2MorseExcitedStatesAndCalculator( 

204 H2MorseExcitedStatesCalculator, H2MorseExcitedStates): 

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

206 

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

208 if isinstance(calculator, str): 

209 exlist = H2MorseExcitedStates.read(calculator, nstates) 

210 else: 

211 atoms = calculator.atoms 

212 atoms.calc = calculator 

213 excalc = H2MorseExcitedStatesCalculator(nstates) 

214 exlist = excalc.calculate(atoms) 

215 H2MorseExcitedStates.__init__(self, nstates=nstates) 

216 for ex in exlist: 

217 self.append(ex)