Reaction path search program, reaction path search system, and reaction path search method

Abstract

To provide a reaction path exploring program, reaction path exploring system and reaction path exploring method that are easy to suppress increase in calculation costs in accordance with increase in number of atoms in an exploration of chemical reaction paths.SOLUTION: A reaction path exploring program is configured to: as to each of a plurality of fragment pairs made of atoms in a structure targeted for an exploration of a reaction path, calculate at least any of a two-story differential coefficient b and three-story differential coefficient a of E(Q) indicative of a potential energy surface (PES) of the structure in Q=Q0 corresponding to a first equilibrium state serving as one of an equilibrium state where the E(Q) takes a minimum value; and, as to the fragment pair higher in an order of priority based on at least a size of any of the three-story differential coefficient a and the two-story differential coefficient b out of the plurality of fragment pairs, implement a reaction path exploring in transition to a second equilibrium serving as other equilibrium from the first equilibrium prior to other fragment pair lower in the order of priority than the fragment pair higher in the order of priority.SELECTED DRAWING: Figure 1

Description

[Technical field]

[0001] The present invention relates to a reaction path searching program, a reaction path searching system, and a reaction path searching method for searching chemical reaction paths. [Background technology]

[0002] One method for searching for chemical reaction paths for a given structure is the artificial force induced reaction method (AFIR method) described in Non-Patent Document 1. The AFIR method is a technique in which a force is applied to a fragment pair consisting of a combination of multiple atoms in a structure to induce deformation of the structure and cause a transition to a different structure. At this time, the corresponding reaction path is obtained from the coordinate change that the system undergoes during the transition process. By systematically applying this operation to various fragment pairs in the structure, reaction paths to various different structures are automatically searched for. [Prior art documents] [Non-patent literature]

[0003] [Non-Patent Document 1] Osamu Maeda, Yu Harabuchi, "Exploring paths of chemical transformations in molecular and periodic systems: An approach utilizing force", Wiley Interdisciplinary Reviews: Computational Molecular Science, 11, e1538 Summary of the Invention [Problem to be solved by the invention]

[0004] The number of fragment pairs to be considered in a system consisting of a certain structure increases in proportion to the square of the number of atoms N contained in the system. In the AFIR method, the search for reaction paths is performed for each fragment pair. Therefore, if the number of atoms in the structure to be searched increases, the number of combinations of fragment pairs to be calculated increases dramatically, and the calculation cost for searching reaction paths may become enormous.

[0005] An object of the present invention is to provide a reaction path searching program, a reaction path searching system, and a reaction path searching method that can easily suppress an increase in calculation cost according to an increase in the number of atoms. [Means for solving the problem]

[0006] The reaction path searching program according to the present invention is a program that causes one or more computers to function as a reaction path searching system that calculates a reaction path in a structure consisting of a plurality of atoms as a structural change represented by a positional relationship between the plurality of atoms, and the reaction path searching system calculates a reaction path in a structure consisting of a plurality of atoms as a structural change represented by a positional relationship between the plurality of atoms, the reaction path searching system calculating a reaction path in a structure consisting of a plurality of atoms as a structural change represented by a positional relationship between the plurality of atoms, the reaction path searching system calculating a reaction path in a structure consisting of a plurality of atoms as a structural change represented by a positional relationship between the plurality of atoms, the reaction path searching system calculating a reaction path in a structure consisting of a plurality of atoms as a structural change represented by a positional relationship between the plurality of atoms, the reaction path searching system calculating a reaction path in a structure consisting of a plurality of atoms as a structural change represented by a positional relationship between the plurality of atoms st Let R s and R t are the covalent bond radii of the s-th atom and the t-th atom, respectively, ρ is either +1 or -1, and p and α are constant values, the function F of Formula 1 calculated based on Formula 2 AFIRThe system includes a structural change calculation means for calculating the structural change in which (Q) takes a minimum value, a differential coefficient calculation means for calculating, for each of a plurality of fragment pairs having different atomic combinations, at least one of the second derivative b and the third derivative a of E(Q) at the position of the plurality of atoms corresponding to a first equilibrium state, which is one of the equilibrium states in which E(Q) takes a minimum value, and an inter-equilibrium state change calculation means for causing the structural change calculation means to calculate, for a fragment pair among the plurality of fragment pairs, a change in structure in the transition from the first equilibrium state to a second equilibrium state, which is another of the equilibrium states, for a fragment pair having a high priority based on the magnitude of at least one of a and b calculated by the differential coefficient calculation means, prior to other fragment pairs having a lower priority.

[0007] (Formula 1) TIFF0007655553000001.tif12170

[0008] (Formula 2) TIFF0007655553000002.tif12170

[0009] A reaction path searching system according to another aspect of the present invention is a system for calculating a reaction path in a structure consisting of a plurality of atoms as a structural change represented by a positional relationship between the plurality of atoms, and for a fragment pair consisting of fragment A consisting of N1 atoms (N1: natural number) extracted from the plurality of atoms and fragment B consisting of fragment A and another N2 atoms (N2: natural number) extracted from the plurality of atoms, a potential energy in a geometric parameter Q indicating the positions of the plurality of atoms is defined as E(Q), a distance between an sth atom (s: natural number satisfying s≦N1) belonging to fragment A and a tth atom (t: natural number satisfying t≦N2) belonging to fragment B is defined as r st Let R s and R tare the covalent bond radii of the s-th atom and the t-th atom, respectively, ρ is either +1 or -1, and p and α are constant values, the function F of Formula 1 calculated based on Formula 2 AFIR The system includes a structural change calculation means for calculating the structural change in which (Q) takes a minimum value, a differential coefficient calculation means for calculating, for each of a plurality of fragment pairs having different atomic combinations, at least one of the second derivative b and the third derivative a of E(Q) at the position of the plurality of atoms corresponding to a first equilibrium state, which is one of the equilibrium states in which E(Q) takes a minimum value, and an inter-equilibrium state change calculation means for causing the structural change calculation means to calculate, for a fragment pair among the plurality of fragment pairs, a change in structure in the transition from the first equilibrium state to a second equilibrium state, which is another of the equilibrium states, for a fragment pair having a high priority based on the magnitude of at least one of a and b calculated by the differential coefficient calculation means, prior to other fragment pairs having a lower priority.

