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

To clarify the nanosized mechanism of good lubrication in copper disulfide (Cu 2 S) crystal which is used as a sliding material, atomistic modeling of Cu 2 S 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 Cu 2 S are carefully estimated by firstprinciple calculations. Then, approximating these interactions, we originally construct a conventional interatomic potential function of Cu 2 S crystal in its hexagonal phase. By using this potential function, we performMD simulation of Cu 2 S crystal which is subjected to shear loading parallel to the basal plane. We compare results obtained by different conditions of sliding directions. Unlike ordinary hexagonalmetallic crystals, it is found that the easy-glide direction does not always show small shear stress for Cu 2 S 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 Cu 2 S crystal exhibits somewhat liquid-like (amorphous) behavior in sliding condition and shear resistance increases with increase of sliding speed.


Introduction
In recent years, using lead (Pb) in metallic material is strictly prohibited due to risk to human health.That is why leadfree 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 (brakepads) 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, MoS 2 , 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 (Cu 2 S) and they are now proceeding to its industrial use [3].Akamatsu et al. are also exploring the fabricating method of Cu 2 S alloy by the motivation of avoiding the risk to human health [4].

Journal of Materials
In the field of mineralogy, Cu 2 S crystal is well known as a mineral called chalcocite.Recently, since the Cu 2 S 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 Cu 2 S 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 Cu 2 S (Cu : S = 2 : 1).Interestingly, this Cu 2 S compound exhibits very complicated solid phase transition behavior [8,9].Around the standard temperature, the unit structure of Cu 2 S 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 Cu 2 S 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 MoS 2 system which is famous as good solid lubricant, where the S atom makes firm bonding with metal atoms [10].
Other than chalcocite (Cu 2 S), there are some well-known natural minerals based on copper alloy [11,12], such as bornite (Cu 5 FeS 4 ), chalcopyrite (CuFeS 2 ), covellite (CuS), and digenite (Cu 9 S 5 ).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, Cu 2 S crystal is still a main compound, so it should be focused on particularly.Thus, in this paper, we will study the atomistic structure of Cu 2 S 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 Cu 2 S.So far, there have been some studies of copper sulfides by ab 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 Cu 2 S 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,  and , assuming the HC configuration.
In order to understand the dynamics of sliding or friction occurring in Cu 2 S crystal, molecular dynamics (MD) simulation is suitable.Lattice constants  and  and groundstate energies which are obtained by ab initio calculations are used in the determination of interatomic potential for MD simulation.In the MD simulation, an existing interatomic potential function of MoS 2 [10] is referred to for its function or functional form.This seems reasonable because those two crystals of metal sulfide (Cu 2 S and MoS 2 ) have many common features, as described above.
Thus, we simulate atomistically the sliding of Cu 2 S 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 Cu 2 S and show some theories needed for the modeling and the calculation.Then, the essence of the ab 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 in ab initio calculation and MD results of sliding simulation are shown and their discussion is made.Finally, the conclusion of this paper will be made.

Theory and Method
2.1.Atomic Structure of Cu 2 S Crystal.As pure compounds of copper sulfide, there are covellite (CuS) and chalcocite (Cu 2 S), and the latter is focused on now.For Cu 2 S, 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 Cu 2 S material here is assumed to be in HC.HC is also taken in MoS 2 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 Cu 2 S and MoS 2 , will be good at lubrication.

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 Cu 2 S 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 Cu 2 S 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 Cu 2 S 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 Cu 2 S in HC is shown in Table 1. Figure 1 shows initial position of atoms in Cu 2 S crystal unit.
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 (, ) are changed as follows.Keeping a certain ratio of /, the lengths (, ) are varied from compressive state to tensile one.The deviation of / from the initial value is +4∼+10%.The volumetric change by the change of (, ) 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 −18 J = 13.60 eV.Table 2 summarizes calculation conditions and parameters used in the present study.The number of -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 = 46552 -points are large enough and this will be used through the present study.

