From Gradient Elasticity to Gradient Interatomic Potentials: The Case-Study of Gradient London Potential


Motivated by the special theory of gradient elasticity (GradEla), a proposal is advanced for extending it to construct gradient models for interatomic potentials, commonly used in atomistic simulations. Our focus is on London’s quantum mechanical potential which is an analytical expression valid until a certain characteristic distance where “attractive” molecular interactions change character and become “repulsive” and cannot be described by the classical form of London’s potential. It turns out that the suggested internal length gradient (ILG) generalization of London’s potential generates both an “attractive” and a “repulsive” branch, and by adjusting the corresponding gradient parameters, the behavior of the empirical Lennard-Jones potentials is theoretically captured.

Share and Cite:

Parisis, K. , Shuang, F. , Wang, B. , Hu, P. , Giannakoudakis, A. and Konstantinidis, A. (2020) From Gradient Elasticity to Gradient Interatomic Potentials: The Case-Study of Gradient London Potential. Journal of Applied Mathematics and Physics, 8, 1826-1837. doi: 10.4236/jamp.2020.89137.

1. Introduction

A simplified version of gradient elasticity theory (GradEla) introducing an extra gradient term (the Laplacian of Hookean stress) in the classical law of linear elasticity [1], has been shown to dispense with various difficulties encountered in the past. Among the advantages stemming from the use of the GradEla model were the removal of singularities from dislocation lines and crack tips, as well as the possibility of conveniently interpreting elastic size effects (e.g. [2] [3] and references quoted therein).

The main feature which made GradEla especially robust was the observation that, under certain conditions [4], solutions of GradEla can be obtained in terms of existing solutions of classical elasticity by solving a non-homogeneous Helmholtz equation. This observation enabled to obtain explicit and easy-to-use non-singular solutions for dislocations and disclinations, as well as for cracks. An account of these developments can be found in a number of previous articles on the subject [5] [6] [7] [8].

In a more recent article [9], the gradient approach was extended to modify the classical Newton’s Law of gravitation, leading to an unexpected preliminary result: the possibility of interpreting the “strong force” of subatomic elementary particles on the basis of a gradient generalization of Newton’s gravitational law. In a related presentation in a soft matter symposium at the University of Florida/Gainesville [10], the question was raised [11] whether or not such a gradient enhancement for the gravitational potential can be extended to modify interatomic potentials used for multiscale simulations in solid state and soft matter calculations. This subject is currently under investigation by the Florida—Thessaloniki groups and a preliminary encouraging result is reported herein.

In Section 2, a brief review of GradEla and its implications in revisiting classical dislocation and fracture mechanics is outlined. In Section 3, a brief review of London’s potential is provided. It is noted that London’s potential is derived on the basis of quantum mechanical calculations, but it describes only the “attractive” (1/r6) interaction of the atoms/molecules considered. It holds up to a critical distance r0, after which the interaction becomes “repulsive”, such that particle “collapse” is prevented and the matter remains intact. Various empirical modifications of the London’s potential have been introduced to model both the “attractive” and “repulsive” branches of the interatomic potential. Among them, a most popular one is the Lennard-Jones potential where an opposite sign (1/r12) term is added in the classical form of London’s potential, such that an “equilibrium” minimum is obtained and the phenomenological constants multiplying the aforementioned power terms are adjusted from experimental data on macroscopic properties of the system at hand. The forms of other similar type of empirical potentials used in the literature are also listed in this section. In Section 4, the gradient generalization of London’s quantum mechanical potential is presented in analogy to the gradient extension of Newton’s gravitational potential [9]. It is shown that this generalization results in a “modified” London’s potential containing both an “attractive” and a “repulsive” branch. By adjusting the new phenomenological parameter characterizing the effect of the gradient (Laplacian) term, the behavior of various types of empirical interatomic potentials, such as the Lennard-Jones which we focus on, can be recovered. Finally in Section 5, conclusions and comments on future work are briefly discussed.

2. Review of GradEla

The classical Hooke’s Law of linear elasticity reads

σ = λ ( t r ε ) 1 + 2 G ε , (1)

where σ is the Hookean stress, ε is the linear strain and ( λ , G ) are the usual Lamé constants. A strain gradient generalization of Equation (1) can be obtained by assuming a nonlocal integral expression for the strain ε of the form

