JP2017049805A - Method, device and program for calculating transportation coefficient - Google Patents

Abstract

PROBLEM TO BE SOLVED: To provide a method for calculating a transportation coefficient which reduces computation cost and improves calculation accuracy.SOLUTION: The method for calculating a transportation coefficient includes: a step ST3 for calculating a behavior of a molecule in a balanced state under the predetermined analysis temperature and analysis pressure on the basis of molecular kinetics computation using molecule model data and calculating (k) pieces of time-series data representing a physical quantity corresponding to the transportation coefficient from a time point tto a time point t; a step ST4 for dividing a period from the time point tto the time point tinto (m) (m<k) pieces of groups and calculating (m) pieces of values G(t) [i=1 to m] of auto-correlation functions on the basis of an average value of the physical quantities calculated for each group; steps (ST5 to ST8) for approximating the (m) pieces of values G(t) of auto-correlation functions with an approximation formula including a KWW (Kohlausch-Williams-Watts) function represented by exp{-(t/τ)}, and determining a parameter including τ and β; and a step ST9 for calculating the transportation coefficient by performing time integrating of the approximation formula using the determined parameter.SELECTED DRAWING: Figure 2

Description

  The present invention relates to a method, an apparatus, and a program for calculating transport coefficients such as a viscosity coefficient, a diffusion coefficient, and a heat conduction coefficient using molecular dynamics simulation.

  According to the fluctuation-dissipation theorem, the transport coefficients such as the diffusion coefficient, viscosity coefficient, and heat conduction coefficient of the material to be analyzed can be obtained by time-integrating the autocorrelation function of the physical quantity obtained from the molecular dynamics calculation. It has been known. For example, Patent Document 1 gives an example of obtaining a heat conduction coefficient as a transport coefficient. In Patent Document 1, molecular dynamics calculation is performed using model data in which a substance to be analyzed is reproduced at an atomic level, and the position of a molecule in an equilibrium state and its behavior are calculated. Based on the calculation result, time series data representing the time change of the heat flux is calculated, and the heat conduction coefficient is calculated from the autocorrelation function of the heat flux.

  Patent Document 2 proposes a means for removing this by paying attention to the fact that a contribution that does not actually exist is mixed into the molecular dynamics calculation in calculating the transport coefficient. The post-removal autocorrelation function is Laplace transformed into a fitting function, and the transport coefficient is calculated by time integration of the fitting function expressed by an exponential function (see paragraph 0036).

JP 2010-139500 A JP 2014-106554 A

The transport coefficient is given by an integral value of a time correlation function of a physical quantity. For example, assuming the shear stress Pxy, the viscosity coefficient η is expressed by the following formula (1). k _B is the Boltzmann constant, V is the volume, and T is the temperature. <...> indicates an autocorrelation function.

  As shown in the equation (1), in order to calculate the transport coefficient, it is necessary to integrate the value of the autocorrelation function from time 0 to infinity ∞. However, since infinite time calculation is impossible in reality, it is necessary to terminate the calculation in a finite time. As shown in FIG. 4, the autocorrelation function fluctuates in a long-time region, and this fluctuation affects the accuracy of calculating the value of the transport coefficient. As shown in FIG. 7, when the calculation of the transport coefficient proceeds, it converges to a substantially constant value at the portion indicated by a square in the figure, and then deteriorates due to the influence of fluctuation. Therefore, in order to reduce the influence of fluctuation of the autocorrelation function in the long time region, as shown in FIG. 7, the calculation of the value of the transport coefficient (in this example, the viscosity coefficient) is continued, and the value converges to a substantially constant value. Then, the calculation is continued until it can be confirmed that the condition is deteriorated, and the calculation is terminated. It is conceivable that the average value of the part that has converged to a certain value to some extent is adopted as the value used as the transport coefficient.

  However, in the above method, in order to reduce the influence of fluctuation of the autocorrelation function in the long time region, the integral calculation must be continued until it can be confirmed that the value of the transport coefficient converges to some extent and then deteriorates. In other words, a sufficiently long calculation time is required and the calculation cost is increased.

  Further, as in Patent Document 2, when the value of the autocorrelation function is Laplace transformed and fitting is performed using a fitting function expressed by an exponential function, a short time region is caused due to the shape of the exponential function. It was found that the value of can not be expressed with an exponential function with high accuracy, and as a result, the calculation accuracy of the transport coefficient is affected.

  The present invention has been made paying attention to such a problem, and an object thereof is to provide a method, an apparatus, and a program for calculating a transportation coefficient with reduced calculation cost and improved calculation accuracy. That is.

  In order to achieve the above object, the present invention takes the following measures.

