JP2006285866A - Simulation method and program - Google Patents

Abstract

<P>PROBLEM TO BE SOLVED: To provide a new simulation method using molecular kinetics. <P>SOLUTION: In a considered particle system, the number of particles is N, mass of each the particle is m, and interaction potential is expressed by a product of an interaction coefficient ε and a dimensionless function f representing inter-particle distance dependence. The number N' of particles, mass m', and an interaction coefficient ε' are respectively found by conversion expressions N'=N/α<SP>d</SP>, m'=mα<SP>2</SP>/α<SP>γ</SP>, and ε'=εα<SP>γ</SP>by use of a spatial dimension number d, a first transfer factor α larger than 1, and a second transfer factor γ not less than 0 and not more than d. Molecular kinetics calculation is executed to a particle system, wherein the number of the particles is N'; the mass of each the particle is m'; and the interaction potential is expressed by ε'f. <P>COPYRIGHT: (C)2007,JPO&INPIT

Description

  The present invention relates to a simulation method and a program, and more particularly to a simulation method using molecular dynamics and a program for causing a computer to execute the simulation method.

  Computer simulations using molecular dynamics have been performed. In molecular dynamics, the equations of motion of the particles that make up the system to be simulated are numerically solved. As the number of particles in the system to be simulated increases, the amount of calculation required increases. With the computing power of current computers, for example, it is only possible to perform a simulation of a system having about 100,000 particles.

In order to reduce the amount of calculation required for the simulation, for example, Non-Patent Document 1 attempts to apply a renormalization technique to molecular dynamics.
Inamura, Takezawa, Shamoto, “Variable scale molecular dynamics simulation based on renormalization technique”, Transactions of the Japan Society of Mechanical Engineers (A), 1997, 63, 608, p. 202-207

  As described in “3. Numerical calculation method” of Non-Patent Document 1, the simulation method of Non-Patent Document 1 cannot positively evaluate the temperature.

  An object of the present invention is to provide a novel simulation method using molecular dynamics and a program for executing the method on a computer.

  Another object of the present invention is to provide a simulation method capable of positively evaluating the temperature based on the renormalization group technique, and a program for executing the method on a computer.

According to one aspect of the present invention, (a) a non-dimensionalization function f that includes N particles, each particle has a mass m, and the interaction potential energy between the particles represents the dependence on the interparticle distance. And the particle system S that can be expressed by the product εf of the interaction coefficient ε, the dimension number d of the space in which the particle system S is arranged, the first renormalization factor α greater than 1, and 0 or more Using the second renormalization factor γ equal to or less than d, the number N ′ of renormalized particles is obtained by the conversion equation N ′ = N / α ^d , and the mass m ′ of the renormalized particles is determined by the conversion equation m. ′ = Mα ² / α ^{γ, the step} of determining the renormalized interaction coefficient ε ′ by the conversion equation ε ′ = εα ^γ , and (b) the number of particles N ′ transferred Each particle has a mass m ′ of renormalized particles, and the interaction potential energy between the particles is And Motoka function f, the particle system S'represented by the product ε'f the interaction coefficient ε'was renormalized, a simulation method and a step of performing a molecular dynamics calculation is provided.

  The number N ′ of particles brought in is smaller than the number N of particles. For this reason, the molecular dynamics calculation for the particle system S ′ can be executed with a smaller amount of calculation than when the molecular dynamics calculation is performed for the particle system S.

  First, the molecular dynamics (Molecular Dynamics; MD) will be briefly described. It consists of N particles (for example, atoms), and the Hamiltonian H is

と表される粒子系について考える。ここで、各粒子の質量がｍであり、粒子ｊの運動量ベクトルがｐ_ｊであり、粒子ｊの位置ベクトル（位置座標）がｑ_ｊであり、粒子間の相互作用ポテンシャルエネルギがφである。

Consider the particle system expressed as Here, the mass of each particle is m, the momentum vector of the particle j is p _j , the position vector (position coordinate) of the particle _j is q _j , and the interaction potential energy between the particles is φ.

  Equation of motion for each particle by substituting Hamiltonian H into Hamilton's canonical equation

及び

as well as

が得られる。ここで式（３）において、ｖ_ｊは粒子ｊの速度ベクトルである。分子動力学では、粒子系を構成する各粒子に対し、式（２）及び（３）で表される運動方程式を数値的に積分して解くことにより、各時刻における各粒子の運動量ベクトル（または速度ベクトル）及び位置ベクトルが求められる。なお、数値積分には、Ｖｅｒｌｅｔ法が用いられることが多い。Ｖｅｒｌｅｔ法については、例えば、Ｊ．Ｍ．Ｔｈｉｊｓｓｅｎ，“Ｃｏｍｐｕｔａｔｉｏｎａｌ Ｐｈｙｓｉｃｓ”，Ｃａｍｂｒｉｄｇｅ Ｕｎｉｖｅｒｓｉｔｙ Ｐｒｅｓｓ（１９９９）のｐ．１７５に説明されている。分子動力学計算で得られた各粒子の位置、速度に基づき、系の種々の物理量を算出することができる。

Is obtained. Here, in Expression (3), v _j is a velocity vector of the particle j. In molecular dynamics, the momentum vector (or each particle at each time) (or by solving numerically integrating the equations of motion represented by equations (2) and (3) for each particle constituting the particle system. Velocity vector) and position vector. Note that the Verlet method is often used for numerical integration. Regarding the Verlet method, for example, J. et al. M.M. Thijssen, “Computational Physics”, Cambridge University Press (1999), p. 175. Various physical quantities of the system can be calculated based on the position and speed of each particle obtained by molecular dynamics calculation.

  Next, molecular dynamics applying the renormalization group technique (hereinafter referred to as renormalization group molecular dynamics) will be conceptually described.

  In the renormalization group molecular dynamics, a desired system S for which a physical quantity is desired is associated with a system having a smaller number of particles than the system S (hereinafter referred to as a renormalized system S ′). Next, molecular dynamics calculation is performed on the renormalized system S ′. Then, the calculation result for the renormalized system S ′ is associated with the desired system S. As a result, it is possible to calculate the physical quantity for the desired system S with a smaller amount of calculation than when the molecular dynamics calculation for the desired system S is directly executed. A conversion rule for associating a physical quantity (for example, the number of particles, the mass of each particle, etc.) in the desired system S with a physical quantity in the renormalized system S ′ is called a scale conversion law.

  Next, the scale conversion rule found by the present inventors will be described. It is assumed that a desired particle system S for which a physical quantity is to be calculated is composed of N particles. As shown in equation (1), the total Hamiltonian H of the particle system S is

と表されるとする。ある粒子ｊのハミルトニアンＨ_ｊは、

