Monkey Mind

Mol2Mol in REINVENT
I recently attended UK QSAR event and Keishi Kohara from AstraZeneca gave a great talk on REINVENT4, using Mol2Mol transformer. I was vaguely aware of the option but had never used it. From Keishi’s talk it looked a lot simpler than the other “whack-a-molecule” options of iteratively editing input TOML files. Here is the relevant description in the paper.

Molecular optimization [22, 52]. A molecule is supplied to the generator as restraint. The generator will find a second molecule within a defined similarity. Depending on the similarity radius the molecule will be relatively similar to the supplied molecule but, importantly, the scaffold can change within the limits of the given similarity. (Mol2Mol)

Basically you just pick a pre-trained prior, a list of project molecules along with a few details and sample the model.

Install

Install is really easy, switch runtime to gpu/tpu, check what version of cuda is running (in this case 12.5). Set up conda for colab.
```
!nvcc --version
!pip install -q condacolab
import condacolab
condacolab.install_miniforge()
```
```
!git clone https://github.com/MolecularAI/REINVENT4.git
%cd REINVENT4
!python install.py cu125
```
That’s it, I wish installing other python packages was so easy! Next step is writing a toml file. We’ll create analogues around OSIMERTINIB where the smiles is saved in “osimertinib.smi”.
```
# REINVENT4 TOML input example for sampling
#

run_type = "sampling"
device = "cuda:0"  # set torch device e.g. "cpu"
json_out_config = "_sampling.json"  # write this TOML to JSON

[parameters]

## Mol2Mol: find molecules similar to the provided molecules
model_file = "priors/mol2mol_medium_similarity.prior"
smiles_file = "osimertinib.smi"  # 1 compound per line
sample_strategy = "beamsearch"  # multinomial or beamsearch (deterministic)
temperature = 1.0 # temperature in multinomial sampling
#tb_logdir = "tb_logs"  # name of the TensorBoard logging directory

output_file = 'sampling.csv'  # sampled SMILES and NLL in CSV format

num_smiles = 100000  # number of SMILES to be sampled, 1 per input SMILES
unique_molecules = true  # if true remove all duplicatesd canonicalize smiles
randomize_smiles = true # if true shuffle atoms in SMILES randomly
```
Running the code is very simple
```
!reinvent -l sampling.log sampling.toml
```
7th May 2025

BDE with orca compound script

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.

26th Mar 2025

compound, ORCA

QSAR with Chemprop
Chemprop is my favourite QSAR tool. It’s packed full of features, easy to use, and you can do some fairly sophisticated model building. Here are some notes on running predictions with evidential uncertainty.

Training data

I’ll use the in vogue Biogen dataset from Fang et al. 2023 which can be conveniently accessed via https://polarishub.io . As it’s a toy model I’ll assume the curated data is fine and won’t spend time finding a good data split, tailoring the descriptors, model parameters etc.
```
pip install git+https://github.com/chemprop/chemprop.git
pip install polaris-lib
pip install -U "ray[data,train,tune,serve]"
```
We’ll first generate the splits:
```
import pandas as pd
from sklearn.model_selection import train_test_split
import sys

df = pd.read_csv(sys.argv[1])

# Step 1: Split into training (80%) and remaining data (20%)
train_df, remaining_df = train_test_split(df, test_size=0.2)

# Step 2: Split remaining data into external validation (10%), and calibration (10%)
external_val_df, calibration_df = train_test_split(remaining_df, test_size=0.5)


# reset indices
train_df = train_df.reset_index(drop=True)
external_val_df = external_val_df.reset_index(drop=True)
calibration_df = calibration_df.reset_index(drop=True)

# Write the concatenated DataFrame to a CSV file
train_df.to_csv("training.csv", index=False)
external_val_df.to_csv("ext_val.csv", index=False)
calibration_df.to_csv("cal.csv",index=False)
```
Model building
```
chemprop hpopt --data-path training.csv -s "SMILES" --target-columns "AC_LOG SOLUBILITY PH 6.8 (ug/mL)" --task-type regression-evidential --search-parameter-keywords depth ffn_num_layers ffn_hidden_dim message_hidden_dim dropout --hpopt-save-dir hpoptv1 --multi-hot-atom-featurizer-mode ORGANIC -l evidential
```
Hyperparameter optimization is straightforward but time-consuming; if you’ve got the hardware, it can be sped up with --raytune-use-gpu --raytune-num-cpus 12. I found this very slow to do on an old MacBook Air M1 (>24 hours and still not complete), but Colab came to the rescue, and a free tier TPU instance worked really well (ca. 1 hr). The --multi-hot-atom-featurizer-mode flag with ORGANIC or, more recently, in v2.1.2, RIGR (Resonance Invariant Graph Representation) is more targeted and thus gives better results. At the time of writing, RIGR is a very new addition; so far, in my limited experience, the ORGANIC flag seems to slightly edge out the RIGR flag, but I expect this to be case dependent. From the pre-print, it does sound like the RIGR flag is particularly helpful for small data sets.

