Crystal structure determination by dual-space iterative algorithms

After the first structure solutions which used symmetry arguments and trial-and-error methods, the Patterson method became the first systematically used approach to structure solution. When the statistical relationships between the reflection intensities and their phases were discovered by Cochran, Sayre, Karle, Hauptman and many others in the 1950s and 1960s, a rich field of direct methods was developed [see e.g. Giacovazzo (1998) for an overview of the subject]. The continuous development and growing power of direct methods made them a leading tool for ab initio structure solution, and their dominance seemed to be incontestable.

The development of crystallographic methods received an important new impulse with the advent of powerful desktop computers. Suddenly, computationally expensive approaches became available to everybody, and methods could be developed that make heavy use of computationally demanding techniques, such as Fourier transformation. The fruits of this revolution are numerous. Direct methods were combined with density modifications in direct space to produce the `Shake and Bake' (SnB) method (Weeks et al., 1993). A flavour of direct methods based on Patterson-function arguments was developed (Rius, 1993) and later transformed into an algorithm cycling between direct and Fourier space (Rius et al., 2007; Rius & Frontera, 2008; Rius, 2012). The `revenge of the Patterson function' was announced (Burla et al., 2004, 2006), combining the superposition minimum function method (Buerger, 1959) with new analysis and computer power, and applying the method to ab initio phasing of macromolecular crystals. Yet another algorithm, called VLD (Burla et al., 2010, 2011a,b, 2012b), is based on an iterative application of difference-Fourier synthesis and real-space density modification techniques.

Next to all these algorithms, a class of so-called dual-space algorithms emerged. This name might be somewhat confusing, because many if not all of the algorithms mentioned in the previous paragraph are in certain sense dual-space. However, in the dual-space algorithms sensu stricto, neither of the two spaces plays a dominant role – neither the operation in direct space nor in the Fourier space is, even in principle, capable of solving the structure alone.

The iterative dual-space phasing algorithms have gained considerable interest in the crystallographic community over the last few years. Among them the best known is the charge flipping algorithm (CFA, Oszlányi & Sütő, 2004). While not the first published algorithm of this kind, it sparked considerable interest in iterative dual-space methods for structure solution, and has marked the beginning of a broad interest in and development of this field in crystallography. Another algorithm that gained considerable popularity thanks to its implementation in the computer program SHELXT is the algorithm called intrinsic phasing (IPA; Sheldrick, 2015).

This chapter describes the dual-space iterative algorithms, summarizes their development and gives an overview of applications of these algorithms to crystallographic problems. It is an updated and extended version of the review by Palatinus (2013), which summarized the field up to then. The text is organized as follows. First, a general overview of dual-space algorithms in phase retrieval is presented with a focus on algorithms relevant for crystallography. Then a detailed description of the charge flipping algorithm and its variants is provided, followed by sections on two important special topics: symmetry and missing data. Finally, applications of the algorithms are described in sections devoted to general structure solution, modulated structures, powder diffraction data and macromolecules. This chapter deals only with algorithms used for ab initio phasing from data with atomic resolution. An important part of the application of dual-space algorithms concerns the phasing of macromolecular structures from data that do not reach atomic resolution. Although the basic principles are similar, the low resolution of macromolecular data poses specific challenges and requires additional information to be used in the phasing process, like the assumption of an existence of molecular envelope or non-crystallographic symmetry. This part is not covered in this text and the reader is referred to the review by Millane & Lo (2013) or the chapter of International Tables for Crystallography Volume C by Millane (2023).

The problem of phase retrieval is omnipresent in various fields of physics and engineering (Millane, 1990). The problem is usually formulated in the general space of square-integrable functions, but for our purposes let us limit the definition to a discretely sampled signal in Euclidean space Eⁿ of dimension n (n = 3 for a typical crystal structure). Let ρ = {ρ_i, i = 1…N_p} be a (generally complex-valued) function sampled on a discrete n-dimensional grid comprising N_p pixels. Let be its (discrete) Fourier transform. Let M = {h_j, j = 1…N_h} be the independent set of vectors h_j, for which Fourier magnitudes are known from experiment. For the most general case of phase retrieval of a complex object, all vectors h_j are independent. If ρ is a real-valued function, then the Fourier coefficients F(h) and F(−h) must be complex conjugate, and only half of all h_j are independent. Further reduction of the number of independent vectors h_j may occur if crystal symmetry is applied.

With these definitions, the phase retrieval problem can be formulated as: find ρ or an approximation to ρ, given the set of known Fourier magnitudes (magnitude constraint). This problem does not have a unique solution, as the constraint ratio Ω, i.e. the ratio of the amount of information available and amount of information needed, equals 1/2 (Elser & Millane, 2008; Millane & Arnal, 2015; Elser et al., 2018). Therefore, additional information about ρ is necessary. This additional information can be the support (i.e. a subset of ρ_i is assumed to be zero), positivity (ρ_i > 0 for all i), atomicity (the signal in ρ is composed of a set of discrete peaks) or any other piece of information. Since this additional information is in all practical applications defined in the direct space of ρ, not in its Fourier space, it will be denoted as the direct-space constraint. The set of all ρ that fulfil the constraint is called the constraint set. Let us denote by S the set of all ρ matching this direct-space constraint, and by R the set of all ρ such that , i.e. matching the magnitude constraint. Then the phase retrieval problem can be simply stated as: find some ρ from S ∩ R. If such ρ exists, the problem is called consistent. If the intersection of S and R is empty, the problem does not have a solution and is called inconsistent. In such a case, it may still be useful to search for ρ such that the sum of its distances to the nearest points in S and R is minimal.

