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