Mod:Hunt Research Group/plot potenital energy surface
Jump to navigation
Jump to search
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()