It is assumed that The Hamiltonian H _j of a particle j is

と表される。

It is expressed.

  For the renormalization discussion, the interaction potential energy φ represents the dependence on the interparticle distance, the dimensionless function f, the interaction strength, and the coefficient ε (which has the energy dimension) It is expressed as follows by the product of the interaction coefficient ε).

例えば、不活性原子に対しては、

For example, for inert atoms:

という相互作用ポテンシャルが用いられる。ここでｒは原子間距離である。相互作用係数ε及び定数σが、原子の種類に対応する値を取る。

The interaction potential is used. Here, r is an interatomic distance. The interaction coefficient ε and the constant σ take values corresponding to the types of atoms.

  Next, a method for deriving the renormalized Hamiltonian H ′ from the Hamiltonian H will be described. As will be shown in detail below, the renormalized Hamiltonian H ′ performs part of the integration in the partition function Z (β) for the particle system S, coarse-grains the Hamiltonian, and then rescales the integration variable. It is obtained with.

  The partition function Z (β) for the canonical ensemble at a fixed number of particles is

と表される。ここで、係数βは、系の温度Ｔとボルツマン定数ｋ_Ｂとによりβ＝１／（ｋ_ＢＴ）と定義される。ｄΓ_Ｎは位相空間内の体積要素であり、より詳細には、

It is expressed. Here, the coefficient β is defined as β = 1 / (k _B T) by the system temperature T and the Boltzmann constant k _B. dΓ _N is a volume element in the phase space, and more specifically,

と表される。ここで、Ｗ_Ｎは、Ｎ！ｈ^３Ｎ（ｈはプランク定数）である。また、Ｄ_ｐ ^Ｎ及びＤ_ｑ ^Ｎをそれぞれ、

It is expressed. Here, W _N is N! h ^3N (h is Planck's constant). Also, D _p ^N and D _q ^N are respectively

It is defined as

It is difficult to discuss the renormalization of potential other than the harmonic oscillator. Therefore, the interaction potential with saddle points is divided into three regions. Since the interparticle distances r → ∞ and r → 0 are fixed points, the potential does not change during the renormalization transformation. Therefore, it is assumed that the renormalization transformation in the vicinity of the saddle point holds for all the interparticle distances r. Develop the tail around the saddle point and leave it to the second fluctuation. If the saddle point position is r ₀ and the relative displacement is δu, the first-order term of δu disappears,

を得る。従って、粒子ｊについて、鞍点近傍を記述するハミルトニアンは、

Get. Therefore, for particle j, the Hamiltonian describing the vicinity of the saddle point is

となる。ここで、ｕ_ｉ及びｕ_ｊは、それぞれ、粒子ｉ及びｊの鞍点からの変位を表す。φ″は、粒子間距離ｒについての２階微分である。

It becomes. Here, u _i and u _j represent the displacements of the particles i and j from the saddle point, respectively. φ ″ is the second derivative with respect to the interparticle distance r.

  Next, coarse graining of the interparticle interaction portion of the partition function Z (β) will be described. First, after considering coarse graining in a particle system arranged on a one-dimensional chain, the system is extended to a particle system arranged on a simple cubic lattice.

As shown in FIG. 1A, particles i and j are arranged closest to each other at lattice points on the one-dimensional chain, and particles j and k are arranged closest to each other. ing. Particle j is arranged between particles i and k. The interaction between the adjacent particles (two arrest neighbor coupling) is shown by a solid line. Consider erasing particles j and coarsening. If the interaction involving the particle j is written out and integration is performed on the displacement u _j of the particle j in the middle between the particles i and k, the particle i and k after the influence of the particle j is renormalized. For the interaction potential φ ′ (q _i −q _k ),

という式を得る。

The following formula is obtained.

  Since the lattice constant becomes α (= 2) times by this coarse graining, it can be expressed in the same form as the original potential function by a variable u ′ scale-converted by the equation u ′ = u / α. In other words, the following similarity law is obtained.

ただし、ここで、式（１２）の２行目へ移るに際して、鞍点近傍以外においても同じ繰り込み変換が成立すると仮定した。αを、第１の繰り込み因子と呼ぶこととする。

However, here, when moving to the second line of the equation (12), it is assumed that the same renormalization conversion is established except in the vicinity of the saddle point. α is referred to as a first renormalization factor.

  Coarse graining in the d-dimensional simple cubic lattice can be realized by a method of potential moving. As for Potential Moving, for example, Leo P.I. Kadanoff, “STATISTICAL PHYSICS, Statics, Dynamics and Renormalization”, Chapter 14 of World Scientific (1999).

  With reference to FIG. 1 (B), the concept of Potential Moving will be described. In Potential Moving, a multidimensional lattice is reduced to a one-dimensional chain. FIG. 1B shows the case of a two-dimensional lattice.

  As shown in the upper diagram of FIG. 1B, particles are arranged on lattice points of a two-dimensional square lattice. The interaction between the most adjacent particles (inter-adjacent interaction) is indicated by a solid line. The direction parallel to the side with the lattice is defined as the X direction, and the direction orthogonal to the X direction is defined as the Y direction.

  As shown in the middle diagram of FIG. 1B, it is considered that every other particle arranged in the X direction is integrated (lattice points for summation). Consider integrating the particles above). A particle to be integrated (erased) is referred to as an integrated particle. The interaction between the nearest neighbors of the integrand particles is indicated by a wavy line. If there is no interaction between the nearest neighbors of the integrand particles, the integrand particles can be regarded as a one-dimensional chain extending in the X direction.

  In Potential Moving, as shown in the lower diagram of FIG. 1B, the interaction between the nearest neighbors of the integrable particles is distributed to the particles arranged on both sides of the integrand particles in the X direction. The interparticle interaction (double coupling) to which the distributed interaction is added is indicated by a double line. The interaction indicated by the double line is twice as strong as the original nearest neighbor interaction. In this way, the integrand particle can be converted from a two-dimensional lattice to a one-dimensional chain. Therefore, the above-described coarse graining method for the one-dimensional chain can be performed on the integrand particles. For example, in the case of a three-dimensional lattice, the same procedure as described above is repeated for the Y and Z directions. In this way, coarse graining in a multidimensional lattice is performed.

  By using the method of Potential Moving, coarse-graining of the interaction part of the partition function Z (β) in the d-dimensional simple cubic lattice,

という結果が導かれる。ここで、繰り込まれた粒子数Ｎ´及び第１の繰り込み因子αが、

The result is derived. Here, the number N ′ of renormalized particles and the first renormalization factor α are

という式で表される。ｎは正整数であり、粗視化の回数を表す。

It is expressed by the formula. n is a positive integer and represents the number of coarse-graining.

