source: trunk/GSASIIstruct.py @ 819

Last change on this file since 819 was 819, checked in by vondreele, 9 years ago

restraint fixes

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