The Rigaku Journal
Vol. 6/ No. 1 / 1989

CONTRIBUTED PAPERSRIETAN: A SOFTWARE PACKAGE FOR THE RIETVELDANALYSIS AND SIMULATION OF X-RAY AND NEUTRONDIFFRACTION PATTERNSFUJIO IZUMINational Institute for Research in Inorganic Materials, 1 - 1 Namiki, Tsukuba, Ibaraki 305, Japan

A computer software package, RIETAN, has been developed for the application of the Rietveld method
to angle-dispersive X-ray and neutron powder data. It incorporates large data bases storing information
about space groups, atomic scattering factors, anomalous dispersion corrections, coherent scattering
lengths, etc.; all space groups can be handled without additional programming. A profile shape function has
been implemented which is the sum of a Gauss function and a Lorentz function with unequal peak heights
and widths at one-half peak height. Three different algorithms for nonlinear least squares have been
employed: Gauss-Newton, modified Marquardt, and conjugate direction methods. The method of leastsquares can be readily changed one after another in the course of refinement. These three methods are
combined with an exterior penalty function method so that objective functions may be minimized under
nonlinear equality and inequality constraints. One can not only plot the observed and calculated intensities
after the refinement butalso simulate diffraction patterns of pure substances and mixtures from user-supplied
profile and structure parameters. A new version of RIETAN has been written to make it possible to a nalyze
time-of-flight neutron data taken on a high-resolution powder diffractometer at KENS.

1. Introduction

The Rietveld method of structure refinement withpowder diffraction patterns was originally devised forthe analysis of fixed-wavelength (angle-dispersive)
neuron data [1] [2] but now constitutes a major
breakthrough for the usefulness of X-ray powder data[3]. This powerful method has greatly extended theamount of structural detail which can be obtained
routinely from powder diffraction patterns. It is applied to the investigation of crystal structures whensingle crystals cannot be prepared at all or when theproperties of interest of the single-crystal form differ
from those of the polycrystalline form, for example, incatalysts, sensors, and bioceramics. However, several
problems limit the straightforward application of themethod: restricted information contained in powder
diffraction patterns, difficulty in correction of preferred orientation, inadequacy of profile shape functions so far used, flat and local minima, narrow convergence range, and inaccurate standard deviations[4].
A computer program, XPD, for the Rietveld analysis of X-ray powder data [5] has been extensively
modified to overcome these difficulties to a considerable extent and to permit the refinement with neutron

10

powder data as well [6]. This new version, called
RIETAN (RIETveld ANalysis), plus XPD are distributed to more than 60 research groups, domestic andforeign. It incorporates almost all convenient featuresof a standard program coded by Wiles and Young [7],for example, handling of anisotropic thermal parameters and the two wavelengths in the Kα doublet,
single-pass operation, multiphase capability, refinement.of lattice and preferred-orientation parameters,
and R factors in various shells. In the description that
follows, emphasis is therefore placed on features thatare not implemented in their program.
2. Hardware and OS Requirements
The RIETAN system has been designed for userconvenience and for use by those who have access toany of standard large computers. It can be run on
FACOM, HITAC, ACOS, VAX, and UNIX (4.2BSD) machines without any modifications. The
virtual storage necessary to process a standard job isless than 2.5 Mbytes. Laser-beam printers with densities larger than 300 dpi are desirable to obtain Rietveld refinement patterns (Fig. 1) with good quality.
In future, RIETAN will be run on personal computers (CPU: 80386 microprocessor) equipped with
80387 co-processors and operating under OS/2, Ver.
The Rigaku Journal

Fig. 1 X-Ray Rietveld refinement of fluorapatite (Cu Kα). The solid line is the best-fit profile, and points superimposed on it
are raw data. The differences between the observed and calculated intensites, ∆y=yi(o)-yi(c), are shown by points appearing
at the bottom. Tick marks below the profile indicate the positions of all allowed Kα 1 and Kα2 peaks.

1.1 or above. Disk files storing X-ray intensity data
written in a RIETAN format can be read directly by
the program without any conversion of the files.
Although it takes longer times to obtain the final
solution on the personal computers than on standard
large ones, this OS/2 version must be fairly practicaland easy to use owing to advanced features of theOS/2 operating system such as hierarchical file system, powerful full-screen editor, overlapping-type
window system, and multi-task capability.
3. Programming Considerations
The RIETAN program has been written in
FORTRAN 77 in conformity with ANSI X3.9-1978
standards for portability, quite independently of the
program developed by Rietveld [1]. In order tofacilitate reading the code and user-desired changes,
the whole program is segmented into 1 main programand a little more than 100 subprograms, and many
com-ment lines were sprinkled liberally throughout
ca. 9000 lines of the program. FORTRAN 77 features
which are not supported in FORTRAN IV are fullyVol. 6 No. 1 1989

