source: trunk/GSASIIstrIO.py @ 960

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

check in CIF dev snapshot

File size: 108.5 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):
1517    'needs a doc string'
1518   
1519    def PrintSize(hapData):
1520        if hapData[0] in ['isotropic','uniaxial']:
1521            line = '\n Size model    : %9s'%(hapData[0])
1522            line += ' equatorial:'+'%12.3f'%(hapData[1][0])+' Refine? '+str(hapData[2][0])
1523            if hapData[0] == 'uniaxial':
1524                line += ' axial:'+'%12.3f'%(hapData[1][1])+' Refine? '+str(hapData[2][1])
1525            line += '\n\t LG mixing coeff.: %12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
1526            print >>pFile,line
1527        else:
1528            print >>pFile,'\n Size model    : %s'%(hapData[0])+ \
1529                '\n\t LG mixing coeff.:%12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
1530            Snames = ['S11','S22','S33','S12','S13','S23']
1531            ptlbls = ' names :'
1532            ptstr =  ' values:'
1533            varstr = ' refine:'
1534            for i,name in enumerate(Snames):
1535                ptlbls += '%12s' % (name)
1536                ptstr += '%12.6f' % (hapData[4][i])
1537                varstr += '%12s' % (str(hapData[5][i]))
1538            print >>pFile,ptlbls
1539            print >>pFile,ptstr
1540            print >>pFile,varstr
1541       
1542    def PrintMuStrain(hapData,SGData):
1543        if hapData[0] in ['isotropic','uniaxial']:
1544            line = '\n Mustrain model: %9s'%(hapData[0])
1545            line += ' equatorial:'+'%12.1f'%(hapData[1][0])+' Refine? '+str(hapData[2][0])
1546            if hapData[0] == 'uniaxial':
1547                line += ' axial:'+'%12.1f'%(hapData[1][1])+' Refine? '+str(hapData[2][1])
1548            line +='\n\t LG mixing coeff.:%12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
1549            print >>pFile,line
1550        else:
1551            print >>pFile,'\n Mustrain model: %s'%(hapData[0])+ \
1552                '\n\t LG mixing coeff.:%12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
1553            Snames = G2spc.MustrainNames(SGData)
1554            ptlbls = ' names :'
1555            ptstr =  ' values:'
1556            varstr = ' refine:'
1557            for i,name in enumerate(Snames):
1558                ptlbls += '%12s' % (name)
1559                ptstr += '%12.6f' % (hapData[4][i])
1560                varstr += '%12s' % (str(hapData[5][i]))
1561            print >>pFile,ptlbls
1562            print >>pFile,ptstr
1563            print >>pFile,varstr
1564
1565    def PrintHStrain(hapData,SGData):
1566        print >>pFile,'\n Hydrostatic/elastic strain: '
1567        Hsnames = G2spc.HStrainNames(SGData)
1568        ptlbls = ' names :'
1569        ptstr =  ' values:'
1570        varstr = ' refine:'
1571        for i,name in enumerate(Hsnames):
1572            ptlbls += '%12s' % (name)
1573            ptstr += '%12.6f' % (hapData[0][i])
1574            varstr += '%12s' % (str(hapData[1][i]))
1575        print >>pFile,ptlbls
1576        print >>pFile,ptstr
1577        print >>pFile,varstr
1578
1579    def PrintSHPO(hapData):
1580        print >>pFile,'\n Spherical harmonics preferred orientation: Order:' + \
1581            str(hapData[4])+' Refine? '+str(hapData[2])
1582        ptlbls = ' names :'
1583        ptstr =  ' values:'
1584        for item in hapData[5]:
1585            ptlbls += '%12s'%(item)
1586            ptstr += '%12.3f'%(hapData[5][item]) 
1587        print >>pFile,ptlbls
1588        print >>pFile,ptstr
1589   
1590    def PrintBabinet(hapData):
1591        print >>pFile,'\n Babinet form factor modification: '
1592        ptlbls = ' names :'
1593        ptstr =  ' values:'
1594        varstr = ' refine:'
1595        for name in ['BabA','BabU']:
1596            ptlbls += '%12s' % (name)
1597            ptstr += '%12.6f' % (hapData[name][0])
1598            varstr += '%12s' % (str(hapData[name][1]))
1599        print >>pFile,ptlbls
1600        print >>pFile,ptstr
1601        print >>pFile,varstr
1602       
1603   
1604    hapDict = {}
1605    hapVary = []
1606    controlDict = {}
1607    poType = {}
1608    poAxes = {}
1609    spAxes = {}
1610    spType = {}
1611   
1612    for phase in Phases:
1613        HistoPhase = Phases[phase]['Histograms']
1614        SGData = Phases[phase]['General']['SGData']
1615        cell = Phases[phase]['General']['Cell'][1:7]
1616        A = G2lat.cell2A(cell)
1617        pId = Phases[phase]['pId']
1618        histoList = HistoPhase.keys()
1619        histoList.sort()
1620        for histogram in histoList:
1621            try:
1622                Histogram = Histograms[histogram]
1623            except KeyError:                       
1624                #skip if histogram not included e.g. in a sequential refinement
1625                continue
1626            hapData = HistoPhase[histogram]
1627            hId = Histogram['hId']
1628            if 'PWDR' in histogram:
1629                limits = Histogram['Limits'][1]
1630                inst = Histogram['Instrument Parameters'][0]
1631                Zero = inst['Zero'][1]
1632                if 'C' in inst['Type'][1]:
1633                    try:
1634                        wave = inst['Lam'][1]
1635                    except KeyError:
1636                        wave = inst['Lam1'][1]
1637                    dmin = wave/(2.0*sind(limits[1]/2.0))
1638                pfx = str(pId)+':'+str(hId)+':'
1639                for item in ['Scale','Extinction']:
1640                    hapDict[pfx+item] = hapData[item][0]
1641                    if hapData[item][1]:
1642                        hapVary.append(pfx+item)
1643                names = G2spc.HStrainNames(SGData)
1644                for i,name in enumerate(names):
1645                    hapDict[pfx+name] = hapData['HStrain'][0][i]
1646                    if hapData['HStrain'][1][i]:
1647                        hapVary.append(pfx+name)
1648                controlDict[pfx+'poType'] = hapData['Pref.Ori.'][0]
1649                if hapData['Pref.Ori.'][0] == 'MD':
1650                    hapDict[pfx+'MD'] = hapData['Pref.Ori.'][1]
1651                    controlDict[pfx+'MDAxis'] = hapData['Pref.Ori.'][3]
1652                    if hapData['Pref.Ori.'][2]:
1653                        hapVary.append(pfx+'MD')
1654                else:                           #'SH' spherical harmonics
1655                    controlDict[pfx+'SHord'] = hapData['Pref.Ori.'][4]
1656                    controlDict[pfx+'SHncof'] = len(hapData['Pref.Ori.'][5])
1657                    for item in hapData['Pref.Ori.'][5]:
1658                        hapDict[pfx+item] = hapData['Pref.Ori.'][5][item]
1659                        if hapData['Pref.Ori.'][2]:
1660                            hapVary.append(pfx+item)
1661                for item in ['Mustrain','Size']:
1662                    controlDict[pfx+item+'Type'] = hapData[item][0]
1663                    hapDict[pfx+item+';mx'] = hapData[item][1][2]
1664                    if hapData[item][2][2]:
1665                        hapVary.append(pfx+item+';mx')
1666                    if hapData[item][0] in ['isotropic','uniaxial']:
1667                        hapDict[pfx+item+';i'] = hapData[item][1][0]
1668                        if hapData[item][2][0]:
1669                            hapVary.append(pfx+item+';i')
1670                        if hapData[item][0] == 'uniaxial':
1671                            controlDict[pfx+item+'Axis'] = hapData[item][3]
1672                            hapDict[pfx+item+';a'] = hapData[item][1][1]
1673                            if hapData[item][2][1]:
1674                                hapVary.append(pfx+item+';a')
1675                    else:       #generalized for mustrain or ellipsoidal for size
1676                        Nterms = len(hapData[item][4])
1677                        if item == 'Mustrain':
1678                            names = G2spc.MustrainNames(SGData)
1679                            pwrs = []
1680                            for name in names:
1681                                h,k,l = name[1:]
1682                                pwrs.append([int(h),int(k),int(l)])
1683                            controlDict[pfx+'MuPwrs'] = pwrs
1684                        for i in range(Nterms):
1685                            sfx = ':'+str(i)
1686                            hapDict[pfx+item+sfx] = hapData[item][4][i]
1687                            if hapData[item][5][i]:
1688                                hapVary.append(pfx+item+sfx)
1689                for bab in ['BabA','BabU']:
1690                    hapDict[pfx+bab] = hapData['Babinet'][bab][0]
1691                    if hapData['Babinet'][bab][1]:
1692                        hapVary.append(pfx+bab)
1693                               
1694                if Print: 
1695                    print >>pFile,'\n Phase: ',phase,' in histogram: ',histogram
1696                    print >>pFile,135*'-'
1697                    print >>pFile,' Phase fraction  : %10.4f'%(hapData['Scale'][0]),' Refine?',hapData['Scale'][1]
1698                    print >>pFile,' Extinction coeff: %10.4f'%(hapData['Extinction'][0]),' Refine?',hapData['Extinction'][1]
1699                    if hapData['Pref.Ori.'][0] == 'MD':
1700                        Ax = hapData['Pref.Ori.'][3]
1701                        print >>pFile,' March-Dollase PO: %10.4f'%(hapData['Pref.Ori.'][1]),' Refine?',hapData['Pref.Ori.'][2], \
1702                            ' Axis: %d %d %d'%(Ax[0],Ax[1],Ax[2])
1703                    else: #'SH' for spherical harmonics
1704                        PrintSHPO(hapData['Pref.Ori.'])
1705                    PrintSize(hapData['Size'])
1706                    PrintMuStrain(hapData['Mustrain'],SGData)
1707                    PrintHStrain(hapData['HStrain'],SGData)
1708                    if hapData['Babinet']['BabA'][0]:
1709                        PrintBabinet(hapData['Babinet'])
1710                HKLd = np.array(G2lat.GenHLaue(dmin,SGData,A))
1711                refList = []
1712                for h,k,l,d in HKLd:
1713                    ext,mul,Uniq,phi = G2spc.GenHKLf([h,k,l],SGData)
1714                    mul *= 2      # for powder overlap of Friedel pairs
1715                    if ext:
1716                        continue
1717                    if 'C' in inst['Type'][0]:
1718                        pos = 2.0*asind(wave/(2.0*d))+Zero
1719                        if limits[0] < pos < limits[1]:
1720                            refList.append([h,k,l,mul,d,pos,0.0,0.0,0.0,0.0,0.0,Uniq,phi,0.0,{}])
1721                            #last item should contain form factor values by atom type
1722                    else:
1723                        raise ValueError 
1724                Histogram['Reflection Lists'][phase] = refList
1725            elif 'HKLF' in histogram:
1726                inst = Histogram['Instrument Parameters'][0]
1727                hId = Histogram['hId']
1728                hfx = ':%d:'%(hId)
1729                for item in inst:
1730                    hapDict[hfx+item] = inst[item][1]
1731                pfx = str(pId)+':'+str(hId)+':'
1732                hapDict[pfx+'Scale'] = hapData['Scale'][0]
1733                if hapData['Scale'][1]:
1734                    hapVary.append(pfx+'Scale')
1735                               
1736                extApprox,extType,extParms = hapData['Extinction']
1737                controlDict[pfx+'EType'] = extType
1738                controlDict[pfx+'EApprox'] = extApprox
1739                controlDict[pfx+'Tbar'] = extParms['Tbar']
1740                controlDict[pfx+'Cos2TM'] = extParms['Cos2TM']
1741                if 'Primary' in extType:
1742                    Ekey = ['Ep',]
1743                elif 'I & II' in extType:
1744                    Ekey = ['Eg','Es']
1745                elif 'Secondary Type II' == extType:
1746                    Ekey = ['Es',]
1747                elif 'Secondary Type I' == extType:
1748                    Ekey = ['Eg',]
1749                else:   #'None'
1750                    Ekey = []
1751                for eKey in Ekey:
1752                    hapDict[pfx+eKey] = extParms[eKey][0]
1753                    if extParms[eKey][1]:
1754                        hapVary.append(pfx+eKey)
1755                for bab in ['BabA','BabU']:
1756                    hapDict[pfx+bab] = hapData['Babinet'][bab][0]
1757                    if hapData['Babinet'][bab][1]:
1758                        hapVary.append(pfx+bab)
1759                if Print: 
1760                    print >>pFile,'\n Phase: ',phase,' in histogram: ',histogram
1761                    print >>pFile,135*'-'
1762                    print >>pFile,' Scale factor     : %10.4f'%(hapData['Scale'][0]),' Refine?',hapData['Scale'][1]
1763                    if extType != 'None':
1764                        print >>pFile,' Extinction  Type: %15s'%(extType),' approx: %10s'%(extApprox),' tbar: %6.3f'%(extParms['Tbar'])
1765                        text = ' Parameters       :'
1766                        for eKey in Ekey:
1767                            text += ' %4s : %10.3e Refine? '%(eKey,extParms[eKey][0])+str(extParms[eKey][1])
1768                        print >>pFile,text
1769                    if hapData['Babinet']['BabA'][0]:
1770                        PrintBabinet(hapData['Babinet'])
1771                Histogram['Reflection Lists'] = phase       
1772               
1773    return hapVary,hapDict,controlDict
1774   
1775def SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms,Print=True,pFile=None):
1776    'needs a doc string'
1777   
1778    def PrintSizeAndSig(hapData,sizeSig):
1779        line = '\n Size model:     %9s'%(hapData[0])
1780        refine = False
1781        if hapData[0] in ['isotropic','uniaxial']:
1782            line += ' equatorial:%12.4f'%(hapData[1][0])
1783            if sizeSig[0][0]:
1784                line += ', sig:%8.4f'%(sizeSig[0][0])
1785                refine = True
1786            if hapData[0] == 'uniaxial':
1787                line += ' axial:%12.4f'%(hapData[1][1])
1788                if sizeSig[0][1]:
1789                    refine = True
1790                    line += ', sig:%8.4f'%(sizeSig[0][1])
1791            line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
1792            if sizeSig[0][2]:
1793                refine = True
1794                line += ', sig:%8.4f'%(sizeSig[0][2])
1795            if refine:
1796                print >>pFile,line
1797        else:
1798            line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
1799            if sizeSig[0][2]:
1800                refine = True
1801                line += ', sig:%8.4f'%(sizeSig[0][2])
1802            Snames = ['S11','S22','S33','S12','S13','S23']
1803            ptlbls = ' name  :'
1804            ptstr =  ' value :'
1805            sigstr = ' sig   :'
1806            for i,name in enumerate(Snames):
1807                ptlbls += '%12s' % (name)
1808                ptstr += '%12.6f' % (hapData[4][i])
1809                if sizeSig[1][i]:
1810                    refine = True
1811                    sigstr += '%12.6f' % (sizeSig[1][i])
1812                else:
1813                    sigstr += 12*' '
1814            if refine:
1815                print >>pFile,line
1816                print >>pFile,ptlbls
1817                print >>pFile,ptstr
1818                print >>pFile,sigstr
1819       
1820    def PrintMuStrainAndSig(hapData,mustrainSig,SGData):
1821        line = '\n Mustrain model: %9s'%(hapData[0])
1822        refine = False
1823        if hapData[0] in ['isotropic','uniaxial']:
1824            line += ' equatorial:%12.1f'%(hapData[1][0])
1825            if mustrainSig[0][0]:
1826                line += ', sig:%8.1f'%(mustrainSig[0][0])
1827                refine = True
1828            if hapData[0] == 'uniaxial':
1829                line += ' axial:%12.1f'%(hapData[1][1])
1830                if mustrainSig[0][1]:
1831                     line += ', sig:%8.1f'%(mustrainSig[0][1])
1832            line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
1833            if mustrainSig[0][2]:
1834                refine = True
1835                line += ', sig:%8.3f'%(mustrainSig[0][2])
1836            if refine:
1837                print >>pFile,line
1838        else:
1839            line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
1840            if mustrainSig[0][2]:
1841                refine = True
1842                line += ', sig:%8.3f'%(mustrainSig[0][2])
1843            Snames = G2spc.MustrainNames(SGData)
1844            ptlbls = ' name  :'
1845            ptstr =  ' value :'
1846            sigstr = ' sig   :'
1847            for i,name in enumerate(Snames):
1848                ptlbls += '%12s' % (name)
1849                ptstr += '%12.6f' % (hapData[4][i])
1850                if mustrainSig[1][i]:
1851                    refine = True
1852                    sigstr += '%12.6f' % (mustrainSig[1][i])
1853                else:
1854                    sigstr += 12*' '
1855            if refine:
1856                print >>pFile,line
1857                print >>pFile,ptlbls
1858                print >>pFile,ptstr
1859                print >>pFile,sigstr
1860           
1861    def PrintHStrainAndSig(hapData,strainSig,SGData):
1862        Hsnames = G2spc.HStrainNames(SGData)
1863        ptlbls = ' name  :'
1864        ptstr =  ' value :'
1865        sigstr = ' sig   :'
1866        refine = False
1867        for i,name in enumerate(Hsnames):
1868            ptlbls += '%12s' % (name)
1869            ptstr += '%12.6g' % (hapData[0][i])
1870            if name in strainSig:
1871                refine = True
1872                sigstr += '%12.6g' % (strainSig[name])
1873            else:
1874                sigstr += 12*' '
1875        if refine:
1876            print >>pFile,'\n Hydrostatic/elastic strain: '
1877            print >>pFile,ptlbls
1878            print >>pFile,ptstr
1879            print >>pFile,sigstr
1880       
1881    def PrintSHPOAndSig(pfx,hapData,POsig):
1882        print >>pFile,'\n Spherical harmonics preferred orientation: Order:'+str(hapData[4])
1883        ptlbls = ' names :'
1884        ptstr =  ' values:'
1885        sigstr = ' sig   :'
1886        for item in hapData[5]:
1887            ptlbls += '%12s'%(item)
1888            ptstr += '%12.3f'%(hapData[5][item])
1889            if pfx+item in POsig:
1890                sigstr += '%12.3f'%(POsig[pfx+item])
1891            else:
1892                sigstr += 12*' ' 
1893        print >>pFile,ptlbls
1894        print >>pFile,ptstr
1895        print >>pFile,sigstr
1896       
1897    def PrintExtAndSig(pfx,hapData,ScalExtSig):
1898        print >>pFile,'\n Single crystal extinction: Type: ',hapData[0],' Approx: ',hapData[1]
1899        text = ''
1900        for item in hapData[2]:
1901            if pfx+item in ScalExtSig:
1902                text += '       %s: '%(item)
1903                text += '%12.2e'%(hapData[2][item][0])
1904                if pfx+item in ScalExtSig:
1905                    text += ' sig: %12.2e'%(ScalExtSig[pfx+item])
1906        print >>pFile,text       
1907       
1908    def PrintBabinetAndSig(pfx,hapData,BabSig):
1909        print >>pFile,'\n Babinet form factor modification: '
1910        ptlbls = ' names :'
1911        ptstr =  ' values:'
1912        sigstr = ' sig   :'
1913        for item in hapData:
1914            ptlbls += '%12s'%(item)
1915            ptstr += '%12.3f'%(hapData[item][0])
1916            if pfx+item in BabSig:
1917                sigstr += '%12.3f'%(BabSig[pfx+item])
1918            else:
1919                sigstr += 12*' ' 
1920        print >>pFile,ptlbls
1921        print >>pFile,ptstr
1922        print >>pFile,sigstr
1923   
1924    PhFrExtPOSig = {}
1925    SizeMuStrSig = {}
1926    ScalExtSig = {}
1927    BabSig = {}
1928    wtFrSum = {}
1929    for phase in Phases:
1930        HistoPhase = Phases[phase]['Histograms']
1931        General = Phases[phase]['General']
1932        SGData = General['SGData']
1933        pId = Phases[phase]['pId']
1934        histoList = HistoPhase.keys()
1935        histoList.sort()
1936        for histogram in histoList:
1937            try:
1938                Histogram = Histograms[histogram]
1939            except KeyError:                       
1940                #skip if histogram not included e.g. in a sequential refinement
1941                continue
1942            hapData = HistoPhase[histogram]
1943            hId = Histogram['hId']
1944            pfx = str(pId)+':'+str(hId)+':'
1945            if hId not in wtFrSum:
1946                wtFrSum[hId] = 0.
1947            if 'PWDR' in histogram:
1948                for item in ['Scale','Extinction']:
1949                    hapData[item][0] = parmDict[pfx+item]
1950                    if pfx+item in sigDict:
1951                        PhFrExtPOSig.update({pfx+item:sigDict[pfx+item],})
1952                wtFrSum[hId] += hapData['Scale'][0]*General['Mass']
1953                if hapData['Pref.Ori.'][0] == 'MD':
1954                    hapData['Pref.Ori.'][1] = parmDict[pfx+'MD']
1955                    if pfx+'MD' in sigDict:
1956                        PhFrExtPOSig.update({pfx+'MD':sigDict[pfx+'MD'],})
1957                else:                           #'SH' spherical harmonics
1958                    for item in hapData['Pref.Ori.'][5]:
1959                        hapData['Pref.Ori.'][5][item] = parmDict[pfx+item]
1960                        if pfx+item in sigDict:
1961                            PhFrExtPOSig.update({pfx+item:sigDict[pfx+item],})
1962                SizeMuStrSig.update({pfx+'Mustrain':[[0,0,0],[0 for i in range(len(hapData['Mustrain'][4]))]],
1963                    pfx+'Size':[[0,0,0],[0 for i in range(len(hapData['Size'][4]))]],
1964                    pfx+'HStrain':{}})                 
1965                for item in ['Mustrain','Size']:
1966                    hapData[item][1][2] = parmDict[pfx+item+';mx']
1967                    hapData[item][1][2] = min(1.,max(0.1,hapData[item][1][2]))
1968                    if pfx+item+';mx' in sigDict:
1969                        SizeMuStrSig[pfx+item][0][2] = sigDict[pfx+item+';mx']
1970                    if hapData[item][0] in ['isotropic','uniaxial']:                   
1971                        hapData[item][1][0] = parmDict[pfx+item+';i']
1972                        if item == 'Size':
1973                            hapData[item][1][0] = min(10.,max(0.001,hapData[item][1][0]))
1974                        if pfx+item+';i' in sigDict: 
1975                            SizeMuStrSig[pfx+item][0][0] = sigDict[pfx+item+';i']
1976                        if hapData[item][0] == 'uniaxial':
1977                            hapData[item][1][1] = parmDict[pfx+item+';a']
1978                            if item == 'Size':
1979                                hapData[item][1][1] = min(10.,max(0.001,hapData[item][1][1]))                       
1980                            if pfx+item+';a' in sigDict:
1981                                SizeMuStrSig[pfx+item][0][1] = sigDict[pfx+item+';a']
1982                    else:       #generalized for mustrain or ellipsoidal for size
1983                        Nterms = len(hapData[item][4])
1984                        for i in range(Nterms):
1985                            sfx = ':'+str(i)
1986                            hapData[item][4][i] = parmDict[pfx+item+sfx]
1987                            if pfx+item+sfx in sigDict:
1988                                SizeMuStrSig[pfx+item][1][i] = sigDict[pfx+item+sfx]
1989                names = G2spc.HStrainNames(SGData)
1990                for i,name in enumerate(names):
1991                    hapData['HStrain'][0][i] = parmDict[pfx+name]
1992                    if pfx+name in sigDict:
1993                        SizeMuStrSig[pfx+'HStrain'][name] = sigDict[pfx+name]
1994                for name in ['BabA','BabU']:
1995                    hapData['Babinet'][name][0] = parmDict[pfx+name]
1996                    if pfx+name in sigDict:
1997                        BabSig[pfx+name] = sigDict[pfx+name]               
1998               
1999            elif 'HKLF' in histogram:
2000                for item in ['Scale',]:
2001                    if parmDict.get(pfx+item):
2002                        hapData[item][0] = parmDict[pfx+item]
2003                        if pfx+item in sigDict:
2004                            ScalExtSig[pfx+item] = sigDict[pfx+item]
2005                for item in ['Ep','Eg','Es']:
2006                    if parmDict.get(pfx+item):
2007                        hapData['Extinction'][2][item][0] = parmDict[pfx+item]
2008                        if pfx+item in sigDict:
2009                            ScalExtSig[pfx+item] = sigDict[pfx+item]
2010                for name in ['BabA','BabU']:
2011                    hapData['Babinet'][name][0] = parmDict[pfx+name]
2012                    if pfx+name in sigDict:
2013                        BabSig[pfx+name] = sigDict[pfx+name]               
2014
2015    if Print:
2016        for phase in Phases:
2017            HistoPhase = Phases[phase]['Histograms']
2018            General = Phases[phase]['General']
2019            SGData = General['SGData']
2020            pId = Phases[phase]['pId']
2021            histoList = HistoPhase.keys()
2022            histoList.sort()
2023            for histogram in histoList:
2024                try:
2025                    Histogram = Histograms[histogram]
2026                except KeyError:                       
2027                    #skip if histogram not included e.g. in a sequential refinement
2028                    continue
2029                print >>pFile,'\n Phase: ',phase,' in histogram: ',histogram
2030                print >>pFile,130*'-'
2031                hapData = HistoPhase[histogram]
2032                hId = Histogram['hId']
2033                pfx = str(pId)+':'+str(hId)+':'
2034                if 'PWDR' in histogram:
2035                    print >>pFile,' Final refinement RF, RF^2 = %.2f%%, %.2f%% on %d reflections'   \
2036                        %(Histogram[pfx+'Rf'],Histogram[pfx+'Rf^2'],Histogram[pfx+'Nref'])
2037               
2038                    if pfx+'Scale' in PhFrExtPOSig:
2039                        wtFr = hapData['Scale'][0]*General['Mass']/wtFrSum[hId]
2040                        sigwtFr = PhFrExtPOSig[pfx+'Scale']*wtFr/hapData['Scale'][0]
2041                        print >>pFile,' Phase fraction  : %10.5f, sig %10.5f Weight fraction  : %8.5f, sig %10.5f' \
2042                            %(hapData['Scale'][0],PhFrExtPOSig[pfx+'Scale'],wtFr,sigwtFr)
2043                    if pfx+'Extinction' in PhFrExtPOSig:
2044                        print >>pFile,' Extinction coeff: %10.4f, sig %10.4f'%(hapData['Extinction'][0],PhFrExtPOSig[pfx+'Extinction'])
2045                    if hapData['Pref.Ori.'][0] == 'MD':
2046                        if pfx+'MD' in PhFrExtPOSig:
2047                            print >>pFile,' March-Dollase PO: %10.4f, sig %10.4f'%(hapData['Pref.Ori.'][1],PhFrExtPOSig[pfx+'MD'])
2048                    else:
2049                        PrintSHPOAndSig(pfx,hapData['Pref.Ori.'],PhFrExtPOSig)
2050                    PrintSizeAndSig(hapData['Size'],SizeMuStrSig[pfx+'Size'])
2051                    PrintMuStrainAndSig(hapData['Mustrain'],SizeMuStrSig[pfx+'Mustrain'],SGData)
2052                    PrintHStrainAndSig(hapData['HStrain'],SizeMuStrSig[pfx+'HStrain'],SGData)
2053                    if len(BabSig):
2054                        PrintBabinetAndSig(pfx,hapData['Babinet'],BabSig)
2055                   
2056                elif 'HKLF' in histogram:
2057                    print >>pFile,' Final refinement RF, RF^2 = %.2f%%, %.2f%% on %d reflections'   \
2058                        %(Histogram[pfx+'Rf'],Histogram[pfx+'Rf^2'],Histogram[pfx+'Nref'])
2059                    print >>pFile,' HKLF histogram weight factor = ','%.3f'%(Histogram['wtFactor'])
2060                    if pfx+'Scale' in ScalExtSig:
2061                        print >>pFile,' Scale factor : %10.4f, sig %10.4f'%(hapData['Scale'][0],ScalExtSig[pfx+'Scale'])
2062                    if hapData['Extinction'][0] != 'None':
2063                        PrintExtAndSig(pfx,hapData['Extinction'],ScalExtSig)
2064                    if len(BabSig):
2065                        PrintBabinetAndSig(pfx,hapData['Babinet'],BabSig)
2066
2067################################################################################
2068##### Histogram data
2069################################################################################       
2070                   
2071def GetHistogramData(Histograms,Print=True,pFile=None):
2072    'needs a doc string'
2073   
2074    def GetBackgroundParms(hId,Background):
2075        Back = Background[0]
2076        DebyePeaks = Background[1]
2077        bakType,bakFlag = Back[:2]
2078        backVals = Back[3:]
2079        backNames = [':'+str(hId)+':Back:'+str(i) for i in range(len(backVals))]
2080        backDict = dict(zip(backNames,backVals))
2081        backVary = []
2082        if bakFlag:
2083            backVary = backNames
2084        backDict[':'+str(hId)+':nDebye'] = DebyePeaks['nDebye']
2085        backDict[':'+str(hId)+':nPeaks'] = DebyePeaks['nPeaks']
2086        debyeDict = {}
2087        debyeList = []
2088        for i in range(DebyePeaks['nDebye']):
2089            debyeNames = [':'+str(hId)+':DebyeA:'+str(i),':'+str(hId)+':DebyeR:'+str(i),':'+str(hId)+':DebyeU:'+str(i)]
2090            debyeDict.update(dict(zip(debyeNames,DebyePeaks['debyeTerms'][i][::2])))
2091            debyeList += zip(debyeNames,DebyePeaks['debyeTerms'][i][1::2])
2092        debyeVary = []
2093        for item in debyeList:
2094            if item[1]:
2095                debyeVary.append(item[0])
2096        backDict.update(debyeDict)
2097        backVary += debyeVary
2098        peakDict = {}
2099        peakList = []
2100        for i in range(DebyePeaks['nPeaks']):
2101            peakNames = [':'+str(hId)+':BkPkpos:'+str(i),':'+str(hId)+ \
2102                ':BkPkint:'+str(i),':'+str(hId)+':BkPksig:'+str(i),':'+str(hId)+':BkPkgam:'+str(i)]
2103            peakDict.update(dict(zip(peakNames,DebyePeaks['peaksList'][i][::2])))
2104            peakList += zip(peakNames,DebyePeaks['peaksList'][i][1::2])
2105        peakVary = []
2106        for item in peakList:
2107            if item[1]:
2108                peakVary.append(item[0])
2109        backDict.update(peakDict)
2110        backVary += peakVary
2111        return bakType,backDict,backVary           
2112       
2113    def GetInstParms(hId,Inst):     
2114        dataType = Inst['Type'][0]
2115        instDict = {}
2116        insVary = []
2117        pfx = ':'+str(hId)+':'
2118        insKeys = Inst.keys()
2119        insKeys.sort()
2120        for item in insKeys:
2121            insName = pfx+item
2122            instDict[insName] = Inst[item][1]
2123            if Inst[item][2]:
2124                insVary.append(insName)
2125#        instDict[pfx+'X'] = max(instDict[pfx+'X'],0.001)
2126#        instDict[pfx+'Y'] = max(instDict[pfx+'Y'],0.001)
2127        instDict[pfx+'SH/L'] = max(instDict[pfx+'SH/L'],0.0005)
2128        return dataType,instDict,insVary
2129       
2130    def GetSampleParms(hId,Sample):
2131        sampVary = []
2132        hfx = ':'+str(hId)+':'       
2133        sampDict = {hfx+'Gonio. radius':Sample['Gonio. radius'],hfx+'Omega':Sample['Omega'],
2134            hfx+'Chi':Sample['Chi'],hfx+'Phi':Sample['Phi']}
2135        Type = Sample['Type']
2136        if 'Bragg' in Type:             #Bragg-Brentano
2137            for item in ['Scale','Shift','Transparency']:       #surface roughness?, diffuse scattering?
2138                sampDict[hfx+item] = Sample[item][0]
2139                if Sample[item][1]:
2140                    sampVary.append(hfx+item)
2141        elif 'Debye' in Type:        #Debye-Scherrer
2142            for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
2143                sampDict[hfx+item] = Sample[item][0]
2144                if Sample[item][1]:
2145                    sampVary.append(hfx+item)
2146        return Type,sampDict,sampVary
2147       
2148    def PrintBackground(Background):
2149        Back = Background[0]
2150        DebyePeaks = Background[1]
2151        print >>pFile,'\n Background function: ',Back[0],' Refine?',bool(Back[1])
2152        line = ' Coefficients: '
2153        for i,back in enumerate(Back[3:]):
2154            line += '%10.3f'%(back)
2155            if i and not i%10:
2156                line += '\n'+15*' '
2157        print >>pFile,line
2158        if DebyePeaks['nDebye']:
2159            print >>pFile,'\n Debye diffuse scattering coefficients'
2160            parms = ['DebyeA','DebyeR','DebyeU']
2161            line = ' names :  '
2162            for parm in parms:
2163                line += '%8s refine?'%(parm)
2164            print >>pFile,line
2165            for j,term in enumerate(DebyePeaks['debyeTerms']):
2166                line = ' term'+'%2d'%(j)+':'
2167                for i in range(3):
2168                    line += '%10.3f %5s'%(term[2*i],bool(term[2*i+1]))                   
2169                print >>pFile,line
2170        if DebyePeaks['nPeaks']:
2171            print >>pFile,'\n Single peak coefficients'
2172            parms =    ['BkPkpos','BkPkint','BkPksig','BkPkgam']
2173            line = ' names :  '
2174            for parm in parms:
2175                line += '%8s refine?'%(parm)
2176            print >>pFile,line
2177            for j,term in enumerate(DebyePeaks['peaksList']):
2178                line = ' peak'+'%2d'%(j)+':'
2179                for i in range(4):
2180                    line += '%10.3f %5s'%(term[2*i],bool(term[2*i+1]))                   
2181                print >>pFile,line
2182       
2183    def PrintInstParms(Inst):
2184        print >>pFile,'\n Instrument Parameters:'
2185        ptlbls = ' name  :'
2186        ptstr =  ' value :'
2187        varstr = ' refine:'
2188        insKeys = Inst.keys()
2189        insKeys.sort()
2190        for item in insKeys:
2191            if item != 'Type':
2192                ptlbls += '%12s' % (item)
2193                ptstr += '%12.6f' % (Inst[item][1])
2194                if item in ['Lam1','Lam2','Azimuth']:
2195                    varstr += 12*' '
2196                else:
2197                    varstr += '%12s' % (str(bool(Inst[item][2])))
2198        print >>pFile,ptlbls
2199        print >>pFile,ptstr
2200        print >>pFile,varstr
2201       
2202    def PrintSampleParms(Sample):
2203        print >>pFile,'\n Sample Parameters:'
2204        print >>pFile,' Goniometer omega = %.2f, chi = %.2f, phi = %.2f'% \
2205            (Sample['Omega'],Sample['Chi'],Sample['Phi'])
2206        ptlbls = ' name  :'
2207        ptstr =  ' value :'
2208        varstr = ' refine:'
2209        if 'Bragg' in Sample['Type']:
2210            for item in ['Scale','Shift','Transparency']:
2211                ptlbls += '%14s'%(item)
2212                ptstr += '%14.4f'%(Sample[item][0])
2213                varstr += '%14s'%(str(bool(Sample[item][1])))
2214           
2215        elif 'Debye' in Type:        #Debye-Scherrer
2216            for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
2217                ptlbls += '%14s'%(item)
2218                ptstr += '%14.4f'%(Sample[item][0])
2219                varstr += '%14s'%(str(bool(Sample[item][1])))
2220
2221        print >>pFile,ptlbls
2222        print >>pFile,ptstr
2223        print >>pFile,varstr
2224       
2225    histDict = {}
2226    histVary = []
2227    controlDict = {}
2228    histoList = Histograms.keys()
2229    histoList.sort()
2230    for histogram in histoList:
2231        if 'PWDR' in histogram:
2232            Histogram = Histograms[histogram]
2233            hId = Histogram['hId']
2234            pfx = ':'+str(hId)+':'
2235            controlDict[pfx+'wtFactor'] = Histogram['wtFactor']
2236            controlDict[pfx+'Limits'] = Histogram['Limits'][1]
2237           
2238            Background = Histogram['Background']
2239            Type,bakDict,bakVary = GetBackgroundParms(hId,Background)
2240            controlDict[pfx+'bakType'] = Type
2241            histDict.update(bakDict)
2242            histVary += bakVary
2243           
2244            Inst = Histogram['Instrument Parameters'][0]
2245            Type,instDict,insVary = GetInstParms(hId,Inst)
2246            controlDict[pfx+'histType'] = Type
2247            if pfx+'Lam1' in instDict:
2248                controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam1']
2249            else:
2250                controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam']           
2251            histDict.update(instDict)
2252            histVary += insVary
2253           
2254            Sample = Histogram['Sample Parameters']
2255            Type,sampDict,sampVary = GetSampleParms(hId,Sample)
2256            controlDict[pfx+'instType'] = Type
2257            histDict.update(sampDict)
2258            histVary += sampVary
2259           
2260   
2261            if Print: 
2262                print >>pFile,'\n Histogram: ',histogram,' histogram Id: ',hId
2263                print >>pFile,135*'-'
2264                Units = {'C':' deg','T':' msec'}
2265                units = Units[controlDict[pfx+'histType'][2]]
2266                Limits = controlDict[pfx+'Limits']
2267                print >>pFile,' Instrument type: ',Sample['Type']
2268                print >>pFile,' Histogram limits: %8.2f%s to %8.2f%s'%(Limits[0],units,Limits[1],units)     
2269                PrintSampleParms(Sample)
2270                PrintInstParms(Inst)
2271                PrintBackground(Background)
2272        elif 'HKLF' in histogram:
2273            Histogram = Histograms[histogram]
2274            hId = Histogram['hId']
2275            pfx = ':'+str(hId)+':'
2276            controlDict[pfx+'wtFactor'] =Histogram['wtFactor']
2277            Inst = Histogram['Instrument Parameters'][0]
2278            controlDict[pfx+'histType'] = Inst['Type'][0]
2279            histDict[pfx+'Lam'] = Inst['Lam'][1]
2280            controlDict[pfx+'keV'] = 12.397639/histDict[pfx+'Lam']                   
2281    return histVary,histDict,controlDict
2282   
2283def SetHistogramData(parmDict,sigDict,Histograms,Print=True,pFile=None):
2284    'needs a doc string'
2285   
2286    def SetBackgroundParms(pfx,Background,parmDict,sigDict):
2287        Back = Background[0]
2288        DebyePeaks = Background[1]
2289        lenBack = len(Back[3:])
2290        backSig = [0 for i in range(lenBack+3*DebyePeaks['nDebye']+4*DebyePeaks['nPeaks'])]
2291        for i in range(lenBack):
2292            Back[3+i] = parmDict[pfx+'Back:'+str(i)]
2293            if pfx+'Back:'+str(i) in sigDict:
2294                backSig[i] = sigDict[pfx+'Back:'+str(i)]
2295        if DebyePeaks['nDebye']:
2296            for i in range(DebyePeaks['nDebye']):
2297                names = [pfx+'DebyeA:'+str(i),pfx+'DebyeR:'+str(i),pfx+'DebyeU:'+str(i)]
2298                for j,name in enumerate(names):
2299                    DebyePeaks['debyeTerms'][i][2*j] = parmDict[name]
2300                    if name in sigDict:
2301                        backSig[lenBack+3*i+j] = sigDict[name]           
2302        if DebyePeaks['nPeaks']:
2303            for i in range(DebyePeaks['nPeaks']):
2304                names = [pfx+'BkPkpos:'+str(i),pfx+'BkPkint:'+str(i),
2305                    pfx+'BkPksig:'+str(i),pfx+'BkPkgam:'+str(i)]
2306                for j,name in enumerate(names):
2307                    DebyePeaks['peaksList'][i][2*j] = parmDict[name]
2308                    if name in sigDict:
2309                        backSig[lenBack+3*DebyePeaks['nDebye']+4*i+j] = sigDict[name]
2310        return backSig
2311       
2312    def SetInstParms(pfx,Inst,parmDict,sigDict):
2313        instSig = {}
2314        insKeys = Inst.keys()
2315        insKeys.sort()
2316        for item in insKeys:
2317            insName = pfx+item
2318            Inst[item][1] = parmDict[insName]
2319            if insName in sigDict:
2320                instSig[item] = sigDict[insName]
2321            else:
2322                instSig[item] = 0
2323        return instSig
2324       
2325    def SetSampleParms(pfx,Sample,parmDict,sigDict):
2326        if 'Bragg' in Sample['Type']:             #Bragg-Brentano
2327            sampSig = [0 for i in range(3)]
2328            for i,item in enumerate(['Scale','Shift','Transparency']):       #surface roughness?, diffuse scattering?
2329                Sample[item][0] = parmDict[pfx+item]
2330                if pfx+item in sigDict:
2331                    sampSig[i] = sigDict[pfx+item]
2332        elif 'Debye' in Sample['Type']:        #Debye-Scherrer
2333            sampSig = [0 for i in range(4)]
2334            for i,item in enumerate(['Scale','Absorption','DisplaceX','DisplaceY']):
2335                Sample[item][0] = parmDict[pfx+item]
2336                if pfx+item in sigDict:
2337                    sampSig[i] = sigDict[pfx+item]
2338        return sampSig
2339       
2340    def PrintBackgroundSig(Background,backSig):
2341        Back = Background[0]
2342        DebyePeaks = Background[1]
2343        lenBack = len(Back[3:])
2344        valstr = ' value : '
2345        sigstr = ' sig   : '
2346        refine = False
2347        for i,back in enumerate(Back[3:]):
2348            valstr += '%10.4g'%(back)
2349            if Back[1]:
2350                refine = True
2351                sigstr += '%10.4g'%(backSig[i])
2352            else:
2353                sigstr += 10*' '
2354        if refine:
2355            print >>pFile,'\n Background function: ',Back[0]
2356            print >>pFile,valstr
2357            print >>pFile,sigstr
2358        if DebyePeaks['nDebye']:
2359            ifAny = False
2360            ptfmt = "%12.3f"
2361            names =  ' names :'
2362            ptstr =  ' values:'
2363            sigstr = ' esds  :'
2364            for item in sigDict:
2365                if 'Debye' in item:
2366                    ifAny = True
2367                    names += '%12s'%(item)
2368                    ptstr += ptfmt%(parmDict[item])
2369                    sigstr += ptfmt%(sigDict[item])
2370            if ifAny:
2371                print >>pFile,'\n Debye diffuse scattering coefficients'
2372                print >>pFile,names
2373                print >>pFile,ptstr
2374                print >>pFile,sigstr
2375        if DebyePeaks['nPeaks']:
2376            ifAny = False
2377            ptfmt = "%14.3f"
2378            names =  ' names :'
2379            ptstr =  ' values:'
2380            sigstr = ' esds  :'
2381            for item in sigDict:
2382                if 'BkPk' in item:
2383                    ifAny = True
2384                    names += '%14s'%(item)
2385                    ptstr += ptfmt%(parmDict[item])
2386                    sigstr += ptfmt%(sigDict[item])
2387            if ifAny:
2388                print >>pFile,'\n Single peak coefficients'
2389                print >>pFile,names
2390                print >>pFile,ptstr
2391                print >>pFile,sigstr
2392       
2393    def PrintInstParmsSig(Inst,instSig):
2394        ptlbls = ' names :'
2395        ptstr =  ' value :'
2396        sigstr = ' sig   :'
2397        refine = False
2398        insKeys = instSig.keys()
2399        insKeys.sort()
2400        for name in insKeys:
2401            if name not in  ['Type','Lam1','Lam2','Azimuth']:
2402                ptlbls += '%12s' % (name)
2403                ptstr += '%12.6f' % (Inst[name][1])
2404                if instSig[name]:
2405                    refine = True
2406                    sigstr += '%12.6f' % (instSig[name])
2407                else:
2408                    sigstr += 12*' '
2409        if refine:
2410            print >>pFile,'\n Instrument Parameters:'
2411            print >>pFile,ptlbls
2412            print >>pFile,ptstr
2413            print >>pFile,sigstr
2414       
2415    def PrintSampleParmsSig(Sample,sampleSig):
2416        ptlbls = ' names :'
2417        ptstr =  ' values:'
2418        sigstr = ' sig   :'
2419        refine = False
2420        if 'Bragg' in Sample['Type']:
2421            for i,item in enumerate(['Scale','Shift','Transparency']):
2422                ptlbls += '%14s'%(item)
2423                ptstr += '%14.4f'%(Sample[item][0])
2424                if sampleSig[i]:
2425                    refine = True
2426                    sigstr += '%14.4f'%(sampleSig[i])
2427                else:
2428                    sigstr += 14*' '
2429           
2430        elif 'Debye' in Sample['Type']:        #Debye-Scherrer
2431            for i,item in enumerate(['Scale','Absorption','DisplaceX','DisplaceY']):
2432                ptlbls += '%14s'%(item)
2433                ptstr += '%14.4f'%(Sample[item][0])
2434                if sampleSig[i]:
2435                    refine = True
2436                    sigstr += '%14.4f'%(sampleSig[i])
2437                else:
2438                    sigstr += 14*' '
2439
2440        if refine:
2441            print >>pFile,'\n Sample Parameters:'
2442            print >>pFile,ptlbls
2443            print >>pFile,ptstr
2444            print >>pFile,sigstr
2445       
2446    histoList = Histograms.keys()
2447    histoList.sort()
2448    for histogram in histoList:
2449        if 'PWDR' in histogram:
2450            Histogram = Histograms[histogram]
2451            hId = Histogram['hId']
2452            pfx = ':'+str(hId)+':'
2453            Background = Histogram['Background']
2454            backSig = SetBackgroundParms(pfx,Background,parmDict,sigDict)
2455           
2456            Inst = Histogram['Instrument Parameters'][0]
2457            instSig = SetInstParms(pfx,Inst,parmDict,sigDict)
2458       
2459            Sample = Histogram['Sample Parameters']
2460            sampSig = SetSampleParms(pfx,Sample,parmDict,sigDict)
2461
2462            print >>pFile,'\n Histogram: ',histogram,' histogram Id: ',hId
2463            print >>pFile,135*'-'
2464            print >>pFile,' Final refinement wR = %.2f%% on %d observations in this histogram'%(Histogram['wR'],Histogram['Nobs'])
2465            print >>pFile,' PWDR histogram weight factor = '+'%.3f'%(Histogram['wtFactor'])
2466            if Print:
2467                print >>pFile,' Instrument type: ',Sample['Type']
2468                PrintSampleParmsSig(Sample,sampSig)
2469                PrintInstParmsSig(Inst,instSig)
2470                PrintBackgroundSig(Background,backSig)
2471               
Note: See TracBrowser for help on using the repository browser.