Explicit-Solute Implicit-Solvent Molecular Simulation with Binary Level-Set, Adaptive-Mobility, and GPU

We set the solvation system region to be $\Omega=(-L,L)^{3}$ for some $L>0$. The side length $L$ is chosen to be large enough so that the region $\Omega$ includes all the solute atoms ${\mathbf{r}}_{i}$ $(1\leq i\leq N)$ whose geometrical center can be shifted to the origin, if necessary; cf. Figure 1. This region $\Omega$ is also our computational box. We cover it by a uniform finite-difference grid of size $h.$ A solute-solvent interface $\Gamma=\partial\Omega_{\rm m}$ is approximated by a binary level-set function $\phi$ that is defined on all the grid cells with $\phi=-1$ and $\phi=+1$ on cells interior and exterior to $\Gamma$, respectively.

We discretize the VISM solvation free-energy functional (2.2) with all the solute atoms fixed at ${\mathbf{r}}_{i}$ $(1\leq i\leq N).$ Let us first rewrite this functional as

where

Approximation of the surface energy. The surface area of the solute-solvent interface $\Gamma$ can be expressed as [31]

where

is a constant, $B(\mathbf{0},A)$ for any $A>0$ is the ball centered at the origin $\mathbf{0}$ with radius $A$, and $y_{3}$ is the third component of the position vector $\mathbf{y}$. The kernel function $K=K(\mathbf{x})$ ($\mathbf{x}\in\mathbb{R}^{3}$) is chosen to be non-negative, compactly supported in the closure of the unit ball $B_{1}(\mathbf{0})$ of ${\mathbb{R}}^{3}$, and spherically symmetric (i.e., it is a function of $|{\mathbf{x}}|$). In our implementation, we set $K(\mathbf{x})=\sin^{2}(\pi|\mathbf{x}|)$ if $|\mathbf{x}|\leq 1$ and $0$ elsewhere. The small parameter $\delta>0$ is the rescaled kernel radius, defined as the radius of the ball of the support of $K(|{\mathbf{x}}|/\delta).$

Discretization of the VISM free energy in the solvation region. By employing the center-point numerical integration rule, one can discretize the double-integral in (3.3) with an optimal choice $\delta=\lambda\sqrt{h},$ where $\lambda>0$ is a constant. Consequently, we obtain the following expression of an approximation of the surface energy, the first term in (3.1) [31]:

where ${\mathbf{x}}_{j}\in\Omega_{\rm m}$ and ${\mathbf{x}}_{k}\in\Omega_{\rm w}$ are the centers of grid cells in $\Omega_{\rm m}$ and $\Omega_{\rm w}$, respectively. The second term in (3.1) can be approximated by the center-point integration rule:

Therefore, these approximations and (3.1) lead to the approximation

The last integral can be written analytically as iterated integrals using the spherical coordinates and evaluated by one-dimensional numerical quadrature; cf. [5].

Flipping grid cells to decrease the energy. Given a solute-solvent interface defined by a binary level-set function, we relax its VISM free energy (3.4) by flipping the grid cells, i.e., by changing the sign of the binary level-set function on the grid cells. The flipping is only done for grid cells that around the interface. This is because that any grid cells centered at ${\mathbf{x}}_{j}\in\Omega_{\rm m}$ and ${\mathbf{x}}_{k}\in\Omega_{\rm w}$ with $|{\mathbf{x}}_{j}-{\mathbf{x}}_{k}|>\lambda\sqrt{h}$ do not contribute to the first term in (3.4).

We pick up a grid cell that is immediate next to be the interface, flip its sign, and calculate the change of the approximate energy based on (3.4). Note that the last term in (3.4) does not change if we flip a grid cell. If the cell is centered at ${\mathbf{x}}_{j}\in\Omega_{\rm m}$, so the sign of the cell is $-1$, then the flip of the cell leads to the change of energy

Otherwise, if the cell is centered at ${\mathbf{x}}_{j}\in\Omega_{\rm w}$, so the sign of the cell is $+1$, then the flip of the cell leads to the change of energy

