source: branch/MPbranch/GSASIIstruct.py

Last change on this file was 518, checked in by toby, 10 years ago

oops, remove single processor debug statement

  • Property svn:keywords set to Date Author Revision URL Id
File size: 147.9 KB
Line 
1#GSASIIstructure - structure computation routines
2########### SVN repository information ###################
3# $Date: 2012-03-09 21:26:45 +0000 (Fri, 09 Mar 2012) $
4# $Author: toby $
5# $Revision: 518 $
6# $URL: branch/MPbranch/GSASIIstruct.py $
7# $Id: GSASIIstruct.py 518 2012-03-09 21:26:45Z 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    return Icorr       
2004   
2005def GetIntensityDerv(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict):
2006    dIdsh = 1./parmDict[hfx+'Scale']
2007    dIdsp = 1./parmDict[phfx+'Scale']
2008    if 'X' in parmDict[hfx+'Type']:
2009        pola,dIdPola = G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5],parmDict[hfx+'Azimuth'])
2010        dIdPola /= pola
2011    else:       #'N'
2012        dIdPola = 0.0
2013    POcorr,dIdPO = GetPrefOriDerv(refl,G,g,phfx,hfx,SGData,calcControls,parmDict)
2014    for iPO in dIdPO:
2015        dIdPO[iPO] /= POcorr
2016    dFdODF = {}
2017    dFdSA = [0,0,0]
2018    if pfx+'SHorder' in parmDict:
2019        odfCor,dFdODF,dFdSA = SHTXcalDerv(refl,g,pfx,hfx,SGData,calcControls,parmDict)
2020        for iSH in dFdODF:
2021            dFdODF[iSH] /= odfCor
2022        for i in range(3):
2023            dFdSA[i] /= odfCor
2024    return dIdsh,dIdsp,dIdPola,dIdPO,dFdODF,dFdSA
2025       
2026def GetSampleGam(refl,wave,G,GB,phfx,calcControls,parmDict):
2027    costh = cosd(refl[5]/2.)
2028    #crystallite size
2029    if calcControls[phfx+'SizeType'] == 'isotropic':
2030        gam = 1.8*wave/(np.pi*parmDict[phfx+'Size:i']*costh)
2031    elif calcControls[phfx+'SizeType'] == 'uniaxial':
2032        H = np.array(refl[:3])
2033        P = np.array(calcControls[phfx+'SizeAxis'])
2034        cosP,sinP = G2lat.CosSinAngle(H,P,G)
2035        gam = (1.8*wave/np.pi)/(parmDict[phfx+'Size:i']*parmDict[phfx+'Size:a']*costh)
2036        gam *= np.sqrt((sinP*parmDict[phfx+'Size:a'])**2+(cosP*parmDict[phfx+'Size:i'])**2)
2037    else:           #ellipsoidal crystallites
2038        Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
2039        H = np.array(refl[:3])
2040        lenR = G2pwd.ellipseSize(H,Sij,GB)
2041        gam = 1.8*wave/(np.pi*costh*lenR)
2042    #microstrain               
2043    if calcControls[phfx+'MustrainType'] == 'isotropic':
2044        gam += 0.018*parmDict[phfx+'Mustrain:i']*tand(refl[5]/2.)/np.pi
2045    elif calcControls[phfx+'MustrainType'] == 'uniaxial':
2046        H = np.array(refl[:3])
2047        P = np.array(calcControls[phfx+'MustrainAxis'])
2048        cosP,sinP = G2lat.CosSinAngle(H,P,G)
2049        Si = parmDict[phfx+'Mustrain:i']
2050        Sa = parmDict[phfx+'Mustrain:a']
2051        gam += 0.018*Si*Sa*tand(refl[5]/2.)/(np.pi*np.sqrt((Si*cosP)**2+(Sa*sinP)**2))
2052    else:       #generalized - P.W. Stephens model
2053        pwrs = calcControls[phfx+'MuPwrs']
2054        sum = 0
2055        for i,pwr in enumerate(pwrs):
2056            sum += parmDict[phfx+'Mustrain:'+str(i)]*refl[0]**pwr[0]*refl[1]**pwr[1]*refl[2]**pwr[2]
2057        gam += 0.018*refl[4]**2*tand(refl[5]/2.)*sum           
2058    return gam
2059       
2060def GetSampleGamDerv(refl,wave,G,GB,phfx,calcControls,parmDict):
2061    gamDict = {}
2062    costh = cosd(refl[5]/2.)
2063    tanth = tand(refl[5]/2.)
2064    #crystallite size derivatives
2065    if calcControls[phfx+'SizeType'] == 'isotropic':
2066        gamDict[phfx+'Size:i'] = -1.80*wave/(np.pi*costh)
2067    elif calcControls[phfx+'SizeType'] == 'uniaxial':
2068        H = np.array(refl[:3])
2069        P = np.array(calcControls[phfx+'SizeAxis'])
2070        cosP,sinP = G2lat.CosSinAngle(H,P,G)
2071        Si = parmDict[phfx+'Size:i']
2072        Sa = parmDict[phfx+'Size:a']
2073        gami = (1.8*wave/np.pi)/(Si*Sa)
2074        sqtrm = np.sqrt((sinP*Sa)**2+(cosP*Si)**2)
2075        gam = gami*sqtrm/costh           
2076        gamDict[phfx+'Size:i'] = gami*Si*cosP**2/(sqtrm*costh)-gam/Si
2077        gamDict[phfx+'Size:a'] = gami*Sa*sinP**2/(sqtrm*costh)-gam/Sa         
2078    else:           #ellipsoidal crystallites
2079        const = 1.8*wave/(np.pi*costh)
2080        Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
2081        H = np.array(refl[:3])
2082        R,dRdS = G2pwd.ellipseSizeDerv(H,Sij,GB)
2083        for i,item in enumerate([phfx+'Size:%d'%(j) for j in range(6)]):
2084            gamDict[item] = -(const/R**2)*dRdS[i]
2085    #microstrain derivatives               
2086    if calcControls[phfx+'MustrainType'] == 'isotropic':
2087        gamDict[phfx+'Mustrain:i'] =  0.018*tanth/np.pi           
2088    elif calcControls[phfx+'MustrainType'] == 'uniaxial':
2089        H = np.array(refl[:3])
2090        P = np.array(calcControls[phfx+'MustrainAxis'])
2091        cosP,sinP = G2lat.CosSinAngle(H,P,G)
2092        Si = parmDict[phfx+'Mustrain:i']
2093        Sa = parmDict[phfx+'Mustrain:a']
2094        gami = 0.018*Si*Sa*tanth/np.pi
2095        sqtrm = np.sqrt((Si*cosP)**2+(Sa*sinP)**2)
2096        gam = gami/sqtrm
2097        gamDict[phfx+'Mustrain:i'] = gam/Si-gami*Si*cosP**2/sqtrm**3
2098        gamDict[phfx+'Mustrain:a'] = gam/Sa-gami*Sa*sinP**2/sqtrm**3
2099    else:       #generalized - P.W. Stephens model
2100        pwrs = calcControls[phfx+'MuPwrs']
2101        const = 0.018*refl[4]**2*tanth
2102        for i,pwr in enumerate(pwrs):
2103            gamDict[phfx+'Mustrain:'+str(i)] = const*refl[0]**pwr[0]*refl[1]**pwr[1]*refl[2]**pwr[2]
2104    return gamDict
2105       
2106def GetReflPos(refl,wave,G,hfx,calcControls,parmDict):
2107    h,k,l = refl[:3]
2108    dsq = 1./G2lat.calc_rDsq2(np.array([h,k,l]),G)
2109    d = np.sqrt(dsq)
2110
2111    refl[4] = d
2112    pos = 2.0*asind(wave/(2.0*d))+parmDict[hfx+'Zero']
2113    const = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius'])                  #shifts in microns
2114    if 'Bragg' in calcControls[hfx+'instType']:
2115        pos -= const*(4.*parmDict[hfx+'Shift']*cosd(pos/2.0)+ \
2116            parmDict[hfx+'Transparency']*sind(pos)*100.0)            #trans(=1/mueff) in cm
2117    else:               #Debye-Scherrer - simple but maybe not right
2118        pos -= const*(parmDict[hfx+'DisplaceX']*cosd(pos)+parmDict[hfx+'DisplaceY']*sind(pos))
2119    return pos
2120
2121def GetReflPosDerv(refl,wave,A,hfx,calcControls,parmDict):
2122    dpr = 180./np.pi
2123    h,k,l = refl[:3]
2124    dstsq = G2lat.calc_rDsq(np.array([h,k,l]),A)
2125    dst = np.sqrt(dstsq)
2126    pos = refl[5]-parmDict[hfx+'Zero']
2127    const = dpr/np.sqrt(1.0-wave**2*dstsq/4.0)
2128    dpdw = const*dst
2129    dpdA = np.array([h**2,k**2,l**2,h*k,h*l,k*l])
2130    dpdA *= const*wave/(2.0*dst)
2131    dpdZ = 1.0
2132    const = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius'])                  #shifts in microns
2133    if 'Bragg' in calcControls[hfx+'instType']:
2134        dpdSh = -4.*const*cosd(pos/2.0)
2135        dpdTr = -const*sind(pos)*100.0
2136        return dpdA,dpdw,dpdZ,dpdSh,dpdTr,0.,0.
2137    else:               #Debye-Scherrer - simple but maybe not right
2138        dpdXd = -const*cosd(pos)
2139        dpdYd = -const*sind(pos)
2140        return dpdA,dpdw,dpdZ,0.,0.,dpdXd,dpdYd
2141           
2142def GetHStrainShift(refl,SGData,phfx,parmDict):
2143    laue = SGData['SGLaue']
2144    uniq = SGData['SGUniq']
2145    h,k,l = refl[:3]
2146    if laue in ['m3','m3m']:
2147        Dij = parmDict[phfx+'D11']*(h**2+k**2+l**2)+ \
2148            refl[4]**2*parmDict[phfx+'eA']*((h*k)**2+(h*l)**2+(k*l)**2)/(h**2+k**2+l**2)**2
2149    elif laue in ['6/m','6/mmm','3m1','31m','3']:
2150        Dij = parmDict[phfx+'D11']*(h**2+k**2+h*k)+parmDict[phfx+'D33']*l**2
2151    elif laue in ['3R','3mR']:
2152        Dij = parmDict[phfx+'D11']*(h**2+k**2+l**2)+parmDict[phfx+'D12']*(h*k+h*l+k*l)
2153    elif laue in ['4/m','4/mmm']:
2154        Dij = parmDict[phfx+'D11']*(h**2+k**2)+parmDict[phfx+'D33']*l**2
2155    elif laue in ['mmm']:
2156        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2
2157    elif laue in ['2/m']:
2158        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2
2159        if uniq == 'a':
2160            Dij += parmDict[phfx+'D23']*k*l
2161        elif uniq == 'b':
2162            Dij += parmDict[phfx+'D13']*h*l
2163        elif uniq == 'c':
2164            Dij += parmDict[phfx+'D12']*h*k
2165    else:
2166        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2+ \
2167            parmDict[phfx+'D12']*h*k+parmDict[phfx+'D13']*h*l+parmDict[phfx+'D23']*k*l
2168    return Dij*refl[4]**2*tand(refl[5]/2.0)
2169           
2170def GetHStrainShiftDerv(refl,SGData,phfx):
2171    laue = SGData['SGLaue']
2172    uniq = SGData['SGUniq']
2173    h,k,l = refl[:3]
2174    if laue in ['m3','m3m']:
2175        dDijDict = {phfx+'D11':h**2+k**2+l**2,
2176            phfx+'eA':((h*k)**2+(h*l)**2+(k*l)**2)/(h**2+k**2+l**2)**2}
2177    elif laue in ['6/m','6/mmm','3m1','31m','3']:
2178        dDijDict = {phfx+'D11':h**2+k**2+h*k,phfx+'D33':l**2}
2179    elif laue in ['3R','3mR']:
2180        dDijDict = {phfx+'D11':h**2+k**2+l**2,phfx+'D12':h*k+h*l+k*l}
2181    elif laue in ['4/m','4/mmm']:
2182        dDijDict = {phfx+'D11':h**2+k**2,phfx+'D33':l**2}
2183    elif laue in ['mmm']:
2184        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2}
2185    elif laue in ['2/m']:
2186        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2}
2187        if uniq == 'a':
2188            dDijDict[phfx+'D23'] = k*l
2189        elif uniq == 'b':
2190            dDijDict[phfx+'D13'] = h*l
2191        elif uniq == 'c':
2192            dDijDict[phfx+'D12'] = h*k
2193            names.append()
2194    else:
2195        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2,
2196            phfx+'D12':h*k,phfx+'D13':h*l,phfx+'D23':k*l}
2197    for item in dDijDict:
2198        dDijDict[item] *= refl[4]**2*tand(refl[5]/2.0)
2199    return dDijDict
2200   
2201def GetFprime(controlDict,Histograms):
2202    FFtables = controlDict['FFtables']
2203    if not FFtables:
2204        return
2205    histoList = Histograms.keys()
2206    histoList.sort()
2207    for histogram in histoList:
2208        if 'PWDR' in histogram[:4]:
2209            Histogram = Histograms[histogram]
2210            hId = Histogram['hId']
2211            hfx = ':%d:'%(hId)
2212            keV = controlDict[hfx+'keV']
2213            for El in FFtables:
2214                Orbs = G2el.GetXsectionCoeff(El.split('+')[0].split('-')[0])
2215                FP,FPP,Mu = G2el.FPcalc(Orbs, keV)
2216                FFtables[El][hfx+'FP'] = FP
2217                FFtables[El][hfx+'FPP'] = FPP               
2218           
2219def GetFobsSq(Histograms,Phases,parmDict,calcControls):
2220    histoList = Histograms.keys()
2221    histoList.sort()
2222    for histogram in histoList:
2223        if 'PWDR' in histogram[:4]:
2224            Histogram = Histograms[histogram]
2225            hId = Histogram['hId']
2226            hfx = ':%d:'%(hId)
2227            Limits = calcControls[hfx+'Limits']
2228            shl = max(parmDict[hfx+'SH/L'],0.002)
2229            Ka2 = False
2230            kRatio = 0.0
2231            if hfx+'Lam1' in parmDict.keys():
2232                Ka2 = True
2233                lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
2234                kRatio = parmDict[hfx+'I(L2)/I(L1)']
2235            x,y,w,yc,yb,yd = Histogram['Data']
2236            ymb = np.array(y-yb)
2237            ycmb = np.array(yc-yb)
2238            ratio = np.where(ycmb!=0.,ymb/ycmb,0.0)           
2239            xB = np.searchsorted(x,Limits[0])
2240            xF = np.searchsorted(x,Limits[1])
2241            refLists = Histogram['Reflection Lists']
2242            for phase in refLists:
2243                Phase = Phases[phase]
2244                pId = Phase['pId']
2245                phfx = '%d:%d:'%(pId,hId)
2246                refList = refLists[phase]
2247                sumFo = 0.0
2248                sumdF = 0.0
2249                sumFosq = 0.0
2250                sumdFsq = 0.0
2251                for refl in refList:
2252                    if 'C' in calcControls[hfx+'histType']:
2253                        yp = np.zeros_like(yb)
2254                        Wd,fmin,fmax = G2pwd.getWidths(refl[5],refl[6],refl[7],shl)
2255                        iBeg = np.searchsorted(x[xB:xF],refl[5]-fmin)
2256                        iFin = np.searchsorted(x[xB:xF],refl[5]+fmax)
2257                        iFin2 = iFin
2258                        yp[iBeg:iFin] = refl[13]*refl[9]*G2pwd.getFCJVoigt3(refl[5],refl[6],refl[7],shl,x[iBeg:iFin])    #>90% of time spent here                           
2259                        if Ka2:
2260                            pos2 = refl[5]+lamRatio*tand(refl[5]/2.0)       # + 360/pi * Dlam/lam * tan(th)
2261                            Wd,fmin,fmax = G2pwd.getWidths(pos2,refl[6],refl[7],shl)
2262                            iBeg2 = np.searchsorted(x,pos2-fmin)
2263                            iFin2 = np.searchsorted(x,pos2+fmax)
2264                            yp[iBeg2:iFin2] += refl[13]*refl[9]*kRatio*G2pwd.getFCJVoigt3(pos2,refl[6],refl[7],shl,x[iBeg2:iFin2])        #and here
2265                        refl[8] = np.sum(np.where(ratio[iBeg:iFin2]>0.,yp[iBeg:iFin2]*ratio[iBeg:iFin2]/(refl[13]*(1.+kRatio)),0.0))
2266                    elif 'T' in calcControls[hfx+'histType']:
2267                        print 'TOF Undefined at present'
2268                        raise Exception    #no TOF yet
2269                    Fo = np.sqrt(np.abs(refl[8]))
2270                    Fc = np.sqrt(np.abs(refl[9]))
2271                    sumFo += Fo
2272                    sumFosq += refl[8]**2
2273                    sumdF += np.abs(Fo-Fc)
2274                    sumdFsq += (refl[8]-refl[9])**2
2275                Histogram[phfx+'Rf'] = min(100.,(sumdF/sumFo)*100.)
2276                Histogram[phfx+'Rf^2'] = min(100.,np.sqrt(sumdFsq/sumFosq)*100.)
2277                Histogram[phfx+'Nref'] = len(refList)
2278               
2279def ComputeReflectionProfile(arg):
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    (RefNum,refList,x,Phase,calcControls,wave,G,g,GB,Vst,
2294     parmDict,pfx,phfx,hfx,Ka2,lamRatio,kRatio) = arg
2295    yc = np.zeros_like(x)
2296    iBegO = None
2297    SGData = Phase['General']['SGData']
2298    shl = max(parmDict[hfx+'SH/L'],0.002)
2299    if 'C' in calcControls[hfx+'histType']:
2300        for refl in refList:
2301            refl[5] = GetReflPos(refl,wave,G,hfx,calcControls,parmDict)         #corrected reflection position
2302            Lorenz = 1./(2.*sind(refl[5]/2.)**2*cosd(refl[5]/2.))           #Lorentz correction
2303            refl[5] += GetHStrainShift(refl,SGData,phfx,parmDict)               #apply hydrostatic strain shift
2304            refl[6:8] = GetReflSIgGam(refl,wave,G,GB,hfx,phfx,calcControls,parmDict)    #peak sig & gam
2305            refl[13] = Vst*Lorenz * GetIntensityCorr(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)    #puts corrections in refl[13]
2306            Wd,fmin,fmax = G2pwd.getWidths(refl[5],refl[6],refl[7],shl)
2307            iBeg = np.searchsorted(x,refl[5]-fmin)
2308            iFin = np.searchsorted(x,refl[5]+fmax)
2309            if not iBeg+iFin:       #peak below low limit - skip peak
2310                continue
2311            elif not iBeg-iFin:     #peak above high limit - done
2312                break
2313            yc[iBeg:iFin] += refl[13]*refl[9]*G2pwd.getFCJVoigt3(refl[5],refl[6],refl[7],shl,x[iBeg:iFin])    #>90% of time spent here
2314            if iBegO is None:
2315                iBegO = iBeg
2316                iFinO = iFin
2317            else:
2318                iBegO = min(iBegO,iBeg)
2319                iFinO = max(iFinO,iFin)
2320            if Ka2:
2321                pos2 = refl[5]+lamRatio*tand(refl[5]/2.0)       # + 360/pi * Dlam/lam * tan(th)
2322                Wd,fmin,fmax = G2pwd.getWidths(pos2,refl[6],refl[7],shl)
2323                iBeg2 = np.searchsorted(x,pos2-fmin)
2324                iFin2 = np.searchsorted(x,pos2+fmax)
2325                if not iBeg2+iFin2:       #peak below low limit - skip peak
2326                    continue
2327                elif not iBeg2-iFin2:     #peak above high limit - done
2328                    continue
2329                yc[iBeg2:iFin2] += refl[13]*refl[9]*kRatio*G2pwd.getFCJVoigt3(pos2,refl[6],refl[7],shl,x[iBeg:iFin])        #and here
2330                iBegO = min(iBegO,iBeg2)
2331                iFinO = max(iFinO,iFin2)
2332    elif 'T' in calcControls[hfx+'histType']:
2333        print 'TOF Undefined at present'
2334        raise Exception    #no TOF yet
2335
2336    if iBegO is None: return None # no reflections were computed
2337    return RefNum,refList,iBegO,iFinO,yc[iBegO:iFinO]
2338
2339def getPowderProfile(parmDict,x,varylist,Histogram,Phases,calcControls,pawleyLookup):               
2340    if calcControls['max Hprocess'] <= 1 and calcControls['max Rprocess'] > 1:
2341        mpPool = mp.Pool()
2342        #print 'multiprocess function call'
2343    #else:
2344        #print 'single thread function call'
2345    hId = Histogram['hId']
2346    hfx = ':%d:'%(hId)
2347    bakType = calcControls[hfx+'bakType']
2348    yb = G2pwd.getBackground(hfx,parmDict,bakType,x)
2349    yc = np.zeros_like(yb)
2350       
2351    if 'C' in calcControls[hfx+'histType']:   
2352        shl = max(parmDict[hfx+'SH/L'],0.002)
2353        if hfx+'Lam1' in parmDict.keys():
2354            wave = parmDict[hfx+'Lam1']
2355            Ka2 = True
2356            lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
2357            kRatio = parmDict[hfx+'I(L2)/I(L1)']
2358        else:
2359            wave = parmDict[hfx+'Lam']
2360            Ka2 = False
2361            lamRatio = None
2362            kRatio = None
2363    else:
2364        print 'TOF Undefined at present'
2365        raise ValueError
2366    for phase in Histogram['Reflection Lists']:
2367        refList = Histogram['Reflection Lists'][phase]
2368        Phase = Phases[phase]
2369        pId = Phase['pId']
2370        pfx = '%d::'%(pId)
2371        phfx = '%d:%d:'%(pId,hId)
2372        hfx = ':%d:'%(hId)
2373        SGData = Phase['General']['SGData']
2374        A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
2375        G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
2376        GA,GB = G2lat.Gmat2AB(G)    #Orthogonalization matricies
2377        Vst = np.sqrt(nl.det(G))    #V*
2378        if 'Pawley' not in Phase['General']['Type']:
2379            refList = StructureFactor(refList,G,hfx,pfx,SGData,calcControls,parmDict)
2380        else:
2381            for refl in refList:
2382                h,k,l = refl[:3]
2383                try:
2384                    refl[9] = abs(parmDict[pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])])
2385                except KeyError:
2386                    #print ' ***Error %d,%d,%d missing from Pawley reflection list ***'%(h,k,l)
2387                    continue
2388
2389        chunk = 200
2390        if (calcControls['max Hprocess'] > 1 or
2391            calcControls['max Rprocess'] == 1 or
2392            mp.cpu_count() < 2 or len(refList) <  2 * chunk):
2393            res = ComputeReflectionProfile(
2394                (0,refList,x,Phase,calcControls,wave,G,g,GB,Vst,
2395                 parmDict,pfx,phfx,hfx,Ka2,lamRatio,kRatio))
2396            if not res: continue
2397            RefNum,refls,iBeg,iFin,yct = res
2398            yc[iBeg:iFin] += yct
2399        else:
2400            argList = []
2401            for i in range(0,len(refList),chunk):
2402                argList.append(
2403                    (i,refList[i:i+chunk],x,Phase,calcControls,wave,G,g,GB,Vst,
2404                     parmDict,pfx,phfx,hfx,Ka2,lamRatio,kRatio))
2405            for res in mpPool.imap_unordered(ComputeReflectionProfile,argList):
2406                if not res: continue
2407                RefNum,refls,iBeg,iFin,yct = res
2408                yc[iBeg:iFin] += yct
2409                for i,refl in enumerate(refls):
2410                    refList[RefNum+i] = refl
2411            mpPool.close()
2412            mpPool.join()
2413    return yc,yb
2414
2415def ComputeReflectionDerivative(arg):
2416    def cellVaryDerv(pfx,SGData,dpdA): 
2417        if SGData['SGLaue'] in ['-1',]:
2418            return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],
2419                [pfx+'A3',dpdA[3]],[pfx+'A4',dpdA[4]],[pfx+'A5',dpdA[5]]]
2420        elif SGData['SGLaue'] in ['2/m',]:
2421            if SGData['SGUniq'] == 'a':
2422                return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A3',dpdA[3]]]
2423            elif SGData['SGUniq'] == 'b':
2424                return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A4',dpdA[4]]]
2425            else:
2426                return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A5',dpdA[5]]]
2427        elif SGData['SGLaue'] in ['mmm',]:
2428            return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]]]
2429        elif SGData['SGLaue'] in ['4/m','4/mmm']:
2430#            return [[pfx+'A0',dpdA[0]+dpdA[1]],[pfx+'A2',dpdA[2]]]
2431            return [[pfx+'A0',dpdA[0]],[pfx+'A2',dpdA[2]]]
2432        elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
2433#            return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[3]],[pfx+'A2',dpdA[2]]]
2434            return [[pfx+'A0',dpdA[0]],[pfx+'A2',dpdA[2]]]
2435        elif SGData['SGLaue'] in ['3R', '3mR']:
2436            return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[2]],[pfx+'A3',dpdA[3]+dpdA[4]+dpdA[5]]]                       
2437        elif SGData['SGLaue'] in ['m3m','m3']:
2438#            return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[2]]]
2439            return [[pfx+'A0',dpdA[0]]]
2440
2441    (RefNum,refList,x,dependentVars,varylist,dFdvDict,
2442     Phase,calcControls,wave,G,g,GB,A,
2443     parmDict,pfx,phfx,hfx,Ka2,lamRatio,kRatio) = arg
2444    dMdv = np.zeros(shape=(len(varylist),len(x)))
2445    depDerivDict = {}
2446    for j in dependentVars:
2447        depDerivDict[j] = np.zeros(shape=(len(x)))
2448    SGData = Phase['General']['SGData']
2449    shl = max(parmDict[hfx+'SH/L'],0.002)
2450    dx = x[1]-x[0]
2451    iBegO = None
2452    for irefS,refl in enumerate(refList):
2453        iref = irefS + RefNum       
2454        if 'C' in calcControls[hfx+'histType']:        #CW powder
2455            dIdsh,dIdsp,dIdpola,dIdPO,dFdODF,dFdSA = GetIntensityDerv(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
2456            Wd,fmin,fmax = G2pwd.getWidths(refl[5],refl[6],refl[7],shl)
2457            iBeg = np.searchsorted(x,refl[5]-fmin)
2458            iFin = np.searchsorted(x,refl[5]+fmax)
2459            if iBegO is None:
2460                iBegO = iBeg
2461                iFinO = iFin
2462            else:
2463                iBegO = min(iBegO,iBeg)
2464                iFinO = max(iFinO,iFin)
2465            if not iBeg+iFin:       #peak below low limit - skip peak
2466                continue
2467                #return None
2468            elif not iBeg-iFin:     #peak above high limit - done
2469                break
2470                #return None
2471            pos = refl[5]
2472            tanth = tand(pos/2.0)
2473            costh = cosd(pos/2.0)
2474            lenBF = iFin-iBeg
2475            dMdpk = np.zeros(shape=(6,lenBF))
2476            dMdipk = G2pwd.getdFCJVoigt3(refl[5],refl[6],refl[7],shl,x[iBeg:iFin])
2477            for i in range(1,5):
2478                dMdpk[i] += 100.*dx*refl[13]*refl[9]*dMdipk[i]
2479            dMdpk[0] += 100.*dx*refl[13]*refl[9]*dMdipk[0]
2480            dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':np.zeros_like(dMdpk[0])}
2481            if Ka2:
2482                pos2 = refl[5]+lamRatio*tanth       # + 360/pi * Dlam/lam * tan(th)
2483                kdelt = int((pos2-refl[5])/dx)               
2484                iBeg2 = min(lenX,iBeg+kdelt)
2485                iFin2 = min(lenX,iFin+kdelt)
2486                if iBeg2-iFin2:
2487                    iBegO = min(iBegO,iBeg2)
2488                    iFinO = max(iFinO,iFin2)
2489                    lenBF2 = iFin2-iBeg2
2490                    dMdpk2 = np.zeros(shape=(6,lenBF2))
2491                    dMdipk2 = G2pwd.getdFCJVoigt3(pos2,refl[6],refl[7],shl,x[iBeg2:iFin2])
2492                    for i in range(1,5):
2493                        dMdpk2[i] = 100.*dx*refl[13]*refl[9]*kRatio*dMdipk2[i]
2494                    dMdpk2[0] = 100.*dx*refl[13]*refl[9]*kRatio*dMdipk2[0]
2495                    dMdpk2[5] = 100.*dx*refl[13]*dMdipk2[0]
2496                    dervDict2 = {'int':dMdpk2[0],'pos':dMdpk2[1],'sig':dMdpk2[2],'gam':dMdpk2[3],'shl':dMdpk2[4],'L1/L2':dMdpk2[5]*refl[9]}
2497            if 'Pawley' in Phase['General']['Type']:
2498                    try:
2499                        idx = varylist.index(pfx+'PWLref:'+str(pawleyLookup[pfx+'%d,%d,%d'%refl[:3]]))
2500                        dMdv[idx][iBeg:iFin] = dervDict['int']/refl[9]
2501                        if Ka2:
2502                            dMdv[idx][iBeg2:iFin2] = dervDict2['int']/refl[9]
2503                        # Assuming Pawley variables not in constraints
2504                    except ValueError:
2505                        pass
2506            dpdA,dpdw,dpdZ,dpdSh,dpdTr,dpdX,dpdY = GetReflPosDerv(refl,wave,A,hfx,calcControls,parmDict)
2507            names = {hfx+'Scale':[dIdsh,'int'],hfx+'Polariz.':[dIdpola,'int'],phfx+'Scale':[dIdsp,'int'],
2508                    hfx+'U':[tanth**2,'sig'],hfx+'V':[tanth,'sig'],hfx+'W':[1.0,'sig'],
2509                    hfx+'X':[1.0/costh,'gam'],hfx+'Y':[tanth,'gam'],hfx+'SH/L':[1.0,'shl'],
2510                    hfx+'I(L2)/I(L1)':[1.0,'L1/L2'],hfx+'Zero':[dpdZ,'pos'],hfx+'Lam':[dpdw,'pos'],
2511                    hfx+'Shift':[dpdSh,'pos'],hfx+'Transparency':[dpdTr,'pos'],hfx+'DisplaceX':[dpdX,'pos'],
2512                    hfx+'DisplaceY':[dpdY,'pos'],}
2513            for name in names:
2514                item = names[name]
2515                if name in varylist:
2516                    dMdv[varylist.index(name)][iBeg:iFin] += item[0]*dervDict[item[1]]
2517                    if Ka2:
2518                        dMdv[varylist.index(name)][iBeg2:iFin2] += item[0]*dervDict2[item[1]]
2519                elif name in dependentVars:
2520                    if Ka2:
2521                        depDerivDict[name][iBeg2:iFin2] += item[0]*dervDict2[item[1]]
2522                    depDerivDict[name][iBeg:iFin] += item[0]*dervDict[item[1]]
2523            for iPO in dIdPO:
2524                if iPO in varylist:
2525                    dMdv[varylist.index(iPO)][iBeg:iFin] += dIdPO[iPO]*dervDict['int']
2526                    if Ka2:
2527                        dMdv[varylist.index(iPO)][iBeg2:iFin2] += dIdPO[iPO]*dervDict2['int']
2528                elif iPO in dependentVars:
2529                    depDerivDict[iPO][iBeg:iFin] += dIdPO[iPO]*dervDict['int']
2530                    if Ka2:
2531                        depDerivDict[iPO][iBeg2:iFin2] += dIdPO[iPO]*dervDict2['int']
2532            for i,name in enumerate(['omega','chi','phi']):
2533                aname = pfx+'SH '+name
2534                if aname in varylist:
2535                    dMdv[varylist.index(aname)][iBeg:iFin] += dFdSA[i]*dervDict['int']
2536                    if Ka2:
2537                        dMdv[varylist.index(aname)][iBeg2:iFin2] += dFdSA[i]*dervDict2['int']
2538                elif aname in dependentVars:
2539                    depDerivDict[aname][iBeg:iFin] += dFdSA[i]*dervDict['int']
2540                    if Ka2:
2541                        depDerivDict[aname][iBeg2:iFin2] += dFdSA[i]*dervDict2['int']
2542            for iSH in dFdODF:
2543                if iSH in varylist:
2544                    dMdv[varylist.index(iSH)][iBeg:iFin] += dFdODF[iSH]*dervDict['int']
2545                    if Ka2:
2546                        dMdv[varylist.index(iSH)][iBeg2:iFin2] += dFdODF[iSH]*dervDict2['int']
2547                elif iSH in dependentVars:
2548                    depDerivDict[iSH][iBeg:iFin] += dFdODF[iSH]*dervDict['int']
2549                    if Ka2:
2550                        depDerivDict[iSH][iBeg2:iFin2] += dFdODF[iSH]*dervDict2['int']
2551            cellDervNames = cellVaryDerv(pfx,SGData,dpdA)
2552            for name,dpdA in cellDervNames:
2553                if name in varylist:
2554                    dMdv[varylist.index(name)][iBeg:iFin] += dpdA*dervDict['pos']
2555                    if Ka2:
2556                        dMdv[varylist.index(name)][iBeg2:iFin2] += dpdA*dervDict2['pos']
2557                elif name in dependentVars:
2558                    depDerivDict[name][iBeg:iFin] += dpdA*dervDict['pos']
2559                    if Ka2:
2560                        depDerivDict[name][iBeg2:iFin2] += dpdA*dervDict2['pos']
2561            dDijDict = GetHStrainShiftDerv(refl,SGData,phfx)
2562            for name in dDijDict:
2563                if name in varylist:
2564                    dMdv[varylist.index(name)][iBeg:iFin] += dDijDict[name]*dervDict['pos']
2565                    if Ka2:
2566                        dMdv[varylist.index(name)][iBeg2:iFin2] += dDijDict[name]*dervDict2['pos']
2567                elif name in dependentVars:
2568                    depDerivDict[name][iBeg:iFin] += dDijDict[name]*dervDict['pos']
2569                    if Ka2:
2570                        depDerivDict[name][iBeg2:iFin2] += dDijDict[name]*dervDict2['pos']
2571            gamDict = GetSampleGamDerv(refl,wave,G,GB,phfx,calcControls,parmDict)
2572            for name in gamDict:
2573                if name in varylist:
2574                    dMdv[varylist.index(name)][iBeg:iFin] += gamDict[name]*dervDict['gam']
2575                    if Ka2:
2576                        dMdv[varylist.index(name)][iBeg2:iFin2] += gamDict[name]*dervDict2['gam']
2577                elif name in dependentVars:
2578                    depDerivDict[name][iBeg:iFin] += gamDict[name]*dervDict['gam']
2579                    if Ka2:
2580                        depDerivDict[name][iBeg2:iFin2] += gamDict[name]*dervDict2['gam']
2581                                               
2582        elif 'T' in calcControls[hfx+'histType']:
2583            print 'TOF Undefined at present'
2584            raise Exception    #no TOF yet
2585        #do atom derivatives -  for F,X & U so far             
2586        corr = dervDict['int']/refl[9]
2587        if Ka2:
2588            corr2 = dervDict2['int']/refl[9]
2589        for name in varylist+dependentVars:
2590            try:
2591                aname = name.split(pfx)[1][:2]
2592                if aname not in ['Af','dA','AU']: continue # skip anything not an atom param
2593            except IndexError:
2594                continue
2595            if name in varylist:
2596                dMdv[varylist.index(name)][iBeg:iFin] += dFdvDict[name][iref]*corr
2597                if Ka2:
2598                    dMdv[varylist.index(name)][iBeg2:iFin2] += dFdvDict[name][iref]*corr2
2599            elif name in dependentVars:
2600                depDerivDict[name][iBeg:iFin] += dFdvDict[name][iref]*corr
2601                if Ka2:
2602                    depDerivDict[name][iBeg2:iFin2] += dFdvDict[name][iref]*corr2
2603        # end loop over reflections
2604    return depDerivDict,dMdv[:,iBegO:iFinO],iBegO,iFinO
2605       
2606def getPowderProfileDerv(parmDict,x,varylist,Histogram,Phases,calcControls,pawleyLookup):   
2607    if calcControls['max Hprocess'] <= 1 and calcControls['max Rprocess'] > 1:
2608        mpPool = mp.Pool()
2609        #print 'mp derivs'
2610    #else:
2611        #print 'single thread derivs'
2612    # create a list of dependent variables and set up a dictionary to hold their derivatives
2613    dependentVars = G2mv.GetDependentVars()
2614    depDerivDict = {}
2615    for j in dependentVars:
2616        depDerivDict[j] = np.zeros(shape=(len(x)))
2617    #print 'dependent vars',dependentVars
2618    lenX = len(x)               
2619    hId = Histogram['hId']
2620    hfx = ':%d:'%(hId)
2621    bakType = calcControls[hfx+'bakType']
2622    dMdv = np.zeros(shape=(len(varylist),len(x)))
2623    dMdb,dMddb = G2pwd.getBackgroundDerv(hfx,parmDict,bakType,x)
2624    if hfx+'Back:0' in varylist: # for now assume that Back:x vars to not appear in constraints
2625        bBpos =varylist.index(hfx+'Back:0')
2626        dMdv[bBpos:bBpos+len(dMdb)] = dMdb
2627    names = [hfx+'DebyeA',hfx+'DebyeR',hfx+'DebyeU']
2628    for name in varylist:
2629        if 'Debye' in name:
2630            id = int(name.split(':')[-1])
2631            parm = name[:int(name.rindex(':'))]
2632            ip = names.index(parm)
2633            dMdv[varylist.index(name)] = dMddb[3*id+ip]
2634    if 'C' in calcControls[hfx+'histType']:   
2635        dx = x[1]-x[0]
2636        shl = max(parmDict[hfx+'SH/L'],0.002)
2637        if hfx+'Lam1' in parmDict.keys():
2638            wave = parmDict[hfx+'Lam1']
2639            Ka2 = True
2640            lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
2641            kRatio = parmDict[hfx+'I(L2)/I(L1)']
2642        else:
2643            wave = parmDict[hfx+'Lam']
2644            Ka2 = False
2645            lamRatio = None
2646            kRatio = None
2647    else:
2648        print 'TOF Undefined at present'
2649        raise ValueError
2650    for phase in Histogram['Reflection Lists']:
2651        refList = Histogram['Reflection Lists'][phase]
2652        Phase = Phases[phase]
2653        SGData = Phase['General']['SGData']
2654        pId = Phase['pId']
2655        pfx = '%d::'%(pId)
2656        phfx = '%d:%d:'%(pId,hId)
2657        A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
2658        G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
2659        GA,GB = G2lat.Gmat2AB(G)    #Orthogonalization matricies
2660        if 'Pawley' not in Phase['General']['Type']:
2661            dFdvDict = StructureFactorDerv(refList,G,hfx,pfx,SGData,calcControls,parmDict)
2662        else:
2663            for iref,refl in enumerate(refList):
2664                try:
2665                    refl[9] = abs(parmDict[pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%refl[:3]])])
2666                except KeyError:
2667                    #print ' ***Error %d,%d,%d missing from Pawley reflection list ***'%refl[:3]
2668                    continue
2669        chunk = 200
2670        if (calcControls['max Hprocess'] > 1 or
2671            calcControls['max Rprocess'] == 1 or
2672            mp.cpu_count() < 2 or len(refList) <  2 * chunk):
2673            TdepDerivDict, TdMdv, iBegO, iFinO = ComputeReflectionDerivative(
2674                (0,refList,x,dependentVars,varylist,dFdvDict,
2675                 Phase,calcControls,wave,G,g,GB,A,
2676                 parmDict,pfx,phfx,hfx,Ka2,lamRatio,kRatio))
2677            dMdv[:,iBegO:iFinO] += TdMdv
2678            for j in dependentVars:         
2679                depDerivDict[j] += TdepDerivDict[j]
2680        else:
2681            argList = []
2682            for iref in range(0,len(refList),chunk):
2683                argList.append(
2684                    (iref,refList[iref:iref+chunk],x,dependentVars,varylist,dFdvDict,
2685                     Phase,calcControls,wave,G,g,GB,A,
2686                     parmDict,pfx,phfx,hfx,Ka2,lamRatio,kRatio))
2687            for res in mpPool.imap_unordered(ComputeReflectionDerivative,argList):
2688                if res is None: continue
2689                TdepDerivDict, TdMdv, iBegO, iFinO = res
2690                dMdv[:,iBegO:iFinO] += TdMdv
2691                for j in dependentVars:
2692                    depDerivDict[j] += TdepDerivDict[j]
2693
2694    # now process derivatives in constraints
2695    G2mv.Dict2Deriv(varylist,depDerivDict,dMdv)
2696    return dMdv
2697
2698def dervRefine(values,HistoPhases,parmdict,varylist,calcControls,pawleyLookup,dlg):
2699    parmdict.update(zip(varylist,values))
2700    G2mv.Dict2Map(parmdict,varylist)
2701    Histograms,Phases = HistoPhases
2702    nvar = len(varylist)
2703    dMdv = np.empty(0)
2704    histoList = Histograms.keys()
2705    histoList.sort()
2706    for histogram in histoList:
2707        if 'PWDR' in histogram[:4]:
2708            Histogram = Histograms[histogram]
2709            hId = Histogram['hId']
2710            hfx = ':%d:'%(hId)
2711            Limits = calcControls[hfx+'Limits']
2712            x,y,w,yc,yb,yd = Histogram['Data']
2713            xB = np.searchsorted(x,Limits[0])
2714            xF = np.searchsorted(x,Limits[1])
2715            dMdvh = np.sqrt(w[xB:xF])*getPowderProfileDerv(parmdict,x[xB:xF],
2716                varylist,Histogram,Phases,calcControls,pawleyLookup)
2717            if len(dMdv):
2718                dMdv = np.concatenate((dMdv.T,dMdvh.T)).T
2719            else:
2720                dMdv = dMdvh
2721    return dMdv
2722
2723def ComputePowderHessian(args):
2724    G2lat.InitTexture()     #initialize fortran arrays for spherical harmonics
2725    Histogram,parmdict,varylist,Phases,calcControls,pawleyLookup = args
2726    hId = Histogram['hId']
2727    hfx = ':%d:'%(hId)
2728    Limits = calcControls[hfx+'Limits']
2729    x,y,w,yc,yb,yd = Histogram['Data']
2730    dy = y-yc
2731    xB = np.searchsorted(x,Limits[0])
2732    xF = np.searchsorted(x,Limits[1])
2733    dMdvh = np.sqrt(w[xB:xF])*getPowderProfileDerv(
2734        parmdict,x[xB:xF],
2735        varylist,Histogram,Phases,calcControls,pawleyLookup)
2736    Vec = np.sum(dMdvh*np.sqrt(w[xB:xF])*dy[xB:xF],axis=1)
2737    Hess = np.inner(dMdvh,dMdvh)
2738    return Vec,Hess
2739
2740#def HessRefine(values,HistoPhases,parmdict,varylist,calcControls,pawleyLookup,dlg):
2741#    parmdict.update(zip(varylist,values))
2742#    G2mv.Dict2Map(parmdict,varylist)
2743#    Histograms,Phases = HistoPhases
2744#    nvar = len(varylist)
2745#    Hess = np.empty(0)
2746#    histoList = Histograms.keys()
2747#    histoList.sort()
2748#    for histogram in histoList:
2749#        if 'PWDR' in histogram[:4]:
2750#            Histogram = Histograms[histogram]
2751#            hId = Histogram['hId']
2752#            hfx = ':%d:'%(hId)
2753#            Limits = calcControls[hfx+'Limits']
2754#            x,y,w,yc,yb,yd = Histogram['Data']
2755#            dy = y-yc
2756#            xB = np.searchsorted(x,Limits[0])
2757#            xF = np.searchsorted(x,Limits[1])
2758#            dMdvh = np.sqrt(w[xB:xF])*getPowderProfileDerv(parmdict,x[xB:xF],
2759#                varylist,Histogram,Phases,calcControls,pawleyLookup)
2760#            if dlg:
2761#                dlg.Update(Histogram['wRp'],newmsg='Hessian for histogram %d Rwp=%8.3f%s'%(hId,Histogram['wRp'],'%'))[0]
2762#            if len(Hess):
2763#                Vec += np.sum(dMdvh*np.sqrt(w[xB:xF])*dy[xB:xF],axis=1)
2764#                Hess += np.inner(dMdvh,dMdvh)
2765#            else:
2766#                Vec = np.sum(dMdvh*np.sqrt(w[xB:xF])*dy[xB:xF],axis=1)
2767#                Hess = np.inner(dMdvh,dMdvh)
2768#    return Vec,Hess
2769
2770def HessRefine(values,HistoPhases,parmdict,varylist,calcControls,pawleyLookup,dlg):
2771    parmdict.update(zip(varylist,values))
2772    G2mv.Dict2Map(parmdict,varylist)
2773    Histograms,Phases = HistoPhases
2774    nvar = len(varylist)
2775    HessSum = np.zeros((nvar,nvar))
2776    VecSum = np.zeros(nvar)
2777    histoList = Histograms.keys()
2778    histoList.sort()
2779    argList = []
2780    MaxProcess = calcControls['max Hprocess']
2781    for histogram in histoList:
2782        if 'PWDR' in histogram[:4]:
2783            Histogram = Histograms[histogram]
2784            argList.append([
2785                Histogram,parmdict,varylist,Phases,
2786                calcControls,pawleyLookup])
2787    if MaxProcess > 1:
2788        mpPool = mp.Pool()
2789        for Vec,Hess in mpPool.imap_unordered(ComputePowderHessian,argList):
2790            VecSum += Vec
2791            HessSum += Hess
2792        mpPool.close()
2793        mpPool.join()
2794    else:
2795        for arg in argList:
2796            Vec,Hess = ComputePowderHessian(arg)
2797            VecSum += Vec
2798            HessSum += Hess
2799    return VecSum,HessSum
2800
2801def ComputePowderProfile(args):
2802    G2lat.InitTexture()     #initialize fortran arrays for spherical harmonics
2803    Histogram,parmdict,varylist,Phases,calcControls,pawleyLookup,histkey = args
2804    hId = Histogram['hId']
2805    hfx = ':%d:'%(hId)
2806    x,y,w,yc,yb,yd = Histogram['Data']
2807    Limits = calcControls[hfx+'Limits']
2808    xB = np.searchsorted(x,Limits[0])
2809    xF = np.searchsorted(x,Limits[1])
2810    yc,yb = getPowderProfile(parmdict,x[xB:xF],
2811                            varylist,Histogram,Phases,calcControls,
2812                            pawleyLookup)
2813    return xB,xF,yc,yb,Histogram['Reflection Lists'],histkey
2814
2815#def errRefine(values,HistoPhases,parmdict,varylist,calcControls,pawleyLookup,dlg):       
2816#    parmdict.update(zip(varylist,values))
2817#    Values2Dict(parmdict, varylist, values)
2818#    G2mv.Dict2Map(parmdict,varylist)
2819#    Histograms,Phases = HistoPhases
2820#    M = np.empty(0)
2821#    sumwYo = 0
2822#    Nobs = 0
2823#    histoList = Histograms.keys()
2824#    histoList.sort()
2825#    for histogram in histoList:
2826#        if 'PWDR' in histogram[:4]:
2827#            Histogram = Histograms[histogram]
2828#            hId = Histogram['hId']
2829#            hfx = ':%d:'%(hId)
2830#            Limits = calcControls[hfx+'Limits']
2831#            x,y,w,yc,yb,yd = Histogram['Data']
2832#            yc *= 0.0                           #zero full calcd profiles
2833#            yb *= 0.0
2834#            yd *= 0.0
2835#            xB = np.searchsorted(x,Limits[0])
2836#            xF = np.searchsorted(x,Limits[1])
2837#            Histogram['Nobs'] = xF-xB
2838#            Nobs += Histogram['Nobs']
2839#            Histogram['sumwYo'] = np.sum(w[xB:xF]*y[xB:xF]**2)
2840#            sumwYo += Histogram['sumwYo']
2841#            yc[xB:xF],yb[xB:xF] = getPowderProfile(parmdict,x[xB:xF],
2842#                varylist,Histogram,Phases,calcControls,pawleyLookup)
2843#            yc[xB:xF] += yb[xB:xF]
2844#            yd[xB:xF] = y[xB:xF]-yc[xB:xF]
2845#            Histogram['sumwYd'] = np.sum(np.sqrt(w[xB:xF])*(yd[xB:xF]))
2846#            wdy = -np.sqrt(w[xB:xF])*(yd[xB:xF])
2847#            Histogram['wRp'] = min(100.,np.sqrt(np.sum(wdy**2)/Histogram['sumwYo'])*100.)
2848#            if dlg:
2849#                dlg.Update(Histogram['wRp'],newmsg='For histogram %d Rwp=%8.3f%s'%(hId,Histogram['wRp'],'%'))[0]
2850#            M = np.concatenate((M,wdy))
2851#    Histograms['sumwYo'] = sumwYo
2852#    Histograms['Nobs'] = Nobs
2853#    Rwp = min(100.,np.sqrt(np.sum(M**2)/sumwYo)*100.)
2854#    if dlg:
2855#        GoOn = dlg.Update(Rwp,newmsg='%s%8.3f%s'%('Powder profile Rwp =',Rwp,'%'))[0]
2856#        if not GoOn:
2857#            parmDict['saved values'] = values
2858#            raise Exception         #Abort!!
2859#    return M
2860#   
2861def errRefine(values,HistoPhases,parmdict,varylist,calcControls,pawleyLookup,dlg):       
2862    parmdict.update(zip(varylist,values))
2863    Values2Dict(parmdict, varylist, values)
2864    G2mv.Dict2Map(parmdict,varylist)
2865    Histograms,Phases = HistoPhases
2866    M = np.empty(0)
2867    sumwYo = 0
2868    Nobs = 0
2869    histoList = Histograms.keys()
2870    histoList.sort()
2871    argList = []
2872    MaxProcess = calcControls['max Hprocess']
2873    for histogram in histoList:
2874        if 'PWDR' in histogram[:4]:
2875            Histogram = Histograms[histogram]
2876            argList.append(
2877                [Histogram,parmdict,varylist,Phases,calcControls,pawleyLookup,histogram]
2878                )
2879    if MaxProcess > 1: 
2880        mpPool = mp.Pool()
2881        for (xB,xF,ycSect,ybSect,RL,histkey
2882             ) in mpPool.imap_unordered(ComputePowderProfile,argList):
2883            Histogram = Histograms[histkey]
2884            Histogram['Reflection Lists'] = RL
2885            x,y,w,yc,yb,yd = Histogram['Data']
2886            yc *= 0.0                           #zero full calcd profiles
2887            yb *= 0.0
2888            yd *= 0.0
2889            Histogram['Nobs'] = xF-xB
2890            Nobs += Histogram['Nobs']
2891            Histogram['sumwYo'] = np.sum(w[xB:xF]*y[xB:xF]**2)
2892            sumwYo += Histogram['sumwYo']
2893            yc[xB:xF] = ycSect
2894            yb[xB:xF] = ybSect
2895            yc[xB:xF] += yb[xB:xF]
2896            yd[xB:xF] = y[xB:xF]-yc[xB:xF]
2897            Histogram['sumwYd'] = np.sum(np.sqrt(w[xB:xF])*(yd[xB:xF]))
2898            wdy = -np.sqrt(w[xB:xF])*(yd[xB:xF])
2899            Histogram['wRp'] = min(100.,np.sqrt(np.sum(wdy**2)/Histogram['sumwYo'])*100.)
2900            M = np.concatenate((M,wdy))
2901        mpPool.close()
2902        mpPool.join()
2903    else:
2904        for arg in argList:
2905            xB,xF,ycSect,ybSect,RL,histkey = ComputePowderProfile(arg)
2906            Histogram = arg[0]
2907            hId = Histogram['hId']
2908            x,y,w,yc,yb,yd = Histogram['Data']
2909            yc *= 0.0                           #zero full calcd profiles
2910            yb *= 0.0
2911            yd *= 0.0
2912            Histogram['Nobs'] = xF-xB
2913            Nobs += Histogram['Nobs']
2914            Histogram['sumwYo'] = np.sum(w[xB:xF]*y[xB:xF]**2)
2915            sumwYo += Histogram['sumwYo']
2916           
2917            yc[xB:xF] = ycSect
2918            yb[xB:xF] = ybSect
2919            yc[xB:xF] += yb[xB:xF]
2920            yd[xB:xF] = y[xB:xF]-yc[xB:xF]
2921            Histogram['sumwYd'] = np.sum(np.sqrt(w[xB:xF])*(yd[xB:xF]))
2922            wdy = -np.sqrt(w[xB:xF])*(yd[xB:xF])
2923            Histogram['wRp'] = min(100.,np.sqrt(np.sum(wdy**2)/Histogram['sumwYo'])*100.)
2924            if dlg:
2925                dlg.Update(Histogram['wRp'],newmsg='For histogram %d Rwp=%8.3f%s'%(hId,Histogram['wRp'],'%'))[0]
2926            M = np.concatenate((M,wdy))
2927
2928    Histograms['sumwYo'] = sumwYo
2929    Histograms['Nobs'] = Nobs
2930    Rwp = min(100.,np.sqrt(np.sum(M**2)/sumwYo)*100.)
2931    if dlg:
2932        GoOn = dlg.Update(Rwp,newmsg='%s%8.3f%s'%('Powder profile Rwp =',Rwp,'%'))[0]
2933        if not GoOn:
2934            parmDict['saved values'] = values
2935            raise Exception         #Abort!!
2936    return M       
2937                   
2938def Refine(GPXfile,dlg):
2939    G2lat.InitTexture()     #initialize fortran arrays for spherical harmonics
2940    ShowBanner()
2941    varyList = []
2942    parmDict = {}
2943    G2mv.InitVars()   
2944    Controls = GetControls(GPXfile)
2945    ShowControls(Controls)
2946    calcControls = {}
2947    calcControls.update(Controls)           
2948    constrDict,constrFlag,fixedList = GetConstraints(GPXfile)
2949    Histograms,Phases = GetUsedHistogramsAndPhases(GPXfile)
2950    if not Phases:
2951        print ' *** ERROR - you have no histograms to refine! ***'
2952        print ' *** Refine aborted ***'
2953        raise Exception
2954    if not Histograms:
2955        print ' *** ERROR - you have no data to refine with! ***'
2956        print ' *** Refine aborted ***'
2957        raise Exception       
2958    Natoms,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables = GetPhaseData(Phases)
2959    calcControls['Natoms'] = Natoms
2960    calcControls['FFtables'] = FFtables
2961    calcControls['BLtables'] = BLtables
2962    hapVary,hapDict,controlDict = GetHistogramPhaseData(Phases,Histograms)
2963    calcControls.update(controlDict)
2964    histVary,histDict,controlDict = GetHistogramData(Histograms)
2965    calcControls.update(controlDict)
2966    varyList = phaseVary+hapVary+histVary
2967    parmDict.update(phaseDict)
2968    parmDict.update(hapDict)
2969    parmDict.update(histDict)
2970    GetFprime(calcControls,Histograms)
2971    # do constraint processing
2972    try:
2973        groups,parmlist = G2mv.GroupConstraints(constrDict)
2974        G2mv.GenerateConstraints(groups,parmlist,varyList,constrDict,constrFlag,fixedList)
2975    except:
2976        print ' *** ERROR - your constraints are internally inconsistent ***'
2977        raise Exception(' *** Refine aborted ***')
2978    # check to see which generated parameters are fully varied
2979    msg = G2mv.SetVaryFlags(varyList)
2980    if msg:
2981        print ' *** ERROR - you have not set the refine flags for constraints consistently! ***'
2982        print msg
2983        raise Exception(' *** Refine aborted ***')
2984    G2mv.Map2Dict(parmDict,varyList)
2985#    print G2mv.VarRemapShow(varyList)
2986
2987    while True:
2988        begin = time.time()
2989        values =  np.array(Dict2Values(parmDict, varyList))
2990        Ftol = Controls['min dM/M']
2991        Factor = Controls['shift factor']
2992        maxCyc = Controls['max cyc']
2993        if 'Jacobian' in Controls['deriv type']:           
2994            result = so.leastsq(errRefine,values,Dfun=dervRefine,full_output=True,
2995                ftol=Ftol,col_deriv=True,factor=Factor,
2996                args=([Histograms,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg))
2997            ncyc = int(result[2]['nfev']/2)
2998        elif 'Hessian' in Controls['deriv type']:
2999            result = G2mth.HessianLSQ(errRefine,values,Hess=HessRefine,ftol=Ftol,maxcyc=maxCyc,
3000                args=([Histograms,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg))
3001            ncyc = result[2]['num cyc']+1                           
3002        else:           #'numeric'
3003            result = so.leastsq(errRefine,values,full_output=True,ftol=Ftol,epsfcn=1.e-8,factor=Factor,
3004                args=([Histograms,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg))
3005            ncyc = int(result[2]['nfev']/len(varyList))
3006#        table = dict(zip(varyList,zip(values,result[0],(result[0]-values))))
3007#        for item in table: print item,table[item]               #useful debug - are things shifting?
3008        runtime = time.time()-begin
3009        chisq = np.sum(result[2]['fvec']**2)
3010        Values2Dict(parmDict, varyList, result[0])
3011        G2mv.Dict2Map(parmDict,varyList)
3012       
3013        Rwp = np.sqrt(chisq/Histograms['sumwYo'])*100.      #to %
3014        GOF = chisq/(Histograms['Nobs']-len(varyList))
3015        print '\n Refinement results:'
3016        print 135*'-'
3017        print ' Number of function calls:',result[2]['nfev'],' Number of observations: ',Histograms['Nobs'],' Number of parameters: ',len(varyList)
3018        print ' Refinement time = %8.3fs, %8.3fs/cycle, for %d cycles'%(runtime,runtime/ncyc,ncyc)
3019        print ' wRp = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f'%(Rwp,chisq,GOF)
3020        try:
3021            covMatrix = result[1]*GOF
3022            sig = np.sqrt(np.diag(covMatrix))
3023            if np.any(np.isnan(sig)):
3024                print '*** Least squares aborted - some invalid esds possible ***'
3025#            table = dict(zip(varyList,zip(values,result[0],(result[0]-values)/sig)))
3026#            for item in table: print item,table[item]               #useful debug - are things shifting?
3027            break                   #refinement succeeded - finish up!
3028        except TypeError:          #result[1] is None on singular matrix
3029            print '**** Refinement failed - singular matrix ****'
3030            if 'Hessian' in Controls['deriv type']:
3031                for i in result[2]['psing'].reverse():
3032                        print 'Removing parameter: ',varyList[i]
3033                        del(varyList[i])                   
3034            else:
3035                Ipvt = result[2]['ipvt']
3036                for i,ipvt in enumerate(Ipvt):
3037                    if not np.sum(result[2]['fjac'],axis=1)[i]:
3038                        print 'Removing parameter: ',varyList[ipvt-1]
3039                        del(varyList[ipvt-1])
3040                        break
3041
3042#    print 'dependentParmList: ',G2mv.dependentParmList
3043#    print 'arrayList: ',G2mv.arrayList
3044#    print 'invarrayList: ',G2mv.invarrayList
3045#    print 'indParmList: ',G2mv.indParmList
3046#    print 'fixedDict: ',G2mv.fixedDict
3047#    print 'test1'
3048    GetFobsSq(Histograms,Phases,parmDict,calcControls)
3049#    print 'test2'
3050    sigDict = dict(zip(varyList,sig))
3051    newCellDict = GetNewCellParms(parmDict,varyList)
3052    newAtomDict = ApplyXYZshifts(parmDict,varyList)
3053    covData = {'variables':result[0],'varyList':varyList,'sig':sig,
3054        'covMatrix':covMatrix,'title':GPXfile,'newAtomDict':newAtomDict,'newCellDict':newCellDict}
3055    # add the uncertainties into the esd dictionary (sigDict)
3056    sigDict.update(G2mv.ComputeDepESD(covMatrix,varyList,parmDict))
3057    SetPhaseData(parmDict,sigDict,Phases,covData)
3058    SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms)
3059    SetHistogramData(parmDict,sigDict,Histograms)
3060    G2mv.PrintIndependentVars(parmDict,varyList,sigDict)
3061    SetUsedHistogramsAndPhases(GPXfile,Histograms,Phases,covData)
3062   
3063#for testing purposes!!!
3064#    import cPickle
3065#    file = open('structTestdata.dat','wb')
3066#    cPickle.dump(parmDict,file,1)
3067#    cPickle.dump(varyList,file,1)
3068#    for histogram in Histograms:
3069#        if 'PWDR' in histogram[:4]:
3070#            Histogram = Histograms[histogram]
3071#    cPickle.dump(Histogram,file,1)
3072#    cPickle.dump(Phases,file,1)
3073#    cPickle.dump(calcControls,file,1)
3074#    cPickle.dump(pawleyLookup,file,1)
3075#    file.close()
3076
3077    ShowControls(Controls)
3078    print 'Last Refinement time = %8.3fs'%(runtime)
3079    if dlg:
3080        return Rwp
3081
3082def SeqRefine(GPXfile,dlg):
3083    G2lat.InitTexture()     #initialize fortran arrays for spherical harmonics
3084    ShowBanner()
3085    print ' Sequential Refinement'
3086    G2mv.InitVars()   
3087    Controls = GetControls(GPXfile)
3088    ShowControls(Controls)           
3089    constrDict,constrFlag,fixedList = GetConstraints(GPXfile)
3090    Histograms,Phases = GetUsedHistogramsAndPhases(GPXfile)
3091    if not Phases:
3092        print ' *** ERROR - you have no histograms to refine! ***'
3093        print ' *** Refine aborted ***'
3094        raise Exception
3095    if not Histograms:
3096        print ' *** ERROR - you have no data to refine with! ***'
3097        print ' *** Refine aborted ***'
3098        raise Exception
3099    Natoms,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables = GetPhaseData(Phases,False)
3100    if 'Seq Data' in Controls:
3101        histNames = Controls['Seq Data']
3102    else:
3103        histNames = GetHistogramNames(GPXfile,['PWDR',])
3104    if 'Reverse Seq' in Controls:
3105        if Controls['Reverse Seq']:
3106            histNames.reverse()
3107    SeqResult = {'histNames':histNames}
3108    makeBack = True
3109    for ihst,histogram in enumerate(histNames):
3110        ifPrint = False
3111        if dlg:
3112            dlg.SetTitle('Residual for histogram '+str(ihst))
3113        calcControls = {}
3114        calcControls['Natoms'] = Natoms
3115        calcControls['FFtables'] = FFtables
3116        calcControls['BLtables'] = BLtables
3117        varyList = []
3118        parmDict = {}
3119        Histo = {histogram:Histograms[histogram],}
3120        hapVary,hapDict,controlDict = GetHistogramPhaseData(Phases,Histo,False)
3121        calcControls.update(controlDict)
3122        histVary,histDict,controlDict = GetHistogramData(Histo,False)
3123        calcControls.update(controlDict)
3124        varyList = phaseVary+hapVary+histVary
3125        if not ihst:
3126            saveVaryList = varyList[:]
3127            for i,item in enumerate(saveVaryList):
3128                items = item.split(':')
3129                if items[1]:
3130                    items[1] = ''
3131                item = ':'.join(items)
3132                saveVaryList[i] = item
3133            SeqResult['varyList'] = saveVaryList
3134        else:
3135            newVaryList = varyList[:]
3136            for i,item in enumerate(newVaryList):
3137                items = item.split(':')
3138                if items[1]:
3139                    items[1] = ''
3140                item = ':'.join(items)
3141                newVaryList[i] = item
3142            if newVaryList != SeqResult['varyList']:
3143                print newVaryList
3144                print SeqResult['varyList']
3145                print '**** ERROR - variable list for this histogram does not match previous'
3146                raise Exception
3147        parmDict.update(phaseDict)
3148        parmDict.update(hapDict)
3149        parmDict.update(histDict)
3150        GetFprime(calcControls,Histo)
3151        constrDict,constrFlag,fixedList = G2mv.InputParse([])        #constraints go here?
3152        groups,parmlist = G2mv.GroupConstraints(constrDict)
3153        G2mv.GenerateConstraints(groups,parmlist,varyList,constrDict,constrFlag,fixedList)
3154        G2mv.Map2Dict(parmDict,varyList)
3155   
3156        while True:
3157            begin = time.time()
3158            values =  np.array(Dict2Values(parmDict, varyList))
3159            Ftol = Controls['min dM/M']
3160            Factor = Controls['shift factor']
3161            maxCyc = Controls['max cyc']
3162
3163            if 'Jacobian' in Controls['deriv type']:           
3164                result = so.leastsq(errRefine,values,Dfun=dervRefine,full_output=True,
3165                    ftol=Ftol,col_deriv=True,factor=Factor,
3166                    args=([Histo,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg))
3167                ncyc = int(result[2]['nfev']/2)
3168            elif 'Hessian' in Controls['deriv type']:
3169                result = G2mth.HessianLSQ(errRefine,values,Hess=HessRefine,ftol=Ftol,maxcyc=maxCyc,
3170                    args=([Histo,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg))
3171                ncyc = result[2]['num cyc']+1                           
3172            else:           #'numeric'
3173                result = so.leastsq(errRefine,values,full_output=True,ftol=Ftol,epsfcn=1.e-8,factor=Factor,
3174                    args=([Histo,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg))
3175                ncyc = int(result[2]['nfev']/len(varyList))
3176
3177
3178
3179            runtime = time.time()-begin
3180            chisq = np.sum(result[2]['fvec']**2)
3181            Values2Dict(parmDict, varyList, result[0])
3182            G2mv.Dict2Map(parmDict,varyList)
3183           
3184            Rwp = np.sqrt(chisq/Histo['sumwYo'])*100.      #to %
3185            GOF = chisq/(Histo['Nobs']-len(varyList))
3186            print '\n Refinement results for histogram: v'+histogram
3187            print 135*'-'
3188            print ' Number of function calls:',result[2]['nfev'],' Number of observations: ',Histo['Nobs'],' Number of parameters: ',len(varyList)
3189            print ' Refinement time = %8.3fs, %8.3fs/cycle, for %d cycles'%(runtime,runtime/ncyc,ncyc)
3190            print ' wRp = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f'%(Rwp,chisq,GOF)
3191            try:
3192                covMatrix = result[1]*GOF
3193                sig = np.sqrt(np.diag(covMatrix))
3194                if np.any(np.isnan(sig)):
3195                    print '*** Least squares aborted - some invalid esds possible ***'
3196                    ifPrint = True
3197                break                   #refinement succeeded - finish up!
3198            except TypeError:          #result[1] is None on singular matrix
3199                print '**** Refinement failed - singular matrix ****'
3200                if 'Hessian' in Controls['deriv type']:
3201                    for i in result[2]['psing'].reverse():
3202                            print 'Removing parameter: ',varyList[i]
3203                            del(varyList[i])                   
3204                else:
3205                    Ipvt = result[2]['ipvt']
3206                    for i,ipvt in enumerate(Ipvt):
3207                        if not np.sum(result[2]['fjac'],axis=1)[i]:
3208                            print 'Removing parameter: ',varyList[ipvt-1]
3209                            del(varyList[ipvt-1])
3210                            break
3211   
3212        GetFobsSq(Histo,Phases,parmDict,calcControls)
3213        sigDict = dict(zip(varyList,sig))
3214        newCellDict = GetNewCellParms(parmDict,varyList)
3215        newAtomDict = ApplyXYZshifts(parmDict,varyList)
3216        covData = {'variables':result[0],'varyList':varyList,'sig':sig,
3217            'covMatrix':covMatrix,'title':histogram,'newAtomDict':newAtomDict,'newCellDict':newCellDict}
3218        SetHistogramPhaseData(parmDict,sigDict,Phases,Histo,ifPrint)
3219        SetHistogramData(parmDict,sigDict,Histo,ifPrint)
3220        SeqResult[histogram] = covData
3221        SetUsedHistogramsAndPhases(GPXfile,Histo,Phases,covData,makeBack)
3222        makeBack = False
3223    SetSeqResult(GPXfile,Histograms,SeqResult)
3224
3225def DistAngle(DisAglCtls,DisAglData):
3226    import numpy.ma as ma
3227   
3228    def ShowBanner(name):
3229        print 80*'*'
3230        print '   Interatomic Distances and Angles for phase '+name
3231        print 80*'*','\n'
3232
3233    ShowBanner(DisAglCtls['Name'])
3234    SGData = DisAglData['SGData']
3235    SGtext = G2spc.SGPrint(SGData)
3236    for line in SGtext: print line
3237    Cell = DisAglData['Cell']
3238   
3239    Amat,Bmat = G2lat.cell2AB(Cell[:6])
3240    covData = {}
3241    if 'covData' in DisAglData:   
3242        covData = DisAglData['covData']
3243        covMatrix = covData['covMatrix']
3244        varyList = covData['varyList']
3245        pfx = str(DisAglData['pId'])+'::'
3246        A = G2lat.cell2A(Cell[:6])
3247        cellSig = getCellEsd(pfx,SGData,A,covData)
3248        names = [' a = ',' b = ',' c = ',' alpha = ',' beta = ',' gamma = ',' Volume = ']
3249        valEsd = [G2mth.ValEsd(Cell[i],cellSig[i],True) for i in range(7)]
3250        line = '\n Unit cell:'
3251        for name,vals in zip(names,valEsd):
3252            line += name+vals 
3253        print line
3254    else: 
3255        print '\n Unit cell: a = ','%.5f'%(Cell[0]),' b = ','%.5f'%(Cell[1]),' c = ','%.5f'%(Cell[2]), \
3256            ' alpha = ','%.3f'%(Cell[3]),' beta = ','%.3f'%(Cell[4]),' gamma = ', \
3257            '%.3f'%(Cell[5]),' volume = ','%.3f'%(Cell[6])
3258    Factor = DisAglCtls['Factors']
3259    Radii = dict(zip(DisAglCtls['AtomTypes'],zip(DisAglCtls['BondRadii'],DisAglCtls['AngleRadii'])))
3260    Units = np.array([                   #is there a nicer way to make this?
3261        [-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],
3262        [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],
3263        [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]])
3264    origAtoms = DisAglData['OrigAtoms']
3265    targAtoms = DisAglData['TargAtoms']
3266    for Oatom in origAtoms:
3267        OxyzNames = ''
3268        IndBlist = []
3269        Dist = []
3270        Vect = []
3271        VectA = []
3272        angles = []
3273        for Tatom in targAtoms:
3274            Xvcov = []
3275            TxyzNames = ''
3276            if 'covData' in DisAglData:
3277                OxyzNames = [pfx+'dAx:%d'%(Oatom[0]),pfx+'dAy:%d'%(Oatom[0]),pfx+'dAz:%d'%(Oatom[0])]
3278                TxyzNames = [pfx+'dAx:%d'%(Tatom[0]),pfx+'dAy:%d'%(Tatom[0]),pfx+'dAz:%d'%(Tatom[0])]
3279                Xvcov = G2mth.getVCov(OxyzNames+TxyzNames,varyList,covMatrix)
3280            result = G2spc.GenAtom(Tatom[3:6],SGData,False,Move=False)
3281            BsumR = (Radii[Oatom[2]][0]+Radii[Tatom[2]][0])*Factor[0]
3282            AsumR = (Radii[Oatom[2]][1]+Radii[Tatom[2]][1])*Factor[1]
3283            for Txyz,Top,Tunit in result:
3284                Dx = (Txyz-np.array(Oatom[3:6]))+Units
3285                dx = np.inner(Amat,Dx)
3286                dist = ma.masked_less(np.sqrt(np.sum(dx**2,axis=0)),0.5)
3287                IndB = ma.nonzero(ma.masked_greater(dist-BsumR,0.))
3288                if np.any(IndB):
3289                    for indb in IndB:
3290                        for i in range(len(indb)):
3291                            if str(dx.T[indb][i]) not in IndBlist:
3292                                IndBlist.append(str(dx.T[indb][i]))
3293                                unit = Units[indb][i]
3294                                tunit = '[%2d%2d%2d]'%(unit[0]+Tunit[0],unit[1]+Tunit[1],unit[2]+Tunit[2])
3295                                pdpx = G2mth.getDistDerv(Oatom[3:6],Tatom[3:6],Amat,unit,Top,SGData)
3296                                sig = 0.0
3297                                if len(Xvcov):
3298                                    sig = np.sqrt(np.inner(pdpx,np.inner(Xvcov,pdpx)))
3299                                Dist.append([Oatom[1],Tatom[1],tunit,Top,ma.getdata(dist[indb])[i],sig])
3300                                if (Dist[-1][-1]-AsumR) <= 0.:
3301                                    Vect.append(dx.T[indb][i]/Dist[-1][-2])
3302                                    VectA.append([OxyzNames,np.array(Oatom[3:6]),TxyzNames,np.array(Tatom[3:6]),unit,Top])
3303                                else:
3304                                    Vect.append([0.,0.,0.])
3305                                    VectA.append([])
3306        Vect = np.array(Vect)
3307        angles = np.zeros((len(Vect),len(Vect)))
3308        angsig = np.zeros((len(Vect),len(Vect)))
3309        for i,veca in enumerate(Vect):
3310            if np.any(veca):
3311                for j,vecb in enumerate(Vect):
3312                    if np.any(vecb):
3313                        angles[i][j],angsig[i][j] = G2mth.getAngSig(VectA[i],VectA[j],Amat,SGData,covData)
3314        line = ''
3315        for i,x in enumerate(Oatom[3:6]):
3316            if len(Xvcov):
3317                line += '%12s'%(G2mth.ValEsd(x,np.sqrt(Xvcov[i][i]),True))
3318            else:
3319                line += '%12.5f'%(x)
3320        print '\n Distances & angles for ',Oatom[1],' at ',line
3321        print 80*'*'
3322        line = ''
3323        for dist in Dist[:-1]:
3324            line += '%12s'%(dist[1].center(12))
3325        print '  To       cell +(sym. op.)      dist.  ',line
3326        for i,dist in enumerate(Dist):
3327            line = ''
3328            for j,angle in enumerate(angles[i][0:i]):
3329                sig = angsig[i][j]
3330                if angle:
3331                    if sig:
3332                        line += '%12s'%(G2mth.ValEsd(angle,sig,True).center(12))
3333                    else:
3334                        val = '%.3f'%(angle)
3335                        line += '%12s'%(val.center(12))
3336                else:
3337                    line += 12*' '
3338            if dist[5]:            #sig exists!
3339                val = G2mth.ValEsd(dist[4],dist[5])
3340            else:
3341                val = '%8.4f'%(dist[4])
3342            print %8s%10s+(%4d) %12s'%(dist[1].ljust(8),dist[2].ljust(10),dist[3],val.center(12)),line
3343
3344def DisAglTor(DATData):
3345    SGData = DATData['SGData']
3346    Cell = DATData['Cell']
3347   
3348    Amat,Bmat = G2lat.cell2AB(Cell[:6])
3349    covData = {}
3350    pfx = ''
3351    if 'covData' in DATData:   
3352        covData = DATData['covData']
3353        covMatrix = covData['covMatrix']
3354        varyList = covData['varyList']
3355        pfx = str(DATData['pId'])+'::'
3356    Datoms = []
3357    Oatoms = []
3358    for i,atom in enumerate(DATData['Datoms']):
3359        symop = atom[-1].split('+')
3360        if len(symop) == 1:
3361            symop.append('0,0,0')       
3362        symop[0] = int(symop[0])
3363        symop[1] = eval(symop[1])
3364        atom.append(symop)
3365        Datoms.append(atom)
3366        oatom = DATData['Oatoms'][i]
3367        names = ['','','']
3368        if pfx:
3369            names = [pfx+'dAx:'+str(oatom[0]),pfx+'dAy:'+str(oatom[0]),pfx+'dAz:'+str(oatom[0])]
3370        oatom += [names,]
3371        Oatoms.append(oatom)
3372    atmSeq = [atom[1]+'('+atom[-2]+')' for atom in Datoms]
3373    if DATData['Natoms'] == 4:  #torsion
3374        Tors,sig = G2mth.GetDATSig(Oatoms,Datoms,Amat,SGData,covData)
3375        print ' Torsion angle for '+DATData['Name']+' atom sequence: ',atmSeq,'=',G2mth.ValEsd(Tors,sig)
3376        print ' NB: Atom sequence determined by selection order'
3377        return      # done with torsion
3378    elif DATData['Natoms'] == 3:  #angle
3379        Ang,sig = G2mth.GetDATSig(Oatoms,Datoms,Amat,SGData,covData)
3380        print ' Angle in '+DATData['Name']+' for atom sequence: ',atmSeq,'=',G2mth.ValEsd(Ang,sig)
3381        print ' NB: Atom sequence determined by selection order'
3382    else:   #2 atoms - distance
3383        Dist,sig = G2mth.GetDATSig(Oatoms,Datoms,Amat,SGData,covData)
3384        print ' Distance in '+DATData['Name']+' for atom sequence: ',atmSeq,'=',G2mth.ValEsd(Dist,sig)
3385               
3386def BestPlane(PlaneData):
3387
3388    def ShowBanner(name):
3389        print 80*'*'
3390        print '   Best plane result for phase '+name
3391        print 80*'*','\n'
3392
3393    ShowBanner(PlaneData['Name'])
3394
3395    Cell = PlaneData['Cell']   
3396    Amat,Bmat = G2lat.cell2AB(Cell[:6])       
3397    Atoms = PlaneData['Atoms']
3398    sumXYZ = np.zeros(3)
3399    XYZ = []
3400    Natoms = len(Atoms)
3401    for atom in Atoms:
3402        xyz = np.array(atom[3:6])
3403        XYZ.append(xyz)
3404        sumXYZ += xyz
3405    sumXYZ /= Natoms
3406    XYZ = np.array(XYZ)-sumXYZ
3407    XYZ = np.inner(Amat,XYZ).T
3408    Zmat = np.zeros((3,3))
3409    for i,xyz in enumerate(XYZ):
3410        Zmat += np.outer(xyz.T,xyz)
3411    print ' Selected atoms centered at %10.5f %10.5f %10.5f'%(sumXYZ[0],sumXYZ[1],sumXYZ[2])
3412    Evec,Emat = nl.eig(Zmat)
3413    Evec = np.sqrt(Evec)/(Natoms-3)
3414    Order = np.argsort(Evec)
3415    XYZ = np.inner(XYZ,Emat.T).T
3416    XYZ = np.array([XYZ[Order[2]],XYZ[Order[1]],XYZ[Order[0]]]).T
3417    print ' Atoms in Cartesian best plane coordinates:'
3418    print ' Name         X         Y         Z'
3419    for i,xyz in enumerate(XYZ):
3420        print ' %6s%10.3f%10.3f%10.3f'%(Atoms[i][1].ljust(6),xyz[0],xyz[1],xyz[2])
3421    print '\n Best plane RMS X =%8.3f, Y =%8.3f, Z =%8.3f'%(Evec[Order[2]],Evec[Order[1]],Evec[Order[0]])   
3422
3423def main():
3424    arg = sys.argv
3425    if len(arg) > 1:
3426        GPXfile = arg[1]
3427        if not ospath.exists(GPXfile):
3428            print 'ERROR - ',GPXfile," doesn't exist!"
3429            exit()
3430        GPXpath = ospath.dirname(arg[1])
3431        Refine(GPXfile,None)
3432    else:
3433        print 'ERROR - missing filename'
3434        exit()
3435    print "Done"
3436         
3437if __name__ == '__main__':
3438    main()
Note: See TracBrowser for help on using the repository browser.