source: trunk/GSASIIstruct.py @ 489

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

fix to speed up scat fac calculations & skip texture calcs if order=0

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