JP5308285B2 - Simulation system and simulation program - Google Patents

Abstract

<P>PROBLEM TO BE SOLVED: To efficiently simulate the behavior of a layered compound. <P>SOLUTION: Each layer of a layered compound is coarse-grained as a rigid sheet being a dynamic unit to determine internal action caused by the interlayer interaction with a rigid sheet of an adjacent layer and internal oscillation (S1-S3) with respect to the rigid sheet of each coarse-grained layer. After this, the behavior of a layered compound formed of a plurality of rigid sheets is computed (S5). <P>COPYRIGHT: (C)2011,JPO&amp;INPIT

Description

  The present invention relates to a simulation system that simulates the behavior of a layered compound such as graphite and mica.

  Certain materials in layered compounds such as graphite and mica may exhibit low friction properties. Mica is not very low friction, but solids such as graphite and molybdenum disulfide have a low coefficient of friction. These compounds are used as solid lubricants in parts of satellites that cannot use oil, such as in outer space, or the compounds that form crystals are mixed into automobile engine oils as additives and crystallized on the sliding surface. And low friction.

  The reason why the frictional force is low in the layered compound is that the interlayer interaction is a crystal like graphite due to van der Waals force or an ionic crystal like molybdenum disulfide because the slip occurs between the layers. It is explained that there is.

  If such a layered compound can be analyzed in detail, structural design such as what layered compound should be used according to the application can be performed.

  In this connection, Patent Document 1 discloses a Brownian dynamics simulation in which a macromolecule (string micelle) in a fluid is coarse-grained. Moreover, the macroscopic coarse-graining technique is disclosed in Patent Documents 2 and 3 and the like. Furthermore, there exists patent document 4 etc. about the molecular simulation regarding a tribo.

Japanese Patent No. 03980600 JP 2002-80406 A JP 2006-338449 A JP-A-8-35831

  In Patent Documents 1 to 3, since an interaction having a two-dimensional periodicity is not defined as an interlayer interaction, the motion of the layered compound cannot be practically handled. In addition, because it is a method for liquids, it handles systems that are isotropic due to the internal vibrations of each coarse-grained particle, such as the diffusion coefficient, and reproduces the characteristics of layered compounds with high anisotropic effects. Can not. Therefore, the method described in Patent Document 1 cannot analyze the behavior of the layered compound.

  In Patent Document 4, in principle, the friction behavior of a layered compound can be handled by an all-atom model. However, when trying to analyze the realistic behavior of layered compounds, the number of symmetric atoms is very large. Therefore, if an all-atom model is used, enormous calculation time is required.

  The present invention coarse-grains each layer of a layered compound as a rigid body sheet as a mechanical unit, and the coarse-grained rigid sheet of each layer exhibits an internal interaction due to interlayer interaction with an adjacent rigid body sheet and internal vibration. After that, the behavior of the layered compound composed of a plurality of layers of rigid sheets is calculated.

  Further, the interlayer interaction includes a mutual potential, and the mutual potential is preferably determined by microscopic simulation at an electron level and an atomic level.

  In addition, it is preferable that the interaction potential has at least one of a two-dimensional periodicity and a terminal effect.

  The internal action is preferably determined at a microscopic level at the atomic level.

  Furthermore, it is preferable to simulate the behavior of the layered compound by applying pressure and shear to the layered compound sandwiched between two solid atomic layers to obtain a mechanical response including shear stress.

  Further, the present invention is a program for performing such a simulation.

  According to the present invention, it is possible to efficiently calculate the behavior of the layered compound such as the mechanical properties such as the frictional force and the thermal stability. This makes it possible to perform molecular design of a layered compound exhibiting desired molecular characteristics.