After calculating $\Delta(G_{\rm solv})_{j}$ by flipping grid cells near the interface, we put $\Delta(G_{\rm solv})_{j}$ in a Min-Heap. We flip the grid cell with the smallest $\Delta(G_{\rm solv})_{j}$ in the heap if $\Delta(G_{\rm solv})_{j}<0$. With the new interface, we add energy changes for new grid cells near the new interface to the heap, delete old grid cells in the heap which are not near the new interface, update $\Delta(G_{\rm solv})_{j}$ for old grid cells near the new interface in the heap, then we sort the Min-Heap. This flipping process stops until all $\Delta(G_{\rm solv})_{j}\geq 0$, indicating the energy reaches a minimum.

Initial surfaces. Our method for minimizing the VISM free-energy functional is of steepest descent type. It starts with an initial surface and iteratively moves it with the free energy decreased in each step of the iteration. The final free-energy minimizing surface is a local minimizer of the functional and depends on the initial surface. Different initial surfaces can lead to different (meta)stable equilibrium conformations that are of interest; see Figure 1. In order to capture multiple local minimizers, we often use two types of initial surfaces. One is a loose initial surface that can be a large sphere enclosing all the solute atoms. The other is a tight initial surface that wraps up all the solute atoms tightly with vdW radii. Such a surface is the zero level-set of the continuous function

where $d_{i}>0$ is the vdW radius of the $i$th solute atom located at $\mathbf{r}_{i}$ $(i=1,\cdots,N).$ The binary level-set function for the surface can then be obtained by setting its value at the center of a grid cell to be the sign of $\varphi$-value at that center. Figure 2 shows a tight surface constructed by both a continuous and the corresponding binary level-set function.

With a fixed solute-solvent interface $\Gamma$, we minimize the free energy $G[\Gamma,{\mathbf{R}}]$ defined in (2.1), (2.2), and (2.4) as a function of ${\mathbf{R}}=({\mathbf{r}}_{1},\dots,{\mathbf{r}}_{N})$ by solving for a steady-state solution to the system of the gradient descent equations

where $({\bf a})_{l}$ denotes the $l$th component of a vector $\bf a\in{\mathbb{R}}^{3}$ $(1\leq l\leq 3)$ and all $M_{nl}>0$ are constants (called mobility constants). The formula of the gradient $\nabla_{\mathbf{r}_{n}}G[\Gamma,\mathbf{R}]$ is given in (A.1) in Appendix. We use the forward Euler method to solve these equations iteratively with a fixed time step.

Due to the complex molecular interactions of an underlying system, the gradient of $G[\Gamma,{\mathbf{R}}]$ can vary significantly with solute atoms and with different components. If the mobility constant $M_{nl}$ is too large, then the motion of the particle $\mathbf{r}_{n}$ in its gradient descent direction may possibly overshoot and increase the free energy. If $M_{nl}$ is too small, the free energy may decrease very slowly. To improve the stability and efficiency, we therefore adaptively change the mobility constants $M_{nl}$ in each step of iteration based on the magnitude of gradient and the total free-energy value.

In our implementation, the mobility constants are chosen to be $\omega M$, where the “base” mobility $\omega$ is updated in each interaction step of the forward Euler method and $M$ is also adjustable. The adjustment of $M$ is based on the decrease or increase of the energy with a high energy threshold. It is controlled by a relative energy value $\delta G^{*}>0$ depending on iteration steps, and a shrinking parameter $\alpha\in(0,1)$ that shrinks $M$ if the energy increases too much and too often tracked by a counting number $N_{\rm cnt}$ which has a threshold value $N_{\rm cnt}^{*}$. Initially, we set $M=1$ and $N_{\rm cnt}=0.$

Given all the atomic positions after a forward Euler iteration step, with some $M>0$ and $0\leq N_{\rm cnt}<N_{\rm cnt}^{*}$, and the corresponding energy $G[\Gamma,{\mathbf{R}}]$ calculated and denoted $G_{\rm old}.$ We calculate all the gradient of $G[\Gamma,{\mathbf{R}}]$ at those positions and set

