Crystal Structure ModelAssembly Program Using the Monte Carlo and the Rfactor Method
Hiroyuki MIURA and Takeshi KIKUCHI
Return
1 Introduction
The Rietveld method is an effective tool to determine the precise crystal structure of powder crystalline materials. It refines crystal structure parameters such as atomic coordinates and temperature factors. However, the Rietveld method can not find structure models by itself, since it requires a starting model. Relatively correct initial data must be prepared to solve the nonlinear equation. In other words, the structure model of the target material must be provided. When no structural information about the material is available, it is very difficult to make a structure model. In a previous paper, we proposed a calculation program based on a packing analysis method [1]. This program searches suitable atom arrangements in a given unit cell from a simple rule stating that "interatomic distances must be larger than the sum of each atomic radius." This method was useful to find the atomic arrangements of simple inorganic materials. However, in the case of complicated materials, it could not finish the calculation within reasonable time because the program must check all possible arrangements. We here propose a new method based on the combination of the Monte Carlo method and the Rfactor method which can be applied to a crystal structure with several independent atoms. This program makes crystal structure models from XRD data, cell constants, a space group, a chemical formula, and a Znumber in a reasonable amount of time using a personal computer.
2 How to construct structure models
Indexed powder diffraction data, from which cell constants are calculated, are indispensable to start crystal structure analysis. Candidates of the space group are narrowed down to several groups by cell constants and extinction of diffraction. Whenever available, information about the symmetry of the diffraction pattern from TEM is useful. A chemical formula and specific gravity lead to a number of formulae in a unit cell (Znumber). Optimizing isotropic temperature factors is not necessary because intensities are not significantly affected by temperature factors. What is necessary is to find the correct positions of independent atoms in a unit cell. Our proposed method makes structure models as follows. First, the program arranges atomic coordinates in an asymmetric unit at random and calculates the theoretical Xray powder diffraction pattern on the basis of the space group, the cell constants, and the isotropic temperature factors (B=1.0). Then, the program calculates the Rfactor of this model from the theoretical and the observed Xray diffraction data. Based on the Rfactors, the structure models are selected and stored. In the calculation process, the following points should be considered.
2. 1 Axis setting
In a tetragonal or a hexagonal system, for example, caxes must be 4fold or 6fold, respectively, and are uniquely determined from powder diffraction data. However, in some cases, the direction of the unit cell axes of an orthorhombic system is not determined uniquely. In this case, one axis is named in three ways as a, b, or caxis and six cell settings are possible. If a space group has glide planes, a distinction must be observed in specific directions. For example, forsterite, whose space group is Pnma, shows the following extinction, viz. (1)0kl:k+l=2n, (2)hk0:h=2n, (3)h00:h=2n, (4)0k0:k=2n, and (5)00l:l=2n. The observed peak of d=3.880A, indexed as 210 from the cell constants a=10.195, b=5.981, and c=4.756, satisfies the abovementioned extinction law. If we indexed by other cell constants, a=4.756, b=5.981, c=10.195, the index of the same diffraction peak would become 021, which does not satisfy the extinction law, and the calculated intensity of this peak would be zero. In this case, it is easy to eliminate wrong axis settings. However, when it is not possible to determine the correct axis setting by the extinction law, all possible arrangements must be tried.
2. 2 Wyckoff position
If there are many probable Wyckoff positions for one atom, one position must be selected from all candidates. Let us take forsterite as an example: it has the formula, Mg_{2}SiO_{4}, its Znumber is 4, and it belongs to space group Pnma. There are eight magnesium (Mg) atoms, four silicon (Si) atoms, and sixteen oxygen (O) atoms in a unit cell. Wyckoff positions of Pnma are 4a, 4b, 4c, and 8d [2]. As there are eight Mg atoms in a unit cell, Mg atoms make it possible to locate two different Wyckoff positions of multiplicity 4 (described as 4+4) or one 8d(8) position. Si atoms must be located at the positions of multiplicity 4(4). Oxygen atoms must be located by combinations such as (8+8), (4+4+8), or (4+4+4+4). The total number of combinations is 2×1×3=6 sets. All these combinations must be checked.
2. 3 Asymmetric unit
The search area must be restricted to an asymmetric unit, and positions of all other symmetrically dependent atoms are calculated by symmetry operations. Sometimes, a pair of models brought into coincidence by translation, rotation, or reflection to the other is found. Two models which are related by translation share the same structure with a different origin. When the space group only has the symmetry operation of the first kind, a pair of models related by mirror or inverse is produced. They are an enantiomorph of the other. Restricting the search area to an asymmetric unit of the Chesire group is effective to reduce the total calculation time and eliminate duplicates[3].
2. 4 Discrimination function
Diffraction intensities mainly depend on atomic coordinates. When the number of diffraction peaks is not enough in comparison to the number of atomic coordinates which we can move as variables, the obtained results frequently converge to a local minimum. Many peaks overlap in powder data, and it is not easy to increase the observed diffraction peaks. We introduced a discriminant to avoid falling to a local minimum. The worst structure models show unrealistic atomic arrangements such as improper localization of cations, unreasonable coordination numbers, or heavy distortion in coordination polyhedral. Using the following criteria, incorrect models are eliminated: (1) Pauling's 1st rule, which states that the interatomic distances must be larger than the sum of each atomic radius, must be applied; (2) Crystal lattice energy must be reasonably small.
2. 5 Calculation process
A structure modelassembly program was then built taking the above points into account. After users indicate the possible Wyckoff position arrangements, the program searches the suitable structure model automatically by the following process:
[preparation]
(1) Produce all possible indexes and sort them by dvalue;
(2) Calculate the multiplicity of diffraction and the Lp factors of each index;
[1st stage] Search in an asymmetric unit
(3) Choose Wyckoff positions for each independent atom at random;
(4) Set the variable atomic coordinates by a random number;
If the interatomic distances are larger than the sum of each atomic radius, the model is accepted;
(5) Calculate the theoretical XRD pattern and Rfactor;
(6) Sort the obtained models by the Rfactor in small order and store the top 300 models;
(7) When the Rfactor becomes smaller than a value fixed in advance, the process proceeds to the 2nd stage;
[2nd stage] Search the atomic position in a restricted area.
(8) Select one structural model made in the first stage and move each atom around the initial positions at random. If the Rfactor becomes smaller than that of the initial model, the newly obtained model is adopted;
(9) Repeat the above calculation and, if the Rfactors converge, proceed to the next stage;
[3rd stage]
(10) Calculate the lattice energy and eliminate incorrect structural models;
If necessary, repeat the above cycle two or three times;
(11) Print out the obtained structure models.
3 Example for calculation
Table 1 Input data for brookite
Formula  TiO_{2}

