source: trunk/GSASIIstruct.py @ 457

Last change on this file since 457 was 457, checked in by vondreele, 10 years ago

begin distance angle calcs
move gpxfile routines from GSASIIstruct.py to GSASIIIO.py
move getVCov & ValEsd? to GSASIImath.py
add some text to help/gsasII.html

  • Property svn:keywords set to Date Author Revision URL Id
File size: 119.2 KB
Line 
1#GSASIIstructure - structure computation routines
2########### SVN repository information ###################
3# $Date: 2012-01-25 15:24:20 +0000 (Wed, 25 Jan 2012) $
4# $Author: vondreele $
5# $Revision: 457 $
6# $URL: trunk/GSASIIstruct.py $
7# $Id: GSASIIstruct.py 457 2012-01-25 15:24:20Z vondreele $
8########### SVN repository information ###################
9import sys
10import numpy as np
11import numpy.linalg as nl
12import time
13import math
14import GSASIIpath
15import GSASIIIO as G2IO
16import GSASIIElem as G2el
17import GSASIIlattice as G2lat
18import GSASIIspc as G2spc
19import GSASIIpwd as G2pwd
20import GSASIImapvars as G2mv
21import GSASIImath as G2mth
22import scipy.optimize as so
23
24sind = lambda x: np.sin(x*np.pi/180.)
25cosd = lambda x: np.cos(x*np.pi/180.)
26tand = lambda x: np.tan(x*np.pi/180.)
27asind = lambda x: 180.*np.arcsin(x)/np.pi
28atan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
29
30
31def ShowBanner():
32    print 80*'*'
33    print '   General Structure Analysis System-II Crystal Structure Refinement'
34    print '     by Robert B. Von Dreele, Argonne National Laboratory(C), 2010'
35    print ' This product includes software developed by the UChicago Argonne, LLC,' 
36    print '            as Operator of Argonne National Laboratory.'
37    print 80*'*','\n'
38
39def ShowControls(Controls):
40    print ' Least squares controls:'
41    print ' Derivative type: ',Controls['deriv type']
42    print ' Minimum delta-M/M for convergence: ','%.2g'%(Controls['min dM/M'])
43    print ' Initial shift factor: ','%.3f'%(Controls['shift factor'])
44   
45def GetFFtable(General):
46    ''' returns a dictionary of form factor data for atom types found in General
47    input:
48        General = dictionary of phase info.; includes AtomTypes
49    return:
50        FFtable = dictionary of form factor data; key is atom type
51    '''
52    atomTypes = General['AtomTypes']
53    FFtable = {}
54    for El in atomTypes:
55        FFs = G2el.GetFormFactorCoeff(El.split('+')[0].split('-')[0])
56        for item in FFs:
57            if item['Symbol'] == El.upper():
58                FFtable[El] = item
59    return FFtable
60   
61def GetBLtable(General):
62    ''' returns a dictionary of neutron scattering length data for atom types & isotopes found in General
63    input:
64        General = dictionary of phase info.; includes AtomTypes & Isotopes
65    return:
66        BLtable = dictionary of scattering length data; key is atom type
67    '''
68    atomTypes = General['AtomTypes']
69    BLtable = {}
70    isotopes = General['Isotopes']
71    isotope = General['Isotope']
72    for El in atomTypes:
73        BLtable[El] = [isotope[El],isotopes[El][isotope[El]]]
74    return BLtable
75       
76def GetPawleyConstr(SGLaue,PawleyRef,pawleyVary):
77    if SGLaue in ['-1','2/m','mmm']:
78        return                      #no Pawley symmetry required constraints
79    for i,varyI in enumerate(pawleyVary):
80        refI = int(varyI.split(':')[-1])
81        ih,ik,il = PawleyRef[refI][:3]
82        for varyJ in pawleyVary[0:i]:
83            refJ = int(varyJ.split(':')[-1])
84            jh,jk,jl = PawleyRef[refJ][:3]
85            if SGLaue in ['4/m','4/mmm']:
86                isum = ih**2+ik**2
87                jsum = jh**2+jk**2
88                if abs(il) == abs(jl) and isum == jsum:
89                    G2mv.StoreEquivalence(varyJ,(varyI,))
90            elif SGLaue in ['3R','3mR']:
91                isum = ih**2+ik**2+il**2
92                jsum = jh**2+jk**2*jl**2
93                isum2 = ih*ik+ih*il+ik*il
94                jsum2 = jh*jk+jh*jl+jk*jl
95                if isum == jsum and isum2 == jsum2:
96                    G2mv.StoreEquivalence(varyJ,(varyI,))
97            elif SGLaue in ['3','3m1','31m','6/m','6/mmm']:
98                isum = ih**2+ik**2+ih*ik
99                jsum = jh**2+jk**2+jh*jk
100                if abs(il) == abs(jl) and isum == jsum:
101                    G2mv.StoreEquivalence(varyJ,(varyI,))
102            elif SGLaue in ['m3','m3m']:
103                isum = ih**2+ik**2+il**2
104                jsum = jh**2+jk**2+jl**2
105                if isum == jsum:
106                    G2mv.StoreEquivalence(varyJ,(varyI,))
107                   
108def cellVary(pfx,SGData): 
109    if SGData['SGLaue'] in ['-1',]:
110        return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3',pfx+'A4',pfx+'A5']
111    elif SGData['SGLaue'] in ['2/m',]:
112        if SGData['SGUniq'] == 'a':
113            return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3']
114        elif SGData['SGUniq'] == 'b':
115            return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A4']
116        else:
117            return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A5']
118    elif SGData['SGLaue'] in ['mmm',]:
119        return [pfx+'A0',pfx+'A1',pfx+'A2']
120    elif SGData['SGLaue'] in ['4/m','4/mmm']:
121        return [pfx+'A0',pfx+'A2']
122    elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
123        return [pfx+'A0',pfx+'A2']
124    elif SGData['SGLaue'] in ['3R', '3mR']:
125        return [pfx+'A0',pfx+'A3']                       
126    elif SGData['SGLaue'] in ['m3m','m3']:
127        return [pfx+'A0',]
128                   
129def GetPhaseData(PhaseData,Print=True):
130           
131    def PrintFFtable(FFtable):
132        print '\n X-ray scattering factors:'
133        print '   Symbol     fa                                      fb                                      fc'
134        print 99*'-'
135        for Ename in FFtable:
136            ffdata = FFtable[Ename]
137            fa = ffdata['fa']
138            fb = ffdata['fb']
139            print ' %8s %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f' %  \
140                (Ename.ljust(8),fa[0],fa[1],fa[2],fa[3],fb[0],fb[1],fb[2],fb[3],ffdata['fc'])
141               
142    def PrintBLtable(BLtable):
143        print '\n Neutron scattering factors:'
144        print '   Symbol   isotope       mass       b       resonant terms'
145        print 99*'-'
146        for Ename in BLtable:
147            bldata = BLtable[Ename]
148            isotope = bldata[0]
149            mass = bldata[1][0]
150            blen = bldata[1][1]
151            bres = []
152            if len(bldata[1]) > 2:
153                bres = bldata[1][2:]
154            line = ' %8s%11s %10.3f %8.3f'%(Ename.ljust(8),isotope.center(11),mass,blen)
155            for item in bres:
156                line += '%10.5g'%(item)
157            print line
158               
159    def PrintAtoms(General,Atoms):
160        print '\n Atoms:'
161        line = '   name    type  refine?   x         y         z    '+ \
162            '  frac site sym  mult I/A   Uiso     U11     U22     U33     U12     U13     U23'
163        if General['Type'] == 'magnetic':
164            line += '   Mx     My     Mz'
165        elif General['Type'] == 'macromolecular':
166            line = ' res no  residue  chain '+line
167        print line
168        if General['Type'] == 'nuclear':
169            print 135*'-'
170            for i,at in enumerate(Atoms):
171                line = '%7s'%(at[0])+'%7s'%(at[1])+'%7s'%(at[2])+'%10.5f'%(at[3])+'%10.5f'%(at[4])+ \
172                    '%10.5f'%(at[5])+'%8.3f'%(at[6])+'%7s'%(at[7])+'%5d'%(at[8])+'%5s'%(at[9])
173                if at[9] == 'I':
174                    line += '%8.4f'%(at[10])+48*' '
175                else:
176                    line += 8*' '
177                    for j in range(6):
178                        line += '%8.4f'%(at[11+j])
179                print line
180       
181    def PrintTexture(textureData):
182        topstr = '\n Spherical harmonics texture: Order:' + \
183            str(textureData['Order'])
184        if textureData['Order']:
185            print topstr+' Refine? '+str(textureData['SH Coeff'][0])
186        else:
187            print topstr
188            return
189        names = ['omega','chi','phi']
190        line = '\n'
191        for name in names:
192            line += ' SH '+name+':'+'%12.4f'%(textureData['Sample '+name][1])+' Refine? '+str(textureData['Sample '+name][0])
193        print line
194        print '\n Texture coefficients:'
195        ptlbls = ' names :'
196        ptstr =  ' values:'
197        SHcoeff = textureData['SH Coeff'][1]
198        for item in SHcoeff:
199            ptlbls += '%12s'%(item)
200            ptstr += '%12.4f'%(SHcoeff[item]) 
201        print ptlbls
202        print ptstr   
203       
204    if Print: print ' Phases:'
205    phaseVary = []
206    phaseDict = {}
207    phaseConstr = {}
208    pawleyLookup = {}
209    FFtables = {}                   #scattering factors - xrays
210    BLtables = {}                   # neutrons
211    Natoms = {}
212    AtMults = {}
213    AtIA = {}
214    shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
215    SamSym = dict(zip(shModels,['0','-1','2/m','mmm']))
216    for name in PhaseData:
217        General = PhaseData[name]['General']
218        pId = PhaseData[name]['pId']
219        pfx = str(pId)+'::'
220        FFtable = GetFFtable(General)
221        BLtable = GetBLtable(General)
222        FFtables.update(FFtable)
223        BLtables.update(BLtable)
224        Atoms = PhaseData[name]['Atoms']
225        try:
226            PawleyRef = PhaseData[name]['Pawley ref']
227        except KeyError:
228            PawleyRef = []
229        SGData = General['SGData']
230        SGtext = G2spc.SGPrint(SGData)
231        cell = General['Cell']
232        A = G2lat.cell2A(cell[1:7])
233        phaseDict.update({pfx+'A0':A[0],pfx+'A1':A[1],pfx+'A2':A[2],pfx+'A3':A[3],pfx+'A4':A[4],pfx+'A5':A[5]})
234        if cell[0]:
235            phaseVary += cellVary(pfx,SGData)
236        Natoms[pfx] = 0
237        if Atoms:
238            if General['Type'] == 'nuclear':
239                Natoms[pfx] = len(Atoms)
240                for i,at in enumerate(Atoms):
241                    phaseDict.update({pfx+'Atype:'+str(i):at[1],pfx+'Afrac:'+str(i):at[6],pfx+'Amul:'+str(i):at[8],
242                        pfx+'Ax:'+str(i):at[3],pfx+'Ay:'+str(i):at[4],pfx+'Az:'+str(i):at[5],
243                        pfx+'dAx:'+str(i):0.,pfx+'dAy:'+str(i):0.,pfx+'dAz:'+str(i):0.,         #refined shifts for x,y,z
244                        pfx+'AI/A:'+str(i):at[9],})
245                    if at[9] == 'I':
246                        phaseDict[pfx+'AUiso:'+str(i)] = at[10]
247                    else:
248                        phaseDict.update({pfx+'AU11:'+str(i):at[11],pfx+'AU22:'+str(i):at[12],pfx+'AU33:'+str(i):at[13],
249                            pfx+'AU12:'+str(i):at[14],pfx+'AU13:'+str(i):at[15],pfx+'AU23:'+str(i):at[16]})
250                    if 'F' in at[2]:
251                        phaseVary.append(pfx+'Afrac:'+str(i))
252                    if 'X' in at[2]:
253                        xId,xCoef = G2spc.GetCSxinel(at[7])
254                        delnames = [pfx+'dAx:'+str(i),pfx+'dAy:'+str(i),pfx+'dAz:'+str(i)]
255                        for j in range(3):
256                            if xId[j] > 0:                               
257                                phaseVary.append(delnames[j])
258                                for k in range(j):
259                                    if xId[j] == xId[k]:
260                                        G2mv.StoreEquivalence(delnames[k],((delnames[j],xCoef[j]),)) 
261                    if 'U' in at[2]:
262                        if at[9] == 'I':
263                            phaseVary.append(pfx+'AUiso:'+str(i))
264                        else:
265                            uId,uCoef = G2spc.GetCSuinel(at[7])[:2]
266                            names = [pfx+'AU11:'+str(i),pfx+'AU22:'+str(i),pfx+'AU33:'+str(i),
267                                pfx+'AU12:'+str(i),pfx+'AU13:'+str(i),pfx+'AU23:'+str(i)]
268                            for j in range(6):
269                                if uId[j] > 0:                               
270                                    phaseVary.append(names[j])
271                                    for k in range(j):
272                                        if uId[j] == uId[k]:
273                                            G2mv.StoreEquivalence(names[k],((names[j],uCoef[j]),))
274#            elif General['Type'] == 'magnetic':
275#            elif General['Type'] == 'macromolecular':
276
277               
278            if 'SH Texture' in General:
279                textureData = General['SH Texture']
280                phaseDict[pfx+'SHmodel'] = SamSym[textureData['Model']]
281                phaseDict[pfx+'SHorder'] = textureData['Order']
282                for name in ['omega','chi','phi']:
283                    phaseDict[pfx+'SH '+name] = textureData['Sample '+name][1]
284                    if textureData['Sample '+name][0]:
285                        phaseVary.append(pfx+'SH '+name)
286                for name in textureData['SH Coeff'][1]:
287                    phaseDict[pfx+name] = textureData['SH Coeff'][1][name]
288                    if textureData['SH Coeff'][0]:
289                        phaseVary.append(pfx+name)
290               
291            if Print:
292                print '\n Phase name: ',General['Name']
293                print 135*'-'
294                PrintFFtable(FFtable)
295                PrintBLtable(BLtable)
296                print ''
297                for line in SGtext: print line
298                PrintAtoms(General,Atoms)
299                print '\n Unit cell: a =','%.5f'%(cell[1]),' b =','%.5f'%(cell[2]),' c =','%.5f'%(cell[3]), \
300                    ' alpha =','%.3f'%(cell[4]),' beta =','%.3f'%(cell[5]),' gamma =', \
301                    '%.3f'%(cell[6]),' volume =','%.3f'%(cell[7]),' Refine?',cell[0]
302                if 'SH Texture' in General:
303                    PrintTexture(textureData)
304                   
305        elif PawleyRef:
306            pawleyVary = []
307            for i,refl in enumerate(PawleyRef):
308                phaseDict[pfx+'PWLref:'+str(i)] = refl[6]
309                pawleyLookup[pfx+'%d,%d,%d'%(refl[0],refl[1],refl[2])] = i
310                if refl[5]:
311                    pawleyVary.append(pfx+'PWLref:'+str(i))
312            GetPawleyConstr(SGData['SGLaue'],PawleyRef,pawleyVary)      #does G2mv.StoreEquivalence
313            phaseVary += pawleyVary
314               
315    return Natoms,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables
316   
317def cellFill(pfx,SGData,parmDict,sigDict): 
318    if SGData['SGLaue'] in ['-1',]:
319        A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
320            parmDict[pfx+'A3'],parmDict[pfx+'A4'],parmDict[pfx+'A5']]
321        sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
322            sigDict[pfx+'A3'],sigDict[pfx+'A4'],sigDict[pfx+'A5']]
323    elif SGData['SGLaue'] in ['2/m',]:
324        if SGData['SGUniq'] == 'a':
325            A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
326                parmDict[pfx+'A3'],0,0]
327            sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
328                sigDict[pfx+'A3'],0,0]
329        elif SGData['SGUniq'] == 'b':
330            A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
331                0,parmDict[pfx+'A4'],0]
332            sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
333                0,sigDict[pfx+'A4'],0]
334        else:
335            A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
336                0,0,parmDict[pfx+'A5']]
337            sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
338                0,0,sigDict[pfx+'A5']]
339    elif SGData['SGLaue'] in ['mmm',]:
340        A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],0,0,0]
341        sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],0,0,0]
342    elif SGData['SGLaue'] in ['4/m','4/mmm']:
343        A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A2'],0,0,0]
344        sigA = [sigDict[pfx+'A0'],0,sigDict[pfx+'A2'],0,0,0]
345    elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
346        A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A2'],
347            parmDict[pfx+'A0'],0,0]
348        sigA = [sigDict[pfx+'A0'],0,sigDict[pfx+'A2'],0,0,0]
349    elif SGData['SGLaue'] in ['3R', '3mR']:
350        A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A0'],
351            parmDict[pfx+'A3'],parmDict[pfx+'A3'],parmDict[pfx+'A3']]
352        sigA = [sigDict[pfx+'A0'],0,0,sigDict[pfx+'A3'],0,0]
353    elif SGData['SGLaue'] in ['m3m','m3']:
354        A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A0'],0,0,0]
355        sigA = [sigDict[pfx+'A0'],0,0,0,0,0]
356    return A,sigA
357       
358def getCellEsd(pfx,SGData,A,covData):
359    dpr = 180./np.pi
360    rVsq = G2lat.calc_rVsq(A)
361    G,g = G2lat.A2Gmat(A)       #get recip. & real metric tensors
362    cell = np.array(G2lat.Gmat2cell(g))   #real cell
363    cellst = np.array(G2lat.Gmat2cell(G)) #recip. cell
364    scos = cosd(cellst[3:6])
365    ssin = sind(cellst[3:6])
366    scot = scos/ssin
367    rcos = cosd(cell[3:6])
368    rsin = sind(cell[3:6])
369    rcot = rcos/rsin
370    RMnames = [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3',pfx+'A4',pfx+'A5']
371    varyList = covData['varyList']
372    covMatrix = covData['covMatrix']
373    vcov = G2mth.getVCov(RMnames,varyList,covMatrix)
374    Ax = np.array(A)
375    Ax[3:] /= 2.
376    drVdA = np.array([Ax[1]*Ax[2]-Ax[5]**2,Ax[0]*Ax[2]-Ax[4]**2,Ax[0]*Ax[1]-Ax[3]**2,
377        Ax[4]*Ax[5]-Ax[2]*Ax[3],Ax[3]*Ax[5]-Ax[1]*Ax[4],Ax[3]*Ax[4]-Ax[0]*Ax[5]])
378    srcvlsq = np.inner(drVdA,np.inner(vcov,drVdA.T))
379    Vol = 1/np.sqrt(rVsq)
380    sigVol = Vol**3*np.sqrt(srcvlsq)/2.
381    R123 = Ax[0]*Ax[1]*Ax[2]
382    dsasdg = np.zeros((3,6))
383    dadg = np.zeros((6,6))
384    for i0 in range(3):         #0  1   2
385        i1 = (i0+1)%3           #1  2   0
386        i2 = (i1+1)%3           #2  0   1
387        i3 = 5-i2               #3  5   4
388        i4 = 5-i1               #4  3   5
389        i5 = 5-i0               #5  4   3
390        dsasdg[i0][i1] = 0.5*scot[i0]*scos[i0]/Ax[i1]
391        dsasdg[i0][i2] = 0.5*scot[i0]*scos[i0]/Ax[i2]
392        dsasdg[i0][i5] = -scot[i0]/np.sqrt(Ax[i1]*Ax[i2])
393        denmsq = Ax[i0]*(R123-Ax[i1]*Ax[i4]**2-Ax[i2]*Ax[i3]**2+(Ax[i4]*Ax[i3])**2)
394        denom = np.sqrt(denmsq)
395        dadg[i5][i0] = -Ax[i5]/denom-rcos[i0]/denmsq*(R123-0.5*Ax[i1]*Ax[i4]**2-0.5*Ax[i2]*Ax[i3]**2)
396        dadg[i5][i1] = -0.5*rcos[i0]/denmsq*(Ax[i0]**2*Ax[i2]-Ax[i0]*Ax[i4]**2)
397        dadg[i5][i2] = -0.5*rcos[i0]/denmsq*(Ax[i0]**2*Ax[i1]-Ax[i0]*Ax[i3]**2)
398        dadg[i5][i3] = Ax[i4]/denom+rcos[i0]/denmsq*(Ax[i0]*Ax[i2]*Ax[i3]-Ax[i3]*Ax[i4]**2)
399        dadg[i5][i4] = Ax[i3]/denom+rcos[i0]/denmsq*(Ax[i0]*Ax[i1]*Ax[i4]-Ax[i3]**2*Ax[i4])
400        dadg[i5][i5] = -Ax[i0]/denom
401    for i0 in range(3):
402        i1 = (i0+1)%3
403        i2 = (i1+1)%3
404        i3 = 5-i2
405        for ij in range(6):
406            dadg[i0][ij] = cell[i0]*(rcot[i2]*dadg[i3][ij]/rsin[i2]-dsasdg[i1][ij]/ssin[i1])
407            if ij == i0:
408                dadg[i0][ij] = dadg[i0][ij]-0.5*cell[i0]/Ax[i0]
409            dadg[i3][ij] = -dadg[i3][ij]*rsin[2-i0]*dpr
410    sigMat = np.inner(dadg,np.inner(vcov,dadg.T))
411    var = np.diag(sigMat)
412    CS = np.where(var>0.,np.sqrt(var),0.)
413    cellSig = [CS[0],CS[1],CS[2],CS[5],CS[4],CS[3],sigVol]  #exchange sig(alp) & sig(gam) to get in right order
414    return cellSig           
415   
416def SetPhaseData(parmDict,sigDict,Phases,covData):
417   
418    def PrintAtomsAndSig(General,Atoms,atomsSig):
419        print '\n Atoms:'
420        line = '   name      x         y         z      frac   Uiso     U11     U22     U33     U12     U13     U23'
421        if General['Type'] == 'magnetic':
422            line += '   Mx     My     Mz'
423        elif General['Type'] == 'macromolecular':
424            line = ' res no  residue  chain '+line
425        print line
426        if General['Type'] == 'nuclear':
427            print 135*'-'
428            fmt = {0:'%7s',1:'%7s',3:'%10.5f',4:'%10.5f',5:'%10.5f',6:'%8.3f',10:'%8.5f',
429                11:'%8.5f',12:'%8.5f',13:'%8.5f',14:'%8.5f',15:'%8.5f',16:'%8.5f'}
430            noFXsig = {3:[10*' ','%10s'],4:[10*' ','%10s'],5:[10*' ','%10s'],6:[8*' ','%8s']}
431            for i,at in enumerate(Atoms):
432                name = fmt[0]%(at[0])+fmt[1]%(at[1])+':'
433                valstr = ' values:'
434                sigstr = ' sig   :'
435                for ind in [3,4,5,6]:
436                    sigind = str(i)+':'+str(ind)
437                    valstr += fmt[ind]%(at[ind])                   
438                    if sigind in atomsSig:
439                        sigstr += fmt[ind]%(atomsSig[sigind])
440                    else:
441                        sigstr += noFXsig[ind][1]%(noFXsig[ind][0])
442                if at[9] == 'I':
443                    valstr += fmt[10]%(at[10])
444                    if str(i)+':10' in atomsSig:
445                        sigstr += fmt[10]%(atomsSig[str(i)+':10'])
446                    else:
447                        sigstr += 8*' '
448                else:
449                    valstr += 8*' '
450                    sigstr += 8*' '
451                    for ind in [11,12,13,14,15,16]:
452                        sigind = str(i)+':'+str(ind)
453                        valstr += fmt[ind]%(at[ind])
454                        if sigind in atomsSig:                       
455                            sigstr += fmt[ind]%(atomsSig[sigind])
456                        else:
457                            sigstr += 8*' '
458                print name
459                print valstr
460                print sigstr
461               
462    def PrintSHtextureAndSig(textureData,SHtextureSig):
463        print '\n Spherical harmonics texture: Order:' + str(textureData['Order'])
464        names = ['omega','chi','phi']
465        namstr = '  names :'
466        ptstr =  '  values:'
467        sigstr = '  esds  :'
468        for name in names:
469            namstr += '%12s'%(name)
470            ptstr += '%12.3f'%(textureData['Sample '+name][1])
471            if 'Sample '+name in SHtextureSig:
472                sigstr += '%12.3f'%(SHtextureSig['Sample '+name])
473            else:
474                sigstr += 12*' '
475        print namstr
476        print ptstr
477        print sigstr
478        print '\n Texture coefficients:'
479        namstr = '  names :'
480        ptstr =  '  values:'
481        sigstr = '  esds  :'
482        SHcoeff = textureData['SH Coeff'][1]
483        for name in SHcoeff:
484            namstr += '%12s'%(name)
485            ptstr += '%12.3f'%(SHcoeff[name])
486            if name in SHtextureSig:
487                sigstr += '%12.3f'%(SHtextureSig[name])
488            else:
489                sigstr += 12*' '
490        print namstr
491        print ptstr
492        print sigstr
493       
494           
495    print '\n Phases:'
496    for phase in Phases:
497        print ' Result for phase: ',phase
498        Phase = Phases[phase]
499        General = Phase['General']
500        SGData = General['SGData']
501        Atoms = Phase['Atoms']
502        cell = General['Cell']
503        pId = Phase['pId']
504        pfx = str(pId)+'::'
505        if cell[0]:
506            A,sigA = cellFill(pfx,SGData,parmDict,sigDict)
507            cellSig = getCellEsd(pfx,SGData,A,covData)  #includes sigVol
508            print ' Reciprocal metric tensor: '
509            ptfmt = "%15.9f"
510            names = ['A11','A22','A33','A12','A13','A23']
511            namstr = '  names :'
512            ptstr =  '  values:'
513            sigstr = '  esds  :'
514            for name,a,siga in zip(names,A,sigA):
515                namstr += '%15s'%(name)
516                ptstr += ptfmt%(a)
517                if siga:
518                    sigstr += ptfmt%(siga)
519                else:
520                    sigstr += 15*' '
521            print namstr
522            print ptstr
523            print sigstr
524            cell[1:7] = G2lat.A2cell(A)
525            cell[7] = G2lat.calc_V(A)
526            print ' New unit cell:'
527            ptfmt = ["%12.6f","%12.6f","%12.6f","%12.4f","%12.4f","%12.4f","%12.3f"]
528            names = ['a','b','c','alpha','beta','gamma','Volume']
529            namstr = '  names :'
530            ptstr =  '  values:'
531            sigstr = '  esds  :'
532            for name,fmt,a,siga in zip(names,ptfmt,cell[1:8],cellSig):
533                namstr += '%12s'%(name)
534                ptstr += fmt%(a)
535                if siga:
536                    sigstr += fmt%(siga)
537                else:
538                    sigstr += 12*' '
539            print namstr
540            print ptstr
541            print sigstr
542           
543        if 'Pawley' in Phase['General']['Type']:
544            pawleyRef = Phase['Pawley ref']
545            for i,refl in enumerate(pawleyRef):
546                key = pfx+'PWLref:'+str(i)
547                refl[6] = abs(parmDict[key])        #suppress negative Fsq
548                if key in sigDict:
549                    refl[7] = sigDict[key]
550                else:
551                    refl[7] = 0
552        else:
553            atomsSig = {}
554            if General['Type'] == 'nuclear':
555                for i,at in enumerate(Atoms):
556                    names = {3:pfx+'Ax:'+str(i),4:pfx+'Ay:'+str(i),5:pfx+'Az:'+str(i),6:pfx+'Afrac:'+str(i),
557                        10:pfx+'AUiso:'+str(i),11:pfx+'AU11:'+str(i),12:pfx+'AU22:'+str(i),13:pfx+'AU33:'+str(i),
558                        14:pfx+'AU12:'+str(i),15:pfx+'AU13:'+str(i),16:pfx+'AU23:'+str(i)}
559                    for ind in [3,4,5,6]:
560                        at[ind] = parmDict[names[ind]]
561                        if ind in [3,4,5]:
562                            name = names[ind].replace('A','dA')
563                        else:
564                            name = names[ind]
565                        if name in sigDict:
566                            atomsSig[str(i)+':'+str(ind)] = sigDict[name]
567                    if at[9] == 'I':
568                        at[10] = parmDict[names[10]]
569                        if names[10] in sigDict:
570                            atomsSig[str(i)+':10'] = sigDict[names[10]]
571                    else:
572                        for ind in [11,12,13,14,15,16]:
573                            at[ind] = parmDict[names[ind]]
574                            if names[ind] in sigDict:
575                                atomsSig[str(i)+':'+str(ind)] = sigDict[names[ind]]
576            PrintAtomsAndSig(General,Atoms,atomsSig)
577       
578        if 'SH Texture' in General:
579            textureData = General['SH Texture']   
580            if textureData['Order']:
581                SHtextureSig = {}
582                for name in ['omega','chi','phi']:
583                    aname = pfx+'SH '+name
584                    textureData['Sample '+name][1] = parmDict[aname]
585                    if aname in sigDict:
586                        SHtextureSig['Sample '+name] = sigDict[aname]
587                for name in textureData['SH Coeff'][1]:
588                    aname = pfx+name
589                    textureData['SH Coeff'][1][name] = parmDict[aname]
590                    if aname in sigDict:
591                        SHtextureSig[name] = sigDict[aname]
592                PrintSHtextureAndSig(textureData,SHtextureSig)
593
594def GetHistogramPhaseData(Phases,Histograms,Print=True):
595   
596    def PrintSize(hapData):
597        if hapData[0] in ['isotropic','uniaxial']:
598            line = '\n Size model    : %9s'%(hapData[0])
599            line += ' equatorial:'+'%12.3f'%(hapData[1][0])+' Refine? '+str(hapData[2][0])
600            if hapData[0] == 'uniaxial':
601                line += ' axial:'+'%12.3f'%(hapData[1][1])+' Refine? '+str(hapData[2][1])
602            print line
603        else:
604            print '\n Size model    : %s'%(hapData[0])
605            Snames = ['S11','S22','S33','S12','S13','S23']
606            ptlbls = ' names :'
607            ptstr =  ' values:'
608            varstr = ' refine:'
609            for i,name in enumerate(Snames):
610                ptlbls += '%12s' % (name)
611                ptstr += '%12.6f' % (hapData[4][i])
612                varstr += '%12s' % (str(hapData[5][i]))
613            print ptlbls
614            print ptstr
615            print varstr
616       
617    def PrintMuStrain(hapData,SGData):
618        if hapData[0] in ['isotropic','uniaxial']:
619            line = '\n Mustrain model: %9s'%(hapData[0])
620            line += ' equatorial:'+'%12.1f'%(hapData[1][0])+' Refine? '+str(hapData[2][0])
621            if hapData[0] == 'uniaxial':
622                line += ' axial:'+'%12.1f'%(hapData[1][1])+' Refine? '+str(hapData[2][1])
623            print line
624        else:
625            print '\n Mustrain model: %s'%(hapData[0])
626            Snames = G2spc.MustrainNames(SGData)
627            ptlbls = ' names :'
628            ptstr =  ' values:'
629            varstr = ' refine:'
630            for i,name in enumerate(Snames):
631                ptlbls += '%12s' % (name)
632                ptstr += '%12.6f' % (hapData[4][i])
633                varstr += '%12s' % (str(hapData[5][i]))
634            print ptlbls
635            print ptstr
636            print varstr
637
638    def PrintHStrain(hapData,SGData):
639        print '\n Hydrostatic/elastic strain: '
640        Hsnames = G2spc.HStrainNames(SGData)
641        ptlbls = ' names :'
642        ptstr =  ' values:'
643        varstr = ' refine:'
644        for i,name in enumerate(Hsnames):
645            ptlbls += '%12s' % (name)
646            ptstr += '%12.6f' % (hapData[0][i])
647            varstr += '%12s' % (str(hapData[1][i]))
648        print ptlbls
649        print ptstr
650        print varstr
651
652    def PrintSHPO(hapData):
653        print '\n Spherical harmonics preferred orientation: Order:' + \
654            str(hapData[4])+' Refine? '+str(hapData[2])
655        ptlbls = ' names :'
656        ptstr =  ' values:'
657        for item in hapData[5]:
658            ptlbls += '%12s'%(item)
659            ptstr += '%12.3f'%(hapData[5][item]) 
660        print ptlbls
661        print ptstr
662   
663    hapDict = {}
664    hapVary = []
665    controlDict = {}
666    poType = {}
667    poAxes = {}
668    spAxes = {}
669    spType = {}
670   
671    for phase in Phases:
672        HistoPhase = Phases[phase]['Histograms']
673        SGData = Phases[phase]['General']['SGData']
674        cell = Phases[phase]['General']['Cell'][1:7]
675        A = G2lat.cell2A(cell)
676        pId = Phases[phase]['pId']
677        histoList = HistoPhase.keys()
678        histoList.sort()
679        for histogram in histoList:
680            try:
681                Histogram = Histograms[histogram]
682            except KeyError:                       
683                #skip if histogram not included e.g. in a sequential refinement
684                continue
685            hapData = HistoPhase[histogram]
686            hId = Histogram['hId']
687            limits = Histogram['Limits'][1]
688            inst = Histogram['Instrument Parameters']
689            inst = dict(zip(inst[3],inst[1]))
690            Zero = inst['Zero']
691            if 'C' in inst['Type']:
692                try:
693                    wave = inst['Lam']
694                except KeyError:
695                    wave = inst['Lam1']
696                dmin = wave/(2.0*sind(limits[1]/2.0))
697            pfx = str(pId)+':'+str(hId)+':'
698            for item in ['Scale','Extinction']:
699                hapDict[pfx+item] = hapData[item][0]
700                if hapData[item][1]:
701                    hapVary.append(pfx+item)
702            names = G2spc.HStrainNames(SGData)
703            for i,name in enumerate(names):
704                hapDict[pfx+name] = hapData['HStrain'][0][i]
705                if hapData['HStrain'][1][i]:
706                    hapVary.append(pfx+name)
707            controlDict[pfx+'poType'] = hapData['Pref.Ori.'][0]
708            if hapData['Pref.Ori.'][0] == 'MD':
709                hapDict[pfx+'MD'] = hapData['Pref.Ori.'][1]
710                controlDict[pfx+'MDAxis'] = hapData['Pref.Ori.'][3]
711                if hapData['Pref.Ori.'][2]:
712                    hapVary.append(pfx+'MD')
713            else:                           #'SH' spherical harmonics
714                controlDict[pfx+'SHord'] = hapData['Pref.Ori.'][4]
715                controlDict[pfx+'SHncof'] = len(hapData['Pref.Ori.'][5])
716                for item in hapData['Pref.Ori.'][5]:
717                    hapDict[pfx+item] = hapData['Pref.Ori.'][5][item]
718                    if hapData['Pref.Ori.'][2]:
719                        hapVary.append(pfx+item)
720            for item in ['Mustrain','Size']:
721                controlDict[pfx+item+'Type'] = hapData[item][0]
722                if hapData[item][0] in ['isotropic','uniaxial']:
723                    hapDict[pfx+item+':i'] = hapData[item][1][0]
724                    if hapData[item][2][0]:
725                        hapVary.append(pfx+item+':i')
726                    if hapData[item][0] == 'uniaxial':
727                        controlDict[pfx+item+'Axis'] = hapData[item][3]
728                        hapDict[pfx+item+':a'] = hapData[item][1][1]
729                        if hapData[item][2][1]:
730                            hapVary.append(pfx+item+':a')
731                else:       #generalized for mustrain or ellipsoidal for size
732                    if item == 'Mustrain':
733                        names = G2spc.MustrainNames(SGData)
734                        pwrs = []
735                        for name in names:
736                            h,k,l = name[1:]
737                            pwrs.append([int(h),int(k),int(l)])
738                        controlDict[pfx+'MuPwrs'] = pwrs
739                    for i in range(len(hapData[item][4])):
740                        sfx = ':'+str(i)
741                        hapDict[pfx+item+sfx] = hapData[item][4][i]
742                        if hapData[item][5][i]:
743                            hapVary.append(pfx+item+sfx)
744                           
745            if Print: 
746                print '\n Phase: ',phase,' in histogram: ',histogram
747                print 135*'-'
748                print ' Phase fraction  : %10.4f'%(hapData['Scale'][0]),' Refine?',hapData['Scale'][1]
749                print ' Extinction coeff: %10.4f'%(hapData['Extinction'][0]),' Refine?',hapData['Extinction'][1]
750                if hapData['Pref.Ori.'][0] == 'MD':
751                    Ax = hapData['Pref.Ori.'][3]
752                    print ' March-Dollase PO: %10.4f'%(hapData['Pref.Ori.'][1]),' Refine?',hapData['Pref.Ori.'][2], \
753                        ' Axis: %d %d %d'%(Ax[0],Ax[1],Ax[2])
754                else: #'SH' for spherical harmonics
755                    PrintSHPO(hapData['Pref.Ori.'])
756                PrintSize(hapData['Size'])
757                PrintMuStrain(hapData['Mustrain'],SGData)
758                PrintHStrain(hapData['HStrain'],SGData)
759            HKLd = np.array(G2lat.GenHLaue(dmin,SGData,A))
760            refList = []
761            for h,k,l,d in HKLd:
762                ext,mul,Uniq,phi = G2spc.GenHKLf([h,k,l],SGData)
763                if ext:
764                    continue
765                if 'C' in inst['Type']:
766                    pos = 2.0*asind(wave/(2.0*d))+Zero
767                    if limits[0] < pos < limits[1]:
768                        refList.append([h,k,l,mul,d,pos,0.0,0.0,0.0,0.0,0.0,Uniq,phi,0.0])
769                else:
770                    raise ValueError 
771            Histogram['Reflection Lists'][phase] = refList
772    return hapVary,hapDict,controlDict
773   
774def SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms,Print=True):
775   
776    def PrintSizeAndSig(hapData,sizeSig):
777        line = '\n Size model:     %9s'%(hapData[0])
778        refine = False
779        if hapData[0] in ['isotropic','uniaxial']:
780            line += ' equatorial:%12.3f'%(hapData[1][0])
781            if sizeSig[0][0]:
782                line += ', sig: %8.3f'%(sizeSig[0][0])
783                refine = True
784            if hapData[0] == 'uniaxial':
785                line += ' axial:%12.3f'%(hapData[1][1])
786                if sizeSig[0][1]:
787                    refine = True
788                    line += ', sig: %8.3f'%(sizeSig[0][1])
789            if refine:
790                print line
791        else:
792            Snames = ['S11','S22','S33','S12','S13','S23']
793            ptlbls = ' name  :'
794            ptstr =  ' value :'
795            sigstr = ' sig   :'
796            for i,name in enumerate(Snames):
797                ptlbls += '%12s' % (name)
798                ptstr += '%12.6f' % (hapData[4][i])
799                if sizeSig[1][i]:
800                    refine = True
801                    sigstr += '%12.6f' % (sizeSig[1][i])
802                else:
803                    sigstr += 12*' '
804            if refine:
805                print line
806                print ptlbls
807                print ptstr
808                print sigstr
809       
810    def PrintMuStrainAndSig(hapData,mustrainSig,SGData):
811        line = '\n Mustrain model: %9s'%(hapData[0])
812        refine = False
813        if hapData[0] in ['isotropic','uniaxial']:
814            line += ' equatorial:%12.1f'%(hapData[1][0])
815            if mustrainSig[0][0]:
816                line += ', sig: %8.1f'%(mustrainSig[0][0])
817                refine = True
818            if hapData[0] == 'uniaxial':
819                line += ' axial:%12.1f'%(hapData[1][1])
820                if mustrainSig[0][1]:
821                     line += ', sig: %8.1f'%(mustrainSig[0][1])
822            if refine:
823                print line
824        else:
825            Snames = G2spc.MustrainNames(SGData)
826            ptlbls = ' name  :'
827            ptstr =  ' value :'
828            sigstr = ' sig   :'
829            for i,name in enumerate(Snames):
830                ptlbls += '%12s' % (name)
831                ptstr += '%12.6f' % (hapData[4][i])
832                if mustrainSig[1][i]:
833                    refine = True
834                    sigstr += '%12.6f' % (mustrainSig[1][i])
835                else:
836                    sigstr += 12*' '
837            if refine:
838                print line
839                print ptlbls
840                print ptstr
841                print sigstr
842           
843    def PrintHStrainAndSig(hapData,strainSig,SGData):
844        Hsnames = G2spc.HStrainNames(SGData)
845        ptlbls = ' name  :'
846        ptstr =  ' value :'
847        sigstr = ' sig   :'
848        refine = False
849        for i,name in enumerate(Hsnames):
850            ptlbls += '%12s' % (name)
851            ptstr += '%12.6g' % (hapData[0][i])
852            if name in strainSig:
853                refine = True
854                sigstr += '%12.6g' % (strainSig[name])
855            else:
856                sigstr += 12*' '
857        if refine:
858            print '\n Hydrostatic/elastic strain: '
859            print ptlbls
860            print ptstr
861            print sigstr
862       
863    def PrintSHPOAndSig(hapData,POsig):
864        print '\n Spherical harmonics preferred orientation: Order:'+str(hapData[4])
865        ptlbls = ' names :'
866        ptstr =  ' values:'
867        sigstr = ' sig   :'
868        for item in hapData[5]:
869            ptlbls += '%12s'%(item)
870            ptstr += '%12.3f'%(hapData[5][item])
871            if item in POsig:
872                sigstr += '%12.3f'%(POsig[item])
873            else:
874                sigstr += 12*' ' 
875        print ptlbls
876        print ptstr
877        print sigstr
878   
879    for phase in Phases:
880        HistoPhase = Phases[phase]['Histograms']
881        SGData = Phases[phase]['General']['SGData']
882        pId = Phases[phase]['pId']
883        histoList = HistoPhase.keys()
884        histoList.sort()
885        for histogram in histoList:
886            try:
887                Histogram = Histograms[histogram]
888            except KeyError:                       
889                #skip if histogram not included e.g. in a sequential refinement
890                continue
891            print '\n Phase: ',phase,' in histogram: ',histogram
892            print 130*'-'
893            hapData = HistoPhase[histogram]
894            hId = Histogram['hId']
895            pfx = str(pId)+':'+str(hId)+':'
896            print ' Final refinement RF, RF^2 = %.2f%%, %.2f%% on %d reflections'   \
897                %(Histogram[pfx+'Rf'],Histogram[pfx+'Rf^2'],Histogram[pfx+'Nref'])
898           
899            PhFrExtPOSig = {}
900            for item in ['Scale','Extinction']:
901                hapData[item][0] = parmDict[pfx+item]
902                if pfx+item in sigDict:
903                    PhFrExtPOSig[item] = sigDict[pfx+item]
904            if hapData['Pref.Ori.'][0] == 'MD':
905                hapData['Pref.Ori.'][1] = parmDict[pfx+'MD']
906                if pfx+'MD' in sigDict:
907                    PhFrExtPOSig['MD'] = sigDict[pfx+'MD']
908            else:                           #'SH' spherical harmonics
909                for item in hapData['Pref.Ori.'][5]:
910                    hapData['Pref.Ori.'][5][item] = parmDict[pfx+item]
911                    if pfx+item in sigDict:
912                        PhFrExtPOSig[item] = sigDict[pfx+item]
913            if Print:
914                if 'Scale' in PhFrExtPOSig:
915                    print ' Phase fraction  : %10.4f, sig %10.4f'%(hapData['Scale'][0],PhFrExtPOSig['Scale'])
916                if 'Extinction' in PhFrExtPOSig:
917                    print ' Extinction coeff: %10.4f, sig %10.4f'%(hapData['Extinction'][0],PhFrExtPOSig['Extinction'])
918                if hapData['Pref.Ori.'][0] == 'MD':
919                    if 'MD' in PhFrExtPOSig:
920                        print ' March-Dollase PO: %10.4f, sig %10.4f'%(hapData['Pref.Ori.'][1],PhFrExtPOSig['MD'])
921                else:
922                    PrintSHPOAndSig(hapData['Pref.Ori.'],PhFrExtPOSig)
923            SizeMuStrSig = {'Mustrain':[[0,0],[0 for i in range(len(hapData['Mustrain'][4]))]],
924                'Size':[[0,0],[0 for i in range(len(hapData['Size'][4]))]],
925                'HStrain':{}}                 
926            for item in ['Mustrain','Size']:
927                if hapData[item][0] in ['isotropic','uniaxial']:                   
928                    hapData[item][1][0] = parmDict[pfx+item+':i']
929                    if item == 'Size':
930                        hapData[item][1][0] = min(10.,max(0.001,hapData[item][1][0]))
931                    if pfx+item+':i' in sigDict: 
932                        SizeMuStrSig[item][0][0] = sigDict[pfx+item+':i']
933                    if hapData[item][0] == 'uniaxial':
934                        hapData[item][1][1] = parmDict[pfx+item+':a']
935                        if item == 'Size':
936                            hapData[item][1][1] = min(10.,max(0.001,hapData[item][1][1]))                       
937                        if pfx+item+':a' in sigDict:
938                            SizeMuStrSig[item][0][1] = sigDict[pfx+item+':a']
939                else:       #generalized for mustrain or ellipsoidal for size
940                    for i in range(len(hapData[item][4])):
941                        sfx = ':'+str(i)
942                        hapData[item][4][i] = parmDict[pfx+item+sfx]
943                        if pfx+item+sfx in sigDict:
944                            SizeMuStrSig[item][1][i] = sigDict[pfx+item+sfx]
945            names = G2spc.HStrainNames(SGData)
946            for i,name in enumerate(names):
947                hapData['HStrain'][0][i] = parmDict[pfx+name]
948                if pfx+name in sigDict:
949                    SizeMuStrSig['HStrain'][name] = sigDict[pfx+name]
950            if Print:
951                PrintSizeAndSig(hapData['Size'],SizeMuStrSig['Size'])
952                PrintMuStrainAndSig(hapData['Mustrain'],SizeMuStrSig['Mustrain'],SGData)
953                PrintHStrainAndSig(hapData['HStrain'],SizeMuStrSig['HStrain'],SGData)
954   
955def GetHistogramData(Histograms,Print=True):
956   
957    def GetBackgroundParms(hId,Background):
958        Back = Background[0]
959        Debye = Background[1]
960        bakType,bakFlag = Back[:2]
961        backVals = Back[3:]
962        backNames = [':'+str(hId)+':Back:'+str(i) for i in range(len(backVals))]
963        backDict = dict(zip(backNames,backVals))
964        backVary = []
965        if bakFlag:
966            backVary = backNames
967        backDict[':'+str(hId)+':nDebye'] = Debye['nDebye']
968        debyeDict = {}
969        debyeList = []
970        for i in range(Debye['nDebye']):
971            debyeNames = [':'+str(hId)+':DebyeA:'+str(i),':'+str(hId)+':DebyeR:'+str(i),':'+str(hId)+':DebyeU:'+str(i)]
972            debyeDict.update(dict(zip(debyeNames,Debye['debyeTerms'][i][::2])))
973            debyeList += zip(debyeNames,Debye['debyeTerms'][i][1::2])
974        debyeVary = []
975        for item in debyeList:
976            if item[1]:
977                debyeVary.append(item[0])
978        backDict.update(debyeDict)
979        backVary += debyeVary   
980        return bakType,backDict,backVary           
981       
982    def GetInstParms(hId,Inst):
983        insVals,insFlags,insNames = Inst[1:4]
984        dataType = insVals[0]
985        instDict = {}
986        insVary = []
987        pfx = ':'+str(hId)+':'
988        for i,flag in enumerate(insFlags):
989            insName = pfx+insNames[i]
990            instDict[insName] = insVals[i]
991            if flag:
992                insVary.append(insName)
993        instDict[pfx+'X'] = max(instDict[pfx+'X'],0.001)
994        instDict[pfx+'Y'] = max(instDict[pfx+'Y'],0.001)
995        instDict[pfx+'SH/L'] = max(instDict[pfx+'SH/L'],0.0005)
996        return dataType,instDict,insVary
997       
998    def GetSampleParms(hId,Sample):
999        sampVary = []
1000        hfx = ':'+str(hId)+':'       
1001        sampDict = {hfx+'Gonio. radius':Sample['Gonio. radius'],hfx+'Omega':Sample['Omega'],
1002            hfx+'Chi':Sample['Chi'],hfx+'Phi':Sample['Phi']}
1003        Type = Sample['Type']
1004        if 'Bragg' in Type:             #Bragg-Brentano
1005            for item in ['Scale','Shift','Transparency']:       #surface roughness?, diffuse scattering?
1006                sampDict[hfx+item] = Sample[item][0]
1007                if Sample[item][1]:
1008                    sampVary.append(hfx+item)
1009        elif 'Debye' in Type:        #Debye-Scherrer
1010            for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
1011                sampDict[hfx+item] = Sample[item][0]
1012                if Sample[item][1]:
1013                    sampVary.append(hfx+item)
1014        return Type,sampDict,sampVary
1015       
1016    def PrintBackground(Background):
1017        Back = Background[0]
1018        Debye = Background[1]
1019        print '\n Background function: ',Back[0],' Refine?',bool(Back[1])
1020        line = ' Coefficients: '
1021        for i,back in enumerate(Back[3:]):
1022            line += '%10.3f'%(back)
1023            if i and not i%10:
1024                line += '\n'+15*' '
1025        print line
1026        if Debye['nDebye']:
1027            print '\n Debye diffuse scattering coefficients'
1028            parms = ['DebyeA','DebyeR','DebyeU']
1029            line = ' names :'
1030            for parm in parms:
1031                line += '%16s'%(parm)
1032            print line
1033            for j,term in enumerate(Debye['debyeTerms']):
1034                line = ' term'+'%2d'%(j)+':'
1035                for i in range(3):
1036                    line += '%10.4g %5s'%(term[2*i],bool(term[2*i+1]))                   
1037                print line
1038       
1039    def PrintInstParms(Inst):
1040        print '\n Instrument Parameters:'
1041        ptlbls = ' name  :'
1042        ptstr =  ' value :'
1043        varstr = ' refine:'
1044        instNames = Inst[3][1:]
1045        for i,name in enumerate(instNames):
1046            ptlbls += '%12s' % (name)
1047            ptstr += '%12.6f' % (Inst[1][i+1])
1048            if name in ['Lam1','Lam2','Azimuth']:
1049                varstr += 12*' '
1050            else:
1051                varstr += '%12s' % (str(bool(Inst[2][i+1])))
1052        print ptlbls
1053        print ptstr
1054        print varstr
1055       
1056    def PrintSampleParms(Sample):
1057        print '\n Sample Parameters:'
1058        print ' Goniometer omega = %.2f, chi = %.2f, phi = %.2f'% \
1059            (Sample['Omega'],Sample['Chi'],Sample['Phi'])
1060        ptlbls = ' name  :'
1061        ptstr =  ' value :'
1062        varstr = ' refine:'
1063        if 'Bragg' in Sample['Type']:
1064            for item in ['Scale','Shift','Transparency']:
1065                ptlbls += '%14s'%(item)
1066                ptstr += '%14.4f'%(Sample[item][0])
1067                varstr += '%14s'%(str(bool(Sample[item][1])))
1068           
1069        elif 'Debye' in Type:        #Debye-Scherrer
1070            for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
1071                ptlbls += '%14s'%(item)
1072                ptstr += '%14.4f'%(Sample[item][0])
1073                varstr += '%14s'%(str(bool(Sample[item][1])))
1074
1075        print ptlbls
1076        print ptstr
1077        print varstr
1078       
1079
1080    histDict = {}
1081    histVary = []
1082    controlDict = {}
1083    histoList = Histograms.keys()
1084    histoList.sort()
1085    for histogram in histoList:
1086        Histogram = Histograms[histogram]
1087        hId = Histogram['hId']
1088        pfx = ':'+str(hId)+':'
1089        controlDict[pfx+'Limits'] = Histogram['Limits'][1]
1090       
1091        Background = Histogram['Background']
1092        Type,bakDict,bakVary = GetBackgroundParms(hId,Background)
1093        controlDict[pfx+'bakType'] = Type
1094        histDict.update(bakDict)
1095        histVary += bakVary
1096       
1097        Inst = Histogram['Instrument Parameters']
1098        Type,instDict,insVary = GetInstParms(hId,Inst)
1099        controlDict[pfx+'histType'] = Type
1100        if pfx+'Lam1' in instDict:
1101            controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam1']
1102        else:
1103            controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam']           
1104        histDict.update(instDict)
1105        histVary += insVary
1106       
1107        Sample = Histogram['Sample Parameters']
1108        Type,sampDict,sampVary = GetSampleParms(hId,Sample)
1109        controlDict[pfx+'instType'] = Type
1110        histDict.update(sampDict)
1111        histVary += sampVary
1112
1113        if Print: 
1114            print '\n Histogram: ',histogram,' histogram Id: ',hId
1115            print 135*'-'
1116            Units = {'C':' deg','T':' msec'}
1117            units = Units[controlDict[pfx+'histType'][2]]
1118            Limits = controlDict[pfx+'Limits']
1119            print ' Instrument type: ',Sample['Type']
1120            print ' Histogram limits: %8.2f%s to %8.2f%s'%(Limits[0],units,Limits[1],units)     
1121            PrintSampleParms(Sample)
1122            PrintInstParms(Inst)
1123            PrintBackground(Background)
1124       
1125    return histVary,histDict,controlDict
1126   
1127def SetHistogramData(parmDict,sigDict,Histograms,Print=True):
1128   
1129    def SetBackgroundParms(pfx,Background,parmDict,sigDict):
1130        Back = Background[0]
1131        Debye = Background[1]
1132        lenBack = len(Back[3:])
1133        backSig = [0 for i in range(lenBack+3*Debye['nDebye'])]
1134        for i in range(lenBack):
1135            Back[3+i] = parmDict[pfx+'Back:'+str(i)]
1136            if pfx+'Back:'+str(i) in sigDict:
1137                backSig[i] = sigDict[pfx+'Back:'+str(i)]
1138        if Debye['nDebye']:
1139            for i in range(Debye['nDebye']):
1140                names = [pfx+'DebyeA:'+str(i),pfx+'DebyeR:'+str(i),pfx+'DebyeU:'+str(i)]
1141                for j,name in enumerate(names):
1142                    Debye['debyeTerms'][i][2*j] = parmDict[name]
1143                    if name in sigDict:
1144                        backSig[lenBack+3*i+j] = sigDict[name]           
1145        return backSig
1146       
1147    def SetInstParms(pfx,Inst,parmDict,sigDict):
1148        insVals,insFlags,insNames = Inst[1:4]
1149        instSig = [0 for i in range(len(insVals))]
1150        for i,flag in enumerate(insFlags):
1151            insName = pfx+insNames[i]
1152            insVals[i] = parmDict[insName]
1153            if insName in sigDict:
1154                instSig[i] = sigDict[insName]
1155        return instSig
1156       
1157    def SetSampleParms(pfx,Sample,parmDict,sigDict):
1158        if 'Bragg' in Sample['Type']:             #Bragg-Brentano
1159            sampSig = [0 for i in range(3)]
1160            for i,item in enumerate(['Scale','Shift','Transparency']):       #surface roughness?, diffuse scattering?
1161                Sample[item][0] = parmDict[pfx+item]
1162                if pfx+item in sigDict:
1163                    sampSig[i] = sigDict[pfx+item]
1164        elif 'Debye' in Sample['Type']:        #Debye-Scherrer
1165            sampSig = [0 for i in range(4)]
1166            for i,item in enumerate(['Scale','Absorption','DisplaceX','DisplaceY']):
1167                Sample[item][0] = parmDict[pfx+item]
1168                if pfx+item in sigDict:
1169                    sampSig[i] = sigDict[pfx+item]
1170        return sampSig
1171       
1172    def PrintBackgroundSig(Background,backSig):
1173        Back = Background[0]
1174        Debye = Background[1]
1175        lenBack = len(Back[3:])
1176        valstr = ' value : '
1177        sigstr = ' sig   : '
1178        refine = False
1179        for i,back in enumerate(Back[3:]):
1180            valstr += '%10.4g'%(back)
1181            if Back[1]:
1182                refine = True
1183                sigstr += '%10.4g'%(backSig[i])
1184            else:
1185                sigstr += 10*' '
1186        if refine:
1187            print '\n Background function: ',Back[0]
1188            print valstr
1189            print sigstr
1190        if Debye['nDebye']:
1191            ifAny = False
1192            ptfmt = "%12.5f"
1193            names =  ' names :'
1194            ptstr =  ' values:'
1195            sigstr = ' esds  :'
1196            for item in sigDict:
1197                if 'Debye' in item:
1198                    ifAny = True
1199                    names += '%12s'%(item)
1200                    ptstr += ptfmt%(parmDict[item])
1201                    sigstr += ptfmt%(sigDict[item])
1202            if ifAny:
1203                print '\n Debye diffuse scattering coefficients'
1204                print names
1205                print ptstr
1206                print sigstr
1207       
1208    def PrintInstParmsSig(Inst,instSig):
1209        ptlbls = ' names :'
1210        ptstr =  ' value :'
1211        sigstr = ' sig   :'
1212        instNames = Inst[3][1:]
1213        refine = False
1214        for i,name in enumerate(instNames):
1215            ptlbls += '%12s' % (name)
1216            ptstr += '%12.6f' % (Inst[1][i+1])
1217            if instSig[i+1]:
1218                refine = True
1219                sigstr += '%12.6f' % (instSig[i+1])
1220            else:
1221                sigstr += 12*' '
1222        if refine:
1223            print '\n Instrument Parameters:'
1224            print ptlbls
1225            print ptstr
1226            print sigstr
1227       
1228    def PrintSampleParmsSig(Sample,sampleSig):
1229        ptlbls = ' names :'
1230        ptstr =  ' values:'
1231        sigstr = ' sig   :'
1232        refine = False
1233        if 'Bragg' in Sample['Type']:
1234            for i,item in enumerate(['Scale','Shift','Transparency']):
1235                ptlbls += '%14s'%(item)
1236                ptstr += '%14.4f'%(Sample[item][0])
1237                if sampleSig[i]:
1238                    refine = True
1239                    sigstr += '%14.4f'%(sampleSig[i])
1240                else:
1241                    sigstr += 14*' '
1242           
1243        elif 'Debye' in Sample['Type']:        #Debye-Scherrer
1244            for i,item in enumerate(['Scale','Absorption','DisplaceX','DisplaceY']):
1245                ptlbls += '%14s'%(item)
1246                ptstr += '%14.4f'%(Sample[item][0])
1247                if sampleSig[i]:
1248                    refine = True
1249                    sigstr += '%14.4f'%(sampleSig[i])
1250                else:
1251                    sigstr += 14*' '
1252
1253        if refine:
1254            print '\n Sample Parameters:'
1255            print ptlbls
1256            print ptstr
1257            print sigstr
1258       
1259    histoList = Histograms.keys()
1260    histoList.sort()
1261    for histogram in histoList:
1262        if 'PWDR' in histogram:
1263            Histogram = Histograms[histogram]
1264            hId = Histogram['hId']
1265            pfx = ':'+str(hId)+':'
1266            Background = Histogram['Background']
1267            backSig = SetBackgroundParms(pfx,Background,parmDict,sigDict)
1268           
1269            Inst = Histogram['Instrument Parameters']
1270            instSig = SetInstParms(pfx,Inst,parmDict,sigDict)
1271       
1272            Sample = Histogram['Sample Parameters']
1273            sampSig = SetSampleParms(pfx,Sample,parmDict,sigDict)
1274
1275            print '\n Histogram: ',histogram,' histogram Id: ',hId
1276            print 135*'-'
1277            print ' Final refinement wRp = %.2f%% on %d observations in this histogram'%(Histogram['wRp'],Histogram['Nobs'])
1278            if Print:
1279                print ' Instrument type: ',Sample['Type']
1280                PrintSampleParmsSig(Sample,sampSig)
1281                PrintInstParmsSig(Inst,instSig)
1282                PrintBackgroundSig(Background,backSig)
1283
1284def GetAtomFXU(pfx,FFtables,BLtables,calcControls,parmDict):
1285    Natoms = calcControls['Natoms'][pfx]
1286    Tdata = Natoms*[' ',]
1287    Mdata = np.zeros(Natoms)
1288    IAdata = Natoms*[' ',]
1289    Fdata = np.zeros(Natoms)
1290    FFdata = []
1291    BLdata = []
1292    Xdata = np.zeros((3,Natoms))
1293    dXdata = np.zeros((3,Natoms))
1294    Uisodata = np.zeros(Natoms)
1295    Uijdata = np.zeros((6,Natoms))
1296    keys = {'Atype:':Tdata,'Amul:':Mdata,'Afrac:':Fdata,'AI/A:':IAdata,
1297        'dAx:':dXdata[0],'dAy:':dXdata[1],'dAz:':dXdata[2],
1298        'Ax:':Xdata[0],'Ay:':Xdata[1],'Az:':Xdata[2],'AUiso:':Uisodata,
1299        'AU11:':Uijdata[0],'AU22:':Uijdata[1],'AU33:':Uijdata[2],
1300        'AU12:':Uijdata[3],'AU13:':Uijdata[4],'AU23:':Uijdata[5]}
1301    for iatm in range(Natoms):
1302        for key in keys:
1303            parm = pfx+key+str(iatm)
1304            if parm in parmDict:
1305                keys[key][iatm] = parmDict[parm]
1306        FFdata.append(FFtables[Tdata[iatm]])
1307        BLdata.append(BLtables[Tdata[iatm]][1])
1308    return FFdata,BLdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata
1309   
1310def StructureFactor(refList,G,hfx,pfx,SGData,calcControls,parmDict):
1311    ''' Compute structure factors for all h,k,l for phase
1312    input:
1313        refList: [ref] where each ref = h,k,l,m,d,...,[equiv h,k,l],phase[equiv]
1314        G:      reciprocal metric tensor
1315        pfx:    phase id string
1316        SGData: space group info. dictionary output from SpcGroup
1317        calcControls:
1318        ParmDict:
1319    puts result F^2 in each ref[8] in refList
1320    '''       
1321    twopi = 2.0*np.pi
1322    twopisq = 2.0*np.pi**2
1323    ast = np.sqrt(np.diag(G))
1324    Mast = twopisq*np.multiply.outer(ast,ast)
1325    FFtables = calcControls['FFtables']
1326    BLtables = calcControls['BLtables']
1327    FFdata,BLdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,FFtables,BLtables,calcControls,parmDict)
1328    if 'N' in parmDict[hfx+'Type']:
1329        FP,FPP = G2el.BlenRes(BLdata,parmDict[hfx+'Lam'])
1330    else:
1331        FP = np.array([El[hfx+'FP'] for El in FFdata])
1332        FPP = np.array([El[hfx+'FPP'] for El in FFdata])
1333    maxPos = len(SGData['SGOps'])
1334    Uij = np.array(G2lat.U6toUij(Uijdata))
1335    bij = Mast*Uij.T
1336    for refl in refList:
1337        fbs = np.array([0,0])
1338        H = refl[:3]
1339        SQ = 1./(2.*refl[4])**2
1340        if 'N' in parmDict[hfx+'Type']:
1341            FF = np.array([El[1] for El in BLdata])
1342        else:       #'X'
1343            FF = np.array([G2el.ScatFac(El,SQ)[0] for El in FFdata])
1344        SQfactor = 4.0*SQ*twopisq
1345        Uniq = refl[11]
1346        phi = refl[12]
1347        phase = twopi*(np.inner(Uniq,(dXdata.T+Xdata.T))+phi[:,np.newaxis])
1348        sinp = np.sin(phase)
1349        cosp = np.cos(phase)
1350        occ = Mdata*Fdata/len(Uniq)
1351        biso = -SQfactor*Uisodata
1352        Tiso = np.where(biso<1.,np.exp(biso),1.0)
1353        HbH = np.array([-np.inner(h,np.inner(bij,h)) for h in Uniq])
1354        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
1355        Tcorr = Tiso*Tuij
1356        fa = np.array([(FF+FP)*occ*cosp*Tcorr,-FPP*occ*sinp*Tcorr])
1357        fas = np.sum(np.sum(fa,axis=1),axis=1)        #real
1358        if not SGData['SGInv']:
1359            fb = np.array([(FF+FP)*occ*sinp*Tcorr,FPP*occ*cosp*Tcorr])
1360            fbs = np.sum(np.sum(fb,axis=1),axis=1)
1361        fasq = fas**2
1362        fbsq = fbs**2        #imaginary
1363        refl[9] = np.sum(fasq)+np.sum(fbsq)
1364        refl[10] = atan2d(fbs[0],fas[0])
1365    return refList
1366   
1367def StructureFactorDerv(refList,G,hfx,pfx,SGData,calcControls,parmDict):
1368    twopi = 2.0*np.pi
1369    twopisq = 2.0*np.pi**2
1370    ast = np.sqrt(np.diag(G))
1371    Mast = twopisq*np.multiply.outer(ast,ast)
1372    FFtables = calcControls['FFtables']
1373    BLtables = calcControls['BLtables']
1374    FFdata,BLdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,FFtables,BLtables,calcControls,parmDict)
1375    if 'N' in parmDict[hfx+'Type']:
1376        FP = 0.
1377        FPP = 0.
1378    else:
1379        FP = np.array([El[hfx+'FP'] for El in FFdata])
1380        FPP = np.array([El[hfx+'FPP'] for El in FFdata])
1381    maxPos = len(SGData['SGOps'])       
1382    Uij = np.array(G2lat.U6toUij(Uijdata))
1383    bij = Mast*Uij.T
1384    dFdvDict = {}
1385    dFdfr = np.zeros((len(refList),len(Mdata)))
1386    dFdx = np.zeros((len(refList),len(Mdata),3))
1387    dFdui = np.zeros((len(refList),len(Mdata)))
1388    dFdua = np.zeros((len(refList),len(Mdata),6))
1389    for iref,refl in enumerate(refList):
1390        H = np.array(refl[:3])
1391        SQ = 1./(2.*refl[4])**2          # or (sin(theta)/lambda)**2
1392        if 'N' in parmDict[hfx+'Type']:
1393            FF = np.array([El[1] for El in BLdata])
1394        else:       #'X'
1395            FF = np.array([G2el.ScatFac(El,SQ)[0] for El in FFdata])
1396        SQfactor = 8.0*SQ*np.pi**2
1397        Uniq = refl[11]
1398        phi = refl[12]
1399        phase = twopi*(np.inner((dXdata.T+Xdata.T),Uniq)+phi[np.newaxis,:])
1400        sinp = np.sin(phase)
1401        cosp = np.cos(phase)
1402        occ = Mdata*Fdata/len(Uniq)
1403        biso = -SQfactor*Uisodata
1404        Tiso = np.where(biso<1.,np.exp(biso),1.0)
1405#        HbH = np.array([-np.inner(h,np.inner(bij,h)) for h in Uniq])
1406        HbH = -np.inner(H,np.inner(bij,H))
1407        Hij = np.array([Mast*np.multiply.outer(U,U) for U in Uniq])
1408        Hij = np.array([G2lat.UijtoU6(Uij) for Uij in Hij])
1409        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
1410        Tcorr = Tiso*Tuij
1411        fot = (FF+FP)*occ*Tcorr
1412        fotp = FPP*occ*Tcorr
1413        fa = np.array([fot[:,np.newaxis]*cosp,fotp[:,np.newaxis]*cosp])       #non positions
1414        fb = np.array([fot[:,np.newaxis]*sinp,-fotp[:,np.newaxis]*sinp])
1415       
1416        fas = np.sum(np.sum(fa,axis=1),axis=1)
1417        fbs = np.sum(np.sum(fb,axis=1),axis=1)
1418        fax = np.array([-fot[:,np.newaxis]*sinp,-fotp[:,np.newaxis]*sinp])   #positions
1419        fbx = np.array([fot[:,np.newaxis]*cosp,-fot[:,np.newaxis]*cosp])
1420        #sum below is over Uniq
1421        dfadfr = np.sum(fa/occ[:,np.newaxis],axis=2)
1422        dfadx = np.sum(twopi*Uniq*fax[:,:,:,np.newaxis],axis=2)
1423        dfadui = np.sum(-SQfactor*fa,axis=2)
1424        dfadua = np.sum(-Hij*fa[:,:,:,np.newaxis],axis=2)
1425        #NB: the above have been checked against PA(1:10,1:2) in strfctr.for     
1426        dFdfr[iref] = 2.*(fas[0]*dfadfr[0]+fas[1]*dfadfr[1])*Mdata/len(Uniq)
1427        dFdx[iref] = 2.*(fas[0]*dfadx[0]+fas[1]*dfadx[1])
1428        dFdui[iref] = 2.*(fas[0]*dfadui[0]+fas[1]*dfadui[1])
1429        dFdua[iref] = 2.*(fas[0]*dfadua[0]+fas[1]*dfadua[1])
1430        if not SGData['SGInv']:
1431            dfbdfr = np.sum(fb/occ[:,np.newaxis],axis=2)        #problem here if occ=0 for some atom
1432            dfbdx = np.sum(twopi*Uniq*fbx[:,:,:,np.newaxis],axis=2)         
1433            dfbdui = np.sum(-SQfactor*fb,axis=2)
1434            dfbdua = np.sum(-Hij*fb[:,:,:,np.newaxis],axis=2)
1435            dFdfr[iref] += 2.*(fbs[0]*dfbdfr[0]-fbs[1]*dfbdfr[1])*Mdata/len(Uniq)
1436            dFdx[iref] += 2.*(fbs[0]*dfbdx[0]+fbs[1]*dfbdx[1])
1437            dFdui[iref] += 2.*(fbs[0]*dfbdui[0]-fbs[1]*dfbdui[1])
1438            dFdua[iref] += 2.*(fbs[0]*dfbdua[0]+fbs[1]*dfbdua[1])
1439        #loop over atoms - each dict entry is list of derivatives for all the reflections
1440    for i in range(len(Mdata)):     
1441        dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i]
1442        dFdvDict[pfx+'dAx:'+str(i)] = dFdx.T[0][i]
1443        dFdvDict[pfx+'dAy:'+str(i)] = dFdx.T[1][i]
1444        dFdvDict[pfx+'dAz:'+str(i)] = dFdx.T[2][i]
1445        dFdvDict[pfx+'AUiso:'+str(i)] = dFdui.T[i]
1446        dFdvDict[pfx+'AU11:'+str(i)] = dFdua.T[0][i]
1447        dFdvDict[pfx+'AU22:'+str(i)] = dFdua.T[1][i]
1448        dFdvDict[pfx+'AU33:'+str(i)] = dFdua.T[2][i]
1449        dFdvDict[pfx+'AU12:'+str(i)] = 2.*dFdua.T[3][i]
1450        dFdvDict[pfx+'AU13:'+str(i)] = 2.*dFdua.T[4][i]
1451        dFdvDict[pfx+'AU23:'+str(i)] = 2.*dFdua.T[5][i]
1452    return dFdvDict
1453       
1454def Dict2Values(parmdict, varylist):
1455    '''Use before call to leastsq to setup list of values for the parameters
1456    in parmdict, as selected by key in varylist'''
1457    return [parmdict[key] for key in varylist] 
1458   
1459def Values2Dict(parmdict, varylist, values):
1460    ''' Use after call to leastsq to update the parameter dictionary with
1461    values corresponding to keys in varylist'''
1462    parmdict.update(zip(varylist,values))
1463   
1464def GetNewCellParms(parmDict,varyList):
1465    newCellDict = {}
1466    Ddict = dict(zip(['D11','D22','D33','D12','D13','D23'],['A'+str(i) for i in range(6)]))
1467    for item in varyList:
1468        keys = item.split(':')
1469        if keys[2] in Ddict:
1470            key = keys[0]+'::'+Ddict[keys[2]]
1471            parm = keys[0]+'::'+keys[2]
1472            newCellDict[parm] = [key,parmDict[key]+parmDict[item]]
1473    return newCellDict
1474   
1475def ApplyXYZshifts(parmDict,varyList):
1476    ''' takes atom x,y,z shift and applies it to corresponding atom x,y,z value
1477        input:
1478            parmDict - parameter dictionary
1479            varyList - list of variables
1480        returns:
1481            newAtomDict - dictitemionary of new atomic coordinate names & values;
1482                key is parameter shift name
1483    '''
1484    newAtomDict = {}
1485    for item in parmDict:
1486        if 'dA' in item:
1487            parm = ''.join(item.split('d'))
1488            parmDict[parm] += parmDict[item]
1489            newAtomDict[item] = [parm,parmDict[parm]]
1490    return newAtomDict
1491   
1492def SHTXcal(refl,g,pfx,hfx,SGData,calcControls,parmDict):
1493    IFCoup = 'Bragg' in calcControls[hfx+'instType']
1494    odfCor = 1.0
1495    H = refl[:3]
1496    cell = G2lat.Gmat2cell(g)
1497    Sangls = [parmDict[pfx+'SH omega'],parmDict[pfx+'SH chi'],parmDict[pfx+'SH phi']]
1498    Gangls = [parmDict[hfx+'Omega'],parmDict[hfx+'Chi'],parmDict[hfx+'Phi'],parmDict[hfx+'Azimuth']]
1499    phi,beta = G2lat.CrsAng(H,cell,SGData)
1500    psi,gam,x,x = G2lat.SamAng(refl[5]/2.,Gangls,Sangls,IFCoup) #ignore 2 sets of angle derivs.
1501    SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],parmDict[pfx+'SHmodel'],parmDict[pfx+'SHorder'])
1502    for item in SHnames:
1503        L,M,N = eval(item.strip('C'))
1504        Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
1505        Ksl,x,x = G2lat.GetKsl(L,M,parmDict[pfx+'SHmodel'],psi,gam)
1506        Lnorm = G2lat.Lnorm(L)
1507        odfCor += parmDict[pfx+item]*Lnorm*Kcl*Ksl
1508    return odfCor
1509   
1510def SHTXcalDerv(refl,g,pfx,hfx,SGData,calcControls,parmDict):
1511    FORPI = 12.5663706143592
1512    IFCoup = 'Bragg' in calcControls[hfx+'instType']
1513    odfCor = 1.0
1514    dFdODF = {}
1515    dFdSA = [0,0,0]
1516    H = refl[:3]
1517    cell = G2lat.Gmat2cell(g)
1518    Sangls = [parmDict[pfx+'SH omega'],parmDict[pfx+'SH chi'],parmDict[pfx+'SH phi']]
1519    Gangls = [parmDict[hfx+'Omega'],parmDict[hfx+'Chi'],parmDict[hfx+'Phi'],parmDict[hfx+'Azimuth']]
1520    phi,beta = G2lat.CrsAng(H,cell,SGData)
1521    psi,gam,dPSdA,dGMdA = G2lat.SamAng(refl[5]/2.,Gangls,Sangls,IFCoup)
1522    SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],parmDict[pfx+'SHmodel'],parmDict[pfx+'SHorder'])
1523    for item in SHnames:
1524        L,M,N = eval(item.strip('C'))
1525        Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
1526        Ksl,dKsdp,dKsdg = G2lat.GetKsl(L,M,parmDict[pfx+'SHmodel'],psi,gam)
1527        Lnorm = G2lat.Lnorm(L)
1528        odfCor += parmDict[pfx+item]*Lnorm*Kcl*Ksl
1529        dFdODF[pfx+item] = Lnorm*Kcl*Ksl
1530        for i in range(3):
1531            dFdSA[i] += parmDict[pfx+item]*Lnorm*Kcl*(dKsdp*dPSdA[i]+dKsdg*dGMdA[i])
1532    return odfCor,dFdODF,dFdSA
1533   
1534def SHPOcal(refl,g,phfx,hfx,SGData,calcControls,parmDict):
1535    odfCor = 1.0
1536    H = refl[:3]
1537    cell = G2lat.Gmat2cell(g)
1538    Sangl = [0.,0.,0.]
1539    if 'Bragg' in calcControls[hfx+'instType']:
1540        Gangls = [0.,90.,0.,parmDict[hfx+'Azimuth']]
1541        IFCoup = True
1542    else:
1543        Gangls = [0.,0.,0.,parmDict[hfx+'Azimuth']]
1544        IFCoup = False
1545    phi,beta = G2lat.CrsAng(H,cell,SGData)
1546    psi,gam,x,x = G2lat.SamAng(refl[5]/2.,Gangls,Sangl,IFCoup) #ignore 2 sets of angle derivs.
1547    SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],'0',calcControls[phfx+'SHord'],False)
1548    for item in SHnames:
1549        L,N = eval(item.strip('C'))
1550        Kcsl,Lnorm = G2lat.GetKclKsl(L,N,SGData['SGLaue'],psi,phi,beta) 
1551        odfCor += parmDict[phfx+item]*Lnorm*Kcsl
1552    return odfCor
1553   
1554def SHPOcalDerv(refl,g,phfx,hfx,SGData,calcControls,parmDict):
1555    FORPI = 12.5663706143592
1556    odfCor = 1.0
1557    dFdODF = {}
1558    H = refl[:3]
1559    cell = G2lat.Gmat2cell(g)
1560    Sangl = [0.,0.,0.]
1561    if 'Bragg' in calcControls[hfx+'instType']:
1562        Gangls = [0.,90.,0.,parmDict[hfx+'Azimuth']]
1563        IFCoup = True
1564    else:
1565        Gangls = [0.,0.,0.,parmDict[hfx+'Azimuth']]
1566        IFCoup = False
1567    phi,beta = G2lat.CrsAng(H,cell,SGData)
1568    psi,gam,x,x = G2lat.SamAng(refl[5]/2.,Gangls,Sangl,IFCoup) #ignore 2 sets of angle derivs.
1569    SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],'0',calcControls[phfx+'SHord'],False)
1570    for item in SHnames:
1571        L,N = eval(item.strip('C'))
1572        Kcsl,Lnorm = G2lat.GetKclKsl(L,N,SGData['SGLaue'],psi,phi,beta) 
1573        odfCor += parmDict[phfx+item]*Lnorm*Kcsl
1574        dFdODF[phfx+item] = Kcsl*Lnorm
1575    return odfCor,dFdODF
1576   
1577def GetPrefOri(refl,G,g,phfx,hfx,SGData,calcControls,parmDict):
1578    if calcControls[phfx+'poType'] == 'MD':
1579        MD = parmDict[phfx+'MD']
1580        MDAxis = calcControls[phfx+'MDAxis']
1581        sumMD = 0
1582        for H in refl[11]:           
1583            cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G)
1584            A = 1.0/np.sqrt((MD*cosP)**2+sinP**2/MD)
1585            sumMD += A**3
1586        POcorr = sumMD/len(refl[11])
1587    else:   #spherical harmonics
1588        POcorr = SHPOcal(refl,g,phfx,hfx,SGData,calcControls,parmDict)
1589    return POcorr
1590   
1591def GetPrefOriDerv(refl,G,g,phfx,hfx,SGData,calcControls,parmDict):
1592    POderv = {}
1593    if calcControls[phfx+'poType'] == 'MD':
1594        MD = parmDict[phfx+'MD']
1595        MDAxis = calcControls[phfx+'MDAxis']
1596        sumMD = 0
1597        sumdMD = 0
1598        for H in refl[11]:           
1599            cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G)
1600            A = 1.0/np.sqrt((MD*cosP)**2+sinP**2/MD)
1601            sumMD += A**3
1602            sumdMD -= (1.5*A**5)*(2.0*MD*cosP**2-(sinP/MD)**2)
1603        POcorr = sumMD/len(refl[11])
1604        POderv[phfx+'MD'] = sumdMD/len(refl[11])
1605    else:   #spherical harmonics
1606        POcorr,POderv = SHPOcalDerv(refl,g,phfx,hfx,SGData,calcControls,parmDict)
1607    return POcorr,POderv
1608   
1609def GetIntensityCorr(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict):
1610    Icorr = parmDict[phfx+'Scale']*parmDict[hfx+'Scale']*refl[3]               #scale*multiplicity
1611    if 'X' in parmDict[hfx+'Type']:
1612        Icorr *= G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5],parmDict[hfx+'Azimuth'])[0]
1613    Icorr *= GetPrefOri(refl,G,g,phfx,hfx,SGData,calcControls,parmDict)
1614    if pfx+'SHorder' in parmDict:
1615        Icorr *= SHTXcal(refl,g,pfx,hfx,SGData,calcControls,parmDict)
1616    refl[13] = Icorr       
1617   
1618def GetIntensityDerv(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict):
1619    dIdsh = 1./parmDict[hfx+'Scale']
1620    dIdsp = 1./parmDict[phfx+'Scale']
1621    if 'X' in parmDict[hfx+'Type']:
1622        pola,dIdPola = G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5],parmDict[hfx+'Azimuth'])
1623        dIdPola /= pola
1624    else:       #'N'
1625        dIdPola = 0.0
1626    POcorr,dIdPO = GetPrefOriDerv(refl,G,g,phfx,hfx,SGData,calcControls,parmDict)
1627    for iPO in dIdPO:
1628        dIdPO[iPO] /= POcorr
1629    dFdODF = {}
1630    dFdSA = [0,0,0]
1631    if pfx+'SHorder' in parmDict:
1632        odfCor,dFdODF,dFdSA = SHTXcalDerv(refl,g,pfx,hfx,SGData,calcControls,parmDict)
1633        for iSH in dFdODF:
1634            dFdODF[iSH] /= odfCor
1635        for i in range(3):
1636            dFdSA[i] /= odfCor
1637    return dIdsh,dIdsp,dIdPola,dIdPO,dFdODF,dFdSA
1638       
1639def GetSampleGam(refl,wave,G,GB,phfx,calcControls,parmDict):
1640    costh = cosd(refl[5]/2.)
1641    #crystallite size
1642    if calcControls[phfx+'SizeType'] == 'isotropic':
1643        gam = 1.8*wave/(np.pi*parmDict[phfx+'Size:i']*costh)
1644    elif calcControls[phfx+'SizeType'] == 'uniaxial':
1645        H = np.array(refl[:3])
1646        P = np.array(calcControls[phfx+'SizeAxis'])
1647        cosP,sinP = G2lat.CosSinAngle(H,P,G)
1648        gam = (1.8*wave/np.pi)/(parmDict[phfx+'Size:i']*parmDict[phfx+'Size:a']*costh)
1649        gam *= np.sqrt((sinP*parmDict[phfx+'Size:a'])**2+(cosP*parmDict[phfx+'Size:i'])**2)
1650    else:           #ellipsoidal crystallites
1651        Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
1652        H = np.array(refl[:3])
1653        lenR = G2pwd.ellipseSize(H,Sij,GB)
1654        gam = 1.8*wave/(np.pi*costh*lenR)
1655    #microstrain               
1656    if calcControls[phfx+'MustrainType'] == 'isotropic':
1657        gam += 0.018*parmDict[phfx+'Mustrain:i']*tand(refl[5]/2.)/np.pi
1658    elif calcControls[phfx+'MustrainType'] == 'uniaxial':
1659        H = np.array(refl[:3])
1660        P = np.array(calcControls[phfx+'MustrainAxis'])
1661        cosP,sinP = G2lat.CosSinAngle(H,P,G)
1662        Si = parmDict[phfx+'Mustrain:i']
1663        Sa = parmDict[phfx+'Mustrain:a']
1664        gam += 0.018*Si*Sa*tand(refl[5]/2.)/(np.pi*np.sqrt((Si*cosP)**2+(Sa*sinP)**2))
1665    else:       #generalized - P.W. Stephens model
1666        pwrs = calcControls[phfx+'MuPwrs']
1667        sum = 0
1668        for i,pwr in enumerate(pwrs):
1669            sum += parmDict[phfx+'Mustrain:'+str(i)]*refl[0]**pwr[0]*refl[1]**pwr[1]*refl[2]**pwr[2]
1670        gam += 0.018*refl[4]**2*tand(refl[5]/2.)*sum           
1671    return gam
1672       
1673def GetSampleGamDerv(refl,wave,G,GB,phfx,calcControls,parmDict):
1674    gamDict = {}
1675    costh = cosd(refl[5]/2.)
1676    tanth = tand(refl[5]/2.)
1677    #crystallite size derivatives
1678    if calcControls[phfx+'SizeType'] == 'isotropic':
1679        gamDict[phfx+'Size:i'] = -1.80*wave/(np.pi*costh)
1680    elif calcControls[phfx+'SizeType'] == 'uniaxial':
1681        H = np.array(refl[:3])
1682        P = np.array(calcControls[phfx+'SizeAxis'])
1683        cosP,sinP = G2lat.CosSinAngle(H,P,G)
1684        Si = parmDict[phfx+'Size:i']
1685        Sa = parmDict[phfx+'Size:a']
1686        gami = (1.8*wave/np.pi)/(Si*Sa)
1687        sqtrm = np.sqrt((sinP*Sa)**2+(cosP*Si)**2)
1688        gam = gami*sqtrm/costh           
1689        gamDict[phfx+'Size:i'] = gami*Si*cosP**2/(sqtrm*costh)-gam/Si
1690        gamDict[phfx+'Size:a'] = gami*Sa*sinP**2/(sqtrm*costh)-gam/Sa         
1691    else:           #ellipsoidal crystallites
1692        const = 1.8*wave/(np.pi*costh)
1693        Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
1694        H = np.array(refl[:3])
1695        R,dRdS = G2pwd.ellipseSizeDerv(H,Sij,GB)
1696        for i,item in enumerate([phfx+'Size:%d'%(j) for j in range(6)]):
1697            gamDict[item] = -(const/R**2)*dRdS[i]
1698    #microstrain derivatives               
1699    if calcControls[phfx+'MustrainType'] == 'isotropic':
1700        gamDict[phfx+'Mustrain:i'] =  0.018*tanth/np.pi           
1701    elif calcControls[phfx+'MustrainType'] == 'uniaxial':
1702        H = np.array(refl[:3])
1703        P = np.array(calcControls[phfx+'MustrainAxis'])
1704        cosP,sinP = G2lat.CosSinAngle(H,P,G)
1705        Si = parmDict[phfx+'Mustrain:i']
1706        Sa = parmDict[phfx+'Mustrain:a']
1707        gami = 0.018*Si*Sa*tanth/np.pi
1708        sqtrm = np.sqrt((Si*cosP)**2+(Sa*sinP)**2)
1709        gam = gami/sqtrm
1710        gamDict[phfx+'Mustrain:i'] = gam/Si-gami*Si*cosP**2/sqtrm**3
1711        gamDict[phfx+'Mustrain:a'] = gam/Sa-gami*Sa*sinP**2/sqtrm**3
1712    else:       #generalized - P.W. Stephens model
1713        pwrs = calcControls[phfx+'MuPwrs']
1714        const = 0.018*refl[4]**2*tanth
1715        for i,pwr in enumerate(pwrs):
1716            gamDict[phfx+'Mustrain:'+str(i)] = const*refl[0]**pwr[0]*refl[1]**pwr[1]*refl[2]**pwr[2]
1717    return gamDict
1718       
1719def GetReflPos(refl,wave,G,hfx,calcControls,parmDict):
1720    h,k,l = refl[:3]
1721    dsq = 1./G2lat.calc_rDsq2(np.array([h,k,l]),G)
1722    d = np.sqrt(dsq)
1723
1724    refl[4] = d
1725    pos = 2.0*asind(wave/(2.0*d))+parmDict[hfx+'Zero']
1726    const = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius'])                  #shifts in microns
1727    if 'Bragg' in calcControls[hfx+'instType']:
1728        pos -= const*(4.*parmDict[hfx+'Shift']*cosd(pos/2.0)+ \
1729            parmDict[hfx+'Transparency']*sind(pos)*100.0)            #trans(=1/mueff) in cm
1730    else:               #Debye-Scherrer - simple but maybe not right
1731        pos -= const*(parmDict[hfx+'DisplaceX']*cosd(pos)+parmDict[hfx+'DisplaceY']*sind(pos))
1732    return pos
1733
1734def GetReflPosDerv(refl,wave,A,hfx,calcControls,parmDict):
1735    dpr = 180./np.pi
1736    h,k,l = refl[:3]
1737    dstsq = G2lat.calc_rDsq(np.array([h,k,l]),A)
1738    dst = np.sqrt(dstsq)
1739    pos = refl[5]-parmDict[hfx+'Zero']
1740    const = dpr/np.sqrt(1.0-wave**2*dstsq/4.0)
1741    dpdw = const*dst
1742    dpdA = np.array([h**2,k**2,l**2,h*k,h*l,k*l])
1743    dpdA *= const*wave/(2.0*dst)
1744    dpdZ = 1.0
1745    const = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius'])                  #shifts in microns
1746    if 'Bragg' in calcControls[hfx+'instType']:
1747        dpdSh = -4.*const*cosd(pos/2.0)
1748        dpdTr = -const*sind(pos)*100.0
1749        return dpdA,dpdw,dpdZ,dpdSh,dpdTr,0.,0.
1750    else:               #Debye-Scherrer - simple but maybe not right
1751        dpdXd = -const*cosd(pos)
1752        dpdYd = -const*sind(pos)
1753        return dpdA,dpdw,dpdZ,0.,0.,dpdXd,dpdYd
1754           
1755def GetHStrainShift(refl,SGData,phfx,parmDict):
1756    laue = SGData['SGLaue']
1757    uniq = SGData['SGUniq']
1758    h,k,l = refl[:3]
1759    if laue in ['m3','m3m']:
1760        Dij = parmDict[phfx+'D11']*(h**2+k**2+l**2)+ \
1761            refl[4]**2*parmDict[phfx+'eA']*((h*k)**2+(h*l)**2+(k*l)**2)/(h**2+k**2+l**2)**2
1762    elif laue in ['6/m','6/mmm','3m1','31m','3']:
1763        Dij = parmDict[phfx+'D11']*(h**2+k**2+h*k)+parmDict[phfx+'D33']*l**2
1764    elif laue in ['3R','3mR']:
1765        Dij = parmDict[phfx+'D11']*(h**2+k**2+l**2)+parmDict[phfx+'D12']*(h*k+h*l+k*l)
1766    elif laue in ['4/m','4/mmm']:
1767        Dij = parmDict[phfx+'D11']*(h**2+k**2)+parmDict[phfx+'D33']*l**2
1768    elif laue in ['mmm']:
1769        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2
1770    elif laue in ['2/m']:
1771        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2
1772        if uniq == 'a':
1773            Dij += parmDict[phfx+'D23']*k*l
1774        elif uniq == 'b':
1775            Dij += parmDict[phfx+'D13']*h*l
1776        elif uniq == 'c':
1777            Dij += parmDict[phfx+'D12']*h*k
1778    else:
1779        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2+ \
1780            parmDict[phfx+'D12']*h*k+parmDict[phfx+'D13']*h*l+parmDict[phfx+'D23']*k*l
1781    return Dij*refl[4]**2*tand(refl[5]/2.0)
1782           
1783def GetHStrainShiftDerv(refl,SGData,phfx):
1784    laue = SGData['SGLaue']
1785    uniq = SGData['SGUniq']
1786    h,k,l = refl[:3]
1787    if laue in ['m3','m3m']:
1788        dDijDict = {phfx+'D11':h**2+k**2+l**2,
1789            phfx+'eA':((h*k)**2+(h*l)**2+(k*l)**2)/(h**2+k**2+l**2)**2}
1790    elif laue in ['6/m','6/mmm','3m1','31m','3']:
1791        dDijDict = {phfx+'D11':h**2+k**2+h*k,phfx+'D33':l**2}
1792    elif laue in ['3R','3mR']:
1793        dDijDict = {phfx+'D11':h**2+k**2+l**2,phfx+'D12':h*k+h*l+k*l}
1794    elif laue in ['4/m','4/mmm']:
1795        dDijDict = {phfx+'D11':h**2+k**2,phfx+'D33':l**2}
1796    elif laue in ['mmm']:
1797        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2}
1798    elif laue in ['2/m']:
1799        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2}
1800        if uniq == 'a':
1801            dDijDict[phfx+'D23'] = k*l
1802        elif uniq == 'b':
1803            dDijDict[phfx+'D13'] = h*l
1804        elif uniq == 'c':
1805            dDijDict[phfx+'D12'] = h*k
1806            names.append()
1807    else:
1808        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2,
1809            phfx+'D12':h*k,phfx+'D13':h*l,phfx+'D23':k*l}
1810    for item in dDijDict:
1811        dDijDict[item] *= refl[4]**2*tand(refl[5]/2.0)
1812    return dDijDict
1813   
1814def GetFprime(controlDict,Histograms):
1815    FFtables = controlDict['FFtables']
1816    if not FFtables:
1817        return
1818    histoList = Histograms.keys()
1819    histoList.sort()
1820    for histogram in histoList:
1821        if 'PWDR' in histogram[:4]:
1822            Histogram = Histograms[histogram]
1823            hId = Histogram['hId']
1824            hfx = ':%d:'%(hId)
1825            keV = controlDict[hfx+'keV']
1826            for El in FFtables:
1827                Orbs = G2el.GetXsectionCoeff(El.split('+')[0].split('-')[0])
1828                FP,FPP,Mu = G2el.FPcalc(Orbs, keV)
1829                FFtables[El][hfx+'FP'] = FP
1830                FFtables[El][hfx+'FPP'] = FPP               
1831           
1832def getPowderProfile(parmDict,x,varylist,Histogram,Phases,calcControls,pawleyLookup):
1833   
1834    def GetReflSIgGam(refl,wave,G,GB,hfx,phfx,calcControls,parmDict):
1835        U = parmDict[hfx+'U']
1836        V = parmDict[hfx+'V']
1837        W = parmDict[hfx+'W']
1838        X = parmDict[hfx+'X']
1839        Y = parmDict[hfx+'Y']
1840        tanPos = tand(refl[5]/2.0)
1841        sig = U*tanPos**2+V*tanPos+W        #save peak sigma
1842        sig = max(0.001,sig)
1843        gam = X/cosd(refl[5]/2.0)+Y*tanPos+GetSampleGam(refl,wave,G,GB,phfx,calcControls,parmDict) #save peak gamma
1844        gam = max(0.001,gam)
1845        return sig,gam
1846               
1847    hId = Histogram['hId']
1848    hfx = ':%d:'%(hId)
1849    bakType = calcControls[hfx+'bakType']
1850    yb = G2pwd.getBackground(hfx,parmDict,bakType,x)
1851    yc = np.zeros_like(yb)
1852       
1853    if 'C' in calcControls[hfx+'histType']:   
1854        shl = max(parmDict[hfx+'SH/L'],0.002)
1855        Ka2 = False
1856        if hfx+'Lam1' in parmDict.keys():
1857            wave = parmDict[hfx+'Lam1']
1858            Ka2 = True
1859            lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
1860            kRatio = parmDict[hfx+'I(L2)/I(L1)']
1861        else:
1862            wave = parmDict[hfx+'Lam']
1863    else:
1864        print 'TOF Undefined at present'
1865        raise ValueError
1866    for phase in Histogram['Reflection Lists']:
1867        refList = Histogram['Reflection Lists'][phase]
1868        Phase = Phases[phase]
1869        pId = Phase['pId']
1870        pfx = '%d::'%(pId)
1871        phfx = '%d:%d:'%(pId,hId)
1872        hfx = ':%d:'%(hId)
1873        SGData = Phase['General']['SGData']
1874        A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
1875        G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
1876        GA,GB = G2lat.Gmat2AB(G)    #Orthogonalization matricies
1877        Vst = np.sqrt(nl.det(G))    #V*
1878        if 'Pawley' not in Phase['General']['Type']:
1879            refList = StructureFactor(refList,G,hfx,pfx,SGData,calcControls,parmDict)
1880        for refl in refList:
1881            if 'C' in calcControls[hfx+'histType']:
1882                h,k,l = refl[:3]
1883                refl[5] = GetReflPos(refl,wave,G,hfx,calcControls,parmDict)         #corrected reflection position
1884                Lorenz = 1./(2.*sind(refl[5]/2.)**2*cosd(refl[5]/2.))           #Lorentz correction
1885                refl[5] += GetHStrainShift(refl,SGData,phfx,parmDict)               #apply hydrostatic strain shift
1886                refl[6:8] = GetReflSIgGam(refl,wave,G,GB,hfx,phfx,calcControls,parmDict)    #peak sig & gam
1887                GetIntensityCorr(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)    #puts corrections in refl[13]
1888                refl[13] *= Vst*Lorenz
1889                if 'Pawley' in Phase['General']['Type']:
1890                    try:
1891                        refl[9] = abs(parmDict[pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])])
1892                    except KeyError:
1893#                        print ' ***Error %d,%d,%d missing from Pawley reflection list ***'%(h,k,l)
1894                        continue
1895                Wd,fmin,fmax = G2pwd.getWidths(refl[5],refl[6],refl[7],shl)
1896                iBeg = np.searchsorted(x,refl[5]-fmin)
1897                iFin = np.searchsorted(x,refl[5]+fmax)
1898                if not iBeg+iFin:       #peak below low limit - skip peak
1899                    continue
1900                elif not iBeg-iFin:     #peak above high limit - done
1901                    break
1902                yc[iBeg:iFin] += refl[13]*refl[9]*G2pwd.getFCJVoigt3(refl[5],refl[6],refl[7],shl,x[iBeg:iFin])    #>90% of time spent here
1903                if Ka2:
1904                    pos2 = refl[5]+lamRatio*tand(refl[5]/2.0)       # + 360/pi * Dlam/lam * tan(th)
1905                    Wd,fmin,fmax = G2pwd.getWidths(pos2,refl[6],refl[7],shl)
1906                    iBeg = np.searchsorted(x,pos2-fmin)
1907                    iFin = np.searchsorted(x,pos2+fmax)
1908                    if not iBeg+iFin:       #peak below low limit - skip peak
1909                        continue
1910                    elif not iBeg-iFin:     #peak above high limit - done
1911                        return yc,yb
1912                    yc[iBeg:iFin] += refl[13]*refl[9]*kRatio*G2pwd.getFCJVoigt3(pos2,refl[6],refl[7],shl,x[iBeg:iFin])        #and here
1913            elif 'T' in calcControls[hfx+'histType']:
1914                print 'TOF Undefined at present'
1915                raise Exception    #no TOF yet
1916    return yc,yb
1917   
1918def GetFobsSq(Histograms,Phases,parmDict,calcControls):
1919    histoList = Histograms.keys()
1920    histoList.sort()
1921    for histogram in histoList:
1922        if 'PWDR' in histogram[:4]:
1923            Histogram = Histograms[histogram]
1924            hId = Histogram['hId']
1925            hfx = ':%d:'%(hId)
1926            Limits = calcControls[hfx+'Limits']
1927            shl = max(parmDict[hfx+'SH/L'],0.002)
1928            Ka2 = False
1929            kRatio = 0.0
1930            if hfx+'Lam1' in parmDict.keys():
1931                Ka2 = True
1932                lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
1933                kRatio = parmDict[hfx+'I(L2)/I(L1)']
1934            x,y,w,yc,yb,yd = Histogram['Data']
1935            ymb = np.array(y-yb)
1936            ycmb = np.array(yc-yb)
1937            ratio = np.where(ycmb!=0.,ymb/ycmb,0.0)           
1938            xB = np.searchsorted(x,Limits[0])
1939            xF = np.searchsorted(x,Limits[1])
1940            refLists = Histogram['Reflection Lists']
1941            for phase in refLists:
1942                Phase = Phases[phase]
1943                pId = Phase['pId']
1944                phfx = '%d:%d:'%(pId,hId)
1945                refList = refLists[phase]
1946                sumFo = 0.0
1947                sumdF = 0.0
1948                sumFosq = 0.0
1949                sumdFsq = 0.0
1950                for refl in refList:
1951                    if 'C' in calcControls[hfx+'histType']:
1952                        yp = np.zeros_like(yb)
1953                        Wd,fmin,fmax = G2pwd.getWidths(refl[5],refl[6],refl[7],shl)
1954                        iBeg = np.searchsorted(x[xB:xF],refl[5]-fmin)
1955                        iFin = np.searchsorted(x[xB:xF],refl[5]+fmax)
1956                        iFin2 = iFin
1957                        yp[iBeg:iFin] = refl[13]*refl[9]*G2pwd.getFCJVoigt3(refl[5],refl[6],refl[7],shl,x[iBeg:iFin])    #>90% of time spent here                           
1958                        if Ka2:
1959                            pos2 = refl[5]+lamRatio*tand(refl[5]/2.0)       # + 360/pi * Dlam/lam * tan(th)
1960                            Wd,fmin,fmax = G2pwd.getWidths(pos2,refl[6],refl[7],shl)
1961                            iBeg2 = np.searchsorted(x,pos2-fmin)
1962                            iFin2 = np.searchsorted(x,pos2+fmax)
1963                            yp[iBeg2:iFin2] += refl[13]*refl[9]*kRatio*G2pwd.getFCJVoigt3(pos2,refl[6],refl[7],shl,x[iBeg2:iFin2])        #and here
1964                        refl[8] = np.sum(np.where(ratio[iBeg:iFin2]>0.,yp[iBeg:iFin2]*ratio[iBeg:iFin2]/(refl[13]*(1.+kRatio)),0.0))
1965                    elif 'T' in calcControls[hfx+'histType']:
1966                        print 'TOF Undefined at present'
1967                        raise Exception    #no TOF yet
1968                    Fo = np.sqrt(np.abs(refl[8]))
1969                    Fc = np.sqrt(np.abs(refl[9]))
1970                    sumFo += Fo
1971                    sumFosq += refl[8]**2
1972                    sumdF += np.abs(Fo-Fc)
1973                    sumdFsq += (refl[8]-refl[9])**2
1974                Histogram[phfx+'Rf'] = min(100.,(sumdF/sumFo)*100.)
1975                Histogram[phfx+'Rf^2'] = min(100.,np.sqrt(sumdFsq/sumFosq)*100.)
1976                Histogram[phfx+'Nref'] = len(refList)
1977               
1978def getPowderProfileDerv(parmDict,x,varylist,Histogram,Phases,calcControls,pawleyLookup):
1979   
1980    def cellVaryDerv(pfx,SGData,dpdA): 
1981        if SGData['SGLaue'] in ['-1',]:
1982            return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],
1983                [pfx+'A3',dpdA[3]],[pfx+'A4',dpdA[4]],[pfx+'A5',dpdA[5]]]
1984        elif SGData['SGLaue'] in ['2/m',]:
1985            if SGData['SGUniq'] == 'a':
1986                return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A3',dpdA[3]]]
1987            elif SGData['SGUniq'] == 'b':
1988                return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A4',dpdA[4]]]
1989            else:
1990                return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A5',dpdA[5]]]
1991        elif SGData['SGLaue'] in ['mmm',]:
1992            return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]]]
1993        elif SGData['SGLaue'] in ['4/m','4/mmm']:
1994#            return [[pfx+'A0',dpdA[0]+dpdA[1]],[pfx+'A2',dpdA[2]]]
1995            return [[pfx+'A0',dpdA[0]],[pfx+'A2',dpdA[2]]]
1996        elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
1997#            return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[3]],[pfx+'A2',dpdA[2]]]
1998            return [[pfx+'A0',dpdA[0]],[pfx+'A2',dpdA[2]]]
1999        elif SGData['SGLaue'] in ['3R', '3mR']:
2000            return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[2]],[pfx+'A3',dpdA[3]+dpdA[4]+dpdA[5]]]                       
2001        elif SGData['SGLaue'] in ['m3m','m3']:
2002#            return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[2]]]
2003            return [[pfx+'A0',dpdA[0]]]
2004    # create a list of dependent variables and set up a dictionary to hold their derivatives
2005    dependentVars = G2mv.GetDependentVars()
2006    depDerivDict = {}
2007    for j in dependentVars:
2008        depDerivDict[j] = np.zeros(shape=(len(x)))
2009    #print 'dependent vars',dependentVars
2010    lenX = len(x)               
2011    hId = Histogram['hId']
2012    hfx = ':%d:'%(hId)
2013    bakType = calcControls[hfx+'bakType']
2014    dMdv = np.zeros(shape=(len(varylist),len(x)))
2015    dMdb,dMddb = G2pwd.getBackgroundDerv(hfx,parmDict,bakType,x)
2016    if hfx+'Back:0' in varylist: # for now assume that Back:x vars to not appear in constraints
2017        bBpos =varylist.index(hfx+'Back:0')
2018        dMdv[bBpos:bBpos+len(dMdb)] = dMdb
2019    names = [hfx+'DebyeA',hfx+'DebyeR',hfx+'DebyeU']
2020    for name in varylist:
2021        if 'Debye' in name:
2022            id = int(name.split(':')[-1])
2023            parm = name[:int(name.rindex(':'))]
2024            ip = names.index(parm)
2025            dMdv[varylist.index(name)] = dMddb[3*id+ip]
2026    if 'C' in calcControls[hfx+'histType']:   
2027        dx = x[1]-x[0]
2028        shl = max(parmDict[hfx+'SH/L'],0.002)
2029        Ka2 = False
2030        if hfx+'Lam1' in parmDict.keys():
2031            wave = parmDict[hfx+'Lam1']
2032            Ka2 = True
2033            lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
2034            kRatio = parmDict[hfx+'I(L2)/I(L1)']
2035        else:
2036            wave = parmDict[hfx+'Lam']
2037    else:
2038        print 'TOF Undefined at present'
2039        raise ValueError
2040    for phase in Histogram['Reflection Lists']:
2041        refList = Histogram['Reflection Lists'][phase]
2042        Phase = Phases[phase]
2043        SGData = Phase['General']['SGData']
2044        pId = Phase['pId']
2045        pfx = '%d::'%(pId)
2046        phfx = '%d:%d:'%(pId,hId)
2047        A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
2048        G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
2049        GA,GB = G2lat.Gmat2AB(G)    #Orthogonalization matricies
2050        if 'Pawley' not in Phase['General']['Type']:
2051            dFdvDict = StructureFactorDerv(refList,G,hfx,pfx,SGData,calcControls,parmDict)
2052        for iref,refl in enumerate(refList):
2053            if 'C' in calcControls[hfx+'histType']:        #CW powder
2054                h,k,l = refl[:3]
2055                dIdsh,dIdsp,dIdpola,dIdPO,dFdODF,dFdSA = GetIntensityDerv(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
2056                if 'Pawley' in Phase['General']['Type']:
2057                    try:
2058                        refl[9] = abs(parmDict[pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])])
2059                    except KeyError:
2060#                        print ' ***Error %d,%d,%d missing from Pawley reflection list ***'%(h,k,l)
2061                        continue
2062                Wd,fmin,fmax = G2pwd.getWidths(refl[5],refl[6],refl[7],shl)
2063                iBeg = np.searchsorted(x,refl[5]-fmin)
2064                iFin = np.searchsorted(x,refl[5]+fmax)
2065                if not iBeg+iFin:       #peak below low limit - skip peak
2066                    continue
2067                elif not iBeg-iFin:     #peak above high limit - done
2068                    break
2069                pos = refl[5]
2070                tanth = tand(pos/2.0)
2071                costh = cosd(pos/2.0)
2072                dMdpk = np.zeros(shape=(6,len(x)))
2073                dMdipk = G2pwd.getdFCJVoigt3(refl[5],refl[6],refl[7],shl,x[iBeg:iFin])
2074                for i in range(1,5):
2075                    dMdpk[i][iBeg:iFin] += 100.*dx*refl[13]*refl[9]*dMdipk[i]
2076                dMdpk[0][iBeg:iFin] += 100.*dx*refl[13]*refl[9]*dMdipk[0]
2077                dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4]}
2078                if Ka2:
2079                    pos2 = refl[5]+lamRatio*tanth       # + 360/pi * Dlam/lam * tan(th)
2080                    kdelt = int((pos2-refl[5])/dx)               
2081                    iBeg2 = min(lenX,iBeg+kdelt)
2082                    iFin2 = min(lenX,iFin+kdelt)
2083                    if iBeg2-iFin2:
2084                        dMdipk2 = G2pwd.getdFCJVoigt3(pos2,refl[6],refl[7],shl,x[iBeg2:iFin2])
2085                        for i in range(1,5):
2086                            dMdpk[i][iBeg2:iFin2] += 100.*dx*refl[13]*refl[9]*kRatio*dMdipk2[i]
2087                        dMdpk[0][iBeg2:iFin2] += 100.*dx*refl[13]*refl[9]*kRatio*dMdipk2[0]
2088                        dMdpk[5][iBeg2:iFin2] += 100.*dx*refl[13]*dMdipk2[0]
2089                        dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':dMdpk[5]*refl[9]}
2090                if 'Pawley' in Phase['General']['Type']:
2091                    try:
2092                        idx = varylist.index(pfx+'PWLref:'+str(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)]))
2093                        dMdv[idx] = dervDict['int']/refl[9]
2094                        # Assuming Pawley variables not in constraints
2095                    except ValueError:
2096                        pass
2097                dpdA,dpdw,dpdZ,dpdSh,dpdTr,dpdX,dpdY = GetReflPosDerv(refl,wave,A,hfx,calcControls,parmDict)
2098                names = {hfx+'Scale':[dIdsh,'int'],hfx+'Polariz.':[dIdpola,'int'],phfx+'Scale':[dIdsp,'int'],
2099                    hfx+'U':[tanth**2,'sig'],hfx+'V':[tanth,'sig'],hfx+'W':[1.0,'sig'],
2100                    hfx+'X':[1.0/costh,'gam'],hfx+'Y':[tanth,'gam'],hfx+'SH/L':[1.0,'shl'],
2101                    hfx+'I(L2)/I(L1)':[1.0,'L1/L2'],hfx+'Zero':[dpdZ,'pos'],hfx+'Lam':[dpdw,'pos'],
2102                    hfx+'Shift':[dpdSh,'pos'],hfx+'Transparency':[dpdTr,'pos'],hfx+'DisplaceX':[dpdX,'pos'],
2103                    hfx+'DisplaceY':[dpdY,'pos'],}
2104                for name in names:
2105                    item = names[name]
2106                    if name in varylist:
2107                        dMdv[varylist.index(name)] += item[0]*dervDict[item[1]]
2108                    elif name in dependentVars:
2109                        depDerivDict[name] += item[0]*dervDict[item[1]]
2110
2111                for iPO in dIdPO:
2112                    if iPO in varylist:
2113                        dMdv[varylist.index(iPO)] += dIdPO[iPO]*dervDict['int']
2114                    elif iPO in dependentVars:
2115                        depDerivDict[iPO] = dIdPO[iPO]*dervDict['int']
2116
2117                for i,name in enumerate(['omega','chi','phi']):
2118                    aname = pfx+'SH '+name
2119                    if aname in varylist:
2120                        dMdv[varylist.index(aname)] += dFdSA[i]*dervDict['int']
2121                    elif aname in dependentVars:
2122                        depDerivDict[aname] += dFdSA[i]*dervDict['int']
2123                for iSH in dFdODF:
2124                    if iSH in varylist:
2125                        dMdv[varylist.index(iSH)] += dFdODF[iSH]*dervDict['int']
2126                    elif iSH in dependentVars:
2127                        depDerivDict[iSH] += dFdODF[iSH]*dervDict['int']
2128                cellDervNames = cellVaryDerv(pfx,SGData,dpdA)
2129                for name,dpdA in cellDervNames:
2130                    if name in varylist:
2131                        dMdv[varylist.index(name)] += dpdA*dervDict['pos']
2132                    elif name in dependentVars:
2133                        depDerivDict[name] += dpdA*dervDict['pos']
2134                dDijDict = GetHStrainShiftDerv(refl,SGData,phfx)
2135                for name in dDijDict:
2136                    if name in varylist:
2137                        dMdv[varylist.index(name)] += dDijDict[name]*dervDict['pos']
2138                    elif name in dependentVars:
2139                        depDerivDict[name] += dDijDict[name]*dervDict['pos']
2140                gamDict = GetSampleGamDerv(refl,wave,G,GB,phfx,calcControls,parmDict)
2141                for name in gamDict:
2142                    if name in varylist:
2143                        dMdv[varylist.index(name)] += gamDict[name]*dervDict['gam']
2144                    elif name in dependentVars:
2145                        depDerivDict[name] += gamDict[name]*dervDict['gam']
2146                                               
2147            elif 'T' in calcControls[hfx+'histType']:
2148                print 'TOF Undefined at present'
2149                raise Exception    #no TOF yet
2150            #do atom derivatives -  for F,X & U so far             
2151            corr = dervDict['int']/refl[9]
2152            for name in varylist+dependentVars:
2153                try:
2154                    aname = name.split(pfx)[1][:2]
2155                    if aname not in ['Af','dA','AU']: continue # skip anything not an atom param
2156                except IndexError:
2157                    continue
2158                if name in varylist:
2159                    dMdv[varylist.index(name)] += dFdvDict[name][iref]*corr
2160                elif name in dependentVars:
2161                    depDerivDict[name] += dFdvDict[name][iref]*corr
2162    # now process derivatives in constraints
2163    G2mv.Dict2Deriv(varylist,depDerivDict,dMdv)
2164    return dMdv
2165
2166def dervRefine(values,HistoPhases,parmdict,varylist,calcControls,pawleyLookup,dlg):
2167    parmdict.update(zip(varylist,values))
2168    G2mv.Dict2Map(parmdict,varylist)
2169    Histograms,Phases = HistoPhases
2170    nvar = len(varylist)
2171    dMdv = np.empty(0)
2172    histoList = Histograms.keys()
2173    histoList.sort()
2174    for histogram in histoList:
2175        if 'PWDR' in histogram[:4]:
2176            Histogram = Histograms[histogram]
2177            hId = Histogram['hId']
2178            hfx = ':%d:'%(hId)
2179            Limits = calcControls[hfx+'Limits']
2180            x,y,w,yc,yb,yd = Histogram['Data']
2181            xB = np.searchsorted(x,Limits[0])
2182            xF = np.searchsorted(x,Limits[1])
2183            dMdvh = np.sqrt(w[xB:xF])*getPowderProfileDerv(parmdict,x[xB:xF],
2184                varylist,Histogram,Phases,calcControls,pawleyLookup)
2185            if len(dMdv):
2186                dMdv = np.concatenate((dMdv.T,dMdvh.T)).T
2187            else:
2188                dMdv = dMdvh
2189    return dMdv
2190
2191def HessRefine(values,HistoPhases,parmdict,varylist,calcControls,pawleyLookup,dlg):
2192    parmdict.update(zip(varylist,values))
2193    G2mv.Dict2Map(parmdict,varylist)
2194    Histograms,Phases = HistoPhases
2195    nvar = len(varylist)
2196    Hess = np.empty(0)
2197    histoList = Histograms.keys()
2198    histoList.sort()
2199    for histogram in histoList:
2200        if 'PWDR' in histogram[:4]:
2201            Histogram = Histograms[histogram]
2202            hId = Histogram['hId']
2203            hfx = ':%d:'%(hId)
2204            Limits = calcControls[hfx+'Limits']
2205            x,y,w,yc,yb,yd = Histogram['Data']
2206            dy = y-yc
2207            xB = np.searchsorted(x,Limits[0])
2208            xF = np.searchsorted(x,Limits[1])
2209            dMdvh = np.sqrt(w[xB:xF])*getPowderProfileDerv(parmdict,x[xB:xF],
2210                varylist,Histogram,Phases,calcControls,pawleyLookup)
2211            if dlg:
2212                dlg.Update(Histogram['wRp'],newmsg='Hessian for histogram %d Rwp=%8.3f%s'%(hId,Histogram['wRp'],'%'))[0]
2213            if len(Hess):
2214                Vec += np.sum(dMdvh*np.sqrt(w[xB:xF])*dy[xB:xF],axis=1)
2215                Hess += np.inner(dMdvh,dMdvh)
2216            else:
2217                Vec = np.sum(dMdvh*np.sqrt(w[xB:xF])*dy[xB:xF],axis=1)
2218                Hess = np.inner(dMdvh,dMdvh)
2219    return Vec,Hess
2220
2221def errRefine(values,HistoPhases,parmdict,varylist,calcControls,pawleyLookup,dlg):       
2222    parmdict.update(zip(varylist,values))
2223    Values2Dict(parmdict, varylist, values)
2224    G2mv.Dict2Map(parmdict,varylist)
2225    Histograms,Phases = HistoPhases
2226    M = np.empty(0)
2227    sumwYo = 0
2228    Nobs = 0
2229    histoList = Histograms.keys()
2230    histoList.sort()
2231    for histogram in histoList:
2232        if 'PWDR' in histogram[:4]:
2233            Histogram = Histograms[histogram]
2234            hId = Histogram['hId']
2235            hfx = ':%d:'%(hId)
2236            Limits = calcControls[hfx+'Limits']
2237            x,y,w,yc,yb,yd = Histogram['Data']
2238            yc *= 0.0                           #zero full calcd profiles
2239            yb *= 0.0
2240            yd *= 0.0
2241            xB = np.searchsorted(x,Limits[0])
2242            xF = np.searchsorted(x,Limits[1])
2243            Histogram['Nobs'] = xF-xB
2244            Nobs += Histogram['Nobs']
2245            Histogram['sumwYo'] = np.sum(w[xB:xF]*y[xB:xF]**2)
2246            sumwYo += Histogram['sumwYo']
2247            yc[xB:xF],yb[xB:xF] = getPowderProfile(parmdict,x[xB:xF],
2248                varylist,Histogram,Phases,calcControls,pawleyLookup)
2249            yc[xB:xF] += yb[xB:xF]
2250            yd[xB:xF] = y[xB:xF]-yc[xB:xF]
2251            Histogram['sumwYd'] = np.sum(np.sqrt(w[xB:xF])*(yd[xB:xF]))
2252            wdy = -np.sqrt(w[xB:xF])*(yd[xB:xF])
2253            Histogram['wRp'] = min(100.,np.sqrt(np.sum(wdy**2)/Histogram['sumwYo'])*100.)
2254            if dlg:
2255                dlg.Update(Histogram['wRp'],newmsg='For histogram %d Rwp=%8.3f%s'%(hId,Histogram['wRp'],'%'))[0]
2256            M = np.concatenate((M,wdy))
2257    Histograms['sumwYo'] = sumwYo
2258    Histograms['Nobs'] = Nobs
2259    Rwp = min(100.,np.sqrt(np.sum(M**2)/sumwYo)*100.)
2260    if dlg:
2261        GoOn = dlg.Update(Rwp,newmsg='%s%8.3f%s'%('Powder profile Rwp =',Rwp,'%'))[0]
2262        if not GoOn:
2263            parmDict['saved values'] = values
2264            raise Exception         #Abort!!
2265    return M
2266   
2267                   
2268def Refine(GPXfile,dlg):
2269    import pytexture as ptx
2270    ptx.pyqlmninit()            #initialize fortran arrays for spherical harmonics
2271   
2272    ShowBanner()
2273    varyList = []
2274    parmDict = {}
2275    calcControls = {}
2276    G2mv.InitVars()   
2277    Controls = G2IO.GetControls(GPXfile)
2278    ShowControls(Controls)           
2279    constrDict,constrFlag,fixedList = G2IO.GetConstraints(GPXfile)
2280    Histograms,Phases = G2IO.GetUsedHistogramsAndPhases(GPXfile)
2281    if not Phases:
2282        print ' *** ERROR - you have no histograms to refine! ***'
2283        print ' *** Refine aborted ***'
2284        raise Exception
2285    if not Histograms:
2286        print ' *** ERROR - you have no data to refine with! ***'
2287        print ' *** Refine aborted ***'
2288        raise Exception       
2289    Natoms,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables = GetPhaseData(Phases)
2290    calcControls['Natoms'] = Natoms
2291    calcControls['FFtables'] = FFtables
2292    calcControls['BLtables'] = BLtables
2293    hapVary,hapDict,controlDict = GetHistogramPhaseData(Phases,Histograms)
2294    calcControls.update(controlDict)
2295    histVary,histDict,controlDict = GetHistogramData(Histograms)
2296    calcControls.update(controlDict)
2297    varyList = phaseVary+hapVary+histVary
2298    parmDict.update(phaseDict)
2299    parmDict.update(hapDict)
2300    parmDict.update(histDict)
2301    GetFprime(calcControls,Histograms)
2302    # do constraint processing
2303    try:
2304        groups,parmlist = G2mv.GroupConstraints(constrDict)
2305        G2mv.GenerateConstraints(groups,parmlist,varyList,constrDict,constrFlag,fixedList)
2306    except:
2307        print ' *** ERROR - your constraints are internally inconsistent ***'
2308        raise Exception(' *** Refine aborted ***')
2309    # check to see which generated parameters are fully varied
2310    msg = G2mv.SetVaryFlags(varyList)
2311    if msg:
2312        print ' *** ERROR - you have not set the refine flags for constraints consistently! ***'
2313        print msg
2314        raise Exception(' *** Refine aborted ***')
2315    G2mv.Map2Dict(parmDict,varyList)
2316#    print G2mv.VarRemapShow(varyList)
2317
2318    while True:
2319        begin = time.time()
2320        values =  np.array(Dict2Values(parmDict, varyList))
2321        Ftol = Controls['min dM/M']
2322        Factor = Controls['shift factor']
2323        maxCyc = Controls['max cyc']
2324        if 'Jacobian' in Controls['deriv type']:           
2325            result = so.leastsq(errRefine,values,Dfun=dervRefine,full_output=True,
2326                ftol=Ftol,col_deriv=True,factor=Factor,
2327                args=([Histograms,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg))
2328            ncyc = int(result[2]['nfev']/2)
2329        elif 'Hessian' in Controls['deriv type']:
2330            result = G2mth.HessianLSQ(errRefine,values,Hess=HessRefine,ftol=Ftol,maxcyc=maxCyc,
2331                args=([Histograms,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg))
2332            ncyc = result[2]['num cyc']+1                           
2333        else:           #'numeric'
2334            result = so.leastsq(errRefine,values,full_output=True,ftol=Ftol,epsfcn=1.e-8,factor=Factor,
2335                args=([Histograms,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg))
2336            ncyc = int(result[2]['nfev']/len(varyList))
2337#        table = dict(zip(varyList,zip(values,result[0],(result[0]-values))))
2338#        for item in table: print item,table[item]               #useful debug - are things shifting?
2339        runtime = time.time()-begin
2340        chisq = np.sum(result[2]['fvec']**2)
2341        Values2Dict(parmDict, varyList, result[0])
2342        G2mv.Dict2Map(parmDict,varyList)
2343       
2344        Rwp = np.sqrt(chisq/Histograms['sumwYo'])*100.      #to %
2345        GOF = chisq/(Histograms['Nobs']-len(varyList))
2346        print '\n Refinement results:'
2347        print 135*'-'
2348        print ' Number of function calls:',result[2]['nfev'],' Number of observations: ',Histograms['Nobs'],' Number of parameters: ',len(varyList)
2349        print ' Refinement time = %8.3fs, %8.3fs/cycle, for %d cycles'%(runtime,runtime/ncyc,ncyc)
2350        print ' wRp = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f'%(Rwp,chisq,GOF)
2351        try:
2352            covMatrix = result[1]*GOF
2353            sig = np.sqrt(np.diag(covMatrix))
2354            if np.any(np.isnan(sig)):
2355                print '*** Least squares aborted - some invalid esds possible ***'
2356#            table = dict(zip(varyList,zip(values,result[0],(result[0]-values)/sig)))
2357#            for item in table: print item,table[item]               #useful debug - are things shifting?
2358            break                   #refinement succeeded - finish up!
2359        except TypeError:          #result[1] is None on singular matrix
2360            print '**** Refinement failed - singular matrix ****'
2361            if 'Hessian' in Controls['deriv type']:
2362                for i in result[2]['psing'].reverse():
2363                        print 'Removing parameter: ',varyList[i]
2364                        del(varyList[i])                   
2365            else:
2366                Ipvt = result[2]['ipvt']
2367                for i,ipvt in enumerate(Ipvt):
2368                    if not np.sum(result[2]['fjac'],axis=1)[i]:
2369                        print 'Removing parameter: ',varyList[ipvt-1]
2370                        del(varyList[ipvt-1])
2371                        break
2372
2373#    print 'dependentParmList: ',G2mv.dependentParmList
2374#    print 'arrayList: ',G2mv.arrayList
2375#    print 'invarrayList: ',G2mv.invarrayList
2376#    print 'indParmList: ',G2mv.indParmList
2377#    print 'fixedDict: ',G2mv.fixedDict
2378#    print 'test1'
2379    GetFobsSq(Histograms,Phases,parmDict,calcControls)
2380#    print 'test2'
2381    sigDict = dict(zip(varyList,sig))
2382    newCellDict = GetNewCellParms(parmDict,varyList)
2383    newAtomDict = ApplyXYZshifts(parmDict,varyList)
2384    covData = {'variables':result[0],'varyList':varyList,'sig':sig,
2385        'covMatrix':covMatrix,'title':GPXfile,'newAtomDict':newAtomDict,'newCellDict':newCellDict}
2386    # add the uncertainties into the esd dictionary (sigDict)
2387    sigDict.update(G2mv.ComputeDepESD(covMatrix,varyList,parmDict))
2388    SetPhaseData(parmDict,sigDict,Phases,covData)
2389    SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms)
2390    SetHistogramData(parmDict,sigDict,Histograms)
2391    G2mv.PrintIndependentVars(parmDict,varyList,sigDict)
2392    G2IO.SetUsedHistogramsAndPhases(GPXfile,Histograms,Phases,covData)
2393   
2394#for testing purposes!!!
2395#    import cPickle
2396#    file = open('structTestdata.dat','wb')
2397#    cPickle.dump(parmDict,file,1)
2398#    cPickle.dump(varyList,file,1)
2399#    for histogram in Histograms:
2400#        if 'PWDR' in histogram[:4]:
2401#            Histogram = Histograms[histogram]
2402#    cPickle.dump(Histogram,file,1)
2403#    cPickle.dump(Phases,file,1)
2404#    cPickle.dump(calcControls,file,1)
2405#    cPickle.dump(pawleyLookup,file,1)
2406#    file.close()
2407
2408    if dlg:
2409        return Rwp
2410
2411def SeqRefine(GPXfile,dlg):
2412    import pytexture as ptx
2413    ptx.pyqlmninit()            #initialize fortran arrays for spherical harmonics
2414   
2415    ShowBanner()
2416    print ' Sequential Refinement'
2417    G2mv.InitVars()   
2418    Controls = G2IO.GetControls(GPXfile)
2419    ShowControls(Controls)           
2420    constrDict,constrFlag,fixedList = G2IO.GetConstraints(GPXfile)
2421    Histograms,Phases = G2IO.GetUsedHistogramsAndPhases(GPXfile)
2422    if not Phases:
2423        print ' *** ERROR - you have no histograms to refine! ***'
2424        print ' *** Refine aborted ***'
2425        raise Exception
2426    if not Histograms:
2427        print ' *** ERROR - you have no data to refine with! ***'
2428        print ' *** Refine aborted ***'
2429        raise Exception
2430    Natoms,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables = GetPhaseData(Phases,False)
2431    if 'Seq Data' in Controls:
2432        histNames = Controls['Seq Data']
2433    else:
2434        histNames = G2IO.GetHistogramNames(GPXfile,['PWDR',])
2435    if 'Reverse Seq' in Controls:
2436        if Controls['Reverse Seq']:
2437            histNames.reverse()
2438    SeqResult = {'histNames':histNames}
2439    makeBack = True
2440    for ihst,histogram in enumerate(histNames):
2441        ifPrint = False
2442        if dlg:
2443            dlg.SetTitle('Residual for histogram '+str(ihst))
2444        calcControls = {}
2445        calcControls['Natoms'] = Natoms
2446        calcControls['FFtables'] = FFtables
2447        calcControls['BLtables'] = BLtables
2448        varyList = []
2449        parmDict = {}
2450        Histo = {histogram:Histograms[histogram],}
2451        hapVary,hapDict,controlDict = GetHistogramPhaseData(Phases,Histo,False)
2452        calcControls.update(controlDict)
2453        histVary,histDict,controlDict = GetHistogramData(Histo,False)
2454        calcControls.update(controlDict)
2455        varyList = phaseVary+hapVary+histVary
2456        if not ihst:
2457            saveVaryList = varyList[:]
2458            for i,item in enumerate(saveVaryList):
2459                items = item.split(':')
2460                if items[1]:
2461                    items[1] = ''
2462                item = ':'.join(items)
2463                saveVaryList[i] = item
2464            SeqResult['varyList'] = saveVaryList
2465        else:
2466            newVaryList = varyList[:]
2467            for i,item in enumerate(newVaryList):
2468                items = item.split(':')
2469                if items[1]:
2470                    items[1] = ''
2471                item = ':'.join(items)
2472                newVaryList[i] = item
2473            if newVaryList != SeqResult['varyList']:
2474                print newVaryList
2475                print SeqResult['varyList']
2476                print '**** ERROR - variable list for this histogram does not match previous'
2477                raise Exception
2478        parmDict.update(phaseDict)
2479        parmDict.update(hapDict)
2480        parmDict.update(histDict)
2481        GetFprime(calcControls,Histo)
2482        constrDict,constrFlag,fixedList = G2mv.InputParse([])        #constraints go here?
2483        groups,parmlist = G2mv.GroupConstraints(constrDict)
2484        G2mv.GenerateConstraints(groups,parmlist,varyList,constrDict,constrFlag,fixedList)
2485        G2mv.Map2Dict(parmDict,varyList)
2486   
2487        while True:
2488            begin = time.time()
2489            values =  np.array(Dict2Values(parmDict, varyList))
2490            Ftol = Controls['min dM/M']
2491            Factor = Controls['shift factor']
2492            maxCyc = Controls['max cyc']
2493
2494            if 'Jacobian' in Controls['deriv type']:           
2495                result = so.leastsq(errRefine,values,Dfun=dervRefine,full_output=True,
2496                    ftol=Ftol,col_deriv=True,factor=Factor,
2497                    args=([Histograms,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg))
2498                ncyc = int(result[2]['nfev']/2)
2499            elif 'Hessian' in Controls['deriv type']:
2500                result = G2mth.HessianLSQ(errRefine,values,Hess=HessRefine,ftol=Ftol,maxcyc=maxCyc,
2501                    args=([Histograms,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg))
2502                ncyc = result[2]['num cyc']+1                           
2503            else:           #'numeric'
2504                result = so.leastsq(errRefine,values,full_output=True,ftol=Ftol,epsfcn=1.e-8,factor=Factor,
2505                    args=([Histograms,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg))
2506                ncyc = int(result[2]['nfev']/len(varyList))
2507
2508
2509
2510            runtime = time.time()-begin
2511            chisq = np.sum(result[2]['fvec']**2)
2512            Values2Dict(parmDict, varyList, result[0])
2513            G2mv.Dict2Map(parmDict,varyList)
2514           
2515            Rwp = np.sqrt(chisq/Histo['sumwYo'])*100.      #to %
2516            GOF = chisq/(Histo['Nobs']-len(varyList))
2517            print '\n Refinement results for histogram: v'+histogram
2518            print 135*'-'
2519            print ' Number of function calls:',result[2]['nfev'],' Number of observations: ',Histo['Nobs'],' Number of parameters: ',len(varyList)
2520            print ' Refinement time = %8.3fs, %8.3fs/cycle, for %d cycles'%(runtime,runtime/ncyc,ncyc)
2521            print ' wRp = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f'%(Rwp,chisq,GOF)
2522            try:
2523                covMatrix = result[1]*GOF
2524                sig = np.sqrt(np.diag(covMatrix))
2525                if np.any(np.isnan(sig)):
2526                    print '*** Least squares aborted - some invalid esds possible ***'
2527                    ifPrint = True
2528                break                   #refinement succeeded - finish up!
2529            except TypeError:          #result[1] is None on singular matrix
2530                print '**** Refinement failed - singular matrix ****'
2531                if 'Hessian' in Controls['deriv type']:
2532                    for i in result[2]['psing'].reverse():
2533                            print 'Removing parameter: ',varyList[i]
2534                            del(varyList[i])                   
2535                else:
2536                    Ipvt = result[2]['ipvt']
2537                    for i,ipvt in enumerate(Ipvt):
2538                        if not np.sum(result[2]['fjac'],axis=1)[i]:
2539                            print 'Removing parameter: ',varyList[ipvt-1]
2540                            del(varyList[ipvt-1])
2541                            break
2542   
2543        GetFobsSq(Histo,Phases,parmDict,calcControls)
2544        sigDict = dict(zip(varyList,sig))
2545        newCellDict = GetNewCellParms(parmDict,varyList)
2546        newAtomDict = ApplyXYZshifts(parmDict,varyList)
2547        covData = {'variables':result[0],'varyList':varyList,'sig':sig,
2548            'covMatrix':covMatrix,'title':histogram,'newAtomDict':newAtomDict,'newCellDict':newCellDict}
2549        SetHistogramPhaseData(parmDict,sigDict,Phases,Histo,ifPrint)
2550        SetHistogramData(parmDict,sigDict,Histo,ifPrint)
2551        SeqResult[histogram] = covData
2552        G2IO.SetUsedHistogramsAndPhases(GPXfile,Histo,Phases,covData,makeBack)
2553        makeBack = False
2554    G2IO.SetSeqResult(GPXfile,Histograms,SeqResult)
2555
2556def DistAngle(DisAglData):
2557    import numpy.ma as ma
2558   
2559    def ShowBanner(name):
2560        print 80*'*'
2561        print '   Interatomic Distances and Angles for phase '+name
2562        print 80*'*','\n'
2563
2564    ShowBanner(DisAglData['Name'])
2565    SGData = DisAglData['SGData']
2566    SGtext = G2spc.SGPrint(SGData)
2567    for line in SGtext: print line
2568    Cell = DisAglData['Cell']
2569   
2570    Amat,Bmat = G2lat.cell2AB(Cell[:6])
2571    if 'covData' in DisAglData:   
2572        covData = DisAglData['covData']
2573        pfx = str(DisAglData['pId'])+'::'
2574        A = G2lat.cell2A(Cell[:6])
2575        cellSig = getCellEsd(pfx,SGData,A,covData)
2576        names = [' a = ',' b = ',' c = ',' alpha = ',' beta = ',' gamma = ',' Volume = ']
2577        valEsd = [G2mth.ValEsd(Cell[i],cellSig[i],True) for i in range(7)]
2578        line = '\n Unit cell:'
2579        for name,vals in zip(names,valEsd):
2580            line += name+vals 
2581        print line
2582    else:
2583        print '\n Unit cell: a = ','%.5f'%(Cell[0]),' b = ','%.5f'%(Cell[1]),' c = ','%.5f'%(Cell[2]), \
2584            ' alpha = ','%.3f'%(Cell[3]),' beta = ','%.3f'%(Cell[4]),' gamma = ', \
2585            '%.3f'%(Cell[5]),' volume = ','%.3f'%(Cell[6])
2586    Factor = DisAglData['Factors']
2587    Radii = dict(zip(DisAglData['AtomTypes'],zip(DisAglData['BondRadii'],DisAglData['AngleRadii'])))
2588    Units = np.array([                   #is there a nicer way to make this?
2589        [-1,-1,-1],[-1,-1,0],[-1,-1,1],[-1,0,-1],[-1,0,0],[-1,0,1],[-1,1,-1],[-1,1,0],[-1,1,1],
2590        [0,-1,-1],[0,-1,0],[0,-1,1],[0,0,-1],[0,0,0],[0,0,1],[0,1,-1],[0,1,0],[0,1,1],
2591        [1,-1,-1],[1,-1,0],[1,-1,1],[1,0,-1],[1,0,0],[1,0,1],[1,1,-1],[1,1,0],[1,1,1]])
2592    origAtoms = zip(DisAglData['OrigIndx'],DisAglData['OrigAtoms'])
2593    targAtoms = DisAglData['TargAtoms']
2594    for Oatom in origAtoms:
2595        Dist = []
2596        Vect = []
2597        for Tatom in targAtoms:
2598            result = G2spc.GenAtom(Tatom[2:5],SGData,False)
2599            BsumR = (Radii[Oatom[1][1]][0]+Radii[Tatom[1]][0])*Factor[0]
2600            AsumR = (Radii[Oatom[1][1]][1]+Radii[Tatom[1]][1])*Factor[1]
2601            for Txyz,Top,Tunit in result:
2602                Dx = (Txyz-np.array(Oatom[1][2:5]))+Units
2603                dx = np.inner(Amat,Dx)
2604                dist = ma.masked_less(np.sqrt(np.sum(dx**2,axis=0)),0.5)
2605                IndB = ma.nonzero(ma.masked_greater(dist-BsumR,0.))
2606                IndA = ma.nonzero(ma.masked_greater(dist-AsumR,0.))
2607                if np.any(IndB):
2608                    for indb in IndB:
2609                        Dist.append([Oatom[1][0],Tatom[0],Tunit,Top,ma.getdata(dist[indb])])
2610                        Vect += dx.T[indb]
2611       
2612        for dist in Dist:
2613            print dist
2614           
2615        for vect in Vect:
2616            print vect
2617           
2618                 
2619
2620#    def FindBonds():                    #uses numpy & masks - very fast even for proteins!
2621#        import numpy.ma as ma
2622#        atomData = data['Atoms']
2623#        generalData = data['General']
2624#        Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
2625#        radii = generalData['BondRadii']
2626#        atomTypes = generalData['AtomTypes']
2627#        try:
2628#            indH = atomTypes.index('H')
2629#            radii[indH] = 0.5
2630#        except:
2631#            pass           
2632#        for atom in atomData:
2633#            atom[-2] = []               #clear out old bond/angle tables
2634#            atom[-1] = []
2635#        Indx = range(len(atomData))
2636#        Atoms = []
2637#        Radii = []
2638#        for atom in atomData:
2639#            Atoms.append(np.array(atom[cx:cx+3]))
2640#            try:
2641#                if not hydro and atom[ct] == 'H':
2642#                    Radii.append(0.0)
2643#                else:
2644#                    Radii.append(radii[atomTypes.index(atom[ct])])
2645#            except ValueError:          #changed atom type!
2646#                Radii.append(0.20)
2647#        Atoms = np.array(Atoms)
2648#        Radii = np.array(Radii)
2649#        IASR = zip(Indx,Atoms,Radii)
2650#        for atomA in IASR:
2651#                Dx = Atoms-atomA[1]
2652#                dist = ma.masked_less(np.sqrt(np.sum(np.inner(Amat,Dx)**2,axis=0)),0.5) #gets rid of self & disorder "bonds" < 0.5A
2653#                sumR = atomA[3]+Radii
2654#                IndB = ma.nonzero(ma.masked_greater(dist-radiusFactor*sumR,0.))                 #get indices of bonded atoms
2655#                i = atomA[0]
2656#                for j in IndB[0]:
2657#                    if j > i:
2658#                        atomData[i][-2].append(Dx[j]*Radii[i]/sumR[j])
2659#                        atomData[j][-2].append(-Dx[j]*Radii[j]/sumR[j])
2660#                       
2661
2662def main():
2663    arg = sys.argv
2664    if len(arg) > 1:
2665        GPXfile = arg[1]
2666        if not ospath.exists(GPXfile):
2667            print 'ERROR - ',GPXfile," doesn't exist!"
2668            exit()
2669        GPXpath = ospath.dirname(arg[1])
2670        Refine(GPXfile,None)
2671    else:
2672        print 'ERROR - missing filename'
2673        exit()
2674    print "Done"
2675         
2676if __name__ == '__main__':
2677    main()
Note: See TracBrowser for help on using the repository browser.