A Posteriori Analysis and Efficient Refinement Strategies for the Poisson-Boltzmann Equation

TL;DR

This paper employs adjoint-based a posteriori analysis to quantify errors in solving the Poisson-Boltzmann equation and introduces efficient refinement strategies to improve the accuracy of solvation free energy calculations.

Contribution

It presents a novel application of a posteriori error analysis to the PBE and develops new refinement strategies for finite element solutions.

Findings

Accurate error quantification of solvation free energy

Identification of key error sources in PBE solutions

Development of effective adaptive refinement methods

Abstract

The Poisson-Boltzmann equation (PBE) models the electrostatic interactions of charged bodies such as molecules and proteins in an electrolyte solvent. The PBE is a challenging equation to solve numerically due to the presence of singularities, discontinuous coefficients and boundary conditions. Hence, there is often large error in the numerical solution of the PBE that needs to be quantified. In this work, we use adjoint based a posteriori analysis to accurately quantify the error in an important quantity of interest, the solvation free energy, for the finite element solution of the PBE. We identify various sources of error and propose novel refinement strategies based on a posteriori error estimates.

Tables

Table 1. Table 1: The discretizations to be refined based on the dominant source of error. The refinement may be uniform or adaptive.

Dominant source	Discretization to refine
$E^r$	Refine $𝒯$
$E^m$	Refine ${𝒯}_m$
$E^{\Gamma }$	Refine simplices containing $\Gamma \cap {\Omega }_m$
$E^{\partial \Omega }$	Refine simplices containing $\partial \Omega \cap \Omega$

Table 2. Table 2: Born ion: Error estimate, effectivity ratio and error contributions for the linearized PBE ( 16 ) with N ​ ( u ) = u 𝑁 𝑢 𝑢 N(u)=u . N 𝑁 N is the number of vertices in 𝒯 𝒯 \mathcal{T} . The terms E r superscript 𝐸 𝑟 E^{r} , E m superscript 𝐸 𝑚 E^{m} , E Γ superscript 𝐸 Γ E^{\Gamma} and E ∂ Ω superscript 𝐸 Ω E^{\partial\Omega} are defined in Theorem 1 .

$N$	Est. Err.	${\gamma }_{\mathrm{eff}}$	$E^r$	$E^m$	$E^{\Gamma }$	$E^{\partial \Omega }$
6718	-1.14	1.05	2.05e-01	5.26e-09	-1.34e+00	4.86e-04
52014	-0.248	1.04	1.19e-01	2.04e-06	-3.67e-01	1.36e-04

Table 3. Table 3: Born ion: Error estimate, effectivity ratio and error contributions for the nonlinear PBE ( 16 ) with N ​ ( u ) = sinh ⁡ ( u ) 𝑁 𝑢 𝑢 N(u)=\sinh(u) . The terms E r superscript 𝐸 𝑟 E^{r} , E m superscript 𝐸 𝑚 E^{m} , E Γ superscript 𝐸 Γ E^{\Gamma} and E ∂ Ω superscript 𝐸 Ω E^{\partial\Omega} are defined in Theorem 2 .

$N$	Est. Err.	${\gamma }_{\mathrm{eff}}$	$E^r$	$E^m$	$E^{\Gamma }$	$E^{\partial \Omega }$
6718	-1.16	1.05	1.88e-01	5.13e-09	-1.34e+00	4.24e-04
52014	-0.254	1.04	1.13e-01	2.69e-06	-3.67e-01	1.18e-04

Table 4. Table 4: Methanol: Error estimate, effectivity ratio and error contributions for the linearized PBE ( 16 ) with N ​ ( u ) = u 𝑁 𝑢 𝑢 N(u)=u .

$N$	Est. Err.	${\gamma }_{\mathrm{eff}}$	$E^r$	$E^m$	$E^{\Gamma }$	$E^{\partial \Omega }$
11769	-0.924	1.03	4.27e-01	4.67e-06	-1.35e+00	2.43e-06
90264	-0.231	1.01	1.49e-01	9.34e-06	-3.81e-01	6.87e-07

Table 5. Table 5: Methanol: Error estimate, effectivity ratio and error contributions for the nonlinear PBE ( 16 ) with N ​ ( u ) = sinh ⁡ ( u ) 𝑁 𝑢 𝑢 N(u)=\sinh(u) .

$N$	Est. Err.	${\gamma }_{\mathrm{eff}}$	$E^r$	$E^m$	$E^{\Gamma }$	$E^{\partial \Omega }$
11769	-0.924	1.02	4.27e-01	3.31e-06	-1.35e+00	2.41e-06
90264	-0.232	1.01	1.49e-01	8.77e-06	-3.81e-01	6.81e-07

Table 6. Table 6: Born ion: Error contributions and refinement for the linearized PBE.

It.	$N$	Error	$E^r$	$E^m$	$E^{\Gamma }$	$E^{\partial \Omega }$
0	6718	-1.14	2.05e-01	5.26e-09	-1.34e+00	4.86e-04
1	15541	-0.214	1.53e-01	9.75e-07	-3.68e-01	4.78e-04
2	41760	0.0124	1.03e-01	2.30e-07	-9.07e-02	4.78e-04
3 ($\Gamma$)	141855	0.0621	8.41e-02	1.38e-07	-2.25e-02	4.78e-04
3-Uniform	323084	0.00853	3.09e-02	1.29e-07	-2.25e-02	1.35e-04

Table 7. Table 7: Uniform Contribution Refinement strategy defined in § 4 applied to the linearized PBE for the Born ion. It. refers to the refinement iteration or level.

It.	$N$	Est. Err.	${\gamma }_{\mathrm{eff}}$	$E^r$	$E^m$	$E^{\Gamma }$	$E^{\partial \Omega }$
0	6718	-1.14	1.05	2.05e-01	5.26e-09	-1.34e+00	4.86e-04
1	15541	-0.214	1.06	1.53e-01	9.75e-07	-3.68e-01	4.78e-04
2	41760	0.0124	–	1.03e-01	2.30e-07	-9.07e-02	4.78e-04

Table 8. Table 8: Uniform Contribution Refinement strategy defined in § 4 applied to the nonlinear PBE for the Born ion.

It.	$N$	Est. Err.	${\gamma }_{\mathrm{eff}}$	$E^r$	$E^m$	$E^{\Gamma }$	$E^{\partial \Omega }$
0	6718	-1.16	1.05	1.88e-01	5.13e-09	-1.34e+00	4.29e-04
1	15541	-0.224	1.05	1.44e-01	5.23e-07	-3.68e-01	4.17e-04
2	41760	3.49e-3	-	9.38e-02	4.86e-07	-9.07e-02	4.16e-04

Table 9. Table 9: Uniform Contribution Refinement strategy defined in § 4 applied to the linearized PBE for the Methanol.

It.	$N$	Est. Err.	${\gamma }_{\mathrm{eff}}$	$E^r$	$E^m$	$E^{\Gamma }$	$E^{\partial \Omega }$
0	11769	-0.924	1.03	4.27e-01	4.67e-06	-1.35e+00	2.43e-06
1	20283	-0.275	1.02	1.06e-01	3.07e-07	-3.82e-01	2.40e-06
2	46212	-0.094	0.998	1.20e-03	7.72e-07	-9.52e-02	2.40e-06

Table 10. Table 10: Uniform Contribution Refinement strategy defined in § 4 applied to the nonlinear PBE for Methanol.

It.	$N$	Est. Err.	${\gamma }_{\mathrm{eff}}$	$E^r$	$E^m$	$E^{\Gamma }$	$E^{\partial \Omega }$
0	11769	-0.924	1.03	4.27e-01	3.31e-06	-1.35e+00	2.41e-06
1	20283	-0.276	1.02	1.06e-01	8.55e-07	-3.82e-01	2.38e-06
2	46212	-0.0944	0.996	8.65e-04	4.38e-07	-9.52e-02	2.38e-06

Table 11. Table 11: Adaptive Contribution Refinement strategy defined in § 4 applied to the linearized PBE for the Born ion.

It.	$N$	Est. Err.	${\gamma }_{\mathrm{eff}}$	$E^r$	$E^m$	$E^{\Gamma }$	$E^{\partial \Omega }$
0	6718	-1.14	1.05	2.05e-01	5.26e-09	-1.34e+00	4.86e-04
1	9476	-1.07	1.03	1.50e-01	1.07e-07	-1.22e+00	4.79e-04
2	14273	-0.521	1.02	1.24e-01	3.08e-07	-6.45e-01	4.78e-04
3	21293	-0.235	1.02	1.07e-01	3.88e-07	-3.42e-01	4.78e-04
4	33320	-0.0895	1.02	9.68e-02	2.72e-07	-1.87e-01	4.78e-04
5	60908	-0.00803	0.98	8.22e-02	2.61e-07	-9.07e-02	4.77e-04
6	112597	0.0112	1.13	5.90e-02	1.75e-07	-4.83e-02	4.77e-04
7	206897	0.0183	1.11	4.31e-02	2.87e-07	-2.53e-02	4.75e-04

Table 12. Table 12: Adaptive Contribution Refinement strategy defined in § 4 applied to the nonlinear PBE for the Born ion.

