Kinisot.py 文件源码

def __init__(self, file, isomer, *args):
      # Frequencies in waveunmbers
      self.frequency_wn = []

      # Extract the Force constants from a g09 logfile and generate the
      # mass-weighted Hessian matrix in Hartree/(amu Bohr^2)
      mw_hessmat = read_hess(file, isomer)

      # Convert from atomic units - a bit ugly
      unit_conversion = ENERGY_AU / (BOHR_RADIUS**2 * ATOMIC_MASS_UNIT) / ((SPEED_OF_LIGHT * 2 * np.pi)**2)
      eigs = np.linalg.eigvalsh(mw_hessmat * unit_conversion)
      freqs = [ np.copysign(np.sqrt(np.abs(freq)),freq) for freq in eigs ]

      # 5 or 6 small normal modes will be removed (depending on whether the molecule is linear or non-linear)
      if is_linear(file) == 'linear': trans_rot_modes = 5
      else: trans_rot_modes = 6

      # Keep a single imaginary frequency. It should be larger than the predefined cut-off
      if np.abs(freqs[0]) > freq_cutoff:
         self.im_frequency_wn = -1.0 * freqs[0]
         trans_rot_modes = trans_rot_modes + 1
      for freq in freqs[trans_rot_modes:]: self.frequency_wn.append(freq)

      # Calculate the excitation factor (EXC), the ZPE (ZPE) and Teller-Redlich product factor (PF)
      # returns a 1D-array of all terms
      self.PF = calc_product_factor(self.frequency_wn, freq_scale_factor)
      self.ZPE = calc_zpe_factor(self.frequency_wn, temperature, freq_scale_factor)
      self.EXC = calc_excitation_factor(self.frequency_wn, temperature, freq_scale_factor)