JP7443712B2 - How to search for diffusion routes - Google Patents

Description

The present invention relates to a diffusion route search method.

With the aim of further improving the performance of various materials, the search for new materials and the selection and search for substitution elements that replace some of the elements constituting substances are actively being carried out.

For example, in Patent Document 1, as a method for selecting a substitution element for a lithium composite oxide, when nickel in a lithium composite oxide LiMO ₂ (M includes nickel) is replaced with a candidate element, the candidate element is accommodated in the nickel site. a solid solubility determination step for determining whether the
Is the adsorption energy Ea of water molecules on the crystal surface when nickel in the lithium composite oxide is replaced with the candidate element smaller than the adsorption energy Eb of water molecules on the crystal surface of the lithium composite oxide before replacement with the candidate element? an adsorption energy determination step for determining the
Substitution in which the candidate element is substituted for nickel in the lithium composite oxide when it is determined in the solid solubility determination step that the candidate element is accommodated in the nickel site, and in the adsorption energy determination step it is determined that Ea<Eb. Disclosed is a method for selecting a substitution element for a lithium composite oxide, which includes a selection step of selecting an element.

Japanese Patent Application Publication No. 2018-140917

As mentioned above, when selecting and searching for new materials and substitution elements, it is important to consider the routes through which atoms that affect the desired reaction, function, etc. move within the crystal of the desired substance. , it is preferable to accurately understand how the particles will spread.

However, the elements inside the crystal are often tightly packed, and there appear to be only small gaps between the atoms that make up the crystal. For this reason, when we estimate the size of a diffused atom, which is the atom whose diffusion path we want to investigate, using the atomic radius or van der Waals radius, we can compare the sizes of multiple gaps in the crystal and determine the diffusion path of the diffused atom. It was difficult to identify.

Furthermore, due to electrostatic interactions between the diffusing atoms and surrounding atoms, the diffusing atoms do not necessarily pass through a wide space within the crystal. Therefore, as described above, it has been difficult to specify the diffusion path of the diffusing atoms based on the atomic radius of the diffusing atoms and the size of the gaps within the crystal.

In view of the problems of the above-mentioned conventional techniques, one aspect of the present invention aims to provide a method for searching for diffusion paths that can efficiently search for diffusion paths of atoms within a crystal.

According to one aspect of the present invention to solve the above problems,
an initial structure setting step of setting an initial structure in which a plurality of atoms other than diffusion atoms are arranged in a crystal used for diffusion route search;
a calculation step of setting a three-dimensional lattice in the initial structure and performing molecular dynamics calculations on the structure in which diffused atoms are arranged at lattice points of the lattice;
a repeated step of changing the lattice points at which the diffused atoms are placed and performing molecular dynamics calculations;
Based on the calculation results of the calculation step and the repetition step, the existence probability of the diffused atom at a plurality of positions in the crystal is determined, and the free energy of the diffused atom at the plurality of positions is determined from the existence probability. a free energy calculation step;
The present invention provides a method for searching for a diffusion path, comprising: a diffusion path searching step of searching for a diffusion path of the diffusing atom based on the free energy at the plurality of positions of the diffusing atom obtained in the free energy calculating step.

According to one aspect of the present invention, it is possible to provide a method for searching for diffusion paths that can efficiently search for diffusion paths of atoms within a crystal.

FIG. 2 is a schematic diagram showing the initial structure of LiMn ₂ O ₄ in which the diffusion path was searched in Example 1. FIG. FIG. 2 is a schematic diagram showing the diffusion path of Li atoms in LiMn ₂ O ₄ determined in Example 1. FIG. 3 is a schematic diagram showing the diffusion path of Li atoms in LiMn ₂ O ₄ determined in Comparative Example 1.

Hereinafter, modes for carrying out the present invention will be described. However, the present invention is not limited to the following embodiments, and various modifications and variations may be made to the following embodiments without departing from the scope of the present invention. Substitutions can be added.

