Multilevel B-Spline Repulsive Energy in Nanomodeling of Graphenes

Quantum energies which are used in applications are usually composed of repulsive and attractive terms. The objective of this study is to use an accurate and efficient fitting of the repulsive energy instead of using standard parametrizations. The investigation is based on Density Functional Theory and Tight Binding simulations. Our objective is not only to capture the values of the repulsive terms but also to efficiently reproduce the elastic properties and the forces. The elasticity values determine the rigidity of a material when some traction or load is applied on it. The pairpotential is based on an exponential term corrected by B-spline terms. In order to accelerate the computations, one uses a hierarchical optimization for the B-splines on different levels. Carbon graphenes constitute the configurations used in the simulations. We report on some results to show the efficiency of the B-splines on different levels.


Introduction
Nanotechnology is a very important field which has emerged in the last decades and developed very quickly in several directions.It has important applications in various disciplines including aircraft, automobile, electronic and medical engineerings.Nanomaterials admit several important properties which can be exploited in applications.For instance, electric conductivity of nanomaterials is applied in electronic components so that the materials conduct electricity more efficiently than diamonds.Thermal resistivity of nanomaterials can be used to reduce or accelerate heat conduction.It also has a good thermic property so that materials can be designed to resist heat at a very high intensity.Graphene has obtained a significant attention from scientists in the last decades for several reasons.Its material property can be controlled for that it can become a stronger material than steel.The objective in this paper is to use an accurate and efficient fitting of the repulsive energy instead of using a standard parametrization.Many approaches have been proposed to represent empirical estimation of repulsive terms.Before presenting our method, let us briefly survey some traditional repulsive methods.Molecular dynamics employing the Lennard-Jones potential have been well understood so far.It is based upon the well-known potential ( ) which is decomposed into attractive and repulsive components.Since it is only expressed in terms of the inter-atomic distances, it is easy to handle.Due to the simple expression of the potential, it can be differentiated easily and it is not computationally expensive to evaluate.Another important parameterization is the Harrison parametrization: The screened Harrison parametrization is an improvement of the former one as provided by ( ) ( ) where c r controls the range of the interaction and µ is some parameters obtained by experiments or a fitting process.Sawada parameterization uses the expression ( ) ( ) ( ) The most currently used parametrization is the GSP parameterization (Goodwin-Skinner-Pettifor) which is expressed as where n , e n and c r are fitting parameters.Several other methods have been also suggested to achieve some desired properties.Some approaches use certain combinations of known ones.
Our motivation is to generate a system which is both accurate and fairly inexpensive to evaluate.We are interested in graphenes and its properties including energy, force and elastic stress.Geometrically, graphenes admit a honeycomb pattern in form of repeated organized hexagons as illustrated in Figure 1 For the generation of the unit cell, one needs a translational vector T perpendicular to the chiral vector C .Let d designate the greatest common divisor of n and m .Define : 3d In the following sections, scaling a graphene amounts to enlarging the unit cells by scaling its primitive vectors.The coordinates of the centers of carbon atoms in the unit cell provided as fractional coordinates within [ ] 0,1 remain unchanged.We are interested in the property of the graphenes as they are confined or stretched as illustrated in Figure 1(c) where we consider a graphene of chirality ( ) 2,1 .For significantly confined graphenes, the repulsive energy is very large.For extremely stretched ones, the repulsive energy vanishes so that the energies are the sum of the energies of all constituting atoms.

