Charge Flipping in GSAS-II - sucrose


In this exercise you will use GSAS-II to solve the structure of sucrose from powder diffraction data via charge flipping. This structure was originally solved from single crystal data in the space group P21. It is a difficult exercise to complete as the charge flipping does not always yield a solution and to even have a chance the Pawley refinement must be done with particular care. We have given a route that sometimes works.

If you have not done so already, start GSAS-II and open the GSAS-II project file (I called it sucrose.gpx) you saved from the third part of the peak fitting-indexing exercise. If you didn’t do that exercise, do it now.

Step 1: Setup for Pawley refinement

There are a number of steps that must be done in preparation for a Pawley refinement after having completed the unit cell indexing and new phase creation. These cover two things: one is to prepare the powder pattern (reset limits, etc.) and the other is to prepare the phase for the refinement.

1.      Select Limits in the GSAS-II data tree for the PWDR data set and expand the plot so that the region near d=1.0Å is readily seen. This is ~24deg 2Q. There is an area relatively clear of peaks that makes a good place to end the calculations with a reflection set of dmin ~1Å resolution needed for charge flipping. Enter 23.55 for Tmax as the upper limit and make a note of the exact d-spacing (1.0125Å) using the mouse cursor on the plot.

2.      Select Instrument Parameters and uncheck all Refine flags and do Operations/Reset profile to recover the default values. The peak fitting done earlier is over a much more limited range of 2Q giving values that might not extrapolate well over the wider range to be used in the Pawley refinement. We will refine them again later.

3.      Select Sample Parameters and uncheck refinement of Histogram scale factor. While you are here set Goniometer radius to 1000.

4.      Go to Phases/sucrose to find the Pawley controls. Check the Do Pawley refinement box and enter the d-spacing (1.0125) into the Pawley dmin box.

When done the General tab should look like

Next find the Data tab; it will be empty except for a single line of text. Do Edit/Add powder histograms; a dialog box will appear

Select either choice; the desired data set will be added to this phase for analysis. The data tab now shows the new data set.

I’ve stretched the window slightly to see the entire contents. This is the location for all the phase dependent parameters for this histogram. Notice that it includes phase fraction, size & mustrain as well as preferred orientation, elastic strain and extinction corrections. The Babinet correction is intended for protein work where a significant region of the structure is disordered solvent. There are buttons for plotting size and mustrain surfaces

as well as preferred orientation correction curves.

Next find the Pawley reflections tab; it will be empty except for column headings. Do Operations/Pawley create; this makes the reflection set over the range covered by setting Pawley dmin. The table should list reflections 0-780. Select and check the refine column using the same technique you used for the peak list.

This completes the Pawley refinement setup; be careful that you didn’t skip a step as it might not work correctly if you did.

Step 2: Pawley refinement

To start the refinement, do Calculate/Refine on the main GSAS-II data tree window. Before it begins a backup of the project file is made; the name will be sucrose.bck0.gpx. It can be used to recover from a bad refinement. A progress dialog box will appear showing the residual as the refinement proceeds. Notice that there was a singular matrix warning on the console window and the last reflection (781) was removed from the refinement because it was off the end of the used part of the pattern. This is the GSAS-II singular matrix recovery mechanism; you can go to the Pawley reflections tab and uncheck the last reflection. Then repeat Calculate/Refine to make sure it has essentially minimized. My refinement completed with Rwp ~26.3%; a dialog box appears asking if you wish to load the result. Press OK. To see the plot select the PWDR line in the GSAS-II data tree (I’ve adjusted the width and height).

It is evident from examining the plot that some improvement can be had by refining the lattice and profile parameters and perhaps the sample X position. Go to Phases/sucrose in the GSAS-II data tree and check the Refine unit cell box. Now select Instrument parameters for the PWDR data set in the GSAS-II data tree and check U, V, W, X & Y. The go to Sample Parameters and check Sample X displ.; make the Goniometer radius 1000 if you didn’t do this earlier. Repeat Calculate/Refine until converged; I obtained an improved residual of ~9.05% with a much better fit at the high angle part of the pattern but still not fit well for the lower part.

