Mod:Hunt Research Group/plot potenital energy surface

Example PES produced using plot_PES.py

UPDATED 15/08/25

This script plots a 1D potential energy surface directly from a .log 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: link

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

x_label = "Dihedral Angle (degrees)"
plot_title = "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 starting_coordinate input_filename_1 input_filename_2
• where you list all of the files that make up your full surface
• 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

You can find an example .log file here: link (starting dihedral angle = 105 degrees)

IMPORTANT NOTES:

• Starting coordinate is the initial bond length/angle/dihedral angle you started the scan from
• Your input files are your .log files from your Gaussian scan
• 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 required modules and libraries
import sys
import os
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import PES_custom as pes

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

# process the input from the user
dir = os.getcwd()

if len(sys.argv) < 3:
    # checks to see that the user has supplied enough arguments when calling the script
    print("To run the script please type: python PES_plotter.py starting_coordinate input_file_name_1 input_file_name_2", 
         " \( and so on for as many .log files as you have)")
    sys.exit()
    
data_files = []
for i in range(len(sys.argv)-2):
    data_files.append(str(sys.argv[i+2])+'.log')


# setting up 
PES_data = pd.DataFrame(columns=['Coordinates','Energy'])
starting_coord = float(sys.argv[1])
supported_scans = {'B':'bond_length', 'A':'bond_angle', 'D':'dihedral_angle'}
angle_scans = ["A", "D"]
input_nextline = False


# defining search strings
input_str = "The following ModRedundant input section has been read:"
step_no_str = "Step number"
energy_str = "SCF Done:"


for data_file in data_files:
    print("Gathering data from", data_file)
    file = open(data_file,'r')
    line = file.readline()
    energys = np.array([])
    iterations = np.array([], dtype=int)
    step_numbers = np.array([], dtype=int)
    
    # Steps through each line in the file, checking if any search strings appear
    while line:
        if input_nextline:
            # Gathers baseline information about the PES
             temp = line.rstrip().split()
             scan_type = temp[0]
             no_steps = int(float(temp[-2]))
             no_structures = no_steps + 1
             step_size = float(temp[-1])
             input_nextline = False
             
        
             if scan_type 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")
        
        if input_str in line:
             input_nextline = True
             
        if step_no_str in line:
            # Identifies which scan and optimisation step we are up to
            temp = line.rstrip().split()
            iterations = np.append(iterations,int(float(temp[2])))
            step_numbers = np.append(step_numbers,int(float(temp[-4])))
            
        if energy_str in line:
            # Records the energy of current step
            temp = line.rstrip().split()
            energys = np.append(energys,float(temp[4]))
        
        # Moves to the next line
        line = file.readline()
    
    print(no_structures, "structures identified, recording optimised structure energies")
    
    # Selects the final structure for each scan step and records its energy and x-coordinate
    for i in range(no_structures):
        coord = starting_coord - i * step_size
        mask = step_numbers != i + 1
        step_iterations = iterations.copy()
        step_iterations[mask] = 0
        energy_index = list(step_iterations).index(np.max(step_iterations))
        eng = energys[energy_index]
        df = pd.DataFrame([[coord, eng]], columns=['Coordinates','Energy'])
        PES_data = pd.concat([PES_data, df], ignore_index=True)
    print("Structure energies identified")
        
        
# If True, adjusts the scan angles to be centered about 0    
if pes.centre_on_zero and scan_type in angle_scans: 
    print("Adjusting angles to range between -180 degrees and 180 degrees")
    PES_data["Coordinates"] = PES_data["Coordinates"].apply(center_angles)
    

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

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


# Plots the potential energy surface
print("Plotting PES...")
ax = plt.axes()    
ax.plot(sorted_PES_data["Coordinates"], 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)

# Saves the PES as an image if user has specified, otherwise plots the figure.
if pes.save_fig:
    print("Plotted PES saved as", pes.fig_filename)
    plt.savefig(pes.fig_filename)
else: plt.show()


sys.exit()