search

Electrostatics / apbs-pdb2pqr

APBS - software for biomolecular electrostatics and solvation
http://www.poissonboltzmann.org/
127 stars 62 forks source link

pdb2pka results are platform dependent #377

Closed kmonson closed 8 years ago

kmonson commented 9 years ago

there is order dependent behavior in pdb2pka. The order that python dictionaries iterated over is platform dependent. This is affecting calculated background, desolvation, and interaction energies which in turn affects the resulting titration curves. Graphcut in NOT affected by this.

pka values have been seen to move by up to 0.6 ph.

kmonson commented 9 years ago

I speculate that is has to do with the fact that we debump the protein if needed.

kmonson commented 9 years ago

My original thought that it has something to do with the order things are iterated over was not correct. There are differences between pc and osx but they do not come into play.

The real reason appears to be a subtle difference in float handling which is having a significant effect on the debumping results in rare cases.

If the debumper cannot find a dihedral angle that is not bumped it chooses the angle that scores the best. Score is determined by the number of bumps and amount of overlap of any bumped atoms. Lowest score wins. The order dihedral angles are checked is not affected by platform.

When trying to debump CYS A 17 on my PC this happens: Best score of 1.79 at angle -34.01. New conflict set: ['HO']

On my Mac this happens at the same place in the run: Best score of 1.79 at angle 295.99. New conflict set: ['HO']

This affects the results of subsequent apbs runs. Note the difference in the resulting desolvation values.

PC: Vpbe_ctor: Using max ion radius (2 A) for exclusion function Debye length: 7.92928 A Current memory usage: 69.364 MB total, 133.521 MB high water Total electrostatic energy = 2.778674936966E+000 kJ/mol

Writing potential to potential0.dx

Vpbe_ctor: Using max ion radius (2 A) for exclusion function Debye length: 7.92928 A Current memory usage: 69.364 MB total, 133.521 MB high water Total electrostatic energy = 2.778674936966E+000 kJ/mol

Writing potential to potential1.dx

Vpbe_ctor: Using max ion radius (2 A) for exclusion function Debye length: 7.92928 A Current memory usage: 69.364 MB total, 133.521 MB high water Total electrostatic energy = 2.778674936966E+000 kJ/mol Writing potential to potential2.dx No energy arrays to destroy. Destroying multigrid structures. Destroying 1 molecules --------> Calculating self energy for residue CYS 17 state CTR02t in the protein Parsing input file D:\workspaces\github\apbs-pdb2pqr\pdb2pqr\pdb2pka_output\workspace\working.in...

Preparing to run 3 PBE calculations.

Vpbe_ctor: Using max ion radius (2 A) for exclusion function Debye length: 7.92928 A Current memory usage: 71.389 MB total, 133.521 MB high water Total electrostatic energy = 3.978953249938E+000 kJ/mol

Writing potential to potential0.dx

Vpbe_ctor: Using max ion radius (2 A) for exclusion function Debye length: 7.92928 A Current memory usage: 71.389 MB total, 134.181 MB high water Total electrostatic energy = 3.978953249938E+000 kJ/mol

Writing potential to potential1.dx

Vpbe_ctor: Using max ion radius (2 A) for exclusion function Debye length: 7.92928 A Current memory usage: 71.389 MB total, 134.181 MB high water Total electrostatic energy = 3.978953249938E+000 kJ/mol Writing potential to potential2.dx No energy arrays to destroy. Destroying multigrid structures. Destroying 1 molecules Desolvation for CYS 17 in state CTR02t is 0.484

OSX: Vpbe_ctor: Using max ion radius (2 A) for exclusion function Debye length: 7.92928 A Current memory usage: 69.357 MB total, 133.521 MB high water Total electrostatic energy = 2.692006121749E+00 kJ/mol

Writing potential to potential0.dx

Vpbe_ctor: Using max ion radius (2 A) for exclusion function Debye length: 7.92928 A Current memory usage: 69.357 MB total, 133.521 MB high water Total electrostatic energy = 2.692006121749E+00 kJ/mol

Writing potential to potential1.dx

