WO2018198273A1 - シミュレーション方法、シミュレーションプログラム、及びシミュレーション装置 - Google Patents

シミュレーション方法、シミュレーションプログラム、及びシミュレーション装置 Download PDF

Info

Publication number
WO2018198273A1
WO2018198273A1 PCT/JP2017/016774 JP2017016774W WO2018198273A1 WO 2018198273 A1 WO2018198273 A1 WO 2018198273A1 JP 2017016774 W JP2017016774 W JP 2017016774W WO 2018198273 A1 WO2018198273 A1 WO 2018198273A1
Authority
WO
WIPO (PCT)
Prior art keywords
particle system
represented
potential
monomer particles
particle
Prior art date
Application number
PCT/JP2017/016774
Other languages
English (en)
French (fr)
Inventor
義崇 小林
大路 市嶋
Original Assignee
住友重機械工業株式会社
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by 住友重機械工業株式会社 filed Critical 住友重機械工業株式会社
Priority to PCT/JP2017/016774 priority Critical patent/WO2018198273A1/ja
Publication of WO2018198273A1 publication Critical patent/WO2018198273A1/ja
Priority to US16/450,341 priority patent/US11164662B2/en

Links

Images

Classifications

    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16C—COMPUTATIONAL CHEMISTRY; CHEMOINFORMATICS; COMPUTATIONAL MATERIALS SCIENCE
    • G16C10/00—Computational theoretical chemistry, i.e. ICT specially adapted for theoretical aspects of quantum chemistry, molecular mechanics, molecular dynamics or the like
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16C—COMPUTATIONAL CHEMISTRY; CHEMOINFORMATICS; COMPUTATIONAL MATERIALS SCIENCE
    • G16C60/00—Computational materials science, i.e. ICT specially adapted for investigating the physical or chemical properties of materials or phenomena associated with their design, synthesis, processing, characterisation or utilisation
    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16Z—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS, NOT OTHERWISE PROVIDED FOR
    • G16Z99/00—Subject matter not provided for in other main groups of this subclass
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/10—Numerical modelling

