Theory of Local Sampling: An Effective Sampling Technique in Monte Carlo and in Molecular Dynamics Simulations with Achieved Reduction in Computation Time Up to 67%

1. Introduction

Molecular simulations provide valuable insights at atomic/molecular level of phenomenon like multiphase systems, dissolutions, mixing systems [1]. These insights help in modeling systems with understanding from atomic/molecular level. Both MD and MC simulations collect atomic configurations that can lead to evaluate thermodynamic properties where MD simulations collect configurations over time steps and in MC simulations collect configurations by making random moves which can accept/reject the move based on acceptance rule.

Monte Carlo (MC) [2] and molecular dynamics (MD) [3] simulations were first developed in 1950s. Liquid Argon MD simulation was first developed by Rahman [4] in 1964, using Lennard-Jones potentials for Argon atoms. Many studies [5,6,7] leads to the development of MD simulation studies followed with algorithmic improvements. NVT and NPT ensembles in MD simulations were developed later by Andersen [8], Parrinello and Rahman [9,10] followed by development of Nosé-Hoover thermostat [11,12,13]. In MC simulations, collecting configurations by making random moves with accepted/rejected moves with acceptance rule technique was first introduced in the early 1950s at Los Alamos [14]. The Metropolis algorithm introduced random moves with accepted/rejected with probability so that low-energy states are sampled with probability [14]. Hastings [15] in 1970 later developed the algorithm called as Metropolis-Hastings method that was widely used today. Later developments in MC algorithms like Gibbs ensemble Monte Carlo method [16] and configurational-bias Monte Carlo method [17] together with modern high performance computing power have made MC simulations used for studying vapor/liquid systems and biological systems [18].

Calculation of nonbonded interactions has been a bottleneck for conventional molecular simulations (for MD and for MC simulations) that are limiting the system size and collecting the volume of atomic configurations. For all atom model, MD simulations last up to several nanoseconds (ns) with maximum number of atoms in the order of ${10}^6$ [1]. To perform molecular simulations at higher time scales from nanoseconds to microseconds, development of new theories in molecular simulations are required. With increasing demand for energy supplies, new developments in MC and in MD simulations are needed, that reduce computation cost for each MD or MC step that further reduce cost for maintaining high performance computing facilities.

In this study, Theory of Local Sampling (TLS) was introduced and applied to molecular simulations to reduce computation time in MC and in MD simulations. Theory of Local Sampling (TLS) was first introduced in 2019, in my previous studies (as future work) [19], with noticed reduction in computation time in MC and in MD simulations. TLS can also be applied to other fields where properties are calculated using global regions, like Computational Fluid Dynamics, Stock Market, ab initio molecular simulations, and quantum calculations. Another technique introduced in this study was using curve fitting methods to predict force/acceleration acting on atoms in MD simulations using previous time steps instead of calculating force/acceleration from conventional methods, which reduces computational load. Higher-order function approximations and curve fitting techniques [20] can be implemented to reduce significant computation time.

In the remainder of this paper, working principle of Theory of Local Sampling (TLS) will be introduced in the next section followed by "Simulation Details" of simulations performed in this study. Detailed analysis were presented in "Results and Discussion" section regarding implementation of TLS in MC and in MD simulations, followed by Conclusion section.

2. Introduction to Theory of Local Sampling (TLS)

Theory of Local Sampling (TLS) is an efficient sampling technique in calculating nonbonded interactions that reduces total computation time in MC and in MD simulations. Conventional way of calculating nonbonded interaction for a particular atom includes calculating pairwise van Der Waals and Coulomb interactions for all neighboring atoms with in cutoff radius and adding tail correction beyond cutoff radius. In Theory of Local Sampling where, nonbonded interactions will be calculated by using selective local region around the atom, that carries and contains information of global region. Instead of considering global/all region around an atom, considering local regions was more optimized way for calculating nonbonded interactions, in MC and in MD simulations, as shown in Figure 1. Local region in TLS can be a sphere around an atom, with inner radius (x₁) and outer radius (x₂) (called as "Bin Width" in this study). In MC simulation, considering local region, that carries information of global region that will reflect in calculating "acceptance rule". Similarly, in MD simulations, calculating force/acceleration acting on an atom based on local region gives good agreement with conventional MD simulations, as reported in this study. TLS technique helps in reducing total computation time in calculating nonbonded interactions for the systems at each simulation step. The number of arrows in Figure 1 which represents pairwise potential calculations required in MC and MD simulations, decreased significantly between conventional/global sampling and local/optimized sampling.

