WO2012056562A1 - 原子構造等計算装置，方法，およびプログラム - Google Patents

原子構造等計算装置，方法，およびプログラム Download PDF

Info

Publication number
WO2012056562A1
WO2012056562A1 PCT/JP2010/069304 JP2010069304W WO2012056562A1 WO 2012056562 A1 WO2012056562 A1 WO 2012056562A1 JP 2010069304 W JP2010069304 W JP 2010069304W WO 2012056562 A1 WO2012056562 A1 WO 2012056562A1
Authority
WO
WIPO (PCT)
Prior art keywords
calculation
region
atomic
atom
substance
Prior art date
Application number
PCT/JP2010/069304
Other languages
English (en)
French (fr)
Inventor
健太郎 高井
Original Assignee
富士通株式会社
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by 富士通株式会社 filed Critical 富士通株式会社
Priority to PCT/JP2010/069304 priority Critical patent/WO2012056562A1/ja
Publication of WO2012056562A1 publication Critical patent/WO2012056562A1/ja

Links

Images

Classifications

    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16C—COMPUTATIONAL CHEMISTRY; CHEMOINFORMATICS; COMPUTATIONAL MATERIALS SCIENCE
    • G16C10/00—Computational theoretical chemistry, i.e. ICT specially adapted for theoretical aspects of quantum chemistry, molecular mechanics, molecular dynamics or the like