With the hyper-parameters identified we are ready to train the model. This can be done with the following command. Number of replicates is analogous to ‘outer loop’ of nested cross validation but at a lower cost, suitable for deep learning applications. Again if the hardware allows the training can be speeded up with --accelerator gpu. Depending on the task it important to specify the appropriate loss function (here evidential), apparently the method choice should pick this up but I am not sure it always does. With evidential regression the --evidential-regularization is model dependent, default is 0.0 and the default recommended by Soleimany et al. (2021) is 0.2
```
chemprop train --data-path training.csv -s "SMILES" --target-columns "LOG SOLUBILITY PH 6.8 (ug/mL)" --task-type regression-evidential --config-path ./hpoptv1/best_config.toml --multi-hot-atom-featurizer-mode ORGANIC --num-replicates 10 -o model_evidential --aggregation norm --epochs 50 -l evidential --evidential-regularization 0.2
```
Final step is get the predictions on the external validation set together with the uncertainty estimation which is calibrated with conformal-regression.
```
chemprop predict -s "SMILES" --test-path ext_val.csv --model-path model_evidential --preds-path preds.csv --multi-hot-atom-featurizer-mode ORGANIC --uncertainty-method evidential-total --calibration-method conformal-regression --cal-path cal.csv
```
The final results don’t look too impressive, but maybe not unsurprising; it’s a very unsophisticated model, the data range is pretty small with logS, it’s a tricky end point to predict, the experimental values can vary a lot depending on sample, the selected data isn’t well distributed over the experimental range and furthermore the applicability domain calculations suggest test set looks pretty different to train. Uncertainty estimations don’t seem to rescue the model either, or look calibrated… I’ll have to more flattering dataset next time but at least it look realistic. Next step is to see if it can be competitive TabPFN.
15th Mar 2025

chemprop, QSAR
GenChem with MoLeR and Molecular Swarm Optimisation
I came across MoLeR a few years ago at a conference in a talk by Krzysztof Maziarz, I liked the fact that it was fast to train, it used graphs rather than text to generate molecules and you could define a scaffold which is typically what you do in almost all projects.

There are loads of generative chemistry programs out there, but I thought MoLeR was particularly well put together and easy to use. For practical applications it was coupled with a molecular swarm optimisation (MSO). Both packages are very good but I had to do a bit of tweaking to get them to play nicely together. Code is available here . Compared to AZ’s Reinvent, I find the MSO + MoLeR combo easier to use, modify and gives less “extreme creative solutions” at least in my hands.

Step one is define is to train a model. This is very straightforward with the instructions on the MoLeR git.

Step two is to define the multidimensional box the swarm(s) will explore. This can be the chemistry used to build the MoLeR model, it can also represents a way of guiding the swarm by selecting a subset of chemistry. Here I encode the smiles with model and return the max and min of the encoding.
```
python define_max_min_of_model.py <model> <smiles>
```
Step three define some objectives for the swarm, initial smiles can be anything, here I use benzene but it can also be a swarm of different molecules. The scaffold is enforced and is present in every molecule. The MPO is a simple JSON dictionary. More details are on Robin Winter’s GitHub.
```
from mso.optimiser import BasePSOptimizer
import numpy as np
from molecule_generation import VaeWrapper
from mso.objectives.swarm_functions import swarm_wt
from mso.objectives.scoring import SwarmScoringFunction, ScoringFunction
from mso.utils import read_model
import time
import os
import sys


init_smiles = "c1ccccc1" # SMILES representation
scaffold = "OC=O"
mwt_desirability = [{"x" : 180, "y" : 0}, {"x":200, "y": 1}, {"x":400, "y": 1.0}, {"x":450, "y": 0.0}]
scoring_functions = [SwarmScoringFunction(func=swarm_wt, name="mwt", desirability = mwt_desirability)]
```
Step four is running the calculation. Here we have two swarms each with 500 probe molecules, the simple scoring function, the box dimensions, the model and the scaffold. Finally I run the model for 20 epochs and save the results to a csv.
```
model_dir, x_max, x_min = read_model(sys.argv[1])

with VaeWrapper(model_dir) as model:
    opt = BasePSOptimizer.from_query(init_smiles=init_smiles,num_part=500,num_swarms=2,inference_model=model,scoring_functions=scoring_functions, x_max=x_max, x_min=x_min, scaffold=scaffold)
    start_time = time.time()
    try:
        opt.run(20)
        timestr = time.strftime("%Y%m%d-%H%M%S")
        print("--- %s seconds ---" % (time.time() - start_time))
        try:
            opt.best_solutions.to_csv("best_solutions_" + timestr + ".csv", index=False)
            opt.best_fitness_history.to_csv("best_history_" + timestr + ".csv", index=False)
        except:
            pass
    except:
        opt.best_solutions.to_csv("best_solutions_" + "error" + ".csv", index=False)
        opt.best_fitness_history.to_csv("best_history_" + "error" + ".csv", index=False)
```
25th Feb 2025

