pKa inconsistencies - ChemAxon Forum Archive

User 870ab5b546

04-06-2007 16:52:43

Hi,

I'm measuring the most acidic and basic pKa's of each atom in the following compound, which is already in the "canonical resonance form" (JChem 3.2.6):

Code:

<?xml version="1.0" ?>

<MDocument>

<MChemicalStruct>

<molecule molID="m1">

<atomArray

atomID="a1 a2 a3 a4 a5 a6 a7 a8 a9 a10 a11 a12 a13 a14 a15 a16 a17"

elementType="C C C C C C O C O C C C C C C O H"

formalCharge="0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0"

x2="5.050965750132235 3.717296409083975 3.717296409083975 5.050965750132235 6.384635091180492 6.384635091180492 5.050965750132235 7.7183142130085285 7.7183142130085285 9.051993334836563 9.051993334836563 7.162519405430373 8.496198527258409 9.829877649086445 11.16355677091448 10.599877649086444 10.764975441456599"

y2="4.382371934332493 3.6123549937734527 2.0723211126553718 1.3023041720963313 2.0723211126553718 3.6123549937734527 5.9223719343324905 4.382354993773452 1.3023211126553722 3.6123549937734527 2.072354993773452 5.958567862126545 6.728567862126544 5.958567862126545 6.728567862126544 4.624888740298509 8.21609363461171"

/>

<bondArray>

<bond atomRefs2="a1 a2" order="1" />

<bond atomRefs2="a1 a6" order="1" />

<bond atomRefs2="a1 a7" order="2" />

<bond atomRefs2="a2 a3" order="1" />

<bond atomRefs2="a3 a4" order="1" />

<bond atomRefs2="a4 a5" order="1" />

<bond atomRefs2="a5 a6" order="1" />

<bond atomRefs2="a5 a9" order="2" />

<bond atomRefs2="a6 a8" order="1" />

<bond atomRefs2="a6 a12" order="1" />

<bond atomRefs2="a8 a10" order="1" />

<bond atomRefs2="a10 a11" order="2" />

<bond atomRefs2="a12 a13" order="1" />

<bond atomRefs2="a13 a14" order="1" />

<bond atomRefs2="a14 a15" order="1" />

<bond atomRefs2="a14 a16" order="2" />

<bond atomRefs2="a15 a17" order="1" />

</bondArray>

</molecule>

</MChemicalStruct>

</MDocument>

When I submit this compound to the pKa calculator in Marvin 4.1.8, it correctly gives me a basic pKa of 19.2 for the anionic C atom. When I call JChem 3.2.6 for the acidic and basic pKas of the same molecule, however, it gives me a basic pKa of 29.98131480199701 for the anionic C atom. I can't figure out why I am getting the wrong answer from JChem. Can you help???

Here's the code. The first part is in this jsp page, and MolFunctions.pKapKbAtom() is in a Java file.

Code:

String substrate = request.getParameter("substrate");

Molecule origMolecule = null;

Molecule molecule = null;

String[] pKas;

String[] pKbs;

double[] molpKas = {0, 0, 0};

double[] molpKbs = {0, 0, 0};

int[] molpKaAtms = {0, 0, 0};

int[] molpKbAtms = {0, 0, 0};

int numAtoms = 0;

int origNumAtoms = 0;

if (substrate != null) {

origMolecule = MolImporter.importMol(substrate);

molecule = MolImporter.importMol(substrate);

origNumAtoms = molecule.getAtomCount();

molecule = MolFunctions.stripMetals(molecule);

numAtoms = molecule.getAtomCount();

pKaPlugin plugin = new pKaPlugin();

plugin.setMaxIons(8); // default 8

plugin.setBasicpKaLowerLimit(-15);

plugin.setAcidicpKaUpperLimit(50);

pKas = new String[numAtoms];

pKbs = new String[numAtoms];

molecule.aromatize(MoleculeGraph.AROM_BASIC);

plugin.setMolecule(molecule);

// new to JChem 3.2.1; makes where pKa and pKb are

// stored NOT dependent on charge of atom

plugin.setpKaPrefixType(pKaPlugin.DYNAMICpKaPREFIX);

plugin.setModel(pKaPlugin.MODEL_LARGE);

boolean runOK = plugin.run();

for (int atmIdx = 0; atmIdx < numAtoms; atmIdx++) {

pKas[atmIdx] = "";

pKbs[atmIdx] = "";

MolAtom atm = molecule.getAtom(atmIdx);

if (atm.getSymbol().equals("H") && atm.getCharge() == 0)

continue;

ArrayList<String> doneThat = new ArrayList<String>();

double[] bothpKas = MolFunctions.pKapKbAtom(molecule,

atmIdx, plugin, doneThat);

if (!Double.isNaN(bothpKas[0])) {

System.out.println("acidic pKa for atom " + (atmIdx + 1)

+ ": " + bothpKas[0]);

pKas[atmIdx] = Utils.truncateFloat((float) bothpKas[0], 1);

} else System.out.println("no acidic pKa for atom " + (atmIdx + 1));

if (!Double.isNaN(bothpKas[1])) {

System.out.println("basic pKa for atom " + (atmIdx + 1)

+ ": " + bothpKas[1]);

pKbs[atmIdx] = Utils.truncateFloat((float) bothpKas[1], 1);

} else System.out.println("no basic pKa for atom " + (atmIdx + 1));

} // for each atom in the molecule

plugin.getSortedValues(pKaPlugin.ACIDIC, molpKas, molpKaAtms);

plugin.getSortedValues(pKaPlugin.BASIC, molpKbs, molpKbAtms);

} else { // substrate is null

pKas = new String[0];

pKbs = new String[0];

} // if substrate is/is not null

