# Vision-Augmented Molecular Dynamics Simulation of Nanoindentation

Authors: Rajab Al-Sayegh; Charalampos Makatsoris
Journal: Journal of Nanomaterials
 (2015)

---

## Abstract

We present a user-friendly vision-augmented technique to carry out atomic simulation using hand gestures. The system is novel in its concept as it enables the user to directly manipulate the atomic structures on the screen, in 3D space using hand gestures, allowing the exploration and visualisation of molecular interactions at different relative conformations. The hand gestures are used to pick and place atoms on the screen allowing thereby the ease of carrying out molecular dynamics simulation in a more efficient way. The end result is that users with limited expertise in developing molecular structures can now do so easily and intuitively by the use of body gestures to interact with the simulator to study the system in question. The proposed system was tested by simulating the crystal anisotropy of crystalline silicon during nanoindentation. A long-range (Screened bond order) Tersoff potential energy function was used during the simulation which revealed the value of hardness and elastic modulus being similar to what has been found previously from the experiments. We anticipate that our proposed system will open up new horizons to the current methods on how an MD simulation is designed and executed.


---

## Body

## 1. Introduction

Computational simulation approaches, such as molecular dynamics (MD), Monte Carlo, or coarse grained MD, are widely employed by chemists and materials scientists to predict and understand the behaviour of materials. These approaches have a high computational cost as they tend to get confined within local regions of the energy landscape spending many CPU during a simulation run. This research is focused on the development of a virtual environment which offers an entirely new way to explore complex energy landscapes by allowing the user to directly interact with the simulated molecular system [1]. Our approach entails the use of computer vision and specifically gesture recognition to enable the user to interact with the simulated system exploiting the numerous degrees of freedom available with hand gestures in 3D space, breaking away from the constraints of traditional methods of interaction involving a high level language based inputs which became the motivation for this paper. This approach progresses the current state of the art via a step change in understanding of the complex interactions as well as allowing the rapid design of novel molecular structures using an intuitive way of interaction with a simulator by means of hand gestures. We demonstrate our approach and describe our vision-augmented molecular dynamics simulation (VAMDS) system on the nanoindentation of a silicon crystal silicon substrate. The paper is organised as follows. Section 2 briefly reviews earlier efforts in design and simulation systems driven by haptic devices. Section 3 presents our proposed VAMDS methodology. In Section 4, the nanoindentation study using the VAMDS system is described. Finally in Section 5, we conclude by presenting the simulation results obtained from the system exploring crystal anisotropy of silicon during its nanoindentation.


## 2. Literature Review

The work in this paper was especially inspired from the pioneering work of Mistry et al. [2] at MIT, in an attempt to help the disabled people produced a wearable gesture interface using a small camera and a projector. The purpose of this device was to interact with real world objects. The system recognised freehand gestures like hand waving, multitouch systems like shapes with both hands, and iconic gestures like in air drawing. They also developed a LCD/CRT display which provides different contents for different viewers wearing 3D shutter glasses simultaneously. The system can switch between transparent and opaque as when one viewer is transparent, other is opaque at frame rate of 120 Hz. This is a good technique to add more viewers at higher frame rate but it requires wearing of special glasses. Mistry et al. [3] also developed a hand free input interface for controlling and commanding a robot. The glaze and blink input directs the robot to move an object from one location to other. The “Blinkbot” replaces traditional interactive devices such as mouse, keyboard, and speech based input systems which are not very good at directional input in space. With the Blinkbot wearable controller, a user provides input by a gaze to select an object and actions are triggered by a blink. The action is to move the object to a desired location by the robot. This is a novel technique which provides highly interactive sense of touch and can be used to help disabled people. They presented a technique to establish the link between computer aided design (CAD) tools and handmade pencil and paper based designs. Taylor and Johnson [4] proposed a haptic interface in which parameters of the computer simulations are controlled by physical means like robots. The demonstration of tangible simulation is performed by three body gravitational simulation. Ricci et al. [5] used haptic devices for computer aided drug design. Haptic driven simulators provided the user an opportunity to drive and control the simulation, increased interactivity, and combines human knowledge with the computational power. Dreher et al. [6] used a framework composed of selector, powerful multicore processors and visualiser to successfully connect 3D virtual reality environment with the molecular dynamics simulator/Gromacs where forces applied to haptic devices were used to select atoms for steering large scale MD simulations. Thus, a number of research efforts are evident concerning automation of the various systems. To this end, we progress the state of the art by utilizing a haptic device to carry out the atomistic simulations.The use of haptic device in molecular dynamics system is a novel proposition and hence is the major focus of this research. Molecular dynamics simulation with haptic devices provides more flexibility of picking and choosing various atoms and offers much flexibility to the inexperienced users by providing an opportunity to deal with the atoms in a very interactive way. Beside many other tasks, the use of haptic devices in conjunction with molecular dynamics will assist in the following:(1)
Make the design of nanostructured materials and processes via atomistic simulation more accessible to inexperienced users.(2)
Enable the design of novel materials, using an intuitive approach.(3)
Introduce a new methodology that will allow the design and execution of MD simulations in a more interactive manner with the expert user in the loop.(4)
Provide a new method that allows for the in depth understanding of atomic scale processes with novel visualisation techniques.Inspired from this, the previous study of the authors [1] was focused on developing a virtual environment for use with the Framework Rigidity Optimised Dynamic Algorithm (FRODA) [8] where the system made use of hand gestures to manipulate various chemical structures on the computer screen to control the simulation. In this paper, we built on this development to test a case study by performing nanoindentation of a silicon substrate by exploring its anisotropy by using a long-range Tersoff potential energy function [7] on a high performance computing (HPC) system.


