There is a review of Burette here

it can be use in four ways

• Finder previews: select a molecular file in Finder and press Space.
• Desktop workspace: open Burette directly to inspect files in tabs, browse project folders, search commands and structures, sketch molecules, review collections, and send files to external chemistry tools.
• iPhone preview app: build the BuretteMobile Xcode target from source, then open molecular documents from Files or another iOS document provider.
• Public molecular plugin: preview an authorized molecular attachment or a public PDB entry directly in the Burette workspace through ChatGPT or Codex.

Download and installation instructions are here https://avogadro.cc/install/index.html

This challenge focuses on predicting how small molecules bind to the Zika virus NS2B-NS3 protease. We invite the molecular modelling community to predict protein–ligand complexes for 400 compounds, spanning small fragments to larger, lead-like molecules. Our aim is not simply to produce a leaderboard, but to understand where methods succeed, where they fail, and how they can be improved

Having just updated the docking scoring Jupyter notebook I thought it might be interesting to provide a version for the blind challenge. This notebook could perhaps be modified to try different docking algorithms or alternative scoring functions.

As the website highlights structural studies have shown that the ZIKV NS2B-NS3 protease can adopt several conformations involving rearrangement of the NS2B cofactor: a closed conformation, an open conformation, and a super-open conformation. The open and super-open conformations are considered proteolytically inactive, as the C-terminal part of NS2B that contributes to the substrate-binding site is unstructured in these conformations.

The focus been on the closed conformation of NS2B-NS3, using a bicistronic vector to co-express the NS2B cofactor and the NS3 protease. This optimised system gave us crystals in the P4₃22 space group with only one NS2B-NS3 complex in the asymmetric unit and an accessible active site, ideal for fragment screening and ligand soaking (PDB code: 29HN). This structure has no ligand bound, but there are several other structures in the PDB that contain a bound ligand including 7H1J and 9RM4. All three structures were imported into MOE and the structures prepared. They were then superimposed as shown below.

A close up view of the binding site is shown below.

And the key ligand interactions are shown here.

This ligand was copied into the binding site of the superimposed structure 29HN and minimised. This was then exported in pdb format and used for the docking studies.

The notebook used for docking can be seen here, this version actually only docks a trial 10 compounds to test the process. In this html version the sliders to scroll through the poses don’t work but you can download the actual notebook at the bottom of the page.

After docking the poses are scored using AEV-Plig and the best scoring poses for each ligand is exported as both a pdb file and as a PyMOL session file.

As you can see from the plot below there is little correlation between Smina docking score and the AEV-plig score.

The notebook and associated files can all be downloaded from GitHub https://github.com/drc007/zika_docking

This work lead to the above publication which has had over 18,500 total views and been downloaded almost 5000 times.

There is also a lab manual and the old Jupyter notebook all of which are freely available on GitHub

https://github.com/UCL/Open_Docking_Lab_Handbook

Over the years I’ve updated the notebook to fix bugs and add features that students at different universities have requested, the result is a project that become a little unwieldy and I thought I’d spend a few weeks (with Claude) cleaning it all up and making it a lot more efficient.

This notebook does not include the preparation of the target protein for docking, but here are a few of the things that need to be addressed.

–Alternates, Residues with alternate locations and/or ambiguous sequence identities (choose highest occupancy) but remember proteins are not static.

–Termini, Protein chain C- or N-termini which need to be charged or capped, or if DNA the terminal PO4 may only have three oxygens bonded to the phosphorous and an additional oxygen needs to be added. Sometimes loops are very disordered and appear as a breaks in the chain, it may be possible to use a loop library to model a replacement.

–Hydrogens, usually not visible and so need to be added/checked, particularly check hydrogens on heteroatoms, especially active site residues where the local environment may influence pKa 

–Ligand, Novel ligands in particular need checking to confirm atoms and bond orders are correct (PDB files do not contain bond information)

–Conformation, check that torsions are reasonable and there are no clashes.

–Charge, It with worth checking the charge on all ionisable groups.

–It can be difficult to be certain of the position of nitrogens in His or the primary amide in Asn, Gln.

The new README for the docking and scoring is shown below

Docking_scoring

A Jupyter notebook for docking ligands into a protein with smina and rescoring the poses with the AEV-PLIG binding-affinity model. It covers conformer generation, binding-site preparation, rigid or flexible docking, pose viewing, affinity prediction using AEV-plig and export in a variety of formats.

The example is the SARS-CoV-2 main protease (Mpro) fragment structure 5RHD. A fully run copy of the notebook is in DockingScoring.html (download it and open it in a browser to see the outputs and interactive 3D views) and, for printing or quick viewing on GitHub, DockingScoring.pdf.

Workflow

1. Ligands. Load an SDF of 3D structures, check that the names are unique, and browse the molecules with a slider.
   1b. Ionisation states. Optionally predict each ligand’s ionisation state at a chosen pH with QupKake (micro-pKa prediction), and dock the ligands in those states.