Definitions

  • the present invention relates to a calculation processing technique of an atomic structure and its electronic state (nucleus position and electron distribution) in a material design support application such as designing a stable structure of metal, semiconductor, living body and the like.
  • the one-electron approximation is an approximation of the state of a multi-electron system by regarding the multi-electron system as a collection of electrons without interaction in an effective potential and occupying the electrons in order from the one-electron state with the lowest energy. It is a method to express.
  • the states of the multi-electron system shown in the following equations (2) and (3) can be expressed by the Schrödinger equation in which the electronic states are integrated.
  • ⁇ eff (r) is the effective potential
  • m is the mass of the electron
  • ⁇ n is the eigenvalue of the electron in the nth state when arranged in order of decreasing energy
  • ⁇ n (r) is the order of decreasing energy. It is a wave function of electrons in the nth state when arranged.
  • the electronic state can be calculated by using the one-electron approximation. If one-electron approximation is used, the energy of the whole system is expressed by the following equation (4).
  • the first item is an electron energy
  • the second item E rep is a term for correcting the interaction energy between the nucleus and the nucleus and the electron energy of the first item.
  • F ( ⁇ n ) represents a Fermi distribution function.
  • the second item E rep there is the following format (C. Molteni, L. Colombo, and L. Miglio, “Condens. Matter (6)”, 1994, 5243).
  • r ij is a distance between atoms of the site i and the site j
  • U 1 , U 2 , r 0 , and r D are parameters that express the characteristics unique to the atom.
  • the force acting on the atom can be calculated by differentiating E tot which is the energy of the entire system with respect to the position RA of the atom, as shown in the following equation (5).
  • molecular dynamics calculation can be performed by moving atoms according to the equation of motion from the forces acting on the obtained atoms.
  • the calculation method for obtaining the eigenvalue and wave function of electrons has the disadvantage that the calculation amount is proportional to the cube of the number of atoms, so that it takes too much time to calculate a large system. Therefore, it was difficult to adopt this method when a reduction in calculation time was required, such as a simulation comparing with experimental data.
  • N the number of atoms to be simulated
  • ⁇ i ⁇ is a localized basis
  • i is an atomic site
  • ⁇ is an atomic orbital of a certain atomic site
  • C n i ⁇ is a coefficient of ⁇ i ⁇ .
  • equation (2) when the wave function is expanded on a certain local basis, the problem of solving the Schrödinger equation (differential equation) in equation (2) can be replaced with the problem of obtaining the coefficient c n i ⁇ of the local basis (eigenvalue problem). it can. Specifically, ⁇ j ⁇ * (r ⁇ j ) is multiplied from the left on both sides of equation (2), and after integrating equation (6), integration is performed over the entire space. That is, the following equation (8) is obtained.
  • H j ⁇ , i ⁇ , S j ⁇ , i ⁇ in the equation (8) are defined as the following equation (9).
  • the component relating to the electron energy can be represented as the following equation (14).
  • ⁇ i ⁇ , j ⁇ represents a matrix element of the energy density matrix.
  • ⁇ and ⁇ are quantities that do not depend on the number of atoms.
  • the density operator and the energy density operator localization i.e. the matrix element of the localized basis phi i.alpha and consisting phi Jbeta density matrix apart over a certain distance " ⁇ i ⁇ , j ⁇ and ⁇ i ⁇ , j ⁇ is 0
  • ⁇ i ⁇ , j ⁇ and ⁇ i ⁇ , j ⁇ is 0
  • the number of sums of j is constant regardless of the number of atoms.
  • the amount of calculation for obtaining energy and force is O (N). That is, structural optimization and molecular dynamics calculation can be performed using the energy and force obtained by the above O (N) method.
  • DC method divides the entire system into a plurality of regions, finds the density matrix of the divided regions, finds the physical quantities of the divided regions using the density matrix of the divided regions, and finally at each region. This is a method to obtain the physical quantities of the entire system by integrating the obtained physical quantities.
  • the molecule is divided into regions in real space, the total energy of the molecule is divided into contributions for each region and contributions between regions, and the contribution of a region with each divided energy term is changed to other molecules including that region.
  • a method has been proposed in which a desired approximation level is selected from the approximation levels prepared in advance and calculated with the contribution of the region.
  • FIG. 1 is a diagram showing an example of divided areas and buffer areas.
  • Fig. 1 is a diagram showing the atomic space of a molecular system in a pseudo plane.
  • the divided area 1 is an area divided for application of the DC method.
  • the buffer area 2 is an area adjacent to the divided area 1.
  • the area obtained by combining the divided area 1 and the buffer area 2 is an area (calculation area) required for obtaining the density matrix.
  • I in the equation (16) is an index of the divided area.
  • ⁇ ij (I) is defined as the following equations (17) and (18).
  • localized basis phi i.alpha and phi Jbeta is, when both are in the divided region 1 shown in FIG. 1, multiplied by a factor 1, if only on the divided area 1 one, multiplied by a factor 0.5, both If it is not in the divided area 1, the density matrix is set to 0.
  • the calculation area composed of the divided area 1 and the buffer area 2 shown in FIG. 1 is extracted, and the localized basis included in the extracted calculation area is ( 9)
  • the matrix component is obtained as shown in the equation.
  • the density matrix component of the divided area 1 is expressed by the following equation (21).
  • the Fermi energy that appears in the Fermi distribution function must be obtained.
  • the Fermi energy can be obtained from the condition shown by the following equation (22) that “the number N of electrons in the entire system is constant”.
  • ⁇ F Fermi energy
  • ⁇ i ⁇ , j ⁇ (I) is the following equation (26).
  • the density matrix is obtained using the local basis included in the calculation region, and the number of electrons in each divided region 1 is calculated from the density matrix to determine the Fermi energy. Thereafter, the physical quantity of each divided region 1 is calculated, and the physical quantity of the entire system can be calculated by obtaining the sum.
  • Matrix elements ⁇ i ⁇ , j ⁇ localized basis phi i.alpha and consisting phi Jbeta density matrix, when localized basis is in the vicinity has a value. Therefore, in order to obtain the matrix element of the density matrix at the end of the divided area 1 (immediately inside the area), the buffer area 2 must be prepared.
  • the size of the buffer area 2 must be set so that the physical quantity to be obtained can be obtained with sufficient accuracy by the calculation system.
  • the size of the buffer area 2 has been manually calculated based on an empirical rule.
  • the number of divisions increases in proportion to the number of atoms when the number of atoms is increased.
  • the coefficient a is proportional to the cube of the number of atoms n of “divided region 1 + buffer region 2”.
  • the coefficient a increases when the buffer area 2 is wide, a calculation method capable of reducing the calculation time by making the buffer area 2 as narrow as possible is desired. In particular, in order to perform simulation while actually comparing with experimental data, it is necessary to further reduce the calculation time.
  • An object of the present invention is to provide a technique capable of reducing the amount of calculation when calculating an atomic structure and an electronic state by a calculation method in which the calculation amount is proportional to the number of atoms and dividing the atomic space into regions. That is.
  • An atomic structure calculation device disclosed as one aspect of the present invention includes a calculation condition acquisition unit that acquires initial structure information of a substance, a region division unit that divides the atomic space of the substance into regions, and virtual atom parameters. For each of the divided areas of the substance, as information indicating the calculation area including the divided area, parameters and end positions of the virtual atoms that end with respect to the area, and the divided areas and the virtual atoms are terminated.
  • the atomic structure or electronic state of the substance by a calculation region determination unit that determines an atomic layer that becomes a buffer region with the atomic layer to be processed, and a calculation process that applies a divide and conquer method to a calculation process in which the calculation amount is proportional to the number of atoms
  • a calculation processing unit that performs calculation using the calculation area.
  • the above atomic structure calculation apparatus by terminating the atoms at the end of the divided region into virtual atoms, the surface electronic state which does not exist at the end of the region can be eliminated. Can be calculated using a narrow buffer area. Therefore, it is possible to reduce the calculation amount of the calculation of the atomic structure and electronic state.
  • FIG. 1 is a diagram illustrating a hardware configuration example of an atomic structure calculation device 10.
  • the atomic structure calculation device 10 is a device that calculates the atomic structure or electronic state of a target substance by dividing the atomic space into regions as shown in FIG. 1 in a calculation method in which the calculation amount is proportional to the number of atoms. is there.
  • FIG. 2 is a diagram illustrating a configuration example in one embodiment of the atomic structure calculation apparatus 10.
  • the atomic structure calculation device 10 includes a calculation condition acquisition unit 11, a region division unit 12, a calculation region determination unit 13, a calculation processing unit 14, and a virtual atom information storage unit 15.
  • the calculation condition acquisition unit 11 acquires calculation conditions and initial atomic structure information of a target substance.
  • Calculation conditions include, for example, simulation temperature, simulation pressure, designation of potential of atomic species to be simulated, number of time steps of molecular dynamics, and calculation accuracy by O (N) method.
  • the calculation accuracy is set by ⁇ E when judging by energy E, and by ⁇ F when judging by force F acting on the atoms.
  • the initial atomic structure information includes, for example, a unit cell vector, lattice constant, number of atoms, atomic species, atomic coordinates, initial atomic velocity, and the like.
  • the region dividing unit 12 divides the entire atomic space of the target substance into rectangular parallelepipeds in order to execute the DC method.
  • the calculation area determination unit 13 determines a calculation area including the divided area 1 and the buffer area 2 adjacent to each of the areas (divided areas) 1 divided by the area dividing unit 12. Specifically, the calculation region determination unit 13 uses the parameters of the virtual atom 3 that terminates at the atoms at the end of the region, the termination position, and the atomic layer of the buffer region 2 as information indicating the calculation region for each divided region 1. decide.
  • FIG. 3 is a diagram showing an example of termination of the virtual atom 3.
  • each atom terminates in a virtual atom 3 in the atomic layer at the end of the calculation region with respect to the divided region 1 in the atomic space of the target substance.
  • the virtual atom 3 is a case where it terminates at a hydrogen atom type atom having a valence electron number of one.
  • the calculation region determination unit 13 obtains the calculation result of the calculation processing unit 14, that is, the energy or the force acting on the atoms in the case of the calculation region determined for the divided region 1 satisfying a predetermined calculation accuracy.
  • the atomic layer of the buffer region 2 when the number of atomic layers included in the calculation region is minimized can be obtained.
  • the calculation area determination unit 13 determines information indicating the area in such a case (parameters of the terminating virtual atom 3, termination position, atomic layer of the buffer area 2) as simulation information.
  • the calculation processing unit 14 determines the atomic structure or the electronic state of the calculation region (the parameter of the virtual atom 3 to be terminated, the termination position, the atomic layer of the buffer region) for the divided region 1 determined by the calculation region determination unit 13. Calculate the energy or the force acting on the atom.
  • the calculation processing unit 14 performs the above-described calculation processing on the calculation regions having various sizes determined by the calculation region determination unit 13 for one divided region 1.
  • the virtual atom information storage unit 15 stores virtual atom information including the parameters of the virtual atom 3 and interatomic distances indicating stable bond distances between the virtual atom 3 and all atomic species.
  • a group 1 element such as hydrogen (H) and lithium (Li)
  • a group 17 group such as fluorine (F) and chlorine (Cl)
  • a parameter indicating the element type is stored.
  • FIG. 4 is a diagram showing an example of a simulation process flow using a calculation result by the atomic structure calculation device 10.
  • Step S1 In the atomic structure calculation device 10, the calculation condition acquisition unit 11 acquires the calculation conditions for the simulation and the initial atomic structure information of the target substance.
  • Step S2 The region dividing unit 12 converts the entire system (total atomic space) of the target substance into a rectangular parallelepiped as a region (divided region) for executing the calculation process in the O (N) method by applying the DC method.
  • Step S3 The calculation region determination unit 13 refers to the information in the virtual atom information storage unit 15, and for each divided region 1, as information indicating the calculation region, the parameters of the virtual atom 3 used for termination, the termination position, The atomic layer of the buffer region 2 is determined.
  • At least one parameter of the virtual atom 3 is prepared in the atomic structure calculation device 10.
  • a parameter of the virtual atom 3 for example, a hydrogen atom type parameter is used.
  • the hydrogen atom type parameter has a single valence electron and has a fully symmetric potential.
  • the parameters of the virtual atom 3 for example, a form of interaction with the partner atom s, p, d,.
  • the atomic structure calculation device 10 may include a virtual atom information storage unit 15.
  • the parameters of the virtual atom 3 are determined from the parameters indicating the element format stored in the virtual atom information storage unit 15. Is selected.
  • the terminal position is determined by the terminal direction and the distance between atoms.
  • the termination direction includes the position of the atom at the end of the region (the atom terminating in the virtual atom 3), and the position of the atom bonded to the atom before dividing the region. It is set in the direction on the straight line.
  • the distance between atoms is the distance at which the energy of the atom and the virtual atom 3 is stably bonded in the direction of the end of the atom at the end of the region.
  • CH bond length of methane (CH4) is 1.1 ⁇
  • Si—H bond length of silane (SiH 4) is 1.5 ⁇
  • Si silicon
  • the virtual atom 3 is hydrogen
  • Step S4 The calculation processing unit 14 calculates the atomic structure and electronic state by the O (N) method for the calculation region determined in step S3. More specifically, the calculation processing unit 14 performs calculation using the above-described equations (16) to (26) as an example.
  • step S3 each time the calculation area determination unit 13 changes the number of atomic layers in the buffer area 2 for one divided area 1, the calculation process part 14 is determined in step S3. Repeat the calculation area processing. Then, the calculation area determination unit 13 is a calculation area in the case where the result of calculation processing of the atomic structure and the electronic state is obtained satisfying the given calculation accuracy and the size of the buffer area 2 is minimized. (Parameters of the virtual atom 3, terminal position, atomic layer of the buffer region 2) are adopted as information for the calculation region for simulation, and the calculation result in this case is passed to the simulation device 20.
  • Step S5 The simulation device 20 performs a simulation using the calculation result obtained from the atomic structure calculation device 10, and outputs the simulation result.
  • FIG. 6 is a diagram showing a more detailed processing flow example of the processing in step S3.
  • step S3 the information indicating the calculation area for each divided area 1, the parameter of the terminating virtual atom 3, the terminal position, and the atomic layer of the buffer area 2 are determined by the following process.
  • Expression (27) is a table expression in which the jump integration of Expression (19) is divided for each type of atomic orbital bond (ss ⁇ , sp ⁇ , pp ⁇ , pp ⁇ , etc.), and Expression (28) is Expression (4). This represents the interaction between a pair of ions in the ion-ion interaction.
  • ⁇ represents the type of atomic orbital bond
  • r represents the interatomic distance.
  • r 0 , h ⁇ 0 , r D , U 1 , U 2 are parameters.
  • Equation (29) is an expression obtained by dividing the jump integration of the equation (19) for each type of atomic orbital coupling (ss ⁇ , sp ⁇ , pp ⁇ , pp ⁇ , etc.).
  • Equation (30) expresses the interaction between a set of ions and ions among the ion-ion interactions of Equation (4).
  • ⁇ represents the type of atomic orbital bond
  • r represents the interatomic distance.
  • the terminal position is the direction of the straight line between the virtual atom 3 and the atom that terminates, and the position of the atom that was bonded before the region splitting. Distance).
  • Step S12 The calculation area determination unit 13 is for calculation when the size of the buffer area 2 is given in advance (for example, when the atomic layer to be the buffer area 2 is sufficiently enlarged by adding the atomic layer). Calculate the energy or force acting on an atom for a region.
  • the calculation region determination unit 13 further increases one atomic layer to the atomic layer determined in step S10, and calculates the energy or force acting on the atom while keeping the parameters and termination distance of the virtual atom 3 the same. Is the energy or force acting on the atoms when the buffer region 2 is sufficiently large.
  • Step S13 The calculation region determination unit 13 calculates energy or force acting on the atom by changing the parameter of the virtual atom 3 and the terminal distance. Further, the calculation region determination unit 13 sets A ′ as the minimum absolute value of the difference between the calculated energy or force acting on the atom and the energy or force acting on the atom when the atomic layer is sufficiently large. calculate.
  • the calculation region determination unit 13 changes the parameters of the virtual atom 3 of the hydrogen atom type, r 0 , h ⁇ 0 , r D , U 1 , U 2 , and the distance between the atom at the end of the region and the virtual atom 3,
  • the absolute value of the difference in energy or force acting on the atoms when the buffer region is sufficiently wide is calculated and compared with the calculation accuracy ⁇ E and ⁇ F.
  • examples of ranges for changing parameters are as follows.
  • Parameter r 0 The range is from (hydrogen atom type parameter ⁇ 0.1 cm) to (hydrogen atom type parameter +0.1 cm), and the step size is 0.01 cm.
  • Parameter h ⁇ 0 The range is from (hydrogen atom type parameter ⁇ 1 eV) to (hydrogen atom type parameter +1 eV), and the step size is 0.1 eV.
  • Parameter r D (Hydrogen atom type parameter -0.1 cm) to (Hydrogen atom type parameter +0.1 cm), with a step size of 0.01 cm.
  • Parameter U 1 The range is from (hydrogen atom type parameter ⁇ 1 eV) to (hydrogen atom type parameter +1 eV), and the step size is 0.1 eV.
  • Parameter U 2 (Hydrogen atom type parameter ⁇ 1 eV) to (Hydrogen atom type parameter + 1 eV), and the step size is 0.1 eV.
  • Step S14 The calculation area determination unit 13 determines whether the absolute minimum value A ′ has a predetermined calculation accuracy and is smaller than ⁇ E when judged by energy or ⁇ F when judged by the force acting on the atoms. To do.
  • step S14 If the determined minimum value A ′ is smaller than the calculation accuracy ⁇ E or ⁇ F (Y in step S14), the process proceeds to step S15. If the minimum value A ′ is not smaller than the calculation accuracy ⁇ E or ⁇ F (step S14) N of S14), the process proceeds to step S19.
  • Step S15 The calculation region determination unit 13 reduces the atomic layer in the buffer region by one layer, and calculates energy or force acting on the atoms with the same parameters and termination positions of the virtual atom 3 as in the process of step S13. Further, the calculation area determination unit 13 calculates the minimum value A ′′ of the absolute value of the difference between the calculated energy or force acting on the atom and the energy or force acting on the atom when the atomic layer is sufficiently large. To do.
  • Step S16 It is determined whether the minimum value A ′′ is smaller than the calculation accuracy ⁇ E or ⁇ F.
  • step S16 If the minimum value A ′′ is smaller than the calculation accuracy ⁇ E or ⁇ F (Y in step S16), the process proceeds to step S17. If the minimum value A ′′ is not smaller than the calculation accuracy ⁇ E or ⁇ F (N in step S16). ), Go to step S18.
  • Step S18 The calculation area determination unit 13 adopts the parameters of the virtual atom 3 in the process of step S13, the end position, and the atomic layer of the buffer area, and ends the process.
  • Step S19 The calculation area determination unit 13 checks whether the processing counter i is greater than 0 (i> 0). If the process counter i is not greater than 0 (N in step S19), the process proceeds to step S110. If the process counter i is greater than 0 (Y in step S19), the process proceeds to step S111.
  • Step S110 The calculation region determination unit 13 increases the atomic layer by one, and returns to the process of step S13.
  • Step S111 The calculation region determination unit 13 adopts the parameters of the virtual atom 3, the terminal position, and the atomic layer of the buffer region 2 immediately before entering the processing of Step S13, and ends the processing.
  • the parameters of the virtual atom 3, the end distance, and the atomic layer which are processing accuracy satisfying the given processing accuracy ⁇ E and ⁇ F and indicate the calculation area when the atomic layer of the buffer region 2 is reduced, are obtained. Obtainable.
  • FIG. 7 is a diagram showing an example of a processing flow for determining an atomic layer in the buffer region 2 when a hydrogen atom parameter is used as the virtual atom 3.
  • Step S20 The calculation region determination unit 13 determines the terminal position (bonding direction and terminal distance) using the parameter of the hydrogen atom as the virtual atom 3, and the atom in the atomic layer i at the end of the buffer region 2 with respect to the divided region 1 Terminate.
  • Step S22 The calculation processing unit 14 obtains energy and force (B i ) acting on the atoms for the calculation region determined in step S20.
  • Step S23 The calculation area determination unit 13 adds one atomic layer.
  • Step S25 Determine the bond direction and end distance in the same way as before adding the atomic layer.
  • Step S26 The calculation processing unit 14 obtains energy and force acting on the atoms (B i + 1 ) for the calculation region determined in step S25.
  • Step S27 The calculation area determination unit 13 calculates the absolute value of the difference between B i + 1 and B i, and the calculated absolute value is based on the condition (calculation accuracy) given to the energy or the force acting on the atom. Determine if it is small.
  • Step S28 If the absolute value obtained in step S27 satisfies the condition (Y in step S28), the calculation area determination unit 13 proceeds to the process of step S29, and if the condition is not satisfied (in step S28) N), the process returns to step S24.
  • Step S29 The calculation area determination unit 13 adopts i as the number of atomic layers.
  • FIG. 8 is a diagram illustrating a hardware configuration example of the atomic structure calculation device 10.
  • the atomic structure calculation device 10 can be implemented by a computer 100 including a CPU 101, a temporary storage device (DRAM, Flash Memory, etc.) 102, and a persistent storage device (HDD, Flash Memory, etc.) 103.
  • a computer 100 including a CPU 101, a temporary storage device (DRAM, Flash Memory, etc.) 102, and a persistent storage device (HDD, Flash Memory, etc.) 103.
  • the atomic structure calculation device 10 can be implemented by a program executable by the computer 100.
  • a program describing the processing content of the functions that the atomic structure calculation apparatus 10 should have is provided.
  • the computer 100 executes the provided program, the processing functions of the atomic structure calculation device 10 described above are realized on the computer 100.
  • the calculation condition acquisition unit 11, the region division unit 12, the calculation region determination unit 13, and the calculation processing unit 14 of the atomic structure calculation device 10 can be configured by a program, and the virtual atom information storage unit 15 is permanently stored. It can be constituted by the sexual storage device 103.
  • the computer 100 can also read a program directly from a portable recording medium and execute processing according to the program.
  • the computer 100 includes a network interface and can sequentially execute processing according to the received program every time the program is transferred from the server computer.
  • this program can be recorded on a recording medium readable by the computer 100.
  • the atomic structure calculation apparatus 10 terminates the atoms at the end of the calculation region in the virtual atom 3 in the atomic structure / electronic state calculation method based on the quantum theory using the DC method.
  • the buffer area 2 can be reduced to reduce the amount of calculation.
  • the atomic structure calculation device 10 can perform high-speed calculation of the atomic structure and electronic state, the number of atoms that can be handled can be expanded.
  • the number of atoms in the divided region and the buffer region conventionally needs to be 216 atoms. According to this, only 162 atoms are required, and the calculation speed is improved by about 2.4 times.
  • regions can be performed independently for every division area. Therefore, it is easy to speed up each divided region by parallel calculation.
  • the disclosed high-speed calculation method using the DC method executed by the atomic structure calculation device 10 can be applied to either first-principles calculations or semi-empirical quantum mechanical calculations.
  • the disclosed calculation apparatus 10 for atomic structure can improve the development efficiency by comparing the experimental results with the simulation results in the development of new materials and devices. Specifically, it takes less time to obtain a guideline for determining the temperature and pressure to obtain the target material and the composition ratio of the material from the simulation results, thereby improving the development speed. be able to.
  • the disclosed atomic structure calculation device 10 contributes to the reduction of development time and cost, and the reduction of the environmental load due to development.

