MOPACでSMILESから酸解離定数pKaを計算する

SMILESで分子構造を入力してMOPACで酸解離定数pKaを計算します。

入力ファイル生成＆計算
計算結果

入力ファイル生成＆計算

必要なものは

Python
RDKit
MOPAC

以上です。Pythonを書きます。コードはこちらです。

from rdkit import Chem
from rdkit.Chem import AllChem

smiles = "CC(=O)O" #SMILES文字列
mol = Chem.MolFromSmiles(smiles)
mol = Chem.AddHs(mol)
AllChem.EmbedMolecule(mol)
AllChem.UFFOptimizeMolecule(mol)

with open("calc_pka.mop", "w") as f:
    f.write("PM6 PKA\n\n\n")
    for i in range(mol.GetNumAtoms()):
        atom = mol.GetAtomWithIdx(i)
        pos = mol.GetConformer().GetAtomPosition(i)
        f.write(f'{atom.GetSymbol():2} {pos.x:10.6f} {pos.y:10.6f} {pos.z:10.6f}\n')

smiles = の文字列を任意の化合物にしてください。このコードを実行してcalc_pka.mopファイルが作成されます。ファイルの中身はこちらです。

PM6 PKA


C   -0.971205  -0.073151  -0.214097
C    0.447528   0.283481   0.074439
O    0.793540   1.493515   0.131202
O    1.377703  -0.719664   0.327015
H   -1.395125   0.627051  -0.964537
H   -1.565508  -0.012197   0.721147
H   -1.028522  -1.106555  -0.615685
H    2.341589  -0.492479   0.540516

MOPACのpKa計算はPM6のみ対応。キーワードにPKAを書くだけです。これにより、分子を構造最適化してその構造でのPKAを計算してくれます。PKAキーワードの詳細はこちら。

PKA — MOPAC documentation

計算結果

計算は0.15秒で終わりました。速いですね。

計算後にcalc_pka.outファイルが出力されます。その中にpKaの結果が書かれています。

  :
  :
     GEOMETRY OPTIMISED USING EIGENVECTOR FOLLOWING (EF).     
     SCF FIELD WAS ACHIEVED                                   


                              PM6 CALCULATION
                                                       MOPAC v22.1.0 Windows
                                                       Fri Aug 29 23:23:20 2025




          FINAL HEAT OF FORMATION =       -101.22323 KCAL/MOL =    -423.51799 KJ/MOL


          COSMO AREA              =         91.89 SQUARE ANGSTROMS
          COSMO VOLUME            =         73.77 CUBIC ANGSTROMS

                          pKa for hydroxyl hydrogens

                      Sorted by pKa       Sorted by atom
                      Atom      pKa       Atom      pKa
                         8      4.960        8      4.960
 

          GRADIENT NORM           =          0.96029          =       0.33951 PER ATOM
          IONIZATION POTENTIAL    =         11.413263 EV
          HOMO LUMO ENERGIES (EV) =        -11.413  0.513
          NO. OF FILLED LEVELS    =         12
          MOLECULAR WEIGHT        =         60.0524         POINT GROUP:  C1  
  :
  :

このように、酢酸分子のOHの水素イオンが放出される際の酸解離定数は4.960と計算されました。実際は4.76なので、いい精度だと思います。