What you see is that the background has a substantial diffuse scattering feature at ~5deg and there are some peak shape issues (which may be a consequence of the poor background fit). We’d like the fit to be better before trying charge flipping. First select Background, set Peaks in background: No. peaks to 1 from the pull down box; the window changes to


Enter 5.0 for pos, 100000 for int and 1000 for sig (as appropriate for broad diffuse peaks) and then check all three refine boxes. Also set No. coeff. pulldown to 5 for the chebyschev part of the background. The window should look like

Go to the Controls item in the GSAS-II data tree and set the Max cycles to 5. Now do Calculate/Refine until the refinement converges. After a set of refinement cycles go to the Pawley reflections tab and check the intensities. If there are many negatives amongst the highest angle reflections do Pawley update, this replaces all the values with those in the extracted from the profile as listed in the Reflection list for the PWDR data set each multiplied by 0.5. Additionally their refinement flags are turned off. Do Calculate/Refine; this will adjust all the reflections in the vicinity of the reset ones. You may repeat this process to suppress any new negatives that may appear. Then set all the Pawley refinement flags back on (excluding any that lead to a singular matrix) and repeat the refinement. After refinement, mine converged to ~5.3% with very few negatives and and the plot looked like

Let’s see if this is good enough for charge flipping.

Step 3: Charge flipping

A few steps are needed to set up for a charge flipping run. Go to Phases/sucrose and notice the Charge flip controls in the middle of the General tab.

1.      Select Reflection set from to be PWDR 11bmb_8716.fxye Bank 1 (the only choice!). This makes the reflection set to be used those in Reflection Lists under the PWDR data set; these should match the Pawley reflection set.

2.      Set the Peak cutoff (in Fourier map controls) to 20%. This controls the map search routine.

For powder data it is not usually advisable to select a Normalizing element; if selected a periodic table appears – select an element or ‘None’. The Resolution defines the map spacing and thus the extent of the reflection set; the observed reflection set is zero filled to this limit for charge flipping. The k-factor and k-Max determine the lower threshold and upper limit for charge flipping, respectively. The latter is used to avoid “Uranium solutions” in the result. Set this to 10-12 for equal atom problems and larger for heavy atom ones. For sucrose set it to 12. When ready select Compute/Charge flipping; a progress bar dialog will appear that tracks the residual found in each charge flipping cycle. It will quickly drop to some level (42-44%) and then for a successful charge flipping run perhaps drop again into the 20-25% range after a considerable length of time. You should save the project if you get an apparently good solution. Many times it will not drop below ~40%; just try Compute/Charge flipping again. If several tries fail to yield a solution, go back to the last step of the Pawley refinement and just repeat it (that is do Pawley update and Calculate/Refine). This will generate a different set of intensities for the highly overlapped part of the powder pattern; perhaps they will give a successful charge flip solution. Unless otherwise interrupted, GSAS-II will try 10000 times for a charge flip solution (for sucrose this takes 5-10min on my laptop). Press Cancel to stop the charge flipping; a peak list and a drawing of the solution will appear. One of my runs gave


The result looks very promising! Charge flipping operates without regard to symmetry and GSAS-II attempts to locate the map with respect to the symmetry elements from the computed reflection phases. (NB: charge flipping also will give a solution without regard to handedness; this solution may be either the correct one for sucrose or its inverse.) Errors in these phases can disrupt this process so that the map and peak positions are offset from their true positions. To discover this and fix it, rotate the drawing around to check that atom positions are related via the 21 screw axis (through the multicolored cross in the center of the cell and parallel to the y-axis – the green cell edge). Avoid motions with the right mouse button pressed, this moves the cross. To recover, press the ‘C’ key to reset the cross to the cell center. I’ve placed the symbol for a 21 axis on my drawing above to show where I think it should be. If you find an offset to the map, orient it so the offset is horizontal or vertical and then press the L, R, U or D keys to shift the map and peaks as needed. The shift is in resolution units (0.5Å) so the fix can only be approximate. Repeat as needed. Notice that the peak positions in the table are updated as the offset is changed. In this case, I needed just a few steps to correct the map offset. NB: you can use these keys to move the result anywhere you wish, e.g. to select a different origin for the structure but pay attention to the location of symmetry elements. After adjustment, my solution looked like