We then move the atomic positions in one time step by the forward Euler iteration for solving (3.7) with all $M_{nl}=\omega M$. We also calculate the energy for the new atomic positions and denote it by $G_{\rm new}.$ If $G_{\rm new}<G_{\rm old}$, then we accept the moved solute atomic positions, do not change the value of $M$, increase $N_{\rm cnt}$ by $1$, and continue the iteration. If $G_{\rm new}>G_{\rm old}$, then we choose the threshold $\delta G^{*}$ to be some fraction of $G_{\rm old}$ and consider two cases. If $G_{\rm new}\geq G_{\rm old}+\delta G^{*}$, then we do not update the atomic positions but shrink $M$ to $M:=\alpha M$ reset $N_{\rm cnt}=0$ and start over with the next iteration step. If $G_{\rm new}<G_{\rm old}+\delta G^{*}$, then we check the counting number $N_{\rm cnt}$ and consider two cases:

1. (1)

   If $N_{\rm cnt}{\geq}N_{\rm cnt}^{*}$, which means that the same $M$ has been used for $N_{\rm cnt}^{*}$ times, then we do not accept the new positions, shrink $M$ to $M:=\alpha M$, and reset $N_{\rm cnt}=0.$

2. (2)

   If $N_{\rm cnt}<N_{\rm cnt}^{*}$, we accept the new atomic positions, keep the same $M$, and increase $N_{\rm cnt}$ by $1$.

Note that $N_{\rm cnt}$ is the number of steps where the same $M$ is used consecutively. In our implementation, we choose $\alpha=1/\sqrt{2}$, $N_{\rm cnt}^{*}=5$, and $\delta G^{*}$ to be $5\%$ of $G_{\rm old}$.

1. Step 1.

   Initialization. Input all the model parameters from (2.2)–(2.4). In particular, position $N$ solute atoms with the center of geometry at the origin by a coordinate translation if necessary. Set the computational box $\Omega=(-L,L)^{3}$ and discretize the box with a uniform finite-difference grid. Set an initial binary level-set function $\phi^{(0)}$ to define the initial solute-solvent interface $\Gamma^{0}$, solute region $\Omega_{\rm m}^{0}$, and solvent region $\Omega_{\rm w}^{0}$.

   Set counter $N_{\rm cnt}=0$, count tolerance $N_{\rm cnt}^{*}=5$, an initial uniform mobility constant $M=1$, and the time step $dt=1$. Set the initial iteration number $k=0$ and $k_{\rm max}$$=3000$. Set the shrinking parameter $\alpha=1/\sqrt{2}$. Set the error tolerance ${\rm Tol}_{1}=1e$-$4$ for the gradient and ${\rm Tol}_{2}=1e$-$5$ for atomic positions update.

2. Step 2.

   Get the optimal solute-solvent interface $\Gamma^{k+1}$ with atomic positions ${\mathbf{R}}^{k}$ by the binary level-set method.

   1. Step 2.1.

      Calculate by (3.5) and (3.6) the change of solvation free energy $\Delta(G_{\rm solv})_{j}$ on all grid cells next to the interface $\Gamma^{k}$. Sort solvation energy change $\Delta(G_{\rm solv})_{j}$ in a Min-Heap.

   2. Step 2.2.

      Flipping Process:
      $\mathbf{while}$ (Smallest $\Delta(G_{\rm solv})_{j}<0$) $\mathbf{do}$
      $~~~~\mathbf{Flip:}$ flip the corresponding grid cell with smallest $\Delta(G_{\rm solv})_{j}$.
      $~~~~\mathbf{Update:}$ check for new grid cells next to the new interface, calculate $\Delta(G_{\rm solv})_{j}$          for new grid cells, and update $\Delta(G_{\rm solv})_{j}$ for old grid cells in the heap.
      $~~~~\mathbf{Sort:}$ sort solvation energy change $\Delta(G_{\rm solv})_{j}$ in a Min-Heap.
      $\mathbf{endwhile}$

   3. Step 2.3.

      Define $\Gamma^{k+1}$ to be the optimal solute-solvent interface from the flipping process.

