# Molecular Dynamics Study on Lubrication Mechanism in Crystalline Structure between Copper and Sulfur

Authors: Ken-ichi Saitoh; Tomohiro Sato; Masanori Takuma; Yoshimasa Takahashi; Ryuketsu Chin
Journal: Journal of Materials
 (2015)

---

## Abstract

To clarify the nanosized mechanism of good lubrication in copper disulfide (Cu2S) crystal which is used as a sliding material, atomistic modeling of Cu2S is conducted and molecular dynamics (MD) simulations are performed in this paper. The interatomic interaction between atoms and crystalline structure in the phase of hexagonal crystal of Cu2S are carefully estimated by first-principle calculations. Then, approximating these interactions, we originally construct a conventional interatomic potential function of Cu2S crystal in its hexagonal phase. By using this potential function, we perform MD simulation of Cu2S crystal which is subjected to shear loading parallel to the basal plane. We compare results obtained by different conditions of sliding directions. Unlike ordinary hexagonal metallic crystals, it is found that the easy-glide direction does not always show small shear stress for Cu2S crystal. Besides, it is found that shearing velocity affects largely the magnitude of averaged shear stress. Generally speaking, higher velocity results in higher resistance against shear deformation. As a result, it is understood that Cu2S crystal exhibits somewhat liquid-like (amorphous) behavior in sliding condition and shear resistance increases with increase of sliding speed.


---

## Body

## 1. Introduction

In recent years, using lead (Pb) in metallic material is strictly prohibited due to risk to human health. That is why lead-free alloy is now being sought for mechanical materials. It has been recognized that bronze, which is composed of copper (Cu) and tin (Sn) atoms, forms a material with good performance of lubrication by the addition of Pb atoms. Therefore, Cu-Sn alloys are widely used for sliding mechanical elements, such as brake friction materials (brake-pads) or mechanical bearing materials. Broadly speaking, it is understood that Pb solids basically have low melting point and they easily tend to be fluidic under high temperature and high pressure condition, as in the environment inside the sliding equipment. One example of substitution of Pb is bismuth (Bi), but it is regarded as one of rare metals and it is hard to obtain much. Therefore, it requires us to clear the hurdle concerning both cost and safety. Recently, it is reported that bronze with the addition of sulfur (S) atoms shows very good performance of lubrication in sliding material, and it has already been patented and been actually developed as an engineering material for industrial use [1]. Historically, it has been well known that S atoms provide a nice machinability (in cutting) of copper alloy, because they work as “tip-breaker” in cutting process [2]. Besides, S atom is well accepted as a component of solid lubricant like Pb atom. For example, molybdenum disulfide, MoS2, including many S atoms, exhibits high performance as the solid lubricant. Therefore, alloys made of Cu and S may provide good sliding and friction performance. In those alloys, S atom will play a role of the long-standing Pb atom instead. Nevertheless, the lubricating mechanism in Cu-S system has not been clarified yet and we should uncover it by studying further.Hirai et al. have succeeded in developing a bronze made with microscopic dispersion of copper sulfide (Cu2S) and they are now proceeding to its industrial use [3]. Akamatsu et al. are also exploring the fabricating method of Cu2S alloy by the motivation of avoiding the risk to human health [4]. In the field of mineralogy, Cu2S crystal is well known as a mineral called chalcocite. Recently, since the Cu2S structure has a unique electronic behavior as for ionic conductivity, it is planned to be used for switching device applications, such as a solar power plant material [5] and a superconductivity material [6, 7]. In the context of unique sliding and friction behavior of Cu2S crystal already found experimentally, the present authors have studied it by atomistic modeling. One of stable stoichiometries of copper and sulfide system is found at Cu2S (Cu : S = 2 : 1). Interestingly, this Cu2S compound exhibits very complicated solid phase transition behavior [8, 9]. Around the standard temperature, the unit structure of Cu2S is orthorhombic, but it transforms into hexagonal symmetry above 105°C (degrees Celsius). Moreover, in higher temperature above 460°C, it is transformed into cubic crystal.When copper alloys are used for sliding and friction materials, they are generally placed in a severe condition of high temperature beyond several hundred °C, and therefore the behavior of hexagonal crystal (HC) would be worth studying extensively. In Cu2S crystal, consequently, one S atom has two neighbor Cu atoms. So, they tend to organize a triangular conformation and keep their firm interatomic bonding when they are stabilized. The same mechanism concerning S atom is also suggested in MoS2 system which is famous as good solid lubricant, where the S atom makes firm bonding with metal atoms [10].Other than chalcocite (Cu2S), there are some well-known natural minerals based on copper alloy [11, 12], such as bornite (Cu5FeS4), chalcopyrite (CuFeS2), covellite (CuS), and digenite (Cu9S5). All of these alloys are interesting as for interrelationship in their fabrication, and all of them may be applied to new functional material in the future. Among them, however, Cu2S crystal is still a main compound, so it should be focused on particularly. Thus, in this paper, we will study the atomistic structure of Cu2S in HC and discuss its atomistic behavior in sliding and friction process.In the present study, aiming at unempirical approach, the first-principle (ab initio) calculation is utilized to identify the crystal structure of Cu2S. So far, there have been some studies of copper sulfides byab initio method, which are mainly employing density functional theory (DFT) [5, 13, 14]. The electronic structure has been reproduced, but unfortunately there still remains ambiguity about atomic position of coppers in the unit crystal, and so further modeling and discussion will be required. Therefore, in this study, by using conventional DFT software package Wien2k [15], the optimization of the unit structure of the Cu2S in HC is carried out based on Kashida et al.’s structural model [13]. As a result, we obtain its energetically stable structure and we determine its basic crystalline information, that is, lattice constants as for lengths, a and c, assuming the HC configuration.In order to understand the dynamics of sliding or friction occurring in Cu2S crystal, molecular dynamics (MD) simulation is suitable. Lattice constants a and c and ground-state energies which are obtained byab initio calculations are used in the determination of interatomic potential for MD simulation. In the MD simulation, an existing interatomic potential function of MoS2 [10] is referred to for its function or functional form. This seems reasonable because those two crystals of metal sulfide (Cu2S and MoS2) have many common features, as described above.Thus, we simulate atomistically the sliding of Cu2S crystal by using MD method with an original and nonempirical interatomic potential. The computation model is built by stacking layers on slip plane (basal plane) of HC lattice and it is at the first place compressed and then shear deformation is applied. Here, the compression is crucial since in the actual situation the friction behavior in sliding is to take place under severe contact condition with large compressive stress (pressure on the contacting surface). We will focus on the dependence on crystalline direction of sliding or on sliding velocity in such severe condition.This paper is organized as follows. Following the Introduction, we describe the basic information of crystal structure of Cu2S and show some theories needed for the modeling and the calculation. Then, the essence of theab initio calculation is shown, and the procedure to implement the interatomic potential is briefly presented. Thereafter, the MD model used in sliding simulation is explained. The summary of results inab initio calculation and MD results of sliding simulation are shown and their discussion is made. Finally, the conclusion of this paper will be made.


## 2. Theory and Method

### 2.1. Atomic Structure of Cu2S Crystal

As pure compounds of copper sulfide, there are covellite (CuS) and chalcocite (Cu2S), and the latter is focused on now. For Cu2S, orthorhombic crystal is obtained in room temperature, but, beyond 105°C, it transforms into hexagonal crystal (HC). Furthermore, in higher temperature than 460°C, it changes to cubic crystal [8]. Copper alloys are often subject to very high pressure and relatively high temperature during operation, such as in the surface layer of mechanical bearing. Thus, the present simulation should focus on the temperature range around several hundred kelvins. Accordingly, the Cu2S material here is assumed to be in HC. HC is also taken in MoS2 and carbon allotrope (graphite or nanometer-thickness graphene), and both of them are well-known solid lubricants. Among conventional metals, however, HC metal such as magnesium (Mg) is relatively inferior in deformation due to the lack of crystalline slip systems, when compared with face-centered cubic (fcc) crystals such as aluminum (Al) and Cu. Only the basal slip of hexagonal crystal is always well activated in Mg crystal, and so its slip is highly direction dependent. It is guessed that sliding motion in HC structure will prefer such direction-dependent feature. This leads to the fact that compounds with HC structure, such as Cu2S and MoS2, will be good at lubrication.


### 2.2. First-Principle Calculation

