Trying to recreate "by hand" Born solvation calculation - problem with solvent radii?

[lvotapka]

Describe the bug A clear and concise description of what the bug is.

I'm performing a Born solvation calculation with pen and paper for testing purposes of other software, and I'm trying to check the results of my calculations against APBS results. However, I'm getting a discrepancy between the two unless I use an interior dielectric that neglects the radius of the solvent ions.

To Reproduce Steps to reproduce the behavior:

Here is a Python script that will compute the Born solvation of a single chloride ion of charge -1 and radius 2.27 Angstroms:

COULOMBS_PER_E = 1.602177e-19
METERS_PER_ANGSTROM = 1e-10
PICOSECONDS_PER_SECOND = 1e12
AVOGADROS_NUMBER = 6.022e23

pi = 3.1415926
q1 = -1.0 * COULOMBS_PER_E
solute_radius = 2.27 * METERS_PER_ANGSTROM
#solvent_radius = 0.0 * METERS_PER_ANGSTROM
solvent_radius = 1.5 * METERS_PER_ANGSTROM
a = solute_radius + solvent_radius
eps_0 = 8.854188e-12 # C^2/Jm
eps_r = 78

deltaG_1_1 = (q1**2 / (8 * pi * eps_0 * a)) * (1.0/eps_r - 1.0) * AVOGADROS_NUMBER / 1000

print("deltaG_1_1:", deltaG_1_1, "kJ/mol") 

Here is an APBS script that should compute the solvation energy of the same system (no other ions, just a dielectric):

read
  mol pqr chloride.pqr
end

elec name solv
  mg-manual
  dime 129 129 129 
  glen 16 16 16 
  gcent 3 4 5 
  mol 1
  npbe
  bcfl sdh

  pdie 1
  sdie 78
  srfm mol
  chgm spl2
  sdens 10.0
  srad 1.5
  temp 298
  calcenergy total
end
elec name ref
  mg-manual
  dime 129 129 129 
  glen 16 16 16 
  gcent 3 4 5 
  mol 1
  npbe
  bcfl sdh

  pdie 1
  sdie 1
  srfm mol
  chgm spl2
  sdens 10.0
  srad 1.5
  temp 298
  calcenergy total
end
print elecEnergy solv - ref end
quit

chloride.pqr.zip

Expected behavior A clear and concise description of what you expected to happen.

The Python script gives a solvation energy of -181.8 kJ/mol

APBS gives a solvation energy of -302.6 kJ/mol

It's very possible that something is wrong with the Python script, but, incidentally, when I change the script's "solvent_radius" to 0, the values match:

With a solvent radius of 0, the Python script gives a solvation energy of -302.1 kJ

When I change the "srad" parameter in the APBS input script to a small value, say, "0.01", then I get essentially the same solvation energy:

APBS with a solvent radius (srad) of 0.01 Angstroms gives a solvation energy of -302.2 kJ/mol

But the Born solvation formula would predict a very different solvation energy. This makes me wonder if there's a bug in APBS the way it handles the dielectric boundary, and that perhaps the solvent radius is being neglected. I observe a similar problem for different arguments of the "srfm" parameter, by the way.

Screenshots If applicable, add screenshots to help explain your problem.

Desktop (please complete the following information):

• OS: [e.g. iOS]
• Browser [e.g. chrome, safari]
• Version [e.g. 22]

APBS version 3.0 Ubuntu Linux 18.04

Smartphone (please complete the following information):

• Device: [e.g. iPhone6]
• OS: [e.g. iOS8.1]
• Browser [e.g. stock browser, safari]
• Version [e.g. 22]

Additional context Add any other context about the problem here.

It's very possible that I'm making a silly mistake somewhere that is causing all of this. But if so I haven't found it yet...

[jbardhan]

Hi Lane,

First, thank you for the detailed description of what you've done to check the behavior and what you were expecting.

As I read your post, I think there may be a misunderstanding of the Born solvation model used by APBS and most other continuum dielectric models for biomolecule solvation.

The dielectric boundary of a solute is given by the molecular surface (also known as the solvent-excluded surface) or some approximation to it. The molecular (aka solvent-excluded surface) is defined as the points of "closest approach" that the surface of a solvent probe sphere can get to the union of solute atoms.

The Born expression that you are using in the Python script is a statement that you expect the dielectric boundary of the solute to be given by the solvent-accessible surface (SAS), which is defined as the points of closest approach that the center of the solvent probe sphere can get to the same union. An equivalent definition is that the SAS is defined by "the union of spheres defined by the atom positions, and each sphere's radius is given by the atom radius plus the solvent radius".

It's possible that I have misunderstood what you're seeking. Does this help?

Best regards, Jay

[jbardhan]

This link provides a graphic that may be helpful. https://www.cgl.ucsf.edu/chimerax/docs/user/commands/surface.html

[lvotapka]

Hi Jay,

Thanks for your answer.

Yes, what you are saying makes perfect sense.

I'm taking a look at the APBS manual entry about srfm: https://apbs.readthedocs.io/en/latest/using/input/old/elec/srfm.html#elecsrfm, and I do see that it's describing what you just mentioned about the definition of the dielectric boundary at the molecular surface (not the solvent/ion accessible surface), so I see now that yes the program is behaving as expected.

But I admit I wouldn't have expected the dielectric boundary to be defined by the surface of the solvent probe, but rather the center of the probe. (After all, the dipole is at the center of the solvent molecule.)

[jbardhan]

Glad to be helpful! I appreciate why this confounds expectations. One way it was rationalized to me was, essentially, the water is modeled as a structureless continuum of point dipoles, so anywhere the interior of a water molecule (modeled as a sphere) can go, you should have the solvent dielectric.

[lvotapka]

Interesting...

They must have decided that using the molecular surface (not the solvent accessible surface) as the dielectric boundary was a best practice.

Hypothetically, just curious, if one wanted to use the solvent accessible surface (a surface defined by solute atom radius + solvent atom radius) for the dielectric boundary in APBS, how is the easiest way to do that?

[jbardhan]

The easiest way I know how to do it is the following 1) add solvent probe radius to each atom's radius in the PQR 2) for the APBS parameter probe radius, set it to 0.0

There may be alternatives!

[lvotapka]

I see. Another way might be to use APBS to generate the dielectric field DX using one PQR with large radii, input that to second APBS calculation that uses a PQR with normal radii - that way one could capture the correct ion-accessible surface for each of the ions with different radii.

Thanks Jay!