Next, coarse graining of the kinetic energy portion of the partition function Z (β) will be described. First, consider a one-dimensional chain. The particle pair j-1 and particle j, which are the particle pairs closest to each other, are coarse-grained into one particle. The concept of the Kadanoff block spin method is applied to the momentum, and the momentum g _j of the coarse-grained particle composed of the particle j and the particle j−1 is defined as follows.

また、重み関数Ｔ｛Ｐ，ｐ｝を以下の様に定義する。

Further, the weight function T {P, p} is defined as follows.

ここで、δはＤｉｒａｃのδ-関数である。明らかに、

Here, δ is a Dirac δ-function. clearly,

を満たす。この条件から、分配関数の運動エネルギ部分Ｚ_ｐへ重み関数Ｔ｛Ｐ，ｐ｝を挿入でき、

Meet. From this condition, the weight function T {P, p} can be inserted into the kinetic energy part Z _p of the partition function,

という式を得る。δ-関数を積分表示すれば、分配関数の運動エネルギ部分Ｚ_ｐは、以下の様に表せる。

The following formula is obtained. If the δ-function is displayed in an integral manner, the kinetic energy portion Z _p of the partition function can be expressed as follows.

この積分は、以下の公式を使えば容易に実行できる。

This integration can be easily performed using the following formula:

ここで、＜Ａ，Ｂ＞は、ベクトルＡとＢとの内積を表す。この公式を用いると、分配関数の運動エネルギ部分Ｚ_ｐは、

Here, <A, B> represents an inner product of the vectors A and B. Using this formula, the kinetic energy part Z _p of the partition function is

と表すことができる。

It can be expressed as.

  This coarse graining increases the kinetic energy by α (= 2) times. Similar results can be obtained when the total kinetic energy of two particles is divided into the kinetic energy of the center of gravity and the kinetic energy of relative motion, and integration is performed on the relative momentum.

Next, coarse graining in a d-dimensional simple cubic lattice will be considered. First, when coarse-graining is performed in the direction of the crystal axis a, the kinetic energy E _kin possessed by the particles becomes α times. Subsequently, when coarse-grained in the b-axis direction, the kinetic energy E _kin α of the particles becomes α times (E _kin α), that is, E _kin α ² . Therefore, in the d-dimensional lattice, E _kin α ^d ,

という関係式が得られる。ｎ回の粗視化では、α＝２^ｎとなることは自明である。

Is obtained. It is self-evident that α = 2 ^{n in} n times of coarse graining.

  From the above discussion, the partition function Z (β) for all coarse-grained Hamiltonians, except for the coefficients that do not affect the result,

と表すことができる。

It can be expressed as.

  Next, in order to make the Hamiltonian isomorphic even in the coarse-grained system, a “retracted” variable is introduced. The renormalized position vector q ′, the renormalized momentum vector p ′, the renormalized mass m ′, the renormalized interaction coefficient ε ′, and the renormalized coefficient β ′ as follows: Define The momentum vector is scaled so that the volume element in the phase space is invariant.

ここで、γは、０以上ｄ以下である。γを、第２の繰り込み因子と呼ぶこととする。

Here, γ is 0 or more and d or less. Let γ be called the second renormalization factor.

  By defining the renormalized variable in this way, the renormalized Hamiltonian H ′ having the same shape as the Hamiltonian H is obtained. The Hamiltonian H 'brought over is

と表される。

It is expressed.

  From Equation (18), the distribution function is as follows.

ここで、位相空間の体積素は、

Here, the volume element of the phase space is

である。

It is.

  The equation of motion in the system described by the renormalized Hamiltonian H ′ (the renormalized particle system S ′) is a canonical equation.

及び

as well as

から与えられ、

Given by

及び

as well as

となる。ここで、繰り込まれた時間ｔ´が、繰り込まれた変数ｑ´、ｐ´、ｍ´及びε´に整合するように、

It becomes. Here, the renormalized time t ′ matches the recursive variables q ′, p ′, m ′ and ε ′.

という式で定められる。運動方程式（２４）及び（２５）に基づいて、繰り込まれた粒子系Ｓ´に対する分子動力学計算を実行することができる。

It is determined by the formula. Based on the equations of motion (24) and (25), a molecular dynamics calculation can be performed on the renormalized particle system S ′.

  The above considerations are summarized. The particle system S and the brought-in particle system S ′ are associated by a scale conversion rule expressed by the following conversion formula.

  The number N of particles in the particle system S and the number N ′ of particles carried in the particle system S ′ are converted into a conversion formula.

により対応付けられる。粒子系Ｓのある粒子の位置ベクトルｑと、繰り込まれた粒子系Ｓ´のある粒子の位置ベクトルｑ´とは、変換式

Are associated with each other. The position vector q of the particle with the particle system S and the position vector q ′ of the particle with the renormalized particle system S ′ are converted to

により対応付けられる。粒子系Ｓのある粒子の運動量ベクトルｐと、繰り込まれた粒子系Ｓ´のある粒子の運動量ベクトルｐ´とは、変換式

Are associated with each other. The momentum vector p of a particle with a particle system S and the momentum vector p ′ of a particle with a renormalized particle system S ′ are converted into a conversion formula:

により対応付けられる。粒子系Ｓの各粒子の質量ｍと、繰り込まれた粒子系Ｓ´の各粒子の質量ｍ´とは、変換式

Are associated with each other. The mass m of each particle of the particle system S and the mass m ′ of each particle of the brought-in particle system S ′ are converted into equations

により対応付けられる。粒子系Ｓにおける相互作用係数εと、繰り込まれた粒子系Ｓ´における相互作用係数ε´とは、変換式

Are associated with each other. The interaction coefficient ε in the particle system S and the interaction coefficient ε ′ in the renormalized particle system S ′

により対応付けられる。

Are associated with each other.

  The time scale t when solving the equation of motion for the particle system S (a certain time interval t in the molecular dynamics calculation assumed to be performed for the particle system S) and the motion for the renormalized particle system S ′. The time scale t ′ (a time interval t ′ in the molecular dynamics calculation performed for the renormalized particle system S ′) in solving the equation is a conversion formula

により対応付けられる。

Are associated with each other.

  Note that the velocity vector v of the particle having the particle system S and the velocity vector v ′ of the particle having the regenerated particle system S ′ are converted because the velocity vector v is expressed as v = p / m. formula

により対応付けられる。

Are associated with each other.

Here, the number of dimensions of the space in which the particle system S is arranged is d. The first renormalization factor α is 2 ⁿ (n is a positive integer), and the second renormalization factor γ is 0 or more and d or less.

  When the value of the second renormalization factor γ is large, for example equal to d, a very large renormalized interaction coefficient ε ′ and a small renormalized mass m ′ appear in the numerical calculation. The second renormalization factor γ is preferably set to 0 for numerical calculation.

  When a one-dimensional chain spring-mass system having a mass m and a spring constant s is considered, a dispersion relation is obtained from the renormalized Hamiltonian as follows.

