Hi Dr. Baker,

I responded too hastily to the last email, particularly regarding some of the workflow questions.  I've found the coulomb utility you  referred to.  Adding the fields would be useful, although because I'm using 0-charge dummy atoms to map the potential from grid points to the atom bond vector, potentials would be more useful.  Looking at the code, this should be rather trivial for me to add myself.

As far as step 3, I still could use some clarification.  Namely, after taking the difference of the .dx files, is there a utility that will map the reaction field .dx file to the atom coordinates in the same manner as the 'write atompot' keyword does in an APBS run?  I've tried reading the pot .dx file in a mg-dummy run, but the output file has 0 charge on all atom positions.

The more I think about your proposed work flow, the more I'm interested in giving it a try.  It appears to take care of many of the quirks that have been troubling me about my current setup.

All the best,
Andrew Ritchie
The University of Texas at Austin
Department of Chemistry and Biochemistry
The Webb Group

On Fri, May 3, 2013 at 4:02 PM, andrew ritchie <drew.w.ritchie@gmail.com> wrote:
Hi Dr. Baker,

Thanks for the reply. 

When changing the box size, I'm increasing glen.  We have not been changing fglen or cglen.  I'm copy/pasting an example of one of our input files at the bottom of this email.    In the behavior I described, I changed the dime such that glen/(dime-1) remained the same.  We're using mg-manual.  

We have been using cubic splines for the initial mapping. 

As far as your proposed work flow, I would just like to make sure I understand how exactly this is done.  
Steps 1 and 2 are perfectly clear.  
For step 3, I would take the difference .dx map, load that with my .pqr file in a mg-dummy calculation with the 'write atompot flat' keywords to map this back to my dummy atoms which I differentiate as before, correct?  
Step 4, when you refer to the analytic coulomb potential, are you referring to pot_prot or an all-atom E~Q/r^2 calculation?  Sorry if that sounds like a dumb question, but I'm not able to find any utilities listed on the user guide site called 'coulomb'.

Once I have clarification on step 4, I'll give this workflow a shot.  

Thanks for your time,
Andrew



----------------------------------------------------------------------------------------------------------------

read 

    mol pqr example.pqr 

end 

 

# stage 4 

# based on suggestions from Gernot Kierseritzky and stage4 described 

# in stages-logic.txt; changed to lpbe for boundary and convergence 

# issues when we go to finer grids, and spl2 because of fewer off- 

# the-grid warnings. 

elec 

    mg-manual 

 

    dime 97 97 97 

    glen 240.0 240.0 240.0 # 2.5 Angstroms 

 

    gcent mol 1 

    mol 1 

    lpbe 

    bcfl sdh # faster than mdh and probably good enough if the coarse grid is large enough 

 

    ion charge 1 conc 0.150 radius 2.0 

    ion charge -1 conc 0.150 radius 2.0 

    pdie 2.0 

    sdie 78.0 

    chgm spl2 

    srfm mol 

    srad 1.4 

    sdens 10.0 

    temp 300.0 

    calcenergy total 

    calcforce no 

end 

# stage2 

#      193      193      193 

#   10.000   10.000   10.000 

elec 

    mg-manual 

 

    glen   10.000   10.000   10.000 

    dime      193      193      193 

    gcent   85.560   58.480   39.830 

    mol 1 

    lpbe 

    bcfl focus 

 

    ion charge 1 conc 0.150 radius 2.0 

    ion charge -1 conc 0.150 radius 2.0 

    pdie 2.0 

    sdie 78.0 

    chgm spl2 

    srfm mol 

    srad 1.4 

    sdens 10.0 

    temp 300.0 

    calcenergy total 

    calcforce no 

    write atompot flat example_10_193 

end 

quit




On Fri, May 3, 2013 at 3:09 PM, Baker, Nathan <Nathan.Baker@pnnl.gov> wrote:
Hi Andrew --

Thanks for your email.  I'll try answer your questions in order and will add a few questions of my own.