3. Step 3.

   Update the solute atomic positions by solving the system of equations (3.7).

   1. Step 3.1.

      Calculate by (A.1) the gradient $\nabla_{\mathbf{r}_{n}}G[\Gamma^{k+1},\mathbf{R}^{k}]$ for all $n=1,\cdots,N.$

   2. Step 3.2.

      Test the convergence. If absolute values of $\nabla_{\mathbf{r}_{n}}G[\Gamma^{k+1},\mathbf{R}^{k}]$ for each solute atom in each coordinate $<\rm{Tol}_{1}$, then stop the algorithm.

   3. Step 3.3.

      Update positions of all moving solute atoms according to equation (3.7) and (3.8). Calculate the total free energy change $\delta G[\Gamma^{k+1},\mathbf{R}^{k+1}]=G[\Gamma^{k+1},\mathbf{R}^{k+1}]-G[\Gamma^{k+1},\mathbf{R}^{k}]$. Set $\delta G^{*}=5\%G[\Gamma^{k+1},\mathbf{R}^{k}]$.

   4. Step 3.4.

      Check the total free energy change:

      $\mathbf{if}$ ($\delta G[\Gamma^{k+1},\mathbf{R}^{k+1}]>\delta G^{*})$ or $(~N_{\rm cnt}\geq N_{\rm cnt}^{*}$ and $0<\delta G[\Gamma^{k+1},\mathbf{R}^{k+1}]\leq\delta G^{*}$)

      Put all moving solute atoms back to $\mathbf{R}^{k}$.

      $M=\alpha M$, $N_{\rm cnt}=0$, go to Step 3.3.

      $\mathbf{else}$
      $N_{\rm cnt}=N_{\rm cnt}+1;$
      $\mathbf{endif}$

4. Step 4.

   Calculate the absolute error $\Delta G[\Gamma^{k+1},\mathbf{R}^{k+1}]_{\rm abs}$ and relative error $\Delta G[\Gamma^{k+1},\mathbf{R}^{k+1}]_{\rm rel}$.

5. Step 5.

   Test the convergence. If either $\Delta G[\Gamma^{k+1},\mathbf{R}^{k+1}]_{\rm rel}<\rm{Tol}_{2}$ or $\Delta G[\Gamma^{k+1},\mathbf{R}^{k+1}]_{\rm abs}<\rm{Tol}_{2}$ stays continuously for 100 steps, or if the number of iterations reaches to $k_{\rm max}$, stop the algorithm. Otherwise, go to Step 2.

We remark that there are two error tolerances and stopping criteria: One is a small tolerance $1$e-$4$ for the gradient descent for updating the solute atomic positions. The other stop criterion is a small tolerance of relative difference or absolute difference of the total free energy. In our experiments, we set the relative difference stop criterion to be $1$e-$5$, and absolute difference stop criteria to be $1$e-$5$. To avoid the situation that the free energy functional decreases slowly because of small mobility factor $M_{nl}$, we determine that the system stops only when the relative difference stop criterion or the absolute difference stop criterion is satisfied continuously for 100 steps.

We consider two molecular systems each with fixed solute atomic positions, and apply the binary level-set method with GPU implementation to minimize the solvation free-energy functional (2.2) of solute-solvent interfaces $\Gamma$ with all atomic positions ${\mathbf{r}}_{i}$ $(1\leq i\leq N)$ fixed. Both systems have been studied extensively with continuous level-set method and CPU computations [29, 32, 34]. Here we show the qualitative accuracy, and efficiency, of our new algorithm and implementation.

We first investigate the VISM surfaces of the two plates at different distances with the like-charge (+0.2e, +0.2e) with different initial configurations. Figure 3 shows a few snapshots of stable 3D equilibrium solute-solvent surfaces of the two parallel charged plates system obtained by the binary level set VISM calculations with loose or tight initial interface at $d=9\,\mathring{A}$, $d=13\,\mathring{A}$, and $d=16\,\mathring{A}$. In the top row of Figure 3, with the loose initial interface, a stable capillary bubble remains between the two charged parallel plates at $d=9\,\mathring{A}$ and $d=13\,\mathring{A}$, and the bubble becomes tighter along the enlarging distance. At $d=16\,\mathring{A}$, the bubble disappears. Comparatively, with the tight initial interface, the equilibrium state is wet at $d=9\,\mathring{A}$, $d=13\,\mathring{A}$, and $d=16\,\mathring{A}$.

