Fast Computation of Electrostatic Interactions for a Charged Polymer with Applied Field

Using a hybrid simulation approach that combines a finite difference method with a Brownian dynamics, we investigated the motion of charged polymers. Owing to the fact that polymer-solution systems often contain a large number of particles and the charged polymer chains are in a state of random motion, it is a time-consuming task to calculate the electrostatic interaction of the system. Accordingly, we propose a new strategy to shorten the CPU time by reducing the iteration area. Our simulation results illustrate the effect of preset parameters on CPU time and accuracy, and demonstrate the feasibility of the “local iteration” method. Importantly, we find that the increase in the number of charged beads has no significant influence on the time of global iterations and local iterations. For a number of 80 × 80 × 80 grids, when the relative error is controlled below 1.5%, the computational efficiency is increased by 8.7 times in the case that contains 500 charged beads. In addition, for a number of 100 × 100 × 100 grids with 100 charged beads, the computational efficiency can be increased up to 12 times. Our work provides new insights for the optimization of iterative algorithms in special problems.


INTRODUCTION
The diffusion of charged polymers is often the rate-determining step in many biological processes. Stochastic methods (such as Brownian dynamics) and continuum methods (such as finite difference method) are two main computational methods for studying diffusion. [1−6] In the process of diffusion, electrostatic force plays an important role in polymer interactions. At present, electrostatic properties of polymers have been extensively studied in the interdisciplinary areas of biomolecule, computational mathematics, physics, and chemistry. There are many methods for calculating electrostatic interactions. However, the existing methods (such as particle-mesh Ewald (PME), Ewald summation, particle-particle particle-mesh (P3M), Coulomb force formula, etc.) all require a lot of CPU time and use large computer memory. [7−11] When a huge three-dimensional simulation system is taken into account, it is unsuitable to use these methods to calculate electrostatic forces owing to the serious challenge in computing time. Therefore, an efficient method is urgently needed for saving computing time. The numerical me-thod is more advantageous, because it sacrifices some accuracy to get a great increase in the calculation speed. At the same time, the numerical method has a weaker dependence on computer memory. [12] Hence, through combining the advantages of traditional methods and numerical methods to simulate the motion of charged polymers is an effective and accurate way, where the numerical methods is mainly used to handle the electrostatic interaction of charged polymers by solving the Poisson equation. [13] This is important because many diffusing molecules in biological systems are charged. The key to solving Poisson equation is how to figure out the potential distribution of electrostatic field quickly and accurately by combining the boundary conditions, and provide accurate force analysis for the subsequent simulation of charged polymer motion. Traditional numerical simulations of three-dimensional problems tend to be computationally intensive because of the requirements on the memory and CPU time to obtain solution with desirable accuracy. Even the most advanced computers may not be able to handle them directly. Thus, in electrostatics, calculating the potential distribution of a given charge density function is crucial. In solving such three-dimensional problems, finite difference method (FDM) is a relatively simple and convenient method in many numerical methods, which has the mature theory and the optional accuracy. [14] When dealing with the regular area, it is easier to program and im-plement on computer than the finite element method (FEM) and the finite volume method (FVM). [15] The principle of FDM is to discretize Poisson's equation. [16−20] Most of the discretized Poisson equations will give a series of linear equations with a large number of unknowns and diagonally dominant coefficient matrices. By using Jacobi, Gauss-Seidel, successive over relaxation (SOR), or Multigrid iterative methods, however, the unknown solutions at the corresponding grid nodes can be approximated iteratively.
Since the iterative algorithm takes up too much CPU time, almost all optimizations are for iteration algorithms. Kamboh et al. have proposed the parallel Jacobi method (PJ), which can reduce the computing time for large data moderately. [21−23] It is superior to the sequential Jacobi (SJ) and the sequential Gauss-Seidel (SGS). Scheduling Relaxation Jacobian (SRJ) is another method for classical optimization, which greatly speeds up convergence and reduces computation time, and it has certain advantages in dealing with complex boundaries. [24−31] In addition, there are also some studies on Multigrid method, for solving 3D Poisson equation with high accuracy, which bases on Richardson extrapolation method. [32−36] The fourth-order compact difference discretization scheme and V-cycle algorithm of multigrid are used to solve the Dirichlet boundary value problem of 3D Poisson equation. Finally, the expensive computation of 3D Poisson equation is solved very well. This method is characterized by high accuracy and some speed improvement.
Nonetheless, in Brownian dynamics, we need to calculate millions or even tens of millions of time steps, which means that the calculation of electrostatic interactions will take substantial time even the aforementioned methods are fast in calculation. Therefore, in the process of simulating the movement of charged polymer chains, it is necessary for us to find a new way to save time and ensure accuracy for millions or even more calculations.
In this study, we present "local iteration" to accelerate the simulation of charged polymers's motion, which possesses the irrelevant local iteration process and global iteration process. The simulation model and corresponding optimization algorithm are presented in the second section. And in the third section, we systematically discuss the effect of the parameters in "local iteration" on the simulation efficiency and accuracy, which further verifies the feasibility of our optimization algorithm. Finally, comprehensive summarization of the processes and results of the simulation test are presented in the last section, as well as some prospections.