The Potential Function of Cu 2 S for MD.
As mentioned above, copper disulfide, Cu 2 S, 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, MoS 2 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 MoS 2 [10] can be utilized as the reference to construct the function needed in the present study for Cu 2 S. Though MoS 2 and Cu 2 S 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 threebody interaction is crucial for stabilization of both crystals.To begin with, we apply the potential forms used in the MoS 2 study, and then our results obtained by first-principle calculation are properly fitted to a new potential function for Cu 2 S. 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 Cu 2 S 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 threebody potential.Parameters for those potential functions are summarized in Table 3: In those equations,   or   is electronic charge of each atom, and the values  S = −2.0 and  Cu = +1.0 are employed for S and Cu atoms, respectively. 0 = 1 [kcal/ Å] = 4.186 [kJ/ Å] is the parameter which determines the stiffness of each ionic sphere and is transferred from the study of MoS 2 .In the Morse potential, the equilibrium distance  0 and bond energy   are fitted to those in the optimized structure which is obtained in our ab initio calculation.In the angular potential, equilibrium three-body angle  0 is also determined by our ab initio calculation.But, it is difficult to derive the spring constant   just from our ab initio calculation, so the same value as in MoS 2 is applied here.Of course, the value of   influences rigidity.But its function form is basically harmonic and we found that the effect of spring constant on sliding is not so crucial.

Molecular Dynamics
Model for Cu 2 S Sliding. Figure 2 shows the MD model for sliding simulation.In the first place, the optimized atomic configuration obtained by ab initio calculation is prepared and then is thermally equilibrated at a finite temperature.Then, the structure is compressed in  direction as shown in Figure 2.After that, the atomic velocity in the  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 Cu 2 S 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.It is supposed that crystalline slip plane of Cu 2 S plays an important role in the slip deformation.As stated above, Cu 2 S 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 10⟩.We will call it (b)-direction hereafter.The second easiest slip direction should be ⟨0110⟩, 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 Cu 2 S 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.

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.

Crystal Structure Analysis by Radial Distribution Function (RDF).
Usually, solid crystal has a three-dimensional order.The Cu 2 S crystal used in the present study exhibits HC structure.One of the methods to identify the crystalline order is radial distribution function (RDF) () [18].Briefly speaking, the function () 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.

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, atomicscale 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 ().When any selected atom  conducts homogeneous deformation with regard to the relative position to another neighbor atom , strain at the interatomic space between  and  can be presented by their interatomic distance (, ) and interatomic vector r(, ) (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  is expressed as where the summation is composed of the work done by atomic motion, which is called "virial."Ω() in the denominator is the volume supposedly occupied by the atom , so ( 5) is regarded as the energy density around that atom.
In the present study, Ω() 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, , , and  (, , ).In practice, the combination of virials of each pair, (, ), (, ), and (, ), is used as the net virial of the triplet.As a result, the formulation of atomic stress comes to the expression where the first term corresponds to pair interaction  pair () and the second term is for three-body interaction.In this expression, Ω(),  pair (, ),  angle (, ), (, ), and   (, ) (,  = 1, 2, 3) are atomic volume, pairwise potential function, three-body potential function, interatomic distance (between  and ), 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].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  0 ,   and the energy parameter   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.

Results and Discussion
The optimized crystalline structure of Cu 2 S 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   0 .Finally, they are  0 = 98.966 ∘ and   0 = 2.532 nm, respectively.
Additionally, the BMH term requires ionic radii  Cu and  S for each element.The radius of Cu atom is supposed to be equivalent to the equilibrium length  0 in Morse potential above; the radius  S = 1.834 nm for S is already available from the previous MD study of MoS 2 [10].
Thus, MD potential parameters to reproduce the HC structure of Cu 2 S are obtained by mostly FP calculation and they are summarized in Table 4.