Definitions

  • the present invention relates to a simulation method, simulation program, and simulation apparatus based on molecular dynamics applying renormalization groups.
  • Patent Document 1 discloses a simulation method using a renormalization transformation method in order to reduce a necessary calculation amount.
  • the renormalization conversion method disclosed in Patent Document 1 will be described.
  • the number of particles of the particle system S to be simulated is represented by N, the mass of each particle is represented by m, and the interaction potential between the particles is represented by ⁇ (r).
  • r is the distance between the particles.
  • the interaction potential ⁇ (r) is represented by the product of the interaction coefficient ⁇ and the function f (r).
  • the interaction coefficient ⁇ represents the strength of interaction and has an energy dimension.
  • the function f (r) represents the dependence on the interparticle distance and is dimensionless.
  • a first renormalization factor ⁇ , a second renormalization factor ⁇ , and a third renormalization factor ⁇ are determined.
  • the first renormalization factor ⁇ is greater than 1.
  • the second renormalization factor ⁇ is not less than 0 and not more than the dimension number d of the space.
  • the third renormalization factor ⁇ is 0 or more.
  • Molecular dynamics calculation is performed on the particle system S ′ subjected to the renormalization transformation.
  • the position vector of each particle obtained by molecular dynamics calculation is represented by q ′, and the momentum vector is represented by p ′.
  • the distance between the monomers of the particle system S and the particle system S ′ is represented by r, the number of dimensions in space is represented by d,
  • the mass of the monomer of the particle system S is represented by m, the number of monomer particles constituting one polymer is represented by Nm, the number of polymers is represented by Np, the coordinates of the monomer particles are represented by q, and between the monomer particles As interaction potential ⁇ (r), Between any monomer particles, And the potential U 0 (r) is
  • the finite stretch nonlinear elastic potential U ch (r: ⁇ , ⁇ ) having the parameters ⁇ , ⁇ used as parameters as the potential U 0 (r) is the potential U 0 (r Added to Is applied,
  • the mass of the monomer particles in the particle system S ′ is represented by m ′
  • the number of monomer particles constituting one polymer is represented by Nm ′
  • the number of polymers is represented by Np ′
  • the coordinates of the monomer particles are represented by q ′.
  • a computer program for executing the simulation a recording medium storing the computer program, and a simulation apparatus are provided.
  • FIG. 1 is a flowchart of a simulation method according to an embodiment.
  • FIG. 2 is a schematic diagram showing a state in which particles i, particles j, and particles k are arranged one-dimensionally in this order.
  • FIG. 3 is a graph showing the calculation result of Expression (27) and the numerical integration result of Expression (28) when the interaction potential ⁇ (r) is a Leonard-Jones type potential.
  • FIGS. 4A to 4C are diagrams for explaining a potential moving method for reducing a two-dimensional lattice to a one-dimensional lattice.
  • FIG. 5 is a graph showing a calculation result of 2 ⁇ ((q i ⁇ q k ) / 2) and a numerical integration result of ⁇ .
  • FIG. 1 is a flowchart of a simulation method according to an embodiment.
  • FIG. 2 is a schematic diagram showing a state in which particles i, particles j, and particles k are arranged one-dimensionally in this order.
  • FIG. 3
  • FIG. 6 is a schematic diagram showing an example in which linear polymers are regularly arranged in a three-dimensional direction.
  • FIG. 7 is a graph showing the results of determining the strain rate dependence of the viscosity of the polymer material by simulation.
  • FIG. 8 is a block diagram of a simulation apparatus according to the embodiment.
  • the momentum vector p of each particle at each time is solved by numerically integrating and solving the equations of motion represented by equations (7) and (8) for each particle constituting the particle system.
  • j and the position vector q j are determined.
  • the Verlet method is often used for numerical integration. The Verlet method is described, for example, on page 175 of JMThijssen, “Computational Physics”, Cambridge University Press (1999).
  • Various physical quantities of the particle system can be calculated based on the momentum vector and position vector of each particle obtained by molecular dynamics calculation.
  • renormalization group molecular dynamics applying the renormalization group technique (hereinafter referred to as renormalization group molecular dynamics) will be explained conceptually.
  • the particle system S to be simulated corresponds to a particle system S ′ (hereinafter referred to as a renormalized particle system S ′) composed of a smaller number of particles than the particle system S. Put it on.
  • a molecular dynamics calculation is performed on the renormalized particle system S ′.
  • the calculation result for the renormalized particle system S ′ is associated with the particle system S to be simulated. Thereby, the amount of calculation can be reduced as compared with the case where the molecular dynamics calculation is directly executed on the particle system S to be simulated.
  • a conversion rule for associating a physical quantity (for example, the number of particles, a mass of a particle, etc.) in the particle system S to be simulated with a physical quantity in the renormalized particle system S ′ is called renormalization measurement.
  • FIG. 1 shows a flowchart of a simulation method according to the embodiment.
  • a renormalized particle system S ' is defined by executing renormalization conversion processing on the particle system S to be simulated based on renormalization conversion measurement.
  • the particle system S is composed of a plurality of polymers, and each of the plurality of polymers is composed of a plurality of monomer particles connected in a one-dimensional manner.
  • ⁇ (r) between the monomer particles of the particle system S between any monomer particles, Applies.
  • r represents the distance from each particle constituting the particle system S.
  • the potential U 0 (r) is basically represented by the following equation. here, f represents a dimensionless function, and ⁇ and ⁇ are parameters that characterize the monomer particles.
  • the parameter ⁇ has a dimension of energy and is called an interaction coefficient.
  • the parameter ⁇ has a distance dimension and depends on the size of the particle.
  • a finite stretch nonlinear elastic potential U ch (r: ⁇ , ⁇ ) is added to the potential U 0 (r). Applies.
  • the finite stretch nonlinear elastic potential U ch (r: ⁇ , ⁇ ) includes parameters that depend on parameters ⁇ and ⁇ that define the potential U 0 (r).
  • a Leonard-Jones type potential can be applied as the potential U 0 (r).
  • the potential U 0 (r) when applying the potential of Lennard-Jones type of shifting, the potential U 0 (r), for example can be defined by the following equation.
  • step S1 the following renormalization measurement is applied.
  • d represents the number of dimensions of the space in which the particle system S to be simulated is arranged.
  • ⁇ is a renormalization factor that depends on the number of renormalizations.
  • the renormalization factor ⁇ is expressed by the following equation.
  • the interaction potential ⁇ ′ (r) applied between any monomer particles in the particle system S ′ is defined by the following equation.
  • the number of particles becomes 1 / ⁇ d times and the mass of the particles becomes ⁇ d times. For this reason, the total mass of the particle system S and the total mass of the particle system S ′ are the same. Furthermore, since the position vector q ′ of the monomer particle in the particle system S ′ is the same as the position vector q of the monomer particle in the particle system S, the macroscopic dimension between the particle system S and the revolved particle system S ′ is Are the same. Since the total mass and size of the particle system are unchanged before and after the renormalization conversion treatment, the density of the particle system is also unchanged.
  • step S2 initial conditions for simulation are set.
  • the initial conditions include initial values of the position vector q j and the momentum vector p j of each particle.
  • the momentum vector p j is set based on the temperature T ′ of the renormalized particle system S ′.
  • T the temperature of the particle system S to be simulated
  • the following renormalization measurement is applied with respect to the temperature.
  • step S3 molecular dynamics calculation is performed for the renormalized particle system S ′.
  • the equation of motion is expressed by the following equation.
  • N ′ represents the total number of monomer particles in the particle system S ′
  • N ′ Nm ′ ⁇ Np ′.
  • the correspondence relationship between the position vector q ′ of each particle of the particle system S ′ that has been carried in and the position vector q of the particle system S to be simulated is as shown in Expression (13).
  • the momentum vector p ′ in the renormalized particle system S ′ and the momentum vector p in the particle system S to be simulated have the following relationship.
  • step S4 the simulation result is output.
  • the position vector q ′ and the momentum vector p ′ may be directly output as numerical values, or based on the position vector q ′, the distribution in the space of a plurality of particles in the particle system S ′ may be displayed as an image. Also good.
  • the Hamiltonian H ′ of the particle system S ′ In order to obtain the Hamiltonian H ′ of the particle system S ′, a part of the integration of the partition function Z ( ⁇ ) with respect to the particle system S is executed, and the Hamiltonian H ′ is obtained by coarse-graining the Hamiltonian. .
  • the distribution function Z ( ⁇ ) for a canonical ensemble with a fixed number of particles is expressed by the following equation.
  • k B is a Boltzmann constant
  • d ⁇ N is a volume element in the phase space.
  • d [gamma] N can be expressed by the following equation.
  • h is a Planck constant.
  • W N is determined so that the quantum sum of all states matches the integral over the phase space.
  • particles i, particles j, and particles k are arranged in this order in a one-dimensional manner.
  • Coarse coarsening of the interaction potential can be performed by writing out the interaction involving the particle j and executing integration with respect to the position coordinate of the particle j in the middle between the particle i and the particle k.
  • the method of potential moving can be used. The method of potential moving is described in Leo P. Kadanoff, STATISTICAL PHYSICS Static, Dynamics and Renormalization, Chap. 14, World Scientific (1999).
  • the interaction potential ⁇ tilde that incorporates the contribution from the second or more adjacent particles can be expressed by the following equation.
  • a represents the interparticle distance in an equilibrium state.
  • the interparticle distance a in the equilibrium state is equal to the distance at which the interaction potential ⁇ (r) is minimized.
  • the position vector q j of the particle j can be represented by one-dimensional coordinates q j .
  • the particles nearest is expressed by the following equation.
  • Equation (25) If integration is performed on the integration variable q j , the following equation is obtained using equation (25).
  • r a denotes the diameter of a particle
  • z (q i -q k) and P (q i -q k) is expressed by the following equation.
  • the integration region is limited to the inner region of the cage potential.
  • equation (29) Substituting equation (29) into equation (28) for numerical integration. When substituting equation (29) into equation (28), equation (24) is used.
  • FIG. 3 shows the calculation result of Expression (27) and the numerical integration result of Expression (28) when the interaction potential ⁇ (r) is a Leonard-Jones type potential.
  • Particles i, particle j, particles k, particles l, and the particle m position coordinates of the particle k when disposed one-dimensionally in this order is represented by q k.
  • the horizontal axis in FIG. 3 represents q k / 2
  • the vertical axis represents P (q i ⁇ q k ) P (q k ⁇ q m ) and z (q i ⁇ q k ) z (q k ⁇ q m ) is represented on a logarithmic scale.
  • ⁇ / k B T 2.0.
  • the probability p (q k ) that the particle k exists at the position coordinate q k can be approximated as follows.
  • the interaction potential of a multidimensional particle system can be realized by the potential moving method.
  • particles are arranged at positions of lattice points of a two-dimensional square lattice.
  • the solid line shows the interaction between nearest particles (nearest neighbor coupling).
  • One direction in which the particles are arranged is defined as an x direction, and a direction perpendicular to the direction is defined as a y direction.
  • the nearest neighbor coupling between integrable particles is distributed to particles arranged on both sides in the x direction.
  • the interparticle interaction (double coupling) to which the distributed interaction is added is represented by a double line.
  • the double coupling represented by the double line is twice as strong as the original nearest neighbor coupling.
  • the procedure for converting the two-dimensional lattice into a one-dimensional chain may be repeated for the three directions of the x direction, the y direction, and the z direction. In this way, it is possible to perform coarse graining of the particle system constituting the multidimensional lattice.
  • Coarse graining of the interaction potential of a particle system constituting a multidimensional lattice having a dimension number d is expressed by the following equation.
  • ⁇ i, j> means that a sum is taken between the nearest grids.
  • equation (33) and equation (34) By substituting equation (33) and equation (34) into equation (22), the following equation is obtained except for coefficients that do not affect the result.
  • the coarse-grained Hamiltonian (the Hamiltonian of the renormalized particle system S ′) H ′ is represented by the following equation.
  • a list of coupling constants when coarse-graining the Hamiltonian is represented by K.
  • the coupling constant list K is expressed as follows.
  • the renormalization transformation R is defined as follows.
  • the Hamiltonian Hn that has been subjected to n renormalization transformations is expressed by the following equation.
  • the momentum vector p j ′ in the renormalized particle system S ′ is expressed by the following equation.
  • Equation (13) the renormalization transformation measurement shown in Equation (13) is derived.
  • the cage potential for the particle j formed by the particles i and k is the same as in the equation (25). It can be expressed as.
  • the integration variable q j is integrated from the position where the monomer particle j contacts one monomer particle i to the position where it contacts the other monomer particle k. Denoting the diameter of the monomer particles r a, the following equation is obtained.
  • the interaction potential after renormalization is a superposition of a function of 2 ⁇ ((qi ⁇ qk) / 2) and ⁇ .
  • the finite elongation nonlinear elastic potential U ch (r: ⁇ , ⁇ )
  • the following potential known as a cremagres model can be adopted.
  • the parameters k and R 0 are determined depending on the parameters ⁇ and ⁇ appearing in the equation (10).
  • the parameters k and R 0 are defined by the following equations.
  • the interaction potential ⁇ (r) is It is expressed.
  • the potential U 0 (r) is as shown in the above equation (12). 2 ⁇ ((q i ⁇ q k ) / 2) in the equation (44) was actually calculated, and ⁇ was obtained by numerical integration.
  • FIG. 5 shows the calculation result of 2 ⁇ ((q i ⁇ q k ) / 2) and the numerical integration result of ⁇ .
  • the horizontal axis represents (q i ⁇ q k ) / 2, and the vertical axis represents energy in the unit “J”.
  • the solid line in FIG. 5 indicates the calculation result of 2 ⁇ ((q i ⁇ q k ) / 2), and the broken line indicates the numerical integration result of ⁇ . From the calculation result of FIG. 5, it can be seen that the change amount of ⁇ is sufficiently smaller than the change amount of 2 ⁇ ((q i ⁇ q k ) / 2). Accordingly, it can be approximated that ⁇ is a constant with respect to 2 ⁇ ((q i ⁇ q k ) / 2).
  • Formula (51) may be applied when performing the simulation.
  • the variable r of the interaction potential ⁇ (r) of the particle system S is replaced with r / ⁇ to obtain the potential ⁇ (r / ⁇ ), and the result is multiplied by ⁇ d, whereby the particle system S ′ Can be numerically calculated.
  • FIG. 6 shows a schematic diagram of a state in which linear polymers are regularly arranged in a three-dimensional direction. Monomer particles constituting each of the polymers are connected in the y direction. A plurality of polymers are arranged in the x, y, and z directions.
  • every other monomer particle aligned in the y direction is erased.
  • the monomer particles arranged on the surface indicated by the arrow Py are erased.
  • the number of monomer particles in the polymer is halved.
  • every other polymer arranged in the x direction is erased.
  • monomer particles arranged on the surface indicated by the arrow Px are removed.
  • the number of polymers is halved.
  • every other polymer arranged in the z direction is erased. For example, monomer particles arranged on the surface indicated by the arrow Pz are removed. This further doubles the number of polymers.
  • ⁇ , ⁇ , k, and R 0 are parameters appearing in the equations (12) and (45).
  • the parameters k and R 0 are determined by ⁇ and ⁇ of the Leonard-Jones type potential (Equation (12)).
  • K ′ and R 0 ′ in the particle system S ′ after renormalization are determined by ⁇ ′ and R 0 ′ after renormalization transformation.
  • Figure 7 shows the simulation results.
  • the horizontal axis represents the strain rate
  • the vertical axis represents the viscosity.
  • the strain rate and the viscosity are expressed in a non-dimensional manner.
  • the round symbol in FIG. 7 shows the result of the simulation for the particle system S that is not subjected to the renormalization transformation
  • the triangular symbol shows the result of the simulation for the particle system S ′ that has been subjected to the renormalization transformation. It can be seen that the viscosity obtained for the particle system S ′ after renormalization is on a curve connecting the viscosities obtained for the particle system S before renormalization. Thus, the strain rate dependence of the viscosity is maintained even in the particle system S ′ subjected to renormalization.
  • the simulation results shown in FIG. 7 show that it is possible to perform macroscopic calculation of a polymer material using the simulation method using the renormalization transformation rule according to the above embodiment.
  • the simulation method according to the above embodiment is realized by a computer executing a computer program.
  • the computer program is provided, for example, in a state recorded in a data recording medium.
  • the computer program can be provided via a telecommunication line.
  • FIG. 8 shows a block diagram of the simulation apparatus 25 according to the embodiment.
  • the simulation device 25 includes an input device 21, a processing unit 10, and an output device 22.
  • the processing unit 10 includes an input unit 11, a calculation unit 12, an output unit 13, and a storage device 14.
  • the various conditions of the simulation are input to the input unit 11 through the input device 21.
  • the dimension number d of the space which is a parameter characterizing the particle system S to be simulated
  • the mass m of the monomer particles the number Nm of monomer particles constituting one polymer
  • the value of the number Np of polymers and the particle system S initial conditions are entered,
  • the storage device 14 stores a computer program for executing a simulation.
  • the computing unit 12 executes a simulation process by executing a computer program stored in the storage device 14 based on the data input to the input unit 11.
  • the simulation result is output to the output device 22 via the output unit 13.
  • the calculation unit 12 calculates the parameters of the renormalized particle system S ′ by performing renormalization conversion processing on the input parameters of the particle system S. A molecular dynamics calculation is performed on the particle system S ′ using the calculated parameters. Since the molecular dynamics calculation is performed on the renormalized particle system S ′, the calculation time can be shortened.

Landscapes

  • Engineering & Computer Science (AREA)
  • Theoretical Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Computing Systems (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Health & Medical Sciences (AREA)
  • General Health & Medical Sciences (AREA)
  • Spectroscopy & Molecular Physics (AREA)
  • General Physics & Mathematics (AREA)
  • General Engineering & Computer Science (AREA)
  • Geometry (AREA)
  • Evolutionary Computation (AREA)
  • Computer Hardware Design (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

くりこみ回数に依存するくりこみ因子λに基づいて、シミュレーション対象の粒子系Ｓに対してくりこみ変換処理を行う。くりこまれた粒子系Ｓ'について分子動力学計算を実行することにより、粒子系Ｓ'の粒子の位置ベクトル及び運動量ベクトルを求める。粒子系Ｓは、各々が一次元状に繋がった複数のモノマー粒子からなる複数のポリマーを含む。ポリマー内の隣接するモノマー粒子間には、有限伸長非線形弾性ポテンシャルが適用される。

Description

シミュレーション方法、シミュレーションプログラム、及びシミュレーション装置
　本発明は、くりこみ群を応用した分子動力学によるシミュレーション方法、シミュレーションプログラム、及びシミュレーション装置に関する。
　分子動力学を用いたコンピュータシミュレーションが行われている。分子動力学では、シミュレーションの対象となる系を構成する粒子の運動方程式が数値的に解かれる。シミュレーション対象の系に含まれる粒子数が増加すれば、必要な計算量が増大する。現行のコンピュータの演算能力でシミュレーション可能な系の粒子数は、通常、数十万個程度である。
　下記の特許文献１に、必要な計算量を減らすためにくりこみ変換の手法を用いたシミュレーション方法が開示されている。以下、特許文献１に開示されたくりこみ変換の手法について説明する。
　シミュレーション対象の粒子系Ｓの粒子数をＮ、各粒子の質量をｍ、粒子間の相互作用ポテンシャルをφ（ｒ）で表す。ここで、ｒは、粒子間の距離である。相互作用ポテンシャルφ（ｒ）が、相互作用係数εと関数ｆ（ｒ）との積で表される。相互作用係数εは、相互作用の強さを表し、エネルギの次元を持つ。関数ｆ（ｒ）は、粒子間距離に対する依存性を表し、無次元化されている。
　第１のくりこみ因子α、第２のくりこみ因子γ、及び第３のくりこみ因子δが決定される。第１のくりこみ因子αは１より大きい。第２のくりこみ因子γは０以上であり、かつ空間の次元数ｄ以下である。第３のくりこみ因子δは０以上である。くりこみの回数をｎで表したとき、第１のくりこみ因子αは、
　α＝２ｎ
と表される。
　くりこみ手法を用いてくりこみ変換された粒子系Ｓ’の粒子数をＮ’、各粒子の質量をｍ’、相互作用係数をε’で表す。以下の変換式を用いて、くりこみ変換された粒子系Ｓ’の粒子数Ｎ’、質量ｍ’、及び相互作用係数ε’を求める。
　Ｎ’＝Ｎ／αｄ
　ｍ’＝ｍαδ－γ
　ε’＝εαγ
　くりこみ変換された粒子系Ｓ’について、分子動力学計算を行なう。分子動力学計算で求められた各粒子の位置ベクトルをｑ’、運動量ベクトルをｐ’で表す。元の粒子系Ｓの各粒子の位置ベクトルｑ、運動量ベクトルｐは、以下の式を用いて算出することができる。
　ｑ＝ｑ’α
　ｐ＝ｐ’／αδ／２
特許第５２４１４６８号公報
　粒子系Ｓが単原子の集合体である場合、特許文献１に開示されたくりこみ変換により、くりこまれた粒子系Ｓ’に対して分子動力学計算を行うことができる。多数のポリマー分子からなる粒子系Ｓに対して分子動力学計算を行う手法として、クレマーグレストモデル（Kremer Grest Model）が提案されている。クレマーグレストモデルにくりこみ変換処理を適用してシミュレーションを行う手法は確立されていない。
　本発明の目的は、ポリマー分子からなる粒子系に対してくりこみ群の考え方を適用し、ポリマー分子の挙動をシミュレーションする方法を提供することである。本発明の他の目的は、このシミュレーション方法を実行するコンピュータプログラムを提供することである。本発明のさらに他の目的は、このシミュレーション方法を実行するシミュレーション装置を提供することである。
　本発明の一観点によると、
　くりこみ回数に依存するくりこみ因子λに基づいて、各々が一次元状に繋がった複数のモノマー粒子からなる複数のポリマーで構成されたシミュレーション対象の粒子系Ｓに対してくりこみ変換処理を行う工程と、
　くりこまれた粒子系Ｓ’について分子動力学計算を実行することにより、くりこまれた粒子系Ｓ’のモノマー粒子の位置ベクトル及び運動量ベクトルを求める工程と
を有し、
　粒子系Ｓ及び粒子系Ｓ’のモノマー間の距離をｒで表し、空間の次元数をｄで表し、
　粒子系Ｓのモノマーの質量をｍで表し、１つのポリマーを構成しているモノマー粒子の数をＮｍで表し、ポリマーの数をＮｐで表し、モノマー粒子の座標をｑで表し、モノマー粒子間の相互作用ポテンシャルφ（ｒ）として、
　任意のモノマー粒子間においては、
Figure JPOXMLDOC01-appb-M000017
 が適用され、ポテンシャルＵ０（ｒ）は、エネルギの次元を持つパラメータεと、長さの次元を持つパラメータσと、無次元関数ｆを用いて
Figure JPOXMLDOC01-appb-M000018
 と表すことができ、
　同一ポリマー内の隣接するモノマー粒子間においては、パラメータとしてポテンシャルＵ０（ｒ）で用いられたパタメータε、σを持つ有限伸長非線形弾性ポテンシャルＵｃｈ（ｒ：ε，σ）がポテンシャルＵ０（ｒ）に追加されて
Figure JPOXMLDOC01-appb-M000019
 が適用され、
　粒子系Ｓ’のモノマー粒子の質量をｍ’で表し、１つのポリマーを構成しているモノマー粒子の数をＮｍ’で表し、ポリマーの数をＮｐ’で表し、モノマー粒子の座標をｑ’で表したとき、下記の変換則
Figure JPOXMLDOC01-appb-M000020
 を適用し、
　粒子系ＳにおけるポテンシャルＵ０（ｒ）、Ｕｃｈ（ｒ：ε，σ）に対応する粒子系Ｓ’におけるポテンシャルＵ０’（ｒ）、Ｕｃｈ’（ｒ：ε，σ）として、
Figure JPOXMLDOC01-appb-M000021
 を用いて、粒子系Ｓ’について前記分子動力学計算を実行するシミュレーション方法が提供される。
　本発明の他の観点によると、上記シミュレーションを実行するコンピュータプログラム、コンピュータプログラムを記録した記録媒体、及びシミュレーション装置が提供される。
　複数のモノマー粒子が鎖状に繋がったポリマーからなる粒子系にくりこみ変換を適用して、分子動力学計算を行うことができる。これにより、計算時間を短縮することができる。
図１は、実施例によるシミュレーション方法のフローチャートである。 図２は、粒子ｉ、粒子ｊ、及び粒子ｋがこの順番に一次元状に並んで配置されている状態を示す模式図である。 図３は、相互作用ポテンシャルφ（ｒ）がレナード・ジョーンズ型ポテンシャルである場合の式（２７）の算出結果、及び式（２８）の数値積分結果を示すグラフである。 図４Ａ～図４Ｃは、二次元格子を一次元格子に帰着させるポテンシャルムービングの方法について説明するための図である。 図５は、２φ（（ｑｉ－ｑｋ）／２）の計算結果、及びδφの数値積分結果を示すグラフである。 図６は、直線状のポリマーが３次元方向に規則正しく配置されている例を示す模式図である。 図７は、高分子材料の粘度のひずみ速度依存性をシミュレーションによって求めた結果を示すグラフである。 図８は、実施例によるシミュレーション装置のブロック図である。
　本願発明の実施例で適用される分子動力学（Molecular Dynamics）について簡単に説明する。Ｎ個の粒子（例えば原子）からなり、ハミルトニアンＨが下記の式で表される粒子系について考える。
Figure JPOXMLDOC01-appb-M000022
 ここで、
ｍは粒子の質量を表し、
φは粒子間の相互作用ポテンシャルを表し、
ベクトルｐｊは粒子の運動量ベクトルを表し、
ベクトルｑｊは粒子の位置ベクトル（位置座標）を表す。
　ハミルトニアンＨを、ハミルトニアンの正準方程式に代入することにより、粒子ｊに対する下記の運動方程式が得られる。
Figure JPOXMLDOC01-appb-M000023
Figure JPOXMLDOC01-appb-M000024
　分子動力学では、粒子系を構成する各粒子に対し、式（７）及び式（８）で表される運動方程式を数値的に積分して解くことにより、各時刻における各粒子の運動量ベクトルｐｊ及び位置ベクトルｑｊが求められる。数値積分には、Verlet法が用いられることが多い。Verlet法については、例えば、J.M.Thijssen,「Computational Physics」,Cambridge University Press (1999)の１７５頁に説明されている。分子動力学計算で得られた各粒子の運動量ベクトル及び位置ベクトルに基づいて、粒子系の種々の物理量を算出することができる。
　次に、くりこみ群の手法を応用した分子動力学（以下、くりこみ群分子動力学という。）について概念的に説明する。
　くりこみ群分子動力学では、シミュレーション対象の粒子系Ｓを、粒子系Ｓの粒子数よりも少ない数の粒子で構成された粒子系Ｓ’（以下、くりこまれた粒子系Ｓ’という。）に対応づける。次に、くりこまれた粒子系Ｓ’に対して分子動力学計算を実行する。くりこまれた粒子系Ｓ’に対する計算結果を、シミュレーション対象の粒子系Ｓに対応付ける。これにより、シミュレーション対象の粒子系Ｓに対して分子動力学計算を直接的に実行する場合よりも計算量を少なくすることができる。シミュレーション対象の粒子系Ｓにおける物理量（例えば、粒子数、粒子の質量等）と、くりこまれた粒子系Ｓ’における物理量とを相互に対応付けるための変換則を、くりこみ変換測という。
　図１に、実施例によるシミュレーション方法のフローチャートを示す。ステップＳ１において、シミュレーション対象の粒子系Ｓに対して、くりこみ変換測に基づいてくりこみ変換処理を実行することにより、くりこまれた粒子系Ｓ’を定義する。粒子系Ｓは、複数のポリマーで構成され、複数のポリマーの各々は、一次元状に繋がった複数のモノマー粒子で構成される。
　粒子系Ｓのモノマー粒子間の相互作用ポテンシャルφ（ｒ）として、任意のモノマー粒子間においては、
Figure JPOXMLDOC01-appb-M000025
 が適用される。ここで、ｒは、粒子系Ｓを構成する各粒子からの距離を表す。
　ポテンシャルＵ０（ｒ）は、基本的に下記の式で表される。
Figure JPOXMLDOC01-appb-M000026
 ここで、
ｆは無次元化された関数を表し、ε及びσは、モノマー粒子を特徴づけるパラメータである。パラメータεはエネルギの次元を持ち、相互作用係数と呼ばれる。パラメータσは距離の次元を持ち、粒子の大きさに依存する。
　同一ポリマー内の隣接するモノマー粒子間においては、ポテンシャルＵ０（ｒ）に有限伸長非線形弾性ポテンシャルＵｃｈ（ｒ：ε，σ）が追加されて
Figure JPOXMLDOC01-appb-M000027
 が適用される。有限伸長非線形弾性ポテンシャルＵｃｈ（ｒ：ε，σ）は、ポテンシャルＵ０（ｒ）を定義するパラメータε及びσに依存するパラメータを含む。
　ポテンシャルＵ０（ｒ）として、例えばレナード・ジョーンズ型のポテンシャルを適用することができる。ポテンシャルＵ０（ｒ）に、シフトさせたレナード・ジョーンズ型のポテンシャルを適用した場合、ポテンシャルＵ０（ｒ）は、例えば下記の式で定義することができる。
Figure JPOXMLDOC01-appb-M000028
　ポテンシャルＵ０（ｒ）として、モース型ポテンシャルを適用することも可能である。さらに、ポテンシャルＵ０（ｒ）として、これらの関数をテーラー展開した近似関数を用いることも可能である。本実施例においては、ポテンシャルＵ０（ｒ）としてレナード・ジョーンズ型のポテンシャルを適用することとして説明を続ける。
　くりこみ変換処理により、シミュレーション対象の粒子系Ｓのモノマー粒子の位置ベクトルｑ、モノマー粒子の質量ｍ、ポリマー内のモノマー粒子の個数Ｎｍ、ポリマーの数Ｎｐ、ポテンシャルＵ０（ｒ）を定義するパラメータε、及びσが、それぞれ繰り込まれた粒子系Ｓ’のモノマー粒子の位置ベクトルｑ’、モノマー粒子の質量ｍ’、ポリマー内のモノマー粒子の個数Ｎｍ’、ポリマーの数Ｎｐ’、ポテンシャルＵ０’（ｒ）を定義するパラメータε’、及びσ’に変換される。ステップＳ１のくりこみ変換処理では、下記のくりこみ変換測が適用される。
Figure JPOXMLDOC01-appb-M000029
 ここで、ｄはシミュレーション対象の粒子系Ｓが配置された空間の次元数を表す。λは繰り込み回数に依存するくりこみ因子である。くりこみ回数がｎ回のとき、くりこみ因子λは以下の式で表される。
Figure JPOXMLDOC01-appb-M000030
　粒子系ＳにおけるポテンシャルＵ０（ｒ）及びＵｃｈ（ｒ：ε，σ）に対応する粒子系Ｓ’におけるポテンシャルＵ０’（ｒ）及びＵｃｈ’（ｒ：ε，σ）として、下記のポテンシャルを用いる。
Figure JPOXMLDOC01-appb-M000031
　粒子系Ｓ’の任意のモノマー粒子間に適用される相互作用ポテンシャルφ’（ｒ）は、下記の式で定義される。
Figure JPOXMLDOC01-appb-M000032
　同一ポリマー内の隣接するモノマー粒子間においては、ポテンシャルＵ０’（ｒ）に有限伸長非線形弾性ポテンシャルＵｃｈ’（ｒ：ε，σ）が追加されて、相互作用ポテンシャルφ’として、
Figure JPOXMLDOC01-appb-M000033
 が適用される。
　上述のくりこみ変換測を適用することにより、粒子の個数が１／λｄ倍になり、粒子の質量がλｄ倍になる。このため、粒子系Ｓの全体の質量と、粒子系Ｓ’の全体の質量とは同一である。さらに、粒子系Ｓ’におけるモノマー粒子の位置ベクトルｑ’が、粒子系Ｓにおけるモノマー粒子の位置ベクトルｑと同一であるため、粒子系Ｓとくりこまれた粒子系Ｓ’との巨視的な寸法は同一である。くりこみ変換処理の前後で、粒子系の全質量及び寸法が不変であるため、粒子系の密度も不変である。
　次に、ステップＳ２において、シミュレーションの初期条件を設定する。初期条件には、各粒子の位置ベクトルｑｊ及び運動量ベクトルｐｊの初期値が含まれる。運動量ベクトルｐｊは、くりこまれた粒子系Ｓ’の温度Ｔ’に基づいて設定される。シミュレーション対象の粒子系Ｓの温度をＴで表したとき、温度に関して下記のくりこみ変換測が適用される。
Figure JPOXMLDOC01-appb-M000034
　次に、ステップＳ３において、くりこまれた粒子系Ｓ’に関して分子動力学計算を実行する。運動方程式は、下記の式で表される。
Figure JPOXMLDOC01-appb-M000035
 　ここで、Ｎ’は粒子系Ｓ’のモノマー粒子の全個数を表しており、Ｎ’＝Ｎｍ’×Ｎｐ’である。
　この運動方程式を、数値的に積分して解く。これにより、くりこまれた粒子系Ｓ’の各粒子の位置ベクトルｑ’及び運動量ベクトルｐ’の時刻歴が求まる。
　くりこまれた粒子系Ｓ’の各粒子の位置ベクトルｑ’と、シミュレーション対象の粒子系Ｓの位置ベクトルｑとの対応関係は、式（１３）に示したとおりである。くりこまれた粒子系Ｓ’における運動量ベクトルｐ’と、シミュレーション対象の粒子系Ｓにおける運動量ベクトルｐとは、下記の関係を有する。
Figure JPOXMLDOC01-appb-M000036
　ステップＳ４において、シミュレーション結果を出力する。例えば、位置ベクトルｑ’及び運動量ベクトルｐ’を、そのまま数値で出力してもよいし、位置ベクトルｑ’に基づいて、粒子系Ｓ’の複数の粒子の空間内の分布を画像として表示してもよい。
　［φ（ｒ）＝Ｕ０（ｒ）の場合について考察］
　次に、相互作用ポテンシャルφ（ｒ）＝Ｕ０（ｒ）の場合について、くりこみ変換則を導き出す。まず、くりこまれた粒子系Ｓ’のハミルトニアンＨ’が
Figure JPOXMLDOC01-appb-M000037
 となることを示す。くりこみ前の粒子系Ｓの粒子数をＮで表し、くりこまれた粒子系Ｓ’の粒子数をＮ’で表す。このハミルトニアンＨ’を正準方程式に代入して運動方程式を得ることができる。ハミルトニアンＨ’のパラメータから、式（１３）のくりこみ変換則が導出される。
　くりこまれた粒子系Ｓ’のハミルトニアンＨ’を得るためには、粒子系Ｓに対する分配関数Ｚ（β）の積分の一部を実行し、ハミルトニアンを粗視化することによりハミルトニアンＨ’が得られる。
　粒子数が一定の正準集団（Canonical ensemble）に対する分配関数Ｚ（β）は、以下の式で表される。
Figure JPOXMLDOC01-appb-M000038
 ここで、ｋＢはボルツマン定数であり、ｄΓＮは位相空間内の体積要素である。ｄΓＮは下記の式で表される。
Figure JPOXMLDOC01-appb-M000039
 ここで、ｈはプランク定数である。ＷＮは、すべての状態の量子論的な和と、位相空間に亘る積分とが一致するように決められる。
　まず、粒子間の相互作用ポテンシャルの粗視化について説明し、次に運動エネルギの粗視化について説明する。続いて、相互作用ポテンシャルの粗視化、及び運動エネルギの粗視化を踏まえて、くりこみ変換測を定義する。
　［モノマー粒子間の相互作用ポテンシャルφ＝Ｕ０の粗視化］
　まず、一次元鎖状に配置された粒子系における相互作用ポテンシャルの粗視化について説明する。その後、単純立方格子状に配置された粒子における相互作用ポテンシャルについて説明する。
　図２に示すように、粒子ｉ、粒子ｊ、及び粒子ｋがこの順番に一次元状に並んで配置されている。粒子ｊが関与する相互作用を書き出し、粒子ｉと粒子ｋとの中間にある粒子ｊの位置座標について積分を実行することにより、相互作用ポテンシャルの粗視化を行うことができる。まず、第二近接以上の粒子からの寄与を取り込むために、ポテンシャルムービング（Potential Moving）の方法を用いることができる。ポテンシャルムービングの方法については、Leo P. Kadanoff, STATISTICAL PHYSICS Static, Dynamics and Renormalization, Chap. 14, World Scientific (1999)に解説されている。
　第二近接以上の粒子からの寄与を取り込んだ相互作用ポテンシャルφティルダは、下記の式で表すことができる。
Figure JPOXMLDOC01-appb-M000040
 ここで、ａは平衡状態における粒子間距離を表す。平衡状態における粒子間距離ａは、相互作用ポテンシャルφ（ｒ）が最小となる距離に等しい。
　複数の粒子が一次元状に配置されているため、粒子ｊの位置ベクトルｑｊは、一次元座標ｑｊで表すことができる。粒子ｊの位置をｑｊで表すと、粒子ｊに対して、最近接の粒子によって作られるケージポテンシャルは、下記の式で表される。
Figure JPOXMLDOC01-appb-M000041
　積分変数ｑｊについて積分を実行すれば、式（２５）を用いて下記の式が得られる。
Figure JPOXMLDOC01-appb-M000042
 ここで、ｒａは粒子の直径を表し、ｚ（ｑｉ－ｑｋ）及びＰ（ｑｉ－ｑｋ）は、下記の式で表される。
Figure JPOXMLDOC01-appb-M000043
Figure JPOXMLDOC01-appb-M000044
 積分領域はケージポテンシャルの内部領域に制限される。
　次に、ｚ（ｑｉ－ｑｋ）を具体的に計算する。相互作用ポテンシャルφがレナード・ジョーンズ型ポテンシャルの場合、φ（２ｎ）は、下記の式で表される。
Figure JPOXMLDOC01-appb-M000045
　式（２９）を式（２８）に代入して数値積分を行う。式（２８）に式（２９）を代入する際には、式（２４）を用いる。
　図３に、相互作用ポテンシャルφ（ｒ）がレナード・ジョーンズ型ポテンシャルである場合の式（２７）の算出結果、及び式（２８）の数値積分結果を示す。粒子ｉ、粒子ｊ、粒子ｋ、粒子ｌ、及び粒子ｍがこの順番に一次元状に配置されている場合の粒子ｋの位置座標がｑｋで表される。図３の横軸は、ｑｋ／２を表し、縦軸は、Ｐ（ｑｉ－ｑｋ）Ｐ（ｑｋ－ｑｍ）、及びｚ（ｑｉ－ｑｋ）ｚ（ｑｋ－ｑｍ）を、対数目盛で表す。図３に至る数値計算において、ε／ｋＢＴ＝２．０とした。
　図３から、Ｐ（ｑｉ－ｑｋ）Ｐ（ｑｋ－ｑｍ）の変化に比べて、ｚ（ｑｉ－ｑｋ）ｚ（ｑｋ－ｑｍ）の変化が緩やかであることがわかる。このため、Ｐ（ｑｉ－ｑｋ）Ｐ（ｑｋ－ｑｍ）に対してｚ（ｑｉ－ｑｋ）ｚ（ｑｋ－ｑｍ）は、ほぼ定数であると近似できる。
　粒子ｋが位置座標ｑｋに存在する確率ｐ（ｑｋ）は、以下のように近似することができる。
Figure JPOXMLDOC01-appb-M000046
　従って、下記の式が導出される。
Figure JPOXMLDOC01-appb-M000047
　ここまでは、一次元状に複数の粒子が配置された粒子系の相互作用ポテンシャルの粗視化について説明した。多次元の粒子系の相互作用ポテンシャルは、ポテンシャルムービングの方法により実現できる。
　図４Ａ～図４Ｃを参照して、二次元格子を一次元格子に帰着させるポテンシャルムービングの方法について説明する。
　図４Ａに示すように、二次元正方格子の格子点の位置に粒子が配置されている。最近接する粒子間の相互作用（ニアレストネイバーカップリング）を実線で示す。粒子が配列する１つの方向をｘ方向とし、それに直交する方向をｙ方向とする。
　図４Ｂに示すように、ｘ方向に並んだ粒子のうち、１つ置きに配置された粒子（中空の円で表された粒子）の変位について積分を実行することを考える。積分の対象となる粒子（消去したい粒子）を被積分粒子と呼ぶこととする。
　図４Ｃに示すように、ポテンシャルムービングの方法では、被積分粒子同士のニアレストネイバーカップリングを、ｘ方向に関して両隣に配置された粒子に振り分ける。振り分けられた相互作用が加算された粒子間相互作用（ダブルカップリング）を二重線で表す。二重線で表されたダブルカップリングは、元のニアレストネイバーカップリングの２倍の強さになる。この手法を用いて、二次元格子を一次元鎖に変換することができる。一次元鎖の粒子系においては、図２を参照して説明した手法により相互作用ポテンシャルの粗視化を行うことができる。シミュレーション対象の粒子系が三次元格子を構成する場合には、二次元格子を一次元鎖に変換する手順をｘ方向、ｙ方向、及びｚ方向の３方向に関して繰り返せばよい。このようにして、多次元格子を構成する粒子系の粗視化を行うことができる。
　次元数ｄの多次元格子を構成する粒子系の相互作用ポテンシャルの粗視化は、下記の式で表される。
Figure JPOXMLDOC01-appb-M000048
 ここで、＜ｉ，ｊ＞は、最近接格子間で和を取ることを意味する。
　最近接格子間の和を、すべての相互作用の和に変更すると、下記の式が得られる。
Figure JPOXMLDOC01-appb-M000049
　［運動エネルギの粗視化］
　次に、運動エネルギの粗視化について説明する。運動エネルギについては、容易に積分を実行することができ、下記の式が導出される。
Figure JPOXMLDOC01-appb-M000050
 式（３４）の導出に際し、下記の式が利用される。ここで、運動量ベクトルｐｊ ２は、ベクトルの内積を意味する。
Figure JPOXMLDOC01-appb-M000051
　［くりこみ変換測の導出］
　次に、上述の相互作用ポテンシャルの粗視化、及び運動エネルギの粗視化から導出されるくりこみ変換測について説明する。
　式（２２）に、式（３３）及び式（３４）を代入することにより、結果に影響しない係数を除いて下記の式が得られる。
Figure JPOXMLDOC01-appb-M000052
　式（３６）から、粗視化されたハミルトニアン（くりこまれた粒子系Ｓ’のハミルトニアン）Ｈ’は、以下の式で表される。
Figure JPOXMLDOC01-appb-M000053
　ハミルトニアンを粗視化する際の結合定数のリストをＫで表す。結合定数のリストＫは下記のように表される。
Figure JPOXMLDOC01-appb-M000054
　くりこみ変換Ｒを以下のように定義する。
Figure JPOXMLDOC01-appb-M000055
　ｎ回のくりこみ変換実行後の結合係数のリストＫｎは、下記の式で表される。
Figure JPOXMLDOC01-appb-M000056
　従って、ｎ回のくりこみ変換が行われたハミルトニアンＨｎは、下記の式で表される。
Figure JPOXMLDOC01-appb-M000057
 ここで、くりこまれた粒子系Ｓ’における運動量ベクトルｐｊ’は、下記の式で表される。
Figure JPOXMLDOC01-appb-M000058
　式（２１）と式（４１）とから、式（１３）に示したくりこみ変換測が導き出される。
　［モノマー粒子間の相互作用ポテンシャルφ＝Ｕ０＋Ｕｃｈの場合について考察］
　次に、ポリマー内の隣接するモノマー粒子間における相互作用ポテンシャルφ（ｒ）に関するくりこみ変換則の導出について説明する。
　図２において、粒子ｉ、ｋが作る粒子ｊに対するケージポテンシャルは、式（２５）と同様に、
Figure JPOXMLDOC01-appb-M000059
 と表すことができる。積分変数ｑｊについて、モノマー粒子ｊが一方のモノマー粒子ｉに接触した位置から、他方のモノマー粒子ｋに接触する位置まで積分する。モノマー粒子の直径をｒａで表すと、下記の式が得られる。
Figure JPOXMLDOC01-appb-M000060
　上述の式のように、くりこみ後の相互作用ポテンシャルは、２φ（（ｑｉ－ｑｋ）／２）とδφとの関数の重ね合わせとなる。
　有限伸長非線形弾性ポテンシャルＵｃｈ（ｒ：ε、σ）として、クレマーグレストモデルとして知られている以下のポテンシャルを採用することができる。
Figure JPOXMLDOC01-appb-M000061
 　ここで、パラメータｋ及びＲ０は、式（１０）に現れているパラメータε及びσに依存して決定される。一例として、パラメータｋ及びＲ０は下記の式で定義される。
Figure JPOXMLDOC01-appb-M000062
　式（４６）の関係を適用した場合、有限伸長非線形弾性ポテンシャルＵｃｈ（ｒ：ε、σ）も、式（１０）に示したポテンシャルＵ０（ｒ）と同様に、下記の式で表すことができる。ここで、ｆｃｈは、無次元化された関数である。
Figure JPOXMLDOC01-appb-M000063
　相互作用ポテンシャルφ（ｒ）は、
Figure JPOXMLDOC01-appb-M000064
 と表される。ポテンシャルＵ０（ｒ）は、上述の式（１２）で示したとおりである。式（４４）の２φ（（ｑｉ－ｑｋ）／２）を実際に計算し、δφを数値積分により求めた。
　図５に、２φ（（ｑｉ－ｑｋ）／２）の計算結果、及びδφの数値積分結果を示す。横軸は（ｑｉ－ｑｋ）／２を表し、縦軸はエネルギを単位「Ｊ」で表す。図５の実線は２φ（（ｑｉ－ｑｋ）／２）の計算結果を示し、破線はδφの数値積分結果を示す。図５の計算結果から、δφの変化量は、２φ（（ｑｉ－ｑｋ）／２）の変化量に比べて十分小さいことがわかる。従って、２φ（（ｑｉ－ｑｋ）／２）に対してδφは定数であると近似することができる。
　上述の考察から、下記の式が導出される。
Figure JPOXMLDOC01-appb-M000065
 　
　この式は、相互作用ポテンシャルφ＝Ｕ０の粗視化の説明で示した式（３１）と等価である。従って、相互作用ポテンシャルφ＝Ｕ０＋Ｕｃｈの粗視化も、相互作用ポテンシャルφ＝Ｕ０の粗視化と同様に取り扱うことができる。
　従って、くりこまれた粒子系Ｓ’における相互作用ポテンシャルφ’（ｒ）は、以下の式で表される。
Figure JPOXMLDOC01-appb-M000066
　上式に、式（１３）のくりこみ変換則を適用すると、次の式が得られる。
Figure JPOXMLDOC01-appb-M000067
　シミュレーションを行うときに、式（５１）を適用してもよい。この場合には、粒子系Ｓの相互作用ポテンシャルφ（ｒ）の変数ｒをｒ／λに置き換えてポテンシャルφ（ｒ／λ）を求め、その結果をλｄ倍することにより、粒子系Ｓ’について数値計算を行うことができる。
　［ポリマー数、ポリマー内のモノマー粒子の個数のくりこみ変換則］
　次に、図６を参照して、ポリマーの数Ｎｐ、及びポリマー内のモノマー粒子の個数Ｎｍのくりこみ変換則について説明する。
　図６に、直線状のポリマーが３次元方向に規則正しく配置されている状態の模式図を示す。ポリマーの各々を構成するモノマー粒子がｙ方向に繋がっている。複数のポリマーが、ｘ、ｙ、ｚ方向に配置されている。
　ｙ方向に関して１回のくりこみを行うと、ｙ方向に並ぶモノマー粒子が１つおきに消去される。例えば、矢印Ｐｙで示された面に配置されたモノマー粒子が消去される。これにより、ポリマー内のモノマー粒子の個数が１／２倍になる。次に、ｘ方向に関して１回のくりこみを行うと、ｘ方向に並ぶポリマーが１つおきに消去される。例えば、矢印Ｐｘで示された面に配置されたモノマー粒子が除去される。これにより、ポリマーの数が１／２倍になる。さらに、ｚ方向に関して１回のくりこみを行うと、ｚ方向に並ぶポリマーが１つおきに消去される。例えば、矢印Ｐｚで示された面に配置されたモノマー粒子が除去される。これにより、ポリマーの数がさらに１／２倍になる。
　上述のように、３次元方向に１回のくりこみを行うと、モノマー粒子の個数が１／２倍になり、ポリマーの数が１／４倍になる。くりこみ回数をｎ回とすると、式（１３）に示したＮｍ及びＮｐのくりこみ変換則が得られる。
　上記実施例によるシミュレーション方法、及びくりこみ変換を行わないシミュレーション方法を用いて、高分子に特有の現象である粘度のひずみ速度依存性を求めた。計算手法として、ＳＬＬＯＤ法を用いた。シミュレーション対象のパラメータは下記の通りである。
ｍ＝４２．３［ｇ／ｍｏｌ］
ε／ｋＢ＝４４３．０［Ｋ］
σ＝５．３×１０－１０［ｍ］
ｋ＝６．５３×１０－１［Ｎ／ｍ］
Ｒ０＝７．９５×１０－１０［ｍ］
Ｔ＝４４３．０［Ｋ］
Ｎｍ’＝６４
Ｎｐ’＝１０００
　ここで、ｋＢはボルツマン定数である。ε、σ、ｋ、及びＲ０は、式（１２）及び式（４５）に現れているパラメータである。パラメータｋ及びＲ０は、レナード・ジョーンズ型ポテンシャル（式（１２））のε及びσによって決定される。くりこみ後の粒子系Ｓ’におけるｋ’及びＲ０’は、くりこみ変換後のε’及びＲ０’により決定される。
　図７に、シミュレーション結果を示す。横軸はひずみ速度を表し、縦軸は粘度を表す。図７では、ひずみ速度及び粘度が無次元化されて表されている。図７中の丸記号はくりこみ変換を行わない粒子系Ｓについてシミュレーションを行った結果を示す、三角記号は、くりこみを行った粒子系Ｓ’についてシミュレーションを行った結果を示す。くりこみ後の粒子系Ｓ’について求めた粘度が、くりこみ前の粒子系Ｓについて求めた粘度を連ねる曲線に乗っていることがわかる。このように、くりこみを行った粒子系Ｓ’においても、粘度のひずみ速度依存性が保たれている。
　図７に示したシミュレーション結果により、上記実施例によるくりこみ変換則を用いたシミュレーション方法を用いて、高分子材料のマクロな計算を行うことが可能であることが示された。
　上記実施例によるシミュレーション方法は、コンピュータがコンピュータプログラムを実行することにより実現される。コンピュータプログラムは、例えばデータ記録媒体に記録された状態で提供される。または、コンピュータプログラムは、電気通信回線を介して提供することも可能である。
　図８に、実施例によるシミュレーション装置２５のブロック図を示す。シミュレーション装置２５は、入力装置２１、処理ユニット１０、及び出力装置２２を含む。処理ユニット１０は、入力部１１、演算部１２、出力部１３、及び記憶装置１４を含む。
　入力装置２１を通して入力部１１にシミュレーションの種々の条件が入力される。例えば、シミュレーション対象の粒子系Ｓを特徴づけるパラメータである空間の次元数ｄ、モノマー粒子の質量ｍ、１つのポリマーを構成しているモノマー粒子の数Ｎｍ、ポリマーの数Ｎｐの値、及び粒子系Ｓの初期条件が入力される、
　記憶装置１４に、シミュレーションを実行するコンピュータプログラムが格納されている。演算部１２は、入力部１１に入力されたデータに基づいて、記憶装置１４に格納されているコンピュータプログラムを実行することにより、シミュレーション処理を実行する。シミュレーション結果が、出力部１３を介して出力装置２２に出力される。
　演算部１２は、入力された粒子系Ｓのパラメータに対してくりこみ変換処理を行うことにより、くりこまれた粒子系Ｓ’のパラメータを算出する。算出されたパラメータを用いて、粒子系Ｓ’について分子動力学計算を行う。くりこまれた粒子系Ｓ’について分子動力学計算を行うため、計算時間を短縮することができる。
　以上実施例に沿って本発明を説明したが、本発明はこれらに制限されるものではない。例えば、種々の変更、改良、組み合わせ等が可能なことは当業者に自明であろう。
１０　処理ユニット
１１　入力部
１２　演算部
１３　出力部
１４　記憶装置
２１　入力装置
２２　出力装置
２５　シミュレーション装置

Claims (7)

  1. 　くりこみ回数に依存するくりこみ因子λに基づいて、各々が一次元状に繋がった複数のモノマー粒子からなる複数のポリマーで構成されたシミュレーション対象の粒子系Ｓに対してくりこみ変換処理を行う工程と、
    　くりこまれた粒子系Ｓ’について分子動力学計算を実行することにより、くりこまれた粒子系Ｓ’のモノマー粒子の位置ベクトル及び運動量ベクトルを求める工程と
    を有し、
    　粒子系Ｓ及び粒子系Ｓ’のモノマー間の距離をｒで表し、空間の次元数をｄで表し、
    　粒子系Ｓのモノマーの質量をｍで表し、１つのポリマーを構成しているモノマー粒子の数をＮｍで表し、ポリマーの数をＮｐで表し、モノマー粒子の座標をｑで表し、モノマー粒子間の相互作用ポテンシャルφ（ｒ）として、
    　任意のモノマー粒子間においては、
    Figure JPOXMLDOC01-appb-M000001
    が適用され、ポテンシャルＵ０（ｒ）は、エネルギの次元を持つパラメータεと、長さの次元を持つパラメータσと、無次元関数ｆを用いて
    Figure JPOXMLDOC01-appb-M000002
    と表すことができ、
    　同一ポリマー内の隣接するモノマー粒子間においては、パラメータとしてポテンシャルＵ０（ｒ）で用いられたパタメータε、σを持つ有限伸長非線形弾性ポテンシャルＵｃｈ（ｒ：ε，σ）がポテンシャルＵ０（ｒ）に追加されて
    Figure JPOXMLDOC01-appb-M000003
    が適用され、
    　粒子系Ｓ’のモノマー粒子の質量をｍ’で表し、１つのポリマーを構成しているモノマー粒子の数をＮｍ’で表し、ポリマーの数をＮｐ’で表し、モノマー粒子の座標をｑ’で表したとき、下記の変換則
    Figure JPOXMLDOC01-appb-M000004
    を適用し、
    　粒子系ＳにおけるポテンシャルＵ０（ｒ）、Ｕｃｈ（ｒ：ε，σ）に対応する粒子系Ｓ’におけるポテンシャルＵ０’（ｒ）、Ｕｃｈ’（ｒ：ε，σ）として、
    Figure JPOXMLDOC01-appb-M000005
    を用いて、粒子系Ｓ’について前記分子動力学計算を実行するシミュレーション方法。
  2. 　粒子系ＳにおけるポテンシャルＵ０（ｒ）は、レナード・ジョーンズ型ポテンシャルの形で表される請求項１に記載のシミュレーション方法。
  3. 　有限伸長非線形弾性ポテンシャルＵｃｈ（ε，σ，ｒ）は、パラメータε、σに依存するパラメータｋ及びＲ０を用いて、
    Figure JPOXMLDOC01-appb-M000006
    と表される請求項１または２に記載のシミュレーション方法。
  4. 　請求項１乃至３のいずれか１項に記載のシミュレーション方法をコンピュータで実行するためのコンピュータプログラム。
  5. 　請求項４に記載のコンピュータプログラムが、コンピュータで読み取り可能に記録された記録媒体。
  6. 　各々が一次元状に繋がった複数のモノマー粒子からなる複数のポリマーで構成されたシミュレーション対象の粒子系Ｓを特徴づけるパラメータである空間の次元数ｄ、モノマー粒子の質量ｍ、１つのポリマーを構成しているモノマー粒子の数Ｎｍ、ポリマーの数Ｎｐの値、及び粒子系Ｓの初期条件が入力される入力装置と、
    　前記入力装置から入力されたデータに基づいてシミュレーション処理を実行する処理ユニットと
    を有し、
    　前記処理ユニットは、
    　粒子系Ｓに対して、くりこみ回数に依存するくりこみ因子λに基づいてくりこみ変換を行う処理と、
    　くりこまれた粒子系Ｓ’について分子動力学計算を実行することにより、くりこまれた粒子系Ｓ’のモノマー粒子の位置及び運動量を求める処理と
    を実行し、
    　前記処理ユニットは、
    　粒子系Ｓ及び粒子系Ｓ’のモノマー間の距離をｒで表し、
    　粒子系Ｓのモノマー粒子の座標をｑで表し、モノマー粒子間の相互作用ポテンシャルφ（ｒ）として、
    　任意のモノマー粒子間においては、
    Figure JPOXMLDOC01-appb-M000007
    が適用され、ポテンシャルＵ０（ｒ）は、エネルギの次元を持つパラメータεと、長さの次元を持つパラメータσと、無次元関数ｆを用いて
    Figure JPOXMLDOC01-appb-M000008
    と表すことができ、
    　同一ポリマー内の隣接するモノマー粒子間においては、パラメータとしてポテンシャルＵ０（ｒ）で用いられたパタメータε、σを持つ有限伸長非線形弾性ポテンシャルＵｃｈ（ｒ：ε，σ）がポテンシャルＵ０（ｒ）に追加されて
    Figure JPOXMLDOC01-appb-M000009
    が適用され、
    　粒子系Ｓ’のモノマー粒子の質量をｍ’で表し、１つのポリマーを構成しているモノマー粒子の数をＮｍ’で表し、ポリマーの数をＮｐ’で表し、モノマー粒子の座標をｑ’で表したとき、下記の変換則
    Figure JPOXMLDOC01-appb-M000010
    を適用し、
    　粒子系ＳにおけるポテンシャルＵ０（ｒ）、Ｕｃｈ（ε，σ，ｒ）に対応する粒子系Ｓ’におけるポテンシャルＵ０’（ｒ）、Ｕｃｈ’（ε，σ，ｒ）として、
    Figure JPOXMLDOC01-appb-M000011
    を用いて、粒子系Ｓ’について前記分子動力学計算を実行するシミュレーション装置。
  7. 　くりこみ回数に依存するくりこみ因子λに基づいて、各々が一次元状に繋がった複数のモノマー粒子からなる複数のポリマーで構成されたシミュレーション対象の粒子系Ｓに対してくりこみ変換処理を行う工程と、
    　くりこまれた粒子系Ｓ’について分子動力学計算を実行することにより、くりこまれた粒子系Ｓ’のモノマー粒子の位置ベクトル及び運動量ベクトルを求める工程と
    を有し、
    　粒子系Ｓ及び粒子系Ｓ’のモノマー間の距離をｒで表し、空間の次元数をｄで表し、
    　粒子系Ｓのモノマーの質量をｍで表し、１つのポリマーを構成しているモノマー粒子の数をＮｍで表し、ポリマーの数をＮｐで表し、モノマー粒子の座標をｑで表し、モノマー粒子間の相互作用ポテンシャルφ（ｒ）として、
    　任意のモノマー粒子間においては、
    Figure JPOXMLDOC01-appb-M000012
    が適用され、
    　同一ポリマー内の隣接するモノマー粒子間においては、パラメータとしてポテンシャルＵ０（ｒ）で用いられたパタメータε、σを持つ有限伸長非線形弾性ポテンシャルＵｃｈ（ｒ：ε，σ）がポテンシャルＵ０（ｒ）に追加されて
    Figure JPOXMLDOC01-appb-M000013
    が適用され、
    　ポテンシャルφ（ｒ）は、エネルギの次元を持つパラメータεと、長さの次元を持つパラメータσと、無次元関数ｆを用いて
    Figure JPOXMLDOC01-appb-M000014
    と表すことができ、
    　粒子系Ｓ’のモノマー粒子の質量をｍ’で表し、１つのポリマーを構成しているモノマー粒子の数をＮｍ’で表し、ポリマーの数をＮｐ’で表し、モノマー粒子の座標をｑ’で表したとき、下記の変換則
    Figure JPOXMLDOC01-appb-M000015
    を適用し、
    　粒子系Ｓにおけるポテンシャルφ（ｒ）に対応する粒子系Ｓ’におけるポテンシャルφ’（ｒ）として、
    Figure JPOXMLDOC01-appb-M000016
    を用いて、粒子系Ｓ’について前記分子動力学計算を実行するシミュレーション方法。
PCT/JP2017/016774 2017-04-27 2017-04-27 シミュレーション方法、シミュレーションプログラム、及びシミュレーション装置 WO2018198273A1 (ja)

Priority Applications (2)

Application Number Priority Date Filing Date Title
PCT/JP2017/016774 WO2018198273A1 (ja) 2017-04-27 2017-04-27 シミュレーション方法、シミュレーションプログラム、及びシミュレーション装置
US16/450,341 US11164662B2 (en) 2017-04-27 2019-06-24 Simulation method, simulation program, and simulation device

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
PCT/JP2017/016774 WO2018198273A1 (ja) 2017-04-27 2017-04-27 シミュレーション方法、シミュレーションプログラム、及びシミュレーション装置

Related Child Applications (1)

Application Number Title Priority Date Filing Date
US16/450,341 Continuation US11164662B2 (en) 2017-04-27 2019-06-24 Simulation method, simulation program, and simulation device

Publications (1)

Publication Number Publication Date
WO2018198273A1 true WO2018198273A1 (ja) 2018-11-01

Family

ID=63919783

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/JP2017/016774 WO2018198273A1 (ja) 2017-04-27 2017-04-27 シミュレーション方法、シミュレーションプログラム、及びシミュレーション装置

Country Status (2)

Country Link
US (1) US11164662B2 (ja)
WO (1) WO2018198273A1 (ja)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP7122209B2 (ja) * 2018-10-01 2022-08-19 住友重機械工業株式会社 シミュレーション装置、シミュレーション方法、及びプログラム

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JPS524168A (en) 1975-05-21 1977-01-13 Du Pont Contacts for integrated circuits
JPS5241468B2 (ja) 1971-09-29 1977-10-18
JP2010146368A (ja) * 2008-12-19 2010-07-01 Sumitomo Heavy Ind Ltd シミュレーション方法及びプログラム
JP2016218767A (ja) * 2015-05-21 2016-12-22 住友重機械工業株式会社 シミュレーション方法、シミュレーションプログラム、及びシミュレーション装置
JP2017027217A (ja) * 2015-07-17 2017-02-02 住友重機械工業株式会社 粉粒体のシミュレーション方法及びシミュレーション装置

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JPS5241468B2 (ja) 1971-09-29 1977-10-18
JPS524168A (en) 1975-05-21 1977-01-13 Du Pont Contacts for integrated circuits
JP2010146368A (ja) * 2008-12-19 2010-07-01 Sumitomo Heavy Ind Ltd シミュレーション方法及びプログラム
JP2016218767A (ja) * 2015-05-21 2016-12-22 住友重機械工業株式会社 シミュレーション方法、シミュレーションプログラム、及びシミュレーション装置
JP2017027217A (ja) * 2015-07-17 2017-02-02 住友重機械工業株式会社 粉粒体のシミュレーション方法及びシミュレーション装置

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
LEO P. KADANOFF: "STATISTICAL PHYSICS Static, Dynamics and Renormalization", 1999, CAMBRIDGE UNIVERSITY PRESS, pages: 175

Also Published As

Publication number Publication date
US20190311786A1 (en) 2019-10-10
US11164662B2 (en) 2021-11-02

Similar Documents

Publication Publication Date Title
Dey et al. Metamodel based high-fidelity stochastic analysis of composite laminates: A concise review with critical comparative assessment
Zhang et al. Analysis of three-dimensional curved beams using isogeometric approach
JP6444260B2 (ja) シミュレーション方法、シミュレーションプログラム、及びシミュレーション装置
Whalen et al. Toward reusable surrogate models: Graph-based transfer learning on trusses
Haider et al. Parallel implementation of k-exact finite volume reconstruction on unstructured grids
Kumar et al. Synthesis of c0 path-generating contact-aided compliant mechanisms using the material mask overlay method
Hauser et al. State-space representations of deep neural networks
JP2020119189A (ja) 流体解析システム、流体解析方法、および流体解析プログラム
Ting Lin et al. Reliability-based multidisciplinary design optimization using probabilistic gradient-based transformation method
JP2018163396A (ja) 区分線形近似関数生成装置および方法
Macedo et al. What is shape? Characterizing particle morphology with genetic algorithms and deep generative models
JP6651254B2 (ja) シミュレーション方法、シミュレーションプログラム、及びシミュレーション装置
Liu et al. Cluster-based optimization of cellular materials and structures for crashworthiness
JP5782008B2 (ja) 非線形構造解析計算装置、非線形構造解析計算方法及び非線形構造解析計算プログラム
WO2018198273A1 (ja) シミュレーション方法、シミュレーションプログラム、及びシミュレーション装置
de Frias et al. A multiscale mass scaling approach for explicit time integration using proper orthogonal decomposition
JP6472346B2 (ja) 粉粒体のシミュレーション方法、プログラム、記録媒体、及びシミュレーション装置
Milewski et al. In search of optimal acceleration approach to iterative solution methods of simultaneous algebraic equations
JP5523364B2 (ja) 解析装置
JP2014013555A (ja) ニューラルネットワーク設計方法、フィッティング方法、及びプログラム
JP6067596B2 (ja) ペアリング演算装置、マルチペアリング演算装置、プログラム
Wei et al. Hybrid Polynomial Chaos Expansion and proper generalized decomposition approach for uncertainty quantification problems in the frame of elasticity
Xu et al. A new DROS-extreme learning machine with differential vector-KPCA approach for real-time fault recognition of nonlinear processes
Manteuffel et al. Projection multilevel methods for quasilinear elliptic partial differential equations: numerical results
US20230281267A1 (en) Calculation device, calculation program, recording medium, and calculation method

Legal Events

Date Code Title Description
121 Ep: the epo has been informed by wipo that ep was designated in this application

Ref document number: 17907295

Country of ref document: EP

Kind code of ref document: A1

122 Ep: pct application non-entry in european phase

Ref document number: 17907295

Country of ref document: EP

Kind code of ref document: A1

NENP Non-entry into the national phase

Ref country code: DE

NENP Non-entry into the national phase

Ref country code: JP