Cu atoms are metal and include metallic bonds, but, on the other hand, S atoms are nonmetal and show large electronegativity. Accordingly, bond between atoms inside Cu2S seems very complicated. Possibly, ionic bond and metallic bond would be nonlinearly interrelated. In considering the structure and bonding state of solid crystal, it is helpful to reproduce and observe the detailed electron density distribution. Besides, those results can be taken into account in the construction of accurate interatomic potential for MD simulations. First-principle (FP) calculation (in Latin words,ab initio calculation) refers to all quantum-mechanical computational chemistry methods and is sometimes called “band calculation.” It obtains electron density around atomic nuclei based on the periodic crystal structure of constituent atoms. Nowadays, a number of software packages of FP method are available in many researches. Among those software packages, our choice is WIEN2k package [15], which is based on density functional theory (DFT). The characteristic of WIEN2k is that it uses full-potential formulation, instead of muffin-tin approximation, and that it uses the linearized augmented plane wave basis set with local orbitals (LAPW + lo). As for correlation exchange term, it can apply either local density approximation (LDA) or generalized gradient approximation (GGA). The detailed theory and usage are not the main focus of this study and should be omitted here.While all the S atoms in Cu2S are simply located at the lattice point of HC, it is suggested that Cu atoms seem to have large mobility inside the unit cell, especially in high temperature, keeping their symmetric atomic configuration. Since actual position of Cu atoms in Cu2S is quite complicated and is, in fact, giving controversy, the present study employs existing crystal model used in previous computation of Kashida et al. [13] It is confirmed that the model agrees with experimental fact as well as theoretical one. Initial atomic position of Cu2S in HC is shown in Table 1. Figure 1 shows initial position of atoms in Cu2S crystal unit.Table 1
Cu2S atomic coordinates and crystal data.
(a)
Element
Coordinate (relative to unit cell)
x
y
z
S
1/3
2/3
1/4
Cu(1)
0
0
1/4
Cu(2)
1/3
2/3
0.578
(b)
Property
Unit
Value
Lattice constant(1)a
nm
0.389
Lattice constant(2)c
nm
0.688
Lattice constant(3)α
deg.
90.0
Lattice constant(4)β
deg.
90.0
Lattice constant(5)γ
deg.
120.0
Atoms per unit cell
—
3
Space group
—
P63Figure 1
Crystal structure of Cu2S (hexagonal crystal unit: 1/3 of conventional hexagonal prism).In order to obtain the minimum (ground-state) energy, the atomic configuration is adjusted via the function of optimization in the software. Lattice constants for angle (α,β,γ) are fixed so that the HC lattice should be maintained. On the other hand, lattice constants for lengths (a,c) are changed as follows. Keeping a certain ratio of c/a, the lengths (a,c) are varied from compressive state to tensile one. The deviation of c/a from the initial value is +4~+10%. The volumetric change by the change of (a,c) ranges within −2~+6%. The exchange correlation energy of PBE-GGA is employed. Judgement of the energy convergence is done at 0.0001 Ry, where 1 Ry = 2.180×10-18J = 13.60 eV. Table 2 summarizes calculation conditions and parameters used in the present study. The number of k-points is varied, and approximately above 20000 we may say that the total energy will be converged and shows sufficient accuracy. Therefore, we recognize that 46×46×22=46552k-points are large enough and this will be used through the present study.Table 2
Calculation condition for structural optimization of Cu2S crystal.
Exchange correlation potential
—
PBE-GGA
Energy to separate core from valence states
Ry∗
−6.0
Energy convergence criteria
Ry∗
0.0001
Change ofc/a ratio
Rate of change
%
+4~+10
Increment of rate (forc/a)
%
+2
Volume change
Rate of change
%
−2~+6
Increment of rate (forc/a)
%
+1
∗Energy unit: 1 Ry = 2.180×10-18 J = 13.60 eV.


### 2.3. The Potential Function of Cu2S for MD

As mentioned above, copper disulfide, Cu2S, has a variety of crystalline structures depending on the temperature. The structure is orthorhombic in standard temperature, HC above 105°C, and cubic further above 460°C. Here, as HC is assumed and the composition ratio is Cu : S = 2 : 1, one S atom should correspond to two Cu atoms. Therefore, those three atoms compose a structural unit, and the three-body bond is formed as a function of their angle. By the way, MoS2 has HC structure as well. The structural feature common to those HC lattices is that their slip plane is only limited to the basal plane (0001). This is completely different from fcc or bcc, which is accompanying many slip systems. In fact, the slip of HC lattice takes place only when a certain condition has been fulfilled. Thus, the potential function which was proposed for MoS2 [10] can be utilized as the reference to construct the function needed in the present study for Cu2S. Though MoS2 and Cu2S have the contrary ratio as to S atom, a lot of common features are likely to be found in these two crystals. For example, angular-dependent three-body interaction is crucial for stabilization of both crystals. To begin with, we apply the potential forms used in the MoS2 study, and then our results obtained by first-principle calculation are properly fitted to a new potential function for Cu2S. We insist that since the fitting data is obtained theoretically only by DFT calculations, the new potential function will be an unempirical (ab initio) one.The detail of the new potential function for Cu2S is as follows. The potential function is a summation of three distinct potential functions, as shown in (1). Ionic and metallic bonds are represented by Born-Mayer-Huggins (BMH) potential, ϕBMH, and Morse potential, ϕMorse, as expressed in (2) and (3), respectively. Angular-dependent potential for the bond of Cu-S interaction, ϕangle, is represented as (4), which is three-body potential. Parameters for those potential functions are summarized in Table 3:(1)ϕ=ϕBMH+ϕMorse+ϕangle,(2)ϕBMH=zizje2r+f0bi+bjexp⁡ai+aj-rbi+bj+cicjr6,(3)ϕMorse=Dijexp⁡-2βijrij-r0-2exp⁡-βijrij-r0,(4)ϕangle=Hθθ-θ02.In those equations, zi or zj is electronic charge of each atom, and the values zS=-2.0 and zCu=+1.0 are employed for S and Cu atoms, respectively. f0 = 1 [kcal/Å] = 4.186 [kJ/Å] is the parameter which determines the stiffness of each ionic sphere and is transferred from the study of MoS2. In the Morse potential, the equilibrium distance r0 and bond energy Dij are fitted to those in the optimized structure which is obtained in ourab initio calculation. In the angular potential, equilibrium three-body angle θ0 is also determined by ourab initio calculation. But, it is difficult to derive the spring constant Hθ just from ourab initio calculation, so the same value as in MoS2 is applied here. Of course, the value of Hθ influences rigidity. But its function form is basically harmonic and we found that the effect of spring constant on sliding is not so crucial.Table 3
Interatomic potential parameters for Cu2S crystal.
(a)
B-M-H potential
Atoms
a
i
b
i
c
i
Unit
[Å]
[Å]
[kJ0.5⋅Å3/mol0.5]
Cu
0.696
0.085
0.000
S
1.831
0.085
70.000
(b)
Morse potential
Application
D
i
j
β
i
j
r
0
Unit
[kJ/mol]
[1/Å]
[Å]
Cu-S
152.076
1.800
2.532
(c)
Angular potential
Application
H
θ
θ
0
Unit
[kJ/mol]
[deg.]
Cu-S-Cu
196.648
98.943
S-Cu-S
196.648
98.943


### 2.4. Molecular Dynamics Model for Cu2S Sliding

Figure2 shows the MD model for sliding simulation. In the first place, the optimized atomic configuration obtained byab initio calculation is prepared and then is thermally equilibrated at a finite temperature. Then, the structure is compressed in z direction as shown in Figure 2. After that, the atomic velocity in the x direction of two special regions (shown with yellow color in the figure) is constrained so that shear deformation is applied to the whole region in the direction parallel to the slip plane. These processes are supposed to mimic a Cu2S crystal in sliding condition. This kind of dynamic sliding simulation has been often studied for tribology-system as found in literatures so far [16], and its advantage is that we can observe directly the sliding behavior of atomic system.Figure 2
Abstract of MD sliding model of Cu2S crystal.It is supposed that crystalline slip plane of Cu2S plays an important role in the slip deformation. As stated above, Cu2S crystal unit has the HC conformation and the primary slip plane should be (0001) in Miller’s indices. The primary slip direction should be the shortest interatomic distance existing on the (0001) plane, that is, 21¯1¯0. We will call it (b)-direction hereafter. The second easiest slip direction should be 01¯10, which we call (a)-direction hereafter. In the present study, we configure these two possible slip directions for sliding simulation and compare these results.Other factors affecting sliding behavior of Cu2S would be velocity and temperature. Sliding velocities are varied so that corresponding shear velocities are 0.1, 1.0, and 5.0 m/s. The reader should feel that these velocities seem quite higher than used in actual sliding machine or equipment. However, in the MD simulation, time increment has to be very short such as 1 fs (femtosecond) or less, so such large deformation rate is inevitable.


