## -*- coding: utf-8 -*- # # Copyright (c) 2018, the cclib development team # # This file is part of cclib (http://cclib.github.io) and is distributed under # the terms of the BSD 3-Clause License. """Parser for ADF output files""" from __future__ import print_function import itertools import re import numpy from cclib.parser import logfileparser from cclib.parser import utils class ADF(logfileparser.Logfile): """An ADF log file""" def __init__(self, *args, **kwargs): # Call the __init__ method of the superclass super(ADF, self).__init__(logname="ADF", *args, **kwargs) def __str__(self): """Return a string representation of the object.""" return "ADF log file %s" % (self.filename) def __repr__(self): """Return a representation of the object.""" return 'ADF("%s")' % (self.filename) def normalisesym(self, label): """Use standard symmetry labels instead of ADF labels. To normalise: (1) any periods are removed (except in the case of greek letters) (2) XXX is replaced by X, and a " added. (3) XX is replaced by X, and a ' added. (4) The greek letters Sigma, Pi, Delta and Phi are replaced by their lowercase equivalent. """ greeks = ['Sigma', 'Pi', 'Delta', 'Phi'] for greek in greeks: if label.startswith(greek): return label.lower() ans = label.replace(".", "") if ans[1:3] == "''": temp = ans[0] + '"' ans = temp l = len(ans) if l > 1 and ans[0] == ans[1]: # Python only tests the second condition if the first is true if l > 2 and ans[1] == ans[2]: ans = ans.replace(ans[0]*3, ans[0]) + '"' else: ans = ans.replace(ans[0]*2, ans[0]) + "'" return ans def normalisedegenerates(self, label, num, ndict=None): """Generate a string used for matching degenerate orbital labels To normalise: (1) if label is E or T, return label:num (2) if label is P or D, look up in dict, and return answer """ if not ndict: ndict = { 'P': {0: "P:x", 1: "P:y", 2: "P:z"}, 'D': {0: "D:z2", 1: "D:x2-y2", 2: "D:xy", 3: "D:xz", 4: "D:yz"} } if label in ndict: if num in ndict[label]: return ndict[label][num] else: return "%s:%i" % (label, num+1) else: return "%s:%i" % (label, num+1) def before_parsing(self): # Used to avoid extracting the final geometry twice in a GeoOpt self.NOTFOUND, self.GETLAST, self.NOMORE = list(range(3)) self.finalgeometry = self.NOTFOUND # Used for calculating the scftarget (variables names taken from the ADF manual) self.accint = self.SCFconv = self.sconv2 = None # keep track of nosym and unrestricted case to parse Energies since it doens't have an all Irreps section self.nosymflag = False self.unrestrictedflag = False SCFCNV, SCFCNV2 = list(range(2)) # used to index self.scftargets[] maxelem, norm = list(range(2)) # used to index scf.values def extract(self, inputfile, line): """Extract information from the file object inputfile.""" # If a file contains multiple calculations, currently we want to print a warning # and skip to the end of the file, since cclib parses only the main system, which # is usually the largest. Here we test this by checking if scftargets has already # been parsed when another INPUT FILE segment is found, although this might # not always be the best indicator. if line.strip() == "(INPUT FILE)" and hasattr(self, "scftargets"): self.logger.warning("Skipping remaining calculations") inputfile.seek(0, 2) return # We also want to check to make sure we aren't parsing "Create" jobs, # which normally come before the calculation we actually want to parse. if line.strip() == "(INPUT FILE)": while True: self.updateprogress(inputfile, "Unsupported Information", self.fupdate) line = next(inputfile) if line.strip() == "(INPUT FILE)" else None if line and not line[:6] in ("Create", "create"): break line = next(inputfile) version_searchstr = "Amsterdam Density Functional (ADF)" if version_searchstr in line: startidx = line.index(version_searchstr) + len(version_searchstr) trimmed_line = line[startidx:].strip()[:-1] # The package version is normally a year with revision # number (such as 2013.01), but it may also be a random # string (such as a version control branch name). match = re.search(r"([\d\.]{4,7})", trimmed_line) if match: package_version = match.groups()[0] self.metadata["package_version"] = package_version else: # This isn't as well-defined, but the field shouldn't # be left empty. self.metadata["package_version"] = trimmed_line.strip() # In ADF 2014.01, there are (INPUT FILE) messages, so we need to use just # the lines that start with 'Create' and run until the title or something # else we are sure is is the calculation proper. It would be good to combine # this with the previous block, if possible. if line[:6] == "Create": while line[:5] != "title" and "NO TITLE" not in line: line = inputfile.next() if line[1:10] == "Symmetry:": info = line.split() if info[1] == "NOSYM": self.nosymflag = True # Use this to read the subspecies of irreducible representations. # It will be a list, with each element representing one irrep. if line.strip() == "Irreducible Representations, including subspecies": self.skip_line(inputfile, 'dashes') self.irreps = [] line = next(inputfile) while line.strip() != "": self.irreps.append(line.split()) line = next(inputfile) if line[4:13] == 'Molecule:': info = line.split() if info[1] == 'UNrestricted': self.unrestrictedflag = True if line[1:6] == "ATOMS": # Find the number of atoms and their atomic numbers # Also extract the starting coordinates (for a GeoOpt anyway) # and the atommasses (previously called vibmasses) self.updateprogress(inputfile, "Attributes", self.cupdate) self.atomcoords = [] self.skip_lines(inputfile, ['header1', 'header2', 'header3']) atomnos = [] atommasses = [] atomcoords = [] coreelectrons = [] line = next(inputfile) while len(line) > 2: # ensure that we are reading no blank lines info = line.split() element = info[1].split('.')[0] atomnos.append(self.table.number[element]) atomcoords.append(list(map(float, info[2:5]))) coreelectrons.append(int(float(info[5]) - float(info[6]))) atommasses.append(float(info[7])) line = next(inputfile) self.atomcoords.append(atomcoords) self.set_attribute('natom', len(atomnos)) self.set_attribute('atomnos', atomnos) self.set_attribute('atommasses', atommasses) self.set_attribute('coreelectrons', coreelectrons) if line[1:10] == "FRAGMENTS": header = next(inputfile) self.frags = [] self.fragnames = [] line = next(inputfile) while len(line) > 2: # ensure that we are reading no blank lines info = line.split() if len(info) == 7: # fragment name is listed here self.fragnames.append("%s_%s" % (info[1], info[0])) self.frags.append([]) self.frags[-1].append(int(info[2]) - 1) elif len(info) == 5: # add atoms into last fragment self.frags[-1].append(int(info[0]) - 1) line = next(inputfile) # Extract charge if line[1:11] == "Net Charge": charge = int(line.split()[2]) self.set_attribute('charge', charge) line = next(inputfile) if len(line.strip()): # Spin polar: 1 (Spin_A minus Spin_B electrons) # (Not sure about this for higher multiplicities) mult = int(line.split()[2]) + 1 else: mult = 1 self.set_attribute('mult', mult) if line[1:22] == "S C F U P D A T E S": # find targets for SCF convergence if not hasattr(self, "scftargets"): self.scftargets = [] self.skip_lines(inputfile, ['e', 'b', 'numbers']) line = next(inputfile) self.SCFconv = float(line.split()[-1]) line = next(inputfile) self.sconv2 = float(line.split()[-1]) # In ADF 2013, the default numerical integration method is fuzzy cells, # although it used to be Voronoi polyhedra. Both methods apparently set # the accint parameter, although the latter does so indirectly, based on # a 'grid quality' setting. This is translated into accint using a # dictionary with values taken from the documentation. if "Numerical Integration : Voronoi Polyhedra (Te Velde)" in line: self.integration_method = "voronoi_polyhedra" if line[1:27] == 'General Accuracy Parameter': # Need to know the accuracy of the integration grid to # calculate the scftarget...note that it changes with time self.accint = float(line.split()[-1]) if "Numerical Integration : Fuzzy Cells (Becke)" in line: self.integration_method = 'fuzzy_cells' if line[1:19] == "Becke grid quality": self.grid_quality = line.split()[-1] quality2accint = { 'BASIC': 2.0, 'NORMAL': 4.0, 'GOOD': 6.0, 'VERYGOOD': 8.0, 'EXCELLENT': 10.0, } self.accint = quality2accint[self.grid_quality] # Half of the atomic orbital overlap matrix is printed since it is symmetric, # but this requires "PRINT Smat" to be in the input. There are extra blank lines # at the end of the block, which are used to terminate the parsing. # # ====== smat # # column 1 2 3 4 # row # 1 1.00000000000000E+00 # 2 2.43370854175315E-01 1.00000000000000E+00 # 3 0.00000000000000E+00 0.00000000000000E+00 1.00000000000000E+00 # ... # if "====== smat" in line: # Initialize the matrix with Nones so we can easily check all has been parsed. overlaps = [[None] * self.nbasis for i in range(self.nbasis)] self.skip_line(inputfile, 'blank') line = inputfile.next() while line.strip(): colline = line assert colline.split()[0] == "column" columns = [int(i) for i in colline.split()[1:]] rowline = inputfile.next() assert rowline.strip() == "row" line = inputfile.next() while line.strip(): i = int(line.split()[0]) vals = [float(col) for col in line.split()[1:]] for j, o in enumerate(vals): k = columns[j] overlaps[k-1][i-1] = o overlaps[i-1][k-1] = o line = inputfile.next() line = inputfile.next() # Now all values should be parsed, and so no Nones remaining. assert all([all([x is not None for x in ao]) for ao in overlaps]) self.set_attribute('aooverlaps', overlaps) if line[1:11] == "CYCLE 1": self.updateprogress(inputfile, "QM convergence", self.fupdate) newlist = [] line = next(inputfile) if not hasattr(self, "geovalues"): # This is the first SCF cycle self.scftargets.append([self.sconv2*10, self.sconv2]) elif self.finalgeometry in [self.GETLAST, self.NOMORE]: # This is the final SCF cycle self.scftargets.append([self.SCFconv*10, self.SCFconv]) else: # This is an intermediate SCF cycle in a geometry optimization, # in which case the SCF convergence target needs to be derived # from the accint parameter. For Voronoi polyhedra integration, # accint is printed and parsed. For fuzzy cells, it can be inferred # from the grid quality setting, as is done somewhere above. if self.accint: oldscftst = self.scftargets[-1][1] grdmax = self.geovalues[-1][1] scftst = max(self.SCFconv, min(oldscftst, grdmax/30, 10**(-self.accint))) self.scftargets.append([scftst*10, scftst]) while line.find("SCF CONVERGED") == -1 and line.find("SCF not fully converged, result acceptable") == -1 and line.find("SCF NOT CONVERGED") == -1: if line[4:12] == "SCF test": if not hasattr(self, "scfvalues"): self.scfvalues = [] info = line.split() newlist.append([float(info[4]), abs(float(info[6]))]) try: line = next(inputfile) except StopIteration: # EOF reached? self.logger.warning("SCF did not converge, so attributes may be missing") break if line.find("SCF not fully converged, result acceptable") > 0: self.logger.warning("SCF not fully converged, results acceptable") if line.find("SCF NOT CONVERGED") > 0: self.logger.warning("SCF did not converge! moenergies and mocoeffs are unreliable") if hasattr(self, "scfvalues"): self.scfvalues.append(newlist) # Parse SCF energy for SP calcs from bonding energy decomposition section. # It seems ADF does not print it earlier for SP calculations. # Geometry optimization runs also print this, and we want to parse it # for them, too, even if it repeats the last "Geometry Convergence Tests" # section (but it's usually a bit different). if line[:21] == "Total Bonding Energy:": if not hasattr(self, "scfenergies"): self.scfenergies = [] energy = utils.convertor(float(line.split()[3]), "hartree", "eV") self.scfenergies.append(energy) if line[51:65] == "Final Geometry": self.finalgeometry = self.GETLAST # Get the coordinates from each step of the GeoOpt. if line[1:24] == "Coordinates (Cartesian)" and self.finalgeometry in [self.NOTFOUND, self.GETLAST]: self.skip_lines(inputfile, ['e', 'b', 'title', 'title', 'd']) atomcoords = [] line = next(inputfile) while list(set(line.strip())) != ['-']: atomcoords.append(list(map(float, line.split()[5:8]))) line = next(inputfile) if not hasattr(self, "atomcoords"): self.atomcoords = [] self.atomcoords.append(atomcoords) # Don't get any more coordinates in this case. # KML: I think we could combine this with optdone (see below). if self.finalgeometry == self.GETLAST: self.finalgeometry = self.NOMORE # There have been some changes in the format of the geometry convergence information, # and this is how it is printed in older versions (2007.01 unit tests). # # ========================== # Geometry Convergence Tests # ========================== # # Energy old : -5.14170647 # new : -5.15951374 # # Convergence tests: # (Energies in hartree, Gradients in hartree/angstr or radian, Lengths in angstrom, Angles in degrees) # # Item Value Criterion Conv. Ratio # ------------------------------------------------------------------------- # change in energy -0.01780727 0.00100000 NO 0.00346330 # gradient max 0.03219530 0.01000000 NO 0.30402650 # gradient rms 0.00858685 0.00666667 NO 0.27221261 # cart. step max 0.07674971 0.01000000 NO 0.75559435 # cart. step rms 0.02132310 0.00666667 NO 0.55335378 # if line[1:27] == 'Geometry Convergence Tests': if not hasattr(self, "geotargets"): self.geovalues = [] self.geotargets = numpy.array([0.0, 0.0, 0.0, 0.0, 0.0], "d") if not hasattr(self, "scfenergies"): self.scfenergies = [] self.skip_lines(inputfile, ['e', 'b']) energies_old = next(inputfile) energies_new = next(inputfile) self.scfenergies.append(utils.convertor(float(energies_new.split()[-1]), "hartree", "eV")) self.skip_lines(inputfile, ['b', 'convergence', 'units', 'b', 'header', 'd']) values = [] for i in range(5): temp = next(inputfile).split() self.geotargets[i] = float(temp[-3]) values.append(float(temp[-4])) self.geovalues.append(values) # This is to make geometry optimization always have the optdone attribute, # even if it is to be empty for unconverged runs. if not hasattr(self, 'optdone'): self.optdone = [] # After the test, there is a message if the search is converged: # # *************************************************************************************************** # Geometry CONVERGED # *************************************************************************************************** # if line.strip() == "Geometry CONVERGED": self.skip_line(inputfile, 'stars') self.optdone.append(len(self.geovalues) - 1) # Here is the corresponding geometry convergence info from the 2013.01 unit test. # Note that the step number is given, which it will be prudent to use in an assertion. # #---------------------------------------------------------------------- #Geometry Convergence after Step 3 (Hartree/Angstrom,Angstrom) #---------------------------------------------------------------------- #current energy -5.16274478 Hartree #energy change -0.00237544 0.00100000 F #constrained gradient max 0.00884999 0.00100000 F #constrained gradient rms 0.00249569 0.00066667 F #gradient max 0.00884999 #gradient rms 0.00249569 #cart. step max 0.03331296 0.01000000 F #cart. step rms 0.00844037 0.00666667 F if line[:31] == "Geometry Convergence after Step": stepno = int(line.split()[4]) # This is to make geometry optimization always have the optdone attribute, # even if it is to be empty for unconverged runs. if not hasattr(self, 'optdone'): self.optdone = [] # The convergence message is inline in this block, not later as it was before. if "** CONVERGED **" in line: if not hasattr(self, 'optdone'): self.optdone = [] self.optdone.append(len(self.geovalues) - 1) self.skip_line(inputfile, 'dashes') current_energy = next(inputfile) energy_change = next(inputfile) constrained_gradient_max = next(inputfile) constrained_gradient_rms = next(inputfile) gradient_max = next(inputfile) gradient_rms = next(inputfile) cart_step_max = next(inputfile) cart_step_rms = next(inputfile) if not hasattr(self, "scfenergies"): self.scfenergies = [] energy = utils.convertor(float(current_energy.split()[-2]), "hartree", "eV") self.scfenergies.append(energy) if not hasattr(self, "geotargets"): self.geotargets = numpy.array([0.0, 0.0, 0.0, 0.0, 0.0], "d") self.geotargets[0] = float(energy_change.split()[-2]) self.geotargets[1] = float(constrained_gradient_max.split()[-2]) self.geotargets[2] = float(constrained_gradient_rms.split()[-2]) self.geotargets[3] = float(cart_step_max.split()[-2]) self.geotargets[4] = float(cart_step_rms.split()[-2]) if not hasattr(self, "geovalues"): self.geovalues = [] self.geovalues.append([]) self.geovalues[-1].append(float(energy_change.split()[-3])) self.geovalues[-1].append(float(constrained_gradient_max.split()[-3])) self.geovalues[-1].append(float(constrained_gradient_rms.split()[-3])) self.geovalues[-1].append(float(cart_step_max.split()[-3])) self.geovalues[-1].append(float(cart_step_rms.split()[-3])) if line.find('Orbital Energies, per Irrep and Spin') > 0 and not hasattr(self, "mosyms") and self.nosymflag and not self.unrestrictedflag: #Extracting orbital symmetries and energies, homos for nosym case #Should only be for restricted case because there is a better text block for unrestricted and nosym self.mosyms = [[]] self.moenergies = [[]] self.skip_lines(inputfile, ['e', 'header', 'd', 'label']) line = next(inputfile) info = line.split() if not info[0] == '1': self.logger.warning("MO info up to #%s is missing" % info[0]) #handle case where MO information up to a certain orbital are missing while int(info[0]) - 1 != len(self.moenergies[0]): self.moenergies[0].append(99999) self.mosyms[0].append('A') homoA = None while len(line) > 10: info = line.split() self.mosyms[0].append('A') self.moenergies[0].append(utils.convertor(float(info[2]), 'hartree', 'eV')) if info[1] == '0.000' and not hasattr(self, 'homos'): self.set_attribute('homos', [len(self.moenergies[0]) - 2]) line = next(inputfile) self.moenergies = [numpy.array(self.moenergies[0], "d")] if line[1:29] == 'Orbital Energies, both Spins' and not hasattr(self, "mosyms") and self.nosymflag and self.unrestrictedflag: #Extracting orbital symmetries and energies, homos for nosym case #should only be here if unrestricted and nosym self.mosyms = [[], []] moenergies = [[], []] self.skip_lines(inputfile, ['d', 'b', 'header', 'd']) homoa = 0 homob = None line = next(inputfile) while len(line) > 5: info = line.split() if info[2] == 'A': self.mosyms[0].append('A') moenergies[0].append(utils.convertor(float(info[4]), 'hartree', 'eV')) if info[3] != '0.00': homoa = len(moenergies[0]) - 1 elif info[2] == 'B': self.mosyms[1].append('A') moenergies[1].append(utils.convertor(float(info[4]), 'hartree', 'eV')) if info[3] != '0.00': homob = len(moenergies[1]) - 1 else: print(("Error reading line: %s" % line)) line = next(inputfile) self.moenergies = [numpy.array(x, "d") for x in moenergies] self.set_attribute('homos', [homoa, homob]) # Extracting orbital symmetries and energies, homos. if line[1:29] == 'Orbital Energies, all Irreps' and not hasattr(self, "mosyms"): self.symlist = {} self.mosyms = [[]] self.moenergies = [[]] self.skip_lines(inputfile, ['e', 'b', 'header', 'd']) homoa = None homob = None #multiple = {'E':2, 'T':3, 'P':3, 'D':5} # The above is set if there are no special irreps names = [irrep[0].split(':')[0] for irrep in self.irreps] counts = [len(irrep) for irrep in self.irreps] multiple = dict(list(zip(names, counts))) irrepspecies = {} for n in range(len(names)): indices = list(range(counts[n])) subspecies = self.irreps[n] irrepspecies[names[n]] = dict(list(zip(indices, subspecies))) line = next(inputfile) while line.strip(): info = line.split() if len(info) == 5: # this is restricted #count = multiple.get(info[0][0],1) count = multiple.get(info[0], 1) for repeat in range(count): # i.e. add E's twice, T's thrice self.mosyms[0].append(self.normalisesym(info[0])) self.moenergies[0].append(utils.convertor(float(info[3]), 'hartree', 'eV')) sym = info[0] if count > 1: # add additional sym label sym = self.normalisedegenerates(info[0], repeat, ndict=irrepspecies) try: self.symlist[sym][0].append(len(self.moenergies[0])-1) except KeyError: self.symlist[sym] = [[]] self.symlist[sym][0].append(len(self.moenergies[0])-1) if info[2] == '0.00' and not hasattr(self, 'homos'): self.homos = [len(self.moenergies[0]) - (count + 1)] # count, because need to handle degenerate cases line = next(inputfile) elif len(info) == 6: # this is unrestricted if len(self.moenergies) < 2: # if we don't have space, create it self.moenergies.append([]) self.mosyms.append([]) # count = multiple.get(info[0][0], 1) count = multiple.get(info[0], 1) if info[2] == 'A': for repeat in range(count): # i.e. add E's twice, T's thrice self.mosyms[0].append(self.normalisesym(info[0])) self.moenergies[0].append(utils.convertor(float(info[4]), 'hartree', 'eV')) sym = info[0] if count > 1: # add additional sym label sym = self.normalisedegenerates(info[0], repeat) try: self.symlist[sym][0].append(len(self.moenergies[0])-1) except KeyError: self.symlist[sym] = [[], []] self.symlist[sym][0].append(len(self.moenergies[0])-1) if info[3] == '0.00' and homoa is None: homoa = len(self.moenergies[0]) - (count + 1) # count because degenerate cases need to be handled if info[2] == 'B': for repeat in range(count): # i.e. add E's twice, T's thrice self.mosyms[1].append(self.normalisesym(info[0])) self.moenergies[1].append(utils.convertor(float(info[4]), 'hartree', 'eV')) sym = info[0] if count > 1: # add additional sym label sym = self.normalisedegenerates(info[0], repeat) try: self.symlist[sym][1].append(len(self.moenergies[1])-1) except KeyError: self.symlist[sym] = [[], []] self.symlist[sym][1].append(len(self.moenergies[1])-1) if info[3] == '0.00' and homob is None: homob = len(self.moenergies[1]) - (count + 1) line = next(inputfile) else: # different number of lines print(("Error", info)) if len(info) == 6: # still unrestricted, despite being out of loop self.set_attribute('homos', [homoa, homob]) self.moenergies = [numpy.array(x, "d") for x in self.moenergies] # Section on extracting vibdisps # Also contains vibfreqs, but these are extracted in the # following section (see below) if line[1:28] == "Vibrations and Normal Modes": self.vibdisps = [] self.skip_lines(inputfile, ['e', 'b', 'header', 'header', 'b', 'b']) freqs = next(inputfile) while freqs.strip() != "": minus = next(inputfile) p = [[], [], []] for i in range(len(self.atomnos)): broken = list(map(float, next(inputfile).split()[1:])) for j in range(0, len(broken), 3): p[j//3].append(broken[j:j+3]) self.vibdisps.extend(p[:(len(broken)//3)]) self.skip_lines(inputfile, ['b', 'b']) freqs = next(inputfile) self.vibdisps = numpy.array(self.vibdisps, "d") if line[1:24] == "List of All Frequencies": # Start of the IR/Raman frequency section self.updateprogress(inputfile, "Frequency information", self.fupdate) # self.vibsyms = [] # Need to look into this a bit more self.vibirs = [] self.vibfreqs = [] for i in range(8): line = next(inputfile) line = next(inputfile).strip() while line: temp = line.split() self.vibfreqs.append(float(temp[0])) self.vibirs.append(float(temp[2])) # or is it temp[1]? line = next(inputfile).strip() self.vibfreqs = numpy.array(self.vibfreqs, "d") self.vibirs = numpy.array(self.vibirs, "d") if hasattr(self, "vibramans"): self.vibramans = numpy.array(self.vibramans, "d") #******************************************************************************************************************8 #delete this after new implementation using smat, eigvec print,eprint? # Extract the number of basis sets if line[1:49] == "Total nr. of (C)SFOs (summation over all irreps)": nbasis = int(line.split(":")[1].split()[0]) self.set_attribute('nbasis', nbasis) # now that we're here, let's extract aonames self.fonames = [] self.start_indeces = {} self.atombasis = [[] for frag in self.frags] # parse atombasis in the case of trivial SFOs self.skip_line(inputfile, 'blank') note = next(inputfile) symoffset = 0 self.skip_line(inputfile, 'blank') line = next(inputfile) if len(line) > 2: # fix for ADF2006.01 as it has another note self.skip_line(inputfile, 'blank') line = next(inputfile) self.skip_line(inputfile, 'blank') self.nosymreps = [] while len(self.fonames) < self.nbasis: symline = next(inputfile) sym = symline.split()[1] line = next(inputfile) num = int(line.split(':')[1].split()[0]) self.nosymreps.append(num) #read until line "--------..." is found while line.find('-----') < 0: line = next(inputfile) line = next(inputfile) # the start of the first SFO while len(self.fonames) < symoffset + num: info = line.split() #index0 index1 occ2 energy3/4 fragname5 coeff6 orbnum7 orbname8 fragname9 if not sym in list(self.start_indeces.keys()): #have we already set the start index for this symmetry? self.start_indeces[sym] = int(info[1]) orbname = info[8] orbital = info[7] + orbname.replace(":", "") fragname = info[5] frag = fragname + info[9] coeff = float(info[6]) # parse atombasis only in the case that all coefficients are 1 # and delete it otherwise if hasattr(self, 'atombasis'): if coeff == 1.: ibas = int(info[0]) - 1 ifrag = int(info[9]) - 1 iat = self.frags[ifrag][0] self.atombasis[iat].append(ibas) else: del self.atombasis line = next(inputfile) while line.strip() and not line[:7].strip(): # while it's the same SFO # i.e. while not completely blank, but blank at the start info = line[43:].split() if len(info) > 0: # len(info)==0 for the second line of dvb_ir.adfout frag += "+" + fragname + info[-1] coeff = float(info[-4]) if coeff < 0: orbital += '-' + info[-3] + info[-2].replace(":", "") else: orbital += '+' + info[-3] + info[-2].replace(":", "") line = next(inputfile) # At this point, we are either at the start of the next SFO or at # a blank line...the end self.fonames.append("%s_%s" % (frag, orbital)) symoffset += num # blankline blankline next(inputfile) next(inputfile) if line[1:32] == "S F O P O P U L A T I O N S ,": #Extract overlap matrix # self.fooverlaps = numpy.zeros((self.nbasis, self.nbasis), "d") symoffset = 0 for nosymrep in self.nosymreps: line = next(inputfile) while line.find('===') < 10: # look for the symmetry labels line = next(inputfile) self.skip_lines(inputfile, ['b', 'b']) text = next(inputfile) if text[13:20] != "Overlap": # verify this has overlap info break self.skip_lines(inputfile, ['b', 'col', 'row']) if not hasattr(self, "fooverlaps"): # make sure there is a matrix to store this self.fooverlaps = numpy.zeros((self.nbasis, self.nbasis), "d") base = 0 while base < nosymrep: # have we read all the columns? for i in range(nosymrep - base): self.updateprogress(inputfile, "Overlap", self.fupdate) line = next(inputfile) parts = line.split()[1:] for j in range(len(parts)): k = float(parts[j]) self.fooverlaps[base + symoffset + j, base + symoffset + i] = k self.fooverlaps[base + symoffset + i, base + symoffset + j] = k #blank, blank, column for i in range(3): next(inputfile) base += 4 symoffset += nosymrep base = 0 # The commented code below makes the atombasis attribute based on the BAS function in ADF, # but this is probably not so useful, since SFOs are used to build MOs in ADF. # if line[1:54] == "BAS: List of all Elementary Cartesian Basis Functions": # # self.atombasis = [] # # # There will be some text, followed by a line: # # (power of) X Y Z R Alpha on Atom # while not line[1:11] == "(power of)": # line = inputfile.next() # dashes = inputfile.next() # blank = inputfile.next() # line = inputfile.next() # # There will be two blank lines when there are no more atom types. # while line.strip() != "": # atoms = [int(i)-1 for i in line.split()[1:]] # for n in range(len(atoms)): # self.atombasis.append([]) # dashes = inputfile.next() # line = inputfile.next() # while line.strip() != "": # indices = [int(i)-1 for i in line.split()[5:]] # for i in range(len(indices)): # self.atombasis[atoms[i]].append(indices[i]) # line = inputfile.next() # line = inputfile.next() if line[48:67] == "SFO MO coefficients": self.mocoeffs = [numpy.zeros((self.nbasis, self.nbasis), "d")] spin = 0 symoffset = 0 lastrow = 0 # Section ends with "1" at beggining of a line. while line[0] != "1": line = next(inputfile) # If spin is specified, then there will be two coefficient matrices. if line.strip() == "***** SPIN 1 *****": self.mocoeffs = [numpy.zeros((self.nbasis, self.nbasis), "d"), numpy.zeros((self.nbasis, self.nbasis), "d")] # Bump up the spin. if line.strip() == "***** SPIN 2 *****": spin = 1 symoffset = 0 lastrow = 0 # Next symmetry. if line.strip()[:4] == "=== ": sym = line.split()[1] if self.nosymflag: aolist = list(range(self.nbasis)) else: aolist = self.symlist[sym][spin] # Add to the symmetry offset of AO ordering. symoffset += lastrow # Blocks with coefficient always start with "MOs :". if line[1:6] == "MOs :": # Next line has the MO index contributed to. monumbers = [int(n) for n in line[6:].split()] self.skip_lines(inputfile, ['occup', 'label']) # The table can end with a blank line or "1". row = 0 line = next(inputfile) while not line.strip() in ["", "1"]: info = line.split() if int(info[0]) < self.start_indeces[sym]: #check to make sure we aren't parsing CFs line = next(inputfile) continue self.updateprogress(inputfile, "Coefficients", self.fupdate) row += 1 coeffs = [float(x) for x in info[1:]] moindices = [aolist[n-1] for n in monumbers] # The AO index is 1 less than the row. aoindex = symoffset + row - 1 for i in range(len(monumbers)): self.mocoeffs[spin][moindices[i], aoindex] = coeffs[i] line = next(inputfile) lastrow = row # ************************************************************************** # * * # * Final excitation energies from Davidson algorithm * # * * # ************************************************************************** # # Number of loops in Davidson routine = 20 # Number of matrix-vector multiplications = 24 # Type of excitations = SINGLET-SINGLET # # Symmetry B.u # # ... several blocks ... # # Normal termination of EXCITATION program part if line[4:53] == "Final excitation energies from Davidson algorithm": while line[1:9] != "Symmetry" and "Normal termination" not in line: line = next(inputfile) symm = self.normalisesym(line.split()[1]) # Excitation energies E in a.u. and eV, dE wrt prev. cycle, # oscillator strengths f in a.u. # # no. E/a.u. E/eV f dE/a.u. # ----------------------------------------------------- # 1 0.17084 4.6488 0.16526E-01 0.28E-08 # ... while line.split() != ['no.', 'E/a.u.', 'E/eV', 'f', 'dE/a.u.'] and "Normal termination" not in line: line = next(inputfile) self.skip_line(inputfile, 'dashes') etenergies = [] etoscs = [] etsyms = [] line = next(inputfile) while len(line) > 2: info = line.split() etenergies.append(utils.convertor(float(info[2]), "eV", "wavenumber")) etoscs.append(float(info[3])) etsyms.append(symm) line = next(inputfile) # There is another section before this, with transition dipole moments, # but this should just skip past it. while line[1:53] != "Major MO -> MO transitions for the above excitations": line = next(inputfile) # Note that here, and later, the number of blank lines can vary between # version of ADF (extra lines are seen in 2013.01 unit tests, for example). self.skip_line(inputfile, 'blank') excitation_occupied = next(inputfile) header = next(inputfile) while not header.strip(): header = next(inputfile) header2 = next(inputfile) x_y_z = next(inputfile) line = next(inputfile) while not line.strip(): line = next(inputfile) # Before we start handeling transitions, we need to create mosyms # with indices; only restricted calcs are possible in ADF. counts = {} syms = [] for mosym in self.mosyms[0]: if list(counts.keys()).count(mosym) == 0: counts[mosym] = 1 else: counts[mosym] += 1 syms.append(str(counts[mosym]) + mosym) etsecs = [] printed_warning = False for i in range(len(etenergies)): etsec = [] info = line.split() while len(info) > 0: match = re.search('[^0-9]', info[1]) index1 = int(info[1][:match.start(0)]) text = info[1][match.start(0):] symtext = text[0].upper() + text[1:] sym1 = str(index1) + self.normalisesym(symtext) match = re.search('[^0-9]', info[3]) index2 = int(info[3][:match.start(0)]) text = info[3][match.start(0):] symtext = text[0].upper() + text[1:] sym2 = str(index2) + self.normalisesym(symtext) try: index1 = syms.index(sym1) except ValueError: if not printed_warning: self.logger.warning("Etsecs are not accurate!") printed_warning = True try: index2 = syms.index(sym2) except ValueError: if not printed_warning: self.logger.warning("Etsecs are not accurate!") printed_warning = True etsec.append([(index1, 0), (index2, 0), float(info[4])]) line = next(inputfile) info = line.split() etsecs.append(etsec) # Again, the number of blank lines between transition can vary. line = next(inputfile) while not line.strip(): line = next(inputfile) if not hasattr(self, "etenergies"): self.etenergies = etenergies else: self.etenergies += etenergies if not hasattr(self, "etoscs"): self.etoscs = etoscs else: self.etoscs += etoscs if not hasattr(self, "etsyms"): self.etsyms = etsyms else: self.etsyms += etsyms if not hasattr(self, "etsecs"): self.etsecs = etsecs else: self.etsecs += etsecs if "M U L L I K E N P O P U L A T I O N S" in line: if not hasattr(self, "atomcharges"): self.atomcharges = {} while line[1:5] != "Atom": line = next(inputfile) self.skip_line(inputfile, 'dashes') mulliken = [] line = next(inputfile) while line.strip(): mulliken.append(float(line.split()[2])) line = next(inputfile) self.atomcharges["mulliken"] = mulliken # Dipole moment is always printed after a point calculation, # and the reference point for this is always the origin (0,0,0) # and not necessarily the center of mass, as explained on the # ADF user mailing list (see cclib/cclib#113 for details). # # ============= # Dipole Moment *** (Debye) *** # ============= # # Vector : 0.00000000 0.00000000 0.00000000 # Magnitude: 0.00000000 # if line.strip()[:13] == "Dipole Moment": self.skip_line(inputfile, 'equals') # There is not always a blank line here, for example when the dipole and quadrupole # moments are printed after the multipole derived atomic charges. Still, to the best # of my knowledge (KML) the values are still in Debye. line = next(inputfile) if not line.strip(): line = next(inputfile) assert line.split()[0] == "Vector" dipole = [float(d) for d in line.split()[-3:]] reference = [0.0, 0.0, 0.0] if not hasattr(self, 'moments'): self.moments = [reference, dipole] else: try: assert self.moments[1] == dipole except AssertionError: self.logger.warning('Overwriting previous multipole moments with new values') self.moments = [reference, dipole] # Molecular response properties. if line.strip()[1:-1].strip() == "RESPONSE program part": while line.strip() != "Normal termination of RESPONSE program part": if "THE DIPOLE-DIPOLE POLARIZABILITY TENSOR:" in line: if not hasattr(self, 'polarizabilities'): self.polarizabilities = [] polarizability = numpy.empty(shape=(3, 3)) self.skip_lines(inputfile, ['b', 'FREQUENCY', 'coordinates']) # Ordering of rows/columns is Y, Z, X. ordering = [1, 2, 0] indices = list(itertools.product(ordering, ordering)) for i in range(3): tokens = next(inputfile).split() for j in range(3): polarizability[indices[(i*3)+j]] = tokens[j] self.polarizabilities.append(polarizability) line = next(inputfile) if line[:24] == ' Buffered I/O statistics': self.metadata['success'] = True