That is, the method of calculating the transport coefficient of the present invention is as follows:
Based on molecular dynamics calculation using atomic model data representing the substance to be analyzed, the behavior of the molecule in an equilibrium state under a predetermined analysis temperature and pressure is calculated, and the physical quantity corresponding to the transport coefficient is calculated at the time t. calculating a k number of time-series data representing ₁ to time t _k,
Dividing the time point t ₁ to time point t _k into m groups (m <k), and calculating m autocorrelation function values based on an average value of physical quantities calculated for each group;
The values of the m autocorrelation functions are approximated by an approximate expression including a KWW function (Kohlausch-Williams-Watts) represented by exp {− (t / τ) ^β }, and parameters including τ and β are determined. And steps to
Calculating the transport coefficient by time-integrating the approximate expression using the determined parameter;
including.

The apparatus for calculating the transport coefficient of the present invention is:
Based on molecular dynamics calculation using atomic model data representing the substance to be analyzed, the behavior of the molecule in an equilibrium state under a predetermined analysis temperature and pressure is calculated, and the physical quantity corresponding to the transport coefficient is calculated at the time t. a physical quantity calculation unit for calculating the k pieces of time-series data representing ₁ to time t _k,
Dividing the time t _k from the time t ₁ to a group of m (m <k), and the autocorrelation function calculation unit for calculating a value of m of the autocorrelation function based on the average value of the calculated physical quantity for each group ,
The values of the m autocorrelation functions are approximated by an approximate expression including a KWW function (Kohlausch-Williams-Watts) represented by exp {− (t / τ) ^β }, and parameters including τ and β are determined. An approximation to
A transport coefficient calculation unit that calculates the transport coefficient by time-integrating the approximate expression using the determined parameter;
Is provided.

In this way, the time t ₁ to the time tk are _divided into m groups (m <k), and the value of the autocorrelation function is calculated based on the average value of the physical quantity calculated for each group. The fluctuation in the time domain can be suppressed, and it is not necessary to calculate more time than necessary to confirm the influence of the fluctuation, so that the calculation cost can be reduced. On the other hand, by taking the average value of physical quantities, the number of time-series data decreases, so the autocorrelation function on the linear scale is unknown, it becomes an autocorrelation function on the log scale, and an approximation with an exponential function is required. . Therefore, since the KWW function is used as an approximate expression, a short time region that cannot be expressed by a simple exponential function can be appropriately expressed, and the calculation accuracy can be improved. Therefore, the calculation cost can be reduced and the calculation accuracy can be improved.

The block diagram which shows typically the apparatus which calculates the transport coefficient of this invention. The flowchart which shows the method of calculating the transport coefficient of this invention. The figure which shows the time change of the relaxation elastic modulus computed by the Multiple-tau correlation method. The figure which shows the value of the autocorrelation function calculated from k time series data, without using a Multiple-tau correlation method. Explanatory drawing regarding the approximate expression which superimposed the exponential function. Explanatory drawing regarding the approximate expression which piled up the KWW function. Explanatory drawing regarding the method of calculating | requiring a transport coefficient by the conventional averaging.

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

[Device for calculating transport coefficient]
The apparatus of the present embodiment is an apparatus that calculates a transport coefficient of a substance to be analyzed under a predetermined analysis temperature and pressure using molecular dynamics simulation. Examples of the transport coefficient include a viscosity coefficient, a diffusion coefficient, and a thermal conductivity. In the present embodiment, an example of calculating a viscosity coefficient will be described.

  As shown in FIG. 1, the apparatus 1 includes an initial setting unit 10, a molecular dynamics calculation unit 11, a physical quantity calculation unit 12, an autocorrelation function calculation unit 13, an approximation unit 14, and a transport coefficient calculation unit 15. Have. These units 10 to 15 are realized by the cooperation of software and hardware by the CPU executing a processing routine (not shown) stored in advance in an information processing apparatus such as a personal computer having a CPU, memory, various interfaces, and the like. Is done.

  The initial setting unit 10 shown in FIG. 1 accepts an operation from a user via a known operation unit such as a keyboard and a mouse, and performs molecular dynamics calculation such as an atomic model representing an analysis target substance, an analysis temperature, and an analysis pressure. Various settings such as necessary settings are executed, and these setting values are stored in the memory.

