Molecular Dynamics for Two-Body Potential from Unobserved Gaussian Regression

A molecular dynamic simulation requires the local potential energy per atom because it provides the force applied to each atom. On the other hand, an electronic structure computation provides only the global potential energy of an atomic system. We propose a stochastic fitting for deducing a potential from a total energy which can be incorporated inside a molecular dynamic program. The objective of that fitting process is twofold. First, the total energy can be reproduced with a sufficient accuracy. In addition, we want to separate the potential energy into atomic contributions. The only inputs in the stochastic regression are the total energies of the atomic systems. As for the applications, we examine the performance of the method with the help of a mathematical model. Afterward, we use it for quantum applications in which we compare the direct method and the unobserved fitting with respect to the energy conservations. For a molecular dynamic application, we examine the atomic cluster formations during freezing when the proposed potential is used. Although the method is generally applicable, we restrict in this paper to total energies which are reproduced from two-body potentials for the molecular dynamic simulation.


Introduction
Multiscale computation (Tadmor & Miller, 2011) is the sequence of simulations from nanoscale quantum mechanics till macroscopic simulation.This sequence usually consists of electronic structure, molecular dynamics (MD), mesoscopic simulation, quasi-continuum, till FEM/CAD (Finite Element Method/Computer Aided Design).The most important step in standard MD packages (Allen & Tildesley, 1989;Ciesla et al., 2007;Kress et al., 1998) is the determination of the local forces which are derived from the potential energy.There are numerous empirical potentials including Lennard-Jones, Finnis-Sinclair, ReaxFF(Reactive Force Field), REBO(Reactive Empirical Bond Order) and BOP(Bond Order Potential) among others.In Figure 1(b) we illustrate a molecular dynamic simulation of a graphene fracture using REBO potentials where the colors indicate the potential energy per atom.In this document, we would like to consider the part of multiscale computation between electronic structure such as DFT and molecular dynamics which is graphically illustrated in Figure 1(a).The main problem in bridging that gap is the separation of atomic contributions of the potential energy.In fact, a quantum packet provides only the total energy of a given atomic configuration whereas a molecular dynamic package requires the atomic contribution in order that the force applied to each atom can be computed at every molecular dynamic step.We are interested in molecular dynamics where the potential is obtained from a prefiting process which consists in splitting the potential energy into atomic local energies.We propose a stochastic method whose only observed inputs are the total energy while the objective is to reconstruct the local energy per atom.The method presented here is based upon machine learning which uses the concept of learning on the fly.It means in general that when an atom or a configuration passes in one state, a certain computation is performed.Upon passing on that state for the next times, the computation is not repeated any more but experiences from previous computation is used.The machine learning process generates a structure (Seeger et al., 2003) which enables the efficient and fast access of such a methodology.Our motivation is to generate a system which is both accurate and fairly inexpensive to evaluate.The proposed method is based on stochastic method where the expectancy and the variance of the prediction need to be computed.One utilizes the Matern correlation function for the assembly of the covariance matrices.The stochastic posterior depends only on the input observed total energies.As a consequence, it can be computed once for all and it is looked up during future computations.For an MD simulation, what is really required is not the energy but the force derived from that energy.Hence, the derivative of the prediction will be computed analytically.That is important because using Finite Difference is known to diverge the energy conservation in MD.During the simulation, one needs to solve a nonlinear optimization pertaining to the logarithm marginal likelihood.We present a method for avoiding the singularity at the diagonal of the covariance matrix.
Our method can be applied to all sorts of quantum information but we mainly focus on the atomic energy in this paper.
The practical implementation is divided into two applications.The first one is the investigation of the method by using some mathematical system.In that case, we consider some exact solution and we reconstruct the internal local energy function by using our theoretical method.The purpose of that test is to examine the performance of the method by comparing the exact solution with the solution provided by the fitting.The local energies are expressed in term of the Coulomb matrix.As for the second application which is based on quantum systems, this paper presents only a preliminary work because evaluating DFT is a very time-consuming process.As a consequence, we compare an MD using Lennard-Jones and an MD where the potential is obtained from unobserved Gaussian regression.In order to gauge the simulation quality, we investigate the radial distribution function (RDF) for various values of the temperature by using the potential which was obtained from the unobserved fitting.The general tendency of the RDF curve is to start from zero values followed by a sudden peak which in turn is followed by a slow approximation of the unit value.In general, the height of the peak depends on the temperature.For very cold temperature, the peak is remarkably high while it is barely observable when the temperature is very warm.To validate the use of the atomic separation approach, we investigate the results of the MD simulation at freezing temperatures.The formation of clusters of atoms which features the phenomenon of freezing temperatures as described in (Chen, 2013;Yamaguchi & Maruyama, 1998) will be practically reconstructed.The adjustment of the temperature value is regulated by using the method of kinetic energy scaling.By using the unobserved Gaussian regression, we achieve total energy conservation as one observes the energy values at many MD steps.
In addition, we will compare the energy values by using the direct method and the proposed atomic separation method.