2. Conformers. Generate conformers with RDKit srETKDGv3 (or ETKDGv3, ETKDG, ETDG), pruned by RMSD, MMFF94-minimised and filtered by an energy window. Undefined stereocentres can optionally be enumerated.
3. Binding site. Either:
   • load a protein–ligand complex PDB and pick the ligand that defines the site, or
   • supply prepared files: a protein without the ligand plus the extracted ligand.
4. Images. Ray-traced PyMOL images of the whole protein with the ligand bound, and of the ligand in the binding site.
5. Docking. An smina input panel for rigid docking, flexible side chains, or --minimize.
6. Results. Step through the poses in 3D with a slider.
7. AEV-PLIG. Predict the pK of every pose, and score the ligand from the PDB in its crystal pose as a reference.
8. Export. Write the best pose of each ligand as PDB files and as a PyMOL session.

Setup

Run these commands from the repository folder.

1. Conda environment and Jupyter kernel

conda env create -f environment.yml

conda run -n docking_scoring python -m ipykernel install --user --name docking_scoring --display-name "Python (docking_scoring)"

environment.yml is portable (macOS and Linux) and installs current versions. docking_scoring.yml pins the exact versions this was developed and tested with on an Apple-silicon Mac. It has never been tested under Windows and the best option might be to try Windows Subsystem for Linux .

Install PyTorch from conda-forge, not pip. The pip wheel bundles its own OpenMP library, which clashes with conda’s and crashes on import with OMP: Error #15.

2. AEV-PLIG trained models

The trained models (about 300 MB) are not stored in this repository. Download them from the official AEV-PLIG repository:

conda run -n docking_scoring python aev_plig/download_models.py

The script fetches the 10-model model_GATv2Net_ligsim90_fep_benchmark ensemble and its scaler from a fixed upstream commit, checks each file’s checksum, and saves them to aev_plig/output/trained_models/.

3. smina

smina is not included in environment.yml: the conda-forge smina package (2020 build) needs an old version of Boost that forces an old RDKit. Install it in one of two ways:

• A separate conda environment, which the notebook finds automatically: conda create -n smina -c conda-forge smina
• The prebuilt binary from https://sourceforge.net/projects/smina/. Make it executable (chmod +x), then put it on your PATH or set the SMINA environment variable to its full path.

The notebook looks for smina in this order:

1. the SMINA environment variable;
2. smina, smina.static or smina.osx.12 in the notebook’s environment or on the PATH;
3. a conda environment called smina;
4. /usr/local/smina.osx.12 and other common locations.

You can also type the path into the smina panel.

4. Run

Open DockingScoring.ipynb in JupyterLab with the Python (docking_scoring) kernel.

• Each step has a panel with a Run or Apply button. Running all cells sets the panels up but does not start the long jobs.
• To run every step with the default settings, set AUTO_RUN = True in the first code cell.

Project layout

DockingScoring.ipynb         the notebook
DockingScoring.html          a fully run copy of the notebook (static)
DockingScoring.pdf           the same, as a PDF (A4 landscape)
docking_tools.py             all helper code (panels, viewers, docking, scoring, export)
environment.yml              conda environment (portable)
docking_scoring.yml          exact tested environment (macOS arm64)
aev_plig/                    AEV-PLIG prediction code (models: aev_plig/download_models.py)
example/                     5RHD example inputs; results are written here when you run the notebook

Generated results are not tracked in git (see .gitignore). Run the notebook to recreate them, or see DockingScoring.html.

Inputs

Ligand SDF

• 3D structures, ideally with hydrogens.
• Each record needs a name on its first line.
• smina keeps only that name, and the results are grouped by it, so different structures must not share a name.

The notebook checks for duplicate names. If it finds any, it offers two fixes:

• use an SD field with a unique value for every molecule (e.g. ID, IDNUMBER, Record_ID) as the ligand ID, or
• assign unique names automatically (name, name_2, …).

The original name is kept in the original_name property, and the input file is not modified.

Binding site

Complex structure mode (PDB or mmCIF)

• Reads .pdb and mmCIF (.cif, .mmcif, optionally gzipped). mmCIF is the PDB’s standard format, and is read with gemmi.
• Lists all ligands, with chemical names where the file has them, and selects the largest automatically. The panel says which one was chosen; pick another from the list if the file holds several ligands.
• Hides waters, solvents and crystallisation additives (e.g. DMSO, glycerol, PEG, sulfate), ions, and terminal caps such as NME. Tick Also list solvents, ions, additives and caps to show them.
• Uses only the first alternate location (A) of the protein and ligand.
• Bond orders for the extracted ligand come from, in order of preference:
  1. a SMILES template, if you paste one. Its charges are kept, so it can set the protonation state;
  2. the ligand’s hydrogens and formal charges (RDKit), when hydrogens are present;
  3. the ligand’s entry in the PDB Chemical Component Dictionary, fetched from RCSB, for files without hydrogens such as RCSB X-ray structures. It is used only if its heavy atoms match, and gives the neutral molecule;
  4. Open Babel perception from the geometry.
