source: Tutorials/PythonScript/Scripting.htm @ 3131

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

set as html; add link to data

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