QM 002: Comparing PCM and APBS

I compare the solvation energy of the Born ion using "quantum mechanics" and electrostatics. The system of interest is an $H^+$ ion, so the reason for the quotation marks around quantum mechanics is that there is no electron, so how can there be any quantum mechanics? Nevermind, I use the term because the PCM model implemented in the GAMESS program is compared with the solution to the Poisson-Boltzmann equation as can be found using the Adaptive Poisson-Boltzmann Solver (APBS) program.
The free energy of solution can be obtained from the equation

$$G_{solution} = E_{solute} + \frac{1}{2} \int_s d\vec{r} \sigma (\vec{r}) V(\vec{r})$$

First, the solution obtained from GAMESS/PCM is presented. The first term on the right side in the equation for the free energy in solution turns out to be exacltly zero because there is no electron in the system. Next, the second term on the right side is evaluated using the below PCM input for GAMESS.

1  $BASIS gbasis=STO ngauss=2 $end
2  $contrl mult=1 scftyp=rhf icharg=1 $end
3  $pcm solvnt=water $end
4  $DATA
5
6 C1
7  H 1.0 0.0 0.0 0.0
8  $END

At the end of the calculation, the following is printed to the log file:

 ----------------------------------------------
-------   RESULTS OF PCM CALCULATION   -------
----------------------------------------------

FREE ENERGY IN SOLVENT =                            =       -0.1813647781 A.U.
INTERNAL ENERGY IN SOLVENT =                        =        0.0000000000 A.U.
DELTA INTERNAL ENERGY =                             =        0.0000000000 A.U.
ELECTROSTATIC INTERACTION                           =       -0.1813647781 A.U.
PIEROTTI CAVITATION ENERGY                          =        0.0000000000 A.U.
DISPERSION FREE ENERGY                              =        0.0000000000 A.U.
REPULSION FREE ENERGY                               =        0.0000000000 A.U.
TOTAL INTERACTION (DELTA + ES + CAV + DISP + REP)   =       -0.1813647781 A.U.
TOTAL FREE ENERGY IN SOLVENT                        =       -0.1813647781 A.U.

The solvation energy (the energy released when transfering the ion from vacuum into the solvent) is now calculated as

$$(-0.1813 - 0.0)Ha*2526.5kJ/(mol*Ha)*\frac{1}{2} = -229.03kJ/mol.$$

Next, we us APBS to solve the problem. This requires a similar approach, involving the calculation of the gas phase and solution phase electrostatic energy. The below input for APBS carries out both calculations in the same run.

# READ IN MOLECULES
# Input copied from poissonboltzmann.org
read
mol pqr ion.pqr
end

# ...................................................................
# Electrostatics calculation on the solvated state
elec name solv        
mg-manual             # Specify the mode for APBS to run
dime 97 97 97         # The grid dimensions
#nlev 4               # Multigrid level parameter, not available in APBS 1.3
grid 0.33 0.33 0.33   # Grid spacing
gcent mol 1           # Center the grid on molecule 1
mol 1                 # Perform the calculation on molecule 1
lpbe                  # Solve the linearized Poisson-Boltzmann equation
bcfl mdh              # Use all multipole moments when calculating the potential
pdie 1.0              # Solute dielectric
sdie 78.54            # Solvent dielectric
chgm spl2             # Spline-based discretization of the delta functions
srfm mol              # Molecular surface definition
srad 1.4              # Solvent probe radius (for molecular surface)
swin 0.3              # Solvent surface spline window (not used here)
sdens 10.0            # Sphere density for accessibility object
temp 298.15           # Temperature
calcenergy total      # Calculate energies
calcforce no          # Do not calculate forces
end

# ...................................................................
# Calculate potential for reference (vacuum) state
elec name ref
mg-manual
dime 97 97 97
#nlev 4
grid 0.33 0.33 0.33
gcent mol 1
mol 1
lpbe
bcfl mdh
pdie 1.0
sdie 1.0
chgm spl2
srfm mol
srad 1.4
swin 0.3
sdens 10.0
temp 298.15
calcenergy total
calcforce no
end

# ...................................................................
# Calculate solvation energy
print elecEnergy solv - ref end
quit

The ion is defined by a PDB formatted input file:

ATOM      1  I   ION     1       0.000   0.000  0.000  1.00  3.0

The charge of the ion is defined in the second to last column, and the radius in the last one. The calculated energy of solvation is the energy of the solvated ion minus the energy of the gas phase ion:

print energy 1 (solv) - 2 (ref) end
Local net energy (PE 0) = -2.295878286327E+02 kJ/mol
Global net ELEC energy = -2.295878286327E+02 kJ/mol

As can be seen, the values from PCM and Poisson-Boltzmann equation are very close.

The problem can also be solved analytically. It involves the calculation of an integral over the surface of a sphere, which can be done as follows. First, the expression for the surface charge densitity $\sigma (r)$ needs to be evaluated. From the Poisson equation, it follows that $\sigma (\vec{r}) = -\frac{1}{4\pi}\left(1- \frac{1}{\epsilon}\right)Q\frac{1}{r^2}$ and using a potential $V(r) = Q\frac{1}{r}$, the integral to solve is thus

$$-\frac{1}{2}\int_s d\vec{r} \sigma (\vec{r}) V(\vec{r}) = -\frac{1}{8\pi} \left(1- \frac{1}{\epsilon}\right)Q^2\int_s d\vec{r} \frac{1}{r^3}$$

After realizing that the integration over the surface of the sphere requires the insertion of the Jacobi-Determinant $r^2 \sin \theta$, the new integration limits can be applied, the integral becomes independent of radius (setting $r \rightarrow R$) and one obtaines

$$-\frac{1}{8\pi}\left(1-\frac{1}{\epsilon}\right)Q^2\frac{1}{R^3}R^2\int_0^{2\pi} d\phi \int_0^{\pi} d\theta \sin \theta$$

which can be evaluated to

$$-\frac{1}{2}\left(1-\frac{1}{\epsilon}\right)\frac{Q^2}{R}.$$

From this it is seen that the solvation energy is depending more on the charge of the ion, than on its radius or the dielectric constant.