• The protein file keeps the HELIX and SHEET records, but other HET groups and waters are removed. Metal ions are kept unless you untick that option.
• Output files are named after the input, e.g. 5RHD_protein.pdb. mmCIF input adds _cif (5RHD_cif_protein.pdb), so it never overwrites files made from a PDB file of the same entry, which may use a different coordinate frame.

Example: example/5RHD.pdb was prepared in MOE (hydrogens added, piperazine protonated, ligand named LIG). example/5RHD.cif is the RCSB file: its ligand US7 has no hydrogens and is taken as the neutral molecule. Paste CS(=O)(=O)c1ccc(N2CC[NH2+]CC2)cc1 as the SMILES template for the protonated form. Using the CIF gave the same best-pose AEV-PLIG pK for FromXRay (5.31) as the PDB file.

Prepared-files mode uses your protein and ligand PDB files as they are.

Note: The two modes can give different smina scores for the same ligand. For example, proteinNoligand.pdb keeps the waters and DMS molecules, which add contacts. Compare runs that used the same mode.

Outputs

Docked poses are stored as SDF files. The best pose of each ligand is also exported as PDB files (ligand only, and protein–ligand complex) and as a PyMOL session file.

Every output is written next to the input it was derived from.

All prediction results go to a predictions/ folder next to the docked SDF:

Copies of the docked SDF, receptor, reference ligand, conformers and images are also placed in predictions/, so the folder is a self-contained record.

Search box

smina docks only inside the search box: the extent of the autobox ligand (all atoms, hydrogens included) plus autobox_add (default 4 Å) on every side. Atoms outside the box are heavily penalised during the search, and smina gives no warning when a ligand is too large. The ligand is simply squeezed: it curls up, or groups that would reach a neighbouring sub-pocket point elsewhere.

Before docking, the smina panel shows the box size and compares each ligand’s longest heavy-atom span, taken from the conformer file, with the smallest box side:

When any ligand is tight or longer, the panel suggests an autobox_add. It also warns above Vina’s recommended 27,000 ų.

After docking, smina’s output reports how many poses have a heavy atom within 1 Å of a box face. The summary and best-poses tables add box_margin, the closest heavy-atom distance to a face in Å (negative means outside), and at_box_edge. A best pose at the box edge means the box may have limited that ligand: re-dock with a larger autobox_add, or set the box explicitly with --center_x/y/z --size_x/y/z under “extra arguments”.

Example (5RHD fragment + 4 Å = 16.8 × 13.2 × 15.5 Å):

• Phenyl (11.4 Å) and iProAmine (11.9 Å) are tight.
• 38 of 297 poses come within 1 Å of the edge, and 3 have an atom slightly outside.
• iProAmine’s best AEV-PLIG pose is 0.13 Å outside the box.

The box is stored in each run record (<output>.run.json).

Ligand efficiency

The summary and best-poses tables include ligand efficiency (LE), the binding free energy per heavy atom, in kcal/mol per heavy atom:

An LE of about 0.3 or more is the usual guide for a good starting point. Small fragments such as the 5RHD series (15–21 heavy atoms) typically score 0.4–0.5 on LE_AEV_PLIG.

Note: After flexible docking, LE_smina is inflated, because the flexible smina score appears to include the flexed side chains. Use LE_AEV_PLIG to compare across docking modes.

Ionisation states (QupKake)