ここで、角周波数がωであり、波数がｋである。平衡状態における原子間距離がａである。なお、第２の繰り込み因子γは０とした。式（２６）は、長波長の極限において、

Here, the angular frequency is ω and the wave number is k. The interatomic distance in the equilibrium state is a. The second renormalization factor γ was 0. Equation (26) is in the limit of long wavelengths,

となり、第１の繰り込み因子αには無関係になる。

And is irrelevant to the first renormalization factor α.

  Next, an example of computer simulation based on renormalization group molecular dynamics will be described.

  First, a simulation in which heat is caused to flow into aluminum to raise the temperature and calculate the melting point and latent heat will be described. In the first simulation, for comparison, a conventional molecular dynamics calculation (without applying renormalization) was performed on a system consisting of 16000 aluminum atoms. Aluminum atoms are arranged in a three-dimensional space, and the dimension number d is three.

In the second and third simulations, molecular dynamics calculation was performed on a system incorporating the system corresponding to the first simulation. The first renormalization factor α is 2 in the second simulation and 2 ² (= 4) in the third simulation. From (conversion equation N), the number of system particles is 2000 in the second simulation, and the number of system particles is 250 in the third simulation. The second renormalization factor γ was set to zero. The conventional method (first simulation) corresponds to the case where the first renormalization factor α is set to 1.

  As the interatomic potential in the first simulation (that is, the interparticle potential in the original system), the Morse potential was used. This potential is

と表される。ここで、相互作用係数εは１．９２×１０^−２１[Ｊ]であり、定数Ａは２．３５×１０^１０[ｍ^−１]であり、定数ｒ_０は２．８６×１０^−１０[ｍ]である。粒子間距離がｒである。

It is expressed. Here, the interaction coefficient ε is 1.92 × 10 ⁻²¹ [J], the constant A is 2.35 × 10 ¹⁰ [m ⁻¹ ], and the constant r ₀ is 2.86 × 10 ⁻¹⁰ [J]. m]. The interparticle distance is r.

As the interparticle potential in the second and third simulations (that is, the interparticle potential in the renormalized system), the same one as in the equation (28) is used. Renormalized interaction coefficients in the second and third simulations are respectively obtained from (conversion equation ε). Since the second renormalization factor γ is 0, the interaction coefficients in the second and third simulations are both equal to the interaction coefficient ε in the first simulation. Constant r ₀ and constant A are equal to those of the first simulation.

The mass of aluminum atoms in the first simulation (that is, the mass of particles in the original system) is 4.48 × 10 ⁻²⁶ [kg]. The mass in which the particles are brought into the second and third simulations is obtained from (conversion equation m), respectively. The mass in the second simulation is 17.9 × 10 ⁻²⁶ [kg] obtained by multiplying the mass in the first simulation by 2 ² (= 4), and the mass in the third simulation is the same as in the first simulation. It becomes 71.7 * 10 < ^-26 [kg] which multiplied 4 < ² > (= 16) of the mass of.

Next, the initial arrangement of the particles will be described. Initially, the particles are placed on lattice points of a face-centered cubic lattice. In the first simulation, 40 layers in which 20 particles in the width direction and 20 particles in the depth direction are arranged overlap each other. Of the original system, 1 / alpha in the width direction, the depth direction 1 / alpha, and the portion of 1 / alpha height direction (1 / α ³ parts of the original system), corresponding to the renormalized system To do. In the second simulation, 20 layers in which 10 particles are arranged in the width direction and 10 particles in the depth direction are stacked, and in the third simulation, 5 particles are arranged in the width direction and 5 particles in the depth direction. The layer in which the particles are arranged is arranged so as to overlap 10 layers.

  Next, the time increment of numerical integration will be described. The time step width in the original system is Δt, and the particle speed is v. The time increment Δt is, for example,

を満たすように定められる（「＜＜」は、式の左辺が右辺より充分小さいという意）。ここでｒ_０は粒子間ポテンシャルにおける定数であり、粒子間距離の目安を与える。繰り込まれた系における時間刻み幅をΔｔ´とし、粒子の速さをｖ´とする。（変換式ｖ）よりｖ´＝ｖ／αであるので、繰り込まれた系における時間刻み幅Δｔ´は、例えば、αΔｔに設定される。

(“<<” means that the left side of the expression is sufficiently smaller than the right side). Here, r ₀ is a constant in the interparticle potential and gives a measure of the interparticle distance. Let the time increment in the renormalized system be Δt ′ and the particle speed be v ′. Since v ′ = v / α from (conversion equation v), the time step width Δt ′ in the renormalized system is set to αΔt, for example.

  The time interval of numerical integration was set to 5.0 [fs] in the first simulation, 10 [fs] in the second simulation, and 20 [fs] in the third simulation.

  Next, the initial speed given to the particles will be described. First, the relationship between the temperature of the original system and the temperature of the retreated system will be described. The temperature T of the original system is

という式で定義される。ここで、ｋ_Ｂはボルツマン定数であり、ｖ_ｉは粒子ｉの速度ベクトルであり、ｖ_ｇは粒子系の重心の速度ベクトルである。温度は、つまり、速度の“揺らぎ”である。繰り込まれた系の温度Ｔ´は、

It is defined by the expression. Here, k _B is the Boltzmann constant, v _i is the velocity vector of the particle i, and _vg is the velocity vector of the center of gravity of the particle system. Temperature is the “fluctuation” of speed. The temperature T ′ of the system brought in is

という式から定まる。

It is determined from the formula.

In the original system, the kinetic energy of a particle with mass m and velocity v is ½ mv ² . In the renormalized system, the kinetic energy of a particle with mass m ′ and speed v ′ is ½ m′v ′ ² . From (conversion equation m) and (conversion equation v),

という関係が得られる。よって、元の系の温度Ｔと、繰り込まれた系の温度Ｔ´とは等しくなる。第１〜第３のシミュレーションにおいて、粒子系の温度が３００[Ｋ]となるように、初期速度が与えられる。

The relationship is obtained. Therefore, the temperature T of the original system is equal to the temperature T ′ of the system that has been brought in. In the first to third simulations, the initial velocity is given so that the temperature of the particle system is 300 [K].

  Next, a method for causing heat to flow into the particle system will be described. In the first to third simulations, kinetic energy was given to the particles arranged on the bottom surface of the particle system to simulate the inflow of heat into the system. The surface into which heat flows will be referred to as a heating surface.