Vpbe_ctor: Using max ion radius (2 A) for exclusion function Debye length: 7.92928 A Current memory usage: 69.357 MB total, 133.521 MB high water Total electrostatic energy = 2.692006121749E+00 kJ/mol Writing potential to potential2.dx No energy arrays to destroy. Destroying multigrid structures. Destroying 1 molecules --------> Calculating self energy for residue CYS 17 state CTR02t in the protein Parsing input file /Users/d3y382/workspaces/apbs-pdb2pqr/pdb2pqr/pdb2pka_output/workspace/working.in...

Preparing to run 3 PBE calculations.

Vpbe_ctor: Using max ion radius (2 A) for exclusion function Debye length: 7.92928 A Current memory usage: 71.389 MB total, 133.521 MB high water Total electrostatic energy = 4.013413119747E+00 kJ/mol

Writing potential to potential0.dx

Vpbe_ctor: Using max ion radius (2 A) for exclusion function Debye length: 7.92928 A Current memory usage: 71.389 MB total, 134.185 MB high water Total electrostatic energy = 4.013413119747E+00 kJ/mol

Writing potential to potential1.dx

Vpbe_ctor: Using max ion radius (2 A) for exclusion function Debye length: 7.92928 A Current memory usage: 71.389 MB total, 134.185 MB high water Total electrostatic energy = 4.013413119747E+00 kJ/mol Writing potential to potential2.dx No energy arrays to destroy. Destroying multigrid structures. Destroying 1 molecules Desolvation for CYS 17 in state CTR02t is 0.533

The problem cascades through to other results in the run.

kmonson commented 9 years ago

Possible fix:

Stop debumping the protein during the pdb2pka run.

Pros: Consistent cross platform results.

Cons: Some residues will no longer have valid curves as they may be stuck in a bumped state.

Originally pdb2pka did not do any debumping. I'm going to do some test runs over the weekend to see if it makes the pc and mac results match up better.

kmonson commented 9 years ago

PC vs MAC debumping still in place. pc_vs_osx

sobolevnrm commented 9 years ago

I would not stop debumping the protein. Can you provide %g or %1.15E output for the numbers so we can determine how large the floating point differences are?

Thanks,

Nathan

On Friday, October 30, 2015, kmonson notifications@github.com wrote:

Possible fix:

Stop debumping the protein during the pdb2pka run.

Pros: Consistent cross platform results.

Cons: Some residues will no longer have valid curves as they may be stuck in a bumped state.

— Reply to this email directly or view it on GitHub https://github.com/Electrostatics/apbs-pdb2pqr/issues/377#issuecomment-152651950 .

sobolevnrm commented 9 years ago

Hi Kyle -

Can you provide more significant digits on the score so we can determine how large the floating point error is?

Thanks,

Nathan

On Friday, October 30, 2015, kmonson notifications@github.com wrote:

My original thought that it has something to do with the order things are iterated over was not correct. There are differences between pc and osx but they do not come into play.

The real reason appears to be a subtle difference in float handling which is having a significant effect on the debumping results in rare cases.

If the debumper cannot find a dihedral angle that is not bumped it chooses the angle that scores the best. Score is determined by the number of bumps and amount of overlap of any bumped atoms. Lowest score wins. The order dihedral angles are checked is not affected by platform.

When trying to debump CYS A 17 on my PC this happens: Best score of 1.79 at angle -34.01. New conflict set: ['HO']

On my Mac this happens at the same place in the run: Best score of 1.79 at angle 295.99. New conflict set: ['HO']

This affects the results of subsequent apbs runs. Note the difference in the resulting desolvation values.

PC: Vpbe_ctor: Using max ion radius (2 A) for exclusion function Debye length: 7.92928 A Current memory usage: 69.364 MB total, 133.521 MB high water Total electrostatic energy = 2.778674936966E+000 kJ/mol Writing potential to potential0.dx

Vpbe_ctor: Using max ion radius (2 A) for exclusion function Debye length: 7.92928 A Current memory usage: 69.364 MB total, 133.521 MB high water Total electrostatic energy = 2.778674936966E+000 kJ/mol Writing potential to potential1.dx

