source: trunk/GSASIIstrIO.py @ 2474

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

fix spin inversion plot & some refine list issues

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 143.5 KB
Line 
1# -*- coding: utf-8 -*-
2########### SVN repository information ###################
3# $Date: 2016-09-20 20:35:05 +0000 (Tue, 20 Sep 2016) $
4# $Author: vondreele $
5# $Revision: 2474 $
6# $URL: trunk/GSASIIstrIO.py $
7# $Id: GSASIIstrIO.py 2474 2016-09-20 20:35:05Z 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: 2474 $")
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        AtInfo = dict(zip(General['AtomTypes'],General['Lande g']))
941        print >>pFile,'\n Magnetic moments:'
942        line = '   name    type  refine?  Mx        My        Mz    '
943        print >>pFile,line
944        print >>pFile,135*'-'
945        for i,at in enumerate(Atoms):
946            if AtInfo[at[ct]]:
947                line = '%7s'%(at[ct-1])+'%7s'%(at[ct])+'%7s'%(at[ct+1])+'%10.5f'%(at[cmx])+'%10.5f'%(at[cmx+1])+ \
948                    '%10.5f'%(at[cmx+2])
949                print >>pFile,line
950       
951               
952    def PrintWaves(General,Atoms):
953        cx,ct,cs,cia = General['AtomPtrs']
954        print >>pFile,'\n Modulation waves'
955        names = {'Sfrac':['Fsin','Fcos','Fzero','Fwid'],'Spos':['Xsin','Ysin','Zsin','Xcos','Ycos','Zcos','Tmin','Tmax','Xmax','Ymax','Zmax'],
956            'Sadp':['U11sin','U22sin','U33sin','U12sin','U13sin','U23sin','U11cos','U22cos',
957            'U33cos','U12cos','U13cos','U23cos'],'Smag':['MXsin','MYsin','MZsin','MXcos','MYcos','MZcos']}
958        print >>pFile,135*'-'
959        for i,at in enumerate(Atoms):
960            AtomSS = at[-1]['SS1']
961            for Stype in ['Sfrac','Spos','Sadp','Smag']:
962                Waves = AtomSS[Stype]
963                if len(Waves):
964                    print >>pFile,' atom: %s, site sym: %s, type: %s wave type: %s:'    \
965                        %(at[ct-1],at[cs],Stype,AtomSS['waveType'])
966                for iw,wave in enumerate(Waves):                   
967                    line = ''
968                    if AtomSS['waveType'] in ['Block','ZigZag'] and Stype == 'Spos' and not iw:
969                        for item in names[Stype][6:]:
970                            line += '%8s '%(item)                       
971                    else:
972                        if Stype == 'Spos':
973                            for item in names[Stype][:6]:
974                                line += '%8s '%(item)
975                        else:
976                            for item in names[Stype]:
977                                line += '%8s '%(item)
978                    print >>pFile,line
979                    line = ''
980                    for item in wave[0]:
981                        line += '%8.4f '%(item)
982                    line += ' Refine? '+str(wave[1])
983                    print >>pFile,line
984       
985    def PrintTexture(textureData):
986        topstr = '\n Spherical harmonics texture: Order:' + \
987            str(textureData['Order'])
988        if textureData['Order']:
989            print >>pFile,topstr+' Refine? '+str(textureData['SH Coeff'][0])
990        else:
991            print >>pFile,topstr
992            return
993        names = ['omega','chi','phi']
994        line = '\n'
995        for name in names:
996            line += ' SH '+name+':'+'%12.4f'%(textureData['Sample '+name][1])+' Refine? '+str(textureData['Sample '+name][0])
997        print >>pFile,line
998        print >>pFile,'\n Texture coefficients:'
999        SHcoeff = textureData['SH Coeff'][1]
1000        SHkeys = SHcoeff.keys()
1001        nCoeff = len(SHcoeff)
1002        nBlock = nCoeff/10+1
1003        iBeg = 0
1004        iFin = min(iBeg+10,nCoeff)
1005        for block in range(nBlock):
1006            ptlbls = ' names :'
1007            ptstr =  ' values:'
1008            for item in SHkeys[iBeg:iFin]:
1009                ptlbls += '%12s'%(item)
1010                ptstr += '%12.4f'%(SHcoeff[item]) 
1011            print >>pFile,ptlbls
1012            print >>pFile,ptstr
1013            iBeg += 10
1014            iFin = min(iBeg+10,nCoeff)
1015       
1016    def MakeRBParms(rbKey,phaseVary,phaseDict):
1017        rbid = str(rbids.index(RB['RBId']))
1018        pfxRB = pfx+'RB'+rbKey+'P'
1019        pstr = ['x','y','z']
1020        ostr = ['a','i','j','k']
1021        for i in range(3):
1022            name = pfxRB+pstr[i]+':'+str(iRB)+':'+rbid
1023            phaseDict[name] = RB['Orig'][0][i]
1024            if RB['Orig'][1]:
1025                phaseVary += [name,]
1026        pfxRB = pfx+'RB'+rbKey+'O'
1027        for i in range(4):
1028            name = pfxRB+ostr[i]+':'+str(iRB)+':'+rbid
1029            phaseDict[name] = RB['Orient'][0][i]
1030            if RB['Orient'][1] == 'AV' and i:
1031                phaseVary += [name,]
1032            elif RB['Orient'][1] == 'A' and not i:
1033                phaseVary += [name,]
1034           
1035    def MakeRBThermals(rbKey,phaseVary,phaseDict):
1036        rbid = str(rbids.index(RB['RBId']))
1037        tlstr = ['11','22','33','12','13','23']
1038        sstr = ['12','13','21','23','31','32','AA','BB']
1039        if 'T' in RB['ThermalMotion'][0]:
1040            pfxRB = pfx+'RB'+rbKey+'T'
1041            for i in range(6):
1042                name = pfxRB+tlstr[i]+':'+str(iRB)+':'+rbid
1043                phaseDict[name] = RB['ThermalMotion'][1][i]
1044                if RB['ThermalMotion'][2][i]:
1045                    phaseVary += [name,]
1046        if 'L' in RB['ThermalMotion'][0]:
1047            pfxRB = pfx+'RB'+rbKey+'L'
1048            for i in range(6):
1049                name = pfxRB+tlstr[i]+':'+str(iRB)+':'+rbid
1050                phaseDict[name] = RB['ThermalMotion'][1][i+6]
1051                if RB['ThermalMotion'][2][i+6]:
1052                    phaseVary += [name,]
1053        if 'S' in RB['ThermalMotion'][0]:
1054            pfxRB = pfx+'RB'+rbKey+'S'
1055            for i in range(8):
1056                name = pfxRB+sstr[i]+':'+str(iRB)+':'+rbid
1057                phaseDict[name] = RB['ThermalMotion'][1][i+12]
1058                if RB['ThermalMotion'][2][i+12]:
1059                    phaseVary += [name,]
1060        if 'U' in RB['ThermalMotion'][0]:
1061            name = pfx+'RB'+rbKey+'U:'+str(iRB)+':'+rbid
1062            phaseDict[name] = RB['ThermalMotion'][1][0]
1063            if RB['ThermalMotion'][2][0]:
1064                phaseVary += [name,]
1065               
1066    def MakeRBTorsions(rbKey,phaseVary,phaseDict):
1067        rbid = str(rbids.index(RB['RBId']))
1068        pfxRB = pfx+'RB'+rbKey+'Tr;'
1069        for i,torsion in enumerate(RB['Torsions']):
1070            name = pfxRB+str(i)+':'+str(iRB)+':'+rbid
1071            phaseDict[name] = torsion[0]
1072            if torsion[1]:
1073                phaseVary += [name,]
1074                   
1075    if Print:
1076        print  >>pFile,'\n Phases:'
1077    phaseVary = []
1078    phaseDict = {}
1079    phaseConstr = {}
1080    pawleyLookup = {}
1081    FFtables = {}                   #scattering factors - xrays
1082    MFtables = {}                   #Mag. form factors
1083    BLtables = {}                   # neutrons
1084    Natoms = {}
1085    AtMults = {}
1086    AtIA = {}
1087    maxSSwave = {}
1088    shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
1089    SamSym = dict(zip(shModels,['0','-1','2/m','mmm']))
1090    atomIndx = {}
1091    for name in PhaseData:
1092        General = PhaseData[name]['General']
1093        pId = PhaseData[name]['pId']
1094        pfx = str(pId)+'::'
1095        FFtable = G2el.GetFFtable(General['AtomTypes'])
1096        BLtable = G2el.GetBLtable(General)
1097        FFtables.update(FFtable)
1098        BLtables.update(BLtable)
1099        if General['Type'] == 'magnetic':
1100            MFtable = G2el.GetMFtable(General['AtomTypes'])
1101            MFtables.update(MFtable)
1102        Atoms = PhaseData[name]['Atoms']
1103        if Atoms and not General.get('doPawley'):
1104            cx,ct,cs,cia = General['AtomPtrs']
1105            AtLookup = G2mth.FillAtomLookUp(Atoms,cia+8)
1106        PawleyRef = PhaseData[name].get('Pawley ref',[])
1107        SGData = General['SGData']
1108        SGtext,SGtable = G2spc.SGPrint(SGData)
1109        cell = General['Cell']
1110        A = G2lat.cell2A(cell[1:7])
1111        phaseDict.update({pfx+'A0':A[0],pfx+'A1':A[1],pfx+'A2':A[2],
1112            pfx+'A3':A[3],pfx+'A4':A[4],pfx+'A5':A[5],pfx+'Vol':G2lat.calc_V(A)})
1113        if cell[0]:
1114            phaseVary += cellVary(pfx,SGData)       #also fills in symmetry required constraints
1115        SSGtext = []    #no superstructure
1116        im = 0
1117        if General.get('Modulated',False):
1118            im = 1      #refl offset
1119            Vec,vRef,maxH = General['SuperVec']
1120            phaseDict.update({pfx+'mV0':Vec[0],pfx+'mV1':Vec[1],pfx+'mV2':Vec[2]})
1121            SSGData = General['SSGData']
1122            SSGtext,SSGtable = G2spc.SSGPrint(SGData,SSGData)
1123            if vRef:
1124                phaseVary += modVary(pfx,SSGData)       
1125        resRBData = PhaseData[name]['RBModels'].get('Residue',[])
1126        if resRBData:
1127            rbids = rbIds['Residue']    #NB: used in the MakeRB routines
1128            for iRB,RB in enumerate(resRBData):
1129                MakeRBParms('R',phaseVary,phaseDict)
1130                MakeRBThermals('R',phaseVary,phaseDict)
1131                MakeRBTorsions('R',phaseVary,phaseDict)
1132       
1133        vecRBData = PhaseData[name]['RBModels'].get('Vector',[])
1134        if vecRBData:
1135            rbids = rbIds['Vector']    #NB: used in the MakeRB routines
1136            for iRB,RB in enumerate(vecRBData):
1137                MakeRBParms('V',phaseVary,phaseDict)
1138                MakeRBThermals('V',phaseVary,phaseDict)
1139                   
1140        Natoms[pfx] = 0
1141        maxSSwave[pfx] = {'Sfrac':0,'Spos':0,'Sadp':0,'Smag':0}
1142        if Atoms and not General.get('doPawley'):
1143            cx,ct,cs,cia = General['AtomPtrs']
1144            Natoms[pfx] = len(Atoms)
1145            for i,at in enumerate(Atoms):
1146                atomIndx[at[cia+8]] = [pfx,i]      #lookup table for restraints
1147                phaseDict.update({pfx+'Atype:'+str(i):at[ct],pfx+'Afrac:'+str(i):at[cx+3],pfx+'Amul:'+str(i):at[cs+1],
1148                    pfx+'Ax:'+str(i):at[cx],pfx+'Ay:'+str(i):at[cx+1],pfx+'Az:'+str(i):at[cx+2],
1149                    pfx+'dAx:'+str(i):0.,pfx+'dAy:'+str(i):0.,pfx+'dAz:'+str(i):0.,         #refined shifts for x,y,z
1150                    pfx+'AI/A:'+str(i):at[cia],})
1151                if at[cia] == 'I':
1152                    phaseDict[pfx+'AUiso:'+str(i)] = at[cia+1]
1153                else:
1154                    phaseDict.update({pfx+'AU11:'+str(i):at[cia+2],pfx+'AU22:'+str(i):at[cia+3],pfx+'AU33:'+str(i):at[cia+4],
1155                        pfx+'AU12:'+str(i):at[cia+5],pfx+'AU13:'+str(i):at[cia+6],pfx+'AU23:'+str(i):at[cia+7]})
1156                if General['Type'] == 'magnetic':
1157                    phaseDict.update({pfx+'AMx:'+str(i):at[cx+4],pfx+'AMy:'+str(i):at[cx+5],pfx+'AMz:'+str(i):at[cx+6]})
1158                if 'F' in at[ct+1]:
1159                    phaseVary.append(pfx+'Afrac:'+str(i))
1160                if 'X' in at[ct+1]:
1161                    try:    #patch for sytsym name changes
1162                        xId,xCoef = G2spc.GetCSxinel(at[cs])
1163                    except KeyError:
1164                        Sytsym = G2spc.SytSym(at[cx:cx+3],SGData)[0]
1165                        at[cs] = Sytsym
1166                        xId,xCoef = G2spc.GetCSxinel(at[cs])
1167                    xId,xCoef = G2spc.GetCSxinel(at[cs])
1168                    names = [pfx+'dAx:'+str(i),pfx+'dAy:'+str(i),pfx+'dAz:'+str(i)]
1169                    equivs = [[],[],[]]
1170                    for j in range(3):
1171                        if xId[j] > 0:                               
1172                            phaseVary.append(names[j])
1173                            equivs[xId[j]-1].append([names[j],xCoef[j]])
1174                    for equiv in equivs:
1175                        if len(equiv) > 1:
1176                            name = equiv[0][0]
1177                            coef = equiv[0][1]
1178                            for eqv in equiv[1:]:
1179                                eqv[1] /= coef
1180                                G2mv.StoreEquivalence(name,(eqv,))
1181                if 'U' in at[ct+1]:
1182                    if at[cia] == 'I':
1183                        phaseVary.append(pfx+'AUiso:'+str(i))
1184                    else:
1185                        try:    #patch for sytsym name changes
1186                            uId,uCoef = G2spc.GetCSuinel(at[cs])[:2]
1187                        except KeyError:
1188                            Sytsym = G2spc.SytSym(at[cx:cx+3],SGData)[0]
1189                            at[cs] = Sytsym
1190                            uId,uCoef = G2spc.GetCSuinel(at[cs])[:2]
1191                        uId,uCoef = G2spc.GetCSuinel(at[cs])[:2]
1192                        names = [pfx+'AU11:'+str(i),pfx+'AU22:'+str(i),pfx+'AU33:'+str(i),
1193                            pfx+'AU12:'+str(i),pfx+'AU13:'+str(i),pfx+'AU23:'+str(i)]
1194                        equivs = [[],[],[],[],[],[]]
1195                        for j in range(6):
1196                            if uId[j] > 0:                               
1197                                phaseVary.append(names[j])
1198                                equivs[uId[j]-1].append([names[j],uCoef[j]])
1199                        for equiv in equivs:
1200                            if len(equiv) > 1:
1201                                name = equiv[0][0]
1202                                coef = equiv[0][1]
1203                                for eqv in equiv[1:]:
1204                                    eqv[1] /= coef
1205                                G2mv.StoreEquivalence(name,equiv[1:])
1206                if 'M' in at[ct+1]:
1207                    pass    #magnetic moment vary here
1208                if General.get('Modulated',False):
1209                    AtomSS = at[-1]['SS1']
1210                    waveType = AtomSS['waveType']
1211                    phaseDict[pfx+'waveType:'+str(i)] = waveType
1212                    for Stype in ['Sfrac','Spos','Sadp','Smag']:
1213                        Waves = AtomSS[Stype]
1214                        nx = 0
1215                        for iw,wave in enumerate(Waves):
1216                            if not iw:
1217                                if waveType in ['ZigZag','Block']:
1218                                    nx = 1
1219                                CSI = G2spc.GetSSfxuinel(waveType,1,at[cx:cx+3],SGData,SSGData)
1220                            else:
1221                                CSI = G2spc.GetSSfxuinel('Fourier',iw+1-nx,at[cx:cx+3],SGData,SSGData)
1222                            uId,uCoef = CSI[Stype]
1223                            stiw = str(i)+':'+str(iw)
1224                            if Stype == 'Spos':
1225                                if waveType in ['ZigZag','Block',] and not iw:
1226                                    names = [pfx+'Tmin:'+stiw,pfx+'Tmax:'+stiw,pfx+'Xmax:'+stiw,pfx+'Ymax:'+stiw,pfx+'Zmax:'+stiw]
1227                                    equivs = [[],[], [],[],[]]
1228                                else:
1229                                    names = [pfx+'Xsin:'+stiw,pfx+'Ysin:'+stiw,pfx+'Zsin:'+stiw,
1230                                        pfx+'Xcos:'+stiw,pfx+'Ycos:'+stiw,pfx+'Zcos:'+stiw]
1231                                    equivs = [[],[],[], [],[],[]]
1232                            elif Stype == 'Sadp':
1233                                names = [pfx+'U11sin:'+stiw,pfx+'U22sin:'+stiw,pfx+'U33sin:'+stiw,
1234                                    pfx+'U12sin:'+stiw,pfx+'U13sin:'+stiw,pfx+'U23sin:'+stiw,
1235                                    pfx+'U11cos:'+stiw,pfx+'U22cos:'+stiw,pfx+'U33cos:'+stiw,
1236                                    pfx+'U12cos:'+stiw,pfx+'U13cos:'+stiw,pfx+'U23cos:'+stiw]
1237                                equivs = [[],[],[],[],[],[], [],[],[],[],[],[]]
1238                            elif Stype == 'Sfrac':
1239                                equivs = [[],[]]
1240                                if 'Crenel' in waveType and not iw:
1241                                    names = [pfx+'Fzero:'+stiw,pfx+'Fwid:'+stiw]
1242                                else:
1243                                    names = [pfx+'Fsin:'+stiw,pfx+'Fcos:'+stiw]
1244                            elif Stype == 'Smag':
1245                                equivs = [[],[],[], [],[],[]]
1246                                names = [pfx+'MXsin:'+stiw,pfx+'MYsin:'+stiw,pfx+'MZsin:'+stiw,
1247                                    pfx+'MXcos:'+stiw,pfx+'MYcos:'+stiw,pfx+'MZcos:'+stiw]
1248                            phaseDict.update(dict(zip(names,wave[0])))
1249                            if wave[1]: #what do we do here for multiple terms in modulation constraints?
1250                                for j in range(len(equivs)):
1251                                    if uId[j][0] > 0:                               
1252                                        phaseVary.append(names[j])
1253                                        equivs[uId[j][0]-1].append([names[j],uCoef[j][0]])
1254                                for equiv in equivs:
1255                                    if len(equiv) > 1:
1256                                        name = equiv[0][0]
1257                                        coef = equiv[0][1]
1258                                        for eqv in equiv[1:]:
1259                                            eqv[1] /= coef
1260                                        G2mv.StoreEquivalence(name,equiv[1:])
1261                            maxSSwave[pfx][Stype] = max(maxSSwave[pfx][Stype],iw+1)
1262            textureData = General['SH Texture']
1263            if textureData['Order'] and not seqRef:
1264                phaseDict[pfx+'SHorder'] = textureData['Order']
1265                phaseDict[pfx+'SHmodel'] = SamSym[textureData['Model']]
1266                for item in ['omega','chi','phi']:
1267                    phaseDict[pfx+'SH '+item] = textureData['Sample '+item][1]
1268                    if textureData['Sample '+item][0]:
1269                        phaseVary.append(pfx+'SH '+item)
1270                for item in textureData['SH Coeff'][1]:
1271                    phaseDict[pfx+item] = textureData['SH Coeff'][1][item]
1272                    if textureData['SH Coeff'][0]:
1273                        phaseVary.append(pfx+item)
1274               
1275            if Print:
1276                print >>pFile,'\n Phase name: ',General['Name']
1277                print >>pFile,135*'-'
1278                PrintFFtable(FFtable)
1279                PrintBLtable(BLtable)
1280                if General['Type'] == 'magnetic':
1281                    PrintMFtable(MFtable)
1282                print >>pFile,''
1283                #how do we print magnetic symmetry table? TBD
1284                if len(SSGtext):    #if superstructure
1285                    for line in SSGtext: print >>pFile,line
1286                    if len(SSGtable):
1287                        for item in SSGtable:
1288                            line = ' %s '%(item)
1289                            print >>pFile,line   
1290                    else:
1291                        print >>pFile,' ( 1)    %s'%(SSGtable[0])
1292                else:
1293                    for line in SGtext: print >>pFile,line
1294                    if len(SGtable):
1295                        for item in SGtable:
1296                            line = ' %s '%(item)
1297                            print >>pFile,line   
1298                    else:
1299                        print >>pFile,' ( 1)    %s'%(SGtable[0])
1300                PrintRBObjects(resRBData,vecRBData)
1301                PrintAtoms(General,Atoms)
1302                if General['Type'] == 'magnetic':
1303                    PrintMoments(General,Atoms)
1304                if General.get('Modulated',False):
1305                    PrintWaves(General,Atoms)
1306                print >>pFile,'\n Unit cell: a = %.5f'%(cell[1]),' b = %.5f'%(cell[2]),' c = %.5f'%(cell[3]), \
1307                    ' alpha = %.3f'%(cell[4]),' beta = %.3f'%(cell[5]),' gamma = %.3f'%(cell[6]), \
1308                    ' volume = %.3f'%(cell[7]),' Refine?',cell[0]
1309                if len(SSGtext):    #if superstructure
1310                    print >>pFile,'\n Modulation vector: mV0 = %.4f'%(Vec[0]),' mV1 = %.4f'%(Vec[1]),   \
1311                        ' mV2 = %.4f'%(Vec[2]),' max mod. index = %d'%(maxH),' Refine?',vRef
1312                if not seqRef:
1313                    PrintTexture(textureData)
1314                if name in RestraintDict:
1315                    PrintRestraints(cell[1:7],SGData,General['AtomPtrs'],Atoms,AtLookup,
1316                        textureData,RestraintDict[name],pFile)
1317                   
1318        elif PawleyRef:
1319            if Print:
1320                print >>pFile,'\n Phase name: ',General['Name']
1321                print >>pFile,135*'-'
1322                print >>pFile,''
1323                if len(SSGtext):    #if superstructure
1324                    for line in SSGtext: print >>pFile,line
1325                    if len(SSGtable):
1326                        for item in SSGtable:
1327                            line = ' %s '%(item)
1328                            print >>pFile,line   
1329                    else:
1330                        print >>pFile,' ( 1)    %s'%(SSGtable[0])
1331                else:
1332                    for line in SGtext: print >>pFile,line
1333                    if len(SGtable):
1334                        for item in SGtable:
1335                            line = ' %s '%(item)
1336                            print >>pFile,line   
1337                    else:
1338                        print >>pFile,' ( 1)    %s'%(SGtable[0])
1339                print >>pFile,'\n Unit cell: a = %.5f'%(cell[1]),' b = %.5f'%(cell[2]),' c = %.5f'%(cell[3]), \
1340                    ' alpha = %.3f'%(cell[4]),' beta = %.3f'%(cell[5]),' gamma = %.3f'%(cell[6]), \
1341                    ' volume = %.3f'%(cell[7]),' Refine?',cell[0]
1342                if len(SSGtext):    #if superstructure
1343                    print >>pFile,'\n Modulation vector: mV0 = %.4f'%(Vec[0]),' mV1 = %.4f'%(Vec[1]),   \
1344                        ' mV2 = %.4f'%(Vec[2]),' max mod. index = %d'%(maxH),' Refine?',vRef
1345            pawleyVary = []
1346            for i,refl in enumerate(PawleyRef):
1347                phaseDict[pfx+'PWLref:'+str(i)] = refl[6+im]
1348                if im:
1349                    pawleyLookup[pfx+'%d,%d,%d,%d'%(refl[0],refl[1],refl[2],refl[3])] = i
1350                else:
1351                    pawleyLookup[pfx+'%d,%d,%d'%(refl[0],refl[1],refl[2])] = i
1352                if refl[5+im]:
1353                    pawleyVary.append(pfx+'PWLref:'+str(i))
1354            GetPawleyConstr(SGData['SGLaue'],PawleyRef,im,pawleyVary)      #does G2mv.StoreEquivalence
1355            phaseVary += pawleyVary
1356               
1357    return Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables,maxSSwave
1358   
1359def cellFill(pfx,SGData,parmDict,sigDict): 
1360    '''Returns the filled-out reciprocal cell (A) terms and their uncertainties
1361    from the parameter and sig dictionaries.
1362
1363    :param str pfx: parameter prefix ("n::", where n is a phase number)
1364    :param dict SGdata: a symmetry object
1365    :param dict parmDict: a dictionary of parameters
1366    :param dict sigDict:  a dictionary of uncertainties on parameters
1367
1368    :returns: A,sigA where each is a list of six terms with the A terms
1369    '''
1370    if SGData['SGLaue'] in ['-1',]:
1371        A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
1372            parmDict[pfx+'A3'],parmDict[pfx+'A4'],parmDict[pfx+'A5']]
1373    elif SGData['SGLaue'] in ['2/m',]:
1374        if SGData['SGUniq'] == 'a':
1375            A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
1376                0,0,parmDict[pfx+'A5']]
1377        elif SGData['SGUniq'] == 'b':
1378            A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
1379                0,parmDict[pfx+'A4'],0]
1380        else:
1381            A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
1382                parmDict[pfx+'A3'],0,0]
1383    elif SGData['SGLaue'] in ['mmm',]:
1384        A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],0,0,0]
1385    elif SGData['SGLaue'] in ['4/m','4/mmm']:
1386        A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A2'],0,0,0]
1387    elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
1388        A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A2'],
1389            parmDict[pfx+'A0'],0,0]
1390    elif SGData['SGLaue'] in ['3R', '3mR']:
1391        A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A0'],
1392            parmDict[pfx+'A3'],parmDict[pfx+'A3'],parmDict[pfx+'A3']]
1393    elif SGData['SGLaue'] in ['m3m','m3']:
1394        A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A0'],0,0,0]
1395
1396    try:
1397        if SGData['SGLaue'] in ['-1',]:
1398            sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
1399                sigDict[pfx+'A3'],sigDict[pfx+'A4'],sigDict[pfx+'A5']]
1400        elif SGData['SGLaue'] in ['2/m',]:
1401            if SGData['SGUniq'] == 'a':
1402                sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
1403                    0,0,sigDict[pfx+'A5']]
1404            elif SGData['SGUniq'] == 'b':
1405                sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
1406                    0,sigDict[pfx+'A4'],0]
1407            else:
1408                sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
1409                    sigDict[pfx+'A3'],0,0]
1410        elif SGData['SGLaue'] in ['mmm',]:
1411            sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],0,0,0]
1412        elif SGData['SGLaue'] in ['4/m','4/mmm']:
1413            sigA = [sigDict[pfx+'A0'],0,sigDict[pfx+'A2'],0,0,0]
1414        elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
1415            sigA = [sigDict[pfx+'A0'],0,sigDict[pfx+'A2'],0,0,0]
1416        elif SGData['SGLaue'] in ['3R', '3mR']:
1417            sigA = [sigDict[pfx+'A0'],0,0,sigDict[pfx+'A3'],0,0]
1418        elif SGData['SGLaue'] in ['m3m','m3']:
1419            sigA = [sigDict[pfx+'A0'],0,0,0,0,0]
1420    except KeyError:
1421        sigA = [0,0,0,0,0,0]
1422    return A,sigA
1423       
1424def PrintRestraints(cell,SGData,AtPtrs,Atoms,AtLookup,textureData,phaseRest,pFile):
1425    'needs a doc string'
1426    if phaseRest:
1427        Amat = G2lat.cell2AB(cell)[0]
1428        cx,ct,cs = AtPtrs[:3]
1429        names = [['Bond','Bonds'],['Angle','Angles'],['Plane','Planes'],
1430            ['Chiral','Volumes'],['Torsion','Torsions'],['Rama','Ramas'],
1431            ['ChemComp','Sites'],['Texture','HKLs']]
1432        for name,rest in names:
1433            itemRest = phaseRest[name]
1434            if itemRest[rest] and itemRest['Use']:
1435                print >>pFile,'\n  %s %10.3f Use: %s'%(name+' restraint weight factor',itemRest['wtFactor'],str(itemRest['Use']))
1436                if name in ['Bond','Angle','Plane','Chiral']:
1437                    print >>pFile,'     calc       obs      sig   delt/sig  atoms(symOp): '
1438                    for indx,ops,obs,esd in itemRest[rest]:
1439                        try:
1440                            AtNames = G2mth.GetAtomItemsById(Atoms,AtLookup,indx,ct-1)
1441                            AtName = ''
1442                            for i,Aname in enumerate(AtNames):
1443                                AtName += Aname
1444                                if ops[i] == '1':
1445                                    AtName += '-'
1446                                else:
1447                                    AtName += '+('+ops[i]+')-'
1448                            XYZ = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cx,3))
1449                            XYZ = G2mth.getSyXYZ(XYZ,ops,SGData)
1450                            if name == 'Bond':
1451                                calc = G2mth.getRestDist(XYZ,Amat)
1452                            elif name == 'Angle':
1453                                calc = G2mth.getRestAngle(XYZ,Amat)
1454                            elif name == 'Plane':
1455                                calc = G2mth.getRestPlane(XYZ,Amat)
1456                            elif name == 'Chiral':
1457                                calc = G2mth.getRestChiral(XYZ,Amat)
1458                            print >>pFile,' %9.3f %9.3f %8.3f %8.3f   %s'%(calc,obs,esd,(obs-calc)/esd,AtName[:-1])
1459                        except KeyError:
1460                            del itemRest[rest]
1461                elif name in ['Torsion','Rama']:
1462                    print >>pFile,'  atoms(symOp)  calc  obs  sig  delt/sig  torsions: '
1463                    coeffDict = itemRest['Coeff']
1464                    for indx,ops,cofName,esd in itemRest[rest]:
1465                        AtNames = G2mth.GetAtomItemsById(Atoms,AtLookup,indx,ct-1)
1466                        AtName = ''
1467                        for i,Aname in enumerate(AtNames):
1468                            AtName += Aname+'+('+ops[i]+')-'
1469                        XYZ = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cx,3))
1470                        XYZ = G2mth.getSyXYZ(XYZ,ops,SGData)
1471                        if name == 'Torsion':
1472                            tor = G2mth.getRestTorsion(XYZ,Amat)
1473                            restr,calc = G2mth.calcTorsionEnergy(tor,coeffDict[cofName])
1474                            print >>pFile,' %8.3f %8.3f %.3f %8.3f %8.3f %s'%(calc,obs,esd,(obs-calc)/esd,tor,AtName[:-1])
1475                        else:
1476                            phi,psi = G2mth.getRestRama(XYZ,Amat)
1477                            restr,calc = G2mth.calcRamaEnergy(phi,psi,coeffDict[cofName])                               
1478                            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])
1479                elif name == 'ChemComp':
1480                    print >>pFile,'     atoms   mul*frac  factor     prod'
1481                    for indx,factors,obs,esd in itemRest[rest]:
1482                        try:
1483                            atoms = G2mth.GetAtomItemsById(Atoms,AtLookup,indx,ct-1)
1484                            mul = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cs+1))
1485                            frac = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cs-1))
1486                            mulfrac = mul*frac
1487                            calcs = mul*frac*factors
1488                            for iatm,[atom,mf,fr,clc] in enumerate(zip(atoms,mulfrac,factors,calcs)):
1489                                print >>pFile,' %10s %8.3f %8.3f %8.3f'%(atom,mf,fr,clc)
1490                            print >>pFile,' Sum:                   calc: %8.3f obs: %8.3f esd: %8.3f'%(np.sum(calcs),obs,esd)
1491                        except KeyError:
1492                            del itemRest[rest]
1493                elif name == 'Texture' and textureData['Order']:
1494                    Start = False
1495                    SHCoef = textureData['SH Coeff'][1]
1496                    shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
1497                    SamSym = dict(zip(shModels,['0','-1','2/m','mmm']))
1498                    print '    HKL  grid  neg esd  sum neg  num neg use unit?  unit esd  '
1499                    for hkl,grid,esd1,ifesd2,esd2 in itemRest[rest]:
1500                        PH = np.array(hkl)
1501                        phi,beta = G2lat.CrsAng(np.array(hkl),cell,SGData)
1502                        ODFln = G2lat.Flnh(Start,SHCoef,phi,beta,SGData)
1503                        R,P,Z = G2mth.getRestPolefig(ODFln,SamSym[textureData['Model']],grid)
1504                        Z = ma.masked_greater(Z,0.0)
1505                        num = ma.count(Z)
1506                        sum = 0
1507                        if num:
1508                            sum = np.sum(Z)
1509                        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)
1510       
1511def getCellEsd(pfx,SGData,A,covData):
1512    'needs a doc string'
1513    dpr = 180./np.pi
1514    rVsq = G2lat.calc_rVsq(A)
1515    G,g = G2lat.A2Gmat(A)       #get recip. & real metric tensors
1516    cell = np.array(G2lat.Gmat2cell(g))   #real cell
1517    cellst = np.array(G2lat.Gmat2cell(G)) #recip. cell
1518    scos = cosd(cellst[3:6])
1519    ssin = sind(cellst[3:6])
1520    scot = scos/ssin
1521    rcos = cosd(cell[3:6])
1522    rsin = sind(cell[3:6])
1523    rcot = rcos/rsin
1524    RMnames = [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3',pfx+'A4',pfx+'A5']
1525    varyList = covData['varyList']
1526    covMatrix = covData['covMatrix']
1527    vcov = G2mth.getVCov(RMnames,varyList,covMatrix)
1528    Ax = np.array(A)
1529    Ax[3:] /= 2.
1530    drVdA = np.array([
1531        Ax[1]*Ax[2]-Ax[5]**2,
1532        Ax[0]*Ax[2]-Ax[4]**2,
1533        Ax[0]*Ax[1]-Ax[3]**2,
1534        Ax[4]*Ax[5]-Ax[2]*Ax[3],
1535        Ax[3]*Ax[5]-Ax[1]*Ax[4],
1536        Ax[3]*Ax[4]-Ax[0]*Ax[5]])
1537    srcvlsq = np.inner(drVdA,np.inner(drVdA,vcov))
1538    Vol = 1/np.sqrt(rVsq)
1539    sigVol = Vol**3*np.sqrt(srcvlsq)/2.         #ok - checks with GSAS
1540    R123 = Ax[0]*Ax[1]*Ax[2]
1541    dsasdg = np.zeros((3,6))
1542    dadg = np.zeros((6,6))
1543    for i0 in range(3):         #0  1   2
1544        i1 = (i0+1)%3           #1  2   0
1545        i2 = (i1+1)%3           #2  0   1
1546        i3 = 5-i2               #3  5   4
1547        i4 = 5-i1               #4  3   5
1548        i5 = 5-i0               #5  4   3
1549        dsasdg[i0][i1] = 0.5*scot[i0]*scos[i0]/Ax[i1]
1550        dsasdg[i0][i2] = 0.5*scot[i0]*scos[i0]/Ax[i2]
1551        dsasdg[i0][i5] = -scot[i0]/np.sqrt(Ax[i1]*Ax[i2])
1552        denmsq = Ax[i0]*(R123-Ax[i1]*Ax[i4]**2-Ax[i2]*Ax[i3]**2+(Ax[i4]*Ax[i3])**2)
1553        denom = np.sqrt(denmsq)
1554        dadg[i5][i0] = -Ax[i5]/denom-rcos[i0]/denmsq*(R123-0.5*Ax[i1]*Ax[i4]**2-0.5*Ax[i2]*Ax[i3]**2)
1555        dadg[i5][i1] = -0.5*rcos[i0]/denmsq*(Ax[i0]**2*Ax[i2]-Ax[i0]*Ax[i4]**2)
1556        dadg[i5][i2] = -0.5*rcos[i0]/denmsq*(Ax[i0]**2*Ax[i1]-Ax[i0]*Ax[i3]**2)
1557        dadg[i5][i3] = Ax[i4]/denom+rcos[i0]/denmsq*(Ax[i0]*Ax[i2]*Ax[i3]-Ax[i3]*Ax[i4]**2)
1558        dadg[i5][i4] = Ax[i3]/denom+rcos[i0]/denmsq*(Ax[i0]*Ax[i1]*Ax[i4]-Ax[i3]**2*Ax[i4])
1559        dadg[i5][i5] = -Ax[i0]/denom
1560    for i0 in range(3):
1561        i1 = (i0+1)%3
1562        i2 = (i1+1)%3
1563        i3 = 5-i2
1564        for ij in range(6):
1565            dadg[i0][ij] = cell[i0]*(rcot[i2]*dadg[i3][ij]/rsin[i2]-dsasdg[i1][ij]/ssin[i1])
1566            if ij == i0:
1567                dadg[i0][ij] = dadg[i0][ij]-0.5*cell[i0]/Ax[i0]
1568            dadg[i3][ij] = -dadg[i3][ij]*rsin[2-i0]*dpr
1569    sigMat = np.inner(dadg,np.inner(dadg,vcov))
1570    var = np.diag(sigMat)
1571    CS = np.where(var>0.,np.sqrt(var),0.)
1572    return [CS[0],CS[1],CS[2],CS[5],CS[4],CS[3],sigVol]
1573   
1574def SetPhaseData(parmDict,sigDict,Phases,RBIds,covData,RestraintDict=None,pFile=None):
1575    'needs a doc string'
1576   
1577    def PrintAtomsAndSig(General,Atoms,atomsSig):
1578        print >>pFile,'\n Atoms:'
1579        line = '   name      x         y         z      frac   Uiso     U11     U22     U33     U12     U13     U23'
1580        if General['Type'] == 'macromolecular':
1581            line = ' res no residue chain '+line
1582        cx,ct,cs,cia = General['AtomPtrs']
1583        print >>pFile,line
1584        print >>pFile,135*'-'
1585        fmt = {0:'%7s',ct:'%7s',cx:'%10.5f',cx+1:'%10.5f',cx+2:'%10.5f',cx+3:'%8.3f',cia+1:'%8.5f',
1586            cia+2:'%8.5f',cia+3:'%8.5f',cia+4:'%8.5f',cia+5:'%8.5f',cia+6:'%8.5f',cia+7:'%8.5f'}
1587        noFXsig = {cx:[10*' ','%10s'],cx+1:[10*' ','%10s'],cx+2:[10*' ','%10s'],cx+3:[8*' ','%8s']}
1588        for atyp in General['AtomTypes']:       #zero composition
1589            General['NoAtoms'][atyp] = 0.
1590        for i,at in enumerate(Atoms):
1591            General['NoAtoms'][at[ct]] += at[cx+3]*float(at[cx+5])     #new composition
1592            if General['Type'] == 'macromolecular':
1593                name = ' %s %s %s %s:'%(at[0],at[1],at[2],at[3])
1594                valstr = ' values:          '
1595                sigstr = ' sig   :          '
1596            else:
1597                name = fmt[0]%(at[ct-1])+fmt[1]%(at[ct])+':'
1598                valstr = ' values:'
1599                sigstr = ' sig   :'
1600            for ind in range(cx,cx+4):
1601                sigind = str(i)+':'+str(ind)
1602                valstr += fmt[ind]%(at[ind])                   
1603                if sigind in atomsSig:
1604                    sigstr += fmt[ind]%(atomsSig[sigind])
1605                else:
1606                    sigstr += noFXsig[ind][1]%(noFXsig[ind][0])
1607            if at[cia] == 'I':
1608                valstr += fmt[cia+1]%(at[cia+1])
1609                if '%d:%d'%(i,cia+1) in atomsSig:
1610                    sigstr += fmt[cia+1]%(atomsSig['%d:%d'%(i,cia+1)])
1611                else:
1612                    sigstr += 8*' '
1613            else:
1614                valstr += 8*' '
1615                sigstr += 8*' '
1616                for ind in range(cia+2,cia+8):
1617                    sigind = str(i)+':'+str(ind)
1618                    valstr += fmt[ind]%(at[ind])
1619                    if sigind in atomsSig:                       
1620                        sigstr += fmt[ind]%(atomsSig[sigind])
1621                    else:
1622                        sigstr += 8*' '
1623            print >>pFile,name
1624            print >>pFile,valstr
1625            print >>pFile,sigstr
1626           
1627    def PrintMomentsAndSig(General,Atoms,atomsSig):
1628        print >>pFile,'\n Magnetic Moments:'    #add magnitude & angle, etc.? TBD
1629        line = '   name      Mx        My        Mz'
1630        cx,ct,cs,cia = General['AtomPtrs']
1631        cmx = cx+4
1632        AtInfo = dict(zip(General['AtomTypes'],General['Lande g']))
1633        print >>pFile,line
1634        print >>pFile,135*'-'
1635        fmt = {0:'%7s',ct:'%7s',cmx:'%10.5f',cmx+1:'%10.5f',cmx+2:'%10.5f'}
1636        noFXsig = {cmx:[10*' ','%10s'],cmx+1:[10*' ','%10s'],cmx+2:[10*' ','%10s']}
1637        for i,at in enumerate(Atoms):
1638            if AtInfo[at[ct]]:
1639                name = fmt[0]%(at[ct-1])+fmt[1]%(at[ct])+':'
1640                valstr = ' values:'
1641                sigstr = ' sig   :'
1642                for ind in range(cmx,cmx+3):
1643                    sigind = str(i)+':'+str(ind)
1644                    valstr += fmt[ind]%(at[ind])                   
1645                    if sigind in atomsSig:
1646                        sigstr += fmt[ind]%(atomsSig[sigind])
1647                    else:
1648                        sigstr += noFXsig[ind][1]%(noFXsig[ind][0])
1649                print >>pFile,name
1650                print >>pFile,valstr
1651                print >>pFile,sigstr
1652           
1653    def PrintWavesAndSig(General,Atoms,wavesSig):
1654        cx,ct,cs,cia = General['AtomPtrs']
1655        print >>pFile,'\n Modulation waves'
1656        names = {'Sfrac':['Fsin','Fcos','Fzero','Fwid'],'Spos':['Xsin','Ysin','Zsin','Xcos','Ycos','Zcos','Tmin','Tmax','Xmax','Ymax','Zmax'],
1657            'Sadp':['U11sin','U22sin','U33sin','U12sin','U13sin','U23sin','U11cos','U22cos',
1658            'U33cos','U12cos','U13cos','U23cos'],'Smag':['MXsin','MYsin','MZsin','MXcos','MYcos','MZcos']}
1659        print >>pFile,135*'-'
1660        for i,at in enumerate(Atoms):
1661            AtomSS = at[-1]['SS1']
1662            waveType = AtomSS['waveType']
1663            for Stype in ['Sfrac','Spos','Sadp','Smag']:
1664                Waves = AtomSS[Stype]
1665                if len(Waves):
1666                    print >>pFile,' atom: %s, site sym: %s, type: %s wave type: %s:'    \
1667                        %(at[ct-1],at[cs],Stype,waveType)
1668                    line = ''
1669                    for iw,wave in enumerate(Waves):
1670                        stiw = ':'+str(i)+':'+str(iw)
1671                        namstr = '  names :'
1672                        valstr = '  values:'
1673                        sigstr = '  esds  :'
1674                        if Stype == 'Spos':
1675                            nt = 6
1676                            ot = 0
1677                            if waveType in ['ZigZag','Block',] and not iw:
1678                                nt = 5
1679                                ot = 6
1680                            for j in range(nt):
1681                                name = names['Spos'][j+ot]
1682                                namstr += '%12s'%(name)
1683                                valstr += '%12.4f'%(wave[0][j])
1684                                if name+stiw in wavesSig:
1685                                    sigstr += '%12.4f'%(wavesSig[name+stiw])
1686                                else:
1687                                    sigstr += 12*' '
1688                        elif Stype == 'Sfrac':
1689                            ot = 0
1690                            if 'Crenel' in waveType and not iw:
1691                                ot = 2
1692                            for j in range(2):
1693                                name = names['Sfrac'][j+ot]
1694                                namstr += '%12s'%(names['Sfrac'][j+ot])
1695                                valstr += '%12.4f'%(wave[0][j])
1696                                if name+stiw in wavesSig:
1697                                    sigstr += '%12.4f'%(wavesSig[name+stiw])
1698                                else:
1699                                    sigstr += 12*' '
1700                        elif Stype == 'Sadp':
1701                            for j in range(12):
1702                                name = names['Sadp'][j]
1703                                namstr += '%10s'%(names['Sadp'][j])
1704                                valstr += '%10.6f'%(wave[0][j])
1705                                if name+stiw in wavesSig:
1706                                    sigstr += '%10.6f'%(wavesSig[name+stiw])
1707                                else:
1708                                    sigstr += 10*' '
1709                        elif Stype == 'Smag':
1710                            for j in range(6):
1711                                name = names['Smag'][j]
1712                                namstr += '%12s'%(names['Smag'][j])
1713                                valstr += '%12.4f'%(wave[0][j])
1714                                if name+stiw in wavesSig:
1715                                    sigstr += '%12.4f'%(wavesSig[name+stiw])
1716                                else:
1717                                    sigstr += 12*' '
1718                               
1719                    print >>pFile,namstr
1720                    print >>pFile,valstr
1721                    print >>pFile,sigstr
1722       
1723               
1724    def PrintRBObjPOAndSig(rbfx,rbsx):
1725        namstr = '  names :'
1726        valstr = '  values:'
1727        sigstr = '  esds  :'
1728        for i,px in enumerate(['Px:','Py:','Pz:']):
1729            name = pfx+rbfx+px+rbsx
1730            namstr += '%12s'%('Pos '+px[1])
1731            valstr += '%12.5f'%(parmDict[name])
1732            if name in sigDict:
1733                sigstr += '%12.5f'%(sigDict[name])
1734            else:
1735                sigstr += 12*' '
1736        for i,po in enumerate(['Oa:','Oi:','Oj:','Ok:']):
1737            name = pfx+rbfx+po+rbsx
1738            namstr += '%12s'%('Ori '+po[1])
1739            valstr += '%12.5f'%(parmDict[name])
1740            if name in sigDict:
1741                sigstr += '%12.5f'%(sigDict[name])
1742            else:
1743                sigstr += 12*' '
1744        print >>pFile,namstr
1745        print >>pFile,valstr
1746        print >>pFile,sigstr
1747       
1748    def PrintRBObjTLSAndSig(rbfx,rbsx,TLS):
1749        namstr = '  names :'
1750        valstr = '  values:'
1751        sigstr = '  esds  :'
1752        if 'N' not in TLS:
1753            print >>pFile,' Thermal motion:'
1754        if 'T' in TLS:
1755            for i,pt in enumerate(['T11:','T22:','T33:','T12:','T13:','T23:']):
1756                name = pfx+rbfx+pt+rbsx
1757                namstr += '%12s'%(pt[:3])
1758                valstr += '%12.5f'%(parmDict[name])
1759                if name in sigDict:
1760                    sigstr += '%12.5f'%(sigDict[name])
1761                else:
1762                    sigstr += 12*' '
1763            print >>pFile,namstr
1764            print >>pFile,valstr
1765            print >>pFile,sigstr
1766        if 'L' in TLS:
1767            namstr = '  names :'
1768            valstr = '  values:'
1769            sigstr = '  esds  :'
1770            for i,pt in enumerate(['L11:','L22:','L33:','L12:','L13:','L23:']):
1771                name = pfx+rbfx+pt+rbsx
1772                namstr += '%12s'%(pt[:3])
1773                valstr += '%12.3f'%(parmDict[name])
1774                if name in sigDict:
1775                    sigstr += '%12.3f'%(sigDict[name])
1776                else:
1777                    sigstr += 12*' '
1778            print >>pFile,namstr
1779            print >>pFile,valstr
1780            print >>pFile,sigstr
1781        if 'S' in TLS:
1782            namstr = '  names :'
1783            valstr = '  values:'
1784            sigstr = '  esds  :'
1785            for i,pt in enumerate(['S12:','S13:','S21:','S23:','S31:','S32:','SAA:','SBB:']):
1786                name = pfx+rbfx+pt+rbsx
1787                namstr += '%12s'%(pt[:3])
1788                valstr += '%12.4f'%(parmDict[name])
1789                if name in sigDict:
1790                    sigstr += '%12.4f'%(sigDict[name])
1791                else:
1792                    sigstr += 12*' '
1793            print >>pFile,namstr
1794            print >>pFile,valstr
1795            print >>pFile,sigstr
1796        if 'U' in TLS:
1797            name = pfx+rbfx+'U:'+rbsx
1798            namstr = '  names :'+'%12s'%('Uiso')
1799            valstr = '  values:'+'%12.5f'%(parmDict[name])
1800            if name in sigDict:
1801                sigstr = '  esds  :'+'%12.5f'%(sigDict[name])
1802            else:
1803                sigstr = '  esds  :'+12*' '
1804            print >>pFile,namstr
1805            print >>pFile,valstr
1806            print >>pFile,sigstr
1807       
1808    def PrintRBObjTorAndSig(rbsx):
1809        namstr = '  names :'
1810        valstr = '  values:'
1811        sigstr = '  esds  :'
1812        nTors = len(RBObj['Torsions'])
1813        if nTors:
1814            print >>pFile,' Torsions:'
1815            for it in range(nTors):
1816                name = pfx+'RBRTr;'+str(it)+':'+rbsx
1817                namstr += '%12s'%('Tor'+str(it))
1818                valstr += '%12.4f'%(parmDict[name])
1819                if name in sigDict:
1820                    sigstr += '%12.4f'%(sigDict[name])
1821            print >>pFile,namstr
1822            print >>pFile,valstr
1823            print >>pFile,sigstr
1824               
1825    def PrintSHtextureAndSig(textureData,SHtextureSig):
1826        print >>pFile,'\n Spherical harmonics texture: Order:' + str(textureData['Order'])
1827        names = ['omega','chi','phi']
1828        namstr = '  names :'
1829        ptstr =  '  values:'
1830        sigstr = '  esds  :'
1831        for name in names:
1832            namstr += '%12s'%(name)
1833            ptstr += '%12.3f'%(textureData['Sample '+name][1])
1834            if 'Sample '+name in SHtextureSig:
1835                sigstr += '%12.3f'%(SHtextureSig['Sample '+name])
1836            else:
1837                sigstr += 12*' '
1838        print >>pFile,namstr
1839        print >>pFile,ptstr
1840        print >>pFile,sigstr
1841        print >>pFile,'\n Texture coefficients:'
1842        SHcoeff = textureData['SH Coeff'][1]
1843        SHkeys = SHcoeff.keys()
1844        nCoeff = len(SHcoeff)
1845        nBlock = nCoeff/10+1
1846        iBeg = 0
1847        iFin = min(iBeg+10,nCoeff)
1848        for block in range(nBlock):
1849            namstr = '  names :'
1850            ptstr =  '  values:'
1851            sigstr = '  esds  :'
1852            for name in SHkeys[iBeg:iFin]:
1853                namstr += '%12s'%(name)
1854                ptstr += '%12.3f'%(SHcoeff[name])
1855                if name in SHtextureSig:
1856                    sigstr += '%12.3f'%(SHtextureSig[name])
1857                else:
1858                    sigstr += 12*' '
1859            print >>pFile,namstr
1860            print >>pFile,ptstr
1861            print >>pFile,sigstr
1862            iBeg += 10
1863            iFin = min(iBeg+10,nCoeff)
1864           
1865    print >>pFile,'\n Phases:'
1866    for phase in Phases:
1867        print >>pFile,' Result for phase: ',phase
1868        Phase = Phases[phase]
1869        General = Phase['General']
1870        SGData = General['SGData']
1871        Atoms = Phase['Atoms']
1872        if Atoms and not General.get('doPawley'):
1873            cx,ct,cs,cia = General['AtomPtrs']
1874            AtLookup = G2mth.FillAtomLookUp(Atoms,cia+8)
1875        cell = General['Cell']
1876        pId = Phase['pId']
1877        pfx = str(pId)+'::'
1878        if cell[0]:
1879            A,sigA = cellFill(pfx,SGData,parmDict,sigDict)
1880            cellSig = getCellEsd(pfx,SGData,A,covData)  #includes sigVol
1881            print >>pFile,' Reciprocal metric tensor: '
1882            ptfmt = "%15.9f"
1883            names = ['A11','A22','A33','A12','A13','A23']
1884            namstr = '  names :'
1885            ptstr =  '  values:'
1886            sigstr = '  esds  :'
1887            for name,a,siga in zip(names,A,sigA):
1888                namstr += '%15s'%(name)
1889                ptstr += ptfmt%(a)
1890                if siga:
1891                    sigstr += ptfmt%(siga)
1892                else:
1893                    sigstr += 15*' '
1894            print >>pFile,namstr
1895            print >>pFile,ptstr
1896            print >>pFile,sigstr
1897            cell[1:7] = G2lat.A2cell(A)
1898            cell[7] = G2lat.calc_V(A)
1899            print >>pFile,' New unit cell:'
1900            ptfmt = ["%12.6f","%12.6f","%12.6f","%12.4f","%12.4f","%12.4f","%12.3f"]
1901            names = ['a','b','c','alpha','beta','gamma','Volume']
1902            namstr = '  names :'
1903            ptstr =  '  values:'
1904            sigstr = '  esds  :'
1905            for name,fmt,a,siga in zip(names,ptfmt,cell[1:8],cellSig):
1906                namstr += '%12s'%(name)
1907                ptstr += fmt%(a)
1908                if siga:
1909                    sigstr += fmt%(siga)
1910                else:
1911                    sigstr += 12*' '
1912            print >>pFile,namstr
1913            print >>pFile,ptstr
1914            print >>pFile,sigstr
1915        ik = 6  #for Pawley stuff below
1916        if General.get('Modulated',False):
1917            ik = 7
1918            Vec,vRef,maxH = General['SuperVec']
1919            if vRef:
1920                print >>pFile,' New modulation vector:'
1921                namstr = '  names :'
1922                ptstr =  '  values:'
1923                sigstr = '  esds  :'
1924                for var in ['mV0','mV1','mV2']:
1925                    namstr += '%12s'%(pfx+var)
1926                    ptstr += '%12.6f'%(parmDict[pfx+var])
1927                    if pfx+var in sigDict:
1928                        sigstr += '%12.6f'%(sigDict[pfx+var])
1929                    else:
1930                        sigstr += 12*' '
1931                print >>pFile,namstr
1932                print >>pFile,ptstr
1933                print >>pFile,sigstr
1934           
1935        General['Mass'] = 0.
1936        if Phase['General'].get('doPawley'):
1937            pawleyRef = Phase['Pawley ref']
1938            for i,refl in enumerate(pawleyRef):
1939                key = pfx+'PWLref:'+str(i)
1940                refl[ik] = parmDict[key]
1941                if key in sigDict:
1942                    refl[ik+1] = sigDict[key]
1943                else:
1944                    refl[ik+1] = 0
1945        else:
1946            VRBIds = RBIds['Vector']
1947            RRBIds = RBIds['Residue']
1948            RBModels = Phase['RBModels']
1949            for irb,RBObj in enumerate(RBModels.get('Vector',[])):
1950                jrb = VRBIds.index(RBObj['RBId'])
1951                rbsx = str(irb)+':'+str(jrb)
1952                print >>pFile,' Vector rigid body parameters:'
1953                PrintRBObjPOAndSig('RBV',rbsx)
1954                PrintRBObjTLSAndSig('RBV',rbsx,RBObj['ThermalMotion'][0])
1955            for irb,RBObj in enumerate(RBModels.get('Residue',[])):
1956                jrb = RRBIds.index(RBObj['RBId'])
1957                rbsx = str(irb)+':'+str(jrb)
1958                print >>pFile,' Residue rigid body parameters:'
1959                PrintRBObjPOAndSig('RBR',rbsx)
1960                PrintRBObjTLSAndSig('RBR',rbsx,RBObj['ThermalMotion'][0])
1961                PrintRBObjTorAndSig(rbsx)
1962            atomsSig = {}
1963            wavesSig = {}
1964            cx,ct,cs,cia = General['AtomPtrs']
1965            for i,at in enumerate(Atoms):
1966                names = {cx:pfx+'Ax:'+str(i),cx+1:pfx+'Ay:'+str(i),cx+2:pfx+'Az:'+str(i),cx+3:pfx+'Afrac:'+str(i),
1967                    cia+1:pfx+'AUiso:'+str(i),cia+2:pfx+'AU11:'+str(i),cia+3:pfx+'AU22:'+str(i),cia+4:pfx+'AU33:'+str(i),
1968                    cia+5:pfx+'AU12:'+str(i),cia+6:pfx+'AU13:'+str(i),cia+7:pfx+'AU23:'+str(i)}
1969                for ind in range(cx,cx+4):
1970                    at[ind] = parmDict[names[ind]]
1971                    if ind in range(cx,cx+3):
1972                        name = names[ind].replace('A','dA')
1973                    else:
1974                        name = names[ind]
1975                    if name in sigDict:
1976                        atomsSig[str(i)+':'+str(ind)] = sigDict[name]
1977                if at[cia] == 'I':
1978                    at[cia+1] = parmDict[names[cia+1]]
1979                    if names[cia+1] in sigDict:
1980                        atomsSig['%d:%d'%(i,cia+1)] = sigDict[names[cia+1]]
1981                else:
1982                    for ind in range(cia+2,cia+8):
1983                        at[ind] = parmDict[names[ind]]
1984                        if names[ind] in sigDict:
1985                            atomsSig[str(i)+':'+str(ind)] = sigDict[names[ind]]
1986                ind = General['AtomTypes'].index(at[ct])
1987                General['Mass'] += General['AtomMass'][ind]*at[cx+3]*at[cx+5]
1988                if General.get('Modulated',False):
1989                    AtomSS = at[-1]['SS1']
1990                    waveType = AtomSS['waveType']
1991                    for Stype in ['Sfrac','Spos','Sadp','Smag']:
1992                        Waves = AtomSS[Stype]
1993                        for iw,wave in enumerate(Waves):
1994                            stiw = str(i)+':'+str(iw)
1995                            if Stype == 'Spos':
1996                                if waveType in ['ZigZag','Block',] and not iw:
1997                                    names = ['Tmin:'+stiw,'Tmax:'+stiw,'Xmax:'+stiw,'Ymax:'+stiw,'Zmax:'+stiw]
1998                                else:
1999                                    names = ['Xsin:'+stiw,'Ysin:'+stiw,'Zsin:'+stiw,
2000                                        'Xcos:'+stiw,'Ycos:'+stiw,'Zcos:'+stiw]
2001                            elif Stype == 'Sadp':
2002                                names = ['U11sin:'+stiw,'U22sin:'+stiw,'U33sin:'+stiw,
2003                                    'U12sin:'+stiw,'U13sin:'+stiw,'U23sin:'+stiw,
2004                                    'U11cos:'+stiw,'U22cos:'+stiw,'U33cos:'+stiw,
2005                                    'U12cos:'+stiw,'U13cos:'+stiw,'U23cos:'+stiw]
2006                            elif Stype == 'Sfrac':
2007                                if 'Crenel' in waveType and not iw:
2008                                    names = ['Fzero:'+stiw,'Fwid:'+stiw]
2009                                else:
2010                                    names = ['Fsin:'+stiw,'Fcos:'+stiw]
2011                            elif Stype == 'Smag':
2012                                names = ['MXsin:'+stiw,'MYsin:'+stiw,'MZsin:'+stiw,
2013                                    'MXcos:'+stiw,'MYcos:'+stiw,'MZcos:'+stiw]
2014                            for iname,name in enumerate(names):
2015                                AtomSS[Stype][iw][0][iname] = parmDict[pfx+name]
2016                                if pfx+name in sigDict:
2017                                    wavesSig[name] = sigDict[pfx+name]
2018                   
2019            PrintAtomsAndSig(General,Atoms,atomsSig)
2020            if General['Type'] == 'magnetic':
2021                PrintMomentsAndSig(General,Atoms,atomsSig)
2022            if General.get('Modulated',False):
2023                PrintWavesAndSig(General,Atoms,wavesSig)
2024           
2025       
2026        textureData = General['SH Texture']   
2027        if textureData['Order']:
2028            SHtextureSig = {}
2029            for name in ['omega','chi','phi']:
2030                aname = pfx+'SH '+name
2031                textureData['Sample '+name][1] = parmDict[aname]
2032                if aname in sigDict:
2033                    SHtextureSig['Sample '+name] = sigDict[aname]
2034            for name in textureData['SH Coeff'][1]:
2035                aname = pfx+name
2036                textureData['SH Coeff'][1][name] = parmDict[aname]
2037                if aname in sigDict:
2038                    SHtextureSig[name] = sigDict[aname]
2039            PrintSHtextureAndSig(textureData,SHtextureSig)
2040        if phase in RestraintDict and not Phase['General'].get('doPawley'):
2041            PrintRestraints(cell[1:7],SGData,General['AtomPtrs'],Atoms,AtLookup,
2042                textureData,RestraintDict[phase],pFile)
2043                   
2044################################################################################
2045##### Histogram & Phase data
2046################################################################################       
2047                   
2048def GetHistogramPhaseData(Phases,Histograms,Print=True,pFile=None,resetRefList=True):
2049    '''Loads the HAP histogram/phase information into dicts
2050
2051    :param dict Phases: phase information
2052    :param dict Histograms: Histogram information
2053    :param bool Print: prints information as it is read
2054    :param file pFile: file object to print to (the default, None causes printing to the console)
2055    :param bool resetRefList: Should the contents of the Reflection List be initialized
2056      on loading. The default, True, initializes the Reflection List as it is loaded.
2057
2058    :returns: (hapVary,hapDict,controlDict)
2059      * hapVary: list of refined variables
2060      * hapDict: dict with refined variables and their values
2061      * controlDict: dict with fixed parameters
2062    '''
2063   
2064    def PrintSize(hapData):
2065        if hapData[0] in ['isotropic','uniaxial']:
2066            line = '\n Size model    : %9s'%(hapData[0])
2067            line += ' equatorial:'+'%12.5f'%(hapData[1][0])+' Refine? '+str(hapData[2][0])
2068            if hapData[0] == 'uniaxial':
2069                line += ' axial:'+'%12.5f'%(hapData[1][1])+' Refine? '+str(hapData[2][1])
2070            line += '\n\t LG mixing coeff.: %12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
2071            print >>pFile,line
2072        else:
2073            print >>pFile,'\n Size model    : %s'%(hapData[0])+ \
2074                '\n\t LG mixing coeff.:%12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
2075            Snames = ['S11','S22','S33','S12','S13','S23']
2076            ptlbls = ' names :'
2077            ptstr =  ' values:'
2078            varstr = ' refine:'
2079            for i,name in enumerate(Snames):
2080                ptlbls += '%12s' % (name)
2081                ptstr += '%12.6f' % (hapData[4][i])
2082                varstr += '%12s' % (str(hapData[5][i]))
2083            print >>pFile,ptlbls
2084            print >>pFile,ptstr
2085            print >>pFile,varstr
2086       
2087    def PrintMuStrain(hapData,SGData):
2088        if hapData[0] in ['isotropic','uniaxial']:
2089            line = '\n Mustrain model: %9s'%(hapData[0])
2090            line += ' equatorial:'+'%12.1f'%(hapData[1][0])+' Refine? '+str(hapData[2][0])
2091            if hapData[0] == 'uniaxial':
2092                line += ' axial:'+'%12.1f'%(hapData[1][1])+' Refine? '+str(hapData[2][1])
2093            line +='\n\t LG mixing coeff.:%12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
2094            print >>pFile,line
2095        else:
2096            print >>pFile,'\n Mustrain model: %s'%(hapData[0])+ \
2097                '\n\t LG mixing coeff.:%12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
2098            Snames = G2spc.MustrainNames(SGData)
2099            ptlbls = ' names :'
2100            ptstr =  ' values:'
2101            varstr = ' refine:'
2102            for i,name in enumerate(Snames):
2103                ptlbls += '%12s' % (name)
2104                ptstr += '%12.6f' % (hapData[4][i])
2105                varstr += '%12s' % (str(hapData[5][i]))
2106            print >>pFile,ptlbls
2107            print >>pFile,ptstr
2108            print >>pFile,varstr
2109
2110    def PrintHStrain(hapData,SGData):
2111        print >>pFile,'\n Hydrostatic/elastic strain: '
2112        Hsnames = G2spc.HStrainNames(SGData)
2113        ptlbls = ' names :'
2114        ptstr =  ' values:'
2115        varstr = ' refine:'
2116        for i,name in enumerate(Hsnames):
2117            ptlbls += '%12s' % (name)
2118            ptstr += '%12.4g' % (hapData[0][i])
2119            varstr += '%12s' % (str(hapData[1][i]))
2120        print >>pFile,ptlbls
2121        print >>pFile,ptstr
2122        print >>pFile,varstr
2123
2124    def PrintSHPO(hapData):
2125        print >>pFile,'\n Spherical harmonics preferred orientation: Order:' + \
2126            str(hapData[4])+' Refine? '+str(hapData[2])
2127        ptlbls = ' names :'
2128        ptstr =  ' values:'
2129        for item in hapData[5]:
2130            ptlbls += '%12s'%(item)
2131            ptstr += '%12.3f'%(hapData[5][item]) 
2132        print >>pFile,ptlbls
2133        print >>pFile,ptstr
2134   
2135    def PrintBabinet(hapData):
2136        print >>pFile,'\n Babinet form factor modification: '
2137        ptlbls = ' names :'
2138        ptstr =  ' values:'
2139        varstr = ' refine:'
2140        for name in ['BabA','BabU']:
2141            ptlbls += '%12s' % (name)
2142            ptstr += '%12.6f' % (hapData[name][0])
2143            varstr += '%12s' % (str(hapData[name][1]))
2144        print >>pFile,ptlbls
2145        print >>pFile,ptstr
2146        print >>pFile,varstr
2147       
2148    hapDict = {}
2149    hapVary = []
2150    controlDict = {}
2151    poType = {}
2152    poAxes = {}
2153    spAxes = {}
2154    spType = {}
2155   
2156    for phase in Phases:
2157        HistoPhase = Phases[phase]['Histograms']
2158        SGData = Phases[phase]['General']['SGData']
2159        cell = Phases[phase]['General']['Cell'][1:7]
2160        A = G2lat.cell2A(cell)
2161        if Phases[phase]['General'].get('Modulated',False):
2162            SSGData = Phases[phase]['General']['SSGData']
2163            Vec,x,maxH = Phases[phase]['General']['SuperVec']
2164        pId = Phases[phase]['pId']
2165        histoList = HistoPhase.keys()
2166        histoList.sort()
2167        for histogram in histoList:
2168            try:
2169                Histogram = Histograms[histogram]
2170            except KeyError:                       
2171                #skip if histogram not included e.g. in a sequential refinement
2172                continue
2173            hapData = HistoPhase[histogram]
2174            hId = Histogram['hId']
2175            if 'PWDR' in histogram:
2176                limits = Histogram['Limits'][1]
2177                inst = Histogram['Instrument Parameters'][0]
2178                Zero = inst['Zero'][0]
2179                if 'C' in inst['Type'][1]:
2180                    try:
2181                        wave = inst['Lam'][1]
2182                    except KeyError:
2183                        wave = inst['Lam1'][1]
2184                    dmin = wave/(2.0*sind(limits[1]/2.0))
2185                elif 'T' in inst['Type'][0]:
2186                    dmin = limits[0]/inst['difC'][1]
2187                if Phases[phase]['General']['Type'] == 'magnetic':
2188                    dmin = max(dmin,Phases[phase]['General']['MagDmin'])
2189                pfx = str(pId)+':'+str(hId)+':'
2190                for item in ['Scale','Extinction']:
2191                    hapDict[pfx+item] = hapData[item][0]
2192                    if hapData[item][1]:
2193                        hapVary.append(pfx+item)
2194                names = G2spc.HStrainNames(SGData)
2195                HSvals = []
2196                for i,name in enumerate(names):
2197                    hapDict[pfx+name] = hapData['HStrain'][0][i]
2198                    HSvals.append(hapDict[pfx+name])
2199                    if hapData['HStrain'][1][i]:
2200                        hapVary.append(pfx+name)
2201                DIJS = G2spc.HStrainVals(HSvals,SGData)
2202                controlDict[pfx+'poType'] = hapData['Pref.Ori.'][0]
2203                if hapData['Pref.Ori.'][0] == 'MD':
2204                    hapDict[pfx+'MD'] = hapData['Pref.Ori.'][1]
2205                    controlDict[pfx+'MDAxis'] = hapData['Pref.Ori.'][3]
2206                    if hapData['Pref.Ori.'][2]:
2207                        hapVary.append(pfx+'MD')
2208                else:                           #'SH' spherical harmonics
2209                    controlDict[pfx+'SHord'] = hapData['Pref.Ori.'][4]
2210                    controlDict[pfx+'SHncof'] = len(hapData['Pref.Ori.'][5])
2211                    controlDict[pfx+'SHnames'] = G2lat.GenSHCoeff(SGData['SGLaue'],'0',controlDict[pfx+'SHord'],False)
2212                    controlDict[pfx+'SHhkl'] = []
2213                    try: #patch for old Pref.Ori. items
2214                        controlDict[pfx+'SHtoler'] = 0.1
2215                        if hapData['Pref.Ori.'][6][0] != '':
2216                            controlDict[pfx+'SHhkl'] = [eval(a.replace(' ',',')) for a in hapData['Pref.Ori.'][6]]
2217                        controlDict[pfx+'SHtoler'] = hapData['Pref.Ori.'][7]
2218                    except IndexError:
2219                        pass
2220                    for item in hapData['Pref.Ori.'][5]:
2221                        hapDict[pfx+item] = hapData['Pref.Ori.'][5][item]
2222                        if hapData['Pref.Ori.'][2]:
2223                            hapVary.append(pfx+item)
2224                for item in ['Mustrain','Size']:
2225                    controlDict[pfx+item+'Type'] = hapData[item][0]
2226                    hapDict[pfx+item+';mx'] = hapData[item][1][2]
2227                    if hapData[item][2][2]:
2228                        hapVary.append(pfx+item+';mx')
2229                    if hapData[item][0] in ['isotropic','uniaxial']:
2230                        hapDict[pfx+item+';i'] = hapData[item][1][0]
2231                        if hapData[item][2][0]:
2232                            hapVary.append(pfx+item+';i')
2233                        if hapData[item][0] == 'uniaxial':
2234                            controlDict[pfx+item+'Axis'] = hapData[item][3]
2235                            hapDict[pfx+item+';a'] = hapData[item][1][1]
2236                            if hapData[item][2][1]:
2237                                hapVary.append(pfx+item+';a')
2238                    else:       #generalized for mustrain or ellipsoidal for size
2239                        Nterms = len(hapData[item][4])
2240                        if item == 'Mustrain':
2241                            names = G2spc.MustrainNames(SGData)
2242                            pwrs = []
2243                            for name in names:
2244                                h,k,l = name[1:]
2245                                pwrs.append([int(h),int(k),int(l)])
2246                            controlDict[pfx+'MuPwrs'] = pwrs
2247                        for i in range(Nterms):
2248                            sfx = ';'+str(i)
2249                            hapDict[pfx+item+sfx] = hapData[item][4][i]
2250                            if hapData[item][5][i]:
2251                                hapVary.append(pfx+item+sfx)
2252                for bab in ['BabA','BabU']:
2253                    hapDict[pfx+bab] = hapData['Babinet'][bab][0]
2254                    if hapData['Babinet'][bab][1]:
2255                        hapVary.append(pfx+bab)
2256                               
2257                if Print: 
2258                    print >>pFile,'\n Phase: ',phase,' in histogram: ',histogram
2259                    print >>pFile,135*'-'
2260                    print >>pFile,' Phase fraction  : %10.4f'%(hapData['Scale'][0]),' Refine?',hapData['Scale'][1]
2261                    print >>pFile,' Extinction coeff: %10.4f'%(hapData['Extinction'][0]),' Refine?',hapData['Extinction'][1]
2262                    if hapData['Pref.Ori.'][0] == 'MD':
2263                        Ax = hapData['Pref.Ori.'][3]
2264                        print >>pFile,' March-Dollase PO: %10.4f'%(hapData['Pref.Ori.'][1]),' Refine?',hapData['Pref.Ori.'][2], \
2265                            ' Axis: %d %d %d'%(Ax[0],Ax[1],Ax[2])
2266                    else: #'SH' for spherical harmonics
2267                        PrintSHPO(hapData['Pref.Ori.'])
2268                        print >>pFile,' Penalty hkl list: '+str(controlDict[pfx+'SHhkl'])+' tolerance: %.2f'%(controlDict[pfx+'SHtoler'])
2269                    PrintSize(hapData['Size'])
2270                    PrintMuStrain(hapData['Mustrain'],SGData)
2271                    PrintHStrain(hapData['HStrain'],SGData)
2272                    if hapData['Babinet']['BabA'][0]:
2273                        PrintBabinet(hapData['Babinet'])                       
2274                if resetRefList:
2275                    refList = []
2276                    Uniq = []
2277                    Phi = []
2278                    if Phases[phase]['General'].get('Modulated',False):
2279                        ifSuper = True
2280                        HKLd = np.array(G2lat.GenSSHLaue(dmin,SGData,SSGData,Vec,maxH,A))
2281                        HKLd = G2mth.sortArray(HKLd,4,reverse=True)
2282                        for h,k,l,m,d in HKLd:
2283                            ext,mul,uniq,phi = G2spc.GenHKLf([h,k,l],SGData)
2284                            mul *= 2      # for powder overlap of Friedel pairs
2285                            if m or not ext:
2286                                if 'C' in inst['Type'][0]:
2287                                    pos = G2lat.Dsp2pos(inst,d)
2288                                    if limits[0] < pos < limits[1]:
2289                                        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])
2290                                        #... sig,gam,fotsq,fctsq, phase,icorr,prfo,abs,ext
2291                                        Uniq.append(uniq)
2292                                        Phi.append(phi)
2293                                elif 'T' in inst['Type'][0]:
2294                                    pos = G2lat.Dsp2pos(inst,d)
2295                                    if limits[0] < pos < limits[1]:
2296                                        wave = inst['difC'][1]*d/(252.816*inst['fltPath'][0])
2297                                        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])
2298                                        # ... sig,gam,fotsq,fctsq, phase,icorr,alp,bet,wave, prfo,abs,ext
2299                                        Uniq.append(uniq)
2300                                        Phi.append(phi)
2301                    else:
2302                        ifSuper = False
2303                        HKLd = np.array(G2lat.GenHLaue(dmin,SGData,A))
2304                        HKLd = G2mth.sortArray(HKLd,3,reverse=True)
2305                        for h,k,l,d in HKLd:
2306                            ext,mul,uniq,phi = G2spc.GenHKLf([h,k,l],SGData)
2307                            mul *= 2      # for powder overlap of Friedel pairs
2308                            if ext:
2309                                continue
2310                            if 'C' in inst['Type'][0]:
2311                                pos = G2lat.Dsp2pos(inst,d)
2312                                if limits[0] < pos < limits[1]:
2313                                    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])
2314                                    #... sig,gam,fotsq,fctsq, phase,icorr,prfo,abs,ext
2315                                    Uniq.append(uniq)
2316                                    Phi.append(phi)
2317                            elif 'T' in inst['Type'][0]:
2318                                pos = G2lat.Dsp2pos(inst,d)
2319                                if limits[0] < pos < limits[1]:
2320                                    wave = inst['difC'][1]*d/(252.816*inst['fltPath'][0])
2321                                    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])
2322                                    # ... sig,gam,fotsq,fctsq, phase,icorr,alp,bet,wave, prfo,abs,ext
2323                                    Uniq.append(uniq)
2324                                    Phi.append(phi)
2325                    Histogram['Reflection Lists'][phase] = {'RefList':np.array(refList),'FF':{},'Type':inst['Type'][0],'Super':ifSuper}
2326            elif 'HKLF' in histogram:
2327                inst = Histogram['Instrument Parameters'][0]
2328                hId = Histogram['hId']
2329                hfx = ':%d:'%(hId)
2330                for item in inst:
2331                    if type(inst) is not list and item != 'Type': continue # skip over non-refined items (such as InstName)
2332                    hapDict[hfx+item] = inst[item][1]
2333                pfx = str(pId)+':'+str(hId)+':'
2334                hapDict[pfx+'Scale'] = hapData['Scale'][0]
2335                if hapData['Scale'][1]:
2336                    hapVary.append(pfx+'Scale')
2337                               
2338                extApprox,extType,extParms = hapData['Extinction']
2339                controlDict[pfx+'EType'] = extType
2340                controlDict[pfx+'EApprox'] = extApprox
2341                if 'C' in inst['Type'][0]:
2342                    controlDict[pfx+'Tbar'] = extParms['Tbar']
2343                    controlDict[pfx+'Cos2TM'] = extParms['Cos2TM']
2344                if 'Primary' in extType:
2345                    Ekey = ['Ep',]
2346                elif 'I & II' in extType:
2347                    Ekey = ['Eg','Es']
2348                elif 'Secondary Type II' == extType:
2349                    Ekey = ['Es',]
2350                elif 'Secondary Type I' == extType:
2351                    Ekey = ['Eg',]
2352                else:   #'None'
2353                    Ekey = []
2354                for eKey in Ekey:
2355                    hapDict[pfx+eKey] = extParms[eKey][0]
2356                    if extParms[eKey][1]:
2357                        hapVary.append(pfx+eKey)
2358                for bab in ['BabA','BabU']:
2359                    hapDict[pfx+bab] = hapData['Babinet'][bab][0]
2360                    if hapData['Babinet'][bab][1]:
2361                        hapVary.append(pfx+bab)
2362                Twins = hapData.get('Twins',[[np.array([[1,0,0],[0,1,0],[0,0,1]]),[1.0,False,0]],])
2363                if len(Twins) == 1:
2364                    hapDict[pfx+'Flack'] = hapData.get('Flack',[0.,False])[0]
2365                    if hapData.get('Flack',[0,False])[1]:
2366                        hapVary.append(pfx+'Flack')
2367                sumTwFr = 0.
2368                controlDict[pfx+'TwinLaw'] = []
2369                controlDict[pfx+'TwinInv'] = []
2370                NTL = 0           
2371                for it,twin in enumerate(Twins):
2372                    if 'bool' in str(type(twin[0])):
2373                        controlDict[pfx+'TwinInv'].append(twin[0])
2374                        controlDict[pfx+'TwinLaw'].append(np.zeros((3,3)))
2375                    else:
2376                        NTL += 1
2377                        controlDict[pfx+'TwinInv'].append(False)
2378                        controlDict[pfx+'TwinLaw'].append(twin[0])
2379                    if it:
2380                        hapDict[pfx+'TwinFr:'+str(it)] = twin[1]
2381                        sumTwFr += twin[1]
2382                    else:
2383                        hapDict[pfx+'TwinFr:0'] = twin[1][0]
2384                        controlDict[pfx+'TwinNMN'] = twin[1][1]
2385                    if Twins[0][1][1]:
2386                        hapVary.append(pfx+'TwinFr:'+str(it))
2387                controlDict[pfx+'NTL'] = NTL
2388                #need to make constraint on TwinFr
2389                controlDict[pfx+'TwinLaw'] = np.array(controlDict[pfx+'TwinLaw'])
2390                if len(Twins) > 1:    #force sum to unity
2391                    hapDict[pfx+'TwinFr:0'] = 1.-sumTwFr
2392                if Print: 
2393                    print >>pFile,'\n Phase: ',phase,' in histogram: ',histogram
2394                    print >>pFile,135*'-'
2395                    print >>pFile,' Scale factor     : %10.4f'%(hapData['Scale'][0]),' Refine?',hapData['Scale'][1]
2396                    if extType != 'None':
2397                        print >>pFile,' Extinction  Type: %15s'%(extType),' approx: %10s'%(extApprox)
2398                        text = ' Parameters       :'
2399                        for eKey in Ekey:
2400                            text += ' %4s : %10.3e Refine? '%(eKey,extParms[eKey][0])+str(extParms[eKey][1])
2401                        print >>pFile,text
2402                    if hapData['Babinet']['BabA'][0]:
2403                        PrintBabinet(hapData['Babinet'])
2404                    if not SGData['SGInv'] and len(Twins) == 1:
2405                        print >>pFile,' Flack parameter: %10.3f'%(hapData['Flack'][0]),' Refine?',hapData['Flack'][1]
2406                    if len(Twins) > 1:
2407                        for it,twin in enumerate(Twins):
2408                            if 'bool' in str(type(twin[0])):
2409                                print >>pFile,' Nonmerohedral twin fr.: %5.3f Inv? %s Refine? '%(hapDict[pfx+'TwinFr:'+str(it)],str(controlDict[pfx+'TwinInv'][it])),Twins[0][1][1] 
2410                            else:
2411                                print >>pFile,' Twin law: %s'%(str(twin[0]).replace('\n',',')),' Twin fr.: %5.3f Refine? '%(hapDict[pfx+'TwinFr:'+str(it)]),Twins[0][1][1] 
2412                       
2413                Histogram['Reflection Lists'] = phase       
2414               
2415    return hapVary,hapDict,controlDict
2416   
2417def SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms,FFtables,Print=True,pFile=None):
2418    'needs a doc string'
2419   
2420    def PrintSizeAndSig(hapData,sizeSig):
2421        line = '\n Size model:     %9s'%(hapData[0])
2422        refine = False
2423        if hapData[0] in ['isotropic','uniaxial']:
2424            line += ' equatorial:%12.5f'%(hapData[1][0])
2425            if sizeSig[0][0]:
2426                line += ', sig:%8.4f'%(sizeSig[0][0])
2427                refine = True
2428            if hapData[0] == 'uniaxial':
2429                line += ' axial:%12.4f'%(hapData[1][1])
2430                if sizeSig[0][1]:
2431                    refine = True
2432                    line += ', sig:%8.4f'%(sizeSig[0][1])
2433            line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
2434            if sizeSig[0][2]:
2435                refine = True
2436                line += ', sig:%8.4f'%(sizeSig[0][2])
2437            if refine:
2438                print >>pFile,line
2439        else:
2440            line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
2441            if sizeSig[0][2]:
2442                refine = True
2443                line += ', sig:%8.4f'%(sizeSig[0][2])
2444            Snames = ['S11','S22','S33','S12','S13','S23']
2445            ptlbls = ' name  :'
2446            ptstr =  ' value :'
2447            sigstr = ' sig   :'
2448            for i,name in enumerate(Snames):
2449                ptlbls += '%12s' % (name)
2450                ptstr += '%12.6f' % (hapData[4][i])
2451                if sizeSig[1][i]:
2452                    refine = True
2453                    sigstr += '%12.6f' % (sizeSig[1][i])
2454                else:
2455                    sigstr += 12*' '
2456            if refine:
2457                print >>pFile,line
2458                print >>pFile,ptlbls
2459                print >>pFile,ptstr
2460                print >>pFile,sigstr
2461       
2462    def PrintMuStrainAndSig(hapData,mustrainSig,SGData):
2463        line = '\n Mustrain model: %9s'%(hapData[0])
2464        refine = False
2465        if hapData[0] in ['isotropic','uniaxial']:
2466            line += ' equatorial:%12.1f'%(hapData[1][0])
2467            if mustrainSig[0][0]:
2468                line += ', sig:%8.1f'%(mustrainSig[0][0])
2469                refine = True
2470            if hapData[0] == 'uniaxial':
2471                line += ' axial:%12.1f'%(hapData[1][1])
2472                if mustrainSig[0][1]:
2473                     line += ', sig:%8.1f'%(mustrainSig[0][1])
2474            line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
2475            if mustrainSig[0][2]:
2476                refine = True
2477                line += ', sig:%8.3f'%(mustrainSig[0][2])
2478            if refine:
2479                print >>pFile,line
2480        else:
2481            line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
2482            if mustrainSig[0][2]:
2483                refine = True
2484                line += ', sig:%8.3f'%(mustrainSig[0][2])
2485            Snames = G2spc.MustrainNames(SGData)
2486            ptlbls = ' name  :'
2487            ptstr =  ' value :'
2488            sigstr = ' sig   :'
2489            for i,name in enumerate(Snames):
2490                ptlbls += '%12s' % (name)
2491                ptstr += '%12.1f' % (hapData[4][i])
2492                if mustrainSig[1][i]:
2493                    refine = True
2494                    sigstr += '%12.1f' % (mustrainSig[1][i])
2495                else:
2496                    sigstr += 12*' '
2497            if refine:
2498                print >>pFile,line
2499                print >>pFile,ptlbls
2500                print >>pFile,ptstr
2501                print >>pFile,sigstr
2502           
2503    def PrintHStrainAndSig(hapData,strainSig,SGData):
2504        Hsnames = G2spc.HStrainNames(SGData)
2505        ptlbls = ' name  :'
2506        ptstr =  ' value :'
2507        sigstr = ' sig   :'
2508        refine = False
2509        for i,name in enumerate(Hsnames):
2510            ptlbls += '%12s' % (name)
2511            ptstr += '%12.4g' % (hapData[0][i])
2512            if name in strainSig:
2513                refine = True
2514                sigstr += '%12.4g' % (strainSig[name])
2515            else:
2516                sigstr += 12*' '
2517        if refine:
2518            print >>pFile,'\n Hydrostatic/elastic strain: '
2519            print >>pFile,ptlbls
2520            print >>pFile,ptstr
2521            print >>pFile,sigstr
2522       
2523    def PrintSHPOAndSig(pfx,hapData,POsig):
2524        print >>pFile,'\n Spherical harmonics preferred orientation: Order:'+str(hapData[4])
2525        ptlbls = ' names :'
2526        ptstr =  ' values:'
2527        sigstr = ' sig   :'
2528        for item in hapData[5]:
2529            ptlbls += '%12s'%(item)
2530            ptstr += '%12.3f'%(hapData[5][item])
2531            if pfx+item in POsig:
2532                sigstr += '%12.3f'%(POsig[pfx+item])
2533            else:
2534                sigstr += 12*' ' 
2535        print >>pFile,ptlbls
2536        print >>pFile,ptstr
2537        print >>pFile,sigstr
2538       
2539    def PrintExtAndSig(pfx,hapData,ScalExtSig):
2540        print >>pFile,'\n Single crystal extinction: Type: ',hapData[0],' Approx: ',hapData[1]
2541        text = ''
2542        for item in hapData[2]:
2543            if pfx+item in ScalExtSig:
2544                text += '       %s: '%(item)
2545                text += '%12.2e'%(hapData[2][item][0])
2546                if pfx+item in ScalExtSig:
2547                    text += ' sig: %12.2e'%(ScalExtSig[pfx+item])
2548        print >>pFile,text       
2549       
2550    def PrintBabinetAndSig(pfx,hapData,BabSig):
2551        print >>pFile,'\n Babinet form factor modification: '
2552        ptlbls = ' names :'
2553        ptstr =  ' values:'
2554        sigstr = ' sig   :'
2555        for item in hapData:
2556            ptlbls += '%12s'%(item)
2557            ptstr += '%12.3f'%(hapData[item][0])
2558            if pfx+item in BabSig:
2559                sigstr += '%12.3f'%(BabSig[pfx+item])
2560            else:
2561                sigstr += 12*' ' 
2562        print >>pFile,ptlbls
2563        print >>pFile,ptstr
2564        print >>pFile,sigstr
2565       
2566    def PrintTwinsAndSig(pfx,twinData,TwinSig):
2567        print >>pFile,'\n Twin Law fractions : '
2568        ptlbls = ' names :'
2569        ptstr =  ' values:'
2570        sigstr = ' sig   :'
2571        for it,item in enumerate(twinData):
2572            ptlbls += '     twin #%d'%(it)
2573            if it:
2574                ptstr += '%12.3f'%(item[1])
2575            else:
2576                ptstr += '%12.3f'%(item[1][0])
2577            if pfx+'TwinFr:'+str(it) in TwinSig:
2578                sigstr += '%12.3f'%(TwinSig[pfx+'TwinFr:'+str(it)])
2579            else:
2580                sigstr += 12*' ' 
2581        print >>pFile,ptlbls
2582        print >>pFile,ptstr
2583        print >>pFile,sigstr
2584       
2585   
2586    PhFrExtPOSig = {}
2587    SizeMuStrSig = {}
2588    ScalExtSig = {}
2589    BabSig = {}
2590    TwinFrSig = {}
2591    wtFrSum = {}
2592    for phase in Phases:
2593        HistoPhase = Phases[phase]['Histograms']
2594        General = Phases[phase]['General']
2595        SGData = General['SGData']
2596        pId = Phases[phase]['pId']
2597        histoList = HistoPhase.keys()
2598        histoList.sort()
2599        for histogram in histoList:
2600            try:
2601                Histogram = Histograms[histogram]
2602            except KeyError:                       
2603                #skip if histogram not included e.g. in a sequential refinement
2604                continue
2605            hapData = HistoPhase[histogram]
2606            hId = Histogram['hId']
2607            pfx = str(pId)+':'+str(hId)+':'
2608            if hId not in wtFrSum:
2609                wtFrSum[hId] = 0.
2610            if 'PWDR' in histogram:
2611                for item in ['Scale','Extinction']:
2612                    hapData[item][0] = parmDict[pfx+item]
2613                    if pfx+item in sigDict:
2614                        PhFrExtPOSig.update({pfx+item:sigDict[pfx+item],})
2615                wtFrSum[hId] += hapData['Scale'][0]*General['Mass']
2616                if hapData['Pref.Ori.'][0] == 'MD':
2617                    hapData['Pref.Ori.'][1] = parmDict[pfx+'MD']
2618                    if pfx+'MD' in sigDict:
2619                        PhFrExtPOSig.update({pfx+'MD':sigDict[pfx+'MD'],})
2620                else:                           #'SH' spherical harmonics
2621                    for item in hapData['Pref.Ori.'][5]:
2622                        hapData['Pref.Ori.'][5][item] = parmDict[pfx+item]
2623                        if pfx+item in sigDict:
2624                            PhFrExtPOSig.update({pfx+item:sigDict[pfx+item],})
2625                SizeMuStrSig.update({pfx+'Mustrain':[[0,0,0],[0 for i in range(len(hapData['Mustrain'][4]))]],
2626                    pfx+'Size':[[0,0,0],[0 for i in range(len(hapData['Size'][4]))]],
2627                    pfx+'HStrain':{}})                 
2628                for item in ['Mustrain','Size']:
2629                    hapData[item][1][2] = parmDict[pfx+item+';mx']
2630#                    hapData[item][1][2] = min(1.,max(0.,hapData[item][1][2]))
2631                    if pfx+item+';mx' in sigDict:
2632                        SizeMuStrSig[pfx+item][0][2] = sigDict[pfx+item+';mx']
2633                    if hapData[item][0] in ['isotropic','uniaxial']:                   
2634                        hapData[item][1][0] = parmDict[pfx+item+';i']
2635                        if item == 'Size':
2636                            hapData[item][1][0] = min(10.,max(0.001,hapData[item][1][0]))
2637                        if pfx+item+';i' in sigDict: 
2638                            SizeMuStrSig[pfx+item][0][0] = sigDict[pfx+item+';i']
2639                        if hapData[item][0] == 'uniaxial':
2640                            hapData[item][1][1] = parmDict[pfx+item+';a']
2641                            if item == 'Size':
2642                                hapData[item][1][1] = min(10.,max(0.001,hapData[item][1][1]))                       
2643                            if pfx+item+';a' in sigDict:
2644                                SizeMuStrSig[pfx+item][0][1] = sigDict[pfx+item+';a']
2645                    else:       #generalized for mustrain or ellipsoidal for size
2646                        Nterms = len(hapData[item][4])
2647                        for i in range(Nterms):
2648                            sfx = ';'+str(i)
2649                            hapData[item][4][i] = parmDict[pfx+item+sfx]
2650                            if pfx+item+sfx in sigDict:
2651                                SizeMuStrSig[pfx+item][1][i] = sigDict[pfx+item+sfx]
2652                names = G2spc.HStrainNames(SGData)
2653                for i,name in enumerate(names):
2654                    hapData['HStrain'][0][i] = parmDict[pfx+name]
2655                    if pfx+name in sigDict:
2656                        SizeMuStrSig[pfx+'HStrain'][name] = sigDict[pfx+name]
2657                for name in ['BabA','BabU']:
2658                    hapData['Babinet'][name][0] = parmDict[pfx+name]
2659                    if pfx+name in sigDict:
2660                        BabSig[pfx+name] = sigDict[pfx+name]               
2661               
2662            elif 'HKLF' in histogram:
2663                for item in ['Scale','Flack']:
2664                    if parmDict.get(pfx+item):
2665                        hapData[item][0] = parmDict[pfx+item]
2666                        if pfx+item in sigDict:
2667                            ScalExtSig[pfx+item] = sigDict[pfx+item]
2668                for item in ['Ep','Eg','Es']:
2669                    if parmDict.get(pfx+item):
2670                        hapData['Extinction'][2][item][0] = parmDict[pfx+item]
2671                        if pfx+item in sigDict:
2672                            ScalExtSig[pfx+item] = sigDict[pfx+item]
2673                for item in ['BabA','BabU']:
2674                    hapData['Babinet'][item][0] = parmDict[pfx+item]
2675                    if pfx+item in sigDict:
2676                        BabSig[pfx+item] = sigDict[pfx+item]
2677                if 'Twins' in hapData:
2678                    it = 1
2679                    sumTwFr = 0.
2680                    while True:
2681                        try:
2682                            hapData['Twins'][it][1] = parmDict[pfx+'TwinFr:'+str(it)]
2683                            if pfx+'TwinFr:'+str(it) in sigDict:
2684                                TwinFrSig[pfx+'TwinFr:'+str(it)] = sigDict[pfx+'TwinFr:'+str(it)]
2685                            if it:
2686                                sumTwFr += hapData['Twins'][it][1]
2687                            it += 1
2688                        except KeyError:
2689                            break
2690                    hapData['Twins'][0][1][0] = 1.-sumTwFr
2691
2692    if Print:
2693        for phase in Phases:
2694            HistoPhase = Phases[phase]['Histograms']
2695            General = Phases[phase]['General']
2696            SGData = General['SGData']
2697            pId = Phases[phase]['pId']
2698            histoList = HistoPhase.keys()
2699            histoList.sort()
2700            for histogram in histoList:
2701                try:
2702                    Histogram = Histograms[histogram]
2703                except KeyError:                       
2704                    #skip if histogram not included e.g. in a sequential refinement
2705                    continue
2706                print >>pFile,'\n Phase: ',phase,' in histogram: ',histogram
2707                print >>pFile,130*'-'
2708                hapData = HistoPhase[histogram]
2709                hId = Histogram['hId']
2710                Histogram['Residuals'][str(pId)+'::Name'] = phase
2711                pfx = str(pId)+':'+str(hId)+':'
2712                hfx = ':%s:'%(hId)
2713                if 'PWDR' in histogram:
2714                    print >>pFile,' Final refinement RF, RF^2 = %.2f%%, %.2f%% on %d reflections'   \
2715                        %(Histogram['Residuals'][pfx+'Rf'],Histogram['Residuals'][pfx+'Rf^2'],Histogram['Residuals'][pfx+'Nref'])
2716                    print >>pFile,' Bragg intensity sum = %.3g'%(Histogram['Residuals'][pfx+'sumInt'])
2717               
2718                    if pfx+'Scale' in PhFrExtPOSig:
2719                        wtFr = hapData['Scale'][0]*General['Mass']/wtFrSum[hId]
2720                        sigwtFr = PhFrExtPOSig[pfx+'Scale']*wtFr/hapData['Scale'][0]
2721                        print >>pFile,' Phase fraction  : %10.5f, sig %10.5f Weight fraction  : %8.5f, sig %10.5f' \
2722                            %(hapData['Scale'][0],PhFrExtPOSig[pfx+'Scale'],wtFr,sigwtFr)
2723                    if pfx+'Extinction' in PhFrExtPOSig:
2724                        print >>pFile,' Extinction coeff: %10.4f, sig %10.4f'%(hapData['Extinction'][0],PhFrExtPOSig[pfx+'Extinction'])
2725                    if hapData['Pref.Ori.'][0] == 'MD':
2726                        if pfx+'MD' in PhFrExtPOSig:
2727                            print >>pFile,' March-Dollase PO: %10.4f, sig %10.4f'%(hapData['Pref.Ori.'][1],PhFrExtPOSig[pfx+'MD'])
2728                    else:
2729                        PrintSHPOAndSig(pfx,hapData['Pref.Ori.'],PhFrExtPOSig)
2730                    PrintSizeAndSig(hapData['Size'],SizeMuStrSig[pfx+'Size'])
2731                    PrintMuStrainAndSig(hapData['Mustrain'],SizeMuStrSig[pfx+'Mustrain'],SGData)
2732                    PrintHStrainAndSig(hapData['HStrain'],SizeMuStrSig[pfx+'HStrain'],SGData)
2733                    if len(BabSig):
2734                        PrintBabinetAndSig(pfx,hapData['Babinet'],BabSig)
2735                   
2736                elif 'HKLF' in histogram:
2737                    Inst = Histogram['Instrument Parameters'][0]
2738                    print >>pFile,' Final refinement RF, RF^2 = %.2f%%, %.2f%% on %d reflections (%d user rejected, %d sp.gp.extinct)'   \
2739                        %(Histogram['Residuals'][pfx+'Rf'],Histogram['Residuals'][pfx+'Rf^2'],Histogram['Residuals'][pfx+'Nref'],
2740                        Histogram['Residuals'][pfx+'Nrej'],Histogram['Residuals'][pfx+'Next'])
2741                    if FFtables != None and 'N' not in Inst['Type'][0]:
2742                        PrintFprime(FFtables,hfx,pFile)
2743                    print >>pFile,' HKLF histogram weight factor = ','%.3f'%(Histogram['wtFactor'])
2744                    if pfx+'Scale' in ScalExtSig:
2745                        print >>pFile,' Scale factor : %10.4f, sig %10.4f'%(hapData['Scale'][0],ScalExtSig[pfx+'Scale'])
2746                    if hapData['Extinction'][0] != 'None':
2747                        PrintExtAndSig(pfx,hapData['Extinction'],ScalExtSig)
2748                    if len(BabSig):
2749                        PrintBabinetAndSig(pfx,hapData['Babinet'],BabSig)
2750                    if pfx+'Flack' in ScalExtSig:
2751                        print >>pFile,' Flack parameter : %10.3f, sig %10.3f'%(hapData['Flack'][0],ScalExtSig[pfx+'Flack'])
2752                    if len(TwinFrSig):
2753                        PrintTwinsAndSig(pfx,hapData['Twins'],TwinFrSig)
2754
2755################################################################################
2756##### Histogram data
2757################################################################################       
2758                   
2759def GetHistogramData(Histograms,Print=True,pFile=None):
2760    'needs a doc string'
2761   
2762    def GetBackgroundParms(hId,Background):
2763        Back = Background[0]
2764        DebyePeaks = Background[1]
2765        bakType,bakFlag = Back[:2]
2766        backVals = Back[3:]
2767        backNames = [':'+str(hId)+':Back;'+str(i) for i in range(len(backVals))]
2768        backDict = dict(zip(backNames,backVals))
2769        backVary = []
2770        if bakFlag:
2771            backVary = backNames
2772        backDict[':'+str(hId)+':nDebye'] = DebyePeaks['nDebye']
2773        backDict[':'+str(hId)+':nPeaks'] = DebyePeaks['nPeaks']
2774        debyeDict = {}
2775        debyeList = []
2776        for i in range(DebyePeaks['nDebye']):
2777            debyeNames = [':'+str(hId)+':DebyeA;'+str(i),':'+str(hId)+':DebyeR;'+str(i),':'+str(hId)+':DebyeU;'+str(i)]
2778            debyeDict.update(dict(zip(debyeNames,DebyePeaks['debyeTerms'][i][::2])))
2779            debyeList += zip(debyeNames,DebyePeaks['debyeTerms'][i][1::2])
2780        debyeVary = []
2781        for item in debyeList:
2782            if item[1]:
2783                debyeVary.append(item[0])
2784        backDict.update(debyeDict)
2785        backVary += debyeVary
2786        peakDict = {}
2787        peakList = []
2788        for i in range(DebyePeaks['nPeaks']):
2789            peakNames = [':'+str(hId)+':BkPkpos;'+str(i),':'+str(hId)+ \
2790                ':BkPkint;'+str(i),':'+str(hId)+':BkPksig;'+str(i),':'+str(hId)+':BkPkgam;'+str(i)]
2791            peakDict.update(dict(zip(peakNames,DebyePeaks['peaksList'][i][::2])))
2792            peakList += zip(peakNames,DebyePeaks['peaksList'][i][1::2])
2793        peakVary = []
2794        for item in peakList:
2795            if item[1]:
2796                peakVary.append(item[0])
2797        backDict.update(peakDict)
2798        backVary += peakVary
2799        return bakType,backDict,backVary           
2800       
2801    def GetInstParms(hId,Inst):     
2802        dataType = Inst['Type'][0]
2803        instDict = {}
2804        insVary = []
2805        pfx = ':'+str(hId)+':'
2806        insKeys = Inst.keys()
2807        insKeys.sort()
2808        for item in insKeys:
2809            insName = pfx+item
2810            instDict[insName] = Inst[item][1]
2811            if len(Inst[item]) > 2 and Inst[item][2]:
2812                insVary.append(insName)
2813        if 'C' in dataType:
2814            instDict[pfx+'SH/L'] = max(instDict[pfx+'SH/L'],0.0005)
2815        return dataType,instDict,insVary
2816       
2817    def GetSampleParms(hId,Sample):
2818        sampVary = []
2819        hfx = ':'+str(hId)+':'       
2820        sampDict = {hfx+'Gonio. radius':Sample['Gonio. radius'],hfx+'Omega':Sample['Omega'],
2821            hfx+'Chi':Sample['Chi'],hfx+'Phi':Sample['Phi']}
2822        for key in ('Temperature','Pressure','FreePrm1','FreePrm2','FreePrm3'):
2823            if key in Sample:
2824                sampDict[hfx+key] = Sample[key]
2825        Type = Sample['Type']
2826        if 'Bragg' in Type:             #Bragg-Brentano
2827            for item in ['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']:
2828                sampDict[hfx+item] = Sample[item][0]
2829                if Sample[item][1]:
2830                    sampVary.append(hfx+item)
2831        elif 'Debye' in Type:        #Debye-Scherrer
2832            for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
2833                sampDict[hfx+item] = Sample[item][0]
2834                if Sample[item][1]:
2835                    sampVary.append(hfx+item)
2836        return Type,sampDict,sampVary
2837       
2838    def PrintBackground(Background):
2839        Back = Background[0]
2840        DebyePeaks = Background[1]
2841        print >>pFile,'\n Background function: ',Back[0],' Refine?',bool(Back[1])
2842        line = ' Coefficients: '
2843        for i,back in enumerate(Back[3:]):
2844            line += '%10.3f'%(back)
2845            if i and not i%10:
2846                line += '\n'+15*' '
2847        print >>pFile,line
2848        if DebyePeaks['nDebye']:
2849            print >>pFile,'\n Debye diffuse scattering coefficients'
2850            parms = ['DebyeA','DebyeR','DebyeU']
2851            line = ' names :  '
2852            for parm in parms:
2853                line += '%8s refine?'%(parm)
2854            print >>pFile,line
2855            for j,term in enumerate(DebyePeaks['debyeTerms']):
2856                line = ' term'+'%2d'%(j)+':'
2857                for i in range(3):
2858                    line += '%10.3f %5s'%(term[2*i],bool(term[2*i+1]))                   
2859                print >>pFile,line
2860        if DebyePeaks['nPeaks']:
2861            print >>pFile,'\n Single peak coefficients'
2862            parms =    ['BkPkpos','BkPkint','BkPksig','BkPkgam']
2863            line = ' names :  '
2864            for parm in parms:
2865                line += '%8s refine?'%(parm)
2866            print >>pFile,line
2867            for j,term in enumerate(DebyePeaks['peaksList']):
2868                line = ' peak'+'%2d'%(j)+':'
2869                for i in range(4):
2870                    line += '%12.3f %5s'%(term[2*i],bool(term[2*i+1]))                   
2871                print >>pFile,line
2872       
2873    def PrintInstParms(Inst):
2874        print >>pFile,'\n Instrument Parameters:'
2875        insKeys = Inst.keys()
2876        insKeys.sort()
2877        iBeg = 0
2878        Ok = True
2879        while Ok:
2880            ptlbls = ' name  :'
2881            ptstr =  ' value :'
2882            varstr = ' refine:'
2883            iFin = min(iBeg+9,len(insKeys))
2884            for item in insKeys[iBeg:iFin]:
2885                if item not in ['Type','Source']:
2886                    ptlbls += '%12s' % (item)
2887                    ptstr += '%12.6f' % (Inst[item][1])
2888                    if item in ['Lam1','Lam2','Azimuth','fltPath','2-theta',]:
2889                        varstr += 12*' '
2890                    else:
2891                        varstr += '%12s' % (str(bool(Inst[item][2])))
2892            print >>pFile,ptlbls
2893            print >>pFile,ptstr
2894            print >>pFile,varstr
2895            iBeg = iFin
2896            if iBeg == len(insKeys):
2897                Ok = False
2898            else:
2899                print >>pFile,'\n'
2900       
2901    def PrintSampleParms(Sample):
2902        print >>pFile,'\n Sample Parameters:'
2903        print >>pFile,' Goniometer omega = %.2f, chi = %.2f, phi = %.2f'% \
2904            (Sample['Omega'],Sample['Chi'],Sample['Phi'])
2905        ptlbls = ' name  :'
2906        ptstr =  ' value :'
2907        varstr = ' refine:'
2908        if 'Bragg' in Sample['Type']:
2909            for item in ['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']:
2910                ptlbls += '%14s'%(item)
2911                ptstr += '%14.4f'%(Sample[item][0])
2912                varstr += '%14s'%(str(bool(Sample[item][1])))
2913           
2914        elif 'Debye' in Type:        #Debye-Scherrer
2915            for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
2916                ptlbls += '%14s'%(item)
2917                ptstr += '%14.4f'%(Sample[item][0])
2918                varstr += '%14s'%(str(bool(Sample[item][1])))
2919
2920        print >>pFile,ptlbls
2921        print >>pFile,ptstr
2922        print >>pFile,varstr
2923       
2924    histDict = {}
2925    histVary = []
2926    controlDict = {}
2927    histoList = Histograms.keys()
2928    histoList.sort()
2929    for histogram in histoList:
2930        if 'PWDR' in histogram:
2931            Histogram = Histograms[histogram]
2932            hId = Histogram['hId']
2933            pfx = ':'+str(hId)+':'
2934            controlDict[pfx+'wtFactor'] = Histogram['wtFactor']
2935            controlDict[pfx+'Limits'] = Histogram['Limits'][1]
2936            controlDict[pfx+'Exclude'] = Histogram['Limits'][2:]
2937            for excl in controlDict[pfx+'Exclude']:
2938                Histogram['Data'][0] = ma.masked_inside(Histogram['Data'][0],excl[0],excl[1])
2939            if controlDict[pfx+'Exclude']:
2940                ma.mask_rows(Histogram['Data'])
2941            Background = Histogram['Background']
2942            Type,bakDict,bakVary = GetBackgroundParms(hId,Background)
2943            controlDict[pfx+'bakType'] = Type
2944            histDict.update(bakDict)
2945            histVary += bakVary
2946           
2947            Inst = Histogram['Instrument Parameters'][0]
2948            Type,instDict,insVary = GetInstParms(hId,Inst)
2949            controlDict[pfx+'histType'] = Type
2950            if 'XC' in Type:
2951                if pfx+'Lam1' in instDict:
2952                    controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam1']
2953                else:
2954                    controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam']           
2955            histDict.update(instDict)
2956            histVary += insVary
2957           
2958            Sample = Histogram['Sample Parameters']
2959            Type,sampDict,sampVary = GetSampleParms(hId,Sample)
2960            controlDict[pfx+'instType'] = Type
2961            histDict.update(sampDict)
2962            histVary += sampVary
2963           
2964   
2965            if Print: 
2966                print >>pFile,'\n Histogram: ',histogram,' histogram Id: ',hId
2967                print >>pFile,135*'-'
2968                Units = {'C':' deg','T':' msec'}
2969                units = Units[controlDict[pfx+'histType'][2]]
2970                Limits = controlDict[pfx+'Limits']
2971                print >>pFile,' Instrument type: ',Sample['Type']
2972                print >>pFile,' Histogram limits: %8.2f%s to %8.2f%s'%(Limits[0],units,Limits[1],units)
2973                if len(controlDict[pfx+'Exclude']):
2974                    excls = controlDict[pfx+'Exclude']
2975                    for excl in excls:
2976                        print >>pFile,' Excluded region:  %8.2f%s to %8.2f%s'%(excl[0],units,excl[1],units)   
2977                PrintSampleParms(Sample)
2978                PrintInstParms(Inst)
2979                PrintBackground(Background)
2980        elif 'HKLF' in histogram:
2981            Histogram = Histograms[histogram]
2982            hId = Histogram['hId']
2983            pfx = ':'+str(hId)+':'
2984            controlDict[pfx+'wtFactor'] = Histogram['wtFactor']
2985            Inst = Histogram['Instrument Parameters'][0]
2986            controlDict[pfx+'histType'] = Inst['Type'][0]
2987            if 'X' in Inst['Type'][0]:
2988                histDict[pfx+'Lam'] = Inst['Lam'][1]
2989                controlDict[pfx+'keV'] = 12.397639/histDict[pfx+'Lam']                   
2990    return histVary,histDict,controlDict
2991   
2992def SetHistogramData(parmDict,sigDict,Histograms,FFtables,Print=True,pFile=None):
2993    'needs a doc string'
2994   
2995    def SetBackgroundParms(pfx,Background,parmDict,sigDict):
2996        Back = Background[0]
2997        DebyePeaks = Background[1]
2998        lenBack = len(Back[3:])
2999        backSig = [0 for i in range(lenBack+3*DebyePeaks['nDebye']+4*DebyePeaks['nPeaks'])]
3000        for i in range(lenBack):
3001            Back[3+i] = parmDict[pfx+'Back;'+str(i)]
3002            if pfx+'Back;'+str(i) in sigDict:
3003                backSig[i] = sigDict[pfx+'Back;'+str(i)]
3004        if DebyePeaks['nDebye']:
3005            for i in range(DebyePeaks['nDebye']):
3006                names = [pfx+'DebyeA;'+str(i),pfx+'DebyeR;'+str(i),pfx+'DebyeU;'+str(i)]
3007                for j,name in enumerate(names):
3008                    DebyePeaks['debyeTerms'][i][2*j] = parmDict[name]
3009                    if name in sigDict:
3010                        backSig[lenBack+3*i+j] = sigDict[name]           
3011        if DebyePeaks['nPeaks']:
3012            for i in range(DebyePeaks['nPeaks']):
3013                names = [pfx+'BkPkpos;'+str(i),pfx+'BkPkint;'+str(i),
3014                    pfx+'BkPksig;'+str(i),pfx+'BkPkgam;'+str(i)]
3015                for j,name in enumerate(names):
3016                    DebyePeaks['peaksList'][i][2*j] = parmDict[name]
3017                    if name in sigDict:
3018                        backSig[lenBack+3*DebyePeaks['nDebye']+4*i+j] = sigDict[name]
3019        return backSig
3020       
3021    def SetInstParms(pfx,Inst,parmDict,sigDict):
3022        instSig = {}
3023        insKeys = Inst.keys()
3024        insKeys.sort()
3025        for item in insKeys:
3026            insName = pfx+item
3027            Inst[item][1] = parmDict[insName]
3028            if insName in sigDict:
3029                instSig[item] = sigDict[insName]
3030            else:
3031                instSig[item] = 0
3032        return instSig
3033       
3034    def SetSampleParms(pfx,Sample,parmDict,sigDict):
3035        if 'Bragg' in Sample['Type']:             #Bragg-Brentano
3036            sampSig = [0 for i in range(5)]
3037            for i,item in enumerate(['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']):
3038                Sample[item][0] = parmDict[pfx+item]
3039                if pfx+item in sigDict:
3040                    sampSig[i] = sigDict[pfx+item]
3041        elif 'Debye' in Sample['Type']:        #Debye-Scherrer
3042            sampSig = [0 for i in range(4)]
3043            for i,item in enumerate(['Scale','Absorption','DisplaceX','DisplaceY']):
3044                Sample[item][0] = parmDict[pfx+item]
3045                if pfx+item in sigDict:
3046                    sampSig[i] = sigDict[pfx+item]
3047        return sampSig
3048       
3049    def PrintBackgroundSig(Background,backSig):
3050        Back = Background[0]
3051        DebyePeaks = Background[1]
3052        lenBack = len(Back[3:])
3053        valstr = ' value : '
3054        sigstr = ' sig   : '
3055        refine = False
3056        for i,back in enumerate(Back[3:]):
3057            valstr += '%10.4g'%(back)
3058            if Back[1]:
3059                refine = True
3060                sigstr += '%10.4g'%(backSig[i])
3061            else:
3062                sigstr += 10*' '
3063        if refine:
3064            print >>pFile,'\n Background function: ',Back[0]
3065            print >>pFile,valstr
3066            print >>pFile,sigstr
3067        if DebyePeaks['nDebye']:
3068            ifAny = False
3069            ptfmt = "%12.3f"
3070            names =  ' names :'
3071            ptstr =  ' values:'
3072            sigstr = ' esds  :'
3073            for item in sigDict:
3074                if 'Debye' in item:
3075                    ifAny = True
3076                    names += '%12s'%(item)
3077                    ptstr += ptfmt%(parmDict[item])
3078                    sigstr += ptfmt%(sigDict[item])
3079            if ifAny:
3080                print >>pFile,'\n Debye diffuse scattering coefficients'
3081                print >>pFile,names
3082                print >>pFile,ptstr
3083                print >>pFile,sigstr
3084        if DebyePeaks['nPeaks']:
3085            print >>pFile,'\n Single peak coefficients:'
3086            parms =    ['BkPkpos','BkPkint','BkPksig','BkPkgam']
3087            line = ' peak no. '
3088            for parm in parms:
3089                line += '%14s%12s'%(parm.center(14),'esd'.center(12))
3090            print >>pFile,line
3091            for ip in range(DebyePeaks['nPeaks']):
3092                ptstr = ' %4d '%(ip)
3093                for parm in parms:
3094                    fmt = '%14.3f'
3095                    efmt = '%12.3f'
3096                    if parm == 'BkPkpos':
3097                        fmt = '%14.4f'
3098                        efmt = '%12.4f'
3099                    name = pfx+parm+';%d'%(ip)
3100                    ptstr += fmt%(parmDict[name])
3101                    if name in sigDict:
3102                        ptstr += efmt%(sigDict[name])
3103                    else:
3104                        ptstr += 12*' '
3105                print >>pFile,ptstr
3106        sumBk = np.array(Histogram['sumBk'])
3107        print >>pFile,' Background sums: empirical %.3g, Debye %.3g, peaks %.3g, Total %.3g'    \
3108            %(sumBk[0],sumBk[1],sumBk[2],np.sum(sumBk))
3109       
3110    def PrintInstParmsSig(Inst,instSig):
3111        refine = False
3112        insKeys = instSig.keys()
3113        insKeys.sort()
3114        iBeg = 0
3115        Ok = True
3116        while Ok:
3117            ptlbls = ' names :'
3118            ptstr =  ' value :'
3119            sigstr = ' sig   :'
3120            iFin = min(iBeg+9,len(insKeys))
3121            for name in insKeys[iBeg:iFin]:
3122                if name not in  ['Type','Lam1','Lam2','Azimuth','Source','fltPath']:
3123                    ptlbls += '%12s' % (name)
3124                    ptstr += '%12.6f' % (Inst[name][1])
3125                    if instSig[name]:
3126                        refine = True
3127                        sigstr += '%12.6f' % (instSig[name])
3128                    else:
3129                        sigstr += 12*' '
3130            if refine:
3131                print >>pFile,'\n Instrument Parameters:'
3132                print >>pFile,ptlbls
3133                print >>pFile,ptstr
3134                print >>pFile,sigstr
3135            iBeg = iFin
3136            if iBeg == len(insKeys):
3137                Ok = False
3138       
3139    def PrintSampleParmsSig(Sample,sampleSig):
3140        ptlbls = ' names :'
3141        ptstr =  ' values:'
3142        sigstr = ' sig   :'
3143        refine = False
3144        if 'Bragg' in Sample['Type']:
3145            for i,item in enumerate(['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']):
3146                ptlbls += '%14s'%(item)
3147                ptstr += '%14.4f'%(Sample[item][0])
3148                if sampleSig[i]:
3149                    refine = True
3150                    sigstr += '%14.4f'%(sampleSig[i])
3151                else:
3152                    sigstr += 14*' '
3153           
3154        elif 'Debye' in Sample['Type']:        #Debye-Scherrer
3155            for i,item in enumerate(['Scale','Absorption','DisplaceX','DisplaceY']):
3156                ptlbls += '%14s'%(item)
3157                ptstr += '%14.4f'%(Sample[item][0])
3158                if sampleSig[i]:
3159                    refine = True
3160                    sigstr += '%14.4f'%(sampleSig[i])
3161                else:
3162                    sigstr += 14*' '
3163
3164        if refine:
3165            print >>pFile,'\n Sample Parameters:'
3166            print >>pFile,ptlbls
3167            print >>pFile,ptstr
3168            print >>pFile,sigstr
3169       
3170    histoList = Histograms.keys()
3171    histoList.sort()
3172    for histogram in histoList:
3173        if 'PWDR' in histogram:
3174            Histogram = Histograms[histogram]
3175            hId = Histogram['hId']
3176            pfx = ':'+str(hId)+':'
3177            Background = Histogram['Background']
3178            backSig = SetBackgroundParms(pfx,Background,parmDict,sigDict)
3179           
3180            Inst = Histogram['Instrument Parameters'][0]
3181            instSig = SetInstParms(pfx,Inst,parmDict,sigDict)
3182       
3183            Sample = Histogram['Sample Parameters']
3184            sampSig = SetSampleParms(pfx,Sample,parmDict,sigDict)
3185
3186            print >>pFile,'\n Histogram: ',histogram,' histogram Id: ',hId
3187            print >>pFile,135*'-'
3188            print >>pFile,' PWDR histogram weight factor = '+'%.3f'%(Histogram['wtFactor'])
3189            print >>pFile,' Final refinement wR = %.2f%% on %d observations in this histogram'%(Histogram['Residuals']['wR'],Histogram['Residuals']['Nobs'])
3190            print >>pFile,' Other residuals: R = %.2f%%, Rb = %.2f%%, wRb = %.2f%% wRmin = %.2f%%'% \
3191                (Histogram['Residuals']['R'],Histogram['Residuals']['Rb'],Histogram['Residuals']['wRb'],Histogram['Residuals']['wRmin'])
3192            if Print:
3193                print >>pFile,' Instrument type: ',Sample['Type']
3194                if FFtables != None and 'N' not in Inst['Type'][0]:
3195                    PrintFprime(FFtables,pfx,pFile)
3196                PrintSampleParmsSig(Sample,sampSig)
3197                PrintInstParmsSig(Inst,instSig)
3198                PrintBackgroundSig(Background,backSig)
3199               
Note: See TracBrowser for help on using the repository browser.