Potential Energy w.r.t. Nuclei Coordinates
Consider a system composed of n atoms R = {r i } n i=1 and m electrons X = {x i } m i=1 .The main problem in electronic structure computation involves the general Hamiltonian (Tadmor & Miller, 2011;Hill et al., 2005;Saad et al., 2010) operator where is the Planck constant divided by 2π, q e is the electronic charge, m e and M k represent the electron mass and the mass of the k-th nucleus while the charge of the j-th nucleus is Z j = |q e |n j in which n j is the atomic number.
The wave function Ψ is the eigenfunction corresponding to the smallest eigenvalue of the eigenvalue problem The kinetic energies of the nucleus and electron are The third summations of (1) provide the nucleon-electron Coulomb attraction energy.The two last summations are the inter-nuclei and inter-electron repulsion energies.We shall define the atom-electron interaction, the internucleon operator and the inter-electron interaction operators V AE , V AA and V EE as follows where Thus, one has We use the Born-Oppenheimer or adiabatic approximation assumption stating that the mass and the volume of the atoms are very large in comparison to those of the electrons.In fact, an electron is 2000 times lighter than any nucleus such that the atoms move comparatively slower than the electrons.Thus, electrons are supposed to follow the movements of the atoms in the sense that As a consequence, we treat the time-independent Hamiltonian operator with respect to the a set of nuclei r i which are supposed to be stationary.That is to say, for each given R = (r 1 , ..., r n ), the initial electronic structure is reduced to the expression Thus, the eigenproblem corresponding to the above Born-Oppenheimer Hamiltonian operator becomes in which the ground state energy E R and the reduced wave function ψ R depend on the given nuclear coordinates R = (r 1 , ..., r n ).Due to the independence of Φ on X = {x i } in (9), one deduces from (11) that By using (2) and ( 9) in the last equality, one obtains By assuming the independence of ψ R (X) on the atomic kinetic term such as Thus, one obtains the Schrödinger equation with respect to the nuclear coordinates only in which For the molecular dynamic applications, one uses V pot which is the effective inter-atomic potential or the potential energy surface.With the help of the quantum assumptions above, that last equation takes only the nuclei coordinates r 1 , ..., r n of the atoms as variables.Since it is generally very expensive to evaluate the function one uses a simplification or an approximation of E pot .It is impossible to express that function analytically for large MD systems which consist of multi-millions atoms.Efficient empirical approximations of (20) for extremely large systems include: Finnis-Sinclair, EAM (Embedded Atom Method) and BOP (Bond Order Potential).Exact solutions for the above equation are only known for very few special cases.Hence, numerical methods must be used for the general cases.The most direct method is to consider the domain of computation Ω ⊂ R 3 which is supposed to be sufficiently larger than the nuclei cloud {r j } n j=1 such that the above operator has neglecting influence beyond it.Every electron x i ∈ Ω admits a spin σ i which is up (σ i := +1/2) or down (σ i := −1/2) so that one has as variables It is possible to have all spins up or all spins down such as ( (x 1 , σ), ..., (x m , σ) ) .Usually, the electron spins are used implicitly and they are omitted in the notations.It is possible to consider the spins as part of the variables or to fix them beforehand.It is also possible to have spinless computations which have their limitations (Ciarlet et al., 2003).For the eigenvalues, one has is the set of permutations over {1, ..., m} and sgn(P) designates the signature of a permutation P. The antisymmetrization operator A has the properties that it commutes with the Hamiltonian operator and that for an m × m matrix M one has A(MΦ) = det(M)A(Φ).The Hartree-Fock approach is the variational formulation on the Hamiltonian operator (10) where the trial functions are antisymmetric functions.Some simplifications of the stationary Hamilton operators have already been proposed.For the DFT(Density Functional Theory), one solves a set of equations for each electron.The similarity of the solutions is then derived from the theory of Kohn-Sham (Kohn & Sham, 1965).The Kohn-Sham formalism consists in replacing the complicated initial problem into several ones.
where V eff is the effective potential energy which depends implicitly on the total electron density ρ The problem is then reduced from dimensions 3m to m sets of 3D smaller problems.The influence of one electron with respect to the other electrons is measured by the total electron density.
These approaches enable the treatment of Hamiltonian problem even for an electronic structure admitting a large number of particles on a single desktop.The eigenvalue problem in ( 22) is nonlinear because its variational operator ⟩ depends on ρ which in turn depends on ψ i .It is solved by using a sequence of the linear eigenvalue problems SCF (Self Consistent Field).The effective potential is constituted of the Hartree potential V H , the exchange correlation potential V XC and the external electrostatic field such as in which the Hartree potential is the inverse of the Poisson operator such as ∆V H (x) = −4πρ(x).The main feature of DFT is that one has to approximate the potential by using some correction terms known as exchange-correlation potential (Perdew & Wang, 1992;Perdew & Zunger, 1981).That is usually done by LDA (Local Density Approximation) or GGA (Generalized Gradient Approximation).Those expressions contain some parameters which are optimized or obtained from some experimental measurements in term of the Weigner-Seitz radius.Analytic expressions of the correlation energy are only known in a few special cases which mainly consist of the high and low density limits.The external electrostatic field potential V ext is provided by the kernel The above exchange-correlation potential is related to the exchange-correlation energy by V XC = δE XC /δρ where one expresses E XC = E X + E C as the exchange and the correlation parts.In term of the exchange-correlation energy density ε XC one has For the local density approximation (LDA), the exchange energy density is expressed as Once the solution E i to ( 22) for all i = 1, ..., m becomes known, the Khon-Sham approach uses the approximation to E of ( 1) by The main improvement from LDA to GGA is that the exchange-correlation energy does not depend only on the total electron density but also on its gradient such as In the computations, one represents Ψ as a set of basis which are usually plane waves (Wimmer et al., 1981), Finite Element Method (Toda et al., 2010;Suryanarayana et al., 2010) or wavelets (Genovese et al., 2011).