GenChem, MoLeR, MSO

Model Uncertainty part2: Applicability Domains

Two simple approaches I like are the sum of distance-weighted contributions (SDC) and calculating the Mahalanobis distance from the training set.

Sum of distance-weighted contributions (SDC) is the easiest to calculate;

Where d_i is the distance between the ith molecule in the training set and the test molecule, a is adjustable parameter that modulates the distance penalty (typically 3) for molecule further away from the training set. Note that d is tanimoto distance, i.e. one minus tanimoto similarity.

def calc_sdc(ref_fps, query_fps):
    d = pariwise_distances(ref_fps, query_fps, metric="jaccard",n_job=-1, force_all_finite=True)
    sdc = np.sum(np.exp(-3*d/(1-d)),axis=0)

    return sdc

Fingerprints are ideally calculated by the same method as used to build the model.

Mahalanobis distance is definitely more difficult to say, and slightly trickier to calculate. It’s simply a measure of the distance between a point and a distribution . There are good implications for in python, Julia , etc. but something I thought would be useful to do would be to use the Facebook A.I. similarity search (FAISS) tool as it is super fast, easy to use, and so enabling application to large libraries.

def create_vector(vector):
    """
    Transform a vector to follow a unit Gaussian distribution using Mahalanobis transformation.
    """
    vector = np.ascontiguousarray(vector)
    xc = vector - vector.mean(0)
    cov = np.dot(xc.T, xc) / xc.shape[0]
    cov += 2 * np.eye(cov.shape[0])
    L = np.linalg.cholesky(cov)
    mahalanobis_transform = np.linalg.inv(L)
    return np.dot(vector, mahalanobis_transform.T)

With the vectors created from this function we create an in memory database of the transformed training vectors which we will then search.

def create_index(input_vector):
    """
    Create a FAISS index for a given set of vectors.
    """
    vector_dimension = input_vector.shape[1]
    index = faiss.IndexFlatL2(vector_dimension)
    index.add(input_vector)
    print(f"Index created with {index.ntotal} vectors")
    return index

Finally we search the indexed database and write the distances to file.

def main():
    if len(sys.argv) != 3:
        print("Usage: python AD_calc_mahalanobis.py <training_csv> <predict_csv>")
        sys.exit(1)

    training_df = pd.read_csv(sys.argv[1])
    predict_df = pd.read_csv(sys.argv[2])

    training_smiles = training_df["smiles"].tolist()
    predict_smiles = predict_df["smiles"].tolist()

    train_fps = smiles_to_fingerprint_array(training_smiles)
    predict_fps = smiles_to_fingerprint_array(predict_smiles)
    ###### SDC #############################################################


    sdc = calc_sdc(predict_fps, train_fps)

    predict_df["scd"] = sdc


    ######## mahalanobisDistance #################################

    array_train = create_vector(train_fps)
    array_predict = create_vector(predict_fps)


    index = create_index(array_train)
    distances, _ = index.search(array_predict, 1)

    predict_df["mahalanobisDistance"] = distances.flatten()

    ############# Write output ####################################
    predict_df.to_csv("test_cmps_with_distance.csv", index=False)
    print("Results saved to test_cmps_with_distance.csv")

if __name__ == "__main__":
    main()

To sanity check everything we will plot the SDC of the external test set, calibration set and a 10000 random smiles generated using MoLeR sampling a chembl derived model provided. The box-plot was generate with AlgebraOfGraphics.

