source: trunk/GSASIIstrIO.py @ 2473

Last change on this file since 2473 was 2473, checked in by vondreele, 5 years ago

work on magnetic structures - import from EXP, plotting, LS refine I/O, mag. form factors, etc.

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