Abstract
Disordered solids are ubiquitous in engineering and everyday use. Although research has made considerable progress in the last decades, our understanding of the mechanics of these materials is, at best, in an embryonic state. Since the nature of disorder complicates the realization of physically meaningful continuum-mechanical models, particle-based molecular descriptions provide a powerful alternative. This paper reviews the numerical realization of classical molecular dynamics from an engineer’s perspective, starting with selecting potential functions, boundary conditions, time integration, and thermodynamic ensembles. Then, we discuss the concept of the potential energy landscape and the computational realization of the most suitable minimization methods. Subsequently, we discuss the algorithms necessary to numerically generate disordered materials, considering their thermodynamic properties and structural identification. We comprehensively and critically review computational methods and strategies available to mimic disordered materials on a molecular level and discuss some intriguing phenomena that are, to date, mostly ignored when applying models based on continuum-mechanical frameworks. We present the crucial difference between the shear response of a crystalline and a disordered structure. In this context, we elaborate on why it is beneficial to use an overdamped, athermal description to disentangle the complex deformation mechanics of disordered solids and comprehensively discuss the theory of the mechanics of disordered materials, including the problems of prediction and reversibility. Furthermore, we examine the fracture process on the nanoscale and investigate the response behavior to more complex deformation protocols. Finally, we provide critical conclusions, including challenges and future perspectives for engineers.
Similar content being viewed by others
1 Introduction
Window glass, golf club heads, toothpaste, gravel, soap bubble rafts, and the earth’s crust are typical examples of disordered solids. Intriguingly, the general idea of such a class of solids has been around for millennia. For example, the discovery of ancient glassy artifacts in Egypt shows that these materials were already recognized by early civilizations around 12000 BC [1]. However, it was only a few decades ago when scientists acknowledged their high potential for engineering applications and design and recognized the necessity to study their unique, rich and fascinating dynamical properties.
General mechanical phenomena associated with disordered solids, such as the fracture process, occur on scales that range over ten orders of magnitude [2], making both the understanding and the computational modeling challenging. For example, the complex dynamics of earthquake faults [3] are linked very intimately with the nanoscale mechanics of disordered materials and their elastoplastic behavior on the nanoscale [4, 5] since the phenomenon of disorder inherently governs the complex responses of both these systems. Out of the extensive collection of disordered solids, we pay attention to the material behavior on the nanoscale in this paper. However, we note that some mechanical phenomena reviewed and simulated here also appear on larger scales.
As disordered materials continuously gain research attention, we briefly summarize important review papers in the field and how our paper positions itself within the broad body of literature. One field of application that exhibits the problem of disorder constitutes the investigation, development, and design of bulk metallic glasses. Several review articles have been published during the last fifteen years. Schuh et al. [6], Trexler et al. [7], and Egami et al. [8] reviewed the deformation and fracture mechanics of metallic glasses from a theoretical and experimental point of view, also with an emphasis on bridging to macroscopic scales. Chen [9] summarized and discussed the structure-property relationship of metallic glasses, with particular emphasis on the strength and ductility of the material on the microscopic scale, while Wang [10] reviewed the elastic properties of bulk metallic glasses, with particular emphasis on experimental studies. Furthermore, Greer et al. [11] presented an overview of shear banding in metallic glasses. Later, Hufnagel et al. [12] published a follow-up review of the previous paper by Schuh et al. [6], presenting an overview of theoretical and experimental research progress in bulk metallic glasses. Thus, this paper does not present the structure-property relationship of bulk metallic glasses and quantitative studies regarding their elastic properties, toughness, ductility, etc.
Falk and Langer [14] published an overview of the shear transformation theory, a mean-field theory for disordered solids, which will, therefore, not be the focus of this paper. Another review article was published by Rodney et al. [15], who presented a survey of modeling disordered materials on multiple scales. Furthermore, Nicolas et al. [4] and Tanguy [5] reviewed the elasto-plastic behavior of disordered materials.
Molecular modeling of disordered solids on the nanoscale requires a toolbox of non-classical techniques that often go beyond classical molecular dynamics. A detailed overview of the relevant computational strategies is missing from the literature. Thus, this paper reviews the latest developments in state-of-the-art computational strategies that simulate the mechanics and thermodynamics of disordered solids using a molecular mechanical framework. Furthermore, we attached importance to reviewing and elaborating on the rich dynamical and thermodynamical phenomena of disordered materials under mechanical deformation. The paper may also serve as a step-by-step instruction for implementing a research tool that includes the state-of-the-art approaches for the mechanical and thermodynamical investigation of glasses and amorphous solids, with a view to stimulating research and building macroscopic, continuum-mechanical material models. In this context, we obtained all numerical results of this paper using our in-house code implemented in Julia [16] to cover all the possible major obstacles that may arise when diving into the topic of disordered solids.
1.1 Outline
In Sects. 2, 3, 4, and 5, we discuss the essential minimum, dealing with the theoretical background and implementation necessary to model disordered structures on the nanoscale. Thus, these sections summarize the basic features of our in-house tool to perform state-of-the-art computational modeling of disordered materials using molecular mechanics.
We start our review with a relatively short introduction to modeling particles on an atomic scale in Sect. 2. This way, the mathematical essentials necessary to model the structure of disordered materials are presented. In Sect. 2.1, we introduce selected pair-potential functions mainly used for densely packed disordered systems and multi-body potentials primarily used for materials that reveal individual network structures. Furthermore, in Sect. 3, we introduce the implementation of boundary conditions in rectangular cells, and subsequently, we present the realization of boundary conditions in tilted cell geometries. Finally, in Sect. 4, we discuss the kinematics of remote deformation and the mechanical loading of tilted cell geometries. We introduce the affine deformation mapping from the reference to current cell geometries, including several strain measures, which one may interpret as the equivalent of the deformation gradient and corresponding strain definitions from continuum mechanics.
Section 5 introduces the general concept of the thermodynamics of disordered solids. To be able to extract the response of the material, especially the part of the response that deviates from the affine deformation introduced in Sect. 4, we discuss microscopic measures, i.e., atomic displacements, atomic deformation gradients, and atomic strains in Sect. 5.1. Furthermore, Sect. 5.2 discusses thermodynamic measures, i.e., kinetic energy, temperature, and stress. Subsequently, Sect. 5.3 presents the most popular time integration scheme used in molecular dynamics to date, represented by the Verlet algorithm. Finally, in Sects. 5.4, 5.5, 5.6, we review a selection of essential dynamical formulations, i.e., the microcanonical ensemble, the canonical ensemble, the isobaric-isothermal ensemble, and their numerical realization using the Verlet algorithm.
In Sect. 6, we present a selection of general numerical algorithms that significantly increase the efficiency of the methods described in the preceding sections. In particular, we introduce the concept of Verlet lists by defining a cutoff radius and present the implementation of linked cell lists.
In Sect. 7, we introduce the concept of the potential energy landscape of disordered systems and the meaning of inherent structures. We review essential minimization algorithms to generate inherent structures, i.e., the steepest descent algorithm in Sect. 7.1, and the conjugate gradient algorithm in Sect. 7.2, including a detailed numerical description of implementing these algorithms in a molecular mechanical framework.
In Sect. 8, we review the computational strategies for generating disordered systems. This approach includes the melting-quenching technique in Sect. 8.1, to obtain glasses in a thermodynamically consistent way, and Monte Carlo bond-switching procedures in Sect. 8.2, to generate two-dimensional glasses and amorphous network materials athermally. We generate two categories of disordered materials using these two techniques. The first category is a dense random-close-packing model material as a benchmark system that should feature the main properties of metallic glass. The second category is a continuous random network model that should feature the main properties of oxide glasses. We numerically generate three examples of disordered solids. Firstly, a Lennard-Jones glass is generated as an example of a random-close-packing model using the melting-quenching approach. This model material has proven to feature some main properties of metallic glasses and has been extensively applied in the literature. Secondly, bulk silica glass samples are generated using the melting-quenching approach and a Stillinger-Weber type potential function. Thirdly, a two-dimensional model network material is presented based on the topological information of a real two-dimensional silica polymorph extracted from experimental data in the literature and prepared with the appropriate relaxation processes using an appropriate pair potential. To generate this material, we perform a Monte Carlo bond-switching algorithm, introducing a sequence of topological flip transformations.
Section 9 presents a comprehensive overview of the deformation mechanics of disordered solids using the athermal quasistatic deformation method and related computational techniques. After introducing the general idea of the athermal quasistatic deformation method, an over-damped zero temperature procedure to model the deformation behavior of disordered solids at low temperatures, we discuss the linearization of the potential energy landscape around a particular reference configuration in Sect. 9.1, introducing the Hessian and its numerical implementation. Subsequently, in Sect. 9.2, we discuss the deformation behavior of solids subjected to athermal quasistatic simple shear. In particular, we demonstrate the fundamental difference in the shear response of an ordered, i.e., a hexagonally packed system in Sect. 9.2.1, and the shear response of a binary Lennard-Jones glass system in Sect. 9.2.2. We present a detailed theoretical framework for describing the non-affinity of the elastic response in Sect. 9.2.3 and the approach to catastrophic inelastic events in Sect. 9.2.4. Subsequently, the transition into cascading events that may be interpreted as nano shear bands is outlined in Sect. 9.2.5. Furthermore, we demonstrate the picture of elementary shear events in network glasses in Sect. 9.2.6. Thereafter, in Sect. 9.3, we critically examine approaches to predict inelastic events in disordered systems and their reversibility upon un- and reverse-loading in Sect. 9.4. We complete this part of the paper by discussing deformation protocols other than simple shear in Sect. 9.5, and provide some mechanical interpretations considering the relevant literature. In this context, in Sect. 9.5.1, uniaxial and simple tensile protocols are presented for densely packed systems and network glasses. In Sect. 9.5.2, we discuss the intriguing mechanical behavior of network glasses under pressure and shear deformation under pressure, presenting a substantial collection of numerical examples.
In Sect. 10, we discuss the statistical classification of the mechanical response of disordered solids and establish the connection with the concept of self-organized criticality.
Finally, in Sect. 11, we present the conclusions and recommendations for possible directions of future research.
1.2 Definitions of Ordered and Disordered States
Before diving into molecular modeling of disordered materials, we summarize a few definitions and distinctions of the types of solids that appear in this paper. At first, we note that we draw on the, to our knowledge, latest definitions by Zanotto and Mauro [17], although we also note here that scholars do not entirely agree on this matter [20,21,22]. Fortunately, the particular topics of disagreement do not interfere with the numerical approaches reviewed in this paper.
Crystals (crystalline solids, ordered materials) are highly ordered structures, as they can be obtained by periodically multiplying an elementary unit cell in space. However, following the laws of thermodynamics, a structure is never perfectly ordered since defects, vacancies, or grain boundaries are always present [23]. Increasing the number of defects in a crystalline structure, the solid may transform into a state whereby it exhibits such a high number of defects that the concept of a regular cell structure must be abandoned, and the material must be categorized as disordered or non-crystalline. Overall, disordered solids can be divided into two groups: glasses and amorphous solids.
Glasses are non-crystalline solids (disordered materials) that exhibit a glass transition. One possibility to generate glasses is to cool down a melt fast enough so that its structure does not have enough time to crystallize. During the cooling process, the structure visits the state of a supercooled liquid that exists between the melting and the glass transition temperature. The atomic structure of a glass is similar to that of its supercooled liquid. Glass is a metastable state of matter that exists on human timescales, while its ultimate fate is to solidify, i.e., crystallize at geological time scales [17].
Amorphous solids are non-crystalline solids (disordered materials) that are not glasses. Amorphous derives from the Greek word “\(\alpha \mu o\rho \phi o\varsigma\)”, which means shapeless. They do not exhibit a glass transition temperature and cannot be generated by quenching a melt [24]. We note here that it is common practice in the literature to use the notation “amorphous solid” as the generic term for non-crystalline solids, cf. [4, 25]. However, most of these works mainly investigate purely mechanical deformation, performed at temperatures, far below the glass transition temperature. In those cases, these differences turn out to be somewhat secondary.
We summarize those notations, initially defined by Gupta [24] and Zanotto and Mauro [17], by dividing a solid into ordered (crystalline) and disordered (non-crystalline) solids. Disordered solids may further be subdivided into amorphous solids and glasses, whereas this classification can be seen as mutually exclusive, i.e., an amorphous solid is never a glass and vice-versa.
2 Molecular Mechanical Description of Particles and Their Interaction
This section briefly discusses the mechanical interaction of molecular particles confined to a finite-sized triclinic cell geometry. Parts of this section are seen as an introduction that allows non-experts to dive into the computational modeling of disordered materials without extensively studying additional literature, such as Leimkuhler and Matthews [26] or Lelièvre and Stoltz [27].
Therefore, numerical methods and algorithms that are relatively straightforward in classical molecular mechanics will be discussed in a nutshell. In contrast, selected strategies may be discussed in more detail if crucial for the problems presented later.
2.1 Inter-atomic Potentials and Forces
In order to tackle the mechanics of molecular systems, we start by localizing N particles or atoms in a three-dimensional space and defining their positions as:
where every vector \(\varvec{r}_i\) denotes a distance vector from the fixed origin of a Cartesian coordinate system to the atomic position. Every atom can be displaced in all three directions, leading to a total number of 3N degrees of freedom. Within the framework of particular methods and algorithms used to describe molecular systems, it may sometimes be advantageous to use an alternative description of the atomic configuration introducing a concatenated position vector, written as:
All remaining quantities of interest, depending on position, such as velocities, accelerations, and forces, will be denoted in terms of these two representations in this paper.
One quantifies the interactions between all atoms via the definition of the total potential energy of the system, which is a scalar quantity that depends only on the atomic positions. One evaluates the total potential energy of the configuration by the sum of 1- to n-body terms via:
where \(m,j,k=1,\dots ,N\). The first term considers the influence of an external field that affects the particles, e.g., through the introduction of box walls. This term is neglected for the demonstrations in this review, as mainly periodic boundary conditions are applied. The second term in Eq. (3) is the two-body term considering the interaction of all possible pairs of atoms in the configuration, leading to a sum of \(N(N-1)\) possible pair potentials \(\mathcal {U}\left( \varvec{r}_m,\varvec{r}_j \right)\). The condition \(j>m\) in the two-body contribution of Eq. (3) ensures that the pairs in the configuration are not double counted. The two-body potential term constitutes the most crucial contribution in molecular mechanics so that, in many cases, all other higher-order terms are neglected. We will, however, also discuss models of disordered network materials, where three-body contributions are relevant, as presented in the third term in Eq. (3). Equivalently to the second term, one must ensure that every interaction triplet \(\mathcal {U}\left( \varvec{r}_m,\varvec{r}_j,\varvec{r}_k \right)\) occurs only once in the summation, introducing the conditions \(j>m\) and by \(k>j\). Otherwise, one would count every triplet six times during the summation of the three-body contribution. Since four-body terms are somewhat irrelevant for describing disordered solids, we neglect all contributions with more than three members.
The force \(\varvec{f}_i\) acting on the ith atom is defined as the negative gradient of the potential energy with respect to its atomic position:
whereas the gradient operator with respect to the ith atomic position is defined as:
Equivalently, one defines the concatenated force vector \(\varvec{f}\) acting on the whole configuration as the negative gradient of the potential energy with respect to the concatenated position coordinate \(\varvec{r}\) as:
whereas the gradient operator is defined as:
where all three directions of every particle lead to a system with 3N degrees of freedom.
Considering Eq. (3) and neglecting external forces on the system and all inter-particle contributions involving more than three atoms, one obtains the force on the ith particle as:
where the first term is the force resulting from the two-body contribution, and the second term is the force resulting from the three-body contribution of the potential.
2.1.1 Two-Body Contributions
In molecular mechanics, one assumes that the potential between two atoms is invariant with respect to translation and rotation. Therefore, one concludes that the potential between two atoms depends only on the interatomic absolute distance:
We note that the pair potential includes the collective information about the atom types but not the information of which atom number within the pair belongs to which atom type. Under these circumstances, one can write:
The two-body contribution of the force on atom i in Eq. (8) is not equal to zero if and only if \(m=i\) or \(j=i\). Consequently, the force on atom i due to the two-body contributions results in:
From the result in Eq. (11), every particle interacts with every other particle in terms of pair-interaction forces. As a consequence, the resulting two-body force \(\varvec{f}_i^{(2)}\) on the ith atom is equal to the sum of its interaction forces \(\varvec{f}_{ij}^{(2)}\) with every other atom j in the ensemble except itself:
For illustrative purposes of the interatomic forces associated with the contribution of the two-body potential, we present a dynamic system of four particles in a three-dimensional reference system.
Illustration of all two-body interactions of a system of four atoms in a fixed three-dimensional Cartesian coordinate system
Indicated by a color code and their size, one can see that the particles may differ regarding their mass and their general physical/mechanical properties. Every ith atom with mass \(m_i\) is subjected to a force vector \(\varvec{f}_i^{(2)}\). Defining a particular amount of mass for every particle and assuming that the consideration of only two-body contributions is sufficient to provide a physically meaningful description of a particular solid, one may apply the momentum principle in differential form, i.e., Newton’s second law of motion, for each of the atoms separately:
which leads to the simplest description of a molecular system.
Taking advantage of the symmetry properties in Eqs. (9) and (10) we define the interatomic distance vector as \(\varvec{r}_{ij}=\varvec{r}_i-\varvec{r}_j\) and its respective magnitude as \(r_{ij}=\left\| \varvec{r}_{ij} \right\|\). Following these definitions, one writes:
with
Considering the relationships from Eq. (16) one can merge Eqs. (14) and (15) to the essential property of pair potentials, written as:
This property is also illustrated in Fig. 1. It is commonly known as Newton’s third law of motion, roughly summarized as to every action there is always an opposed reaction. Thus, every pair potential \(\mathcal {U}_2\left( \varvec{r}_i,\varvec{r}_j \right)\) results in a force pair \(\varvec{f}_{ij}^{(2)}\) and \(-\varvec{f}_{ij}^{(2)}\). From a numeric implementation point of view one sets all forces to zero and, then, performs a double loop with the counters i and j over all atoms i and all partner atoms \(j>i\), and adds the force \(\varvec{f}_{ij}^{(2)}\) to the force vector \(\varvec{f}_i^{(2)}\) on atom i, while subtracting the force \(\varvec{f}_{ij}^{(2)}\) from the force vector \(\varvec{f}_{j}^{(2)}\).
Considering the findings from Eqs. (14)–(16), one writes the interatomic force due to one pair potential as:
where \(\nicefrac {\varvec{r}_{ij}}{r_{ij}}\) defines the unit direction from atom j to i; while the derivative of the pair potential with respect to the radial distance \(\mathcal {U}_2'(r_{ij})\) leads to the scalar value of the interatomic pair force contribution.
We summarize that the direction of a two-body force contribution is defined by the distance vector from atom j to i, and the corresponding magnitude is defined by the negative derivative of the pair potential with respect to the inter-atomic distance \(r_{ij}\).
Inter-atomic force pairs due to a the pair potential and b the parts \(h^{(1)}\left( \varvec{r}_{ij},\varvec{r}_{ik} \right)\), \(h^{(2)}\left( \varvec{r}_{jk},\varvec{r}_{ji} \right)\), \(h^{(3)}\left( \varvec{r}_{ki},\varvec{r}_{kj} \right)\) of the three-body potential
2.1.2 Three-Body Contributions
We now pay attention to the third term in Eq. (8) and define all possible triplets of atoms that are—equivalently to the two-body contribution—symmetrical under permutation of the atom indices and invariant under translation and rotation. In doing so, the objective is to stabilize particular bond angles within the atomic triplets in order to enable a more realistic description of particular materials, such as silicon, which is, besides oxygen, one of the indispensable ingredients of silica glass and, therefore, a crucial element in some implementations covered in this paper. The force \(\varvec{f}_i^{(3)}\) on atom i resulting from the three body contributions, as presented in term two of Eq. (8), is only not equal to zero if and only if \(i=m\) or \(i=j\) or \(i=k\). Taking advantage of the property of permutation of the atom indices, the three-body force on atom i is rewritten as:
Thus, atom i experiences the force \(\varvec{f}_i^{(3)}\) as a result of all possible three-body terms of which it is a part.
Taking advantage of the invariance with respect to translation and rotation every contribution \(\mathcal {U}_3\left( \varvec{r}_i, \varvec{r}_j, \varvec{r}_k \right)\) involving three atoms may be subdivided, so that the contribution is rewritten as the sum of three triplet terms (Fig. 2):
For the sake of compactness in the formulations to come, we denote \(h^{(1)}:=h^{(1)}\left( \varvec{r}_{ij},\varvec{r}_{ik} \right)\), \(h^{(2)}:=h^{(2)}\left( \varvec{r}_{jk},\varvec{r}_{ji} \right)\), and \(h^{(3)}:=h^{(3)}\left( \varvec{r}_{ki},\varvec{r}_{kj} \right)\).
From an implementation point of view, one starts by setting all atomic forces that originate from the three-body contributions to zero. Considering Eq. (20), one concludes that, in order to compute the force that acts on atom i due to one three-body contribution \(\mathcal {U}_{3}\left( \varvec{r}_i,\varvec{r}_j,\varvec{r}_k\right)\), one must compute the gradient of all three triplets \(h^{(1)}\), \(h^{(2)}\), and \(h^{(3)}\) with respect to the atomic position \(\varvec{r}_i\) and repeat this for every atom i in the configuration, as also indicated by the sum in the last line of Eq. (19). However, by repeating this procedure for every atom in the configuration, one realizes that a significant portion of the gradients is evaluated multiple times. Thus, to avoid repetitive evaluations of the gradients, one considers only the contributions for which \(i>j>k\) applies. To identify the gradients ignored through this assumption, we investigate the relation of the gradients in one triplet with respect to the atomic positions. For instance, we present the gradient of the triplet \(h^{(1)}\) with respect to atom i as:
Considering that
one can rewrite the gradient of \(h^{(1)}\) with respect to atom i as:
In Eq. (22), \(\varvec{I}_3\) denotes the three-dimensional identity matrix. It can easily be shown that the property in Eq. (23) can be equivalently applied to the two remaining triplets \(h^{(2)}\) and \(h^{(3)}\), which leads to the more general statement:
where \(\alpha =1,2,3\). The mechanical interpretation of this crucial finding may be seen as the well-known static force equilibrium that applies at all time instants in every three-body triplet, written as:
where \(\varvec{f}_{\kappa }^{(3,\alpha )}= -\nabla _{\varvec{r}_{\kappa }}h^{(\alpha )}\). Since every triplet is in a state of force equilibrium, one calculates only the two gradients at the right-hand side of Eq. (24) and uses them to evaluate the third gradient. Thus, this approach allows one to save \(\nicefrac {1}{3}\) of the total evaluation of all three-body contributions in the configuration.
Since the three-body terms are invariant with respect to translation and rotation, one can redefine the triplets as functions dependent on the interatomic distances and the enclosed bond angle as:
where we denote \(\tilde{h}^{(1)}:=\tilde{h}^{(1)}\left( r_{ij},r_{ik},\vartheta _{jik}\right)\), \(\tilde{h}^{(2)}:=\tilde{h}^{(2)}\left( r_{jk},r_{ji},\vartheta _{kji}\right)\), and \(\tilde{h}^{(3)}:=\tilde{h}^{(3)}\left( r_{ki},r_{kj},\vartheta _{ikj}\right)\) for the sake of compactness.
Physically meaningful molecular simulations live and die by choosing an appropriate potential function. Although the literature is constantly enriched by modified and improved potentials, we qualitatively and quantitatively present only a few potential functions that are particularly relevant to studying the mechanical behavior of disordered materials. The following subsection presents a model potential for densely packed disordered systems, such as metallic glasses. Furthermore, we present two potential types for bulk network glasses, such as silica glass, and a Yukawa-type potential for model network materials, i.e., systems that correspond with Zachariasen’s random network model in two dimensions, in Appendix 10.
2.2 A Lennard-Jones Potential for Binary 2D Glasses
One may imagine the interaction of two spherical particles by (i) a prominent resistance to compression; therefore, leading to a strong repelling force at small distances, and by (ii) some distance range of attraction so that particles may seek their equilibrium and, therefore, their force-free configuration. The first potential function that fulfills such conditions so that they can physically/mechanically be described in a rather simple way is the Lennard-Jones (LJ) potential [28, 29]. This pair potential is defined as a high-order \(\left( \frac{1}{r} \right)\) rational function:
where the term of order 12 is responsible for the repelling part, while the lower term of order 6 is responsible for the attractive part. Physically interpreted, the first term models the non-bonded overlap between the electron clouds, while the second part models the van der Waals interaction due to electron correlations [30]. The solid blue plot in Fig. 3 presents the LJ potential function, while the red solid plot shows its derivative, i.e., the principal ingredient for the interatomic forces, as discussed earlier in this section.
Function of the Lennard-Jones potential and its derivative
Dimensionless Lennard-Jones units are applied in most papers dealing with such potential types. Distances are then measured in \(\sigma\), which is \(r_{ij}\) at the zero crossing of the potential. The potential energy is measured in \(\epsilon\), which is the absolute value of the minimum of the potential, and masses are measured in Dalton [30]. The time is measured as: \(T^* = \sqrt{\nicefrac {m\sigma ^2}{\epsilon }}\). Both parameters \(\sigma\) and \(\epsilon\) depend on the combination of types of the respective atom pair. Numerous numerical studies in the literature indicate that a particular proportion of larger (heavier) and smaller (lighter) particles are considered, to enable a better understanding of the dynamics of disordered systems, as will be discussed later in this paper. In the case of binary materials using atomic types \(\bar{\alpha }\) and \(\bar{\beta }\), three combinations must be provided for each of the parameters \(\epsilon\) and \(\sigma\) to model all possible interactions: \(\bar{\alpha }-\bar{\alpha }\), \(\bar{\beta }-\bar{\beta }\), \(\bar{\alpha }-\bar{\beta }\).
We will employ this potential to model densely packed disordered systems as a qualitative two-dimensional mechanical model for the structure of a metallic glass. As mentioned in Sect. 1, quantitative studies for metallic glasses are not reviewed in this paper, since there is already a wide variety of literature available.
Additionally, we present relevant potential functions describing the behavior of network glasses implemented for the numerical investigations in this paper in Appendix 10.
3 Boundary Conditions
Most studies investigating the mechanical behavior of disordered solids apply periodic boundary conditions [30]. Using this framework, every point in the sample can be considered a center point. Thus, for investigations of mechanical phenomena, such models excel compared to those that include boundary walls or pulling areas since they generally prevent unwanted localization effects.
Applying periodic boundary conditions, the sample is periodically duplicated in all three directions in a three-dimensional space. Consequently, the distance calculation would include not only the atoms in the original cell but also the atoms from all surrounding cells duplicated, which would be an infinite number of pair interactions. Consequently, in the case of short-range potentials, one uses the concept of the minimum image convention [31], illustrated in Fig. 4. In this figure, a rectangular original box geometry O is surrounded by its periodic duplicates, denoted as A–H, arranged clockwise. The general idea is that every atom in the original box O is located at the center of its region, which has the same size and shape as the geometry of the box. In this illustration, a quadratic box area is translated so that atom 4 is located at its center, presented by the dotted lines. Instead of calculating the distances with all remaining infinite number atoms, one considers only the minimum distance with neighbor atoms in the dotted box geometry. In this particular example, atom 4, highlighted in blue, has 5 neighbor atoms, i.e., 3, 5, \(1_\text {A}\), \(2_\text {B}\), and \(6_\text {H}\). In this manner, every atom still has \(N-1\) neighbors. It is noteworthy that long-range interactions, such as electrostatic forces whose pair-potentials decrease proportionally to \(\frac{1}{n}\), must additionally be considered when incorporating periodic boundary conditions using strategies such as the Ewald method, the tree-code approach, or the fast-multipole method. For more information regarding the treatment of long-range interactions in periodic cells, we refer the reader to the relevant literature, cf. [30, 32,33,34,35].
Periodic boundary conditions and minimum image criterion
For the practical realization of boundary conditions, let us first consider a rectangular cell geometry containing a particular atomic configuration for the purpose of illustration. Suppose the projection of the interatomic distance-vector \(\varvec{r}_{ij}\) in one basis direction is larger than half the box length in this direction. In that case, the component of this distance vector is replaced by its periodic duplicate, considering positive or negative distance directions. A two-dimensional illustration of periodicity in \(x_1\)-direction is illustrated in Fig. 5. The two remaining directions are equivalent, while all combinations of the three directions of periodic interactions are possible.
Two-dimensional illustration of a system with periodic boundary conditions in \(x_1\)-direction in a rectangular cell geometry
However, cell geometries are, in general, not rectangular. Nevertheless, we have to ensure periodicity for our simulations, as initially proposed by Lees and Edwards [36]. We provide a two-dimensional example of such a tilted box in Fig. 6. The box or simulation cell is defined as:
where \(\varvec{h}_1,\varvec{h}_2,\varvec{h}_3 \in \mathbb {R}^3\) denote the bravais vectors defining the cell boundaries [37]. The bravais tensor \(\varvec{H} \in \mathbb {R}^{3\times 3}\) describes, by construction, a parallelepiped. The first vector \(\varvec{h}_1\) is collinear with the \(x_1\)-axis, and the plane spanned by \(\varvec{h}_1\) and \(\varvec{h}_2\) remains in the \(x_1\)-\(x_2\)-plane [38]. Obviously, the equivalent two-dimensional cell is defined by the two bravais vectors \(\varvec{h}_1\) and \(\varvec{h}_2\), where \(\varvec{h}_1\) is collinear with the \(x_1\) axis. To realize periodic boundary conditions in tilted cells, the system is periodically duplicated over the bravais vectors, as shown in Fig. 6. A convenient way of dealing with periodic parallelepiped cells is to take advantage of the transformation of an atomic position \(\varvec{r}_i\) in \(\varvec{H}\) into an artificial position \(\varvec{\textsf{r}}_i\) in a unit cell, whose basis consists of the three unit vectors of the Cartesian coordinate system. As firstly proposed by Ray and Rahman [39], this transformation reads as follows:
so that the detection of the distance vectors of the periodic duplicates is performed in a unit box \(\varvec{I}_3\) and is transformed back into the physical configuration \(\varvec{H}\). Thus, the practical realization of periodic boundary conditions remains equivalent to the rectangular case when considering the transformation in Eq. (31).
Furthermore, one must consider the effect that atoms may leave the box boundaries and move into their periodic duplicate, which results in the atoms entering the simulation box on the opposite side. Implementation-wise, this is not an issue as long as they do not leave the first duplicated box. However, it makes sense to let all atoms that leave the box geometrically enter at the corresponding opposite side.
Two-dimensional illustrations of a system with periodic boundary conditions of a tilted cell geometry
The introduction of periodic boundary conditions has a crucial influence on the potential energy of the system since the particles interact over their boundaries with their periodic duplicates. Following this line of thought, we must consider these contributions; however, one must also consider the additional information of the shape of the periodic cell. Thus, we extend the definition of the potential energy in a periodic simulation cell to:
4 Affine Deformation Kinematics
Before macroscopically deforming a system, one starts with defining the geometry of a reference cell, written as:
where, equivalently to Eq. (30), one denotes the bravais vectors in the reference configuration of the simulation cell as \(\hat{\varvec{h}}_1, \hat{\varvec{h}}_2, \hat{\varvec{h}}_3 \in \mathbb {R}^3\) so that \(\varvec{H} \in \mathbb {R}^{3\times 3}\). In this reference configuration, every particle is assigned to its reference position \(\hat{\varvec{r}}_{i0}\). External deformation is superimposed to the configuration by altering the bravais vectors of the simulation cell and, consequently, changing the boundary conditions from the reference configuration \(\hat{\varvec{H}}\) to the current configuration \(\varvec{H}\). This approach includes the geometric correction of the box boundaries and the definition of the corresponding affine displacement field that must be subjected to the entire atomic configuration, as qualitatively illustrated by the two-dimensional sketch in Fig. 7.
Qualitative illustration of an affine deformation of a triclinic cell geometry in two dimensions
Considering Eq. (31), every particle can be mapped to its position in the unit cell in the reference configuration by \(\varvec{\textsf{r}}_i=\hat{\varvec{H}}^{-1}\hat{\varvec{r}}_i\) and in the current configuration by \(\varvec{\textsf{r}}_i=\varvec{H}^{-1}\varvec{r}_i\). Following this line of reasoning, every particle can be mapped from the reference configuration \(\hat{\varvec{H}}\) to the current configuration \(\varvec{H}\) [39, 40]:
This mapping may be seen as the discrete equivalent to the deformation gradient from continuum-mechanical formulations [41]. Thus, we follow this notation and introduce the deformation gradient as a linear mapping from the reference to the current cell geometry as:
which is a tensor of order two. We note that a rigorous comparison with the continuum-mechanically motivated deformation gradient is challenging, since the models dealt with in this paper are discrete by definition. Nevertheless, where possible, we emphasize the connections between the discrete and continuum descriptions.
The mappings discussed in this section are purely affine since they are used for subjecting external deformation to the deformation cell. Depending on the complexity of the material, a configuration may respond in such a manner that it deviates from the purely affine deformation. In this case, the mapping of the displaced particles back into the reference configuration via \(\hat{\varvec{r}}_i = \varvec{F}^{-1}\varvec{r}_i\) leads to an atomic configuration in the reference cell \(\hat{\varvec{H}}\), where the position vectors \(\hat{\varvec{r}}_i\) may have varied from their initial positions. The mechanics of such processes and the extraction of the portion that deviates from affine deformation, i.e., the so-called non-affine deformation, is discussed in Sect. 9.
Although one may evaluate a respective remote strain measure directly using the deformation gradient, we chose to perform a small exercise to extract the remote strain by following the change in length of a vector that points from one coordinate in the cell to a neighboring point in the cell. This way, we follow a distance vector located in the simulation cell during deformation. We will take advantage of this additional exercise later in Sect. 5.1, where we will discuss the local, non-affine strain response of an atomic structure. In the reference and the current configuration, this distance vector pointing from a neighboring to a center point is defined as:
The position vectors \(\hat{\varvec{s}}\) and \(\hat{\varvec{s}}'\), as well as \(\varvec{s}\) and \(\varvec{s}'\), denote a center and a neighboring point in the reference and the current configuration, respectively. Inserting the mapping from the reference to the current configuration in terms of the deformation gradient above, one can easily derive the relation:
which reminds us that the affine mapping of distance vectors between two positions in the simulation cell is invariant with respect to the choice of these positions and is, therefore, constant throughout the deformation cell. One measures the change in the quadratic length of the distance vector \(\varvec{d}\) during deformation from the reference to the current configuration. One obtains:
The term \(\varvec{E} = \frac{1}{2} \left( \varvec{C} - \varvec{I}_3 \right)\) is referred to as the Green–Lagrangian (or also Green-St. Venant) strain tensor. The tensor \(\varvec{C}=\varvec{F}^T\varvec{F}\) is the right Cauchy-Green deformation tensor (\(\varvec{F}\) is on the right), also referred to as Green deformation tensor. Both strain measures are second-order tensors. Looking at Eq. (38), it is obvious that \(\varvec{E}\) and \(\varvec{C}\) are both associated with the distance vector \(\hat{\varvec{d}}\), leading to strain measures that are associated with the reference configuration. We again make use of the change in the quadratic length of the distance vector. However, this time, we use the inverse mapping \(\hat{\varvec{d}}=\varvec{F}^{-1}\varvec{d}\), so that one obtains:
The term \(\varvec{A} = \frac{1}{2} \left( \varvec{B}^{-1} - \varvec{I}_3 \right)\) is referred to as the Euler-Almansi strain tensor. The tensor \(\varvec{B}=\varvec{F}\varvec{F}^T\) is referred to as the left Cauchy-Green deformation tensor (\(\varvec{F}\) is on the left side), also referred to as the Finger deformation tensor. Looking at Eq. (39), it is obvious that both tensors \(\varvec{A}\) and \(\varvec{B}\) are associated with the distance vector \(\varvec{d}\), leading to strain measures that are associated with the current configuration. For further information regarding strain tensors and their relationships, we refer to Holzapfel [41]. We further note that applying logarithmic strain formulations, such as the Hencky strain, can be advantageous in particular situations. Regarding the Hencky strain, we refer to the relevant literature [42,43,44]. This paper will use the Green–Lagrangian framework to extract the remote strains during mechanical deformation.
All deformation tensors above depend only on the deformation gradient \(\varvec{F}\). Therefore, in the case of a macroscopic box deformation state that is performed purely affine, the corresponding macroscopic remote strain state is constant and quantified by one strain tensor only. At this point, we note again that molecular systems may not respond in an affine manner, leading to complex local strain responses, as discussed later in this paper.
When performing molecular mechanics simulations, the macroscopic strain constitutes one essential tensorial object for quantifying the structure-property relationships of complex materials. One incrementally alters the strain state by manipulating the cell bravais vectors and observes the structural response, e.g., measuring the potential energy or the corresponding stress tensor often as average quantities evaluated over all particles in the entire cell. It becomes evident that displacement- or strain-controlled protocols are convenient and mainly accepted as the standard approach to investigate molecular mechanical investigations of disordered materials since softening effects and fracture processes can be represented straightforwardly. However, we note that force-controlled approaches have also been successfully applied to study the behavior of disordered materials [45].
We illustrate three two-dimensional box deformation examples by manipulating the set of bravais vectors, and present the corresponding affine displacement fields in Fig. 8.
Three two-dimensional examples of cell deformation types starting from a reference cell subjected to an initial shear angle \(\hat{\gamma }:=\hat{\gamma }_{12}\) and the corresponding affine displacement fields; a stretch of the box by stretching the \(\varvec{h}_1\) vector; b stretch of the box by stretching the \(\varvec{h}_2\) vector; c engineering shear \(\gamma :=\gamma _{12}\) by manipulating the \(x_1\)-component of the \(\varvec{h}_2\) vector
As a first deformation example, we show the effect of the elongation of the first bravais vector \(\hat{\varvec{h}}_1\) by multiplying it by a factor of \(1+\alpha\), as shown on the left side of Fig. 8a. To remain general, we set the two-dimensional reference cell, defined by the bravais letters \(\hat{\varvec{h}}_1\) and \(\hat{\varvec{h}}_2\), to a configuration that is initially sheared by a reference tilt angle \(\hat{\gamma }:=\hat{\gamma }_{12}\), as depicted in Fig. 8.
As a second deformation example, Fig. 8b presents the equivalent elongation of the reference cell by multiplying the second bravais vector \(\hat{\varvec{h}}_2\), by a factor of \(1+\beta\).
As a third deformation example, we present the engineering (simple) shear strain deformation mode in Fig. 8c. When performing simple shear, one fixes the bottom of the cell while horizontally deforming the top of the box. Manipulating the second bravais vector \(\hat{\varvec{h}}_2\) simple shear is realized by adding a value of \(\gamma \hat{h}_{2,2}\) to the component \(\hat{h}_{2,1}\), where we quantify the shear angle as \(\gamma :=\gamma _{12}\).
The three illustrations on the right side of Fig. 8 reveal the affine displacement fields corresponding to the respective bravais vector manipulations. The affine deformation fields constitute a linear displacement interpolation between the cell boundaries defined by the bravais vectors. Thus, the strain corresponding to the affine displacement field is constant throughout the cell.
A combination of different deformation modes allows for an infinite variety of possible deformation protocols, such as volumetric compression and combinations of shear and volumetric deformation for the respective studies of the structure-property relationships of complex materials. However, as also stated in the literature, although it can be easily realized, it is shown that simple shear deformation is not that “simple” when it comes to the response of the material [46]. The deformation gradient of simple shear is obtained by the superposition of the pure shear contribution, i.e., the deformation of the material that results in a measurable stress reaction and the contribution of pure rotation. Consequently, in two dimensions, one writes:
Figure 9a shows the tilting deformation of a configuration box by \(\gamma\) during simple shear.
The kinematics of shear deformation; a simple shear is a combination of pure shear and rotation; b pure shear deformation on a rectangular box due to a perpendicular push-pull combination
The square and inscribed circle become a parallelogram with an inscribed ellipse. Thus, the material experiences a push-pull combination that increases nonlinearly with increasing shear angle \(\gamma\), highlighted in blue in Fig. 9a. Below, we investigate the kinematics of simple shear in more detail. In order to find the principal strain states of simple shear and the corresponding directions, one solves the eigenvalue problem of the Green–Lagrangian strain tensor of simple shear. However, we note that this procedure may be performed with the help of any other strain measure discussed above. The Green–Lagrangian strain tensor of simple shear in two dimensions is written as:
which is, by definition, symmetric and positive definite, if \(\gamma\) is non-zero. Solving the eigenvalue problem \((\varvec{E}-\varvec{I}_3\varepsilon ) \varvec{\lambda } = \varvec{0}\), one obtains two eigenvalues in two dimensions. These are referred to as the principal strains \(\varepsilon _1\) and \(\varepsilon _2\), and one obtains two corresponding eigenvectors \(\varvec{\lambda }_1\) and \(\varvec{\lambda }_2\), defining the directions of the principal strains, highlighted in red in Fig. 9a. Solving the eigenvalue problem, the Green–Lagrangian strain tensor may be decomposed as:
where the principal strains are evaluated as:
The transformation matrix \(\varvec{T}=[\varvec{\lambda }_1,\varvec{\lambda }_2]\) is assembled by the orthonormal eigenvectors belonging to the corresponding eigenvalues \(\varepsilon _1\) and \(\varepsilon _2\), and it transforms the strain state from the Cartesian coordinate system \(\varvec{e}_i\) into the principal strain coordinates \(\varvec{\lambda }_i\). Since we rotate from the coordinate system \(\varvec{e}_i\) into the principal strain state, one may assemble the rotation matrix not only by the eigenvectors depending on the shear angle \(\gamma\), but also by the orthogonal rotation depending on an angle \(\alpha\) defined between \(\varvec{e}_1\) and \(\varvec{\lambda }_1\) or equivalently between \(\varvec{e}_2\) and \(\varvec{\lambda }_2\):
Using simple trigonometric relations, the relationship between the incline angle \(\alpha\) and the tilting angle \(\gamma\), as visualized in Fig. 9a, can be conveniently defined by the relationship [44]:
When performing molecular simulations, it can be convenient to induce a pure shear deformation protocol into the directions of \(\varvec{e}_1\) and \(\varvec{e}_2\) by subjecting the material to the push-pull combinations and preserving a constant cell volume, as presented in Fig. 9b.
In a purely elastic framework, simple and pure shear can be compared to this kinematical framework. However, we emphasize that a transformation from one to the other shear strain state is not feasible as soon as irreversible effects occur during the shear deformation protocol.
We note that configurations can also be subjected to mechanical deformation by using pulling layers. Using this type of boundary condition, particular atom groups are frozen and subjected to specific displacement histories. This strategy is rather unpopular for investigating the mechanical behavior of disordered materials because of unwanted local effects near the pulling areas. However, it can be beneficial when investigating the influence of the loading direction on particular anisotropic materials [47]. In this case, it must be ensured that the local area of investigation is far enough away from the pulling area to minimize its local influence.
We furthermore note that the complexity of the material response does not necessarily originate only from the type of the spacial deformation picture but also the deformation history, such as loading- and unloading curves. The high complexity and unpredictability of such problems, particularly related to the mechanics of disordered materials, are discussed in more detail in Sect. 9.3.
5 Molecular Dynamics
This section summarizes the essential minimum information required to understand and implement classical molecular dynamics simulations. We refer to the relevant literature for more detailed information on statistical mechanics and classical molecular dynamics simulations [30, 31, 48]. We will use the formulations presented in this section, to prepare glasses using the melting quenching method, which will be introduced later in Sect. 8.1.
The major goal of molecular dynamics is to evaluate the response trajectories of the particles while trying to restrain the system to particular geometric and thermodynamic boundaries. To this end, we first summarize some important geometric and thermodynamic measures in Sects. 5.1 and 5.2, before we proceed with the numeric response evaluation in Sect. 5.3 and outline some important ensembles in molecular dynamics in Sects. 5.4, 5.5 and 5.6.
5.1 Local Geometric Measures
During a molecular simulation, every particle travels on an arbitrarily complex motion curve which is generally confined to the cell \(\varvec{H}\) of the configuration that may also change its geometry with time. Defining a reference cell \(\hat{\varvec{H}}\) and the corresponding reference configuration \(\hat{\varvec{r}}_i\) at a distinct molecular state, we introduce the first geometric response measure, which is the time-dependent particle displacement vector, defined as:
where it is possible to concatenate the atomic displacements to a vector \(\varvec{u} \in \mathbb {R}^{3N}\). Thus, the displacement field describes the total deformation response of every particle in the configuration at any time. It must be noted that the periodic boundary conditions must be taken into account when applying Eq. (47).
In Sect. 3, we introduced the boundary conditions of a molecular cell and Sect. 4 how to deform the whole cell by manipulating the three Bravais vectors that define the cell boundaries, introducing one global deformation gradient \(\varvec{F}\). However, as also mentioned in Sect. 4, the system does, in general, not respond in a purely affine manner, so that a portion of the displacements responds non-affinely even if thermal vibrations are neglected. This property is illustratively visualized in Fig. 10. This deviation from the affine response is always true for disordered systems. However, for non-centrosymmetric crystalline structures, the structural response also deviates from the affine deformation [49,50,51,52].
Qualitative illustration of a non-affine material response due to remote cell deformation
Consequently, one may quantify local deformation gradients that deviate from the macroscopic remote deformation gradient \(\varvec{F}\). Following this thought, we introduce a local deformation gradient allocated to every atomic position. One way of performing that is to extract the movement of every center particle i with respect to all its surrounding “neighboring” particles. Since this objective attempts to include continuum mechanical ideas in the discrete molecular framework, we apply the definition of the distance vector from a center to a neighboring point in Eq. (36). In this way, \(\hat{\varvec{s}}\) and \(\varvec{s}\) now point to a center particle i, denoted by the position vectors \(\hat{\varvec{r}}_i\) and \(\varvec{r}_i\) in the reference and current configuration. Equivalently, \(\hat{\varvec{s}}'\) and \(\varvec{s}'\) now point to a neighboring particle, denoted by the positions \(\hat{\varvec{r}}_j\) and \(\varvec{r}_j\) in the reference and current configuration. The local deformation gradient is supposed to simultaneously map all distance vectors from the reference configuration \(\hat{\varvec{d}}_{ij}:=\hat{\varvec{r}}_i-\hat{\varvec{r}}_j\) to the current configuration \(\varvec{d}_{ij}:=\varvec{r}_i-\varvec{r}_j\).
In order to evaluate the local deformation gradient, one assembles a matrix \(\hat{\varvec{Y}}\), containing all atomic distance vectors in the reference configuration, written as:
and a matrix \(\varvec{Y}\), containing all atomic distance vectors in the current configuration, written as:
One aims to find a local linear deformation mapping \(\varvec{F}_i\) such that \(\mathcal {J}\left( \varvec{F}_i\right) = \left\| \varvec{Y} - \hat{\varvec{Y}}\varvec{F}_i \right\| ^2\) is minimized. The first order optimality condition for the minimum of the convex functional \(\mathcal {J}\left( \varvec{F}_i\right)\) is provided by:
meaning that the local deformation gradient can be conveniently evaluated by:
It is worth mentioning that the periodic boundary conditions must be considered when evaluating the inter-particle distance vectors. The average of all local deformation gradients over all particles in the cell corresponds to the macroscopic deformation gradient \(\varvec{F}\), which can directly be evaluated by the bravais vectors of the simulation cell. From an implementation perspective, one may not consider all possible pairs in the configuration but introduce a radial cutoff value around the center atom i, so that the respective neighbor atom is only considered if it is situated within this cutoff distance. More details regarding the cutoff-radius are provided in Sect. 4.
Local strain measures, discussed in, such as the local Green–Lagrangian strain tensor \(\varvec{E}_{i}\), can directly be evaluated from the local deformation gradient \(\varvec{F}_i\) as follows:
Other local strain definitions can be evaluated equivalently using \(\varvec{F}_i\).
It is worth noting that Scagnetti et al. [53] did not use the deformation gradient but the displacement gradient as a basis for measuring the local strain. In doing so, they numerically differentiated the trajectories of the atomic displacement vectors to evaluate the Green–Lagrangian strain tensor.
Evaluating the local (atomic) deformation gradients \(\varvec{F}_i\) and the corresponding local (atomic) strain tensors \(\varvec{E}_i\) for all atoms (\(i=1\dots N\)) provides discrete tensorial fields describing the deformation kinematics of the current atomic configuration confined to the cell in \(\varvec{H}\), compared to the reference atomic configuration confined to the cell \(\hat{\varvec{H}}\). However, extracting such local geometric measures requires knowledge of atomic positions in the current configuration. Thus, one must apply an appropriate solution algorithm to evaluate the response from the reference to the current state. In that regard, selected solution algorithms are discussed in Sects. 5.3 and 9.
5.2 Thermodynamic Measures: Kinetic Energy, Temperature, Stress
We start with the definition of the kinetic energy of an ensemble [31], which is, on the one hand, the vibration energy of all particles:
Considering the generalized equipartition theorem [31], the kinetic energy is, on the other hand, uniformly distributed to every particle:
In this equation, d denotes the dimension of the problem. In our case, we have \(d=2\) or \(d=3\). Eq. (54) states that every degree of freedom in the system is associated with an average energy of \(\frac{1}{2}k_B T\) kgm\(^{2}\)s\(^{-2}\). The factor \(k_B=1.38066\times 10^{-23}\) kgm\(^{2}\)s\(^{-2}\)K\(^{-1}\) denotes the Boltzmann constant. In the case of dimensionless units, \(k_B\) is set equal to one, which allows one to circumvent tedious conversions into appropriate units and concentrate solely on the dynamics of the system. The expression in Eq. (53) is equated with Eq. (54), leading to the measure of the temperature, which is a positive quantity that is proportional to the vibrational energy of the system:
As the particles vibrate nonlinearly in the ensemble at any temperature greater than zero, the time evolution of the temperature fluctuates and is, particularly for smaller systems, a noisy function. However, with an increasing number of particles, this noise becomes smaller and averages out.
The stress evaluation in microscopic atomic scale mechanics may also be seen as an attempt to connect with the formulations of continuum mechanics. It originates from the formulation of the virial, and considers both a term that relates to the kinetic energy and a term dependent on the atomic positions [54, 55]. In order to measure the stress, we first introduce the kinetic energy in terms of a second-order tensor [56,57,58]:
Furthermore, the virial tensor \(\varvec{W}\) depends on the atomic positions and forces, and is defined via:
Neglecting all many-body interactions, one may simplify \(\varvec{W}\) to a virial that considers only pair potentials, \(\varvec{W}^{(2)}\). It can be written as:
This formulation allows for a demonstrative mechanical interpretation. The virial is the sum of all interatomic forces projected into the respective coordinate directions. Thus, even if the system remains in static equilibrium, where all forces \(\varvec{f}_i\) are equal to zero, the components of the virial tensor are not equal to zero, as the atomic interaction forces \(\varvec{f}_{ij}\) do generally not vanish. Finally, the pressure tensor is calculated taking into account both the kinetic energy tensor (56) and the virial tensor (57):
The Cauchy stress tensor, as also defined in continuum-mechanical formulations, is then equal to the negative pressure tensor: \(\varvec{\sigma } = -\varvec{P}\). In the absence of thermal vibrations, i.e., in the athermal limit, the Cauchy stress tensor is equal to the negative virial tensor multiplied by the factor \(\frac{2}{V}\).
The tensorial quantities discussed above are obviously stress measures that are averaged over all particles in the cell. However, one may also evaluate atomic level stresses with respect to every particle i surrounded by its local surrounding volume \(V_i\), cf. [53, 59,60,61].
5.3 Time Integration
Eq. (13) is a second-order nonlinear ordinary differential equation of motion. Due to the high number of atomic interactions, molecular systems require an efficient time-stepping method that approximates the dynamics over a large number of timesteps with a required level of accuracy. The molecular systems are highly nonlinear and, therefore, may exhibit sensitive dependence on perturbations, such as slightly different initial conditions. An overview of various integration schemes for molecular dynamics simulations is provided by Mesirov et al. [62]. Implicit time integration schemes build mostly on iteratively solving the linearized system. The high dimension of molecular systems leads to huge matrices, and the high number of interatomic interactions leads to a broad bandwidth of the matrices. Thus, implicit time integration schemes are rather unpopular in molecular modeling. Counting on explicit schemes, the most widely used method for integrating the equations of motions of molecular systems has been introduced by Verlet [63], which is still the predominant time integration scheme for such systems [30].
For the sake of compactness, we define the time derivatives of the quantities of interest as \(\dot{\varvec{r}}_i := \frac{d\varvec{r}_i}{dt}\), \(\ddot{\varvec{r}}_i := \frac{d^2\varvec{r}_i}{dt^2}\), et cetera. To obtain the Verlet integration scheme, one starts from the Taylor forward expansion of the position \(\varvec{r}_i\), from the time instant t to the time instant \(t + \Delta t\):
where \(\Delta t\) is the timestep length, and \(O(\Delta t^4)\) denotes all the remaining higher-order terms of at least fourth order. Equivalently, the Taylor backward expansion of the position \(\varvec{r}_i\), from the time instant t to the time instant \(t - \Delta t\), reads as follows:
The summation of Eqs. (60) and (61) leads to:
One observes that the odd order terms are no longer included in Eq. (62). Furthermore, the velocities \(\dot{\varvec{r}}_i\) are not needed to obtain the new position \(\varvec{r}_i(t + \Delta t)\). However, the velocities may be needed to determine crucial system quantities, such as temperature. For this purpose, the velocity \(\dot{\varvec{r}}_i(t)\) of each atom can be derived by the subtraction of Eqs. (60) and (61), which leads to:
It has to be noted that the position \(\varvec{r}_i(t - \Delta t)\) must be stored during forward stepping, but for the initial step, \(\varvec{r}_i(t - \Delta t)\) is unknown. Thus, the Verlet algorithm is not self-starting. Before performing the initial step, one must evaluate the artificial position \(\varvec{r}_i(-\Delta t)\) which allows one to initiate the algorithm together with the initial conditions. Moreover, the velocities are less accurate than the positions calculated via the basic Verlet algorithm. Whereas Eq. (62) is correct except for errors of order \(\Delta t^4\), the velocities calculated from Eq. (63) exhibit errors of order \(\Delta t^2\) [31].
The shortcomings of the basic Verlet algorithm can be resolved by slight modifications, such that the velocities appear explicitly. This so-called Velocity Verlet algorithm is derived by replacing t by \(t + \Delta t\), therefore shifting the formulation by the length of \(\Delta t\) in Eq. (61), such that:
The summation of Eqs. (60) and (64) leads to the velocity update:
Furthermore, the position update is obtained by inserting Eq. (65) into Eq. (64):
Eqs. (65) and (66) are the representations of the exact trajectories. One defines the numerical approximation of the exact dynamics by assuming dynamic equilibrium only at discrete time instants \(t_n\) and approximating the solution trajectory in between these time instants by neglecting the higher order terms \(O(\Delta t^3)\) from Eqs. (65) and (66). The numeric approximation is written as:
where \(\dot{\varvec{r}}_i^{(n+\frac{1}{2})}\) defines a half-step velocity term. The Velocity Verlet algorithm may then be implemented following a four-step procedure:
This procedure includes the position update in the second step. Based on the new positions, the accelerations are evaluated using the ordinary differential equation describing the dynamic system at step \(n+1\).
This velocity version of the Verlet algorithm is, to date, the most attractive and used algorithm due to its numerical stability, accuracy, simplicity, and time-reversibility, satisfying the conservation laws of energy and momentum [30, 31]. In the following subsection, we review the implementation of selected ensembles used for the melting-quenching method to obtain thermodynamic descriptions of disordered structures.
5.4 Microcanonical Ensemble
We start by transforming Newton’s second law of motion (13) from 3N ordinary differential equations into its canonical equivalent, consisting of \(2\times 3N\) equations:
where \(\varvec{p}_i\) is the momentum vector that the ith atom generates. The atoms experience an undamped motion in this micro-canonical ensemble [63], where the complexity of such motion mainly depends on the corresponding reaction force to which each atom is subjected. Applying the Velocity Verlet stepping algorithm (69) to the micro-canonical ensemble is straightforward, whereas the forces in step three can be directly evaluated from the second line in Eq. (70). The numerical realization of the microcanonical ensemble, using the Velocity Verlet algorithm, is written as:
Molecular dynamics simulations based on the microcanonical ensemble (70) occur on the condition that the number of atoms N, the system volume V, and the total energy E are conserved. However, this so-called micro-canonical, or NVE ensemble cannot reproduce real conditions, since it would be useful to perform simulations under a constant temperature T and constant pressure P. In the canonical or NVT ensemble, the system exchanges heat with its environment to achieve constant temperature. Furthermore, the so-called isothermal-isobaric or NPT ensemble can be used to conserve both pressure and temperature. Several techniques have been developed to control the temperature and pressure, cf. [31, 48, 64, 65]. In the following subsection, we apply and explain the technique also presented in Melchionna et al. [66].
5.5 Canonical Ensemble
One of the most frequently used algorithms for constant-temperature Molecular dynamics simulations was developed by Nosé and Hoover [67, 68]—the Nosé–Hoover thermostat. They extended the NVE ensemble to the NVT ensemble, where N, V, and T stand for a constant number of particles, volume, and temperature, respectively. This method makes use of an additional imaginary particle. In this way, the Lagrangian functional of the system, consisting of N atoms, is extended by an artificial degree of freedom, representing a heat bath or a thermostat [48], which is coupled to every atom and modifies its velocity. Essentially, this is realized by a friction parameter \(\zeta\), which may become positive or negative, to accelerate or decelerate the atoms and achieve dynamic states that oscillate around the desired temperature. The formulations in Eq. (70) are extended to the canonical ensemble, written as:
where the parameter \(\nu _T\) represents the thermostating rate, while T and \(T_d\) are the instantaneous system temperature and the desired target temperature, respectively. The equations of motion may again be integrated numerically using the Velocity Verlet algorithm, realized in four steps:
It must be noted that it can be beneficial to further extend single thermostats, as described above, to multi-degree-of-freedom thermostats, referred to as Nosé–Hoover chains. This may help to damp unphysical oscillations resulting in more stable simulations [69]. We note that other popular methods exist in the literature to control the temperature of molecular systems, such as the Andersen thermostat. For further information, the reader is referred to the relevant literature, cf. [48, 70].
5.6 Isobaric-isothermal Ensembles
The equations of motion of the canonical ensemble in (71) can be further extended to realize simulations under NPT conditions, whereby N, P, and T stand for a constant number of particles, constant pressure, and constant temperature, respectively. To realize such conditions, the system is allowed to change its volume. According to the Nosé–Hoover approach, this can be achieved by altering the atomic positions and velocities. Thus, in the isobaric-isothermal ensemble, the equations of motion take the form:
In this formulation, \(\varvec{r}_{\text {com}}=\frac{\sum _i m_i \varvec{r}_i}{\sum _j m_j}\) describes the center-of-mass of the system. The parameter \(\mu _P\) is a strain-rate factor that balances the deviations of the instantaneous system pressure P from the desired external pressure \(P_d\). The parameter \(\nu _p\) represents the barostating rate, \(k_B\) is the Boltzmann constant, d is the dimensionality of the system, and V is the system volume [66].
The equations of the isobaric-isothermal ensemble (72) describe volumetric changes of a system subjected to hydrostatic pressure, i.e., all bravais vectors of the simulation cell are adjusted by the same rate. This approach may yield satisfying results for fluids, since they reveal no—or very small—shear resistance. For solids, however, this may be insufficient. Thus, to simulate a box of varying shapes, the hydrostatic pressure P is replaced by the pressure tensor \(\textbf{P}\). Moreover, the scalar parameter \(\mu _P\) is replaced by a tensor \(\varvec{\mu }_P\). Then, the diagonal elements of \(\varvec{\mu }_P\) adjust the internal system pressure to the external hydrostatic pressure \(P_d\), and the non-diagonal elements of \(\varvec{\mu }_P\) eliminate the fluctuations of the non-diagonal elements of the pressure tensor \(\textbf{P}\). The equations describing the case of the hydrostatic pressure (72) are extended to:
where \(\textbf{I}_3\) is the second-order identity tensor. The trace of \(\varvec{\mu }_P = \varvec{H} \varvec{H}^{-1}\) has the property \(\mathrm {tr(\varvec{\mu }_P) = {\dot{V}}/V}\) [66]. The numeric integration procedure of the coupled set of equations of motion (72) and (73) can be efficiently realized using a Liouville-based approach that preserves the correct phase space volume element. For a detailed description of this approach, we refer to Tuckerman et al. [71].
6 Verlet and Linked Cell Lists
The most time-consuming processes in algorithms used in molecular simulations constitute the evaluation of the energy and the forces acting on each particle since every atom interacts with every other atom within the ensemble in terms of their shared potentials. Taking into account only pair potentials, \(N (N-1)\) interaction calculations are required. However, taking into consideration Newton’s third law \(\varvec{f}_{ij} = -\varvec{f}_{ji}\), as emphasized in Eq. (17), reduces the number of calculations to \(\nicefrac {N(N-1)}{2}\). Consequently, the computational effort scales quadratically with an increasing number of particles N, making simulations with only a few thousand atoms painfully slow.
6.1 Verlet Lists
Looking at the potential in Fig. 3, one can see that the interaction converges rather quickly to zero for larger interatomic distances. Therefore, one can conclude that the contributions of far-distant atomic pairs and triplets can be neglected. A cutoff radius \(r_c\) can be introduced that defines this circular interaction area, as indicated in Fig. 11. Noteworthy, taking into account the minimum image criterion, discussed in Sect. 3, the cutoff radius must not be larger than half the smallest box length.
Introduction of a cutoff radius \(r_c\) for a reduced evaluation of the atomic interactions and application of a cutoff-skin \(r_s\) to estimate the required re-evaluation of the Verlet list
Introducing a cutoff leads to a discontinuous jump in the functions of the pair potentials and their derivatives at the cutoff distance. While this approach is less dramatic when performing classical molecular dynamics simulations, where the system is integrated over time with an appropriate stepping scheme, as discussed in Sect. 5.3, it has a catastrophic influence on the convergence behavior of minimization strategies discussed later in Sect. 7. This issue can be overcome by ensuring that the potential energy function is zero and has zero slope at the cutoff distance by shifting and tilting the potential [72]. The shifted and tilted potential reads:
Every center atom is assigned to a neighbor list, called the Verlet list, which contains all neighbor atoms within the cutoff distance \(r_c\) [30, 63]. For the evaluation of the atomic force and the potential energy of every center atom, one only considers the atoms that are included in the neighbor list, while all remaining atoms outside the cutoff radius are neglected. The neighbor lists are updated after a certain number of calculation steps. The approach relies on the slowly changing atomic neighborhood, which implies that the neighbor lists remain valid throughout several calculation steps [30]. In this way, this approach avoids the time-consuming and unnecessary computational distance updates with atoms that are not contained in the neighbor list [31], in order that the speed of the program can be significantly improved. The number of required distance calculations using the Verlet list approach is evaluated by:
where \(\rho _n\) represents the atomic number density.
At the first step of the simulation, the initial construction of every list for every center atom requires evaluating the distances to all the remaining \((N-1)\) atoms, i.e., a full update. During the simulations, atoms outside the cutoff region may travel into the radial area, and atoms inside the radial area may have left the cutoff region. So, when should one perform such an expensive full update again to re-evaluate the Verlet lists? In order to stay on the safe side, one introduces a skin cutoff distance \(r_s\) resulting in a larger sphere of radius \(r_c + r_s\), and monitors the absolute value of the highest possible travel distance an atom can cover. The atoms inside the extended sphere of radius \(r_c + r_s\) are included in the Verlet neighbor list, while all atoms outside the sphere are neglected. The algorithm is quite straightforward. After one full update, one sets a maximum travel distance value \(d_{max}\) to zero. The atom that traveled the furthest within one calculation step needs to be identified, and its travel distance is added to the value \(d_{max}\). This value is then monitored after every calculation step. A new full update is required when an atom from outside could penetrate the cutoff skin in the radial direction. However, as the center atom also travels within every calculation step, the condition for the next distance update is: \(d_{max} \le \nicefrac {r_s}{2}\). The efficiency of this algorithm depends on the careful choice of the size of the skin. The larger the skin, the larger the number of atoms within the neighbor list, but the smaller the required number of full-distance updates. It should be noted here that updating the Verlet neighbor lists, even though performed every \(N_m\)th calculation step only, is still of order \(O(N^2)\).
6.2 Linked Cell Lists
Let us now minimize the time needed to replace a full update based on the simulation domain subdivision [73]. These methods partition the simulation box into a regular lattice of appropriately chosen \(N_{cell,x} \times N_{cell,y} \times N_{cell,z}\) smaller cells, to which the atoms are assigned according to their position relative to the corresponding cell [31, 73]. At first, we require the edge length of the cells to be greater than the potential cutoff radius \(r_c\). A two-dimensional schematic representation of this so-called linked cell list method is provided in Fig. 12a.
Linked cell-list method combined with the Verlet list method; a schematic representation of the combined Verlet and the linked cell list method with the cell edge lengths \(l > r_c+r_s\); b schematic representation of the combined Verlet and modified cell-linked method with the cell edge lengths \(l<r_c+r_s\). Applying the modified cell-linked method, the atoms contained in the light blue cells are excluded from the neighbor distance calculation
We note that this procedure is not limited to rectangular box shapes. In the case of parallelepiped boxes, the subdivision cells must also follow the equivalent parallelepiped shape. Furthermore, the method does not require the cells to have the same edge lengths, although it is advantageous if the edge lengths are relatively similar. From an implementation point of view, the approach builds on the property that monitoring if an atom, which is situated in a confined cell, is significantly less expensive than calculating the squared distances to all other atoms. As the atoms are assigned to the cells based on their current position, interactions between a center atom assigned to a particular cell are only possible with a neighboring atom if situated either in the same cell or in nearest-neighbor cells, i.e., the linked cells [30]. Thus, the force evaluation occurs via looping only over the neighbor atoms situated in the adjacent linked cells rather than over all remaining atoms [73]. Using the linked cell list method, the computational cost to calculate the interatomic distances is noticeably reduced to the order O(N). In three dimensions, the number of required interatomic distance calculations using the linked cell list method is evaluated by:
where \(\rho _n\) represents the atomic number density. This approach can be further optimized by decreasing the size of the cells, as presented in Fig. 12b. In this figure, we manipulated the linked cells from Fig. 12a so that every cell is again subdivided into \(2\times 2\) cells. In this case, the edge length is \(l<r_c+r_s\). Looking at the same configuration as in Fig. 12b, one realizes that cells emerge that do not have to be included in the update, as indicated in light blue in Fig. 12b so that the atoms in the light blue cells can be neglected. However, this modified, linked cell list method requires storing information regarding the neighboring cells and the second neighboring cells. This procedure requires one to increase the number of cell neighbors from \(3^3-1=26\) to \(5^3-1=124\).
6.3 Verlet and Linked Cell Lists
Using a combination of the Verlet and the linked cell list method significantly increases the efficiency of molecular simulations. In the case of a simple two-body potential, as presented in Eq. (29), it can be shown that using the Verlet list approach, the number of interatomic distance calculations is 16 times less than using the linked cell list. Thus, using the linked cell list method to build the Verlet lists was proposed in the literature. In doing so, the linked cell lists eliminate the disadvantage of the Verlet approach, which is the process of updating the neighbor lists and scales quadratically with the number of atoms, i.e., \(O(N^2)\) [48]. Figure 12 also shows a schematic representation of the combined Verlet and cell-linked methods. This figure considers only atoms in the blue cells to build and update the Verlet lists of the atoms in the center (red) cell. It is worth mentioning that, to better approximate the area of interaction, i.e., the potential cutoff sphere of radius \(r_c+r_s\), a slight modification of the conventional linked cell list method was proposed in the literature [73].
7 Locating Minima in the Potential Energy Landscape
Dependent on the respective temperature, solids have the property of remaining mostly in a local basin of a highly complex potential energy landscape (PEL), which depends on all atomic positions and the current cell shape \(\mathcal {U}(\varvec{r},\varvec{H})\), as discussed in Sects. 2.1 and 3. Figure 13a depicts a one-dimensional illustration of the PEL of a glass. Investigating glassy structures, the PEL can be divided into a two-level system consisting of global and local oscillations depending on the configuration coordinate. Globally, the PEL visits metabasins, while locally the system may be trapped in a local basin, as depicted by the red ball in Fig. 13a, b. Jumping from one metabasin into the other is unlikely for a solid at ambient temperature. One must heat the system, applying, e.g., the canonical ensemble, to create a melt and eventually cool it down again. Such a process is also referred to as \(\alpha\)-relaxation. However, under ambient conditions, the influence of the temperature in such a glassy solid can then be interpreted as a nonlinear random vibration of the configuration around the local minimum, as shown in Fig. 13b. With increasing time, thermal activation of the configuration so that it may fall into another basin of the PEL will—independently of the present magnitude of a positive temperature—finally occur, as depicted in Fig. 13c. Observing this illustration, it is easily imaginable that a jump from A to B is more likely than a jump from A to C. This indicates the already highly complex behavior of amorphous solids if no external mechanical deformation is applied at all. One can imagine that only statistical approaches make sense when it comes to the thermo-dynamical investigation of atomic structures and disordered solids, in particular [30].
One dimensional, schematic representation of the potential energy landscape (PEL) of a glass; a Global view of the PEL with metabasins, the \(\alpha\)-relaxation refers to a transition from one into another metabasin; b Local view of one basin (inherent structure), the structure oscillates around the local minimum (thermal vibration); c \(\beta\)-relaxation through thermal activation
When running the canonical ensemble at very low temperatures, the atoms will slightly oscillate around the inherent structure, i.e., the structure belonging to a local minimum of the PEL, at which the total forces of each atom vanish. Decreasing the temperature close to zero, one can imagine that the atoms will almost come to rest at this configuration. Finding the nearest local minimum of the PEL from a random configuration may be achieved by numerous algorithms of a constantly expanding minimization toolbox, which may be classified into two major groups, namely, line search type and trust-region types of algorithms [74]. Based on their superb effectiveness, simplicity, and applicability to molecular systems, we will concentrate on line search types of algorithms, particularly the two most popular methods: the steepest descent algorithm and the conjugate gradient algorithm. Since minimization algorithms play a crucial role in the athermal quasistatic deformation method, which will be discussed later in Sect. 9, we focus on introducing these two candidate methodologies in more detail, with particular emphasis on their application to molecular systems.
From an arbitrary atomic configuration at iteration step k, which is defined by the vector \(\varvec{r}_k\), one introduces a step length \(\alpha\) into a chosen searching direction \(\varvec{p}_k\), leading to the configuration at iteration step \(k+1\):
In this formulation, we choose the representation of the 3N-dimensional concatenated vector space, i.e., the three coordinates of all N atoms strung together, as defined in Eq. (2). An optimal search length for \(\alpha\) follows the one-dimensional minimization problem:
so that the value of the potential energy of the next configuration \(\varvec{r}_{k+1}\) is lower than the value of the potential energy of the current configuration \(\varvec{r}_k\). Finding an exact \(\alpha\) that fulfills the condition in Eq. (78) turns out to be costly. Instead, one proceeds step-wise downwards until the minimum is approximately found. Dependent on the applied algorithm, one may also adapt the searching direction \(\varvec{p}_k\) during the iterative search procedure. The variable \(k \in \mathbb {N}\) defines the running number of the iteration steps. Two quantities in Eq. (77) are a priori unknown: the searching direction \(\varvec{p}_k\) and the step length \(\alpha\). However, the search direction may vary during the iteration and is defined based on the specific algorithm used. The step length \(\alpha\) may a priori be chosen carefully by experience. If it is too small, there will not be a sufficiently substantial reduction in \(\mathcal {U}\). If too large, the algorithm diverges, leading to unphysical configurations. However, in general cases, it may be necessary to introduce certain ranges for \(\alpha\) using inequalities such as the Wolf conditions, as presented in the following subsections.
7.1 The Steepest Descent Algorithm
The most natural choice for \(\varvec{p}_k\) may be the direction towards the steepest descent of the potential energy, which is achieved by the negative gradient of the potential energy at the particular configuration \(-\nabla \mathcal {U}\left( \varvec{r}_k\right)\), which, in turn, is equal to the total force on every atom, here strung together in one 3N-dimensional vector \(\varvec{f}_k\). The configuration is updated iteratively and is evaluated as:
Thus, one iteration step must include the update of the atomic distances, the evaluation of the atomic forces, and the subsequent update of the atomic positions. The steepest descent algorithm guarantees convergence if a sufficiently small \(\alpha\) is chosen. However, although the numerical effort per step is relatively cheap, it may become excruciatingly slow to find the local minima of molecular systems, as the number of the required calculation steps is considerably high. Most convergence criteria are naturally defined by the absolute value of the difference in potential energies:
or the absolute norm of the force difference:
Figure 14 illustrates the convergence rate of the steepest descent method for a two-dimensional PEL.
Convergence behavior of the steepest descent algorithm
One observes a zigzag behavior, often leading to unacceptably slow convergence rates. We note that the zigzag behavior depends on the conditioning of the Hessian (introduced later in Sect. 9.1), i.e., the ratio of its largest and smallest eigenvalues. However, we emphasize that the steepest descent method can always be a backup option if faster minimization algorithms fail.
7.2 The Conjugate Gradient Algorithm
We start to introduce the group of conjugate gradient algorithms by firstly investigating conjugate directions for linear systems \(\textbf{A}\varvec{x}=\varvec{b}\), whereas \(\textbf{A}\) is a \(N\times N\) positive definite matrix [75]. One may introduce the equivalent minimization problem that results in the same solution \(\varvec{x}\):
If a step-wise procedure is applied, the gradient of this function at step k, \(\nabla \Phi (\varvec{x}_k)\), equals the residuum \(\varvec{r}_k = \textbf{A}\varvec{x} - \varvec{b}\). The basic idea is to define a set of N vectors that are conjugate with respect to the matrix \(\textbf{A}\):
This equation implies the necessary condition that the conjugate set of vectors \(\left[ \varvec{p}_1, \varvec{p}_2, \dots ,\varvec{p}_N \right]\) is also linearly independent. Subsequently, the minimum can be found by stepping consecutively into the conjugate directions via \(\varvec{x}_{k+1} = \varvec{x}_k + \alpha _k\varvec{p}_k\). For a linear system with a positive definite matrix \(\textbf{A}\) in Eq. (82), one needs a maximum of N calculation steps, whereby the step length is evaluated as:
However, this maximum number is not significant in practice since the algorithm converges, in general, much earlier. In theory, the best choice for such a set of conjugate directions is defined by the eigenvectors of \(\textbf{A}\), as they are mutually orthogonal and conjugate with respect to \(\textbf{A}\) [74]. An illustration of this minimization procedure for a two-dimensional quadratic PEL, realized by the set of the two eigenvectors, is depicted in Fig. 15.
Visualization of the set of two conjugate minimization directions \(\left[ \varvec{p}_1, \varvec{p}_2 \right]\) for a two-dimensional quadratic problem
Having defined conjugate directions, the conjugate gradient method is a special form of the conjugate direction method, where the conjugate vectors \(\varvec{p}_k\) are not known a priori. The new conjugate vector \(\varvec{p}_k\) is computed only by a linear combination of the negative gradient of the potential, i.e., the total \(3N\times 1\) force vector due to the current configuration, and the previous conjugate direction:
where the factor \(\beta _{k+1}\) is evaluated by:
Since molecular systems are highly nonlinear, we are not further presenting the algorithm for the linear system but directly addressing the general nonlinear problem. The extension to nonlinear systems concerns only the evaluation of the factor \(\beta _{k+1}\). We discuss the two most important extensions. The first possibility is the direct nonlinear equivalent of Eq. (86), presented by Fletcher and Reeves [74]. It is written as:
For nonlinear versions of the conjugate gradient method, it can happen that not every step is a descent step. Therefore, the step length \(\alpha _k\) must be chosen so that the strong Wolfe conditions [74] are satisfied:
where \(0<c_1<0.5\) and \(0<c_2<0.5\). The second possibility is an extension to Eq. (87), presented by Polak and Ribière [74]. It is written as:
Numerical studies show that the Polak-Ribière version tends to be numerically more stable than the Fletcher-Reeves version. However, the strong Wolfe conditions must be extended by a third condition:
We present a detailed step-by-step instruction of the Polak-Ribière version of the conjugate gradient method for molecular dynamics, which was implemented for the examples in this paper, in Algorithm 1.
8 Generation and Topological Identification
Two strategies for the generation of disordered materials are presented in this section. The first method is the melting-quenching approach, i.e., to cool down a melt until a solid emerges. This approach is a numerical reproduction of the natural glass fabrication process. It is, therefore, exclusively a strategy that allows for the generation of glassy structures, which is the first subgroup of disordered solids, as discussed in Sect. 1.2. The second method is the Monte Carlo bond switching algorithm, i.e., an athermal method that induces a sequence of local manipulations to an initially crystalline lattice. Using this approach, the result may be a glass or an amorphous solid, depending on the atomic composition and type of material.
8.1 The Melting-Quenching Procedure
Glasses are non-crystalline solids that exhibit a glass transition [17, 76, 77]. Thermodynamically, the configurations, i.e., the underlying inherent structures, of a melt continuously change, visiting wide sectors of the PEL at very low time scales, as presented by the red mark in Fig. 16a. They are generated by quenching the melt fast enough to prohibit a crystallization, i.e., to trap the system in a local minimum of the PEL, as depicted by the blue mark in Fig. 16a. During this intentionally fast quenching process, the glass visits the state of a supercooled liquid until it finally ends up in its solid glassy state, as depicted in Fig. 16b. The supercooled liquid state exists between the glass transition temperature \(T_g\) and the melting temperature \(T_m\). The structure of the glass solid is similar to that of its parent supercooled liquid state. The latest definitions of glass suggest that glasses are typical solids on a human time scale. However, their ultimate fate is to relax into their crystalline state, i.e., towards the absolute minimum of the PEL [17].
a Quenching a continuously changing configuration (glass melt) into an inherent structure within the PEL; b schematic enthalpy versus temperature plot showing four different states: liquid, supercooled liquid, glass, crystal
From the definition above, the natural motivation to numerically generate a glassy structure is to follow the process of cooling down the molecular system with the help of a heat bath using the NVT ensemble, as discussed earlier in Sect. 5.5 and elaborated in Eq. (72). One possibility for implementing such a numerical procedure is depicted in Fig. 17a. One starts from an initial configuration, either crystalline or random, at a temperature above \(T_m\), equilibrating the glass melt long enough until the initial state does not affect the structure of the system [19]. The equilibrating time depends mainly on the choice of the potential and the temperature.
a Melting quenching approach to generate a set of glass samples; b cutout of the equilibration and the subsequent quenching process of a binary Lennard Jones glass realized by the NVT ensemble and the Verlet integration scheme: integration time steps versus temperature T and damping parameter \(\zeta\); c the result of the melting quenching procedure: a binary Lennard-Jones glass sample; d bulk silica glass sample, visualized using Visual Molecular Dynamics [18]
Subsequently, the ensemble is further equilibrated for a long period, whereas melted configurations are extracted at equidistant periods. The equidistant periods must be chosen large enough to be considered independent from each other. This is a delicate process since the chosen temperature significantly influences the required equilibration time. For example, we found that an equilibration time of 25 picoseconds is sufficient to generate an independent silica glass configuration. Every extracted state is then quenched at appropriately chosen cooling rates, either to ambient or zero temperature, as shown in Fig. 17a. Such procedures are mainly realized using the canonical or the isobaric-isothermal ensemble integrated using the Verlet time-stepping scheme.
This paper presents two melting-quenching benchmark examples: a binary two-dimensional Lennard-Jones glass as a theoretical, qualitative model of densely packed glassy systems, as well as a three-dimensional binary vitreous silica structure as an example of a network glass.
Binary Lennard-Jones glass samples As our first benchmark model, we chose to generate a two-dimensional binary Lennard-Jones glass using the potential function presented in Eq. (29) via the melting-quenching procedure. For this example, we applied a unitless formulation as discussed in Sect. 2.2. The two-component system has a quasi-crystalline ground state and can easily be quenched into a stable glassy state. The system consists of atoms of two different sizes, the small atom type (S) and the large atom type (L), resulting in the zero-energy distances, written as:
and in the bond strengths, written as:
The number ratio of the particles is chosen as:
which is equal to half of the numerical value of the golden ratio, which is approximately 1.618. The system is energetically favorable for ten small atoms to surround one large atom or, vice-versa, five large atoms to surround one small atom in order to avoid local crystallization. This binary model system undergoes something like a glass transition when quenched [78]. In our simulation, we defined the mass of a large atom type to be twice as large as the mass of a small atom type. Time stepping was realized using the Velocity Verlet integration scheme. A system of 900 randomly placed atoms, i.e. 498 small and 402 large atoms, was equilibrated for \(10^6\) time steps at a temperature of \(1.5\,T_m\) using the canonical ensemble. Subsequently, the binary system was quenched for \(2\times 10^3\) steps, during which the ensemble was forced to cool down to the (almost) zero temperature state \(T_0=1.0\times 10^{-5}\,T_m\). The period of the last \(10^5\) integration steps and the subsequent \(2\times 10^3\) quenching steps are shown in Fig. 17b. The system oscillates around the target temperature function, whereas the size of the oscillations depends on the choice of the parameter \(\nu _T\). The magenta plot presents the corresponding time evolution of the damping parameter \(\zeta\) in Fig. 17b. During the quenching period, the damping parameter considerably increases, slowing down the Brownian motion of the particles to almost zero. The quenched result of the binary Lennard-Jones glass is depicted in Fig. 17c. It was used for the mechanical and thermodynamical investigations of several studies as the first theoretical benchmark model for a binary metallic glass [79,80,81,82].
Bulk silica samples To obtain our second benchmark model from the family of network glasses, we applied the melting-quenching procedure for a set of bulk silica glass samples. Here, the more sophisticated Stillinger-Weber type potential (157) was used, which consists of two- and three-body contributions. We note that the purely two-body Best–Kramer–Santen type of potential, presented in Eq. (156), results in similar atomic configurations and is a comparably popular choice in the literature [83]. Units are measured according to the definitions in Appendix A. One Si-atom and two O-atoms are positioned in an elementary unit cuboid box. The length of this box is chosen as 3.565 Å so that the material has the same density as observed from experiments, i.e., 2.2 g/cm\(^3\) [19]. When positioning the atoms within the elementary box, one must carefully choose the coordinates so that the interatomic distances do not get too small. As an example, the chosen coordinates of \(\left[ 1.0, 1.5, 1.75 \right]\) Å, \(\left[ 0.5, 3.0, 0.5 \right]\) Å, and \(\left[ 3.0, 3.0, 3.0 \right]\) Å regarding the Si-atom and the two O-atoms, respectively, constitute an appropriate initial configuration and lead to a stable time integration process. This elementary box is replicated in the three coordinate directions. The randomization of the system originates from the Boltzmann distribution that is used to define the initial conditions for the atomic velocities. Using a time step of 0.5 fs, a system of ten replications in every direction, which corresponds to an ensemble of \(3\times 10^3\) atoms, was equilibrated for one ns at a temperature of \(8\times 10^3\) K, using the NVT ensemble, as presented in Eq. (72). Subsequently, the melted configuration was quenched to a very low temperature of 0.01 K using a quenching rate of 50 K/ps. We note that such a quenching rate is many orders of magnitude larger than those performable in a laboratory. In reality, the quenching process is realized by heat exchange through the surfaces instead of a damping parameter affecting every atom of the ensemble simultaneously. Defining a physically meaningful quenching rate is still a current topic of scholarly discussion [84]. In this regard, many studies accept the numerical incapability to describe macroscopic time scales using molecular dynamics models and focus on generating a physically meaningful configuration and comparing it to experimental results. Intriguingly, a machine-learning approach has lately been developed to generate a silica glass structure [85]. Comparing the results with Stillinger-Weber-type potentials, this alternative approach is highly efficient while providing the essential features of a network glass structure.
Bidault et al. [86] have shown that, during the quenching process, the silica melt reveals a structural anisotropy, and the melt shows non-Newtonian and strain-rate dependent properties when subjected to external loading.
The quenched result of one sample is depicted in Fig. 17d. One obtains a complex disordered network structure fundamentally different from the binary Lennard-Jones glass in Fig. 17c. In this figure, the red balls represent the Si-atoms, while the green balls represent the O-atoms. Vitreous silica is assembled by SiO\(_4\) tetrahedra that build complex networks of various shapes and sizes. One may define three levels to classify a silica glass network structure.
The first level concerns the intra-tetrahedral structure, and gives information on the architecture and the shape of the SiO\(_4\) tetrahedra. The second level concerns the inter-tetrahedral structure and gives information on the arrangement of the SiO\(_4\) tetrahedra with their neighboring tetrahedra. Both levels can be quantified statistically using the concept of the radial distribution function (RDF) and the bond angle distributions (BADs). These two statistical quantities can also be extracted experimentally and are mainly used to verify the numerical studies of vitreous silica [87, 88].
The RDF \(g_{\alpha \beta }\) represents the probability density of atoms of two different types, \(\alpha\), and \(\beta\), at a certain distance r. The number of particles in a spherical shell of a thickness \(\Delta r\) is equal to the radial distribution multiplied by the shell volume and the particle density:
The quantity \(c_\beta\) defines the concentration of atom type \(\beta\), and \(\rho\) defines the total number density.
The BAD represents the second identification level and constitutes the probability of occurrence of a particular angle spanned by the two vectors \(\varvec{r}_{\beta \alpha }\) and \(\varvec{r}_{\beta \gamma }\) within every possible triplet that includes the atoms of type \(\alpha\), \(\beta\), and \(\gamma\). Practically, the angles are only evaluated between nearest-neighbor atoms of the respective type [19].
The third level constitutes the ring structure of the corner-sharing SiO\(_4\) tetrahedra, which has become an accepted measure for the medium-range order of network materials, such as vitreous silica. One may count the number of SiO\(_4\) tetrahedra, i.e., the number of Si-atoms, of every ring in the ensemble and evaluate a statistic of rings embedded in the network structure. A ring that is optimal in terms of potential energy consists of six corner-sharing tetrahedra, i.e., six Si-atoms. Using the melting-quenching procedure, applying such a low cooling rate so that the configuration will finally fall into the absolute minimum of the PEL, i.e., a crystalline state, is numerically impossible. Thus, the system will always be trapped in a local minimum of the PEL, so that the network also consists of rings that reveal a smaller and a larger number of Si-atoms than six. The broader the histogram of the ring statistics, the larger the level of heterogeneity, i.e., the larger the number of smaller and larger rings [89, 90].
a Radial distribution functions; b bond angle distributions; c visualization of a six and a twelve-membered ring, reprinted from Ebrahem et al. [89] with permission from Elsevier; d ring statistics
Vollmayr et al. [91] showed that the choice of the quenching rate does influence the structure of the silica network. Their study includes the short-range order of the network structure, i.e., the geometry of the tetrahedral SO\(_4\) units (represented by the RDF) and the orientation angle between these units (represented by the BAD). Later, it was shown that the medium-range order, quantified by the statistic of rings (a six- and a twelve-membered ring are exemplarily visualized in Fig. 18c), can be influenced more highly when changing the quenching rate. The higher the cooling velocity, the larger the level of heterogeneity, i.e., the broader the histogram in Fig. 18d [89].
While the statistical identification of the RDFs and the BADs is relatively straightforward, evaluating the statistics of rings of a network structure, such as silica glass, becomes rather challenging and computationally expensive. King [92] was one of the first to demonstrate the complexity of the possible ring structures in a network glass, such as vitreous silica. He pointed out that the rings in the crystalline polymorph, \(\alpha\)-cristobalite, consist only of six-membered corner-sharing SiO\(_4\) tetrahedral ring structures. The possibilities of rings in a vitreous structure are almost infinite when considering them as closed paths in the network of corner-sharing SiO\(_4\) tetrahedra. Thus, he pointed out that one could consider the definition of the shortest path to select the physically relevant rings only. Thorpe [93] emphasized the relevance of the ring formations qualitatively, discussing the local rigidity levels of vitreous networks. He pointed out that rings with fewer than six members are rather rigid, while rings with seven or more members are rather floppy. These thoughts lead to a complex spatial distribution of local flexibilities that does not only depend on one local ring but also its neighborhood rings. Guttmann [94] provided a clarification of a shortest-path criterion, which was first mentioned by King [92], and stated that the ring search procedure is a straightforward, distinct problem as it eliminates all closed paths that exhibit “shortcuts” across them. Below, we briefly review the essential minimum data from the field of graph theory [95] necessary to extract a set of physically relevant rings from an atomic network by realizing the shortest path criterion, which was first introduced by Franzblau [96].
By definition, a graph (or a network) \(G=(V,E)\) consists of a set of vertices V and a set of edges E, naturally describing atoms and bonds between atoms in our context. It is noteworthy that the structural information of the network topology of silica glass can be identified by placing the vertices at the coordinates of the Si atoms. The edges describe the set of all possible Si pairs in the sample. If an O-bridge links a respective Si atom pair, they constitute an adjacent tetrahedral pair that shares a corner. In other words, the network structure of silica glass is equivalently defined by considering the linked Si-atoms only. Defining two vertices, y, and z, in a network, a path is a chain of subsequent edges linked together. To ensure the existence of such a chain, every vertex in the path must be a part of at least two edges. Thus, it is possible to define a chain of edges, i.e., paths, by listing the subsequent vertices along the path. The number of edges in a path defines its length. Naturally, one can define a ring as a path that returns to its starting point, i.e., the start and end vertices are identical. Proceeding along this closed path, the size of a ring may then be defined by the number of edges, k, contained in this closed path, which is identical to the number of vertices. Proceeding one edge forward and the same edge back again is impossible in vitreous silica or at least very unlikely. In other words, the structure of a two-membered ring is impossible, as the smallest possible closed path found in vitreous silica contains at least three vertices. On the one hand, the distance between two vertices, y, and z, is the minimum number of subsequent edges, k, out of all possible paths from y to z. On the other hand, the diameter between two vertices, y, and z, is the maximum possible distance out of all possible paths from y to z. The diameter of a ring is \(\frac{k}{2}\) if the ring contains an even number of vertices, and it is \(\frac{k-1}{2}\) if the ring contains an odd number of vertices.
A disproportionally high number of closed paths can be found in a graph/network. However, not all of those closed paths are relevant from a mechanical point of view. Having defined the required tools from the graph theory above, we can extract the set of “primitive rings”. Franzblau [96] and Hobbs et al. [97] define a primitive ring simply as a closed path that can not be divided into any sub-paths. Thus, we present the following definition of a primitive ring, according to Franzblau [96] (Fig. 19).
Definition Given a graph G and a ring R in G, R is a shortest path ring if R contains a shortest path for every possible vertex pair on the ring. We emphasize this statement with the help of the simple network example, depicted in Fig. 20.
Minimalist example for extracting possible paths for the ring search algorithm and sorting out primitive rings
Out of six possible rings (closed paths) one can count the following four primitive rings in this network: abcfea, abdfea, bdfcb, and bgdb. One already observes that primitive rings are not entirely exclusive as they may share sub-paths. The ring abgdfea does not count as a primitive ring as one can find a shortcut bd that divides this closed path into smaller closed paths abdfea and bgdb. Equivalently, the ring bgdfcb does not count as a primitive ring as one can find the shortcut bd that divides this closed path into smaller closed paths bdfcb and bgdb.
Implementation-wise, the ring search algorithm consists of two procedures. The first procedure is to find the collection of all closed paths in the network. In the case of vitreous silica, one starts to form a randomly chosen Si-atom, from which six exit-incoming possibilities for the closed-path search are considered, as presented in the top right-hand illustration in Fig. 20 in blue, green, magenta, orange, red, and yellow color.
Evaluation of the medium range order: finding closed paths in the network structure considering the all income-exit possibilities of a vitreous silica network
Following one of these six income-exit possibilities, one proceeds from Si-atom to Si-atom over the corresponding O-bridges. The path possibilities multiply by a factor of three with every Si-atom, as indicated by the blue path in Fig. 20. Realizing a search algorithm up to n-membered rings, one obtains a maximum of \(3^{n-1}\) possible closed path candidates. This approach is performed for all six income-exit possibilities. All path candidates that do not enter over the incoming oxygen bridge are erased from the candidate list so that an initial set of closed paths is found for every Si-atom (vertex) in the system. One essential restriction that must be considered when performing the closed path search for periodic systems is ensuring that the candidates’ circumference of k vertices is smaller than the box length. If this criterion is not fulfilled, one detects closed paths over the periodic boundary conditions that are not rings. The second procedure is to extract all primitive rings by excluding all rings that contain shortcuts when applying the primitive ring criterion defined above. However, it is not necessary to check if all possible combinations of pair distances are sub-paths of the ring candidate. For the sake of efficiency, it is sufficient to check if every diameter of the ring is a shortest path. Therefore, one has to start from one vertex in two different directions, while the shortest path to the opposite vertices must be a sub-path of the ring candidate. It is noteworthy that, in the case of an even ring, one has one opposite vertex, while, in the case of an odd ring, one has two opposite vertices, which is in line with the above-defined diameter of a ring. For a k-membered even ring candidate, this must be repeated \(\frac{k}{2}\) times, and for a k-membered odd ring candidate, this must be repeated \(\frac{k-1}{2}\) times [89]. If the shortest paths found are all sub-paths of the ring candidate, i.e., if no shortcuts can be found, the candidate ring is indeed a primitive ring and can be added to the already evaluated primitive rings list. As the total ring search algorithm is performed for every vertex, one finds every k-membered primitive ring k times. Thus, all multiple primitive rings must be erased from the entire list as a final step.
Evaluating the medium-range order in terms of the ring statistics is significantly more expensive than identifying the short-range order realized by the RDF and the BAD. Especially if large rings occur in the network structure, the ring search algorithm becomes disproportionally expensive. Thus, the realization of the ring search algorithm has been further optimized in the literature [89, 98,99,100].
8.2 Monte Carlo Bond Switching Algorithms and the 2D Zachariasen Model
Analogously to the motivation of creating a two-dimensional binary Lennard-Jones model glass for theoretical investigations, as shown in Fig. 17c, one may also choose to create a two-dimensional version of vitreous silica, as shown in Fig. 17d. Firstly, one obvious advantage of investigating a purely two-dimensional model is the possibility of visualizing—one may even call it imaging—which allows for better classification and understanding of complex mechanical phenomena. Secondly, abandoning one dimension leads to significantly computationally cheaper models than three-dimensional models of similar sizes. Thirdly, using such purely theoretical model materials, one can bother less about the quantitative correctness or the level of accuracy of the model compared to its experimental counterpart, which allows one to mainly focus on the complex phenomena of the highly nonlinear chaotic physical and mechanical behavior of disordered solids. Adapting the two-dimensional observed phenomena to their more realistic three-dimensional parent models is often a technical task.
By generating a two-dimensional analogy of vitreous silica, we first want to introduce the random network model, introduced almost ninety years ago in the pioneering work of William Holder Zachariasen [101]. With the help of a two-dimensional cartoon model of a crystalline and a vitreous state of silica, as shown in the left and the right subplot of Fig. 21a, he correctly argued that vitreous silica does not have a crystalline arrangement of rings; instead, it has a disordered arrangement of rings of various shapes and sizes, which he described by the random network theory. The latter consists of three simple but essential rules, all of which are consistent with our numerical vitreous bulk silica sample in Fig. 17d: (i) one O-atom is not linked to more than two Si-atoms; (ii) the number of oxygen atoms surrounding an Si-atom is, in general, smaller or equal to four (in the case of an ideal two-dimensional network, it is equal to three); (iii) the tetrahedra, or, in two dimensions, triangles, share corners but no edges, i.e., rings with less than three Si members are not possible or at least very improbable.
Two-dimensional models of vitreous silica; a Zachariasen’s two-dimensional analogy of the crystalline and the vitreous state of silica glass, the figures are taken from the paper by Zachariasen [101] and reused with permission from ACS; the numerically relaxed network structure of an experimentally imaged 2D silica material from the crystal to the vitreous state [104] (b) and the vitreous state [104] (c); note the intriguing similarity with Zachariasen’s cartoon models [101] (a)
While the experimental extraction of RDFs and BADs has been possible for half a century [87], the rough microscopic surface of silica glass has definitively proved the network theory in terms of microscopic imaging impossible. However, eighty years after Zachariasen’s theory, the generation of a flat 2D silica material enabled such an imaging procedure [102,103,104,105] using transmission electron microscopy. Figures 21b and 21c present images of such a real two-dimensional silica glass material. To present the corresponding atomic models in those two figures, we took the atomic coordinates from Lichtenstein et al. [103]. We relaxed the ensemble to zero temperature into the nearest local basin by minimizing the potential energy using the Yukawa-type potential, as discussed in Sect. A, and the conjugate gradient algorithm, as discussed in Sect. 7.2. These figures are numerical reproductions from real images, so periodic boundary conditions are not realizable. Therefore, the position vectors of the atoms belonging to the thin border layers are frozen and excluded from the position updates during the minimization procedure. Figure 21b shows both the crystalline and the vitreous state and the transition between the two states, while Fig. 21c shows a vitreous network structure. Although this material is different from bulk silica, as it is a mirrored bilayer linked over the Si-atoms by oxygen bonds [106], it gives an idea about the natural statistic of rings. Therefore, the architecture of the medium-range order of network materials, such as silicate glasses [107,108,109]. Not only the mechanical and physical properties of two-dimensional layered silica have been the focus of interest [110,111,112,113], but the extracted picture of the silica network also serves as a benchmark of two-dimensional networks in nature [114,115,116,117]. Thus, it serves as a role model for physically meaningful two-dimensional network glasses, used to investigate and further understand the thermodynamics and deformation mechanics of network glasses, in general, as will be discussed in detail in Sect. 9. We note that also other experimentally imaged two-dimensional network structures exist in the literature. In this context, the authors acknowledge the papers by Huang et al. [118, 119], whose images could also be used as a benchmark model for the numerical generation of networks.
The numerical generation of such materials has lately been performed by molecular dynamics and purely athermal procedures. The latter was performed by a random series of local topological manipulations. We will review both strategies below.
The thermodynamic process was performed by an adaptation of the melting-quenching process, as shown in Fig. 17a. Roy et al. [115] created two-dimensional silica melts of samples between 80 and 500 atoms and integrated the canonical ensemble in time using the Nosé–Hoover thermostat for a large number of time steps. At equidistant time intervals, they minimized the structure using the conjugate gradient algorithm to find the closest basin, i.e., the melt’s underlying inherent structure. For this procedure, they adapted the Yukawa-potential function to fit the model material SiO\(_{1.5}\), i.e., the microscopic picture of the network structures imaged by Lichtenstein et al. [103], as shown in Fig. 21b, c. The network topology can be tuned by varying the temperature from which the melt is minimized. The higher the temperature, the higher the level of heterogeneity of the network structure. This can be adapted until the statistics of rings of the numerically generated sample is comparable to the one of the real vitreous material, presented in Fig. 21c. As elegant as this procedure is, it has one disadvantage, as also claimed by the authors. For obvious statistical reasons, the number of samples that fulfill the requirements of Zachariasen’s cartoon model and the topological constraints in terms of coordination number extracted from the images, i.e., every Si-atom has three neighbors, and every O-atom has two neighbors, are very small compared to a large set of minimized inherent structures. With an increasing number of atoms, the possibility of finding such a “perfect” sample that topologically matches the role models decreases significantly, leading to a numerical threshold regarding the size of the sample to date. One concludes that the generation of a whole set of samples of appropriate sizes used for mechanical investigations turns out to be an unfeasible task using the thermodynamic quenching procedure. A numerically more efficient procedure must be found that allows for the generation of Zachariasen network structures.
The alternative strategy to obtain such a network structure is realized via the Monte Carlo bond switching algorithm, which is generally a randomized series of topological flip transformations in the network structure [120]. One starts from a purely hexagonal lattice and performs subsequent flip transformations of a bonded triangle SiO\(_3\) pair. The first flip transformation is depicted in Fig. 22, where one pair of SiO\(_3\) triangles, linked by an O-bridge, is selected and highlighted in red. This triangle is rotated, followed by minimization of the potential energy so that the configuration can relax towards its nearest energy basin. The switched result is depicted in Fig. 22b.
Topological flip transformation of a 6-6-6-6 configuration within a pristine honeycomb lattice structure to a 5-7-5-7 configuration (Stone-Wales defect)
It can easily be seen that every Si-atom still has three covalent O-neighbors, and each O-atom still has two Si-neighbors. Therefore, the atomic structure is still fully-coordinated. Initially, the triplet pair is surrounded by four six-membered rings, i.e., a 6-6-6-6 configuration. After the switch, the transformed bond is surrounded by two five-membered and two seven-membered rings, a 5-7-5-7 configuration. This particular case is often referred to as Stone-Wales (SW) defect [121]. It is evident that the flipping rotation angle is exactly \(\nicefrac {\pi }{2}\) in the particular case of a purely hexagonal lattice that is subjected to one flip transformation. To obtain a Zachariasen random network, one must perform several subsequent random switches in the network structure—a Monte Carlo approach. One switch is then considered a Monte Carlo move. Firstly, one randomly chooses an adjacent triangle pair and geometrically rotates it around its center. After this purely geometrical transformation, the molecular system is minimized using the Yukawa-type potential from Eq. (162). Since, at this point, the choice of the triangle pair to be flipped is purely random, it is possible that very unlikely or unphysical configurations might emerge. Thus, it makes sense not to accept every randomly proposed switch by introducing the Metropolis-Hastings acceptance condition [122, 123], here slightly adapted to the problem at hand according to Kumar et al. [124]:
where \(\mathcal {U}_b\) and \(\mathcal {U}_a\) denote the potential energy before and after the switch, \(k_B\) is the Bolzmann constant and T is the temperature. In general, the lower the T is, the higher the overall probability of acceptance. Therefore, the bond switching procedure can be seen as a Monte Carlo Markov Chain method, where sampling is realized by the Metropolis-Hastings algorithm. In summary, the following conditions must be fulfilled to accept a bond switch finally:
-
(i)
acceptance through the Metropolis-Hastings condition in Eq. (96);
-
(ii)
the configuration remains fully coordinated; and
-
(iii)
the ring sizes remain in an acceptable range in accordance with experimental images (in the case of 2D silica, four to ten-membered rings are allowed).
We emphasize that, throughout the procedure, only one interaction potential is generally needed. However, in order to save computational effort, one may perform a local relaxation near the switching center, check the consistency, and then relax the whole sample. In doing so, it may be advantageous to use several interaction potentials consecutively during one Monte Carlo move. In particular, Kumar et al. [124] used the Keating potential to relax the configuration locally around the switching center, i.e., only the four rings involved. In this way, the switching is efficiently repeated until one obtains an accepted state according to conditions (i) and (ii). Subsequently, they relaxed the whole configuration and, finally, accepted or rejected the switch based on condition (iii). Thereby, it is also possible to get rid of the oxygen bonds and only model the network using Lennard-Jones, Keating, or Tersoff potential functions [124, 125]. Following this procedure, archetypal two-dimensional glass models that follow Zachariasen’s framework have been generated using the Monte Carlo bond switching algorithm and experimental inputs as benchmarks [124], as presented in Fig. 21b, c.
As discussed above, the rotation procedure during one switch is a geometrical manipulation, also proposed by numerous studies before. Since this manipulation is purely geometrical before one relaxes the molecular structure, one may argue that, strictly speaking, one Monte Carlo move may not be thermodynamically consistent. Experimentally, one bond switch can be induced by locally exciting the material with electron beams, as presented by Huang et al. [119]. To apply a switching procedure that is thermodynamically consistent, one may try to simulate the experimental process by locally increasing the kinetic energy in an NVT ensemble until a flip transformation occurs and then relaxing the structure. This should be repeated for the whole sequence of switches until a molecular structure is obtained that reveals the desired level of disorder. This strategy is numerically unfeasible and has not been pursued further in the literature. The accepted approach is, therefore, to mainly concentrate on altering the network topology and follow consistency checks after every Monte Carlo move.
On this occasion, we take the opportunity to note that such configurations are ubiquitous in several types of 2D materials, whose pristine configurations reveal the structure of a honeycomb lattice, opening new possibilities for the application of 2D networks. The most important example here constitutes the family of graphene monolayer structures, firstly generated by Geim and Novoselov [126,127,128], which exhibit extraordinary mechanical properties [129, 130]. Graphene monolayers have C-atoms instead of Si-atoms and, obviously, no bridging atoms, and the lattice size is smaller by a factor of about two [125]. The topological defects significantly alter the mechanical and electronic properties of the 2D material [131,132,133,134]. Lately, it has even been possible to generate graphene glass structures that fulfill Zachariasen’s continuous network theory experimentally [135, 136] and numerically [137]. However, regarding the rich body of the structure-property relationship of this and related 2D materials, we refer to the extensive literature that deals with these particular type of materials and their promising future in electro-mechanical applications [138,139,140,141,142,143,144,145,146]. The motivation here is rather to extract and generate realistic vitreous and/or amorphous network configurations for two-dimensional thermo-dynamical and purely athermal mechanical investigations.
The dual Monte Carlo bond switching algorithm follows the same idea as the classical Monte Carlo bond switching algorithm. The main difference is that the switching procedure is performed on a dual lattice, which carries the information of the ring structure. This is particularly convenient since a flip transformation is seen as a topological defect that only alters the ring structure but not the coordination number. Below, we present the dual Monte Carlo bond switching algorithm that allows one to control both the level of disorder in terms of the ring statistics and the ring-neighborhood statistics. The algorithm starts from the crystalline lattice, as depicted by the blue structure in Fig. 23a. Every honeycomb ring is assigned a dual node placed in its center, as depicted by the black nodes in Fig. 23a. Thus, the initially assigned integer value to every dual node is six. Every dual node is linked by a dual bond with a corresponding dual neighboring node if the overlying rings of the physical lattice share a bond, i.e., if the two rings are adjacent. This situation leads to a reduced representation of the network structure, represented by a dual network of tessellating equilateral triangles, where every dual node has six neighbors in case of a pristine, hexagonal configuration.
The physical lattice is described by the Yukawa potential for two-dimensional silica model structures, as discussed in Appendix A. Conversely, the dual lattice can be described by a harmonic dual pair potential, or referred to as ring-ring potential by Morley and Wilson [147], which depends on the types of the dual node pair and which is equivalent to the sizes of the overlying two adjacent rings. The harmonic dual potential is written as:
In this potential function \(k_h\), denotes a factor for the strength of the interaction, which we set to 1.0 for the demonstrations in this paper. The value \(r_{ij}^{(0)}\) denotes the equilibrium distance of two dual nodes, i.e., rings, of sizes \(n_j\) and \(n_j\). It is evaluated by:
where \(r_{bl}\) refers to the bond length of the physical structure, i.e., the equilibrium distance of every triangle center.
Instead of switching the physical bond consisting of two SiO\(_3\) triangles by \(\nicefrac {\pi }{2}\), one switches the dual bond, as presented by the solid red line in Fig. 23a. This is done by redefining the dual node connections, as presented by the red dashed line in Fig. 23a. The integer type of the dual nodes belonging to the bond before the switch is reduced by one, while the integer type of the dual nodes belonging to the dual bond after the switch is increased by one. Subsequently, the system is relaxed using the harmonic dual potential. The corresponding overlying Si atoms of the physical network structure are obtained by placing them in the dual triangle centers, and the O-bridges between triangle pairs are positioned to link neighboring triangles. Additionally, the physical structure is relaxed using the Yukawa potential to obtain the final network, as discussed in Appendix A.
For illustration, we track a sequence of the first three flips of such a crystalline-to-vitreous transformation. Starting the procedure, a 6-6-6-6 configuration in Fig. 23a is topologically transformed into a 5-7-5-7 configuration in Fig. 23b. This includes the minimization of the dual potential in Eq. (97) on the dual network after the switch and the subsequent minimization of the physical lattice using the Yukawa potential (162). This first flip is the realization of the SW defect [121] in a two-dimensional silica structure.
Three subsequent topological flip transformations using the dual version of the bond switching algorithm; a hexagonal structure and the first selected dual bond of a 6-6-6-6 configuration; b switched dual bond to 5-7-5-7 configuration followed by relaxation of the dual network and the physical network; c and d switch of a 5-6-6-7 to a 4-7-5-8 configuration; e and f switch of a 6-6-7-6 to a 5-7-6-7 configuration
Starting from the configuration in Fig. 23b, one defines the next dual bond, as depicted by the red line in Fig. 23c. One then topologically transforms the respective dual bond structure by redefining the dual node connections of the corresponding dual bond, as shown by the dashed red line in Fig. 23c. A local 5-6-6-7 configuration was topologically transformed into a 4-7-5-8 configuration, followed by the minimization approach on the dual and the physical lattice, as discussed above. The result of this flip is shown in Fig. 23d. Finally, we present a third topological transformation in Fig. 23e, f by switching a local 6-6-7-6 into a 5-7-6-7 configuration, followed by the minimization, equivalently to the two flip procedures before. Choosing a sequence of dual bonds randomly can finally lead to a physically meaningful network structure, as also microscopically imaged by Lichtenstein et al. [104] (see Fig. 21c) since the relaxation of the overlying network structure using the Yukawa potential (162) provides the final physical check of the configuration. If the network is not fully coordinated, i.e., if every Si-atom is not surrounded by three O-atoms or every O-atom is not surrounded by two Si-atoms, the system is not consistent with Zachariasen’s network theory or experimental observations. This means that the result of the latest flip transformation must be abandoned, and one goes back to the configuration before the flip. In order to save significant computational time, one may perform the dual bond switching algorithm purely on the dual lattice. However, performing a random sequence of switches will most probably not lead to a physically meaningful network structure. For instance, we subjected a \(15\times 18\) cell hexagonal structure to a sequence of 200 random switches, where only the dual structure is minimized using the harmonic ring-ring potential. The result is presented in Fig. 24a.
a Network after 200 random switches; b Network after 200 switches controlling the ring neighborhood statistics
As shown in this figure, a purely random sequence of switches may lead to unphysical clusters of small and large rings. One possibility to avoid that can be to use an empirical investigation, introduced by Aboav and Weaire [148,149,150]. The Aboav Weaire law can be written as follows:
This empirically observed relation evaluates the mean ring size \(m_n\) that is naturally found around all rings of size n in a particular network structure. The magnitudes \(\mu _n\) and \(\sigma _n\) denote the mean and the variance of the ring statistics. Formula (99) can then be used to control the appearance of the mean ring size around a center ring, which one may interpret as ring neighborhood statistics. Physically meaningful ring configurations, according to Aboav and Weaire [149], reveal values of around \(\alpha _n = 0.3\). This means that, in a natural silica network structure, it is more probable to find small rings surrounding large rings and vice-versa, while clusters of small or large rings are rather improbable. We extracted the ring neighborhood statistics of the measured configuration in Fig. 21c and evaluated Aboav-Weaire parameters of \(\alpha _4 = 0.31\), \(\alpha _5 = 0.36\), \(\alpha _7 = 0.29\), \(\alpha _8 = 0.31\), and \(\alpha _9 = 0.32\) for center ring sizes of four, five, seven, eight and nine corner-sharing triangles (Si-atoms), respectively. In order to create physically meaningful networks in terms of the ring neighborhood statistics, one might choose target values for \(\alpha _4^{(t)}\) and \(\alpha _5^{(t)}\) and evaluate the difference \(\vartheta _\alpha\) between the target and the actual values of the configuration via:
where \(c_1\) and \(c_2\) are constant weighting factors. One may interpret this difference as an objective function that must be minimized during the switching procedure in order to obtain the topological target statistics. The constant weighting factors \(c_1\) and \(c_2\) define the importance of the objective function \(\theta _\alpha\). For the examples presented in this paper, values between 0.1 and 10.0 are feasible. The numerical results in this paper are produced by choosing a value of 2.0 for the weighting factors \(c_1\) and \(c_2\). If the change in \(\vartheta _\alpha\) during the kth switch \(\Delta \vartheta ^{(k)} = \vartheta _{\alpha }^{(k+1)} - \vartheta _{\alpha }^{(k)}\) is negative, i.e., if the step is downward, the network configuration, represented by \(\alpha _4\) and \(\alpha _5\), approaches the statistics of the target configuration so that the kth switch is accepted. If \(\Delta \vartheta\) is upward, the switch is rejected. Following this procedure may not lead to satisfying results, as the system may be trapped in a local minimum far away from the statistics of the actual target configuration. Thus, the Metropolis-Hastings condition from Eq. (96) is adapted and, additionally, included so that every uphill step is accepted with the probability:
where the chance of upward acceptance is adjusted by choosing an appropriate \(\bar{T}\). In our simulations, we found that a value between \(10^{-3}\) and \(10^{-4}\) for the quantity \(\bar{T}\) leads to satisfying results. The network configuration after 200 switches, obtained by controlling the ring neighborhood statistics using the Aboav-Weaire parameter \(\alpha\), is depicted in Fig. 24b. One observes the even distribution of small and large rings of this artificial network structure and its similarity to the real measured configuration in Fig. 21c, while the configuration in Fig. 24a does not represent a physically meaningful network structure when compared to the measurement.
The dual lattice formulation has two significant advantages. Firstly, it allows one to perform a switching move more accurately since switches in an amorphous network structure are generally not rotations of precisely 90 degrees, but they may vary depending on the level of disorder of the system. Secondly, it allows performing a certain number of switches on the dual lattice only without relaxing the physical lattice using the Yukawa-type potential. On the one hand, this considerably increases the efficiency of the algorithm, allowing one to generate significantly larger samples, which is a necessity for meaningful mechanical investigations. On the other hand, one may argue that the procedure is, strictly speaking, not thermodynamically consistent. However, the physical check comes from the minimization of the physical lattice a priori to the switching sequence on the dual lattice and, additionally, from the Aboav-Weaire law that must be satisfied after every switch.
To also control the level of heterogeneity of the network, the overall statistics of rings must be considered in an additional objective function \(\vartheta _p\). The microscopic pictures in Fig. 21b and c reveal ring sizes between four and nine. We present the ring statistics, extracted from this figure, in the gray histogram in Fig. 25a.
a Experimentally extracted ring statistics from Fig. 21c and the corresponding log-normal fit; b tuning the extracted ring statistics—from crystalline to vitreous; c generated sample with a low target heterogeneity after 2000 switching attempts; d generated sample with a higher target heterogeneity after 2000 switching attempts
The peak of the histogram is six, i.e., rings with six members are found most frequently in the network, which is the energetically favorable configuration. Furthermore, one observes more five-membered than seven-membered rings, which results in an asymmetric histogram, which may be fitted best by a logarithmic distribution function, as also discussed by Büchner et al. [107]. The logarithmic distribution, as presented by the black line in Fig. 25a, follows the form:
where \(x_n\) is the continuous counterpart to the ring size integer value. Moreover, \(\mu ^{(m)}\) and \(\sigma ^{(m)}\) denote the measured target mean and target standard deviation that can be extracted from the experimental topology from Fig. 21c. The level of the target heterogeneity of the network can then be manipulated by changing the value of the extracted standard deviation of the distribution. A higher standard deviation leads to an increase in small- and large-membered rings and a decrease in the six-membered rings, as demonstrated in Fig. 25b. After tuning the continuous log-normal distribution, artificial asymmetric histograms \(p_n^{(t)}\) are extracted that serve as target distributions for the ring statistic objective function:
where the magnitude \(c_3\) serves as a weighting factor, chosen as 1.0 for the calculations in this paper. The higher \(c_1\) and \(c_2\) compared to the factor \(c_3\), the more emphasis on the ring neighborhood statistics during the switching sequence. Vice-versa, the higher \(c_3\) compared to the \(c_1\) and \(c_2\), the more emphasis on the overall ring statistics. Clearly, unphysical choices that significantly deviate from the chosen factors here may result in non-converging network configurations.
The total objective function to be considered during the dual Monte Carlo bond switching algorithm is then written as:
Figure 25c, d are the result of \(2\times 10^3\) dual switching attempts using target Aboav-Weaire parameters of 0.3 for both \(\alpha _4\) and \(\alpha _5\). However, to generate the sample in Fig. 25c, we used a target variance \(0.6\,\sigma ^{(m)}\), i.e., a lower target heterogeneity, while to generate the sample in Fig. 25d, we used a target variance of \(1.4\,\sigma ^{(m)}\), i.e., a higher target heterogeneity compared to the value extracted from the measurement. We note that it was necessary to relax the physical network after every tenth successful dual switch, to check if the network topology is stable and the configuration remains in accordance with the experimental observations and Zachariasen’s requirements.
Notably, the number of experimentally imaged 2D silica samples to be found in the literature is rather limited. Thus, an experimental extraction of different levels of heterogeneity is not possible to date. However, the bond switching algorithm shows that the networks smaller than \(0.6\,\sigma ^{(m)}\) and larger than \(1.4\,\sigma ^{(m)}\) do not converge regarding the objective function \(\vartheta\) in Eq. (104), which gives us an idea about the network topologies that are physically meaningful.
9 Mechanical Deformation
Born and Kármán [151] already pointed out in their pioneering paper that the vibrational properties govern the mechanical response behavior of a particular molecular structure. Indeed, numerical model structures employing inter-particle interactions, as discussed in Sect. 2.1, result in complex statistical atomic vibrations and diffuse response trajectories. However, it may not always be necessary to follow these noisy statistical trajectories to understand the underlying deformation mechanics of the system. Athermal, quasistatic strategies have been shown to clearly demonstrate the mechanical response behavior of disordered solids if the temperature range of investigation is far enough away from the transition temperature \(T_g\).
Consequently, in order to investigate the pure structure-property relationship of materials, it has become popular to study materials at the athermal limit of zero temperature [80, 81]. In this situation, the system is overdamped, the velocity of each particle vanishes, and the Newtonian set of equations of motion (13) results in a static equilibrium governed by a highly nonlinear PEL:
The sum of all interatomic forces \(\varvec{f}_{ij}\) vanishes for every atom i. The configuration stands still in a local minimum of the PEL, and no thermal vibration is considered.
We emphasize that such a situation may be approximated by performing isothermal and isothermal-isobaric simulations at very low temperatures [79]. In doing so, one controls the cell size and shape while performing the chosen time integration scheme. The cell borders are deformed incrementally, and the atomic configuration is subjected to the corresponding affine displacement field during every integration step. The deformation velocity should be as low as possible to avoid deviating too much from a macroscopic deformation rate applicable within a laboratory. Considering the integration time step lengths, e.g., for vitreous silica 0.5 fs using the Vashishta- or Born–Meyer–Hudgins-type of potential, this seems to be an unrealizable task from an engineering perspective. An accepted strategy in this context seems to be an adaptive procedure whereby the deformation velocity is step-wisely decreased until no further change or feasible response quantities are observed. Furthermore, in the case of atomic rearrangements, i.e., sudden topological changes, such as nano cracks, unwanted dynamic wave effects occur, during which the temperature temporally increases and vibration energy travels through the material. These dynamic effects complicate the investigation of the structure-property relationship of the material.
The problems above do not occur when using the athermal quasistatic (AQS) deformation protocol. Instead of a time integration scheme, one uses an athermal approach applying a sequence of incrementally increasing the strain and then minimizing the potential energy. Thus, local basins can only be escaped by externally induced strain, not by thermal activation. This approach is a physically meaningful description of glassy structures concerning temperatures much lower than the glass transition temperature \(T_g\). Therefore, the AQS protocol allows a pure mechanical description of the disordered structure without thermal disturbances.
During one AQS step, the shape of the PEL slightly alters. Adding subsequent AQS steps, one obtains an entire AQS deformation protocol during which the configuration, indicated by the red circle in Fig. 26b, remains continuously in a local basin of the PEL.
Qualitative illustration of the mechanical deformation-induced change in the PEL
Since external deformation alters the shape of the PEL, a point may appear that allows the system to drop down into an adjacent basin, as indicated on the right-hand side of Fig. 26. As a result of this drop, the system experiences a catastrophic, discontinuous rearrangement event.
9.1 Linearization
Using the AQS deformation protocol, the system remains in a local minimum of the potential energy, except when the minimum disappears and the system drops into an adjacent basin. However, for this subsection, we consider the configuration to be in a fixed deformation state so that \(\varvec{H}=\hat{\varvec{H}}\). In this way, one can define the atomic coordinates referring to this local energy minimum as reference positions \(\varvec{r}=\hat{\varvec{r}}\).
Linearization of the molecular system around a minimum of the potential energy landscape; a illustrative cutout of the PEL with one configurational coordinate (black line) and the linearization around the basin of the PEL (blue line); linearized oscillation (blue ball and blue dashed lines) versus the nonlinear response (blacked dashed line); b practical implementation of the first order force approximation of the potential energy basin—the Hessian. One recognizes the severe deviation from the second order approximation, especially, far away from the minimum in a
Any change in atomic position that results in the system moving away from this local minimum is defined by the concatenated displacement vector \(\varvec{u}\). Brownian motion that does not lead to any activation such that a \(\beta\)- or even an \(\alpha\)-relaxation process occurs, as shown in Fig. 13, may be interpreted as a nonlinear oscillation in \(\varvec{u}\) around the atomic position of the local minimum \(\varvec{r}\), as indicated in the one-dimensional illustration of Fig. 27a. However, as also indicated in this figure, if the displacement exceeds a critical value, the system may jump over the adjacent maximum and drop down into a nearby minimum of the PEL.
In the concatenated coordinates, the set of equations of motion (13) is written in matrix form as:
where \(\varvec{M}=\text {diag}\left[ m_k \right] \in \mathbb {R}^{3N\times 3N}\), and \(\varvec{f}\in \mathbb {R}^{3N}\) denotes the concatenated force that corresponds to the displacement vector \(\varvec{u}\).
The Taylor expansion of the PEL around the minimum \(\varvec{r}\) is written as:
Neglecting all terms higher than order two and considering that the gradient in the local basin \(\nabla _{\varvec{r}}\mathcal {U}\left( \varvec{r},\varvec{H} \right)\) vanishes leads to a second order approximation of the PEL:
where \(\mathcal {U}^{(2)}\left( \varvec{u},\varvec{H} \right) \approx \mathcal {U}\left( \varvec{r},\varvec{H} \right)\) for small displacements \(\varvec{u}\). Consequently, the increase in potential energy \(\Delta \mathcal {U}^{(2)} \left( \varvec{u},\varvec{H} \right) = \mathcal {U}^{(2)}\left( \varvec{u},\varvec{H} \right) - \mathcal {U}\left( \varvec{r} \right)\) is represented by a quadratic PEL around the local minimum at \(\varvec{r}\), as indicated by the blue function in Fig. 27a, and depends only on the displacement coordinate \(\varvec{u}\). One realizes that for small oscillations, the linearized system behaves very similarly to the nonlinear system, while the error becomes significant in the case of large displacements, particularly in case of possible jumps into adjacent basins of the PEL. The second order derivatives at the configuration \(\varvec{r}\) of the PEL \(\frac{\partial ^2\mathcal {U}\left( \varvec{r},\varvec{H} \right) }{\partial \varvec{r}^2}\) are referred to as the Hessian \(\varvec{\mathcal {H}}\in \mathbb {R}^{3N\times 3N}\), which is written as:
Eq. (109) reveals that the Hessian is symmetrical with respect to both i and j, as well as \(\alpha\) and \(\beta\).
The concatenated force vector is the gradient of the PEL with respect to the corresponding position vector. However, using the second-order approximation of the PEL, as presented in Eq. (108), one performs a transformation from the position coordinate \(\varvec{r}\) to the displacement coordinate \(\varvec{u}\). Thus, the second-order approximation of the concatenated force can be evaluated by the second-order approximation of the PEL depending on the displacement vector:
Taking into account Eq. (106) and the result from Eq. (110), this leads to the linearized matrix set of equations of motion of the molecular system:
The Hessian matrix containing the second derivatives of the potential energy is then equivalent to the stiffness matrix in structural dynamics, and finite element analysis [152]. This linear set of equation of motion obeys the solution ansatz function \(\varvec{u} = \varvec{\varphi }e^{i\lambda t}\). Inserting the solution function into Eq. (111) leads to the eigenvalue equation:
or, equivalently:
Eqs. (112) and (113) have 3N non-trivial solutions for the eigenvectors \(\varvec{\varphi }_m\) belonging to the corresponding eigenvalues \(\lambda _m\). The matrix \(\varvec{D}=\varvec{M}^{-1}\varvec{\mathcal {H}}\) is referred to as the dynamical matrix and is equal to the Hessian matrix regarding its off-diagonal components. As \(\varvec{M}\) is a diagonal matrix, the inversion is cheaply performed by inverting its diagonal entries and \(\textbf{D}\) remains symmetric implying existence of 3N real eigenvalues (counting multiplicities) that we enumerate in increasing order \(\lambda _1 \le \lambda _2 \le \cdots \le \lambda _{3N}\). We collect the eigenvalues in the spectral matrix \(\mathbf {\Omega }=\text{ diag }\left[ \lambda _k \right]\). Subsequently, one can evaluate the corresponding eigenvectors that assemble the eigenmatrix \(\varvec{\Phi }\), such that \(\varvec{\Phi } = \left[ \varvec{\varphi }_1, \varvec{\varphi }_2, \cdots , \varvec{\varphi }_{3N} \right]\). For our calculations, the eigenvectors are orthonormalized with respect to \(\varvec{M}\). Thus, the following properties hold true:
9.1.1 Implementation of the Hessian
The most straightforward implementation of the Hessian is a subsequent evaluation of pair interaction matrices \(\varvec{\mathcal {H}}_{ij}\), whereby one such interaction of atoms i and j embedded in the configuration is written as:
When evaluating the set of interaction matrices one takes advantage of the symmetry properties, i.e., \(\varvec{\mathcal {H}}_{ij} = \varvec{\mathcal {H}}_{ji}\). For the purpose of illustration, we present the implementation of the interaction matrix for the three-dimensional space:
However, one can take advantage of \(\varvec{f}_i=-\nabla _{\varvec{r}_i}\mathcal {U}\left( \varvec{r},\varvec{H}\right)\) since analytically evaluating the second derivatives of the potential energy can be quite cumbersome. As soon as the atomic forces are evaluated, the Hessian can be expressed as the gradient of the atomic force vector \(\varvec{\mathcal {H}}_{ij}=-\nabla _{\varvec{r}_j}\varvec{f}_i\), which is implemented as:
Based on this approach, one may implement the Hessian using the finite difference approach, as lately employed by Bonfanti et al. [153]. As also visualized by a two-dimensional representation in Fig. 27c, one locates a pair of atoms i and j. Considering the symmetry properties, it is sufficient to choose j, always greater or equal to i. The purpose of the interaction matrix is to measure the change in force on the atom i caused by small position changes of the atom j while keeping all other atomic positions frozen. As an example, we demonstrate the component in the first row and second column \(\frac{\partial f_{i,1}}{\partial x_{j,2}}\) of \(\varvec{\mathcal {H}}_{ij}\). Firstly, atom j is displaced by a small value \(\delta\) into the positive \(x_2\)-direction. The interatomic distances are updated, and the \(x_1\)-component of the force on atom i is evaluated \(f_{i,1}^{+x_2}\). Secondly, atom j is displaced by \(2\delta\) in the negative \(x_2\)-direction. Then, the interatomic distances are again updated, and the \(x_1\)-component of the force on atom i is evaluated \(f_{i,1}^{-x_2}\). The component of the pair interaction matrix is then approximated by:
The value \(\delta\) should neither be fixed too small, so that one does not fall below the level of machine accuracy, nor too large so that the values of the Hessian become inaccurate. For the examples in this paper, we observed that, choosing a value of \(10^{-6}\) for \(\delta\), leads to accurate and stable numerical results. The eigenvalues \(\lambda _k\) and corresponding eigenvectors \(\varvec{\varphi }_k\) are evaluated using the Lanczos algorithm, which is available in the majority of computer libraries [154]. In many cases, it is sufficient to evaluate only a smaller set of the lowest eigenvalues and corresponding eigenvectors, as will be discussed in the following sections. In Fig. 28, we present the visualization of the lowest eigenvectors of the Lennard-Jones glass sample, depicted in Fig. 17c.
Linearization around the relaxed structure; eigenvalue decomposition of the Hessian; plots a–f represent mode 1 to mode 6, respectively. We note that the first two modes also called the Goldstone modes, represent the rigid body motion due to the first two zero eigenvalues of periodic systems in two dimensions
The Hessian of a system with periodic boundary conditions is positive semi-definite [155]. Thus, the first two lowest eigenvalues in a two-dimensional system are zero. The two corresponding eigenvectors are referred to as the Goldstone modes [156], and are depicted in Fig. 28a, b. These modes derive from the fact that the sample can be translated in two linearly independent directions in two-dimensional systems with periodic boundary conditions without subjecting any forces to the system. Thus, the atoms leave the box at one end while entering on the other, as discussed in Sect. 3. The next four eigenvectors, depicted in Fig. 28c–f, represent the four lowest modes that correspond to the vibration patterns which assemble the global oscillation movement of the linearized system (111). One may interpret the eigenvalue \(\lambda _k\) as a pseudo force magnitude, necessary to deform the system into the corresponding eigenmode state, cf. Eq. (112). This perception already underlines the importance of eigenmodes during deformation. They provide information about the soft patterns of material deformation and are, therefore, an indicator of material response behavior—an issue we will discuss in more detail in the following section.
9.2 Shear Deformation Response
We start our investigations with the AQS protocol of simple shear deformation on two illustrative two-dimensional numerical examples to show the fundamental mechanical differences between ordered and disordered solids. We illustrate the rich dynamics of deformed glassy structures using the second example, and we refer to the relevant literature.
During simple shear, the cell shape, described by the Bravais tensor \(\varvec{H}:=\varvec{H}\left( \gamma \right)\), depends only on one scalar parameter \(\gamma\). Recall that, for the sake of abbreviation, we define the shear angle as \(\gamma :=\gamma _{12}\). Consequently, the Bravais tensor in the reference configuration is defined as \(\hat{\varvec{H}}:=\varvec{H}\left( \hat{\gamma }\right)\), and the potential energy also depends only on the atomic position and one scalar parameter describing the shape change, i.e., the shear angle:
Before discussing the shear deformation mechanics of disordered materials in detail, we summarize and refine the geometrical background and the required transformation properties. Figure 29 presents a two-dimensional illustration of the configuration in the reference domain \(\hat{\Omega }\) and the respective configuration in the current domain \(\Omega _\gamma\).
Simple shear mappings and configurations
The deformation gradient for simple shear maps a reference configuration in the domain to a current configuration \(\varvec{F}_{\gamma }:\hat{\Omega }\rightarrow \Omega _\gamma\), as discussed in Sect. 4. In Fig. 29, we present one particle during three different states of interest, two states in the reference domain on the left-hand side and one state in the current domain on the right-hand side. The first state is the initial position prior to the subjected mechanical loading \(\hat{\varvec{r}}_{0i}\in \hat{\Omega }\). The second state \(\varvec{r}_i\) is the total outcome of the affine mapping \(\varvec{F}_{\gamma }\) applied to \(\varvec{r}_{0i}\) and the additional material response that deviates from the affine mapping. This additional material response is effectively caused by the displacements resulting from the minimization procedure of the PEL, also referred to as non-affine displacement. In the following section, we will, among others, present an example where the non-affine shear response is equal to zero and another example where the non-affine shear response is not equal to zero. The third state is the current atomic position \(\varvec{r}_i\) mapped back to the reference domain using the inverse mapping of \(\varvec{F}_{\gamma }\). As we will see later, this “artificial” state allows us to abandon the affine motion of the particle and study the non-affine part only. We summarize the transformation properties of simple shear for one particle as:
Equivalently to the definition in Eq. (2), all particle positions of the three states can conveniently be described by \(\hat{\varvec{r}}_0\), \(\varvec{r}\) and \(\hat{\varvec{r}}\), regarding state one, two and three respectively. To map concatenated quantities we define the global deformation gradient \(\bar{\varvec{F}}_\gamma : \hat{\Omega }^N \rightarrow \Omega ^N\) so that it is written as follows:
This mapping can be easily realized so that it is consistent with the definition of the concatenated coordinate by defining a block-wise diagonal matrix, written as:
where all off-diagonal block elements remain zero.
9.2.1 Simple Shear Deformation of a Hexagonally Packed Lattice
The first example shows the simple shear response of a simple crystalline structure. Using only one particle type, we modeled the system by employing the Lennard-Jones potential, as presented in Eq. (29). We placed 20 rows containing 20 atoms in a box with periodic boundary conditions to prepare the structure. The distance between every atom and its nearest neighbor equals the equilibrium distance. All forces vanish in this configuration since the system is located in a global minimum of the PEL. The result is a perfectly ordered, hexagonally packed lattice, as shown in Fig. 30a.
Deformation mechanism of a crystalline solid; a crystalline structure of a hexagonally packed model material implemented by the Lennard-Jones potential, as presented in Eq. (29); b sample box subjected to a linearly increasing shear deformation protocol; c plastic shear slip event; d non-affine displacement field of the plastic slip event; e stress-strain curve of the linearly increasing shear deformation protocol with the rearrangement event indicated by the corresponding stress-drop in red; the red particle in a–c is subjected to a numerical disturbance by displacing it horizontally to the right by a value of \(1.0e-6\), which allows the nucleation of a clean inelastic stress drop event
The mechanical shear behavior of such a defectless system is well known, as elaborately discussed in the literature [157], and ends up in a defectless flow effect for higher strains. However, as every atom has the same neighboring environment, we introduced an artificial numerical defect at the position of one particle by horizontal displacing it by a small value. The corresponding particle is highlighted in red in Fig. 30a. This small modification allows one to locally initiate deformation-induced atomic rearrangements within the crystalline environment.
Using the simple shear deformation protocol, we subject the linear transformation mapping to the box boundaries and the corresponding affine deformation to every particle. However, the response of the material is evaluated by:
which is equal to the affine mapping. Thus, we emphasize that the non-affine elastic response of the hexagonally packed structure to simple shear deformation (124) is zero. The entire chosen shear deformation protocol for this example, is a linear step-wise increase in shear strain using the transformation mapping in Eq. (124), applying a step size of \(\Delta \gamma = 10^{-3}\). For low shear angles \(\gamma\), one observes a linear, reversible stress-strain relation. We present the deformed shear state at a relatively small shear strain value in Fig. 30b and the corresponding stress-strain curve with the perfectly linear behavior for the initial strain regime in Fig. 30e. It is noteworthy that, in the case of a perfect defect-free configuration, the AQS protocol may not be the best choice since the total deformation of each particle coincides with the affine field of the simple shear deformation as long as the system response remains reversible.
To locally initiate a possible clean catastrophic event, one introduces a simple artificial dislocation in the form of a displacement disturbance to one particle in the system before every AQS step. As mentioned above, we horizontally disturbed the tenth particle in the tenth row by a value of \(10^{-6}\), highlighted in red in Fig. 30a–c. The numerical disturbance is sufficiently small not to noticeably alter the mechanical response behavior but large enough to favor the jump over one particular possible saddle-node bifurcation that may emerge in the PEL of the system when increasing the shear strain. Thus, during every elastic AQS step, the task of the minimization algorithm is to push the disturbing atom back into the deformed crystalline formation. After the initial elastic strain behavior, one observes a distinct, sudden stress drop, as highlighted by the red line in Fig. 30e. This is accompanied by a microscopic rearrangement event, as shown in Fig. 30c. To emphasize this structural rearrangement, we did not project back any atoms that left the box boundaries after this AQS step to their corresponding duplicate. One observes that the row that holds the disturbing atom slips abruptly into the shearing direction by one atomic pair distance. As shown in Fig. 30d, the non-affine displacement field does not vanish during this event. The non-affine displacement vanishes exactly half the box length away from the slipping line and linearly increases towards the slipping line from both sides. After the slip, the structure of the material jumps topologically back into its original configuration. Therefore, neither volume nor density changes locally after the event. Energy dissipates, and the stress-strain behavior remains parallel to the initial stress-strain line until the next stress drop occurs. Thus, the shear stress-strain response of the crystalline material is a recurrence of equidistant stress drops that correspond to catastrophic slip events and parallel elastic branches in between them. The interatomic distance determines the size of one elementary slip event. Thus, one may easily imagine that an increase in the number of atoms leads to an effect whereby the slip event becomes smaller compared to the sample size. The size of the stress drops decreases so that, in the asymptotic limit, one obtains a linear elastic and perfectly plastic shear stress behavior following the assumptions proposed by basic continuum-mechanical material descriptions [158]. Overall, inelastic shear rearrangements in crystalline solids emerge from local dislocations or material defects that govern the plastic flow behavior [80, 159].
9.2.2 Simple Shear Deformation of a Lennard-Jones Glass Structure
We start the investigation of glassy materials with a numerical illustration, with which we will demonstrate the progress of the relevant literature. The numerical example presents the shear response of a disordered structure, represented by a binary Lennard-Jones glass. The numerical generation of this material is, amongst other glass types, elaborately discussed in Sect. 8.1. Starting from a physically possible random configuration of 498 large and 402 small particles, we let the configuration relax at a temperature significantly higher than \(T_m\) for \(10^6\) integration steps using the NVT ensemble. Subsequently, the system is quenched to zero temperature infinitely fast, to obtain its corresponding inherent structure, as shown in Fig. 31a. As discussed in Sect. 8.1, the system does not have the time to rearrange into a crystalline structure. Additionally, we note that this material is also inherently frustrated, so that crystallization is not possible [78, 160, 161].
Deformation mechanism of a disordered solid; a generated sample of a disordered material, represented by a binary Lennard Jones glass; b sample box subjected to a linearly increasing shear deformation protocol; c stress-strain relation of the linearly increasing shear deformation protocol; d enlargement of a local region before the shear rearrangement; e enlargement of the same local region after the shear rearrangement
Equivalently to the crystalline example in Sect. 9.2.1, we subjected the configuration to a stepwise linear increase in shear deformation using the AQS protocol, as depicted in Fig. 31b, and measured the stress response using Eq. (58). However, we note that the material response to the shear deformation map (124) is not equal to the affine motion only. We also expect a deviation from the affine term, written as:
The deviation of the affine deformation response caused by the structural complexity of the material \(\varvec{x}_i = \left( x_{i1}, x_{i2} \right) ^T\) is referred to as a non-affine part of deformation. However, one also recognizes the additive split of affine and non-affine displacements. Later in this section, we will see that the picture of non-affine motion is vigorously linked with the rich and complex behavior of disordered materials.
The stress-strain response of the example at hand is presented in Fig. 31c. While the stress response of the crystalline structure seems relatively straightforward, the shear stress-strain curve of the disordered structure reveals a highly complex, even chaotic, response behavior. As a result of this, we substantially introduce the definition by Jin et al. [162], who classified the shear response of disordered solids and glasses into three states:
-
(i)
For small loading, the system behaves purely elastic, where the stress increases smoothly with the strain—the elastic regime;
-
(ii)
With increasing \(\gamma\), the stress-strain curve becomes jerky. Every incremental increase in strain may result in a stress drop—the marginally stable regime; and
-
(iii)
Further increase in \(\gamma\) results in yielding, which manifests in a jerky flow consisting of a sequence of elastic branches irregularly intersected by stress drops that strongly vary in their magnitudes. During this state, the system experiences structural rearrangements that may cover large parts of the cell volume—the limit of existence of the solid.
In his pioneering work, Spaepen [163] demonstrated that fracture occurs along planes of shear bands, in which the material softens locally and deforms much faster than in the rest of the sample. Performing a pure shear experiment, as illustrated in Fig. 9b, the shear plane forms at a \(45^{\circ }\) angle with respect to the principle strain axes and is weakened by some shear-induced structural change. Thus, macroscopic flow occurs due to a finite number of individual atomic jumps that are somehow arranged in the shear plane. However, Spaepen argued in his paper that the atoms could only jump if there were enough free atomic volume in the nearest neighborhood. Every jump occurs from one into another minimum in the PEL, with the free volume being the basis for his flow formulation. Overall, these atomic jumps are biased in the direction of the external forces on the system. Linking the free volume with a particular region that experiences a shear-induced structural change may be seen as a pioneering attempt to provide an a priori quantification, i.e., a prediction of such localized events. However, we note that recently it was shown that the local density, which highly correlates with the free volume, is a relatively poor local quantity of prediction that allows one to find spots that are susceptible to rearrangements [164]. In fact, more promising prediction quantities have been proposed in the literature, as we will present later in Sect. 9.3.
Argon and Kuo [165] published another pioneering work. Intriguingly, they performed bubble experiments to elucidate the microscopic deformation picture of metallic glasses. Shortly after, Argon and Shi [166, 167] provided extensive numerical investigations using an inter-bubble potential. Their findings are in accordance with the statements by Spaepen [163], which are that shear-induced deformation results in very localized shear transformations with the size of a few bubbles only, which can qualitatively be translated to shear transformations of disordered solids [168]. In their pioneering work, they reported two types of shear transformations in their bubble experiments. The first type is an eye-catching event, where a local row of a few atoms slides against a parallel situated row along the direction shear plane. They called this type a disk-shaped shear transformation. The second type is a rearrangement along the principal axes of the simple shear deformation, which they call diffusive shear transformation events. Such an event is the local material response governed by the global principal directions of shear strain we elaborated on earlier in Eqs. (40)–(46). As we will discuss later in this paper, the second type is an elementary rearrangement picture for inelastic deformation events in disordered solids. Shear transformations were later experimentally studied in colloidal glass systems [169,170,171]. Intriguingly, the behavior of disordered solids under shear deformation can essentially be captured in such models.
To further investigate shear transformation events, we now inspect one stress drop in more detail, which is highlighted in red in Fig. 31c. As observed in the studies mentioned above, such a stress drop corresponds to a strongly localized rearrangement event that affects one or more groups of a relatively small number of atoms. We denoted the location of such an atom group by the blue window in Fig. 31b and enlarged this window before the stress drop event in Fig. 31d and after the stress-drop event in Fig. 31e. In agreement with the essence of the study of Argon and Kuo [165], one observes the microscopic picture of a sudden local—what they called diffusive—shape change, marked by a red circular area in Fig. 31d, which is biased to deform into the principal direction of shear, marked by the ellipsoidal area in Fig. 31e. In other words, the global simple shear deformation protocol results in a sequence of rearrangements, each occurring around a center of a relatively small number of atoms. The rearrangement center is then biased to being pulled along the first principal axis while pushed inwards along the second principal axis of simple shear deformation. The elastic displacement response around the rearrangement center follows the picture of a quadrupolar shape, and it is in good agreement with the elastic response of an ellipsoidal inclusion presented by Eshelby [172].
Identifying spots in disordered solids that are prone to rearrangements, so-called shear transformation zones [79], turns out to be a challenging task, not least because the highly nonlinear nature of such systems complicates prediction attempts, as discussed later in this section. The question arises whether such local rearrangement spots exist in the material a priori to the external loading protocol and could be seen as microstructural defects, such as dislocations in crystals. Based on the assumption that local spots preexist in the material, constitutive equations were developed that incorporate local entropic fluctuations in the volume of shear transformation zones [14, 79, 173, 174]. The literature did, however, not entirely agree on this matter, as Gendelman et al. [175] demonstrated that shear transformation zones might change when altering the shearing direction. Thus, the discussion regarding the existence of shear transformation zones depends on the configuration and the deformation protocol, which can significantly complicate an a priori identification of such local spots in the material. In this context, Barbot et al. [176] have shown that shear transformation zones possess weak planes. In this way, they show that protocol-dependent plastic activity does not rule out the existence of such zones a priori encoded in the material. Nevertheless, the local stress thresholds are very sensitive to the loading direction.
9.2.3 The Nonlinear Elastic Simple Shear Response
The total displacement field of the crystalline lattice response is equal to the affine displacement field, and the atomic forces are zero directly after the affine deformation step that corresponds to the shear increment \(\Delta \gamma\), as we were subjecting our deformation protocol to a Bravais lattice structure. However, subjecting the disordered configuration, which is not a Bravais crystal, to shear strain, the forces on the system do not vanish. We need a correction so that the system remains in an equilibrium position, i.e., the non-affine correction \(\varvec{x}_i\), as presented in Eq. (125).
Taking into account the definition of the potential energy for simple shear deformation protocols in Eq. (120), one introduces a potential energy measure with respect to the reference domain as follows:
where we considered the mapping from the actual into the reference configuration using the concatenated deformation gradient \(\bar{\varvec{F}}_\gamma\), as defined in Eq. (123).
Since the mapping of the microscopic coordinate back into the reference cell by applying the inverse transformation \(\hat{\varvec{r}} = \varvec{F}_\gamma ^{-1}\varvec{r}\) is equivalent to subtracting the affine part from the actual coordinate we will refer to \(\hat{\varvec{r}}\) as the non-affine coordinate.
Such a division of the mechanics in affine and non-affine parts raises the question of how much disorder influences the overall mechanical elastic response since disorder seems to be somewhat closely related to non-affinity. Looking at the stress-strain response of the shear deformation example of the disordered solid in Fig. 31c, one may expect—at first sight—that these elastic branches may display similar dynamics to the elastic crystalline solid, and the difference may only occur when approaching catastrophic events. However, in what follows, we will show that this is not the case, and we will unravel the dynamics of elasticity of disordered solids based on the findings of Maloney and Lemaître [80,81,82] and Lerner [155].
Before we study the mechanics of disorder in more detail, we note that the derivative of an arbitrary function f dependent on \(\hat{\varvec{r}}\) and \(\gamma\) is written as:
since \(\hat{\varvec{r}}\) itself depends on \(\gamma\). On the right-hand side, one detects the change in the non-affine coordinate with respect to the shear angle gamma \(\frac{d\hat{\varvec{r}}(\gamma )}{d\gamma }\). This term can be interpreted as a non-affine deformation velocity. Firstly, we determine static equilibrium, which means that the concatenated force vector vanishes at every deformation state \(\gamma\):
It is worth recalling at this point that \(\varvec{r}\) and \(\hat{\varvec{r}}\) implicitly depend on \(\gamma\) in Eq. (128). However, we omit this dependence in the majority of the equations for the sake of compactness. Since \(\bar{\varvec{F}}_\gamma ^T\) is invertible, one can define the force mapped back to the reference configuration, which must also vanish, as follows:
During an elastic response segment in the quasistatic framework, the system remains continuously in a minimum of the PEL so that both the force vectors \(\varvec{f}\) and \(\hat{\varvec{f}}\) remain at zero with varying shear angle. For the change in \(\hat{\varvec{f}}\) with varying shear angle, represented in the reference configuration, one obtains the following:
where \(\hat{\varvec{\mathcal {H}}}(\gamma )\) denotes the second order derivatives with respect to the reference coordinates \(\frac{\partial \hat{\mathcal {U}}(\hat{\varvec{r}},\gamma )}{\partial \hat{\varvec{r}}\partial \hat{\varvec{r}}}\). This \(3N\times 3N\) matrix is written as follows:
where it holds true that:
One concludes that \(\hat{\varvec{\mathcal {H}}} \rightarrow \varvec{\mathcal {H}}\) for \(\gamma \rightarrow \hat{\gamma }\) since the deformation gradient becomes the identity tensor. Thus, in this case, the derivatives with respect to the non-affine coordinates are equivalent to the derivatives with respect to the microscopic coordinates. The first term in Eq. (130) represents the change in the force in the configuration if only affine deformation is applied and the system is not yet corrected by the non-affine part. Thus, we define this non-affine correction force as:
This correction force may also be interpreted as a virtual force on the system due to an affine disturbance. The change in the force due to the affine deformation, \(\tilde{\varvec{f}}(\gamma )\), can therefore be used as a mechanical measure for the level of disorder in a system. Thus, the second part of Eq. (130) is responsible so that the system remains in a state of zero force. Moving the second part in Eq. (130) to the right-hand side, the mechanical force equilibrium of disordered elasticity is written as:
which reveals that the change in force on the system caused by disorder is equal to the Hessian, i.e., a tangential stiffness operator depending on the shear strain \(\gamma\), mapped to the reference configuration and contracted by the non-affine velocity field. In the limit of \(\gamma \rightarrow \hat{\gamma }\) the mechanical force equilibrium of disordered elasticity (134) simplifies to:
When increasing the strain step-wisely in finite shear strain increments, the non-affine velocity field \(\frac{d\hat{\varvec{r}}(\hat{\gamma })}{d\gamma }\) leads to a certain finite displacement after each deformation step. Thus, in this discrete context, one may also call it a non-affine displacement field. Consequently, the representation of the deformation mechanics in the reference configuration allows one to represent the pure non-affine motion only. In order to invert the Hessian in Eq. (134), we have to remove the eigenmodes that belong to the zero eigenvalues, i.e., the Goldstone modes describing rigid body motion (134), as shown in Fig. 28a, b.
After discussing the equation of disordered microscopic motion, we aim to determine the influence of disorder on some macroscopic measures of interest. In doing so, the shear stress is, by definition, evaluated as the change in potential energy with the shear strain. Considering the derivation law in Eq. (127), the shear stress is evaluated as:
The term \(\nabla _{\varvec{r}}\mathcal {U}\left( \varvec{r},\gamma \right)\) vanishes as the system is continuously located in a minimum of the PEL where it remains in a state of zero forces. Thus, we cannot extract any further information regarding the influence of structural disorder on the stress-strain behavior. However, one may differentiate this relation once more with respect to the shear angle to obtain the shear modulus, which is, by definition, the change in the shear stress with respect to the shear strain:
where, in the last two lines, we considered the force equilibrium for disordered elasticity, as presented in Eqs. (133) and (134). Thus, in the limit of \(\gamma \rightarrow \hat{\gamma }\) the shear modulus is written as:
The Hessian is, by definition, positive semidefinite. Thus, the second term in Eq. (137) will always be greater or equal to zero, depending on the level of disorder. This property reveals that disorder will always reduce and never increase the magnitude of the shear modulus. Thus, Zaccone et al. [177] developed an analytical constitutive stress-strain relation for metallic glasses based on the non-affine displacement fields. Their model uses only a few fitting parameters while considering the microscopic mechanics.
After disentangling the elastic theory of the mechanics of amorphous solids, we present the non-affine displacement field of the elastic response at a relatively small shear angle \(\gamma\) in Fig. 32a. Applying the transformation in Eq. (31) to the configuration, we map the atomic positions to the unit cell via \(\varvec{\textsf{r}}_i = \varvec{H}^{-1}\hat{\varvec{r}}_i\), to illustrate the non-affine displacement vector field in a unit square, as depicted in Fig. 32a. Inspection of this vector field shows again the complexity of the elastic deformation procedure of disordered materials, which reminds us of the picture of vortex-like structures [178, 179].
Non-affine displacement fields during a an elastic deformation step and during b an inelastic rearrangement event, i.e., a shear transformation (STZ). The red circle indicates the rearrangement center
9.2.4 Approaching Catastrophes
As the curvature of the potential energy landscape becomes flatter and flatter, as indicated by the cartoon plot in Fig. 26, the second term in Eq. (137) increases, while the shear stiffness in one or more localized areas in the material continuously decreases so that the shear modulus, as presented in Eq. (137), decreases and, eventually, becomes negative. The system destabilizes in one direction. Thus, at this point, the system behaves very distinctively and very localized in terms of one catastrophic displacement shape [82]. The non-affine displacement field suddenly jumps by several magnitudes. In short, the microscopic configuration experiences a sudden change when the local minimum transforms into a first-order saddle node, as depicted in Fig. 26. Figure 32b shows a typical picture of the non-affine displacement field of an elementary catastrophic event taken from the numerical example, presented in Fig. 31. For the purpose of visualization, we used a higher scaling factor for Fig. 31b than for Fig. 31a. Qualitatively in accordance with the microscopic changes of another event, microscopically captured in Fig. 31d, e, one observes the displacement field of elementary shear. In accordance with Maloney and Lemaitre [82], we defined the center of rearrangement by identifying the position of the atom that experiences the maximum displacement, indicated by a red circle in Fig. 32b. From a qualitative inspection of this figure, one can say that the statistical expectation of the deformation picture of such events is that the material is pulled into the first principal direction \(\varvec{\lambda }_1\) while being pushed into the second principal direction \(\varvec{\lambda }_2\) of simple shear deformation. For the evaluation of the principal strain directions of the simple shear deformation mode, we refer to Eqs. (40)–(46). The nature of disordered materials is responsible for single events deviating from the principal directions of elementary shear dependent on the particular local situation and its microscopic surroundings.
We recall the concept of modal analysis, which describes motion in terms of a superposition of structurally dependent eigenvectors. Interpreting a local spot that will be the center of a catastrophic rearrangement as “soft”, as it experiences a distinct, sudden local decrease in shear stiffness, motivates one to monitor this situation using the eigenvalue analysis. Such soft spots are responsible for vibration modes that need a very small amount of energy to be activated, which is the case for the vibration modes belonging to the lowest set of eigenvalues. Thus, the motivation is to inspect the change in the eigenvalues and corresponding eigenmodes when approaching a sudden local rearrangement. First, we define the dynamical matrix in the reference configuration as:
Taking into account that \(\lambda _k = \varvec{\varphi }_k^T\hat{\varvec{D}}\varvec{\varphi }_k\) and the derivation law from Eq. (127), the change in the kth eigenvalue with respect to shear deformation is written as:
Due to normality conditions, only the first term in this equation remains different from zero. Applying the derivation law from Eq. (127) and considering the equilibrium equation for disordered systems (134) one may rewrite Eq. (140) as:
where \(\odot\) refers to a contraction of a third-order tensor into a second-order tensor, represented by a matrix, so that the result of this equation is a scalar. The inverse of the Hessian is now rewritten in its eigendecomposed form as:
We note that we have only \(2N-2\) degrees of freedom in two dimensions here instead of 2N since we did not consider the two Goldstone modes belonging to the first two zero eigenvalues in this sum. We further note that \(\tilde{\lambda }_p\) and \(\tilde{\varvec{\varphi }}_p\) are the pth eigenvalue and eigenvector of the Hessian, which differ from the eigenvectors and eigenvalues of the dynamical matrix if the masses of the particles are not all set to one. Inserting relation (134) into Eq. (141) and taking into consideration the eigendecomposition of the inverse of the Hessian, one can rewrite the change of the eigenvalues with respect to the shear angle (141) as:
The first thing that comes to mind is the eigenvalue \(\tilde{\lambda }_p\) in the denominator. Only the terms in this sum with a small denominator are relevant, i.e., the terms containing the lowest eigenvalues. Therefore, near a catastrophic event, the system falls into a motion pattern mostly governed by the eigenvectors belonging to the lowest eigenvalues. Thus, we consider Eq. (143) only for the first non-zero eigenvalue, i.e., in the case of two dimensions, \(k=p=3\). The first summation term, i.e., \(p=3\), in Eq. (143) becomes more and more prominent when the strain \(\gamma\) approaches the strain value of the rearrangement \(\gamma _c\) (\(\tilde{\lambda }_3\ll 1\)):
where the scalar c is evaluated as:
This differential equation (144) can be trivially solved. The missing boundary condition comes from the assumption that no additional energy at the rearrangement strain \(\gamma _c\) is required to push the system over the saddle-node: \(\lambda _3(\gamma =\gamma _c) = 0\). Consequently, the solution of Eq. (144) is written as:
We summarize the dynamical behavior when the system approaches a saddle-node as follows:
Considering Eq. (133) and relation (148), one may approximate the shear modulus when approaching \(\gamma _c\):
leading to an evolution of:
as \(\gamma \rightarrow \gamma _c\). Equivalently, the shear stress evolution takes the form of:
as \(\gamma \rightarrow \gamma _c\).
Intriguingly, the dynamics above provide general statements regarding the mechanics of disordered solids and are mainly independent of particular potential functions that may be used for numeric implementations and demonstrations. Numerical model systems, such as Lennard-Jones glasses, do not aim to capture a real material quantitatively. Nevertheless, they are very well able to represent the underlying dynamics of disordered materials. Thus, to underline the mathematical model of a disordered system approaching a catastrophic event, we now carefully inspect the system’s dynamics shortly before the particular event, shown in Fig. 32b. This requires us to decrease the deformation step size for the shear increment \(\Delta \gamma\) from a value of \(1.667\times 10^{-5}\) to \(3.333\times 10^{-7}\) and to also decrease the accuracy threshold \(\varvec{f}^T\varvec{f}\) for the conjugate gradient algorithm from \(10^{-7}\) to \(10^{-9}\). A visualization of the stress-strain curve, enlarged before the rearrangement from Fig. 32b, is presented in the top plot of Fig. 33.
Approaching a catastrophe; a stress-strain behavior during softening; b evolution of the first five non-zero eigenvalues of the Hessian approaching the catastrophe
Further enlarging the peak of this curve reveals a progressive decrease in the slope of the stress-strain curve a few steps before the actual event, which also corresponds to the decrease in the shear modulus, as predicted in Eq. (150). During the step directly before the event, the material’s resistance is almost zero as the slope becomes horizontal, and the curvature in the local basin into the rearrangement direction is almost zero. To quantify the curvature, we evaluated the Hessian matrix \(\varvec{H}\) in the range of 250 shear increments before and 250 shear increments after the event. The evolutions of the first five non-zero eigenvalues during this deformation range are depicted in the bottom plot of Fig. 33. As discussed in detail in Eqs. (140)–(146), one clearly observes that the first eigenvalue continuously decreases long before it can be detected investigating the stress-strain curve and, finally, the system reveals singular behavior and approaches zero shortly before the rearrangement event in terms of a square-root singularity \(\lambda _1 \propto \sqrt{\gamma _c - \gamma }\) (\(\gamma < \gamma _c\)). [81, 82, 180]. Therefore, an observation of the eigenvalue is a much better indicator for anticipating catastrophic events in disordered materials than only an inspection of the stress-strain relationship. From the bottom plot of Fig. 33, it becomes apparent that the first mode mainly governs the dynamics of the upcoming rearrangement. This phenomenon confirms our assumption: soft spots are governed by “soft” eigenmodes, whose corresponding eigenvalues follow a square-root singular behavior when approaching the rearrangement strain.
Linearization directly before the stress drop; a–f linear modes sorted from one to six, respectively (the red mark indicates the rearrangement location); note that the first two modes describe rigid body motion with their corresponding zero eigenvalues
This fact also becomes clear when studying the corresponding microscopic deformation picture of the first six modes one step before rearrangement. For the sake of consistency, we also mapped the mode deformation pictures to the unit box, as defined in Eq. (31). The first two eigenvectors in Fig. 34a, b refer to the Goldstone modes, which describe rigid body motion without any material deformation. Both Goldstone modes are extracted from the Hessian in such a manner that an inversion is possible. Their corresponding eigenvalues are zero, and they are not presented in the bottom plots of Fig. 33. Much more interesting is the third mode, which corresponds to the third eigenvalue that eventually approaches zero in terms of the square root singularity. Comparison of the mode deformation picture in Fig. 34c with the non-affine displacements of the subsequent rearrangement event in Fig. 32b reveals the strong correlation already evident from a qualitative inspection. However, one observes slight differences caused by the fact that, in this particular case, the second deformation mode, i.e., eigenvector four, including the Goldstone modes, shown in Fig. 34d, seems to have a small contribution to the rearrangement as well. This effect is also confirmed when observing the evolution of the corresponding eigenvalue \(\lambda _4\) at the bottom plot of Fig. 33, which noticeably decreases when approaching the strain state of rearrangement. The displacement pictures of modes five and six are depicted in Fig. 34e, f. However, inspecting the evolution of the corresponding eigenvalues, they do not noticeably contribute to the dynamics of the future rearrangement.
9.2.5 From Elementary Inelastic Events to Elementary Shear Bands
Visual inspection of the rearrangement events in the example in Fig. 31 shows that not all events are directly recognizable as a direct response to elementary shear, such as the rearrangement in Fig. 32b already discussed. In particular, the non-affine deformation picture of larger events, such as the red-marked stress drop in Fig. 31c, leads to a more complex superposition of spatially distributed rearrangements. However, it was shown that separating these superpositions is rarely possible [80]. Once the saddle-node bifurcation emerges, the system drops down, and a local rearrangement induces a cascade of one or more subsequent rearrangements [80]. We present a simple approach, initially introduced by Maloney and Lemaitre [80], that separates the single rearrangements appearing during a cascade, i.e., during a particularly large event. Following their strategy, we perform one AQS step directly before the red stress drop in Fig. 31c, indicated by the first black circle. Instead of the conjugate gradient algorithm, we use the steepest descent method, discussed in Sect. 7.1, to find the next basin. Interpreting the progress of the steepest descent algorithm as time, the time derivative of the potential energy is then precisely the sum of the squares of the forces in the system. In other words, the higher the forces during the algorithm, the more substantial the progress made toward the next energy basin.
Investigation of subsequent subevents by following the minimization step history of a catastrophic cascade event
We present the minimization progress using the steepest descent method in Fig. 35. The top plot of this figure shows the evolution of the potential energy, while the bottom plot reveals the corresponding force evolution. This figure shows that the potential energy does not decrease steadily. It drops rapidly during distinct single steep gradients. Between those drops, the algorithm stagnates in the form of long energy plateaus. After a comparably small initial rearrangement effect during the first few hundred steps (not shown in this figure), the progress stagnates. The potential energy reaches its first plateau, indicated by configuration A, highlighted in red. Subsequently, the system drops into three other plateaus, indicated by configurations B, C, and D, until the basin is finally reached. During every drop, the sum of the squared forces increases significantly, while during the plateaus, the forces almost vanish. This evolution suggests that the cascade event can be divided into particular sub-events. In the present example, this results in sub-events A\(\longrightarrow\)B, B\(\longrightarrow\)C, and C\(\longrightarrow\)D, as shown in the three non-affine rearrangement plots at the bottom of Fig. 35. We emphasize that the evolution of the potential energy is constantly decreasing, which is the mathematical necessity of the steepest descent algorithm. Thus, even during the long plateaus, the slope is always negative, reminding us that these sub-events do not separate basins of the PEL. Thus, Maloney and Lemaitre describe the situation as if the system is almost in equilibrium when resting on such a plateau and refer to them as quasi-basins [80].
As it is accepted in the literature that elementary shear banding results from local structural softening in the material [164, 181, 182], one may argue that cascading, which is the interaction of two or more sub-events, may constitute the initial mechanism of elementary shear band formations. Figure 36 depicts a possible simplified mechanism.
Cascade mechanism according to Maloney and Lemaître [82]; a virgin material; b development of a rearrangement that follows the shape of quadrupolar pattern; c development of a cascade by triggering other shear transformation zones
Figure 36a depicts the virgin cartoon material subjected to elementary shear by the pull-push combination into the principal directions \(\varvec{\lambda }_1\) and \(\varvec{\lambda }_2\). Drawing on the concept of shear transformation zones by Falk and Langer [79], there are particular spots, indicated by the gray-filled circles in this figure, prone to rearrangement. Figure 36b indicates the occurrence of an elementary rearrangement, whereby the flow mechanism of the material coincides with the push-pull combination of elementary shear, following the displacement field of a quadrupolar shape with a \(\cos {2\Theta }\)-symmetry pattern [183]. We color-coded the displacement field qualitatively and indicated the zones of maximum shear stress by the two crossing red lines, i.e., where the flow changes its sign. The rearrangement event activates another shear transformation zone, which triggers and interacts with the first rearrangement center, as shown in Fig. 36c. One might interpret this extremely simplified visualization as the mechanism of an elementary shear band that develops in the plane of maximum shear between the rearrangement centers during the cascade. Dasgupta et al. [184, 185] confirm and extend this theory. When the stress increases beyond a particular yield strain, cascades occur in quadrupolar-like rearrangements, aligned at an angle of \(45^{\circ }\) with respect to the principal stress axis. In line with this theory, Lemaître and Caroli [186] pointed out that randomly located irreversible rearrangements are coupled via their surrounding quadrupolar displacement fields. Hinkle et al. [187] applied a coarse-graining methodology to an atomistic model of bulk metallic glass, based on which they were able to develop a well-informed, predictive continuum description of the plasticity in shear bands using the shear transformation zone theory. They developed an optimal coarse-graining length scale, so that shear banding is captured best compared to the molecular dynamics model.
While Hieronymus-Schmidt et al. [188] experimentally confirm that shear bands evolve around a \(45^{\circ }\) angle with respect to the principal stress axes, their study does not align with the shear banding mechanism in Fig. 36, and, consequently, with the orientation of the quadrupoles earlier proposed by Maloney and Lemaître [82]. Figure 37 presents their elementary shear band formation mechanism.
Shear band formation as an alignment of rearrangement events that follow the deformation pictures of Eshelby quadrupoles according to Hieronymus-Schmidt et al. [188]
One sees that the push-pull directions of the quadrupoles are not aligned with the push-pull directions of elementary shear; rather they are alternately rotated by \(45^{\circ }\). However, their shear band model was strengthened by identifying periodic density fluctuations in experimental observations.
We emphasize that shear bands do not always occur at an inclination angle of \(45^{\circ }\). Depending on the loading conditions, the angle can vary between \(40^{\circ }\) and \(50^{\circ }\) [189]. However, theoretically, angles between \(30^{\circ }\) and \(60^{\circ }\) with respect to the principal shear directions are possible [190]. Intriguingly, the angles observed from shear band formation in a macroscopic framework differ from the angle in microscopic correlations. In this context, Karimi and Barrat [183] demonstrated on a mesoscopic model that the varying angle of shear bands can be interpreted as local plastic events in the form of Eshelby transformations combined with a local Mohr-Coulomb yielding criterion that incorporates a pressure-sensitive effect during the propagation of plasticity.
9.2.6 Shear-Induced Elementary Events in Network Glasses
The AQS protocol is identical for Lennard Jones and other glass systems. However, the mechanical response of network glasses includes particular effects that are worth reviewing and might be the origin of the complex fracture behavior of silica glass-based materials. As elaborated in detail in Sect. 8.1, network glasses are an assemblage of tetrahedral units in three dimensions or as triangle units in two dimensions that build complex rings of various shapes and sizes so that this distinct network of directional covalent bonds affects the microscopic picture of rearrangement events. We underline this statement with an illustrative numerical example. We subject the two-dimensional network glass sample from Fig. 25d to the AQS simple shear protocol. After observing the shear response and the first inelastic rearrangement event, one localizes the zone that will first experience a shear-induced transformation, the first shear transformation zone, during the AQS protocol. We present the local atomic configuration of this zone in the stress-free sample in a circular cutout in Fig. 38a.
Circular cutout of a two-dimensional network glass sample around an STZ; a STZ of the stress-free configuration; b STZ shortly before the event; c STZ shortly after the event
The microscopic pictures before and after the rearrangement are depicted in Fig. 38b and c, respectively. Compared to the densely packed Lennard-Jones glass, the network material follows a substantial structural elastic orientation with respect to the principal axes of the simple shear deformation. Roughly speaking, the rings are stretched into the first principal axis while squeezed into the direction of the second principal axis. The rearrangement occurs in breaking a covalent Si–O bond, accompanied by a certain amount of stress release. We emphasize that this is just one typical example out of the infinitely rich repertoire of elementary shear transformation pictures in network glasses. Bonfanti et al. [153] classified elementary events in bulk silica into two types. The first event type is a sudden change in the shape of a local group of rings, while all covalent bonds stay unaltered, i.e., shape-changing events. The second event type is a local breakage of one or more bonds and, potentially, the forming of new bonds. However, both event types have a quadrupolar non-affine system response comparable to the one in Lennard-Jones glasses, shown in Fig. 32b. Consequently, two types of events are observed in three-dimensional bulk silica models [153]. However, to the best of the authors’ knowledge, only bond-breaking and bond-forming events were observed in the literature regarding the shear deformation response of two-dimensional network models [114, 191]. An open question to date is whether these elementary events in silica glass can be again divided into sub-events, as shown in Fig. 35, which could decode elementary inelastic unit events, such as the breakage or the formation of one bond, in silica glass.
Motivated by the picture of this quadrupolar-type of event, Albaret et al. [192] demonstrated on an amorphous silicon model that it is possible to fit the atomic displacement fields on the continuous Eshelby solution, which can, in turn, be used to quantitatively reproduce stress-strain curves—a strategy that may be useful for future mesoscopic models of disordered materials. Furthermore, Boioli et al. [193] characterized shear transformations by comparing them with the continuous Eshelby deformation picture. Using a Stillinger-Weber-type potential, they also pointed out that bond character and direction greatly influence the shear transformation picture.
9.3 Prediction of Inelastic Shear Rearrangement Events
The idea of shear transformation zones preexisting in the material provides the basis for mean-field approaches, such as the shear transformation zone theory, that may describe the essential features of a disordered material subjected to remote deformation [14]. However, the concept of the shear transformation zone theory builds on the characterization of local structural identification parameters, such as free volume, elastic moduli, or local atomic structure. While, for crystals, the identification of zones prone to rearrangements is trivial, as they are localized at or around the positions of dislocations, regarding disordered systems, the local identification based on microscopic structural quantities turns out to be challenging [194, 195].
Whereas thermodynamic molecular simulations are statistical in nature and, therefore, require a set of response statistics, the dynamics of systems in the athermal quasistatic limit are, on the contrary, purely deterministic. If the same simulation is repeated, one will obtain precisely the same result. Generally, this can be seen as an encouraging fact since local zones of rearrangements and their triggering strain must then somehow be a priori decoded in the dynamics. As mentioned before, not only the microscopic configurations themselves may contain the rich dynamics of amorphous solids but also the external loading history to which they may be subjected [175]. The issue of prediction is caused by the highly nonlinear, even chaotic nature of such problems and, therefore, by the lack of closed-form solutions available in the literature. It may even be questioned if such a solution could ever exist. It is also worth mentioning that the PEL can significantly change its shape due to deformation, while only navigation through its nearest vicinity is, to date, possible using numerically expensive methodologies. Therefore, prediction strategies build instead on harvesting potential candidates, based mainly on configurational information, and are limited to events that may happen in the nearer future.
However, some promising structural correlations have been developed in the last decades. In this context, the main connections between local structural indicators and the appearance of local zones susceptible to rearrangements were elaborately summarized by Richard et al. [196]. They categorized the overall group of structural indicators into five subgroups, which they classify as (i) purely structural indicators, (ii) machine learning-enhanced structural indicators, (iii) indicators investigating linear modes, (iv) linear response indicators, and (v) indicators involving the local nonlinear structural response. We will discuss some main prediction strategies—not necessarily in the above-classified sequence—below. However, for a direct comparison of structural predictors, including numerical studies, we refer the reader to the summary of Richard et al. [196].
The oldest identification quantity that comes to mind is the local free volume, as first suggested by Spaepen [163]. However, as mentioned before, this is a relatively poor prediction indicator. Nevertheless, the local Voronoi cell anisotropy, which enriches the free volume by directional information, showed high correlations with zones susceptible to rearrangements [197]. Therefore, it may be worthwhile to investigate further predictors that originate exclusively from atomic positions since they are computationally cheap.
Ding et al. [198] introduced the local “flexibility volume”, the local squared vibrational displacement extracted from multiple NVE runs around an equilibrium position multiplied by the local atomic volume. As stated by the authors, this prediction quantity does not only represent the local space free for dilatation and relaxation, i.e., the atomic free volume, it also adds information about the local configurational dynamic response. Introducing the “local flexibility volume,” they present prediction maps that correlate well with shear transformation zones in metallic glasses.
The first strategy that comes to mind to predict shear transformation zones, which does not depend on structural indicators, is to monitor the local shape of the potential energy basin that the system is situated in during deformation. Such an investigation is possible by linearizing the problem by evaluating the Hessian of the system. This method is, therefore, an indicator from a particular configuration \(\varvec{r}\) without investigating the nearer neighborhood of the PEL at a fixed deformation state. As indicated above, it is possible to predict rearrangement events in the nearer vicinity of the present configuration by the identification of the vibration modes that correspond to the lowest eigenvalues, also called soft vibrational modes [195, 200,201,202]. When approaching a shear transformation, the eigenvalue corresponding to the lowest vibrational mode approaches zero in terms of a square-root singularity (146), and the non-affine displacement field of the rearrangement event correlates highly with the soft vibrational mode directly before the event. However, comparing the displacement pictures of the eigenmodes of the Hessian shortly before the rearrangement in Fig. 34 with the modes that correspond to the initial configuration in Fig. 28 reveals the complex nonlinear nature of disordered materials. Except for the Goldstone modes that describe rigid body motion and the mode in Fig. 34e, which describes a global material displacement, modes three to six differ fundamentally. Through a series of stress drops and the corresponding atomic rearrangements, the configuration experiences a sequence of microscopic transformations, which makes a prediction of the whole trajectory from an inspection of only the initial configuration almost impossible.
The lowest eigenmodes of a configuration are a state-dependent quantity, which complicates a prediction when changing the deformation protocol, as critically pointed out by Gendelman et al. [175]. Thus, when predicting a local spot prone to rearrangement, one must acknowledge that shear transformation zones are protocol- and state-dependent. Xu et al. [203] proposed a promising prediction approach that considers these circumstances. In particular, when using the vibrational coordinate \(q^{(1)}\) for a prediction of the rearrangement displacements as \(\varvec{x}^{(1)} = \varvec{\Phi ^{(1)}}q^{(1)}\) the rearrangement direction of the prediction of the shear transformation zone is somehow fixed. If the shear direction changes shortly before the rearrangement, the same inelastic event is not observed, and the corresponding vibration mode does not predict the rearrangement. This phenomenon is also in line with the critical comments by Gendelman et al. [175]. However, one solution to deal with this problem, according to Xu et al. [203], is to introduce the second and third derivative of the potential energy with respect to the modal coordinate \(q^{(1)}\) at the configuration shortly before the rearrangement. Based on this, they evaluate the triggering strain for different shear orientations and the softest shear orientation of the lower mode. As a result, the prediction approach uses atomic nonaffinity based on the anisotropy of the shear stress derivative and the coordinate of the lowest vibration mode \(q^{(1)}\), varying the orientations. For a prediction far away from the instability, they suggest using a whole set of lower vibration modes for the approach. Furthermore, Baggioloi et al. [199] showed in a very recent paper that an investigation of the non-affine displacement fields during the deformation protocol leads to a successful prediction of shear transformation zones and shear banding, quantitatively characterized by the Burgers vector [23].
A different but promising prediction strategy was proposed by Patinet et al. [47] and further extended by Barbot et al. [176]. The strategy is to define a scanning region, as shown in Fig. 39a, which confines an atomic subsample much smaller than the actual glass sample.
The scanning window may either proceed from atomic position to atomic position or move along an equidistant grid covering the entire sample. Every local visit corresponds to a series of AQS shear deformation protocols on the subsample. In doing so, the subsample is rotated at different angles, whereby the rotation angle \(\alpha\) varies between \(0^{\circ }\) and \(180^{\circ }\). The atomic configuration inside the subsample is divided into two groups. During the AQS protocol, both groups are subjected to the affine deformation field. However, only the atoms in group I can move freely during deformation, while the atoms in group II are frozen, as they are excluded from the position update. They extracted the minimum shear stress of the subsample projected into the loading direction of the glass sample. This local quantification predicts shear transformation zones and could even serve as an ordered prediction sequence of when those shear transformations could occur. Ruan et al. [204] have applied the local yield stress method to a three-dimensional glass showing the local anisotropy of the material, which can be applied to support continuum-mechanical constitutive models. This method was also employed to enable further understanding of the origins of rejuvenation [164] and polarization [205] and to pave the way from atomistic simulations to elastoplastic models [206, 207].
Lately, the activation-relaxation technique [208,209,210,211] has been used as a quite successful prediction strategy. The activation-relaxation technique allows one to travel across the nearer vicinity around the inherent structure in the PEL. Motivated by a well-known problem in real life, where one has to conquer a pass between two mountaintops to get from one valley to another, the activation-relaxation technique aims at carrying the system from one minimum of the PEL to the other by passing over a saddle node. A visualization of this procedure applied to an extremely simplified comic PEL is presented in Fig. 39b. Starting from a local minimum, the activation-relaxation technique proceeds in two steps: (i) the activation step, where the system moves towards a saddle point of the PEL (denoted as S in Fig. 39b); and (ii) the relaxation step, where the system is pushed over the saddle node and relaxed into the new minimum. The activation step (i) starts by displacing one or more atoms in a random direction. From this point on, the system is iteratively pushed out of its minimum by the application of an artificial force \(\varvec{g}\) that is step-wisely adapted:
In this approach, \(\varvec{f}\in \mathbb {R}^{3N}\) is the force on the system due to the introduced perturbation vector \(\tilde{\varvec{x}}\in \mathbb {R}^{3N}\) that moves a selected atom group out of the minimum, while \(\alpha ^{(ar)}\) denotes a control parameter. The force \(\varvec{f}_c\) always points away from the minimum. From the point of view of practical implementation, the activation step (i) should be divided into two substeps: the stabilizing phase and the convergence phase. In the first substep, the system must escape the local minimum. Using an entirely random displacement vector, one observes a significant correlation between the random displacement and the lower vibration modes in such a manner that the system tends to fall back into its initial minimum. Thus, one can neglect the contribution of the lower eigenmodes by setting the initial displacement vector as a linear combination of the higher eigenmodes and escaping the minimum step-wisely until a certain threshold is reached. Rodney and Schuh [45], who used the activation-relaxation technique for Monte Carlo simulations of a Lennard Jones glass, set the threshold to a minimum curvature to become smaller than a negative threshold value, for example. Subsequently, the algorithm enters the second part of the activation step, i.e., the convergence phase, following the direction of \(\varvec{g}\) until the system approaches the saddle node. This is the case when \(\varvec{f} \cdot \frac{\tilde{\varvec{x}}}{\left| \tilde{\varvec{x}} \right| }\) changes its sign from a negative to a positive value. Once the system is lifted over a saddle node, the subsequent relaxation step can be realized by applying an appropriate minimization algorithm, such as the conjugate gradient method, as discussed in Sect. 7.2. While the relaxation step is straightforward, the activation step needs careful attention to ensure convergence to a saddle node. The algorithm is somewhat sensitive to the input parameters, and it can often result in unsuccessful searches and fail during the destabilizing phase if \(\alpha ^{(ar)}\) is set at a value that is too small.
Mechanical deformation results in the minimum dissolving into a first-order saddle node and the system dropping into an adjacent minimum. One may “artificially” excite such a rearrangement event using the activation relaxation technique, neglecting any external deformation, thereby investigating the system only in its initial state. Such a prediction was performed by Xu et al. [212] and applied to a bulk metallic glass. Out of a metallic glass sample with \(10^4\) atoms, they harvested 20 candidates for rearrangement events out of \(2\times 10^5\) searches. This result is not surprising. As indicated in Fig. 39b, the PEL is a high-dimensional, highly complex shape consisting of numerous minima and maxima. Therefore, an inspection using the activation-relaxation technique will always result in multiple rearrangement candidates. Furthermore, they included the stress gradient from the corresponding mode of the predicted rearrangement events, which they realized by following the minimum energy path from the initial minimum to the adjacent minima using the nudged elastic band method [213, 214]. This strategy allowed them to include protocol-dependent information in the prediction procedure. However, the activation-relaxation technique could not predict all of the events, which also stimulated them to unload the system after each rearrangement. They observed that rearrangements that recovered during unloading could not be predicted from the initial state. In contrast, rearrangements that recovered only after unloading and further reverse loading into the negative strain range could be predicted. In other words, prediction using the PEL is only possible if the corresponding strain state is part of a hysteresis confined by a triggering and a recovering strain.
Lately, the complexity of the prediction problem has encouraged researchers to add data-driven algorithms into the toolbox of prediction strategies. The objective is to feed machine-learning models, such as neural networks, by response quantities, such as positions or non-affine displacement fields, and train the network models to return prediction quantities, such as the rearrangement position and strain. Cubuk et al. [215] spotted candidates susceptible to undergoing atomic-scale rearrangements on two disordered systems, i.e., frictional granular packings and two-dimensional Lennard-Jones glasses. In particular, they classified the particles as “soft”, i.e., susceptible to rearrangements; or “hard”, i.e., not susceptible to rearrangements, using support vector machines [216]. Wang et al. [217] applied a tree-boosting algorithm to find atoms with strong resistance or high compliance to thermal \(\beta\) activation, which they numerically realized by employing the activation relaxation technique [211]. They fed the machine learning algorithm by a local interstice measure [218], defining the susceptibility of the material subjected to shear loading as the output value. Fan et al. [219] proposed local structural representation to predict stress-induced shear transformations taking into account the loading direction. In doing so, they introduced a spatial density map to characterize the local atomic structure, which results in two- or three-dimensional images, depending on the system’s dimensionality. They used several neural network architectures, i.e., convolutional neural networks, graph neural networks, and support vector machines, to process the images and return the prediction output. The apparent advantage of such methods is the possibility of capturing unexpected combinations of features, which seems impossible when using conventional strategies. At the same time, one must always admit that machine learning methods are mostly black boxes fitting massive parameter sets on training data sets. After comparing the predictions with the training data set, one can show that machine learning methods work for data sets they did not see before to a certain extent and physically interpret results. Thus, the reliability of such algorithms depends significantly on a carefully chosen, representative training data set. However, a general problem is the apparent lack of physical interpretation of the fitting parameters and, therefore, a missing physically motivated description. Furthermore, it must be emphasized that rare outliers not contained in the fitting parameter scale may correspond to important cascade events, so machine learning-enhanced predictions may fail entirely.
9.4 Recovering of Rearrangement Events
The arrows in the illustrative PEL in Fig. 39b reveal several thermally activated events using the activation-relaxation technique. One may identify an event as a journey from the inherent structure in minimum \(\textcircled {1}\) to the inherent structure in minimum \(\textcircled {2}\). It is clear that once the system rests in \(\textcircled {2}\), one may also jump back into the original position \(\textcircled {1}\). The event back into the original position corresponds exactly to the reverse deformation of the initial event in the configuration space. As a result, recovering thermally activated events is trivial from the perspective of the PEL. It may even be possible for the configuration to jump from the initial state \(\textcircled {1}\) to \(\textcircled {2}\), from \(\textcircled {2}\) to \(\textcircled {3}\), and from \(\textcircled {3}\) back again into the initial state \(\textcircled {1}\). Such a detour through the PEL could be arbitrarily extended.
However, mechanical deformation continuously alters the shape of the PEL. The distances between the maxima and minima may undergo significant changes, and minima can transform into first-order saddle nodes and vanish. If the configuration lies in such a minimum, it experiences a rearrangement event, as extensively discussed in Sects. 9.2.4 and 9.2.5. When loading a configuration until a certain stress drop, followed by subsequent unloading, it may be possible that the system recovers upon unloading. The rearrangement at the recovering strain is then the reverse event at the recovering strain to the rearrangement at the triggering strain. Such a situation is depicted in Fig. 40a.
Two examples of loading until the first rearrangement, followed by subsequent unloading and reverse loading; a a rearrangement event that recovers upon unloading; b a rearrangement event that does not recover upon unloading and reverse loading—a recovering strain due to further reverse loading may not happen at all
As also pointed out by Lundberg et al. [220] and mentioned above, an event may not recover upon unloading and reverse loading, as indicated in Fig. 40b. In this case, a thorough investigation of the PEL using the activation-relaxation technique starting from \(\textcircled {1}\) will find a nearby minimum corresponding to the inherent structure \(\textcircled {2}\). Using the method of mechanical deformation, one may load, unload, and reverse load the system until an arbitrary rearrangement occurs that is not identified as the reverse event to the triggering rearrangement. Such a rearrangement can change the shape of the PEL significantly and, therefore, also the minimum \(\textcircled {1}\) and its neighborhood so that the system may not likely fall into this minimum again. As also pointed out by Xu et al. [212], the situation in Fig. 40a is different. The minimum \(\textcircled {2'}\) cannot be found via the activation-relaxation technique starting from \(\textcircled {1}\), as it does not exist yet. One must load the system until a particular state \(\textcircled {1'}\) from which the minimum \(\textcircled {2'}\) can be accessed via activation-relaxation. The strain range in which \(\textcircled {1'}\) must be chosen is a priori unknown, and continuous activation-relaxation testing, while increasing the strain or mechanical cyclic loading, is necessary to access \(\textcircled {2'}\). Regarding the denomination of these two types of inelastic events, we draw on the definition of Xu et al. [212], who denote rearrangements that recover upon unloading as anelastic events and rearrangements that do not recover upon unloading as plastic events.
The highly nonlinear nature of molecular systems combined with the complex dynamics of amorphous plasticity makes predictions of the respective types of inelastic events almost impossible from just navigating through the PEL at certain strain states. Thus, Bamer et al. [221] proposed a data-driven approach using a polynomial regression technique and presented their results on bulk silica glass. This method allows for predicting or classifying the types of events while harvesting physical information in terms of polynomial coefficients. However, such data-driven methodologies are strongly linked to the material’s properties, and generalizations to other disordered materials may be challenging.
9.5 Other Deformation Protocols
Although local shear constitutes the most dominant deformation protocol inducing material damage, deformation states may be defined by subjecting the cell geometry to any arbitrary deformation gradient in such a manner that elementary events are not limited to happening during shear transformation only. In particular, this also concerns materials that experience local deformation modes emerging from elementary processes that also include volume changes, such as silica glass. In this context, fracture simulations often provide beneficial information to enable a better understanding of the material response behavior better. Thus, this section will discuss essential features from simulations following other deformation protocols than shear loading.
9.5.1 Tensile Deformation
Generally, tensile deformation protocols start from rectangular cell shapes, although they are not restricted to this requirement. One may distinguish between two different types of tensile deformation protocols: (i) uniaxial tensile deformation, which is the alteration of one diagonal element of the deformation gradient \(\varvec{F}\) that generally results in a multi-dimensional corresponding stress response, and (ii) simple tension, which is a combination of an increase in the strain in the axial direction followed by a strain relaxation in the perpendicular direction. Regarding the second variant, the deformation gradient is not predefined since the level of relaxation in the perpendicular direction depends on the Poisson’s ratio of the material and generally alters throughout the protocol. Thus, the component of the stress tensor in the perpendicular direction of deformation vanishes in this case. Regarding either of the two methods, the mathematical mechanism is equivalent. The PEL experiences a deformation-induced growing shift in its shape until the system reaches a saddle node and drops into an adjacent minimum. From a mechanical point of view, these deformation types are fundamentally different from shear since they may favor the nucleation of voids caused by a volume change in the simulation box. Eventually, the material in the simulation box, consisting of a set of linked atomic particles, can end up in a collection of two pieces, each of which consists of a set of linked atomic particles. In other words, the material is broken and cannot take any more tensile stress.
Tensile simulations of Lennard Jones glasses are rarely found in the literature. They were performed to investigate the fracture process of disordered solids and the elementary growth of nano-cracks in the material [222, 223]. The lack of such simulations is mainly caused by the fact that metallic glasses fail in the form of cascades of shear events. Uniaxial tensile simulations are mechanically more complex as they produce a combination of tensile, transverse, and shear stress in the material. We prepared a simulation by performing a uniaxial tensile test on the same sample subjected to simple shear, as shown in Fig. 31. In particular, we present the uniaxial tensile test in Fig. 41 and the simple tensile test in Fig. 42.
Uniaxial tension of a a two-dimensional Lennard Jones glass
Simple tension of a two-dimensional Lennard Jones glass
Uniaxial tension of a two-dimensional network glass
Simple tension of a two-dimensional network glass
The tensile test in Fig. 41 presents an effect that may macroscopically be defined as mode I fracture [224]. When performing molecular simulations with many atoms, local dislocations in the form of a notch are often introduced. This approach allows one to locally induce the nucleation of the fracture process [225]. One may interpret them as downsized laboratory samples, which allow one to study the crack propagation phenomenon in detail. However, at the same time, this procedure has the disadvantage of artificially introducing a geometrical disturbance in the form of a notch that favors localized potential spots of rearrangements around the notch and therefore hinders a detailed investigation of the crack nucleation in general. In the simulations in this paper, we do not apply such crack initiators and let the crack form itself due to the amorphous nature of the structure. The “softest” spots experience isolated rearrangements, whereas the bi-axially increasing stress-strain state favors the formation of nano-voids instead of shear transformations. The interaction of these nano-voids leads to accumulating a distinct crack with increasing tensile strain. The non-affine displacement fields on the right side correspond to the stress drops highlighted in red in the stress-strain curves in the middle column of Fig. 41.
The simple tension response of the densely packed glass, presented in Fig. 42, shows a response behavior that reveals, in this case, a significantly larger plastic capacity before the material experiences considerable damage effects. The relaxation into the transverse direction of the simple tension protocol allows for the occurrence of atomic scale rearrangements that are more similar to shear transformations than to the rearrangements of the uniaxial response since the void formations due to axial movement are more likely to be closed through the relaxation phenomenon into the transverse direction. However, we emphasize that a quantitative investigation of the mechanical properties of densely packed materials, such as bulk metallic glasses, depends on the rich variety of the possible atomic composition, and we refer to the relevant literature [12].
However, the investigation of network glasses reveals more intriguing mechanical phenomena when investigating the mode I fracture process. Thus, we subject the two-dimensional network material from Fig. 25c to both uniaxial and simple tensile deformation protocols. The results are presented in Figs. 43 and 44.
From the uniaxial tensile simulation, one observes that the fracture process is a successive occurrence of bond-breaking events in the principal tensile stress direction. In both uniaxial and simple tensile loading cases, voids emerge through bond breaking, and smaller rings break and dissolve into larger rings—an effect that has also been observed in bulk silica glass samples [13]. With increasing strain, these larger rings eventually merge into a distinct mode I nano-crack and the formation of two pieces of material in the simulation box. At this stage, the material is broken and cannot take any further tensile stress. Another interesting phenomenon observed in the results of these two tensile examples is the effect of softening with increasing tensile strain. With every bond that breaks, stress drops occur within elementary rearrangement events. We also note a decrease in material stiffness when inspecting the elastic stress-strain branches with increasing tensile strain. This phenomenon is caused by the fact that the number of active bonds in the first principal direction decreases with the progressing fracture process.
The literature investigating the fracture phenomenon of network glasses under tensile simulations is quite extensive. The fracture process of silicate glasses depends on the chemical composition and additives, such as sodium or aluminum, that alter the network structure and greatly influence the mechanical behavior. For further information in this regard, we refer to the relevant literature [226,227,228,229]. As already mentioned in Sect. 1, we concentrate on binary glasses that exhibit networks of SiO\(_3\) triangles or SiO\(_4\) tetrahedra in two and three dimensions, respectively. Intriguingly, the response behavior in Figs. 43 and 44 reveal a different fracture behavior than one would expect from the everyday use of silicate glasses on the macro scale, which is brittle. However, while surface cracks mainly govern the macroscopic brittle fracture behavior [230, 231], it was experimentally proven that silica glass exhibits ductile fracture behavior on the nanometer scale [232]. In particular, nanoindentation results suggest plastic deformation in the mechanical response of silicate glasses [233], and it is expected that the silica glass significantly plastifies in the vicinity of the crack tip during fracture [232]. However, the truth of the latter statement has resulted in controversy [234, 235].
From a purely mechanical deformation perspective, one may interpret – in a simplified framework—the dynamics of the tensile crackling of a network glass as the loading procedure of two rigid plates linked through a set of truss bars arranged in parallel. An engineering analogy is depicted in Fig. 45a. Such descriptions are well-known as fiber bundle models in the statistical mechanics community, and although these models are simplified, they can display the non-trivial effects of the fracture phenomenon [236]. Every truss bar behaves in a purely elastic manner until its force threshold is reached. At this moment, it fails and cannot carry any increase in load. This is the activation of a tensile or simple tensile soft spot, which one may interpret as an elementary tensile event, equivalent to a shear transformation [14], as shown in Fig. 31d, e. Another mechanical interpretation would be the concept of statical indeterminacy. In the case of every bar having the same stiffness and tensile strength, one would observe a sudden fracture process, after which the material cannot take any tensile load. As a result of the variety in the truss strengths, fracture occurs gradually during the loading process, and total material failure may be indicated long before it finally happens. In addition, the softening effect is considered in such a simplified model as the evolution of the material stiffness depends on the decreasing number of intact truss bars.
Influence of the network heterogeneity of a Zachariasen network glass on the ductility of the material; a truss analogy of fracture; b and c mean stress-strain result of 100 uniaxial and simple tension simulations, adapted from Ebrahem et al. [238]
It was shown that ductility depends greatly on the level of heterogeneity of the silica network structure [89, 100, 237, 238]. We refer to the discussions in Sect. 8, where the generation of silica glass structures is discussed, taking into account, the level of heterogeneity. The higher the level of heterogeneity, the lower the tensile strength but, the higher the ductility of the material. This effect was demonstrated on both tensile simulations at ambient temperature using the canonical ensemble [89] and athermal quasistatic tensile simulation protocols [238]. Again, this effect can be interpreted with the help of the simplified mechanical model in Fig. 45a. On the one hand, a high level of heterogeneity refers to larger deviations in the force thresholds of the truss bars. One observes the rupture of the first truss bars long before all trusses finally break. On the other hand, a very low heterogeneity refers to smaller deviations in the force thresholds of the truss bars, and the others break shortly after the first truss breaks. If all trusses have the same strength, one observes a purely brittle fracture phenomenon, which is observed when crystalline network structures are subjected to tensile loading [239, 240]. This effect has also been observed by Pedone et al. [241], who performed uniaxial and simple tensile tests on bulk silica using the Berendsen thermostat at ambient temperature on both crystalline and vitreous phases of silica glass. Chowdhury et al. [242] investigated the mechanical behavior of silica glasses using reactive force fields. They also reported that varying the quenching rate alters the structure of silica and, consequently, their mechanical response behavior when subjected to tensile loading.
Figure 45b and c show the mean stress-strain responses of uniaxial and simple tension simulations of a set of one hundred model network glass samples generated using the dual Monte Carlo bond switching algorithm. These results are adapted from Ebrahem et al. [238]. For this series of molecular simulations, the standard deviation \(\sigma {(m)}\) of the logarithmic distribution was fitted to the statistics extracted from pictures of real two-dimensional silica and defined as the reference heterogeneity \(s^2_m\), from which five levels of the target heterogeneity were evaluated: \(0.6\,s^2_m\), \(0.8\,s^2_m\), \(1.0\,s^2_m\), \(1.2\,s^2_m\), \(1.4\,s^2_m\). For both deformation protocols, one observes a clear brittle-to-ductile transition when increasing the level of heterogeneity of the material. This effect has also been observed for bulk silica [89, 100] and amorphous graphene monolayers [137].
9.5.2 Pressure and Shear Deformation Under Pressure
Densification of the two-dimensional model network glass sample
These deformation protocols are highly relevant for network materials that densify when subjected to external pressure. Furthermore, network glasses may not only densify when subjected to external pressure, but they also show anomalous shear behavior that varies with the external pressure [243,244,245].
Firstly, we qualitatively demonstrate the densification process by subjecting our two-dimensional silica sample from Fig. 25c to external volumetric deformation using the AQS protocol realized by step wisely altering the two cell vectors of the two-dimensional simulation box. The result is depicted in Fig. 46. We emphasize that, for convenience, we plot the absolute value of the volumetric strain \(|\varepsilon _v|= |\frac{\varepsilon _{11} + \varepsilon _{22}}{2}|\), which is why the curve increases to the right, although the volumetric strain is negative. The volumetric stress-strain response results in a highly nonlinear curve that has a considerable increase in its slope with an increase in the absolute value of the volumetric deformation. During the densification process, the material also experiences sudden stress drops, enlarged for the first two drops in Fig. 46, where the network structure collapses locally, and the material experiences distinct rearrangements that may not follow the typical quadrupolar shapes. After higher densification strains, over-coordinated local configurations occur, as highlighted by the red circles in Fig. 46. At very high densification strains, the atomic configuration reassembles itself into a more ordered structure, as there is less space for local density fluctuations caused by the disorder of the network material. The observations from this two-dimensional model can be equivalently transformed into a three-dimensional structure. Jin et al. [246] observed from molecular dynamics simulations that the Si–O coordination of bulk silica glass changes from four to six when densified and the bond angle \(\vartheta _{\text{O}{-}\text{Si}{-}\text{O}}\) from 109\(^\circ\) to 90\(^\circ\). The structure of bulk silica glass changes from corner-sharing SiO\(_4\) tetrahedra to corner and edge-sharing SiO\(_6\) octahedra. In the two-dimensional model framework, the network changes from SiO\(_3\) triangles to SiO\(_4\) quadrangles, as shown in Fig. 46. However, Tian et al. [247] showed that porous silica does not experience over-coordination when subjected to volumetric loading, but the medium-range order alters while the saturation increases. Furthermore, Yuan et al. [248] showed that densified silica glass reveals higher ductility after unloading when subsequently subjected to tensile loading. This effect can be attributed to the energy dissipated when over-coordinated Si–O bonds release during tensile loading, providing additional plastic capacity.
Results of a two-dimensional network glass subjected to shear deformation under pressure (adapted from Bamer et al. [114]); top: simple shear response and simple shear response under pressure; middle and bottom: mean simple shear response under external pressure evaluating the results of one hundred samples
We finish this section by reviewing the anomalous yielding behavior of network glasses. Intriguingly, the shear response behavior alters significantly with increasing pressure and reveals surprising fluctuations regarding shear modulus, elastic limits, material strengths, and yielding plateau [249, 250]. We qualitatively discuss these properties by adapting the results from Bamer et al. [114]. In doing so, one hundred two-dimensional model network samples were generated using the dual Monte Carlo bond switching algorithm, as discussed in Sect. 8. Every sample was densified to a maximum value of 15% in 2.5% steps and then subjected to simple shear deformation using the AQS protocol. The top plot in Fig. 47 shows the simple shear response of one sample when subjected to zero pressure and the simple shear response of this sample subjected to 15% densification. One immediately realizes that the material behaves in a fundamentally different manner in both the elastic and plastic range. To further investigate this effect, we present the middle and bottom plots, where the mean response curves of the one hundred samples subjected to different levels of densification followed by simple shear deformation are presented. One observes a decrease in shear modulus until around 5% densification. Further increasing the pressure leads to an increase in the shear modulus. However, the elastic range decreases when the pressure increases further from 5%. At 15% densification, the material plastifies very early and gradually approaches a yielding plateau. These findings are qualitatively in line with Schill et al. [250], who investigated the anomalous yielding behavior of bulk silica by performing molecular dynamics simulations using the isobaric isothermal ensemble. Furthermore, Bamer et al. [114] investigated the deformation picture of elementary events of the two-dimensional network glass subjected to shear and pressure. Intriguingly, a shear transformation at zero pressure, as indicated in Fig. 38, results, on average, in a stronger outward flow and a rather weak perpendicular inward flow since it often occurs during a bond-breaking event in the first principal direction and the formation of a nano void. However, the external pressure forces the material to flow into the free volume in the second principal direction of simple shear, as indicated in Fig. 9. They applied the divergence theorem to the atomic flow field to quantify the relation of the in- to outward atomic flow. They defined the center of the rearrangement by spotting the atom with the highest non-affine displacement. The position of this atom serves as the center of a radial region \(\Omega\) to quantify the in-to-outward atomic flow by a magnitude \(f(\Omega )\). Consequently, a quantification of the atomic flow is provided by:
where \(\varvec{n}\) denotes the normal unit vector in the outward direction. A negative value for \(f(\Omega )\) indicates that more material enters the region \(\Omega\) in the second principal direction of shear compared to the amount of material leaving \(\Omega\) in the first principal direction of shear. Vice versa, a positive value for \(f(\Omega )\) indicates that more material leaves the region \(\Omega\) in the first direction of principal shear compared to the amount of material entering \(\Omega\) in the second principal direction of shear, i.e., the formation of nano-voids. If \(f(\Omega )\) is equal to zero, no change in the local material density occurs during the rearrangement event. Intriguingly, Bamer et al. [114] showed that, with increasing pressure, the average in-to-outward atomic flow decreases from an initial positive value, crosses the zero mark, and further decreases until a negative minimum. With an even further increase in pressure, the average in-to-outward atomic flow increases again until \(f(\Omega )\) vanishes again, i.e., no average local change in density is observed.
At this point, we note that the mechanics of silicate glasses, particularly their brittle fracture behavior on larger scales and how to advance their mechanical behavior, are the subject of extensive studies with an experimental focus. In this context, we refer to the relevant literature, cf. [251,252,253,254,255,256,257,258].
10 Statistical Classification
In the athermal quasistatic limit, the deformation mechanics of disordered systems, such as the Lennard Jones glasses subjected to shear deformation history, visualized in Fig. 48a, are governed by elastic response branches intersected by sudden stress drops, as discussed elaborately in Sect. 9. This results in a category of materials that respond to a slow-driving mechanical load in terms of intermittent dynamics, characterized by sudden stress drops of different sizes.
Stress drop statistics of a Lennard Jones glass subjected to simple shear; a stress-strain curve and extracted stress drops; b stress drop statistics; note that the stress drop statistics follow a power law distribution
As also discussed in Sect. 9.3, the prediction of the occurrence of a stress drop turns out to be challenging and is the focus of current research. The stress drop magnitudes are proportional to the energy released. Therefore, how far the system drops down in a highly complex PEL that continuously changes during deformation is highly relevant. At first sight, one may conclude that the stress-strain response of the system is in a state of deterministic chaos. However, instead of pursuing prediction strategies, as discussed in Sect. 9.3, one may follow a statistical approach to classify and understand the underlying dynamics in terms of statistical rules. Thus, we subjected 100 glass samples to simple shear deformation. We collected all stress drops and evaluated the probability of occurrence of their absolute magnitudes. This probability is presented as a histogram on a logarithmic scale and marked with quadratic marks in Fig. 48b. Small stress drops occur very often, while large ones occur rarely. We fitted the histogram to an exponential distribution function extended by a cutoff term. The function is written as:
where \(\kappa _1\) and \(\kappa _2\) are constants to be evaluated during the fitting procedure. The cutoff term becomes less pronounced with increasing sample size. Thus, it can be shown that the stress drop statistics are, overall, proportional to a power-law distribution \(P\left\{ \Delta \tau \right\} \propto \Delta \tau ^{-\kappa }\). Thus, one may conclude that the deformation mechanics of disordered systems are not entirely irregular but follow a particular pattern. The fact that the probability density function extracted is proportional to an exponential distribution indicates that the system may be categorized as a part of the framework of self-organized criticality [259,260,261], which describes a dynamical phenomenon initially proposed by Bak et al. [262]. Self-organized criticality has later been shown to apply to a broad range of other phenomena observed in nature, for example, the appearance of forest fires and earthquakes, neuroscience, ecology, etc. [263, 264]. Following this theory, the disordered structure is continuously in a critical state of equilibrium. It organizes itself in terms of catastrophic rearrangement events responsible for the jerky flow in Fig. 48a. The original example of self-organized criticality is the sandpile model by Bak et al. [262], and it reveals a scaling parameter of \(\kappa = \nicefrac {3}{2}\). However, this scaling parameter varies when fitting the stress drop statistics of molecular models [265, 266]. In this context, Liu et al. [267] have proposed that the rate dependence influences the statistics of stress drops, and they provide an extensive study regarding the exponents of the fitting distributions.
The phenomenon of self-organized criticality has also experimentally been proven by Sun et al. [268] by characterizing plastic yielding in metallic glasses as self-organized critical. Furthermore, it has been shown that the crack growth phenomenon in disordered solids appears as self-organized critical and, to some extent, predictable scaling laws [269].
We conclude that the mechanics of disordered systems follow power-law statistics and are, therefore, in a critical state during deformation. At this key point, one may argue that such systems are fractal [270], which implies that the statistics of the mechanics of disordered systems are a scale-invariant phenomenon. In this context, we note that such scale-free distributions have also been observed in elastoplastic models, describing the mechanical behavior of disordered solids under shear, cf. [271]. It is noteworthy that the critical exponent changes with the deformation state [265].
11 Conclusions and Recommendations
Although extensive research has been done for decades, the complex and rich mechanical behavior of disordered materials continues to be poorly understood. This paper reviews and discusses the mechanics of disordered solids using molecular mechanical descriptions. Although the molecular descriptions are simplified by investigating the mechanical response of over-damped structures in the athermal quasistatic limit, the complexity that originates from the structural disorder may challenge researchers for decades to come.
The elastic response of disordered systems does not obey the Cauchy-Born rule, but it results in a diffusive displacement field reminiscent of turbulent fluid flow fields. The inelastic response is governed by sudden, localized rearrangement events that are the size of only a few particles, followed by an elastic response that reveals the strain field of a quadrupolar shape. These very localized phenomena support the idea that disordered solids have defective zones prone to experience atomic-scale rearrangements. In the case of densely packed systems, the local structure transforms primarily due to shear deformation, leading to the definition of a shear transformation zone [79, 163, 165,166,167,168, 170]. In the case of network materials, this localized process may be more complex and might also be induced by volumetric loading or combinations of both volumetric and shear loading.
The existing research gaps in this area are numerous. Obviously, filling those gaps is an ambitious objective and would lead to a better understanding of the mechanics of disordered solids. Thus, our thoughts presented in the following may be seen as recommendations with a personal note.
11.1 Predicting Localized Events
The pre-existence of zones susceptible to atomic scale rearrangements has been the subject of an ongoing discussion during the last decade. Some researchers have disagreed regarding the pre-existence of such rearrangement zones since they argued that such zones might change with the deformation protocol applied, such as the direction of the shear deformation imposed [175]. However, this directional aspect has been elaborately discussed by Barbot et al. [176], who showed that shear transformation zones have slip orientations. Thus, according to them, the change of loading direction “turns on” and “turns off” zones susceptible to experience rearrangements. Consequently, we argue that taking into account external deformation protocols does not prohibit the prediction of rearrangement zones but makes it much more complex and challenging. We conclude that, in an athermal quasistatic framework, such rearrangement regions pre-exist in the material since, from a mathematical perspective, the mechanical problem remains deterministic, although the nature of the response turns out to be highly complex and even chaotic.
Since the existence of shear transformation zones seems to have been accepted by the communities investigating disordered solids, the last decade has been occupied by predicting their location and triggering strain, i.e., where and when they may occur. A forecast of the mechanical response of a disordered solid on the nanoscale based on its properties taken from the undeformed state seems to be a Holy Grail of the mechanics of disordered materials. Predicting strategies utilize input quantities that range from structural indicators, such as free volume [163]; packing [272]; short-range order [182]; indicators originating from the linearization of the system, such as elastic moduli [194] or low vibrational soft modes [200,201,202, 273,274,275,276,277]; and more refined indicators, such as flexibility volumes [198]; nonlinear plastic modes [155] and machine learning approaches [215, 278]. However, it turns out that, to date, the local yield stress predicts shear transformation zones with the highest goodness while also considering the directional information [47, 176] and therefore constitutes one of the most promising approaches. However, one downside of this approach is that the local identifications are realized through the generation of a whole set of molecular simulations for every scanning grid point, so the identification of the yield stress field may become quite cumbersome for larger scanning grids. Thus, the application of surrogate models, such as model order reduction or machine learning-based approaches that include the local atomic configuration, may prove to be of great added value to the numerical realization of the local yield stress method.
The prediction of shear transformations mainly focuses on densely packed Lennard-Jones systems and metallic glass models in the literature. However, an adaptation of this idea to network glasses, such as silica glass, is, to date, missing from the literature. In this paper, we have shown that a network glass behaves fundamentally differently than a densely packed system since it shows damage effects due to shear deformation while it plastifies under pressure. Firstly, it would be of immense interest if strategies like the local yield stress method could be directly transferable to network glasses. The two-dimensional silica model glass presented in this paper seems to be an excellent candidate for such investigations. Secondly, applying the local yield stress to densified network glasses could provide valuable insights into the local effect of the pressure on the yielding behavior of such materials. The literature suggests that silica glass reveals two different shear transformation pictures: a shape-changing and a bond-breaking event [153, 221]. Regarding our simplified model network material in two dimensions, we observed only bond-breaking events and events that also include bond-forming processes. Considering this local network information, it would be intriguing if it is possible to develop a unified theory of fracture and plasticity. The bond rearrangement in the material and the prediction goodness from Patinet’s local probing method [47, 176] may provide valuable information for further use in a multiscale model that unifies plasticity and fracture.
11.2 From Prediction to Multiscale Modeling
The transportation of microscopic mechanics into larger scales constitutes another Holy Grail in modeling disordered solids. In this context, very recently, Castellanos et al. [206] developed elasto-plastic models quantified from atomistic Lennard-Jones glasses. They extended these models using the local yield stress method to extract the local mechanical response and maintain information regarding the microscale physics [207]. A realization of this approach to develop mesoscale models of network glasses would greatly enrich the mechanical understanding of glass fracture. The objective is to transport atomistic-informed local behavior to better model the combination of damage and plasticity around the crack tip of network glasses. This approach requires the successfully adapted local probing method, such as the above-mentioned local yield stress method, to predict local rearrangement events in network materials.
However, scanning a whole molecular cell to extract local yield stress quantities turns out to be challenging in terms of computational effort since it requires numerous local molecular mechanical investigations. To realize such a computational burden, various tools are available in the literature, including but not limited to machine learning methodologies and generic algorithms, model order reduction methods, et cetera. New methodologies to numerically speed up mechanical investigations of disordered solids would enrich current local scanning approaches in this context.
An intriguing mesoscale framework was published by Rycroft and Bouchbinder [279, 280], who presented a numerical realization of a low-temperature shear transformation zone model [14, 281,282,283,284,285]. In doing so, they were able to model an elastoplastic crack-tip instability using an Eulerian framework. It would be intriguing to extend the low-temperature shear transformation zone theory and apply it to network glasses. Such mesoscale simulations would then be performed on the continuum level. Pursuing such an objective may not be straightforward since network glasses experience damage when subjected to shear and plastify when subjected to pressure, as discussed elaborately in this paper.
11.3 Reduced Order Models and Metamodeling
The computational efficiency of molecular modeling has continuously improved in the last decades, as discussed in Section 6. In particular, parallelization, which was not the focus of this paper, allows one to investigate large models when calculating the problem on many processors simultaneously. In this context, one might argue that computational strategies to make molecular calculations efficient are pretty advanced. However, reduced-order models, mainly targeted toward describing the mechanical behavior of disordered solids, are, to date, absent from the literature. The thermal, high-frequency vibration is hard to capture in a reduced-order model. However, approaching catastrophes, the athermal quasistatic description of disordered solids shows very localized distinct response patterns. In particular, the material responds mainly in terms of the lowest vibration modes, also called soft modes [202, 274]. This property opens up promising perspectives for transformations into low-order subspaces. In this context, the main challenge may be to track changes in the low-vibrational soft modes [286]. As a result, the challenge is to predict the changes in the PEL of the highly nonlinear system during the deformation protocol.
Another promising strategy may be the application of machine learning to generate surrogate models in order to describe the complex stress-strain behavior of disordered systems. In particular, smaller configurations that occur when performing local probing, such as the local yield stress method [47, 176], may be identifiable by neural networks, for example. The configuration coordinates serve as input parameters, while the desired prediction quantity is the output parameter. Since it was shown that neural networks could describe very complex relationships, the first applications of neural networks were performed in relation to the prediction of inelastic events [215, 278]. It is noteworthy that the variety of possible machine learning approaches turns out to be extensively rich, and it is continuously growing. In this context, physics-informed neural networks have lately been presented, whereby the objective function to be minimized during the training algorithm is extended by physical information from the partial differential equations or the equilibrium equations, such as residual forces in the molecular system [287].
11.4 Statistical Models
Numerous mechanical processes of disordered solids on the nanoscale show power-law distributions and therefore follow self-organized critical behavior, which suggests that they are scale-invariant. One may use this property to generate a probabilistic pendant to the more deterministic approaches suggested in Sect. 11.2. Self-organized criticality could offer a bridge across the scales to enable the development of physically motivated engineering models of disordered solids on the meso- and macroscale in order to describe their complex yielding and fracture behavior. However, using purely continuum-mechanical frameworks, representing the typical serrated flow picture that governs the concept of self-organized criticality, may not be possible. One may introduce elastoplastic models representing this statistical framework due to stick-slip friction descriptions and incorporate these models into a continuum-mechanical framework. The focus should then be on capturing statistical behavior and reproducing such phenomena qualitatively on coarser scales [4]. Comparisons and validations with experimental results are necessary to provide reliable material models of disordered systems.
Concerning a macroscopic description of fracture, incorporating statistical information on microscopic slip motion into phase-field models could lead to a new class of physics-informed multiscale fracture models. In this context, we note that although both types of materials, i.e., densely-packed systems, such as metallic glasses, and network materials, such as silicate glasses, follow similar concepts of disorder, their mechanical behavior on the nanoscale and the macroscale turns out to be fundamentally different. Therefore, it may still be crucial to transfer statistical information and local response identification from the nano- to the macroscale when performing statistical multiscale approaches of such materials.
All things considered, a profound understanding of the disorder in materials and structures will foster the engineering process of structural and electro-mechanical industrial components and provide valuable insights into phenomena on different lengths- and time scales, such as the mechanics of tectonic plates.
References
Varshneya AK, Mauro JC (2019) Fundamentals of inorganic glasses, 3rd edn. Elsevier Inc., Amsterdam
Carlson JM, Langer JS, Shaw BE (1994) Dynamics of earthquake faults. Rev Mod Phys 66:657–671
Kawamura H, Hatano T, Kato N, Biswas S, Chakrabarti K (2012) Statistical physics of fracture, friction, and earthquakes. Rev Mod Phys 84:839–884
Nicolas A, Ferrero E, Martens K, Barrat J (2018) Deformation and flow of amorphous solids: insights from elastoplastic models. Rev Mod Phys 90:045006
Tanguy A (2021) Elasto-plastic behavior of amorphous materials: a brief review. C R Physique. https://doi.org/10.5802/crphys.49
Schuh C, Hufnagel T, Ramamurty U (2007) Mechanical behavior of amorphous alloys. Acta Mater 55:4067–4109
Trexler MM, Thadhani NN (2010) Mechanical properties of bulk metallic glasses. Prog Mater Sci 55:759–839
Egami T, Iwashita T, Dmowski W (2013) Mechanical properties of metallic glasses. Metals 3:77–113
Chen M (2008) Mechanical behavior of metallic glasses: microscopic understanding of strength and ductility. Annu Rev Matter Res 38:445–469
Wang WH (2012) The elastic properties, elastic models and elastic perspectives of metallic glasses. Prog Mater Sci 57:487–656
Greer AL, Cheng YQ, Ma E (2013) Shear bands in metallic glasses. Mater Sci Eng R Rep 74:71–132
Hufnagel T, Schuh C, Falk M (2016) Deformation of metallic glasses: recent developments in theory, simulations, and experiments. Acta Mater 109:375–393
Muralidharan K, Oh K, Deymier P, Runge K, Simmons J (2007) Molecular dynamics simulation of atomic-level brittle fracture mechanisms in amorphous silica. J Mater Sci 42:4159–4169
Falk ML, Langer JS (2011) Deformation and failure of amorphous, solidlike materials. Annu Rev Condens Matter Phys 2:353–373
Rodney D, Tanguy A, Vandembroucq D (2011) Modeling the mechanics of amorphous solids at different length scale and time scale. Model Simul Mater Sci Eng 19:083001
Bezanson J, Edelman A, Karpinski S, Shah VB (2017) Julia: a fresh approach to numerical computing. SIAM Rev Soc Ind Appl Math 59:65–98
Zanotto E, Mauro J (2017) The glassy state of matter: its definition and ultimate fate. J Non-Cryst Solids 471:490–495
Humphrey W, Dalke A, Schulten K (1996) VMD: visual molecular dynamics. J Mol Graph 14:33–38
Vashishta P, Kalia R, Rino J, Ebbsjö I (1990) Interaction potential for \(\text{ SiO}_{2}\): a molecular-dynamics study of structural correlations. Phys Rev B 41:12197–12210
Schmelzer JWP, Tropin TV (2018) Glass transition. Crystallization of glass-forming melts, and entropy. Entropy 20:103
Zanotto ED, Mauro JC (2018) Comment on “Glass transition, crystallization of glass-forming melts, and entropy’’ Entropy 2018, 20, 103. Entropy 20:703
Schmelzer JWP, Tropin TV (2018) Reply to “Comment on Glass transition. Crystallization of glass-forming melts, and entropy” by Zanotto and Mauro. Entropy 20:704
Raabe D (1998) Computational materials science: the simulation of materials. Microstructures and properties, Wiley, Weinheim
Gupta PK (1996) Non-crystalline solids: glasses and amorphous solids. J Non-Cryst Solids 195:158–164
Cheng YT, Johnson WL (1987) Disordered materials: a survey of amorphous solids. Science 235:997–1002
Leimkuhler B, Matthews C (2015) Molecular dynamics. Springer, Switzerland, With Deterministic and Stochastic Numerical Methods
Lelièvre T, Stoltz G (2016) Partial differential equations and stochastic methods in molecular dynamics. Acta Numer 25:681–880
Jones JE (1997) On the determination of molecular fields.–I. From the variation of the viscosity of a gas with temperature. Proc R Soc Lond A 106:441–462
Jones JE (1997) On the determination of molecular fields.–II. From the equation of state of a gas. Proc R Soc Lond A 106:463–477
Rappaport DC (2004) The art of molecular dynamics simulation, 2nd edn. Cambridge University Press, Cambridge
Allen MP, Tildesley DJ (1989) Computer simulation of liquids. Oxford University Press, New York
Ewald PP (1921) Die Berechnung optischer und elektrostatischer Gitterpotentiale. Ann Phys 369:253–287
Darden T, York D, Pedersen L (1993) Particle mesh Ewald: An N.log(N) method for Ewald sums in large systems. J Chem Phys 98:10089
Essmann U, Perera L, Berkowitz ML (1995) A smooth particle mesh Ewald method. J Chem Phys 103:8577
Stamm B, Lagardère L, Pollack È, Maday Y, Piquemal JP (2018) A coherent derivation of the Ewald summation for arbitrary orders of multipoles: the self-terms. J Chem Phys 149:124103
Lees AW, Edwards SF (1972) The computer study of transport processes under extreme conditions. J Phys C 5:1921–1929
Lemaître A, Maloney C (2006) Sum rules for the quasi-static and visco-elastic response of disordered solids at zero temperature. J Stat Phys 123:415–453
Kobayashi H, Yamamoto R (2011) Implementation of Lees-Edwards periodic boundary conditions for direct numerical simulations of particle dispersions under shear flow. J Chem Phys 134:064110
Ray JR, Moody MC, Rahman A (1985) Molecular dynamics calculation of elastic constants for a crystalline system in equilibrium. Phys Rev B 32:733–735
Ray JR, Rahman A (1984) Statistical ensembles and molecular dynamics studies of anisotropic solids. J Chem Phys 80:4423–4428
Holzapfel GA (2001) Nonlinear solid mechanics, a continuum approach for engineering. Wiley, LTD, Chichester
Baẑant ZP (1998) Easy-to-compute tensors with symmetric inverse approximating Hencky finite strain and its rate. J Eng Mater Technol 120:131–136
Shrivastava S, Ghosh C, Jonas JJ (2012) A comparison of the von Mises and Hencky equivalent strains for use in simple shear experiments. Philos Mag 92:779–786
Jonas JJ, Ghosh C, Basabe V, Shrivastava S (2012) The Hencky equivalent strain and its inapplicability to the interpretation of torsion testing experiments. Philos Mag 92:2313–2328
Rondey D, Schuh CA (2009) Yield stress in metallic glasses: the jamming-unjamming transition studied through Monte Carlo simulations based on the activation-relaxation technique. Phys Rev B 80:184203
Destrade M, Murphy JG, Saccomandi G (2012) Simple shear is not so simple. Int J Non Linear Mech 47:210–214
Patinet S, Vandembroucq D, Falk ML (2016) Connecting local yield stresses with plastic activity in amorphous solids. Phys Rev Lett 117:045501
Frenkel D, Smit B (2001) Understanding molecular simulation: from algorithms to application. Academic Press, San Diego
Milkus R, Zaccone A (2016) Local inversion-symmetry breaking controls the boson peak in glasses and crystals. Phys Rev B 93:094204
Milkus R, Zaccone A (2017) Atomic-scale origin of dynamic viscoelastic response and creep in disordered solids. Phys Rev E 95:023001
Cui B, Zaccone A, Rodney D (2019) Nonaffine lattice dynamics with the Ewald method reveals strongly nonaffine elasticity of \(\alpha\)-quartz. J Chem Phys 151:224509
Focks T, Bamer F, Markert B, Wu Z, Stamm B (2022) Displacement field splitting of defective hexagonal lattices. Phys Rev B 106:014105
Scagnetti PA, Nagem RJ, Sandri GVH, Bifano TG (1996) Stress and strain analysis in molecular dynamics simulation of solids. ASME J Appl Mech 63:450–452
Clausius R (1870) On a mechanical theorem applicable to heat. Philos Mag 40:122–127
Tsai DH (1979) The virial theorem and stress calculation in molecular dynamics. J Chem Phys 70:1375–1382
Hansen J, McDonald I (1986) Theory of simple liquids. Academic Press, London
Plimpton S (1995) Fast parallel algorithms for short-range molecular dynamics. J Comput Phys 117:1–17
Thompson A, Plimpton S, Mattson W (2009) General formulation of pressure and stress tensor for arbitrary many-body interaction potentials under periodic boundary conditions. J Chem Phys 131:154107–154113
Vitek V, Egami T (1987) Atomic level stresses in solids and liquids. Phys Stat Sol B 144:145–156
Cheung KS, Yip S (1991) Atomic-level stress in an inhomogeneous system. J Appl Phys 70:5688–5690
Yang J, Komvopoulos K (2020) A stress analysis method for molecular dynamics systems. Int J Solids Struct 193:98–105
Mesirov J, Schulten K, Sumners D (1996) Mathematical approaches to biomolecular structure and dynamics. Springer, New York
Verlet L (1967) Computer “Experiments’’ on classical fluids. I. Thermodynamical properties of Lennard-Jones molecules. Phys Rev 159:98–103
DiTolla F, Ronchetti M (1993) Applicability of Nosé isothermal reversible dynamics. Phys Rev E 48:1726
Hünenberger P (2005) Thermostat algorithms for molecular dynamics simulations. Adv Polym Sci 173:105–149
Melchionna S, Ciccotti G, Holian B (1993) Hoover NPT dynamics for systems varying in shape and size. Mol Phys 78:533–544
Nosé S (1984) A unified formulation of the constant temperature molecular dynamics methods. J Chem Phys 81:511
Hoover W (1985) Canonical dynamics: equilibrium phase-space distributions. Phys Rev A 31:1695
Martyna G, Klein M, Tuckerman M (1992) Nosé-Hoover chains: the canonical ensemble via continuous dynamics. J Chem Phys 97:2635–2643
Tanaka H, Nakanishi K, Watanabe N (1983) Constant temperature molecular dynamics calculation on Lennard-Jones fluid and its application to water. J Chem Phys 78:2626
Tuckerman ME, Alejandre J, López-Rendón R, Jochim AL, Martyna GJ (2006) A Liouville-operator derived measure-preserving integrator for molecular dynamics simulations in the isothermal-isobaric ensemble. J Phys A 39:5629–5651
Branicio P, Rino J, Gan C, Tsuzuki H (2009) Interaction potential for indium phosphide: a molecular dynamics and first-principles study of the elastic constants, generalized stacking fault and surface energies. J Phys 21:095002
Mattson W, Rice B (1999) Near-Neighbor calculations using a modified cell-linked list method. Comput Phys Commun 119:135–148
Nocedal J, Wright S (2006) Numerical optimization. Springer, New York
Hestenes MR, Stiefel E (1952) Methods of conjugate gradients for solving linear systems. J Res Natl Bur Stand 49:409–436
Mauro JC, Zanotto ED (2014) Two centuries of glass research: historical trends, current status, and grand challenges for the future. Int J Appl Glass Sci 5:313–327
Mauro J, Smedskjaer M (2014) Statistical mechanics of glass. J Non-Cryst Solids 396:41–53
Lancon F, Billard L, Chaudhari P (1986) Thermodynamical properties of a two-dimensional quasi-crystal from molecular dynamics calculations. Europhys Lett 2:625–629
Falk M, Langer J (1998) Dynamics of viscoplastic deformation in amorphous solids. Phys Rev E 57:7192
Maloney C, Lemaître A (2004) Universal breakdown of elasticity at the onset of material failure. Phys Rev Lett 93:195501
Maloney C, Lemaître A (2004) Subextensive scaling in the athermal, quasistatic limit of amorphous matter in plastic shear flow. Phys Rev Lett 93:016001
Maloney C, Lemaître A (2006) Amorphous systems in athermal quasistatic shear. Phys Rev E 74:016118
Beest BV, Kramer G, van Santen R (1990) Force fields for silicas and aluminophosphates based on ab initio calculations. Phys Rev Lett 64:1955–1958
Smedskjaer M, Bauchy M, Mauro J, Sylwester J, Bockowski M (2015) Unique effects of thermal and pressure histories on glass hardness: structural and topological origin. J Chem Phys 143:164505
Comin M, Levis LJ (2019) Deep-learning approach to the structure of amorphous silicon. Phys Rev B 100:094107
Bidault X, Chaussedent S, Blanc W, Neuville D (2016) Deformation of silica glass studied by molecular dynamics: structural origin of the anisotropy and non-Newtonian behavior. J Non-Cryst Solids 433:38–44
Mozzi R, Warren B (1969) The structure of vitreous silica. J Appl Crystallogr 2:164–172
Pettifer R, Dupree R, Farnan I, Sternberg U (1988) NMR determinations Si-O-Si bond angle distributions in silica. J Non-Cryst Solids 106:408–412
Ebrahem F, Bamer F, Markert B (2018) The influence of the network topology on the deformation and fracture behaviour of silica glass: a molecular dynamics study. Comput Mater Sci 149:162–169
Koziatek P, Barrat J, Rodney D (2015) Short- and medium-range orders in as-quenched and deformed SiO\(_2\) glasses: an atomistic study. J Non-Cryst Solids 414:7–15
Vollmayr K, Kob W, Binder K (1996) Cooling-rate effects in amorphous silica: a computer-simulation study. Phys Rev B 54:15808
King S (1967) Ring Configurations in a random network model of vitreous silica. Nature 213:1112–1113
Thorpe M (1983) Continuous deformations in random networks. J Non-Cryst Solids 57:355–370
Guttman L (1990) Ring structure of the crystalline and amorphous forms of silicon dioxide. J Non-Cryst Solids 116:145–147
Bondy JA, Murty USR (1976) Graph theory with applications. Macmillan, London
Franzblau D (1991) Computation of ring statistics for network models of solids. Phys Rev B 44:4925–4930
Hobbs L, Jesurum C, Pulim V, Berger B (1998) Local topology of silica networks. Philos Mag A 78:679–711
Yuan X, Cormack A (2002) Efficient algorithm for primitive ring statistics in topological networks. Comput Mater Sci 24:343–360
Roux SL, Jund P (2010) Ring statistics analysis of topological networks: new approach and application to amorphous GeS2 and SiO2 systems. Comput Mater Sci 49:70–83
Bamer F, Ebrahem F, Markert B (2019) Plasticity in vitreous silica induced by cyclic tension considering rate-dependence: role of the network topology. J Non-Cryst Solids 503:176–181
Zachariasen W (1932) The atomic arrangement in glass. J Am Chem Soc 54:3841–3851
Lichtenstein L, Büchner C, Yang B, Shaikhutdinov S, Heyde M, Sierka M, Włodarczyk R, Sauer J, Freund H (2012) The atomic structure of a metal-supported vitreous thin silica film. Angew Chem Int Ed 51:404–407
Lichtenstein L, Heyde M, Freund H (2012) Atomic arrangement in two-dimensional silica: from crystalline to vitreous structures. J Phys Chem C 116:20426–20432
Lichtenstein L, Heyde M, Freund H (2012) Crystalline-vitreous interface in two dimensional silica. Phys Rev Lett 109:106101
Heyde M, Shaikhutdinov S, Freund H (2012) Two-dimensional silica: crystalline and vitreous. Chem Phys Lett 550:1–7
Shaikhutdinov S, Freund H (2013) Ultrathin silica films on metals: the long and winding road to understanding the atomic structure. Adv Mater 25:49–67
Büchner C, Liu L, Stuckenholz S, Burson K, Lichtenstein L, Heyde M, Gao H, Freund H (2016) Building block analysis of 2D amorphous networks reveals medium range correlation. J Non-Cryst Solids 435:40–47
Büchner C, Lichtenstein L, Yu X, Boscoboinik J, Yang B, Kaden W, Heyde M, Shaikhutdinov S, Włodarczyk R, Sierka M, Sauer J, Freund H (2014) Ultrathin silica films: the atomic structure of two-dimensional crystals and glasses. Chem Eur J 20:9176
Büchner C, Heyde M (2017) Two-dimensional silica opens new perspectives. Prog Surf Sci 92:341–374
Gao E, Xie B, Xu Z (2016) Two-dimensional silica: structural, mechanical properties, and strain-induced band gap tuning. J Appl Phys 119:014301
Gao Z, Dong X, Li N, Rien J (2017) Novel two-dimensional silicon dioxide with in-plane negative Poisson’s ratio. Nano Lett 17:772
Safaei S, Tavakoli R, Jafary-Zadeh M (2018) Molecular dynamics study of two dimensional silicon dioxides with in-plane negative Poisson’s ratio. Comput Mater Sci 153:258
Zhang J (2018) Phase-dependent mechanical properties of two-dimensional silica films: a molecular dynamics study. Comput Mater Sci 142:7
Bamer F, Ebrahem F, Markert B (2020) Elementary plastic events in a Zachariasen glass under shear and pressure. Materialia 9:100556
Roy P, Heyde M, Heuer A (2018) Modelling the atomic arrangement of amorphous 2D silica: a network analysis. Phys Chem Chem Phys 20:14725–14739
Roy P, Heuer A (2019) Ring stastitics in 2D silica: effective temperatures in equilibrium. Phys Rev Lett 122:016104
Roy P, Heuer A (2019) Relating local structures, energies, and occurrence probabilities in a two-dimensional silica network. J Phys 31:225703
Huang P, Kurasch S, Srivastava A, Skakalova V, Kotakoski J, Krasheninnikov A, Hovden R, Mao Q, Meyer J, Smet J, Muller K, Kaiser U (2012) Direct imaging of a two-dimensional silica glass on graphene. Nano Lett 12:1081–1086
Huang PY, Kurasch S, Alden JS, Shekhawat A, Alemi AA, McEuen PL, Sethna JP, Kaiser U, Muller DA (2013) Imaging atomic rearrangements in two-dimensional silica glass: watching silica’s dance. Science 342:224–227
Godréche C, Kostov I, Yekutieli I (1992) Topological correlations in cellular structures and planar graph theory. Phys Rev Lett 69:2674–2677
Heyde M (2013) Structure and motion of a 2D glass. Science 342:201–202
Metropolis N, Rosenbluth AW, Rosenbluth MN, Teller AH, Teller E (1953) Equation of state calculations by fast computing machines. J Chem Phys 21:1087–1092
Hastings WK (1970) Monte Carlo sampling methods using Markov chains and their applications. Biometrika 57:97–109
Kumar A, Wilson M, Thorpe M (2012) Amorphous graphene: a realization of Zachariasen’s glass. J Phys 24:485003
Wilson M, Kumar A, Sherrington D, Thorpe M (2013) Modelling vitreous silica bilayers. Phys Rev B 87:214108
Novoselov K, Geim A, Morozov S, Jiang D, Zhang Y, Dubonos S, Grigorieva I, Firsov A (2004) Electric field effect in atomically thin carbon films. Science 306:666–669
Meyer J, Geim A, Katsnelson M, Novoselov K, Booth T, Roth S (2007) The structure of suspended graphene sheets. Nature 446:60–63
Geim A, Novoselov K (2007) The rise of graphene. Nature Mater. 6:183–191
Akinwande D, Brennan C, Bunch J, Egberts P, Felts J, Gao H, Huang R, Kim J, Li T, Li Y, Liechti K, Lu N, Park H, Reed E, Wang P, Yakobson B, Zhang T, Zhang Y, Zhou Y, Zhu Y (2017) A review on mechanics and mechanical properties of 2D materials–graphene and beyond. Extreme Mech. Lett. 13:42–72
Grima J, Winczewski S, Mizzi L, Grech M, Cauchi R, Gatt R, Attard D, Wojciechowski K, Rybicki J (2015) Tailoring graphene to achieve negative Poisson’s ratio properties. Adv Mat 27:1455–1459
Banhart F, Kotakoski J, Krasheninnikov A (2010) Structural defects in graphene. ASC. Nano 5:26–41
He L, Guo S, Lei J, Sha Z, Liu Z (2014) The effect of Stone-Thrower-Wales defects on mechanical properties of graphene sheets–a molecular dynamics study. Carbon 75:124–132
Rajasekaran G, Parashar A (2016) Molecular dynamics study on the mechanical response and failure behaviour of graphene: perfomance enhancement via 5-7-7-5 defects. RSC Adv 6:26361
Hossain M, Ahmed Z, Silverman B, Khawaja M, Calderon J, Rutten A, Tse S (2018) Anisotropic toughness and strength in graphene and its atomistic origin. J Mech Phys Solids 110:118–136
Joo W, Lee J, Jang Y, Kang S, Kwon Y, Chung J, Lee S, Kim C, Kim T, Yang C, Kim U, Choi B, Whang D, Hwang S (2017) Realization of continuous Zachariasen carbon monolayer. Sci Adv 3:e1601821
Toh C, Zhang H, Lin J, Mayorov A, Wang Y, Orofeo C, Ferry D, Andersen H, Kakenov N, Guo Z, Abidi I, Sims H, Suenaga K, Pantelides S, Özyilmaz B (2020) Synthesis and properties of free-standing monolayer amorphous carbon. Nature 577:199–203
Ebrahem F, Stratmann J, Stoffel M, Markert B, Bamer F (2020) Continuous Zachariasen carbon monolayers under tensile deformation: Insights from molecular dynamics simulations. Extreme Mech. Lett. 38:100744
Kotakoski J, Krasheninnikov A, Kaiser U, Meyer J (2011) From point defects in graphene to two-dimensional amorphous carbon. Phys Rev Lett 106:105505
Xu M, Liang T, Shi M, Chen H (2013) Graphene-like two-dimensional materials. Chem Rev 113:3766–3798
Partovi-Azar P, Jand S, Namiranian A, Rafii-Tabar H (2013) Electronic features induced by Stone-Wales defects in zigzag and chiral carbon nanotubes. Comput Mater Sci 79:82–86
Fu Y, Ragab T, Basaran C (2016) The effect of Stone-Wales defects on the mechanical behavior of graphene nano-ribbons. Comput Mater Sci 124:142–150
Wang S, Fan Z, Cui Y, Zhang S, Yang B, Chen H (2017) Fracture behaviors of brittle to ductile 2D carbon structures under uniaxial tensile stress. Carbon 111:486–492
Yang G, Li L, Lee W, Ng M (2018) Structure of graphene and its disorders: a review. Sci Technol Adv Mater 19:613–648
Ansari R, Ajori S, Motevalli B (2012) Mechanical properties of defective single-layered graphene sheets via molecular dynamics simulation. Superlattices Microstruct 51:274
Openov L, Podlivaev A (2015) Interaction of the Stone-Wales defects in graphene. Phys Solid State 57:1477
Wang S, Yang B, Yuan J, Si Y, Chen H (2015) Large-scale molecular simulations on the mechanical response and failure behavior of a defective graphene: cases of 5-8-5 defects. Sci Rep 5:14957
Morley D, Wilson M (2018) Controlling disorder in two-dimensional networks. J Phys 30:50LT02
Aboav D (1970) The arrangement of grains in a polycrystal. Metallography 3:383–390
Aboav D (1980) The arrangement of cells in a net. Metallography 13:43–58
Weaire D (1974) Some remarks on the arrangement of grains in a polycrystal. Metallography 7:157–160
Born M, Kármán T (1912) Über Schwingungen in Raumgittern. Phys. Zeit. 8:297–309
Bathe KJ (2006) Finite element procedures, 2nd edn. Prentice Hall, Pearson Education, Watertown
Bonfanti S, Guerra R, Mondal C, Procaccia I, Zapperi S (2019) Elementary plastic events in amorphous silica. Phys Rev E 100:060602(R)
Lanczos C (1988) Applied analysis. Dover, New York
Lerner E (2016) Micromechanics of nonlinear plastic modes. Phys Rev E 93:053004
Ghysels A, Van Neck D, Van Speybroeck V, Verstraelen T, Waroquier M (2007) Vibrational modes in partially optimized molecular systems. J Chem Phys 126:224102
Ashby MF (1972) A first report on deformation-mechanism maps. Acta Metall 20:887–897
Jones RM (2009) Deformation theory of plasticity. Bull Ridge Publishing, Blacksburg
Li S, Ren B, Minaki H (2014) Multiscale crystal defect dynamics: a dual-lattice process zone model. Philos Mag 94:1414–1450
Lancon F, Billard L, Laugier J, Chamberot A (1982) Simulation of a reproducible model of metallic glasses by hardsphere relaxation. J Phys F 12:259–268
Lancon F, Billard L, Chamberot A (1984) Structural description of a metallic glass model. J Phys F 14:579–591
Jin Y, Urbani P, Zamponi F, Yoshino H (2018) A stability-reversibility map unifies elasticity, plasticity, yielding, jamming, in hard sphere glasses. Sci Adv 4:eaat6387
Spaepen F (1977) A microscopic mechanism for steady state inhomogeneous flow in metallic glasses. Acta Metall 25:407–415
Barbot A, Lerbinger M, Lemaïtre A, Vandembroucq D, Patinet S (2020) Rejuvenation and shear banding in model amorphous solids. Phys Rev E 101:033001
Argon A, Kuo H (1979) Plastic flow in a disordered bubble raft (an analogy of a metallic glass). Mater Sci Eng 39:101–109
Argon A, Shi L (1982) Analysis of plastic flow in amorphous soap bubble raft by the use of an inter-bubble potential. Philos Mag A 46:275–294
Argon A, Shi L (1983) Development of visco-plastic deformation in metallic glasses. Acta Metall 31:499–507
Argon A (1979) Plastic deformation in metallic glasses. Acta Metall 27:47–58
Weeks ER, Crocker JC, Levitt AC, Schofield A, Weitz DA (2000) Three-dimensional direct imaging of structural relaxation near the colloidal glass transition. Science 287:627–631
Schall P, Weitz D, Spaepen F (2007) Structural rearrangements that govern flow in colloidal glasses. Science 318:1895–1899
Jensen K, Weitz D, Spaepen F (2014) Local shear transformations in deformed and quiescent hard-sphere colloidal glasses. Phys Rev E 90:042305
Eshelby J (1957) The determination of the elastic field of an ellipsoidal inclusion, and related problems. Proc R Soc A 241:376–396
Falk ML (1999) Molecular-dynamics study of ductile and brittle fracture in model noncrystalline solids. Phys Rev B 60:7062–7070
Langer JS (2006) Shear-transformation-zone theory of deformation in metallic glasses. Scr Mater 54:375–379
Gendelman O, Jaiswal P, Procaccia I, Gupta B, Zylberg J (2013) Shear Transformation Zones: state determined or protocol dependent? Europhys Lett 88:022310
Barbot A, Lerbinger M, Hernandez-Garcia A, García-García R, Falk ML, Vandembroucq D, Patinet S (2018) Local yield stress statistics in model amorphous solids. Phys Rev E 97:033001
Zaccone A, Schall P, Ternetjev E (2014) Microscopic origin of nonlinear nonaffine deformation in bulk metallic glasses. Phys Rev B 90:140203(R)
Wittmer JP, Tanguy A, Barrat JL, Lewis L (2002) Vibrations of amorphous, nanometric structures: when does continuum theory apply? Europhys Lett 57:423–429
Tanguy A, Wittmer JP, Leonforte F, Barrat JL (2002) Continuum limit of amorphous elastic bodies: a finite-size study of low-frequency harmonic vibration. Phys Rev B 66:174205
Dasgupta R, Karmakar S, Procaccia I (2012) Universality of the plastic instability in strained amorphous solids. Phys Rev Lett 108:075701
Shi Y, Falk M (2006) Does metallic glass have a backbone? The role of percolating short range order in strength and failure. Scr Mater 54:381–386
Shi Y, Falk M (2007) Stress-induced structural transformation and shear banding during simulated nanoindentation of a metallic glass. Acta Mater 55:4317–4324
Karimi K, Barrat JL (2018) Correlation and shear bands in a plastically deformed granular medium. Sci Rep 8:4021
Dasgupta R, Hentschel HGE, Procaccia I (2012) Microscopic mechanism of shear bands in amorphous solids. Phys Rev Lett 109:255502
Dasgupta R, Hentschel HGE, Procaccia I (2013) Yield strain in shear banding amorphous solids. Phys Rev E 87:022810
Lemaitre A, Caroli C (2007) Plastic response of a two-dimensional amorphous solid to quasistatic shear: transverse particle diffusion and phenomenology of dissipative events. Phys Rev E 76:036104
Hinkle AR, Rycroft CH, Shields MD, Falk ML (2017) Coarse graining atomistic simulations of plastically deforming amorphous solids. Phys Rev E 95:053001
Hieronymus-Schmidt V, Rösner H, Wilde G, Zaccone A (2017) Shear banding in metallic glasses described by alignments of Eshelby quadrupoles. Phys Rev B 95:134111
Ott RT, Sansoz F, Jiao T, Warner D, Molinari JF, Ramesh KT, Hufnagel TC, Fan C (2006) Yield criteria and strain-rate behavior of Zr\(_{57.4}\)Cu\(_{16.4}\)Ni\(_{8.2}\)Ta\(_{8}\)Al\(_{10}\) metallic-glass-matrix composites. Metall Mater Trans A 3251–3258
Ashwin J, Gendelman O, Procaccia I, Shor C (2013) Yield strain and shear-band direction in amorphous solids under two-dimensional uniaxial loading. Phys Rev E 88:022310
Ebrahem F, Bamer F, Markert B (2020) Origin of reversible and irreversible atomic-scale rearrangements in a model two-dimensional network glass. Phys Rev E 102:033006
Albaret T, Tanguy A, Boioli F, Rodney D (2016) Mapping between atomistic simulations and Eshelby inclusions in the shear deformation of an amorphous silicon model. Phys Rev E 93:053002
Boioli F, Albaret T, Rodney D (2017) Shear transformation distribution and activation in glasses at the atomic scale. Phys Rev E 95:033005
Tsamados M, Tanguy A, Goldenberg C, Barrat JL (2009) Local elasticity map and plasticity in a model Lennard-Jones glass. Phys Rev E 80:026112
Jacks RL, Dunleavy AJ, Royall CP (2014) Information-theoretic measurements of coupling between structure and dynamics in glass formers. Phys Rev Lett 113:095703
Richard D, Ozawa M, Patinet S, Stanifer E, Shang B, Ridout SA, Xu B, Zhang G, Morse PK, Barrat JL, Berthier L, Falk ML, Guan P, Liu AJ, Martens K, Sastry S, Vandembroucq D, Lerner E, Manning ML (2020) Predicting plasticity in disordered solids from structural indicators. Phys Rev Mater 4:113609
Rieser JM, Goodrich CP, Liu AJ, Durian DJ (2016) Divergence of Voronoi cell anisotropy vector: a threshold-free characterization of local structure in amorphous materials. Phys Rev Lett 116:088001
Ding J, Cheng YQ, Sheng H, Asta M, Richie RO, Ma E (2017) Universal structural parameter to quantitatively predict metallic glass properties. Nat Commun 7:13733
Baggioli M, Kriuchevskyi I, Sirk TW, Zaccone A (2021) Plasticity in amorphous solids is mediated by topological defects in the displacement field. Phys Rev Lett 127:015501
Manning ML, Liu AJ (2011) Vibrational modes identify soft spots in a sheared disordered packing. Phys Rev Lett 107:108302
Rottler J, Schoenholz SS, Liu AJ (2014) Predicting plasticity with soft vibrational modes: from dislocations to glass. Phys Rev E 89:042304
Mosayebi M, Ilg P, Widmer-Cooper A, Del Gado E (2014) Soft modes and nonaffine rearrangements in the inherent structures of supercooled liquids. Phys Rev Lett 112:105503
Xu B, Falk ML, Patinet S, Guan P (2021) Atomic nonaffinity as a predictor of plasticity in amorphous solids. Phys. Rev. Mater. 5:025603
Ruan D, Patinet S, Falk ML (2022) Predicting plastic events and quantifying the local yield surface in 3D model glasses. J. Mech. Phys. Sol. 158:104671
Patinet S, Barbot A, Lerbinger M, Vandembroucq D, Lemaître A (2020) Origin of the Bauschinger effect in amorphous solids. Phys Rev Lett 124:205503
Castellanos DF, Roux S, Patinet S (2021) Insights from the quantitative calibration of an elasto-plastic model from a Lennard-Jones atomic glass. C R Phys 22:135–162
Castellanos DF, Roux S, Patinet S. History dependent plasticity of glass: a mapping between atomistic and elasto-plastic models. Acta Mater. https://doi.org/10.2139/ssrn.4015207
Barkema GT, Mousseau N (1996) Event-based relaxation of continuous disordered systems. Phys Rev Lett 77:4358–4361
Mousseau N, Barkema GT (1998) Traveling through potential energy landscapes of disordered materials: the activation-relaxation technique. Phys Rev E 57:2419–2424
Cances E, Legoll F, Marinica MC, Minoukadeh K, Willaime F (2009) Some improvements of the activation-relaxation technique method for finding transition pathways on potential energy surfaces. J Chem Phys 130:114711
Fan Y, Iwashita T, Egami T (2017) Energy landscape-driven non-equilibrium evolution of inherent structure in disordered material. Nat Commun 8:15417
Xu B, Falk M, Li J, Kong L (2018) Predicting shear transformation events in metallic glasses. Phys Rev Lett 120:125503
Mills G, Jónsson H (1994) Quantum and thermal effects in H2 dissociative adsorption: evaluation of free energy barriers in multidimensional quantum systems. Phys Rev Lett 72:1124–1127
Xu B, Falk M, Li J, Kong L (2017) Strain-dependent activation energy of shear transformation in metallic glasses. Phys Rev B 95:144201
Cubuk ED, Schoenholz SS, Rieser JM, Malone BD, Rottler J, Durian DJ, Kaxiras E, Liu AJ (2015) Identifying structural flow defects in disordered solids using machine-learning methods. Phys Rev Lett 114:108001
Ma Y, Guo G (2014) Support vector machines applications. Springer, Berlin
Wang Q, Ding J, Zhang L, Podryabinkin E, Shapeev A, Ma E (2020) Predicting the propensity for thermally activated \(\beta\)-events in metallic glasses via interpretable machine learning. NPJ Comput Mater 6:194
Wang Q, Jain A (2019) A transferable machine-learning framework linking interstice distribution and plastic heterogeneity in metallic glasses. Nat Commun 10:5537
Fan Z, Ma E (2021) Predicting orientation-dependent plastic susceptibility from static structure in amorphous solids via deep learning. Nat Commun 12:1506
Lundberg M, Krishan K, Xu N, O’Hern CS, Dennin M (2008) Reversible plastic events in amorphous materials. Phys Rev E 77:041505
Bamer F, Alshabab SS, Paul A, Ebrahem F, Markert B, Stamm B (2021) Data-driven classification of elementary rearrangement events in silica glass. Scr Mater 205:114179
Gagnon G, Patton J, Lacks D (2001) Energy landscape view of fracture and avalanches in disordered materials. Phys Rev E 64:051508
Lorenz CD, Stevens MJ (2003) Fracture behavior of Lennard-Jones glasses. Phys Rev E 68:021802
Anderson T (2005) Fracture mechanics: fundamentals and applications. Taylor & Francis Group, Boca Raton
Rountree C, Kalia R, Lidorikis E, Nakano A, Brutzel LV, Vashishta P (2002) Atomistic aspects of crack propagation in brittle materials: multimillion atom molecular dynamics simulations. Annu Rev Mater Res 32:377–400
Rountree C (2017) Recent progress to understand stress corrosion cracking in sodium borosilicate glasses: linking the chemical composition to structural, physical and fracture properties. J Phys D Appl Phys 50:343002
Pedone A (2009) Properties calculations of silica-based glasses by atomistic simulations techniques: a review. J Phys Chem C 113:20773–20784
Bauchy M (2014) Structural, vibrational, and elastic properties of a calcium alumina silicate glass from molecular dynamics simulations: the role of the potential. J Chem Phys 141:024507
Pedone A, Menziani M, Cormack A (2015) Dynamics of fracture in silica and soda-silicate glasses: from bulk materials to nanowires. J Phys Chem C 119:25499–25507
Smedskjaer MM, Bauchy M, Mauro JC, Sylvester MR, Bockovsky M (2015) Unique effects of thermal and pressure histories on glass hardness: structural and topological origin. J Chem Phys 143:164505
Yu Y, Anoop Krishnan NM, Smedskjaer MM, Sant G, Bauchy M (2018) The hydrophilic-to-hydrophobic transition of is driven by the topology of its surface. J Chem Phys 148:074503
Célarié F, Padres S, Bonamy D, Ferrero L, Bouchaud E, Guillot C, Marlièère C (2003) glass breaks like metal, but at the nanometer scale. Phys Rev Lett 90:075504
Lacroix R, Kermouche G, Teisseire J, Barthel E (2012) Plastic deformation and residual stresses in amorphous silica pillars under uniaxial loading. Acta Mater 60:5555–5566
Luo J, Wang J, Bitzek E, Huang Y, Zheng H, Tong L, Yang Q, Li J, Mao S (2015) Size-dependent brittle-to-ductile transition in silica glass nanofibers. Nano Lett 16:105–103
Guin JP, Wiederhorn SM (2004) Fracture of silicate glasses: ductile or brittle. Phys Rev Lett 92:215502
Patinet S, Vandembroucq D, Hansen A (2014) Cracks in random brittle solids. Eur. Phys. J.: Spec. Top. 223:2339–2351
Wang B, Yu Y, Wang M, Mauro J, Bauchy M (2016) Nanoductility in silicate glasses is driven by topological heterogeneity. Phys Rev B 93:064202
Ebrahem F, Bamer F, Markert B (2020) Vitreous 2D silica under tension: from brittle to ductile behaviour. Mater Sci Eng A 780:139189
Bamer F, Ebrahem F, Markert B (2019) Athermal mechanical analysis of Stone-Wales defects in two-dimensional silica. Comput Mater Sci 163:301–307
Ebrahem F, Bamer F, Markert B (2020) Stone-Wales defect interaction in quasistatically deformed 2D silica. J Mater Sci 55:3470–3483
Pedone A, Malavasi G, Menziani M, Serge U, AN C (2008) Molecular dynamics studies of stress-strain behavior of silica glass under a tensile load. Chem Mater 20:4356–4366
Chowdhury S, Haque B, Gillespie J (2016) Molecular dynamics simulations of the structure and mechanical properties of silica glass using ReaxFF. J Mater Sci 51:10139–10159
Kondo K, Lio S, Sawaoka A (1981) Nonlinear pressure dependence of the elastic moduli of fused quartz up to 3 GPa. J Appl Phys 52:2826–2831
Meade C, Jeanloz R (1988) Effect of a coordination change on the strength of amorphous \(\text{ SiO}_2\). Science 241:1072–1074
Mantisi B, Kermouche G, Barthel E, Tanguy A (2016) Impact of pressure on plastic yield in amorphous solids with open structure. Phys Rev E 93:033001
Jin W, Kalia R, Vashishta P (1994) Structural transformation in densified silica glass: a molecular-dynamics study. Phys Rev B 50:118–133
Tian Y, Du J, Hu D, Zheng W, Han W (2018) Densification effects on porous silica: A molecular dynamics study. Scr Mater 149:58–61
Yuan F, Huang L (2014) Brittle to ductile transition in densified silica glass. Sci Rep 4:5035
Mantisi B, Tanguy A, Kermouche G, Barthel E (2012) Atomic response of a model silica glass under shear and pressure. Eur Phys J B 85:304
Schill W, Heyden S, Contil S, Ortiz M (2018) The anomalous yield behavior of fused silica glass. J Mech Phys Solids 113:105–125
Wondraczek L, Mauro JC, Eckert J, Kühn U, Horbach J, Deubener J, Rouxel T (2011) Towards ultrastrong glasses. Adv Mater 23:4578–4586
Limbach R, Rodrigues BP, Wondraczek L (2014) Strain-rate sensitivity of glasses. J Non-Cryst Solids 404:124–134
Limbach R, Winterstein-Beckmann A, Dellith J, Möncke D, Wondraczek L (2015) Plasticity, crack initiation and defect resistance in alkali-borosilicate glasses: from normal to anomalous behavior. J Non-Cryst Solids 417:15–27
Kapoor S, Wondraczek L, Smedskjaer MM (2017) Pressure-induced densification of oxide glasses at the glass transition. Front Mater 4:1
Benzine O, Bruns S, Pan Z, Durst K, Wondraczek L (2018) Local deformation of glasses is mediated by rigidity fluctuation on nanometer scale. Adv Sci 5:1800916
Sawamura S, Wondraczek L (2018) Scratch hardness of glass. Phys. Rev. Mater. 2:092601
Wondraczek L, Bouchbinder E, Ehrlicher A, Mauro JC, Sajzew R, Smedskjaer MM (2022) Advancing the mechanical performance of glasses: perspectives and challenges. Adv Mater 34:2109029
Rodrigues BP, To T, Smedskjaer MM, Wondraczek L (2022) Mechanical properties of oxide glasses. Rev Mineral Geochem 87:229–281
Salerno KM, Maloney CE, Robbins MO (2012) Avalanches in strained amorphous solids: does inertia destroy critical behavior? Phys Rev Lett 109:105703
Salerno KM, Robbins MO (2013) Effect of inertia on sheared disordered solids: critical scaling of avalanches in two and three dimensions. Phys Rev E 88:062206
Clemmer JT, Salerno KM, Robbins MO (2021) Criticality in sheared, disordered solids I Rate effects in stress and diffusion. Phys Rev E 103:042605
Bak P, Tang C, Wiesenfeld K (1987) Self-organized criticality: an explanation of \(1/f\) noise. Phys Rev Lett 59:381–384
Hergarten S (2002) Self-organized criticality in earth systems. Springer, Berlin
Wattkins NW, Pruessner G, Chapman SC, Crosby NB, Jensen HJ (2016) 25 years of self-organized criticality: concepts and controversies. Space Sci Rev 198:3–44
Ferrero EE, Jagla EA (2019) Elastic interfaces on disordered substrates: from mean-field depinning to yielding. Phys Rev Lett 123:218002
Zhang D, Dahmen K, Ostoja-Starzewski M (2017) Scaling of slip avalanches in sheared amorphous materials based on large-scale atomistic simulations. Phys Rev E 95:032902
Liu C, Ferrero EE, Puosi F, Barrat JL, Martens K (2016) Driving rate dependence of avalanche statistics and shapes at the yielding transition. Phys Rev Lett 116:065501
Sun BA, Yu HB, Jiaon W, Bai HY, Zhao DQ, Wang WH (2010) Plasticity of ductile metallic glasses: a self-organized critical state. Phys Rev Lett 105:035501
Bonamy D, Santucci S, Ponson L (2008) Crackling dynamics in material failure as the signature of a self-organized dynamic phase transition. Phys Rev Lett 101:045501
Bak P, Creutz M (1994) Fractals and self-organized criticality. In: Bunde A, Havlin S (eds) Fractals in science. Springer, Berlin
Talamali M, Petäjä V, Vandembroucq D, Roux S (2011) Avalanches, precursors, and finite-size fluctuations in a mesoscopic model of amorphous plasticity. Phys Rev E 84:016115
Jack RL, Dunleavy AJ, Royall CP (2014) Information-theoretic measurements of coupling between structure and dynamics in glass formers. Phys Rev Lett 113:095703
Chen K, Manning ML, Yunker PJ, Ellenbroek WG, Zhang Z, Liu AJ, Yodh AG (2011) Measurement of correlations between low-frequency vibrational modes and particle rearrangements in quasi-two-dimensional colloidal glasses. Phys Rev Lett 107:108301
Schoenholz JJ, Liu AJ, Riggleman RA, Rottler J (2014) Understanding plastic deformation in thermal glasses from single-soft-spot dynamics. Phys Rev X 4:031014
Smessart A, Rottler J (2014) Structural relaxation in glassy polymers predicted by soft modes: a quantitative analysis. Soft Matter 10:8533–8541
Tanguy A, Mantisi B, Tsamados M (2010) Vibrational modes as a predictor for plasticity in a model glass. Europhys Lett 90:16004
Ghosh A, Chikkadi V, Shall P, Bonn D (2011) Connecting structural relaxation with the low frequency modes in a hard-sphere colloidal glass. Phys Rev Lett 107:108302
Schoenholz SS, Cubuk ED, Sussman DM, Kaxiras E, Liu AJ (2016) A structural approach to relaxation in glassy liquids. Nat Phys 12:469–471
Rycroft CH, Gibou F (2012) Simulations of a stretching bar using a plasticity model from the shear transformation zone theory. J Comput Phys 231:2155–2179
Rycroft CH, Bouchbinder E (2012) Fracture toughness of metallic glasses: annealing-induced embrittlement. Phys Rev Lett 109:194301
Bouchbinder E, Langer JS, Procaccia I (2007) Athermal shear-transformation-zone theory of amorphous plastic deformation. I. Basic principles. Phys Rev E 75:036107
Bouchbinder E, Langer JS, Procaccia I (2007) Athermal shear-transformation-zone theory of amorphous plastic deformation. II. Analysis of simulated amorphous silicon. Phys Rev E 75:036108
Bouchbinder E, Langer JS (2009) Nonequilibrium thermodynamics of driven amorphous materials. I. Internal degrees of freedom and volume deformation. Phys Rev E 80:031131
Bouchbinder E, Langer JS (2009) Nonequilibrium thermodynamics of driven amorphous materials. II. Effective-temperature theory. Phys Rev E 80:031132
Bouchbinder E, Langer JS (2009) Nonequilibrium thermodynamics of driven amorphous materials. III. Shear-transformation-zone plasticity. Phys Rev E 80:031133
Yang J, Duan J, Wang YJ, Jiang MQ (2020) Complexity of plastic instability in amorphous solids: Insights from spatiotemporal evolution of vibrational modes. Eur Phys J E 43:56
Raissi M, Perdikaris P, Karniadakis GE (2019) Physics-informed neural networks: a deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. J Comput Phys 378:686–707
Pedone A, Malavasi G, Menziani M, Cormack A, Segre U (2006) A new self-consistent empirical interatomic potential model for oxides, silicates, and silica-based glasses. J Phys Chem B 110:11780
Carré A, Horbach J, Ispas S, Kob W (2008) New fitting scheme to obtain effective potential from Car-Parrinello molecular-dynamics simulations. Europhys Lett 82:17001
Yu Y, Wang B, Wang M, Sant G, Bauchy M (2016) Revisiting silica with ReaxFF: towards improved predictions of glass structure and properties via reactive molecular dynamics. J Non-Cryst Solids 443:148–154
Tersoff J (1988) New empirical approach for the structure and energy of covalent systems. Phys Rev B 37:6991–7000
van Beest BWH, Kramer GJ, van Santen RA (1990) Force fields for silicas and alumina phosphates based on ab initio calculations. Phys Rev Lett 64:1955–1958
Zhang L, Sun S, Jahanshahi S (2001) Molecular dynamics simulations of silicate slags and slag-solid interfaces. J Non-Cryst Solids 282:24–29
Stillinger FH, Weber TA (1985) Computer simulation of local order in condensed phases of silicon. Phys Rev B 31:5262–5271
Rino J, Ebbsjo I, Kalia R, Nakano A, Vashishta P (1993) Structure of rings in vitreous \(\text{ SiO}_{2}\). Phys Rev B 47:3053–3062
Rountree C, Kalia R, Lidorikis E, Nakano A, Brutzel LV, Vashishta P (2002) Atomistic aspects of crack propagation in brittle materials: multimillion atom molecular dynamics simulations. Annu Rev Mater Res 32:377–400
Vashishta P, Kalia R, Nakano A (1999) Large-scale atomistic simulations of dynamic fracture. Comput Sci Eng 1:56–65
Kalia R, Nakano A, Vashishta P, Rountree C, Brutzel LV, Ogata S (2003) Multiresolution atomistic simulations of dynamic fracture in nanostructured ceramics and glasses. Int J Fract 121:71–79
Rountree C, Prades S, Bonamy D, Bouchaud E, Kalia R, Guillot C (2007) A unified study of crack propagation in amorphous silica: using experiments and simulations. J Alloys Compd 434–435:60–63
Rountree C, Bonamy D, Padres S, Dalmas D, Guillot C, Elisabeth B (2010) Fracture in glass via molecular dynamics simulations and atomic force microscopy experiments. Phys Chem Glass 51:127–132
Rountree C, Vandembroucq D, Tamali M, Bouchaud E, Roux S (2009) Plasticity-induced structural anisotropy of silica glass. Phys Rev Lett 102:195501
Vashishta P, Kalia R, Nakano A, Rino J (2007) Interaction potential for silicon carbide: a molecular dynamics study of elastic constants and vibrational density of states for crystalline and amorphous silicon carbide. J Appl Phys 101:103515
Vashishta P, Kalia R, Nakano A, Li W, Ebbsjö I (1997) Molecular dynamics methods and large-scale simulations of amorphous materials. In: Thorpe MF, Mitkova MI (eds) Amorphous insulators and semiconductors, vol 23. Kluwer Academic Publishers, pp 151–213
Méndez-Maldonado G, González-Melchor M, Alejandre J (2012) Phase equilibria and interfacial properties of two-dimensional Yukawa fluids. Condens Matter Phys 15:23002
Méndez-Maldonado G, González-Melchor M, Alejandre J (2012) Liquid-vapor phase diagram and cluster formation of two-dimensional ionic fluids. J Chem Phys 137:054711
Morley D, Wilson M (2019) Constructing bilayers with tuneable ring statistics and topologies. Mol Phys 117:3148–3157
Acknowledgements
This research has been funded by the Federal Ministry of Education and Research (BMBF) and the Ministry of Culture and Science of the German State of North Rhine-Westphalia (MKW) under the Excellence Strategy of the Federal Government and the Länder. Additionally, the first author thanks Michael Falk, Christopher Rycroft, Michael Shields, and Thomas Hardin for fruitful discussions.
Funding
Open Access funding enabled and organized by Projekt DEAL.
Author information
Authors and Affiliations
Corresponding author
Ethics declarations
Conflict of interest
On behalf of all authors, the corresponding author states that there is no conflict of interest.
Additional information
Publisher's Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Appendix A: Potentials for Network Glass Structures
Appendix A: Potentials for Network Glass Structures
In order to provide a realistic model of a network glass, such as bulk silica, various potentials were developed that not only include purely two-body but also combined two- and three-body interactions, cf. [19, 83, 288,289,290,291]. With regard to the numerical demonstrations in this paper, we will, however, concentrate on the two most essential potentials used in this context and present them in a nutshell.
1.1 The Best–Kramer–Santen Potential for Silica Glass
Van Best, Kramer, and van Santen presented a purely two-body interaction potential [292] to model the structure of silica glass efficiently. Starting from ab-initio calculations, they derived a bulk force field, providing an accurate description of the interactions of Si and O. The proposed pair potential has the form:
where the Coulomb interaction is described by the effective atomic charges \(q_{\text {Si}}=2.4\) and \(q_{\text {O}}=-1.2\). The parameters A, B, and C in the second and third terms form a Buckingham-type potential accounting for the covalent contribution. In order to reproduce the tetrahedral building units of silica, it is sufficient to allow for only two short-range interactions, i.e., the Si–O bond and the non-bonded O–O interaction. The parameters of the BKS potential are provided in Table 1.
We note that another popular force field used for the mechanical descriptions of silicate glasses, which also follows the formulation of purely two-body contributions, is represented by the Born-Meyer-Hudgins potential. However, the atomic charges are not described effectively in this case, but they are divided by a constant. We refer to Zhang et al. [293] for further reading.
1.2 The Vashishta Potential for Silica Glass
Silica glass consists of a binary mixture of two constituents, whereas \(\nicefrac {2}{3}\) of all particles are oxygen (O) atoms, and \(\nicefrac {1}{3}\) of all particles are silicon (Si) atoms. The potential proposed here includes two- and three-body contributions [19]. The three-body part of the potential constitutes an adaptation of the original three-body formulation that was initially developed by Stillinger and Weber [294] for pure silicon. However, the Vashishta potential describes the interaction of Si and O atoms and has extensively been applied in mechanical studies and validated by comparing numerical demonstrations and experimental results [295,296,297,298,299,300,301,302,303]. The two-body part \(\tilde{\mathcal {U}}_2\) includes the effects of steric repulsion, electronic polarisabilities, and Coulombic charge transfer and is written in the form of a pair-potential as follows:
The parameters \(H^{*}_{ij}\) and \(\eta _{ij}\) denote the repulsive strength and the exponent, while \(\alpha ^{*}_i\), \(Z_i\), and \(r_{4s}\) are the electronic polarizability, effective charge, and decay length, respectively.
Vashishta used the description of the three-body part \(\mathcal {U}_3\left( \varvec{r}_i,\varvec{r}_j,\varvec{r}_k\right)\), using the model proposed by Stillinger and Weber [294]. The three-body contribution is applied to favor bond angles in particular triplet combinations. Following the model by Vashishta [19], the potential of the triplets follows a quadratic form, written as:
where
if \(r_{ij},r_{ik} \le r_{3c}\) and zero otherwise. The remaining two terms \(\Lambda _{kji}\) and \(\Lambda _{ikj}\) belonging to the triplets \(\hat{h}^{(2)}\) and \(\hat{h}^{(3)}\) are defined equivalently. The cutoff distance for the three-body contributions \(r_{3c}\) should be fixed at such a value that takes into consideration the nearest neighbors of all possible atom types during the evaluation. All remaining atoms situated outside this radial distance are neglected for the three-body term. Exemplary, the parameter \(\vartheta _{jik,0}\) denotes the equilibrium bond angle enclosed by the vectors \(\varvec{r}_{ij}\) and \(\varvec{r}_{ik}\). All possible combinations of the B- and the \(\vartheta\)-constants represent the strength of the three-body part and the angle between the two respective triplet vectors. Since silica glass consists of two constituents, each of the three atomic coordinates \(\varvec{r}_i\), \(\varvec{r}_j\), and \(\varvec{r}_k\), defining the three-body contribution \(\mathcal {U}_3\left( \varvec{r}_i,\varvec{r}_j,\varvec{r}_k\right)\), may either direct to an Si or an O atom. Following this logic, there are eight (\(2^3\)) possible triplet combinations, i.e., Si–Si–Si, Si–Si–O, Si–O–O, O–Si–Si, O–O–Si, O–O-O, O–Si–O, and Si–O–Si. However, it is sufficient to only consider the last two triplets [19] with the respective equilibrium bond angles \(\vartheta _{\text{Si}{-}\text{O}{-}\text{Si},0}\) and the \(\vartheta _{\text{O}{-}\text{Si}{-}\text{O},0}\), as indicated in Fig. 49. This results in a significant saving of computational effort.
The two relevant triplets of the three-body contribution of the Vashishta potential: a the O–Si–O triplet; b the Si–O–Si triplet
The corresponding parameters of the potential are presented in Table 2.
1.3 Yukawa Potential for 2D Network Materials
We present this potential as the counterpart to the Lennard-Jones potential for two-dimensional model glasses since the objective is to generate a physically meaningful network structure variant of a model 2D glass in this paper. As we show in Sect. 8, this potential captures all information regarding the network topology, following the definitions of Zachariasen’s random network theory [101]. Thereby, one may investigate the in-plane mechanical behavior of silica-like networks. In doing so, one implements a Yukawa-type potential function that allows for the network formation of silica in two dimensions [115, 304, 305]. The functional form of the potential reads:
where \(r_{ij}\), \(a_{ij}\), \(q_{ij}\) and \(\kappa\) denote the interatomic distance, the repulsive strength, the ionic charge, and the screening parameter, respectively. For more details about the parameters, the readers are referred to Roy et al. [115]. This potential function successfully reproduces the topological measures of a two-dimensional representation of silica glass, represented by the short-range order and the ring statistics, as discussed in Section 8. Both topological identification measures align closely with experimental studies [116, 117, 306]. We furthermore note that this potential is dimensionless by definition. Therefore, all quantities, except length scales, are expressed in dimensionless units.
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
About this article
Cite this article
Bamer, F., Ebrahem, F., Markert, B. et al. Molecular Mechanics of Disordered Solids. Arch Computat Methods Eng 30, 2105–2180 (2023). https://doi.org/10.1007/s11831-022-09861-1
Received:
Accepted:
Published:
Issue Date:
DOI: https://doi.org/10.1007/s11831-022-09861-1