plot_en_angle_gaussian_scan.py

#!/usr/bin/env python3
"""
Receives a Gaussian's ".log" of a scan generated by plot_eff_tors and print the energy vs. torsional curve.

Author: Henrique Musseli Cezar
Date: APR/2018
"""

import argparse
import matplotlib as mpl
# Force matplotlib to not use any Xwindows backend.
mpl.use('Agg')
import matplotlib.pyplot as plt
from distutils.spawn import find_executable

def parse_en_log_gaussian(fname):
  died = []
  ener = []
  conv = True
  with open(fname, 'r') as f:
    while True:
      line = f.readline()
      if line.strip().startswith("#"):
        break

    # if MP2 get the energy corrected by the perturbation
    if "MP2" in line.upper():
      mp2 = True
    else:
      mp2 = False

    for line in f:
      if "Convergence criterion not met" in line:
        conv = False
        
      if " dihedral =" in line:
        angle = line.split("=")[1].strip()

      # get energy and convert to kcal/mol
      if not mp2 and "SCF Done: " in line:
        if conv:
          en = float(line.split()[4])*627.509
          died.append(float(angle.split()[0]))
          ener.append(en)
        else:
          conv = True
      elif mp2 and "EUMP2" in line:
        en = float(line.split()[5].replace("D","E"))*627.509
        died.append(float(angle.split()[0]))
        ener.append(en)

  return sorted(died), [x for _,x in sorted(zip(died,ener))]


if __name__ == '__main__':
  parser = argparse.ArgumentParser(description='Receives a Gaussians ".log" of a scan generated by plot_eff_tors and print the energy vs. torsional curve.')
  parser.add_argument("logfile", help="Gaussian's .log file")
  args = parser.parse_args()

  died, ener = parse_en_log_gaussian(args.logfile)

  # print to screen
  for ang, en in zip(died,ener):
    print("%f\t%f" % (ang, en))

  # plot it
  if find_executable('latex') and find_executable('dvipng'):
    mpl.rcParams.update({'font.size':18, 'text.usetex':True, 'font.family':'serif', 'ytick.major.pad':4})
  else:
    mpl.rcParams.update({'font.size':18, 'font.family':'serif', 'ytick.major.pad':4})

  plt.plot(died,ener)
  plt.xlim([-180,180])
  plt.xticks([-180,-120,-60,0,60,120,180])

  # some example options
  # plt.xlabel(r"$\phi$ ($^\circ$)")
  # plt.ylabel(r"Estimated effective torsional (kcal/mol)")
  # plt.xlim([-180.0,180.0])
  # plt.xticks([-180,-120,-60,0,60,120,180])

  plt.savefig("scan_gaussian.eps", bbox_inches='tight')