source: trunk/GSASIIstruct.py @ 494

Last change on this file since 494 was 494, checked in by vondreele, 11 years ago

brians version

  • Property svn:keywords set to Date Author Revision URL Id
File size: 140.8 KB
Line 
1#GSASIIstructure - structure computation routines
2
3########### SVN repository information ###################
4
5# $Date: 2012-02-24 20:25:49 +0000 (Fri, 24 Feb 2012) $
6
7# $Author: vondreele $
8
9# $Revision: 494 $
10
11# $URL: trunk/GSASIIstruct.py $
12
13# $Id: GSASIIstruct.py 494 2012-02-24 20:25:49Z vondreele $
14
15########### SVN repository information ###################
16
17MaxThreads = 1
18
19MaxProcess = 2
20
21
22
23import sys
24
25import os
26
27import os.path as ospath
28
29import time
30
31import math
32
33import cPickle
34
35import multiprocessing as mp
36
37import numpy as np
38
39import numpy.linalg as nl
40
41import scipy.optimize as so
42
43import GSASIIpath
44
45import GSASIIElem as G2el
46
47import GSASIIlattice as G2lat
48
49import GSASIIspc as G2spc
50
51import GSASIIpwd as G2pwd
52
53import GSASIImapvars as G2mv
54
55import GSASIImath as G2mth
56
57#import buggerwx
58
59try:
60
61    import mkl
62
63except:
64
65    print "MKL module not present"
66
67
68
69
70
71sind = lambda x: np.sin(x*np.pi/180.)
72
73cosd = lambda x: np.cos(x*np.pi/180.)
74
75tand = lambda x: np.tan(x*np.pi/180.)
76
77asind = lambda x: 180.*np.arcsin(x)/np.pi
78
79acosd = lambda x: 180.*np.arccos(x)/np.pi
80
81atan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
82
83
84
85
86
87def GetControls(GPXfile):
88
89    ''' Returns dictionary of control items found in GSASII gpx file
90
91    input:
92
93        GPXfile = .gpx full file name
94
95    return:
96
97        Controls = dictionary of control items
98
99    '''
100
101    Controls = {'deriv type':'analytical','min dM/M':0.0001,'shift factor':1.}
102
103    file = open(GPXfile,'rb')
104
105    while True:
106
107        try:
108
109            data = cPickle.load(file)
110
111        except EOFError:
112
113            break
114
115        datum = data[0]
116
117        if datum[0] == 'Controls':
118
119            Controls.update(datum[1])
120
121    file.close()
122
123    return Controls
124
125   
126
127def GetConstraints(GPXfile):
128
129    constList = []
130
131    file = open(GPXfile,'rb')
132
133    while True:
134
135        try:
136
137            data = cPickle.load(file)
138
139        except EOFError:
140
141            break
142
143        datum = data[0]
144
145        if datum[0] == 'Constraints':
146
147            constDict = datum[1]
148
149            for item in constDict:
150
151                constList += constDict[item]
152
153    file.close()
154
155    constDict = []
156
157    constFlag = []
158
159    fixedList = []
160
161    for item in constList:
162
163        if item[-2]:
164
165            fixedList.append(str(item[-2]))
166
167        else:
168
169            fixedList.append('0')
170
171        if item[-1]:
172
173            constFlag.append(['VARY'])
174
175        else:
176
177            constFlag.append([])
178
179        itemDict = {}
180
181        for term in item[:-2]:
182
183            itemDict[term[1]] = term[0]
184
185        constDict.append(itemDict)
186
187    return constDict,constFlag,fixedList
188
189   
190
191def GetPhaseNames(GPXfile):
192
193    ''' Returns a list of phase names found under 'Phases' in GSASII gpx file
194
195    input:
196
197        GPXfile = gpx full file name
198
199    return:
200
201        PhaseNames = list of phase names
202
203    '''
204
205    file = open(GPXfile,'rb')
206
207    PhaseNames = []
208
209    while True:
210
211        try:
212
213            data = cPickle.load(file)
214
215        except EOFError:
216
217            break
218
219        datum = data[0]
220
221        if 'Phases' == datum[0]:
222
223            for datus in data[1:]:
224
225                PhaseNames.append(datus[0])
226
227    file.close()
228
229    return PhaseNames
230
231
232
233def GetAllPhaseData(GPXfile,PhaseName):
234
235    ''' Returns the entire dictionary for PhaseName from GSASII gpx file
236
237    input:
238
239        GPXfile = gpx full file name
240
241        PhaseName = phase name
242
243    return:
244
245        phase dictionary
246
247    '''       
248
249    file = open(GPXfile,'rb')
250
251    General = {}
252
253    Atoms = []
254
255    while True:
256
257        try:
258
259            data = cPickle.load(file)
260
261        except EOFError:
262
263            break
264
265        datum = data[0]
266
267        if 'Phases' == datum[0]:
268
269            for datus in data[1:]:
270
271                if datus[0] == PhaseName:
272
273                    break
274
275    file.close()
276
277    return datus[1]
278
279   
280
281def GetHistograms(GPXfile,hNames):
282
283    """ Returns a dictionary of histograms found in GSASII gpx file
284
285    input:
286
287        GPXfile = .gpx full file name
288
289        hNames = list of histogram names
290
291    return:
292
293        Histograms = dictionary of histograms (types = PWDR & HKLF)
294
295    """
296
297    file = open(GPXfile,'rb')
298
299    Histograms = {}
300
301    while True:
302
303        try:
304
305            data = cPickle.load(file)
306
307        except EOFError:
308
309            break
310
311        datum = data[0]
312
313        hist = datum[0]
314
315        if hist in hNames:
316
317            if 'PWDR' in hist[:4]:
318
319                PWDRdata = {}
320
321                PWDRdata['Data'] = datum[1][1]          #powder data arrays
322
323                PWDRdata[data[2][0]] = data[2][1]       #Limits
324
325                PWDRdata[data[3][0]] = data[3][1]       #Background
326
327                PWDRdata[data[4][0]] = data[4][1]       #Instrument parameters
328
329                PWDRdata[data[5][0]] = data[5][1]       #Sample parameters
330
331                try:
332
333                    PWDRdata[data[9][0]] = data[9][1]       #Reflection lists might be missing
334
335                except IndexError:
336
337                    PWDRdata['Reflection lists'] = {}
338
339   
340
341                Histograms[hist] = PWDRdata
342
343            elif 'HKLF' in hist[:4]:
344
345                HKLFdata = []
346
347                datum = data[0]
348
349                HKLFdata = datum[1:][0]
350
351                Histograms[hist] = HKLFdata           
352
353    file.close()
354
355    return Histograms
356
357   
358
359def GetHistogramNames(GPXfile,hType):
360
361    """ Returns a list of histogram names found in GSASII gpx file
362
363    input:
364
365        GPXfile = .gpx full file name
366
367        hType = list ['PWDR','HKLF']
368
369    return:
370
371        HistogramNames = list of histogram names (types = PWDR & HKLF)
372
373    """
374
375    file = open(GPXfile,'rb')
376
377    HistogramNames = []
378
379    while True:
380
381        try:
382
383            data = cPickle.load(file)
384
385        except EOFError:
386
387            break
388
389        datum = data[0]
390
391        if datum[0][:4] in hType:
392
393            HistogramNames.append(datum[0])
394
395    file.close()
396
397    return HistogramNames
398
399   
400
401def GetUsedHistogramsAndPhases(GPXfile):
402
403    ''' Returns all histograms that are found in any phase
404
405    and any phase that uses a histogram
406
407    input:
408
409        GPXfile = .gpx full file name
410
411    return:
412
413        Histograms = dictionary of histograms as {name:data,...}
414
415        Phases = dictionary of phases that use histograms
416
417    '''
418
419    phaseNames = GetPhaseNames(GPXfile)
420
421    histoList = GetHistogramNames(GPXfile,['PWDR','HKLF'])
422
423    allHistograms = GetHistograms(GPXfile,histoList)
424
425    phaseData = {}
426
427    for name in phaseNames: 
428
429        phaseData[name] =  GetAllPhaseData(GPXfile,name)
430
431    Histograms = {}
432
433    Phases = {}
434
435    for phase in phaseData:
436
437        Phase = phaseData[phase]
438
439        if Phase['Histograms']:
440
441            if phase not in Phases:
442
443                pId = phaseNames.index(phase)
444
445                Phase['pId'] = pId
446
447                Phases[phase] = Phase
448
449            for hist in Phase['Histograms']:
450
451                if hist not in Histograms:
452
453                    Histograms[hist] = allHistograms[hist]
454
455                    #future restraint, etc. histograms here           
456
457                    hId = histoList.index(hist)
458
459                    Histograms[hist]['hId'] = hId
460
461    return Histograms,Phases
462
463   
464
465def GPXBackup(GPXfile,makeBack=True):
466
467    import distutils.file_util as dfu
468
469    GPXpath,GPXname = ospath.split(GPXfile)
470
471    if GPXpath == '': GPXpath = '.'
472
473    Name = ospath.splitext(GPXname)[0]
474
475    files = os.listdir(GPXpath)
476
477    last = 0
478
479    for name in files:
480
481        name = name.split('.')
482
483        if len(name) >= 3 and name[0] == Name and 'bak' in name[-2]:
484
485            if makeBack:
486
487                last = max(last,int(name[-2].strip('bak'))+1)
488
489            else:
490
491                last = max(last,int(name[-2].strip('bak')))
492
493    GPXback = ospath.join(GPXpath,GPXname.rstrip('.'.join(name[-2:]))+'.bak'+str(last)+'.gpx')
494
495    dfu.copy_file(GPXfile,GPXback)
496
497    return GPXback
498
499       
500
501def SetUsedHistogramsAndPhases(GPXfile,Histograms,Phases,CovData,makeBack=True):
502
503    ''' Updates gpxfile from all histograms that are found in any phase
504
505    and any phase that used a histogram
506
507    input:
508
509        GPXfile = .gpx full file name
510
511        Histograms = dictionary of histograms as {name:data,...}
512
513        Phases = dictionary of phases that use histograms
514
515        CovData = dictionary of refined variables, varyList, & covariance matrix
516
517        makeBack = True if new backup of .gpx file is to be made; else use the last one made
518
519    '''
520
521                       
522
523    GPXback = GPXBackup(GPXfile,makeBack)
524
525    print '\n',135*'-'
526
527    print 'Read from file:',GPXback
528
529    print 'Save to file  :',GPXfile
530
531    infile = open(GPXback,'rb')
532
533    outfile = open(GPXfile,'wb')
534
535    while True:
536
537        try:
538
539            data = cPickle.load(infile)
540
541        except EOFError:
542
543            break
544
545        datum = data[0]
546
547#        print 'read: ',datum[0]
548
549        if datum[0] == 'Phases':
550
551            for iphase in range(len(data)):
552
553                if data[iphase][0] in Phases:
554
555                    phaseName = data[iphase][0]
556
557                    data[iphase][1].update(Phases[phaseName])
558
559        elif datum[0] == 'Covariance':
560
561            data[0][1] = CovData
562
563        try:
564
565            histogram = Histograms[datum[0]]
566
567#            print 'found ',datum[0]
568
569            data[0][1][1] = histogram['Data']
570
571            for datus in data[1:]:
572
573#                print '    read: ',datus[0]
574
575                if datus[0] in ['Background','Instrument Parameters','Sample Parameters','Reflection Lists']:
576
577                    datus[1] = histogram[datus[0]]
578
579        except KeyError:
580
581            pass
582
583                               
584
585        cPickle.dump(data,outfile,1)
586
587    infile.close()
588
589    outfile.close()
590
591    print 'GPX file save successful'
592
593   
594
595def SetSeqResult(GPXfile,Histograms,SeqResult):
596
597    GPXback = GPXBackup(GPXfile)
598
599    print '\n',135*'-'
600
601    print 'Read from file:',GPXback
602
603    print 'Save to file  :',GPXfile
604
605    infile = open(GPXback,'rb')
606
607    outfile = open(GPXfile,'wb')
608
609    while True:
610
611        try:
612
613            data = cPickle.load(infile)
614
615        except EOFError:
616
617            break
618
619        datum = data[0]
620
621        if datum[0] == 'Sequental results':
622
623            data[0][1] = SeqResult
624
625        try:
626
627            histogram = Histograms[datum[0]]
628
629            data[0][1][1] = histogram['Data']
630
631            for datus in data[1:]:
632
633                if datus[0] in ['Background','Instrument Parameters','Sample Parameters','Reflection Lists']:
634
635                    datus[1] = histogram[datus[0]]
636
637        except KeyError:
638
639            pass
640
641                               
642
643        cPickle.dump(data,outfile,1)
644
645    infile.close()
646
647    outfile.close()
648
649    print 'GPX file save successful'
650
651                       
652
653def GetPWDRdata(GPXfile,PWDRname):
654
655    ''' Returns powder data from GSASII gpx file
656
657    input:
658
659        GPXfile = .gpx full file name
660
661        PWDRname = powder histogram name as obtained from GetHistogramNames
662
663    return:
664
665        PWDRdata = powder data dictionary with:
666
667            Data - powder data arrays, Limits, Instrument Parameters, Sample Parameters
668
669       
670
671    '''
672
673    file = open(GPXfile,'rb')
674
675    PWDRdata = {}
676
677    while True:
678
679        try:
680
681            data = cPickle.load(file)
682
683        except EOFError:
684
685            break
686
687        datum = data[0]
688
689        if datum[0] == PWDRname:
690
691            PWDRdata['Data'] = datum[1][1]          #powder data arrays
692
693            PWDRdata[data[2][0]] = data[2][1]       #Limits
694
695            PWDRdata[data[3][0]] = data[3][1]       #Background
696
697            PWDRdata[data[4][0]] = data[4][1]       #Instrument parameters
698
699            PWDRdata[data[5][0]] = data[5][1]       #Sample parameters
700
701            try:
702
703                PWDRdata[data[9][0]] = data[9][1]       #Reflection lists might be missing
704
705            except IndexError:
706
707                PWDRdata['Reflection lists'] = {}
708
709    file.close()
710
711    return PWDRdata
712
713   
714
715def GetHKLFdata(GPXfile,HKLFname):
716
717    ''' Returns single crystal data from GSASII gpx file
718
719    input:
720
721        GPXfile = .gpx full file name
722
723        HKLFname = single crystal histogram name as obtained from GetHistogramNames
724
725    return:
726
727        HKLFdata = single crystal data list of reflections: for each reflection:
728
729            HKLF = [np.array([h,k,l]),FoSq,sigFoSq,FcSq,Fcp,Fcpp,phase]
730
731    '''
732
733    file = open(GPXfile,'rb')
734
735    HKLFdata = []
736
737    while True:
738
739        try:
740
741            data = cPickle.load(file)
742
743        except EOFError:
744
745            break
746
747        datum = data[0]
748
749        if datum[0] == HKLFname:
750
751            HKLFdata = datum[1:][0]
752
753    file.close()
754
755    return HKLFdata
756
757   
758
759def ShowBanner():
760
761    print 80*'*'
762
763    print '   General Structure Analysis System-II Crystal Structure Refinement'
764
765    print '     by Robert B. Von Dreele, Argonne National Laboratory(C), 2010'
766
767    print ' This product includes software developed by the UChicago Argonne, LLC,' 
768
769    print '            as Operator of Argonne National Laboratory.'
770
771    print 80*'*','\n'
772
773
774
775def ShowControls(Controls):
776
777    print ' Least squares controls:'
778
779    print ' Derivative type: ',Controls['deriv type']
780
781    print ' Minimum delta-M/M for convergence: ','%.2g'%(Controls['min dM/M'])
782
783    print ' Initial shift factor: ','%.3f'%(Controls['shift factor'])
784
785   
786
787def GetFFtable(General):
788
789    ''' returns a dictionary of form factor data for atom types found in General
790
791    input:
792
793        General = dictionary of phase info.; includes AtomTypes
794
795    return:
796
797        FFtable = dictionary of form factor data; key is atom type
798
799    '''
800
801    atomTypes = General['AtomTypes']
802
803    FFtable = {}
804
805    for El in atomTypes:
806
807        FFs = G2el.GetFormFactorCoeff(El.split('+')[0].split('-')[0])
808
809        for item in FFs:
810
811            if item['Symbol'] == El.upper():
812
813                FFtable[El] = item
814
815    return FFtable
816
817   
818
819def GetBLtable(General):
820
821    ''' returns a dictionary of neutron scattering length data for atom types & isotopes found in General
822
823    input:
824
825        General = dictionary of phase info.; includes AtomTypes & Isotopes
826
827    return:
828
829        BLtable = dictionary of scattering length data; key is atom type
830
831    '''
832
833    atomTypes = General['AtomTypes']
834
835    BLtable = {}
836
837    isotopes = General['Isotopes']
838
839    isotope = General['Isotope']
840
841    for El in atomTypes:
842
843        BLtable[El] = [isotope[El],isotopes[El][isotope[El]]]
844
845    return BLtable
846
847       
848
849def GetPawleyConstr(SGLaue,PawleyRef,pawleyVary):
850
851    if SGLaue in ['-1','2/m','mmm']:
852
853        return                      #no Pawley symmetry required constraints
854
855    for i,varyI in enumerate(pawleyVary):
856
857        refI = int(varyI.split(':')[-1])
858
859        ih,ik,il = PawleyRef[refI][:3]
860
861        for varyJ in pawleyVary[0:i]:
862
863            refJ = int(varyJ.split(':')[-1])
864
865            jh,jk,jl = PawleyRef[refJ][:3]
866
867            if SGLaue in ['4/m','4/mmm']:
868
869                isum = ih**2+ik**2
870
871                jsum = jh**2+jk**2
872
873                if abs(il) == abs(jl) and isum == jsum:
874
875                    G2mv.StoreEquivalence(varyJ,(varyI,))
876
877            elif SGLaue in ['3R','3mR']:
878
879                isum = ih**2+ik**2+il**2
880
881                jsum = jh**2+jk**2*jl**2
882
883                isum2 = ih*ik+ih*il+ik*il
884
885                jsum2 = jh*jk+jh*jl+jk*jl
886
887                if isum == jsum and isum2 == jsum2:
888
889                    G2mv.StoreEquivalence(varyJ,(varyI,))
890
891            elif SGLaue in ['3','3m1','31m','6/m','6/mmm']:
892
893                isum = ih**2+ik**2+ih*ik
894
895                jsum = jh**2+jk**2+jh*jk
896
897                if abs(il) == abs(jl) and isum == jsum:
898
899                    G2mv.StoreEquivalence(varyJ,(varyI,))
900
901            elif SGLaue in ['m3','m3m']:
902
903                isum = ih**2+ik**2+il**2
904
905                jsum = jh**2+jk**2+jl**2
906
907                if isum == jsum:
908
909                    G2mv.StoreEquivalence(varyJ,(varyI,))
910
911                   
912
913def cellVary(pfx,SGData): 
914
915    if SGData['SGLaue'] in ['-1',]:
916
917        return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3',pfx+'A4',pfx+'A5']
918
919    elif SGData['SGLaue'] in ['2/m',]:
920
921        if SGData['SGUniq'] == 'a':
922
923            return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3']
924
925        elif SGData['SGUniq'] == 'b':
926
927            return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A4']
928
929        else:
930
931            return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A5']
932
933    elif SGData['SGLaue'] in ['mmm',]:
934
935        return [pfx+'A0',pfx+'A1',pfx+'A2']
936
937    elif SGData['SGLaue'] in ['4/m','4/mmm']:
938
939        return [pfx+'A0',pfx+'A2']
940
941    elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
942
943        return [pfx+'A0',pfx+'A2']
944
945    elif SGData['SGLaue'] in ['3R', '3mR']:
946
947        return [pfx+'A0',pfx+'A3']                       
948
949    elif SGData['SGLaue'] in ['m3m','m3']:
950
951        return [pfx+'A0',]
952
953                   
954
955def GetPhaseData(PhaseData,Print=True):
956
957           
958
959    def PrintFFtable(FFtable):
960
961        print '\n X-ray scattering factors:'
962
963        print '   Symbol     fa                                      fb                                      fc'
964
965        print 99*'-'
966
967        for Ename in FFtable:
968
969            ffdata = FFtable[Ename]
970
971            fa = ffdata['fa']
972
973            fb = ffdata['fb']
974
975            print ' %8s %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f' %  \
976
977                (Ename.ljust(8),fa[0],fa[1],fa[2],fa[3],fb[0],fb[1],fb[2],fb[3],ffdata['fc'])
978
979               
980
981    def PrintBLtable(BLtable):
982
983        print '\n Neutron scattering factors:'
984
985        print '   Symbol   isotope       mass       b       resonant terms'
986
987        print 99*'-'
988
989        for Ename in BLtable:
990
991            bldata = BLtable[Ename]
992
993            isotope = bldata[0]
994
995            mass = bldata[1][0]
996
997            blen = bldata[1][1]
998
999            bres = []
1000
1001            if len(bldata[1]) > 2:
1002
1003                bres = bldata[1][2:]
1004
1005            line = ' %8s%11s %10.3f %8.3f'%(Ename.ljust(8),isotope.center(11),mass,blen)
1006
1007            for item in bres:
1008
1009                line += '%10.5g'%(item)
1010
1011            print line
1012
1013               
1014
1015    def PrintAtoms(General,Atoms):
1016
1017        print '\n Atoms:'
1018
1019        line = '   name    type  refine?   x         y         z    '+ \
1020
1021            '  frac site sym  mult I/A   Uiso     U11     U22     U33     U12     U13     U23'
1022
1023        if General['Type'] == 'magnetic':
1024
1025            line += '   Mx     My     Mz'
1026
1027        elif General['Type'] == 'macromolecular':
1028
1029            line = ' res no  residue  chain '+line
1030
1031        print line
1032
1033        if General['Type'] == 'nuclear':
1034
1035            print 135*'-'
1036
1037            for i,at in enumerate(Atoms):
1038
1039                line = '%7s'%(at[0])+'%7s'%(at[1])+'%7s'%(at[2])+'%10.5f'%(at[3])+'%10.5f'%(at[4])+ \
1040
1041                    '%10.5f'%(at[5])+'%8.3f'%(at[6])+'%7s'%(at[7])+'%5d'%(at[8])+'%5s'%(at[9])
1042
1043                if at[9] == 'I':
1044
1045                    line += '%8.4f'%(at[10])+48*' '
1046
1047                else:
1048
1049                    line += 8*' '
1050
1051                    for j in range(6):
1052
1053                        line += '%8.4f'%(at[11+j])
1054
1055                print line
1056
1057       
1058
1059    def PrintTexture(textureData):
1060
1061        topstr = '\n Spherical harmonics texture: Order:' + \
1062
1063            str(textureData['Order'])
1064
1065        if textureData['Order']:
1066
1067            print topstr+' Refine? '+str(textureData['SH Coeff'][0])
1068
1069        else:
1070
1071            print topstr
1072
1073            return
1074
1075        names = ['omega','chi','phi']
1076
1077        line = '\n'
1078
1079        for name in names:
1080
1081            line += ' SH '+name+':'+'%12.4f'%(textureData['Sample '+name][1])+' Refine? '+str(textureData['Sample '+name][0])
1082
1083        print line
1084
1085        print '\n Texture coefficients:'
1086
1087        ptlbls = ' names :'
1088
1089        ptstr =  ' values:'
1090
1091        SHcoeff = textureData['SH Coeff'][1]
1092
1093        for item in SHcoeff:
1094
1095            ptlbls += '%12s'%(item)
1096
1097            ptstr += '%12.4f'%(SHcoeff[item]) 
1098
1099        print ptlbls
1100
1101        print ptstr   
1102
1103       
1104
1105    if Print: print ' Phases:'
1106
1107    phaseVary = []
1108
1109    phaseDict = {}
1110
1111    phaseConstr = {}
1112
1113    pawleyLookup = {}
1114
1115    FFtables = {}                   #scattering factors - xrays
1116
1117    BLtables = {}                   # neutrons
1118
1119    Natoms = {}
1120
1121    AtMults = {}
1122
1123    AtIA = {}
1124
1125    shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
1126
1127    SamSym = dict(zip(shModels,['0','-1','2/m','mmm']))
1128
1129    for name in PhaseData:
1130
1131        General = PhaseData[name]['General']
1132
1133        pId = PhaseData[name]['pId']
1134
1135        pfx = str(pId)+'::'
1136
1137        FFtable = GetFFtable(General)
1138
1139        BLtable = GetBLtable(General)
1140
1141        FFtables.update(FFtable)
1142
1143        BLtables.update(BLtable)
1144
1145        Atoms = PhaseData[name]['Atoms']
1146
1147        try:
1148
1149            PawleyRef = PhaseData[name]['Pawley ref']
1150
1151        except KeyError:
1152
1153            PawleyRef = []
1154
1155        SGData = General['SGData']
1156
1157        SGtext = G2spc.SGPrint(SGData)
1158
1159        cell = General['Cell']
1160
1161        A = G2lat.cell2A(cell[1:7])
1162
1163        phaseDict.update({pfx+'A0':A[0],pfx+'A1':A[1],pfx+'A2':A[2],pfx+'A3':A[3],pfx+'A4':A[4],pfx+'A5':A[5]})
1164
1165        if cell[0]:
1166
1167            phaseVary += cellVary(pfx,SGData)
1168
1169        Natoms[pfx] = 0
1170
1171        if Atoms:
1172
1173            if General['Type'] == 'nuclear':
1174
1175                Natoms[pfx] = len(Atoms)
1176
1177                for i,at in enumerate(Atoms):
1178
1179                    phaseDict.update({pfx+'Atype:'+str(i):at[1],pfx+'Afrac:'+str(i):at[6],pfx+'Amul:'+str(i):at[8],
1180
1181                        pfx+'Ax:'+str(i):at[3],pfx+'Ay:'+str(i):at[4],pfx+'Az:'+str(i):at[5],
1182
1183                        pfx+'dAx:'+str(i):0.,pfx+'dAy:'+str(i):0.,pfx+'dAz:'+str(i):0.,         #refined shifts for x,y,z
1184
1185                        pfx+'AI/A:'+str(i):at[9],})
1186
1187                    if at[9] == 'I':
1188
1189                        phaseDict[pfx+'AUiso:'+str(i)] = at[10]
1190
1191                    else:
1192
1193                        phaseDict.update({pfx+'AU11:'+str(i):at[11],pfx+'AU22:'+str(i):at[12],pfx+'AU33:'+str(i):at[13],
1194
1195                            pfx+'AU12:'+str(i):at[14],pfx+'AU13:'+str(i):at[15],pfx+'AU23:'+str(i):at[16]})
1196
1197                    if 'F' in at[2]:
1198
1199                        phaseVary.append(pfx+'Afrac:'+str(i))
1200
1201                    if 'X' in at[2]:
1202
1203                        xId,xCoef = G2spc.GetCSxinel(at[7])
1204
1205                        delnames = [pfx+'dAx:'+str(i),pfx+'dAy:'+str(i),pfx+'dAz:'+str(i)]
1206
1207                        for j in range(3):
1208
1209                            if xId[j] > 0:                               
1210
1211                                phaseVary.append(delnames[j])
1212
1213                                for k in range(j):
1214
1215                                    if xId[j] == xId[k]:
1216
1217                                        G2mv.StoreEquivalence(delnames[k],((delnames[j],xCoef[j]),)) 
1218
1219                    if 'U' in at[2]:
1220
1221                        if at[9] == 'I':
1222
1223                            phaseVary.append(pfx+'AUiso:'+str(i))
1224
1225                        else:
1226
1227                            uId,uCoef = G2spc.GetCSuinel(at[7])[:2]
1228
1229                            names = [pfx+'AU11:'+str(i),pfx+'AU22:'+str(i),pfx+'AU33:'+str(i),
1230
1231                                pfx+'AU12:'+str(i),pfx+'AU13:'+str(i),pfx+'AU23:'+str(i)]
1232
1233                            for j in range(6):
1234
1235                                if uId[j] > 0:                               
1236
1237                                    phaseVary.append(names[j])
1238
1239                                    for k in range(j):
1240
1241                                        if uId[j] == uId[k]:
1242
1243                                            G2mv.StoreEquivalence(names[k],((names[j],uCoef[j]),))
1244
1245#            elif General['Type'] == 'magnetic':
1246
1247#            elif General['Type'] == 'macromolecular':
1248
1249
1250
1251               
1252
1253            if 'SH Texture' in General:
1254
1255                textureData = General['SH Texture']
1256
1257                phaseDict[pfx+'SHmodel'] = SamSym[textureData['Model']]
1258
1259                phaseDict[pfx+'SHorder'] = textureData['Order']
1260
1261                for name in ['omega','chi','phi']:
1262
1263                    phaseDict[pfx+'SH '+name] = textureData['Sample '+name][1]
1264
1265                    if textureData['Sample '+name][0]:
1266
1267                        phaseVary.append(pfx+'SH '+name)
1268
1269                for name in textureData['SH Coeff'][1]:
1270
1271                    phaseDict[pfx+name] = textureData['SH Coeff'][1][name]
1272
1273                    if textureData['SH Coeff'][0]:
1274
1275                        phaseVary.append(pfx+name)
1276
1277               
1278
1279            if Print:
1280
1281                print '\n Phase name: ',General['Name']
1282
1283                print 135*'-'
1284
1285                PrintFFtable(FFtable)
1286
1287                PrintBLtable(BLtable)
1288
1289                print ''
1290
1291                for line in SGtext: print line
1292
1293                PrintAtoms(General,Atoms)
1294
1295                print '\n Unit cell: a =','%.5f'%(cell[1]),' b =','%.5f'%(cell[2]),' c =','%.5f'%(cell[3]), \
1296
1297                    ' alpha =','%.3f'%(cell[4]),' beta =','%.3f'%(cell[5]),' gamma =', \
1298
1299                    '%.3f'%(cell[6]),' volume =','%.3f'%(cell[7]),' Refine?',cell[0]
1300
1301                if 'SH Texture' in General:
1302
1303                    PrintTexture(textureData)
1304
1305                   
1306
1307        elif PawleyRef:
1308
1309            pawleyVary = []
1310
1311            for i,refl in enumerate(PawleyRef):
1312
1313                phaseDict[pfx+'PWLref:'+str(i)] = refl[6]
1314
1315                pawleyLookup[pfx+'%d,%d,%d'%(refl[0],refl[1],refl[2])] = i
1316
1317                if refl[5]:
1318
1319                    pawleyVary.append(pfx+'PWLref:'+str(i))
1320
1321            GetPawleyConstr(SGData['SGLaue'],PawleyRef,pawleyVary)      #does G2mv.StoreEquivalence
1322
1323            phaseVary += pawleyVary
1324
1325               
1326
1327    return Natoms,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables
1328
1329   
1330
1331def cellFill(pfx,SGData,parmDict,sigDict): 
1332
1333    if SGData['SGLaue'] in ['-1',]:
1334
1335        A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
1336
1337            parmDict[pfx+'A3'],parmDict[pfx+'A4'],parmDict[pfx+'A5']]
1338
1339        sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
1340
1341            sigDict[pfx+'A3'],sigDict[pfx+'A4'],sigDict[pfx+'A5']]
1342
1343    elif SGData['SGLaue'] in ['2/m',]:
1344
1345        if SGData['SGUniq'] == 'a':
1346
1347            A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
1348
1349                parmDict[pfx+'A3'],0,0]
1350
1351            sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
1352
1353                sigDict[pfx+'A3'],0,0]
1354
1355        elif SGData['SGUniq'] == 'b':
1356
1357            A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
1358
1359                0,parmDict[pfx+'A4'],0]
1360
1361            sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
1362
1363                0,sigDict[pfx+'A4'],0]
1364
1365        else:
1366
1367            A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
1368
1369                0,0,parmDict[pfx+'A5']]
1370
1371            sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
1372
1373                0,0,sigDict[pfx+'A5']]
1374
1375    elif SGData['SGLaue'] in ['mmm',]:
1376
1377        A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],0,0,0]
1378
1379        sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],0,0,0]
1380
1381    elif SGData['SGLaue'] in ['4/m','4/mmm']:
1382
1383        A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A2'],0,0,0]
1384
1385        sigA = [sigDict[pfx+'A0'],0,sigDict[pfx+'A2'],0,0,0]
1386
1387    elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
1388
1389        A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A2'],
1390
1391            parmDict[pfx+'A0'],0,0]
1392
1393        sigA = [sigDict[pfx+'A0'],0,sigDict[pfx+'A2'],0,0,0]
1394
1395    elif SGData['SGLaue'] in ['3R', '3mR']:
1396
1397        A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A0'],
1398
1399            parmDict[pfx+'A3'],parmDict[pfx+'A3'],parmDict[pfx+'A3']]
1400
1401        sigA = [sigDict[pfx+'A0'],0,0,sigDict[pfx+'A3'],0,0]
1402
1403    elif SGData['SGLaue'] in ['m3m','m3']:
1404
1405        A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A0'],0,0,0]
1406
1407        sigA = [sigDict[pfx+'A0'],0,0,0,0,0]
1408
1409    return A,sigA
1410
1411       
1412
1413def getCellEsd(pfx,SGData,A,covData):
1414
1415    dpr = 180./np.pi
1416
1417    rVsq = G2lat.calc_rVsq(A)
1418
1419    G,g = G2lat.A2Gmat(A)       #get recip. & real metric tensors
1420
1421    cell = np.array(G2lat.Gmat2cell(g))   #real cell
1422
1423    cellst = np.array(G2lat.Gmat2cell(G)) #recip. cell
1424
1425    scos = cosd(cellst[3:6])
1426
1427    ssin = sind(cellst[3:6])
1428
1429    scot = scos/ssin
1430
1431    rcos = cosd(cell[3:6])
1432
1433    rsin = sind(cell[3:6])
1434
1435    rcot = rcos/rsin
1436
1437    RMnames = [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3',pfx+'A4',pfx+'A5']
1438
1439    varyList = covData['varyList']
1440
1441    covMatrix = covData['covMatrix']
1442
1443    vcov = G2mth.getVCov(RMnames,varyList,covMatrix)
1444
1445    Ax = np.array(A)
1446
1447    Ax[3:] /= 2.
1448
1449    drVdA = np.array([Ax[1]*Ax[2]-Ax[5]**2,Ax[0]*Ax[2]-Ax[4]**2,Ax[0]*Ax[1]-Ax[3]**2,
1450
1451        Ax[4]*Ax[5]-Ax[2]*Ax[3],Ax[3]*Ax[5]-Ax[1]*Ax[4],Ax[3]*Ax[4]-Ax[0]*Ax[5]])
1452
1453    srcvlsq = np.inner(drVdA,np.inner(vcov,drVdA.T))
1454
1455    Vol = 1/np.sqrt(rVsq)
1456
1457    sigVol = Vol**3*np.sqrt(srcvlsq)/2.
1458
1459    R123 = Ax[0]*Ax[1]*Ax[2]
1460
1461    dsasdg = np.zeros((3,6))
1462
1463    dadg = np.zeros((6,6))
1464
1465    for i0 in range(3):         #0  1   2
1466
1467        i1 = (i0+1)%3           #1  2   0
1468
1469        i2 = (i1+1)%3           #2  0   1
1470
1471        i3 = 5-i2               #3  5   4
1472
1473        i4 = 5-i1               #4  3   5
1474
1475        i5 = 5-i0               #5  4   3
1476
1477        dsasdg[i0][i1] = 0.5*scot[i0]*scos[i0]/Ax[i1]
1478
1479        dsasdg[i0][i2] = 0.5*scot[i0]*scos[i0]/Ax[i2]
1480
1481        dsasdg[i0][i5] = -scot[i0]/np.sqrt(Ax[i1]*Ax[i2])
1482
1483        denmsq = Ax[i0]*(R123-Ax[i1]*Ax[i4]**2-Ax[i2]*Ax[i3]**2+(Ax[i4]*Ax[i3])**2)
1484
1485        denom = np.sqrt(denmsq)
1486
1487        dadg[i5][i0] = -Ax[i5]/denom-rcos[i0]/denmsq*(R123-0.5*Ax[i1]*Ax[i4]**2-0.5*Ax[i2]*Ax[i3]**2)
1488
1489        dadg[i5][i1] = -0.5*rcos[i0]/denmsq*(Ax[i0]**2*Ax[i2]-Ax[i0]*Ax[i4]**2)
1490
1491        dadg[i5][i2] = -0.5*rcos[i0]/denmsq*(Ax[i0]**2*Ax[i1]-Ax[i0]*Ax[i3]**2)
1492
1493        dadg[i5][i3] = Ax[i4]/denom+rcos[i0]/denmsq*(Ax[i0]*Ax[i2]*Ax[i3]-Ax[i3]*Ax[i4]**2)
1494
1495        dadg[i5][i4] = Ax[i3]/denom+rcos[i0]/denmsq*(Ax[i0]*Ax[i1]*Ax[i4]-Ax[i3]**2*Ax[i4])
1496
1497        dadg[i5][i5] = -Ax[i0]/denom
1498
1499    for i0 in range(3):
1500
1501        i1 = (i0+1)%3
1502
1503        i2 = (i1+1)%3
1504
1505        i3 = 5-i2
1506
1507        for ij in range(6):
1508
1509            dadg[i0][ij] = cell[i0]*(rcot[i2]*dadg[i3][ij]/rsin[i2]-dsasdg[i1][ij]/ssin[i1])
1510
1511            if ij == i0:
1512
1513                dadg[i0][ij] = dadg[i0][ij]-0.5*cell[i0]/Ax[i0]
1514
1515            dadg[i3][ij] = -dadg[i3][ij]*rsin[2-i0]*dpr
1516
1517    sigMat = np.inner(dadg,np.inner(vcov,dadg.T))
1518
1519    var = np.diag(sigMat)
1520
1521    CS = np.where(var>0.,np.sqrt(var),0.)
1522
1523    cellSig = [CS[0],CS[1],CS[2],CS[5],CS[4],CS[3],sigVol]  #exchange sig(alp) & sig(gam) to get in right order
1524
1525    return cellSig           
1526
1527   
1528
1529def SetPhaseData(parmDict,sigDict,Phases,covData):
1530
1531   
1532
1533    def PrintAtomsAndSig(General,Atoms,atomsSig):
1534
1535        print '\n Atoms:'
1536
1537        line = '   name      x         y         z      frac   Uiso     U11     U22     U33     U12     U13     U23'
1538
1539        if General['Type'] == 'magnetic':
1540
1541            line += '   Mx     My     Mz'
1542
1543        elif General['Type'] == 'macromolecular':
1544
1545            line = ' res no  residue  chain '+line
1546
1547        print line
1548
1549        if General['Type'] == 'nuclear':
1550
1551            print 135*'-'
1552
1553            fmt = {0:'%7s',1:'%7s',3:'%10.5f',4:'%10.5f',5:'%10.5f',6:'%8.3f',10:'%8.5f',
1554
1555                11:'%8.5f',12:'%8.5f',13:'%8.5f',14:'%8.5f',15:'%8.5f',16:'%8.5f'}
1556
1557            noFXsig = {3:[10*' ','%10s'],4:[10*' ','%10s'],5:[10*' ','%10s'],6:[8*' ','%8s']}
1558
1559            for i,at in enumerate(Atoms):
1560
1561                name = fmt[0]%(at[0])+fmt[1]%(at[1])+':'
1562
1563                valstr = ' values:'
1564
1565                sigstr = ' sig   :'
1566
1567                for ind in [3,4,5,6]:
1568
1569                    sigind = str(i)+':'+str(ind)
1570
1571                    valstr += fmt[ind]%(at[ind])                   
1572
1573                    if sigind in atomsSig:
1574
1575                        sigstr += fmt[ind]%(atomsSig[sigind])
1576
1577                    else:
1578
1579                        sigstr += noFXsig[ind][1]%(noFXsig[ind][0])
1580
1581                if at[9] == 'I':
1582
1583                    valstr += fmt[10]%(at[10])
1584
1585                    if str(i)+':10' in atomsSig:
1586
1587                        sigstr += fmt[10]%(atomsSig[str(i)+':10'])
1588
1589                    else:
1590
1591                        sigstr += 8*' '
1592
1593                else:
1594
1595                    valstr += 8*' '
1596
1597                    sigstr += 8*' '
1598
1599                    for ind in [11,12,13,14,15,16]:
1600
1601                        sigind = str(i)+':'+str(ind)
1602
1603                        valstr += fmt[ind]%(at[ind])
1604
1605                        if sigind in atomsSig:                       
1606
1607                            sigstr += fmt[ind]%(atomsSig[sigind])
1608
1609                        else:
1610
1611                            sigstr += 8*' '
1612
1613                print name
1614
1615                print valstr
1616
1617                print sigstr
1618
1619               
1620
1621    def PrintSHtextureAndSig(textureData,SHtextureSig):
1622
1623        print '\n Spherical harmonics texture: Order:' + str(textureData['Order'])
1624
1625        names = ['omega','chi','phi']
1626
1627        namstr = '  names :'
1628
1629        ptstr =  '  values:'
1630
1631        sigstr = '  esds  :'
1632
1633        for name in names:
1634
1635            namstr += '%12s'%(name)
1636
1637            ptstr += '%12.3f'%(textureData['Sample '+name][1])
1638
1639            if 'Sample '+name in SHtextureSig:
1640
1641                sigstr += '%12.3f'%(SHtextureSig['Sample '+name])
1642
1643            else:
1644
1645                sigstr += 12*' '
1646
1647        print namstr
1648
1649        print ptstr
1650
1651        print sigstr
1652
1653        print '\n Texture coefficients:'
1654
1655        namstr = '  names :'
1656
1657        ptstr =  '  values:'
1658
1659        sigstr = '  esds  :'
1660
1661        SHcoeff = textureData['SH Coeff'][1]
1662
1663        for name in SHcoeff:
1664
1665            namstr += '%12s'%(name)
1666
1667            ptstr += '%12.3f'%(SHcoeff[name])
1668
1669            if name in SHtextureSig:
1670
1671                sigstr += '%12.3f'%(SHtextureSig[name])
1672
1673            else:
1674
1675                sigstr += 12*' '
1676
1677        print namstr
1678
1679        print ptstr
1680
1681        print sigstr
1682
1683       
1684
1685           
1686
1687    print '\n Phases:'
1688
1689    for phase in Phases:
1690
1691        print ' Result for phase: ',phase
1692
1693        Phase = Phases[phase]
1694
1695        General = Phase['General']
1696
1697        SGData = General['SGData']
1698
1699        Atoms = Phase['Atoms']
1700
1701        cell = General['Cell']
1702
1703        pId = Phase['pId']
1704
1705        pfx = str(pId)+'::'
1706
1707        if cell[0]:
1708
1709            A,sigA = cellFill(pfx,SGData,parmDict,sigDict)
1710
1711            cellSig = getCellEsd(pfx,SGData,A,covData)  #includes sigVol
1712
1713            print ' Reciprocal metric tensor: '
1714
1715            ptfmt = "%15.9f"
1716
1717            names = ['A11','A22','A33','A12','A13','A23']
1718
1719            namstr = '  names :'
1720
1721            ptstr =  '  values:'
1722
1723            sigstr = '  esds  :'
1724
1725            for name,a,siga in zip(names,A,sigA):
1726
1727                namstr += '%15s'%(name)
1728
1729                ptstr += ptfmt%(a)
1730
1731                if siga:
1732
1733                    sigstr += ptfmt%(siga)
1734
1735                else:
1736
1737                    sigstr += 15*' '
1738
1739            print namstr
1740
1741            print ptstr
1742
1743            print sigstr
1744
1745            cell[1:7] = G2lat.A2cell(A)
1746
1747            cell[7] = G2lat.calc_V(A)
1748
1749            print ' New unit cell:'
1750
1751            ptfmt = ["%12.6f","%12.6f","%12.6f","%12.4f","%12.4f","%12.4f","%12.3f"]
1752
1753            names = ['a','b','c','alpha','beta','gamma','Volume']
1754
1755            namstr = '  names :'
1756
1757            ptstr =  '  values:'
1758
1759            sigstr = '  esds  :'
1760
1761            for name,fmt,a,siga in zip(names,ptfmt,cell[1:8],cellSig):
1762
1763                namstr += '%12s'%(name)
1764
1765                ptstr += fmt%(a)
1766
1767                if siga:
1768
1769                    sigstr += fmt%(siga)
1770
1771                else:
1772
1773                    sigstr += 12*' '
1774
1775            print namstr
1776
1777            print ptstr
1778
1779            print sigstr
1780
1781           
1782
1783        if 'Pawley' in Phase['General']['Type']:
1784
1785            pawleyRef = Phase['Pawley ref']
1786
1787            for i,refl in enumerate(pawleyRef):
1788
1789                key = pfx+'PWLref:'+str(i)
1790
1791                refl[6] = abs(parmDict[key])        #suppress negative Fsq
1792
1793                if key in sigDict:
1794
1795                    refl[7] = sigDict[key]
1796
1797                else:
1798
1799                    refl[7] = 0
1800
1801        else:
1802
1803            atomsSig = {}
1804
1805            if General['Type'] == 'nuclear':
1806
1807                for i,at in enumerate(Atoms):
1808
1809                    names = {3:pfx+'Ax:'+str(i),4:pfx+'Ay:'+str(i),5:pfx+'Az:'+str(i),6:pfx+'Afrac:'+str(i),
1810
1811                        10:pfx+'AUiso:'+str(i),11:pfx+'AU11:'+str(i),12:pfx+'AU22:'+str(i),13:pfx+'AU33:'+str(i),
1812
1813                        14:pfx+'AU12:'+str(i),15:pfx+'AU13:'+str(i),16:pfx+'AU23:'+str(i)}
1814
1815                    for ind in [3,4,5,6]:
1816
1817                        at[ind] = parmDict[names[ind]]
1818
1819                        if ind in [3,4,5]:
1820
1821                            name = names[ind].replace('A','dA')
1822
1823                        else:
1824
1825                            name = names[ind]
1826
1827                        if name in sigDict:
1828
1829                            atomsSig[str(i)+':'+str(ind)] = sigDict[name]
1830
1831                    if at[9] == 'I':
1832
1833                        at[10] = parmDict[names[10]]
1834
1835                        if names[10] in sigDict:
1836
1837                            atomsSig[str(i)+':10'] = sigDict[names[10]]
1838
1839                    else:
1840
1841                        for ind in [11,12,13,14,15,16]:
1842
1843                            at[ind] = parmDict[names[ind]]
1844
1845                            if names[ind] in sigDict:
1846
1847                                atomsSig[str(i)+':'+str(ind)] = sigDict[names[ind]]
1848
1849            PrintAtomsAndSig(General,Atoms,atomsSig)
1850
1851       
1852
1853        if 'SH Texture' in General:
1854
1855            textureData = General['SH Texture']   
1856
1857            if textureData['Order']:
1858
1859                SHtextureSig = {}
1860
1861                for name in ['omega','chi','phi']:
1862
1863                    aname = pfx+'SH '+name
1864
1865                    textureData['Sample '+name][1] = parmDict[aname]
1866
1867                    if aname in sigDict:
1868
1869                        SHtextureSig['Sample '+name] = sigDict[aname]
1870
1871                for name in textureData['SH Coeff'][1]:
1872
1873                    aname = pfx+name
1874
1875                    textureData['SH Coeff'][1][name] = parmDict[aname]
1876
1877                    if aname in sigDict:
1878
1879                        SHtextureSig[name] = sigDict[aname]
1880
1881                PrintSHtextureAndSig(textureData,SHtextureSig)
1882
1883
1884
1885def GetHistogramPhaseData(Phases,Histograms,Print=True):
1886
1887   
1888
1889    def PrintSize(hapData):
1890
1891        if hapData[0] in ['isotropic','uniaxial']:
1892
1893            line = '\n Size model    : %9s'%(hapData[0])
1894
1895            line += ' equatorial:'+'%12.3f'%(hapData[1][0])+' Refine? '+str(hapData[2][0])
1896
1897            if hapData[0] == 'uniaxial':
1898
1899                line += ' axial:'+'%12.3f'%(hapData[1][1])+' Refine? '+str(hapData[2][1])
1900
1901            print line
1902
1903        else:
1904
1905            print '\n Size model    : %s'%(hapData[0])
1906
1907            Snames = ['S11','S22','S33','S12','S13','S23']
1908
1909            ptlbls = ' names :'
1910
1911            ptstr =  ' values:'
1912
1913            varstr = ' refine:'
1914
1915            for i,name in enumerate(Snames):
1916
1917                ptlbls += '%12s' % (name)
1918
1919                ptstr += '%12.6f' % (hapData[4][i])
1920
1921                varstr += '%12s' % (str(hapData[5][i]))
1922
1923            print ptlbls
1924
1925            print ptstr
1926
1927            print varstr
1928
1929       
1930
1931    def PrintMuStrain(hapData,SGData):
1932
1933        if hapData[0] in ['isotropic','uniaxial']:
1934
1935            line = '\n Mustrain model: %9s'%(hapData[0])
1936
1937            line += ' equatorial:'+'%12.1f'%(hapData[1][0])+' Refine? '+str(hapData[2][0])
1938
1939            if hapData[0] == 'uniaxial':
1940
1941                line += ' axial:'+'%12.1f'%(hapData[1][1])+' Refine? '+str(hapData[2][1])
1942
1943            print line
1944
1945        else:
1946
1947            print '\n Mustrain model: %s'%(hapData[0])
1948
1949            Snames = G2spc.MustrainNames(SGData)
1950
1951            ptlbls = ' names :'
1952
1953            ptstr =  ' values:'
1954
1955            varstr = ' refine:'
1956
1957            for i,name in enumerate(Snames):
1958
1959                ptlbls += '%12s' % (name)
1960
1961                ptstr += '%12.6f' % (hapData[4][i])
1962
1963                varstr += '%12s' % (str(hapData[5][i]))
1964
1965            print ptlbls
1966
1967            print ptstr
1968
1969            print varstr
1970
1971
1972
1973    def PrintHStrain(hapData,SGData):
1974
1975        print '\n Hydrostatic/elastic strain: '
1976
1977        Hsnames = G2spc.HStrainNames(SGData)
1978
1979        ptlbls = ' names :'
1980
1981        ptstr =  ' values:'
1982
1983        varstr = ' refine:'
1984
1985        for i,name in enumerate(Hsnames):
1986
1987            ptlbls += '%12s' % (name)
1988
1989            ptstr += '%12.6f' % (hapData[0][i])
1990
1991            varstr += '%12s' % (str(hapData[1][i]))
1992
1993        print ptlbls
1994
1995        print ptstr
1996
1997        print varstr
1998
1999
2000
2001    def PrintSHPO(hapData):
2002
2003        print '\n Spherical harmonics preferred orientation: Order:' + \
2004
2005            str(hapData[4])+' Refine? '+str(hapData[2])
2006
2007        ptlbls = ' names :'
2008
2009        ptstr =  ' values:'
2010
2011        for item in hapData[5]:
2012
2013            ptlbls += '%12s'%(item)
2014
2015            ptstr += '%12.3f'%(hapData[5][item]) 
2016
2017        print ptlbls
2018
2019        print ptstr
2020
2021   
2022
2023    hapDict = {}
2024
2025    hapVary = []
2026
2027    controlDict = {}
2028
2029    poType = {}
2030
2031    poAxes = {}
2032
2033    spAxes = {}
2034
2035    spType = {}
2036
2037   
2038
2039    for phase in Phases:
2040
2041        HistoPhase = Phases[phase]['Histograms']
2042
2043        SGData = Phases[phase]['General']['SGData']
2044
2045        cell = Phases[phase]['General']['Cell'][1:7]
2046
2047        A = G2lat.cell2A(cell)
2048
2049        pId = Phases[phase]['pId']
2050
2051        histoList = HistoPhase.keys()
2052
2053        histoList.sort()
2054
2055        for histogram in histoList:
2056
2057            try:
2058
2059                Histogram = Histograms[histogram]
2060
2061            except KeyError:                       
2062
2063                #skip if histogram not included e.g. in a sequential refinement
2064
2065                continue
2066
2067            hapData = HistoPhase[histogram]
2068
2069            hId = Histogram['hId']
2070
2071            limits = Histogram['Limits'][1]
2072
2073            inst = Histogram['Instrument Parameters']
2074
2075            inst = dict(zip(inst[3],inst[1]))
2076
2077            Zero = inst['Zero']
2078
2079            if 'C' in inst['Type']:
2080
2081                try:
2082
2083                    wave = inst['Lam']
2084
2085                except KeyError:
2086
2087                    wave = inst['Lam1']
2088
2089                dmin = wave/(2.0*sind(limits[1]/2.0))
2090
2091            pfx = str(pId)+':'+str(hId)+':'
2092
2093            for item in ['Scale','Extinction']:
2094
2095                hapDict[pfx+item] = hapData[item][0]
2096
2097                if hapData[item][1]:
2098
2099                    hapVary.append(pfx+item)
2100
2101            names = G2spc.HStrainNames(SGData)
2102
2103            for i,name in enumerate(names):
2104
2105                hapDict[pfx+name] = hapData['HStrain'][0][i]
2106
2107                if hapData['HStrain'][1][i]:
2108
2109                    hapVary.append(pfx+name)
2110
2111            controlDict[pfx+'poType'] = hapData['Pref.Ori.'][0]
2112
2113            if hapData['Pref.Ori.'][0] == 'MD':
2114
2115                hapDict[pfx+'MD'] = hapData['Pref.Ori.'][1]
2116
2117                controlDict[pfx+'MDAxis'] = hapData['Pref.Ori.'][3]
2118
2119                if hapData['Pref.Ori.'][2]:
2120
2121                    hapVary.append(pfx+'MD')
2122
2123            else:                           #'SH' spherical harmonics
2124
2125                controlDict[pfx+'SHord'] = hapData['Pref.Ori.'][4]
2126
2127                controlDict[pfx+'SHncof'] = len(hapData['Pref.Ori.'][5])
2128
2129                for item in hapData['Pref.Ori.'][5]:
2130
2131                    hapDict[pfx+item] = hapData['Pref.Ori.'][5][item]
2132
2133                    if hapData['Pref.Ori.'][2]:
2134
2135                        hapVary.append(pfx+item)
2136
2137            for item in ['Mustrain','Size']:
2138
2139                controlDict[pfx+item+'Type'] = hapData[item][0]
2140
2141                if hapData[item][0] in ['isotropic','uniaxial']:
2142
2143                    hapDict[pfx+item+':i'] = hapData[item][1][0]
2144
2145                    if hapData[item][2][0]:
2146
2147                        hapVary.append(pfx+item+':i')
2148
2149                    if hapData[item][0] == 'uniaxial':
2150
2151                        controlDict[pfx+item+'Axis'] = hapData[item][3]
2152
2153                        hapDict[pfx+item+':a'] = hapData[item][1][1]
2154
2155                        if hapData[item][2][1]:
2156
2157                            hapVary.append(pfx+item+':a')
2158
2159                else:       #generalized for mustrain or ellipsoidal for size
2160
2161                    if item == 'Mustrain':
2162
2163                        names = G2spc.MustrainNames(SGData)
2164
2165                        pwrs = []
2166
2167                        for name in names:
2168
2169                            h,k,l = name[1:]
2170
2171                            pwrs.append([int(h),int(k),int(l)])
2172
2173                        controlDict[pfx+'MuPwrs'] = pwrs
2174
2175                    for i in range(len(hapData[item][4])):
2176
2177                        sfx = ':'+str(i)
2178
2179                        hapDict[pfx+item+sfx] = hapData[item][4][i]
2180
2181                        if hapData[item][5][i]:
2182
2183                            hapVary.append(pfx+item+sfx)
2184
2185                           
2186
2187            if Print: 
2188
2189                print '\n Phase: ',phase,' in histogram: ',histogram
2190
2191                print 135*'-'
2192
2193                print ' Phase fraction  : %10.4f'%(hapData['Scale'][0]),' Refine?',hapData['Scale'][1]
2194
2195                print ' Extinction coeff: %10.4f'%(hapData['Extinction'][0]),' Refine?',hapData['Extinction'][1]
2196
2197                if hapData['Pref.Ori.'][0] == 'MD':
2198
2199                    Ax = hapData['Pref.Ori.'][3]
2200
2201                    print ' March-Dollase PO: %10.4f'%(hapData['Pref.Ori.'][1]),' Refine?',hapData['Pref.Ori.'][2], \
2202
2203                        ' Axis: %d %d %d'%(Ax[0],Ax[1],Ax[2])
2204
2205                else: #'SH' for spherical harmonics
2206
2207                    PrintSHPO(hapData['Pref.Ori.'])
2208
2209                PrintSize(hapData['Size'])
2210
2211                PrintMuStrain(hapData['Mustrain'],SGData)
2212
2213                PrintHStrain(hapData['HStrain'],SGData)
2214
2215            HKLd = np.array(G2lat.GenHLaue(dmin,SGData,A))
2216
2217            refList = []
2218
2219            for h,k,l,d in HKLd:
2220
2221                ext,mul,Uniq,phi = G2spc.GenHKLf([h,k,l],SGData)
2222
2223                if ext:
2224
2225                    continue
2226
2227                if 'C' in inst['Type']:
2228
2229                    pos = 2.0*asind(wave/(2.0*d))+Zero
2230
2231                    if limits[0] < pos < limits[1]:
2232
2233                        refList.append([h,k,l,mul,d,pos,0.0,0.0,0.0,0.0,0.0,Uniq,phi,0.0])
2234
2235                else:
2236
2237                    raise ValueError 
2238
2239            Histogram['Reflection Lists'][phase] = refList
2240
2241    return hapVary,hapDict,controlDict
2242
2243   
2244
2245def SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms,Print=True):
2246
2247   
2248
2249    def PrintSizeAndSig(hapData,sizeSig):
2250
2251        line = '\n Size model:     %9s'%(hapData[0])
2252
2253        refine = False
2254
2255        if hapData[0] in ['isotropic','uniaxial']:
2256
2257            line += ' equatorial:%12.3f'%(hapData[1][0])
2258
2259            if sizeSig[0][0]:
2260
2261                line += ', sig: %8.3f'%(sizeSig[0][0])
2262
2263                refine = True
2264
2265            if hapData[0] == 'uniaxial':
2266
2267                line += ' axial:%12.3f'%(hapData[1][1])
2268
2269                if sizeSig[0][1]:
2270
2271                    refine = True
2272
2273                    line += ', sig: %8.3f'%(sizeSig[0][1])
2274
2275            if refine:
2276
2277                print line
2278
2279        else:
2280
2281            Snames = ['S11','S22','S33','S12','S13','S23']
2282
2283            ptlbls = ' name  :'
2284
2285            ptstr =  ' value :'
2286
2287            sigstr = ' sig   :'
2288
2289            for i,name in enumerate(Snames):
2290
2291                ptlbls += '%12s' % (name)
2292
2293                ptstr += '%12.6f' % (hapData[4][i])
2294
2295                if sizeSig[1][i]:
2296
2297                    refine = True
2298
2299                    sigstr += '%12.6f' % (sizeSig[1][i])
2300
2301                else:
2302
2303                    sigstr += 12*' '
2304
2305            if refine:
2306
2307                print line
2308
2309                print ptlbls
2310
2311                print ptstr
2312
2313                print sigstr
2314
2315       
2316
2317    def PrintMuStrainAndSig(hapData,mustrainSig,SGData):
2318
2319        line = '\n Mustrain model: %9s'%(hapData[0])
2320
2321        refine = False
2322
2323        if hapData[0] in ['isotropic','uniaxial']:
2324
2325            line += ' equatorial:%12.1f'%(hapData[1][0])
2326
2327            if mustrainSig[0][0]:
2328
2329                line += ', sig: %8.1f'%(mustrainSig[0][0])
2330
2331                refine = True
2332
2333            if hapData[0] == 'uniaxial':
2334
2335                line += ' axial:%12.1f'%(hapData[1][1])
2336
2337                if mustrainSig[0][1]:
2338
2339                     line += ', sig: %8.1f'%(mustrainSig[0][1])
2340
2341            if refine:
2342
2343                print line
2344
2345        else:
2346
2347            Snames = G2spc.MustrainNames(SGData)
2348
2349            ptlbls = ' name  :'
2350
2351            ptstr =  ' value :'
2352
2353            sigstr = ' sig   :'
2354
2355            for i,name in enumerate(Snames):
2356
2357                ptlbls += '%12s' % (name)
2358
2359                ptstr += '%12.6f' % (hapData[4][i])
2360
2361                if mustrainSig[1][i]:
2362
2363                    refine = True
2364
2365                    sigstr += '%12.6f' % (mustrainSig[1][i])
2366
2367                else:
2368
2369                    sigstr += 12*' '
2370
2371            if refine:
2372
2373                print line
2374
2375                print ptlbls
2376
2377                print ptstr
2378
2379                print sigstr
2380
2381           
2382
2383    def PrintHStrainAndSig(hapData,strainSig,SGData):
2384
2385        Hsnames = G2spc.HStrainNames(SGData)
2386
2387        ptlbls = ' name  :'
2388
2389        ptstr =  ' value :'
2390
2391        sigstr = ' sig   :'
2392
2393        refine = False
2394
2395        for i,name in enumerate(Hsnames):
2396
2397            ptlbls += '%12s' % (name)
2398
2399            ptstr += '%12.6g' % (hapData[0][i])
2400
2401            if name in strainSig:
2402
2403                refine = True
2404
2405                sigstr += '%12.6g' % (strainSig[name])
2406
2407            else:
2408
2409                sigstr += 12*' '
2410
2411        if refine:
2412
2413            print '\n Hydrostatic/elastic strain: '
2414
2415            print ptlbls
2416
2417            print ptstr
2418
2419            print sigstr
2420
2421       
2422
2423    def PrintSHPOAndSig(hapData,POsig):
2424
2425        print '\n Spherical harmonics preferred orientation: Order:'+str(hapData[4])
2426
2427        ptlbls = ' names :'
2428
2429        ptstr =  ' values:'
2430
2431        sigstr = ' sig   :'
2432
2433        for item in hapData[5]:
2434
2435            ptlbls += '%12s'%(item)
2436
2437            ptstr += '%12.3f'%(hapData[5][item])
2438
2439            if item in POsig:
2440
2441                sigstr += '%12.3f'%(POsig[item])
2442
2443            else:
2444
2445                sigstr += 12*' ' 
2446
2447        print ptlbls
2448
2449        print ptstr
2450
2451        print sigstr
2452
2453   
2454
2455    for phase in Phases:
2456
2457        HistoPhase = Phases[phase]['Histograms']
2458
2459        SGData = Phases[phase]['General']['SGData']
2460
2461        pId = Phases[phase]['pId']
2462
2463        histoList = HistoPhase.keys()
2464
2465        histoList.sort()
2466
2467        for histogram in histoList:
2468
2469            try:
2470
2471                Histogram = Histograms[histogram]
2472
2473            except KeyError:                       
2474
2475                #skip if histogram not included e.g. in a sequential refinement
2476
2477                continue
2478
2479            print '\n Phase: ',phase,' in histogram: ',histogram
2480
2481            print 130*'-'
2482
2483            hapData = HistoPhase[histogram]
2484
2485            hId = Histogram['hId']
2486
2487            pfx = str(pId)+':'+str(hId)+':'
2488
2489            print ' Final refinement RF, RF^2 = %.2f%%, %.2f%% on %d reflections'   \
2490
2491                %(Histogram[pfx+'Rf'],Histogram[pfx+'Rf^2'],Histogram[pfx+'Nref'])
2492
2493           
2494
2495            PhFrExtPOSig = {}
2496
2497            for item in ['Scale','Extinction']:
2498
2499                hapData[item][0] = parmDict[pfx+item]
2500
2501                if pfx+item in sigDict:
2502
2503                    PhFrExtPOSig[item] = sigDict[pfx+item]
2504
2505            if hapData['Pref.Ori.'][0] == 'MD':
2506
2507                hapData['Pref.Ori.'][1] = parmDict[pfx+'MD']
2508
2509                if pfx+'MD' in sigDict:
2510
2511                    PhFrExtPOSig['MD'] = sigDict[pfx+'MD']
2512
2513            else:                           #'SH' spherical harmonics
2514
2515                for item in hapData['Pref.Ori.'][5]:
2516
2517                    hapData['Pref.Ori.'][5][item] = parmDict[pfx+item]
2518
2519                    if pfx+item in sigDict:
2520
2521                        PhFrExtPOSig[item] = sigDict[pfx+item]
2522
2523            if Print:
2524
2525                if 'Scale' in PhFrExtPOSig:
2526
2527                    print ' Phase fraction  : %10.4f, sig %10.4f'%(hapData['Scale'][0],PhFrExtPOSig['Scale'])
2528
2529                if 'Extinction' in PhFrExtPOSig:
2530
2531                    print ' Extinction coeff: %10.4f, sig %10.4f'%(hapData['Extinction'][0],PhFrExtPOSig['Extinction'])
2532
2533                if hapData['Pref.Ori.'][0] == 'MD':
2534
2535                    if 'MD' in PhFrExtPOSig:
2536
2537                        print ' March-Dollase PO: %10.4f, sig %10.4f'%(hapData['Pref.Ori.'][1],PhFrExtPOSig['MD'])
2538
2539                else:
2540
2541                    PrintSHPOAndSig(hapData['Pref.Ori.'],PhFrExtPOSig)
2542
2543            SizeMuStrSig = {'Mustrain':[[0,0],[0 for i in range(len(hapData['Mustrain'][4]))]],
2544
2545                'Size':[[0,0],[0 for i in range(len(hapData['Size'][4]))]],
2546
2547                'HStrain':{}}                 
2548
2549            for item in ['Mustrain','Size']:
2550
2551                if hapData[item][0] in ['isotropic','uniaxial']:                   
2552
2553                    hapData[item][1][0] = parmDict[pfx+item+':i']
2554
2555                    if item == 'Size':
2556
2557                        hapData[item][1][0] = min(10.,max(0.001,hapData[item][1][0]))
2558
2559                    if pfx+item+':i' in sigDict: 
2560
2561                        SizeMuStrSig[item][0][0] = sigDict[pfx+item+':i']
2562
2563                    if hapData[item][0] == 'uniaxial':
2564
2565                        hapData[item][1][1] = parmDict[pfx+item+':a']
2566
2567                        if item == 'Size':
2568
2569                            hapData[item][1][1] = min(10.,max(0.001,hapData[item][1][1]))                       
2570
2571                        if pfx+item+':a' in sigDict:
2572
2573                            SizeMuStrSig[item][0][1] = sigDict[pfx+item+':a']
2574
2575                else:       #generalized for mustrain or ellipsoidal for size
2576
2577                    for i in range(len(hapData[item][4])):
2578
2579                        sfx = ':'+str(i)
2580
2581                        hapData[item][4][i] = parmDict[pfx+item+sfx]
2582
2583                        if pfx+item+sfx in sigDict:
2584
2585                            SizeMuStrSig[item][1][i] = sigDict[pfx+item+sfx]
2586
2587            names = G2spc.HStrainNames(SGData)
2588
2589            for i,name in enumerate(names):
2590
2591                hapData['HStrain'][0][i] = parmDict[pfx+name]
2592
2593                if pfx+name in sigDict:
2594
2595                    SizeMuStrSig['HStrain'][name] = sigDict[pfx+name]
2596
2597            if Print:
2598
2599                PrintSizeAndSig(hapData['Size'],SizeMuStrSig['Size'])
2600
2601                PrintMuStrainAndSig(hapData['Mustrain'],SizeMuStrSig['Mustrain'],SGData)
2602
2603                PrintHStrainAndSig(hapData['HStrain'],SizeMuStrSig['HStrain'],SGData)
2604
2605   
2606
2607def GetHistogramData(Histograms,Print=True):
2608
2609   
2610
2611    def GetBackgroundParms(hId,Background):
2612
2613        Back = Background[0]
2614
2615        Debye = Background[1]
2616
2617        bakType,bakFlag = Back[:2]
2618
2619        backVals = Back[3:]
2620
2621        backNames = [':'+str(hId)+':Back:'+str(i) for i in range(len(backVals))]
2622
2623        backDict = dict(zip(backNames,backVals))
2624
2625        backVary = []
2626
2627        if bakFlag:
2628
2629            backVary = backNames
2630
2631        backDict[':'+str(hId)+':nDebye'] = Debye['nDebye']
2632
2633        debyeDict = {}
2634
2635        debyeList = []
2636
2637        for i in range(Debye['nDebye']):
2638
2639            debyeNames = [':'+str(hId)+':DebyeA:'+str(i),':'+str(hId)+':DebyeR:'+str(i),':'+str(hId)+':DebyeU:'+str(i)]
2640
2641            debyeDict.update(dict(zip(debyeNames,Debye['debyeTerms'][i][::2])))
2642
2643            debyeList += zip(debyeNames,Debye['debyeTerms'][i][1::2])
2644
2645        debyeVary = []
2646
2647        for item in debyeList:
2648
2649            if item[1]:
2650
2651                debyeVary.append(item[0])
2652
2653        backDict.update(debyeDict)
2654
2655        backVary += debyeVary   
2656
2657        return bakType,backDict,backVary           
2658
2659       
2660
2661    def GetInstParms(hId,Inst):
2662
2663        insVals,insFlags,insNames = Inst[1:4]
2664
2665        dataType = insVals[0]
2666
2667        instDict = {}
2668
2669        insVary = []
2670
2671        pfx = ':'+str(hId)+':'
2672
2673        for i,flag in enumerate(insFlags):
2674
2675            insName = pfx+insNames[i]
2676
2677            instDict[insName] = insVals[i]
2678
2679            if flag:
2680
2681                insVary.append(insName)
2682
2683        instDict[pfx+'X'] = max(instDict[pfx+'X'],0.001)
2684
2685        instDict[pfx+'Y'] = max(instDict[pfx+'Y'],0.001)
2686
2687        instDict[pfx+'SH/L'] = max(instDict[pfx+'SH/L'],0.0005)
2688
2689        return dataType,instDict,insVary
2690
2691       
2692
2693    def GetSampleParms(hId,Sample):
2694
2695        sampVary = []
2696
2697        hfx = ':'+str(hId)+':'       
2698
2699        sampDict = {hfx+'Gonio. radius':Sample['Gonio. radius'],hfx+'Omega':Sample['Omega'],
2700
2701            hfx+'Chi':Sample['Chi'],hfx+'Phi':Sample['Phi']}
2702
2703        Type = Sample['Type']
2704
2705        if 'Bragg' in Type:             #Bragg-Brentano
2706
2707            for item in ['Scale','Shift','Transparency']:       #surface roughness?, diffuse scattering?
2708
2709                sampDict[hfx+item] = Sample[item][0]
2710
2711                if Sample[item][1]:
2712
2713                    sampVary.append(hfx+item)
2714
2715        elif 'Debye' in Type:        #Debye-Scherrer
2716
2717            for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
2718
2719                sampDict[hfx+item] = Sample[item][0]
2720
2721                if Sample[item][1]:
2722
2723                    sampVary.append(hfx+item)
2724
2725        return Type,sampDict,sampVary
2726
2727       
2728
2729    def PrintBackground(Background):
2730
2731        Back = Background[0]
2732
2733        Debye = Background[1]
2734
2735        print '\n Background function: ',Back[0],' Refine?',bool(Back[1])
2736
2737        line = ' Coefficients: '
2738
2739        for i,back in enumerate(Back[3:]):
2740
2741            line += '%10.3f'%(back)
2742
2743            if i and not i%10:
2744
2745                line += '\n'+15*' '
2746
2747        print line
2748
2749        if Debye['nDebye']:
2750
2751            print '\n Debye diffuse scattering coefficients'
2752
2753            parms = ['DebyeA','DebyeR','DebyeU']
2754
2755            line = ' names :'
2756
2757            for parm in parms:
2758
2759                line += '%16s'%(parm)
2760
2761            print line
2762
2763            for j,term in enumerate(Debye['debyeTerms']):
2764
2765                line = ' term'+'%2d'%(j)+':'
2766
2767                for i in range(3):
2768
2769                    line += '%10.4g %5s'%(term[2*i],bool(term[2*i+1]))                   
2770
2771                print line
2772
2773       
2774
2775    def PrintInstParms(Inst):
2776
2777        print '\n Instrument Parameters:'
2778
2779        ptlbls = ' name  :'
2780
2781        ptstr =  ' value :'
2782
2783        varstr = ' refine:'
2784
2785        instNames = Inst[3][1:]
2786
2787        for i,name in enumerate(instNames):
2788
2789            ptlbls += '%12s' % (name)
2790
2791            ptstr += '%12.6f' % (Inst[1][i+1])
2792
2793            if name in ['Lam1','Lam2','Azimuth']:
2794
2795                varstr += 12*' '
2796
2797            else:
2798
2799                varstr += '%12s' % (str(bool(Inst[2][i+1])))
2800
2801        print ptlbls
2802
2803        print ptstr
2804
2805        print varstr
2806
2807       
2808
2809    def PrintSampleParms(Sample):
2810
2811        print '\n Sample Parameters:'
2812
2813        print ' Goniometer omega = %.2f, chi = %.2f, phi = %.2f'% \
2814
2815            (Sample['Omega'],Sample['Chi'],Sample['Phi'])
2816
2817        ptlbls = ' name  :'
2818
2819        ptstr =  ' value :'
2820
2821        varstr = ' refine:'
2822
2823        if 'Bragg' in Sample['Type']:
2824
2825            for item in ['Scale','Shift','Transparency']:
2826
2827                ptlbls += '%14s'%(item)
2828
2829                ptstr += '%14.4f'%(Sample[item][0])
2830
2831                varstr += '%14s'%(str(bool(Sample[item][1])))
2832
2833           
2834
2835        elif 'Debye' in Type:        #Debye-Scherrer
2836
2837            for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
2838
2839                ptlbls += '%14s'%(item)
2840
2841                ptstr += '%14.4f'%(Sample[item][0])
2842
2843                varstr += '%14s'%(str(bool(Sample[item][1])))
2844
2845
2846
2847        print ptlbls
2848
2849        print ptstr
2850
2851        print varstr
2852
2853       
2854
2855
2856
2857    histDict = {}
2858
2859    histVary = []
2860
2861    controlDict = {}
2862
2863    histoList = Histograms.keys()
2864
2865    histoList.sort()
2866
2867    for histogram in histoList:
2868
2869        Histogram = Histograms[histogram]
2870
2871        hId = Histogram['hId']
2872
2873        pfx = ':'+str(hId)+':'
2874
2875        controlDict[pfx+'Limits'] = Histogram['Limits'][1]
2876
2877       
2878
2879        Background = Histogram['Background']
2880
2881        Type,bakDict,bakVary = GetBackgroundParms(hId,Background)
2882
2883        controlDict[pfx+'bakType'] = Type
2884
2885        histDict.update(bakDict)
2886
2887        histVary += bakVary
2888
2889       
2890
2891        Inst = Histogram['Instrument Parameters']
2892
2893        Type,instDict,insVary = GetInstParms(hId,Inst)
2894
2895        controlDict[pfx+'histType'] = Type
2896
2897        if pfx+'Lam1' in instDict:
2898
2899            controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam1']
2900
2901        else:
2902
2903            controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam']           
2904
2905        histDict.update(instDict)
2906
2907        histVary += insVary
2908
2909       
2910
2911        Sample = Histogram['Sample Parameters']
2912
2913        Type,sampDict,sampVary = GetSampleParms(hId,Sample)
2914
2915        controlDict[pfx+'instType'] = Type
2916
2917        histDict.update(sampDict)
2918
2919        histVary += sampVary
2920
2921
2922
2923        if Print: 
2924
2925            print '\n Histogram: ',histogram,' histogram Id: ',hId
2926
2927            print 135*'-'
2928
2929            Units = {'C':' deg','T':' msec'}
2930
2931            units = Units[controlDict[pfx+'histType'][2]]
2932
2933            Limits = controlDict[pfx+'Limits']
2934
2935            print ' Instrument type: ',Sample['Type']
2936
2937            print ' Histogram limits: %8.2f%s to %8.2f%s'%(Limits[0],units,Limits[1],units)     
2938
2939            PrintSampleParms(Sample)
2940
2941            PrintInstParms(Inst)
2942
2943            PrintBackground(Background)
2944
2945       
2946
2947    return histVary,histDict,controlDict
2948
2949   
2950
2951def SetHistogramData(parmDict,sigDict,Histograms,Print=True):
2952
2953   
2954
2955    def SetBackgroundParms(pfx,Background,parmDict,sigDict):
2956
2957        Back = Background[0]
2958
2959        Debye = Background[1]
2960
2961        lenBack = len(Back[3:])
2962
2963        backSig = [0 for i in range(lenBack+3*Debye['nDebye'])]
2964
2965        for i in range(lenBack):
2966
2967            Back[3+i] = parmDict[pfx+'Back:'+str(i)]
2968
2969            if pfx+'Back:'+str(i) in sigDict:
2970
2971                backSig[i] = sigDict[pfx+'Back:'+str(i)]
2972
2973        if Debye['nDebye']:
2974
2975            for i in range(Debye['nDebye']):
2976
2977                names = [pfx+'DebyeA:'+str(i),pfx+'DebyeR:'+str(i),pfx+'DebyeU:'+str(i)]
2978
2979                for j,name in enumerate(names):
2980
2981                    Debye['debyeTerms'][i][2*j] = parmDict[name]
2982
2983                    if name in sigDict:
2984
2985                        backSig[lenBack+3*i+j] = sigDict[name]           
2986
2987        return backSig
2988
2989       
2990
2991    def SetInstParms(pfx,Inst,parmDict,sigDict):
2992
2993        insVals,insFlags,insNames = Inst[1:4]
2994
2995        instSig = [0 for i in range(len(insVals))]
2996
2997        for i,flag in enumerate(insFlags):
2998
2999            insName = pfx+insNames[i]
3000
3001            insVals[i] = parmDict[insName]
3002
3003            if insName in sigDict:
3004
3005                instSig[i] = sigDict[insName]
3006
3007        return instSig
3008
3009       
3010
3011    def SetSampleParms(pfx,Sample,parmDict,sigDict):
3012
3013        if 'Bragg' in Sample['Type']:             #Bragg-Brentano
3014
3015            sampSig = [0 for i in range(3)]
3016
3017            for i,item in enumerate(['Scale','Shift','Transparency']):       #surface roughness?, diffuse scattering?
3018
3019                Sample[item][0] = parmDict[pfx+item]
3020
3021                if pfx+item in sigDict:
3022
3023                    sampSig[i] = sigDict[pfx+item]
3024
3025        elif 'Debye' in Sample['Type']:        #Debye-Scherrer
3026
3027            sampSig = [0 for i in range(4)]
3028
3029            for i,item in enumerate(['Scale','Absorption','DisplaceX','DisplaceY']):
3030
3031                Sample[item][0] = parmDict[pfx+item]
3032
3033                if pfx+item in sigDict:
3034
3035                    sampSig[i] = sigDict[pfx+item]
3036
3037        return sampSig
3038
3039       
3040
3041    def PrintBackgroundSig(Background,backSig):
3042
3043        Back = Background[0]
3044
3045        Debye = Background[1]
3046
3047        lenBack = len(Back[3:])
3048
3049        valstr = ' value : '
3050
3051        sigstr = ' sig   : '
3052
3053        refine = False
3054
3055        for i,back in enumerate(Back[3:]):
3056
3057            valstr += '%10.4g'%(back)
3058
3059            if Back[1]:
3060
3061                refine = True
3062
3063                sigstr += '%10.4g'%(backSig[i])
3064
3065            else:
3066
3067                sigstr += 10*' '
3068
3069        if refine:
3070
3071            print '\n Background function: ',Back[0]
3072
3073            print valstr
3074
3075            print sigstr
3076
3077        if Debye['nDebye']:
3078
3079            ifAny = False
3080
3081            ptfmt = "%12.5f"
3082
3083            names =  ' names :'
3084
3085            ptstr =  ' values:'
3086
3087            sigstr = ' esds  :'
3088
3089            for item in sigDict:
3090
3091                if 'Debye' in item:
3092
3093                    ifAny = True
3094
3095                    names += '%12s'%(item)
3096
3097                    ptstr += ptfmt%(parmDict[item])
3098
3099                    sigstr += ptfmt%(sigDict[item])
3100
3101            if ifAny:
3102
3103                print '\n Debye diffuse scattering coefficients'
3104
3105                print names
3106
3107                print ptstr
3108
3109                print sigstr
3110
3111       
3112
3113    def PrintInstParmsSig(Inst,instSig):
3114
3115        ptlbls = ' names :'
3116
3117        ptstr =  ' value :'
3118
3119        sigstr = ' sig   :'
3120
3121        instNames = Inst[3][1:]
3122
3123        refine = False
3124
3125        for i,name in enumerate(instNames):
3126
3127            ptlbls += '%12s' % (name)
3128
3129            ptstr += '%12.6f' % (Inst[1][i+1])
3130
3131            if instSig[i+1]:
3132
3133                refine = True
3134
3135                sigstr += '%12.6f' % (instSig[i+1])
3136
3137            else:
3138
3139                sigstr += 12*' '
3140
3141        if refine:
3142
3143            print '\n Instrument Parameters:'
3144
3145            print ptlbls
3146
3147            print ptstr
3148
3149            print sigstr
3150
3151       
3152
3153    def PrintSampleParmsSig(Sample,sampleSig):
3154
3155        ptlbls = ' names :'
3156
3157        ptstr =  ' values:'
3158
3159        sigstr = ' sig   :'
3160
3161        refine = False
3162
3163        if 'Bragg' in Sample['Type']:
3164
3165            for i,item in enumerate(['Scale','Shift','Transparency']):
3166
3167                ptlbls += '%14s'%(item)
3168
3169                ptstr += '%14.4f'%(Sample[item][0])
3170
3171                if sampleSig[i]:
3172
3173                    refine = True
3174
3175                    sigstr += '%14.4f'%(sampleSig[i])
3176
3177                else:
3178
3179                    sigstr += 14*' '
3180
3181           
3182
3183        elif 'Debye' in Sample['Type']:        #Debye-Scherrer
3184
3185            for i,item in enumerate(['Scale','Absorption','DisplaceX','DisplaceY']):
3186
3187                ptlbls += '%14s'%(item)
3188
3189                ptstr += '%14.4f'%(Sample[item][0])
3190
3191                if sampleSig[i]:
3192
3193                    refine = True
3194
3195                    sigstr += '%14.4f'%(sampleSig[i])
3196
3197                else:
3198
3199                    sigstr += 14*' '
3200
3201
3202
3203        if refine:
3204
3205            print '\n Sample Parameters:'
3206
3207            print ptlbls
3208
3209            print ptstr
3210
3211            print sigstr
3212
3213       
3214
3215    histoList = Histograms.keys()
3216
3217    histoList.sort()
3218
3219    for histogram in histoList:
3220
3221        if 'PWDR' in histogram:
3222
3223            Histogram = Histograms[histogram]
3224
3225            hId = Histogram['hId']
3226
3227            pfx = ':'+str(hId)+':'
3228
3229            Background = Histogram['Background']
3230
3231            backSig = SetBackgroundParms(pfx,Background,parmDict,sigDict)
3232
3233           
3234
3235            Inst = Histogram['Instrument Parameters']
3236
3237            instSig = SetInstParms(pfx,Inst,parmDict,sigDict)
3238
3239       
3240
3241            Sample = Histogram['Sample Parameters']
3242
3243            sampSig = SetSampleParms(pfx,Sample,parmDict,sigDict)
3244
3245
3246
3247            print '\n Histogram: ',histogram,' histogram Id: ',hId
3248
3249            print 135*'-'
3250
3251            print ' Final refinement wRp = %.2f%% on %d observations in this histogram'%(Histogram['wRp'],Histogram['Nobs'])
3252
3253            if Print:
3254
3255                print ' Instrument type: ',Sample['Type']
3256
3257                PrintSampleParmsSig(Sample,sampSig)
3258
3259                PrintInstParmsSig(Inst,instSig)
3260
3261                PrintBackgroundSig(Background,backSig)
3262
3263
3264
3265def GetAtomFXU(pfx,FFtables,BLtables,calcControls,parmDict):
3266
3267    Natoms = calcControls['Natoms'][pfx]
3268
3269    Tdata = Natoms*[' ',]
3270
3271    Mdata = np.zeros(Natoms)
3272
3273    IAdata = Natoms*[' ',]
3274
3275    Fdata = np.zeros(Natoms)
3276
3277    FFdata = []
3278
3279    BLdata = []
3280
3281    Xdata = np.zeros((3,Natoms))
3282
3283    dXdata = np.zeros((3,Natoms))
3284
3285    Uisodata = np.zeros(Natoms)
3286
3287    Uijdata = np.zeros((6,Natoms))
3288
3289    keys = {'Atype:':Tdata,'Amul:':Mdata,'Afrac:':Fdata,'AI/A:':IAdata,
3290
3291        'dAx:':dXdata[0],'dAy:':dXdata[1],'dAz:':dXdata[2],
3292
3293        'Ax:':Xdata[0],'Ay:':Xdata[1],'Az:':Xdata[2],'AUiso:':Uisodata,
3294
3295        'AU11:':Uijdata[0],'AU22:':Uijdata[1],'AU33:':Uijdata[2],
3296
3297        'AU12:':Uijdata[3],'AU13:':Uijdata[4],'AU23:':Uijdata[5]}
3298
3299    for iatm in range(Natoms):
3300
3301        for key in keys:
3302
3303            parm = pfx+key+str(iatm)
3304
3305            if parm in parmDict:
3306
3307                keys[key][iatm] = parmDict[parm]
3308
3309        FFdata.append(FFtables[Tdata[iatm]])
3310
3311        BLdata.append(BLtables[Tdata[iatm]][1])
3312
3313    return FFdata,BLdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata
3314
3315   
3316
3317def StructureFactor(refList,G,hfx,pfx,SGData,calcControls,parmDict):
3318
3319    ''' Compute structure factors for all h,k,l for phase
3320
3321    input:
3322
3323        refList: [ref] where each ref = h,k,l,m,d,...,[equiv h,k,l],phase[equiv]
3324
3325        G:      reciprocal metric tensor
3326
3327        pfx:    phase id string
3328
3329        SGData: space group info. dictionary output from SpcGroup
3330
3331        calcControls:
3332
3333        ParmDict:
3334
3335    puts result F^2 in each ref[8] in refList
3336
3337    '''       
3338
3339    twopi = 2.0*np.pi
3340
3341    twopisq = 2.0*np.pi**2
3342
3343    ast = np.sqrt(np.diag(G))
3344
3345    Mast = twopisq*np.multiply.outer(ast,ast)
3346
3347    FFtables = calcControls['FFtables']
3348
3349    BLtables = calcControls['BLtables']
3350
3351    FFdata,BLdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,FFtables,BLtables,calcControls,parmDict)
3352
3353    if 'N' in parmDict[hfx+'Type']:
3354
3355        FP,FPP = G2el.BlenRes(BLdata,parmDict[hfx+'Lam'])
3356
3357    else:
3358
3359        FP = np.array([El[hfx+'FP'] for El in FFdata])
3360
3361        FPP = np.array([El[hfx+'FPP'] for El in FFdata])
3362
3363    maxPos = len(SGData['SGOps'])
3364
3365    Uij = np.array(G2lat.U6toUij(Uijdata))
3366
3367    bij = Mast*Uij.T
3368
3369    for refl in refList:
3370
3371        fbs = np.array([0,0])
3372
3373        H = refl[:3]
3374
3375        SQ = 1./(2.*refl[4])**2
3376
3377        if 'N' in parmDict[hfx+'Type']:
3378
3379            FF = np.array([El[1] for El in BLdata])
3380
3381        else:       #'X'
3382
3383            FF = np.array([G2el.ScatFac(El,SQ)[0] for El in FFdata])
3384
3385        SQfactor = 4.0*SQ*twopisq
3386
3387        Uniq = refl[11]
3388
3389        phi = refl[12]
3390
3391        phase = twopi*(np.inner(Uniq,(dXdata.T+Xdata.T))+phi[:,np.newaxis])
3392
3393        sinp = np.sin(phase)
3394
3395        cosp = np.cos(phase)
3396
3397        occ = Mdata*Fdata/len(Uniq)
3398
3399        biso = -SQfactor*Uisodata
3400
3401        Tiso = np.where(biso<1.,np.exp(biso),1.0)
3402
3403        HbH = np.array([-np.inner(h,np.inner(bij,h)) for h in Uniq])
3404
3405        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
3406
3407        Tcorr = Tiso*Tuij
3408
3409        fa = np.array([(FF+FP)*occ*cosp*Tcorr,-FPP*occ*sinp*Tcorr])
3410
3411        fas = np.sum(np.sum(fa,axis=1),axis=1)        #real
3412
3413        if not SGData['SGInv']:
3414
3415            fb = np.array([(FF+FP)*occ*sinp*Tcorr,FPP*occ*cosp*Tcorr])
3416
3417            fbs = np.sum(np.sum(fb,axis=1),axis=1)
3418
3419        fasq = fas**2
3420
3421        fbsq = fbs**2        #imaginary
3422
3423        refl[9] = np.sum(fasq)+np.sum(fbsq)
3424
3425        refl[10] = atan2d(fbs[0],fas[0])
3426
3427    return refList
3428
3429   
3430
3431def StructureFactorDerv(refList,G,hfx,pfx,SGData,calcControls,parmDict):
3432
3433    twopi = 2.0*np.pi
3434
3435    twopisq = 2.0*np.pi**2
3436
3437    ast = np.sqrt(np.diag(G))
3438
3439    Mast = twopisq*np.multiply.outer(ast,ast)
3440
3441    FFtables = calcControls['FFtables']
3442
3443    BLtables = calcControls['BLtables']
3444
3445    FFdata,BLdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,FFtables,BLtables,calcControls,parmDict)
3446
3447    if 'N' in parmDict[hfx+'Type']:
3448
3449        FP = 0.
3450
3451        FPP = 0.
3452
3453    else:
3454
3455        FP = np.array([El[hfx+'FP'] for El in FFdata])
3456
3457        FPP = np.array([El[hfx+'FPP'] for El in FFdata])
3458
3459    maxPos = len(SGData['SGOps'])       
3460
3461    Uij = np.array(G2lat.U6toUij(Uijdata))
3462
3463    bij = Mast*Uij.T
3464
3465    dFdvDict = {}
3466
3467    dFdfr = np.zeros((len(refList),len(Mdata)))
3468
3469    dFdx = np.zeros((len(refList),len(Mdata),3))
3470
3471    dFdui = np.zeros((len(refList),len(Mdata)))
3472
3473    dFdua = np.zeros((len(refList),len(Mdata),6))
3474
3475    for iref,refl in enumerate(refList):
3476
3477        H = np.array(refl[:3])
3478
3479        SQ = 1./(2.*refl[4])**2          # or (sin(theta)/lambda)**2
3480
3481        if 'N' in parmDict[hfx+'Type']:
3482
3483            FF = np.array([El[1] for El in BLdata])
3484
3485        else:       #'X'
3486
3487            FF = np.array([G2el.ScatFac(El,SQ)[0] for El in FFdata])
3488
3489        SQfactor = 8.0*SQ*np.pi**2
3490
3491        Uniq = refl[11]
3492
3493        phi = refl[12]
3494
3495        phase = twopi*(np.inner((dXdata.T+Xdata.T),Uniq)+phi[np.newaxis,:])
3496
3497        sinp = np.sin(phase)
3498
3499        cosp = np.cos(phase)
3500
3501        occ = Mdata*Fdata/len(Uniq)
3502
3503        biso = -SQfactor*Uisodata
3504
3505        Tiso = np.where(biso<1.,np.exp(biso),1.0)
3506
3507#        HbH = np.array([-np.inner(h,np.inner(bij,h)) for h in Uniq])
3508
3509        HbH = -np.inner(H,np.inner(bij,H))
3510
3511        Hij = np.array([Mast*np.multiply.outer(U,U) for U in Uniq])
3512
3513        Hij = np.array([G2lat.UijtoU6(Uij) for Uij in Hij])
3514
3515        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
3516
3517        Tcorr = Tiso*Tuij
3518
3519        fot = (FF+FP)*occ*Tcorr
3520
3521        fotp = FPP*occ*Tcorr
3522
3523        fa = np.array([fot[:,np.newaxis]*cosp,fotp[:,np.newaxis]*cosp])       #non positions
3524
3525        fb = np.array([fot[:,np.newaxis]*sinp,-fotp[:,np.newaxis]*sinp])
3526
3527       
3528
3529        fas = np.sum(np.sum(fa,axis=1),axis=1)
3530
3531        fbs = np.sum(np.sum(fb,axis=1),axis=1)
3532
3533        fax = np.array([-fot[:,np.newaxis]*sinp,-fotp[:,np.newaxis]*sinp])   #positions
3534
3535        fbx = np.array([fot[:,np.newaxis]*cosp,-fot[:,np.newaxis]*cosp])
3536
3537        #sum below is over Uniq
3538
3539        dfadfr = np.sum(fa/occ[:,np.newaxis],axis=2)
3540
3541        dfadx = np.sum(twopi*Uniq*fax[:,:,:,np.newaxis],axis=2)
3542
3543        dfadui = np.sum(-SQfactor*fa,axis=2)
3544
3545        dfadua = np.sum(-Hij*fa[:,:,:,np.newaxis],axis=2)
3546
3547        #NB: the above have been checked against PA(1:10,1:2) in strfctr.for     
3548
3549        dFdfr[iref] = 2.*(fas[0]*dfadfr[0]+fas[1]*dfadfr[1])*Mdata/len(Uniq)
3550
3551        dFdx[iref] = 2.*(fas[0]*dfadx[0]+fas[1]*dfadx[1])
3552
3553        dFdui[iref] = 2.*(fas[0]*dfadui[0]+fas[1]*dfadui[1])
3554
3555        dFdua[iref] = 2.*(fas[0]*dfadua[0]+fas[1]*dfadua[1])
3556
3557        if not SGData['SGInv']:
3558
3559            dfbdfr = np.sum(fb/occ[:,np.newaxis],axis=2)        #problem here if occ=0 for some atom
3560
3561            dfbdx = np.sum(twopi*Uniq*fbx[:,:,:,np.newaxis],axis=2)         
3562
3563            dfbdui = np.sum(-SQfactor*fb,axis=2)
3564
3565            dfbdua = np.sum(-Hij*fb[:,:,:,np.newaxis],axis=2)
3566
3567            dFdfr[iref] += 2.*(fbs[0]*dfbdfr[0]-fbs[1]*dfbdfr[1])*Mdata/len(Uniq)
3568
3569            dFdx[iref] += 2.*(fbs[0]*dfbdx[0]+fbs[1]*dfbdx[1])
3570
3571            dFdui[iref] += 2.*(fbs[0]*dfbdui[0]-fbs[1]*dfbdui[1])
3572
3573            dFdua[iref] += 2.*(fbs[0]*dfbdua[0]+fbs[1]*dfbdua[1])
3574
3575        #loop over atoms - each dict entry is list of derivatives for all the reflections
3576
3577    for i in range(len(Mdata)):     
3578
3579        dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i]
3580
3581        dFdvDict[pfx+'dAx:'+str(i)] = dFdx.T[0][i]
3582
3583        dFdvDict[pfx+'dAy:'+str(i)] = dFdx.T[1][i]
3584
3585        dFdvDict[pfx+'dAz:'+str(i)] = dFdx.T[2][i]
3586
3587        dFdvDict[pfx+'AUiso:'+str(i)] = dFdui.T[i]
3588
3589        dFdvDict[pfx+'AU11:'+str(i)] = dFdua.T[0][i]
3590
3591        dFdvDict[pfx+'AU22:'+str(i)] = dFdua.T[1][i]
3592
3593        dFdvDict[pfx+'AU33:'+str(i)] = dFdua.T[2][i]
3594
3595        dFdvDict[pfx+'AU12:'+str(i)] = 2.*dFdua.T[3][i]
3596
3597        dFdvDict[pfx+'AU13:'+str(i)] = 2.*dFdua.T[4][i]
3598
3599        dFdvDict[pfx+'AU23:'+str(i)] = 2.*dFdua.T[5][i]
3600
3601    return dFdvDict
3602
3603       
3604
3605def Dict2Values(parmdict, varylist):
3606
3607    '''Use before call to leastsq to setup list of values for the parameters
3608
3609    in parmdict, as selected by key in varylist'''
3610
3611    return [parmdict[key] for key in varylist] 
3612
3613   
3614
3615def Values2Dict(parmdict, varylist, values):
3616
3617    ''' Use after call to leastsq to update the parameter dictionary with
3618
3619    values corresponding to keys in varylist'''
3620
3621    parmdict.update(zip(varylist,values))
3622
3623   
3624
3625def GetNewCellParms(parmDict,varyList):
3626
3627    newCellDict = {}
3628
3629    Ddict = dict(zip(['D11','D22','D33','D12','D13','D23'],['A'+str(i) for i in range(6)]))
3630
3631    for item in varyList:
3632
3633        keys = item.split(':')
3634
3635        if keys[2] in Ddict:
3636
3637            key = keys[0]+'::'+Ddict[keys[2]]
3638
3639            parm = keys[0]+'::'+keys[2]
3640
3641            newCellDict[parm] = [key,parmDict[key]+parmDict[item]]
3642
3643    return newCellDict
3644
3645   
3646
3647def ApplyXYZshifts(parmDict,varyList):
3648
3649    ''' takes atom x,y,z shift and applies it to corresponding atom x,y,z value
3650
3651        input:
3652
3653            parmDict - parameter dictionary
3654
3655            varyList - list of variables
3656
3657        returns:
3658
3659            newAtomDict - dictitemionary of new atomic coordinate names & values;
3660
3661                key is parameter shift name
3662
3663    '''
3664
3665    newAtomDict = {}
3666
3667    for item in parmDict:
3668
3669        if 'dA' in item:
3670
3671            parm = ''.join(item.split('d'))
3672
3673            parmDict[parm] += parmDict[item]
3674
3675            newAtomDict[item] = [parm,parmDict[parm]]
3676
3677    return newAtomDict
3678
3679   
3680
3681def SHTXcal(refl,g,pfx,hfx,SGData,calcControls,parmDict):
3682
3683    IFCoup = 'Bragg' in calcControls[hfx+'instType']
3684
3685    odfCor = 1.0
3686
3687    H = refl[:3]
3688
3689    cell = G2lat.Gmat2cell(g)
3690
3691    Sangls = [parmDict[pfx+'SH omega'],parmDict[pfx+'SH chi'],parmDict[pfx+'SH phi']]
3692
3693    Gangls = [parmDict[hfx+'Omega'],parmDict[hfx+'Chi'],parmDict[hfx+'Phi'],parmDict[hfx+'Azimuth']]
3694
3695    phi,beta = G2lat.CrsAng(H,cell,SGData)
3696
3697    psi,gam,x,x = G2lat.SamAng(refl[5]/2.,Gangls,Sangls,IFCoup) #ignore 2 sets of angle derivs.
3698
3699    SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],parmDict[pfx+'SHmodel'],parmDict[pfx+'SHorder'])
3700
3701    for item in SHnames:
3702
3703        L,M,N = eval(item.strip('C'))
3704
3705        Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
3706
3707        Ksl,x,x = G2lat.GetKsl(L,M,parmDict[pfx+'SHmodel'],psi,gam)
3708
3709        Lnorm = G2lat.Lnorm(L)
3710
3711        odfCor += parmDict[pfx+item]*Lnorm*Kcl*Ksl
3712
3713    return odfCor
3714
3715   
3716
3717def SHTXcalDerv(refl,g,pfx,hfx,SGData,calcControls,parmDict):
3718
3719    FORPI = 12.5663706143592
3720
3721    IFCoup = 'Bragg' in calcControls[hfx+'instType']
3722
3723    odfCor = 1.0
3724
3725    dFdODF = {}
3726
3727    dFdSA = [0,0,0]
3728
3729    H = refl[:3]
3730
3731    cell = G2lat.Gmat2cell(g)
3732
3733    Sangls = [parmDict[pfx+'SH omega'],parmDict[pfx+'SH chi'],parmDict[pfx+'SH phi']]
3734
3735    Gangls = [parmDict[hfx+'Omega'],parmDict[hfx+'Chi'],parmDict[hfx+'Phi'],parmDict[hfx+'Azimuth']]
3736
3737    phi,beta = G2lat.CrsAng(H,cell,SGData)
3738
3739    psi,gam,dPSdA,dGMdA = G2lat.SamAng(refl[5]/2.,Gangls,Sangls,IFCoup)
3740
3741    SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],parmDict[pfx+'SHmodel'],parmDict[pfx+'SHorder'])
3742
3743    for item in SHnames:
3744
3745        L,M,N = eval(item.strip('C'))
3746
3747        Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
3748
3749        Ksl,dKsdp,dKsdg = G2lat.GetKsl(L,M,parmDict[pfx+'SHmodel'],psi,gam)
3750
3751        Lnorm = G2lat.Lnorm(L)
3752
3753        odfCor += parmDict[pfx+item]*Lnorm*Kcl*Ksl
3754
3755        dFdODF[pfx+item] = Lnorm*Kcl*Ksl
3756
3757        for i in range(3):
3758
3759            dFdSA[i] += parmDict[pfx+item]*Lnorm*Kcl*(dKsdp*dPSdA[i]+dKsdg*dGMdA[i])
3760
3761    return odfCor,dFdODF,dFdSA
3762
3763   
3764
3765def SHPOcal(refl,g,phfx,hfx,SGData,calcControls,parmDict):
3766
3767    odfCor = 1.0
3768
3769    H = refl[:3]
3770
3771    cell = G2lat.Gmat2cell(g)
3772
3773    Sangl = [0.,0.,0.]
3774
3775    if 'Bragg' in calcControls[hfx+'instType']:
3776
3777        Gangls = [0.,90.,0.,parmDict[hfx+'Azimuth']]
3778
3779        IFCoup = True
3780
3781    else:
3782
3783        Gangls = [0.,0.,0.,parmDict[hfx+'Azimuth']]
3784
3785        IFCoup = False
3786
3787    phi,beta = G2lat.CrsAng(H,cell,SGData)
3788
3789    psi,gam,x,x = G2lat.SamAng(refl[5]/2.,Gangls,Sangl,IFCoup) #ignore 2 sets of angle derivs.
3790
3791    SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],'0',calcControls[phfx+'SHord'],False)
3792
3793    for item in SHnames:
3794
3795        L,N = eval(item.strip('C'))
3796
3797        Kcsl,Lnorm = G2lat.GetKclKsl(L,N,SGData['SGLaue'],psi,phi,beta) 
3798
3799        odfCor += parmDict[phfx+item]*Lnorm*Kcsl
3800
3801    return odfCor
3802
3803   
3804
3805def SHPOcalDerv(refl,g,phfx,hfx,SGData,calcControls,parmDict):
3806
3807    FORPI = 12.5663706143592
3808
3809    odfCor = 1.0
3810
3811    dFdODF = {}
3812
3813    H = refl[:3]
3814
3815    cell = G2lat.Gmat2cell(g)
3816
3817    Sangl = [0.,0.,0.]
3818
3819    if 'Bragg' in calcControls[hfx+'instType']:
3820
3821        Gangls = [0.,90.,0.,parmDict[hfx+'Azimuth']]
3822
3823        IFCoup = True
3824
3825    else:
3826
3827        Gangls = [0.,0.,0.,parmDict[hfx+'Azimuth']]
3828
3829        IFCoup = False
3830
3831    phi,beta = G2lat.CrsAng(H,cell,SGData)
3832
3833    psi,gam,x,x = G2lat.SamAng(refl[5]/2.,Gangls,Sangl,IFCoup) #ignore 2 sets of angle derivs.
3834
3835    SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],'0',calcControls[phfx+'SHord'],False)
3836
3837    for item in SHnames:
3838
3839        L,N = eval(item.strip('C'))
3840
3841        Kcsl,Lnorm = G2lat.GetKclKsl(L,N,SGData['SGLaue'],psi,phi,beta) 
3842
3843        odfCor += parmDict[phfx+item]*Lnorm*Kcsl
3844
3845        dFdODF[phfx+item] = Kcsl*Lnorm
3846
3847    return odfCor,dFdODF
3848
3849   
3850
3851def GetPrefOri(refl,G,g,phfx,hfx,SGData,calcControls,parmDict):
3852
3853    if calcControls[phfx+'poType'] == 'MD':
3854
3855        MD = parmDict[phfx+'MD']
3856
3857        MDAxis = calcControls[phfx+'MDAxis']
3858
3859        sumMD = 0
3860
3861        for H in refl[11]:           
3862
3863            cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G)
3864
3865            A = 1.0/np.sqrt((MD*cosP)**2+sinP**2/MD)
3866
3867            sumMD += A**3
3868
3869        POcorr = sumMD/len(refl[11])
3870
3871    else:   #spherical harmonics
3872
3873        POcorr = SHPOcal(refl,g,phfx,hfx,SGData,calcControls,parmDict)
3874
3875    return POcorr
3876
3877   
3878
3879def GetPrefOriDerv(refl,G,g,phfx,hfx,SGData,calcControls,parmDict):
3880
3881    POderv = {}
3882
3883    if calcControls[phfx+'poType'] == 'MD':
3884
3885        MD = parmDict[phfx+'MD']
3886
3887        MDAxis = calcControls[phfx+'MDAxis']
3888
3889        sumMD = 0
3890
3891        sumdMD = 0
3892
3893        for H in refl[11]:           
3894
3895            cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G)
3896
3897            A = 1.0/np.sqrt((MD*cosP)**2+sinP**2/MD)
3898
3899            sumMD += A**3
3900
3901            sumdMD -= (1.5*A**5)*(2.0*MD*cosP**2-(sinP/MD)**2)
3902
3903        POcorr = sumMD/len(refl[11])
3904
3905        POderv[phfx+'MD'] = sumdMD/len(refl[11])
3906
3907    else:   #spherical harmonics
3908
3909        POcorr,POderv = SHPOcalDerv(refl,g,phfx,hfx,SGData,calcControls,parmDict)
3910
3911    return POcorr,POderv
3912
3913   
3914
3915def GetIntensityCorr(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict):
3916
3917    Icorr = parmDict[phfx+'Scale']*parmDict[hfx+'Scale']*refl[3]               #scale*multiplicity
3918
3919    if 'X' in parmDict[hfx+'Type']:
3920
3921        Icorr *= G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5],parmDict[hfx+'Azimuth'])[0]
3922
3923    Icorr *= GetPrefOri(refl,G,g,phfx,hfx,SGData,calcControls,parmDict)
3924
3925    if pfx+'SHorder' in parmDict:
3926
3927        Icorr *= SHTXcal(refl,g,pfx,hfx,SGData,calcControls,parmDict)
3928
3929    refl[13] = Icorr       
3930
3931   
3932
3933def GetIntensityDerv(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict):
3934
3935    dIdsh = 1./parmDict[hfx+'Scale']
3936
3937    dIdsp = 1./parmDict[phfx+'Scale']
3938
3939    if 'X' in parmDict[hfx+'Type']:
3940
3941        pola,dIdPola = G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5],parmDict[hfx+'Azimuth'])
3942
3943        dIdPola /= pola
3944
3945    else:       #'N'
3946
3947        dIdPola = 0.0
3948
3949    POcorr,dIdPO = GetPrefOriDerv(refl,G,g,phfx,hfx,SGData,calcControls,parmDict)
3950
3951    for iPO in dIdPO:
3952
3953        dIdPO[iPO] /= POcorr
3954
3955    dFdODF = {}
3956
3957    dFdSA = [0,0,0]
3958
3959    if pfx+'SHorder' in parmDict:
3960
3961        odfCor,dFdODF,dFdSA = SHTXcalDerv(refl,g,pfx,hfx,SGData,calcControls,parmDict)
3962
3963        for iSH in dFdODF:
3964
3965            dFdODF[iSH] /= odfCor
3966
3967        for i in range(3):
3968
3969            dFdSA[i] /= odfCor
3970
3971    return dIdsh,dIdsp,dIdPola,dIdPO,dFdODF,dFdSA
3972
3973       
3974
3975def GetSampleGam(refl,wave,G,GB,phfx,calcControls,parmDict):
3976
3977    costh = cosd(refl[5]/2.)
3978
3979    #crystallite size
3980
3981    if calcControls[phfx+'SizeType'] == 'isotropic':
3982
3983        gam = 1.8*wave/(np.pi*parmDict[phfx+'Size:i']*costh)
3984
3985    elif calcControls[phfx+'SizeType'] == 'uniaxial':
3986
3987        H = np.array(refl[:3])
3988
3989        P = np.array(calcControls[phfx+'SizeAxis'])
3990
3991        cosP,sinP = G2lat.CosSinAngle(H,P,G)
3992
3993        gam = (1.8*wave/np.pi)/(parmDict[phfx+'Size:i']*parmDict[phfx+'Size:a']*costh)
3994
3995        gam *= np.sqrt((sinP*parmDict[phfx+'Size:a'])**2+(cosP*parmDict[phfx+'Size:i'])**2)
3996
3997    else:           #ellipsoidal crystallites
3998
3999        Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
4000
4001        H = np.array(refl[:3])
4002
4003        lenR = G2pwd.ellipseSize(H,Sij,GB)
4004
4005        gam = 1.8*wave/(np.pi*costh*lenR)
4006
4007    #microstrain               
4008
4009    if calcControls[phfx+'MustrainType'] == 'isotropic':
4010
4011        gam += 0.018*parmDict[phfx+'Mustrain:i']*tand(refl[5]/2.)/np.pi
4012
4013    elif calcControls[phfx+'MustrainType'] == 'uniaxial':
4014
4015        H = np.array(refl[:3])
4016
4017        P = np.array(calcControls[phfx+'MustrainAxis'])
4018
4019        cosP,sinP = G2lat.CosSinAngle(H,P,G)
4020
4021        Si = parmDict[phfx+'Mustrain:i']
4022
4023        Sa = parmDict[phfx+'Mustrain:a']
4024
4025        gam += 0.018*Si*Sa*tand(refl[5]/2.)/(np.pi*np.sqrt((Si*cosP)**2+(Sa*sinP)**2))
4026
4027    else:       #generalized - P.W. Stephens model
4028
4029        pwrs = calcControls[phfx+'MuPwrs']
4030
4031        sum = 0
4032
4033        for i,pwr in enumerate(pwrs):
4034
4035            sum += parmDict[phfx+'Mustrain:'+str(i)]*refl[0]**pwr[0]*refl[1]**pwr[1]*refl[2]**pwr[2]
4036
4037        gam += 0.018*refl[4]**2*tand(refl[5]/2.)*sum           
4038
4039    return gam
4040
4041       
4042
4043def GetSampleGamDerv(refl,wave,G,GB,phfx,calcControls,parmDict):
4044
4045    gamDict = {}
4046
4047    costh = cosd(refl[5]/2.)
4048
4049    tanth = tand(refl[5]/2.)
4050
4051    #crystallite size derivatives
4052
4053    if calcControls[phfx+'SizeType'] == 'isotropic':
4054
4055        gamDict[phfx+'Size:i'] = -1.80*wave/(np.pi*costh)
4056
4057    elif calcControls[phfx+'SizeType'] == 'uniaxial':
4058
4059        H = np.array(refl[:3])
4060
4061        P = np.array(calcControls[phfx+'SizeAxis'])
4062
4063        cosP,sinP = G2lat.CosSinAngle(H,P,G)
4064
4065        Si = parmDict[phfx+'Size:i']
4066
4067        Sa = parmDict[phfx+'Size:a']
4068
4069        gami = (1.8*wave/np.pi)/(Si*Sa)
4070
4071        sqtrm = np.sqrt((sinP*Sa)**2+(cosP*Si)**2)
4072
4073        gam = gami*sqtrm/costh           
4074
4075        gamDict[phfx+'Size:i'] = gami*Si*cosP**2/(sqtrm*costh)-gam/Si
4076
4077        gamDict[phfx+'Size:a'] = gami*Sa*sinP**2/(sqtrm*costh)-gam/Sa         
4078
4079    else:           #ellipsoidal crystallites
4080
4081        const = 1.8*wave/(np.pi*costh)
4082
4083        Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)]
4084
4085        H = np.array(refl[:3])
4086
4087        R,dRdS = G2pwd.ellipseSizeDerv(H,Sij,GB)
4088
4089        for i,item in enumerate([phfx+'Size:%d'%(j) for j in range(6)]):
4090
4091            gamDict[item] = -(const/R**2)*dRdS[i]
4092
4093    #microstrain derivatives               
4094
4095    if calcControls[phfx+'MustrainType'] == 'isotropic':
4096
4097        gamDict[phfx+'Mustrain:i'] =  0.018*tanth/np.pi           
4098
4099    elif calcControls[phfx+'MustrainType'] == 'uniaxial':
4100
4101        H = np.array(refl[:3])
4102
4103        P = np.array(calcControls[phfx+'MustrainAxis'])
4104
4105        cosP,sinP = G2lat.CosSinAngle(H,P,G)
4106
4107        Si = parmDict[phfx+'Mustrain:i']
4108
4109        Sa = parmDict[phfx+'Mustrain:a']
4110
4111        gami = 0.018*Si*Sa*tanth/np.pi
4112
4113        sqtrm = np.sqrt((Si*cosP)**2+(Sa*sinP)**2)
4114
4115        gam = gami/sqtrm
4116
4117        gamDict[phfx+'Mustrain:i'] = gam/Si-gami*Si*cosP**2/sqtrm**3
4118
4119        gamDict[phfx+'Mustrain:a'] = gam/Sa-gami*Sa*sinP**2/sqtrm**3
4120
4121    else:       #generalized - P.W. Stephens model
4122
4123        pwrs = calcControls[phfx+'MuPwrs']
4124
4125        const = 0.018*refl[4]**2*tanth
4126
4127        for i,pwr in enumerate(pwrs):
4128
4129            gamDict[phfx+'Mustrain:'+str(i)] = const*refl[0]**pwr[0]*refl[1]**pwr[1]*refl[2]**pwr[2]
4130
4131    return gamDict
4132
4133       
4134
4135def GetReflPos(refl,wave,G,hfx,calcControls,parmDict):
4136
4137    h,k,l = refl[:3]
4138
4139    dsq = 1./G2lat.calc_rDsq2(np.array([h,k,l]),G)
4140
4141    d = np.sqrt(dsq)
4142
4143
4144
4145    refl[4] = d
4146
4147    pos = 2.0*asind(wave/(2.0*d))+parmDict[hfx+'Zero']
4148
4149    const = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius'])                  #shifts in microns
4150
4151    if 'Bragg' in calcControls[hfx+'instType']:
4152
4153        pos -= const*(4.*parmDict[hfx+'Shift']*cosd(pos/2.0)+ \
4154
4155            parmDict[hfx+'Transparency']*sind(pos)*100.0)            #trans(=1/mueff) in cm
4156
4157    else:               #Debye-Scherrer - simple but maybe not right
4158
4159        pos -= const*(parmDict[hfx+'DisplaceX']*cosd(pos)+parmDict[hfx+'DisplaceY']*sind(pos))
4160
4161    return pos
4162
4163
4164
4165def GetReflPosDerv(refl,wave,A,hfx,calcControls,parmDict):
4166
4167    dpr = 180./np.pi
4168
4169    h,k,l = refl[:3]
4170
4171    dstsq = G2lat.calc_rDsq(np.array([h,k,l]),A)
4172
4173    dst = np.sqrt(dstsq)
4174
4175    pos = refl[5]-parmDict[hfx+'Zero']
4176
4177    const = dpr/np.sqrt(1.0-wave**2*dstsq/4.0)
4178
4179    dpdw = const*dst
4180
4181    dpdA = np.array([h**2,k**2,l**2,h*k,h*l,k*l])
4182
4183    dpdA *= const*wave/(2.0*dst)
4184
4185    dpdZ = 1.0
4186
4187    const = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius'])                  #shifts in microns
4188
4189    if 'Bragg' in calcControls[hfx+'instType']:
4190
4191        dpdSh = -4.*const*cosd(pos/2.0)
4192
4193        dpdTr = -const*sind(pos)*100.0
4194
4195        return dpdA,dpdw,dpdZ,dpdSh,dpdTr,0.,0.
4196
4197    else:               #Debye-Scherrer - simple but maybe not right
4198
4199        dpdXd = -const*cosd(pos)
4200
4201        dpdYd = -const*sind(pos)
4202
4203        return dpdA,dpdw,dpdZ,0.,0.,dpdXd,dpdYd
4204
4205           
4206
4207def GetHStrainShift(refl,SGData,phfx,parmDict):
4208
4209    laue = SGData['SGLaue']
4210
4211    uniq = SGData['SGUniq']
4212
4213    h,k,l = refl[:3]
4214
4215    if laue in ['m3','m3m']:
4216
4217        Dij = parmDict[phfx+'D11']*(h**2+k**2+l**2)+ \
4218
4219            refl[4]**2*parmDict[phfx+'eA']*((h*k)**2+(h*l)**2+(k*l)**2)/(h**2+k**2+l**2)**2
4220
4221    elif laue in ['6/m','6/mmm','3m1','31m','3']:
4222
4223        Dij = parmDict[phfx+'D11']*(h**2+k**2+h*k)+parmDict[phfx+'D33']*l**2
4224
4225    elif laue in ['3R','3mR']:
4226
4227        Dij = parmDict[phfx+'D11']*(h**2+k**2+l**2)+parmDict[phfx+'D12']*(h*k+h*l+k*l)
4228
4229    elif laue in ['4/m','4/mmm']:
4230
4231        Dij = parmDict[phfx+'D11']*(h**2+k**2)+parmDict[phfx+'D33']*l**2
4232
4233    elif laue in ['mmm']:
4234
4235        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2
4236
4237    elif laue in ['2/m']:
4238
4239        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2
4240
4241        if uniq == 'a':
4242
4243            Dij += parmDict[phfx+'D23']*k*l
4244
4245        elif uniq == 'b':
4246
4247            Dij += parmDict[phfx+'D13']*h*l
4248
4249        elif uniq == 'c':
4250
4251            Dij += parmDict[phfx+'D12']*h*k
4252
4253    else:
4254
4255        Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2+ \
4256
4257            parmDict[phfx+'D12']*h*k+parmDict[phfx+'D13']*h*l+parmDict[phfx+'D23']*k*l
4258
4259    return Dij*refl[4]**2*tand(refl[5]/2.0)
4260
4261           
4262
4263def GetHStrainShiftDerv(refl,SGData,phfx):
4264
4265    laue = SGData['SGLaue']
4266
4267    uniq = SGData['SGUniq']
4268
4269    h,k,l = refl[:3]
4270
4271    if laue in ['m3','m3m']:
4272
4273        dDijDict = {phfx+'D11':h**2+k**2+l**2,
4274
4275            phfx+'eA':((h*k)**2+(h*l)**2+(k*l)**2)/(h**2+k**2+l**2)**2}
4276
4277    elif laue in ['6/m','6/mmm','3m1','31m','3']:
4278
4279        dDijDict = {phfx+'D11':h**2+k**2+h*k,phfx+'D33':l**2}
4280
4281    elif laue in ['3R','3mR']:
4282
4283        dDijDict = {phfx+'D11':h**2+k**2+l**2,phfx+'D12':h*k+h*l+k*l}
4284
4285    elif laue in ['4/m','4/mmm']:
4286
4287        dDijDict = {phfx+'D11':h**2+k**2,phfx+'D33':l**2}
4288
4289    elif laue in ['mmm']:
4290
4291        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2}
4292
4293    elif laue in ['2/m']:
4294
4295        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2}
4296
4297        if uniq == 'a':
4298
4299            dDijDict[phfx+'D23'] = k*l
4300
4301        elif uniq == 'b':
4302
4303            dDijDict[phfx+'D13'] = h*l
4304
4305        elif uniq == 'c':
4306
4307            dDijDict[phfx+'D12'] = h*k
4308
4309            names.append()
4310
4311    else:
4312
4313        dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2,
4314
4315            phfx+'D12':h*k,phfx+'D13':h*l,phfx+'D23':k*l}
4316
4317    for item in dDijDict:
4318
4319        dDijDict[item] *= refl[4]**2*tand(refl[5]/2.0)
4320
4321    return dDijDict
4322
4323   
4324
4325def GetFprime(controlDict,Histograms):
4326
4327    FFtables = controlDict['FFtables']
4328
4329    if not FFtables:
4330
4331        return
4332
4333    histoList = Histograms.keys()
4334
4335    histoList.sort()
4336
4337    for histogram in histoList:
4338
4339        if 'PWDR' in histogram[:4]:
4340
4341            Histogram = Histograms[histogram]
4342
4343            hId = Histogram['hId']
4344
4345            hfx = ':%d:'%(hId)
4346
4347            keV = controlDict[hfx+'keV']
4348
4349            for El in FFtables:
4350
4351                Orbs = G2el.GetXsectionCoeff(El.split('+')[0].split('-')[0])
4352
4353                FP,FPP,Mu = G2el.FPcalc(Orbs, keV)
4354
4355                FFtables[El][hfx+'FP'] = FP
4356
4357                FFtables[El][hfx+'FPP'] = FPP               
4358
4359           
4360
4361def getPowderProfile(parmDict,x,varylist,Histogram,Phases,calcControls,pawleyLookup):
4362
4363   
4364
4365    def GetReflSIgGam(refl,wave,G,GB,hfx,phfx,calcControls,parmDict):
4366
4367        U = parmDict[hfx+'U']
4368
4369        V = parmDict[hfx+'V']
4370
4371        W = parmDict[hfx+'W']
4372
4373        X = parmDict[hfx+'X']
4374
4375        Y = parmDict[hfx+'Y']
4376
4377        tanPos = tand(refl[5]/2.0)
4378
4379        sig = U*tanPos**2+V*tanPos+W        #save peak sigma
4380
4381        sig = max(0.001,sig)
4382
4383        gam = X/cosd(refl[5]/2.0)+Y*tanPos+GetSampleGam(refl,wave,G,GB,phfx,calcControls,parmDict) #save peak gamma
4384
4385        gam = max(0.001,gam)
4386
4387        return sig,gam
4388
4389               
4390
4391    hId = Histogram['hId']
4392
4393    hfx = ':%d:'%(hId)
4394
4395    bakType = calcControls[hfx+'bakType']
4396
4397    yb = G2pwd.getBackground(hfx,parmDict,bakType,x)
4398
4399    yc = np.zeros_like(yb)
4400
4401       
4402
4403    if 'C' in calcControls[hfx+'histType']:   
4404
4405        shl = max(parmDict[hfx+'SH/L'],0.002)
4406
4407        Ka2 = False
4408
4409        if hfx+'Lam1' in parmDict.keys():
4410
4411            wave = parmDict[hfx+'Lam1']
4412
4413            Ka2 = True
4414
4415            lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
4416
4417            kRatio = parmDict[hfx+'I(L2)/I(L1)']
4418
4419        else:
4420
4421            wave = parmDict[hfx+'Lam']
4422
4423    else:
4424
4425        print 'TOF Undefined at present'
4426
4427        raise ValueError
4428
4429    for phase in Histogram['Reflection Lists']:
4430
4431        refList = Histogram['Reflection Lists'][phase]
4432
4433        Phase = Phases[phase]
4434
4435        pId = Phase['pId']
4436
4437        pfx = '%d::'%(pId)
4438
4439        phfx = '%d:%d:'%(pId,hId)
4440
4441        hfx = ':%d:'%(hId)
4442
4443        SGData = Phase['General']['SGData']
4444
4445        A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
4446
4447        G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
4448
4449        GA,GB = G2lat.Gmat2AB(G)    #Orthogonalization matricies
4450
4451        Vst = np.sqrt(nl.det(G))    #V*
4452
4453        if 'Pawley' not in Phase['General']['Type']:
4454
4455            refList = StructureFactor(refList,G,hfx,pfx,SGData,calcControls,parmDict)
4456
4457        for refl in refList:
4458
4459            if 'C' in calcControls[hfx+'histType']:
4460
4461                h,k,l = refl[:3]
4462
4463                refl[5] = GetReflPos(refl,wave,G,hfx,calcControls,parmDict)         #corrected reflection position
4464
4465                Lorenz = 1./(2.*sind(refl[5]/2.)**2*cosd(refl[5]/2.))           #Lorentz correction
4466
4467                refl[5] += GetHStrainShift(refl,SGData,phfx,parmDict)               #apply hydrostatic strain shift
4468
4469                refl[6:8] = GetReflSIgGam(refl,wave,G,GB,hfx,phfx,calcControls,parmDict)    #peak sig & gam
4470
4471                GetIntensityCorr(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)    #puts corrections in refl[13]
4472
4473                refl[13] *= Vst*Lorenz
4474
4475                if 'Pawley' in Phase['General']['Type']:
4476
4477                    try:
4478
4479                        refl[9] = abs(parmDict[pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])])
4480
4481                    except KeyError:
4482
4483#                        print ' ***Error %d,%d,%d missing from Pawley reflection list ***'%(h,k,l)
4484
4485                        continue
4486
4487                Wd,fmin,fmax = G2pwd.getWidths(refl[5],refl[6],refl[7],shl)
4488
4489                iBeg = np.searchsorted(x,refl[5]-fmin)
4490
4491                iFin = np.searchsorted(x,refl[5]+fmax)
4492
4493                if not iBeg+iFin:       #peak below low limit - skip peak
4494
4495                    continue
4496
4497                elif not iBeg-iFin:     #peak above high limit - done
4498
4499                    break
4500
4501                yc[iBeg:iFin] += refl[13]*refl[9]*G2pwd.getFCJVoigt3(refl[5],refl[6],refl[7],shl,x[iBeg:iFin])    #>90% of time spent here
4502
4503                if Ka2:
4504
4505                    pos2 = refl[5]+lamRatio*tand(refl[5]/2.0)       # + 360/pi * Dlam/lam * tan(th)
4506
4507                    Wd,fmin,fmax = G2pwd.getWidths(pos2,refl[6],refl[7],shl)
4508
4509                    iBeg = np.searchsorted(x,pos2-fmin)
4510
4511                    iFin = np.searchsorted(x,pos2+fmax)
4512
4513                    if not iBeg+iFin:       #peak below low limit - skip peak
4514
4515                        continue
4516
4517                    elif not iBeg-iFin:     #peak above high limit - done
4518
4519                        return yc,yb
4520
4521                    yc[iBeg:iFin] += refl[13]*refl[9]*kRatio*G2pwd.getFCJVoigt3(pos2,refl[6],refl[7],shl,x[iBeg:iFin])        #and here
4522
4523            elif 'T' in calcControls[hfx+'histType']:
4524
4525                print 'TOF Undefined at present'
4526
4527                raise Exception    #no TOF yet
4528
4529    return yc,yb
4530
4531   
4532
4533def GetFobsSq(Histograms,Phases,parmDict,calcControls):
4534
4535    histoList = Histograms.keys()
4536
4537    histoList.sort()
4538
4539    for histogram in histoList:
4540
4541        if 'PWDR' in histogram[:4]:
4542
4543            Histogram = Histograms[histogram]
4544
4545            hId = Histogram['hId']
4546
4547            hfx = ':%d:'%(hId)
4548
4549            Limits = calcControls[hfx+'Limits']
4550
4551            shl = max(parmDict[hfx+'SH/L'],0.002)
4552
4553            Ka2 = False
4554
4555            kRatio = 0.0
4556
4557            if hfx+'Lam1' in parmDict.keys():
4558
4559                Ka2 = True
4560
4561                lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
4562
4563                kRatio = parmDict[hfx+'I(L2)/I(L1)']
4564
4565            x,y,w,yc,yb,yd = Histogram['Data']
4566
4567            ymb = np.array(y-yb)
4568
4569            ycmb = np.array(yc-yb)
4570
4571            ratio = np.where(ycmb!=0.,ymb/ycmb,0.0)           
4572
4573            xB = np.searchsorted(x,Limits[0])
4574
4575            xF = np.searchsorted(x,Limits[1])
4576
4577            refLists = Histogram['Reflection Lists']
4578
4579            for phase in refLists:
4580
4581                Phase = Phases[phase]
4582
4583                pId = Phase['pId']
4584
4585                phfx = '%d:%d:'%(pId,hId)
4586
4587                refList = refLists[phase]
4588
4589                sumFo = 0.0
4590
4591                sumdF = 0.0
4592
4593                sumFosq = 0.0
4594
4595                sumdFsq = 0.0
4596
4597                for refl in refList:
4598
4599                    if 'C' in calcControls[hfx+'histType']:
4600
4601                        yp = np.zeros_like(yb)
4602
4603                        Wd,fmin,fmax = G2pwd.getWidths(refl[5],refl[6],refl[7],shl)
4604
4605                        iBeg = np.searchsorted(x[xB:xF],refl[5]-fmin)
4606
4607                        iFin = np.searchsorted(x[xB:xF],refl[5]+fmax)
4608
4609                        iFin2 = iFin
4610
4611                        yp[iBeg:iFin] = refl[13]*refl[9]*G2pwd.getFCJVoigt3(refl[5],refl[6],refl[7],shl,x[iBeg:iFin])    #>90% of time spent here                           
4612
4613                        if Ka2:
4614
4615                            pos2 = refl[5]+lamRatio*tand(refl[5]/2.0)       # + 360/pi * Dlam/lam * tan(th)
4616
4617                            Wd,fmin,fmax = G2pwd.getWidths(pos2,refl[6],refl[7],shl)
4618
4619                            iBeg2 = np.searchsorted(x,pos2-fmin)
4620
4621                            iFin2 = np.searchsorted(x,pos2+fmax)
4622
4623                            yp[iBeg2:iFin2] += refl[13]*refl[9]*kRatio*G2pwd.getFCJVoigt3(pos2,refl[6],refl[7],shl,x[iBeg2:iFin2])        #and here
4624
4625                        refl[8] = np.sum(np.where(ratio[iBeg:iFin2]>0.,yp[iBeg:iFin2]*ratio[iBeg:iFin2]/(refl[13]*(1.+kRatio)),0.0))
4626
4627                    elif 'T' in calcControls[hfx+'histType']:
4628
4629                        print 'TOF Undefined at present'
4630
4631                        raise Exception    #no TOF yet
4632
4633                    Fo = np.sqrt(np.abs(refl[8]))
4634
4635                    Fc = np.sqrt(np.abs(refl[9]))
4636
4637                    sumFo += Fo
4638
4639                    sumFosq += refl[8]**2
4640
4641                    sumdF += np.abs(Fo-Fc)
4642
4643                    sumdFsq += (refl[8]-refl[9])**2
4644
4645                Histogram[phfx+'Rf'] = min(100.,(sumdF/sumFo)*100.)
4646
4647                Histogram[phfx+'Rf^2'] = min(100.,np.sqrt(sumdFsq/sumFosq)*100.)
4648
4649                Histogram[phfx+'Nref'] = len(refList)
4650
4651               
4652
4653def getPowderProfileDerv(parmDict,x,varylist,Histogram,Phases,calcControls,pawleyLookup):
4654
4655   
4656
4657    def cellVaryDerv(pfx,SGData,dpdA): 
4658
4659        if SGData['SGLaue'] in ['-1',]:
4660
4661            return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],
4662
4663                [pfx+'A3',dpdA[3]],[pfx+'A4',dpdA[4]],[pfx+'A5',dpdA[5]]]
4664
4665        elif SGData['SGLaue'] in ['2/m',]:
4666
4667            if SGData['SGUniq'] == 'a':
4668
4669                return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A3',dpdA[3]]]
4670
4671            elif SGData['SGUniq'] == 'b':
4672
4673                return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A4',dpdA[4]]]
4674
4675            else:
4676
4677                return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A5',dpdA[5]]]
4678
4679        elif SGData['SGLaue'] in ['mmm',]:
4680
4681            return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]]]
4682
4683        elif SGData['SGLaue'] in ['4/m','4/mmm']:
4684
4685#            return [[pfx+'A0',dpdA[0]+dpdA[1]],[pfx+'A2',dpdA[2]]]
4686
4687            return [[pfx+'A0',dpdA[0]],[pfx+'A2',dpdA[2]]]
4688
4689        elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
4690
4691#            return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[3]],[pfx+'A2',dpdA[2]]]
4692
4693            return [[pfx+'A0',dpdA[0]],[pfx+'A2',dpdA[2]]]
4694
4695        elif SGData['SGLaue'] in ['3R', '3mR']:
4696
4697            return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[2]],[pfx+'A3',dpdA[3]+dpdA[4]+dpdA[5]]]                       
4698
4699        elif SGData['SGLaue'] in ['m3m','m3']:
4700
4701#            return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[2]]]
4702
4703            return [[pfx+'A0',dpdA[0]]]
4704
4705    # create a list of dependent variables and set up a dictionary to hold their derivatives
4706
4707    dependentVars = G2mv.GetDependentVars()
4708
4709    depDerivDict = {}
4710
4711    for j in dependentVars:
4712
4713        depDerivDict[j] = np.zeros(shape=(len(x)))
4714
4715    #print 'dependent vars',dependentVars
4716
4717    lenX = len(x)               
4718
4719    hId = Histogram['hId']
4720
4721    hfx = ':%d:'%(hId)
4722
4723    bakType = calcControls[hfx+'bakType']
4724
4725    dMdv = np.zeros(shape=(len(varylist),len(x)))
4726
4727    dMdb,dMddb = G2pwd.getBackgroundDerv(hfx,parmDict,bakType,x)
4728
4729    if hfx+'Back:0' in varylist: # for now assume that Back:x vars to not appear in constraints
4730
4731        bBpos =varylist.index(hfx+'Back:0')
4732
4733        dMdv[bBpos:bBpos+len(dMdb)] = dMdb
4734
4735    names = [hfx+'DebyeA',hfx+'DebyeR',hfx+'DebyeU']
4736
4737    for name in varylist:
4738
4739        if 'Debye' in name:
4740
4741            id = int(name.split(':')[-1])
4742
4743            parm = name[:int(name.rindex(':'))]
4744
4745            ip = names.index(parm)
4746
4747            dMdv[varylist.index(name)] = dMddb[3*id+ip]
4748
4749    if 'C' in calcControls[hfx+'histType']:   
4750
4751        dx = x[1]-x[0]
4752
4753        shl = max(parmDict[hfx+'SH/L'],0.002)
4754
4755        Ka2 = False
4756
4757        if hfx+'Lam1' in parmDict.keys():
4758
4759            wave = parmDict[hfx+'Lam1']
4760
4761            Ka2 = True
4762
4763            lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1'])
4764
4765            kRatio = parmDict[hfx+'I(L2)/I(L1)']
4766
4767        else:
4768
4769            wave = parmDict[hfx+'Lam']
4770
4771    else:
4772
4773        print 'TOF Undefined at present'
4774
4775        raise ValueError
4776
4777    for phase in Histogram['Reflection Lists']:
4778
4779        refList = Histogram['Reflection Lists'][phase]
4780
4781        Phase = Phases[phase]
4782
4783        SGData = Phase['General']['SGData']
4784
4785        pId = Phase['pId']
4786
4787        pfx = '%d::'%(pId)
4788
4789        phfx = '%d:%d:'%(pId,hId)
4790
4791        A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
4792
4793        G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
4794
4795        GA,GB = G2lat.Gmat2AB(G)    #Orthogonalization matricies
4796
4797        if 'Pawley' not in Phase['General']['Type']:
4798
4799            dFdvDict = StructureFactorDerv(refList,G,hfx,pfx,SGData,calcControls,parmDict)
4800
4801        for iref,refl in enumerate(refList):
4802
4803            if 'C' in calcControls[hfx+'histType']:        #CW powder
4804
4805                h,k,l = refl[:3]
4806
4807                dIdsh,dIdsp,dIdpola,dIdPO,dFdODF,dFdSA = GetIntensityDerv(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
4808
4809                if 'Pawley' in Phase['General']['Type']:
4810
4811                    try:
4812
4813                        refl[9] = abs(parmDict[pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])])
4814
4815                    except KeyError:
4816
4817#                        print ' ***Error %d,%d,%d missing from Pawley reflection list ***'%(h,k,l)
4818
4819                        continue
4820
4821                Wd,fmin,fmax = G2pwd.getWidths(refl[5],refl[6],refl[7],shl)
4822
4823                iBeg = np.searchsorted(x,refl[5]-fmin)
4824
4825                iFin = np.searchsorted(x,refl[5]+fmax)
4826
4827                if not iBeg+iFin:       #peak below low limit - skip peak
4828
4829                    continue
4830
4831                elif not iBeg-iFin:     #peak above high limit - done
4832
4833                    break
4834
4835                pos = refl[5]
4836
4837                tanth = tand(pos/2.0)
4838
4839                costh = cosd(pos/2.0)
4840
4841                dMdpk = np.zeros(shape=(6,len(x)))
4842
4843                dMdipk = G2pwd.getdFCJVoigt3(refl[5],refl[6],refl[7],shl,x[iBeg:iFin])
4844
4845                for i in range(1,5):
4846
4847                    dMdpk[i][iBeg:iFin] += 100.*dx*refl[13]*refl[9]*dMdipk[i]
4848
4849                dMdpk[0][iBeg:iFin] += 100.*dx*refl[13]*refl[9]*dMdipk[0]
4850
4851                dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4]}
4852
4853                if Ka2:
4854
4855                    pos2 = refl[5]+lamRatio*tanth       # + 360/pi * Dlam/lam * tan(th)
4856
4857                    kdelt = int((pos2-refl[5])/dx)               
4858
4859                    iBeg2 = min(lenX,iBeg+kdelt)
4860
4861                    iFin2 = min(lenX,iFin+kdelt)
4862
4863                    if iBeg2-iFin2:
4864
4865                        dMdipk2 = G2pwd.getdFCJVoigt3(pos2,refl[6],refl[7],shl,x[iBeg2:iFin2])
4866
4867                        for i in range(1,5):
4868
4869                            dMdpk[i][iBeg2:iFin2] += 100.*dx*refl[13]*refl[9]*kRatio*dMdipk2[i]
4870
4871                        dMdpk[0][iBeg2:iFin2] += 100.*dx*refl[13]*refl[9]*kRatio*dMdipk2[0]
4872
4873                        dMdpk[5][iBeg2:iFin2] += 100.*dx*refl[13]*dMdipk2[0]
4874
4875                        dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':dMdpk[5]*refl[9]}
4876
4877                if 'Pawley' in Phase['General']['Type']:
4878
4879                    try:
4880
4881                        idx = varylist.index(pfx+'PWLref:'+str(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)]))
4882
4883                        dMdv[idx] = dervDict['int']/refl[9]
4884
4885                        # Assuming Pawley variables not in constraints
4886
4887                    except ValueError:
4888
4889                        pass
4890
4891                dpdA,dpdw,dpdZ,dpdSh,dpdTr,dpdX,dpdY = GetReflPosDerv(refl,wave,A,hfx,calcControls,parmDict)
4892
4893                names = {hfx+'Scale':[dIdsh,'int'],hfx+'Polariz.':[dIdpola,'int'],phfx+'Scale':[dIdsp,'int'],
4894
4895                    hfx+'U':[tanth**2,'sig'],hfx+'V':[tanth,'sig'],hfx+'W':[1.0,'sig'],
4896
4897                    hfx+'X':[1.0/costh,'gam'],hfx+'Y':[tanth,'gam'],hfx+'SH/L':[1.0,'shl'],
4898
4899                    hfx+'I(L2)/I(L1)':[1.0,'L1/L2'],hfx+'Zero':[dpdZ,'pos'],hfx+'Lam':[dpdw,'pos'],
4900
4901                    hfx+'Shift':[dpdSh,'pos'],hfx+'Transparency':[dpdTr,'pos'],hfx+'DisplaceX':[dpdX,'pos'],
4902
4903                    hfx+'DisplaceY':[dpdY,'pos'],}
4904
4905                for name in names:
4906
4907                    item = names[name]
4908
4909                    if name in varylist:
4910
4911                        dMdv[varylist.index(name)] += item[0]*dervDict[item[1]]
4912
4913                    elif name in dependentVars:
4914
4915                        depDerivDict[name] += item[0]*dervDict[item[1]]
4916
4917
4918
4919                for iPO in dIdPO:
4920
4921                    if iPO in varylist:
4922
4923                        dMdv[varylist.index(iPO)] += dIdPO[iPO]*dervDict['int']
4924
4925                    elif iPO in dependentVars:
4926
4927                        depDerivDict[iPO] = dIdPO[iPO]*dervDict['int']
4928
4929
4930
4931                for i,name in enumerate(['omega','chi','phi']):
4932
4933                    aname = pfx+'SH '+name
4934
4935                    if aname in varylist:
4936
4937                        dMdv[varylist.index(aname)] += dFdSA[i]*dervDict['int']
4938
4939                    elif aname in dependentVars:
4940
4941                        depDerivDict[aname] += dFdSA[i]*dervDict['int']
4942
4943                for iSH in dFdODF:
4944
4945                    if iSH in varylist:
4946
4947                        dMdv[varylist.index(iSH)] += dFdODF[iSH]*dervDict['int']
4948
4949                    elif iSH in dependentVars:
4950
4951                        depDerivDict[iSH] += dFdODF[iSH]*dervDict['int']
4952
4953                cellDervNames = cellVaryDerv(pfx,SGData,dpdA)
4954
4955                for name,dpdA in cellDervNames:
4956
4957                    if name in varylist:
4958
4959                        dMdv[varylist.index(name)] += dpdA*dervDict['pos']
4960
4961                    elif name in dependentVars:
4962
4963                        depDerivDict[name] += dpdA*dervDict['pos']
4964
4965                dDijDict = GetHStrainShiftDerv(refl,SGData,phfx)
4966
4967                for name in dDijDict:
4968
4969                    if name in varylist:
4970
4971                        dMdv[varylist.index(name)] += dDijDict[name]*dervDict['pos']
4972
4973                    elif name in dependentVars:
4974
4975                        depDerivDict[name] += dDijDict[name]*dervDict['pos']
4976
4977                gamDict = GetSampleGamDerv(refl,wave,G,GB,phfx,calcControls,parmDict)
4978
4979                for name in gamDict:
4980
4981                    if name in varylist:
4982
4983                        dMdv[varylist.index(name)] += gamDict[name]*dervDict['gam']
4984
4985                    elif name in dependentVars:
4986
4987                        depDerivDict[name] += gamDict[name]*dervDict['gam']
4988
4989                                               
4990
4991            elif 'T' in calcControls[hfx+'histType']:
4992
4993                print 'TOF Undefined at present'
4994
4995                raise Exception    #no TOF yet
4996
4997            #do atom derivatives -  for F,X & U so far             
4998
4999            corr = dervDict['int']/refl[9]
5000
5001            for name in varylist+dependentVars:
5002
5003                try:
5004
5005                    aname = name.split(pfx)[1][:2]
5006
5007                    if aname not in ['Af','dA','AU']: continue # skip anything not an atom param
5008
5009                except IndexError:
5010
5011                    continue
5012
5013                if name in varylist:
5014
5015                    dMdv[varylist.index(name)] += dFdvDict[name][iref]*corr
5016
5017                elif name in dependentVars:
5018
5019                    depDerivDict[name] += dFdvDict[name][iref]*corr
5020
5021    # now process derivatives in constraints
5022
5023    G2mv.Dict2Deriv(varylist,depDerivDict,dMdv)
5024
5025    return dMdv
5026
5027
5028
5029def dervRefine(values,HistoPhases,parmdict,varylist,calcControls,pawleyLookup,dlg):
5030
5031    parmdict.update(zip(varylist,values))
5032
5033    G2mv.Dict2Map(parmdict,varylist)
5034
5035    Histograms,Phases = HistoPhases
5036
5037    nvar = len(varylist)
5038
5039    dMdv = np.empty(0)
5040
5041    histoList = Histograms.keys()
5042
5043    histoList.sort()
5044
5045    for histogram in histoList:
5046
5047        if 'PWDR' in histogram[:4]:
5048
5049            Histogram = Histograms[histogram]
5050
5051            hId = Histogram['hId']
5052
5053            hfx = ':%d:'%(hId)
5054
5055            Limits = calcControls[hfx+'Limits']
5056
5057            x,y,w,yc,yb,yd = Histogram['Data']
5058
5059            xB = np.searchsorted(x,Limits[0])
5060
5061            xF = np.searchsorted(x,Limits[1])
5062
5063            dMdvh = np.sqrt(w[xB:xF])*getPowderProfileDerv(parmdict,x[xB:xF],
5064
5065                varylist,Histogram,Phases,calcControls,pawleyLookup)
5066
5067            if len(dMdv):
5068
5069                dMdv = np.concatenate((dMdv.T,dMdvh.T)).T
5070
5071            else:
5072
5073                dMdv = dMdvh
5074
5075    return dMdv
5076
5077
5078
5079def ComputePowderHessian(args):
5080
5081    Histogram,parmdict,varylist,Phases,calcControls,pawleyLookup = args
5082
5083    hId = Histogram['hId']
5084
5085    hfx = ':%d:'%(hId)
5086
5087    Limits = calcControls[hfx+'Limits']
5088
5089    x,y,w,yc,yb,yd = Histogram['Data']
5090
5091    dy = y-yc
5092
5093    xB = np.searchsorted(x,Limits[0])
5094
5095    xF = np.searchsorted(x,Limits[1])
5096
5097    dMdvh = np.sqrt(w[xB:xF])*getPowderProfileDerv(
5098
5099        parmdict,x[xB:xF],
5100
5101        varylist,Histogram,Phases,calcControls,pawleyLookup)
5102
5103    Vec = np.sum(dMdvh*np.sqrt(w[xB:xF])*dy[xB:xF],axis=1)
5104
5105    Hess = np.inner(dMdvh,dMdvh)
5106
5107    return Vec,Hess
5108
5109
5110
5111def HessRefine(values,HistoPhases,parmdict,varylist,calcControls,pawleyLookup,dlg):
5112
5113    parmdict.update(zip(varylist,values))
5114
5115    G2mv.Dict2Map(parmdict,varylist)
5116
5117    Histograms,Phases = HistoPhases
5118
5119    nvar = len(varylist)
5120
5121    HessSum = np.zeros((nvar,nvar))
5122
5123    VecSum = np.zeros(nvar)
5124
5125    histoList = Histograms.keys()
5126
5127    histoList.sort()
5128
5129    argList = []
5130
5131    for histogram in histoList:
5132
5133        if 'PWDR' in histogram[:4]:
5134
5135            Histogram = Histograms[histogram]
5136
5137            argList.append([
5138
5139                Histogram,parmdict,varylist,Phases,
5140
5141                calcControls,pawleyLookup])
5142
5143    if MaxProcess > 1: 
5144
5145        mpPool = mp.Pool(processes=MaxProcess)
5146
5147        #results = mpPool.map(ComputePowderHessian,argList)
5148
5149        #for Vec,Hess in results:
5150
5151        for Vec,Hess in mpPool.imap_unordered(ComputePowderHessian,argList):
5152
5153            VecSum += Vec
5154
5155            HessSum += Hess
5156
5157    else:
5158
5159        for arg in argList:
5160
5161            Vec,Hess = ComputePowderHessian(arg)
5162
5163            VecSum += Vec
5164
5165            HessSum += Hess
5166
5167    return VecSum,HessSum
5168
5169
5170
5171def ComputePowderProfile(args):
5172
5173    Histogram,parmdict,varylist,Phases,calcControls,pawleyLookup = args
5174
5175    hId = Histogram['hId']
5176
5177    print 'ComputePowderProfile',hId
5178
5179    hfx = ':%d:'%(hId)
5180
5181    x,y,w,yc,yb,yd = Histogram['Data']
5182
5183    Limits = calcControls[hfx+'Limits']
5184
5185    xB = np.searchsorted(x,Limits[0])
5186
5187    xF = np.searchsorted(x,Limits[1])
5188
5189    yc,yb = getPowderProfile(parmdict,x[xB:xF],
5190
5191                            varylist,Histogram,Phases,calcControls,
5192
5193                            pawleyLookup)
5194
5195    return xB,xF,yc,yb,Histogram['Reflection Lists']
5196
5197
5198
5199def errRefine(values,HistoPhases,parmdict,varylist,calcControls,pawleyLookup,dlg):       
5200
5201    parmdict.update(zip(varylist,values))
5202
5203    Values2Dict(parmdict, varylist, values)
5204
5205    G2mv.Dict2Map(parmdict,varylist)
5206
5207    Histograms,Phases = HistoPhases
5208
5209    M = np.empty(0)
5210
5211    sumwYo = 0
5212
5213    Nobs = 0
5214
5215    histoList = Histograms.keys()
5216
5217    histoList.sort()
5218
5219    argList = []
5220
5221    for histogram in histoList:
5222
5223        if 'PWDR' in histogram[:4]:
5224
5225            Histogram = Histograms[histogram]
5226
5227            argList.append(
5228
5229                [Histogram,parmdict,varylist,Phases,calcControls,pawleyLookup]
5230
5231                )
5232
5233    if MaxProcess > 1: 
5234
5235        mpPool = mp.Pool(processes=MaxProcess)
5236
5237        #results = mpPool.map(ComputePowderProfile,argList)
5238
5239        #for arg,res in zip(argList,results):
5240
5241        #results = mpPool.map(ComputePowderProfile,argList)
5242
5243        #for i,res in enumerate(results):
5244
5245        for i,res in enumerate(
5246
5247            mpPool.imap_unordered(ComputePowderProfile,argList)
5248
5249            ):
5250
5251            print 'process',i
5252
5253            xB,xF,ycSect,ybSect,RL = res
5254
5255            Histogram = argList[i][0]
5256
5257            Histogram['Reflection Lists'] = RL
5258
5259            x,y,w,yc,yb,yd = Histogram['Data']
5260
5261            yc *= 0.0                           #zero full calcd profiles
5262
5263            yb *= 0.0
5264
5265            yd *= 0.0
5266
5267            Histogram['Nobs'] = xF-xB
5268
5269            Nobs += Histogram['Nobs']
5270
5271            Histogram['sumwYo'] = np.sum(w[xB:xF]*y[xB:xF]**2)
5272
5273            sumwYo += Histogram['sumwYo']
5274
5275           
5276
5277            yc[xB:xF] = ycSect
5278
5279            yb[xB:xF] = ybSect
5280
5281            yc[xB:xF] += yb[xB:xF]
5282
5283            yd[xB:xF] = y[xB:xF]-yc[xB:xF]
5284
5285            Histogram['sumwYd'] = np.sum(np.sqrt(w[xB:xF])*(yd[xB:xF]))
5286
5287            wdy = -np.sqrt(w[xB:xF])*(yd[xB:xF])
5288
5289            Histogram['wRp'] = min(100.,np.sqrt(np.sum(wdy**2)/Histogram['sumwYo'])*100.)
5290
5291            M = np.concatenate((M,wdy))
5292
5293    else:
5294
5295        for arg in argList:
5296
5297            xB,xF,ycSect,ybSect,RL = ComputePowderProfile(arg)
5298
5299            Histogram = arg[0]
5300
5301            x,y,w,yc,yb,yd = Histogram['Data']
5302
5303            yc *= 0.0                           #zero full calcd profiles
5304
5305            yb *= 0.0
5306
5307            yd *= 0.0
5308
5309            Histogram['Nobs'] = xF-xB
5310
5311            Nobs += Histogram['Nobs']
5312
5313            Histogram['sumwYo'] = np.sum(w[xB:xF]*y[xB:xF]**2)
5314
5315            sumwYo += Histogram['sumwYo']
5316
5317           
5318
5319            yc[xB:xF] = ycSect
5320
5321            yb[xB:xF] = ybSect
5322
5323            yc[xB:xF] += yb[xB:xF]
5324
5325            yd[xB:xF] = y[xB:xF]-yc[xB:xF]
5326
5327            Histogram['sumwYd'] = np.sum(np.sqrt(w[xB:xF])*(yd[xB:xF]))
5328
5329            wdy = -np.sqrt(w[xB:xF])*(yd[xB:xF])
5330
5331            Histogram['wRp'] = min(100.,np.sqrt(np.sum(wdy**2)/Histogram['sumwYo'])*100.)
5332
5333            if dlg:
5334
5335                dlg.Update(Histogram['wRp'],newmsg='For histogram %d Rwp=%8.3f%s'%(hId,Histogram['wRp'],'%'))[0]
5336
5337            M = np.concatenate((M,wdy))
5338
5339
5340
5341    Histograms['sumwYo'] = sumwYo
5342
5343    Histograms['Nobs'] = Nobs
5344
5345    Rwp = min(100.,np.sqrt(np.sum(M**2)/sumwYo)*100.)
5346
5347    if dlg:
5348
5349        GoOn = dlg.Update(Rwp,newmsg='%s%8.3f%s'%('Powder profile Rwp =',Rwp,'%'))[0]
5350
5351        if not GoOn:
5352
5353            parmDict['saved values'] = values
5354
5355            raise Exception         #Abort!!
5356
5357    return M
5358
5359   
5360
5361                   
5362
5363def Refine(GPXfile,dlg):
5364
5365    import pytexture as ptx
5366
5367    ptx.pyqlmninit()            #initialize fortran arrays for spherical harmonics
5368
5369   
5370
5371    ShowBanner()
5372
5373    varyList = []
5374
5375    parmDict = {}
5376
5377    calcControls = {}
5378
5379    G2mv.InitVars()   
5380
5381    Controls = GetControls(GPXfile)
5382
5383    ShowControls(Controls)           
5384
5385    constrDict,constrFlag,fixedList = GetConstraints(GPXfile)
5386
5387    Histograms,Phases = GetUsedHistogramsAndPhases(GPXfile)
5388
5389    if not Phases:
5390
5391        print ' *** ERROR - you have no histograms to refine! ***'
5392
5393        print ' *** Refine aborted ***'
5394
5395        raise Exception
5396
5397    if not Histograms:
5398
5399        print ' *** ERROR - you have no data to refine with! ***'
5400
5401        print ' *** Refine aborted ***'
5402
5403        raise Exception       
5404
5405    Natoms,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables = GetPhaseData(Phases)
5406
5407    calcControls['Natoms'] = Natoms
5408
5409    calcControls['FFtables'] = FFtables
5410
5411    calcControls['BLtables'] = BLtables
5412
5413    hapVary,hapDict,controlDict = GetHistogramPhaseData(Phases,Histograms)
5414
5415    calcControls.update(controlDict)
5416
5417    histVary,histDict,controlDict = GetHistogramData(Histograms)
5418
5419    calcControls.update(controlDict)
5420
5421    varyList = phaseVary+hapVary+histVary
5422
5423    parmDict.update(phaseDict)
5424
5425    parmDict.update(hapDict)
5426
5427    parmDict.update(histDict)
5428
5429    GetFprime(calcControls,Histograms)
5430
5431    # do constraint processing
5432
5433    try:
5434
5435        groups,parmlist = G2mv.GroupConstraints(constrDict)
5436
5437        G2mv.GenerateConstraints(groups,parmlist,varyList,constrDict,constrFlag,fixedList)
5438
5439    except:
5440
5441        print ' *** ERROR - your constraints are internally inconsistent ***'
5442
5443        raise Exception(' *** Refine aborted ***')
5444
5445    # check to see which generated parameters are fully varied
5446
5447    msg = G2mv.SetVaryFlags(varyList)
5448
5449    if msg:
5450
5451        print ' *** ERROR - you have not set the refine flags for constraints consistently! ***'
5452
5453        print msg
5454
5455        raise Exception(' *** Refine aborted ***')
5456
5457    G2mv.Map2Dict(parmDict,varyList)
5458
5459#    print G2mv.VarRemapShow(varyList)
5460
5461
5462
5463    while True:
5464
5465        begin = time.time()
5466
5467        values =  np.array(Dict2Values(parmDict, varyList))
5468
5469        Ftol = Controls['min dM/M']
5470
5471        Factor = Controls['shift factor']
5472
5473        maxCyc = Controls['max cyc']
5474
5475        if 'Jacobian' in Controls['deriv type']:           
5476
5477            result = so.leastsq(errRefine,values,Dfun=dervRefine,full_output=True,
5478
5479                ftol=Ftol,col_deriv=True,factor=Factor,
5480
5481                args=([Histograms,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg))
5482
5483            ncyc = int(result[2]['nfev']/2)
5484
5485        elif 'Hessian' in Controls['deriv type']:
5486
5487            result = G2mth.HessianLSQ(errRefine,values,Hess=HessRefine,ftol=Ftol,maxcyc=maxCyc,
5488
5489                args=([Histograms,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg))
5490
5491            ncyc = result[2]['num cyc']+1                           
5492
5493        else:           #'numeric'
5494
5495            result = so.leastsq(errRefine,values,full_output=True,ftol=Ftol,epsfcn=1.e-8,factor=Factor,
5496
5497                args=([Histograms,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg))
5498
5499            ncyc = int(result[2]['nfev']/len(varyList))
5500
5501#        table = dict(zip(varyList,zip(values,result[0],(result[0]-values))))
5502
5503#        for item in table: print item,table[item]               #useful debug - are things shifting?
5504
5505        runtime = time.time()-begin
5506
5507        chisq = np.sum(result[2]['fvec']**2)
5508
5509        Values2Dict(parmDict, varyList, result[0])
5510
5511        G2mv.Dict2Map(parmDict,varyList)
5512
5513       
5514
5515        Rwp = np.sqrt(chisq/Histograms['sumwYo'])*100.      #to %
5516
5517        GOF = chisq/(Histograms['Nobs']-len(varyList))
5518
5519        print '\n Refinement results:'
5520
5521        print 135*'-'
5522
5523        print ' Number of function calls:',result[2]['nfev'],' Number of observations: ',Histograms['Nobs'],' Number of parameters: ',len(varyList)
5524
5525        print ' Refinement time = %8.3fs, %8.3fs/cycle, for %d cycles'%(runtime,runtime/ncyc,ncyc)
5526
5527        print ' wRp = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f'%(Rwp,chisq,GOF)
5528
5529        try:
5530
5531            covMatrix = result[1]*GOF
5532
5533            sig = np.sqrt(np.diag(covMatrix))
5534
5535            if np.any(np.isnan(sig)):
5536
5537                print '*** Least squares aborted - some invalid esds possible ***'
5538
5539#            table = dict(zip(varyList,zip(values,result[0],(result[0]-values)/sig)))
5540
5541#            for item in table: print item,table[item]               #useful debug - are things shifting?
5542
5543            break                   #refinement succeeded - finish up!
5544
5545        except TypeError:          #result[1] is None on singular matrix
5546
5547            print '**** Refinement failed - singular matrix ****'
5548
5549            if 'Hessian' in Controls['deriv type']:
5550
5551                for i in result[2]['psing'].reverse():
5552
5553                        print 'Removing parameter: ',varyList[i]
5554
5555                        del(varyList[i])                   
5556
5557            else:
5558
5559                Ipvt = result[2]['ipvt']
5560
5561                for i,ipvt in enumerate(Ipvt):
5562
5563                    if not np.sum(result[2]['fjac'],axis=1)[i]:
5564
5565                        print 'Removing parameter: ',varyList[ipvt-1]
5566
5567                        del(varyList[ipvt-1])
5568
5569                        break
5570
5571
5572
5573#    print 'dependentParmList: ',G2mv.dependentParmList
5574
5575#    print 'arrayList: ',G2mv.arrayList
5576
5577#    print 'invarrayList: ',G2mv.invarrayList
5578
5579#    print 'indParmList: ',G2mv.indParmList
5580
5581#    print 'fixedDict: ',G2mv.fixedDict
5582
5583#    print 'test1'
5584
5585    GetFobsSq(Histograms,Phases,parmDict,calcControls)
5586
5587#    print 'test2'
5588
5589    sigDict = dict(zip(varyList,sig))
5590
5591    newCellDict = GetNewCellParms(parmDict,varyList)
5592
5593    newAtomDict = ApplyXYZshifts(parmDict,varyList)
5594
5595    covData = {'variables':result[0],'varyList':varyList,'sig':sig,
5596
5597        'covMatrix':covMatrix,'title':GPXfile,'newAtomDict':newAtomDict,'newCellDict':newCellDict}
5598
5599    # add the uncertainties into the esd dictionary (sigDict)
5600
5601    sigDict.update(G2mv.ComputeDepESD(covMatrix,varyList,parmDict))
5602
5603    SetPhaseData(parmDict,sigDict,Phases,covData)
5604
5605    SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms)
5606
5607    SetHistogramData(parmDict,sigDict,Histograms)
5608
5609    G2mv.PrintIndependentVars(parmDict,varyList,sigDict)
5610
5611    SetUsedHistogramsAndPhases(GPXfile,Histograms,Phases,covData)
5612
5613   
5614
5615#for testing purposes!!!
5616
5617#    import cPickle
5618
5619#    file = open('structTestdata.dat','wb')
5620
5621#    cPickle.dump(parmDict,file,1)
5622
5623#    cPickle.dump(varyList,file,1)
5624
5625#    for histogram in Histograms:
5626
5627#        if 'PWDR' in histogram[:4]:
5628
5629#            Histogram = Histograms[histogram]
5630
5631#    cPickle.dump(Histogram,file,1)
5632
5633#    cPickle.dump(Phases,file,1)
5634
5635#    cPickle.dump(calcControls,file,1)
5636
5637#    cPickle.dump(pawleyLookup,file,1)
5638
5639#    file.close()
5640
5641
5642
5643    if dlg:
5644
5645        return Rwp
5646
5647
5648
5649def SeqRefine(GPXfile,dlg):
5650
5651    import pytexture as ptx
5652
5653    ptx.pyqlmninit()            #initialize fortran arrays for spherical harmonics
5654
5655   
5656
5657    ShowBanner()
5658
5659    print ' Sequential Refinement'
5660
5661    G2mv.InitVars()   
5662
5663    Controls = GetControls(GPXfile)
5664
5665    ShowControls(Controls)           
5666
5667    constrDict,constrFlag,fixedList = GetConstraints(GPXfile)
5668
5669    Histograms,Phases = GetUsedHistogramsAndPhases(GPXfile)
5670
5671    if not Phases:
5672
5673        print ' *** ERROR - you have no histograms to refine! ***'
5674
5675        print ' *** Refine aborted ***'
5676
5677        raise Exception
5678
5679    if not Histograms:
5680
5681        print ' *** ERROR - you have no data to refine with! ***'
5682
5683        print ' *** Refine aborted ***'
5684
5685        raise Exception
5686
5687    Natoms,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables = GetPhaseData(Phases,False)
5688
5689    if 'Seq Data' in Controls:
5690
5691        histNames = Controls['Seq Data']
5692
5693    else:
5694
5695        histNames = GetHistogramNames(GPXfile,['PWDR',])
5696
5697    if 'Reverse Seq' in Controls:
5698
5699        if Controls['Reverse Seq']:
5700
5701            histNames.reverse()
5702
5703    SeqResult = {'histNames':histNames}
5704
5705    makeBack = True
5706
5707    for ihst,histogram in enumerate(histNames):
5708
5709        ifPrint = False
5710
5711        if dlg:
5712
5713            dlg.SetTitle('Residual for histogram '+str(ihst))
5714
5715        calcControls = {}
5716
5717        calcControls['Natoms'] = Natoms
5718
5719        calcControls['FFtables'] = FFtables
5720
5721        calcControls['BLtables'] = BLtables
5722
5723        varyList = []
5724
5725        parmDict = {}
5726
5727        Histo = {histogram:Histograms[histogram],}
5728
5729        hapVary,hapDict,controlDict = GetHistogramPhaseData(Phases,Histo,False)
5730
5731        calcControls.update(controlDict)
5732
5733        histVary,histDict,controlDict = GetHistogramData(Histo,False)
5734
5735        calcControls.update(controlDict)
5736
5737        varyList = phaseVary+hapVary+histVary
5738
5739        if not ihst:
5740
5741            saveVaryList = varyList[:]
5742
5743            for i,item in enumerate(saveVaryList):
5744
5745                items = item.split(':')
5746
5747                if items[1]:
5748
5749                    items[1] = ''
5750
5751                item = ':'.join(items)
5752
5753                saveVaryList[i] = item
5754
5755            SeqResult['varyList'] = saveVaryList
5756
5757        else:
5758
5759            newVaryList = varyList[:]
5760
5761            for i,item in enumerate(newVaryList):
5762
5763                items = item.split(':')
5764
5765                if items[1]:
5766
5767                    items[1] = ''
5768
5769                item = ':'.join(items)
5770
5771                newVaryList[i] = item
5772
5773            if newVaryList != SeqResult['varyList']:
5774
5775                print newVaryList
5776
5777                print SeqResult['varyList']
5778
5779                print '**** ERROR - variable list for this histogram does not match previous'
5780
5781                raise Exception
5782
5783        parmDict.update(phaseDict)
5784
5785        parmDict.update(hapDict)
5786
5787        parmDict.update(histDict)
5788
5789        GetFprime(calcControls,Histo)
5790
5791        constrDict,constrFlag,fixedList = G2mv.InputParse([])        #constraints go here?
5792
5793        groups,parmlist = G2mv.GroupConstraints(constrDict)
5794
5795        G2mv.GenerateConstraints(groups,parmlist,varyList,constrDict,constrFlag,fixedList)
5796
5797        G2mv.Map2Dict(parmDict,varyList)
5798
5799   
5800
5801        while True:
5802
5803            begin = time.time()
5804
5805            values =  np.array(Dict2Values(parmDict, varyList))
5806
5807            Ftol = Controls['min dM/M']
5808
5809            Factor = Controls['shift factor']
5810
5811            maxCyc = Controls['max cyc']
5812
5813
5814
5815            if 'Jacobian' in Controls['deriv type']:           
5816
5817                result = so.leastsq(errRefine,values,Dfun=dervRefine,full_output=True,
5818
5819                    ftol=Ftol,col_deriv=True,factor=Factor,
5820
5821                    args=([Histo,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg))
5822
5823                ncyc = int(result[2]['nfev']/2)
5824
5825            elif 'Hessian' in Controls['deriv type']:
5826
5827                result = G2mth.HessianLSQ(errRefine,values,Hess=HessRefine,ftol=Ftol,maxcyc=maxCyc,
5828
5829                    args=([Histo,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg))
5830
5831                ncyc = result[2]['num cyc']+1                           
5832
5833            else:           #'numeric'
5834
5835                result = so.leastsq(errRefine,values,full_output=True,ftol=Ftol,epsfcn=1.e-8,factor=Factor,
5836
5837                    args=([Histo,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg))
5838
5839                ncyc = int(result[2]['nfev']/len(varyList))
5840
5841
5842
5843
5844
5845
5846
5847            runtime = time.time()-begin
5848
5849            chisq = np.sum(result[2]['fvec']**2)
5850
5851            Values2Dict(parmDict, varyList, result[0])
5852
5853            G2mv.Dict2Map(parmDict,varyList)
5854
5855           
5856
5857            Rwp = np.sqrt(chisq/Histo['sumwYo'])*100.      #to %
5858
5859            GOF = chisq/(Histo['Nobs']-len(varyList))
5860
5861            print '\n Refinement results for histogram: v'+histogram
5862
5863            print 135*'-'
5864
5865            print ' Number of function calls:',result[2]['nfev'],' Number of observations: ',Histo['Nobs'],' Number of parameters: ',len(varyList)
5866
5867            print ' Refinement time = %8.3fs, %8.3fs/cycle, for %d cycles'%(runtime,runtime/ncyc,ncyc)
5868
5869            print ' wRp = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f'%(Rwp,chisq,GOF)
5870
5871            try:
5872
5873                covMatrix = result[1]*GOF
5874
5875                sig = np.sqrt(np.diag(covMatrix))
5876
5877                if np.any(np.isnan(sig)):
5878
5879                    print '*** Least squares aborted - some invalid esds possible ***'
5880
5881                    ifPrint = True
5882
5883                break                   #refinement succeeded - finish up!
5884
5885            except TypeError:          #result[1] is None on singular matrix
5886
5887                print '**** Refinement failed - singular matrix ****'
5888
5889                if 'Hessian' in Controls['deriv type']:
5890
5891                    for i in result[2]['psing'].reverse():
5892
5893                            print 'Removing parameter: ',varyList[i]
5894
5895                            del(varyList[i])                   
5896
5897                else:
5898
5899                    Ipvt = result[2]['ipvt']
5900
5901                    for i,ipvt in enumerate(Ipvt):
5902
5903                        if not np.sum(result[2]['fjac'],axis=1)[i]:
5904
5905                            print 'Removing parameter: ',varyList[ipvt-1]
5906
5907                            del(varyList[ipvt-1])
5908
5909                            break
5910
5911   
5912
5913        GetFobsSq(Histo,Phases,parmDict,calcControls)
5914
5915        sigDict = dict(zip(varyList,sig))
5916
5917        newCellDict = GetNewCellParms(parmDict,varyList)
5918
5919        newAtomDict = ApplyXYZshifts(parmDict,varyList)
5920
5921        covData = {'variables':result[0],'varyList':varyList,'sig':sig,
5922
5923            'covMatrix':covMatrix,'title':histogram,'newAtomDict':newAtomDict,'newCellDict':newCellDict}
5924
5925        SetHistogramPhaseData(parmDict,sigDict,Phases,Histo,ifPrint)
5926
5927        SetHistogramData(parmDict,sigDict,Histo,ifPrint)
5928
5929        SeqResult[histogram] = covData
5930
5931        SetUsedHistogramsAndPhases(GPXfile,Histo,Phases,covData,makeBack)
5932
5933        makeBack = False
5934
5935    SetSeqResult(GPXfile,Histograms,SeqResult)
5936
5937
5938
5939def DistAngle(DisAglCtls,DisAglData):
5940
5941    import numpy.ma as ma
5942
5943   
5944
5945    def ShowBanner(name):
5946
5947        print 80*'*'
5948
5949        print '   Interatomic Distances and Angles for phase '+name
5950
5951        print 80*'*','\n'
5952
5953
5954
5955    ShowBanner(DisAglCtls['Name'])
5956
5957    SGData = DisAglData['SGData']
5958
5959    SGtext = G2spc.SGPrint(SGData)
5960
5961    for line in SGtext: print line
5962
5963    Cell = DisAglData['Cell']
5964
5965   
5966
5967    Amat,Bmat = G2lat.cell2AB(Cell[:6])
5968
5969    covData = {}
5970
5971    if 'covData' in DisAglData:   
5972
5973        covData = DisAglData['covData']
5974
5975        covMatrix = covData['covMatrix']
5976
5977        varyList = covData['varyList']
5978
5979        pfx = str(DisAglData['pId'])+'::'
5980
5981        A = G2lat.cell2A(Cell[:6])
5982
5983        cellSig = getCellEsd(pfx,SGData,A,covData)
5984
5985        names = [' a = ',' b = ',' c = ',' alpha = ',' beta = ',' gamma = ',' Volume = ']
5986
5987        valEsd = [G2mth.ValEsd(Cell[i],cellSig[i],True) for i in range(7)]
5988
5989        line = '\n Unit cell:'
5990
5991        for name,vals in zip(names,valEsd):
5992
5993            line += name+vals 
5994
5995        print line
5996
5997    else: 
5998
5999        print '\n Unit cell: a = ','%.5f'%(Cell[0]),' b = ','%.5f'%(Cell[1]),' c = ','%.5f'%(Cell[2]), \
6000
6001            ' alpha = ','%.3f'%(Cell[3]),' beta = ','%.3f'%(Cell[4]),' gamma = ', \
6002
6003            '%.3f'%(Cell[5]),' volume = ','%.3f'%(Cell[6])
6004
6005    Factor = DisAglCtls['Factors']
6006
6007    Radii = dict(zip(DisAglCtls['AtomTypes'],zip(DisAglCtls['BondRadii'],DisAglCtls['AngleRadii'])))
6008
6009    Units = np.array([                   #is there a nicer way to make this?
6010
6011        [-1,-1,-1],[-1,-1,0],[-1,-1,1],[-1,0,-1],[-1,0,0],[-1,0,1],[-1,1,-1],[-1,1,0],[-1,1,1],
6012
6013        [0,-1,-1],[0,-1,0],[0,-1,1],[0,0,-1],[0,0,0],[0,0,1],[0,1,-1],[0,1,0],[0,1,1],
6014
6015        [1,-1,-1],[1,-1,0],[1,-1,1],[1,0,-1],[1,0,0],[1,0,1],[1,1,-1],[1,1,0],[1,1,1]])
6016
6017    origAtoms = DisAglData['OrigAtoms']
6018
6019    targAtoms = DisAglData['TargAtoms']
6020
6021    for Oatom in origAtoms:
6022
6023        OxyzNames = ''
6024
6025        IndBlist = []
6026
6027        Dist = []
6028
6029        Vect = []
6030
6031        VectA = []
6032
6033        angles = []
6034
6035        for Tatom in targAtoms:
6036
6037            Xvcov = []
6038
6039            TxyzNames = ''
6040
6041            if 'covData' in DisAglData:
6042
6043                OxyzNames = [pfx+'dAx:%d'%(Oatom[0]),pfx+'dAy:%d'%(Oatom[0]),pfx+'dAz:%d'%(Oatom[0])]
6044
6045                TxyzNames = [pfx+'dAx:%d'%(Tatom[0]),pfx+'dAy:%d'%(Tatom[0]),pfx+'dAz:%d'%(Tatom[0])]
6046
6047                Xvcov = G2mth.getVCov(OxyzNames+TxyzNames,varyList,covMatrix)
6048
6049            result = G2spc.GenAtom(Tatom[3:6],SGData,False,Move=False)
6050
6051            BsumR = (Radii[Oatom[2]][0]+Radii[Tatom[2]][0])*Factor[0]
6052
6053            AsumR = (Radii[Oatom[2]][1]+Radii[Tatom[2]][1])*Factor[1]
6054
6055            for Txyz,Top,Tunit in result:
6056
6057                Dx = (Txyz-np.array(Oatom[3:6]))+Units
6058
6059                dx = np.inner(Amat,Dx)
6060
6061                dist = ma.masked_less(np.sqrt(np.sum(dx**2,axis=0)),0.5)
6062
6063                IndB = ma.nonzero(ma.masked_greater(dist-BsumR,0.))
6064
6065                if np.any(IndB):
6066
6067                    for indb in IndB:
6068
6069                        for i in range(len(indb)):
6070
6071                            if str(dx.T[indb][i]) not in IndBlist:
6072
6073                                IndBlist.append(str(dx.T[indb][i]))
6074
6075                                unit = Units[indb][i]
6076
6077                                tunit = '[%2d%2d%2d]'%(unit[0]+Tunit[0],unit[1]+Tunit[1],unit[2]+Tunit[2])
6078
6079                                pdpx = G2mth.getDistDerv(Oatom[3:6],Tatom[3:6],Amat,unit,Top,SGData)
6080
6081                                sig = 0.0
6082
6083                                if len(Xvcov):
6084
6085                                    sig = np.sqrt(np.inner(pdpx,np.inner(Xvcov,pdpx)))
6086
6087                                Dist.append([Oatom[1],Tatom[1],tunit,Top,ma.getdata(dist[indb])[i],sig])
6088
6089                                if (Dist[-1][-1]-AsumR) <= 0.:
6090
6091                                    Vect.append(dx.T[indb][i]/Dist[-1][-2])
6092
6093                                    VectA.append([OxyzNames,np.array(Oatom[3:6]),TxyzNames,np.array(Tatom[3:6]),unit,Top])
6094
6095                                else:
6096
6097                                    Vect.append([0.,0.,0.])
6098
6099                                    VectA.append([])
6100
6101        Vect = np.array(Vect)
6102
6103        angles = np.zeros((len(Vect),len(Vect)))
6104
6105        angsig = np.zeros((len(Vect),len(Vect)))
6106
6107        for i,veca in enumerate(Vect):
6108
6109            if np.any(veca):
6110
6111                for j,vecb in enumerate(Vect):
6112
6113                    if np.any(vecb):
6114
6115                        angles[i][j],angsig[i][j] = G2mth.getAngSig(VectA[i],VectA[j],Amat,SGData,covData)
6116
6117        line = ''
6118
6119        for i,x in enumerate(Oatom[3:6]):
6120
6121            if len(Xvcov):
6122
6123                line += '%12s'%(G2mth.ValEsd(x,np.sqrt(Xvcov[i][i]),True))
6124
6125            else:
6126
6127                line += '%12.5f'%(x)
6128
6129        print '\n Distances & angles for ',Oatom[1],' at ',line
6130
6131        print 80*'*'
6132
6133        line = ''
6134
6135        for dist in Dist[:-1]:
6136
6137            line += '%12s'%(dist[1].center(12))
6138
6139        print '  To       cell +(sym. op.)      dist.  ',line
6140
6141        for i,dist in enumerate(Dist):
6142
6143            line = ''
6144
6145            for j,angle in enumerate(angles[i][0:i]):
6146
6147                sig = angsig[i][j]
6148
6149                if angle:
6150
6151                    if sig:
6152
6153                        line += '%12s'%(G2mth.ValEsd(angle,sig,True).center(12))
6154
6155                    else:
6156
6157                        val = '%.3f'%(angle)
6158
6159                        line += '%12s'%(val.center(12))
6160
6161                else:
6162
6163                    line += 12*' '
6164
6165            if dist[5]:            #sig exists!
6166
6167                val = G2mth.ValEsd(dist[4],dist[5])
6168
6169            else:
6170
6171                val = '%8.4f'%(dist[4])
6172
6173            print %8s%10s+(%4d) %12s'%(dist[1].ljust(8),dist[2].ljust(10),dist[3],val.center(12)),line
6174
6175
6176
6177def Torsion(TorsionData):
6178
6179
6180
6181    def ShowBanner(name):
6182
6183        print 80*'*'
6184
6185        print '   Torsion angle for phase '+name
6186
6187        print 80*'*'
6188
6189
6190
6191    ShowBanner(TorsionData['Name'])
6192
6193    SGData = TorsionData['SGData']
6194
6195    Cell = TorsionData['Cell']
6196
6197   
6198
6199    Amat,Bmat = G2lat.cell2AB(Cell[:6])
6200
6201    covData = {}
6202
6203    pfx = ''
6204
6205    if 'covData' in TorsionData:   
6206
6207        covData = TorsionData['covData']
6208
6209        covMatrix = covData['covMatrix']
6210
6211        varyList = covData['varyList']
6212
6213        pfx = str(TorsionData['pId'])+'::'
6214
6215    #find one end of 4 atom string - involved in longest distance
6216
6217    dist = {}
6218
6219    for i,X1 in enumerate(TorsionData['Datoms']):
6220
6221        for j,X2 in enumerate(TorsionData['Datoms'][:i]):
6222
6223            dist[np.sqrt(np.sum(np.inner(Amat,np.array(X2[3:6])-np.array(X1[3:6]))**2))] = [i,j]
6224
6225    sortdist = dist.keys()
6226
6227    sortdist.sort()
6228
6229    end = dist[sortdist[-1]][0]
6230
6231    #order atoms in distance from end - defines sequence of atoms for the torsion angle
6232
6233    dist = {}
6234
6235    X1 = TorsionData['Datoms'][end]
6236
6237    for i,X2 in enumerate(TorsionData['Datoms']):
6238
6239        dist[np.sqrt(np.sum(np.inner(Amat,np.array(X2[3:6])-np.array(X1[3:6]))**2))] = i
6240
6241    sortdist = dist.keys()
6242
6243    sortdist.sort()
6244
6245    Datoms = []
6246
6247    Oatoms = []
6248
6249    for d in sortdist:
6250
6251        atom = TorsionData['Datoms'][dist[d]]
6252
6253        symop = atom[-1].split('+')
6254
6255        if len(symop) == 1:
6256
6257            symop.append('0,0,0')       
6258
6259        symop[0] = int(symop[0])
6260
6261        symop[1] = eval(symop[1])
6262
6263        atom[-1] = symop
6264
6265        Datoms.append(atom)
6266
6267        oatom = TorsionData['Oatoms'][dist[d]]
6268
6269        names = ['','','']
6270
6271        if pfx:
6272
6273            names = [pfx+'dAx:'+str(oatom[0]),pfx+'dAy:'+str(oatom[0]),pfx+'dAz:'+str(oatom[0])]
6274
6275        oatom += [names,]
6276
6277        Oatoms.append(oatom)
6278
6279    Tors,sig = G2mth.GetTorsionSig(Oatoms,Datoms,Amat,SGData,covData)
6280
6281    print ' Torsion angle for atom sequence: ',[Datoms[i][1] for i in range(4)],'=',G2mth.ValEsd(Tors,sig)
6282
6283    print ' NB: Atom sequence determined by interatomic distances'
6284
6285               
6286
6287def BestPlane(PlaneData):
6288
6289
6290
6291    def ShowBanner(name):
6292
6293        print 80*'*'
6294
6295        print '   Best plane result for phase '+name
6296
6297        print 80*'*','\n'
6298
6299
6300
6301    ShowBanner(PlaneData['Name'])
6302
6303
6304
6305    Cell = PlaneData['Cell']   
6306
6307    Amat,Bmat = G2lat.cell2AB(Cell[:6])       
6308
6309    Atoms = PlaneData['Atoms']
6310
6311    sumXYZ = np.zeros(3)
6312
6313    XYZ = []
6314
6315    Natoms = len(Atoms)
6316
6317    for atom in Atoms:
6318
6319        xyz = np.array(atom[3:6])
6320
6321        XYZ.append(xyz)
6322
6323        sumXYZ += xyz
6324
6325    sumXYZ /= Natoms
6326
6327    XYZ = np.array(XYZ)-sumXYZ
6328
6329    XYZ = np.inner(Amat,XYZ).T
6330
6331    Zmat = np.zeros((3,3))
6332
6333    for i,xyz in enumerate(XYZ):
6334
6335        Zmat += np.outer(xyz.T,xyz)
6336
6337    print ' Selected atoms centered at %10.5f %10.5f %10.5f'%(sumXYZ[0],sumXYZ[1],sumXYZ[2])
6338
6339    Evec,Emat = nl.eig(Zmat)
6340
6341    Evec = np.sqrt(Evec)/(Natoms-3)
6342
6343    Order = np.argsort(Evec)
6344
6345    XYZ = np.inner(XYZ,Emat.T).T
6346
6347    XYZ = np.array([XYZ[Order[2]],XYZ[Order[1]],XYZ[Order[0]]]).T
6348
6349    print ' Atoms in Cartesian best plane coordinates:'
6350
6351    print ' Name         X         Y         Z'
6352
6353    for i,xyz in enumerate(XYZ):
6354
6355        print ' %6s%10.3f%10.3f%10.3f'%(Atoms[i][1].ljust(6),xyz[0],xyz[1],xyz[2])
6356
6357    print '\n Best plane RMS X =%8.3f, Y =%8.3f, Z =%8.3f'%(Evec[Order[2]],Evec[Order[1]],Evec[Order[0]])   
6358
6359
6360
6361           
6362
6363def main():
6364
6365    arg = sys.argv
6366
6367    try: 
6368
6369        mkl.set_num_threads(MaxThreads)
6370
6371        print "Max threads ",mkl.get_max_threads()
6372
6373        ActualThreads = mkl.get_max_threads()
6374
6375    except:
6376
6377        print "MKL module not present"
6378
6379        ActualThreads = 1
6380
6381    if len(arg) > 1:
6382
6383        GPXfile = arg[1]
6384
6385        if not ospath.exists(GPXfile):
6386
6387            print 'ERROR - ',GPXfile," doesn't exist!"
6388
6389            exit()
6390
6391        GPXpath = ospath.dirname(arg[1])
6392
6393        Refine(GPXfile,None)
6394
6395    else:
6396
6397        print 'ERROR - missing filename'
6398
6399        exit()
6400
6401    print "Requested threads ",MaxThreads
6402
6403    print "Actual threads ",ActualThreads
6404
6405    print "Max subprocesses ",MaxProcess
6406
6407    print "Done"
6408
6409         
6410
6411if __name__ == '__main__':
6412
6413    main()
6414
Note: See TracBrowser for help on using the repository browser.