The molecular dynamics calculation unit 11 shown in FIG. 1 uses various parameters such as the atomic model data, analysis temperature, and analysis pressure set by the initial setting unit 10, and each molecule in an equilibrium state based on the molecular dynamics calculation of the atomic model. Calculate time-series data related to the position of the model. Along with this, the physical quantity calculation unit 12 shown in FIG. 1 calculates the k pieces of time-series data representing a physical quantity corresponding to the transport coefficients on the basis of the calculation results of the molecular dynamics calculation part 11 from the time t ₁ to time t _k . k is the number at the time of reproducing the behavior of the atomic model, and is also the length of the analysis time. The length of the analysis time is preset by the initial setting unit 10. The calculation time becomes longer as k is larger. As shown in FIG. 1, the calculated k time-series data are stored in a memory. When the transport coefficient is a viscosity coefficient, it becomes a shear stress tensor, when the transport coefficient is a diffusion coefficient, it becomes a velocity, and when the transport coefficient is a heat conduction coefficient, it becomes a heat flux.

In this embodiment, according to the rocking dissipation theorem, the transport coefficient is given as an integral value of the time correlation function of the physical quantity. If the shear stress Pxy, the viscosity coefficient η is expressed by the following equation (1), if the velocity v, the diffusion coefficient D is expressed by the following equation (2), and if the heat flux density J, the heat conduction coefficient. λ is expressed by the following formula (3). Here, k _B is Boltzmann's constant, V is volume, and T is temperature. <...> indicates an autocorrelation function.

Autocorrelation function calculation unit 13 shown in FIG. 1, as shown in FIG. 1, divided time t _k from the time t ₁ to a group of m (m <k), calculates the average value of the physical quantity for each group Based on the average value, m autocorrelation function values G (t _i ) are calculated. However, i = 1 to m. In the present embodiment, m autocorrelation function values (for example, stress relaxation time series data) are calculated using the Multiple-tau correlation method. Of course, the method is not limited to the multiple-tau correlation method as long as the physical quantity is averaged by some method. For example, the maximum entropy method, the block average method, etc. are mentioned. In addition, although the average value is calculated from the time-series data (k pieces) of the physical quantity and the m autocorrelation function values are calculated in the memory of FIG. 1, this is schematically shown. It is explanatory drawing for doing, and is not exact.

1 approximates the autocorrelation function value G (t _i ) calculated by the autocorrelation function calculation unit 13 to a KWW function (Kohlausch-Williams) expressed by exp {− (t / τ) ^β }. Approximation with an approximate expression including -Watts) and determining parameters including τ and β.

In the present embodiment, the determination is made as follows.
The approximate expression is expressed by Expression (4) which is a superposition expansion expression of N KWW functions.

First, the approximation unit 14 determines parameters a _i , τ _i , and β _i (i = 1 to N) by the least square method with N = 1. Specifically, the parameters a _i , τ _i , β _i (i = 1 to N) so that the sum of square errors Δ ² with respect to the stress relaxation time series data (m autocorrelation function values) is minimized. To decide. The total sum Δ ² is expressed by equation (5).

Next, the approximation unit 14 determines that an error between the determined parameters a _i , τ _i , β _i and the approximate expression (4) and the values of the m autocorrelation functions is equal to or less than a predetermined allowable value. Until it becomes, N is increased by 1, and the determination of the parameters a _i , τ _i , and β _i by the least square method is repeated.
Specifically, the length N of the expansion formula starts from N = 1 and increases by 1 until the formula (6) is satisfied.

Equation (6) indicates that the degree of coincidence R ² only needs to be 99.9% or more. The error is 0.1% or less. A value 0.999 corresponding to the predetermined allowable value is set in advance via the initial setting unit 10. In the present embodiment, N = 3 and the condition of Expression (6) is satisfied.

The transport coefficient calculation unit 15 shown in FIG. 1 calculates the transport coefficient by time-integrating the approximate expression using the determined parameters. If it is a viscosity coefficient, it represents with Formula (7). Γ represents a gamma function. The parameters a _i , τ _i , β _i , and N are determined by the approximation unit 14.

  In order to show the effect of this embodiment, the viscosity coefficient of liquid styrene was calculated.

