1# -*- coding: utf-8 -*- 

2# 

3# Copyright (c) 2017, the cclib development team 

4# 

5# This file is part of cclib (http://cclib.github.io) and is distributed under 

6# the terms of the BSD 3-Clause License. 

7 

8"""Parser for Psi3 output files.""" 

9 

10import numpy 

11 

12from cclib.parser import logfileparser 

13from cclib.parser import utils 

14 

15 

16class Psi3(logfileparser.Logfile): 

17 """A Psi3 log file.""" 

18 

19 def __init__(self, *args, **kwargs): 

20 

21 # Call the __init__ method of the superclass 

22 super(Psi3, self).__init__(logname="Psi3", *args, **kwargs) 

23 

24 def __str__(self): 

25 """Return a string representation of the object.""" 

26 return "Psi3 log file %s" % (self.filename) 

27 

28 def __repr__(self): 

29 """Return a representation of the object.""" 

30 return 'Psi3("%s")' % (self.filename) 

31 

32 def normalisesym(self, label): 

33 """Psi3 does not require normalizing symmetry labels.""" 

34 return label 

35 

36 def extract(self, inputfile, line): 

37 """Extract information from the file object inputfile.""" 

38 

39 if "Version" in line: 

40 self.metadata["package_version"] = ''.join(line.split()[1:]).lower() 

41 

42 # Psi3 prints the coordinates in several configurations, and we will parse the 

43 # the canonical coordinates system in Angstroms as the first coordinate set, 

44 # although it is actually somewhere later in the input, after basis set, etc. 

45 # We can also get or verify the number of atoms and atomic numbers from this block. 

46 if line.strip() == "-Geometry in the canonical coordinate system (Angstrom):": 

47 

48 self.skip_lines(inputfile, ['header', 'd']) 

49 

50 coords = [] 

51 numbers = [] 

52 line = next(inputfile) 

53 while line.strip(): 

54 

55 tokens = line.split() 

56 

57 element = tokens[0] 

58 numbers.append(self.table.number[element]) 

59 

60 x = float(tokens[1]) 

61 y = float(tokens[2]) 

62 z = float(tokens[3]) 

63 coords.append([x, y, z]) 

64 

65 line = next(inputfile) 

66 

67 self.set_attribute('natom', len(coords)) 

68 self.set_attribute('atomnos', numbers) 

69 

70 if not hasattr(self, 'atomcoords'): 

71 self.atomcoords = [] 

72 self.atomcoords.append(coords) 

73 

74 if line.strip() == '-SYMMETRY INFORMATION:': 

75 line = next(inputfile) 

76 while line.strip(): 

77 if "Number of atoms" in line: 

78 self.set_attribute('natom', int(line.split()[-1])) 

79 line = next(inputfile) 

80 if line.strip() == "-BASIS SET INFORMATION:": 

81 line = next(inputfile) 

82 while line.strip(): 

83 if "Number of SO" in line: 

84 self.set_attribute('nbasis', int(line.split()[-1])) 

85 line = next(inputfile) 

86 

87 # In Psi3, the section with the contraction scheme can be used to infer atombasis. 

88 if line.strip() == "-Contraction Scheme:": 

89 

90 self.skip_lines(inputfile, ['header', 'd']) 

91 

92 indices = [] 

93 line = next(inputfile) 

94 while line.strip(): 

95 shells = line.split('//')[-1] 

96 expression = shells.strip().replace(' ', '+') 

97 expression = expression.replace('s', '*1') 

98 expression = expression.replace('p', '*3') 

99 expression = expression.replace('d', '*6') 

100 nfuncs = eval(expression) 

101 if len(indices) == 0: 

102 indices.append(range(nfuncs)) 

103 else: 

104 start = indices[-1][-1] + 1 

105 indices.append(range(start, start+nfuncs)) 

106 line = next(inputfile) 

107 

108 self.set_attribute('atombasis', indices) 

109 

110 if line.strip() == "CINTS: An integrals program written in C": 

111 

112 self.skip_lines(inputfile, ['authors', 'd', 'b', 'b']) 

113 

114 line = next(inputfile) 

115 assert line.strip() == "-OPTIONS:" 

116 while line.strip(): 

117 line = next(inputfile) 

118 

119 line = next(inputfile) 

120 assert line.strip() == "-CALCULATION CONSTANTS:" 

121 while line.strip(): 

122 if "Number of atoms" in line: 

123 natom = int(line.split()[-1]) 

124 self.set_attribute('natom', natom) 

125 if "Number of symmetry orbitals" in line: 

126 nbasis = int(line.split()[-1]) 

127 self.set_attribute('nbasis', nbasis) 