//----------------------------------------------------------------------

// pKapKbAtom

//----------------------------------------------------------------------

// gets the smallest pKa and largest pKb of an atom in a molecule

// The plugin has already calculated the pKas of the molecule before being

// sent here.

// doneThat stores resonance structures already examined;

// prevents infinite loops in recursion.

public static double[] pKapKbAtom(Molecule molecule, int atmIdx,

pKaPlugin plugin, ArrayList<String> doneThat) {

println("Calculating pKa and pKb of atom " + (atmIdx + 1) + "...");

MolAtom atm = molecule.getAtom(atmIdx);

boolean bears_H = (atm.getImplicitHcount()

+ atm.getExplicitHcount() > 0);

// needs to bear H to be acidic

double[] pKaAcidic = (bears_H ?

plugin.getpKaValues(atmIdx, pKaPlugin.ACIDIC)

: null);

double[] pKaBasic =

plugin.getpKaValues(atmIdx, pKaPlugin.BASIC);

double[] retVal = new double[2];

retVal[0] = (pKaAcidic != null ? pKaAcidic[0] : Double.NaN);

retVal[1] = (pKaBasic != null ? pKaBasic[0] : Double.NaN);

// if atom is C and plugin has failed to give value, guess

if (atm.getSymbol().equals("C")) {

if (Double.isNaN(retVal[0]) && atm.getCharge() == 0

&& atm.getImplicitHcount()

+ atm.getExplicitHcount() > 0) {

retVal[0] = getCpKa(molecule, atmIdx);

println(" guessing at pKaAcidic of CH ... value is "

+ retVal[0]);

} // C is negative, pKaBasic not calculated

if (Double.isNaN(retVal[1]) && atm.getCharge() == -1) {

retVal[1] = getCpKa(molecule, atmIdx);

println(" guessing at pKaBasic of C(-) ... value is "

+ retVal[1]);

} // C is negative, pKaBasic not calculated

} // if it's a C atom

// see if there's a resonance form with a charge < 0 on this atom

ArrayList<Integer> beenThereIndices = new ArrayList<Integer>();

beenThereIndices.add(atmIdx); // prevents infinite loops in recursion

doneThat.add(molecule.toFormat("mrv")); // prevents infinite loops in recursion

Molecule resMol = resStructWithNegChg(molecule, atmIdx,

beenThereIndices, doneThat);

if (resMol != null) {

println("Found a resonance form with a negative charge on "

+ "this atom; recursively calling pKapKbAtom() for:\n"

+ resMol.toFormat("mrv"));

pKaPlugin newPlugin = new pKaPlugin();

newPlugin.setpKaPrefixType(pKaPlugin.DYNAMICpKaPREFIX);

newPlugin.setModel(pKaPlugin.MODEL_LARGE);

newPlugin.setMaxIons(8);

newPlugin.setBasicpKaLowerLimit(-15);

newPlugin.setAcidicpKaUpperLimit(50);

resMol.calcHybridization();

boolean runOK = false;

try {

newPlugin.setMolecule(resMol);

runOK = newPlugin.run();

double[] newpKapKb = pKapKbAtom(resMol, atmIdx, newPlugin, doneThat);

println("Resonance form gave new pKa = " + newpKapKb[0]

+ ", new pKb = " + newpKapKb[1]);

if (Double.isNaN(retVal[1]) || newpKapKb[1] > retVal[1]) {

println("Changing old pKb " + retVal[1]

+ " to new pKb " + newpKapKb[1]);

retVal[1] = newpKapKb[1];

} // if the resonance structure is more basic than the original

} catch (PluginException e) {

System.out.println("pKaPlugin exception for:\n"

+ resMol.toFormat("mrv"));

}

} // if there's a resonance form with a negative charge on this atom

println("For atom " + (atmIdx + 1) + ", returning pKa = "

+ retVal[0] + ", pKb = " + retVal[1]);

return retVal;

} // pKapKbAtom(Molecule, int, pKaPlugin, ArrayList<String>)