### 2.5. Analyzing MD Results

In the present study, the MD results are mainly analyzed by radial distribution function (RDF) and atomic shear stress. Strict formalization of stress tensor with regard to three-body term would require further discussion at this stage [17], so we use here an approximate and practial treatment. These analyzing methods for MD results are explained below.


#### 2.5.1. Crystal Structure Analysis by Radial Distribution Function (RDF)

Usually, solid crystal has a three-dimensional order. The Cu2S crystal used in the present study exhibits HC structure. One of the methods to identify the crystalline order is radial distribution function (RDF) g(r) [18]. Briefly speaking, the function g(r) gives the number density of neighbor atoms which are found in a certain distance from each atom. Its distribution shows some peaks at some specific distances inherent to the crystalline order. When the temperature arises or when the crystalline order is partially lost due to outbreak of lattice defects, the peak value of RDF generally reduces and the shape at the peak becomes broadened. Thus, by observing the change of the peak position and the shape of the curve of RDF, the deterioration of the crystalline order can be detected and occurrence of lattice defects can be guessed.


#### 2.5.2. Atomic Stress Analysis (Approximate Treatment including Three-Body Terms)

Since our MD model is in solid state, the behavior of sliding and friction is reflected by a shear stress averaged over the specimen. Therefore, atomic-scale stress is being analysed. But we note that stress (or strain) is basically the concept in continuum mechanics and is not essential for any atomic system. The facts that atomic positions are discrete and that interatomic potential has nonlocal nature will conflict with the local theory of continuum. Accordingly, we need to formulate atomic stress together with some approximation.Suppose that total potential energy inside the system comprises all the superposition of pairwise interactionsϕ(r). When any selected atom i conducts homogeneous deformation with regard to the relative position to another neighbor atom j, strain at the interatomic space between i and j can be presented by their interatomic distance r(i,j) and interatomic vector r(i,j) (this assumes that there is no slip, no phase transition, nor long atomic diffusion over the lattice). Then, atomic stress tensor occurring at the atom i is expressed as(5)σi=12Ωi∑j∂ϕi,j∂ri,jri,j⊗ri,jri,j,where the summation is composed of the work done by atomic motion, which is called “virial.” Ω(i) in the denominator is the volume supposedly occupied by the atom i, so (5) is regarded as the energy density around that atom. In the present study, Ω(i) is assumed to be a constant, which has been estimated in the reference structure (which has been estimated in the undeformed lattice). Additionally, vector product (r⊗r) means the dyadic between interatomic vectors, so (5) has the form of the second-order (symmetric) tensor and its components are six. Thus, if an interatomic potential of MD contains just pairwise term, the virial can be straightforwardly converted to atomic stress. But if it includes three-body terms depending on the angle as in the present interatomic potential, they can never strictly be decomposed into virials of each pair of atoms. However, we can apply an approximate formulation to a triplet, i, j, and k (i,j,k). In practice, the combination of virials of each pair, (i,j), (j,k), and (k,i), is used as the net virial of the triplet. As a result, the formulation of atomic stress comes to the expression(6)σαβi=12Ωi∑j∂ϕpairi,j∂ri,jrαi,jrβi,jri,j+13Ωi∑j,k∂ϕanglei,j∂ri,jrαi,jrβi,jri,j+∂ϕanglej,k∂rj,krαj,krβj,krj,k+∂ϕanglek,i∂rk,irαk,irβk,irk,i,where the first term corresponds to pair interaction ϕpair(r) and the second term is for three-body interaction. In this expression, Ω(i), ϕpair(i,j), ϕangle(i,j), r(i,j), and rα(i,j)(α,β=1,2,3) are atomic volume, pairwise potential function, three-body potential function, interatomic distance (between i and j), and its vector components. In comparing all the components in (6) during actual MD simulation, it is understood that three-body term is relatively large so that it is not negligible. This expression of atomic stress including three-body term is also utilized in widely used MD software (e.g., LAMMPS) [19]. Besides, this estimation of atomic stress was found appropriate when we applied it to another MD study using the three-body-type Tersoff potential (for Si-Ge system) [20].


## 2.1. Atomic Structure of Cu2S Crystal

As pure compounds of copper sulfide, there are covellite (CuS) and chalcocite (Cu2S), and the latter is focused on now. For Cu2S, orthorhombic crystal is obtained in room temperature, but, beyond 105°C, it transforms into hexagonal crystal (HC). Furthermore, in higher temperature than 460°C, it changes to cubic crystal [8]. Copper alloys are often subject to very high pressure and relatively high temperature during operation, such as in the surface layer of mechanical bearing. Thus, the present simulation should focus on the temperature range around several hundred kelvins. Accordingly, the Cu2S material here is assumed to be in HC. HC is also taken in MoS2 and carbon allotrope (graphite or nanometer-thickness graphene), and both of them are well-known solid lubricants. Among conventional metals, however, HC metal such as magnesium (Mg) is relatively inferior in deformation due to the lack of crystalline slip systems, when compared with face-centered cubic (fcc) crystals such as aluminum (Al) and Cu. Only the basal slip of hexagonal crystal is always well activated in Mg crystal, and so its slip is highly direction dependent. It is guessed that sliding motion in HC structure will prefer such direction-dependent feature. This leads to the fact that compounds with HC structure, such as Cu2S and MoS2, will be good at lubrication.


## 2.2. First-Principle Calculation

Cu atoms are metal and include metallic bonds, but, on the other hand, S atoms are nonmetal and show large electronegativity. Accordingly, bond between atoms inside Cu2S seems very complicated. Possibly, ionic bond and metallic bond would be nonlinearly interrelated. In considering the structure and bonding state of solid crystal, it is helpful to reproduce and observe the detailed electron density distribution. Besides, those results can be taken into account in the construction of accurate interatomic potential for MD simulations. First-principle (FP) calculation (in Latin words,ab initio calculation) refers to all quantum-mechanical computational chemistry methods and is sometimes called “band calculation.” It obtains electron density around atomic nuclei based on the periodic crystal structure of constituent atoms. Nowadays, a number of software packages of FP method are available in many researches. Among those software packages, our choice is WIEN2k package [15], which is based on density functional theory (DFT). The characteristic of WIEN2k is that it uses full-potential formulation, instead of muffin-tin approximation, and that it uses the linearized augmented plane wave basis set with local orbitals (LAPW + lo). As for correlation exchange term, it can apply either local density approximation (LDA) or generalized gradient approximation (GGA). The detailed theory and usage are not the main focus of this study and should be omitted here.While all the S atoms in Cu2S are simply located at the lattice point of HC, it is suggested that Cu atoms seem to have large mobility inside the unit cell, especially in high temperature, keeping their symmetric atomic configuration. Since actual position of Cu atoms in Cu2S is quite complicated and is, in fact, giving controversy, the present study employs existing crystal model used in previous computation of Kashida et al. [13] It is confirmed that the model agrees with experimental fact as well as theoretical one. Initial atomic position of Cu2S in HC is shown in Table 1. Figure 1 shows initial position of atoms in Cu2S crystal unit.Table 1
Cu2S atomic coordinates and crystal data.
(a)
Element
Coordinate (relative to unit cell)
x
y
z
S
1/3
2/3
1/4
Cu(1)
0
0
1/4
Cu(2)
1/3
2/3
0.578
(b)
Property
Unit
Value
Lattice constant(1)a
nm
0.389
Lattice constant(2)c
nm
0.688
Lattice constant(3)α
deg.
90.0
Lattice constant(4)β
deg.
90.0
Lattice constant(5)γ
deg.
120.0
Atoms per unit cell
—
3
Space group
—
P63Figure 1
Crystal structure of Cu2S (hexagonal crystal unit: 1/3 of conventional hexagonal prism).In order to obtain the minimum (ground-state) energy, the atomic configuration is adjusted via the function of optimization in the software. Lattice constants for angle (α,β,γ) are fixed so that the HC lattice should be maintained. On the other hand, lattice constants for lengths (a,c) are changed as follows. Keeping a certain ratio of c/a, the lengths (a,c) are varied from compressive state to tensile one. The deviation of c/a from the initial value is +4~+10%. The volumetric change by the change of (a,c) ranges within −2~+6%. The exchange correlation energy of PBE-GGA is employed. Judgement of the energy convergence is done at 0.0001 Ry, where 1 Ry = 2.180×10-18J = 13.60 eV. Table 2 summarizes calculation conditions and parameters used in the present study. The number of k-points is varied, and approximately above 20000 we may say that the total energy will be converged and shows sufficient accuracy. Therefore, we recognize that 46×46×22=46552k-points are large enough and this will be used through the present study.Table 2
Calculation condition for structural optimization of Cu2S crystal.
Exchange correlation potential
—
PBE-GGA
Energy to separate core from valence states
Ry∗
−6.0
Energy convergence criteria
Ry∗
0.0001
Change ofc/a ratio
Rate of change
%
+4~+10
Increment of rate (forc/a)
%
+2
Volume change
Rate of change
%
−2~+6
Increment of rate (forc/a)
%
+1
∗Energy unit: 1 Ry = 2.180×10-18 J = 13.60 eV.