The diffusion route search method of this embodiment can include the following steps.
An initial structure setting step for setting an initial structure in which multiple atoms other than diffusion atoms are arranged, which is included in the crystal used for diffusion route search.
A calculation process in which a three-dimensional lattice is set as the initial structure and a molecular dynamics calculation is performed on a structure in which diffused atoms are placed at the lattice points of the lattice.
An iterative process of changing the lattice points where diffused atoms are placed and performing molecular dynamics calculations.
A free energy calculation step of calculating the existence probability of the diffused atom at multiple positions in the crystal based on the calculation results of the calculation step and the repetition step, and calculating the free energy of the diffused atom at each of the multiple positions from the existence probability.
Diffusion path searching step for searching the diffusion path of the diffusing atoms from the free energies at multiple positions of the diffusing atoms obtained in the free energy calculation step.

The inventor of the present invention has conducted extensive studies on an efficient method for searching for the diffusion path of atoms within a crystal.

Molecular dynamics method is a method for analyzing the movement and energy of atoms and molecules. For this reason, it may be possible to use molecular dynamics calculations as a method of searching for the diffusion paths of atoms within the crystal. However, there is a large difference between the time scale of the diffusion process and the time scale of molecular dynamics calculations, and calculating the diffusion process using molecular dynamics calculations requires a very long calculation time and is extremely expensive. Become. Furthermore, since molecular dynamics calculations depend on initial conditions, calculations from a small number of initial conditions do not necessarily lead to the diffusion path. Therefore, it is difficult to search for diffusion paths of diffusing atoms within a crystal using only molecular dynamics calculations.

When a diffusing atom whose diffusion path is to be investigated is located at a coordinate position within the diffusion path, the free energy of the diffusing atom at that coordinate position is greater than the free energy of the diffusing atom when the diffusing atom is outside the diffusion path. will also become smaller. Therefore, the inventor of the present invention focused on the fact that the diffusion path of the diffusing atoms can be searched by searching for a path where the free energy of the diffusing atoms becomes small.

Therefore, the inventor of the present invention further investigated a method for efficiently searching for a path in which the free energy of the diffusing atoms becomes small. The probability ρ(r) that a diffused atom exists at a specific position r and the free energy F(r) at that position have a relationship expressed by the following equation (1).

F(r)=-k _B Tln(ρ(r))+C...(1)
In addition, k _B in the above formula (1) represents Boltzmann's constant, T represents temperature (K), and C represents a constant.

Therefore, if the probability of existence of a diffusing atom at each position within the crystal is determined, the free energy at each position within the crystal can be determined, and the diffusion path can be searched based on the determined free energy.

For example, by using a large number of results of short-time molecular dynamics calculations, it is possible to collect a large number of samples in which the positions of diffused atoms within the crystal have changed. Therefore, it was discovered that the existence probability of diffusing atoms in a crystal can be calculated using the results of such molecular dynamics calculations, which reduces the cost of implementing the molecular dynamics method and allows efficient derivation of diffusion paths. , completed the present invention.

Each step will be explained below.
(Initial structure setting process)
In the initial structure setting step, it is possible to set an initial structure in which a plurality of atoms other than the diffusion atoms are arranged, which are included in the crystal used for the diffusion route search. That is, it is possible to set the initial coordinates of a plurality of atoms other than the diffusion atoms included in the crystal used for diffusion route search.

In the initial structure setting step, there is no particular limitation on the specific method for setting the positions of a plurality of atoms other than the diffusion atoms included in the crystal used for the diffusion route search. For example, the atomic arrangement of each atom can be set based on the crystal structure of the crystal, which has been experimentally determined or disclosed in literature, etc., and can be used as an initial structure.
(calculation process)
In the calculation step, a three-dimensional lattice is set in the initial structure whose positions were set in the initial structure setting step, and molecular dynamics calculations can be performed on the structure in which diffused atoms are arranged at the lattice points of the lattice.