1) FIG. 3 is a diagram in which the change in relaxation elastic modulus with time is calculated by the method of the present embodiment (Multiple-tau correlation method). FIG. 4 is a diagram in which the value of the autocorrelation function is calculated from k pieces of time-series data without using the Multiple-tau correlation method. In FIG. 4, fluctuations that occur in the long time region exist, but in FIG. 3, it can be seen that the fluctuations are reduced. Incidentally, the relaxation modulus is because it is calculated by multiplying the V / k B _T of the value of the autocorrelation function, be multiplied to V / k B _T of the value of FIG. 4, FIG. 3 and FIG. 4 is a comparable .

  2) FIG. 5 is a diagram showing a value obtained by approximating the time change (stress relaxation) of the relaxation elastic modulus calculated by the method of the present embodiment (Multiple-tau correlation method) with an approximate expression in which exponential functions are superimposed. . The exponential function is exp {− (t / τ)}. The number of exponential superpositions was three. FIG. 6 is a diagram in which the time change (stress relaxation) of the relaxation modulus calculated by the method of the present embodiment (Multiple-tau correlation method) is approximated by an approximate expression in which the KWW function is superimposed. The number of superpositions of the KWW function was three. In FIG. 5, it can be seen that in the short-time region, the error between the approximate expression and the stress relaxation value is large, and the approximate expression is not appropriate. In FIG. 6, the error between the approximate expression and the stress relaxation value is smaller than in FIG. 5, and it can be said that the approximate expression is appropriate.

  3) FIG. 7 is a diagram in which the viscosity coefficient is calculated by time integration of Equation (1). As shown in the figure, the calculation is continued until it can be confirmed that the value of the viscosity coefficient has converged to a substantially constant value and then deteriorated, and the calculation is terminated. The calculation time is 5. Calculated to E + 04. For the viscosity coefficient, the average value of the portion (square in the figure) where the value is considered to have converged was adopted.

4) The calculated value and calculation time of the viscosity coefficient of liquid styrene are as follows.
Experimental value: 696 μPaS
7. Method using average value shown in FIG. 7: 707 μPaS; calculation time E + 04
Approximate expression using exponential function shown in FIG. E + 04
Approximate expression using KWW function shown in FIG. 6: 692 μPaS; E + 04
As described above, the conventional method of FIG. 7 can calculate the transport coefficient with a certain degree of accuracy with respect to the experimental value, but has a demerit that the calculation cost is high.
The method shown in FIG. 5 is preferable because the calculation cost is low, but the physical quantity in the short-time region cannot be expressed well by the approximate expression, so the calculation accuracy of the transport coefficient is poor.
The method shown in FIG. 6 has a low calculation cost, and the calculation accuracy of the transport coefficient is higher than other methods. This is considered to be because fluctuations in the long-time region are reduced by the effect of averaging by the Multiple-tau correlation method, and furthermore, the value in the short-time region can be appropriately expressed by an approximate expression.
Therefore, it can be understood that the method of the present embodiment can achieve both reduction in calculation cost and improvement in calculation accuracy. In this embodiment, the calculation error is reduced from 5% to 1%, and the calculation time is reduced to about 1/5.

[Method of calculating transport coefficient]
A method for calculating the transport coefficient using the apparatus 1 will be described with reference to FIG.

  First, in step ST1, the initial setting unit 10 shown in FIG. 1 executes various settings such as an atomic model representing a substance to be analyzed, settings necessary for molecular dynamics calculation such as analysis temperature and analysis pressure, and the like. Store the value in memory.

In the next step ST2, the molecular dynamics calculation unit 11 shown in FIG. 1 uses time-series data relating to the position of each molecular model in an equilibrium state based on the molecular dynamics calculation of the atomic model using various preset parameters. calculate. Along with this, the physical quantity calculation unit 12 shown in FIG. 1, at step ST3, the molecular dynamics calculation time t _k from the time t ₁ (shear stress tensor in the present embodiment) physical quantity corresponding to the basis transport coefficients on the calculation result of K time-series data expressed up to are calculated.

In the next step ST4, the autocorrelation function calculating section 13 shown in FIG. 1 divides the time t _k from the time point t ₁ to a group of m (m <k), based on the average value of the calculated physical quantity for each group The value of m autocorrelation functions is calculated. In the present embodiment, the Multiple-tau correlation method is used, but the present invention is not limited to this. For example, the maximum entropy method or the block average method can be used.

In steps ST5 to ST8, the approximating unit 14 shown in FIG. 1 converts the m autocorrelation function values G (t _i ) into a KWW function (Kohlausch-Williams) expressed by exp {− (t / τ) ^β }. Approximation with an approximate expression including -Watts) and determining parameters including τ and β.