## 2.3. The Potential Function of Cu2S for MD

As mentioned above, copper disulfide, Cu2S, has a variety of crystalline structures depending on the temperature. The structure is orthorhombic in standard temperature, HC above 105°C, and cubic further above 460°C. Here, as HC is assumed and the composition ratio is Cu : S = 2 : 1, one S atom should correspond to two Cu atoms. Therefore, those three atoms compose a structural unit, and the three-body bond is formed as a function of their angle. By the way, MoS2 has HC structure as well. The structural feature common to those HC lattices is that their slip plane is only limited to the basal plane (0001). This is completely different from fcc or bcc, which is accompanying many slip systems. In fact, the slip of HC lattice takes place only when a certain condition has been fulfilled. Thus, the potential function which was proposed for MoS2 [10] can be utilized as the reference to construct the function needed in the present study for Cu2S. Though MoS2 and Cu2S have the contrary ratio as to S atom, a lot of common features are likely to be found in these two crystals. For example, angular-dependent three-body interaction is crucial for stabilization of both crystals. To begin with, we apply the potential forms used in the MoS2 study, and then our results obtained by first-principle calculation are properly fitted to a new potential function for Cu2S. We insist that since the fitting data is obtained theoretically only by DFT calculations, the new potential function will be an unempirical (ab initio) one.The detail of the new potential function for Cu2S is as follows. The potential function is a summation of three distinct potential functions, as shown in (1). Ionic and metallic bonds are represented by Born-Mayer-Huggins (BMH) potential, ϕBMH, and Morse potential, ϕMorse, as expressed in (2) and (3), respectively. Angular-dependent potential for the bond of Cu-S interaction, ϕangle, is represented as (4), which is three-body potential. Parameters for those potential functions are summarized in Table 3:(1)ϕ=ϕBMH+ϕMorse+ϕangle,(2)ϕBMH=zizje2r+f0bi+bjexp⁡ai+aj-rbi+bj+cicjr6,(3)ϕMorse=Dijexp⁡-2βijrij-r0-2exp⁡-βijrij-r0,(4)ϕangle=Hθθ-θ02.In those equations, zi or zj is electronic charge of each atom, and the values zS=-2.0 and zCu=+1.0 are employed for S and Cu atoms, respectively. f0 = 1 [kcal/Å] = 4.186 [kJ/Å] is the parameter which determines the stiffness of each ionic sphere and is transferred from the study of MoS2. In the Morse potential, the equilibrium distance r0 and bond energy Dij are fitted to those in the optimized structure which is obtained in ourab initio calculation. In the angular potential, equilibrium three-body angle θ0 is also determined by ourab initio calculation. But, it is difficult to derive the spring constant Hθ just from ourab initio calculation, so the same value as in MoS2 is applied here. Of course, the value of Hθ influences rigidity. But its function form is basically harmonic and we found that the effect of spring constant on sliding is not so crucial.Table 3
Interatomic potential parameters for Cu2S crystal.
(a)
B-M-H potential
Atoms
a
i
b
i
c
i
Unit
[Å]
[Å]
[kJ0.5⋅Å3/mol0.5]
Cu
0.696
0.085
0.000
S
1.831
0.085
70.000
(b)
Morse potential
Application
D
i
j
β
i
j
r
0
Unit
[kJ/mol]
[1/Å]
[Å]
Cu-S
152.076
1.800
2.532
(c)
Angular potential
Application
H
θ
θ
0
Unit
[kJ/mol]
[deg.]
Cu-S-Cu
196.648
98.943
S-Cu-S
196.648
98.943


## 2.4. Molecular Dynamics Model for Cu2S Sliding

Figure2 shows the MD model for sliding simulation. In the first place, the optimized atomic configuration obtained byab initio calculation is prepared and then is thermally equilibrated at a finite temperature. Then, the structure is compressed in z direction as shown in Figure 2. After that, the atomic velocity in the x direction of two special regions (shown with yellow color in the figure) is constrained so that shear deformation is applied to the whole region in the direction parallel to the slip plane. These processes are supposed to mimic a Cu2S crystal in sliding condition. This kind of dynamic sliding simulation has been often studied for tribology-system as found in literatures so far [16], and its advantage is that we can observe directly the sliding behavior of atomic system.Figure 2
Abstract of MD sliding model of Cu2S crystal.It is supposed that crystalline slip plane of Cu2S plays an important role in the slip deformation. As stated above, Cu2S crystal unit has the HC conformation and the primary slip plane should be (0001) in Miller’s indices. The primary slip direction should be the shortest interatomic distance existing on the (0001) plane, that is, 21¯1¯0. We will call it (b)-direction hereafter. The second easiest slip direction should be 01¯10, which we call (a)-direction hereafter. In the present study, we configure these two possible slip directions for sliding simulation and compare these results.Other factors affecting sliding behavior of Cu2S would be velocity and temperature. Sliding velocities are varied so that corresponding shear velocities are 0.1, 1.0, and 5.0 m/s. The reader should feel that these velocities seem quite higher than used in actual sliding machine or equipment. However, in the MD simulation, time increment has to be very short such as 1 fs (femtosecond) or less, so such large deformation rate is inevitable.


## 2.5. Analyzing MD Results

In the present study, the MD results are mainly analyzed by radial distribution function (RDF) and atomic shear stress. Strict formalization of stress tensor with regard to three-body term would require further discussion at this stage [17], so we use here an approximate and practial treatment. These analyzing methods for MD results are explained below.


### 2.5.1. Crystal Structure Analysis by Radial Distribution Function (RDF)

Usually, solid crystal has a three-dimensional order. The Cu2S crystal used in the present study exhibits HC structure. One of the methods to identify the crystalline order is radial distribution function (RDF) g(r) [18]. Briefly speaking, the function g(r) gives the number density of neighbor atoms which are found in a certain distance from each atom. Its distribution shows some peaks at some specific distances inherent to the crystalline order. When the temperature arises or when the crystalline order is partially lost due to outbreak of lattice defects, the peak value of RDF generally reduces and the shape at the peak becomes broadened. Thus, by observing the change of the peak position and the shape of the curve of RDF, the deterioration of the crystalline order can be detected and occurrence of lattice defects can be guessed.


### 2.5.2. Atomic Stress Analysis (Approximate Treatment including Three-Body Terms)