Step 4: Obtain solution from charge flip result

The next step is to extract the unique peak positions from the map and make atoms out of them. The map peak list is initially sorted by peak magnitude; you can change the sorting to be along any of the column headings in the table by a single click on the desired heading. The dzero column gives the distance of each peak from the origin. The menu items under Map peaks give you several tools to aid in the peak selection process; we will just use a couple in this example. I have chosen to proceed by sorting the atoms by z, selecting all and then finding the unique ones.

1.      Do a single click on the z column heading. This will sort the peaks in increasing z-axis coordinate.

2.      Do a single click on the upper left (empty) box of the table. All entries will be colored grey indicating their selection.

3.      Do Map peaks/Unique peaks. That should select (grey) 23 out of 46 peak positions in the table and all should be near the top of the list. Given the formula for sucrose is C12H22O11 that would be all the C & O atoms! They will be colored green in the plot. If more are selected then check the map offset – it may need a slight correction.

The drawing shows the selected atoms in green.

Notice that in my case one atom is in a neighboring molecule, one can deselect it and then pick the desired equivalent. Hold the Shift key down and use the right mouse button to pick/unpick atoms in the drawing. NB: the left mouse button will clear all the selections. The other two “lone” atoms belong to molecules in neighboring unit cells; we’ll deal with them later.

These are now ready to turn into atoms. Do Map peaks/Move peaks. The drawing will show white balls for each of the selected positions and the Atoms table will now have 23 new H atoms.

The atom names indicate the peak magnitude to aid identification (M100 for the largest in the Map peaks table which may not be one of those moved to the Atoms table). You can hand sort them by selecting a row with the Alt key (or Shift/Ctrl keys) down and then selecting a position in the table with the Alt key still down. On Linux machines the Alt key may have been hijacked by the operating system; use Shift/Ctrl instead to move atoms. After sorting by magnitude, my list is now

Notice that in this case the atom list is neatly sorted out into two groups with a clear gap between M79 and M64. The higher ones are likely to be the 11 O-atoms and the rest are the 12 C-atoms. There may be a bit of ambiguity in the middle, but just assign the top 11 as O and the rest as C and see what happens. You can do this by selecting a block of atoms, go to Edit/Modify atom parameters and select Type from the dialog box; a periodic table will appear. Select the desired element and they will be changed in the list and drawing.

The original charge flipping solution is still visible and the atoms are now identified, I’ve changed the drawing to show balls and sticks and labeled the atoms with their names, this is done in the Draw atoms tab. Double click the appropriate column heading and make the selection from the dialog box. 

Now one can safely remove the map peaks (do Map peaks/Clear peaks in the Map peaks tab) and clear the charge flip map (do Clear map in the General tab). The atom positions I obtained have two atoms (O3 & O4) not connected to the sucrose molecule because they belong to neighboring molecules. We need to discover the equivalent ones that belong to this molecule. Follow the steps:

1.      Go to the Draw Atoms tab and select them (hold the Ctrl key down and select the appropriate rows).

2.      Do Edit/Fill unit cell; new atoms will appear at end of list and in the drawing.

3.      Select all the atoms in the Draw Atoms table by a double click on upper left corner box of the table.

4.      Do Edit/Fill CN-sphere; new atoms will appear that are bonded to the molecule.

5.      Select these new atoms that are in the molecule; they will probably be outside the unit cell box. Do this on the drawing with Shift left for the first one & Shift right for the rest, they will show in green.