using Colors
using AlgebraOfGraphics, CairoMakie
using DataFrames, CSV

update_theme!(
           fontsize=16,
           markersize=8,
           Axis=(title=L"Box Plot SDC",))

# define colour scheme 
Sun_Orange_Bright = RGB(252/255, 229/255, 0)
Plant_Green_Bright = RGB(80/255, 230/255, 0)
Air_Blue_Bright = RGB(0, 230/255, 190/255)
Energy_Pink_Bright = RGB(255/255, 0, 100/255)

#read data
df_1 = CSV.read(ARGS[1], DataFrame)
df_2 = CSV.read(ARGS[2], DataFrame)
df_3 = CSV.read(ARGS[3], DataFrame)


y1 = df_1[!, :scd]
y2 = df_2[!, :scd]
y3 = df_3[!, :scd]

x = vcat(fill("test", length(y1)), fill("cal", length(y2)), fill("random", length(y3)))  # Group labels
y = vcat(y1, y2, y3)  # Combine columns into one vector

# Create a named tuple for AlgebraOfGraphics
df = (; x, y)
plt = data(df) * mapping(:x=>L"x", :y =>L"y", color=:x) * visual(BoxPlot) * data(df) * mapping(:x=>L"x", :y =>L"y", color=:x) * visual(BoxPlot)

colors = [Plant_Green_Bright, Air_Blue_Bright, Energy_Pink_Bright]

fg = draw(plt, scales(Color = (; palette = colors)))

save("BoxPlot.svg",fg, px_per_unit = 3)

As expected test set and calibration set look similar, whereas the random smiles looks more different.

15th Feb 2025

AD, AlgebraOfGraphics, BoxPlot, SDC, uncertainty

Extracting uncertainty estimations from raw data.

I recently came across this blog article from Dr. Lewis Martin on Bayesian IC₅₀ fitting. The determination of uncertainty from raw data is something that is of interest to me. In the article Lewis use PyMC3, which is pretty popular, the analogues package in Julia is Turing so we will use that.

Imports

StaticArray.jl is mainly there to speed up the code, there maybe be better ways of doing this with threading or using a GPU.

using AlgebraOfGraphics, CairoMakie
using DataFrames
using  Turing, Distributions
using StaticArrays

Generate Fake Data

First define the concentration range for the x-axis.

# x points are in units of µM - they describe the 
# test ligand concentrations across a 9pt CRC
concs = [0.01, 0.032, 0.1, 0.32, 1, 3.2, 10, 32, 100]
# Convert concentrations to log10 Molar units
log_concs = -log10.(concs ./ 1e6)

Make up some pIC50 data for the y-axis using a sigmoid function.

function sigmoid(x; middle, slope, top, bottom)
    return top / (1 + ℯ^(-(x - middle) * slope) - bottom)
end

sigmoid_values = sigmoid.(log_concs; middle=6.1, slope=2, top=100, bottom=0)

We’ll check the data by plotting the data, here we will use Algebra of Graphics

# plot initial data.
df = (; log_concs, sigmoid_values, )
xy = data(df) * mapping(:log_concs => L"concentrations", :sigmoid_values => L"Activity")
layers = visual(Scatter)
fg = draw(layers * xy)
save("figure1.png", fg, px_per_unit = 3)

Next we’ll re-plot the data with duplicates and noise.

#add some duplicates and noise.
concs = repeat(log_concs, inner=2)        
responses = repeat(sigmoid_values, inner=2)  
# Add random noise (Normal distribution with mean 0, std dev 5)
responses .+= randn(length(responses)) .* 5

With the data in hand we ready for modelling, first we define a model. Priors for top and bottom are inherently set to 100 and 0. For the variance we use the inverse gamma, the alternative would be the truncated normal err ~ Truncated(Normal(0, 20), 0, Inf). After this we build the model and then do a “No U–Turn Sampling” with 10,000 datapoints.

# Define the Turing model
@model function sigmoid_model(concs, responses)
    concs = SVector{length(concs)}(concs)
    responses = SVector{length(responses)}(responses)
    # Priors
    s² ~ InverseGamma(0.5, 20)
    midpoint ~ Normal(5, sqrt(s²))
    slope ~ Normal(0, sqrt(s²))
    bottom ~ Normal(0, sqrt(s²))
    top ~ Normal(100, sqrt(s²))
    
    # Sigmoid function: 
    proportion = @. top / (1 + exp(-(concs + midpoint) * slope)) - bottom
    # Likelihood (observed data with Gaussian noise)
    responses ~ MvNormal(proportion, s²)
