source: trunk/GSASIIstruct.py @ 391

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

set default atom frac=1 for cif files
store covMatrix not normalized by diagonal
make esds on cell parameters

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