Main principle involved in TLS is any small change in global region of a system will reflect in the considered local region and can model the whole system by correlating between global and local region, as shown in Figure 1. TLS is considered to be efficient sampling technique where properties can be modeled by using local sample than sampled globally in a conventional way. Other types of local regions can be considered like a cube or cylindrical ring around the atom, or wavy surface of cos/sin function around the atom, depending on the system, that can give more optimized local sampling in TLS.

Theory of Local Sampling in molecular simulations can be implemented as follows:

• Bin width around an atom (with friction in RDF plots, at the end of bin width).
• Bin width around an atom with varying bin size randomly at each time step (without friction in RDF plots).
• Randomly chosen samples around atom with or without bin width around an atom.
• Considering multiple local regions instead on one local region.

Correlation between local sampling and global sampling can be a captured by modeling the amount of error estimated between global and local region property values. Neural networks and artificial intelligence (AI) are some of the tools that can be integrated in TLS in molecular simulation that will take simulation to higher time scale from nanoseconds to microseconds in MD simulations. With depleting fossil fuel reserves and increase in global warming effects in the past few decades, efficient methods in saving energy is very much needed for power requirements of high performance computing facilities to run molecular simulations for bigger systems.

In this study, results of Theory of Local Sampling (TLS) were reported for MC simulations, MD Simulations, "Curve Fitting in MD simulations" and "TLS & Curve Fitting MD simulations" with significant reduction in computation time. By applying TLS in MC simulations, reduction in computation time was achieved up to 37% and in MD simulations, up to 67% reduction was achieved, as shown in the following sections.

Another approach in molecular simulations to reduce computation time is using "Local Function" instead of conventional potential functions for calculating potential acting on an atom. Idea behind this approach is local function for potential calculation requires less computational load, where conventional potential functions demands higher computational load. Later exact force/potential acting on atom can be deduced by extrapolating/correlating the local function with conventional potential functions. TLS sampling techniques can be applied to systems where properties of single atom/point were estimated by global sampling around the atom/point considered. Systems like ab initio molecular simulations, quantum simulations, Computational Fluid Dynamics (CFD), and dynamics of stock market, where global properties can be estimated by considered local population/local sampling and scaling/correlating back to global values.

3. Simulation Details

Argon system of 500 atoms was used to study the implentation of Theory of Local Sampling (TLS) in Monte Carlo and in Molecular Dynamics simulations. MC and MD simulation codes for Argon atoms were developed using fortran language. 500 Argon atoms were considered in a simulation box with density of 1400 kg/m³. Lennard Jones parameters for argon atoms for $\sigma$ as 3.4 Åand $ϵ$ as 1.65*10⁻²¹ Joules were considered in this study [21]. For nonbonded interactions, cutoff radius of 14 Åwas used and tail correction was applied to atoms beyond the cutoff. Periodic boundary conditions were applied in all directions in a simulation box using Minimum Image Convention algorithm. For MD simulations, Velocity Verlet integrator [22] was used with initial velocities given to argon atoms using Gaussian distribution at 86 K. Velocity rescaling algorithm [23] was implemented at temperature of 86 K. Dimensionless time step of 0.01 was used in MD simulations. Bin width (x₁* - x₂*) for TLS is defined as a local region around an atom with inner radius of x₁* and outer radius of x₂*, making it a sphere region around an atom, as shown in Figure 1. Radial distribution function (RDF) plots and all other variables reported in this study are dimensionless variables (x* = r/$\sigma$) with dimensionless bin widths (x₁* - x₂*) that are considered for local regions in TLS. For MC simulations, dimensionless distance (dx* = dr/$\sigma$) of 0.1 was used to move the atoms randomly at each MC cycle (trail move). Both MC and MD simulations were performed for 100,000 steps to calculate the required computation time. Computation time for both MC and for MD simulations were reported in this study using single processor. Average properties like potential energy, kinetic energy, and temperature of the argon system were calculated for the last 2000 steps for MC and MD simulations with dimensionless values reported in this study. For "Curve Fitting Algorithm in MD simulation", computational load for calculating potential energies at fourth step is excluded.