In the calculation process, a structure in which diffused atoms are arranged in the initial structure can be created as described above. At this time, it is preferable to create the initial structure so that there are many variations in the positions of the diffused atoms included in the structure obtained after the calculation step or the repeating step described below. For this reason, a three-dimensional lattice corresponding to the initial structure is virtually set up for the initial structure so that the positions of the diffusing atoms in the system containing the diffusing atoms used for molecular dynamics calculations are scattered, and the lattice points of the diffusing atoms are Preferably, the atoms are arranged. The size, shape, etc. of the lattice are not particularly limited as long as they can be selected so that the positions of the diffused atoms are dispersed as described above.

Molecular dynamics calculation is a computer simulation method of the physical movement of atoms, and by numerically solving Newton's equation of motion, it is possible to determine the evolution of the position of atoms over time.

In molecular dynamics calculations, information about atoms and interactions between atoms is represented by a functional form for describing potential energy and its parameter set (force field).

The type of force field used in the molecular dynamics calculation in the calculation process is not particularly limited, and various force fields can be used. For example, metal/alloy systems include EAM and MEAM, inorganic compound systems include Buckingham, BKS, Clay-FF, CVFF_aug, etc., semiconductor systems include Tersoff, and organic compound systems include PCFF, Compass, MMFF, OPLS-AA, AMBER, CHARMM, etc. UFF etc. can be used. In addition, forces selected from existing force fields such as polarization force fields such as X-Pol, AMBER polarization force field, and CHARMM polarization force field, and reaction force fields such as ReaxFF, as well as force fields created by yourself as necessary. You can use the field.

Existing force fields may not specify the charge of the target atom. In that case, RESP (Restrained ElectroStatic Potential) charges, AM1-BCC (Bond Charge Correction) charges, etc. can also be used.

The program (software) used for molecular dynamics calculations is not particularly limited, but may be selected from existing programs such as LAMMPS, DL_POLY, Gromacs (Groningen Machine for Chemical Simulations), AMBER, CHARMM, NAMD, or self-made programs. A program can be used.

The setting environment for performing molecular dynamics calculations can be, for example, vacuum or periodic boundary conditions if a solvent is included.

The numerical integration method for solving Newton's equation of motion when performing molecular dynamics calculations is not particularly limited, but may be selected from, for example, the Berlet method, the velocity Berlet method, the Leap-frog method, the predictor-modifier method, etc. Other methods can be used.

The time span for performing molecular dynamics calculations is not particularly limited. However, in the diffusion route search method of this embodiment, it is preferable to perform a plurality of short-time molecular dynamics calculations and collect a plurality of samples in which the positions of diffusing atoms within the crystal have changed. Therefore, it is not necessary to perform long-time molecular dynamics calculations, and it is preferable to select the time width so that calculation costs can be suppressed. The time width for performing the molecular dynamics calculation can be, for example, 0.5 fs or more and 2 fs or less.

Further, the temperature control method is not particularly limited, but may be selected from, for example, the rate scaling method, the Nose-Hoover hot bath method, the Nose-Hoover chain method, the Berendsen hot bath method, the Andersen hot bath method, the Langevin kinetic method, etc. Other methods can be used.

The method of controlling the pressure under periodic boundary conditions is not particularly limited either, but for example, a method selected from the Berendsen method, the Parinello-Rahman method, etc. can be used.

A cutoff method can be used to calculate long-range interactions such as electrostatic interactions and van der Waals interactions. In particular, the Particle-Mesh Ewald method, the multipole expansion method, etc. can be used to calculate electrostatic interactions under periodic boundary conditions.

Molecular dynamics calculation in the calculation process uses, for example, a CPU (Central Processing Unit), a RAM (Random Access Memory), various storage media such as a hard disk, an output device such as a display, an input device such as a keyboard, various peripheral devices, etc. It can be carried out using a conventional computer system equipped with a computer system. Note that examples of the computer system include a network server, a workstation, a personal computer, and the like.

