## Re: [Apbs-users] apbs-users Digest, Vol 42, Issue 4

 Re: [Apbs-users] apbs-users Digest, Vol 42, Issue 4 From: J.Dziedzic - 2009-12-14 10:36:00  Thank you very much, Gernot Kieseritzky for a detailed explanation of how to subtract the grid artefact out. This method works fine for the trivial dimer case I have presented earlier, yet I'm seeing weird results for the next case on my list -- a collection of three point charges that are the atomic cores of H2O. Of course I have no idea what the solvation energy is for such a system of isolated cores, yet I am worried that your approach seems to claim its total electrostatic energy with sdie=80, pdie=1 is negative. The input file is copied from your example: read mol pqr test.pqr end elec name test_inhom mg-auto dime 65 65 65 nlev 4 cglen 50 50 50 fglen 12 12 12 cgcent mol 1 fgcent mol 1 mol 1 lpbe bcfl mdh pdie 1.0 sdie 80.0 chgm spl2 srfm smol swin 0.3 srad 1.4 sdens 10.0 temp 298.15 calcenergy total calcforce no end elec name test_hom mg-auto dime 65 65 65 nlev 4 cglen 50 50 50 fglen 12 12 12 cgcent mol 1 fgcent mol 1 mol 1 lpbe bcfl mdh pdie 1.0 sdie 1.0 chgm spl2 srfm smol swin 0.3 srad 1.4 sdens 10.0 temp 298.15 calcenergy total calcforce no end print elecEnergy test_inhom - test_hom end quit ... with test.pqr now being REMARK 1 Cores of H2O ATOM 1 O XXX 1 12.000 12.957 12.000 6.00 1.50 ATOM 2 H XXX 1 12.000 12.000 12.000 1.00 1.20 ATOM 3 H XXX 1 12.927 13.197 12.000 1.00 1.20 The Coulomb energy with sdie=pdie=1 is E1=18334 kJ/mol as calculated by the coulomb.c utility. Calculation by hand gives a similar result of 18326 kJ/mol. Below I'm listing the \Delta{}G values produced by your script along with the value of the total electrostatic energy in sdie=1, pdie=\epsilon, obtained after \Delta{}G to the E1 value calculated above. \epsilon \Delta{}G E in dielectric of \epsilon 80 -27119 -8785 50 -26905 -8571 30 -26526 -8192 10 -24650 -6316 5 -21865 -3531 4 -20480 -2146 3 -18181 153 2 -13606 4728 1 0 18334 What is the reason for my getting negative energies? Could that be that the dielectric is polarizing sooo much? On a side note -- I thought of a simple approach to remove the grid artefact, can you comment on the feasibility? 1) Output the charge distribution to an OpenDX file. 2) Locate all 3x3x3 pockets of isolated charge floating in a sea of zero charge. 3) Calculate the Coulomb sum for the interaction of these 27 charges. Since they will always be in a cavity, epsilon will be that of pdiel and constant. Repeat for all 3x3x3 charge pockets. 4) Subtract the obtained result from what APBS produces. This assumes that all point charges are separated by at least one grid point of zero charge and that chgm spl2 is used. Is this feasible? thank you, - J. 

 Re: [Apbs-users] apbs-users Digest, Vol 42, Issue 4 From: J.Dziedzic - 2009-12-14 10:36:00  Thank you very much, Gernot Kieseritzky for a detailed explanation of how to subtract the grid artefact out. This method works fine for the trivial dimer case I have presented earlier, yet I'm seeing weird results for the next case on my list -- a collection of three point charges that are the atomic cores of H2O. Of course I have no idea what the solvation energy is for such a system of isolated cores, yet I am worried that your approach seems to claim its total electrostatic energy with sdie=80, pdie=1 is negative. The input file is copied from your example: read mol pqr test.pqr end elec name test_inhom mg-auto dime 65 65 65 nlev 4 cglen 50 50 50 fglen 12 12 12 cgcent mol 1 fgcent mol 1 mol 1 lpbe bcfl mdh pdie 1.0 sdie 80.0 chgm spl2 srfm smol swin 0.3 srad 1.4 sdens 10.0 temp 298.15 calcenergy total calcforce no end elec name test_hom mg-auto dime 65 65 65 nlev 4 cglen 50 50 50 fglen 12 12 12 cgcent mol 1 fgcent mol 1 mol 1 lpbe bcfl mdh pdie 1.0 sdie 1.0 chgm spl2 srfm smol swin 0.3 srad 1.4 sdens 10.0 temp 298.15 calcenergy total calcforce no end print elecEnergy test_inhom - test_hom end quit ... with test.pqr now being REMARK 1 Cores of H2O ATOM 1 O XXX 1 12.000 12.957 12.000 6.00 1.50 ATOM 2 H XXX 1 12.000 12.000 12.000 1.00 1.20 ATOM 3 H XXX 1 12.927 13.197 12.000 1.00 1.20 The Coulomb energy with sdie=pdie=1 is E1=18334 kJ/mol as calculated by the coulomb.c utility. Calculation by hand gives a similar result of 18326 kJ/mol. Below I'm listing the \Delta{}G values produced by your script along with the value of the total electrostatic energy in sdie=1, pdie=\epsilon, obtained after \Delta{}G to the E1 value calculated above. \epsilon \Delta{}G E in dielectric of \epsilon 80 -27119 -8785 50 -26905 -8571 30 -26526 -8192 10 -24650 -6316 5 -21865 -3531 4 -20480 -2146 3 -18181 153 2 -13606 4728 1 0 18334 What is the reason for my getting negative energies? Could that be that the dielectric is polarizing sooo much? On a side note -- I thought of a simple approach to remove the grid artefact, can you comment on the feasibility? 1) Output the charge distribution to an OpenDX file. 2) Locate all 3x3x3 pockets of isolated charge floating in a sea of zero charge. 3) Calculate the Coulomb sum for the interaction of these 27 charges. Since they will always be in a cavity, epsilon will be that of pdiel and constant. Repeat for all 3x3x3 charge pockets. 4) Subtract the obtained result from what APBS produces. This assumes that all point charges are separated by at least one grid point of zero charge and that chgm spl2 is used. Is this feasible? thank you, - J. 
 Re: [Apbs-users] apbs-users Digest, Vol 42, Issue 4 From: Gernot Kieseritzky - 2009-12-15 13:59:20 Hi! On Mon, 2009-12-14 at 10:35 +0000, J.Dziedzic wrote: > Below I'm listing the \Delta{}G values produced by your script > along with the value of the total electrostatic energy in sdie=1, > pdie=\epsilon, obtained after \Delta{}G to the E1 value calculated above. Shouldn't it be sdie=\epsilon and pdie=1? Otherwise \Delta{}G would be the solvation energy of set of point charges in a small \epsilon=80 bubble surrounded by vacuum. > \epsilon \Delta{}G E in dielectric of \epsilon > 80 -27119 -8785 > 50 -26905 -8571 > 30 -26526 -8192 > 10 -24650 -6316 > 5 -21865 -3531 > 4 -20480 -2146 > 3 -18181 153 > 2 -13606 4728 > 1 0 18334 > > What is the reason for my getting negative energies? Could that > be that the dielectric is polarizing sooo much? Well, it would suggest that the solvation energy in absolute values is larger than the Coulomb repulsion of your system. Sounds kind of counter intuitive with point charges of +6 and +1 at close distance. But I have the feeling you mixed up sdie and pdie (see above). > On a side note -- I thought of a simple approach to remove the grid > artefact, can you comment on the feasibility? > > 1) Output the charge distribution to an OpenDX file. > 2) Locate all 3x3x3 pockets of isolated charge floating in a sea of > zero charge. > 3) Calculate the Coulomb sum for the interaction of these 27 charges. > Since they will always be in a cavity, epsilon will be that of pdiel > and constant. Repeat for all 3x3x3 charge pockets. > 4) Subtract the obtained result from what APBS produces. > > This assumes that all point charges are separated by at least one > grid point of zero charge and that chgm spl2 is used. > > Is this feasible? Maybe. But in the end it's easier (and faster?) to simply subtract the homogenious APBS result. If electrostatics is really time-limiting and the grid artefact creates a serious problem in your application you should consider using a Generalized Born approach. The currently fastest implemention comes from the group of Prof. Caflisch: http://www3.interscience.wiley.com/journal/116327343/abstract?CRETRY=1&SRETRY=0 and is part of the latest CHARMM version 35. Best regards, Gernot Kieseritzky 
 Re: [Apbs-users] apbs-users Digest, Vol 42, Issue 4 From: J.Dziedzic - 2009-12-15 13:59:01 Gernot Kieseritzky wrote: > Hi! > > On Mon, 2009-12-14 at 10:35 +0000, J.Dziedzic wrote: >> Below I'm listing the \Delta{}G values produced by your script >> along with the value of the total electrostatic energy in sdie=1, >> pdie=\epsilon, obtained after \Delta{}G to the E1 value calculated above. > > Shouldn't it be sdie=\epsilon and pdie=1? Otherwise \Delta{}G would be > the solvation energy of set of point charges in a small \epsilon=80 > bubble surrounded by vacuum. Sorry, a typo on my part. I got it right in the input file, so the calculation was for pdie=1.0 and varying sdie, as it should be. >> \epsilon \Delta{}G E in dielectric of \epsilon >> 80 -27119 -8785 >> 50 -26905 -8571 >> 30 -26526 -8192 >> 10 -24650 -6316 >> 5 -21865 -3531 >> 4 -20480 -2146 >> 3 -18181 153 >> 2 -13606 4728 >> 1 0 18334 >> >> What is the reason for my getting negative energies? Could that >> be that the dielectric is polarizing sooo much? > > Well, it would suggest that the solvation energy in absolute values is > larger than the Coulomb repulsion of your system. Sounds kind of counter > intuitive with point charges of +6 and +1 at close distance. But I have > the feeling you mixed up sdie and pdie (see above). I mixed them up only when typing the email, the input file was fine. >> On a side note -- I thought of a simple approach to remove the grid >> artefact, can you comment on the feasibility? >> >> 1) Output the charge distribution to an OpenDX file. >> 2) Locate all 3x3x3 pockets of isolated charge floating in a sea of >> zero charge. >> 3) Calculate the Coulomb sum for the interaction of these 27 charges. >> Since they will always be in a cavity, epsilon will be that of pdiel >> and constant. Repeat for all 3x3x3 charge pockets. >> 4) Subtract the obtained result from what APBS produces. >> >> This assumes that all point charges are separated by at least one >> grid point of zero charge and that chgm spl2 is used. >> >> Is this feasible? > > Maybe. I have tried this approach. It doesn't really work, because there is a shift in the potential by a constant, I believe. APBS uses the q-phi integration and my approach uses a sum over pairs of q's. When phi is shifted, the results do not match. Is there any way to determine the constant term in the potential? > But in the end it's easier (and faster?) to simply subtract the > homogenious APBS result. I agree. I would have used the approach you have suggested, but I was worried about the negative results I got for the three point charges. So, do you think this is a valid result, just counterintuitive? I am also worried that with this technique the result will depend on the precise dimensions of the cavity. > If electrostatics is really time-limiting and > the grid artefact creates a serious problem in your application you > should consider using a Generalized Born approach. The currently fastest > implemention comes from the group of Prof. Caflisch: > > http://www3.interscience.wiley.com/journal/116327343/abstract?CRETRY=1&SRETRY=0 > > and is part of the latest CHARMM version 35. Thank you. The electrostatic calculation is not a time bottleneck in my case, far from it. It's just that I need to be able to get reliable values for electrostatic energy for a set of point charges embedded in an inhomogenous dielectric. I have another doubt as well. When a point charge is spilled onto the grid, the grid artefact occurs because all of these grid points are now in the potential generated by their own collection of charges. This artefact is meant to cancel out when calculating \Delta{}G, because a similar scenario occurs in the inhomogeneous dielectric case. However, as the boundary conditions are the same for both cases, the addition of the dielectric around the cavity now implies a subtle change in the potential in the cavity to accomodate the fact, that the same boundary conditions must be satisfied. Would this not change the grid artefact energy subtly, leading to the canceling-out being not 100% thanks, - Jacek Dziedzic 
 Re: [Apbs-users] apbs-users Digest, Vol 42, Issue 4 From: Gernot Kieseritzky - 2009-12-15 17:37:30 Hi! On Tue, 2009-12-15 at 13:58 +0000, J.Dziedzic wrote: > >> \epsilon \Delta{}G E in dielectric of \epsilon > >> 80 -27119 -8785 > >> 50 -26905 -8571 > >> 30 -26526 -8192 > >> 10 -24650 -6316 > >> 5 -21865 -3531 > >> 4 -20480 -2146 > >> 3 -18181 153 > >> 2 -13606 4728 > >> 1 0 18334 > >> > >> What is the reason for my getting negative energies? Could that > >> be that the dielectric is polarizing sooo much? > > > > Well, it would suggest that the solvation energy in absolute values is > > larger than the Coulomb repulsion of your system. Sounds kind of counter > > intuitive with point charges of +6 and +1 at close distance. But I have > > the feeling you mixed up sdie and pdie (see above). > > I mixed them up only when typing the email, the input file was fine. I repeated your calculation at sdie = 80 and I think the values are fine (whether they are physically reasonable is another story). The very low solvation energy of -27119 is due to the high charge density creating a strong reaction field. Setting the partial charges to that of TIP3 water from the CHARMM force field roughly reproduces its solvation energy of about 10 kcal/mol. So there's no grid artefact at work here. In your file 'test.pqr' you have assigned a charge +6 corresponding to the oxygen nucleus, and +1 corresponding to the protons. So the result somehow suggests that such an arrangement of naked nuclei would be stable in bulk water due to solvation. I guess, that's non-sense, because forces other than electrostatics would prevent this. In reality, of course, you have electrons around the nuclei -- so you would have to add the electrostatic energy due to the electrons to get the full picture. Best regards, Gernot Kieseritzky 
 Re: [Apbs-users] apbs-users Digest, Vol 42, Issue 4 From: J.Dziedzic - 2009-12-15 17:43:21 Gernot Kieseritzky wrote: > Hi! > [...] > In your file 'test.pqr' you have assigned a charge +6 corresponding to > the oxygen nucleus, and +1 corresponding to the protons. So the result > somehow suggests that such an arrangement of naked nuclei would be > stable in bulk water due to solvation. I guess, that's non-sense, > because forces other than electrostatics would prevent this. In reality, > of course, you have electrons around the nuclei -- so you would have to > add the electrostatic energy due to the electrons to get the full > picture. Yes, of course. The arrangement I've suggested represents the naked cores of the oxygen and hydrogen atoms, as seen from the point of view of a DFT-with-pseudopotential approach. The electrons are treated separately, they are represented by a continuous distribution of charge density, that integrates to a charge of -8, discretized on a grid. This part of my APBS calculation is progressing without doubts, that's why I have omitted it. thank you, - Jacek Dziedzic