Monkey Mind

Most of the time for Orca QM calculations I’ll use Weasel FACCT’s rather fantastic workflow tool but occasionally there are workflows which either don’t exist or don’t quite do exactly what I want. After attending the Virtual Winter School on Computational Chemistry and getting lessons on Dimitrios Liakos on compound scripts I’ve been mean to put the theory into practice. The talks on line here , if your new to Orca Frank’s and Anneke Dittmer’s talks are very good.

The example here is looking a C-H bond dissociation energies, which are quite useful to study sites of metabolism. Typically I use ML to get a rough idea of which bonds to look at, then follow up with QM. The writing of multiple orca scripts was getting a bit tedious, so to facilitate the calculations I’ve used Julia to write a compound script.

julia CH_bde_calc.jl <file.xyz> <index>

There are no imports as it’s just using base module.

filename = ARGS[1]
index = parse(Int, ARGS[2])
outfilename = "radical"

# Function to extract atomic coordinates from an XYZ file
function split_xyz(filename::String)
    file_contents = readlines(filename)
    atoms = String[]
    pattern = r"([A-Za-z]{1,2})\s+([\d.-]+)\s+([\d.-]+)\s+([\d.-]+)"
    for line in file_contents
        if occursin(pattern, line)
            push!(atoms, line)
        end
    end
    return atoms
end

# Get atom list from the XYZ file
arr = split_xyz(filename)

# Remove the atom at the specified index
new_arr = copy(arr)
deleteat!(new_arr, index)

I’ve just used PBE0 as the functional for speed. Variable are defined first, then there are three jobs; calculation of energy of the whole molecule then calculations of the fragment energies. Finally we read the property files and use the three energies to calculate the BDE.

# Define functional
functional = "PBE0"

# Convert arrays to newline-separated strings
parent_xyz = join(arr, "\n")
radical_xyz = join(new_arr, "\n")

# Output filenames
basename1 = "$(outfilename)_Compound_1"
basename2 = "$(outfilename)_Compound_2"
basename3 = "$(outfilename)_Compound_3"

# Prepare script for ORCA compound mode
script = """
%Compound
  # ------------------------------------
  # First job: Get energy of parent molecule
  # ------------------------------------
  Variable HartToKcal = 627.5096080305927;  # Hartree to kcal/mol conversion factor
  Variable molec = 0.0;
  Variable frag = 0.0;
  Variable Hatom = 0.0;
  Variable basename1 = "$basename1";
  Variable basename2 = "$basename2";
  Variable basename3 = "$basename3";
  Variable energyfrag, bde, deltaG;
  Variable myProperty = "DFT_Total_en";

  New_Step
    !$functional
    ! opt anfreq
    *xyz 0 1
$parent_xyz
    *
  Step_End

  # ------------------------------------
  # Second job: Get energy of radical
  # ------------------------------------
  New_Step
    !$functional
    ! opt anfreq
    *xyz 0 2
$radical_xyz
    *
  Step_End

  # ------------------------------------
  # Third job: Get energy of hydrogen atom
  # ------------------------------------
  New_Step
    !$functional
    *xyz 0 2
    H 0.0 0.0 0.0
    *
  Step_End

  # Read calculated properties from ORCA output
  molec.ReadProperty(propertyName=myProperty, filename=basename1);
  frag.ReadProperty(propertyName=myProperty, filename=basename2);
  Hatom.ReadProperty(propertyName=myProperty, filename=basename3);

  # Compute bond dissociation energy
  energyfrag = Hatom + frag;
  deltaG = energyfrag - molec;
  bde = HartToKcal * deltaG;

  print("BDE: %.12lf kcal/mol\\n", bde);
EndRun
"""

# Write script to output file
open("$outfilename.cmp", "w") do file
    write(file, script)
end

Using tetrachloroethane as a test BDE comes out at 100.52 kcal/mol which seems about right ball park.