It.	$N$	Est. Err.	${\gamma }_{\mathrm{eff}}$	$E^r$	$E^m$	$E^{\Gamma }$	$E^{\partial \Omega }$
0	6718	-1.16	1.05	1.88e-01	5.13e-09	-1.34e+00	4.29e-04
1	9476	-1.08	1.02	1.38e-01	8.05e-08	-1.22e+00	4.19e-04
2	14273	-0.531	1.02	1.14e-01	4.24e-07	-6.45e-01	4.17e-04
3	21293	-0.244	1.02	9.77e-02	7.15e-08	-3.42e-01	4.16e-04
4	33320	-0.0983	1.01	8.80e-02	3.05e-07	-1.87e-01	4.16e-04
5	54112	-0.021	0.98	7.94e-02	-7.17e-08	-1.01e-01	4.16e-04
6	106120	-0.00159	-	4.96e-02	3.95e-08	-5.16e-02	4.13e-04
7	200323	0.0109	1.22	3.66e-02	1.31e-07	-2.61e-02	4.10e-04

Table 13. Table 13: Adaptive Contribution Refinement strategy defined in § 4 applied to the linearized PBE for Methanol.

It.	$N$	Est. Err.	${\gamma }_{\mathrm{eff}}$	$E^r$	$E^m$	$E^{\Gamma }$	$E^{\partial \Omega }$
0	11769	-0.924	1.03	4.27e-01	4.67e-06	-1.35e+00	2.43e-06
1	12380	-0.564	1.04	1.98e-01	1.55e-06	-7.62e-01	2.42e-06
2	13852	-0.391	1.03	8.72e-02	1.56e-06	-4.78e-01	2.41e-06
3	17353	-0.275	1.02	1.21e-02	-1.17e-06	-2.87e-01	2.40e-06
4	22732	-0.206	1.01	-2.94e-02	-1.81e-06	-1.77e-01	2.40e-06
5	33019	-0.156	0.99	-5.16e-02	-5.96e-07	-1.04e-01	2.40e-06
6	50784	-0.127	0.98	-6.92e-02	-1.08e-06	-5.80e-02	2.40e-06
7	86224	0.000197	-	3.48e-02	2.20e-07	-3.46e-02	2.40e-06

Table 14. Table 14: Adaptive Contribution Refinement strategy defined in § 4 applied to the nonlinear PBE for methanol.

It.	$N$	Est. Err.	${\gamma }_{\mathrm{eff}}$	$E^r$	$E^m$	$E^{\Gamma }$	$E^{\partial \Omega }$
0	11769	-0.924	1.03	4.27e-01	3.31e-06	-1.35e+00	2.41e-06
1	12380	-0.564	1.04	1.97e-01	8.87e-07	-7.62e-01	2.40e-06
2	13852	-0.392	1.03	8.68e-02	4.96e-06	-4.78e-01	2.39e-06
3	17353	-0.275	1.02	1.17e-02	-1.42e-06	-2.87e-01	2.38e-06
4	22732	-0.206	1.00	-2.97e-02	-1.47e-06	-1.77e-01	2.38e-06
5	33036	-0.156	0.99	-5.19e-02	-8.73e-07	-1.04e-01	2.38e-06
6	50796	-0.127	0.98	-6.95e-02	-1.01e-06	-5.80e-02	2.38e-06
7	86276	1.67e-05	-	3.46e-02	-8.72e-07	-3.45e-02	2.38e-06

Table 15. Table 15: Classical Refinement strategy defined in § 4 applied to the linearized PBE for the Born ion.

It.	$N$	Est. Err.	${\gamma }_{\mathrm{eff}}$	$E^r$	$E^m$	$E^{\Gamma }$	$E^{\partial \Omega }$
0	6718	-1.14	1.05	2.05e-01	5.26e-09	-1.34e+00	4.86e-04
1	9481	-1.07	1.03	1.50e-01	5.62e-08	-1.22e+00	4.79e-04
2	14299	-0.512	1.02	1.24e-01	-4.16e-08	-6.36e-01	4.78e-04
3	21638	-0.23	1.02	1.06e-01	2.82e-07	-3.37e-01	4.78e-04
4	33959	-0.0868	1.02	9.58e-02	1.83e-07	-1.83e-01	4.78e-04
5	55762	-0.0116	1.00	8.68e-02	9.50e-08	-9.89e-02	4.78e-04
6	95162	0.0275	1.03	8.10e-02	1.28e-08	-5.40e-02	4.78e-04
7	165202	0.0487	1.02	7.75e-02	9.00e-08	-2.94e-02	4.78e-04

Table 16. Table 16: Classical Refinement strategy defined in § 4 applied to the nonlinear PBE for the Born ion.

It.	$N$	Est. Err.	${\gamma }_{\mathrm{eff}}$	$E^r$	$E^m$	$E^{\Gamma }$	$E^{\partial \Omega }$
0	6718	-1.16	1.05	1.88e-01	5.13e-09	-1.34e+00	4.29e-04
1	9481	-1.08	1.02	1.38e-01	1.06e-07	-1.22e+00	4.19e-04
2	14299	-0.521	1.02	1.14e-01	-3.29e-08	-6.36e-01	4.17e-04
3	21638	-0.239	1.02	9.70e-02	1.84e-07	-3.37e-01	4.16e-04
4	33959	-0.0955	1.01	8.71e-02	1.99e-07	-1.83e-01	4.16e-04
5	55762	-0.02	0.98	7.85e-02	1.07e-07	-9.89e-02	4.16e-04
6	95162	0.0192	1.06	7.28e-02	2.90e-07	-5.40e-02	4.16e-04
7	165192	0.0404	1.04	6.93e-02	4.39e-08	-2.94e-02	4.16e-04

Table 17. Table 17: Classical Refinement strategy defined in § 4 applied to the linearized PBE for Methanol.

It.	$N$	Est. Err.	${\gamma }_{\mathrm{eff}}$	$E^r$	$E^m$	$E^{\Gamma }$	$E^{\partial \Omega }$
0	11769	-0.924	1.03	4.27e-01	4.67e-06	-1.35e+00	2.43e-06
1	12637	-0.613	1.03	1.25e-01	7.97e-07	-7.38e-01	2.41e-06
2	14798	-0.352	1.02	5.26e-02	4.42e-07	-4.04e-01	2.40e-06
3	20782	-0.267	1.01	-5.43e-02	-1.61e-06	-2.13e-01	2.40e-06
4	34556	-0.131	0.99	-3.32e-02	6.86e-08	-9.74e-02	2.40e-06
5	63134	-0.0388	0.94	1.25e-02	-2.33e-07	-5.12e-02	2.40e-06
6	118378	-0.0098	-	1.77e-02	-1.35e-07	-2.75e-02	2.40e-06

Table 18. Table 18: Classical Refinement strategy defined in § 4 applied to the nonlinear PBE for methanol.

It.	$N$	Est. Err.	${\gamma }_{\mathrm{eff}}$	$E^r$	$E^m$	$E^{\Gamma }$	$E^{\partial \Omega }$
0	11769	-0.924	1.03	4.27e-01	3.31e-06	-1.35e+00	2.41e-06
1	12637	-0.613	1.03	1.24e-01	2.05e-06	-7.38e-01	2.40e-06
2	14798	-0.352	1.02	5.22e-02	-1.86e-07	-4.04e-01	2.38e-06
3	20782	-0.267	1.01	-5.47e-02	-3.94e-07	-2.13e-01	2.38e-06
4	34538	-0.131	0.98	-3.35e-02	-3.19e-07	-9.75e-02	2.38e-06
5	63070	-0.0391	0.94	1.22e-02	1.45e-07	-5.13e-02	2.38e-06
6	118208	-0.01	-	1.75e-02	5.36e-08	-2.75e-02	2.38e-06

Table 19. Table 19: Uniform Contribution Refinement strategy applied to linearized PBE for the setup described in § 5.3.5 .

It.	$N$	Est. Err.	${\gamma }_{\mathrm{eff}}$	$E^r$	$E^m$	$E^{\Gamma }$	$E^{\partial \Omega }$
0	6718	-128	1.04	6.81e+00	1.92e-07	-1.35e+02	3.76e-14
1	15541	-28.2	1.04	8.51e+00	7.96e-05	-3.68e+01	4.10e-15
2	41760	-5.12	1.01	3.95e+00	1.94e-05	-9.07e+00	3.14e-15
3	141855	-0.0159	-	2.24e+00	1.66e-05	-2.25e+00	3.12e-15

Table 20. Table 20: Uniform Contribution Refinement strategy applied to nonlinear PBE for the setup described in § 5.3.5 .

It.	$N$	Est. Err.	${\gamma }_{\mathrm{eff}}$	$E^r$	$E^m$	$E^{\Gamma }$	$E^{\partial \Omega }$
0	6718	-131	1.04	3.72e+00	3.27e-07	-1.35e+02	2.25e-15
1	15541	-29.4	1.04	7.36e+00	7.20e-05	-3.68e+01	2.46e-16
2	41760	-6.24	1.01	2.83e+00	2.36e-05	-9.07e+00	1.88e-16
3	141855	-1.18	-	1.07e+00	-5.90e-07	-2.25e+00	1.83e-16

Table 21. Table 21: Adaptive Contribution Refinement strategy applied to linearized PBE for the setup described in § 5.3.5 .