4. Results and Discussions

4.1. Theory of Local Sampling in Monte Carlo Simulations

TLS was applied to MC simulations (known as "TLS MC") for Argon atoms, for calculating potential acting on an atom by considering local region instead of global region (all region around the atom) with various bin width, as shown in Table 1. Bin width (x₁* - x₂*) were considered with minimum bin width of (0.5 - 1.5) to maximum bin width of (0.5 - 2.0). For calculating "acceptance rule", to accept or reject a random move of an atom, potential energies of new trail move (E₂) and potential energy before trail move (E₁) were calculated in the local region. The differences of conventional potential energies values, in acceptance rule, will be reflected (in lower magnitude) in the considered local region, which is the basic principle of Theory of Local Sampling. Maximum of 37% of reduction in computation time was achieved for bin width of (0.5 - 1.5) for "TLS MC" simulation when compared with conventional MC simulation, as shown in Table 1. Conventional MC simulation takes 579.6 seconds to complete 100,000 MC cycles, where as "TLS MC" simulation takes a minimum of 363.8 seconds for the same Argon system, with reduction in computation time 37%. Percentage of reduction in computation time decreases with increase in bin width of local region from (0.5 - 1.5) to (0.5 - 2.0) for "TLS MC" simulation.

Potential energies calculated for conventional MC, "TLS MC" simulation with bin width of (0.5 - 1.5) and (0.5 - 2.0) was shown in Figure 2. It was observed that atomic configurations of "TLS MC" were different from conventional MC simulation as potential energies having unique atomic configurations over MC cycles. Observation of Figure 2 shows that "TLS MC" simulation with bin width of (0.5 - 2.0) and conventional MC simulation have merged potential values after 1300 MC cycles. Average PE* for bin (0.5 - 1.5) is lower than conventional MC simulation as shown in Figure 2 and Table 1.

Radial distribution function plots were considered for comparing atomic configuration distribution at equilibrium stage for conventional MC and "TLS MC" simulation, as shown in Figure 3. For "TLS MC" simulation, RDF plots shows a friction or discontinuity in atom distribution at the end of bin width, as shown in Figure 3 for bin widths of (0.5 - 1.5) and (0.5 - 2.0). "TLS MC" simulation follows maximum and minimum peaks of first, second, and third solvation shells, when compares to conventional MC simulations. One of the reason for having friction in atoms distribution around local region bin width is "Acceptance rule" is considered only for the local region and atoms in remaining regions have to adjust itself, creating a shock wave or friction at the end of outer radius of bin width (x₂*). To reduce the shock or friction in atom distribution at the end of bin width, outer radius bin width size (x₂*) was varied randomly at the magnitude of 0.4 units at each trail move of atoms in each MC cycle. Variation of outer radius bin width size at each trail move reduces shock in atom distribution and atoms will adjust itself at all regions instead of only at local region bin width. Table 2 shows "TLS MC" simulation with varying bin width randomly from (0.5 to 1.3-1.7) to (0.5 to 2.0-2.4), with maximum reduction in computation time of 36%. Figure 4 show RDF plots comparison for conventional MC simulation to "TLS MC" simulation with varying bin width size, with no friction in RDF curves unlike in Figure 3. It is clear that, varying bin width at each trail move of an atom, reduces shock/friction in equilibrium atomic configurations as shown in Figure 3 and Figure 4.

4.2. Theory of Local Sampling in Molecular Dynamics Simulations

Theory of local sampling (TLS) was applied to MD simulation (TLS MD) of Argon atoms for calculating force acting on an atom that further used in calculating acceleration of an atom that was used in Velocity Verlet Algorithm at each time step. For "TLS MD" simulations, local regions were selected with inner and outer radius (x₁* - x₂*) (bin width) for each atom, for calculating force acting using Lennard-Jones potential. Conventional MD simulations for 500 Argon atoms were performed and compared with "TLS MD" as shown in Table 3. Maximum of 60% reduction in computation time was achieved for "TLS MD" with bin width of (0.75 - 1.25) when compared with conventional MD simulation. Conventional MD simulation takes 324.1 seconds to complete for 100,000 steps using single processor, while "TLS MD" simulation takes minimum of 130.6 seconds (60% reduction in total computation time) for the same system. For bin width (0.75 - 1.25), average potential energy slightly differs between conventional MD and "TLS MD" with values -2996.5 and -2936.3 (dimensionless potential energies, PE*), respectively. As bin width increases from (0.75 - 1.25) to (0.75 - 2.2), average potential energy decreases and came closer to conventional MD simulation energies, as shown in Table 3. Temperature of "TLS MD" and conventional MD maintain same consistency irrespective of bin width chosen for "TLS MD", as velocity rescale thermostat was used in both cases. Reduction in computation time decreases with increase in bin width size, from (0.75 - 1.25) to (0.75 - 2.2), with minimum values achieved for 50% reduction. Minimum potential energy of the Argon system was achieved in "TLS MD" of -3009.8 when compared to conventional MD simulation (-2996.5), indicating more optimized phase space configurations were achieved in "TLS MD" than in conventional MD simulation.