Specifically, in step ST5, the approximation unit 14 sets N = 1. In the next step ST6, the approximating unit 14 approximates the m autocorrelation function values G (t _i ) by the equation (4) using the least square method, and parameters a _i , τ _i , β _i ( i = 1 to N) is determined. In the next step ST7, the approximating unit 14 determines that an error between the determined parameters a _i , τ _i , β _i and the value calculated by the approximate expression (4) and the values of the m autocorrelation functions is a predetermined tolerance. Determine if it is less than or equal to the value. In step ST7, N is increased by 1 (ST8) until it is determined that the error is equal to or smaller than a predetermined allowable value (ST7: NO), and step ST6 is repeated. Of course, instead of adopting the method of verifying by gradually increasing N in this way, an approximate expression in which a plurality of KWW functions are superimposed may be defined.

If it is determined in step ST7 that the error is equal to or less than a predetermined allowable value (ST7: YES), in the next step ST9, the transport coefficient calculation unit 15 uses the determined parameters a _i , τ _i , β _i . The approximate expression (4) is integrated using time to calculate the transport coefficient. Note that the expression obtained by integrating the approximate expression (4) with time is (7).

As described above, the method for calculating the transport coefficient of this embodiment is as follows.
Based on molecular dynamics calculation using atomic model data representing the substance to be analyzed, the behavior of the molecule in an equilibrium state under a predetermined analysis temperature and pressure is calculated, and the physical quantity corresponding to the transport coefficient is calculated at the time t. a step ST3 for calculating the k pieces of time-series data representing ₁ to time _{t k,}
The time point t ₁ to the time point tk are _divided into m groups (m <k), and m autocorrelation function values G (t _i ) [i = 1 to 1) based on the average value of physical quantities calculated for each group. m] for calculating ST],
The value G (t _i ) of m autocorrelation functions is approximated by an approximate expression including a KWW function (Kohlausch-Williams-Watts) represented by exp {− (t / τ) ^β }. Determining parameters to include (ST5-8);
Step ST9 for calculating the transport coefficient by time integrating the approximate expression using the determined parameters;
including.

The apparatus for calculating the transport coefficient of this embodiment is
Based on molecular dynamics calculation using atomic model data representing the substance to be analyzed, the behavior of the molecule in an equilibrium state under a predetermined analysis temperature and pressure is calculated, and the physical quantity corresponding to the transport coefficient is calculated at the time t. a physical quantity calculation unit 12 for calculating the k pieces of time-series data representing ₁ to time t _k,
The time points t ₁ to t _k are _divided into m groups (m <k), and m autocorrelation function values G (t _i ) [i = 1] based on the average value of physical quantities calculated for each group. ~ M], an autocorrelation function calculation unit 13;
The m autocorrelation function values G (t _i ) are approximated by an approximate expression including a KWW function (Kohlausch-Williams-Watts) represented by exp {− (t / τ) ^β }, and τ and β An approximation unit 14 for determining parameters including
A transport coefficient calculation unit 15 that calculates the transport coefficient by time-integrating the approximate expression using the determined parameters;
Is provided.

In this way, the time point t ₁ to the time point tk are _divided into m groups (m <k), and the autocorrelation function value G (t _i ) [i = 1] based on the average value of the physical quantities calculated for each group. ˜m] is calculated, it is possible to suppress fluctuations in the long-time region by averaging, and it is not necessary to calculate more time than necessary in order to confirm the influence of fluctuations, so that the calculation cost can be reduced. On the other hand, by taking the average value of physical quantities, the number of time-series data decreases, so the autocorrelation function on the linear scale is unknown, it becomes an autocorrelation function on the log scale, and an approximation with an exponential function is required. . Therefore, since the KWW function is used as an approximate expression, a short time region that cannot be expressed by a simple exponential function can be appropriately expressed, and the calculation accuracy can be improved. Therefore, the calculation cost can be reduced and the calculation accuracy can be improved.

In the present embodiment, the autocorrelation function value G (t _i ) is calculated using a multiple-tau correlation method, a maximum entropy method, or a block average method.