end
end

model = sigmoid_model(concs, responses)
chain = sample(model, NUTS(), 10_000; tune=3000, target_accept=0.9)

Finally plot the data as density plot, which looks good mean is approximately 6.1 and standard deviation is approximately 0.1.

df = DataFrame(chain)
df.midpoint = abs.(df.midpoint)
xy = data(df) * mapping(:midpoint => L"midpoint")
layers=AlgebraOfGraphics.density()
fg = draw(layers * xy)
save("figure3.png", fg, px_per_unit = 3)

1st Jan 2025

Bayesian, Julia, Turing

carboxylic acid pKa calculations
As we know, pKa measurements are an important physical property, and much effort has been devoted to their calculation. I often feel that the devil is in the detail, which I think some QSPR models can gloss over, as such I was drawn to this paper by Gerald Monard group which uses atomic QM charges to predict the pK_a of a range of simple acids with good accuracy (R² = 0.955). The authors do a reasonable job exploring solvent models, DFT methods and charge models.

Detailed here are my attempts to reproduce their work with the tools I have available. Let’s break it down into the different parts: a) QM calculations b) extracting the charges c) Model evaluation

Code can be found at this GitHub

a) QM calculations

The Gaussian input file can be found in the papers supporting info, however for various reasons ORCA is by far my favoured quantum chemistry software package (and free for non-commercial use). As such I have translated the Gaussian .com file into a Orca .inp file.
```
--Link1--
%chk=CAS_62-23-7.chk
# opt nosymm freq M06L/6-311G(d,p) SCRF=(SMD,solvent=water) pop=nbo IOP(6/80=1)

CAS_62-23-7

-1   1
    C   -2.9342860000       -0.0000210000       -0.0000290000     
    O   -3.4980180000       -1.1241640000        0.0014510000     
    O   -3.4980380000        1.1241130000       -0.0013430000     
    C   -1.4110770000       -0.0000090000       -0.0000170000     
    C   -0.7017960000       -1.2046730000        0.0021260000     
    C   -0.7018130000        1.2046660000       -0.0021430000     
    C    0.6827570000       -1.2173540000        0.0021950000     
    C    0.6827400000        1.2173660000       -0.0021760000     
    C    1.3557730000        0.0000110000        0.0000190000     
    H    1.2409380000        2.1448510000       -0.0037910000     
    H   -1.2544120000        2.1372820000       -0.0037420000     
    H   -1.2543830000       -2.1372960000        0.0037130000     
    H    1.2409670000       -2.1448320000        0.0038240000     
    N    2.8177720000        0.0000200000        0.0000390000     
    O    3.4015070000        1.0813120000       -0.0012430000     
    O    3.4015190000       -1.0812670000        0.0012640000     
```
If you are well dosed with caffeine, you will have noticed that I have neglected the to include the “nbo” keyword in the transcription, while there is a simple interface with ORCA the NBO program, it requires a licence. The omission does compromise the predictive power, as the authors do get a better correlation with experimental results when the natural population analysis atomic charges are used rather than the Löwdin or Mulliken charges. As an alternative I will use the Minimal Basis Iterative Stockholder (MBIS) [Ayers2016] charges which weren’t looked in the paper.
```
!B3LYP D4 DEF2-SVP OPT FREQ MBIS
%PAL
NPROCS 12 
END

%Method
WriteJSONPropertyfile True
End

%cpcm SMD TRUE
SMDSOLVENT "WATER"
end

* xyz -1   1
    C   -2.9342860000       -0.0000210000       -0.0000290000     
    O   -3.4980180000       -1.1241640000        0.0014510000     
    O   -3.4980380000        1.1241130000       -0.0013430000     
    C   -1.4110770000       -0.0000090000       -0.0000170000     
    C   -0.7017960000       -1.2046730000        0.0021260000     
    C   -0.7018130000        1.2046660000       -0.0021430000     
    C    0.6827570000       -1.2173540000        0.0021950000     
    C    0.6827400000        1.2173660000       -0.0021760000     
    C    1.3557730000        0.0000110000        0.0000190000     
    H    1.2409380000        2.1448510000       -0.0037910000     
    H   -1.2544120000        2.1372820000       -0.0037420000     
    H   -1.2543830000       -2.1372960000        0.0037130000     
    H    1.2409670000       -2.1448320000        0.0038240000     
    N    2.8177720000        0.0000200000        0.0000390000     
    O    3.4015070000        1.0813120000       -0.0012430000     
    O    3.4015190000       -1.0812670000        0.0012640000
*     
```
We run the calculations as normal, but given there are a few we just use a bash script to launch the calculations. Calculations are super speedy with ORCA6.0.1 and once we are done we have a xzy file and a property file in a json format (among many others).
```
for file in $(ls *.inp); do
   orca $file
done
```
b) Processing the output files