ε ( r ) = 1 V V G ε ( | r r | ) ε ( r ) d V , (2)

where V is the elementary volume considered under the macroscopic strain ε at the point r ; ε ( r ) is the local microscopic strain at each particular point r within the volume V, and G ε denotes a corresponding kernel describing the “weighted” effect on the microscopic strain ε ( r ) on the macroscopic strain ε ( r ) . By taking the Fourier transform of Equation (2), expanding in Taylor series in the Fourier space, and inverting we can replace the classical strain ε in Equation (1) with its gradient counterpart ε l ε 2 2 ε , so that Equation (1) is replaced by the GradEla constitutive equation of the form

σ = λ ( t r ε ) 1 + 2 G ε c 2 [ λ ( t r ε ) 1 + 2 G ε ] , (3)

where the new phenomenological gradient coefficient c = l ε 2 is the square of an internal length characterizing the inhomogeneity of the underlying (micro/nano) structure of the non-classical elastic material at hand. In fact, it turns out that l ε is given by the relation

l ε 2 = 1 2 | d 2 G ˜ ε ( 0 ) / d k 2 | ,

where G ˜ ε denotes the Fourier transform of the kernel G ε and k = | k | is the wave vector.In general, the sign in front of Equation (3) can be either positive or negative since the outlined mathematical procedure gives c = s i g n ( l ε ) l ε 2 ; l ε = ± 1 2 | d 2 G ˜ ε ( 0 ) / d k 2 | . Even though some early authors have used Equation (3) with the “+” sign in front of the Laplacian term, stability reasons require the “−” sign to be used. A simpler and more direct way to arrive at Equation (3) is toset G ε = 1 in Equation (2), so that ε ( r ) is the average strain tensor over the elementary spherical volume V centered at r . By performing a Taylor expansion(up to the second order) of the local strain ε ( r ) in physical space around the point r , integrating over the elementary volume V, and subsequently inverting, we arrive at the same result where now l ε 2 R 2 / 10 with R denoting the radius of V. In other words, the internal length l ε is directly related in this case with the size of the elementary volume at hand which, in the case of a polycrystal, can be identified with an average grain size.

It is noted that the second term in Equation (3) enhancing the classical elasticity law of Equation (1) is indeed the Laplacian of Hookean stress. It also turns out [4] that solutions of the usual equilibrium equation

d i v σ = 0 or σ i j , j = 0 , (4)

can be obtained in terms of classical elasticity solutions for infinite domains or finite domains under certain conditions on the boundaries. In fact, for traction boundary conditions, the displacement field u and the strain field ε of GradEla turn out to be given in terms of solutions of the following inhomogeneous Helmholtz equation

u c 2 u = u 0 ε c 2 ε = ε 0 , (5)

where the source terms ( u 0 , ε 0 ) are the solutions of a corresponding boundary-value problem based on classicallinear elasticity, and the symbol 2 denotes as usual the Laplace operator 2 ( ) = ( ) , i i .

An argument for the stress σ in Equation (1) similar to that employed for the strain ε through Equation (2) can lead to the following expression between the gradient stress σ and the classical stress σ 0

σ c 2 σ = σ 0 , (6)

where it was assumed, for simplicity, that l ε 2 = l σ 2 c .

In applying the above approach to revisit dislocation mechanics, i.e. by using the Ru-Aifantis formalism recapitulating below:

· Gradient Constitutive Equation: σ = λ ( t r ε ) 1 + 2 G ε c 2 [ λ ( t r ε ) 1 + 2 G ε ] ,