utilized: block IF (block IF, ELSE, ELSE IF, and
END IF statements), character type, character substring, list-directed READ statement, internal data set
input/output statement, and PARAMETER statement.
The use of block IF makes it possible to reduce
the number of GO TO statements to a great extent andfacilitates understanding of the algorithms of the
program. Except for lines in which only a character
variable is input, all data are input in free format: data
placed in a line are separated by a space or a comma;'/' is placed as the last character of the line if default
values are assigned for the rest of data to be input.
Maximum numbers of reflections, step-scanned
intensity data, total parameters, refinable parameters,chemical species, etc. are declared in PARAMETER
statements and used as declarator subscripts of many
arrays. Accordingly, the dimensions of these arrays in
the whole text can be changed very easily and without
mistakes by using a text editor.

11

4. Data Bases
The following data stored in three sequential files
are input during the execution of RIETAN.
1) Laue group numbers, presence/absence of centersof symmetry, Hermann-Mauguin space groupsymbols, symmetry conditions, and coordinates ofgeneral equivalent positions described for all thesettings of 230 space groups in International
Tables, Vol. 1 [8].
2) Corresponding descriptions in International
Tables, Vol. A [9].
3) Coefficients for analytic approximations of
scattering factors and anomalous dispersion
corrections listed in International Tables, Vol. IV[10], coherent scattering lengths, incoherent
scattering cross sections, and absorption crosssections [11], and atomic weights.
One merely enters names of constituent atoms
and space group and setting numbers in conformity
with International Tables, Vol. 1 [8] or Vol. A [9].
After symmetry conditions have been input directlyfrom File 1) or 2), indices of possible reflections are
generated from them with KDRREF subprograms
written by Cooper, Rouse and Sayers [12]; no preparatory program is needed to generate unique reflection sets as in the original program of Rietveld [1]. If
non-standard settings are required which are not described in International Tables [8][9], symmetry
conditions may be specified directly in the input data.
A further option is also provided to allow the suppression of specified reflection type; this facility isused when atoms are at special positions only. Afterall possible reflections have been generated, theirmultiplicities, interplanar spacings, diffraction angles,hs, h·ts [13], etc. are calculated and then sorted onincreasing scattering angle, 2θ.
As described above, this program calls for only a
restricted amount of input data in free format owing to
the presence of the data bases and use of list-directed
READ statements. It is, therefore, especially suitedfor those who are not familiar with crystallography or
FORMAT specification codes in FORTRAN.
5. Model Function
The observed intensity, yi(o), at a particular step,i, is modeled by a calculated intensity, y i(c):yi(c)=∑ s|Fk|2mkPkL(θk)G(∆θik)+yib(c)

(1)

k

with
|Fk|2 = |Fk(cryst.)|2 + |Fk(magn.)|2

12

(2)

∆θik = θi - θk
(3)
In Eqs. (1), (2), and (3), k=reflection number, s =scale factor, Fk = structure factor, mk = multiplicity, Pk= correction factor for preferred orientation, L(θk) =Lorentz and polarization factors,* G(∆θik) = profile
shape function, y ib(c) = background, F k(cryst.) =crystal-structure factor, Fk(magn.) = magneticstructure factor,** θi = scattering angle at the ith step,θk = Bragg angle. The summation in Eq. (1) must be
carried out over all reflections contribut-ing to the net
intensity at the ith step. The diffraction intensity ismultiplied by the scale factor, s, prior to thesummation because RIETAN can deal with mixturesof two or more phases; s has, of course, a constantvalue for all the reflections of a phase.
6. Structure Factors
Several crystal-structure parameters are contained in Fk(cryst.): fractional coordinates (xj, yj, zj),
occupation factors (nj), overall thermal parameters(Q), individual isotropic thermal parameters (Bj), andanisotropic thermal parameters (β11j, β22j, β33j, β12j,
β13j, β23j), where j is the site number. Either isotropic
or anisotropic thermal motion can be arbitrarily
assigned to each atom.
The CHARACTER data of the coordinates ofgeneral equivalent points read in from the data baseare converted into a set of rotation matrices and translation vectors. The program then automatically selects symmetry operations required for each site in an
asymmetric unit by checking whether or not generalequivalent positions calculated for the site overlap
with each other; only these essential sym-metry
operations are used for the subsequent structure factor
calculations. This procedure greatly reduces the times
required for the calculation of structure factors andtheir derivatives especially when dealing with crystals
with high symmetry. Furthermore, occupation factors
can be input and refined directly without considering
the site multiplicities arising because one or more
symmetry elements act at the site.
The present version can handle magnetic
structures with collinear spin arrangements. The magnetic moment, µj, of each site, the angle between thespin direction and the unique axis of the lattice(tetragonal, hexagonal, and rhombohedral), Φ, andthe angles between the spin direction and three axes
(orthorhombic), Φa, Φb, and Φc, are directly refined* The polarization factor is unnecessary in neutron andsynchrotron X-ray diffraction.
** Fk(magn.) = 0 in X-ray diffraction.