It is a figure explaining the friction mechanism of a layered compound. It is a figure explaining coarse-graining of a graphene sheet and 6 degrees of freedom. It is a figure which shows the multiscale analysis scheme of a layered compound. It is a figure which shows the molecular structure of a graphite. It is a figure which shows the Lennard-Jones potential between carbon atoms. It is a figure which shows the potential (Ca and Cb are the lattice points on between a and b) between a carbon atom and a graphene sheet. It is a figure which shows the potential (x, y space) between a carbon atom and a graphene sheet. It is a figure which shows the potential (x, y space) between a carbon atom and a graphene sheet. It is a figure which shows the model of simulation. It is a figure which shows the displacement of the x direction of the sheet | seat in a thermal equilibrium state. It is a figure which shows the displacement of the y direction of the sheet | seat in a thermal equilibrium state. It is a figure which shows the fluctuation | variation of the total energy U of the sheet | seat in a thermal equilibrium state. It is a figure which shows the "transfer piece" model on HOPG (highly oriented graphite). It is a figure which shows the displacement of the x direction of the sheet | seat under sliding motion in an asymmetrical system. It is a figure which shows the displacement of the y direction of the sheet | seat under sliding motion in an asymmetrical system. It is a figure which shows the time-dependent change of the total energy U of the sheet | seat under the sliding motion in an asymmetrical system. It is a figure which shows the displacement of the x and y direction of the sheet | seat under sliding motion in an asymmetrical system. It is a flowchart of calculation.

  Hereinafter, embodiments of the present invention will be described with reference to the drawings.

"Configuration of the simulation system"
The present embodiment relates to a simulation system for performing simulation using a model for a layered compound and a program thereof, and a program for performing simulation using the model is installed in a storage unit of a computer by CD, DVD, communication, etc. Is executed and various conditions and data are input to perform a simulation for obtaining physical properties such as a friction phenomenon and a friction coefficient. In addition, it is also suitable to use what is already sold for the simulation program itself, and to set the model in the simulation program by the program of this embodiment.

  The simulation is performed on an actual layered compound such as graphite or mica, and the structure of a functional molecule composed of the layered compound is designed.

"Behavior analysis of layered compounds"
For example, let us consider a state in which a steel ball having a radius of several millimeters is in contact with a graphite base on which infinitely long plates (graphene [1]) are stacked to make a relative sliding motion. At this time, it is unlikely that slip will occur between the graphene layers (if the interlayer friction force per unit area is finite, the force required for interlayer slip is infinite). In an actual substance, even highly oriented graphite has a curvature of the order of μm, so that it is geometrically difficult for shear to occur inside the underlying solid.

  Therefore, from a practical point of view, as shown in FIG. 1, transfer of graphite fragments occurs from the base graphite to the steel ball side, and sliding friction occurs between the transfer piece and the base graphite. Should be considered. The steel balls used for mechanical testing have a surface roughness of the order of sub-μm, so it is thought that there is a transfer piece at the tip of the steel ball. Some studies have been observed by Electron Microscope (TEM) [2]. Then, the mechanism of friction between the transfer piece and the base graphite surface is questionable.

  In this embodiment, whether there is a slip between the transfer piece and the base, or whether there is a slip between the transfer pieces, and the relationship between the shape and size of the transfer piece. Handle.

  Regarding the low-friction mechanism of graphite, there is an example in which the behavior between two layers is analyzed by the All-Atom Molecular Dynamics (MD) method [3], but the tip of Atomic Force Microscope (AFM) Even the piece of graphite transferred to the substrate actually consists of hundreds of thousands of atoms, so it is difficult to handle the dynamics of such a whole system by the all-atom molecular dynamics method. Therefore, in this embodiment, the transfer piece is coarse-grained to simplify the calculation.

  By the way, there are many modes in the internal vibration of the graphene sheet composing the transfer piece. If each graphene sheet is regarded as one rigid body, the fluctuation of the electron cloud due to this internal vibration is van der Waals force. It can be regarded as the fluctuation force exerted on the mass point of the adjacent rigid body as the fluctuation of the vibration. Further, this rigid body has six degrees of freedom of translation 3 and rotation 3. Of these, three important motions are translation and four yaw angles which are rotations about the axis in the layer direction. The sliding motion of the transfer piece can be described as the motion of a stack of graphene sheet rigid bodies having a Brownian motion in four degrees of freedom (see FIG. 2).

  In this embodiment, this is calculated by the Monte Carlo Brownian Dynamics (MCBD) method [4] to simulate the frictional behavior of the layered compound. The MCBD method was proposed by Kikuchi et al. In 1991 as a method for calculating colloidal solutions, but as described below, it is equivalent to Langevin dynamics and a phenomenon in which thermal fluctuations are on the same order as the interparticle potential. Is an effective coarse-grained simulation method.

