source: trunk/GSASIIstruct.py @ 403

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

speedup GetUsedHistogramsAndPhases?

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