The Rigaku Journal

by using equations given by Shirane [14]. Implementation of magnetic scattering for non-collinear
spin arrangements is roundly desired; this will be
undertaken some day.
7. Preferred-Orientation Function
Preferred orientation is corrected by the following
two-parameter function [15] [16]:Pk =p1 + (1 - pl) exp(-p2φk)(4)
with• Flat-plate sample
φk =ψk (plate crystals)
π/2 - ψk (needle crystals)• Cylindrical sample
Φk=
π/2 - ψk (plate crystals)
ψk (needle crystals)
where ψ k is the acute angle between the reciprocal
vector for the kth reflection and preferred-orientation
vector, hpa*2 + kpb *2 + ipc*2. The direction of the preferred-orientation vector corresponds to that of anormal to the cleavage plane in the plate crystal, and tothat of the extension axis in the needle-shaped crystal.
Equation (4) includes the preferred-orientation
function proposed by Rietveld [1] as a special case (p1= 0). One may be forced to fix p1 at 0 or anotherconstant value between 0 and 1 because the standard
deviations of p1 and p2 often become extraordinarily
large on refinement of both p1 and p2.
Equation (4) is merely an empirical function, andneither p1 nor p2 has physical meaning. Neutron diffraction is highly recommended when dealing withvery orientation-prone substances because it gives
much smaller errors arising from preferred orientation
than X-ray diffraction [17].
8. Profile Shape Function
The profile shape function to approximate thepeak shape of each reflection can be represented by
the product of a symmetric profile shape function,g(∆θik), and an asymmetric correction, a(∆θik) [1]:G(∆θik) = g(∆θik)a(∆θik)(5)This asymmetry correction is applied to intensity datain the low 2θ region.
8.1 Symmetric Profile Shape Function

Symmetric profile shape functions used earlier in
Rietveld analysis for the X-ray case were Gaussian[1], Lorentzian [18], modified Lorentzian [19], and intermediate Lorentzian [20]. At present, better symmetric profile shape functions, i.e., pseudo-Voigt [21]Vol. 6 No. 1 1989

and Pearson VII functions [22], are widely used inplace of the above four functions [23].
The present program adopts a modified pseudoVoigt function, which is a linear combination of a
Gauss function and a Lorentz function with unequalpeak heights and full widths at half-maximum intensity (FWHM):2 2∆θ ik  





g ( ∆θ ik ) = C  γ exp −4 ln 2
 H k (G )  



 

(6)
−1
2



 2∆θ ik   
+ (1 − γ )1 + 4
 H k ( L)   


with
−1

 π 1 2π(1 − γ )H k ( L) 
C = 
 γH k (G ) + (7)2 4 ln 2 [

H k (G ) = U ( tan θ k − c s ) +V ( tan θ k − c s ) + W
2

H k ( L) = H k (G ) δ

]

12

(8)

(9)

In the above equations, C =normalizing constant, γ =fraction of the Gaussian component, Hk(G) = FWHM
of the Gaussian component, Hk(L) = FWHM of theLorentzian component, δ =ratio of H k(G) to H k(L), cs= 0 or 0.6. Equation (8) represents the dependence ofHk(G) on θk; U, V, and W are referred to as FWHM
parameters. When correlations among the FWHM
parameters are very high, cs should be set at 0.6 [23].
Equation (6) contains five refinable profile parameters: U, V, W, γ, and δ. With appropriate choices ofγ, this profile shape function can vary from Gaussian(γ= 1) to Lorentzian (γ = 0). Neutron data are usually
analyzed by fixing γ at 1 because the Gauss function
fits neutron diffraction profiles very well [1].
The difference between the conventional pseudoVoigt function [21] [23] and Eq. (6) is that δ is fixed at
one in the former whereas δ is a variable parameter inthe latter. The use of Eq. (6) leads to R factors a little
lower than those obtained with the conventional
pseudo-Voigt function. However, computational experience with refinements using Eq. (6) has shownthat the correlation coefficient between γ and δ isusually very high; one may sometimes be obliged tofix δ at an appropriate value.
8.2 Asymmetric Correction

In the present version, an empirical equation described by Rietveld [1] is used as a(∆θik):13