"Method of simulation"
"MCBD method"
The molecular dynamics (MD) method and the Monte Carlo (MC) method are used as representative methods of computer simulation for handling multi-particle systems. The MD method tracks the trajectory of particles by simultaneously solving and solving equations of motion for each particle. On the other hand, the MC method calculates the configuration integral in statistical mechanics using random numbers. Each method has many variations depending on the degree of coarse-graining and the required physical quantity, but the difference is that the basis of the algorithm is Newtonian mechanics and the latter is statistical mechanics.

  The Brownian motion is macroscopically described by the Fokker-Planck diffusion equation, but microscopically, there are many theories that debate from the Langevin equation [5]. The Langevin equation takes a form that includes a term due to the fluctuation force due to the thermal motion of solvent molecules, which is added to the equation of motion stochastically. As a simulation method, the MD method using random numbers for the fluctuation term of the Langevin equation has been developed [6]. On the other hand, a representative method of the MC method is the method of Metropolis et al. [7]. This is based on the concept of a Markov chain that repeats a transition in which the transition probability at time n is determined only by the state at time n-1. In the normal Metropolis MC method, there is no clear correspondence between the state transition interval and the physical “time”.

  By the way, in the Brownian motion, it is assumed that the particle movement is Markovian, so in the Metropolis MC method, if the correspondence with time can be taken in the process of state transition, this method is also used in the Brownian motion. Should be a simulation method.

  The Metropolis MC method assumes an equilibrium system, and the Markov chain generated in the transition process is introduced to obtain the configuration integral in the thermal equilibrium state. Therefore, it cannot be directly applied to simulation of Brownian motion in non-equilibrium systems including transient processes. However, since the transition process in Brownian motion is also Markov, if a certain condition is added to the method of giving the transition probability, the Markov chain realized by the Metropolis MC method is a Brown particle corresponding to the physical "time". It should be the trajectory of [4].

In the MC step, “time” t is introduced as a label for identifying the state of the system, and the probability distribution function for the arrangement X at “time” t is P (X, t). At this time, since sampling assumes a Markov process, the change in “time” of P (X, t) is determined by the following master equation [5].

Here, W (X | X ′) is a transition probability per unit “time” regarding the transition from the arrangement X ′ to X. The master equation can be solved approximately by developing the parameters representing the size of the system. However, it is assumed that W (X | X ′) has the following canonical form.

Here, x ′ = X ′ / Ω and r = XX ′. At this time, the jump efficiency_α _{ν, λ} (x) is defined by the following equation.

When α ₁ , ₀ , this expansion gives the following equation as the first approximation of the nonlinear Fokker-Planck equation.

However, τ = ω ⁻² t. In the case of Metropolis MC, ω can be taken as the reciprocal of the maximum displacement width in the “time” interval Δt. Therefore, the first term on the right side of equation (2) is
Is calculated. For example, when the derivative U ′ (x) of x of the potential energy U (x) of the system is negative, for sufficiently large Ω
It is. Therefore, the jump efficiency is
Equation (4) results in the following diffusion equation:
Here, the diffusion constant D is defined by the following equation.

That is, by making the maximum displacement width Ω ⁻¹ in the Metropolis MC method small enough that the force acting on the particles can be considered constant for the MC “time” interval Δt, the Metropolis MC method becomes a method for simulating Brownian motion. It has been shown. At this time, MC “time” corresponds to physical time via equation (7).

  The diffusion coefficient can be determined using experimental values or by a more microscopic (atomic level) simulation method. Examples of using experimental values include low molecular ions in a polymer electrolyte solution. Yes [8].

  The advantage of this method over the method based on the MD method is that it is not necessary to calculate the forces acting on the particles, only the change in potential energy U (x) need be calculated. In addition to the thermal equilibrium value as compared with the general MC method, a change with time of the target physical quantity is also given. For this purpose, integration over a large number of trajectories is required.

  FIG. 3 shows an example of a multi-scale analysis scheme centering on the present MCBD method for handling the dynamics of layered compounds. As described above, the physicochemical reality of simulation in the MCBD method is to make the potential energy U (x) of the interlayer interaction and the diffusion coefficient D of the sheet realistic.

  As an example in the present invention, U (x) uses a literature value proposed by Girifalco et al. Well-known as van der Waals potential between carbon atoms of graphene (graphite, carbon nanotube, C60, etc.) [9]. As for the diffusion constant D, if there is an experimental value with high reliability, it may be used, or it can be determined by performing molecular dynamics of all atoms at a more microscopic level. If the force field parameters necessary for molecular dynamics of all atoms do not exist or if reliability is a problem, further force field parameters can be determined quantum-chemically by ab initio molecular orbital calculations.

  By using this multi-scale analysis scheme, it is possible to handle the dynamics of arbitrary layered compounds, including not only van der Waals systems such as graphite and h-BN, but also ionic systems such as molybdenum disulfide and mica. It becomes. Therefore, this method can be applied to the analysis and development of molybdenum-based layered compounds that are well known as engine oil additives and solid lubricants.