Specifically, for example, a program for the molecular dynamics calculation described above is stored in a storage medium, and the program is executed by a CPU, and an input device such as a keyboard or the like is stored in a storage medium such as a RAM. This can be achieved by reading the initial structure and conditions input from .

By performing molecular dynamics calculations, the positions of the diffused atoms and the multiple atoms set in the initial structure setting process can be changed over time.
(repetitive process)
In the iterative process, molecular dynamics calculations can be performed by changing the lattice points at which the diffusing atoms are placed. Specifically, among the lattice points of the three-dimensional lattice set in the calculation process, diffuse atoms are placed at lattice points that are different from the lattice points on which molecular dynamics calculations have already been performed, and the process is performed in the same way as in the calculation process. Molecular dynamics calculations can be performed using In the iterative process, for example, diffusing atoms can be placed at all lattice points of a set three-dimensional lattice, and the process can be repeated until molecular dynamics calculations are performed. Alternatively, for example, diffusion atoms may be placed at predetermined lattice points among the lattice points of a three-dimensional lattice, and the calculations may be repeated until the molecular dynamics calculation is performed.
(Free energy calculation process)
In the free energy calculation step, first, based on the calculation results of the calculation step and the repetition step, the existence probabilities of diffused atoms at a plurality of positions within the crystal can be determined. Then, the free energies of the diffused atoms at the plurality of positions can be calculated from the existence probabilities.

Now consider a system that is in equilibrium at a certain temperature. This system assumes that all structural spaces can be divided into discrete states. If we observe the time evolution of this system in discrete time, the series of evolution of the system can be represented by the transitions between the states visited by the system. The observed state can be regarded as one state of a stochastic process in which time is discretized. If this process is described by a Markov chain, the probability of observing each state at a certain time depends only on the previous state. For a process with L states, the process is characterized by an L×L transition matrix T(τ) that depends only on the unit time τ, which is the interval of observation times.

The probability vector occupying L states at time t is denoted by p(t). If the probability vector representing the initial state is given by p(0), the probability vector after time nτ is expressed by the following equation (2).

p(nτ) = T(nτ) p(0) = [T(τ)] ⁿ p(0)...(2)
The transition matrix T(τ) has an eigenvalue of μ _k (τ) and a corresponding eigenvector μ _k . Each eigenvalue has a time scale τ _k =−τ/ln (μ _k (τ)) called the implemented time scale, and the associated eigenvector provides information on state transitions according to this time scale.

If this system can be described as a Markov process, then when the transition matrix T is obtained, the relationship between the state p(t) at time t and the state p(t+τ) at time (t+τ) is p(t+τ)= It can be expressed as Tp(t).

When the system reaches an equilibrium state, the above equation becomes π=Tπ. Here, π represents a probability distribution in an equilibrium state, and Σ _i π _i =1. The above equation indicates that the probability distribution in the equilibrium state corresponds to an eigenvector of T with an eigenvalue of 1.

The free energy in the equilibrium state is obtained from the probability distribution in the equilibrium state, as shown in equation (3) below.

なお、上記式中のＦ_ｉは状態ｉの自由エネルギーを、ｋ_Ｂはボルツマン定数を、Ｔは系の絶対温度を、π_ｉは状態ｉの確率をそれぞれ表している。

In the above formula, F _i represents the free energy of state i, k _B represents the Boltzmann constant, T represents the absolute temperature of the system, and π _i represents the probability of state i.

However, in the above equation, the free energy of the value with the largest probability distribution is set to 0.

Therefore, in the free energy calculation step of the diffusion route search method of this embodiment, first, the structure of the system containing the diffusing atoms calculated in the calculation step and the repetition step is extracted at every arbitrary time interval Δt and clustered. However, Δt≦τ. For example, Δt can be set to 0.5 ps or more and 5 ps or less. Although the specific clustering method is not particularly limited, for example, the K-Means method, the K-Centers method, etc. can be used. By performing clustering, the structure of a system containing diffused atoms can be classified into multiple states.