With this formulation it becomes obvious that the phase retrieval problem is a specific case of the general constraint satisfaction problem, where the two (or more) constraints need not be specifically the magnitude and direct-space constraint. A special, mathematically more easily tractable case is the convex feasibility problem, where all the constraints are convex. A constraint set A is convex if for any two elements of the constraint set the following statement is valid: The convexity of the constraints allows important conclusions to be made about the properties and convergence of algorithms proposed for the solution of the convex feasibility problem. Algorithms exist that always converge to a solution of this problem. Unfortunately, the magnitude constraint central to the phase retrieval problem is obviously non-convex, and the results of the convex optimization theory cannot be carried over to the phase retrieval problem. Nevertheless, it is useful to compare the algorithms developed in these two frameworks. Analysis of the convex versions of the algorithms gives valuable insight into the relationships between various algorithms proposed in phase retrieval.

The formulation of constraints plays an important role in the specification of the algorithms. Let us summarize the specific forms of the magnitude and direct-space constraints encountered in crystallographic phase retrieval. The basic form of the magnitude constraint set was defined above as the set of all ρ such that for all h_j ∈ M. The recipe to transform an arbitrary set of values of ρ into a set that belongs to the constraint set is called a projection operator or projector, and the operation itself is called projection. In order to be well defined, the projection is required to be distance-minimizing, i.e. the change of ρ needed to bring it onto the constraint set must be the smallest possible. Operations that do not satisfy this condition, but still transform ρ to a point in the constraint set are called pseudoprojections.

The basic projector onto the magnitude constraint set can be defined as follows: ρ^new is then obtained as an inverse Fourier transform of . In other words, the new ρ is obtained from the old by replacing the known Fourier magnitudes with the observed ones, keeping all phases and the unknown Fourier magnitudes intact.¹ However, alternative variants of this basic scheme have been proposed and used. Elser (2003b) proposed placing an upper bound on the magnitude of unknown Fourier coefficients:

In crystallography, suitable can be conveniently estimated from the Wilson statistics as the expected mean Fourier amplitude for structure factors with the given |h_j|. For normalized Fourier magnitudes (E values), F^bound = 1 for all j. If c goes to infinity, this operator transforms to (2). Another special instance of this operator is the case c = 0:

Yet another variant has been devised for practical applications. The Fourier magnitudes |F^o| are known only for associated scattering vectors up to a certain length . The sphere of known |F^o| is denoted the resolution sphere. In practical applications, the set M of scattering vectors with known Fourier magnitudes rarely covers all vectors in the resolution sphere. For example, the forward scattering intensity, F(000), is never measured. Such missing data inside the resolution sphere require a treatment different from unmeasured data at high resolution. For this purpose, a combination of projector (3) for data inside the resolution sphere and (4) outside the resolution sphere is useful:

Again, the constant c can be set infinite, in which case the second condition is never met and the rule reduces to leaving all unmeasured magnitudes inside the resolution sphere unchanged, and setting everything outside the resolution sphere to zero.

Equation (5) has a very important special case, namely c = ∞ and so small that the only coefficient with is the coefficient F(0). This case corresponds essentially to equation (4), but instead of constraining F(0) to zero, it is left unchanged by the projector. This form was used in the original article on the CFA (Oszlányi & Sütő, 2004):

Among direct-space constraints, the most studied constraint in phase retrieval is the support constraint. This constraint can be applied if some part of ρ is known or assumed to be zero. This is a relevant constraint in single-particle imaging and in macromolecular crystallography. However, in most crystallographic applications the distribution of scattering density in the unit cell is usually unknown, and no support constraint can be defined. In the absence of a known support, the positivity of the electron density can still be exploited, and the positivity constraint can be conveniently defined. The corresponding projector is very simple:

Such an operation has been frequently used in macromolecular crystallography as a part of phase extension and refinement procedures (solvent flattening; Wang, 1985). For ab initio crystal-structure determination, however, if was found that a more aggressive density modification technique is needed:

Introduction of the free parameter δ gives the algorithms more freedom to find a balance between the perturbing and stabilizing effect of the operation. Such an operation is the basis of some of the electron-density modification (EDM) techniques used in direct methods (e.g. Giacovazzo & Siliqi, 1997). In the context of ab initio structure solution by dual-space methods, it was first proposed by Shiono & Woolfson (1992), and it is in the background of the charge-flipping operation. As noted e.g. in Oszlányi & Sütő (2008), this operation is not distance-minimizing, i.e. it is not a projection in the strict sense, but a pseudoprojection. It is important to note that while projections onto convex constraint sets are uniquely defined, the number of possible pseudoprojections is in general infinite. The choice of a particular pseudoprojection instead of a projection is therefore to a large extent arbitrary. Its choice is, typically, not based on any fundamental theory, but on the practical observation that that particular choice performs well in practice. For the sake of simplicity, we will not make an explicit distinction between true distance-minimizing projections and pseudoprojections in this text, where this is not explicitly needed.