"Interaction between graphene sheets"
The structure of graphite at the molecular level is shown in FIG. In an AB stack structure in a thermal equilibrium state, the distance between carbon atoms la is 1.42 mm, and the distance between sheets lz is 3.35 mm. The interaction between graphene sheets was based on literature values [9] of van der Waals potential between carbon atoms by Girifalc et al. At this potential, the AB stack structure, which is a stable structure, is also reproduced, and the peeling stress between the graphene sheets is also reproduced. By introducing this potential as an interaction between rigid bodies, it is possible to capture the anisotropy of the interlayer interaction and the energy topography according to the direction of motion.

"Carbon interatomic potential"
The interatomic interaction Uatom-atom (ra) is expressed as the following Lennard-Jones type function.

Here, ra is an interatomic distance, and A and B are constants. This function can be rewritten using the index ε representing the depth of the potential and σ related to the equilibrium distance.

ε and the equilibrium distance ra0 are expressed as follows.

The following values were used as potentials between carbon atoms on graphene. A = 15.2 eV × Å ⁶ , B = 24.1 × 10 −3 eV × Å ¹² , r 0 = 3.83 Å, ε = 2.89 meV. FIG. 5 shows the interatomic potential equation (8) expressed as above.

"Carbon atom-graphene potential"
Next, the potential between carbon atom and graphene is determined by the following equation based on the equation (8).

Here, in the sheet orthogonal coordinate system, x is a slip direction, y is a direction perpendicular to the slip, and the sum of the right side is a carbon atom group having a honeycomb structure of one side la arranged on a plane z = 0, and arranged at z = 1z. Let it be an interaction with a carbon atom. In the simulation, in the coordinate system having the principal axis vectors a = ((√3 / 2) la, (1/2) la), b = (0, la) in consideration of the symmetry of the system, the directions a and b are used. The Uatom-sheet is calculated in advance on the lattice points divided into 100 (Fig. 6), and the sheet potential is returned to the orthogonal coordinate system when the sheet moves to calculate the inter-sheet potential. The Uatom-sheet at each lattice point was accumulated from the vicinity of the atom to the far side until the difference ΔUatom-sheet <10 ⁻⁶ eV. FIG. 7 shows Uatom-sheet (x, y) in the orthogonal coordinate system.

"Graphene-graphene potential"
The atoms arranged in two opposing honeycomb structures are geometrically distinguished into two types depending on the periodicity, so the potential Usheet-sheet between graphene sheets uses Uatom-sheet (x, y),
Can be calculated by Here, Δx and Δy are differences in coordinates in the x and y directions between the sheets, and Nx, N′x, Ny, and N′y are the numbers of atoms in the x and y directions of the opposite sheets. Equation (13) makes it possible to incorporate end effects in the interaction of sheets of finite size. The potential surface of the Usheet-sheet in the system of Nx = N′x = 100 and Ny = N′y = 100 is shown in FIG. 8, and cross sections on the Δx and Δy axes are shown. When Δx = 0 and Δy = 0, this corresponds to the AB structure and is the most stable, and it is shown that the sheet-to-sheet interaction is weakened due to the sheets being displaced.

When one side of the sheet is a highly oriented pyrolytic graphite (HOPG) with no end, the potential Usheet-plate is
Is represented by Here, x and y are the coordinates of the sheet, and Nx and Ny are the number of atoms in the x and y directions of the sheet. In this interaction, there is no change in the baseline of the potential surface shown in FIG. 8, and it becomes a periodic function.

