Get 20M+ Full-Text Papers For Less Than $1.50/day. Start a 14-Day Trial for You or Your Team.

Learn More →

An Adaptive Local Iterative Method for Fast Calculation of Electrostatic Interactions between Charged Polymers in Dielectric Inhomogeneous System

An Adaptive Local Iterative Method for Fast Calculation of Electrostatic Interactions between... IntroductionElectrostatic interactions between the charged polymer and biomolecule‐solvent environment play an important role in many biological processes, since many diffusing molecules of interest in biological systems carry charges.[1–7] It is well known that when the dielectric value in an electrolyte solution is represented by the segmental dielectric function, this belongs to the research category of dielectric inhomogeneity, which helps us understand the structure, function and dynamics of biomolecules, etc. In the experiment investigating phase separation of block copolymers, when the dielectric inhomogeneity is considered, the ions added to the solution tend to be enriched in the high dielectric constant components, thereby promoting phase separation, which is beyond the original cognition.[8–11] Furthermore, a particularly interesting aspect of this work is the consideration of the dependence of the dielectric coefficient on the electric field strength in electrolyte solutions, since the solvent dipole is highly oriented under a large electric field and thus affects the dielectric distribution.[12,13] The heterogeneity of the dielectric coefficient dramatically affects the structure of the bilayer region and charge transfer. Therefore, understanding the dielectric properties of solvents is essential for the accurate characterization of biomolecules in molecular level research, which can reasonably explain the electrophoresis of charged colloids in high‐valent salt solutions presenting charge reversal and the ability of ions to promote the phase separation of block copolymers.In such a dielectric inhomogeneous system, a key challenge is how to solve electrostatic interactions quickly. For instance, the biomolecule‐solvent two medium model is a classical dielectric inhomogeneity system.[14] Obviously, since electrolyte solution systems contain biomolecules that usually dissociate a huge number of charged particles, calculating the electrostatic interactions of this system is an extremely time‐consuming task. A priori, for an N‐particle system, using traditional methods[15–18] (e.g., Ewald summation, particle–particle particle‐mesh, Coulomb force formula, etc.) to calculate these interactions would seem to require O(N2)$O({N}^{2})$ operations. We also need to recognize that in systems governed by the Poisson equation under dielectric inhomogeneity, the electrostatic forces on atom are not simply the product of electrostatic field E at the atom and the atomic charge q. However, these computational efforts may consume a lot of computer memory, especially when sufficient statistical sampling requires hundreds of thousands of simulations or a single simulation for extremely long time. Therefore, there is an urgent need for an efficient method to save computing time. In the past decades, these challenges have been addressed by classical numerical methods (e.g., finite difference method (FDM), finite element method (FEM)) as well as through popular linear or nonlinear iterative methods (e.g., successive over‐relaxation, conjugate gradient method, and multigrid method).[19–22] Among them, FDM is widely favored by researchers due to its theoretical simplicity and convenience in dealing with regular boundaries, while optimization efficiency and accuracy are two core problems for FDM research. A large number of excellent works have appeared one after another in this field.[23–28] In this study, we develop a hybrid approach to combine the advantages of Brownian dynamics and finite difference methods in the same simulation.[29,30] In this way, the calculation of electrostatic interactions is converted into the problem of solving the 3D Poisson equation under a dielectric inhomogeneous system in numerical techniques. The partial differential equations of this system are well defined at boundary and interface conditions. It can be solved approximately by the finite difference method as a sparse linear matrix problem. Consequently, it is feasible to adopt it for the electrostatic calculations of large‐scale biomolecular systems. In addition, another challenge for us is to deal with the molecular interface problem accurately and efficiently to achieve a smooth extension of the potential values over different dielectric regions. The main difficulties include three points: 1) the complex geometric configuration of the molecular surface; 2) the dielectric coefficient and electric field are discontinuous across the molecular interface; 3) the existence of singular sources. Recently, Egan[31] introduced a second‐order exact geometric discretization of point charge sources to solve the Poisson equation. In order to deal with the boundary conditions of the Poisson problem on irregular domain, Towers[32] developed a central difference discretization of the Dirac‐delta function representing the singular source at the interface. Moreover, some works have been dedicated to the study of the regularization in the setting of smoothly varying dielectric functions on diffuse interfaces and complex structures, of which four regularization methods based on Green's function are the most prevalent.The first approach decomposes the numerical solution into two fractions over the entire computational domain: the singular component (Green's function) and the residual component (the regular component). Despite the lack of accuracy, this scheme is one of the pioneering works in the field and many other subsequence works are motivated by it.[33,34] Chern et al.[35] proposed the second method: a regular format separation for the solution, which treats the regular and singular parts of the solution separately. This decomposition method has been applied to the subsequent finite element methods of Lu and Zhou et al.[36,37] The third method was proposed and explored in Xie's work.[38] This scheme divides the electrolyte solution into three parts over the entire region, which is applied to solvents that contain any number of ionic species. However, Xie's and Chern's schemes are both of finite element aspect. The fourth method, matched interface and boundary (MIB) technique based on finite difference, was proposed by Zhou and Wei et al.[39–45] The MIB method to smoothly extend the both sides of the interface by means of fictitious domains, then obtains a central finite difference scheme with high‐order accuracy repeatedly utilizing the interface conditions.However, in order to numerically solve the dynamic processes of electrostatic forces into the typical Brownian dynamics simulations of molecular conjugation or dissociation, which may involve millions or even tens of millions of time steps, the existing numerical methods are difficult to load such a large calculation amount. This means that it is necessary to find a new approach to avoid expensive computational costs and to ensure the accuracy of millions or even more calculations. In our previous work, we have proposed an optimization strategy that substantially improves the iteration speed and have used the method to calculate electrostatic interactions in the dielectric homogeneous system with excellent accuracy.[46]In this paper, the motion of charged polymers in a two‐dielectric model was simulated using finite difference scheme and simple yet effective interface treatment. We propose an iteration optimization strategy for large charged systems and further optimize the algorithm according to the shape of the polymer, which greatly improves the computational efficiency and saves a substantial CPU time. The rest of this paper is organized as follows. Section 2 presents the simulation model and the calculation scheme of interface problem, and proposes the corresponding algorithmic optimization strategy. In Section 3, we first verify the viability of the SOR global iteration method based on the FDM for the calculation of electrostatic forces. In addition, we are devoted to discussing the effects of the parameters in the “local iteration” on the computational speed and accuracy. The computational efficiency and adaptability of the banded local iteration method are further validated extensively. Finally, we draw some conclusions, and provide brief overview and some prospections in Section 4.Simulation Model and MethodsThe Biomolecule‐Solvent and Polymer ModelWe construct a coarse‐grained computational model for the movement of the charged polymers in a solvated biomolecular system, which can be described by the Poisson–Boltzmann equation (PBE) or viewed as a Poisson equation carrying a series of singular charge sources. As shown in Figure 1, the computational domain Ω is separated by the molecular surface Γ into the solute (biomolecule) region Ωm$\Omega_{\text{m}}$ and the solvent region Ωs$\Omega_{\text{s}}$, that is, Γ = Ωm∩Ωs$\Omega_{\text{m}}\cap \Omega_{\text{s}}$, Ω = Ωm⋃Ωs$\Omega_{\text{m}}\bigcup \Omega_{\text{s}}$. The boundary of the entire domain is ∂Ω$\partial \Omega$. In continuum models, biomolecules are usually approximated as homogeneous dielectrics, and the relative permittivity is taken to be between 1 and 20 (here set it as 1). In our simulations, in order to capture some important physical properties of biological systems more accurately and to further improve the quality of the classical Poisson model in electrostatic models, the dielectric value of the solvent region is replaced by a dielectric function dependent on the electric field strength.[47] The dielectric coefficient of the whole computational domain is a piecewise constant function given as1ε(r)=εm,r∈Ωmn2+απN0(n2+2)μν4πε0ELβ(n2+2)μνEkBT,r∈Ωs$$\begin{eqnarray} \varepsilon (r)=\left\{ \def\eqcellsep{&}\begin{array}{ll} {\varepsilon}_{\text{m}},& r\in \Omega_{\text{m}}\\ [6pt] {n}^{2}+\displaystyle\frac{\alpha \pi {N}_{0}({n}^{2}+2){\mu}_{\nu}}{4\pi {\varepsilon}_{0}E}L\left[\frac{\beta \left[({n}^{2}+2)\right]{\mu}_{\nu}E}{{k}_{\text{B}}T}\right],& r\in \Omega_{\text{s}} \end{array} \right.\nonumber\\ \end{eqnarray}$$where εm is the biomolecule region dielectric constant (the interface dielectric constants are assigned as the harmonic average of the solvent and solute dielectric constants), n is the optical refractive index, N0 is the number of molecules per unit volume, μν${\mu}_{\nu}$ is the dipole moment of the water molecule, ε0 is the vacuum dielectric constant, and E is the field strength. L(x)$L(x)$ denotes the usual Langevin function, α and β are numerical factors, T and kB represent the absolute temperature and Boltzmann constant, respectively. In our numerical implementation, the molecular parameters at T=300K$T=300\ \text{K}$ have values as following, n=1.33$n=1.33$, N0=3.334×1028m−3${N}_{0}=3.334\ensuremath{\times{}}{10}^{28}\ {\text{m}}^{-3}$, μν=2.1D${\mu}_{\nu}=2.1\ \text{D}$, α=28/373$\alpha =28/3\sqrt{73}$, β=73/6$\beta =\sqrt{73}/6$.1FigureSchematic illustration of the computational domain Ω: biomolecule‐solvent model. The molecular surface Γ divides the biomolecule region Ωm$\Omega_{\text{m}}$ and the solvent region Ωs$\Omega_{\text{s}}$.In the biomolecular region, the charge distribution is usually fixed because the charged particles are bound by chemical bonds and are not easily affected by the surrounding electrostatic field. For the solvent region, however, the charged particles (Na+${\text{Na}}^{+}$, Cl−${\text{Cl}}^{-}$, K+${\text{K}}^{+}$, Mg2+${\text{Mg}}^{2+}$, etc.) in the solution will change their spatial positions frequently due to the influence of the surrounding electrostatic field and the collision effect of the particles. Therefore, we use a probabilistic approach to quantitatively describe the distribution of charged particles in solvent. By Gauss's law, the system of equations for the electrostatic potential u(r)$u(r)$ of the model is expressed as2εm∇2u(r)=−ρm(r),r∈Ωm$$\begin{equation} {\varepsilon}_{\text{m}}{\nabla}^{2}u(r)=-{\rho}_{\text{m}}(r),\quad r\in \Omega_{\text{m}} \end{equation}$$3∇·(εs(r)∇u(r))=−ρs(r),r∈Ωs$$\begin{equation} \nabla \cdot ({\varepsilon}_{\text{s}}(r)\nabla u(r))=-{\rho}_{\text{s}}(r),\quad r\in \Omega_{\text{s}} \end{equation}$$4[u]Γ=0,r∈Γ$$\begin{equation} {[u]}_{\Gamma}=0,\quad r\in \Gamma \end{equation}$$5[εuξ]Γ=0,r∈Γ$$\begin{equation} {[\varepsilon {u}_{\xi}]}_{\Gamma}=0,\quad r\in \Gamma \end{equation}$$6u(r)=ub(r),r∈∂Ω$$\begin{equation} u(r)={u}_{b}(r),\quad r\in \partial \Omega \end{equation}$$where ρm(r)${\rho}_{\text{m}}(r)$ can be expressed by the sum of a series of Dirac‐Delta functions, which are redistributed to their neighboring grid points with appropriate interpolation formats. Thus the second‐order or higher‐order numerical accuracy at Cartesian grid points is realized. ρs(r)${\rho}_{\text{s}}(r)$ denotes the charge density distribution function in solvent region. [·]Γ${[\cdot ]}_{\Gamma}$ denotes the jump on Γ of the value of the function from Ωm$\Omega_{\text{m}}$ to Ωs$\Omega_{\text{s}}$. Equations (2) and (3) both represent dimensionless model equations. For example, the original equation of Equation (2) is εm∇2Φ=−ρm(r),r∈Ω${\varepsilon}_{\text{m}}{\nabla}^{2}\Phi=-{\rho}_{\text{m}}(r),r\in \Omega$. By the formula 1V=1J/C$1V=1J/C$, the unit of eckBT$\frac{{e}_{c}}{{k}_{\text{B}}T}$ is found to be 1/V$1/V$. Thus, by setting u=eckBTΦ$u=\frac{{e}_{\mathrm{c}}}{{k}_{\text{B}}T}\Phi$. Clearly, both ρ(r)$\rho (r)$ and u have unit m−2${\text{m}}^{-}2$ due to C2/(FJ)=1${C}^{2}/(\textit{FJ})=1$. Hence, the both sides of each equation of (2) have the same unit m−2${\text{m}}^{-2}$, which can be removed from the both sides of each equation to make (2) become dimensionless. The solution u of the Poisson model is called the dimensionless electrostatic potential. When u is found, the electrostatic potential Φ can be simply said to have a value of u in units kBT/ec${k}_{\text{B}}T/{e}_{\mathrm{c}}$. Equations (4) and (5) are the interface jump conditions across the molecular surface with ξ as the outer normal direction, where uξ=∂u∂ξ${u}_{\xi}=\frac{\partial u}{\partial \xi}$ is the directional derivative in ξ. Mathematically, the uniqueness of the solution is ensured by jump conditions. Physically, these conditions are necessary for the continuity of the electrostatic potential and its fluxes across the interface. Here, for grid‐based approach, the computational domain Ω is still assumed as a cube. As shown in Figure 2, we fix the electric potential on the left side of the simulation box as U and further get a Dirichlet boundary condition imposed on ∂Ω$\partial \Omega$.2Figure3D schematic of a charged polymer diffusing under dielectric inhomogeneity, where the dimensions of the simulation box are (L,L,L)$(L,L,L)$ and the left surface of the electric potential is fixed to U and the others to 0.We put a polymer into the domain Ω, which is modeled as a conventional bead‐spring chain consisting of N beads connected via the finitely extensible nonlinear elastic (FENE) potential[48]7UFENE(r)=−0.5kbb2ln(1−r2/b2)$$\begin{equation} {U}_{\text{FENE}}(r)=-0.5{k}^{b}{b}^{2}\text{ln}(1-{r}^{2}/{b}^{2}) \end{equation}$$the excluded volume interactions between the beads are considered through the Lennard–Jones (LJ) potential[49–51]8ULJ(r)=4ε[(σ/r)12−(σ/r)6]+ε,r≤21/6σ$$\begin{equation} {U}_{\textit{LJ}}(r)=4\varepsilon [{(\sigma /r)}^{12}-{(\sigma /r)}^{6}]+\varepsilon ,\quad r\le {2}^{1/6}\sigma \end{equation}$$where r and b are the bead–bead distance and the maximum bead–bead separation distance between consecutive beads along the polymer chain, kb${k}^{b}$ denotes the spring constant. ε represents the strength of the LJ potential, σ is used as the unit of length in this work and is the collision diameter of the LJ potential. Additionally, the charge of beads in the polymer is assumed to be q. The Poisson equation is solved with the interface and boundary potential conditions to compute the electrostatic interactions between charged beads. Hence, the Brownian dynamics can be applied to simulate the diffusion of charged polymers, including a noise term and a friction term. After specifying the interaction potentials, the BD equation of motion for particle i is9r⃗̇i(t)=−1μ▽Ui+1μf⃗i(t)$$\begin{equation} {\dot{\vec{r}}}_{i}(t)=-\frac{1}{\mu}\bigtriangledown {U}_{i}+\frac{1}{\mu}{\vec{f}}_{i}(t) \end{equation}$$where Ui${U}_{i}$ denotes the sum of all the interaction potentials discussed above, μ is the friction coefficient. f⃗i(t)${\vec{f}}_{i}(t)$ represents the random force associated with μ through the fluctuation dissipation theorem ⟨f⃗i(t)f⃗j(t′)⟩=2μkBTδijδ(t−t′)I$ \langle {\vec{f}}_{i}(t){\vec{f}}_{j}({t}^{{}^{\prime}}) \rangle =2\mu {k}_{\text{B}}T{\delta}_{\textit{ij}}\delta (t-{t}^{{}^{\prime}})\mathbf{I}$, where I is the unit tensor. The equation of motion of the charged polymer is integrated based on the time step Δt$\Delta t$. In our model, we present our results in reduced units, which result in a unit of the time as t=t*(ε/mσ2)1/2$t={t}^{\ast}{(\varepsilon /m{\sigma}^{2})}^{1/2}$, the charge as q=q*/(4πε0σε)1/2$q={q}^{\ast}/{(4\pi {\varepsilon}_{0}\sigma \varepsilon )}^{1/2}$, and the electric field as E=E*(4πε0σε)1/2σ/ε$E={E}^{\ast}{(4\pi {\varepsilon}_{0}\sigma \varepsilon )}^{1/2}\sigma /\varepsilon $. Where σ, ε and m denote the units of length, energy and mass, respectively. We assume that σ = ε = m = 1 and the temperature is kBT${k}_{\text{B}}T$ = 1.0, and the friction coefficient is set as μ=0.5Δt−1$\mu =0.5 \Delta {t}^{-1}$. Additionally, the dimensionless parameters of the polymer are fixed as kb=30${k}^{b}=30$ and b=1.5$b=1.5$, the charge of beads in the polymer is q=10$q=10$, and the time step is chosen as Δt=0.0002$\Delta t=0.0002$. We fix the size of the simulation box as L=80σ$L=80\sigma $, the electric potential of left surface as U=4$U=4$.Improved Interface Computation SchemeThis work mainly focuses on the fast solution of electrostatic interactions between charged particles in polymers with an applied electric field under the dielectric inhomogeneous system. While considering a uniform Cartesian mesh division of the domain Ω, the standard finite difference format loses its designed convergence around the interface. Hence the interface jump condition must be used to recover the accuracy. In order to achieve a smooth extension of the electric potential values, it is necessary to derive the computational expressions of the molecular interface, which is crucial for the situations where the interface intersects or does not intersect with the grid nodes. For this purpose, we divide all grid points into two categories: regular and irregular grid points. An irregular grid point is defined as a node in the standard finite difference format containing grid points that span the interface, that is, at least one of its six neighboring points is from different region. Otherwise, the points far from the interface Γ are regular points, performing the standard 7‐point central difference to discretize Equations (2) and (3). As shown in Figure 3, the blue and orange solid dots represent irregular points from the solute and solvent regions, respectively. The green and pink dots represent interface points that intersect and do not intersect with grid nodes, respectively. In our scheme, the difference format of irregular points is derived first, and then the computational format of interface points obtained by discretizing the interface conditions is substituted to modify the finite difference approximation of irregular points.3FigureThe interface value calculation scheme illustrated in a 2D setting with outside irregular points in blue and inside irregular points in orange. Green and pink dots indicate interface points that intersect and do not intersect with grid nodes, respectively. The light blue area denotes the solute area, the purple area denotes the irregular point distribution area, and the rest is the solvent area.For the electrostatic equations in 3D space with irregular dielectric interfaces, there are two interface conditions at any direction along the interface. In this way, the continuity of the potential flux Equation (5) is directed toward the normal direction ξ. At a given intersection point, it is convenient to introduce a coordinate transformation from (ξ, η, ζ) to (x, y, z), where η and ζ are the tangential and subnormal directions, respectively. Thus, the flux continuity condition Equations (4) and (5) are projected locally into the Cartesian grid10∂u∂η=(−uxssinθ+uyscosθ)−(−uxmsinθ+uymcosθ),$$\begin{equation} \left[\frac{\partial u}{\partial \eta}\right]=(-{u}_{x}^{\mathrm{s}}\sin\theta +{u}_{y}^{\mathrm{s}}\cos\theta )-(-{u}_{x}^{\mathrm{m}}\sin\theta +{u}_{y}^{\mathrm{m}}\cos\theta ), \end{equation}$$11∂u∂ζ=(−uxssinψcosθ−uyssinψsinθ+uzscosψ)−(−uxmsinψcosθ−uymsinψsinθ+uzmcosψ),$$\begin{eqnarray} \left[\frac{\partial u}{\partial \zeta}\right] &=& (-{u}_{x}^{\mathrm{s}}\sin\psi \cos\theta -{u}_{y}^{\mathrm{s}}\sin\psi \sin\theta +{u}_{z}^{\mathrm{s}}\cos\psi )\nonumber\\ &&-\ (-{u}_{x}^{\mathrm{m}}\sin\psi \cos\theta -{u}_{y}^{\mathrm{m}}\sin\psi \sin\theta +{u}_{z}^{\mathrm{m}}\cos\psi ), \end{eqnarray}$$12ε∂u∂ξ=εs(uxscosψcosθ+uyscosψsinθ+uzssinψ)−εm(uxmcosψcosθ+uymcosψsinθ+uzmsinψ),$$\begin{eqnarray} \left[\varepsilon \frac{\partial u}{\partial \xi}\right] &=& {\varepsilon}^{\mathrm{s}}({u}_{x}^{\mathrm{s}}\cos\psi \cos\theta +{u}_{y}^{\mathrm{s}}\cos\psi \sin\theta +{u}_{z}^{\mathrm{s}}\sin\psi )\nonumber\\ &&-\ {\varepsilon}^{\mathrm{m}}({u}_{x}^{\mathrm{m}}\cos\psi \cos\theta +{u}_{y}^{\mathrm{m}}\cos\psi \sin\theta +{u}_{z}^{\mathrm{m}}\sin\psi ), \end{eqnarray}$$where ψ denotes the angle between ξ and its projection on x–y plane while θ is the angle between this projection and x‐direction. uxs${u}_{x}^{\mathrm{s}}$ and uxm${u}_{x}^{\mathrm{m}}$ denote the first order derivatives in the x direction from different regions, respectively. In a uniform grid, the interface either intersects or does not intersect with the grid points. We first discuss the case where the interface intersects the grid points. At this time, all the irregular points use the successive over‐relaxation (SOR) iteration scheme with standard 7‐point difference13ui,j,k(n+1)=ui,j,k(n)+λ16ui+1,j,k(n)+ui−1,k(n+1)+ui,j+1,k(n)+ui,j−1,k(n+1)+ui,j,k+1(n)+ui,j,k−1(n+1)+h2ρi,j,kεi,j,k−ui,j,k(n)$$\begin{equation} {u}_{i,j,k}^{(n+1)}={u}_{i,j,k}^{(n)}+\lambda \left[\frac{1}{6}\left({u}_{i+1,j,k}^{(n)}+{u}_{i-1,k}^{(n+1)}+{u}_{i,j+1,k}^{(n)}+{u}_{i,j-1,k}^{(n+1)}+{u}_{i,j,k+1}^{(n)}+{u}_{i,j,k-1}^{(n+1)}+{h}^{2}\frac{{\rho}_{i,j,k}}{{\varepsilon}_{i,j,k}}\right)-{u}_{i,j,k}^{(n)}\right] \end{equation}$$where i, j and k denote the node position in the discretized grid, h is the grid spacing in the x, y, and z directions, the superscript n and n+1 represent the current iteration step and next iteration step, respectively. λ is the relaxation factor. It is possible to estimate numerically the optimal value of the parameter λ in the SOR method in some cases.[28] For instance, the optimal relaxation factor can be obtained according to the dimensions of the computational domain for both solvent and solute, that is, λoptimal=2/[1+sin(πh/Lmin)]${\lambda}_{\text{optimal}}=2/[1+\text{sin}(\pi h/{L}_{\text{min}})]$, where Lmin${L}_{\text{min}}$ is the minimum grid number of the simulation box in three directions. For the irregular points ui,j,k${u}_{i,j,k}$ in Equation (13), we assume that one of its six adjacent nodes, ui,j+1,k${u}_{i,j+1,k}$, is an interface point, then we need to combine Equations (5) and (12) to derive the interface point formula rather than use its value directly,14εscosψcosθui,j+1,k−ui−1,j+1,kh+cosψsinθui,j+1,k−ui,j,kh+sinψui,j+1,k−ui,j+1,k−1h=εmcosψcosθui+1,j+1,k−ui,j+1,kh+cosψsinθui,j+2,k−ui,j+1,kh+sinψui,j+1,k+1−ui,j+1,kh$$\begin{align} & {\varepsilon}^{\mathrm{s}} \left(\cos\psi \cos\theta \frac{{u}_{i,j+1,k}-{u}_{i-1,j+1,k}}{h}+\cos\psi \sin\theta \frac{{u}_{i,j+1,k}-{u}_{i,j,k}}{h}+\sin\psi \frac{{u}_{i,j+1,k}-{u}_{i,j+1,k-1}}{h}\right)\nonumber\\ & \quad ={\varepsilon}^{\mathrm{m}}\left(\cos\psi \cos\theta \frac{{u}_{i+1,j+1,k}-{u}_{i,j+1,k}}{h}+\cos\psi \sin\theta \frac{{u}_{i,j+2,k}-{u}_{i,j+1,k}}{h}+\sin\psi \frac{{u}_{i,j+1,k+1}-{u}_{i,j+1,k}}{h}\right) \end{align}$$according to the above equation, we will replace the interface point ui,j+1,k${u}_{i,j+1,k}$ in Equation (13) with the following formula15ui,j+1,k=(εi+1,j+1,kmui+1,j+1,k+εi−1,j+1,ksui−1,j+1,k)T1+(εi,j+2,kmui,j+2,k+εi,j,ksui,j,k)T2+(εi,j+1,k+1mui,j+1,k+1+εi,j+1,k−1sui,j+1,k−1)T3(εi+1,j+1,km+εi−1,j+1,ks)T1+(εi,j+2,km+εi,j,ks)T2+(εi,j+1,k+1m+εi,j+1,k−1s)T3$$\begin{equation} {u}_{i,j+1,k}=\frac{({\varepsilon}_{i+1,j+1,k}^{\text{m}}{u}_{i+1,j+1,k}+{\varepsilon}_{i-1,j+1,k}^{\text{s}}{u}_{i-1,j+1,k}){T}_{1}+({\varepsilon}_{i,j+2,k}^{\text{m}}{u}_{i,j+2,k}+{\varepsilon}_{i,j,k}^{s}{u}_{i,j,k}){T}_{2}+({\varepsilon}_{i,j+1,k+1}^{\text{m}}{u}_{i,j+1,k+1}+{\varepsilon}_{i,j+1,k-1}^{\text{s}}{u}_{i,j+1,k-1}){T}_{3}}{({\varepsilon}_{i+1,j+1,k}^{\text{m}}+{\varepsilon}_{i-1,j+1,k}^{s}){T}_{1}+({\varepsilon}_{i,j+2,k}^{\text{m}}+{\varepsilon}_{i,j,k}^{s}){T}_{2}+({\varepsilon}_{i,j+1,k+1}^{\text{m}}+{\varepsilon}_{i,j+1,k-1}^{\text{s}}){T}_{3}} \end{equation}$$with T1=cosψcosθ${T}_{1}=\cos\psi \cos\theta $, T2=cosψsinθ${T}_{2}=\cos\psi \sin\theta $, and T3=sinψ${T}_{3}=\sin\psi $.In addition, in the situation where the interface does not intersect with the grid points, the standard 7‐point difference scheme is no longer used for irregular points. We combine the mesh division and the geometry of biomolecule to identify the nodes involving interface in six directions, further deriving a 7‐point difference scheme containing the interface points. We use a second‐order accuracy approximation for ∂2u∂x2$\frac{{\partial}^{2}u}{\partial {x}^{2}}$, ∂2u∂y2$\frac{{\partial}^{2}u}{\partial {y}^{2}}$, and ∂2u∂z2$\frac{{\partial}^{2}u}{\partial {z}^{2}}$, which means that an additional point is needed in the direction involving the interface point to compensate for the loss of accuracy due to the interface. As shown in Figure 4, ui,j,k${u}_{i,j,k}$ has an interface point ub,j,k${u}_{b,j,k}$ among the difference points in the x‐direction, then an additional difference point ui+2,j,k${u}_{i+2,j,k}$ should be added. Using Taylor expansion and after algebraic operations, ∂2u∂x2$\frac{{\partial}^{2}u}{\partial {x}^{2}}$ can be discretized by ub,j,k${u}_{b,j,k}$, ui,j,k${u}_{i,j,k}$, ui+1,j,k${u}_{i+1,j,k}$, and ui+2,j,k${u}_{i+2,j,k}$, and its second‐order accuracy expression is16∂2u∂x2|ui,j,k=1hα−1α+2ui+2,j,k+2(2−α)α+1ui+1,j,k+α−3αui,j,k+6α(α+1)(α+2)ub,j,k$$\begin{eqnarray} \frac{{\partial}^{2}u}{\partial {x}^{2}}{|}_{{u}_{i,j,k}}&=&\frac{1}{h}\left[\frac{\alpha -1}{\alpha +2}{u}_{i+2,j,k}+\frac{2(2-\alpha )}{\alpha +1}{u}_{i+1,j,k}+\frac{\alpha -3}{\alpha}{u}_{i,j,k}\right.\nonumber\\ &&\left.+\ \frac{6}{\alpha (\alpha +1)(\alpha +2)}{u}_{b,j,k}\right] \end{eqnarray}$$4FigureLocal topology around irregular point (i,j,k)$(i,j,k)$. Take the x‐direction as an example, the interface crosses the x‐mesh line at the point B(b,j,k)$B(b,j,k)$ between (i,j,k)$(i,j,k)$ and (i−1,j,k)$(i-1,j,k)$, which are two irregular points. Two auxiliary lines (marked with dashed line) through B, one on the x–y plane and the other on the x–z plane, are drawn in order to derive the format for the calculation of the interface point B. Two auxiliary points (b,j,k+1)$(b,j,k+1)$ and (b,j,k+2)$(b,j,k+2)$ are placed on the auxiliary line on the x–z plane to facilitate the discretization of uzs${u}_{z}^{\text{s}}$; Additionally, two auxiliary points (b,j−1,k)$(b,j-1,k)$ and (b,j−2,k)$(b,j-2,k)$ are placed on the auxiliary line on the x–y plane to facilitate the discretization of uys${u}_{y}^{\text{s}}$. αh$\alpha h$ and βh$\beta h$ represent the interface point distances of (i,j,k)$(i,j,k)$ to the x and y directions, respectively.In order to obtain a second‐order accuracy approximation for ∂2u∂y2$\frac{{\partial}^{2}u}{\partial {y}^{2}}$ and ∂2u∂z2$\frac{{\partial}^{2}u}{\partial {z}^{2}}$, a similar process must be performed. In Figure 4, taking the x and y directions involving the difference points as an example. After discretizing the Poisson equation, we will obtain the SOR iteration format of ui,j,k${u}_{i,j,k}$17ui,j,k(n+1)=ui,j,k(n)+λαββ(α−3)+α(β−3)−2αβ6α(α+1)(α+2)ub,j,k(n+1)+2(2−α)1+αui+1,j,k(n)+α−12+αui+2,j,k(n)+6β(β+1)(β+2)ui,f,k(n+1)+2(2−β)1+βui,j+1,k(n)+β−12+βui,j+2,k(n)+ui,j,k+1(n)+ui,j,k−1(n+1)+h2ρi,j,kεi,j,k−ui,j,k(n)$$\begin{eqnarray} {u}_{i,j,k}^{(n+1)}& =&{u}_{i,j,k}^{(n)}+\lambda \left\{\frac{\alpha \beta}{\beta (\alpha -3)+\alpha (\beta -3)-2\alpha \beta}\right.\nonumber\\ && \left.\left[\frac{6}{\alpha (\alpha +1)(\alpha +2)}{u}_{b,j,k}^{(n+1)}+\frac{2(2-\alpha )}{1+\alpha}{u}_{i+1,j,k}^{(n)}+\frac{\alpha -1}{2+\alpha}{u}_{i+2,j,k}^{(n)}\right.\right.\nonumber\\ && \left.\left.+\ \frac{6}{\beta (\beta +1)(\beta +2)}{u}_{i,f,k}^{(n+1)}+\frac{2(2-\beta )}{1+\beta}{u}_{i,j+1,k}^{(n)}+\frac{\beta -1}{2+\beta}{u}_{i,j+2,k}^{(n)}\right.\right.\nonumber\\ && \left.\left.+\ {u}_{i,j,k+1}^{(n)}+ {u}_{i,j,k-1}^{(n+1)}+{h}^{2}\frac{{\rho}_{i,j,k}}{{\varepsilon}_{i,j,k}}\right]-{u}_{i,j,k}^{(n)}\right\} \end{eqnarray}$$where α and β denote the distance ratios of ui,j,k${u}_{i,j,k}$ to the interface points in the x and y directions, respectively. Furthermore, the modified finite difference approximation at irregular points maintains second‐order accuracy, provided that the interface point expressions are accurately estimated. Therefore, the interface points ub,j,k${u}_{b,j,k}$ and ui,f,k${u}_{i,f,k}$ have their special expressions in the above computational format. Inspired by Zhou,[41] we use a linear combination of the interface jump conditions to eliminate any two first‐order derivatives that are not easily approximated, which two are eliminated depends on the molecular geometry. A specific example is discussed to further illustrate the expression for the interface point B according to Figure 4. We will keep uxm${u}_{x}^{\text{m}}$ and uxs${u}_{x}^{\text{s}}$ to generate the jump condition along the x‐direction. In the y‐direction, it is more difficult to calculate uym${u}_{y}^{\text{m}}$ at (b,j,k)$(b,j,k)$ than uys${u}_{y}^{\text{s}}$. In the z‐direction, computing uzm${u}_{z}^{\text{m}}$ at (b,j,k)$(b,j,k)$ is more difficult than computing uzs${u}_{z}^{\text{s}}$. Therefore, uym${u}_{y}^{\text{m}}$ and uzm${u}_{z}^{\text{m}}$ will be eliminated by Equations (10)–(12) in interface format18ε∂u∂ξ−εmtanψ∂u∂ζ+Cη∂u∂η=Cxsuxs+Cxmuxm+Cysuys+Czsuzs$$\begin{eqnarray} &&\left[\varepsilon \frac{\partial u}{\partial \xi}\right]-{\varepsilon}^{\text{m}}\tan\psi \left[\frac{\partial u}{\partial \zeta}\right]+{C}_{\eta}\left[\frac{\partial u}{\partial \eta}\right]={C}_{x}^{\text{s}}{u}_{x}^{\text{s}}+{C}_{x}^{\text{m}}{u}_{x}^{\text{m}}\nonumber\\ &&\quad +\ {C}_{y}^{\text{s}}{u}_{y}^{\text{s}}+{C}_{z}^{\text{s}}{u}_{z}^{\text{s}} \end{eqnarray}$$where the superscripts m$\text{m}$ and s$\text{s}$ represent the molecular region and the solvent region, respectively. ux${u}_{x}$, uy${u}_{y}$, and uz${u}_{z}$ are the derivative of the electric potential, which are defined in the normal direction ξ at the dielectric interface. Then, through approximating the remaining four first‐order derivatives, we obtain the expression for ub,j,k${u}_{b,j,k}$19ub,j,k=(1−α)Cxsui−1,j,k−αCxmui,j,k+(α+α2)Cysub,j−1,k−(α+α2)Czsub,j,k+1(1−α)Cxs−αCxm+(α+α2)Cys−(α+α2)Czs$$\begin{eqnarray} {{ \def\eqcellsep{&}\begin{array}{l} {u}_{b,j,k}=\displaystyle\frac{(1-\alpha ){C}_{x}^{\text{s}}{u}_{i-1,j,k}-\alpha {C}_{x}^{\text{m}}{u}_{i,j,k}+(\alpha +{\alpha}^{2}){C}_{y}^{\text{s}}{u}_{b,j-1,k}-(\alpha +{\alpha}^{2}){C}_{z}^{\text{s}}{u}_{b,j,k+1}}{(1-\alpha ){C}_{x}^{\text{s}}-\alpha {C}_{x}^{\text{m}}+(\alpha +{\alpha}^{2}){C}_{y}^{\text{s}}-(\alpha +{\alpha}^{2}){C}_{z}^{\text{s}}} \end{array} }}\nonumber\hspace*{-10pt}\\ \end{eqnarray}$$where Cxs=εscosψcosζ+εmtanψsinψcosζ+εmsinζtanζ/cosψ${C}_{x}^{\text{s}}={\varepsilon}^{\text{s}}\cos\psi \cos\zeta +{\varepsilon}^{\text{m}}\tan\psi \sin\psi \cos\zeta +{\varepsilon}^{\text{m}}\sin\zeta \tan \zeta /\cos\psi $, Cxm=−εm/(cosψcosζ)${C}_{x}^{\text{m}}=-{\varepsilon}^{\text{m}}/(\cos\psi \cos\zeta )$, Cys=(εs+εm)cosψsinζ${C}_{y}^{\text{s}}=({\varepsilon}^{\text{s}}+{\varepsilon}^{\text{m}})\cos\psi \sin\zeta $, Czs=(εs−εm)sinψ${C}_{z}^{\text{s}}=({\varepsilon}^{\text{s}}-{\varepsilon}^{\text{m}})\sin\psi $. Note that the auxiliary points ub,j−1,k${u}_{b,j-1,k}$ and ub,j,k+1${u}_{b,j,k+1}$ are not grid nodes in the above equation, which can be calculated by interpolating points on the same horizontal line in the same region. For example, we interpolate ub,j−1,k${u}_{b,j-1,k}$ with ui−2,j−1,k${u}_{i-2,j-1,k}$, ui−1,j−1,k${u}_{i-1,j-1,k}$ and ui,j−1,k${u}_{i,j-1,k}$, and then substitute these interpolations back into the above interface point expression. In addition, these solutions can be extended to more complicated macromolecular (solute) shapes.For the discretized Poisson equation, a successive over‐relaxation iteration (SOR) stabilized solver is adopted to solve the linear system. The maximum numerical error is measured, in terms of the discrete L∞${L}_{\infty}$ norm20eu=maxi=1,2,⋯,N|uexa(ri)−unum(ri)|$$\begin{equation} {e}_{u}=\underset{i=1,2,\text{\ensuremath{\cdots}},N}{\max}|{u}_{\text{exa}}({r}_{i})-{u}_{\text{num}}({r}_{i})| \end{equation}$$where ri${r}_{i}$ is the regular grid points in solvent and solute region (use standard 7‐point differential scheme) or the irregular grid points near the interface (use new interface scheme), and the total number is denoted as N. Here, uexa(ri)${u}_{\text{exa}}({r}_{i})$ and unum(ri)${u}_{\text{num}}({r}_{i})$ are the exact solution and the numerical solution, respectively. We are interested in examining the performance of the proposed method on the interface by comparing it with the standard 7‐point differential scheme.We conducted a number of numerical experiments to examine the performance of the proposed new computational scheme near the interface, with particular attention to the order of convergence. As shown in Table 1, the grid spacing h=0.01$h=0.01$ is set as the exact solution (there is no exact solution since the right‐hand side function of the model equation is obtained from the random motion of the polymer and has no fixed form). Through solving Equations (2) and (3), we report the electrostatic potential error eu${e}_{u}$ and its convergence order regarding the grid refinement. As can be seen from Table 1, the new computational format achieves very similar results to the traditional format in terms of electric potential. The numerical results indicate that the convergence in electric potential error eu${e}_{u}$ is second order. As far as the potential value is concerned, the results obtained at coarse grid (e.g., h=1.0$h=1.0$) are very close to that at the exact solution (i.e., fine grid h=0.01$h=0.01$).1TableSolving the Poisson equation: error and order of convergence for the standard 7‐point difference scheme and the new interface schemeStandard 7‐point schemeImproved interface schemeheu${e}_{u}$Ordereu${e}_{u}$Order19.8e−047.3e−041/25.7e−042.072.6e−042.221/41.2e−041.988.9e−051.961/83.8e−052.279.6e−062.451/167.7e−062.241.1e−062.31Proposed and Optimized Local Iteration StrategyAbove interface treatment scheme lays the basic framework for solving elliptic interface problems. Solving the differential discretized system of equations using SOR iteration indeed has certain advantages. While in Brownian dynamics, one usually need to compute millions or even tens of millions of time steps, which means that the calculation of electrostatic interactions will take a lot of time. Actually, the charged particles in polymer tend to occupy only a small region in simulation box, and the spatial distribution of electrostatic potential changes little due to the small perturbation in dielectric and charge distributions at each step during such simulations. Especially when an electric field was applied, we find that the value of electric potential hardly changes within a certain time interval for grid nodes that are far from the charged polymer. This result suggests that in a short simulation time, one can save most of the iteration time if we only calculate the electric potential of grids near the charged polymer. The electric potential for other grids is fixed and updated after several steps. Overall, we can save substantial insignificant iteration time by using local iteration instead of global iteration while without much compromise in accuracy. Three important parameters need to be introduced here to ensure that the local iteration method can improve the computational efficiency while also guaranteeing a certain accuracy. As shown in Figure 5, we first locate the spatial location of the beads in polymer and enclose them with a minimal cubic basal box (A1B1C1D1${A}_{1}{B}_{1}{C}_{1}{D}_{1}$). Then, from six sides of the basal box, we extend outward the grid of Lw${L}_{\text{w}}$ units to form a “local box.” Once the “local box” exceeds the boundary of the domain ∂Ω$\partial \Omega$ in a certain direction, it will not further extend. Additionally, the parameter dstep${d}_{\text{step}}$ is defined as the step interval between the first global iteration and the next global iteration. We also set a warning value (Ld), which is considered as the threshold for updating the local iteration area. When the nearest distance between the basal box and the “local box” is equal to the warning value, the “local box” will be re‐divided based on current position of basal box.5Figure2D view of the basal box (A1B1C1D1${A}_{1}{B}_{1}{C}_{1}{D}_{1}$) and “local box” (A2B2C2D2${A}_{2}{B}_{2}{C}_{2}{D}_{2}$) under dielectric inhomogeneity. Lw denotes the length of the base box extending outward, and Ld is the warning value between the basal box and the local box. Further, A1B1C1D1${A}_{1}{B}_{1}{C}_{1}{D}_{1}$ is defined as the local basal box, and the green highlighted area indicates the banded local basal box. After several steps, the polymer is moved to virtual box (A1′B1′C1′D1′${A}_{1}^{{}^{\prime}}{B}_{1}^{{}^{\prime}}{C}_{1}^{{}^{\prime}}{D}_{1}^{{}^{\prime}}$).In fact, we can further optimize the local iteration method for polymers with special shapes. If the polymer is clustered on a certain plane (e.g., x–y, x–z, y–z), then local iteration can maximize its advantages. However, when the polymer is diagonal or dispersed in 3D space, the efficiency of the speedup is dramatically reduced. To address the adverse effects of these shape characteristics, we propose an adaptive banded local iteration strategy and expect that this approach will save a significant amount of computational expense and further improve computational efficiency. As shown in Figure 5, the banded iteration area is marked in green. The previous local iteration method is called traditional local iteration (TLI) method,[46] and the new adaptive local iteration strategy is called banded local iteration (BLI) method. Importantly, it is quite common for polymers to overlap domains in molecular simulations. In response, our algorithm is also able to accurately locate the banded local iterative regions.In our model, a new interface format is derived using a Cartesian grid, a standard finite difference scheme, and the lowest order interface jump condition. Then, the format is combined with a local iteration method to solve the elliptic equations with discontinuity coefficients and relatively smooth interfaces to simulate the movement of the polymer. We are more concerned with the performance of the local iteration algorithm under dielectric inhomogeneity rather than with how to handle interface problems with singularities. Since there has been a lot of work on the explicit interface treatment of geometric singularities, which is not the focus of this paper. Therefore, we only calculate two molecular examples with different geometrical features and classify them into two categories: molecular interfaces and grid nodes all intersect (called prism example) and partially intersect (the solute shape similar to peanut, called ellipsoid example). In different molecular configurations, we simulated three cases (as shown in Figure 6) in which the same shape of polymer was put at three different initial positions: the solvent region, the intra‐molecular region, and the junction of solvent and biomolecule, respectively. Finally, we enumerate two polymers with special shapes, and compare the numerical results of traditional local iteration and banded local iteration in two different molecular configuration computational domains.6FigureThe polymer is placed in three different positions in the computational domain, cases 1–3 corresponding to the solvent region, the internal region of the biomolecule (i.e., solute region) and the junction between solvent and molecule, respectively.Results and DiscussionIn the numerical experiments, all codes were written in Matlab 2017b and run on PC with Intel Core i5‐7500 CPU (3.40 GHz) and 8 GB RAM. At first, independent effects of three parameters (Lw${L}_{\text{w}}$, dstep${d}_{\text{step}}$, and Ld${L}_{\text{d}}$) on the computational speed and accuracy are briefly discussed. Additionally, we mainly focus on three aspects of simulations: the polymer motion at different positions in computational domain, the intersection of molecular interface with grid points (i.e., molecular shape), and the shape of the polymers. In order to make the results clearer and more universal, we will use the relative values, that is, RLw=2×Lw/L${\text{RL}}_{\text{w}}=2\ensuremath{\times{}}{L}_{\text{w}}/L$, RS=dstep/Step$\text{RS}={d}_{\text{step}}/\text{Step}$, RLd=Ld/Lw${\text{RL}}_{\text{d}}={L}_{\text{d}}/{L}_{\text{w}}$. For example, a local box is formed by extending Lw units in six directions from the center of the basal box, thus RLw${\text{RL}}_{\text{w}}$ portrays the relative size of the local box well. When RLw=100%${\text{RL}}_{\text{w}}=100\%$, whether the “local box” completely covers the simulation box depends on the initial position of the charged beads in a polymer.Benchmarked with SOR global iteration solution, both terms of accuracy and computational efficiency are investigated for the local iteration method through 3D simulation. In terms of accuracy, we report the electric potential error and the error in the ultimate position of charged particles because the final objective is to simulate the motion of charged particles, which needs an accurate calculation of the spatial potential distribution and ultimate position of the polymer. Therefore, the relative error of potential and relative position error of charged particles are defined as follows21REU=ΘRMS(Uglobal−Ulocal)/Uglobal$$\begin{equation} {\text{RE}}_{\text{U}}={\Theta}_{\text{RMS}}\left[({U}_{\text{global}}-{U}_{\text{local}})/{U}_{\text{global}}\right] \end{equation}$$22REPos=ΘRMS|Posglobal−Poslocal|/|Posglobal−Posinitial|$$\begin{equation} {\text{RE}}_{\text{Pos}}={\Theta}_{\text{RMS}}\left[|{\text{Pos}}_{\text{global}}-{\text{Pos}}_{\textit{local}}|/|{\text{Pos}}_{\text{global}}-{\text{Pos}}_{\text{initial}}|\right] \end{equation}$$where ΘRMS${\Theta}_{\text{RMS}}$ denotes the root mean square sign, the overall similarity of the polymer conformation is checked by calculating the RMS (root mean square) deviations of all charged particle positions and spatial electric potential. Uglobal${U}_{\text{global}}$ and Ulocal${U}_{\text{local}}$ represent the potential distribution at each charged particles in the global iteration and the local iteration, respectively. Posinitial${\text{Pos}}_{\text{initial}}$, Posglobal${\text{Pos}}_{\text{global}}$, and Poslocal${\text{Pos}}_{\text{local}}$ denote the 3D spatial coordinates of charged particles in the initial, global iteration and local iteration, respectively, and the displacement distance of each particle is expressed by the absolute value of the difference between them. In terms of computational speed, the CPU times of global iteration (denoted as tG) and local iteration (denoted as tL) were recorded and compared. In previous work,[46] we chose a simulation time of 10 000 time steps, which was long enough to test the stability of the simulation results. However, the global iteration simulation under dielectric inhomogeneity takes a lot of CPU time, so that we set the total step number of motion Step=1000$\text{Step}=1000$ in all simulations.In order to verify the validity of the SOR global iteration method based on a finite difference scheme (with a grid spacing h=1$h=1$), we compare it with existing conventional algorithms for calculating electrostatic forces. In particular, the Coulomb formulation, an exact method for solving electrostatic interactions, will be used to evaluate the validity of our proposed numerical method. During the simulation, we use the global iteration method and Coulomb's formula in the section of solving for the electrostatic force respectively, while keeping the settings consistent in all other sections, and then record the CPU time and numerical results such as particle position and spatial potential. We display the CPU time as a function of the number of charged particles (Num$\text{Num}$) for these two methods in Figure 7. The comparison results indicates that as the number of charged particles increases, the CPU time of these two methods both exhibits an increase. However, the SOR method is different from the Coulomb formula for solving the electrostatic force directly, which is insensitive to the increases in Num$\text{Num}$. Obviously, as the number of particles increases, the SOR global iteration only increases the workload of interpolating more charges to nearby nodes, which can even be ignored compared with time‐consuming iteration loop. On the contrary, with increasing number of particles, the CPU time of the Coulomb formula increases significantly. The numerical results indicate that the SOR method is more suitable for large‐scale charged particle system. This conclusion is supported by the results in the inset of Figure 7, where the computational efficiency of the global iteration gradually outperforms the Coulomb formula when Num>400$\text{Num}&gt;400$. As Num$\text{Num}$ continues to increase, the computational efficiency advantage of the global iteration becomes more pronounced, and it is nine times more efficient than the Coulomb formula when Num=1500$\text{Num}=1500$. In general, the calculation complexity of the Coulomb formula is O(N2)$O({N}^{2})$, while the calculation complexity of a single iteration of the SOR global iteration method with optimal relaxation parameters is between O(NlnN)$O(N\text{ln}N)$ and O(NlgN)$O(N\text{lg}N)$. Theoretically, there is a significant difference in computational effort between the two methods, and Figure 7 should show a corresponding CPU time difference. In fact, the iteration method requires a loop of iterations until the convergence criterion is satisfied. We assume that a complete loop process requires the implementation of n iterations, then the total computational effort to compute the single electrostatic force needs to be multiplied by a factor n. In addition, using the results of the Coulomb formula as a benchmark, Table 2 shows the root mean square deviation of the charged particle position and electric potential for global iteration method. We find that as NUM increases, both REU${\text{RE}}_{\mathrm{U}}$ and REPos${\text{RE}}_{\text{Pos}}$ increase. In the REU${\text{RE}}_{\mathrm{U}}$ term, as Num$\text{Num}$ increases from 100 to 1500, the error increases by only about 1.12e−3. In the REPos${\text{RE}}_{\text{Pos}}$ term, although the error increases from 6.66e−4 to 1.05e−3, we can observe that the increase in error is not significant when Num>1000$\text{Num}&gt;1000$. Overall, the two errors are controlled to be below 3.16e−3 and 1.05e−3, respectively, both in the 1e−3 magnitude, which is an acceptable error range for the experiment. Therefore, it is a feasible scheme to apply a global SOR iteration method based on a finite difference format to solve the electrostatic force and thus simulate charged polymer motion.2TableThe root mean square deviation of the charged particle position and electric potential10030050010001500REU${\text{RE}}_{\text{U}}$2.04e−32.42e−32.67e−33.02e−33.16e−3REPos${\text{RE}}_{\text{Pos}}$6.66e−48.20e−49.01e−49.92e−41.05e−37FigureCPU time t as a function of the number of charged particles (Num$\text{Num}$). The inset shows the relative percentage of CPU time for the global iteration method and the Coulomb formula at different numbers of particles.In Figures 8 and 9, parts (a) and (b), we display CPU time t and the relative error as a function of RLw${\text{RL}}_{\mathrm{w}}$ (the length ratio of “local box” to “global box”), respectively, where the relative error is divided into two parts with the relative error of electrostatic potential (REU${\text{RE}}_{\mathrm{U}}$) and the relative error of the final position of the beads (REPos${\text{RE}}_{\text{Pos}}$). Meanwhile, the dimension of simulation box was set as L=80$L=80$ with grid spacing h=Δx=Δy=Δz=1$h= \Delta x= \Delta y= \Delta z=1$, which means that the total number of grid nodes is 81×81×81$81\ensuremath{\times{}}81\ensuremath{\times{}}81$. The number of charged particles in polymer is set as Num=500$\text{Num}=500$ (the FDM method has significant advantages for handling electrostatic interactions containing large numbers of charged particles), the relative distance between charged beads and boundary RLd=80%${\text{RL}}_{\mathrm{d}}=80\%$. And the value of relative global iteration step interval RS$\text{RS}$ is set as 10% and 25% in simulations to explore the impact of update frequency of global iteration on the computational speed and accuracy. After verification, we find that the increasing RLd${\text{RL}}_{\mathrm{d}}$ is beneficial to reduce REU${\text{RE}}_{\mathrm{U}}$, while the CPU time only slightly increases at the same time. Therefore, we uniformly set RLd${\text{RL}}_{\mathrm{d}}$ in all the examples to 80%. Under these conditions, we calculate tG, tL, REU${\text{RE}}_{\mathrm{U}}$, and REPos${\text{RE}}_{\text{Pos}}$ of the two numerical examples of prismatic and ellipsoidal biomolecules, in which the numerical results of different examples are represented by different colored lines. Based on these two examples, we constructed three cases as shown in Figures 8 and 9.8FigureThe case of polymer in solvent region: a) CPU time t and b) the relative error of electrostatic potential REU${\text{RE}}_{\text{U}}$ as a function of the length ratio between “local box” and “global box” (RLw${\text{RL}}_{\text{w}}$). The relative distance between charged polymer and boundary is set as RLd${\text{RL}}_{\text{d}}$ = 80%. a) RS$\text{a) RS}$ represents relative global iteration step interval, and tG represents the global iteration time. b) The relative error of the final position of the charged particles (REPos${\text{RE}}_{\text{Pos}}$) as a function of the relative size of local box (RLw${\text{RL}}_{\mathrm{w}}$).9FigureThe case of polymer in internal region of the solute: a) CPU time t and b) the relative error of electrostatic potential REU${\text{RE}}_{\text{U}}$ as a function of the length ratio between “local box” and “global box” (RLw${\text{RL}}_{\text{w}}$). The relative distance between charged polymer and boundary is set as RLd${\text{RL}}_{\text{d}}$ = 80%. RS$\text{RS}$ represents relative global iteration step interval, and tG in (a) represents the global iteration time. The inset of (b) shows the relative error of the final position of the charged particles (REPos${\text{RE}}_{\text{Pos}}$) as a function of the relative size of local box (RLw${\text{RL}}_{\mathrm{w}}$).Initially, as shown in Figure 6, we put the polymer chain in solvent area and the charged particles of the polymer chain are under electrostatic interaction, elastic force and repulsive force to obtain an equilibrium conformation (case 1). For the solvent region, real‐time updating of the dielectric value everywhere in the solvent area may spend more computing costs because the dielectric value is not a constant but as a function of the electric field strength, which implies that the calculation of electrostatic interaction will consume more CPU time. However, these calculation costs are necessary since the dielectric response can accurately capture certain characteristics of the biological process. As shown in Figure 8a, with the decrease of RLw${\text{RL}}_{\mathrm{w}}$, the CPU time decreases. Obviously, reducing the local iteration domain means saving a part of the iteration time. In addition, the computational expense of iteration increases with the decrease of RS$\text{RS}$. In other words, as the RS$\textit{RS}$ decreases, the frequency of global iteration increases and the cumulative error decreases, which further leads to the increase of total CPU time and the decrease of REU${\text{RE}}_{\mathrm{U}}$. These results imply that a smaller RLw${\text{RL}}_{\mathrm{w}}$ and a larger RS$\text{RS}$ are beneficial to improve calculation efficiency. Under the same conditions, the ellipsoid group reported significantly more CPU time than the prismatic group. There are two main reasons: in the program of ellipsoid example, it is necessary to open up an independent space to store the information of interface values that do not intersect with nodes. Otherwise, the calculation format of these interface points is more complicated. As the local box continually expands, we find that the tL of both prismatic and ellipsoid examples are approximately equal to tG when RLw=90%${\text{RL}}_{\text{w}}=90\%$, because the smallest basal box occupies part of the computational domain. It's worth noting that during RLw${\text{RL}}_{\mathrm{w}}$ increasing from 35% to 60%, a significant increase in the CPU time was recorded. According to the preliminary analysis, when RLw<35%${\text{RL}}_{\text{w}}&lt;35\%$, the local iteration is only performed in the solvent area. With the successive increase of RLw${\text{RL}}_{\mathrm{w}}$, the local iteration range further covers the molecular (solute) area. As a result, the computational complexity of nodes near the molecular interface significantly increases the CPU time. When RLw>60%${\text{RL}}_{\text{w}}&gt;60\%$, the local iteration area completely covers the entire biomolecule, and this trend stabilizes. Besides, the results of REPos${\text{RE}}_{\text{Pos}}$ illustrate the relative error of the polymer translocation, where the REPos${\text{RE}}_{\text{Pos}}$ and REU${\text{RE}}_{\mathrm{U}}$ share a similar change trend. As shown in Figure 8b, REU${\text{RE}}_{\mathrm{U}}$ and REPos${\text{RE}}_{\text{Pos}}$ decrease with the increase of RLw${\text{RL}}_{\mathrm{w}}$. These results indicate that choosing a larger RLw${\text{RL}}_{\mathrm{w}}$ gets closer to the global solution. We find that the REU${\text{RE}}_{\mathrm{U}}$ has dropped by an order of magnitude when RLw≈18%${\text{RL}}_{\text{w}}\approx 18\%$. Moreover, we also find that there is a difference in REU${\text{RE}}_{\mathrm{U}}$ between the two examples. The reason for this phenomenon is that the ellipsoid calculation example compensate the loss of accuracy by repeatedly using the interface conditions and adding iteration nodes near the interface. Importantly, when RLw=7.5%${\text{RL}}_{\text{w}}=7.5\%$ and RS=25%$\text{RS}=25\%$, the single most striking observation to emerge from the data comparison was the calculation speed of the two examples is increased by 20 and 17 times respectively, and the REU${\text{RE}}_{\mathrm{U}}$ is controlled below 3.2e−4.Besides, we also put the polymer chain in solute area and observed its motion and the performance of local iteration method (case 2). The results in Figure 9a show that CPU time decreases with the decrease of RLw${\text{RL}}_{\mathrm{w}}$. At the same time, a decreasing RS$\text{RS}$ also induces the decrease of CPU time under the same RLw${\text{RL}}_{\mathrm{w}}$. As shown in Figure 9b, with the increase of RLw${\text{RL}}_{\mathrm{w}}$, the REU${\text{RE}}_{\mathrm{U}}$ first decreases rapidly, then remains unchanged. In the inset of Figure 9b shows the relative position error (REPos${\text{RE}}_{\text{Pos}}$) as a function of RLw${\text{RL}}_{\mathrm{w}}$, which decreases as RLw${\text{RL}}_{\mathrm{w}}$ increases. In fact, if RLw${\text{RL}}_{\mathrm{w}}$ is large enough the error will be almost negligible. These situations are very similar to case 1, which is exactly what we expected. Consequently, we can get a conclusion that the local iteration algorithm is not affected by the spatial position of the polymer, which reflects its stability. Interestingly, we find that when RLw<50%${\text{RL}}_{\text{w}}&lt;50\%$, cases 1 and 2 presented large differences in CPU time. Theoretically, as long as the local box does not contact with the molecular interface, the node values in the local iterative domain are calculated by the standard 7‐point center difference scheme. However, because the molecular volume is only 1.5% of the computational domain and the large number of charged beads in polymer located inside the biomolecule has occupied almost the entire biomolecule (Num=500$\text{Num}=500$), it drives the local box to contact the molecular interface prematurely (in case 2). Therefore, the computational amount of case 2 is much higher than that of case 1 even in the initial stage of RLw${\text{RL}}_{\mathrm{w}}$ increase. In a way, the acceleration effect of local iteration in case 2 is not as good as that in cases 1 and 3.In fact, as shown in Figure 6, setting the polymer chain at the junction of solute and solvent area is another important situation for exploring the molecular behavior on singlemolecule level, which takes charged polymer intersected vesicles as prototype (case 3). Fast calculation of electrostatic potential provides convenience for applications of vesicles in pharmaceutical and physiological processes. As shown in Figure 10a, with the increase of RLw${\text{RL}}_{\mathrm{w}}$, the CPU time first increases rapidly, and then the rate of increase becomes slower. Similarly, Figure 10b shows that REU${\text{RE}}_{\mathrm{U}}$ decreases with the decrease in RS$\text{RS}$ for the same RLw${\text{RL}}_{\mathrm{w}}$, which indicates that increasing the frequency of global iteration can significantly reduce the accumulation of errors. When RLw=37%${\text{RL}}_{\text{w}}=37\%$, the CPU time of the ellipsoid example gradually draws a significant gap with the prism example. It is also at this stage that the growth in CPU time becomes stable. Overall, the result of case 3 is very similar to the other two cases, which is consistent with our expectations. While, the CPU time of case 3 first increases rapidly and then increases slowly, opposite to cases 1. Through analysis, we find that this has a great relationship with the region contained in local iteration box. When the polymer spans two regions, no matter how small the local box, the initial calculation will always involve interface nodes, leading to more difficult iterative solution. Additionally, in the process of polymer movement, free moving ions are distributed on the left, while fixed charges are distributed in the solute region on the right, which further increases the computational cost. Even though the calculation speed decreases slightly, when RLw=7.5%${\text{RL}}_{\text{w}}=7.5\%$ and RS=25%$\text{RS}=25\%$, the computation efficiency still increased by 8.2 times compared to the global iteration in the ellipsoid numerical example.10FigureThe case where the polymer is located at the solute‐solvent junction: a) CPU time t and b) the relative error of electrostatic potential REU${\text{RE}}_{\text{U}}$ as a function of the length ratio between “local box” and “global box” (RLw${\text{RL}}_{\text{w}}$). The relative distance between charged polymer and boundary is set as RLd${\text{RL}}_{\text{d}}$ = 80%. RS$\text{RS}$ represents relative global iteration step interval, and tG in (a) represents the global iteration time. The inset of (b) shows the relative error of the final position of the charged particles (REPos${\text{RE}}_{\text{Pos}}$) as a function of the relative size of local box (RLw${\text{RL}}_{\mathrm{w}}$).In order to further investigate our analysis in detail, we combine the results of the three position cases of the two calculation examples and shown them in Figures 11a and 11b, respectively. The left axis and right axis denote CPU time t and the relative potential error (REU${\text{RE}}_{\mathrm{U}}$) as a function of RLw${\text{RL}}_{\mathrm{w}}$, respectively. First, with the increase of RLw${\text{RL}}_{\mathrm{w}}$, the CPU time increases. Meanwhile, the decrease of RS$\text{RS}$ is also an important factor for the increase of CPU time. We find that tG is almost equal for the three cases in the same numerical example (e.g., ellipsoid example is about 41000s$41\ 000\ \text{s}$), this indicates that the global iteration time is independent of the polymer placement and is only related to the molecular conformation. Furthermore, by comparing the simulation results of Figure 11a,b, it clearly implies that tG of the ellipsoid example is always large than the prismatic example under same conditions. Accordingly, it is further verified that the CPU time of the global iteration only depends on the calculation format and has no relationship with the position of the polymer. We know that the potential value around the molecular interface is the major calculation consumption in the dielectric inhomogeneous system. It also explains that the CPU time of the dielectric inhomogeneity is much more than the dielectric homogeneity in global iteration. When the curves of all cases are plot in figure, we can clearly observe some differences, the overall trend of the connection between RLw${\text{RL}}_{\mathrm{w}}$ and CPU time is very similar, but the computational efficiency changes at different RLw${\text{RL}}_{\mathrm{w}}$ (note that the difference mentioned above is not related to the numerical example, but only to the spatial location of the polymer, thus the two diagrams are described uniformly): i) for case 1, the computational efficiency is greatly improved in the initial stage, and the CPU time increases slowly. When RLw>35%${\text{RL}}_{\text{w}}&gt;35\%$, the CPU time increases suddenly and rapidly, ii) in case 2, the rate of the change in computational speed is different from that of case 1 (RLw=25%${\text{RL}}_{\text{w}}=25\%$), because the solute region is smaller than the solvent region. When the local iteration domain is expanding, case 2 covers the interface area earlier than case 1. Undoubtedly, the calculation of the potential value around the interface occupies most of the CPU time. Hence, the CPU time grows rapidly when RLw=7.5--25%${\text{RL}}_{\text{w}}=7.5\text{--}25\%$ and then increases slowly. iii) As expectation, this is further verified in case 3, which shows a similar phenomenon to case 2. In the RLw<37%${\text{RL}}_{\text{w}}&lt;37\%$ stage, the CPU time increases rapidly, which is not only related to the current position of the polymer (the junction between solute and solvent), but also has a certain relationship with the direction of polymer movement. The polymer moves toward the interior of the biomolecule under the action of the applied electric field force and the electrostatic force in the electrolyte solution, so that the speed of the local box covering the molecular region becomes faster. Until the critical value 37% is exceeded, the molecular region is completely covered. As RLw${\text{RL}}_{\mathrm{w}}$ increases continuously, the local box expands the number of nodes in the solvent region, so that the increase rate of CPU time gradually slows down.11Figurea,b) The left axis and right axis denote CPU time t and the relative potential error REU${\text{RE}}_{\mathrm{U}}$ as a function of the relative size of local box (RLw${\text{RL}}_{\mathrm{w}}$), respectively. (a) Calculates an ellipsoid numerical example, (b) calculates a prism numerical example. Marked the CPU time and error for the equal conditions with a same color line, where the solid represents the CPU time and the hollow represents the REU${\text{RE}}_{\mathrm{U}}$, and the intersection point is defined as the optimal parameter point for this group (mark with solid circle).Additionally, as shown in Figure 11a,b, with the increase of RLw${\text{RL}}_{\mathrm{w}}$ and the decrease of RS$\text{RS}$, the value of REU${\text{RE}}_{\mathrm{U}}$ decreases significantly. However, the RLw${\text{RL}}_{\text{w}}$–REU${\text{RE}}_{\text{U}}$ function curves of the three cases seem to have some differences. In the discussion of interface point calculation format, we construct a high‐precision discrete difference scheme via adding differential nodes. Therefore, we consider that the major error originates from the solute and solvent region (7‐point center difference scheme) rather than the region around the interface (high‐precision difference scheme). And we also infer that this is opposite to the change trend of CPU time. Numerical results prove that the initial error of case 3 is the largest, because the polymer is originally located at the solute–solvent junction. The electric potential in the region of major error source is not updated for a period of time, which will lead to an error propagation effect and increase REU${\text{RE}}_{\mathrm{U}}$. For cases 1 and 2, the smaller RLw${\text{RL}}_{\mathrm{w}}$ makes the local box only work in the solvent or solute region, so that error propagation in the solvent or solute region is blocked. Theoretically, according to current simulation data, we can explore the optimal parameter selection for each numerical example, that is, the intersection of the same set of curves, which is marked as solid circle symbol in Figure 11. Obviously, the optimum point does not have the best acceleration effect, but it is the most cost‐effective situation, which is a good balance between computational efficiency and accuracy. For example, in case 1 of the ellipsoid example (Figure 11a), the optimal parameters of the ellipsoid example are RLw=28%${\text{RL}}_{\text{w}}=28\%$, RS=25%$\text{RS}=25\%$, RLd=80%${\text{RL}}_{\text{d}}=80\%$. Certainly, if test more parameter combinations, you may get a better parameter group. In addition, in both Figure 11a,b under the same case, even different RS$\text{RS}$ are selected, the optimal parameters almost remain within a certain interval. For instance, in case 1 of Figure 11b, RLw=15--20%${\text{RL}}_{\text{w}}=15\text{--}20\%$ is the optimal choice, while in case 3, the optimal RLw${\text{RL}}_{\mathrm{w}}$ takes the value between 19% and 23%.As noted above, we now examine the performance of the local iteration method for simulating the motion of charged polymer under two different shape of solute molecular models. Generally, one find that local iteration performs well when the type of the polymer is clustered or planar. However, once the shape of the polymer is scattered or diagonal, the effect of local iteration will be greatly reduced. This causes certain challenges to the universality of our algorithm. In addition, the computational difficulties caused by the special iteration scheme of the nodes around molecular interface bring another challenge to the total CPU time. Therefore, in order to improve the local iteration method, we learn from the original optimization ideas of the traditional local iteration (TLI) to further reduce unimportant iterations time. We do this by constructing a compact base box that tightly wraps the polymer. That is, the configuration of the initial base box is the same as the polymer shape, which minimizes the unnecessary iteration area. Another advantage of this approach is that even the shape of the polymer is constantly changing as it moves, the adaptive banded local iteration (BLI) can locate the current shape of the polymer in real time, and expand RLw${\text{RL}}_{\mathrm{w}}$ units in six directions to form local box. Note that we cancel the warning value RLd${\text{RL}}_{\mathrm{d}}$ in the BLI method, because many numerical experiments show that it does not spend much time to re‐divide the local box frequently.Finally, we simulate two polymers of different shapes under two models with different molecular structures, using a banded local iterative approach. Figure 12a,b display the CPU time (t) as a function of the relative length of the local box (RLw${\text{RL}}_{\mathrm{w}}$) at the relative global iteration step interval RS=25%$\text{RS}=25\%$. With the decrease in RLw${\text{RL}}_{\mathrm{w}}$, the iteration domain reduce, which results in a decrease in CPU time. On the other hand, REU${\text{RE}}_{\mathrm{U}}$ decreases with increasing RLw${\text{RL}}_{\mathrm{w}}$. In terms of error, when RLw<25%${\text{RL}}_{\text{w}}&lt;25\%$, the REU${\text{RE}}_{\mathrm{U}}$ of the ellipsoid numerical example is always lower than that of the prism example. This result coincides with our aim to construct a high‐precision differential scheme by adding differential nodes at the disjoint interface points. Surprisingly, compared to global iteration, BLI method has better acceleration efficiency than TLI method, accompanied by a decrease in accuracy. As shown in Figure 12a, we first consider a diagonal type of polymer (3D), in which the position of the polymer is at the intersection of the solute and solvent. The effects of BLI method on CPU time and accuracy in different numerical examples are evaluated. In this situation, we only increase the simulation efficiency of the TLI method by about three times compared to that of the global iteration due to the limitations of its strategy and polymer shape. As from the inset of Figure 12a, we can imagine the diagonal polymer as a 3D diagonal matrix with a large number of “zero elements,” and try to get faster computation speed by cutting these unnecessary nodes. In Figure 12a, we observe that the computation speed of BLI method is ≈25 times faster than the global iteration (tG${t}_{G}$ is about 34000--42000s$34000\text{--}42000\ \text{s}$). Now, the TLI method takes almost nine times as much CPU time as the BLI method. We also focus on the variation in error. Under the same conditions, the error of the BLI method is higher than that of the TLI method, but is still in the same order of magnitude (10−4), which demonstrates that the BLI has a higher benefit. Moreover, we wanted to verify the performance of the BLI method with a polymer of scattered type, which is also put in the solute region (similar to case 2). We can see the powerful acceleration advantage of BLI method in Figure 12b. When RLw=7.5%${\text{RL}}_{\text{w}}=7.5\%$, the computational efficiency improves 29 times over the TLI method, reaching a peak. At the same time, REU${\text{RE}}_{\mathrm{U}}$ increases from 2.48e−4 to 5.15e−4. We find that the computational efficiency of BLI method can be increased by at least 3.3 times, while the TLI method can only improve the computational efficiency at most 2.8 times, which shows that the BLI method indeed improves the computational efficiency. Otherwise, the results demonstrate that BLI method not only improve the simulation efficiency for polymers with special shapes, but also for case 2 in Figure 9 (where the performance of TLI method is not excellent, only 5.6 times at most). By reducing the iteration area, the local box contacts the interface more slowly and results in improved simulation efficiency.12Figurea,b) The left axis and right axis denote CPU time t and the relative potential error REU${\text{RE}}_{\mathrm{U}}$ as a function of the relative size of local box (RLw${\text{RL}}_{\mathrm{w}}$), respectively. Marked the CPU time and error for the equal conditions with a same color line, where the solid represents the CPU time and the hollow represents the REU${\text{RE}}_{\mathrm{U}}$, and the intersection point is defined as the optimal parameter point for this group (mark with solid circle). The inset is sketched a simple 3D coarse‐grained model of the polymer. (a) simulates a diagonal type polymer, and set it at the junction of solute and solvent, (b) simulates a scattered type polymer, and put it on the solute region.We further extract the optimal parameters from two examples, which can be seen from Figure 12a,b. The current numerical results show that the optimal banded local box size is RLw=34%${\text{RL}}_{\text{w}}=34\%$ in the prismatic example and RLw=32%${\text{RL}}_{\text{w}}=32\%$ in the ellipsoid example for simulating the motion of diagonal polymers. Furthermore, when the shape of the polymer is scattered, the optimal band local box relative size is selected as 34% and 30%, respectively. By comparing and analyzing Figure 12a,b, we find that for the polymer of different shapes, if use the same local iteration algorithm, the optimal parameter ranges are very similar. When the TLI method is used, the optimal RLw${\text{RL}}_{\mathrm{w}}$ is ≈30%$\approx 30\%$ to 34%. However, for the BLI method, setting RLw=10--15%${\text{RL}}_{\mathrm{w}}=10\text{--}15\%$ can achieve the best balance between computational precision and computational cost. This simulation result implies that the distribution interval of the optimal parameters appeared to be unaffected by the shape of the polymer and is more related to the type of local iteration method.ConclusionsIn this work, using a hybrid simulation approach that combines a finite difference method and a Brownian dynamics model (FD–BD) for the polymer, we study the diffusion of charged polymers in the dielectric inhomogeneity system. We believe that further optimization of the hybrid FE‐BD approach might include the following. Perhaps the main difficulty is need to solve an equation based on a discontinuous medium in the entire domain, that is, the Poisson equation with variable dielectric coefficient. We propose a new computational scheme based on FDM with second‐order accuracy and effectively handles with the interface problem posed by such equations, achieving a smooth transition of the potential value from one region to another. The numerical results show that the improved scheme has a second order of convergence. In addition, we validate the validity of the global SOR iteration method based on the improved interface scheme. It is an excellent numerical tool for simulating charged polymer motion since its computational efficiency is better than the traditional Coulomb formula. However, in simulations that involve millions or even tens of millions of time steps, such computational efficiency is far from adequate. In order to avoid the iterative section to occupy too much CPU resources, the local iteration method is used to accelerate this process. In the simulations, we modeled numerical examples of two different biomolecule configurations (ellipsoid and prism) and put the polymer of same shape in the solvent region, the solute (biomolecule) region, and the solute–solvent interface, respectively. As expected, we find that both the decrease in RLw${\text{RL}}_{\mathrm{w}}$ and the increase in RS$\text{RS}$ can significantly reduce CPU time and thus improve computational efficiency, while the corresponding cost is the increase in error. We also find that although the overall trends of computational efficiency is very similar for the three cases, they exhibit slight differences in the inflection points of the magnitude changes: as local box continues to expand, computational efficiency is expressed as the solvent region > the interface > the internal region of solute in the early stages, while in the late stage, the internal region of solute > the interface > the solvent region. In general, the variation in the magnitude of computational efficiency depends on the rate at which the local box covers the molecular interface. In addition, our simulation results demonstrate that the error of local iteration method increases as the distance between the polymer and the molecular interface decreases.We further discuss the optimal parameter sets in the numerical experiments of this paper, which are chosen by combining both computational speed and accuracy for the sake of reasonableness. For example, RLw=28%${\text{RL}}_{\text{w}}=28\%$, RS=25%$\text{RS}=25\%$, RLd=80%${\text{RL}}_{\text{d}}=80\%$ are the optimal set of parameters for case 1 in the ellipsoid numerical example, where the computational efficiency is improved by 4.5 times compared to the global iteration, and REU${\text{RE}}_{\mathrm{U}}$ is controlled below 9.5e−5. Importantly, in all cases of the ellipsoidal example, RLw${\text{RL}}_{\mathrm{w}}$ is basically concentrated between 10% and 28% to ensure the best results of simulation, while in the prismatic example it takes values from 12% to 23%. Accordingly, the optimal parameters are determined based on different examples. In this way, the local iteration method effectively solves the computational complexity of global iteration in the dielectric inhomogeneous system, which saves a significant amount of CPU time.However, we observe that the traditional local iteration (TLI) method does not have significant acceleration effect for the polymers of two particular shapes. Therefore, we propose the banded local iteration (BLI) method, which is better adapted to the effects of deformation of the polymer during motion. Our results indicate that the BLI method additionally improves the computational efficiency over the TLI method. In particular, when the polymer shape is scattered, the BLI method improves the computational efficiency by a factor of 29 over the global iteration, which is much higher than that of the TLI method (4.5 times). We hope this work can be helpful for accelerating the solution of electrostatic interaction forces under dielectric inhomogeneity, and provide theoretical guidance for the optimization of iterations for polymer chains with different geometries in large electrostatic systems (containing a large number of charged particles or huge computational domain). Furthermore, we intend to study the performance comparison of local iteration method and P3M, Ewald summation in high concentration salt solution in our further work. In future studies, it would be very fascinating to understand the motion of charged particles under the dielectric inhomogeneous system from a different angle, for example, by searching the energy minimum of the system to achieve a mutual equilibrium, which will help us to solve more complex molecular with singularities.AcknowledgementsThis work is supported by the National Natural Science Foundation of China (Grant No. 12071363 and 11701433), the Science and Technology Major Project of Hubei Province under Grant 2021AAA010, and the National Key Research and Development Program of China (No. 2020YFA0714200).Conflict of InterestThe authors declare no conflict of interest.Data Availability StatementThe data that support the findings of this study are available from the corresponding author upon reasonable request.D. L. Ermak, J. A. McCammon, J. Chem. Phys. 1978, 69, 1352.B. A. Luty, J. A. McCammon, H. X. Zhou, J. Chem. Phys. 1992, 97, 5682.H. X. Zhou, J. Chem. Phys. 1990, 92, 3092.S. R. McGuffee, A. H. Elcock, Plos Comput. Biol. 2010, 6, e1000694.S. H. Northrup, S. A. Allison, J. A. McCammon, J. Chem. Phys. 1984, 80, 1517.S. H. Northrup, J. O. Boles, J. C. Reynolds, J. Phys. Chem. 1987, 91, 5991.S. Lu, A. P. Michailova, J. J. Saucerman, Y. Cheng, Z. Yu, T. H. Kaiser, W. W. Li, R. E. Bank, M. J. Holst, J. A. McCammon, IEEE Eng. Med. 2009, 28, 46.I. Nakamura, Z. G. Wang, Soft Matter. 2012, 8, 9356.I. Gunkel, T. Thurn‐Albrecht, Macromolecules 2012, 45, 283.C.‐L. Ren, I. Nakamura, Z.‐G. Wang, Macromolecules 2016, 49, 425.A. A. Teran, N. P. Balsara, J. Phys. Chem. B 2014, 118, 4.A. Abrashkin, D. Andelman, H. Orland, Phys. Rev. Lett. 2007, 99, 077801.A. Hildebrandt, R. Blossey, S. Rjasanow, O. Kohlbacher, H.‐P. Lenhof, Phys. Rev. Lett. 2004, 93, 108104.B. Lu, Y. Zhou, M. Holst, J. McCammon, Commun. Comput. Phys. 2008, 3, 973.A. Y. Toukmaji, J. A. Board Jr, Comput. Phys. Commun. 1996, 95, 73.M. Deserno, C. Holm, J. Chem. Phys. 1998, 109, 7694.T. Darden, D. York, L. Pedersen, J. Chem. Phys. 1993, 98, 10089.V. Ballenegger, J. J. Cerda, C. Holm, J. Chem. Theory Comput. 2012, 8, 936.H. M. Antia, Numerical Methods for Scientists and Engineers, Third Edition, Springer, 2012.G. D. Smith, Numerical Solution of Partial Differential Equations: Finite Difference Methods, Oxford University Press, Oxford 1985.D. Young, T. Am. Math. Soc. 1954, 76, 92.J. Zhang, J. Comput. Phys. 1998, 143, 449.X. I. A. Yang, R. Mittal, J. Comput. Phys. 2014, 274, 695.T. Fukuchi, AIP Adv. 2014, 4, 055312.W. F. Spotz, G. F. Carey, Numer. Meth., Part D E 2015, 12, 235.G. Sutmann, B. Steffen, J. Comput. Appl. Math. 2006, 187, 142.H. Wang, Y. Zhang, X. Ma, J. Qiu, Y. Liang, Comput. Math. Appl. 2016, 71, 1843.J. E. Adsuara, I. Cordero‐Carrion, P. Cerda‐Duran, M. A. Aloy, J. Comput. Phys. 2016, 321, 369.G. A. Huber, Y. Miao, S. Zhou, B. Li, J. A. McCammon, J. Chem. Phys. 2016, 144.P. Bauler, G. A. Huber, J. A. McCammon, J. Chem. Phys. 2012, 136, 164107.R. Egan, F. Gibou, J. Comput. Phys. 2017, 346, 71.J. D. Towers, J. Comput. Phys. 2018, 361, 424.M. K. Gilson, M. E. Davis, B. A. Luty, J. A. McCammon, J. Phys. Chem. 1993, 97, 3591.Z. Zhou, P. Payne, M. Vasquez, N. Kuhn, M. Levitt, J. Comput. Chem. 1996, 17, 1344.I.‐L. Chern, J.‐G. Liu, W.‐C. Wang, Methods. Appl. Anal. 2003, 10, 309.H. Li, B. Lu, J. Chem. Phys. 2014, 141, 07B609_1.X. Liu, Y. Qiao, B. Lu, SIAM J. Appl. Math. 2018, 78, 1131.D. Xie, J. Comput. Phys. 2014, 275, 294.Y. Zhou, S. Zhao, M. Feig, G.‐W. Wei, J. Comput. Phys. 2006, 213, 1.Y. Zhou, G.‐W. Wei, J. Comput. Phys. 2006, 219, 228.Y. Zhou, M. Feig, G.‐W. Wei, J. Comput. Chem. 2008, 29, 87.K. Xia, M. Zhan, G.‐W. Wei, J. Comput. Phys. 2011, 230, 4588.W. Geng, S. Zhao, J. Comput. Phys. 2017, 351, 25.H. Feng, G. Long, S. Zhao, J. Comput. Appl. 2019, 361, 426.A. Lee, W. Geng, S. Zhao, J. Comput. Phys. 2021, 426, 109958.H. Lin, Z.‐T. Lei, M.‐M. Ding, H.‐J. Wang, T.‐F. Shi, Chinses. J. Polym. Sci. 2020, 38, 665.F. Booth, J. Chem. Phys. 1951, 19, 391.H. Warner, Ind. Eng. Chem. Fund. 1972, 11, 379.D. Frenkel, B. Smit, Understanding Molecular Simulation, Academic Press, New York 2002.M. Allen, D. Tildesley, Computer Simulation of Molecular Liquids, Clarendon Press, Oxford 1987.G. C. Maitland, Intermolecular Forces: Their Origin and Determination, Oxford University Press, Oxford 1981. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Advanced Theory and Simulations Wiley