128 line = next(inputfile) 

129 

130 if line.strip() == "CSCF3.0: An SCF program written in C": 

131 

132 self.skip_lines(inputfile, ['b', 'authors', 'b', 'd', 'b', 

133 'mult', 'mult_comment', 'b']) 

134 

135 line = next(inputfile) 

136 while line.strip(): 

137 if line.split()[0] == "multiplicity": 

138 mult = int(line.split()[-1]) 

139 self.set_attribute('mult', mult) 

140 if line.split()[0] == "charge": 

141 charge = int(line.split()[-1]) 

142 self.set_attribute('charge', charge) 

143 if line.split()[0] == "convergence": 

144 conv = float(line.split()[-1]) 

145 if line.split()[0] == "reference": 

146 self.reference = line.split()[-1] 

147 line = next(inputfile) 

148 

149 if not hasattr(self, 'scftargets'): 

150 self.scftargets = [] 

151 self.scftargets.append([conv]) 

152 

153 # ==> Iterations <== 

154 

155 # Psi3 converges just the density elements, although it reports in the iterations 

156 # changes in the energy as well as the DIIS error. 

157 psi3_iterations_header = "iter total energy delta E delta P diiser" 

158 if line.strip() == psi3_iterations_header: 

159 

160 if not hasattr(self, 'scfvalues'): 

161 self.scfvalues = [] 

162 self.scfvalues.append([]) 

163 

164 line = next(inputfile) 

165 while line.strip(): 

166 ddensity = float(line.split()[-2]) 

167 self.scfvalues[-1].append([ddensity]) 

168 line = next(inputfile) 

169 

170 # This section, from which we parse molecular orbital symmetries and 

171 # orbital energies, is quite similar for both Psi3 and Psi4, and in fact 

172 # the format for orbtials is the same, although the headers and spacers 

173 # are a bit different. Let's try to get both parsed with one code block. 

174 # 

175 # Here is how the block looks like for Psi4: 

176 # 

177 # Orbital Energies (a.u.) 

178 # ----------------------- 

179 # 

180 # Doubly Occupied: 

181 # 

182 # 1Bu -11.040586 1Ag -11.040524 2Bu -11.031589 

183 # 2Ag -11.031589 3Bu -11.028950 3Ag -11.028820 

184 # (...) 

185 # 15Ag -0.415620 1Bg -0.376962 2Au -0.315126 

186 # 2Bg -0.278361 3Bg -0.222189 

187 # 

188 # Virtual: 

189 # 

190 # 3Au 0.198995 4Au 0.268517 4Bg 0.308826 

191 # 5Au 0.397078 5Bg 0.521759 16Ag 0.565017 

192 # (...) 

193 # 24Ag 0.990287 24Bu 1.027266 25Ag 1.107702 

194 # 25Bu 1.124938 

195 # 

196 # The case is different in the trigger string. 

197 if "orbital energies (a.u.)" in line.lower(): 

198 

199 self.moenergies = [[]] 

200 self.mosyms = [[]] 

201 

202 self.skip_line(inputfile, 'blank') 

203 

204 occupied = next(inputfile) 

205 if self.reference[0:2] == 'RO' or self.reference[0:1] == 'R': 

206 assert 'doubly occupied' in occupied.lower() 

207 elif self.reference[0:1] == 'U': 

208 assert 'alpha occupied' in occupied.lower() 

209 

210 # Parse the occupied MO symmetries and energies. 

211 self._parse_mosyms_moenergies(inputfile, 0) 

212 

213 # The last orbital energy here represents the HOMO. 

214 self.homos = [len(self.moenergies[0])-1] 

215 # For a restricted open-shell calculation, this is the 

216 # beta HOMO, and we assume the singly-occupied orbitals 

217 # are all alpha, which are handled next. 

218 if self.reference[0:2] == 'RO': 

219 self.homos.append(self.homos[0]) 

220 

221 self.skip_line(inputfile, 'blank') 

222 

223 unoccupied = next(inputfile) 

224 if self.reference[0:2] == 'RO': 

225 assert unoccupied.strip() == 'Singly Occupied:' 

226 elif self.reference[0:1] == 'R': 

227 assert unoccupied.strip() == 'Unoccupied orbitals' 

228 elif self.reference[0:1] == 'U': 

229 assert unoccupied.strip() == 'Alpha Virtual:' 

230 

231 # Parse the unoccupied MO symmetries and energies. 

232 self._parse_mosyms_moenergies(inputfile, 0) 

233 

234 # Here is where we handle the Beta or Singly occupied orbitals. 

235 if self.reference[0:1] == 'U': 

236 self.mosyms.append([]) 