First, the heat inflow method (that is, the heat inflow method in the original system) in the first simulation will be described. It is assumed that kinetic energy ΔE is given for each particle on the heating surface during the time interval Δt. The change in the momentum vector of the particle is obtained from the following equation (32). Here, the momentum vector at a certain time t is p _t, the momentum vector at time t + Delta] t is p _{t + Delta] t.}

運動エネルギΔＥによって、運動量ベクトルがλ倍されるとする。つまり、

Assume that the momentum vector is multiplied by λ by the kinetic energy ΔE. That means

と表されるとする。式（３３）を式（３２）に代入することにより、λは、

It is assumed that By substituting equation (33) into equation (32), λ becomes

と表される。時間刻み幅Δｔごとに、式（３３）により運動量ベクトルを更新することにより、時間刻み幅Δｔ当たりに運動エネルギΔＥが流入するような熱の流入を模擬できる。

It is expressed. By updating the momentum vector by the equation (33) for each time step width Δt, it is possible to simulate the inflow of heat in which the kinetic energy ΔE flows per time step width Δt.

When the number of particles on the heating surfaces and N _h, the area of the heating surface and S _h, power f _h which flows per unit area,

と表される。第１のシミュレーションでは、加熱面上の粒子数Ｎ_ｈが４００個であり、加熱面の面積Ｓ_ｈが２．８３×１０^−１７[ｍ^２]であり、運動エネルギΔＥが４．１４×１０^−２４[Ｊ]であり、時間刻み幅Δｔが５．０[ｆｓ]であり、パワーｆ_ｈが１．１７×１０^１０[Ｗ／ｍ^２]である。

It is expressed. In the first simulation, the particle number _{N h} of the heating surface is 400, the area _{S h} of the heating surface is ^{2.83 × 10 -17 [m 2]} , the kinetic energy ΔE is 4.14 × 10 ⁻²⁴ [J], the time step Δt is 5.0 [fs], and the power f _h is 1.17 × 10 ¹⁰ [W / m ² ].

In the second and third simulations (in the renormalized system), the power f _h ′ flowing in per unit area is

という式で与えられる。運動量ベクトルを更新するための係数λ´は、

It is given by the formula. The coefficient λ ′ for updating the momentum vector is

という式で与えられる。なお式（３１）から、与えられる運動エネルギΔＥは元の系と繰り込まれた系とで等しくなる。

It is given by the formula. Note that from equation (31), the given kinetic energy ΔE is equal between the original system and the renormalized system.

Next, simulation results will be described with reference to FIG. In the graph shown in FIG. 2, the lower graph represents the result of the first simulation (conventional molecular dynamics calculation), and the middle graph represents the second simulation (renormalization group molecular dynamics with the renormalization factor α being 2). The upper graph represents the result of the third simulation (renormalization group molecular dynamics calculation with the renormalization factor α being 4). In the three graphs, the vertical axis indicates the temperature (unit [K]), and the horizontal axis indicates the amount of heat that flows in, that is, the energy that flows (unit 10 ⁻¹⁶ [J]).

First, the result of the first simulation will be described with reference to the lower graph. When energy is not flowing into the system, the temperature of the system is 300 [K]. As the incoming energy increases, the temperature of the system increases. When the inflow energy reaches about 5 × 10 ⁻¹⁶ [J], the temperature temporarily stops at about 850 [K], and the temperature is almost constant until the inflow energy reaches about 9 × 10 ⁻¹⁶ [J]. It becomes. This temperature indicates the melting point. When the inflow energy exceeds about 9 × 10 ⁻¹⁶ [J], the temperature rises again.

  Energy that flows into the system during a period of constant temperature corresponds to latent heat. The latent heat obtained from this simulation is 14.7 [kJ / mol]. In addition, melting | fusing point of aluminum calculated | required from experiment is 933 [K], and latent heat is 11 [kJ / mol]. The difference between the value obtained by the simulation and the experimental value is thought to be due to the accuracy of the interatomic potential used.

Next, the results of the second and third simulations will be described with reference to the middle and upper graphs. The results of the second and third simulations show the same tendency as the results of the first simulation. That is, the temperature tends to be substantially constant around 850 [K] during the period from about 5 × 10 ⁻¹⁶ [J] to about 9 × 10 ⁻¹⁶ [J]. Thereby, it can be determined from the second and third simulations that the melting point is about 850 [K] and the latent heat is 14.7 [kJ / mol].

  The energy flowing into the particle system can be obtained by integrating the inflow power for the simulated period. Assuming that the energy that has flowed into the system that has been brought in is E ′, the energy E that has flowed into the original system is

という式から求められる。なお、粒子１つ当たりの流入エネルギは、元の系と繰り込まれた系とで等しい。