a(∆θik) = 1-( sin(∆θik)(2∆θik)2 )/tan(θk)(10)
a is the asymmetry parameter to be refined. As θk isincreased, a(∆θik) moves on toward 1.
8.3 Another Representation for the Dependence of
FWHM on Scattering Angle

Thompson, Cox, and Hastings [24] reported thatthe pseudo-Voigt function can be represented by only
Hk(G) and Hk(L) and that the variation of Hk(G) andHk(L) with Bragg angle can be approximated closely
by the functions Vtanθk and Xcosθk which respectively represent the contributions from instrumental
resolution and particle size broadening. Thus, theprofile shape function contains only two profile
parameters: V and X. Their formulation is worth
noting because it is soundly based on physical reality.
The effectiveness of this function in Rietveld analysis
will be tested in the near future.
9. Background Corrections
Background corrections are evaluated as follows.
(i) The background function, yib(c), is used which
is linear in six refinable parameters b0-b5:
j 2θ − θ max − θ min y ib ( c) = ∑ b j  i (11)j =0
 θ max − θ min
where 2θmax and 2θmin are maximum and minimum
diffraction angle, respectively.
(ii) Four pairs of smoothed values of operatorselected points in the pattern are fitted with a powerseries polynomial of degree 3, and the value of thispolynomial corresponding to the given value of 2θi is
calculated.
The use of Eq. (11) is preferable except for thecase dealing with samples showing very simple
diffraction patterns, because peaks overlap to a great
extent in the high 2θ region. In Eq. (11), 2θi is normalized between -1 and 1 to reduce the correlations
between b 0 - b5.
5

10. Summary of Refinable Parameters
Table 1 shows variable parameters contained inyi(c).
Lattice parameters are contained in θk in Eq. (1).
However, it is not lattice parameters but elements ofthe metric tensor for the reciprocal lattice (a*2, b* 2,
c*2, b*c*cosα*, c*a*cosβ*, a*b*cosγ*) that are
refined in the program [1]; these elements areconverted into lattice parameters upon printing ofparameter values.

14

Table 1 Parameters refined in Rietveld analysis.

• Global parameters
Parameter to correct for the zero-point shift of a goniometer
1) Zero-point error: Z
Parameters to calculate the background2) Background parameters b0, b1, b2, b3, b4, b5• Phase-dependent parameters
Parameters to adjust integrated intensities
3) Scale factor: s4) Preferred-orientation parameters: p1, p2
Profile parameters
5) FWHM parameters: U, V, W
6) Gaussian fraction: γ7) FWHM ratio: δ8) Asymmetery parameter: A
Parameters to determine peak positions
9) Lattice parameters: a, b, c, α, β, γCrystal-structure parameters
10) Fractional coordinates: xj, yj, zj11) Occupation factors: nj12) Overall thermal parameter: Q13) Isotropic thermal parameters: Bj
14) Anisotropic thermal parameters: β11j, β22j,
β 33j, β 12j, β 13j, β 23j
Magnetic-structure parameters
15) Magnetic moments in Bohr magneton: µj16) Angle between the spin direction and the unique axis*: Φ17) Angles between the spin direction and the three axes**:Φa, Φb, Φc* Tetragonal, hexagonal and rhombohedral systems** Orthorhombic systems

11. Constrained Minimization Procedures
11.1 Linear Equality Constraints

In the RIETAN program, simple linear equality
constraints are explicitly solved for selected variables, and those variables eliminated from the problem as independent variables. Linear constraints are
input in just the same form as assignment state-ments
in FORTRAN, e.g., ‘A(60) = A(9)’, ‘A(29) = 1.0A(24)’, and ‘A(41)=O.95-0.5*A(31)’. In these expressions, array A stores both variable and fixedparameters contained in yi(c). When dealing with
mixtures of two or more phases, common profile
parameters are usually assigned to each phase byintroducing linear equality constraints.

The Rigaku Journal

11.2 Refinement Based on Information about the
Crystal Structure

Powder diffraction patterns contain far poorer
information than single-crystal intensity data mainly
owing to the overlap of equivalent and neighboringreflections. Accordingly, a priori information on the
crystal structure should be introduced into the program as constraints when reflections overlap very
extensively or when a number of structure parameters
have to be refined because many atoms are containedin an asymmetric unit. Most programs for Rietveld
analysis are able to apply only equality constraints
such as those required to keep certain bond lengths orbond angles fixed or equal [25]. True bond lengths or
angles are not exactly equal to expected ones.
Therefore, the introduction of inequality constraints[26] is preferred, permitting the imposition of ranges
expected for interatomic distances and bond angles aswell as thermal parameters as boundary conditions.

11.3 Nonlinear Constraints