Since our MD model is in solid state, the behavior of sliding and friction is reflected by a shear stress averaged over the specimen. Therefore, atomic-scale stress is being analysed. But we note that stress (or strain) is basically the concept in continuum mechanics and is not essential for any atomic system. The facts that atomic positions are discrete and that interatomic potential has nonlocal nature will conflict with the local theory of continuum. Accordingly, we need to formulate atomic stress together with some approximation.Suppose that total potential energy inside the system comprises all the superposition of pairwise interactionsϕ(r). When any selected atom i conducts homogeneous deformation with regard to the relative position to another neighbor atom j, strain at the interatomic space between i and j can be presented by their interatomic distance r(i,j) and interatomic vector r(i,j) (this assumes that there is no slip, no phase transition, nor long atomic diffusion over the lattice). Then, atomic stress tensor occurring at the atom i is expressed as(5)σi=12Ωi∑j∂ϕi,j∂ri,jri,j⊗ri,jri,j,where the summation is composed of the work done by atomic motion, which is called “virial.” Ω(i) in the denominator is the volume supposedly occupied by the atom i, so (5) is regarded as the energy density around that atom. In the present study, Ω(i) is assumed to be a constant, which has been estimated in the reference structure (which has been estimated in the undeformed lattice). Additionally, vector product (r⊗r) means the dyadic between interatomic vectors, so (5) has the form of the second-order (symmetric) tensor and its components are six. Thus, if an interatomic potential of MD contains just pairwise term, the virial can be straightforwardly converted to atomic stress. But if it includes three-body terms depending on the angle as in the present interatomic potential, they can never strictly be decomposed into virials of each pair of atoms. However, we can apply an approximate formulation to a triplet, i, j, and k (i,j,k). In practice, the combination of virials of each pair, (i,j), (j,k), and (k,i), is used as the net virial of the triplet. As a result, the formulation of atomic stress comes to the expression(6)σαβi=12Ωi∑j∂ϕpairi,j∂ri,jrαi,jrβi,jri,j+13Ωi∑j,k∂ϕanglei,j∂ri,jrαi,jrβi,jri,j+∂ϕanglej,k∂rj,krαj,krβj,krj,k+∂ϕanglek,i∂rk,irαk,irβk,irk,i,where the first term corresponds to pair interaction ϕpair(r) and the second term is for three-body interaction. In this expression, Ω(i), ϕpair(i,j), ϕangle(i,j), r(i,j), and rα(i,j)(α,β=1,2,3) are atomic volume, pairwise potential function, three-body potential function, interatomic distance (between i and j), and its vector components. In comparing all the components in (6) during actual MD simulation, it is understood that three-body term is relatively large so that it is not negligible. This expression of atomic stress including three-body term is also utilized in widely used MD software (e.g., LAMMPS) [19]. Besides, this estimation of atomic stress was found appropriate when we applied it to another MD study using the three-body-type Tersoff potential (for Si-Ge system) [20].


## 2.5.1. Crystal Structure Analysis by Radial Distribution Function (RDF)

Usually, solid crystal has a three-dimensional order. The Cu2S crystal used in the present study exhibits HC structure. One of the methods to identify the crystalline order is radial distribution function (RDF) g(r) [18]. Briefly speaking, the function g(r) gives the number density of neighbor atoms which are found in a certain distance from each atom. Its distribution shows some peaks at some specific distances inherent to the crystalline order. When the temperature arises or when the crystalline order is partially lost due to outbreak of lattice defects, the peak value of RDF generally reduces and the shape at the peak becomes broadened. Thus, by observing the change of the peak position and the shape of the curve of RDF, the deterioration of the crystalline order can be detected and occurrence of lattice defects can be guessed.


## 2.5.2. Atomic Stress Analysis (Approximate Treatment including Three-Body Terms)

Since our MD model is in solid state, the behavior of sliding and friction is reflected by a shear stress averaged over the specimen. Therefore, atomic-scale stress is being analysed. But we note that stress (or strain) is basically the concept in continuum mechanics and is not essential for any atomic system. The facts that atomic positions are discrete and that interatomic potential has nonlocal nature will conflict with the local theory of continuum. Accordingly, we need to formulate atomic stress together with some approximation.Suppose that total potential energy inside the system comprises all the superposition of pairwise interactionsϕ(r). When any selected atom i conducts homogeneous deformation with regard to the relative position to another neighbor atom j, strain at the interatomic space between i and j can be presented by their interatomic distance r(i,j) and interatomic vector r(i,j) (this assumes that there is no slip, no phase transition, nor long atomic diffusion over the lattice). Then, atomic stress tensor occurring at the atom i is expressed as(5)σi=12Ωi∑j∂ϕi,j∂ri,jri,j⊗ri,jri,j,where the summation is composed of the work done by atomic motion, which is called “virial.” Ω(i) in the denominator is the volume supposedly occupied by the atom i, so (5) is regarded as the energy density around that atom. In the present study, Ω(i) is assumed to be a constant, which has been estimated in the reference structure (which has been estimated in the undeformed lattice). Additionally, vector product (r⊗r) means the dyadic between interatomic vectors, so (5) has the form of the second-order (symmetric) tensor and its components are six. Thus, if an interatomic potential of MD contains just pairwise term, the virial can be straightforwardly converted to atomic stress. But if it includes three-body terms depending on the angle as in the present interatomic potential, they can never strictly be decomposed into virials of each pair of atoms. However, we can apply an approximate formulation to a triplet, i, j, and k (i,j,k). In practice, the combination of virials of each pair, (i,j), (j,k), and (k,i), is used as the net virial of the triplet. As a result, the formulation of atomic stress comes to the expression(6)σαβi=12Ωi∑j∂ϕpairi,j∂ri,jrαi,jrβi,jri,j+13Ωi∑j,k∂ϕanglei,j∂ri,jrαi,jrβi,jri,j+∂ϕanglej,k∂rj,krαj,krβj,krj,k+∂ϕanglek,i∂rk,irαk,irβk,irk,i,where the first term corresponds to pair interaction ϕpair(r) and the second term is for three-body interaction. In this expression, Ω(i), ϕpair(i,j), ϕangle(i,j), r(i,j), and rα(i,j)(α,β=1,2,3) are atomic volume, pairwise potential function, three-body potential function, interatomic distance (between i and j), and its vector components. In comparing all the components in (6) during actual MD simulation, it is understood that three-body term is relatively large so that it is not negligible. This expression of atomic stress including three-body term is also utilized in widely used MD software (e.g., LAMMPS) [19]. Besides, this estimation of atomic stress was found appropriate when we applied it to another MD study using the three-body-type Tersoff potential (for Si-Ge system) [20].


## 3. Results and Discussion

### 3.1. Optimization of Lattice Constant of Cu2S Crystal by First-Principle Calculation

The relation between volume and energy of Cu2S crystal unit is obtained as in Figure 3 by first-principle (FP) calculation including structural optimization. For the variety of the lattice constant ratio c/a, the lowest energy is obtained at the deviation Δc/a = 8%. Moreover, under that condition, volume expansion at +3% shows the smallest energy. From this optimized structure, therefore, lattice constants of the HC structure are a = 0.383 nm and c = 0.731 nm, for which c/a is 1.91. The length for c and the ratio c/a seem relatively larger than the HC structure of usual metals. It means that layers of sulfur (S) and copper (Cu) atoms are remarkably separated from each other in the Cu2S crystal. In the present paper, these optimized lattice constants are used as values to configure atomic positions in the specimen of Cu2S crystal.Figure 3
Volume-energy relation of Cu2S crystal structure.


### 3.2. Construction of Interatomic Potential Function for MD Calculations

It is experimentally and theoretically understood that interaction between Cu and S atoms is generally strong. So, the Cu-S interaction can be expressed by Morse potential function as shown by (3). The optimized configuration of Cu-S dimer is obtained by the FP calculation by us [21]. From this, we determine both the equilibrium distance parameters r0, βij and the energy parameter Dij needed in Morse potential. Figure 4 shows the relation between interatomic distance and its energy, which is obtained by the FP calculation. The obtained parameters for Morse potential is shown in Table 4.Table 4
Potential parameters fitted by FP calculations.
(a)
Morse potential (Cu-S)
Kind of pair
D
i
j
β
i
j
r
0
Unit
kJ/mole
1/nm
nm
Cu-S
298.916
1.500
2.029
(Mo-S)
(152.076)
(1.800)
(2.532)
(b)
Angular potential (S-Cu-S, Cu-S-Cu)
Combination
θ
0
Unit
deg.
S-Cu-S
98.966
Cu-S-Cu
98.966
(c)
B-M-H ionic potential (Cu-S)
Kind of pair
r
0
′
a
C
u
a
S
Unit
nm
nm
nm
2.5319
0.70052
1.8314∗
a
S
∗
has been already obtained for MoS2 potential [10].Figure 4
Distance-energy relation of Cu-S interaction obtained by FP calculation (Mo-S interaction is also shown for comparison).The optimized crystalline structure of Cu2S which has been obtained in the previous section can be used to adjust an angular-dependent three-body potential parameters shown in (4). In this process of fitting, we assume that Cu atoms are located at the averaged position between tetragonal and octahedral sites as shown in Figure 5, since there is possibility that Cu atoms may be located on both sites. In practice, we just need equilibrium three-body angle θ0 and pairwise distance r0′. Finally, they are θ0=98.966° and r0′ = 2.532 nm, respectively.Figure 5
Schematic of the way to average two types of coordination for Cu atoms in Cu2S crystal.Additionally, the BMH term requires ionic radiiaCu and aS for each element. The radius of Cu atom is supposed to be equivalent to the equilibrium length r0 in Morse potential above; the radius aS=1.834 nm for S is already available from the previous MD study of MoS2 [10].Thus, MD potential parameters to reproduce the HC structure of Cu2S are obtained by mostly FP calculation and they are summarized in Table 4.


