1 | # -*- coding: utf-8 -*- |
---|
2 | ########### SVN repository information ################### |
---|
3 | # $Date: 2023-02-05 15:40:17 +0000 (Sun, 05 Feb 2023) $ |
---|
4 | # $Author: vondreele $ |
---|
5 | # $Revision: 5494 $ |
---|
6 | # $URL: trunk/GSASIIstrIO.py $ |
---|
7 | # $Id: GSASIIstrIO.py 5494 2023-02-05 15:40:17Z vondreele $ |
---|
8 | ########### SVN repository information ################### |
---|
9 | ''' |
---|
10 | *GSASIIstrIO: structure I/O routines* |
---|
11 | ------------------------------------- |
---|
12 | |
---|
13 | Contains routines for reading from GPX files and printing to the .LST file. |
---|
14 | Used for refinements and in G2scriptable. |
---|
15 | |
---|
16 | Should not contain any wxpython references as this should be able to be used |
---|
17 | in non-GUI settings. |
---|
18 | |
---|
19 | ''' |
---|
20 | from __future__ import division, print_function |
---|
21 | import platform |
---|
22 | import re |
---|
23 | import os |
---|
24 | import os.path as ospath |
---|
25 | import time |
---|
26 | import math |
---|
27 | import copy |
---|
28 | if '2' in platform.python_version_tuple()[0]: |
---|
29 | import cPickle |
---|
30 | else: |
---|
31 | import pickle as cPickle |
---|
32 | import numpy as np |
---|
33 | import numpy.ma as ma |
---|
34 | import GSASIIpath |
---|
35 | GSASIIpath.SetVersionNumber("$Revision: 5494 $") |
---|
36 | import GSASIIElem as G2el |
---|
37 | import GSASIIlattice as G2lat |
---|
38 | import GSASIIspc as G2spc |
---|
39 | import GSASIIobj as G2obj |
---|
40 | import GSASIImapvars as G2mv |
---|
41 | import GSASIImath as G2mth |
---|
42 | import GSASIIstrMath as G2stMth |
---|
43 | import GSASIIfiles as G2fil |
---|
44 | import GSASIIpy3 as G2py3 |
---|
45 | |
---|
46 | sind = lambda x: np.sin(x*np.pi/180.) |
---|
47 | cosd = lambda x: np.cos(x*np.pi/180.) |
---|
48 | tand = lambda x: np.tan(x*np.pi/180.) |
---|
49 | asind = lambda x: 180.*np.arcsin(x)/np.pi |
---|
50 | acosd = lambda x: 180.*np.arccos(x)/np.pi |
---|
51 | atan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi |
---|
52 | |
---|
53 | ateln2 = 8.0*math.log(2.0) |
---|
54 | |
---|
55 | #=============================================================================== |
---|
56 | # Support for GPX file reading |
---|
57 | #=============================================================================== |
---|
58 | def cPickleLoad(fp): |
---|
59 | if '2' in platform.python_version_tuple()[0]: |
---|
60 | return cPickle.load(fp) |
---|
61 | else: |
---|
62 | return cPickle.load(fp,encoding='latin-1') |
---|
63 | |
---|
64 | gpxIndex = {}; gpxNamelist = []; gpxSize = -1 |
---|
65 | '''Global variables used in :func:`IndexGPX` to see if file has changed |
---|
66 | (gpxSize) and to index where to find each 1st-level tree item in the file. |
---|
67 | ''' |
---|
68 | |
---|
69 | def GetFullGPX(GPXfile): |
---|
70 | ''' Returns complete contents of GSASII gpx file. |
---|
71 | Used in :func:`GSASIIscriptable.LoadDictFromProjFile`. |
---|
72 | |
---|
73 | :param str GPXfile: full .gpx file name |
---|
74 | :returns: Project,nameList, where |
---|
75 | |
---|
76 | * Project (dict) is a representation of gpx file following the GSAS-II |
---|
77 | tree structure for each item: key = tree name (e.g. 'Controls', |
---|
78 | 'Restraints', etc.), data is dict |
---|
79 | * nameList (list) has names of main tree entries & subentries used to reconstruct project file |
---|
80 | ''' |
---|
81 | return IndexGPX(GPXfile,read=True) |
---|
82 | |
---|
83 | def IndexGPX(GPXfile,read=False): |
---|
84 | '''Create an index to a GPX file, optionally the file into memory. |
---|
85 | The byte size of the GPX file is saved. If this routine is called |
---|
86 | again, and if this size does not change, indexing is not repeated |
---|
87 | since it is assumed the file has not changed (this can be overriden |
---|
88 | by setting read=True). |
---|
89 | |
---|
90 | :param str GPXfile: full .gpx file name |
---|
91 | :returns: Project,nameList if read=, where |
---|
92 | |
---|
93 | * Project (dict) is a representation of gpx file following the GSAS-II |
---|
94 | tree structure for each item: key = tree name (e.g. 'Controls', |
---|
95 | 'Restraints', etc.), data is dict |
---|
96 | * nameList (list) has names of main tree entries & subentries used to reconstruct project file |
---|
97 | ''' |
---|
98 | global gpxSize |
---|
99 | if gpxSize == os.path.getsize(GPXfile) and not read: |
---|
100 | return |
---|
101 | global gpxIndex |
---|
102 | gpxIndex = {} |
---|
103 | global gpxNamelist |
---|
104 | gpxNamelist = [] |
---|
105 | if GSASIIpath.GetConfigValue('debug'): print("DBG: Indexing GPX file") |
---|
106 | gpxSize = os.path.getsize(GPXfile) |
---|
107 | fp = open(GPXfile,'rb') |
---|
108 | Project = {} |
---|
109 | try: |
---|
110 | while True: |
---|
111 | pos = fp.tell() |
---|
112 | data = cPickleLoad(fp) |
---|
113 | datum = data[0] |
---|
114 | gpxIndex[datum[0]] = pos |
---|
115 | if read: Project[datum[0]] = {'data':datum[1]} |
---|
116 | gpxNamelist.append([datum[0],]) |
---|
117 | for datus in data[1:]: |
---|
118 | if read: Project[datum[0]][datus[0]] = datus[1] |
---|
119 | gpxNamelist[-1].append(datus[0]) |
---|
120 | # print('project load successful') |
---|
121 | except EOFError: |
---|
122 | pass |
---|
123 | except Exception as msg: |
---|
124 | G2fil.G2Print('Read Error:',msg) |
---|
125 | raise Exception("Error reading file "+str(GPXfile)+". This is not a GSAS-II .gpx file") |
---|
126 | finally: |
---|
127 | fp.close() |
---|
128 | if read: return Project,gpxNamelist |
---|
129 | |
---|
130 | def GetControls(GPXfile): |
---|
131 | ''' Returns dictionary of control items found in GSASII gpx file |
---|
132 | |
---|
133 | :param str GPXfile: full .gpx file name |
---|
134 | :return: dictionary of control items |
---|
135 | ''' |
---|
136 | Controls = copy.deepcopy(G2obj.DefaultControls) |
---|
137 | IndexGPX(GPXfile) |
---|
138 | pos = gpxIndex.get('Controls') |
---|
139 | if pos is None: |
---|
140 | G2fil.G2Print('Warning: Controls not found in gpx file {}'.format(GPXfile)) |
---|
141 | return Controls |
---|
142 | fp = open(GPXfile,'rb') |
---|
143 | fp.seek(pos) |
---|
144 | datum = cPickleLoad(fp)[0] |
---|
145 | fp.close() |
---|
146 | Controls.update(datum[1]) |
---|
147 | return Controls |
---|
148 | |
---|
149 | def ReadConstraints(GPXfile, seqHist=None): |
---|
150 | '''Read the constraints from the GPX file and interpret them |
---|
151 | |
---|
152 | called in :func:`ReadCheckConstraints`, :func:`GSASIIstrMain.Refine` |
---|
153 | and :func:`GSASIIstrMain.SeqRefine`. |
---|
154 | ''' |
---|
155 | IndexGPX(GPXfile) |
---|
156 | fl = open(GPXfile,'rb') |
---|
157 | pos = gpxIndex.get('Constraints') |
---|
158 | if pos is None: |
---|
159 | raise Exception("No constraints in GPX file") |
---|
160 | fl.seek(pos) |
---|
161 | ConstraintsItem = cPickleLoad(fl)[0] |
---|
162 | seqmode = 'use-all' |
---|
163 | if seqHist is not None: |
---|
164 | seqmode = ConstraintsItem[1].get('_seqmode','wildcards-only') |
---|
165 | fl.close() |
---|
166 | constList = [] |
---|
167 | for item in ConstraintsItem[1]: |
---|
168 | if item.startswith('_'): continue |
---|
169 | constList += ConstraintsItem[1][item] |
---|
170 | constrDict,fixedList,ignored = G2mv.ProcessConstraints(constList,seqmode,seqHist) |
---|
171 | #if ignored: |
---|
172 | # G2fil.G2Print ('Warning: {} Constraints were rejected. Was a constrained phase, histogram or atom deleted?'.format(ignored)) |
---|
173 | return constrDict,fixedList |
---|
174 | |
---|
175 | def ReadCheckConstraints(GPXfile, seqHist=None,Histograms=None,Phases=None): |
---|
176 | '''Load constraints and related info and return any error or warning messages |
---|
177 | This is done from the GPX file rather than the tree. |
---|
178 | |
---|
179 | :param str GPXfile: specifies the path to a .gpx file. |
---|
180 | :param str seqHist: specifies a histogram to be loaded for |
---|
181 | a sequential refinement. If None (default) all are loaded. |
---|
182 | :param dict Histograms: output from :func:`GetUsedHistogramsAndPhases`, |
---|
183 | can optionally be supplied to save time for sequential refinements |
---|
184 | :param dict Phases: output from :func:`GetUsedHistogramsAndPhases`, can |
---|
185 | optionally be supplied to save time for sequential refinements |
---|
186 | ''' |
---|
187 | G2mv.InitVars() # init constraints |
---|
188 | # get variables |
---|
189 | if Histograms is None or Phases is None: |
---|
190 | Histograms,Phases = GetUsedHistogramsAndPhases(GPXfile) |
---|
191 | if not Phases: |
---|
192 | return 'Error: No phases or no histograms in phases!','' |
---|
193 | if not Histograms: |
---|
194 | return 'Error: no diffraction data','' |
---|
195 | if seqHist: |
---|
196 | Histograms = {seqHist:Histograms[seqHist]} # sequential fit: only need one histogram |
---|
197 | hId = Histograms[seqHist]['hId'] |
---|
198 | else: |
---|
199 | hId = None |
---|
200 | constrDict,fixedList = ReadConstraints(GPXfile, hId) # load user constraints from file (uses ProcessConstraints) |
---|
201 | parmDict = {} |
---|
202 | # generate symmetry constraints to check for conflicts |
---|
203 | rigidbodyDict = GetRigidBodies(GPXfile) |
---|
204 | rbIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[],'Spin':[]}) |
---|
205 | rbVary,rbDict = GetRigidBodyModels(rigidbodyDict,Print=False) |
---|
206 | parmDict.update(rbDict) |
---|
207 | (Natoms, atomIndx, phaseVary,phaseDict, pawleyLookup,FFtables,EFtables,BLtables,MFtables,maxSSwave) = \ |
---|
208 | GetPhaseData(Phases,RestraintDict=None,seqHistName=seqHist,rbIds=rbIds,Print=False) # generates atom symmetry constraints |
---|
209 | parmDict.update(phaseDict) |
---|
210 | hapVary,hapDict,controlDict = GetHistogramPhaseData(Phases,Histograms,Print=False,resetRefList=False) |
---|
211 | parmDict.update(hapDict) |
---|
212 | histVary,histDict,controlDict = GetHistogramData(Histograms,Print=False) |
---|
213 | parmDict.update(histDict) |
---|
214 | varyList = rbVary+phaseVary+hapVary+histVary |
---|
215 | msg = G2mv.EvaluateMultipliers(constrDict,phaseDict,hapDict,histDict) |
---|
216 | if msg: |
---|
217 | return 'Unable to interpret multiplier(s): '+msg,'' |
---|
218 | errmsg,warnmsg,groups,parmlist = G2mv.GenerateConstraints(varyList,constrDict,fixedList,parmDict, |
---|
219 | seqHistNum=hId) |
---|
220 | G2mv.Map2Dict(parmDict,varyList) # changes varyList |
---|
221 | return errmsg, warnmsg |
---|
222 | |
---|
223 | def makeTwinFrConstr(Phases,Histograms,hapVary): |
---|
224 | TwConstr = [] |
---|
225 | TwFixed = [] |
---|
226 | for Phase in Phases: |
---|
227 | pId = Phases[Phase]['pId'] |
---|
228 | for Histogram in Phases[Phase]['Histograms']: |
---|
229 | try: |
---|
230 | hId = Histograms[Histogram]['hId'] |
---|
231 | phfx = '%d:%d:'%(pId,hId) |
---|
232 | if phfx+'TwinFr:0' in hapVary: |
---|
233 | TwFixed.append('1.0') #constraint value |
---|
234 | nTwin = len(Phases[Phase]['Histograms'][Histogram]['Twins']) |
---|
235 | TwConstr.append({phfx+'TwinFr:'+str(i):'1.0' for i in range(nTwin)}) |
---|
236 | except KeyError: #unused histograms? |
---|
237 | pass |
---|
238 | return TwConstr,TwFixed |
---|
239 | |
---|
240 | def GetRestraints(GPXfile): |
---|
241 | '''Read the restraints from the GPX file. |
---|
242 | Throws an exception if not found in the .GPX file |
---|
243 | ''' |
---|
244 | IndexGPX(GPXfile) |
---|
245 | fl = open(GPXfile,'rb') |
---|
246 | pos = gpxIndex.get('Restraints') |
---|
247 | if pos is None: |
---|
248 | raise Exception("No Restraints in GPX file") |
---|
249 | fl.seek(pos) |
---|
250 | datum = cPickleLoad(fl)[0] |
---|
251 | fl.close() |
---|
252 | return datum[1] |
---|
253 | |
---|
254 | def GetRigidBodies(GPXfile): |
---|
255 | '''Read the rigid body models from the GPX file |
---|
256 | ''' |
---|
257 | IndexGPX(GPXfile) |
---|
258 | fl = open(GPXfile,'rb') |
---|
259 | pos = gpxIndex.get('Rigid bodies') |
---|
260 | if pos is None: |
---|
261 | raise Exception("No Rigid bodies in GPX file") |
---|
262 | fl.seek(pos) |
---|
263 | datum = cPickleLoad(fl)[0] |
---|
264 | fl.close() |
---|
265 | return datum[1] |
---|
266 | |
---|
267 | def GetFprime(controlDict,Histograms): |
---|
268 | 'Needs a doc string' |
---|
269 | FFtables = controlDict['FFtables'] |
---|
270 | if not FFtables: |
---|
271 | return |
---|
272 | histoList = list(Histograms.keys()) |
---|
273 | histoList.sort() |
---|
274 | for histogram in histoList: |
---|
275 | if histogram[:4] in ['PWDR','HKLF']: |
---|
276 | Histogram = Histograms[histogram] |
---|
277 | hId = Histogram['hId'] |
---|
278 | hfx = ':%d:'%(hId) |
---|
279 | if 'E' in controlDict[hfx+'histType']: |
---|
280 | for El in FFtables: |
---|
281 | FFtables[El][hfx+'FP'] = 0. |
---|
282 | FFtables[El][hfx+'FPP'] = 0. |
---|
283 | elif 'X' in controlDict[hfx+'histType']: |
---|
284 | keV = controlDict[hfx+'keV'] |
---|
285 | for El in FFtables: |
---|
286 | Orbs = G2el.GetXsectionCoeff(El.split('+')[0].split('-')[0]) |
---|
287 | FP,FPP,Mu = G2el.FPcalc(Orbs, keV) |
---|
288 | FFtables[El][hfx+'FP'] = FP |
---|
289 | FFtables[El][hfx+'FPP'] = FPP |
---|
290 | |
---|
291 | def PrintFprime(FFtables,pfx,pFile): |
---|
292 | pFile.write('\n Resonant form factors:(ref: D.T. Cromer & D.A. Liberman (1981), Acta Cryst. A37, 267-268.)\n') |
---|
293 | Elstr = ' Element:' |
---|
294 | FPstr = " f' :" |
---|
295 | FPPstr = ' f" :' |
---|
296 | for El in FFtables: |
---|
297 | Elstr += ' %8s'%(El) |
---|
298 | FPstr += ' %8.3f'%(FFtables[El][pfx+'FP']) |
---|
299 | FPPstr += ' %8.3f'%(FFtables[El][pfx+'FPP']) |
---|
300 | pFile.write(Elstr+'\n') |
---|
301 | pFile.write(FPstr+'\n') |
---|
302 | pFile.write(FPPstr+'\n') |
---|
303 | |
---|
304 | def PrintBlength(BLtables,wave,pFile): |
---|
305 | pFile.write('\n Resonant neutron scattering lengths:\n') |
---|
306 | Elstr = ' Element:' |
---|
307 | FPstr = " b' :" |
---|
308 | FPPstr = ' b" :' |
---|
309 | for El in BLtables: |
---|
310 | BP,BPP = G2el.BlenResCW([El,],BLtables,wave) |
---|
311 | Elstr += ' %8s'%(El) |
---|
312 | FPstr += ' %8.3f'%(BP) |
---|
313 | FPPstr += ' %8.3f'%(BPP) |
---|
314 | pFile.write(Elstr+'\n') |
---|
315 | pFile.write(FPstr+'\n') |
---|
316 | pFile.write(FPPstr+'\n') |
---|
317 | |
---|
318 | def PrintISOmodes(pFile,Phases,parmDict,sigDict): |
---|
319 | '''Prints the values for the ISODISTORT modes into the project's |
---|
320 | .lst file after a refinement. |
---|
321 | ''' |
---|
322 | for phase in Phases: |
---|
323 | if 'ISODISTORT' not in Phases[phase]: continue |
---|
324 | data = Phases[phase] |
---|
325 | ISO = data['ISODISTORT'] |
---|
326 | atNames = [atom[0] for atom in data['Atoms']] |
---|
327 | |
---|
328 | if 'G2VarList' in ISO: |
---|
329 | deltaList = [] |
---|
330 | notfound = [] |
---|
331 | for gv,Ilbl in zip(ISO['G2VarList'],ISO['IsoVarList']): |
---|
332 | dvar = gv.varname() |
---|
333 | var = dvar.replace('::dA','::A') |
---|
334 | atnum = atNames.index(Ilbl[:Ilbl.rfind('_')]) |
---|
335 | v = Ilbl[Ilbl.rfind('_')+1:] |
---|
336 | pval = ISO['G2parentCoords'][atnum][['dx','dy','dz'].index(v)] |
---|
337 | if var in parmDict: |
---|
338 | cval = parmDict[var] |
---|
339 | else: |
---|
340 | notfound.append(var) |
---|
341 | continue |
---|
342 | deltaList.append(cval-pval) |
---|
343 | if notfound: |
---|
344 | msg = 'PrintISOmodes warning: Atom parameters ' |
---|
345 | for i,v in enumerate(notfound): |
---|
346 | if i == len(notfound)-1: |
---|
347 | msg += ' & ' |
---|
348 | elif i != 0: |
---|
349 | msg += ', ' |
---|
350 | msg += v |
---|
351 | print(msg,'not found') |
---|
352 | print(' skipping computation for modes:') |
---|
353 | for i,j in zip(ISO['IsoModeList'],ISO['G2ModeList']): |
---|
354 | print(' ',i,'({})'.format(j)) |
---|
355 | continue |
---|
356 | modeVals = np.inner(ISO['Var2ModeMatrix'],deltaList) |
---|
357 | |
---|
358 | pFile.write('\n ISODISTORT Displacive Modes for phase {}\n'.format(data['General'].get('Name',''))) |
---|
359 | l = str(max([len(i) for i in ISO['IsoModeList']])+3) |
---|
360 | fmt = ' {:'+l+'}{}' |
---|
361 | for varid,[var,val,norm,G2mode] in enumerate(zip( |
---|
362 | ISO['IsoModeList'],modeVals,ISO['NormList'],ISO['G2ModeList'] )): |
---|
363 | try: |
---|
364 | value = G2mth.ValEsd(val/norm,-0.001) |
---|
365 | item = str(G2mode).replace('::','::nv-') |
---|
366 | if item in sigDict: |
---|
367 | ISO['modeDispl'][varid] = val/norm |
---|
368 | value = G2mth.ValEsd(val/norm,sigDict[item]/norm) |
---|
369 | except TypeError: |
---|
370 | value = '?' |
---|
371 | pFile.write(fmt.format(var,value)+'\n') |
---|
372 | |
---|
373 | if 'G2OccVarList' in ISO: #untested - probably wrong |
---|
374 | deltaOccList = [] |
---|
375 | notfound = [] |
---|
376 | for gv,Ilbl in zip(ISO['G2OccVarList'],ISO['OccVarList']): |
---|
377 | var = gv.varname() |
---|
378 | albl = Ilbl[:Ilbl.rfind('_')] |
---|
379 | pval = ISO['BaseOcc'][albl] |
---|
380 | if var in parmDict: |
---|
381 | cval = parmDict[var] |
---|
382 | else: |
---|
383 | notfound.append(var) |
---|
384 | continue |
---|
385 | deltaOccList.append(cval-pval) |
---|
386 | if notfound: |
---|
387 | msg = 'PrintISOmodes warning: Atom parameters ' |
---|
388 | for i,v in enumerate(notfound): |
---|
389 | if i == len(notfound)-1: |
---|
390 | msg += ' & ' |
---|
391 | elif i != 0: |
---|
392 | msg += ', ' |
---|
393 | msg += v |
---|
394 | print(msg,'not found') |
---|
395 | print(' skipping computation for modes:') |
---|
396 | for i,j in zip(ISO['OccVarList'],ISO['G2OccVarList']): |
---|
397 | print(' ',i,'({})'.format(j)) |
---|
398 | continue |
---|
399 | modeOccVals = np.inner(ISO['Var2OccMatrix'],deltaOccList) |
---|
400 | |
---|
401 | pFile.write('\n ISODISTORT Occupancy Modes for phase {}\n'.format(data['General'].get('Name',''))) |
---|
402 | l = str(max([len(i) for i in ISO['OccModeList']])+3) |
---|
403 | fmt = ' {:'+l+'}{}' |
---|
404 | for var,val,norm,G2mode in zip( |
---|
405 | ISO['OccModeList'],modeOccVals,ISO['OccNormList'],ISO['G2OccModeList'] ): |
---|
406 | try: |
---|
407 | value = G2py3.FormatSigFigs(val/norm) |
---|
408 | if str(G2mode) in sigDict: |
---|
409 | value = G2mth.ValEsd(val/norm,sigDict[str(G2mode)]/norm) |
---|
410 | except TypeError: |
---|
411 | value = '?' |
---|
412 | pFile.write(fmt.format(var,value)+'\n') |
---|
413 | |
---|
414 | def PrintIndependentVars(parmDict,varyList,sigDict,PrintAll=False,pFile=None): |
---|
415 | '''Print the values and uncertainties on the independent parameters''' |
---|
416 | printlist = [] |
---|
417 | mvs = G2mv.GetIndependentVars() |
---|
418 | for i,name in enumerate(mvs): |
---|
419 | if PrintAll or name in varyList: |
---|
420 | sig = sigDict.get(name) |
---|
421 | printlist.append([name,parmDict[name],sig]) |
---|
422 | if len(printlist) == 0: return |
---|
423 | s1 = '' |
---|
424 | s2 = '' |
---|
425 | s3 = '' |
---|
426 | pFile.write(130*'-'+'\n') |
---|
427 | pFile.write("Parameters generated by constraints\n") |
---|
428 | printlist.append(3*[None]) |
---|
429 | for name,val,esd in printlist: |
---|
430 | if len(s1) > 120 or name is None: |
---|
431 | pFile.write(''+'\n') |
---|
432 | pFile.write(s1+'\n') |
---|
433 | pFile.write(s2+'\n') |
---|
434 | pFile.write(s3+'\n') |
---|
435 | s1 = '' |
---|
436 | if name is None: break |
---|
437 | if s1 == "": |
---|
438 | s1 = ' name :' |
---|
439 | s2 = ' value :' |
---|
440 | s3 = ' sig : ' |
---|
441 | wdt = len(name)+1 |
---|
442 | s1 += ('%15s' % (name)).rjust(wdt) |
---|
443 | s2 += ('%15.5f' % (val)).center(wdt) |
---|
444 | if esd is None: |
---|
445 | s3 += ('%15s' % ('n/a')).center(wdt) |
---|
446 | else: |
---|
447 | s3 += fmtESD(name,sigDict,'f',15,5).center(wdt) |
---|
448 | |
---|
449 | def GetPhaseNames(GPXfile): |
---|
450 | ''' Returns a list of phase names found under 'Phases' in GSASII gpx file |
---|
451 | |
---|
452 | :param str GPXfile: full .gpx file name |
---|
453 | :return: list of phase names |
---|
454 | ''' |
---|
455 | IndexGPX(GPXfile) |
---|
456 | fl = open(GPXfile,'rb') |
---|
457 | pos = gpxIndex.get('Phases') |
---|
458 | if pos is None: |
---|
459 | raise Exception("No Phases in GPX file") |
---|
460 | fl.seek(pos) |
---|
461 | data = cPickleLoad(fl) |
---|
462 | fl.close() |
---|
463 | return [datus[0] for datus in data[1:]] |
---|
464 | |
---|
465 | def GetAllPhaseData(GPXfile,PhaseName): |
---|
466 | ''' Returns the entire dictionary for PhaseName from GSASII gpx file |
---|
467 | |
---|
468 | :param str GPXfile: full .gpx file name |
---|
469 | :param str PhaseName: phase name |
---|
470 | :return: phase dictionary or None if PhaseName is not present |
---|
471 | ''' |
---|
472 | IndexGPX(GPXfile) |
---|
473 | fl = open(GPXfile,'rb') |
---|
474 | pos = gpxIndex.get('Phases') |
---|
475 | if pos is None: |
---|
476 | raise Exception("No Phases in GPX file") |
---|
477 | fl.seek(pos) |
---|
478 | data = cPickleLoad(fl) |
---|
479 | fl.close() |
---|
480 | |
---|
481 | for datus in data[1:]: |
---|
482 | if datus[0] == PhaseName: |
---|
483 | return datus[1] |
---|
484 | |
---|
485 | def GetHistograms(GPXfile,hNames): |
---|
486 | """ Returns a dictionary of histograms found in GSASII gpx file |
---|
487 | |
---|
488 | :param str GPXfile: full .gpx file name |
---|
489 | :param str hNames: list of histogram names |
---|
490 | :return: dictionary of histograms (types = PWDR & HKLF) |
---|
491 | |
---|
492 | """ |
---|
493 | IndexGPX(GPXfile) |
---|
494 | fl = open(GPXfile,'rb') |
---|
495 | Histograms = {} |
---|
496 | for hist in hNames: |
---|
497 | pos = gpxIndex.get(hist) |
---|
498 | if pos is None: |
---|
499 | raise Exception("Histogram {} not found in GPX file".format(hist)) |
---|
500 | fl.seek(pos) |
---|
501 | data = cPickleLoad(fl) |
---|
502 | datum = data[0] |
---|
503 | if 'PWDR' in hist[:4]: |
---|
504 | PWDRdata = {} |
---|
505 | PWDRdata.update(datum[1][0]) #weight factor |
---|
506 | PWDRdata['Data'] = ma.array(ma.getdata(datum[1][1])) #masked powder data arrays/clear previous masks |
---|
507 | PWDRdata[data[2][0]] = data[2][1] #Limits & excluded regions (if any) |
---|
508 | PWDRdata[data[3][0]] = data[3][1] #Background |
---|
509 | PWDRdata[data[4][0]] = data[4][1] #Instrument parameters |
---|
510 | PWDRdata[data[5][0]] = data[5][1] #Sample parameters |
---|
511 | try: |
---|
512 | PWDRdata[data[9][0]] = data[9][1] #Reflection lists might be missing |
---|
513 | except IndexError: |
---|
514 | PWDRdata['Reflection Lists'] = {} |
---|
515 | PWDRdata['Residuals'] = {} |
---|
516 | Histograms[hist] = PWDRdata |
---|
517 | elif 'HKLF' in hist[:4]: |
---|
518 | HKLFdata = {} |
---|
519 | HKLFdata.update(datum[1][0]) #weight factor |
---|
520 | #patch |
---|
521 | if 'list' in str(type(datum[1][1])): |
---|
522 | #if isinstance(datum[1][1],list): |
---|
523 | RefData = {'RefList':[],'FF':{}} |
---|
524 | for ref in datum[1][1]: |
---|
525 | RefData['RefList'].append(ref[:11]+[ref[13],]) |
---|
526 | RefData['RefList'] = np.array(RefData['RefList']) |
---|
527 | datum[1][1] = RefData |
---|
528 | #end patch |
---|
529 | datum[1][1]['FF'] = {} |
---|
530 | HKLFdata['Data'] = datum[1][1] |
---|
531 | HKLFdata['Instrument Parameters'] = dict(data)['Instrument Parameters'] |
---|
532 | HKLFdata['Reflection Lists'] = None |
---|
533 | HKLFdata['Residuals'] = {} |
---|
534 | Histograms[hist] = HKLFdata |
---|
535 | fl.close() |
---|
536 | return Histograms |
---|
537 | |
---|
538 | def GetHistogramNames(GPXfile,hTypes): |
---|
539 | """ Returns a list of histogram names found in a GSAS-II .gpx file that |
---|
540 | match specifed histogram types. Names are returned in the order they |
---|
541 | appear in the file. |
---|
542 | |
---|
543 | :param str GPXfile: full .gpx file name |
---|
544 | :param str hTypes: list of histogram types |
---|
545 | :return: list of histogram names (types = PWDR & HKLF) |
---|
546 | |
---|
547 | """ |
---|
548 | IndexGPX(GPXfile) |
---|
549 | return [n[0] for n in gpxNamelist if n[0][:4] in hTypes] |
---|
550 | |
---|
551 | def GetUsedHistogramsAndPhases(GPXfile): |
---|
552 | ''' Returns all histograms that are found in any phase |
---|
553 | and any phase that uses a histogram. This also |
---|
554 | assigns numbers to used phases and histograms by the |
---|
555 | order they appear in the file. |
---|
556 | |
---|
557 | :param str GPXfile: full .gpx file name |
---|
558 | :returns: (Histograms,Phases) |
---|
559 | |
---|
560 | * Histograms = dictionary of histograms as {name:data,...} |
---|
561 | * Phases = dictionary of phases that use histograms |
---|
562 | |
---|
563 | ''' |
---|
564 | phaseNames = GetPhaseNames(GPXfile) |
---|
565 | histoList = GetHistogramNames(GPXfile,['PWDR','HKLF']) |
---|
566 | allHistograms = GetHistograms(GPXfile,histoList) |
---|
567 | phaseData = {} |
---|
568 | for name in phaseNames: |
---|
569 | phaseData[name] = GetAllPhaseData(GPXfile,name) |
---|
570 | Histograms = {} |
---|
571 | Phases = {} |
---|
572 | for phase in phaseData: |
---|
573 | Phase = phaseData[phase] |
---|
574 | if Phase['General']['Type'] == 'faulted': continue #don't use faulted phases! |
---|
575 | if Phase['Histograms']: |
---|
576 | for hist in Phase['Histograms']: |
---|
577 | if 'Use' not in Phase['Histograms'][hist]: #patch |
---|
578 | Phase['Histograms'][hist]['Use'] = True |
---|
579 | if Phase['Histograms'][hist]['Use'] and phase not in Phases: |
---|
580 | pId = phaseNames.index(phase) |
---|
581 | Phase['pId'] = pId |
---|
582 | Phases[phase] = Phase |
---|
583 | if hist not in Histograms and Phase['Histograms'][hist]['Use']: |
---|
584 | try: |
---|
585 | Histograms[hist] = allHistograms[hist] |
---|
586 | hId = histoList.index(hist) |
---|
587 | Histograms[hist]['hId'] = hId |
---|
588 | except KeyError: # would happen if a referenced histogram were |
---|
589 | # renamed or deleted |
---|
590 | G2fil.G2Print('Warning: For phase "'+phase+ |
---|
591 | '" unresolved reference to histogram "'+hist+'"') |
---|
592 | # load the fix background info into the histograms |
---|
593 | for hist in Histograms: |
---|
594 | if 'Background' not in Histograms[hist]: continue |
---|
595 | fixedBkg = Histograms[hist]['Background'][1].get('background PWDR') |
---|
596 | if fixedBkg: |
---|
597 | if not fixedBkg[0]: continue |
---|
598 | # patch: add refinement flag, if needed |
---|
599 | if len(fixedBkg) == 2: fixedBkg += [False] |
---|
600 | h = Histograms[hist]['Background'][1] |
---|
601 | try: |
---|
602 | Limits = Histograms[hist]['Limits'][1] |
---|
603 | x = Histograms[hist]['Data'][0] |
---|
604 | xB = np.searchsorted(x,Limits[0]) |
---|
605 | xF = np.searchsorted(x,Limits[1])+1 |
---|
606 | h['fixback'] = allHistograms[fixedBkg[0]]['Data'][1][xB:xF] |
---|
607 | except KeyError: # would happen if a referenced histogram were renamed or deleted |
---|
608 | G2fil.G2Print('Warning: For hist "{}" unresolved background reference to hist "{}"' |
---|
609 | .format(hist,fixedBkg[0])) |
---|
610 | G2obj.IndexAllIds(Histograms=Histograms,Phases=Phases) |
---|
611 | return Histograms,Phases |
---|
612 | |
---|
613 | def getBackupName(GPXfile,makeBack): |
---|
614 | ''' |
---|
615 | Get the name for the backup .gpx file name |
---|
616 | |
---|
617 | :param str GPXfile: full .gpx file name |
---|
618 | :param bool makeBack: if True the name of a new file is returned, if |
---|
619 | False the name of the last file that exists is returned |
---|
620 | :returns: the name of a backup file |
---|
621 | |
---|
622 | ''' |
---|
623 | GPXpath,GPXname = ospath.split(GPXfile) |
---|
624 | if GPXpath == '': GPXpath = '.' |
---|
625 | Name = ospath.splitext(GPXname)[0] |
---|
626 | files = os.listdir(GPXpath) |
---|
627 | last = 0 |
---|
628 | for name in files: |
---|
629 | name = name.split('.') |
---|
630 | if len(name) == 3 and name[0] == Name and 'bak' in name[1]: |
---|
631 | if makeBack: |
---|
632 | last = max(last,int(name[1].strip('bak'))+1) |
---|
633 | else: |
---|
634 | last = max(last,int(name[1].strip('bak'))) |
---|
635 | GPXback = ospath.join(GPXpath,ospath.splitext(GPXname)[0]+'.bak'+str(last)+'.gpx') |
---|
636 | return GPXback |
---|
637 | |
---|
638 | def GPXBackup(GPXfile,makeBack=True): |
---|
639 | ''' |
---|
640 | makes a backup of the specified .gpx file |
---|
641 | |
---|
642 | :param str GPXfile: full .gpx file name |
---|
643 | :param bool makeBack: if True (default), the backup is written to |
---|
644 | a new file; if False, the last backup is overwritten |
---|
645 | :returns: the name of the backup file that was written |
---|
646 | ''' |
---|
647 | import distutils.file_util as dfu |
---|
648 | GPXback = getBackupName(GPXfile,makeBack) |
---|
649 | tries = 0 |
---|
650 | while True: |
---|
651 | try: |
---|
652 | dfu.copy_file(GPXfile,GPXback) |
---|
653 | break |
---|
654 | except: |
---|
655 | tries += 1 |
---|
656 | if tries > 10: |
---|
657 | return GPXfile #failed! |
---|
658 | time.sleep(1) #just wait a second! |
---|
659 | return GPXback |
---|
660 | |
---|
661 | def SetUsedHistogramsAndPhases(GPXfile,Histograms,Phases,RigidBodies,CovData,parmFrozenList,makeBack=True): |
---|
662 | ''' Updates gpxfile from all histograms that are found in any phase |
---|
663 | and any phase that used a histogram. Also updates rigid body definitions. |
---|
664 | This is used for non-sequential fits, but not for sequential fitting. |
---|
665 | |
---|
666 | :param str GPXfile: full .gpx file name |
---|
667 | :param dict Histograms: dictionary of histograms as {name:data,...} |
---|
668 | :param dict Phases: dictionary of phases that use histograms |
---|
669 | :param dict RigidBodies: dictionary of rigid bodies |
---|
670 | :param dict CovData: dictionary of refined variables, varyList, & covariance matrix |
---|
671 | :param list parmFrozenList: list of parameters (as str) that are frozen |
---|
672 | due to limits; converted to :class:`GSASIIobj.G2VarObj` objects. |
---|
673 | :param bool makeBack: True if new backup of .gpx file is to be made; else |
---|
674 | use the last one made |
---|
675 | ''' |
---|
676 | |
---|
677 | import distutils.file_util as dfu |
---|
678 | GPXback = GPXBackup(GPXfile,makeBack) |
---|
679 | G2fil.G2Print ('Read from file:'+GPXback) |
---|
680 | G2fil.G2Print ('Save to file :'+GPXfile) |
---|
681 | infile = open(GPXback,'rb') |
---|
682 | outfile = open(GPXfile,'wb') |
---|
683 | while True: |
---|
684 | try: |
---|
685 | data = cPickleLoad(infile) |
---|
686 | except EOFError: |
---|
687 | break |
---|
688 | datum = data[0] |
---|
689 | # print 'read: ',datum[0] |
---|
690 | if datum[0] == 'Phases': |
---|
691 | for iphase in range(len(data)): |
---|
692 | if data[iphase][0] in Phases: |
---|
693 | phaseName = data[iphase][0] |
---|
694 | data[iphase][1].update(Phases[phaseName]) |
---|
695 | elif datum[0] == 'Covariance': |
---|
696 | data[0][1] = CovData |
---|
697 | elif datum[0] == 'Rigid bodies': |
---|
698 | data[0][1] = RigidBodies |
---|
699 | elif datum[0] == 'Controls': |
---|
700 | Controls = data[0][1] |
---|
701 | if 'parmFrozen' not in Controls: |
---|
702 | Controls['parmFrozen'] = {} |
---|
703 | Controls['parmFrozen']['FrozenList'] = [i if type(i) is G2obj.G2VarObj |
---|
704 | else G2obj.G2VarObj(i) for i in parmFrozenList] |
---|
705 | try: |
---|
706 | histogram = Histograms[datum[0]] |
---|
707 | # print 'found ',datum[0] |
---|
708 | data[0][1][1] = histogram['Data'] |
---|
709 | data[0][1][0].update(histogram['Residuals']) |
---|
710 | for datus in data[1:]: |
---|
711 | # print ' read: ',datus[0] |
---|
712 | if datus[0] in ['Instrument Parameters','Sample Parameters','Reflection Lists']: |
---|
713 | datus[1] = histogram[datus[0]] |
---|
714 | if datus[0] == 'Background': # remove fixed background from file |
---|
715 | d1 = {key:histogram['Background'][1][key] |
---|
716 | for key in histogram['Background'][1] |
---|
717 | if not key.startswith('_fixed')} |
---|
718 | datus[1] = copy.deepcopy(histogram['Background']) |
---|
719 | datus[1][1] = d1 |
---|
720 | except KeyError: |
---|
721 | pass |
---|
722 | try: |
---|
723 | cPickle.dump(data,outfile,1) |
---|
724 | except AttributeError: |
---|
725 | G2fil.G2Print ('ERROR - bad data in least squares result') |
---|
726 | infile.close() |
---|
727 | outfile.close() |
---|
728 | dfu.copy_file(GPXback,GPXfile) |
---|
729 | G2fil.G2Print ('GPX file save failed - old version retained',mode='error') |
---|
730 | return |
---|
731 | |
---|
732 | infile.close() |
---|
733 | outfile.close() |
---|
734 | |
---|
735 | G2fil.G2Print ('GPX file save successful') |
---|
736 | |
---|
737 | def GetSeqResult(GPXfile): |
---|
738 | ''' |
---|
739 | Returns the sequential results table information from a GPX file. |
---|
740 | Called at the beginning of :meth:`GSASIIstrMain.SeqRefine` |
---|
741 | |
---|
742 | :param str GPXfile: full .gpx file name |
---|
743 | :returns: a dict containing the sequential results table |
---|
744 | ''' |
---|
745 | IndexGPX(GPXfile) |
---|
746 | pos = gpxIndex.get('Sequential results') |
---|
747 | if pos is None: |
---|
748 | return {} |
---|
749 | fl = open(GPXfile,'rb') |
---|
750 | fl.seek(pos) |
---|
751 | datum = cPickleLoad(fl)[0] |
---|
752 | fl.close() |
---|
753 | return datum[1] |
---|
754 | |
---|
755 | def SetupSeqSavePhases(GPXfile): |
---|
756 | '''Initialize the files used to save intermediate results from |
---|
757 | sequential fits. |
---|
758 | ''' |
---|
759 | IndexGPX(GPXfile) |
---|
760 | # load initial Phase results from GPX |
---|
761 | fl = open(GPXfile,'rb') |
---|
762 | pos = gpxIndex.get('Phases') |
---|
763 | if pos is None: |
---|
764 | raise Exception("No Phases in GPX file") |
---|
765 | fl.seek(pos) |
---|
766 | data = cPickleLoad(fl) |
---|
767 | fl.close() |
---|
768 | # create GPX-like file to store latest Phase info; init with start vals |
---|
769 | GPXphase = os.path.splitext(GPXfile)[0]+'.seqPhase' |
---|
770 | fp = open(GPXphase,'wb') |
---|
771 | cPickle.dump(data,fp,1) |
---|
772 | fp.close() |
---|
773 | # create empty file for histogram info |
---|
774 | GPXhist = os.path.splitext(GPXfile)[0]+'.seqHist' |
---|
775 | fp = open(GPXhist,'wb') |
---|
776 | fp.close() |
---|
777 | |
---|
778 | def SaveUpdatedHistogramsAndPhases(GPXfile,Histograms,Phases,RigidBodies,CovData,parmFrozen): |
---|
779 | ''' |
---|
780 | Save phase and histogram information into "pseudo-gpx" files. The phase |
---|
781 | information is overwritten each time this is called, but histogram information is |
---|
782 | appended after each sequential step. |
---|
783 | |
---|
784 | :param str GPXfile: full .gpx file name |
---|
785 | :param dict Histograms: dictionary of histograms as {name:data,...} |
---|
786 | :param dict Phases: dictionary of phases that use histograms |
---|
787 | :param dict RigidBodies: dictionary of rigid bodies |
---|
788 | :param dict CovData: dictionary of refined variables, varyList, & covariance matrix |
---|
789 | :param dict parmFrozen: dict with frozen parameters for all phases |
---|
790 | and histograms (specified as str values) |
---|
791 | ''' |
---|
792 | |
---|
793 | GPXphase = os.path.splitext(GPXfile)[0]+'.seqPhase' |
---|
794 | fp = open(GPXphase,'rb') |
---|
795 | data = cPickleLoad(fp) # first block in file should be Phases |
---|
796 | if data[0][0] != 'Phases': |
---|
797 | raise Exception('Unexpected block in {} file. How did this happen?' |
---|
798 | .format(GPXphase)) |
---|
799 | fp.close() |
---|
800 | # update previous phase info |
---|
801 | for datum in data[1:]: |
---|
802 | if datum[0] in Phases: |
---|
803 | datum[1].update(Phases[datum[0]]) |
---|
804 | # save latest Phase/refinement info |
---|
805 | fp = open(GPXphase,'wb') |
---|
806 | cPickle.dump(data,fp,1) |
---|
807 | cPickle.dump([['Covariance',CovData]],fp,1) |
---|
808 | cPickle.dump([['Rigid bodies',RigidBodies]],fp,1) |
---|
809 | cPickle.dump([['parmFrozen',parmFrozen]],fp,1) |
---|
810 | fp.close() |
---|
811 | # create an entry that looks like a PWDR tree item |
---|
812 | for key in Histograms: |
---|
813 | if key.startswith('PWDR '): |
---|
814 | break |
---|
815 | else: |
---|
816 | raise Exception('No PWDR entry in Histogram dict!') |
---|
817 | histname = key |
---|
818 | hist = copy.deepcopy(Histograms[key]) |
---|
819 | xfer_dict = {'Index Peak List': [[], []], |
---|
820 | 'Comments': [], |
---|
821 | 'Unit Cells List': [], |
---|
822 | 'Peak List': {'peaks': [], 'sigDict': {}}, |
---|
823 | } |
---|
824 | histData = hist['Data'] |
---|
825 | del hist['Data'] |
---|
826 | for key in ('Limits','Background','Instrument Parameters', |
---|
827 | 'Sample Parameters','Reflection Lists'): |
---|
828 | xfer_dict[key] = hist[key] |
---|
829 | if key == 'Background': # remove fixed background from file |
---|
830 | xfer_dict['Background'][1] = {k:hist['Background'][1][k] |
---|
831 | for k in hist['Background'][1] |
---|
832 | if not k.startswith('_fixed')} |
---|
833 | del hist[key] |
---|
834 | # xform into a gpx-type entry |
---|
835 | data = [] |
---|
836 | data.append([histname,[hist,histData,histname]]) |
---|
837 | for key in ['Comments','Limits','Background','Instrument Parameters', |
---|
838 | 'Sample Parameters','Peak List','Index Peak List', |
---|
839 | 'Unit Cells List','Reflection Lists']: |
---|
840 | data.append([key,xfer_dict[key]]) |
---|
841 | # append histogram to histogram info |
---|
842 | GPXhist = os.path.splitext(GPXfile)[0]+'.seqHist' |
---|
843 | fp = open(GPXhist,'ab') |
---|
844 | cPickle.dump(data,fp,1) |
---|
845 | fp.close() |
---|
846 | return |
---|
847 | |
---|
848 | def SetSeqResult(GPXfile,Histograms,SeqResult): |
---|
849 | ''' |
---|
850 | Places the sequential results information into a GPX file |
---|
851 | after a refinement has been completed. |
---|
852 | Called at the end of :meth:`GSASIIstrMain.SeqRefine` |
---|
853 | |
---|
854 | :param str GPXfile: full .gpx file name |
---|
855 | ''' |
---|
856 | GPXback = GPXBackup(GPXfile) |
---|
857 | G2fil.G2Print ('Read from file:'+GPXback) |
---|
858 | G2fil.G2Print ('Save to file :'+GPXfile) |
---|
859 | GPXphase = os.path.splitext(GPXfile)[0]+'.seqPhase' |
---|
860 | fp = open(GPXphase,'rb') |
---|
861 | data = cPickleLoad(fp) # first block in file should be Phases |
---|
862 | if data[0][0] != 'Phases': |
---|
863 | raise Exception('Unexpected block in {} file. How did this happen?'.format(GPXphase)) |
---|
864 | Phases = {} |
---|
865 | for name,vals in data[1:]: |
---|
866 | Phases[name] = vals |
---|
867 | name,CovData = cPickleLoad(fp)[0] # 2nd block in file should be Covariance |
---|
868 | name,RigidBodies = cPickleLoad(fp)[0] # 3rd block in file should be Rigid Bodies |
---|
869 | name,parmFrozenDict = cPickleLoad(fp)[0] # 4th block in file should be frozen parameters |
---|
870 | fp.close() |
---|
871 | GPXhist = os.path.splitext(GPXfile)[0]+'.seqHist' |
---|
872 | hist = open(GPXhist,'rb') |
---|
873 | # build an index to the GPXhist file |
---|
874 | histIndex = {} |
---|
875 | while True: |
---|
876 | loc = hist.tell() |
---|
877 | try: |
---|
878 | datum = cPickleLoad(hist)[0] |
---|
879 | except EOFError: |
---|
880 | break |
---|
881 | histIndex[datum[0]] = loc |
---|
882 | |
---|
883 | infile = open(GPXback,'rb') |
---|
884 | outfile = open(GPXfile,'wb') |
---|
885 | while True: |
---|
886 | try: |
---|
887 | data = cPickleLoad(infile) |
---|
888 | except EOFError: |
---|
889 | break |
---|
890 | datum = data[0] |
---|
891 | if datum[0] == 'Sequential results': |
---|
892 | data[0][1] = SeqResult |
---|
893 | elif datum[0] == 'Phases': |
---|
894 | for pdata in data[1:]: |
---|
895 | if pdata[0] in Phases: |
---|
896 | pdata[1].update(Phases[pdata[0]]) |
---|
897 | elif datum[0] == 'Covariance': |
---|
898 | data[0][1] = CovData |
---|
899 | elif datum[0] == 'Rigid bodies': |
---|
900 | data[0][1] = RigidBodies |
---|
901 | elif datum[0] == 'Controls': # reset the Copy Next flag after a sequential fit |
---|
902 | Controls = data[0][1] |
---|
903 | Controls['Copy2Next'] = False |
---|
904 | for key in parmFrozenDict: |
---|
905 | Controls['parmFrozen'][key] = [ |
---|
906 | i if type(i) is G2obj.G2VarObj |
---|
907 | else G2obj.G2VarObj(i) |
---|
908 | for i in parmFrozenDict[key]] |
---|
909 | elif datum[0] in histIndex: |
---|
910 | hist.seek(histIndex[datum[0]]) |
---|
911 | hdata = cPickleLoad(hist) |
---|
912 | if data[0][0] != hdata[0][0]: |
---|
913 | G2fil.G2Print('Error! Updating {} with {}'.format(data[0][0],hdata[0][0])) |
---|
914 | data[0] = hdata[0] |
---|
915 | xferItems = ['Background','Instrument Parameters','Sample Parameters','Reflection Lists'] |
---|
916 | hItems = {name:j+1 for j,(name,val) in enumerate(hdata[1:]) if name in xferItems} |
---|
917 | for j,(name,val) in enumerate(data[1:]): |
---|
918 | if name not in xferItems: continue |
---|
919 | data[j+1][1] = hdata[hItems[name]][1] |
---|
920 | cPickle.dump(data,outfile,1) |
---|
921 | hist.close() |
---|
922 | infile.close() |
---|
923 | outfile.close() |
---|
924 | # clean up tmp files |
---|
925 | try: |
---|
926 | os.remove(GPXphase) |
---|
927 | except: |
---|
928 | G2fil.G2Print('Warning: unable to delete {}'.format(GPXphase)) |
---|
929 | try: |
---|
930 | os.remove(GPXhist) |
---|
931 | except: |
---|
932 | G2fil.G2Print('Warning: unable to delete {}'.format(GPXhist)) |
---|
933 | G2fil.G2Print ('GPX file merge completed') |
---|
934 | |
---|
935 | #============================================================================== |
---|
936 | # Refinement routines |
---|
937 | #============================================================================== |
---|
938 | def ShowBanner(pFile=None): |
---|
939 | 'Print authorship, copyright and citation notice' |
---|
940 | pFile.write(80*'*'+'\n') |
---|
941 | pFile.write(' General Structure Analysis System-II Crystal Structure Refinement\n') |
---|
942 | pFile.write(' by Robert B. Von Dreele & Brian H. Toby\n') |
---|
943 | pFile.write(' Argonne National Laboratory(C), 2010\n') |
---|
944 | pFile.write(' This product includes software developed by the UChicago Argonne, LLC,\n') |
---|
945 | pFile.write(' as Operator of Argonne National Laboratory.\n') |
---|
946 | pFile.write(' Please cite:\n') |
---|
947 | pFile.write(' B.H. Toby & R.B. Von Dreele, J. Appl. Cryst. 46, 544-549 (2013)\n') |
---|
948 | |
---|
949 | pFile.write(80*'*'+'\n') |
---|
950 | |
---|
951 | def ShowControls(Controls,pFile=None,SeqRef=False,preFrozenCount=0): |
---|
952 | 'Print controls information' |
---|
953 | pFile.write(' Least squares controls:\n') |
---|
954 | pFile.write(' Refinement type: %s\n'%Controls['deriv type']) |
---|
955 | if 'Hessian' in Controls['deriv type']: |
---|
956 | pFile.write(' Maximum number of cycles: %d\n'%Controls['max cyc']) |
---|
957 | else: |
---|
958 | pFile.write(' Minimum delta-M/M for convergence: %.2g\n'%(Controls['min dM/M'])) |
---|
959 | pFile.write(' Regularize hydrogens (if any): %s\n'%Controls.get('HatomFix',False)) |
---|
960 | pFile.write(' Initial shift factor: %.3f\n'%(Controls['shift factor'])) |
---|
961 | if SeqRef: |
---|
962 | pFile.write(' Sequential refinement controls:\n') |
---|
963 | pFile.write(' Copy of histogram results to next: %s\n'%(Controls['Copy2Next'])) |
---|
964 | pFile.write(' Process histograms in reverse order: %s\n'%(Controls['Reverse Seq'])) |
---|
965 | if preFrozenCount: |
---|
966 | pFile.write('\n Starting refinement with {} Frozen variables\n\n'.format(preFrozenCount)) |
---|
967 | |
---|
968 | def GetPawleyConstr(SGLaue,PawleyRef,im,pawleyVary): |
---|
969 | 'needs a doc string' |
---|
970 | # if SGLaue in ['-1','2/m','mmm']: |
---|
971 | # return #no Pawley symmetry required constraints |
---|
972 | eqvDict = {} |
---|
973 | for i,varyI in enumerate(pawleyVary): |
---|
974 | eqvDict[varyI] = [] |
---|
975 | refI = int(varyI.split(':')[-1]) |
---|
976 | ih,ik,il = PawleyRef[refI][:3] |
---|
977 | dspI = PawleyRef[refI][4+im] |
---|
978 | for varyJ in pawleyVary[i+1:]: |
---|
979 | refJ = int(varyJ.split(':')[-1]) |
---|
980 | jh,jk,jl = PawleyRef[refJ][:3] |
---|
981 | dspJ = PawleyRef[refJ][4+im] |
---|
982 | if SGLaue in ['4/m','4/mmm']: |
---|
983 | isum = ih**2+ik**2 |
---|
984 | jsum = jh**2+jk**2 |
---|
985 | if abs(il) == abs(jl) and isum == jsum: |
---|
986 | eqvDict[varyI].append(varyJ) |
---|
987 | elif SGLaue in ['3R','3mR']: |
---|
988 | isum = ih**2+ik**2+il**2 |
---|
989 | jsum = jh**2+jk**2+jl**2 |
---|
990 | isum2 = ih*ik+ih*il+ik*il |
---|
991 | jsum2 = jh*jk+jh*jl+jk*jl |
---|
992 | if isum == jsum and isum2 == jsum2: |
---|
993 | eqvDict[varyI].append(varyJ) |
---|
994 | elif SGLaue in ['3','3m1','31m','6/m','6/mmm']: |
---|
995 | isum = ih**2+ik**2+ih*ik |
---|
996 | jsum = jh**2+jk**2+jh*jk |
---|
997 | if abs(il) == abs(jl) and isum == jsum: |
---|
998 | eqvDict[varyI].append(varyJ) |
---|
999 | elif SGLaue in ['m3','m3m']: |
---|
1000 | isum = ih**2+ik**2+il**2 |
---|
1001 | jsum = jh**2+jk**2+jl**2 |
---|
1002 | if isum == jsum: |
---|
1003 | eqvDict[varyI].append(varyJ) |
---|
1004 | elif abs(dspI-dspJ)/dspI < 1.e-4: |
---|
1005 | eqvDict[varyI].append(varyJ) |
---|
1006 | for item in pawleyVary: |
---|
1007 | if eqvDict[item]: |
---|
1008 | for item2 in pawleyVary: |
---|
1009 | if item2 in eqvDict[item]: |
---|
1010 | eqvDict[item2] = [] |
---|
1011 | G2mv.StoreEquivalence(item,eqvDict[item]) |
---|
1012 | |
---|
1013 | def cellVary(pfx,SGData): |
---|
1014 | '''Creates equivalences for a phase based on the Laue class. |
---|
1015 | Returns a list of A tensor terms that are non-zero. |
---|
1016 | ''' |
---|
1017 | if SGData['SGLaue'] in ['-1',]: |
---|
1018 | return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3',pfx+'A4',pfx+'A5'] |
---|
1019 | elif SGData['SGLaue'] in ['2/m',]: |
---|
1020 | if SGData['SGUniq'] == 'a': |
---|
1021 | return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A5'] |
---|
1022 | elif SGData['SGUniq'] == 'b': |
---|
1023 | return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A4'] |
---|
1024 | else: |
---|
1025 | return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3'] |
---|
1026 | elif SGData['SGLaue'] in ['mmm',]: |
---|
1027 | return [pfx+'A0',pfx+'A1',pfx+'A2'] |
---|
1028 | elif SGData['SGLaue'] in ['4/m','4/mmm']: |
---|
1029 | G2mv.StoreEquivalence(pfx+'A0',(pfx+'A1',)) |
---|
1030 | return [pfx+'A0',pfx+'A1',pfx+'A2'] |
---|
1031 | elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']: |
---|
1032 | G2mv.StoreEquivalence(pfx+'A0',(pfx+'A1',pfx+'A3',)) |
---|
1033 | return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3'] |
---|
1034 | elif SGData['SGLaue'] in ['3R', '3mR']: |
---|
1035 | G2mv.StoreEquivalence(pfx+'A0',(pfx+'A1',pfx+'A2',)) |
---|
1036 | G2mv.StoreEquivalence(pfx+'A3',(pfx+'A4',pfx+'A5',)) |
---|
1037 | return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3',pfx+'A4',pfx+'A5'] |
---|
1038 | elif SGData['SGLaue'] in ['m3m','m3']: |
---|
1039 | G2mv.StoreEquivalence(pfx+'A0',(pfx+'A1',pfx+'A2',)) |
---|
1040 | return [pfx+'A0',pfx+'A1',pfx+'A2'] |
---|
1041 | |
---|
1042 | def modVary(pfx,SSGData): |
---|
1043 | vary = [] |
---|
1044 | for i,item in enumerate(SSGData['modSymb']): |
---|
1045 | if item in ['a','b','g']: |
---|
1046 | vary.append(pfx+'mV%d'%(i)) |
---|
1047 | return vary |
---|
1048 | |
---|
1049 | ################################################################################ |
---|
1050 | ##### Rigid Body Models and not General.get('doPawley') |
---|
1051 | ################################################################################ |
---|
1052 | |
---|
1053 | def GetRigidBodyModels(rigidbodyDict,Print=True,pFile=None): |
---|
1054 | '''Get Rigid body info from tree entry and print it to .LST file |
---|
1055 | Adds variables and dict items for vector RBs, but for Residue bodies |
---|
1056 | this is done in :func:`GetPhaseData`. |
---|
1057 | ''' |
---|
1058 | def PrintSpnRBModel(RBModel): |
---|
1059 | pFile.write('Spinning RB name: %s atom type: %s, no.atoms: %d, No. times used: %d\n'% |
---|
1060 | (RBModel['RBname'],RBModel['atType'],RBModel['Natoms'],RBModel['useCount'])) |
---|
1061 | for i in WriteSpnRBModel(RBModel): |
---|
1062 | pFile.write(i) |
---|
1063 | |
---|
1064 | def PrintResRBModel(RBModel): |
---|
1065 | pFile.write('Residue RB name: %s no.atoms: %d, No. times used: %d\n'% |
---|
1066 | (RBModel['RBname'],len(RBModel['rbTypes']),RBModel['useCount'])) |
---|
1067 | for i in WriteResRBModel(RBModel): |
---|
1068 | pFile.write(i) |
---|
1069 | |
---|
1070 | def PrintVecRBModel(RBModel): |
---|
1071 | pFile.write('Vector RB name: %s no.atoms: %d, No. times used: %d\n'% |
---|
1072 | (RBModel['RBname'],len(RBModel['rbTypes']),RBModel['useCount'])) |
---|
1073 | for i in WriteVecRBModel(RBModel): |
---|
1074 | pFile.write(i) |
---|
1075 | pFile.write('Orientation defined by: atom %s -> atom %s & atom %s -> atom %s\n'% |
---|
1076 | (RBModel['rbRef'][0],RBModel['rbRef'][1],RBModel['rbRef'][0],RBModel['rbRef'][2])) |
---|
1077 | |
---|
1078 | if Print and pFile is None: raise Exception("specify pFile or Print=False") |
---|
1079 | rbVary = [] |
---|
1080 | rbDict = {} |
---|
1081 | rbIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[],'Spin':[]}) |
---|
1082 | if len(rbIds.get('Spin',{})): |
---|
1083 | for irb,item in enumerate(rbIds['Spin']): |
---|
1084 | if rigidbodyDict['Spin'][item]['useCount']: |
---|
1085 | RBradius = rigidbodyDict['Spin'][item]['radius'] |
---|
1086 | pid = '::RBS;0:'+str(irb) |
---|
1087 | rbDict[pid] = RBradius[0] |
---|
1088 | if RBradius[1]: |
---|
1089 | rbVary.append(pid) |
---|
1090 | if Print: |
---|
1091 | pFile.write('\nSpinning rigid body model:\n') |
---|
1092 | PrintSpnRBModel(rigidbodyDict['Spin'][item]) |
---|
1093 | |
---|
1094 | if len(rbIds['Vector']): |
---|
1095 | for irb,item in enumerate(rbIds['Vector']): |
---|
1096 | if rigidbodyDict['Vector'][item]['useCount']: |
---|
1097 | RBmags = rigidbodyDict['Vector'][item]['VectMag'] |
---|
1098 | RBrefs = rigidbodyDict['Vector'][item]['VectRef'] |
---|
1099 | for i,[mag,ref] in enumerate(zip(RBmags,RBrefs)): |
---|
1100 | pid = '::RBV;'+str(i)+':'+str(irb) |
---|
1101 | rbDict[pid] = mag |
---|
1102 | if ref: |
---|
1103 | rbVary.append(pid) |
---|
1104 | if Print: |
---|
1105 | pFile.write('\nVector rigid body model:\n') |
---|
1106 | PrintVecRBModel(rigidbodyDict['Vector'][item]) |
---|
1107 | if Print: |
---|
1108 | if len(rbIds['Residue']): |
---|
1109 | for item in rbIds['Residue']: |
---|
1110 | if rigidbodyDict['Residue'][item]['useCount']: |
---|
1111 | pFile.write('\nResidue rigid body model:\n') |
---|
1112 | PrintResRBModel(rigidbodyDict['Residue'][item]) |
---|
1113 | return rbVary,rbDict |
---|
1114 | |
---|
1115 | def SetRigidBodyModels(parmDict,sigDict,rigidbodyDict,pFile=None): |
---|
1116 | 'needs a doc string' |
---|
1117 | |
---|
1118 | def PrintRBVectandSig(VectRB,VectSig): |
---|
1119 | pFile.write('\n Rigid body vector magnitudes for %s:\n'%VectRB['RBname']) |
---|
1120 | namstr = ' names :' |
---|
1121 | valstr = ' values:' |
---|
1122 | sigstr = ' esds :' |
---|
1123 | for i,[val,sig] in enumerate(zip(VectRB['VectMag'],VectSig)): |
---|
1124 | namstr += '%12s'%('Vect '+str(i)) |
---|
1125 | valstr += '%12.4f'%(val) |
---|
1126 | if sig: |
---|
1127 | sigstr += '%12.4f'%(sig) |
---|
1128 | else: |
---|
1129 | sigstr += 12*' ' |
---|
1130 | pFile.write(namstr+'\n') |
---|
1131 | pFile.write(valstr+'\n') |
---|
1132 | pFile.write(sigstr+'\n') |
---|
1133 | |
---|
1134 | def PrintRBSpnandSig(SpinRB,SpinSig): |
---|
1135 | pFile.write('\n Spinning radius for %s:\n'%SpinRB['RBname']) |
---|
1136 | namstr = ' names :' |
---|
1137 | valstr = ' values:' |
---|
1138 | sigstr = ' esds :' |
---|
1139 | for i,[val,sig] in enumerate(zip(SpinRB['radius'],SpinSig)): |
---|
1140 | namstr += '%12s'%('Spin '+str(i)) |
---|
1141 | valstr += '%12.4f'%(val) |
---|
1142 | if sig: |
---|
1143 | sigstr += '%12.4f'%(sig) |
---|
1144 | else: |
---|
1145 | sigstr += 12*' ' |
---|
1146 | pFile.write(namstr+'\n') |
---|
1147 | pFile.write(valstr+'\n') |
---|
1148 | pFile.write(sigstr+'\n') |
---|
1149 | |
---|
1150 | RBIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[],'Spin':[]}) #these are lists of rbIds |
---|
1151 | if not RBIds['Vector'] and not RBIds['Spin']: |
---|
1152 | return |
---|
1153 | for irb,item in enumerate(RBIds['Vector']): |
---|
1154 | if rigidbodyDict['Vector'][item]['useCount']: |
---|
1155 | VectSig = [] |
---|
1156 | RBmags = rigidbodyDict['Vector'][item]['VectMag'] |
---|
1157 | for i,mag in enumerate(RBmags): |
---|
1158 | name = '::RBV;'+str(i)+':'+str(irb) |
---|
1159 | if name in sigDict: |
---|
1160 | VectSig.append(sigDict[name]) |
---|
1161 | PrintRBVectandSig(rigidbodyDict['Vector'][item],VectSig) |
---|
1162 | for irb,item in enumerate(RBIds['Spin']): |
---|
1163 | if rigidbodyDict['Spin'][item]['useCount']: |
---|
1164 | SpinSig = [] |
---|
1165 | name = '::RBS;'+str(0)+':'+str(irb) |
---|
1166 | if name in sigDict: |
---|
1167 | SpinSig.append(sigDict[name]) |
---|
1168 | |
---|
1169 | |
---|
1170 | PrintRBSpnandSig(rigidbodyDict['Spin'][item],SpinSig) |
---|
1171 | |
---|
1172 | ################################################################################ |
---|
1173 | ##### Phase data |
---|
1174 | ################################################################################ |
---|
1175 | def GetPhaseData(PhaseData,RestraintDict={},rbIds={},Print=True,pFile=None, |
---|
1176 | seqHistName=None,symHold=None): |
---|
1177 | '''Setup the phase information for a structural refinement, used for |
---|
1178 | regular and sequential refinements, optionally printing information |
---|
1179 | to the .lst file (if Print is True) |
---|
1180 | |
---|
1181 | :param dict PhaseData: the contents of the Phase tree item (may be read from |
---|
1182 | .gpx file) with information on all phases |
---|
1183 | :param dict RestraintDict: an optional dict with restraint information |
---|
1184 | :param dict rbIds: an optional dict with rigid body information |
---|
1185 | :param bool Print: a flag that determines if information will be formatted and |
---|
1186 | printed to the .lst file |
---|
1187 | :param file pFile: a file object (created by open) where print information is sent |
---|
1188 | when Print is True |
---|
1189 | :param str seqHistName: will be None, except for sequential fits. For sequential |
---|
1190 | fits, this can be the name of the current histogram or 'All'. If a histogram |
---|
1191 | name is supplied, only the phases used in the current histogram are loaded. |
---|
1192 | If 'All' is specified, all phases are loaded (used for error checking). |
---|
1193 | :param list symHold: if not None (None is the default) the names of parameters |
---|
1194 | held due to symmetry are placed in this list |
---|
1195 | :returns: lots of stuff: Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup, |
---|
1196 | FFtables,EFtables,BLtables,MFtables,maxSSwave (see code for details). |
---|
1197 | ''' |
---|
1198 | |
---|
1199 | def PrintFFtable(FFtable): |
---|
1200 | pFile.write('\n X-ray scattering factors:\n') |
---|
1201 | pFile.write(' Symbol fa fb fc\n') |
---|
1202 | pFile.write(99*'-'+'\n') |
---|
1203 | for Ename in FFtable: |
---|
1204 | ffdata = FFtable[Ename] |
---|
1205 | fa = ffdata['fa'] |
---|
1206 | fb = ffdata['fb'] |
---|
1207 | pFile.write(' %8s %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f\n'% |
---|
1208 | (Ename.ljust(8),fa[0],fa[1],fa[2],fa[3],fb[0],fb[1],fb[2],fb[3],ffdata['fc'])) |
---|
1209 | |
---|
1210 | def PrintEFtable(EFtable): |
---|
1211 | pFile.write('\n Electron scattering factors:\n') |
---|
1212 | pFile.write(' Symbol fa fb\n') |
---|
1213 | pFile.write(99*'-'+'\n') |
---|
1214 | for Ename in EFtable: |
---|
1215 | efdata = EFtable[Ename] |
---|
1216 | fa = efdata['fa'] |
---|
1217 | fb = efdata['fb'] |
---|
1218 | pFile.write(' %8s %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f\n'% |
---|
1219 | (Ename.ljust(8),fa[0],fa[1],fa[2],fa[3],fa[4],fb[0],fb[1],fb[2],fb[3],fb[4])) |
---|
1220 | |
---|
1221 | def PrintMFtable(MFtable): |
---|
1222 | pFile.write('\n <j0> Magnetic scattering factors:\n') |
---|
1223 | pFile.write(' Symbol mfa mfb mfc\n') |
---|
1224 | pFile.write(99*'-'+'\n') |
---|
1225 | for Ename in MFtable: |
---|
1226 | mfdata = MFtable[Ename] |
---|
1227 | fa = mfdata['mfa'] |
---|
1228 | fb = mfdata['mfb'] |
---|
1229 | pFile.write(' %8s %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f\n'% |
---|
1230 | (Ename.ljust(8),fa[0],fa[1],fa[2],fa[3],fb[0],fb[1],fb[2],fb[3],mfdata['mfc'])) |
---|
1231 | pFile.write('\n <j2> Magnetic scattering factors:\n') |
---|
1232 | pFile.write(' Symbol nfa nfb nfc\n') |
---|
1233 | pFile.write(99*'-'+'\n') |
---|
1234 | for Ename in MFtable: |
---|
1235 | mfdata = MFtable[Ename] |
---|
1236 | fa = mfdata['nfa'] |
---|
1237 | fb = mfdata['nfb'] |
---|
1238 | pFile.write(' %8s %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f\n'% |
---|
1239 | (Ename.ljust(8),fa[0],fa[1],fa[2],fa[3],fb[0],fb[1],fb[2],fb[3],mfdata['nfc'])) |
---|
1240 | |
---|
1241 | def PrintBLtable(BLtable): |
---|
1242 | pFile.write('\n Neutron scattering factors:\n') |
---|
1243 | pFile.write(' Symbol isotope mass b resonant terms\n') |
---|
1244 | pFile.write(99*'-'+'\n') |
---|
1245 | for Ename in BLtable: |
---|
1246 | bldata = BLtable[Ename] |
---|
1247 | isotope = bldata[0] |
---|
1248 | mass = bldata[1]['Mass'] |
---|
1249 | if 'BW-LS' in bldata[1]: |
---|
1250 | bres = bldata[1]['BW-LS'] |
---|
1251 | blen = 0 |
---|
1252 | else: |
---|
1253 | blen = bldata[1]['SL'][0] |
---|
1254 | bres = [] |
---|
1255 | line = ' %8s%11s %10.3f %8.3f'%(Ename.ljust(8),isotope.center(11),mass,blen) |
---|
1256 | for item in bres: |
---|
1257 | line += '%10.5g'%(item) |
---|
1258 | pFile.write(line+'\n') |
---|
1259 | |
---|
1260 | def PrintRBObjects(resRBData,vecRBData,spnRBData): |
---|
1261 | |
---|
1262 | def PrintRBThermals(): |
---|
1263 | tlstr = ['11','22','33','12','13','23'] |
---|
1264 | sstr = ['12','13','21','23','31','32','AA','BB'] |
---|
1265 | TLS = RB['ThermalMotion'][1] |
---|
1266 | TLSvar = RB['ThermalMotion'][2] |
---|
1267 | if 'T' in RB['ThermalMotion'][0]: |
---|
1268 | pFile.write('TLS data\n') |
---|
1269 | text = '' |
---|
1270 | for i in range(6): |
---|
1271 | text += 'T'+tlstr[i]+' %8.4f %s '%(TLS[i],str(TLSvar[i])[0]) |
---|
1272 | pFile.write(text+'\n') |
---|
1273 | if 'L' in RB['ThermalMotion'][0]: |
---|
1274 | text = '' |
---|
1275 | for i in range(6,12): |
---|
1276 | text += 'L'+tlstr[i-6]+' %8.2f %s '%(TLS[i],str(TLSvar[i])[0]) |
---|
1277 | pFile.write(text+'\n') |
---|
1278 | if 'S' in RB['ThermalMotion'][0]: |
---|
1279 | text = '' |
---|
1280 | for i in range(12,20): |
---|
1281 | text += 'S'+sstr[i-12]+' %8.3f %s '%(TLS[i],str(TLSvar[i])[0]) |
---|
1282 | pFile.write(text+'\n') |
---|
1283 | if 'U' in RB['ThermalMotion'][0]: |
---|
1284 | pFile.write('Uiso data\n') |
---|
1285 | text = 'Uiso'+' %10.3f %s'%(TLS[0],str(TLSvar[0])[0]) |
---|
1286 | pFile.write(text+'\n') |
---|
1287 | |
---|
1288 | if len(resRBData): |
---|
1289 | for RB in resRBData: |
---|
1290 | Oxyz = RB['Orig'][0] |
---|
1291 | Qrijk = RB['Orient'][0] |
---|
1292 | Angle = 2.0*acosd(Qrijk[0]) |
---|
1293 | pFile.write('\nRBObject %s at %10.4f %10.4f %10.4f Refine? %s\n'% |
---|
1294 | (RB['RBname'],Oxyz[0],Oxyz[1],Oxyz[2],RB['Orig'][1])) |
---|
1295 | pFile.write('Orientation angle,vector: %10.3f %10.4f %10.4f %10.4f Refine? %s\n'% |
---|
1296 | (Angle,Qrijk[1],Qrijk[2],Qrijk[3],RB['Orient'][1])) |
---|
1297 | pFile.write('Atom site frac: %10.3f Refine? %s\n'%(RB['AtomFrac'][0],RB['AtomFrac'][1])) |
---|
1298 | Torsions = RB['Torsions'] |
---|
1299 | if len(Torsions): |
---|
1300 | text = 'Torsions: ' |
---|
1301 | for torsion in Torsions: |
---|
1302 | text += '%10.4f Refine? %s'%(torsion[0],torsion[1]) |
---|
1303 | pFile.write(text+'\n') |
---|
1304 | PrintRBThermals() |
---|
1305 | |
---|
1306 | if len(vecRBData): |
---|
1307 | for RB in vecRBData: |
---|
1308 | Oxyz = RB['Orig'][0] |
---|
1309 | Qrijk = RB['Orient'][0] |
---|
1310 | Angle = 2.0*acosd(Qrijk[0]) |
---|
1311 | pFile.write('\nRBObject %s at %10.4f %10.4f %10.4f Refine? %s\n'% |
---|
1312 | (RB['RBname'],Oxyz[0],Oxyz[1],Oxyz[2],RB['Orig'][1])) |
---|
1313 | pFile.write('Orientation angle,vector: %10.3f %10.4f %10.4f %10.4f Refine? %s\n'% |
---|
1314 | (Angle,Qrijk[1],Qrijk[2],Qrijk[3],RB['Orient'][1])) |
---|
1315 | pFile.write('Atom site frac: %10.3f Refine? %s\n'%(RB['AtomFrac'][0],RB['AtomFrac'][1])) |
---|
1316 | PrintRBThermals() |
---|
1317 | |
---|
1318 | if len(spnRBData): #TODO: fix this |
---|
1319 | for RB in spnRBData: |
---|
1320 | Oxyz = RB['Orig'][0] |
---|
1321 | Qrijk = RB['Orient'][0] |
---|
1322 | Angle = 2.0*acosd(Qrijk[0]) |
---|
1323 | pFile.write('\nRBObject %s at %10.4f %10.4f %10.4f Refine? %s\n'% |
---|
1324 | (RB['RBname'][0],Oxyz[0],Oxyz[1],Oxyz[2],RB['Orig'][1])) |
---|
1325 | pFile.write('Orientation angle,vector: %10.3f %10.4f %10.4f %10.4f Refine? %s\n'% |
---|
1326 | (Angle,Qrijk[1],Qrijk[2],Qrijk[3],RB['Orient'][1])) |
---|
1327 | pFile.write('Bessel/Spherical Harmonics coefficients; symmetry required sign shown\n') |
---|
1328 | for ish,SHC in enumerate(RB['SHC']): |
---|
1329 | if len(SHC): |
---|
1330 | pFile.write('Spin shell (%s) shell no: %d\n'%(RB['RBname'][ish],ish)) |
---|
1331 | SHkeys = list(SHC.keys()) |
---|
1332 | nCoeff = len(SHC) |
---|
1333 | nBlock = nCoeff//6+1 |
---|
1334 | iBeg = 0 |
---|
1335 | iFin = min(iBeg+6,nCoeff) |
---|
1336 | for block in range(nBlock): |
---|
1337 | ptlbls = ' names :' |
---|
1338 | ptstr = ' values:' |
---|
1339 | ptref = ' refine:' |
---|
1340 | for item in SHkeys[iBeg:iFin]: |
---|
1341 | ptlbls += '%12s'%(item) |
---|
1342 | ptstr += '%9.4f*%2d'%(SHC[item][0],SHC[item][1]) |
---|
1343 | ptref += '%12s'%(SHC[item][2]) |
---|
1344 | pFile.write(ptlbls+'\n') |
---|
1345 | pFile.write(ptstr+'\n') |
---|
1346 | pFile.write(ptref+'\n') |
---|
1347 | iBeg += 6 |
---|
1348 | iFin = min(iBeg+6,nCoeff) |
---|
1349 | |
---|
1350 | def PrintAtoms(General,Atoms): |
---|
1351 | cx,ct,cs,cia = General['AtomPtrs'] |
---|
1352 | pFile.write('\n Atoms:\n') |
---|
1353 | line = ' name type refine? x y z '+ \ |
---|
1354 | ' frac site sym mult I/A Uiso U11 U22 U33 U12 U13 U23' |
---|
1355 | if General['Type'] == 'macromolecular': |
---|
1356 | line = ' res no residue chain'+line |
---|
1357 | pFile.write(line+'\n') |
---|
1358 | if General['Type'] in ['nuclear','magnetic','faulted',]: |
---|
1359 | pFile.write(135*'-'+'\n') |
---|
1360 | for i,at in enumerate(Atoms): |
---|
1361 | line = '%7s'%(at[ct-1])+'%7s'%(at[ct])+'%7s'%(at[ct+1])+'%10.5f'%(at[cx])+'%10.5f'%(at[cx+1])+ \ |
---|
1362 | '%10.5f'%(at[cx+2])+'%8.3f'%(at[cx+3])+'%7s'%(at[cs])+'%5d'%(at[cs+1])+'%5s'%(at[cia]) |
---|
1363 | if at[cia] == 'I': |
---|
1364 | line += '%8.5f'%(at[cia+1])+48*' ' |
---|
1365 | else: |
---|
1366 | line += 8*' ' |
---|
1367 | for j in range(6): |
---|
1368 | line += '%8.5f'%(at[cia+2+j]) |
---|
1369 | pFile.write(line+'\n') |
---|
1370 | elif General['Type'] == 'macromolecular': |
---|
1371 | pFile.write(135*'-'+'\n') |
---|
1372 | for i,at in enumerate(Atoms): |
---|
1373 | line = '%7s'%(at[0])+'%7s'%(at[1])+'%7s'%(at[2])+'%7s'%(at[ct-1])+'%7s'%(at[ct])+'%7s'%(at[ct+1])+'%10.5f'%(at[cx])+'%10.5f'%(at[cx+1])+ \ |
---|
1374 | '%10.5f'%(at[cx+2])+'%8.3f'%(at[cx+3])+'%7s'%(at[cs])+'%5d'%(at[cs+1])+'%5s'%(at[cia]) |
---|
1375 | if at[cia] == 'I': |
---|
1376 | line += '%8.4f'%(at[cia+1])+48*' ' |
---|
1377 | else: |
---|
1378 | line += 8*' ' |
---|
1379 | for j in range(6): |
---|
1380 | line += '%8.4f'%(at[cia+2+j]) |
---|
1381 | pFile.write(line+'\n') |
---|
1382 | |
---|
1383 | def PrintMoments(General,Atoms): |
---|
1384 | cx,ct,cs,cia = General['AtomPtrs'] |
---|
1385 | cmx = cx+4 |
---|
1386 | AtInfo = dict(zip(General['AtomTypes'],General['Lande g'])) |
---|
1387 | pFile.write('\n Magnetic moments:\n') |
---|
1388 | line = ' name type refine? Mx My Mz ' |
---|
1389 | pFile.write(line+'\n') |
---|
1390 | pFile.write(135*'-'+'\n') |
---|
1391 | for i,at in enumerate(Atoms): |
---|
1392 | if AtInfo[at[ct]]: |
---|
1393 | line = '%7s'%(at[ct-1])+'%7s'%(at[ct])+'%7s'%(at[ct+1])+'%10.5f'%(at[cmx])+'%10.5f'%(at[cmx+1])+ \ |
---|
1394 | '%10.5f'%(at[cmx+2]) |
---|
1395 | pFile.write(line+'\n') |
---|
1396 | |
---|
1397 | |
---|
1398 | def PrintWaves(General,Atoms): |
---|
1399 | cx,ct,cs,cia = General['AtomPtrs'] |
---|
1400 | pFile.write('\n Modulation waves\n') |
---|
1401 | names = {'Sfrac':['Fsin','Fcos','Fzero','Fwid'],'Spos':['Xsin','Ysin','Zsin','Xcos','Ycos','Zcos','Tmin','Tmax','Xmax','Ymax','Zmax'], |
---|
1402 | 'Sadp':['U11sin','U22sin','U33sin','U12sin','U13sin','U23sin','U11cos','U22cos', |
---|
1403 | 'U33cos','U12cos','U13cos','U23cos'],'Smag':['MXsin','MYsin','MZsin','MXcos','MYcos','MZcos']} |
---|
1404 | pFile.write(135*'-'+'\n') |
---|
1405 | for i,at in enumerate(Atoms): |
---|
1406 | AtomSS = at[-1]['SS1'] |
---|
1407 | for Stype in ['Sfrac','Spos','Sadp','Smag']: |
---|
1408 | Waves = AtomSS[Stype] |
---|
1409 | if len(Waves): |
---|
1410 | pFile.write(' atom: %s, site sym: %s, type: %s wave type: %s:\n'% |
---|
1411 | (at[ct-1],at[cs],Stype,Waves[0])) |
---|
1412 | else: |
---|
1413 | continue |
---|
1414 | for iw,wave in enumerate(Waves[1:]): |
---|
1415 | line = '' |
---|
1416 | if Waves[0] in ['Block','ZigZag'] and Stype == 'Spos' and not iw: |
---|
1417 | for item in names[Stype][6:]: |
---|
1418 | line += '%8s '%(item) |
---|
1419 | else: |
---|
1420 | if Stype == 'Spos': |
---|
1421 | for item in names[Stype][:6]: |
---|
1422 | line += '%8s '%(item) |
---|
1423 | else: |
---|
1424 | for item in names[Stype]: |
---|
1425 | line += '%8s '%(item) |
---|
1426 | pFile.write(line+'\n') |
---|
1427 | line = '' |
---|
1428 | for item in wave[0]: |
---|
1429 | line += '%8.4f '%(item) |
---|
1430 | line += ' Refine? '+str(wave[1]) |
---|
1431 | pFile.write(line+'\n') |
---|
1432 | |
---|
1433 | def PrintTexture(textureData): |
---|
1434 | topstr = '\n Spherical harmonics texture: Order:' + \ |
---|
1435 | str(textureData['Order']) |
---|
1436 | if textureData['Order']: |
---|
1437 | pFile.write('%s Refine? %s\n'%(topstr,textureData['SH Coeff'][0])) |
---|
1438 | else: |
---|
1439 | pFile.write(topstr+'\n') |
---|
1440 | return |
---|
1441 | names = ['omega','chi','phi'] |
---|
1442 | line = '\n' |
---|
1443 | for name in names: |
---|
1444 | line += ' SH '+name+':'+'%12.4f'%(textureData['Sample '+name][1])+' Refine? '+str(textureData['Sample '+name][0]) |
---|
1445 | pFile.write(line+'\n') |
---|
1446 | pFile.write('\n Texture coefficients:\n') |
---|
1447 | SHcoeff = textureData['SH Coeff'][1] |
---|
1448 | SHkeys = list(SHcoeff.keys()) |
---|
1449 | nCoeff = len(SHcoeff) |
---|
1450 | nBlock = nCoeff//10+1 |
---|
1451 | iBeg = 0 |
---|
1452 | iFin = min(iBeg+10,nCoeff) |
---|
1453 | for block in range(nBlock): |
---|
1454 | ptlbls = ' names :' |
---|
1455 | ptstr = ' values:' |
---|
1456 | for item in SHkeys[iBeg:iFin]: |
---|
1457 | ptlbls += '%12s'%(item) |
---|
1458 | ptstr += '%12.4f'%(SHcoeff[item]) |
---|
1459 | pFile.write(ptlbls+'\n') |
---|
1460 | pFile.write(ptstr+'\n') |
---|
1461 | iBeg += 10 |
---|
1462 | iFin = min(iBeg+10,nCoeff) |
---|
1463 | |
---|
1464 | def MakeRBParms(rbKey,phaseVary,phaseDict): |
---|
1465 | #### patch 2/24/21 BHT: new param, AtomFrac in RB |
---|
1466 | if 'AtomFrac' not in RB and rbKey != 'S': raise Exception('out of date RB: edit in RB Models') |
---|
1467 | # end patch |
---|
1468 | if rbKey == 'S': |
---|
1469 | rbid = str(rbids.index(RB['RBId'][0])) #for spin RBs |
---|
1470 | else: |
---|
1471 | rbid = str(rbids.index(RB['RBId'])) #others |
---|
1472 | pfxRB = pfx+'RB'+rbKey+'P' |
---|
1473 | pstr = ['x','y','z'] |
---|
1474 | ostr = ['a','i','j','k'] |
---|
1475 | Sytsym = G2spc.SytSym(RB['Orig'][0],SGData)[0] |
---|
1476 | xId,xCoef = G2spc.GetCSxinel(Sytsym)[:2] # gen origin site sym |
---|
1477 | equivs = {1:[],2:[],3:[]} |
---|
1478 | for i in range(3): |
---|
1479 | name = pfxRB+pstr[i]+':'+str(iRB)+':'+rbid |
---|
1480 | phaseDict[name] = RB['Orig'][0][i] |
---|
1481 | if RB['Orig'][1]: |
---|
1482 | if xId[i] > 0: |
---|
1483 | phaseVary += [name,] |
---|
1484 | equivs[xId[i]].append([name,xCoef[i]]) |
---|
1485 | else: |
---|
1486 | if symHold is not None: #variable is held due to symmetry |
---|
1487 | symHold.append(name) |
---|
1488 | G2mv.StoreHold(name,'In rigid body') |
---|
1489 | for equiv in equivs: |
---|
1490 | if len(equivs[equiv]) > 1: |
---|
1491 | name = equivs[equiv][0][0] |
---|
1492 | coef = equivs[equiv][0][1] |
---|
1493 | for eqv in equivs[equiv][1:]: |
---|
1494 | eqv[1] /= coef |
---|
1495 | G2mv.StoreEquivalence(name,(eqv,)) |
---|
1496 | pfxRB = pfx+'RB'+rbKey+'O' |
---|
1497 | A,V = G2mth.Q2AV(RB['Orig'][0]) |
---|
1498 | fixAxis = [0, np.abs(V).argmax()+1] |
---|
1499 | for i in range(4): |
---|
1500 | name = pfxRB+ostr[i]+':'+str(iRB)+':'+rbid |
---|
1501 | phaseDict[name] = RB['Orient'][0][i] |
---|
1502 | if RB['Orient'][1] == 'AV' and i: |
---|
1503 | phaseVary += [name,] |
---|
1504 | elif RB['Orient'][1] == 'A' and not i: |
---|
1505 | phaseVary += [name,] |
---|
1506 | elif RB['Orient'][1] == 'V' and i not in fixAxis: |
---|
1507 | phaseVary += [name,] |
---|
1508 | if rbKey != 'S': |
---|
1509 | name = pfx+'RB'+rbKey+'f:'+str(iRB)+':'+rbid |
---|
1510 | phaseDict[name] = RB['AtomFrac'][0] |
---|
1511 | if RB['AtomFrac'][1]: |
---|
1512 | phaseVary += [name,] |
---|
1513 | else: |
---|
1514 | name = pfx+'RB'+rbKey+'AtNo:'+str(iRB)+':'+rbid |
---|
1515 | phaseDict[name] = atomIndx[RB['Ids'][0]][1] |
---|
1516 | name = pfx+'RB'+rbKey+'symAxis:'+str(iRB)+':'+rbid |
---|
1517 | phaseDict[name] = RB['symAxis'] |
---|
1518 | name = pfx+'RB'+rbKey+'SytSym:'+str(iRB)+':'+rbid |
---|
1519 | phaseDict[name] = RB['SytSym'] |
---|
1520 | for ish,atype in enumerate(RB['atType']): |
---|
1521 | name = '%sRBSAtType;%d:%d:%s'%(pfx,ish,iRB,rbid) |
---|
1522 | phaseDict[name] = RB['atType'][ish] |
---|
1523 | name = '%sRBSNatoms;%d:%d:%s'%(pfx,ish,iRB,rbid) |
---|
1524 | phaseDict[name] = RB['Natoms'][ish] |
---|
1525 | |
---|
1526 | def MakeRBThermals(rbKey,phaseVary,phaseDict): |
---|
1527 | rbid = str(rbids.index(RB['RBId'])) |
---|
1528 | tlstr = ['11','22','33','12','13','23'] |
---|
1529 | sstr = ['12','13','21','23','31','32','AA','BB'] |
---|
1530 | if 'T' in RB['ThermalMotion'][0]: |
---|
1531 | pfxRB = pfx+'RB'+rbKey+'T' |
---|
1532 | for i in range(6): |
---|
1533 | name = pfxRB+tlstr[i]+':'+str(iRB)+':'+rbid |
---|
1534 | phaseDict[name] = RB['ThermalMotion'][1][i] |
---|
1535 | if RB['ThermalMotion'][2][i]: |
---|
1536 | phaseVary += [name,] |
---|
1537 | if 'L' in RB['ThermalMotion'][0]: |
---|
1538 | pfxRB = pfx+'RB'+rbKey+'L' |
---|
1539 | for i in range(6): |
---|
1540 | name = pfxRB+tlstr[i]+':'+str(iRB)+':'+rbid |
---|
1541 | phaseDict[name] = RB['ThermalMotion'][1][i+6] |
---|
1542 | if RB['ThermalMotion'][2][i+6]: |
---|
1543 | phaseVary += [name,] |
---|
1544 | if 'S' in RB['ThermalMotion'][0]: |
---|
1545 | pfxRB = pfx+'RB'+rbKey+'S' |
---|
1546 | for i in range(8): |
---|
1547 | name = pfxRB+sstr[i]+':'+str(iRB)+':'+rbid |
---|
1548 | phaseDict[name] = RB['ThermalMotion'][1][i+12] |
---|
1549 | if RB['ThermalMotion'][2][i+12]: |
---|
1550 | phaseVary += [name,] |
---|
1551 | if 'U' in RB['ThermalMotion'][0]: |
---|
1552 | name = pfx+'RB'+rbKey+'U:'+str(iRB)+':'+rbid |
---|
1553 | phaseDict[name] = RB['ThermalMotion'][1][0] |
---|
1554 | if RB['ThermalMotion'][2][0]: |
---|
1555 | phaseVary += [name,] |
---|
1556 | |
---|
1557 | def MakeRBTorsions(rbKey,phaseVary,phaseDict): |
---|
1558 | rbid = str(rbids.index(RB['RBId'])) |
---|
1559 | pfxRB = pfx+'RB'+rbKey+'Tr;' |
---|
1560 | for i,torsion in enumerate(RB['Torsions']): |
---|
1561 | name = pfxRB+str(i)+':'+str(iRB)+':'+rbid |
---|
1562 | phaseDict[name] = torsion[0] |
---|
1563 | if torsion[1]: |
---|
1564 | phaseVary += [name,] |
---|
1565 | |
---|
1566 | def MakeRBSphHarm(rbKey,phaseVary,phaseDict): |
---|
1567 | for ish,Shcof in enumerate(RB['SHC']): |
---|
1568 | if not len(Shcof): |
---|
1569 | continue |
---|
1570 | shid = str(rbids.index(RB['RBId'][ish])) |
---|
1571 | pfxRB = '%sRBSSh;%d;'%(pfx,ish) |
---|
1572 | for i,shcof in enumerate(Shcof): |
---|
1573 | SHcof = Shcof[shcof] |
---|
1574 | name = pfxRB+shcof.strip('+').strip('-')+':'+str(iRB)+':'+shid #patch; remove sign |
---|
1575 | phaseDict[name] = SHcof[0]*SHcof[1] #apply sign p |
---|
1576 | if SHcof[2]: |
---|
1577 | phaseVary += [name,] |
---|
1578 | |
---|
1579 | if Print and pFile is None: raise Exception("specify pFile or Print=False") |
---|
1580 | if Print: |
---|
1581 | pFile.write('\n Phases:\n') |
---|
1582 | phaseVary = [] |
---|
1583 | phaseDict = {} |
---|
1584 | pawleyLookup = {} |
---|
1585 | FFtables = {} #scattering factors - xrays |
---|
1586 | EFtables = {} #scattering factors - electrons |
---|
1587 | MFtables = {} #Mag. form factors |
---|
1588 | BLtables = {} # neutrons |
---|
1589 | Natoms = {} |
---|
1590 | maxSSwave = {} |
---|
1591 | shModels = ['cylindrical','none','shear - 2/m','rolling - mmm'] |
---|
1592 | SamSym = dict(zip(shModels,['0','-1','2/m','mmm'])) |
---|
1593 | atomIndx = {} |
---|
1594 | for name in PhaseData: |
---|
1595 | if seqHistName is not None and seqHistName != 'All': # sequential: load only used phases |
---|
1596 | if seqHistName not in PhaseData[name]['Histograms']: continue |
---|
1597 | if not PhaseData[name]['Histograms'][seqHistName]['Use']: continue |
---|
1598 | General = PhaseData[name]['General'] |
---|
1599 | pId = PhaseData[name]['pId'] |
---|
1600 | pfx = str(pId)+'::' |
---|
1601 | FFtable = G2el.GetFFtable(General['AtomTypes']) |
---|
1602 | EFtable = G2el.GetEFFtable(General['AtomTypes']) |
---|
1603 | BLtable = G2el.GetBLtable(General) |
---|
1604 | FFtables.update(FFtable) |
---|
1605 | EFtables.update(EFtable) |
---|
1606 | BLtables.update(BLtable) |
---|
1607 | phaseDict[pfx+'isMag'] = False |
---|
1608 | SGData = General['SGData'] |
---|
1609 | SGtext,SGtable = G2spc.SGPrint(SGData) |
---|
1610 | if General['Type'] == 'magnetic': |
---|
1611 | MFtable = G2el.GetMFtable(General['AtomTypes'],General['Lande g']) |
---|
1612 | MFtables.update(MFtable) |
---|
1613 | phaseDict[pfx+'isMag'] = True |
---|
1614 | SpnFlp = SGData['SpnFlp'] |
---|
1615 | Atoms = PhaseData[name]['Atoms'] |
---|
1616 | if Atoms and not General.get('doPawley'): |
---|
1617 | cx,ct,cs,cia = General['AtomPtrs'] |
---|
1618 | AtLookup = G2mth.FillAtomLookUp(Atoms,cia+8) |
---|
1619 | PawleyRef = PhaseData[name].get('Pawley ref',[]) |
---|
1620 | cell = General['Cell'] |
---|
1621 | A = G2lat.cell2A(cell[1:7]) |
---|
1622 | phaseDict.update({pfx+'A0':A[0],pfx+'A1':A[1],pfx+'A2':A[2], |
---|
1623 | pfx+'A3':A[3],pfx+'A4':A[4],pfx+'A5':A[5],pfx+'Vol':G2lat.calc_V(A)}) |
---|
1624 | if cell[0]: |
---|
1625 | phaseVary += cellVary(pfx,SGData) #also fills in symmetry required constraints |
---|
1626 | SSGtext = [] #no superstructure |
---|
1627 | im = 0 |
---|
1628 | if General.get('Modulated',False): |
---|
1629 | im = 1 #refl offset |
---|
1630 | Vec,vRef,maxH = General['SuperVec'] |
---|
1631 | phaseDict.update({pfx+'mV0':Vec[0],pfx+'mV1':Vec[1],pfx+'mV2':Vec[2]}) |
---|
1632 | SSGData = General['SSGData'] |
---|
1633 | SSGtext,SSGtable = G2spc.SSGPrint(SGData,SSGData) |
---|
1634 | if vRef: |
---|
1635 | phaseVary += modVary(pfx,SSGData) |
---|
1636 | if Atoms and not General.get('doPawley'): |
---|
1637 | cia = General['AtomPtrs'][3] |
---|
1638 | for i,at in enumerate(Atoms): |
---|
1639 | atomIndx[at[cia+8]] = [pfx,i] #lookup table for restraints |
---|
1640 | resRBData = PhaseData[name]['RBModels'].get('Residue',[]) |
---|
1641 | if resRBData: |
---|
1642 | rbids = rbIds['Residue'] #NB: used in the MakeRB routines |
---|
1643 | for iRB,RB in enumerate(resRBData): |
---|
1644 | MakeRBParms('R',phaseVary,phaseDict) |
---|
1645 | MakeRBThermals('R',phaseVary,phaseDict) |
---|
1646 | MakeRBTorsions('R',phaseVary,phaseDict) |
---|
1647 | |
---|
1648 | vecRBData = PhaseData[name]['RBModels'].get('Vector',[]) |
---|
1649 | if vecRBData: |
---|
1650 | rbids = rbIds['Vector'] #NB: used in the MakeRB routines |
---|
1651 | for iRB,RB in enumerate(vecRBData): |
---|
1652 | MakeRBParms('V',phaseVary,phaseDict) |
---|
1653 | MakeRBThermals('V',phaseVary,phaseDict) |
---|
1654 | |
---|
1655 | spnRBData = PhaseData[name]['RBModels'].get('Spin',[]) |
---|
1656 | if spnRBData: |
---|
1657 | rbids = rbIds['Spin'] #NB: used in the MakeRB routines |
---|
1658 | for iRB,RB in enumerate(spnRBData): |
---|
1659 | MakeRBParms('S',phaseVary,phaseDict) |
---|
1660 | MakeRBSphHarm('S',phaseVary,phaseDict) |
---|
1661 | |
---|
1662 | Natoms[pfx] = 0 |
---|
1663 | maxSSwave[pfx] = {'Sfrac':0,'Spos':0,'Sadp':0,'Smag':0} |
---|
1664 | if Atoms and not General.get('doPawley'): |
---|
1665 | cx,ct,cs,cia = General['AtomPtrs'] |
---|
1666 | Natoms[pfx] = len(Atoms) |
---|
1667 | for i,at in enumerate(Atoms): |
---|
1668 | # atomIndx[at[cia+8]] = [pfx,i] #lookup table for restraints |
---|
1669 | phaseDict.update({pfx+'Atype:'+str(i):at[ct],pfx+'Afrac:'+str(i):at[cx+3],pfx+'Amul:'+str(i):at[cs+1], |
---|
1670 | pfx+'Ax:'+str(i):at[cx],pfx+'Ay:'+str(i):at[cx+1],pfx+'Az:'+str(i):at[cx+2], |
---|
1671 | pfx+'dAx:'+str(i):0.,pfx+'dAy:'+str(i):0.,pfx+'dAz:'+str(i):0., #refined shifts for x,y,z |
---|
1672 | pfx+'AI/A:'+str(i):at[cia],}) |
---|
1673 | if at[cia] == 'I': |
---|
1674 | phaseDict[pfx+'AUiso:'+str(i)] = at[cia+1] |
---|
1675 | else: |
---|
1676 | phaseDict.update({pfx+'AU11:'+str(i):at[cia+2],pfx+'AU22:'+str(i):at[cia+3],pfx+'AU33:'+str(i):at[cia+4], |
---|
1677 | pfx+'AU12:'+str(i):at[cia+5],pfx+'AU13:'+str(i):at[cia+6],pfx+'AU23:'+str(i):at[cia+7]}) |
---|
1678 | if General['Type'] == 'magnetic': |
---|
1679 | phaseDict.update({pfx+'AMx:'+str(i):at[cx+4],pfx+'AMy:'+str(i):at[cx+5],pfx+'AMz:'+str(i):at[cx+6]}) |
---|
1680 | if 'F' in at[ct+1]: |
---|
1681 | phaseVary.append(pfx+'Afrac:'+str(i)) |
---|
1682 | if 'X' in at[ct+1]: |
---|
1683 | try: #patch for sytsym name changes |
---|
1684 | xId,xCoef = G2spc.GetCSxinel(at[cs])[:2] |
---|
1685 | except KeyError: |
---|
1686 | Sytsym = G2spc.SytSym(at[cx:cx+3],SGData)[0] |
---|
1687 | at[cs] = Sytsym |
---|
1688 | xId,xCoef = G2spc.GetCSxinel(at[cs])[:2] |
---|
1689 | names = [pfx+'dAx:'+str(i),pfx+'dAy:'+str(i),pfx+'dAz:'+str(i)] |
---|
1690 | equivs = {1:[],2:[],3:[]} |
---|
1691 | for j in range(3): |
---|
1692 | if xId[j] > 0: |
---|
1693 | phaseVary.append(names[j]) |
---|
1694 | equivs[xId[j]].append([names[j],xCoef[j]]) |
---|
1695 | else: |
---|
1696 | if symHold is not None: #variable is held due to symmetry |
---|
1697 | symHold.append(names[j]) |
---|
1698 | G2mv.StoreHold(names[j],'Fixed by symmetry') |
---|
1699 | for equiv in equivs: |
---|
1700 | if len(equivs[equiv]) > 1: |
---|
1701 | name = equivs[equiv][0][0] |
---|
1702 | coef = equivs[equiv][0][1] |
---|
1703 | for eqv in equivs[equiv][1:]: |
---|
1704 | eqv[1] /= coef |
---|
1705 | G2mv.StoreEquivalence(name,(eqv,)) |
---|
1706 | if 'U' in at[ct+1]: |
---|
1707 | if at[cia] == 'I': |
---|
1708 | phaseVary.append(pfx+'AUiso:'+str(i)) |
---|
1709 | else: |
---|
1710 | try: #patch for sytsym name changes |
---|
1711 | uId,uCoef = G2spc.GetCSuinel(at[cs])[:2] |
---|
1712 | except KeyError: |
---|
1713 | Sytsym = G2spc.SytSym(at[cx:cx+3],SGData)[0] |
---|
1714 | at[cs] = Sytsym |
---|
1715 | uId,uCoef = G2spc.GetCSuinel(at[cs])[:2] |
---|
1716 | names = [pfx+'AU11:'+str(i),pfx+'AU22:'+str(i),pfx+'AU33:'+str(i), |
---|
1717 | pfx+'AU12:'+str(i),pfx+'AU13:'+str(i),pfx+'AU23:'+str(i)] |
---|
1718 | equivs = {1:[],2:[],3:[],4:[],5:[],6:[]} |
---|
1719 | for j in range(6): |
---|
1720 | if uId[j] > 0: |
---|
1721 | phaseVary.append(names[j]) |
---|
1722 | equivs[uId[j]].append([names[j],uCoef[j]]) |
---|
1723 | for equiv in equivs: |
---|
1724 | if len(equivs[equiv]) > 1: |
---|
1725 | name = equivs[equiv][0][0] |
---|
1726 | coef = equivs[equiv][0][1] |
---|
1727 | for eqv in equivs[equiv][1:]: |
---|
1728 | eqv[1] /= coef |
---|
1729 | G2mv.StoreEquivalence(name,(eqv,)) |
---|
1730 | if 'M' in at[ct+1]: |
---|
1731 | SytSym,Mul,Nop,dupDir = G2spc.SytSym(at[cx:cx+3],SGData) |
---|
1732 | mId,mCoef = G2spc.GetCSpqinel(SpnFlp,dupDir) |
---|
1733 | names = [pfx+'AMx:'+str(i),pfx+'AMy:'+str(i),pfx+'AMz:'+str(i)] |
---|
1734 | equivs = {1:[],2:[],3:[]} |
---|
1735 | for j in range(3): |
---|
1736 | if mId[j] > 0: |
---|
1737 | phaseVary.append(names[j]) |
---|
1738 | equivs[mId[j]].append([names[j],mCoef[j]]) |
---|
1739 | for equiv in equivs: |
---|
1740 | if len(equivs[equiv]) > 1: |
---|
1741 | name = equivs[equiv][0][0] |
---|
1742 | coef = equivs[equiv][0][1] |
---|
1743 | for eqv in equivs[equiv][1:]: |
---|
1744 | eqv[1] /= coef |
---|
1745 | G2mv.StoreEquivalence(name,(eqv,)) |
---|
1746 | if General.get('Modulated',False): |
---|
1747 | AtomSS = at[-1]['SS1'] |
---|
1748 | for Stype in ['Sfrac','Spos','Sadp','Smag']: |
---|
1749 | Waves = AtomSS[Stype] |
---|
1750 | if len(Waves): |
---|
1751 | waveType = Waves[0] |
---|
1752 | else: |
---|
1753 | continue |
---|
1754 | phaseDict[pfx+Stype[1].upper()+'waveType:'+str(i)] = waveType |
---|
1755 | nx = 0 |
---|
1756 | for iw,wave in enumerate(Waves[1:]): |
---|
1757 | if not iw: |
---|
1758 | if waveType in ['ZigZag','Block']: |
---|
1759 | nx = 1 |
---|
1760 | CSI = G2spc.GetSSfxuinel(waveType,Stype,1,at[cx:cx+3],SGData,SSGData) |
---|
1761 | else: |
---|
1762 | CSI = G2spc.GetSSfxuinel('Fourier',Stype,iw+1-nx,at[cx:cx+3],SGData,SSGData) |
---|
1763 | uId,uCoef = CSI[0] |
---|
1764 | stiw = str(i)+':'+str(iw) |
---|
1765 | if Stype == 'Spos': |
---|
1766 | if waveType in ['ZigZag','Block',] and not iw: |
---|
1767 | names = [pfx+'Tmin:'+stiw,pfx+'Tmax:'+stiw,pfx+'Xmax:'+stiw,pfx+'Ymax:'+stiw,pfx+'Zmax:'+stiw] |
---|
1768 | equivs = {1:[],2:[], 3:[],4:[],5:[]} |
---|
1769 | else: |
---|
1770 | names = [pfx+'Xsin:'+stiw,pfx+'Ysin:'+stiw,pfx+'Zsin:'+stiw, |
---|
1771 | pfx+'Xcos:'+stiw,pfx+'Ycos:'+stiw,pfx+'Zcos:'+stiw] |
---|
1772 | equivs = {1:[],2:[],3:[], 4:[],5:[],6:[]} |
---|
1773 | elif Stype == 'Sadp': |
---|
1774 | names = [pfx+'U11sin:'+stiw,pfx+'U22sin:'+stiw,pfx+'U33sin:'+stiw, |
---|
1775 | pfx+'U12sin:'+stiw,pfx+'U13sin:'+stiw,pfx+'U23sin:'+stiw, |
---|
1776 | pfx+'U11cos:'+stiw,pfx+'U22cos:'+stiw,pfx+'U33cos:'+stiw, |
---|
1777 | pfx+'U12cos:'+stiw,pfx+'U13cos:'+stiw,pfx+'U23cos:'+stiw] |
---|
1778 | equivs = {1:[],2:[],3:[],4:[],5:[],6:[], 7:[],8:[],9:[],10:[],11:[],12:[]} |
---|
1779 | elif Stype == 'Sfrac': |
---|
1780 | equivs = {1:[],2:[]} |
---|
1781 | if 'Crenel' in waveType and not iw: |
---|
1782 | names = [pfx+'Fzero:'+stiw,pfx+'Fwid:'+stiw] |
---|
1783 | else: |
---|
1784 | names = [pfx+'Fsin:'+stiw,pfx+'Fcos:'+stiw] |
---|
1785 | elif Stype == 'Smag': |
---|
1786 | equivs = {1:[],2:[],3:[], 4:[],5:[],6:[]} |
---|
1787 | names = [pfx+'MXsin:'+stiw,pfx+'MYsin:'+stiw,pfx+'MZsin:'+stiw, |
---|
1788 | pfx+'MXcos:'+stiw,pfx+'MYcos:'+stiw,pfx+'MZcos:'+stiw] |
---|
1789 | phaseDict.update(dict(zip(names,wave[0]))) |
---|
1790 | if wave[1]: #what do we do here for multiple terms in modulation constraints? |
---|
1791 | for j in range(len(equivs)): |
---|
1792 | if uId[j][0] > 0: |
---|
1793 | phaseVary.append(names[j]) |
---|
1794 | equivs[uId[j][0]].append([names[j],uCoef[j][0]]) |
---|
1795 | for equiv in equivs: |
---|
1796 | if len(equivs[equiv]) > 1: |
---|
1797 | name = equivs[equiv][0][0] |
---|
1798 | coef = equivs[equiv][0][1] |
---|
1799 | for eqv in equivs[equiv][1:]: |
---|
1800 | eqv[1] /= coef |
---|
1801 | G2mv.StoreEquivalence(name,(eqv,)) |
---|
1802 | maxSSwave[pfx][Stype] = max(maxSSwave[pfx][Stype],iw+1) |
---|
1803 | textureData = General['SH Texture'] |
---|
1804 | if textureData['Order']: # and seqHistName is not None: |
---|
1805 | phaseDict[pfx+'SHorder'] = textureData['Order'] |
---|
1806 | phaseDict[pfx+'SHmodel'] = SamSym[textureData['Model']] |
---|
1807 | for item in ['omega','chi','phi']: |
---|
1808 | phaseDict[pfx+'SH '+item] = textureData['Sample '+item][1] |
---|
1809 | if textureData['Sample '+item][0]: |
---|
1810 | phaseVary.append(pfx+'SH '+item) |
---|
1811 | for item in textureData['SH Coeff'][1]: |
---|
1812 | phaseDict[pfx+item] = textureData['SH Coeff'][1][item] |
---|
1813 | if textureData['SH Coeff'][0]: |
---|
1814 | phaseVary.append(pfx+item) |
---|
1815 | |
---|
1816 | if Print: |
---|
1817 | pFile.write('\n Phase name: %s\n'%General['Name']) |
---|
1818 | pFile.write(135*'='+'\n') |
---|
1819 | PrintFFtable(FFtable) |
---|
1820 | PrintEFtable(EFtable) |
---|
1821 | PrintBLtable(BLtable) |
---|
1822 | if General['Type'] == 'magnetic': |
---|
1823 | PrintMFtable(MFtable) |
---|
1824 | pFile.write('\n') |
---|
1825 | #how do we print magnetic symmetry table? TBD |
---|
1826 | if len(SSGtext): #if superstructure |
---|
1827 | for line in SSGtext: pFile.write(line+'\n') |
---|
1828 | if len(SSGtable): |
---|
1829 | for item in SSGtable: |
---|
1830 | line = ' %s '%(item) |
---|
1831 | pFile.write(line+'\n') |
---|
1832 | else: |
---|
1833 | pFile.write(' ( 1) %s\n'%(SSGtable[0])) |
---|
1834 | else: |
---|
1835 | for line in SGtext: pFile.write(line+'\n') |
---|
1836 | if len(SGtable): |
---|
1837 | for item in SGtable: |
---|
1838 | line = ' %s '%(item) |
---|
1839 | pFile.write(line+'\n') |
---|
1840 | else: |
---|
1841 | pFile.write(' ( 1) %s\n'%(SGtable[0])) |
---|
1842 | PrintRBObjects(resRBData,vecRBData,spnRBData) |
---|
1843 | PrintAtoms(General,Atoms) |
---|
1844 | if General['Type'] == 'magnetic': |
---|
1845 | PrintMoments(General,Atoms) |
---|
1846 | if General.get('Modulated',False): |
---|
1847 | PrintWaves(General,Atoms) |
---|
1848 | pFile.write('\n Unit cell: a = %.5f b = %.5f c = %.5f alpha = %.3f beta = %.3f gamma = %.3f volume = %.3f Refine? %s\n'% |
---|
1849 | (cell[1],cell[2],cell[3],cell[4],cell[5],cell[6],cell[7],cell[0])) |
---|
1850 | if len(SSGtext): #if superstructure |
---|
1851 | pFile.write('\n Modulation vector: mV0 = %.4f mV1 = %.4f mV2 = %.4f max mod. index = %d Refine? %s\n'% |
---|
1852 | (Vec[0],Vec[1],Vec[2],maxH,vRef)) |
---|
1853 | if seqHistName is not None: |
---|
1854 | PrintTexture(textureData) |
---|
1855 | if name in RestraintDict: |
---|
1856 | PrintRestraints(cell[1:7],SGData,General['AtomPtrs'],Atoms,AtLookup, |
---|
1857 | textureData,RestraintDict[name],pFile) |
---|
1858 | |
---|
1859 | elif PawleyRef: |
---|
1860 | if Print: |
---|
1861 | pFile.write('\n Phase name: %s\n'%General['Name']) |
---|
1862 | pFile.write(135*'='+'\n') |
---|
1863 | pFile.write('\n') |
---|
1864 | if len(SSGtext): #if superstructure |
---|
1865 | for line in SSGtext: pFile.write(line+'\n') |
---|
1866 | if len(SSGtable): |
---|
1867 | for item in SSGtable: |
---|
1868 | line = ' %s '%(item) |
---|
1869 | pFile.write(line+'\n') |
---|
1870 | else: |
---|
1871 | pFile.write(' ( 1) %s\n'%SSGtable[0]) |
---|
1872 | else: |
---|
1873 | for line in SGtext: pFile.write(line+'\n') |
---|
1874 | if len(SGtable): |
---|
1875 | for item in SGtable: |
---|
1876 | line = ' %s '%(item) |
---|
1877 | pFile.write(line+'\n') |
---|
1878 | else: |
---|
1879 | pFile.write(' ( 1) %s\n'%(SGtable[0])) |
---|
1880 | pFile.write('\n Unit cell: a = %.5f b = %.5f c = %.5f alpha = %.3f beta = %.3f gamma = %.3f volume = %.3f Refine? %s\n'% |
---|
1881 | (cell[1],cell[2],cell[3],cell[4],cell[5],cell[6],cell[7],cell[0])) |
---|
1882 | if len(SSGtext): #if superstructure |
---|
1883 | pFile.write('\n Modulation vector: mV0 = %.4f mV1 = %.4f mV2 = %.4f max mod. index = %d Refine? %s\n'% |
---|
1884 | (Vec[0],Vec[1],Vec[2],maxH,vRef)) |
---|
1885 | pawleyVary = [] |
---|
1886 | for i,refl in enumerate(PawleyRef): |
---|
1887 | phaseDict[pfx+'PWLref:'+str(i)] = refl[6+im] |
---|
1888 | if im: |
---|
1889 | pawleyLookup[pfx+'%d,%d,%d,%d'%(refl[0],refl[1],refl[2],refl[3])] = i |
---|
1890 | else: |
---|
1891 | pawleyLookup[pfx+'%d,%d,%d'%(refl[0],refl[1],refl[2])] = i |
---|
1892 | if refl[5+im]: |
---|
1893 | pawleyVary.append(pfx+'PWLref:'+str(i)) |
---|
1894 | GetPawleyConstr(SGData['SGLaue'],PawleyRef,im,pawleyVary) #does G2mv.StoreEquivalence |
---|
1895 | phaseVary += pawleyVary |
---|
1896 | |
---|
1897 | return Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,FFtables,EFtables,BLtables,MFtables,maxSSwave |
---|
1898 | |
---|
1899 | def cellFill(pfx,SGData,parmDict,sigDict): |
---|
1900 | '''Returns the filled-out reciprocal cell (A) terms and their uncertainties |
---|
1901 | from the parameter and sig dictionaries. |
---|
1902 | |
---|
1903 | :param str pfx: parameter prefix ("n::", where n is a phase number) |
---|
1904 | :param dict SGdata: a symmetry object |
---|
1905 | :param dict parmDict: a dictionary of parameters |
---|
1906 | :param dict sigDict: a dictionary of uncertainties on parameters |
---|
1907 | |
---|
1908 | :returns: A,sigA where each is a list of six terms with the A terms |
---|
1909 | ''' |
---|
1910 | if SGData['SGLaue'] in ['-1',]: |
---|
1911 | A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'], |
---|
1912 | parmDict[pfx+'A3'],parmDict[pfx+'A4'],parmDict[pfx+'A5']] |
---|
1913 | elif SGData['SGLaue'] in ['2/m',]: |
---|
1914 | if SGData['SGUniq'] == 'a': |
---|
1915 | A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'], |
---|
1916 | 0,0,parmDict[pfx+'A5']] |
---|
1917 | elif SGData['SGUniq'] == 'b': |
---|
1918 | A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'], |
---|
1919 | 0,parmDict[pfx+'A4'],0] |
---|
1920 | else: |
---|
1921 | A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'], |
---|
1922 | parmDict[pfx+'A3'],0,0] |
---|
1923 | elif SGData['SGLaue'] in ['mmm',]: |
---|
1924 | A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],0,0,0] |
---|
1925 | elif SGData['SGLaue'] in ['4/m','4/mmm']: |
---|
1926 | A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A2'],0,0,0] |
---|
1927 | elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']: |
---|
1928 | A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A2'], |
---|
1929 | parmDict[pfx+'A0'],0,0] |
---|
1930 | elif SGData['SGLaue'] in ['3R', '3mR']: |
---|
1931 | A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A0'], |
---|
1932 | parmDict[pfx+'A3'],parmDict[pfx+'A3'],parmDict[pfx+'A3']] |
---|
1933 | elif SGData['SGLaue'] in ['m3m','m3']: |
---|
1934 | A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A0'],0,0,0] |
---|
1935 | |
---|
1936 | try: |
---|
1937 | if SGData['SGLaue'] in ['-1',]: |
---|
1938 | sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'], |
---|
1939 | sigDict[pfx+'A3'],sigDict[pfx+'A4'],sigDict[pfx+'A5']] |
---|
1940 | elif SGData['SGLaue'] in ['2/m',]: |
---|
1941 | if SGData['SGUniq'] == 'a': |
---|
1942 | sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'], |
---|
1943 | 0,0,sigDict[pfx+'A5']] |
---|
1944 | elif SGData['SGUniq'] == 'b': |
---|
1945 | sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'], |
---|
1946 | 0,sigDict[pfx+'A4'],0] |
---|
1947 | else: |
---|
1948 | sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'], |
---|
1949 | sigDict[pfx+'A3'],0,0] |
---|
1950 | elif SGData['SGLaue'] in ['mmm',]: |
---|
1951 | sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],0,0,0] |
---|
1952 | elif SGData['SGLaue'] in ['4/m','4/mmm']: |
---|
1953 | sigA = [sigDict[pfx+'A0'],0,sigDict[pfx+'A2'],0,0,0] |
---|
1954 | elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']: |
---|
1955 | sigA = [sigDict[pfx+'A0'],0,sigDict[pfx+'A2'],0,0,0] |
---|
1956 | elif SGData['SGLaue'] in ['3R', '3mR']: |
---|
1957 | sigA = [sigDict[pfx+'A0'],0,0,sigDict[pfx+'A3'],0,0] |
---|
1958 | elif SGData['SGLaue'] in ['m3m','m3']: |
---|
1959 | sigA = [sigDict[pfx+'A0'],0,0,0,0,0] |
---|
1960 | except KeyError: |
---|
1961 | sigA = [0,0,0,0,0,0] |
---|
1962 | return A,sigA |
---|
1963 | |
---|
1964 | def PrintRestraints(cell,SGData,AtPtrs,Atoms,AtLookup,textureData,phaseRest,pFile): |
---|
1965 | 'needs a doc string' |
---|
1966 | if phaseRest: |
---|
1967 | Amat = G2lat.cell2AB(cell)[0] |
---|
1968 | cx,ct,cs = AtPtrs[:3] |
---|
1969 | names = [['Bond','Bonds'],['Angle','Angles'],['Plane','Planes'], |
---|
1970 | ['Chiral','Volumes'],['Torsion','Torsions'],['Rama','Ramas'], |
---|
1971 | ['ChemComp','Sites'],['Texture','HKLs']] |
---|
1972 | for name,rest in names: |
---|
1973 | itemRest = phaseRest[name] |
---|
1974 | if rest in itemRest and itemRest[rest] and itemRest['Use']: |
---|
1975 | pFile.write('\n %s restraint weight factor %10.3f Use: %s\n'%(name,itemRest['wtFactor'],str(itemRest['Use']))) |
---|
1976 | if name in ['Bond','Angle','Plane','Chiral']: |
---|
1977 | pFile.write(' calc obs sig delt/sig atoms(symOp): \n') |
---|
1978 | for indx,ops,obs,esd in itemRest[rest]: |
---|
1979 | try: |
---|
1980 | AtNames = G2mth.GetAtomItemsById(Atoms,AtLookup,indx,ct-1) |
---|
1981 | AtName = '' |
---|
1982 | for i,Aname in enumerate(AtNames): |
---|
1983 | AtName += Aname |
---|
1984 | if ops[i] == '1': |
---|
1985 | AtName += '-' |
---|
1986 | else: |
---|
1987 | AtName += '+('+ops[i]+')-' |
---|
1988 | XYZ = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cx,3)) |
---|
1989 | XYZ = G2mth.getSyXYZ(XYZ,ops,SGData) |
---|
1990 | if name == 'Bond': |
---|
1991 | calc = G2mth.getRestDist(XYZ,Amat) |
---|
1992 | elif name == 'Angle': |
---|
1993 | calc = G2mth.getRestAngle(XYZ,Amat) |
---|
1994 | elif name == 'Plane': |
---|
1995 | calc = G2mth.getRestPlane(XYZ,Amat) |
---|
1996 | elif name == 'Chiral': |
---|
1997 | calc = G2mth.getRestChiral(XYZ,Amat) |
---|
1998 | pFile.write(' %9.3f %9.3f %8.3f %8.3f %s\n'%(calc,obs,esd,(obs-calc)/esd,AtName[:-1])) |
---|
1999 | except KeyError: |
---|
2000 | del itemRest[rest] |
---|
2001 | elif name in ['Torsion','Rama']: |
---|
2002 | pFile.write(' atoms(symOp) calc obs sig delt/sig torsions: \n') |
---|
2003 | coeffDict = itemRest['Coeff'] |
---|
2004 | for indx,ops,cofName,esd in itemRest[rest]: |
---|
2005 | AtNames = G2mth.GetAtomItemsById(Atoms,AtLookup,indx,ct-1) |
---|
2006 | AtName = '' |
---|
2007 | for i,Aname in enumerate(AtNames): |
---|
2008 | AtName += Aname+'+('+ops[i]+')-' |
---|
2009 | XYZ = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cx,3)) |
---|
2010 | XYZ = G2mth.getSyXYZ(XYZ,ops,SGData) |
---|
2011 | if name == 'Torsion': |
---|
2012 | tor = G2mth.getRestTorsion(XYZ,Amat) |
---|
2013 | restr,calc = G2mth.calcTorsionEnergy(tor,coeffDict[cofName]) |
---|
2014 | pFile.write(' %8.3f %8.3f %.3f %8.3f %8.3f %s\n'%(calc,obs,esd,(obs-calc)/esd,tor,AtName[:-1])) |
---|
2015 | else: |
---|
2016 | phi,psi = G2mth.getRestRama(XYZ,Amat) |
---|
2017 | restr,calc = G2mth.calcRamaEnergy(phi,psi,coeffDict[cofName]) |
---|
2018 | pFile.write(' %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %s\n'%(calc,obs,esd,(obs-calc)/esd,phi,psi,AtName[:-1])) |
---|
2019 | elif name == 'ChemComp': |
---|
2020 | pFile.write(' atoms mul*frac factor prod\n') |
---|
2021 | for indx,factors,obs,esd in itemRest[rest]: |
---|
2022 | try: |
---|
2023 | atoms = G2mth.GetAtomItemsById(Atoms,AtLookup,indx,ct-1) |
---|
2024 | mul = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cs+1)) |
---|
2025 | frac = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cs-1)) |
---|
2026 | mulfrac = mul*frac |
---|
2027 | calcs = mul*frac*factors |
---|
2028 | for iatm,[atom,mf,fr,clc] in enumerate(zip(atoms,mulfrac,factors,calcs)): |
---|
2029 | pFile.write(' %10s %8.3f %8.3f %8.3f\n'%(atom,mf,fr,clc)) |
---|
2030 | pFile.write(' Sum: calc: %8.3f obs: %8.3f esd: %8.3f\n'%(np.sum(calcs),obs,esd)) |
---|
2031 | except KeyError: |
---|
2032 | del itemRest[rest] |
---|
2033 | elif name == 'Texture' and textureData['Order']: |
---|
2034 | SHCoef = textureData['SH Coeff'][1] |
---|
2035 | shModels = ['cylindrical','none','shear - 2/m','rolling - mmm'] |
---|
2036 | SamSym = dict(zip(shModels,['0','-1','2/m','mmm'])) |
---|
2037 | pFile.write (' HKL grid neg esd sum neg num neg use unit? unit esd \n') |
---|
2038 | for hkl,grid,esd1,ifesd2,esd2 in itemRest[rest]: |
---|
2039 | phi,beta = G2lat.CrsAng(np.array(hkl),cell,SGData) |
---|
2040 | ODFln = G2lat.Flnh(SHCoef,phi,beta,SGData) |
---|
2041 | R,P,Z = G2mth.getRestPolefig(ODFln,SamSym[textureData['Model']],grid) |
---|
2042 | Z = ma.masked_greater(Z,0.0) |
---|
2043 | num = ma.count(Z) |
---|
2044 | sum = 0 |
---|
2045 | if num: |
---|
2046 | sum = np.sum(Z) |
---|
2047 | pFile.write (' %d %d %d %d %8.3f %8.3f %8d %s %8.3f\n'%(hkl[0],hkl[1],hkl[2],grid,esd1,sum,num,str(ifesd2),esd2)) |
---|
2048 | |
---|
2049 | def getCellEsd(pfx,SGData,A,covData,unique=False): |
---|
2050 | '''Compute the standard uncertainty on cell parameters |
---|
2051 | |
---|
2052 | :param str pfx: prefix of form p\\:\\: |
---|
2053 | :param SGdata: space group information |
---|
2054 | :param list A: Reciprocal cell Ai terms |
---|
2055 | :param dict covData: covariance tree item |
---|
2056 | :param bool unique: when True, only directly refined parameters |
---|
2057 | (a in cubic, a & alpha in rhombohedral cells) are assigned |
---|
2058 | positive s.u. values. Used for CIF generation. |
---|
2059 | ''' |
---|
2060 | rVsq = G2lat.calc_rVsq(A) |
---|
2061 | G,g = G2lat.A2Gmat(A) #get recip. & real metric tensors |
---|
2062 | RMnames = [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3',pfx+'A4',pfx+'A5'] |
---|
2063 | varyList = covData['varyList'] |
---|
2064 | covMatrix = covData['covMatrix'] |
---|
2065 | if len(covMatrix): |
---|
2066 | vcov = G2mth.getVCov(RMnames,varyList,covMatrix) |
---|
2067 | if SGData['SGLaue'] in ['3', '3m1', '31m', '6/m', '6/mmm']: |
---|
2068 | vcov[1,1] = vcov[3,3] = vcov[0,1] = vcov[1,0] = vcov[0,0] |
---|
2069 | vcov[1,3] = vcov[3,1] = vcov[0,3] = vcov[3,0] = vcov[0,0] |
---|
2070 | vcov[1,2] = vcov[2,1] = vcov[2,3] = vcov[3,2] = vcov[0,2] |
---|
2071 | elif SGData['SGLaue'] in ['m3','m3m']: |
---|
2072 | vcov[0:3,0:3] = vcov[0,0] |
---|
2073 | elif SGData['SGLaue'] in ['4/m', '4/mmm']: |
---|
2074 | vcov[0:2,0:2] = vcov[0,0] |
---|
2075 | vcov[1,2] = vcov[2,1] = vcov[0,2] |
---|
2076 | elif SGData['SGLaue'] in ['3R','3mR']: |
---|
2077 | vcov[0:3,0:3] = vcov[0,0] |
---|
2078 | # vcov[4,4] = vcov[5,5] = vcov[3,3] |
---|
2079 | vcov[3:6,3:6] = vcov[3,3] |
---|
2080 | vcov[0:3,3:6] = vcov[0,3] |
---|
2081 | vcov[3:6,0:3] = vcov[3,0] |
---|
2082 | else: |
---|
2083 | vcov = np.eye(6) |
---|
2084 | delt = 1.e-9 |
---|
2085 | drVdA = np.zeros(6) |
---|
2086 | for i in range(6): |
---|
2087 | A[i] += delt |
---|
2088 | drVdA[i] = G2lat.calc_rVsq(A) |
---|
2089 | A[i] -= 2*delt |
---|
2090 | drVdA[i] -= G2lat.calc_rVsq(A) |
---|
2091 | A[i] += delt |
---|
2092 | drVdA /= 2.*delt |
---|
2093 | srcvlsq = np.inner(drVdA,np.inner(drVdA,vcov)) |
---|
2094 | Vol = 1/np.sqrt(rVsq) |
---|
2095 | sigVol = Vol**3*np.sqrt(srcvlsq)/2. #ok - checks with GSAS |
---|
2096 | |
---|
2097 | dcdA = np.zeros((6,6)) |
---|
2098 | for i in range(6): |
---|
2099 | pdcdA =np.zeros(6) |
---|
2100 | A[i] += delt |
---|
2101 | pdcdA += G2lat.A2cell(A) |
---|
2102 | A[i] -= 2*delt |
---|
2103 | pdcdA -= G2lat.A2cell(A) |
---|
2104 | A[i] += delt |
---|
2105 | dcdA[i] = pdcdA/(2.*delt) |
---|
2106 | |
---|
2107 | sigMat = np.inner(dcdA,np.inner(dcdA,vcov)) |
---|
2108 | var = np.diag(sigMat) |
---|
2109 | CS = np.where(var>0.,np.sqrt(var),0.) |
---|
2110 | if SGData['SGLaue'] in ['3', '3m1', '31m', '6/m', '6/mmm','m3','m3m','4/m','4/mmm']: |
---|
2111 | CS[3:6] = 0.0 |
---|
2112 | # show s.u. values only for the unique values |
---|
2113 | if not unique: |
---|
2114 | pass |
---|
2115 | elif SGData['SGLaue'] in ['3', '3m1', '31m', '6/m', '6/mmm','4/m', '4/mmm']: |
---|
2116 | CS[1] = -CS[1] |
---|
2117 | elif SGData['SGLaue'] in ['m3','m3m']: |
---|
2118 | CS[1] = -CS[1] |
---|
2119 | CS[2] = -CS[2] |
---|
2120 | elif SGData['SGLaue'] in ['3R','3mR']: |
---|
2121 | CS[1] = -CS[1] |
---|
2122 | CS[2] = -CS[2] |
---|
2123 | CS[4] = -CS[4] |
---|
2124 | CS[3] = -CS[3] |
---|
2125 | return [CS[0],CS[1],CS[2],CS[5],CS[4],CS[3],sigVol] |
---|
2126 | |
---|
2127 | def getCellSU(pId,hId,SGData,parmDict,covData): |
---|
2128 | '''Compute the unit cell parameters and standard uncertainties |
---|
2129 | where lattice parameters and Hstrain (Dij) may be refined. This is |
---|
2130 | called only for generation of CIFs. |
---|
2131 | |
---|
2132 | :param pId: phase index |
---|
2133 | :param hId: histogram index |
---|
2134 | :param SGdata: space group information |
---|
2135 | :param dict parmDict: parameter dict, must have all non-zero Dij and Ai terms |
---|
2136 | :param dict covData: covariance tree item |
---|
2137 | ''' |
---|
2138 | |
---|
2139 | Dnames = ['{}:{}:D{}'.format(pId,hId,i) for i in ['11','22','33','12','13','23']] |
---|
2140 | Anames = ['{}::A{}'.format(pId,i) for i in range(6)] |
---|
2141 | Ai = [parmDict[i] for i in Anames] |
---|
2142 | Dij = [parmDict.get(i,0.) for i in Dnames] |
---|
2143 | A = np.array(Ai) + np.array(Dij) |
---|
2144 | cell = list(G2lat.A2cell(A)) + [G2lat.calc_V(A)] |
---|
2145 | rVsq = G2lat.calc_rVsq(A) |
---|
2146 | G,g = G2lat.A2Gmat(A) #get recip. & real metric tensors |
---|
2147 | varyList = covData['varyList'] |
---|
2148 | covMatrix = covData['covMatrix'] |
---|
2149 | if len(covMatrix): |
---|
2150 | vcov = G2mth.getVCov(Anames+Dnames,varyList,covMatrix) |
---|
2151 | for i in [0,6]: |
---|
2152 | for j in [0,6]: |
---|
2153 | if SGData['SGLaue'] in ['3', '3m1', '31m', '6/m', '6/mmm']: |
---|
2154 | vcov[1+i,1+j] = vcov[3+i,3+j] = vcov[i,1+j] = vcov[1+i,j] = vcov[i,j] |
---|
2155 | vcov[1+i,3+j] = vcov[3+i,1+j] = vcov[i,3+j] = vcov[3+i,j] = vcov[i,j] |
---|
2156 | vcov[1+i,2+j] = vcov[2+i,1+j] = vcov[2+i,3+j] = vcov[3+i,2+j] = vcov[i,2+j] |
---|
2157 | elif SGData['SGLaue'] in ['m3','m3m']: |
---|
2158 | vcov[i:3+i,j:3+j] = vcov[i,j] |
---|
2159 | elif SGData['SGLaue'] in ['4/m', '4/mmm']: |
---|
2160 | vcov[i:2+i,j:2+j] = vcov[i,j] |
---|
2161 | vcov[1+i,2+j] = vcov[2+i,1+j] = vcov[i,2+j] |
---|
2162 | elif SGData['SGLaue'] in ['3R','3mR']: |
---|
2163 | vcov[i:3+j,i:3+j] = vcov[i,j] |
---|
2164 | # vcov[4,4] = vcov[5,5] = vcov[3,3] |
---|
2165 | vcov[3+i:6+i,3+j:6+j] = vcov[3,3+j] |
---|
2166 | vcov[i:3+i,3+j:6+j] = vcov[i,3+j] |
---|
2167 | vcov[3+i:6+i,j:3+j] = vcov[3+i,j] |
---|
2168 | else: |
---|
2169 | vcov = np.eye(12) |
---|
2170 | delt = 1.e-9 |
---|
2171 | drVdA = np.zeros(12) |
---|
2172 | for i in range(12): |
---|
2173 | A[i%6] += delt |
---|
2174 | drVdA[i] = G2lat.calc_rVsq(A) |
---|
2175 | A[i%6] -= 2*delt |
---|
2176 | drVdA[i] -= G2lat.calc_rVsq(A) |
---|
2177 | A[i%6] += delt |
---|
2178 | drVdA /= 2.*delt |
---|
2179 | srcvlsq = np.inner(drVdA,np.inner(drVdA,vcov)) |
---|
2180 | Vol = 1/np.sqrt(rVsq) |
---|
2181 | sigVol = Vol**3*np.sqrt(srcvlsq)/2. #ok - checks with GSAS |
---|
2182 | |
---|
2183 | dcdA = np.zeros((12,12)) |
---|
2184 | for i in range(12): |
---|
2185 | pdcdA =np.zeros(12) |
---|
2186 | A[i%6] += delt |
---|
2187 | pdcdA += G2lat.A2cell(A)+G2lat.A2cell(A) |
---|
2188 | A[i%6] -= 2*delt |
---|
2189 | pdcdA -= G2lat.A2cell(A)+G2lat.A2cell(A) |
---|
2190 | A[i%6] += delt |
---|
2191 | dcdA[i] = pdcdA/(2.*delt) |
---|
2192 | sigMat = np.inner(dcdA,np.inner(dcdA,vcov)) |
---|
2193 | var = np.diag(sigMat) |
---|
2194 | CS = np.where(var>0.,np.sqrt(var),0.) |
---|
2195 | if SGData['SGLaue'] in ['3', '3m1', '31m', '6/m', '6/mmm','m3','m3m','4/m','4/mmm']: |
---|
2196 | CS[3:6] = 0.0 |
---|
2197 | # show s.u. values only for the unique values |
---|
2198 | if SGData['SGLaue'] in ['3', '3m1', '31m', '6/m', '6/mmm','4/m', '4/mmm']: |
---|
2199 | CS[1] = -CS[1] |
---|
2200 | elif SGData['SGLaue'] in ['m3','m3m']: |
---|
2201 | CS[1] = -CS[1] |
---|
2202 | CS[2] = -CS[2] |
---|
2203 | elif SGData['SGLaue'] in ['3R','3mR']: |
---|
2204 | CS[1] = -CS[1] |
---|
2205 | CS[2] = -CS[2] |
---|
2206 | CS[4] = -CS[4] |
---|
2207 | CS[3] = -CS[3] |
---|
2208 | return cell,[CS[0],CS[1],CS[2],CS[5],CS[4],CS[3],sigVol] |
---|
2209 | |
---|
2210 | def SetPhaseData(parmDict,sigDict,Phases,RBIds,covData,RestraintDict=None,pFile=None): |
---|
2211 | '''Called after a refinement to transfer parameters from the parameter dict to |
---|
2212 | the phase(s) information read from a GPX file. Also prints values to the .lst file |
---|
2213 | ''' |
---|
2214 | |
---|
2215 | def PrintAtomsAndSig(General,Atoms,sigDict,sigKey): |
---|
2216 | pFile.write('\n Atoms:\n') |
---|
2217 | line = ' name x y z frac Uiso U11 U22 U33 U12 U13 U23' |
---|
2218 | if General['Type'] == 'macromolecular': |
---|
2219 | line = ' res no residue chain '+line |
---|
2220 | cx,ct,cs,cia = General['AtomPtrs'] |
---|
2221 | pFile.write(line+'\n') |
---|
2222 | pFile.write(135*'-'+'\n') |
---|
2223 | fmt = {0:'%7s',ct:'%7s',cx:'%10.5f',cx+1:'%10.5f',cx+2:'%10.5f',cx+3:'%8.3f',cia+1:'%8.5f', |
---|
2224 | cia+2:'%8.5f',cia+3:'%8.5f',cia+4:'%8.5f',cia+5:'%8.5f',cia+6:'%8.5f',cia+7:'%8.5f'} |
---|
2225 | #noFXsig = {cx:[10*' ','%10s'],cx+1:[10*' ','%10s'],cx+2:[10*' ','%10s'],cx+3:[8*' ','%8s']} |
---|
2226 | for atyp in General['AtomTypes']: #zero composition |
---|
2227 | General['NoAtoms'][atyp] = 0. |
---|
2228 | for i,at in enumerate(Atoms): |
---|
2229 | General['NoAtoms'][at[ct]] += at[cx+3]*float(at[cx+5]) #new composition |
---|
2230 | if General['Type'] == 'macromolecular': |
---|
2231 | name = ' %s %s %s %s:'%(at[0],at[1],at[2],at[3]) |
---|
2232 | valstr = ' values: ' |
---|
2233 | sigstr = ' sig : ' |
---|
2234 | else: |
---|
2235 | name = fmt[0]%(at[ct-1])+fmt[1]%(at[ct])+':' |
---|
2236 | valstr = ' values:' |
---|
2237 | sigstr = ' sig : ' |
---|
2238 | for ind in range(cx,cx+4): |
---|
2239 | sigind = str(i)+':'+str(ind) |
---|
2240 | valstr += fmt[ind]%(at[ind]) |
---|
2241 | # if sigind in atomsSig: |
---|
2242 | # sigstr += fmt[ind]%(atomsSig[sigind]) |
---|
2243 | # else: |
---|
2244 | # sigstr += noFXsig[ind][1]%(noFXsig[ind][0]) |
---|
2245 | sigstr += fmtESD(sigKey[sigind],sigDict,fmt[ind]) |
---|
2246 | if at[cia] == 'I': |
---|
2247 | valstr += fmt[cia+1]%(at[cia+1]) |
---|
2248 | # if '%d:%d'%(i,cia+1) in atomsSig: |
---|
2249 | # sigstr += fmt[cia+1]%(atomsSig['%d:%d'%(i,cia+1)]) |
---|
2250 | # else: |
---|
2251 | # sigstr += 8*' ' |
---|
2252 | sigstr += fmtESD(sigKey['%d:%d'%(i,cia+1)],sigDict,fmt[cia+1]) |
---|
2253 | else: |
---|
2254 | valstr += 8*' ' |
---|
2255 | sigstr += 8*' ' |
---|
2256 | for ind in range(cia+2,cia+8): |
---|
2257 | sigind = str(i)+':'+str(ind) |
---|
2258 | valstr += fmt[ind]%(at[ind]) |
---|
2259 | # if sigind in atomsSig: |
---|
2260 | # sigstr += fmt[ind]%(atomsSig[sigind]) |
---|
2261 | # else: |
---|
2262 | # sigstr += 8*' ' |
---|
2263 | sigstr += fmtESD(sigKey[sigind],sigDict,fmt[ind]) |
---|
2264 | pFile.write(name+'\n') |
---|
2265 | pFile.write(valstr+'\n') |
---|
2266 | pFile.write(sigstr+'\n') |
---|
2267 | |
---|
2268 | def PrintMomentsAndSig(General,Atoms,atomsSig): |
---|
2269 | cell = General['Cell'][1:7] |
---|
2270 | G = G2lat.fillgmat(cell) |
---|
2271 | ast = np.sqrt(np.diag(G)) |
---|
2272 | GS = G/np.outer(ast,ast) |
---|
2273 | pFile.write('\n Magnetic Moments:\n') #add magnitude & angle, etc.? TBD |
---|
2274 | line = ' name Mx My Mz |Mag|' |
---|
2275 | cx,ct,cs,cia = General['AtomPtrs'] |
---|
2276 | cmx = cx+4 |
---|
2277 | AtInfo = dict(zip(General['AtomTypes'],General['Lande g'])) |
---|
2278 | pFile.write(line+'\n') |
---|
2279 | pFile.write(135*'-'+'\n') |
---|
2280 | fmt = {0:'%7s',ct:'%7s',cmx:'%10.3f',cmx+1:'%10.3f',cmx+2:'%10.3f'} |
---|
2281 | noFXsig = {cmx:[10*' ','%10s'],cmx+1:[10*' ','%10s'],cmx+2:[10*' ','%10s']} |
---|
2282 | for i,at in enumerate(Atoms): |
---|
2283 | if AtInfo[at[ct]]: |
---|
2284 | name = fmt[0]%(at[ct-1])+fmt[1]%(at[ct])+':' |
---|
2285 | valstr = ' values:' |
---|
2286 | sigstr = ' sig :' |
---|
2287 | for ind in range(cmx,cmx+3): |
---|
2288 | sigind = str(i)+':'+str(ind) |
---|
2289 | valstr += fmt[ind]%(at[ind]) |
---|
2290 | if sigind in atomsSig: |
---|
2291 | sigstr += fmt[ind]%(atomsSig[sigind]) |
---|
2292 | else: |
---|
2293 | sigstr += noFXsig[ind][1]%(noFXsig[ind][0]) |
---|
2294 | mag = np.array(at[cmx:cmx+3]) |
---|
2295 | Mag = np.sqrt(np.inner(mag,np.inner(mag,GS))) |
---|
2296 | valstr += '%10.3f'%Mag |
---|
2297 | sigstr += 10*' ' |
---|
2298 | pFile.write(name+'\n') |
---|
2299 | pFile.write(valstr+'\n') |
---|
2300 | pFile.write(sigstr+'\n') |
---|
2301 | |
---|
2302 | def PrintWavesAndSig(General,Atoms,wavesSig): |
---|
2303 | cx,ct,cs,cia = General['AtomPtrs'] |
---|
2304 | pFile.write('\n Modulation waves\n') |
---|
2305 | names = {'Sfrac':['Fsin','Fcos','Fzero','Fwid'],'Spos':['Xsin','Ysin','Zsin','Xcos','Ycos','Zcos','Tmin','Tmax','Xmax','Ymax','Zmax'], |
---|
2306 | 'Sadp':['U11sin','U22sin','U33sin','U12sin','U13sin','U23sin','U11cos','U22cos', |
---|
2307 | 'U33cos','U12cos','U13cos','U23cos'],'Smag':['MXsin','MYsin','MZsin','MXcos','MYcos','MZcos']} |
---|
2308 | pFile.write(135*'-'+'\n') |
---|
2309 | for i,at in enumerate(Atoms): |
---|
2310 | AtomSS = at[-1]['SS1'] |
---|
2311 | for Stype in ['Sfrac','Spos','Sadp','Smag']: |
---|
2312 | Waves = AtomSS[Stype] |
---|
2313 | if len(Waves) > 1: |
---|
2314 | waveType = Waves[0] |
---|
2315 | else: |
---|
2316 | continue |
---|
2317 | if len(Waves): |
---|
2318 | pFile.write(' atom: %s, site sym: %s, type: %s wave type: %s:\n'% |
---|
2319 | (at[ct-1],at[cs],Stype,waveType)) |
---|
2320 | for iw,wave in enumerate(Waves[1:]): |
---|
2321 | stiw = ':'+str(i)+':'+str(iw) |
---|
2322 | namstr = ' names :' |
---|
2323 | valstr = ' values:' |
---|
2324 | sigstr = ' esds :' |
---|
2325 | if Stype == 'Spos': |
---|
2326 | nt = 6 |
---|
2327 | ot = 0 |
---|
2328 | if waveType in ['ZigZag','Block',] and not iw: |
---|
2329 | nt = 5 |
---|
2330 | ot = 6 |
---|
2331 | for j in range(nt): |
---|
2332 | name = names['Spos'][j+ot] |
---|
2333 | namstr += '%12s'%(name) |
---|
2334 | valstr += '%12.4f'%(wave[0][j]) |
---|
2335 | if name+stiw in wavesSig: |
---|
2336 | sigstr += '%12.4f'%(wavesSig[name+stiw]) |
---|
2337 | else: |
---|
2338 | sigstr += 12*' ' |
---|
2339 | elif Stype == 'Sfrac': |
---|
2340 | ot = 0 |
---|
2341 | if 'Crenel' in waveType and not iw: |
---|
2342 | ot = 2 |
---|
2343 | for j in range(2): |
---|
2344 | name = names['Sfrac'][j+ot] |
---|
2345 | namstr += '%12s'%(names['Sfrac'][j+ot]) |
---|
2346 | valstr += '%12.4f'%(wave[0][j]) |
---|
2347 | if name+stiw in wavesSig: |
---|
2348 | sigstr += '%12.4f'%(wavesSig[name+stiw]) |
---|
2349 | else: |
---|
2350 | sigstr += 12*' ' |
---|
2351 | elif Stype == 'Sadp': |
---|
2352 | for j in range(12): |
---|
2353 | name = names['Sadp'][j] |
---|
2354 | namstr += '%10s'%(names['Sadp'][j]) |
---|
2355 | valstr += '%10.6f'%(wave[0][j]) |
---|
2356 | if name+stiw in wavesSig: |
---|
2357 | sigstr += '%10.6f'%(wavesSig[name+stiw]) |
---|
2358 | else: |
---|
2359 | sigstr += 10*' ' |
---|
2360 | elif Stype == 'Smag': |
---|
2361 | for j in range(6): |
---|
2362 | name = names['Smag'][j] |
---|
2363 | namstr += '%12s'%(names['Smag'][j]) |
---|
2364 | valstr += '%12.4f'%(wave[0][j]) |
---|
2365 | if name+stiw in wavesSig: |
---|
2366 | sigstr += '%12.4f'%(wavesSig[name+stiw]) |
---|
2367 | else: |
---|
2368 | sigstr += 12*' ' |
---|
2369 | |
---|
2370 | pFile.write(namstr+'\n') |
---|
2371 | pFile.write(valstr+'\n') |
---|
2372 | pFile.write(sigstr+'\n') |
---|
2373 | |
---|
2374 | |
---|
2375 | def PrintRBObjPOAndSig(rbfx,rbsx): |
---|
2376 | for i in WriteRBObjPOAndSig(pfx,rbfx,rbsx,parmDict,sigDict): |
---|
2377 | pFile.write(i+'\n') |
---|
2378 | |
---|
2379 | def PrintRBObjTLSAndSig(rbfx,rbsx,TLS): |
---|
2380 | for i in WriteRBObjTLSAndSig(pfx,rbfx,rbsx,TLS,parmDict,sigDict): |
---|
2381 | pFile.write(i) |
---|
2382 | |
---|
2383 | def PrintRBObjSHCAndSig(rbfx,SHC,rbsx): |
---|
2384 | for i in WriteRBObjSHCAndSig(pfx,rbfx,rbsx,parmDict,sigDict,SHC): |
---|
2385 | pFile.write(i) |
---|
2386 | |
---|
2387 | def PrintRBObjTorAndSig(rbsx): |
---|
2388 | nTors = len(RBObj['Torsions']) |
---|
2389 | if nTors: |
---|
2390 | for i in WriteRBObjTorAndSig(pfx,rbsx,parmDict,sigDict,nTors): |
---|
2391 | pFile.write(i) |
---|
2392 | |
---|
2393 | def PrintSHtextureAndSig(textureData,SHtextureSig): |
---|
2394 | Tindx = 1.0 |
---|
2395 | Tvar = 0.0 |
---|
2396 | pFile.write('\n Spherical harmonics texture: Order: %d\n'%textureData['Order']) |
---|
2397 | names = ['omega','chi','phi'] |
---|
2398 | namstr = ' names :' |
---|
2399 | ptstr = ' values:' |
---|
2400 | sigstr = ' esds :' |
---|
2401 | for name in names: |
---|
2402 | namstr += '%12s'%(name) |
---|
2403 | ptstr += '%12.3f'%(textureData['Sample '+name][1]) |
---|
2404 | if 'Sample '+name in SHtextureSig: |
---|
2405 | sigstr += '%12.3f'%(SHtextureSig['Sample '+name]) |
---|
2406 | else: |
---|
2407 | sigstr += 12*' ' |
---|
2408 | pFile.write(namstr+'\n') |
---|
2409 | pFile.write(ptstr+'\n') |
---|
2410 | pFile.write(sigstr+'\n') |
---|
2411 | pFile.write('\n Texture coefficients:\n') |
---|
2412 | SHcoeff = textureData['SH Coeff'][1] |
---|
2413 | SHkeys = list(SHcoeff.keys()) |
---|
2414 | nCoeff = len(SHcoeff) |
---|
2415 | nBlock = nCoeff//10+1 |
---|
2416 | iBeg = 0 |
---|
2417 | iFin = min(iBeg+10,nCoeff) |
---|
2418 | for block in range(nBlock): |
---|
2419 | namstr = ' names :' |
---|
2420 | ptstr = ' values:' |
---|
2421 | sigstr = ' esds :' |
---|
2422 | for name in SHkeys[iBeg:iFin]: |
---|
2423 | namstr += '%12s'%(name) |
---|
2424 | ptstr += '%12.3f'%(SHcoeff[name]) |
---|
2425 | l = 2.0*eval(name.strip('C'))[0]+1 |
---|
2426 | Tindx += SHcoeff[name]**2/l |
---|
2427 | if name in SHtextureSig: |
---|
2428 | Tvar += (2.*SHcoeff[name]*SHtextureSig[name]/l)**2 |
---|
2429 | sigstr += '%12.3f'%(SHtextureSig[name]) |
---|
2430 | else: |
---|
2431 | sigstr += 12*' ' |
---|
2432 | pFile.write(namstr+'\n') |
---|
2433 | pFile.write(ptstr+'\n') |
---|
2434 | pFile.write(sigstr+'\n') |
---|
2435 | iBeg += 10 |
---|
2436 | iFin = min(iBeg+10,nCoeff) |
---|
2437 | pFile.write(' Texture index J = %.3f(%d)'%(Tindx,int(1000*np.sqrt(Tvar)))) |
---|
2438 | |
---|
2439 | ########################################################################## |
---|
2440 | # SetPhaseData starts here |
---|
2441 | if pFile: pFile.write('\n Phases:\n') |
---|
2442 | for phase in Phases: |
---|
2443 | if pFile: pFile.write(' Result for phase: %s\n'%phase) |
---|
2444 | if pFile: pFile.write(135*'='+'\n') |
---|
2445 | Phase = Phases[phase] |
---|
2446 | General = Phase['General'] |
---|
2447 | SGData = General['SGData'] |
---|
2448 | Atoms = Phase['Atoms'] |
---|
2449 | AtLookup = [] |
---|
2450 | if Atoms and not General.get('doPawley'): |
---|
2451 | cx,ct,cs,cia = General['AtomPtrs'] |
---|
2452 | AtLookup = G2mth.FillAtomLookUp(Atoms,cia+8) |
---|
2453 | cell = General['Cell'] |
---|
2454 | pId = Phase['pId'] |
---|
2455 | pfx = str(pId)+'::' |
---|
2456 | if cell[0]: |
---|
2457 | A,sigA = cellFill(pfx,SGData,parmDict,sigDict) |
---|
2458 | cellSig = getCellEsd(pfx,SGData,A,covData,unique=True) #includes sigVol |
---|
2459 | if pFile: pFile.write(' Reciprocal metric tensor: \n') |
---|
2460 | ptfmt = "%15.9f" |
---|
2461 | names = ['A11','A22','A33','A12','A13','A23'] |
---|
2462 | namstr = ' names :' |
---|
2463 | ptstr = ' values:' |
---|
2464 | sigstr = ' esds :' |
---|
2465 | for name,a,siga in zip(names,A,sigA): |
---|
2466 | namstr += '%15s'%(name) |
---|
2467 | ptstr += ptfmt%(a) |
---|
2468 | if siga: |
---|
2469 | sigstr += ptfmt%(siga) |
---|
2470 | else: |
---|
2471 | sigstr += 15*' ' |
---|
2472 | if pFile: pFile.write(namstr+'\n') |
---|
2473 | if pFile: pFile.write(ptstr+'\n') |
---|
2474 | if pFile: pFile.write(sigstr+'\n') |
---|
2475 | cell[1:7] = G2lat.A2cell(A) |
---|
2476 | cell[7] = G2lat.calc_V(A) |
---|
2477 | if pFile: pFile.write(' New unit cell:\n') |
---|
2478 | ptfmt = ["%12.6f","%12.6f","%12.6f","%12.4f","%12.4f","%12.4f","%12.3f"] |
---|
2479 | names = ['a','b','c','alpha','beta','gamma','Volume'] |
---|
2480 | namstr = ' names :' |
---|
2481 | ptstr = ' values:' |
---|
2482 | sigstr = ' esds :' |
---|
2483 | for name,fmt,a,siga in zip(names,ptfmt,cell[1:8],cellSig): |
---|
2484 | namstr += '%12s'%(name) |
---|
2485 | ptstr += fmt%(a) |
---|
2486 | if siga and siga > 0: |
---|
2487 | sigstr += fmt%(siga) |
---|
2488 | else: |
---|
2489 | sigstr += 12*' ' |
---|
2490 | if pFile: pFile.write(namstr+'\n') |
---|
2491 | if pFile: pFile.write(ptstr+'\n') |
---|
2492 | if pFile: pFile.write(sigstr+'\n') |
---|
2493 | ik = 6 #for Pawley stuff below |
---|
2494 | if General.get('Modulated',False): |
---|
2495 | ik = 7 |
---|
2496 | Vec,vRef,maxH = General['SuperVec'] |
---|
2497 | if vRef: |
---|
2498 | if pFile: pFile.write(' New modulation vector:\n') |
---|
2499 | namstr = ' names :' |
---|
2500 | ptstr = ' values:' |
---|
2501 | sigstr = ' esds :' |
---|
2502 | for iv,var in enumerate(['mV0','mV1','mV2']): |
---|
2503 | namstr += '%12s'%(pfx+var) |
---|
2504 | ptstr += '%12.6f'%(parmDict[pfx+var]) |
---|
2505 | if pfx+var in sigDict: |
---|
2506 | Vec[iv] = parmDict[pfx+var] |
---|
2507 | sigstr += '%12.6f'%(sigDict[pfx+var]) |
---|
2508 | else: |
---|
2509 | sigstr += 12*' ' |
---|
2510 | if pFile: pFile.write(namstr+'\n') |
---|
2511 | if pFile: pFile.write(ptstr+'\n') |
---|
2512 | if pFile: pFile.write(sigstr+'\n') |
---|
2513 | |
---|
2514 | General['Mass'] = 0. |
---|
2515 | if Phase['General'].get('doPawley'): |
---|
2516 | pawleyRef = Phase['Pawley ref'] |
---|
2517 | for i,refl in enumerate(pawleyRef): |
---|
2518 | key = pfx+'PWLref:'+str(i) |
---|
2519 | refl[ik] = parmDict[key] |
---|
2520 | if key in sigDict: |
---|
2521 | refl[ik+1] = sigDict[key] |
---|
2522 | else: |
---|
2523 | refl[ik+1] = 0 |
---|
2524 | else: |
---|
2525 | VRBIds = RBIds['Vector'] |
---|
2526 | RRBIds = RBIds['Residue'] |
---|
2527 | SRBIds = RBIds['Spin'] |
---|
2528 | RBModels = Phase['RBModels'] |
---|
2529 | if pFile: |
---|
2530 | for irb,RBObj in enumerate(RBModels.get('Vector',[])): |
---|
2531 | jrb = VRBIds.index(RBObj['RBId']) |
---|
2532 | rbsx = str(irb)+':'+str(jrb) |
---|
2533 | pFile.write(' Vector rigid body parameters:\n') |
---|
2534 | PrintRBObjPOAndSig('RBV',rbsx) |
---|
2535 | PrintRBObjTLSAndSig('RBV',rbsx,RBObj['ThermalMotion'][0]) |
---|
2536 | for irb,RBObj in enumerate(RBModels.get('Residue',[])): |
---|
2537 | jrb = RRBIds.index(RBObj['RBId']) |
---|
2538 | rbsx = str(irb)+':'+str(jrb) |
---|
2539 | pFile.write(' Residue rigid body parameters:\n') |
---|
2540 | PrintRBObjPOAndSig('RBR',rbsx) |
---|
2541 | PrintRBObjTLSAndSig('RBR',rbsx,RBObj['ThermalMotion'][0]) |
---|
2542 | PrintRBObjTorAndSig(rbsx) |
---|
2543 | for irb,RBObj in enumerate(RBModels.get('Spin',[])): |
---|
2544 | for ish in range(len(RBObj['RBId'])): #shell no. |
---|
2545 | jrb = SRBIds.index(RBObj['RBId'][ish]) #spin rb no. |
---|
2546 | rbsx = str(irb)+':'+str(jrb) |
---|
2547 | pFile.write(' Spin rigid body parameters:\n') |
---|
2548 | PrintRBObjPOAndSig('RBS',rbsx) |
---|
2549 | rbsx = ':%d:%d'%(irb,jrb) |
---|
2550 | PrintRBObjSHCAndSig('RBSSh;%d;'%ish,RBObj['SHC'][ish],rbsx) |
---|
2551 | atomsSig = {} |
---|
2552 | sigKey = {} |
---|
2553 | wavesSig = {} |
---|
2554 | cx,ct,cs,cia = General['AtomPtrs'] |
---|
2555 | for i,at in enumerate(Atoms): |
---|
2556 | names = {cx:pfx+'Ax:'+str(i),cx+1:pfx+'Ay:'+str(i),cx+2:pfx+'Az:'+str(i),cx+3:pfx+'Afrac:'+str(i), |
---|
2557 | cia+1:pfx+'AUiso:'+str(i),cia+2:pfx+'AU11:'+str(i),cia+3:pfx+'AU22:'+str(i),cia+4:pfx+'AU33:'+str(i), |
---|
2558 | cia+5:pfx+'AU12:'+str(i),cia+6:pfx+'AU13:'+str(i),cia+7:pfx+'AU23:'+str(i), |
---|
2559 | cx+4:pfx+'AMx:'+str(i),cx+5:pfx+'AMy:'+str(i),cx+6:pfx+'AMz:'+str(i)} |
---|
2560 | for ind in range(cx,cx+4): |
---|
2561 | at[ind] = parmDict[names[ind]] |
---|
2562 | if ind in range(cx,cx+3): |
---|
2563 | name = names[ind].replace('A','dA') |
---|
2564 | else: |
---|
2565 | name = names[ind] |
---|
2566 | sigKey[str(i)+':'+str(ind)] = name |
---|
2567 | if name in sigDict: |
---|
2568 | atomsSig[str(i)+':'+str(ind)] = sigDict[name] |
---|
2569 | if at[cia] == 'I': |
---|
2570 | at[cia+1] = parmDict[names[cia+1]] |
---|
2571 | sigKey['%d:%d'%(i,cia+1)] = names[cia+1] |
---|
2572 | if names[cia+1] in sigDict: |
---|
2573 | atomsSig['%d:%d'%(i,cia+1)] = sigDict[names[cia+1]] |
---|
2574 | else: |
---|
2575 | for ind in range(cia+2,cia+8): |
---|
2576 | at[ind] = parmDict[names[ind]] |
---|
2577 | sigKey[str(i)+':'+str(ind)] = names[ind] |
---|
2578 | if names[ind] in sigDict: |
---|
2579 | atomsSig[str(i)+':'+str(ind)] = sigDict[names[ind]] |
---|
2580 | if General['Type'] == 'magnetic': |
---|
2581 | for ind in range(cx+4,cx+7): |
---|
2582 | at[ind] = parmDict[names[ind]] |
---|
2583 | sigKey[str(i)+':'+str(ind)] = names[ind] |
---|
2584 | if names[ind] in sigDict: |
---|
2585 | atomsSig[str(i)+':'+str(ind)] = sigDict[names[ind]] |
---|
2586 | ind = General['AtomTypes'].index(at[ct]) |
---|
2587 | General['Mass'] += General['AtomMass'][ind]*at[cx+3]*at[cx+5] |
---|
2588 | if General.get('Modulated',False): |
---|
2589 | AtomSS = at[-1]['SS1'] |
---|
2590 | for Stype in ['Sfrac','Spos','Sadp','Smag']: |
---|
2591 | Waves = AtomSS[Stype] |
---|
2592 | if len(Waves): |
---|
2593 | waveType = Waves[0] |
---|
2594 | else: |
---|
2595 | continue |
---|
2596 | for iw,wave in enumerate(Waves[1:]): |
---|
2597 | stiw = str(i)+':'+str(iw) |
---|
2598 | if Stype == 'Spos': |
---|
2599 | if waveType in ['ZigZag','Block',] and not iw: |
---|
2600 | names = ['Tmin:'+stiw,'Tmax:'+stiw,'Xmax:'+stiw,'Ymax:'+stiw,'Zmax:'+stiw] |
---|
2601 | else: |
---|
2602 | names = ['Xsin:'+stiw,'Ysin:'+stiw,'Zsin:'+stiw, |
---|
2603 | 'Xcos:'+stiw,'Ycos:'+stiw,'Zcos:'+stiw] |
---|
2604 | elif Stype == 'Sadp': |
---|
2605 | names = ['U11sin:'+stiw,'U22sin:'+stiw,'U33sin:'+stiw, |
---|
2606 | 'U12sin:'+stiw,'U13sin:'+stiw,'U23sin:'+stiw, |
---|
2607 | 'U11cos:'+stiw,'U22cos:'+stiw,'U33cos:'+stiw, |
---|
2608 | 'U12cos:'+stiw,'U13cos:'+stiw,'U23cos:'+stiw] |
---|
2609 | elif Stype == 'Sfrac': |
---|
2610 | if 'Crenel' in waveType and not iw: |
---|
2611 | names = ['Fzero:'+stiw,'Fwid:'+stiw] |
---|
2612 | else: |
---|
2613 | names = ['Fsin:'+stiw,'Fcos:'+stiw] |
---|
2614 | elif Stype == 'Smag': |
---|
2615 | names = ['MXsin:'+stiw,'MYsin:'+stiw,'MZsin:'+stiw, |
---|
2616 | 'MXcos:'+stiw,'MYcos:'+stiw,'MZcos:'+stiw] |
---|
2617 | for iname,name in enumerate(names): |
---|
2618 | sigKey[name] = pfx+name |
---|
2619 | AtomSS[Stype][iw+1][0][iname] = parmDict[pfx+name] |
---|
2620 | if pfx+name in sigDict: |
---|
2621 | wavesSig[name] = sigDict[pfx+name] |
---|
2622 | |
---|
2623 | if pFile: PrintAtomsAndSig(General,Atoms,sigDict,sigKey) |
---|
2624 | if pFile and General['Type'] == 'magnetic': |
---|
2625 | PrintMomentsAndSig(General,Atoms,atomsSig) |
---|
2626 | if pFile and General.get('Modulated',False): |
---|
2627 | PrintWavesAndSig(General,Atoms,wavesSig) |
---|
2628 | |
---|
2629 | density = G2mth.getDensity(General)[0] |
---|
2630 | if pFile: pFile.write('\n Density: {:.4f} g/cm**3\n'.format(density)) |
---|
2631 | |
---|
2632 | |
---|
2633 | textureData = General['SH Texture'] |
---|
2634 | if textureData['Order']: |
---|
2635 | SHtextureSig = {} |
---|
2636 | for name in ['omega','chi','phi']: |
---|
2637 | aname = pfx+'SH '+name |
---|
2638 | textureData['Sample '+name][1] = parmDict[aname] |
---|
2639 | if aname in sigDict: |
---|
2640 | SHtextureSig['Sample '+name] = sigDict[aname] |
---|
2641 | for name in textureData['SH Coeff'][1]: |
---|
2642 | aname = pfx+name |
---|
2643 | textureData['SH Coeff'][1][name] = parmDict[aname] |
---|
2644 | if aname in sigDict: |
---|
2645 | SHtextureSig[name] = sigDict[aname] |
---|
2646 | PrintSHtextureAndSig(textureData,SHtextureSig) |
---|
2647 | if phase in RestraintDict and not Phase['General'].get('doPawley'): |
---|
2648 | PrintRestraints(cell[1:7],SGData,General['AtomPtrs'],Atoms,AtLookup, |
---|
2649 | textureData,RestraintDict[phase],pFile) |
---|
2650 | |
---|
2651 | def SetISOmodes(parmDict,sigDict,Phases,pFile=None): |
---|
2652 | '''After a refinement, sets the values for the ISODISTORT modes into |
---|
2653 | the parameter and s.u. dicts. |
---|
2654 | Also, in the case of a non-sequential refinement, prints them into |
---|
2655 | the project's .lst file. |
---|
2656 | |
---|
2657 | :param dict parmDict: parameter dict |
---|
2658 | :param dict sigDict: s.u. (uncertainty) dict |
---|
2659 | :param dict Phases: Phase info from tree/.gpx |
---|
2660 | :param file pFile: file for .lst info or None (None for sequential fits) |
---|
2661 | ''' |
---|
2662 | for phase in Phases: |
---|
2663 | if 'ISODISTORT' not in Phases[phase]: continue |
---|
2664 | data = Phases[phase] |
---|
2665 | ISO = data['ISODISTORT'] |
---|
2666 | atNames = [atom[0] for atom in data['Atoms']] |
---|
2667 | |
---|
2668 | if 'G2VarList' in ISO: |
---|
2669 | deltaList = [] |
---|
2670 | notfound = [] |
---|
2671 | for gv,Ilbl in zip(ISO['G2VarList'],ISO['IsoVarList']): |
---|
2672 | dvar = gv.varname() |
---|
2673 | var = dvar.replace('::dA','::A') |
---|
2674 | atnum = atNames.index(Ilbl[:Ilbl.rfind('_')]) |
---|
2675 | v = Ilbl[Ilbl.rfind('_')+1:] |
---|
2676 | pval = ISO['G2parentCoords'][atnum][['dx','dy','dz'].index(v)] |
---|
2677 | if var in parmDict: |
---|
2678 | cval = parmDict[var] |
---|
2679 | else: |
---|
2680 | notfound.append(var) |
---|
2681 | continue |
---|
2682 | deltaList.append(cval-pval) |
---|
2683 | |
---|
2684 | if notfound and pFile: |
---|
2685 | msg = 'SetISOmodes warning: Atom parameters ' |
---|
2686 | for i,v in enumerate(notfound): |
---|
2687 | if i == len(notfound)-1: |
---|
2688 | msg += ' & ' |
---|
2689 | elif i != 0: |
---|
2690 | msg += ', ' |
---|
2691 | msg += v |
---|
2692 | print(msg,'not found') |
---|
2693 | print(' skipping computation for modes:') |
---|
2694 | for i,j in zip(ISO['IsoModeList'],ISO['G2ModeList']): |
---|
2695 | print(' ',i,'({})'.format(j)) |
---|
2696 | continue |
---|
2697 | elif notfound: |
---|
2698 | continue |
---|
2699 | modeVals = np.inner(ISO['Var2ModeMatrix'],deltaList) |
---|
2700 | |
---|
2701 | if pFile: |
---|
2702 | pFile.write('\n ISODISTORT Displacive Modes for phase {}\n'.format( |
---|
2703 | data['General'].get('Name',''))) |
---|
2704 | |
---|
2705 | l = str(max([len(i) for i in ISO['IsoModeList']])+3) |
---|
2706 | fmt = ' {:'+l+'}{}' |
---|
2707 | for varid,[var,val,norm,G2mode] in enumerate(zip( |
---|
2708 | ISO['IsoModeList'],modeVals,ISO['NormList'],ISO['G2ModeList'] )): |
---|
2709 | try: |
---|
2710 | value = G2mth.ValEsd(val/norm,-0.001) |
---|
2711 | item = str(G2mode).replace('::','::nv-') |
---|
2712 | parmDict[str(G2mode)] = val/norm |
---|
2713 | if item in sigDict: |
---|
2714 | ISO['modeDispl'][varid] = val/norm |
---|
2715 | value = G2mth.ValEsd(val/norm,sigDict[item]/norm) |
---|
2716 | sigDict[str(G2mode)] = sigDict[item]/norm |
---|
2717 | except TypeError: |
---|
2718 | value = '?' |
---|
2719 | if pFile: |
---|
2720 | pFile.write(fmt.format(var,value)+'\n') |
---|
2721 | |
---|
2722 | if 'G2OccVarList' in ISO: #untested - probably wrong |
---|
2723 | deltaOccList = [] |
---|
2724 | notfound = [] |
---|
2725 | for gv,Ilbl in zip(ISO['G2OccVarList'],ISO['OccVarList']): |
---|
2726 | var = gv.varname() |
---|
2727 | albl = Ilbl[:Ilbl.rfind('_')] |
---|
2728 | pval = ISO['BaseOcc'][albl] |
---|
2729 | if var in parmDict: |
---|
2730 | cval = parmDict[var] |
---|
2731 | else: |
---|
2732 | notfound.append(var) |
---|
2733 | continue |
---|
2734 | deltaOccList.append(cval-pval) |
---|
2735 | |
---|
2736 | if notfound and pFile: |
---|
2737 | msg = 'SetISOmodes warning: Atom parameters ' |
---|
2738 | for i,v in enumerate(notfound): |
---|
2739 | if i == len(notfound)-1: |
---|
2740 | msg += ' & ' |
---|
2741 | elif i != 0: |
---|
2742 | msg += ', ' |
---|
2743 | msg += v |
---|
2744 | print(msg,'not found') |
---|
2745 | print(' skipping computation for modes:') |
---|
2746 | for i,j in zip(ISO['IsoModeList'],ISO['G2ModeList']): |
---|
2747 | print(' ',i,'({})'.format(j)) |
---|
2748 | continue |
---|
2749 | elif notfound: |
---|
2750 | continue |
---|
2751 | modeOccVals = np.inner(ISO['Var2OccMatrix'],deltaOccList) |
---|
2752 | |
---|
2753 | if pFile: |
---|
2754 | pFile.write('\n ISODISTORT Occupancy Modes for phase {}\n'.format(data['General'].get('Name',''))) |
---|
2755 | l = str(max([len(i) for i in ISO['OccModeList']])+3) |
---|
2756 | fmt = ' {:'+l+'}{}' |
---|
2757 | for var,val,norm,G2mode in zip( |
---|
2758 | ISO['OccModeList'],modeOccVals,ISO['OccNormList'],ISO['G2OccModeList'] ): |
---|
2759 | try: |
---|
2760 | value = G2py3.FormatSigFigs(val/norm) |
---|
2761 | if str(G2mode) in sigDict: |
---|
2762 | value = G2mth.ValEsd(val/norm,sigDict[str(G2mode)]/norm) |
---|
2763 | except TypeError: |
---|
2764 | value = '?' |
---|
2765 | if pFile: |
---|
2766 | pFile.write(fmt.format(var,value)+'\n') |
---|
2767 | |
---|
2768 | ################################################################################ |
---|
2769 | ##### Histogram & Phase data |
---|
2770 | ################################################################################ |
---|
2771 | |
---|
2772 | def GetHistogramPhaseData(Phases,Histograms,Controls={},Print=True,pFile=None,resetRefList=True): |
---|
2773 | '''Loads the HAP histogram/phase information into dicts |
---|
2774 | |
---|
2775 | :param dict Phases: phase information |
---|
2776 | :param dict Histograms: Histogram information |
---|
2777 | :param bool Print: prints information as it is read |
---|
2778 | :param file pFile: file object to print to (the default, None causes printing to the console) |
---|
2779 | :param bool resetRefList: Should the contents of the Reflection List be initialized |
---|
2780 | on loading. The default, True, initializes the Reflection List as it is loaded. |
---|
2781 | |
---|
2782 | :returns: (hapVary,hapDict,controlDict) |
---|
2783 | * hapVary: list of refined variables |
---|
2784 | * hapDict: dict with refined variables and their values |
---|
2785 | * controlDict: dict with fixed parameters |
---|
2786 | ''' |
---|
2787 | |
---|
2788 | def PrintSize(hapData): |
---|
2789 | if hapData[0] in ['isotropic','uniaxial']: |
---|
2790 | line = '\n Size model : %9s'%(hapData[0]) |
---|
2791 | line += ' equatorial:'+'%12.3f'%(hapData[1][0])+' Refine? '+str(hapData[2][0]) |
---|
2792 | if hapData[0] == 'uniaxial': |
---|
2793 | line += ' axial:'+'%12.3f'%(hapData[1][1])+' Refine? '+str(hapData[2][1]) |
---|
2794 | line += '\n\t LG mixing coeff.: %12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2]) |
---|
2795 | pFile.write(line+'\n') |
---|
2796 | else: |
---|
2797 | pFile.write('\n Size model : %s\n\t LG mixing coeff.:%12.4f Refine? %s\n'% |
---|
2798 | (hapData[0],hapData[1][2],hapData[2][2])) |
---|
2799 | Snames = ['S11','S22','S33','S12','S13','S23'] |
---|
2800 | ptlbls = ' names :' |
---|
2801 | ptstr = ' values:' |
---|
2802 | varstr = ' refine:' |
---|
2803 | for i,name in enumerate(Snames): |
---|
2804 | ptlbls += '%12s' % (name) |
---|
2805 | ptstr += '%12.3f' % (hapData[4][i]) |
---|
2806 | varstr += '%12s' % (str(hapData[5][i])) |
---|
2807 | pFile.write(ptlbls+'\n') |
---|
2808 | pFile.write(ptstr+'\n') |
---|
2809 | pFile.write(varstr+'\n') |
---|
2810 | |
---|
2811 | def PrintMuStrain(hapData,SGData): |
---|
2812 | if hapData[0] in ['isotropic','uniaxial']: |
---|
2813 | line = '\n Mustrain model: %9s'%(hapData[0]) |
---|
2814 | line += ' equatorial:'+'%12.1f'%(hapData[1][0])+' Refine? '+str(hapData[2][0]) |
---|
2815 | if hapData[0] == 'uniaxial': |
---|
2816 | line += ' axial:'+'%12.1f'%(hapData[1][1])+' Refine? '+str(hapData[2][1]) |
---|
2817 | line +='\n\t LG mixing coeff.:%12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2]) |
---|
2818 | pFile.write(line+'\n') |
---|
2819 | else: |
---|
2820 | pFile.write('\n Mustrain model: %s\n\t LG mixing coeff.:%12.4f Refine? %s\n'% |
---|
2821 | (hapData[0],hapData[1][2],hapData[2][2])) |
---|
2822 | Snames = G2spc.MustrainNames(SGData) |
---|
2823 | ptlbls = ' names :' |
---|
2824 | ptstr = ' values:' |
---|
2825 | varstr = ' refine:' |
---|
2826 | for i,name in enumerate(Snames): |
---|
2827 | ptlbls += '%12s' % (name) |
---|
2828 | ptstr += '%12.1f' % (hapData[4][i]) |
---|
2829 | varstr += '%12s' % (str(hapData[5][i])) |
---|
2830 | pFile.write(ptlbls+'\n') |
---|
2831 | pFile.write(ptstr+'\n') |
---|
2832 | pFile.write(varstr+'\n') |
---|
2833 | |
---|
2834 | def PrintHStrain(hapData,SGData): |
---|
2835 | pFile.write('\n Hydrostatic/elastic strain:\n') |
---|
2836 | Hsnames = G2spc.HStrainNames(SGData) |
---|
2837 | ptlbls = ' names :' |
---|
2838 | ptstr = ' values:' |
---|
2839 | varstr = ' refine:' |
---|
2840 | for i,name in enumerate(Hsnames): |
---|
2841 | ptlbls += '%12s' % (name) |
---|
2842 | ptstr += '%12.4g' % (hapData[0][i]) |
---|
2843 | varstr += '%12s' % (str(hapData[1][i])) |
---|
2844 | pFile.write(ptlbls+'\n') |
---|
2845 | pFile.write(ptstr+'\n') |
---|
2846 | pFile.write(varstr+'\n') |
---|
2847 | |
---|
2848 | def PrintSHPO(hapData): |
---|
2849 | pFile.write('\n Spherical harmonics preferred orientation: Order: %d Refine? %s\n'%(hapData[4],hapData[2])) |
---|
2850 | ptlbls = ' names :' |
---|
2851 | ptstr = ' values:' |
---|
2852 | for item in hapData[5]: |
---|
2853 | ptlbls += '%12s'%(item) |
---|
2854 | ptstr += '%12.3f'%(hapData[5][item]) |
---|
2855 | pFile.write(ptlbls+'\n') |
---|
2856 | pFile.write(ptstr+'\n') |
---|
2857 | |
---|
2858 | def PrintBabinet(hapData): |
---|
2859 | pFile.write('\n Babinet form factor modification:\n') |
---|
2860 | ptlbls = ' names :' |
---|
2861 | ptstr = ' values:' |
---|
2862 | varstr = ' refine:' |
---|
2863 | for name in ['BabA','BabU']: |
---|
2864 | ptlbls += '%12s' % (name) |
---|
2865 | ptstr += '%12.6f' % (hapData[name][0]) |
---|
2866 | varstr += '%12s' % (str(hapData[name][1])) |
---|
2867 | pFile.write(ptlbls+'\n') |
---|
2868 | pFile.write(ptstr+'\n') |
---|
2869 | pFile.write(varstr+'\n') |
---|
2870 | |
---|
2871 | hapDict = {} |
---|
2872 | hapVary = [] |
---|
2873 | controlDict = {} |
---|
2874 | |
---|
2875 | for phase in Phases: |
---|
2876 | HistoPhase = Phases[phase]['Histograms'] |
---|
2877 | SGData = Phases[phase]['General']['SGData'] |
---|
2878 | cell = Phases[phase]['General']['Cell'][1:7] |
---|
2879 | A = G2lat.cell2A(cell) |
---|
2880 | if Phases[phase]['General'].get('Modulated',False): |
---|
2881 | SSGData = Phases[phase]['General']['SSGData'] |
---|
2882 | Vec,x,maxH = Phases[phase]['General']['SuperVec'] |
---|
2883 | pId = Phases[phase]['pId'] |
---|
2884 | for histogram in Histograms: |
---|
2885 | if histogram not in HistoPhase and phase in Histograms[histogram]['Reflection Lists']: |
---|
2886 | #remove previously created reflection list if histogram is removed from phase |
---|
2887 | #print("removing ",phase,"from",histogram) |
---|
2888 | del Histograms[histogram]['Reflection Lists'][phase] |
---|
2889 | histoList = list(HistoPhase.keys()) |
---|
2890 | histoList.sort() |
---|
2891 | for histogram in histoList: |
---|
2892 | try: |
---|
2893 | Histogram = Histograms[histogram] |
---|
2894 | except KeyError: |
---|
2895 | #skip if histogram not included e.g. in a sequential refinement |
---|
2896 | continue |
---|
2897 | if not HistoPhase[histogram]['Use']: #remove previously created & now unused phase reflection list |
---|
2898 | if phase in Histograms[histogram]['Reflection Lists']: |
---|
2899 | del Histograms[histogram]['Reflection Lists'][phase] |
---|
2900 | continue |
---|
2901 | hapData = HistoPhase[histogram] |
---|
2902 | hId = Histogram['hId'] |
---|
2903 | if 'PWDR' in histogram: |
---|
2904 | limits = Histogram['Limits'][1] |
---|
2905 | inst = Histogram['Instrument Parameters'][0] #TODO - grab table here if present |
---|
2906 | if 'C' in inst['Type'][1]: |
---|
2907 | try: |
---|
2908 | wave = inst['Lam'][1] |
---|
2909 | except KeyError: |
---|
2910 | wave = inst['Lam1'][1] |
---|
2911 | dmin = wave/(2.0*sind(limits[1]/2.0)) |
---|
2912 | elif 'T' in inst['Type'][0]: |
---|
2913 | dmin = G2lat.Pos2dsp(inst,limits[1]) |
---|
2914 | dmin = limits[0]/inst['difC'][1] |
---|
2915 | elif 'E' in inst['Type'][0]: |
---|
2916 | dmin = G2lat.Pos2dsp(inst,limits[1]) |
---|
2917 | else: |
---|
2918 | wave = inst['Lam'][1] |
---|
2919 | dmin = wave/(2.0*sind(limits[1]/2.0)) |
---|
2920 | pfx = str(pId)+':'+str(hId)+':' |
---|
2921 | if Phases[phase]['General']['doPawley']: |
---|
2922 | hapDict[pfx+'LeBail'] = False #Pawley supercedes LeBail |
---|
2923 | Tmin = G2lat.Dsp2pos(inst,dmin) |
---|
2924 | if 'T' in inst['Type'][1]: |
---|
2925 | limits[0] = max(limits[0],Tmin) |
---|
2926 | else: |
---|
2927 | limits[1] = min(limits[1],Tmin) |
---|
2928 | else: |
---|
2929 | hapDict[pfx+'LeBail'] = hapData.get('LeBail',False) |
---|
2930 | if Phases[phase]['General']['Type'] == 'magnetic': |
---|
2931 | dmin = max(dmin,Phases[phase]['General'].get('MagDmin',0.)) |
---|
2932 | for item in ['Scale','Extinction']: |
---|
2933 | hapDict[pfx+item] = hapData[item][0] |
---|
2934 | # if hapData[item][1] and not hapDict[pfx+'LeBail']: |
---|
2935 | if hapData[item][1]: |
---|
2936 | hapVary.append(pfx+item) |
---|
2937 | names = G2spc.HStrainNames(SGData) |
---|
2938 | HSvals = [] |
---|
2939 | for i,name in enumerate(names): |
---|
2940 | hapDict[pfx+name] = hapData['HStrain'][0][i] |
---|
2941 | HSvals.append(hapDict[pfx+name]) |
---|
2942 | if hapData['HStrain'][1][i]: |
---|
2943 | hapVary.append(pfx+name) |
---|
2944 | # Acorr = G2lat.AplusDij(A,hapData['HStrain'][0],SGData) |
---|
2945 | Acorr = A |
---|
2946 | #adjust A for the Dij rith here? |
---|
2947 | if 'Layer Disp' in hapData: |
---|
2948 | hapDict[pfx+'LayerDisp'] = hapData['Layer Disp'][0] |
---|
2949 | if hapData['Layer Disp'][1]: |
---|
2950 | hapVary.append(pfx+'LayerDisp') |
---|
2951 | else: |
---|
2952 | hapDict[pfx+'LayerDisp'] = 0.0 |
---|
2953 | if 'E' not in inst['Type'][0]: |
---|
2954 | controlDict[pfx+'poType'] = hapData['Pref.Ori.'][0] |
---|
2955 | if hapData['Pref.Ori.'][0] == 'MD': |
---|
2956 | hapDict[pfx+'MD'] = hapData['Pref.Ori.'][1] |
---|
2957 | controlDict[pfx+'MDAxis'] = hapData['Pref.Ori.'][3] |
---|
2958 | if hapData['Pref.Ori.'][2]: # and not hapDict[pfx+'LeBail']: |
---|
2959 | hapVary.append(pfx+'MD') |
---|
2960 | else: #'SH' spherical harmonics |
---|
2961 | controlDict[pfx+'SHord'] = hapData['Pref.Ori.'][4] |
---|
2962 | controlDict[pfx+'SHncof'] = len(hapData['Pref.Ori.'][5]) |
---|
2963 | controlDict[pfx+'SHnames'] = G2lat.GenSHCoeff(SGData['SGLaue'],'0',controlDict[pfx+'SHord'],False) |
---|
2964 | controlDict[pfx+'SHhkl'] = [] |
---|
2965 | try: #patch for old Pref.Ori. items |
---|
2966 | controlDict[pfx+'SHtoler'] = 0.1 |
---|
2967 | if hapData['Pref.Ori.'][6][0] != '': |
---|
2968 | controlDict[pfx+'SHhkl'] = [eval(a.replace(' ',',')) for a in hapData['Pref.Ori.'][6]] |
---|
2969 | controlDict[pfx+'SHtoler'] = hapData['Pref.Ori.'][7] |
---|
2970 | except IndexError: |
---|
2971 | pass |
---|
2972 | for item in hapData['Pref.Ori.'][5]: |
---|
2973 | hapDict[pfx+item] = hapData['Pref.Ori.'][5][item] |
---|
2974 | if hapData['Pref.Ori.'][2]: # and not hapDict[pfx+'LeBail']: |
---|
2975 | hapVary.append(pfx+item) |
---|
2976 | for item in ['Mustrain','Size']: |
---|
2977 | controlDict[pfx+item+'Type'] = hapData[item][0] |
---|
2978 | hapDict[pfx+item+';mx'] = hapData[item][1][2] |
---|
2979 | if hapData[item][2][2]: |
---|
2980 | hapVary.append(pfx+item+';mx') |
---|
2981 | if hapData[item][0] in ['isotropic','uniaxial']: |
---|
2982 | hapDict[pfx+item+';i'] = hapData[item][1][0] |
---|
2983 | if hapData[item][2][0]: |
---|
2984 | hapVary.append(pfx+item+';i') |
---|
2985 | if hapData[item][0] == 'uniaxial': |
---|
2986 | controlDict[pfx+item+'Axis'] = hapData[item][3] |
---|
2987 | hapDict[pfx+item+';a'] = hapData[item][1][1] |
---|
2988 | if hapData[item][2][1]: |
---|
2989 | hapVary.append(pfx+item+';a') |
---|
2990 | else: #generalized for mustrain or ellipsoidal for size |
---|
2991 | Nterms = len(hapData[item][4]) |
---|
2992 | if item == 'Mustrain': |
---|
2993 | names = G2spc.MustrainNames(SGData) |
---|
2994 | pwrs = [] |
---|
2995 | for name in names: |
---|
2996 | h,k,l = name[1:] |
---|
2997 | pwrs.append([int(h),int(k),int(l)]) |
---|
2998 | controlDict[pfx+'MuPwrs'] = pwrs |
---|
2999 | for i in range(Nterms): |
---|
3000 | sfx = ';'+str(i) |
---|
3001 | hapDict[pfx+item+sfx] = hapData[item][4][i] |
---|
3002 | if hapData[item][5][i]: |
---|
3003 | hapVary.append(pfx+item+sfx) |
---|
3004 | if Phases[phase]['General']['Type'] != 'magnetic': |
---|
3005 | for bab in ['BabA','BabU']: |
---|
3006 | hapDict[pfx+bab] = hapData['Babinet'][bab][0] |
---|
3007 | if hapData['Babinet'][bab][1]: # and not hapDict[pfx+'LeBail']: |
---|
3008 | hapVary.append(pfx+bab) |
---|
3009 | |
---|
3010 | if Print: |
---|
3011 | pFile.write('\n Phase: %s in histogram: %s\n'%(phase,histogram)) |
---|
3012 | pFile.write(135*'='+'\n') |
---|
3013 | if hapDict.get(pfx+'LeBail'): |
---|
3014 | pFile.write(' Perform LeBail extraction\n') |
---|
3015 | elif 'E' not in inst['Type'][0]: |
---|
3016 | pFile.write(' Phase fraction : %10.4g Refine? %s\n'%(hapData['Scale'][0],hapData['Scale'][1])) |
---|
3017 | pFile.write(' Extinction coeff: %10.4f Refine? %s\n'%(hapData['Extinction'][0],hapData['Extinction'][1])) |
---|
3018 | if hapData['Pref.Ori.'][0] == 'MD': |
---|
3019 | Ax = hapData['Pref.Ori.'][3] |
---|
3020 | pFile.write(' March-Dollase PO: %10.4f Refine? %s Axis: %d %d %d\n'% |
---|
3021 | (hapData['Pref.Ori.'][1],hapData['Pref.Ori.'][2],Ax[0],Ax[1],Ax[2])) |
---|
3022 | else: #'SH' for spherical harmonics |
---|
3023 | PrintSHPO(hapData['Pref.Ori.']) |
---|
3024 | pFile.write(' Penalty hkl list: %s tolerance: %.2f\n'%(controlDict[pfx+'SHhkl'],controlDict[pfx+'SHtoler'])) |
---|
3025 | if 'E' not in inst['Type'][0]: |
---|
3026 | PrintSize(hapData['Size']) |
---|
3027 | PrintMuStrain(hapData['Mustrain'],SGData) |
---|
3028 | PrintHStrain(hapData['HStrain'],SGData) |
---|
3029 | if 'Layer Disp' in hapData: |
---|
3030 | pFile.write(' Layer Displacement: %10.3f Refine? %s\n'%(hapData['Layer Disp'][0],hapData['Layer Disp'][1])) |
---|
3031 | if Phases[phase]['General']['Type'] != 'magnetic'and 'E' not in inst['Type'][0]: |
---|
3032 | if hapData['Babinet']['BabA'][0]: |
---|
3033 | PrintBabinet(hapData['Babinet']) |
---|
3034 | if resetRefList and (not hapDict.get(pfx+'LeBail') or (hapData.get('LeBail',False) and Controls.get('newLeBail',False))): |
---|
3035 | Scale = Histogram['Sample Parameters']['Scale'][0] #for initializing reflection structure factors. |
---|
3036 | StartI = hapData['Scale'][0]*Scale |
---|
3037 | refList = [] |
---|
3038 | useExt = 'magnetic' in Phases[phase]['General']['Type'] and 'N' in inst['Type'][0] |
---|
3039 | if Phases[phase]['General'].get('Modulated',False): |
---|
3040 | ifSuper = True |
---|
3041 | HKLd = np.array(G2lat.GenSSHLaue(dmin,SGData,SSGData,Vec,maxH,Acorr)) |
---|
3042 | HKLd = G2mth.sortArray(HKLd,4,reverse=True) |
---|
3043 | for h,k,l,m,d in HKLd: |
---|
3044 | ext,mul,uniq,phi = G2spc.GenHKLf([h,k,l],SGData) |
---|
3045 | mul *= 2 # for powder overlap of Friedel pairs |
---|
3046 | if m or not ext or useExt: |
---|
3047 | if 'C' in inst['Type'][0]: |
---|
3048 | pos = G2lat.Dsp2pos(inst,d) |
---|
3049 | if limits[0] < pos < limits[1]: |
---|
3050 | refList.append([h,k,l,m,mul,d, pos,0.0,0.0,0.0,StartI, 0.0,0.0,1.0,1.0,1.0]) |
---|
3051 | #... sig,gam,fotsq,fctsq, phase,icorr,prfo,abs,ext |
---|
3052 | elif 'T' in inst['Type'][0]: |
---|
3053 | pos = G2lat.Dsp2pos(inst,d) |
---|
3054 | if limits[0] < pos < limits[1]: |
---|
3055 | wave = inst['difC'][1]*d/(252.816*inst['fltPath'][0]) |
---|
3056 | refList.append([h,k,l,m,mul,d, pos,0.0,0.0,0.0,StartI, 0.0,0.0,0.0,0.0,wave, 1.0,1.0,1.0]) |
---|
3057 | # ... sig,gam,fotsq,fctsq, phase,icorr,alp,bet,wave, prfo,abs,ext |
---|
3058 | #TODO - if tabulated put alp & bet in here |
---|
3059 | elif 'B' in inst['Type'][0]: |
---|
3060 | pos = G2lat.Dsp2pos(inst,d) |
---|
3061 | if limits[0] < pos < limits[1]: |
---|
3062 | refList.append([h,k,l,m,mul,d, pos,0.0,0.0,0.0,StartI, 0.0,0.0,0.0,0.0, 1.0,1.0,1.0]) |
---|
3063 | # ... sig,gam,fotsq,fctsq, phase,icorr,alp,bet, prfo,abs,ext |
---|
3064 | else: |
---|
3065 | ifSuper = False |
---|
3066 | HKLd = np.array(G2lat.GenHLaue(dmin,SGData,Acorr)) |
---|
3067 | HKLd = G2mth.sortArray(HKLd,3,reverse=True) |
---|
3068 | for h,k,l,d in HKLd: |
---|
3069 | ext,mul,uniq,phi = G2spc.GenHKLf([h,k,l],SGData) |
---|
3070 | if ext and 'N' in inst['Type'][0] and 'MagSpGrp' in SGData: |
---|
3071 | ext = G2spc.checkMagextc([h,k,l],SGData) |
---|
3072 | mul *= 2 # for powder overlap of Friedel pairs |
---|
3073 | if ext and not useExt: |
---|
3074 | continue |
---|
3075 | if 'C' in inst['Type'][0]: |
---|
3076 | pos = G2lat.Dsp2pos(inst,d) |
---|
3077 | if limits[0] < pos < limits[1]: |
---|
3078 | refList.append([h,k,l,mul,d, pos,0.0,0.0,0.0,StartI, 0.0,0.0,1.0,1.0,1.0]) |
---|
3079 | #... sig,gam,fotsq,fctsq, phase,icorr,prfo,abs,ext |
---|
3080 | elif 'T' in inst['Type'][0]: |
---|
3081 | pos = G2lat.Dsp2pos(inst,d) |
---|
3082 | if limits[0] < pos < limits[1]: |
---|
3083 | wave = inst['difC'][1]*d/(252.816*inst['fltPath'][0]) |
---|
3084 | refList.append([h,k,l,mul,d, pos,0.0,0.0,0.0,StartI, 0.0,0.0,0.0,0.0,wave, 1.0,1.0,1.0]) |
---|
3085 | # ... sig,gam,fotsq,fctsq, phase,icorr,alp,bet,wave, prfo,abs,ext |
---|
3086 | elif 'B' in inst['Type'][0]: |
---|
3087 | pos = G2lat.Dsp2pos(inst,d) |
---|
3088 | if limits[0] < pos < limits[1]: |
---|
3089 | refList.append([h,k,l,mul,d, pos,0.0,0.0,0.0,StartI, 0.0,0.0,0.0,0.0, 1.0,1.0,1.0]) |
---|
3090 | # ... sig,gam,fotsq,fctsq, phase,icorr,alp,bet, prfo,abs,ext |
---|
3091 | elif 'E' in inst['Type'][0]: |
---|
3092 | pos = G2lat.Dsp2pos(inst,d) |
---|
3093 | if limits[0] < pos < limits[1]: |
---|
3094 | refList.append([h,k,l,mul,d, pos,0.0,0.0,0.0,StartI, 0.0,0.0]) |
---|
3095 | # ... sig,gam,fotsq,fctsq, phase,icorr |
---|
3096 | Histogram['Reflection Lists'][phase] = {'RefList':np.array(refList),'FF':{},'Type':inst['Type'][0],'Super':ifSuper} |
---|
3097 | elif 'HKLF' in histogram: |
---|
3098 | inst = Histogram['Instrument Parameters'][0] |
---|
3099 | hId = Histogram['hId'] |
---|
3100 | hfx = ':%d:'%(hId) |
---|
3101 | for item in inst: |
---|
3102 | if type(inst) is not list and item != 'Type': continue # skip over non-refined items (such as InstName) |
---|
3103 | hapDict[hfx+item] = inst[item][1] |
---|
3104 | pfx = str(pId)+':'+str(hId)+':' |
---|
3105 | hapDict[pfx+'Scale'] = hapData['Scale'][0] |
---|
3106 | if hapData['Scale'][1]: |
---|
3107 | hapVary.append(pfx+'Scale') |
---|
3108 | |
---|
3109 | extApprox,extType,extParms = hapData['Extinction'] |
---|
3110 | controlDict[pfx+'EType'] = extType |
---|
3111 | controlDict[pfx+'EApprox'] = extApprox |
---|
3112 | if 'C' in inst['Type'][0]: |
---|
3113 | controlDict[pfx+'Tbar'] = extParms['Tbar'] |
---|
3114 | controlDict[pfx+'Cos2TM'] = extParms['Cos2TM'] |
---|
3115 | if 'Primary' in extType: |
---|
3116 | Ekey = ['Ep',] |
---|
3117 | elif 'I & II' in extType: |
---|
3118 | Ekey = ['Eg','Es'] |
---|
3119 | elif 'Secondary Type II' == extType: |
---|
3120 | Ekey = ['Es',] |
---|
3121 | elif 'Secondary Type I' == extType: |
---|
3122 | Ekey = ['Eg',] |
---|
3123 | else: #'None' |
---|
3124 | Ekey = [] |
---|
3125 | for eKey in Ekey: |
---|
3126 | hapDict[pfx+eKey] = extParms[eKey][0] |
---|
3127 | if extParms[eKey][1]: |
---|
3128 | hapVary.append(pfx+eKey) |
---|
3129 | for bab in ['BabA','BabU']: |
---|
3130 | hapDict[pfx+bab] = hapData['Babinet'][bab][0] |
---|
3131 | if hapData['Babinet'][bab][1]: |
---|
3132 | hapVary.append(pfx+bab) |
---|
3133 | Twins = hapData.get('Twins',[[np.array([[1,0,0],[0,1,0],[0,0,1]]),[1.0,False,0]],]) |
---|
3134 | if len(Twins) == 1: |
---|
3135 | hapDict[pfx+'Flack'] = hapData.get('Flack',[0.,False])[0] |
---|
3136 | if hapData.get('Flack',[0,False])[1]: |
---|
3137 | hapVary.append(pfx+'Flack') |
---|
3138 | sumTwFr = 0. |
---|
3139 | controlDict[pfx+'TwinLaw'] = [] |
---|
3140 | controlDict[pfx+'TwinInv'] = [] |
---|
3141 | NTL = 0 |
---|
3142 | for it,twin in enumerate(Twins): |
---|
3143 | if 'bool' in str(type(twin[0])): |
---|
3144 | controlDict[pfx+'TwinInv'].append(twin[0]) |
---|
3145 | controlDict[pfx+'TwinLaw'].append(np.zeros((3,3))) |
---|
3146 | else: |
---|
3147 | NTL += 1 |
---|
3148 | controlDict[pfx+'TwinInv'].append(False) |
---|
3149 | controlDict[pfx+'TwinLaw'].append(twin[0]) |
---|
3150 | if it: |
---|
3151 | hapDict[pfx+'TwinFr:'+str(it)] = twin[1] |
---|
3152 | sumTwFr += twin[1] |
---|
3153 | else: |
---|
3154 | hapDict[pfx+'TwinFr:0'] = twin[1][0] |
---|
3155 | controlDict[pfx+'TwinNMN'] = twin[1][1] |
---|
3156 | if Twins[0][1][1]: |
---|
3157 | hapVary.append(pfx+'TwinFr:'+str(it)) |
---|
3158 | controlDict[pfx+'NTL'] = NTL |
---|
3159 | #need to make constraint on TwinFr |
---|
3160 | controlDict[pfx+'TwinLaw'] = np.array(controlDict[pfx+'TwinLaw']) |
---|
3161 | if len(Twins) > 1: #force sum to unity |
---|
3162 | hapDict[pfx+'TwinFr:0'] = 1.-sumTwFr |
---|
3163 | if Print: |
---|
3164 | pFile.write('\n Phase: %s in histogram: %s\n'%(phase,histogram)) |
---|
3165 | pFile.write(135*'='+'\n') |
---|
3166 | pFile.write(' Scale factor : %10.4g Refine? %s\n'%(hapData['Scale'][0],hapData['Scale'][1])) |
---|
3167 | if extType != 'None': |
---|
3168 | pFile.write(' Extinction Type: %15s approx: %10s\n'%(extType,extApprox)) |
---|
3169 | text = ' Parameters :' |
---|
3170 | for eKey in Ekey: |
---|
3171 | text += ' %4s : %10.3e Refine? '%(eKey,extParms[eKey][0])+str(extParms[eKey][1]) |
---|
3172 | pFile.write(text+'\n') |
---|
3173 | if hapData['Babinet']['BabA'][0]: |
---|
3174 | PrintBabinet(hapData['Babinet']) |
---|
3175 | if not SGData['SGInv'] and len(Twins) == 1: |
---|
3176 | pFile.write(' Flack parameter: %10.3f Refine? %s\n'%(hapData['Flack'][0],hapData['Flack'][1])) |
---|
3177 | if len(Twins) > 1: |
---|
3178 | for it,twin in enumerate(Twins): |
---|
3179 | if 'bool' in str(type(twin[0])): |
---|
3180 | pFile.write(' Nonmerohedral twin fr.: %5.3f Inv? %s Refine? %s\n'% |
---|
3181 | (hapDict[pfx+'TwinFr:'+str(it)],str(controlDict[pfx+'TwinInv'][it]),str(Twins[0][1][1]))) |
---|
3182 | else: |
---|
3183 | pFile.write(' Twin law: %s Twin fr.: %5.3f Refine? %s\n'% |
---|
3184 | (str(twin[0]).replace('\n',','),hapDict[pfx+'TwinFr:'+str(it)],str(Twins[0][1][1]))) |
---|
3185 | |
---|
3186 | Histogram['Reflection Lists'] = phase |
---|
3187 | |
---|
3188 | return hapVary,hapDict,controlDict |
---|
3189 | |
---|
3190 | def SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms,calcControls,Print=True,pFile=None, |
---|
3191 | covMatrix=[],varyList=[]): |
---|
3192 | '''Updates parmDict with HAP results from refinement and prints a summary if Print is True |
---|
3193 | ''' |
---|
3194 | |
---|
3195 | def PrintSizeAndSig(hapData,sizeSig): |
---|
3196 | line = '\n Size model: %9s'%(hapData[0]) |
---|
3197 | refine = False |
---|
3198 | if hapData[0] in ['isotropic','uniaxial']: |
---|
3199 | line += ' equatorial:%12.5f'%(hapData[1][0]) |
---|
3200 | if sizeSig[0][0]: |
---|
3201 | line += ', sig:%8.4f'%(sizeSig[0][0]) |
---|
3202 | refine = True |
---|
3203 | if hapData[0] == 'uniaxial': |
---|
3204 | line += ' axial:%12.4f'%(hapData[1][1]) |
---|
3205 | if sizeSig[0][1]: |
---|
3206 | refine = True |
---|
3207 | line += ', sig:%8.4f'%(sizeSig[0][1]) |
---|
3208 | line += ' LG mix coeff.:%12.4f'%(hapData[1][2]) |
---|
3209 | if sizeSig[0][2]: |
---|
3210 | refine = True |
---|
3211 | line += ', sig:%8.4f'%(sizeSig[0][2]) |
---|
3212 | if refine: |
---|
3213 | pFile.write(line+'\n') |
---|
3214 | else: |
---|
3215 | line += ' LG mix coeff.:%12.4f'%(hapData[1][2]) |
---|
3216 | if sizeSig[0][2]: |
---|
3217 | refine = True |
---|
3218 | line += ', sig:%8.4f'%(sizeSig[0][2]) |
---|
3219 | Snames = ['S11','S22','S33','S12','S13','S23'] |
---|
3220 | ptlbls = ' name :' |
---|
3221 | ptstr = ' value :' |
---|
3222 | sigstr = ' sig :' |
---|
3223 | for i,name in enumerate(Snames): |
---|
3224 | ptlbls += '%12s' % (name) |
---|
3225 | ptstr += '%12.6f' % (hapData[4][i]) |
---|
3226 | if sizeSig[1][i]: |
---|
3227 | refine = True |
---|
3228 | sigstr += '%12.6f' % (sizeSig[1][i]) |
---|
3229 | else: |
---|
3230 | sigstr += 12*' ' |
---|
3231 | if refine: |
---|
3232 | pFile.write(line+'\n') |
---|
3233 | pFile.write(ptlbls+'\n') |
---|
3234 | pFile.write(ptstr+'\n') |
---|
3235 | pFile.write(sigstr+'\n') |
---|
3236 | |
---|
3237 | def PrintMuStrainAndSig(hapData,mustrainSig,SGData): |
---|
3238 | line = '\n Mustrain model: %9s\n'%(hapData[0]) |
---|
3239 | refine = False |
---|
3240 | if hapData[0] in ['isotropic','uniaxial']: |
---|
3241 | line += ' equatorial:%12.1f'%(hapData[1][0]) |
---|
3242 | if mustrainSig[0][0]: |
---|
3243 | line += ', sig:%8.1f'%(mustrainSig[0][0]) |
---|
3244 | refine = True |
---|
3245 | if hapData[0] == 'uniaxial': |
---|
3246 | line += ' axial:%12.1f'%(hapData[1][1]) |
---|
3247 | if mustrainSig[0][1]: |
---|
3248 | line += ', sig:%8.1f'%(mustrainSig[0][1]) |
---|
3249 | line += ' LG mix coeff.:%12.4f'%(hapData[1][2]) |
---|
3250 | if mustrainSig[0][2]: |
---|
3251 | refine = True |
---|
3252 | line += ', sig:%8.3f'%(mustrainSig[0][2]) |
---|
3253 | if refine: |
---|
3254 | pFile.write(line+'\n') |
---|
3255 | else: |
---|
3256 | line += ' LG mix coeff.:%12.4f'%(hapData[1][2]) |
---|
3257 | if mustrainSig[0][2]: |
---|
3258 | refine = True |
---|
3259 | line += ', sig:%8.3f'%(mustrainSig[0][2]) |
---|
3260 | Snames = G2spc.MustrainNames(SGData) |
---|
3261 | ptlbls = ' name :' |
---|
3262 | ptstr = ' value :' |
---|
3263 | sigstr = ' sig :' |
---|
3264 | for i,name in enumerate(Snames): |
---|
3265 | ptlbls += '%12s' % (name) |
---|
3266 | ptstr += '%12.1f' % (hapData[4][i]) |
---|
3267 | if mustrainSig[1][i]: |
---|
3268 | refine = True |
---|
3269 | sigstr += '%12.1f' % (mustrainSig[1][i]) |
---|
3270 | else: |
---|
3271 | sigstr += 12*' ' |
---|
3272 | if refine: |
---|
3273 | pFile.write(line+'\n') |
---|
3274 | pFile.write(ptlbls+'\n') |
---|
3275 | pFile.write(ptstr+'\n') |
---|
3276 | pFile.write(sigstr+'\n') |
---|
3277 | |
---|
3278 | def PrintHStrainAndSig(hapData,strainSig,SGData): |
---|
3279 | Hsnames = G2spc.HStrainNames(SGData) |
---|
3280 | ptlbls = ' name :' |
---|
3281 | ptstr = ' value :' |
---|
3282 | sigstr = ' sig :' |
---|
3283 | refine = False |
---|
3284 | for i,name in enumerate(Hsnames): |
---|
3285 | ptlbls += '%12s' % (name) |
---|
3286 | ptstr += '%12.4g' % (hapData[0][i]) |
---|
3287 | if name in strainSig: |
---|
3288 | refine = True |
---|
3289 | sigstr += '%12.4g' % (strainSig[name]) |
---|
3290 | else: |
---|
3291 | sigstr += 12*' ' |
---|
3292 | if refine: |
---|
3293 | pFile.write('\n Hydrostatic/elastic strain:\n') |
---|
3294 | pFile.write(ptlbls+'\n') |
---|
3295 | pFile.write(ptstr+'\n') |
---|
3296 | pFile.write(sigstr+'\n') |
---|
3297 | |
---|
3298 | def PrintSHPOAndSig(pfx,hapData,POsig): |
---|
3299 | Tindx = 1.0 |
---|
3300 | Tvar = 0.0 |
---|
3301 | pFile.write('\n Spherical harmonics preferred orientation: Order: %d\n'%hapData[4]) |
---|
3302 | ptlbls = ' names :' |
---|
3303 | ptstr = ' values:' |
---|
3304 | sigstr = ' sig :' |
---|
3305 | for item in hapData[5]: |
---|
3306 | ptlbls += '%12s'%(item) |
---|
3307 | ptstr += '%12.3f'%(hapData[5][item]) |
---|
3308 | l = 2.0*eval(item.strip('C'))[0]+1 |
---|
3309 | Tindx += hapData[5][item]**2/l |
---|
3310 | if pfx+item in POsig: |
---|
3311 | Tvar += (2.*hapData[5][item]*POsig[pfx+item]/l)**2 |
---|
3312 | sigstr += '%12.3f'%(POsig[pfx+item]) |
---|
3313 | else: |
---|
3314 | sigstr += 12*' ' |
---|
3315 | pFile.write(ptlbls+'\n') |
---|
3316 | pFile.write(ptstr+'\n') |
---|
3317 | pFile.write(sigstr+'\n') |
---|
3318 | pFile.write('\n Texture index J = %.3f(%d)\n'%(Tindx,int(1000*np.sqrt(Tvar)))) |
---|
3319 | |
---|
3320 | def PrintExtAndSig(pfx,hapData,ScalExtSig): |
---|
3321 | pFile.write('\n Single crystal extinction: Type: %s Approx: %s\n'%(hapData[0],hapData[1])) |
---|
3322 | text = '' |
---|
3323 | for item in hapData[2]: |
---|
3324 | if pfx+item in ScalExtSig: |
---|
3325 | text += ' %s: '%(item) |
---|
3326 | text += '%12.2e'%(hapData[2][item][0]) |
---|
3327 | if pfx+item in ScalExtSig: |
---|
3328 | text += ' sig: %12.2e'%(ScalExtSig[pfx+item]) |
---|
3329 | pFile.write(text+'\n') |
---|
3330 | |
---|
3331 | def PrintBabinetAndSig(pfx,hapData,BabSig): |
---|
3332 | pFile.write('\n Babinet form factor modification:\n') |
---|
3333 | ptlbls = ' names :' |
---|
3334 | ptstr = ' values:' |
---|
3335 | sigstr = ' sig :' |
---|
3336 | for item in hapData: |
---|
3337 | ptlbls += '%12s'%(item) |
---|
3338 | ptstr += '%12.3f'%(hapData[item][0]) |
---|
3339 | if pfx+item in BabSig: |
---|
3340 | sigstr += '%12.3f'%(BabSig[pfx+item]) |
---|
3341 | else: |
---|
3342 | sigstr += 12*' ' |
---|
3343 | pFile.write(ptlbls+'\n') |
---|
3344 | pFile.write(ptstr+'\n') |
---|
3345 | pFile.write(sigstr+'\n') |
---|
3346 | |
---|
3347 | def PrintTwinsAndSig(pfx,twinData,TwinSig): |
---|
3348 | pFile.write('\n Twin Law fractions :\n') |
---|
3349 | ptlbls = ' names :' |
---|
3350 | ptstr = ' values:' |
---|
3351 | sigstr = ' sig :' |
---|
3352 | for it,item in enumerate(twinData): |
---|
3353 | ptlbls += ' twin #%d'%(it) |
---|
3354 | if it: |
---|
3355 | ptstr += '%12.3f'%(item[1]) |
---|
3356 | else: |
---|
3357 | ptstr += '%12.3f'%(item[1][0]) |
---|
3358 | if pfx+'TwinFr:'+str(it) in TwinSig: |
---|
3359 | sigstr += '%12.3f'%(TwinSig[pfx+'TwinFr:'+str(it)]) |
---|
3360 | else: |
---|
3361 | sigstr += 12*' ' |
---|
3362 | pFile.write(ptlbls+'\n') |
---|
3363 | pFile.write(ptstr+'\n') |
---|
3364 | pFile.write(sigstr+'\n') |
---|
3365 | |
---|
3366 | # global PhFrExtPOSig # this is not used externally anymore. Remove? |
---|
3367 | PhFrExtPOSig = {} |
---|
3368 | SizeMuStrSig = {} |
---|
3369 | ScalExtSig = {} |
---|
3370 | BabSig = {} |
---|
3371 | TwinFrSig = {} |
---|
3372 | # wtFrSum = {} |
---|
3373 | for phase in Phases: |
---|
3374 | HistoPhase = Phases[phase]['Histograms'] |
---|
3375 | General = Phases[phase]['General'] |
---|
3376 | SGData = General['SGData'] |
---|
3377 | pId = Phases[phase]['pId'] |
---|
3378 | histoList = list(HistoPhase.keys()) |
---|
3379 | histoList.sort() |
---|
3380 | for histogram in histoList: |
---|
3381 | try: |
---|
3382 | Histogram = Histograms[histogram] |
---|
3383 | except KeyError: |
---|
3384 | #skip if histogram not included e.g. in a sequential refinement |
---|
3385 | continue |
---|
3386 | if not Phases[phase]['Histograms'][histogram]['Use']: |
---|
3387 | #skip if phase absent from this histogram |
---|
3388 | continue |
---|
3389 | hapData = HistoPhase[histogram] |
---|
3390 | hId = Histogram['hId'] |
---|
3391 | pfx = str(pId)+':'+str(hId)+':' |
---|
3392 | # if hId not in wtFrSum: |
---|
3393 | # wtFrSum[hId] = 0. |
---|
3394 | if 'PWDR' in histogram: |
---|
3395 | inst = Histogram['Instrument Parameters'][0] #TODO - grab table here if present |
---|
3396 | if 'E' not in inst['Type'][0]: |
---|
3397 | parmDict[pfx+'Scale'] = max(1.e-12,parmDict[pfx+'Scale']) |
---|
3398 | for item in ['Scale','Extinction']: |
---|
3399 | hapData[item][0] = parmDict[pfx+item] |
---|
3400 | if pfx+item in sigDict and not parmDict.get(pfx+'LeBail'): |
---|
3401 | PhFrExtPOSig.update({pfx+item:sigDict[pfx+item],}) |
---|
3402 | if hapData['Pref.Ori.'][0] == 'MD': |
---|
3403 | hapData['Pref.Ori.'][1] = parmDict[pfx+'MD'] |
---|
3404 | if pfx+'MD' in sigDict and not parmDict.get(pfx+'LeBail'): |
---|
3405 | PhFrExtPOSig.update({pfx+'MD':sigDict[pfx+'MD'],}) |
---|
3406 | else: #'SH' spherical harmonics |
---|
3407 | for item in hapData['Pref.Ori.'][5]: |
---|
3408 | hapData['Pref.Ori.'][5][item] = parmDict[pfx+item] |
---|
3409 | if pfx+item in sigDict and not parmDict.get(pfx+'LeBail'): |
---|
3410 | PhFrExtPOSig.update({pfx+item:sigDict[pfx+item],}) |
---|
3411 | SizeMuStrSig.update({pfx+'Mustrain':[[0,0,0],[0 for i in range(len(hapData['Mustrain'][4]))]], |
---|
3412 | pfx+'Size':[[0,0,0],[0 for i in range(len(hapData['Size'][4]))]],pfx+'HStrain':{}}) |
---|
3413 | for item in ['Mustrain','Size']: |
---|
3414 | hapData[item][1][2] = parmDict[pfx+item+';mx'] |
---|
3415 | # hapData[item][1][2] = min(1.,max(0.,hapData[item][1][2])) |
---|
3416 | if pfx+item+';mx' in sigDict: |
---|
3417 | SizeMuStrSig[pfx+item][0][2] = sigDict[pfx+item+';mx'] |
---|
3418 | if hapData[item][0] in ['isotropic','uniaxial']: |
---|
3419 | hapData[item][1][0] = parmDict[pfx+item+';i'] |
---|
3420 | if item == 'Size': |
---|
3421 | hapData[item][1][0] = min(10.,max(0.001,hapData[item][1][0])) |
---|
3422 | if pfx+item+';i' in sigDict: |
---|
3423 | SizeMuStrSig[pfx+item][0][0] = sigDict[pfx+item+';i'] |
---|
3424 | if hapData[item][0] == 'uniaxial': |
---|
3425 | hapData[item][1][1] = parmDict[pfx+item+';a'] |
---|
3426 | if item == 'Size': |
---|
3427 | hapData[item][1][1] = min(10.,max(0.001,hapData[item][1][1])) |
---|
3428 | if pfx+item+';a' in sigDict: |
---|
3429 | SizeMuStrSig[pfx+item][0][1] = sigDict[pfx+item+';a'] |
---|
3430 | else: #generalized for mustrain or ellipsoidal for size |
---|
3431 | Nterms = len(hapData[item][4]) |
---|
3432 | for i in range(Nterms): |
---|
3433 | sfx = ';'+str(i) |
---|
3434 | hapData[item][4][i] = parmDict[pfx+item+sfx] |
---|
3435 | if pfx+item+sfx in sigDict: |
---|
3436 | SizeMuStrSig[pfx+item][1][i] = sigDict[pfx+item+sfx] |
---|
3437 | else: |
---|
3438 | SizeMuStrSig.update({pfx+'HStrain':{}}) |
---|
3439 | names = G2spc.HStrainNames(SGData) |
---|
3440 | for i,name in enumerate(names): |
---|
3441 | hapData['HStrain'][0][i] = parmDict[pfx+name] |
---|
3442 | if pfx+name in sigDict: |
---|
3443 | SizeMuStrSig[pfx+'HStrain'][name] = sigDict[pfx+name] |
---|
3444 | if 'Layer Disp' in hapData: |
---|
3445 | hapData['Layer Disp'][0] = parmDict[pfx+'LayerDisp'] |
---|
3446 | if pfx+'LayerDisp' in sigDict: |
---|
3447 | SizeMuStrSig[pfx+'LayerDisp'] = sigDict[pfx+'LayerDisp'] |
---|
3448 | if Phases[phase]['General']['Type'] != 'magnetic' and 'E' not in inst['Type'][0]: |
---|
3449 | for name in ['BabA','BabU']: |
---|
3450 | hapData['Babinet'][name][0] = parmDict[pfx+name] |
---|
3451 | if pfx+name in sigDict and not parmDict.get(pfx+'LeBail'): |
---|
3452 | BabSig[pfx+name] = sigDict[pfx+name] |
---|
3453 | |
---|
3454 | elif 'HKLF' in histogram: |
---|
3455 | for item in ['Scale','Flack']: |
---|
3456 | if parmDict.get(pfx+item): |
---|
3457 | hapData[item][0] = parmDict[pfx+item] |
---|
3458 | if pfx+item in sigDict: |
---|
3459 | ScalExtSig[pfx+item] = sigDict[pfx+item] |
---|
3460 | for item in ['Ep','Eg','Es']: |
---|
3461 | if parmDict.get(pfx+item): |
---|
3462 | hapData['Extinction'][2][item][0] = parmDict[pfx+item] |
---|
3463 | if pfx+item in sigDict: |
---|
3464 | ScalExtSig[pfx+item] = sigDict[pfx+item] |
---|
3465 | for item in ['BabA','BabU']: |
---|
3466 | hapData['Babinet'][item][0] = parmDict[pfx+item] |
---|
3467 | if pfx+item in sigDict: |
---|
3468 | BabSig[pfx+item] = sigDict[pfx+item] |
---|
3469 | if 'Twins' in hapData: |
---|
3470 | it = 1 |
---|
3471 | sumTwFr = 0. |
---|
3472 | while True: |
---|
3473 | try: |
---|
3474 | hapData['Twins'][it][1] = parmDict[pfx+'TwinFr:'+str(it)] |
---|
3475 | if pfx+'TwinFr:'+str(it) in sigDict: |
---|
3476 | TwinFrSig[pfx+'TwinFr:'+str(it)] = sigDict[pfx+'TwinFr:'+str(it)] |
---|
3477 | if it: |
---|
3478 | sumTwFr += hapData['Twins'][it][1] |
---|
3479 | it += 1 |
---|
3480 | except KeyError: |
---|
3481 | break |
---|
3482 | hapData['Twins'][0][1][0] = 1.-sumTwFr |
---|
3483 | # HAP parameters updated, now compute derived quantities |
---|
3484 | for hist in Phases[phase]['Histograms']: |
---|
3485 | try: |
---|
3486 | Histogram = Histograms[hist] |
---|
3487 | except KeyError: #skip if histogram not included e.g. in a sequential refinement |
---|
3488 | continue |
---|
3489 | vDict,sDict = G2stMth.calcMassFracs(varyList,covMatrix, |
---|
3490 | Phases,hist,Histograms[hist]['hId']) |
---|
3491 | parmDict.update(vDict) |
---|
3492 | sigDict.update(sDict) |
---|
3493 | |
---|
3494 | if Print: |
---|
3495 | for phase in Phases: |
---|
3496 | HistoPhase = Phases[phase]['Histograms'] |
---|
3497 | General = Phases[phase]['General'] |
---|
3498 | SGData = General['SGData'] |
---|
3499 | pId = Phases[phase]['pId'] |
---|
3500 | histoList = list(HistoPhase.keys()) |
---|
3501 | histoList.sort() |
---|
3502 | for histogram in histoList: |
---|
3503 | try: |
---|
3504 | Histogram = Histograms[histogram] |
---|
3505 | except KeyError: |
---|
3506 | #skip if histogram not included e.g. in a sequential refinement |
---|
3507 | continue |
---|
3508 | hapData = HistoPhase[histogram] |
---|
3509 | hId = Histogram['hId'] |
---|
3510 | Histogram['Residuals'][str(pId)+'::Name'] = phase |
---|
3511 | pfx = str(pId)+':'+str(hId)+':' |
---|
3512 | hfx = ':%s:'%(hId) |
---|
3513 | if pfx+'Nref' not in Histogram['Residuals']: #skip not used phase in histogram |
---|
3514 | continue |
---|
3515 | pFile.write('\n Phase: %s in histogram: %s\n'%(phase,histogram)) |
---|
3516 | pFile.write(135*'='+'\n') |
---|
3517 | Inst = Histogram['Instrument Parameters'][0] |
---|
3518 | if 'PWDR' in histogram: |
---|
3519 | pFile.write(' Final refinement RF, RF^2 = %.2f%%, %.2f%% on %d reflections\n'% |
---|
3520 | (Histogram['Residuals'][pfx+'Rf'],Histogram['Residuals'][pfx+'Rf^2'],Histogram['Residuals'][pfx+'Nref'])) |
---|
3521 | pFile.write(' Durbin-Watson statistic = %.3f\n'%(Histogram['Residuals']['Durbin-Watson'])) |
---|
3522 | pFile.write(' Bragg intensity sum = %.3g\n'%(Histogram['Residuals'][pfx+'sumInt'])) |
---|
3523 | |
---|
3524 | if parmDict.get(pfx+'LeBail') or 'E' in Inst['Type'][0]: |
---|
3525 | pFile.write(' Performed LeBail extraction for phase %s in histogram %s\n'%(phase,histogram)) |
---|
3526 | elif 'E' not in Inst['Type'][0]: |
---|
3527 | var = pfx + 'WgtFrac' |
---|
3528 | if pfx+'Scale' in PhFrExtPOSig or var in sigDict: |
---|
3529 | wtFr = parmDict.get(var,0) |
---|
3530 | sigwtFr = sigDict.get(var,0) |
---|
3531 | pFile.write(' Phase fraction : %10.5g, sig %10.5g Weight fraction : %8.5f, sig %10.5f\n'% |
---|
3532 | (hapData['Scale'][0],PhFrExtPOSig.get(pfx+'Scale',0),wtFr,sigwtFr)) |
---|
3533 | if pfx+'Extinction' in PhFrExtPOSig: |
---|
3534 | pFile.write(' Extinction coeff: %10.4f, sig %10.4f\n'%(hapData['Extinction'][0],PhFrExtPOSig[pfx+'Extinction'])) |
---|
3535 | if hapData['Pref.Ori.'][0] == 'MD': |
---|
3536 | if pfx+'MD' in PhFrExtPOSig: |
---|
3537 | pFile.write(' March-Dollase PO: %10.4f, sig %10.4f\n'%(hapData['Pref.Ori.'][1],PhFrExtPOSig[pfx+'MD'])) |
---|
3538 | else: |
---|
3539 | PrintSHPOAndSig(pfx,hapData['Pref.Ori.'],PhFrExtPOSig) |
---|
3540 | if 'E' not in Inst['Type'][0]: |
---|
3541 | PrintSizeAndSig(hapData['Size'],SizeMuStrSig[pfx+'Size']) |
---|
3542 | PrintMuStrainAndSig(hapData['Mustrain'],SizeMuStrSig[pfx+'Mustrain'],SGData) |
---|
3543 | PrintHStrainAndSig(hapData['HStrain'],SizeMuStrSig[pfx+'HStrain'],SGData) |
---|
3544 | if pfx+'LayerDisp' in SizeMuStrSig: |
---|
3545 | pFile.write(' Layer displacement : %10.3f, sig %10.3f\n'%(hapData['Layer Disp'][0],SizeMuStrSig[pfx+'LayerDisp'])) |
---|
3546 | if Phases[phase]['General']['Type'] != 'magnetic' and not parmDict.get(pfx+'LeBail') and 'E' not in Inst['Type'][0]: |
---|
3547 | if len(BabSig): |
---|
3548 | PrintBabinetAndSig(pfx,hapData['Babinet'],BabSig) |
---|
3549 | |
---|
3550 | elif 'HKLF' in histogram: |
---|
3551 | pFile.write(' Final refinement RF, RF^2 = %.2f%%, %.2f%% on %d reflections (%d user rejected, %d sp.gp.extinct)\n'% |
---|
3552 | (Histogram['Residuals'][pfx+'Rf'],Histogram['Residuals'][pfx+'Rf^2'],Histogram['Residuals'][pfx+'Nref'], |
---|
3553 | Histogram['Residuals'][pfx+'Nrej'],Histogram['Residuals'][pfx+'Next'])) |
---|
3554 | if calcControls != None: #skipped in seqRefine |
---|
3555 | if 'X'in Inst['Type'][0]: |
---|
3556 | PrintFprime(calcControls['FFtables'],hfx,pFile) |
---|
3557 | elif 'NC' in Inst['Type'][0]: |
---|
3558 | PrintBlength(calcControls['BLtables'],Inst['Lam'][1],pFile) |
---|
3559 | pFile.write(' HKLF histogram weight factor = %.3f\n'%(Histogram['wtFactor'])) |
---|
3560 | if pfx+'Scale' in ScalExtSig: |
---|
3561 | pFile.write(' Scale factor : %10.4g, sig %10.4g\n'%(hapData['Scale'][0],ScalExtSig[pfx+'Scale'])) |
---|
3562 | if hapData['Extinction'][0] != 'None': |
---|
3563 | PrintExtAndSig(pfx,hapData['Extinction'],ScalExtSig) |
---|
3564 | if len(BabSig): |
---|
3565 | PrintBabinetAndSig(pfx,hapData['Babinet'],BabSig) |
---|
3566 | if pfx+'Flack' in ScalExtSig: |
---|
3567 | pFile.write(' Flack parameter : %10.3f, sig %10.3f\n'%(hapData['Flack'][0],ScalExtSig[pfx+'Flack'])) |
---|
3568 | if len(TwinFrSig): |
---|
3569 | PrintTwinsAndSig(pfx,hapData['Twins'],TwinFrSig) |
---|
3570 | |
---|
3571 | ################################################################################ |
---|
3572 | ##### Histogram data |
---|
3573 | ################################################################################ |
---|
3574 | |
---|
3575 | def GetHistogramData(Histograms,Print=True,pFile=None): |
---|
3576 | 'needs a doc string' |
---|
3577 | |
---|
3578 | def GetBackgroundParms(hId,Background): |
---|
3579 | Back = Background[0] |
---|
3580 | DebyePeaks = Background[1] |
---|
3581 | bakType,bakFlag = Back[:2] |
---|
3582 | backVals = Back[3:] |
---|
3583 | backNames = [':'+str(hId)+':Back;'+str(i) for i in range(len(backVals))] |
---|
3584 | backDict = dict(zip(backNames,backVals)) |
---|
3585 | backVary = [] |
---|
3586 | if bakFlag: |
---|
3587 | backVary = backNames |
---|
3588 | backDict[':'+str(hId)+':nDebye'] = DebyePeaks['nDebye'] |
---|
3589 | backDict[':'+str(hId)+':nPeaks'] = DebyePeaks['nPeaks'] |
---|
3590 | debyeDict = {} |
---|
3591 | debyeList = [] |
---|
3592 | for i in range(DebyePeaks['nDebye']): |
---|
3593 | debyeNames = [':'+str(hId)+':DebyeA;'+str(i),':'+str(hId)+':DebyeR;'+str(i),':'+str(hId)+':DebyeU;'+str(i)] |
---|
3594 | debyeDict.update(dict(zip(debyeNames,DebyePeaks['debyeTerms'][i][::2]))) |
---|
3595 | debyeList += zip(debyeNames,DebyePeaks['debyeTerms'][i][1::2]) |
---|
3596 | debyeVary = [] |
---|
3597 | for item in debyeList: |
---|
3598 | if item[1]: |
---|
3599 | debyeVary.append(item[0]) |
---|
3600 | backDict.update(debyeDict) |
---|
3601 | backVary += debyeVary |
---|
3602 | peakDict = {} |
---|
3603 | peakList = [] |
---|
3604 | for i in range(DebyePeaks['nPeaks']): |
---|
3605 | peakNames = [':'+str(hId)+':BkPkpos;'+str(i),':'+str(hId)+ \ |
---|
3606 | ':BkPkint;'+str(i),':'+str(hId)+':BkPksig;'+str(i),':'+str(hId)+':BkPkgam;'+str(i)] |
---|
3607 | peakDict.update(dict(zip(peakNames,DebyePeaks['peaksList'][i][::2]))) |
---|
3608 | peakList += zip(peakNames,DebyePeaks['peaksList'][i][1::2]) |
---|
3609 | peakVary = [] |
---|
3610 | for item in peakList: |
---|
3611 | if item[1]: |
---|
3612 | peakVary.append(item[0]) |
---|
3613 | backDict.update(peakDict) |
---|
3614 | backVary += peakVary |
---|
3615 | if 'background PWDR' in Background[1]: |
---|
3616 | backDict[':'+str(hId)+':Back File'] = Background[1]['background PWDR'][0] |
---|
3617 | backDict[':'+str(hId)+':BF mult'] = Background[1]['background PWDR'][1] |
---|
3618 | try: |
---|
3619 | if Background[1]['background PWDR'][0] and Background[1]['background PWDR'][2]: |
---|
3620 | backVary.append(':'+str(hId)+':BF mult') |
---|
3621 | except IndexError: # old version without refine flag |
---|
3622 | pass |
---|
3623 | return bakType,backDict,backVary |
---|
3624 | |
---|
3625 | def GetInstParms(hId,Inst): |
---|
3626 | #patch |
---|
3627 | dataType = Inst['Type'][0] |
---|
3628 | if 'Z' not in Inst and 'E' not in dataType: |
---|
3629 | Inst['Z'] = [0.0,0.0,False] |
---|
3630 | instDict = {} |
---|
3631 | insVary = [] |
---|
3632 | pfx = ':'+str(hId)+':' |
---|
3633 | insKeys = list(Inst.keys()) |
---|
3634 | insKeys.sort() |
---|
3635 | for item in insKeys: |
---|
3636 | if 'Azimuth' in item: |
---|
3637 | continue |
---|
3638 | insName = pfx+item |
---|
3639 | instDict[insName] = Inst[item][1] |
---|
3640 | if len(Inst[item]) > 2 and Inst[item][2]: |
---|
3641 | insVary.append(insName) |
---|
3642 | if 'C' in dataType: |
---|
3643 | instDict[pfx+'SH/L'] = max(instDict[pfx+'SH/L'],0.0005) |
---|
3644 | elif 'T' in dataType: #trap zero alp, bet coeff. |
---|
3645 | if not instDict[pfx+'alpha']: |
---|
3646 | instDict[pfx+'alpha'] = 1.0 |
---|
3647 | if not instDict[pfx+'beta-0'] and not instDict[pfx+'beta-1']: |
---|
3648 | instDict[pfx+'beta-1'] = 1.0 |
---|
3649 | elif 'B' in dataType: #trap zero alp, bet coeff. |
---|
3650 | if not instDict[pfx+'alpha-0'] and not instDict[pfx+'alpha-1']: |
---|
3651 | instDict[pfx+'alpha-1'] = 1.0 |
---|
3652 | if not instDict[pfx+'beta-0'] and not instDict[pfx+'beta-1']: |
---|
3653 | instDict[pfx+'beta-1'] = 1.0 |
---|
3654 | elif 'E' in dataType: |
---|
3655 | pass |
---|
3656 | return dataType,instDict,insVary |
---|
3657 | |
---|
3658 | def GetSampleParms(hId,Sample): |
---|
3659 | sampVary = [] |
---|
3660 | hfx = ':'+str(hId)+':' |
---|
3661 | sampDict = {hfx+'Gonio. radius':Sample['Gonio. radius'],hfx+'Omega':Sample['Omega'], |
---|
3662 | hfx+'Chi':Sample['Chi'],hfx+'Phi':Sample['Phi'],hfx+'Azimuth':Sample['Azimuth']} |
---|
3663 | for key in ('Temperature','Pressure','FreePrm1','FreePrm2','FreePrm3'): |
---|
3664 | if key in Sample: |
---|
3665 | sampDict[hfx+key] = Sample[key] |
---|
3666 | Type = Sample['Type'] |
---|
3667 | if 'Bragg' in Type: #Bragg-Brentano |
---|
3668 | for item in ['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']: |
---|
3669 | sampDict[hfx+item] = Sample[item][0] |
---|
3670 | if Sample[item][1]: |
---|
3671 | sampVary.append(hfx+item) |
---|
3672 | elif 'Debye' in Type: #Debye-Scherrer |
---|
3673 | for item in ['Scale','Absorption','DisplaceX','DisplaceY']: |
---|
3674 | sampDict[hfx+item] = Sample[item][0] |
---|
3675 | if Sample[item][1]: |
---|
3676 | sampVary.append(hfx+item) |
---|
3677 | return Type,sampDict,sampVary |
---|
3678 | |
---|
3679 | def PrintBackground(Background): |
---|
3680 | Back = Background[0] |
---|
3681 | DebyePeaks = Background[1] |
---|
3682 | pFile.write('\n Background function: %s Refine? %s\n'%(Back[0],Back[1])) |
---|
3683 | line = ' Coefficients: ' |
---|
3684 | for i,back in enumerate(Back[3:]): |
---|
3685 | line += '%10.3f'%(back) |
---|
3686 | if i and not i%10: |
---|
3687 | line += '\n'+15*' ' |
---|
3688 | pFile.write(line+'\n') |
---|
3689 | if DebyePeaks['nDebye']: |
---|
3690 | pFile.write('\n Debye diffuse scattering coefficients\n') |
---|
3691 | parms = ['DebyeA','DebyeR','DebyeU'] |
---|
3692 | line = ' names : ' |
---|
3693 | for parm in parms: |
---|
3694 | line += '%8s refine?'%(parm) |
---|
3695 | pFile.write(line+'\n') |
---|
3696 | for j,term in enumerate(DebyePeaks['debyeTerms']): |
---|
3697 | line = ' term'+'%2d'%(j)+':' |
---|
3698 | for i in range(3): |
---|
3699 | line += '%10.3f %5s'%(term[2*i],bool(term[2*i+1])) |
---|
3700 | pFile.write(line+'\n') |
---|
3701 | if DebyePeaks['nPeaks']: |
---|
3702 | pFile.write('\n Single peak coefficients\n') |
---|
3703 | parms = ['BkPkpos','BkPkint','BkPksig','BkPkgam'] |
---|
3704 | line = ' names : ' |
---|
3705 | for parm in parms: |
---|
3706 | line += '%8s refine?'%(parm) |
---|
3707 | pFile.write(line+'\n') |
---|
3708 | for j,term in enumerate(DebyePeaks['peaksList']): |
---|
3709 | line = ' peak'+'%2d'%(j)+':' |
---|
3710 | for i in range(4): |
---|
3711 | line += '%12.3f %5s'%(term[2*i],bool(term[2*i+1])) |
---|
3712 | pFile.write(line+'\n') |
---|
3713 | if 'background PWDR' in DebyePeaks: |
---|
3714 | try: |
---|
3715 | pFile.write(' Fixed background file: %s; mult: %.3f Refine? %s\n'%(DebyePeaks['background PWDR'][0], |
---|
3716 | DebyePeaks['background PWDR'][1],DebyePeaks['background PWDR'][2])) |
---|
3717 | except IndexError: #old version without refine flag |
---|
3718 | pass |
---|
3719 | |
---|
3720 | def PrintInstParms(Inst): |
---|
3721 | pFile.write('\n Instrument Parameters:\n') |
---|
3722 | insKeys = list(Inst.keys()) |
---|
3723 | insKeys.sort() |
---|
3724 | iBeg = 0 |
---|
3725 | Ok = True |
---|
3726 | while Ok: |
---|
3727 | ptlbls = ' name :' |
---|
3728 | ptstr = ' value :' |
---|
3729 | varstr = ' refine:' |
---|
3730 | iFin = min(iBeg+9,len(insKeys)) |
---|
3731 | for item in insKeys[iBeg:iFin]: |
---|
3732 | if item not in ['Type','Source','Bank']: |
---|
3733 | ptlbls += '%12s' % (item) |
---|
3734 | ptstr += '%12.6f' % (Inst[item][1]) |
---|
3735 | if item in ['Lam1','Lam2','Azimuth','fltPath']: |
---|
3736 | varstr += 12*' ' |
---|
3737 | else: |
---|
3738 | varstr += '%12s' % (str(bool(Inst[item][2]))) |
---|
3739 | pFile.write(ptlbls+'\n') |
---|
3740 | pFile.write(ptstr+'\n') |
---|
3741 | pFile.write(varstr+'\n') |
---|
3742 | iBeg = iFin |
---|
3743 | if iBeg == len(insKeys): |
---|
3744 | Ok = False |
---|
3745 | else: |
---|
3746 | pFile.write('\n') |
---|
3747 | |
---|
3748 | def PrintSampleParms(Sample): |
---|
3749 | pFile.write('\n Sample Parameters:\n') |
---|
3750 | pFile.write(' Goniometer omega = %.2f, chi = %.2f, phi = %.2f\n'% |
---|
3751 | (Sample['Omega'],Sample['Chi'],Sample['Phi'])) |
---|
3752 | ptlbls = ' name :' |
---|
3753 | ptstr = ' value :' |
---|
3754 | varstr = ' refine:' |
---|
3755 | if 'Bragg' in Sample['Type']: |
---|
3756 | for item in ['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']: |
---|
3757 | ptlbls += '%14s'%(item) |
---|
3758 | ptstr += '%14.4f'%(Sample[item][0]) |
---|
3759 | varstr += '%14s'%(str(bool(Sample[item][1]))) |
---|
3760 | |
---|
3761 | elif 'Debye' in Type: #Debye-Scherrer |
---|
3762 | for item in ['Scale','Absorption','DisplaceX','DisplaceY']: |
---|
3763 | ptlbls += '%14s'%(item) |
---|
3764 | ptstr += '%14.4f'%(Sample[item][0]) |
---|
3765 | varstr += '%14s'%(str(bool(Sample[item][1]))) |
---|
3766 | |
---|
3767 | pFile.write(ptlbls+'\n') |
---|
3768 | pFile.write(ptstr+'\n') |
---|
3769 | pFile.write(varstr+'\n') |
---|
3770 | |
---|
3771 | histDict = {} |
---|
3772 | histVary = [] |
---|
3773 | controlDict = {} |
---|
3774 | histoList = list(Histograms.keys()) |
---|
3775 | histoList.sort() |
---|
3776 | for histogram in histoList: |
---|
3777 | if 'PWDR' in histogram: |
---|
3778 | Histogram = Histograms[histogram] |
---|
3779 | hId = Histogram['hId'] |
---|
3780 | pfx = ':'+str(hId)+':' |
---|
3781 | controlDict[pfx+'wtFactor'] = Histogram['wtFactor'] |
---|
3782 | controlDict[pfx+'Limits'] = Histogram['Limits'][1] |
---|
3783 | controlDict[pfx+'Exclude'] = Histogram['Limits'][2:] |
---|
3784 | for excl in controlDict[pfx+'Exclude']: |
---|
3785 | Histogram['Data'][0] = ma.masked_inside(Histogram['Data'][0],excl[0],excl[1]) |
---|
3786 | if controlDict[pfx+'Exclude']: |
---|
3787 | ma.mask_rows(Histogram['Data']) |
---|
3788 | Background = Histogram['Background'] |
---|
3789 | Type,bakDict,bakVary = GetBackgroundParms(hId,Background) |
---|
3790 | controlDict[pfx+'bakType'] = Type |
---|
3791 | histDict.update(bakDict) |
---|
3792 | histVary += bakVary |
---|
3793 | |
---|
3794 | Inst = Histogram['Instrument Parameters'] #TODO ? ignores tabulated alp,bet & delt for TOF |
---|
3795 | if 'T' in Type and len(Inst[1]): #patch - back-to-back exponential contribution to TOF line shape is removed |
---|
3796 | G2fil.G2Print ('Warning: tabulated profile coefficients are ignored') |
---|
3797 | Type,instDict,insVary = GetInstParms(hId,Inst[0]) |
---|
3798 | controlDict[pfx+'histType'] = Type |
---|
3799 | if 'XC' in Type or 'XB' in Type: |
---|
3800 | if pfx+'Lam1' in instDict: |
---|
3801 | controlDict[pfx+'keV'] = G2mth.wavekE(instDict[pfx+'Lam1']) |
---|
3802 | else: |
---|
3803 | controlDict[pfx+'keV'] = G2mth.wavekE(instDict[pfx+'Lam']) |
---|
3804 | histDict.update(instDict) |
---|
3805 | histVary += insVary |
---|
3806 | |
---|
3807 | Sample = Histogram['Sample Parameters'] |
---|
3808 | Type,sampDict,sampVary = GetSampleParms(hId,Sample) |
---|
3809 | controlDict[pfx+'instType'] = Type |
---|
3810 | histDict.update(sampDict) |
---|
3811 | histVary += sampVary |
---|
3812 | |
---|
3813 | |
---|
3814 | if Print: |
---|
3815 | pFile.write('\n Histogram: %s histogram Id: %d\n'%(histogram,hId)) |
---|
3816 | pFile.write(135*'='+'\n') |
---|
3817 | Units = {'C':' deg','T':' msec','B':' deg','E':'keV'} |
---|
3818 | units = Units[controlDict[pfx+'histType'][2]] |
---|
3819 | Limits = controlDict[pfx+'Limits'] |
---|
3820 | pFile.write(' Instrument type: %s\n'%Sample['Type']) |
---|
3821 | pFile.write(' Histogram limits: %8.2f%s to %8.2f%s\n'%(Limits[0],units,Limits[1],units)) |
---|
3822 | if len(controlDict[pfx+'Exclude']): |
---|
3823 | excls = controlDict[pfx+'Exclude'] |
---|
3824 | for excl in excls: |
---|
3825 | pFile.write(' Excluded region: %8.2f%s to %8.2f%s\n'%(excl[0],units,excl[1],units)) |
---|
3826 | PrintSampleParms(Sample) |
---|
3827 | PrintInstParms(Inst[0]) |
---|
3828 | PrintBackground(Background) |
---|
3829 | elif 'HKLF' in histogram: |
---|
3830 | Histogram = Histograms[histogram] |
---|
3831 | hId = Histogram['hId'] |
---|
3832 | pfx = ':'+str(hId)+':' |
---|
3833 | controlDict[pfx+'wtFactor'] = Histogram['wtFactor'] |
---|
3834 | Inst = Histogram['Instrument Parameters'][0] |
---|
3835 | controlDict[pfx+'histType'] = Inst['Type'][1] |
---|
3836 | if 'X' in Inst['Type'][1]: |
---|
3837 | histDict[pfx+'Lam'] = Inst['Lam'][1] |
---|
3838 | controlDict[pfx+'keV'] = G2mth.wavekE(histDict[pfx+'Lam']) |
---|
3839 | elif 'SEC' in Inst['Type'][1]: |
---|
3840 | histDict[pfx+'Lam'] = Inst['Lam'][1] |
---|
3841 | elif 'NC' in Inst['Type'][1] or 'NB' in Inst['Type'][1]: |
---|
3842 | histDict[pfx+'Lam'] = Inst['Lam'][1] |
---|
3843 | return histVary,histDict,controlDict |
---|
3844 | |
---|
3845 | def SetHistogramData(parmDict,sigDict,Histograms,calcControls,Print=True,pFile=None,seq=False): |
---|
3846 | 'Shows histogram data after a refinement' |
---|
3847 | |
---|
3848 | def SetBackgroundParms(pfx,Background,parmDict,sigDict): |
---|
3849 | Back = Background[0] |
---|
3850 | DebyePeaks = Background[1] |
---|
3851 | lenBack = len(Back[3:]) |
---|
3852 | backSig = [0 for i in range(lenBack+3*DebyePeaks['nDebye']+4*DebyePeaks['nPeaks'])] |
---|
3853 | for i in range(lenBack): |
---|
3854 | Back[3+i] = parmDict[pfx+'Back;'+str(i)] |
---|
3855 | if pfx+'Back;'+str(i) in sigDict: |
---|
3856 | backSig[i] = sigDict[pfx+'Back;'+str(i)] |
---|
3857 | if DebyePeaks['nDebye']: |
---|
3858 | for i in range(DebyePeaks['nDebye']): |
---|
3859 | names = [pfx+'DebyeA;'+str(i),pfx+'DebyeR;'+str(i),pfx+'DebyeU;'+str(i)] |
---|
3860 | for j,name in enumerate(names): |
---|
3861 | DebyePeaks['debyeTerms'][i][2*j] = parmDict[name] |
---|
3862 | if name in sigDict: |
---|
3863 | backSig[lenBack+3*i+j] = sigDict[name] |
---|
3864 | if DebyePeaks['nPeaks']: |
---|
3865 | for i in range(DebyePeaks['nPeaks']): |
---|
3866 | names = [pfx+'BkPkpos;'+str(i),pfx+'BkPkint;'+str(i), |
---|
3867 | pfx+'BkPksig;'+str(i),pfx+'BkPkgam;'+str(i)] |
---|
3868 | for j,name in enumerate(names): |
---|
3869 | DebyePeaks['peaksList'][i][2*j] = parmDict[name] |
---|
3870 | if name in sigDict: |
---|
3871 | backSig[lenBack+3*DebyePeaks['nDebye']+4*i+j] = sigDict[name] |
---|
3872 | if pfx+'BF mult' in sigDict: |
---|
3873 | DebyePeaks['background PWDR'][1] = parmDict[pfx+'BF mult'] |
---|
3874 | backSig.append(sigDict[pfx+'BF mult']) |
---|
3875 | |
---|
3876 | return backSig |
---|
3877 | |
---|
3878 | def SetInstParms(pfx,Inst,parmDict,sigDict): |
---|
3879 | instSig = {} |
---|
3880 | insKeys = list(Inst.keys()) |
---|
3881 | insKeys.sort() |
---|
3882 | for item in insKeys: |
---|
3883 | insName = pfx+item |
---|
3884 | Inst[item][1] = parmDict[insName] |
---|
3885 | if insName in sigDict: |
---|
3886 | instSig[item] = sigDict[insName] |
---|
3887 | else: |
---|
3888 | instSig[item] = 0 |
---|
3889 | return instSig |
---|
3890 | |
---|
3891 | def SetSampleParms(pfx,Sample,parmDict,sigDict): |
---|
3892 | if 'Bragg' in Sample['Type']: #Bragg-Brentano |
---|
3893 | sampSig = [0 for i in range(5)] |
---|
3894 | for i,item in enumerate(['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']): |
---|
3895 | Sample[item][0] = parmDict[pfx+item] |
---|
3896 | if pfx+item in sigDict: |
---|
3897 | sampSig[i] = sigDict[pfx+item] |
---|
3898 | elif 'Debye' in Sample['Type']: #Debye-Scherrer |
---|
3899 | sampSig = [0 for i in range(4)] |
---|
3900 | for i,item in enumerate(['Scale','Absorption','DisplaceX','DisplaceY']): |
---|
3901 | Sample[item][0] = parmDict[pfx+item] |
---|
3902 | if pfx+item in sigDict: |
---|
3903 | sampSig[i] = sigDict[pfx+item] |
---|
3904 | return sampSig |
---|
3905 | |
---|
3906 | def PrintBackgroundSig(Background,backSig): |
---|
3907 | Back = Background[0] |
---|
3908 | DebyePeaks = Background[1] |
---|
3909 | valstr = ' value : ' |
---|
3910 | sigstr = ' sig : ' |
---|
3911 | refine = False |
---|
3912 | for i,back in enumerate(Back[3:]): |
---|
3913 | valstr += '%10.4g'%(back) |
---|
3914 | if Back[1]: |
---|
3915 | refine = True |
---|
3916 | sigstr += '%10.4g'%(backSig[i]) |
---|
3917 | else: |
---|
3918 | sigstr += 10*' ' |
---|
3919 | if refine: |
---|
3920 | pFile.write('\n Background function: %s\n'%Back[0]) |
---|
3921 | pFile.write(valstr+'\n') |
---|
3922 | pFile.write(sigstr+'\n') |
---|
3923 | if DebyePeaks['nDebye']: |
---|
3924 | ifAny = False |
---|
3925 | ptfmt = "%12.3f" |
---|
3926 | names = ' names :' |
---|
3927 | ptstr = ' values:' |
---|
3928 | sigstr = ' esds :' |
---|
3929 | for item in sigDict: |
---|
3930 | if 'Debye' in item: |
---|
3931 | ifAny = True |
---|
3932 | names += '%12s'%(item) |
---|
3933 | ptstr += ptfmt%(parmDict[item]) |
---|
3934 | sigstr += ptfmt%(sigDict[item]) |
---|
3935 | if ifAny: |
---|
3936 | pFile.write('\n Debye diffuse scattering coefficients\n') |
---|
3937 | pFile.write(names+'\n') |
---|
3938 | pFile.write(ptstr+'\n') |
---|
3939 | pFile.write(sigstr+'\n') |
---|
3940 | if DebyePeaks['nPeaks']: |
---|
3941 | pFile.write('\n Single peak coefficients:\n') |
---|
3942 | parms = ['BkPkpos','BkPkint','BkPksig','BkPkgam'] |
---|
3943 | line = ' peak no. ' |
---|
3944 | for parm in parms: |
---|
3945 | line += '%14s%12s'%(parm.center(14),'esd'.center(12)) |
---|
3946 | pFile.write(line+'\n') |
---|
3947 | for ip in range(DebyePeaks['nPeaks']): |
---|
3948 | ptstr = ' %4d '%(ip) |
---|
3949 | for parm in parms: |
---|
3950 | fmt = '%14.3f' |
---|
3951 | efmt = '%12.3f' |
---|
3952 | if parm == 'BkPkpos': |
---|
3953 | fmt = '%14.4f' |
---|
3954 | efmt = '%12.4f' |
---|
3955 | name = pfx+parm+';%d'%(ip) |
---|
3956 | ptstr += fmt%(parmDict[ |
---|