SIMULATION MODEL AND METHODS
As shown in Fig. 1, in a simulation box, we fix the left face of electric potential (blue) as U and the others (gray) as 0, and put into a polymer, which is modelled as a conventional beadspring chain composed of N successive Lennard-Jones (LJ) [37−39] beads and connected through the finitely extensible nonlinear elastic (FENE) potential, [40] U LJ (r) = 4ε[(σ/r) 12 and where σ is the diameter of a bead, ε gives the strength of the LJ κ b potential, represents the spring constant, and b denotes the maximum separation distance between the consecutive beads along the polymer chain.
In addition, we set the charge of beads in the polymer as q and solve the Poisson equation with the boundary electric potential to calculate the electrostatic interactions between charged beads. Therefore, we can employ the Brownian dynamics to simulate the charged polymer diffusion, which includes a noise and a friction term. Having specified the interaction potentials, the BD equation of motion for bead i is where U i is the sum of all the interaction potentials discussed above, μ denotes the friction coefficient, and is the random force related to μ by the fluctuation dissipation theorem , where I is the unit tensor. We integrate the equations of charged polymer motions against a time step Dt. In our model, σ, ε, and m denote the units of length, energy and mass, respectively, which result in a unit of the time as t = t*(ε/m/σ 2 ) 1/2 , the charge as q = q * /(4πε 0 σε) 1/2 , and the electric field as E = E * (4πε 0 σε) 1/2 σ/ε. We present our results in reduced units in which σ = ε = m = 1 and the temperature is k B T = 1.0, and the friction coefficient is set as μ = 0.5τ -1 . Therefore, the dimensionless parameters of the polymer are fixed as = 30 and b = 1.5, the charge of beads in the polymer is q = 10, and the time step is chosen as Dt = 0.0002. We fix the size of the simulation box as L = 80, the electric potential of left surface as U = 4, and the relative dielectric constant ε r = 80.
In this work, we focus on the rapid solution of electrostatic interactions between charged beads in a polymer with the external electric field, which corresponds to solving a threedimensional Poisson equation with: and which can be represented by the sum of a series of Dirac-Delta functions, u is the electric potential, ρ is the charge density as a function of position, ε 0 and ε r are the vacuum permittivity and the relative permittivity of the solvent, respectively, α repres- ents a potential function of the boundary node, and Ω denotes a continuous domain in 3D space with boundary conditions prescribed on its boundary. As shown in Fig. 2, in order to obtain the numerical solution of u i,j,k , we use the traditional 7-point central difference to discretize Eq. (4), where i, j, and k represent the node location in the discretized mesh and Δx, Δy, and Δz represent the mesh spacing in the x, y, and z directions (Δx = Δy = Δz = 1). The charge density of each node is estimated by a linear interpolation scheme from charged beads in the polymer. After discretization, Eq. (6) represents the linear equations, which can be expressed by a standard matrix form Au = B. Thus, we can obtain u i,j,k as: where ρ i,j,k is the charge density distribution function. In this work, we use Successive Over-Relaxation (SOR) method for solving the above equations, which is one of the most efficient iterative methods. Compared with Gauss-Seidel and Jacobi iteration methods, this method shows great advantage in calculation speed for large linear equations, whose iterative format can be written as: where λ represents the relaxation factor and n denotes the number of iteration steps. Obviously, when λ = 1, the SOR is the Gauss-Seidel iteration. The optimal relaxation factor can be obtained by [30] λ optimal = 2/[1 + sin( where Mig is the minimum number of grid nodes in three directions. For 80 × 80 × 80 grids, Mig = 81, λ optimal = 2/{1 + sin[π/(81 -1)]} ≈ 1.924. Thus, we can obtain the potential of each nodes, and further calculate the electric field in the position of charged beads by the same linear interpolation scheme, which means that we can calculate the electrostatic interactions between charged bead in the polymer with the external electric field. In order to compare the CPU time among SOR, Gauss-Seidel, and Jacobi iteration methods, we display the CPU time for these methods in Fig. 3, which includes the direct solution of Coulomb force between charged beads. The comparison results show that with the increase of charged beads, the CPU time of these methods all exhibits an increase. The SOR method is different from direct solution of Coulomb force, which is insensitive to the increase of the number of charged beads. For current grid numbers (80 × 80 × 80), we find that when the number of charge beads Num is about 75, the CPU time for SOR is almost equal to that of direct solution of Coulomb force, which indicates that the SOR method is more suitable for large scale charged bead system.
In fact, in Brownian dynamics, we often need to calculate millions or even tens of millions of time steps, which means that the calculation of electrostatic interactions will take substantial time. If we only consider one charged polymer in solvents with uniform dielectric constant, within a certain time interval, the effect of the motion of charged polymer on the potential of grids far from the charged polymer can be neglected, especially in the presence of the external electric field. Thus, we can save most of the iteration time for grids far from the charged polymer, whose value of potential has almost no change within a certain time interval. This means that we can use "local iteration" instead of "global iteration" on the premise of guaranteeing calculation accuracy to save a lot of unimportant iteration time. Specially, in the simulation box, the charged beads in the polymer often occupy only a small area.
Thus, in a short simulation time, we only calculate the potential of adjacent grids and the potential of other grids is fixed, which will be solved after several steps. For example, in Fig. 4, in order to better verify the "local iteration" method, we put the charged polymer in the middle of the simulation box. We first define a smallest "box", which just surrounds the charged polymer. Then we expand L w grids outwards and form a "local box". When the "local box" exceeds the global box in a certain direction, it will not be further expanded. In addition, we set d step as the interval of steps for each global iteration and a warning value (L d ). When the nearest distance between the smallest box and the "local box" equals the warning value, we carry out the global iteration and re-divide the "local box". Therefore, in this work, we mainly study the effects of three parameters on the calculation accuracy and speed, which include L w , d step , and L d . We expect that this method can greatly improve the calculation speed while ensuring the accuracy.

RESULTS AND DISCUSSION
In the numerical experiment, all the codes use double precision arithmetic and run on a personal desktop with Intel Core i7-6700 CPU (3.40 GHz) and 16 GB RAM. We mainly study the independent effects of three parameters (L w , d step , and L d ) on the calculation accuracy and speed. The effect of the number of charged beads in a polymer (denoted as Num) on the calculation is also discussed. Finally, the applicability of the local iteration method is verified by the different numbers of grids. Based on the comparison between the results of local iteration and global iteration, the advantages of local iteration are explained from two aspects: on the one hand, the relative potential error (RE U ) of each charged beads between global iteration and local iteration is calculated through: where θ RMS denotes the root mean square sign, U global and U local represent the potential at each charged beads in the global iteration and the local iteration, respectively. On the other hand, the CPU time of global iteration (denoted as G-t) and local iteration (denoted as L-t) is compared. In order to make L w , d step , L d , and other parameters more universal, we use the relative size for illustration. For example, RL w = L w /L, RS = d step /Step, RL d = L d /L w . Because the charged beads occupy some grids, we ex-tend L w grids in six directions. When RL w = 50%, whether the "local box" completely covers the simulation box depends on the initial position of the charged beads in a polymer. Firstly, we discuss the influence of the value of L w in local iteration on the accuracy and CPU time, which are compared with the result of global iteration. We set the grid size as 80 × 80 × 80, the number of charged beads in a polymer Num = 100, total step number of motion Step = 10000, the relative distance between charged polymer and boundary RL d = 80%, and the relative global iteration step interval RS = 0.25%, 1%, 2.5%, 5%. Under these conditions, we calculate L-t, G-t, and RE U . As shown in Fig. 5(a), with the decrease of RL w , the CPU time decreases, which implies an increase in the simulation efficiency. However, the decrease of RL w has no effect on the global time G-t, because G-t only depends on the grid size and is not affected by the parameters in "local iteration" (such as L w , d step ). Otherwise, with the same RL w , an increasing RS also induces the decrease of CPU time. On the other hand, as shown in Fig. 5(b), both the decrease of RL w and the increase of RS, which are of advantages to the computational efficiency, induce the increase of relative potential error RE U . Obviously, with the increase of RL w , the scope of local iteration is also enlarged, which implies the increase of CPU time consumed by local iteration and further results in the increase of total CPU time. At the same time, enlarged local iteration also makes the solutions approximate that of the global iteration (i.e. a smaller RE U ). Besides, under the same total iteration step, the number of global iterations increases with RS decrease. Hence, a smaller RS indicates more computational consumption in global iteration and fewer cumulative errors, which further induces the increase of total CPU time and the reduction of relative potential error RE U . Note that when RL w is beyond 47%, L-t approaches G-t and RE U approximates zero. This is because RL w > 47% implies that the "box" of local iteration is nearly as large as simulation box, resulting in the similar computational consumption and accuracy. Both the results of CPU time and RE U demonstrate the availability of our "local iteration": when RL w = 6.25% and RS = 5%, the simulation efficiency is increased by 12.5 times compared to that of the traditional global iteration; meanwhile, RE U is merely 1.43%, which is in the acceptable range. Even one wants a greater accuracy, RE U ≈ 0.5% when RS = 0.25%, the smaller RL w (6.25%) also ensures a sufficient efficiency (8 times than global iteration).
In order to further observe and verify the effect of d step on accuracy and speed, we use RS to explain more intuitively. Similarly, under the grid number of 80 × 80 × 80, we set Num = 100, the length ratio of "local box" to "global box" RL w = 6.25%, the relative distance between charged polymer and boundary RL d = 80%, and Step = 10000, 15000, 20000, 30000. Then, we calculated L-t, G-t, and RE U . As shown in Fig. 6(a), with the increase of RS, the CPU time decreases, which implies an increase in the simulation efficiency. In addition, we find that an increasing Step induces the increase of CPU time under the same RS. Otherwise, as shown in Fig. 6(b), the increase of RS, which is of advantage to the computational efficiency, would induce the increase of relative potential error RE U . We know that with the increase of RS, the number of global iterations decreases, which implies the decrease of total CPU time. However, decreasing the value of RS can make the solution approximate that of the global iteration. We find that when RS = 0.25%, the computational efficiency of Step = 10000 and Step = 30000 is increased by 10.1 times and 11.3 times, respectively. The result indicates that with the increase of Step, the computation efficiency of local iteration increases more obviously. In other words, the large number of steps provide more speedup than that of the small number of steps under same conditions. That 1% ≤ RS ≤ 2.5% is the best under these Step is clearly exhibited in Fig. 6(b). This is due to the enough large number of global iterations even under large Step. As a result, RE U is lower than 1.35%, and the computational efficiency can be increased by at least 9.8 times. Even one wants a greater accuracy, we can choose a larger RL w to reduce RE U .
Moreover, when the charged polymer moves towards the boundary of the "box" (the following "box" denotes "local box"), we also set a "warning line" -L d . We still select the number of 80 × 80 × 80 grids, under Num = 100, Step = 10000, RS = 1%, and RL w = 6.25%, 10%, 12.5%, 15%; CPU time and RE U corresponding to different RL d are individually calculated. Owing to different values of L w , the value of L d also varies correspondingly. For example, when RL w = 6.25% (L w = 5), L d can only take 4, 3, 2, or 1. According to Fig. 7(a), we find that the size of L d has almost no effect on the CPU time, because L d is only a parameter for re-dividing the "box" and it does not increase the amount of calculations like L w and d step . As shown in Fig. 7(b), RE U decreases with the increase of RL d . The result indicates that the accuracy increases with the increase of update frequency of "box". This is because the effect of the charged polymer on the potential of external node of the "box" weakens with the increase of the distance between charged polymer and boundary of the "box" and RE U thereby increases. In addition, we find that the RE U /RL d curve becomes steeper with the increase of RL w . In theory, RE U decreases with the increase of RL w . However, we find that the accuracy of RL w = 15% is not as good as that of RL w = 12.5% when RL d value is very small. This phenomenon is caused by the fact that the moving distance for updating the "box" increases with the increase of RL w when RL d is the same, Step = 10000 Step = 15000 Step = 20000 Step = 30000 Step = 10000 Step = 15000 Step = 20000 Step = 30000 Fig. 6 (a) CPU time t and (b) the relative error of electrostatic potential RE U as a function of the relative global iteration step interval RS. Grid size 80 × 80 × 80, the number of charged beads in a polymer Num = 100, the length ratio of "local box" to "global box" RL w = 6.25%, the relative distance between charged polymer and boundary RL d = 80%.
leading to the lower frequency of updating the "box", and the cumulative error further becomes too large. Such phenomenon is especially obvious when RL d ≤ 20%. Then, aiming at the uncertainty of the number of charged beads in polymer Num in practical problems, we choose Num = 100, 200, 300, and 500 under the condition that the number of grids is 80 × 80 × 80, Step = 10000, RS = 1%, and RL d = 80%. As shown in Fig. 8(a), the CPU time decreases with the decrease of RL w , which is similar to Fig. 5(a). Moreover, we find that with the same RL w , the increase of Num leads to the increase of CPU time for both local iteration and global iteration, while G-t of Num = 500 is only 1123 s more than G-t of Num = 100. There are two reasons why the CPU time of large number of charged beads (L Num ) is slightly longer than that of small number of charged beads (S Num ). Firstly, the number of grids occupied by the charged polymer increases with the increasing number of charged beads in a polymer. Therefore, the local iteration range of L Num is larger than that of S Num even though RL w is the same. Secondly, L Num will consume more CPU time in the process of charge interpolation. However, the CPU time consumed by the interpolation calculation process is much less than the iterative calculation pro-cess. As shown in Fig. 8(b), RE U increases slightly with the increase of Num. The reason for the difference in error is that the influence of L Num on the external node potential of the "box" is greater than that of S Num . We find that when RL w is large enough, the influence of L Num on the potential of external nodes is gradually weakened. We also observe that RE U of L Num and S Num approach more obviously around RL w = 28.75%. When RL w is between 36% and 50%, RE U of the four cases of Num is infinitely close. Of course, if we want to analyze a huge number of charged beads, we should expand the "local box" appropriately.
Finally, we use the "local iteration" method for testing the effect of different numbers of grids. Under the conditions of Num = 100, Step = 10000, RS = 1%, and RL d = 80%, we test the G-t and L-t of 100 × 100 × 100, 80 × 80 × 80, 100 × 50 × 50, and 50 × 50 × 50 grids and find that the values of G-t are 46514.8, 23252.4, 6158.3, and 2373.9 s, respectively. As shown in Fig. 9(a), the CPU time decreases with the decrease of RL w , which is similar to Fig. 5(a)   same RL w , the increasing number of grids leads to the increase of CPU time. As shown in Fig. 9(b), with L w = 5, the ratio of L-t/G-t decreases with the increase of grid number, which indicates that the efficiency of acceleration increases with the increase of grid number. The grid of 100 × 100 × 100 increases computational efficiency by 12 times more than the other cases. Our conclusion is that the "local iteration" is more suitable to computationally intensive problems. As shown in Fig. 9(c), the decrease of RL w , which is of advantage to the computational efficiency, induces the increase of RE U . Distinctly, the number of 100 × 100 × 100 grids is much larger than that of other sizes, and the external node increases as the "box" becomes smaller. Therefore, cumulative error caused by external nodes of the "box" leads to the larger RE U . However, we can still control RE U of the huge number of grids such as 100 × 100 × 100 below 3.5% with RS = 1%. In addition, we also find that RE U of 50 × 50 × 50 is even greater than that of 80 × 80 × 80 when RL w ≤ 18%. For smaller simulation boxes, the effect of charge on nodes outside the local iteration is greater, so RL w > 18% is a better choice for the small number of grids such as 50 × 50 × 50 from the Fig. 9(c).

CONCLUSIONS
In this work, basing on the SOR iteration algorithm in the finite difference method, we present a "local iteration" method to accelerate the simulation of three-dimensional diffusion of charged polymers in solution. The effects of three parameters (L w , d step , and L d ) in the "local iteration" method on the computational efficiency and accuracy are further tested. Our simulation results demonstrate that both the decrease of RL w and RS could dramatically increase the computational efficiency, while the corresponding price is the decrease of accuracy. Otherwise, the increase of RL d could enhance the computational accuracy and has no influence on the calculative expense. Therefore, in order to reduce the effect of charged beads on the potential of external nodes of the "box", we should make RL d as large as possible. Accordingly, for a simulation system with Num = 100 and grid size 80 × 80 × 80, the "local iteration" method could present an approximate 12-fold calculating speed in comparison with that of the global iteration method, where the relative potential error (RE U ) is controlled below 1.5%. Moreover, the increased complexity of the system (i.e., the number of charged beads) has little influence on the acceleration effect: when the same accuracy is ensured, the "local iteration" method could also reach 8-fold calculating speed for a system with Num = 500 and grid size 80 × 80 × 80. However, the enlargement of simulation system will increase the relative potential error. In order to ensure the calculating accuracy, a larger RL w should be used (RL w ≥ 25%), while the acceleration effect will decrease to only 3-fold (the relative error is controlled below 2%). Hence, the "local iteration" method needs the optimization for special problems, which includes the grid number and three parameters. This could significantly accelerate the simulation processes and guarantee precision control, thereby efficiently accelerating the simulation processes. In practical problems, we may also encounter inhomogeneous media. Therefore, it is very meaningful to further optimize different types of three-dimensional Poisson equation algorithm.   The relative error of electrostatic potential RE U as a function of the length ratio of "local box" to "global box" RL w . Simulation box with dimensions of (100, 100, 100), (80, 80, 80), (100, 50, 50), (50, 50, 50), grid spacing Δx = Δy = Δz = 1. The number of charged beads in a polymer Num = 100, the total step number of motion Step = 10000, the relative global iteration step interval RS = 1%, the relative distance between charged polymer and boundary RL d = 80%.