Operation (8) lends itself to a modification, which is not possible with (7), namelyThis constraint does not impose positivity on ρ, but only eliminates low-density regions. It can be understood as a `dynamical support' constraint, where, unlike in the standard support constraint, the support is newly identified at every iteration cycle as the region with small density. This operation is at the basis of the band-flipping approach (Oszlányi & Sütő, 2007). An asymmetric interval of the eliminated region is also possible and useful if the negative and positive signals in the density are of unequal strength (Whitfield & Coelho, 2016):

Another direct-space constraint used in crystallography is the atomicity constraint. This constraint is less easily expressed by a simple recipe, but, in rough terms, the corresponding projector consists of selecting a prescribed number of peaks in ρ, and setting to zero all pixels outside these peaks (Elser, 2003b; Feng, 2012). A very special type of the atomicity constraint is used in the IPA. The corresponding operation consists of locating the peaks in ρ, placing a Gaussian with unit volume at the positions of the peaks, and multiplying ρ by the mask formed by these Gaussians (Sheldrick, 2015). Although this operation is not a projection, as the resulting density does not have a specific property that would define it, it is closely related to an atomicity projection because the result of the operation is a set of discrete, sharp peaks, and all negative densities are set to zero.

Whatever the exact constraint and the recipe for the projector, it can be symbolically denoted as P, and the transformation of the image then as ρ^new = Pρ^old, where P is an operator acting on ρ. Such operators can be combined to yield more complicated transformations. A particularly important combination of projections is a so-called overprojection,where I is the identity operator. Geometrically, such overprojection means that the shift from ρ to Pρ is continued in the same direction by a fraction γ of the original distance from ρ to Pρ. The special case of γ = 1 is called a reflector, and will be denoted R without an explicit superscript:

We are now equipped to review the different algorithms suggested in the literature for phase retrieval, and specifically for crystal structure solution. In some works, the iterative phase retrieval algorithms are described as explicit schemes for pixelwise obtaining ρ of cycle n + 1 from ρ of cycle n. Such recipes do not explicitly separate the combination of projectors acting on ρ from the exact definition of the projector, and sometimes even cannot be expressed in the form of a combination of projectors acting on ρ. Another approach to describe the algorithms is to define the iteration scheme in terms of operators acting on ρ⁽ⁿ⁾ to obtain ρ⁽ⁿ⁺¹⁾, and define the exact form of the operators separately. Where possible, we will adopt this second approach, and we will point out cases where this approach leads to difficulties. Explicit flowcharts of the most important algorithms are then summarized in Fig. 1.

Figure 1 Explicit schemes of the most important dual-space algorithms. Projections (6) and (8) were used for PM and PD, respectively. τ, ω, φ and ψ are intermediate images, , , and are their Fourier transforms, respectively. The schemes correspond to equations (14) (CFA), (18) (HIO), (23) (AAR) and (20) (DM).

The basic algorithm is the alternating application of the magnitude and direct-space projections. Expressed as an iteration scheme, this algorithm can be written as Here P_M denotes the magnitude projector, and P_D the direct-space projector. As for all algorithms that will be presented in this section, the iteration typically starts from a random starting image, but that is not strictly necessary, and starting from a non-random starting point changes these algorithms from phase retrieval to phase refinement or phase extension algorithms. This simple algorithm is known in the phase retrieval community as the Gerchberg–Saxton (Gerchberg & Saxton, 1972) or error-reduction algorithm (Fienup, 1982), and as the POCS (projections onto convex sets) for convex feasibility problems (Censor & Zenios, 1997). In crystallography, this algorithm was used in conjunction with projection (8) under the name low-density elimination (LDE; Shiono & Woolfson, 1992). This method was developed as a phase extension method for macromolecular crystallography, but the authors added a short comment stating that one of the test structures could be solved ab initio from random phases. This seems to be the first published record of a crystal structure solved ab initio by dual-space methods, although this possibility was considered much earlier, for example in the ground-breaking paper by Fienup (1982). A detailed account of the performance of LDE in ab initio solution is given in Matsugaki & Shiono (2001).

The Gerchberg–Saxton/error-reduction/LDE algorithm is known to be prone to stagnation at false solutions. One way of reducing the risk of stagnation is to replace the projectors by reflectors. Replacing the direct-space projector P_D by the corresponding reflector yields this algorithm: This is the iteration scheme of the basic CFA (Oszlányi & Sütő, 2004), with R_D being the reflector of (8) and P_M being the magnitude projection (6). As noted by Wu et al. (2004b), this algorithm is a special case of Fienup's output–output algorithm [Fienup (1982), equation (43) with β = 2 and with γ being the set of all pixels, where ρ_i < δ].

A symmetric counterpart of this algorithm is the following scheme: Here the reflector is applied in reciprocal space. This scheme is used in the IPA (Sheldrick, 2015). However, in the presented algorithm the recommended (and default) value of γ in R_M is 2 and not 1. The phase retrieval algorithm of Feng (2012) also resembles very much this type of algorithm, although the operator in Fourier space uses a special version of the Fourier coefficients, yielding the operatorThe structure of this operator resembles a reflector, but it is not a reflector in a strict sense.

Another logical extension of the iteration scheme is to replace both projectors by reflectors: It turns out that this scheme is difficult to use because of its instability. The perturbation induced by the alternating reflectors is too strong, and the solutions are not stable. This scheme has been made to work only by replacing the magnitude reflector by a special `partial reflector' [see Section 3 and equation (29)].

