source: trunk/GSASIIstrIO.py @ 2495

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

some cleanup of Data tab
fix (!) crash after Refine for Costraints & Rigid Bodies - DELETED windows problem
remove Babinet & Flack from magnetic structures
remove magnetic from StructureFactorDeriv2
make new StructureFactorDerivMag? & work in that

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