Mod:Hunt Research Group/plot potenital energy surface

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"
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("Energy (kJ/mol)")
ax.set_title(pes.plot_title)

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


sys.exit()