We now examine the potential of mean force (PMF) of the two-plate system with respect to the plate-plate separation distance $d$. This is the VISM free-energy value as a function of $d$, with an additive constant such that the free energy is $0$ at the infinite plate-plate separation. For a given $d$, we may have two VISM free-energy minimizing solute-solvent interfaces corresponding to a tight and a loose initial surface, respectively. We denote by $G_{\rm VISM}^{\rm pmf}(d)$ the corresponding minimum free energy of one of the two branches, and denote by $G_{\rm geom}^{\rm pmf}(d),$ $G_{\rm vdW}^{\rm pmf}(d),$ and $G_{\rm elec}^{\rm pmf}(d)$ the components of the PMF, corresponding to the first, second, and third terms in (2.2), respectively. Precise definition is given in Appendix.

Figure 4 displays the bimodal behavior and hysteresis of the two different PMF branches stemming from the equilibria of wet and dry states, i.e., the VISM free-energy minimizing surfaces corresponding to initial tight and loose surfaces. Atomic charges considered here are (+0.2e,-0,2e) and (+0,2e, +0,2e), respectively. We can see that like-charged and oppositely charged plates give different free-energy branches and hysteresis. For the like-charged cases in Figure 4, a strong hysteresis is presented for $8\lesssim d\lesssim 15\mathring{A}$. For the oppositely charged plates, strong hysteresis is presented for $6\lesssim d\lesssim 8\mathring{A}$.

In Table 2, we show a comparison of the calculation speed and different parts of the free energy using the binary level-set method between GPU single precision code and CPU double precision code of the two charged parallel plates system. Three different grid sizes are shown. We can observe that the results for free-energy estimates from CPU double precision code and GPU single precision code are nearly the same. The improvement of speed using GPU can be obtained by comparing the time. The cost time of CPU code is around 5 times of the GPU code with three different grid sizes.

The functional unit of this protein is a homo-octamer, and each subunit consists of two domains. We set up a series of configurations where the two domains are increasingly separated from $d=0$ to $d=20\,\mathring{A}$ apart, perpendicular to their interface. The domain separation $d$ is chosen here to be the reaction coordinate. Note that the zero domain separation corresponds to the native configuration in the crystal structure.

Three pairs of stable equilibrium solute-solvent interfaces of BphC at $d=8\,\mathring{A}$, $d=12\,\mathring{A}$, and $d=16\,\mathring{A}$ with tight or loose initial interfaces are presented in Figure 5. The top row is with the loose initial interfaces, and the bottom row is with the tight initial interfaces. We observe that the equilibria of the loose initial interface wrap the two domains of BphC at $d=8\,\mathring{A}$ and $d=12\,\mathring{A}$, and the equilibria surface becomes tighter along the increasing distance. At $d=16\,\mathring{A}$, the interfaces of two domains separate, changing to the wet state. In contract, with the tight initial interface, all equilibria states are wet.

In Figure 6, different parts of the PMF of BphC with respect to the separation of two domains, from $d=0$ to $d=20\mathring{A}$ obtained by our binary level set calculations using tight and loose initial surfaces are displayed. These PMFs exclude the solute-solute vdW interactions. We observe the bimodal hydration behavior: the branches of different parts of the PMF of BphC between $4\,\mathring{A}$ and $14\,\mathring{A}$, indicating that initial interfaces can strongly affect the PMF of BphC.

In Table 3, we show a comparison of the calculation speed and different parts of the free energy with the binary level-set method between GPU single precision code and CPU double precision code of the BphC. Three different grid sizes are shown. We observe that the results from CPU double precision code and GPU single precision code are nicely consistent. With the grid size of $50^{3}$, $100^{3}$, and $200^{3}$, the time cost of CPU code is around 15 times, 78 times, and 216 times correspondingly to the time cost of the GPU code.

