# -*- coding: utf-8 -*- # # Copyright (c) 2018, the cclib development team # # This file is part of cclib (http://cclib.github.io) and is distributed under # the terms of the BSD 3-Clause License. """Parser for Gaussian output files""" from __future__ import print_function import re import numpy from cclib.parser import data from cclib.parser import logfileparser from cclib.parser import utils class Gaussian(logfileparser.Logfile): """A Gaussian 98/03 log file.""" def __init__(self, *args, **kwargs): # Call the __init__ method of the superclass super(Gaussian, self).__init__(logname="Gaussian", *args, **kwargs) def __str__(self): """Return a string representation of the object.""" return "Gaussian log file %s" % (self.filename) def __repr__(self): """Return a representation of the object.""" return 'Gaussian("%s")' % (self.filename) def normalisesym(self, label): """Use standard symmetry labels instead of Gaussian labels. To normalise: (1) If label is one of [SG, PI, PHI, DLTA], replace by [sigma, pi, phi, delta] (2) replace any G or U by their lowercase equivalent """ # note: DLT must come after DLTA greek = [('SG', 'sigma'), ('PI', 'pi'), ('PHI', 'phi'), ('DLTA', 'delta'), ('DLT', 'delta')] for k, v in greek: if label.startswith(k): tmp = label[len(k):] label = v if tmp: label = v + "." + tmp ans = label.replace("U", "u").replace("G", "g") return ans def before_parsing(self): # Used to index self.scftargets[]. SCFRMS, SCFMAX, SCFENERGY = list(range(3)) # Extract only well-formed numbers in scientific notation. self.re_scinot = re.compile('(\w*)=\s*(-?\d\.\d{2}D[+-]\d{2})') # Extract only well-formed numbers in traditional # floating-point format. self.re_float = re.compile('(\w*-?\w*)=\s*(-?\d+\.\d{10,})') # Flag for identifying Coupled Cluster runs. self.coupledcluster = False # Fragment number for counterpoise or fragment guess calculations # (normally zero). self.counterpoise = 0 # Flag for identifying ONIOM calculations. self.oniom = False # Flag for identifying BOMD calculations. # These calculations have a back-integration algorithm so that not all # geometries should be kept. # We also add a "time" attribute to the parser. self.BOMD = False # Do we have high-precision polarizabilities printed from a # dedicated `polar` job? If so, avoid duplicate parsing. self.hp_polarizabilities = False def after_parsing(self): # Correct the percent values in the etsecs in the case of # a restricted calculation. The following has the # effect of including each transition twice. if hasattr(self, "etsecs") and len(self.homos) == 1: new_etsecs = [[(x[0], x[1], x[2] * numpy.sqrt(2)) for x in etsec] for etsec in self.etsecs] self.etsecs = new_etsecs if hasattr(self, "scanenergies"): self.scancoords = [] self.scancoords = self.atomcoords if (hasattr(self, 'enthalpy') and hasattr(self, 'temperature') and hasattr(self, 'freeenergy')): self.set_attribute('entropy', (self.enthalpy - self.freeenergy) / self.temperature) # This bit is needed in order to trim coordinates that are printed a second time # at the end of geometry optimizations. Note that we need to do this for both atomcoords # and inputcoords. The reason is that normally a standard orientation is printed and that # is what we parse into atomcoords, but inputcoords stores the input (unmodified) coordinates # and that is copied over to atomcoords if no standard oritentation was printed, which happens # for example for jobs with no symmetry. This last step, however, is now generic for all parsers. # Perhaps then this part should also be generic code... # Regression that tests this: Gaussian03/cyclopropenyl.rhf.g03.cut.log if hasattr(self, 'optstatus') and len(self.optstatus) > 0: last_point = len(self.optstatus) - 1 if hasattr(self, 'atomcoords'): self.atomcoords = self.atomcoords[:last_point + 1] if hasattr(self, 'inputcoords'): self.inputcoords = self.inputcoords[:last_point + 1] # If we parsed high-precision vibrational displacements, overwrite # lower-precision displacements in self.vibdisps if hasattr(self, 'vibdispshp'): self.vibdisps = self.vibdispshp del self.vibdispshp if hasattr(self, 'time'): self.time = [self.time[i] for i in sorted(self.time.keys())] if hasattr(self, 'energies_BOMD'): self.set_attribute('scfenergies', [self.energies_BOMD[i] for i in sorted(self.energies_BOMD.keys())]) if hasattr(self, 'atomcoords_BOMD'): self.atomcoords= \ [self.atomcoords_BOMD[i] for i in sorted(self.atomcoords_BOMD.keys())] def extract(self, inputfile, line): """Extract information from the file object inputfile.""" # Extract the version number: "Gaussian 09, Revision D.01" # becomes "09revisionD.01". if line.strip() == "Cite this work as:": line = inputfile.next() tokens = line.split() self.metadata["package_version"] = ''.join([ tokens[1][:-1], 'revision', tokens[-1][:-1], ]) if line.strip().startswith("Link1: Proceeding to internal job step number"): self.new_internal_job() # This block contains some general information as well as coordinates, # which could be parsed in the future: # # Symbolic Z-matrix: # Charge = 0 Multiplicity = 1 # C 0.73465 0. 0. # C 1.93465 0. 0. # C # ... # # It also lists fragments, if there are any, which is potentially valuable: # # Symbolic Z-matrix: # Charge = 0 Multiplicity = 1 in supermolecule # Charge = 0 Multiplicity = 1 in fragment 1. # Charge = 0 Multiplicity = 1 in fragment 2. # B(Fragment=1) 0.06457 -0.0279 0.01364 # H(Fragment=1) 0.03117 -0.02317 1.21604 # ... # # Note, however, that currently we only parse information for the whole system # or supermolecule as Gaussian calls it. if line.strip() == "Symbolic Z-matrix:": self.updateprogress(inputfile, "Symbolic Z-matrix", self.fupdate) line = inputfile.next() while line.split()[0] == 'Charge': # For the supermolecule, we can parse the charge and multicplicity. regex = ".*=(.*)Mul.*=\s*-?(\d+).*" match = re.match(regex, line) assert match, "Something unusual about the line: '%s'" % line self.set_attribute('charge', int(match.groups()[0])) self.set_attribute('mult', int(match.groups()[1])) if line.split()[-2] == "fragment": self.nfragments = int(line.split()[-1].strip('.')) if line.strip()[-13:] == "model system.": self.nmodels = getattr(self, 'nmodels', 0) + 1 line = inputfile.next() # The remaining part will allow us to get the atom count. # When coordinates are given, there is a blank line at the end, but if # there is a Z-matrix here, there will also be variables and we need to # stop at those to get the right atom count. # Also, in older versions there is bo blank line (G98 regressions), # so we need to watch out for leaving the link. natom = 0 while line.split() and not "Variables" in line and not "Leave Link" in line: natom += 1 line = inputfile.next() self.set_attribute('natom', natom) # Continuing from above, there is not always a symbolic matrix, for example # if the Z-matrix was in the input file. In such cases, try to match the # line and get at the charge and multiplicity. # # Charge = 0 Multiplicity = 1 in supermolecule # Charge = 0 Multiplicity = 1 in fragment 1. # Charge = 0 Multiplicity = 1 in fragment 2. if line[1:7] == 'Charge' and line.find("Multiplicity") >= 0: self.updateprogress(inputfile, "Charge and Multiplicity", self.fupdate) if line.split()[-1] == "supermolecule" or not "fragment" in line: regex = ".*=(.*)Mul.*=\s*-?(\d+).*" match = re.match(regex, line) assert match, "Something unusual about the line: '%s'" % line self.set_attribute('charge', int(match.groups()[0])) self.set_attribute('mult', int(match.groups()[1])) if line.split()[-2] == "fragment": self.nfragments = int(line.split()[-1].strip('.')) # Number of atoms is also explicitely printed after the above. if line[1:8] == "NAtoms=": self.updateprogress(inputfile, "Attributes", self.fupdate) natom = int(re.search('NAtoms=\s*(\d+)', line).group(1)) self.set_attribute('natom', natom) # Basis set name if line[1:15] == "Standard basis": self.metadata["basis_set"] = line.split()[2] # Dipole moment # e.g. from G09 # Dipole moment (field-independent basis, Debye): # X= 0.0000 Y= 0.0000 Z= 0.0930 # e.g. from G03 # X= 0.0000 Y= 0.0000 Z= -1.6735 Tot= 1.6735 # need the "field independent" part - ONIOM and other calc use diff formats if line[1:39] == "Dipole moment (field-independent basis": self.updateprogress(inputfile, "Dipole and Higher Moments", self.fupdate) self.reference = [0.0, 0.0, 0.0] self.moments = [self.reference] tokens = inputfile.next().split() # split - dipole would need to be *huge* to fail a split # and G03 and G09 use different spacing if len(tokens) >= 6: dipole = (float(tokens[1]), float(tokens[3]), float(tokens[5])) if not hasattr(self, 'moments'): self.moments = [self.reference, dipole] else: self.moments.append(dipole) if line[1:43] == "Quadrupole moment (field-independent basis": # e.g. (g09) # Quadrupole moment (field-independent basis, Debye-Ang): # XX= -6.1213 YY= -4.2950 ZZ= -5.4175 # XY= 0.0000 XZ= 0.0000 YZ= 0.0000 # or from g03 # XX= -6.1213 YY= -4.2950 ZZ= -5.4175 quadrupole = {} for j in range(2): # two rows line = inputfile.next() if line[22] == '=': # g03 file for i in (1, 18, 35): quadrupole[line[i:i+4]] = float(line[i+5:i+16]) else: for i in (1, 27, 53): quadrupole[line[i:i+4]] = float(line[i+5:i+25]) lex = sorted(quadrupole.keys()) quadrupole = [quadrupole[key] for key in lex] if not hasattr(self, 'moments') or len(self.moments) < 2: self.logger.warning("Found quadrupole moments but no previous dipole") self.reference = [0.0, 0.0, 0.0] self.moments = [self.reference, None, quadrupole] else: if len(self.moments) == 2: self.moments.append(quadrupole) else: assert self.moments[2] == quadrupole if line[1:41] == "Octapole moment (field-independent basis": # e.g. # Octapole moment (field-independent basis, Debye-Ang**2): # XXX= 0.0000 YYY= 0.0000 ZZZ= -0.1457 XYY= 0.0000 # XXY= 0.0000 XXZ= 0.0136 XZZ= 0.0000 YZZ= 0.0000 # YYZ= -0.5848 XYZ= 0.0000 octapole = {} for j in range(2): # two rows line = inputfile.next() if line[22] == '=': # g03 file for i in (1, 18, 35, 52): octapole[line[i:i+4]] = float(line[i+5:i+16]) else: for i in (1, 27, 53, 79): octapole[line[i:i+4]] = float(line[i+5:i+25]) # last line only 2 moments line = inputfile.next() if line[22] == '=': # g03 file for i in (1, 18): octapole[line[i:i+4]] = float(line[i+5:i+16]) else: for i in (1, 27): octapole[line[i:i+4]] = float(line[i+5:i+25]) lex = sorted(octapole.keys()) octapole = [octapole[key] for key in lex] if not hasattr(self, 'moments') or len(self.moments) < 3: self.logger.warning("Found octapole moments but no previous dipole or quadrupole") self.reference = [0.0, 0.0, 0.0] self.moments = [self.reference, None, None, octapole] else: if len(self.moments) == 3: self.moments.append(octapole) else: assert self.moments[3] == octapole if line[1:20] == "Hexadecapole moment": # e.g. # Hexadecapole moment (field-independent basis, Debye-Ang**3): # XXXX= -3.2614 YYYY= -6.8264 ZZZZ= -4.9965 XXXY= 0.0000 # XXXZ= 0.0000 YYYX= 0.0000 YYYZ= 0.0000 ZZZX= 0.0000 # ZZZY= 0.0000 XXYY= -1.8585 XXZZ= -1.4123 YYZZ= -1.7504 # XXYZ= 0.0000 YYXZ= 0.0000 ZZXY= 0.0000 hexadecapole = {} # read three lines worth of 4 moments per line for j in range(3): line = inputfile.next() if line[22] == '=': # g03 file for i in (1, 18, 35, 52): hexadecapole[line[i:i+4]] = float(line[i+5:i+16]) else: for i in (1, 27, 53, 79): hexadecapole[line[i:i+4]] = float(line[i+5:i+25]) # last line only 3 moments line = inputfile.next() if line[22] == '=': # g03 file for i in (1, 18, 35): hexadecapole[line[i:i+4]] = float(line[i+5:i+16]) else: for i in (1, 27, 53): hexadecapole[line[i:i+4]] = float(line[i+5:i+25]) lex = sorted(hexadecapole.keys()) hexadecapole = [hexadecapole[key] for key in lex] if not hasattr(self, 'moments') or len(self.moments) < 4: self.reference = [0.0, 0.0, 0.0] self.moments = [self.reference, None, None, None, hexadecapole] else: if len(self.moments) == 4: self.append_attribute("moments", hexadecapole) else: try: numpy.testing.assert_equal(self.moments[4], hexadecapole) except AssertionError: self.logger.warning("Attribute hexadecapole changed value (%s -> %s)" % (self.moments[4], hexadecapole)) self.append_attribute("moments", hexadecapole) # Catch message about completed optimization. if line[1:23] == "Optimization completed": if not hasattr(self, 'optdone'): self.optdone = [] self.optdone.append(len(self.geovalues) - 1) assert hasattr(self, "optstatus") and len(self.optstatus) > 0 self.optstatus[-1] += data.ccData.OPT_DONE # Catch message about stopped optimization (not converged). if line[1:21] == "Optimization stopped": if not hasattr(self, "optdone"): self.optdone = [] assert hasattr(self, "optstatus") and len(self.optstatus) > 0 self.optstatus[-1] += data.ccData.OPT_UNCONVERGED # Extract the atomic numbers and coordinates from the input orientation, # in the event the standard orientation isn't available. # Don't extract from Input or Z-matrix orientation in a BOMD run, as only # the final geometry should be kept but extract inputatoms. if line.find("Input orientation") > -1 or line.find("Z-Matrix orientation") > -1: # If this is a counterpoise calculation, this output means that # the supermolecule is now being considered, so we can set: self.counterpoise = 0 self.updateprogress(inputfile, "Attributes", self.cupdate) if not self.BOMD and not hasattr(self, "inputcoords"): self.inputcoords = [] self.inputatoms = [] self.skip_lines(inputfile, ['d', 'cols', 'cols', 'd']) atomcoords = [] line = next(inputfile) while list(set(line.strip())) != ["-"]: broken = line.split() self.inputatoms.append(int(broken[1])) atomcoords.append(list(map(float, broken[3:6]))) line = next(inputfile) if not self.BOMD: self.inputcoords.append(atomcoords) self.set_attribute('atomnos', self.inputatoms) self.set_attribute('natom', len(self.inputatoms)) if self.BOMD and line.startswith(' Summary information for step'): # We keep time and energies_BOMD and coordinates in a dictionary # because steps can be recalculated, and we need to overwite the # previous data broken = line.split() step = int(broken[-1]) line = next(inputfile) broken = line.split() if not hasattr(self, "time"): self.set_attribute('time', {step:float(broken[-1])}) else: self.time[step] = float(broken[-1]) line = next(inputfile) broken = line.split(';')[1].split() ene = utils.convertor(self.float(broken[-1]), "hartree", "eV") if not hasattr(self, "energies_BOMD"): self.set_attribute('energies_BOMD', {step:ene}) else: self.energies_BOMD[step] = ene self.updateprogress(inputfile, "Attributes", self.cupdate) if not hasattr(self, "atomcoords_BOMD"): self.atomcoords_BOMD = {} #self.inputatoms = [] self.skip_lines(inputfile, ['EKin', 'Angular', 'JX', 'Total', 'Total', 'Cartesian']) atomcoords = [] line = next(inputfile) while not "MW cartesian" in line: broken = line.split() atomcoords.append(list(map(self.float, (broken[3], broken[5], broken[7])))) # self.inputatoms.append(int(broken[1])) line = next(inputfile) self.atomcoords_BOMD[step] = atomcoords #self.set_attribute('atomnos', self.inputatoms) #self.set_attribute('natom', len(self.inputatoms)) # Extract the atomic masses. # Typical section: # Isotopes and Nuclear Properties: #(Nuclear quadrupole moments (NQMom) in fm**2, nuclear magnetic moments (NMagM) # in nuclear magnetons) # # Atom 1 2 3 4 5 6 7 8 9 10 # IAtWgt= 12 12 12 12 12 1 1 1 12 12 # AtmWgt= 12.0000000 12.0000000 12.0000000 12.0000000 12.0000000 1.0078250 1.0078250 1.0078250 12.0000000 12.0000000 # NucSpn= 0 0 0 0 0 1 1 1 0 0 # AtZEff= -3.6000000 -3.6000000 -3.6000000 -3.6000000 -3.6000000 -1.0000000 -1.0000000 -1.0000000 -3.6000000 -3.6000000 # NQMom= 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 # NMagM= 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 2.7928460 2.7928460 2.7928460 0.0000000 0.0000000 # ... with blank lines dividing blocks of ten, and Leave Link 101 at the end. # This is generally parsed before coordinates, so atomnos is not defined. # Note that in Gaussian03 the comments are not there yet and the labels are different. if line.strip() == "Isotopes and Nuclear Properties:": if not hasattr(self, "atommasses"): self.atommasses = [] line = next(inputfile) while line[1:16] != "Leave Link 101": if line[1:8] == "AtmWgt=": self.atommasses.extend(list(map(float, line.split()[1:]))) line = next(inputfile) # Extract the atomic numbers and coordinates of the atoms. if line.strip() == "Standard orientation:": self.updateprogress(inputfile, "Attributes", self.cupdate) # If this is a counterpoise calculation, this output means that # the supermolecule is now being considered, so we can set: self.counterpoise = 0 if not hasattr(self, "atomcoords"): self.atomcoords = [] self.skip_lines(inputfile, ['d', 'cols', 'cols', 'd']) atomnos = [] atomcoords = [] line = next(inputfile) while list(set(line.strip())) != ["-"]: broken = line.split() atomnos.append(int(broken[1])) atomcoords.append(list(map(float, broken[-3:]))) line = next(inputfile) self.atomcoords.append(atomcoords) self.set_attribute('natom', len(atomnos)) self.set_attribute('atomnos', atomnos) # This is a bit of a hack for regression Gaussian09/BH3_fragment_guess.pop_minimal.log # to skip output for all fragments, assuming the supermolecule is always printed first. # Eventually we want to make this more general, or even better parse the output for # all fragment, but that will happen in a newer version of cclib. if line[1:16] == "Fragment guess:" and getattr(self, 'nfragments', 0) > 1: if not "full" in line: inputfile.seek(0, 2) # Another hack for regression Gaussian03/ortho_prod_prod_freq.log, which is an ONIOM job. # Basically for now we stop parsing after the output for the real system, because # currently we don't support changes in system size or fragments in cclib. When we do, # we will want to parse the model systems, too, and that is what nmodels could track. if "ONIOM: generating point" in line and line.strip()[-13:] == 'model system.' and getattr(self, 'nmodels', 0) > 0: inputfile.seek(0, 2) # With the gfinput keyword, the atomic basis set functios are: # # AO basis set in the form of general basis input (Overlap normalization): # 1 0 # S 3 1.00 0.000000000000 # 0.7161683735D+02 0.1543289673D+00 # 0.1304509632D+02 0.5353281423D+00 # 0.3530512160D+01 0.4446345422D+00 # SP 3 1.00 0.000000000000 # 0.2941249355D+01 -0.9996722919D-01 0.1559162750D+00 # 0.6834830964D+00 0.3995128261D+00 0.6076837186D+00 # 0.2222899159D+00 0.7001154689D+00 0.3919573931D+00 # **** # 2 0 # S 3 1.00 0.000000000000 # 0.7161683735D+02 0.1543289673D+00 # ... # # The same is also printed when the gfprint keyword is used, but the # interstitial lines differ and there are no stars between atoms: # # AO basis set (Overlap normalization): # Atom C1 Shell 1 S 3 bf 1 - 1 0.509245180608 -2.664678875191 0.000000000000 # 0.7161683735D+02 0.1543289673D+00 # 0.1304509632D+02 0.5353281423D+00 # 0.3530512160D+01 0.4446345422D+00 # Atom C1 Shell 2 SP 3 bf 2 - 5 0.509245180608 -2.664678875191 0.000000000000 # 0.2941249355D+01 -0.9996722919D-01 0.1559162750D+00 # ... #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) if line[1:13] == "AO basis set" and not self.oniom: self.gbasis = [] # For counterpoise fragment calcualtions, skip these lines. if self.counterpoise != 0: return atom_line = inputfile.next() self.gfprint = atom_line.split()[0] == "Atom" self.gfinput = not self.gfprint # Note how the shell information is on a separate line for gfinput, # whereas for gfprint it is on the same line as atom information. if self.gfinput: shell_line = inputfile.next() shell = [] while len(self.gbasis) < self.natom: if self.gfprint: cols = atom_line.split() subshells = cols[4] ngauss = int(cols[5]) else: cols = shell_line.split() subshells = cols[0] ngauss = int(cols[1]) parameters = [] for ig in range(ngauss): line = inputfile.next() parameters.append(list(map(self.float, line.split()))) for iss, ss in enumerate(subshells): contractions = [] for param in parameters: exponent = param[0] coefficient = param[iss+1] contractions.append((exponent, coefficient)) subshell = (ss, contractions) shell.append(subshell) if self.gfprint: line = inputfile.next() if line.split()[0] == "Atom": atomnum = int(re.sub(r"\D", "", line.split()[1])) if atomnum == len(self.gbasis) + 2: self.gbasis.append(shell) shell = [] atom_line = line else: self.gbasis.append(shell) else: line = inputfile.next() if line.strip() == "****": self.gbasis.append(shell) shell = [] atom_line = inputfile.next() shell_line = inputfile.next() else: shell_line = line # Find the targets for SCF convergence (QM calcs). # Not for BOMD as targets are not available in the summary if not self.BOMD and line[1:44] == 'Requested convergence on RMS density matrix': if not hasattr(self, "scftargets"): self.scftargets = [] # The following can happen with ONIOM which are mixed SCF # and semi-empirical if type(self.scftargets) == type(numpy.array([])): self.scftargets = [] scftargets = [] # The RMS density matrix. scftargets.append(self.float(line.split('=')[1].split()[0])) line = next(inputfile) # The MAX density matrix. scftargets.append(self.float(line.strip().split('=')[1][:-1])) line = next(inputfile) # For G03, there's also the energy (not for G98). if line[1:10] == "Requested": scftargets.append(self.float(line.strip().split('=')[1][:-1])) self.scftargets.append(scftargets) # Extract SCF convergence information (QM calcs). if line[1:10] == 'Cycle 1': if not hasattr(self, "scfvalues"): self.scfvalues = [] scfvalues = [] line = next(inputfile) while line.find("SCF Done") == -1: self.updateprogress(inputfile, "QM convergence", self.fupdate) if line.find(' E=') == 0: self.logger.debug(line) # RMSDP=3.74D-06 MaxDP=7.27D-05 DE=-1.73D-07 OVMax= 3.67D-05 # or # RMSDP=1.13D-05 MaxDP=1.08D-04 OVMax= 1.66D-04 if line.find(" RMSDP") == 0: # Fields of interest: # RMSDP # MaxDP # (DE) -> Only add the energy if it's a target criteria matches = self.re_scinot.findall(line) matches = { match[0]: self.float(match[1]) for match in matches } scfvalues_step = [ matches.get('RMSDP', numpy.nan), matches.get('MaxDP', numpy.nan) ] if len(self.scftargets[0]) == 3: scfvalues_step.append(matches.get('DE', numpy.nan)) scfvalues.append(scfvalues_step) try: line = next(inputfile) # May be interupted by EOF. except StopIteration: self.logger.warning('File terminated before end of last SCF!') break self.scfvalues.append(scfvalues) # Extract SCF convergence information (AM1, INDO and other semi-empirical calcs). # The output (for AM1) looks like this: # Ext34=T Pulay=F Camp-King=F BShift= 0.00D+00 # It= 1 PL= 0.103D+01 DiagD=T ESCF= 31.564733 Diff= 0.272D+02 RMSDP= 0.152D+00. # It= 2 PL= 0.114D+00 DiagD=T ESCF= 7.265370 Diff=-0.243D+02 RMSDP= 0.589D-02. # ... # It= 11 PL= 0.184D-04 DiagD=F ESCF= 4.687669 Diff= 0.260D-05 RMSDP= 0.134D-05. # It= 12 PL= 0.105D-04 DiagD=F ESCF= 4.687669 Diff=-0.686D-07 RMSDP= 0.215D-05. # 4-point extrapolation. # It= 13 PL= 0.110D-05 DiagD=F ESCF= 4.687669 Diff=-0.111D-06 RMSDP= 0.653D-07. # Energy= 0.172272018655 NIter= 14. if line[1:4] == 'It=': scftargets = numpy.array([1E-7], "d") # This is the target value for the rms scfvalues = [[]] while line.find(" Energy") == -1: self.updateprogress(inputfile, "AM1 Convergence") if line[1:4] == "It=": parts = line.strip().split() scfvalues[0].append(self.float(parts[-1][:-1])) line = next(inputfile) # If an AM1 or INDO guess is used (Guess=INDO in the input, for example), # this will be printed after a single iteration, so that is the line # that should trigger a break from this loop. At least that's what we see # for regression Gaussian/Gaussian09/guessIndo_modified_ALT.out if line[:14] == " Initial guess": break # Attach the attributes to the object Only after the energy is found . if line.find(" Energy") == 0: self.scftargets = scftargets self.scfvalues = scfvalues # Note: this needs to follow the section where 'SCF Done' is used # to terminate a loop when extracting SCF convergence information. if not self.BOMD and line[1:9] == 'SCF Done': t1 = line.split()[2] if t1 == 'E(RHF)': self.metadata["methods"].append("HF") else: self.metadata["methods"].append("DFT") self.metadata["functional"] = t1[t1.index("(") + 2:t1.rindex(")")] if not hasattr(self, "scfenergies"): self.scfenergies = [] self.scfenergies.append(utils.convertor(self.float(line.split()[4]), "hartree", "eV")) # gmagoon 5/27/09: added scfenergies reading for PM3 case # Example line: " Energy= -0.077520562724 NIter= 14." # See regression Gaussian03/QVGXLLKOCUKJST-UHFFFAOYAJmult3Fixed.out if line[1:8] == 'Energy=': if not hasattr(self, "scfenergies"): self.scfenergies = [] self.scfenergies.append(utils.convertor(self.float(line.split()[1]), "hartree", "eV")) # Total energies after Moller-Plesset corrections. # Second order correction is always first, so its first occurance # triggers creation of mpenergies (list of lists of energies). # Further MP2 corrections are appended as found. # # Example MP2 output line: # E2 = -0.9505918144D+00 EUMP2 = -0.28670924198852D+03 # Warning! this output line is subtly different for MP3/4/5 runs if "EUMP2" in line[27:34]: self.metadata["methods"].append("MP2") if not hasattr(self, "mpenergies"): self.mpenergies = [] self.mpenergies.append([]) mp2energy = self.float(line.split("=")[2]) self.mpenergies[-1].append(utils.convertor(mp2energy, "hartree", "eV")) # Example MP3 output line: # E3= -0.10518801D-01 EUMP3= -0.75012800924D+02 if line[34:39] == "EUMP3": self.metadata["methods"].append("MP3") mp3energy = self.float(line.split("=")[2]) self.mpenergies[-1].append(utils.convertor(mp3energy, "hartree", "eV")) # Example MP4 output lines: # E4(DQ)= -0.31002157D-02 UMP4(DQ)= -0.75015901139D+02 # E4(SDQ)= -0.32127241D-02 UMP4(SDQ)= -0.75016013648D+02 # E4(SDTQ)= -0.32671209D-02 UMP4(SDTQ)= -0.75016068045D+02 # Energy for most substitutions is used only (SDTQ by default) if line[34:42] == "UMP4(DQ)": self.metadata["methods"].append("MP4") mp4energy = self.float(line.split("=")[2]) line = next(inputfile) if line[34:43] == "UMP4(SDQ)": mp4energy = self.float(line.split("=")[2]) line = next(inputfile) if line[34:44] == "UMP4(SDTQ)": mp4energy = self.float(line.split("=")[2]) self.mpenergies[-1].append(utils.convertor(mp4energy, "hartree", "eV")) # Example MP5 output line: # DEMP5 = -0.11048812312D-02 MP5 = -0.75017172926D+02 if line[29:32] == "MP5": self.metadata["methods"].append("MP5") mp5energy = self.float(line.split("=")[2]) self.mpenergies[-1].append(utils.convertor(mp5energy, "hartree", "eV")) # Total energies after Coupled Cluster corrections. # Second order MBPT energies (MP2) are also calculated for these runs, # but the output is the same as when parsing for mpenergies. # Read the consecutive correlated energies # but append only the last one to ccenergies. # Only the highest level energy is appended - ex. CCSD(T), not CCSD. if line[1:10] == "DE(Corr)=" and line[27:35] == "E(CORR)=": self.metadata["methods"].append("CCSD") self.ccenergy = self.float(line.split()[3]) if line[1:10] == "T5(CCSD)=": line = next(inputfile) if line[1:9] == "CCSD(T)=": self.metadata["methods"].append("CCSD-T") self.ccenergy = self.float(line.split()[1]) if line[12:53] == "Population analysis using the SCF density": if hasattr(self, "ccenergy"): if not hasattr(self, "ccenergies"): self.ccenergies = [] self.ccenergies.append(utils.convertor(self.ccenergy, "hartree", "eV")) del self.ccenergy # Find step number for current optimization/IRC # Matches "Step number 123", "Pt XX Step number 123" and "PtXXX Step number 123" if " Step number" in line: step = int(line.split()[line.split().index('Step') + 2]) if not hasattr(self, "optstatus"): self.optstatus = [] self.optstatus.append(data.ccData.OPT_UNKNOWN) if step == 1: self.optstatus[-1] += data.ccData.OPT_NEW # Geometry convergence information. if line[49:59] == 'Converged?': if not hasattr(self, "geotargets"): self.geovalues = [] self.geotargets = numpy.array([0.0, 0.0, 0.0, 0.0], "d") newlist = [0]*4 for i in range(4): line = next(inputfile) self.logger.debug(line) parts = line.split() try: value = self.float(parts[2]) except ValueError: self.logger.error("Problem parsing the value for geometry optimisation: %s is not a number." % parts[2]) else: newlist[i] = value self.geotargets[i] = self.float(parts[3]) self.geovalues.append(newlist) # Gradients. # Read in the cartesian energy gradients (forces) from a block like this: # ------------------------------------------------------------------- # Center Atomic Forces (Hartrees/Bohr) # Number Number X Y Z # ------------------------------------------------------------------- # 1 1 -0.012534744 -0.021754635 -0.008346094 # 2 6 0.018984731 0.032948887 -0.038003451 # 3 1 -0.002133484 -0.006226040 0.023174772 # 4 1 -0.004316502 -0.004968213 0.023174772 # -2 -0.001830728 -0.000743108 -0.000196625 # ------------------------------------------------------------------ # # The "-2" line is for a dummy atom # # Then optimization is done in internal coordinates, Gaussian also # print the forces in internal coordinates, which can be produced from # the above. This block looks like this: # Variable Old X -DE/DX Delta X Delta X Delta X New X # (Linear) (Quad) (Total) # ch 2.05980 0.01260 0.00000 0.01134 0.01134 2.07114 # hch 1.75406 0.09547 0.00000 0.24861 0.24861 2.00267 # hchh 2.09614 0.01261 0.00000 0.16875 0.16875 2.26489 # Item Value Threshold Converged? # We could get the gradients in BOMD, but it is more complex because # they are not in the summary, and they are not as relevant as for # an optimization if not self.BOMD and line[37:43] == "Forces": if not hasattr(self, "grads"): self.grads = [] self.skip_lines(inputfile, ['header', 'd']) forces = [] line = next(inputfile) while list(set(line.strip())) != ['-']: tmpforces = [] for N in range(3): # Fx, Fy, Fz force = line[23+N*15:38+N*15] if force.startswith("*"): force = "NaN" tmpforces.append(float(force)) forces.append(tmpforces) line = next(inputfile) self.grads.append(forces) #Extract PES scan data #Summary of the potential surface scan: # N A SCF #---- --------- ----------- # 1 109.0000 -76.43373 # 2 119.0000 -76.43011 # 3 129.0000 -76.42311 # 4 139.0000 -76.41398 # 5 149.0000 -76.40420 # 6 159.0000 -76.39541 # 7 169.0000 -76.38916 # 8 179.0000 -76.38664 # 9 189.0000 -76.38833 # 10 199.0000 -76.39391 # 11 209.0000 -76.40231 #---- --------- ----------- if "Summary of the potential surface scan:" in line: scanenergies = [] scanparm = [] colmnames = next(inputfile) hyphens = next(inputfile) line = next(inputfile) while line != hyphens: broken = line.split() scanenergies.append(float(broken[-1])) scanparm.append(map(float, broken[1:-1])) line = next(inputfile) if not hasattr(self, "scanenergies"): self.scanenergies = [] self.scanenergies = scanenergies if not hasattr(self, "scanparm"): self.scanparm = [] self.scanparm = scanparm if not hasattr(self, "scannames"): self.scannames = colmnames.split()[1:-1] # Orbital symmetries. if line[1:20] == 'Orbital symmetries:' and not hasattr(self, "mosyms"): # For counterpoise fragments, skip these lines. if self.counterpoise != 0: return self.updateprogress(inputfile, "MO Symmetries", self.fupdate) self.mosyms = [[]] line = next(inputfile) unres = False if line.find("Alpha Orbitals") == 1: unres = True line = next(inputfile) i = 0 while len(line) > 18 and line[17] == '(': if line.find('Virtual') >= 0: self.homos = [i - 1] parts = line[17:].split() for x in parts: self.mosyms[0].append(self.normalisesym(x.strip('()'))) i += 1 line = next(inputfile) if unres: line = next(inputfile) # Repeat with beta orbital information i = 0 self.mosyms.append([]) while len(line) > 18 and line[17] == '(': if line.find('Virtual') >= 0: # Here we consider beta # If there was also an alpha virtual orbital, # we will store two indices in the array # Otherwise there is no alpha virtual orbital, # only beta virtual orbitals, and we initialize # the array with one element. See the regression # QVGXLLKOCUKJST-UHFFFAOYAJmult3Fixed.out # donated by Gregory Magoon (gmagoon). if (hasattr(self, "homos")): # Extend the array to two elements # 'HOMO' indexes the HOMO in the arrays self.homos.append(i-1) else: # 'HOMO' indexes the HOMO in the arrays self.homos = [i - 1] parts = line[17:].split() for x in parts: self.mosyms[1].append(self.normalisesym(x.strip('()'))) i += 1 line = next(inputfile) # Some calculations won't explicitely print the number of basis sets used, # and will occasionally drop some without warning. We can infer the number, # however, from the MO symmetries printed here. Specifically, this fixes # regression Gaussian/Gaussian09/dvb_sp_terse.log (#23 on github). self.set_attribute('nmo', len(self.mosyms[-1])) # Alpha/Beta electron eigenvalues. if line[1:6] == "Alpha" and line.find("eigenvalues") >= 0: # For counterpoise fragments, skip these lines. if self.counterpoise != 0: return # For ONIOM calcs, ignore this section in order to bypass assertion failure. if self.oniom: return self.updateprogress(inputfile, "Eigenvalues", self.fupdate) self.moenergies = [[]] HOMO = -2 while line.find('Alpha') == 1: if line.split()[1] == "virt." and HOMO == -2: # If there aren't any symmetries, this is a good way to find the HOMO. HOMO = len(self.moenergies[0])-1 self.homos = [HOMO] # Convert to floats and append to moenergies, but sometimes Gaussian # doesn't print correctly so test for ValueError (bug 1756789). part = line[28:] i = 0 while i*10+4 < len(part): s = part[i*10:(i+1)*10] try: x = self.float(s) except ValueError: x = numpy.nan self.moenergies[0].append(utils.convertor(x, "hartree", "eV")) i += 1 line = next(inputfile) # If, at this point, self.homos is unset, then there were not # any alpha virtual orbitals if not hasattr(self, "homos"): HOMO = len(self.moenergies[0])-1 self.homos = [HOMO] if line.find('Beta') == 2: self.moenergies.append([]) HOMO = -2 while line.find('Beta') == 2: if line.split()[1] == "virt." and HOMO == -2: # If there aren't any symmetries, this is a good way to find the HOMO. # Also, check for consistency if homos was already parsed. HOMO = len(self.moenergies[1])-1 self.homos.append(HOMO) part = line[28:] i = 0 while i*10+4 < len(part): x = part[i*10:(i+1)*10] self.moenergies[1].append(utils.convertor(self.float(x), "hartree", "eV")) i += 1 line = next(inputfile) self.moenergies = [numpy.array(x, "d") for x in self.moenergies] # Start of the IR/Raman frequency section. # Caution is advised here, as additional frequency blocks # can be printed by Gaussian (with slightly different formats), # often doubling the information printed. # See, for a non-standard exmaple, regression Gaussian98/test_H2.log # If either the Gaussian freq=hpmodes keyword or IOP(7/33=1) is used, # an extra frequency block with higher-precision vibdisps is # printed before the normal frequency block. # Note that the code parses only the vibsyms and vibdisps # from the high-precision block, but parses vibsyms, vibfreqs, # vibramans and vibirs from the normal block. vibsyms parsed # from the high-precision block are discarded and replaced by those # from the normal block while the high-precision vibdisps, if present, # are used to overwrite default-precision vibdisps at the end of the parse. if line[1:14] == "Harmonic freq": # This matches in both freq block types self.updateprogress(inputfile, "Frequency Information", self.fupdate) # The whole block should not have any blank lines. while line.strip() != "": # The line with indices if line[1:15].strip() == "" and line[15:60].split()[0].isdigit(): freqbase = int(line[15:60].split()[0]) if freqbase == 1 and hasattr(self, 'vibsyms'): # we are coming accross duplicated information. # We might be be parsing a default-precision block having # already parsed (only) vibsyms and displacements from # the high-precision block, or might be encountering # a second low-precision block (see e.g. 25DMF_HRANH.log # regression). self.vibsyms = [] if hasattr(self, "vibirs"): self.vibirs = [] if hasattr(self, 'vibfreqs'): self.vibfreqs = [] if hasattr(self, 'vibramans'): self.vibramans = [] if hasattr(self, 'vibdisps'): self.vibdisps = [] # Lines with symmetries and symm. indices begin with whitespace. if line[1:15].strip() == "" and not line[15:60].split()[0].isdigit(): if not hasattr(self, 'vibsyms'): self.vibsyms = [] syms = line.split() self.vibsyms.extend(syms) if line[1:15] == "Frequencies --": # note: matches low-precision block, and if not hasattr(self, 'vibfreqs'): self.vibfreqs = [] freqs = [self.float(f) for f in line[15:].split()] self.vibfreqs.extend(freqs) if line[1:15] == "IR Inten --": # note: matches only low-precision block if not hasattr(self, 'vibirs'): self.vibirs = [] irs = [] for ir in line[15:].split(): try: irs.append(self.float(ir)) except ValueError: irs.append(self.float('nan')) self.vibirs.extend(irs) if line[1:15] == "Raman Activ --": # note: matches only low-precision block if not hasattr(self, 'vibramans'): self.vibramans = [] ramans = [] for raman in line[15:].split(): try: ramans.append(self.float(raman)) except ValueError: ramans.append(self.float('nan')) self.vibramans.extend(ramans) # Block with (default-precision) displacements should start with this. # 1 2 3 # A A A # Frequencies -- 370.7936 370.7987 618.0103 # Red. masses -- 2.3022 2.3023 1.9355 # Frc consts -- 0.1865 0.1865 0.4355 # IR Inten -- 0.0000 0.0000 0.0000 # Atom AN X Y Z X Y Z X Y Z # 1 6 0.00 0.00 -0.04 0.00 0.00 0.19 0.00 0.00 0.12 # 2 6 0.00 0.00 0.19 0.00 0.00 -0.06 0.00 0.00 -0.12 if line.strip().split()[0:3] == ["Atom", "AN", "X"]: if not hasattr(self, 'vibdisps'): self.vibdisps = [] disps = [] for n in range(self.natom): line = next(inputfile) numbers = [float(s) for s in line[10:].split()] N = len(numbers) // 3 if not disps: for n in range(N): disps.append([]) for n in range(N): disps[n].append(numbers[3*n:3*n+3]) self.vibdisps.extend(disps) # Block with high-precision (freq=hpmodes) displacements should start with this. # 1 2 3 4 5 # A A A A A # Frequencies --- 370.7936 370.7987 618.0103 647.7864 647.7895 # Reduced masses --- 2.3022 2.3023 1.9355 6.4600 6.4600 # Force constants --- 0.1865 0.1865 0.4355 1.5971 1.5972 # IR Intensities --- 0.0000 0.0000 0.0000 0.0000 0.0000 # Coord Atom Element: # 1 1 6 0.00000 0.00000 0.00000 -0.18677 0.05592 # 2 1 6 0.00000 0.00000 0.00000 0.28440 0.21550 # 3 1 6 -0.04497 0.19296 0.11859 0.00000 0.00000 # 1 2 6 0.00000 0.00000 0.00000 0.03243 0.37351 # 2 2 6 0.00000 0.00000 0.00000 0.14503 -0.06117 # 3 2 6 0.18959 -0.05753 -0.11859 0.00000 0.00000 if line.strip().split()[0:3] == ["Coord", "Atom", "Element:"]: # Wait until very end of parsing to assign vibdispshp to self.vibdisps # as otherwise the higher precision displacements will be overwritten # by low precision displacements which are printed further down file if not hasattr(self, 'vibdispshp'): self.vibdispshp = [] disps = [] for n in range(3*self.natom): line = next(inputfile) numbers = [float(s) for s in line[16:].split()] atomindex = int(line[4:10])-1 # atom index, starting at zero numbermodes = len(numbers) if not disps: for mode in range(numbermodes): # For each mode, make list of list [atom][coord_index] disps.append([[] for x in range(0, self.natom)]) for mode in range(numbermodes): disps[mode][atomindex].append(numbers[mode]) self.vibdispshp.extend(disps) line = next(inputfile) # Electronic transitions. if line[1:14] == "Excited State": if not hasattr(self, "etenergies"): self.etenergies = [] self.etoscs = [] self.etsyms = [] self.etsecs = [] # Need to deal with lines like: # (restricted calc) # Excited State 1: Singlet-BU 5.3351 eV 232.39 nm f=0.1695 # (unrestricted calc) (first excited state is 2!) # Excited State 2: ?Spin -A 0.1222 eV 10148.75 nm f=0.0000 # (Gaussian 09 ZINDO) # Excited State 1: Singlet-?Sym 2.5938 eV 478.01 nm f=0.0000 =0.000 p = re.compile(":(?P.*?)(?P-?\d*\.\d*) eV") groups = p.search(line).groups() self.etenergies.append(utils.convertor(self.float(groups[1]), "eV", "wavenumber")) self.etoscs.append(self.float(line.split("f=")[-1].split()[0])) self.etsyms.append(groups[0].strip()) line = next(inputfile) p = re.compile("(\d+)") CIScontrib = [] while line.find(" ->") >= 0: # This is a contribution to the transition parts = line.split("->") self.logger.debug(parts) # Has to deal with lines like: # 32 -> 38 0.04990 # 35A -> 45A 0.01921 frommoindex = 0 # For restricted or alpha unrestricted fromMO = parts[0].strip() if fromMO[-1] == "B": frommoindex = 1 # For beta unrestricted fromMO = int(p.match(fromMO).group())-1 # subtract 1 so that it is an index into moenergies t = parts[1].split() tomoindex = 0 toMO = t[0] if toMO[-1] == "B": tomoindex = 1 toMO = int(p.match(toMO).group())-1 # subtract 1 so that it is an index into moenergies percent = self.float(t[1]) # For restricted calculations, the percentage will be corrected # after parsing (see after_parsing() above). CIScontrib.append([(fromMO, frommoindex), (toMO, tomoindex), percent]) line = next(inputfile) self.etsecs.append(CIScontrib) # Circular dichroism data (different for G03 vs G09) # # G03 # # ## <0|r|b> * (Au), Rotatory Strengths (R) in # ## cgs (10**-40 erg-esu-cm/Gauss) # ## state X Y Z R(length) # ## 1 0.0006 0.0096 -0.0082 -0.4568 # ## 2 0.0251 -0.0025 0.0002 -5.3846 # ## 3 0.0168 0.4204 -0.3707 -15.6580 # ## 4 0.0721 0.9196 -0.9775 -3.3553 # # G09 # # ## 1/2[<0|r|b>* + (<0|rxdel|b>*)*] # ## Rotatory Strengths (R) in cgs (10**-40 erg-esu-cm/Gauss) # ## state XX YY ZZ R(length) R(au) # ## 1 -0.3893 -6.7546 5.7736 -0.4568 -0.0010 # ## 2 -17.7437 1.7335 -0.1435 -5.3845 -0.0114 # ## 3 -11.8655 -297.2604 262.1519 -15.6580 -0.0332 if line[1:52] == "<0|r|b> * (Au), Rotatory Strengths (R)" or \ line[1:50] == "1/2[<0|r|b>* + (<0|rxdel|b>*)*]": self.etrotats = [] self.skip_lines(inputfile, ['units']) headers = next(inputfile) Ncolms = len(headers.split()) line = next(inputfile) parts = line.strip().split() while len(parts) == Ncolms: try: R = self.float(parts[4]) except ValueError: # nan or -nan if there is no first excited state # (for unrestricted calculations) pass else: self.etrotats.append(R) line = next(inputfile) temp = line.strip().split() parts = line.strip().split() self.etrotats = numpy.array(self.etrotats, "d") # Number of basis sets functions. # Has to deal with lines like: # NBasis = 434 NAE= 97 NBE= 97 NFC= 34 NFV= 0 # and... # NBasis = 148 MinDer = 0 MaxDer = 0 # Although the former is in every file, it doesn't occur before # the overlap matrix is printed. if line[1:7] == "NBasis" or line[4:10] == "NBasis": # For counterpoise fragment, skip these lines. if self.counterpoise != 0: return # For ONIOM calcs, ignore this section in order to bypass assertion failure. if self.oniom: return # If nbasis was already parsed, check if it changed. If it did, issue a warning. # In the future, we will probably want to have nbasis, as well as nmo below, # as a list so that we don't need to pick one value when it changes. nbasis = int(line.split('=')[1].split()[0]) if hasattr(self, "nbasis"): try: assert nbasis == self.nbasis except AssertionError: self.logger.warning("Number of basis functions (nbasis) has changed from %i to %i" % (self.nbasis, nbasis)) self.nbasis = nbasis # Number of linearly-independent basis sets. if line[1:7] == "NBsUse": # For counterpoise fragment, skip these lines. if self.counterpoise != 0: return # For ONIOM calcs, ignore this section in order to bypass assertion failure. if self.oniom: return nmo = int(line.split('=')[1].split()[0]) self.set_attribute('nmo', nmo) # For AM1 calculations, set nbasis by a second method, # as nmo may not always be explicitly stated. if line[7:22] == "basis functions, ": nbasis = int(line.split()[0]) self.set_attribute('nbasis', nbasis) # Molecular orbital overlap matrix. # Has to deal with lines such as: # *** Overlap *** # ****** Overlap ****** # Note that Gaussian sometimes drops basis functions, # causing the overlap matrix as parsed below to not be # symmetric (which is a problem for population analyses, etc.) if line[1:4] == "***" and (line[5:12] == "Overlap" or line[8:15] == "Overlap"): # Ensure that this is the main calc and not a fragment if self.counterpoise != 0: return self.aooverlaps = numpy.zeros((self.nbasis, self.nbasis), "d") # Overlap integrals for basis fn#1 are in aooverlaps[0] base = 0 colmNames = next(inputfile) while base < self.nbasis: self.updateprogress(inputfile, "Overlap", self.fupdate) for i in range(self.nbasis-base): # Fewer lines this time line = next(inputfile) parts = line.split() for j in range(len(parts)-1): # Some lines are longer than others k = float(parts[j+1].replace("D", "E")) self.aooverlaps[base+j, i+base] = k self.aooverlaps[i+base, base+j] = k base += 5 colmNames = next(inputfile) self.aooverlaps = numpy.array(self.aooverlaps, "d") # Molecular orbital coefficients (mocoeffs). # Essentially only produced for SCF calculations. # This is also the place where aonames and atombasis are parsed. if line[5:35] == "Molecular Orbital Coefficients" or line[5:41] == "Alpha Molecular Orbital Coefficients" or line[5:40] == "Beta Molecular Orbital Coefficients": # If counterpoise fragment, return without parsing orbital info if self.counterpoise != 0: return # Skip this for ONIOM calcs if self.oniom: return if line[5:40] == "Beta Molecular Orbital Coefficients": beta = True if self.popregular: return # This was continue before refactoring the parsers. #continue # Not going to extract mocoeffs # Need to add an extra array to self.mocoeffs self.mocoeffs.append(numpy.zeros((self.nmo, self.nbasis), "d")) else: beta = False self.aonames = [] self.atombasis = [] mocoeffs = [numpy.zeros((self.nmo, self.nbasis), "d")] base = 0 self.popregular = False for base in range(0, self.nmo, 5): self.updateprogress(inputfile, "Coefficients", self.fupdate) colmNames = next(inputfile) if not colmNames.split(): self.logger.warning("Molecular coefficients header found but no coefficients.") break if base == 0 and int(colmNames.split()[0]) != 1: # Implies that this is a POP=REGULAR calculation # and so, only aonames (not mocoeffs) will be extracted self.popregular = True symmetries = next(inputfile) eigenvalues = next(inputfile) for i in range(self.nbasis): line = next(inputfile) if i == 0: # Find location of the start of the basis function name start_of_basis_fn_name = line.find(line.split()[3]) - 1 if base == 0 and not beta: # Just do this the first time 'round parts = line[:start_of_basis_fn_name].split() if len(parts) > 1: # New atom if i > 0: self.atombasis.append(atombasis) atombasis = [] atomname = "%s%s" % (parts[2], parts[1]) orbital = line[start_of_basis_fn_name:20].strip() self.aonames.append("%s_%s" % (atomname, orbital)) atombasis.append(i) part = line[21:].replace("D", "E").rstrip() temp = [] for j in range(0, len(part), 10): temp.append(float(part[j:j+10])) if beta: self.mocoeffs[1][base:base + len(part) // 10, i] = temp else: mocoeffs[0][base:base + len(part) // 10, i] = temp if base == 0 and not beta: # Do the last update of atombasis self.atombasis.append(atombasis) if self.popregular: # We now have aonames, so no need to continue break if not self.popregular and not beta: self.mocoeffs = mocoeffs # Natural orbital coefficients (nocoeffs) and occupation numbers (nooccnos), # which are respectively define the eigenvectors and eigenvalues of the # diagonalized one-electron density matrix. These orbitals are formed after # configuration interaction (CI) calculations, but not only. Similarly to mocoeffs, # we can parse and check aonames and atombasis here. # # Natural Orbital Coefficients: # 1 2 3 4 5 # Eigenvalues -- 2.01580 2.00363 2.00000 2.00000 1.00000 # 1 1 O 1S 0.00000 -0.15731 -0.28062 0.97330 0.00000 # 2 2S 0.00000 0.75440 0.57746 0.07245 0.00000 # ... # def natural_orbital_single_spin_parsing(inputfile, updateprogress_title): coeffs = numpy.zeros((self.nmo, self.nbasis), "d") occnos = [] aonames = [] atombasis = [] for base in range(0, self.nmo, 5): self.updateprogress(inputfile, updateprogress_title, self.fupdate) colmNames = next(inputfile) eigenvalues = next(inputfile) occnos.extend(map(float, eigenvalues.split()[2:])) for i in range(self.nbasis): line = next(inputfile) # Just do this the first time 'round. if base == 0: # Changed below from :12 to :11 to deal with Elmar Neumann's example. parts = line[:11].split() # New atom. if len(parts) > 1: if i > 0: atombasis.append(basisonatom) basisonatom = [] atomname = "%s%s" % (parts[2], parts[1]) orbital = line[11:20].strip() aonames.append("%s_%s" % (atomname, orbital)) basisonatom.append(i) part = line[21:].replace("D", "E").rstrip() temp = [] for j in range(0, len(part), 10): temp.append(float(part[j:j+10])) coeffs[base:base + len(part) // 10, i] = temp # Do the last update of atombasis. if base == 0: atombasis.append(basisonatom) return occnos, coeffs, aonames, atombasis if line[5:33] == "Natural Orbital Coefficients": updateprogress_title = "Natural orbitals" nooccnos, nocoeffs, aonames, atombasis = natural_orbital_single_spin_parsing(inputfile, updateprogress_title) self.set_attribute("nocoeffs", nocoeffs) self.set_attribute("nooccnos", nooccnos) self.set_attribute("atombasis", atombasis) self.set_attribute("aonames", aonames) # Natural spin orbital coefficients (nsocoeffs) and occupation numbers (nsooccnos) # Parsed attributes are similar to the natural orbitals above except # the natural spin orbitals and occupation numbers are the eigenvalues # and eigenvectors of the one particles spin density matrices # Alpha Natural Orbital Coefficients: # 1 2 3 4 5 # Eigenvalues -- 1.00000 1.00000 0.99615 0.99320 0.99107 # 1 1 O 1S 0.70425 0.70600 -0.16844 -0.14996 -0.00000 # 2 2S 0.01499 0.01209 0.36089 0.34940 -0.00000 # ... # Beta Natural Orbital Coefficients: # 1 2 3 4 5 # Eigenvalues -- 1.00000 1.00000 0.99429 0.98790 0.98506 # 1 1 O 1S 0.70822 0.70798 -0.15316 -0.13458 0.00465 # 2 2S 0.00521 0.00532 0.33837 0.33189 -0.01301 # 3 3S -0.02542 -0.00841 0.28649 0.53224 0.18902 # ... if line[5:39] == "Alpha Natural Orbital Coefficients": updateprogress_title = "Natural Spin orbitals (alpha)" nsooccnos, nsocoeffs, aonames, atombasis = natural_orbital_single_spin_parsing(inputfile, updateprogress_title) if self.unified_no_nso: self.append_attribute("nocoeffs", nsocoeffs) self.append_attribute("nooccnos", nsooccnos) else: self.append_attribute("nsocoeffs", nsocoeffs) self.append_attribute("nsooccnos", nsooccnos) self.set_attribute("atombasis", atombasis) self.set_attribute("aonames", aonames) if line[5:38] == "Beta Natural Orbital Coefficients": updateprogress_title = "Natural Spin orbitals (beta)" nsooccnos, nsocoeffs, aonames, atombasis = natural_orbital_single_spin_parsing(inputfile, updateprogress_title) if self.unified_no_nso: self.append_attribute("nocoeffs", nsocoeffs) self.append_attribute("nooccnos", nsooccnos) else: self.append_attribute("nsocoeffs", nsocoeffs) self.append_attribute("nsooccnos", nsooccnos) self.set_attribute("atombasis", atombasis) self.set_attribute("aonames", aonames) # For FREQ=Anharm, extract anharmonicity constants if line[1:40] == "X matrix of Anharmonic Constants (cm-1)": Nvibs = len(self.vibfreqs) self.vibanharms = numpy.zeros((Nvibs, Nvibs), "d") base = 0 colmNames = next(inputfile) while base < Nvibs: for i in range(Nvibs-base): # Fewer lines this time line = next(inputfile) parts = line.split() for j in range(len(parts)-1): # Some lines are longer than others k = float(parts[j+1].replace("D", "E")) self.vibanharms[base+j, i+base] = k self.vibanharms[i+base, base+j] = k base += 5 colmNames = next(inputfile) # Pseudopotential charges. if line.find("Pseudopotential Parameters") > -1: self.skip_lines(inputfile, ['e', 'label1', 'label2', 'e']) line = next(inputfile) if line.find("Centers:") < 0: return # This was continue before parser refactoring. # continue # Needs to handle code like the following: # # Center Atomic Valence Angular Power Coordinates # Number Number Electrons Momentum of R Exponent Coefficient X Y Z # =================================================================================================================================== # Centers: 1 # Centers: 16 # Centers: 21 24 # Centers: 99100101102 # 1 44 16 -4.012684 -0.696698 0.006750 # F and up # 0 554.3796303 -0.05152700 centers = [] while line.find("Centers:") >= 0: temp = line[10:] for i in range(0, len(temp)-3, 3): centers.append(int(temp[i:i+3])) line = next(inputfile) centers.sort() # Not always in increasing order self.coreelectrons = numpy.zeros(self.natom, "i") for center in centers: front = line[:10].strip() while not (front and int(front) == center): line = next(inputfile) front = line[:10].strip() info = line.split() self.coreelectrons[center-1] = int(info[1]) - int(info[2]) line = next(inputfile) # This will be printed for counterpoise calcualtions only. # To prevent crashing, we need to know which fragment is being considered. # Other information is also printed in lines that start like this. if line[1:14] == 'Counterpoise:': if line[42:50] == "fragment": self.counterpoise = int(line[51:54]) # This will be printed only during ONIOM calcs; use it to set a flag # that will allow assertion failures to be bypassed in the code. if line[1:7] == "ONIOM:": self.oniom = True # This will be printed only during BOMD calcs; if line.startswith(" INPUT DATA FOR L118"): self.BOMD = True # Atomic charges are straightforward to parse, although the header # has changed over time somewhat. # # Mulliken charges: # 1 # 1 C -0.004513 # 2 C -0.077156 # ... # Sum of Mulliken charges = 0.00000 # Mulliken charges with hydrogens summed into heavy atoms: # 1 # 1 C -0.004513 # 2 C 0.002063 # ... # if line[1:25] == "Mulliken atomic charges:" or line[1:18] == "Mulliken charges:" or \ line[1:23] == "Lowdin Atomic Charges:" or line[1:16] == "Lowdin charges:": if not hasattr(self, "atomcharges"): self.atomcharges = {} ones = next(inputfile) charges = [] nline = next(inputfile) while not "Sum of" in nline: charges.append(float(nline.split()[2])) nline = next(inputfile) if "Mulliken" in line: self.atomcharges["mulliken"] = charges else: self.atomcharges["lowdin"] = charges if line.strip() == "Natural Population": if not hasattr(self, 'atomcharges'): self.atomcharges = {} line1 = next(inputfile) line2 = next(inputfile) if line1.split()[0] == 'Natural' and line2.split()[2] == 'Charge': dashes = next(inputfile) charges = [] for i in range(self.natom): nline = next(inputfile) charges.append(float(nline.split()[2])) self.atomcharges["natural"] = charges #Extract Thermochemistry #Temperature 298.150 Kelvin. Pressure 1.00000 Atm. #Zero-point correction= 0.342233 (Hartree/ #Thermal correction to Energy= 0. #Thermal correction to Enthalpy= 0. #Thermal correction to Gibbs Free Energy= 0.302940 #Sum of electronic and zero-point Energies= -563.649744 #Sum of electronic and thermal Energies= -563.636699 #Sum of electronic and thermal Enthalpies= -563.635755 #Sum of electronic and thermal Free Energies= -563.689037 if "Sum of electronic and thermal Enthalpies" in line: self.set_attribute('enthalpy', float(line.split()[6])) if "Sum of electronic and thermal Free Energies=" in line: self.set_attribute('freeenergy', float(line.split()[7])) if line[1:13] == "Temperature ": self.set_attribute('temperature', float(line.split()[1])) self.set_attribute('pressure', float(line.split()[4])) # Static polarizability (from `polar`), lower triangular # matrix. if line[1:26] == "SCF Polarizability for W=": self.hp_polarizabilities = True if not hasattr(self, 'polarizabilities'): self.polarizabilities = [] polarizability = numpy.zeros(shape=(3, 3)) self.skip_line(inputfile, 'directions') for i in range(3): line = next(inputfile) polarizability[i, :i+1] = [self.float(x) for x in line.split()[1:]] polarizability = utils.symmetrize(polarizability, use_triangle='lower') self.polarizabilities.append(polarizability) # Static polarizability (from `freq`), lower triangular matrix. if line[1:16] == "Polarizability=": self.hp_polarizabilities = True if not hasattr(self, 'polarizabilities'): self.polarizabilities = [] polarizability = numpy.zeros(shape=(3, 3)) polarizability_list = [] polarizability_list.extend([line[16:31], line[31:46], line[46:61]]) line = next(inputfile) polarizability_list.extend([line[16:31], line[31:46], line[46:61]]) indices = numpy.tril_indices(3) polarizability[indices] = [self.float(x) for x in polarizability_list] polarizability = utils.symmetrize(polarizability, use_triangle='lower') self.polarizabilities.append(polarizability) # Static polarizability, compressed into a single line from # terse printing. # Order is XX, YX, YY, ZX, ZY, ZZ (lower triangle). if line[2:23] == "Exact polarizability:": if not self.hp_polarizabilities: if not hasattr(self, 'polarizabilities'): self.polarizabilities = [] polarizability = numpy.zeros(shape=(3, 3)) indices = numpy.tril_indices(3) polarizability[indices] = [self.float(x) for x in [line[23:31], line[31:39], line[39:47], line[47:55], line[55:63], line[63:71]]] polarizability = utils.symmetrize(polarizability, use_triangle='lower') self.polarizabilities.append(polarizability) # IRC Computation convergence checks. # # -------- Sample extract for IRC step -------- # # IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC # ------------------------------------------------------------------------ # INPUT DATA FOR L123 # ------------------------------------------------------------------------ # GENERAL PARAMETERS: # Follow reaction path in both directions. # Maximum points per path = 200 # Step size = 0.100 bohr # Integration scheme = HPC # Redo corrector integration= Yes # DWI Weight Power = 2 # DWI will use Hessian update vectors when possible. # Max correction cycles = 50 # Initial Hessian = CalcFC # Hessian evaluation = Analytic every 5 predictor steps # = Analytic every 5 corrector steps # Hessian updating method = Bofill # ------------------------------------------------------------------------ # IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC # # -------- Sample extract for converged step -------- # IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC # Pt 1 Step number 1 out of a maximum of 50 # Modified Bulirsch-Stoer Extrapolation Cycles: # EPS = 0.000010000000000 # Maximum DWI energy std dev = 0.000000595 at pt 1 # Maximum DWI gradient std dev = 0.135684493 at pt 2 # CORRECTOR INTEGRATION CONVERGENCE: # Recorrection delta-x convergence threshold: 0.010000 # Delta-x Convergence Met # Point Number: 1 Path Number: 1 # CHANGE IN THE REACTION COORDINATE = 0.16730 # NET REACTION COORDINATE UP TO THIS POINT = 0.16730 # # OF POINTS ALONG THE PATH = 1 # # OF STEPS = 1 # # Calculating another point on the path. # IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC # # -------- Sample extract for unconverged intermediate step -------- # IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC # Error in corrector energy = -0.0000457166 # Magnitude of corrector gradient = 0.0140183779 # Magnitude of analytic gradient = 0.0144021969 # Magnitude of difference = 0.0078709968 # Angle between gradients (degrees)= 32.1199 # Pt 40 Step number 2 out of a maximum of 20 # Modified Bulirsch-Stoer Extrapolation Cycles: # EPS = 0.000010000000000 # Maximum DWI energy std dev = 0.000007300 at pt 31 # Maximum DWI gradient std dev = 0.085197906 at pt 59 # CORRECTOR INTEGRATION CONVERGENCE: # Recorrection delta-x convergence threshold: 0.010000 # Delta-x Convergence NOT Met # IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC-IRC if line[1:20] == "INPUT DATA FOR L123": # First IRC step if not hasattr(self, "optstatus"): self.optstatus = [] self.optstatus.append(data.ccData.OPT_NEW) if line[3:22] == "Delta-x Convergence": if line[23:30] == "NOT Met": self.optstatus[-1] += data.ccData.OPT_UNCONVERGED elif line[23:26] == "Met": self.optstatus[-1] += data.ccData.OPT_DONE if not hasattr(self, 'optdone'): self.optdone = [] self.optdone.append(len(self.optstatus) - 1) if line[:31] == ' Normal termination of Gaussian': self.metadata['success'] = True