It is calculated from the formula. Note that the inflow energy per particle is the same for the original system and the renormalized system.

  The temperature of the system brought in is obtained from the equation (30). From the above discussion, the temperature of the original system is equal to the temperature of the renormalized system.

  In the above-described simulation, the inflow energy and temperature in the original system were obtained, but other physical quantities in the original system, such as the position, momentum, and velocity of the particles, are also in accordance with the scale conversion law as necessary. Can be based on.

  As described above, in the renormalization group molecular dynamics, the physical quantity expected to be obtained by the molecular dynamics calculation for the system before renormalization (the system corresponding to the first simulation of the above example) (in the above example, (Melting point and latent heat) can be calculated by performing molecular dynamics calculation for a system having a smaller number of particles than the system before retraction. Various calculations required for simulation by renormalization group molecular dynamics can be executed by a computer by a program.

  In the simulation of the above embodiment, aluminum has a face-centered cubic structure. The scale conversion law is derived from the consideration for a particle system arranged on a simple cubic lattice, but it is also effective for other crystal structures that are not limited to a simple cubic structure.

Next, a description will be given of a simulation in which heat is caused to flow into silicon to raise the temperature, and the melting point and latent heat are calculated. Silicon takes a diamond structure. Similar to the simulation for aluminum, three simulations were performed. The first simulation is based on a conventional molecular dynamics method performed for comparison. In the second and third simulations, molecular dynamics calculation was performed on a system incorporating the system corresponding to the first simulation. The renormalization factor α is 2 in the second simulation and 2 ² (= 4) in the third simulation. The number of particles in the system is 8192 in the first simulation, 1024 in the second simulation, and 128 in the third simulation.

  As the interatomic potential in the first simulation (that is, the interparticle potential in the original system), SW (Stillinger-Weber) potential was used. This potential is

と表される。ここで、

It is expressed. here,

であり、相互作用係数εは、２．１６７[ｅＶ]であり、σは２．０９５[Å]であり、定数Ａ、Ｂ、ｐ、ｑ、ａ、λ、及びγは、それぞれ、７．０５、０．６０２、４．０、０．０、１．８、２１．０、及び１．２である。

The interaction coefficient ε is 2.167 [eV], σ is 2.095 [Å], and the constants A, B, p, q, a, λ, and γ are 7. 05, 0.602, 4.0, 0.0, 1.8, 21.0, and 1.2.

Further, as shown in FIG. 3A, the angle formed by the combination of the particle i and the particle k and the combination of the particle i and the particle j is θ _ijk , and the interparticle distance from the particle i to the particle j is r _ij . The renormalized interaction coefficient in the interparticle potential in the second and third simulations is obtained by the same method as in the simulation for aluminum described above.

The mass of silicon atoms (that is, the mass of particles in the original system) in the first simulation was 4.67 × 10 ⁻²⁶ [kg]. The mass with which the particles are brought into the second and third simulations is obtained by the same method as the simulation for the above-described aluminum.

  Next, simulation results will be described with reference to FIG. In the graph shown in FIG. 3B, the lower graph represents the result of the first simulation (conventional molecular dynamics calculation), and the middle graph represents the second simulation (renormalization factor α being 2). The results of the molecular dynamics calculation are shown, and the upper graph shows the results of the third simulation (renormalization group molecular dynamics calculation with the renormalization factor α being 4).

First, the result of the first simulation will be described with reference to the lower graph. During the period until the inflow energy reaches about 1.3 × 10 ⁻¹⁶ [J] to about 4.6 × 10 ⁻¹⁶ [J], the temperature becomes substantially constant at about 1750 [K]. The latent heat is determined to be 25.6 [kJ / mol]. The melting point of silicon obtained from the experiment is 1687 [K], and the latent heat is 50.6 [kJ / mol]. The difference between the value obtained by the simulation and the experimental value is thought to be due to the accuracy of the interatomic potential.

  Next, the results of the second and third simulations will be described with reference to the middle and upper graphs. The results of the second and third simulations show almost the same tendency as the results of the first simulation. Thereby, based on the 2nd and 3rd simulation, melting | fusing point and latent heat close | similar to the value obtained by the 1st simulation can be obtained.

Next, with reference to FIG. 4, the simulation which calculated | required the dispersion | distribution relationship by renormalization group molecular dynamics is demonstrated. A wave number in the [111] direction was calculated by applying vibration in parallel to the [111] direction of aluminum. The atoms were connected by springs. The spring constant s was 40.3 [N / m]. The spring constant s is obtained by approximating the minimum vicinity of the interatomic potential of aluminum with a quadratic function. The mass m of the aluminum atoms was 4.48 × 10 ⁻²⁶ [kg]. The interatomic distance a was 2.86 × 10 ⁻¹⁰ [m]. The exact solution of the dispersion relationship (analytical solution) is ω as the angular frequency of excitation and k as the wave number in the [111] direction.

という式で与えられる。なお、この厳密解は、式（２６）においてα＝１として求めることができる。
グラフ中の実線Ｌが、厳密解を示す。点Ｐ１が、従来の分子動力学計算（繰り込み因子α＝１の場合に対応）の結果を示す。点Ｐ２〜点Ｐ５が、それぞれ、繰り込み因子αを、２^７、２^１５、２^２１、及び２^２８とした場合の結果を示す。繰り込み因子αを大きくすることは、元の系（繰り込まれる前の系）のサイズを大きくすることに対応する。繰り込み因子αを大きくすることにより、小さい波数を得ることができる。繰り込み群分子動力学で得られた結果は、厳密解と良く一致している。このように、繰り込み群分子動力学を用いれば、従来の分子動力学で実行するのが困難な長波長域の計算（連続体とみなせるようなマクロな系に対する計算）を、精度良く行うことができる。

It is given by the formula. This exact solution can be obtained as α = 1 in equation (26).
A solid line L in the graph indicates an exact solution. Point P1 shows the result of the conventional molecular dynamics calculation (corresponding to the renormalization factor α = 1). Points P2 to P5 show the results when the renormalization factor α is 2 ⁷ , 2 ¹⁵ , 2 ²¹ , and 2 ²⁸ , respectively. Increasing the renormalization factor α corresponds to increasing the size of the original system (system before the renormalization). A small wave number can be obtained by increasing the renormalization factor α. The results obtained with renormalization group molecular dynamics are in good agreement with the exact solution. In this way, using renormalization group molecular dynamics, it is possible to accurately calculate long wavelength regions (calculations for macroscopic systems that can be regarded as continuum) that are difficult to perform with conventional molecular dynamics. it can.

  An XYZ orthogonal coordinate system is considered in the space where the particle system is arranged. In the above-described embodiment, the renormalized position vector q ′ using the renormalization factor α from the original system position vector q is

と求められた。これは、Ｘ、Ｙ及びＺのすべての方向に関する空間的なスケール変換が、共通の繰り込み因子αにより行われることを示す。繰り込み因子αを、Ｘ、Ｙ及びＺの各方向について独立に設定することも可能である。つまり、異方性を有するスケール変換を行うことも可能である。

I was asked. This indicates that the spatial scaling for all directions of X, Y and Z is performed by a common renormalization factor α. It is also possible to set the renormalization factor α independently in each of the X, Y, and Z directions. That is, it is possible to perform scale conversion having anisotropy.

When the second renormalization factor γ is 0, the scale conversion rule having anisotropy is represented by the following conversion equation. The renormalization factor in the X direction and alpha _X, the renormalization factor in the Y direction and alpha _Y, the renormalization factor Z direction is alpha _Z. Each of the renormalization factors α _{X to} α _Z is, for example, 2 ⁿ (n is a positive integer).

  The number N of particles in the particle system S and the number N ′ of particles carried in the particle system S ′ are converted into a conversion formula.

により対応付けられる。粒子系Ｓのある粒子の位置ベクトルｑの各成分と、繰り込まれた粒子系Ｓ´のある粒子の位置ベクトルｑ´の各成分とは、変換式

Are associated with each other. Each component of the position vector q of the particle having the particle system S and each component of the position vector q ′ of the particle having the regenerated particle system S ′ are converted into a conversion formula.

により対応付けられる。ここで、ηは、ベクトルのＸ，Ｙ、Ｚのいずれかの成分を示す（以下、（変換式ｐａ）、（変換式ｖａ）、及び（変換式ｍａ）においても同様である。）。粒子系Ｓのある粒子の運動量ベクトルｐの各成分と、繰り込まれた粒子系Ｓ´のある粒子の運動量ベクトルｐ´の各成分とは、変換式

Are associated with each other. Here, η represents one of the X, Y, and Z components of the vector (hereinafter, the same applies to (conversion equation pa), (conversion equation va), and (conversion equation ma)). Each component of the momentum vector p of the particle with the particle system S and each component of the momentum vector p ′ of the particle with the retrained particle system S ′

により対応付けられる。粒子系Ｓのある粒子の速度ベクトルｖの各成分と、繰り込まれた粒子系Ｓ´のある粒子の速度ベクトルｖ´の各成分とは、変換式

Are associated with each other. Each component of the velocity vector v of the particle having the particle system S and each component of the velocity vector v 'of the particle having the regenerated particle system S ′ are converted into a conversion formula.

により対応付けられる。

Are associated with each other.

Since a renormalization factor is set for each direction, the particle mass is also scaled for each direction. And the mass m of the particles of the particle system S, of each particle of the renormalized particle system S', X direction of the mass _{m'X, Y} direction of the mass _m'Y, and the Z direction of the mass _m'Z is , Conversion formula