Unobserved Information and Gaussian Regression
In this section, we will describe the main points about the unobserved Gaussian regression.We are going beyond the usual reconstruction of a function f by using the observations {y i } such that y i = f (x i ) + ε i because y i is unobserved in our situation.During the data fitting, we do not reconstruct a global function with respect to the atom coordinates.Instead, we use implicit variables q of dimension D. That is to say, for every atomic system R having n atoms we have the new variables q α (R) = [q α 1 (R), ..., q α D (R)] ∈ R D for each atom x α ∈ R. In this theoretical description, q α (R) is supposed to be a general mapping depending on the atom coordinates.But for the applications, they might have some physical significations such as the Coulomb matrix which we will encounter toward the end of this paper.We consider many observations involving atomic systems {R i } of number N. For each observation, we have a corresponding total energy E i tot such that in which ε i are measurement imperfections.The observed data are {E i tot } which represent the total energy for quantum application while the unknown is the function E loc .We sometimes use the term total energy to refer to the sum of all local energies.It is not to be confused with the usual total energy which is the sum of the kinetic energy and the potential energy.The quantum motivation for this kind of reconstruction is that we are interested in the local energies E loc in applications including molecular dynamics, whereas it is impossible to obtain some samples {E i loc } from an electronic structure packet.The only values that a quantum simulator provides are a set of total energies {E i tot }.It is possible to provide gradient values ∇E tot and Hessian values ∂ 2 E tot /∂x p ∂x q but we do not treat that yet in this current article because using the gradients and Hessians is computationally very intensive.Fitting using gradient and higher order derivatives will be reported in a subsequent paper.For the next stochastic deduction, we will use the following Matern correlation function where (ν, ℓ) are positive hyper-parameters and K ν is the modified Bessel function.In dimension D, the spectral density of the Matern function is If ν tends to infinity, one obtains the squared exponential case.For the special case of half-integers such as ν = p + 1/2, one obtains a product of an exponential and a polynomial of order p such as In the particular cases where ν = 3/2 and ν = 5/2, one has The number of atoms in the i-th observation is n i for i = 1, ..., N and we denote M := ∑ N i=1 n i .Let us denote by Q the following of implicit variable . . . . . . . . . . . .
which is a matrix admitting the size (M×D) where D is the dimension representing the number of implicit variables q per atom.Let L be a binary matrix of size (M × N) having entries 0 or 1 such that only the Coulomb entries in the i-th observation are relevant for the i-th row of L T .Denote Q = [q i ] M i=1 so that one has the relation E i loc := E loc (q i ) where {E i loc } M i=1 are unobserved samples.As a consequence, one has the following covariance relation for noise-free observations For a noisy data E i tot = E tot (R i ) + ε where the additional noise follows the Gaussian distribution ε ∼ N(0, Introduce the next vectors of observed and unobserved inputs In the presence of noise, the marginal likelihood is the integral of the likelihood and the prior such as and one has no noise, then the marginal likelihood is which implies by taking the logarithm marginal likelihood In the presence of noise ε ∼ N(0, σ 2 n ) by using (36), the log marginal likelihood becomes The first term is the complexity term.The last term − N 2 log(2π) is the normalization term.For a given test value Q * , the predicted estimation of the local energy follows the Gaussian distribution as That is to say, for Q and E tot , the predictive distribution for a covariate vector Q * is Gaussian having the following expectancy and variance where those expressions are dependent on some set of hyper-parameters.The main objective is to theoretically optimize the marginal likelihood.In practice, using the log marginal likelihood is more efficient to implement.Now, we would like to investigate the derivative of the log marginal likelihood with respect to a hyper-parameter θ.Let the value of r for given q p ∈ R D and q s ∈ R D be in which λ 1 ,..., λ D are hyper-parameters which should be determined by using the log marginal likelihood.One has θ = log(ℓ) or θ = log(λ i ) where i = 1, ..., D. In order to apply an efficient nonlinear optimization to the former functional, we need the partial derivative of the inverse ] −1 with respect to a hyper-parameter θ.From the fact that ∂ ∂θ one deduces ∂ ∂θ which cannot be simplified any further as L is not invertible.Similarly, one has the relation By using the Matern function k from ( 29), the partial derivative of a covariance matrix entry with respect to a hyperparameter θ is , where r(q p , q s ) = This is singular in the situation where q p = q s which occurs on the diagonal entries of the covariance matrix.In fact, one has In order to suppress the singularity, one utilizes R(q p , q s ) := r 2 (q p , q s ) such that ∂ λ µ R(q p , q s ) is regular everywhere.
In addition, k(R) where where q MATERN is a polynomial.Therefore, one has which is regular for all values of r(q p , q s ).The partial derivatives of the above estimation with respect to components of q can be expressed analytically because we have used linear entities.