Vpbe_ctor: Using max ion radius (2 A) for exclusion function Debye length: 7.92928 A Current memory usage: 69.364 MB total, 133.521 MB high water Total electrostatic energy = 2.778674936966E+000 kJ/mol Writing potential to potential2.dx No energy arrays to destroy. Destroying multigrid structures. Destroying 1 molecules --------> Calculating self energy for residue CYS 17 state CTR02t in the protein Parsing input file D:\workspaces\github\apbs-pdb2pqr\pdb2pqr\pdb2pka_output\workspace\working.in... Preparing to run 3 PBE calculations.

Vpbe_ctor: Using max ion radius (2 A) for exclusion function Debye length: 7.92928 A Current memory usage: 71.389 MB total, 133.521 MB high water Total electrostatic energy = 3.978953249938E+000 kJ/mol Writing potential to potential0.dx

Vpbe_ctor: Using max ion radius (2 A) for exclusion function Debye length: 7.92928 A Current memory usage: 71.389 MB total, 134.181 MB high water Total electrostatic energy = 3.978953249938E+000 kJ/mol Writing potential to potential1.dx

Vpbe_ctor: Using max ion radius (2 A) for exclusion function Debye length: 7.92928 A Current memory usage: 71.389 MB total, 134.181 MB high water Total electrostatic energy = 3.978953249938E+000 kJ/mol Writing potential to potential2.dx No energy arrays to destroy. Destroying multigrid structures. Destroying 1 molecules Desolvation for CYS 17 in state CTR02t is 0.484

OSX: Vpbe_ctor: Using max ion radius (2 A) for exclusion function Debye length: 7.92928 A Current memory usage: 69.357 MB total, 133.521 MB high water Total electrostatic energy = 2.692006121749E+00 kJ/mol Writing potential to potential0.dx

Vpbe_ctor: Using max ion radius (2 A) for exclusion function Debye length: 7.92928 A Current memory usage: 69.357 MB total, 133.521 MB high water Total electrostatic energy = 2.692006121749E+00 kJ/mol Writing potential to potential1.dx

Vpbe_ctor: Using max ion radius (2 A) for exclusion function Debye length: 7.92928 A Current memory usage: 69.357 MB total, 133.521 MB high water Total electrostatic energy = 2.692006121749E+00 kJ/mol Writing potential to potential2.dx No energy arrays to destroy. Destroying multigrid structures. Destroying 1 molecules --------> Calculating self energy for residue CYS 17 state CTR02t in the protein Parsing input file /Users/d3y382/workspaces/apbs-pdb2pqr/pdb2pqr/pdb2pka_output/workspace/working.in... Preparing to run 3 PBE calculations.

Vpbe_ctor: Using max ion radius (2 A) for exclusion function Debye length: 7.92928 A Current memory usage: 71.389 MB total, 133.521 MB high water Total electrostatic energy = 4.013413119747E+00 kJ/mol Writing potential to potential0.dx

Vpbe_ctor: Using max ion radius (2 A) for exclusion function Debye length: 7.92928 A Current memory usage: 71.389 MB total, 134.185 MB high water Total electrostatic energy = 4.013413119747E+00 kJ/mol Writing potential to potential1.dx

Vpbe_ctor: Using max ion radius (2 A) for exclusion function Debye length: 7.92928 A Current memory usage: 71.389 MB total, 134.185 MB high water Total electrostatic energy = 4.013413119747E+00 kJ/mol Writing potential to potential2.dx No energy arrays to destroy. Destroying multigrid structures. Destroying 1 molecules Desolvation for CYS 17 in state CTR02t is 0.533

The problem cascades through to other results in the run.

— Reply to this email directly or view it on GitHub https://github.com/Electrostatics/apbs-pdb2pqr/issues/377#issuecomment-152649974 .

kmonson commented 9 years ago

Can do.

I used repr(value) as it gives "the shortest string that when converted to float will return the original value; in case of ambiguity, the last digit is rounded to the closest value." That way we can see any difference at all.

PC: Checking if we must debump any residues... CYS A 17 HO is too close to CYS A 17 OXT Starting to debump CYS A 17... Debumping cutoffs: 2.0 for heavy-heavy, 1.5 for hydrogen-heavy, and 1.0 for hydrogen-hydrogen. Using dihedral angle number 0 to debump the residue. Best score of 1.7914065591240647 at angle -34.01276582806889. New conflict set: ['HO']