により対応付けられる。繰り込まれた系における運動方程式（２５）で、方向ごとに、対応する質量が用いられる。

Are associated with each other. In the equation of motion (25) in the renormalized system, the corresponding mass is used for each direction.

  The interaction coefficient ε in the particle system S and the interaction coefficient ε ′ in the renormalized particle system S ′

により対応付けられ、粒子系Ｓにおける時間スケールｔと、繰り込まれた粒子系Ｓ´におけるスケールｔ´とは、変換式

The time scale t in the particle system S and the scale t ′ in the renormalized particle system S ′ are converted by the conversion formula

により対応付けられる。

Are associated with each other.

For example, consider a simulation for a thin film. The thin film has a larger dimension in the direction parallel to the film surface than the dimension in the film thickness direction. For this reason. For the direction parallel to the film thickness, the renormalization factor may be larger than that in the film thickness direction. For example, if the film thickness direction is the Z direction, the renormalization factors α _X and α _Y in the direction parallel to the film surface may be set larger than the renormalization factor α _Z in the film thickness direction.

Note that the second scaling parameter γ is included in (conversion equation m), (conversion equation ε), and (conversion equation t) in the form of α ^γ . When γ = 0 is set, α ^γ is simply “1”. Thus, with respect to particle mass, interaction coefficient, and time, without using the second scaling parameter γ (without introducing a second scaling parameter γ)

というスケール変換を施したとしても、第２のスケーリングパラメータγを０と設定したのと同様なスケール変換が行われる。なお、γ＝０のとき、速度ベクトルについては、

Even if the scale conversion is performed, the same scale conversion as when the second scaling parameter γ is set to 0 is performed. When γ = 0, the velocity vector is

というスケール変換が行われる。

Scale conversion is performed.

If the renormalization factors α, α _X , α _Y , and α _Z are each greater than 1, the number of particles in the renormalized system is smaller than the number of particles in the original system.

  If necessary, the renormalization factor (for example, α) can be made different between a part of the system to be simulated and another part.

  In the above-described examples, the melting point, latent heat, and dispersion relationship were obtained as examples of renormalization group molecular dynamics. Renormalization group molecular dynamics can be used for simulation of various physical phenomena, but not limited to these. Since renormalization group molecular dynamics uses molecular dynamics, it can be applied to various problems that can be handled by molecular dynamics. For example, problems with flow, structure, heat, reaction, etc., as well as problems with mechanism, wear and lubrication can be addressed. Renormalization group molecular dynamics makes it possible to obtain physical quantities for systems with a larger number of particles than systems that actually perform molecular dynamics calculations. Therefore, if the renormalization group molecular dynamics discovered by the inventors of the present application are used, it is possible to perform a highly accurate simulation with respect to a scale problem that has been difficult to perform in the conventional simulation. In addition, as a method of molecular dynamics calculation performed in renormalization group molecular dynamics, various well-known things can be used.

  In Non-Patent Document 1 described in the “Background Art” column of this specification, an attempt is made to apply the renormalization technique to molecular dynamics. However, in Non-Patent Document 1, only the interaction part of the Hamiltonian is renormalized, and the kinetic energy part is not renormalized. In Non-Patent Document 1, a simple average is taken with respect to the momentum (or speed) instead of renormalization. For this reason, it is limited to static analysis, the temperature cannot be taken into account explicitly, and an artificial method for forcibly causing energy dissipation (see “3. Numerical calculation method” in Non-Patent Document 1). (23), (24), and (25)) described in (1) are required. For this reason, there are problems such as lack of reliability when applied to dynamic problems involving calculation of finite temperature, vibration / noise, and phase transition. The scale conversion rule disclosed in Non-Patent Document 1 is as follows.

ここで、粒子の位置ベクトルがｑであり、運動量ベクトルがｐであり、質量がｍである。粒子間ポテンシャルの相互作用係数がεであり、時間がｔである。繰り込まれた系の物理量にプライム「´」を付す。αが繰り込み因子である。

Here, the position vector of the particle is q, the momentum vector is p, and the mass is m. The interaction coefficient of the interparticle potential is ε, and the time is t. The prime “′” is added to the physical quantity of the system that has been transferred. α is a renormalization factor.

  In the process of deriving the scale conversion law of the present inventors, the kinetic energy part of the Hamiltonian is also renormalized. This leads to an appropriate scale conversion law for the momentum vector (or velocity vector). Therefore, it is possible to explicitly consider the temperature in the simulation.

  Although the present invention has been described with reference to the embodiments, the present invention is not limited thereto. It will be apparent to those skilled in the art that various modifications, improvements, combinations, and the like can be made.

FIG. 1A and FIG. 1B are diagrams for explaining a method of deriving a scale conversion rule. It is a graph which shows the result of the simulation with respect to aluminum using renormalization group molecular dynamics. FIG. 3A is a diagram for explaining the SW potential, and FIG. 3B is a graph showing the result of simulation for silicon using renormalization group molecular dynamics. It is a graph which shows the result of the simulation which calculated | required the dispersion | distribution relationship using renormalization group molecular dynamics.

Explanation of symbols

α First renormalization factor N Number of particles N ′ Number of renormalized particles

Claims (11)