Sliding and Friction
Behavior of Cu 2 S 3.3.1.Dependency on Crystalline Orientation.The dependence on crystalline orientation in sliding behavior of Cu 2 S is as follows.Figure 6 shows the atomic configurations at  = 600 ps with sliding velocity V = 0.5 m/s, comparing between models which are sliding in (a)-direction (⟨0110⟩) sliding and in (b)-direction (⟨21 10⟩).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.
The result of RDF analysis is shown in Figure 7, 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.8nm 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.
Just when the resulting sliding distance becomes 1.0 nm, is the RDF obtained as shown in Figure 8 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.
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 This stress deviation   reflects a sudden increase of friction due to roughness of slide planes.The time transition of   is obtained as shown in Figure 9.We summarize the dependency of sliding direction as follows.When the shear velocity is small, (b)-direction (⟨0110⟩) sliding is easier than (a)-direction (⟨0110⟩).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 Cu 2 S 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.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,   linearly increases with the same gradient.This fact means that the material deforms in elastic manner for small sliding and deformation.However, after that,   decreases suddenly.The larger the sliding velocity is, the larger the averaged value of   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   .For (b)-direction, Figure 10(B) shows that the value of   increases with increase of sliding velocity.It is concluded that shear stress of this material depends largely on sliding velocity.The relation between   and sliding velocity V is summarized in Table 5.  ave , time average of   , 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 ⟨0110⟩(0001) (i.e., (b)-direction), it does not show small resistance in the sliding.Thus, it is concluded that Cu 2 S crystal does surely show sliding motion, but it is not by usual crystalline slip (glide) mechanism of ordinary metals.

Dependency on Sliding Velocity.
On 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 Figure 11, 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.

Conclusion
In this paper, we perform a molecular dynamics (MD) study on copper sulfide material (Cu 2 S) under the sliding

Figure 2 :
Figure 2: Abstract of MD sliding model of Cu 2 S crystal.

Figure 3 :
Figure 3: Volume-energy relation of Cu 2 S crystal structure.

Figure 4 :
Figure 4: Distance-energy relation of Cu-S interaction obtained by FP calculation (Mo-S interaction is also shown for comparison).

Figure 5 :
Figure 5: Schematic of the way to average two types of coordination for Cu atoms in Cu 2 S crystal.

Figure 6 :
Figure 6: Comparison of instantaneous atomic configurations obtained at  = 600 ps with V = 0.5 m/s (the picture is drawn onto - plane).

Figure 9 (
A) shows the relation between the time and the stress deviation   .  is averaged over the whole specimen.Those graphs show that   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   is changed as shown in Figures9(B) and 9(C).As shown in Figure9(B), for 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 10 (
A) shows the relation between shear distance and shear stress deviation   .

Figure 8 :
Figure 8: The comparison of RDF at the same sliding distance (1.0 nm) between (a)-and (b)-direction sliding.

Table 5 :
Relation between averaged shear stress  ave and sliding velocity V.Sliding velocity V m/sAveraged shear stress  ave GPa (A) (a)-direction sliding (B) (b)

Table 1 :
Cu 2 S atomic coordinates and crystal data.

Table 3 :
Interatomic potential parameters for Cu 2 S crystal.
The relation between volume and energy of Cu 2 S crystal unit is obtained as in Figure3by firstprinciple (FP) calculation including structural optimization.For the variety of the lattice constant ratio /, the lowest energy is obtained at the deviation Δ/ = 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  = 0.383 nm and  = 0.731 nm, for which / is 1.91.The length for  and the ratio / 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 Cu 2 S crystal.In the present paper, these optimized lattice constants are used as values to configure atomic positions in the specimen of Cu 2 S crystal.
3.1.Optimization of Lattice Constant of Cu 2 S Crystal byFirst-Principle Calculation.3.2.Construction of Interatomic Potential Function for MDCalculations.It is experimentally and theoretically understood that interaction between Cu and S atoms is generally

Table 4 :
Potential parameters fitted by FP calculations.