Coverage for cclib/parser/gaussianparser.py : 96%
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
1# -*- coding: utf-8 -*-
2#
3# Copyright (c) 2020, 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.
8"""Parser for Gaussian output files"""
11import re
13import numpy
15from cclib.parser import data
16from cclib.parser import logfileparser
17from cclib.parser import utils
20class Gaussian(logfileparser.Logfile):
21 """A Gaussian 98/03 log file."""
23 def __init__(self, *args, **kwargs):
25 # Call the __init__ method of the superclass
26 super(Gaussian, self).__init__(logname="Gaussian", *args, **kwargs)
28 def __str__(self):
29 """Return a string representation of the object."""
30 return "Gaussian log file %s" % (self.filename)
32 def __repr__(self):
33 """Return a representation of the object."""
34 return 'Gaussian("%s")' % (self.filename)
36 def normalisesym(self, label):
37 """Use standard symmetry labels instead of Gaussian labels.
39 To normalise:
40 (1) If label is one of [SG, PI, PHI, DLTA], replace by [sigma, pi, phi, delta]
41 (2) replace any G or U by their lowercase equivalent
42 """
43 # note: DLT must come after DLTA
44 greek = [('SG', 'sigma'), ('PI', 'pi'), ('PHI', 'phi'),
45 ('DLTA', 'delta'), ('DLT', 'delta')]
46 for k, v in greek:
47 if label.startswith(k):
48 tmp = label[len(k):]
49 label = v
50 if tmp:
51 label = v + "." + tmp
53 ans = label.replace("U", "u").replace("G", "g")
54 return ans
56 # Use to map from the usual year suffixed to full years so package
57 # versions can be sorted properly after parsing with
58 # `packaging.parse.version`.
59 YEAR_SUFFIXES_TO_YEARS = {
60 '70': '1970',
61 '76': '1976',
62 '80': '1980',
63 '82': '1982',
64 '86': '1986',
65 '88': '1988',
66 '90': '1990',
67 '92': '1992',
68 '94': '1994',
69 '98': '1998',
70 '03': '2003',
71 '09': '2009',
72 '16': '2016',
73 }
75 def before_parsing(self):
77 # Extract only well-formed numbers in scientific notation.
78 self.re_scinot = re.compile(r'(\w*)=\s*(-?\d\.\d{2}D[+-]\d{2})')
79 # Extract only well-formed numbers in traditional
80 # floating-point format.
81 self.re_float = re.compile(r'(\w*-?\w*)=\s*(-?\d+\.\d{10,})')
83 # Flag for identifying Coupled Cluster runs.
84 self.coupledcluster = False
86 # Fragment number for counterpoise or fragment guess calculations
87 # (normally zero).
88 self.counterpoise = 0
90 # Flag for identifying ONIOM calculations.
91 self.oniom = False
93 # Flag for identifying BOMD calculations.
94 # These calculations have a back-integration algorithm so that not all
95 # geometries should be kept.
96 # We also add a "time" attribute to the parser.
97 self.BOMD = False
99 # Do we have high-precision polarizabilities printed from a
100 # dedicated `polar` job? If so, avoid duplicate parsing.
101 self.hp_polarizabilities = False
103 def after_parsing(self):
104 # atomcoords are parsed as a list of lists but it should be an array
105 if hasattr(self, "atomcoords"):
106 self.atomcoords = numpy.array(self.atomcoords)
108 # Correct the percent values in the etsecs in the case of
109 # a restricted calculation. The following has the
110 # effect of including each transition twice.
111 if hasattr(self, "etsecs") and len(self.homos) == 1:
112 new_etsecs = [[(x[0], x[1], x[2] * numpy.sqrt(2)) for x in etsec]
113 for etsec in self.etsecs]
114 self.etsecs = new_etsecs
116 if hasattr(self, "scanenergies"):
117 self.scancoords = []
118 if hasattr(self, 'optstatus') and hasattr(self, 'atomcoords'):
119 converged_indexes = [x for x, y in enumerate(self.optstatus) if y & data.ccData.OPT_DONE > 0]
120 self.scancoords = self.atomcoords[converged_indexes,:,:]
121 elif hasattr(self, 'atomcoords'):
122 self.scancoords = self.atomcoords
124 if (hasattr(self, 'enthalpy') and hasattr(self, 'temperature')
125 and hasattr(self, 'freeenergy')):
126 self.set_attribute('entropy', (self.enthalpy - self.freeenergy) / self.temperature)
128 # This bit is needed in order to trim coordinates that are printed a second time
129 # at the end of geometry optimizations. Note that we need to do this for both atomcoords
130 # and inputcoords. The reason is that normally a standard orientation is printed and that
131 # is what we parse into atomcoords, but inputcoords stores the input (unmodified) coordinates
132 # and that is copied over to atomcoords if no standard oritentation was printed, which happens
133 # for example for jobs with no symmetry. This last step, however, is now generic for all parsers.
134 # Perhaps then this part should also be generic code...
135 # Regression that tests this: Gaussian03/cyclopropenyl.rhf.g03.cut.log
136 if hasattr(self, 'optstatus') and len(self.optstatus) > 0:
137 last_point = len(self.optstatus) - 1
138 if hasattr(self, 'atomcoords'):
139 self.atomcoords = self.atomcoords[:last_point + 1]
140 if hasattr(self, 'inputcoords'):
141 self.inputcoords = self.inputcoords[:last_point + 1]
143 # If we parsed high-precision vibrational displacements, overwrite
144 # lower-precision displacements in self.vibdisps
145 if hasattr(self, 'vibdispshp'):
146 self.vibdisps = self.vibdispshp
147 del self.vibdispshp
148 if hasattr(self, 'time'):
149 self.time = [self.time[i] for i in sorted(self.time.keys())]
150 if hasattr(self, 'energies_BOMD'):
151 self.set_attribute('scfenergies',
152 [self.energies_BOMD[i] for i in sorted(self.energies_BOMD.keys())])
153 if hasattr(self, 'atomcoords_BOMD'):
154 self.atomcoords= \
155 [self.atomcoords_BOMD[i] for i in sorted(self.atomcoords_BOMD.keys())]
157 # Gaussian prints 'forces' in input orientation unlike other values such as 'moments' or 'vibdisp'.
158 # Therefore, we convert 'grads' to the values in standard orientation with rotation matrix.
159 if hasattr(self, 'grads') and hasattr(self, 'inputcoords') and hasattr(self, 'atomcoords'):
160 grads_std = []
161 for grad, inputcoord, atomcoord in zip(self.grads, self.inputcoords, self.atomcoords):
162 rotation = utils.get_rotation(numpy.array(inputcoord), numpy.array(atomcoord))
163 grads_std.append(rotation.apply(grad))
164 self.set_attribute('grads', numpy.array(grads_std))
166 def extract(self, inputfile, line):
167 """Extract information from the file object inputfile."""
169 # Extract the version number: "Gaussian 09, Revision D.01"
170 # becomes "09revisionD.01".
171 if line.strip() == "Cite this work as:":
172 tokens = next(inputfile).split()
173 self.metadata["legacy_package_version"] = ''.join([
174 tokens[1][:-1],
175 'revision',
176 tokens[-1][:-1],
177 ])
179 # Extract the version number: "Gaussian 98: x86-Linux-G98RevA.11.3
180 # 5-Feb-2002" becomes "1998+A.11.3", and "Gaussian 16:
181 # ES64L-G16RevA.03 25-Dec-2016" becomes "2016+A.03".
182 if "Gaussian, Inc.," in line:
183 self.skip_lines(inputfile, ["b", "s"])
184 _, _, platform_full_version, compile_date = next(inputfile).split()
185 run_date = next(inputfile).strip()
186 platform_full_version_tokens = platform_full_version.split("-")
187 full_version = platform_full_version_tokens[-1]
188 platform = "-".join(platform_full_version_tokens[:-1])
189 year_suffix = full_version[1:3]
190 revision = full_version[6:]
191 self.metadata["package_version"] = "{}+{}".format(
192 self.YEAR_SUFFIXES_TO_YEARS[year_suffix], revision
193 )
194 self.metadata["platform"] = platform
196 if line.strip().startswith("Link1: Proceeding to internal job step number"):
197 self.new_internal_job()
199 # This block contains some general information as well as coordinates,
200 # which could be parsed in the future:
201 #
202 # Symbolic Z-matrix:
203 # Charge = 0 Multiplicity = 1
204 # C 0.73465 0. 0.
205 # C 1.93465 0. 0.
206 # C
207 # ...
208 #
209 # It also lists fragments, if there are any, which is potentially valuable:
210 #
211 # Symbolic Z-matrix:
212 # Charge = 0 Multiplicity = 1 in supermolecule
213 # Charge = 0 Multiplicity = 1 in fragment 1.
214 # Charge = 0 Multiplicity = 1 in fragment 2.
215 # B(Fragment=1) 0.06457 -0.0279 0.01364
216 # H(Fragment=1) 0.03117 -0.02317 1.21604
217 # ...
218 #
219 # Note, however, that currently we only parse information for the whole system
220 # or supermolecule as Gaussian calls it.
221 if line.strip() == "Symbolic Z-matrix:":
223 self.updateprogress(inputfile, "Symbolic Z-matrix", self.fupdate)
225 line = inputfile.next()
226 while line.split()[0] == 'Charge':
228 # For the supermolecule, we can parse the charge and multicplicity.
229 regex = r".*=(.*)Mul.*=\s*-?(\d+).*"
230 match = re.match(regex, line)
231 assert match, "Something unusual about the line: '%s'" % line
233 self.set_attribute('charge', int(match.groups()[0]))
234 self.set_attribute('mult', int(match.groups()[1]))
236 if line.split()[-2] == "fragment":
237 self.nfragments = int(line.split()[-1].strip('.'))
239 if line.strip()[-13:] == "model system.":
240 self.nmodels = getattr(self, 'nmodels', 0) + 1
242 line = inputfile.next()
244 # The remaining part will allow us to get the atom count.
245 # When coordinates are given, there is a blank line at the end, but if
246 # there is a Z-matrix here, there will also be variables and we need to
247 # stop at those to get the right atom count.
248 # Also, in older versions there is bo blank line (G98 regressions),
249 # so we need to watch out for leaving the link.
250 natom = 0
251 while line.split() and not "Variables" in line and not "Leave Link" in line:
252 natom += 1
253 line = inputfile.next()
254 self.set_attribute('natom', natom)
256 # Continuing from above, there is not always a symbolic matrix, for example
257 # if the Z-matrix was in the input file. In such cases, try to match the
258 # line and get at the charge and multiplicity.
259 #
260 # Charge = 0 Multiplicity = 1 in supermolecule
261 # Charge = 0 Multiplicity = 1 in fragment 1.
262 # Charge = 0 Multiplicity = 1 in fragment 2.
263 if line[1:7] == 'Charge' and line.find("Multiplicity") >= 0:
265 self.updateprogress(inputfile, "Charge and Multiplicity", self.fupdate)
267 if line.split()[-1] == "supermolecule" or not "fragment" in line:
269 regex = r".*=(.*)Mul.*=\s*-?(\d+).*"
270 match = re.match(regex, line)
271 assert match, "Something unusual about the line: '%s'" % line
273 self.set_attribute('charge', int(match.groups()[0]))
274 self.set_attribute('mult', int(match.groups()[1]))
276 if line.split()[-2] == "fragment":
277 self.nfragments = int(line.split()[-1].strip('.'))
279 # Number of atoms is also explicitely printed after the above.
280 if line[1:8] == "NAtoms=":
282 self.updateprogress(inputfile, "Attributes", self.fupdate)
284 natom = int(re.search(r'NAtoms=\s*(\d+)', line).group(1))
285 self.set_attribute('natom', natom)
287 # Necessary for `if line.strip().split()[0:3] == ["Atom", "AN", "X"]:` block
288 if not hasattr(self, 'nqmf'):
289 match = re.search('NQMF=\s*(\d+)', line)
290 if match is not None:
291 nqmf = int(match.group(1))
292 if nqmf > 0:
293 self.set_attribute('nqmf', nqmf)
295 # Basis set name
296 if line[1:15] == "Standard basis":
297 self.metadata["basis_set"] = line.split()[2]
299 # Dipole moment
300 # e.g. from G09
301 # Dipole moment (field-independent basis, Debye):
302 # X= 0.0000 Y= 0.0000 Z= 0.0930
303 # e.g. from G03
304 # X= 0.0000 Y= 0.0000 Z= -1.6735 Tot= 1.6735
305 # need the "field independent" part - ONIOM and other calc use diff formats
306 if line[1:39] == "Dipole moment (field-independent basis":
308 self.updateprogress(inputfile, "Dipole and Higher Moments", self.fupdate)
310 self.reference = [0.0, 0.0, 0.0]
311 self.moments = [self.reference]
313 tokens = inputfile.next().split()
314 # split - dipole would need to be *huge* to fail a split
315 # and G03 and G09 use different spacing
316 if len(tokens) >= 6:
317 dipole = (float(tokens[1]), float(tokens[3]), float(tokens[5]))
319 if not hasattr(self, 'moments'):
320 self.moments = [self.reference, dipole]
321 else:
322 self.moments.append(dipole)
324 if line[1:43] == "Quadrupole moment (field-independent basis":
325 # e.g. (g09)
326 # Quadrupole moment (field-independent basis, Debye-Ang):
327 # XX= -6.1213 YY= -4.2950 ZZ= -5.4175
328 # XY= 0.0000 XZ= 0.0000 YZ= 0.0000
329 # or from g03
330 # XX= -6.1213 YY= -4.2950 ZZ= -5.4175
331 quadrupole = {}
332 for j in range(2): # two rows
333 line = inputfile.next()
334 if line[22] == '=': # g03 file
335 for i in (1, 18, 35):
336 quadrupole[line[i:i+4]] = float(line[i+5:i+16])
337 else:
338 for i in (1, 27, 53):
339 quadrupole[line[i:i+4]] = float(line[i+5:i+25])
341 lex = sorted(quadrupole.keys())
342 quadrupole = [quadrupole[key] for key in lex]
344 if not hasattr(self, 'moments') or len(self.moments) < 2:
345 self.logger.warning("Found quadrupole moments but no previous dipole")
346 self.reference = [0.0, 0.0, 0.0]
347 self.moments = [self.reference, None, quadrupole]
348 else:
349 if len(self.moments) == 2:
350 self.moments.append(quadrupole)
351 else:
352 assert self.moments[2] == quadrupole
354 if line[1:41] == "Octapole moment (field-independent basis":
355 # e.g.
356 # Octapole moment (field-independent basis, Debye-Ang**2):
357 # XXX= 0.0000 YYY= 0.0000 ZZZ= -0.1457 XYY= 0.0000
358 # XXY= 0.0000 XXZ= 0.0136 XZZ= 0.0000 YZZ= 0.0000
359 # YYZ= -0.5848 XYZ= 0.0000
360 octapole = {}
361 for j in range(2): # two rows
362 line = inputfile.next()
363 if line[22] == '=': # g03 file
364 for i in (1, 18, 35, 52):
365 octapole[line[i:i+4]] = float(line[i+5:i+16])
366 else:
367 for i in (1, 27, 53, 79):
368 octapole[line[i:i+4]] = float(line[i+5:i+25])
370 # last line only 2 moments
371 line = inputfile.next()
372 if line[22] == '=': # g03 file
373 for i in (1, 18):
374 octapole[line[i:i+4]] = float(line[i+5:i+16])
375 else:
376 for i in (1, 27):
377 octapole[line[i:i+4]] = float(line[i+5:i+25])
379 lex = sorted(octapole.keys())
380 octapole = [octapole[key] for key in lex]
382 if not hasattr(self, 'moments') or len(self.moments) < 3:
383 self.logger.warning("Found octapole moments but no previous dipole or quadrupole")
384 self.reference = [0.0, 0.0, 0.0]
385 self.moments = [self.reference, None, None, octapole]
386 else:
387 if len(self.moments) == 3:
388 self.moments.append(octapole)
389 else:
390 assert self.moments[3] == octapole
392 if line[1:20] == "Hexadecapole moment":
393 # e.g.
394 # Hexadecapole moment (field-independent basis, Debye-Ang**3):
395 # XXXX= -3.2614 YYYY= -6.8264 ZZZZ= -4.9965 XXXY= 0.0000
396 # XXXZ= 0.0000 YYYX= 0.0000 YYYZ= 0.0000 ZZZX= 0.0000
397 # ZZZY= 0.0000 XXYY= -1.8585 XXZZ= -1.4123 YYZZ= -1.7504
398 # XXYZ= 0.0000 YYXZ= 0.0000 ZZXY= 0.0000
399 hexadecapole = {}
400 # read three lines worth of 4 moments per line
401 for j in range(3):
402 line = inputfile.next()
403 if line[22] == '=': # g03 file
404 for i in (1, 18, 35, 52):
405 hexadecapole[line[i:i+4]] = float(line[i+5:i+16])
406 else:
407 for i in (1, 27, 53, 79):
408 hexadecapole[line[i:i+4]] = float(line[i+5:i+25])
410 # last line only 3 moments
411 line = inputfile.next()
412 if line[22] == '=': # g03 file
413 for i in (1, 18, 35):
414 hexadecapole[line[i:i+4]] = float(line[i+5:i+16])
415 else:
416 for i in (1, 27, 53):
417 hexadecapole[line[i:i+4]] = float(line[i+5:i+25])
419 lex = sorted(hexadecapole.keys())
420 hexadecapole = [hexadecapole[key] for key in lex]
422 if not hasattr(self, 'moments') or len(self.moments) < 4:
423 self.reference = [0.0, 0.0, 0.0]
424 self.moments = [self.reference, None, None, None, hexadecapole]
425 else:
426 if len(self.moments) == 4:
427 self.append_attribute("moments", hexadecapole)
428 else:
429 try:
430 numpy.testing.assert_equal(self.moments[4], hexadecapole)
431 except AssertionError:
432 self.logger.warning("Attribute hexadecapole changed value (%s -> %s)" % (self.moments[4], hexadecapole))
433 self.append_attribute("moments", hexadecapole)
435 # Catch message about completed optimization.
436 if line[1:23] == "Optimization completed":
438 if not hasattr(self, 'optdone'):
439 self.optdone = []
440 self.optdone.append(len(self.geovalues) - 1)
442 assert hasattr(self, "optstatus") and len(self.optstatus) > 0
443 self.optstatus[-1] += data.ccData.OPT_DONE
445 # Catch message about stopped optimization (not converged).
446 if line[1:21] == "Optimization stopped":
448 if not hasattr(self, "optdone"):
449 self.optdone = []
451 assert hasattr(self, "optstatus") and len(self.optstatus) > 0
452 self.optstatus[-1] += data.ccData.OPT_UNCONVERGED
454 # Extract the atomic numbers and coordinates from the input orientation,
455 # in the event the standard orientation isn't available.
456 # Don't extract from Input or Z-matrix orientation in a BOMD run, as only
457 # the final geometry should be kept but extract inputatoms.
458 # We also use "inputcoords" to convert "grads" from input orientation
459 # to standard orientation
460 if line.find("Input orientation") > -1 or line.find("Z-Matrix orientation") > -1:
462 # If this is a counterpoise calculation, this output means that
463 # the supermolecule is now being considered, so we can set:
464 self.counterpoise = 0
466 self.updateprogress(inputfile, "Attributes", self.cupdate)
468 if not self.BOMD and not hasattr(self, "inputcoords"):
469 self.inputcoords = []
470 self.inputatoms = []
472 self.skip_lines(inputfile, ['d', 'cols', 'cols', 'd'])
474 atomcoords = []
475 line = next(inputfile)
476 while list(set(line.strip())) != ["-"]:
477 broken = line.split()
478 atomno = int(broken[1])
479 # Atom with atomno -1 only appears on "Z-Matrix orientation", and excluded on
480 # "Input orientation" or "Standard orientation".
481 # We remove this line to keep shape consistency of "atomcoords" and "inputcoords",
482 # so that we can convert "grads" from input orientation to standard orientaion
483 # with rotation matrix calculated from "atomcoords" and "inputcoords"
484 if atomno != -1:
485 self.inputatoms.append(atomno)
486 atomcoords.append(list(map(float, broken[3:6])))
487 line = next(inputfile)
489 if not self.BOMD: self.inputcoords.append(atomcoords)
491 self.set_attribute('atomnos', numpy.array(self.inputatoms))
492 self.set_attribute('natom', len(self.inputatoms))
494 if self.BOMD and line.startswith(' Summary information for step'):
496 # We keep time and energies_BOMD and coordinates in a dictionary
497 # because steps can be recalculated, and we need to overwite the
498 # previous data
500 broken = line.split()
501 step = int(broken[-1])
502 line = next(inputfile)
503 broken = line.split()
504 if not hasattr(self, "time"):
505 self.set_attribute('time', {step:float(broken[-1])})
506 else:
507 self.time[step] = float(broken[-1])
509 line = next(inputfile)
510 broken = line.split(';')[1].split()
511 ene = utils.convertor(utils.float(broken[-1]), "hartree", "eV")
512 if not hasattr(self, "energies_BOMD"):
513 self.set_attribute('energies_BOMD', {step:ene})
514 else:
515 self.energies_BOMD[step] = ene
517 self.updateprogress(inputfile, "Attributes", self.cupdate)
519 if not hasattr(self, "atomcoords_BOMD"):
520 self.atomcoords_BOMD = {}
521 #self.inputatoms = []
523 self.skip_lines(inputfile, ['EKin', 'Angular', 'JX', 'Total', 'Total', 'Cartesian'])
525 atomcoords = []
526 line = next(inputfile)
527 while not "MW cartesian" in line:
528 broken = line.split()
529 atomcoords.append(list(map(utils.float, (broken[3], broken[5], broken[7]))))
530 # self.inputatoms.append(int(broken[1]))
531 line = next(inputfile)
533 self.atomcoords_BOMD[step] = atomcoords
535 #self.set_attribute('atomnos', self.inputatoms)
536 #self.set_attribute('natom', len(self.inputatoms))
539 # Extract the atomic masses.
540 # Typical section:
541 # Isotopes and Nuclear Properties:
542 #(Nuclear quadrupole moments (NQMom) in fm**2, nuclear magnetic moments (NMagM)
543 # in nuclear magnetons)
544 #
545 # Atom 1 2 3 4 5 6 7 8 9 10
546 # IAtWgt= 12 12 12 12 12 1 1 1 12 12
547 # AtmWgt= 12.0000000 12.0000000 12.0000000 12.0000000 12.0000000 1.0078250 1.0078250 1.0078250 12.0000000 12.0000000
548 # NucSpn= 0 0 0 0 0 1 1 1 0 0
549 # AtZEff= -3.6000000 -3.6000000 -3.6000000 -3.6000000 -3.6000000 -1.0000000 -1.0000000 -1.0000000 -3.6000000 -3.6000000
550 # NQMom= 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
551 # NMagM= 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 2.7928460 2.7928460 2.7928460 0.0000000 0.0000000
552 # ... with blank lines dividing blocks of ten, and Leave Link 101 at the end.
553 # This is generally parsed before coordinates, so atomnos is not defined.
554 # Note that in Gaussian03 the comments are not there yet and the labels are different.
555 if line.strip() == "Isotopes and Nuclear Properties:":
557 if not hasattr(self, "atommasses"):
558 self.atommasses = []
560 line = next(inputfile)
561 while line[1:16] != "Leave Link 101":
562 if line[1:8] == "AtmWgt=":
563 self.atommasses.extend(list(map(float, line.split()[1:])))
564 line = next(inputfile)
566 # Extract the atomic numbers and coordinates of the atoms.
567 if line.strip() == "Standard orientation:":
569 self.updateprogress(inputfile, "Attributes", self.cupdate)
571 # If this is a counterpoise calculation, this output means that
572 # the supermolecule is now being considered, so we can set:
573 self.counterpoise = 0
575 if not hasattr(self, "atomcoords"):
576 self.atomcoords = []
578 self.skip_lines(inputfile, ['d', 'cols', 'cols', 'd'])
580 atomnos = []
581 atomcoords = []
582 line = next(inputfile)
583 while list(set(line.strip())) != ["-"]:
584 broken = line.split()
585 atomnos.append(int(broken[1]))
586 atomcoords.append(list(map(float, broken[-3:])))
587 line = next(inputfile)
588 self.atomcoords.append(atomcoords)
590 self.set_attribute('natom', len(atomnos))
591 self.set_attribute('atomnos', atomnos)
593 # This is a bit of a hack for regression Gaussian09/BH3_fragment_guess.pop_minimal.log
594 # to skip output for all fragments, assuming the supermolecule is always printed first.
595 # Eventually we want to make this more general, or even better parse the output for
596 # all fragment, but that will happen in a newer version of cclib.
597 if line[1:16] == "Fragment guess:" and getattr(self, 'nfragments', 0) > 1:
598 if not "full" in line:
599 inputfile.seek(0, 2)
601 # Another hack for regression Gaussian03/ortho_prod_prod_freq.log, which is an ONIOM job.
602 # Basically for now we stop parsing after the output for the real system, because
603 # currently we don't support changes in system size or fragments in cclib. When we do,
604 # we will want to parse the model systems, too, and that is what nmodels could track.
605 if "ONIOM: generating point" in line and line.strip()[-13:] == 'model system.' and getattr(self, 'nmodels', 0) > 0:
606 inputfile.seek(0, 2)
608 # With the gfinput keyword, the atomic basis set functios are:
609 #
610 # AO basis set in the form of general basis input (Overlap normalization):
611 # 1 0
612 # S 3 1.00 0.000000000000
613 # 0.7161683735D+02 0.1543289673D+00
614 # 0.1304509632D+02 0.5353281423D+00
615 # 0.3530512160D+01 0.4446345422D+00
616 # SP 3 1.00 0.000000000000
617 # 0.2941249355D+01 -0.9996722919D-01 0.1559162750D+00
618 # 0.6834830964D+00 0.3995128261D+00 0.6076837186D+00
619 # 0.2222899159D+00 0.7001154689D+00 0.3919573931D+00
620 # ****
621 # 2 0
622 # S 3 1.00 0.000000000000
623 # 0.7161683735D+02 0.1543289673D+00
624 # ...
625 #
626 # The same is also printed when the gfprint keyword is used, but the
627 # interstitial lines differ and there are no stars between atoms:
628 #
629 # AO basis set (Overlap normalization):
630 # Atom C1 Shell 1 S 3 bf 1 - 1 0.509245180608 -2.664678875191 0.000000000000
631 # 0.7161683735D+02 0.1543289673D+00
632 # 0.1304509632D+02 0.5353281423D+00
633 # 0.3530512160D+01 0.4446345422D+00
634 # Atom C1 Shell 2 SP 3 bf 2 - 5 0.509245180608 -2.664678875191 0.000000000000
635 # 0.2941249355D+01 -0.9996722919D-01 0.1559162750D+00
636 # ...
638 #ONIOM calculations result basis sets reported for atoms that are not in order of atom number which breaks this code (line 390 relies on atoms coming in order)
639 if line[1:13] == "AO basis set" and not self.oniom:
641 self.gbasis = []
643 # For counterpoise fragment calcualtions, skip these lines.
644 if self.counterpoise != 0:
645 return
647 atom_line = inputfile.next()
648 self.gfprint = atom_line.split()[0] == "Atom"
649 self.gfinput = not self.gfprint
651 # Note how the shell information is on a separate line for gfinput,
652 # whereas for gfprint it is on the same line as atom information.
653 if self.gfinput:
654 shell_line = inputfile.next()
656 shell = []
657 while len(self.gbasis) < self.natom:
659 if self.gfprint:
660 cols = atom_line.split()
661 subshells = cols[4]
662 ngauss = int(cols[5])
663 else:
664 cols = shell_line.split()
665 subshells = cols[0]
666 ngauss = int(cols[1])
668 parameters = []
669 for ig in range(ngauss):
670 line = inputfile.next()
671 parameters.append(list(map(utils.float, line.split())))
672 for iss, ss in enumerate(subshells):
673 contractions = []
674 for param in parameters:
675 exponent = param[0]
676 coefficient = param[iss+1]
677 contractions.append((exponent, coefficient))
678 subshell = (ss, contractions)
679 shell.append(subshell)
681 if self.gfprint:
682 line = inputfile.next()
683 if line.split()[0] == "Atom":
684 atomnum = int(re.sub(r"\D", "", line.split()[1]))
685 if atomnum == len(self.gbasis) + 2:
686 self.gbasis.append(shell)
687 shell = []
688 atom_line = line
689 else:
690 self.gbasis.append(shell)
691 else:
692 line = inputfile.next()
693 if line.strip() == "****":
694 self.gbasis.append(shell)
695 shell = []
696 atom_line = inputfile.next()
697 shell_line = inputfile.next()
698 else:
699 shell_line = line
701 if "Dispersion energy=" in line:
702 dispersion = utils.convertor(float(line.split()[-2]), "hartree", "eV")
703 self.append_attribute("dispersionenergies", dispersion)
705 # Find the targets for SCF convergence (QM calcs).
706 # Not for BOMD as targets are not available in the summary
707 if not self.BOMD and line[1:44] == 'Requested convergence on RMS density matrix':
709 if not hasattr(self, "scftargets"):
710 self.scftargets = []
711 # The following can happen with ONIOM which are mixed SCF
712 # and semi-empirical
713 if type(self.scftargets) == type(numpy.array([])):
714 self.scftargets = []
716 scftargets = []
717 # The RMS density matrix.
718 scftargets.append(utils.float(line.split('=')[1].split()[0]))
719 line = next(inputfile)
720 # The MAX density matrix.
721 scftargets.append(utils.float(line.strip().split('=')[1][:-1]))
722 line = next(inputfile)
723 # For G03, there's also the energy (not for G98).
724 if line[1:10] == "Requested":
725 scftargets.append(utils.float(line.strip().split('=')[1][:-1]))
727 self.scftargets.append(scftargets)
729 # Extract SCF convergence information (QM calcs).
730 if line[1:10] == 'Cycle 1':
732 if not hasattr(self, "scfvalues"):
733 self.scfvalues = []
735 scfvalues = []
736 line = next(inputfile)
737 while line.find("SCF Done") == -1:
739 self.updateprogress(inputfile, "QM convergence", self.fupdate)
741 if line.find(' E=') == 0:
742 self.logger.debug(line)
744 # RMSDP=3.74D-06 MaxDP=7.27D-05 DE=-1.73D-07 OVMax= 3.67D-05
745 # or
746 # RMSDP=1.13D-05 MaxDP=1.08D-04 OVMax= 1.66D-04
747 if line.find(" RMSDP") == 0:
749 # Fields of interest:
750 # RMSDP
751 # MaxDP
752 # (DE) -> Only add the energy if it's a target criteria
754 matches = self.re_scinot.findall(line)
755 matches = {
756 match[0]: utils.float(match[1])
757 for match in matches
758 }
759 scfvalues_step = [
760 matches.get('RMSDP', numpy.nan),
761 matches.get('MaxDP', numpy.nan)
762 ]
763 if hasattr(self, "scftargets") and len(self.scftargets[0]) == 3:
764 scfvalues_step.append(matches.get('DE', numpy.nan))
765 scfvalues.append(scfvalues_step)
767 try:
768 line = next(inputfile)
769 # May be interupted by EOF.
770 except StopIteration:
771 self.logger.warning('File terminated before end of last SCF!')
772 break
774 self.scfvalues.append(scfvalues)
776 # Extract SCF convergence information (AM1, INDO and other semi-empirical calcs).
777 # The output (for AM1) looks like this:
778 # Ext34=T Pulay=F Camp-King=F BShift= 0.00D+00
779 # It= 1 PL= 0.103D+01 DiagD=T ESCF= 31.564733 Diff= 0.272D+02 RMSDP= 0.152D+00.
780 # It= 2 PL= 0.114D+00 DiagD=T ESCF= 7.265370 Diff=-0.243D+02 RMSDP= 0.589D-02.
781 # ...
782 # It= 11 PL= 0.184D-04 DiagD=F ESCF= 4.687669 Diff= 0.260D-05 RMSDP= 0.134D-05.
783 # It= 12 PL= 0.105D-04 DiagD=F ESCF= 4.687669 Diff=-0.686D-07 RMSDP= 0.215D-05.
784 # 4-point extrapolation.
785 # It= 13 PL= 0.110D-05 DiagD=F ESCF= 4.687669 Diff=-0.111D-06 RMSDP= 0.653D-07.
786 # Energy= 0.172272018655 NIter= 14.
787 if line[1:4] == 'It=':
789 scftargets = numpy.array([1E-7], "d") # This is the target value for the rms
790 scfvalues = [[]]
792 while line.find(" Energy") == -1:
794 self.updateprogress(inputfile, "AM1 Convergence")
796 if line[1:4] == "It=":
797 parts = line.strip().split()
798 scfvalues[0].append(utils.float(parts[-1][:-1]))
800 line = next(inputfile)
802 # If an AM1 or INDO guess is used (Guess=INDO in the input, for example),
803 # this will be printed after a single iteration, so that is the line
804 # that should trigger a break from this loop. At least that's what we see
805 # for regression Gaussian/Gaussian09/guessIndo_modified_ALT.out
806 if line[:14] == " Initial guess":
807 break
809 # Attach the attributes to the object Only after the energy is found .
810 if line.find(" Energy") == 0:
811 self.scftargets = scftargets
812 self.scfvalues = scfvalues
814 # Note: this needs to follow the section where 'SCF Done' is used
815 # to terminate a loop when extracting SCF convergence information.
816 if not self.BOMD and line[1:9] == 'SCF Done':
818 t1 = line.split()[2]
819 if t1 == 'E(RHF)':
820 self.metadata["methods"].append("HF")
821 else:
822 self.metadata["methods"].append("DFT")
823 self.metadata["functional"] = t1[t1.index("(") + 2:t1.rindex(")")]
825 if not hasattr(self, "scfenergies"):
826 self.scfenergies = []
828 self.scfenergies.append(utils.convertor(utils.float(line.split()[4]), "hartree", "eV"))
829 # gmagoon 5/27/09: added scfenergies reading for PM3 case
830 # Example line: " Energy= -0.077520562724 NIter= 14."
831 # See regression Gaussian03/QVGXLLKOCUKJST-UHFFFAOYAJmult3Fixed.out
832 if line[1:8] == 'Energy=':
833 if not hasattr(self, "scfenergies"):
834 self.scfenergies = []
835 self.scfenergies.append(utils.convertor(utils.float(line.split()[1]), "hartree", "eV"))
837 # Total energies after Moller-Plesset corrections.
838 # Second order correction is always first, so its first occurance
839 # triggers creation of mpenergies (list of lists of energies).
840 # Further MP2 corrections are appended as found.
841 #
842 # Example MP2 output line:
843 # E2 = -0.9505918144D+00 EUMP2 = -0.28670924198852D+03
844 # Warning! this output line is subtly different for MP3/4/5 runs
845 if "EUMP2" in line[27:34]:
846 self.metadata["methods"].append("MP2")
848 if not hasattr(self, "mpenergies"):
849 self.mpenergies = []
850 self.mpenergies.append([])
851 mp2energy = utils.float(line.split("=")[2])
852 self.mpenergies[-1].append(utils.convertor(mp2energy, "hartree", "eV"))
854 # Example MP3 output line:
855 # E3= -0.10518801D-01 EUMP3= -0.75012800924D+02
856 if line[34:39] == "EUMP3":
857 self.metadata["methods"].append("MP3")
859 mp3energy = utils.float(line.split("=")[2])
860 self.mpenergies[-1].append(utils.convertor(mp3energy, "hartree", "eV"))
862 # Example MP4 output lines:
863 # E4(DQ)= -0.31002157D-02 UMP4(DQ)= -0.75015901139D+02
864 # E4(SDQ)= -0.32127241D-02 UMP4(SDQ)= -0.75016013648D+02
865 # E4(SDTQ)= -0.32671209D-02 UMP4(SDTQ)= -0.75016068045D+02
866 # Energy for most substitutions is used only (SDTQ by default)
867 if line[34:42] == "UMP4(DQ)":
868 self.metadata["methods"].append("MP4")
870 mp4energy = utils.float(line.split("=")[2])
871 line = next(inputfile)
872 if line[34:43] == "UMP4(SDQ)":
873 mp4energy = utils.float(line.split("=")[2])
874 line = next(inputfile)
875 if line[34:44] == "UMP4(SDTQ)":
876 mp4energy = utils.float(line.split("=")[2])
877 self.mpenergies[-1].append(utils.convertor(mp4energy, "hartree", "eV"))
879 # Example MP5 output line:
880 # DEMP5 = -0.11048812312D-02 MP5 = -0.75017172926D+02
881 if line[29:32] == "MP5":
882 self.metadata["methods"].append("MP5")
883 mp5energy = utils.float(line.split("=")[2])
884 self.mpenergies[-1].append(utils.convertor(mp5energy, "hartree", "eV"))
886 # Total energies after Coupled Cluster corrections.
887 # Second order MBPT energies (MP2) are also calculated for these runs,
888 # but the output is the same as when parsing for mpenergies.
889 # Read the consecutive correlated energies
890 # but append only the last one to ccenergies.
891 # Only the highest level energy is appended - ex. CCSD(T), not CCSD.
892 if line[1:10] == "DE(Corr)=" and line[27:35] == "E(CORR)=":
893 self.metadata["methods"].append("CCSD")
894 self.ccenergy = utils.float(line.split()[3])
895 if line[1:14] == "T1 Diagnostic":
896 self.metadata["t1_diagnostic"] = utils.float(line.split()[-1])
897 if line[1:10] == "T5(CCSD)=":
898 line = next(inputfile)
899 if line[1:9] == "CCSD(T)=":
900 self.metadata["methods"].append("CCSD-T")
901 self.ccenergy = utils.float(line.split()[1])
902 if line[12:53] == "Population analysis using the SCF density":
903 if hasattr(self, "ccenergy"):
904 if not hasattr(self, "ccenergies"):
905 self.ccenergies = []
906 self.ccenergies.append(utils.convertor(self.ccenergy, "hartree", "eV"))
907 del self.ccenergy
908 # Find step number for current optimization/IRC
909 # Matches "Step number 123", "Pt XX Step number 123" and "PtXXX Step number 123"
910 if " Step number" in line:
911 step = int(line.split()[line.split().index('Step') + 2])
912 if not hasattr(self, "optstatus"):
913 self.optstatus = []
914 self.optstatus.append(data.ccData.OPT_UNKNOWN)
915 if step == 1:
916 self.optstatus[-1] += data.ccData.OPT_NEW
918 # Geometry convergence information.
919 if line[49:59] == 'Converged?':
921 if not hasattr(self, "geotargets"):
922 self.geovalues = []
923 self.geotargets = numpy.array([0.0, 0.0, 0.0, 0.0], "d")
924 allconverged = True
925 newlist = [0]*4
926 for i in range(4):
927 line = next(inputfile)
928 self.logger.debug(line)
929 parts = line.split()
930 if "NO" in parts[-1]:
931 allconverged = False
932 try:
933 value = utils.float(parts[2])
934 except ValueError:
935 self.logger.error("Problem parsing the value for geometry optimisation: %s is not a number." % parts[2])
936 else:
937 newlist[i] = value
938 self.geotargets[i] = utils.float(parts[3])
939 # reset some parameters that are printed each iteration if the
940 # optimization has not yet converged. For example, etenergies
941 # (Issue #889) and similar properties are only reported for the
942 # final step of an optimization.
943 if not allconverged:
944 for reset_attr in ["etenergies", "etoscs", "etsyms", "etsecs", "etdips", "etveldips", "etmagdips"]:
945 if hasattr(self, reset_attr):
946 setattr(self, reset_attr, [])
948 self.geovalues.append(newlist)
950 # Gradients.
951 # Read in the cartesian energy gradients (forces) from a block like this:
952 # -------------------------------------------------------------------
953 # Center Atomic Forces (Hartrees/Bohr)
954 # Number Number X Y Z
955 # -------------------------------------------------------------------
956 # 1 1 -0.012534744 -0.021754635 -0.008346094
957 # 2 6 0.018984731 0.032948887 -0.038003451
958 # 3 1 -0.002133484 -0.006226040 0.023174772
959 # 4 1 -0.004316502 -0.004968213 0.023174772
960 # -2 -0.001830728 -0.000743108 -0.000196625
961 # ------------------------------------------------------------------
962 #
963 # The "-2" line is for a dummy atom
964 #
965 # Then optimization is done in internal coordinates, Gaussian also
966 # print the forces in internal coordinates, which can be produced from
967 # the above. This block looks like this:
968 # Variable Old X -DE/DX Delta X Delta X Delta X New X
969 # (Linear) (Quad) (Total)
970 # ch 2.05980 0.01260 0.00000 0.01134 0.01134 2.07114
971 # hch 1.75406 0.09547 0.00000 0.24861 0.24861 2.00267
972 # hchh 2.09614 0.01261 0.00000 0.16875 0.16875 2.26489
973 # Item Value Threshold Converged?
975 # We could get the gradients in BOMD, but it is more complex because
976 # they are not in the summary, and they are not as relevant as for
977 # an optimization
978 if not self.BOMD and line[37:43] == "Forces":
980 if not hasattr(self, "grads"):
981 self.grads = []
983 self.skip_lines(inputfile, ['header', 'd'])
985 forces = []
986 line = next(inputfile)
987 while list(set(line.strip())) != ['-']:
988 tmpforces = []
989 for N in range(3): # Fx, Fy, Fz
990 force = line[23+N*15:38+N*15]
991 if force.startswith("*"):
992 force = "NaN"
993 tmpforces.append(float(force))
994 forces.append(tmpforces)
995 line = next(inputfile)
996 self.grads.append(forces)
998 if "Number of optimizations in scan" in line:
999 self.scan_length = int(line.split()[-1])
1001 # All PES scans have a list of initial parameters from which we
1002 # can get the names and more.
1003 #
1004 # ----------------------------
1005 # ! Initial Parameters !
1006 # ! (Angstroms and Degrees) !
1007 # -------------------------- --------------------------
1008 # ! Name Definition Value Derivative Info. !
1009 # --------------------------------------------------------------------------------
1010 # ! R1 R(1,2) 1.4212 estimate D2E/DX2 !
1011 # ! R2 R(1,14) 1.4976 estimate D2E/DX2 !
1012 # ...
1013 if "Initial Parameters" in line:
1014 self.scannames_all = []
1015 self.scannames_scanned = []
1016 self.skip_lines(inputfile, ['units', 'd', 'header', 'd'])
1017 line = next(inputfile)
1018 while line.strip()[0] == '!':
1019 name = line.split()[1]
1020 definition = line.split()[2]
1021 self.scannames_all.append(name)
1022 if line.split()[4] == 'Scan':
1023 self.scannames_scanned.append(name)
1024 self.append_attribute('scannames', definition)
1025 line = next(inputfile)
1027 # Extract unrelaxed PES scan data, which looks something like:
1028 #
1029 # Summary of the potential surface scan:
1030 # N A SCF
1031 # ---- --------- -----------
1032 # 1 109.0000 -76.43373
1033 # 2 119.0000 -76.43011
1034 # 3 129.0000 -76.42311
1035 # 4 139.0000 -76.41398
1036 # 5 149.0000 -76.40420
1037 # 6 159.0000 -76.39541
1038 # 7 169.0000 -76.38916
1039 # 8 179.0000 -76.38664
1040 # 9 189.0000 -76.38833
1041 # 10 199.0000 -76.39391
1042 # 11 209.0000 -76.40231
1043 # ---- --------- -----------
1044 if "Summary of the potential surface scan:" in line:
1046 colmnames = next(inputfile)
1047 if not hasattr(self, "scannames"):
1048 self.set_attribute("scannames", colmnames.split()[1:-1])
1050 hyphens = next(inputfile)
1051 line = next(inputfile)
1052 scanparm = [[] for _ in range(len(self.scannames))]
1053 while line != hyphens:
1054 broken = line.split()
1055 self.append_attribute('scanenergies', (utils.convertor(float(broken[-1]), "hartree", "eV")))
1056 for idx,p in enumerate(broken[1:-1]):
1057 scanparm[idx].append(float(p))
1058 #self.append_attribute('scanparm', [float(p) for p in broken[1:-1]])
1059 line = next(inputfile)
1060 self.set_attribute('scanparm', scanparm)
1062 # Extract relaxed (optimized) PES scan data, for which the form
1063 # of the output is transposed:
1064 #
1065 # Summary of Optimized Potential Surface Scan (add -382.0 to energies):
1066 # 1 2 3 4 5
1067 # Eigenvalues -- -0.30827 -0.30695 -0.30265 -0.29955 -0.30260
1068 # R1 1.42115 1.42152 1.42162 1.42070 1.42071
1069 # R2 1.49761 1.49787 1.49855 1.49901 1.49858
1070 # R3 1.42245 1.42185 1.42062 1.42048 1.42147
1071 # R4 1.40217 1.40236 1.40306 1.40441 1.40412
1072 # ...
1073 if "Summary of Optimized Potential Surface Scan" in line:
1074 # The base energy separation is version dependent, and first
1075 # appears in Gaussian16.
1076 base_energy = 0.0
1077 if "add" in line and "to energies" in line:
1078 base_energy = float(line.split()[-3])
1080 scanenergies = []
1081 scanparm = [[] for _ in range(len(self.scannames))]
1082 while len(scanenergies) != self.scan_length:
1083 line = next(inputfile)
1084 indices = [int(i) for i in line.split()]
1085 widths = [10]*len(indices)
1086 splitter = utils.WidthSplitter(widths)
1088 line = next(inputfile)
1089 eigenvalues_in_line = line[21:].rstrip()
1090 assert len(eigenvalues_in_line) == sum(widths)
1091 cols = list(splitter.split(eigenvalues_in_line))
1092 try:
1093 eigenvalues = [float(e) for e in cols]
1094 eigenvalues = [base_energy + e for e in eigenvalues]
1095 except ValueError:
1096 eigenvalues = [numpy.nan for _ in cols]
1097 assert len(eigenvalues) == len(indices)
1098 eigenvalues = [utils.convertor(e, "hartree", "eV") for e in eigenvalues]
1099 scanenergies.extend(eigenvalues)
1101 for _, name in enumerate(self.scannames_all):
1102 line = next(inputfile)
1103 assert line.split()[0] == name
1104 if name in self.scannames_scanned:
1105 iname = self.scannames_scanned.index(name)
1106 params_in_line = line[21:].rstrip()
1107 assert len(params_in_line) == sum(widths)
1108 params = [float(v) for v in splitter.split(params_in_line)]
1109 scanparm[iname].extend(params)
1111 self.set_attribute('scanenergies', scanenergies)
1112 self.set_attribute('scanparm', scanparm)
1114 # Orbital symmetries.
1115 if line[1:20] == 'Orbital symmetries:' and not hasattr(self, "mosyms"):
1117 # For counterpoise fragments, skip these lines.
1118 if self.counterpoise != 0:
1119 return
1121 self.updateprogress(inputfile, "MO Symmetries", self.fupdate)
1123 self.mosyms = [[]]
1124 line = next(inputfile)
1125 unres = False
1126 if line.find("Alpha Orbitals") == 1:
1127 unres = True
1128 line = next(inputfile)
1129 i = 0
1130 while len(line) > 18 and line[17] == '(':
1131 if line.find('Virtual') >= 0:
1132 self.homos = [i - 1]
1133 parts = line[17:].split()
1134 for x in parts:
1135 self.mosyms[0].append(self.normalisesym(x.strip('()')))
1136 i += 1
1137 line = next(inputfile)
1138 if unres:
1139 line = next(inputfile)
1140 # Repeat with beta orbital information
1141 i = 0
1142 self.mosyms.append([])
1143 while len(line) > 18 and line[17] == '(':
1144 if line.find('Virtual') >= 0:
1145 # Here we consider beta
1146 # If there was also an alpha virtual orbital,
1147 # we will store two indices in the array
1148 # Otherwise there is no alpha virtual orbital,
1149 # only beta virtual orbitals, and we initialize
1150 # the array with one element. See the regression
1151 # QVGXLLKOCUKJST-UHFFFAOYAJmult3Fixed.out
1152 # donated by Gregory Magoon (gmagoon).
1153 if (hasattr(self, "homos")):
1154 # Extend the array to two elements
1155 # 'HOMO' indexes the HOMO in the arrays
1156 self.homos.append(i-1)
1157 else:
1158 # 'HOMO' indexes the HOMO in the arrays
1159 self.homos = [i - 1]
1160 parts = line[17:].split()
1161 for x in parts:
1162 self.mosyms[1].append(self.normalisesym(x.strip('()')))
1163 i += 1
1164 line = next(inputfile)
1166 # Some calculations won't explicitely print the number of basis sets used,
1167 # and will occasionally drop some without warning. We can infer the number,
1168 # however, from the MO symmetries printed here. Specifically, this fixes
1169 # regression Gaussian/Gaussian09/dvb_sp_terse.log (#23 on github).
1170 self.set_attribute('nmo', len(self.mosyms[-1]))
1172 # Alpha/Beta electron eigenvalues.
1173 if line[1:6] == "Alpha" and line.find("eigenvalues") >= 0:
1175 # For counterpoise fragments, skip these lines.
1176 if self.counterpoise != 0:
1177 return
1179 # For ONIOM calcs, ignore this section in order to bypass assertion failure.
1180 if self.oniom:
1181 return
1183 self.updateprogress(inputfile, "Eigenvalues", self.fupdate)
1184 self.moenergies = [[]]
1185 HOMO = -2
1187 while line.find('Alpha') == 1:
1188 if line.split()[1] == "virt." and HOMO == -2:
1190 # If there aren't any symmetries, this is a good way to find the HOMO.
1191 HOMO = len(self.moenergies[0])-1
1192 self.homos = [HOMO]
1194 # Convert to floats and append to moenergies, but sometimes Gaussian
1195 # doesn't print correctly so test for ValueError (bug 1756789).
1196 part = line[28:]
1197 i = 0
1198 while i*10+4 < len(part):
1199 s = part[i*10:(i+1)*10]
1200 try:
1201 x = utils.float(s)
1202 except ValueError:
1203 x = numpy.nan
1204 self.moenergies[0].append(utils.convertor(x, "hartree", "eV"))
1205 i += 1
1206 line = next(inputfile)
1208 # If, at this point, self.homos is unset, then there were not
1209 # any alpha virtual orbitals
1210 if not hasattr(self, "homos"):
1211 HOMO = len(self.moenergies[0])-1
1212 self.homos = [HOMO]
1214 if line.find('Beta') == 2:
1215 self.moenergies.append([])
1217 HOMO = -2
1218 while line.find('Beta') == 2:
1219 if line.split()[1] == "virt." and HOMO == -2:
1221 # If there aren't any symmetries, this is a good way to find the HOMO.
1222 # Also, check for consistency if homos was already parsed.
1223 HOMO = len(self.moenergies[1])-1
1224 self.homos.append(HOMO)
1226 part = line[28:]
1227 i = 0
1228 while i*10+4 < len(part):
1229 x = part[i*10:(i+1)*10]
1230 self.moenergies[1].append(utils.convertor(utils.float(x), "hartree", "eV"))
1231 i += 1
1232 line = next(inputfile)
1234 self.moenergies = [numpy.array(x, "d") for x in self.moenergies]
1236 # Start of the IR/Raman frequency section.
1237 # Caution is advised here, as additional frequency blocks
1238 # can be printed by Gaussian (with slightly different formats),
1239 # often doubling the information printed.
1240 # See, for a non-standard exmaple, regression Gaussian98/test_H2.log
1241 # If either the Gaussian freq=hpmodes keyword or IOP(7/33=1) is used,
1242 # an extra frequency block with higher-precision vibdisps is
1243 # printed before the normal frequency block.
1244 # Note that the code parses only the vibsyms and vibdisps
1245 # from the high-precision block, but parses vibsyms, vibfreqs, vibfconsts,
1246 # vibramans, vibrmasses and vibirs from the normal block. vibsyms parsed
1247 # from the high-precision block are discarded and replaced by those
1248 # from the normal block while the high-precision vibdisps, if present,
1249 # are used to overwrite default-precision vibdisps at the end of the parse.
1250 if line[1:14] == "Harmonic freq": # This matches in both freq block types
1252 self.updateprogress(inputfile, "Frequency Information", self.fupdate)
1254 # The whole block should not have any blank lines.
1255 while line.strip() != "":
1257 # The line with indices
1258 if line[1:15].strip() == "" and line[15:60].split()[0].isdigit():
1259 freqbase = int(line[15:60].split()[0])
1260 if freqbase == 1 and hasattr(self, 'vibsyms'):
1261 # we are coming accross duplicated information.
1262 # We might be be parsing a default-precision block having
1263 # already parsed (only) vibsyms and displacements from
1264 # the high-precision block, or might be encountering
1265 # a second low-precision block (see e.g. 25DMF_HRANH.log
1266 # regression).
1267 self.vibsyms = []
1268 if hasattr(self, "vibirs"):
1269 self.vibirs = []
1270 if hasattr(self, 'vibfreqs'):
1271 self.vibfreqs = []
1272 if hasattr(self, 'vibramans'):
1273 self.vibramans = []
1274 if hasattr(self, "vibrmasses"):
1275 self.vibrmasses = []
1276 if hasattr(self, "vibfconsts"):
1277 self.vibfconsts = []
1278 if hasattr(self, 'vibdisps'):
1279 self.vibdisps = []
1281 # Lines with symmetries and symm. indices begin with whitespace.
1282 if line[1:15].strip() == "" and not line[15:60].split()[0].isdigit():
1284 if not hasattr(self, 'vibsyms'):
1285 self.vibsyms = []
1286 syms = line.split()
1287 self.vibsyms.extend(syms)
1289 if line[1:15] == "Frequencies --": # note: matches low-precision block, and
1291 if not hasattr(self, 'vibfreqs'):
1292 self.vibfreqs = []
1294 freqs = [utils.float(f) for f in line[15:].split()]
1295 self.vibfreqs.extend(freqs)
1297 if line[1:15] == "Red. masses --": # note: matches only low-precision block
1299 if not hasattr(self, 'vibrmasses'):
1300 self.vibrmasses = []
1302 rmasses = [utils.float(f) for f in line[15:].split()]
1303 self.vibrmasses.extend(rmasses)
1305 if line[1:15] == "Frc consts --": # note: matches only low-precision block
1307 if not hasattr(self, 'vibfconsts'):
1308 self.vibfconsts = []
1310 fconsts = [utils.float(f) for f in line[15:].split()]
1311 self.vibfconsts.extend(fconsts)
1313 if line[1:15] == "IR Inten --": # note: matches only low-precision block
1315 if not hasattr(self, 'vibirs'):
1316 self.vibirs = []
1318 irs = []
1319 for ir in line[15:].split():
1320 try:
1321 irs.append(utils.float(ir))
1322 except ValueError:
1323 irs.append(utils.float('nan'))
1324 self.vibirs.extend(irs)
1326 if line[1:15] == "Raman Activ --": # note: matches only low-precision block
1328 if not hasattr(self, 'vibramans'):
1329 self.vibramans = []
1331 ramans = []
1332 for raman in line[15:].split():
1333 try:
1334 ramans.append(utils.float(raman))
1335 except ValueError:
1336 ramans.append(utils.float('nan'))
1338 self.vibramans.extend(ramans)
1340 # Block with (default-precision) displacements should start with this.
1341 # 1 2 3
1342 # A A A
1343 # Frequencies -- 370.7936 370.7987 618.0103
1344 # Red. masses -- 2.3022 2.3023 1.9355
1345 # Frc consts -- 0.1865 0.1865 0.4355
1346 # IR Inten -- 0.0000 0.0000 0.0000
1347 # Atom AN X Y Z X Y Z X Y Z
1348 # 1 6 0.00 0.00 -0.04 0.00 0.00 0.19 0.00 0.00 0.12
1349 # 2 6 0.00 0.00 0.19 0.00 0.00 -0.06 0.00 0.00 -0.12
1351 if line.strip().split()[0:3] == ["Atom", "AN", "X"]:
1352 if not hasattr(self, 'vibdisps'):
1353 self.vibdisps = []
1354 disps = []
1355 if not hasattr(self, 'nqmf'):
1356 self.set_attribute('nqmf', self.natom)
1357 for n in range(self.nqmf):
1358 line = next(inputfile)
1359 numbers = [float(s) for s in line[10:].split()]
1360 N = len(numbers) // 3
1361 if not disps:
1362 for n in range(N):
1363 disps.append([])
1364 for n in range(N):
1365 disps[n].append(numbers[3*n:3*n+3])
1366 self.vibdisps.extend(disps)
1368 # Block with high-precision (freq=hpmodes) displacements should start with this.
1369 # 1 2 3 4 5
1370 # A A A A A
1371 # Frequencies --- 370.7936 370.7987 618.0103 647.7864 647.7895
1372 # Reduced masses --- 2.3022 2.3023 1.9355 6.4600 6.4600
1373 # Force constants --- 0.1865 0.1865 0.4355 1.5971 1.5972
1374 # IR Intensities --- 0.0000 0.0000 0.0000 0.0000 0.0000
1375 # Coord Atom Element:
1376 # 1 1 6 0.00000 0.00000 0.00000 -0.18677 0.05592
1377 # 2 1 6 0.00000 0.00000 0.00000 0.28440 0.21550
1378 # 3 1 6 -0.04497 0.19296 0.11859 0.00000 0.00000
1379 # 1 2 6 0.00000 0.00000 0.00000 0.03243 0.37351
1380 # 2 2 6 0.00000 0.00000 0.00000 0.14503 -0.06117
1381 # 3 2 6 0.18959 -0.05753 -0.11859 0.00000 0.00000
1382 if line.strip().split()[0:3] == ["Coord", "Atom", "Element:"]:
1383 # Wait until very end of parsing to assign vibdispshp to self.vibdisps
1384 # as otherwise the higher precision displacements will be overwritten
1385 # by low precision displacements which are printed further down file
1386 if not hasattr(self, 'vibdispshp'):
1387 self.vibdispshp = []
1389 disps = []
1390 for n in range(3*self.natom):
1391 line = next(inputfile)
1392 numbers = [float(s) for s in line[16:].split()]
1393 atomindex = int(line[4:10])-1 # atom index, starting at zero
1394 numbermodes = len(numbers)
1396 if not disps:
1397 for mode in range(numbermodes):
1398 # For each mode, make list of list [atom][coord_index]
1399 disps.append([[] for x in range(0, self.natom)])
1400 for mode in range(numbermodes):
1401 disps[mode][atomindex].append(numbers[mode])
1402 self.vibdispshp.extend(disps)
1404 line = next(inputfile)
1406 # Electronic transitions.
1407 if line[1:14] == "Excited State":
1409 if not hasattr(self, "etenergies"):
1410 self.etenergies = []
1411 self.etoscs = []
1412 self.etsyms = []
1413 self.etsecs = []
1415 # Need to deal with lines like:
1416 # (restricted calc)
1417 # Excited State 1: Singlet-BU 5.3351 eV 232.39 nm f=0.1695
1418 # (unrestricted calc) (first excited state is 2!)
1419 # Excited State 2: ?Spin -A 0.1222 eV 10148.75 nm f=0.0000
1420 # (Gaussian 09 ZINDO)
1421 # Excited State 1: Singlet-?Sym 2.5938 eV 478.01 nm f=0.0000 <S**2>=0.000
1422 p = re.compile(r":(?P<sym>.*?)(?P<energy>-?\d*\.\d*) eV")
1423 groups = p.search(line).groups()
1424 self.etenergies.append(utils.convertor(utils.float(groups[1]), "eV", "wavenumber"))
1425 self.etoscs.append(utils.float(line.split("f=")[-1].split()[0]))
1426 self.etsyms.append(groups[0].strip())
1428 line = next(inputfile)
1430 p = re.compile(r"(\d+)")
1431 CIScontrib = []
1432 while line.find(" ->") >= 0: # This is a contribution to the transition
1433 parts = line.split("->")
1434 self.logger.debug(parts)
1435 # Has to deal with lines like:
1436 # 32 -> 38 0.04990
1437 # 35A -> 45A 0.01921
1438 frommoindex = 0 # For restricted or alpha unrestricted
1439 fromMO = parts[0].strip()
1440 if fromMO[-1] == "B":
1441 frommoindex = 1 # For beta unrestricted
1442 fromMO = int(p.match(fromMO).group())-1 # subtract 1 so that it is an index into moenergies
1444 t = parts[1].split()
1445 tomoindex = 0
1446 toMO = t[0]
1447 if toMO[-1] == "B":
1448 tomoindex = 1
1449 toMO = int(p.match(toMO).group())-1 # subtract 1 so that it is an index into moenergies
1451 percent = utils.float(t[1])
1452 # For restricted calculations, the percentage will be corrected
1453 # after parsing (see after_parsing() above).
1454 CIScontrib.append([(fromMO, frommoindex), (toMO, tomoindex), percent])
1455 line = next(inputfile)
1456 self.etsecs.append(CIScontrib)
1458 # Electronic transition transition-dipole data
1459 #
1460 # Ground to excited state transition electric dipole moments (Au):
1461 # state X Y Z Dip. S. Osc.
1462 # 1 0.0008 -0.0963 0.0000 0.0093 0.0005
1463 # 2 0.0553 0.0002 0.0000 0.0031 0.0002
1464 # 3 0.0019 2.3193 0.0000 5.3790 0.4456
1465 # Ground to excited state transition velocity dipole moments (Au):
1466 # state X Y Z Dip. S. Osc.
1467 # 1 -0.0001 0.0032 0.0000 0.0000 0.0001
1468 # 2 -0.0081 0.0000 0.0000 0.0001 0.0005
1469 # 3 -0.0002 -0.2692 0.0000 0.0725 0.3887
1470 # Ground to excited state transition magnetic dipole moments (Au):
1471 # state X Y Z
1472 # 1 0.0000 0.0000 -0.0003
1473 # 2 0.0000 0.0000 0.0000
1474 # 3 0.0000 0.0000 0.0035
1475 #
1476 # NOTE: In Gaussian03, there were some inconsitancies in the use of
1477 # small / capital letters: e.g.
1478 # Ground to excited state Transition electric dipole moments (Au)
1479 # Ground to excited state transition velocity dipole Moments (Au)
1480 # so to look for a match, we will lower() everything.
1482 if line[1:51].lower() == "ground to excited state transition electric dipole":
1483 if not hasattr(self, "etdips"):
1484 self.etdips = []
1485 self.etveldips = []
1486 self.etmagdips = []
1487 if self.etdips == []:
1488 self.netroot = 0
1489 etrootcount = 0 # to count number of et roots
1491 # now loop over lines reading eteltrdips until we find eteltrdipvel
1492 line = next(inputfile) # state X ...
1493 line = next(inputfile) # 1 -0.0001 ...
1494 while line[1:40].lower() != "ground to excited state transition velo":
1495 self.etdips.append(list(map(float, line.split()[1:4])))
1496 etrootcount += 1
1497 line = next(inputfile)
1498 if not self.netroot:
1499 self.netroot = etrootcount
1501 # now loop over lines reading etveleltrdips until we find
1502 # etmagtrdip
1503 line = next(inputfile) # state X ...
1504 line = next(inputfile) # 1 -0.0001 ...
1505 while line[1:40].lower() != "ground to excited state transition magn":
1506 self.etveldips.append(list(map(float, line.split()[1:4])))
1507 line = next(inputfile)
1509 # now loop over lines while the line starts with at least 3 spaces
1510 line = next(inputfile) # state X ...
1511 line = next(inputfile) # 1 -0.0001 ...
1512 while line[0:3] == " ":
1513 self.etmagdips.append(list(map(float, line.split()[1:4])))
1514 line = next(inputfile)
1516 # Circular dichroism data (different for G03 vs G09)
1517 #
1518 # G03
1519 #
1520 # ## <0|r|b> * <b|rxdel|0> (Au), Rotatory Strengths (R) in
1521 # ## cgs (10**-40 erg-esu-cm/Gauss)
1522 # ## state X Y Z R(length)
1523 # ## 1 0.0006 0.0096 -0.0082 -0.4568
1524 # ## 2 0.0251 -0.0025 0.0002 -5.3846
1525 # ## 3 0.0168 0.4204 -0.3707 -15.6580
1526 # ## 4 0.0721 0.9196 -0.9775 -3.3553
1527 #
1528 # G09
1529 #
1530 # ## 1/2[<0|r|b>*<b|rxdel|0> + (<0|rxdel|b>*<b|r|0>)*]
1531 # ## Rotatory Strengths (R) in cgs (10**-40 erg-esu-cm/Gauss)
1532 # ## state XX YY ZZ R(length) R(au)
1533 # ## 1 -0.3893 -6.7546 5.7736 -0.4568 -0.0010
1534 # ## 2 -17.7437 1.7335 -0.1435 -5.3845 -0.0114
1535 # ## 3 -11.8655 -297.2604 262.1519 -15.6580 -0.0332
1536 if line[1:52] == "<0|r|b> * <b|rxdel|0> (Au), Rotatory Strengths (R)" or \
1537 line[1:50] == "1/2[<0|r|b>*<b|rxdel|0> + (<0|rxdel|b>*<b|r|0>)*]":
1539 self.etrotats = []
1541 self.skip_lines(inputfile, ['units'])
1543 headers = next(inputfile)
1544 Ncolms = len(headers.split())
1545 line = next(inputfile)
1546 parts = line.strip().split()
1547 while len(parts) == Ncolms:
1548 try:
1549 R = utils.float(parts[4])
1550 except ValueError:
1551 # nan or -nan if there is no first excited state
1552 # (for unrestricted calculations)
1553 pass
1554 else:
1555 self.etrotats.append(R)
1556 line = next(inputfile)
1557 temp = line.strip().split()
1558 parts = line.strip().split()
1559 self.etrotats = numpy.array(self.etrotats, "d")
1561 # Number of basis sets functions.
1562 # Has to deal with lines like:
1563 # NBasis = 434 NAE= 97 NBE= 97 NFC= 34 NFV= 0
1564 # and...
1565 # NBasis = 148 MinDer = 0 MaxDer = 0
1566 # Although the former is in every file, it doesn't occur before
1567 # the overlap matrix is printed.
1568 if line[1:7] == "NBasis" or line[4:10] == "NBasis":
1570 # For counterpoise fragment, skip these lines.
1571 if self.counterpoise != 0:
1572 return
1574 # For ONIOM calcs, ignore this section in order to bypass assertion failure.
1575 if self.oniom:
1576 return
1578 # If nbasis was already parsed, check if it changed. If it did, issue a warning.
1579 # In the future, we will probably want to have nbasis, as well as nmo below,
1580 # as a list so that we don't need to pick one value when it changes.
1581 nbasis = int(line.split('=')[1].split()[0])
1582 if hasattr(self, "nbasis"):
1583 try:
1584 assert nbasis == self.nbasis
1585 except AssertionError:
1586 self.logger.warning("Number of basis functions (nbasis) has changed from %i to %i" % (self.nbasis, nbasis))
1587 self.nbasis = nbasis
1589 # Number of linearly-independent basis sets.
1590 if line[1:7] == "NBsUse":
1592 # For counterpoise fragment, skip these lines.
1593 if self.counterpoise != 0:
1594 return
1596 # For ONIOM calcs, ignore this section in order to bypass assertion failure.
1597 if self.oniom:
1598 return
1600 nmo = int(line.split('=')[1].split()[0])
1601 self.set_attribute('nmo', nmo)
1603 # For AM1 calculations, set nbasis by a second method,
1604 # as nmo may not always be explicitly stated.
1605 if line[7:22] == "basis functions, ":
1607 nbasis = int(line.split()[0])
1608 self.set_attribute('nbasis', nbasis)
1610 # Molecular orbital overlap matrix.
1611 # Has to deal with lines such as:
1612 # *** Overlap ***
1613 # ****** Overlap ******
1614 # Note that Gaussian sometimes drops basis functions,
1615 # causing the overlap matrix as parsed below to not be
1616 # symmetric (which is a problem for population analyses, etc.)
1617 if line[1:4] == "***" and (line[5:12] == "Overlap" or line[8:15] == "Overlap"):
1619 # Ensure that this is the main calc and not a fragment
1620 if self.counterpoise != 0:
1621 return
1623 self.aooverlaps = numpy.zeros((self.nbasis, self.nbasis), "d")
1624 # Overlap integrals for basis fn#1 are in aooverlaps[0]
1625 base = 0
1626 colmNames = next(inputfile)
1627 while base < self.nbasis:
1629 self.updateprogress(inputfile, "Overlap", self.fupdate)
1631 for i in range(self.nbasis-base): # Fewer lines this time
1632 line = next(inputfile)
1633 parts = line.split()
1634 for j in range(len(parts)-1): # Some lines are longer than others
1635 k = float(parts[j+1].replace("D", "E"))
1636 self.aooverlaps[base+j, i+base] = k
1637 self.aooverlaps[i+base, base+j] = k
1638 base += 5
1639 colmNames = next(inputfile)
1640 self.aooverlaps = numpy.array(self.aooverlaps, "d")
1642 # Molecular orbital coefficients (mocoeffs).
1643 # Essentially only produced for SCF calculations.
1644 # This is also the place where aonames and atombasis are parsed.
1645 if line[5:35] == "Molecular Orbital Coefficients" or line[5:41] == "Alpha Molecular Orbital Coefficients" or line[5:40] == "Beta Molecular Orbital Coefficients":
1647 # If counterpoise fragment, return without parsing orbital info
1648 if self.counterpoise != 0:
1649 return
1650 # Skip this for ONIOM calcs
1651 if self.oniom:
1652 return
1654 if line[5:40] == "Beta Molecular Orbital Coefficients":
1655 beta = True
1656 if self.popregular:
1657 return
1658 # This was continue before refactoring the parsers.
1659 #continue # Not going to extract mocoeffs
1660 # Need to add an extra array to self.mocoeffs
1661 self.mocoeffs.append(numpy.zeros((self.nmo, self.nbasis), "d"))
1662 else:
1663 beta = False
1664 self.aonames = []
1665 self.atombasis = []
1666 mocoeffs = [numpy.zeros((self.nmo, self.nbasis), "d")]
1668 base = 0
1669 self.popregular = False
1670 for base in range(0, self.nmo, 5):
1672 self.updateprogress(inputfile, "Coefficients", self.fupdate)
1674 colmNames = next(inputfile)
1676 if not colmNames.split():
1677 self.logger.warning("Molecular coefficients header found but no coefficients.")
1678 break
1680 if base == 0 and int(colmNames.split()[0]) != 1:
1681 # Implies that this is a POP=REGULAR calculation
1682 # and so, only aonames (not mocoeffs) will be extracted
1683 self.popregular = True
1684 symmetries = next(inputfile)
1685 eigenvalues = next(inputfile)
1686 for i in range(self.nbasis):
1688 line = next(inputfile)
1689 if i == 0:
1690 # Find location of the start of the basis function name
1691 start_of_basis_fn_name = line.find(line.split()[3]) - 1
1692 if base == 0 and not beta: # Just do this the first time 'round
1693 parts = line[:start_of_basis_fn_name].split()
1694 if len(parts) > 1: # New atom
1695 if i > 0:
1696 self.atombasis.append(atombasis)
1697 atombasis = []
1698 atomname = "%s%s" % (parts[2], parts[1])
1699 orbital = line[start_of_basis_fn_name:20].strip()
1700 self.aonames.append("%s_%s" % (atomname, orbital))
1701 atombasis.append(i)
1703 part = line[21:].replace("D", "E").rstrip()
1704 temp = []
1705 for j in range(0, len(part), 10):
1706 temp.append(float(part[j:j+10]))
1707 if beta:
1708 self.mocoeffs[1][base:base + len(part) // 10, i] = temp
1709 else:
1710 mocoeffs[0][base:base + len(part) // 10, i] = temp
1712 if base == 0 and not beta: # Do the last update of atombasis
1713 self.atombasis.append(atombasis)
1714 if self.popregular:
1715 # We now have aonames, so no need to continue
1716 break
1717 if not self.popregular and not beta:
1718 self.mocoeffs = mocoeffs
1720 # Natural orbital coefficients (nocoeffs) and occupation numbers (nooccnos),
1721 # which are respectively define the eigenvectors and eigenvalues of the
1722 # diagonalized one-electron density matrix. These orbitals are formed after
1723 # configuration interaction (CI) calculations, but not only. Similarly to mocoeffs,
1724 # we can parse and check aonames and atombasis here.
1725 #
1726 # Natural Orbital Coefficients:
1727 # 1 2 3 4 5
1728 # Eigenvalues -- 2.01580 2.00363 2.00000 2.00000 1.00000
1729 # 1 1 O 1S 0.00000 -0.15731 -0.28062 0.97330 0.00000
1730 # 2 2S 0.00000 0.75440 0.57746 0.07245 0.00000
1731 # ...
1732 #
1733 def natural_orbital_single_spin_parsing(inputfile, updateprogress_title):
1734 coeffs = numpy.zeros((self.nmo, self.nbasis), "d")
1735 occnos = []
1736 aonames = []
1737 atombasis = []
1738 for base in range(0, self.nmo, 5):
1739 self.updateprogress(inputfile, updateprogress_title, self.fupdate)
1740 colmNames = next(inputfile)
1741 eigenvalues = next(inputfile)
1742 occnos.extend(map(float, eigenvalues.split()[2:]))
1743 for i in range(self.nbasis):
1744 line = next(inputfile)
1745 # Just do this the first time 'round.
1746 if base == 0:
1747 parts = line[:12].split()
1748 # New atom.
1749 if len(parts) > 1:
1750 if i > 0:
1751 atombasis.append(basisonatom)
1752 basisonatom = []
1753 atomname = "%s%s" % (parts[2], parts[1])
1754 orbital = line[11:20].strip()
1755 aonames.append("%s_%s" % (atomname, orbital))
1756 basisonatom.append(i)
1757 part = line[21:].replace("D", "E").rstrip()
1758 temp = []
1759 for j in range(0, len(part), 10):
1760 temp.append(float(part[j:j+10]))
1761 coeffs[base:base + len(part) // 10, i] = temp
1762 # Do the last update of atombasis.
1763 if base == 0:
1764 atombasis.append(basisonatom)
1765 return occnos, coeffs, aonames, atombasis
1767 if line[5:33] == "Natural Orbital Coefficients":
1768 updateprogress_title = "Natural orbitals"
1769 nooccnos, nocoeffs, aonames, atombasis = natural_orbital_single_spin_parsing(inputfile, updateprogress_title)
1770 self.set_attribute("nocoeffs", nocoeffs)
1771 self.set_attribute("nooccnos", nooccnos)
1772 self.set_attribute("atombasis", atombasis)
1773 self.set_attribute("aonames", aonames)
1775 # Natural spin orbital coefficients (nsocoeffs) and occupation numbers (nsooccnos)
1776 # Parsed attributes are similar to the natural orbitals above except
1777 # the natural spin orbitals and occupation numbers are the eigenvalues
1778 # and eigenvectors of the one particles spin density matrices
1779 # Alpha Natural Orbital Coefficients:
1780 # 1 2 3 4 5
1781 # Eigenvalues -- 1.00000 1.00000 0.99615 0.99320 0.99107
1782 # 1 1 O 1S 0.70425 0.70600 -0.16844 -0.14996 -0.00000
1783 # 2 2S 0.01499 0.01209 0.36089 0.34940 -0.00000
1784 # ...
1785 # Beta Natural Orbital Coefficients:
1786 # 1 2 3 4 5
1787 # Eigenvalues -- 1.00000 1.00000 0.99429 0.98790 0.98506
1788 # 1 1 O 1S 0.70822 0.70798 -0.15316 -0.13458 0.00465
1789 # 2 2S 0.00521 0.00532 0.33837 0.33189 -0.01301
1790 # 3 3S -0.02542 -0.00841 0.28649 0.53224 0.18902
1791 # ...
1793 if line[5:39] == "Alpha Natural Orbital Coefficients":
1794 updateprogress_title = "Natural Spin orbitals (alpha)"
1795 nsooccnos, nsocoeffs, aonames, atombasis = natural_orbital_single_spin_parsing(inputfile, updateprogress_title)
1796 if self.unified_no_nso:
1797 self.append_attribute("nocoeffs", nsocoeffs)
1798 self.append_attribute("nooccnos", nsooccnos)
1799 else:
1800 self.append_attribute("nsocoeffs", nsocoeffs)
1801 self.append_attribute("nsooccnos", nsooccnos)
1802 self.set_attribute("atombasis", atombasis)
1803 self.set_attribute("aonames", aonames)
1804 if line[5:38] == "Beta Natural Orbital Coefficients":
1805 updateprogress_title = "Natural Spin orbitals (beta)"
1806 nsooccnos, nsocoeffs, aonames, atombasis = natural_orbital_single_spin_parsing(inputfile, updateprogress_title)
1807 if self.unified_no_nso:
1808 self.append_attribute("nocoeffs", nsocoeffs)
1809 self.append_attribute("nooccnos", nsooccnos)
1810 else:
1811 self.append_attribute("nsocoeffs", nsocoeffs)
1812 self.append_attribute("nsooccnos", nsooccnos)
1813 self.set_attribute("atombasis", atombasis)
1814 self.set_attribute("aonames", aonames)
1816 # For FREQ=Anharm, extract anharmonicity constants
1817 if line[1:40] == "X matrix of Anharmonic Constants (cm-1)":
1818 Nvibs = len(self.vibfreqs)
1819 self.vibanharms = numpy.zeros((Nvibs, Nvibs), "d")
1821 base = 0
1822 colmNames = next(inputfile)
1823 while base < Nvibs:
1825 for i in range(Nvibs-base): # Fewer lines this time
1826 line = next(inputfile)
1827 parts = line.split()
1828 for j in range(len(parts)-1): # Some lines are longer than others
1829 k = float(parts[j+1].replace("D", "E"))
1830 self.vibanharms[base+j, i+base] = k
1831 self.vibanharms[i+base, base+j] = k
1832 base += 5
1833 colmNames = next(inputfile)
1835 # Pseudopotential charges.
1836 if line.find("Pseudopotential Parameters") > -1:
1838 self.skip_lines(inputfile, ['e', 'label1', 'label2', 'e'])
1840 line = next(inputfile)
1841 if line.find("Centers:") < 0:
1842 return
1843 # This was continue before parser refactoring.
1844 # continue
1846 # Needs to handle code like the following:
1847 #
1848 # Center Atomic Valence Angular Power Coordinates
1849 # Number Number Electrons Momentum of R Exponent Coefficient X Y Z
1850 # ===================================================================================================================================
1851 # Centers: 1
1852 # Centers: 16
1853 # Centers: 21 24
1854 # Centers: 99100101102
1855 # 1 44 16 -4.012684 -0.696698 0.006750
1856 # F and up
1857 # 0 554.3796303 -0.05152700
1858 centers = []
1859 while line.find("Centers:") >= 0:
1860 temp = line[10:]
1861 for i in range(0, len(temp)-3, 3):
1862 centers.append(int(temp[i:i+3]))
1863 line = next(inputfile)
1864 centers.sort() # Not always in increasing order
1866 self.coreelectrons = numpy.zeros(self.natom, "i")
1868 for center in centers:
1869 front = line[:10].strip()
1870 while not (front and int(front) == center):
1871 line = next(inputfile)
1872 front = line[:10].strip()
1873 info = line.split()
1874 self.coreelectrons[center-1] = int(info[1]) - int(info[2])
1875 line = next(inputfile)
1877 # This will be printed for counterpoise calcualtions only.
1878 # To prevent crashing, we need to know which fragment is being considered.
1879 # Other information is also printed in lines that start like this.
1880 if line[1:14] == 'Counterpoise:':
1882 if line[42:50] == "fragment":
1883 self.counterpoise = int(line[51:54])
1885 # This will be printed only during ONIOM calcs; use it to set a flag
1886 # that will allow assertion failures to be bypassed in the code.
1887 if line[1:7] == "ONIOM:":
1888 self.oniom = True
1890 # This will be printed only during BOMD calcs;
1891 if line.startswith(" INPUT DATA FOR L118"):
1892 self.BOMD = True
1894 # Atomic charges are straightforward to parse, although the header
1895 # has changed over time somewhat.
1896 #
1897 # Mulliken charges:
1898 # 1
1899 # 1 C -0.004513
1900 # 2 C -0.077156
1901 # ...
1902 # Sum of Mulliken charges = 0.00000
1903 # Mulliken charges with hydrogens summed into heavy atoms:
1904 # 1
1905 # 1 C -0.004513
1906 # 2 C 0.002063
1907 # ...
1908 #
1909 if line[1:25] == "Mulliken atomic charges:" or line[1:18] == "Mulliken charges:" or \
1910 line[1:23] == "Lowdin Atomic Charges:" or line[1:16] == "Lowdin charges:" or \
1911 line[1:37] == "Mulliken charges and spin densities:" or \
1912 line[1:32] == "Mulliken atomic spin densities:":
1914 has_spin = 'spin densities' in line
1915 has_charges = 'charges' in line
1917 if has_charges and not hasattr(self, "atomcharges"):
1918 self.atomcharges = {}
1920 if has_spin and not hasattr(self, "atomspins"):
1921 self.atomspins = {}
1923 ones = next(inputfile)
1925 charges = []
1926 spins = []
1927 nline = next(inputfile)
1928 while not "Sum of" in nline:
1929 if has_charges:
1930 charges.append(float(nline.split()[2]))
1932 if has_spin and has_charges:
1933 spins.append(float(nline.split()[3]))
1935 if has_spin and not has_charges:
1936 spins.append(float(nline.split()[2]))
1938 nline = next(inputfile)
1940 if "Mulliken" in line:
1941 if has_charges:
1942 self.atomcharges["mulliken"] = charges
1943 if has_spin:
1944 self.atomspins["mulliken"] = spins
1946 elif "Lowdin" in line:
1947 self.atomcharges["lowdin"] = charges
1950 if line.strip() == "Natural Population":
1951 if not hasattr(self, 'atomcharges'):
1952 self.atomcharges = {}
1953 line1 = next(inputfile)
1954 line2 = next(inputfile)
1955 if line1.split()[0] == 'Natural' and line2.split()[2] == 'Charge':
1956 dashes = next(inputfile)
1957 charges = []
1958 for i in range(self.natom):
1959 nline = next(inputfile)
1960 charges.append(float(nline.split()[2]))
1961 self.atomcharges["natural"] = charges
1963 #Extract Thermochemistry
1964 #Temperature 298.150 Kelvin. Pressure 1.00000 Atm.
1965 #Zero-point correction= 0.342233 (Hartree/
1966 #Thermal correction to Energy= 0.
1967 #Thermal correction to Enthalpy= 0.
1968 #Thermal correction to Gibbs Free Energy= 0.302940
1969 #Sum of electronic and zero-point Energies= -563.649744
1970 #Sum of electronic and thermal Energies= -563.636699
1971 #Sum of electronic and thermal Enthalpies= -563.635755
1972 #Sum of electronic and thermal Free Energies= -563.689037
1973 if "Zero-point correction" in line:
1974 self.set_attribute('zpve', float(line.split()[2]))
1975 if "Sum of electronic and thermal Enthalpies" in line:
1976 self.set_attribute('enthalpy', float(line.split()[6]))
1977 if "Sum of electronic and thermal Free Energies=" in line:
1978 self.set_attribute('freeenergy', float(line.split()[7]))
1979 if line[1:13] == "Temperature ":
1980 self.set_attribute('temperature', float(line.split()[1]))
1981 self.set_attribute('pressure', float(line.split()[4]))
1983 # Static polarizability (from `polar`), lower triangular
1984 # matrix.
1985 if line[1:26] == "SCF Polarizability for W=":
1986 self.hp_polarizabilities = True
1987 if not hasattr(self, 'polarizabilities'):
1988 self.polarizabilities = []
1989 polarizability = numpy.zeros(shape=(3, 3))
1990 self.skip_line(inputfile, 'directions')
1991 for i in range(3):
1992 line = next(inputfile)
1993 polarizability[i, :i+1] = [utils.float(x) for x in line.split()[1:]]
1995 polarizability = utils.symmetrize(polarizability, use_triangle='lower')
1996 self.polarizabilities.append(polarizability)
1998 # Static polarizability (from `freq`), lower triangular matrix.
1999 if line[1:16] == "Polarizability=":
2000 self.hp_polarizabilities = True
2001 if not hasattr(self, 'polarizabilities'):
2002 self.polarizabilities = []
2003 polarizability = numpy.zeros(shape=(3, 3))
2004 polarizability_list = []
2005 polarizability_list.extend([line[16:31], line[31:46], line[46:61]])
2006 line = next(inputfile)
2007 polarizability_list.extend([line[16:31], line[31:46], line[46:61]])
2008 indices = numpy.tril_indices(3)
2009 polarizability[indices] = [utils.float(x) for x in polarizability_list]
2010 polarizability = utils.symmetrize(polarizability, use_triangle='lower')
2011 self.polarizabilities.append(polarizability)
2013 # Static polarizability, compressed into a single line from
2014 # terse printing.
2015 # Order is XX, YX, YY, ZX, ZY, ZZ (lower triangle).
2016 if line[2:23] == "Exact polarizability:":
2017 if not self.hp_polarizabilities:
2018 if not hasattr(self, 'polarizabilities'):
2019 self.polarizabilities = []
2020 polarizability = numpy.zeros(shape=(3, 3))
2021 indices = numpy.tril_indices(3)
2022 polarizability[indices] = [utils.float(x) for x in
2023 [line[23:31], line[31:39], line[39:47], line[47:55], line[55:63], line[63:71]]]
2024 polarizability = utils.symmetrize(polarizability, use_triangle='lower')
2025 self.polarizabilities.append(polarizability)
2027 # IRC Computation convergence checks.
2028 #
2029 # -------- Sample extract for IRC step --------
2030 #
2031 # IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC
2032 # ------------------------------------------------------------------------
2033 # INPUT DATA FOR L123
2034 # ------------------------------------------------------------------------
2035 # GENERAL PARAMETERS:
2036 # Follow reaction path in both directions.
2037 # Maximum points per path = 200
2038 # Step size = 0.100 bohr
2039 # Integration scheme = HPC
2040 # Redo corrector integration= Yes
2041 # DWI Weight Power = 2
2042 # DWI will use Hessian update vectors when possible.
2043 # Max correction cycles = 50
2044 # Initial Hessian = CalcFC
2045 # Hessian evaluation = Analytic every 5 predictor steps
2046 # = Analytic every 5 corrector steps
2047 # Hessian updating method = Bofill
2048 # ------------------------------------------------------------------------
2049 # IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC
2050 #
2051 # -------- Sample extract for converged step --------
2052 # IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC
2053 # Pt 1 Step number 1 out of a maximum of 50
2054 # Modified Bulirsch-Stoer Extrapolation Cycles:
2055 # EPS = 0.000010000000000
2056 # Maximum DWI energy std dev = 0.000000595 at pt 1
2057 # Maximum DWI gradient std dev = 0.135684493 at pt 2
2058 # CORRECTOR INTEGRATION CONVERGENCE:
2059 # Recorrection delta-x convergence threshold: 0.010000
2060 # Delta-x Convergence Met
2061 # Point Number: 1 Path Number: 1
2062 # CHANGE IN THE REACTION COORDINATE = 0.16730
2063 # NET REACTION COORDINATE UP TO THIS POINT = 0.16730
2064 # # OF POINTS ALONG THE PATH = 1
2065 # # OF STEPS = 1
2066 #
2067 # Calculating another point on the path.
2068 # IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC
2069 #
2070 # -------- Sample extract for unconverged intermediate step --------
2071 # IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC
2072 # Error in corrector energy = -0.0000457166
2073 # Magnitude of corrector gradient = 0.0140183779
2074 # Magnitude of analytic gradient = 0.0144021969
2075 # Magnitude of difference = 0.0078709968
2076 # Angle between gradients (degrees)= 32.1199
2077 # Pt 40 Step number 2 out of a maximum of 20
2078 # Modified Bulirsch-Stoer Extrapolation Cycles:
2079 # EPS = 0.000010000000000
2080 # Maximum DWI energy std dev = 0.000007300 at pt 31
2081 # Maximum DWI gradient std dev = 0.085197906 at pt 59
2082 # CORRECTOR INTEGRATION CONVERGENCE:
2083 # Recorrection delta-x convergence threshold: 0.010000
2084 # Delta-x Convergence NOT Met
2085 # IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC
2086 if line[1:20] == "INPUT DATA FOR L123": # First IRC step
2087 if not hasattr(self, "optstatus"):
2088 self.optstatus = []
2089 self.optstatus.append(data.ccData.OPT_NEW)
2090 if line[3:22] == "Delta-x Convergence":
2091 if line[23:30] == "NOT Met":
2092 self.optstatus[-1] += data.ccData.OPT_UNCONVERGED
2093 elif line[23:26] == "Met":
2094 self.optstatus[-1] += data.ccData.OPT_DONE
2095 if not hasattr(self, 'optdone'):
2096 self.optdone = []
2097 self.optdone.append(len(self.optstatus) - 1)
2099 if line[:31] == ' Normal termination of Gaussian':
2100 self.metadata['success'] = True