In their pioneering work, Fienup (1982) proposed a set of even more complicated iterations schemes, of which the hybrid input–output (HIO) algorithm proved to be the most successful. The hybrid input–output algorithm was defined as an explicit recipe for pixelwise obtaining ρ^new from ρ^old. Adapted to our notation, this scheme readswhere β is a free parameter of the algorithm, and (P_Dρ⁽ⁿ⁾)_i means the ith pixel of the image P_Dρ⁽ⁿ⁾. This algorithm cannot be expressed in a form of an operator acting on ρ⁽ⁿ⁾ (Bauschke et al., 2003). Only if the direct-space constraint is the support constraint (or another constraint, for which the projector is a linear operator) can the HIO algorithm be expressed as a fixed-point operator (Bauschke et al., 2002, 2003): If the second form of the iteration scheme (19) is combined with the positivity constraints and not the support constraint, yet another algorithm, called hybrid projection reflection (HPR), is obtained (Bauschke et al., 2003). The HIO algorithm, or elements thereof, has been used in crystallographic phase retrieval schemes (Wu et al., 2006; Lei, 2007).

Another algorithm that bears a strong relationship to the HIO algorithm was proposed by Elser (2003b) and named the difference map (DM). It is a three-parameter algorithm defined by the following scheme:

In the original work (Elser, 2003b), the optimal values of parameters γ_M and γ_D were estimated to be equal to β⁻¹ and −β⁻¹, respectively. In subsequent work (Elser, 2003c), a slightly different choice of γ_D was suggested. It can be easily shown (Elser, 2003b) that, for the case of support constraint only, the HIO algorithm is a special case of the DM with γ_M = β⁻¹ and γ_D = −1. This equivalence, however, does not hold for the positivity or atomicity constraint. The difference map was demonstrated to work for structure solution (Elser, 2003a). The specific form of the magnitude constraint was that of equation (3), and the direct-space constraint was the atomicity. For some reason this algorithm has never become very popular in crystallography, while it has found more applications in the phase retrieval of non-periodic objects.

When the special value of β = 1 is used in the second equality of equation (19), we obtain a particularly simple expression:

This algorithm was first proposed and analysed by Douglas & Rachford (1956) for the solution of differential equations, and adapted for convex sets by Lions & Mercier (1979). In the phase retrieval context it was suggested and analysed by Bauschke et al. (2004) under the name averaged alternating reflections (AAR). The AAR algorithm has the interesting property that, under certain circumstances, it is an important special case of both the HIO algorithm and the difference map. More specifically, assuming P_D is a linear operator, the AAR algorithm is the HIO algorithm with β = 1, and the difference map algorithm with β = 1, γ_D = −1 and γ_M = 1. Moreover, if, in addition to P_D, P_M is also assumed to be linear (keep in mind that this is not the case for the magnitude projection), Elser's difference map with the recommended parameters γ_M = β⁻¹ and γ_D = −β⁻¹ becomes a weighted average of two symmetric versions of AAR:

These relationships cannot be carried over directly to phase retrieval, where the magnitude constraint and often also the direct-space constraints are not linear. Nevertheless, they indicate that all these algorithms have a closely related structure.

Scheme (21) has a closely related symmetric counterpart, also obtained from (22) with β = −1: The two schemes are very similar, but they are not the same because the magnitude and direct-space constraints have a very different structure. The latter scheme was used by Oszlányi & Sütő (2011) and shown to be superior to the original CFA, especially when dealing with low-resolution data (Oszlányi & Sütő, 2011; van der Lee, 2013).

So far, we have not considered the consistency of the phase retrieval problem, and we have assumed that the solution exists. However, crystallographic phase retrieval is most often inconsistent. This is caused by the limited resolution of the diffraction data and by the presence of noise in the data. Moreover, if constraint (8) or some of its variants are used, the problem is inconsistent, because we are approximating the true ρ by a function which does not assume values between 0 and δ. The same is true for atomicity constraints, which assume zero density outside of a limited number of pixels. For inconsistent problems, the AAR, HIO, DM and related algorithms have a tendency to diverge from the solution (Marchesini, 2007; Fig. 2). This leads to a frequently observed problem of wandering of the iterates away from the solution. The error-reduction algorithm and the CFA do not diverge, and stay close to the optimal point (i.e. at the place with the shortest distance between the constraint sets), but they suffer from stagnation at local distance minima (Fig. 2c). An interesting algorithm that does not diverge for inconsistent problems, but inherits most of the ability of the AAR-related algorithms to escape from the local minima, is the relaxed alternating averaged reflection (RAAR) algorithm (Luke, 2005). This algorithm is a one-parameter relaxation of the AAR algorithm:

Figure 2 Convergence of selected dual-space algorithms on a two-dimensional example. (a) Two convex constraint sets with intersection. (b) Two non-convex constraint sets with several intersections. (c) Two non-convex constraint sets without intersection (unfeasible problem). All iterations start from the same point in the right-hand part of the plots. Symbols represent the actual iterates, the dotted lines connect consecutive iterates. ER = error-reduction algorithm (13), CFA = charge flipping algorithm (14), AAR = averaged alternating reflections (21), RAAR = relaxed averaged alternating reflections (24), DM = difference map [(20), with γM = β−1, γD = −β−1].

This algorithm has a fixed point, even if the corresponding problem is inconsistent, and as such appears to be a very attractive variant of applications in crystallography. Skubák (2018) applied this algorithm to the solution of heavy-atom substructures from anomalous difference data and found it better performing than the standard CFA. The best value of the parameter β was found empirically to be 0.82. It is possible that in contexts other than heavy-atom substructures the optimal β might be different.

Realistic phase retrieval problems operate in spaces of very large dimension. It is, however, very enlightening to observe the behaviour of the algorithms for a simple, two-pixel problem, which can be represented in a plane. Several of the presented algorithms (ER, the CFA, AAR, RAAR, DM) were used to solve a simple problem in two dimensions, where the two constraint sets are represented by two curves. Fig. 2 shows the results for a convex consistent problem, a non-convex consistent problem with multiple solutions and a non-convex inconsistent problem. Each symbol in the plots shows an iterate of the algorithms. Successive iterates are connected with a line, forming a path. For the convex consistent problem (Fig. 2a), all algorithms converge to the correct solution. For the non-convex consistent problem (Fig. 2b), ER and the CFA stagnate at local minima. However, the CFA is able to avoid some of the local minima, and approaches the solution more than ER. All other algorithms find one of the solutions. The non-convex inconsistent problem is the most challenging. ER and the CFA approach the solution, but stagnate at local minima. Again, the CFA avoids some of the minima that are trapping the ER algorithm. The AAR, HIO and DM algorithms all diverge from the solution. The RAAR algorithm converges close to the solution, and a point very close to the solution would be found by a single application of one of the projections. More details of the behaviour of various algorithms for general phase retrieval problems are discussed in an excellent overview by Marchesini (2007).

Most of the algorithms presented so far can be regarded as special cases of a general, six-parameter iteration scheme of the following form:

Table 1 gives the parameters of this general algorithm that correspond to the algorithms presented in this section. The only algorithm that cannot be represented as a special case of the above scheme is the general form of the HIO algorithm [equation (18)] and the HPR algorithm [equation (19) with the positivity constraint].

Table 1 Dual-space iterative algorithms expressed as instances of the general iteration scheme (25) β is the free parameter of the algorithms.   β1 γM,1 γD,1 β2 γM,2 γD,2 Gerchberg–Saxton/error reduction 1 0 0 0 – – Original CFA 1 0 1 0 – – Intrinsic phasing† 1 γM,1 0 0 – – HIO‡ β β−1 0 −β 0 −1 DM§ β β−1 0 −β 0 −β−1 AAR¶ 1/2 1 1 0 – – AAR, equation (23) 0 – – 1/2 1 1 RAAR β/2 1 1 1 − β 0 −1 †The recommended value of γM,1 is 2. ‡Strictly valid only for support constraint. §With optimal parameters derived from the assumption of locally orthogonal constraint sets. ¶And also HPR with β = 1, equation (21).

This overview of algorithms should provide an idea of the complexity of the field, and point out the similarities, but also the differences, between the algorithms. The flowcharts of the most important algorithms are summarized in Fig. 1. The overview should also make it clear that the potential and performance of various algorithms in a crystallographic context is far from being entirely understood and explored.

In the previous section, various phase retrieval algorithms were reviewed. Although several algorithms have been applied to crystallographic problems, the charge-flipping-based algorithms remain the most studied. This section therefore provides a detailed overview of variants and flavours of the CFA. Many of the findings on the CFA might apply also to other dual-space algorithms.

The first and obvious improvement of the efficiency of the algorithm is not related to the algorithm itself, but to the data employed. The original paper on the CFA (Oszlányi & Sütő, 2004) used as an input the magnitudes of the standard, non-normalized Fourier coefficients. Using normalized Fourier coefficients (the E values) yields sharper maps and thus a much better performance of the algorithm. This general knowledge was quantitatively probed by Oszlányi & Sütő (2008), who showed, on a selected example, a reduction of the number of iteration steps by about two orders of magnitude upon introducing normalized Fourier coefficients. Other dual-space algorithms (Lei, 2007; Feng, 2012) directly employ the normalized coefficients. In the IPA, the magnitudes entering the iteration are a mix between normalized and unnormalized structure-factor amplitudes in the form G = E^qF^1−q, with q being a user-defined parameter with a default value of 0.5.

