source: trunk/GSASIIstruct.py @ 885

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

add 'Global' to constraints
remove a stray debug print
fix U6toUij so it returns a numpy array
fixes to Q2AV & Q2AVdeg
RB changes - get derivatives in

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