JP2007004530A - Molecular simulation method - Google Patents

Abstract

<P>PROBLEM TO BE SOLVED: To provide a molecular simulation method with which the phenomenon analysis of a large-scale system can be highly, accurately and quickly performed. <P>SOLUTION: The electronic states of a model constituted of a first atom group composed of the prescribed number of atoms, a second atom group arranged on the periphery of the first atom group and a third atom group arranged on the periphery of the second atom group are calculated to prepare a potential map indicating between the positions of atoms and energy or force. In a phenomenon analytical procedure, phenomenon analysis is performed by referring to the potential map. Since calculation required for the calculation of atomic energy or the like can be reduced by using the potential map in the molecular simulation method, the calculation time can be shortened. Further since the potential map is prepared by suitably using electronic state calculation while considering the effects of atoms arranged on far positions in average, the phenomenon analysis can be highly and accurately performed. <P>COPYRIGHT: (C)2007,JPO&INPIT

Description

  The present invention relates to a molecular simulation method.

  Conventionally, molecular simulation methods for analyzing phenomena such as phase transition and adsorption / desorption of a system composed of a plurality of atoms are known. For example, Patent Document 1 describes a molecular structure transition simulation method for calculating a point defect diffusion path and barrier energy in crystalline silicon of a semiconductor device. Further, Patent Documents 2 to 6 describe a molecular simulation method using a molecular dynamics method or the like.

JP 2002-260975 A Japanese Patent No. 2950253 Japanese Patent No. 3110310 Japanese Patent No. 2792487 JP-A-5-274277 JP-A-5-151190

  However, in the techniques described in Patent Documents 1 to 6 described above, the scale of the model is reduced because the application range is limited due to the existence of interatomic potential and the accuracy, and the calculation load is enormous. In some cases, only a short time local phenomenon analysis can be performed.

  The present invention has been made to solve such problems, and the object of the present invention is to provide a molecular simulation method capable of performing phenomenon analysis with high accuracy and high speed for a large-scale system. Is to provide.

  In one aspect of the present invention, a first atomic group composed of a predetermined number of atoms, a second atomic group located around a region where the first atomic group exists, the first atomic group, and the first atomic group A potential map that shows the relationship between the position of an atom and energy or force is created by performing electronic state calculations on a model consisting of a third atomic group located around the region where the two atomic group exists. And a phenomenon map analysis procedure for performing a phenomenon analysis with reference to the generated potential map.

  The above-described molecular simulation method is suitably used for analyzing phenomena such as phase transition and adsorption / desorption of a system composed of a plurality of atoms. In the potential map creation procedure, a potential map showing the relationship between the position of an atom and energy or force by performing electronic state calculation on a model composed of a first atomic group, a second atomic group, and a third atomic group. Create In the phenomenon analysis procedure, the phenomenon analysis is performed with reference to the created potential map. Specifically, the first atomic group is composed of a predetermined number of atoms and is generally located at the central portion in the model. The second atomic group is located around a region where the first atomic group exists. Specifically, the second atomic group is located at a position at least a predetermined distance away from the central atom of the first atomic group. On the other hand, the third atomic group is located around a region where the first atomic group and the second atomic group exist. As described above, according to the molecular simulation method described above, the calculation time required for atomic energy can be reduced by using the potential map, so that the calculation time can be shortened. Furthermore, since the potential map is created by appropriately using the electronic state calculation while taking into account the effects of atoms located far away, it is possible to analyze the phenomenon with high accuracy. Therefore, according to the molecular simulation method described above, it is possible to analyze a phenomenon with high accuracy and high speed for a large-scale system.

  In a preferred embodiment, the first atomic group and the third atomic group are fixed, and the second atomic group is moved by the influence of the first atomic group and the third atomic group.

  In this case, since the atoms constituting the first atomic group are located at the discrete points of the potential map, the arrangement positions are fixed. On the other hand, the atoms constituting the second atomic group are moved according to the arrangement position of the first atomic group in the potential map. That is, the atoms constituting the second atomic group are optimized in the arrangement position in the potential map creation procedure. On the other hand, the atoms constituting the third atomic group affect surrounding atoms, but the arrangement positions are not moved.

  In one aspect of the molecular simulation method, the electronic state calculation is performed by band calculation or molecular orbital calculation.

  In this aspect, in the potential map creation procedure, the potential map is created by performing band calculation or molecular orbital calculation. As a result, even if the electronic state changes due to recombination of atoms due to reaction or the like, it can be treated as a classical potential.

  In another aspect of the molecular simulation method, the phenomenon analysis procedure includes a procedure for obtaining energy or force of an atom by performing spline interpolation on the potential map.

  In this aspect, the phenomenon analysis procedure uses the energy (potential energy) or force obtained from the potential map by performing spline interpolation. When the potential map is created from discrete points, the atoms for which the potential energy and the like are to be obtained are not necessarily located on the discrete points. Therefore, the potential energy at a position in the vicinity of this atom is read from the potential map, and a value obtained by spline interpolation of the read potential energy is set as a potential energy to be obtained. Thereby, even if the number of data is reduced by representing the potential map with discrete points, the accuracy of the phenomenon analysis can be ensured.

  Preferred embodiments of the present invention will be described below with reference to the drawings.

