source: trunk/GSASIIstrIO.py @ 1618

Last change on this file since 1618 was 1618, checked in by vondreele, 7 years ago

fix space group print in shelx export
fix errors when data is deleted
add HKLshow to change in space group in index page
more on supersymmetry refinement

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