source: trunk/GSASIIstrIO.py @ 981

Last change on this file since 981 was 981, checked in by toby, 8 years ago

introduce regress option; fix esd printing; more docs; new Mac app with drag & drop for open; control reset of ref list on load

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