source: MPbranch/GSASIIstruct.py @ 514

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

fix mp.Pool creation

  • 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-06 17:01:57 +0000 (Tue, 06 Mar 2012) $
4# $Author: toby $
5# $Revision: 514 $
6# $URL: MPbranch/GSASIIstruct.py $
7# $Id: GSASIIstruct.py 514 2012-03-06 17:01:57Z 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':1,
44                'max Rprocess':2,
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    hId = Histogram['hId']
2343    hfx = ':%d:'%(hId)
2344    bakType = calcControls[hfx+'bakType']
2345    yb = G2pwd.getBackground(hfx,parmDict,bakType,x)
2346    yc = np.zeros_like(yb)
2347       
2348    if 'C' in calcControls[hfx+'histType']:   
2349        shl = max(parmDict[hfx+'SH/L'],0.002)
2350        if hfx+'Lam1' in parmDict.keys():
2351            wave = parmDict[hfx+'Lam1']
2352            Ka2 = True
2353            lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
2354            kRatio = parmDict[hfx+'I(L2)/I(L1)']
2355        else:
2356            wave = parmDict[hfx+'Lam']
2357            Ka2 = False
2358            lamRatio = None
2359            kRatio = None
2360    else:
2361        print 'TOF Undefined at present'
2362        raise ValueError
2363    for phase in Histogram['Reflection Lists']:
2364        refList = Histogram['Reflection Lists'][phase]
2365        Phase = Phases[phase]
2366        pId = Phase['pId']
2367        pfx = '%d::'%(pId)
2368        phfx = '%d:%d:'%(pId,hId)
2369        hfx = ':%d:'%(hId)
2370        SGData = Phase['General']['SGData']
2371        A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
2372        G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
2373        GA,GB = G2lat.Gmat2AB(G)    #Orthogonalization matricies
2374        Vst = np.sqrt(nl.det(G))    #V*
2375        if 'Pawley' not in Phase['General']['Type']:
2376            refList = StructureFactor(refList,G,hfx,pfx,SGData,calcControls,parmDict)
2377        else:
2378            for refl in refList:
2379                h,k,l = refl[:3]
2380                try:
2381                    refl[9] = abs(parmDict[pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])])
2382                except KeyError:
2383                    #print ' ***Error %d,%d,%d missing from Pawley reflection list ***'%(h,k,l)
2384                    continue
2385
2386        chunk = 200
2387        if (calcControls['max Hprocess'] > 1 or
2388            calcControls['max Rprocess'] == 1 or
2389            mp.cpu_count() < 2 or len(refList) <  2 * chunk):
2390            print 'single thread function',hId,len(refList)
2391            res = ComputeReflectionProfile(
2392                (0,refList,x,Phase,calcControls,wave,G,g,GB,Vst,
2393                 parmDict,pfx,phfx,hfx,Ka2,lamRatio,kRatio))
2394            if not res: continue
2395            RefNum,refls,iBeg,iFin,yct = res
2396            yc[iBeg:iFin] += yct
2397        else:
2398            print 'multiprocess function',hId,len(refList)
2399            argList = []
2400            for i in range(0,len(refList),chunk):
2401                argList.append(
2402                    (i,refList[i:i+chunk],x,Phase,calcControls,wave,G,g,GB,Vst,
2403                     parmDict,pfx,phfx,hfx,Ka2,lamRatio,kRatio))
2404            for res in mpPool.imap_unordered(ComputeReflectionProfile,argList):
2405                if not res: continue
2406                RefNum,refls,iBeg,iFin,yct = res
2407                yc[iBeg:iFin] += yct
2408                for i,refl in enumerate(refls):
2409                    refList[RefNum+i] = refl
2410            mpPool.close()
2411            mpPool.join()
2412    return yc,yb
2413
2414def ComputeReflectionDerivative(arg):
2415    def cellVaryDerv(pfx,SGData,dpdA): 
2416        if SGData['SGLaue'] in ['-1',]:
2417            return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],
2418                [pfx+'A3',dpdA[3]],[pfx+'A4',dpdA[4]],[pfx+'A5',dpdA[5]]]
2419        elif SGData['SGLaue'] in ['2/m',]:
2420            if SGData['SGUniq'] == 'a':
2421                return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A3',dpdA[3]]]
2422            elif SGData['SGUniq'] == 'b':
2423                return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A4',dpdA[4]]]
2424            else:
2425                return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A5',dpdA[5]]]
2426        elif SGData['SGLaue'] in ['mmm',]:
2427            return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]]]
2428        elif SGData['SGLaue'] in ['4/m','4/mmm']:
2429#            return [[pfx+'A0',dpdA[0]+dpdA[1]],[pfx+'A2',dpdA[2]]]
2430            return [[pfx+'A0',dpdA[0]],[pfx+'A2',dpdA[2]]]
2431        elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
2432#            return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[3]],[pfx+'A2',dpdA[2]]]
2433            return [[pfx+'A0',dpdA[0]],[pfx+'A2',dpdA[2]]]
2434        elif SGData['SGLaue'] in ['3R', '3mR']:
2435            return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[2]],[pfx+'A3',dpdA[3]+dpdA[4]+dpdA[5]]]                       
2436        elif SGData['SGLaue'] in ['m3m','m3']:
2437#            return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[2]]]
2438            return [[pfx+'A0',dpdA[0]]]
2439
2440    (RefNum,refList,x,dependentVars,varylist,dFdvDict,
2441     Phase,calcControls,wave,G,g,GB,A,
2442     parmDict,pfx,phfx,hfx,Ka2,lamRatio,kRatio) = arg
2443    dMdv = np.zeros(shape=(len(varylist),len(x)))
2444    depDerivDict = {}
2445    for j in dependentVars:
2446        depDerivDict[j] = np.zeros(shape=(len(x)))
2447    SGData = Phase['General']['SGData']
2448    shl = max(parmDict[hfx+'SH/L'],0.002)
2449    dx = x[1]-x[0]
2450    iBegO = None
2451    for irefS,refl in enumerate(refList):
2452        iref = irefS + RefNum       
2453        if 'C' in calcControls[hfx+'histType']:        #CW powder
2454            dIdsh,dIdsp,dIdpola,dIdPO,dFdODF,dFdSA = GetIntensityDerv(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
2455            Wd,fmin,fmax = G2pwd.getWidths(refl[5],refl[6],refl[7],shl)
2456            iBeg = np.searchsorted(x,refl[5]-fmin)
2457            iFin = np.searchsorted(x,refl[5]+fmax)
2458            if iBegO is None:
2459                iBegO = iBeg
2460                iFinO = iFin
2461            else:
2462                iBegO = min(iBegO,iBeg)
2463                iFinO = max(iFinO,iFin)
2464            if not iBeg+iFin:       #peak below low limit - skip peak
2465                continue
2466                #return None
2467            elif not iBeg-iFin:     #peak above high limit - done
2468                break
2469                #return None
2470            pos = refl[5]
2471            tanth = tand(pos/2.0)
2472            costh = cosd(pos/2.0)
2473            lenBF = iFin-iBeg
2474            dMdpk = np.zeros(shape=(6,lenBF))
2475            dMdipk = G2pwd.getdFCJVoigt3(refl[5],refl[6],refl[7],shl,x[iBeg:iFin])
2476            for i in range(1,5):
2477                dMdpk[i] += 100.*dx*refl[13]*refl[9]*dMdipk[i]
2478            dMdpk[0] += 100.*dx*refl[13]*refl[9]*dMdipk[0]
2479            dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':np.zeros_like(dMdpk[0])}
2480            if Ka2:
2481                pos2 = refl[5]+lamRatio*tanth       # + 360/pi * Dlam/lam * tan(th)
2482                kdelt = int((pos2-refl[5])/dx)               
2483                iBeg2 = min(lenX,iBeg+kdelt)
2484                iFin2 = min(lenX,iFin+kdelt)
2485                if iBeg2-iFin2:
2486                    iBegO = min(iBegO,iBeg2)
2487                    iFinO = max(iFinO,iFin2)
2488                    lenBF2 = iFin2-iBeg2
2489                    dMdpk2 = np.zeros(shape=(6,lenBF2))
2490                    dMdipk2 = G2pwd.getdFCJVoigt3(pos2,refl[6],refl[7],shl,x[iBeg2:iFin2])
2491                    for i in range(1,5):
2492                        dMdpk2[i] = 100.*dx*refl[13]*refl[9]*kRatio*dMdipk2[i]
2493                    dMdpk2[0] = 100.*dx*refl[13]*refl[9]*kRatio*dMdipk2[0]
2494                    dMdpk2[5] = 100.*dx*refl[13]*dMdipk2[0]
2495                    dervDict2 = {'int':dMdpk2[0],'pos':dMdpk2[1],'sig':dMdpk2[2],'gam':dMdpk2[3],'shl':dMdpk2[4],'L1/L2':dMdpk2[5]*refl[9]}
2496            if 'Pawley' in Phase['General']['Type']:
2497                    try:
2498                        idx = varylist.index(pfx+'PWLref:'+str(pawleyLookup[pfx+'%d,%d,%d'%refl[:3]]))
2499                        dMdv[idx][iBeg:iFin] = dervDict['int']/refl[9]
2500                        if Ka2:
2501                            dMdv[idx][iBeg2:iFin2] = dervDict2['int']/refl[9]
2502                        # Assuming Pawley variables not in constraints
2503                    except ValueError:
2504                        pass
2505            dpdA,dpdw,dpdZ,dpdSh,dpdTr,dpdX,dpdY = GetReflPosDerv(refl,wave,A,hfx,calcControls,parmDict)
2506            names = {hfx+'Scale':[dIdsh,'int'],hfx+'Polariz.':[dIdpola,'int'],phfx+'Scale':[dIdsp,'int'],
2507                    hfx+'U':[tanth**2,'sig'],hfx+'V':[tanth,'sig'],hfx+'W':[1.0,'sig'],
2508                    hfx+'X':[1.0/costh,'gam'],hfx+'Y':[tanth,'gam'],hfx+'SH/L':[1.0,'shl'],
2509                    hfx+'I(L2)/I(L1)':[1.0,'L1/L2'],hfx+'Zero':[dpdZ,'pos'],hfx+'Lam':[dpdw,'pos'],
2510                    hfx+'Shift':[dpdSh,'pos'],hfx+'Transparency':[dpdTr,'pos'],hfx+'DisplaceX':[dpdX,'pos'],
2511                    hfx+'DisplaceY':[dpdY,'pos'],}
2512            for name in names:
2513                item = names[name]
2514                if name in varylist:
2515                    dMdv[varylist.index(name)][iBeg:iFin] += item[0]*dervDict[item[1]]
2516                    if Ka2:
2517                        dMdv[varylist.index(name)][iBeg2:iFin2] += item[0]*dervDict2[item[1]]
2518                elif name in dependentVars:
2519                    if Ka2:
2520                        depDerivDict[name][iBeg2:iFin2] += item[0]*dervDict2[item[1]]
2521                    depDerivDict[name][iBeg:iFin] += item[0]*dervDict[item[1]]
2522            for iPO in dIdPO:
2523                if iPO in varylist:
2524                    dMdv[varylist.index(iPO)][iBeg:iFin] += dIdPO[iPO]*dervDict['int']
2525                    if Ka2:
2526                        dMdv[varylist.index(iPO)][iBeg2:iFin2] += dIdPO[iPO]*dervDict2['int']
2527                elif iPO in dependentVars:
2528                    depDerivDict[iPO][iBeg:iFin] += dIdPO[iPO]*dervDict['int']
2529                    if Ka2:
2530                        depDerivDict[iPO][iBeg2:iFin2] += dIdPO[iPO]*dervDict2['int']
2531            for i,name in enumerate(['omega','chi','phi']):
2532                aname = pfx+'SH '+name
2533                if aname in varylist:
2534                    dMdv[varylist.index(aname)][iBeg:iFin] += dFdSA[i]*dervDict['int']
2535                    if Ka2:
2536                        dMdv[varylist.index(aname)][iBeg2:iFin2] += dFdSA[i]*dervDict2['int']
2537                elif aname in dependentVars:
2538                    depDerivDict[aname][iBeg:iFin] += dFdSA[i]*dervDict['int']
2539                    if Ka2:
2540                        depDerivDict[aname][iBeg2:iFin2] += dFdSA[i]*dervDict2['int']
2541            for iSH in dFdODF:
2542                if iSH in varylist:
2543                    dMdv[varylist.index(iSH)][iBeg:iFin] += dFdODF[iSH]*dervDict['int']
2544                    if Ka2:
2545                        dMdv[varylist.index(iSH)][iBeg2:iFin2] += dFdODF[iSH]*dervDict2['int']
2546                elif iSH in dependentVars:
2547                    depDerivDict[iSH][iBeg:iFin] += dFdODF[iSH]*dervDict['int']
2548                    if Ka2:
2549                        depDerivDict[iSH][iBeg2:iFin2] += dFdODF[iSH]*dervDict2['int']
2550            cellDervNames = cellVaryDerv(pfx,SGData,dpdA)
2551            for name,dpdA in cellDervNames:
2552                if name in varylist:
2553                    dMdv[varylist.index(name)][iBeg:iFin] += dpdA*dervDict['pos']
2554                    if Ka2:
2555                        dMdv[varylist.index(name)][iBeg2:iFin2] += dpdA*dervDict2['pos']
2556                elif name in dependentVars:
2557                    depDerivDict[name][iBeg:iFin] += dpdA*dervDict['pos']
2558                    if Ka2:
2559                        depDerivDict[name][iBeg2:iFin2] += dpdA*dervDict2['pos']
2560            dDijDict = GetHStrainShiftDerv(refl,SGData,phfx)
2561            for name in dDijDict:
2562                if name in varylist:
2563                    dMdv[varylist.index(name)][iBeg:iFin] += dDijDict[name]*dervDict['pos']
2564                    if Ka2:
2565                        dMdv[varylist.index(name)][iBeg2:iFin2] += dDijDict[name]*dervDict2['pos']
2566                elif name in dependentVars:
2567                    depDerivDict[name][iBeg:iFin] += dDijDict[name]*dervDict['pos']
2568                    if Ka2:
2569                        depDerivDict[name][iBeg2:iFin2] += dDijDict[name]*dervDict2['pos']
2570            gamDict = GetSampleGamDerv(refl,wave,G,GB,phfx,calcControls,parmDict)
2571            for name in gamDict:
2572                if name in varylist:
2573                    dMdv[varylist.index(name)][iBeg:iFin] += gamDict[name]*dervDict['gam']
2574                    if Ka2:
2575                        dMdv[varylist.index(name)][iBeg2:iFin2] += gamDict[name]*dervDict2['gam']
2576                elif name in dependentVars:
2577                    depDerivDict[name][iBeg:iFin] += gamDict[name]*dervDict['gam']
2578                    if Ka2:
2579                        depDerivDict[name][iBeg2:iFin2] += gamDict[name]*dervDict2['gam']
2580                                               
2581        elif 'T' in calcControls[hfx+'histType']:
2582            print 'TOF Undefined at present'
2583            raise Exception    #no TOF yet
2584        #do atom derivatives -  for F,X & U so far             
2585        corr = dervDict['int']/refl[9]
2586        if Ka2:
2587            corr2 = dervDict2['int']/refl[9]
2588        for name in varylist+dependentVars:
2589            try:
2590                aname = name.split(pfx)[1][:2]
2591                if aname not in ['Af','dA','AU']: continue # skip anything not an atom param
2592            except IndexError:
2593                continue
2594            if name in varylist:
2595                dMdv[varylist.index(name)][iBeg:iFin] += dFdvDict[name][iref]*corr
2596                if Ka2:
2597                    dMdv[varylist.index(name)][iBeg2:iFin2] += dFdvDict[name][iref]*corr2
2598            elif name in dependentVars:
2599                depDerivDict[name][iBeg:iFin] += dFdvDict[name][iref]*corr
2600                if Ka2:
2601                    depDerivDict[name][iBeg2:iFin2] += dFdvDict[name][iref]*corr2
2602        # end loop over reflections
2603    return depDerivDict,dMdv[:,iBegO:iFinO],iBegO,iFinO
2604       
2605def getPowderProfileDerv(parmDict,x,varylist,Histogram,Phases,calcControls,pawleyLookup):   
2606    if (calcControls['max Hprocess'] <= 1 and calcControls['max Rprocess'] > 1:
2607        mpPool = mp.Pool()
2608    # create a list of dependent variables and set up a dictionary to hold their derivatives
2609    dependentVars = G2mv.GetDependentVars()
2610    depDerivDict = {}
2611    for j in dependentVars:
2612        depDerivDict[j] = np.zeros(shape=(len(x)))
2613    #print 'dependent vars',dependentVars
2614    lenX = len(x)               
2615    hId = Histogram['hId']
2616    hfx = ':%d:'%(hId)
2617    bakType = calcControls[hfx+'bakType']
2618    dMdv = np.zeros(shape=(len(varylist),len(x)))
2619    dMdb,dMddb = G2pwd.getBackgroundDerv(hfx,parmDict,bakType,x)
2620    if hfx+'Back:0' in varylist: # for now assume that Back:x vars to not appear in constraints
2621        bBpos =varylist.index(hfx+'Back:0')
2622        dMdv[bBpos:bBpos+len(dMdb)] = dMdb
2623    names = [hfx+'DebyeA',hfx+'DebyeR',hfx+'DebyeU']
2624    for name in varylist:
2625        if 'Debye' in name:
2626            id = int(name.split(':')[-1])
2627            parm = name[:int(name.rindex(':'))]
2628            ip = names.index(parm)
2629            dMdv[varylist.index(name)] = dMddb[3*id+ip]
2630    if 'C' in calcControls[hfx+'histType']:   
2631        dx = x[1]-x[0]
2632        shl = max(parmDict[hfx+'SH/L'],0.002)
2633        if hfx+'Lam1' in parmDict.keys():
2634            wave = parmDict[hfx+'Lam1']
2635            Ka2 = True
2636            lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
2637            kRatio = parmDict[hfx+'I(L2)/I(L1)']
2638        else:
2639            wave = parmDict[hfx+'Lam']
2640            Ka2 = False
2641            lamRatio = None
2642            kRatio = None
2643    else:
2644        print 'TOF Undefined at present'
2645        raise ValueError
2646    for phase in Histogram['Reflection Lists']:
2647        refList = Histogram['Reflection Lists'][phase]
2648        Phase = Phases[phase]
2649        SGData = Phase['General']['SGData']
2650        pId = Phase['pId']
2651        pfx = '%d::'%(pId)
2652        phfx = '%d:%d:'%(pId,hId)
2653        A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
2654        G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
2655        GA,GB = G2lat.Gmat2AB(G)    #Orthogonalization matricies
2656        if 'Pawley' not in Phase['General']['Type']:
2657            dFdvDict = StructureFactorDerv(refList,G,hfx,pfx,SGData,calcControls,parmDict)
2658        else:
2659            for iref,refl in enumerate(refList):
2660                try:
2661                    refl[9] = abs(parmDict[pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%refl[:3]])])
2662                except KeyError:
2663                    #print ' ***Error %d,%d,%d missing from Pawley reflection list ***'%refl[:3]
2664                    continue
2665        chunk = 200
2666        if (calcControls['max Hprocess'] > 1 or
2667            calcControls['max Rprocess'] == 1 or
2668            mp.cpu_count() < 2 or len(refList) <  2 * chunk):
2669            print 'single thread derivs'
2670            TdepDerivDict, TdMdv, iBegO, iFinO = ComputeReflectionDerivative(
2671                (0,refList,x,dependentVars,varylist,dFdvDict,
2672                 Phase,calcControls,wave,G,g,GB,A,
2673                 parmDict,pfx,phfx,hfx,Ka2,lamRatio,kRatio))
2674            dMdv[:,iBegO:iFinO] += TdMdv
2675            for j in dependentVars:         
2676                depDerivDict[j] += TdepDerivDict[j]
2677        else:
2678            print 'multiprocess derivs'
2679            argList = []
2680            for iref in range(0,len(refList),chunk):
2681                argList.append(
2682                    (iref,refList[iref:iref+chunk],x,dependentVars,varylist,dFdvDict,
2683                     Phase,calcControls,wave,G,g,GB,A,
2684                     parmDict,pfx,phfx,hfx,Ka2,lamRatio,kRatio))
2685            for res in mpPool.imap_unordered(ComputeReflectionDerivative,argList):
2686                if res is None: continue
2687                TdepDerivDict, TdMdv, iBegO, iFinO = res
2688                dMdv[:,iBegO:iFinO] += TdMdv
2689                for j in dependentVars:
2690                    depDerivDict[j] += TdepDerivDict[j]
2691
2692    # now process derivatives in constraints
2693    G2mv.Dict2Deriv(varylist,depDerivDict,dMdv)
2694    return dMdv
2695
2696def dervRefine(values,HistoPhases,parmdict,varylist,calcControls,pawleyLookup,dlg):
2697    parmdict.update(zip(varylist,values))
2698    G2mv.Dict2Map(parmdict,varylist)
2699    Histograms,Phases = HistoPhases
2700    nvar = len(varylist)
2701    dMdv = np.empty(0)
2702    histoList = Histograms.keys()
2703    histoList.sort()
2704    for histogram in histoList:
2705        if 'PWDR' in histogram[:4]:
2706            Histogram = Histograms[histogram]
2707            hId = Histogram['hId']
2708            hfx = ':%d:'%(hId)
2709            Limits = calcControls[hfx+'Limits']
2710            x,y,w,yc,yb,yd = Histogram['Data']
2711            xB = np.searchsorted(x,Limits[0])
2712            xF = np.searchsorted(x,Limits[1])
2713            dMdvh = np.sqrt(w[xB:xF])*getPowderProfileDerv(parmdict,x[xB:xF],
2714                varylist,Histogram,Phases,calcControls,pawleyLookup)
2715            if len(dMdv):
2716                dMdv = np.concatenate((dMdv.T,dMdvh.T)).T
2717            else:
2718                dMdv = dMdvh
2719    return dMdv
2720
2721def ComputePowderHessian(args):
2722    Histogram,parmdict,varylist,Phases,calcControls,pawleyLookup = args
2723    hId = Histogram['hId']
2724    hfx = ':%d:'%(hId)
2725    Limits = calcControls[hfx+'Limits']
2726    x,y,w,yc,yb,yd = Histogram['Data']
2727    dy = y-yc
2728    xB = np.searchsorted(x,Limits[0])
2729    xF = np.searchsorted(x,Limits[1])
2730    dMdvh = np.sqrt(w[xB:xF])*getPowderProfileDerv(
2731        parmdict,x[xB:xF],
2732        varylist,Histogram,Phases,calcControls,pawleyLookup)
2733    Vec = np.sum(dMdvh*np.sqrt(w[xB:xF])*dy[xB:xF],axis=1)
2734    Hess = np.inner(dMdvh,dMdvh)
2735    return Vec,Hess
2736
2737#def HessRefine(values,HistoPhases,parmdict,varylist,calcControls,pawleyLookup,dlg):
2738#    parmdict.update(zip(varylist,values))
2739#    G2mv.Dict2Map(parmdict,varylist)
2740#    Histograms,Phases = HistoPhases
2741#    nvar = len(varylist)
2742#    Hess = np.empty(0)
2743#    histoList = Histograms.keys()
2744#    histoList.sort()
2745#    for histogram in histoList:
2746#        if 'PWDR' in histogram[:4]:
2747#            Histogram = Histograms[histogram]
2748#            hId = Histogram['hId']
2749#            hfx = ':%d:'%(hId)
2750#            Limits = calcControls[hfx+'Limits']
2751#            x,y,w,yc,yb,yd = Histogram['Data']
2752#            dy = y-yc
2753#            xB = np.searchsorted(x,Limits[0])
2754#            xF = np.searchsorted(x,Limits[1])
2755#            dMdvh = np.sqrt(w[xB:xF])*getPowderProfileDerv(parmdict,x[xB:xF],
2756#                varylist,Histogram,Phases,calcControls,pawleyLookup)
2757#            if dlg:
2758#                dlg.Update(Histogram['wRp'],newmsg='Hessian for histogram %d Rwp=%8.3f%s'%(hId,Histogram['wRp'],'%'))[0]
2759#            if len(Hess):
2760#                Vec += np.sum(dMdvh*np.sqrt(w[xB:xF])*dy[xB:xF],axis=1)
2761#                Hess += np.inner(dMdvh,dMdvh)
2762#            else:
2763#                Vec = np.sum(dMdvh*np.sqrt(w[xB:xF])*dy[xB:xF],axis=1)
2764#                Hess = np.inner(dMdvh,dMdvh)
2765#    return Vec,Hess
2766
2767def HessRefine(values,HistoPhases,parmdict,varylist,calcControls,pawleyLookup,dlg):
2768    parmdict.update(zip(varylist,values))
2769    G2mv.Dict2Map(parmdict,varylist)
2770    Histograms,Phases = HistoPhases
2771    nvar = len(varylist)
2772    HessSum = np.zeros((nvar,nvar))
2773    VecSum = np.zeros(nvar)
2774    histoList = Histograms.keys()
2775    histoList.sort()
2776    argList = []
2777    MaxProcess = calcControls['max Hprocess']
2778    for histogram in histoList:
2779        if 'PWDR' in histogram[:4]:
2780            Histogram = Histograms[histogram]
2781            argList.append([
2782                Histogram,parmdict,varylist,Phases,
2783                calcControls,pawleyLookup])
2784    if MaxProcess > 1:
2785        mpPool = mp.Pool()
2786        for Vec,Hess in mpPool.imap_unordered(ComputePowderHessian,argList):
2787            VecSum += Vec
2788            HessSum += Hess
2789        mpPool.close()
2790        mpPool.join()
2791    else:
2792        for arg in argList:
2793            Vec,Hess = ComputePowderHessian(arg)
2794            VecSum += Vec
2795            HessSum += Hess
2796    return VecSum,HessSum
2797
2798def ComputePowderProfile(args):
2799    Histogram,parmdict,varylist,Phases,calcControls,pawleyLookup,histkey = args
2800    hId = Histogram['hId']
2801    hfx = ':%d:'%(hId)
2802    x,y,w,yc,yb,yd = Histogram['Data']
2803    Limits = calcControls[hfx+'Limits']
2804    xB = np.searchsorted(x,Limits[0])
2805    xF = np.searchsorted(x,Limits[1])
2806    yc,yb = getPowderProfile(parmdict,x[xB:xF],
2807                            varylist,Histogram,Phases,calcControls,
2808                            pawleyLookup)
2809    return xB,xF,yc,yb,Histogram['Reflection Lists'],histkey
2810
2811#def errRefine(values,HistoPhases,parmdict,varylist,calcControls,pawleyLookup,dlg):       
2812#    parmdict.update(zip(varylist,values))
2813#    Values2Dict(parmdict, varylist, values)
2814#    G2mv.Dict2Map(parmdict,varylist)
2815#    Histograms,Phases = HistoPhases
2816#    M = np.empty(0)
2817#    sumwYo = 0
2818#    Nobs = 0
2819#    histoList = Histograms.keys()
2820#    histoList.sort()
2821#    for histogram in histoList:
2822#        if 'PWDR' in histogram[:4]:
2823#            Histogram = Histograms[histogram]
2824#            hId = Histogram['hId']
2825#            hfx = ':%d:'%(hId)
2826#            Limits = calcControls[hfx+'Limits']
2827#            x,y,w,yc,yb,yd = Histogram['Data']
2828#            yc *= 0.0                           #zero full calcd profiles
2829#            yb *= 0.0
2830#            yd *= 0.0
2831#            xB = np.searchsorted(x,Limits[0])
2832#            xF = np.searchsorted(x,Limits[1])
2833#            Histogram['Nobs'] = xF-xB
2834#            Nobs += Histogram['Nobs']
2835#            Histogram['sumwYo'] = np.sum(w[xB:xF]*y[xB:xF]**2)
2836#            sumwYo += Histogram['sumwYo']
2837#            yc[xB:xF],yb[xB:xF] = getPowderProfile(parmdict,x[xB:xF],
2838#                varylist,Histogram,Phases,calcControls,pawleyLookup)
2839#            yc[xB:xF] += yb[xB:xF]
2840#            yd[xB:xF] = y[xB:xF]-yc[xB:xF]
2841#            Histogram['sumwYd'] = np.sum(np.sqrt(w[xB:xF])*(yd[xB:xF]))
2842#            wdy = -np.sqrt(w[xB:xF])*(yd[xB:xF])
2843#            Histogram['wRp'] = min(100.,np.sqrt(np.sum(wdy**2)/Histogram['sumwYo'])*100.)
2844#            if dlg:
2845#                dlg.Update(Histogram['wRp'],newmsg='For histogram %d Rwp=%8.3f%s'%(hId,Histogram['wRp'],'%'))[0]
2846#            M = np.concatenate((M,wdy))
2847#    Histograms['sumwYo'] = sumwYo
2848#    Histograms['Nobs'] = Nobs
2849#    Rwp = min(100.,np.sqrt(np.sum(M**2)/sumwYo)*100.)
2850#    if dlg:
2851#        GoOn = dlg.Update(Rwp,newmsg='%s%8.3f%s'%('Powder profile Rwp =',Rwp,'%'))[0]
2852#        if not GoOn:
2853#            parmDict['saved values'] = values
2854#            raise Exception         #Abort!!
2855#    return M
2856#   
2857def errRefine(values,HistoPhases,parmdict,varylist,calcControls,pawleyLookup,dlg):       
2858    parmdict.update(zip(varylist,values))
2859    Values2Dict(parmdict, varylist, values)
2860    G2mv.Dict2Map(parmdict,varylist)
2861    Histograms,Phases = HistoPhases
2862    M = np.empty(0)
2863    sumwYo = 0
2864    Nobs = 0
2865    histoList = Histograms.keys()
2866    histoList.sort()
2867    argList = []
2868    MaxProcess = calcControls['max Hprocess']
2869    for histogram in histoList:
2870        if 'PWDR' in histogram[:4]:
2871            Histogram = Histograms[histogram]
2872            argList.append(
2873                [Histogram,parmdict,varylist,Phases,calcControls,pawleyLookup,histogram]
2874                )
2875    if MaxProcess > 1: 
2876        mpPool = mp.Pool()
2877        for (xB,xF,ycSect,ybSect,RL,histkey
2878             ) in mpPool.imap_unordered(ComputePowderProfile,argList):
2879            Histogram = Histograms[histkey]
2880            Histogram['Reflection Lists'] = RL
2881            x,y,w,yc,yb,yd = Histogram['Data']
2882            yc *= 0.0                           #zero full calcd profiles
2883            yb *= 0.0
2884            yd *= 0.0
2885            Histogram['Nobs'] = xF-xB
2886            Nobs += Histogram['Nobs']
2887            Histogram['sumwYo'] = np.sum(w[xB:xF]*y[xB:xF]**2)
2888            sumwYo += Histogram['sumwYo']
2889           
2890            yc[xB:xF] = ycSect
2891            yb[xB:xF] = ybSect
2892            yc[xB:xF] += yb[xB:xF]
2893            yd[xB:xF] = y[xB:xF]-yc[xB:xF]
2894            Histogram['sumwYd'] = np.sum(np.sqrt(w[xB:xF])*(yd[xB:xF]))
2895            wdy = -np.sqrt(w[xB:xF])*(yd[xB:xF])
2896            Histogram['wRp'] = min(100.,np.sqrt(np.sum(wdy**2)/Histogram['sumwYo'])*100.)
2897            M = np.concatenate((M,wdy))
2898        mpPool.close()
2899        mpPool.join()
2900    else:
2901        for arg in argList:
2902            xB,xF,ycSect,ybSect,RL,histkey = ComputePowderProfile(arg)
2903            Histogram = arg[0]
2904            hId = Histogram['hId']
2905            x,y,w,yc,yb,yd = Histogram['Data']
2906            yc *= 0.0                           #zero full calcd profiles
2907            yb *= 0.0
2908            yd *= 0.0
2909            Histogram['Nobs'] = xF-xB
2910            Nobs += Histogram['Nobs']
2911            Histogram['sumwYo'] = np.sum(w[xB:xF]*y[xB:xF]**2)
2912            sumwYo += Histogram['sumwYo']
2913           
2914            yc[xB:xF] = ycSect
2915            yb[xB:xF] = ybSect
2916            yc[xB:xF] += yb[xB:xF]
2917            yd[xB:xF] = y[xB:xF]-yc[xB:xF]
2918            Histogram['sumwYd'] = np.sum(np.sqrt(w[xB:xF])*(yd[xB:xF]))
2919            wdy = -np.sqrt(w[xB:xF])*(yd[xB:xF])
2920            Histogram['wRp'] = min(100.,np.sqrt(np.sum(wdy**2)/Histogram['sumwYo'])*100.)
2921            if dlg:
2922                dlg.Update(Histogram['wRp'],newmsg='For histogram %d Rwp=%8.3f%s'%(hId,Histogram['wRp'],'%'))[0]
2923            M = np.concatenate((M,wdy))
2924
2925    Histograms['sumwYo'] = sumwYo
2926    Histograms['Nobs'] = Nobs
2927    Rwp = min(100.,np.sqrt(np.sum(M**2)/sumwYo)*100.)
2928    if dlg:
2929        GoOn = dlg.Update(Rwp,newmsg='%s%8.3f%s'%('Powder profile Rwp =',Rwp,'%'))[0]
2930        if not GoOn:
2931            parmDict['saved values'] = values
2932            raise Exception         #Abort!!
2933    return M       
2934                   
2935def Refine(GPXfile,dlg):
2936    import pytexture as ptx
2937    ptx.pyqlmninit()            #initialize fortran arrays for spherical harmonics
2938   
2939    ShowBanner()
2940    varyList = []
2941    parmDict = {}
2942    G2mv.InitVars()   
2943    Controls = GetControls(GPXfile)
2944    ShowControls(Controls)
2945    calcControls = {}
2946    calcControls.update(Controls)           
2947    constrDict,constrFlag,fixedList = GetConstraints(GPXfile)
2948    Histograms,Phases = GetUsedHistogramsAndPhases(GPXfile)
2949    if not Phases:
2950        print ' *** ERROR - you have no histograms to refine! ***'
2951        print ' *** Refine aborted ***'
2952        raise Exception
2953    if not Histograms:
2954        print ' *** ERROR - you have no data to refine with! ***'
2955        print ' *** Refine aborted ***'
2956        raise Exception       
2957    Natoms,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables = GetPhaseData(Phases)
2958    calcControls['Natoms'] = Natoms
2959    calcControls['FFtables'] = FFtables
2960    calcControls['BLtables'] = BLtables
2961    hapVary,hapDict,controlDict = GetHistogramPhaseData(Phases,Histograms)
2962    calcControls.update(controlDict)
2963    histVary,histDict,controlDict = GetHistogramData(Histograms)
2964    calcControls.update(controlDict)
2965    varyList = phaseVary+hapVary+histVary
2966    parmDict.update(phaseDict)
2967    parmDict.update(hapDict)
2968    parmDict.update(histDict)
2969    GetFprime(calcControls,Histograms)
2970    # do constraint processing
2971    try:
2972        groups,parmlist = G2mv.GroupConstraints(constrDict)
2973        G2mv.GenerateConstraints(groups,parmlist,varyList,constrDict,constrFlag,fixedList)
2974    except:
2975        print ' *** ERROR - your constraints are internally inconsistent ***'
2976        raise Exception(' *** Refine aborted ***')
2977    # check to see which generated parameters are fully varied
2978    msg = G2mv.SetVaryFlags(varyList)
2979    if msg:
2980        print ' *** ERROR - you have not set the refine flags for constraints consistently! ***'
2981        print msg
2982        raise Exception(' *** Refine aborted ***')
2983    G2mv.Map2Dict(parmDict,varyList)
2984#    print G2mv.VarRemapShow(varyList)
2985
2986    while True:
2987        begin = time.time()
2988        values =  np.array(Dict2Values(parmDict, varyList))
2989        Ftol = Controls['min dM/M']
2990        Factor = Controls['shift factor']
2991        maxCyc = Controls['max cyc']
2992        if 'Jacobian' in Controls['deriv type']:           
2993            result = so.leastsq(errRefine,values,Dfun=dervRefine,full_output=True,
2994                ftol=Ftol,col_deriv=True,factor=Factor,
2995                args=([Histograms,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg))
2996            ncyc = int(result[2]['nfev']/2)
2997        elif 'Hessian' in Controls['deriv type']:
2998            result = G2mth.HessianLSQ(errRefine,values,Hess=HessRefine,ftol=Ftol,maxcyc=maxCyc,
2999                args=([Histograms,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg))
3000            ncyc = result[2]['num cyc']+1                           
3001        else:           #'numeric'
3002            result = so.leastsq(errRefine,values,full_output=True,ftol=Ftol,epsfcn=1.e-8,factor=Factor,
3003                args=([Histograms,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg))
3004            ncyc = int(result[2]['nfev']/len(varyList))
3005#        table = dict(zip(varyList,zip(values,result[0],(result[0]-values))))
3006#        for item in table: print item,table[item]               #useful debug - are things shifting?
3007        runtime = time.time()-begin
3008        chisq = np.sum(result[2]['fvec']**2)
3009        Values2Dict(parmDict, varyList, result[0])
3010        G2mv.Dict2Map(parmDict,varyList)
3011       
3012        Rwp = np.sqrt(chisq/Histograms['sumwYo'])*100.      #to %
3013        GOF = chisq/(Histograms['Nobs']-len(varyList))
3014        print '\n Refinement results:'
3015        print 135*'-'
3016        print ' Number of function calls:',result[2]['nfev'],' Number of observations: ',Histograms['Nobs'],' Number of parameters: ',len(varyList)
3017        print ' Refinement time = %8.3fs, %8.3fs/cycle, for %d cycles'%(runtime,runtime/ncyc,ncyc)
3018        print ' wRp = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f'%(Rwp,chisq,GOF)
3019        try:
3020            covMatrix = result[1]*GOF
3021            sig = np.sqrt(np.diag(covMatrix))
3022            if np.any(np.isnan(sig)):
3023                print '*** Least squares aborted - some invalid esds possible ***'
3024#            table = dict(zip(varyList,zip(values,result[0],(result[0]-values)/sig)))
3025#            for item in table: print item,table[item]               #useful debug - are things shifting?
3026            break                   #refinement succeeded - finish up!
3027        except TypeError:          #result[1] is None on singular matrix
3028            print '**** Refinement failed - singular matrix ****'
3029            if 'Hessian' in Controls['deriv type']:
3030                for i in result[2]['psing'].reverse():
3031                        print 'Removing parameter: ',varyList[i]
3032                        del(varyList[i])                   
3033            else:
3034                Ipvt = result[2]['ipvt']
3035                for i,ipvt in enumerate(Ipvt):
3036                    if not np.sum(result[2]['fjac'],axis=1)[i]:
3037                        print 'Removing parameter: ',varyList[ipvt-1]
3038                        del(varyList[ipvt-1])
3039                        break
3040
3041#    print 'dependentParmList: ',G2mv.dependentParmList
3042#    print 'arrayList: ',G2mv.arrayList
3043#    print 'invarrayList: ',G2mv.invarrayList
3044#    print 'indParmList: ',G2mv.indParmList
3045#    print 'fixedDict: ',G2mv.fixedDict
3046#    print 'test1'
3047    GetFobsSq(Histograms,Phases,parmDict,calcControls)
3048#    print 'test2'
3049    sigDict = dict(zip(varyList,sig))
3050    newCellDict = GetNewCellParms(parmDict,varyList)
3051    newAtomDict = ApplyXYZshifts(parmDict,varyList)
3052    covData = {'variables':result[0],'varyList':varyList,'sig':sig,
3053        'covMatrix':covMatrix,'title':GPXfile,'newAtomDict':newAtomDict,'newCellDict':newCellDict}
3054    # add the uncertainties into the esd dictionary (sigDict)
3055    sigDict.update(G2mv.ComputeDepESD(covMatrix,varyList,parmDict))
3056    SetPhaseData(parmDict,sigDict,Phases,covData)
3057    SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms)
3058    SetHistogramData(parmDict,sigDict,Histograms)
3059    G2mv.PrintIndependentVars(parmDict,varyList,sigDict)
3060    SetUsedHistogramsAndPhases(GPXfile,Histograms,Phases,covData)
3061   
3062#for testing purposes!!!
3063#    import cPickle
3064#    file = open('structTestdata.dat','wb')
3065#    cPickle.dump(parmDict,file,1)
3066#    cPickle.dump(varyList,file,1)
3067#    for histogram in Histograms:
3068#        if 'PWDR' in histogram[:4]:
3069#            Histogram = Histograms[histogram]
3070#    cPickle.dump(Histogram,file,1)
3071#    cPickle.dump(Phases,file,1)
3072#    cPickle.dump(calcControls,file,1)
3073#    cPickle.dump(pawleyLookup,file,1)
3074#    file.close()
3075
3076    ShowControls(Controls)
3077    print 'Last Refinement time = %8.3fs'%(runtime)
3078    if dlg:
3079        return Rwp
3080
3081def SeqRefine(GPXfile,dlg):
3082    import pytexture as ptx
3083    ptx.pyqlmninit()            #initialize fortran arrays for spherical harmonics
3084   
3085    ShowBanner()
3086    print ' Sequential Refinement'
3087    G2mv.InitVars()   
3088    Controls = GetControls(GPXfile)
3089    ShowControls(Controls)           
3090    constrDict,constrFlag,fixedList = GetConstraints(GPXfile)
3091    Histograms,Phases = GetUsedHistogramsAndPhases(GPXfile)
3092    if not Phases:
3093        print ' *** ERROR - you have no histograms to refine! ***'
3094        print ' *** Refine aborted ***'
3095        raise Exception
3096    if not Histograms:
3097        print ' *** ERROR - you have no data to refine with! ***'
3098        print ' *** Refine aborted ***'
3099        raise Exception
3100    Natoms,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables = GetPhaseData(Phases,False)
3101    if 'Seq Data' in Controls:
3102        histNames = Controls['Seq Data']
3103    else:
3104        histNames = GetHistogramNames(GPXfile,['PWDR',])
3105    if 'Reverse Seq' in Controls:
3106        if Controls['Reverse Seq']:
3107            histNames.reverse()
3108    SeqResult = {'histNames':histNames}
3109    makeBack = True
3110    for ihst,histogram in enumerate(histNames):
3111        ifPrint = False
3112        if dlg:
3113            dlg.SetTitle('Residual for histogram '+str(ihst))
3114        calcControls = {}
3115        calcControls['Natoms'] = Natoms
3116        calcControls['FFtables'] = FFtables
3117        calcControls['BLtables'] = BLtables
3118        varyList = []
3119        parmDict = {}
3120        Histo = {histogram:Histograms[histogram],}
3121        hapVary,hapDict,controlDict = GetHistogramPhaseData(Phases,Histo,False)
3122        calcControls.update(controlDict)
3123        histVary,histDict,controlDict = GetHistogramData(Histo,False)
3124        calcControls.update(controlDict)
3125        varyList = phaseVary+hapVary+histVary
3126        if not ihst:
3127            saveVaryList = varyList[:]
3128            for i,item in enumerate(saveVaryList):
3129                items = item.split(':')
3130                if items[1]:
3131                    items[1] = ''
3132                item = ':'.join(items)
3133                saveVaryList[i] = item
3134            SeqResult['varyList'] = saveVaryList
3135        else:
3136            newVaryList = varyList[:]
3137            for i,item in enumerate(newVaryList):
3138                items = item.split(':')
3139                if items[1]:
3140                    items[1] = ''
3141                item = ':'.join(items)
3142                newVaryList[i] = item
3143            if newVaryList != SeqResult['varyList']:
3144                print newVaryList
3145                print SeqResult['varyList']
3146                print '**** ERROR - variable list for this histogram does not match previous'
3147                raise Exception
3148        parmDict.update(phaseDict)
3149        parmDict.update(hapDict)
3150        parmDict.update(histDict)
3151        GetFprime(calcControls,Histo)
3152        constrDict,constrFlag,fixedList = G2mv.InputParse([])        #constraints go here?
3153        groups,parmlist = G2mv.GroupConstraints(constrDict)
3154        G2mv.GenerateConstraints(groups,parmlist,varyList,constrDict,constrFlag,fixedList)
3155        G2mv.Map2Dict(parmDict,varyList)
3156   
3157        while True:
3158            begin = time.time()
3159            values =  np.array(Dict2Values(parmDict, varyList))
3160            Ftol = Controls['min dM/M']
3161            Factor = Controls['shift factor']
3162            maxCyc = Controls['max cyc']
3163
3164            if 'Jacobian' in Controls['deriv type']:           
3165                result = so.leastsq(errRefine,values,Dfun=dervRefine,full_output=True,
3166                    ftol=Ftol,col_deriv=True,factor=Factor,
3167                    args=([Histo,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg))
3168                ncyc = int(result[2]['nfev']/2)
3169            elif 'Hessian' in Controls['deriv type']:
3170                result = G2mth.HessianLSQ(errRefine,values,Hess=HessRefine,ftol=Ftol,maxcyc=maxCyc,
3171                    args=([Histo,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg))
3172                ncyc = result[2]['num cyc']+1                           
3173            else:           #'numeric'
3174                result = so.leastsq(errRefine,values,full_output=True,ftol=Ftol,epsfcn=1.e-8,factor=Factor,
3175                    args=([Histo,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg))
3176                ncyc = int(result[2]['nfev']/len(varyList))
3177
3178
3179
3180            runtime = time.time()-begin
3181            chisq = np.sum(result[2]['fvec']**2)
3182            Values2Dict(parmDict, varyList, result[0])
3183            G2mv.Dict2Map(parmDict,varyList)
3184           
3185            Rwp = np.sqrt(chisq/Histo['sumwYo'])*100.      #to %
3186            GOF = chisq/(Histo['Nobs']-len(varyList))
3187            print '\n Refinement results for histogram: v'+histogram
3188            print 135*'-'
3189            print ' Number of function calls:',result[2]['nfev'],' Number of observations: ',Histo['Nobs'],' Number of parameters: ',len(varyList)
3190            print ' Refinement time = %8.3fs, %8.3fs/cycle, for %d cycles'%(runtime,runtime/ncyc,ncyc)
3191            print ' wRp = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f'%(Rwp,chisq,GOF)
3192            try:
3193                covMatrix = result[1]*GOF
3194                sig = np.sqrt(np.diag(covMatrix))
3195                if np.any(np.isnan(sig)):
3196                    print '*** Least squares aborted - some invalid esds possible ***'
3197                    ifPrint = True
3198                break                   #refinement succeeded - finish up!
3199            except TypeError:          #result[1] is None on singular matrix
3200                print '**** Refinement failed - singular matrix ****'
3201                if 'Hessian' in Controls['deriv type']:
3202                    for i in result[2]['psing'].reverse():
3203                            print 'Removing parameter: ',varyList[i]
3204                            del(varyList[i])                   
3205                else:
3206                    Ipvt = result[2]['ipvt']
3207                    for i,ipvt in enumerate(Ipvt):
3208                        if not np.sum(result[2]['fjac'],axis=1)[i]:
3209                            print 'Removing parameter: ',varyList[ipvt-1]
3210                            del(varyList[ipvt-1])
3211                            break
3212   
3213        GetFobsSq(Histo,Phases,parmDict,calcControls)
3214        sigDict = dict(zip(varyList,sig))
3215        newCellDict = GetNewCellParms(parmDict,varyList)
3216        newAtomDict = ApplyXYZshifts(parmDict,varyList)
3217        covData = {'variables':result[0],'varyList':varyList,'sig':sig,
3218            'covMatrix':covMatrix,'title':histogram,'newAtomDict':newAtomDict,'newCellDict':newCellDict}
3219        SetHistogramPhaseData(parmDict,sigDict,Phases,Histo,ifPrint)
3220        SetHistogramData(parmDict,sigDict,Histo,ifPrint)
3221        SeqResult[histogram] = covData
3222        SetUsedHistogramsAndPhases(GPXfile,Histo,Phases,covData,makeBack)
3223        makeBack = False
3224    SetSeqResult(GPXfile,Histograms,SeqResult)
3225
3226def DistAngle(DisAglCtls,DisAglData):
3227    import numpy.ma as ma
3228   
3229    def ShowBanner(name):
3230        print 80*'*'
3231        print '   Interatomic Distances and Angles for phase '+name
3232        print 80*'*','\n'
3233
3234    ShowBanner(DisAglCtls['Name'])
3235    SGData = DisAglData['SGData']
3236    SGtext = G2spc.SGPrint(SGData)
3237    for line in SGtext: print line
3238    Cell = DisAglData['Cell']
3239   
3240    Amat,Bmat = G2lat.cell2AB(Cell[:6])
3241    covData = {}
3242    if 'covData' in DisAglData:   
3243        covData = DisAglData['covData']
3244        covMatrix = covData['covMatrix']
3245        varyList = covData['varyList']
3246        pfx = str(DisAglData['pId'])+'::'
3247        A = G2lat.cell2A(Cell[:6])
3248        cellSig = getCellEsd(pfx,SGData,A,covData)
3249        names = [' a = ',' b = ',' c = ',' alpha = ',' beta = ',' gamma = ',' Volume = ']
3250        valEsd = [G2mth.ValEsd(Cell[i],cellSig[i],True) for i in range(7)]
3251        line = '\n Unit cell:'
3252        for name,vals in zip(names,valEsd):
3253            line += name+vals 
3254        print line
3255    else: 
3256        print '\n Unit cell: a = ','%.5f'%(Cell[0]),' b = ','%.5f'%(Cell[1]),' c = ','%.5f'%(Cell[2]), \
3257            ' alpha = ','%.3f'%(Cell[3]),' beta = ','%.3f'%(Cell[4]),' gamma = ', \
3258            '%.3f'%(Cell[5]),' volume = ','%.3f'%(Cell[6])
3259    Factor = DisAglCtls['Factors']
3260    Radii = dict(zip(DisAglCtls['AtomTypes'],zip(DisAglCtls['BondRadii'],DisAglCtls['AngleRadii'])))
3261    Units = np.array([                   #is there a nicer way to make this?
3262        [-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],
3263        [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],
3264        [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]])
3265    origAtoms = DisAglData['OrigAtoms']
3266    targAtoms = DisAglData['TargAtoms']
3267    for Oatom in origAtoms:
3268        OxyzNames = ''
3269        IndBlist = []
3270        Dist = []
3271        Vect = []
3272        VectA = []
3273        angles = []
3274        for Tatom in targAtoms:
3275            Xvcov = []
3276            TxyzNames = ''
3277            if 'covData' in DisAglData:
3278                OxyzNames = [pfx+'dAx:%d'%(Oatom[0]),pfx+'dAy:%d'%(Oatom[0]),pfx+'dAz:%d'%(Oatom[0])]
3279                TxyzNames = [pfx+'dAx:%d'%(Tatom[0]),pfx+'dAy:%d'%(Tatom[0]),pfx+'dAz:%d'%(Tatom[0])]
3280                Xvcov = G2mth.getVCov(OxyzNames+TxyzNames,varyList,covMatrix)
3281            result = G2spc.GenAtom(Tatom[3:6],SGData,False,Move=False)
3282            BsumR = (Radii[Oatom[2]][0]+Radii[Tatom[2]][0])*Factor[0]
3283            AsumR = (Radii[Oatom[2]][1]+Radii[Tatom[2]][1])*Factor[1]
3284            for Txyz,Top,Tunit in result:
3285                Dx = (Txyz-np.array(Oatom[3:6]))+Units
3286                dx = np.inner(Amat,Dx)
3287                dist = ma.masked_less(np.sqrt(np.sum(dx**2,axis=0)),0.5)
3288                IndB = ma.nonzero(ma.masked_greater(dist-BsumR,0.))
3289                if np.any(IndB):
3290                    for indb in IndB:
3291                        for i in range(len(indb)):
3292                            if str(dx.T[indb][i]) not in IndBlist:
3293                                IndBlist.append(str(dx.T[indb][i]))
3294                                unit = Units[indb][i]
3295                                tunit = '[%2d%2d%2d]'%(unit[0]+Tunit[0],unit[1]+Tunit[1],unit[2]+Tunit[2])
3296                                pdpx = G2mth.getDistDerv(Oatom[3:6],Tatom[3:6],Amat,unit,Top,SGData)
3297                                sig = 0.0
3298                                if len(Xvcov):
3299                                    sig = np.sqrt(np.inner(pdpx,np.inner(Xvcov,pdpx)))
3300                                Dist.append([Oatom[1],Tatom[1],tunit,Top,ma.getdata(dist[indb])[i],sig])
3301                                if (Dist[-1][-1]-AsumR) <= 0.:
3302                                    Vect.append(dx.T[indb][i]/Dist[-1][-2])
3303                                    VectA.append([OxyzNames,np.array(Oatom[3:6]),TxyzNames,np.array(Tatom[3:6]),unit,Top])
3304                                else:
3305                                    Vect.append([0.,0.,0.])
3306                                    VectA.append([])
3307        Vect = np.array(Vect)
3308        angles = np.zeros((len(Vect),len(Vect)))
3309        angsig = np.zeros((len(Vect),len(Vect)))
3310        for i,veca in enumerate(Vect):
3311            if np.any(veca):
3312                for j,vecb in enumerate(Vect):
3313                    if np.any(vecb):
3314                        angles[i][j],angsig[i][j] = G2mth.getAngSig(VectA[i],VectA[j],Amat,SGData,covData)
3315        line = ''
3316        for i,x in enumerate(Oatom[3:6]):
3317            if len(Xvcov):
3318                line += '%12s'%(G2mth.ValEsd(x,np.sqrt(Xvcov[i][i]),True))
3319            else:
3320                line += '%12.5f'%(x)
3321        print '\n Distances & angles for ',Oatom[1],' at ',line
3322        print 80*'*'
3323        line = ''
3324        for dist in Dist[:-1]:
3325            line += '%12s'%(dist[1].center(12))
3326        print '  To       cell +(sym. op.)      dist.  ',line
3327        for i,dist in enumerate(Dist):
3328            line = ''
3329            for j,angle in enumerate(angles[i][0:i]):
3330                sig = angsig[i][j]
3331                if angle:
3332                    if sig:
3333                        line += '%12s'%(G2mth.ValEsd(angle,sig,True).center(12))
3334                    else:
3335                        val = '%.3f'%(angle)
3336                        line += '%12s'%(val.center(12))
3337                else:
3338                    line += 12*' '
3339            if dist[5]:            #sig exists!
3340                val = G2mth.ValEsd(dist[4],dist[5])
3341            else:
3342                val = '%8.4f'%(dist[4])
3343            print %8s%10s+(%4d) %12s'%(dist[1].ljust(8),dist[2].ljust(10),dist[3],val.center(12)),line
3344
3345def DisAglTor(DATData):
3346    SGData = DATData['SGData']
3347    Cell = DATData['Cell']
3348   
3349    Amat,Bmat = G2lat.cell2AB(Cell[:6])
3350    covData = {}
3351    pfx = ''
3352    if 'covData' in DATData:   
3353        covData = DATData['covData']
3354        covMatrix = covData['covMatrix']
3355        varyList = covData['varyList']
3356        pfx = str(DATData['pId'])+'::'
3357    Datoms = []
3358    Oatoms = []
3359    for i,atom in enumerate(DATData['Datoms']):
3360        symop = atom[-1].split('+')
3361        if len(symop) == 1:
3362            symop.append('0,0,0')       
3363        symop[0] = int(symop[0])
3364        symop[1] = eval(symop[1])
3365        atom.append(symop)
3366        Datoms.append(atom)
3367        oatom = DATData['Oatoms'][i]
3368        names = ['','','']
3369        if pfx:
3370            names = [pfx+'dAx:'+str(oatom[0]),pfx+'dAy:'+str(oatom[0]),pfx+'dAz:'+str(oatom[0])]
3371        oatom += [names,]
3372        Oatoms.append(oatom)
3373    atmSeq = [atom[1]+'('+atom[-2]+')' for atom in Datoms]
3374    if DATData['Natoms'] == 4:  #torsion
3375        Tors,sig = G2mth.GetDATSig(Oatoms,Datoms,Amat,SGData,covData)
3376        print ' Torsion angle for '+DATData['Name']+' atom sequence: ',atmSeq,'=',G2mth.ValEsd(Tors,sig)
3377        print ' NB: Atom sequence determined by selection order'
3378        return      # done with torsion
3379    elif DATData['Natoms'] == 3:  #angle
3380        Ang,sig = G2mth.GetDATSig(Oatoms,Datoms,Amat,SGData,covData)
3381        print ' Angle in '+DATData['Name']+' for atom sequence: ',atmSeq,'=',G2mth.ValEsd(Ang,sig)
3382        print ' NB: Atom sequence determined by selection order'
3383    else:   #2 atoms - distance
3384        Dist,sig = G2mth.GetDATSig(Oatoms,Datoms,Amat,SGData,covData)
3385        print ' Distance in '+DATData['Name']+' for atom sequence: ',atmSeq,'=',G2mth.ValEsd(Dist,sig)
3386               
3387def BestPlane(PlaneData):
3388
3389    def ShowBanner(name):
3390        print 80*'*'
3391        print '   Best plane result for phase '+name
3392        print 80*'*','\n'
3393
3394    ShowBanner(PlaneData['Name'])
3395
3396    Cell = PlaneData['Cell']   
3397    Amat,Bmat = G2lat.cell2AB(Cell[:6])       
3398    Atoms = PlaneData['Atoms']
3399    sumXYZ = np.zeros(3)
3400    XYZ = []
3401    Natoms = len(Atoms)
3402    for atom in Atoms:
3403        xyz = np.array(atom[3:6])
3404        XYZ.append(xyz)
3405        sumXYZ += xyz
3406    sumXYZ /= Natoms
3407    XYZ = np.array(XYZ)-sumXYZ
3408    XYZ = np.inner(Amat,XYZ).T
3409    Zmat = np.zeros((3,3))
3410    for i,xyz in enumerate(XYZ):
3411        Zmat += np.outer(xyz.T,xyz)
3412    print ' Selected atoms centered at %10.5f %10.5f %10.5f'%(sumXYZ[0],sumXYZ[1],sumXYZ[2])
3413    Evec,Emat = nl.eig(Zmat)
3414    Evec = np.sqrt(Evec)/(Natoms-3)
3415    Order = np.argsort(Evec)
3416    XYZ = np.inner(XYZ,Emat.T).T
3417    XYZ = np.array([XYZ[Order[2]],XYZ[Order[1]],XYZ[Order[0]]]).T
3418    print ' Atoms in Cartesian best plane coordinates:'
3419    print ' Name         X         Y         Z'
3420    for i,xyz in enumerate(XYZ):
3421        print ' %6s%10.3f%10.3f%10.3f'%(Atoms[i][1].ljust(6),xyz[0],xyz[1],xyz[2])
3422    print '\n Best plane RMS X =%8.3f, Y =%8.3f, Z =%8.3f'%(Evec[Order[2]],Evec[Order[1]],Evec[Order[0]])   
3423
3424def main():
3425    arg = sys.argv
3426    if len(arg) > 1:
3427        GPXfile = arg[1]
3428        if not ospath.exists(GPXfile):
3429            print 'ERROR - ',GPXfile," doesn't exist!"
3430            exit()
3431        GPXpath = ospath.dirname(arg[1])
3432        Refine(GPXfile,None)
3433    else:
3434        print 'ERROR - missing filename'
3435        exit()
3436    print "Done"
3437         
3438if __name__ == '__main__':
3439    main()
Note: See TracBrowser for help on using the repository browser.