source: MPbranch/GSASIIstruct.py @ 496

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

making 1st multiprocess version of GSAS-II

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