Automated pKa calculation¶

import os, math
from scm.plams import Settings, init, finish, CRSJob, config
import matplotlib.pyplot as plt


init()
# suppress plams output
config.log.stdout = 0


class PkaSystem:
    def __init__(
        self,
        comp_protonated: Settings = None,
        comp_deprotonated: Settings = None,
        systype="acid",
    ):
        self.comp_protonated = comp_protonated
        self.comp_deprotonated = comp_deprotonated
        self.type = systype.lower()
        self.g_cp = []
        self.g_cd = []
        self.pka = {}
        self.error = self.check_type()

    def check_type(self):
        if not all(isinstance(x, Settings) for x in [self.comp_protonated, self.comp_deprotonated]):
            print("Incorrect input type for compounds in PkaSystem... must be a plams.Settings instance")
            return False
        return True


class PkaCalculator:
    def __init__(
        self,
        systems=None,
        solv_protonated=None,
        solv_deprotonated=None,
        temp_range=[298.15, 298.15],
        steps=0,
        use_correction=True,
    ):

        self.systems = systems
        self.solv_protonated = solv_protonated
        self.solv_deprotonated = solv_deprotonated
        self.temp_range = temp_range
        self.steps = steps
        self.use_correction = use_correction
        self.R_const = 0.001987  # kcal/(mol K)

        self.g_sp = []
        self.g_sd = []

        if self.steps == 0:
            self.temp_vals = [temp_range[0]]
        else:
            self.temp_vals = [
                self.temp_range[0] * (1 - t_idx / self.steps) + self.temp_range[1] * (t_idx / self.steps)
                for t_idx in range(self.steps + 1)
            ]

        self.calc_G_values()
        self.calc_pkas()

    def calc_G_values(self):

        for temp in self.temp_vals:
            settings = Settings()
            settings.input.property._h = "ACTIVITYCOEF"
            # optionally, change to the COSMOSAC2013 method
            # settings.input.method = 'COSMOSAC2013'
            settings.input.temperature = temp

            compounds = [Settings()] * (2 + 2 * len(systems))
            compounds[0] = self.solv_deprotonated
            compounds[0].frac1 = 1.0
            compounds[1] = self.solv_protonated

            for i, system in enumerate(systems):
                compounds[2 + 2 * i] = system.comp_deprotonated
                compounds[2 + 2 * i + 1] = system.comp_protonated

            settings.input.compound = compounds
            my_job = CRSJob(settings=settings)
            out = my_job.run()
            res = out.get_results()

            g_vals = res["G solute"]
            for i, system in enumerate(systems):
                system.g_cd.append(g_vals[2 + 2 * i])
                system.g_cp.append(g_vals[2 + 2 * i + 1])
            self.g_sd.append(g_vals[0])
            self.g_sp.append(g_vals[1])

    def calc_pkas(self):
        for i, temp in enumerate(self.temp_vals):
            temp_key = round(temp, 3)
            for system in self.systems:
                g_diss = system.g_cd[i] - system.g_cp[i] + self.g_sp[i] - self.g_sd[i]
                if self.use_correction:
                    pka = self.calc_corrected_pka(g_diss, system, temp)
                else:
                    pka = g_diss / (self.R_const * temp * math.log(10)) - 1.74
                system.pka[temp_key] = pka

    def calc_corrected_pka(self, g_diss, system, temp):
        if system.type == "acid":
            return 0.62 * g_diss / (self.R_const * temp * math.log(10)) + 2.1
        else:
            return 0.67 * g_diss / (self.R_const * temp * math.log(10)) - 2.0


if __name__ == "__main__":

    ##################  Note: Ensure to download the coskf_pkas before running the script  ##################
    database_path = os.path.abspath("./coskf_pkas")

    if not os.path.exists(database_path):
        raise OSError(f"The provided path does not exist. Exiting.")

    systems = []

    benzoic_acid = Settings({"_h": os.path.join(database_path, "Benzoic_acid.coskf")})
    benzoic_acid_deprotonated = Settings({"_h": os.path.join(database_path, "conjugate_base_Benzoic_acid.coskf")})
    systems.append(PkaSystem(comp_protonated=benzoic_acid, comp_deprotonated=benzoic_acid_deprotonated, systype="acid"))

    pyridine = Settings({"_h": os.path.join(database_path, "Pyridine.coskf")})
    pyridineH = Settings({"_h": os.path.join(database_path, "conjugate_acid_Pyridine.coskf")})
    systems.append(PkaSystem(comp_protonated=pyridineH, comp_deprotonated=pyridine, systype="base"))

    water = Settings({"_h": os.path.join(database_path, "Water.coskf")})
    h3o = Settings({"_h": os.path.join(database_path, "conjugate_acid_Water.coskf")})

    PkaCalculator(systems, solv_protonated=h3o, solv_deprotonated=water, temp_range=[298.15, 348.15], steps=5)

    for system in systems:
        print(system.pka)

finish()