source: trunk/GSASIIstruct.py @ 792

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

change instrument parameters to dict from lists
chase down effects - got them all?
start on TOF; read TimeMap? style data - not complete

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