OSX: Checking if we must debump any residues... CYS A 17 HO is too close to CYS A 17 OXT Starting to debump CYS A 17... Debumping cutoffs: 2.0 for heavy-heavy, 1.5 for hydrogen-heavy, and 1.0 for hydrogen-hydrogen. Using dihedral angle number 0 to debump the residue. Best score of 1.7914065591240578 at angle 295.98723417193116. New conflict set: ['HO']

The first 14 digits of the score agree.

This residue is debumped 5 times before this. The first 4 the angles agree completely. The one right before the results disagree slightly and looks like this:

PC: Checking if we must debump any residues... CYS A 17 HO is too close to GLY A 16 C CYS A 17 HO is too close to GLY A 16 O CYS A 17 HO is too close to CYS A 17 N Starting to debump CYS A 17... Debumping cutoffs: 2.0 for heavy-heavy, 1.5 for hydrogen-heavy, and 1.0 for hydrogen-hydrogen. Using dihedral angle number 0 to debump the residue. No conflicts found at angle 320.98723238411094 Debumping Successful!

OSX: Checking if we must debump any residues... CYS A 17 HO is too close to GLY A 16 C CYS A 17 HO is too close to GLY A 16 O CYS A 17 HO is too close to CYS A 17 N Starting to debump CYS A 17... Debumping cutoffs: 2.0 for heavy-heavy, 1.5 for hydrogen-heavy, and 1.0 for hydrogen-hydrogen. Using dihedral angle number 0 to debump the residue. No conflicts found at angle 320.987232384111 Debumping Successful!

Back when we improved the debumping code I made sure that there were no longer accumulated error problems. Most likely the difference is the different representations of the same float on each platform.

kmonson commented 9 years ago

To detemine if the results I was getting from 1etm was a fluke I ran the remaining test residues on OSX. Every one failed to a difference degree.

Right now the test fails if for any pH on any titration curve the value at that pH differs from the accepted value by more than Epsilon. Epsilon = 0.01.

1a1p: ERROR: test curve CTR_A_13 has 28 bad points with 1.0020132953 cumulative error. ERROR: test curve ASP_A_6 has 1 bad points with 0.048457910123 cumulative error.

2m0w: ERROR: test curve CTR_A_23 has 32 bad points with 1.82526771798 cumulative error. ERROR: test curve NTR_A_1 has 2 bad points with 0.0237504658597 cumulative error. ERROR: test curve TYR_A_11 has 19 bad points with 0.37689734059 cumulative error. ERROR: test curve LYS_A_23 has 65 bad points with 23.2337772006 cumulative error. ERROR: test curve ASP_A_1 has 30 bad points with 1.51741408547 cumulative error.

2m8f: ERROR: test curve ASP_A_24 has 43 bad points with 5.35997089817 cumulative error. ERROR: test curve ASP_A_9 has 6 bad points with 0.063632561504 cumulative error. ERROR: test curve CTR_A_24 has 45 bad points with 1.74502968042 cumulative error.

I can make up graphs of these if desired.

1etm for comparison: ERROR: test curve CTR_A_17 has 46 bad points with 5.23565460061 cumulative error. ERROR: test curve NTR_A_6 has 5 bad points with 0.0905552842186 cumulative error. ERROR: test curve GLU_A_7 has 2 bad points with 0.027364182049 cumulative error.

sobolevnrm commented 9 years ago

What do you mean by "remaining test residues"?

On Wed, Nov 4, 2015 at 12:43 PM, kmonson notifications@github.com wrote:

To detemine if the results I was getting from 1etm was a fluke I ran the remaining test residues on OSX. Every one failed to a difference degree.

Right now the test fails if for any pH on any titration curve the value at that pH differs from the accepted value by more than Epsilon. Epsilon = 0.01.

1a1p: ERROR: test curve CTR_A_13 has 28 bad points with 1.0020132953 cumulative error. ERROR: test curve ASP_A_6 has 1 bad points with 0.048457910123 cumulative error.

2m0w: ERROR: test curve CTR_A_23 has 32 bad points with 1.82526771798 cumulative error. ERROR: test curve NTR_A_1 has 2 bad points with 0.0237504658597 cumulative error. ERROR: test curve TYR_A_11 has 19 bad points with 0.37689734059 cumulative error. ERROR: test curve LYS_A_23 has 65 bad points with 23.2337772006 cumulative error. ERROR: test curve ASP_A_1 has 30 bad points with 1.51741408547 cumulative error.