## 3. Methodology

This section presents the working scheme of the system. This work made use of the open source codes only. Accordingly, we used LAMMPS [9] to perform the MD simulations while VMD [10], OVITO [11], and JMOL [12] were used to visualize and analyse the atomistic simulation data. The commands captured by the haptic device concerning modified atomic positions are sent to the visualization software, for example, Jmol and OVITO, using gestures recognition. Algorithms were written which in turn help to obtain a description of a newer geometry with modified positions or orientations of the atoms. The modified file is saved in a location which is then fed to the molecular dynamics simulator like LAMMPS for performing the simulation with revised inputs. A schematic diagram of this scheme is shown in Figure 1 while a more functional scheme of this method along with the representation of simulation results at various steps is shown in Figure 2.Figure 1
Vision-augmented molecular dynamics simulation.Figure 2
Integration between computer vision and molecular dynamics simulation.Human gestures recognition is an important task where haptic devices are used. The Kinect device (Xbox technology) introduced by Microsoft in the year 2010 that has the features of gestures recognition, facial recognition, and voice recognition was adapted as a hardware in this work. In addition, Prime Sense Company developed and released an open source motion tracking software called NITE middleware that analyse the data from the hardware [13]. NITE provides an application programming interface (API) which was used in this work to control the hand gestures. NITE algorithm also provides additional features such as scene segmentation, hand point detection and tracking, and full body tracking. As shown in Figures 1 and 2, there are four distinct components of this system: the haptic device (Microsoft Kinect), JMOL (to modify the positions of the atoms), LAMMPS (to run the MD simulation), and OVITO (to prepare input file for LAMMPS and to simultaneously visualize the LAMMPS output). Since the system is based on open source codes, it provides full flexibility to modify the system in any desirable way. Gesture application (App) and Jmol visualiser were connectedvia a TCP/IP protocol (see dotted lines in Figure 1). The Gesture App detects human gestures, converts them into equivalent Jmol commands, and then sends this to Jmol server using TCP/IP protocols; for example, a push gesture for left hand was converted to spin on/off command of the Jmol. Jmol server receives this spin on/off command and performs the spin operation on the atomic structure of the material. Subsequently, OVITO was used to convert XYZ file format (output of Jmol) to the LAMMPS data file format which was sent to HCP server.The entire process was automated so that the HPC server can call the LAMMPS engine to generate the simulation output in the form of a trajectory at the desired time steps. The output files obtained from LAMMPS were programmed to save in the local computer from the HPC server and the process continues until the whole directory containing the LAMMPS output is copied to the local computer. This cycle continues while in the meantime, OVITO interactively make the use visualise each file copied from the HPC server to the local computer. Briefly, the vision-augmented molecular dynamics simulation developed in this work is comprised of the key steps shown in Figure3.Figure 3
Flowchart of the VAMD developed in this work.


### 3.1. Gestures to Jmol Command Conversion

Since we intend to manipulate the position of atoms merely using gestures, we need to convert the gestures into appropriate JMOL commands. Some of these conversions are tabulated in Table1. In order to send the JMOL command to the JMOL Visualizer, we use a Wireshark tool to capture the data on TCP/IP port to identify the exact structure of the command. As an example, in order to send the command zoom in, we need to execute movement Swipe left (left hand) and Swipe right (right hand) in front of haptic device. Gesture application will convert this to the Jmol command format which is “zoom in” and send it to Jmol server via TCP/IP port to do the operation “zoom in” to the molecular structure visualised on Jmol.Table 1
Gestures to JMOL command mapping.
S/number
Hand gestures
Jmol commands
1
Push (right hand)
Left mouse button click
2
Swipe right (right hand)
Right button mouse click
3
Swipe left (right hand)
Set picking drag molecule
4
Swipe up (right hand)
Set picking drag atom
5
Swipe down (right hand)
Set picking drag bond
6
Swipe right (left hand) and Swipe left (right hand)
Zoom out
7
Swipe left (left hand) and Swipe right (right hand)
Zoom in
8
Swipe up (left hand)
Set picking delete atom
9
Swipe down (left hand)
Set picking delete bond
10
Push (left hand)
Spin on/off
11
Circle (right hand)
Write the output from Jmol inXYZ format and send this to OVITOWe prepare this string for each command and send it over TCP/IP to the JMOL Visualizer. Once the string is received by the JMOL server, it executes a particular command and the respective action is taken accordingly. As an example, in the above case, the loaded molecule gets zoomed in. For other commands, we just replace the zoom in field with the corresponding JMOL command.


## 3.1. Gestures to Jmol Command Conversion

Since we intend to manipulate the position of atoms merely using gestures, we need to convert the gestures into appropriate JMOL commands. Some of these conversions are tabulated in Table1. In order to send the JMOL command to the JMOL Visualizer, we use a Wireshark tool to capture the data on TCP/IP port to identify the exact structure of the command. As an example, in order to send the command zoom in, we need to execute movement Swipe left (left hand) and Swipe right (right hand) in front of haptic device. Gesture application will convert this to the Jmol command format which is “zoom in” and send it to Jmol server via TCP/IP port to do the operation “zoom in” to the molecular structure visualised on Jmol.Table 1
Gestures to JMOL command mapping.
S/number
Hand gestures
Jmol commands
1
Push (right hand)
Left mouse button click
2
Swipe right (right hand)
Right button mouse click
3
Swipe left (right hand)
Set picking drag molecule
4
Swipe up (right hand)
Set picking drag atom
5
Swipe down (right hand)
Set picking drag bond
6
Swipe right (left hand) and Swipe left (right hand)
Zoom out
7
Swipe left (left hand) and Swipe right (right hand)
Zoom in
8
Swipe up (left hand)
Set picking delete atom
9
Swipe down (left hand)
Set picking delete bond
10
Push (left hand)
Spin on/off
11
Circle (right hand)
Write the output from Jmol inXYZ format and send this to OVITOWe prepare this string for each command and send it over TCP/IP to the JMOL Visualizer. Once the string is received by the JMOL server, it executes a particular command and the respective action is taken accordingly. As an example, in the above case, the loaded molecule gets zoomed in. For other commands, we just replace the zoom in field with the corresponding JMOL command.