An Adaptive Local Iterative Method for Fast Calculation of Electrostatic Interactions between Charged Polymers in Dielectric Inhomogeneous System

Download PDF
14 pages

Loading next page...
 
/lp/wiley/an-adaptive-local-iterative-method-for-fast-calculation-of-KMtaitwiu0

References (41)

Publisher
Wiley
Copyright
© 2023 Wiley‐VCH GmbH
eISSN
2513-0390
DOI
10.1002/adts.202200776
Publisher site
See Article on Publisher Site

Abstract

IntroductionElectrostatic interactions between the charged polymer and biomolecule‐solvent environment play an important role in many biological processes, since many diffusing molecules of interest in biological systems carry charges.[1–7] It is well known that when the dielectric value in an electrolyte solution is represented by the segmental dielectric function, this belongs to the research category of dielectric inhomogeneity, which helps us understand the structure, function and dynamics of biomolecules, etc. In the experiment investigating phase separation of block copolymers, when the dielectric inhomogeneity is considered, the ions added to the solution tend to be enriched in the high dielectric constant components, thereby promoting phase separation, which is beyond the original cognition.[8–11] Furthermore, a particularly interesting aspect of this work is the consideration of the dependence of the dielectric coefficient on the electric field strength in electrolyte solutions, since the solvent dipole is highly oriented under a large electric field and thus affects the dielectric distribution.[12,13] The heterogeneity of the dielectric coefficient dramatically affects the structure of the bilayer region and charge transfer. Therefore, understanding the dielectric properties of solvents is essential for the accurate characterization of biomolecules in molecular level research, which can reasonably explain the electrophoresis of charged colloids in high‐valent salt solutions presenting charge reversal and the ability of ions to promote the phase separation of block copolymers.In such a dielectric inhomogeneous system, a key challenge is how to solve electrostatic interactions quickly. For instance, the biomolecule‐solvent two medium model is a classical dielectric inhomogeneity system.[14] Obviously, since electrolyte solution systems contain biomolecules that usually dissociate a huge number of charged particles, calculating the electrostatic interactions of this system is an extremely time‐consuming task. A priori, for an N‐particle system, using traditional methods[15–18] (e.g., Ewald summation, particle–particle particle‐mesh, Coulomb force formula, etc.) to calculate these interactions would seem to require O(N2)$O({N}^{2})$ operations. We also need to recognize that in systems governed by the Poisson equation under dielectric inhomogeneity, the electrostatic forces on atom are not simply the product of electrostatic field E at the atom and the atomic charge q. However, these computational efforts may consume a lot of computer memory, especially when sufficient statistical sampling requires hundreds of thousands of simulations or a single simulation for extremely long time. Therefore, there is an urgent need for an efficient method to save computing time. In the past decades, these challenges have been addressed by classical numerical methods (e.g., finite difference method (FDM), finite element method (FEM)) as well as through popular linear or nonlinear iterative methods (e.g., successive over‐relaxation, conjugate gradient method, and multigrid method).[19–22] Among them, FDM is widely favored by researchers due to its theoretical simplicity and convenience in dealing with regular boundaries, while optimization efficiency and accuracy are two core problems for FDM research. A large number of excellent works have appeared one after another in this field.[23–28] In this study, we develop a hybrid approach to combine the advantages of Brownian dynamics and finite difference methods in the same simulation.[29,30] In this way, the calculation of electrostatic interactions is converted into the problem of solving the 3D Poisson equation under a dielectric inhomogeneous system in numerical techniques. The partial differential equations of this system are well defined at boundary and interface conditions. It can be solved approximately by the finite difference method as a sparse linear matrix problem. Consequently, it is feasible to adopt it for the electrostatic calculations of large‐scale biomolecular systems. In addition, another challenge for us is to deal with the molecular interface problem accurately and efficiently to achieve a smooth extension of the potential values over different dielectric regions. The main difficulties include three points: 1) the complex geometric configuration of the molecular surface; 2) the dielectric coefficient and electric field are discontinuous across the molecular interface; 3) the existence of singular sources. Recently, Egan[31] introduced a second‐order exact geometric discretization of point charge sources to solve the Poisson equation. In order to deal with the boundary conditions of the Poisson problem on irregular domain, Towers[32] developed a central difference discretization of the Dirac‐delta function representing the singular source at the interface. Moreover, some works have been dedicated to the study of the regularization in the setting of smoothly varying dielectric functions on diffuse interfaces and complex structures, of which four regularization methods based on Green's function are the most prevalent.The first approach decomposes the numerical solution into two fractions over the entire computational domain: the singular component (Green's function) and the residual component (the regular component). Despite the lack of accuracy, this scheme is one of the pioneering works in the field and many other subsequence works are motivated by it.[33,34] Chern et al.[35] proposed the second method: a regular format separation for the solution, which treats the regular and singular parts of the solution separately. This decomposition method has been applied to the subsequent finite element methods of Lu and Zhou et al.[36,37] The third method was proposed and explored in Xie's work.[38] This scheme divides the electrolyte solution into three parts over the entire region, which is applied to solvents that contain any number of ionic species. However, Xie's and Chern's schemes are both of finite element aspect. The fourth method, matched interface and boundary (MIB) technique based on finite difference, was proposed by Zhou and Wei et al.[39–45] The MIB method to smoothly extend the both sides of the interface by means of fictitious domains, then obtains a central finite difference scheme with high‐order accuracy repeatedly utilizing the interface conditions.However, in order to numerically solve the dynamic processes of electrostatic forces into the typical Brownian dynamics simulations of molecular conjugation or dissociation, which may involve millions or even tens of millions of time steps, the existing numerical methods are difficult to load such a large calculation amount. This means that it is necessary to find a new approach to avoid expensive computational costs and to ensure the accuracy of millions or even more calculations. In our previous work, we have proposed an optimization strategy that substantially improves the iteration speed and have used the method to calculate electrostatic interactions in the dielectric homogeneous system with excellent accuracy.[46]In this paper, the motion of charged polymers in a two‐dielectric model was simulated using finite difference scheme and simple yet effective interface treatment. We propose an iteration optimization strategy for large charged systems and further optimize the algorithm according to the shape of the polymer, which greatly improves the computational efficiency and saves a substantial CPU time. The rest of this paper is organized as follows. Section 2 presents the simulation model and the calculation scheme of interface problem, and proposes the corresponding algorithmic optimization strategy. In Section 3, we first verify the viability of the SOR global iteration method based on the FDM for the calculation of electrostatic forces. In addition, we are devoted to discussing the effects of the parameters in the “local iteration” on the computational speed and accuracy. The computational efficiency and adaptability of the banded local iteration method are further validated extensively. Finally, we draw some conclusions, and provide brief overview and some prospections in Section 4.Simulation Model and MethodsThe Biomolecule‐Solvent and Polymer ModelWe construct a coarse‐grained computational model for the movement of the charged polymers in a solvated biomolecular system, which can be described by the Poisson–Boltzmann equation (PBE) or viewed as a Poisson equation carrying a series of singular charge sources. As shown in Figure 1, the computational domain Ω is separated by the molecular surface Γ into the solute (biomolecule) region Ωm$\Omega_{\text{m}}$ and the solvent region Ωs$\Omega_{\text{s}}$, that is, Γ = Ωm∩Ωs$\Omega_{\text{m}}\cap \Omega_{\text{s}}$, Ω = Ωm⋃Ωs$\Omega_{\text{m}}\bigcup \Omega_{\text{s}}$. The boundary of the entire domain is ∂Ω$\partial \Omega$. In continuum models, biomolecules are usually approximated as homogeneous dielectrics, and the relative permittivity is taken to be between 1 and 20 (here set it as 1). In our simulations, in order to capture some important physical properties of biological systems more accurately and to further improve the quality of the classical Poisson model in electrostatic models, the dielectric value of the solvent region is replaced by a dielectric function dependent on the electric field strength.[47] The dielectric coefficient of the whole computational domain is a piecewise constant function given as1ε(r)=εm,r∈Ωmn2+απN0(n2+2)μν4πε0ELβ(n2+2)μνEkBT,r∈Ωs$$\begin{eqnarray} \varepsilon (r)=\left\{ \def\eqcellsep{&}\begin{array}{ll} {\varepsilon}_{\text{m}},& r\in \Omega_{\text{m}}\\ [6pt] {n}^{2}+\displaystyle\frac{\alpha \pi {N}_{0}({n}^{2}+2){\mu}_{\nu}}{4\pi {\varepsilon}_{0}E}L\left[\frac{\beta \left[({n}^{2}+2)\right]{\mu}_{\nu}E}{{k}_{\text{B}}T}\right],& r\in \Omega_{\text{s}} \end{array} \right.\nonumber\\ \end{eqnarray}$$where εm is the biomolecule region dielectric constant (the interface dielectric constants are assigned as the harmonic average of the solvent and solute dielectric constants), n is the optical refractive index, N0 is the number of molecules per unit volume, μν${\mu}_{\nu}$ is the dipole moment of the water molecule, ε0 is the vacuum dielectric constant, and E is the field strength. L(x)$L(x)$ denotes the usual Langevin function, α and β are numerical factors, T and kB represent the absolute temperature and Boltzmann constant, respectively. In our numerical implementation, the molecular parameters at T=300K$T=300\ \text{K}$ have values as following, n=1.33$n=1.33$, N0=3.334×1028m−3${N}_{0}=3.334\ensuremath{\times{}}{10}^{28}\ {\text{m}}^{-3}$, μν=2.1D${\mu}_{\nu}=2.1\ \text{D}$, α=28/373$\alpha =28/3\sqrt{73}$, β=73/6$\beta =\sqrt{73}/6$.1FigureSchematic illustration of the computational domain Ω: biomolecule‐solvent model. The molecular surface Γ divides the biomolecule region Ωm$\Omega_{\text{m}}$ and the solvent region Ωs$\Omega_{\text{s}}$.In the biomolecular region, the charge distribution is usually fixed because the charged particles are bound by chemical bonds and are not easily affected by the surrounding electrostatic field. For the solvent region, however, the charged particles (Na+${\text{Na}}^{+}$, Cl−${\text{Cl}}^{-}$, K+${\text{K}}^{+}$, Mg2+${\text{Mg}}^{2+}$, etc.) in the solution will change their spatial positions frequently due to the influence of the surrounding electrostatic field and the collision effect of the particles. Therefore, we use a probabilistic approach to quantitatively describe the distribution of charged particles in solvent. By Gauss's law, the system of equations for the electrostatic potential u(r)$u(r)$ of the model is expressed as2εm∇2u(r)=−ρm(r),r∈Ωm$$\begin{equation} {\varepsilon}_{\text{m}}{\nabla}^{2}u(r)=-{\rho}_{\text{m}}(r),\quad r\in \Omega_{\text{m}} \end{equation}$$3∇·(εs(r)∇u(r))=−ρs(r),r∈Ωs$$\begin{equation} \nabla \cdot ({\varepsilon}_{\text{s}}(r)\nabla u(r))=-{\rho}_{\text{s}}(r),\quad r\in \Omega_{\text{s}} \end{equation}$$4[u]Γ=0,r∈Γ$$\begin{equation} {[u]}_{\Gamma}=0,\quad r\in \Gamma \end{equation}$$5[εuξ]Γ=0,r∈Γ$$\begin{equation} {[\varepsilon {u}_{\xi}]}_{\Gamma}=0,\quad r\in \Gamma \end{equation}$$6u(r)=ub(r),r∈∂Ω$$\begin{equation} u(r)={u}_{b}(r),\quad r\in \partial \Omega \end{equation}$$where ρm(r)${\rho}_{\text{m}}(r)$ can be expressed by the sum of a series of Dirac‐Delta functions, which are redistributed to their neighboring grid points with appropriate interpolation formats. Thus the second‐order or higher‐order numerical accuracy at Cartesian grid points is realized. ρs(r)${\rho}_{\text{s}}(r)$ denotes the charge density distribution function in solvent region. [·]Γ${[\cdot ]}_{\Gamma}$ denotes the jump on Γ of the value of the function from Ωm$\Omega_{\text{m}}$ to Ωs$\Omega_{\text{s}}$. Equations (2) and (3) both represent dimensionless model equations. For example, the original equation of Equation (2) is εm∇2Φ=−ρm(r),r∈Ω${\varepsilon}_{\text{m}}{\nabla}^{2}\Phi=-{\rho}_{\text{m}}(r),r\in \Omega$. By the formula 1V=1J/C$1V=1J/C$, the unit of eckBT$\frac{{e}_{c}}{{k}_{\text{B}}T}$ is found to be 1/V$1/V$. Thus, by setting u=eckBTΦ$u=\frac{{e}_{\mathrm{c}}}{{k}_{\text{B}}T}\Phi$. Clearly, both ρ(r)$\rho (r)$ and u have unit m−2${\text{m}}^{-}2$ due to C2/(FJ)=1${C}^{2}/(\textit{FJ})=1$. Hence, the both sides of each equation of (2) have the same unit m−2${\text{m}}^{-2}$, which can be removed from the both sides of each equation to make (2) become dimensionless. The solution u of the Poisson model is called the dimensionless electrostatic potential. When u is found, the electrostatic potential Φ can be simply said to have a value of u in units kBT/ec${k}_{\text{B}}T/{e}_{\mathrm{c}}$. Equations (4) and (5) are the interface jump conditions across the molecular surface with ξ as the outer normal direction, where uξ=∂u∂ξ${u}_{\xi}=\frac{\partial u}{\partial \xi}$ is the directional derivative in ξ. Mathematically, the uniqueness of the solution is ensured by jump conditions. Physically, these conditions are necessary for the continuity of the electrostatic potential and its fluxes across the interface. Here, for grid‐based approach, the computational domain Ω is still assumed as a cube. As shown in Figure 2, we fix the electric potential on the left side of the simulation box as U and further get a Dirichlet boundary condition imposed on ∂Ω$\partial \Omega$.2Figure3D schematic of a charged polymer diffusing under dielectric inhomogeneity, where the dimensions of the simulation box are (L,L,L)$(L,L,L)$ and the left surface of the electric potential is fixed to U and the others to 0.We put a polymer into the domain Ω, which is modeled as a conventional bead‐spring chain consisting of N beads connected via the finitely extensible nonlinear elastic (FENE) potential[48]7UFENE(r)=−0.5kbb2ln(1−r2/b2)$$\begin{equation} {U}_{\text{FENE}}(r)=-0.5{k}^{b}{b}^{2}\text{ln}(1-{r}^{2}/{b}^{2}) \end{equation}$$the excluded volume interactions between the beads are considered through the Lennard–Jones (LJ) potential[49–51]8ULJ(r)=4ε[(σ/r)12−(σ/r)6]+ε,r≤21/6σ$$\begin{equation} {U}_{\textit{LJ}}(r)=4\varepsilon [{(\sigma /r)}^{12}-{(\sigma /r)}^{6}]+\varepsilon ,\quad r\le {2}^{1/6}\sigma \end{equation}$$where r and b are the bead–bead distance and the maximum bead–bead separation distance between consecutive beads along the polymer chain, kb${k}^{b}$ denotes the spring constant. ε represents the strength of the LJ potential, σ is used as the unit of length in this work and is the collision diameter of the LJ potential. Additionally, the charge of beads in the polymer is assumed to be q. The Poisson equation is solved with the interface and boundary potential conditions to compute the electrostatic interactions between charged beads. Hence, the Brownian dynamics can be applied to simulate the diffusion of charged polymers, including a noise term and a friction term. After specifying the interaction potentials, the BD equation of motion for particle i is9r⃗̇i(t)=−1μ▽Ui+1μf⃗i(t)$$\begin{equation} {\dot{\vec{r}}}_{i}(t)=-\frac{1}{\mu}\bigtriangledown {U}_{i}+\frac{1}{\mu}{\vec{f}}_{i}(t) \end{equation}$$where Ui${U}_{i}$ denotes the sum of all the interaction potentials discussed above, μ is the friction coefficient. f⃗i(t)${\vec{f}}_{i}(t)$ represents the random force associated with μ through the fluctuation dissipation theorem ⟨f⃗i(t)f⃗j(t′)⟩=2μkBTδijδ(t−t′)I$ \langle {\vec{f}}_{i}(t){\vec{f}}_{j}({t}^{{}^{\prime}}) \rangle =2\mu {k}_{\text{B}}T{\delta}_{\textit{ij}}\delta (t-{t}^{{}^{\prime}})\mathbf{I}$, where I is the unit tensor. The equation of motion of the charged polymer is integrated based on the time step Δt$\Delta t$. In our model, we present our results in reduced units, which result in a unit of the time as t=t*(ε/mσ2)1/2$t={t}^{\ast}{(\varepsilon /m{\sigma}^{2})}^{1/2}$, the charge as q=q*/(4πε0σε)1/2$q={q}^{\ast}/{(4\pi {\varepsilon}_{0}\sigma \varepsilon )}^{1/2}$, and the electric field as E=E*(4πε0σε)1/2σ/ε$E={E}^{\ast}{(4\pi {\varepsilon}_{0}\sigma \varepsilon )}^{1/2}\sigma /\varepsilon $. Where σ, ε and m denote the units of length, energy and mass, respectively. We assume that σ = ε = m = 1 and the temperature is kBT${k}_{\text{B}}T$ = 1.0, and the friction coefficient is set as μ=0.5Δt−1$\mu =0.5 \Delta {t}^{-1}$. Additionally, the dimensionless parameters of the polymer are fixed as kb=30${k}^{b}=30$ and b=1.5$b=1.5$, the charge of beads in the polymer is q=10$q=10$, and the time step is chosen as Δt=0.0002$\Delta t=0.0002$. We fix the size of the simulation box as L=80σ$L=80\sigma $, the electric potential of left surface as U=4$U=4$.Improved Interface Computation SchemeThis work mainly focuses on the fast solution of electrostatic interactions between charged particles in polymers with an applied electric field under the dielectric inhomogeneous system. While considering a uniform Cartesian mesh division of the domain Ω, the standard finite difference format loses its designed convergence around the interface. Hence the interface jump condition must be used to recover the accuracy. In order to achieve a smooth extension of the electric potential values, it is necessary to derive the computational expressions of the molecular interface, which is crucial for the situations where the interface intersects or does not intersect with the grid nodes. For this purpose, we divide all grid points into two categories: regular and irregular grid points. An irregular grid point is defined as a node in the standard finite difference format containing grid points that span the interface, that is, at least one of its six neighboring points is from different region. Otherwise, the points far from the interface Γ are regular points, performing the standard 7‐point central difference to discretize Equations (2) and (3). As shown in Figure 3, the blue and orange solid dots represent irregular points from the solute and solvent regions, respectively. The green and pink dots represent interface points that intersect and do not intersect with grid nodes, respectively. In our scheme, the difference format of irregular points is derived first, and then the computational format of interface points obtained by discretizing the interface conditions is substituted to modify the finite difference approximation of irregular points.3FigureThe interface value calculation scheme illustrated in a 2D setting with outside irregular points in blue and inside irregular points in orange. Green and pink dots indicate interface points that intersect and do not intersect with grid nodes, respectively. The light blue area denotes the solute area, the purple area denotes the irregular point distribution area, and the rest is the solvent area.For the electrostatic equations in 3D space with irregular dielectric interfaces, there are two interface conditions at any direction along the interface. In this way, the continuity of the potential flux Equation (5) is directed toward the normal direction ξ. At a given intersection point, it is convenient to introduce a coordinate transformation from (ξ, η, ζ) to (x, y, z), where η and ζ are the tangential and subnormal directions, respectively. Thus, the flux continuity condition Equations (4) and (5) are projected locally into the Cartesian grid10∂u∂η=(−uxssinθ+uyscosθ)−(−uxmsinθ+uymcosθ),$$\begin{equation} \left[\frac{\partial u}{\partial \eta}\right]=(-{u}_{x}^{\mathrm{s}}\sin\theta +{u}_{y}^{\mathrm{s}}\cos\theta )-(-{u}_{x}^{\mathrm{m}}\sin\theta +{u}_{y}^{\mathrm{m}}\cos\theta ), \end{equation}$$11∂u∂ζ=(−uxssinψcosθ−uyssinψsinθ+uzscosψ)−(−uxmsinψcosθ−uymsinψsinθ+uzmcosψ),$$\begin{eqnarray} \left[\frac{\partial u}{\partial \zeta}\right] &=& (-{u}_{x}^{\mathrm{s}}\sin\psi \cos\theta -{u}_{y}^{\mathrm{s}}\sin\psi \sin\theta +{u}_{z}^{\mathrm{s}}\cos\psi )\nonumber\\ &&-\ (-{u}_{x}^{\mathrm{m}}\sin\psi \cos\theta -{u}_{y}^{\mathrm{m}}\sin\psi \sin\theta +{u}_{z}^{\mathrm{m}}\cos\psi ), \end{eqnarray}$$12ε∂u∂ξ=εs(uxscosψcosθ+uyscosψsinθ+uzssinψ)−εm(uxmcosψcosθ+uymcosψsinθ+uzmsinψ),$$\begin{eqnarray} \left[\varepsilon \frac{\partial u}{\partial \xi}\right] &=& {\varepsilon}^{\mathrm{s}}({u}_{x}^{\mathrm{s}}\cos\psi \cos\theta +{u}_{y}^{\mathrm{s}}\cos\psi \sin\theta +{u}_{z}^{\mathrm{s}}\sin\psi )\nonumber\\ &&-\ {\varepsilon}^{\mathrm{m}}({u}_{x}^{\mathrm{m}}\cos\psi \cos\theta +{u}_{y}^{\mathrm{m}}\cos\psi \sin\theta +{u}_{z}^{\mathrm{m}}\sin\psi ), \end{eqnarray}$$where ψ denotes the angle between ξ and its projection on x–y plane while θ is the angle between this projection and x‐direction. uxs${u}_{x}^{\mathrm{s}}$ and uxm${u}_{x}^{\mathrm{m}}$ denote the first order derivatives in the x direction from different regions, respectively. In a uniform grid, the interface either intersects or does not intersect with the grid points. We first discuss the case where the interface intersects the grid points. At this time, all the irregular points use the successive over‐relaxation (SOR) iteration scheme with standard 7‐point difference13ui,j,k(n+1)=ui,j,k(n)+λ16ui+1,j,k(n)+ui−1,k(n+1)+ui,j+1,k(n)+ui,j−1,k(n+1)+ui,j,k+1(n)+ui,j,k−1(n+1)+h2ρi,j,kεi,j,k−ui,j,k(n)$$\begin{equation} {u}_{i,j,k}^{(n+1)}={u}_{i,j,k}^{(n)}+\lambda \left[\frac{1}{6}\left({u}_{i+1,j,k}^{(n)}+{u}_{i-1,k}^{(n+1)}+{u}_{i,j+1,k}^{(n)}+{u}_{i,j-1,k}^{(n+1)}+{u}_{i,j,k+1}^{(n)}+{u}_{i,j,k-1}^{(n+1)}+{h}^{2}\frac{{\rho}_{i,j,k}}{{\varepsilon}_{i,j,k}}\right)-{u}_{i,j,k}^{(n)}\right] \end{equation}$$where i, j and k denote the node position in the discretized grid, h is the grid spacing in the x, y, and z directions, the superscript n and n+1 represent the current iteration step and next iteration step, respectively. λ is the relaxation factor. It is possible to estimate numerically the optimal value of the parameter λ in the SOR method in some cases.[28] For instance, the optimal relaxation factor can be obtained according to the dimensions of the computational domain for both solvent and solute, that is, λoptimal=2/[1+sin(πh/Lmin)]${\lambda}_{\text{optimal}}=2/[1+\text{sin}(\pi h/{L}_{\text{min}})]$, where Lmin${L}_{\text{min}}$ is the minimum grid number of the simulation box in three directions. For the irregular points ui,j,k${u}_{i,j,k}$ in Equation (13), we assume that one of its six adjacent nodes, ui,j+1,k${u}_{i,j+1,k}$, is an interface point, then we need to combine Equations (5) and (12) to derive the interface point formula rather than use its value directly,14εscosψcosθui,j+1,k−ui−1,j+1,kh+cosψsinθui,j+1,k−ui,j,kh+sinψui,j+1,k−ui,j+1,k−1h=εmcosψcosθui+1,j+1,k−ui,j+1,kh+cosψsinθui,j+2,k−ui,j+1,kh+sinψui,j+1,k+1−ui,j+1,kh$$\begin{align} & {\varepsilon}^{\mathrm{s}} \left(\cos\psi \cos\theta \frac{{u}_{i,j+1,k}-{u}_{i-1,j+1,k}}{h}+\cos\psi \sin\theta \frac{{u}_{i,j+1,k}-{u}_{i,j,k}}{h}+\sin\psi \frac{{u}_{i,j+1,k}-{u}_{i,j+1,k-1}}{h}\right)\nonumber\\ & \quad ={\varepsilon}^{\mathrm{m}}\left(\cos\psi \cos\theta \frac{{u}_{i+1,j+1,k}-{u}_{i,j+1,k}}{h}+\cos\psi \sin\theta \frac{{u}_{i,j+2,k}-{u}_{i,j+1,k}}{h}+\sin\psi \frac{{u}_{i,j+1,k+1}-{u}_{i,j+1,k}}{h}\right) \end{align}$$according to the above equation, we will replace the interface point ui,j+1,k${u}_{i,j+1,k}$ in Equation (13) with the following formula15ui,j+1,k=(εi+1,j+1,kmui+1,j+1,k+εi−1,j+1,ksui−1,j+1,k)T1+(εi,j+2,kmui,j+2,k+εi,j,ksui,j,k)T2+(εi,j+1,k+1mui,j+1,k+1+εi,j+1,k−1sui,j+1,k−1)T3(εi+1,j+1,km+εi−1,j+1,ks)T1+(εi,j+2,km+εi,j,ks)T2+(εi,j+1,k+1m+εi,j+1,k−1s)T3$$\begin{equation} {u}_{i,j+1,k}=\frac{({\varepsilon}_{i+1,j+1,k}^{\text{m}}{u}_{i+1,j+1,k}+{\varepsilon}_{i-1,j+1,k}^{\text{s}}{u}_{i-1,j+1,k}){T}_{1}+({\varepsilon}_{i,j+2,k}^{\text{m}}{u}_{i,j+2,k}+{\varepsilon}_{i,j,k}^{s}{u}_{i,j,k}){T}_{2}+({\varepsilon}_{i,j+1,k+1}^{\text{m}}{u}_{i,j+1,k+1}+{\varepsilon}_{i,j+1,k-1}^{\text{s}}{u}_{i,j+1,k-1}){T}_{3}}{({\varepsilon}_{i+1,j+1,k}^{\text{m}}+{\varepsilon}_{i-1,j+1,k}^{s}){T}_{1}+({\varepsilon}_{i,j+2,k}^{\text{m}}+{\varepsilon}_{i,j,k}^{s}){T}_{2}+({\varepsilon}_{i,j+1,k+1}^{\text{m}}+{\varepsilon}_{i,j+1,k-1}^{\text{s}}){T}_{3}} \end{equation}$$with T1=cosψcosθ${T}_{1}=\cos\psi \cos\theta $, T2=cosψsinθ${T}_{2}=\cos\psi \sin\theta $, and T3=sinψ${T}_{3}=\sin\psi $.In addition, in the situation where the interface does not intersect with the grid points, the standard 7‐point difference scheme is no longer used for irregular points. We combine the mesh division and the geometry of biomolecule to identify the nodes involving interface in six directions, further deriving a 7‐point difference scheme containing the interface points. We use a second‐order accuracy approximation for ∂2u∂x2$\frac{{\partial}^{2}u}{\partial {x}^{2}}$, ∂2u∂y2$\frac{{\partial}^{2}u}{\partial {y}^{2}}$, and ∂2u∂z2$\frac{{\partial}^{2}u}{\partial {z}^{2}}$, which means that an additional point is needed in the direction involving the interface point to compensate for the loss of accuracy due to the interface. As shown in Figure 4, ui,j,k${u}_{i,j,k}$ has an interface point ub,j,k${u}_{b,j,k}$ among the difference points in the x‐direction, then an additional difference point ui+2,j,k${u}_{i+2,j,k}$ should be added. Using Taylor expansion and after algebraic operations, ∂2u∂x2$\frac{{\partial}^{2}u}{\partial {x}^{2}}$ can be discretized by ub,j,k${u}_{b,j,k}$, ui,j,k${u}_{i,j,k}$, ui+1,j,k${u}_{i+1,j,k}$, and ui+2,j,k${u}_{i+2,j,k}$, and its second‐order accuracy expression is16∂2u∂x2|ui,j,k=1hα−1α+2ui+2,j,k+2(2−α)α+1ui+1,j,k+α−3αui,j,k+6α(α+1)(α+2)ub,j,k$$\begin{eqnarray} \frac{{\partial}^{2}u}{\partial {x}^{2}}{|}_{{u}_{i,j,k}}&=&\frac{1}{h}\left[\frac{\alpha -1}{\alpha +2}{u}_{i+2,j,k}+\frac{2(2-\alpha )}{\alpha +1}{u}_{i+1,j,k}+\frac{\alpha -3}{\alpha}{u}_{i,j,k}\right.\nonumber\\ &&\left.+\ \frac{6}{\alpha (\alpha +1)(\alpha +2)}{u}_{b,j,k}\right] \end{eqnarray}$$4FigureLocal topology around irregular point (i,j,k)$(i,j,k)$. Take the x‐direction as an example, the interface crosses the x‐mesh line at the point B(b,j,k)$B(b,j,k)$ between (i,j,k)$(i,j,k)$ and (i−1,j,k)$(i-1,j,k)$, which are two irregular points. Two auxiliary lines (marked with dashed line) through B, one on the x–y plane and the other on the x–z plane, are drawn in order to derive the format for the calculation of the interface point B. Two auxiliary points (b,j,k+1)$(b,j,k+1)$ and (b,j,k+2)$(b,j,k+2)$ are placed on the auxiliary line on the x–z plane to facilitate the discretization of uzs${u}_{z}^{\text{s}}$; Additionally, two auxiliary points (b,j−1,k)$(b,j-1,k)$ and (b,j−2,k)$(b,j-2,k)$ are placed on the auxiliary line on the x–y plane to facilitate the discretization of uys${u}_{y}^{\text{s}}$. αh$\alpha h$ and βh$\beta h$ represent the interface point distances of (i,j,k)$(i,j,k)$ to the x and y directions, respectively.In order to obtain a second‐order accuracy approximation for ∂2u∂y2$\frac{{\partial}^{2}u}{\partial {y}^{2}}$ and ∂2u∂z2$\frac{{\partial}^{2}u}{\partial {z}^{2}}$, a similar process must be performed. In Figure 4, taking the x and y directions involving the difference points as an example. After discretizing the Poisson equation, we will obtain the SOR iteration format of ui,j,k${u}_{i,j,k}$17ui,j,k(n+1)=ui,j,k(n)+λαββ(α−3)+α(β−3)−2αβ6α(α+1)(α+2)ub,j,k(n+1)+2(2−α)1+αui+1,j,k(n)+α−12+αui+2,j,k(n)+6β(β+1)(β+2)ui,f,k(n+1)+2(2−β)1+βui,j+1,k(n)+β−12+βui,j+2,k(n)+ui,j,k+1(n)+ui,j,k−1(n+1)+h2ρi,j,kεi,j,k−ui,j,k(n)$$\begin{eqnarray} {u}_{i,j,k}^{(n+1)}& =&{u}_{i,j,k}^{(n)}+\lambda \left\{\frac{\alpha \beta}{\beta (\alpha -3)+\alpha (\beta -3)-2\alpha \beta}\right.\nonumber\\ && \left.\left[\frac{6}{\alpha (\alpha +1)(\alpha +2)}{u}_{b,j,k}^{(n+1)}+\frac{2(2-\alpha )}{1+\alpha}{u}_{i+1,j,k}^{(n)}+\frac{\alpha -1}{2+\alpha}{u}_{i+2,j,k}^{(n)}\right.\right.\nonumber\\ && \left.\left.+\ \frac{6}{\beta (\beta +1)(\beta +2)}{u}_{i,f,k}^{(n+1)}+\frac{2(2-\beta )}{1+\beta}{u}_{i,j+1,k}^{(n)}+\frac{\beta -1}{2+\beta}{u}_{i,j+2,k}^{(n)}\right.\right.\nonumber\\ && \left.\left.+\ {u}_{i,j,k+1}^{(n)}+ {u}_{i,j,k-1}^{(n+1)}+{h}^{2}\frac{{\rho}_{i,j,k}}{{\varepsilon}_{i,j,k}}\right]-{u}_{i,j,k}^{(n)}\right\} \end{eqnarray}$$where α and β denote the distance ratios of ui,j,k${u}_{i,j,k}$ to the interface points in the x and y directions, respectively. Furthermore, the modified finite difference approximation at irregular points maintains second‐order accuracy, provided that the interface point expressions are accurately estimated. Therefore, the interface points ub,j,k${u}_{b,j,k}$ and ui,f,k${u}_{i,f,k}$ have their special expressions in the above computational format. Inspired by Zhou,[41] we use a linear combination of the interface jump conditions to eliminate any two first‐order derivatives that are not easily approximated, which two are eliminated depends on the molecular geometry. A specific example is discussed to further illustrate the expression for the interface point B according to Figure 4. We will keep uxm${u}_{x}^{\text{m}}$ and uxs${u}_{x}^{\text{s}}$ to generate the jump condition along the x‐direction. In the y‐direction, it is more difficult to calculate uym${u}_{y}^{\text{m}}$ at (b,j,k)$(b,j,k)$ than uys${u}_{y}^{\text{s}}$. In the z‐direction, computing uzm${u}_{z}^{\text{m}}$ at (b,j,k)$(b,j,k)$ is more difficult than computing uzs${u}_{z}^{\text{s}}$. Therefore, uym${u}_{y}^{\text{m}}$ and uzm${u}_{z}^{\text{m}}$ will be eliminated by Equations (10)–(12) in interface format18ε∂u∂ξ−εmtanψ∂u∂ζ+Cη∂u∂η=Cxsuxs+Cxmuxm+Cysuys+Czsuzs$$\begin{eqnarray} &&\left[\varepsilon \frac{\partial u}{\partial \xi}\right]-{\varepsilon}^{\text{m}}\tan\psi \left[\frac{\partial u}{\partial \zeta}\right]+{C}_{\eta}\left[\frac{\partial u}{\partial \eta}\right]={C}_{x}^{\text{s}}{u}_{x}^{\text{s}}+{C}_{x}^{\text{m}}{u}_{x}^{\text{m}}\nonumber\\ &&\quad +\ {C}_{y}^{\text{s}}{u}_{y}^{\text{s}}+{C}_{z}^{\text{s}}{u}_{z}^{\text{s}} \end{eqnarray}$$where the superscripts m$\text{m}$ and s$\text{s}$ represent the molecular region and the solvent region, respectively. ux${u}_{x}$, uy${u}_{y}$, and uz${u}_{z}$ are the derivative of the electric potential, which are defined in the normal direction ξ at the dielectric interface. Then, through approximating the remaining four first‐order derivatives, we obtain the expression for ub,j,k${u}_{b,j,k}$19ub,j,k=(1−α)Cxsui−1,j,k−αCxmui,j,k+(α+α2)Cysub,j−1,k−(α+α2)Czsub,j,k+1(1−α)Cxs−αCxm+(α+α2)Cys−(α+α2)Czs$$\begin{eqnarray} {{ \def\eqcellsep{&}\begin{array}{l} {u}_{b,j,k}=\displaystyle\frac{(1-\alpha ){C}_{x}^{\text{s}}{u}_{i-1,j,k}-\alpha {C}_{x}^{\text{m}}{u}_{i,j,k}+(\alpha +{\alpha}^{2}){C}_{y}^{\text{s}}{u}_{b,j-1,k}-(\alpha +{\alpha}^{2}){C}_{z}^{\text{s}}{u}_{b,j,k+1}}{(1-\alpha ){C}_{x}^{\text{s}}-\alpha {C}_{x}^{\text{m}}+(\alpha +{\alpha}^{2}){C}_{y}^{\text{s}}-(\alpha +{\alpha}^{2}){C}_{z}^{\text{s}}} \end{array} }}\nonumber\hspace*{-10pt}\\ \end{eqnarray}$$where Cxs=εscosψcosζ+εmtanψsinψcosζ+εmsinζtanζ/cosψ${C}_{x}^{\text{s}}={\varepsilon}^{\text{s}}\cos\psi \cos\zeta +{\varepsilon}^{\text{m}}\tan\psi \sin\psi \cos\zeta +{\varepsilon}^{\text{m}}\sin\zeta \tan \zeta /\cos\psi $, Cxm=−εm/(cosψcosζ)${C}_{x}^{\text{m}}=-{\varepsilon}^{\text{m}}/(\cos\psi \cos\zeta )$, Cys=(εs+εm)cosψsinζ${C}_{y}^{\text{s}}=({\varepsilon}^{\text{s}}+{\varepsilon}^{\text{m}})\cos\psi \sin\zeta $, Czs=(εs−εm)sinψ${C}_{z}^{\text{s}}=({\varepsilon}^{\text{s}}-{\varepsilon}^{\text{m}})\sin\psi $. Note that the auxiliary points ub,j−1,k${u}_{b,j-1,k}$ and ub,j,k+1${u}_{b,j,k+1}$ are not grid nodes in the above equation, which can be calculated by interpolating points on the same horizontal line in the same region. For example, we interpolate ub,j−1,k${u}_{b,j-1,k}$ with ui−2,j−1,k${u}_{i-2,j-1,k}$, ui−1,j−1,k${u}_{i-1,j-1,k}$ and ui,j−1,k${u}_{i,j-1,k}$, and then substitute these interpolations back into the above interface point expression. In addition, these solutions can be extended to more complicated macromolecular (solute) shapes.For the discretized Poisson equation, a successive over‐relaxation iteration (SOR) stabilized solver is adopted to solve the linear system. The maximum numerical error is measured, in terms of the discrete L∞${L}_{\infty}$ norm20eu=maxi=1,2,⋯,N|uexa(ri)−unum(ri)|$$\begin{equation} {e}_{u}=\underset{i=1,2,\text{\ensuremath{\cdots}},N}{\max}|{u}_{\text{exa}}({r}_{i})-{u}_{\text{num}}({r}_{i})| \end{equation}$$where ri${r}_{i}$ is the regular grid points in solvent and solute region (use standard 7‐point differential scheme) or the irregular grid points near the interface (use new interface scheme), and the total number is denoted as N. Here, uexa(ri)${u}_{\text{exa}}({r}_{i})$ and unum(ri)${u}_{\text{num}}({r}_{i})$ are the exact solution and the numerical solution, respectively. We are interested in examining the performance of the proposed method on the interface by comparing it with the standard 7‐point differential scheme.We conducted a number of numerical experiments to examine the performance of the proposed new computational scheme near the interface, with particular attention to the order of convergence. As shown in Table 1, the grid spacing h=0.01$h=0.01$ is set as the exact solution (there is no exact solution since the right‐hand side function of the model equation is obtained from the random motion of the polymer and has no fixed form). Through solving Equations (2) and (3), we report the electrostatic potential error eu${e}_{u}$ and its convergence order regarding the grid refinement. As can be seen from Table 1, the new computational format achieves very similar results to the traditional format in terms of electric potential. The numerical results indicate that the convergence in electric potential error eu${e}_{u}$ is second order. As far as the potential value is concerned, the results obtained at coarse grid (e.g., h=1.0$h=1.0$) are very close to that at the exact solution (i.e., fine grid h=0.01$h=0.01$).1TableSolving the Poisson equation: error and order of convergence for the standard 7‐point difference scheme and the new interface schemeStandard 7‐point schemeImproved interface schemeheu${e}_{u}$Ordereu${e}_{u}$Order19.8e−047.3e−041/25.7e−042.072.6e−042.221/41.2e−041.988.9e−051.961/83.8e−052.279.6e−062.451/167.7e−062.241.1e−062.31Proposed and Optimized Local Iteration StrategyAbove interface treatment scheme lays the basic framework for solving elliptic interface problems. Solving the differential discretized system of equations using SOR iteration indeed has certain advantages. While in Brownian dynamics, one usually need to compute millions or even tens of millions of time steps, which means that the calculation of electrostatic interactions will take a lot of time. Actually, the charged particles in polymer tend to occupy only a small region in simulation box, and the spatial distribution of electrostatic potential changes little due to the small perturbation in dielectric and charge distributions at each step during such simulations. Especially when an electric field was applied, we find that the value of electric potential hardly changes within a certain time interval for grid nodes that are far from the charged polymer. This result suggests that in a short simulation time, one can save most of the iteration time if we only calculate the electric potential of grids near the charged polymer. The electric potential for other grids is fixed and updated after several steps. Overall, we can save substantial insignificant iteration time by using local iteration instead of global iteration while without much compromise in accuracy. Three important parameters need to be introduced here to ensure that the local iteration method can improve the computational efficiency while also guaranteeing a certain accuracy. As shown in Figure 5, we first locate the spatial location of the beads in polymer and enclose them with a minimal cubic basal box (A1B1C1D1${A}_{1}{B}_{1}{C}_{1}{D}_{1}$). Then, from six sides of the basal box, we extend outward the grid of Lw${L}_{\text{w}}$ units to form a “local box.” Once the “local box” exceeds the boundary of the domain ∂Ω$\partial \Omega$ in a certain direction, it will not further extend. Additionally, the parameter dstep${d}_{\text{step}}$ is defined as the step interval between the first global iteration and the next global iteration. We also set a warning value (Ld), which is considered as the threshold for updating the local iteration area. When the nearest distance between the basal box and the “local box” is equal to the warning value, the “local box” will be re‐divided based on current position of basal box.5Figure2D view of the basal box (A1B1C1D1${A}_{1}{B}_{1}{C}_{1}{D}_{1}$) and “local box” (A2B2C2D2${A}_{2}{B}_{2}{C}_{2}{D}_{2}$) under dielectric inhomogeneity. Lw denotes the length of the base box extending outward, and Ld is the warning value between the basal box and the local box. Further, A1B1C1D1${A}_{1}{B}_{1}{C}_{1}{D}_{1}$ is defined as the local basal box, and the green highlighted area indicates the banded local basal box. After several steps, the polymer is moved to virtual box (A1′B1′C1′D1′${A}_{1}^{{}^{\prime}}{B}_{1}^{{}^{\prime}}{C}_{1}^{{}^{\prime}}{D}_{1}^{{}^{\prime}}$).In fact, we can further optimize the local iteration method for polymers with special shapes. If the polymer is clustered on a certain plane (e.g., x–y, x–z, y–z), then local iteration can maximize its advantages. However, when the polymer is diagonal or dispersed in 3D space, the efficiency of the speedup is dramatically reduced. To address the adverse effects of these shape characteristics, we propose an adaptive banded local iteration strategy and expect that this approach will save a significant amount of computational expense and further improve computational efficiency. As shown in Figure 5, the banded iteration area is marked in green. The previous local iteration method is called traditional local iteration (TLI) method,[46] and the new adaptive local iteration strategy is called banded local iteration (BLI) method. Importantly, it is quite common for polymers to overlap domains in molecular simulations. In response, our algorithm is also able to accurately locate the banded local iterative regions.In our model, a new interface format is derived using a Cartesian grid, a standard finite difference scheme, and the lowest order interface jump condition. Then, the format is combined with a local iteration method to solve the elliptic equations with discontinuity coefficients and relatively smooth interfaces to simulate the movement of the polymer. We are more concerned with the performance of the local iteration algorithm under dielectric inhomogeneity rather than with how to handle interface problems with singularities. Since there has been a lot of work on the explicit interface treatment of geometric singularities, which is not the focus of this paper. Therefore, we only calculate two molecular examples with different geometrical features and classify them into two categories: molecular interfaces and grid nodes all intersect (called prism example) and partially intersect (the solute shape similar to peanut, called ellipsoid example). In different molecular configurations, we simulated three cases (as shown in Figure 6) in which the same shape of polymer was put at three different initial positions: the solvent region, the intra‐molecular region, and the junction of solvent and biomolecule, respectively. Finally, we enumerate two polymers with special shapes, and compare the numerical results of traditional local iteration and banded local iteration in two different molecular configuration computational domains.6FigureThe polymer is placed in three different positions in the computational domain, cases 1–3 corresponding to the solvent region, the internal region of the biomolecule (i.e., solute region) and the junction between solvent and molecule, respectively.Results and DiscussionIn the numerical experiments, all codes were written in Matlab 2017b and run on PC with Intel Core i5‐7500 CPU (3.40 GHz) and 8 GB RAM. At first, independent effects of three parameters (Lw${L}_{\text{w}}$, dstep${d}_{\text{step}}$, and Ld${L}_{\text{d}}$) on the computational speed and accuracy are briefly discussed. Additionally, we mainly focus on three aspects of simulations: the polymer motion at different positions in computational domain, the intersection of molecular interface with grid points (i.e., molecular shape), and the shape of the polymers. In order to make the results clearer and more universal, we will use the relative values, that is, RLw=2×Lw/L${\text{RL}}_{\text{w}}=2\ensuremath{\times{}}{L}_{\text{w}}/L$, RS=dstep/Step$\text{RS}={d}_{\text{step}}/\text{Step}$, RLd=Ld/Lw${\text{RL}}_{\text{d}}={L}_{\text{d}}/{L}_{\text{w}}$. For example, a local box is formed by extending Lw units in six directions from the center of the basal box, thus RLw${\text{RL}}_{\text{w}}$ portrays the relative size of the local box well. When RLw=100%${\text{RL}}_{\text{w}}=100\%$, whether the “local box” completely covers the simulation box depends on the initial position of the charged beads in a polymer.Benchmarked with SOR global iteration solution, both terms of accuracy and computational efficiency are investigated for the local iteration method through 3D simulation. In terms of accuracy, we report the electric potential error and the error in the ultimate position of charged particles because the final objective is to simulate the motion of charged particles, which needs an accurate calculation of the spatial potential distribution and ultimate position of the polymer. Therefore, the relative error of potential and relative position error of charged particles are defined as follows21REU=ΘRMS(Uglobal−Ulocal)/Uglobal$$\begin{equation} {\text{RE}}_{\text{U}}={\Theta}_{\text{RMS}}\left[({U}_{\text{global}}-{U}_{\text{local}})/{U}_{\text{global}}\right] \end{equation}$$22REPos=ΘRMS|Posglobal−Poslocal|/|Posglobal−Posinitial|$$\begin{equation} {\text{RE}}_{\text{Pos}}={\Theta}_{\text{RMS}}\left[|{\text{Pos}}_{\text{global}}-{\text{Pos}}_{\textit{local}}|/|{\text{Pos}}_{\text{global}}-{\text{Pos}}_{\text{initial}}|\right] \end{equation}$$where ΘRMS${\Theta}_{\text{RMS}}$ denotes the root mean square sign, the overall similarity of the polymer conformation is checked by calculating the RMS (root mean square) deviations of all charged particle positions and spatial electric potential. Uglobal${U}_{\text{global}}$ and Ulocal${U}_{\text{local}}$ represent the potential distribution at each charged particles in the global iteration and the local iteration, respectively. Posinitial${\text{Pos}}_{\text{initial}}$, Posglobal${\text{Pos}}_{\text{global}}$, and Poslocal${\text{Pos}}_{\text{local}}$ denote the 3D spatial coordinates of charged particles in the initial, global iteration and local iteration, respectively, and the displacement distance of each particle is expressed by the absolute value of the difference between them. In terms of computational speed, the CPU times of global iteration (denoted as tG) and local iteration (denoted as tL) were recorded and compared. In previous work,[46] we chose a simulation time of 10 000 time steps, which was long enough to test the stability of the simulation results. However, the global iteration simulation under dielectric inhomogeneity takes a lot of CPU time, so that we set the total step number of motion Step=1000$\text{Step}=1000$ in all simulations.In order to verify the validity of the SOR global iteration method based on a finite difference scheme (with a grid spacing h=1$h=1$), we compare it with existing conventional algorithms for calculating electrostatic forces. In particular, the Coulomb formulation, an exact method for solving electrostatic interactions, will be used to evaluate the validity of our proposed numerical method. During the simulation, we use the global iteration method and Coulomb's formula in the section of solving for the electrostatic force respectively, while keeping the settings consistent in all other sections, and then record the CPU time and numerical results such as particle position and spatial potential. We display the CPU time as a function of the number of charged particles (Num$\text{Num}$) for these two methods in Figure 7. The comparison results indicates that as the number of charged particles increases, the CPU time of these two methods both exhibits an increase. However, the SOR method is different from the Coulomb formula for solving the electrostatic force directly, which is insensitive to the increases in Num$\text{Num}$. Obviously, as the number of particles increases, the SOR global iteration only increases the workload of interpolating more charges to nearby nodes, which can even be ignored compared with time‐consuming iteration loop. On the contrary, with increasing number of particles, the CPU time of the Coulomb formula increases significantly. The numerical results indicate that the SOR method is more suitable for large‐scale charged particle system. This conclusion is supported by the results in the inset of Figure 7, where the computational efficiency of the global iteration gradually outperforms the Coulomb formula when Num>400$\text{Num}&gt;400$. As Num$\text{Num}$ continues to increase, the computational efficiency advantage of the global iteration becomes more pronounced, and it is nine times more efficient than the Coulomb formula when Num=1500$\text{Num}=1500$. In general, the calculation complexity of the Coulomb formula is O(N2)$O({N}^{2})$, while the calculation complexity of a single iteration of the SOR global iteration method with optimal relaxation parameters is between O(NlnN)$O(N\text{ln}N)$ and O(NlgN)$O(N\text{lg}N)$. Theoretically, there is a significant difference in computational effort between the two methods, and Figure 7 should show a corresponding CPU time difference. In fact, the iteration method requires a loop of iterations until the convergence criterion is satisfied. We assume that a complete loop process requires the implementation of n iterations, then the total computational effort to compute the single electrostatic force needs to be multiplied by a factor n. In addition, using the results of the Coulomb formula as a benchmark, Table 2 shows the root mean square deviation of the charged particle position and electric potential for global iteration method. We find that as NUM increases, both REU${\text{RE}}_{\mathrm{U}}$ and REPos${\text{RE}}_{\text{Pos}}$ increase. In the REU${\text{RE}}_{\mathrm{U}}$ term, as Num$\text{Num}$ increases from 100 to 1500, the error increases by only about 1.12e−3. In the REPos${\text{RE}}_{\text{Pos}}$ term, although the error increases from 6.66e−4 to 1.05e−3, we can observe that the increase in error is not significant when Num>1000$\text{Num}&gt;1000$. Overall, the two errors are controlled to be below 3.16e−3 and 1.05e−3, respectively, both in the 1e−3 magnitude, which is an acceptable error range for the experiment. Therefore, it is a feasible scheme to apply a global SOR iteration method based on a finite difference format to solve the electrostatic force and thus simulate charged polymer motion.2TableThe root mean square deviation of the charged particle position and electric potential10030050010001500REU${\text{RE}}_{\text{U}}$2.04e−32.42e−32.67e−33.02e−33.16e−3REPos${\text{RE}}_{\text{Pos}}$6.66e−48.20e−49.01e−49.92e−41.05e−37FigureCPU time t as a function of the number of charged particles (Num$\text{Num}$). The inset shows the relative percentage of CPU time for the global iteration method and the Coulomb formula at different numbers of particles.In Figures 8 and 9, parts (a) and (b), we display CPU time t and the relative error as a function of RLw${\text{RL}}_{\mathrm{w}}$ (the length ratio of “local box” to “global box”), respectively, where the relative error is divided into two parts with the relative error of electrostatic potential (REU${\text{RE}}_{\mathrm{U}}$) and the relative error of the final position of the beads (REPos${\text{RE}}_{\text{Pos}}$). Meanwhile, the dimension of simulation box was set as L=80$L=80$ with grid spacing h=Δx=Δy=Δz=1$h= \Delta x= \Delta y= \Delta z=1$, which means that the total number of grid nodes is 81×81×81$81\ensuremath{\times{}}81\ensuremath{\times{}}81$. The number of charged particles in polymer is set as Num=500$\text{Num}=500$ (the FDM method has significant advantages for handling electrostatic interactions containing large numbers of charged particles), the relative distance between charged beads and boundary RLd=80%${\text{RL}}_{\mathrm{d}}=80\%$. And the value of relative global iteration step interval RS$\text{RS}$ is set as 10% and 25% in simulations to explore the impact of update frequency of global iteration on the computational speed and accuracy. After verification, we find that the increasing RLd${\text{RL}}_{\mathrm{d}}$ is beneficial to reduce REU${\text{RE}}_{\mathrm{U}}$, while the CPU time only slightly increases at the same time. Therefore, we uniformly set RLd${\text{RL}}_{\mathrm{d}}$ in all the examples to 80%. Under these conditions, we calculate tG, tL, REU${\text{RE}}_{\mathrm{U}}$, and REPos${\text{RE}}_{\text{Pos}}$ of the two numerical examples of prismatic and ellipsoidal biomolecules, in which the numerical results of different examples are represented by different colored lines. Based on these two examples, we constructed three cases as shown in Figures 8 and 9.8FigureThe case of polymer in solvent region: a) CPU time t and b) the relative error of electrostatic potential REU${\text{RE}}_{\text{U}}$ as a function of the length ratio between “local box” and “global box” (RLw${\text{RL}}_{\text{w}}$). The relative distance between charged polymer and boundary is set as RLd${\text{RL}}_{\text{d}}$ = 80%. a) RS$\text{a) RS}$ represents relative global iteration step interval, and tG represents the global iteration time. b) The relative error of the final position of the charged particles (REPos${\text{RE}}_{\text{Pos}}$) as a function of the relative size of local box (RLw${\text{RL}}_{\mathrm{w}}$).9FigureThe case of polymer in internal region of the solute: a) CPU time t and b) the relative error of electrostatic potential REU${\text{RE}}_{\text{U}}$ as a function of the length ratio between “local box” and “global box” (RLw${\text{RL}}_{\text{w}}$). The relative distance between charged polymer and boundary is set as RLd${\text{RL}}_{\text{d}}$ = 80%. RS$\text{RS}$ represents relative global iteration step interval, and tG in (a) represents the global iteration time. The inset of (b) shows the relative error of the final position of the charged particles (REPos${\text{RE}}_{\text{Pos}}$) as a function of the relative size of local box (RLw${\text{RL}}_{\mathrm{w}}$).Initially, as shown in Figure 6, we put the polymer chain in solvent area and the charged particles of the polymer chain are under electrostatic interaction, elastic force and repulsive force to obtain an equilibrium conformation (case 1). For the solvent region, real‐time updating of the dielectric value everywhere in the solvent area may spend more computing costs because the dielectric value is not a constant but as a function of the electric field strength, which implies that the calculation of electrostatic interaction will consume more CPU time. However, these calculation costs are necessary since the dielectric response can accurately capture certain characteristics of the biological process. As shown in Figure 8a, with the decrease of RLw${\text{RL}}_{\mathrm{w}}$, the CPU time decreases. Obviously, reducing the local iteration domain means saving a part of the iteration time. In addition, the computational expense of iteration increases with the decrease of RS$\text{RS}$. In other words, as the RS$\textit{RS}$ decreases, the frequency of global iteration increases and the cumulative error decreases, which further leads to the increase of total CPU time and the decrease of REU${\text{RE}}_{\mathrm{U}}$. These results imply that a smaller RLw${\text{RL}}_{\mathrm{w}}$ and a larger RS$\text{RS}$ are beneficial to improve calculation efficiency. Under the same conditions, the ellipsoid group reported significantly more CPU time than the prismatic group. There are two main reasons: in the program of ellipsoid example, it is necessary to open up an independent space to store the information of interface values that do not intersect with nodes. Otherwise, the calculation format of these interface points is more complicated. As the local box continually expands, we find that the tL of both prismatic and ellipsoid examples are approximately equal to tG when RLw=90%${\text{RL}}_{\text{w}}=90\%$, because the smallest basal box occupies part of the computational domain. It's worth noting that during RLw${\text{RL}}_{\mathrm{w}}$ increasing from 35% to 60%, a significant increase in the CPU time was recorded. According to the preliminary analysis, when RLw<35%${\text{RL}}_{\text{w}}&lt;35\%$, the local iteration is only performed in the solvent area. With the successive increase of RLw${\text{RL}}_{\mathrm{w}}$, the local iteration range further covers the molecular (solute) area. As a result, the computational complexity of nodes near the molecular interface significantly increases the CPU time. When RLw>60%${\text{RL}}_{\text{w}}&gt;60\%$, the local iteration area completely covers the entire biomolecule, and this trend stabilizes. Besides, the results of REPos${\text{RE}}_{\text{Pos}}$ illustrate the relative error of the polymer translocation, where the REPos${\text{RE}}_{\text{Pos}}$ and REU${\text{RE}}_{\mathrm{U}}$ share a similar change trend. As shown in Figure 8b, REU${\text{RE}}_{\mathrm{U}}$ and REPos${\text{RE}}_{\text{Pos}}$ decrease with the increase of RLw${\text{RL}}_{\mathrm{w}}$. These results indicate that choosing a larger RLw${\text{RL}}_{\mathrm{w}}$ gets closer to the global solution. We find that the REU${\text{RE}}_{\mathrm{U}}$ has dropped by an order of magnitude when RLw≈18%${\text{RL}}_{\text{w}}\approx 18\%$. Moreover, we also find that there is a difference in REU${\text{RE}}_{\mathrm{U}}$ between the two examples. The reason for this phenomenon is that the ellipsoid calculation example compensate the loss of accuracy by repeatedly using the interface conditions and adding iteration nodes near the interface. Importantly, when RLw=7.5%${\text{RL}}_{\text{w}}=7.5\%$ and RS=25%$\text{RS}=25\%$, the single most striking observation to emerge from the data comparison was the calculation speed of the two examples is increased by 20 and 17 times respectively, and the REU${\text{RE}}_{\mathrm{U}}$ is controlled below 3.2e−4.Besides, we also put the polymer chain in solute area and observed its motion and the performance of local iteration method (case 2). The results in Figure 9a show that CPU time decreases with the decrease of RLw${\text{RL}}_{\mathrm{w}}$. At the same time, a decreasing RS$\text{RS}$ also induces the decrease of CPU time under the same RLw${\text{RL}}_{\mathrm{w}}$. As shown in Figure 9b, with the increase of RLw${\text{RL}}_{\mathrm{w}}$, the REU${\text{RE}}_{\mathrm{U}}$ first decreases rapidly, then remains unchanged. In the inset of Figure 9b shows the relative position error (REPos${\text{RE}}_{\text{Pos}}$) as a function of RLw${\text{RL}}_{\mathrm{w}}$, which decreases as RLw${\text{RL}}_{\mathrm{w}}$ increases. In fact, if RLw${\text{RL}}_{\mathrm{w}}$ is large enough the error will be almost negligible. These situations are very similar to case 1, which is exactly what we expected. Consequently, we can get a conclusion that the local iteration algorithm is not affected by the spatial position of the polymer, which reflects its stability. Interestingly, we find that when RLw<50%${\text{RL}}_{\text{w}}&lt;50\%$, cases 1 and 2 presented large differences in CPU time. Theoretically, as long as the local box does not contact with the molecular interface, the node values in the local iterative domain are calculated by the standard 7‐point center difference scheme. However, because the molecular volume is only 1.5% of the computational domain and the large number of charged beads in polymer located inside the biomolecule has occupied almost the entire biomolecule (Num=500$\text{Num}=500$), it drives the local box to contact the molecular interface prematurely (in case 2). Therefore, the computational amount of case 2 is much higher than that of case 1 even in the initial stage of RLw${\text{RL}}_{\mathrm{w}}$ increase. In a way, the acceleration effect of local iteration in case 2 is not as good as that in cases 1 and 3.In fact, as shown in Figure 6, setting the polymer chain at the junction of solute and solvent area is another important situation for exploring the molecular behavior on singlemolecule level, which takes charged polymer intersected vesicles as prototype (case 3). Fast calculation of electrostatic potential provides convenience for applications of vesicles in pharmaceutical and physiological processes. As shown in Figure 10a, with the increase of RLw${\text{RL}}_{\mathrm{w}}$, the CPU time first increases rapidly, and then the rate of increase becomes slower. Similarly, Figure 10b shows that REU${\text{RE}}_{\mathrm{U}}$ decreases with the decrease in RS$\text{RS}$ for the same RLw${\text{RL}}_{\mathrm{w}}$, which indicates that increasing the frequency of global iteration can significantly reduce the accumulation of errors. When RLw=37%${\text{RL}}_{\text{w}}=37\%$, the CPU time of the ellipsoid example gradually draws a significant gap with the prism example. It is also at this stage that the growth in CPU time becomes stable. Overall, the result of case 3 is very similar to the other two cases, which is consistent with our expectations. While, the CPU time of case 3 first increases rapidly and then increases slowly, opposite to cases 1. Through analysis, we find that this has a great relationship with the region contained in local iteration box. When the polymer spans two regions, no matter how small the local box, the initial calculation will always involve interface nodes, leading to more difficult iterative solution. Additionally, in the process of polymer movement, free moving ions are distributed on the left, while fixed charges are distributed in the solute region on the right, which further increases the computational cost. Even though the calculation speed decreases slightly, when RLw=7.5%${\text{RL}}_{\text{w}}=7.5\%$ and RS=25%$\text{RS}=25\%$, the computation efficiency still increased by 8.2 times compared to the global iteration in the ellipsoid numerical example.10FigureThe case where the polymer is located at the solute‐solvent junction: a) CPU time t and b) the relative error of electrostatic potential REU${\text{RE}}_{\text{U}}$ as a function of the length ratio between “local box” and “global box” (RLw${\text{RL}}_{\text{w}}$). The relative distance between charged polymer and boundary is set as RLd${\text{RL}}_{\text{d}}$ = 80%. RS$\text{RS}$ represents relative global iteration step interval, and tG in (a) represents the global iteration time. The inset of (b) shows the relative error of the final position of the charged particles (REPos${\text{RE}}_{\text{Pos}}$) as a function of the relative size of local box (RLw${\text{RL}}_{\mathrm{w}}$).In order to further investigate our analysis in detail, we combine the results of the three position cases of the two calculation examples and shown them in Figures 11a and 11b, respectively. The left axis and right axis denote CPU time t and the relative potential error (REU${\text{RE}}_{\mathrm{U}}$) as a function of RLw${\text{RL}}_{\mathrm{w}}$, respectively. First, with the increase of RLw${\text{RL}}_{\mathrm{w}}$, the CPU time increases. Meanwhile, the decrease of RS$\text{RS}$ is also an important factor for the increase of CPU time. We find that tG is almost equal for the three cases in the same numerical example (e.g., ellipsoid example is about 41000s$41\ 000\ \text{s}$), this indicates that the global iteration time is independent of the polymer placement and is only related to the molecular conformation. Furthermore, by comparing the simulation results of Figure 11a,b, it clearly implies that tG of the ellipsoid example is always large than the prismatic example under same conditions. Accordingly, it is further verified that the CPU time of the global iteration only depends on the calculation format and has no relationship with the position of the polymer. We know that the potential value around the molecular interface is the major calculation consumption in the dielectric inhomogeneous system. It also explains that the CPU time of the dielectric inhomogeneity is much more than the dielectric homogeneity in global iteration. When the curves of all cases are plot in figure, we can clearly observe some differences, the overall trend of the connection between RLw${\text{RL}}_{\mathrm{w}}$ and CPU time is very similar, but the computational efficiency changes at different RLw${\text{RL}}_{\mathrm{w}}$ (note that the difference mentioned above is not related to the numerical example, but only to the spatial location of the polymer, thus the two diagrams are described uniformly): i) for case 1, the computational efficiency is greatly improved in the initial stage, and the CPU time increases slowly. When RLw>35%${\text{RL}}_{\text{w}}&gt;35\%$, the CPU time increases suddenly and rapidly, ii) in case 2, the rate of the change in computational speed is different from that of case 1 (RLw=25%${\text{RL}}_{\text{w}}=25\%$), because the solute region is smaller than the solvent region. When the local iteration domain is expanding, case 2 covers the interface area earlier than case 1. Undoubtedly, the calculation of the potential value around the interface occupies most of the CPU time. Hence, the CPU time grows rapidly when RLw=7.5--25%${\text{RL}}_{\text{w}}=7.5\text{--}25\%$ and then increases slowly. iii) As expectation, this is further verified in case 3, which shows a similar phenomenon to case 2. In the RLw<37%${\text{RL}}_{\text{w}}&lt;37\%$ stage, the CPU time increases rapidly, which is not only related to the current position of the polymer (the junction between solute and solvent), but also has a certain relationship with the direction of polymer movement. The polymer moves toward the interior of the biomolecule under the action of the applied electric field force and the electrostatic force in the electrolyte solution, so that the speed of the local box covering the molecular region becomes faster. Until the critical value 37% is exceeded, the molecular region is completely covered. As RLw${\text{RL}}_{\mathrm{w}}$ increases continuously, the local box expands the number of nodes in the solvent region, so that the increase rate of CPU time gradually slows down.11Figurea,b) The left axis and right axis denote CPU time t and the relative potential error REU${\text{RE}}_{\mathrm{U}}$ as a function of the relative size of local box (RLw${\text{RL}}_{\mathrm{w}}$), respectively. (a) Calculates an ellipsoid numerical example, (b) calculates a prism numerical example. Marked the CPU time and error for the equal conditions with a same color line, where the solid represents the CPU time and the hollow represents the REU${\text{RE}}_{\mathrm{U}}$, and the intersection point is defined as the optimal parameter point for this group (mark with solid circle).Additionally, as shown in Figure 11a,b, with the increase of RLw${\text{RL}}_{\mathrm{w}}$ and the decrease of RS$\text{RS}$, the value of REU${\text{RE}}_{\mathrm{U}}$ decreases significantly. However, the RLw${\text{RL}}_{\text{w}}$–REU${\text{RE}}_{\text{U}}$ function curves of the three cases seem to have some differences. In the discussion of interface point calculation format, we construct a high‐precision discrete difference scheme via adding differential nodes. Therefore, we consider that the major error originates from the solute and solvent region (7‐point center difference scheme) rather than the region around the interface (high‐precision difference scheme). And we also infer that this is opposite to the change trend of CPU time. Numerical results prove that the initial error of case 3 is the largest, because the polymer is originally located at the solute–solvent junction. The electric potential in the region of major error source is not updated for a period of time, which will lead to an error propagation effect and increase REU${\text{RE}}_{\mathrm{U}}$. For cases 1 and 2, the smaller RLw${\text{RL}}_{\mathrm{w}}$ makes the local box only work in the solvent or solute region, so that error propagation in the solvent or solute region is blocked. Theoretically, according to current simulation data, we can explore the optimal parameter selection for each numerical example, that is, the intersection of the same set of curves, which is marked as solid circle symbol in Figure 11. Obviously, the optimum point does not have the best acceleration effect, but it is the most cost‐effective situation, which is a good balance between computational efficiency and accuracy. For example, in case 1 of the ellipsoid example (Figure 11a), the optimal parameters of the ellipsoid example are RLw=28%${\text{RL}}_{\text{w}}=28\%$, RS=25%$\text{RS}=25\%$, RLd=80%${\text{RL}}_{\text{d}}=80\%$. Certainly, if test more parameter combinations, you may get a better parameter group. In addition, in both Figure 11a,b under the same case, even different RS$\text{RS}$ are selected, the optimal parameters almost remain within a certain interval. For instance, in case 1 of Figure 11b, RLw=15--20%${\text{RL}}_{\text{w}}=15\text{--}20\%$ is the optimal choice, while in case 3, the optimal RLw${\text{RL}}_{\mathrm{w}}$ takes the value between 19% and 23%.As noted above, we now examine the performance of the local iteration method for simulating the motion of charged polymer under two different shape of solute molecular models. Generally, one find that local iteration performs well when the type of the polymer is clustered or planar. However, once the shape of the polymer is scattered or diagonal, the effect of local iteration will be greatly reduced. This causes certain challenges to the universality of our algorithm. In addition, the computational difficulties caused by the special iteration scheme of the nodes around molecular interface bring another challenge to the total CPU time. Therefore, in order to improve the local iteration method, we learn from the original optimization ideas of the traditional local iteration (TLI) to further reduce unimportant iterations time. We do this by constructing a compact base box that tightly wraps the polymer. That is, the configuration of the initial base box is the same as the polymer shape, which minimizes the unnecessary iteration area. Another advantage of this approach is that even the shape of the polymer is constantly changing as it moves, the adaptive banded local iteration (BLI) can locate the current shape of the polymer in real time, and expand RLw${\text{RL}}_{\mathrm{w}}$ units in six directions to form local box. Note that we cancel the warning value RLd${\text{RL}}_{\mathrm{d}}$ in the BLI method, because many numerical experiments show that it does not spend much time to re‐divide the local box frequently.Finally, we simulate two polymers of different shapes under two models with different molecular structures, using a banded local iterative approach. Figure 12a,b display the CPU time (t) as a function of the relative length of the local box (RLw${\text{RL}}_{\mathrm{w}}$) at the relative global iteration step interval RS=25%$\text{RS}=25\%$. With the decrease in RLw${\text{RL}}_{\mathrm{w}}$, the iteration domain reduce, which results in a decrease in CPU time. On the other hand, REU${\text{RE}}_{\mathrm{U}}$ decreases with increasing RLw${\text{RL}}_{\mathrm{w}}$. In terms of error, when RLw<25%${\text{RL}}_{\text{w}}&lt;25\%$, the REU${\text{RE}}_{\mathrm{U}}$ of the ellipsoid numerical example is always lower than that of the prism example. This result coincides with our aim to construct a high‐precision differential scheme by adding differential nodes at the disjoint interface points. Surprisingly, compared to global iteration, BLI method has better acceleration efficiency than TLI method, accompanied by a decrease in accuracy. As shown in Figure 12a, we first consider a diagonal type of polymer (3D), in which the position of the polymer is at the intersection of the solute and solvent. The effects of BLI method on CPU time and accuracy in different numerical examples are evaluated. In this situation, we only increase the simulation efficiency of the TLI method by about three times compared to that of the global iteration due to the limitations of its strategy and polymer shape. As from the inset of Figure 12a, we can imagine the diagonal polymer as a 3D diagonal matrix with a large number of “zero elements,” and try to get faster computation speed by cutting these unnecessary nodes. In Figure 12a, we observe that the computation speed of BLI method is ≈25 times faster than the global iteration (tG${t}_{G}$ is about 34000--42000s$34000\text{--}42000\ \text{s}$). Now, the TLI method takes almost nine times as much CPU time as the BLI method. We also focus on the variation in error. Under the same conditions, the error of the BLI method is higher than that of the TLI method, but is still in the same order of magnitude (10−4), which demonstrates that the BLI has a higher benefit. Moreover, we wanted to verify the performance of the BLI method with a polymer of scattered type, which is also put in the solute region (similar to case 2). We can see the powerful acceleration advantage of BLI method in Figure 12b. When RLw=7.5%${\text{RL}}_{\text{w}}=7.5\%$, the computational efficiency improves 29 times over the TLI method, reaching a peak. At the same time, REU${\text{RE}}_{\mathrm{U}}$ increases from 2.48e−4 to 5.15e−4. We find that the computational efficiency of BLI method can be increased by at least 3.3 times, while the TLI method can only improve the computational efficiency at most 2.8 times, which shows that the BLI method indeed improves the computational efficiency. Otherwise, the results demonstrate that BLI method not only improve the simulation efficiency for polymers with special shapes, but also for case 2 in Figure 9 (where the performance of TLI method is not excellent, only 5.6 times at most). By reducing the iteration area, the local box contacts the interface more slowly and results in improved simulation efficiency.12Figurea,b) The left axis and right axis denote CPU time t and the relative potential error REU${\text{RE}}_{\mathrm{U}}$ as a function of the relative size of local box (RLw${\text{RL}}_{\mathrm{w}}$), respectively. Marked the CPU time and error for the equal conditions with a same color line, where the solid represents the CPU time and the hollow represents the REU${\text{RE}}_{\mathrm{U}}$, and the intersection point is defined as the optimal parameter point for this group (mark with solid circle). The inset is sketched a simple 3D coarse‐grained model of the polymer. (a) simulates a diagonal type polymer, and set it at the junction of solute and solvent, (b) simulates a scattered type polymer, and put it on the solute region.We further extract the optimal parameters from two examples, which can be seen from Figure 12a,b. The current numerical results show that the optimal banded local box size is RLw=34%${\text{RL}}_{\text{w}}=34\%$ in the prismatic example and RLw=32%${\text{RL}}_{\text{w}}=32\%$ in the ellipsoid example for simulating the motion of diagonal polymers. Furthermore, when the shape of the polymer is scattered, the optimal band local box relative size is selected as 34% and 30%, respectively. By comparing and analyzing Figure 12a,b, we find that for the polymer of different shapes, if use the same local iteration algorithm, the optimal parameter ranges are very similar. When the TLI method is used, the optimal RLw${\text{RL}}_{\mathrm{w}}$ is ≈30%$\approx 30\%$ to 34%. However, for the BLI method, setting RLw=10--15%${\text{RL}}_{\mathrm{w}}=10\text{--}15\%$ can achieve the best balance between computational precision and computational cost. This simulation result implies that the distribution interval of the optimal parameters appeared to be unaffected by the shape of the polymer and is more related to the type of local iteration method.ConclusionsIn this work, using a hybrid simulation approach that combines a finite difference method and a Brownian dynamics model (FD–BD) for the polymer, we study the diffusion of charged polymers in the dielectric inhomogeneity system. We believe that further optimization of the hybrid FE‐BD approach might include the following. Perhaps the main difficulty is need to solve an equation based on a discontinuous medium in the entire domain, that is, the Poisson equation with variable dielectric coefficient. We propose a new computational scheme based on FDM with second‐order accuracy and effectively handles with the interface problem posed by such equations, achieving a smooth transition of the potential value from one region to another. The numerical results show that the improved scheme has a second order of convergence. In addition, we validate the validity of the global SOR iteration method based on the improved interface scheme. It is an excellent numerical tool for simulating charged polymer motion since its computational efficiency is better than the traditional Coulomb formula. However, in simulations that involve millions or even tens of millions of time steps, such computational efficiency is far from adequate. In order to avoid the iterative section to occupy too much CPU resources, the local iteration method is used to accelerate this process. In the simulations, we modeled numerical examples of two different biomolecule configurations (ellipsoid and prism) and put the polymer of same shape in the solvent region, the solute (biomolecule) region, and the solute–solvent interface, respectively. As expected, we find that both the decrease in RLw${\text{RL}}_{\mathrm{w}}$ and the increase in RS$\text{RS}$ can significantly reduce CPU time and thus improve computational efficiency, while the corresponding cost is the increase in error. We also find that although the overall trends of computational efficiency is very similar for the three cases, they exhibit slight differences in the inflection points of the magnitude changes: as local box continues to expand, computational efficiency is expressed as the solvent region > the interface > the internal region of solute in the early stages, while in the late stage, the internal region of solute > the interface > the solvent region. In general, the variation in the magnitude of computational efficiency depends on the rate at which the local box covers the molecular interface. In addition, our simulation results demonstrate that the error of local iteration method increases as the distance between the polymer and the molecular interface decreases.We further discuss the optimal parameter sets in the numerical experiments of this paper, which are chosen by combining both computational speed and accuracy for the sake of reasonableness. For example, RLw=28%${\text{RL}}_{\text{w}}=28\%$, RS=25%$\text{RS}=25\%$, RLd=80%${\text{RL}}_{\text{d}}=80\%$ are the optimal set of parameters for case 1 in the ellipsoid numerical example, where the computational efficiency is improved by 4.5 times compared to the global iteration, and REU${\text{RE}}_{\mathrm{U}}$ is controlled below 9.5e−5. Importantly, in all cases of the ellipsoidal example, RLw${\text{RL}}_{\mathrm{w}}$ is basically concentrated between 10% and 28% to ensure the best results of simulation, while in the prismatic example it takes values from 12% to 23%. Accordingly, the optimal parameters are determined based on different examples. In this way, the local iteration method effectively solves the computational complexity of global iteration in the dielectric inhomogeneous system, which saves a significant amount of CPU time.However, we observe that the traditional local iteration (TLI) method does not have significant acceleration effect for the polymers of two particular shapes. Therefore, we propose the banded local iteration (BLI) method, which is better adapted to the effects of deformation of the polymer during motion. Our results indicate that the BLI method additionally improves the computational efficiency over the TLI method. In particular, when the polymer shape is scattered, the BLI method improves the computational efficiency by a factor of 29 over the global iteration, which is much higher than that of the TLI method (4.5 times). We hope this work can be helpful for accelerating the solution of electrostatic interaction forces under dielectric inhomogeneity, and provide theoretical guidance for the optimization of iterations for polymer chains with different geometries in large electrostatic systems (containing a large number of charged particles or huge computational domain). Furthermore, we intend to study the performance comparison of local iteration method and P3M, Ewald summation in high concentration salt solution in our further work. In future studies, it would be very fascinating to understand the motion of charged particles under the dielectric inhomogeneous system from a different angle, for example, by searching the energy minimum of the system to achieve a mutual equilibrium, which will help us to solve more complex molecular with singularities.AcknowledgementsThis work is supported by the National Natural Science Foundation of China (Grant No. 12071363 and 11701433), the Science and Technology Major Project of Hubei Province under Grant 2021AAA010, and the National Key Research and Development Program of China (No. 2020YFA0714200).Conflict of InterestThe authors declare no conflict of interest.Data Availability StatementThe data that support the findings of this study are available from the corresponding author upon reasonable request.D. L. Ermak, J. A. McCammon, J. Chem. Phys. 1978, 69, 1352.B. A. Luty, J. A. McCammon, H. X. Zhou, J. Chem. Phys. 1992, 97, 5682.H. X. Zhou, J. Chem. Phys. 1990, 92, 3092.S. R. McGuffee, A. H. Elcock, Plos Comput. Biol. 2010, 6, e1000694.S. H. Northrup, S. A. Allison, J. A. McCammon, J. Chem. Phys. 1984, 80, 1517.S. H. Northrup, J. O. Boles, J. C. Reynolds, J. Phys. Chem. 1987, 91, 5991.S. Lu, A. P. Michailova, J. J. Saucerman, Y. Cheng, Z. Yu, T. H. Kaiser, W. W. Li, R. E. Bank, M. J. Holst, J. A. McCammon, IEEE Eng. Med. 2009, 28, 46.I. Nakamura, Z. G. Wang, Soft Matter. 2012, 8, 9356.I. Gunkel, T. Thurn‐Albrecht, Macromolecules 2012, 45, 283.C.‐L. Ren, I. Nakamura, Z.‐G. Wang, Macromolecules 2016, 49, 425.A. A. Teran, N. P. Balsara, J. Phys. Chem. B 2014, 118, 4.A. Abrashkin, D. Andelman, H. Orland, Phys. Rev. Lett. 2007, 99, 077801.A. Hildebrandt, R. Blossey, S. Rjasanow, O. Kohlbacher, H.‐P. Lenhof, Phys. Rev. Lett. 2004, 93, 108104.B. Lu, Y. Zhou, M. Holst, J. McCammon, Commun. Comput. Phys. 2008, 3, 973.A. Y. Toukmaji, J. A. Board Jr, Comput. Phys. Commun. 1996, 95, 73.M. Deserno, C. Holm, J. Chem. Phys. 1998, 109, 7694.T. Darden, D. York, L. Pedersen, J. Chem. Phys. 1993, 98, 10089.V. Ballenegger, J. J. Cerda, C. Holm, J. Chem. Theory Comput. 2012, 8, 936.H. M. Antia, Numerical Methods for Scientists and Engineers, Third Edition, Springer, 2012.G. D. Smith, Numerical Solution of Partial Differential Equations: Finite Difference Methods, Oxford University Press, Oxford 1985.D. Young, T. Am. Math. Soc. 1954, 76, 92.J. Zhang, J. Comput. Phys. 1998, 143, 449.X. I. A. Yang, R. Mittal, J. Comput. Phys. 2014, 274, 695.T. Fukuchi, AIP Adv. 2014, 4, 055312.W. F. Spotz, G. F. Carey, Numer. Meth., Part D E 2015, 12, 235.G. Sutmann, B. Steffen, J. Comput. Appl. Math. 2006, 187, 142.H. Wang, Y. Zhang, X. Ma, J. Qiu, Y. Liang, Comput. Math. Appl. 2016, 71, 1843.J. E. Adsuara, I. Cordero‐Carrion, P. Cerda‐Duran, M. A. Aloy, J. Comput. Phys. 2016, 321, 369.G. A. Huber, Y. Miao, S. Zhou, B. Li, J. A. McCammon, J. Chem. Phys. 2016, 144.P. Bauler, G. A. Huber, J. A. McCammon, J. Chem. Phys. 2012, 136, 164107.R. Egan, F. Gibou, J. Comput. Phys. 2017, 346, 71.J. D. Towers, J. Comput. Phys. 2018, 361, 424.M. K. Gilson, M. E. Davis, B. A. Luty, J. A. McCammon, J. Phys. Chem. 1993, 97, 3591.Z. Zhou, P. Payne, M. Vasquez, N. Kuhn, M. Levitt, J. Comput. Chem. 1996, 17, 1344.I.‐L. Chern, J.‐G. Liu, W.‐C. Wang, Methods. Appl. Anal. 2003, 10, 309.H. Li, B. Lu, J. Chem. Phys. 2014, 141, 07B609_1.X. Liu, Y. Qiao, B. Lu, SIAM J. Appl. Math. 2018, 78, 1131.D. Xie, J. Comput. Phys. 2014, 275, 294.Y. Zhou, S. Zhao, M. Feig, G.‐W. Wei, J. Comput. Phys. 2006, 213, 1.Y. Zhou, G.‐W. Wei, J. Comput. Phys. 2006, 219, 228.Y. Zhou, M. Feig, G.‐W. Wei, J. Comput. Chem. 2008, 29, 87.K. Xia, M. Zhan, G.‐W. Wei, J. Comput. Phys. 2011, 230, 4588.W. Geng, S. Zhao, J. Comput. Phys. 2017, 351, 25.H. Feng, G. Long, S. Zhao, J. Comput. Appl. 2019, 361, 426.A. Lee, W. Geng, S. Zhao, J. Comput. Phys. 2021, 426, 109958.H. Lin, Z.‐T. Lei, M.‐M. Ding, H.‐J. Wang, T.‐F. Shi, Chinses. J. Polym. Sci. 2020, 38, 665.F. Booth, J. Chem. Phys. 1951, 19, 391.H. Warner, Ind. Eng. Chem. Fund. 1972, 11, 379.D. Frenkel, B. Smit, Understanding Molecular Simulation, Academic Press, New York 2002.M. Allen, D. Tildesley, Computer Simulation of Molecular Liquids, Clarendon Press, Oxford 1987.G. C. Maitland, Intermolecular Forces: Their Origin and Determination, Oxford University Press, Oxford 1981.

Journal

Advanced Theory and Simulations – Wiley

Published: Apr 1, 2023

Keywords: banded local iteration; Brownian dynamics; dielectric inhomogeneity; electrostatic interaction; finite difference method

There are no references for this article.