source: trunk/GSASIIstrIO.py @ 942

Last change on this file since 942 was 942, checked in by vondreele, 9 years ago

mods for MC/SA:
moved scat fac routines from GSASIIstrIO.py & GSASIIstrMath.py to GSASIIElem.py

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