· Ru-Aifantis Theorem: u c 2 u = u 0 { ε c 2 ε = ε 0 , σ c 2 σ = σ 0 ,

we obtain the following non-singular solutions for a screw dislocation

ε x z = b z 4 π [ y r 2 + y r c K 1 ( r c ) ] , ε y z = b z 4 π [ x r 2 x r c K 1 ( r c ) ] , (7)

where r denotes as usual the radial coordinate. It is readily noted that these expressions (similar ones hold for the stresses σ x z , σ y z ) give finite (zero) values at the dislocation line (see Figure 1) since for r 0 the modified Bessel function K 1 gives K 1 ( r / c ) c / r ( σ y z , ε y z ) 0 . It is also noted that the self-energy is given by the expression W s = G b z 2 4 π { γ + ln R 2 c } ( γ = 0.577 is Euler’s constant), such that there is no need for an ad hoc dislocation core r 0 assumption. This expression for the self-energy W s of the screw dislocation at hand ( b z is the Burgers vector) holds for an infinite cylinder ( R ) surrounding the dislocation line.

In calculating the self-energy of an edge dislocation W e surrounded by a finite

Figure 1. Distribution of strains/stresses associated with a screw dislocation. The first figure corresponds to the classical singular solution, while the second corresponds to the gradient nonsingular one.

cylinder of radius R, the corresponding expression for an edge dislocation is given by the more complex equation below [5] [12]

W e = b 2 4 π ( 1 ν ) { ln R 2 c + γ + 2 K 0 ( R c ) + 2 c R K 1 ( R c ) 2 c R 2 } , (8)

where b denotes again the Burgers vector, ν is the Poisson’s ratio, ( K 0 , K 1 ) denote the modified Bessel functions. By letting R , we have (as for the case of screw dislocation) the following limiting value W e = b 2 4 π ( 1 ν ) { ln R 2 c + γ + 1 2 } . A Plot of the aforementioned generalized finite expression for the self-energy of an edge dislocation is provided in Figure 2 in comparison with corresponding atomistic simulations based on a Stillinger-Weber potential [12] [13]. This comparison between the GradEla model and atomistic Stillinger-Weber calculations provide the following estimate [12] for the gradient coefficient c = l ε 2 = 0.2 - 2.2 Å 2 (the symbol Å denotes as usual Angstrom units). It is noted that as R all three models (classical elasticity, Stillinger-Weber and GradEla) converge. But as R approaches the dislocation line ( R 0 ) only the GradEla model goes smoothly to zero. The classical elasticity model holds up to distances bounded by the dislocation core, while the atomistic simulations provide results at smaller distances, but not up to the dislocation line.

Figure 2. Plots of the self-energy of an edge dislocation for classical elasticity (green), GradEla (red/blue) and Stillinger-Weber atomistic simulations (dots) for three edge partial dislocation configurations. See [12] for details.

3. London’s Quantum Mechanical Potential

Based on exact quantum mechanical calculations London [14] [15] has arrived at the following forms of the interatomic force F ( = d w / d r ) and interatomic potential w ( = w ( r ) )

F = d w d r ; w = w ( r ) = { 3 α 0 2 h v 4 ( 4 π ε 0 ) 2 1 r 6 = C r 6 ; r σ , ; r < σ , (9)

where C = 3 α 0 2 h v / 4 ( 4 π ε 0 ) 2 , α 0 is the atomic polarizability and ε 0 the vacuum dielectric permittivity. The quantities ( h , v ) denote respectively the Planck constant and the electron orbital frequency. The above analytical relation for the attractive interaction which holds until a critical distance σ , and there is no available a similar expression for the region where the interaction becomes repulsive going to infinity as r 0 . To describe quantitatively “repulsive” interactions for distances r < σ , Lennard-Jones [16] suggested the following modification of London’s potential

w ( r ) = A r 6 + B r 12 , (10)

where A and B are determined by fitting them to obtain through atomistic simulations the measured experimental values of macroscopic properties. Qualitative graphs for F and w ( r ) are provided in Figure 3 [17].

4. Gradient Modification of London’s Potential

Motivated by GradEla and a corresponding generalization of Newton’s gravitational (1/r) potential, we discuss below an analogous gradient modification of London’s quantum mechanical potential. It turns out [10] [11] that the gradient

Figure 3. Plots of the potential w ( r ) and the force F ( r ) for the Lennard-Jones potential [17].

enhanced London’s potential w L G is obtained in terms of its classical counterpart w through the inhomogeneous Helmholtz equation

( 1 l 2 2 ) w L G = w , w = w ( r ) = C r 6 . (11)

The solution of Equation (11) for ( w L G 0 , r ) is given by the expression

w L G ( r ) = A l e r / l r + C 48 l 6 { 4 l 4 r 4 + 2 l 2 r 2 + l r [ e r / l Ei ( r l ) e r / l Ei ( r l ) ] } , (12)

where A is a new integration constant, l is an internal length parameter, and Ei denotes the exponential integral Ei ( x ) = x e t t d t . Near the origin ( r 0 ), it turns out that w L G ( r ) C 12 l 2 r 4 , while at large distances ( r ) it approaches the classical London’s potential, i.e. w L G ( r ) C r 6 for r l . The qualitative behavior of London’s gradient potential, given by Equation (12) is provided in Figure 4(a). In Figure 4(b) the corresponding plots for the force F L G ( r ) ( = d w L G / d r ) are provided.

As an example application of the newly derived gradient potential, we consider the case of the Argon, for which parameter values of the non-gradient counterpart of the potential are available from computer simulations for liquid argon in accordance with experiment (see Table 6.1 of [17] and the data of [18] for the deduced numerical/experimental values). In particular, the ionization potential, designated as ε (in units of Joules or eV), as well as its location r m (in Å), can be estimated as ε = 1.95 × 10 21 J and r m = 0.37 nm respectively. The Lennard-Jones potential then can be uniquely determined from these parameters. For this purpose, Equation (10) is written in the form w L J ( r ) = ε ( ( r m / r ) 12 2 ( r m / r ) 6 ) , where it is evident that the minimum occurs

Figure 4. Qualitative plots of the London gradient potential w L G and the corresponding interaction force F L G . The scaling factors are C l 6 and C l 7 for w L G and F L G respectively.

at r m with w L J ( r m ) = ε and d w L J ( r m ) / d r = 0 . This point determines the transition from “attractive” to “repulsive” branch for distances r < r m . Additionally, the Lennard-Jones potential crosses zero at r = σ = 2 1 / 6 r m = 0.89 r m = 0. 32 nm . The parameters ( ε , r m ) are related with the ( A , B ) of Equation (10) through the relationship A = ε r m 12 = 4 ε σ 12 , B = 2 ε r m 6 = 4 ε σ 6 . The fitted London’s constant is C = 50 × 10 79 J m 6 , which was determined such as the classical London’s potential passes through the experimental potential minimum exactly at r m .

In order to demonstrate the ability of the gradient modification of London’s potential to recover the behavior of the Lennard-Jones potential for the Ar-Ar interaction case, we can adjust the gradient parameters ( A , l , C ) , such as the position of the potential minimum occurs at r m , i.e. w L G ( r m ) = ε , the corresponding potential curves are as close as possible (by minimizing their mean square error) for the attractive branch, and letting free the repulsive branch to approach infinity as r 0 , as dictated by the gradient part of the so modified potential. The obtained parameter values are A = 1.392 × 10 17 J , l = 0.57 Å, and C = 90 × 10 79 J m 6 (see also [12]). Figure 5 shows that the gradient modification of London’s potential fits quite well the Lennard-Jones potential, while both curves have their minima intersect at distance r m . It is noted that the gradient potential has the same asymptotic O ( r 6 ) distances r > r m , in agreement with the classical forms of both Lennard-Jones and London’s potential. It turns out, as expected, that the gradient modified London’s potential exhibits “repulsive”

Figure 5. Quantitative plots of the Ar-Ar fitting for the classical London, the Lennard-Jones and gradient London potential respectively.

branch for r < r m , where the change of slope occurs, in contrast to the classical form of London’s potential which exhibits only an attractive (1/r6) branch.

In a similar way other types of intermolecular or atomistic potentials, such as the Stillinger-Weber potential [19], can be generalized, a task currently being in progress by the Florida-Thessaloniki groups [11].

5. Discussion—Future Directions

A brief review of the robust gradient elasticity (GradEla) model was given first with emphasis on its ability to remove the undesirable dislocation singularities predicted by the classical Hookean elasticity. The removal of the singularity and the size of dislocation core are obtained by properly adjusting the internal length parameter multiplying the extra Laplacian term introduced in the classical elasticity stress-strain relation to account for heterogeneity effects within the elementary material volume considered. The GradEla robustness is due to the observation (Ru-Aifantis theorem) that explicit solutions of gradient elasticity can be obtained in terms of existing classical elasticity solutions through the use of an inhomogeneous Helmholtz equation for which the analytical mathematical results are available. This methodology is extended to revisit London’s quantum mechanical (attractive interaction) potential and generalize it to include both an “attractive” and a “repulsive” branch. By properly adjusting the internal length parameter, the behavior of a variety of empirical interatomic potentials, such as the Lennard-Jones and the Stillinger-Weber potentials, can qualitatively and quantitatively be readily recovered.

To further substantiate the usefulness of the approach specific materials should be considered in detail and the phenomenological parameter(s) of the so-derived gradient London’s potential should be determined for the material system at hand. In doing so, other possible generalizations of the gradient approach may also be pursued. In fact, two such generalizations are currently being considered: One is concerned with a fractional implementation of GradEla and another with the introduction of an additional biharmonic-like operator (or a bi-Laplacian) term to further generalize London’s potential. For the fractional implementation of GradEla, the reader can consult references [20] [21] [22] [23]. For the bi-Laplacian generalization of GradEla, the reader can consult [24] [25].

In concluding, we point out that the connection between average and local quantities through the Laplacian operator (as discussed first for the case of GradEla and adopted subsequently for the case of gradient London potential) was originally pointed out by Maxwell—the physical meaning of Laplacian. More details can be found in [26] where also combined gradient-stochastic models are discussed. The role of stochasticity to GradEla and intermolecular gradient potentials will also be a subject of future studies. In this connection, reference is made to a forthcoming review chapter [27], where a detailed account of generalized gradient interatomic potentials is provided, along with their extension to their fractional counterparts.


All authors acknowledge the support of the H2020-MSCA-RISE project, grant No. 734485 “Fracture Across Scales and Materials, Processes and Disciplines” (FRAMED), as well as the RISE project No. 824022 “Atomistic to Molecular to Bulk Turbulence” (ATM2BT). F. Shuang, B. Wang and P. Hu also gratefully acknowledge the support of the U.S. Department of Energy, Office of Basic Energy Sciences under the grant DE-SC0017715, which made this work possible under the supervision of K. E. Aifantis. The authors are grateful to KEA for suggesting and encouraging FS, BW and PH to work on the project, as well as to Professor Elias C. Aifantis who coordinated the whole effort (see also [27]). In fact, the AUTh-UF interaction was initiated during a visit of E.C. Aifantis to Gainsville, which was partially supported by the grants NSRF 2014-2020 for the project MIS 5005134 “Nano-chemomechanics in Deformation and Fracture: Theory and Applications in LiBs and SGS”, as welll as MIS 5045454 “Material Instabilities, Size Effects and Morphogenesis: Nanomaterials and Brain”.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.


[1] Aifantis, E.C. (1992) On the Role of Gradients in the Localization of Deformation and Fracture. International Journal of Engineering Science, 30, 1279-1299.
[2] Aifantis, E.C. (2011) On the Gradient Approach—Relation to Eringen’s Nonlocal Theory. International Journal of Engineering Science, 49, 1367-1377.
[3] Aifantis, E.C. (2016) Internal Length Gradient (ILG) Material Mechanics across Scales and Disciplines. Advances in Applied Mechanics, 49, 1-110.
[4] Ru, C.Q. and Aifantis, E.C. (1993) A Simple Approach to Solve Boundary-Value Problems in Gradient Elasticity. Acta Mechanica, 101, 59-68.
[5] Yu, M.G. and Aifantis, E.C. (1999) Dislocations and Disclinations in Gradient Elasticity. Physica Status Solidi (B), 214, 245-284.<245::AID-PSSB245>3.0.CO;2-P
[6] Askes, H. and Aifantis, E.C. (2011) Gradient Elasticity in Statics and Dynamics: An Overview of Formulations, Length Scale Identification Procedures, Finite Element Implementations and New Results. International Journal of Solids and Structures, 48, 1962-1990.
[7] Lazar, M., Maugin, G.A. and Aifantis, E.C. (2005) On Dislocations in a Special Class of Generalized Elasticity. Physica Status Solidi (B), 242, 2365-2390.
[8] Lazar, M., Maugin, G.A. and Aifantis, E.C. (2006) Dislocations in Second Strain Gradient Elasticity. International Journal of Solids and Structures, 43, 1787-1817.
[9] Aifantis, E.C. (2020) A Concise Review of Gradient Models in Mechanics and Physics. Frontiers in Physics, 7, 239.
[10] Aifantis, E.C. (2019) Invited Lecture in: Soft Matter Symposium 2019: Celebrating Five Years of the SMS, Organizers: T.E. Angelini, K.E. Aifantis and X. Tang, University of Florida, Gainesville, October 8-10.
[11] Aifantis, K.E. (2019) Current Work in Various Types of Gradient Enhanced Intermolecular and Atomistic Potentials Is on Progress by the Florida-Thessaloniki Groups.
[12] Kioseoglou, J., Dimitrakopulos, G.P., Komninou, Ph., Karakostas, T., Konstantopoulos, I., Avlonitis, M. and Aifantis, E.C. (2006) Analysis of Partial Dislocations in Wurtzite GaN Using Gradient Elasticity. Physica Status Solidi (A), 203, 2161-2166.
[13] Kioseoglou, J., Dimitrakopulos, G.P., Komninou, Ph., Karakostas, T. and Aifantis, E.C. (2008) Dislocation Core Investigation by Geometric Phase Analysis and the Dislocation Density Tensor. Journal of Physics D: Applied Physics, 41, Article ID: 035408.
[14] London, F. (1930) Zur Theorie und Systematik der Molekularkräfte. Zeitschrift für Physik, 63, 245-279.
[15] London, F. (1937) The General Theory of Molecular Forces. Transactions of the Faraday Society, 33, 8-26.
[16] Jones, J.E. (1924) On the Determination of Molecular Fields I. From the Variation of the Viscosity of a Gas with Temperature. Proceedings of the Royal Society of London A, 106, 441-462.
[17] Israelachvili, J.N. (2011) Intermolecular and Surface Forces. 3rd Edition, Academic Press, Cambridge.
[18] Parson, J.M., Siska, P.E. and Lee, Y.T. (1972) Intermolecular Potentials from Crossed-Beam Differential Elastic Scattering Measurements. IV. Ar + Ar. The Journal of Chemical Physics, 56, 1511-1516.
[19] Stillinger, F.H. and Weber, T.A. (1985) Computer Simulation of Local Order in Condensed Phases of Silicon. Physical Review B, 31, 5262-5271.
[20] Tarasov, V.E. and Aifantis, E.C. (2015) Non-Standard Extensions of Gradient Elasticity: Fractional Non-Locality, Memory and Fractality. Communications in Nonlinear Science and Numerical Simulation, 22, 197-227.
[21] Aifantis, E.C. (2019) Fractional Generalizations of Gradient Mechanics. In: Tarasov, V.E., Ed., Handbook of Fractional Calculus with Applications, Volume 4, De Gruyter, Berlin, 241-262.
[22] Tarasov, V.E. and Aifantis, E.C. (2019) On Fractional and Fractal Formulations of Gradient Linear and Nonlinear Elasticity. Acta Mechanica, 230, 2043-2070.
[23] Parisis, K., Konstantopoulos, I. and Aifantis, E.C. (2018) Nonsingular Solutions of GradEla Models for Dislocations: An Extension to Fractional GradEla. Journal of Micromechanics and Molecular Physics, 3, Article ID: 1840013.
[24] Lazar, M., Maugin, G.A. and Aifantis, E.C. (2006) On the Theory of Nonlocal Elasticity of Bi-Helmholtz Type and Some Applications. International Journal of Solids and Structures, 43, 1404-1421.
[25] Aifantis, E.C. (2009) Non-Singular Dislocation Fields. IOP Conference Series: Materials Science and Engineering, 3, Article ID: 0712026.
[26] Aifantis, E.C. (2014) Gradient Material Mechanics: Perspectives and Prospects. Acta Mechanica, 225, 999-1012.
[27] Aifantis, E.C. (2020) Gradient Extension of Classical Material Models: From Nuclear & Condensed Matter Scales to Earth & Cosmological Scales. In: Ghavanloo, E., Fazelzadeh, S.A. and Marotti de Sciarra, F., Eds., Size-Dependent Continuum Mechanics Approaches: Theory and Applications, Springer Tracts in Mechanical Engineering, Springer, in press.

Copyright © 2024 by authors and Scientific Research Publishing Inc.

Creative Commons License

This work and the related PDF file are licensed under a Creative Commons Attribution 4.0 International License.