Boundary condition.  When you increase the box size, which box are you increasing:  cglen or fglen?  Are you increase dime at the same time?  Is this an mg-manual or mg-auto calculation?  Our mg-auto focusing procedure is designed to gradually focus from the coarsest to finest grid, reducing the grid spacing by a controlled amount (https://sourceforge.net/p/apbs/code/ci/master/tree/src/generic/vhal.h#l393) at each step.  You might try changing this VREDFRAC parameter to see how it influences your results.

Charge interpolation.  We use trilinear interpolation to map the potential back onto the atoms.  However, if you used cubic splines to discretize the charges originally, you may get inconsistent results.  My suggestion would be to calculate your field somewhat differently:

 1.  Perform your current calculation on the system with distinct dielectric values for the interior (pdie eps1) and exterior (sdie eps2) of the protein.  Let's call the potential obtained from this calculation pot_prot.
 2.  Perform exactly the same calculation on the system with identical dielectric values for the interior (pdie eps1) and exterior (sdie eps1) of the protein.  Let's call the potential obtained from this calculation pot_ref.
 3.  Calculate the reaction field potential from the difference of the two maps (pot_prot – pot_ref).  You can use the APBS tool dxmath for this.  Differentiate this map as you describe below to obtain the reaction field.
 4.  Add back in the analytic Coulomb potential or field, scaled by the internal dielectric eps1.  You can use the APBS tool coulomb for this potential (we can add the field calculation, too, if it's useful).

The goal of this process is to use the high-accuracy analytic results for the Coulombic field wherever possible (particularly near charge centers) and to use the numerical lower-accuracy results for the reaction field only where absolutely necessary to account for the effects of the dielectric boundary.

Hope this helps,

--
Nathan Baker
Pacific Northwest National Laboratory
Phone: +1-509-375-3997
http://nabaker.me

From: Andrew Ritchie <drew.w.ritchie@gmail.com<mailto:drew.w.ritchie@gmail.com>>
Date: Thursday, May 2, 2013 9:31 AM
To: APBS-USERS mailing list <apbs-users@lists.sourceforge.net<mailto:apbs-users@lists.sourceforge.net>>
Subject: [apbs-users] Questions about grid spacing, boundary conditions, and write atompot

Greetings,

I'm using APBS to calculate electric fields within proteins.  I am doing this is using a two-stage focusing strategy with the multigrid solver.  I create dummy atoms along the vector which I want to know the electric field on and use the 'write atompot' keyword to return the potentials at those dummy atoms.  I then take the negative of the change in potential to find the field.  There are a couple interesting behaviors I have observed and I was wondering if this has been addressed/observed before and quantitatively explained. It should be noted that, although my end-goal is the electric field, I am talking about the potential values here.  I do not see a systematic change in the potentials in such a way as, although the potentials are different, their rate of change are the same.  I am seeing different potentials which is leading to different fields.

My first observation regards the boundary condition.  For the same protein structure with a fixed grid spacing, increasing the box size incrementally, I sometimes see relatively large changes in the potential (from ~7 kbT/e to ~5 kbT/e, so about a 30% change) while other times I see essentially no changes in the potential at my dummy atoms.  Situations in which I see large changes seem to be correlated to a number of large partially charged atoms (typically oxygen) going from outside the second-stage calculation box to inside the second-stage calculation box.  For snapshots I have done this test one, once the second-stage box exceeds ~20 angstroms, I stop seeing these large changes in potential.  This leads me to believe that within some r-distance, the boundary condition inadequately describes the surrounding system.  Is this an artifact of using a very coarse grid spacing for the first stage calculation, (first-stage calculation is is 240x240x240 angstrom box with 97 grid points in each dimension, so a grid spacing of 2.5 angstroms) or does the boundary condition simply not adequately account for large charges outside the box while the box is small enough for 1/r to be significant?

The second observation regards the interpolation of charge at the grid points to atoms (specifically my dummy atoms).  In my system, the dummy atoms are along a nitrile bond vector.  In the previously described observation, I zeroed out the charge on the nitrile atoms to better observe the contribution of the field due to everything else (since the potentials due to the nitrile part of the system dominate due to the small 1/r term).  Here I am zeroing out the charge on everything that is NOT the nitrile so that I can better observe just the nitrile contribution (although, again, the nitrile contributions dominate so zeroing nothing out gives essentially the same results).  What I've seen here are, for a fixed second-stage box size (using a 19 angstrom box), large changes in the potential near the nitrile atoms (changes of ~40%, 20-400 kbT/e within ~.2 angstroms of an atom) and smaller changes ( <0-3 kbT/e when > .2 angstroms from the bond vector).  I do not see this behavior when the charges on the nitrile are zeroed out, meaning this is only present when there are charged atoms near the dummy atoms.  The solution to this seems to be to cluster the dummy atoms closer to the bond mid point, although I've been unable to find the algorithm which maps the potential on the grid back to atoms, so I'm not sure what the resolution limit (if any) relative to the grid spacing is.

Thanks for the assistance,
Andrew Ritchie
Department of Chemistry and Biochemistry
The University of Texas at Austin
The Webb Group



--
-Andrew Ritchie



--
-Andrew Ritchie