[Potential map creation method]
First, a method for creating a potential map according to an embodiment of the present invention will be described.

  FIG. 1 is a diagram for explaining a basic concept of an electronic state calculation model used for creating a potential map.

  This model is composed of an atomic group located approximately in the center of the cube in FIG. 1 (hereinafter referred to as “first atomic group 1”) and an atomic group located around the region where the first atomic group 1 exists. (Hereinafter referred to as “second atomic group 2”) and an atomic group located around the region where the second atomic group 2 exists (hereinafter referred to as “third atomic group 3”). . The first atomic group 1 is configured by the number of atoms about the coordination number of the crystal (four in the example of FIG. 1), and is located at each discrete point in the potential map described later. In other words, a potential map is created with the first atomic group 1 as the center. Further, the second atomic group 2 is constituted by an atomic group located at a distance that is approximately the second or third proximity of the crystal from the central atom of the first atomic group 1. On the other hand, the third atomic group 3 is constituted by an atomic group located at a distance further away from the central atom of the first atomic group 1 than the distance of the second or third proximity.

  In the present embodiment, a potential map is created by performing electronic state calculation on a model constituted by the atomic group as described above. Specifically, since the atoms constituting the first atomic group 1 are positioned at each discrete point of the potential map, the arrangement position is fixed. That is, the first atomic group 1 is not moved by the force of surrounding atoms. On the other hand, the atoms constituting the second atomic group 2 are moved according to the arrangement position of the first atomic group 1 and the like in the potential map. That is, the atoms constituting the second atomic group 2 move by receiving the forces of surrounding atoms and the like, and the arrangement position is optimized in the potential map creation procedure. On the other hand, the atoms constituting the third atomic group 3 affect surrounding atoms, but the arrangement positions are not moved.

  Here, a specific example of the electronic state calculation described above will be described with reference to FIG.

  FIG. 2A is a diagram for explaining band calculation. The first atomic group 1 is located in the region indicated by the broken line 11 (hereinafter, this region is also referred to as “cluster”). Further, the second atomic group 2 is located in the region indicated by the broken line 12 excluding the region indicated by the broken line 11. In this case, the second atomic group 2 forms a crystal lattice as indicated by a broken line 15. On the other hand, the third atomic group 3 is located in the region indicated by the solid line 13 excluding the regions indicated by the broken lines 11 and 12. According to the band calculation, the total energy can be obtained by calculating the force and energy by imposing periodic boundary conditions on various atomic arrangements. In this case, all atoms in the region surrounded by the solid line 13 constitute a supercell in the band calculation. That is, the band calculation is performed by imposing a periodic boundary condition on the unit in which the central portion of the supercell 13 in which the crystal lattice 15 is appropriately stacked is replaced with the cluster 11. A potential map is obtained by mapping the obtained total energy to the atomic arrangement of the first atomic group 1.

  FIG. 2B is a diagram for explaining molecular orbital calculation (MO; MOPAC). Also in FIG. 2B, the first atom group 1, the second atom group 2, and the third atom group 3 are located at the same arrangement positions as in FIG. The molecular orbital calculation is different from the band calculation in that the above unit is regarded as an isolated molecule and the total energy is obtained without using the periodic boundary condition. Also in this case, the potential map is obtained by mapping the obtained total energy to the atomic arrangement of the first atomic group 1.

  Next, a potential map creation method will be described using the flowchart shown in FIG. Here, as an example, a method for creating a potential map for single crystal growth of SiC by the sublimation method will be described. The potential map creation method is executed by a personal computer or the like. Specifically, the personal computer creates a potential map by executing a predetermined program.

  First, in step S101, the number N of atoms included in the cluster (that is, the number of the first atomic group 1) and the cluster size (force partition distance) are set. For example, when a model having a gas and a crystal is adopted, the number N of atoms in the gas phase is “3” or less, and the number N of atoms in the crystal layer is greater than “3” and less than “6”. When the setting in step S101 is completed, the procedure proceeds to step S102. In step S102, the element type included in the system to be simulated is set, and the process proceeds to step S103.

  In step S103, it is determined whether the number of atoms N is greater than “3”. This determination is performed to distinguish between the case where the target of the potential map to be created is the gas phase (N ≦ 3) and the case where the target is the crystal phase (N> 3). When the number of atoms N is larger than “3” (step S103; Yes), the process proceeds to step S120, and when the number of atoms N is “3” or less (step S103; No), the process proceeds to step S104.

  In step S104, “0” is substituted for variable i, and the process proceeds to step S105. In step S105, the interatomic distance r defined in the atomic group in which the number of atoms N is “2” or more is obtained by the arithmetic expression “r = r0 + rd × REAL (i)”. In this case, “r0” indicates the minimum interatomic distance, and “rd” indicates the update interval. “REAL (i)” indicates an arithmetic expression for converting the variable i into a real number. When the procedure of step S105 ends, the process proceeds to step S106. The procedure of step S105 is repeated until the interatomic distance r reaches the maximum interatomic distance.

  In step S106, it is determined whether or not the number N of atoms is “3”. When the number N of atoms is “3” (step S106; Yes), the process proceeds to step S108. In this case, it is necessary to obtain a three-body force acting between three atoms. On the other hand, when the number of atoms N is not “3” (step S106; No), that is, when the number of atoms N is “2”, the process proceeds to step S107. In this case, it is necessary to obtain a two-body force acting between two atoms.