In this section, we conduct numerical experiments on two molecuales that are treated as nonpolar systems (i.e., no charges) to demonstrate the efficiency of our free-energy minimization algorithm.

We design two sets of experiments on the optimization process and equilibria with different initial configurations. In the first set of experiments, Experiment 1.1.a and Experiment 1.1.b, we set the equilibrium bond length $r_{0}=3\,\mathring{A}$. In the second set of experiments, Experiment 1.2.a and Experiment 1.2.b, we set the equilibrium bond length $r_{0}=8\,\mathring{A}$. In each set of experiments, we test two types of initial configurations. In Experiment 1.1.a and Experiment 1.2.a, we place initially the two solute atoms far away from each other so that their distance is much larger than the equilibrium bond length. We place the two solute atoms at positions $(7,0,0)$ and $(-7,0,0)$, respectively. In Experiment 1.1.b and Experiment 1.2.b, we place initially the two solute atoms very close to each other so that their distance is smaller than the equilibrium distance. Specifically, we place the two solute atoms initially at positions $(0.5,0,0)$ and $(-0.5,0,0)$, respectively.

In Figure 7, the minimization processes of Experiment 1.1.a and Experiment 1.1.b are displayed. The red dots represent two atoms, and the blue segment represents the bond. In the top row of Figure 7, we can see that initially surface consists of two disconnected spheres, then two atoms get closer, and spheres merge, until the system reaches an equilibrium state. In the bottom row of Figure 7, initially, the two atoms are very close to each other, then the atoms are pushed apart due to the force from strong bonding energy, the interface is moved accordingly, and then the system reaches to an equilibrium state.

We observe from Table 4 that that atoms have the exact same positions and free energy in the equilibrium from the two experiments 1.1.a and 1.1.b, indicating that the molecular system in the two simulations reached the same equilibrium.

Figure 8 shows the minimization processes of Experiment 1.2.a and Experiment 1.2.b. In the top row of Figure 8, we can see that initially surface consists of two disconnected spheres, then two atoms get closer, until the system reaches an equilibrium state. Comparing with the equilibrium of Experiment 1.1.a, the spheres are not merged due to a larger bond length $8\mathring{A}$. In the bottom row of Figure 8, initially, the two atoms are very close to each other, then the atoms are pushed apart due to the force from strong bonding energy, the interface moves and then splits apart, then the system reaches to an equilibrium state with which the interface consists of two separate spheres. Table 5 shows that the two experiments 1.2.a and 1.2.b reach the same equilibrium.

In the ethane molecule, each atom is connected to other atoms through bonding, bending, and torsion structure. The effect of vdW interaction energy among unbonded pairs of solute atoms is relatively less important when compared with molecular mechanical interactions, so in our simulation experiment of the ethane molecule, we neglect the vdW interaction energy among unbonded pairs of solute atoms.

We design three different initial configurations of the ethane molecule from its equilibrium:

1. $\bullet$

   In Experiment 2.1.a, we stretch all hydrogen-carbon bonds to be $2\mathring{A}$.

2. $\bullet$

   In Experiment 2.1.b, we stretch or shrink all hydrogen-carbon bonds such that hydrogen-carbon bonds have the length of $1\,\mathring{A}$, $1.25\,\mathring{A}$, $1.5\,\mathring{A}$, $1.75\,\mathring{A}$, $2\,\mathring{A}$, and $2.25\,\mathring{A}$.

3. $\bullet$

   In Experiment 2.2, we introduce a small fluctuation, then rotate one set of three hydrogen-carbon bonds 50 degrees with respect to the carbon-carbon bond.

In Table 6, the bond lengths of hydrogen-carbon bonds and the carbon-carbon bond in their equilibrium of Experiment 2.1.a, Experiment 2.1.b, and Experiment 2.2 are compared. It is clear that in the equilibrium, all hydrogen-carbon bonds in three experiments have the same length, which is consistent with the reference length of hydrogen-carbon bond $1.093\,\mathring{A}$. The carbon-carbon bond in each of the three experiments is the same as the reference length $1.508\,\mathring{A}$. This verifies the accuracy of our free-energy minimization algorithm.