[0010] A reaction path searching method according to another aspect of the present invention is a method for calculating a reaction path in a structure consisting of a plurality of atoms as a structural change represented by a positional relationship between the plurality of atoms, the method comprising the steps of: for a fragment pair consisting of fragment A consisting of N1 atoms (N1: a natural number) extracted from the plurality of atoms, and fragment B consisting of fragment A and another N2 atoms (N2: a natural number) extracted from the plurality of atoms, defining potential energy in a geometric parameter Q indicating positions of the plurality of atoms as E(Q), defining a distance between an sth atom (s: a natural number satisfying s≦N1) belonging to fragment A and a tth atom (t: a natural number satisfying t≦N2) belonging to fragment B as r st Let R s and R t are the covalent bond radii of the s-th atom and the t-th atom, respectively, ρ is either +1 or -1, and p and α are constant values, the function F of Formula 1 calculated based on Formula 2 AFIRThe method includes a structural change calculation step of calculating the structural change in which E(Q) takes a minimum value, a differential coefficient calculation step of calculating at least one of the second-order differential coefficient b and the third-order differential coefficient a of E(Q) at the position of the multiple atoms corresponding to a first equilibrium state, which is one of the equilibrium states in which E(Q) takes a minimum value, for each of the multiple fragment pairs having different atomic combinations, and an equilibrium state change calculation step of causing the structural change in the transition from the first equilibrium state to a second equilibrium state, which is another of the equilibrium states, to be calculated in the structural change calculation step for a fragment pair among the multiple fragment pairs having a high priority based on the magnitude of at least one of a and b calculated in the differential coefficient calculation step, prior to other fragment pairs having a lower priority.

[0011] According to the reaction path search program, reaction path search system, and reaction path search method of the present invention, a reaction path in a structure consisting of a plurality of atoms is calculated as a structural change represented by the positional relationship of the atoms. At this time, for fragment pairs with a high priority based on the magnitude of at least one of the second-order differential coefficient b and the third-order differential coefficient a of E(Q) at the position of the atoms corresponding to the first equilibrium state, the structural change is calculated prior to fragment pairs with a low priority. For example, a small second-order differential coefficient b indicates a tendency for a local increase in E(Q) to be small. Therefore, for example, when a fragment pair with a small b is prioritized, a path with a small increase in E(Q) near the first equilibrium state can be obtained. In addition, a small third-order differential coefficient a indicates a tendency for a change from a downward convex to an upward convex in E(Q). Therefore, as another example, when a fragment pair with a small a is prioritized, it is expected that a transition state is easily found near the first equilibrium state. Therefore, by performing calculations on some fragment pairs based on at least the magnitude of the second derivative b and the third derivative a of E(Q) over other fragment pairs, it is more likely that the desired reaction path can be obtained without incurring significant calculation costs. This makes it easier to suppress the increase in calculation costs associated with an increase in the number of atoms.

[0012] In the present invention, the reaction path search system further includes state selection means for selecting any one of N (N: a natural number of 2 or more) equilibrium states including all the second equilibrium states and the initial state obtained based on the past calculation results by the equilibrium state change calculation means as the first equilibrium state, and it is preferable to repeat the process of causing the equilibrium state change calculation means to calculate the change in the structure regarding the fragment pair in the transition from the first equilibrium state selected by the state selection means to the second equilibrium state by the structure change calculation means. According to this, any of the already acquired equilibrium states is used as the first equilibrium state, and the calculation of the change in the structure from the first equilibrium state to the second equilibrium state, that is, the acquisition of the reaction path is repeated. As a result, it is possible to acquire a network composed of the acquired equilibrium states including the initial state and the second equilibrium states acquired in the past, and the reaction paths connecting these equilibrium states.

[0013] In the present invention, the reaction path search system includes rate constant reduction means for obtaining, by repeating m times (m: a natural number satisfying m = N - l) the reduction of a rate constant matrix based on the RCMC method, a rate constant matrix composed of l * l rate constants related to the transition between l (l: a natural number satisfying l < N) super states represented by the weighted sum of the N equilibrium states from an N-row N-column rate constant matrix composed of N * N rate constants related to the transition between the N equilibrium states. For each of m = 1, 2,..., M (M: a natural number satisfying M < N), the state selection means sets each of the N equilibrium states as EQ i (i: a natural number not exceeding N) and each of the l super states as SS j (j: a natural number not exceeding l), sets the population of EQ i after m times of reduction as p i (m), sets the population of SS j after m times of reduction as Q j (m), sets the contribution of EQ j to SS i after m times of reduction as χ ji (m), and EQ iThe relative Gibbs energy of ΔG i The gas constant is R, and the model temperature parameter is T R When this is done, p is calculated based on Equation 3. i (m) and Λ calculated based on Equation 4 i The higher the EQ i It is preferable to select one of the N equilibrium states as the first equilibrium state so that the first equilibrium state is likely to be selected. i Equilibrium states with larger θ are more likely to be selected in calculating reaction paths.

[0014] (Formula 3) TIFF0007655553000003.tif20170

[0015] (Formula 4) TIFF0007655553000004.tif15170

[0016] In the present invention, the state selection means is i Large EQ i The more likely it is to be selected, the more likely it is to be selected. i The number of times that the structural change calculation means calculates the structural change with the first equilibrium state set as n i Small EQ i The number of times that the structural change calculation means has calculated the structural change for the N equilibrium states is n all The larger Λ i It is preferable to select one of the N equilibrium states as the first equilibrium state so that the difference in the number of equilibrium states is unlikely to be reflected in the difference in the likelihood of the equilibrium state being selected. all When is relatively small, Λ i is large and n i The smaller the equilibrium state, the more likely it is to be selected as a calculation target. This allows you to start calculating the structural change with a certain degree of n all Until becomes large, Λ i The equilibrium state where n is large (traffic volume is large) is likely to be selected appropriately. allAs becomes relatively large, Λ i The difference in the priority level is less likely to be reflected in the ease with which an equilibrium state is selected. Therefore, an equilibrium state with a small amount of traffic is also likely to be selected. In this way, with the above configuration, an equilibrium state is selected that is appropriately responsive to the priority level that changes with the progress of the calculation.

[0017] In the present invention, the state selection means defines the total number of times that the structural change calculation means has calculated the structural change for the N equilibrium states as n all Then, EQ i The number of times that the structural change calculation means calculates the structural change with the first equilibrium state being n i Then, ξ i is a real number between 0 and 1, and α and β are real numbers between 0 and 1, respectively. i The loudest EQ i It is preferable to select v as the first equilibrium state. i is Λ i The larger the larger, and n all When is relatively small, n i The smaller the value, the larger the value of n all When becomes relatively large, Λ i The difference is ν i Therefore, the difference in ν i The loudest EQ i By selecting the equilibrium state, an equilibrium state that appropriately corresponds to the priority that changes according to the progress of the calculation is selected.

[0018] (Formula 5) TIFF0007655553000005.tif11170 [Brief description of the drawings]

[0019] [Figure 1] FIG. 2 is a conceptual diagram of a fragment pair consisting of multiple atoms in the AFIR method used in accordance with one embodiment of the present invention. [Diagram 2]FIG. 1 is a conceptual diagram showing a network of reaction pathways obtained based on the AFIR method used in accordance with one embodiment of the present invention. [Diagram 3] FIG. 1 is a block diagram showing a functional configuration of a reaction path searching system according to one embodiment of the present invention. [Figure 4] 2 is a graph conceptually showing the change in PES around the equilibrium state EQ0 when the relative positions of the 0th layer atoms in the fragment pair of FIG. 1 are changed along a line passing through these atoms. [Diagram 5] FIG. 4 is a flow chart showing the flow of a series of processes executed by the reaction path searching system of FIG. DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS

[0020] A reaction path searching system 1 according to one embodiment of the present invention will be described with reference to the drawings. The reaction path searching system 1 searches for a reaction path in a structure such as a molecule composed of a plurality of atoms. The reaction path searching system 1 searches for a reaction path based on the artificial force induced reaction method (AFIR method), in particular, the single component AFIR method (SC-AFIR method). First, an overview of a reaction path searching method based on the SC-AFIR method will be described. Next, the configuration, functions, and other details of the reaction path searching system 1 will be described.

[0021] [Outline of reaction path search method] The AFIR method uses the AFIR function F shown in the following Equations 1 and 2. AFIR In this method, a structural transformation is induced by applying a force between a fragment pair consisting of fragments A and B as shown in FIG.

[0022] (Formula 1) TIFF0007655553000006.tif12170(Formula 2) TIFF0007655553000007.tif12170

[0023] Here, fragment A is a fragment of a structure consisting of N1 atoms (N1: natural number) in the above structure. Fragment B is a fragment of a structure consisting of another N2 atoms (N2: natural number) in the above structure. E(Q) corresponds to the potential energy surface (PES) in the geometric parameter Q indicating the position of the atoms in the above structure. r st is the distance between the sth atom (s: a natural number satisfying s≦N1) belonging to fragment A and the tth atom (t: a natural number satisfying t≦N2) belonging to fragment B. R s and R t are the covalent bond radii of the sth atom and the tth atom, respectively. ρ is set to either +1 or -1. p is set appropriately within the range of 4 to 8 (for example, p=6). α, which indicates the strength of the force acting between fragments, is a constant based on the following formula 1-1. R0 and ε in formula 1-1 are the Ar-Ar Lennard-Jones parameters, that is, 3.8164*10 -10 m and 1.0061 kJ / mol are used. α in Equation 1-1 corresponds to the average force that an Ar-Ar pair experiences when moving from the minimum point to the turning point when they collide directly with each other at the collision energy γ. Therefore, γ is called the model collision energy parameter, and corresponds to the approximate upper limit of the barrier that the force in the AFIR method can overcome. γ can be set appropriately by the user.

[0024] (Formula 1-1) TIFF0007655553000008.tif22170

[0025] The reaction path is obtained by minimizing the AFIR function by a standard quasi-Newton method (see Non-Patent Document 1). This path is called the "AFIR path". Along the AFIR path, approximate equilibrium (EQ) and transition state (TS) structures are obtained as structures where the AFIR function has a local minimum and maximum, respectively. These structures can be used as initial estimates of the actual EQ and TS. However, these energy maxima may deviate significantly from the actual TS. This situation may occur especially when γ is too large. Even in such cases, overlooking these TSs is avoided by performing a relaxation calculation of the AFIR path. For this purpose, any double-ended method that optimizes the entire path points can be used. As a method for this, the local plane update (LUP) method is adopted (J. Chem. Phys. 1991, 94, 751). According to one implementation using the LUP method, the path points are evenly distributed along a certain path and moved to low-energy points in a hyperplane perpendicular to the tangent of the path. The path obtained in this way is called the "LUP path". Here, the two endpoints and all local maxima are directly optimized to their minimum and maximum, respectively. The approximate EQs and TSs on the AFIR or LUP paths are optimized to the real EQs and TSs by standard quasi-Newton methods. Then, starting from all the obtained TSs, a calculation of the intrinsic reaction coordinates (IRCs) is performed. Finally, a normal mode analysis is performed on all the obtained EQs and TSs.

[0026] In the SC-AFIR method, AFIR paths are calculated from EQs (J.Comput.Chem.2014,35,166). The structures corresponding to EQs are molecules, molecular complexes, organometallic complexes, clusters, adsorption structures on surfaces, crystal structures, etc. From a single EQ, various AFIR paths are obtained by minimizing the AFIR function of fragment pairs consisting of various fragments A and B. In the structures corresponding to each EQ, fragment pairs are systematically defined by focusing on atom pairs. That is, as shown in Figure 1, in fragment A and fragment B, two atoms are assigned as 0th layer atoms. Then, the 1st layer atoms connected to the 0th layer atoms and the 2nd layer atoms connected to the 1st layer atoms are added to fragments A and B. In the 1st and 2nd layer atoms in each fragment, atoms whose distance to the atoms in the other fragment is less than the distance between the 0th layer atoms are removed from each fragment. A set of fragment pairs is generated by applying this procedure to all atom pairs except those with very long distances. Then, the AFIR path of each fragment pair is calculated. If AFIR paths are calculated for all fragment pairs, N atoms If is the number of atoms, the cost required to complete a calculation around a single EQ is N atoms 2 The reaction path network shown in Figure 2 is obtained by repeatedly calculating the AFIR path from the already obtained EQ and obtaining new EQs and paths between EQs. In the following, the reaction path network consisting of N EQs (i: a natural number of 2 or more) is called EQ i In this case, EQ1 represents the reactant as the initial structure. Figure 2 shows an example of a reaction path network for N=6.

[0027] Next, we explain the kinetic-based navigation adopted to reduce the computational cost of the reaction path. The number of EQs is the number of atoms, N atoms It has been shown that the function grows exponentially with N atomsUnless is small, it is difficult to apply the above calculation procedure based on the SC-AFIR method to all EQs. On the other hand, the number of EQs that can be dynamically accessed under mild experimental conditions is limited. Therefore, we use the rate constant matrix reduction (RCMC) method, which is a kinetic method, to determine whether an EQ is kinetically accessible under specific experimental conditions, thereby eliminating the search cost for kinetically inaccessible EQs (Chem. Lett. 2019, 48, 47).

[0028] The RCMC method enables kinetic analysis that takes into account the entire reaction pathway network consisting of N EQs by repeating an operation called contraction. An N-by-N rate constant matrix is ​​obtained from the N*N rate constants for transitions between N EQs. In the RCMC method, a smaller n-by-n rate constant matrix (n=NM) is obtained by M contractions. The n-by-n rate constant matrix after M contractions represents the rate constants for transitions between n superstates (SSs). The SSs are expressed as the weighted sum of all EQs. Each contraction reduces the size of the rate constant matrix by one, and all off-diagonal elements of the rate constant matrix are updated based on the quasi-steady state approximation and the Boltzmann distribution preservation condition. The off-diagonal elements of the rate constant matrix obtained by contraction correspond to the rate constants for transitions from one SS to another SS. Contraction is performed using a threshold t that can be set by the user. MAX It applies to all states with a lifetime shorter than t MAX A number of n SSs are obtained that transition between each other on a timescale longer than . A more detailed explanation of the RCMC method is given below.

