-
Notifications
You must be signed in to change notification settings - Fork 6
Ctrmodel
To get started type:
>>>import sxrd
and then
>>>sxrd.start()
in the PDS shell. This should make the python interface structure refinement window pop up.
To start working on your CTR data you need to have several input files prepared. You can read these files into the program using the 'Read Files' Menu. You need:
a datafile: *.dat
a bulk file: *.bul
a surface file: *.sur
a parameter file: *.par
a rigid body file: *.rbf
a bond valence file: *.bvc
Datafile, bulk file, surface file, and parameter file are necessary input, rigid body file and bond valence file are optional. As it's tedious to read them all in separately every time you start the GUI you can save a session in a .ssn file (this is still very basic, it just saves the paths to the input files) and you can read all files at once by loading the session.
The datafile contains a list with:
H K L F Ferr Lb Db
for all your datapoints. Lb denotes the L value of the lowest bragg peak. Note that you need to have an separate (positive) Lb value for the negative side of a rod. Db is the distance between two bragg peaks in L.
After you read in the datafile a list of all your rods should appear on the "Main controls" panel of the window. With each of the rods you have a text box that allows you to specify a rod weight. Rod weights can be used to emphasize or ignore (weight = 0) a certain CTR in the fitting procedure.
The bulk file contains a list of the atoms in the bulk unit cell:
Atom_label X Y Z Uiso
There has to be one line in the bulk file specifying the cell parameters
cell a b c alpha beta gamma delta1 delta2 NLayers
This line is indicated by the leading 'cell'. There have to be 9 numbers after ‘cell’. delta1 delta2 are the corrections for incommensurate surfaces, NLayers is the number of layers in the unit cell. In the text box next to bulk file on the window the cell parameters are displayed.
The surface file contains a list of all atoms in your surface unit cell:
Atom_label X Y Z Ux Uy Uz Uxy Uxz Uyz occupancy
Besides the list of atoms there is a float and a string defining how a parameter is used and defining the name of the parameter used to modify a certain value for each of the values given in the list of atoms (despite Atom_label). Note that parameters for XYZ are delta values while parameters for Ux ... replace the number from the surface list.
That means before every fit X will be replaced by:
X = X + dx_float * dX_parameter
while Ux for example will be replaced by:
Ux = Ux_float * Ux_parameter
The floats can be used to keep the surface symmetry during the fit. They may e.g. enforce that U tensors for atoms which are symmetry equivalent according to a symmetry operation 'R' fulfill the equation:
U' = R * U * RT
All the parameters defined in the surface- (and the rigid body-) file plus zwater, sig_water, Scale, specScale, and beta (roughness value after Robinson) must be listed in the parameter file.
The parameter file contains a list with
Label value min max refine_flag
for all parameters. In the newer version a standard deviation may be given as well after value.
zwater, sig_water, sig_water_bar, d_water, Scale, specScale, and beta are global parameters which cannot be changed. zwater, sig_water, sig_water_bar, and d_water are only needed if the “use bulk water” box in the upper right corner of main control panel is checked. They specify a layered water profile used to resemble bulk water.
zwater - is the distance between origin and the first bulk water layer (fractional coordinate)
sig_water - is the Debye Waller factor of the first bulk water layer (an Uiso like value).
sig_water_bar - is the increase in the DW factor of the bulk water with each consecutive layer
d_water - is the distance between water layers (in Angstroem)
After reading in the parameter file a list of all you parameters appears on the “Parameters” panel. Parameters and their limits can be modified there. With the refine check boxes you can decide whether you want to refine a certain parameter in a fit. You can set a step size and play with the < > buttons for a manual adjustment of single parameters.
The rigid body file contains info as it was used in ROD in the group definition, for each group (or rigid body) you need:
Label Number_of_atoms_in_the_group (Positions of all atoms in the group in the surface file (counting from zero)) definitions of the angular parameters (e.g. chi phi theta)
The first atom after number of atoms is the center of rotation. The atoms in a rigid body have to have the same translational parameters, otherwise the body would not be rigid.
In the bond valence file you can specify clusters of atoms that form a coordination polyhedron around a center atom for which the bond valence sum is calculated. Each bond valence cluster corresponds to one line in the input file containing:
Number_of_atoms_in_cluster center_atom_label center_atom_valence line_number_of_center_in_surface_file center_x_shift center_y_shift neighbor1_atom_label neighbor1_valence pos_of_neigh1_in_surface neigh2_x_shift neigh2_y_shift neighbor2_atom_label ... bond_valence_impact_param1 bond_valence_impact_param2
That means after the number of atoms you first list the center atom and then all neighbors. For each atom in the cluster you need five inputs: label, valence (note: +2 is just 2 (without +), -2 is -2), line in the surface file where this atom is found (integer counting from 0), x_shift and y_shift are integers which are used if the neighboring atom is in a laterally offset unit cell. At the end off each line the two so called bond valence impact parameters are given. These two values 'a' and 'b' specify an exponential function that links the impact of the bond valence on the goodness of fit parameter 'y' to the deviation of the bond valence sum from the valence of the center atom 'x'. The function is:
y = (a*x)**b
If bond valence constraints are used (the box in the lower left corner is checked), bond valence sums 'BVS' for all BV-clusters are calculated after each iteration of the fit . This gives x: x = (BVS - valence)/valence and y = (a*x)**b. The y values off all BV-clusters are summed up and the resulting sum(y) influences the goodness of fit value R:
R_new = R * (1+sum(y))
As an example the two values a = 8.26 and b = 24.16 in the first line of the example file are chosen so that if the bond valence sum is equal to the valence of the center atom (x = 0) then there is no influence on R (y = 0). If BVS is 10% of the valence of the center atom (x = 0.1) R is increased by 1% (y = 0.01). If BVS is 11% of the valence of the center atom (x = 0.11) R is increased by 10% (y = 0.1). That means with these two values a and b we can make sure bond valence has almost no influence on the fit within a tolerance of +/- 10%. Above that value the impact becomes extreme (e.g. 566% increase of R if BVS is 13% off). This model is very flexible and you can specify all kinds of shapes of the 'impact function' by choosing proper impact parameters.
Using bond valence constraints allows you to leave the limits of translational and angular parameters relatively wide open and still end up with a reasonable structure (reasonable concerning bond valence sums).
Even if you do not intend to use bond valence constraints for fitting your data, it's still worth the effort to set up a bond valence file, because it allows you to calculate the bond valences and the relevant distances between atoms with one mouse click on the 'calculate Bond Valence Sums' button.
The 'Plot CTRs' button cam be used to plot the data and the calculated structure factors. if you want to plot additional information check the corresponding boxes.
The 'Plot edens' button creates a plot of the electron density in the surface structure projected onto the surface normal. It's calculated from the input in surface-, rigid body-, and parameter files, but you can imagine it as part of a Fourier Transform of the specular rod.
with the check box you can choose whether you want to use Bond Valence constraints in the fit or not. That means according to the Bond Valence impact parameters, you set for each Bond Valence Cluster, the R value is influenced by the proximity of the calculated Bond Valence sums to the Valences of the center atoms. (see the section about the bond valence file for more info)
A click on the 'calculate Bond Valence Sums' button prints bond valence sums and bond distances for all the Bond Valence clusters to the PDS shell.
The use bulk water check box allows you to choose whether you wish to put bulk water on top of the surface sturucture. Bulk water is simulated as an layered water profile (check the e-density plot to see what it looks like). Make sure 'zwater', 'd_water', 'sig_water', and 'sig_water_bar' are specified on the parameters Panel if you want to use bulk water. 'zwater' is the height of the first water layer above the origin. 'sig_water' is a Uiso like value specifying the width of the gaussian describing the first water layer, 'dwater' is the distance between the water layeres, and 'sig_water_bar' is the increas in the width of the gaussian with each water layer.
You can spread the starting parameters randomly in parameter space if you check the use random parameters box.
The optimization routine follows the Downhill Simplex algorithm. The Simplex is by default formed by the parameter starting values and randomly distributed points in parameter space and then "crawls" along a gradient into a minimum. (In case you checked the use random parameters box, all points are random.) It's pretty robust, even if you are fitting many parameters at a time. Downhill Simplex is a compromise between the ability to probe many points in parameter space as e.g. done by genetic algorithms and the efficiency of Algorithms going straight into a minimum, as e.g. Newton- or Levenberg-Marquardt fitting. The model is optimized to produce a minimal normalized chi**2.
A nice description of Downhill Simplex (or Nelder-Mead) minimization can be found here. Parameters alpha, beta, and gamma (corresponding to alpha, rho, and gamma on wikipedia) define the stepsize of Simplex movements. Parameter delta defines the fraction of parameter space in which the initial points spanning the Simplex are randomly distributed. Ftol, Xtol, and Maxiter are convergence criteria. Ftol is the standard deviation of function values at which the routine stops. Xtol is the proximity of points in parameter space at which the routine stops. Maxiter is the maximum number of iterations allowed.
All parameters for which the refine flag on the parameter panel is set are optimized. To set parameters as equal, you can write the parameter label off one parameter into the value text boxes off other parameters. This allows you e.g. to define anisotropic thermal vibration parameters in the surface file, and nevertheless fit them as isotropic...
Intermediate results are printed continuously to the python interface structure refinement GUI status bar, so you can see how the fit procedes. If there's a new best fit the CTR- and Structure plots are updated. Figure 3 shows how points spanning the Simplex are distributed in parameter space.
Parameter statistics are computed by calculating numerical derivatives for all refined parameters with respect to their influence on all data points. This produces standard deviations of the parameters, displayed on the parameter panel. The dropdown menu allows you to print the scaled sensitivities of the selected parameter to the CTR plot window. Note that parameter standard deviations and chi**2 values depend on the total number of parameters checked. For reliable values you should, at the end of the fitting procedure, check all parameters used for the model adjustment and run the parameter statistics calculation. Parameter statistics for individual parameters can be calculated by clicking on the parameter label of a parameter in the parameter panel. Parameter statistics calculations are based on the math described in the USGS software UCODE documentation.
If you have vtk installed you may click on the Display structure button to obtain a 3D plot of the structure. If you leave this window open it should automatically update if you do any manipulations on the structure. Atom types not defined with RGB values in "tdl\modules\sxrd\atom_styles.py" appear in white. (Note: Display structure as well as Plot edens are optimized for structures with the crystal surface at the origin of the z-scale, all other cell definitions may produce strange plots.)
Use the 'Write Files' Menu to create output files.
write .cif File: produces a crude adaption to a .cif file, containing 2*2 surface unit cells. The advantage is that it includes the Debye Waller parameters and if you set up your crystal structure viewer correctly (display as a molecular structure, no translation, use U ellipsoids ...) you can directly plot the surface structure including anisotropic vibrational parameters.
write parameter File: saves your parameters to a file that can be reread in a next modeling session.
write data File: writes out CTR data plus Fcalc, Fbulk, Frough, and Fwater for further plotting in origin, excel, or whatever.
write edensity File: writes the data in the e-density plot to a file. Which is a simple and handy way to compare different structures.
write surface File: writes out the surface structure considering current parameters and rigid bodies.
save session prodces a file in which the paths to all files loaded are stored. After closing and restarting he GUI this should allow you to quickly get back to the same status by choosing load session in the Read Files menu.
The python interface structure refinement GUI Resonant Data panel allows you to analyze resonant data (RASD or RIDS or RAXR ...).