Automated calculation of pKa values

pka.py

import sys

pKa_ref = 11.4
ref = "heliotridane"

filename = sys.argv[1]
file = open(filename, "r+")
energies = {}
log_files = {}
names = []

for line in file:
    words = line.split()
    log_file  = words[0]
    ion = words[0].split("_")[0]
    name = ion[:-1]
    energy = float(words[1])
    if name not in names:
        names.append(name)
    if ion not in energies:
       energies[ion] = energy
       log_files[ion] = log_file
    elif energy < energies[ion]:
       energies[ion] = energy
       log_files[ion] = log_file

delta_G_ref = delta_G = energies[ref+"+"] - energies[ref+"-"]
for name in names:
    delta_G = energies[name+"+"] - energies[name+"-"]
    pKa = pKa_ref + (delta_G_ref - delta_G)/1.36
    print name, pKa, log_files[name+"+"], log_files[name+"-"]

pm6_opt.header

 $contrl runtyp=optimize icharg=0 $end
 $basis gbasis=pm6-d3h+ $end
 $statpt nstep=500 $end

pm6_sol.header

 $contrl icharg=0 $end
 $basis gbasis=PM6-D3H+ $end
 $pcm solvnt=WATER mxts=15000 smd=.T.  $end
 $tescav mthall=4 ntsall=60 $end

readv2.py

"""
type "python readv2.py *.log" to use
This program finds the last heat of formation in all GAMESS log files
for semiempirical geometry optimizations
"""

import sys

for filename in sys.argv[1:]:

    for line in reversed(open(filename,'r').readlines()):
        words = line.split() # split line into words
        if len(words) < 1:
            continue
        if words[0] == 'HEAT':
            heat = words[-2]
            print filename, heat
            break

smiles.list

heliotridane-a C[C@H]1CC[N@]2CCC[C@@H]12
heliotridane-b C[C@H]1CC[N@@]2CCC[C@@H]12
heliotridane+a C[C@H]1CC[N@H+]2CCC[C@@H]12
heliotridane+b C[C@H]1CC[N@@H+]2CCC[C@@H]12
comp1-a C[C@H]1[C@@H]2[C@@H](C(N(C)C2=O)=O)C3[N@]1CCC3
comp1-b C[C@H]1[C@@H]2[C@@H](C(N(C)C2=O)=O)C3[N@@]1CCC3
comp1+a C[C@H]1[C@@H]2[C@@H](C(N(C)C2=O)=O)C3[N@H+]1CCC3
comp1+b C[C@H]1[C@@H]2[C@@H](C(N(C)C2=O)=O)C3[N@@H+]1CCC3

smiles2coord.py

import subprocess,sys

filename = sys.argv[1]
file = open(filename, "r+")
babel = '/opt/openbabel/openbabel-master/bin/babel'

for line in file:
    words = line.split()
    name = words[0]
    smiles = words[1]
    smiles = smiles.replace("#","%23")
    smiles = smiles.replace("[","%5B")
    smiles = smiles.replace("]","%5D")
    sdffile = name+".sdf"
    xyzfile = name+".xyz"
    print sdffile,xyzfile
    url="https://cactus.nci.nih.gov/chemical/structure/"+smiles+"/sdf"
    subprocess.call(['curl',url, '-o',sdffile])
    subprocess.call([babel, '-isdf', sdffile, '-oxyz', xyzfile])