source: trunk/GSASIIstruct.py @ 917

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

final fixes to RB derivs - all seem OK now

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