Molecular Dynamics of Large Systems
In general, by considering a real valued Lagrange function L(t, x 1 , ..., x n , v 1 , ..., v n ) and several continuously differentiable functions q 1 ,...,q n , one obtains the following Euler-Lagrange equation for each k = 1, ..., n We consider a time-dependent system consisting of n atoms where the k-th atom admits the Cartesian coordinates r k = r k (t) depending on time.By considering the MD Lagrange function where M k is the mass of the k-th ion, the Euler-Lagrange equation ( 51) yields the following equation of motion which is described as an ODE (ordinary differential equation) for each component j = 1, 2, 3.That expression is the second law of motion where [r k,1 (t), rk,2 (t), rk,3 (t)] denotes the acceleration while the force ] is applied on the k-th atom.By considering the kernel V pot of the potential energy, one should normally obtain the relation where Ψ corresponds to the ground state energy.But the Hellman-Feynmann rule states the independence of the force on the derivative of the ground state energy such that one eventually obtains Every MD iteration consists of the following steps.First, one determines the atomic force F k .Then, one has to solve the above equation of motion with the help of any ODE solver.The usual methods for the ODE solver are based on variants of Runge-Kutta or Verlet algorithm.Finally, one updates the atomic position r k .For an expression A = A(t) depending on t, we designate the time average by Since the MD trajectory involves only a large number of discrete values, we utilize the following discrete counterpart for L number of MD iterations In the case of two-body potentials where the potential function V depends only on the inter-atomic distances r i j = ∥r i − r j ∥, let us review some MD thermodynamic expressions.The pressure and the temperature are respectively defined by the expressions where Ω, n and k B are respectively the volume, the number of particles and the Boltzmann constant.Equation ( 59) relates the regulation of the temperature with the kinetic energy.The radial distribution function for nonzero radius r is which quantifies the general distribution of the atoms with respect to their corresponding neighbors.We have in addition One obtains finally the total potential energy PE = 2ρπn ∫ r 2 g(r)V(r)dr which can be used to verify the total energy conservation when the kinetic energy is added to it.

