source: Tutorials/PythonScript/Scripting.htm @ 3132

Last change on this file since 3132 was 3132, checked in by toby, 4 years ago

more minor G2scriptable changes

  • Property svn:mime-type set to text/html
File size: 17.8 KB
Line 
1<!DOCTYPE HTML PUBLIC "-//IETF//DTD HTML//EN">
2<html> <head>
3<title>Scripting GSAS-II</title>
4</head>
5
6<body>
7<h1>Running a GSAS-II Refinement From a Python Script</h1>
8In this training example we create a Python script to duplicate the
9refinement in the
10<A
11href="https://subversion.xray.aps.anl.gov/pyGSAS/Tutorials/CWCombined/Combined%20refinement.htm"
12target="_blank">GSAS-II CW Combined Refinement</A>
13tutorial. This uses the
14<A href="http://gsas-ii.readthedocs.io/en/latest/GSASIIscripts.html" target="_blank">
15GSASIIscriptable module</A>
16to perform the same refinements steps, but without use of the GSAS-II
17graphical user interface.
18
19<h2>Prerequisites</h2>
20This exercise assumes that the reader has reasonable familiarity with
21the Python language.
22Before beginning this exercise, the documentation on the
23<A href="http://gsas-ii.readthedocs.io/en/latest/GSASIIscripts.html"
24target="_blank">
25Scripting Tools in the GSAS-II GSASIIscriptable module</A> should be
26read from here: <A href="http://gsas-ii.readthedocs.io/en/latest/GSASIIscripts.html"
27target="_blank">http://gsas-ii.readthedocs.io/en/latest/GSASIIscripts.html</A>. It
28is also wise to read the <A href="https://subversion.xray.aps.anl.gov/pyGSAS/Tutorials/CWCombined/Combined%20refinement.htm"
29target="_blank">GSAS-II CW Combined Refinement</A> tutorial that this
30exercise is modeled upon, which explains why each refinement step is
31being used.
32<P>
33The exercise can be performed by placing all of the Python commands
34into a script, (which can also be
35<A href="https://subversion.xray.aps.anl.gov/trac/pyGSAS/export/3131/Tutorials/PythonScript/data/example.py">
36downloaded here</A>),
37but a more pedagogical approach will be to enter the
38commands into a Python interpreter. Use of IPython or Jupyter to run
39Python will make this a more pleasant experience.
40
41<h4>0: Load the GSASIIscriptable module</H4>
42<blockquote>
43The first step in writing any Python module is to load the modules
44that will be needed. Here we will need the <tt>os</tt> and
45<tt>sys</tt> modules from the Python standard library and the
46GSASIIscriptable from the GSAS-II code. The location where the
47GSAS-II source code is installed must usually be hard-coded into the
48script as is done in the example below. Note that a common location
49for this will be
50<TT>os.path.expanduser("~/g2conda/GSASII/")</TT>. Thus the script will
51begin with:
52
53<blockquote><textarea rows="3" cols="60" readonly>
54import os,sys
55sys.path.insert(0,'/Users/toby/software/G2/GSASII')
56import GSASIIscriptable as G2sc</textarea></blockquote>
57
58If a <tt>ImportError: No module named GSASIIscriptable</tt> error
59occurs, then the wrong directory location has been supplied for the
60GSAS-II files.
61</blockquote>
62
63<h4>1: Define some other prerequisite code</H4>
64<blockquote>
65To simplify this example, we will define the location where files will
66be written as <tt>workdir</tt> (this directory must exist, but it may
67be empty) and the location where the input files for this exercise may
68be found as <tt>datadir</tt>. (Note that this directory must have the
69the following files: "PBSO4.XRA", "INST_XRY.PRM", "PBSO4.CWN",
70"inst_d1a.prm" and "PbSO4-Wyckoff.cif". These files can be downloaded
71from
72<A href="https://subversion.xray.aps.anl.gov/pyGSAS/Tutorials/PythonScript/data/">
73https://subversion.xray.aps.anl.gov/pyGSAS/Tutorials/PythonScript/data/</A>)
74<P>
75We also define here a short function to display the weighted profile R
76factor for every histogram in a project, <tt>HistStats</tt>. This will
77be <A href="#HistStats">discussed later</A>.
78
79<blockquote><textarea rows="10" cols="70" readonly>
80workdir = "/Users/toby/Scratch/PythonScript"
81datadir = "/Users/toby/software/G2/Tutorials/PythonScript/data"
82
83def HistStats(gpx):
84    '''prints profile rfactors for all histograms'''
85    print(u"*** profile Rwp, "+os.path.split(gpx.filename)[1])
86    for hist in gpx.histograms():
87        print("\t{:20s}: {:.2f}".format(hist.name,hist.get_wR()))
88    print("")
89    gpx.save()</textarea></blockquote>
90
91</blockquote>
92
93<h4>2: Create a GSAS-II project</H4>
94<blockquote>
95The first step in creating a GSASIIscriptable script to to create or
96access a GSAS-II project, which is done with a call to
97<TT>GSASIIscriptable.G2Project()</tt>. This can be done in one of two
98ways, a call with <tt>G2sc.G2Project(filename=</tt><I>file</I><tt>)</tt> creates a new
99(empty) project, while a call with
100<tt>G2sc.G2Project(gpxfile=</tt><I>file</I><tt>)</tt> opens and
101reads an existing project (.gpx) file. Both return a Project wrapper
102object that is used to access a number of methods and variables.
103Note that <TT>GSASIIscriptable.G2Project()</tt> can read .gpx files
104written by both previous <TT>GSASIIscriptable.G2Project()</tt> runs or
105the GSAS-II GUI.
106<P>
107In this example, we create a new project by using the
108<tt>filename=</tt><I>file</I> argument. Note that this will not
109actually create a file until a save operation is done.
110   
111<blockquote><textarea rows="3" cols="70" readonly>
112# create a project with a default project name
113gpx = G2sc.G2Project(filename='PbSO4.gpx')</textarea></blockquote>
114
115</blockquote>
116<h4>3: Add Histograms and Phase to the GSAS-II project</H4>
117<blockquote>
118To add two powder diffraction datasets (histograms) to the
119project we use the <I>project</I><tt>.add_powder_histogram()</tt> method in
120the <TT>GSASIIscriptable</tt> class. The two arguments to
121<tt>.add_powder_histogram()</tt> are the powder dataset and the
122  instrument parameter file. While neither powder data file uses a
123  standard extension, the importer is able to determine the file
124  format anyway. This will not be true for all file formats.
125<UL><LI>
126  Note that <I>project</I><tt>.add_powder_histogram()</tt> returns a
127  powder histogram objects which are saved for later reference. It is
128  also possible to obtain these using <tt>gpx.histograms()</tt>, which
129  returns a list of defined histograms.
130</UL>
131<P>
132Then we add a phase to the project using
133  <I>project</I><tt>.add_phase()</tt>. This specifies a CIF containing the
134  structural information, a name for the phase and specifies that the
135  two histograms are "added" (linked) to the phase.
136  <UL><LI>
137  Note that <I>project</I><tt>.add_phase()</tt>
138  returns a phase object</li>
139  <LI>Also, the previously saved histogram objects are used in the
140  <I>project</I><tt>.add_phase()</tt> call. Note that it is
141  also possible to reference the two histgrams by their numbers in the
142  project (here <tt>histograms=[0,1]</tt>) or by the histogram names
143  (here <BR><TT>histograms=['PWDR PBSO4.XRA Bank 1', 'PWDR PBSO4.CWN
144  Bank 1']</TT>).
145  These three ways to use the <tt>histograms</tt> parameter produce
146  the same result.
147   
148<blockquote><textarea rows="10" cols="70" readonly>
149# setup step 1: add two histograms to the project
150hist1 = gpx.add_powder_histogram(os.path.join(datadir,"PBSO4.XRA"),
151                          os.path.join(datadir,"INST_XRY.PRM"))
152hist2 = gpx.add_powder_histogram(os.path.join(datadir,"PBSO4.CWN"),
153                          os.path.join(datadir,"inst_d1a.prm"))
154# setup step 2: add a phase and link it to the previous histograms
155phase0 = gpx.add_phase(os.path.join(datadir,"PbSO4-Wyckoff.cif"),
156                      phasename="PbSO4",
157                      histograms=[hist1,hist2])</textarea></blockquote>
158
159 
160</blockquote>
161<A name="ChangeCycles">
162<h4>4: Change the number of refinement cycles</H4>
163</A><blockquote>
164While this is not noted in the original tutorial, this exercise will
165  not complete properly if more variables are added to the refinement
166  without converging the refinement (or at least coming close to
167  convergence) at each refinement step. This is best accomplished by
168  increasing the number of least-squares cycles. There at present is
169  no method in the project object that allows this parameter to be
170  set, so this must be done by finding appropriate parameter
171  dictionary entry.
172<P>
173  At this point, the Python command <tt>gpx.data.keys()</tt> shows the names of the
174  entries in the GSAS-II tree; 'Controls' is the tree item
175  corresponding to the section where Least Squares cycles are
176  set. Command <tt>gpx.data['Controls'].keys()</tt> shows that all
177  values are located in an entry labeled 'data' and
178  <tt>gpx.data['Controls']['data'].keys()</tt> shows the entries in
179this section. Examination of
180<BR><tt>gpx.data['Controls']['data']['max cyc']</tt>shows the value 3,
181which is default number of cycles. Thus,
182the number of cycles is changed with this Python command:
183
184<blockquote><textarea rows="3" cols="70" readonly>
185# not in tutorial: increase # of cycles to improve convergence
186gpx.data['Controls']['data']['max cyc'] = 8 # not in API </textarea></blockquote>
187
188</blockquote>
189<h4>5: Set initial variables and refine</H4>
190<blockquote>
191
192In Step 4 of the original tutorial, the refinement is performed with
193the variables that are on by default. These variables are the two
194histogram scale factors and three background parameters for each
195histogram (8 in total). With GSASIIscriptable, the background
196refinement flags are not turned on by default, so a dictionary must be
197created to set these histogram variables. This code:
198<UL>
199<LI>defines a dictionary (refdict0) that will be used to set the number of background coefficients (not really
200needed, since the default is 3) and sets the background refinement
201flag "on".</li>
202<LI>The project is saved under a new name, so that we can later look at
203the results from each step separately.</li>
204<LI>The parameters in refdict0 are set using
205<I>project</i><tt>.do_refinements()</tt> for both histograms in the
206project and the a refinement is run.
207<LI>The weighted profile R-factor values for all histograms in the project
208are printed using the previously defined <tt>HistStats()</tt> function.
209</UL>
210
211<A name="HistStats">
212Function <tt>HistStats()</tt> works by using <tt>gpx.histograms()</tt>
213to iterate over all histograms in the project, setting <tt>hist</tt> to each histogram
214object. Class member <tt>hist.name</tt> provides the histogram name
215and method <tt>hist.get_wR()</tt> looks up the profile R-factor and
216prints them. The function also writes the final refinement results
217into the current project file. </A>
218
219<blockquote><textarea rows="5" cols="75" readonly>
220# tutorial step 4: turn on background refinement (Hist)
221refdict0 = {"set": {"Background": { "no. coeffs": 3, "refine": True }}}
222gpx.save('step4.gpx')
223gpx.do_refinements([refdict0])
224HistStats(gpx)</textarea></blockquote>
225
226The output from this will be:<blockquote><pre>
227 Hessian Levenburg-Marquardt SVD refinement on 8 variables:
228initial chi^2 9.6912e+06
229 Cycle: 0, Time: 1.88s, Chi**2: 6.7609e+05, Lambda: 0.001,  Delta: 0.93
230initial chi^2 6.7609e+05
231 Cycle: 1, Time: 1.84s, Chi**2: 6.7602e+05, Lambda: 0.001,  Delta: 0.000104
232converged
233Found 0 SVD zeros
234Read from file: /Users/toby/software/G2/Tutorials/PythonScript/step4.bak0.gpx
235Save to file  : /Users/toby/software/G2/Tutorials/PythonScript/step4.gpx
236GPX file save successful
237 Refinement results are in file: /Users/toby/software/G2/Tutorials/PythonScript/step4.lst
238 ***** Refinement successful *****
239*** profile Rwp, step4.gpx
240        PWDR PBSO4.XRA Bank 1: 40.88
241        PWDR PBSO4.CWN Bank 1: 18.65
242</pre></blockquote>Note that the Rwp values agree with what is
243 expected from the original tutorial.
244<P>
245<B>Note:</B> that there are several equivalent ways to set the histogram
246parameters using
247<tt>G2PwdrData.set_refinements()</tt>,
248<tt>G2Project.set_refinement()</tt> or
249<tt>my_project.do_refinements()</tt>, as
250<A href="http://gsas-ii.readthedocs.io/en/latest/GSASIIscripts.html#refinement-parameters">
251described here</A>. The <tt>gpx.do_refinements([refdict0])</tt>
252statement above could be replaced with:
253<blockquote><pre>
254gpx.set_refinement({"Background": { "no. coeffs": 3, "refine": True }})
255gpx.do_refinements([{}])
256</pre></blockquote>
257or
258<blockquote><pre> 
259for hist in gpx.histograms():
260    hist.set_refinements({"Background": { "no. coeffs": 3, "refine": True }})
261gpx.do_refinements([{}])
262</pre></blockquote>
263
264</blockquote>
265<h4>6: Add unit cell parameters to refinement</h4>
266<blockquote>
267
268In Step 5 of the original tutorial, the refinement is performed again,
269but the unit cell is refined in the phase.
270
271<blockquote><textarea rows="6" cols="75" readonly>
272# tutorial step 5: add unit cell refinement (Phase)
273gpx.save('step5.gpx')
274refdict1 = {"set": {"Cell": True}} # set the cell flag (for all phases)
275gpx.set_refinement(refdict1)
276gpx.do_refinements([{}])
277HistStats(gpx)</textarea></blockquote>
278
279In this case the <tt>gpx.set_refinement(refdict1)</tt> statement can
280be replaced with:<blockquote><pre>
281phase0.set_refinements({"Cell": True})
282</pre></blockquote>
283
284Note that it is also possible to combine the refinement in the current
285and previous section using
286<blockquote><pre>
287gpx.do_refinements([refdict0,refdict1])
288</pre></blockquote>
289but then the results are not saved in separate project files and it is
290not possible to see the Rwp values after each refinement.
291
292</blockquote>
293<h4>7: Add hydrostatic strain for just one histogram</h4>
294<blockquote>
295
296In Step 6 of the original tutorial, the refinement is performed again
297after adding Histogram/Phase (HAP) parameters that allow different
298lattice constants for the first histogram only.
299
300<blockquote><textarea rows="6" cols="75" readonly>
301# tutorial step 6: add Dij terms (HAP) for histogram 1 only
302gpx.save('step6.gpx')
303refdict2 = {"set": {"HStrain": True}} # set HAP parameters
304gpx.set_refinement(refdict2,phase=phase0,histogram=[hist1])
305gpx.do_refinements([{}]) # refine after setting
306HistStats(gpx)</textarea></blockquote>
307
308Here we cannot use <tt>gpx.do_refinements([refdict2])</tt> because
309that would turn on refinement of the Hstrain terms for both histograms
310(if there were more than one phase, it would be all phases and all
311histograms), but in the above the <tt>gpx.set_refinement(...)</tt>
312statement alternately can be replaced with this:
313<blockquote><pre>
314phase0.set_HAP_refinements({"HStrain": True},histograms=[hist1])
315</pre></blockquote>
316
317</blockquote>
318<h4>8: Add X-ray Sample broadening terms</h4>
319<blockquote>
320The next step in the original tutorial is to treat sample broadening
321by turning on refinement of the "Mustrain" (microstrain) and
322"Size" (Scherrer broadening) terms using an isotropic (single-term)
323model. As described in the
324<A
325href="http://gsas-ii.readthedocs.io/en/latest/GSASIIscripts.html#histogram-and-phase-parameters">
326documentation for Histogram/Phase parameters</A>, type must always be
327specfied even as in this case where it is not being changed from the
328existing default. Again, since these parameters are being set only for
329one histogram, either <tt>phase0.set_HAP_refinements()</tt> or
330<tt>gpx.set_refinement()</tt> must be used.
331
332<blockquote><textarea rows="8" cols="75" readonly>
333# tutorial step 7: add size & strain broadening (HAP) for histogram 1 only
334gpx.save('step7.gpx')
335refdict2 = {"set": {"Mustrain": {"type":"isotropic","refine":True},
336                    "Size":{"type":"isotropic","refine":True},
337                    }}
338gpx.set_refinement(refdict2,phase=phase0,histogram=[hist1])
339gpx.do_refinements([{}]) # refine after setting
340HistStats(gpx)</textarea></blockquote>
341
342</blockquote>
343<h4>9: Add Structural and Sample Displacement Parameters</h4>
344<blockquote>
345The step 8 in the original tutorial is to treat sample displacement
346for each histogram/phase. These parameters are different because of
347the differing diffractometer geometries. We also add refinement of
348sample parameters using <tt>phase0.set_refinements()</tt> to set the
349"X" (coordinate) and "U" (displacement) refinement flags for all
350atoms. The original tutorial
351calls for the diffractometer radius to be changed. This parameter
352cannot be set from any GSASIIscriptable routines, but following a
353similar
354<A href="#ChangeCycles">process, as before</A>, the location for this
355setting can be located in the histogram's 'Sample Parameters' section
356(<tt>hist2.data['Sample Parameters']['Gonio. radius']</tt>).
357
358<blockquote><textarea rows="9" cols="75" readonly>
359# tutorial step 8: add sample parameters & set radius (Hist); refine atom parameters (phase)
360gpx.save('step8.gpx')
361hist1.set_refinements({'Sample Parameters': ['Shift']})
362hist2.set_refinements({'Sample Parameters': ['DisplaceX', 'DisplaceY']})
363hist2.data['Sample Parameters']['Gonio. radius'] = 650. # not in API
364phase0.set_refinements({"Atoms":{"all":"XU"}})
365gpx.do_refinements([{}]) # refine after setting
366HistStats(gpx)</textarea></blockquote>
367
368Note that the settings provided in the
369<tt>phase0.set_refinements()</tt> statement could have been done using
370the <tt>gpx.do_refinements()</tt> call, but not those in the
371<tt>hist</tt><I>X</I><tt>.set_refinements()</tt> calls, since they
372must be different for each histogram.
373
374</blockquote>
375<h4>10: Change Data Limits; Vary Gaussian Profile Terms</h4>
376<blockquote>
377The final step in the original tutorial is to trim the range of data
378used in the refinement to exclude data where no reflections occur and
379where the peaks are cut off at high angle. Also, additional parameters
380are refined here because the supplied instrumental profile parameters
381are not very accurate descriptions for the datasets that are used
382here. It is not possible to refine the Lorentzian x-ray instrumental
383broadening terms, since this broadening is treated as sample
384broadening. The Lorentzian neutron broadening is negligible.
385
386<blockquote><textarea rows="6" cols="75" readonly>
387# tutorial step 9: change data limits & inst. parm refinements (Hist)
388gpx.save('step9.gpx')
389hist1.set_refinements({'Limits': [16.,158.4]})
390hist2.set_refinements({'Limits': [19.,153.]})
391gpx.do_refinements([{"set": {'Instrument Parameters': ['U', 'V', 'W']}}])
392HistStats(gpx)</textarea></blockquote>
393
394Note that the final <tt>gpx.do_refinements()</tt> statement can be
395replaced with calls to
396<tt>hist</tt><I>X</I><tt>.set_refinements()</tt> or
397<tt>gpx.set_refinement()</tt>, such as
398<blockquote><pre>
399hist1.set_refinements({'Instrument Parameters': ['U', 'V', 'W']})
400hist2.set_refinements({'Instrument Parameters': ['U', 'V', 'W']})
401gpx.do_refinements([{}])
402</pre></blockquote>
403
404
405</blockquote><hr>
406<address></address>
407<!-- hhmts start -->Last modified: Thu Oct 12 10:45:51 CDT 2017 <!-- hhmts end -->
408</body> </html>
Note: See TracBrowser for help on using the repository browser.