"Simulation model"
An MCBD simulation model of the graphene transfer piece is shown in FIG. A graphene sheet of Nz layers (Nz layers) is sandwiched between the upper and lower sliders. Each sheet is assumed to be composed of Nx and Ny atoms with respect to the sliding direction x and the perpendicular direction y of the slider, respectively. For example, in the system of (Nx, Ny, Nz) = (100, 100, 20), the number of atoms in the transfer piece is 200,000, and it is very difficult to handle dynamics by molecular dynamics based on the all-atom model. On the other hand, in the MCBD method of the present embodiment, the calculation time is proportional to Nz and does not depend on Nx and Ny. This is because the atom-graphene interaction potential Uatom-sheet having no terminal effect is first determined in the above-described calculation of the interlayer interaction, and then the graphene-graphene interaction potential Usheet- including the terminal effect is determined. This is because the procedure for determining the sheet was taken. The transfer piece sizes in the system of (Nx, Ny, Nz) = (100, 100, 20) are 12.3 nm, 21.3 nm, and 6.4 nm in the (x, y, z) direction, respectively. It is the same scale as the size of the transfer piece attached to the tip of the AFM tip.

  The physical time Δt and the time between MC steps in the MCBD method are related by the diffusion coefficient D and the equation (4). The value of D is an experimental value or a more micro-scale simulation, for example, a finite graphene sheet Can be obtained by performing a molecular dynamics calculation of all atoms for the motion of.

The maximum displacement width was Ω ⁻¹ = 0: 01 mm, the interlayer distance was a very light load condition, a thermal equilibrium value of 3.35 mm, and the sliding speed was 10 ⁻⁶ mm / MCstep. The temperature T is a room temperature of 300 K. In most of the simulations shown below, as an initial arrangement of sheets, x = 0 and y = 0 in which all layers are stacked.

  As described above, when each graphene sheet is regarded as a rigid body, it has degrees of freedom of translation 3 and rotation 3 (right in FIG. 2). As an example, the result of focusing on only x and y motions is shown.

"result"
"Thermal fluctuation of graphite in equilibrium"
10, 11, and 12 plot the fluctuations of x, y and the total energy of the system in the equilibrium state of the system of graphene sheets each having 100 × 100 atoms and 20 layers. The 0th (upper slider), 1, 5, 10, 15, and 20th sheets are shown.

  Here, FIG. 10 shows the displacement of the sheet in the x direction in the thermal equilibrium state. Hereinafter, the numbers in the legend in the figure indicate that the 0th (upper slider), No. 1, 5, 10, 15, and 20 sheets are plotted. Nx = Ny = 100, Nz = 20, degrees of freedom: x, y, end effect: upper and lower, transfer piece shape: rectangular parallelepiped, sliding speed: 0. FIG. 11 shows the displacement in the y direction of the sheet in a thermal equilibrium state. Nx = Ny = 100, Nz = 20, degrees of freedom: x, y, end effect: upper and lower, transfer piece shape: rectangular parallelepiped, sliding speed: 0. FIG. 12 shows the fluctuation of the total energy U of the sheet in the thermal equilibrium state. Nx = Ny = 100, Nz = 20, degrees of freedom: x, y, end effect: upper and lower, transfer piece shape: rectangular parallelepiped, sliding speed: 0.

  The end effect is taken into consideration for both the upper slider and the lower slider, and as shown in FIG. 8, the potential surface has an energy disadvantageous deviation from the equilibrium position (origin) in each sheet. The fluctuation is small in both the x and y directions, staying in the vicinity of the origin, and each sheet does not exceed the convex portion of the potential in FIG.

  When observed in detail, the fluctuations of No. 1 and No. 20 in contact with the upper and lower sliders are suppressed as compared with other sliders.