It.	$N$	Est. Err.	${\gamma }_{\mathrm{eff}}$	$E^r$	$E^m$	$E^{\Gamma }$	$E^{\partial \Omega }$
0	6718	-128	1.04	6.81e+00	1.92e-07	-1.35e+02	3.76e-14
1	9458	-116	1.02	6.51e+00	9.61e-06	-1.22e+02	8.52e-15
2	14164	-59.5	1.02	5.72e+00	7.12e-05	-6.53e+01	4.42e-15
3	21205	-29.9	1.01	4.49e+00	-2.03e-05	-3.44e+01	3.30e-15
4	33141	-15.3	1.01	3.55e+00	2.44e-05	-1.89e+01	3.29e-15
5	53605	-7.32	0.99	2.85e+00	-4.94e-06	-1.02e+01	3.28e-15
6	91792	-3.3	0.95	2.30e+00	2.61e-05	-5.59e+00	3.15e-15
7	160005	-1.12	-	1.92e+00	1.41e-05	-3.04e+00	3.15e-15

Table 22. Table 22: Adaptive Contribution Refinement strategy applied to nonlinear PBE for the setup described in § 5.3.5 .

It.	$N$	Est. Err.	${\gamma }_{\mathrm{eff}}$	$E^r$	$E^m$	$E^{\Gamma }$	$E^{\partial \Omega }$
0	6718	-131	1.04	3.72e+00	3.27e-07	-1.35e+02	2.25e-15
1	9458	-118	1.02	4.23e+00	6.37e-06	-1.22e+02	5.13e-16
2	14164	-60.8	1.02	4.51e+00	1.48e-05	-6.53e+01	2.68e-16
3	21206	-30.9	1.01	3.59e+00	2.14e-06	-3.45e+01	1.98e-16
4	33216	-16.2	1.00	2.63e+00	2.77e-05	-1.88e+01	1.97e-16
5	53748	-8.24	0.98	1.90e+00	-3.05e-06	-1.01e+01	1.98e-16
6	92219	-4.27	0.95	1.30e+00	2.44e-06	-5.56e+00	1.88e-16
7	160740	-2.13	0.90	8.86e-01	-2.27e-06	-3.02e+00	1.87e-16

Table 23. Table 23: Classical Refinement strategy applied to linearized PBE for the setup described in § 5.3.5 .

It.	$N$	Est. Err.	${\gamma }_{\mathrm{eff}}$	$E^r$	$E^m$	$E^{\Gamma }$	$E^{\partial \Omega }$
0	6718	-128	1.04	6.81e+00	1.92e-07	-1.35e+02	3.76e-14
1	9481	-116	1.02	6.51e+00	9.62e-06	-1.22e+02	8.52e-15
2	14299	-57.9	1.02	5.69e+00	7.71e-05	-6.36e+01	4.10e-15
3	21638	-29.2	1.01	4.46e+00	1.38e-05	-3.37e+01	3.60e-15
4	33952	-14.8	1.01	3.49e+00	1.23e-05	-1.83e+01	3.58e-15
5	55747	-7.09	0.98	2.80e+00	2.70e-05	-9.89e+00	3.03e-15
6	95115	-3.16	0.95	2.25e+00	1.60e-05	-5.40e+00	2.91e-15

Table 24. Table 24: Classical Refinement strategy applied to nonlinear PBE for the setup described in § 5.3.5 .

It.	$N$	Est. Err.	${\gamma }_{\mathrm{eff}}$	$E^r$	$E^m$	$E^{\Gamma }$	$E^{\partial \Omega }$
0	6718	-131	1.04	3.72e+00	3.27e-07	-1.35e+02	2.25e-15
1	9477	-118	1.02	4.23e+00	8.49e-06	-1.22e+02	5.12e-16
2	14286	-59.8	1.02	4.49e+00	1.86e-05	-6.43e+01	2.47e-16
3	21365	-30.4	1.01	3.61e+00	2.82e-05	-3.40e+01	2.21e-16
4	33470	-16	1.00	2.60e+00	1.84e-05	-1.86e+01	2.18e-16
5	54450	-8.16	0.98	1.88e+00	6.44e-07	-1.00e+01	1.85e-16
6	93369	-4.22	0.95	1.28e+00	-5.89e-07	-5.50e+00	1.77e-16

Equations