In step S107, the molecular orbital calculation (MO) is used to calculate the interatomic interaction energy V ₂ (i) when the number of atoms N is “2” (hereinafter, the interatomic interaction energy is expressed as “ Potential energy ", or simply" potential "or" energy "). In this case, when the variable i is updated until the interatomic distance r reaches the maximum interatomic distance r0, the interatomic interaction energy V ₂ (i) corresponding to the variable i is repeatedly calculated. To do. When the procedure of step S107 ends, the process proceeds to step S117.

  On the other hand, in step S108, the variable j is set to “0”, and the process proceeds to step S109. In step S109, the third interatomic distance s defined in the atomic group in which the number N of atoms is “3” or more is obtained by an arithmetic expression “s = s0 + sd × REAL (j)”. In this case, “s0” represents the minimum third interatomic distance, “sd” represents the interval for updating, and “REAL (j)” represents an arithmetic expression for converting the variable j into a real number. When the procedure of step S109 ends, the process proceeds to step S110. Note that the procedure of step S109 is repeated until the third interatomic distance r reaches the maximum third interatomic distance.

  In step S110, the variable k is set to “0”, and the process proceeds to step S111. In step S111, an angle t sandwiching the central atom (hereinafter referred to as “interatomic angle”) defined in the atomic group having the number of atoms N of “3” or more is referred to as an arithmetic expression “t = t0 + td × REAL (k ) ”. In this case, “t0” indicates the minimum interatomic angle, “td” indicates the interval for updating, and “REAL (k)” indicates an arithmetic expression for converting the variable k to a real number. When the procedure of step S111 ends, the process proceeds to step S112. Note that the procedure of step S111 is repeated until the interatomic angle t reaches the maximum interatomic angle.

In step S112, intermolecular interaction energy V ₃ (i, j, k) when the number of atoms N is “3” is calculated using molecular orbital calculation (MO). In addition, until the interatomic distance r becomes the maximum interatomic distance from the minimum interatomic distance r0, until the third interatomic distance s becomes the maximum interatomic distance from the minimum third interatomic distance s0, and the interatomic angle When the variables i, j, k are updated to various values until t reaches the maximum interatomic angle t0, the interatomic interaction energy V ₃ corresponding to these variables i, j, k. (I, j, k) is repeatedly calculated. When the procedure of step S112 is completed, the process proceeds to step S113.

  In step S113, the variable k is updated to “k + 1”, and the process proceeds to step S114. In step S114, it is determined whether or not the variable k has been updated. Specifically, it is determined whether or not the interatomic angle t has reached the maximum interatomic angle. When the interatomic angle t is the maximum interatomic angle (step S114; Yes), the process proceeds to step S115. On the other hand, when the interatomic angle t is not the maximum interatomic angle (step S114; No), the process returns to step S111.

  In step S115, the variable j is updated to “j + 1”, and the process proceeds to step S116. In step S116, it is determined whether or not the update of the variable j has been completed. Specifically, it is determined whether the third interatomic distance s has reached the maximum third interatomic distance. If the third interatomic distance s is the maximum third interatomic distance (step S116; Yes), the process proceeds to step S117. On the other hand, if the third interatomic distance s is not the maximum third interatomic distance (step S116; No), the process returns to step S109.

  In step S117, the variable i is updated to “i + 1”, and the process proceeds to step S118. In step S118, it is determined whether or not the update of the variable i has been completed. Specifically, it is determined whether or not the interatomic distance r has reached the maximum interatomic distance. When the interatomic distance r is the maximum interatomic distance (step S118; Yes), the process proceeds to step S123. On the other hand, when the interatomic distance r is not the maximum interatomic distance (step S118; No), the process returns to step S105.

In step S123, the obtained potential map is stored. Specifically, when the number of atoms N is “2”, the interatomic interaction energy V ₂ (i) obtained by substituting various variables i, or the number of atoms N is “3” The interatomic interaction energy V ₃ (i, j, k) obtained by substituting various variables i, j, and k into the potential map. When the above procedure is completed, the flow is exited.

  On the other hand, when the number N of atoms is larger than “3”, that is, when the potential map is created by a crystal, steps S120 to S122 are performed. In step S120, a procedure for updating the interatomic distance r, the third interatomic distance s, the interatomic angle t, and the like is performed as in steps S105, S109, and S111 described above. Then, the process proceeds to step S121.

In step S121, any of the electronic state calculations described above is selected, and the interatomic interaction energy V _N (i, j, k,...) In the number N of atoms is calculated. More specifically, the interatomic interaction energy V _N (i, j, k, ...) is calculated. For example, the maximum number of atoms contained in the cluster is about 5-6. The interatomic interaction energy V _N (i, j, k,...) Can be calculated by electronic state calculation such as band calculation and molecular orbital calculation described above. When the procedure of step S121 ends, the process proceeds to step S122.

  In step S122, it is determined whether or not the variables i, j, k,... Have been updated until parameters such as the interatomic distance r and the third interatomic distance s reach the maximum value. If the update has been completed (step S122; Yes), the process proceeds to step S123. If the update has not been completed (step S122; No), the process returns to step S120.

In step S123, the obtained potential map is stored. Specifically, interatomic interaction energy V _N (i, j, k,...) Obtained by substituting various variables i, j, k,... Is stored as a potential map. When the above procedure is completed, the flow is exited.

  An example of the created potential map is shown in FIG. FIG. 4 shows the potential energy field on the surface of the SiC seed crystal. In addition, “V” in FIG. 4 indicates potential energy. Specifically, FIG. 4A shows a 3C—SiC seed crystal, and FIG. 4B shows a 2H—SiC seed crystal. From this, it can be seen that an expected stable point (grid position) appears on the potential map.

  Thus, the potential map according to the present embodiment is created by appropriately using electronic state calculation while taking into account the effects of atoms located far away on average. Therefore, the potential determined by this potential map is close to the potential obtained in the crystal field to be grown, and can be treated as a classical potential even if the electronic state changes due to recombination of atoms due to reaction or the like. it can.

[Phenomenon analysis method]
Next, a phenomenon analysis method using the above-described potential map will be described with reference to the flowchart shown in FIG. Here, as an example, a case where a single crystal growth of SiC by a sublimation method is simulated will be described. The phenomenon analysis method is executed by a personal computer or the like. Specifically, the personal computer performs a phenomenon analysis by executing a predetermined program.

  First, in step S201, conditions for the phenomenon to be analyzed (simulation conditions) are set. Specifically, the type of gas species, temperature, density (partial pressure), supply conditions, seed crystal size, orientation, temperature conditions, and the like are set. When the setting in step S201 ends, the process proceeds to step S202.

  In step S202, an initial state such as the position and velocity of the atom to be analyzed is set. For example, Si atoms and C atoms are arranged in a lower portion in a simulation box (model unit set in the simulation) in consideration of the crystal system and the surface orientation of the surface. Then, source gas molecules having a density proportional to a predetermined partial pressure are arranged in the upper part of the simulation box. In this case, the speeds of the crystal phase atoms and the gas phase atoms are set so as to have an average speed corresponding to the temperature. In addition, the size of the simulation box, the temperature of the crystal phase, the temperature of the gas phase, the gas partial pressure, the element type, coordinates and speed for each numbered atom are set as the initial state. When the procedure of step S202 is completed, the process proceeds to step S203.

  In step S203, the number, element type, position, and velocity of atoms in the vicinity are recorded as table data (hereinafter referred to as “cluster table”) for all atoms that are the target of the phenomenon analysis. . In this cluster table, atomic arrangement data displayed by polar coordinates in which the central atom is the origin, the closest atom is positioned on the z axis, and the next closest atom is positioned in the xz plane is recorded. Further, the cluster table also records the cluster orientation represented by a transformation matrix that represents the relationship between the polar coordinate system and the coordinate system fixed in space by Euler angles. When the procedure in step S203 is completed, the process proceeds to step S204.

  In step S204, the potential energy is read with reference to the stored potential map. Specifically, first, the atomic arrangement data of the central atom (hereinafter referred to as “target atom”) that is the target of potential energy calculation is read from the cluster table, and the point in the potential map closest to this atomic arrangement data ( Hereinafter, a point closest to the closest point (hereinafter referred to as an “adjacent point”) is determined. Then, the potential energy at the nearest point and the potential energy at the adjacent point are read out. Such a procedure is performed because the target atom is not necessarily positioned on the discrete point because the potential map is defined by the discrete point. When the procedure of step S204 is completed, the process proceeds to step S205.

  In step S205, the following procedure is performed. First, the potential energy of the target atom is spline-interpolated between the potential energy of the closest point determined in step S204 and the potential energy of the adjacent point based on the distance between the target atom's atomic arrangement data and the closest point in the potential map. Ask by. Further, the force acting on the target atom is obtained by obtaining the difference between the potential energy at the nearest point and the potential energy at the adjacent point and performing spline interpolation on the value obtained by decomposing the obtained difference into each component. Thus, by obtaining the potential energy by spline interpolation, the accuracy of the phenomenon analysis can be ensured even if the number of potential map data is reduced. Note that the method of interpolating the potential map when obtaining the potential energy is not limited to spline interpolation. In order to improve the calculation speed, the potential energy at the nearest point may be used as it is without performing interpolation.

  Next, a change in the position and velocity of the target atom is obtained by substituting the force acting on the target atom obtained as described above into an expression obtained by differentiating the motion equation (for example, Berlet's motion equation). Then, the position and speed of the target atom are updated by adding the obtained change in the position and speed of the target atom to the value obtained in the procedure of step S205. The calculation of energy, force, position, and velocity described above is performed for all atoms included in the system.

  Further, in step S205, the temperature and pressure of the gas phase are obtained by calculating the average value and virial of the kinetic energy for all atoms present in the gas phase. At the same time, the temperature of the crystal phase is obtained by calculating the average value of the kinetic energy for all atoms present in the crystal phase. When the procedure in step S205 is completed, the procedure in step S206 is performed.

  In the procedure of step S206, it is determined whether or not the number of steps in which the procedure of steps S203 to S205 has been repeated has reached a predetermined cycle (correction cycle). For example, the correction cycle corresponds to the number of times that the thermodynamic data such as temperature and pressure is repeatedly executed before correction. If the number of steps has reached the correction cycle (step S206; Yes), the process proceeds to step S207. If the number of steps has not reached the correction cycle (step S206; No), the process proceeds to step S203. Return.

  In step S207, the average value of temperature and pressure is updated, and the temperature and pressure are corrected. Specifically, the temperature and pressure of the system are updated in a predetermined period, and the average value of the temperature and pressure is calculated during this update. Before updating, the system temperature and pressure are compared with the set conditions set in steps S201 and S202, and if there is a difference, correction is performed by an appropriate method. When the procedure of step S207 is completed, the process proceeds to step S208.

  In the procedure of step S208, it is determined whether or not the number of steps obtained by repeating the procedure from steps S203 to S207 has reached a predetermined step (final step). This final step corresponds to the time (that is, the number of repeated executions) for which simulation should be performed. If the number of steps has reached the final step (step S208; Yes), the procedure proceeds to step S209. If the number of steps has not reached the final step (step S208; No), the procedure returns to step S203. .

  In step S209, a phenomenon analysis file obtained by performing the above simulation is created. Specifically, a phenomenon analysis file such as a time-dependent change of the snapshot in the simulation box is created based on data such as the position / velocity of atoms recorded at any time during the execution of the simulation. When the above procedure is completed, the flow is exited.

  Thus, in the phenomenon analysis method according to the present embodiment, the phenomenon analysis is performed using a potential map created in advance. Accordingly, it is not necessary to calculate the potential energy of the atom by the Monte Carlo method (MC) or the molecular dynamics method every time the state of the system is changed, so that the calculation time can be shortened. Furthermore, since the potential map is created by appropriately using the electronic state calculation while taking into account the effects of atoms located far away, the phenomenon is analyzed with high accuracy by using the potential map. be able to. As described above, according to the phenomenon analysis method according to the present embodiment, it is possible to perform phenomenon analysis with high accuracy and high speed for a large-scale system.

It is a figure for demonstrating the basic concept of the electronic state calculation model used for potential map preparation. It is a figure for demonstrating the specific example of electronic state calculation. It is a flowchart which shows the preparation method of the potential map which concerns on this embodiment. It is a figure which shows the specific example of a potential map. It is a flowchart which shows the phenomenon analysis method which concerns on this embodiment.

Explanation of symbols

1 First atomic group 2 Second atomic group 3 Third atomic group 11 Cluster 13 Supercell 15 Crystal lattice

Claims (4)

A first atomic group composed of a predetermined number of atoms, a second atomic group located around a region where the first atomic group exists, a region where the first atomic group and the second atomic group exist A potential map creation procedure for creating a potential map showing the relationship between the position of an atom and energy or force by performing electronic state calculation on a model composed of a third group of atoms located around;
A molecular simulation method comprising: a phenomenon analysis procedure for performing a phenomenon analysis with reference to the created potential map. The first atomic group and the third atomic group are fixed;
The molecular simulation method according to claim 1, wherein the second atomic group is moved under the influence of the first atomic group and the third atomic group.   The molecular simulation method according to claim 1, wherein the electronic state calculation is performed by band calculation or molecular orbital calculation. The molecular simulation method according to any one of claims 1 to 3, wherein the phenomenon analysis procedure includes a procedure for obtaining energy or force of an atom by performing spline interpolation on the potential map.

Images