Changeset 494


Ignore:
Timestamp:
Feb 24, 2012 2:25:49 PM (10 years ago)
Author:
vondreele
Message:

brians version

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIIstruct.py

    r493 r494  
    11#GSASIIstructure - structure computation routines
     2
    23########### SVN repository information ###################
     4
    35# $Date$
     6
    47# $Author$
     8
    59# $Revision$
     10
    611# $URL$
     12
    713# $Id$
     14
    815########### SVN repository information ###################
     16
     17MaxThreads = 1
     18
     19MaxProcess = 2
     20
     21
     22
    923import sys
     24
    1025import os
     26
    1127import os.path as ospath
     28
    1229import time
     30
    1331import math
     32
    1433import cPickle
     34
     35import multiprocessing as mp
     36
    1537import numpy as np
     38
    1639import numpy.linalg as nl
     40
    1741import scipy.optimize as so
     42
    1843import GSASIIpath
     44
    1945import GSASIIElem as G2el
     46
    2047import GSASIIlattice as G2lat
     48
    2149import GSASIIspc as G2spc
     50
    2251import GSASIIpwd as G2pwd
     52
    2353import GSASIImapvars as G2mv
     54
    2455import GSASIImath as G2mth
    2556
     57#import buggerwx
     58
     59try:
     60
     61    import mkl
     62
     63except:
     64
     65    print "MKL module not present"
     66
     67
     68
     69
     70
    2671sind = lambda x: np.sin(x*np.pi/180.)
     72
    2773cosd = lambda x: np.cos(x*np.pi/180.)
     74
    2875tand = lambda x: np.tan(x*np.pi/180.)
     76
    2977asind = lambda x: 180.*np.arcsin(x)/np.pi
     78
    3079acosd = lambda x: 180.*np.arccos(x)/np.pi
     80
    3181atan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
    3282
    3383
     84
     85
     86
    3487def GetControls(GPXfile):
     88
    3589    ''' Returns dictionary of control items found in GSASII gpx file
     90
    3691    input:
     92
    3793        GPXfile = .gpx full file name
     94
    3895    return:
     96
    3997        Controls = dictionary of control items
     98
    4099    '''
     100
    41101    Controls = {'deriv type':'analytical','min dM/M':0.0001,'shift factor':1.}
     102
    42103    file = open(GPXfile,'rb')
     104
    43105    while True:
     106
    44107        try:
     108
    45109            data = cPickle.load(file)
     110
    46111        except EOFError:
     112
    47113            break
     114
    48115        datum = data[0]
     116
    49117        if datum[0] == 'Controls':
     118
    50119            Controls.update(datum[1])
     120
    51121    file.close()
     122
    52123    return Controls
     124
    53125   
     126
    54127def GetConstraints(GPXfile):
     128
    55129    constList = []
     130
    56131    file = open(GPXfile,'rb')
     132
    57133    while True:
     134
    58135        try:
     136
    59137            data = cPickle.load(file)
     138
    60139        except EOFError:
     140
    61141            break
     142
    62143        datum = data[0]
     144
    63145        if datum[0] == 'Constraints':
     146
    64147            constDict = datum[1]
     148
    65149            for item in constDict:
     150
    66151                constList += constDict[item]
     152
    67153    file.close()
     154
    68155    constDict = []
     156
    69157    constFlag = []
     158
    70159    fixedList = []
     160
    71161    for item in constList:
     162
    72163        if item[-2]:
     164
    73165            fixedList.append(str(item[-2]))
     166
    74167        else:
     168
    75169            fixedList.append('0')
     170
    76171        if item[-1]:
     172
    77173            constFlag.append(['VARY'])
     174
    78175        else:
     176
    79177            constFlag.append([])
     178
    80179        itemDict = {}
     180
    81181        for term in item[:-2]:
     182
    82183            itemDict[term[1]] = term[0]
     184
    83185        constDict.append(itemDict)
     186
    84187    return constDict,constFlag,fixedList
     188
    85189   
     190
    86191def GetPhaseNames(GPXfile):
     192
    87193    ''' Returns a list of phase names found under 'Phases' in GSASII gpx file
     194
    88195    input:
     196
    89197        GPXfile = gpx full file name
     198
    90199    return:
     200
    91201        PhaseNames = list of phase names
     202
    92203    '''
     204
    93205    file = open(GPXfile,'rb')
     206
    94207    PhaseNames = []
     208
    95209    while True:
     210
    96211        try:
     212
    97213            data = cPickle.load(file)
     214
    98215        except EOFError:
     216
    99217            break
     218
    100219        datum = data[0]
     220
    101221        if 'Phases' == datum[0]:
     222
    102223            for datus in data[1:]:
     224
    103225                PhaseNames.append(datus[0])
     226
    104227    file.close()
     228
    105229    return PhaseNames
    106230
     231
     232
    107233def GetAllPhaseData(GPXfile,PhaseName):
     234
    108235    ''' Returns the entire dictionary for PhaseName from GSASII gpx file
     236
    109237    input:
     238
    110239        GPXfile = gpx full file name
     240
    111241        PhaseName = phase name
     242
    112243    return:
     244
    113245        phase dictionary
     246
    114247    '''       
     248
    115249    file = open(GPXfile,'rb')
     250
    116251    General = {}
     252
    117253    Atoms = []
     254
    118255    while True:
     256
    119257        try:
     258
    120259            data = cPickle.load(file)
     260
    121261        except EOFError:
     262
    122263            break
     264
    123265        datum = data[0]
     266
    124267        if 'Phases' == datum[0]:
     268
    125269            for datus in data[1:]:
     270
    126271                if datus[0] == PhaseName:
     272
    127273                    break
     274
    128275    file.close()
     276
    129277    return datus[1]
     278
    130279   
     280
    131281def GetHistograms(GPXfile,hNames):
     282
    132283    """ Returns a dictionary of histograms found in GSASII gpx file
     284
    133285    input:
     286
    134287        GPXfile = .gpx full file name
     288
    135289        hNames = list of histogram names
     290
    136291    return:
     292
    137293        Histograms = dictionary of histograms (types = PWDR & HKLF)
     294
    138295    """
     296
    139297    file = open(GPXfile,'rb')
     298
    140299    Histograms = {}
     300
    141301    while True:
     302
    142303        try:
     304
    143305            data = cPickle.load(file)
     306
    144307        except EOFError:
     308
    145309            break
     310
    146311        datum = data[0]
     312
    147313        hist = datum[0]
     314
    148315        if hist in hNames:
     316
    149317            if 'PWDR' in hist[:4]:
     318
    150319                PWDRdata = {}
     320
    151321                PWDRdata['Data'] = datum[1][1]          #powder data arrays
     322
    152323                PWDRdata[data[2][0]] = data[2][1]       #Limits
     324
    153325                PWDRdata[data[3][0]] = data[3][1]       #Background
     326
    154327                PWDRdata[data[4][0]] = data[4][1]       #Instrument parameters
     328
    155329                PWDRdata[data[5][0]] = data[5][1]       #Sample parameters
     330
    156331                try:
     332
    157333                    PWDRdata[data[9][0]] = data[9][1]       #Reflection lists might be missing
     334
    158335                except IndexError:
     336
    159337                    PWDRdata['Reflection lists'] = {}
     338
    160339   
     340
    161341                Histograms[hist] = PWDRdata
     342
    162343            elif 'HKLF' in hist[:4]:
     344
    163345                HKLFdata = []
     346
    164347                datum = data[0]
     348
    165349                HKLFdata = datum[1:][0]
     350
    166351                Histograms[hist] = HKLFdata           
     352
    167353    file.close()
     354
    168355    return Histograms
     356
    169357   
     358
    170359def GetHistogramNames(GPXfile,hType):
     360
    171361    """ Returns a list of histogram names found in GSASII gpx file
     362
    172363    input:
     364
    173365        GPXfile = .gpx full file name
     366
    174367        hType = list ['PWDR','HKLF']
     368
    175369    return:
     370
    176371        HistogramNames = list of histogram names (types = PWDR & HKLF)
     372
    177373    """
     374
    178375    file = open(GPXfile,'rb')
     376
    179377    HistogramNames = []
     378
    180379    while True:
     380
    181381        try:
     382
    182383            data = cPickle.load(file)
     384
    183385        except EOFError:
     386
    184387            break
     388
    185389        datum = data[0]
     390
    186391        if datum[0][:4] in hType:
     392
    187393            HistogramNames.append(datum[0])
     394
    188395    file.close()
     396
    189397    return HistogramNames
     398
    190399   
     400
    191401def GetUsedHistogramsAndPhases(GPXfile):
     402
    192403    ''' Returns all histograms that are found in any phase
     404
    193405    and any phase that uses a histogram
     406
    194407    input:
     408
    195409        GPXfile = .gpx full file name
     410
    196411    return:
     412
    197413        Histograms = dictionary of histograms as {name:data,...}
     414
    198415        Phases = dictionary of phases that use histograms
     416
    199417    '''
     418
    200419    phaseNames = GetPhaseNames(GPXfile)
     420
    201421    histoList = GetHistogramNames(GPXfile,['PWDR','HKLF'])
     422
    202423    allHistograms = GetHistograms(GPXfile,histoList)
     424
    203425    phaseData = {}
     426
    204427    for name in phaseNames:
     428
    205429        phaseData[name] =  GetAllPhaseData(GPXfile,name)
     430
    206431    Histograms = {}
     432
    207433    Phases = {}
     434
    208435    for phase in phaseData:
     436
    209437        Phase = phaseData[phase]
     438
    210439        if Phase['Histograms']:
     440
    211441            if phase not in Phases:
     442
    212443                pId = phaseNames.index(phase)
     444
    213445                Phase['pId'] = pId
     446
    214447                Phases[phase] = Phase
     448
    215449            for hist in Phase['Histograms']:
     450
    216451                if hist not in Histograms:
     452
    217453                    Histograms[hist] = allHistograms[hist]
     454
    218455                    #future restraint, etc. histograms here           
     456
    219457                    hId = histoList.index(hist)
     458
    220459                    Histograms[hist]['hId'] = hId
     460
    221461    return Histograms,Phases
     462
    222463   
    223 def getBackupName2(GPXfile,makeBack=True):      #not work correctly
     464
     465def GPXBackup(GPXfile,makeBack=True):
     466
     467    import distutils.file_util as dfu
     468
    224469    GPXpath,GPXname = ospath.split(GPXfile)
     470
    225471    if GPXpath == '': GPXpath = '.'
     472
    226473    Name = ospath.splitext(GPXname)[0]
     474
    227475    files = os.listdir(GPXpath)
     476
    228477    last = 0
     478
    229479    for name in files:
     480
    230481        name = name.split('.')
     482
    231483        if len(name) >= 3 and name[0] == Name and 'bak' in name[-2]:
     484
    232485            if makeBack:
     486
    233487                last = max(last,int(name[-2].strip('bak'))+1)
     488
    234489            else:
     490
    235491                last = max(last,int(name[-2].strip('bak')))
     492
    236493    GPXback = ospath.join(GPXpath,GPXname.rstrip('.'.join(name[-2:]))+'.bak'+str(last)+'.gpx')
     494
     495    dfu.copy_file(GPXfile,GPXback)
     496
    237497    return GPXback
    238498
    239 def getBackupName(GPXfile,makeBack):       #recovered old one
    240     GPXpath,GPXname = ospath.split(GPXfile)
    241     if GPXpath == '': GPXpath = '.'
    242     Name = ospath.splitext(GPXname)[0]
    243     files = os.listdir(GPXpath)
    244     last = 0
    245     for name in files:
    246         name = name.split('.')
    247         if len(name) == 3 and name[0] == Name and 'bak' in name[1]:
    248             if makeBack:
    249                 last = max(last,int(name[1].strip('bak'))+1)
     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
    2501643            else:
    251                 last = max(last,int(name[1].strip('bak')))
    252     GPXback = ospath.join(GPXpath,ospath.splitext(GPXname)[0]+'.bak'+str(last)+'.gpx')
    253     return GPXback   
     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
    2541683       
    255 def GPXBackup(GPXfile,makeBack=True):
    256     import distutils.file_util as dfu
    257     GPXback = getBackupName(GPXfile,makeBack)
    258     dfu.copy_file(GPXfile,GPXback)
    259     return GPXback
    260 
    261 def SetUsedHistogramsAndPhases(GPXfile,Histograms,Phases,CovData,makeBack=True):
    262     ''' Updates gpxfile from all histograms that are found in any phase
    263     and any phase that used a histogram
     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
    2643321    input:
    265         GPXfile = .gpx full file name
    266         Histograms = dictionary of histograms as {name:data,...}
    267         Phases = dictionary of phases that use histograms
    268         CovData = dictionary of refined variables, varyList, & covariance matrix
    269         makeBack = True if new backup of .gpx file is to be made; else use the last one made
     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
    2703663    '''
    271                        
    272     GPXback = GPXBackup(GPXfile,makeBack)
    273     print '\n',135*'-'
    274     print 'Read from file:',GPXback
    275     print 'Save to file  :',GPXfile
    276     infile = open(GPXback,'rb')
    277     outfile = open(GPXfile,'wb')
     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
    2785463    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
    2795529        try:
    280             data = cPickle.load(infile)
    281         except EOFError:
    282             break
    283         datum = data[0]
    284 #        print 'read: ',datum[0]
    285         if datum[0] == 'Phases':
    286             for iphase in range(len(data)):
    287                 if data[iphase][0] in Phases:
    288                     phaseName = data[iphase][0]
    289                     data[iphase][1].update(Phases[phaseName])
    290         elif datum[0] == 'Covariance':
    291             data[0][1] = CovData
    292         try:
    293             histogram = Histograms[datum[0]]
    294 #            print 'found ',datum[0]
    295             data[0][1][1] = histogram['Data']
    296             for datus in data[1:]:
    297 #                print '    read: ',datus[0]
    298                 if datus[0] in ['Background','Instrument Parameters','Sample Parameters','Reflection Lists']:
    299                     datus[1] = histogram[datus[0]]
    300         except KeyError:
    301             pass
    302                                
    303         cPickle.dump(data,outfile,1)
    304     infile.close()
    305     outfile.close()
    306     print 'GPX file save successful'
     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
    3075613   
    308 def SetSeqResult(GPXfile,Histograms,SeqResult):
    309     GPXback = GPXBackup(GPXfile)
    310     print '\n',135*'-'
    311     print 'Read from file:',GPXback
    312     print 'Save to file  :',GPXfile
    313     infile = open(GPXback,'rb')
    314     outfile = open(GPXfile,'wb')
    315     while True:
    316         try:
    317             data = cPickle.load(infile)
    318         except EOFError:
    319             break
    320         datum = data[0]
    321         if datum[0] == 'Sequental results':
    322             data[0][1] = SeqResult
    323         try:
    324             histogram = Histograms[datum[0]]
    325             data[0][1][1] = histogram['Data']
    326             for datus in data[1:]:
    327                 if datus[0] in ['Background','Instrument Parameters','Sample Parameters','Reflection Lists']:
    328                     datus[1] = histogram[datus[0]]
    329         except KeyError:
    330             pass
    331                                
    332         cPickle.dump(data,outfile,1)
    333     infile.close()
    334     outfile.close()
    335     print 'GPX file save successful'
    336                        
    337 def GetPWDRdata(GPXfile,PWDRname):
    338     ''' Returns powder data from GSASII gpx file
    339     input:
    340         GPXfile = .gpx full file name
    341         PWDRname = powder histogram name as obtained from GetHistogramNames
    342     return:
    343         PWDRdata = powder data dictionary with:
    344             Data - powder data arrays, Limits, Instrument Parameters, Sample Parameters
    345        
    346     '''
    347     file = open(GPXfile,'rb')
    348     PWDRdata = {}
    349     while True:
    350         try:
    351             data = cPickle.load(file)
    352         except EOFError:
    353             break
    354         datum = data[0]
    355         if datum[0] == PWDRname:
    356             PWDRdata['Data'] = datum[1][1]          #powder data arrays
    357             PWDRdata[data[2][0]] = data[2][1]       #Limits
    358             PWDRdata[data[3][0]] = data[3][1]       #Background
    359             PWDRdata[data[4][0]] = data[4][1]       #Instrument parameters
    360             PWDRdata[data[5][0]] = data[5][1]       #Sample parameters
     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
    3615871            try:
    362                 PWDRdata[data[9][0]] = data[9][1]       #Reflection lists might be missing
    363             except IndexError:
    364                 PWDRdata['Reflection lists'] = {}
    365     file.close()
    366     return PWDRdata
     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
    3675911   
    368 def GetHKLFdata(GPXfile,HKLFname):
    369     ''' Returns single crystal data from GSASII gpx file
    370     input:
    371         GPXfile = .gpx full file name
    372         HKLFname = single crystal histogram name as obtained from GetHistogramNames
    373     return:
    374         HKLFdata = single crystal data list of reflections: for each reflection:
    375             HKLF = [np.array([h,k,l]),FoSq,sigFoSq,FcSq,Fcp,Fcpp,phase]
    376     '''
    377     file = open(GPXfile,'rb')
    378     HKLFdata = []
    379     while True:
    380         try:
    381             data = cPickle.load(file)
    382         except EOFError:
    383             break
    384         datum = data[0]
    385         if datum[0] == HKLFname:
    386             HKLFdata = datum[1:][0]
    387     file.close()
    388     return HKLFdata
     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
    3895943   
    390 def ShowBanner():
    391     print 80*'*'
    392     print '   General Structure Analysis System-II Crystal Structure Refinement'
    393     print '              by Robert B. Von Dreele & Brian H. Toby'
    394     print '                Argonne National Laboratory(C), 2010'
    395     print ' This product includes software developed by the UChicago Argonne, LLC,'
    396     print '            as Operator of Argonne National Laboratory.'
    397     print 80*'*','\n'
    398 
    399 def ShowControls(Controls):
    400     print ' Least squares controls:'
    401     print ' Refinement type: ',Controls['deriv type']
    402     if 'Hessian' in Controls['deriv type']:
    403         print ' Maximum number of cycles:',Controls['max cyc']
     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
    4046395    else:
    405         print ' Minimum delta-M/M for convergence: ','%.2g'%(Controls['min dM/M'])
    406     print ' Initial shift factor: ','%.3f'%(Controls['shift factor'])
    407    
    408 def GetFFtable(General):
    409     ''' returns a dictionary of form factor data for atom types found in General
    410     input:
    411         General = dictionary of phase info.; includes AtomTypes
    412     return:
    413         FFtable = dictionary of form factor data; key is atom type
    414     '''
    415     atomTypes = General['AtomTypes']
    416     FFtable = {}
    417     for El in atomTypes:
    418         FFs = G2el.GetFormFactorCoeff(El.split('+')[0].split('-')[0])
    419         for item in FFs:
    420             if item['Symbol'] == El.upper():
    421                 FFtable[El] = item
    422     return FFtable
    423    
    424 def GetBLtable(General):
    425     ''' returns a dictionary of neutron scattering length data for atom types & isotopes found in General
    426     input:
    427         General = dictionary of phase info.; includes AtomTypes & Isotopes
    428     return:
    429         BLtable = dictionary of scattering length data; key is atom type
    430     '''
    431     atomTypes = General['AtomTypes']
    432     BLtable = {}
    433     isotopes = General['Isotopes']
    434     isotope = General['Isotope']
    435     for El in atomTypes:
    436         BLtable[El] = [isotope[El],isotopes[El][isotope[El]]]
    437     return BLtable
    438        
    439 def GetPawleyConstr(SGLaue,PawleyRef,pawleyVary):
    440     if SGLaue in ['-1','2/m','mmm']:
    441         return                      #no Pawley symmetry required constraints
    442     for i,varyI in enumerate(pawleyVary):
    443         refI = int(varyI.split(':')[-1])
    444         ih,ik,il = PawleyRef[refI][:3]
    445         for varyJ in pawleyVary[0:i]:
    446             refJ = int(varyJ.split(':')[-1])
    447             jh,jk,jl = PawleyRef[refJ][:3]
    448             if SGLaue in ['4/m','4/mmm']:
    449                 isum = ih**2+ik**2
    450                 jsum = jh**2+jk**2
    451                 if abs(il) == abs(jl) and isum == jsum:
    452                     G2mv.StoreEquivalence(varyJ,(varyI,))
    453             elif SGLaue in ['3R','3mR']:
    454                 isum = ih**2+ik**2+il**2
    455                 jsum = jh**2+jk**2*jl**2
    456                 isum2 = ih*ik+ih*il+ik*il
    457                 jsum2 = jh*jk+jh*jl+jk*jl
    458                 if isum == jsum and isum2 == jsum2:
    459                     G2mv.StoreEquivalence(varyJ,(varyI,))
    460             elif SGLaue in ['3','3m1','31m','6/m','6/mmm']:
    461                 isum = ih**2+ik**2+ih*ik
    462                 jsum = jh**2+jk**2+jh*jk
    463                 if abs(il) == abs(jl) and isum == jsum:
    464                     G2mv.StoreEquivalence(varyJ,(varyI,))
    465             elif SGLaue in ['m3','m3m']:
    466                 isum = ih**2+ik**2+il**2
    467                 jsum = jh**2+jk**2+jl**2
    468                 if isum == jsum:
    469                     G2mv.StoreEquivalence(varyJ,(varyI,))
    470                    
    471 def cellVary(pfx,SGData):
    472     if SGData['SGLaue'] in ['-1',]:
    473         return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3',pfx+'A4',pfx+'A5']
    474     elif SGData['SGLaue'] in ['2/m',]:
    475         if SGData['SGUniq'] == 'a':
    476             return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3']
    477         elif SGData['SGUniq'] == 'b':
    478             return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A4']
    479         else:
    480             return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A5']
    481     elif SGData['SGLaue'] in ['mmm',]:
    482         return [pfx+'A0',pfx+'A1',pfx+'A2']
    483     elif SGData['SGLaue'] in ['4/m','4/mmm']:
    484         return [pfx+'A0',pfx+'A2']
    485     elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
    486         return [pfx+'A0',pfx+'A2']
    487     elif SGData['SGLaue'] in ['3R', '3mR']:
    488         return [pfx+'A0',pfx+'A3']                       
    489     elif SGData['SGLaue'] in ['m3m','m3']:
    490         return [pfx+'A0',]
    491                    
    492 def GetPhaseData(PhaseData,Print=True):
    493            
    494     def PrintFFtable(FFtable):
    495         print '\n X-ray scattering factors:'
    496         print '   Symbol     fa                                      fb                                      fc'
    497         print 99*'-'
    498         for Ename in FFtable:
    499             ffdata = FFtable[Ename]
    500             fa = ffdata['fa']
    501             fb = ffdata['fb']
    502             print ' %8s %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f' %  \
    503                 (Ename.ljust(8),fa[0],fa[1],fa[2],fa[3],fb[0],fb[1],fb[2],fb[3],ffdata['fc'])
    504                
    505     def PrintBLtable(BLtable):
    506         print '\n Neutron scattering factors:'
    507         print '   Symbol   isotope       mass       b       resonant terms'
    508         print 99*'-'
    509         for Ename in BLtable:
    510             bldata = BLtable[Ename]
    511             isotope = bldata[0]
    512             mass = bldata[1][0]
    513             blen = bldata[1][1]
    514             bres = []
    515             if len(bldata[1]) > 2:
    516                 bres = bldata[1][2:]
    517             line = ' %8s%11s %10.3f %8.3f'%(Ename.ljust(8),isotope.center(11),mass,blen)
    518             for item in bres:
    519                 line += '%10.5g'%(item)
    520             print line
    521                
    522     def PrintAtoms(General,Atoms):
    523         print '\n Atoms:'
    524         line = '   name    type  refine?   x         y         z    '+ \
    525             '  frac site sym  mult I/A   Uiso     U11     U22     U33     U12     U13     U23'
    526         if General['Type'] == 'magnetic':
    527             line += '   Mx     My     Mz'
    528         elif General['Type'] == 'macromolecular':
    529             line = ' res no  residue  chain '+line
    530         print line
    531         if General['Type'] == 'nuclear':
    532             print 135*'-'
    533             for i,at in enumerate(Atoms):
    534                 line = '%7s'%(at[0])+'%7s'%(at[1])+'%7s'%(at[2])+'%10.5f'%(at[3])+'%10.5f'%(at[4])+ \
    535                     '%10.5f'%(at[5])+'%8.3f'%(at[6])+'%7s'%(at[7])+'%5d'%(at[8])+'%5s'%(at[9])
    536                 if at[9] == 'I':
    537                     line += '%8.4f'%(at[10])+48*' '
    538                 else:
    539                     line += 8*' '
    540                     for j in range(6):
    541                         line += '%8.4f'%(at[11+j])
    542                 print line
    543        
    544     def PrintTexture(textureData):
    545         topstr = '\n Spherical harmonics texture: Order:' + \
    546             str(textureData['Order'])
    547         if textureData['Order']:
    548             print topstr+' Refine? '+str(textureData['SH Coeff'][0])
    549         else:
    550             print topstr
    551             return
    552         names = ['omega','chi','phi']
    553         line = '\n'
    554         for name in names:
    555             line += ' SH '+name+':'+'%12.4f'%(textureData['Sample '+name][1])+' Refine? '+str(textureData['Sample '+name][0])
    556         print line
    557         print '\n Texture coefficients:'
    558         ptlbls = ' names :'
    559         ptstr =  ' values:'
    560         SHcoeff = textureData['SH Coeff'][1]
    561         for item in SHcoeff:
    562             ptlbls += '%12s'%(item)
    563             ptstr += '%12.4f'%(SHcoeff[item])
    564         print ptlbls
    565         print ptstr   
    566        
    567     if Print: print ' Phases:'
    568     phaseVary = []
    569     phaseDict = {}
    570     phaseConstr = {}
    571     pawleyLookup = {}
    572     FFtables = {}                   #scattering factors - xrays
    573     BLtables = {}                   # neutrons
    574     Natoms = {}
    575     AtMults = {}
    576     AtIA = {}
    577     shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
    578     SamSym = dict(zip(shModels,['0','-1','2/m','mmm']))
    579     for name in PhaseData:
    580         General = PhaseData[name]['General']
    581         pId = PhaseData[name]['pId']
    582         pfx = str(pId)+'::'
    583         FFtable = GetFFtable(General)
    584         BLtable = GetBLtable(General)
    585         FFtables.update(FFtable)
    586         BLtables.update(BLtable)
    587         Atoms = PhaseData[name]['Atoms']
    588         try:
    589             PawleyRef = PhaseData[name]['Pawley ref']
    590         except KeyError:
    591             PawleyRef = []
    592         SGData = General['SGData']
    593         SGtext = G2spc.SGPrint(SGData)
    594         cell = General['Cell']
    595         A = G2lat.cell2A(cell[1:7])
    596         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]})
    597         if cell[0]:
    598             phaseVary += cellVary(pfx,SGData)
    599         Natoms[pfx] = 0
    600         if Atoms:
    601             if General['Type'] == 'nuclear':
    602                 Natoms[pfx] = len(Atoms)
    603                 for i,at in enumerate(Atoms):
    604                     phaseDict.update({pfx+'Atype:'+str(i):at[1],pfx+'Afrac:'+str(i):at[6],pfx+'Amul:'+str(i):at[8],
    605                         pfx+'Ax:'+str(i):at[3],pfx+'Ay:'+str(i):at[4],pfx+'Az:'+str(i):at[5],
    606                         pfx+'dAx:'+str(i):0.,pfx+'dAy:'+str(i):0.,pfx