source: trunk/GSASIIstruct.py @ 811

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

split GSASIIgrid.py into 3 parts; GSASIIconstrGUI.py & GSASIIrestrGUI.py are new

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