import os, math
from scm.plams import *
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: Be sure to add the path to your own AMSCRS directory here ##################
database_path = os.getcwd()
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()