### 3.3. Sliding and Friction Behavior of Cu2S

#### 3.3.1. Dependency on Crystalline Orientation

The dependence on crystalline orientation in sliding behavior of Cu2S is as follows. Figure 6 shows the atomic configurations at t = 600 ps with sliding velocity v = 0.5 m/s, comparing between models which are sliding in (a)-direction (01¯10) sliding and in (b)-direction (21¯1¯0). Both models once collapse and lose their crystalline stacking during sliding. But, after that, the atomic structure recovers its original stacking during long-time sliding. The difference between two models is not found just visually from atomic configurations, as seen in the broken black- and yellow-colored rectangular areas depicted in the figures. Therefore, the RDF analysis will be helpful for identifying the change of crystalline structures. Besides, in order to recognize transition of an atomic force during sliding, the analysis of shear stress will be helpful.Figure 6
Comparison of instantaneous atomic configurations obtained att = 600 ps with v = 0.5 m/s (the picture is drawn onto y-z plane).The result of RDF analysis is shown in Figure7, which is for sliding velocity v = 0.5 m/s and compares results between (a)-direction and (b)-direction models. In these RDF figures, peaks are found at some distances, by which the nature of crystalline structure is confirmed. In particular, when we take a look at peaks found in the range of 0.2~0.8 nm for two models, the steepness of the distribution in (b)-direction is totally stronger than that in (a)-direction. So, the structure obtained by the sliding in (b)-direction tends to retain more local crystalline structure than that obtained in (a)-direction sliding.Figure 7
Time transition of RDF forv = 0.5 m/s (comparison between (a)- and (b)-direction sliding).Just when the resulting sliding distance becomes 1.0 nm, is the RDF obtained as shown in Figure8 for different sliding velocities (v = 1.0 and 5.0 m/s). At the same sliding distance, the crystalline stacking is the same, and so the results may be compared as for different velocity conditions. One peak of RDF is just identified for the 1st neighbor distance, which means the crystalline structure has been lost and the atomic movement shows somewhat of fluidity.Figure 8
The comparison of RDF at the same sliding distance (1.0 nm) between (a)- and (b)-direction sliding.Shear stress averaged over the whole specimen,τtot, represents the resistance to the sliding. The time transition of τtot reflects atomic force state in the sliding. During relaxation process, in fact, a certain value of shear stress already occurs. This is due to atomic rearrangement inside the crystal unit. After the specimen is fully relaxed, we reset the stress value. Let the value at the beginning of sliding be τtot(0). Then, the absolute value of deviation of current value τtot from initial value τtot(0) is defined to be (7)τa=τtot-τtot0.This stress deviation τa reflects a sudden increase of friction due to roughness of slide planes. The time transition of τa is obtained as shown in Figure 9. Figure 9(A) shows the relation between the time and the stress deviation τa. τa is averaged over the whole specimen. Those graphs show that τa for (a)-direction sliding is larger than that in (b)-direction. The downward arrows in the figure display each maximum value. At low velocity, atoms near contact layers feel relatively larger resistance to slide in (a)-direction sliding, while the (b)-direction sliding seems much easier. However, for the cases with higher velocity such as v = 1.0 m/s or 5.0 m/s, the tendency of stress τa is changed as shown in Figures 9(B) and 9(C). As shown in Figure 9(B), for example, for v = 1.0 m/s, the peak value in (b)-direction is sometimes larger than that in (a)-direction, especially during 50–100 ps or 250–300 ps.Figure 9
Comparison of the absolute value of shear stress deviationτa between (a)- and (b)-direction.We summarize the dependency of sliding direction as follows. When the shear velocity is small, (b)-direction (01¯10) sliding is easier than (a)-direction (01¯10). However, once the sliding has completed, such as to 100 or 200 ps from the start of sliding, distinction between sliding directions becomes negligible. The cause of this behavior would be guessed and explained as follows. When the sliding is carried out sufficiently, the crystalline structure of Cu2S has collapsed and the crystalline slip system (plane and direction) no longer determines the sliding behavior. Even if there exists a certain easy-glide plane or direction, the crystal does not always show especially smaller shear stress in sliding.


#### 3.3.2. Dependency on Sliding Velocity

Figure10(A) shows the relation between shear distance and shear stress deviation τa. It compares results for a variety of velocities ranging from 0.5 to 5.0 m/s. For all velocities, and for the sliding distance up to 0.2 nm, τa linearly increases with the same gradient. This fact means that the material deforms in elastic manner for small sliding and deformation. However, after that, τa decreases suddenly. The larger the sliding velocity is, the larger the averaged value of τa is obtained. It is reasonable that time rate of decrease is smaller for higher speed sliding. As the whole time calculated here is seen, the case of v = 5.0 m/s exhibits the largest level of τa. For (b)-direction, Figure 10(B) shows that the value of τa increases with increase of sliding velocity. It is concluded that shear stress of this material depends largely on sliding velocity.Figure 10
Relation between shear distance and the deviation of shear stressτa: comparison among different sliding velocities.The relation betweenτa and sliding velocity v is summarized in Table 5. τave, time average of τa, is also shown at the bottom line. It indicates that τave of (a)-direction is approximately 10% smaller than that of (b)-direction. While the easiest glide direction/plane of the HC structure is 01¯10(0001) (i.e., (b)-direction), it does not show small resistance in the sliding. Thus, it is concluded that Cu2S crystal does surely show sliding motion, but it is not by usual crystalline slip (glide) mechanism of ordinary metals.Table 5
Relation between averaged shear stressτave and sliding velocity v.
Sliding velocityv m/s
Averaged shear stressτave GPa
(A) (a)-direction sliding
(B) (b)-direction sliding
0.5
19.3
22.8
1.0
21.4
25.4
5.0
26.4
28.3
Average
22.4
25.5On the other hand, the fact that shear stress increases with increase of velocity is interesting. It means that this material deforms by sliding with somewhat liquid-like (amorphous) behavior. The deformation rate (i.e., velocity) in such amorphous substance produces shear stress inside due to its viscosity, and so the shear stress depends on sliding (shear) velocity gradient. As shown in Figure11, averaged shear stresses τave summarized in Table 5 are replotted in terms of sliding velocity. As recognized in Figure 11, shear stress is almost proportional to the sliding velocity. In the present simulation condition, system temperature is kept at 10 kelvins, which is quite low. But the atomic mobility at contact layers between materials should rise due to the work done by sliding motion. In observing atomic motion, they show intermittent slip, which is sometimes called “stick and slip motion” in nanotribological consideration. Furthermore, when in Figure 11 we extrapolate approximated straight lines down to v=0, shear stress τave does not return to zero, but it remains at a finite value τ0 = 20~25 GPa. These ideal shear stresses at v=0 correspond to static friction of materials. It also means that τ0 are at least required for layered atoms to exhibit a material flow.Figure 11
Averaged shear stress replotted against sliding velocity ((a)- and (b)-directions (21¯1¯0 and 01¯10)).


## 3.1. Optimization of Lattice Constant of Cu2S Crystal by First-Principle Calculation

The relation between volume and energy of Cu2S crystal unit is obtained as in Figure 3 by first-principle (FP) calculation including structural optimization. For the variety of the lattice constant ratio c/a, the lowest energy is obtained at the deviation Δc/a = 8%. Moreover, under that condition, volume expansion at +3% shows the smallest energy. From this optimized structure, therefore, lattice constants of the HC structure are a = 0.383 nm and c = 0.731 nm, for which c/a is 1.91. The length for c and the ratio c/a seem relatively larger than the HC structure of usual metals. It means that layers of sulfur (S) and copper (Cu) atoms are remarkably separated from each other in the Cu2S crystal. In the present paper, these optimized lattice constants are used as values to configure atomic positions in the specimen of Cu2S crystal.Figure 3
Volume-energy relation of Cu2S crystal structure.


## 3.2. Construction of Interatomic Potential Function for MD Calculations

