source: trunk/GSASIIstrIO.py @ 939

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

fix & cleanup unit tests; add/change doc strings for sphinx; add all G2 py files to sphinx

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