QupKake (Abarbanel & Hutchison, J. Chem. Theory Comput. 2024, https://doi.org/10.1021/acs.jctc.4c00328) predicts micro-pKa values from GFN2-xTB calculations and graph neural networks. Step 1b of the notebook uses it to set the ionisation state of the input ligands and of the ligand from the PDB.

How the state is set:

1. Each ligand is neutralised and submitted to QupKake, which predicts a pKa for every basic and acidic site.
2. At the chosen pH (default 7.4), the most basic site is protonated if its pKa is above the pH, and the most acidic site is deprotonated if its pKa is below it. At most one of each is applied. Other sites that also cross the pH are reported, not applied.
3. Implausible basic sites are ignored (on by default; untick the box to use every site). For this ligand series QupKake reported sulphonamide N, amide N and aryl-amine N, such as the N-aryl nitrogen of an aryl-piperazine, as the most basic sites. Taken at face value, that would protonate a sulphonamide ([NH3+]S(=O)(=O)…).

The ligands in their predicted states are written to <input>_pH7.4.sdf, with a table in <input>_pH7.4.csv. When docking in the predicted states is ticked, that file is docked, so every later output is named <input>_pH7.4_… and never overwrites the results for the input states. The PDB ligand keeps its crystal coordinates and is used as the reference in the final table. QupKake’s own files are in qupkake/ next to the input.

Example result (pH 7.4). With the filter on, QupKake’s states match the MOE-prepared input exactly: the piperazine or piperidine NH is protonated, or the N-isopropyl N for iProAmine. Re-docking therefore reproduced all 297 poses and scores. The neutral US7 from 5RHD.cif is protonated to the same state.

With the filter off, 6 of the 7 ligands were charged on a sulphonamide or aryl-amine N, and the best-pose AEV-PLIG pK changed by −0.27 to +0.35:

Installing QupKake. It needs its own conda environment, named qupkake, with xtb 6.4.1 built from source; see the QupKake README. On Apple silicon, build xtb at -O2, because gfortran 16 at -O3 miscompiles it. The notebook finds the environment automatically and limits each xtb process to one thread. If QupKake is not installed, the input ionisation states are used.

Full details for installing QupKake are available here https://github.com/drc007/qupkake-apple-silicon/tree/main/qupkake-apple-silicon

Reproducibility

Every smina run writes a run record next to its output, <output>.run.json, and the same summary at the top of <output>.log. The record holds:

• the smina version and path, and the full command;
• seed, --cpu, exhaustiveness and every other setting;
• the platform;
• SHA-256 checksums of the receptor, ligand and autobox files and of the output.

AEV-PLIG adds predictions/<docked>_run_info.json, which contains that docking record plus the model, scoring options, library versions, receptor checksum and the PDB-ligand scores.

Seed and thread count. With the tested smina build (master:dc3dfab, based on AutoDock Vina 1.1.2), a fixed --seed gave identical poses whatever the --cpu value:

At --cpu 3 (fewer threads than exhaustiveness runs), flexible docking showed a run-to-run jitter of up to 0.002 kcal/mol in the reported score. That is below the 0.01 kcal/mol precision of smina’s log.

--cpu therefore defaults to the number of cores. Results can still differ with a different smina or Vina build, platform or compiler, so keep the run record with your results.

Flexible docking

• Choose Flexible side chains in the smina panel.
• Automatic selection takes every residue within a set distance of the ligand (default 5 Å).
• Manual selection lets you pick residues from a list, or type them (A:41,A:145).
• Gly, Ala and Pro are excluded, because their side chains have no rotatable bonds.

It is much slower than rigid docking. For the example’s 33 conformers:

• rigid docking took about 40 s;
• flexible docking with 4 residues at exhaustiveness 8 took about 8½ min.

The automatic 5 Å selection picks 13 residues for 5RHD, which would be far slower again.

smina writes the flexed side chains for every pose. They come out in the same order as the poses, but without PDB atom names. The notebook maps them back onto the receptor and rebuilds a complete receptor for each pose, which is used:

• in the 3D viewer, where the flexed side chains are magenta and the original ones grey;
• by AEV-PLIG, which scores each pose against its own receptor;
• in the PDB and PyMOL exports.

Pairing poses with their side chains is verified, not assumed. smina writes each ligand’s poses and side-chain models together, so if it skips a ligand both files lose it and stay in step (tested). Before anything is paired, the notebook checks:

• Names: the docked SDF’s record titles and the models’ names must match. If they don’t, ligand blocks are aligned by name, and a block is paired only when its pose and model counts agree.
• Geometry: a pose that overlaps its own flexed side chains (< 1.5 Å) is rejected.

AEV-PLIG stops with an error on any mismatch, rather than score a pose against the wrong receptor. The viewer shows such poses with the rigid protein and a warning.

Note: The flexible-docking smina score appears to include the flexed side chains, so its values (around −8 to −11 kcal/mol for the example) are not comparable with rigid scores or with the PDB-ligand row. AEV-PLIG pK values are comparable between modes.

AEV-PLIG

aev_plig/ contains the AEV-PLIG prediction code (process_and_predict.py and its modules) from https://github.com/oxpig/AEV-PLIG, licensed under the BSD 3-Clause licence in aev_plig/LICENSE.txt. The trained model_GATv2Net_ligsim90_fep_benchmark ensemble (10 models) is downloaded separately (see Setup).

Changes for the current library versions:

• torchani 2.9 changed its SpeciesConverter, so the species index is now computed directly as atomic number − 1 (equivalent).
• pkg_resources is replaced by importlib.metadata.
• Model weights are loaded with torch.load(..., weights_only=False), which torch ≥ 2.6 requires.
• A new --output_dir option controls where the predictions are written.
• Dataset paths are handled with os.path.splitext.
• RDKit’s deprecated GetExplicitValence() is replaced.

Regression check: on AEV-PLIG’s example dataset, all 241 predictions from the patched code in the new environment match the original code in a Python 3.8 / torch 2.4 environment to within 2×10⁻⁶.

Hydrogens

Before scoring, all hydrogens are added to each ligand pose, because AEV-PLIG was trained on fully protonated ligands and one of its atom features counts attached hydrogens. smina writes polar hydrogens only. Scoring the same poses with every hydrogen removed changes the predicted pK by up to about 1 unit.

Citation

Valsson et al., “Narrowing the gap between machine learning scoring functions and free energy perturbation using augmented data”, Communications Chemistry (2025), https://doi.org/10.1038/s42004-025-01428-y. The licence is in aev_plig/LICENSE.txt.

Display options

• Protein colouring uses the RasMol secondary-structure scheme by default: helix magenta-red, sheet yellow, coil and turn white. Coil is drawn light grey so it shows on a white background.
  • Set dt.PROTEIN_COLOURING to "rainbow" or "grey" for other schemes.
  • Set dt.RASMOL_COIL = "#FFFFFF" for strict RasMol white coil.
• Pose viewers have toggles for the binding-site residues, the ligand from the PDB, and (after flexible docking) the original side chains.
• Binding-site image has a toggle to hide the side chains.

Known limitations

• smina silently stops at an unreadable SDF record. Every later ligand is skipped, the exit code is still 0, and the log says nothing. The smina panel warns before docking if the ligand file has unreadable records. After docking it compares the input ligands with the docked ones, reports any ligand with no poses, and records them in the run record (ligands_missing).
• Running the notebook from the command line with jupyter nbconvert --execute stalls on the slider cells; this is an nbclient issue. All the steps still complete and their files are written. Use JupyterLab interactively.
• Turns are coloured as coil unless the PDB file has TURN records, because neither 3Dmol.js nor PyMOL assigns turns.
• AEV-PLIG accepts only the elements C, N, O, S, P, F, Cl, Br, I and B. Ligands with any other element are skipped, and the notebook reports how many were rejected.

Licence

The notebook, docking_tools.py and documentation are released under the MIT licence (see LICENSE) and are available on GitHub https://github.com/drc007/docking_scoring.

The bundled AEV-PLIG code in aev_plig/ keeps its own BSD 3-Clause licence (see aev_plig/LICENSE.txt). The 5RHD structure comes from the Protein Data Bank (https://www.rcsb.org/structure/5RHD).

High-throughput experimentation (HTE) accelerates molecular discovery but struggles to quantify yields without pure analyte standards. The current need to purify standards increases cost, complexity, and cycle time, limiting the scalability of new molecular discovery. Here, we introduce an integrated “detect-predict” framework that leverages artificial intelligence (AI) for analyte standard-free yield quantification.

Paper is here https://doi.org/10.1021/jacsau.6c00812

Code https://github.com/xtalpi-inc/lcms2yield

Safari Extension, it allows you so select a drug name in a web page (for example), right-click and from the dropdown menu there is the option to choose either PubChem or ChEMBL as shown below. Compound Lookup turns any compound name on a web page into a structure and a set of database identifiers, without leaving the page.

For Safari https://macinchem.org/2026/09/29/safari-extension/

For Chrome https://macinchem.org/2026/09/29/google-chrome-extension/

• The bundled xtb-641 binary is Linux x86-64 and can’t run on macOS.
• conda-forge xtb is rejected by QupKake on import (“Conda version of xTB is currently not supported”). The conda-forge build also fails at runtime with Library not loaded: @rpath/libgfortran.5.dylib.
• A source build of xtb 6.4.1 with current compilers builds and runs, but at the default -O3 it is silently miscompiled: geometry optimizations climb in energy and never converge (details under Pitfalls).

This guide gives a setup that works, tested on macOS 26.6 (Mac Studio, M-series, 24 cores) in October 2026.

Throughout this guide, QK is the directory the QupKake repo is cloned into:

export QK=$HOME/Projects/QupKake

1. Clone QupKake

git clone https://github.com/Shualdon/QupKake.git "$QK"

2. Create the conda environment

Use conda-forge only. The pytorch channel named in the repo’s environment.yaml no longer publishes builds. The same environment also supplies the Fortran toolchain for xtb, so the binary links against the env’s own libgfortran and openblas.

conda create -n qupkake -c conda-forge --override-channels -y \
  python=3.9 pip "pytorch=2.1" torchvision torchmetrics "pytorch-lightning=2.1" \
  pandas numpy scipy rdkit matplotlib jupyterlab ipykernel \
  meson ninja pkg-config gfortran openblas liblapack c-compiler
conda activate qupkake

pip install torch_geometric==2.4
pip install "$QK"

Torch is pinned to 2.1, as in the repo’s environment.yaml, to match the shipped model checkpoints.

3. Build xtb 6.4.1 from source

3a. Download and verify

mkdir -p "$QK/xtb-src" && cd "$QK/xtb-src"
curl -LO https://github.com/grimme-lab/xtb/releases/download/v6.4.1/xtb-6.4.1.tar.xz
curl -LO https://github.com/grimme-lab/xtb/releases/download/v6.4.1/xtb-6.4.1.tar.xz.sha256
shasum -a 256 xtb-6.4.1.tar.xz; cat xtb-6.4.1.tar.xz.sha256   # must match:
# 15cfd8743dd963c54e0fee12da5ee34a08dde78625cf8a0a21eee1eaa652beae
tar xf xtb-6.4.1.tar.xz && cd xtb-6.4.1

3b. Patch one format string

gfortran 16 rejects a malformed format descriptor at runtime: Fortran runtime error: Missing comma between descriptors. xtb hits it during every --opt run. The patch adds the missing commas:

sed -i '' '624s/(1x,"("f7.2"%)")/(1x,"(",f7.2,"%)")/' src/optimizer.f90
sed -n 624p src/optimizer.f90
# expected:       write(env%unit,'(1x,"(",f7.2,"%)")')         (depred-echng)/echng*100

3c. Configure and build at -O2 (not the default -O3)

FC=gfortran CC=cc meson setup build-o2 \
  --prefix="$QK/xtb-6.4.1" \
  --buildtype=debugoptimized \
  -Dla_backend=openblas -Dopenmp=true
ninja -C build-o2 install

• debugoptimized means -O2 -g. Do not use --buildtype=release. It gives -O3, which gfortran 16 miscompiles.
• la_backend=openblas: the default is mkl-static, which doesn’t exist on Apple Silicon.
• meson prints WARNING: FC and CC are not from the same vendor. That’s expected and harmless.

The installed binary finds libgfortran, libopenblas and libomp through rpaths into the conda env. To confirm:

otool -l "$QK/xtb-6.4.1/bin/xtb" | grep -A2 LC_RPATH   # should list .../envs/qupkake/lib

3d. Verify with a real optimization

A --version check or a single-point energy is not enough: the miscompiled -O3 binary passes both. Run a full geometry optimization, in gas phase and in ALPB water (the settings QupKake uses):

cd "$(mktemp -d)"
python -c "
from rdkit import Chem; from rdkit.Chem import AllChem
m = Chem.AddHs(Chem.MolFromSmiles('CC(=O)Oc1ccccc1C(=O)O')); AllChem.EmbedMolecule(m, randomSeed=42)
AllChem.MMFFOptimizeMolecule(m); Chem.MolToMolFile(m, 'aspirin.mol')"
export OMP_STACKSIZE=4G
"$QK/xtb-6.4.1/bin/xtb" aspirin.mol --opt > gas.txt 2>&1
"$QK/xtb-6.4.1/bin/xtb" aspirin.mol --opt --alpb water --lmo -P 4 > water.txt 2>&1
grep -E 'OPTIMIZATION CONVERGED|FAILED TO CONVERGE|TOTAL ENERGY|abnormal' gas.txt water.txt

Expected result: both runs print GEOMETRY OPTIMIZATION CONVERGED after about 20–25 iterations, with total energies of about −39.632 Eh (gas) and −39.652 Eh (water). A broken build gives FAILED TO CONVERGE or abnormal termination instead, and its energy rises at every optimization step.

Shell gotcha when scripting tests: in zsh (the macOS default shell), an unquoted $args holding "--opt --alpb water" does not word-split. xtb receives one garbage argument, quietly ignores it and runs a gas-phase single point. Use bash, or ${=args} in zsh.

4. Point QupKake at xtb

QupKake reads XTBPATH once, at import time (qupkake/__init__.py), and it must be the path to the executable. Set it, along with OpenMP settings, in conda activation hooks so conda activate qupkake sets everything:

mkdir -p "$CONDA_PREFIX/etc/conda/activate.d" "$CONDA_PREFIX/etc/conda/deactivate.d"

cat > "$CONDA_PREFIX/etc/conda/activate.d/xtb.sh" <<EOF
export _QUPKAKE_OLD_OMP_STACKSIZE="\${OMP_STACKSIZE:-}"
export _QUPKAKE_OLD_OMP_NUM_THREADS="\${OMP_NUM_THREADS:-}"
export XTBPATH=$QK/xtb-6.4.1/bin/xtb
export XTBHOME=$QK/xtb-6.4.1/share/xtb
export OMP_STACKSIZE=4G
export OMP_NUM_THREADS=1
EOF

cat > "$CONDA_PREFIX/etc/conda/deactivate.d/xtb.sh" <<'EOF'
unset XTBPATH XTBHOME
if [ -n "$_QUPKAKE_OLD_OMP_STACKSIZE" ]; then export OMP_STACKSIZE="$_QUPKAKE_OLD_OMP_STACKSIZE"; else unset OMP_STACKSIZE; fi
if [ -n "$_QUPKAKE_OLD_OMP_NUM_THREADS" ]; then export OMP_NUM_THREADS="$_QUPKAKE_OLD_OMP_NUM_THREADS"; else unset OMP_NUM_THREADS; fi
unset _QUPKAKE_OLD_OMP_STACKSIZE _QUPKAKE_OLD_OMP_NUM_THREADS
EOF

conda deactivate && conda activate qupkake && echo "$XTBPATH"

xtb 6.4.1 has its GFN2 parameters compiled in, so XTBHOME is optional. It’s set anyway for completeness.

5. Jupyter kernel

Jupyter kernels don’t run conda activation hooks, so put the variables in the kernel spec. The kernel’s PATH must also include conda: on import QupKake runs conda list xtb, and if conda is missing the import crashes with a confusing NameError: name 'XTB_LOCATION' is not defined.

python -m ipykernel install --user --name qupkake --display-name "Python (qupkake)"
python - <<EOF
import json, os
p = os.path.expanduser("~/Library/Jupyter/kernels/qupkake/kernel.json")
k = json.load(open(p))
k["env"] = {
    "XTBPATH": "$QK/xtb-6.4.1/bin/xtb",
    "XTBHOME": "$QK/xtb-6.4.1/share/xtb",
    "OMP_STACKSIZE": "4G",
    "OMP_NUM_THREADS": "1",
    "PATH": "$CONDA_PREFIX/bin:$(dirname "$CONDA_EXE")/../condabin:/usr/bin:/bin:/usr/sbin:/sbin",
}
json.dump(k, open(p, "w"), indent=1)
EOF

A worked example notebook is in QupKake_usage.ipynb. It covers a single SMILES, a CSV batch, tautomers, site drawings and calling QupKake from Python.

6. Verify QupKake end to end

cd "$(mktemp -d)"
python -c "import qupkake; print(qupkake.XTB_LOCATION)"     # source-built xtb path
qupkake smiles "CC(=O)Oc1ccccc1C(=O)O" -n aspirin -r test
grep -A1 "<pka" test/output/qupkake_output.sdf

Expected result: it runs in about 6 s and predicts an acidic pKa of about 3.58 for aspirin’s carboxylic acid (experimental ≈ 3.5), plus a basic site at about 1.5.

Two harmless messages appear on every import:

• CondaValueError: No packages match 'xtb' — QupKake checking for a conda xtb, which is deliberately absent.
• GPU available but not used — QupKake runs on CPU.

Pitfalls

-O3 miscompilation (the main trap)

With gfortran 16 at -O3, xtb 6.4.1 runs, prints the right version and gives correct single-point energies. In geometry optimizations, though, the energy rises at every step until the SCF fails. QupKake reports this as:

RuntimeError: Problem running xtb on molecule
...
TypeError: 'NoneType' object is not iterable      (in predict_sites)

and the input SDF in <root>/raw/ is truncated to 0 bytes. -O2 and -O0 give identical, correct results. The cause wasn’t investigated further: a gfortran 16 optimizer bug or undefined behaviour in the 2021 xtb code.

-mp N oversubscribes the CPU

QupKake passes its worker count to every xtb call as -P N. -mp 16 therefore launches 16 xtb jobs with 16 threads each. For batch runs, point XTBPATH at a small wrapper that forces one thread per xtb process:

#!/bin/bash
# xtb1: force "-P 1" whatever QupKake asks for
args=(); skip=0
for a in "$@"; do
  if [ $skip = 1 ]; then args+=("1"); skip=0; continue; fi
  [ "$a" = "-P" ] && skip=1
  args+=("$a")
done
exec "$QK/xtb-6.4.1/bin/xtb" "${args[@]}"

Use it with XTBPATH=/path/to/xtb1 qupkake file mols.csv -mp 22 .... Without -mp, QupKake still passes -P <all cores>, so the wrapper helps there too.

Input structures

• Give 2D SDFs as SMILES (CSV), not SDF. For SDF input QupKake passes the coordinates to xtb unchanged, with no hydrogens added. For CSV input it adds hydrogens and embeds in 3D first.
• One unembeddable molecule kills the whole file. RDKit’s Bad Conformer Id aborts the run. Pre-check with qupkake.cli.embed_molecule and remove failures.
• Output pKa values are sometimes the string tensor(4.21). Strip that before converting to float.
• Tautomer and charge state matter. QupKake predicts micro-pKa of the structure as drawn, so a different tautomer gives different values. On a 4.6k-compound set:
  • RDKit’s canonical tautomer (rdMolStandardize.TautomerEnumerator.Canonicalize) changed the most acidic pKa by >1 unit for 18% of the 1,900 molecules it altered.
  • The registered structures agreed better with QupKake’s own xtb-based tautomer search (-t) than the RDKit canonical ones did.

-t (tautomer search)

• -t runs one xtb optimization per RDKit-enumerated tautomer: a median of about 12 per molecule, up to 700+ for some. At about 4 CPU-seconds each, budget accordingly.
• If xtb fails on any tautomer, QupKake aborts the molecule (KeyError: 'totalenergy'). A drop-in launcher that skips failed tautomers is in runs/ForPKA_taut/qupkake_t.py (use python qupkake_t.py file ... -t in place of qupkake file ... -t).
• The GFN2-xTB/ALPB energy ranking sometimes favours implausible imidic-acid (C(OH)=N) and enol tautomers. Check its picks by eye.

Batch runs

For thousands of molecules, run many small single-worker QupKake jobs through xargs -P <cores>, not one big -mp run. QupKake splits a file into equal-count blocks per worker, so a single slow molecule stalls the whole run. Example scripts are in runs/ (ForPKA*/run_all.sh, run_job.sh, combine.py). Throughput on 22 cores is about 4.5 h per 4,600 drug-like molecules without -t.

Rebuilding xtb

conda activate qupkake
cd "$QK/xtb-src/xtb-6.4.1"
ninja -C build-o2 install

Then re-run the check in 3d.

I’ve uploaded details to GitHub https://github.com/drc007/qupkake-apple-silicon/tree/main/qupkake-apple-silicon

Citation

Abarbanel, O. D.; Hutchison, G. R. QupKake: Integrating Machine Learning and Quantum Chemistry for micro-pKa Predictions. J. Chem. Theory Comput. 2024. doi:10.1021/acs.jctc.4c00328

xtb: Bannwarth, C.; Ehlert, S.; Grimme, S. J. Chem. Theory Comput. 2019, 15, 1652. github.com/grimme-lab/xtb

His website also describes several other useful scripts http://www.tsurumi.yokohama-cu.ac.jp/pdl/PDL/PYTHON/index.html

Details here https://www.ukri.org/opportunity/frontier-ai-sme-champions-phase-1/

I thought it might be interesting to compare the performance of various clustering approaches on Apple Silicon (Apple M2 Ultra Mac with 192 GB of RAM.).

These four Jupyter notebooks cluster a set of molecules by structural similarity, using RDKit Morgan fingerprints (radius 3, 2048 bits) and Tanimoto similarity.

All four notebooks have the same inputs, outputs and layout, so you can run a dataset through any of them and compare the results. BitBIRCH and the two Butina notebooks group molecules by a similarity threshold. k-means instead uses a number of clusters you choose. All notebook use the same conda enviroment. Passing a square matrix to Butina.ClusterData needs a recent RDKit, so environment.yml requires RDKit 2026.03 or newer.

Choosing a threshold

Diverse sets, such as approved drugs or random screening compounds, give mostly single-molecule clusters at 0.65. Around 0.35–0.45 gives useful families. Some reference points:

Performance and memory

Timings are from an Apple M2 Ultra Mac with 192 GB of RAM.

All four notebooks were run on the same two sets: 149,396 random ZINC compounds and the 4,629-molecule ERK set. The clustering time is the clustering step alone. The notebook time includes loading the file (about 50 s for 150k), fingerprints and export. “Agreement” is the adjusted Rand index against RDKit Butina at the same threshold (1 = identical partition, 0 = no better than chance).

150k molecules:

An optional t-SNE map adds about 4 minutes at 150k, and about 20 s for the ERK set, in every notebook.

ERK set (4,629 molecules): every method clusters in 0.2–0.6 s, apart from k-means (2 s at k = 41, 14 s at k = 250, plus 23 s for the k scan). Peak RAM is under 1.2 GB. The agreement with RDKit Butina is 0.95 for MLX Butina at both thresholds, and 0.60 (at 0.65) and 0.34 (at 0.35) for BitBIRCH.

What this shows:

• MLX Butina gives Butina results at BitBIRCH cost. At 150k it is about 6× faster than RDKit Butina and uses about 20× less memory, because it keeps only the neighbour pairs instead of the full distance matrix. Its GPU memory peak was only 0.16 GB.
• MLX and RDKit Butina differ only in the greedy step. The neighbour lists are identical: every pair above the threshold matches RDKit’s Tanimoto exactly. mlxmolkit recounts neighbours after each cluster (RDKit’s reordering=True) and breaks ties the other way. RDKit’s notebook defaults to reordering=False. The difference is small at 0.65 (0.94–0.95 agreement, the same largest clusters). It grows at 0.35 on 150k (0.76, 21,362 against 23,064 clusters), where recounting lets later clusters absorb more molecules.
• BitBIRCH is fast but partitions differently. It makes more, smaller clusters than Butina at the same threshold, and agrees less with it at low thresholds.
• k-means isn’t suited to diverse sets. MiniBatchKMeans on 150k diverse molecules produces one huge catch-all cluster (19k molecules at k = 103, 44k at k = 250) alongside many small ones. Use it to split a set into a fixed number of groups, not to find families.

MLX Butina check. Before clustering, the MLX notebook checks that a Metal GPU is available. It estimates memory from the neighbour density in a sample, which matters most at low thresholds on sets of close analogues. It also times two trial runs to give a rough, usually high, estimate of the run time.

Butina memory check. Butina needs every pairwise distance, and RDKit holds them in a full n × n matrix. The notebook builds this matrix as float32 in blocks, which uses about 5× less memory than RDKit’s usual flat distance list. Before clustering starts, a check runs that:

• estimates the peak memory, including the matrix and neighbour lists
• compares it with the RAM available at that moment
• prints the largest set size that would fit, and an estimate of the run time
• stops with an error if the estimate is more than MEMORY_FRACTION of the available RAM