Comparison of potential energy (PE) for "TLS MD" with conventional MD was shown in Figure 5, with bin width of (0.75 - 1.25) and (0.75 - 1.8). "TLS MD" phase space has its unique movement over MD time steps when compared to conventional MD simulation. It was observed that bin width of (0.75 - 1.8) have closer potential energies with conventional MD simulation. Average PE for bin width (0.75 - 1.25) is higher than conventional MD simulation, as shown in Figure 5 and Table 3.

Radial distribution function (RDF) for conventional MD and "TLS MD" were considered to study the atomic distributions after reaching equilibrium. Figure 6 shows RDF plots for conventional MD and "TLS MD" with bin width of (0.75 - 1.25) and (0.75 - 1.8). It was observed that "TLS MD" shows good agreement with conventional MD in first, second, and third solvation shells. For "TLS MD" with bin width of (0.75 - 1.25) have higher peak in first solvation shell which reflects in average potential energy and in PE plot, in Table 3 and Figure 5, respectively. Same happens at the minimum of first solvation shell for "TLS MD" with bin width of (0.75 - 1.25) when compared with conventional MD simulation, as shown in Figure 6. Minimum deviations were found for second, third, and fourth solvation shells between "TLS MD" and conventional MD simulations. Bin width less than (0.75 - 1.25) gives RDF that deviate from RDF of conventional MD simulations, giving different atomic configurations and analysis of such configurations is beyond the scope of this study.

By comparing average energies, average temperature (Table 3), PE plot (Figure 5), and RDF plot (Figure 6), "TLS MD" shows good agreement with conventional MD simulations with reduction in computation time from 50% - 60%. Correlation and extrapolation of force/potential values of local region to global region was not considered in this study.

4.3. Curve Fitting Algorithm for Acceleration Prediction in MD Simulations

Another approach in reducing computation time in MD simulations is predicting/extrapolating properties of future time steps, from previous time steps using function approximations like curve fitting techniques [20], as shown in Figure 7. It was observed that force/acceleration time profiles on a particular atom in MD simulations shows smooth curves with joining splines with quadratic/cubic in nature as shown in Figure 8. Figure 8 shows acceleration time profile in x-axis (a_x) for an individual Argon atom in conventional MD simulation for Argon system. It shows that force/acceleration at intermediate steps can be predicted using previous time steps using quadratic/cubic equations. This technique of apply curve fitting techniques show reduction in overall computation time, as intermediate steps are predicted with free of cost in terms of computational load. In this study, acceleration for each Argon atom was predicted at every 4^th time step from the previous 3 time steps using quadratic function approximation. For every 4^th MD time steps, one MD step is predicting instead of calculating force/acceleration in conventional way which reduces computational load in MD simulation. It was expected to reduce computation time of 25%, since for every 4 time steps, one time step is free of cost. Higher order function approximation [20] for force/acceleration time profiles can be considered to predict two or more time steps ahead using previous time steps, that reduces computation time more than 25% and this approach is beyond the scope of this study.

In this study, acceleration was calculated at every 4^th time step by assuming quadratic equation (Equation 1) for acceleration over time. To find the constant of a, b, and c, previous three time steps are considered as -1, 0, 1, time values (Equations 2, 3, and 4), for simplicity. Fourth time step will be calculated as per Equation 5. Constants of a, b, and c were calculated from Equations 6, 7, and 8, using previous acceleration values (a_c1, a_c2, and a_c3). Equation 9 shows acceleration value at every 4^th time step (a_c4) in terms of acceleration values from previous time steps.