"Slip dynamics"
Consider the state of the transfer piece in the experiment. Since the transfer piece is considered to be attached to the upper slider that moves integrally with the steel ball, it is natural to assume that the upper slider is a sheet having a finite area, that is, a terminal effect. Since the lower slider is highly oriented graphite that is expected to have high flatness, it is considered appropriate to place it as a periodic function represented by the formula (14). FIG. 13 shows the relationship between the asymmetric vertical slider.

  The result of the shear simulation is shown in FIG. FIG. 14 shows the displacement in the x direction of the sheet under sliding motion in an asymmetric system, Nx = Ny = 100, Nz = 20, degrees of freedom: x, y, end effect: upper, transfer piece shape: rectangular parallelepiped, Slip: This is a condition. As the upper slider moves, you can see the whole transfer piece following.

  Further, when attention is paid to the motion of the twentieth layer in contact with the lower slider (HOPG side), it can be seen that it periodically vibrates with the motion following the upper slider.

  FIG. 15 shows the displacement in the y direction in the simulation, and FIG. 16 shows the change in energy. 15 and 16 also have the conditions that Nx = Ny = 100, Nz = 20, degrees of freedom: x, y, end effect: top, transfer piece shape: rectangular parallelepiped, slip: existence. Also in the y direction, it can be read that the twentieth layer vibrates periodically. In addition, the energy vibrates at a period twice that of the 20th layer, suggesting that the frictional expression is closely related to the interaction between the lowermost layer of the transfer piece and the HOPG surface.

  FIG. 17 shows the displacement of the 20th layer in the x and y planes. This FIG. 17 also has the conditions that Nx = Ny = 100, Nz = 20, degrees of freedom: x, y, end effect: top, transfer piece shape: rectangular parallelepiped, slip: existence. Compared with the potential surface of FIG. 8, it can be seen that the twentieth layer is moving along the equipotential surface.

  From the above, it can be seen that in the asymmetric system imitating the steel ball-transfer piece-HOPG surface, the transfer piece moves together with the steel ball and slips by movement along the equipotential surface of HOPG.

  From the viewpoint of the transfer piece, the transfer piece has no energy loss regardless of the x-axis direction with respect to the lower slider, but the upper slider has a finite effect because it has an end effect. It is more energetically advantageous to follow the slider. It can be seen that it is important that the slider has a finite size.

  As described above, the state of interlayer slip under various conditions can be captured by the MCBD calculation. The base was graphene represented by a periodic potential, and the upper part of the transfer piece was a slider that moved at a constant speed. In a transfer piece having a size of about 100 atoms (12 nm) on one side, the transfer piece moves together with the slider, and friction occurs between the lowermost layer and the base.

  In the sliding friction of the layered compound, it is understood that the energy topography between sheets at such a realistic sliding speed and temperature is very important.

  In the present embodiment, the degrees of freedom of x and y have been described, but it is naturally possible to incorporate the simulator into the z direction and the yaw angle. At the same time, the diffusion coefficient D can be determined by determining the total atomic molecular dynamics.

  By using the simulator of this embodiment, if an appropriate transfer piece size effective for low friction is found, direct observation, improvement of surface properties that actively promote transfer, etc. are considered. It is done. Further, it is possible to search for a new low-friction material by using the multi-scale analysis scheme shown in FIG.

"Processing Procedures and Effects"
FIG. 18 shows the overall processing flow of calculation in this embodiment. First, a candidate for a layered compound (layered molecule) is selected (S1). As described above, there are various layered compounds such as graphite, mica, and molybdenum disulfide, and the substance to be calculated is specified.

  Based on the literature value, experimental value, or quantum calculation result of the target substance, the interlayer interaction is determined and input (S2). That is, if a layered compound is specified, a basic configuration such as an atomic arrangement of the layered compound can be obtained from literature. Therefore, the interaction such as the mutual potential between layers is determined from literature values, experimental values, and the like. In addition, as for the interaction such as the mutual potential, it is often the case that the literature value cannot be used as it is, and the layer-to-layer (sheet-to-sheet) interaction as described above is determined based on the basic layer compound configuration. Is preferred. For example, the value of equation (13) is determined.

  In particular, according to this embodiment, before calculating the sheet-to-sheet interaction (mutual potential), the atom-to-sheet interaction is calculated without considering the end effect, and the end effect is included based on the result. Calculate sheet-to-sheet interaction. Accordingly, the sheet-to-sheet calculation with the end effect effective is facilitated. And by such a calculation, the interaction between layers can be prescribed | regulated as a characteristic of one sheet | seat.

  Next, parameters (diffusion coefficient, swinging force, etc.) indicating the action caused by internal vibration of each layer are determined and input based on literature values, experimental values, or calculation results at the atomic level (S3). The diffusion coefficient D and the like can also be determined by calculation based on basic data as described above. In addition, if it can be determined directly by literatures, it is used. For example, the value of equation (7) is determined thereby.

  As described above, since the diffusion coefficient and the like can also be described as the layer characteristics, each layer can be coarse-grained as a rigid sheet in which the mutual potential between the layers and the diffusion coefficient are specified.

  The initial conditions such as the initial arrangement of each layer and the external field are input (S4). As described in the section of “Simulation Model”, the maximum displacement width, interlayer distance, temperature, sheet arrangement, and the like are determined.

  Then, structural calculation or dynamic calculation is executed (S5). For example, the calculations shown in the sections of “thermal fluctuation of graphite in equilibrium” and “slip dynamics” are executed.

  Such a simulation is performed by applying a force on the coarse-grained rigid sheet or moving the upper one. Therefore, the calculation can be performed with a relatively simple calculation, and the mutual potential and diffusion coefficient can be calculated in consideration of the influence of the atomic level.

  When such calculation is completed, an execution result is output (S6). In other words, the behavior of the target layered compound under various conditions is simulated, and as a result, the target layered compound is evaluated and the structure of the layered compound appropriate for the specified condition is determined. it can.

  Conventionally, when the number of layered compound pieces is Nxy per sheet, and the whole is composed of Nz layers, it is necessary to handle Nxy × Nz particles in the molecular simulation based on the all-atom model, which requires enormous calculation time. . As an example, the graphite transfer piece attached to the tip on the atomic force microscope has Nxy of several tens of thousands and Nz of several tens, and it requires several months of calculation by a large computer to calculate one physical quantity. . According to this embodiment, the influence of the vibration of the atoms constituting each layer of the layered compound can be transferred to the interaction and diffusion coefficient in each layer unit, limiting the internal degrees of freedom, thereby handling particles The number is reduced to Nz, and Nxy times calculation is speeded up.

  In addition, since the interlayer interaction and the like are based on microscopic calculations, the characteristics of each compound can be reproduced. Thereby, the molecular design of a novel layered compound can be performed. Furthermore, in this method, since the system anisotropy is taken such that the diffusion coefficient in the thickness direction of the layer is small and the diffusion coefficient in the direction perpendicular thereto is large, it is possible to realize a calculation reflecting the characteristics of the layered compound.

