source: trunk/GSASIIstrIO.py @ 2520

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

remove unused & commented Srr fctr & deriv routines
Mag moment derivatives semireasonable but not close
transform nucl - mag ok if no off diag terms in transform matrix

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