{'fraction_non_alphanumeric': 0.054499608397486105, 'fraction_numerical': 0.018605895050526275, 'mean_word_length': 3.887042569081404, 'pattern_counts': {'":': 0, '<': 12, '<?xml version=': 0, '>': 16, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 15, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A new approach to constrained total variation solvation models and the study of solute-solvent interface profiles

In the past decade, variational implicit solvation models (VISM) have achieved great success in solvation energy predictions. However, all existing VISMs in literature lack the uniqueness of an energy minimizing solute-solvent interface and thus prevent us from studying many important properties of the interface profile. To overcome this difficulty, we introduce a new constrained VISM and conduct a rigorous analysis of the model. Existence, uniqueness and regularity of the energy minimizing interface has been studied. A necessary condition for the formation of a sharp solute-solvent interface has been derived. Moreover, we develop a novel approach to the variational analysis of the constrained model, which provides a complete answer to a question in our previous work [55]. Model validation and numerical implementation have been demonstrated by using several common biomolecular modeling tasks. Numerical simulations show that the solvation energies calculated from our new model match the experimental data very well.


Introduction
The description of the complex interactions between the solute and solvent plays an important role in essentially all chemical and biological processes. Solute-solvent interactions are typically described by solvation energies (or closely related quantities): the free energy of transferring the solute (e.g. macromolecules including proteins, DNA, RNA) from the vacuum to a solvent environment of interest (e.g. water at a certain ionic strength). There are two major approaches for solvation energy analysis, i.e., explicit solvent models and implicit solvent models [47]. Explicit models, treating solvent as individual molecules, are too computationally expensive for large solute-solvent systems, such as the solvation of macromolecules in ionic environments; in contrast, implicit models, by averaging the effect of solvent phase as continuum media [5,6,9,10,15,31,46], are much more efficient and thus are able to handle much larger systems [6,20,32,36,37,40,49,61].
Central in the description of the solvation energy in implicit solvent models is an interface separating the discrete solute and the continuum solvent domains. All of the physical properties of interest, including electrostatic free energies, biomolecular surface areas, molecular cavitation volumes and pK a values are very sensitive to the interface definition [26,59,63]. Variational implicit solvation models (VISM) stand out as a successful approach to compute the disposition of an interface separating the solute and the solvent [8,16,17,21,22,28,28,65,71]. In a VISM, the desired interface profile is obtained by minimizing a solvation energy functional coupling the discrete description of solute and the continuum description of solvent.
Despite of their initial successes in solvation energy calculations, sharp solute-solvent interface models suffer from several drawbacks. Firstly, from a physical point of view, there should be a smooth transition region, in which atoms of solute and solvent are mixed. In principle, an isolated molecule can be analyzed by the first principle -a quantum mechanical description of the wave function or density distribution of all the electrons and nuclei. However, such a description is computationally intractable for large biomolecules. Under physiological conditions, biomolecules are in a non-isolated environment, and are interacting with solvent molecules and/or other biomolecules. Therefore, their wave functions overlap spatially, so do their electron density distributions. Secondly, from an analytic point of view, the presence of geometric singularities is inevitable in many conventional VISMs. It makes the underlying model lack stability and differentiability, which generates an intrinsic difficulty in the rigorous analysis of the model. Thirdly, from a computational point of view, these surface configurations produce fundamental difficulty in the simulation of the governing partial differential equations (PDEs), like the Poisson-Boltzmann (PB) equation. Those considerations motivate the use of the diffuse solvent-solute interface definition.
Among all effort to ameliorate the solvent-solute interface definition, arguably, one of the most extensively used models is the total variation based model (TVBVISM), cf. [17,27,64,[66][67][68]. The main idea of TVBVISM is based on a transition parameter u : Ω → [0, 1] such that u takes value 1 in the solute and 0 in the solvent region. More precisely, the following total solvation free energy was proposed in terms of u: Here the constant γ > 0 is the surface tension. By the coarea formula for a Lipschitz function u : Ω → [0, 1], where H 2 stands for the 2-dimensional Hausdorff measure. Hence, the total variation term Du (Ω) represents the mean surface area of a family of isosurfaces Ω ∩ u −1 (t). See [66] for more detail. According to this geometric interpretation, γ Du (Ω), measures the disruption of intermolecular and/or intramolecular bonds during the solvation process.
The constant P h is the hydrodynamic pressure. In a previous work [55], we proposed a novel physical interpretation of the characteristic function u so that u(x) represents the volume ratio of the solute at x ∈ Ω. Therefore, Ω P h u dx is the mechanical work of creating the biomolecular size vacuum in the solvent. ρ s is the constant solvent bulk density, and U vdW (x) is the attractive portion of the Van der Waals potential at point x. It represents the attractive dispersion effects near the solute-solvent interface and has been shown by Wagoner and Baker [63] to play a crucial role in accurate nonpolar solvation analysis. The first three terms are usually termed the nonpolar portion of the solvation free energy.
The second and third lines of (1) are usually called the polar portion of the solvation free energy, in which ψ is the electrostatic potential. ρ m is an L ∞ -approximation of the density of molecular charges; m and s are the dielectric constants of the solute molecule and the solvent, respectively, with 0 < m s . q j is the charge of ion species j = 1, 2, · · · , N c ; and c ∞ j is the bulk concentration of the j-th ionic species. Finally, β = 1/k B T , where k B is the Boltzmann constant and T is the absolute temperature. For notational brevity, throughout this paper, we put (1.2) Numerical simulations show that diffuse-interface models can significantly improve the accuracy and efficiency of solvation energy computation [8,16,17,21,22,28,28,45,65,71]. In contrast, on a theoretic level, there are several open questions concerning model (1).
First, the uniqueness of a minimizer is unknown for (1). Indeed, most of the solvation energy functionals, regardless of sharp or diffuse interfaces, only predict local minimizers, cf. [8,16,17,21,22,28,45,65,71]. As a consequence, solutions of the corresponding Euler-Lagrange equations may not correctly depict the energy minimizing interface profile. In contrast, any minimizer of (1) is global. However, lacking strict convexity, (1) may admit multiple global minimizers. This prevents us from studying many properties of the interface profile, e.g. the size of the set of discontinuities. These observations motivate us to introduce strict convexity into model (1) by including a new parameter p = 2N 2N −1 with N ∈ N so that u p (x) represents the volume ratio of the solute at x ∈ Ω. It is important to notice that the geometric meaning of the term Du (Ω) remains the same as in the original model (1). We will establish the existence, uniqueness and regularity of the global minimizer of the modified model, see (2.3).
Second, the natural admissible space to minimize (1) is the space of BV −functions. Therefore, it is possible that model (1) is minimized by the characteristic function of a set of finite perimete. This corresponds to a sharp solute-solvent interface, an unrealistic situation as discussed before. Nevertheless, it is mathematically impossible to exclude such situations in model (1) due to the lack of uniqueness of a minimizer. Based on the modified model, this work provides a partial answer to the question why the solvation free energy is not minimized by a sharp interface. More precisely, we show that a necessary condition for a nonpolar molecule to have a sharp energy-minimizing interface is that the mean curvature of its Van Der Waals surface is everywhere nonpositive, which is unrealistic for almost all real-world biomolecules. To the best of our knowledge, our work is the first to give a mathematical explanation of such phenomenon.
Third, the physical meaning of the characteristic function u enforces two biological constraints: (1) u needs to be 1 for the pure solute region and 0 in the pure solvent area, and (2) as a volume ratio function, it must satisfy that 0 ≤ u ≤ 1. This leads to a constrained total variation model (2.3), which is a non-differentiable functional with a two-sided obstacle. It is known that the Euler-Lagrange equations of similar functionals with simpler structure and without obstacle, e.g. Rudin-Osher-Fatemi models, were formally derived by using the 1−Laplacian operator [54]. With the presence of the obstacle, on a heuristic level with sufficiently smooth minimizer u and energy functional, one expects the corresponding first variations with respect to u to take the form of a variational inequality, or equivalently, of a 1−Laplacian type equation involving a measure supported on the coincidence sets {u = 0} and {u = 1}. Unfortunately, both the functional (2.3) and the minimizer u lack the required smoothness. This casts a shadow over the study of the first variations of the constrained total variation model, not even formally. In [55], we proposed a novel approach to the variational analysis of such constrained VISM via approximation by a sequence of q-energy type functionals. This approach was applied to the numerical study of the nonpolar energy in our previous work [55]. Using a similar idea and the new volume ratio function u p , we will rigorously derive the variational formulas of the new total energy functional.
The rest of the paper is organized as follows. A list of the main theorems is stated at the end of the introduction. In Section 2, we state the precise definition of our new model. In Section 3, we study a family of perturbed Poisson-Boltzmann equations. These equations will be used in Sections 4 and 6. Section 4 is devoted to the validation of the model, in which we prove the existence and uniqueness of a minimizer and the continuous dependence of the solvation energy on the biological constraints. In Section 5, a necessary condition for the formation of a sharp solute-solvent interface is derived. The argument heavily relies on the tools from nonsmooth convex analysis. In Section 6, we conduct a variational analysis of our new model by means of an approximation argument. Base on this analysis, our model, including its solvation energy and solute-solvent interface predictions, is studied through numerical simulations. For the readers' convenience, we include two appendices at the end of this article, one on BV −functions and the other on nonsmooth convex analysis.
For the reader's convenience, we will give a list of the main theoretic results here: • Theorem 4.1: the existence and uniqueness of global minimizers of the total solvation energy; • Theorem 4.2: the continuous dependence of the solvation energy on the biological constraints; • Theorem 5.10: a necessary condition for the formation of a sharp solute-solvent interface; • Theorem 6.3: the theoretic basis of the numerical simulations.

Solvation Free Energy Functional
2.1. Notations. In this article, we use x = (x 1 , x 2 , · · · , x N ) to denote the coordinates in R N . S N −1 denotes the (N − 1)−sphere in R N . Given two vectors u, v ∈ R N , u · v is their inner products.
Given U ⊆ R N , U stands for the closure of U . The topological boundary of U is denoted by ∂U . Given two domains U and Ω in R N , U ⊂⊂ Ω means that U ⊂ Ω.
For any two Banach spaces X, Y , the notation X → Y means that X is continuously embedded in Y . Given a sequence {u k } ∞ k=1 = (u 1 , u 2 , · · · ) in X, u k u in X means that u k converge weakly to some u ∈ X.
Given 1 ≤ p ≤ ∞, let p be its Hölder conjugate. L p (U, X) is the set of all X-valued p−integrable (Lebesgue) measurable functions defined on U , whose norm is denoted by · p . The notation X is sometimes omitted when its choice is clear from the context. W k,p (U ) stands for the Sobolev space consisting of functions whose weak derivatives up to k−th power belong to L p (U ). Additionally, H 1 (U ) = W 1,2 (U ).
Given two sets A and B, A ⊆ B and A ⊂ B mean that A is a subset and a proper subset of B, respectively. Finally, we denote by L N and H N −1 the N −dimensional Lebesgue measure and the (N − 1)−dimensional Hausdorff measure, respectively.

An Experimental Based Domain Decomposition.
Let Ω ⊆ R 3 be a bounded and connected Lipschitz domain composed of three disjoint subdomains: • Ω m : solute (molecular) region; • Ω s : solvent region; • Ω t : solute-solvent mixing region.
We further assume that ∂Ω ⊂ ∂Ω s and ∂Ω m ⊂ ∂Ω t . Let Σ 1 = ∂Ω m be a smoothed Van Der Waals surface enclosing the pure solute region and be the smoothed solvent accessible surface outside which is the pure solvent domain. Suppose that Σ 1 ∩Σ 0 = ∅ and Ω m , Ω s are non-empty. In addition, we assume that Σ i , i = 0, 1, are embedded closed Lipschitz surfaces. In this article, a closed surface always means one that is compact, without boundary and embedded in R 3 . We further assume that the solute region Ω m contains N a solute atoms located at x 1 , · · · , x Na ; and there are N c ion species outside Ω m . Finally, for notational brevity, we put Ω w = Ω \ Ω s . A picture illustration of the domain definition and decomposition can be found in Figure 3

A Novel Solvation Energy Functional.
As an improvement of the previous differential geometric based solvation model [17,55], we study a novel solvation free energy, whose nonpolar portion is defined as 2N −1 for some integer N > 1 and λ, P h > 0. Note that 1 < p < 3 2 . Since 3 2 = 1 * is the Sobolev dual of 1, we have BV (Ω) → L p (Ω). Here u : Ω → R represents a characterizing function of the solute such that u p (x) is the volume ratio at position x ∈ Ω (as shown in Figure 3). As such, the physical constraints u(x) ∈ [0, 1] for a.a. x ∈ Ω (2.1) and u = 1 a.e. in Ω m and u = 0 a.e. in Ω s (2.2) need to be imposed. Note that U vdW (x) can be formulated by i U att i (x) in which U att i (x) represents the attractive part of Lennard-Jones potential [17,63]. To this end, the L-J potential can be divided into attractive U att i and repulsive U rep i in different ways. Here we take a Weeks-Chandler-Andersen (WCA) decomposition based on the original WCA theory [42]: with parameters is of energy and σ is of length.
We choose Ω m in such a way that there exist balls B(x i , τ ) with i = 1, · · · , N a and τ > 0 such that The polar portion of the solvation free energy is defined as Here (u) = u p m + (1 − u p ) s is the dielectric constant of the solvent/solute mixture. ρ m is supported in Ω m . In addition, the neutral condition holds for some ψ ∞ ∈ W 1,∞ (Ω). Therefore given u ∈ BV (Ω) satisfying (2.3), ψ = ψ(u) is determined via the elliptic boundary value problem (2.3).
With the above observations, the minimization problem can be restated as to minimize in the admissible space

A Family of Perturbed Poisson-Boltzmann Equation
In this section, we study a sequence of functionals associated with the polar free energy, which will be used in the numerical simulations in Section 6.
Proposition 3.1. Given any u ∈ Y k , k = 0, 1, · · · , there exists a unique ψ u ∈ A such that Moreover, ψ u is the unique weak solution to (3). Further, ψ u satisfies In particular, the constant C 0 is independent of Ω m , Ω s , u and k.
Proof. Analogous problems have been studied in the literature on various Poisson-Boltzmann type equations, cf. [22,44,45,55]. In order to show the determining factors of the constant C 0 in (3.1), we will, nevertheless, state a brief proof.
The constant M 0 depends only on Ω, s , m , q 1 and ψ ∞ . DefineG k u : By the direct method of calculus of variation and the strict convexity ofG k u (·), there exists a global minimizer Let ψ u =ψ u +ψ u . From the above equality, we learn that ψ u minimizes G k u (·) in Y k . Then following Steps (iii) and (iv) in the proof of [55, Proposition 2.2], we can show that ψ u ∞ + ψ u H 1 ≤ M 1 for some constant M 1 depending only on M 0 . We can take C 0 = M 0 + M 1 .
The above proposition immediately gives the following crucial estimates. For every k and u ∈ Y k , where ψ u is the solution to (3). The constant C 1 is independent of Ω m , Ω s , k and the choice of u.
Proposition 3.2. Let u k ∈ Y k , k = 0, 1, · · · , be such that If, in addition, u k ∈ Y and ψ k ∈ A satisfies G 0 u k ( ψ k ) = min w∈A G 0 u k (w). Then Proof. We will only prove (3.2). The proof for (3.2) is similar.

This establishes (3). It follows from the Poincaré inequality that
then can be shown by using (3) and the dominated convergence theorem.

Properties of Global Minimizers
The following theorem on the existence and uniqueness of a minimizer of I(·) can be proved essentially in the same way as [  To show the robustness of the model (2.3), one need to answer the question whether the solvation energy I(u min ) depends continuously on Ω m and Ω s in a suitable topology? The answer to the above question is affirmative. We will present the proof of a partial result in this subsection. Due to the length of this article, a complete answer will be presented in a subsequent paper.
Assume that { Ω m;n } ∞ n=1 and { Ω s;n } ∞ n=1 are two sequences of Lipschitz subdomains such that We consider the sequence of energy functionals I n (·) defined by replacing Ω m and Ω s by Ω m;n and Ω s;n in I(·), respectively. The corresponding admissible spaces are in Ω and u = 1 a.e. in Ω m;n and u = 0 a.e. in Ω s;n }.
Then for each n, there is a unique minimizer u n of I n (·) in Y n . Moreover, lim n→∞ I n (u n ) = I(u min ).
Proof. The existence and uniqueness of a minimizer of I n (·) in Y n for each n follows from Theorem 4.1.
Observe that u min ∈ Y n for all n. Thus This implies that where C 1 is the constant in (3). Therefore, u n BV is uniformly bounded with respect to n. Proposition A.2 implies that there exists a subsequence, not relabelled, and some u ∈ BV (Ω) such that u n → u in L 1 (Ω). From Propositions A.3, Propositions 3.2 and the dominated convergence theorem, we infer that This proves the convergence assertion.
where Per(Ω m ; Ω) is the perimeter of Ω m in Ω, see Appendix A, and G ele (Ω m ) is the electrostatic free energy.
In the classic Poisson-Boltzmann theory, it is defined by cf. [2,15,23,43,56,69,70]. The electrostatic potential ψ solves the classic sharp-interface Poisson-Boltzmann equation: The following corollary shows that (4) is in some sense the limiting case of our diffuse interface model. Remark 4.4. In a subsequent paper, we will show that, under mild regularity assumption on Σ 1 and Σ 0 , the conditions Ω m;n ⊆ Ω m and Ω s;n ⊆ Ω s in Theorem 4.2 can be relaxed.

How to Exclude the Formation of Sharp Interfaces?
In Theorem 4.1, we have shown that there is a unique characterizing function u min ∈ BV (Ω) minimizing (2.3) in Y. However, since BV −functions allow jump discontinuities, a natural question to ask is whether the minimizing energy state is achieved by a sharp interface between the solute and solvent regions, or equivalently, whether the characterizing function u min is the characteristic function of a set of finite perimeter.
To simplify the analysis, we will focus on the nonpolar portion of the solvation energy, i.e. (2.3). Motived by the idea in [12][13][14], we will show that when the mean curvature of Σ 0 is positive at some point, the energy minimizing state is never achieved by a sharp interface. See Theorem 5.10.

Necessary Conditions for the Minimizer of Nonpolar Energy.
Throughout this section, we assume that Ω t = ∅. First consider the minimization problem of the nonpolar energy 3)}. One will show that the minimizer u min of (5.1) automatically satisfies Constraint (2.3). The reason to exclude (2.3) in the definition of the admissible space is due to the following consideration. Any subdifferential of I np (·) with Constraint (2.3) contains a function which may be discontinuous along ∂{u min = 1} and ∂{u min = 0}. This will prevent us from establishing the continuity of u min in these two sets.
Proof. Note that X is closed and convex in BV (Ω). Based on the strict convexity, lower semicontinuity of I np and the direct method of Calculus of Variation, we can readily establish the existence and uniqueness of a global minimizer u min . If Direct computations show that I np ( u min ) < I np (u min ). A contradiction. Therefore, 0 ≤ u min ≤ 1 a.e. in Ω.
Next, we derive necessary conditions for the minimizer of (5.1). We will use tools from non-smooth analysis, c.f. [24,25,29], to derive the subdifferential of (5.1). However, very little is known about the dual space of BV (Ω). To overcome this difficulty and tackle the Constraint (2.3), we will consider I np as a functional defined on L p (Ω) and include two extra terms. Define where Tru is the trace of u on ∂Ω and K = {u ∈ L p (Ω) : u = 1 in Ω m , and u = 0 in Ω s a.e.} and I K is the indicator function of K. In addition, we put and The latter is Lipschitz continuous in L p (Ω). It is understood that So, dom(E 1 ) = BV (Ω) and dom(I K ) = K. Using these notations, we can restate Problem (5.1) as to minimize a functional E np : Direct computations show that u min minimizes (5.1) in X iff it minimizes E np (·) in L p (Ω).
Note that K is closed and convex in L p (Ω). This implies that I K is convex and lower semicontinuous. What is more, by the definition of subdifferentials, for every u ∈ K, u * ∈ ∂I K (u) iff Here ·, · is the duality pairing between L p (Ω) and Then in Ω t . This is also the sufficient condition of u * ∈ ∂I K (u). Indeed, given any u * ∈ L p (Ω) with u * = 0 a.e. in Ω t , for any v ∈ K, To sum up, a function u * ∈ L p (Ω) belongs to ∂I K (u) iff u * = 0 in Ω t .
To compute ∂E 1 (u), we define . Given any u ∈ BV (Ω) and z ∈ X ∞ p , there exists a Radon measure, denoted by z · Du, such that for any φ ∈ C ∞ 0 (Ω), with a little abuse of notation, The measure z · Du is absolutely continuous with respect to |Du|. By the Radon-Nikodym Theorem, there is a |Du|-measurable function θ(z, Du) s.t.
Next, Proposition B.1 implies that for any u ∈ L p (Ω), Because of the lack of continuity of E 1 and I K , in general, we can only conclude that ∂E 1 (u) + ∂I K (u) ⊆ ∂(E 1 + I K )(u). In order to compute ∂(E 1 + I K )(u), we will use Propositions B.3. It suffices to verify the closed linear space condition. An easy computation shows that which is obviously a linear subspace of L p (Ω). We learn from Propositions A.3 and A.6 that dom(E 1 ) − dom(I K ) is closed. Now Proposition B.3 immediately implies that From the definition of subdifferential and (5.1), we learn that More precisely, this means that there is some z ∈ X ∞ p with z ∞ = 1 satisfying (5.1) and w ∈ L p (Ω) with In particular, it holds that

5.2.
Regularity of the Minimizer u min . As in the previous subsection, u min is the minimizer of (5.1) in where the minimum is taken in the set Proof. Take z as in (5.1 Given any φ ∈ C 0 (Ω), by [3, Proposition 2.7(i) and Formula (2.15)], we have On the other hand, by the coarea formula (A), holds for all F ∈ M. Combining with (5.1), we thus deduce that Therefore, (5.2) holds for t. We thus deduce that the assertion holds for any t ∈ [0, 1).
Remark 5.3. The existence of a minimizer of (5.2) can be proved by using the classical method of Calculus of Variation for every t ∈ [0, 1). See [38] for a related problem.
Lemma 5.4. Let t < t. If F t and F t are minimizers of (5.2) with t and t , respectively, then F t ⊆ F t .
Proof. We clearly have But t < t. This implies that F t ⊆ F t .
It follows from Proposition 5.2 and Lemma 5.4 that However, there are only countably many t such that L 3 ({u = t}) > 0. This implies that This completes the proof.
Proof. For any x ∈ ∂E t ∩ Ω t , for sufficiently small r > 0, the ball B(x, r) is contained in Ω t . For any local perturbation of E t in B(x, r), i.e. a set F of finite perimeter such that by Hölder inequality for any δ ∈ (0, 1). Note that the constant C in the above inequality is independent of the position of x. Hence E t ∩ Ω t is almost minimal in Ω t in the sense of [60, Definition 1.5]. Therefore, [60, Theorem 1.9] implies that the singular set of E t is contained in Σ 0 ∪Σ 1 and ∂E t \(Σ 0 ∪Σ 1 ) is a C 1 -hypersurface. Then the assertion follows from the standard regularity theorem of non-parametric minimizing surfaces, see [35] for example. For the reader's convenience, we will state a proof here. For every x 0 ∈ ∂E t \ (Σ 0 ∪ Σ 1 ), denote by H x0 the tangent plane of ∂E t at x 0 . Near x 0 , we can rewrite the coordinates in the form x = (y, z), where y is the coordinates in H and z is the coordinate in the normal direction of H. We use the convention z = y = 0 at x 0 . For sufficiently small r > 0, let U r = B(x, r)∩H x0 . Build a cylinder C r = U r ×(−r, r) ⊂⊂ Ω t in (y, z)−coordinates centered at x 0 . Inside C r , we can express ∂E t as the graph of a C 1 −function v: See Figure 2. Then By the fundamental theorem of Calculus, v solves By choosing r > 0 sufficiently small, one can infer from [34,Theorems 16.10] that v ∈ C 2 (U r ). The remaining regularity follows from a bootstrapping argument, cf. [34, Theorems 6.13 and 6.17].
Remark 5.7. If we assume, in addition, that Σ i ∈ C 1,1 for i = 0, 1, then following the argument in [60, Section 1.14(iv)], one can show that the singular set of E t is empty and ∂E t ∈ C 1,1 . Since this fact will not be used below, to keep the article in a reasonable length, we will not provide a rigorous proof here.
Proposition 5.8. The jump set, J umin , of u min is contained in Σ 0 ∪ Σ 1 .
Proof. The proof follows the idea in [13,Theorem 3.4]. By (A), it suffices to show that for any t 1 < t 2 ∈ [0, 1) and t 1 , t 2 ∈ Q, it holds . By Proposition 5.6, both ∂E t1 and ∂E t2 are regular in a neighbourhood of x 0 . From the fact E t2 ⊆ E t1 , we deduce that the tangent space of E t2 and E t1 at x 0 agree. Denote the tangent space by H x0 . We define the coordinates in the form x = (y, z) and the cylinder C r = (−r, r) × U r as in the previous proof. Then we can express E ti with i = 1, 2 as graphs over U r as Since t 2 > t 1 , v i (0) = 0, ∇ y v i (0) = 0, by choosing r > 0 small enough, we have for all y ∈ U r . This implies that In view of the boundary condition v 1 ≥ v 2 on ∂U r , we infer from [34, Theorem 10.1] that v 2 < v 1 in U r , which contradicts v 1 (x 0 ) = v 2 (x 0 ). Therefore, (∂E t1 ∩ ∂E t2 ) \ (Σ 0 ∪ Σ 1 ) = ∅. Remark 5.9. In particular, Proposition 5.8 implies that u ∈ C(Ω t ).

Necessary Conditions for the Formation of a Sharp Interface.
In this section, we first consider the case that Ω t is connected. In order to state the main theorem of this section, we define the orientations of Σ i in such a way that • the outer normal of Σ 1 points into Ω t , and • the outer normal of Σ 0 points into Ω s .
With these conventions, a sphere of radius R > 0 has constant mean curvature −1/R. Theorem 5.10. Suppose that Ω t is connected and Σ i , for i = 0, 1, are C 2 −closed surfaces. Let κ be the mean curvature of Σ 1 . If κ(p) > 0 for some p ∈ Σ 1 , then there is no sharp solute-solvent interface, that is, the minimizer u min of (2.3) is not the characteristic function of a set E of finite perimeter with Ω m ⊆ E ⊆ Ω w .
Proof. Assume, to the contrary, that there exists a set E of finite perimeter such that Ω m ⊆ E ⊆ Ω w and χ E minimizes (2.3).
(3) Let κ be the mean curvature of Σ 0 . If E = Ω w , then following the above argument, we conclude that is a necessary condition for E = Ω w However, it is well known that there is no closed hypersurface with everywhere positive mean curvature in R 3 . Therefore, E = Ω w Remark 5.11. The mean curvature condition κ(p) > 0 for some p ∈ Σ 1 is satisfied by almost all macromolecules. This explains why diffuse interfaces are indeed more realistic in real-world solvation processes. It is equally important to point out that the mean curvature condition is in some sense "stable". Recall that the Hausdorff metric on compact subsets K ⊂ R n , n ∈ N, is defined by Given a closed surface Σ in R 3 , its second normal bundle is given by where ∇ Σ is the surface gradient defined by for some r > 0. Here P Σ (p) = I − ν Σ (p) ⊗ ν Σ (p). Denote by M the set of all connected closed surfaces in R 3 . Equipped with the metric d H , M is a Banach manifold, cf. [51,52]. If a connected component, M 1 , of Σ 1 satisfies the condition in Theorem 5.10, then any Σ ∈ M that is sufficiently close to M 1 with respect to the metric d H satisfies the same condition.
Remark 5.12. The connectedness condition of Ω t was used in the proof of Theorem 5.10. It is well-known that cavities may appear inside macromolecules, which corresponds to the situation of disconnected Ω m . In the case of N cavities inside Ω m , Ω t consists of N + 1 connected components. More precisely, where Γ j are C 2 −closed and connected hypersurfaces and Γ j , j = 1, · · · , N , is the boundary of the j−th cavity. Correspondingly, where Ω t,j are the connected components of Ω t and Ω t,j , j = 1, · · · , N , is the j−th cavity bounded by Γ j and Ω t,0 ∩ Ω s = ∅. See Figure 3 for a picture illustration of a solute molecule with one cavity inside. To make the convention of the mean curvature consistent, we define the orientation of Γ j in the following way: • the outer normal of Γ 0 points into Ω b,0 ; • for j = 1, · · · , N , the outer normal of Γ j points into Ω m .
Under these conventions, we can follow the proof of Theorem 5.10 and show that Γ j (j = 1, · · · , N ) is a sharp interface iff Γ j has everywhere positive mean curvature, which is impossible. Therefore, none of the cavities can be purely occupied by the solvent.

Numerical Simulations
The non-differentiable structure of (2.3) and the Constraints (2.3) and (2.3) generate an essential difficulty in the numerical simulations of (2.3). This motives us to study a sequence of approximation problems.

An Approximation Problem.
Recall the definition of {q k } ∞ k=1 from Section 3. We introduce a family of perturbed solvation free energy functionals where ψ ∈ A satisfies (3). We will seek a minimizer of I k (·) in Y k , c.f. (3). For notational brevity, we term the second line of (6.1) I p,k (u, ψ).
Let u min be a minimizer of (2.3) in Y and ψ min = ψ umin be the solution of (2.3) with u = u min .
To prepare for the main result of this section, we introduce Ω j,k := {x ∈ Ω : dis(x, Ω j ) < 1/k}, j ∈ {m, s} and X k := {u ∈ Y : u ≡ 1 in Ω m,k and u ≡ 0 in Ω s,k }, k ∈ N, and quote the following two lemmas from [55].
The theoretic basis of the numerical simulation is the following theorem. Theorem 6.3. For each k = 1, 2, · · · , there exists a unique u min,k ∈ Y k ∩Y such that I k (u min,k ) = min If, in addition, and as k → ∞ u min,k → u min in L r (Ω) for all r ∈ [1, ∞) and ψ min,k → ψ min in H 1 (Ω), where ψ min,k = ψ u min,k is the solution to (3) with u = u min,k .
Proof. (i) The existence and uniqueness of a minimizer u min,k ∈ Y k of I k (·) for each k can be proved in the same manner as in Theorem 4.1.
(ii) We will show that u min,k is a global minimizer of I k (·) iff (u min,k , ψ min,k ) is a saddle point of Here C 0 is the constant in (3.1). Proposition 3.1 shows that ψ min,k ∈ A . Denote by S k the set of all saddle points of L k . Recall that (u 0 , ψ 0 ) ∈ S k iff Note that Y k and A are closed and convex in W 1,q k (Ω) and H 1 (Ω), respectively. Moreover, It follows from the direct method of Calculus of Variation that the infimum is achieved. Therefore, By [29, Proposition VI.1.2], (u min,k , ψ min,k ) ∈ S k . Conversely, if (u 0 , ψ 0 ) ∈ S k , then (6.1) and Proposition 3.1 show that ψ 0 is the solution of (3) with u = u 0 . What is more, since (6.1) still holds true if we replace (u min,k , ψ min,k ) by (u 0 , ψ 0 ), we infer that u 0 = u min,k .
A contradiction. Hence, u min,k ∈ Y.
where C 0 is the constant in Proposition 3.1 and C 2 is independent of k and v. This yields that where C 3 = Ω\Ωm ρ s U vdW dx. We thus infer from (6.1) that for some C 4 independent of k. Proposition A.2 implies that there exists a subsequence of {u min,k } ∞ k=1 , not relabelled, converging to some u 0 ∈ Y in L 1 (Ω). The Riesz-Thorin interpolation theorem then implies that

Then it follows from Propositions A.3 and 3.2 that
On the other hand, we define x ∈ Ω s,n u 0 (x), elsewhere.
We will show that lim sup k→∞ I k (u min,k ) ≤ I(w n ). (6.6) Lemma 6.1 implies that we can find a sequence {w n,i } ∞ i=1 such that w n,i ∈ C ∞ (Ω) ∩ Y k for all k and w n,i → w n in L 1 (Ω) and Dw n,i (Ω) → Dw n (Ω) as i → ∞.
Since u min,k minimizes I k (·) in Y k , we have Pushing k → ∞, the dominated convergence theorem and Proposition 3.  By the uniqueness of a global minimizer of I(·), we conclude that u 0 = u min .

Variation of Solvation Free Energy.
Motivated by Theorem 6.3, we will study the numerical simulations of the approximating functional (6.1). As the first step, we will derive the variational formulas of (6.1) at u min,k . Recall that u min,k minimizes (6.1) in Y k iff (u min,k , ψ min,k ) is a saddle point of (6.1) in Y k × A , where ψ min,k solves (3) with u = u min,k . This means that L k (u min,k , ψ min,k ) = min u∈Y k L k (u, ψ min,k ).
Therefore, we can verify that u min,k satisfies In view of (3), (u min,k , ψ min,k ) solves the following elliptic system Remark 6.4. The approach in this section actually gives a solution to the variational analysis of (1) with Constraints 2.3 and 2.3, which provides a complete answer to a question in our previous work [55].
6.3. Computational methods. This section presents the computational methods and algorithms for the solution of the coupled system (6.2) and its associated parameterization process. The solution of (6.2) provides a physically sound "diffuse solute-solvent interface profile" u and the electrostatic potential ψ, and thereby the calculation of the total solvation free energy.
While solving for u and ψ, the surface evolution equation and the perturbed PB equation cannot be decoupled and thus need to be solved simultaneously. In the following, we first describe in more detail about the solution methods for each equation and their discretized formulations. Then the scheme for the convergence of two coupled equations is presented as well as a simple parameterization approach for optimal parameter values.

The perturbed Poisson-Boltzmann equation.
For the solution of perturbed PB (PPB) equation, we adopted the finite difference scheme. Thanks to the continuous dielectric function, an accurate solution can be achieved with a standard second-order center difference scheme. Specifically, for a solvent without salt, the PPB equation can be simplified to a perturbed Poisson equation. If the position (x i , y j , z k ) is represented by the pixel (i, j, k), its discretized form becomes where the uniform grid spacing h is applied at x, y and z directions, and (i + 1 2 , j, k) = (u(x i + 1 2 h, y j , z k )), q(i, j, k) is used to describe the fractional charge at grid point (x i , y j , z k ). The fractional charge is calculated by the second-order interpolation (trilinear) of the charge density ρ m . Then a standard linear algebraic equation system is resulted from the the discretized perturbed Poisson equation in the form of AX = B, in which X is the targeted solution. Matrix A is the discretization matrix and B is the source term according to the discrete charges.
The boundary condition of PPB equation is computed via the summation of electrostatic potential contributions of individual atom charges [33]. The resulted linear system can be solved by various linear solvers (like biconjugate gradient in this study) together with pre-conditioners for potential acceleration. 0 can be used for the initial guess of the solution and convergence tolerance is set as a small number such as 10 −6 . It has been shown that the designed PB solver is capable of delivering second-order accuracy [17].
6.3.2. The surface evolution equation. The solution of the surface evolution equation can be attained via the following parabolic PDE as done in earlier work [7,17]. Computationally, the equation (6.3.2) can be expanded into a form as follows.
In particular, the time-dependent derivative is carried out by explicit Euler scheme. Note that other implicit schemes can be designed to improve the solution efficiency and will be pursued later. The first and second order spatial derivatives are handled by finite difference schemes [17]. To impose the domain decomposition in (6.2), we let u be fixed as one in the pure solute area Ω m and as zero in the pure solvent region Ω s . Here the pure solute area is numerically defined to be enclosed by a smoothed Van Der Waals surface (vdW) and the the pure solvent region is the area outside a smoothed solvent accessible surface (SAS). The initial value of u in between Ω m and Ω s can be set between 0 and 1.

6.3.3.
Coupling of the perturbed Poisson Boltzmann and surface evolution equations. In principle, the surface evolution equation needs to be solved simultaneously with the perturbed PB equation until the solution process reaches a self-consistency. To speed up the whole iterative process, electrostatic potential ψ is updated after a number of time steps (i.e., 10 to 100 steps) evolution of the parabolic surface equation [17].
Moreover, a simple relaxation algorithm is adopted to guarantee the convergence of the iterative process as follows [17] : where u new and u old are the new and old u profile values from current and previous steps, respectively. ψ old and ψ new denote previous and new electrostatic potentials, respectively. α = 0.5 and α = 0.5 are set in our calculation.
In addition, a simple cutoff strategy is conducted to apply Constraint (2.3) and to avoid possible numerical errors: u > 1. The cutoff checkup is carried out every time step or several steps during the solution of surface evolution equation.
Finally, to reduces the total iteration number and save the computational time significantly, first of all, one may start the iterative process with an initial u from solving Eq. (6.3.2) without the electrostatic potential term. Second, one may take the prior potential ψ as a good guess for the next resulted linear system in the PPB solution. That will make the PPB solver converge faster. 6.3.4. Parameterization. There are some parameter values that need to be determined for real numerical simulations of solvation free energy. They include solvent density ρ s , the solvent radius σ s ,γ, P h and so on. Since most of the parameters are involved in nonpolar solvation energy, a previous simple parameter fitting strategy is adopted here [19,55]. In particular, on the one side, some parameter values are fixed or given such as: ρ s =0.03341/Å 3 ; solvent radius σ s =0.65Å ; radii of solute atoms like σ c =1.87Å . On the other side, some are considered as fitting parameters like γ, P h , and well depth parameters is where i denotes different atom types. The following iterative procedure is used to obtain the optimal fitting parameter values: Step 0: An initial guess of fitting parameters and a trial set of molecules are determined with their existing informaion such as atomic coordinates, radii, and experimental data of solvation free energies.
Step 1: For individual j-th molecule, j = 1, · · · N m where N m is the total number of molecules in the trial set, the coupled system (6.2) is solved until self-consistency is reached to find the quasi-steady state solution of u j and ψ j with latest parameter values. Note that if the trial set is nonpolar, one only needs to solve the surface evolution equation without a driven potential from the electrostatic field. Then the fitting process is exactly the same as our previous paper [55].
Step 2: Electrostatic solvation energy I j p,q k is calculated for each molecules using the profile of ψ j .
Step 3: A non-negative least squares algorithm is used to update all non-negative parameters P h , γ, and is with a minimization problem where I j,exp q k is the existing experimental data of solvation free energies in the literature.
Step 4: The iterative loop from Step 1 to Step 3 is repeated until all fitting parameters converge to a certain set of values within a pre-set tolerance.
6.4. Simulation Results. In this section, both nonpolar and polar molecules are taken for the numerical simulation and model validation. Nonpolar molecules are simulated first to justify the usage of u p which represents the volume ratio of solute. That may minimize modeling uncertainties from solvent-solute electrostatic interactions. It is followed by the calculation of polar molecules to demonstrate the potential of current proposed model for the prediction of polar solvation energies.
6.4.1. Nonpolar molecules. To validate the current constrained variational model, we start with a set of 11 alkanes as a calibration set for numerical implementation of model solution and the associated parameterizaton process. First of all, two parameters N and q k need to be pre-determined for each simulation. It turns out that optimal fitting parameters are uniquely computed for a set of arbitrary N > 1 and q k , where p = 2N 2N −1 , and q k → 1 + . For instance, when N = 40 and q k = 1.00001, the calculated optimal fitting parameters are the following: γ = 0.0746 kcal/(molÅ 2 ), P h = 0.0090 kcal/(molÅ 3 ) and cs = 0.486 kcal/mol, and hs = 0.00 kcal/mol. Note that hs and cs are well depth parameters of the hydrogen and carbon, respectively. Moreover, it is shown that the current model is able to reproduce the total solvation free energies of 11 alkanes very well (see Table 1). The root mean square (RMS) error of 11 alkenes is 0.109 kcal/mol. For the nonpolar solvation free energy, the repulsive and attractive parts of solvation free energy are also calculated for detailed comparisons with others in the literature. Note that the first two terms of (2.3) are considered as the repulsive part of solvation free energy.
Next, it is interesting to see whether the model parameter N or equivalently p = 2N 2N −1 , which is introduced in the volume ratio of solute u p , plays an important role in the solvation free energy calculation and prediction. For this purpose, different N values are chosen for the set of 11 alkanes while fixing all other simulation  Table 2). Moreover, with q = 1.00001 and N = 40, a predictive study is conducted for a set of 11 alkene compounds which was also used before [19,53,55]. The assumed similar solvent environment allows one to apply the above-obtained optimized parameters of 11 alkanes here because of the fact that both nonpolar sets only possess two types of atoms (C and H). It turns out that the numerical prediction of the current model matches the experimental data well as shown in Table 3. The RMS error of 11 alkenes is 0.21 kcal/mol. Furthermore, we have theoretically proved that total solvation energies converge to the case of q k = 1 when q k → 1 + . Numerically, the convergence can be demonstrated as follow: choosing a set of molecules like the above alkene compounds and fixing all other numerical settings, one allows the value of q k to approach 1 by creating a sequence of q k (q k = 1.01, 1.001, 1.0001, 1.00001, 1.000001). Then the total solvation free energy of each molecule is computed. Table 4 illustrates the convergence of total solvation free energies for all eleven alkenes.
Remark that regarding the numerical calculation of solvation free energy for nonpolar molecules, the currently computed results are almost the same as the previous constrained solvation model [55] when N is large enough. The similarity can be explained by the fact that pu p−1 → 1 for 0 < u < 1 when p = 2N 2N −1 → 1 with N → ∞. The introduction of u p as solute volume ratio enables us to derive the system (6.2) from proposed constrained total solvation energy model (2.3). It has been a theoretical advance from our previous constrained model in which a PDE was derived only for nonpolar energy functional due to the complex two-obstacle problem [55].
In this section, the model potential and validation are demonstrated numerically for polar molecules. To the end, a challenging set of 17 compounds is chosen. The challenge arises partially due to strong solventsolute interactions caused by polyfunctional or interacting polar groups. Actually, its challenge can be seen quantitatively. For instance, using an explicit solvent model, Nicholls et al. obtained the root mean square error (RMS) as 1.71 ± 0.05 kcal/mol via [50]. With an improved multiscale model equipped with self-consistent quantum charge density by Chen et al [18], RMS was still around 1.50 kcal/mol.
For the current simulation, the structure data of the set of 17 molecules is taken from the supporting information of the paper of Nicholls et al [50] as we did before. The dielectric constants are slightly adjusted.
In the solute region m ≈ 1, while s ≤ 80 for the solvent region. For this 17 set, different well-depth parameters is need to be optimized based on the above-described simple parameterization scheme. It is shown that the computed solvation free energy is quite comparable with the experimental data. The root mean square error can be improved to 1.107 kcal/mol (See table 5) when m = 1.15 and s = 70. In addition, it is found that almost identical simulation results are obtained for large enough N . In other words, model parameter value N does not play an important role for the solvation energy prediction while it obviously benefits the theoretical derivation and the proof for current constrained variational model. The minor effect of different N values can be found in Table 6.

Conclusions
Variational implicit solvation models (VISM) with diffuse solvent-solute interface definition have been considered as a successful approach to compute the disposition of an interface separating the solute and the solvent. It has been shown numerically that variational diffuse-interface solvation models can significantly improve the accuracy and efficiency of solvation energy computation. However, there are several open questions concerning those models at a theoretic level. In particular, all existing VISMs in literature lack the uniqueness of an energy minimizing solute-solvent interface and thus prevent us from studying many important properties of the interface profile.
Therefore, by introducing a new volume ratio function u p , in this work, we have developed a novel constrained VISM based on a promising previously-proposed total variation based model (TVBVISM). Existence, uniqueness and regularity of the energy minimizing solute-solvent interface have been studied. Moreover, with the assistance of the precise depiction of the interface profile, this work provides a partial answer to the question why the solvation free energy is not minimized by a sharp solute-solvent interface. It turns out that when the mean curvature of Σ 0 is positive at some point, the energy minimizing state is never achieved by a sharp interface.
In addition, for the variational analysis of the new model and for the numerical computation of the solvation energy, a novel approach has been proposed to overcome the essential difficulty generated by the involved constraints in the model. Specifically, the variational formulas of the new energy functional can be rigorously derived via the introduction of the new volume ratio function u p together with an approximation technique by a sequence of q-energy type functionals. This is another advance from our previous work in which only the numerical study of nonpolar energy can be conducted for a constrained VISM. Model validation and numerical implementation have been demonstrated by using several common biomolecular modeling tasks. Numerical simulations show that the solvation energies calculated from our new model match the experimental data very well.
For the future work, we will provide a complete proof for the continuous dependence of the solvation free energy on the surfaces Ω m and Ω s in a suitable topology. Numerically, based on the derived elliptic system, we intend to further improve the accuracy and efficiency of the solvation energy prediction via refined parameterization schemes. Moreover, analysis of the current and potential numerical schemes like convergence will be a topic for future study.
Appendix A. BV-functions In Appendix A, we will introduce some notations and preliminaries of BV −functions. The main reference is [1,30]. Let Ω ⊂ R N be open. By the structure theorem of BV −functions, for every f ∈ BV (Ω), there exist Radon measure µ and a µ−measurable function σ : Ω → R N such that • |σ(x)| = 1 a.e. and • Ω f divz dx = − Ω z · σ dµ for all z ∈ C ∞ 0 (Ω; R N ). We write |Df | for the measure µ.
Sobolev embedding also holds for functions of bounded variations: The embedding is compact when 1 ≤ p < 1 * .
Then there exists a subsequence, not relabelled, such that f n → f in L 1 (Ω) for some f ∈ BV (Ω). An Lebesgue measurable set E ⊂ R N is said to have finite perimeter in Ω if χ E ∈ BV (Ω).
Definition A.4. Let E be of finite perimeter in Ω. We call the reduced boundary ∂E * the collection of all points x ∈ supp|Dχ E | ∩ Ω such that the limit ν E (x) := − lim r→0 + Dχ E (B(x, r)) Dχ E (B(x, r)) exists in R N and satisfies |ν E |(x) = 1 a.e.. The function ν E : ∂E * → S N −1 is called the generalized outer normal to E. ∂E \ ∂E * is called the singular set of E. In particular, we have such that where ν is the outer unit normal on ∂Ω. It is understood that the measure on ∂Ω is H N −1 . The function Trf , which is uniquely defined H N −1 a.e. on ∂Ω, is called the trace of f on ∂Ω.
Proposition A.6. Let Ω be bounded and Lipschitz. Assume that f 1 ∈ BV (Ω) and f 2 ∈ BV (R N \ Ω). Define Given f ∈ L 1 loc (Ω), we say that f has an approximate limit at x ∈ Ω if there exists z ∈ R such that lim r→0 + 1 |B(x, r)| B(x,r) |u(y) − z| dy = 0. (A. 2) The set of points where this does not hold is called the approximate discontinuity set of f , and it is denoted by S f . By Lebesgue differentiation theorem, L N (S f ) = 0. z is uniquely determined via (A) and is denoted byf (x). f is said to be approximately continuous at x if x / ∈ S f and f (x) =f (x).
We say f ∈ L 1 loc (Ω) has an approximate jump point at x ∈ Ω if there exist a = b ∈ R and µ ∈ S N −1 such that a = b and Then u is approximately continuous at x ∈ Ω iff u * (x) = u * (x).

Appendix B. Tools from convex analysis
In Appendix B, we will state some useful tools from Convex Analysis. Interested readers may refer to the books [29,58] for more details.
Let X be a Banach space with norm · . Throughout, we assume that f : X → R ∪ {±∞} is convex and lower semicontinuous (l.s.c.) function. Its effective domain is defined by is f is said to be proper if it nowhere takes value −∞ and is not identically equal to +∞ on X.
Given any subset U ⊂ X, its indicator function I U is defined by I U (x) = 0 when x ∈ U ∞ when x ∈ X \ U.
We denote by X * the topological dual of X and ·, · the duality pairing. When f is proper, the subdifferential of f at u ∈ dom(f ) is the set of all u * ∈ X * such that and is denoted by ∂f (u). Each element of ∂f (u) is called a subdifferential of f at u. When ∂f (u) = ∅, we say that f is subdifferentiable at u.
The relationship between subdifferentiability and Gâteaux-differentiability is described by the following proposition.
By the definition of the subdifferential, it is obvious that However, the converse is not always true. We list below several cases where the converse holds.