source: trunk/GSASIIstruct.py @ 400

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

new default for Vcov contour plot - RdYlGn?
faster cleanup on changing/reloading projects
cleanup data delete
implement sample parameter copy
improve Vcov plotting routine
implement plot of vcov from seq refinements

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