By classifying the structure of a system containing diffused atoms by clustering, it is possible to estimate the unit time τ and count the transitions between the classified states. From the obtained number of transitions, the previously described transition matrix T(τ) can be estimated using, for example, the maximum likelihood method.

From the obtained transition matrix T(τ), the eigenvalues and eigenvectors of the transition matrix T(τ) are determined. The eigenvector at this time represents the existence probability of each classified state. That is, by determining the eigenvectors, the existence probabilities of the diffused atoms at multiple positions within the crystal can be determined. Then, using the calculated eigenvector, that is, the existence probability, the free energy for the position of the diffused atom in each classified state can be calculated using the above-mentioned equation (3).
(Diffusion route search process)
In the diffusion path searching step, the diffusion path of the diffusing atom can be searched from the position of the diffusing atom in each state and the corresponding free energy obtained in the free energy calculation step. That is, the diffusion path of the diffusing atoms can be searched from the free energies at multiple positions of the diffusing atoms obtained in the free energy calculation step.

In the diffusion path searching step, the coordinates of the diffusing atoms in each classified state and the free energy calculated in the free energy calculation step in each state can be mapped, that is, three-dimensionally plotted. From the obtained map, a continuous region with small free energy can be identified as the diffusion route of the diffusing atoms.

Note that if an appropriate diffusion path cannot be found, the setting of the unit time τ used in the calculation process or the repetition process may be inappropriate. Therefore, if an appropriate diffusion path cannot be found, the unit time τ can be changed and the free energy calculation step and the diffusion path searching step described above can be performed again.

According to the diffusion route search method of the present embodiment described above, an initial structure is used in which the positions of diffused atoms are changed throughout the crystal. For this reason, the diffused atoms can reach spatial regions that are difficult to reach by thermal fluctuations based on calculations, so it is possible to increase the spatial search efficiency and search for diffusion paths without overlooking them.

Furthermore, since it is not necessary to perform molecular dynamics calculations for a long time, calculation costs can be suppressed.

Hereinafter, the present invention will be described in more detail with reference to Examples. However, the present invention is not limited to the following examples.
[Example 1]
The following procedure was used to search for the diffusion route of Li atoms in LiMn ₂ O ₄ .
(Initial structure setting process)
The initial structure of LiMn ₂ O ₄ was set. Specifically, as shown in FIG. 1, a structure 10 of LiMn ₂ O ₄ in which lithium atoms 11, manganese atoms 12, and oxygen atoms 13 were arranged was set in the cell. Then, the lithium atom 111 arranged in the center of the cell, which becomes a diffusion atom, was deleted to obtain an initial structure. In addition, in FIG. 1, atoms with the same hatching indicate atoms of the same type.
(calculation process)
A lattice with a lattice spacing of 0.5 Å, that is, a three-dimensional lattice in which a plurality of cubic cells each having a side length of 0.5 Å are arranged, was set in the initial structure obtained in the initial structure setting step. Then, a lithium atom, which is a diffusion atom, was arranged at any one lattice point of the lattice.

Next, molecular dynamics calculations were performed on the structure in which lithium atoms, which are diffusion atoms, are arranged.

Molecular dynamics calculations were performed using LAMMPS as software and a force field developed by Ishizawa et al. of Nagoya University. Then, the coordinates of each atom were entered to set the environment in the crystal.

In addition, the velocity Berlet method was used as a velocity calculation method when performing molecular dynamics calculations, and the time width was set to 1 fs. The Nose-Hoover chain method was used as a temperature control method, and the set temperature was 300K.

The Particle-Mesh Ewald method was used to calculate long-range interactions.

