Example: pKa prediction (PLAMS)¶
This example should be executed using PLAMS.
from scm.plams.interfaces.molecule.rdkit import from_smiles
import numpy as np
import multiprocessing
# In this example we compute pKa (acid dissociation constant) using MOPAC for a set of
# molecules. The molecules are defined using smiles strings, and are converted to xyz
# structures using the plams-rdkit interface.
# Important note: the predicted pKa strongly depend on the molecule's conformer.
# Here we use the lowest conformer predicted by rdkit's UFF.
# The difference between the values computed here and the results on the
# MOPAC website (ref_mopac_pKa) is due to different conformers
# Data taken from the online MOPAC manual: http://openmopac.net/manual/ (only a sub set)
data_tmp = [
# Molecule name smiles exp_pKa ref_mopac_pKa (from mopac's website)
["1-Naphthoic_acid", "C1=CC=C2C(=C1)C=CC=C2C(=O)O", 3.69, 4.35],
["2,2,2-Trichloroethanol", "C(C(Cl)(Cl)Cl)O", 12.02, 12.22],
["2,2,2-Trifluoroethanol", "C(C(F)(F)F)O", 12.40, 12.27],
["2,2-Dimethylpropionic_acid", "CC(C)(C)C(=O)O", 5.03, 5.23],
["2,3,4,6-Tetrachlorophenol", "C1=C(C(=C(C(=C1Cl)Cl)Cl)O)Cl", 7.10, 6.08],
["Acetic_acid", "CC(=O)O", 4.76, 5.00],
["Acrylic_acid", "C=CC(=O)O", 4.25, 4.65],
["Benzoid_acid", "C1=CC=C(C=C1)C(=O)O", 4.20, 4.30],
["Citric_acid", "C(C(=O)O)C(CC(=O)O)(C(=O)O)O", 3.13, 2.56],
["Ethanol", "CCO", 16.00, 16.37],
["Formic_acid", "C(=O)O", 3.77, 3.77],
["Glycine", "C(C(=O)O)N", 2.35, 2.53],
["Isoleucine", "CCC(C)C(C(=O)O)N", 2.32, 2.48],
["Methanol", "CO", 15.54, 15.23],
["o-Nitrophenol", "C1=CC=C(C(=C1)[N+](=O)[O-])O", 7.17, 7.52],
["Pentachlorophenol", "C1(=C(C(=C(C(=C1Cl)Cl)Cl)Cl)Cl)O", 4.90, 5.55],
["Phenol", "C1=CC=C(C=C1)O", 10.00, 9.71],
["Pyruvic_acid", "CC(=O)C(=O)O", 2.50, 2.85],
["T-Butanol", "CC(C)(C)O", 17.00, 16.25],
["Terephthalic_acid", "C1=CC(=CC=C1C(=O)O)C(=O)O", 3.51, 3.59],
["Valine", "CC(C)C(C(=O)O)N", 2.29, 2.61],
["Water", "O", 15.74, 15.75],
]
# Turn data_tmp into a dictionary:
systems = [{"name": d[0], "smiles": d[1], "exp_pKa": d[2], "ref_mopac_pKa": d[3]} for d in data_tmp]
# Create the molecules from the smiles using rdkit:
molecules = []
for system in systems:
# Compute 30 conformers, optimize with UFF and pick the lowest in energy.
mol = from_smiles(system["smiles"], nconfs=30, forcefield="uff")[0]
mol.properties.name = system["name"]
mol.properties.exp_pKa = system["exp_pKa"]
mol.properties.ref_mopac_pKa = system["ref_mopac_pKa"]
molecules.append(mol)
# MOPAC input:
s = Settings()
s.runscript.nproc = 1 # serial calculation
s.input.ams.Task = "GeometryOptimization"
s.input.ams.GeometryOptimization.Convergence.Step = 1.0e-3
s.input.ams.GeometryOptimization.Convergence.Gradients = 1.0e-5
s.input.mopac.model = "PM6"
s.input.mopac.properties.pKa = "Yes"
# Set up and run jobs:
jobs = MultiJob(children=[AMSJob(name=mol.properties.name, molecule=mol, settings=s) for mol in molecules])
jr = JobRunner(parallel=True, maxjobs=multiprocessing.cpu_count()) # run jobs in parallel
jobs.run(jobrunner=jr)
# Collect results:
for i, mol in enumerate(molecules):
pKaValues = jobs.children[i].results.readrkf("Properties", "pKaValues", file="mopac")
mol.properties.calc_pKa = np.mean(pKaValues) # If there is more than one pKa, take the average value
# Print results in a table:
print("Results:\n")
print("| {:28} | {:8} | {:8} | {:8} | {:8} |".format("Molecule", "exp pKa", "calc pKa", "ref", "calc-exp"))
for mol in molecules:
print(
"| {:28} | {:>8.2f} | {:>8.4f} | {:>8.2f} | {:>8.2f} |".format(
mol.properties.name,
mol.properties.exp_pKa,
mol.properties.calc_pKa,
mol.properties.ref_mopac_pKa,
mol.properties.calc_pKa - mol.properties.exp_pKa,
)
)
print("")
errors = [mol.properties.calc_pKa - mol.properties.exp_pKa for mol in molecules]
print("Mean signed error : {:4.2f}".format(np.mean(errors)))
print("Mean unsigned error: {:4.2f}".format(np.mean([abs(e) for e in errors])))
print("Root mean square error: {:4.2f}".format(np.sqrt(np.mean([e**2 for e in errors]))))
print("Done")