Renormalization group theory of molecular dynamics

Large scale computation by molecular dynamics (MD) method is often challenging or even impractical due to its computational cost, in spite of its wide applications in a variety of fields. Although the recent advancement in parallel computing and introduction of coarse-graining methods have enabled large scale calculations, macroscopic analyses are still not realizable. Here, we present renormalized molecular dynamics (RMD), a renormalization group of MD in thermal equilibrium derived by using the Migdal–Kadanoff approximation. The RMD method improves the computational efficiency drastically while retaining the advantage of MD. The computational efficiency is improved by a factor of 2n(D+1)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$2^{n(D+1)}$$\end{document} over conventional MD where D is the spatial dimension and n is the number of applied renormalization transforms. We verify RMD by conducting two simulations; melting of an aluminum slab and collision of aluminum spheres. Both problems show that the expectation values of physical quantities are in good agreement after the renormalization, whereas the consumption time is reduced as expected. To observe behavior of RMD near the critical point, the critical exponent of the Lennard-Jones potential is extracted by calculating specific heat on the mesoscale. The critical exponent is obtained as ν=0.63±0.01\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\nu =0.63\pm 0.01$$\end{document}. In addition, the renormalization group of dissipative particle dynamics (DPD) is derived. Renormalized DPD is equivalent to RMD in isothermal systems under the condition such that Deborah number De≪1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$De\ll 1$$\end{document}.