Landscapes

  • Engineering & Computer Science (AREA)
  • Computing Systems (AREA)
  • Theoretical Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Health & Medical Sciences (AREA)
  • General Health & Medical Sciences (AREA)
  • Spectroscopy & Molecular Physics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

　開示する原子構造等計算装置１０は，物質の初期構造情報を取得する計算条件取得部１１と，物質の原子空間を領域に分割する領域分割部１２と，仮想原子のパラメータを用意し，分割領域各々に対して，領域の端の原子に仮想原子を終端させて，分割領域と仮想原子が終端する原子層とのバッファ領域でなる計算用領域を決定する計算用領域決定部１３と，計算量が原子数に比例する計算処理に分割統治法を適用した計算処理によって，物質の原子構造または電子状態を求める計算処理部１４とを備え，物質の原子構造または電子状態の計算量を軽減する。

Description

原子構造等計算装置，方法，およびプログラム
　本発明は，金属，半導体，生体等の安定構造の設計等の材料設計支援アプリケーションにおける，原子構造とその電子状態（原子核の位置と電子分布）の計算処理技術に関する。
　トランジスタ，太陽電池等の新素材の設計やナノデバイスの開発を行う際に，量子論にもとづくシミュレーションを用いる傾向が年々高まっている。デバイスの微細化が進む現在において，目的にあった材料の開発には，ナノスケールでの物理現象を理解することが必須だからである。新材料のナノスケールでの開発には，構造最適化シミュレーション，分子動力学シミュレーションを，実験データと比較しながら行う必要がある。そのためには，短時間でシミュレーションできるアプリケーションが必要となる。
　量子論にもとづいて物質の原子構造・電子状態などのシミュレーションを行うためには，以下の（１）式に示す，ミクロなスケールでの支配方程式であるシュレディンガー方程式を解く必要がある。
Figure JPOXMLDOC01-appb-M000001
　現実の物質には多くの電子が存在している。多くの電子が含まれている系（以下，多電子系という）を計算する１つの手法として，一電子近似がある。一電子近似とは，多電子系をある有効ポテンシャル中における相互作用のない電子の集まりとみなして，電子をエネルギーの低い一電子状態から順番に占有させることによって，多電子系の状態を近似的に表す方法である。一電子近似を用いることにより，以下の（２），（３）式に示す，多電子系の状態を，電子状態は一体であるシュレディンガー方程式によって表すことができる。
Figure JPOXMLDOC01-appb-M000002
　ここで，υeff(r)は有効ポテンシャル，ｍは電子の質量，εnはエネルギーが低い順番に並べたときのｎ番目の状態の電子の固有値，ψn(r)はエネルギーが低い順番に並べたときのｎ番目の状態の電子の波動関数である。このように，一電子近似を用いることによって，電子状態を計算できる。
一電子近似を用いると全系のエネルギーは，以下の（４）式となる。
Figure JPOXMLDOC01-appb-M000003
　ここで，第１項目は電子エネルギー，第２項目Ｅrepは原子核と原子核の相互作用エネルギーと第１項目の電子エネルギーを補正する項である。また，ｆ(εn)はフェルミ分布関数を表す。例えば，第２項目Ｅrepの例として，以下の形式がある（C. Molteni, L. Colombo, and L. Miglio, "Condens. Matter (6)", 1994, 5243）。
Figure JPOXMLDOC01-appb-M000004
　ここで，ｒijは，サイトｉとサイトｊの原子間の距離，Ｕ1，Ｕ２，ｒ０，ｒDは，原子固有の性質を表現するパラメータである。
　原子に働く力は，以下の（５）式に示すように，全系のエネルギーであるＥtotを原子の位置ＲAで微分して計算することができる。
Figure JPOXMLDOC01-appb-M000005
　このようにして得た原子に働く力から，原子に働く力が，０，または，実際のシミュレーションにおいてパラメータとして入力した値以下になるように原子を動かすことによって，安定な原子配置を計算することができる（構造最適化計算）。
　また，求めた原子に働く力から，運動方程式に従って原子を動かすことによって，分子動力学計算を行うことができる。
　上述のシュレディンガー方程式について，電子の固有値と波動関数を求める計算方法は，計算量が原子数の３乗に比例するため，大きな系を対象とする場合には計算時間がかかりすぎるという欠点がある。したがって，実験データと比較しながらのシミュレーションのように，計算時間の短縮が求められている場合に採用することは難しかった。
　そのため，必要な物理量だけを短時間で求める方法として，計算量が原子数に比例するオーダーN法（Ｏ（Ｎ）法，Nはシミュレーションする原子数）がある。
　量子論を用いて原子構造・電子状態を求めるＯ（Ｎ）法では，一電子の固有値と波動関数を求めずに，必要な物理量だけを求める。全系のエネルギーや原子に働く力をＯ（Ｎ）法を用いて求める場合には，密度行列を導入して計算を行う。
　以下に，密度行列を用いたＯ（Ｎ）法を説明する。
　電子の波動関数ψnを局在基底で展開すると以下の（６）式のように表せる。
Figure JPOXMLDOC01-appb-M000006
　ここで，φiαは局在基底，ｉは原子のサイト，αはある原子サイトの原子軌道を表す。また，ｃｎ iαはφiαの係数である。（６）式を，ディラックが発明したブラ-ケットの記法を用いると下記の（７）式のように表すことができる。
Figure JPOXMLDOC01-appb-M000007
　このように，波動関数をある局在基底で展開すると，（２）式のシュレディンガー方程式（微分方程式）を解く問題は，局在基底の係数ｃｎ iαを求める問題（固有値問題）に置き換えることができる。具体的には，（２）式の両辺に左からφjβ *(ｒ－τj)を掛け，（６）式を代入した後に全空間で積分する。すなわち，以下の（８）式のようになる。
Figure JPOXMLDOC01-appb-M000008
　ここで，（８）式のＨjβ，iα，Ｓjβ，iαを，以下の（９）式のように定義する。
　（９）式の右辺はブラ-ケットの記法を用いて表したものである。
Figure JPOXMLDOC01-appb-M000009
　ここで，密度演算子ρ＾（式では，ρの上部に＾）を導入する。密度演算子は以下のように書き表すことができる。
Figure JPOXMLDOC01-appb-M000010
　したがって，（１０）式に（７）式を代入すると，以下の式（１１）となり，
ここで，ρiα，jβは，以下の式（１２）のように表すことができる。
Figure JPOXMLDOC01-appb-M000011
　すると，電子エネルギーを表す（４）式の右辺第１項は，（７）式を代入して（１２）式を用いることにより，以下の式（１３）のように，密度行列の行列要素を用いて表すことができる。
Figure JPOXMLDOC01-appb-M000012
　また，（５）式で表した原子に働く力の項の中で，電子エネルギーに関する成分は，以下の式（１４）のように表すことができる。
Figure JPOXMLDOC01-appb-M000013
　ここで，式（１５）のように，ωiα，jβは，エネルギー密度行列の行列要素を表す。
Figure JPOXMLDOC01-appb-M000014
　（１３）式と（１４）式の４つの和のインデックスのうち，α，βは，原子数に依存しない量である。また，密度演算子とエネルギー密度演算子の局在性，すなわち「ある距離以上に離れた局在基底φiαとφjβからなる密度行列の行列要素ρiα，jβとωiα，jβは，０とみなせる」（W. Kohn, "Int. J. Quart. Chem. 56", 1995, 299）という物理的性質より，ｊの和の数は原子数によらず一定となる。
　したがって，エネルギーや力を求めるための計算量はＯ（Ｎ）となる。すなわち，上記のＯ（Ｎ）法で求めたエネルギーや力を用いて，構造最適化や分子動力学計算を行うことができる。
　したがって，エネルギーや原子に働く力をより精度良くポイントは，「波動関数を求めずに，近傍に存在する２つの局在基底φiαとφjβからなる密度行列の行列要素ρiα，jβをいかに正しく，かつ，計算量がＯ（Ｎ）となるように求めるか」ということになる。そのような方法として数多くの手法が提案されている。一手法として，分割統治法（Divide-and-Conquer法）が提案されている（W. Yang, "Phys. Rev. Lett. 66", 1991, 1438）。
　分割統治法（ＤＣ法）は，全系を複数の領域に分割し，分割した領域の密度行列を求め，分割した領域の密度行列を用いて分割した領域の物理量を求め，最後に各領域で求めた物理量を統合して，系全体の物理量を求める手法である。
　なお，分子を実空間で領域に分割し，分子の全エネルギーを領域ごとの寄与，領域間の寄与に分割し，分割された各エネルギー項のある領域の寄与を，その領域を含む他の分子の領域の寄与で置き換え，予め用意した近似レベルのうちの所望の近似レベルを選択して計算する手法が提案されている。
特開２０００－１１０２５号公報
　量子論にもとづいて物質の原子構造・電子状態などのシミュレーションを行うために，上述の（１）式のシュレディンガー方程式を解く場合に，オーダーN法（Ｏ（Ｎ）法）を用いることによって，一電子の固有値と波動関数を求めることは不要となり，必要な物理量だけを求めればよい。
　しかし，大規模な分子系を対象とする場合には，以下のような問題がある。
　Ｏ（Ｎ）法によりエネルギーや原子に働く力を正しく求めるためにＤＣ法を採用する場合，分割した領域の密度行列を求める際は，分割した領域内の原子だけではなく，分割した領域に隣接する領域（バッファ領域）の原子も用いて計算する必要がある。
　図１は，分割した領域とバッファ領域の例を示す図である。
　図１は，分子系の原子空間を擬似的に平面で表した図である。分割領域１は，ＤＣ法の適用のために分割された領域である。バッファ領域２は，分割領域１と隣接する領域である。分割領域１とバッファ領域２とを合わせた領域が，密度行列を求める際に要する領域（計算用領域）となる。
　各分割領域１について，具体的には，以下のように物理量が計算される。
　密度行列の行列要素ρiα，jβを，下記の（１６）式のように書き表す。
Figure JPOXMLDOC01-appb-M000015
　ここで，（１６）式のIは，分割した領域のインデックスである。今，Ρij (I)を，以下の（１７）式，（１８）式のように定義する。
Figure JPOXMLDOC01-appb-M000016
　すなわち，局在基底φiαとφjβが，両方とも図１に示す分割領域１にある場合は，係数１を掛け，片方のみ分割領域１にある場合は，係数０．５を掛け，両方とも分割領域１にない場合は，密度行列を０とする。
　そして，分割領域１の密度行列の成分を求めるために，図１に示す分割領域１とバッファ領域２でなる計算用領域を抜き出し，抜き出した計算用領域に含まれる局在基底を用いて，（９）式のように行列成分を求める。
Figure JPOXMLDOC01-appb-M000017
　そして，（１９）式を用いて各計算用領域の固有値問題を解く。
Figure JPOXMLDOC01-appb-M000018
　各計算用領域の固有値問題を解いて求めた波動関数を用いると，分割領域１の密度行列の成分は，以下の（２１）式となる。
Figure JPOXMLDOC01-appb-M000019
　ここで注意すべき点は，分割領域１に対して固有値問題を解いたので，（２１）式の密度行列の要素は，近似した物理量という点である。ただし，計算用領域として全系を含むようにとれば，密度行列の要素は厳密な物理量となる。
　エネルギーなどの物理量を計算するためには，フェルミ分布関数に現れるフェルミエネルギーを求めなければならない。フェルミエネルギーは，「全系の電子数Nが一定」という下記の（２２）式で示す条件から求めることができる。
Figure JPOXMLDOC01-appb-M000020
　ここで，フェルミ分布関数ｆ（εn (I)）は，以下の（２３）式と表される。εＦがフェルミエネルギーである。
Figure JPOXMLDOC01-appb-M000021
　フェルミエネルギーを決定する際，全系の電子数を求める必要があるので，各分割領域１の電子数の総和を求める必要がある。（２１）式と（２２）式で決定したフェルミエネルギーを用いると，既出の（１３）式は下記の（２４）式のように表すことができる。
Figure JPOXMLDOC01-appb-M000022
　また，既出の（１４）式は，（２５）式と表すことができる。
Figure JPOXMLDOC01-appb-M000023
　ここで，ωiα,jβ (I)は，以下の（２６）式である。
Figure JPOXMLDOC01-appb-M000024
　以上のように，計算用領域内に含まれる局在基底を用いて密度行列を求め，密度行列から各分割領域１の電子数を計算してフェルミエネルギーを決定する。その後，各分割領域１の物理量を計算して，その和を求めることによって全系の物理量が計算できる。
　局在基底φiαとφjβからなる密度行列の行列要素Ρiα，jβは，局在基底が近傍の場合には値を持つ。そのため，分割領域１の端（領域のすぐ内側）の密度行列の行列要素を求めるために，バッファ領域２を用意しなければならない。
　さらに，バッファ領域２の広さは，計算する系によって，求めるべき物理量が十分な精度で求めるように設定しなければならない。従来では，バッファ領域２の広さは，経験則にもとづいて手動で計算していた。
　Ｏ（Ｎ）法は，原子数を増やすと分割数が原子数に比例して増加するので，計算量が３乗の計算回数が比例して増えることから呼ばれているように，式で表すと，「計算量はオーダーＮ：計算量＝ａ×Ｎ（Ｎは系全体の原子数）」となる。係数ａは，「分割領域１＋バッファ領域２」の原子数ｎの３乗に比例する。
　したがって，バッファ領域２が広いと係数ａが大きくなるため，できる限りバッファ領域２を狭くして計算時間を短縮できる計算手法が所望されている。特に，実際に実験データと比較しながらシミュレーションを行うためには，計算時間をさらに短縮する必要がある。
　本発明の目的は，計算量が原子数に比例する計算手法において原子空間を領域に分割して計算する処理によって，原子構造・電子状態を計算する場合に，計算量を軽減できる技術を提供することである。
　本発明の一態様として開示する原子構造等計算装置は，物質の初期構造情報を取得する計算条件取得部と，前記物質の原子空間を領域に分割する領域分割部と，仮想原子のパラメータを用意し，前記物質の分割領域各々に対して，前記分割領域を含む計算用領域を示す情報として，前記領域に対して終端する仮想原子のパラメータおよび終端位置，ならびに前記分割領域と前記仮想原子が終端する原子層とのバッファ領域となる原子層を決定する計算用領域決定部と，計算量が原子数に比例する計算処理に分割統治法を適用した計算処理によって，前記物質の原子構造または電子状態を求める場合に，前記計算用領域を用いて計算を行う計算処理部とを備える。
　上記の原子構造等計算装置によれば，分割した領域の端の原子を仮想原子に終端させることによって，領域の端の本来存在しない表面電子状態をなくすことができるので，分割した領域の密度行列を狭いバッファ領域を用いて計算することができる。よって，原子構造・電子状態の計算の計算量軽減を実現することができる。
分割領域とバッファ領域の例を示す図である。 開示する原子構造等計算装置の一実施形態における構成例を示す図である。 仮想原子の終端例を示す図である。 原子構造等計算装置による計算結果を用いるシミュレーション処理フロー例を示す図である。 終端方向と終端位置との例を示す図である。 ステップＳ３の処理の，より詳細な処理フロー例を示す図である。 仮想原子として水素原子のパラメータを用いた場合のバッファ領域の原子層を決定する処理フロー例を示す図である。 原子構造等計算装置１０のハードウェア構成例を示す図である。
　以下に，本発明の一態様として開示される原子構造等計算装置の一実施形態を説明する。
　原子構造等計算装置１０は，ターゲットとする物質の原子構造または電子状態を，計算量が原子数に比例する計算手法において，図１に示すように，原子空間を領域分割して計算する装置である。
　図２は，原子構造等計算装置１０の一実施形態における構成例を示す図である。
　原子構造等計算装置１０は，計算条件取得部１１，領域分割部１２，計算用領域決定部１３，計算処理部１４，および仮想原子情報記憶部１５を備える。
　計算条件取得部１１は，ターゲットとする物質の計算条件，初期原子構造情報を取得する。
　計算条件は，例えば，シミュレーション温度，シミュレーション圧力，シミュレーションする原子種のポテンシャルの指定，分子動力学のタイムステップ数，Ｏ（Ｎ）法による計算精度である。計算精度は，エネルギーＥで判断する場合にはΔＥで，原子に働く力Ｆで判断する場合はΔＦで設定される。
　初期原子構造情報は，例えば，単位格子ベクトル，格子定数，原子数，原子種，原子座標，初期原子速度などである。
　領域分割部１２は，ＤＣ法を実行するために，ターゲットとする物質の全系の原子空間を，直方体に分割する。
　計算用領域決定部１３は，領域分割部１２で分割された領域（分割領域）１各々に対して，分割領域１と隣接するバッファ領域２とでなる計算用領域を決定する。具体的には，計算用領域決定部１３は，各分割領域１に対する計算用領域を示す情報として，領域の端の原子に終端する仮想原子３のパラメータ，終端位置，バッファ領域２の原子層を決定する。
　図３は，仮想原子３の終端例を示す図である。
　図３に示すように，ターゲットとする物質の原子空間の分割領域１に対する計算用領域の端の原子層で各原子が仮想原子３に終端する。図３では，仮想原子３は，価電子数が１である水素原子型の原子に終端する場合である。
　さらに，計算用領域決定部１３は，計算処理部１４の計算結果，すなわち，分割領域１について決定した計算用領域の場合のエネルギーまたは原子に働く力が，予め与えられた計算精度を満たして求まり，かつ，計算用領域が含む原子層数が最小となる場合のバッファ領域２の原子層を求めることができる。計算用領域決定部１３は，このような場合の領域を示す情報（終端する仮想原子３のパラメータ，終端位置，バッファ領域２の原子層）をシミュレーション用情報に決定する。
　計算処理部１４は，計算用領域決定部１３が決定した分割領域１に対する計算用領域（終端する仮想原子３のパラメータ，終端位置，バッファ領域の原子層）の場合について，原子構造または電子状態の計算処理を行い，エネルギーもしくは原子に働く力を求める。計算処理部１４は，１つの分割領域１について，計算用領域決定部１３が決定した様々な大きさの計算用領域に対して上記の計算処理を行う。
　仮想原子情報記憶部１５は，仮想原子３のパラメータと，仮想原子３と全原子種との安定的結合距離を示す原子間距離を含む仮想原子情報を記憶する。
　例えば，仮想原子情報記憶部１５には，仮想原子３のパラメータとして，水素（Ｈ），リチウム（Ｌｉ）等の第１族の元素，フッ素（Ｆ），塩素（Ｃｌ）等の第１７族の元素の形式を示すパラメータが記憶される。
　さらに，仮想原子情報として，第１族の水素（Ｈ），リチウム（Ｌｉ），…の各元素，または，第１７族のフッ素（Ｆ），塩素（Ｃｌ），…の各元素の各々について，全原子種との安定的な結合距離を示す原子間距離が記憶される。
　図４は，原子構造等計算装置１０による計算結果を用いるシミュレーション処理フロー例を示す図である。
　ステップＳ１：　原子構造等計算装置１０において，計算条件取得部１１は，シミュレーションに対する計算条件，ターゲットとなる物質の初期原子構造情報を取得する。
　ステップＳ２：　領域分割部１２は，ターゲットとなる物質の全系（全原子空間）を，Ｏ（Ｎ）法における計算処理をＤＣ法を適用して実行するための領域（分割領域）として直方体に分割する，
　ステップＳ３：　計算用領域決定部１３は，仮想原子情報記憶部１５の情報を参照して，各分割領域１について，計算用領域を示す情報として，終端に用いる仮想原子３のパラメータ，終端位置，バッファ領域２の原子層を決定する。
　また，原子構造等計算装置１０には，少なくとも１つの仮想原子３のパラメータが，用意される。仮想原子３のパラメータとして，例えば，水素原子型のパラメータが用いられる。水素原子型のパラメータは，価電子数が１つであり，全対称のポテンシャルを持つような形式を示すパラメータである。仮想原子３のパラメータは，例えば，水素原子型のｓ軌道と結合する相手の原子ｓ，ｐ，ｄ，…，軌道との相互作用の形式を用いる。
　または，原子構造等計算装置１０は，仮想原子情報記憶部１５を備えていてもよく，この場合に，仮想原子情報記憶部１５に記憶された元素の形式を示すパラメータから，仮想原子３のパラメータが選択される。
　終端位置は，終端方向と原子間の距離で定まる。
　終端方向は，図５（Ａ）に示すように，領域の端の原子（仮想原子３に終端する原子）の位置と，その領域を分割する前に当該原子と結合していた原子の位置との直線上の方向に設定される。
　原子間の距離は，図５（Ｂ）に示すように，領域の端の原子の終端方向で，当該原子と仮想原子３のエネルギーが安定的に結合する距離とする。例えば，メタン（ＣＨ４）のＣ－Ｈの結合長は，１．１Åであるので，領域の端の原子が炭素（Ｃ）であり仮想原子３が水素（Ｈ）であるペアの場合には，結合距離を１．１Åとする。また，シラン（ＳｉＨ４）のＳｉ－Ｈの結合長は，１．５Åであるので，領域の端の原子がケイ素（Ｓｉ）であり仮想原子３が水素（Ｈ）である場合には，結合距離を１．５Åとする。
　バッファ領域２の原子層は，Ｏ（Ｎ）法による原子構造，電子状態の計算において所定の計算制度ΔＥ，ΔＦを満たす十分に広いバッファ領域２の大きさが与えられている場合に，その大きさの範囲内の任意の原子層数をとる。
　ステップＳ４：　計算処理部１４は，ステップＳ３で決定された計算用領域の場合について，Ｏ（Ｎ）法による原子構造，電子状態の計算を行う。より詳しくは，計算処理部１４は，一例として，前述の（１６）～（２６）式を用いて計算を行う。
　ステップＳ３の処理において，計算用領域決定部１３が，１つの分割領域１について，バッファ領域２の原子層数を変える度に，ステップＳ４の処理で，計算処理部１４は，ステップＳ３で決定された計算用領域の処理を繰り返す。そして，計算用領域決定部１３は，与えられた計算精度を満たして原子構造，電子状態の計算処理の結果が求まる場合であって，バッファ領域２の大きさが最小となる場合の計算用領域を示す情報（仮想原子３のパラメータ，終端位置，バッファ領域２の原子層）を，シミュレーション用の計算用領域の情報として採用し，この場合の計算結果をシミュレーション装置２０へ渡す。
　ステップＳ５：　シミュレーション装置２０は，原子構造等計算装置１０から得た計算結果を用いてシミュレーションを行い，シミュレーション結果を出力する。
　図６は，ステップＳ３の処理の，より詳細な処理フロー例を示す図である。
　ステップＳ３の処理において，各分割領域１に対する計算用領域を示す情報，終端する仮想原子３のパラメータ，終端位置，バッファ領域２の原子層は，以下の処理により決定される。
　ステップＳ１０：　計算用領域決定部１３は，ターゲットとなる物質の分割領域１について，仮想原子３のパラメータ＝水素原子型のパラメータ，終端距離＝領域の端の原子と水素原子との原子間距離，終端方向＝分割の前の結合方向を条件として，バッファ領域２の原子層を決定する。
　仮想原子３である水素原子を表すための表式として，例えばタイトバインディングの場合には以下の（２７）式，（２８）式が知られている。
Figure JPOXMLDOC01-appb-M000025
　（２７）式が，（１９）式の飛び移り積分を原子軌道の結合の種類（ｓｓσ，ｓｐσ，ｐｐσ，ｐｐπなど）ごとに分けた表式であり，（２８）式が，（４）式のイオン－イオン相互作用のうち，一組のイオンとイオンの相互作用を表したものである。ここで，αは原子軌道の結合の種類を表し，ｒが原子間距離を表す。そして，ｒ０，ｈα０，ｒＤ，Ｕ１，Ｕ２がパラメータである。
　また，他の表式として，以下の（２９）式，（３０）式が知られている。
Figure JPOXMLDOC01-appb-M000026
　（２７）式および（２８）式と同様に，（２９）式が，（１９）式の飛び移り積分を原子軌道の結合の種類（ｓｓσ，ｓｐσ，ｐｐσ，ｐｐπなど）ごとに分けた表式であり，（３０）式が，（４）式のイオン－イオン相互作用のうち，一組のイオンとイオンの相互作用を表したものである。ここで，αは原子軌道の結合の種類を表し，ｒが原子間距離を表す。
　終端位置は，仮想原子３と終端する原子と，領域分割前に結合していた原子の位置との直線上の方向で，仮想原子３と終端する原子と水素との安定的な結合長（終端距離）の位置とする。
　ステップＳ１１：　計算用領域決定部１３は，処理カウンタｉ＝０とする。
　ステップＳ１２：　計算用領域決定部１３は，予め与えられたバッファ領域２の大きさの場合（例えば，原子層を加えてバッファ領域２となる原子層を十分に大きくした状態の場合）の計算用領域に対する，エネルギーまたは原子に働く力を計算する。
　計算用領域決定部１３は，ステップＳ１０で決定した原子層に，さらに１層原子層を増やし，仮想原子３のパラメータ，終端距離を同じままにして，エネルギーまたは原子に働く力を計算し，これを，バッファ領域２が十分に大きい場合のエネルギーまたは原子に働く力とする。
　ステップＳ１３：　計算用領域決定部１３は，仮想原子３のパラメータと，終端距離とを変化させて，エネルギーまたは原子に働く力を計算する。さらに，計算用領域決定部１３は，計算したエネルギーまたは原子に働く力と，原子層を十分に大きくした状態の場合のエネルギーまたは原子に働く力との差の絶対値の最小値をＡ′を計算する。
　計算用領域決定部１３は，水素原子型の仮想原子３のパラメータ，ｒ０，ｈα０，ｒＤ，Ｕ１，Ｕ２と，領域の端の原子と仮想原子３との距離を変えて，バッファ領域が十分に広い場合のエネルギーもしくは原子に働く力の差の絶対値を計算して，計算精度ΔＥ，ΔＦと比較していく。この場合に，パラメータを変化させる範囲の例は，以下のとおりとする。
　パラメータｒ０：（水素原子型のパラメータ－０．１Å）～（水素原子型のパラメータ＋０．１Å）の範囲とし，刻み幅は０．０１Åとする。
　パラメータｈα０：（水素原子型のパラメータ－１eV）～（水素原子型のパラメータ＋１eV）の範囲とし，刻み幅は０．１eVとする。
　パラメータｒＤ：（水素原子型のパラメータ－０．１Å）～（水素原子型のパラメータ＋０．１Å）の範囲とし，刻み幅は０．０１Åとする。
　パラメータＵ１：（水素原子型のパラメータ－１eV）～（水素原子型のパラメータ＋１eV）の範囲とし，刻み幅は０．１eVとする。
　パラメータＵ２：（水素原子型のパラメータ－１eV）～（水素原子型のパラメータ＋１eV）の範囲とし，刻み幅は０．１eVとする。
　領域の端の原子と仮想原子との距離：（水素原子型のパラメータ－０．１Å）～（水素原子型のパラメータ＋０．１Å）の範囲とし，刻み幅は０．０１Åとする。
　ステップＳ１４：　計算用領域決定部１３は，絶対値の最小値Ａ′が，所定の計算精度であってエネルギーで判断する場合のΔＥまたは原子に働く力で判断する場合のΔＦより小さいかを判定する。
　求めた最小値Ａ′が，計算精度ΔＥまたはΔＦより小さい場合には（ステップＳ１４のＹ），ステップＳ１５へ進み，最小値をＡ′が，計算精度ΔＥまたはΔＦより小さくない場合には（ステップＳ１４のＮ），ステップＳ１９へ進む。
　ステップＳ１５：　計算用領域決定部１３は，バッファ領域の原子層を１層少なくし，ステップＳ１３の処理と同一の仮想原子３のパラメータと終端位置とで，エネルギーまたは原子に働く力を計算する。さらに，計算用領域決定部１３は，計算したエネルギーまたは原子に働く力と，原子層を十分に大きくした状態の場合のエネルギーまたは原子に働く力との差の絶対値の最小値Ａ″を計算する。
　ステップＳ１６：　最小値Ａ″が，計算精度ΔＥまたはΔＦより小さいかを判定する。
　最小値Ａ″が，計算精度ΔＥまたはΔＦより小さい場合には（ステップＳ１６のＹ），ステップＳ１７へ進み，最小値Ａ″が，計算精度ΔＥまたはΔＦより小さくない場合には（ステップＳ１６のＮ），ステップＳ１８へ進む。
　ステップＳ１７：　計算用領域決定部１３は，原子層を１層減らした状態を保持し，処理カウンタｉ＝ｉ＋１とする。
　ステップＳ１８：　計算用領域決定部１３は，ステップＳ１３の処理での仮想原子３のパラメータ，終端位置，バッファ領域の原子層を採用して，処理を終了する。
　ステップＳ１９：　計算用領域決定部１３は，処理カウンタｉが０より大きいか（ｉ＞０）を調べる。処理カウンタｉが０より大きくない場合には（ステップＳ１９のＮ），ステップＳ１１０へ進み，処理カウンタｉが０より大きい場合には（ステップＳ１９のＹ），ステップＳ１１１へ進む。
　ステップＳ１１０：　計算用領域決定部１３は，原子層を１層増やして，ステップＳ１３の処理へ戻る。
　ステップＳ１１１：　計算用領域決定部１３は，ステップＳ１３の処理へ入る直前の状態での仮想原子３のパラメータ，終端位置，バッファ領域２の原子層を採用して，処理を終了する。
　以上の処理により，与えられた処理精度ΔＥ，ΔＦを満たす処理精度であって，バッファ領域２の原子層を少なくした場合の計算用領域を示す，仮想原子３のパラメータ，終端距離，原子層を得ることができる。
　図７は，仮想原子３として水素原子のパラメータを用いた場合のバッファ領域２の原子層を決定する処理フロー例を示す図である。
　ステップＳ２０：　計算用領域決定部１３は，仮想原子３として水素原子のパラメータを用いて，終端位置（結合方向と終端距離）定め，分割領域１についてのバッファ領域２の端の原子層ｉの原子を終端する。
　ステップＳ２１：　原子層の数ｉ＝０とする。
　ステップＳ２２：　計算処理部１４は，ステップＳ２０で決定した計算用領域について，エネルギー，原子に働く力（Ｂｉ）を求める。
　ステップＳ２３：　計算用領域決定部１３は，原子層を１層追加する。
　ステップＳ２４：　原子層の数ｉ＝ｉ＋１とする。
　ステップＳ２５：　原子層を追加する前と同じ方法で，結合方向と終端距離を定める。
　ステップＳ２６：　計算処理部１４は，ステップＳ２５で決定した計算用領域について，エネルギー，原子に働く力（Ｂｉ＋１）を求める。
　ステップＳ２７：　計算用領域決定部１３は，Ｂｉ＋１とＢｉとの差の絶対値を求め，求めた絶対値が，エネルギーまたは原子に働く力に対して与えられている条件（計算精度）より小さいかを判定する。
　ステップＳ２８：　計算用領域決定部１３は，ステップＳ２７で求めた絶対値が，条件を満たしていれば（ステップＳ２８のＹ），ステップＳ２９の処理へ進み，条件を満たしていなければ（ステップＳ２８のＮ），ステップＳ２４の処理へ戻る。
　ステップＳ２９：　計算用領域決定部１３は，原子層の数として，ｉを採用する。
　図８は，原子構造等計算装置１０のハードウェア構成例を示す図である。
　原子構造等計算装置１０は，ＣＰＵ１０１，一時記憶装置（ＤＲＡＭ，Ｆｌａｓｈ　Ｍｅｍｏｒｙ等）１０２，永続性記憶装置（ＨＤＤ，Ｆｌａｓｈ　Ｍｅｍｏｒｙ等）１０３を備えるコンピュータ１００によって実施することができる。
　また，原子構造等計算装置１０は，コンピュータ１００が実行可能なプログラムによって実施することができる。この場合に，原子構造等計算装置１０が有すべき機能の処理内容を記述したプログラムが提供される。提供されたプログラムをコンピュータ１００が実行することによって，上記説明した原子構造等計算装置１０の処理機能がコンピュータ１００上で実現される。
　すなわち，原子構造等計算装置１０の計算条件取得部１１，領域分割部１２，計算用領域決定部１３，計算処理部１４は，プログラムで構成することができ，仮想原子情報記憶部１５は，永続性記憶装置１０３で構成することができる。
　なお，コンピュータ１００は，可搬型記録媒体から直接プログラムを読み取り，そのプログラムに従った処理を実行することもできる。また，コンピュータ１００は，ネットワークインターフェースを備えて，サーバコンピュータからプログラムが転送されるごとに，逐次，受け取ったプログラムに従った処理を実行することもできる。
　さらに，このプログラムは，コンピュータ１００で読み取り可能な記録媒体に記録しておくことができる。
　上記説明したように，原子構造等計算装置１０は，ＤＣ法を用いた量子論にもとづく原子構造・電子状態の計算処理手法において，計算用の領域の端の原子を仮想原子３に終端させるため，バッファ領域２を小さくして計算量を軽減することができる。
　よって，原子構造，電子状態の計算を高速化させることができ，よって，実験結果とシミュレーション結果とを照らし合わせた材料開発を，より効果的に行うことができる。
　また，原子構造等計算装置１０は，原子構造・電子状態の高速演算が可能になるため，取り扱える原子数を拡張することができる。具体的な例としては，シリコンのダイヤモンド結晶において分割領域に８原子含む場合に，従来，分割領域とバッファ領域の原子数は２１６原子が必要であったが，開示する原子構造等計算装置１０によれば，１６２原子で済み，計算速度は約２．４倍向上する。
　さらに，以下の２点の効果がある。
　・　開示する原子構造等計算装置１０によれば，計算用領域の一部をなすバッファ領域の決定処理は，分割領域ごとに独立して行うことができる。よって，各分割領域について並列計算による高速化が容易になる。
　・　開示する原子構造等計算装置１０が実行するＤＣ手法を用いる計算処理の高速化手法は，第一原理計算または半経験的量子力学計算のいずれにも適応することができる。
　したがって，開示する原子構造等計算装置１０により，新材料やデバイスの開発において，実験結果とシミュレーション結果とを照らし合わせた開発効率を向上させることができる。具体的には，目的とする材料を得るための温度や圧力，材料の組成比を決定するための指針をシミュレーション結果から得る際に，従来に比べて短時間で済むため，開発速度を向上させることができる。
　よって，開示する原子構造等計算装置１０は，開発時間と経費の削減，開発による環境負荷の低減に寄与する。
　１０　原子構造等計算装置
　１１　計算条件取得部
　１２　領域分割部
　１３　計算用領域決定部
　１４　計算処理部
　１５　仮想原子情報記憶部
　２０　シミュレーション装置
 

