JP4733307B2 - Simulation device - Google Patents

Description

[0001]
BACKGROUND OF THE INVENTION
The present invention relates to a simulation apparatus useful for development of drugs and the like.
[0002]
[Prior art]
The Tsallis distribution [T] is a noncanonically extended generalized canonical distribution suitable for describing nonadditiveness, which is a feature of systems with long-range interactions and long-term memory effects. In recent years, it has attracted attention as a theoretical system of statistical mechanics. In the following description, a continuous Tsallis distribution in which variables x and p or their function E (x, p) takes a continuous value is an object. Such a continuous Tsallis distribution is defined by the following probability density (note that this definition is not unique and various modifications are made).
[0003]
[Expression 2]
[0004]
Where x = (x ₁ , ..., x _n ) Is a set of (generalized) coordinates, p = (p ₁ , ..., p _n ) Represents a set of momentums (conjugate to it). E (x, p) = U (x) + K (p) is total energy, U is potential energy, and K is kinetic energy.
[0005]
[Equation 3]
[0006]
It is. (The mass was all 1). , C is a normalization factor, β = 1 / kT, T is a parameter corresponding to temperature (hereinafter, temperature), and k is a Boltzmann constant. q is a real parameter, ρ depending on the range of the value of E _Tsallis Is given in the range defined as a real number. Boltzmann-Gyps (BG) density traditionally used to characterize the canonical distribution in the q → 1 limit
[0007]
[Expression 4]
[0008]
Is equal to This is one reason why the distribution defined in (1.1) is called a generalized canonical distribution.
Regarding the physical and chemical properties of a substance, the information usually obtained in the laboratory can be said to be the information of the system in the heat bath under a certain temperature. The physical system to which attention is focused as an application target of the following description is a system represented by a set of atoms or molecules of several to several million atoms. Of particular interest is a system represented by a (biological) polymer that is the target of a drug and molecules that constitute the environment in which it works. A typical example of the potential energy function in such a system is as follows:
[0009]
[Equation 5]
[0010]
The first term represents the potential energy determined by the distance between each atom, and the second term represents the sum of potential energies representing the equilibrium values such as the bond length and bond angle of atoms in the molecule.
[0011]
The BG distribution has been well studied as a distribution representing the state of the physical system in the heat bath. In particular, the dynamic properties of the polymer described in (1.3) and the system constituting its environment can be well explained by thermodynamic quantities such as free energy based on the BG distribution. It has been revealed by research. For example, the stable three-dimensional structure taken by a protein that is a biopolymer has the lowest free energy based on the BG distribution among various three-dimensional structures that can construct a polypeptide chain having a unique amino acid sequence that defines each protein. The state is considered selected. In addition, the affinity between a drug and a receptor protein is quantitatively explained by the difference in free energy between a system that exists as a complex in which both are bound and a system in which both are discrete.
[0012]
However, the system represented by the above (1.3) is generally complicated in terms of energy, and it has been difficult for the conventional technique to accurately obtain the BG distribution. Furthermore, the Tsallis distribution (1.1) is power-law, and the energy when these distributions are viewed as a function of the energy E, as is clear compared with the exponential BG distribution (1.2). In the case of a Tsallis distribution with q> 1, the dependence on E decreases more gradually than the BG distribution as the energy E increases, and is a set of states having a wide range of energy values. That is, when sampling is performed in the phase space based on the Tsallis distribution and a set of (x, p) pairs is created, the phase space much wider than the BG distribution is sampled in the Tsallis distribution with q> 1. Will be. Therefore, when calculating the thermodynamic quantity by numerical calculation, if a calculation is performed based on this sampling data, a wider phase space is sampled when the Tsallis distribution is used, so a better result is obtained. It is done.
In fact, there is an example where the concept of Tsallis distribution is used for effective sampling [AS]. If effective sampling can be realized by realizing the Tsallis distribution and the BG distribution can be accurately reconstructed from the Tsallis distribution, a stable three-dimensional structure of the protein can be predicted from the genome sequence, Applications to drug design on computers that estimate affinity are expected.
[0013]
Furthermore, the Tsallis distribution itself has begun to be considered as a candidate for a distribution representing the state of an actual physical system in the heat bath, like the conventional BG distribution [PP, AR1, AR2]. Therefore, by realizing the Tsallis distribution on a computer, it is possible to obtain information on the substance in the heat bath with a computer. Furthermore, for real physical phenomena that cannot be handled well by BG additive statistical mechanics, There is a report that it can be explained [A]. Therefore, by realizing the Tsallis distribution on a computer, it is expected to obtain completely new knowledge about the behavior of the physical system in the heat bath. The target system can be of various levels:
・ From microscopic [BCM] at the elementary particle level such as the transverse momentum distribution of hadrons and jets generated by electron and positron scattering,
・ Large scale systems such as the velocity distribution of galactic clusters [L et al]
An example in which these features can be described by the Tsallis distribution has been reported. Basically, the target is a system in which an energy function is defined, and the effect of the Tsallis distribution is expected especially in a system having a long-range interaction or a long-time memory effect. The first term of equation (1.3) is generally a long-range interaction. By realizing the Tsallis distribution in these systems, the dynamic behavior and structure formation of complex (biological system) polymers are completely eliminated. We can expect new knowledge and drug design from a new perspective.
[0014]
The following types of conventional methods for realizing the Tsallis distribution are known:
(1) q = 1 Tsallis distribution (ie, BG distribution) using both (x, p) variables
(2) Tsallis distribution of q ≧ 1 using x variable
(3) Tsallis distribution with q ≦ 1 using both (x, p) variables
Among the methods for realizing each, the following are typical deterministic methods.
[0015]
(Method 1) Nos é and Hoover propose an ODE (Ordinary Differential Equation) for a physical system with an energy function defined, and the BG distribution of temperature T can be obtained simply by adding one degree of freedom. Shown that it is feasible [N, H]. As a result, the simulation of the physical system in the heat bath was made realistic, and the knowledge about the canonical (Boltzmann-Gyps) group obtained by the simulation was dramatically increased. The Nos é -Hoover (NH) method is widely used in the field of molecular dynamics (MD). This equation is called the Nos é -Hoover equation and is expressed below.
[0016]
[Formula 6]
[0017]
Here, U is a potential function, and Q is a positive parameter.
(Method 2) Andricioaei and Straub proposed a method to realize the Tsallis distribution, applied to chemical systems, and demonstrated their effectiveness [AS]. In this approach aimed at effective sampling, a method based on Newton's equation was proposed, considering only the coordinate variable x as the Tsallis distribution.
[0018]
(Method 3) Plastino and Anteneodo's method [PA] is a method for generating a Tsallis distribution. However, unlike Method 2, it considers both (x, p) variables. This method supports mainly when 0 <q ≦ 1. (Velocity distributions of the galactic clusters have been reported as phenomena where analysis using a distribution of q <1 is effective.)
On the other hand, as a conventional method for performing sampling in the BG distribution, there is a method in which the BG distribution and a Tsallis distribution that is an extension thereof are realized. Of these, the above method 1-3 can be cited as a representative method as a deterministic method.
[0019]
[Problems to be solved by the invention]
As a conventional method for realizing the Tsallis distribution for the above (1)-(3), there is a method 1-3, but a method for realizing a Tsallis distribution of q ≧ 1 using both (x, p) variables is still successful. Not. The present invention achieves this distribution.
[0020]
Next, problems of the conventional method when performing sampling in the BG distribution will be described. (Method 1) This method is often used as a method for generating a BG distribution, but the setting of the parameter Q is not simple [N2, PHV]. Since Q → ∞ approaches a pure Newton equation with constant friction, if this is too large, the effect of temperature control on temperature T is reduced. Conversely, it is known that if Q is too small, the nonlinearity of the friction term becomes too large and the numerical integration algorithm fails. Further, depending on the system to be handled, the point (x, p) cannot be sufficiently sampled, and the thermodynamic quantity may not be calculated correctly.
[0021]
Various methods other than the NH method are known as methods for generating the BG distribution [N2]. However, in the method of realizing the BG distribution using the deterministic equation using both (x, p) variables, it is essentially difficult to exceed the high potential barrier with respect to the target temperature. That is, the barrier U (x _w ) Increases, the probability P (x _w , P) ∝exp [−E (x _w , P) / kT] ≦ exp [−U (x _w ) / KT] decreases exponentially, so the barrier point x _w To reach, a statistically significant amount of very rare events need to be sampled. Accordingly, in order to perform a considerable number of samplings corresponding to such extremely rare phenomenon, an essentially long calculation time is required.
[0022]
(Method 2) This is a deterministic method for realizing the Tsallis distribution considering only the x variable. This is basically a method of using a Newton equation using a potential U ′ modified for a Tsallis distribution having x as a variable with respect to a potential U representing a system. The effect of reducing the force acting on the barrier of U ′ compared to U can be expected, and in fact, the sampling effect is improved as compared with the conventional method [AS, PW1, PW2, HO]. In application to chemical systems, attempts have been made to use appropriate temperature controls [PW1, PW2]. These represent the Tsallis distribution with only x variables, ie
[0023]
[Expression 7]
[0024]
It is an attempt to realize in such a form. Note that U ′ is
[0025]
[Equation 8]
[0026]
As defined in.
However, in the basic equations used in those attempts, the law of conservation of energy holds, so that U ′ (x)> E (x (0), p (0)) against the barrier of U ′ Cannot be exceeded, and a trap to local energy minima can occur. Therefore, it is trapped inside such a barrier, and sufficient sampling cannot be achieved. This is a problem depending on the initial value (x (0), p (0)), and there is room for improvement. In application, there are reports that clearly indicate the method of temperature control and reports that do not. However, in any report, whether or not (3.1) is truly realized is not clarified except in a simple case without a barrier (one-dimensional harmonic oscillator).
[0027]
(Method 3) When 0 <q ≦ 1, which is mainly supported by this method, E is a finite value E ₀ Above ≡1 / (1-q) β,
[0028]
[Equation 9]
[0029]
Is defined. However, as described in Method 1, this is disadvantageous for sampling (note that [PA] is not intended for sampling). High energy realization probability is small or zero (U (x) ≧ E _c This is because it becomes difficult to cross the barrier of the potential U. this is,(
x, p) -means that sufficient sampling of the space cannot be performed and correct thermodynamic functions cannot be calculated. On the other hand, when q> 1, the problem of overflow in the numerical solution is pointed out. The details are not disclosed.
[0030]
The present invention aims to improve the sampling capability as compared with the conventional method, and this is an implementation of a method for realizing the Tsallis distribution of q ≧ 1 using both the (x, p) variables described above. This is achieved naturally. This achievement makes it possible to obtain an accurate value of the thermodynamic quantity in the BG distribution.
[0031]
Thereby, the efficiency and precision of the drug design by computer simulation can be achieved.
An object of the present invention is to provide a simulation apparatus useful for simulation in drug design and the like.
[0032]
[Means for Solving the Problems]
The simulation apparatus of the present invention is a simulation apparatus for accurately calculating a thermodynamic quantity of a physical system or a chemical system, and has a parameter setting means for setting parameters necessary for the simulation and energy compared to the Boltzmann-Gyps distribution. Based on the parameters, a generalized distribution setting unit that sets a generalized distribution that gradually approaches 0, a differential equation setting unit that sets a differential equation that reproduces the generalized distribution, and Physical quantity calculating means for integrating the differential equation, sampling along the obtained trajectory, and calculating a long-term average of the physical quantity based on the sampling.
[0033]
According to the present invention, when a thermodynamic quantity is calculated using a Boltzmann-Gyps distribution, sampling is performed using not a Boltzmann-Gyps distribution but a generalized distribution that is less asymptotic to zero. Therefore, even when the shape of the energy function has a minimum value, sampling can be made effective by overcoming the barrier surrounding the minimum value. Therefore, the possibility that sampling is trapped by the minimum value of the energy function is reduced, and a more accurate thermodynamic quantity can be calculated.
[0034]
DETAILED DESCRIPTION OF THE INVENTION
In the following description of the embodiment of the present invention, the Tsallis distribution is assumed. However, the present invention is not limited to the Tsallis distribution, and the energy E when the BG distribution is viewed as a function of the energy E. The present invention can be similarly applied even if the distribution has a function whose distribution changes as it approaches 0 later than the exponential function.
[0035]
First, in the embodiment of the present invention, a Tsallis distribution of q ≧ 1 using both (x, p) variables is realized by a deterministic algorithm, specifically, an ordinary differential equation. The long-term limit of the staying time rate in the minute region Δ (x, p) in the phase space of the solution of the ordinary differential equation is expressed by ρ of (1, 1). _Tsallis An ordinary differential equation is defined so as to be proportional to (x, p). Here, Δ (x, p) means a minute region including an arbitrary point (x, p). In the embodiment of the present invention, the following ordinary differential equation is used.
[0036]
[Expression 10]
[0037]
Here, the symbol D represents a differential operator. Ζ is a control variable for generating a Tsallis distribution appropriately, and ρ _z (Ζ) is an arbitrary function given by the user.
[0038]
As will be described later, when this ordinary differential equation is solved, the solution changes with time in appearance frequency proportional to the Tsallis distribution density, and the Tsallis distribution is realized (Equation (5.4)). In addition, the expected value of the physical quantity in the Tsallis distribution is also obtained by sampling a point in the phase space along the orbit of this ordinary differential equation and averaging this for a long time by sampling in the direction of time evolution according to the ordinary differential equation. Can be calculated. It is shown that the long-time average is not usually dependent on the initial value.
[0039]
Thus, to obtain the expected value of the physical quantity A in the BG distribution by sampling that reproduces the Tsallis distribution,
[0040]
[Expression 11]
[0041]
It can be calculated by obtaining the ratio of the two long-time averages as follows. Here, the long-time average is represented by a superscript line (hereinafter the same). Where ρ _Tsallis And ρ _BG Are (1.1) and
[0042]
[Expression 12]
[0043]
Is given by Here, β ′ is 1 / kT ′, but does not have to be the same value as β in equation (1.1). That is, the temperature T ′ of the BG distribution may be different from the temperature T of the Tsallis distribution. Although the Tsallis distribution is obtained, the thermodynamic quantity (long-term average of the physical quantity) is obtained from the BG distribution because the physical or chemical of the target system is calculated from the thermodynamic quantity calculated using the BG distribution. This is because a method for estimating the physical properties has been established.
[0044]
To show that equation (4.1) is extensible, we explain that this is a special case of an equation. The equation is R ^{2n + 1} It is defined for a density function ρ arbitrarily given in the region Γ of (2n + 1 dimensional space (corresponding to phase space)):
[0045]
[Formula 13]
[0046]
here,
[0047]
[Expression 14]
[0048]
The function was defined in. Outer 1 is
[0049]
[Outside 1]
[0050]
X of the function Θ _i Ingredient, p _i Represents a partial differential at the point ω≡ (x, p, ζ) of the component and ζ component. In equation (4.2), the Liuville equation,
[0051]
[Expression 15]
[0052]
Can be confirmed (X represents the right side of equation (4.2)). Assuming that ρ satisfies the appropriate mathematical condition, the long-term limit of the orbital residence time rate for the region A ⊂Γ is almost all with respect to the measure ρdω, according to the Lieuville theorem and Berkov's ergodic theorem [M]. Exists for the initial value of (the long-time average long-time limit along the trajectory starting from a point in the phase space converges to a finite value for most points in the phase space. To do). If the system is ergodic with respect to this measure,
[0053]
[Expression 16]
[0054]
It becomes. That is, the realization probability density of the state ω∈Γ is proportional to the given ρ. Equation (4.1) is expressed by Equations (4.2) and (4.3)
[0055]
[Expression 17]
[0056]
Is obtained by substituting Thus, the realization probability density of (x, p) is defined as ρ _Tsallis It is intended to be proportional to (x, p).
It should be noted that a density other than the expression (1.1) may be used to define the Tsallis distribution [TMP], which can also be handled by the present invention (ρ in the expression (4.5) _Tsallis The definition of (x, p) may be changed) and is not limited to the form of (1.1). By the way, usually it is not the formula (1.1),
[0057]
[Expression 18]
[0058]
Is often defined using. However, when the expected value of the physical quantity is obtained, it is often used in the form of equation (1.1) obtained by raising the above to the qth power, as proposed in [TMP]. Since this embodiment is intended to obtain the expected value of the physical quantity, (1.1) is adopted as its definition for convenience. E (x, p) can also be handled as a more general function (ρ in equation (4.5)). _Tsallis (This corresponds to the definition content of (x, p) being changed). In addition, X satisfying the equation (4.4) has a trivial arbitrary property (X can be multiplied by a constant). Therefore, using this, for example, when Θ in the equation (4.3) has a physical dimension, (Θ can be scaled by a constant ε with this same dimension,
[0059]
[Equation 19]
[0060]
Considering the equation (4.2) for, and multiplying the right side by ε, the equation can be expressed by Θ having dimensions). This is also true for equation (4.1).
Under mathematical additional conditions, for equation (4.1), the long-term limit of function f exists for almost all initial values for the measure ρdω, according to Berkov's Ergod Theorem [M]. here,
[0061]
[Expression 20]
[0062]
It is. If the system is ergodic with respect to this measure,
[0063]
[Expression 21]
[0064]
In particular, the long-term average value of the physical quantity A as a function of x and p is the expected value of A (weight ρ in the Tsallis distribution). _Tsallis The spatial average). That is,
[0065]
[Expression 22]
[0066]
Similarly, for the small region Δ (x, p), the following holds:
[0067]
[Expression 23]
[0068]
The left side is the long-term limit of the frequency of trajectory visits to Δ (x, p). Therefore, the realization probability density of the point (x, p) is ρ _Tsallis Is proportional to
In the computer simulation, the value on the left side of (5.3) is set to a sufficiently large time step N.
[0069]
[Expression 24]
[0070]
Evaluate with approximate values such as
Further, from the equation (5.3), the following equation is established:
[0071]
[Expression 25]
[0072]
Thereby, the expected value (right side) of the physical quantity A in the BG distribution can be calculated by averaging the two long sides on the left side.
1 to 20 are diagrams showing the effectiveness of the method of the present embodiment by comparison between numerical calculation and theoretical values.
[0073]
First, an example in which a Tsallis distribution of q ≧ 1 using both (x, p) variables is realized. Next, an example in which the sampling capability in the BG distribution is improved over the conventional method will be shown. And the example of verification of the sampling capability in BG distribution in different temperature is shown.
[0074]
As a first example of realizing the Tsallis distribution, the potential function U is converted to a one-dimensional harmonic oscillator.
[0075]
[Equation 26]
[0076]
These are shown in FIGS. The temperature T was 1.0. For simplicity, all quantities are represented dimensionlessly. The initial values of ODE are x (0) = 0.1, p (0) = 0.0, ζ (0) = 0.0, and numerical integration is performed by the fourth-order Runge-Kutta method. Calculated for the value of. Table 1 shows the parameter values used in the simulation.
[0077]
[Table 1]
[0078]
Δt represents a time step size, and N represents a simulation time step. ρ when q is 3.0 (case 1) _Tsallis FIG. 1 and FIG. 2 show simulation values and theoretical values (equation (1.1)) obtained from the histogram. These are in good agreement. Similarly, the results when the value of q is 2.0 (case 2) are shown in FIGS. It can be seen that the simulation values and the theoretical values agree well.
[0079]
A second example of realizing the Tsallis distribution is the one in which the potential function U is a double well oscillator type (FIG. 3) (shown in FIGS. 5 to 9):
[0080]
[Expression 27]
[0081]
The potential function is expressed by the above equation, and the diagram is shown in FIG.
For a temperature T = 1.0, the height of the barrier is 10.0. The results when q is 3.0 (case 4) and 2.0 (case 5) are shown in FIGS. As in the case of the harmonic oscillator, the simulation value and the theoretical value are in good agreement.
[0082]
In order to verify the sampling ability in the BG distribution, the result obtained using the equation (5.6) in the simulation of case 4 is shown and compared with the result obtained by the conventional method. First, the density ρ of the BG distribution with T = 1.0 _BG FIG. 10 and FIG. 11 show simulation values and theoretical values (equation (1.2)) obtained from the histogram. These are in good agreement. Next, a simulation result (equation (5.5)) of an expected value in a BG distribution of energy E (x, p) as a physical quantity is shown (FIG. 12). The theoretical value of the expected energy is a definite integral
[0083]
[Expression 28]
[0084]
Is estimated to be 1.0248. Therefore, it can be seen that the simulation result converges to the correct value. As a comparison, the result of calculating the thermodynamic function in the same manner as the conventional method of Method 1 is shown. The temperature T, initial value, numerical integration method, and time increment are all the same as in case 4. The value of parameter 1 / 2Q in equation (2.1) is 10 ^-2 10 ^-1 ... 10 ^Four 10 ^Five , And 2.0 × 10 ^Five I went 9 ways. In all cases, the temperature of the system
[0085]
[Expression 29]
[0086]
However, the coordinate value x of the solution was localized in the valley region on the left side of the potential including the initial coordinate 0.1, and no movement occurred on the right side. Some of the simulation results are shown in FIGS. Here, FIG. 13 shows the parameter 1 / 2Q = 10. ^-2 , Initial value x = 0.1, FIG. 14 has the same initial value, parameter 1 / 2Q = 10, and FIG. 15 has the same initial value, parameter 1 / 2Q = 10 ^Three 16, the initial values are the same, and the parameter 1 / 2Q = 10 ^Five 17, the initial values are the same, and the parameter 1 / 2Q = 2.0 × 10 ^Five 18, the initial value is x = 0.9, and the parameter 1 / 2Q = 2.0 × 10 ^Five It is said.
In Method 1, which is the conventional method, the BG distribution is not correctly generated as a result of localizing the sampling to one well of the double well. Therefore, the thermodynamic function value is not calculated correctly. Parameter 1 / 2Q = 2.0 × 10 ^Five FIG. 17 shows the simulation result of the expected energy value for. Sampling is inadequate and does not converge to the correct value. As described above, if the value of the parameter 1 / Q is further reduced, the Newton equation is approached, so that the ability to generate the BG distribution is reduced, and these results are not improved. Conversely, it is known that if the value of parameter 1 / Q is increased too much, the numerical integration algorithm fails. In fact, even in this calculation, 1 / 2Q = 10 ⁶ Then, the coordinate x caused a numerical overflow, and the numerical integration algorithm broke down. Further, even when the initial coordinate value is set to the right valley region x = 0.9, localization is similarly performed, and sampling is insufficient (1 / 2Q = 2.0 × 10 6). ^Five This case is shown in FIG.
[0087]
In order to verify the sampling capability in the BG distribution at different temperatures, the same conditions as in the above case 4 were used, but the calculation was performed according to the present embodiment for the BG distribution having a temperature 2.0 that is twice as high. Density ρ _BG 19 and 20 show simulation values and theoretical values obtained from the histogram. You can see these agreements.
[0088]
FIG. 21 is a block configuration diagram of the simulation apparatus according to the embodiment of the present invention.
The simulation apparatus 1 includes an input unit 10, a calculation unit 11, and an output unit 12. Input data required for calculation is input to the input unit 10. The input data input to the input unit 10 is sent to the calculation unit 11. The calculation unit 11 uses the parameter values included in the input data to reproduce the Tsallis distribution as shown in the description of the embodiment (4.1), or the energy more slowly than the BG distribution. Reproduce a distribution that is asymptotic to 0 (which is “slow” asymptotic behavior of a function that is slower than an exponential function as a function of energy E, for example, asymptotic behavior is governed by a finite power of E) Integrate the differential equation by a sequential numerical method such as Runge-Kutta method. Then, sampling is performed from a point on the orbit obtained as a result of integrating the differential equation, and a physical quantity value and an expected value are calculated using the sampled point.
[0089]
The calculated physical quantity value or expected value is passed to the output unit 12 as output data and presented to the user by a display device (not shown).
FIG. 22 is a diagram illustrating an example of input data.
[0090]
As shown in the figure, as an example of input data, the number of degrees of freedom n and the potential function No (in this case, a plurality of different potential functions are registered in the calculation unit and selected from these functions). Note that it is also possible to input the function form directly and let the simulation device interpret it)), simulation conditions, boundary conditions, Tsallis distribution parameters, ρ _z There are a function parameter, a BG distribution temperature list, a frequency calculation parameter, and the like.
[0091]
The simulation conditions include the number of simulation steps, a time step size, and an output designation parameter. Further, the Tsallis distribution parameter includes q, temperature T, or β. ρ _z As the function parameter, there are a function No in the same meaning as the potential function No, a setting coefficient, and the like. The BG distribution temperature list includes the number of temperatures of the BG distribution to be calculated and the value of the temperature of the BG distribution to be calculated. The frequency calculation parameter includes a discrete width when the distribution sampling is formed as a histogram, a domain parameter for controlling the domain, and the like.
[0092]
FIG. 23 is a diagram illustrating an example of output data.
Output data includes time changes in the values of the state variables x and p solved by integrating differential equations, time changes in other variable values (ζ, potential energy, etc.), and expected values of physical quantities (time average). Change, time change (and final value) of expected value (time average) of physical quantity in BG distribution at each temperature (energy, specific heat, etc.), frequency of potential energy value (realization probability), frequency of kinetic energy value (realization probability) , Tsallis distribution density (density, peripheral distribution), BG distribution density (density, peripheral distribution), and the like. Among the above data, the thermodynamic quantities are the expected value of the physical quantity and the expected value of the physical quantity in the BG distribution at each temperature. Other data can be used to evaluate the accuracy of the simulation.
[0093]
24 and 25 are flowcharts showing the overall processing flow of the embodiment of the present invention.
First, in step S1, data input is performed. Next, constant definition is performed in step S2. In step S3, preparation for output to a file for outputting the calculation result is performed. In step S4, initial settings of variables are defined. In this case, how to set is selected between automatic generation (step S6) and reading from an external file (step S7) in step S5. When initialization of variables is completed, preparation for frequency calculation is performed in step S8. In step S9, numerical integration is started.
[0094]
In step S10, t is set to t + Δt, and numerical integration is executed in step S11. In executing the numerical calculation, an equation (4.1 in the case of the embodiment of the present invention) of an equation for solving the results of the z-density value calculation unit (step S16), the potential function calculation unit (step S17), and the like. The right side is calculated (step S18), and numerical integration is performed using this result. In this case, it is possible to select or add an integration method. Further, the calculation result of the potential function calculation unit is passed to the energy calculation unit in step S15 and used for energy calculation. In step S12, a boundary condition for numerical integration is processed. In step S13, energy calculation is performed. In step S14, the density of the distribution is calculated. At this time, the distribution includes a Tsallis distribution and a BG distribution. The density calculation uses data from the calculation part of each distribution.
[0095]
In step S19, a numerical integration check is performed based on the data from the z-density value calculation unit. In step S20, a numerical error check is performed. If there is an error, the numerical integration is terminated. If there is no error, in step S21, frequency calculation is performed to determine the frequency of potential and momentum energy and the density of the Tsallis distribution. In step S22, the expected value of the physical quantity is calculated. Examples of physical quantities include kinetic energy, total energy, energy in BG distribution, specific heat in BG distribution, and the like.
[0096]
In step S23, the calculation result is output as output data. In step S24, it is determined whether or not a calculation end condition predetermined by the user is satisfied. If not satisfied, the process returns to step S10 to repeat the calculation. If satisfied, the numerical integration is terminated. Finally, after the completion of numerical integration, as a check of the output data, the theoretical value of Tsallis distribution density (peripheral distribution) and BG distribution density (peripheral distribution) is calculated when possible, and compared with the output data. Examine the accuracy of and verify that an appropriate distribution has been generated. Thereby, it is evaluated whether or not the calculated physical quantity gives a correct value.
[0097]
FIG. 26 is a flowchart showing the flow of processing of the input unit.
First, in step S30, data is input. In step S31, constants for numerical integration are defined, and in step S32, file output preparation is performed. In step S33, the initial value of the variable is defined. At this time, it is selected whether the initial value of the variable is automatically generated or read from an external file. Alternatively, the user may input manually. In step S34, preparation for density / frequency calculation is performed, and the process proceeds to the calculation unit.
[0098]
FIG. 27 is a flowchart showing the flow of processing of the calculation unit.
In step S40, numerical integration is started. In step S41, the time parameter t is advanced by t + Δt and Δt. In step S42, numerical integration is executed. At this time, in the z-density value calculation unit, a function selected and added and the result of calculation using an appropriate function are substituted into the right side of the equation together with the calculation result of the potential function calculation unit, and calculation is performed. To do. The result is used for numerical integration calculation in step S42. In numerical integration calculations, it is convenient to be able to select and add calculation methods. In step S43, the boundary condition is processed, and in step S44, the energy calculation using the data from the potential function calculation unit is performed in the energy calculation unit, and the energy value is calculated using the energy calculation. In step S45, the density ratio is calculated. In step S46, numerical integration is checked. In step S47, it is determined whether there is a numerical error. If there is an error, numerical integration is terminated. If there is no error, frequency calculation is performed in step S48, and physical quantity expected value calculation is performed in step S49. In step S50, data is output, and in step S51, the calculation end condition is determined. If the end condition is not satisfied, the process returns to step S41 to repeat the calculation, but the end condition is satisfied. Ends the numerical integration.
[0099]
FIG. 28 is a flowchart showing the flow of processing of the output unit.
In step S60, the output data calculated by the calculation unit is confirmed. That is, if possible, the theoretical value of the Tsallis distribution density and the peripheral distribution of the BG distribution density is calculated, and in step S61, the output data is visualized to determine whether the user has performed a correct simulation. Output data.
[0100]
FIG. 29 is a diagram for explaining the calculation contents of the potential function value calculation unit.
When x (t) is updated as an integration result of the differential equation, the following calculation is performed. That is, the potential function value, the potential function partial function (Force Function value is known to indicate the force acting on the system), the adjustment of the energy origin, and if necessary, the solution Add a potential function to constrain. Here, what kind of potential function is used can be conveniently selected from a plurality of candidates and can be added.
[0101]
FIG. 30 is a diagram illustrating a hardware environment of a computer when the simulation apparatus according to the embodiment of the present invention is realized by a program.
The CPU 21 executes the processing according to the embodiment of the present invention by executing the program stored in the ROM 22 and the RAM 23 connected by the bus 20. The program copied to the RAM 23 is stored in a storage device 27 such as a hard disk, or is stored in a portable recording medium 29 such as a floppy disk, CD-ROM, or DVD.
[0102]
The program stored in the portable recording medium 29 is read by the reading device 28 and copied to the RAM 23 via the bus 20 or stored in the storage device 27 and installed in the computer.
[0103]
The input / output device 30 includes a display, a keyboard, a mouse, a template, and the like, and is used for transmitting input from the user to the CPU 21 and presenting the calculation result of the CPU 21 to the user.
[0104]
The communication interface 24 can be connected to the information provider 26 via the network 25 to download the program so that the CPU 21 can execute it. Alternatively, the program that the information provider 26 has may be executed in a network environment.
<References>
-[A] Jun Abe, Mathematical Sciences, No.439, No.440, No.441; No.442 (2000)
・ [AR1] S.Abe and AKRajagopal, Nonuniqueness of canonical ensemble theory arising from microcanonical basis, Phys. Lett. A272,341 (2000)
・ [AR2] S. Abe and A, K, Rajagopal, Microcanonical foundation for systems with power-law distributions, J. Phys. A33, 8733 (2000)
・ [AS] I. Andricioaei and JEStraub, On Monte Carlo and molecular dynamics method inspired by Tsallis statistics: Methodology, optimization and applications to atomic clusters, J. Chem. Phys. 107, 9117 (1997)
・ [BCM] I.Bediaga, EMFCurado and J.Miranda, A nonextensive thermodynamical equilibrium approach in e + e_! Hadrons, Physica A 286,156 (2000)
・ [H] WGHoover, Phys. Rev. A 31 1695 (1985)
・ [HO] UHEHansmann and Y. Okamoto, Tackling the protein folding problem by ageneralized-ensemble approach with Tsallis statistics, in NonextensiveStatistical Mechanics and Thermodynamics, eds.SRA Salinas and C. Tsallis, Braz. J.Phys. 29, 187 (1999 )
・ [L et al] A. Lavagno, G. Kaniadakis, M. Rego-Monteiro, P. Quarati and C. Tsallis, Nonextensive thermostatistical approach of the peculiar velocity function of galaxy clusters, Astrophysical Letters and Communications 35, 449 (1989)
・ [M] R. Mane, "Ergodic Theory and Differentiable Dynamics," p.35, Springer-Verlag (1987)
・ [N] S. Nos &eacute;, J. Chem. Phys. 81, 511 (1984)
・ [N2] S. Nos &eacute;, Prog. Theor.Phys. Suppl. 103, 1 (1991)
・ [PA] ARPlastino and C. Anteneodo, A dynamical thermostatting approach to non-extensive canonical ensembles, Ann. Phys. 255, 250 (1997)
・ [PHV] HAPosh, WG Hoover, and FjVesely, Phys. Rev. A 33, 4253 (1986)
・ [PP] ARPlastino and A. Plastino, From Gibbs microcanonical ensemble to Tsallisgeneralized canonical distribution, Phys. Lett. A 193, 140 (1994)
・ [PW1] Y. Pakand S. Wang, Folding of a 16_residue helical peptide using molecular dynamicssimulation with Tsallis effective potential, J. Chem. Phys. 111,4356 (1999)
・ [PW2] Y. Pakand S. Wang, J. Phys. Chem. B 2000, 104, 354 (2000)
・ [T] C. Tsallis, Possble generalization of Boltzmann-Gibbs statistics, J. Stat. Phys. 52, 479 (1988)
・ [TMP] C. Tsallis, RSMendes and AR Plastino, The role of constraints within generalized nonextensivestatistics, Physica A 261,534 (1998)
[0105]
【The invention's effect】
The present invention has the following effects.
It is possible to generate a Tsallis distribution with q ≧ 1 using both (x, p) variables. This has not been successful with conventional methods.
-Generation of BG distribution at any temperature is possible.
・ The thermodynamic quantity can be calculated by generating the BG distribution.
・ Since the sampling capability is improved over the conventional method, the correct value can be calculated even when the conventional method cannot calculate the correct distribution or thermodynamic quantity value due to insufficient sampling.
[Brief description of the drawings]
FIG. 1 is a diagram (part 1) showing the effectiveness of a technique according to the present embodiment by comparison between numerical calculation and theoretical value.
FIG. 2 is a diagram (part 2) showing the effectiveness of the method of the present embodiment by comparison between numerical calculation and theoretical value;
FIG. 3 is a diagram (part 3) illustrating the effectiveness of the technique of the present embodiment by comparison between numerical calculation and theoretical value;
FIG. 4 is a diagram (part 4) illustrating the effectiveness of the method of the present embodiment by comparison between numerical calculation and theoretical value;
FIG. 5 is a diagram (No. 5) showing the effectiveness of the technique of the present embodiment by comparison between numerical calculation and theoretical value;
FIG. 6 is a diagram (No. 6) showing the effectiveness of the method of the present embodiment by comparison between numerical calculation and theoretical value;
FIG. 7 is a diagram (No. 7) showing the effectiveness of the method of the present embodiment by comparison between numerical calculation and theoretical value;
FIG. 8 is a diagram (No. 8) showing the effectiveness of the technique of the present embodiment by comparison between numerical calculation and theoretical value;
FIG. 9 is a diagram (No. 9) showing the effectiveness of the technique of the present embodiment by comparison between numerical calculation and theoretical value;
FIG. 10 is a diagram (No. 10) showing the effectiveness of the technique of the present embodiment by comparison between numerical calculation and theoretical value;
FIG. 11 is a diagram (No. 11) showing the effectiveness of the technique of the present embodiment by comparison between numerical calculation and theoretical value;
FIG. 12 is a diagram (No. 12) showing the effectiveness of the technique of the present embodiment by comparison between numerical calculation and theoretical value;
FIG. 13 is a diagram (No. 13) showing the effectiveness of the technique of the present embodiment by comparison between numerical calculation and theoretical value;
FIG. 14 is a diagram (No. 14) showing the effectiveness of the method of the present embodiment by comparison between numerical calculation and theoretical value;
FIG. 15 is a view (No. 15) showing the effectiveness of the technique of the present embodiment by comparison between numerical calculation and theoretical value;
FIG. 16 is a view (No. 16) showing the effectiveness of the technique of the present embodiment by comparison between numerical calculation and theoretical value;
FIG. 17 is a view (No. 17) showing the effectiveness of the method of the present embodiment by comparison between numerical calculation and theoretical value;
FIG. 18 is a diagram (No. 18) showing the effectiveness of the method of the present embodiment by comparison between numerical calculation and theoretical value;
FIG. 19 is a diagram (No. 19) showing the effectiveness of the technique of the present embodiment by comparison between numerical calculation and theoretical value;
FIG. 20 is a diagram (No. 20) showing the effectiveness of the technique of the present embodiment by comparison between numerical calculation and theoretical value;
FIG. 21 is a block configuration diagram of a simulation apparatus according to an embodiment of the present invention.
FIG. 22 is a diagram showing an example of input data.
FIG. 23 is a diagram illustrating an example of output data.
FIG. 24 is a flowchart (No. 1) showing the overall processing flow of the embodiment of the present invention;
FIG. 25 is a flowchart (No. 2) showing the overall processing flow of the embodiment of the present invention;
FIG. 26 is a flowchart illustrating a process flow of an input unit.
FIG. 27 is a flowchart illustrating a processing flow of a calculation unit.
FIG. 28 is a flowchart illustrating a processing flow of the output unit.
FIG. 29 is a diagram illustrating the calculation contents of a potential function value calculation unit.
FIG. 30 is a diagram illustrating a hardware environment of a computer when the simulation apparatus according to the embodiment of the present invention is realized by a program.
[Explanation of symbols]
10 Input section
11 Calculator
12 Output section

