source: trunk/GSASIIstrIO.py @ 1615

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

more work on incommensurate wave input to LS
fix FillAtomLookup? problem in G2restrGUI

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