source: trunk/GSASIIstruct.py @ 411

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

fix very old (!) bug in psvfcj.for
implement neutron scattering lengths in GSASII including the dozen anomalous scattering isotopes
isotope choice is in General
so GSASII will now do neutron CW Rietveld refinements
some cleanup of the constraints GUI
remove varyList from GSASIImapvars.py globals

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