source: Tutorials/PythonScript/Scripting.htm @ 3130

Last change on this file since 3130 was 3130, checked in by toby, 5 years ago

updates for G2Scriptable and associated tutorial

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