There are a vast number of modern applications where phase transition or interface plays an important role such as laser annealing and 3D printing. It has long been desired to perform computational experiments of such macroscopic, multiphysics problems, which are often challenging for continuum approaches. There has also been a strong interest in a method capable of seamlessly simulating from microscale, mesoscale to even a continuum regime. The only method that satisfies these demands in principle is molecular dynamics (MD) 1,2 , and it has been adopted in a wide variety of fields such as nanostructures 3 and biochemistry 4,5 .
A major shortcoming of MD is that it requires a vast amount of computational resources especially for solving large systems. Recently, advancement in hardware, parallel computing methods as well as graphical processing units (GP-GPUs) have drastically increased the simulation scale achievable [6][7][8][9] , but a majority of MD studies are still limited to the orders of micrometers and nanoseconds, and macroscopic analyses are often even impractical. Several solutions have been developed to overcome this limitation. In the fields of biochemistry and biophysics, common techniques are coarse-graining [10][11][12] and enhanced sampling methods [13][14][15] . Another approach, often adopted for nanofluids and nanostructures, is hybrid methods of atomic-continuum domains 16,17 .
Here we propose a new approach to expand the applicability of conventional MD. P. W. Anderson suggested that calculations in arbitrary scale are possible by applying the renormalization group (RNG) theory to transport phenomenon 18 [pp. 212]. If the Hamiltonian describing MD is at a fixed point, following the RNG transformation, a calculation in arbitrary scale can be conducted with the same number of particles. This idea coarse-grains the entire domain uniformly so that a drastic improvement on the computational efficiency may be achieved in mesoscale or even macroscopic systems, while retaining the aforementioned advantages of MD. The RNG of elastic bodies has been constructed for coarse-grained MD [19][20][21] [Sect. 13.5]. Faccioli et. al adopted the idea of RNG to decouple the short-time fluctuations and the long-term dynamics of the molecular system described by the Langevin equation 22,23 (i.e. dissipative dynamics), but the same idea cannot be applied to MD since energy of a particle is not conserved. To the best of our knowledge, the RNG for MD has not yet been constructed.
In the current work, we derive the RNG by using the Migdal-Kadanoff approximation 21,24 [Sect. 16.6.4]. The obtained RNG for MD, or renormalized molecular dynamics (RMD), is verified by two test problems; melting of an aluminum slab and collision of aluminum spheres. We then discuss the computational efficiency of RMD over conventional MD and the similarity of physical phenomena after the renormalization. Furthermore, to observe behavior of RMD near the critical point, we calculate specific heat at constant volume of the Lennard-Jones fluid Coarse-graining of an atomic chain. We start coarse-graining from the Hamiltonian of an atomic chain of length L. A distribution function of the Canonical ensemble with N particles is given as: where h is the Planck constant, β = 1/(k B T) , T is temperature and k B is the Boltzmann constant. The coarsegraining of the atomic chain is obtained by removing a degree of freedom of a particle j locating at midpoint of particles i and k. The potential on the particle j formed by the nearest neighbors can be written as: Taylor series expansion on r i −r k 2 gives: Here, x j and φ (2l) are defined as: We define z 2l as a root of φ (2l) (r) . For the Morse and LJ potentials, φ (2l) (r = z 2l ) are given as: Therefore, having φ (2l) (z 2l ) = 0 yields: Next step is the integration for variables p j and r j . Using Eq. (4) gives: for LJ potential 2 increases, a sign of φ (2l) (r i,k /2) changes from positive to negative at z 2l , then φ (2l) (r i,k /2) , after taking minimum values, asymptotically approaches to zero. In the range of φ (2l) (r i,k /2) ≤ 0 , the integral Eq. (10) diverges at a limit of �L/(2σ ) → ∞ . It is required to have v f ≤ L 3 (i.e. the collective entropy), so φ (2l) (r i,k /2) should be adjusted to zero if φ (2l) (r i,k /2) ≤ 0 . By this manipulation, a leading term in the range of φ (2l) ≤ 0 would become a term of ∼ φ (2(l+1)) x 2(l+1) j .
As an example, we discuss an approximation when the summation is cut off at l = 2 . Eq. (10) can be written as: We have φ (2) > 0 in the range of 0 < |r ij | 2 < z 2 . Ignoring a term of x 4 j , we obtain: Similarly, we have φ (2) < 0 and φ (4) > 0 in the range of z 2 < |r ij | 2 < z 4 . Thus, φ (2) is adjusted to zero and we have: By the approximation described above, at an arbitrary r i,j , the integral can be rewritten as: Here, z 0 = 0 and always z 2(l−1) < z 2l . Large x j does not contribute to the integration because φ (2l) ≤ 0 does not exist. Therefore, by extending the integration range as �L/(2σ ) → ∞ , Eq. (11) is reduced to the gamma function. From definition of the gamma function, Ŵ(l −1 ) = l ∞ 0 e −x l dx , the integration becomes: By substituting Eqs. (9) and (13) into Eq. (3), the distribution function of the atomic chain after the coarsegraining can be obtained as: where F ′ is the Helmholtz free energy after the coarse-graining. Hereafter, an apostrophe indicates the coarsegrained property. The second term of the right hand side of Eq. (16) is given by the Stirring approximation: The entropy S ′ and the internal energy E ′ of the atomic chain after the coarsegraining are: www.nature.com/scientificreports/ where s(r) is a step function derived from s(r) = T ∂ ∂T log v f such that: By the decimation N → N/2 , the coarse-grained Hamiltonian of the atomic chain is represented as: The last term of the right hand side of Eq. (18), 3 2 k B T , is removed from Eq. (20) since it disappears by the normalization.
Coarse-graining of the D-dimensional Hamiltonian. In the previous section, we derived the onedimensional coarse-graining of the atomic chain. Next, the result is extended to the D-dimensional Hamiltonian by the Migdal-Kadanoff approximation 21,24 [Sect. 16.6.4].
First, we discuss the coarse-graining of a two-dimensional Hamiltonian. We assume that the time average of the atomic positions forms a simple (i.e. square, oblique, rectangular and hexagonal) lattice. The lattice with the time averaged positions of atoms maintains the periodicity of a simple lattice, so that the Migdal-Kadanoff approximation can be applied to atomic systems. For convenience, we adopted the square lattice, which is one of the simple lattice, in the following discussion including figures. Nevertheless, the coarse-graining procedure on this section does not lose its generality.
In the square lattice configuration, a particle at point r i,j has four nearest neighbors at r i+1,j , r i−1,j , r i,j+1 , r i,j−1 and four next nearest neighbors at r i+1,j−1 , r i−1,j−1 , r i−1,j+1 , r i+1,j+1 . The particle at r i,j can move inside an area of L 2 . The atomic arrangement is shown on Fig. 1a. A basic idea of using the Migdal-Kadanoff approximation is to remove an interaction and distribute it into adjacent interactions, so that atomic chains appear and one-dimensional coarse-graining is applicable. For example, let us focus on φ(r i+1,j − r i,j ) , one of the nearest neighbor interactions on the particle at r i,j . By the Migdal-Kadanoff approximation, the interaction (17) www.nature.com/scientificreports/ is split into half (each has strength of 1/2ǫ ) and distributed to adjacent interactions φ(r i+1,j+1 − r i,j+1 ) and φ(r i+1,j−1 − r i,j−1 ) . This process is applied to all four nearest neighbor interactions on the particle at r i,j . The same approximation is repeated to adjacent blocks, which contain particles at r i+2,j , r i−2,j , r i,j+2 and r i,j−2 . Therefore, By these manipulations, the Hamiltonian can be rewritten as: Schematics of the square lattice and the Migdal-Kadanoff approximation is shown on Fig. 1a, b. An integration about p i,j is trivial, and the particles at r i,j−1 , r i−1,j , r i,j+1 , r i+1,j can be considered as particles on the atomic chain. Applying the result of the coarse-graining of the atomic chain Eq. (20) gives: Finally, by extending the same procedure to the D-dimension, the coarse-graining of the D-dimensional Hamiltonian is obtained as: Derivation of the renormalization transforms. A list of coupling constants, denoted as K , is given as: There exist two possible renormalization transforms such that an action −βH is at a fixed point. One of such transforms introduces scaling of the space, L → L ′ = 2 −1 L . By this transform, the coupling constants are scaled as: Another transform retains the spatial size L invariant, so K is scaled as: Mass density ρ = Nm/L D is at a fixed point for both transforms. However, R a results in the Young modulus Y and the speed of sound c s relevant, since Y ∼ ǫ/(σ 2 r o ) and c s ∼ ǫr 2 o /(mσ 2 ) ( ∼ √ k B T/m for gas) 25 [Sect. 3.2], whereas R b yields both Y and c s invariant (i.e. at a fixed point). Therefore, we adopted R b for the current work. After the n renormalization transforms, the coupling constants K n and the total potential n can be written as: www.nature.com/scientificreports/ It should be noted that the 2l-th spatial derivative of s in the interaction, s (2l) , disappears when calculating v f . The D-dimensional Hamiltonian after the n renormalization transforms H n is given as: where = 2 n is a scaling factor and r i,k = |r i − r k | . < k, i > denotes taking a summation of k on the nearest neighbor particles around the particle i. The second term of the right hand side of Eq. (30) has a contribution of (k Finally, the equation of motion of the i-th particle is obtained by substituting the renormalized Hamiltonian Eqs. (29) and (30) Table 1.

Renormalization of Langevin dynamics.
We discuss an application of the RNG to dissipative dynamics. Renormalization group of the dissipative dynamics is derived. The damping coefficient ξ can be expressed as ξ ∼ ηr o , so it is scaled as ξ → 2 ξ by the renormalization. Hence, the RNG for the dissipative dynamics system R d (K ) can be constructed by adding the transform ξ → 2 ξ into the renormalization transforms R b (K ) as: Renormalization of dissipative particle dynamics (DPD) 26 is derived as follows:  www.nature.com/scientificreports/ Here, e ij = r i,j /r i,j ,w D is a weighting function and ζ ij (t) is δ-correlated Gaussian noise. It can be easily verified that it follows the fluctuation dissipation theorem after the renormalization. In DPD, φ(r i,j ) of Eq. (30) is defined as an interaction potential acting between particles in the dissipative dynamics system. To have the RNG constructed, it is required that the interaction potential as an identity element does not cause divergence on calculation of a free volume v f (Eq. (10)). For example, the Morse and LJ type potentials are renormalizable. A power-law repulsive potential φ = ǫ(σ/(r i,j − r o )) m is also renormalizable because φ (2l) (r i,j ) > 0 and s(r i,j ) ≈ 1 2 . However, the attractive potential φ = −ǫ(σ/(r i,j − r o )) m can not be renormalized.
Since the Deborah number De of DPD is scaled as De ′ = De , De is at a fixed point only in the limitation of De → 0 or De → ∞ . Therefore, De of renormalized DPD cannot be the same as that of the real system. On the other hand, renormalized DPD is equivalent to RMD in isothermal systems, provided that De ≪ 1 and φ is an interatomic potential.

Results
Melting of an aluminum slab. In this section, we present verifications of RMD method by conducting two test problems. For the first problem, we considered melting of aluminum by simulating a slab of aluminum with heat added from its bottom. On Fig. 2a, temperature profile of the original system ( = 2 0 ) is plotted along E in , total energy input to the system. A presence of a plateau between two linear lines indicates the latent heat of melting. Three linear curves were fitted to the linear lines and the plateau of the profile by the least squares method to determine the melting temperature and the latent heat. These values were found to be 833.0 K and 467.3 kJ/kg, respectively. A horizontal dashed line on Fig. 2a indicates the melting temperature, and vertical dashed lines are endpoints of the plateau.
Simulations were repeated in renormalized systems at scaling factors of = 2 1 and 2 2 . The temperature profiles on Fig. 2b, c show the similar plateaus at the same temperature and energy input as the case of = 2 0 , except that larger fluctuations are present due to the fewer particles in the renormalized systems. The similarity between the original and the renormalized systems, such as crystal structures and progress of melting surfaces, can be visually observed on snapshots shown on Fig. 3 while the particles are drastically coarsened.

Collision of aluminum spheres.
For the next verification, we carried out a simulation of a collision of two aluminum spheres. The simulations were conducted with rather large scaling factors, = 2 22 , 2 23 and 2 24 to investigate applicability of RMD to macroscopic systems. An aluminum sphere with radius of 34.3 mm and velocity of 50.0 m/s was collided to another stationary sphere with identical size and temperature, and behavior of the collision was monitored through sphere velocity and system kinetic energy. Figure 4a is a plot of the velocity of the moving sphere versus time after contact, and Fig. 4b shows profiles of the system kinetic energy. The test conditions and the simulated contact durations are tabulated on Table 2. For reference, a theoretical contact duration of elastic spheres 27 [Sect. [1][2][3][4][5][6][7][8][9] is calculated to be 93.24 × 10 −6 s.
The kinetic energy and velocity profiles as well as the contact durations of three cases are agreed well whilst the number of particles are drastically reduced as the system is renormalized. This indicates that the similarity of the phenomena is retained after the renormalization in the same way as the melting of the slab, therefore suggests that the macroscopic MD calculation with only a few hundred thousands particles is achievable. It should be also noted that the total energy of each system is conserved throughout the calculations.  www.nature.com/scientificreports/ N ′ ∼ N/ D , t c is transformed as t ′ c = t c / D+1 . Therefore, RMD achieves a higher efficiency by a factor of D+1 over conventional MD. Figure 5 shows the total consumption time versus the scaling factors. The plots show a clear indication that the computational efficiency improves accordingly as the system is renormalized.   www.nature.com/scientificreports/ It should be noted, however, that a term s(r) resulted from the renormalization is neglected in the current test problems. If the term is not negligible, such as at a high temperature condition k B T ≫ ǫ or in the lower dimensions D < 3 , s(r) needs to be included to calculations so that it would result in additional computational cost. In MD calculations, the largest portion of the computational cost is calculations of forces. A simple estimation suggests that, in the worst case scenario, the practical computational efficiency would be approximated as a factor of 2 n(D+1)−1 , because only the calculations of forces are doubled. DPD contains a time constant τ = m/ξ . The time step of DPD is also restricted by �t ′ ≪ τ ′ and τ ′ = τ . This implies that the time step t is scaled as t ′ = t , therefore it yields the same computational efficiency as RMD.

Critical behavior in the liquid-vapor region on mesoscale.
We investigate the applicability of the renormalized Hamiltonian at a fixed point to phase transition near the critical point. To do so, we observe the behavior of specific heat at constant volume C V in the near-critical region by the finite size scaling 28,29 [Sect. V.2.2]. When scaling of the distribution function by the renormalization is invariant near the critical point, following the finite size scaling theory, C V and T can be scaled as: where T c is critical temperature, ν is critical exponent and L is domain size (side length of a cubic domain). Although �(x) is a function of a variable x, a form of the function �(x) of x is not explicitly known. However, �(x) obtained from different system sizes can be collapsed onto a single curve if values of T c and ν are valid.  www.nature.com/scientificreports/ We calculated C ′ V of argon near the critical temperature in three different mesoscale sizes, L = 11σ ′ , 16.5σ ′ and 22σ ′ by RMD. The scaling factor is set to = 2 8 . �(x) and x are calculated by Eqs. (35) and (36) from the resulting C ′ V and its temperature dependence. We explored several values of the critical exponent ν for the critical temperature, T ′ c / 3 = 157.2K [30][31][32] and located that �(x) of the different mesoscale sizes approximately collapse onto a single curve at ν = 0.63 ± 0.04 as shown Fig. 6. Figure 6 is a plot of �(x)/L 3−2/ν max versus x/L 1/ν max at ν = 0.59, 0.63 and 0.67, where L max refers to the largest one of three sizes, i.e. L max = 22σ ′ . ν = 0.62 and 0.64 were also explored, and the curves of three sizes were not clearly distinguishable (the plots are omitted). Therefore, we found that the curves of three sizes overlap at ν = 0.63 ± 0.01 . This is close to the value of the existing work ν = 0.630 obtained by the conventional MD calculation of the LJ fluid 33 and by the study of the same universality class (i.e. three dimensional Ising model) 34,35 . More precise extraction of the critical exponent is not a scope of the current work, however, our results demonstrated that RMD is capable of reproducing the critical behavior. It should be mentioned that the critical slowing down is not improved due to the reduction in the number of particles by the renormalization.

Discussion
In this study, we have constructed the RNG of MD by using the Migdal-Kadanoff approximation. The total mass and energy are conserved in the renormalized systems, and the obtained scale transformation rules show that especially Young modulus and speed of sound remain invariant. We conducted two simulations to validate the RMD method. Both problems have demonstrated that the expected values of physical quantities are in good agreement after the renormalization, and the computational efficiency is improved as expected by 2 n(D+1) over conventional MD. Furthermore, to observe behavior of RMD near the critical region, temperature dependence of C V of the LJ fluids was calculated near the critical temperature. The results showed that the critical exponent ν can be extracted by the finite size scaling theory in RMD, and it was confirmed that the obtained ν was agreed with previous work regarding the conventional MD calculation of the LJ fluid and the study of the same universality class (i.e. three dimensional Ising model). However, RMD does not improve the critical slowing down due to the reduction in the number of particles by the renormalization.
In addition, we explored the application of the RNG to DPD, and derived the renormalized DPD equations. To have the RNG constructed for DPD, it is required that an interaction potential as an identity element does not cause divergence on calculation of a free volume v f (Eq. (10)). Renormalized DPD is equivalent to RMD in isothermal systems, provided that an interaction φ is an interatomic potential and the Deborah number De ≪ 1 . We also showed that the computational efficiency of renormalized DPD is similar to that of RMD.
Our work proposes a new approach for MD, which provides a drastic improvement on the computational efficiency for large scale systems while retaining the advantages of conventional MD. For future work, we will also investigate and discuss how the term of s(r) would contribute to the RMD calculations in the conditions such that the term is not negligible. Starting from the current study, which coarse-grains the whole domain uniformly, we would extend RMD to more practical method, such as the derivation of the local group which should be applicable in a similar manner to irregular meshes.

Methods
All the simulations of the current work are carried out by our in-house code written in C. Practically, computational algorithms of RMD is identical to conventional MD, which solve the equations of motion Eqs. (31) and (32), except that the temperature T, the potential parameters σ , r o , ǫ , and the mass m need to be scaled accordingly to the renormalization transforms K n = (2 −nD β, 2 nD ǫ, 2 n r o , 2 n σ , 2 nD m) (Eq. (27)). Second term of the right hand side of Eq. (31) is neglected in the current work. The velocity Verlet scheme is used for the time integration.
Test conditions of the melting phenomenon. The Morse potential is adopted for the potential of alu- is (L x , L y , L z ) = (3.883 × 10 −6 m, 0.971 × 10 −6 m, 0.971 × 10 −6 m) . Note that L x is set to be twice the initial thickness of the slab so that it provides some space for the slab to expand and evaporate on the surface. Mirror boundary conditions are applied in the x-direction to prevent atoms from leaving the domain, and periodic boundary conditions are applied to all other directions.
Four layers of the atoms from the bottom of the domain are treated as a wall for heat transfer (shown as red particles on Fig. 3). The number of wall atoms is N ′ p,wall 3 = 4608 . The heat is added to the slab by increasing the kinetic energy of the atoms of the wall via rescaling the momentum. Only momentum of each particle in the wall is rescaled as p i,rescaled = p i 1 + 2m ′ �e ′ i /p 2 i , where p i is the current momentum of particle i, p i,rescaled is the momentum after rescaling and e ′ i is energy added to the particle at every 4.0 × 10 −13 s. This interval is www.nature.com/scientificreports/ commonly used for all = 2 0 , 2 1 and 2 2 . The energy E ′ = N ′ p,wall e ′ i = 7.68 × 10 − 20 J is the total energy added to the wall. The wall atoms interact with other atoms via the Morse potential, and spring-like forces also act on these atoms. This allow the wall atoms to vibrate but not freely move away from their original positions. The spring-like force is given as F ′ i,spring = −k ′ δr i where k ′ / = 20.0 kg/s 2 is a spring constant and δr i is a displacement of an atom i relative to its original position.
Temperature is initially set as T ′ / 3 = 700 K, and is relaxed at the same temperature for the first 1.0 × 10 −6 s by the Berendsen thermostat. The thermostat is applied at every 4.0 × 10 −13 s, and the time constant of the thermostat is set to 1.0 × 10 −5 s. This interval and duration are commonly used for all = 2 0 , 2 1 and 2 2 . Heat addition is then started and the temperature of the slab is sampled at every 2,000 time steps. The wall particles are not included to the temperature calculation. Note that the total energy added to the system is invariant, E ′ = E , but since the number of the wall atoms are reduced as N ′ p,wall 3 = N p,wall , the energy added to each atom is scaled as �e ′ i / 3 = �e i .
Test conditions of the collision of two sphere. The potential and parameters are identical to the first problem, except that the potential is constant φ(r o ) if r > r o for interactions between particles which belong to different spheres. This provides only repulsive forces between spheres to prevent spheres from sticking together upon contact. Particles are initially placed in fcc configuration inside a spherical region of a given radius R = 3.863 × 10 −2 m. The radius is selected so that particles forming sphere surface are located on the exactly same distance from the center of the sphere at all three cases. Temperature is initially set to be T ′ / 3 = 100 K, and gradually cooled down to 0 K by the weak Berendsen thermostat. The time constant of the thermostat is set to 20.0 t ′ . The thermostat is applied every 2.19 × 10 −6 s, and its duration is 1.47 × 10 −3 s. This interval and duration are commonly used for all = 2 22 , 2 23 and 2 24 .
The obtained sphere is duplicated and placed so that separation distance of sphere centers is 7.49 × 10 −2 m. The spheres are left stationary (without thermostat) for another 3.66 × 10 −4 s for relaxation, then velocities of the particles of one sphere are set to be 50.0 m/s toward the other sphere. The velocity of the sphere is calculated by averaging the velocities of the particles which belong to the sphere.
Test conditions of the specific heat calculation. The Lennard-Jones potential is adopted for the poten- The temperature is initially set at T ′ / 3 = 360 K, then the Nose-Hoover thermostat is turned on at t = 5.12 × 10 −9 s to bring the system to the target temperature. The total energy is sampled from t = 5.12 × 10 −7 to 1.02 × 10 −6 s at every 100 time steps. C ′ V is calculated from the fluctuation of the total energy E ′ under the NVT condition as follows: where indicates the ensemble average. Note that C V is scaled as C ′ V = C V / 3 as shown on Table 1a. Fluctuations of the energy near the critical point is highly sensitive to initial configurations of the system. To obtain statistically more meaningful results, we averaged over the results of five runs with different initial configurations for each test condition.

Data availability
The data that support the findings of this study are available from the corresponding author upon request.

Code availability
An in-house MD computer code was used for all the calculations of this study. The code is available from the corresponding author upon request.