source: trunk/GSASIIstrIO.py @ 3054

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

fixes to MCSA gui routines; trap <3 atom RBs
fix missing AtLookup? error in LeBail? refinements

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