Acta Polytechnica CTU Proceedings


doi:10.14311/APP.2018.15.0057
Acta Polytechnica CTU Proceedings 15:57–62, 2018 © Czech Technical University in Prague, 2018

available online at http://ojs.cvut.cz/ojs/index.php/app

MOLECULAR STATICS SIMULATION OF NANOINDENTATION
USING ADAPTIVE QUASICONTINUUM METHOD

Karel Mikeša,∗, Ondřej Rokošb, Ron H. J. Peerlingsb

a Czech Technical University in Prague, Thákurova 7, 166 29 Prague 6, Czech Republic
b Eindhoven University of Technology, P.O. Box 513, 5600 MB, Eindhoven, The Netherlands
∗ corresponding author: Mikes.Karel.1@fsv.cvut.cz

Abstract. In this work, molecular statics is used to model a nanoindentation test on a two-dimensional
hexagonal lattice. To this end, the QuasiContinuum (QC) method with adaptive propagation of the fully
resolved domain is used to reduce the computational cost required by the full atomistic model. Three
different adaptive mesh refinement criteria are introduced and tested, based on: (i) the Zienkiewicz–
Zhu criterion (used for the deformation gradient), (ii) local atoms’ site energy, and (iii) local lattice
disregistry. Accuracy and efficiency of individual refinement schemes are compared against the full
atomistic model and obtained results are discussed.

Keywords: Molecular statics, quasicontinuum method, adaptivity, nanoindentation.

1. Introduction
Nanoindentation is a commonly used testing proce-
dure applied to small volumes of materials for mea-
suring their micromechanical properties. Typically, a
hard tip (i.e. indenter) with known mechanical prop-
erties is pressed into an examined sample of unknown
mechanical properties. Loading force and penetration
depth of the indenter are recorded during the load-
ing and unloading stages, providing a basis for the
estimation of the unknown mechanical properties.
Numerical models are typically used as a tool for

better understanding the underlying phenomena, and
to obtain detailed information about local mechanisms
occurring below the indenter tip (such as dislocation
nucleation, propagation, and interaction), which di-
rectly influence measured reaction force. To this end,
both the indenter and specimen are typically mod-
elled at the atomistic level using molecular statics
or molecular dynamics, entailing high computational
costs when realistic configurations and dimensions are
used. The QuasiContinuum (QC) method (cf. e.g. [1])
is employed to simplify the full atomistic model, to re-
duce the associated computational costs, and to allow
for modelling of realistic situations.
This paper focuses on the predictive abilities of

an adaptive QC methodology (recalled in Section 3)
in combination with three types of error indica-
tors/estimators for local mesh refinement compared
against the full atomistic simulations. In particular,
(i) the Zienkiewicz–Zhu error estimator (used for the
deformation gradient), (ii) an indicator based on local
atoms’ site energy, and (iii) an estimator based on
local disregistry profiles are tested for a simple two-
dimensional indentation test. The individual defini-
tions are outlined in Section 3.3, whereas the accuracy
and associated computational costs are discussed in
Section 4.

The paper closes with conclusions and recommen-
dations in Section 5.

2. Full atomistic model
Atomistic models based on molecular statics are char-
acterized by an underlying lattice in combination
with an interatomic potential. In this work, a two-
dimensional hexagonal lattice with lattice spacing d0 is
used, as shown in Fig. 1. Individual atom interactions
are described by the Lennard–Jones (LJ) potential,
defined as

φαβ(rαβ) = ε

[( rm

rαβ

)12
− 2

( rm

rαβ

)6
]
, (1)

