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. 7import collections 8 9"""Parser for Turbomole output files.""" 10 11import re 12 13import numpy 14 15from cclib.parser import logfileparser 16from cclib.parser import utils 17from cclib.parser import data 18 19class AtomBasis: 20 def __init__(self, atname, basis_name, inputfile): 21 self.symmetries=[] 22 self.coefficients=[] 23 self.atname=atname 24 self.basis_name=basis_name 25 26 self.parse_basis(inputfile) 27 28 def parse_basis(self, inputfile): 29 i=0 30 line=inputfile.next() 31 32 while(line[0]!="*"): 33 (nbasis_text, symm)=line.split() 34 self.symmetries.append(symm) 35 36 nbasis=int(nbasis_text) 37 coeff_arr=numpy.zeros((nbasis, 2), float) 38 39 for j in range(0, nbasis, 1): 40 line=inputfile.next() 41 (e1_text, e2_text)=line.split() 42 coeff_arr[j][0]=float(e1_text) 43 coeff_arr[j][1]=float(e2_text) 44 45 self.coefficients.append(coeff_arr) 46 line=inputfile.next() 47 48class Turbomole(logfileparser.Logfile): 49 """A Turbomole log file.""" 50 51 def __init__(self, *args, **kwargs): 52 super(Turbomole, self).__init__(logname="Turbomole", *args, **kwargs) 53 54 # Flag for whether this calc is DFT. 55 self.is_DFT = False 56 57 # A Regex that we use to extract version info. 58 self.version_regex = re.compile(r"TURBOMOLE(?: rev\.)? V([\d]+)[.-]([\d]+)(?:[.-]([\d]))?(?: \( ?([0-9A-z]+) ?\))?") 59 60 # A list of previous lines to allow look-behind functionality. 61 self.last_lines = collections.deque([""] * 10, 10) 62 63 def __str__(self): 64 """Return a string representation of the object.""" 65 return "Turbomole output file %s" % (self.filename) 66 67 def __repr__(self): 68 """Return a representation of the object.""" 69 return 'Turbomole("%s")' % (self.filename) 70 71 def normalisesym(self, label): 72 """Normalise the symmetries used by Turbomole. 73 74 The labels are standardized except for the first character being lowercase. 75 """ 76 # TODO more work could be required, but we don't have any logfiles 77 # with non-C1 symmetry. 78 return label.capitalize() 79 80 def before_parsing(self): 81 82 self.periodic_table = utils.PeriodicTable() 83 84 @staticmethod 85 def split_molines(inline): 86 """Splits the lines containing mocoeffs (each of length 20) 87 and converts them to float correctly. 88 """ 89 line = inline.replace("D", "E") 90 f1 = line[0:20] 91 f2 = line[20:40] 92 f3 = line[40:60] 93 f4 = line[60:80] 94 95 if(len(f4) > 1): 96 return [float(f1), float(f2), float(f3), float(f4)] 97 if(len(f3) > 1): 98 return [float(f1), float(f2), float(f3)] 99 if(len(f2) > 1): 100 return [float(f1), float(f2)] 101 if(len(f1) > 1): 102 return [float(f1)] 103 104 def extract(self, inputfile, line): 105 """Extract information from the file object inputfile.""" 106 107 ## This information is in the control file. 108 # $rundimensions 109 # dim(fock,dens)=1860 110 # natoms=20 111 # nshell=40 112 # nbf(CAO)=60 113 # nbf(AO)=60 114 # dim(trafo[SAO<-->AO/CAO])=60 115 # rhfshells=1 116 if line[3:10]=="natoms=": 117 self.natom=int(line[10:]) 118 119 if line[3:11] == "nbf(AO)=": 120 nmo = int(line.split('=')[1]) 121 self.set_attribute('nbasis', nmo) 122 123 # The DFT functional. 124 # This information is printed by dscf but not in an easily parsable format, so we'll take it from the control file instead... 125 # Additionally, turbomole stores functional names in lower case. This looks odd, so we'll convert to uppercase (?) 126 # We are parsing this section from the control file. 127 if line[3:13] == "functional": 128 self.metadata['functional'] = line.split()[1].upper() 129 self.is_DFT = True 130 131 # Information about DFT is also printed by dscf in the main .log file. 132 # We don't parse this at the moment, but it is still important to know 133 # if we are using DFT for parsing later on (to distinguish between DFT 134 # Vs HF etc). If we don't have the control file available, we will 135 # need this check: 136 if "density functional" in line: 137 self.is_DFT = True 138 139 # Extract the version number and optionally the build number. 140 version_match = self.version_regex.search(line) 141 if version_match: 142 # We only combine the parts of the version string we actually got. 143 # Our regex ensures we have at least the first two parts (x.y), and optionally a third part (x.y.z). 144 # This line could look like any of the following: 145 # - TURBOMOLE rev. V7.4.1 (2bfdd732) 146 # - TURBOMOLE V7.2 ( 21471 ) 11 Oct 2017 at 17:04:51 147 # - TURBOMOLE V5-9-0 29 Nov 2006 at 22:06:41 148 version = ".".join([version_part for version_part in version_match.group(1,2,3) if version_part is not None]) 149 150 # We may also have a build ID. 151 build_id = version_match.group(4) 152 153 self.metadata["legacy_package_version"] = version 154 self.metadata["package_version"] = "{}.r{}".format(version, build_id) if build_id is not None else version 155 156 # We have entered a new module (sub program); reset our success flag. 157 self.metadata['success'] = False 158 159 ## Orbital occupation info from dscf. 160 # orbitals $scfmo will be written to file mos 161 # 162 # irrep 1a 2a 3a 4a 5a 163 # eigenvalues H -20.25992 -1.24314 -0.57053 -0.46144 -0.39295 164 # eV -551.3047 -33.8279 -15.5250 -12.5564 -10.6929 165 # occupation 2.0000 2.0000 2.0000 2.0000 2.0000 166 # 167 # irrep 6a 7a 168 # eigenvalues H 0.55091 0.64409 169 # eV 14.9910 17.5268 170 ## Or 171 # orbitals $uhfmo_beta will be written to file beta 172 # 173 # orbitals $uhfmo_alpha will be written to file alpha 174 # 175 # alpha: 176 # 177 # irrep 31a 32a 33a 34a 35a 178 # eigenvalues H -0.47570 -0.46573 -0.40741 -0.39213 -0.35411 179 # eV -12.9446 -12.6733 -11.0862 -10.6705 -9.6358 180 # occupation 1.0000 1.0000 1.0000 1.0000 1.0000 181 # 182 # irrep 36a 37a 38a 39a 40a 183 # eigenvalues H -0.18634 -0.10035 -0.09666 -0.02740 0.06072 184 # eV -5.0705 -2.7306 -2.6303 -0.7455 1.6522 185 # 186 # beta: 187 # 188 # irrep 30a 31a 32a 33a 34a 189 # eigenvalues H -0.49118 -0.47348 -0.44470 -0.39020 -0.37919 190 # eV -13.3658 -12.8842 -12.1009 -10.6181 -10.3184 191 # occupation 1.0000 1.0000 1.0000 1.0000 1.0000 192 # 193 # irrep 35a 36a 37a 38a 39a 194 # eigenvalues H -0.28091 -0.15088 -0.09343 -0.07531 -0.00688 195 # eV -7.6440 -4.1058 -2.5424 -2.0493 -0.1873 196 # There's no point parsing HOMO/LUMO if we don't already have orbitals. 197 if hasattr(self, "mosyms"): 198 if "will be written to file mos" in line: 199 orbitals, line = self.parse_dscf_orbitals(inputfile, line) 200 self.set_attribute("homos", [self.determine_homo(self.mosyms[0], orbitals)]) 201 202 if "alpha:" in line: 203 orbitals, line = self.parse_dscf_orbitals(inputfile, line) 204 if len(orbitals) > 0: 205 homo = self.determine_homo(self.mosyms[0], orbitals) 206 if not hasattr(self, "homos"): 207 self.set_attribute('homos', [homo]) 208 else: 209 self.homos[0] = homo 210 211 if "beta:" in line: 212 orbitals, line = self.parse_dscf_orbitals(inputfile, line) 213 if len(orbitals) > 0: 214 homo = self.determine_homo(self.mosyms[1], orbitals) 215 if not hasattr(self, "homos"): 216 self.set_attribute('homos', [homo]) 217 elif len(self.homos) == 1: 218 self.homos.append(homo) 219 else: 220 self.homos[1] = homo 221 222 # Molecular charge and dipole info from dscf. 223 # ============================================================================== 224 # electrostatic moments 225 # ============================================================================== 226 # 227 # reference point for electrostatic moments: 0.00000 0.00000 0.00000 228 # 229 # 230 # nuc elec -> total 231 # ------------------------------------------------------------------------------ 232 # charge 233 # ------------------------------------------------------------------------------ 234 # 70.000000 -70.000000 -0.000000 235 # 236 # ------------------------------------------------------------------------------ 237 # dipole moment 238 # ------------------------------------------------------------------------------ 239 # x -0.000000 0.000000 -0.000000 240 # y 0.001384 -0.001340 0.000044 241 # z -0.000000 0.000000 -0.000000 242 # 243 # | dipole moment | = 0.0000 a.u. = 0.0001 debye 244 # 245 # ------------------------------------------------------------------------------ 246 # quadrupole moment 247 # ------------------------------------------------------------------------------ 248 # xx 1499.472650 -1537.186938 -37.714287 249 # yy 0.000002 -43.588053 -43.588051 250 # zz 244.507989 -281.855472 -37.347483 251 # xy -0.000000 0.000000 -0.000000 252 # xz -5.011477 5.064680 0.053203 253 # yz -0.000000 0.000000 0.000000 254 # 255 # 1/3 trace= -39.549940 256 # anisotropy= 6.066190 257 if "reference point for electrostatic moments:" in line: 258 # This indicates the start of a new dipole section. 259 # Safe to overwrite any old dipoles. 260 parts = line.split() 261 self.moments = [[ 262 utils.convertor(float(parts[-3]), "bohr", "Angstrom"), 263 utils.convertor(float(parts[-2]), "bohr", "Angstrom"), 264 utils.convertor(float(parts[-1]), "bohr", "Angstrom") 265 ]] 266 267 if "nuc elec -> total" in line: 268 line = next(inputfile) 269 line = next(inputfile) 270 if "charge" in line: 271 line = next(inputfile) 272 line = next(inputfile) 273 274 total_charge = float(line.split()[2]) 275 total_charge_int = round(total_charge) 276 277 # Check we won't loose information converting to int. 278 if total_charge != total_charge_int: 279 self.logger.warning("Converting non integer total charge '{}' to integer".format(total_charge)) 280 281 # Set regardless. 282 self.set_attribute("charge", total_charge_int) 283 284 if line.strip() == "dipole moment": 285 line = next(inputfile) 286 if set(line.strip()) == {"-"}: 287 line = next(inputfile) 288 x_coord = utils.convertor(float(line.split()[-1]), "ebohr", "Debye") 289 line = next(inputfile) 290 y_coord = utils.convertor(float(line.split()[-1]), "ebohr", "Debye") 291 line = next(inputfile) 292 z_coord = utils.convertor(float(line.split()[-1]), "ebohr", "Debye") 293 294 # Assume 0,0,0 as origin if not given. 295 if not hasattr(self, "moments"): 296 self.moments = [[0,0,0]] 297 self.moments.append([x_coord, y_coord, z_coord]) 298 299 if line.strip() == "quadrupole moment": 300 line = next(inputfile) 301 if set(line.strip()) == {"-"}: 302 line = next(inputfile) 303 xx_coord = utils.convertor(float(line.split()[-1]), "ebohr2", "Buckingham") 304 line = next(inputfile) 305 yy_coord = utils.convertor(float(line.split()[-1]), "ebohr2", "Buckingham") 306 line = next(inputfile) 307 zz_coord = utils.convertor(float(line.split()[-1]), "ebohr2", "Buckingham") 308 line = next(inputfile) 309 xy_coord = utils.convertor(float(line.split()[-1]), "ebohr2", "Buckingham") 310 line = next(inputfile) 311 xz_coord = utils.convertor(float(line.split()[-1]), "ebohr2", "Buckingham") 312 line = next(inputfile) 313 yz_coord = utils.convertor(float(line.split()[-1]), "ebohr2", "Buckingham") 314 self.moments.append([xx_coord, xy_coord, xz_coord, yy_coord, yz_coord, zz_coord]) 315 316 ## Basis set info from dscf. 317 # +--------------------------------------------------+ 318 # | basis set information | 319 # +--------------------------------------------------+ 320 # 321 # we will work with the 1s 3p 5d 7f 9g ... basis set 322 # ...i.e. with spherical basis functions... 323 # 324 # type atoms prim cont basis 325 # --------------------------------------------------------------------------- 326 # o 1 15 5 sto-3g hondo [2s1p|6s3p] 327 # h 2 3 1 sto-3g hondo [1s|3s] 328 # --------------------------------------------------------------------------- 329 if "type atoms prim cont basis" in line: 330 line = next(inputfile) 331 line = next(inputfile) 332 basis_sets = [] 333 while set(line.strip()) != {"-"}: 334 basis_sets.append(" ".join(line.split()[4:-1])) 335 line = next(inputfile) 336 337 # Turbomole gives us the basis set for each atom, but we're only interested if the same basis set is used throughout (for now). 338 if len(set(basis_sets)) == 1: 339 self.metadata["basis_set"] = list(set(basis_sets))[0] 340 341 # Coordinates and gradients from statpt. 342 # ************************************************************************* 343 # ATOM CARTESIAN COORDINATES 344 # 1 c -2.67642286424381 0.00038527796998 -0.44566112589039 345 # 2 c -1.69183504162310 0.00075591605173 2.05352357416443 346 # 3 c 0.92311977253458 0.00052122921704 2.48381308900370 347 # 4 c 2.67642286424491 0.00038528140768 0.44566112589851 348 # 5 c 1.69183504161656 0.00075592009652 -2.05352357415432 349 # 6 h -2.98983238801665 0.00149374762706 3.67174843088014 350 # 7 h 1.64279730150941 0.00038458181684 4.43158185858240 351 # 8 h 2.98983238800216 0.00149376082387 -3.67174843087059 352 # 9 c 5.44975417206469 -0.00039526372012 1.01184691031725 353 # 10 c 7.34299179214000 -0.00071986894317 -0.68271530533475 354 # 11 h 7.01332300460381 -0.00029648805084 -2.72785997302809 355 # 12 h 5.91422629512709 -0.00052260839470 3.03917498747826 356 # 13 h 9.32456490822245 -0.00139866158379 -0.07769877084777 357 # 14 c -5.44975417205833 -0.00039526324514 -1.01184691032156 358 # 15 h -5.91422629512079 -0.00052260600240 -3.03917498748360 359 # 16 c -7.34299179213657 -0.00071986698481 0.68271530532303 360 # 17 h -7.01332300460557 -0.00029648750734 2.72785997301777 361 # 18 h -9.32456490821458 -0.00139865699894 0.07769877082649 362 # 19 c -0.92311977253569 0.00052122610063 -2.48381308899233 363 # 20 h -1.64279730151053 0.00038457006746 -4.43158185856859 364 # ************************************************************************* 365 # ATOM CARTESIAN GRADIENTS 366 # 1 c 0.00001122250995 0.00000667056223 0.00001828623085 367 # 2 c 0.00000193381649 -0.00001184358472 0.00001337145426 368 # 3 c -0.00000578603102 0.00000187474419 0.00000987215044 369 # 4 c -0.00001122250794 0.00000667041002 -0.00001828623169 370 # 5 c -0.00000193381900 -0.00001184350870 -0.00001337145668 371 # 6 h -0.00001402399615 0.00001821929428 -0.00000611193875 372 # 7 h 0.00001742260704 -0.00000462667818 -0.00000195847724 373 # 8 h 0.00001402399610 0.00001821957100 0.00000611193919 374 # 9 c -0.00006557131494 -0.00002136001589 0.00002475454295 375 # 10 c 0.00006109842779 -0.00000007565809 -0.00004223963891 376 # 11 h -0.00002879535900 0.00000574108155 -0.00000337279221 377 # 12 h 0.00004754086240 0.00000863480152 -0.00000658955451 378 # 13 h -0.00000914653785 -0.00000320753859 0.00003066647865 379 # 14 c 0.00006557131503 -0.00002135978825 -0.00002475454202 380 # 15 h -0.00004754086227 0.00000863481272 0.00000658955373 381 # 16 c -0.00006109842855 -0.00000007566316 0.00004223963942 382 # 17 h 0.00002879535887 0.00000574109116 0.00000337279251 383 # 18 h 0.00000914653868 -0.00000320756462 -0.00003066647871 384 # 19 c 0.00000578603127 0.00000187466275 -0.00000987214898 385 # 20 h -0.00001742260682 -0.00000462696679 0.00000195847629 386 # ************************************************************************* 387 # 388 # In older versions (5.9) 389 # ATOM CARTESIAN GRADIENTS 390 # 1 c -0.00006249573129 -0.00001452041025 -0.00002380577094 391 # 2 c 0.00007390978833 -0.00000084052225 0.00000584172986 392 # 3 c -0.00000362984325 0.00000466314052 -0.00002733252443 393 # 4 c 0.00006249783071 -0.00001452158676 0.00002380177091 394 # 5 c -0.00007390342509 -0.00000084270316 -0.00000584066725 395 # 6 h -0.00005950589946 0.00000095087455 0.00005137853319 396 # 7 h -0.00000645073463 -0.00000179899463 0.00003328307243 397 # 8 h 0.00005950123260 0.00000095090879 -0.00005137678599 398 # 9 c -0.00004885046058 -0.00002593955699 -0.00004887642413 399 # 10 c 0.00000015222157 0.00004027444288 0.00003204866033 400 # 11 h 0.00001239310045 0.00001308636889 -0.00000092060218 401 # 12 h 0.00001561988506 -0.00001608132438 0.00002717833564 402 # 13 h -0.00000616830304 0.00000016851704 -0.00003103539353 403 # 14 c 0.00004884621752 -0.00002593547764 0.00004887866674 404 # 15 h -0.00001561780942 -0.00001607888272 -0.00002718198410 405 # 16 c -0.00000015329549 0.00004026946593 -0.00003204885735 406 # 17 h -0.00001239289077 0.00001308531162 0.00000092210161 407 # 18 h 0.00000616838338 0.00000016858902 0.00003103523777 408 # 19 c 0.00000362187006 0.00000466096401 0.00002733645543 409 # 20 h 0.00000645786975 -0.00000179869042 -0.00003328556910 410 # Optint: norm of internal gradient 0.00029710 411 # 412 if "CARTESIAN GRADIENTS" in line: 413 atomcoords = [] 414 atomnos = [] 415 line = next(inputfile) 416 while len(line.split()) == 5: 417 atomnos.append(self.periodic_table.number[line.split()[-4].capitalize()]) 418 atomcoords.append([utils.convertor(float(x), "bohr", "Angstrom") 419 for x in line.split()[-3:]]) 420 line = next(inputfile) 421 422 self.set_attribute('atomnos', atomnos) 423 self.set_attribute('natom', len(atomcoords)) 424 self.append_attribute("grads", atomcoords) 425 426 ## Atomic coordinates in job.last: 427 # +--------------------------------------------------+ 428 # | Atomic coordinate, charge and isotop information | 429 # +--------------------------------------------------+ 430 # 431 # 432 # atomic coordinates atom shells charge pseudo isotop 433 # -2.69176330 -0.00007129 -0.44712612 c 3 6.000 0 0 434 # -1.69851645 -0.00007332 2.06488947 c 3 6.000 0 0 435 # 0.92683848 -0.00007460 2.49592179 c 3 6.000 0 0 436 # 2.69176331 -0.00007127 0.44712612 c 3 6.000 0 0 437 # 1.69851645 -0.00007331 -2.06488947 c 3 6.000 0 0 438 #... 439 # -7.04373606 0.00092244 2.74543891 h 1 1.000 0 0 440 # -9.36352819 0.00017229 0.07445322 h 1 1.000 0 0 441 # -0.92683849 -0.00007461 -2.49592179 c 3 6.000 0 0 442 # -1.65164853 -0.00009927 -4.45456858 h 1 1.000 0 0 443 # 444 # During a jobex driven optimisation, the above geometry section will be printed at the start of each module (subprogram). 445 # We only want one coord entry per step, so we need to be careful not to add the same geometry multiple times. 446 # 447 # For 'large' molecules (probably >= 100 atoms), certain turbomole modules (eg, STATP) will only print a truncated 448 # coordinate section here, with some lines being cut and replaced with 3 lines consisting only of whitespace and 449 # full stops (.): 450 # 5.98651205 -2.33496770 -0.80565619 c 6.000 0 451 # 9.67800538 -0.33507410 1.19593467 c 6.000 0 452 # 6.76781184 0.69969618 5.52907043 c 6.000 0 453 # . . . . . . 454 # . . . . . . 455 # . . . . . . 456 # 21.85796486 5.25815448 5.52575969 h 1.000 0 457 # 16.70506498 7.08657699 13.19684279 h 1.000 0 458 # 25.54043040 7.25266316 7.52378680 h 1.000 0 459 # 460 # This is very annoying, but not disastrous as it is likely the coordinates from this step will have been already printed 461 # from another module, so we can just ignore sections like this. 462 try: 463 if 'Atomic coordinate, charge and isotop information' in line: 464 while 'atomic coordinates' not in line: 465 line = next(inputfile) 466 467 atomcoords = [] 468 atomnos = [] 469 line = next(inputfile) 470 while len(line.split()) >= 3: 471 atomnos.append(self.periodic_table.number[line.split()[3].capitalize()]) 472 atomcoords.append([utils.convertor(float(x), "bohr", "Angstrom") 473 for x in line.split()[:3]]) 474 line = next(inputfile) 475 476 # Check the coordinates we just parsed are not already in atomcoords. 477 if not hasattr(self, 'atomcoords') or self.atomcoords[-1] != atomcoords: 478 self.append_attribute('atomcoords', atomcoords) 479 self.set_attribute('atomnos', atomnos) 480 self.set_attribute('natom', len(atomcoords)) 481 except Exception as e: 482 # Something went wrong, check to see if the coord line we parsed was garbage, or if something else happened. 483 if set(line.strip()).issubset({' ', '.'}): 484 # This line is a truncated coord line, we can skip this coord section. 485 pass 486 else: 487 # Panic. 488 raise e 489 490 # Flag that indicates we are doing an opt. 491 if "OPTIMIZATION CYCLE" in line: 492 self.append_attribute("optstatus", data.ccData.OPT_UNKNOWN) 493 494 if "OPTIMIZATION CYCLE 1" in line: 495 # This is the start of the opt. 496 self.optstatus[-1] += data.ccData.OPT_NEW 497 498 # Optimisation convergence criteria using statpt. 499 # 500 # ****************************************************************** 501 # CONVERGENCE INFORMATION 502 # 503 # Converged? Value Criterion 504 # Energy change no 0.0000011 0.0000010 505 # RMS of displacement yes 0.0001152 0.0005000 506 # RMS of gradient yes 0.0000548 0.0005000 507 # MAX displacement yes 0.0001409 0.0010000 508 # MAX gradient yes 0.0000670 0.0010000 509 # ****************************************************************** 510 if "CONVERGENCE INFORMATION" in line: 511 # This is an optimisation. 512 # Skip lines. 513 line = next(inputfile) 514 line = next(inputfile) 515 line = next(inputfile) 516 517 convergence = [] 518 geovalues = [] 519 geotargets = [] 520 521 # There are a variable number of criteria. 522 while len(line.split()) > 3: 523 parts = line.split() 524 # lower() is for (unnecessary?) future proofing... 525 converged = parts[-3].lower() == "yes" 526 value = float(parts[-2]) 527 criterion = float(parts[-1]) 528 529 # TODO: possibly require some unit conversions? 530 convergence.append(converged) 531 geovalues.append(value) 532 geotargets.append(criterion) 533 534 # Next. 535 line = next(inputfile) 536 537 self.set_attribute("geotargets", geotargets) 538 self.append_attribute("geovalues", geovalues) 539 540 if all(convergence): 541 # This iteration has converged. 542 self.append_attribute("optdone", len(self.geovalues) -1) 543 self.optstatus[-1] += data.ccData.OPT_DONE 544 else: 545 # Not converged. 546 if not hasattr(self, 'optdone'): 547 self.set_attribute("optdone", []) 548 #self.optstatus[-1] += data.ccData.OPT_UNCONVERGED 549 550 551 # Frequency values in aoforce.out 552 # mode 7 8 9 10 11 12 553 # 554 # frequency 53.33 88.32 146.85 171.70 251.75 289.44 555 # 556 # symmetry a a a a a a 557 # 558 # IR YES YES YES YES YES YES 559 # |dDIP/dQ| (a.u.) 0.0002 0.0000 0.0005 0.0004 0.0000 0.0000 560 # intensity (km/mol) 0.05 0.00 0.39 0.28 0.00 0.00 561 # intensity ( % ) 0.05 0.00 0.40 0.28 0.00 0.00 562 # 563 # RAMAN YES YES YES YES YES YES 564 # 565 # 1 c x 0.00000 0.00001 0.00000 -0.01968 -0.04257 0.00001 566 # y -0.08246 -0.08792 0.02675 -0.00010 0.00000 0.17930 567 # z 0.00001 0.00003 0.00004 -0.10350 0.11992 -0.00003 568 # .... 569 # 570 # reduced mass(g/mol) 3.315 2.518 2.061 3.358 3.191 2.323 571 572 if 'NORMAL MODES and VIBRATIONAL FREQUENCIES (cm**(-1))' in line: 573 has_raman = False 574 while line[7:11] != 'mode': 575 line = next(inputfile) 576 if line.startswith(" differential RAMAN cross sections"): 577 has_raman = True 578 vibfreqs, vibsyms, vibirs, vibdisps, vibrmasses = [], [], [], [], [] 579 while 'all done ****' not in line: 580 581 if line.strip().startswith('mode'): 582 self.skip_line(inputfile, 'b') 583 line = next(inputfile) 584 assert line.strip().startswith('frequency') 585 freqs = [float(i.replace('i', '-')) for i in line.split()[1:]] 586 vibfreqs.extend(freqs) 587 self.skip_lines(inputfile, ['b']) 588 line = next(inputfile) 589 assert line.strip().startswith('symmetry') 590 syms = [self.normalisesym(sym) for sym in line.split()[1:]] 591 vibsyms.extend(syms) 592 593 self.skip_lines(inputfile, ['b', 'IR', 'dDIP/dQ']) 594 line = next(inputfile) 595 assert line.strip().startswith('intensity (km/mol)') 596 irs = [utils.float(f) for f in line.split()[2:]] 597 vibirs.extend(irs) 598 599 self.skip_lines(inputfile, ['intensity %', 'b', 'RAMAN']) 600 if has_raman: 601 self.skip_lines( 602 inputfile, 603 ['(par,par)', '(ort,ort)', '(ort,unpol)', 'depol. ratio'] 604 ) 605 line = next(inputfile) 606 assert not line.strip() 607 line = next(inputfile) 608 x, y, z = [], [], [] 609 atomcounter = 0 610 while atomcounter < self.natom: 611 x.append([float(i) for i in line.split()[3:]]) 612 line = next(inputfile) 613 y.append([float(i) for i in line.split()[1:]]) 614 line = next(inputfile) 615 z.append([float(i) for i in line.split()[1:]]) 616 line = next(inputfile) 617 atomcounter += 1 618 619 for j in range(len(x[0])): 620 disps = [] 621 for i in range(len(x)): 622 disps.append([x[i][j], y[i][j], z[i][j]]) 623 vibdisps.append(disps) 624 625 line = next(inputfile) 626 assert line.startswith('reduced mass(g/mol)') 627 rmasses = [utils.float(f) for f in line.split()[2:]] 628 vibrmasses.extend(rmasses) 629 630 if "zero point VIBRATIONAL energy" in line: 631 self.set_attribute("zpve", float(line.split()[6])) 632 633 line = next(inputfile) 634 635 self.set_attribute('vibfreqs', vibfreqs) 636 self.set_attribute('vibsyms', vibsyms) 637 self.set_attribute('vibirs', vibirs) 638 self.set_attribute('vibdisps', vibdisps) 639 self.set_attribute('vibrmasses', vibrmasses) 640 641 # In this section we are parsing mocoeffs and moenergies from 642 # the files like: mos, alpha and beta. 643 # $scfmo scfconv=6 format(4d20.14) 644 # # SCF total energy is -382.3457535740 a.u. 645 # # 646 # 1 a eigenvalue=-.97461484059799D+01 nsaos=60 647 # 0.69876828353937D+000.32405121159405D-010.87670894913921D-03-.85232349313288D-07 648 # 0.19361534257922D-04-.23841194890166D-01-.81711001390807D-020.13626356942047D-02 649 # ... 650 # ... 651 # $end 652 if (line.startswith('$scfmo') or line.startswith('$uhfmo')) and line.find('scfconv') > 0: 653 if line.strip().startswith('$uhfmo_alpha'): 654 self.unrestricted = True 655 656 # Need to skip the first line to start with lines starting with '#'. 657 line = next(inputfile) 658 while line.strip().startswith('#') and not line.find('eigenvalue') > 0: 659 line = next(inputfile) 660 661 moirreps = [] 662 moenergies = [] 663 mocoeffs = [] 664 mosyms = [] 665 666 while not line.strip().startswith('$'): 667 number, sym = line.split()[:2] 668 number = int(number) 669 sym = self.normalisesym(sym) 670 671 info = re.match(r".*eigenvalue=(?P<moenergy>[0-9D\.+-]{20})\s+nsaos=(?P<count>\d+).*", line) 672 eigenvalue = utils.float(info.group('moenergy')) 673 orbital_energy = utils.convertor(eigenvalue, 'hartree', 'eV') 674 675 moenergies.append(orbital_energy) 676 mosyms.append(sym) 677 moirreps.append((number, sym)) 678 679 single_coeffs = [] 680 nsaos = int(info.group('count')) 681 682 while(len(single_coeffs) < nsaos): 683 line = next(inputfile) 684 single_coeffs.extend(Turbomole.split_molines(line)) 685 686 mocoeffs.append(single_coeffs) 687 line = next(inputfile) 688 689 max_nsaos = max([len(i) for i in mocoeffs]) 690 for i in mocoeffs: 691 while len(i) < max_nsaos: 692 i.append(numpy.nan) 693 694 # We now need to sort our orbitals (because Turbomole groups them by symm). 695 mos = list(zip(moenergies, mocoeffs, mosyms, moirreps)) 696 mos.sort(key = lambda mo: mo[0]) 697 moenergies, mocoeffs, mosyms, moirreps = zip(*mos) 698 699 # MO irreps is not actually recognised as a cclib attribute, 700 # but we may need this info to parse other sections. 701 self.append_attribute("moirreps", moirreps) 702 self.append_attribute("moenergies", moenergies) 703 self.append_attribute("mocoeffs", mocoeffs) 704 self.append_attribute("mosyms", mosyms) 705 self.set_attribute("nmo", len(moenergies)) 706 707 # Parsing the scfenergies, scfvalues and scftargets from job.last file. 708 # scf convergence criterion : increment of total energy < .1000000D-05 709 # and increment of one-electron energy < .1000000D-02 710 # 711 # ... 712 # ... 713 # current damping : 0.700 714 # ITERATION ENERGY 1e-ENERGY 2e-ENERGY NORM[dD(SAO)] TOL 715 # 1 -382.34543727790 -1396.8009423 570.56292464 0.000D+00 0.556D-09 716 # Exc = -57.835278090846 N = 69.997494722 717 # max. resid. norm for Fia-block= 2.782D-05 for orbital 33a 718 # ... 719 # ... 720 # current damping : 0.750 721 # ITERATION ENERGY 1e-ENERGY 2e-ENERGY NORM[dD(SAO)] TOL 722 # 3 -382.34575357399 -1396.8009739 570.56263988 0.117D-03 0.319D-09 723 # Exc = -57.835593208072 N = 69.999813370 724 # max. resid. norm for Fia-block= 7.932D-06 for orbital 33a 725 # max. resid. fock norm = 8.105D-06 for orbital 33a 726 # 727 # convergence criteria satisfied after 3 iterations 728 # 729 # 730 # ------------------------------------------ 731 # | total energy = -382.34575357399 | 732 # ------------------------------------------ 733 # : kinetic energy = 375.67398458525 : 734 # : potential energy = -758.01973815924 : 735 # : virial theorem = 1.98255043001 : 736 # : wavefunction norm = 1.00000000000 : 737 # .......................................... 738 if 'scf convergence criterion' in line: 739 total_energy_threshold = utils.float(line.split()[-1]) 740 one_electron_energy_threshold = utils.float(next(inputfile).split()[-1]) 741 scftargets = [total_energy_threshold, one_electron_energy_threshold] 742 self.append_attribute('scftargets', scftargets) 743 744 # Get ready for SCF iteration energies. 745 self.iter_energy = [] 746 self.iter_one_elec_energy = [] 747 748 if 'ITERATION ENERGY' in line: 749 while 'convergence criteria satisfied' not in line: 750 # nbasis 751 # Doesn't appear to do anything, number of basis functions does not appear in this section? 752 #if "number of basis functions" in line: 753 # self.set_attribute("nbasis", int(line.split()[-1])) 754 755 if 'ITERATION ENERGY' in line: 756 line = next(inputfile) 757 info = line.split() 758 self.iter_energy.append(utils.float(info[1])) 759 self.iter_one_elec_energy.append(utils.float(info[2])) 760 line = next(inputfile) 761 762 assert len(self.iter_energy) == len(self.iter_one_elec_energy), \ 763 'Different number of values found for total energy and one electron energy.' 764 scfvalues = [[x - y, a - b] for x, y, a, b in 765 zip(self.iter_energy[1:], self.iter_energy[:-1], self.iter_one_elec_energy[1:], self.iter_one_elec_energy[:-1])] 766 self.append_attribute('scfvalues', scfvalues) 767 while 'total energy' not in line: 768 line = next(inputfile) 769 770 scfenergy = utils.convertor(utils.float(line.split()[4]), 'hartree', 'eV') 771 self.append_attribute('scfenergies', scfenergy) 772 773 # We need to determine whether this is a HF or DFT energy for metadata. 774 if self.is_DFT: 775 self.metadata['methods'].append("DFT") 776 else: 777 self.metadata['methods'].append("HF") 778 779 # ********************************************************************** 780 # * * 781 # * RHF energy : -74.9644564256 * 782 # * MP2 correlation energy (doubles) : -0.0365225363 * 783 # * * 784 # * Final MP2 energy : -75.0009789619 * 785 # ... 786 # * Norm of MP1 T2 amplitudes : 0.0673494687 * 787 # * * 788 # ********************************************************************** 789 # OR 790 # ********************************************************************** 791 # * * 792 # * RHF energy : -74.9644564256 * 793 # * correlation energy : -0.0507799360 * 794 # * * 795 # * Final CCSD energy : -75.0152363616 * 796 # * * 797 # * D1 diagnostic : 0.0132 * 798 # * * 799 # ********************************************************************** 800 if 'Final CCSD energy' in line or 'Final CC2 energy' in line: 801 self.append_attribute( 802 'ccenergies', 803 utils.convertor(utils.float(line.split()[5]), 'hartree', 'eV') 804 ) 805 self.metadata['methods'].append("CCSD") 806 807 # Look for MP energies. 808 for mp_level in range(2,6): 809 if "Final MP{} energy".format(mp_level) in line: 810 mpenergy = utils.convertor(utils.float(line.split()[5]), 'hartree', 'eV') 811 if mp_level == 2: 812 self.append_attribute('mpenergies', [mpenergy]) 813 else: 814 self.mpenergies[-1].append(mpenergy) 815 self.metadata['methods'].append("MP{}".format(mp_level)) 816 817 # ***************************************************** 818 # * * 819 # * SCF-energy : -74.49827196840999 * 820 # * MP2-energy : -0.19254365976227 * 821 # * total : -74.69081562817226 * 822 # * * 823 # * (MP2-energy evaluated from T2 amplitudes) * 824 # * * 825 # ***************************************************** 826 if 'MP2-energy' in line: 827 line = next(inputfile) 828 if 'total' in line: 829 mp2energy = [utils.convertor(utils.float(line.split()[3]), 'hartree', 'eV')] 830 self.append_attribute('mpenergies', mp2energy) 831 self.metadata['methods'].append("MP2") 832 833 # Support for the now outdated (?) rimp2 834 # ------------------------------------------------ 835 # Method : MP2 836 # Total Energy : -75.0009789796 837 # ------------------------------------------------ 838 if "Method : MP2" in line: 839 line = next(inputfile) 840 mp2energy = [utils.convertor(utils.float(line.split()[3]), 'hartree', 'eV')] 841 self.append_attribute('mpenergies', mp2energy) 842 self.metadata['methods'].append("MP2") 843 844 845 # Excited state info from escf. 846 # 1 singlet a excitation 847 # 848 # 849 # Total energy: -112.9060086800507 850 # 851 # Excitation energy: 0.3111453714493050 852 # 853 # Excitation energy / eV: 8.466699968968847 854 # 855 # Excitation energy / nm: 146.4375074832943 856 # 857 # Excitation energy / cm^(-1): 68288.51549285055 858 # 859 # 860 # Oscillator strength: 861 # 862 # velocity representation: 0.1011193926229111 863 # 864 # length representation: 0.9461905579062439E-01 865 # 866 # mixed representation: 0.9781524141002400E-01 867 # 868 # 869 # Rotatory strength: 870 # 871 # velocity representation: 0.1282566570726007E-10 872 # 873 # velocity rep. / 10^(-40)erg*cm^3: 0.8285997783054736E-06 874 # 875 # length representation: 0.1242930160103103E-10 876 # 877 # length rep. / 10^(-40)erg*cm^3: 0.8029927479925187E-06 878 # 879 # 880 # Dominant contributions: 881 # 882 # occ. orbital energy / eV virt. orbital energy / eV |coeff.|^2*100 883 # 7 a -10.76 9 a -0.78 96.7 884 # 885 ## For UHF: 886 # occ. orbital energy / eV virt. orbital energy / eV |coeff.|^2*100 887 # 7 a alpha -15.12 12 a alpha 4.74 24.5 888 # 7 a beta -15.12 12 a beta 4.74 24.5 889 if "Excitation energy:" in line and "excitation" in self.last_lines[-5].split(): 890 # The irrep of the state is a few lines back. 891 symm_parts = self.last_lines[-5].split() 892 # We don't always have the multiplicity available. 893 if len(symm_parts) < 4: 894 # No mult. 895 mult = "???" 896 else: 897 mult = symm_parts[1].capitalize() 898 899 symmetry = "{}-{}".format(mult, symm_parts[-2].capitalize()) 900 self.append_attribute("etsyms", symmetry) 901 902 # Energy should be in cm-1... 903 energy = utils.convertor(utils.float(line.split()[-1]), 'hartree', 'wavenumber') 904 self.append_attribute("etenergies", energy) 905 906 # Oscillator strength. 907 while "length representation:" not in line: 908 line = next(inputfile) 909 oscillator_strength = utils.float(line.split()[-1]) 910 self.append_attribute("etoscs", oscillator_strength) 911 912 line = next(inputfile) 913 914 # Rotatory strength. 915 while "length representation:" not in line: 916 line = next(inputfile) 917 rotatory_strength = utils.float(line.split()[-1]) 918 self.append_attribute("etrotats", rotatory_strength) 919 920 while "Dominant contributions:" not in line: 921 line = next(inputfile) 922 line = next(inputfile) 923 line = next(inputfile) 924 line = next(inputfile) 925 926 # We can't get transitions if we don't have any orbitals. 927 if hasattr(self, 'moirreps'): 928 transitions = [] 929 930 while len(line.split()) > 0: 931 parts = line.split() 932 933 # Get alpha/beta, deleting "alpha" or "beta" from the line so 934 # our indexes are aligned for both RHF and UHF. 935 if parts[2] == "alpha": 936 start_AB = 0 937 parts.pop(2) 938 if parts[2] == "beta": 939 start_AB = 1 940 parts.pop(2) 941 else: 942 start_AB = 0 943 944 # Determine our start orbital. 945 start_MO_irrep = (int(parts[0]), self.normalisesym(parts[1])) 946 start_MO = self.moirreps[start_AB].index(start_MO_irrep) 947 948 if parts[5] == "beta": 949 end_AB = 1 950 else: 951 end_AB = 0 952 953 # And end orbital. 954 end_MO_irrep = (int(parts[3]), self.normalisesym(parts[4])) 955 end_MO = self.moirreps[end_AB].index(end_MO_irrep) 956 957 # Finally, get our coefficient. 958 # Sadly, Turbomole only prints |coeff.|^2*100, so we can't determine whether 959 # the coefficient is negative or positive... 960 coeff = (utils.float(parts[-1]) /100) ** 0.5 961 962 # Add to list. 963 transitions.append(( 964 (start_MO, start_AB), 965 (end_MO, end_AB), 966 coeff 967 )) 968 969 # Go again. 970 line = next(inputfile) 971 972 # Print a warning because we're missing the sign (+/-) of etsecs. 973 self.logger.warning("The sign (+/-) of the coefficient for singly-excited configurations is not available") 974 self.append_attribute("etsecs", transitions) 975 976 # Transition dipoles. 977 while "Electric transition dipole moment (length rep.):" not in line: 978 line = next(inputfile) 979 line = next(inputfile) 980 981 # Unlike PDM, we leave TDMs in units of au (based on the implementation for Gaussian). 982 line = next(inputfile) 983 tdm_x = float(line.split()[1]) 984 line = next(inputfile) 985 tdm_y = float(line.split()[1]) 986 line = next(inputfile) 987 tdm_z = float(line.split()[1]) 988 989 self.append_attribute("etdips", [tdm_x, tdm_y, tdm_z]) 990 991 # Excitation energies with ricc2. 992 # +================================================================================+ 993 # | sym | multi | state | CC2 excitation energies | %t1 | %t2 | 994 # | | | +----------------------------------------+--------+--------+ 995 # | | | | Hartree | eV | cm-1 | % | % | 996 # +================================================================================+ 997 # | a | 1 | 1 | 0.3257471 | 8.86403 | 71493.234 | 94.89 | 5.11 | 998 # | a | 1 | 2 | 0.3257471 | 8.86403 | 71493.232 | 94.89 | 5.11 | 999 # | a | 1 | 3 | 0.3864541 | 10.51595 | 84816.880 | 95.05 | 4.95 | 1000 # | a | 1 | 4 | 0.3938325 | 10.71673 | 86436.241 | 94.09 | 5.91 | 1001 # | a | 1 | 5 | 0.3938325 | 10.71673 | 86436.241 | 94.09 | 5.91 | 1002 # | a | 1 | 6 | 0.4122171 | 11.21700 | 90471.202 | 93.32 | 6.68 | 1003 # | a | 1 | 7 | 0.4351537 | 11.84114 | 95505.195 | 94.16 | 5.84 | 1004 # | a | 1 | 8 | 0.4454969 | 12.12259 | 97775.278 | 94.26 | 5.74 | 1005 # | a | 1 | 9 | 0.4454969 | 12.12259 | 97775.273 | 94.26 | 5.74 | 1006 # | a | 1 | 10 | 0.5036295 | 13.70446 | 110533.909 | 93.51 | 6.49 | 1007 # +================================================================================+ 1008 # TODO: Perhaps need a more general way to look for lines like this? 1009 if "| sym | multi | state | CC2 excitation energies | %t1 | %t2 |" in line \ 1010 or "| sym | multi | state | ADC(2) excitation energies | %t1 | %t2 |" in line: 1011 self.skip_lines(inputfile, ["divider", "units", "divider"]) 1012 line = next(inputfile) 1013 1014 # Reset in case we've parsed this section before for some reason... 1015 for attr in ("etenergies", "etsyms", "etoscs", "etsecs"): 1016 if hasattr(self, attr): 1017 delattr(self, attr) 1018 1019 while len(line.split("|")) > 1: 1020 # Split. 1021 parts = [part.strip() for part in line.split("|")] 1022 1023 if parts[2] == "1": 1024 mult = "Singlet" 1025 elif parts[2] == "2": 1026 mult = "Doublet" 1027 elif parts[2] == "3": 1028 mult = "Triplet" 1029 elif parts[2] == "4": 1030 mult = "Quartet" 1031 elif parts[2] == "-": 1032 mult = "???" 1033 else: 1034 mult = parts[2] 1035 1036 symmetry = "{}-{}".format(mult, parts[1].capitalize()) 1037 self.append_attribute("etsyms", symmetry) 1038 1039 energy = utils.float(parts[6]) 1040 self.append_attribute("etenergies", energy) 1041 1042 # Go again. 1043 line = next(inputfile) 1044 1045 # Singly excited configurations are printed separately with ricc2. 1046 # +=======================================================================+ 1047 # | type: RE0 symmetry: a state: 1 | 1048 # +-----------------------+-----------------------+-----------------------+ 1049 # | occ. orb. index spin | vir. orb. index spin | coeff/|amp| % | 1050 # +=======================+=======================+=======================+ 1051 # | 7 a 7 | 9 a 9 | 0.65683 43.1 | 1052 # | 7 a 7 | 13 a 13 | 0.47926 23.0 | 1053 # | 7 a 7 | 12 a 12 | -0.44814 20.1 | 1054 # | 7 a 7 | 10 a 10 | -0.24445 6.0 | 1055 # | 7 a 7 | 16 a 16 | -0.14810 2.2 | 1056 # | 4 a 4 | 13 a 13 | 0.11053 1.2 | 1057 # +=======================+=======================+=======================+ 1058 # norm of printed elements: 0.95585 1059 # for UHF: 1060 # +=======================================================================+ 1061 # | type: RE0 symmetry: a state: 1 | 1062 # +-----------------------+-----------------------+-----------------------+ 1063 # | occ. orb. index spin | vir. orb. index spin | coeff/|amp| % | 1064 # +=======================+=======================+=======================+ 1065 # | 7 a 7 (b) | 12 a 12 (b) | 0.49335 24.3 | 1066 # | 7 a 7 (a) | 12 a 12 (a) | -0.49335 24.3 | 1067 # | 7 a 7 (a) | 9 a 9 (a) | 0.42257 17.9 | 1068 # | 7 a 7 (b) | 9 a 9 (b) | -0.42257 17.9 | 1069 # | 7 a 7 (b) | 13 a 13 (b) | 0.20077 4.0 | 1070 # | 7 a 7 (a) | 13 a 13 (a) | -0.20077 4.0 | 1071 # | 7 a 7 (a) | 15 a 15 (a) | 0.11607 1.3 | 1072 # | 7 a 7 (b) | 15 a 15 (b) | -0.11607 1.3 | 1073 # +=======================+=======================+=======================+ 1074 if "| occ. orb. index spin | vir. orb. index spin | coeff/|amp| % |" in line: 1075 self.skip_lines(inputfile, ["divider"]) 1076 line = next(inputfile) 1077 1078 transitions = [] 1079 1080 while len(line.split()) > 1: 1081 parts = line.split() 1082 1083 # Determine our start orbital. 1084 # Unlike the HF/DFT level excited states printed by escf, ricc2 prints 1085 # the orbital index directly. 1086 start_MO = int(parts[3]) -1 1087 1088 # Get alpha/beta, deleting "alpha" or "beta" from the line so our 1089 # indexes are aligned for both RHF and UHF. 1090 if parts[4] == "(a)": 1091 start_AB = 0 1092 parts.pop(4) 1093 if parts[4] == "(b)": 1094 start_AB = 1 1095 parts.pop(4) 1096 else: 1097 start_AB = 0 1098 1099 # And end orbital. 1100 end_MO = int(parts[7]) -1 1101 1102 if parts[8] == "(b)": 1103 end_AB = 1 1104 else: 1105 end_AB = 0 1106 1107 # Finally, get our coefficient. 1108 # Once again, ricc2 beats escf and we have both the coefficient and % 1109 # available. 1110 coeff = utils.float(parts[-3]) 1111 1112 # Add to list. 1113 transitions.append(( 1114 (start_MO, start_AB), 1115 (end_MO, end_AB), 1116 coeff 1117 )) 1118 1119 # Go again. 1120 line = next(inputfile) 1121 1122 self.append_attribute("etsecs", transitions) 1123 1124 # Oscillator strengths are also printed separately with ricc2, 1125 # and may be missing entirely. 1126 # 1127 # oscillator strength (length gauge) : 0.09614727 1128 # 1129 if "oscillator strength (length gauge) :" in line: 1130 self.append_attribute("etoscs", utils.float(line.split()[-1])) 1131 1132 1133 # All done for this loop. 1134 # Keep track of last lines. 1135 self.last_lines.append(line) 1136 1137 if ": all done ****" in line: 1138 # End of module, set success flag. 1139 self.metadata['success'] = True 1140 1141 1142 1143 def split_irrep(self, irrep): 1144 """ 1145 Split a Turbomole irrep into number and symmetry. 1146 """ 1147 rematch = re.match(r"^([0-9]+)(.+)$", irrep) 1148 number = int(rematch.group(1)) 1149 sym = self.normalisesym(rematch.group(2)) 1150 return (number, sym) 1151 1152 def parse_dscf_orbitals(self, inputfile, line): 1153 """ 1154 Extract orbital occupation and energies from a dscf logfile. 1155 1156 Returns 1157 ------- 1158 tuple 1159 a two membered tuple where the first element is a list of dictionaries of the the orbitals parsed, while the second is the line on which parsing should continue. 1160 """ 1161 ## Orbital occupation info from dscf. 1162 # orbitals $scfmo will be written to file mos 1163 # 1164 # irrep 1a 2a 3a 4a 5a 1165 # eigenvalues H -20.25992 -1.24314 -0.57053 -0.46144 -0.39295 1166 # eV -551.3047 -33.8279 -15.5250 -12.5564 -10.6929 1167 # occupation 2.0000 2.0000 2.0000 2.0000 2.0000 1168 # ... 1169 # irrep 6a 7a 1170 # eigenvalues H 0.55091 0.64409 1171 # eV 14.9910 17.5268 1172 ## Or 1173 # orbitals $uhfmo_beta will be written to file beta 1174 # 1175 # orbitals $uhfmo_alpha will be written to file alpha 1176 # 1177 # alpha: 1178 # 1179 # irrep 31a 32a 33a 34a 35a 1180 # eigenvalues H -0.47570 -0.46573 -0.40741 -0.39213 -0.35411 1181 # eV -12.9446 -12.6733 -11.0862 -10.6705 -9.6358 1182 # occupation 1.0000 1.0000 1.0000 1.0000 1.0000 1183 # ... 1184 # irrep 36a 37a 38a 39a 40a 1185 # eigenvalues H -0.18634 -0.10035 -0.09666 -0.02740 0.06072 1186 # eV -5.0705 -2.7306 -2.6303 -0.7455 1.6522 1187 # 1188 # beta: 1189 # 1190 # irrep 30a 31a 32a 33a 34a 1191 # eigenvalues H -0.49118 -0.47348 -0.44470 -0.39020 -0.37919 1192 # eV -13.3658 -12.8842 -12.1009 -10.6181 -10.3184 1193 # occupation 1.0000 1.0000 1.0000 1.0000 1.0000 1194 # ... 1195 # irrep 35a 36a 37a 38a 39a 1196 # eigenvalues H -0.28091 -0.15088 -0.09343 -0.07531 -0.00688 1197 # eV -7.6440 -4.1058 -2.5424 -2.0493 -0.1873 1198 # Skip blank line. 1199 line = next(inputfile) 1200 1201 orbitals = [] 1202 while True: 1203 irreps = [] 1204 energies_hartree = [] 1205 energies_eV = [] 1206 occupations = [] 1207 1208 # MO index 1209 line = next(inputfile) 1210 1211 # Check we're still in the right section. 1212 if "irrep" not in line: 1213 # All done. 1214 break 1215 else: 1216 # Turbomole lists orbitals of different symmetry separately. 1217 irreps = line.split()[1:] 1218 1219 # Energy in H. 1220 line = next(inputfile) 1221 energies_hartree = [float(energy) for energy in line.split()[2:]] 1222 1223 # Energy in eV. 1224 line = next(inputfile) 1225 energies_eV = [float(energy) for energy in line.split()[1:]] 1226 1227 # Occupation. 1228 # This line will be missing if the orbitals are virtual (unoccupied). 1229 line = next(inputfile) 1230 if "occupation" in line: 1231 occupations = [float(occupation) for occupation in line.split()[1:]] 1232 line = next(inputfile) 1233 1234 # If we have any missing occupations, fill with 0 1235 occupations.extend([0.0] * (len(irreps) - len(occupations))) 1236 1237 # Add to list. 1238 orbitals.extend([ 1239 {'irrep': irrep, 'energy_H': energy_H, 'energy_eV': energy_eV, 'occupancy': occupation} 1240 for irrep, energy_H, energy_eV, occupation 1241 in zip(irreps, energies_hartree, energies_eV, occupations) 1242 ]) 1243 1244 return orbitals, line 1245 1246 1247 def determine_homo(self, mosyms, dscf_mos): 1248 """ 1249 Determine the highest occupied molecular orbital. 1250 1251 Returns 1252 ------- 1253 int 1254 the index of the HOMO. 1255 """ 1256 # First, we need a mapping between each irrep and its index in mosyms etc. 1257 symm_count = {} 1258 irreps = [] 1259 1260 for symm in mosyms: 1261 try: 1262 symm_count[symm] += 1 1263 except KeyError: 1264 symm_count[symm] = 1 1265 1266 irreps.append((symm_count[symm], symm)) 1267 1268 # We also have occupancy info from dscf (but probably only for those orbitals close to the HOMO/LUMO gap). 1269 # We now need to determine which orbital with occupancy is highest. 1270 homo = 0 1271 for mo in dscf_mos: 1272 if mo['occupancy'] > 0: 1273 # This orbital has electrons, determine its position out of all orbitals. 1274 index = irreps.index(self.split_irrep(mo['irrep'])) 1275 if index > homo: 1276 homo = index 1277 1278 return homo 1279 1280 def deleting_modes(self, vibfreqs, vibdisps, vibirs, vibrmasses): 1281 """Deleting frequencies relating to translations or rotations""" 1282 i = 0 1283 while i < len(vibfreqs): 1284 if vibfreqs[i] == 0.0: 1285 # Deleting frequencies that have value 0 since they 1286 # do not correspond to vibrations. 1287 del vibfreqs[i], vibdisps[i], vibirs[i], vibrmasses[i] 1288 i -= 1 1289 i += 1 1290 1291 def after_parsing(self): 1292 if hasattr(self, 'vibfreqs'): 1293 self.deleting_modes(self.vibfreqs, self.vibdisps, self.vibirs, self.vibrmasses) 1294 1295 # Try and determine our multiplicity from our orbitals. 1296 if hasattr(self, "homos"): 1297 if len(self.homos) == 1: 1298 # If we are restricted (len(homos) == 1); assume we have to be singlet. 1299 self.set_attribute("mult", 1) 1300 else: 1301 # Unrestricted, the difference in homos should tell us the no. of unpaired e-. 1302 unpaired_e = abs(self.homos[0] - self.homos[1]) 1303 self.set_attribute("mult", unpaired_e +1) 1304 1305 # Set a flag if we stopped part way through an opt. 1306 if hasattr(self, "optstatus") and self.optstatus[-1] != data.ccData.OPT_DONE: 1307 self.optstatus[-1] += data.ccData.OPT_UNCONVERGED 1308 1309 1310class OldTurbomole(logfileparser.Logfile): 1311 """A Turbomole output file. Code is outdated and is not being used.""" 1312 1313 def __init__(self, *args): 1314 # Call the __init__ method of the superclass 1315 super(Turbomole, self).__init__(logname="Turbomole", *args) 1316 1317 def __str__(self): 1318 """Return a string representation of the object.""" 1319 return "Turbomole output file %s" % (self.filename) 1320 1321 def __repr__(self): 1322 """Return a representation of the object.""" 1323 return 'Turbomole("%s")' % (self.filename) 1324 1325 def atlist(self, atstr): 1326 # turn atstr from atoms section into array 1327 1328 fields=atstr.split(',') 1329 list=[] 1330 for f in fields: 1331 if(f.find('-')!=-1): 1332 rangefields=f.split('-') 1333 start=int(rangefields[0]) 1334 end=int(rangefields[1]) 1335 1336 for j in range(start, end+1, 1): 1337 list.append(j-1) 1338 else: 1339 list.append(int(f)-1) 1340 return(list) 1341 1342 def normalisesym(self, label): 1343 """Normalise the symmetries used by Turbomole.""" 1344 return ans 1345 1346 def before_parsing(self): 1347 self.geoopt = False # Is this a GeoOpt? Needed for SCF targets/values. 1348 1349 def split_molines(self, inline): 1350 line=inline.replace("D", "E") 1351 f1=line[0:20] 1352 f2=line[20:40] 1353 f3=line[40:60] 1354 f4=line[60:80] 1355 1356 if(len(f4)>1): 1357 return( (float(f1), float(f2), float(f3), float(f4)) ) 1358 if(len(f3)>1): 1359 return( (float(f1), float(f2), float(f3)) ) 1360 if(len(f2)>1): 1361 return( (float(f1), float(f2)) ) 1362 if(len(f1)>1): 1363 return([float(f1)]) 1364 return 1365 1366 def extract(self, inputfile, line): 1367 """Extract information from the file object inputfile.""" 1368 1369 if line[3:11]=="nbf(AO)=": 1370 nmo=int(line[11:]) 1371 self.nbasis=nmo 1372 self.nmo=nmo 1373 if line[3:9]=="nshell": 1374 temp=line.split('=') 1375 homos=int(temp[1]) 1376 1377 if line[0:6] == "$basis": 1378 print("Found basis") 1379 self.basis_lib=[] 1380 line = inputfile.next() 1381 line = inputfile.next() 1382 1383 while line[0] != '*' and line[0] != '$': 1384 temp=line.split() 1385 line = inputfile.next() 1386 while line[0]=="#": 1387 line = inputfile.next() 1388 self.basis_lib.append(AtomBasis(temp[0], temp[1], inputfile)) 1389 line = inputfile.next() 1390 if line == "$ecp\n": 1391 self.ecp_lib=[] 1392 1393 line = inputfile.next() 1394 line = inputfile.next() 1395 1396 while line[0] != '*' and line[0] != '$': 1397 fields=line.split() 1398 atname=fields[0] 1399 ecpname=fields[1] 1400 line = inputfile.next() 1401 line = inputfile.next() 1402 fields=line.split() 1403 ncore = int(fields[2]) 1404 1405 while line[0] != '*': 1406 line = inputfile.next() 1407 self.ecp_lib.append([atname, ecpname, ncore]) 1408 1409 if line[0:6] == "$coord": 1410 if line[0:11] == "$coordinate": 1411# print "Breaking" 1412 return 1413 1414# print "Found coords" 1415 self.atomcoords = [] 1416 self.atomnos = [] 1417 atomcoords = [] 1418 atomnos = [] 1419 1420 line = inputfile.next() 1421 if line[0:5] == "$user": 1422# print "Breaking" 1423 return 1424 1425 while line[0] != "$": 1426 temp = line.split() 1427 atsym=temp[3].capitalize() 1428 atomnos.append(self.table.number[atsym]) 1429 atomcoords.append([utils.convertor(float(x), "bohr", "Angstrom") 1430 for x in temp[0:3]]) 1431 line = inputfile.next() 1432 self.atomcoords.append(atomcoords) 1433 self.atomnos = numpy.array(atomnos, "i") 1434 1435 if line[14:32] == "atomic coordinates": 1436 atomcoords = [] 1437 atomnos = [] 1438 1439 line = inputfile.next() 1440 1441 while len(line) > 2: 1442 temp = line.split() 1443 atsym = temp[3].capitalize() 1444 atomnos.append(self.table.number[atsym]) 1445 atomcoords.append([utils.convertor(float(x), "bohr", "Angstrom") 1446 for x in temp[0:3]]) 1447 line = inputfile.next() 1448 1449 if not hasattr(self,"atomcoords"): 1450 self.atomcoords = [] 1451 1452 self.atomcoords.append(atomcoords) 1453 self.atomnos = numpy.array(atomnos, "i") 1454 1455 if line[0:6] == "$atoms": 1456 print("parsing atoms") 1457 line = inputfile.next() 1458 self.atomlist=[] 1459 while line[0]!="$": 1460 temp=line.split() 1461 at=temp[0] 1462 atnosstr=temp[1] 1463 while atnosstr[-1] == ",": 1464 line = inputfile.next() 1465 temp=line.split() 1466 atnosstr=atnosstr+temp[0] 1467# print "Debug:", atnosstr 1468 atlist=self.atlist(atnosstr) 1469 1470 line = inputfile.next() 1471 1472 temp=line.split() 1473# print "Debug basisname (temp):",temp 1474 basisname=temp[2] 1475 ecpname='' 1476 line = inputfile.next() 1477 while(line.find('jbas')!=-1 or line.find('ecp')!=-1 or 1478 line.find('jkbas')!=-1): 1479 if line.find('ecp')!=-1: 1480 temp=line.split() 1481 ecpname=temp[2] 1482 line = inputfile.next() 1483 1484 self.atomlist.append( (at, basisname, ecpname, atlist)) 1485 1486# I have no idea what this does, so "comment" out 1487 if line[3:10]=="natoms=": 1488# if 0: 1489 1490 self.natom=int(line[10:]) 1491 1492 basistable=[] 1493 1494 for i in range(0, self.natom, 1): 1495 for j in range(0, len(self.atomlist), 1): 1496 for k in range(0, len(self.atomlist[j][3]), 1): 1497 if self.atomlist[j][3][k]==i: 1498 basistable.append((self.atomlist[j][0], 1499 self.atomlist[j][1], 1500 self.atomlist[j][2])) 1501 self.aonames=[] 1502 counter=1 1503 for a, b, c in basistable: 1504 ncore=0 1505 if len(c) > 0: 1506 for i in range(0, len(self.ecp_lib), 1): 1507 if self.ecp_lib[i][0]==a and \ 1508 self.ecp_lib[i][1]==c: 1509 ncore=self.ecp_lib[i][2] 1510 1511 for i in range(0, len(self.basis_lib), 1): 1512 if self.basis_lib[i].atname==a and self.basis_lib[i].basis_name==b: 1513 pa=a.capitalize() 1514 basis=self.basis_lib[i] 1515 1516 s_counter=1 1517 p_counter=2 1518 d_counter=3 1519 f_counter=4 1520 g_counter=5 1521# this is a really ugly piece of code to assign the right labels to 1522# basis functions on atoms with an ecp 1523 if ncore == 2: 1524 s_counter=2 1525 elif ncore == 10: 1526 s_counter=3 1527 p_counter=3 1528 elif ncore == 18: 1529 s_counter=4 1530 p_counter=4 1531 elif ncore == 28: 1532 s_counter=4 1533 p_counter=4 1534 d_counter=4 1535 elif ncore == 36: 1536 s_counter=5 1537 p_counter=5 1538 d_counter=5 1539 elif ncore == 46: 1540 s_counter=5 1541 p_counter=5 1542 d_counter=6 1543 1544 for j in range(0, len(basis.symmetries), 1): 1545 if basis.symmetries[j]=='s': 1546 self.aonames.append("%s%d_%d%s" % \ 1547 (pa, counter, s_counter, "S")) 1548 s_counter=s_counter+1 1549 elif basis.symmetries[j]=='p': 1550 self.aonames.append("%s%d_%d%s" % \ 1551 (pa, counter, p_counter, "PX")) 1552 self.aonames.append("%s%d_%d%s" % \ 1553 (pa, counter, p_counter, "PY")) 1554 self.aonames.append("%s%d_%d%s" % \ 1555 (pa, counter, p_counter, "PZ")) 1556 p_counter=p_counter+1 1557 elif basis.symmetries[j]=='d': 1558 self.aonames.append("%s%d_%d%s" % \ 1559 (pa, counter, d_counter, "D 0")) 1560 self.aonames.append("%s%d_%d%s" % \ 1561 (pa, counter, d_counter, "D+1")) 1562 self.aonames.append("%s%d_%d%s" % \ 1563 (pa, counter, d_counter, "D-1")) 1564 self.aonames.append("%s%d_%d%s" % \ 1565 (pa, counter, d_counter, "D+2")) 1566 self.aonames.append("%s%d_%d%s" % \ 1567 (pa, counter, d_counter, "D-2")) 1568 d_counter=d_counter+1 1569 elif basis.symmetries[j]=='f': 1570 self.aonames.append("%s%d_%d%s" % \ 1571 (pa, counter, f_counter, "F 0")) 1572 self.aonames.append("%s%d_%d%s" % \ 1573 (pa, counter, f_counter, "F+1")) 1574 self.aonames.append("%s%d_%d%s" % \ 1575 (pa, counter, f_counter, "F-1")) 1576 self.aonames.append("%s%d_%d%s" % \ 1577 (pa, counter, f_counter, "F+2")) 1578 self.aonames.append("%s%d_%d%s" % \ 1579 (pa, counter, f_counter, "F-2")) 1580 self.aonames.append("%s%d_%d%s" % \ 1581 (pa, counter, f_counter, "F+3")) 1582 self.aonames.append("%s%d_%d%s" % \ 1583 (pa, counter, f_counter, "F-3")) 1584 elif basis.symmetries[j]=='g': 1585 self.aonames.append("%s%d_%d%s" % \ 1586 (pa, counter, f_counter, "G 0")) 1587 self.aonames.append("%s%d_%d%s" % \ 1588 (pa, counter, f_counter, "G+1")) 1589 self.aonames.append("%s%d_%d%s" % \ 1590 (pa, counter, f_counter, "G-1")) 1591 self.aonames.append("%s%d_%d%s" % \ 1592 (pa, counter, g_counter, "G+2")) 1593 self.aonames.append("%s%d_%d%s" % \ 1594 (pa, counter, g_counter, "G-2")) 1595 self.aonames.append("%s%d_%d%s" % \ 1596 (pa, counter, g_counter, "G+3")) 1597 self.aonames.append("%s%d_%d%s" % \ 1598 (pa, counter, g_counter, "G-3")) 1599 self.aonames.append("%s%d_%d%s" % \ 1600 (pa, counter, g_counter, "G+4")) 1601 self.aonames.append("%s%d_%d%s" % \ 1602 (pa, counter, g_counter, "G-4")) 1603 break 1604 counter=counter+1 1605 1606 if line=="$closed shells\n": 1607 line = inputfile.next() 1608 temp = line.split() 1609 occs = int(temp[1][2:]) 1610 self.homos = numpy.array([occs-1], "i") 1611 1612 if line == "$alpha shells\n": 1613 line = inputfile.next() 1614 temp = line.split() 1615 occ_a = int(temp[1][2:]) 1616 line = inputfile.next() # should be $beta shells 1617 line = inputfile.next() # the beta occs 1618 temp = line.split() 1619 occ_b = int(temp[1][2:]) 1620 self.homos = numpy.array([occ_a-1,occ_b-1], "i") 1621 1622 if line[12:24]=="OVERLAP(CAO)": 1623 line = inputfile.next() 1624 line = inputfile.next() 1625 overlaparray=[] 1626 self.aooverlaps=numpy.zeros( (self.nbasis, self.nbasis), "d") 1627 while line != " ----------------------\n": 1628 temp=line.split() 1629 overlaparray.extend(map(float, temp)) 1630 line = inputfile.next() 1631 counter=0 1632 1633 for i in range(0, self.nbasis, 1): 1634 for j in range(0, i+1, 1): 1635 self.aooverlaps[i][j]=overlaparray[counter] 1636 self.aooverlaps[j][i]=overlaparray[counter] 1637 counter=counter+1 1638 1639 if ( line[0:6] == "$scfmo" or line[0:12] == "$uhfmo_alpha" ) and line.find("scf") > 0: 1640 temp = line.split() 1641 1642 if temp[1][0:7] == "scfdump": 1643# self.logger.warning("SCF not converged?") 1644 print("SCF not converged?!") 1645 1646 if line[0:12] == "$uhfmo_alpha": # if unrestricted, create flag saying so 1647 unrestricted = 1 1648 else: 1649 unrestricted = 0 1650 1651 self.moenergies=[] 1652 self.mocoeffs=[] 1653 1654 for spin in range(unrestricted + 1): # make sure we cover all instances 1655 title = inputfile.next() 1656 while(title[0] == "#"): 1657 title = inputfile.next() 1658 1659# mocoeffs = numpy.zeros((self.nbasis, self.nbasis), "d") 1660 moenergies = [] 1661 moarray=[] 1662 1663 if spin == 1 and title[0:11] == "$uhfmo_beta": 1664 title = inputfile.next() 1665 while title[0] == "#": 1666 title = inputfile.next() 1667 1668 while(title[0] != '$'): 1669 temp=title.split() 1670 1671 orb_symm=temp[1] 1672 1673 try: 1674 energy = float(temp[2][11:].replace("D", "E")) 1675 except ValueError: 1676 print(spin, ": ", title) 1677 1678 orb_en = utils.convertor(energy,"hartree","eV") 1679 1680 moenergies.append(orb_en) 1681 single_mo = [] 1682 1683 while(len(single_mo)<self.nbasis): 1684 self.updateprogress(inputfile, "Coefficients", self.cupdate) 1685 title = inputfile.next() 1686 lines_coeffs=self.split_molines(title) 1687 single_mo.extend(lines_coeffs) 1688 1689 moarray.append(single_mo) 1690 title = inputfile.next() 1691 1692# for i in range(0, len(moarray), 1): 1693# for j in range(0, self.nbasis, 1): 1694# try: 1695# mocoeffs[i][j]=moarray[i][j] 1696# except IndexError: 1697# print "Index Error in mocoeffs.", spin, i, j 1698# break 1699 1700 mocoeffs = numpy.array(moarray,"d") 1701 self.mocoeffs.append(mocoeffs) 1702 self.moenergies.append(moenergies) 1703 1704 if line[26:49] == "a o f o r c e - program": 1705 self.vibirs = [] 1706 self.vibfreqs = [] 1707 self.vibsyms = [] 1708 self.vibdisps = [] 1709 1710# while line[3:31] != "**** force : all done ****": 1711 1712 if line[12:26] == "ATOMIC WEIGHTS": 1713#begin parsing atomic weights 1714 self.vibmasses=[] 1715 line=inputfile.next() # lines ======= 1716 line=inputfile.next() # notes 1717 line=inputfile.next() # start reading 1718 temp=line.split() 1719 while(len(temp) > 0): 1720 self.vibmasses.append(float(temp[2])) 1721 line=inputfile.next() 1722 temp=line.split() 1723 1724 if line[5:14] == "frequency": 1725 if not hasattr(self,"vibfreqs"): 1726 self.vibfreqs = [] 1727 self.vibfreqs = [] 1728 self.vibsyms = [] 1729 self.vibdisps = [] 1730 self.vibirs = [] 1731 1732 temp=line.replace("i","-").split() 1733 1734 freqs = [utils.float(f) for f in temp[1:]] 1735 self.vibfreqs.extend(freqs) 1736 1737 line=inputfile.next() 1738 line=inputfile.next() 1739 1740 syms=line.split() 1741 self.vibsyms.extend(syms[1:]) 1742 1743 line=inputfile.next() 1744 line=inputfile.next() 1745 line=inputfile.next() 1746 line=inputfile.next() 1747 1748 temp=line.split() 1749 irs = [utils.float(f) for f in temp[2:]] 1750 self.vibirs.extend(irs) 1751 1752 line=inputfile.next() 1753 line=inputfile.next() 1754 line=inputfile.next() 1755 line=inputfile.next() 1756 1757 x=[] 1758 y=[] 1759 z=[] 1760 1761 line=inputfile.next() 1762 while len(line) > 1: 1763 temp=line.split() 1764 x.append(map(float, temp[3:])) 1765 1766 line=inputfile.next() 1767 temp=line.split() 1768 y.append(map(float, temp[1:])) 1769 1770 line=inputfile.next() 1771 temp=line.split() 1772 z.append(map(float, temp[1:])) 1773 line=inputfile.next() 1774 1775# build xyz vectors for each mode 1776 1777 for i in range(0, len(x[0]), 1): 1778 disp=[] 1779 for j in range(0, len(x), 1): 1780 disp.append( [x[j][i], y[j][i], z[j][i]]) 1781 self.vibdisps.append(disp) 1782 1783# line=inputfile.next() 1784 1785 def after_parsing(self): 1786 1787 # delete all frequencies that correspond to translations or rotations 1788 1789 if hasattr(self,"vibfreqs"): 1790 i = 0 1791 while i < len(self.vibfreqs): 1792 if self.vibfreqs[i]==0.0: 1793 del self.vibfreqs[i] 1794 del self.vibdisps[i] 1795 del self.vibirs[i] 1796 del self.vibsyms[i] 1797 i -= 1 1798 i += 1 1799 1800