Here's the log output:

Code:

Calculating pKa and pKb of atom 1...

resStructWithNegChg: beenThere = {1}

For atom 1, returning pKa = NaN, pKb = NaN

stage 3, cpd C13H17O3, atom C1 has no acidic pKa, no basic pKa

Calculating pKa and pKb of atom 2...

resStructWithNegChg: beenThere = {2}

For atom 2, returning pKa = 18.46324186445971, pKb = NaN

stage 3, cpd C13H17O3, atom C2 has acidic pKa 18.46324186445971, no basic pKa

Calculating pKa and pKb of atom 3...

Entering getCpKa for C atom whose pKa is not calculated by JChem.

getCpKa: Hybridization is sp3

guessing at pKaAcidic of CH ... value is 51.0

resStructWithNegChg: beenThere = {3}

For atom 3, returning pKa = 51.0, pKb = NaN

stage 3, cpd C13H17O3, atom C3 has acidic pKa 51.0, no basic pKa

Calculating pKa and pKb of atom 4...

resStructWithNegChg: beenThere = {4}

For atom 4, returning pKa = 19.3351664642691, pKb = NaN

stage 3, cpd C13H17O3, atom C4 has acidic pKa 19.3351664642691, no basic pKa

Calculating pKa and pKb of atom 5...

resStructWithNegChg: beenThere = {5}

For atom 5, returning pKa = NaN, pKb = NaN

stage 3, cpd C13H17O3, atom C5 has no acidic pKa, no basic pKa

Calculating pKa and pKb of atom 6...

resStructWithNegChg: beenThere = {6}

For atom 6, returning pKa = NaN, pKb = NaN

stage 3, cpd C13H17O3, atom C6 has no acidic pKa, no basic pKa

Calculating pKa and pKb of atom 7...

resStructWithNegChg: beenThere = {7}

resStructWithNegChg: beenThere = {7, 1}

resStructWithNegChg: beenThere = {7, 1, 1}

For atom 7, returning pKa = NaN, pKb = -7.562260074973646

stage 3, cpd C13H17O3, atom O7 has no acidic pKa, basic pKa -7.562260074973646

Calculating pKa and pKb of atom 8...

Entering getCpKa for C atom whose pKa is not calculated by JChem.

getCpKa: Hybridization is sp3

getCpKa: allylic or propargylic

guessing at pKaAcidic of CH ... value is 41.0

resStructWithNegChg: beenThere = {8}

For atom 8, returning pKa = 41.0, pKb = NaN

stage 3, cpd C13H17O3, atom C8 has acidic pKa 41.0, no basic pKa

Calculating pKa and pKb of atom 9...

resStructWithNegChg: beenThere = {9}

resStructWithNegChg: beenThere = {9, 5}

resStructWithNegChg: beenThere = {9, 5, 5}

For atom 9, returning pKa = NaN, pKb = -8.333727472295461

stage 3, cpd C13H17O3, atom O9 has no acidic pKa, basic pKa -8.333727472295461

Calculating pKa and pKb of atom 10...

Entering getCpKa for C atom whose pKa is not calculated by JChem.

getCpKa: Hybridization is sp2

guessing at pKaAcidic of CH ... value is 43.0

resStructWithNegChg: beenThere = {10}

For atom 10, returning pKa = 43.0, pKb = NaN

stage 3, cpd C13H17O3, atom C10 has acidic pKa 43.0, no basic pKa

Calculating pKa and pKb of atom 11...

Entering getCpKa for C atom whose pKa is not calculated by JChem.

getCpKa: Hybridization is sp2

guessing at pKaAcidic of CH ... value is 43.0

resStructWithNegChg: beenThere = {11}

resStructWithNegChg: beenThere = {11, 10}

For atom 11, returning pKa = 43.0, pKb = NaN

stage 3, cpd C13H17O3, atom C11 has acidic pKa 43.0, no basic pKa

Calculating pKa and pKb of atom 12...

Entering getCpKa for C atom whose pKa is not calculated by JChem.

