source: trunk/GSASIIstruct.py @ 915

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

RB now seem OK, SC now seem OK for F & F2 refinements.

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