Practical Results
In this section, we would like to give some details about the results of the practical implementation of the formerly described theoretical method.We will present results related to mathematical models as well as to quantum application for molecular dynamics.The implementation uses C/C++, BLAS/LAPACK and NLOPT.The BLAS packet is used for the fast vector operations.We use LAPACK for the linear operations such as Cholesky factorization and dense matrix solvers.The code is a very developed version of the matlab implementation provided in (Rasmussen & Williams, 2006).The new additional enhancements from the matlab version consist of the following features.First, using NLOPT provides a lot of improvements as compared to the original matlab nonlinear conjugate gradients.That can be observed when both the number of points and the dimension become large in the nonlinear optimization of the hyper-parameters.Furthermore, the gradient of the kernel-based approximation can be evaluated by using analytical expressions instead of using finite difference.The initial guess of the hyper-parameters is provided by the users.One can consider the determination of the final hyper-parameters as an unconstrained nonlinear optimization.We use NLOPT for the nonlinear operations (Nocedal, 1999) for the optimal hyper-parameters in the log marginal likelihood (41).NLOPT supports diverse nonlinear optimization operations (Johnson, 2010) in which local optimizers are involved.A local one searches only inside a neighborhood of a certain provided starting initial guess.The optimizations are performed by using derivative-free or gradient-based algorithm which are available in NLOPT.Derivative-free algorithms include BOBYQA (Bond Optimization BY Quadratic Approximation), COBYLA (Constrained Optimization BY Linear Approximation), NEWUOA (NEW Unconstrained Optimization Algorithm).Gradient-based methods include MMA (Method of Moving Asymptotes) and LBFGB (Limited memory Broyden-Fletcher-Goldfarb-Shanno).
As a first application, we consider the reconstruction of a mathematical model where the exact solution is explicitly known.In an atomic system R constituting of n atoms (r 1 , r 2 , ...,r n ), the Coulomb matrix [ q α,β ] admits the following entries for α, β = 1, ..., n: in which Z α is the nuclear charge of the α-th atom.In our situation, we consider only homogeneous systems consisting of one element type such that Z α ≡ Z is a constant.We fix any positive number D which will then be the implicit dimension.That is to say, for each atom r α , the implicit variables ] will be the largest Coulomb entries involving the atom r α .Except for the Coulomb entry q i = 0.5Z 2.4  i , the largest Coulomb values correspond to the closest atom neighbors.The purpose of the first application is to reconstruct the local energy E loc which is a D-variate function.In order to be able to compare the reconstructed function with the original one, we consider an exact function which maps an atomic system R = (r 1 ,..., r n ) to , where (64) It is possible to evaluate the errors with respect to the local energy or the global one.As in the real situation, the above exact function is unchanged by reordering the sequence of the atoms and the Coulomb matrix.The reconstruction of the implicit function E loc from the total energy values follows the theoretical prediction from section 2.2.The inputs are the set of total energies {E i } in which each total energy E i corresponds to one atomic system R i .We do not put any restriction on the size of the atomic systems for that two different atomic systems R i and R j may contain different numbers of atoms and that D is not necessarily the whole n i .On the other hand, the size D of the Coulomb values per atom has to remain constant during the whole simulation.Otherwise, it would be impossible to assemble the covariance matrices.For our simulation, each coordinate of the atom x α ∈ R is randomly defined.Note that the implicit size is much larger than the number of atomic systems so that much larger capacity of the computer memories are required than in the usual fitting (Rasmussen & Williams, 2006).We examine the accuracies for each group of atomic systems categorized according to their sizes n i = 12, ..., 15.In our case, the implicit dimension and the atomic number are D = 4 and Z = 5 respectively.We observe that we have generally a good accuracy improvement when the number N of the atomic systems increases irrespective of the size n i of each atomic system because the value of closest neighbors D is fixed.
We would like now to investigate the influence of the implicit dimension D on the accuracy of the local energies.
As in the former test from Table , we take again 5000 testing atomic systems.We vary the value of the dimension in the range D = 4 till D = 7 which controls the number of Coulomb values that are considered for each atom.The results of the test are collected in Table where it can be observed that the dimension affects the accuracy.A higher dimension is somewhat more difficult to approximate than a lower dimension.Nonetheless, it can be observed that the error reductions are somehow comparable for all dimensions.As the number of the atomic systems grows, the improvements of the accuracy are relatively of the same size for all considered dimensions.In order to observe the error reduction more clearly, we would like to analyze the accuracy of the total energy.In Figure 1(c), we display the error plots in function of the number of training atom systems.Both the axis for the number of atomic systems and the axis for the error in the total energy are logarithmically scaled.We examine two different cases admitting dimensions D = 4 and D = 6 which correspond respectively to the blue curve with triangular markers and the red curve with square markers.The results exhibit that the performance is efficient for both test cases.The case D = 6 is somewhat more difficult to approximate than the case D = 4 but both situations show that a good accuracy is possible to be reached as more and more atomic systems are added into the training configurations.The purpose of the following test is to apply the potential energies obtained from the fitting process into a molecular dynamic model.We focus only on a two-body potential for the computation of the potential energy.In fact, we utilize the following Lennard-Jones potential energy (Tadmor & Miller, 2011) for the molecular dynamics: We employ reduced units so that all involved entities are dimensionless.For example, the temperature T is factored into T k B /ε.In addition, ε and σ admit unit values in our case while the method is still applicable to more general parameters.This corresponds to the general theory in (26) where the unknown function is E LJ .We replace only the sum in ( 26) by a sum with two indices which can be organized lexicographically in order to obtain a sum with a single index.The implicit variables in this case are q α,β := ∥r α − r β ∥ so that we have the dimension D = 1 for the two-body potential.The atomic forces at each MD step are obtained by taking the analytical partial derivatives of the Gaussian fitting which can be computed explicitly.Taking finite difference from the potential energy is known to possess a diverging energy conservation.The purpose of the test is to compare the MD simulation using the direct potential energy on the one hand and the one using the Gaussian fitting on the other.The total energy which is the sum of the kinetic energy and the potential energy should be kept approximately constant during the whole molecular dynamic steps.We want to compare the results pertaining to the energy conservations.The curves for the total energies are shown in Figure 1(d) for several MD steps.One can observe that the total energies for the two methods align well.The closeup shows that at each MD steps, the two results differ somewhat but the general energy conservations remain practically unchanged.
We want to examine the influence of the temperature on the whole molecular dynamic simulation.The adjustment of the temperatures is regulated by the scaling method of the kinetic energies which is a well-known method for regulating temperatures.In fact, one executes the molecular dynamic simulation for sufficient number of steps until the equilibrium is reached before one applies a constant scaling factor to the kinetic energy.An illustration of that situation is depicted on Figure 1(g) where 200 MD steps are executed to reach the equilibrium.Three different temperature adjustments are displayed there.One of them corresponds to a heating adjustment while the other two are used for cooling.These temperatures are respectively T = 0.78T equil , T = 0.62T equil and T = 1.28T equil .
When the desired cooling or heating temperature values are attained, one applies again a normal MD simulation.
In our case, the whole simulation consists of 2000 molecular dynamics steps.We want now to investigate the radial distribution function which was met in (60).We use the local potential energies which have been separated on atoms by using the former method.The radial distribution function quantifies the general distribution of each atom with respect to its neighbors.The values have been averaged with respect to the whole molecular dynamic trajectories.The results of the RDF tests are depicted on Figure 1(h) where the horizontal axis describes the average radius from an atom.We consider four situations which correspond to four different thermal values.In general for all values of the temperature, the RDF value is zero until the ideal radius is reached at the position where the RDF attains its peak.That is followed by some values which alternate in the neighborhood of the unity.The heights of the peak are very much affected by the temperature.The peak values are 1.7, 2.3, 2.8 and 4.9 for the temperature T = 285.71Tequil , T = 2.1930T equil , T = 0.2917T equil and T = 0.0604T equil respectively.By using the proposed potential, one obtains the phenomenon of atom clustering in the case of freezing temperatures.By executing sufficiently many molecular dynamic steps, the atoms automatically form clusters by accumulating close to one another.In Figure 3, we can observe some instances of such cluster formations.Atoms which are close enough to one another are depicted with the same colors.That is another validation that the proposed Gaussian potential energies can efficiently simulate molecular dynamics.This clustering phenomenon can also be confirmed by the radial distribution functions in Figure 1(h) because a high peak averagely reflects the phenomenon of dense neighbors next to each atom.