Claims (7)

  1. 　物質の初期構造情報を取得する計算条件取得部と，
    　前記物質の原子空間を領域に分割する領域分割部と，
    　仮想原子のパラメータを用意し，前記物質の分割領域各々に対して，前記分割領域を含む計算用領域を示す情報として，前記領域に対して終端する仮想原子のパラメータおよび終端位置，ならびに前記分割領域と前記仮想原子が終端する原子層とのバッファ領域となる原子層を決定する計算用領域決定部と，
    　計算量が原子数に比例する計算処理に分割統治法を適用した計算処理によって，前記物質の原子構造または電子状態を求める場合に，前記計算用領域を用いて計算を行う計算処理部とを備える
    　ことを特徴とする原子構造等計算装置。
  2. 　前記仮想原子のパラメータは，価電子数が１であり全対称のポテンシャルを有する形式を示すものである
    　ことを特徴とする請求項１に記載の原子構造等計算装置。
  3. 　前記仮想原子のパラメータの決定は，第１族または第１７族に属するいずれかの元素のパラメータを初期値に用いる
    　ことを特徴とする請求項１または請求項２に記載の原子構造等計算装置。
  4. 　前記仮想原子のパラメータとともに，前記仮想原子と全原子種との安定的結合を生じる原子間距離を設定した仮想原子情報を記憶する仮想原子情報記憶部を備えて，
    　前記計算用領域決定部は，前記仮想原子の終端位置を定める前記仮想原子と前記仮想原子に終端する原子との距離の決定に，前記仮想原子情報にもとづいて，前記仮想原子と終端する原子種に対応する原子間距離を初期値に用いる
    　ことを特徴とする請求項１ないし請求項３のいずれか一項に記載の原子構造等計算装置。
  5. 　前記計算処理部は， 
    　前記計算用領域決定部は，前記物質に対応する計算精度を取得し，前記計算処理部による計算結果が前記計算精度を満たし，かつ，前記バッファ領域の原子層の数が最小の場合の計算用領域を示す，前記計算用領域に対して終端する仮想原子のパラメータおよび終端位置，ならびに前記バッファ領域となる原子層を，シミュレーション用情報として決定する
    　ことを特徴とする請求項１ないし請求項４のいずれか一項に記載の原子構造等計算装置。
  6. 　物質の原子構造または電子状態を計算する方法であって，
    　コンピュータが，
    　物質の初期構造情報を取得する処理ステップと，
    　前記物質の原子空間を領域に分割する処理ステップと，
    　仮想原子のパラメータを用意し，前記物質の分割領域各々に対して，前記分割領域を含む計算用領域を示す情報として，前記領域に対して終端する仮想原子のパラメータおよび終端位置，ならびに前記分割領域と前記仮想原子が終端する原子層とのバッファ領域となる原子層を決定する処理ステップと，
    　計算量が原子数に比例する計算処理に分割統治法を適用した計算処理によって，前記物質の原子構造または電子状態を求める場合に，前記計算用領域を用いて計算を行う処理ステップとを，実行する
    　ことを特徴とする原子構造等計算方法。
  7. 　コンピュータに，物質の原子構造または電子状態を計算する処理を実行させるためのプログラムであって，
    　物質の初期構造情報を取得する処理と，
    　前記物質の原子空間を領域に分割する処理と，
    　仮想原子のパラメータを用意し，前記物質の分割領域各々に対して，前記分割領域を含む計算用領域を示す情報として，前記領域に対して終端する仮想原子のパラメータおよび終端位置，ならびに前記分割領域と前記仮想原子が終端する原子層とのバッファ領域となる原子層を決定する処理と，
    　計算量が原子数に比例する計算処理に分割統治法を適用した計算処理によって，前記物質の原子構造または電子状態を求める場合に，前記計算用領域を用いて計算を行う処理とを，実行させるための
    　原子構造等計算プログラム。
     