Constrained nonlinear programming problems
are solved by an exterior penalty function method[27]. Rietveld analysis under nonlinear constraints isformally stated as
Minimize: f(x)= ∑ w i [yi(o)-y i(c)]2(12)
isubject to p linear and/or nonlinear inequality
constraints
gn(x) ≥ 0 n= 1, 2, . . . ,P
(13)
and q linear and/or nonlinear equality constraints
hn(x) = 0 n=1,2, . . .,Q(14)where x is the vector of variable parameters, and w i (=1/yi(o)) is the weighting based on counting statistics.
The exterior penalty function method transforms
a constrained optimization problem into a sequenceof unconstrained optimizations for K=0, 1, 2, given
byMinimize: F(x, t(K) = f(x) +Fig. 2 Function subprogran CON written for the introduction of 12 constraints in the Rietveld anaylsis of ammonioleucite.

Vol. 6 No. 1 1989

15

p
q
K 

t ( )  ∑ H (g n ( x ))g n2 ( x ) + ∑ hn2 ( x ) (15)n =1
 n =1
where t (K) is a strictly increasing sequence of positive
numbers, and H is the Heaviside operator such that
H(a) = 0 for a ≥ 0 and H(a) = 1 for a < 0. The second
and third terms in Eq. (15) prevent the x vector fromgetting too far away from the feasible region. A
nonfeasible sequence of points generated by thisalgorithm may yield a feasible solution in the limit,
i.e., t(K) = ∞. Because Eq. (15) has a sum-of-squares
form, conventional algorithms for nonlinear leastsquares computations, which will be described inchapter 12, may be employed without any modifications.
The algorithm of the exterior penalty function
method consists of the following four steps:
1) Set K at 0. Give the initial values of x and t(0).
2) Refine a set of parameters x which minimize F(x,T(K)).
3) If the second and third terms in Eq. (15) are reduced
to nil, stop the calculation since the current valuesof x are the solution.
4) Add 1 to K Increase t(K) and return to step 2).
As Fig. 2 shows, the constraints are introduced
into the program by means of a separately written
function subprogram named CON, the purpose ofwhich is to evaluate H(gn(x))gn(x) and hn(x) whenprovided with the x vector. Those partial derivativesof the functionals of constraints with respect toparameters which are used in least-squares calculations are approximated by centered differences not
only to save preparation times to formulate analyticalderivatives but also to avoid human errors.
The above method has been successfully applied
to the refinement of the structure of yttriumcontaining α-sialon [28]. In that work, X-ray powder
data were supplemented with eight inequality constraints that restrict the (Si, Al)-(N, O) bond lengths
within a reasonable range: 0. 175 ± 0.005 nm.

12. Nonlinear Least-Squares Procedures
Almost all computer codes for Rietveld
refinement employ the Gauss-Newton algorithm tofind parameters which minimize the weighted sumof-squares of the residuals apart from the progressive
programs of Baerlocher [26], and Howard and Snyder[29]. However, when applied to Rietveld analysis
without any improvements, this method suffers
disadvantages that the range of convergence is notvery great and that the refinements often converge tolocal minima [4]. Since none of the algorithms has
16

proved to be so superior that it can be classified as auniversal panacea for nonlinear least-squares solutions, it is impolitic to persist in only one method ofleast squares.
Three different methods are available fornonlinear least-squares fitting: Gauss-Newton, modified Marquardt, and conjugate direction methods. Allof them are designed to give stable convergence. The
Gauss-Newton and Marquardt methods use derivatives of yi(c) with respect to x, whereas the conjugate
direction method is based on a direct-search algorithm. Since the values of initial parameters in the,input file can be replaced with those of final parameters after the refinement, it is very easy to change themethod of least squares one after another and continue the calculation from the last point of the previous refinement.
In each of these methods, variable parameters in
an input file can be updated after they have beenrefined. Thus, one can carry out Rietveld analysis
again with parameter values obtained in the previous
refinement.
12.1 Gauss-Newton Method

This algorithm calculates changes in parametersat each iterative step, ∆x, from the normal equation:
M∆X=N
(16)
with
M jk = ∑ wi
i

+t

(K )

∂yi ( c) ∂y i ( c)∂ x j ∂ xk

 p
∂g n ∂g n q∂h n ∂ h n ()
Hg+ ∑ H( hn )
∑
n

∂
x
∂
x
∂ x j ∂x k 
n=1
 n=1
j
k

N j = ∑ wi [ y i ( o ) − y i ( c)]i−t(K ) 

(17)∂y i ( c)
∂x j

q
∂g n
∂h 
+ ∑ hn n  ∑ H ( g n )g n
∂x j n =1 ∂x j  n =1(18)p