Conclusion
We presented an approach which approximates the total energy and which decomposes it into atomic contributions.
The approach is based on stochastic method where the expectancy and the variance are computed.The main purpose is to determine the local energy function while the only observed inputs are the total energy.The method has been tested by using both a mathematical model and a molecular dynamic simulation.For the mathematical model, the accuracy between the exact solution and the fitting results has been investigated while we use two-body potentials for the molecular dynamic model.With the help of the new method, we simulated atom clustering in case of freezing temperatures.In addition, the radial distribution function has been examined for various values of the temperatures.

Figure 1 .
Figure 1.Motivation: (a)Local energy computation and molecular dynamic simulation (b)Potential energy in a graphene fracture Figure 2. (a)Error in the total energy (b)Energy conservation: Lennard-Jones vs. unobserved fitting potential

Figure 4 .
Figure 4. Using the unobserved Gauss reconstruction in reduced units (dimensionless): (a)Temperature adjustment and fluctuation (b)Radial distributed functions for different temperatures

Table 1 .
Errors in term of number of atomic systems with four Coulomb coordinates per atomAs a first test for examining the accuracy, we consider several training atomic systems where the number of atoms in each system R i varies from n i = 12 to n i = 15 which is chosen by some random number generator.To test the accuracy, we consider additional 5000 atomic systems whose sizes also vary from 12 to 15.The results of the computations are collected in Table which contains the accuracy of the local energy E loc .It shows the influence of the number N of training atomic systems and the implicit size which is the size of the involved covariance matrices.

Table 2 .
Influence of the implicit dimension (Number of the closest neighbors