As described in the previous section, the original version of the CFA employed the iteration scheme (14) with the magnitude projection (6) and direct-space projection (8). Owing to its crucial role in the CFA, we give here the corresponding reflector of (8) explicitly:

This is the so-called charge-flipping operation, which gave the algorithm its name. The exact form of the magnitude constraint (6) is important. For example, replacing the constraint by the closely related form (2) has a devastating effect on the efficiency of the algorithm. The zero Fourier coefficient F(0) also deserves special attention. This coefficient is not known experimentally. In the original formulation of the CFA [ equation (6)] it was left unconstrained. The variant with constraining F(0) to zero [i.e. using projector (4)] was also sometimes used (Palatinus, 2004; Coelho, 2007a; Zhou & Harris, 2008). However, leaving the F(0) coefficient unconstrained seems to be the most efficient approach. If left unconstrained, the abrupt drop of the value of F(0) can also be used as an indicator of the convergence of the iteration (Oszlányi & Sütő, 2004).

The parameter δ is the single free parameter of the basic algorithm. Its value on an absolute scale differs from one problem to another. However, it was shown (Oszlányi & Sütő, 2008) that δ can be conveniently expressed in terms of the standard deviation of the reconstructed density: where σ(ρ) is the standard deviation of the distribution of the density values. The optimal value of k_ed was shown to be most often between 0.9 and 1.3.

Soon after the publication of the algorithm, first applications and improvements of the algorithm appeared. Wu et al. (2004a), along with the first application to experimental data, proposed to replace the constant δ in the charge-flipping operation by a dynamical δ determined newly every cycle so that a constant number of pixels is flipped. While this modification does not seem to have an important effect for most structure-solution problems, it appears to have a stabilizing effect for problems where the solution is less stable.

Naturally, most efforts concentrated on the improvement of the phasing power of the algorithm. These attempts focused either on the modification of the constraints or of the iteration scheme. The first class involves the so-called π-half variant (Oszlányi & Sütő, 2005), where the magnitude constraint is modified to [cf. equation (6)] The threshold between weak and strong reflections is selected so that a certain fraction of the reflections – typically 20–30% – are treated as weak. This modification dramatically improves the performance of the algorithm.

Another improvement of the operation on the Fourier magnitudes was the replacement of the simple magnitude projection by this operation (Oszlányi & Sütő, 2008): Here is the phase of the coefficient . It can be regarded as a standard F + ΔF Fourier synthesis used commonly in macromolecular crystallography. This operator resembles a reflector, but a true reflector would have to change the sign of all unobserved Fourier coefficients (). The authors recommend that a limit is imposed on the change of |F_j| so that , where W is a new free parameter of the algorithm. In some sense this variant changes the iteration scheme from the standard CFA scheme (14) to (17). The damping introduced by the parameter W and by not applying the operation to unobserved reflections is necessary to ensure that this algorithm does not suffer from the instability typical for the iteration scheme (17). This modification provided similar or somewhat better results than the π-half variant on a test organic structure (Oszlányi & Sütő, 2008). A very similar reciprocal-space operation is used in the IPA.

The modifications described so far aimed at the improvement of the efficiency of the algorithm. A variant called band flipping (Oszlányi & Sütő, 2007) instead aims at lifting the requirement of the positivity of the direct-space constraint. It employs the `dynamical support' constraint (9) instead of the standard constraint (8). The dynamical support constraint does not enforce positivity of ρ. The action of the corresponding reflector (the band-flipping operator) is to change the sign of all pixels with − δ < ρ_i < δ. This variant has a weaker phasing power than the standard variant, but is applicable to neutron scattering densities with negative scatterers (Oszlányi & Sütő, 2007), or to the reconstruction of difference electron densities, and hence to the solution of superstructures (Palatinus et al., 2011). An asymmetric interval for the flipping proved to be useful in the solution of structures from neutron diffraction (Eq. 10 of Whitfield & Coelho, 2016).

Oszlányi & Sütő (2011) combined the charge-flipping operation with the AAR iteration scheme (23). It was shown that the AAR iteration scheme clearly outperforms the standard charge-flipping scheme. The improvement is especially striking for low-resolution data. However, the tests were performed on synthetic, noise-free data. The advantage of the AAR algorithm over the standard CFA is less striking for experimental, noisy data, owing to the tendency of the AAR scheme to diverge for inconsistent problems (Marchesini, 2007). The advantage for low-resolution data, however, remains visible even in experimental data (van der Lee, 2013).

The challenge of the dual-space algorithms is finding the balance between the perturbation strength and stability of the solution at convergence. One possibility to influence this balance is to design a moderately perturbing algorithm, but add additional perturbation in the projection step, i.e. replacing the perfect projection operation with an operation that adds additional perturbation or suppresses certain unwanted features in the density. A prominent example of this strategy are so-called omit maps. Omit maps are frequently used in macromolecular crystallography to remove model bias. Oszlányi & Süto (2016) analysed the effect of resetting parts of ρ to 0 at every nth iteration (n between 1 and 10). The most efficient of the tested variants turned out to be dividing the unit cell by a randomly oriented and randomly positioned plane into two equally large halves, and setting to zero one of the halves. This apparently drastic perturbation improves the efficiency of the CFA and turns also other, less well performing variants like LDE into viable phasing algorithms. The random omit procedure is also an integral part of the IPA. In the stage of construction of the mask consisting of Gaussians placed at the peaks of ρ, a small number of these Gaussians are removed at random, creating an additional perturbation in the iteration.

Another example is the modification introduced by Coelho (2007a), who proposed replacing the charge-flipping operation (26) by the recipeThis operation has the effect of damping the highest maxima in the density and thus avoiding false solutions consisting of a single very large peak (known as a `uranium atom solution').

