Hi Dr. Baker,
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----------------------------------------------------------------------------------------------------------------
mol pqr example.pqr
# 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.
dime 97 97 97
glen 240.0 240.0 240.0 # 2.5 Angstroms
gcent mol 1
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
# 193 193 193
# 10.000 10.000 10.000
glen 10.000 10.000 10.000
dime 193 193 193
gcent 85.560 58.480 39.830
ion charge 1 conc 0.150 radius 2.0
ion charge -1 conc 0.150 radius 2.0
write atompot flat example_10_193
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,
Pacific Northwest National Laboratory
From: Andrew Ritchie <firstname.lastname@example.org<mailto:email@example.com>>
Date: Thursday, May 2, 2013 9:31 AM
To: APBS-USERS mailing list <firstname.lastname@example.org<mailto:email@example.com>>
Subject: [apbs-users] Questions about grid spacing, boundary conditions, and write atompot
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,
Department of Chemistry and Biochemistry
The University of Texas at Austin
The Webb Group