2m8f: ERROR: test curve ASP_A_24 has 43 bad points with 5.35997089817 cumulative error. ERROR: test curve ASP_A_9 has 6 bad points with 0.063632561504 cumulative error. ERROR: test curve CTR_A_24 has 45 bad points with 1.74502968042 cumulative error.

I can make up graphs of these if desired.

1etm for comparison: ERROR: test curve CTR_A_17 has 46 bad points with 5.23565460061 cumulative error. ERROR: test curve NTR_A_6 has 5 bad points with 0.0905552842186 cumulative error. ERROR: test curve GLU_A_7 has 2 bad points with 0.027364182049 cumulative error.

— Reply to this email directly or view it on GitHub https://github.com/Electrostatics/apbs-pdb2pqr/issues/377#issuecomment-153857270 .

kmonson commented 9 years ago

Sorry, I meant test proteins.

The test suite run 1etm, 1a1p, 2m0wm and 2m8f and compares with previously run results. There were picked because they run quickly and cover the cases of no HIS, 1 HIS, and 2 HIS.

keith923 commented 9 years ago

Here's the results from my Mac, note that some are the same and others are not:

1a1p

ERROR: test curve CTR_A_13 has 28 bad points with 1.0020132953 cumulative error. ERROR: test curve ASP_A_6 has 1 bad points with 0.048457910123 cumulative error.

1etm

ERROR: test curve CTR_A_17 has 46 bad points with 5.23565460061 cumulative error. ERROR: test curve NTR_A_6 has 5 bad points with 0.0905552842186 cumulative error. ERROR: test curve GLU_A_7 has 2 bad points with 0.027364182049 cumulative error.

2m0w

ERROR: test curve CTR_A_23 has 38 bad points with 3.41658097439 cumulative error. ERROR: test curve NTR_A_1 has 29 bad points with 0.840778347842 cumulative error. ERROR: test curve TYR_A_11 has 9 bad points with 0.167473715115 cumulative error. ERROR: test curve LYS_A_23 has 66 bad points with 25.3832202564 cumulative error. ERROR: test curve ASP_A_1 has 20 bad points with 0.470800100885 cumulative error.

2m8f

ERROR: test curve ASP_A_24 has 43 bad points with 5.35997089817 cumulative error. ERROR: test curve ASP_A_9 has 6 bad points with 0.063632561504 cumulative error. ERROR: test curve CTR_A_24 has 45 bad points with 1.74502968042 cumulative error.

keith923 commented 9 years ago

And these are the results from running on constance.

1a1p

ERROR: test curve CTR_A_13 has 46 bad points with 7.30081274411 cumulative error. ERROR: test curve ASP_A_6 has 7 bad points with 0.182040046079 cumulative error. ERROR: test curve NTR_A_1 has 11 bad points with 0.148747489576 cumulative error.

1etm

ERROR: test curve NTR_A_6 has 30 bad points with 1.56128164154 cumulative error. ERROR: test curve CTR_A_17 has 48 bad points with 6.98104981128 cumulative error. ERROR: test curve GLU_A_7 has 9 bad points with 0.111403729853 cumulative error.

2m0w

ERROR: test curve CTR_A_23 has 20 bad points with 0.437979681589 cumulative error. ERROR: test curve NTR_A_1 has 6 bad points with 0.065629377701 cumulative error. ERROR: test curve LYS_A_23 has 42 bad points with 5.19827524462 cumulative error. ERROR: test curve ASP_A_1 has 35 bad points with 2.32554255331 cumulative error.

2m8f

ERROR: test curve HIS_A_19 has 2 bad points with 0.031265639929 cumulative error. ERROR: test curve ASP_A_24 has 86 bad points with 34.348151658 cumulative error. ERROR: test curve TYR_A_15 has 5 bad points with 0.083392058539 cumulative error. ERROR: test curve CTR_A_24 has 104 bad points with 61.9038870089 cumulative error. ERROR: test curve ASP_A_9 has 23 bad points with 0.603528007126 cumulative error.

kmonson commented 8 years ago

TL;DR

