source: trunk/GSASIIstruct.py @ 395

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

Add goniometer omega, chi & phi to sample data
put SH texture in General
fix phase delete to remove it from reflection lists as well
continue development of constraints/restraints GUI
fixes to texture computations, GUI & least squares refinement

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