Figure 9 displays the snapshots of minimization process of Experiment 2.2. The red dots represent carbon atoms, the blue dots represent the hydrogen atoms, the light blue segments and green segments represent the hydrogen-carbon bonds, and the black segment represents the carbon-carbon bond. It is captured that during the relaxation, the set of three hydrogen-carbon bonds rotated back to their equilibrium, and all hydrogen-carbon bonds have the same length. Free energy of steady state in Experiment 2.1.a is 3.692 $k_{B}T$, the free energy of steady state in Experiment 2.1.b is 3.685 $k_{B}T$, and the free energy of steady state in Experiment 2.2 is 3.687 $k_{B}T$. Thus, the three experiments get to the same equilibria. We remark that the free energy here does not include the vdW interaction energy among unbonded pairs of solute atoms.

In Figure 10, we plot the free energy vs. iteration steps in numerical computations for Experiment 2.1.a and Experiment 2.2. In Experiment 2.1.a, initially, the free energy is very large, that is because the stretch or shrink of the initial bond length causes a large value of the bonding energy. We can see that the free energy decays very fast in the first few steps, which is caused by the dominant force from the bonding energy. It takes around 1000 steps to adjust positions of solute atoms in ethane molecule to reach the equilibrium. In contract, the initial free energy of Experiment 2.2 in Figure 10 is a relatively small value, as we only introduced a small fluctuation to the initial atomic positions and the rotation of a set of three hydrogen-carbon bonds did not cause large free energy change. We observe that the rate of free-energy change at around the 3rd step and 200th computational step becomes slower and slower, which indicates that our free-energy minimization algorithm was adjusting the suitable mobility factor $M$ during the minimization process. Although the initial free energy is small, it takes more than 1400 steps to rotate the hydrogen-carbon bonds back to the right position and reach the equilibrium.

In this section, we choose a biologically important and realistic system, the p53/MDM2 protein complex, and investigate the binding behavior using our free-energy minimization model and algorithm. To make the calculations of molecular movement easier, the receptor protein MDM2 here is fixed.

We use the CHARMM36 force field [15, 2, 21, 20] for our VISM simulations for the binding of p53 and MDM2.

Table 7 shows the solvation free energy and its components obtained by our VISM simulations for the bound complex p53/MDM2, and the computational times of the simulations with the GPU single precision and the CPU double precision, respectively.

During those simulations, we relax the relative difference stop criterion and absolute difference stop criteria to be $1$e-$3$, just for the efficiency comparison of the GPU code and the CPU code. We set the initial configuration of p53/MDM2 to be a tight initial interface for atomic position with small fluctuations of p53/MDM2 in the bound complex, where the bound complex is taken from the Protein Data Bank (PDB code: 1ycr.pdb). It can be observed in Table 7 that the difference between GPU code with single precision data type and CPU with double precision data type is less than 1%, but the cost time of CPU is more than 78 times, and 307 times slower than the cost time of GPU with number of grid points $50^{3}$, and $100^{3}$, correspondingly.

We further investigate the binding behavior of p53 and MDM2 with our free-energy minimization method and algorithm. Note that, here we use a grid size $h=0.5\,\mathring{A}$, the relative difference stop criterion and absolute difference stop criterion are $1$e-$5$. We construct the initial configuration with a tight initial interface by pulling p53 away from the MDM2 pocket in the bound complex along the line passing through the geometrical centers.

In Figure 11, a few snapshots of numerical results of p53/MDM2 are displayed, showing the minimization process. The position of the protein MDM2 is fixed, positions of p53 atoms are adjusted by the free-energy minimization process. We color each piece of surface according to whether its closet solute atom comes from MDM2 (red) or p53 (blue) to show the relative positions of MDM2 and p53. In the initial configuration, the protein p53 and the receptor protein MDM2 are separated as we can observe a hole between them; that is a small region filled with water. In the process, the relative positions of p53 and MDM2 are adjusted, p53 and MDM2 become closer and closer. In the equilibrium, the hole disappears, and the two proteins are combined together.