Quantum Repulsive Representation
We consider 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.The atoms move comparatively slower than the electrons.Thus, we treat the time-independent Hamiltonian operator with respect to the a set of nuclei which are supposed to be stationary: x x (6) where the coordinates of the i -th electron are denoted by , , . The above formula describes the kinetic energy along with the nuclear-electron interaction and the inter-electron interaction.Several simplifications of the stationary Hamilton operators have already been proposed.Our proposed potential energy uses two of such simplifications which we survey below.
For the DFT(Density Functional Theory), one solves one equation for each electron.The Kohn-Sham formalism [1] consists in replacing the complicated single problem into several simpler ones.For each 1, where eff V is the effective potential energy which depends implicitly on the total electron density ( ) ( ) x x .The problem is then reduced from dimensions 3 e N to e N sets of smaller problems of dimension 3D .The influence of one electron with respect to the other electron is measured by the total electron density.These approaches enable the treatment of Hamiltonian problem even for an electronic structure having a large number of particles on a single desktop.The eigenvalue problem in (7) is nonlinear because its variational [2] [3] operator ( ) x Ψ 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 H V , the exchange correlation potential XC V and the external electrostatic field such as in which the Hartree potential is the inverse of the Poisson operator such as For its evaluation, either a Poisson problem is solved or one convolves with the Green fundamental solution such as x x x x .The main feature of DFT is that one has to approximate the potential by using some correction terms known as exchange-correlation potential [4] [5].That is usually done by LDA (Local Density Approximation) or GGA (Generalized Gradient Approximation).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 ext V is provided by the kernel where For the local density approximation (LDA), the exchange energy density is expressed as ( ) ( ) ( ) [ ] ( ) Analytic values of the correlation energy density are only known for some extreme cases.For the high density limit, the exchange correlation energy density is approximated by when the Weigner-Seitz radius s r is very small.For the low density limit where s r is very large, one has ( ) ( ) ( ) ( ) For other values of s r , some interpolation of those extreme values is considered.For example, by using the VWN-approximation (Vosko, Wilk, Nusair) as in [6], one has Once the solutions i E to (7) become known for all 1, , e i N =  , the Khon-Sham approach uses the approximation to E of ( 6) 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 [ ]( ) x .As a second simplification, we survey the semi-empirical (SE) method using Hueckle method.Consider the spherical coordinates ( ) r θ φ such that ( ) The atomic orbitals sharp (s), principal (p), diffuse (d) and fundamental (f) correspond to linear combinations of ( ) respectively.The basis functions centered at the origin are [7] defined by where the radial function The parameters = − a a a .In fact, the overlap matrix entries can be expanded as , , = ∑ a where one uses the inter-atomic distance ij i j d = − a a .Computing the integrals by quadrature is too expensive.
One stores the expansion coefficients ( ) C α β a .The value of ( ) are stored in Slater-Koster table.It does not depend on the coordinates of i a and j a but only on the inter-distance (see [7] for similar discussion).For the Hamiltonian of the SE Hückle, the on-site situation with respect to the center i a is ( ) where E α is approximately the eigenenergy for index α .The off-site term is 2 The value of x .This partial differential equation needs to be solved for every evaluation of the Hartree term ( )( ) In the Atomistix Toolkit package [7], that is solved by a fast multigrid solver.The coefficient ( ) ε x is a dielectric coefficient [8] and As a matter of fact, the SE empirical method is much more efficient than the DFT method in term of computational speed.But the DFT computation produces much more accurate results.As a consequence, one searches a certain correction term for the SE method in such a way that the resulting method keeps the efficiency of the SE method while approximating the quality of the DFT approach.The ultimate objective is thus to find a repulsive term to add to the SE energy as described below.We want to generate a repulsive term which conserves most of the properties from the DFT computation.For a configuration { } 1 , we intend to conserve the energy ( ) . In addition, we are also interested in approximating the forces.For each atom ( ) , , , the corresponding force is ( ) , , In addition, we focus also on the elastic property of the graphenes [9].In general, this property determines the rigidity of a graphene when a traction is applied on it.The strain tensor which is is represented in the longitudinal, transversal and normal components.The stress σ is also represented in a similar tensor way.The strain is related to the displacement u having components i u by ( ) The correlation between the strain ε , stress σ and displacements u is governed by some elasticity equation [9].Practically, the stress contains implicitly some property of the second derivatives of the energy for the reason that it is the derivative of the energy with respect to strains which are functions of the gradients of the displacements ( ) For a set  of graphene configurations, the ideal objective functional for the nonlinear optimization is . uniformly.Afterwards, one refines gradually in the vicinity of the optimal scaling factor ( ) opt λ c of the configuration c .The principal objective for that construction is to accumulate many points in the neighborhood of the optimum ( ) opt λ c .The determination of the stress is computationally more intensive compared to the computation of the energies.That situation holds even for the case of semi-empirical Hueckle method.The computation of stress for the DFT case is even more intensive but it needs only be done once and stored during the whole optimization.As a consequence, one needs only to handle elastic properties at a few positions in the course of the optimization computation.Otherwise, the whole optimization execution would be too slow since the evaluation of the objective functional would be very intensive.For example, the stress is only applied in the neighborhood of the minimal energy in our computation.Generally, λ ∈ Λ c is of the same importance.The vicinity of the optimal scaling factor opt λ is more valuable because the equilibrium takes place there.As a consequence, one introduces some positive weights to the scaling factors.For our implementation, we used some Gaussian functions centered at the optimal value added by some minimal shift shift > 0 δ such as ( ) The purpose of shift δ is to prevent the value of λ ω from being practically zero when λ is far from opt λ .
In our case, we have taken the parameter values to be 0.005 c σ = = and shift 1 δ = .Since the objective function ( 10) is very intensive to evaluate, we use in practice its simplification where the forces are provided by finite difference of the energy.Now we would like to describe the parameters with respect to which the nonlinear optimization is performed.The semi-empirical energy with zero repulsive term 0 SE E behaves as a pure attractive energy.That is, in order to obtain an energy comparable to the DFT energy, one appends a .
That is, the function  acts pairwise on the carbon atoms with nuclei coordinates i a and j a such that ( ) ( ) a a .In other words, the whole process amounts to replacing the repulsive term of the SE energy by an optimal potential energy ( ) ∑  a a .We search for the optimal pair potential function in the form ( ) ( ) in which N designates B-spline basis functions such that we obtain an energy that behaves very similarly to the DFT in term of energy, force and stress.
In the expression (12), the function e bt a − captures the general behavior of the pair potential function  .The role of the is to correct the small imperfection produced by the principal function e bt a − .Some cut-off radius C r is used so that the pair potential function  vanishes beyond that value.In our situation, a cut-off radius of about 4.0 C r = Å suffices completely.In order to obtain zero value and derivative , we insert a short transition function next to the cut-off radius C r .One extends  next to the cut-off radius C r by a polynomial so that one obtains a smooth transition toward zero.Since the unknown pair potential function is partly expressed in B-spline basis as in ( 12 . One defines the B-spline basis functions for 0, , , , , , in which the truncated power functions ( ) We only focus on B-splines which are internally uniform: except for the boundary multiple knots, all knot entries i ζ are uniformly spaced.The integer k controls the smoothness of the B-spline for which the resulting function admits an overall smoothness of 2 k −  so that the case 1 k = corresponds to discontinuous piecewise constant functions.The integer n controls the number of B-spline functions.In Figure 3(a), we see some illustration of B-spline bases on an internally uniform knot sequence.Figure 3(b) displays an instance of a B-spline curve defined on [ ] 0,1 .In Figure 3(c), the knot sequence has been refined uniformly by increasing n to 2n while keeping 2 k = .That is achieved by introducing a new knot entry between every two knots of the B-spline in the former Figure 3(b).For our application, we insert several knots at once so that the new knot sequence is again internally uniform.A new knot entry is inserted between two consecutive old ones.The evaluation of B-spline functions is not calculated by using the above definition but rather by means of the de-Boor algorithm.Since the knot sequence is internally uniform, we use the notation ).We will describe next the procedure of inserting new knots into existing ones.That is important when one needs to increase the degree of freedom in the pair potential function in (12).The principal objective is to efficiently express a function defined on the coarse knot sequence in term of B-splines on a fine one.Consider two knot sequences ( ) ) where : are evaluated by using the recurrence In our simulation, we took 3 k = which corresponds to continuously differentiable pair potentials.In Figure 4, we observe some illustration of such knot insertions.Not only the two B-spline functions admit the same image but their parametrizations from their interval of definition [ ] 0,1 are completely identical.The whole process of the determination of the repulsive energy is performed in increasing levels as follows.First, one determines the optimal value for e bt a − without the B-spline part in (12) by using a global optimizer.
Then, one fixes the resulting optimal values of ( ) , a b during the subsequent computation.Second, one searches the optimal B-spline ( ) where 4 n = by starting a local optimization with the initial guess , , . Now, one repeats the following steps iteratively.Inject the optimal value ( )