In the present embodiment, the approximate expression is expressed by Expression (4),
N = 1, parameters a _i , τ _i , β _i (i = 1 to N) are determined by the least square method,
Until the error between the determined parameters a _i , τ _i , β _i and the value G ′ (t) calculated by the approximate expression and the values of the m autocorrelation functions is equal to or less than a predetermined allowable value, N Is increased by 1 and the determination of the parameters a _i , τ _i , and β _i by the least square method is repeated.

  According to this method, the number of expansions of the approximate expression is increased so that the error is less than or equal to a predetermined allowable value. Therefore, the accuracy of the approximation satisfies the target allowable value, and as a result, the accuracy of the transport coefficient is further improved. It becomes possible.

  The program of the present embodiment is a program that causes a computer to execute each step constituting the above method. By executing this program, it is possible to obtain the operational effects of the above method. In other words, it can be said that the above calculation method is used.

  As mentioned above, although embodiment of this invention was described based on drawing, it should be thought that a specific structure is not limited to these embodiment. The scope of the present invention is shown not only by the above description of the embodiments but also by the scope of claims for patent, and further includes all modifications within the meaning and scope equivalent to the scope of claims for patent.

  For example, each of the units 10 to 15 illustrated in FIG. 1 is realized by executing a predetermined program by a CPU of a computer, but each unit may be configured by a dedicated memory or a dedicated circuit.

  The structure employed in each of the above embodiments can be employed in any other embodiment. The specific configuration of each unit is not limited to the above-described embodiment, and various modifications can be made without departing from the spirit of the present invention.

DESCRIPTION OF SYMBOLS 11 ... Molecular dynamics calculation part 12 ... Physical quantity calculation part 13 ... Autocorrelation function calculation part 14 ... Approximation part 15 ... Transport coefficient calculation part

Claims (7)

Based on molecular dynamics calculation using atomic model data representing the substance to be analyzed, the behavior of the molecule in an equilibrium state under a predetermined analysis temperature and pressure is calculated, and the physical quantity corresponding to the transport coefficient is calculated at the time t. calculating a k number of time-series data representing ₁ to time t _k,
Dividing the time point t ₁ to time point t _k into m groups (m <k), and calculating m autocorrelation function values based on an average value of physical quantities calculated for each group;
The values of the m autocorrelation functions are approximated by an approximate expression including a KWW function (Kohlausch-Williams-Watts) represented by exp {− (t / τ) ^β }, and parameters including τ and β are determined. And steps to
Calculating the transport coefficient by time-integrating the approximate expression using the determined parameter;
A method for calculating a transportation coefficient including:   The method according to claim 1, wherein the value of the autocorrelation function is calculated using a multiple-tau correlation method, a maximum entropy method, or a block average method. The approximate expression is expressed by Expression (4),

N = 1, parameters a _i , τ _i , β _i (i = 1 to N) are determined by the least square method,
Until the error between the determined parameters a _i , τ _i , β _i and the value G ′ (t) calculated by the approximate expression and the values of the m autocorrelation functions is equal to or less than a predetermined allowable value, N The method according to claim 1, wherein the parameter a _i , τ _i , and β _i are repeatedly determined by increasing the value by 1 and the least square method. Based on molecular dynamics calculation using atomic model data representing the substance to be analyzed, the behavior of the molecule in an equilibrium state under a predetermined analysis temperature and pressure is calculated, and the physical quantity corresponding to the transport coefficient is calculated at the time t. a physical quantity calculation unit for calculating the k pieces of time-series data representing ₁ to time t _k,
Dividing the time t _k from the time t ₁ to a group of m (m <k), and the autocorrelation function calculation unit for calculating a value of m of the autocorrelation function based on the average value of the calculated physical quantity for each group ,
The values of the m autocorrelation functions are approximated by an approximate expression including a KWW function (Kohlausch-Williams-Watts) represented by exp {− (t / τ) ^β }, and parameters including τ and β are determined. An approximation to
A transport coefficient calculation unit that calculates the transport coefficient by time-integrating the approximate expression using the determined parameter;
An apparatus for calculating a transport coefficient.   The apparatus according to claim 4, wherein the value of the autocorrelation function is calculated using a multiple-tau correlation method, a maximum entropy method, or a block average method. The approximate expression is expressed by Expression (4),

The approximation unit determines parameters a _i , τ _i , and β _i (i = 1 to N) by the least square method with N = 1,
Until the error between the determined parameters a _i , τ _i , β _i and the value G ′ (t) calculated by the approximate expression and the values of the m autocorrelation functions is equal to or less than a predetermined allowable value, N 6 is increased by 1, and the determination of the parameters a _i , τ _i , and β _i by the least square method is repeated.   The program which makes a computer perform the method in any one of Claims 1-3.