where rαβ = ‖rβ−rα‖`2 denotes the distance between
two atoms α and β, rm denotes the distance at which
the interaction energy reaches its minimum, and ε is
the energy well depth.
The total potential energy associated with the en-

tire atomic structure is computed as a sum over all
interactions, i.e.

E(r) = 1
2

NAtm∑
α,β; α 6=β

φαβ(rαβ), (2)

where NAtm represents the number of atoms, and r
is a column storing their positions.

Because evaluation of the interatomic potential for
all pairwise combinations is computationally expen-
sive, and because long-distance interactions have neg-
ligible contributions to the total potential energy, a
cut-off radius rcut is considered [1], beyond which inter-
actions are neglected. Such a simplification introduces
a discontinuity of φαβ at rcut, which is removed by
subtracting a linear function to assure zero value and
zero slope of φαβ at rcut. As sketched in Fig. 1, a

57

http://dx.doi.org/10.14311/APP.2018.15.0057
http://ojs.cvut.cz/ojs/index.php/app


K. Mikeš, O. Rokoš, R. H. J. Peerlings Acta Polytechnica CTU Proceedings

d0

d0 θ = 120◦

rcut = 2.5d0

h0 = d0
√
3/2

Figure 1. A geometry of hexagonal lattice and corre-
sponding cut-off radius (dashed line).

cut-off radius rcut = 2.5 d0 is employed to provide
next-to-nearest interactions.

In order to find a stress-free configuration, an initial
relaxation is carried out on an ideal periodic lattice
with spacing d0 = rm, which results in a reduced lat-
tice spacing d0 = 0.9917496 rm used for constructing
the initial system.
The geometry of employed indentation test is

sketched in Fig. 2. The specimen domain is of the size
128d0 × 128h0, contains 16, 862 atoms, and considers
atoms near the bottom and both vertical edges as
fixed, whereas the top edge is a free surface. The
flat indenter is modelled at the atomistic level using
the same hexagonal lattice as used for the specimen,
but having infinite stiffness. Its geometry is specified
through a width 11d0 at the tip and two surfaces in-
clined by 60◦, as shown in Fig. 3 (left). The positions
of all indenter atoms are prescribed in 80 uniform
loading, and 80 uniform unloading steps, achieving
the maximum indentation depth 8d0. The interaction
strength between atoms of the indenter and the tested
material is reduced by a factor of 0.55 (compared
to the atoms of the tested material) to prevent tear-
ing of the indented specimen during the unloading
stage. The total potential energy of the atomistic
system E(r) is minimized at each time step using the
trust-region algorithm; for further details see e.g. [2].

3. Quasicontinuum Method
The Quasicontinuum (QC) method is a concurrent
multiscale technique introduced in [3]. The key idea
consists in combining the accurate but expensive atom-
istic description only in regions of high interest with
a cheap continuum approximation elsewhere. The
specimen domain therefore is divided into two parts:
(i) the fully-resolved region, in which the full non-
local atomistic model is used, seamlessly coupled with
(ii) a coarse-grained continuum region, in which inter-
polation through triangular elements along with an
efficient summation scheme is introduced.

x
y

12
8h

0

128d0

Figure 2. A sketch of the tested setup: specimen
being indented (grey), and indenter (black).

Figure 3. Detail of indenter area at the begin-
ning of the loading process (left), and after unloading
(right).

3.1. Interpolation
The first step in the QC reduction is interpolation,
which introduces the so-called repatoms through which
the kinematic behaviour of the entire system is recon-
structed according to

r = Φrrep, (3)

where rrep is a column storing positions of all
repatoms, interpolated through an interpolation ma-
trix Φ associated with the adopted triangulation. The
number of repatoms is typically much smaller com-
pared to the number of all atoms, reducing thus the
computational effort required. Two example trian-
gulations associated with different sizes of the fully-
resolved regions are shown in Fig. 4.

3.2. Summation
In the second QC reduction step, a so-called summa-
tion rule is introduced to avoid the necessity of visiting
all atoms when assembling the total potential energy
in Eq. (2). To this end, the site energies of all atoms
situated inside a triangular element are approximated
by the energy of only a few, or even one sampling
atom and its corresponding weight factor wα, i.e.

E =
NSampAtm∑

α

wαφα. (4)

58



vol. 15/2018 Molecular Statics Simulation of Nanoindentation

0 20 40 60 80 100 120 140

-60

-40

-20

0

20

40

60

0 20 40 60 80 100 120 140

-60

-40

-20

0

20

40

60

Figure 4. Initial triangulation with small (left)
and large (right) fully-resolved region. Repatoms are
shown as black dots, interpolation elements as blue
triangles, sampling atoms as red dots, and remaining
atoms as grey dots.

In Eq. (4), φα is the site energy of a sampling atom
α, defined as

φα = 1
2

NAtm∑
β; α6=β
rαβ<rcut

φαβ(rαβ). (5)

Although multiple possibilities are available in the
literature to select sampling atoms and their corre-
sponding weight factors, the central summation rule
introduced in [4] is used hereafter (cf. Fig. 4).

3.3. Adaptivity
The area of high interest (i.e. the fully-resolved region
and hence also the associated triangulation) can adap-
tively evolve at each time increment or iteration to
accommodate dislocation movements, while retaining
QC efficiency.

To this end, the fully-resolved region is sequentially
updated as follows. At each time increment, the sys-
tem is equilibrated for a fixed fully-resolved region.
Using a selected refinement criterion (detailed below),
the triangulation is checked and refined if required.
Any additional atoms are added as repatoms, the inter-
polation mesh is updated, and equilibrium is restored.
Such a procedure is repeated until the mesh refine-
ment criterion is satisfied for all elements, proceeding
subsequently to a new load increment.
In total three refinement criteria are tested:
(i) The Zienkiewicz–Zhu error estimator (ZZ), as
introduced in [5], used for the deformation gradient.
That is, projection of the deformation gradient is
used to estimate the local error inside each element.
A threshold value ZZtr is specified to determine
elements to be refined, and atoms located in those
elements to become repatoms.

(ii) The error indicator based on a local energy crite-
rion uses the local atoms’ site energy to determine
elements to be refined. First, the site energy of each
sampling atom is tested for the following condition

φα ≥ Etr, (6)

where Etr is a selected critical energy threshold.
Sampling atoms satisfying the condition of Eq. (6),

as well as their neighbours within radius rref , are
labelled as critical atoms. Interpolation elements
that contain at least one critical atom are refined, i.e.
all atoms in such elements are added as repatoms.

(iii) The error indicator based on disregistry profile
works in a similar way as the local energy criterion.
The only difference consists in the condition for the
selection of critical atoms, now specified as

NAtm∑
β; α6=β
rαβ<rcut

||rβ − rαt ||`2 ≥ Dtr, (7)