Computer Implementation
In this section, we report on some practical results of the formerly proposed method.The implementation of the method was realized by combining ATK, NLOPT and python.The ATK (Atomistix ToolKit) has some GUI extension well known as VNL (Virtual NanoLab).We use NLOPT for the diverse nonlinear optimization operations [10] in which both global and local optimizers are involved.A global optimizer searches for the best parameters among all possibilities while a local one searches only inside a local neighborhood of a certain provided starting initial guess.For the global optimizer, we use the NLOPT option GN-CRS2-LM standing for Controlled Random Search with Local Mutation.The local optimizers are performed by using BOBYQA algorithm which is an efficient gradient-free method available in NLOPT.In order to facilitate the combinations of options, we implemented several python classes.The class for the reference configurations organizes the graphene structures to be used together with their respective optimization weights.There is also a class for the optimization parameters specifying the property such as orders and levels of B-splines as well as the abortion criterion.It controls the contribution of the energy, force and stress ( ) , , E F S µ µ µ in the optimization functional (10).The construction of the sets ( ) Λ c as well as the interval for the range of interest has been equally supported by some python classes.In order to save computations, one needs to precompute and store the data for the DFT as well as the semi-empirical with zero pair potential.
As a first test, we consider multiple computations for different configurations of graphenes.The configuration n  is based upon the first index n of the chirality parameters ( ) , n m where m is allowed to vary.That is, each configuration n  is composed of all graphenes admitting chirality ( ) , n m such that 0 m n ≤ ≤ .In Figure 4(a), we observe some comparisons for graphenes in n  where 1 n = .Most values align on the diagonal which implies the agreement between the outcomes provided by the DFT and the SE methods.Similar tests for graphenes where 2 n = and 3 n = are depicted respectively in Figure 4(b), Figure 4(c).The resulting SE energies do not exactly provide the same results as the DFT but the current SE energies should be more reliable in comparison to the empirical potential in (1)-( 5) which contain very few parameters.In addition, the speed of computation is much faster for the currently presented SE than the one for DFT.In fact, the execution time of the SE in comparison to DFT has a speed of factor 10 or more.Due to that acceleration gain, the method is in many aspects good to attain efficiency.If the accuracy is not satisfactory, then one has to use the direct DFT with the cost of much more computing time.
As a further test, we investigate the decrease of the objective function with regard to the B-spline levels.We consider again the three configurations n  above for 1, 2,3 n = . The results of such tests are displayed in Table 1 where the initial line describes the SE with zero pair potential (PP).The next one is the SE with exponential pair potential ( ) e bt t a − =  without B-splines.The following ones are the pair potentials with more and more B-splines as in (12).The error barely drops down after level 4 for all graphene configurations n  .In fact, the minimal value of the function in (10) is not always zero.As a consequence, one cannot expect an arbitrarily accurate approximation.As a next test, we consider the complex band structures for using the DFT and SE computations whose results are respectively displayed in Figure 5(a), Figure 5(b) for the graphene with chirality ( ) 1, 0 .The plots depict band lines which are not shown as continuous curves but as sets of sampling points.The points which are purely real and explicitly complex are depicted in red and green respectively.In order to provide more validation for the efficiency of the proposed method, some comparison of the elastic properties was performed when computed by means of the DFT and SE methods.In Figure 6, we observe the elastic properties corresponding to the two methods.In general, the stress tensor σ is presented in three directions similar to (9).Nevertheless, we omit the normal components of the stress tensor σ in this particular

