source: trunk/GSASIIstrIO.py @ 946

Last change on this file since 946 was 946, checked in by toby, 9 years ago

update self-docs, start work on constraints object

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