$$a_c=a{t}^2+bt+c$$

(1)

$$\mathrm{at}\mathrm{t}=-1,a_{c1}=a-b+c$$

(2)

$$\mathrm{at}\mathrm{t}=0,a_{c2}=c$$

(3)

$$\mathrm{at}\mathrm{t}=1,a_{c3}=a+b+c$$

(4)

$$\mathrm{at}\mathrm{t}=2,a_{c(4,now)}=4a+2b+c$$

(5)

from equations (2), (3), and (4), constants a, b, c were obtained as follows:

$$a={\displaystyle \frac{a_{c1}+a_{c3}-2a_{c2}}{2}}$$

(6)

$$b={\displaystyle \frac{a_{c3}-a_{c2}}{2}}$$

(7)

$$c=a_{c2}$$

(8)

Substituting the constants (a, b, and c) in Equation (5) gives acceleration for the present time step (a_c4,now) as follows:

$$a_{c4,now}=3a_{c3}-3a_{c2}+a_{c1}$$

(9)

Table 4 shows reduction in computation time in apply curve fitting algorithm in MD simulation (called as "Curve Fitting MD") with conventional MD simulation. As predicted, 24% reduction in computation time was achieved with potential energies, kinetic energies and temperature average values shows similar to conventional MD simulations (Table 4).

Figure 9 compares potential energy (PE*) profile for "Curve Fitting MD" and conventional MD for the first 2000 time steps of Argon atoms. During initial time steps, "Curve Fitting MD" follows same as conventional MD profile and later departs by have its own phase space configurations over time. Both MD simulations have average potential values as same with slight variations having their own phase space configurations over time.

RDF plots comparison shows similar atomic distribution for first, second, and third solvation shells for "Curve fitting MD" and conventional MD simulations, as shown in Figure 10.

4.4. Combining TLS & Curve Fitting Algorithm in MD Simulations

Further reduction in computation time in MD simulations can be achieved by combining "TLS MD" and "Curve Fitting MD", where force/potential can be calculated for the first 3 time steps using "TLS MD" and the 4^th time step will be predicted using "Curve Fitting MD" technique. Another 10% of overall reduction in computation time will be expected, as 25% reduction in computation time of "TLS MD" will be achieved due to "Curve Fitting MD". Every 4^th time step of acceleration value is predicted using previous 3 time steps whose acceleration was calculated using local sampling. Table 5 shows reduction in computation time for "TLS & Curve Fitting MD" when compared with conventional MD simulations along with calculated potential energy, kinetic energy, and temperature average values. Maximum percentage of reduction in achieved by 67% for bin width of (0.75 - 1.25) with minimum value of 61% for bin width of (0.75 - 2.2). Potential energies and kinetic energies shows good agreement with conventional MD simulations as shown in Table 5.

RDF plots comparison for "TLS & Curve Fitting MD" and conventional MD simulation were analyzed for bin widths of (0.75 - 1.25) and (0.75 - 1.8), as shown in Figure 11. RDF plot shows similar behavior as "TLS MD" (Figure 6) and inclusion of "Curve Fitting MD" in "TLS MD" has lesser impact in distribution of atomic configurations after reaching equilibrium.

5. Conclusion

Theory of local sampling (TLS) in molecular simulations was introduced in this study with reported reduction in computation time in MC and in MD simulations. Reduction in computation time up to 60% was achieved in Theory of Local Sampling Molecular Dynamics (TLS MD) and in "TLS MC" simulations up to 37% reduction was achieved. "Curve fitting" techniques was implemented in MD simulations ("Curve Fitting MD") for predicting force/acceleration acting on atom at every 4^th step using previous 3 time steps, with achieved reduction in computation time of 24%. Further, combining "TLS MD" and "Curve Fitting MD" reduces more computational time with maximum reduction of 67% was achieved when compared to conventional MD simulation. TLS can be applied to other systems with principle of calculating local sample and scaling back to global sampling values for systems like ab initio molecular simulations, quantum calculations, Computational Fluid Dynamics (CFD), and dynamics of stock/equity market. "TLS MD" using random local sampling, modeling of error estimation with conventional sampling for MD simulations will be another approach to reduce computational load in MD simulations. Neural networks and artificial intelligence (AI) are some of the tools that can be integrated in TLS in molecular simulation that will reduce computation time to many folds.