Conclusion
A method was presented to determine the optimal pair potential for the repulsive quantum energy.We concentrated on configurations which are constituted of carbon The method was based upon hierarchical B-splines layered on different levels.The principal objective function consists of terms involving not only energies but also forces and elastic stresses.Several computer results validate the reliability of the newly proposed method as compared to outcomes from Density Functional Theory.
(a).They are controlled by the chirality which is a couple of integers ( ) , n m so that 0 m n ≤ ≤ .In the case n m = , one has an armchair graphene while 0 m = corresponds to the case of a zigzag graphene as in Figure1(b).Suppose 3 a designates the carbon bond length of the graphene.Define a and b the directive vectors of the honeycomb describing a 2D-lattice so that

Figure 1 .
Figure 1.(a) Nanosheet: chiral vector C and translational vector T ; (b) Armchair and Zigzag graphene; (c) Ground state energy for graphenes.Confining and enlarging the volume of the graphenes.
The above exchange-correlation potential is related to the exchange-correlation energy by XC XC and the correlation parts.In term of the exchange-correlation energy density XC ε one has

Λ
c are sets of scaling factors with respect to the reference configuration ∈  c for the energy, force and stress respectively.Now, we show the construction of prescribes the range of interest.That interval contains the optimal factor geometry optimization.The construction is performed in several steps as depicted in Figure2.As a first step, one refines the interval [ ]
), we recall briefly some important properties of a B-spline setting.It is in fact a very flexible way of representing piecewise polynomials on any interval of definition [ ] , a b .Consider two integers , n k such that 1 n k ≥ ≥ .Suppose the interval [ ] , a b is subdivided by the knot sequence For both knot sequences, the smoothness index k is conserved intact.The following discrete B-splines enable the expression of a coarse basis i N ζ as a linear combination of the fine basis

Figure 3 .
Figure 3. (a) B-spline bases; (b) Original B-spline; (c) a finer B-spline which has the same parametrization as the original B-spline.
above knot insertion technique where n is increased into 2n .Apply then a local optimization with respect to ( )

Figure 6 (
Figure 6(b) where one observes the alignments of DFT and SE elasticities.The two results for ( ) 1, 0 and For the nondiagonal values or off-site cases centered at two different atoms i a and j a , the entries are computed by a Slater-Koster table lookup where the values are functions of the interatomic coordinates ij

Table 1 .
Errors at each B-spline level.case of the graphene configurations which are planar.In addition, the stress components LT