Things got a bit trickier here, I would be interested in hearing suggestions for a more elegant solution.

The first (Julia) script takes the outputted coordinate file (.xyz) and the property file (.json), uses the old & dependable obabel to write a mol2 file from the xyz, the script then update partial charges with those found in the ORCA json output. A key snippet is shown below.
```
function main(args)
    isempty(args) && error("Please provide a file name as an argument")
    file = args[1]
    base, _ = splitext(file)
    
    run(`obabel $base.xyz -omol2 -O $base.mol2`)
    atoms, bonds = read_mol2("$base.mol2")
    data = read_json("$base.property.json")
    new_charges = first.(data["Geometry_1"]["Mulliken_Population_Analysis"]["ATOMICCHARGES"])
    
    atoms_with_charges = add_charges(atoms, new_charges)
    output_file = "$base.mol2"
    write_mol2(output_file, atoms_with_charges, bonds)
end
```
The second script, grudgingly in python as I couldn’t workout an easy way of doing this in Julia, reads the newly modified mol2 file, locates the carboxylate with a SMARTS pattern and returns the partial charges for just the carboxylate group. Final we pivot the data and output a CSV file.

c) How did we do.

We did ok for quite a simple model, R² = 0.81 compared to R²=0.955 from the authors, the coefficient of determination is brought down by a couple of outliers (acetopyruvic acid and pyruvic acid) where the proton can H-bond to the adjacent oxygen. Pleasingly MBIS charges worked pretty well, and maybe unsurprisingly better than Mulliken charges with all their known problems. Maybe NPA charges would work better. Thanks again to the authors Zeynep Pinar Haslak, Sabrina Zareb, Ilknur Dogan, Viktorya Aviyente and Gerald Monard.
17th Nov 2024

ORCA, pKa
Making a pretty scatter plot with Algebra of Graphics
Inspired by a rather useful blog by the venerable Pat Walters on how to make pretty scatter plots using seaborn. I thought it would be useful to me to document how I could do this in Julia using Algebra of Graphics. What particular appeals to me about this package is how quickly it can construct beautiful (and complex) figures.

We will use the same data Pat used. First import the modules and tweak the theme.
```
using AlgebraOfGraphics, CairoMakie
using DataFrames
using CSV

update_theme!(
           fontsize=16,
           markersize=8,
           Axis=(title="hERG Model Performance",))
```
Next, we will read in the data, and extract out the logP, experimental and calculated IC50 values. It’s much better to work in pIC50 values so we convert these (it also makes plotting massively easier). A few things to quicker note for those are new to Julia; i) μ is just inputed as \mu ii) The .* 1e-6 is multiplying every element in the list “x” by 1e-6, similarly for -log10 function. iii) Julia uses quotes use ” ” rather than ‘ ‘ for stings.
```
df = CSV.read("seaborn_scatterplot_example.csv", DataFrame)
expt_IC50_μM = df[!, "Experimental IC50(uM)"]
predicted_IC50_μM = df[!, "Predicted IC50(uM)"]
logP = df[!,"LogP"]

function convert_pI50(x)
    _IC50_M = x .* 1e-6
    return -log10.(_IC50_M)
end

x = convert_pI50(expt_IC50_μM)
y = convert_pI50(predicted_IC50_μM)
```
It would be more elegant to use the input data frame and modify it in the mapping step, but I think this is harder to read. So we reassign the data frame with just what we want to plot. After some renaming of the x and y axis using LaTeX. We then compute a linear fit of y ~ 1 + x , and combine this with the scatter plot. The result is a very clean plot with minimal coding.
```
#take replace the dataframe with the modifed data
df = (; x, y, logP)
xy = data(df) * mapping(:x => L"experimental $pIC_{50}$", :y => L"predicted $pIC_{50}$", color=:logP => L"$logP_{o/w}$")
layers = linear() + visual(Scatter)
fg = draw(layers * xy)

save("figure.png", fg, px_per_unit = 3) 
```
16th Nov 2024

Imports

Generate Fake Data

a) QM calculations

b) Processing the output files

c) How did we do.