Molecular dynamics calculations were performed for 500 ps under the above conditions, and the positions of all atoms were saved every 0.1 ps.
(repetitive process)
Among the lattice points of the three-dimensional lattice set in the calculation process, a lithium atom as a diffusion atom was already placed, and one lithium atom was placed at a lattice point different from the point where the molecular dynamics calculation was performed. Then, molecular dynamics calculations were performed on the structure in which the lithium atoms were arranged in the same manner as in the calculation process.

Furthermore, except for changing the lattice points where lithium atoms, which are diffusion atoms, are placed, the procedure was repeated in the same manner as described above, until lithium atoms were placed at all lattice points and molecular dynamics calculations were performed. .
(Free energy calculation process)
In the free energy calculation process, first, the structure of the system containing diffused atoms recorded in the calculation process and the repetition process is extracted so that the interval is 0.1 ps, and clustered using the K-Means method to calculate the structure of the system containing diffused atoms. was classified.

After classifying the structure of the system containing diffused atoms by clustering, the unit time τ was set to 10 ps, transitions between the classified states were counted, and the transition matrix T(τ) was estimated using the maximum likelihood method.

From the obtained transition matrix T(τ), the eigenvalues and eigenvectors of the transition matrix T(τ) were determined. Then, using the calculated eigenvectors, the free energy for the position of the diffused atom in each classified state was calculated using the above-mentioned equation (3).
(Diffusion route search process)
The coordinates of the diffused atoms in each classified state and the free energy calculated in the free energy calculation process for each state are mapped, that is, three-dimensionally plotted, and from the resulting map, a continuous region with small free energy is searched. , as a diffusion route.

FIG. 2 shows the diffusion paths of the diffused atoms obtained in the diffusion path search step.

As shown in FIG. 2, diffusion paths 211 to 214 were found as diffusion paths for the lithium atoms 111. These diffusion paths are the diffusion paths for the lithium atoms 111, which are diffusion atoms placed in the center of the cell (see FIG. 1). ) and lithium atoms arranged around the lithium atom 111.

These diffusion paths 211 to 214 match the previously reported diffusion paths of lithium atoms in LiMn ₂ O ₄ , and the method of searching for diffusion paths used in this example is based on actual phenomena. I was able to confirm that.
[Comparative example 1]
(Initial structure setting process)
As in Example 1, the initial structure of LiMn ₂ O ₄ was set. Specifically, as shown in FIG. 1, a structure 10 of LiMn ₂ O ₄ in which lithium atoms 11, manganese atoms 12, and oxygen atoms 13 were arranged was set in the cell. In this comparative example, structure 10 in which lithium atoms 111, which are diffusion atoms, were also arranged was used as the initial structure.
(Diffusion route search process)
A diffusion route was searched based on the position of each atom and the atomic radius of the lithium atom 111, which is a diffusing atom.

As a result, the diffusion route 31 could be searched as shown in FIG.

As is clear from the comparison with the results of the above examples, it was not possible to search all the diffusion routes in this comparative example. In addition, since the diffusion route is searched based on the position of each atom and the atomic radius of the lithium atom 111, it takes a lot of time to search for the diffusion route, and it is difficult to search for the diffusion route efficiently as in Example 1. I couldn't.

Claims (1)

an initial structure setting step of setting an initial structure in which a plurality of atoms other than diffusion atoms are arranged in a crystal used for diffusion route search;
a calculation step of setting a three-dimensional lattice in the initial structure and performing molecular dynamics calculations on the structure in which diffused atoms are arranged at lattice points of the lattice;
a repeated step of changing the lattice points at which the diffused atoms are placed and performing molecular dynamics calculations;
Based on the calculation results of the calculation step and the repetition step, the existence probability of the diffused atom at a plurality of positions in the crystal is determined, and the free energy of the diffused atom at the plurality of positions is determined from the existence probability. a free energy calculation step;
A method for searching a diffusion path, comprising: a diffusion path searching step of searching for a diffusion path of the diffusing atom from the free energy at the plurality of positions of the diffusing atom obtained in the free energy calculating step.