source: trunk/GSASIIstruct.py @ 889

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

add torsion angle derivatives in RB
small fixes to other RB stuff

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 218.4 KB
Line 
1# -*- coding: utf-8 -*-
2#GSASIIstructure - structure computation routines
3########### SVN repository information ###################
4# $Date: 2013-04-19 20:38:07 +0000 (Fri, 19 Apr 2013) $
5# $Author: vondreele $
6# $Revision: 889 $
7# $URL: trunk/GSASIIstruct.py $
8# $Id: GSASIIstruct.py 889 2013-04-19 20:38:07Z vondreele $
9########### SVN repository information ###################
10import sys
11import os
12import os.path as ospath
13import time
14import math
15import copy
16import random
17import cPickle
18import numpy as np
19import numpy.ma as ma
20import numpy.linalg as nl
21import scipy.optimize as so
22import GSASIIpath
23GSASIIpath.SetVersionNumber("$Revision: 889 $")
24import GSASIIElem as G2el
25import GSASIIlattice as G2lat
26import GSASIIspc as G2spc
27import GSASIIpwd as G2pwd
28import GSASIImapvars as G2mv
29import GSASIImath as G2mth
30
31sind = lambda x: np.sin(x*np.pi/180.)
32cosd = lambda x: np.cos(x*np.pi/180.)
33tand = lambda x: np.tan(x*np.pi/180.)
34asind = lambda x: 180.*np.arcsin(x)/np.pi
35acosd = lambda x: 180.*np.arccos(x)/np.pi
36atan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
37   
38ateln2 = 8.0*math.log(2.0)
39DEBUG = False       #only for powders!
40
41def GetControls(GPXfile):
42    ''' Returns dictionary of control items found in GSASII gpx file
43    input:
44        GPXfile = .gpx full file name
45    return:
46        Controls = dictionary of control items
47    '''
48    Controls = {'deriv type':'analytic Hessian','max cyc':3,'max Hprocess':1,
49        'max Rprocess':1,'min dM/M':0.0001,'shift factor':1.}
50    fl = open(GPXfile,'rb')
51    while True:
52        try:
53            data = cPickle.load(fl)
54        except EOFError:
55            break
56        datum = data[0]
57        if datum[0] == 'Controls':
58            Controls.update(datum[1])
59    fl.close()
60    return Controls
61   
62def GetConstraints(GPXfile):
63    '''Read the constraints from the GPX file and interpret them
64    '''
65    constList = []
66    fl = open(GPXfile,'rb')
67    while True:
68        try:
69            data = cPickle.load(fl)
70        except EOFError:
71            break
72        datum = data[0]
73        if datum[0] == 'Constraints':
74            constDict = datum[1]
75            for item in constDict:
76                constList += constDict[item]
77    fl.close()
78    constDict,fixedList,ignored = ProcessConstraints(constList)
79    if ignored:
80        print ignored,'old-style Constraints were rejected'
81    return constDict,fixedList
82   
83def ProcessConstraints(constList):
84    "interpret constraints"
85    constDict = []
86    fixedList = []
87    ignored = 0
88    for item in constList:
89        if item[-1] == 'h':
90            # process a hold
91            fixedList.append('0')
92            constDict.append({item[0][1]:0.0})
93        elif item[-1] == 'f':
94            # process a new variable
95            fixedList.append(None)
96            constDict.append({})
97            for term in item[:-3]:
98                constDict[-1][term[1]] = term[0]
99            #constFlag[-1] = ['Vary']
100        elif item[-1] == 'c': 
101            # process a contraint relationship
102            fixedList.append(str(item[-3]))
103            constDict.append({})
104            for term in item[:-3]:
105                constDict[-1][term[1]] = term[0]
106            #constFlag[-1] = ['VaryFree']
107        elif item[-1] == 'e':
108            # process an equivalence
109            firstmult = None
110            eqlist = []
111            for term in item[:-3]:
112                if term[0] == 0: term[0] = 1.0
113                if firstmult is None:
114                    firstmult,firstvar = term
115                else:
116                    eqlist.append([term[1],firstmult/term[0]])
117            G2mv.StoreEquivalence(firstvar,eqlist)
118        else:
119            ignored += 1
120    return constDict,fixedList,ignored
121
122def CheckConstraints(GPXfile):
123    '''Load constraints and related info and return any error or warning messages'''
124    # init constraints
125    G2mv.InitVars()   
126    # get variables
127    Histograms,Phases = GetUsedHistogramsAndPhases(GPXfile)
128    if not Phases:
129        return 'Error: No Phases!',''
130    if not Histograms:
131        return 'Error: no diffraction data',''
132    rigidbodyDict = GetRigidBodies(GPXfile)
133    rbIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]})
134    rbVary,rbDict = GetRigidBodyModels(rigidbodyDict,Print=False)
135    Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables = GetPhaseData(Phases,RestraintDict=None,rbIds=rbIds,Print=False)
136    hapVary,hapDict,controlDict = GetHistogramPhaseData(Phases,Histograms,Print=False)
137    histVary,histDict,controlDict = GetHistogramData(Histograms,Print=False)
138    varyList = rbVary+phaseVary+hapVary+histVary
139    constrDict,fixedList = GetConstraints(GPXfile)
140    return G2mv.CheckConstraints(varyList,constrDict,fixedList)
141   
142def GetRestraints(GPXfile):
143    '''Read the restraints from the GPX file
144    '''
145    fl = open(GPXfile,'rb')
146    while True:
147        try:
148            data = cPickle.load(fl)
149        except EOFError:
150            break
151        datum = data[0]
152        if datum[0] == 'Restraints':
153            restraintDict = datum[1]
154    fl.close()
155    return restraintDict
156   
157def GetRigidBodies(GPXfile):
158    '''Read the rigid body models from the GPX file
159    '''
160    fl = open(GPXfile,'rb')
161    while True:
162        try:
163            data = cPickle.load(fl)
164        except EOFError:
165            break
166        datum = data[0]
167        if datum[0] == 'Rigid bodies':
168            rigidbodyDict = datum[1]
169    fl.close()
170    return rigidbodyDict
171       
172def GetPhaseNames(GPXfile):
173    ''' Returns a list of phase names found under 'Phases' in GSASII gpx file
174    input:
175        GPXfile = gpx full file name
176    return:
177        PhaseNames = list of phase names
178    '''
179    fl = open(GPXfile,'rb')
180    PhaseNames = []
181    while True:
182        try:
183            data = cPickle.load(fl)
184        except EOFError:
185            break
186        datum = data[0]
187        if 'Phases' == datum[0]:
188            for datus in data[1:]:
189                PhaseNames.append(datus[0])
190    fl.close()
191    return PhaseNames
192
193def GetAllPhaseData(GPXfile,PhaseName):
194    ''' Returns the entire dictionary for PhaseName from GSASII gpx file
195    input:
196        GPXfile = gpx full file name
197        PhaseName = phase name
198    return:
199        phase dictionary
200    '''       
201    fl = open(GPXfile,'rb')
202    General = {}
203    Atoms = []
204    while True:
205        try:
206            data = cPickle.load(fl)
207        except EOFError:
208            break
209        datum = data[0]
210        if 'Phases' == datum[0]:
211            for datus in data[1:]:
212                if datus[0] == PhaseName:
213                    break
214    fl.close()
215    return datus[1]
216   
217def GetHistograms(GPXfile,hNames):
218    """ Returns a dictionary of histograms found in GSASII gpx file
219    input:
220        GPXfile = .gpx full file name
221        hNames = list of histogram names
222    return:
223        Histograms = dictionary of histograms (types = PWDR & HKLF)
224    """
225    fl = open(GPXfile,'rb')
226    Histograms = {}
227    while True:
228        try:
229            data = cPickle.load(fl)
230        except EOFError:
231            break
232        datum = data[0]
233        hist = datum[0]
234        if hist in hNames:
235            if 'PWDR' in hist[:4]:
236                PWDRdata = {}
237                try:
238                    PWDRdata.update(datum[1][0])        #weight factor
239                except ValueError:
240                    PWDRdata['wtFactor'] = 1.0          #patch
241                PWDRdata['Data'] = datum[1][1]          #powder data arrays
242                PWDRdata[data[2][0]] = data[2][1]       #Limits
243                PWDRdata[data[3][0]] = data[3][1]       #Background
244                PWDRdata[data[4][0]] = data[4][1]       #Instrument parameters
245                PWDRdata[data[5][0]] = data[5][1]       #Sample parameters
246                try:
247                    PWDRdata[data[9][0]] = data[9][1]       #Reflection lists might be missing
248                except IndexError:
249                    PWDRdata['Reflection Lists'] = {}
250   
251                Histograms[hist] = PWDRdata
252            elif 'HKLF' in hist[:4]:
253                HKLFdata = {}
254                try:
255                    HKLFdata.update(datum[1][0])        #weight factor
256                except ValueError:
257                    HKLFdata['wtFactor'] = 1.0          #patch
258                HKLFdata['Data'] = datum[1][1]
259                HKLFdata[data[1][0]] = data[1][1]       #Instrument parameters
260                HKLFdata['Reflection Lists'] = None
261                Histograms[hist] = HKLFdata           
262    fl.close()
263    return Histograms
264   
265def GetHistogramNames(GPXfile,hType):
266    """ Returns a list of histogram names found in GSASII gpx file
267    input:
268        GPXfile = .gpx full file name
269        hType = list ['PWDR','HKLF']
270    return:
271        HistogramNames = list of histogram names (types = PWDR & HKLF)
272    """
273    fl = open(GPXfile,'rb')
274    HistogramNames = []
275    while True:
276        try:
277            data = cPickle.load(fl)
278        except EOFError:
279            break
280        datum = data[0]
281        if datum[0][:4] in hType:
282            HistogramNames.append(datum[0])
283    fl.close()
284    return HistogramNames
285   
286def GetUsedHistogramsAndPhases(GPXfile):
287    ''' Returns all histograms that are found in any phase
288    and any phase that uses a histogram
289    input:
290        GPXfile = .gpx full file name
291    return:
292        Histograms = dictionary of histograms as {name:data,...}
293        Phases = dictionary of phases that use histograms
294    '''
295    phaseNames = GetPhaseNames(GPXfile)
296    histoList = GetHistogramNames(GPXfile,['PWDR','HKLF'])
297    allHistograms = GetHistograms(GPXfile,histoList)
298    phaseData = {}
299    for name in phaseNames: 
300        phaseData[name] =  GetAllPhaseData(GPXfile,name)
301    Histograms = {}
302    Phases = {}
303    for phase in phaseData:
304        Phase = phaseData[phase]
305        if Phase['Histograms']:
306            if phase not in Phases:
307                pId = phaseNames.index(phase)
308                Phase['pId'] = pId
309                Phases[phase] = Phase
310            for hist in Phase['Histograms']:
311                if 'Use' not in Phase['Histograms'][hist]:      #patch
312                    Phase['Histograms'][hist]['Use'] = True         
313                if hist not in Histograms and Phase['Histograms'][hist]['Use']:
314                    Histograms[hist] = allHistograms[hist]
315                    hId = histoList.index(hist)
316                    Histograms[hist]['hId'] = hId
317    return Histograms,Phases
318   
319def getBackupName(GPXfile,makeBack):
320    GPXpath,GPXname = ospath.split(GPXfile)
321    if GPXpath == '': GPXpath = '.'
322    Name = ospath.splitext(GPXname)[0]
323    files = os.listdir(GPXpath)
324    last = 0
325    for name in files:
326        name = name.split('.')
327        if len(name) == 3 and name[0] == Name and 'bak' in name[1]:
328            if makeBack:
329                last = max(last,int(name[1].strip('bak'))+1)
330            else:
331                last = max(last,int(name[1].strip('bak')))
332    GPXback = ospath.join(GPXpath,ospath.splitext(GPXname)[0]+'.bak'+str(last)+'.gpx')
333    return GPXback   
334       
335def GPXBackup(GPXfile,makeBack=True):
336    import distutils.file_util as dfu
337    GPXback = getBackupName(GPXfile,makeBack)
338    dfu.copy_file(GPXfile,GPXback)
339    return GPXback
340
341def SetUsedHistogramsAndPhases(GPXfile,Histograms,Phases,RigidBodies,CovData,makeBack=True):
342    ''' Updates gpxfile from all histograms that are found in any phase
343    and any phase that used a histogram. Also updates rigid body definitions.
344    input:
345        GPXfile = .gpx full file name
346        Histograms = dictionary of histograms as {name:data,...}
347        Phases = dictionary of phases that use histograms
348        RigidBodies = dictionary of rigid bodies
349        CovData = dictionary of refined variables, varyList, & covariance matrix
350        makeBack = True if new backup of .gpx file is to be made; else use the last one made
351    '''
352                       
353    GPXback = GPXBackup(GPXfile,makeBack)
354    print 'Read from file:',GPXback
355    print 'Save to file  :',GPXfile
356    infile = open(GPXback,'rb')
357    outfile = open(GPXfile,'wb')
358    while True:
359        try:
360            data = cPickle.load(infile)
361        except EOFError:
362            break
363        datum = data[0]
364#        print 'read: ',datum[0]
365        if datum[0] == 'Phases':
366            for iphase in range(len(data)):
367                if data[iphase][0] in Phases:
368                    phaseName = data[iphase][0]
369                    data[iphase][1].update(Phases[phaseName])
370        elif datum[0] == 'Covariance':
371            data[0][1] = CovData
372        elif datum[0] == 'Rigid bodies':
373            data[0][1] = RigidBodies
374        try:
375            histogram = Histograms[datum[0]]
376#            print 'found ',datum[0]
377            data[0][1][1] = histogram['Data']
378            for datus in data[1:]:
379#                print '    read: ',datus[0]
380                if datus[0] in ['Background','Instrument Parameters','Sample Parameters','Reflection Lists']:
381                    datus[1] = histogram[datus[0]]
382        except KeyError:
383            pass
384                               
385        cPickle.dump(data,outfile,1)
386    infile.close()
387    outfile.close()
388    print 'GPX file save successful'
389   
390def SetSeqResult(GPXfile,Histograms,SeqResult):
391    GPXback = GPXBackup(GPXfile)
392    print 'Read from file:',GPXback
393    print 'Save to file  :',GPXfile
394    infile = open(GPXback,'rb')
395    outfile = open(GPXfile,'wb')
396    while True:
397        try:
398            data = cPickle.load(infile)
399        except EOFError:
400            break
401        datum = data[0]
402        if datum[0] == 'Sequential results':
403            data[0][1] = SeqResult
404        try:
405            histogram = Histograms[datum[0]]
406            data[0][1][1] = histogram['Data']
407            for datus in data[1:]:
408                if datus[0] in ['Background','Instrument Parameters','Sample Parameters','Reflection Lists']:
409                    datus[1] = histogram[datus[0]]
410        except KeyError:
411            pass
412                               
413        cPickle.dump(data,outfile,1)
414    infile.close()
415    outfile.close()
416    print 'GPX file save successful'
417                       
418def ShowBanner(pFile=None):
419    print >>pFile,80*'*'
420    print >>pFile,'   General Structure Analysis System-II Crystal Structure Refinement'
421    print >>pFile,'              by Robert B. Von Dreele & Brian H. Toby'
422    print >>pFile,'                Argonne National Laboratory(C), 2010'
423    print >>pFile,' This product includes software developed by the UChicago Argonne, LLC,' 
424    print >>pFile,'            as Operator of Argonne National Laboratory.'
425    print >>pFile,'                          Please cite:'
426    print >>pFile,'   B.H. Toby & R.B. Von Dreele, J. Appl. Cryst. 46, 544-549 (2013)'
427
428    print >>pFile,80*'*','\n'
429
430def ShowControls(Controls,pFile=None):
431    print >>pFile,' Least squares controls:'
432    print >>pFile,' Refinement type: ',Controls['deriv type']
433    if 'Hessian' in Controls['deriv type']:
434        print >>pFile,' Maximum number of cycles:',Controls['max cyc']
435    else:
436        print >>pFile,' Minimum delta-M/M for convergence: ','%.2g'%(Controls['min dM/M'])
437    print >>pFile,' Initial shift factor: ','%.3f'%(Controls['shift factor'])
438   
439def GetFFtable(General):
440    ''' returns a dictionary of form factor data for atom types found in General
441    input:
442        General = dictionary of phase info.; includes AtomTypes
443    return:
444        FFtable = dictionary of form factor data; key is atom type
445    '''
446    atomTypes = General['AtomTypes']
447    FFtable = {}
448    for El in atomTypes:
449        FFs = G2el.GetFormFactorCoeff(El.split('+')[0].split('-')[0])
450        for item in FFs:
451            if item['Symbol'] == El.upper():
452                FFtable[El] = item
453    return FFtable
454   
455def GetBLtable(General):
456    ''' returns a dictionary of neutron scattering length data for atom types & isotopes found in General
457    input:
458        General = dictionary of phase info.; includes AtomTypes & Isotopes
459    return:
460        BLtable = dictionary of scattering length data; key is atom type
461    '''
462    atomTypes = General['AtomTypes']
463    BLtable = {}
464    isotopes = General['Isotopes']
465    isotope = General['Isotope']
466    for El in atomTypes:
467        BLtable[El] = [isotope[El],isotopes[El][isotope[El]]]
468    return BLtable
469       
470def GetPawleyConstr(SGLaue,PawleyRef,pawleyVary):
471#    if SGLaue in ['-1','2/m','mmm']:
472#        return                      #no Pawley symmetry required constraints
473    eqvDict = {}
474    for i,varyI in enumerate(pawleyVary):
475        eqvDict[varyI] = []
476        refI = int(varyI.split(':')[-1])
477        ih,ik,il = PawleyRef[refI][:3]
478        dspI = PawleyRef[refI][4]
479        for varyJ in pawleyVary[i+1:]:
480            refJ = int(varyJ.split(':')[-1])
481            jh,jk,jl = PawleyRef[refJ][:3]
482            dspJ = PawleyRef[refJ][4]
483            if SGLaue in ['4/m','4/mmm']:
484                isum = ih**2+ik**2
485                jsum = jh**2+jk**2
486                if abs(il) == abs(jl) and isum == jsum:
487                    eqvDict[varyI].append(varyJ) 
488            elif SGLaue in ['3R','3mR']:
489                isum = ih**2+ik**2+il**2
490                jsum = jh**2+jk**2*jl**2
491                isum2 = ih*ik+ih*il+ik*il
492                jsum2 = jh*jk+jh*jl+jk*jl
493                if isum == jsum and isum2 == jsum2:
494                    eqvDict[varyI].append(varyJ) 
495            elif SGLaue in ['3','3m1','31m','6/m','6/mmm']:
496                isum = ih**2+ik**2+ih*ik
497                jsum = jh**2+jk**2+jh*jk
498                if abs(il) == abs(jl) and isum == jsum:
499                    eqvDict[varyI].append(varyJ) 
500            elif SGLaue in ['m3','m3m']:
501                isum = ih**2+ik**2+il**2
502                jsum = jh**2+jk**2+jl**2
503                if isum == jsum:
504                    eqvDict[varyI].append(varyJ)
505            elif abs(dspI-dspJ)/dspI < 1.e-4:
506                eqvDict[varyI].append(varyJ) 
507    for item in pawleyVary:
508        if eqvDict[item]:
509            for item2 in pawleyVary:
510                if item2 in eqvDict[item]:
511                    eqvDict[item2] = []
512            G2mv.StoreEquivalence(item,eqvDict[item])
513                   
514def cellVary(pfx,SGData): 
515    if SGData['SGLaue'] in ['-1',]:
516        return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3',pfx+'A4',pfx+'A5']
517    elif SGData['SGLaue'] in ['2/m',]:
518        if SGData['SGUniq'] == 'a':
519            return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3']
520        elif SGData['SGUniq'] == 'b':
521            return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A4']
522        else:
523            return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A5']
524    elif SGData['SGLaue'] in ['mmm',]:
525        return [pfx+'A0',pfx+'A1',pfx+'A2']
526    elif SGData['SGLaue'] in ['4/m','4/mmm']:
527        return [pfx+'A0',pfx+'A2']
528    elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
529        return [pfx+'A0',pfx+'A2']
530    elif SGData['SGLaue'] in ['3R', '3mR']:
531        return [pfx+'A0',pfx+'A3']                       
532    elif SGData['SGLaue'] in ['m3m','m3']:
533        return [pfx+'A0',]
534       
535################################################################################
536##### Rigid Body Models
537################################################################################
538       
539def GetRigidBodyModels(rigidbodyDict,Print=True,pFile=None):
540   
541    def PrintResRBModel(RBModel):
542        atNames = RBModel['atNames']
543        rbRef = RBModel['rbRef']
544        rbSeq = RBModel['rbSeq']
545        print >>pFile,'Residue RB name: ',RBModel['RBname'],' no.atoms: ',len(RBModel['rbTypes']), \
546            'No. times used: ',RBModel['useCount']
547        print >>pFile,'    At name       x          y          z'
548        for name,xyz in zip(atNames,RBModel['rbXYZ']):
549            print >>pFile,%8s %10.4f %10.4f %10.4f'%(name,xyz[0],xyz[1],xyz[2])
550        print >>pFile,'Orientation defined by:',atNames[rbRef[0]],' -> ',atNames[rbRef[1]], \
551            ' & ',atNames[rbRef[0]],' -> ',atNames[rbRef[2]]
552        if rbSeq:
553            for i,rbseq in enumerate(rbSeq):
554                print >>pFile,'Torsion sequence ',i,' Bond: '+atNames[rbseq[0]],' - ', \
555                    atNames[rbseq[1]],' riding: ',[atNames[i] for i in rbseq[3]]
556       
557    def PrintVecRBModel(RBModel):
558        rbRef = RBModel['rbRef']
559        atTypes = RBModel['rbTypes']
560        print >>pFile,'Vector RB name: ',RBModel['RBname'],' no.atoms: ',len(RBModel['rbTypes']), \
561            'No. times used: ',RBModel['useCount']
562        for i in range(len(RBModel['VectMag'])):
563            print >>pFile,'Vector no.: ',i,' Magnitude: ', \
564                '%8.4f'%(RBModel['VectMag'][i]),' Refine? ',RBModel['VectRef'][i]
565            print >>pFile,'  No. Type     vx         vy         vz'
566            for j,[name,xyz] in enumerate(zip(atTypes,RBModel['rbVect'][i])):
567                print >>pFile,%d   %2s %10.4f %10.4f %10.4f'%(j,name,xyz[0],xyz[1],xyz[2])
568        print >>pFile,'  No. Type      x          y          z'
569        for i,[name,xyz] in enumerate(zip(atTypes,RBModel['rbXYZ'])):
570            print >>pFile,%d   %2s %10.4f %10.4f %10.4f'%(i,name,xyz[0],xyz[1],xyz[2])
571        print >>pFile,'Orientation defined by: atom ',rbRef[0],' -> atom ',rbRef[1], \
572            ' & atom ',rbRef[0],' -> atom ',rbRef[2]
573    rbVary = []
574    rbDict = {}
575    rbIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]})
576    if len(rbIds['Vector']):
577        for irb,item in enumerate(rbIds['Vector']):
578            if rigidbodyDict['Vector'][item]['useCount']:
579                RBmags = rigidbodyDict['Vector'][item]['VectMag']
580                RBrefs = rigidbodyDict['Vector'][item]['VectRef']
581                for i,[mag,ref] in enumerate(zip(RBmags,RBrefs)):
582                    pid = '::RBV;'+str(i)+':'+str(irb)
583                    rbDict[pid] = mag
584                    if ref:
585                        rbVary.append(pid)
586                if Print:
587                    print >>pFile,'\nVector rigid body model:'
588                    PrintVecRBModel(rigidbodyDict['Vector'][item])
589    if len(rbIds['Residue']):
590        for item in rbIds['Residue']:
591            if rigidbodyDict['Residue'][item]['useCount']:
592                if Print:
593                    print >>pFile,'\nResidue rigid body model:'
594                    PrintResRBModel(rigidbodyDict['Residue'][item])
595    return rbVary,rbDict
596   
597def SetRigidBodyModels(parmDict,sigDict,rigidbodyDict,pFile=None):
598   
599    def PrintRBVectandSig(VectRB,VectSig):
600        print >>pFile,'\n Rigid body vector magnitudes for '+VectRB['RBname']+':'
601        namstr = '  names :'
602        valstr = '  values:'
603        sigstr = '  esds  :'
604        for i,[val,sig] in enumerate(zip(VectRB['VectMag'],VectSig)):
605            namstr += '%12s'%('Vect '+str(i))
606            valstr += '%12.4f'%(val)
607            if sig:
608                sigstr += '%12.4f'%(sig)
609            else:
610                sigstr += 12*' '
611        print >>pFile,namstr
612        print >>pFile,valstr
613        print >>pFile,sigstr       
614       
615    RBIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]})  #these are lists of rbIds
616    if not RBIds['Vector']:
617        return
618    for irb,item in enumerate(RBIds['Vector']):
619        if rigidbodyDict['Vector'][item]['useCount']:
620            VectSig = []
621            RBmags = rigidbodyDict['Vector'][item]['VectMag']
622            for i,mag in enumerate(RBmags):
623                name = '::RBV;'+str(i)+':'+str(irb)
624                mag = parmDict[name]
625                if name in sigDict:
626                    VectSig.append(sigDict[name])
627            PrintRBVectandSig(rigidbodyDict['Vector'][item],VectSig)   
628       
629def ApplyRBModels(parmDict,Phases,rigidbodyDict,Update=False):
630    ''' Takes RB info from RBModels in Phase and RB data in rigidbodyDict along with
631    current RB values in parmDict & modifies atom contents (xyz & Uij) of parmDict
632    '''
633    atxIds = ['Ax:','Ay:','Az:']
634    atuIds = ['AU11:','AU22:','AU33:','AU12:','AU13:','AU23:']
635    RBIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]})  #these are lists of rbIds
636    if not RBIds['Vector'] and not RBIds['Residue']:
637        return
638    VRBIds = RBIds['Vector']
639    RRBIds = RBIds['Residue']
640    if Update:
641        RBData = rigidbodyDict
642    else:
643        RBData = copy.deepcopy(rigidbodyDict)     # don't mess with original!
644    if RBIds['Vector']:                       # first update the vector magnitudes
645        VRBData = RBData['Vector']
646        for i,rbId in enumerate(VRBIds):
647            for j in range(len(VRBData[rbId]['VectMag'])):
648                name = '::RBV;'+str(j)+':'+str(i)
649                VRBData[rbId]['VectMag'][j] = parmDict[name]
650    for phase in Phases:
651        Phase = Phases[phase]
652        General = Phase['General']
653        cell = General['Cell'][1:7]
654        Amat,Bmat = G2lat.cell2AB(cell)
655        AtLookup = G2mth.FillAtomLookUp(Phase['Atoms'])
656        pfx = str(Phase['pId'])+'::'
657        if Update:
658            RBModels = Phase['RBModels']
659        else:
660            RBModels =  copy.deepcopy(Phase['RBModels']) # again don't mess with original!
661        for irb,RBObj in enumerate(RBModels.get('Vector',[])):
662            jrb = VRBIds.index(RBObj['RBId'])
663            rbsx = str(irb)+':'+str(jrb)
664            for i,px in enumerate(['RBVPx:','RBVPy:','RBVPz:']):
665                RBObj['Orig'][0][i] = parmDict[pfx+px+rbsx]
666            for i,po in enumerate(['RBVOa:','RBVOi:','RBVOj:','RBVOk:']):
667                RBObj['Orient'][0][i] = parmDict[pfx+po+rbsx]
668            TLS = RBObj['ThermalMotion']
669            if 'T' in TLS[0]:
670                for i,pt in enumerate(['RBVT11:','RBVT22:','RBVT33:','RBVT12:','RBVT13:','RBVT23:']):
671                    TLS[1][i] = parmDict[pfx+pt+rbsx]
672            if 'L' in TLS[0]:
673                for i,pt in enumerate(['RBVL11:','RBVL22:','RBVL33:','RBVL12:','RBVL13:','RBVL23:']):
674                    TLS[1][i+6] = parmDict[pfx+pt+rbsx]
675            if 'S' in TLS[0]:
676                for i,pt in enumerate(['RBVS12:','RBVS13:','RBVS21:','RBVS23:','RBVS31:','RBVS32:','RBVSAA:','RBVSBB:']):
677                    TLS[1][i+12] = parmDict[pfx+pt+rbsx]
678            if 'U' in TLS[0]:
679                TLS[1][0] = parmDict[pfx+'RBVU:'+rbsx]
680            XYZ,Cart = G2mth.UpdateRBXYZ(Bmat,RBObj,RBData,'Vector')
681            UIJ = G2mth.UpdateRBUIJ(Bmat,Cart,RBObj)
682            for i,x in enumerate(XYZ):
683                atId = RBObj['Ids'][i]
684                for j in [0,1,2]:
685                    parmDict[pfx+atxIds[j]+str(AtLookup[atId])] = x[j]
686                if UIJ[i][0] == 'A':
687                    for j in range(6):
688                        parmDict[pfx+atuIds[j]+str(AtLookup[atId])] = UIJ[i][j+2]
689                elif UIJ[i][0] == 'I':
690                    parmDict[pfx+'AUiso:'+str(AtLookup[atId])] = UIJ[i][1]
691           
692        for irb,RBObj in enumerate(RBModels.get('Residue',[])):
693            jrb = RRBIds.index(RBObj['RBId'])
694            rbsx = str(irb)+':'+str(jrb)
695            for i,px in enumerate(['RBRPx:','RBRPy:','RBRPz:']):
696                RBObj['Orig'][0][i] = parmDict[pfx+px+rbsx]
697            for i,po in enumerate(['RBROa:','RBROi:','RBROj:','RBROk:']):
698                RBObj['Orient'][0][i] = parmDict[pfx+po+rbsx]               
699            TLS = RBObj['ThermalMotion']
700            if 'T' in TLS[0]:
701                for i,pt in enumerate(['RBRT11:','RBRT22:','RBRT33:','RBRT12:','RBRT13:','RBRT23:']):
702                    RBObj['ThermalMotion'][1][i] = parmDict[pfx+pt+rbsx]
703            if 'L' in TLS[0]:
704                for i,pt in enumerate(['RBRL11:','RBRL22:','RBRL33:','RBRL12:','RBRL13:','RBRL23:']):
705                    RBObj['ThermalMotion'][1][i+6] = parmDict[pfx+pt+rbsx]
706            if 'S' in TLS[0]:
707                for i,pt in enumerate(['RBRS12:','RBRS13:','RBRS21:','RBRS23:','RBRS31:','RBRS32:','RBRSAA:','RBRSBB:']):
708                    RBObj['ThermalMotion'][1][i+12] = parmDict[pfx+pt+rbsx]
709            if 'U' in TLS[0]:
710                RBObj['ThermalMotion'][1][0] = parmDict[pfx+'RBRU:'+rbsx]
711            for itors,tors in enumerate(RBObj['Torsions']):
712                tors[0] = parmDict[pfx+'RBRTr;'+str(itors)+':'+rbsx]
713            XYZ,Cart = G2mth.UpdateRBXYZ(Bmat,RBObj,RBData,'Residue')
714            UIJ = G2mth.UpdateRBUIJ(Bmat,Cart,RBObj)
715            for i,x in enumerate(XYZ):
716                atId = RBObj['Ids'][i]
717                for j in [0,1,2]:
718                    parmDict[pfx+atxIds[j]+str(AtLookup[atId])] = x[j]
719                if UIJ[i][0] == 'A':
720                    for j in range(6):
721                        parmDict[pfx+atuIds[j]+str(AtLookup[atId])] = UIJ[i][j+2]
722                elif UIJ[i][0] == 'I':
723                    parmDict[pfx+'AUiso:'+str(AtLookup[atId])] = UIJ[i][1]
724                   
725def ApplyRBModelDervs(dFdvDict,parmDict,rigidbodyDict,Phase):
726    atxIds = ['dAx:','dAy:','dAz:']
727    atuIds = ['AU11:','AU22:','AU33:','AU12:','AU13:','AU23:']
728    TIds = ['T11:','T22:','T33:','T12:','T13:','T23:']
729    LIds = ['L11:','L22:','L33:','L12:','L13:','L23:']
730    SIds = ['S12:','S13:','S21:','S23:','S31:','S32:','SAA:','SBB:']
731    PIds = ['Px:','Py:','Pz:']
732    OIds = ['Oa:','Oi:','Oj:','Ok:']
733    RBIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]})  #these are lists of rbIds
734    if not RBIds['Vector'] and not RBIds['Residue']:
735        return
736    VRBIds = RBIds['Vector']
737    RRBIds = RBIds['Residue']
738    RBData = rigidbodyDict
739    for item in parmDict:
740        if 'RB' in item:
741            dFdvDict[item] = 0.        #NB: this is a vector which is no. refl. long & must be filled!
742    General = Phase['General']
743    cell = General['Cell'][1:7]
744    Amat,Bmat = G2lat.cell2AB(cell)
745    rpd = np.pi/180.
746    rpd2 = rpd**2
747    g = np.inner(Bmat,Bmat.T)
748    gvec = np.sqrt(np.array([g[0][0]**2,g[1][1]**2,g[2][2]**2,
749        g[0][0]*g[1][1],g[0][0]*g[2][2],g[1][1]*g[2][2]]))
750    AtLookup = G2mth.FillAtomLookUp(Phase['Atoms'])
751    pfx = str(Phase['pId'])+'::'
752    RBModels =  Phase['RBModels']
753    dx = 0.0001
754    for irb,RBObj in enumerate(RBModels.get('Vector',[])):
755        VModel = RBData['Vector'][RBObj['RBId']]
756        Q = RBObj['Orient'][0]
757        Pos = RBObj['Orig'][0]
758        jrb = VRBIds.index(RBObj['RBId'])
759        rbsx = str(irb)+':'+str(jrb)
760        dXdv = []
761        for iv in range(len(VModel['VectMag'])):
762            dXdv.append(np.inner(Bmat.T,VModel['rbVect'][iv]).T)
763        XYZ,Cart = G2mth.UpdateRBXYZ(Bmat,RBObj,RBData,'Vector')
764        for ia,atId in enumerate(RBObj['Ids']):
765            atNum = AtLookup[atId]
766            for iv in range(len(VModel['VectMag'])):
767                for ix in [0,1,2]:
768                    dFdvDict['::RBV;'+str(iv)+':'+str(jrb)] += dXdv[iv][ia][ix]*dFdvDict[pfx+atxIds[ix]+str(atNum)]
769            for i,name in enumerate(['RBVPx:','RBVPy:','RBVPz:']):
770                dFdvDict[pfx+name+rbsx] += dFdvDict[pfx+atxIds[i]+str(atNum)]
771            for iv in range(4):
772                Q[iv] -= dx
773                XYZ1,Cart1 = G2mth.UpdateRBXYZ(Bmat,RBObj,RBData,'Vector')
774                Q[iv] += 2.*dx
775                XYZ2,Cart2 = G2mth.UpdateRBXYZ(Bmat,RBObj,RBData,'Vector')
776                Q[iv] -= dx
777                dXdO = (XYZ2[ia]-XYZ1[ia])/(2.*dx)
778                for ix in [0,1,2]:
779                    dFdvDict[pfx+'RBV'+OIds[iv]+rbsx] += dXdO[ix]*dFdvDict[pfx+atxIds[ix]+str(atNum)]
780            X = G2mth.prodQVQ(Q,Cart[ia])
781            dFdu = np.array([dFdvDict[pfx+Uid+str(AtLookup[atId])] for Uid in atuIds]).T*gvec
782            dFdu = G2lat.U6toUij(dFdu.T)
783            dFdu = np.tensordot(Amat.T,np.tensordot(Amat,dFdu,([1,0])),([0,1]))
784            dFdu = G2lat.UijtoU6(dFdu)
785            atNum = AtLookup[atId]
786            if 'T' in RBObj['ThermalMotion'][0]:
787                for i,name in enumerate(['RBVT11:','RBVT22:','RBVT33:','RBVT12:','RBVT13:','RBVT23:']):
788                    dFdvDict[pfx+name+rbsx] += dFdu[i]
789            if 'L' in RBObj['ThermalMotion'][0]:
790                dFdvDict[pfx+'RBVL11:'+rbsx] += rpd2*(dFdu[1]*X[2]**2+dFdu[2]*X[1]**2-dFdu[5]*X[1]*X[2])
791                dFdvDict[pfx+'RBVL22:'+rbsx] += rpd2*(dFdu[0]*X[2]**2+dFdu[2]*X[0]**2-dFdu[4]*X[0]*X[2])
792                dFdvDict[pfx+'RBVL33:'+rbsx] += rpd2*(dFdu[0]*X[1]**2+dFdu[1]*X[0]**2-dFdu[3]*X[0]*X[1])
793                dFdvDict[pfx+'RBVL12:'+rbsx] += rpd2*(-dFdu[3]*X[2]**2-2.*dFdu[2]*X[0]*X[1]+
794                    dFdu[4]*X[1]*X[2]+dFdu[5]*X[0]*X[2])
795                dFdvDict[pfx+'RBVL13:'+rbsx] += rpd2*(dFdu[0]*X[1]**2-2.*dFdu[1]*X[0]*X[2]+
796                    dFdu[3]*X[1]*X[2]+dFdu[5]*X[0]*X[1])
797                dFdvDict[pfx+'RBVL23:'+rbsx] += rpd2*(dFdu[0]*X[1]**2-2.*dFdu[0]*X[1]*X[2]+
798                    dFdu[3]*X[0]*X[2]+dFdu[4]*X[0]*X[1])
799            if 'S' in RBObj['ThermalMotion'][0]:
800                dFdvDict[pfx+'RBVS12:'+rbsx] += rpd*(dFdu[5]*X[1]-2.*dFdu[1]*X[2])
801                dFdvDict[pfx+'RBVS13:'+rbsx] += rpd*(-dFdu[5]*X[2]+2.*dFdu[2]*X[1])
802                dFdvDict[pfx+'RBVS21:'+rbsx] += rpd*(-dFdu[4]*X[0]+2.*dFdu[0]*X[2])
803                dFdvDict[pfx+'RBVS23:'+rbsx] += rpd*(dFdu[4]*X[2]-2.*dFdu[2]*X[0])
804                dFdvDict[pfx+'RBVS31:'+rbsx] += rpd*(dFdu[3]*X[0]-2.*dFdu[0]*X[1])
805                dFdvDict[pfx+'RBVS32:'+rbsx] += rpd*(-dFdu[3]*X[1]+2.*dFdu[1]*X[0])
806                dFdvDict[pfx+'RBVSAA:'+rbsx] += rpd*(dFdu[4]*X[1]-dFdu[3]*X[2])
807                dFdvDict[pfx+'RBVSBB:'+rbsx] += rpd*(dFdu[5]*X[0]-dFdu[3]*X[2])
808            if 'U' in RBObj['ThermalMotion'][0]:
809                dFdvDict[pfx+'RBVU:'+rbsx] += dFdvDict[pfx+'AUiso:'+str(AtLookup[atId])]
810
811
812    for irb,RBObj in enumerate(RBModels.get('Residue',[])):
813        Q = RBObj['Orient'][0]
814        Pos = RBObj['Orig'][0]
815        jrb = RRBIds.index(RBObj['RBId'])
816        torData = RBData['Residue'][RBObj['RBId']]['rbSeq']
817        rbsx = str(irb)+':'+str(jrb)
818        XYZ,Cart = G2mth.UpdateRBXYZ(Bmat,RBObj,RBData,'Residue')
819        for itors,tors in enumerate(RBObj['Torsions']):
820            tname = pfx+'RBRTr;'+str(itors)+':'+rbsx           
821            orId,pvId = torData[itors][:2]
822            pivotVec = Cart[pvId]-Cart[orId]
823            QA = G2mth.AVdeg2Q(0.001,pivotVec)
824            QB = G2mth.AVdeg2Q(-0.001,pivotVec)
825            for ir in torData[itors][3]:
826                atNum = AtLookup[RBObj['Ids'][ir]]
827                rVec = Cart[ir]-Cart[pvId]
828                dRdT = np.inner(Bmat,(G2mth.prodQVQ(QA,rVec)-G2mth.prodQVQ(QB,rVec)))/.002
829                for ix in [0,1,2]:
830                    dFdvDict[tname] += dRdT[ix]*dFdvDict[pfx+atxIds[ix]+str(atNum)]
831        for ia,atId in enumerate(RBObj['Ids']):
832            atNum = AtLookup[atId]
833            dx = 0.0001
834            for i,name in enumerate(['RBRPx:','RBRPy:','RBRPz:']):
835                dFdvDict[pfx+name+rbsx] += dFdvDict[pfx+atxIds[i]+str(atNum)]
836            for iv in range(4):
837                Q[iv] -= dx
838                XYZ1,Cart1 = G2mth.UpdateRBXYZ(Bmat,RBObj,RBData,'Residue')
839                Q[iv] += 2.*dx
840                XYZ2,Cart2 = G2mth.UpdateRBXYZ(Bmat,RBObj,RBData,'Residue')
841                Q[iv] -= dx
842                dXdO = (XYZ2[ia]-XYZ1[ia])/(2.*dx)
843                for ix in [0,1,2]:
844                    dFdvDict[pfx+'RBR'+OIds[iv]+rbsx] += dXdO[ix]*dFdvDict[pfx+atxIds[ix]+str(atNum)]
845            X = G2mth.prodQVQ(Q,Cart[ia])
846            dFdu = np.array([dFdvDict[pfx+Uid+str(AtLookup[atId])] for Uid in atuIds]).T*gvec
847            dFdu = G2lat.U6toUij(dFdu.T)
848            dFdu = np.tensordot(Amat.T,np.tensordot(Amat,dFdu,([1,0])),([0,1]))
849            dFdu = G2lat.UijtoU6(dFdu)
850            atNum = AtLookup[atId]
851            if 'T' in RBObj['ThermalMotion'][0]:
852                for i,name in enumerate(['RBRT11:','RBRT22:','RBRT33:','RBRT12:','RBRT13:','RBRT23:']):
853                    dFdvDict[pfx+name+rbsx] += dFdu[i]
854            if 'L' in RBObj['ThermalMotion'][0]:
855                dFdvDict[pfx+'RBRL11:'+rbsx] += rpd2*(dFdu[1]*X[2]**2+dFdu[2]*X[1]**2-dFdu[5]*X[1]*X[2])
856                dFdvDict[pfx+'RBRL22:'+rbsx] += rpd2*(dFdu[0]*X[2]**2+dFdu[2]*X[0]**2-dFdu[4]*X[0]*X[2])
857                dFdvDict[pfx+'RBRL33:'+rbsx] += rpd2*(dFdu[0]*X[1]**2+dFdu[1]*X[0]**2-dFdu[3]*X[0]*X[1])
858                dFdvDict[pfx+'RBRL12:'+rbsx] += rpd2*(-dFdu[3]*X[2]**2-2.*dFdu[2]*X[0]*X[1]+
859                    dFdu[4]*X[1]*X[2]+dFdu[5]*X[0]*X[2])
860                dFdvDict[pfx+'RBRL13:'+rbsx] += rpd2*(dFdu[0]*X[1]**2-2.*dFdu[1]*X[0]*X[2]+
861                    dFdu[3]*X[1]*X[2]+dFdu[5]*X[0]*X[1])
862                dFdvDict[pfx+'RBRL23:'+rbsx] += rpd2*(dFdu[0]*X[1]**2-2.*dFdu[0]*X[1]*X[2]+
863                    dFdu[3]*X[0]*X[2]+dFdu[4]*X[0]*X[1])
864            if 'S' in RBObj['ThermalMotion'][0]:
865                dFdvDict[pfx+'RBRS12:'+rbsx] += rpd*(dFdu[5]*X[1]-2.*dFdu[1]*X[2])
866                dFdvDict[pfx+'RBRS13:'+rbsx] += rpd*(-dFdu[5]*X[2]+2.*dFdu[2]*X[1])
867                dFdvDict[pfx+'RBRS21:'+rbsx] += rpd*(-dFdu[4]*X[0]+2.*dFdu[0]*X[2])
868                dFdvDict[pfx+'RBRS23:'+rbsx] += rpd*(dFdu[4]*X[2]-2.*dFdu[2]*X[0])
869                dFdvDict[pfx+'RBRS31:'+rbsx] += rpd*(dFdu[3]*X[0]-2.*dFdu[0]*X[1])
870                dFdvDict[pfx+'RBRS32:'+rbsx] += rpd*(-dFdu[3]*X[1]+2.*dFdu[1]*X[0])
871                dFdvDict[pfx+'RBRSAA:'+rbsx] += rpd*(dFdu[4]*X[1]-dFdu[3]*X[2])
872                dFdvDict[pfx+'RBRSBB:'+rbsx] += rpd*(dFdu[5]*X[0]-dFdu[3]*X[2])
873            if 'U' in RBObj['ThermalMotion'][0]:
874                dFdvDict[pfx+'RBRU:'+rbsx] += dFdvDict[pfx+'AUiso:'+str(AtLookup[atId])]
875   
876       
877################################################################################
878##### Phase data
879################################################################################       
880                   
881def GetPhaseData(PhaseData,RestraintDict={},rbIds={},Print=True,pFile=None):
882           
883    def PrintFFtable(FFtable):
884        print >>pFile,'\n X-ray scattering factors:'
885        print >>pFile,'   Symbol     fa                                      fb                                      fc'
886        print >>pFile,99*'-'
887        for Ename in FFtable:
888            ffdata = FFtable[Ename]
889            fa = ffdata['fa']
890            fb = ffdata['fb']
891            print >>pFile,' %8s %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f' %  \
892                (Ename.ljust(8),fa[0],fa[1],fa[2],fa[3],fb[0],fb[1],fb[2],fb[3],ffdata['fc'])
893               
894    def PrintBLtable(BLtable):
895        print >>pFile,'\n Neutron scattering factors:'
896        print >>pFile,'   Symbol   isotope       mass       b       resonant terms'
897        print >>pFile,99*'-'
898        for Ename in BLtable:
899            bldata = BLtable[Ename]
900            isotope = bldata[0]
901            mass = bldata[1][0]
902            blen = bldata[1][1]
903            bres = []
904            if len(bldata[1]) > 2:
905                bres = bldata[1][2:]
906            line = ' %8s%11s %10.3f %8.3f'%(Ename.ljust(8),isotope.center(11),mass,blen)
907            for item in bres:
908                line += '%10.5g'%(item)
909            print >>pFile,line
910           
911    def PrintRBObjects(resRBData,vecRBData):
912       
913        def PrintRBThermals():
914            tlstr = ['11','22','33','12','13','23']
915            sstr = ['12','13','21','23','31','32','AA','BB']
916            TLS = RB['ThermalMotion'][1]
917            TLSvar = RB['ThermalMotion'][2]
918            if 'T' in RB['ThermalMotion'][0]:
919                print >>pFile,'TLS data'
920                text = ''
921                for i in range(6):
922                    text += 'T'+tlstr[i]+' %8.4f %s '%(TLS[i],str(TLSvar[i])[0])
923                print >>pFile,text
924                if 'L' in RB['ThermalMotion'][0]: 
925                    text = ''
926                    for i in range(6,12):
927                        text += 'L'+tlstr[i-6]+' %8.2f %s '%(TLS[i],str(TLSvar[i])[0])
928                    print >>pFile,text
929                if 'S' in RB['ThermalMotion'][0]:
930                    text = ''
931                    for i in range(12,20):
932                        text += 'S'+sstr[i-12]+' %8.3f %s '%(TLS[i],str(TLSvar[i])[0])
933                    print >>pFile,text
934            if 'U' in RB['ThermalMotion'][0]:
935                print >>pFile,'Uiso data'
936                text = 'Uiso'+' %10.3f %s'%(TLS[0],str(TLSvar[0])[0])           
937           
938        if len(resRBData):
939            for RB in resRBData:
940                Oxyz = RB['Orig'][0]
941                Qrijk = RB['Orient'][0]
942                Angle = 2.0*acosd(Qrijk[0])
943                print >>pFile,'\nRBObject ',RB['RBname'],' at ',      \
944                    '%10.4f %10.4f %10.4f'%(Oxyz[0],Oxyz[1],Oxyz[2]),' Refine?',RB['Orig'][1]
945                print >>pFile,'Orientation angle,vector:',      \
946                    '%10.3f %10.4f %10.4f %10.4f'%(Angle,Qrijk[1],Qrijk[2],Qrijk[3]),' Refine? ',RB['Orient'][1]
947                Torsions = RB['Torsions']
948                if len(Torsions):
949                    text = 'Torsions: '
950                    for torsion in Torsions:
951                        text += '%10.4f Refine? %s'%(torsion[0],torsion[1])
952                    print >>pFile,text
953                PrintRBThermals()
954        if len(vecRBData):
955            for RB in vecRBData:
956                Oxyz = RB['Orig'][0]
957                Qrijk = RB['Orient'][0]
958                Angle = 2.0*acosd(Qrijk[0])
959                print >>pFile,'\nRBObject ',RB['RBname'],' at ',      \
960                    '%10.4f %10.4f %10.4f'%(Oxyz[0],Oxyz[1],Oxyz[2]),' Refine?',RB['Orig'][1]           
961                print >>pFile,'Orientation angle,vector:',      \
962                    '%10.3f %10.4f %10.4f %10.4f'%(Angle,Qrijk[1],Qrijk[2],Qrijk[3]),' Refine? ',RB['Orient'][1]
963                PrintRBThermals()
964               
965    def PrintAtoms(General,Atoms):
966        cx,ct,cs,cia = General['AtomPtrs']
967        print >>pFile,'\n Atoms:'
968        line = '   name    type  refine?   x         y         z    '+ \
969            '  frac site sym  mult I/A   Uiso     U11     U22     U33     U12     U13     U23'
970        if General['Type'] == 'magnetic':
971            line += '   Mx     My     Mz'
972        elif General['Type'] == 'macromolecular':
973            line = ' res no residue chain'+line
974        print >>pFile,line
975        if General['Type'] == 'nuclear':
976            print >>pFile,135*'-'
977            for i,at in enumerate(Atoms):
978                line = '%7s'%(at[ct-1])+'%7s'%(at[ct])+'%7s'%(at[ct+1])+'%10.5f'%(at[cx])+'%10.5f'%(at[cx+1])+ \
979                    '%10.5f'%(at[cx+2])+'%8.3f'%(at[cx+3])+'%7s'%(at[cs])+'%5d'%(at[cs+1])+'%5s'%(at[cia])
980                if at[cia] == 'I':
981                    line += '%8.4f'%(at[cia+1])+48*' '
982                else:
983                    line += 8*' '
984                    for j in range(6):
985                        line += '%8.4f'%(at[cia+1+j])
986                print >>pFile,line
987        elif General['Type'] == 'macromolecular':
988            print >>pFile,135*'-'           
989            for i,at in enumerate(Atoms):
990                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])+ \
991                    '%10.5f'%(at[cx+2])+'%8.3f'%(at[cx+3])+'%7s'%(at[cs])+'%5d'%(at[cs+1])+'%5s'%(at[cia])
992                if at[cia] == 'I':
993                    line += '%8.4f'%(at[cia+1])+48*' '
994                else:
995                    line += 8*' '
996                    for j in range(6):
997                        line += '%8.4f'%(at[cia+1+j])
998                print >>pFile,line
999       
1000    def PrintTexture(textureData):
1001        topstr = '\n Spherical harmonics texture: Order:' + \
1002            str(textureData['Order'])
1003        if textureData['Order']:
1004            print >>pFile,topstr+' Refine? '+str(textureData['SH Coeff'][0])
1005        else:
1006            print >>pFile,topstr
1007            return
1008        names = ['omega','chi','phi']
1009        line = '\n'
1010        for name in names:
1011            line += ' SH '+name+':'+'%12.4f'%(textureData['Sample '+name][1])+' Refine? '+str(textureData['Sample '+name][0])
1012        print >>pFile,line
1013        print >>pFile,'\n Texture coefficients:'
1014        ptlbls = ' names :'
1015        ptstr =  ' values:'
1016        SHcoeff = textureData['SH Coeff'][1]
1017        for item in SHcoeff:
1018            ptlbls += '%12s'%(item)
1019            ptstr += '%12.4f'%(SHcoeff[item]) 
1020        print >>pFile,ptlbls
1021        print >>pFile,ptstr
1022       
1023    def MakeRBParms(rbKey,phaseVary,phaseDict):
1024        rbid = str(rbids.index(RB['RBId']))
1025        pfxRB = pfx+'RB'+rbKey+'P'
1026        pstr = ['x','y','z']
1027        ostr = ['a','i','j','k']
1028        for i in range(3):
1029            name = pfxRB+pstr[i]+':'+str(iRB)+':'+rbid
1030            phaseDict[name] = RB['Orig'][0][i]
1031            if RB['Orig'][1]:
1032                phaseVary += [name,]
1033        pfxRB = pfx+'RB'+rbKey+'O'
1034        for i in range(4):
1035            name = pfxRB+ostr[i]+':'+str(iRB)+':'+rbid
1036            phaseDict[name] = RB['Orient'][0][i]
1037            if RB['Orient'][1] == 'AV' and i:
1038                phaseVary += [name,]
1039            elif RB['Orient'][1] == 'A' and not i:
1040                phaseVary += [name,]
1041           
1042    def MakeRBThermals(rbKey,phaseVary,phaseDict):
1043        rbid = str(rbids.index(RB['RBId']))
1044        tlstr = ['11','22','33','12','13','23']
1045        sstr = ['12','13','21','23','31','32','AA','BB']
1046        if 'T' in RB['ThermalMotion'][0]:
1047            pfxRB = pfx+'RB'+rbKey+'T'
1048            for i in range(6):
1049                name = pfxRB+tlstr[i]+':'+str(iRB)+':'+rbid
1050                phaseDict[name] = RB['ThermalMotion'][1][i]
1051                if RB['ThermalMotion'][2][i]:
1052                    phaseVary += [name,]
1053        if 'L' in RB['ThermalMotion'][0]:
1054            pfxRB = pfx+'RB'+rbKey+'L'
1055            for i in range(6):
1056                name = pfxRB+tlstr[i]+':'+str(iRB)+':'+rbid
1057                phaseDict[name] = RB['ThermalMotion'][1][i+6]
1058                if RB['ThermalMotion'][2][i+6]:
1059                    phaseVary += [name,]
1060        if 'S' in RB['ThermalMotion'][0]:
1061            pfxRB = pfx+'RB'+rbKey+'S'
1062            for i in range(8):
1063                name = pfxRB+sstr[i]+':'+str(iRB)+':'+rbid
1064                phaseDict[name] = RB['ThermalMotion'][1][i+12]
1065                if RB['ThermalMotion'][2][i+12]:
1066                    phaseVary += [name,]
1067        if 'U' in RB['ThermalMotion'][0]:
1068            name = pfx+'RB'+rbKey+'U:'+str(iRB)+':'+rbid
1069            phaseDict[name] = RB['ThermalMotion'][1][0]
1070            if RB['ThermalMotion'][2][0]:
1071                phaseVary += [name,]
1072               
1073    def MakeRBTorsions(rbKey,phaseVary,phaseDict):
1074        rbid = str(rbids.index(RB['RBId']))
1075        pfxRB = pfx+'RB'+rbKey+'Tr;'
1076        for i,torsion in enumerate(RB['Torsions']):
1077            name = pfxRB+str(i)+':'+str(iRB)+':'+rbid
1078            phaseDict[name] = torsion[0]
1079            if torsion[1]:
1080                phaseVary += [name,]
1081                   
1082    if Print:
1083        print  >>pFile,'\n Phases:'
1084    phaseVary = []
1085    phaseDict = {}
1086    phaseConstr = {}
1087    pawleyLookup = {}
1088    FFtables = {}                   #scattering factors - xrays
1089    BLtables = {}                   # neutrons
1090    Natoms = {}
1091    AtMults = {}
1092    AtIA = {}
1093    shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
1094    SamSym = dict(zip(shModels,['0','-1','2/m','mmm']))
1095    atomIndx = {}
1096    for name in PhaseData:
1097        General = PhaseData[name]['General']
1098        pId = PhaseData[name]['pId']
1099        pfx = str(pId)+'::'
1100        FFtable = GetFFtable(General)
1101        BLtable = GetBLtable(General)
1102        FFtables.update(FFtable)
1103        BLtables.update(BLtable)
1104        Atoms = PhaseData[name]['Atoms']
1105        AtLookup = G2mth.FillAtomLookUp(Atoms)
1106        PawleyRef = PhaseData[name].get('Pawley ref',[])
1107        SGData = General['SGData']
1108        SGtext = G2spc.SGPrint(SGData)
1109        cell = General['Cell']
1110        A = G2lat.cell2A(cell[1:7])
1111        phaseDict.update({pfx+'A0':A[0],pfx+'A1':A[1],pfx+'A2':A[2],
1112            pfx+'A3':A[3],pfx+'A4':A[4],pfx+'A5':A[5],pfx+'Vol':G2lat.calc_V(A)})
1113        if cell[0]:
1114            phaseVary += cellVary(pfx,SGData)
1115        resRBData = PhaseData[name]['RBModels'].get('Residue',[])
1116        if resRBData:
1117            rbids = rbIds['Residue']    #NB: used in the MakeRB routines
1118            for iRB,RB in enumerate(resRBData):
1119                MakeRBParms('R',phaseVary,phaseDict)
1120                MakeRBThermals('R',phaseVary,phaseDict)
1121                MakeRBTorsions('R',phaseVary,phaseDict)
1122       
1123        vecRBData = PhaseData[name]['RBModels'].get('Vector',[])
1124        if vecRBData:
1125            rbids = rbIds['Vector']    #NB: used in the MakeRB routines
1126            for iRB,RB in enumerate(vecRBData):
1127                MakeRBParms('V',phaseVary,phaseDict)
1128                MakeRBThermals('V',phaseVary,phaseDict)
1129                   
1130        Natoms[pfx] = 0
1131        if Atoms and not General.get('doPawley'):
1132            cx,ct,cs,cia = General['AtomPtrs']
1133            if General['Type'] in ['nuclear','macromolecular']:
1134                Natoms[pfx] = len(Atoms)
1135                for i,at in enumerate(Atoms):
1136                    atomIndx[at[-1]] = [pfx,i]      #lookup table for restraints
1137                    phaseDict.update({pfx+'Atype:'+str(i):at[ct],pfx+'Afrac:'+str(i):at[cx+3],pfx+'Amul:'+str(i):at[cs+1],
1138                        pfx+'Ax:'+str(i):at[cx],pfx+'Ay:'+str(i):at[cx+1],pfx+'Az:'+str(i):at[cx+2],
1139                        pfx+'dAx:'+str(i):0.,pfx+'dAy:'+str(i):0.,pfx+'dAz:'+str(i):0.,         #refined shifts for x,y,z
1140                        pfx+'AI/A:'+str(i):at[cia],})
1141                    if at[cia] == 'I':
1142                        phaseDict[pfx+'AUiso:'+str(i)] = at[cia+1]
1143                    else:
1144                        phaseDict.update({pfx+'AU11:'+str(i):at[cia+2],pfx+'AU22:'+str(i):at[cia+3],pfx+'AU33:'+str(i):at[cia+4],
1145                            pfx+'AU12:'+str(i):at[cia+5],pfx+'AU13:'+str(i):at[cia+6],pfx+'AU23:'+str(i):at[cia+7]})
1146                    if 'F' in at[ct+1]:
1147                        phaseVary.append(pfx+'Afrac:'+str(i))
1148                    if 'X' in at[ct+1]:
1149                        xId,xCoef = G2spc.GetCSxinel(at[cs])
1150                        names = [pfx+'dAx:'+str(i),pfx+'dAy:'+str(i),pfx+'dAz:'+str(i)]
1151                        equivs = [[],[],[]]
1152                        for j in range(3):
1153                            if xId[j] > 0:                               
1154                                phaseVary.append(names[j])
1155                                equivs[xId[j]-1].append([names[j],xCoef[j]])
1156                        for equiv in equivs:
1157                            if len(equiv) > 1:
1158                                name = equiv[0][0]
1159                                for eqv in equiv[1:]:
1160                                    G2mv.StoreEquivalence(name,(eqv,))
1161                    if 'U' in at[ct+1]:
1162                        if at[9] == 'I':
1163                            phaseVary.append(pfx+'AUiso:'+str(i))
1164                        else:
1165                            uId,uCoef = G2spc.GetCSuinel(at[cs])[:2]
1166                            names = [pfx+'AU11:'+str(i),pfx+'AU22:'+str(i),pfx+'AU33:'+str(i),
1167                                pfx+'AU12:'+str(i),pfx+'AU13:'+str(i),pfx+'AU23:'+str(i)]
1168                            equivs = [[],[],[],[],[],[]]
1169                            for j in range(6):
1170                                if uId[j] > 0:                               
1171                                    phaseVary.append(names[j])
1172                                    equivs[uId[j]-1].append([names[j],uCoef[j]])
1173                            for equiv in equivs:
1174                                if len(equiv) > 1:
1175                                    name = equiv[0][0]
1176                                    for eqv in equiv[1:]:
1177                                        G2mv.StoreEquivalence(name,(eqv,))
1178#            elif General['Type'] == 'magnetic':
1179#            elif General['Type'] == 'macromolecular':
1180            textureData = General['SH Texture']
1181            if textureData['Order']:
1182                phaseDict[pfx+'SHorder'] = textureData['Order']
1183                phaseDict[pfx+'SHmodel'] = SamSym[textureData['Model']]
1184                for item in ['omega','chi','phi']:
1185                    phaseDict[pfx+'SH '+item] = textureData['Sample '+item][1]
1186                    if textureData['Sample '+item][0]:
1187                        phaseVary.append(pfx+'SH '+item)
1188                for item in textureData['SH Coeff'][1]:
1189                    phaseDict[pfx+item] = textureData['SH Coeff'][1][item]
1190                    if textureData['SH Coeff'][0]:
1191                        phaseVary.append(pfx+item)
1192               
1193            if Print:
1194                print >>pFile,'\n Phase name: ',General['Name']
1195                print >>pFile,135*'-'
1196                PrintFFtable(FFtable)
1197                PrintBLtable(BLtable)
1198                print >>pFile,''
1199                for line in SGtext: print >>pFile,line
1200                PrintRBObjects(resRBData,vecRBData)
1201                PrintAtoms(General,Atoms)
1202                print >>pFile,'\n Unit cell: a =','%.5f'%(cell[1]),' b =','%.5f'%(cell[2]),' c =','%.5f'%(cell[3]), \
1203                    ' alpha =','%.3f'%(cell[4]),' beta =','%.3f'%(cell[5]),' gamma =', \
1204                    '%.3f'%(cell[6]),' volume =','%.3f'%(cell[7]),' Refine?',cell[0]
1205                PrintTexture(textureData)
1206                if name in RestraintDict:
1207                    PrintRestraints(cell[1:7],SGData,General['AtomPtrs'],Atoms,AtLookup,
1208                        textureData,RestraintDict[name],pFile)
1209                   
1210        elif PawleyRef:
1211            pawleyVary = []
1212            for i,refl in enumerate(PawleyRef):
1213                phaseDict[pfx+'PWLref:'+str(i)] = refl[6]
1214                pawleyLookup[pfx+'%d,%d,%d'%(refl[0],refl[1],refl[2])] = i
1215                if refl[5]:
1216                    pawleyVary.append(pfx+'PWLref:'+str(i))
1217            GetPawleyConstr(SGData['SGLaue'],PawleyRef,pawleyVary)      #does G2mv.StoreEquivalence
1218            phaseVary += pawleyVary
1219               
1220    return Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables
1221   
1222def cellFill(pfx,SGData,parmDict,sigDict): 
1223    if SGData['SGLaue'] in ['-1',]:
1224        A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
1225            parmDict[pfx+'A3'],parmDict[pfx+'A4'],parmDict[pfx+'A5']]
1226        sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
1227            sigDict[pfx+'A3'],sigDict[pfx+'A4'],sigDict[pfx+'A5']]
1228    elif SGData['SGLaue'] in ['2/m',]:
1229        if SGData['SGUniq'] == 'a':
1230            A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
1231                parmDict[pfx+'A3'],0,0]
1232            sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
1233                sigDict[pfx+'A3'],0,0]
1234        elif SGData['SGUniq'] == 'b':
1235            A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
1236                0,parmDict[pfx+'A4'],0]
1237            sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
1238                0,sigDict[pfx+'A4'],0]
1239        else:
1240            A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
1241                0,0,parmDict[pfx+'A5']]
1242            sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
1243                0,0,sigDict[pfx+'A5']]
1244    elif SGData['SGLaue'] in ['mmm',]:
1245        A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],0,0,0]
1246        sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],0,0,0]
1247    elif SGData['SGLaue'] in ['4/m','4/mmm']:
1248        A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A2'],0,0,0]
1249        sigA = [sigDict[pfx+'A0'],0,sigDict[pfx+'A2'],0,0,0]
1250    elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
1251        A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A2'],
1252            parmDict[pfx+'A0'],0,0]
1253        sigA = [sigDict[pfx+'A0'],0,sigDict[pfx+'A2'],0,0,0]
1254    elif SGData['SGLaue'] in ['3R', '3mR']:
1255        A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A0'],
1256            parmDict[pfx+'A3'],parmDict[pfx+'A3'],parmDict[pfx+'A3']]
1257        sigA = [sigDict[pfx+'A0'],0,0,sigDict[pfx+'A3'],0,0]
1258    elif SGData['SGLaue'] in ['m3m','m3']:
1259        A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A0'],0,0,0]
1260        sigA = [sigDict[pfx+'A0'],0,0,0,0,0]
1261    return A,sigA
1262       
1263def PrintRestraints(cell,SGData,AtPtrs,Atoms,AtLookup,textureData,phaseRest,pFile):
1264    if phaseRest:
1265        Amat = G2lat.cell2AB(cell)[0]
1266        cx,ct,cs = AtPtrs[:3]
1267        names = [['Bond','Bonds'],['Angle','Angles'],['Plane','Planes'],
1268            ['Chiral','Volumes'],['Torsion','Torsions'],['Rama','Ramas'],
1269            ['ChemComp','Sites'],['Texture','HKLs']]
1270        for name,rest in names:
1271            itemRest = phaseRest[name]
1272            if itemRest[rest] and itemRest['Use']:
1273                print >>pFile,'\n  %s %10.3f Use: %s'%(name+' restraint weight factor',itemRest['wtFactor'],str(itemRest['Use']))
1274                if name in ['Bond','Angle','Plane','Chiral']:
1275                    print >>pFile,'     calc       obs      sig   delt/sig  atoms(symOp): '
1276                    for indx,ops,obs,esd in itemRest[rest]:
1277                        try:
1278                            AtNames = G2mth.GetAtomItemsById(Atoms,AtLookup,indx,ct-1)
1279                            AtName = ''
1280                            for i,Aname in enumerate(AtNames):
1281                                AtName += Aname
1282                                if ops[i] == '1':
1283                                    AtName += '-'
1284                                else:
1285                                    AtName += '+('+ops[i]+')-'
1286                            XYZ = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cx,3))
1287                            XYZ = G2mth.getSyXYZ(XYZ,ops,SGData)
1288                            if name == 'Bond':
1289                                calc = G2mth.getRestDist(XYZ,Amat)
1290                            elif name == 'Angle':
1291                                calc = G2mth.getRestAngle(XYZ,Amat)
1292                            elif name == 'Plane':
1293                                calc = G2mth.getRestPlane(XYZ,Amat)
1294                            elif name == 'Chiral':
1295                                calc = G2mth.getRestChiral(XYZ,Amat)
1296                            print >>pFile,' %9.3f %9.3f %8.3f %8.3f   %s'%(calc,obs,esd,(obs-calc)/esd,AtName[:-1])
1297                        except KeyError:
1298                            del itemRest[rest]
1299                elif name in ['Torsion','Rama']:
1300                    print >>pFile,'  atoms(symOp)  calc  obs  sig  delt/sig  torsions: '
1301                    coeffDict = itemRest['Coeff']
1302                    for indx,ops,cofName,esd in itemRest[rest]:
1303                        AtNames = G2mth.GetAtomItemsById(Atoms,AtLookup,indx,ct-1)
1304                        AtName = ''
1305                        for i,Aname in enumerate(AtNames):
1306                            AtName += Aname+'+('+ops[i]+')-'
1307                        XYZ = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cx,3))
1308                        XYZ = G2mth.getSyXYZ(XYZ,ops,SGData)
1309                        if name == 'Torsion':
1310                            tor = G2mth.getRestTorsion(XYZ,Amat)
1311                            restr,calc = G2mth.calcTorsionEnergy(tor,coeffDict[cofName])
1312                            print >>pFile,' %8.3f %8.3f %.3f %8.3f %8.3f %s'%(calc,obs,esd,(obs-calc)/esd,tor,AtName[:-1])
1313                        else:
1314                            phi,psi = G2mth.getRestRama(XYZ,Amat)
1315                            restr,calc = G2mth.calcRamaEnergy(phi,psi,coeffDict[cofName])                               
1316                            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])
1317                elif name == 'ChemComp':
1318                    print >>pFile,'     atoms   mul*frac  factor     prod'
1319                    for indx,factors,obs,esd in itemRest[rest]:
1320                        try:
1321                            atoms = G2mth.GetAtomItemsById(Atoms,AtLookup,indx,ct-1)
1322                            mul = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cs+1))
1323                            frac = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cs-1))
1324                            mulfrac = mul*frac
1325                            calcs = mul*frac*factors
1326                            for iatm,[atom,mf,fr,clc] in enumerate(zip(atoms,mulfrac,factors,calcs)):
1327                                print >>pFile,' %10s %8.3f %8.3f %8.3f'%(atom,mf,fr,clc)
1328                            print >>pFile,' Sum:                   calc: %8.3f obs: %8.3f esd: %8.3f'%(np.sum(calcs),obs,esd)
1329                        except KeyError:
1330                            del itemRest[rest]
1331                elif name == 'Texture' and textureData['Order']:
1332                    Start = False
1333                    SHCoef = textureData['SH Coeff'][1]
1334                    shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
1335                    SamSym = dict(zip(shModels,['0','-1','2/m','mmm']))
1336                    print '    HKL  grid  neg esd  sum neg  num neg use unit?  unit esd  '
1337                    for hkl,grid,esd1,ifesd2,esd2 in itemRest[rest]:
1338                        PH = np.array(hkl)
1339                        phi,beta = G2lat.CrsAng(np.array(hkl),cell,SGData)
1340                        ODFln = G2lat.Flnh(Start,SHCoef,phi,beta,SGData)
1341                        R,P,Z = G2mth.getRestPolefig(ODFln,SamSym[textureData['Model']],grid)
1342                        Z = ma.masked_greater(Z,0.0)
1343                        num = ma.count(Z)
1344                        sum = 0
1345                        if num:
1346                            sum = np.sum(Z)
1347                        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)
1348       
1349def getCellEsd(pfx,SGData,A,covData):
1350    dpr = 180./np.pi
1351    rVsq = G2lat.calc_rVsq(A)
1352    G,g = G2lat.A2Gmat(A)       #get recip. & real metric tensors
1353    cell = np.array(G2lat.Gmat2cell(g))   #real cell
1354    cellst = np.array(G2lat.Gmat2cell(G)) #recip. cell
1355    scos = cosd(cellst[3:6])
1356    ssin = sind(cellst[3:6])
1357    scot = scos/ssin
1358    rcos = cosd(cell[3:6])
1359    rsin = sind(cell[3:6])
1360    rcot = rcos/rsin
1361    RMnames = [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3',pfx+'A4',pfx+'A5']
1362    varyList = covData['varyList']
1363    covMatrix = covData['covMatrix']
1364    vcov = G2mth.getVCov(RMnames,varyList,covMatrix)
1365    Ax = np.array(A)
1366    Ax[3:] /= 2.
1367    drVdA = np.array([Ax[1]*Ax[2]-Ax[5]**2,Ax[0]*Ax[2]-Ax[4]**2,Ax[0]*Ax[1]-Ax[3]**2,
1368        Ax[4]*Ax[5]-Ax[2]*Ax[3],Ax[3]*Ax[5]-Ax[1]*Ax[4],Ax[3]*Ax[4]-Ax[0]*Ax[5]])
1369    srcvlsq = np.inner(drVdA,np.inner(vcov,drVdA.T))
1370    Vol = 1/np.sqrt(rVsq)
1371    sigVol = Vol**3*np.sqrt(srcvlsq)/2.
1372    R123 = Ax[0]*Ax[1]*Ax[2]
1373    dsasdg = np.zeros((3,6))
1374    dadg = np.zeros((6,6))
1375    for i0 in range(3):         #0  1   2
1376        i1 = (i0+1)%3           #1  2   0
1377        i2 = (i1+1)%3           #2  0   1
1378        i3 = 5-i2               #3  5   4
1379        i4 = 5-i1               #4  3   5
1380        i5 = 5-i0               #5  4   3
1381        dsasdg[i0][i1] = 0.5*scot[i0]*scos[i0]/Ax[i1]
1382        dsasdg[i0][i2] = 0.5*scot[i0]*scos[i0]/Ax[i2]
1383        dsasdg[i0][i5] = -scot[i0]/np.sqrt(Ax[i1]*Ax[i2])
1384        denmsq = Ax[i0]*(R123-Ax[i1]*Ax[i4]**2-Ax[i2]*Ax[i3]**2+(Ax[i4]*Ax[i3])**2)
1385        denom = np.sqrt(denmsq)
1386        dadg[i5][i0] = -Ax[i5]/denom-rcos[i0]/denmsq*(R123-0.5*Ax[i1]*Ax[i4]**2-0.5*Ax[i2]*Ax[i3]**2)
1387        dadg[i5][i1] = -0.5*rcos[i0]/denmsq*(Ax[i0]**2*Ax[i2]-Ax[i0]*Ax[i4]**2)
1388        dadg[i5][i2] = -0.5*rcos[i0]/denmsq*(Ax[i0]**2*Ax[i1]-Ax[i0]*Ax[i3]**2)
1389        dadg[i5][i3] = Ax[i4]/denom+rcos[i0]/denmsq*(Ax[i0]*Ax[i2]*Ax[i3]-Ax[i3]*Ax[i4]**2)
1390        dadg[i5][i4] = Ax[i3]/denom+rcos[i0]/denmsq*(Ax[i0]*Ax[i1]*Ax[i4]-Ax[i3]**2*Ax[i4])
1391        dadg[i5][i5] = -Ax[i0]/denom
1392    for i0 in range(3):
1393        i1 = (i0+1)%3
1394        i2 = (i1+1)%3
1395        i3 = 5-i2
1396        for ij in range(6):
1397            dadg[i0][ij] = cell[i0]*(rcot[i2]*dadg[i3][ij]/rsin[i2]-dsasdg[i1][ij]/ssin[i1])
1398            if ij == i0:
1399                dadg[i0][ij] = dadg[i0][ij]-0.5*cell[i0]/Ax[i0]
1400            dadg[i3][ij] = -dadg[i3][ij]*rsin[2-i0]*dpr
1401    sigMat = np.inner(dadg,np.inner(vcov,dadg.T))
1402    var = np.diag(sigMat)
1403    CS = np.where(var>0.,np.sqrt(var),0.)
1404    cellSig = [CS[0],CS[1],CS[2],CS[5],CS[4],CS[3],sigVol]  #exchange sig(alp) & sig(gam) to get in right order
1405    return cellSig           
1406   
1407def SetPhaseData(parmDict,sigDict,Phases,RBIds,covData,RestraintDict=None,pFile=None):
1408   
1409    def PrintAtomsAndSig(General,Atoms,atomsSig):
1410        print >>pFile,'\n Atoms:'
1411        line = '   name      x         y         z      frac   Uiso     U11     U22     U33     U12     U13     U23'
1412        if General['Type'] == 'magnetic':
1413            line += '   Mx     My     Mz'
1414        elif General['Type'] == 'macromolecular':
1415            line = ' res no  residue  chain '+line
1416        print >>pFile,line
1417        if General['Type'] == 'nuclear':
1418            print >>pFile,135*'-'
1419            fmt = {0:'%7s',1:'%7s',3:'%10.5f',4:'%10.5f',5:'%10.5f',6:'%8.3f',10:'%8.5f',
1420                11:'%8.5f',12:'%8.5f',13:'%8.5f',14:'%8.5f',15:'%8.5f',16:'%8.5f'}
1421            noFXsig = {3:[10*' ','%10s'],4:[10*' ','%10s'],5:[10*' ','%10s'],6:[8*' ','%8s']}
1422            for atyp in General['AtomTypes']:       #zero composition
1423                General['NoAtoms'][atyp] = 0.
1424            for i,at in enumerate(Atoms):
1425                General['NoAtoms'][at[1]] += at[6]*float(at[8])     #new composition
1426                name = fmt[0]%(at[0])+fmt[1]%(at[1])+':'
1427                valstr = ' values:'
1428                sigstr = ' sig   :'
1429                for ind in [3,4,5,6]:
1430                    sigind = str(i)+':'+str(ind)
1431                    valstr += fmt[ind]%(at[ind])                   
1432                    if sigind in atomsSig:
1433                        sigstr += fmt[ind]%(atomsSig[sigind])
1434                    else:
1435                        sigstr += noFXsig[ind][1]%(noFXsig[ind][0])
1436                if at[9] == 'I':
1437                    valstr += fmt[10]%(at[10])
1438                    if str(i)+':10' in atomsSig:
1439                        sigstr += fmt[10]%(atomsSig[str(i)+':10'])
1440                    else:
1441                        sigstr += 8*' '
1442                else:
1443                    valstr += 8*' '
1444                    sigstr += 8*' '
1445                    for ind in [11,12,13,14,15,16]:
1446                        sigind = str(i)+':'+str(ind)
1447                        valstr += fmt[ind]%(at[ind])
1448                        if sigind in atomsSig:                       
1449                            sigstr += fmt[ind]%(atomsSig[sigind])
1450                        else:
1451                            sigstr += 8*' '
1452                print >>pFile,name
1453                print >>pFile,valstr
1454                print >>pFile,sigstr
1455               
1456    def PrintRBObjPOAndSig(rbfx,rbsx):
1457        namstr = '  names :'
1458        valstr = '  values:'
1459        sigstr = '  esds  :'
1460        for i,px in enumerate(['Px:','Py:','Pz:']):
1461            name = pfx+rbfx+px+rbsx
1462            namstr += '%12s'%('Pos '+px[1])
1463            valstr += '%12.5f'%(parmDict[name])
1464            if name in sigDict:
1465                sigstr += '%12.5f'%(sigDict[name])
1466            else:
1467                sigstr += 12*' '
1468        for i,po in enumerate(['Oa:','Oi:','Oj:','Ok:']):
1469            name = pfx+rbfx+po+rbsx
1470            namstr += '%12s'%('Ori '+po[1])
1471            valstr += '%12.5f'%(parmDict[name])
1472            if name in sigDict:
1473                sigstr += '%12.5f'%(sigDict[name])
1474            else:
1475                sigstr += 12*' '
1476        print >>pFile,namstr
1477        print >>pFile,valstr
1478        print >>pFile,sigstr
1479       
1480    def PrintRBObjTLSAndSig(rbfx,rbsx,TLS):
1481        namstr = '  names :'
1482        valstr = '  values:'
1483        sigstr = '  esds  :'
1484        if 'N' not in TLS:
1485            print >>pFile,' Thermal motion:'
1486        if 'T' in TLS:
1487            for i,pt in enumerate(['T11:','T22:','T33:','T12:','T13:','T23:']):
1488                name = pfx+rbfx+pt+rbsx
1489                namstr += '%12s'%(pt[:3])
1490                valstr += '%12.5f'%(parmDict[name])
1491                if name in sigDict:
1492                    sigstr += '%12.5f'%(sigDict[name])
1493                else:
1494                    sigstr += 12*' '
1495            print >>pFile,namstr
1496            print >>pFile,valstr
1497            print >>pFile,sigstr
1498        if 'L' in TLS:
1499            namstr = '  names :'
1500            valstr = '  values:'
1501            sigstr = '  esds  :'
1502            for i,pt in enumerate(['L11:','L22:','L33:','L12:','L13:','L23:']):
1503                name = pfx+rbfx+pt+rbsx
1504                namstr += '%12s'%(pt[:3])
1505                valstr += '%12.3f'%(parmDict[name])
1506                if name in sigDict:
1507                    sigstr += '%12.3f'%(sigDict[name])
1508                else:
1509                    sigstr += 12*' '
1510            print >>pFile,namstr
1511            print >>pFile,valstr
1512            print >>pFile,sigstr
1513        if 'S' in TLS:
1514            namstr = '  names :'
1515            valstr = '  values:'
1516            sigstr = '  esds  :'
1517            for i,pt in enumerate(['S12:','S13:','S21:','S23:','S31:','S32:','SAA:','SBB:']):
1518                name = pfx+rbfx+pt+rbsx
1519                namstr += '%12s'%(pt[:3])
1520                valstr += '%12.4f'%(parmDict[name])
1521                if name in sigDict:
1522                    sigstr += '%12.4f'%(sigDict[name])
1523                else:
1524                    sigstr += 12*' '
1525            print >>pFile,namstr
1526            print >>pFile,valstr
1527            print >>pFile,sigstr
1528        if 'U' in TLS:
1529            name = pfx+rbfx+'U:'+rbsx
1530            namstr = '  names :'+'%12s'%('Uiso')
1531            valstr = '  values:'+'%12.5f'%(parmDict[name])
1532            if name in sigDict:
1533                sigstr = '  esds  :'+'%12.5f'%(sigDict[name])
1534            else:
1535                sigstr = '  esds  :'+12*' '
1536            print >>pFile,namstr
1537            print >>pFile,valstr
1538            print >>pFile,sigstr
1539       
1540    def PrintRBObjTorAndSig(rbsx):
1541        namstr = '  names :'
1542        valstr = '  values:'
1543        sigstr = '  esds  :'
1544        nTors = len(RBObj['Torsions'])
1545        if nTors:
1546            print >>pFile,' Torsions:'
1547            for it in range(nTors):
1548                name = pfx+'RBRTr;'+str(it)+':'+rbsx
1549                namstr += '%12s'%('Tor'+str(it))
1550                valstr += '%12.4f'%(parmDict[name])
1551                if name in sigDict:
1552                    sigstr += '%12.4f'%(sigDict[name])
1553            print >>pFile,namstr
1554            print >>pFile,valstr
1555            print >>pFile,sigstr
1556               
1557    def PrintSHtextureAndSig(textureData,SHtextureSig):
1558        print >>pFile,'\n Spherical harmonics texture: Order:' + str(textureData['Order'])
1559        names = ['omega','chi','phi']
1560        namstr = '  names :'
1561        ptstr =  '  values:'
1562        sigstr = '  esds  :'
1563        for name in names:
1564            namstr += '%12s'%(name)
1565            ptstr += '%12.3f'%(textureData['Sample '+name][1])
1566            if 'Sample '+name in SHtextureSig:
1567                sigstr += '%12.3f'%(SHtextureSig['Sample '+name])
1568            else:
1569                sigstr += 12*' '
1570        print >>pFile,namstr
1571        print >>pFile,ptstr
1572        print >>pFile,sigstr
1573        print >>pFile,'\n Texture coefficients:'
1574        namstr = '  names :'
1575        ptstr =  '  values:'
1576        sigstr = '  esds  :'
1577        SHcoeff = textureData['SH Coeff'][1]
1578        for name in SHcoeff:
1579            namstr += '%12s'%(name)
1580            ptstr += '%12.3f'%(SHcoeff[name])
1581            if name in SHtextureSig:
1582                sigstr += '%12.3f'%(SHtextureSig[name])
1583            else:
1584                sigstr += 12*' '
1585        print >>pFile,namstr
1586        print >>pFile,ptstr
1587        print >>pFile,sigstr
1588           
1589    print >>pFile,'\n Phases:'
1590    for phase in Phases:
1591        print >>pFile,' Result for phase: ',phase
1592        Phase = Phases[phase]
1593        General = Phase['General']
1594        SGData = General['SGData']
1595        Atoms = Phase['Atoms']
1596        AtLookup = G2mth.FillAtomLookUp(Atoms)
1597        cell = General['Cell']
1598        pId = Phase['pId']
1599        pfx = str(pId)+'::'
1600        if cell[0]:
1601            A,sigA = cellFill(pfx,SGData,parmDict,sigDict)
1602            cellSig = getCellEsd(pfx,SGData,A,covData)  #includes sigVol
1603            print >>pFile,' Reciprocal metric tensor: '
1604            ptfmt = "%15.9f"
1605            names = ['A11','A22','A33','A12','A13','A23']
1606            namstr = '  names :'
1607            ptstr =  '  values:'
1608            sigstr = '  esds  :'
1609            for name,a,siga in zip(names,A,sigA):
1610                namstr += '%15s'%(name)
1611                ptstr += ptfmt%(a)
1612                if siga:
1613                    sigstr += ptfmt%(siga)
1614                else:
1615                    sigstr += 15*' '
1616            print >>pFile,namstr
1617            print >>pFile,ptstr
1618            print >>pFile,sigstr
1619            cell[1:7] = G2lat.A2cell(A)
1620            cell[7] = G2lat.calc_V(A)
1621            print >>pFile,' New unit cell:'
1622            ptfmt = ["%12.6f","%12.6f","%12.6f","%12.4f","%12.4f","%12.4f","%12.3f"]
1623            names = ['a','b','c','alpha','beta','gamma','Volume']
1624            namstr = '  names :'
1625            ptstr =  '  values:'
1626            sigstr = '  esds  :'
1627            for name,fmt,a,siga in zip(names,ptfmt,cell[1:8],cellSig):
1628                namstr += '%12s'%(name)
1629                ptstr += fmt%(a)
1630                if siga:
1631                    sigstr += fmt%(siga)
1632                else:
1633                    sigstr += 12*' '
1634            print >>pFile,namstr
1635            print >>pFile,ptstr
1636            print >>pFile,sigstr
1637           
1638        General['Mass'] = 0.
1639        if Phase['General'].get('doPawley'):
1640            pawleyRef = Phase['Pawley ref']
1641            for i,refl in enumerate(pawleyRef):
1642                key = pfx+'PWLref:'+str(i)
1643                refl[6] = parmDict[key]
1644                if key in sigDict:
1645                    refl[7] = sigDict[key]
1646                else:
1647                    refl[7] = 0
1648        else:
1649            VRBIds = RBIds['Vector']
1650            RRBIds = RBIds['Residue']
1651            RBModels = Phase['RBModels']
1652            for irb,RBObj in enumerate(RBModels.get('Vector',[])):
1653                jrb = VRBIds.index(RBObj['RBId'])
1654                rbsx = str(irb)+':'+str(jrb)
1655                print >>pFile,' Vector rigid body parameters:'
1656                PrintRBObjPOAndSig('RBV',rbsx)
1657                PrintRBObjTLSAndSig('RBV',rbsx,RBObj['ThermalMotion'][0])
1658            for irb,RBObj in enumerate(RBModels.get('Residue',[])):
1659                jrb = RRBIds.index(RBObj['RBId'])
1660                rbsx = str(irb)+':'+str(jrb)
1661                print >>pFile,' Residue rigid body parameters:'
1662                PrintRBObjPOAndSig('RBR',rbsx)
1663                PrintRBObjTLSAndSig('RBR',rbsx,RBObj['ThermalMotion'][0])
1664                PrintRBObjTorAndSig(rbsx)
1665            atomsSig = {}
1666            if General['Type'] == 'nuclear':        #this needs macromolecular variant!
1667                for i,at in enumerate(Atoms):
1668                    names = {3:pfx+'Ax:'+str(i),4:pfx+'Ay:'+str(i),5:pfx+'Az:'+str(i),6:pfx+'Afrac:'+str(i),
1669                        10:pfx+'AUiso:'+str(i),11:pfx+'AU11:'+str(i),12:pfx+'AU22:'+str(i),13:pfx+'AU33:'+str(i),
1670                        14:pfx+'AU12:'+str(i),15:pfx+'AU13:'+str(i),16:pfx+'AU23:'+str(i)}
1671                    for ind in [3,4,5,6]:
1672                        at[ind] = parmDict[names[ind]]
1673                        if ind in [3,4,5]:
1674                            name = names[ind].replace('A','dA')
1675                        else:
1676                            name = names[ind]
1677                        if name in sigDict:
1678                            atomsSig[str(i)+':'+str(ind)] = sigDict[name]
1679                    if at[9] == 'I':
1680                        at[10] = parmDict[names[10]]
1681                        if names[10] in sigDict:
1682                            atomsSig[str(i)+':10'] = sigDict[names[10]]
1683                    else:
1684                        for ind in [11,12,13,14,15,16]:
1685                            at[ind] = parmDict[names[ind]]
1686                            if names[ind] in sigDict:
1687                                atomsSig[str(i)+':'+str(ind)] = sigDict[names[ind]]
1688                    ind = General['AtomTypes'].index(at[1])
1689                    General['Mass'] += General['AtomMass'][ind]*at[6]*at[8]
1690            PrintAtomsAndSig(General,Atoms,atomsSig)
1691       
1692        textureData = General['SH Texture']   
1693        if textureData['Order']:
1694            SHtextureSig = {}
1695            for name in ['omega','chi','phi']:
1696                aname = pfx+'SH '+name
1697                textureData['Sample '+name][1] = parmDict[aname]
1698                if aname in sigDict:
1699                    SHtextureSig['Sample '+name] = sigDict[aname]
1700            for name in textureData['SH Coeff'][1]:
1701                aname = pfx+name
1702                textureData['SH Coeff'][1][name] = parmDict[aname]
1703                if aname in sigDict:
1704                    SHtextureSig[name] = sigDict[aname]
1705            PrintSHtextureAndSig(textureData,SHtextureSig)
1706        if phase in RestraintDict:
1707            PrintRestraints(cell[1:7],SGData,General['AtomPtrs'],Atoms,AtLookup,
1708                textureData,RestraintDict[phase],pFile)
1709                   
1710################################################################################
1711##### Histogram & Phase data
1712################################################################################       
1713                   
1714def GetHistogramPhaseData(Phases,Histograms,Print=True,pFile=None):
1715   
1716    def PrintSize(hapData):
1717        if hapData[0] in ['isotropic','uniaxial']:
1718            line = '\n Size model    : %9s'%(hapData[0])
1719            line += ' equatorial:'+'%12.3f'%(hapData[1][0])+' Refine? '+str(hapData[2][0])
1720            if hapData[0] == 'uniaxial':
1721                line += ' axial:'+'%12.3f'%(hapData[1][1])+' Refine? '+str(hapData[2][1])
1722            line += '\n\t LG mixing coeff.: %12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
1723            print >>pFile,line
1724        else:
1725            print >>pFile,'\n Size model    : %s'%(hapData[0])+ \
1726                '\n\t LG mixing coeff.:%12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
1727            Snames = ['S11','S22','S33','S12','S13','S23']
1728            ptlbls = ' names :'
1729            ptstr =  ' values:'
1730            varstr = ' refine:'
1731            for i,name in enumerate(Snames):
1732                ptlbls += '%12s' % (name)
1733                ptstr += '%12.6f' % (hapData[4][i])
1734                varstr += '%12s' % (str(hapData[5][i]))
1735            print >>pFile,ptlbls
1736            print >>pFile,ptstr
1737            print >>pFile,varstr
1738       
1739    def PrintMuStrain(hapData,SGData):
1740        if hapData[0] in ['isotropic','uniaxial']:
1741            line = '\n Mustrain model: %9s'%(hapData[0])
1742            line += ' equatorial:'+'%12.1f'%(hapData[1][0])+' Refine? '+str(hapData[2][0])
1743            if hapData[0] == 'uniaxial':
1744                line += ' axial:'+'%12.1f'%(hapData[1][1])+' Refine? '+str(hapData[2][1])
1745            line +='\n\t LG mixing coeff.:%12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
1746            print >>pFile,line
1747        else:
1748            print >>pFile,'\n Mustrain model: %s'%(hapData[0])+ \
1749                '\n\t LG mixing coeff.:%12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
1750            Snames = G2spc.MustrainNames(SGData)
1751            ptlbls = ' names :'
1752            ptstr =  ' values:'
1753            varstr = ' refine:'
1754            for i,name in enumerate(Snames):
1755                ptlbls += '%12s' % (name)
1756                ptstr += '%12.6f' % (hapData[4][i])
1757                varstr += '%12s' % (str(hapData[5][i]))
1758            print >>pFile,ptlbls
1759            print >>pFile,ptstr
1760            print >>pFile,varstr
1761
1762    def PrintHStrain(hapData,SGData):
1763        print >>pFile,'\n Hydrostatic/elastic strain: '
1764        Hsnames = G2spc.HStrainNames(SGData)
1765        ptlbls = ' names :'
1766        ptstr =  ' values:'
1767        varstr = ' refine:'
1768        for i,name in enumerate(Hsnames):
1769            ptlbls += '%12s' % (name)
1770            ptstr += '%12.6f' % (hapData[0][i])
1771            varstr += '%12s' % (str(hapData[1][i]))
1772        print >>pFile,ptlbls
1773        print >>pFile,ptstr
1774        print >>pFile,varstr
1775
1776    def PrintSHPO(hapData):
1777        print >>pFile,'\n Spherical harmonics preferred orientation: Order:' + \
1778            str(hapData[4])+' Refine? '+str(hapData[2])
1779        ptlbls = ' names :'
1780        ptstr =  ' values:'
1781        for item in hapData[5]:
1782            ptlbls += '%12s'%(item)
1783            ptstr += '%12.3f'%(hapData[5][item]) 
1784        print >>pFile,ptlbls
1785        print >>pFile,ptstr
1786   
1787    def PrintBabinet(hapData):
1788        print >>pFile,'\n Babinet form factor modification: '
1789        ptlbls = ' names :'
1790        ptstr =  ' values:'
1791        varstr = ' refine:'
1792        for name in ['BabA','BabU']:
1793            ptlbls += '%12s' % (name)
1794            ptstr += '%12.6f' % (hapData[name][0])
1795            varstr += '%12s' % (str(hapData[name][1]))
1796        print >>pFile,ptlbls
1797        print >>pFile,ptstr
1798        print >>pFile,varstr
1799       
1800   
1801    hapDict = {}
1802    hapVary = []
1803    controlDict = {}
1804    poType = {}
1805    poAxes = {}
1806    spAxes = {}
1807    spType = {}
1808   
1809    for phase in Phases:
1810        HistoPhase = Phases[phase]['Histograms']
1811        SGData = Phases[phase]['General']['SGData']
1812        cell = Phases[phase]['General']['Cell'][1:7]
1813        A = G2lat.cell2A(cell)
1814        pId = Phases[phase]['pId']
1815        histoList = HistoPhase.keys()
1816        histoList.sort()
1817        for histogram in histoList:
1818            try:
1819                Histogram = Histograms[histogram]
1820            except KeyError:                       
1821                #skip if histogram not included e.g. in a sequential refinement
1822                continue
1823            hapData = HistoPhase[histogram]
1824            hId = Histogram['hId']
1825            if 'PWDR' in histogram:
1826                limits = Histogram['Limits'][1]
1827                inst = Histogram['Instrument Parameters'][0]
1828                Zero = inst['Zero'][1]
1829                if 'C' in inst['Type'][1]:
1830                    try:
1831                        wave = inst['Lam'][1]
1832                    except KeyError:
1833                        wave = inst['Lam1'][1]
1834                    dmin = wave/(2.0*sind(limits[1]/2.0))
1835                pfx = str(pId)+':'+str(hId)+':'
1836                for item in ['Scale','Extinction']:
1837                    hapDict[pfx+item] = hapData[item][0]
1838                    if hapData[item][1]:
1839                        hapVary.append(pfx+item)
1840                names = G2spc.HStrainNames(SGData)
1841                for i,name in enumerate(names):
1842                    hapDict[pfx+name] = hapData['HStrain'][0][i]
1843                    if hapData['HStrain'][1][i]:
1844                        hapVary.append(pfx+name)
1845                controlDict[pfx+'poType'] = hapData['Pref.Ori.'][0]
1846                if hapData['Pref.Ori.'][0] == 'MD':
1847                    hapDict[pfx+'MD'] = hapData['Pref.Ori.'][1]
1848                    controlDict[pfx+'MDAxis'] = hapData['Pref.Ori.'][3]
1849                    if hapData['Pref.Ori.'][2]:
1850                        hapVary.append(pfx+'MD')
1851                else:                           #'SH' spherical harmonics
1852                    controlDict[pfx+'SHord'] = hapData['Pref.Ori.'][4]
1853                    controlDict[pfx+'SHncof'] = len(hapData['Pref.Ori.'][5])
1854                    for item in hapData['Pref.Ori.'][5]:
1855                        hapDict[pfx+item] = hapData['Pref.Ori.'][5][item]
1856                        if hapData['Pref.Ori.'][2]:
1857                            hapVary.append(pfx+item)
1858                for item in ['Mustrain','Size']:
1859                    controlDict[pfx+item+'Type'] = hapData[item][0]
1860                    hapDict[pfx+item+';mx'] = hapData[item][1][2]
1861                    if hapData[item][2][2]:
1862                        hapVary.append(pfx+item+';mx')
1863                    if hapData[item][0] in ['isotropic','uniaxial']:
1864                        hapDict[pfx+item+';i'] = hapData[item][1][0]
1865                        if hapData[item][2][0]:
1866                            hapVary.append(pfx+item+';i')
1867                        if hapData[item][0] == 'uniaxial':
1868                            controlDict[pfx+item+'Axis'] = hapData[item][3]
1869                            hapDict[pfx+item+';a'] = hapData[item][1][1]
1870                            if hapData[item][2][1]:
1871                                hapVary.append(pfx+item+';a')
1872                    else:       #generalized for mustrain or ellipsoidal for size
1873                        Nterms = len(hapData[item][4])
1874                        if item == 'Mustrain':
1875                            names = G2spc.MustrainNames(SGData)
1876                            pwrs = []
1877                            for name in names:
1878                                h,k,l = name[1:]
1879                                pwrs.append([int(h),int(k),int(l)])
1880                            controlDict[pfx+'MuPwrs'] = pwrs
1881                        for i in range(Nterms):
1882                            sfx = ':'+str(i)
1883                            hapDict[pfx+item+sfx] = hapData[item][4][i]
1884                            if hapData[item][5][i]:
1885                                hapVary.append(pfx+item+sfx)
1886                for bab in ['BabA','BabU']:
1887                    hapDict[pfx+bab] = hapData['Babinet'][bab][0]
1888                    if hapData['Babinet'][bab][1]:
1889                        hapVary.append(pfx+bab)
1890                               
1891                if Print: 
1892                    print >>pFile,'\n Phase: ',phase,' in histogram: ',histogram
1893                    print >>pFile,135*'-'
1894                    print >>pFile,' Phase fraction  : %10.4f'%(hapData['Scale'][0]),' Refine?',hapData['Scale'][1]
1895                    print >>pFile,' Extinction coeff: %10.4f'%(hapData['Extinction'][0]),' Refine?',hapData['Extinction'][1]
1896                    if hapData['Pref.Ori.'][0] == 'MD':
1897                        Ax = hapData['Pref.Ori.'][3]
1898                        print >>pFile,' March-Dollase PO: %10.4f'%(hapData['Pref.Ori.'][1]),' Refine?',hapData['Pref.Ori.'][2], \
1899                            ' Axis: %d %d %d'%(Ax[0],Ax[1],Ax[2])
1900                    else: #'SH' for spherical harmonics
1901                        PrintSHPO(hapData['Pref.Ori.'])
1902                    PrintSize(hapData['Size'])
1903                    PrintMuStrain(hapData['Mustrain'],SGData)
1904                    PrintHStrain(hapData['HStrain'],SGData)
1905                    if hapData['Babinet']['BabA'][0]:
1906                        PrintBabinet(hapData['Babinet'])
1907                HKLd = np.array(G2lat.GenHLaue(dmin,SGData,A))
1908                refList = []
1909                for h,k,l,d in HKLd:
1910                    ext,mul,Uniq,phi = G2spc.GenHKLf([h,k,l],SGData)
1911                    mul *= 2      # for powder overlap of Friedel pairs
1912                    if ext:
1913                        continue
1914                    if 'C' in inst['Type'][0]:
1915                        pos = 2.0*asind(wave/(2.0*d))+Zero
1916                        if limits[0] < pos < limits[1]:
1917                            refList.append([h,k,l,mul,d,pos,0.0,0.0,0.0,0.0,0.0,Uniq,phi,0.0,{}])
1918                            #last item should contain form factor values by atom type
1919                    else:
1920                        raise ValueError 
1921                Histogram['Reflection Lists'][phase] = refList
1922            elif 'HKLF' in histogram:
1923                inst = Histogram['Instrument Parameters'][0]
1924                hId = Histogram['hId']
1925                hfx = ':%d:'%(hId)
1926                for item in inst:
1927                    hapDict[hfx+item] = inst[item][1]
1928                pfx = str(pId)+':'+str(hId)+':'
1929                hapDict[pfx+'Scale'] = hapData['Scale'][0]
1930                if hapData['Scale'][1]:
1931                    hapVary.append(pfx+'Scale')
1932                               
1933                extApprox,extType,extParms = hapData['Extinction']
1934                controlDict[pfx+'EType'] = extType
1935                controlDict[pfx+'EApprox'] = extApprox
1936                controlDict[pfx+'Tbar'] = extParms['Tbar']
1937                controlDict[pfx+'Cos2TM'] = extParms['Cos2TM']
1938                if 'Primary' in extType:
1939                    Ekey = ['Ep',]
1940                elif 'I & II' in extType:
1941                    Ekey = ['Eg','Es']
1942                elif 'Secondary Type II' == extType:
1943                    Ekey = ['Es',]
1944                elif 'Secondary Type I' == extType:
1945                    Ekey = ['Eg',]
1946                else:   #'None'
1947                    Ekey = []
1948                for eKey in Ekey:
1949                    hapDict[pfx+eKey] = extParms[eKey][0]
1950                    if extParms[eKey][1]:
1951                        hapVary.append(pfx+eKey)
1952                for bab in ['BabA','BabU']:
1953                    hapDict[pfx+bab] = hapData['Babinet'][bab][0]
1954                    if hapData['Babinet'][bab][1]:
1955                        hapVary.append(pfx+bab)
1956                if Print: 
1957                    print >>pFile,'\n Phase: ',phase,' in histogram: ',histogram
1958                    print >>pFile,135*'-'
1959                    print >>pFile,' Scale factor     : %10.4f'%(hapData['Scale'][0]),' Refine?',hapData['Scale'][1]
1960                    if extType != 'None':
1961                        print >>pFile,' Extinction  Type: %15s'%(extType),' approx: %10s'%(extApprox),' tbar: %6.3f'%(extParms['Tbar'])
1962                        text = ' Parameters       :'
1963                        for eKey in Ekey:
1964                            text += ' %4s : %10.3e Refine? '%(eKey,extParms[eKey][0])+str(extParms[eKey][1])
1965                        print >>pFile,text
1966                    if hapData['Babinet']['BabA'][0]:
1967                        PrintBabinet(hapData['Babinet'])
1968                Histogram['Reflection Lists'] = phase       
1969               
1970    return hapVary,hapDict,controlDict
1971   
1972def SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms,Print=True,pFile=None):
1973   
1974    def PrintSizeAndSig(hapData,sizeSig):
1975        line = '\n Size model:     %9s'%(hapData[0])
1976        refine = False
1977        if hapData[0] in ['isotropic','uniaxial']:
1978            line += ' equatorial:%12.4f'%(hapData[1][0])
1979            if sizeSig[0][0]:
1980                line += ', sig:%8.4f'%(sizeSig[0][0])
1981                refine = True
1982            if hapData[0] == 'uniaxial':
1983                line += ' axial:%12.4f'%(hapData[1][1])
1984                if sizeSig[0][1]:
1985                    refine = True
1986                    line += ', sig:%8.4f'%(sizeSig[0][1])
1987            line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
1988            if sizeSig[0][2]:
1989                refine = True
1990                line += ', sig:%8.4f'%(sizeSig[0][2])
1991            if refine:
1992                print >>pFile,line
1993        else:
1994            line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
1995            if sizeSig[0][2]:
1996                refine = True
1997                line += ', sig:%8.4f'%(sizeSig[0][2])
1998            Snames = ['S11','S22','S33','S12','S13','S23']
1999            ptlbls = ' name  :'
2000            ptstr =  ' value :'
2001            sigstr = ' sig   :'
2002            for i,name in enumerate(Snames):
2003                ptlbls += '%12s' % (name)
2004                ptstr += '%12.6f' % (hapData[4][i])
2005                if sizeSig[1][i]:
2006                    refine = True
2007                    sigstr += '%12.6f' % (sizeSig[1][i])
2008                else:
2009                    sigstr += 12*' '
2010            if refine:
2011                print >>pFile,line
2012                print >>pFile,ptlbls
2013                print >>pFile,ptstr
2014                print >>pFile,sigstr
2015       
2016    def PrintMuStrainAndSig(hapData,mustrainSig,SGData):
2017        line = '\n Mustrain model: %9s'%(hapData[0])
2018        refine = False
2019        if hapData[0] in ['isotropic','uniaxial']:
2020            line += ' equatorial:%12.1f'%(hapData[1][0])
2021            if mustrainSig[0][0]:
2022                line += ', sig:%8.1f'%(mustrainSig[0][0])
2023                refine = True
2024            if hapData[0] == 'uniaxial':
2025                line += ' axial:%12.1f'%(hapData[1][1])
2026                if mustrainSig[0][1]:
2027                     line += ', sig:%8.1f'%(mustrainSig[0][1])
2028            line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
2029            if mustrainSig[0][2]:
2030                refine = True
2031                line += ', sig:%8.3f'%(mustrainSig[0][2])
2032            if refine:
2033                print >>pFile,line
2034        else:
2035            line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
2036            if mustrainSig[0][2]:
2037                refine = True
2038                line += ', sig:%8.3f'%(mustrainSig[0][2])
2039            Snames = G2spc.MustrainNames(SGData)
2040            ptlbls = ' name  :'
2041            ptstr =  ' value :'
2042            sigstr = ' sig   :'
2043            for i,name in enumerate(Snames):
2044                ptlbls += '%12s' % (name)
2045                ptstr += '%12.6f' % (hapData[4][i])
2046                if mustrainSig[1][i]:
2047                    refine = True
2048                    sigstr += '%12.6f' % (mustrainSig[1][i])
2049                else:
2050                    sigstr += 12*' '
2051            if refine:
2052                print >>pFile,line
2053                print >>pFile,ptlbls
2054                print >>pFile,ptstr
2055                print >>pFile,sigstr
2056           
2057    def PrintHStrainAndSig(hapData,strainSig,SGData):
2058        Hsnames = G2spc.HStrainNames(SGData)
2059        ptlbls = ' name  :'
2060        ptstr =  ' value :'
2061        sigstr = ' sig   :'
2062        refine = False
2063        for i,name in enumerate(Hsnames):
2064            ptlbls += '%12s' % (name)
2065            ptstr += '%12.6g' % (hapData[0][i])
2066            if name in strainSig:
2067                refine = True
2068                sigstr += '%12.6g' % (strainSig[name])
2069            else:
2070                sigstr += 12*' '
2071        if refine:
2072            print >>pFile,'\n Hydrostatic/elastic strain: '
2073            print >>pFile,ptlbls
2074            print >>pFile,ptstr
2075            print >>pFile,sigstr
2076       
2077    def PrintSHPOAndSig(pfx,hapData,POsig):
2078        print >>pFile,'\n Spherical harmonics preferred orientation: Order:'+str(hapData[4])
2079        ptlbls = ' names :'
2080        ptstr =  ' values:'
2081        sigstr = ' sig   :'
2082        for item in hapData[5]:
2083            ptlbls += '%12s'%(item)
2084            ptstr += '%12.3f'%(hapData[5][item])
2085            if pfx+item in POsig:
2086                sigstr += '%12.3f'%(POsig[pfx+item])
2087            else:
2088                sigstr += 12*' ' 
2089        print >>pFile,ptlbls
2090        print >>pFile,ptstr
2091        print >>pFile,sigstr
2092       
2093    def PrintExtAndSig(pfx,hapData,ScalExtSig):
2094        print >>pFile,'\n Single crystal extinction: Type: ',hapData[0],' Approx: ',hapData[1]
2095        text = ''
2096        for item in hapData[2]:
2097            if pfx+item in ScalExtSig:
2098                text += '       %s: '%(item)
2099                text += '%12.2e'%(hapData[2][item][0])
2100                if pfx+item in ScalExtSig:
2101                    text += ' sig: %12.2e'%(ScalExtSig[pfx+item])
2102        print >>pFile,text       
2103       
2104    def PrintBabinetAndSig(pfx,hapData,BabSig):
2105        print >>pFile,'\n Babinet form factor modification: '
2106        ptlbls = ' names :'
2107        ptstr =  ' values:'
2108        sigstr = ' sig   :'
2109        for item in hapData:
2110            ptlbls += '%12s'%(item)
2111            ptstr += '%12.3f'%(hapData[item][0])
2112            if pfx+item in BabSig:
2113                sigstr += '%12.3f'%(BabSig[pfx+item])
2114            else:
2115                sigstr += 12*' ' 
2116        print >>pFile,ptlbls
2117        print >>pFile,ptstr
2118        print >>pFile,sigstr
2119   
2120    PhFrExtPOSig = {}
2121    SizeMuStrSig = {}
2122    ScalExtSig = {}
2123    BabSig = {}
2124    wtFrSum = {}
2125    for phase in Phases:
2126        HistoPhase = Phases[phase]['Histograms']
2127        General = Phases[phase]['General']
2128        SGData = General['SGData']
2129        pId = Phases[phase]['pId']
2130        histoList = HistoPhase.keys()
2131        histoList.sort()
2132        for histogram in histoList:
2133            try:
2134                Histogram = Histograms[histogram]
2135            except KeyError:                       
2136                #skip if histogram not included e.g. in a sequential refinement
2137                continue
2138            hapData = HistoPhase[histogram]
2139            hId = Histogram['hId']
2140            pfx = str(pId)+':'+str(hId)+':'
2141            if hId not in wtFrSum:
2142                wtFrSum[hId] = 0.
2143            if 'PWDR' in histogram:
2144                for item in ['Scale','Extinction']:
2145                    hapData[item][0] = parmDict[pfx+item]
2146                    if pfx+item in sigDict:
2147                        PhFrExtPOSig.update({pfx+item:sigDict[pfx+item],})
2148                wtFrSum[hId] += hapData['Scale'][0]*General['Mass']
2149                if hapData['Pref.Ori.'][0] == 'MD':
2150                    hapData['Pref.Ori.'][1] = parmDict[pfx+'MD']
2151                    if pfx+'MD' in sigDict:
2152                        PhFrExtPOSig.update({pfx+'MD':sigDict[pfx+'MD'],})
2153                else:                           #'SH' spherical harmonics
2154                    for item in hapData['Pref.Ori.'][5]:
2155                        hapData['Pref.Ori.'][5][item] = parmDict[pfx+item]
2156                        if pfx+item in sigDict:
2157                            PhFrExtPOSig.update({pfx+item:sigDict[pfx+item],})
2158                SizeMuStrSig.update({pfx+'Mustrain':[[0,0,0],[0 for i in range(len(hapData['Mustrain'][4]))]],
2159                    pfx+'Size':[[0,0,0],[0 for i in range(len(hapData['Size'][4]))]],
2160                    pfx+'HStrain':{}})                 
2161                for item in ['Mustrain','Size']:
2162                    hapData[item][1][2] = parmDict[pfx+item+';mx']
2163                    hapData[item][1][2] = min(1.,max(0.1,hapData[item][1][2]))
2164                    if pfx+item+';mx' in sigDict:
2165                        SizeMuStrSig[pfx+item][0][2] = sigDict[pfx+item+';mx']
2166                    if hapData[item][0] in ['isotropic','uniaxial']:                   
2167                        hapData[item][1][0] = parmDict[pfx+item+';i']
2168                        if item == 'Size':
2169                            hapData[item][1][0] = min(10.,max(0.001,hapData[item][1][0]))
2170                        if pfx+item+';i' in sigDict: 
2171                            SizeMuStrSig[pfx+item][0][0] = sigDict[pfx+item+';i']
2172                        if hapData[item][0] == 'uniaxial':
2173                            hapData[item][1][1] = parmDict[pfx+item+';a']
2174                            if item == 'Size':
2175                                hapData[item][1][1] = min(10.,max(0.001,hapData[item][1][1]))                       
2176                            if pfx+item+';a' in sigDict:
2177                                SizeMuStrSig[pfx+item][0][1] = sigDict[pfx+item+';a']
2178                    else:       #generalized for mustrain or ellipsoidal for size
2179                        Nterms = len(hapData[item][4])
2180                        for i in range(Nterms):
2181                            sfx = ':'+str(i)
2182                            hapData[item][4][i] = parmDict[pfx+item+sfx]
2183                            if pfx+item+sfx in sigDict:
2184                                SizeMuStrSig[pfx+item][1][i] = sigDict[pfx+item+sfx]
2185                names = G2spc.HStrainNames(SGData)
2186                for i,name in enumerate(names):
2187                    hapData['HStrain'][0][i] = parmDict[pfx+name]
2188                    if pfx+name in sigDict:
2189                        SizeMuStrSig[pfx+'HStrain'][name] = sigDict[pfx+name]
2190                for name in ['BabA','BabU']:
2191                    hapData['Babinet'][name][0] = parmDict[pfx+name]
2192                    if pfx+name in sigDict:
2193                        BabSig[pfx+name] = sigDict[pfx+name]               
2194               
2195            elif 'HKLF' in histogram:
2196                for item in ['Scale',]:
2197                    if parmDict.get(pfx+item):
2198                        hapData[item][0] = parmDict[pfx+item]
2199                        if pfx+item in sigDict:
2200                            ScalExtSig[pfx+item] = sigDict[pfx+item]
2201                for item in ['Ep','Eg','Es']:
2202                    if parmDict.get(pfx+item):
2203                        hapData['Extinction'][2][item][0] = parmDict[pfx+item]
2204                        if pfx+item in sigDict:
2205                            ScalExtSig[pfx+item] = sigDict[pfx+item]
2206                for name in ['BabA','BabU']:
2207                    hapData['Babinet'][name][0] = parmDict[pfx+name]
2208                    if pfx+name in sigDict:
2209                        BabSig[pfx+name] = sigDict[pfx+name]               
2210
2211    if Print:
2212        for phase in Phases:
2213            HistoPhase = Phases[phase]['Histograms']
2214            General = Phases[phase]['General']
2215            SGData = General['SGData']
2216            pId = Phases[phase]['pId']
2217            histoList = HistoPhase.keys()
2218            histoList.sort()
2219            for histogram in histoList:
2220                try:
2221                    Histogram = Histograms[histogram]
2222                except KeyError:                       
2223                    #skip if histogram not included e.g. in a sequential refinement
2224                    continue
2225                print >>pFile,'\n Phase: ',phase,' in histogram: ',histogram
2226                print >>pFile,130*'-'
2227                hapData = HistoPhase[histogram]
2228                hId = Histogram['hId']
2229                pfx = str(pId)+':'+str(hId)+':'
2230                if 'PWDR' in histogram:
2231                    print >>pFile,' Final refinement RF, RF^2 = %.2f%%, %.2f%% on %d reflections'   \
2232                        %(Histogram[pfx+'Rf'],Histogram[pfx+'Rf^2'],Histogram[pfx+'Nref'])
2233               
2234                    if pfx+'Scale' in PhFrExtPOSig:
2235                        wtFr = hapData['Scale'][0]*General['Mass']/wtFrSum[hId]
2236                        sigwtFr = PhFrExtPOSig[pfx+'Scale']*wtFr/hapData['Scale'][0]
2237                        print >>pFile,' Phase fraction  : %10.5f, sig %10.5f Weight fraction  : %8.5f, sig %10.5f' \
2238                            %(hapData['Scale'][0],PhFrExtPOSig[pfx+'Scale'],wtFr,sigwtFr)
2239                    if pfx+'Extinction' in PhFrExtPOSig:
2240                        print >>pFile,' Extinction coeff: %10.4f, sig %10.4f'%(hapData['Extinction'][0],PhFrExtPOSig[pfx+'Extinction'])
2241                    if hapData['Pref.Ori.'][0] == 'MD':
2242                        if pfx+'MD' in PhFrExtPOSig:
2243                            print >>pFile,' March-Dollase PO: %10.4f, sig %10.4f'%(hapData['Pref.Ori.'][1],PhFrExtPOSig[pfx+'MD'])
2244                    else:
2245                        PrintSHPOAndSig(pfx,hapData['Pref.Ori.'],PhFrExtPOSig)
2246                    PrintSizeAndSig(hapData['Size'],SizeMuStrSig[pfx+'Size'])
2247                    PrintMuStrainAndSig(hapData['Mustrain'],SizeMuStrSig[pfx+'Mustrain'],SGData)
2248                    PrintHStrainAndSig(hapData['HStrain'],SizeMuStrSig[pfx+'HStrain'],SGData)
2249                    if len(BabSig):
2250                        PrintBabinetAndSig(pfx,hapData['Babinet'],BabSig)
2251                   
2252                elif 'HKLF' in histogram:
2253                    print >>pFile,' Final refinement RF, RF^2 = %.2f%%, %.2f%% on %d reflections'   \
2254                        %(Histogram[pfx+'Rf'],Histogram[pfx+'Rf^2'],Histogram[pfx+'Nref'])
2255                    print >>pFile,' HKLF histogram weight factor = ','%.3f'%(Histogram['wtFactor'])
2256                    if pfx+'Scale' in ScalExtSig:
2257                        print >>pFile,' Scale factor : %10.4f, sig %10.4f'%(hapData['Scale'][0],ScalExtSig[pfx+'Scale'])
2258                    if hapData['Extinction'][0] != 'None':
2259                        PrintExtAndSig(pfx,hapData['Extinction'],ScalExtSig)
2260                    if len(BabSig):
2261                        PrintBabinetAndSig(pfx,hapData['Babinet'],BabSig)
2262
2263################################################################################
2264##### Histogram data
2265################################################################################       
2266                   
2267def GetHistogramData(Histograms,Print=True,pFile=None):
2268   
2269    def GetBackgroundParms(hId,Background):
2270        Back = Background[0]
2271        DebyePeaks = Background[1]
2272        bakType,bakFlag = Back[:2]
2273        backVals = Back[3:]
2274        backNames = [':'+str(hId)+':Back:'+str(i) for i in range(len(backVals))]
2275        backDict = dict(zip(backNames,backVals))
2276        backVary = []
2277        if bakFlag:
2278            backVary = backNames
2279        backDict[':'+str(hId)+':nDebye'] = DebyePeaks['nDebye']
2280        backDict[':'+str(hId)+':nPeaks'] = DebyePeaks['nPeaks']
2281        debyeDict = {}
2282        debyeList = []
2283        for i in range(DebyePeaks['nDebye']):
2284            debyeNames = [':'+str(hId)+':DebyeA:'+str(i),':'+str(hId)+':DebyeR:'+str(i),':'+str(hId)+':DebyeU:'+str(i)]
2285            debyeDict.update(dict(zip(debyeNames,DebyePeaks['debyeTerms'][i][::2])))
2286            debyeList += zip(debyeNames,DebyePeaks['debyeTerms'][i][1::2])
2287        debyeVary = []
2288        for item in debyeList:
2289            if item[1]:
2290                debyeVary.append(item[0])
2291        backDict.update(debyeDict)
2292        backVary += debyeVary
2293        peakDict = {}
2294        peakList = []
2295        for i in range(DebyePeaks['nPeaks']):
2296            peakNames = [':'+str(hId)+':BkPkpos:'+str(i),':'+str(hId)+ \
2297                ':BkPkint:'+str(i),':'+str(hId)+':BkPksig:'+str(i),':'+str(hId)+':BkPkgam:'+str(i)]
2298            peakDict.update(dict(zip(peakNames,DebyePeaks['peaksList'][i][::2])))
2299            peakList += zip(peakNames,DebyePeaks['peaksList'][i][1::2])
2300        peakVary = []
2301        for item in peakList:
2302            if item[1]:
2303                peakVary.append(item[0])
2304        backDict.update(peakDict)
2305        backVary += peakVary
2306        return bakType,backDict,backVary           
2307       
2308    def GetInstParms(hId,Inst):     
2309        dataType = Inst['Type'][0]
2310        instDict = {}
2311        insVary = []
2312        pfx = ':'+str(hId)+':'
2313        insKeys = Inst.keys()
2314        insKeys.sort()
2315        for item in insKeys:
2316            insName = pfx+item
2317            instDict[insName] = Inst[item][1]
2318            if Inst[item][2]:
2319                insVary.append(insName)
2320#        instDict[pfx+'X'] = max(instDict[pfx+'X'],0.001)
2321#        instDict[pfx+'Y'] = max(instDict[pfx+'Y'],0.001)
2322        instDict[pfx+'SH/L'] = max(instDict[pfx+'SH/L'],0.0005)
2323        return dataType,instDict,insVary
2324       
2325    def GetSampleParms(hId,Sample):
2326        sampVary = []
2327        hfx = ':'+str(hId)+':'       
2328        sampDict = {hfx+'Gonio. radius':Sample['Gonio. radius'],hfx+'Omega':Sample['Omega'],
2329            hfx+'Chi':Sample['Chi'],hfx+'Phi':Sample['Phi']}
2330        Type = Sample['Type']
2331        if 'Bragg' in Type:             #Bragg-Brentano
2332            for item in ['Scale','Shift','Transparency']:       #surface roughness?, diffuse scattering?
2333                sampDict[hfx+item] = Sample[item][0]
2334                if Sample[item][1]:
2335                    sampVary.append(hfx+item)
2336        elif 'Debye' in Type:        #Debye-Scherrer
2337            for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
2338                sampDict[hfx+item] = Sample[item][0]
2339                if Sample[item][1]:
2340                    sampVary.append(hfx+item)
2341        return Type,sampDict,sampVary
2342       
2343    def PrintBackground(Background):
2344        Back = Background[0]
2345        DebyePeaks = Background[1]
2346        print >>pFile,'\n Background function: ',Back[0],' Refine?',bool(Back[1])
2347        line = ' Coefficients: '
2348        for i,back in enumerate(Back[3:]):
2349            line += '%10.3f'%(back)
2350            if i and not i%10:
2351                line += '\n'+15*' '
2352        print >>pFile,line
2353        if DebyePeaks['nDebye']:
2354            print >>pFile,'\n Debye diffuse scattering coefficients'
2355            parms = ['DebyeA','DebyeR','DebyeU']
2356            line = ' names :  '
2357            for parm in parms:
2358                line += '%8s refine?'%(parm)
2359            print >>pFile,line
2360            for j,term in enumerate(DebyePeaks['debyeTerms']):
2361                line = ' term'+'%2d'%(j)+':'
2362                for i in range(3):
2363                    line += '%10.3f %5s'%(term[2*i],bool(term[2*i+1]))                   
2364                print >>pFile,line
2365        if DebyePeaks['nPeaks']:
2366            print >>pFile,'\n Single peak coefficients'
2367            parms =    ['BkPkpos','BkPkint','BkPksig','BkPkgam']
2368            line = ' names :  '
2369            for parm in parms:
2370                line += '%8s refine?'%(parm)
2371            print >>pFile,line
2372            for j,term in enumerate(DebyePeaks['peaksList']):
2373                line = ' peak'+'%2d'%(j)+':'
2374                for i in range(4):
2375                    line += '%10.3f %5s'%(term[2*i],bool(term[2*i+1]))                   
2376                print >>pFile,line
2377       
2378    def PrintInstParms(Inst):
2379        print >>pFile,'\n Instrument Parameters:'
2380        ptlbls = ' name  :'
2381        ptstr =  ' value :'
2382        varstr = ' refine:'
2383        insKeys = Inst.keys()
2384        insKeys.sort()
2385        for item in insKeys:
2386            if item != 'Type':
2387                ptlbls += '%12s' % (item)
2388                ptstr += '%12.6f' % (Inst[item][1])
2389                if item in ['Lam1','Lam2','Azimuth']:
2390                    varstr += 12*' '
2391                else:
2392                    varstr += '%12s' % (str(bool(Inst[item][2])))
2393        print >>pFile,ptlbls
2394        print >>pFile,ptstr
2395        print >>pFile,varstr
2396       
2397    def PrintSampleParms(Sample):
2398        print >>pFile,'\n Sample Parameters:'
2399        print >>pFile,' Goniometer omega = %.2f, chi = %.2f, phi = %.2f'% \
2400            (Sample['Omega'],Sample['Chi'],Sample['Phi'])
2401        ptlbls = ' name  :'
2402        ptstr =  ' value :'
2403        varstr = ' refine:'
2404        if 'Bragg' in Sample['Type']:
2405            for item in ['Scale','Shift','Transparency']:
2406                ptlbls += '%14s'%(item)
2407                ptstr += '%14.4f'%(Sample[item][0])
2408                varstr += '%14s'%(str(bool(Sample[item][1])))
2409           
2410        elif 'Debye' in Type:        #Debye-Scherrer
2411            for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
2412                ptlbls += '%14s'%(item)
2413                ptstr += '%14.4f'%(Sample[item][0])
2414                varstr += '%14s'%(str(bool(Sample[item][1])))
2415
2416        print >>pFile,ptlbls
2417        print >>pFile,ptstr
2418        print >>pFile,varstr
2419       
2420    histDict = {}
2421    histVary = []
2422    controlDict = {}
2423    histoList = Histograms.keys()
2424    histoList.sort()
2425    for histogram in histoList:
2426        if 'PWDR' in histogram:
2427            Histogram = Histograms[histogram]
2428            hId = Histogram['hId']
2429            pfx = ':'+str(hId)+':'
2430            controlDict[pfx+'wtFactor'] = Histogram['wtFactor']
2431            controlDict[pfx+'Limits'] = Histogram['Limits'][1]
2432           
2433            Background = Histogram['Background']
2434            Type,bakDict,bakVary = GetBackgroundParms(hId,Background)
2435            controlDict[pfx+'bakType'] = Type
2436            histDict.update(bakDict)
2437            histVary += bakVary
2438           
2439            Inst = Histogram['Instrument Parameters'][0]
2440            Type,instDict,insVary = GetInstParms(hId,Inst)
2441            controlDict[pfx+'histType'] = Type
2442            if pfx+'Lam1' in instDict:
2443                controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam1']
2444            else:
2445                controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam']           
2446            histDict.update(instDict)
2447            histVary += insVary
2448           
2449            Sample = Histogram['Sample Parameters']
2450            Type,sampDict,sampVary = GetSampleParms(hId,Sample)
2451            controlDict[pfx+'instType'] = Type
2452            histDict.update(sampDict)
2453            histVary += sampVary
2454           
2455   
2456            if Print: 
2457                print >>pFile,'\n Histogram: ',histogram,' histogram Id: ',hId
2458                print >>pFile,135*'-'
2459                Units = {'C':' deg','T':' msec'}
2460                units = Units[controlDict[pfx+'histType'][2]]
2461                Limits = controlDict[pfx+'Limits']
2462                print >>pFile,' Instrument type: ',Sample['Type']
2463                print >>pFile,' Histogram limits: %8.2f%s to %8.2f%s'%(Limits[0],units,Limits[1],units)     
2464                PrintSampleParms(Sample)
2465                PrintInstParms(Inst)
2466                PrintBackground(Background)
2467        elif 'HKLF' in histogram:
2468            Histogram = Histograms[histogram]
2469            hId = Histogram['hId']
2470            pfx = ':'+str(hId)+':'
2471            controlDict[pfx+'wtFactor'] =Histogram['wtFactor']
2472            Inst = Histogram['Instrument Parameters'][0]
2473            controlDict[pfx+'histType'] = Inst['Type'][0]
2474            histDict[pfx+'Lam'] = Inst['Lam'][1]
2475            controlDict[pfx+'keV'] = 12.397639/histDict[pfx+'Lam']                   
2476    return histVary,histDict,controlDict
2477   
2478def SetHistogramData(parmDict,sigDict,Histograms,Print=True,pFile=None):
2479   
2480    def SetBackgroundParms(pfx,Background,parmDict,sigDict):
2481        Back = Background[0]
2482        DebyePeaks = Background[1]
2483        lenBack = len(Back[3:])
2484        backSig = [0 for i in range(lenBack+3*DebyePeaks['nDebye']+4*DebyePeaks['nPeaks'])]
2485        for i in range(lenBack):
2486            Back[3+i] = parmDict[pfx+'Back:'+str(i)]
2487            if pfx+'Back:'+str(i) in sigDict:
2488                backSig[i] = sigDict[pfx+'Back:'+str(i)]
2489        if DebyePeaks['nDebye']:
2490            for i in range(DebyePeaks['nDebye']):
2491                names = [pfx+'DebyeA:'+str(i),pfx+'DebyeR:'+str(i),pfx+'DebyeU:'+str(i)]
2492                for j,name in enumerate(names):
2493                    DebyePeaks['debyeTerms'][i][2*j] = parmDict[name]
2494                    if name in sigDict:
2495                        backSig[lenBack+3*i+j] = sigDict[name]           
2496        if DebyePeaks['nPeaks']:
2497            for i in range(DebyePeaks['nPeaks']):
2498                names = [pfx+'BkPkpos:'+str(i),pfx+'BkPkint:'+str(i),
2499                    pfx+'BkPksig:'+str(i),pfx+'BkPkgam:'+str(i)]
2500                for j,name in enumerate(names):
2501                    DebyePeaks['peaksList'][i][2*j] = parmDict[name]
2502                    if name in sigDict:
2503                        backSig[lenBack+3*DebyePeaks['nDebye']+4*i+j] = sigDict[name]
2504        return backSig
2505       
2506    def SetInstParms(pfx,Inst,parmDict,sigDict):
2507        instSig = {}
2508        insKeys = Inst.keys()
2509        insKeys.sort()
2510        for item in insKeys:
2511            insName = pfx+item
2512            Inst[item][1] = parmDict[insName]
2513            if insName in sigDict:
2514                instSig[item] = sigDict[insName]
2515            else:
2516                instSig[item] = 0
2517        return instSig
2518       
2519    def SetSampleParms(pfx,Sample,parmDict,sigDict):
2520        if 'Bragg' in Sample['Type']:             #Bragg-Brentano
2521            sampSig = [0 for i in range(3)]
2522            for i,item in enumerate(['Scale','Shift','Transparency']):       #surface roughness?, diffuse scattering?
2523                Sample[item][0] = parmDict[pfx+item]
2524                if pfx+item in sigDict:
2525                    sampSig[i] = sigDict[pfx+item]
2526        elif 'Debye' in Sample['Type']:        #Debye-Scherrer
2527            sampSig = [0 for i in range(4)]
2528            for i,item in enumerate(['Scale','Absorption','DisplaceX','DisplaceY']):
2529                Sample[item][0] = parmDict[pfx+item]
2530                if pfx+item in sigDict:
2531                    sampSig[i] = sigDict[pfx+item]
2532        return sampSig
2533       
2534    def PrintBackgroundSig(Background,backSig):
2535        Back = Background[0]
2536        DebyePeaks = Background[1]
2537        lenBack = len(Back[3:])
2538        valstr = ' value : '
2539        sigstr = ' sig   : '
2540        refine = False
2541        for i,back in enumerate(Back[3:]):
2542            valstr += '%10.4g'%(back)
2543            if Back[1]:
2544                refine = True
2545                sigstr += '%10.4g'%(backSig[i])
2546            else:
2547                sigstr += 10*' '
2548        if refine:
2549            print >>pFile,'\n Background function: ',Back[0]
2550            print >>pFile,valstr
2551            print >>pFile,sigstr
2552        if DebyePeaks['nDebye']:
2553            ifAny = False
2554            ptfmt = "%12.3f"
2555            names =  ' names :'
2556            ptstr =  ' values:'
2557            sigstr = ' esds  :'
2558            for item in sigDict:
2559                if 'Debye' in item:
2560                    ifAny = True
2561                    names += '%12s'%(item)
2562                    ptstr += ptfmt%(parmDict[item])
2563                    sigstr += ptfmt%(sigDict[item])
2564            if ifAny:
2565                print >>pFile,'\n Debye diffuse scattering coefficients'
2566                print >>pFile,names
2567                print >>pFile,ptstr
2568                print >>pFile,sigstr
2569        if DebyePeaks['nPeaks']:
2570            ifAny = False
2571            ptfmt = "%14.3f"
2572            names =  ' names :'
2573            ptstr =  ' values:'
2574            sigstr = ' esds  :'
2575            for item in sigDict:
2576                if 'BkPk' in item:
2577                    ifAny = True
2578                    names += '%14s'%(item)
2579                    ptstr += ptfmt%(parmDict[item])
2580                    sigstr += ptfmt%(sigDict[item])
2581            if ifAny:
2582                print >>pFile,'\n Single peak coefficients'
2583                print >>pFile,names
2584                print >>pFile,ptstr
2585                print >>pFile,sigstr
2586       
2587    def PrintInstParmsSig(Inst,instSig):
2588        ptlbls = ' names :'
2589        ptstr =  ' value :'
2590        sigstr = ' sig   :'
2591        refine = False
2592        insKeys = instSig.keys()
2593        insKeys.sort()
2594        for name in insKeys:
2595            if name not in  ['Type','Lam1','Lam2','Azimuth']:
2596                ptlbls += '%12s' % (name)
2597                ptstr += '%12.6f' % (Inst[name][1])
2598                if instSig[name]:
2599                    refine = True
2600                    sigstr += '%12.6f' % (instSig[name])
2601                else:
2602                    sigstr += 12*' '
2603        if refine:
2604            print >>pFile,'\n Instrument Parameters:'
2605            print >>pFile,ptlbls
2606            print >>pFile,ptstr
2607            print >>pFile,sigstr
2608       
2609    def PrintSampleParmsSig(Sample,sampleSig):
2610        ptlbls = ' names :'
2611        ptstr =  ' values:'
2612        sigstr = ' sig   :'
2613        refine = False
2614        if 'Bragg' in Sample['Type']:
2615            for i,item in enumerate(['Scale','Shift','Transparency']):
2616                ptlbls += '%14s'%(item)
2617                ptstr += '%14.4f'%(Sample[item][0])
2618                if sampleSig[i]:
2619                    refine = True
2620                    sigstr += '%14.4f'%(sampleSig[i])
2621                else:
2622                    sigstr += 14*' '
2623           
2624        elif 'Debye' in Sample['Type']:        #Debye-Scherrer
2625            for i,item in enumerate(['Scale','Absorption','DisplaceX','DisplaceY']):
2626                ptlbls += '%14s'%(item)
2627                ptstr += '%14.4f'%(Sample[item][0])
2628                if sampleSig[i]:
2629                    refine = True
2630                    sigstr += '%14.4f'%(sampleSig[i])
2631                else:
2632                    sigstr += 14*' '
2633
2634        if refine:
2635            print >>pFile,'\n Sample Parameters:'
2636            print >>pFile,ptlbls
2637            print >>pFile,ptstr
2638            print >>pFile,sigstr
2639       
2640    histoList = Histograms.keys()
2641    histoList.sort()
2642    for histogram in histoList:
2643        if 'PWDR' in histogram:
2644            Histogram = Histograms[histogram]
2645            hId = Histogram['hId']
2646            pfx = ':'+str(hId)+':'
2647            Background = Histogram['Background']
2648            backSig = SetBackgroundParms(pfx,Background,parmDict,sigDict)
2649           
2650            Inst = Histogram['Instrument Parameters'][0]
2651            instSig = SetInstParms(pfx,Inst,parmDict,sigDict)
2652       
2653            Sample = Histogram['Sample Parameters']
2654            sampSig = SetSampleParms(pfx,Sample,parmDict,sigDict)
2655
2656            print >>pFile,'\n Histogram: ',histogram,' histogram Id: ',hId
2657            print >>pFile,135*'-'
2658            print >>pFile,' Final refinement wR = %.2f%% on %d observations in this histogram'%(Histogram['wR'],Histogram['Nobs'])
2659            print >>pFile,' PWDR histogram weight factor = '+'%.3f'%(Histogram['wtFactor'])
2660            if Print:
2661                print >>pFile,' Instrument type: ',Sample['Type']
2662                PrintSampleParmsSig(Sample,sampSig)
2663                PrintInstParmsSig(Inst,instSig)
2664                PrintBackgroundSig(Background,backSig)
2665               
2666################################################################################
2667##### Penalty & restraint functions
2668################################################################################
2669
2670def penaltyFxn(HistoPhases,parmDict,varyList):
2671    Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
2672    pNames = []
2673    pVals = []
2674    pWt = []
2675    negWt = {}
2676    for phase in Phases:
2677        pId = Phases[phase]['pId']
2678        negWt[pId] = Phases[phase]['General']['Pawley neg wt']
2679        General = Phases[phase]['General']
2680        textureData = General['SH Texture']
2681        SGData = General['SGData']
2682        AtLookup = G2mth.FillAtomLookUp(Phases[phase]['Atoms'])
2683        cell = General['Cell'][1:7]
2684        Amat,Bmat = G2lat.cell2AB(cell)
2685        if phase not in restraintDict:
2686            continue
2687        phaseRest = restraintDict[phase]
2688        names = [['Bond','Bonds'],['Angle','Angles'],['Plane','Planes'],
2689            ['Chiral','Volumes'],['Torsion','Torsions'],['Rama','Ramas'],
2690            ['ChemComp','Sites'],['Texture','HKLs']]
2691        for name,rest in names:
2692            itemRest = phaseRest[name]
2693            if itemRest[rest] and itemRest['Use']:
2694                wt = itemRest['wtFactor']
2695                if name in ['Bond','Angle','Plane','Chiral']:
2696                    for i,[indx,ops,obs,esd] in enumerate(itemRest[rest]):
2697                        pNames.append(str(pId)+':'+name+':'+str(i))
2698                        XYZ = np.array(G2mth.GetAtomCoordsByID(pId,parmDict,AtLookup,indx))
2699                        XYZ = G2mth.getSyXYZ(XYZ,ops,SGData)
2700                        if name == 'Bond':
2701                            calc = G2mth.getRestDist(XYZ,Amat)
2702                        elif name == 'Angle':
2703                            calc = G2mth.getRestAngle(XYZ,Amat)
2704                        elif name == 'Plane':
2705                            calc = G2mth.getRestPlane(XYZ,Amat)
2706                        elif name == 'Chiral':
2707                            calc = G2mth.getRestChiral(XYZ,Amat)
2708                        pVals.append(obs-calc)
2709                        pWt.append(wt/esd**2)
2710                elif name in ['Torsion','Rama']:
2711                    coeffDict = itemRest['Coeff']
2712                    for i,[indx,ops,cofName,esd] in enumerate(itemRest[rest]):
2713                        pNames.append(str(pId)+':'+name+':'+str(i))
2714                        XYZ = np.array(G2mth.GetAtomCoordsByID(pId,parmDict,AtLookup,indx))
2715                        XYZ = G2mth.getSyXYZ(XYZ,ops,SGData)
2716                        if name == 'Torsion':
2717                            tor = G2mth.getRestTorsion(XYZ,Amat)
2718                            restr,calc = G2mth.calcTorsionEnergy(tor,coeffDict[cofName])
2719                        else:
2720                            phi,psi = G2mth.getRestRama(XYZ,Amat)
2721                            restr,calc = G2mth.calcRamaEnergy(phi,psi,coeffDict[cofName])                               
2722                        pVals.append(obs-calc)
2723                        pWt.append(wt/esd**2)
2724                elif name == 'ChemComp':
2725                    for i,[indx,factors,obs,esd] in enumerate(itemRest[rest]):
2726                        pNames.append(str(pId)+':'+name+':'+str(i))
2727                        mul = np.array(G2mth.GetAtomItemsById(Atoms,AtLookUp,indx,cs+1))
2728                        frac = np.array(G2mth.GetAtomItemsById(Atoms,AtLookUp,indx,cs-1))
2729                        calc = np.sum(mul*frac*factors)
2730                        pVals.append(obs-calc)
2731                        pWt.append(wt/esd**2)                   
2732                elif name == 'Texture':
2733                    SHkeys = textureData['SH Coeff'][1].keys()
2734                    SHCoef = G2mth.GetSHCoeff(pId,parmDict,SHkeys)
2735                    shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
2736                    SamSym = dict(zip(shModels,['0','-1','2/m','mmm']))
2737                    for i,[hkl,grid,esd1,ifesd2,esd2] in enumerate(itemRest[rest]):
2738                        PH = np.array(hkl)
2739                        phi,beta = G2lat.CrsAng(np.array(hkl),cell,SGData)
2740                        ODFln = G2lat.Flnh(False,SHCoef,phi,beta,SGData)
2741                        R,P,Z = G2mth.getRestPolefig(ODFln,SamSym[textureData['Model']],grid)
2742                        Z1 = -ma.masked_greater(Z,0.0)
2743                        IndZ1 = np.array(ma.nonzero(Z1))
2744                        for ind in IndZ1.T:
2745                            pNames.append('%d:%s:%d:%.2f:%.2f'%(pId,name,i,R[ind[0],ind[1]],P[ind[0],ind[1]]))
2746                            pVals.append(Z1[ind[0]][ind[1]])
2747                            pWt.append(wt/esd1**2)
2748                        if ifesd2:
2749                            Z2 = 1.-Z
2750                            for ind in np.ndindex(grid,grid):
2751                                pNames.append('%d:%s:%d:%.2f:%.2f'%(pId,name+'-unit',i,R[ind[0],ind[1]],P[ind[0],ind[1]]))
2752                                pVals.append(Z1[ind[0]][ind[1]])
2753                                pWt.append(wt/esd2**2)
2754         
2755    for item in varyList:
2756        if 'PWLref' in item and parmDict[item] < 0.:
2757            pId = int(item.split(':')[0])
2758            if negWt[pId]:
2759                pNames.append(item)
2760                pVals.append(-parmDict[item])
2761                pWt.append(negWt[pId])
2762    pVals = np.array(pVals)
2763    pWt = np.array(pWt)         #should this be np.sqrt?
2764    return pNames,pVals,pWt
2765   
2766def penaltyDeriv(pNames,pVal,HistoPhases,parmDict,varyList):
2767    Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
2768    pDerv = np.zeros((len(varyList),len(pVal)))
2769    for phase in Phases:
2770        if phase not in restraintDict:
2771            continue
2772        pId = Phases[phase]['pId']
2773        General = Phases[phase]['General']
2774        SGData = General['SGData']
2775        AtLookup = G2mth.FillAtomLookUp(Phases[phase]['Atoms'])
2776        cell = General['Cell'][1:7]
2777        Amat,Bmat = G2lat.cell2AB(cell)
2778        textureData = General['SH Texture']
2779
2780        SHkeys = textureData['SH Coeff'][1].keys()
2781        SHCoef = G2mth.GetSHCoeff(pId,parmDict,SHkeys)
2782        shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
2783        SamSym = dict(zip(shModels,['0','-1','2/m','mmm']))
2784        sam = SamSym[textureData['Model']]
2785        phaseRest = restraintDict[phase]
2786        names = {'Bond':'Bonds','Angle':'Angles','Plane':'Planes',
2787            'Chiral':'Volumes','Torsion':'Torsions','Rama':'Ramas',
2788            'ChemComp':'Sites','Texture':'HKLs'}
2789        lasthkl = np.array([0,0,0])
2790        for ip,pName in enumerate(pNames):
2791            pnames = pName.split(':')
2792            if pId == int(pnames[0]):
2793                name = pnames[1]
2794                if not name:        #empty for Pawley restraints; pName has '::' in it
2795                    pDerv[varyList.index(pName)][ip] += 1.
2796                    continue
2797                id = int(pnames[2]) 
2798                itemRest = phaseRest[name]
2799                if name in ['Bond','Angle','Plane','Chiral']:
2800                    indx,ops,obs,esd = itemRest[names[name]][id]
2801                    dNames = []
2802                    for ind in indx:
2803                        dNames += [str(pId)+'::dA'+Xname+':'+str(AtLookup[ind]) for Xname in ['x','y','z']]
2804                    XYZ = np.array(G2mth.GetAtomCoordsByID(pId,parmDict,AtLookup,indx))
2805                    if name == 'Bond':
2806                        deriv = G2mth.getRestDeriv(G2mth.getRestDist,XYZ,Amat,ops,SGData)
2807                    elif name == 'Angle':
2808                        deriv = G2mth.getRestDeriv(G2mth.getRestAngle,XYZ,Amat,ops,SGData)
2809                    elif name == 'Plane':
2810                        deriv = G2mth.getRestDeriv(G2mth.getRestPlane,XYZ,Amat,ops,SGData)
2811                    elif name == 'Chiral':
2812                        deriv = G2mth.getRestDeriv(G2mth.getRestChiral,XYZ,Amat,ops,SGData)
2813                elif name in ['Torsion','Rama']:
2814                    coffDict = itemRest['Coeff']
2815                    indx,ops,cofName,esd = itemRest[names[name]][id]
2816                    dNames = []
2817                    for ind in indx:
2818                        dNames += [str(pId)+'::dA'+Xname+':'+str(AtLookup[ind]) for Xname in ['x','y','z']]
2819                    XYZ = np.array(G2mth.GetAtomCoordsByID(pId,parmDict,AtLookup,indx))
2820                    if name == 'Torsion':
2821                        deriv = G2mth.getTorsionDeriv(XYZ,Amat,coffDict[cofName])
2822                    else:
2823                        deriv = G2mth.getRamaDeriv(XYZ,Amat,coffDict[cofName])
2824                elif name == 'ChemComp':
2825                    indx,factors,obs,esd = itemRest[names[name]][id]
2826                    dNames = []
2827                    for ind in indx:
2828                        dNames += [str(pId)+'::Afrac:'+str(AtLookup[ind])]
2829                        mul = np.array(G2mth.GetAtomItemsById(Atoms,AtLookUp,indx,cs+1))
2830                        deriv = mul*factors
2831                elif 'Texture' in name:
2832                    deriv = []
2833                    dNames = []
2834                    hkl,grid,esd1,ifesd2,esd2 = itemRest[names[name]][id]
2835                    hkl = np.array(hkl)
2836                    if np.any(lasthkl-hkl):
2837                        PH = np.array(hkl)
2838                        phi,beta = G2lat.CrsAng(np.array(hkl),cell,SGData)
2839                        ODFln = G2lat.Flnh(False,SHCoef,phi,beta,SGData)
2840                        lasthkl = copy.copy(hkl)                       
2841                    if 'unit' in name:
2842                        pass
2843                    else:
2844                        gam = float(pnames[3])
2845                        psi = float(pnames[4])
2846                        for SHname in ODFln:
2847                            l,m,n = eval(SHname[1:])
2848                            Ksl = G2lat.GetKsl(l,m,sam,psi,gam)[0]
2849                            dNames += [str(pId)+'::'+SHname]
2850                            deriv.append(-ODFln[SHname][0]*Ksl/SHCoef[SHname])
2851                for dName,drv in zip(dNames,deriv):
2852                    try:
2853                        ind = varyList.index(dName)
2854                        pDerv[ind][ip] += drv
2855                    except ValueError:
2856                        pass
2857    return pDerv
2858
2859################################################################################
2860##### Function & derivative calculations
2861################################################################################       
2862                   
2863def GetAtomFXU(pfx,calcControls,parmDict):
2864    Natoms = calcControls['Natoms'][pfx]
2865    Tdata = Natoms*[' ',]
2866    Mdata = np.zeros(Natoms)
2867    IAdata = Natoms*[' ',]
2868    Fdata = np.zeros(Natoms)
2869    FFdata = []
2870    BLdata = []
2871    Xdata = np.zeros((3,Natoms))
2872    dXdata = np.zeros((3,Natoms))
2873    Uisodata = np.zeros(Natoms)
2874    Uijdata = np.zeros((6,Natoms))
2875    keys = {'Atype:':Tdata,'Amul:':Mdata,'Afrac:':Fdata,'AI/A:':IAdata,
2876        'dAx:':dXdata[0],'dAy:':dXdata[1],'dAz:':dXdata[2],
2877        'Ax:':Xdata[0],'Ay:':Xdata[1],'Az:':Xdata[2],'AUiso:':Uisodata,
2878        'AU11:':Uijdata[0],'AU22:':Uijdata[1],'AU33:':Uijdata[2],
2879        'AU12:':Uijdata[3],'AU13:':Uijdata[4],'AU23:':Uijdata[5]}
2880    for iatm in range(Natoms):
2881        for key in keys:
2882            parm = pfx+key+str(iatm)
2883            if parm in parmDict:
2884                keys[key][iatm] = parmDict[parm]
2885    return Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata
2886   
2887def getFFvalues(FFtables,SQ):
2888    FFvals = {}
2889    for El in FFtables:
2890        FFvals[El] = G2el.ScatFac(FFtables[El],SQ)[0]
2891    return FFvals
2892   
2893def getBLvalues(BLtables):
2894    BLvals = {}
2895    for El in BLtables:
2896        BLvals[El] = BLtables[El][1][1]
2897    return BLvals
2898       
2899def StructureFactor(refList,G,hfx,pfx,SGData,calcControls,parmDict):
2900    ''' Compute structure factors for all h,k,l for phase
2901    input:
2902        refList: [ref] where each ref = h,k,l,m,d,...,[equiv h,k,l],phase[equiv]
2903        G:      reciprocal metric tensor
2904        pfx:    phase id string
2905        SGData: space group info. dictionary output from SpcGroup
2906        calcControls:
2907        ParmDict:
2908    puts result F^2 in each ref[8] in refList
2909    '''       
2910    twopi = 2.0*np.pi
2911    twopisq = 2.0*np.pi**2
2912    phfx = pfx.split(':')[0]+hfx
2913    ast = np.sqrt(np.diag(G))
2914    Mast = twopisq*np.multiply.outer(ast,ast)
2915    FFtables = calcControls['FFtables']
2916    BLtables = calcControls['BLtables']
2917    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
2918    FF = np.zeros(len(Tdata))
2919    if 'N' in parmDict[hfx+'Type']:
2920        FP,FPP = G2el.BlenRes(Tdata,BLtables,parmDict[hfx+'Lam'])
2921    else:
2922        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
2923        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
2924    maxPos = len(SGData['SGOps'])
2925    Uij = np.array(G2lat.U6toUij(Uijdata))
2926    bij = Mast*Uij.T
2927    for refl in refList:
2928        fbs = np.array([0,0])
2929        H = refl[:3]
2930        SQ = 1./(2.*refl[4])**2
2931        SQfactor = 4.0*SQ*twopisq
2932        Bab = parmDict[phfx+'BabA']*np.exp(-parmDict[phfx+'BabU']*SQfactor)
2933        if not len(refl[-1]):                #no form factors
2934            if 'N' in parmDict[hfx+'Type']:
2935                refl[-1] = getBLvalues(BLtables)
2936            else:       #'X'
2937                refl[-1] = getFFvalues(FFtables,SQ)
2938        for i,El in enumerate(Tdata):
2939            FF[i] = refl[-1][El]           
2940        Uniq = refl[11]
2941        phi = refl[12]
2942        phase = twopi*(np.inner(Uniq,(dXdata.T+Xdata.T))+phi[:,np.newaxis])
2943        sinp = np.sin(phase)
2944        cosp = np.cos(phase)
2945        occ = Mdata*Fdata/len(Uniq)
2946        biso = -SQfactor*Uisodata
2947        Tiso = np.where(biso<1.,np.exp(biso),1.0)
2948        HbH = np.array([-np.inner(h,np.inner(bij,h)) for h in Uniq])
2949        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
2950        Tcorr = Tiso*Tuij
2951        fa = np.array([(FF+FP-Bab)*occ*cosp*Tcorr,-FPP*occ*sinp*Tcorr])
2952        fas = np.sum(np.sum(fa,axis=1),axis=1)        #real
2953        if not SGData['SGInv']:
2954            fb = np.array([(FF+FP-Bab)*occ*sinp*Tcorr,FPP*occ*cosp*Tcorr])
2955            fbs = np.sum(np.sum(fb,axis=1),axis=1)
2956        fasq = fas**2
2957        fbsq = fbs**2        #imaginary
2958        refl[9] = np.sum(fasq)+np.sum(fbsq)
2959        refl[10] = atan2d(fbs[0],fas[0])
2960    return refList
2961   
2962def StructureFactorDerv(refList,G,hfx,pfx,SGData,calcControls,parmDict):
2963    twopi = 2.0*np.pi
2964    twopisq = 2.0*np.pi**2
2965    phfx = pfx.split(':')[0]+hfx
2966    ast = np.sqrt(np.diag(G))
2967    Mast = twopisq*np.multiply.outer(ast,ast)
2968    FFtables = calcControls['FFtables']
2969    BLtables = calcControls['BLtables']
2970    Tdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,calcControls,parmDict)
2971    FF = np.zeros(len(Tdata))
2972    if 'N' in parmDict[hfx+'Type']:
2973        FP = 0.
2974        FPP = 0.
2975    else:
2976        FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
2977        FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
2978    maxPos = len(SGData['SGOps'])       
2979    Uij = np.array(G2lat.U6toUij(Uijdata))
2980    bij = Mast*Uij.T
2981    dFdvDict = {}
2982    dFdfr = np.zeros((len(refList),len(Mdata)))
2983    dFdx = np.zeros((len(refList),len(Mdata),3))
2984    dFdui = np.zeros((len(refList),len(Mdata)))
2985    dFdua = np.zeros((len(refList),len(Mdata),6))
2986    dFdbab = np.zeros((len(refList),2))
2987    for iref,refl in enumerate(refList):
2988        H = np.array(refl[:3])
2989        SQ = 1./(2.*refl[4])**2             # or (sin(theta)/lambda)**2
2990        SQfactor = 8.0*SQ*np.pi**2
2991        dBabdA = np.exp(-parmDict[phfx+'BabU']*SQfactor)
2992        Bab = parmDict[phfx+'BabA']*dBabdA
2993        for i,El in enumerate(Tdata):           
2994            FF[i] = refl[-1][El]           
2995        Uniq = refl[11]
2996        phi = refl[12]
2997        phase = twopi*(np.inner((dXdata.T+Xdata.T),Uniq)+phi[np.newaxis,:])
2998        sinp = np.sin(phase)
2999        cosp = np.cos(phase)
3000        occ = Mdata*Fdata/len(Uniq)
3001        biso = -SQfactor*Uisodata
3002        Tiso = np.where(biso<1.,np.exp(biso),1.0)
3003        HbH = -np.inner(H,np.inner(bij,H))
3004        Hij = np.array([Mast*np.multiply.outer(U,U) for U in Uniq])
3005        Hij = np.array([G2lat.UijtoU6(Uij) for Uij in Hij])
3006        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
3007        Tcorr = Tiso*Tuij
3008        fot = (FF+FP-Bab)*occ*Tcorr
3009        fotp = FPP*occ*Tcorr
3010        fa = np.array([fot[:,np.newaxis]*cosp,fotp[:,np.newaxis]*cosp])       #non positions
3011        fb = np.array([fot[:,np.newaxis]*sinp,-fotp[:,np.newaxis]*sinp])
3012       
3013        fas = np.sum(np.sum(fa,axis=1),axis=1)
3014        fbs = np.sum(np.sum(fb,axis=1),axis=1)
3015        fax = np.array([-fot[:,np.newaxis]*sinp,-fotp[:,np.newaxis]*sinp])   #positions
3016        fbx = np.array([fot[:,np.newaxis]*cosp,-fot[:,np.newaxis]*cosp])
3017        #sum below is over Uniq
3018        dfadfr = np.sum(fa/occ[:,np.newaxis],axis=2)
3019        dfadx = np.sum(twopi*Uniq*fax[:,:,:,np.newaxis],axis=2)
3020        dfadui = np.sum(-SQfactor*fa,axis=2)
3021        dfadua = np.sum(-Hij*fa[:,:,:,np.newaxis],axis=2)
3022        dfadba = np.sum(-cosp*(occ*Tcorr)[:,np.newaxis],axis=1)
3023        #NB: the above have been checked against PA(1:10,1:2) in strfctr.for     
3024        dFdfr[iref] = 2.*(fas[0]*dfadfr[0]+fas[1]*dfadfr[1])*Mdata/len(Uniq)
3025        dFdx[iref] = 2.*(fas[0]*dfadx[0]+fas[1]*dfadx[1])
3026        dFdui[iref] = 2.*(fas[0]*dfadui[0]+fas[1]*dfadui[1])
3027        dFdua[iref] = 2.*(fas[0]*dfadua[0]+fas[1]*dfadua[1])
3028        dFdbab[iref] = np.array([np.sum(dfadba*dBabdA),np.sum(-dfadba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T
3029        if not SGData['SGInv']:
3030            dfbdfr = np.sum(fb/occ[:,np.newaxis],axis=2)        #problem here if occ=0 for some atom
3031            dfbdx = np.sum(twopi*Uniq*fbx[:,:,:,np.newaxis],axis=2)         
3032            dfbdui = np.sum(-SQfactor*fb,axis=2)
3033            dfbdua = np.sum(-Hij*fb[:,:,:,np.newaxis],axis=2)
3034            dfbdba = np.sum(-sinp*(occ*Tcorr)[:,np.newaxis],axis=1)
3035            dFdfr[iref] += 2.*(fbs[0]*dfbdfr[0]-fbs[1]*dfbdfr[1])*Mdata/len(Uniq)
3036            dFdx[iref] += 2.*(fbs[0]*dfbdx[0]+fbs[1]*dfbdx[1])
3037            dFdui[iref] += 2.*(fbs[0]*dfbdui[0]-fbs[1]*dfbdui[1])
3038            dFdua[iref] += 2.*(fbs[0]*dfbdua[0]+fbs[1]*dfbdua[1])
3039            dFdbab[iref] += np.array([np.sum(dfbdba*dBabdA),np.sum(-dfbdba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T
3040        #loop over atoms - each dict entry is list of derivatives for all the reflections
3041    for i in range(len(Mdata)):     
3042        dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i]
3043        dFdvDict[pfx+'dAx:'+str(i)] = dFdx.T[0][i]
3044        dFdvDict[pfx+'dAy:'+str(i)] = dFdx.T[1][i]
3045        dFdvDict[pfx+'dAz:'+str(i)] = dFdx.T[2][i]
3046        dFdvDict[pfx+'AUiso:'+str(i)] = dFdui.T[i]
3047        dFdvDict[pfx+'AU11:'+str(i)] = dFdua.T[0][i]
3048        dFdvDict[pfx+'AU22:'+str(i)] = dFdua.T[1][i]
3049        dFdvDict[pfx+'AU33:'+str(i)] = dFdua.T[2][i]
3050        dFdvDict[pfx+'AU12:'+str(i)] = 2.*dFdua.T[3][i]
3051        dFdvDict[pfx+'AU13:'+str(i)] = 2.*dFdua.T[4][i]
3052        dFdvDict[pfx+'AU23:'+str(i)] = 2.*dFdua.T[5][i]
3053        dFdvDict[pfx+'BabA'] = dFdbab.T[0]
3054        dFdvDict[pfx+'BabU'] = dFdbab.T[1]
3055    return dFdvDict
3056   
3057def SCExtinction(ref,phfx,hfx,pfx,calcControls,parmDict,varyList):
3058    ''' Single crystal extinction function; puts correction in ref[13] and returns
3059    corrections needed for derivatives
3060    '''
3061    ref[13] = 1.0
3062    dervCor = 1.0
3063    dervDict = {}
3064    if calcControls[phfx+'EType'] != 'None':
3065        cos2T = 1.0-0.5*(parmDict[hfx+'Lam']/ref[4])**2         #cos(2theta)
3066        if 'SXC' in parmDict[hfx+'Type']:
3067            AV = 7.9406e5/parmDict[pfx+'Vol']**2
3068            PL = np.sqrt(1.0-cos2T**2)/parmDict[hfx+'Lam']
3069            P12 = (calcControls[phfx+'Cos2TM']+cos2T**4)/(calcControls[phfx+'Cos2TM']+cos2T**2)
3070        elif 'SNT' in parmDict[hfx+'Type']:
3071            AV = 1.e7/parmDict[pfx+'Vol']**2
3072            PL = 1./(4.*refl[4]**2)
3073            P12 = 1.0
3074        elif 'SNC' in parmDict[hfx+'Type']:
3075            AV = 1.e7/parmDict[pfx+'Vol']**2
3076            PL = np.sqrt(1.0-cos2T**2)/parmDict[hfx+'Lam']
3077            P12 = 1.0
3078           
3079        PLZ = AV*P12*parmDict[hfx+'Lam']**2*ref[7]
3080        if 'Primary' in calcControls[phfx+'EType']:
3081            PLZ *= 1.5
3082        else:
3083            PLZ *= calcControls[phfx+'Tbar']
3084                       
3085        if 'Primary' in calcControls[phfx+'EType']:
3086            PSIG = parmDict[phfx+'Ep']
3087        elif 'I & II' in calcControls[phfx+'EType']:
3088            PSIG = parmDict[phfx+'Eg']/np.sqrt(1.+(parmDict[phfx+'Es']*PL/parmDict[phfx+'Eg'])**2)
3089        elif 'Type II' in calcControls[phfx+'EType']:
3090            PSIG = parmDict[phfx+'Es']
3091        else:       # 'Secondary Type I'
3092            PSIG = parmDict[phfx+'Eg']/PL
3093           
3094        AG = 0.58+0.48*cos2T+0.24*cos2T**2
3095        AL = 0.025+0.285*cos2T
3096        BG = 0.02-0.025*cos2T
3097        BL = 0.15-0.2*(0.75-cos2T)**2
3098        if cos2T < 0.:
3099            BL = -0.45*cos2T
3100        CG = 2.
3101        CL = 2.
3102        PF = PLZ*PSIG
3103       
3104        if 'Gaussian' in calcControls[phfx+'EApprox']:
3105            PF4 = 1.+CG*PF+AG*PF**2/(1.+BG*PF)
3106            extCor = np.sqrt(PF4)
3107            PF3 = 0.5*(CG+2.*AG*PF/(1.+BG*PF)-AG*PF**2*BG/(1.+BG*PF)**2)/(PF4*extCor)
3108        else:
3109            PF4 = 1.+CL*PF+AL*PF**2/(1.+BL*PF)
3110            extCor = np.sqrt(PF4)
3111            PF3 = 0.5*(CL+2.*AL*PF/(1.+BL*PF)-AL*PF**2*BL/(1.+BL*PF)**2)/(PF4*extCor)
3112
3113        dervCor = (1.+PF)*PF3
3114        if 'Primary' in calcControls[phfx+'EType'] and phfx+'Ep' in varyList:
3115            dervDict[phfx+'Ep'] = -ref[7]*PLZ*PF3
3116        if 'II' in calcControls[phfx+'EType'] and phfx+'Es' in varyList:
3117            dervDict[phfx+'Es'] = -ref[7]*PLZ*PF3*(PSIG/parmDict[phfx+'Es'])**3
3118        if 'I' in calcControls[phfx+'EType'] and phfx+'Eg' in varyList:
3119            dervDict[phfx+'Eg'] = -ref[7]*PLZ*PF3*(PSIG/parmDict[phfx+'Eg'])**3*PL**2
3120               
3121        ref[13] = 1./extCor
3122    return dervCor,dervDict
3123       
3124   
3125def Dict2Values(parmdict, varylist):
3126    '''Use before call to leastsq to setup list of values for the parameters
3127    in parmdict, as selected by key in varylist'''
3128    return [parmdict[key] for key in varylist] 
3129   
3130def Values2Dict(parmdict, varylist, values):
3131    ''' Use after call to leastsq to update the parameter dictionary with
3132    values corresponding to keys in varylist'''
3133    parmdict.update(zip(varylist,values))
3134   
3135def GetNewCellParms(parmDict,varyList):
3136    newCellDict = {}
3137    Anames = ['A'+str(i) for i in range(6)]
3138    Ddict = dict(zip(['D11','D22','D33','D12','D13','D23'],Anames))
3139    for item in varyList:
3140        keys = item.split(':')
3141        if keys[2] in Ddict:
3142            key = keys[0]+'::'+Ddict[keys[2]]       #key is e.g. '0::A0'
3143            parm = keys[0]+'::'+keys[2]             #parm is e.g. '0::D11'
3144            newCellDict[parm] = [key,parmDict[key]+parmDict[item]]
3145    return newCellDict          # is e.g. {'0::D11':A0+D11}
3146   
3147def ApplyXYZshifts(parmDict,varyList):
3148    ''' takes atom x,y,z shift and applies it to corresponding atom x,y,z value
3149        input:
3150            parmDict - parameter dictionary
3151            varyList - list of variables
3152        returns:
3153            newAtomDict - dictionary of new atomic coordinate names & values;
3154                key is parameter shift name
3155    '''
3156    newAtomDict = {}
3157    for item in parmDict:
3158        if 'dA' in item:
3159            parm = ''.join(item.split('d'))
3160            parmDict[parm] += parmDict[item]
3161            newAtomDict[item] = [parm,parmDict[parm]]
3162    return newAtomDict
3163   
3164def SHTXcal(refl,g,pfx,hfx,SGData,calcControls,parmDict):
3165    #Spherical harmonics texture
3166    IFCoup = 'Bragg' in calcControls[hfx+'instType']
3167    odfCor = 1.0
3168    H = refl[:3]
3169    cell = G2lat.Gmat2cell(g)
3170    Sangls = [parmDict[pfx+'SH omega'],parmDict[pfx+'SH chi'],parmDict[pfx+'SH phi']]
3171    Gangls = [parmDict[hfx+'Omega'],parmDict[hfx+'Chi'],parmDict[hfx+'Phi'],parmDict[hfx+'Azimuth']]
3172    phi,beta = G2lat.CrsAng(H,cell,SGData)
3173    psi,gam,x,x = G2lat.SamAng(refl[5]/2.,Gangls,Sangls,IFCoup) #ignore 2 sets of angle derivs.
3174    SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],parmDict[pfx+'SHmodel'],parmDict[pfx+'SHorder'])
3175    for item in SHnames:
3176        L,M,N = eval(item.strip('C'))
3177        Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
3178        Ksl,x,x = G2lat.GetKsl(L,M,parmDict[pfx+'SHmodel'],psi,gam)
3179        Lnorm = G2lat.Lnorm(L)
3180        odfCor += parmDict[pfx+item]*Lnorm*Kcl*Ksl
3181    return odfCor
3182   
3183def SHTXcalDerv(refl,g,pfx,hfx,SGData,calcControls,parmDict):
3184    #Spherical harmonics texture derivatives
3185    FORPI = 4.0*np.pi
3186    IFCoup = 'Bragg' in calcControls[hfx+'instType']
3187    odfCor = 1.0
3188    dFdODF = {}
3189    dFdSA = [0,0,0]
3190    H = refl[:3]
3191    cell = G2lat.Gmat2cell(g)
3192    Sangls = [parmDict[pfx+'SH omega'],parmDict[pfx+'SH chi'],parmDict[pfx+'SH phi']]
3193    Gangls = [parmDict[hfx+'Omega'],parmDict[hfx+'Chi'],parmDict[hfx+'Phi'],parmDict[hfx+'Azimuth']]
3194    phi,beta = G2lat.CrsAng(H,cell,SGData)
3195    psi,gam,dPSdA,dGMdA = G2lat.SamAng(refl[5]/2.,Gangls,Sangls,IFCoup)
3196    SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],parmDict[pfx+'SHmodel'],parmDict[pfx+'SHorder'])
3197    for item in SHnames:
3198        L,M,N = eval(item.strip('C'))
3199        Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
3200        Ksl,dKsdp,dKsdg = G2lat.GetKsl(L,M,parmDict[pfx+'SHmodel'],psi,gam)
3201        Lnorm = G2lat.Lnorm(L)
3202        odfCor += parmDict[pfx+item]*Lnorm*Kcl*Ksl
3203        dFdODF[pfx+item] = Lnorm*Kcl*Ksl
3204        for i in range(3):
3205            dFdSA[i] += parmDict[pfx+item]*Lnorm*Kcl*(dKsdp*dPSdA[i]+dKsdg*dGMdA[i])
3206    return odfCor,dFdODF,dFdSA
3207   
3208def SHPOcal(refl,g,phfx,hfx,SGData,calcControls,parmDict):
3209    #sphericaql harmonics preferred orientation (cylindrical symmetry only)
3210    odfCor = 1.0
3211    H = refl[:3]
3212    cell = G2lat.Gmat2cell(g)
3213    Sangl = [0.,0.,0.]
3214    if 'Bragg' in calcControls[hfx+'instType']:
3215        Gangls = [0.,90.,0.,parmDict[hfx+'Azimuth']]
3216        IFCoup = True
3217    else:
3218        Gangls = [0.,0.,0.,parmDict[hfx+'Azimuth']]
3219        IFCoup = False
3220    phi,beta = G2lat.CrsAng(H,cell,SGData)
3221    psi,gam,x,x = G2lat.SamAng(refl[5]/2.,Gangls,Sangl,IFCoup) #ignore 2 sets of angle derivs.
3222    SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],'0',calcControls[phfx+'SHord'],False)
3223    for item in SHnames:
3224        L,N = eval(item.strip('C'))
3225        Kcsl,Lnorm = G2lat.GetKclKsl(L,N,SGData['SGLaue'],psi,phi,beta) 
3226        odfCor += parmDict[phfx+item]*Lnorm*Kcsl
3227    return odfCor
3228   
3229def SHPOcalDerv(refl,g,phfx,hfx,SGData,calcControls,parmDict):
3230    #spherical harmonics preferred orientation derivatives (cylindrical symmetry only)
3231    FORPI = 12.5663706143592
3232    odfCor = 1.0
3233    dFdODF = {}
3234    H = refl[:3]
3235    cell = G2lat.Gmat2cell(g)
3236    Sangl = [0.,0.,0.]
3237    if 'Bragg' in calcControls[hfx+'instType']:
3238        Gangls = [0.,90.,0.,parmDict[hfx+'Azimuth']]
3239        IFCoup = True
3240    else:
3241        Gangls = [0.,0.,0.,parmDict[hfx+'Azimuth']]
3242        IFCoup = False
3243    phi,beta = G2lat.CrsAng(H,cell,SGData)
3244    psi,gam,x,x = G2lat.SamAng(refl[5]/2.,Gangls,Sangl,IFCoup) #ignore 2 sets of angle derivs.
3245    SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],'0',calcControls[phfx+'SHord'],False)
3246    for item in SHnames:
3247        L,N = eval(item.strip('C'))
3248        Kcsl,Lnorm = G2lat.GetKclKsl(L,N,SGData['SGLaue'],psi,phi,beta) 
3249        odfCor += parmDict[phfx+item]*Lnorm*Kcsl
3250        dFdODF[phfx+item] = Kcsl*Lnorm
3251    return odfCor,dFdODF
3252   
3253def GetPrefOri(refl,G,g,phfx,hfx,SGData,calcControls,parmDict):
3254    POcorr = 1.0
3255    if calcControls[phfx+'poType'] == 'MD':
3256        MD = parmDict[phfx+'MD']
3257        if MD != 1.0:
3258            MDAxis = calcControls[phfx+'MDAxis']
3259            sumMD = 0
3260            for H in refl[11]:           
3261                cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G)
3262                A = 1.0/np.sqrt((MD*cosP)**2+sinP**2/MD)
3263                sumMD += A**3
3264            POcorr = sumMD/len(refl[11])
3265    else:   #spherical harmonics
3266        if calcControls[phfx+'SHord']:
3267            POcorr = SHPOcal(refl,g,phfx,hfx,SGData,calcControls,parmDict)
3268    return POcorr
3269   
3270def GetPrefOriDerv(refl,G,g,phfx,hfx,SGData,calcControls,parmDict):
3271    POcorr = 1.0
3272    POderv = {}
3273    if calcControls[phfx+'poType'] == 'MD':
3274        MD = parmDict[phfx+'MD']
3275        MDAxis = calcControls[phfx+'MDAxis']
3276        sumMD = 0
3277        sumdMD = 0
3278        for H in refl[11]:           
3279            cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G)
3280            A = 1.0/np.sqrt((MD*cosP)**2+sinP**2/MD)
3281            sumMD += A**3
3282            sumdMD -= (1.5*A**5)*(2.0*MD*cosP**2-(sinP/MD)**2)
3283        POcorr = sumMD/len(refl[11])
3284        POderv[phfx+'MD'] = sumdMD/len(refl[11])
3285    else:   #spherical harmonics
3286        if calcControls[phfx+'SHord']:
3287            POcorr,POderv = SHPOcalDerv(refl,g,phfx,hfx,SGData,calcControls,parmDict)
3288    return POcorr,POderv
3289   
3290def GetAbsorb(refl,hfx,calcControls,parmDict):
3291    if 'Debye' in calcControls[hfx+'instType']:
3292        return G2pwd.Absorb('Cylinder',parmDict[hfx+'Absorption'],refl[5],0,0)
3293    else:
3294        return 1.0
3295   
3296def GetAbsorbDerv(refl,hfx,calcControls,parmDict):
3297    if 'Debye' in calcControls[hfx+'instType']:
3298        return G2pwd.AbsorbDerv('Cylinder',parmDict[hfx+'Absorption'],refl[5],0,0)
3299    else:
3300        return 0.0
3301   
3302def GetIntensityCorr(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict):
3303    Icorr = parmDict[phfx+'Scale']*parmDict[hfx+'Scale']*refl[3]               #scale*multiplicity
3304    if 'X' in parmDict[hfx+'Type']:
3305        Icorr *= G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5],parmDict[hfx+'Azimuth'])[0]
3306    Icorr *= GetPrefOri(refl,G,g,phfx,hfx,SGData,calcControls,parmDict)
3307    if pfx+'SHorder' in parmDict:
3308        Icorr *= SHTXcal(refl,g,pfx,hfx,SGData,calcControls,parmDict)
3309    Icorr *= GetAbsorb(refl,hfx,calcControls,parmDict)
3310    refl[13] = Icorr       
3311   
3312def GetIntensityDerv(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict):
3313    dIdsh = 1./parmDict[hfx+'Scale']
3314    dIdsp = 1./parmDict[phfx+'Scale']
3315    if 'X' in parmDict[hfx+'Type']:
3316        pola,dIdPola = G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5],parmDict[hfx+'Azimuth'])
3317        dIdPola /= pola
3318    else:       #'N'
3319        dIdPola = 0.0
3320    POcorr,dIdPO = GetPrefOriDerv(refl,G,g,phfx,hfx,SGData,calcControls,parmDict)
3321    for iPO in dIdPO:
3322        dIdPO[iPO] /= POcorr
3323    dFdODF = {}
3324    dFdSA = [0,0,0]
3325    if pfx+'SHorder' in parmDict:
3326        odfCor,dFdODF,dFdSA = SHTXcalDerv(refl,g,pfx,hfx,SGData,calcControls,parmDict)
3327        for iSH in dFdODF:
3328            dFdODF[iSH] /= odfCor
3329        for i in range(3):
3330            dFdSA[i] /= odfCor
3331    dFdAb = GetAbsorbDerv(refl,hfx,calcControls,parmDict)
3332    return dIdsh,dIdsp,dIdPola,dIdPO,dFdODF,dFdSA,dFdAb
3333       
3334def GetSampleSigGam(refl,wave,G,GB,phfx,calcControls,parmDict):
3335    costh = cosd(refl[5]/2.)
3336    #crystallite size
3337    if calcControls[phfx+'SizeType'] == 'isotropic':
3338        Sgam = 1.8*wave/(np.pi*parmDict[phfx+'Size;i']*costh)
3339    elif calcControls[phfx+'SizeType'] == 'uniaxial':
3340        H = np.array(refl[:3])
3341        P = np.array(calcControls[phfx+'SizeAxis'])
3342        cosP,sinP = G2lat.CosSinAngle(H,P,G)
3343        Sgam = (1.8*wave/np.pi)/(parmDict[phfx+'Size;i']*parmDict[phfx+'Size;a']*costh)
3344        Sgam *= np.sqrt((sinP*parmDict[phfx+'Size;a'])**2+(cosP*parmDict[phfx+'Size;i'])**2)
3345    else:           #ellipsoidal crystallites
3346        Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
3347        H = np.array(refl[:3])
3348        lenR = G2pwd.ellipseSize(H,Sij,GB)
3349        Sgam = 1.8*wave/(np.pi*costh*lenR)
3350    #microstrain               
3351    if calcControls[phfx+'MustrainType'] == 'isotropic':
3352        Mgam = 0.018*parmDict[phfx+'Mustrain;i']*tand(refl[5]/2.)/np.pi
3353    elif calcControls[phfx+'MustrainType'] == 'uniaxial':
3354        H = np.array(refl[:3])
3355        P = np.array(calcControls[phfx+'MustrainAxis'])
3356        cosP,sinP = G2lat.CosSinAngle(H,P,G)
3357        Si = parmDict[phfx+'Mustrain;i']
3358        Sa = parmDict[phfx+'Mustrain;a']
3359        Mgam = 0.018*Si*Sa*tand(refl[5]/2.)/(np.pi*np.sqrt((Si*cosP)**2+(Sa*sinP)**2))
3360    else:       #generalized - P.W. Stephens model
3361        pwrs = calcControls[phfx+'MuPwrs']
3362        sum = 0
3363        for i,pwr in enumerate(pwrs):
3364            sum += parmDict[phfx+'Mustrain:'+str(i)]*refl[0]**pwr[0]*refl[1]**pwr[1]*refl[2]**pwr[2]
3365        Mgam = 0.018*refl[4]**2*tand(refl[5]/2.)*sum
3366    gam = Sgam*parmDict[phfx+'Size;mx']+Mgam*parmDict[phfx+'Mustrain;mx']
3367    sig = (Sgam*(1.-parmDict[phfx+'Size;mx']))**2+(Mgam*(1.-parmDict[phfx+'Mustrain;mx']))**2
3368    sig /= ateln2
3369    return sig,gam
3370       
3371def GetSampleSigGamDerv(refl,wave,G,GB,phfx,calcControls,parmDict):
3372    gamDict = {}
3373    sigDict = {}
3374    costh = cosd(refl[5]/2.)
3375    tanth = tand(refl[5]/2.)
3376    #crystallite size derivatives
3377    if calcControls[phfx+'SizeType'] == 'isotropic':
3378        Sgam = 1.8*wave/(np.pi*parmDict[phfx+'Size;i']*costh)
3379        gamDict[phfx+'Size;i'] = -1.8*wave*parmDict[phfx+'Size;mx']/(np.pi*costh)
3380        sigDict[phfx+'Size;i'] = -3.6*Sgam*wave*(1.-parmDict[phfx+'Size;mx'])**2/(np.pi*costh*ateln2)
3381    elif calcControls[phfx+'SizeType'] == 'uniaxial':
3382        H = np.array(refl[:3])
3383        P = np.array(calcControls[phfx+'SizeAxis'])
3384        cosP,sinP = G2lat.CosSinAngle(H,P,G)
3385        Si = parmDict[phfx+'Size;i']
3386        Sa = parmDict[phfx+'Size;a']
3387        gami = (1.8*wave/np.pi)/(Si*Sa)
3388        sqtrm = np.sqrt((sinP*Sa)**2+(cosP*Si)**2)
3389        Sgam = gami*sqtrm
3390        gam = Sgam/costh
3391        dsi = (gami*Si*cosP**2/(sqtrm*costh)-gam/Si)
3392        dsa = (gami*Sa*sinP**2/(sqtrm*costh)-gam/Sa)
3393        gamDict[phfx+'Size;i'] = dsi*parmDict[phfx+'Size;mx']
3394        gamDict[phfx+'Size;a'] = dsa*parmDict[phfx+'Size;mx']
3395        sigDict[phfx+'Size;i'] = 2.*dsi*Sgam*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
3396        sigDict[phfx+'Size;a'] = 2.*dsa*Sgam*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
3397    else:           #ellipsoidal crystallites
3398        const = 1.8*wave/(np.pi*costh)
3399        Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
3400        H = np.array(refl[:3])
3401        lenR,dRdS = G2pwd.ellipseSizeDerv(H,Sij,GB)
3402        Sgam = 1.8*wave/(np.pi*costh*lenR)
3403        for i,item in enumerate([phfx+'Size:%d'%(j) for j in range(6)]):
3404            gamDict[item] = -(const/lenR**2)*dRdS[i]*parmDict[phfx+'Size;mx']
3405            sigDict[item] = -2.*Sgam*(const/lenR**2)*dRdS[i]*(1.-parmDict[phfx+'Size;mx'])**2/ateln2
3406    gamDict[phfx+'Size;mx'] = Sgam
3407    sigDict[phfx+'Size;mx'] = -2.*Sgam**2*(1.-parmDict[phfx+'Size;mx'])/ateln2
3408           
3409    #microstrain derivatives               
3410    if calcControls[phfx+'MustrainType'] == 'isotropic':
3411        Mgam = 0.018*parmDict[phfx+'Mustrain;i']*tand(refl[5]/2.)/np.pi
3412        gamDict[phfx+'Mustrain;i'] =  0.018*tanth*parmDict[phfx+'Mustrain;mx']/np.pi
3413        sigDict[phfx+'Mustrain;i'] =  0.036*Mgam*tanth*(1.-parmDict[phfx+'Mustrain;mx'])**2/(np.pi*ateln2)       
3414    elif calcControls[phfx+'MustrainType'] == 'uniaxial':
3415        H = np.array(refl[:3])
3416        P = np.array(calcControls[phfx+'MustrainAxis'])
3417        cosP,sinP = G2lat.CosSinAngle(H,P,G)
3418        Si = parmDict[phfx+'Mustrain;i']
3419        Sa = parmDict[phfx+'Mustrain;a']
3420        gami = 0.018*Si*Sa*tanth/np.pi
3421        sqtrm = np.sqrt((Si*cosP)**2+(Sa*sinP)**2)
3422        Mgam = gami/sqtrm
3423        dsi = -gami*Si*cosP**2/sqtrm**3
3424        dsa = -gami*Sa*sinP**2/sqtrm**3
3425        gamDict[phfx+'Mustrain;i'] = (Mgam/Si+dsi)*parmDict[phfx+'Mustrain;mx']
3426        gamDict[phfx+'Mustrain;a'] = (Mgam/Sa+dsa)*parmDict[phfx+'Mustrain;mx']
3427        sigDict[phfx+'Mustrain;i'] = 2*(Mgam/Si+dsi)*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2
3428        sigDict[phfx+'Mustrain;a'] = 2*(Mgam/Sa+dsa)*Mgam*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2       
3429    else:       #generalized - P.W. Stephens model
3430        pwrs = calcControls[phfx+'MuPwrs']
3431        const = 0.018*refl[4]**2*tanth
3432        sum = 0
3433        for i,pwr in enumerate(pwrs):
3434            term = refl[0]**pwr[0]*refl[1]**pwr[1]*refl[2]**pwr[2]
3435            sum += parmDict[phfx+'Mustrain:'+str(i)]*term
3436            gamDict[phfx+'Mustrain:'+str(i)] = const*term*parmDict[phfx+'Mustrain;mx']
3437            sigDict[phfx+'Mustrain:'+str(i)] = \
3438                2.*const*term*(1.-parmDict[phfx+'Mustrain;mx'])**2/ateln2
3439        Mgam = 0.018*refl[4]**2*tand(refl[5]/2.)*sum
3440        for i in range(len(pwrs)):
3441            sigDict[phfx+'Mustrain:'+str(i)] *= Mgam
3442    gamDict[phfx+'Mustrain;mx'] = Mgam
3443    sigDict[phfx+'Mustrain;mx'] = -2.*Mgam**2*(1.-parmDict[phfx+'Mustrain;mx'])/ateln2
3444    return sigDict,gamDict
3445       
3446def GetReflPos(refl,wave,G,hfx,calcControls,parmDict):
3447    h,k,l = refl[:3]
3448    dsq = 1./G2lat.calc_rDsq2(np.array([h,k,l]),G)
3449    d = np.sqrt(dsq)
3450
3451    refl[4] = d
3452    pos = 2.0*asind(wave/(2.0*d))+parmDict[hfx+'Zero']
3453    const = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius'])                  #shifts in microns
3454    if 'Bragg' in calcControls[hfx+'instType']:
3455        pos -= const*(4.*parmDict[hfx+'Shift']*cosd(pos/2.0)+ \
3456            parmDict[hfx+'Transparency']*sind(pos)*100.0)            #trans(=1/mueff) in cm
3457    else:               #Debye-Scherrer - simple but maybe not right
3458        pos -= const*(parmDict[hfx+'DisplaceX']*cosd(pos)+parmDict[hfx+'DisplaceY']*sind(pos))
3459    return pos
3460
3461def GetReflPosDerv(refl,wave,A,hfx,calcControls,parmDict):
3462    dpr = 180./np.pi
3463    h,k,l = refl[:3]
3464    dstsq = G2lat.calc_rDsq(np.array([h,k,l]),A)
3465    dst = np.sqrt(dstsq)
3466    pos = refl[5]-parmDict[hfx+'Zero']
3467    const = dpr/np.sqrt(1.0-wave**2*dstsq/4.0)
3468    dpdw = const*dst
3469    dpdA = np.array([h**2,k**2,l**2,h*k,h*l,k*l])
3470    dpdA *= const*wave/(2.0*dst)
3471    dpdZ = 1.0
3472    const = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius'])                  #shifts in microns
3473    if 'Bragg' in calcControls[hfx+'instType']:
3474        dpdSh = -4.*const*cosd(pos/2.0)
3475        dpdTr = -const*sind(pos)*100.0
3476        return dpdA,dpdw,dpdZ,dpdSh,dpdTr,0.,0.
3477    else:               #Debye-Scherrer - simple but maybe not right
3478        dpdXd = -const*cosd(pos)
3479        dpdYd = -const*sind(pos)
3480        return dpdA,dpdw,dpdZ,0.,0.,dpdXd,dpdYd
3481           
3482def GetHStrainShift(refl,SGData,phfx,parmDict):
3483    laue = SGData['SGLaue']
3484    uniq = SGData['SGUniq']
3485    h,k,l = refl[:3]
3486    if laue in ['m3','m3m']:
3487        Dij = parmDict[phfx+'D11']*(h**2+k**2+l**2)+ \
3488            refl[4]**2*parmDict[phfx+'eA']*((h*k)**2+(h*l)**2+(k*l)**2)/(h**2+k**2+l**2)**2
3489    elif laue in ['6/m','6/mmm','3m1','31m','3']:
3490        Dij = parmDict[phfx+'D11']*(h**2+k**2+h*k)+parmDict[phfx+'D33']*l**2
3491    elif laue in ['3R','3mR']:
3492        Dij = parmDict[phfx+'D11']*(h**2+k**2+l**2)+parmDict[phfx+'D12']*(h*k+h*l+k*l)
3493    elif laue in ['4/m','4/mmm']:
3494        Dij = parmDict[phfx+'D11']*(h**2+k**2)+parmDict[phfx+'D33']*l**2
3495    elif laue in ['mmm']:
3496        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2
3497    elif laue in ['2/m']:
3498        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2
3499        if uniq == 'a':
3500            Dij += parmDict[phfx+'D23']*k*l
3501        elif uniq == 'b':
3502            Dij += parmDict[phfx+'D13']*h*l
3503        elif uniq == 'c':
3504            Dij += parmDict[phfx+'D12']*h*k
3505    else:
3506        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2+ \
3507            parmDict[phfx+'D12']*h*k+parmDict[phfx+'D13']*h*l+parmDict[phfx+'D23']*k*l
3508    return -Dij*refl[4]**2*tand(refl[5]/2.0)
3509           
3510def GetHStrainShiftDerv(refl,SGData,phfx):
3511    laue = SGData['SGLaue']
3512    uniq = SGData['SGUniq']
3513    h,k,l = refl[:3]
3514    if laue in ['m3','m3m']:
3515        dDijDict = {phfx+'D11':h**2+k**2+l**2,
3516            phfx+'eA':refl[4]**2*((h*k)**2+(h*l)**2+(k*l)**2)/(h**2+k**2+l**2)**2}
3517    elif laue in ['6/m','6/mmm','3m1','31m','3']:
3518        dDijDict = {phfx+'D11':h**2+k**2+h*k,phfx+'D33':l**2}
3519    elif laue in ['3R','3mR']:
3520        dDijDict = {phfx+'D11':h**2+k**2+l**2,phfx+'D12':h*k+h*l+k*l}
3521    elif laue in ['4/m','4/mmm']:
3522        dDijDict = {phfx+'D11':h**2+k**2,phfx+'D33':l**2}
3523    elif laue in ['mmm']:
3524        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2}
3525    elif laue in ['2/m']:
3526        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2}
3527        if uniq == 'a':
3528            dDijDict[phfx+'D23'] = k*l
3529        elif uniq == 'b':
3530            dDijDict[phfx+'D13'] = h*l
3531        elif uniq == 'c':
3532            dDijDict[phfx+'D12'] = h*k
3533            names.append()
3534    else:
3535        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2,
3536            phfx+'D12':h*k,phfx+'D13':h*l,phfx+'D23':k*l}
3537    for item in dDijDict:
3538        dDijDict[item] *= -refl[4]**2*tand(refl[5]/2.0)
3539    return dDijDict
3540   
3541def GetFprime(controlDict,Histograms):
3542    FFtables = controlDict['FFtables']
3543    if not FFtables:
3544        return
3545    histoList = Histograms.keys()
3546    histoList.sort()
3547    for histogram in histoList:
3548        if histogram[:4] in ['PWDR','HKLF']:
3549            Histogram = Histograms[histogram]
3550            hId = Histogram['hId']
3551            hfx = ':%d:'%(hId)
3552            keV = controlDict[hfx+'keV']
3553            for El in FFtables:
3554                Orbs = G2el.GetXsectionCoeff(El.split('+')[0].split('-')[0])
3555                FP,FPP,Mu = G2el.FPcalc(Orbs, keV)
3556                FFtables[El][hfx+'FP'] = FP
3557                FFtables[El][hfx+'FPP'] = FPP               
3558           
3559def GetFobsSq(Histograms,Phases,parmDict,calcControls):
3560    histoList = Histograms.keys()
3561    histoList.sort()
3562    for histogram in histoList:
3563        if 'PWDR' in histogram[:4]:
3564            Histogram = Histograms[histogram]
3565            hId = Histogram['hId']
3566            hfx = ':%d:'%(hId)
3567            Limits = calcControls[hfx+'Limits']
3568            shl = max(parmDict[hfx+'SH/L'],0.0005)
3569            Ka2 = False
3570            kRatio = 0.0
3571            if hfx+'Lam1' in parmDict.keys():
3572                Ka2 = True
3573                lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
3574                kRatio = parmDict[hfx+'I(L2)/I(L1)']
3575            x,y,w,yc,yb,yd = Histogram['Data']
3576            xB = np.searchsorted(x,Limits[0])
3577            xF = np.searchsorted(x,Limits[1])
3578            ymb = np.array(y-yb)
3579            ymb = np.where(ymb,ymb,1.0)
3580            ycmb = np.array(yc-yb)
3581            ratio = 1./np.where(ycmb,ycmb/ymb,1.e10)         
3582            refLists = Histogram['Reflection Lists']
3583            for phase in refLists:
3584                Phase = Phases[phase]
3585                pId = Phase['pId']
3586                phfx = '%d:%d:'%(pId,hId)
3587                refList = refLists[phase]
3588                sumFo = 0.0
3589                sumdF = 0.0
3590                sumFosq = 0.0
3591                sumdFsq = 0.0
3592                for refl in refList:
3593                    if 'C' in calcControls[hfx+'histType']:
3594                        yp = np.zeros_like(yb)
3595                        Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5],refl[6],refl[7],shl)
3596                        iBeg = max(xB,np.searchsorted(x,refl[5]-fmin))
3597                        iFin = max(xB,min(np.searchsorted(x,refl[5]+fmax),xF))
3598                        iFin2 = iFin
3599                        yp[iBeg:iFin] = refl[13]*refl[9]*G2pwd.getFCJVoigt3(refl[5],refl[6],refl[7],shl,x[iBeg:iFin])    #>90% of time spent here
3600                        if Ka2:
3601                            pos2 = refl[5]+lamRatio*tand(refl[5]/2.0)       # + 360/pi * Dlam/lam * tan(th)
3602                            Wd,fmin,fmax = G2pwd.getWidthsCW(pos2,refl[6],refl[7],shl)
3603                            iBeg2 = max(xB,np.searchsorted(x,pos2-fmin))
3604                            iFin2 = min(np.searchsorted(x,pos2+fmax),xF)
3605                            if not iBeg2+iFin2:       #peak below low limit - skip peak
3606                                continue
3607                            elif not iBeg2-iFin2:     #peak above high limit - done
3608                                break
3609                            yp[iBeg2:iFin2] += refl[13]*refl[9]*kRatio*G2pwd.getFCJVoigt3(pos2,refl[6],refl[7],shl,x[iBeg2:iFin2])        #and here
3610                        refl[8] = np.sum(np.where(ratio[iBeg:iFin2]>0.,yp[iBeg:iFin2]*ratio[iBeg:iFin2]/(refl[13]*(1.+kRatio)),0.0))
3611                    elif 'T' in calcControls[hfx+'histType']:
3612                        print 'TOF Undefined at present'
3613                        raise Exception    #no TOF yet
3614                    Fo = np.sqrt(np.abs(refl[8]))
3615                    Fc = np.sqrt(np.abs(refl[9]))
3616                    sumFo += Fo
3617                    sumFosq += refl[8]**2
3618                    sumdF += np.abs(Fo-Fc)
3619                    sumdFsq += (refl[8]-refl[9])**2
3620                Histogram[phfx+'Rf'] = min(100.,(sumdF/sumFo)*100.)
3621                Histogram[phfx+'Rf^2'] = min(100.,np.sqrt(sumdFsq/sumFosq)*100.)
3622                Histogram[phfx+'Nref'] = len(refList)
3623               
3624def getPowderProfile(parmDict,x,varylist,Histogram,Phases,calcControls,pawleyLookup):
3625   
3626    def GetReflSigGam(refl,wave,G,GB,hfx,phfx,calcControls,parmDict):
3627        U = parmDict[hfx+'U']
3628        V = parmDict[hfx+'V']
3629        W = parmDict[hfx+'W']
3630        X = parmDict[hfx+'X']
3631        Y = parmDict[hfx+'Y']
3632        tanPos = tand(refl[5]/2.0)
3633        Ssig,Sgam = GetSampleSigGam(refl,wave,G,GB,phfx,calcControls,parmDict)
3634        sig = U*tanPos**2+V*tanPos+W+Ssig     #save peak sigma
3635        sig = max(0.001,sig)
3636        gam = X/cosd(refl[5]/2.0)+Y*tanPos+Sgam     #save peak gamma
3637        gam = max(0.001,gam)
3638        return sig,gam
3639               
3640    hId = Histogram['hId']
3641    hfx = ':%d:'%(hId)
3642    bakType = calcControls[hfx+'bakType']
3643    yb = G2pwd.getBackground(hfx,parmDict,bakType,x)
3644    yc = np.zeros_like(yb)
3645       
3646    if 'C' in calcControls[hfx+'histType']:   
3647        shl = max(parmDict[hfx+'SH/L'],0.002)
3648        Ka2 = False
3649        if hfx+'Lam1' in parmDict.keys():
3650            wave = parmDict[hfx+'Lam1']
3651            Ka2 = True
3652            lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
3653            kRatio = parmDict[hfx+'I(L2)/I(L1)']
3654        else:
3655            wave = parmDict[hfx+'Lam']
3656    else:
3657        print 'TOF Undefined at present'
3658        raise ValueError
3659    for phase in Histogram['Reflection Lists']:
3660        refList = Histogram['Reflection Lists'][phase]
3661        Phase = Phases[phase]
3662        pId = Phase['pId']
3663        pfx = '%d::'%(pId)
3664        phfx = '%d:%d:'%(pId,hId)
3665        hfx = ':%d:'%(hId)
3666        SGData = Phase['General']['SGData']
3667        A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
3668        G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
3669        GA,GB = G2lat.Gmat2AB(G)    #Orthogonalization matricies
3670        Vst = np.sqrt(nl.det(G))    #V*
3671        if not Phase['General'].get('doPawley'):
3672            refList = StructureFactor(refList,G,hfx,pfx,SGData,calcControls,parmDict)
3673        for refl in refList:
3674            if 'C' in calcControls[hfx+'histType']:
3675                h,k,l = refl[:3]
3676                refl[5] = GetReflPos(refl,wave,G,hfx,calcControls,parmDict)         #corrected reflection position
3677                Lorenz = 1./(2.*sind(refl[5]/2.)**2*cosd(refl[5]/2.))           #Lorentz correction
3678                refl[5] += GetHStrainShift(refl,SGData,phfx,parmDict)               #apply hydrostatic strain shift
3679                refl[6:8] = GetReflSigGam(refl,wave,G,GB,hfx,phfx,calcControls,parmDict)    #peak sig & gam
3680                GetIntensityCorr(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)    #puts corrections in refl[13]
3681                refl[13] *= Vst*Lorenz
3682                if Phase['General'].get('doPawley'):
3683                    try:
3684                        pInd =pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])
3685                        refl[9] = parmDict[pInd]
3686                    except KeyError:
3687#                        print ' ***Error %d,%d,%d missing from Pawley reflection list ***'%(h,k,l)
3688                        continue
3689                Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5],refl[6],refl[7],shl)
3690                iBeg = np.searchsorted(x,refl[5]-fmin)
3691                iFin = np.searchsorted(x,refl[5]+fmax)
3692                if not iBeg+iFin:       #peak below low limit - skip peak
3693                    continue
3694                elif not iBeg-iFin:     #peak above high limit - done
3695                    break
3696                yc[iBeg:iFin] += refl[13]*refl[9]*G2pwd.getFCJVoigt3(refl[5],refl[6],refl[7],shl,x[iBeg:iFin])    #>90% of time spent here
3697                if Ka2:
3698                    pos2 = refl[5]+lamRatio*tand(refl[5]/2.0)       # + 360/pi * Dlam/lam * tan(th)
3699                    Wd,fmin,fmax = G2pwd.getWidthsCW(pos2,refl[6],refl[7],shl)
3700                    iBeg = np.searchsorted(x,pos2-fmin)
3701                    iFin = np.searchsorted(x,pos2+fmax)
3702                    if not iBeg+iFin:       #peak below low limit - skip peak
3703                        continue
3704                    elif not iBeg-iFin:     #peak above high limit - done
3705                        return yc,yb
3706                    yc[iBeg:iFin] += refl[13]*refl[9]*kRatio*G2pwd.getFCJVoigt3(pos2,refl[6],refl[7],shl,x[iBeg:iFin])        #and here
3707            elif 'T' in calcControls[hfx+'histType']:
3708                print 'TOF Undefined at present'
3709                raise Exception    #no TOF yet
3710    return yc,yb
3711   
3712def getPowderProfileDerv(parmDict,x,varylist,Histogram,Phases,rigidbodyDict,calcControls,pawleyLookup):
3713   
3714    def cellVaryDerv(pfx,SGData,dpdA): 
3715        if SGData['SGLaue'] in ['-1',]:
3716            return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],
3717                [pfx+'A3',dpdA[3]],[pfx+'A4',dpdA[4]],[pfx+'A5',dpdA[5]]]
3718        elif SGData['SGLaue'] in ['2/m',]:
3719            if SGData['SGUniq'] == 'a':
3720                return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A3',dpdA[3]]]
3721            elif SGData['SGUniq'] == 'b':
3722                return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A4',dpdA[4]]]
3723            else:
3724                return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A5',dpdA[5]]]
3725        elif SGData['SGLaue'] in ['mmm',]:
3726            return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]]]
3727        elif SGData['SGLaue'] in ['4/m','4/mmm']:
3728            return [[pfx+'A0',dpdA[0]],[pfx+'A2',dpdA[2]]]
3729        elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
3730            return [[pfx+'A0',dpdA[0]],[pfx+'A2',dpdA[2]]]
3731        elif SGData['SGLaue'] in ['3R', '3mR']:
3732            return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[2]],[pfx+'A3',dpdA[3]+dpdA[4]+dpdA[5]]]                       
3733        elif SGData['SGLaue'] in ['m3m','m3']:
3734            return [[pfx+'A0',dpdA[0]]]
3735    # create a list of dependent variables and set up a dictionary to hold their derivatives
3736    dependentVars = G2mv.GetDependentVars()
3737    depDerivDict = {}
3738    for j in dependentVars:
3739        depDerivDict[j] = np.zeros(shape=(len(x)))
3740    #print 'dependent vars',dependentVars
3741    lenX = len(x)               
3742    hId = Histogram['hId']
3743    hfx = ':%d:'%(hId)
3744    bakType = calcControls[hfx+'bakType']
3745    dMdv = np.zeros(shape=(len(varylist),len(x)))
3746    dMdb,dMddb,dMdpk = G2pwd.getBackgroundDerv(hfx,parmDict,bakType,x)
3747    if hfx+'Back:0' in varylist: # for now assume that Back:x vars to not appear in constraints
3748        bBpos =varylist.index(hfx+'Back:0')
3749        dMdv[bBpos:bBpos+len(dMdb)] = dMdb
3750    names = [hfx+'DebyeA',hfx+'DebyeR',hfx+'DebyeU']
3751    for name in varylist:
3752        if 'Debye' in name:
3753            id = int(name.split(':')[-1])
3754            parm = name[:int(name.rindex(':'))]
3755            ip = names.index(parm)
3756            dMdv[varylist.index(name)] = dMddb[3*id+ip]
3757    names = [hfx+'BkPkpos',hfx+'BkPkint',hfx+'BkPksig',hfx+'BkPkgam']
3758    for name in varylist:
3759        if 'BkPk' in name:
3760            id = int(name.split(':')[-1])
3761            parm = name[:int(name.rindex(':'))]
3762            ip = names.index(parm)
3763            dMdv[varylist.index(name)] = dMdpk[4*id+ip]
3764    cw = np.diff(x)
3765    cw = np.append(cw,cw[-1])
3766    if 'C' in calcControls[hfx+'histType']:   
3767        shl = max(parmDict[hfx+'SH/L'],0.002)
3768        Ka2 = False
3769        if hfx+'Lam1' in parmDict.keys():
3770            wave = parmDict[hfx+'Lam1']
3771            Ka2 = True
3772            lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
3773            kRatio = parmDict[hfx+'I(L2)/I(L1)']
3774        else:
3775            wave = parmDict[hfx+'Lam']
3776    else:
3777        print 'TOF Undefined at present'
3778        raise ValueError
3779    for phase in Histogram['Reflection Lists']:
3780        refList = Histogram['Reflection Lists'][phase]
3781        Phase = Phases[phase]
3782        SGData = Phase['General']['SGData']
3783        pId = Phase['pId']
3784        pfx = '%d::'%(pId)
3785        phfx = '%d:%d:'%(pId,hId)
3786        A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
3787        G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
3788        GA,GB = G2lat.Gmat2AB(G)    #Orthogonalization matricies
3789        if not Phase['General'].get('doPawley'):
3790            dFdvDict = StructureFactorDerv(refList,G,hfx,pfx,SGData,calcControls,parmDict)
3791            ApplyRBModelDervs(dFdvDict,parmDict,rigidbodyDict,Phase)
3792        for iref,refl in enumerate(refList):
3793            if 'C' in calcControls[hfx+'histType']:        #CW powder
3794                h,k,l = refl[:3]
3795                dIdsh,dIdsp,dIdpola,dIdPO,dFdODF,dFdSA,dFdAb = GetIntensityDerv(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
3796                Wd,fmin,fmax = G2pwd.getWidthsCW(refl[5],refl[6],refl[7],shl)
3797                iBeg = np.searchsorted(x,refl[5]-fmin)
3798                iFin = np.searchsorted(x,refl[5]+fmax)
3799                if not iBeg+iFin:       #peak below low limit - skip peak
3800                    continue
3801                elif not iBeg-iFin:     #peak above high limit - done
3802                    break
3803                pos = refl[5]
3804                tanth = tand(pos/2.0)
3805                costh = cosd(pos/2.0)
3806                lenBF = iFin-iBeg
3807                dMdpk = np.zeros(shape=(6,lenBF))
3808                dMdipk = G2pwd.getdFCJVoigt3(refl[5],refl[6],refl[7],shl,x[iBeg:iFin])
3809                for i in range(5):
3810                    dMdpk[i] += 100.*cw[iBeg:iFin]*refl[13]*refl[9]*dMdipk[i]
3811                dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':np.zeros_like(dMdpk[0])}
3812                if Ka2:
3813                    pos2 = refl[5]+lamRatio*tanth       # + 360/pi * Dlam/lam * tan(th)
3814                    iBeg2 = np.searchsorted(x,pos2-fmin)
3815                    iFin2 = np.searchsorted(x,pos2+fmax)
3816                    if iBeg2-iFin2:
3817                        lenBF2 = iFin2-iBeg2
3818                        dMdpk2 = np.zeros(shape=(6,lenBF2))
3819                        dMdipk2 = G2pwd.getdFCJVoigt3(pos2,refl[6],refl[7],shl,x[iBeg2