It is experimentally and theoretically understood that interaction between Cu and S atoms is generally strong. So, the Cu-S interaction can be expressed by Morse potential function as shown by (3). The optimized configuration of Cu-S dimer is obtained by the FP calculation by us [21]. From this, we determine both the equilibrium distance parameters r0, βij and the energy parameter Dij needed in Morse potential. Figure 4 shows the relation between interatomic distance and its energy, which is obtained by the FP calculation. The obtained parameters for Morse potential is shown in Table 4.Table 4
Potential parameters fitted by FP calculations.
(a)
Morse potential (Cu-S)
Kind of pair
D
i
j
β
i
j
r
0
Unit
kJ/mole
1/nm
nm
Cu-S
298.916
1.500
2.029
(Mo-S)
(152.076)
(1.800)
(2.532)
(b)
Angular potential (S-Cu-S, Cu-S-Cu)
Combination
θ
0
Unit
deg.
S-Cu-S
98.966
Cu-S-Cu
98.966
(c)
B-M-H ionic potential (Cu-S)
Kind of pair
r
0
′
a
C
u
a
S
Unit
nm
nm
nm
2.5319
0.70052
1.8314∗
a
S
∗
has been already obtained for MoS2 potential [10].Figure 4
Distance-energy relation of Cu-S interaction obtained by FP calculation (Mo-S interaction is also shown for comparison).The optimized crystalline structure of Cu2S which has been obtained in the previous section can be used to adjust an angular-dependent three-body potential parameters shown in (4). In this process of fitting, we assume that Cu atoms are located at the averaged position between tetragonal and octahedral sites as shown in Figure 5, since there is possibility that Cu atoms may be located on both sites. In practice, we just need equilibrium three-body angle θ0 and pairwise distance r0′. Finally, they are θ0=98.966° and r0′ = 2.532 nm, respectively.Figure 5
Schematic of the way to average two types of coordination for Cu atoms in Cu2S crystal.Additionally, the BMH term requires ionic radiiaCu and aS for each element. The radius of Cu atom is supposed to be equivalent to the equilibrium length r0 in Morse potential above; the radius aS=1.834 nm for S is already available from the previous MD study of MoS2 [10].Thus, MD potential parameters to reproduce the HC structure of Cu2S are obtained by mostly FP calculation and they are summarized in Table 4.


## 3.3. Sliding and Friction Behavior of Cu2S

### 3.3.1. Dependency on Crystalline Orientation

The dependence on crystalline orientation in sliding behavior of Cu2S is as follows. Figure 6 shows the atomic configurations at t = 600 ps with sliding velocity v = 0.5 m/s, comparing between models which are sliding in (a)-direction (01¯10) sliding and in (b)-direction (21¯1¯0). Both models once collapse and lose their crystalline stacking during sliding. But, after that, the atomic structure recovers its original stacking during long-time sliding. The difference between two models is not found just visually from atomic configurations, as seen in the broken black- and yellow-colored rectangular areas depicted in the figures. Therefore, the RDF analysis will be helpful for identifying the change of crystalline structures. Besides, in order to recognize transition of an atomic force during sliding, the analysis of shear stress will be helpful.Figure 6
Comparison of instantaneous atomic configurations obtained att = 600 ps with v = 0.5 m/s (the picture is drawn onto y-z plane).The result of RDF analysis is shown in Figure7, which is for sliding velocity v = 0.5 m/s and compares results between (a)-direction and (b)-direction models. In these RDF figures, peaks are found at some distances, by which the nature of crystalline structure is confirmed. In particular, when we take a look at peaks found in the range of 0.2~0.8 nm for two models, the steepness of the distribution in (b)-direction is totally stronger than that in (a)-direction. So, the structure obtained by the sliding in (b)-direction tends to retain more local crystalline structure than that obtained in (a)-direction sliding.Figure 7
Time transition of RDF forv = 0.5 m/s (comparison between (a)- and (b)-direction sliding).Just when the resulting sliding distance becomes 1.0 nm, is the RDF obtained as shown in Figure8 for different sliding velocities (v = 1.0 and 5.0 m/s). At the same sliding distance, the crystalline stacking is the same, and so the results may be compared as for different velocity conditions. One peak of RDF is just identified for the 1st neighbor distance, which means the crystalline structure has been lost and the atomic movement shows somewhat of fluidity.Figure 8
The comparison of RDF at the same sliding distance (1.0 nm) between (a)- and (b)-direction sliding.Shear stress averaged over the whole specimen,τtot, represents the resistance to the sliding. The time transition of τtot reflects atomic force state in the sliding. During relaxation process, in fact, a certain value of shear stress already occurs. This is due to atomic rearrangement inside the crystal unit. After the specimen is fully relaxed, we reset the stress value. Let the value at the beginning of sliding be τtot(0). Then, the absolute value of deviation of current value τtot from initial value τtot(0) is defined to be (7)τa=τtot-τtot0.This stress deviation τa reflects a sudden increase of friction due to roughness of slide planes. The time transition of τa is obtained as shown in Figure 9. Figure 9(A) shows the relation between the time and the stress deviation τa. τa is averaged over the whole specimen. Those graphs show that τa for (a)-direction sliding is larger than that in (b)-direction. The downward arrows in the figure display each maximum value. At low velocity, atoms near contact layers feel relatively larger resistance to slide in (a)-direction sliding, while the (b)-direction sliding seems much easier. However, for the cases with higher velocity such as v = 1.0 m/s or 5.0 m/s, the tendency of stress τa is changed as shown in Figures 9(B) and 9(C). As shown in Figure 9(B), for example, for v = 1.0 m/s, the peak value in (b)-direction is sometimes larger than that in (a)-direction, especially during 50–100 ps or 250–300 ps.Figure 9
Comparison of the absolute value of shear stress deviationτa between (a)- and (b)-direction.We summarize the dependency of sliding direction as follows. When the shear velocity is small, (b)-direction (01¯10) sliding is easier than (a)-direction (01¯10). However, once the sliding has completed, such as to 100 or 200 ps from the start of sliding, distinction between sliding directions becomes negligible. The cause of this behavior would be guessed and explained as follows. When the sliding is carried out sufficiently, the crystalline structure of Cu2S has collapsed and the crystalline slip system (plane and direction) no longer determines the sliding behavior. Even if there exists a certain easy-glide plane or direction, the crystal does not always show especially smaller shear stress in sliding.


### 3.3.2. Dependency on Sliding Velocity

Figure10(A) shows the relation between shear distance and shear stress deviation τa. It compares results for a variety of velocities ranging from 0.5 to 5.0 m/s. For all velocities, and for the sliding distance up to 0.2 nm, τa linearly increases with the same gradient. This fact means that the material deforms in elastic manner for small sliding and deformation. However, after that, τa decreases suddenly. The larger the sliding velocity is, the larger the averaged value of τa is obtained. It is reasonable that time rate of decrease is smaller for higher speed sliding. As the whole time calculated here is seen, the case of v = 5.0 m/s exhibits the largest level of τa. For (b)-direction, Figure 10(B) shows that the value of τa increases with increase of sliding velocity. It is concluded that shear stress of this material depends largely on sliding velocity.Figure 10
Relation between shear distance and the deviation of shear stressτa: comparison among different sliding velocities.The relation betweenτa and sliding velocity v is summarized in Table 5. τave, time average of τa, is also shown at the bottom line. It indicates that τave of (a)-direction is approximately 10% smaller than that of (b)-direction. While the easiest glide direction/plane of the HC structure is 01¯10(0001) (i.e., (b)-direction), it does not show small resistance in the sliding. Thus, it is concluded that Cu2S crystal does surely show sliding motion, but it is not by usual crystalline slip (glide) mechanism of ordinary metals.Table 5
Relation between averaged shear stressτave and sliding velocity v.
Sliding velocityv m/s
Averaged shear stressτave GPa
(A) (a)-direction sliding
(B) (b)-direction sliding
0.5
19.3
22.8
1.0
21.4
25.4
5.0
26.4
28.3
Average
22.4
25.5On the other hand, the fact that shear stress increases with increase of velocity is interesting. It means that this material deforms by sliding with somewhat liquid-like (amorphous) behavior. The deformation rate (i.e., velocity) in such amorphous substance produces shear stress inside due to its viscosity, and so the shear stress depends on sliding (shear) velocity gradient. As shown in Figure11, averaged shear stresses τave summarized in Table 5 are replotted in terms of sliding velocity. As recognized in Figure 11, shear stress is almost proportional to the sliding velocity. In the present simulation condition, system temperature is kept at 10 kelvins, which is quite low. But the atomic mobility at contact layers between materials should rise due to the work done by sliding motion. In observing atomic motion, they show intermittent slip, which is sometimes called “stick and slip motion” in nanotribological consideration. Furthermore, when in Figure 11 we extrapolate approximated straight lines down to v=0, shear stress τave does not return to zero, but it remains at a finite value τ0 = 20~25 GPa. These ideal shear stresses at v=0 correspond to static friction of materials. It also means that τ0 are at least required for layered atoms to exhibit a material flow.Figure 11
Averaged shear stress replotted against sliding velocity ((a)- and (b)-directions (21¯1¯0 and 01¯10)).


