source: trunk/GSASIIstrIO.py @ 2546

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

cleanup of all GSASII main routines - remove unused variables, dead code, etc.

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