## 4. MD Simulation of the Nanoindentation of Crystalline Silicon with VAMDS

### 4.1. MD Simulation Methodology

In this section we demonstrate the methodology adopted for simulating nanoindentation of silicon using VAMDS. The aim of the study was to explore the crystal anisotropy of silicon in such a way that the orientation of the crystal can be switched merely by the hand gestures. To this end, the simulation model was first developed by following the standard algorithm, boundary conditions, and ensemble in the same way as has been previously done in other research papers [14–16]. A schematic simulation model used to carry out the simulation is shown in Figure 4. In Figure 4, the atoms in the indenter were kept fixed (diamond indenter was assumed to be an infinite rigid body structure) and, in the three test cases, the crystal orientation of silicon was varied. Goel et al. [14] have provided a glance of other research works where a higher velocity (up to 500 m/s) has been frequently used to overcome the limitations of MD simulation. To achieve more realistic results, we made use of an indentation velocity of 20 m/s as has been recently used on HPC systems [17]. Also, we made use of long-range Tersoff potential energy function [7] for describing silicon and carbon as it provides improved description of silicon.Figure 4
Schematic of the MD simulation.The details of the parameters used to develop the MD simulation model are shown in Table2. In the MD simulation model, to avoid any artificial effect of the temperature (due to thermal fluctuations and thermal vibrations), a low temperature of 10 K was used to equilibrate the sample and perform the nanoindentation. During the MD simulation, the atomic stress tensor (http://lammps.sandia.gov/doc/compute_stress_atom.html) was calculated by considering an elemental atomic volume (1 nm × 1 nm × 1 nm) of silicon in the deformation zone right underneath the indenter.Table 2
Details used for the development of MD simulation model.
Dimensions of the silicon workpiece
Crystal orientation of silicon
Number of atoms
16.296 nm × 10.864 nm × 16.296 nm
(010)
144000
15.9667 nm × 10.7548 nm × 15.9667 nm
(110)
137088
15.9667 nm × 11.29 nm × 16.132 nm
(111)
145152
Equilibrium parameters
Silicon: 5.432 Å and Diamond: 3.5656 Å
Indenter and specimen surface distance
1 nm
Depth of indentation
2 nm
Speed of indentation and retraction
20 m/s = 0.02 nm/ps
Diameter of the indenter
4 nm
Boundary conditions
Periodic inX and Z direction
Potential function used
Screening bond order Tersoff potential function [7]In addition to stress calculation, further analysis was carried out by estimating the contact pressure (pm) underneath the indenter. Contact pressure (pm) is defined as the ratio of the instantaneous normal force (F) on the indenter and the projected area (A) which is calculated as (2π×a×h) for a spherical indenter where a=h(2R-h) is the contact radius of the spherical indenter, R is the radius of the indenter, and h is the instantaneous displacement of the indenter. The average value of the contact pressure (after it attains saturation value) is also referred to as nanoindentation hardness (H) of the material.A popular way of calculation of the elastic modulus (Es) and hardness (H) of the material from the load-displacement (P-h) plot was proposed by Oliver and Pharr [18] during the last decade of 20th century. Their method relies on calculation of the projected contact area by drawing an imaginary line, typically following the power law, on the top 1/3rd part along the unloading curve of the P-h plot. The slope of this curve (S) thus enables to obtain the reduced elastic modulus of the material which can then be used to obtain Es (elastic modulus) of the specimen using the following equations:(1)Er=1βπ2SA,1-υs2Es=1Er-1-υi2Ei.In (1), Er is the reduced elastic modulus (N/m2), S is the slope of the top 1/3rd part of the unloading curve which were obtained as 729 N/m (for (010) oriented Si), 820 N/m (for (110) oriented Si), and 914 N/m (for (111) oriented Si) from the simulation results, β=1 is a constant for spherical indenter, and A is the projected area (m2) which varies with the depth of indentation. For accuracy of computation of Es, other parameters such as υs=0.3459 (Poisson’s ratio of silicon), υi=0.103 (Poisson’s ratio of Diamond), and Ei=1056.48 GPa (Young’s modulus of diamond) were obtained directly from the potential function.


### 4.2. MD Simulation Results

#### 4.2.1. Calibration of the Gesture Application with the MD Simulation Model

Figure5 shows the progressive physical movement of the human hand sensed by the gesture application (left column) and transformation of this physical hand movement into the corresponding displacement of the indenter into the substrate (right hand column). The physical distances of the hand movements were captured through a measurement tape and were compared with the movement of the indenter in the MD simulation model so as to see the factor by which the hand movement transfers into the equivalent displacement of the indenter. The relation between these two movements was plotted and is shown in Figure 6. Knowing this relation was a preliminary step for calibrating the newly proposed VAMDS system.Figure 5
Comparison of physical movement with the movement in the MD simulation model (a) at the beginning of the simulation, (b) when the indenter establishes contact with the substrate and (c) final indentation depth was achieved.
(a)
(b)
(c)Figure 6
Relationship between physical movement of hand gesture and indenter in the MD model.


#### 4.2.2. HPPT of Silicon and Extraction of the Mechanical Properties fromP-h Plot

The ductile behaviour of brittle materials is often attributed to an event recognized as high-pressure phase transformation (HPPT). HPPT enforces brittle materials to become ductile for the short duration of loading which upon unloading eventually transforms to amorphous phase of the brittle material. HPPT of silicon during its nanoindentation has been explained to be driven from the deviatoric stress rather than temperature [14] in the deformation zone. Much of the literature on silicon nanoindentation have reported HPPT to be the primary mechanism governing the plasticity of silicon that causes brittle-ductile transition. There is however a noticeable exception of Mylvaganam and Zhang [19] observed nanotwinning (associated with Si-I to bct-5 phase transformation) in silicon along the 110 direction that stops at Shockley partial dislocation. It was therefore relevant to assert the occurrence of HPPT in this work for which an indicator known as coordination number was used for its ease of implementation.A snapshot from the MD simulation showing the change in coordination number during peak loading of silicon on various orientations is shown in Figure7. The coordination number of 4 (Figure 7) is representative of bulk silicon while coordination number of 1 and 2 on the free surface of silicon workpiece represents the dangling bonds on the Si-terminated nascent surface. During the nanoindentation, the coordination number in the indentation zone was observed to change from 4 to 6 with a corresponding change in the number of nearest neighbour atoms. Literature suggests that this change corresponds to an allotropic transformation of silicon from the diamond cubic (alpha-silicon) to the body centred tetragonal (beta-silicon) structure and is consistent with the previously reported results of Si-I to Si-II during loading. Most of the metastable phases of silicon were observed to form along the 110 direction, which is the direction of slip plane in silicon [20] as is shown in Figure 7 while the (111) planes act as both slip planes and cleavage planes.Figure 7
Change in coordination number at peak indentation depth of 2 nm on three orientations.Figure8 is a P-h plot obtained from the displacement controlled nanoindentation simulation where the penetration depth of the diamond tool was kept fixed at 2 nm. An interesting observation concerning the P-h plot comes from the work of Jang et al. [21] where a good correlation between unloading curve and the phase of the material formed postindentation of silicon was observed. Based on the micro-Raman spectra, the authors proposed that the unloading discontinuity, often called as “pop-out,” corresponds to the formation of metastable Si-XII/Si-III crystalline phases, while the hysteresis, called the “elbow” is associated with the formation of a-Si. In Figure 8, an interesting feature observed from the P-h plot was that the unloading curve for the (110) and the (111) substrate followed power law whereas the unloading curve on the (010) substrate showed an elbow. This characteristic is a testimony to the formation of a-Si upon the retraction of the indenter on the (010) orientation.Figure 8
P
-
h plots of the nanoindentation for all three crystal orientations of silicon.Another interesting feature observed from theP-h plot obtained from the MD simulation was the residual indentation depth (hf) that can also be used to characterize the recovery of the indented surface upon retraction of the indenter. In this case, hf was quantified from the P-h plot (Figure 8) for all the three simulated orientations and was found inevitably different in all the cases; that is, hf (010) was minimum while hf for the (111) orientation was maximum whereas hf (110) was intermittent. This suggests that at atomic scale, the (010) orientation will show higher extent of elastic recovery than the (111) orientation of silicon. This result is somewhat consistent with the earlier work of Shibata et al. [22] where a Schmidt-type slip orientation factor was proposed and the 11-0 direction either on the (100) or on the (111) planes was recognized as more preferred combination to plastically deform the silicon. Lastly, beside a monotonic increase in the loading curve, a strong cohesion between the indenter and the substrate during retraction was also observed.HPPT is highly direction sensitive [23, 24]; that is, pressure or hydrostatic stress required for transformation on the (111) orientation is lower than that required to induce the transformation on the (100) orientation [20]. Literature suggests that the transformation pressure required to cause metallisation of silicon varies between 9 and 18 GPa [25]. It may however be noted that the hydrostatic pressure differs from deviatoric stress [26]. During nanoindentation, silicon undergoes both shear and compression making it different from the case of simple compressive pressure.As it can be seen from Table3, the hydrostatic stress (pressure) [27] in the indentation zone was observed to be maximum on the (110) orientation followed by (010) orientation and it was minimum for the (111) orientation. This correlation was however not the same for the von Mises stress and Tresca stress in the indentation zone.Table 3
Computational results obtained from the nanoindentation simulation.
Silicon
(010) orientation
(110) orientation
(111) orientation
Average hardness (GPa) at indent depth of 2 nm
16.99
18.22
18.48
Elastic modulus (GPa) at indent depth of 2 nm
129
147
167
h
f
/
h
max
⁡ (extent of elastic recovery)
1.17/2=0.585
1.36/2=0.68
1.47/2=0.735
Peak maximum von Mises stress in the indent zone (GPa)
11.07
22.75
16.72
Peak maximum Tresca stress in the indent zone (GPa)
5.91
12.11
9.14
Peak maximum octahedral shear stress in the indent zone (GPa)
5.21
10.72
7.88
Peak hydrostatic stress in the indent zone (GPa)
20.11
21.98
16.62
Peak temperature in the indentation zone (K)
300
419
438
Features of theP-h curve during unloading
Elbow formation corresponding to a-Si formation
No feature
No featureFor the two kinds of deviatoric stresses, the (110) orientation showed maximum values whereas the (010) orientation showed least value. The peak average temperature in the cutting zone was rather high, especially considering the fact that nanoindentation process was carried out at 10 K. At its peak, the (111) orientation showed the maximum temperature of about 438 K followed by the (110) orientation at 419 K whereas the (010) orientation showed the least temperature of 300 K. Finally, the elastic modulus was obtained fromP-h plot using Oliver and Pharr method and accordingly the variation in E and hardness on all the three crystal surfaces with respect to change in the indentation depth is plotted and shown in Figure 9.Figure 9
Variation in the (a) elastic modulus and (b) nanoindentation hardness of silicon with respect to the indentation depth.
(a)
(b)It is interesting to note from Figure9(a) that the elastic modulus at finite indentation depth of few nanometres decrease sharply until attaining a saturated value as the indentation depth increases. At an indentation depth of 2 nm, the elastic modulus of the three simulated orientations was obtained as 129 GPa for the (010) orientation, 147 GPa for the (110) orientation, and 167 GPa for the (111) orientation. These values are well comparable with the experimental values of the elastic modulus of silicon [28, 29] and resembles similar to what has been obtained from MD [30].It has been seen how a new approach is helpful in carrying out atomic simulations. Future work of this research will realise multiscale approach by integrating computer vision (e.g., gesture recognition and artificial intelligence), molecular dynamics (e.g., molecular modelling and material modelling), and intuitive understating of atomic scale. By doing this, we will be able to interact with moleculesvia our senses.


## 4.1. MD Simulation Methodology

In this section we demonstrate the methodology adopted for simulating nanoindentation of silicon using VAMDS. The aim of the study was to explore the crystal anisotropy of silicon in such a way that the orientation of the crystal can be switched merely by the hand gestures. To this end, the simulation model was first developed by following the standard algorithm, boundary conditions, and ensemble in the same way as has been previously done in other research papers [14–16]. A schematic simulation model used to carry out the simulation is shown in Figure 4. In Figure 4, the atoms in the indenter were kept fixed (diamond indenter was assumed to be an infinite rigid body structure) and, in the three test cases, the crystal orientation of silicon was varied. Goel et al. [14] have provided a glance of other research works where a higher velocity (up to 500 m/s) has been frequently used to overcome the limitations of MD simulation. To achieve more realistic results, we made use of an indentation velocity of 20 m/s as has been recently used on HPC systems [17]. Also, we made use of long-range Tersoff potential energy function [7] for describing silicon and carbon as it provides improved description of silicon.Figure 4
Schematic of the MD simulation.The details of the parameters used to develop the MD simulation model are shown in Table2. In the MD simulation model, to avoid any artificial effect of the temperature (due to thermal fluctuations and thermal vibrations), a low temperature of 10 K was used to equilibrate the sample and perform the nanoindentation. During the MD simulation, the atomic stress tensor (http://lammps.sandia.gov/doc/compute_stress_atom.html) was calculated by considering an elemental atomic volume (1 nm × 1 nm × 1 nm) of silicon in the deformation zone right underneath the indenter.Table 2
Details used for the development of MD simulation model.
Dimensions of the silicon workpiece
Crystal orientation of silicon
Number of atoms
16.296 nm × 10.864 nm × 16.296 nm
(010)
144000
15.9667 nm × 10.7548 nm × 15.9667 nm
(110)
137088
15.9667 nm × 11.29 nm × 16.132 nm
(111)
145152
Equilibrium parameters
Silicon: 5.432 Å and Diamond: 3.5656 Å
Indenter and specimen surface distance
1 nm
Depth of indentation
2 nm
Speed of indentation and retraction
20 m/s = 0.02 nm/ps
Diameter of the indenter
4 nm
Boundary conditions
Periodic inX and Z direction
Potential function used
Screening bond order Tersoff potential function [7]In addition to stress calculation, further analysis was carried out by estimating the contact pressure (pm) underneath the indenter. Contact pressure (pm) is defined as the ratio of the instantaneous normal force (F) on the indenter and the projected area (A) which is calculated as (2π×a×h) for a spherical indenter where a=h(2R-h) is the contact radius of the spherical indenter, R is the radius of the indenter, and h is the instantaneous displacement of the indenter. The average value of the contact pressure (after it attains saturation value) is also referred to as nanoindentation hardness (H) of the material.A popular way of calculation of the elastic modulus (Es) and hardness (H) of the material from the load-displacement (P-h) plot was proposed by Oliver and Pharr [18] during the last decade of 20th century. Their method relies on calculation of the projected contact area by drawing an imaginary line, typically following the power law, on the top 1/3rd part along the unloading curve of the P-h plot. The slope of this curve (S) thus enables to obtain the reduced elastic modulus of the material which can then be used to obtain Es (elastic modulus) of the specimen using the following equations:(1)Er=1βπ2SA,1-υs2Es=1Er-1-υi2Ei.In (1), Er is the reduced elastic modulus (N/m2), S is the slope of the top 1/3rd part of the unloading curve which were obtained as 729 N/m (for (010) oriented Si), 820 N/m (for (110) oriented Si), and 914 N/m (for (111) oriented Si) from the simulation results, β=1 is a constant for spherical indenter, and A is the projected area (m2) which varies with the depth of indentation. For accuracy of computation of Es, other parameters such as υs=0.3459 (Poisson’s ratio of silicon), υi=0.103 (Poisson’s ratio of Diamond), and Ei=1056.48 GPa (Young’s modulus of diamond) were obtained directly from the potential function.


## 4.2. MD Simulation Results

### 4.2.1. Calibration of the Gesture Application with the MD Simulation Model

Figure5 shows the progressive physical movement of the human hand sensed by the gesture application (left column) and transformation of this physical hand movement into the corresponding displacement of the indenter into the substrate (right hand column). The physical distances of the hand movements were captured through a measurement tape and were compared with the movement of the indenter in the MD simulation model so as to see the factor by which the hand movement transfers into the equivalent displacement of the indenter. The relation between these two movements was plotted and is shown in Figure 6. Knowing this relation was a preliminary step for calibrating the newly proposed VAMDS system.Figure 5
Comparison of physical movement with the movement in the MD simulation model (a) at the beginning of the simulation, (b) when the indenter establishes contact with the substrate and (c) final indentation depth was achieved.
(a)
(b)
(c)Figure 6
Relationship between physical movement of hand gesture and indenter in the MD model.


### 4.2.2. HPPT of Silicon and Extraction of the Mechanical Properties fromP-h Plot

The ductile behaviour of brittle materials is often attributed to an event recognized as high-pressure phase transformation (HPPT). HPPT enforces brittle materials to become ductile for the short duration of loading which upon unloading eventually transforms to amorphous phase of the brittle material. HPPT of silicon during its nanoindentation has been explained to be driven from the deviatoric stress rather than temperature [14] in the deformation zone. Much of the literature on silicon nanoindentation have reported HPPT to be the primary mechanism governing the plasticity of silicon that causes brittle-ductile transition. There is however a noticeable exception of Mylvaganam and Zhang [19] observed nanotwinning (associated with Si-I to bct-5 phase transformation) in silicon along the 110 direction that stops at Shockley partial dislocation. It was therefore relevant to assert the occurrence of HPPT in this work for which an indicator known as coordination number was used for its ease of implementation.A snapshot from the MD simulation showing the change in coordination number during peak loading of silicon on various orientations is shown in Figure7. The coordination number of 4 (Figure 7) is representative of bulk silicon while coordination number of 1 and 2 on the free surface of silicon workpiece represents the dangling bonds on the Si-terminated nascent surface. During the nanoindentation, the coordination number in the indentation zone was observed to change from 4 to 6 with a corresponding change in the number of nearest neighbour atoms. Literature suggests that this change corresponds to an allotropic transformation of silicon from the diamond cubic (alpha-silicon) to the body centred tetragonal (beta-silicon) structure and is consistent with the previously reported results of Si-I to Si-II during loading. Most of the metastable phases of silicon were observed to form along the 110 direction, which is the direction of slip plane in silicon [20] as is shown in Figure 7 while the (111) planes act as both slip planes and cleavage planes.Figure 7
Change in coordination number at peak indentation depth of 2 nm on three orientations.Figure8 is a P-h plot obtained from the displacement controlled nanoindentation simulation where the penetration depth of the diamond tool was kept fixed at 2 nm. An interesting observation concerning the P-h plot comes from the work of Jang et al. [21] where a good correlation between unloading curve and the phase of the material formed postindentation of silicon was observed. Based on the micro-Raman spectra, the authors proposed that the unloading discontinuity, often called as “pop-out,” corresponds to the formation of metastable Si-XII/Si-III crystalline phases, while the hysteresis, called the “elbow” is associated with the formation of a-Si. In Figure 8, an interesting feature observed from the P-h plot was that the unloading curve for the (110) and the (111) substrate followed power law whereas the unloading curve on the (010) substrate showed an elbow. This characteristic is a testimony to the formation of a-Si upon the retraction of the indenter on the (010) orientation.Figure 8
P
-
h plots of the nanoindentation for all three crystal orientations of silicon.Another interesting feature observed from theP-h plot obtained from the MD simulation was the residual indentation depth (hf) that can also be used to characterize the recovery of the indented surface upon retraction of the indenter. In this case, hf was quantified from the P-h plot (Figure 8) for all the three simulated orientations and was found inevitably different in all the cases; that is, hf (010) was minimum while hf for the (111) orientation was maximum whereas hf (110) was intermittent. This suggests that at atomic scale, the (010) orientation will show higher extent of elastic recovery than the (111) orientation of silicon. This result is somewhat consistent with the earlier work of Shibata et al. [22] where a Schmidt-type slip orientation factor was proposed and the 11-0 direction either on the (100) or on the (111) planes was recognized as more preferred combination to plastically deform the silicon. Lastly, beside a monotonic increase in the loading curve, a strong cohesion between the indenter and the substrate during retraction was also observed.HPPT is highly direction sensitive [23, 24]; that is, pressure or hydrostatic stress required for transformation on the (111) orientation is lower than that required to induce the transformation on the (100) orientation [20]. Literature suggests that the transformation pressure required to cause metallisation of silicon varies between 9 and 18 GPa [25]. It may however be noted that the hydrostatic pressure differs from deviatoric stress [26]. During nanoindentation, silicon undergoes both shear and compression making it different from the case of simple compressive pressure.As it can be seen from Table3, the hydrostatic stress (pressure) [27] in the indentation zone was observed to be maximum on the (110) orientation followed by (010) orientation and it was minimum for the (111) orientation. This correlation was however not the same for the von Mises stress and Tresca stress in the indentation zone.Table 3
Computational results obtained from the nanoindentation simulation.
Silicon
(010) orientation
(110) orientation
(111) orientation
Average hardness (GPa) at indent depth of 2 nm
16.99
18.22
18.48
Elastic modulus (GPa) at indent depth of 2 nm
129
147
167
h
f
/
h
max
⁡ (extent of elastic recovery)
1.17/2=0.585
1.36/2=0.68
1.47/2=0.735
Peak maximum von Mises stress in the indent zone (GPa)
11.07
22.75
16.72
Peak maximum Tresca stress in the indent zone (GPa)
5.91
12.11
9.14
Peak maximum octahedral shear stress in the indent zone (GPa)
5.21
10.72
7.88
Peak hydrostatic stress in the indent zone (GPa)
20.11
21.98
16.62
Peak temperature in the indentation zone (K)
300
419
438
Features of theP-h curve during unloading
Elbow formation corresponding to a-Si formation
No feature
No featureFor the two kinds of deviatoric stresses, the (110) orientation showed maximum values whereas the (010) orientation showed least value. The peak average temperature in the cutting zone was rather high, especially considering the fact that nanoindentation process was carried out at 10 K. At its peak, the (111) orientation showed the maximum temperature of about 438 K followed by the (110) orientation at 419 K whereas the (010) orientation showed the least temperature of 300 K. Finally, the elastic modulus was obtained fromP-h plot using Oliver and Pharr method and accordingly the variation in E and hardness on all the three crystal surfaces with respect to change in the indentation depth is plotted and shown in Figure 9.Figure 9
Variation in the (a) elastic modulus and (b) nanoindentation hardness of silicon with respect to the indentation depth.
(a)
(b)It is interesting to note from Figure9(a) that the elastic modulus at finite indentation depth of few nanometres decrease sharply until attaining a saturated value as the indentation depth increases. At an indentation depth of 2 nm, the elastic modulus of the three simulated orientations was obtained as 129 GPa for the (010) orientation, 147 GPa for the (110) orientation, and 167 GPa for the (111) orientation. These values are well comparable with the experimental values of the elastic modulus of silicon [28, 29] and resembles similar to what has been obtained from MD [30].It has been seen how a new approach is helpful in carrying out atomic simulations. Future work of this research will realise multiscale approach by integrating computer vision (e.g., gesture recognition and artificial intelligence), molecular dynamics (e.g., molecular modelling and material modelling), and intuitive understating of atomic scale. By doing this, we will be able to interact with moleculesvia our senses.


## 4.2.1. Calibration of the Gesture Application with the MD Simulation Model

Figure5 shows the progressive physical movement of the human hand sensed by the gesture application (left column) and transformation of this physical hand movement into the corresponding displacement of the indenter into the substrate (right hand column). The physical distances of the hand movements were captured through a measurement tape and were compared with the movement of the indenter in the MD simulation model so as to see the factor by which the hand movement transfers into the equivalent displacement of the indenter. The relation between these two movements was plotted and is shown in Figure 6. Knowing this relation was a preliminary step for calibrating the newly proposed VAMDS system.Figure 5
Comparison of physical movement with the movement in the MD simulation model (a) at the beginning of the simulation, (b) when the indenter establishes contact with the substrate and (c) final indentation depth was achieved.
(a)
(b)
(c)Figure 6
Relationship between physical movement of hand gesture and indenter in the MD model.


## 4.2.2. HPPT of Silicon and Extraction of the Mechanical Properties fromP-h Plot

The ductile behaviour of brittle materials is often attributed to an event recognized as high-pressure phase transformation (HPPT). HPPT enforces brittle materials to become ductile for the short duration of loading which upon unloading eventually transforms to amorphous phase of the brittle material. HPPT of silicon during its nanoindentation has been explained to be driven from the deviatoric stress rather than temperature [14] in the deformation zone. Much of the literature on silicon nanoindentation have reported HPPT to be the primary mechanism governing the plasticity of silicon that causes brittle-ductile transition. There is however a noticeable exception of Mylvaganam and Zhang [19] observed nanotwinning (associated with Si-I to bct-5 phase transformation) in silicon along the 110 direction that stops at Shockley partial dislocation. It was therefore relevant to assert the occurrence of HPPT in this work for which an indicator known as coordination number was used for its ease of implementation.A snapshot from the MD simulation showing the change in coordination number during peak loading of silicon on various orientations is shown in Figure7. The coordination number of 4 (Figure 7) is representative of bulk silicon while coordination number of 1 and 2 on the free surface of silicon workpiece represents the dangling bonds on the Si-terminated nascent surface. During the nanoindentation, the coordination number in the indentation zone was observed to change from 4 to 6 with a corresponding change in the number of nearest neighbour atoms. Literature suggests that this change corresponds to an allotropic transformation of silicon from the diamond cubic (alpha-silicon) to the body centred tetragonal (beta-silicon) structure and is consistent with the previously reported results of Si-I to Si-II during loading. Most of the metastable phases of silicon were observed to form along the 110 direction, which is the direction of slip plane in silicon [20] as is shown in Figure 7 while the (111) planes act as both slip planes and cleavage planes.Figure 7
Change in coordination number at peak indentation depth of 2 nm on three orientations.Figure8 is a P-h plot obtained from the displacement controlled nanoindentation simulation where the penetration depth of the diamond tool was kept fixed at 2 nm. An interesting observation concerning the P-h plot comes from the work of Jang et al. [21] where a good correlation between unloading curve and the phase of the material formed postindentation of silicon was observed. Based on the micro-Raman spectra, the authors proposed that the unloading discontinuity, often called as “pop-out,” corresponds to the formation of metastable Si-XII/Si-III crystalline phases, while the hysteresis, called the “elbow” is associated with the formation of a-Si. In Figure 8, an interesting feature observed from the P-h plot was that the unloading curve for the (110) and the (111) substrate followed power law whereas the unloading curve on the (010) substrate showed an elbow. This characteristic is a testimony to the formation of a-Si upon the retraction of the indenter on the (010) orientation.Figure 8
P
-
h plots of the nanoindentation for all three crystal orientations of silicon.Another interesting feature observed from theP-h plot obtained from the MD simulation was the residual indentation depth (hf) that can also be used to characterize the recovery of the indented surface upon retraction of the indenter. In this case, hf was quantified from the P-h plot (Figure 8) for all the three simulated orientations and was found inevitably different in all the cases; that is, hf (010) was minimum while hf for the (111) orientation was maximum whereas hf (110) was intermittent. This suggests that at atomic scale, the (010) orientation will show higher extent of elastic recovery than the (111) orientation of silicon. This result is somewhat consistent with the earlier work of Shibata et al. [22] where a Schmidt-type slip orientation factor was proposed and the 11-0 direction either on the (100) or on the (111) planes was recognized as more preferred combination to plastically deform the silicon. Lastly, beside a monotonic increase in the loading curve, a strong cohesion between the indenter and the substrate during retraction was also observed.HPPT is highly direction sensitive [23, 24]; that is, pressure or hydrostatic stress required for transformation on the (111) orientation is lower than that required to induce the transformation on the (100) orientation [20]. Literature suggests that the transformation pressure required to cause metallisation of silicon varies between 9 and 18 GPa [25]. It may however be noted that the hydrostatic pressure differs from deviatoric stress [26]. During nanoindentation, silicon undergoes both shear and compression making it different from the case of simple compressive pressure.As it can be seen from Table3, the hydrostatic stress (pressure) [27] in the indentation zone was observed to be maximum on the (110) orientation followed by (010) orientation and it was minimum for the (111) orientation. This correlation was however not the same for the von Mises stress and Tresca stress in the indentation zone.Table 3
Computational results obtained from the nanoindentation simulation.
Silicon
(010) orientation
(110) orientation
(111) orientation
Average hardness (GPa) at indent depth of 2 nm
16.99
18.22
18.48
Elastic modulus (GPa) at indent depth of 2 nm
129
147
167
h
f
/
h
max
⁡ (extent of elastic recovery)
1.17/2=0.585
1.36/2=0.68
1.47/2=0.735
Peak maximum von Mises stress in the indent zone (GPa)
11.07
22.75
16.72
Peak maximum Tresca stress in the indent zone (GPa)
5.91
12.11
9.14
Peak maximum octahedral shear stress in the indent zone (GPa)
5.21
10.72
7.88
Peak hydrostatic stress in the indent zone (GPa)
20.11
21.98
16.62
Peak temperature in the indentation zone (K)
300
419
438
Features of theP-h curve during unloading
Elbow formation corresponding to a-Si formation
No feature
No featureFor the two kinds of deviatoric stresses, the (110) orientation showed maximum values whereas the (010) orientation showed least value. The peak average temperature in the cutting zone was rather high, especially considering the fact that nanoindentation process was carried out at 10 K. At its peak, the (111) orientation showed the maximum temperature of about 438 K followed by the (110) orientation at 419 K whereas the (010) orientation showed the least temperature of 300 K. Finally, the elastic modulus was obtained fromP-h plot using Oliver and Pharr method and accordingly the variation in E and hardness on all the three crystal surfaces with respect to change in the indentation depth is plotted and shown in Figure 9.Figure 9
Variation in the (a) elastic modulus and (b) nanoindentation hardness of silicon with respect to the indentation depth.
(a)
(b)It is interesting to note from Figure9(a) that the elastic modulus at finite indentation depth of few nanometres decrease sharply until attaining a saturated value as the indentation depth increases. At an indentation depth of 2 nm, the elastic modulus of the three simulated orientations was obtained as 129 GPa for the (010) orientation, 147 GPa for the (110) orientation, and 167 GPa for the (111) orientation. These values are well comparable with the experimental values of the elastic modulus of silicon [28, 29] and resembles similar to what has been obtained from MD [30].It has been seen how a new approach is helpful in carrying out atomic simulations. Future work of this research will realise multiscale approach by integrating computer vision (e.g., gesture recognition and artificial intelligence), molecular dynamics (e.g., molecular modelling and material modelling), and intuitive understating of atomic scale. By doing this, we will be able to interact with moleculesvia our senses.


## 5. Conclusions

This work was aimed at developing a vision-augmented molecular dynamics simulation (VAMDS) method which was tested by performing a case study on nanoindentation of silicon. It was remarkably challenging yet innovative because it will open up new horizons for the development of new knowledge by building a clear base to enable design and creation of novel materials. This was achieved through a systematic design approach integrating the use of a haptic device that fed the data of hand gestures to the simulator. Aside from explaining various aspects of the hardware, design, and development of this vision-augmented technique, the following points can be concluded based on the aforementioned discussions:(1)
The mechanism of ductility in silicon during its nanoindentation and the observed plasticity occurs due to the high pressure phase transformation (HPPT) that is driven by stress rather than temperature. For HPPT to occur in silicon, the von Mises and Tresca stress for the (110) orientation were found maximum while the (010) orientation showed the least values whereas the (111) orientation showed somewhat an intermittent value.(2)
At an indentation depth of 2 nm, it was seen that the unloading curve for the (110) and the (111) substrate follows power law whereas the unloading curve on the (010) substrate showed an elbow which was found related with the formation of a-Si upon the retraction of the indenter on the (010) orientation only. Furthermore, the elastic constants of silicon varied in the order ofE111>E110>E010.(3)
Obtaining the variation in the hardness and elastic modulus using classical way of doing molecular dynamics simulation would have been a relatively tedious task but, with the implementation of the VAMD, this is easily achievable. The use of haptic device as shown in this work reduces the pain of obtaining mechanical propertyviaP-h profile as the simulation run can be restarted from any point merely by the use of hand gestures which is the major novelty of this work.



---

*Source: 857574-2015-06-23.xml*