$$\{\begin{array}{l}\begin{array}{r}-\nabla \cdot (ϵ(x)\nabla \tilde{u}(x))+{\bar{\kappa }}^2(x)\left(\frac{k_{B}T}{e_c}\right)\sinh \left(\frac{e_c\tilde{u}(x)}{k_{B}T}\right)=4\pi \sum _{i=1}^m{Q}_i\delta (x-x_i),\\ \underset{\Vert x\Vert \to \infty }{\lim }\tilde{u}(x)=0.\end{array}\end{array}$$

$$ϵ(x)=\{\begin{array}{ll}{ϵ}_m& x\in {\Omega }_m\\ {ϵ}_s& x\in {\Omega }_s^{\infty }\end{array}\text{and}{\bar{\kappa }}^2(x)=\{\begin{array}{ll}0& x\in {\Omega }_m\\ {\bar{\kappa }}_s^2={ϵ}_s\frac{8\pi N_A{e}_c^2}{1000k_{B}T}{I}_s& x\in {\Omega }_s^{\infty }\end{array}.$$

$$\{\begin{array}{l}\begin{array}{rl}-\nabla \cdot (ϵ(x)\nabla u(x))+{\bar{\kappa }}^2(x)\sinh (u(x))=\frac{4\pi e_c}{k_{B}T}\sum _{i=1}^m{Q}_i\delta (x-x_i),x& \in \Omega ,\\ u(x)=g(x),x& \in \partial \Omega .\end{array}\end{array}$$

$$g=\frac{e_c}{k_{B}T}\sum _{i=1}^m\frac{Q_i}{{ϵ}_s|x-x_i|}\exp \left(\frac{-{\bar{\kappa }}_s|x-x_i|}{\sqrt{{ϵ}_s}}\right).$$

$$u=\{\begin{array}{ll}{u}^s+u^h+u^r& \text{ in }{\Omega }_m,\\ u^r& \text{ in }{\Omega }_s,\end{array}$$

$$\{\begin{array}{l}\begin{array}{r}-\nabla \cdot {ϵ}_m\nabla u^s=\frac{4\pi e_c}{k_{B}T}\sum _{i=1}^m{Q}_i\delta (x-x_i),\\ u^s(\infty )=0.\end{array}\end{array}$$

$$u^s(x)=\frac{e_c}{{ϵ}_m{k}_{B}T}\sum _{i=1}^m\frac{Q_i}{|x-x_i|}.$$

$$\{\begin{array}{l}\begin{array}{rl}{\nabla }^2{u}^h=0& \text{ in }{\Omega }_m,\\ u^h=-u^s& \text{ on }\Gamma .\end{array}\end{array}$$

$$\begin{cases}\begin{aligned} -\nabla\cdot(\epsilon(x)\nabla u^{r}(x))+\bar{\kappa}^{2}(x)\sinh(u^{r}(x))=0\quad\text{ in }\Omega_{m}\cup\Omega_{s},\\ \left\llbracket u^{r}\right\rrbracket_{\Gamma}=0,\\ \left\llbracket\epsilon(x)\frac{\partial{u^{r}}(x)}{\partial\mathbf{n}}\right\rrbracket_{\Gamma}=-\epsilon_{m}\frac{\partial u^{s}+u^{h}}{\partial\mathbf{n}}\\ u=g\quad\text{ on }\partial\Omega,\end{aligned}\end{cases}$$

$$\left\llbracket v(x)\right\rrbracket_{\Gamma}=\lim_{\alpha\rightarrow 0^{+}}v(x+\alpha\mathbf{n})-v(x-\alpha\mathbf{n}),$$

$$\begin{cases}\begin{aligned} -\nabla\cdot(\epsilon(x)\nabla u^{r}(x))+\bar{\kappa}^{2}(x)N(u^{r}(x))=0\quad\text{ in }\Omega_{m}\cup\Omega_{s},\\ \left\llbracket u^{r}\right\rrbracket_{\Gamma}=0,\\ \left\llbracket\epsilon(x)\frac{\partial{u^{r}}(x)}{\partial\mathbf{n}}\right\rrbracket_{\Gamma}=-\epsilon_{m}\frac{\partial u^{s}+u^{h}}{\partial\mathbf{n}},\\ u=g\quad\text{ on }\partial\Omega.\end{aligned}\end{cases}$$

$$N(u^r(x))=\{\begin{array}{l}\begin{array}{rlrl}& \sinh (u^r(x))& & \text{for nonlinear PBE},\\ & u^r(x)& & \text{for linearized PBE}.\end{array}\end{array}$$

$$H_{u^s}^1({\Omega }_m):=\{v\in H^1({\Omega }_m):v(x)=-u^s\text{ on }\Gamma \}$$

$$H_g^1(\Omega ):=\{v\in H^1(\Omega ):v=g\text{ on }\Gamma ,\alpha \le v\le \beta \text{ in }\Omega \}.$$

$$\{\begin{array}{l}\begin{array}{rl}& (ϵ\nabla u^h,\nabla w{)}_m=0\\ & (ϵ\nabla u^r,\nabla v)+({\bar{\kappa }}^{2}N(u^r),v)+⟨{ϵ}_m\frac{\partial u^h}{\partial \mathbf{n}},v{⟩}_{\Gamma }=-⟨{ϵ}_m\frac{\partial u^s}{\partial \mathbf{n}},v{⟩}_{\Gamma }\end{array}\end{array}$$

$$⟨{ϵ}_m\frac{\partial u^h}{\partial \mathbf{n}},v⟩=(ϵ{\nabla }^2{u}^h,v{)}_m+(ϵ\nabla u^h,\nabla v{)}_m$$

$$\{\begin{array}{l}\begin{array}{rl}& (ϵ\nabla u^h,\nabla w{)}_m=0\\ & (ϵ\nabla u^r,\nabla v)+({\bar{\kappa }}^{2}N(u^r),v)+(ϵ\nabla u^h,\nabla v{)}_m=-⟨{ϵ}_m\frac{\partial u^s}{\partial \mathbf{n}},v{⟩}_{\Gamma },\end{array}\end{array}$$

$$\Delta G_{\text{sol}}=\frac{\alpha }{2}\int \sum _{i=1}^m{Q}_i\delta (x-x_i)(u^h(x)+u^r(x))\mathrm{d}x,$$

$$Q(u^h,u^r)={\int }_{{\Omega }_m}\sum _{i=1}^m{Q}_i{\eta }^{-3}\rho ((x-x_i)/\eta )(u^h(x)+u^r(x))\mathrm{d}x=(\psi ,u^h+u^r{)}_m,$$

$$\psi (x)=\sum _{i=1}^m{Q}_i{\eta }^{-3}\rho ((x-x_i)/\eta ),$$

$$\rho (x)=\{\begin{array}{l}c{e}^{(-1/(1-|x{|}^2))}\text{ if }|x|<1,\\ 0\text{otherwise},\end{array}$$

$$\{\begin{array}{l}\begin{array}{rl}& (ϵ\nabla U^h,\nabla w{)}_m=0\\ & (ϵ\nabla U^r,\nabla v)+({\bar{\kappa }}^{2}N(U^r),v)+(ϵ\nabla U^h,\nabla v{)}_m=-⟨{ϵ}_m\frac{\partial u^s}{\partial \mathbf{n}},v{⟩}_{\Gamma },\end{array}\end{array}$$

$$\text{set@color}⟨\text{leavevmode}\text{set@color}⟨\mathcal{L}x,y^{\ast }{\text{set@color}⟩\text{leavevmode}\text{set@color}⟩}_Y=\text{set@color}⟨\text{leavevmode}\text{set@color}⟨x,{\mathcal{L}}^{\ast }{y}^{\ast }{\text{set@color}⟩\text{leavevmode}\text{set@color}⟩}_X,x\in X,y^{\ast }\in Y^{\ast },$$

$$\text{set@color}⟨\text{leavevmode}\text{set@color}⟨\psi ,u-U{\text{set@color}⟩\text{leavevmode}\text{set@color}⟩}_X=\text{set@color}⟨\text{leavevmode}\text{set@color}⟨L^{\ast }\varphi ,u-U{\text{set@color}⟩\text{leavevmode}\text{set@color}⟩}_X=\text{set@color}⟨\text{leavevmode}\text{set@color}⟨\varphi ,Lu-LU{\text{set@color}⟩\text{leavevmode}\text{set@color}⟩}_Y=\text{set@color}⟨\text{leavevmode}\text{set@color}⟨\varphi ,f-LU{\text{set@color}⟩\text{leavevmode}\text{set@color}⟩}_Y.$$

$$\{\begin{array}{l}\begin{array}{rl}& (ϵ\nabla {\varphi }^r,\nabla v)+({\bar{\kappa }}^2{\varphi }^r,v)=(\psi ,v{)}_m,\\ & (ϵ\nabla {\varphi }^h,\nabla w{)}_m+(ϵ\nabla {\varphi }^r,\nabla w{)}_m=(\psi ,w{)}_m,\end{array}\end{array}$$

$$\{\begin{array}{l}\begin{array}{rl}& (ϵ\nabla w^r,\nabla v)+({\bar{\kappa }}^2{w}^r,v)=(\psi ,v{)}_m,\\ & (ϵ\nabla w^h,\nabla w{)}_m+⟨w^h,{ϵ}_m\frac{\partial w}{\partial \mathbf{n}}{⟩}_{\Gamma }=(\psi ,w{)}_m,\end{array}\end{array}$$

$$u^r=u_0^r+u_d^r\text{and}{u}^h=u_0^h+u_d^h$$

$$U^r=U_0^r+U_d^r\text{and}{U}^h=U_0^h+U_d^h$$

$$Q(u^h-U^h,u^r-U^r)=E^r+E^m+E^{\Gamma }+E^{\partial \Omega }+E^{\mathrm{neg}},$$

$$E^r$$

Full text

1 Introduction

Electrostatic interactions play a critical role in determining macroscopic properties of dielectric biomolecular systems, such as solvation free energy and binding affinities [26, 46, 61]. The Poisson-Boltzmann equation (PBE) has been widely used for modeling the electrostatic interactions of charged bodies such as molecules and proteins in electrolyte solvents. The PBE was introduced decades ago [39, 15], and we refer to the classical texts [54, 60] for its derivation.

The focus of this article is robust error estimation and refinement strategies for computing a quantity of interest (QoI), such as the solvation free energy, from the solution of the PBE. The PBE is a challenging equation to solve numerically and numerous computational methods and software packages have been derived for its solution [38, 58, 64, 62, 67, 11, 25, 43, 51, 9, 20, 56, 24, 42, 5, 59, 21, 50, 12, 41]. In this article we follow the approach in [36, 41] to solve the PBE using a three term splitting method which accounts for the well-posedness of the continuum problem as well as avoiding amplification of numerical rounding errors. However, even this method, like all numerical methods, often has significant errors in the computation of the QoI and this error needs to be accurately estimated from computed information for reliable use of the PBE in biophysics, biochemistry, medical and other science and engineering fields [32, 31].

In this article we employ adjoint based a posteriori analysis to accurately quantify the error in a QoI computed from the numerical solution of the PBE. Adjoint based error estimation is widely used for a host of numerical methods including finite elements, finite difference, time integration, multi-scale simulations and inverse problems [30, 29, 32, 1, 7, 8, 10, 37, 13, 16, 22, 57]. The error estimate weights computable residuals of the numerical solution with the solution of an adjoint problem to quantify the accumulation and propagation of error. The resulting estimates have the useful feature that the total error is decomposed as a sum of contributions from various aspects of the discretization and therefore provide insight in to the effect of different choices for the parameters controlling the discretization. Thus, we not only quantify the error using adjoint based a posteriori analysis, we also partition the error to identify contributions from various sources of error. For example, we can identify if the boundary discretization or the interior discretization is the major source of error.

Since the error in the numerical solution of the PBE is often significant, there are a number of adaptive refinement strategies proposed for obtaining accurate solutions of the PBE [12, 42, 41, 63, 2]. Most of the adaptive algorithms are based on controlling the error in global norms and some of the algorithms are shown to be provably convergent [20]. However, if the goal of the numerical computation is accurate approximation of the QoI, then a refinement strategy based on solution residuals weighted by the adjoint information is an appealing option. In this paper we propose refinement strategies based on the relative contribution to the error of a discretization choice. The adjoint based analysis and its partitioning of the error suggests novel refinement strategies for obtaining accurate estimates of the QoI from the numerical solution of the PBE.

Adaptive refinement using adjoint based analysis and optimal multilevel preconditioning for the PBE are developed previously in [2]. However, the analysis and results of this article differ significantly from that paper. The focus of that paper was adaptive refinement for the linearized PBE using the two term splitting [20], whereas we focus on the three term splitting for both the linearized and the nonlinear PBE [36, 41]. Moreover, our aim is to derive accurate error estimates in the QoI. While [2] derived an error estimate for the QoI for the three term splitting, no numerical experiments were performed for the three term splitting. Even for the two term splitting no numerical results indicating the accuracy of the derived estimates were shown and instead the focus was on adaptive refinement. In addition, the adjoint problem for the three term splitting derived in [2] leads to an ill-posed problem as we discuss in §3. In this work, we not only derive an estimate for the three term splitting based on the correct formulation of the adjoint operator, we also decompose the error so that the various sources of error and their relative contributions are also available. Moreover, the error estimate derived in [2] assumes that the continuum and discrete solutions satisfy the boundary conditions exactly. While this assumption may be justifiable for the results in the two term splitting in [2], we point out the importance of the role of boundary condition for the harmonic component of the three term splitting. This boundary condition is defined on the interface between the solvent and molecular regions, and hence impacts the computation of the QoI significantly. Finally, we propose a fundamentally different refinement strategy since the standard goal oriented refinement strategy employed in [2] appears to be sub-optimal. Adaptive refinement for obtaining accurate values of a QoI is a challenging task as the error contributions of an individual element may be positive or negative leading to significant cancellation of error. In [2], the refinement strategy takes the absolute value of error contributions and applies the principle of equidistribution for marking elements for refinement. This strategy ignores the cancellation of error, and hence the resulting adaptive algorithm may have less than desirable convergence properties. This drawback is overcome in [2] by defining a somewhat ad-hoc error indicator. On the other hand, this article decomposes the error into different contributions and use this information to devise adaptive schemes to target the discretization choices which have the most significant effect on error.

The rest of the paper is organized as follows. Section 2 introduces the PBE, its linearized and nonlinear versions, weak forms and a finite element method to solve it. Section 3 performs adjoint based a posteriori analyses for both the linearized and nonlinear PBE. In particular, a well-posed adjoint problem for the three-term PBE and error representations are derived. Section 4 discusses refinement strategies based on a posteriori error estimates. Numerical experiments are presented in Section 5, which illustrate the accuracy of the estimates as well as the efficacy of employing refinement strategies which target specific sources of error. Section 6 presents conclusions.

2 The Poisson-Boltzmann Equation

2.1 The nonlinear Equation and its dimensionless form

The Poisson-Boltzmann equation models the electrostatic activity between molecules in an ionic solvent. In this model, it is assumed that the ions in the solvent are distributed according to the Boltzmann distribution and that the potential of the mean force on a particle is simply the charge of the ion times the electrostatic potential. For a 1:1 electrolyte solvent (e.g. NaCl), the nonlinear Poisson-Boltzmann equation is [35, 4],

[TABLE]

Here, $\tilde{u}$ is the unknown electrostatic potential, $\epsilon$ is the dielectric coefficient, $\bar{\kappa}(x)$ is the modified Debye-Hückel parameter which describes the accessibility of the solvent to the solute, $k_{B}$ is the Boltzmann constant, $e_{c}$ is the charge on a proton and $T$ is the temperature. Moreover, the solute contains a total of $m$ fixed point charges, with the $i$th charge $Q_{i}$ centered at position $x_{i}$. The resulting distribution is a linear combination of Dirac delta functions $\delta(x-x_{i})$.

The domain for the problem $\mathbb{R}^{3}$, is subdivided into a molecular region, $\Omega_{m}$, a solvent region $\Omega_{s}^{\infty}$, and an interface between the two denoted by $\Gamma$. The solute is surrounded by solvent, which is represented as a continuum over the subdomain $\Omega_{s}^{\infty}=\mathbb{R}^{3}\backslash\overline{\Omega}_{m}$. The subdomains for a typical biomolecular solute are shown in Figure 1 which has been adopted from [12].

The dielectric coefficient $\epsilon(x)$ and modified Debye-Hückel parameter $\bar{\kappa}(x)$ are defined on $\Omega_{m}\cup\Omega_{s}^{\infty}$ by the piecewise constant functions

[TABLE]

Here, $\epsilon_{m}$ and $\epsilon_{s}$ are positive constants and $N_{A}$ is Avogadro’s number. The ionic strength $I_{S}$ is a physical parameter which varies depending on the solvent.

Numerical simulations are not feasible on the unbounded solvent domain, $\Omega_{s}^{\infty}$, and hence it is truncated at a finite radius from the “center” of the molecule, to form a bounded domain $\Omega_{s}$. Dirichlet boundary conditions are imposed to capture the asymptotic behavior of the solution on an unbounded domain. Combining this with the change of variables, $u(x)=e_{c}\tilde{u}(x)/k_{B}T$, results in a dimensionless Poisson-Boltzmann equation on the spherical domain $\Omega=\Omega_{m}\cup\Omega_{s}\cup\Gamma$:

[TABLE]

The boundary conditions are prescribed using a linear combination of Helmholtz Green’s functions [12],

[TABLE]

2.2 Weak form based on three term splitting

We denote by $L_{2}(\Omega)$ as the space of square integrable functions, $H^{1}(\Omega)$ as the space of functions having an integrable (weak) derivative, $H^{1}_{0}(\Omega)$ as the subspace of $H^{1}(\Omega)$ of functions satisfying homogeneous Dirichlet boundary conditions (in the sense of the trace operator) and $H^{-1}$ as the dual space of $H^{1}_{0}(\Omega)$. The right hand side of (3) contains $\delta$ functions, which are unbounded linear functionals over the space $H^{1}_{0}$ and hence a well-posed weak form cannot be derived directly from (3). To overcome this problem, two and three term splittings of the PBE have been proposed [68, 20, 41]. The two and three term splitting are equivalent mathematically, however, the three term splitting is numerically more desirable [41]. The three term splitting decomposes the function $u$ as

[TABLE]

where $u^{s}$, $u^{h}$ and $u^{r}$ are the singular, harmonic and regular components respectively. The singular function $u^{s}$ is the solution of the following Poisson equation

[TABLE]

Recognizing that the singular component is the Green’s function of the Laplace operator leads to an analytical expression for $u^{s}$ as

[TABLE]

The harmonic component $u^{h}$ is the solution to

[TABLE]

The regular component $u^{r}$ satisfies

[TABLE]

where the jump at the interface is defined as

[TABLE]

with $\mathbf{n}$ as the unit normal to the interface $\Gamma$, pointing outward from $\Omega_{m}$. The condition involving the jump in the normal derivative of $u^{r}$ arises by substituting (5) in (3), using the definitions of $u^{s}$ and $u^{h}$, and the fact that for the solution $u$ of (3) we have \left\llbracket\epsilon(x)\frac{\partial{u}(x)}{\partial\mathbf{n}}\right\rrbracket=0. Sometimes the nonlinear PBE is linearized by the assumption $\sinh(u)\approx u$ leading to the dimensionless linearized PBE. We can write both the linear and nonlinear versions as

[TABLE]

where

[TABLE]

2.3 Weak forms

We define the affine spaces

[TABLE]

and

[TABLE]

Here $\alpha$ and $\beta$ are positive constants used to control the nonlinear $\sinh$ term, see [41] for details. The weak form for the three term split PBE, (8) and (10), is to find $(u^{h},u^{r})\in H^{1}_{u^{s}}(\Omega_{m})\times H^{1}_{g}(\Omega)$ such that

[TABLE]

for all $(w,v)\in H^{1}_{0}(\Omega_{m})\times H^{1}_{0}(\Omega)$. Here we used the notation $(a,b)=\int_{\Omega}ab\,\mathrm{d}x$, $(a,b)_{m}=\int_{\Omega_{m}}ab\,\mathrm{d}x$ and $\langle a,b\rangle_{\Gamma}=\int_{\Gamma}ab\,\mathrm{d}x$ to represent the standard $L_{2}$ inner products over $\Omega$, $\Omega_{m}$ and $\Gamma$ respectively. The existence and uniqueness of the weak solution is shown in [41]. The weak form (14) is a one-way coupled system; we first solve for $u^{h}$ and then use it to compute $u^{r}$. Now using the Green’s identity

[TABLE]

and (8) in (14) leads to an different weak form: find $(u^{h},u^{r})\in H^{1}_{u^{s}}(\Omega_{m})\times H^{1}_{g}(\Omega)$ such that

[TABLE]

for all $(w,v)\in H^{1}_{0}(\Omega_{m})\times H^{1}_{0}(\Omega)$. The weak forms (14) and (16) are mathematically equivalent, however, the form (16) is amenable to defining the adjoint operator as discussed in §3.

2.4 Quantity of interest: solvation free energy

The QoI may be any bounded linear functional of the weak solution $(u^{h},u^{r})$. An important physical quantity computed from the solution of the PBE is electrostatic free energy of solvation [36],

[TABLE]

where $\alpha=k_{B}T/e_{c}$. Unfortunately, $\Delta G_{\text{sol}}$ is not a bounded linear functional in $H^{1}_{0}(\Omega)$ due to the presence of $\delta$ functions. A common approach is to “mollify” the unbounded functional [2, 7, 1] to obtain a bounded linear functional. We thus define our quantity of interest to be a mollified version of solvation free energy, scaled by $2/\alpha$ for simplicity, as

[TABLE]

where

[TABLE]

$\rho$ is the standard mollifier

[TABLE]

$|x|$ denotes the Euclidean norm of $x\in\mathbb{R}^{3}$ and $c$ is a scaling constant to ensure that $\int_{\mathbb{R}^{3}}\rho(x)=1$. Now, as $\eta\rightarrow 0$, $\eta^{-3}\rho(x/\eta)\rightarrow\delta(x)$. Hence the value of the QoI approaches the value of the (scaled) solvation free energy for small values of $\eta$.

2.5 Finite element method

We discretize $\Omega$ and $\Omega_{m}$ into three dimensional triangulations $\mathcal{T}$ and $\mathcal{T}_{m}$. We assume that the interface $\Gamma$ is polygonal and exactly represented by the triangulation. Although the triangulations $\mathcal{T}$ and $\mathcal{T}_{m}$ may differ in $\Omega_{m}$, they respect the interface $\Gamma$ in the sense that $(\cup_{T_{m}\in\mathcal{T}_{m}}{T}_{m})\cap\Gamma=(\cup_{T\in\mathcal{T}}{T})\cap\Gamma$. Each of these triangulations is arranged in such a way that the union of the elements of $\mathcal{T}$ (resp. $\mathcal{T}_{m}$ ) is $\Omega$ (resp. $\Omega_{m}$) and the intersection of any two elements is either a common edge, node, or is empty. The finite element space consists of continuous piecewise polynomials. We let $V_{h}\subset H_{0}^{1}(\Omega)$ (resp. $V_{h}^{m}\subset H_{0}^{1}(\Omega_{m})$) denote the space of continuous piecewise polynomial functions $v(x)\in\mathbb{R}$ defined on $\mathcal{T}$ (resp. $\mathcal{T}_{m})$. Similarly, we let $V_{h,g}$ (resp. $V_{h,u_{s}}^{m})$ be the affine space of continuous piecewise polynomial functions $v(x)\in\mathbb{R}$ such that $v(x)=g(x)$ for $x$ on $\partial\Omega$ (resp. $v(x)=-u_{s}(x)$ for $x$ on $\Gamma=\partial\Omega_{m}$).

The discrete weak problem is to find $(U^{h},U^{r})\in V_{h,u_{s}}^{m}\times V_{h,g}$ such that

[TABLE]

for all $(w,v)\in V_{h}^{m}\times V_{h}$. Note that throughout this article we use lower case letters for continuum solutions and uppercase letters for discrete solutions.

3 Adjoint based a posteriori analysis

In this section we derive the adjoint equation corresponding to the PBE and then form error representations for both the linearized and nonlinear PBE.

3.1 Abstract definition of adjoint operator and error representation

The adjoint operator $\mathcal{L}^{\ast}:Y^{\ast}\rightarrow X^{\ast}$ of a linear operator $\mathcal{L}:X\rightarrow Y$ between Banach spaces $X,Y$ with dual spaces $X^{\ast},Y^{\ast}$ is defined by the bilinear identity [53, 44, 65],

[TABLE]

where \mathopen{\hbox{\set@color{\langle}}\kern-1.94444pt\leavevmode\hbox{\set@color{\langle}}}\cdot,\cdot\mathclose{\hbox{\set@color{\rangle}}\kern-1.94444pt\leavevmode\hbox{\set@color{\rangle}}}_{S} denotes duality-pairing in the space $S\in\{X,Y\}$. Now, if $Lu=f$ and $L^{\ast}\phi=\psi$, and $U$ is a discrete approximation to $u$, we obtain a representation for the error $(u-U)$ as

[TABLE]

The above abstract error representation is a standard form for all adjoint based error analysis: residual(s) of the discrete solution weighted by the adjoint solution(s). The weighting of the residual by the adjoint solution accounts for the accumulation and cancellation of error in the discrete solution. We remark that the derivation in (23) is similar to the derivation of standard Green’s functions in PDE analysis, and hence adjoint solutions may be thought of as generalized Green’s functions [33].

3.2 A posteriori analysis of the linearized PBE

This section forms an adjoint operator and an error representation for the linearized PBE.

3.2.1 Adjoint operator for the linearized PBE

In the context of the linearized PBE, the duality pairing \mathopen{\hbox{\set@color{\langle}}\kern-1.94444pt\leavevmode\hbox{\set@color{\langle}}}\mathcal{L}x,y^{\ast}\mathclose{\hbox{\set@color{\rangle}}\kern-1.94444pt\leavevmode\hbox{\set@color{\rangle}}}_{Y} is described by the left hand side of (16) with $N(u)=u$. Applying the definition (22) leads to the following adjoint problem: find $(\phi^{h},\phi^{r})\in H^{1}_{0}(\Omega_{m})\times H^{1}_{0}(\Omega)$ such that

[TABLE]

for all $w,v\in H^{1}_{0}(\Omega_{m})\times H^{1}_{0}(\Omega)$. Here $\psi$ arises from the definition of the QoI, see (19). Observe that (24) is also a one way coupled system, similar to (16), however, the direction of coupling is now reversed: we first solve for the component in $\Omega$ and use that in the equation posed on $\Omega_{m}$.

Remark 1**.**

In Section 4.2 of [2], an adjoint to the three term split PBE is defined as: find $(w^{h},w^{r})\in H^{1}_{0}(\Omega_{m})\times H^{1}_{0}(\Omega)$ such that

[TABLE]

*for all $(w,v)\in H^{1}_{0}(\Omega_{m})\times H^{1}_{0}(\Omega)$. However, this is not a well-posed problem as $\langle\epsilon_{m}\frac{\partial w}{\partial\mathbf{n}},w^{h}\rangle_{\Gamma}$ is not continuous in $w$ for all $w\in H^{1}_{0}(\Omega_{m})$. Continuity of $\langle\epsilon_{m}\frac{\partial w}{\partial\mathbf{n}},w^{h}\rangle_{\Gamma}$ requires addition regularity on $w$ e.g. $w\in H^{\frac{3}{2}}(\Omega)$. *

3.2.2 Error representation for the linearized PBE

The effect on approximating the boundary conditions on the interface $\Gamma$ for $U^{h}$ may have significant effect on the accuracy of the method. Hence, we quantify the effect of boundary conditions, both at $\Gamma=\partial\Omega_{m}$ corresponding to the harmonic component $U^{h}$ and at $\partial\Omega$ corresponding to the regular component $U^{r}$. We employ the decompositions

[TABLE]

where $u_{0}^{r}\in H^{1}_{0}(\Omega)$ (resp. $u_{0}^{h}\in H^{1}_{0}(\Omega_{m})$) and ${u}_{d}^{r}\in H^{1}(\Omega)$ (resp. ${u}_{d}^{h}\in H^{1}(\Omega_{m})$ ) such that ${u}_{d}^{r}=g$ on $\partial\Omega$ (resp. ${u}_{d}^{h}=-u^{s}$ on $\partial\Omega_{m}=\Gamma$) . Similarly we have the decompositions,

[TABLE]

where $U^{r}_{0}\in V_{h}$ (resp. $U^{h}_{0}\in V_{h}^{m}$) and ${U}_{d}^{r}\in V_{h,g}$ (resp. ${U}_{d}^{r}\in V_{h,u_{s}}^{m}$). Note that due to the finite dimension of $V_{h,g}$ and $V_{h,u_{s}}^{m}$ and the nature of the boundary conditions $g$ and $-u^{s}$, ${u}_{d}^{r}\neq{U}_{d}^{r}$ and ${u}_{d}^{h}\neq{U}_{d}^{h}$. Moreover, there are infinitely many choices for the functions $u_{d}^{h},u_{d}^{r},U_{d}^{h},U_{d}^{h}$ and we assume a choice is made such that these functions are known. This leads to the following error representation.

Theorem 1**.**

Let $(u^{h},u^{r})$ be the true solutions to the linearized PBE (16) with $N(u^{r})=u^{r}$, $(U^{h},U^{r})$ be the finite element solutions to the discrete weak form (21) and $(\phi^{h},\phi^{r})$ be the solutions to the adjoint weak form (24). Then the error in the QoI (18) is given by

[TABLE]

where

[TABLE]

Proof.

The tuple $(u^{h}-U^{h},u^{r}-U^{r})$ is not in $H^{1}_{0}(\Omega_{m})\times H^{1}_{0}(\Omega)$. However, if we use the decompositions (26) and (27) along with the linearity of the QoI $Q$,

[TABLE]

The tuple $((u_{0}^{h}+U^{h}_{d})-U^{h},(u^{r}_{0}+U^{r}_{d})-U^{r})$ is in $H^{1}_{0}(\Omega_{m})\times H^{1}_{0}(\Omega)$. Hence, setting $w=((u_{0}^{h}+U^{h}_{d})-U^{h})$ and $v=((u^{r}_{0}+U^{r}_{d})-U^{r})$ in the adjoint equation (24) and adding the two equations leads to,

[TABLE]

Substituting $u^{r}_{0}=(u^{r}_{0}+u^{r}_{d})-u^{r}_{d}=u^{r}-u^{r}_{d}$ and similarly $u^{h}_{0}=u^{h}-u^{h}_{d}$ and rearranging,

[TABLE]

Now, since $(u^{h},u^{r})$ is the true solution, it satisfies the weak form (16). Substituting this in (32) and rearranging terms,

[TABLE]

Combining (30) and (33) completes the proof.

∎

In the above theorem $E^{r}$, $E^{m}$, $E^{\Gamma}$, $E^{\partial\Omega}$ and $E^{\mathrm{neg}}$ denote different sources of error. The first four terms $E^{r}$, $E^{m}$, $E^{\Gamma}$ and $E^{\partial\Omega}$ have the form of adjoint weighted residuals and reflect error contributions due to FEM solution of $u^{r}$, FEM solution of $u^{h}$, representation of boundary data for $u^{h}$ and representation of boundary data for $u^{r}$. The term $E^{\mathrm{neg}}$, which is computable since all the functions involved are known, is referred to as the “negligible” component of error as it is typically negligible due to the standard choice of the boundary functions. See §5 for more details on the choice of the boundary functions involved as well as the numerical value of this term.

3.3 A posteriori Analysis of the nonlinear PBE

We now extend the ideas for the linearized PBE to derive an adjoint and error representation for the nonlinear PBE.

3.3.1 Adjoint operator for the nonlinear PBE

The extension of the above approach to the nonlinear PBE is complicated by the fact that there is no unique definition of an adjoint operator corresponding to a nonlinear differential operator. Rather, an adjoint problem useful for the purpose at hand has to be selected. A common choice useful for various kinds of analysis is based on linearization [53, 52]. Defining $z=su+(1-s)U$ and $\overline{\alpha}(x)=\int_{0}^{1}\cosh(z(x))\,\mathrm{d}s$, we observe that

[TABLE]

Then the adjoint corresponding to the nonlinear PBE (16) is: find $(\phi^{h},\phi^{r})\in$    $H^{1}_{0}(\Omega_{m})\times H^{1}_{0}(\Omega)$ such that

[TABLE]

for all $w,v\in H^{1}_{0}(\Omega_{m})\times H^{1}_{0}(\Omega)$. In practice, we cannot compute the linearization $\overline{\alpha}$ since we do not know the true solution $u^{r}$. Instead, the differential operator is typically linearized around the numerical solution, in this case $U^{r}$. The resulting estimate can be shown to converge to the true estimate in the limit of refined discretization [32]. In practice, this approach yields robustly accurate error estimates.

3.3.2 Error representation for the nonlinear PBE

The above adjoint equation leads to the following error representation for the nonlinear PBE.

Theorem 2**.**

Let $(u^{h},u^{r})$ be the true solutions to the nonlinear PBE (16) with $N(u^{r})=\sinh(u^{r})$, $(U^{h},U^{r})$ be the finite element solutions to the discrete weak form (21) and $(\phi^{h},\phi^{r})$ be the solutions to the adjoint weak form (24). Then the error in the QoI (18) is given by,

[TABLE]

where

[TABLE]

Proof.

The proof is similar to the proof of Theorem 1. The difference is in the term $(\bar{\kappa}^{2}\phi^{r},(u^{r}_{0}+U^{r}_{d})-U^{r}))$ in (31) which now becomes,

[TABLE]

where we again made the substitution $u^{r}_{0}=u^{r}-u^{r}_{d}$. Combining the above equation with (34) leads to

[TABLE]

∎

3.3.3 Error representation for the alternate formulation of the PBE

In this article, the focus is on quantifying the error due to the solution of the FEM problem in (21) which corresponds to the solution of (16). However, some existing codes may be based on the discrete solution of the weak form (14). In such a case, the error representation is easily modified as shown in the next theorem.

Theorem 3**.**

Let $(u^{h},u^{r})$ be the true solutions to the (linearized or nonlinear) PBE (16), $(U^{h},U^{r})$ be the finite element solutions to the discrete weak form corresponding to (14) and $(\phi^{h},\phi^{r})$ be the solutions to the adjoint weak form (24). Then the error in the QoI (18) is given by,

[TABLE]

where

[TABLE]

and the remaining terms are the same as in Theorem 2.

Proof.

Adding and subtracting $\langle\epsilon_{m}\frac{\partial U^{h}}{\partial\mathbf{n}},\phi^{r}\rangle_{\Gamma}$ to term $E^{r}$ in Theorem 2 completes the proof. ∎

4 Refinement strategies based on a posteriori error estimates

This section discusses the accuracy of a posteriori error estimates and the potential for obtaining accurate QoI values using the error information to refine the discretization.

4.1 A posteriori error estimates: implementation and accuracy

The error representations (28) and (36) involve analytic adjoint solutions $(\phi^{h},\phi^{r})$ and representation of boundary conditions by functions $(u^{h}_{d},u^{r}_{d})$. In practice, these quantities need to be estimated computationally. As is common in literature for adjoint based a posteriori analysis, the adjoint solutions are approximated in a space $W^{h}$ (resp. $W^{r}$) such that $V^{h}\subset W^{h}$ (resp. $V^{r}\subset W^{r}$) [32, 30, 28, 19, 16, 32, 30, 28, 19, 16, 23, 19, 14, 8]. $W^{h}$ may be obtained by refining the mesh or by increasing the polynomial order. Similarly, the functions $(u^{h}_{d},u^{r}_{d})$ are approximated in $W^{h}$, such that they satisfy the boundary condition exactly on a boundary vertex and are zero on the interior vertices. These approximations lead to error estimates from the error representations (28) and (36). Since the formulas are similar, except that the $(\phi^{h},\phi^{r})$ and $(u^{h}_{d},u^{r}_{d})$ are replaced by their approximations, we avoid re-writing the error estimates explicitly and instead now refer to (28) and (36) as error estimates.

The accuracy of the error estimate is measured by the effectivity ratio defined as

[TABLE]

An accurate error estimator has an effectivity ratio close to one. Since the true solution is not known, we compute a more accurate reference numerical solution using a higher dimensional space for measuring the true error.

4.2 Guiding refinement decisions using error estimates

4.2.1 Error Contributions and cancellation of error

Once the error estimate is in place, its various components $E^{r},E^{m},E^{\Gamma}$ and $E^{\partial\Omega}$ reflect different sources of error. Refinement strategies based on these components can then be derived. For example, if $E^{\Gamma}$ is the dominant component, then simplices in $\mathcal{T}_{m}$ which intersect with the interface $\Gamma$ may be refined to reduce the error. This strategy of refining the mesh is quite different from classical adaptive refinement schemes. One main difference is that, in refining the simplices on the interface to reduce $E^{\Gamma}$ we may use either uniform refinement or an adaptive refinement strategy. The other difference is in the treatment of cancellation of errors which we now discuss.

Classical adaptive refinement schemes form elemental error indicators and refine elements which have the largest value of such indicators [1, 55, 34]. While adaptive refinement often outperforms uniform refinement, its efficiency is somewhat limited for decreasing error in a QoI as the error contributions may be both positive or negative, and hence there is often significant cancellation of error [19]. This is in contrast to reducing error in standard norms which are always positive [40, 5, 12]. In classical adjoint based adaptivity, the absolute value of the elemental error indicators is taken and the principle of equidistribution applied. By taking the absolute value of the elemental error contributions, the cancellation of error due to opposing signs is lost. This phenomenon, along with a novel refinement strategy based on “mesoscale” regions is illustrated for ODEs in [19]. On the other hand, uniform refinement reduces the error predictably in the asymptotic regime and hence it is expected to reduce both the positive and negative elemental contributions equally. Thus, uniform refinement is expected to preserve the cancellation of error and this was observed experimentally in [19]. Uniform refinement is also more predictable in the expected decrease of error. In this article, we outline refinement strategies targeting sources of error as well as those based on elemental error indicators.

The main idea behind targeting sources of error to obtain accurate solutions is to reduce the dominant (in magnitude) source of error $E^{r},E^{m},E^{\Gamma}$ and $E^{\partial\Omega}$. This is accomplished by refining (either uniformly or adaptively) the corresponding discretization as shown in Table 1.

4.2.2 Uniform Contribution Refinement

In the Uniform Contribution Refinement, we choose the dominant component for refinement if it is at least 3 times larger than the next dominant component, or if both the top two dominant components have the same sign, so that the cancellation of error is preserved. If this requirement is not satisfied, the scheme defaults to standard uniform refinement.

4.2.3 Adaptive Contribution Refinement

The Adaptive Contribution Refinement is similar to the standard algorithms for adjoint weighted adaptive algorithms [33, 8, 2]. E.g., if the aim is to reduce the component $E^{r}$, then we define an elemental error indicator based on (37) as

[TABLE]

where $T\in\mathcal{T}$, and the subscripts ${T,\Gamma}$, $T$ and $T,m$ refer to evaluations of the integrals restricted to the element $T$ such that $T\cap\Gamma\neq\phi$, $T\in\Omega$ and $T\cap\Omega_{m}\neq\phi$ respectively. Once a per elemental error estimator is defined, the Dörfler scheme is used for marking the elements for refinement [27]. To preserve the cancellation of errors between different sources of error, all sources which have a total error contribution of at least half the dominant error contribution are selected to be adaptively refined.

4.2.4 Classical Refinement

In the classical adaptive refinement strategy we add up the terms in $E^{r},E^{m},E^{\Gamma},E^{\partial\Omega}$ and $E^{\mathrm{neg}}$ in Theorem 2 so that the error in the QoI for the nonlinear PBE is, $Q(u^{h}-U^{h},u^{r}-U^{r})\equiv E$ is,

[TABLE]

We define projection operators, $\pi_{m}:H^{1}_{0}(\Omega_{m})\rightarrow V_{h}^{m}$ and $\pi_{r}:H^{1}_{0}(\Omega)\rightarrow V_{h}$. From (21) we have,

[TABLE]

Combining (43) with (44) leads to the following elemental error indicator for element $T$

[TABLE]

The elemental error indicator for the linearized PBE is similar except that $\sinh({U}^{r})$ is replaced by $U^{r}$ and $\overline{\alpha}$ by $1$.

5 Numerical experiments

We show the accuracy of the a posteriori error estimates and utilization of the different sources of error to obtain an accurate computation of the QoI for the Born ion and methanol. The values of the constants in the PBE are chosen as $\epsilon_{m}=1$, $\epsilon_{s}=78$ and $\bar{\kappa}^{2}=0.918168$ unless otherwise stated. The value $\bar{\kappa}^{2}=0.918168$ corresponds to an ionic concentration of 0.1 M. These values reflect typical scenarios for PBE simulations [12, 18]. The initial meshes, defining the domains $\Omega_{m}$, $\Omega_{s}$ and the interface $\Gamma$ are generated using GAMer[66]. We use the standard space of continuous piecewise linear polynomials for the solution spaces corresponding to $u^{m}$ and $u^{r}$, that is for spaces $V_{h}^{m}$ and $V_{h}$. The spaces for the adjoint solutions $W^{h}$ and $W^{r}$ are chosen to be continuous piecewise quadratic polynomials. For ease of implementation, we always ensure that $\mathcal{T}_{m}=\mathcal{T}\cap\Omega_{m}$. The QoI (18) requires accurate integration near the points $x_{i}$. This is achieved by refining the cells near $x_{i}$ a few times. The functions $u^{h}_{d},u^{r}_{d},U^{h}_{d},U^{r}_{d},$ are such that they satisfy the boundary condition exactly on a boundary vertex and are zero on the interior vertices. This choice results in the component $E^{\mathrm{neg}}$ being exactly zero which was also verified numerically. Experiments are performed for the Born ion and the methanol molecule. The reference solutions needed for the effectivity ratios are computed using a mesh with 411635 vertices for the Born ion and a mesh with 90264 vertices for the methanol molecule and using continuous piecewise quadratic polynomials for the finite element space. The reference values of the QoI for the linearized PBE for the Born ion and methanol are -276.749875 and -48.477443 respectively. The corresponding values for the nonlinear PBE are -276.825527 and -48.479878. Since the reference solutions themselves have some error, effectivity ratios are only shown for experiments for which the reference solution is relatively accurate. All computations are carried out in the finite element software package DOLFIN from the FEniCS library [48, 49, 3, 47]. The value of $\eta$ for the QoI in (18) was chosen as 0.005. The Dörfler marking parameter is chosen as 0.2. The projection operators $\pi_{m}$ and $\pi_{r}$ were chosen as $L_{2}$ projectors.

5.1 Born Ion

The Born ion consists of a single point charge $Q_{1}$ in the center of a spherical solute domain $\Omega_{m}$ of radius $R$ [45]. The solute is surrounded by a large spherical solvent domain, $\Omega_{s}$, as depicted in Fig. 2a which has been adopted from [18].

Table 2 shows the error estimate, the effectivity ratio and different sources of error for the linearized PBE for two different meshes: an initial mesh of 6718 vertices and a uniformly refined mesh of 52014 vertices. In both cases, the effectivity ratio is close to one, indicating the accuracy of the error estimate. Moreover, we see that uniform refinement decreases all sources of error while preserving their signs, and hence accounts for cancellation of error. Similar results for the nonlinear PBE are shown in Table 3.

5.2 Methanol

We examine the accuracy of the error estimates in the more challenging setting of a methanol molecule, obtained from the APBS software package [6]. The methanol molecule consists of three charged particles representing charge groups: $\text{CH}_{3}$ and H with positive charges of 0.27 and 0.43 respectively, and an O atom with a negative charge of 0.7. The model is depicted in Fig. 2b adopted from [17]. The numerical experiments are performed on two meshes: an initial mesh of 11769 vertices and a uniformly refined mesh of 90264 vertices. The results for the linearized and nonlinear PBE are shown in Tables 4 and 5. The effectivity ratios are again close to $1.0$ and highlight the accuracy and robustness of the error estimate for both cases.

5.3 Refinement strategies

We use the different sources of error identified by $E^{r},E^{m},E^{\Gamma}$ and $E^{\partial\Omega}$ to guide refinement decisions. We first give an example of the effect of refining different discretization components, highlighting the significance of cancellation of error. Finally we present examples based on the Uniform Contribution Refinement, Adaptive Contribution Refinement and Classical Refinement schemes explained in §4.

5.3.1 Effect of refinement decisions on the QoI error

Consider the error information in Table 2 for the coarse mesh with 6718 vertices and error of $-1.14$. The uniformly refined mesh had 52014 vertices and an error of $-0.248$. On examining the different sources of error, we observe that the dominant error contribution is represented by $E^{\Gamma}$. Thus, instead of uniformly refining the mesh, we only refine simplices on the interface $\Gamma$. This refinement strategy is carried out by marking simplices in $\Omega_{m}$ which have one face on the interface, that is, marking simplices $T\in\mathcal{T}$ (and in $\mathcal{T}_{m}$) such that $T\in\Omega_{m}$ and $T\cap\Gamma$ is not the empty set.

The refinement results are shown in Table 6. The “It.” indicates the refinement level or iteration, with [math] indicating the starting coarse mesh. After the interface is refined, we arrive at level 1, corresponding to row having $N=15541$ vertices. Comparing this to the solution obtained by uniform refinement in Table 2, the error is now slightly less while the number of vertices is only 30% of the number of vertices of the uniformly refined mesh. This reflects a significant cost savings in obtaining accurate solutions.

At refinement level 1, once again, $E^{\Gamma}$ is the dominant component and we again refine the interface to reduce the error to $0.0124$. Now, both $E^{r}$ and $E^{\Gamma}$ have the same order of magnitude, but opposite signs. If we still carry on refining the interface to arrive at level 3-$\Gamma$. However, now the error has increased to $0.0621$! This behavior is quite common in numerical simulations, where refining a discretization parameter leads to an increase in the error rather than a reduction. Without the aid of adjoint based estimates, the cause of this increase may be hard to diagnose. The error information at levels 2 and 3 indicate why this increase occurred. The error at level 2 involved cancellation between the terms $E^{r}$ and $E^{\Gamma}$. Refining the interface significantly reduced $E^{\Gamma}$, while having only a marginal effect on $E^{r}$. Thus, there is less cancellation of error and the error increased to $0.0621$. A better option here is uniform refinement, which preserves the cancellation of error between different contributions [19]. The results of applying uniform refinement to Level 2 are shown as level 3-Uniform. The cancellation of error is preserved and the error decreased. Note that both the contribution of $E^{\Gamma}$ is the same for both levels $3(\Gamma)$ and $3-$Uniform, while the contribution of $E^{r}$ only sees a significant decrease at level $3-$Uniform.

5.3.2 Results for Uniform Contribution Refinement

The results of the Uniform Contribution Refinement strategy defined in §4 for the solution of the linear and nonlinear PBE for the Born ion are shown in Tables 7 and  8, while the results for methanol are shown in Tables 9 and 10. Comparing these results to Tables 2, 3, 4 and 5, we observe that the Uniform Contribution Refinement achieves significantly more accurate solutions with a lower computational cost (as measured by the number of vertices in the mesh) for both Born ion and methanol. An interesting observation is at level 2 of Table 8 where $E^{r}$ and $E^{\Gamma}$ have almost the same magnitude but opposite signs. These two sources of error cancel, leading to an unexpectedly low error.

5.3.3 Results for Adaptive Contribution Refinement

The results of the Adaptive Contribution Refinement strategy defined in §4 for the solution of the linear and nonlinear PBE for the Born ion are shown in Tables 11 and 12, while the results for methanol are shown in Tables 13 and 14. Comparing these results to Tables 2, 3, 4 and 5, we observe that the Adaptive Contribution Refinement is almost an order of magnitude more accurate for a uniformly refined mesh having the same number of vertices. Adaptive Contribution Refinement also outperforms the Uniform Contribution Refinement strategy for relatively small values of $N$. A couple of interesting observations are in order. In Table 12 the error decreases up to level 6 after which loss of error cancellation leads to an increase on level 7. In Table 14 we have an unexpectedly low error due to the cancellation between the terms $E^{r}$ and $E^{\Gamma}$.

5.3.4 Results for Classical Refinement

The results of the Classical Refinement strategy defined in §4 for the solution of the linearized and nonlinear PBE for the Born ion are shown in Tables 15 and 16, while the results for methanol are shown in Tables 17 and 18. Comparing these results to Tables 2, 3, 4 and 5, we observe that Classical refinement also performs well compared to uniform refinement. However, its performance is slightly worse than Adaptive Contribution Refinement as illustrated by Tables 11 and 15. In fact, the error at level 7 in table 15 shows almost a doubling of error at level 6. This is explained by observing the behavior of the terms $E^{r}$ and $E^{\Gamma}$, which are the two dominant sources of error, at these levels. Although both terms decrease in magnitude, there is less cancellation of error, leading to an overall increase. A similar increase in the error is observed at level 7 of Table 16.

5.3.5 Experiment illustrating difference between linear and nonlinear PBE results

In this section we perform an experiment to illustrate the difference in the results of the linear and nonlinear PBE solutions. To this end, we again choose the Born ion but now the charge on the ion, $Q_{1}$, is taken to be ten times its value in earlier experiments and also set $\bar{\kappa}^{2}=9.18168$ which is also ten times larger than in earlier experiments. We call this setup the highly charged Born ion. The difference in the computed QoI between the linear and nonlinear PBE for a mesh of 6718 vertices was approximately 59 units. The results for the different adaptive strategies also indicate different behavior between the linearized and nonlinear PBE.

The results for the linear and nonlinear PBE using Uniform Contribution refinement, Adaptive Contribution refinement and Classical Refinement are shown in Tables 19, 20, 21, 22, 23, 24. The results indicate that the Adaptive Contribution Refinement performs better than Classical Refinement for the linearized PBE, while they both perform equally well for the nonlinear PBE. Uniform Contribution Refinement outperforms both Classical Refinement and Adaptive Contribution Refinement.

6 Conclusions

Computing a QoI from the numerical solution of the PBE often has significant error that needs to be quantified. In this article, we develop adjoint based error estimates for this purpose. The adjoint operators are defined by accounting for the coupled nature of the three term split PBE as well as the issues arising due to the regularity of the normal derivative. The resulting error estimates are shown to be accurate, with effectivity ratios close to one. The error is partitioned in such a way that specific sources of error are identified and addressed. Moreover, novel refinement schemes, called Uniform Contribution Refinement and Adaptive Contribution Refinement in this article, utilize the information about the sources of error to arrive at accurate computed values of the QoI.

The effects of interface geometry on the error is an interesting area of future research. The current article is based on the the standard assumption in the PBE literature that the tessellated geometric representation of the interface is the true interface, e.g. as in references [2, 5, 20]. The effect of the geometry, which could be considered a “model form error”, is an interesting topic to explore and the author intends to pursue it in future.