Coverage for cclib/parser/mopacparser.py : 97%
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 MOPAC output files"""
10# Based on parser in RMG-Py by Greg Magoon
11# https://github.com/ReactionMechanismGenerator/RMG-Py/blob/master/external/cclib/parser/mopacparser.py
12# Also parts from Ben Albrecht
13# https://github.com/ben-albrecht/cclib/blob/master/cclib/parser/mopacparser.py
14# Merged and modernized by Geoff Hutchison
16import re
17import math
19import numpy
21from cclib.parser import data
22from cclib.parser import logfileparser
23from cclib.parser import utils
26def symbol2int(symbol):
27 t = utils.PeriodicTable()
28 return t.number[symbol]
30class MOPAC(logfileparser.Logfile):
31 """A MOPAC20XX output file."""
33 def __init__(self, *args, **kwargs):
35 # Call the __init__ method of the superclass
36 super(MOPAC, self).__init__(logname="MOPAC", *args, **kwargs)
38 def __str__(self):
39 """Return a string representation of the object."""
40 return "MOPAC log file %s" % (self.filename)
42 def __repr__(self):
43 """Return a representation of the object."""
44 return 'MOPAC("%s")' % (self.filename)
46 def normalisesym(self, label):
47 """MOPAC does not require normalizing symmetry labels."""
48 return label
50 def before_parsing(self):
51 #TODO
53 # Defaults
54 charge = 0
55 self.set_attribute('charge', charge)
56 mult = 1
57 self.set_attribute('mult', mult)
59 # Keep track of whether or not we're performing an
60 # (un)restricted calculation.
61 self.unrestricted = False
62 self.is_rohf = False
64 # Keep track of 1SCF vs. gopt since gopt is default
65 self.onescf = False
66 self.geomdone = False
68 # Compile the dashes-and-or-spaces-only regex.
69 self.re_dashes_and_spaces = re.compile(r'^[\s-]+$')
71 self.star = ' * '
72 self.stars = ' *******************************************************************************'
74 self.spinstate = {'SINGLET': 1,
75 'DOUBLET': 2,
76 'TRIPLET': 3,
77 'QUARTET': 4,
78 'QUINTET': 5,
79 'SEXTET': 6,
80 'HEPTET': 7,
81 'OCTET': 8,
82 'NONET': 9}
84 def extract(self, inputfile, line):
85 """Extract information from the file object inputfile."""
87 # Extract the package version.
88 if "For non-commercial use only" in line:
89 # Ignore the platorm information for now (the last character).
90 self.metadata["package_version"] = line.split()[8][:-1]
91 # Use the year as the legacy (short) package version.
92 self.skip_lines(
93 inputfile, ["Stewart Computational Chemistry", "s", "s", "s", "s"]
94 )
95 self.metadata["legacy_package_version"] = next(inputfile).split()[1][5:]
97 # Extract the atomic numbers and coordinates from the optimized geometry
98 # note that cartesian coordinates section occurs multiple times in the file, and we want to end up using the last instance
99 # also, note that the section labeled cartesian coordinates doesn't have as many decimal places as the one used here
100 # Example 1 (not used):
101 # CARTESIAN COORDINATES
102 #
103 # NO. ATOM X Y Z
104 #
105 # 1 O 4.7928 -0.8461 0.3641
106 # 2 O 5.8977 -0.3171 0.0092
107 # ...
108 # Example 2 (used):
109 # ATOM CHEMICAL X Y Z
110 # NUMBER SYMBOL (ANGSTROMS) (ANGSTROMS) (ANGSTROMS)
111 #
112 # 1 O 4.79280259 * -0.84610232 * 0.36409474 *
113 # 2 O 5.89768035 * -0.31706418 * 0.00917035 *
114 # ... etc.
115 if line.split() == ["NUMBER", "SYMBOL", "(ANGSTROMS)", "(ANGSTROMS)", "(ANGSTROMS)"]:
117 self.updateprogress(inputfile, "Attributes", self.cupdate)
119 self.inputcoords = []
120 self.inputatoms = []
122 blankline = inputfile.next()
124 atomcoords = []
125 line = inputfile.next()
126 while len(line.split()) > 6:
127 # MOPAC Version 14.019L 64BITS suddenly appends this block with
128 # "CARTESIAN COORDINATES" block with no blank line.
129 tokens = line.split()
130 self.inputatoms.append(symbol2int(tokens[1]))
131 xc = float(tokens[2])
132 yc = float(tokens[4])
133 zc = float(tokens[6])
134 atomcoords.append([xc, yc, zc])
135 line = inputfile.next()
137 self.inputcoords.append(atomcoords)
139 if not hasattr(self, "natom"):
140 self.atomnos = numpy.array(self.inputatoms, 'i')
141 self.natom = len(self.atomnos)
143 if 'CHARGE ON SYSTEM =' in line:
144 charge = int(line.split()[5])
145 self.set_attribute('charge', charge)
147 if 'SPIN STATE DEFINED' in line:
148 # find the multiplicity from the line token (SINGLET, DOUBLET, TRIPLET, etc)
149 mult = self.spinstate[line.split()[1]]
150 self.set_attribute('mult', mult)
152 # Read energy (in kcal/mol, converted to eV)
153 #
154 # FINAL HEAT OF FORMATION = -333.88606 KCAL = -1396.97927 KJ
155 if 'FINAL HEAT OF FORMATION =' in line:
156 if not hasattr(self, "scfenergies"):
157 self.scfenergies = []
158 self.scfenergies.append(utils.convertor(utils.float(line.split()[5]), "kcal/mol", "eV"))
160 # Molecular mass parsing (units will be amu)
161 #
162 # MOLECULAR WEIGHT == 130.1890
163 if line[0:35] == ' MOLECULAR WEIGHT =':
164 self.molmass = utils.float(line.split()[3])
166 #rotational constants
167 #Example:
168 # ROTATIONAL CONSTANTS IN CM(-1)
169 #
170 # A = 0.01757641 B = 0.00739763 C = 0.00712013
171 # could also read in moment of inertia, but this should just differ by a constant: rot cons= h/(8*Pi^2*I)
172 # note that the last occurence of this in the thermochemistry section has reduced precision,
173 # so we will want to use the 2nd to last instance
174 if line[0:40] == ' ROTATIONAL CONSTANTS IN CM(-1)':
175 blankline = inputfile.next()
176 rotinfo = inputfile.next()
177 if not hasattr(self, "rotcons"):
178 self.rotcons = []
179 broken = rotinfo.split()
180 # leave the rotational constants in Hz
181 a = float(broken[2])
182 b = float(broken[5])
183 c = float(broken[8])
184 self.rotcons.append([a, b, c])
186 # Start of the IR/Raman frequency section.
187 # Example:
188 # VIBRATION 1 1A ATOM PAIR ENERGY CONTRIBUTION RADIAL
189 # FREQ. 15.08 C 12 -- C 16 +7.9% (999.0%) 0.0%
190 # T-DIPOLE 0.2028 C 16 -- H 34 +5.8% (999.0%) 28.0%
191 # TRAVEL 0.0240 C 16 -- H 32 +5.6% (999.0%) 35.0%
192 # RED. MASS 1.7712 O 1 -- O 4 +5.2% (999.0%) 0.4%
193 # EFF. MASS7752.8338
194 #
195 # VIBRATION 2 2A ATOM PAIR ENERGY CONTRIBUTION RADIAL
196 # FREQ. 42.22 C 11 -- C 15 +9.0% (985.8%) 0.0%
197 # T-DIPOLE 0.1675 C 15 -- H 31 +6.6% (843.6%) 3.3%
198 # TRAVEL 0.0359 C 15 -- H 29 +6.0% (802.8%) 24.5%
199 # RED. MASS 1.7417 C 13 -- C 17 +5.8% (792.7%) 0.0%
200 # EFF. MASS1242.2114
201 if line[1:10] == 'VIBRATION':
202 self.updateprogress(inputfile, "Frequency Information", self.fupdate)
204 # get the vib symmetry
205 if len(line.split()) >= 3:
206 sym = line.split()[2]
207 if not hasattr(self, 'vibsyms'):
208 self.vibsyms = []
209 self.vibsyms.append(sym)
211 line = inputfile.next()
212 if 'FREQ' in line:
213 if not hasattr(self, 'vibfreqs'):
214 self.vibfreqs = []
215 freq = float(line.split()[1])
216 self.vibfreqs.append(freq)
218 line = inputfile.next()
219 if 'T-DIPOLE' in line:
220 if not hasattr(self, 'vibirs'):
221 self.vibirs = []
222 tdipole = float(line.split()[1])
223 # transform to km/mol
224 self.vibirs.append(math.sqrt(tdipole))
226 line = inputfile.next()
227 if 'TRAVEL' in line:
228 pass
230 line = inputfile.next()
231 if 'RED. MASS' in line:
232 if not hasattr(self, 'vibrmasses'):
233 self.vibrmasses = []
234 rmass = float(line.split()[2])
235 self.vibrmasses.append(rmass)
237 # Orbital eigenvalues, e.g.
238 # ALPHA EIGENVALUES
239 # BETA EIGENVALUES
240 # or just "EIGENVALUES" for closed-shell
241 if 'EIGENVALUES' in line:
242 if not hasattr(self, 'moenergies'):
243 self.moenergies = [] # list of arrays
245 energies = []
246 line = inputfile.next()
247 while len(line.split()) > 0:
248 energies.extend([float(i) for i in line.split()])
249 line = inputfile.next()
250 self.moenergies.append(energies)
252 # todo:
253 # Partial charges and dipole moments
254 # Example:
255 # NET ATOMIC CHARGES
257 if line[:16] == '== MOPAC DONE ==':
258 self.metadata['success'] = True