# -*- coding: utf-8 -*- from __future__ import print_function, absolute_import import ctypes, math, os.path, sys import numpy as np # Importing regardless of relative import try: from .io import * except: from io import * # PHYSICAL CONSTANTS UNITS GAS_CONSTANT = 8.3144621 # J / K / mol PLANCK_CONSTANT = 6.62606957e-34 # J * s BOLTZMANN_CONSTANT = 1.3806488e-23 # J / K SPEED_OF_LIGHT = 2.99792458e10 # cm / s AVOGADRO_CONSTANT = 6.0221415e23 # 1 / mol AMU_to_KG = 1.66053886E-27 # UNIT CONVERSION J_TO_AU = 4.184 * 627.509541 * 1000.0 # UNIT CONVERSION # Symmetry numbers for different point groups pg_sm = {"C1": 1, "Cs": 1, "Ci": 1, "C2": 2, "C3": 3, "C4": 4, "C5": 5, "C6": 6, "C7": 7, "C8": 8, "D2": 4, "D3": 6, "D4": 8, "D5": 10, "D6": 12, "D7": 14, "D8": 16, "C2v": 2, "C3v": 3, "C4v": 4, "C5v": 5, "C6v": 6, "C7v": 7, "C8v": 8, "C2h": 2, "C3h": 3, "C4h": 4, "C5h": 5, "C6h": 6, "C7h": 7, "C8h": 8, "D2h": 4, "D3h": 6, "D4h": 8, "D5h": 10, "D6h": 12, "D7h": 14, "D8h": 16, "D2d": 4, "D3d": 6, "D4d": 8, "D5d": 10, "D6d": 12, "D7d": 14, "D8d": 16, "S4": 4, "S6": 6, "S8": 8, "T": 6, "Th": 12, "Td": 12, "O": 12, "Oh": 24, "Cinfv": 1, "Dinfh": 2, "I": 30, "Ih": 60, "Kh": 1} def sharepath(filename): """ Get absolute pathway to GoodVibes project. Used in finding location of compiled C files used in symmetry corrections. Parameter: filename (str): name of compiled C file, OS specific. Returns: str: absolute path on machine to compiled C file. """ here = os.path.dirname(os.path.abspath(__file__)) return os.path.join(here, 'share', filename) def calc_translational_energy(temperature): """ Translational energy evaluation Calculates the translational energy (J/mol) of an ideal gas. i.e. non-interacting molecules so molar energy = Na * atomic energy. This approximation applies to all energies and entropies computed within. Etrans = 3/2 RT! Parameter: temperature (float): temperature for calculations to be performed at. Returns: float: translational energy of chemical system. """ energy = 1.5 * GAS_CONSTANT * temperature return energy def calc_rotational_energy(zpe, symmno, temperature, linear): """ Rotational energy evaluation Calculates the rotational energy (J/mol) Etrans = 0 (atomic) ; RT (linear); 3/2 RT (non-linear) Parameters: zpe (float): zero point energy of chemical system. symmno (float): symmetry number, used for adding a symmetry correction. temperature (float): temperature for calculations to be performed at. linear (bool): flag for linear molecules, changes how calculation is performed. Returns: float: rotational energy of chemical system. """ if zpe == 0.0: energy = 0.0 elif linear == 1: energy = GAS_CONSTANT * temperature else: energy = 1.5 * GAS_CONSTANT * temperature return energy def calc_vibrational_energy(frequency_wn, temperature, freq_scale_factor, fract_modelsys): """ Vibrational energy evaluation. Calculates the vibrational energy contribution (J/mol). Includes ZPE (0K) and thermal contributions. Evib = R * Sum(0.5 hv/k + (hv/k)/(e^(hv/KT)-1)) Parameters: frequency_wn (list): list of frequencies parsed from file. temperature (float): temperature for calculations to be performed at. freq_scale_factor (float): frequency scaling factor based on level of theory and basis set used. fract_modelsys (list): MM frequency scale factors obtained from ONIOM calculations. Returns: float: vibrational energy of chemical system. """ if fract_modelsys is not False: freq_scale_factor = [freq_scale_factor[0] * fract_modelsys[i] + freq_scale_factor[1] * (1.0 - fract_modelsys[i]) for i in range(len(fract_modelsys))] factor = [(PLANCK_CONSTANT * frequency_wn[i] * SPEED_OF_LIGHT * freq_scale_factor[i]) / (BOLTZMANN_CONSTANT * temperature) for i in range(len(frequency_wn))] else: factor = [(PLANCK_CONSTANT * freq * SPEED_OF_LIGHT * freq_scale_factor) / (BOLTZMANN_CONSTANT * temperature) for freq in frequency_wn] # Error occurs if T is too low when performing math.exp for entry in factor: if entry > math.log(sys.float_info.max): sys.exit("\nx Warning! Temperature may be too low to calculate vibrational energy. Please adjust using the `-t` option and try again.\n") energy = [entry * GAS_CONSTANT * temperature * (0.5 + (1.0 / (math.exp(entry) - 1.0))) for entry in factor] return sum(energy) def calc_zeropoint_energy(frequency_wn, freq_scale_factor, fract_modelsys): """ Vibrational Zero point energy evaluation. Calculates the vibrational ZPE (J/mol) EZPE = Sum(0.5 hv/k) Parameters: frequency_wn (list): list of frequencies parsed from file. freq_scale_factor (float): frequency scaling factor based on level of theory and basis set used. fract_modelsys (list): MM frequency scale factors obtained from ONIOM calculations. Returns: float: zerp point energy of chemical system. """ if fract_modelsys is not False: freq_scale_factor = [freq_scale_factor[0] * fract_modelsys[i] + freq_scale_factor[1] * (1.0 - fract_modelsys[i]) for i in range(len(fract_modelsys))] factor = [(PLANCK_CONSTANT * frequency_wn[i] * SPEED_OF_LIGHT * freq_scale_factor[i]) / (BOLTZMANN_CONSTANT) for i in range(len(frequency_wn))] else: factor = [(PLANCK_CONSTANT * freq * SPEED_OF_LIGHT * freq_scale_factor) / (BOLTZMANN_CONSTANT) for freq in frequency_wn] energy = [0.5 * entry * GAS_CONSTANT for entry in factor] return sum(energy) def get_free_space(solv): """ Computed the amount of accessible free space (ml per L) in solution. Calculates the free space in a litre of bulk solvent, based on Shakhnovich and Whitesides (J. Org. Chem. 1998, 63, 3821-3830). Free space based on accessible to a solute immersed in bulk solvent, i.e. this is the volume not occupied by solvent molecules, calculated using literature values for molarity and B3LYP/6-31G* computed molecular volumes. Parameter: solv (str): solvent used in chemical calculation. Returns: float: accessible free space in solution. """ solvent_list = ["none", "H2O", "toluene", "DMF", "AcOH", "chloroform"] molarity = [1.0, 55.6, 9.4, 12.9, 17.4, 12.5] # mol/l molecular_vol = [1.0, 27.944, 149.070, 77.442, 86.10, 97.0] # Angstrom^3 nsolv = 0 for i in range(0, len(solvent_list)): if solv == solvent_list[i]: nsolv = i solv_molarity = molarity[nsolv] solv_volume = molecular_vol[nsolv] if nsolv > 0: v_free = 8 * ((1E27 / (solv_molarity * AVOGADRO_CONSTANT)) ** 0.333333 - solv_volume ** 0.333333) ** 3 freespace = v_free * solv_molarity * AVOGADRO_CONSTANT * 1E-24 else: freespace = 1000.0 return freespace def calc_translational_entropy(molecular_mass, conc, temperature, solv): """ Translational entropy evaluation. Calculates the translational entropic contribution (J/(mol*K)) of an ideal gas. Needs the molecular mass. Convert mass in amu to kg; conc in mol/l to number per m^3 Strans = R(Ln(2pimkT/h^2)^3/2(1/C)) + 1 + 3/2) Parameters: molecular_mass (float): total molecular mass of chemical system. conc (float): concentration to perform calculations at. temperature (float): temperature for calculations to be performed at. solv (str): solvent used in chemical calculation. Returns: float: translational entropy of chemical system. """ lmda = ((2.0 * math.pi * molecular_mass * AMU_to_KG * BOLTZMANN_CONSTANT * temperature) ** 0.5) / PLANCK_CONSTANT freespace = get_free_space(solv) ndens = conc * 1000 * AVOGADRO_CONSTANT / (freespace / 1000.0) entropy = GAS_CONSTANT * (2.5 + math.log(lmda ** 3 / ndens)) return entropy def calc_electronic_entropy(multiplicity): """ Electronic entropy evaluation. Calculates the electronic entropic contribution (J/(mol*K)) of the molecule Selec = R(Ln(multiplicity) Parameter: multiplicity (int): multiplicity of chemical system. Returns: float: electronic entropy of chemical system. """ entropy = GAS_CONSTANT * (math.log(multiplicity)) return entropy def calc_rotational_entropy(zpe, linear, symmno, rotemp, temperature): """ Rotational entropy evaluation. Calculates the rotational entropy (J/(mol*K)) Strans = 0 (atomic) ; R(Ln(q)+1) (linear); R(Ln(q)+3/2) (non-linear) Parameters: zpe (float): zero point energy of chemical system. linear (bool): flag for linear molecules. symmno (float): symmetry number of chemical system. rotemp (list): list of parsed rotational temperatures of chemical system. temperature (float): temperature for calculations to be performed at. Returns: float: rotational entropy of chemical system. """ if rotemp == [0.0, 0.0, 0.0] or zpe == 0.0: # Monatomic entropy = 0.0 else: if len(rotemp) == 1: # Diatomic or linear molecules linear = 1 qrot = temperature / rotemp[0] elif len(rotemp) == 2: # Possible gaussian problem with linear triatomic linear = 2 else: qrot = math.pi * temperature ** 3 / (rotemp[0] * rotemp[1] * rotemp[2]) qrot = qrot ** 0.5 if linear == 1: entropy = GAS_CONSTANT * (math.log(qrot / symmno) + 1) elif linear == 2: entropy = 0.0 else: entropy = GAS_CONSTANT * (math.log(qrot / symmno) + 1.5) return entropy def calc_rrho_entropy(frequency_wn, temperature, freq_scale_factor, fract_modelsys): """ Rigid rotor harmonic oscillator (RRHO) entropy evaluation - this is the default treatment Entropic contributions (J/(mol*K)) according to a rigid-rotor harmonic-oscillator description for a list of vibrational modes Sv = RSum(hv/(kT(e^(hv/kT)-1) - ln(1-e^(-hv/kT))) Parameters: frequency_wn (list): list of frequencies parsed from file. temperature (float): temperature for calculations to be performed at. freq_scale_factor (float): frequency scaling factor based on level of theory and basis set used. fract_modelsys (list): MM frequency scale factors obtained from ONIOM calculations. Returns: float: RRHO entropy of chemical system. """ if fract_modelsys is not False: freq_scale_factor = [freq_scale_factor[0] * fract_modelsys[i] + freq_scale_factor[1] * (1.0 - fract_modelsys[i]) for i in range(len(fract_modelsys))] factor = [(PLANCK_CONSTANT * frequency_wn[i] * SPEED_OF_LIGHT * freq_scale_factor[i]) / (BOLTZMANN_CONSTANT * temperature) for i in range(len(frequency_wn))] else: factor = [(PLANCK_CONSTANT * freq * SPEED_OF_LIGHT * freq_scale_factor) / (BOLTZMANN_CONSTANT * temperature) for freq in frequency_wn] entropy = [entry * GAS_CONSTANT / (math.exp(entry) - 1) - GAS_CONSTANT * math.log(1 - math.exp(-entry)) for entry in factor] return entropy def calc_qRRHO_energy(frequency_wn, temperature, freq_scale_factor): """ Quasi-rigid rotor harmonic oscillator energy evaluation. Head-Gordon RRHO-vibrational energy contribution (J/mol*K) of vibrational modes described by a rigid-rotor harmonic approximation. V_RRHO = 1/2(Nhv) + RT(hv/kT)e^(-hv/kT)/(1-e^(-hv/kT)) Parameters: frequency_wn (list): list of frequencies parsed from file. temperature (float): temperature for calculations to be performed at. freq_scale_factor (float): frequency scaling factor based on level of theory and basis set used. Returns: float: quasi-RRHO energy of chemical system. """ factor = [PLANCK_CONSTANT * freq * SPEED_OF_LIGHT * freq_scale_factor for freq in frequency_wn] energy = [0.5 * AVOGADRO_CONSTANT * entry + GAS_CONSTANT * temperature * entry / BOLTZMANN_CONSTANT / temperature * math.exp(-entry / BOLTZMANN_CONSTANT / temperature) / (1 - math.exp(-entry / BOLTZMANN_CONSTANT / temperature)) for entry in factor] return energy def calc_freerot_entropy(frequency_wn, temperature, freq_scale_factor, fract_modelsys, file, inertia, roconst): """ Free rotor entropy evaluation. Entropic contributions (J/(mol*K)) according to a free-rotor description for a list of vibrational modes Sr = R(1/2 + 1/2ln((8pi^3u'kT/h^2)) Parameters: frequency_wn (list): list of frequencies parsed from file. temperature (float): temperature for calculations to be performed at. freq_scale_factor (float): frequency scaling factor based on level of theory and basis set used. fract_modelsys (list): MM frequency scale factors obtained from ONIOM calculations. inertia (str): flag for choosing global average moment of inertia for all molecules or computing individually from parsed rotational constants roconst (list): list of parsed rotational constants for computing the average moment of inertia. Returns: float: free rotor entropy of chemical system. """ # This is the average moment of inertia used by Grimme if inertia == "global" or len(roconst) == 0: bav = 1.00e-44 else: av_roconst_ghz = sum(roconst)/len(roconst) #GHz av_roconst_hz = av_roconst_ghz * 1000000000 #Hz av_roconst_s = 1 / av_roconst_hz #s av_roconst = av_roconst_s * PLANCK_CONSTANT #kg m^2 bav = av_roconst if fract_modelsys is not False: freq_scale_factor = [freq_scale_factor[0] * fract_modelsys[i] + freq_scale_factor[1] * (1.0 - fract_modelsys[i]) for i in range(len(fract_modelsys))] mu = [PLANCK_CONSTANT / (8 * math.pi ** 2 * frequency_wn[i] * SPEED_OF_LIGHT * freq_scale_factor[i]) for i in range(len(frequency_wn))] else: mu = [PLANCK_CONSTANT / (8 * math.pi ** 2 * freq * SPEED_OF_LIGHT * freq_scale_factor) for freq in frequency_wn] mu_primed = [entry * bav / (entry + bav) for entry in mu] factor = [8 * math.pi ** 3 * entry * BOLTZMANN_CONSTANT * temperature / PLANCK_CONSTANT ** 2 for entry in mu_primed] entropy = [(0.5 + math.log(entry ** 0.5)) * GAS_CONSTANT for entry in factor] return entropy def calc_damp(frequency_wn, freq_cutoff): """A damping function to interpolate between RRHO and free rotor vibrational entropy values""" alpha = 4 damp = [1 / (1 + (freq_cutoff / entry) ** alpha) for entry in frequency_wn] return damp class calc_bbe: """ The function to compute the "black box" entropy and enthalpy values along with all other thermochemical quantities. Parses energy, program version, frequencies, charge, multiplicity, solvation model, computation time. Computes H, S from partition functions, applying qhasi-harmonic corrections, COSMO-RS solvation corrections, considering frequency scaling factors from detected level of theory/basis set, and optionally ONIOM frequency scaling. Attributes: xyz (getoutData object): contains Cartesian coordinates, atom connectivity. job_type (str): contains information on the type of Gaussian job such as ground or transition state optimization, frequency. roconst (list): list of parsed rotational constants from Gaussian calculations. program (str): program used in chemical computation. version_program (str): program version used in chemical computation. solvation_model (str): solvation model used in chemical computation. file (str): input chemical computation output file. charge (int): overall charge of molecule. empirical_dispersion (str): empirical dispersion model used in computation. multiplicity (int): multiplicity of molecule or chemical system. mult (int): multiplicity of molecule or chemical system. point_group (str): point group of molecule or chemical system used for symmetry corrections. sp_energy (float): single-point energy parsed from output file. sp_program (str): program used for single-point energy calculation. sp_version_program (str): version of program used for single-point energy calculation. sp_solvation_model (str): solvation model used for single-point energy calculation. sp_file (str): single-point energy calculation output file. sp_charge (int): overall charge of molecule in single-point energy calculation. sp_empirical_dispersion (str): empirical dispersion model used in single-point energy computation. sp_multiplicity (int): multiplicity of molecule or chemical system in single-point energy computation. cpu (list): days, hours, mins, secs, msecs of computation time. scf_energy (float): self-consistent field energy. frequency_wn (list): frequencies parsed from chemical computation output file. im_freq (list): imaginary frequencies parsed from chemical computation output file. inverted_freqs (list): frequencies inverted from imaginary to real numbers. zero_point_corr (float): thermal corrections for zero-point energy parsed from file. zpe (float): vibrational zero point energy computed from frequencies. enthalpy (float): enthalpy computed from partition functions. qh_enthalpy (float): enthalpy computed from partition functions, quasi-harmonic corrections applied. entropy (float): entropy of chemical system computed from partition functions. qh_entropy (float): entropy of chemical system computed from partition functions, quasi-harmonic corrections applied. gibbs_free_energy (float): Gibbs free energy of chemical system computed from enthalpy and entropy. qh_gibbs_free_energy (float): Gibbs free energy of chemical system computed from quasi-harmonic enthalpy and/or entropy. cosmo_qhg (float): quasi-harmonic Gibbs free energy with COSMO-RS correction for Gibbs free energy of solvation linear_warning (bool): flag for linear molecules, may be missing a rotational constant. """ def __init__(self, file, QS, QH, s_freq_cutoff, H_FREQ_CUTOFF, temperature, conc, freq_scale_factor, solv, spc, invert, d3_term, ssymm=False, cosmo=None, mm_freq_scale_factor=False,inertia='global',g4=False): # List of frequencies and default values im_freq_cutoff, frequency_wn, im_frequency_wn, rotemp, roconst, linear_mol, link, freqloc, linkmax, symmno, self.cpu, inverted_freqs = 0.0, [], [], [ 0.0, 0.0, 0.0], [0.0, 0.0, 0.0], 0, 0, 0, 0, 1, [0, 0, 0, 0, 0], [] linear_warning = False if mm_freq_scale_factor is False: fract_modelsys = False else: fract_modelsys = [] freq_scale_factor = [freq_scale_factor, mm_freq_scale_factor] self.xyz = getoutData(file) self.job_type = gaussian_jobtype(file) self.roconst = [] # Parse some useful information from the file self.sp_energy, self.program, self.version_program, self.solvation_model, self.file, self.charge, self.empirical_dispersion, self.multiplicity = parse_data( file) with open(file) as f: g_output = f.readlines() self.cosmo_qhg = 0.0 # Read any single point energies if requested if spc != False and spc != 'link': name, ext = os.path.splitext(file) try: self.sp_energy, self.sp_program, self.sp_version_program, self.sp_solvation_model, self.sp_file, self.sp_charge, self.sp_empirical_dispersion, self.sp_multiplicity = parse_data( name + '_' + spc + ext) self.cpu = sp_cpu(name + '_' + spc + ext) except ValueError: self.sp_energy = '!' pass else: self.sp_energy, self.sp_program, self.sp_version_program, self.sp_solvation_model, self.sp_file, self.sp_charge, self.sp_empirical_dispersion, self.sp_multiplicity = parse_data( file) if self.sp_program == 'Gaussian' or self.program == 'Gaussian': # Count number of links for line in g_output: # Only read first link + freq not other link jobs if "Normal termination" in line: linkmax += 1 else: frequency_wn = [] if 'Frequencies --' in line: freqloc = linkmax # Iterate over output if freqloc == 0: freqloc = len(g_output) for i, line in enumerate(g_output): # Link counter if "Normal termination" in line: link += 1 # Reset frequencies if in final freq link if link == freqloc: frequency_wn = [] im_frequency_wn = [] if mm_freq_scale_factor is not False: fract_modelsys = [] # If spc specified will take last Energy from file, otherwise will break after freq calc if not g4: if link > freqloc: break # Iterate over output: look out for low frequencies if line.strip().startswith('Frequencies -- '): if mm_freq_scale_factor is not False: newline = g_output[i + 3] all_freqs = [] for j in range(2,5): try: fr = float(line.strip().split()[j]) all_freqs.append(fr) except IndexError: pass most_low_freq = min(all_freqs) for j in range(2, 5): try: x = float(line.strip().split()[j]) # If given MM freq scale factor fill the fract_modelsys array: if mm_freq_scale_factor is not False: y = float(newline.strip().split()[j]) / 100.0 y = float('{:.6f}'.format(y)) else: y = 1.0 # Only deal with real frequencies if x > 0.00: frequency_wn.append(x) if mm_freq_scale_factor is not False: fract_modelsys.append(y) # Check if we want to make any low lying imaginary frequencies positive elif x < -1 * im_freq_cutoff: if invert is not False: if invert == 'auto': if "TSFreq" in self.job_type: if x == most_low_freq: im_frequency_wn.append(x) else: frequency_wn.append(x * -1.) inverted_freqs.append(x) else: frequency_wn.append(x * -1.) inverted_freqs.append(x) elif x > float(invert): frequency_wn.append(x * -1.) inverted_freqs.append(x) else: im_frequency_wn.append(x) else: im_frequency_wn.append(x) except IndexError: pass # For QM calculations look for SCF energies, last one will be the optimized energy elif line.strip().startswith('SCF Done:'): self.scf_energy = float(line.strip().split()[4]) elif line.strip().startswith('E2('): spe_value = line.strip().split()[-1] self.scf_energy = float(spe_value.replace('D','E')) # For Counterpoise calculations the corrected energy value will be taken elif line.strip().startswith('Counterpoise corrected energy'): self.scf_energy = float(line.strip().split()[4]) # For MP2 calculations replace with EUMP2 elif 'EUMP2 =' in line.strip(): self.scf_energy = float((line.strip().split()[5]).replace('D', 'E')) # For ONIOM calculations use the extrapolated value rather than SCF value elif "ONIOM: extrapolated energy" in line.strip(): self.scf_energy = (float(line.strip().split()[4])) # For G4 calculations look for G4 energies (Gaussian16a bug prints G4(0 K) as DE(HF)) --Brian modified to work for G16c-where bug is fixed. elif line.strip().startswith('G4(0 K)'): self.scf_energy = float(line.strip().split()[2]) self.scf_energy -= self.zero_point_corr #Remove G4 ZPE elif line.strip().startswith('E(ZPE)='): #Overwrite DFT ZPE with G4 ZPE self.zero_point_corr = float(line.strip().split()[1]) # For TD calculations look for SCF energies of the first excited state elif 'E(TD-HF/TD-DFT)' in line.strip(): self.scf_energy = float(line.strip().split()[4]) # For Semi-empirical or Molecular Mechanics calculations elif "Energy= " in line.strip() and "Predicted" not in line.strip() and "Thermal" not in line.strip() and "G4" not in line.strip(): self.scf_energy = (float(line.strip().split()[1])) # Look for thermal corrections, paying attention to point group symmetry elif line.strip().startswith('Zero-point correction='): self.zero_point_corr = float(line.strip().split()[2]) # Grab Multiplicity elif 'Multiplicity' in line.strip(): try: self.mult = int(line.split('=')[-1].strip().split()[0]) except: self.mult = int(line.split()[-1]) # Grab molecular mass elif line.strip().startswith('Molecular mass:'): molecular_mass = float(line.strip().split()[2]) # Grab rational symmetry number elif line.strip().startswith('Rotational symmetry number'): if not ssymm: symmno = int((line.strip().split()[3]).split(".")[0]) # Grab point group elif line.strip().startswith('Full point group'): if line.strip().split()[3] == 'D*H' or line.strip().split()[3] == 'C*V': linear_mol = 1 # Grab rotational constants elif line.strip().startswith('Rotational constants (GHZ):'): try: self.roconst = [float(line.strip().replace(':', ' ').split()[3]), float(line.strip().replace(':', ' ').split()[4]), float(line.strip().replace(':', ' ').split()[5])] except ValueError: if line.strip().find('********'): linear_warning = True self.roconst = [float(line.strip().replace(':', ' ').split()[4]), float(line.strip().replace(':', ' ').split()[5])] # Grab rotational temperatures elif line.strip().startswith('Rotational temperature '): rotemp = [float(line.strip().split()[3])] elif line.strip().startswith('Rotational temperatures'): try: rotemp = [float(line.strip().split()[3]), float(line.strip().split()[4]), float(line.strip().split()[5])] except ValueError: rotemp = None if line.strip().find('********'): linear_warning = True rotemp = [float(line.strip().split()[4]), float(line.strip().split()[5])] if "Job cpu time" in line.strip(): days = int(line.split()[3]) + self.cpu[0] hours = int(line.split()[5]) + self.cpu[1] mins = int(line.split()[7]) + self.cpu[2] secs = 0 + self.cpu[3] msecs = int(float(line.split()[9]) * 1000.0) + self.cpu[4] self.cpu = [days, hours, mins, secs, msecs] if self.sp_program == 'NWChem' or self.program == 'NWChem': print("Parsing NWChem output...") # Iterate for i,line in enumerate(g_output): #scanning for low frequencies... if line.strip().startswith('P.Frequency'): newline=g_output[i+3] for j in range(1,7): try: x = float(line.strip().split()[j]) y = 1.0 # Only deal with real frequencies if x > 0.00: frequency_wn.append(x) if mm_freq_scale_factor is not False: fract_modelsys.append(y) # Check if we want to make any low lying imaginary frequencies positive elif x < -1 * im_freq_cutoff: if invert is not False: if x > float(invert): frequency_wn.append(x * -1.) inverted_freqs.append(x) else: im_frequency_wn.append(x) else: im_frequency_wn.append(x) except IndexError: pass # For QM calculations look for SCF energies, last one will be the optimized energy elif line.strip().startswith('Total DFT energy ='): self.scf_energy = float(line.strip().split()[4]) # Look for thermal corrections, paying attention to point group symmetry elif line.strip().startswith('Zero-Point'): self.zero_point_corr = float(line.strip().split()[8]) # Grab Multiplicity elif 'mult ' in line.strip(): try: self.mult = int(line.split()[1]) except: self.mult = 1 # Grab molecular mass elif line.strip().find('mol. weight') != -1: molecular_mass = float(line.strip().split()[-1][0:-1]) # Grab rational symmetry number elif line.strip().find('symmetry #') != -1: if not ssymm: symmno = int(line.strip().split()[-1][0:-1]) # Grab point group elif line.strip().find('symmetry detected') != -1: if line.strip().split()[0] == 'D*H' or line.strip().split()[0] == 'C*V': linear_mol = 1 # Grab rotational constants (convert cm-1 to GHz) elif line.strip().startswith('A=') or line.strip().startswith('B=') or line.strip().startswith('C=') : print(line.strip().split()[1]) letter=line.strip()[0] h = 0 if letter == 'A': h = 0 elif letter == 'B': h = 1 elif letter == 'C': h = 2 roconst[h]=float(line.strip().split()[1])*29.9792458 rotemp[h]=float(line.strip().split()[4]) if "Total times" in line.strip(): days = 0 hours = 0 mins = 0 secs = line.strip().split()[3][0:-1] msecs = 0 self.cpu = [days,hours,mins,secs,msecs] self.inverted_freqs = inverted_freqs # Skip the calculation if unable to parse the frequencies or zpe from the output file if hasattr(self, "zero_point_corr") and rotemp: cutoffs = [s_freq_cutoff for freq in frequency_wn] # Translational and electronic contributions to the energy and entropy do not depend on frequencies u_trans = calc_translational_energy(temperature) s_trans = calc_translational_entropy(molecular_mass, conc, temperature, solv) s_elec = calc_electronic_entropy(self.mult) # Rotational and Vibrational contributions to the energy entropy if len(frequency_wn) > 0: zpe = calc_zeropoint_energy(frequency_wn, freq_scale_factor, fract_modelsys) u_rot = calc_rotational_energy(self.zero_point_corr, symmno, temperature, linear_mol) u_vib = calc_vibrational_energy(frequency_wn, temperature, freq_scale_factor, fract_modelsys) s_rot = calc_rotational_entropy(self.zero_point_corr, linear_mol, symmno, rotemp, temperature) # Calculate harmonic entropy, free-rotor entropy and damping function for each frequency Svib_rrho = calc_rrho_entropy(frequency_wn, temperature, freq_scale_factor, fract_modelsys) if s_freq_cutoff > 0.0: Svib_rrqho = calc_rrho_entropy(cutoffs, temperature, freq_scale_factor, fract_modelsys) Svib_free_rot = calc_freerot_entropy(frequency_wn, temperature, freq_scale_factor, fract_modelsys,file, inertia, self.roconst) S_damp = calc_damp(frequency_wn, s_freq_cutoff) # check for qh if QH: Uvib_qrrho = calc_qRRHO_energy(frequency_wn, temperature, freq_scale_factor) H_damp = calc_damp(frequency_wn, H_FREQ_CUTOFF) # Compute entropy (cal/mol/K) using the two values and damping function vib_entropy = [] vib_energy = [] for j in range(0, len(frequency_wn)): # Entropy correction if QS == "grimme": vib_entropy.append(Svib_rrho[j] * S_damp[j] + (1 - S_damp[j]) * Svib_free_rot[j]) elif QS == "truhlar": if s_freq_cutoff > 0.0: if frequency_wn[j] > s_freq_cutoff: vib_entropy.append(Svib_rrho[j]) else: vib_entropy.append(Svib_rrqho[j]) else: vib_entropy.append(Svib_rrho[j]) # Enthalpy correction if QH: vib_energy.append(H_damp[j] * Uvib_qrrho[j] + (1 - H_damp[j]) * 0.5 * GAS_CONSTANT * temperature) qh_s_vib, h_s_vib = sum(vib_entropy), sum(Svib_rrho) if QH: qh_u_vib = sum(vib_energy) else: zpe, u_rot, u_vib, qh_u_vib, s_rot, h_s_vib, qh_s_vib = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 # The D3 term is added to the energy term here. If not requested then this term is zero # It is added to the SPC energy if defined (instead of the SCF energy) if spc is False: self.scf_energy += d3_term else: self.sp_energy += d3_term # Add terms (converted to au) to get Free energy - perform separately # for harmonic and quasi-harmonic values out of interest self.enthalpy = self.scf_energy + (u_trans + u_rot + u_vib + GAS_CONSTANT * temperature) / J_TO_AU self.qh_enthalpy = 0.0 if QH: self.qh_enthalpy = self.scf_energy + (u_trans + u_rot + qh_u_vib + GAS_CONSTANT * temperature) / J_TO_AU # Single point correction replaces energy from optimization with single point value if spc is not False: try: self.enthalpy = self.enthalpy - self.scf_energy + self.sp_energy except TypeError: pass if QH: try: self.qh_enthalpy = self.qh_enthalpy - self.scf_energy + self.sp_energy except TypeError: pass self.zpe = zpe / J_TO_AU self.entropy = (s_trans + s_rot + h_s_vib + s_elec) / J_TO_AU self.qh_entropy = (s_trans + s_rot + qh_s_vib + s_elec) / J_TO_AU # Symmetry - entropy correction for molecular symmetry if ssymm: sym_entropy_correction, pgroup = self.sym_correction(file.split('.')[0].replace('/', '_')) self.point_group = pgroup self.entropy += sym_entropy_correction self.qh_entropy += sym_entropy_correction # Calculate Free Energy if QH: self.gibbs_free_energy = self.enthalpy - temperature * self.entropy self.qh_gibbs_free_energy = self.qh_enthalpy - temperature * self.qh_entropy else: self.gibbs_free_energy = self.enthalpy - temperature * self.entropy self.qh_gibbs_free_energy = self.enthalpy - temperature * self.qh_entropy if cosmo: self.cosmo_qhg = self.qh_gibbs_free_energy + cosmo self.im_freq = [] for freq in im_frequency_wn: if freq < -1 * im_freq_cutoff: self.im_freq.append(freq) self.frequency_wn = frequency_wn self.im_frequency_wn = im_frequency_wn self.linear_warning = linear_warning # Get external symmetry number def ex_sym(self, file): coords_string = self.xyz.coords_string() coords = coords_string.encode('utf-8') c_coords = ctypes.c_char_p(coords) # Determine OS with sys.platform to see what compiled symmetry file to use platform = sys.platform if platform.startswith('linux'): # linux - .so file path1 = sharepath('symmetry_linux.so') newlib = 'lib_' + file + '.so' path2 = sharepath(newlib) copy = 'cp ' + path1 + ' ' + path2 os.popen(copy).close() symmetry = ctypes.CDLL(path2) elif platform.startswith('darwin'): # macOS - .dylib file path1 = sharepath('symmetry_mac.dylib') newlib = 'lib_' + file + '.dylib' path2 = sharepath(newlib) copy = 'cp ' + path1 + ' ' + path2 os.popen(copy).close() symmetry = ctypes.CDLL(path2) elif platform.startswith('win'): # windows - .dll file path1 = sharepath('symmetry_windows.dll') newlib = 'lib_' + file + '.dll' path2 = sharepath(newlib) copy = 'copy ' + path1 + ' ' + path2 os.popen(copy).close() symmetry = ctypes.cdll.LoadLibrary(path2) symmetry.symmetry.restype = ctypes.c_char_p pgroup = symmetry.symmetry(c_coords).decode('utf-8') ex_sym = pg_sm.get(pgroup) # Remove file if platform.startswith('linux'): # linux - .so file remove = 'rm ' + path2 os.popen(remove).close() elif platform.startswith('darwin'): # macOS - .dylib file remove = 'rm ' + path2 os.popen(remove).close() elif platform.startswith('win'): # windows - .dll file handle = symmetry._handle del symmetry ctypes.windll.kernel32.FreeLibrary(ctypes.c_void_p(handle)) remove = 'Del /F "' + path2 + '"' os.popen(remove).close() return ex_sym, pgroup def int_sym(self): self.xyz.get_connectivity() cap = [1, 9, 17] neighbor = [5, 6, 7, 8, 14, 15, 16] int_sym = 1 for i, row in enumerate(self.xyz.connectivity): if self.xyz.atom_nums[i] != 6: continue As = np.array(self.xyz.atom_nums)[row] if len(As == 4): neighbors = [x for x in As if x in neighbor] caps = [x for x in As if x in cap] if (len(neighbors) == 1) and (len(set(caps)) == 1): int_sym *= 3 return int_sym def sym_correction(self, file): ex_sym, pgroup = self.ex_sym(file) int_sym = self.int_sym() #override int_sym int_sym = 1 sym_num = ex_sym * int_sym sym_correction = (-GAS_CONSTANT * math.log(sym_num)) / J_TO_AU return sym_correction, pgroup