6.      Look at the Draw Atoms table for the selections and make note of the Sym Op for each one. In my case they were “2+1,0,0” for O(4) and “2+2,-1,1” for O(3). The first number is the symmetry operator number and the 3 numbers after the ‘+’ are unit cell translations.

7.      Go to the Atoms tab.

8.      Select an atom in turn that needs to be moved, e.g. in my case O(3) or O(4).

9.      Do Edit/Transform atoms, select the operator and unit cell translations, then press OK. The atom coordinates will be transformed.

10.  Repeat 8 & 9 for any other atoms to be moved.

This will complete the molecule with the extra atoms. Do Edit/Reload draw atoms to see the completed molecule. I’ve made it a balls & sticks with names as labels.

Save the project file; the next step is to complete the structure by finding H-atoms and then do the final Rietveld refinements.

Step 5: Complete the structure

To complete the structure you will need to perform a preliminary Rietveld refinement. This will make the atom positions more accurate and prepare a set of calculated structure factors needed to generate a difference Fourier map good enough to maybe reveal the H atoms. There are a number of steps needed to set all the controls from the previous Pawley refinement to get a successful start for the Rietveld refinement.

1.      Go to the Instrument Parameters entry in the GSAS-II data tree and uncheck all Refine boxes. You don’t want to start with these when the scale factor is unknown.

2.      Go to Sample Parameters and check the Histogram scale factor box. This is the scale factor for the data set. Uncheck the Sample X displ.

3.      Go to Phases/sucrose General tab and uncheck Refine unit cell and Do Pawley refinement boxes. The latter is most important; if you forget you will just get another Pawley refinement!

Now do Calculate/Refine from the main GSAS-II data tree menu; only scale and background are refined. I got a Rwp ~20%. The plot shows clear intensity differences; the atom positions obviously need refinement.

I’ve shifted & zoomed in to make the differences more obvious.

Go to Phases/sucrose and select the Atoms tab. Note the warning at the bottom of the window; the origin in space group P21 has an arbitrary location along the y-axis. Then double click the refine column heading, a dialog box will appear, select the X-coordinates box and press OK. The atoms table will show the refinement controls

To deal with the warning we need to set a “Hold” on one atom y-coordinate. Select Constraints in the main GSAS-II data tree.

Do an Edit/Add hold and select ‘0::dAy:0 for O(1)’ from the dialog box. This will fix the shift along y for the O(1) atom.

Press OK; the Constraints window will show the hold

Do Calculate/Refine again; my Rwp is now ~10.3%. You can now safely refine lattice parameters (in the General tab for Phases/sucrose), the Instrument parameters U, V, W, X & Y and the Sample X displ. You can also include refinement of atom thermal motion (Uiso) parameters. My refinement improved slightly to Rwp~10.0%. Let’s see if a difference Fourier map will show anything.
Go to the Phases/sucrose General tab; the Fourier map controls are midway down the window.

1.      Select delt-F for the Map type.

2.      Select PWDR 11bmb_8716.fxye Bank 1 for the Reflection set from item (the only choice). The Fourier calculation will use the structure factors shown in the Reflection list item for this PWDR data set.

3.      Do Compute/Fourier map from the General tab window. The map will be computed via fast Fourier techniques and is very fast. A few lines summarizing the map is printed on the console, the structure plot may show a single green dot somewhere (it could be hidden beneath an atom!), this is the highest point in the map.

4.      Go to the Draw Options tab and adjust the Contour level slider until something appears on the map. You can also reduce the Map radius to help clarify the view.

What is drawn here are green dots whose size are proportional to the electron density. There are peaks that might be H-atom positions as well as features that suggest C & O-atoms out of position. Recall that the refinement was just done without any H-atoms; it will make the best fit compensating for the missing H-atoms by displacing the C & O atoms in the direction of the missing H-atoms.

1.      Go to the Phases/sucrose Draw Atoms tab.

2.      Select a C or O atom with possible H-atom density nearby. It will appear green. Do Edit/View point; this will position the cross on the selected atom.