PCT/JP2010/069304 2010-10-29 2010-10-29 原子構造等計算装置，方法，およびプログラム WO2012056562A1 (ja)

Priority Applications (1)

Application Number Priority Date Filing Date Title
PCT/JP2010/069304 WO2012056562A1 (ja) 2010-10-29 2010-10-29 原子構造等計算装置，方法，およびプログラム

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
PCT/JP2010/069304 WO2012056562A1 (ja) 2010-10-29 2010-10-29 原子構造等計算装置，方法，およびプログラム

Publications (1)

Publication Number Publication Date
WO2012056562A1 true WO2012056562A1 (ja) 2012-05-03

Family

ID=45993316

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/JP2010/069304 WO2012056562A1 (ja) 2010-10-29 2010-10-29 原子構造等計算装置，方法，およびプログラム

Country Status (1)

Country Link
WO (1) WO2012056562A1 (ja)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104805846A (zh) * 2015-04-24 2015-07-29 成都理工大学 一种浅表层土质滑坡的危险性划分方法及应用
WO2019014309A1 (en) 2017-07-13 2019-01-17 Verdezyne (Abc), Llc BIOLOGICAL METHODS FOR MODIFYING A CELL CARBON FLOW

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
HIROYUKI KAGESHIMA: "<Kiso-hen> Kiso Koza <Bunshi Simulation no Kiso to Oyo>", FIRST-PRINCIPLES CALCULATION 1, OYO BUTSURI, vol. 75, no. 10, 10 October 2006 (2006-10-10), pages 1258 - 1265 *
MASATO KOBAYASHI ET AL.: "Implementation of Divide-and-Conquer (DC) Electronic Structure Code to GAMESS Program Package", JOURNAL OF COMPUTER CHEMISTRY, vol. 8, no. 1, 2009, JAPAN, pages 1 - 12 *
MASAYUKI HATA ET AL.: "Diffusion and Nucleation Mechanism of Silicon Adatoms on Si (001) Surface", EXTENDED ABSTRACTS, JAPAN SOCIETY OF APPLIED PHYSICS AND RELATED SOCIETIES, DAI 43 KAI, no. 2, 26 March 1996 (1996-03-26), pages 697 *

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104805846A (zh) * 2015-04-24 2015-07-29 成都理工大学 一种浅表层土质滑坡的危险性划分方法及应用
WO2019014309A1 (en) 2017-07-13 2019-01-17 Verdezyne (Abc), Llc BIOLOGICAL METHODS FOR MODIFYING A CELL CARBON FLOW