getCpKa: Hybridization is sp3

guessing at pKaAcidic of CH ... value is 51.0

resStructWithNegChg: beenThere = {12}

For atom 12, returning pKa = 51.0, pKb = NaN

stage 3, cpd C13H17O3, atom C12 has acidic pKa 51.0, no basic pKa

Calculating pKa and pKb of atom 13...

resStructWithNegChg: beenThere = {13}

For atom 13, returning pKa = NaN, pKb = 29.98131480199701

stage 3, cpd C13H17O3, atom C13 has no acidic pKa, basic pKa 29.98131480199701

Calculating pKa and pKb of atom 14...

resStructWithNegChg: beenThere = {14}

For atom 14, returning pKa = NaN, pKb = NaN

stage 3, cpd C13H17O3, atom C14 has no acidic pKa, no basic pKa

Calculating pKa and pKb of atom 15...

resStructWithNegChg: beenThere = {15}

For atom 15, returning pKa = 26.173321502035975, pKb = NaN

stage 3, cpd C13H17O3, atom C15 has acidic pKa 26.173321502035975, no basic pKa

Calculating pKa and pKb of atom 16...

resStructWithNegChg: beenThere = {16}

Found a resonance form:

<?xml version="1.0" ?>

<MDocument>

<MChemicalStruct>

<molecule molID="m1">

<atomArray

atomID="a1 a2 a3 a4 a5 a6 a7 a8 a9 a10 a11 a12 a13 a14 a15 a16 a17"

elementType="C C C C C C O C O C C C C C C O H"

formalCharge="0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0"

x2="5.050965750132235 3.717296409083975 3.717296409083975 5.050965750132235 6.384635091180492 6.384635091180492 5.050965750132235 7.7183142130085285 7.7183142130085285 9.051993334836563 9.051993334836563 7.162519405430373 8.496198527258409 9.829877649086445 11.16355677091448 10.599877649086444 10.764975441456599"

y2="4.382371934332493 3.6123549937734527 2.0723211126553718 1.3023041720963313 2.0723211126553718 3.6123549937734527 5.9223719343324905 4.382354993773452 1.3023211126553722 3.6123549937734527 2.072354993773452 5.958567862126545 6.728567862126544 5.958567862126545 6.728567862126544 4.624888740298509 8.21609363461171"

/>

<bondArray>

<bond atomRefs2="a1 a2" order="1" />

<bond atomRefs2="a1 a6" order="1" />

<bond atomRefs2="a1 a7" order="2" />

<bond atomRefs2="a2 a3" order="1" />

<bond atomRefs2="a3 a4" order="1" />

<bond atomRefs2="a4 a5" order="1" />

<bond atomRefs2="a5 a6" order="1" />

<bond atomRefs2="a5 a9" order="2" />

<bond atomRefs2="a6 a8" order="1" />

<bond atomRefs2="a6 a12" order="1" />

<bond atomRefs2="a8 a10" order="1" />

<bond atomRefs2="a10 a11" order="2" />

<bond atomRefs2="a12 a13" order="1" />

<bond atomRefs2="a13 a14" order="2" />

<bond atomRefs2="a14 a15" order="1" />

<bond atomRefs2="a14 a16" order="1" />

<bond atomRefs2="a15 a17" order="1" />

</bondArray>

</molecule>

</MChemicalStruct>

</MDocument>

Checking against 1 resonance structures whether resonance structure has already been found and checked

haven't already looked at this resonance structure

Found a resonance form with a negative charge on this atom; recursively calling pKapKbAtom() for:

<?xml version="1.0" ?>

<MDocument>

<MChemicalStruct>

<molecule molID="m1">

<atomArray

atomID="a1 a2 a3 a4 a5 a6 a7 a8 a9 a10 a11 a12 a13 a14 a15 a16 a17"

elementType="C C C C C C O C O C C C C C C O H"

formalCharge="0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0"

x2="5.050965750132235 3.717296409083975 3.717296409083975 5.050965750132235 6.384635091180492 6.384635091180492 5.050965750132235 7.7183142130085285 7.7183142130085285 9.051993334836563 9.051993334836563 7.162519405430373 8.496198527258409 9.829877649086445 11.16355677091448 10.599877649086444 10.764975441456599"

y2="4.382371934332493 3.6123549937734527 2.0723211126553718 1.3023041720963313 2.0723211126553718 3.6123549937734527 5.9223719343324905 4.382354993773452 1.3023211126553722 3.6123549937734527 2.072354993773452 5.958567862126545 6.728567862126544 5.958567862126545 6.728567862126544 4.624888740298509 8.21609363461171"

