From: David Minh (daveminh_at_gmail.com)
Date: Tue Jul 24 2012 - 10:36:23 CDT
Clayton Jarrat, a summer REU student, has resolved this problem. See below:
NAMD 2.9 had been calculating negative solvent accessible surface areas
for some molecules when when using the Generalized Born Implicit Solvent
algorithm. It is not immediately obvious when this happens because the
SASA is only used internally, as one of the values used in determining
the total electrostatic force (the "ELECT" value reported in the .out
file). I have done some research, picked apart the source code, and
found the problem.
NAMD uses the Linear Combination of Pairwise Overlaps method to
approximate the surface area of molecules. This method models each atom
as a hard sphere, and geometrically calculates the area of each sphere,
then how much of the area of any given sphere is buried inside another.
>From the paper in which this method was first proposed (Jo¨rg Weiser,
Peter S. Shenkin, W. Clark Still, Approximate Atomic Surfaces from
Linear Combinations of Pairwise Overlaps (LCPO), 1998), the formula used
for this is shown in the attached image. A rigid mathematical
derivation seems excessive here, but the meaning of each term is
simple. Ai is the SASA of an individual atom modeled as a hard sphere
i. The first term is simply the area of sphere i. The second term
calculates from the radii and distance between centers of neighbouring
spheres (possibly two spheres of arbitrary size j and k, for this
example), how much surfaces area of sphere i is buried inside its
neighbours, j and k. However, it is possible for spheres j and k to
overlap each other inside i, and in the second term, the area buried
inside both j and k is subtracted twice, which underestimates Ai. The
third term calculates how much area of i is subtracted twice in this
way, and corrects this. Lastly, it is easy to draw three spheres, i, j,
and k, such that j and k both overlap with i, and with each other, but
none of the overlap of j and k is inside i, regardless of how large the
overlap is between j and k. As the distance between i and its two
neighbours decreases, and also as the distance between j and k
decreases, the percentage of the overlap between j and k that is inside
i increases, and fourth term accounts for this. The constants P1
through P4 are constants which are specific to each possible number of
large neighbours for hybridization of each element.
NAMD uses exactly the formula from this paper, and has, at the end of
the source code file "ComputeLCPO.c", a structure which contains all of
the P1 - P4 values for it will need, which are copied exactly from the
appendix of the aforementioned paper. Obviously, for every atom type,
the area of the hard sphere model buried inside neighbors can not exceed
the total area of that sphere. However, the P1 and P2 values for sp3
hybridized N with one large neighbour are such that this is precisely
what NAMD calculates. The P1 value is an order of magnitude smaller
than the P2 value, which makes no physical sense. NAMD subtracts as
"buried" area more area than there is in total for sp3 N with one
neighbour by nearly a factor of 10.
Weiser, Shenkin, and Still calculate these P values three different ways
for each atom, and the numbers they show in the appendix for the very
same atom's P1 value that were arrived at using the other two methods
are entirely sensible. Each is larger than magnitude, but of opposite
sign, from the corresponding P2 value, and an order of magnitude larger
than this value. If one changes the exponent from -2 to -1, this value
will be very close to what was obtained using the other two methods, but
slightly smaller, which is consistent with P1 values obtained using this
method for every other atom in the table. It seems extremely likely
that this is simply a typo in the exponent. However, it may call the
whole table into question. I have tried to contact the authors of the
paper, but haven't heard back.
In any case, the surface area which NAMD calculates for any molecule
containing an sp3 N with one neighbour will be underestimated, and can
be negative if such N atoms make up a large percentage of the area of
On Jun 15, 2012, at 10:08 AM, David Minh wrote:
> Hi NAMD users,
> Has anyone else seen a negative surface area value?
> I'm using NAMD 2.9 and get a negative surface area for a host-guest system. I see it for the host, and for the host-guest complex, but not for the guest itself. An example is attached. Is this a quirk of the SASA algorithm or is it a bug?
> ETITLE: TS BOND ANGLE DIHED IMPRP ELECT VDW BOUNDARY MISC KINETIC TOTAL TEMP POTENTIAL TOTAL3 TEMPAVG
> Vacuum ENERGY: 0 0.8593 149.6652 34.6591 23.1692 76.6095 -65.2340 0.0000 0.0000 93.0849 312.8133 278.8241 219.7284 312.8224 278.8241
> GB ENERGY: 0 0.8593 149.6652 34.6591 23.1692 -70.1686 -65.2340 0.0000 0.0000 93.0849 166.0352 278.8241 72.9503 166.0377 278.8241
> GBSA ENERGY: 0 0.8593 149.6652 34.6591 23.1692 -77.1064 -65.2340 0.0000 0.0000 93.0849 159.0973 278.8241 66.0124 159.0998 278.8241
> Info: NAMD 2.9 for Linux-x86_64-multicore
> Info: Please visit http://www.ks.uiuc.edu/Research/namd/
> Info: for updates, documentation, and support information.
> Info: Please cite Phillips et al., J. Comp. Chem. 26:1781-1802 (2005)
> Info: in all publications reporting results obtained with NAMD.
> Info: Based on Charm++/Converse 60400 for multicore-linux64-iccstatic
> Info: Built Mon Apr 30 14:00:48 CDT 2012 by jim on naiad.ks.uiuc.edu
> Info: 1 NAMD 2.9 Linux-x86_64-multicore 8 dnb09.chem.duke.edu dm225
> Info: Running on 8 processors, 1 nodes, 1 physical nodes.
> Info: CPU topology information available.
> Info: Charm++/Converse parallel runtime startup completed at 0.00875401 s
> Info: 549.258 MB of memory in use based on /proc/self/stat
> Info: Configuration file is CUC_energy_vac.namd
> Info: Working in the current directory /xtmp/dm225/clusters/DSCR/Cu/1-build/M
This archive was generated by hypermail 2.1.6 : Tue Dec 31 2013 - 23:22:18 CST