237 self.moenergies.append([]) 

238 line = next(inputfile) 

239 assert line.strip() == 'Beta Occupied:' 

240 self.skip_line(inputfile, 'blank') 

241 self._parse_mosyms_moenergies(inputfile, 1) 

242 self.homos.append(len(self.moenergies[1])-1) 

243 line = next(inputfile) 

244 assert line.strip() == 'Beta Virtual:' 

245 self.skip_line(inputfile, 'blank') 

246 self._parse_mosyms_moenergies(inputfile, 1) 

247 elif self.reference[0:2] == 'RO': 

248 line = next(inputfile) 

249 assert line.strip() == 'Virtual:' 

250 self.skip_line(inputfile, 'blank') 

251 self._parse_mosyms_moenergies(inputfile, 0) 

252 

253 # Both Psi3 and Psi4 print the final SCF energy right after 

254 # the orbital energies, but the label is different. Psi4 also 

255 # does DFT, and the label is also different in that case. 

256 if "* SCF total energy" in line: 

257 e = float(line.split()[-1]) 

258 if not hasattr(self, 'scfenergies'): 

259 self.scfenergies = [] 

260 self.scfenergies.append(utils.convertor(e, 'hartree', 'eV')) 

261 

262 # We can also get some higher moments in Psi3, although here the dipole is not printed 

263 # separately and the order is not lexicographical. However, the numbers seem 

264 # kind of strange -- the quadrupole seems to be traceless, although I'm not sure 

265 # whether the standard transformation has been used. So, until we know what kind 

266 # of moment these are and how to make them raw again, we will only parse the dipole. 

267 # 

268 # -------------------------------------------------------------- 

269 # *** Electric multipole moments *** 

270 # -------------------------------------------------------------- 

271 # 

272 # CAUTION : The system has non-vanishing dipole moment, therefore 

273 # quadrupole and higher moments depend on the reference point. 

274 # 

275 # -Coordinates of the reference point (a.u.) : 

276 # x y z 

277 # -------------------- -------------------- -------------------- 

278 # 0.0000000000 0.0000000000 0.0000000000 

279 # 

280 # -Electric dipole moment (expectation values) : 

281 # 

282 # mu(X) = -0.00000 D = -1.26132433e-43 C*m = -0.00000000 a.u. 

283 # mu(Y) = 0.00000 D = 3.97987832e-44 C*m = 0.00000000 a.u. 

284 # mu(Z) = 0.00000 D = 0.00000000e+00 C*m = 0.00000000 a.u. 

285 # |mu| = 0.00000 D = 1.32262368e-43 C*m = 0.00000000 a.u. 

286 # 

287 # -Components of electric quadrupole moment (expectation values) (a.u.) : 

288 # 

289 # Q(XX) = 10.62340220 Q(YY) = 1.11816843 Q(ZZ) = -11.74157063 

290 # Q(XY) = 3.64633112 Q(XZ) = 0.00000000 Q(YZ) = 0.00000000 

291 # 

292 if line.strip() == "*** Electric multipole moments ***": 

293 

294 self.skip_lines(inputfile, ['d', 'b', 'caution1', 'caution2', 'b']) 

295 

296 coordinates = next(inputfile) 

297 assert coordinates.split()[-2] == "(a.u.)" 

298 self.skip_lines(inputfile, ['xyz', 'd']) 

299 line = next(inputfile) 

300 self.origin = numpy.array([float(x) for x in line.split()]) 

301 self.origin = utils.convertor(self.origin, 'bohr', 'Angstrom') 

302 

303 self.skip_line(inputfile, "blank") 

304 line = next(inputfile) 

305 assert "Electric dipole moment" in line 

306 self.skip_line(inputfile, "blank") 

307 

308 # Make sure to use the column that has the value in Debyes. 

309 dipole = [] 

310 for i in range(3): 

311 line = next(inputfile) 

312 dipole.append(float(line.split()[2])) 

313 

314 if not hasattr(self, 'moments'): 

315 self.moments = [self.origin, dipole] 

316 else: 

317 assert self.moments[1] == dipole 

318 

319 def _parse_mosyms_moenergies(self, inputfile, spinidx): 

320 """Parse molecular orbital symmetries and energies from the 

321 'Post-Iterations' section. 

322 """ 

323 line = next(inputfile) 

324 while line.strip(): 

325 for i in range(len(line.split()) // 2): 

326 self.mosyms[spinidx].append(line.split()[i*2][-2:]) 

327 moenergy = utils.convertor(float(line.split()[i*2+1]), "hartree", "eV") 

328 self.moenergies[spinidx].append(moenergy) 

329 line = next(inputfile) 

330 return