Znumber  8

Space Group  Pbca

Cell constants  a=9.182  b=5.456  c=5.143A

Atoms  valence  atomic radius  Wyckoff position

Ti1  4  0.70  8c

O1  2  1.30  8c

O2  2  1.30  8c

For Xray diffraction data used in calculation, see Dobs. and Iobs in Table 2.
Table 2 Obtained structure models for brookite
Rfactor= 0.046

 obtained structure  reference structure*

atom  w.p.  x  y  z  x  y  z

Ti  8c  0.129  0.104  0.138  0.1290  0.0972  0.1371

O  8c  0.006  0.093  0.199  0.0101  0.1486  0.1824

O  8c  0.237  0.161  0.461  0.2304  0.1130  0.4629


Dobs.  Dcal.  Iobs.  Ical.  H  K  L

3.510  3.513  100  100  2  1  0

3.470  3.466  80  81  1  1  1

2.900  2.901  90  90  2  1  1

2.729  2.728  4  3  0  2  0

2.476  2.476  25  23  1  0  2

2.409  2.410  18  18  0  2  1

2.370  2.369  6  8  3  1  1

2.344  2.345  4  4  2  2  0

2.332  2.331  4  3  1  2  1

2.296  2.295  5  3  4  0  0

2.254  2.255  8  9  1  1  2

2.244  2.243  18  18  2  0  2

2.133  2.134  16  15  2  2  1

1.969  1.969  16  21  3  0  2

1.893  1.893  30  30  3  2  1

  1.871    2  0  2  2

1.851  1.852  13  15  3  1  2