In Eqs. (17) and (18), g n(x) and h n(x) are, respectively,
abbreviated as gn and hn. The coefficient matrix M
and the column matrixes ∆x and N are calculated indouble precision owing to the large number of data
points in Rietveld analysis and considerable numerical errors in the least-squares method with the normal
equation.
Although ∆x is evaluated from M-1N in moststructure-refinement programs, there is little torecommend such a technique because of long computation time and low precision. In RIETAN, theCholeski decomposition of the symmetric, positive

The Rigaku Journaldefinite matrix, M, and back-solution for consistentsystems of linear equations are performed. In cases
where M is non-positive definite, the Marquardt
method described in the next section should be used.
A new set of x, x', is readily obtained byx'= x + d∆x
(19)
with
d=2 -n n = 0, 1, 2, 3, 4
The initial value of the damping factor, d, is set at 1 (n= 0). If F(x, t(K)) does not decrease with x', d isdecreased, and x' is calculated again with Eq. (19).
The value of d is automatically adjusted according tojust the same rule adopted in a SALS program [30].
12.2 Modified Marquardt's Method

This method also uses Eqs. (17) and (18) but adds
λ-diag(M) (λ: Marquardt parameter, diag: diagonal
matrix) to M to stabilize the convergence to theminimum [31]:(M+ λ - diag(M))∆x = N(20)
Then, ∆x tends towards the steepest descents
direction as λ becomes larger, while the GaussNewton solution is obtained when λ becomes negligible. The value of λ is adjusted during a series of
iterations by the same technique as described by
Fletcher [32]. Even if the coefficient matrix M is not
positive definite, it can be made computationally
positive definite by choosing the Marquardt parameter
large enough. Modified Marquardt's method is very
effective for dealing with highly nonlinear model

functions, yi(c), or problems in which starting values
for parameters differ markedly from the true ones.
12.3 Conjugate Direction Method

The conjugate direction method [33] is one of themost efficient algorithms to minimize objective functions without calculating derivatives. The minimumof F(x, t(K)) is located by successive unidimensional
searches from an initial point along a set of conjugate
directions generated by the procedure. A combinationof the DSC and Powell algorithms [34] was selected
as a unidimensional optimization routine.
Since the directions for minimization are determined solely from successive evaluations of the objective function, F(x, t(K) ), this procedure is muchslower than the two least-squares methods with derivatives but capable of solving ill-conditioned problems in which very high correlations exist between
parameters. It is mainly used in the late stages ofrefinement to test the prospect of a local minimum
being the global minimum or escape from a local
minimum by using sufficiently large step sizes of line
searches. On the other hand, the Gauss-Newton andMarquardt algorithms can check the convergence tothe global minimum only by using different starting
vectors.
12.4 Incremental Refinements

One usually proceeds in steps in Rietveld
analysis, first refining only one or two parameters andthen gradually letting more and more of theparameters be adjusted in the successive least-squares
refinement cycles [4]. RIETAN requires only a single

Fig. 3 X-ray powder diffraction pattern simulated for the orthorhobmic form of Ba2YCu3O7-y .

Vol. 6 No. 1 198917

job to execute such incremental refinements; that is,variable parameters in each cycle can be appointed
when using the methods of least squares with
derivatives. Troublesome repetition of batch jobs is,therefore, unnecessary in most Rietveld refinements.
Usually, linear parameters, that is, the backgroundparameters and the scale factor are refined in the firstcycle, profile parameters in the second cycle, latticepa-rameters in the third cycle, and subsequently all the.parameters simultaneously. Thus far, almost all
Rietveld refinements have been carried out routinely
with this procedure.
Automatic refinements are also possible in which
the parameters obtained by the successive refinements
described above are further adjusted by the conjugate
direction method to ensure that there are no lower
minima in the vicinity of the the one found by the
initial refinement.
12.5 Standard Deviations

The standard deviations of parameters are calculated from diagonal elements in the inverse of thecoefficient matrix M according to the procedure proposed by Scott [35]. M is neither inverted in each
iteration in the two methods using derivatives nor
calculated in Powell's method. For this reason, thestandard deviations are calculated by selecting theGauss-Newton method and specifying the number of
iterations as zero after the refinement has finished.