3.      Using the mouse, position the multicolored cross directly in the center of the nearby H atom peak. Rotate the drawing around with the left mouse button and shift the cross with the right mouse button. You can zoom in using the roll wheel on the mouse (or else use the Camera distance slider in the Draw Options tab if your mouse lacks a wheel).

4.      Go to the Atoms tab and do Edit/Insert view point, a white ball will appear at the selected position and an H-atom (name = UNK) will be appended to the end of the atom list. Repeat steps 1-4 for as many recognizable H-atom positions you can find; the goal is to find all 22 of them (you probably will find less than that, but you could guess many of the rest). I easily found all but 5 of them.

5.      Change their Type to H, they will be renamed and make then all ball & stick. Do a Calculate/Refine; there should be a drop in Rwp. Then repeat the Compute/Fourier map in the General tab, adjust the Contour level and Map radius as desired. Look for the remaining H-atoms; you may have to guess positions for some of them using known stereochemistry rules. My model looked like

Do a Calculate/Refine to finish this stage in the structure analysis. My Rwp was ~8.02% with all the H-atoms included (but not refined). The fit looks pretty good

Step 6: Final refinement

Looking at the fit, there may be room for improvement and we’d like to refine the H-atom positions. First we need to change the profile model to include crystallite size and mustrain broadening. Go to Instrument Parameters, set X and Y to zero and uncheck the Refine flags for them.

Next go to Phases/sucrose and look at the Data tab; make sure the Show box is checked. Check the boxes for both Cryst. size and microstrain.

Do a Calculate/Refine; the Rwp will start out high but should converge to where it was before this change. I got 8.08% and the Data tab looked like

Next select Generalized for the Mustrain model; this is the Stephens phenomenological model for microstrain which has 9 terms for monoclinic structures. Check all the boxes for it.

Do a Calculate/Refine; the Rwp will start out high but should converge lower than where it was before this change. I got 7.65% and the Data tab looked like. If you check the Mustrain button in the Select plot type box, you will see the interesting mustrain surface plot for ground sucrose.

Before refining H-atom positions it would be wise to create restraints on their bond distances and angles. Locate Restraints in the main GSAS-II data tree; there will be an almost blank window telling you there are no bond restraints for sucrose. Do Edit/Add restraints, a dialog box will appear asking for the origin atom for the bond.

Select all H and press OK. A new dialog box will appear asking for the target atoms for the bonds.

Select all C and all O and a third dialog box appears asking for the bond length


Enter 1.09 and press OK. I got a short list of 8 H-C and H-O restraints, but I expect 22 of them. Change the Search range to 1.5 and repeat Edit/Add restraints. There are now 22 restraints on H-C and H-O bond lengths.

Next select the Angle restraints tab and again there aren’t any. Do Edit/Add restraints, a dialog box will appear asking for the A atom in the angle A-B-C. Select all H. The next dialog asks for the B atom in the A-B-C angle, select all C and all O. The third dialog asks for the angle and has a default of 109.54. This is the tetrahedral angle so just press OK. 36 restraints will be listed; they are listed in order of the B atom. O-atoms should have one angle apiece; C-atoms will have more.

Now go back the Atoms tab in Phases/sucrose, select H from the dialog box by double clicking on the Type column. Then do Edit/Set atom refinement flags. Select X-coordinates; this will set those for all H-atoms. Now do Calculate/Refine, the console will list progress on the “Penalty function” (i.e. the restraints) and the Rwp drops in my case to ~7.05%. One can increase the Restraint weight factors to tighten them up at some expense in Rwp. In the end I chose 2.0 for both Bond and Angle restraints to give a final Rwp ~7.1%. NB: the H-atom positions may not be correct even though the refinement apparently succeeded. You should check that the OH ones all give reasonable H-bonds to neighboring O-atoms. As this is powder data such changes will hardly affect the fit but will be more satisfying. At present such checking is outside the scope of GSAS-II.