Mod:Hunt Research Group/plot potenital energy surface

Example PES produced using plot_PES.py

Extracting Scan Data From GaussView

Currently this script requires you first extract the data from your scans from GaussView. For each calculation making up your overall scan:

• In GaussView open your .log file
• Go to Results then Scan to open the Gaussian formatted PES
• Right click on the plot and select Save Data...
• Save your data as a .txt file

Running plot_PES.py

This script is broken into 2 files. The first is a customisation file that you can edit to change the title, x axis label, plot colour, and save location of the produced figure. You can find a list of colours and their names here: colours

• Copy the below code into a script called "PES_custom.py"

x_label = "Dihedral Angle (degrees)"
plot_title = "Potenital Energy Scan over Rotation of Imidazolium Ethyl Group"
plot_colour = "black"
centre_on_zero = True
save_fig = False
fig_filename = "test_save.png"

The second file is the main script you will call. This needs to be in the same place as the customisation file you just created.

• Copy the below code into a script called "plot_PES.py"
• To execute type: python plot_PES.py scan_type input_filename_1 input_filename_2
• Depending on your customisation file, this will either display a figure with your PES plot or save the plot under the file name you gave in PES_custom.py

IMPORTANT NOTES:

• Replace scan_type with the type of scan that you did (e.g bond_angle, bond_length, dihedral_angle)
• Your input files .txt files you saved from Gaussian
• You can list as many input files as you need, just please omit the file extensions
• You can choose whether or not to save your figure as a .png automatically by setting save_fig = True in PES_custom.py
• If you do not automatically save your figure, the script will only stop running when you close the figure
• Any issues flick me an email at laura.doyle@vuw.ac.nz :)

import sys
import os
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import PES_custom as pes

dir = os.getcwd()

supported_scans = ['bond_length', 'bond_angle', 'dihedral_angle']

if len(sys.argv) < 3:
    print("To run the script please type: python PES_plotter.py scan_type input_file_name_1 input_file_name_2", 
         " \( and so on for as many .txt files as you have)")
    sys.exit()
    
else:
    if str(sys.argv[1]) not in supported_scans:
        print("WARNING: This scan type may not be supported, please double check results.",
              " Supported scan types are:", *supported_scans, sep="\n")

    data_files = []
    for i in range(len(sys.argv)-2):
        data_files.append(str(sys.argv[i+2])+'.txt')
    
scan_type = sys.argv[1] 


def center_angles(angle):
    # Adjust scan angles to be between -180 and 180 degrees
    centered = ((angle + 180) % 360) - 180
    return centered

# Grab the data from the scan data files and combine into a single dataframe
PES_data = pd.DataFrame()
for data_file in data_files:
    df = pd.read_csv(data_file,sep='\s+',skiprows=4,names=["Angle", "Energy"])
    PES_data = pd.concat([PES_data, df], ignore_index=True)
    
# If True, adjusts the scan angles to be centered about 0    
if pes.centre_on_zero and "angle" in scan_type: PES_data["Angle"] = PES_data["Angle"].apply(center_angles)

# Sorts the data by scanning variable
sorted_PES_data = PES_data.sort_values(by="Angle")

# Converts Energy from a.u. to kJ/mol and zeroes the energy to the minimum
energies_kj = sorted_PES_data["Energy"] * 2625
zeroed_energies_kj = energies_kj - np.min(energies_kj)


# Plots the potential energy surface
ax = plt.axes()    
ax.plot(sorted_PES_data["Angle"], zeroed_energies_kj, c=pes.plot_colour, marker="o", linewidth=0.75)
ax.set_xlabel(pes.x_label)
ax.set_ylabel("Relative Energy (kJ/mol)")
ax.set_title(pes.plot_title)

if pes.save_fig:
    plt.savefig(pes.fig_filename)
else: plt.show()


sys.exit()