source: MPbranch/GSASIIstruct.py @ 507

Last change on this file since 507 was 507, checked in by toby, 11 years ago

fix function calc to allow for out of sequence finish of 'threads'; other minor fixes

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