source: trunk/GSASIIstruct.py @ 480

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

missing import in G2struct for ospath

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