References
In the above description, the cited documents are shown below.
[1] AK Geim and KS Novoselov. Aaa. Nature Materials, Vol. 6, p. 183, 2007.
[2] M. Dienwiebel et al. Aaa. Surface Science, Vol. 576, p. 197, 2005.
[3] M. Tsukada N. Sasaki, K. Kobayashi. Atomic-scale friction image of graphite in atomic-force microscopy. Phys. Rev. B., Vol. 54, pp. 2138 {2149, 1996.
[4] K. Kikuchi, M. Yoshida, T. Maekawa, and H. Watanabe. Metropolis monte carlo method as a numerical technique to solve the fokker-planck equation. Chem. Phys. Lett., Vol. 185, pp. 335 {338, 1991.
[5] NG van Kampen. Stochastic processes in physics and chemistry, revised and enlarged edition, chapter 11. North-Holland, Amsterdam, 1992.
[6] DL Ermak and JA McCammon. Brownian dynamics with hydrodynamic interac-tions. J. Chem. Phys., Vol. 69, pp. 1352 {1360, 1978.
[7] N. Metropolis, AW Rosenbluth, MN Rosenbluth, AH Teller, and E. Teller. Equation of state calculations by fast computing machines. J. Chem. Phys., Vol. 21, pp. 1087 {1092, 1953.
[8] M. Yoshida and K. Kikuchi. Metropolis monte carlo brownian dynamics simulation of the ion atmosphere polarization around a rodlike polyion. J. Phys. Chem., Vol. 98, pp. 10303 {10306, 1994.
[9] M. Hodak LA Girifalco and RS Lee. Aaa. Phys. Rev. B., Vol. 62, p. 13104, 2000.

Claims (6)

Coarse each layer of the layered compound as a rigid sheet that is a mechanical unit,
For the rigid sheet of each layer coarse-grained, determine the interlayer interaction with the rigid sheet of the adjacent layer, the internal action due to internal vibration,
After that, a simulation system characterized in that the calculation of the behavior of the layered compound composed of a plurality of layers of rigid sheets is performed. The simulation system according to claim 1,
The interlayer interaction includes a mutual potential, and the mutual potential is determined by microscopic simulation at an electron level and an atom level. The simulation system according to claim 2,
A simulation system, wherein the interaction potential has at least one of a two-dimensional periodicity and a terminal effect. The simulation system according to claim 1,
The simulation system characterized in that the internal action is determined at a microscopic level at an atomic level. In the simulation system according to any one of claims 1 to 4,
Furthermore, a simulation system characterized by simulating the behavior of a layered compound by applying pressure and shear to a layered compound sandwiched between two solid atomic layers to obtain a mechanical response including shear stress. On the computer,
For each layered compound, each layer is coarse-grained as a rigid sheet that is a mechanical unit,
For the coarse sheet of each layer, the interlayer interaction with the rigid sheet of the adjacent layer, the internal action due to internal vibration is determined,
After that, the calculation of the behavior of the layered compound consisting of multiple layers of rigid sheets is performed.
A simulation program characterized by that.