The simplicity of the CFA makes it easy to combine it with other iterations schemes or completely different solution strategies. The CFA was combined with histogram matching in powder diffraction (Section 7.3). Coelho (2007a) combined the CFA with the tangent formula (i.e. classical direct methods) to obtain an algorithm that merges the two worlds. The algorithm proposed by Coelho contains several modifications with respect to the original algorithm, but the principal one is the introduction of the tangent formula in the Fourier-space modification step. Instead of keeping the phases of the Fourier coefficients intact [cf. equation (2)], they are shifted towards phases obtained by the application of the tangent formula. With this approach, a substantial improvement of performance was obtained. This algorithm was implemented in the crystallographic suite TOPAS (Coelho, 2007b).

The potential use of the symmetry information in the dual-space algorithms, especially in the CFA, has been subject to a recurring debate over the years. All major publications on dual-space structure solution methods (Matsugaki & Shiono, 2001; Elser, 2003b; Oszlányi & Sütő, 2004; Sheldrick, 2015) note that symmetry has not been used in their tests, and Elser (2003b) and Oszlányi & Sütő (2004) express the hope that the proper use of symmetry will improve the power of the algorithms. On the one hand, one could argue that the symmetry is already implicitly contained in the data, and adding it explicitly does not provide additional information. On the other hand, the hope for improvement was driven by the observation that by explicitly introducing phase relationships due to the space-group symmetry in the iteration, the number of phases that needs to be determined in the iteration is reduced. This would mean a reduction of the dimension of the search space by a factor approaching the order of the space group, i.e. up to 48 for some cubic space groups.