Where it couldn't debump a residue the debumping code was settling ties for best configuration with floating point error. This could results in dihedral angles differing dramatically after debumping.

PDB2PKA is very sensitive to different debumping results.

It seems that some sort of symmetry between two bumped residues that cannot be fixed appears to be a common occurrence in pdb2pka. This symmetry results in scoring ties for multiple angles tested in the debumping code. If desired I can look this more closely to be sure it is what is going on.

Comparing scores using an Epsilon value and resolving ties with "first in wins" resolves all problems with all 4 test cases. We should test more proteins.

There is no perfect solution to this problem but using an Espilon should fix it most of the time.

Extra Details

When the debumping code fails to fix a bumped residue the chosen angle is the one that scores the lowest (least amount of overlap with atoms outside the residue). The atoms to be considered for calculating the score come from an optimization data structure call "cells". The order Atoms come out of the cells data structure for consideration is somewhat random and is platform dependent.

Trying to fix the problem by changing the order that Atoms are processed in is fixed is a bad idea. Even if the values are added in the same order that is no guarantee of the same result. See http://christian-seiler.de/projekte/fpmath/ under the heading "Problems with internal precision" specifically item 1. Add to that the fact that we have no control over how python itself was compiled and it becomes a fools errand.

The result is that multiple floating point numbers are added together in a different order on different platforms when calculating the bump score. In cases of what should be an exact tie you get platform dependent results.

As mentioned in previous comments above the APBS results are sensitive to the exact configuration of dihedral angles resulting from debumping. Fixing the debumping code to choose the same angle in the case of a tie will fix the problem.

I fixed the problem by introducing an Epsilon value to determine if two angles tied.

The old code looked like this:

if score == 0:
    #Hurray we are done!
   #Do some stuff

# Set the best angle
elif score < bestscore:
    bestscore = score
    bestangle = newangle
    foundImprovement = True

The new code looks like this:

if score == 0:
    #Hurray we are done!
    #Do some stuff

# Set the best angle
elif score < bestscore:
    diff = abs(bestscore - score)
    #Don't update if it's effectively a tie
    if diff > EPSILON:
        bestscore = score
        bestangle = newangle
        foundImprovement = True

The new code only affects cases where we don't find an unbumped state. If the improvement to the score does not exceed Epsilon than the best score and best angle are not updated. First in wins. EPSILON is 0.0000001.

This solution is not perfect. There is still a possibility that a pair of angle scores could be within Epsilon of each other on one platform and not on another. I don't believe that this is resolvable. I haven't seen it happen yet.

I speculate that in the cases we are frequently seeing ties in the bump scoring is because of some sort of symmetry between the two bumped residues. The only evidence I have is that before the fix the resulting angles disagree by 30 degrees and the set of conflicting atoms is not changed for the cases I've looked at. Also, introducing an Epsilon seems to seems to be a reliable fix also indicates some sort of symmetry.

Further study of the pdb state files would be need to confirm this.

sobolevnrm commented 8 years ago

This looks like a reasonable fix but I noticed you reopened the issue. Are there further problems?

kmonson commented 8 years ago

Closing it was unintentional. I wanted to make sure that it works on a few other proteins first and that you were good with the solution before I closed it. I'll test it on a few more proteins today then close it.

kmonson commented 8 years ago

Looks good. Every almost every titration curve is indistinguishable based on the platform. Ones that are only just barely. Every platform would pass the automated tests.

Biggest difference I could find anywhere was 1a1p, windows vs Constance.

ARG A 11 PH 12.9 Win 0.641808385103 Constance 0.641799425436

CTR A 13 PH 6.3 Win 0.409334548484 Constance 0.409298500778

None of this differs enough to fail the automated test suite (values must differ by more than 0.01) and the user would have to select a VERY specific ph value to tease out any platform specific differences in the resulting PQR file. I'm calling this done.

Fixed in f4e1d7873b9fb13b9be07d080e76a6a31aa30036

sobolevnrm commented 8 years ago

Congratulations! On Nov 10, 2015 4:21 PM, "kmonson" notifications@github.com wrote:

Closed #377 https://github.com/Electrostatics/apbs-pdb2pqr/issues/377.

— Reply to this email directly or view it on GitHub https://github.com/Electrostatics/apbs-pdb2pqr/issues/377#event-460773353 .