|
From: Ramon C. <rc...@iq...> - 2014-03-28 08:11:24
|
Dear all,
I'm trying to parse Gaussian 03 IRC calculations and I face two problems:
1) My IRC had some instability at the end of the calculation, therefore some big
numbers could not be properly formatted, and parsing the output gives an error
(see attached file). Fair enough, but it would be good if I could still parse
the file and get all the previous points. Would that be possible? Here is a line
that gives an error:
Z-Matrix orientation:
---------------------------------------------------------------------
Center Atomic Atomic Coordinates (Angstroms)
Number Number Type X Y Z
---------------------------------------------------------------------
1 6 0 42354.982772************93749.516872
2 1 0 42355.937399************93750.035409
3 1 0 42354.724535************93748.877845
4 8 0 42354.612001************93749.055343
5 8 0 42353.940461************93750.217216
---------------------------------------------------------------------
2) If I parse an IRC I get the energy of all the points, either optimized or
not. It would be useful to me to be able to get which points are the final
points of the optimization (for example if data.optdone was a boolean array for
all the points). For that, I have modified gaussianparser.py and data.py to work
this way. I enclose the source files, as they can be useful to others or to the
development team. However I have only tested them with my log file. They need
further testing.
Cheers,
Ramon
--
Ramon Crehuet
Cientific Titular (Assistant Professor)
Institute of Advanced Chemistry of Catalonia
IQAC - CSIC
http://www.iqac.csic.es/qteor
https://twitter.com/rcrehuet
http://ramoncrehuet.wordpress.com/
Tel. +34 934006116
Jordi Girona 18-26
08034 Barcelona (Spain)
|
|
From: Ramon C. <rc...@iq...> - 2014-03-28 10:35:36
Attachments:
gaussianparser.py
|
Hi again, I'll answer my own question, because I've already modified gaussianparser.py so that it deals with missing coordinates. I took the idea from the lines reading the Forces. I just changed a few lines from line 190. I attach the file in case it can be useful to the community. I know I should create a pull request on github, but I'm still a bit of a newby in github. Next time... Ramon PS. This version of gaussianparser.py also includes my previous contribution where optdone attribute is a list of booleans, not a boolean. On 28/03/14 10:20, ccl...@li... wrote: > Dear all, > I'm trying to parse Gaussian 03 IRC calculations and I face two problems: > 1) My IRC had some instability at the end of the calculation, therefore some big > numbers could not be properly formatted, and parsing the output gives an error > (see attached file). Fair enough, but it would be good if I could still parse > the file and get all the previous points. Would that be possible? Here is a line > that gives an error: > > Z-Matrix orientation: > --------------------------------------------------------------------- > Center Atomic Atomic Coordinates (Angstroms) > Number Number Type X Y Z > --------------------------------------------------------------------- > 1 6 0 42354.982772************93749.516872 > 2 1 0 42355.937399************93750.035409 > 3 1 0 42354.724535************93748.877845 > 4 8 0 42354.612001************93749.055343 > 5 8 0 42353.940461************93750.217216 > --------------------------------------------------------------------- > > 2) If I parse an IRC I get the energy of all the points, either optimized or > not. It would be useful to me to be able to get which points are the final > points of the optimization (for example if data.optdone was a boolean array for > all the points). For that, I have modified gaussianparser.py and data.py to work > this way. I enclose the source files, as they can be useful to others or to the > development team. However I have only tested them with my log file. They need > further testing. > > Cheers, > Ramon -- Ramon Crehuet Cientific Titular (Assistant Professor) Institute of Advanced Chemistry of Catalonia IQAC - CSIC http://www.iqac.csic.es/qteor https://twitter.com/rcrehuet http://ramoncrehuet.wordpress.com/ Tel. +34 934006116 Jordi Girona 18-26 08034 Barcelona (Spain) |
|
From: Adam T. <ate...@gm...> - 2014-03-28 14:41:51
|
Hi Ramon, Thank you for alerting us to these issues, and providing an example logfile and patches. Just to clarify: you are willing to place the logfile in the public domain? Also, for cclib v1.3 (6-12 months?), I do plan on changing optdone from a bool to a list of integers corresponding to the energy/atomcoords index of when an optimization has finished. This would make it easy to select the energies and coordinates that correspond to points on a IRC or relaxed PES. Thanks again, Adam On Fri, Mar 28, 2014 at 3:35 AM, Ramon Crehuet <rc...@iq...> wrote: > Hi again, > I'll answer my own question, because I've already modified > gaussianparser.py so that it deals with missing coordinates. I took the > idea from the lines reading the Forces. I just changed a few lines from > line 190. I attach the file in case it can be useful to the community. > I know I should create a pull request on github, but I'm still a bit of a > newby in github. Next time... > Ramon > > PS. This version of gaussianparser.py also includes my previous > contribution where optdone attribute is a list of booleans, not a boolean. > > > > On 28/03/14 10:20, ccl...@li... wrote: > >> Dear all, >> I'm trying to parse Gaussian 03 IRC calculations and I face two problems: >> 1) My IRC had some instability at the end of the calculation, therefore >> some big >> numbers could not be properly formatted, and parsing the output gives an >> error >> (see attached file). Fair enough, but it would be good if I could still >> parse >> the file and get all the previous points. Would that be possible? Here is >> a line >> that gives an error: >> >> Z-Matrix orientation: >> --------------------------------------------------------------------- >> Center Atomic Atomic Coordinates (Angstroms) >> Number Number Type X Y Z >> --------------------------------------------------------------------- >> 1 6 0 42354.982772************93749.516872 >> 2 1 0 42355.937399************93750.035409 >> 3 1 0 42354.724535************93748.877845 >> 4 8 0 42354.612001************93749.055343 >> 5 8 0 42353.940461************93750.217216 >> --------------------------------------------------------------------- >> >> 2) If I parse an IRC I get the energy of all the points, either optimized >> or >> not. It would be useful to me to be able to get which points are the final >> points of the optimization (for example if data.optdone was a boolean >> array for >> all the points). For that, I have modified gaussianparser.py and data.py >> to work >> this way. I enclose the source files, as they can be useful to others or >> to the >> development team. However I have only tested them with my log file. They >> need >> further testing. >> >> Cheers, >> Ramon >> > > -- > Ramon Crehuet > Cientific Titular (Assistant Professor) > Institute of Advanced Chemistry of Catalonia > IQAC - CSIC > http://www.iqac.csic.es/qteor > https://twitter.com/rcrehuet > http://ramoncrehuet.wordpress.com/ > Tel. +34 934006116 > Jordi Girona 18-26 > 08034 Barcelona (Spain) > > > > ------------------------------------------------------------------------------ > > _______________________________________________ > cclib-users mailing list > ccl...@li... > https://lists.sourceforge.net/lists/listinfo/cclib-users > > |
|
From: Karol M. L. <kar...@gm...> - 2014-03-28 13:16:20
|
Hola Ramon, Thanks for the contribution! Could you confirm that you place your logfiles in the public domain, so that we can use them for testing cclib? I will incorporate your changes into the Gaussian parser. Of course, we welcome suggestions and pull requests also on github in the future. We have some instructions on how to contribute in the docs (http://cclib.github.io/development.html) and will add more soon. If you have any questions, feel free to ask. Concerning optdone, you've actually anticipated a change that has already happened, because in the development version optdone is in fact a list. It is a list of integers, however, not Booleans. Each element is an index for atomcorods and scgenergies where the optimization has converged. This will be the default for all versions after v1.3. Best, Karol On Mar 28 2014, Ramon Crehuet wrote: > Hi again, > I'll answer my own question, because I've already modified > gaussianparser.py so that it deals with missing coordinates. I took > the idea from the lines reading the Forces. I just changed a few > lines from line 190. I attach the file in case it can be useful to > the community. > I know I should create a pull request on github, but I'm still a bit > of a newby in github. Next time... > Ramon > > PS. This version of gaussianparser.py also includes my previous > contribution where optdone attribute is a list of booleans, not a > boolean. > > > On 28/03/14 10:20, ccl...@li... wrote: > >Dear all, > >I'm trying to parse Gaussian 03 IRC calculations and I face two problems: > >1) My IRC had some instability at the end of the calculation, therefore some big > >numbers could not be properly formatted, and parsing the output gives an error > >(see attached file). Fair enough, but it would be good if I could still parse > >the file and get all the previous points. Would that be possible? Here is a line > >that gives an error: > > > > Z-Matrix orientation: > > --------------------------------------------------------------------- > > Center Atomic Atomic Coordinates (Angstroms) > > Number Number Type X Y Z > > --------------------------------------------------------------------- > > 1 6 0 42354.982772************93749.516872 > > 2 1 0 42355.937399************93750.035409 > > 3 1 0 42354.724535************93748.877845 > > 4 8 0 42354.612001************93749.055343 > > 5 8 0 42353.940461************93750.217216 > > --------------------------------------------------------------------- > > > >2) If I parse an IRC I get the energy of all the points, either optimized or > >not. It would be useful to me to be able to get which points are the final > >points of the optimization (for example if data.optdone was a boolean array for > >all the points). For that, I have modified gaussianparser.py and data.py to work > >this way. I enclose the source files, as they can be useful to others or to the > >development team. However I have only tested them with my log file. They need > >further testing. > > > >Cheers, > >Ramon > > -- > Ramon Crehuet > Cientific Titular (Assistant Professor) > Institute of Advanced Chemistry of Catalonia > IQAC - CSIC > http://www.iqac.csic.es/qteor > https://twitter.com/rcrehuet > http://ramoncrehuet.wordpress.com/ > Tel. +34 934006116 > Jordi Girona 18-26 > 08034 Barcelona (Spain) > > # This file is part of cclib (http://cclib.sf.net), a library for parsing > # and interpreting the results of computational chemistry packages. > # > # Copyright (C) 2006, the cclib development team > # > # The library is free software, distributed under the terms of > # the GNU Lesser General Public version 2.1 or later. You should have > # received a copy of the license along with cclib. You can also access > # the full license online at http://www.gnu.org/copyleft/lgpl.html. > > import re > > import numpy > > from . import logfileparser > from . 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 > > >>> sym = Gaussian("dummyfile").normalisesym > >>> labels = ['A1', 'AG', 'A1G', "SG", "PI", "PHI", "DLTA", 'DLTU', 'SGG'] > >>> map(sym, labels) > ['A1', 'Ag', 'A1g', 'sigma', 'pi', 'phi', 'delta', 'delta.u', 'sigma.g'] > """ > # 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)) > > # Flag that indicates whether it has reached the end of a geoopt. > self.optfinished = False > > # 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 > > 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, 'enthaply') and hasattr(self, 'temperature') > and hasattr(self, 'freeenergy')): > self.entropy = (self.enthaply - self.freeenergy)/self.temperature > > def extract(self, inputfile, line): > """Extract information from the file object inputfile.""" > > #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] > > #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: > if not hasattr(self, 'enthaply'): > self.enthaply = float(line.split()[6]) > if "Sum of electronic and thermal Free Energies=" in line: > if not hasattr(self, 'freeenergy'): > self.freeenergy = float(line.split()[7]) > if line[1:12] == "Temperature": > if not hasattr(self, 'temperature'): > self.temperature = float(line.split()[1]) > > # Number of atoms. > if line[1:8] == "NAtoms=": > > self.updateprogress(inputfile, "Attributes", self.fupdate) > > natom = int(line.split()[1]) > if not hasattr(self, "natom"): > self.natom = natom > > # Catch message about completed optimization. > if line[1:23] == "Optimization completed": > self.optfinished = True > self.optdone[-1] = True > > # Extract the atomic numbers and coordinates from the input orientation, > # in the event the standard orientation isn't available. > if not self.optfinished and 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 hasattr(self, "inputcoords"): > self.inputcoords = [] > self.inputatoms = [] > > if not hasattr(self, "optdone"): > self.optdone = [] > > hyphens = next(inputfile) > colmNames = next(inputfile) > colmNames = next(inputfile) > hyphens = next(inputfile) > > atomcoords = [] > line = next(inputfile) > while line != hyphens: > tmpcoords= [] > for N in range(3): #X, Y, Z > coord = line[34+N*12:46+N*12] > if coord.startswith("*"): > coord = "NaN" > tmpcoords.append(float(coord)) > broken = line.split() > self.inputatoms.append(int(broken[1])) > atomcoords.append(tmpcoords) > line = next(inputfile) > > self.inputcoords.append(atomcoords) > self.optdone.append(False) > > if not hasattr(self, "atomnos"): > self.atomnos = numpy.array(self.inputatoms, 'i') > self.natom = len(self.atomnos) > > # 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 not self.optfinished and 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 = [] > > if not hasattr(self, "optdone"): > self.optdone = [] > > hyphens = next(inputfile) > colmNames = next(inputfile) > colmNames = next(inputfile) > hyphens = next(inputfile) > > atomnos = [] > atomcoords = [] > line = next(inputfile) > while line != hyphens: > broken = line.split() > atomnos.append(int(broken[1])) > atomcoords.append(list(map(float, broken[-3:]))) > line = next(inputfile) > self.atomcoords.append(atomcoords) > self.optdone.append(False) > > if not hasattr(self, "natom"): > self.atomnos = numpy.array(atomnos, 'i') > self.natom = len(self.atomnos) > > # make sure atomnos is added for the case where natom has already been set > elif not hasattr(self, "atomnos"): > self.atomnos = numpy.array(atomnos, 'i') > > # Find the targets for SCF convergence (QM calcs). > if 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: > > parts = line.split() > newlist = [self.float(x.split('=')[1]) for x in parts[0:2]] > energy = 1.0 > if len(parts) > 4: > energy = parts[2].split('=')[1] > if energy == "": > energy = self.float(parts[3]) > else: > energy = self.float(energy) > if len(self.scftargets[0]) == 3: # Only add the energy if it's a target criteria > newlist.append(energy) > scfvalues.append(newlist) > > try: > line = next(inputfile) > # May be interupted by EOF. > except StopIteration: > break > > self.scfvalues.append(scfvalues) > > # Extract SCF convergence information (AM1 calcs). > if line[1:4] == 'It=': > > self.scftargets = numpy.array([1E-7], "d") # This is the target value for the rms > self.scfvalues = [[]] > > line = next(inputfile) > while line.find(" Energy") == -1: > > if self.progress: > step = inputfile.tell() > if step != oldstep: > self.progress.update(step, "AM1 Convergence") > oldstep = step > > if line[1:4] == "It=": > parts = line.strip().split() > self.scfvalues[0].append(self.float(parts[-1][:-1])) > line = next(inputfile) > > # Note: this needs to follow the section where 'SCF Done' is used > # to terminate a loop when extracting SCF convergence information. > if line[1:9] == 'SCF Done': > > 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]: > > 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": > > 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)": > > 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": > 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.ccenergy = self.float(line.split()[3]) > if line[1:10] == "T5(CCSD)=": > line = next(inputfile) > if line[1:9] == "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 > > # 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? > if line[37:43] == "Forces": > > if not hasattr(self, "grads"): > self.grads = [] > > header = next(inputfile) > dashes = next(inputfile) > line = next(inputfile) > forces = [] > while line != dashes: > 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) > > # Charge and multiplicity. > # If counterpoise correction is used, multiple lines match. > # The first one contains charge/multiplicity of the whole molecule.: > # 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: > > regex = ".*=(.*)Mul.*=\s*-?(\d+).*" > match = re.match(regex, line) > assert match, "Something unusual about the line: '%s'" % line > > if not hasattr(self, "charge"): > self.charge = int(match.groups()[0]) > if not hasattr(self, "mult"): > self.mult = int(match.groups()[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 = numpy.array([i-1], "i") # 'HOMO' indexes the HOMO in the arrays > 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.resize([2]) > self.homos[1] = i-1 > else: > # 'HOMO' indexes the HOMO in the arrays > self.homos = numpy.array([i-1], "i") > parts = line[17:].split() > for x in parts: > self.mosyms[1].append(self.normalisesym(x.strip('()'))) > i += 1 > line = next(inputfile) > > # 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.fupdat... [truncated message content] |