The standard deviations of lattice parameters arecalculated from those of the elements of the metrictensor (cf chapter 10) with formulae given by Kelsey[36].
13. Graphic Data Processing
One can not only plot the observed and calculated
intensities after refinement by the least-squares methods (Fig. 1) but also simulate powder diffraction
patterns of pure substances and mixtures from usersupplied profile, lattice, and structure parameters(Fig. 3). These figures are drawn using an incremental
plotter, a laser beam printer, or a graphic display ininteractive mode. Theoretical diffraction patterns
serve for various semi-qualitative analyses such as
examination of the validity of structural models,order-disorder distribution, degree of preferred orientation, characterization of polytype structures, andestimate of occupation factors in solid solutions.
Since CalComp-compatible graphic subroutine packages are used for portability, RIETAN must be run onother computers without difficulty.
14. Connection with Other Programs
After the refinement has been finished, three files
can be created which store (i) y i(o), y i(c), and 2θk, (ii)final crystal-structure parameters, their standard deviations, final lattice parameters, symmetry operations etc., and (iii) hkl, F’o’, and Fc. Single quotation

Fig. 4 Rietveld refinement profile for TOF neutron diffraction data of La1.9Ca1.1Cu2O6 . Q(=2π/d) is plotted as abscissa and the net
intensity as ordinate. Backgournd was fit as part of the refinement but has been subtracted before plotting.

18

The Rigaku Journalmarks in F’o’ are needed because it is estimatedindirectly from summation of contributions of thepeak to net observed intensities [1]. The above three
files respectively enable (i) drawing Rietveld
refinement patterns (Fig. 1), (ii) calculation of interatomic distances and bond angles with an ORFFE
program [37], calculation of Madelung energies by
Ewald's method, and drawing crystal structureillustrations with an ORTEP 11 program [38], and (iii)
Fourier and D syntheses with F’o’ and Fc data. Fourier
and D syntheses, followed by drawing contour maps,
are very useful for improving incomplete structuralmodels.

One can enjoy the benefits of the common and userfriendly software for the Rietveld analysis of resultsfrom the four categories of experimental technique.
Some time later, the two versions of RIETAN
will be integrated into one program so that structure
parameters can be refined simultaneously with theuse of both X-ray and neutron diffraction data [42].
The different diffraction methods allow the collectionof several, in many respects complementary, data setsfrom one kind of sample. Therefore, Rietveld refinements with the combined intensity data will be very
helpful to obtain detailed information about the structure parameters of whole kinds of atoms.

15. Further Developments of RIETANThe Rietveld method has been successively
applied to time-of-flight (TOF) neutron powder diffraction data with an adaptation of the above program(Fig. 4). These intensity data were measured at a fixed
scattering angle on a High-Resolution Powder
diffractometer, HRP [39], at the pulsed neutron source(KENS) at the National Laboratory for High Energy
Physics (KEK). This TOF neutron diffraction version
retains all the features of the angle-dispersive one. A
new profile shape function optimized for a cold
neutron source (solid methane at 20 K) at KENS was
implemented in it. Details on this version have beenreported elsewhere [40]. Recently, it has been applied
successfully to the structure refinement of a series ofhigh-Tc superconducting oxides [41].
The program written for the fixed-wavelength
case has been modified so as to enable the Rietveld
analysis of synchrotron X-ray data obtained from a
new dedicated powder diffractometer set up at the
Photon Factory at KEK. The flat-plate geometry is
unsuitable for collecting intensity data because of
marked preferred orientation caused by highly parallel
X rays. The Debye-Scherrer geometry with rotation ofglass capillaries containing samples (some-times
mixed with amorphous silica) is usually adopted toreduce the preferred-orientation effect as much as
possible [24].
Thus, RIETAN is now applicable to therefinement of crystal structures from powder dataobtained by four different methods [3]:Acknowledgment
The author wishes to thank Dr. A. Yamamoto for
helpful advice on computing methods in crystallography. X-Ray powder diffraction data used to plot
Fig. 2 were kindly provided by Professor R. A.
Young.

characteristic X-ray source
Fixed-wavelength methods withsynchrotron X-ray source
reactor neutron source

Fixed-angle method with pulsed neutron source
Vol. 6 No. 1 1989