[0029] EQ after m-fold contraction i The population of p i (m). A population is a set of EQs at a given time. i In kinetic-based navigation, the user needs to provide the initial population considering the experimental conditions. For example, if EQ1 is a reactant, the initial population of EQ1 is set to 1 (p1(0)=1), and the populations of all other EQs are set to 0 (p i≠1(0)=0). Then, all EQ i Based on Equations 3 and 4, the index Λ i , the so-called traffic volume, is calculated.

[0030] (Formula 3) TIFF0007655553000009.tif20170

[0031] (Formula 4) TIFF0007655553000010.tif15170

[0032] In formulas 3 and 4, Q (m) j Let each of the n SSs be SS j (j: a natural number satisfying j≦n), the SS after m contractions j The population of χ (m) ji is the SS after m contractions. j EQ to i is the contribution of ΔG i is the EQ obtained by normal mode analysis based on quantum chemical calculations. i is the relative Gibbs energy of the gas constant, T R is the model temperature parameter. Equation 4 is MAX EQ generated within i Population influx and EQ i For example, in a three-state system consisting of a reactant EQ, an intermediate EQ that is sufficiently more stable than the reactant EQ, and a product EQ that is sufficiently more stable than the intermediate EQ, if the experimental condition is that the path from the reactant EQ to the product EQ is open only via the intermediate EQ, then the reactant EQ and the product EQ are expressed as Λ i = 1.0, and the intermediate EQ is Λ i On the other hand, from the experimental conditions, EQ i If you do not have access to Λ i In this embodiment, as described later, such Λ iOne of the N EQs is selected using the value. Then, the AFIR path is calculated for the selected EQ. The AFIR path is calculated for a fragment pair appropriately selected from the fragment pairs that have not yet been processed. T R is a nonzero Λ i Determines how often a high-energy EQ is selected from among EQs with R =3000K or some other high value.

[0033] [System Details] The reaction path searching system 1 includes hardware including one or more computers, and software including a program that causes the hardware to function so as to realize various functions of the reaction path searching system 1. Each computer includes a central processing unit (CPU), a memory device such as a read-only memory (ROM) and a random access memory (RAM), and various interfaces such as an input / output interface. In addition, an input device that accepts user input, an output device such as a display, a recording device such as a hard disk, and the like may be connected to the computer as external devices or may be installed as internal devices. The software includes program data and the like recorded in a memory device or a recording device such as a hard disk. Then, the hardware such as a computer executes various information processing such as arithmetic processing and input / output processing according to the software. Various functions of the reaction path searching system 1 described below are realized by the functions of the hardware according to the software.

[0034] In addition, when the hardware includes multiple computers, the multiple computers may cooperate to execute the processes related to the various functions. There are various modes of such cooperation. For example, the computers may be connected to each other through a communication network such as the Internet, and each computer may share and execute the necessary arithmetic processing while communicating the necessary information with each other. In addition, one computer may be responsible for the user interface function, and the other computer may execute the arithmetic processing based on input information input by a user to the first computer and transmitted through the communication network.

[0035] As shown in Fig. 3, the reaction path searching system 1 includes a storage unit 10, a fragment pair derivation unit 20, a fragment pair ranking acquisition unit 30, a state selection unit 40, and a reaction path calculation unit 50. The storage unit 10 stores data indicating a list of acquired EQs, structure data indicating atomic coordinates, fragment data indicating combinations of atoms in fragments A and B, search count data, and the like. The list of acquired EQs includes EQ1 corresponding to a reactant and other EQs acquired by reaction path searching (EQ i≠1 ). The structure data includes data indicating the coordinates of atoms for each EQ included in the list of acquired EQs. The atomic coordinates represent the type and position of each atom constituting the structure. The position of an atom is indicated, for example, by the three-dimensional position coordinates of the atom (for example, coordinates on a Cartesian coordinate system). The fragment data is data indicating the combination of atoms in a fragment pair (fragment A and B) set in each EQ for the list of acquired EQs. The memory unit 10 stores the fragment data in association with the EQ. When a ranking (priority in the present invention) has been acquired for the fragment pair, the fragment data includes data indicating the ranking of the fragment pair. The level of the ranking indicates the level of priority for performing a reaction path search. The number of searches data is data indicating the number of searches for an EQ. i The number of times that reaction path searches starting from were performed is n i Each EQ in the list of acquired EQs i(i=1, 2, ...). Furthermore, the fragment data includes data indicating whether or not a reaction path search has been performed for each fragment pair. This data indicates that it has not been performed by default. Furthermore, the memory unit 10 stores data indicating various numerical conditions used by each part of the reaction path searching system 1 when performing the following calculations.

[0036] The fragment pair derivation unit 20 extracts combinations of atoms constituting each of the fragments A and B from all atoms constituting the structure based on the above-mentioned method according to the SC-AFIR method. Fragment data indicating the extraction results is stored in the storage unit 10. ... i A combination of atoms constituting a fragment pair is generated for

[0037] The fragment pair ranking acquisition unit 30 refers to the fragment data in the memory unit 10, and from the list of acquired EQs, identifies EQs for which fragment pair rankings have not been acquired (hereinafter referred to as target EQs). Then, it acquires a ranking for each of the target EQs as follows:

[0038] First, the fragment pair ranking acquisition unit 30 derives the second-order differential coefficient b and the third-order differential coefficient a of E(Q) at Q=Q0 for each of the fragment pairs derived by the fragment pair derivation unit 20 for the target EQ. Here, E(Q) indicates the PES as described above. Q0 indicates the position of the atom in the structure corresponding to the target EQ. The second-order differential coefficient b and the third-order differential coefficient a are calculated based on, for example, the slope of the tangent to the curve C at Q=Q0 and Q=Q' shown in FIG. 4 and the value of E(Q). Q'(=Q0+ΔQ) indicates the position of the atom in the structure after at least one of the positions of the pair of 0th layer atoms is minutely moved relative to the other along a straight line connecting the pair of 0th layer atoms so that the distance between the pair of 0th layer atoms in the fragment pair at Q=Q0 is D (see FIG. 1), and the distance changes minutely from D to D+ΔD (for example, ΔD=0.05 angstroms). Curve C is a curve along the PES that shows the change in E(Q) when the relative positions of the pair of 0th layer atoms are changed along the straight line. The fragment pair ranking acquisition unit 30 optimizes Q' or ΔQ to satisfy the following condition. The condition is that, while fixing the distance between the pair of 0th layer atoms to D+ΔD, for a distance matrix whose elements are distances between all atoms, matrix elements other than the elements corresponding to the 0th layer atoms in the distance matrix at Q=Q' reproduce as closely as possible matrix elements other than the elements corresponding to the 0th layer atoms in the distance matrix at Q=Q0. E(Q) is calculated by a cubic function f(x)=a3*x 3 +a2*x 2 +a1*x+a0, and the coefficients of f(x) are found by solving four simultaneous equations consisting of four formulas: f(x(Q0))=E(Q0), f(x(Q'))=E(Q'), f'(x(Q0))=0, and f'(x(Q'))=(slope of tangent line l to curve C at Q=Q') for a0, a1, a2, and a3. Note that x indicates the coordinate axis along the line connecting a pair of 0th layer atoms, f'(x) is the first derivative of f(x) (df / dx), and f''(x) is the second derivative of f(x) (d 2 f / dx 2). Let f(x(Q)) and f'(x(Q)) denote the values ​​of f(x) and f'(x) at the x corresponding to the position Q. Then, use f(x) with the calculated a0, a1, a2, and a3 substituted in to find a and b as b = f''(Q0) and a = f''(Q0). Note that f''(x) is the third derivative (d 3 f / dx 3 ) Note that the function and step of calculating the second-order differential coefficient b and the third-order differential coefficient a by the fragment pair ranking acquisition unit 30 corresponds to the function of the differential coefficient calculation means and the differential coefficient calculation step in the present invention.

[0039] The fragment pair ranking acquisition unit 30 calculates the ζ of Equations 4-1 to 4-4 expressed by at least one of the second-order differential coefficient b and the third-order differential coefficient a acquired as described above. Second , ζ Third , ζ Scaled and ζ Average The ranking of the fragment pairs is determined based on either Second , ζ Third , ζ Scaled and ζ Average Which of the following is used is selected according to a user input via an input device. Second , ζ Third , ζ Scaled and ζ Average The selected one of these is designated as ζ. At this time, the fragment pair ranking acquisition unit 30 sets a higher ranking for fragment pairs with smaller ζ. The acquired ranking is stored in the storage unit 10.

[0040] (Formula 4-1) TIFF0007655553000011.tif6170

[0041] (Formula 4-2) TIFF0007655553000012.tif6170

[0042] (Formula 4-3) TIFF0007655553000013.tif8170

[0043] (Formula 4-4) TIFF0007655553000014.tif8170

[0044] ζ Second The direction where is small can be said to be the direction where the local increase of E(Q) is small. Second By performing reaction path search for fragment pairs with small , it is expected that we can obtain a path with a small increase in E(Q) near EQ. Third In the direction where is small, E(Q) changes from downward convex to upward convex, and it is expected that a transition state will be found. In this case, the height of the maximum point on the cubic function is determined by the balance between a and b. Third becomes especially small when E(Q) suddenly changes from a downward convex to an upward convex in the direction where b is large. Therefore, ζ Third By performing reaction path search for fragment pairs with small ζ, it is expected that we can obtain paths through which hard bonds such as chemical bonds can be rearranged with low barriers. Scaled The smaller is, the lower the maximum point on the cubic function. Therefore, ζ Scaled By performing reaction path search on small fragment pairs, it is expected that it will be possible to obtain low-barrier pathways such as internal rotations and hydrogen bond rearrangements. In many chemical reactions, low-barrier events such as conformational changes are repeated on short time scales such as nanoseconds or microseconds, achieving local thermal equilibrium. So-called rare events that involve high barriers, such as chemical bond rearrangements, proceed on long time scales such as seconds or hours. In multi-step reactions, the entire reaction progresses through repeated local thermal equilibrium and rare events. To reproduce this situation, the ζ Third Small fragment pairs and ζ Scaled We need to calculate both of the smaller fragment pairs. Third and ζ Scaled The geometric mean of Average has been introduced. Average is, Third Or ζScaled becomes smaller for small fragment pairs. Therefore, ζ Average By performing reaction path searches on fragment pairs with small sizes, it is expected that a variety of pathways, such as internal rotation, hydrogen bond rearrangement, and chemical bond rearrangement, can be obtained.

[0045] The state selection unit 40 (state selection means in the present invention) selects EQ, which is the starting point of the reaction path search. Specifically, the state selection unit 40 selects Λ shown in the above-mentioned formula 4. i Calculate the following formula 5: i Then, the state selection unit 40 calculates v i The loudest EQ i The state selection unit 40 selects Λ as the starting point for the next reaction path search. i When calculating, the rate constant matrix contraction calculation based on the RCMC method is performed as described above. The function of the state selection unit 40 corresponds to the function of the rate constant contraction means in the present invention.

[0046] (Formula 5) TIFF0007655553000015.tif11170

[0047] n i is represented by the search count data stored in the storage unit 10. all will be the same as all the EQs up to that point. i n in i where α and β are real numbers between 0 and 1. For example, α=0.5 and β=0.5 may be set. α is the sum of Λ i Small EQ i Even those EQ i β plays a role in ensuring that EQ is selected several times. i It plays a role in making it difficult for the user to select the i is a random real number between 0 and 1. i is ν based on Equation 5 i It is set for each calculation.

[0048] The reaction path calculation unit 50 executes a reaction path search with the EQ selected by the state selection unit 40 as the starting EQ (first equilibrium state in the present invention). The reaction path calculation unit 50 executes a reaction path search for the fragment pairs indicated by the fragment data stored in the memory unit 10 in association with this EQ. At this time, the reaction path calculation unit 50 executes a reaction path search for the fragment pairs with higher rankings before those with lower rankings. In other words, a reaction path search is executed for the fragment pairs with higher rankings before those with lower rankings. For the fragment pairs for which such a reaction path search has been executed, the fragment data stored in the memory unit 10 is updated to indicate that a reaction path search has been executed.

[0049] In the reaction path search, as described above, calculation of the AFIR path and the LUP path, optimization of the transition state, and calculation of the IRC are performed. As a result, a new reaction path is acquired starting from the EQ selected by the state selection unit 40. The reaction path calculation unit 50 refers to the memory contents of the memory unit 10, and when there is a new EQ that is not included in the list of acquired EQs in the EQ (the second equilibrium state in the present invention) through which the acquired reaction path passes or reaches, adds the new EQ to the list of acquired EQs indicated by the memory contents of the memory unit 10. In addition, structural data indicating the coordinates of the atoms corresponding to the new EQ is added to the memory contents of the memory unit 10. For this new EQ, the fragment pair derivation unit 20 extracts a combination of atoms constituting a fragment pair, and the fragment pair ranking acquisition unit 30 acquires a ranking for the fragment pair derived by this. Note that the function and step of calculating the AFIR path by the reaction path calculation unit 50 correspond to the function of the structural change calculation means and the structural change calculation step in the present invention. In addition, the functions and steps of the reaction path calculation unit 50 that calculates a reaction path corresponding to a transition to another EQ by calculating an AFIR path starting from the EQ selected by the state selection unit 40 correspond to the function of the inter-equilibrium state change calculation means and the inter-equilibrium state change calculation step in the present invention.

[0050] A series of steps in the reaction path searching method executed by the reaction path searching system 1 will be described below with reference to Fig. 5. First, for an initial state EQ1 provided by a user, the fragment pair derivation unit 20 and the fragment pair ranking acquisition unit 30 acquire fragment pairs and their rankings (step S1). The structural data in the initial state EQ1 is provided, for example, by importing a data file generated by the user into the system.

[0051] Next, the reaction path calculation unit 50 executes a reaction path search starting from EQ1 for the fragment pair corresponding to the highest ranking among the rankings acquired in step S1 (step S2). Based on the reaction path acquired in this way, a new EQ, EQ2, is added to the list of acquired EQs.

[0052] Next, the state selection unit 40 selects each EQ in the list of acquired EQs. i For (i=1, 2, ...), ν i Next, the fragment pair ranking acquisition unit 30 calculates v i The EQ corresponding to the largest of the i It is determined whether or not a ranking of the fragment pair has been acquired for the EQ (step S4). This determination is made by referring to the fragment data stored in the memory unit 10. If the fragment pair ranking acquisition unit 30 determines that a ranking has been acquired (step S4, Yes), the process of step S6 is executed. i is the target EQ) (step S4, No), the v i The EQ corresponding to the largest of the i The fragment pair ranking acquisition unit 30 acquires the ranking for the fragment pairs related to (step S5).

[0053] Next, the reaction path calculation unit 50 calculates the v iThe EQ corresponding to the largest of the i The contents stored in the storage unit 10 are referred to, and a reaction path search is performed for the highest-ranked fragment pair among those for which a reaction path search has not yet been performed (step S6). If a new EQ that is not in the list of previously obtained EQs is obtained by the reaction path search, data related to the new EQ is added to the contents stored in the storage unit 10.

[0054] Next, the reaction path search system 1 judges whether or not a search condition for the reaction path is satisfied (step S7). The search condition may be, for example, "in the past n calculations of the AFIR path, the traffic volume Λ i The largest m EQs have not been updated. In this case, the natural numbers n and m may be set by the user. If it is determined that the search condition is not satisfied (step S7, No), the process from step S3 is executed. If it is determined that the search condition is satisfied (step S7, Yes), the reaction path searching system 1 ends the series of processes.

[0055] According to the reaction path searching system 1 of the present invention described above, a reaction path search is performed for fragment pairs with high priority based on at least one of the magnitudes of the second-order differential coefficient b and the third-order differential coefficient a of E(Q) at the position of an atom corresponding to a certain equilibrium state EQ, prior to fragment pairs with low priority. Second , ζ Third , ζ Scaled and ζ Average The fragment pairs are ranked based on the ζ selected by the user. Then, reaction path search is performed in descending order of the ranking of the fragment pairs. Second , ζ Third , ζ Scaled and ζ AverageWhen a reaction path search is preferentially performed for fragment pairs with high rankings based on the above, as described above, it becomes easier to obtain a target reaction path, such as a transition state being easily found, a path in which hard bonds such as chemical bonds are recombined with low barriers being easily obtained, a path with low barriers such as internal rotation and hydrogen bond rearrangement being easily obtained, and a variety of paths such as internal rotation, hydrogen bond rearrangement, and chemical bond rearrangement being easily obtained. This increases the possibility of obtaining a target reaction path without incurring too much calculation cost. Therefore, an increase in calculation cost in accordance with an increase in the number of atoms is easily suppressed.

[0056] In addition, in this embodiment, v i The loudest EQ i is selected and its EQ i A reaction path search is performed for ν i As shown in the above formula 5, n all When is relatively small, α^{n i / log(n all )} is Λ i becomes relatively small for Λ i Therefore, the equilibrium state with a larger traffic volume is more likely to be selected as the calculation target. all As becomes relatively large, α^{n i / log(n all )} is Λ i In other words, Λ i The difference is ν i Therefore, the difference in n all When is relatively large, the ease of selecting the equilibrium state increases. i Therefore, an equilibrium state with a small amount of traffic will also be selected. In this way, with the above configuration, an equilibrium state with a high priority is appropriately selected according to the progress of the calculation.

[0057] [Example] An embodiment of the present invention will be described below. In this embodiment, the reaction path search method shown in FIG. 5 according to the above embodiment was carried out for the Strecker reaction, which is the most basic organic synthesis reaction, with some steps being changed. The Strecker reaction is a reaction in which a ketone (or aldehyde), ammonia, and hydrogen cyanide react to produce an aminonitrile. In this embodiment, acetone was used as the ketone. In addition, one extra molecule of water was added to the system in consideration of the efficiency of proton transfer in an aqueous solution. The collision energy parameter γ of the force applied by the SC-AFIR method was set to 500 kJ / mol. A weak force of γ=100 / {N(N-1) / 2} kJ / mol was added to all atom pairs so that the molecules would not be too far apart (N here is the number of atoms included in the system). In order to reduce the amount of calculation, as the above-mentioned process change, optimization of the transition state after calculation of the LUP path and calculation of the IRC were not performed, and the energy maximum point of the LUP path was used as an approximate transition state, and the LUP path was used as a path defining the connection between EQs. In the harmonic analysis for calculating the Gibbs energy, 50 cm -1 The frequency of the mode with the following frequency is 50 cm -1 The electronic structure calculations were performed at the RHF / SV level using Gaussian 16.

[0058] ν i In the calculation of the values, Λ was obtained by performing the RCMC method at three temperatures: 200K, 300K, and 400K. i The largest value was used. Note that Λ was i If the value becomes zero, the corresponding ν i The value was also forced to zero. The reaction time was 1 day, the pressure was 1 atm, and T R was set to 3000K. Second , ζ Third , ζ Scaled and ζ Average In addition to the case where the fragment pairs to be searched for reaction paths were selected based on the rankings given by using , we also performed reaction path searches when the fragment pairs to be searched for were selected randomly. The search was terminated when the product was obtained, and the computational cost required up to that point was calculated using ζ Second , ζThird , ζ Scaled and ζ Average The calculations were compared between the cases where the same initial structure (reactant structure) and product structure were used. Table 1 shows the comparison results. Table 2 shows the x-, y-, and z-coordinates of each atom in the initial structure, and Table 3 shows the x-, y-, and z-coordinates of each atom in the product structure.

[0059] (Table 1) TIFF0007655553000016.tif53170

[0060] (Table 2) TIFF0007655553000017.tif174170

[0061] (Table 3) TIFF0007655553000018.tif174170

[0062] In Table 1, N gradient , N Hessian , N EQ , N path correspond to the number of gradient calculations performed until the product was obtained, the number of Hessian calculations (calculation related to second-order differentials), the number of EQs obtained until the product was obtained, and the number of paths, respectively. Here, it was confirmed that the reaction mechanism in which aminonitrile is produced by nucleophilic addition of cyano anion to iminium ion was obtained in all calculations.

[0063] As can be seen from Table 1, the calculation cost was the highest when the fragment pairs were randomly selected. In other words, it was confirmed that the present invention is effective. Second , ζ Third , ζ Scaled and ζ Average Among them, ζ Averageperformed best. This is believed to be because both the local thermal equilibration pathway and the slow pathway corresponding to rare events were sampled efficiently. When reactant molecules are weakly associated, there are many pathways in which the molecular orientation changes with a very low barrier of less than 10 kJ / mol. Such pathways are called ζ Scaled On the other hand, the pathway in which a water molecule is eliminated from a hemiaminal intermediate to produce an iminium ion has a relatively high barrier of about 100 kJ / mol, and Scaled If we use ζ as an index, the calculation will be postponed. Scaled and ζ Third The geometric mean of Average If we use ζ as an index, Scaled Small fragment pairs and ζ Third Therefore, it is likely that preferential exploration was also performed for recombination pathways with relatively high barriers to chemical bonding, such as the conversion of a hemiaminal intermediate to an iminium ion.

[0064] In this example, the reaction path search was performed on a system consisting of at most 20 atoms (N=20). On the other hand, the number of fragments increases in proportion to the square of N, so random and ζ Average It is expected that the cost ratio to N will also increase rapidly with N. In this example, we performed numerical verification using an organic synthesis reaction as an example, but in processes that do not involve rearrangement of chemical bonds, such as conformational changes in polymers and phase transitions of aggregate structures, ζ Scaled It is expected that it would be better to prioritize only low-barrier routes using ζ as an index. Average and ζ Scaled It is recommended to use them appropriately.

[0065] [Details of the RCMC Act] The details of the RCMC method are explained. Let the set of states in the n-th contraction be N (n) , the rate constant from state i to j is k i→j (n) , the Boltzmann distribution of state i is P i(n) In other words, N (0) is the set of all stable structures (EQs) included in the reaction path network, k i→j (0) is the rate constant for the transition from stable structure i to j, P i (n) corresponds to the Boltzmann distribution of the stable structure i. The reduction operation using these quantities is as follows:

[0066] 1. Upper limit of the rate constant k MAX (=1 / t MAX ) and set n=0. 2. For all elementary processes, k i→j (0) Calculate the rate constant matrix. In this case, if there are two or more elementary processes directly connecting stable structures i and j, k i→j (0) is the sum of the rate constants for them. If there is no elementary process directly connecting stable structures i and j, k i→j (0) is set to zero. k for the elementary process connecting stable structures i and i i→i (0) is also set to zero. 3. Boltzmann distribution of all stable structures P i (0) Request. 4. (Initial) distribution Q of all stable structures i (0) Gives. 5. Chi ii (0) = 1, χ ji (0) = 0 (j ≠ i). The following steps 6 to 14 are executed repeatedly as a loop. 6. Among all pairs of i and j, k i→j (n) Find the one that is largest. 7. If the maximum k obtained in 6 i→j (n) About i→j (n) <k MAX If so, exit the loop. 8. The steady-state approximation is applied to state i. In this case, all rate constants in the rate constant matrix are corrected using Equations 6 and 7. (Equation 6) TIFF0007655553000019.tif9170(Formula 7) TIFF0007655553000020.tif151709. Set all rate constants associated with state i to zero (k´ i→m (n+1) =k´ m→i (n+1) =0). 10.As shown in Equation 8, the Boltzmann distribution of state i is set to zero (P i (n+1) =0), and then distribute it to another state k (k ≠ i) weighted by the rate constant. (Formula 8) TIFF0007655553000021.tif817011.As shown in Equation 9, the distribution of state i is set to zero (Q i (n+1) =0), and then distribute it to another state k (k ≠ i) weighted by the rate constant. (Equation 9) TIFF0007655553000022.tif817012.The contribution of stable structure j to state k is updated using equation 10. (Equation 10) TIFF0007655553000023.tif917013.To reproduce the updated Boltzmann distribution, k´ k→l (n+1) is corrected according to Equation 11. (Formula 11) TIFF0007655553000024.tif1317014. The remaining N after eliminating state i (weighted distribution) (n+1) Let be the set of states and go back to step 6.

[0067] Q in 4 and 11 above i (0) and Q i (n) Q in Equation 3 (m) j In addition, the χ in 5 and 12 above is used as ii (0) , χji (0) (j ≠ i) and χ ji (n) is the χ of formula 3 ji (m) It is used as.

[0068] <Modification> The above is a description of a preferred embodiment of the present invention, but the present invention is not limited to the above-described embodiment, and various modifications are possible within the scope of the means for solving the problems.

[0069] For example, in this embodiment, ζ Second , ζ Third , ζ Scaled and ζ Average A ranking is set for each fragment pair such that the smaller the ζ selected by the user from any one of the above, the higher the ranking. Then, the reaction path search is performed in order from the highest ranked fragment pair. However, other methods may be adopted for setting the ranking and the order in which the reaction path search is performed. For example, a two-stage ranking may be set for each fragment pair based on whether ζ is equal to or less than a predetermined value. Then, after the reaction path search is performed in a random order for the fragment pairs with higher rankings, the search may be performed in a random order for the fragment pairs with lower rankings. In this case, too, by performing a search for the fragment pairs with higher rankings before the fragment pairs with lower rankings, the desired reaction path is more likely to be obtained in advance compared to the case where a search is performed for all fragment pairs in a random order.

[0070] In the above embodiment, the state selection unit 40 selects v i The loudest EQ i is selected as the starting point for reaction path search. However, such v i The starting point for the reaction path search may be selected based on a criterion other than Λ i Equilibrium states with larger n are more likely to be selected. all The larger Λ iThe starting point for the reaction path search should be selected so that differences in the reaction rate are unlikely to result in differences in the likelihood of the equilibrium state being selected. [Explanation of symbols]

[0071] 1. Reaction Path Search System 10 Storage section 20 Fragment pair derivation section 30 Fragment pair ranking acquisition unit 40 Status selection section 50 Reaction Path Calculation Unit

Claims

1. A program that causes one or more computers to function as a reaction path search system that calculates a reaction path in a structure consisting of a plurality of atoms as a change in structure represented by a positional relationship between the plurality of atoms, The reaction path searching system comprises: N extracted from the plurality of atoms １ pieces (N １ A fragment A consisting of atoms (N: natural number) and a fragment A and another N extracted from the plurality of atoms ２ pieces (N ２ With respect to a fragment pair constituted by a fragment A and a fragment B consisting of atoms of s-th (s: s≦N) atoms, the potential energy in a geometric parameter Q indicating the positions of the plurality of atoms is defined as E(Q), and the s-th (s: s≦N) atom belonging to the fragment A is defined as １ (t: a natural number satisfying N ２ The distance between the atom and the atom is r ｓｔ Let R ｓ and R ｔ are the covalent bond radii of the s-th atom and the t-th atom, respectively, ρ is either +1 or -1, and p and α are constant values, the function F of Equation 1 calculated based on Equation 2 is ＡＦＩＲ A structural change calculation means for calculating a structural change that causes (Q) to have a minimum value; (Formula 1) (Formula 2) A differential coefficient calculation means for calculating at least one of a second-order differential coefficient b and a third-order differential coefficient a of E(Q) at a position of the plurality of atoms corresponding to a first equilibrium state, which is one of the equilibrium states in which E(Q) takes a minimum value, for each of the plurality of fragment pairs having different atomic combinations; A reaction path searching program characterized in that it comprises an equilibrium state change calculation means that causes the structural change calculation means to calculate the structural change in the transition from the first equilibrium state to a second equilibrium state, which is another of the equilibrium states, for a fragment pair among the multiple fragment pairs that has a high priority based on the magnitude of at least one of a and b calculated by the differential coefficient calculation means, prior to other fragment pairs that have a lower priority.

2. The reaction path searching system comprises: a state selection means for selecting, as the first equilibrium state, one of N equilibrium states (N: a natural number of 2 or more) consisting of all the second equilibrium states and an initial state obtained based on past calculation results by the inter-equilibrium-state change calculation means, The reaction path searching program according to claim 1, characterized in that the equilibrium state change calculation means repeatedly causes the structural change calculation means to calculate the structural change for the fragment pair in the transition from the first equilibrium state selected by the state selection means to the second equilibrium state.

3. The reaction path searching system comprises: a rate constant reduction means for obtaining a rate constant matrix consisting of l*l rate constants associated with transitions between l superstates (l: a natural number satisfying l<N) represented by a weighted sum of the N equilibrium states from a rate constant matrix of N rows and N columns consisting of N*N rate constants associated with transitions between the N equilibrium states, by repeating reduction of the rate constant matrix based on an RCMC method m times (m: a natural number satisfying m=N-l); The state selection means: For each of m=1, 2, ..., M (M: a natural number satisfying M < N), the rate constant reduction means acquires a contracted rate constant matrix, and based on the result, each of the N equilibrium states is calculated by EQ. ｉ (i: a natural number equal to or less than N), and each of the l superstates is ｊ (j: a natural number less than or equal to l), and the EQ after m contractions ｉ Population of p ｉ (m), and the SS after m contractions is ｊ Population of Q ｊ (m), and the SS after m contractions is ｊ Eq. ｉ The contribution of χ ｊｉ (m), and EQ. ｉ The relative Gibbs energy of ΔG ｉ The gas constant is R, and the model temperature parameter is T Ｒ When this is done, p calculated based on Equation 3 ｉ (m) and Λ calculated based on Equation 4 ｉ The larger the EQ ｉ 3. The reaction path searching program according to claim 2, wherein one of the N equilibrium states is selected as the first equilibrium state such that the first equilibrium state is more likely to be selected. (Formula 3) (Formula 4)

4. The state selection means: Λ ｉ Large EQ ｉ The more likely it is to be selected, and EQ ｉ The number of times that the structural change calculation means calculates the structural change with the first equilibrium state being n ｉ Small EQ ｉ and the total number of times that the structural change calculation means has calculated the structural change for the N equilibrium states is n ａｌｌ The larger Λ ｉ The reaction path search program according to claim 3, characterized in that one of the N equilibrium states is selected as the first equilibrium state so that differences in the N equilibrium states are unlikely to be reflected in differences in the likelihood of the equilibrium state being selected.

5. The state selection means: The total number of times that the structural change calculation means calculates the structural change for the N equilibrium states is n ａｌｌ Then, EQ. ｉ The number of times that the structural change calculation means calculates the structural change with the first equilibrium state set as n ｉ Then, ξ ｉ is a real number between 0 and 1, and α and β are real numbers between 0 and 1, respectively. ｉ The largest EQ ｉ 5. The reaction path searching program according to claim 3, wherein the reaction path searching program selects as the first equilibrium state. (Formula 5)

6. A system for calculating a reaction path in a structure consisting of a plurality of atoms as a change in structure represented by a positional relationship between the plurality of atoms, comprising: N extracted from the plurality of atoms １ pieces (N １ A fragment A consisting of atoms (N: natural number) and a fragment A and another N extracted from the plurality of atoms ２ pieces (N ２ With respect to a fragment pair constituted by a fragment A and a fragment B consisting of atoms of s-th (s: s≦N) atoms, the potential energy in a geometric parameter Q indicating the positions of the plurality of atoms is defined as E(Q), and the s-th (s: s≦N) atom belonging to the fragment A is defined as １ (t: a natural number satisfying N ２ The distance between the atom and the atom is r ｓｔ Let R ｓ and R ｔ are the covalent bond radii of the s-th atom and the t-th atom, respectively, ρ is either +1 or -1, and p and α are constant values, the function F of Equation 1 calculated based on Equation 2 is ＡＦＩＲ A structural change calculation means for calculating a structural change that causes (Q) to have a minimum value; (Formula 1) (Formula 2) A differential coefficient calculation means for calculating at least one of a second-order differential coefficient b and a third-order differential coefficient a of E(Q) at a position of the plurality of atoms corresponding to a first equilibrium state, which is one of the equilibrium states in which E(Q) takes a minimum value, for each of the plurality of fragment pairs having different atomic combinations; A reaction path searching system characterized in that it comprises an equilibrium state change calculation means that causes the structural change calculation means to calculate the structural change in the transition from the first equilibrium state to a second equilibrium state, which is another of the equilibrium states, for a fragment pair among the multiple fragment pairs that has a high priority based on the magnitude of at least one of a and b calculated by the differential coefficient calculation means, prior to other fragment pairs that have a lower priority.

7. A method for calculating a reaction path in a structure consisting of a plurality of atoms as a change in structure represented by a positional relationship between the plurality of atoms, comprising the steps of: N extracted from the plurality of atoms １ pieces (N １ A fragment A consisting of atoms (N: natural number) and a fragment A and another N extracted from the plurality of atoms ２ pieces (N ２ With respect to a fragment pair constituted by a fragment A and a fragment B consisting of atoms of s-th (s: s≦N) atoms, the potential energy in a geometric parameter Q indicating the positions of the plurality of atoms is defined as E(Q), and the s-th (s: s≦N) atom belonging to the fragment A is defined as １ (t: a natural number satisfying N ２ The distance between the atom and the atom is r ｓｔ Let R ｓ and R ｔ are the covalent bond radii of the s-th atom and the t-th atom, respectively, ρ is either +1 or -1, and p and α are constant values, the function F of Equation 1 calculated based on Equation 2 is ＡＦＩＲ A structural change calculation step of calculating a structural change in which (Q) takes a minimum value; (Formula 1) (Formula 2) A differential coefficient calculation step of calculating at least one of a second-order differential coefficient b and a third-order differential coefficient a of E(Q) at a position of the plurality of atoms corresponding to a first equilibrium state, which is one of the equilibrium states in which E(Q) takes a minimum value, for each of the plurality of fragment pairs having different atomic combinations; A reaction path searching method characterized in that it comprises an equilibrium state change calculation step in which, for a fragment pair among the multiple fragment pairs having a high priority based on the magnitude of at least one of a and b calculated in the differential coefficient calculation step, the structural change in the transition from the first equilibrium state to a second equilibrium state, which is another of the equilibrium states, is calculated by the structural change calculation step prior to other fragment pairs having a lower priority.

Images