Similar Documents

Publication Publication Date Title
Momeni et al. Multiscale computational understanding and growth of 2D materials: a review
Perkyns et al. Protein solvation from theory and simulation: Exact treatment of Coulomb interactions in three-dimensional theories
Xiao et al. Solid-state dimer method for calculating solid-solid phase transitions
Moses et al. Density functional study of the adsorption and van der Waals binding of aromatic and conjugated compounds on the basal plane of MoS2
Fraux et al. Note: Heterogeneous ice nucleation on silver-iodide-like surfaces
Steinmann et al. A generalized-gradient approximation exchange hole model for dispersion coefficients
Zervos et al. Two finite-element discretizations for gradient elasticity
Cui et al. A second nearest-neighbor embedded atom method interatomic potential for Li–Si alloys
Ismail et al. Capillary waves at the liquid-vapor interface and the surface tension of water
Jeanmairet et al. Molecular density functional theory of water describing hydrophobicity at short and long length scales
dos Santos et al. Surface tensions and surface potentials of acid solutions
CN106372400B (zh) 构建极化力场的方法及应用、预测药物晶型的方法及系统
Ohwaki et al. Large-scale first-principles molecular dynamics for electrochemical systems with O (N) methods
Aviat et al. The truncated conjugate gradient (TCG), a non-iterative/fixed-cost strategy for computing polarization in molecular dynamics: Fast evaluation of analytical forces
Chakrabarti et al. Fully coupled discrete energy-averaged model for Terfenol-D
WO2012056562A1 (ja) 原子構造等計算装置，方法，およびプログラム
Dubecký et al. Benchmarking fundamental gap of Sc2C (OH) 2 MXene by many-body methods
Scemama et al. An efficient sampling algorithm for variational Monte Carlo
Otsuka et al. Linear-scaling first-principles molecular dynamics of complex biological systems with the Conquest code
Houriez et al. Ion hydration free energies and water surface potential in water nano drops: The cluster pair approximation and the proton hydration Gibbs free energy in solution
Feng et al. Communication: Molecular dynamics simulations of the interfacial structure of alkali metal fluoride solutions
Kalemos et al. Theoretical investigation of the He–I2 (E3Πg) ion-pair state: Ab initio intermolecular potential and vibrational levels
Itakura et al. Atomistic study on the cross-slip process of a screw< a> dislocation in magnesium
Parandekar et al. Applications and assessment of QM: QM electronic embedding using generalized asymmetric Mulliken atomic charges
Di Prinzio et al. Molecular dynamics simulations of tilt grain boundaries in ice

Legal Events

Date Code Title Description
121 Ep: the epo has been informed by wipo that ep was designated in this application

Ref document number: 10858951

Country of ref document: EP

Kind code of ref document: A1

NENP Non-entry into the national phase

Ref country code: DE

122 Ep: pct application non-entry in european phase

Ref document number: 10858951

Country of ref document: EP

Kind code of ref document: A1

NENP Non-entry into the national phase

Ref country code: JP