Interestingly, apart from specific exceptions discussed later, an improvement by using symmetry was never accomplished, and application of the algorithms without any symmetry constraints remains the most efficient approach. Although no mathematically rigorous analysis of the problem has yet been published, an intuitive explanation of this fact is available. If symmetry is imposed on the density, then all features must develop from a random density exactly at the positions determined by symmetry. Without symmetry restrictions (i.e. `in the space group P1'), the structure can develop anywhere in the unit cell, giving the algorithm much more freedom to randomly develop a `seed' of the correct structure, and to converge to a complete solution from that seed. Moreover, the parameter δ of the flipping operation must be set to find a balance between the perturbing effect of the operation and the stability at the solution. If the symmetry is fixed – for simplicity, let us consider just the presence of an inversion centre at the origin – all Fourier-coefficient phases are fixed to either 0 or π. Switching the phases of important reflections from 0 to π requires an extremely strong perturbation of the density in real space. If δ is set so high as to permit such changes, it will be too high to stabilize the iteration at the solution. If δ is smaller, the phases of the most important reflections will be fixed and the iteration will stagnate. Interestingly, a similar observation was also made in the framework of direct methods (Sheldrick & Gould, 1995; Burla et al., 2000) and a procedure called RELAX, which relaxes the symmetry constraints on the structure (Burla et al., 2002), has become a standard part of the structure solution process in the program SIR2011 (Burla et al., 2012a) and later.

An apparent contradiction to the arguments just stated is the method due to Eggeman et al. (2009), which used symmetry-enhanced charge flipping to solve a two-dimensional structure from zonal electron diffraction data. The approach is the following: first run the CFA in P1 for a couple of cycles, and then symmetrize the density regularly every few cycles. This approach stabilized and improved the solution. But the contradiction is only apparent. In the particular case of two-dimensional electron diffraction data, the problem is not to reach convergence. The structure is very simple to solve, but it is difficult to stabilize the solution owing to the limited accuracy and limited amount of data. In such cases, applying the symmetry may indeed help to find the best solution and stabilize it by adding more constraints to the problem. A similar effect can be expected for structures solved from powder diffraction data. Indeed, a possibility to partially or completely impose the symmetry on the current iterate has been implemented in the charge-flipping routine of the program TOPAS Academic (Coelho, 2007b), and is reported to have a stabilizing effect on the structure solution from low-resolution or powder diffraction data (Coelho, 2012). Partially imposing the symmetry on the iterates has also been used by Coelho (2021) in a work focused on the ab initio solution of protein structures by the CFA. In this work, enforcing the symmetry in the iterates during iteration is interpreted as narrowing the searched phase space, but increasing the radius of convergence of the iteration. The narrowed phase-space search is compensated by several other modifications of the basic algorithm that increase the perturbation strength and thus broaden the searched phase space.

The fact that the algorithm performs best without any symmetry restrictions can be turned into an advantage. If the solution is found without symmetry restrictions, then it can be analysed for the presence of symmetry elements, and the most probable space group of the structure can be deduced after the solution (Palatinus & van der Lee, 2008; Sheldrick, 2015). This approach is fundamentally different from the standard space-group determination using the symmetry and systematic absences in diffraction data, because it uses the phased Fourier coefficients and not just their magnitudes. It thus does not suffer from the ambiguities present in the standard approach, and allows, for example, an unambiguous discrimination between space groups differing only by the presence/absence of an inversion centre. This approach, on the other hand, is less sensitive to small deviations from higher symmetry, which can often be reliably revealed in Fourier space. Ideally, the best estimate of the space-group symmetry should be obtained by combining both methods.

Incomplete diffraction data are a severe problem for the structure solution process regardless of the solution method. Several methods have been devised to overcome the problem. The missing coefficients can be extrapolated by imposing positivity on the Patterson map (Karle & Hauptman, 1964; Langs, 1998) or estimated using probabilistic formulae relating the unknown magnitudes either to the experimental observations (David, 1987; Cascarano et al., 1991; Xu & Hauptman, 2000) or to the Fourier magnitudes of a model density (Caliandro et al., 2005a,b, 2009).

The dual-space algorithms are a Fourier-based technique, and thus the problem of an incomplete data set is probably even more critical here than in other methods. In the original formulation of the CFA, the magnitude constraint had the form (6), i.e. all unmeasured Fourier magnitudes except for F(0) were reset to zero. This severely hampers the algorithm's performance, even if only a few important low-order Fourier magnitudes are missing. A partial remedy to the problem is to replace the projection (6) by projection (5), possibly with infinite c. This modification solves the problem of missing low-order data to a large extent. For cases of an extreme amount of missing data, Palatinus et al. (2007) proposed a method based on the optimization of the Patterson function by the maximum entropy method (MEM). The optimization of the Patterson map by MEM leads to an estimation of the missing Fourier magnitudes, which can then be used as experimental data in the charge-flipping iteration. Using this technique, test structures could be solved with more than 50% of the reflections missing inside the resolution sphere. The method is, however, relatively tedious, time consuming and not very efficient in extrapolating the data to higher resolution.

Compared with the standard algorithm with the constraint (5), significant improvement of performance with missing data was reported with the AAR variant (Oszlányi & Sütő, 2011). The published tests show that in some cases the AAR algorithm can solve purely organic, non-centrosymmetric structures from low-resolution data (d_min = 1.5 Å), while the standard algorithm fails already slightly above d_min = 1.2 Å. For centrosymmetric structures, solution from data with d_min = 1.6 Å was easily possible with AAR, while it was very difficult or impossible with the standard CFA. A beneficial effect of the AAR scheme was confirmed on a large set of experimental data sets by van der Lee (2013).

In the IPA, the problem of an incomplete data set is less severe. The replacement of the peaks in the iterate by the idealized Gaussian peaks provides the additional information about the electron density that automatically leads to an effective estimation of the missing Fourier amplitudes. This approach is very efficient, but it may fail in the rare cases when the correct electron density is not composed dominantly from isolated round maxima, for example in cases of severe disorder in the structure.

No modern computing method can hope for widespread usage without publicly available software implementing the method. It is likely that one of the key reasons for the success of the CFA is that a rich collection of such software is available. Quickly after publication, the CFA has become available for users as a module in the crystallographic software package PLATON (Spek, 2003), and in the computer program BayMEM (van Smaalen et al., 2003). The CFA is also implemented in the program TOPAS (Coelho, 2007b). This program contains an implementation with several special features (Coelho, 2007a, 2021). It is focused on structure solution from powder diffraction data and includes the powder CFA with histogram matching (Section 7.3), but can also be used with single-crystal data. Another program for structure solution using dual-space algorithms is Superflip (Palatinus & Chapuis, 2007). This program allows the application of the CFA in arbitrary dimensions, allowing the solution of ordinary periodic structures as well as modulated structures and quasicrystals (see Section 7.2). The program also contains many of the modifications of the charge flipping algorithm described in Section 3 and the symmetry-determination algorithm due to Palatinus & van der Lee (2008). It also contains the general, six-parameter iteration scheme [equation (25)], allowing the application of a wide variety of iterative algorithms (Table 1). The IPA (intrinsic phasing algorithm) is implemented in the popular structure determination program SHELXT (Sheldrick, 2015). This program is optimized for the solution of small-molecule crystal structures, contains its own symmetry determination algorithm and provides an advanced algorithm for the chemical interpretation of the obtained electron density.

Dual-space algorithms are still a relatively young class of phasing algorithm, yet they have already established themselves in the portfolio of tools for the solution of the phase problem. Many variations of the algorithms exist and numerous applications testify to their usefulness. The main problem in the further development of phase retrieval algorithms is that no solid mathematical theory is available that would allow the determination of the perfect algorithm. A lot of insight can be gained from the analogies with the convex optimization theory, but the results are not directly applicable. Moreover, an algorithm is not just an iteration scheme, but a combination of the iteration scheme and the exact definition of the constraints and projections. A successful algorithm is a fine balance between all components, and a small, seemingly unimportant change can result in a change of efficiency by several orders of magnitude.