source: trunk/GSASIIstrIO.py @ 2478

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

work up to magnetic structure factor calcs.

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