Claims (7)

On the computer,
From the input unit, the simulation parameters required for calculating the thermodynamic quantities of physical systems or chemical system from the user, and receives an input of data for setting the differential equation contains at least potential function, the inputted differential Based on the data for setting the equation, when the statistical distribution defined by statistical mechanics is viewed as a function of energy, the change in distribution with increasing energy reproduces a generalized distribution that gradually approaches the zero later than the exponential function. Processing to set the differential equation to
Based on the input parameters, the process of integrating the set differential equations, along the resulting trajectory, to sample, to calculate the long-time average of the physical quantity based on the sampling,
A simulation program characterized by executing   The simulation program according to claim 1, wherein the generalized distribution is a Tsallis distribution.   The simulation program according to claim 2, wherein the differential equation is a simultaneous ordinary differential equation that reproduces a Tsallis distribution. In the simultaneous differential equations, x _i is a generalized coordinate, p _i is a momentum conjugate to the generalized coordinate, U is a potential function, n is an integer, T is a temperature, k is a Boltzmann constant, β = 1 / kT, and E When an energy function, D is a differential operator, q is a constant of 1 or more, ζ is a control variable for appropriately creating a Tsallis distribution, and ρ _z is an arbitrary function given by the user,
The simulation program according to claim 3, wherein the simulation program is given by: A simulation device for calculating a thermodynamic quantity of a physical system or a chemical system,
Accepts user input of parameters necessary for simulation to calculate physical or chemical thermodynamic quantities and data for setting differential equations including at least potential functions , and sets the input differential equations based on data for, when the statistical distribution defined by statistical mechanics viewed as a function of energy, the change in the distribution of as the energy becomes larger, to reproduce the generalized distribution of asymptotic to 0 slower than exponential derivative An input means for setting an equation;
A physical quantity calculating means for integrating the set differential equation based on the input parameters, sampling along the obtained trajectory, and calculating a long-term average of the physical quantities based on the sampling;
A simulation apparatus comprising: In the simulation device,
The input unit accepts input of parameters necessary for simulation to calculate physical or chemical thermodynamic quantities from the user and data for setting a differential equation including at least a potential function. based on the data for setting the equations, when the statistical distribution defined by statistical mechanics viewed as a function of energy, the change in the distribution of as the energy becomes larger, generalized distribution you asymptotic to 0 slower than exponential A process of setting a differential equation that reproduces
Based on the input parameters, the process of integrating the set differential equations, along the resulting trajectory, to sample, to calculate the long-time average of the physical quantity based on the sampling,
Simulation method characterized by to run. On the computer,
The input unit accepts input of parameters necessary for simulation to calculate physical or chemical thermodynamic quantities from the user and data for setting a differential equation including at least a potential function. based on the data for setting the equations, when the statistical distribution defined by statistical mechanics viewed as a function of energy, the change in the distribution of as the energy becomes larger, generalized distribution you asymptotic to 0 slower than exponential A process of setting a differential equation that reproduces
Based on the input parameters, the process of integrating the set differential equations, along the resulting trajectory, to sample, to calculate the long-time average of the physical quantity based on the sampling,
A computer-readable recording medium on which a simulation program is recorded.