References[1] H. M. Rietveld: J. Appl. Crystallogr. 2 (1969) 65.
[2] A. K. Cheetham and J. C. Taylor: J. Solid State Chem. 21(1977) 253.
[3] A. Albinati and B. T. M. Willis: J. Appl. Crystallogr. 15(1982) 361.
[4] R. A. Young: Accuracy in Powder Diffraction, NBS Spec.
Publ. (U.S.) 567, National Bureau of Standards,
Washington (1980), p. 143.
[5] F. Izumi: X-Sen Bunseki no Shinpo (Adv. X-Ray Chem.
Anal. Jpn., in Japanese) No. 14, Agne, Tokyo (1982), p.
43.
[6] F. Izumi: Nippon Kessho Gakkai Shi (J. Crystallogr. Soc.
Jpn., in Japanese) 27 (1985) 23.
[7] D. B. Wiles and R. A. Young: J. Appl. Crystallogr. 14(1981) 149.
[8] International Tablesfor X-Ray Crystallography, Vol. 1,
Kynoch Press, Birmingham (1969).
[9] International Tables for Crystallography, Vol. A, D.
Reidet, Dordrecht (1983).
[10] International Tables for X-Ray Crystallography, Vol. IV,
Kynoch Press, Birmingham (1974).
[11] V. F. Sears: Neutron Scattering (Part A), Methods ofExperimental Physics, Vol. 23, ed. by K. Sk6ld and D. L.
Price, Academic Press, New York (1986), p. 521.
[12] M. J. Cooper, K. D. Rouse and R. Sayers: A Computer
Program for the Generation of Unique Reflection Sets
for All Space Groups, AERE-R8695, AERE Harwell(1977).
[13] J. S. Rollet: Computing Methods in Crystallography, ed.
by J. S. Rollet, Pergamon Press, Oxford (1965), p. 25.
[14] G. Shirane: Acta Crystallogr. 12 (1959) 282.
[15] Y. Sasa and M. Uda: J. Solid State Chem. 18 (1976) 63.

19

[16] H. Toraya and F. Marumo: Mineral. J. 10 (1 981) 21 1.
[17] A. W. Hewat: Accuracy in Powder Diffraction, NBS
Spec. Publ. (U.S.) 567, National Bureau of Standards,
Washington (I 980), P. III.
[18] R. A. Young, P. E. Mackie and R. B. Von Dreele: J.
Appi. Crystallogr. 10 (1977) 262.
[19] E. J. Sonneveld and J. W. Visser: J. Appl. Crystallogr. 8(1975) 1.
[20] C. P. Khattak and D. E. Cox: J. Appl. Crystallogr. 10(1977) 405.
[21] G. K. Wertheim, M. A. Butler, K. W. West and D. N. E.
Buchanan: Rev. Sci. Instrum. 45 (1974) 1369.

[30] T. Nakagawa and Y. Oyanagi: Analysis of Experimental
Data with Least-Squares Methods (in Japanese), Tokyo
Daigaku Shuppankai, Tokyo (1982), p. 98.
[31] D. W. Marquardt: J. Soc. Ind. Appl. Math. 11 (1963)
431.
[32] R. Fletcher: A Modified Marquardt Subroutine for NonLinear Least Squares, AERE-R6799, AERE Harwell(1971).
[33] M. J. D. Powell: Computer J. 7 (1964) 155.
[34] D. M. Himmelblau: Applied Nonlinear Programming,
McGraw-Hill, New York (1972), p. 42.
[35] H. G. Scott: J. Appl. Crystallogr. 16 (1983) 159.

[22] M. M. Hall, Jr., V. G. Veeraraghavan, H. Rubin and P. G.
Winchell: J. Appl. Crystallogr. 10 (1977) 66.
[23] R. A. Young and D. B. Wiles: J. Appl. Crystallogr. 15(1982) 430.

[36] C. H. Kelsey: Mineral. Mag. 33 (1964) 809.
[37] W. R. Busing, K. 0. Martin and H. A. Levy: ORFFE,
Report ORNL-TM-306, Oak Ridge National Laboratory,
Tennessee (1964).

[24] P. Thompson, D. E. Cox and J. B. Hastings: J. Appl.
Crystallogr. 20 (1986) 79.

[38] C. K. Johnson: ORTEP II, Report ORNL-5138, Oak
Ridge National Laboratory, Tennessee (1976).

[25] G. S. Pawley: J. Appl. Crystallogr. 13 (1980) 630.
[26] Ch. Baerlocher: Acta Crystallogr., Sect. A 40 (1984) C368.

[39] N. Watanabe, H. Asano, H. lwasa, S. Satoh, H. Murata,
K. Karahashi, S. Tomiyoshi, F. Izumi and K. Inoue: Jpn.
J. Appl. Phys. 26 (1987) 1164.
[40] F. Izumi, H. Asano, H. Murata and N. Watanabe: J.
Appl. Crystallogr. 20 (1987) 411.

[27] W. I. Zangwill: Management Sci. 13 (1967) 344.
[28] F. Izumi, M. Mitomo and Y. Bando: J. Mater. Sci. 19(1984) 3115.
[29] S. A. Howard and R. L. Snyder: Acta Crystallogr., Sect.
A 40 (1984) C-369.

20[41] F. Izumi and I. Nakai: Nippon Kessho Gakkai Shi (J.
Crystallogr. Soc. Jpn., in Japanese), 29 (1987) 365.
[42] J. K. Maichle, J. lhringer and W. Prandl: J. Appl.
Crystallogr. 21 (1988) 22.

The Rigaku Journal