where Dtr is a selected disregistry threshold, and
||rβ − rαt ||`2 denotes the distance of atom β from
the closest position of its theoretical neighbour-
ing atom α considered in the reference configura-
tion. The positions of either all eighteen nearest-
and next-to-nearest-neighbour atoms or just the six
nearest-neighbour atoms can be used to evaluate
the condition of Eq. (7). All atoms within radius
rref are again labelled as critical atoms.

4. Results
Numerical simulations of the nanoindentation test,
as described in Section 2, are carried out to analyse
the adaptive propagation of the fully-refined region
associated with movements of individual dislocations.
In the full atomistic model, initially all atoms are

displaced elastically until a critical penetration depth
is reached, which triggers nucleation of four disloca-
tions positioned symmetrically (due to the symmetry
of the problem) under both indenter edges. Upon
further loading, these dislocations propagate along
preferred lattice directions until reaching the bound-
ary of the specimen. During unloading, however, when
the last symmetric dislocation pair annihilates, one di-
rection is preferred due to numerical round-off errors,
resulting in a nonsymmetric final shape of the indent
(see Fig. 3 (right)). Results of the full simulation are
used as the reference solution against which the QC
simulations are compared, and are summarized along
with the corresponding performance in Table 1. The
force–displacement diagram is shown in Fig. 6.
In QC simulations, movements of individual dis-

locations can be captured properly only inside the
fully-resolved region. The adaptive algorithms intro-
duced in Section 3.3 are used to this end, providing
results summarized in Fig. 5 for the Zienkiewicz–Zhu
error estimator and loading steps number 28 (pene-
tration of the first symmetric dislocation pair outside
of the original area of high interest), 39 (one step
before nucleation of the second symmetric disloca-
tion pair), 55 (one step before the nucleation of the
third symmetric dislocation pair), and 80 (maximum
penetration depth). The following cases for error es-
timators/indicators in QC methodology have been
tested:

59



K. Mikeš, O. Rokoš, R. H. J. Peerlings Acta Polytechnica CTU Proceedings

Figure 5. Results corresponding to QC simulations with Zienkiewicz–Zhu mesh refinement criterion. Repatoms in
different loading steps are shown, corresponding to: loading step number 28 (top left), 39 (top right), 55 (bottom
left), and 80 (bottom right). The colour of individual atoms corresponds to the local disregistry level. Dashed lines
indicate the reference specimen shape.

Initial/final Iterations Mesh Time
repatoms iterations

Full solution 16862/16862 4980 - 1.00 (1:26:31)

QC fix small 952/952 3549 - 0.061

QF fix large 1900/1900 3560 - 0.099

QC zz 952/6224 8957 242 0.772

QC dis 952/8143 7536 149 0.751

Table 1. Numerical performance corresponding to individual computational models.

60



vol. 15/2018 Molecular Statics Simulation of Nanoindentation

Figure 6. Normalized force–displacement diagrams
corresponding to the different computational models.

Figure 7. The number of repatoms relative to the
number of all atoms as a function of load increment
for the different computational models.

(1) two simulations with fixed fully-resolved regions
of different sizes, see Fig. 4, referred to as QC small
and QC large;

(2) one simulation with an adaptive area of high
interest using the ZZ criterion with ZZtr = 0.015,
referred to as QC zz;

(3) one simulation with an adaptive area of high
interest based on local energy with refinement ra-
dius rref = 5d0 and threshold value corresponding
to 6% change of the initial potential energy (for
a typical internal atom this value corresponds to
Etr = −2.95); and

(4) one simulation with adaptive area of high interest
based on lattice disregistry with threshold value
Dtr = 0.5d0, the nearest-neighbour option, and
refinement radius rref = 5d0, referred to as QC dis.
The simplest QC simulation with a small fixed fully-

resolved area (QC small) provides a significant speed-
up of a factor of 16 compared to the full atomistic
model. The initial elastic behaviour and nucleation
of the first symmetric dislocation pair is captured
accurately, cf. Fig. 6, whereas in later stages (i.e.
once the first symmetric dislocation pair reaches the
boundary of the fully-resolved region), the resulting

behaviour becomes overly stiff. This is clearly visible
in the corresponding force–displacement diagram.
A partial improvement can be achieved by enlarg-

ing the fully-resolved QC region (QC large), which,
nevertheless, suffers from the same shortcomings. The
overestimation of the force–displacement diagram is
less significant, capturing more accurately the unload-
ing branch. The shapes of local force peaks do not
correspond, however, to the exact solution due to
the obstructed dislocation movements. The achieved
speed-up is approximately of the order of 10.

Clearly, only the adaptive schemes are able to cap-
ture properly the dislocation movement under the
tip of the indenter. Note that all adaptive QC ex-
amples are initialized with the fully-resolved area of
the same size as used for the QC small simulation.
The ZZ criterion detects and refines elements near
individual dislocations, allowing for their propagation
throughout the entire specimen. This is reflected in
the corresponding force-displacement curve in Fig. 6,
where a significant improvement in accuracy compared
to the approaches with fixed meshes can be observed.

Finally, the local lattice deflection criterion is used,
which provides (in the case considered) practically
the same results as the energy criterion (not shown).
Generally, both criteria are not equivalent and dif-
ferences can be observed in examples with different
type of atoms, where the energy criterion is more
robust at the interface. However, in case of a sin-
gle crystal with one type of atom, for each value of
disregistry threshold a corresponding value of energy
threshold that provides virtually the same results can
be found. Compared to the ZZ condition, local lattice
deflection provides similar accuracy in the resulting
reaction force, while requiring a slightly higher num-
ber of repatoms (cf. Fig. 7). In spite of this fact, its
overall performance is slightly faster compared to ZZ
(cf. Tab. 1), explained by the fewer number of mesh
iterations needed. Unlike ZZ, the indicator based on
the local disregistry provides two parameters, Dtr and
rref , for controlling the mesh evolution. Whereas Dtr
locates the position inside the dislocation core region
(see colours in Fig. 5), rref reflects the size of the
fully-resolved region considered around that position.
Large values of rref therefore provide fast propagation
in the direction of the dislocation movement, but at
the same time may require too many repatoms in the
perpendicular direction. The value rref = 5d0 proved
to provide a good balance.
The performance of all QC approaches is summa-

rized in Tab. 1, where the number of repatoms, re-
quired number of solver and mesh iterations, and
computational times are reported. All adaptive QC
approaches provide accurate descriptions during the
loading as well as unloading stages, resulting in ac-
ceptable estimates of the corresponding reaction force.
Although the obtained gain in computation time is
relatively small (about 25% compared to the full sim-
ulation), a more substantial speed-up is expected for

61



K. Mikeš, O. Rokoš, R. H. J. Peerlings Acta Polytechnica CTU Proceedings

simulations with larger domains.

5. Conclusions
In this contribution, a QuasiContinuum (QC) method
with different mesh refinement strategies has been
compared against the full molecular statics simulation
for the case of a two-dimensional nanoindentation test.
It has been shown that all three tested adaptive QC
approaches provide acceptable agreement with the
underlying full atomistic simulation.
Two adaptive approaches (based on local energy

and lattice disregistry) provide equivalent results com-
pared to the widely used Zienkiewicz–Zhu error esti-
mator, while performing slightly faster due to fewer
mesh iterations. Overall, adaptive QC approaches
reduce computational time by approximately 25%
compared to the full simulation. For large domains,
more significant speed-ups with similar accuracy can
be expected.

List of symbols
d0 horizontal lattice spacing
h0 vertical lattice spacing
φαβ Lennard–Jones (LJ) potential
ε energy parameter of LJ potential
rm distance parameter of LJ potential
rαβ distance between two atoms α and β
rcut cut-off radius
E potential energy
NAtm number of all atoms
NSampAtm number of sampling atoms
r atom positions
rrep repatom positions
Φ interpolation matrix
φα site energy of an atom α

wα weight factor
Etr energy threshold
Dtr disregistry threshold

Acknowledgements
Financial support of this research from the Grant
Agency of the Czech Technical University in Prague
(SGS18/037/OHK1/1T/11), and from the Czech Science
Foundation (GAČR 17-04150J) is gratefully acknowledged.

References
[1] E. B. Tadmor, R. E. Miller. Modeling materials:

continuum, atomistic and multiscale techniques.
Cambridge University Press, 2011.

[2] A. R. Conn, N. I. Gould, P. L. Toint. Trust region
methods, vol. 1. Siam, 2000.

[3] E. B. Tadmor, M. Ortiz, R. Phillips. Quasicontinuum
analysis of defects in solids. Philosophical magazine A
73(6):1529–1563, 1996.

[4] L. A. A. Beex, R. H. J. Peerlings, M. G. D. Geers.
Central summation in the quasicontinuum method.
Journal of the Mechanics and Physics of Solids
70:242–261, 2014.

[5] O. C. Zienkiewicz, J. Z. Zhu. The superconvergent
patch recovery and a posteriori error estimates. part 2:
Error estimates and adaptivity. International Journal
for Numerical Methods in Engineering 33(7):1365–1382,
1992.

62


	Acta Polytechnica CTU Proceedings 15:57–62, 2018
	1 Introduction
	2 Full atomistic model
	3 Quasicontinuum Method
	3.1 Interpolation
	3.2 Summation
	3.3 Adaptivity

	4 Results
	5 Conclusions
	List of symbols
	Acknowledgements
	References