*:reference structure is derived from Wyckoff[5]
Figure 1. Convergence of atomic coordinates of brookite. Ti, O1 and O2 positions in asymmetric unite of best one hundred models after 3×10^{3,} 3×10^{4}, 3×10^{5} and 3×10^{6} Monte Carlo trials are projected on (001) plane. All atoms concentrate to each circles in the asymmetric unit which shows the real atomic position.
3. 1 Brookite (TiO_{2})
The structure model of brookite was produced by the crystal data described in JCPDS Card No. 291361[4], including the space group, the chemical formula, the Znumber, and the cell constants (Table 1). The space group of brookite reported in JCPDS Card No. 291361, Pcab, was changed to Pbca, the standard lattice setting described in International Tables Vol. A (IT). Cell parameters were changed to a=9.189, b=5.4558, and c=5.1429A. There are eight titanium atoms and sixteen oxygen atoms in a unit cell. As the multiplicity of general position is eight, it is assumed that one independent Ti atom and two independent O atoms exist on 8c positions. In the first stage, 300 models were selected from 600 thousand produced by a random number. The change of the coordinates of independent atoms in an asymmetric unit (0<=x,y,z<=1/2) of the best 100 models was projected on 001 planes (Figure 1). As the number of loops increased, the area in which each atom exists became clear. In the second stage, the search area for the independent atoms was narrowed down to the vicinity of the atoms of the obtained models, and the Rfactors were dramatically reduced (Figure 2). The Rfactor decreased linearly to the logarithm of the Monte Carlo number. If the search area were not reduced, 6×10^{7} Monte Carlo numbers would be required. This process reduces calculation time dramatically as well (1/30). Eight models giving the best Rfactors were obtained. Atoms were situated in mirror symmetry by the (100) plane with x=1/4, the (010) plane with y=1/4, and the (001) plane with Z=1/4. These models were in enantiomorphic relation. Searching all asymmetric units of the space group was not necessary. The Ti atom could be found in the asymmetric units of the Chesire group (0<=x<=1/4, 0<=y<=1/4, 0<=z<=1/4). The finally obtained crystal structural model and the observed and calculated XRD data after 2.58 million loops are listed in Table 2.
Figure 2. Change of Rfactor. Rfactor decrease linearly to the logarithm of Monte Carlo number in the first stage and decrease dramatically after the search range was narrowed.
3. 2 Forsterite (Mg_{2}SiO_{4})
The crystal structure of forsterite was solved from the data of JCPDS 211260. The initial data used for the calculation are listed in Table 3. There were 6 possible arrangements of atom positions, as described earlier. The Rfactor convergence of 6 arrangements was compared from 30 thousand models, and arrangements of Wyckoff positions (4+4) for Mg and (4+4+8) for O were selected for further calculation because they gave the best Rfactor. After 30 million Monte Carlo trials, the program proceeded to the second stage. Although no theoretical parameters were available to indicate the exact timing to start the second stage, this timing was estimated by the experimental parameter of the difference between the Rfactors of the top and the 300th model. If the difference became smaller than 0.1, it could be concluded that the calculation converged. In the second stage, the same process mentioned in the case of brookite was applied. In this stage, the Rfactor decreased dramatically and reached the real structure model. The used XRD data from JCPDS 211260 and the calculated XRD data from the obtained crystal model are listed in Table 4.
Table 3 Input data for forsterite
Formula  Mg_{2}SiO_{4}

Znumber  4

Space Group  Pnma

Cell constants  a = 10.195  b = 5.981  c = 4.756A

Atoms  valence  atomic radius  Wyckoff position

Mg1  2  0.40  4a 4c

Mg2  2  0.40  4b 4c

Si1  4  0.40  4a 4b 4c

O1  2  1.30  4a 4b 4c

O2  2  1.30  4a 4b 4c

O3  2  1.30  8d

For Xray diffraction data used in calculation, see Dobs. and Iobs. in Table 4.
Table 4 Obtained structure model for forsterite
Rfactor= 0.118

 obtained structure  reference structure*

atom  w.p.  x  y  z  x  y  z

Mg  4a  0  0  0  0  0  0

Mg  4c  0.283  0.250  0.000  0.280  0.250  0.013

Si  4c  0.101  0.250  0.443  0.098  0.250  0.433

O1  4c  0.086  0.250  0.783  0.092  0.250  0.790

O2  4c  0.464  0.250  0.211  0.449  0.250  0.219

O3  8d  0.168  0.022  0.272  0.163  0.038  0.277


Dobs.  Dcal.  Iobs.  Ical.  H  K  L

5.100  5.098  20  20  2  0  0

4.310  4.310  1  2  1  0  1

3.880  3.880  64  62  2  1  0

3.720  3.723  22  14  0  1  1

3.500  3.497  15  13  1  1  1

3.480  3.477  14  15  2  0  1

3.010  3.006  6  4  2  1  1

2.991  2.991  18  15  0  2  0

2.764  2.765  62  0  3  0  1

2.579    1     

2.509  2.510  80  65  3  1  1

2.457  2.457  100  100  1  2  1

  2.378    1  0  0  2

2.344  2.345  12  16  4  1  0

2.316  2.316  10  6  1  0  2

2.267  2.267  44  41  2  2  1

2.246  2.245  30  32  3  2  0

2.160  2.155  20  20  2  0  2

  2.103    1  4  1  1

2.030  2.027  5  1  2  1  2

*:reference structure is derived from Wyckoff[5]
4 Conclusion
The combination of the Monte Carlo method and the Rfactor method can solve simple oxide and silicate structures on the basis of indexed powder diffraction data and chemical formulae. This method is useful for Xray crystal structural analysis. Atom positions in a unit cell were searched approximately in the first stage and precisely in a restricted area in the second stage. This process, which reduces the total calculation time, can be carried out with a small personal computer. Its most significant advantage is that it does not require preliminary knowledge about the structure. As machine power increases, our method will become more efficient to solve more complex structures.
This work was supported by a grantinaid for Scientific Research (Project number, 10640458) from the Ministry of Education, Japan, which is gratefully acknowledged.
References
[ 1] Miura, H., Journal of the Ceramic Society of Japan, 105, 536540 (1997).
[ 2] International Tables for Crystallography Vol.A, Ed. by Theo Hahn, D.Reidel Publishing Co., Dordrecht (1987).
[ 3] Sakurai, T., A practical Guide for Xray crystal structure analysis, Syokabo (1983), p144145.
[ 4] Mineral Powder Diffraction File, JCPDS International Center for Diffraction Data (1985).
[ 5] Crystal Structures, Ed. by Wyckoff, R. W., Interscience, New York (1965).
Return