## 3.3.1. Dependency on Crystalline Orientation

The dependence on crystalline orientation in sliding behavior of Cu2S is as follows. Figure 6 shows the atomic configurations at t = 600 ps with sliding velocity v = 0.5 m/s, comparing between models which are sliding in (a)-direction (01¯10) sliding and in (b)-direction (21¯1¯0). Both models once collapse and lose their crystalline stacking during sliding. But, after that, the atomic structure recovers its original stacking during long-time sliding. The difference between two models is not found just visually from atomic configurations, as seen in the broken black- and yellow-colored rectangular areas depicted in the figures. Therefore, the RDF analysis will be helpful for identifying the change of crystalline structures. Besides, in order to recognize transition of an atomic force during sliding, the analysis of shear stress will be helpful.Figure 6
Comparison of instantaneous atomic configurations obtained att = 600 ps with v = 0.5 m/s (the picture is drawn onto y-z plane).The result of RDF analysis is shown in Figure7, which is for sliding velocity v = 0.5 m/s and compares results between (a)-direction and (b)-direction models. In these RDF figures, peaks are found at some distances, by which the nature of crystalline structure is confirmed. In particular, when we take a look at peaks found in the range of 0.2~0.8 nm for two models, the steepness of the distribution in (b)-direction is totally stronger than that in (a)-direction. So, the structure obtained by the sliding in (b)-direction tends to retain more local crystalline structure than that obtained in (a)-direction sliding.Figure 7
Time transition of RDF forv = 0.5 m/s (comparison between (a)- and (b)-direction sliding).Just when the resulting sliding distance becomes 1.0 nm, is the RDF obtained as shown in Figure8 for different sliding velocities (v = 1.0 and 5.0 m/s). At the same sliding distance, the crystalline stacking is the same, and so the results may be compared as for different velocity conditions. One peak of RDF is just identified for the 1st neighbor distance, which means the crystalline structure has been lost and the atomic movement shows somewhat of fluidity.Figure 8
The comparison of RDF at the same sliding distance (1.0 nm) between (a)- and (b)-direction sliding.Shear stress averaged over the whole specimen,τtot, represents the resistance to the sliding. The time transition of τtot reflects atomic force state in the sliding. During relaxation process, in fact, a certain value of shear stress already occurs. This is due to atomic rearrangement inside the crystal unit. After the specimen is fully relaxed, we reset the stress value. Let the value at the beginning of sliding be τtot(0). Then, the absolute value of deviation of current value τtot from initial value τtot(0) is defined to be (7)τa=τtot-τtot0.This stress deviation τa reflects a sudden increase of friction due to roughness of slide planes. The time transition of τa is obtained as shown in Figure 9. Figure 9(A) shows the relation between the time and the stress deviation τa. τa is averaged over the whole specimen. Those graphs show that τa for (a)-direction sliding is larger than that in (b)-direction. The downward arrows in the figure display each maximum value. At low velocity, atoms near contact layers feel relatively larger resistance to slide in (a)-direction sliding, while the (b)-direction sliding seems much easier. However, for the cases with higher velocity such as v = 1.0 m/s or 5.0 m/s, the tendency of stress τa is changed as shown in Figures 9(B) and 9(C). As shown in Figure 9(B), for example, for v = 1.0 m/s, the peak value in (b)-direction is sometimes larger than that in (a)-direction, especially during 50–100 ps or 250–300 ps.Figure 9
Comparison of the absolute value of shear stress deviationτa between (a)- and (b)-direction.We summarize the dependency of sliding direction as follows. When the shear velocity is small, (b)-direction (01¯10) sliding is easier than (a)-direction (01¯10). However, once the sliding has completed, such as to 100 or 200 ps from the start of sliding, distinction between sliding directions becomes negligible. The cause of this behavior would be guessed and explained as follows. When the sliding is carried out sufficiently, the crystalline structure of Cu2S has collapsed and the crystalline slip system (plane and direction) no longer determines the sliding behavior. Even if there exists a certain easy-glide plane or direction, the crystal does not always show especially smaller shear stress in sliding.


## 3.3.2. Dependency on Sliding Velocity

Figure10(A) shows the relation between shear distance and shear stress deviation τa. It compares results for a variety of velocities ranging from 0.5 to 5.0 m/s. For all velocities, and for the sliding distance up to 0.2 nm, τa linearly increases with the same gradient. This fact means that the material deforms in elastic manner for small sliding and deformation. However, after that, τa decreases suddenly. The larger the sliding velocity is, the larger the averaged value of τa is obtained. It is reasonable that time rate of decrease is smaller for higher speed sliding. As the whole time calculated here is seen, the case of v = 5.0 m/s exhibits the largest level of τa. For (b)-direction, Figure 10(B) shows that the value of τa increases with increase of sliding velocity. It is concluded that shear stress of this material depends largely on sliding velocity.Figure 10
Relation between shear distance and the deviation of shear stressτa: comparison among different sliding velocities.The relation betweenτa and sliding velocity v is summarized in Table 5. τave, time average of τa, is also shown at the bottom line. It indicates that τave of (a)-direction is approximately 10% smaller than that of (b)-direction. While the easiest glide direction/plane of the HC structure is 01¯10(0001) (i.e., (b)-direction), it does not show small resistance in the sliding. Thus, it is concluded that Cu2S crystal does surely show sliding motion, but it is not by usual crystalline slip (glide) mechanism of ordinary metals.Table 5
Relation between averaged shear stressτave and sliding velocity v.
Sliding velocityv m/s
Averaged shear stressτave GPa
(A) (a)-direction sliding
(B) (b)-direction sliding
0.5
19.3
22.8
1.0
21.4
25.4
5.0
26.4
28.3
Average
22.4
25.5On the other hand, the fact that shear stress increases with increase of velocity is interesting. It means that this material deforms by sliding with somewhat liquid-like (amorphous) behavior. The deformation rate (i.e., velocity) in such amorphous substance produces shear stress inside due to its viscosity, and so the shear stress depends on sliding (shear) velocity gradient. As shown in Figure11, averaged shear stresses τave summarized in Table 5 are replotted in terms of sliding velocity. As recognized in Figure 11, shear stress is almost proportional to the sliding velocity. In the present simulation condition, system temperature is kept at 10 kelvins, which is quite low. But the atomic mobility at contact layers between materials should rise due to the work done by sliding motion. In observing atomic motion, they show intermittent slip, which is sometimes called “stick and slip motion” in nanotribological consideration. Furthermore, when in Figure 11 we extrapolate approximated straight lines down to v=0, shear stress τave does not return to zero, but it remains at a finite value τ0 = 20~25 GPa. These ideal shear stresses at v=0 correspond to static friction of materials. It also means that τ0 are at least required for layered atoms to exhibit a material flow.Figure 11
Averaged shear stress replotted against sliding velocity ((a)- and (b)-directions (21¯1¯0 and 01¯10)).


## 4. Conclusion

In this paper, we perform a molecular dynamics (MD) study on copper sulfide material (Cu2S) under the sliding motion. The MD simulation includes a new implementation of potential energy function for Cu2S crystal, where some first-principle calculations are utilized. The potential energy function has three-body term. The following results are obtained and we will make conclusion:(1)
First-principle calculations are performed for the hexagonal unit structure of Cu2S crystal. From the optimized structure, original lattice parameters a and c of unit cell and stable configuration of S and Cu atoms are obtained. Based on those results, an interatomic potential function including ionic and Morse terms as well as angular-dependent term is constructed for Cu and S system. It is confirmed that this potential function can show enough stability for the hexagonal crystal structure of Cu2S in MD simulations.(2)
An easy-glide direction21¯1¯0 on (0001) plane generally found for usual hexagonal crystals is not the direction with small resistance to sliding.(3)
Cu2S crystal shows partially liquid-like structure and behavior, in which the shear stress occurring in the material depends on the sliding velocity.In the further study, we should clarify the temperature dependence and dependency on compressive condition during sliding of Cu2S. It is supposed that they are made possible by analyzing more atomic trajectories obtained by MD simulations. It will lead to total understanding of sliding mechanism of Cu2S crystal.



---

*Source: 963257-2015-10-26.xml*