source: trunk/GSASIIstrIO.py @ 1956

Last change on this file since 1956 was 1956, checked in by vondreele, 6 years ago

fix (again) fill unit cell, dynamic SS drawing & GenAtom? issues - now all OK
fixes to get SS structure factors to compute & twins in SS data issues

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