/>

<bondArray>

<bond atomRefs2="a1 a2" order="1" />

<bond atomRefs2="a1 a6" order="1" />

<bond atomRefs2="a1 a7" order="2" />

<bond atomRefs2="a2 a3" order="1" />

<bond atomRefs2="a3 a4" order="1" />

<bond atomRefs2="a4 a5" order="1" />

<bond atomRefs2="a5 a6" order="1" />

<bond atomRefs2="a5 a9" order="2" />

<bond atomRefs2="a6 a8" order="1" />

<bond atomRefs2="a6 a12" order="1" />

<bond atomRefs2="a8 a10" order="1" />

<bond atomRefs2="a10 a11" order="2" />

<bond atomRefs2="a12 a13" order="1" />

<bond atomRefs2="a13 a14" order="2" />

<bond atomRefs2="a14 a15" order="1" />

<bond atomRefs2="a14 a16" order="1" />

<bond atomRefs2="a15 a17" order="1" />

</bondArray>

</molecule>

</MChemicalStruct>

</MDocument>

Calculating pKa and pKb of atom 16...

resStructWithNegChg: beenThere = {16}

For atom 16, returning pKa = NaN, pKb = 10.51728653351482

Resonance form gave new pKa = NaN, new pKb = 10.51728653351482

Changing old pKb -6.996512990109083 to new pKb 10.51728653351482

For atom 16, returning pKa = NaN, pKb = 10.51728653351482

stage 3, cpd C13H17O3, atom O16 has no acidic pKa, basic pKa 10.51728653351482

in stage 3, cpd C13H17O3 has molMostBasic = 29.98131480199701 and molMostAcidic = 18.46324186445971

User 851ac690a0

04-06-2007 18:44:23

Hi,

First of all your code is correct.

Two facts seem to be important.

1.

If the "micro pKa" calculation method used then the calculated pKa of the C(-) = 19.2…I think this was your personal setting in Marvin.

See.micropKa fig.

If the "macro pKa" calculation method used (this is the default) then the calculated pKa of the C(-) = 29.981… See.macropKa fig.

2.

Bug in my code.

The pKa=29.981... is assigned to the C(-) atom. The pKa=19.34 would be a better assignment. This is a bug in my "pKa - atom" assignment algorithm.

I will fix this bug now.

Jozsi

User 870ab5b546

04-06-2007 19:08:08

Is there a way to set the calculation method (micro or macro) from the API?

Of course, dynamic pKas depend on the microspecies, so setting the pKa prefix type to DYNAMICpKaPREFIX should automatically give microstate pKas.

Can you provide me with a JChem build that includes your bug fix? (A new Marvin build is not needed, as Marvin already works properly.)

User 851ac690a0

04-06-2007 19:27:43

Hi,

API functions of the pKaPlugin are given here:

http://www.chemaxon.com/marvin/doc/api/index.html

"micro - macro" calculation type can be adjusted with this function:

public void setMicropKaCalc(boolean micro)

Default: false (macro pKa calculation).

Parameter: micro - is true if micro pKa calculation is required.

Unfortunately I can not provide you with a bug fixed version.

Jozsi

User 870ab5b546

04-06-2007 19:49:22

And unfortunately, setting setMicropKaCalc(true) does not fix the bug. Sigh.

If you can't provide a build with the bugfix, do you have an estimate of when the next JChem version will be released?

User 851ac690a0

04-06-2007 22:32:05

Hi,

Quote:

...setting setMicropKaCalc(true) does not fix the bug

What is the pKa of the C(-) that you obtained with this setting?

Quote:

...when the next JChem version will be released

As far as I know we will make a release in mid. of June.

Jozsi

User 870ab5b546

04-06-2007 23:29:25

Jozsi wrote:

Quote:

...setting setMicropKaCalc(true) does not fix the bug

What is the pKa of the C(-) that you obtained with this setting?

Same as before: 29.98. Adding the line did not change the values returned by getpKaValues().

User 870ab5b546

06-06-2007 02:07:18

You will be relieved to hear, I was mistaken. Adding setMicropKaCalc(true) to the code does indeed fix the problem. I must have been looking at the wrong output or the value for a different atom.

Anyway, I don't need a new release, now that I have a workaround. You may want to add to your documentation for setpKaPrefixType() that when it is set to pKaPlugin.DYNAMICpKaPREFIX, then setMicropKaCalc(true) should also be specified.

User 851ac690a0

06-06-2007 06:34:13

Hi,

OK. Thank you.

Jozsi