source: trunk/GSASIIstrIO.py @ 3340

Last change on this file since 3340 was 3340, checked in by vondreele, 4 years ago

make testDeriv py3 compatible
fix problem with Reload draw atoms - mag moments
add |Mag| to lst output of magnetic moments
correct errors in mag moment structure factor & derivatives

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