(A) including N particles, the mass of each particle is m, and the interaction potential energy between the particles is expressed by a product εf of a non-dimensionalization function f and an interaction coefficient ε representing the dependence on the interparticle distance. When considering a particle system S that can be represented, the dimension number d of the space in which the particle system S is arranged, a first renormalization factor α greater than 1, and a second renormalization factor γ greater than or equal to 0 and less than or equal to d. And the number N ′ of renormalized particles is obtained by the conversion formula N ′ = N / α ^d , and the mass m ′ of the retreated particles is obtained by the conversion formula m ′ = mα ² / α ^γ. , The step of determining the renormalized interaction coefficient ε ′ by the conversion equation ε ′ = εα ^γ ,
(B) including particles having a number N ′ of renormalized particles, each particle having a mass m ′ of the renormalized particles, and the interaction potential energy between the particles is the dimensionless function f And a step of executing molecular dynamics calculation for the particle system S ′ represented by the product ε′f with the renormalized interaction coefficient ε ′. further,
(C) A position vector obtained by molecular dynamics calculation in the step (b) of a particle contained in the particle system S ′ is represented as q ′, and the particle contained in the particle system S ′ The momentum vector obtained by the molecular dynamics calculation in the step (b) is represented as p ′, and the velocity vector obtained by the molecular dynamics calculation in the step (b) of a certain particle included in the particle system S ′ is represented by When expressed as v ′ and a certain time interval in the molecular dynamics calculation in the step (b) of a particle included in the particle system S ′ is expressed as t ′, the first renormalization factor α and the second Calculation using a renormalization factor γ to obtain a position vector q of a certain particle contained in the particle system S by a conversion equation q = q′α, and a momentum vector p of a certain particle contained in the particle system S Total calculated by p = p ′ / α In calculation, a velocity vector v of a particle contained in the particle system S is calculated by a conversion formula v = v′α, and a certain time interval t in the molecular dynamics calculation assumed to be performed on the particle system S is of the calculation for obtaining the conversion equation t = t'α ^γ, the simulation method according to claim 1 comprising the step of performing at least one calculation.   The simulation method according to claim 1, wherein the second renormalization factor γ is zero. The simulation method according to claim 1, wherein the first renormalization factor is 2 ⁿ where n is a positive integer.   The simulation method according to claim 1, wherein the dimension number d is three. (A) including N particles, the mass of each particle is m, and the interaction potential energy between the particles is expressed by a product εf of a non-dimensionalization function f and an interaction coefficient ε representing the dependence on the interparticle distance. When a particle system S that can be expressed is considered, the number of renormalized particles N ′ is converted by using a dimension d of the space in which the particle system S is arranged and a renormalization factor α larger than 1. determined by the equation N'= N / α ^d, the mass m'of renormalized particles, a step of determining the conversion equation m'= m.alpha ^2,
(B) including particles having a number N ′ of renormalized particles, each particle having a mass m ′ of the renormalized particles, and the interaction potential energy between the particles is the dimensionless function f And a step of performing molecular dynamics calculation for the particle system S ′ represented by the product εf with the interaction coefficient ε. further,
(C) A position vector obtained by molecular dynamics calculation in the step (b) of a particle contained in the particle system S ′ is represented as q ′, and the particle contained in the particle system S ′ The momentum vector obtained by the molecular dynamics calculation in the step (b) is represented as p ′, and the velocity vector obtained by the molecular dynamics calculation in the step (b) of a certain particle included in the particle system S ′ is represented by When expressed as v ′, using the renormalization factor α, the position vector q of a certain particle contained in the particle system S is calculated by a conversion equation q = q′α, and the calculation of the certain particle contained in the particle system S At least one of a calculation for obtaining the momentum vector p by the conversion formula p = p ′ / α and a calculation for obtaining the velocity vector v of a particle contained in the particle system S by the conversion formula v = v′α. Have a process to execute Simulation method according to claim 6 that. (D) includes N particles, the mass of each particle is m, and the interaction potential energy between the particles is expressed by the product εf of the non-dimensionalization function f and the interaction coefficient ε representing the dependence on the interparticle distance. When an XYZ rectangular coordinate system is considered in the space where the particle system S can be represented and the space in which the particle system S is arranged, the renormalization factor α _{X in} the X direction and the renormalization factor α _Y in the Y direction, each greater than 1. And the renormalization factor α _Z in the Z direction, the number N ′ of renormalized particles is obtained by the conversion formula N ′ = N / (α _X α _Y α _Z ), and the mass of the renormalized particles in the X direction m _X ′ is obtained by the conversion formula m _X ′ = mα _X ² , and the mass m _Y ′ of the renormalized particle in the Y direction is obtained by the conversion formula m _Y ′ = mα _Y ^{2, and the} retraction in the Z direction is performed. The mass m _Z ′ of the particles is converted into the conversion formula m _Z ′ = mα _Z ^The process determined by ² .
(E) including particles having a number N ′ of renormalized particles, each particle having a mass m ′ of the renormalized particles, and the interaction potential energy between the particles is the dimensionless function f For the particle system S ′ represented by the product εf with the interaction coefficient ε, the renormalized particle mass m _X ′ in the _X direction is used for the X direction motion equation, and the Y direction motion equation is Using molecular mass m _Y ′ renormalized in the Y direction and using the mass m _Z ′ in the Z direction for the equation of motion in the Z direction, and performing a molecular dynamics calculation. Having a simulation method. further,
(F) A position vector of a certain particle contained in the particle system S ′ obtained by molecular dynamics calculation in the step (e) is represented by q ′, and each of the position vector q ′ in the X, Y, and Z directions Are represented as q _X ′, q _Y ′, and q _Z ′, and a momentum vector obtained by molecular dynamics calculation in the step (e) of a particle contained in the particle system S ′ is represented by p ′. And the components in the X, Y, and Z directions of the momentum vector p ′ are represented as p _X ′, p _Y ′, and p _Z ′. The velocity vector obtained by the molecular dynamics calculation in e) is represented by v ′, and the components of the velocity vector v ′ in the X, Y and Z directions are represented by v _X ′, v _Y ′ and v _Z ′. when the X-direction of the renormalization factor alpha _X, Y direction renormalization factor alpha _Y, Using fine Z direction renormalization factor alpha _Z, X of the position vector q of a particle contained in the particles based S, Y and Z directions of the component q _X, q _Y, and q _Z, respectively, conversion formula _{q X = q X'α X,} q Y = q Y'α Y, and _{q Z} = _{q Z'α} Z calculation for obtaining the, X momentum vector p of a particle included in the particle system S, Y and Z The components p _X , p _Y , and p _Z in the respective directions are converted into the transformation formulas p _X = p _X '/ α _X , p _Y = p _Y ' / α _Y , and p _Z = p _Z '/ α _Z , respectively. Calculations to be obtained and components v _X , v _Y , and v _Z in the X, Y, and Z directions of a velocity vector v of a particle included in the particle system S are respectively converted into transformation equations v _X = v _{X ′} α. _{X, v Y = v Y'α} Y, and _v Z = _v of calculation for determining the _{Z'α Z,} small The simulation method according to claim 8, further comprising a step of executing at least one calculation. The simulation method according to claim 8 or 9, wherein at least two of the renormalization factor α _{X in} the X direction, the renormalization factor α _{Y in} the Y direction, and the renormalization factor α _{Z in} the Z direction are different from each other. The program for performing the simulation method in any one of Claims 1-10 with a computer.

Images