Changeset 380


Ignore:
Timestamp:
Oct 3, 2011 3:43:49 PM (10 years ago)
Author:
vondreele
Message:

completed Rietveld refinement version

Location:
trunk
Files:
14 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASII.py

    r373 r380  
    14531453        parmDict = {}
    14541454        Histograms,Phases = self.GetUsedHistogramsAndPhasesfromTree()
    1455         phaseVary,phaseDict,pawleyLookup = G2str.GetPhaseData(Phases,Print=False)       
     1455        Natoms,phaseVary,phaseDict,pawleyLookup,FFtable = G2str.GetPhaseData(Phases,Print=False)       
    14561456        hapVary,hapDict,controlDict = G2str.GetHistogramPhaseData(Phases,Histograms,Print=False)
    14571457        histVary,histDict,controlDict = G2str.GetHistogramData(Histograms,Print=False)
     
    14631463            if parm.split(':')[-1] in ['Azimuth','Gonio. radius','Lam1','Lam2']:
    14641464                parmDict[parm] = [parmDict[parm],' ']
     1465            elif parm.split(':')[-2] in ['Ax','Ay','Az']:
     1466                parmDict[parm] = [parmDict[parm],' ']
    14651467            elif parm in varyList:
    14661468                parmDict[parm] = [parmDict[parm],'True']
     
    14701472        try:
    14711473            if dlg.ShowModal() == wx.ID_OK:
    1472                 print 'do something with changes??'
    1473         finally:
    1474             dlg.Destroy()
    1475 
    1476 
     1474                print 'do something with changes?? - No!'
     1475        finally:
     1476            dlg.Destroy()
    14771477       
    14781478    def OnRefine(self,event):
     
    14801480        #works - but it'd be better if it could restore plots
    14811481        dlg = wx.ProgressDialog('Residual','Powder profile Rwp =',101.0,
    1482             style = wx.PD_ELAPSED_TIME|wx.PD_AUTO_HIDE|wx.PD_REMAINING_TIME|wx.PD_CAN_ABORT)
     1482            style = wx.PD_ELAPSED_TIME|wx.PD_AUTO_HIDE|wx.PD_CAN_ABORT)
    14831483        screenSize = wx.ClientDisplayRect()
    14841484        Size = dlg.GetSize()
  • trunk/GSASIIElem.py

    r327 r380  
    4646        S = FFdata.readline()
    4747        if S[3:5] == ElS:
    48             if S[5:6] != '_':
     48            if S[5:6] != '_' and S[8] not in ['N','M']:
    4949                Z=int(S[:2])
    5050                Symbol = S[3:7].strip()
  • trunk/GSASIIlattice.py

    r379 r380  
    208208def U6toUij(U6):
    209209    """Fill matrix (Uij) from U6 = [U11,U22,U33,U12,U13,U23]
    210     returns
     210    NB: there is a non numpy version in GSASIIspc: U2Uij
     211    input:
     212        U6 - 6 terms of u11,u22,...
     213    returns:
     214        Uij - numpy [3][3] array of uij
    211215    """
    212216    U = np.zeros(shape=(3,3))
     
    216220        [U6[4],  U6[5],  U6[2]]]
    217221    return U
    218        
     222
     223def UijtoU6(U):
     224    """Fill vector [U11,U22,U33,U12,U13,U23] from Uij
     225    NB: there is a non numpy version in GSASIIspc: Uij2U
     226    """
     227    U6 = np.zeros(6)
     228    U6 = [U[0][0],U[1][1],U[2][2],U[0][1],U[0][2],U[1][2]]
     229    return U6
     230
    219231def Uij2betaij(Uij,G):
    220232    """
  • trunk/GSASIIphsGUI.py

    r379 r380  
    485485            3*[wg.GRID_VALUE_FLOAT+':10,5',]+[wg.GRID_VALUE_FLOAT+':10,4', #x,y,z,frac
    486486            wg.GRID_VALUE_STRING,wg.GRID_VALUE_STRING,wg.GRID_VALUE_CHOICE+":I,A",]
    487         Types += 7*[wg.GRID_VALUE_FLOAT+':10,4',]
     487        Types += 7*[wg.GRID_VALUE_FLOAT+':10,5',]
    488488        colLabels = ['Name','Type','refine','x','y','z','frac','site sym','mult','I/A','Uiso','U11','U22','U33','U12','U13','U23']
    489489        if generalData['Type'] == 'magnetic':
     
    742742        for i,atom in enumerate(atomData):
    743743            table.append(atom)
    744             rowLabels.append(str(i+1))
     744            rowLabels.append(str(i))
    745745        atomTable = G2gd.Table(table,rowLabels=rowLabels,colLabels=colLabels,types=Types)
    746746        Atoms.SetTable(atomTable, True)
     
    11881188        for i,atom in enumerate(drawingData['Atoms']):
    11891189            table.append(atom[:colLabels.index('I/A')+1])
    1190             rowLabels.append(str(i+1))
     1190            rowLabels.append(str(i))
    11911191
    11921192        atomTable = G2gd.Table(table,rowLabels=rowLabels,colLabels=colLabels,types=Types)
  • trunk/GSASIIpwdGUI.py

    r373 r380  
    12931293    rowLabels = []
    12941294    refList = []
    1295     for h,k,l,m,d,pos,sig,gam,fo,fc,x in data[self.RefList]:
    1296         refList.append([h,k,l,m,d,pos,sig,gam,fo,fc])
     1295    for h,k,l,m,d,pos,sig,gam,fo,fc,phi,x,x in data[self.RefList]:
     1296        refList.append([h,k,l,m,d,pos,sig,gam,fo,fc,phi])
    12971297    for i in range(len(refList)): rowLabels.append(str(i))
    1298     colLabels = ['H','K','L','mul','d','pos','sig','gam','Fosq','Fcsq',]
    1299     Types = 4*[wg.GRID_VALUE_LONG,]+4*[wg.GRID_VALUE_FLOAT+':10,4',]+2*[wg.GRID_VALUE_FLOAT+':10,2',]
     1298    colLabels = ['H','K','L','mul','d','pos','sig','gam','Fosq','Fcsq','phase',]
     1299    Types = 4*[wg.GRID_VALUE_LONG,]+4*[wg.GRID_VALUE_FLOAT+':10,4',]+2*[wg.GRID_VALUE_FLOAT+':10,2',]+[wg.GRID_VALUE_FLOAT+':10,3',]
    13001300    self.PeakTable = G2gd.Table(refList,rowLabels=rowLabels,colLabels=colLabels,types=Types)
    13011301    self.dataFrame.SetLabel('Reflection List for '+self.RefList)
  • trunk/GSASIIspc.py

    r372 r380  
    322322    h,k,l,f = Uniq
    323323    Uniq=np.array(zip(h[:Nuniq],k[:Nuniq],l[:Nuniq]))
     324    phi = f[:Nuniq]
    324325    Uniq = np.array(Uniq)
    325326   
    326     return iabsnt,2*mulp,Uniq       #include Friedel pairs in powder mulp
    327 
    328        
    329 def GenHKL(HKL,SGData,Friedel=False):          #gives wrong answers!
    330     '''
    331     Generate set of equivalent reflections
    332     input:
    333         HKL - [h,k,l]
    334         SGData - space group data obtained from SpcGroup
    335         Friedel = True to retain Friedel pairs for centrosymmetric case
    336     returns:
    337         iabsnt = True is reflection is forbidden by symmetry
    338         mulp = reflection multiplicity including Fridel pairs
    339         Uniq = numpy array of equivalent hkl in descending order of h,k,l
    340         Phs = numpy array of corresponding phase factors (multiples of 2pi)
    341        
    342     NB: only works on one HKL at a time -
    343         it can be made to do an array of HKLs - not any faster!
    344     '''
    345     iabsnt = False
    346     H = np.array(HKL)
    347     Inv = SGData['SGInv']
    348     Cen = SGData['SGCen'][1:]           #skip P
    349     # do cell centering extinction check
    350     for cen in Cen:
    351         iabsnt |= bool(int(round(np.sum(cen*H*12)))%12)
    352     # get operators & expand if centrosymmetric
    353     Ops = SGData['SGOps']
    354     opM = np.array([op[0] for op in Ops])
    355     opT = np.array([op[1] for op in Ops])
    356     if Inv:
    357         opM = np.concatenate((opM,-opM))
    358         opT = np.concatenate((opT,-opT))
    359     # generate full hkl table and phase factors; find duplicates for multiplicity
    360     eqH = np.inner(opM,H)
    361     HT = np.dot(opT,H)%1.
    362     dup = len(ma.nonzero([ma.allclose(H,hkl,atol=0.001) for hkl in eqH])[0])
    363     mulp = len(eqH)/dup
    364     # generate unique reflection set (with or without Friedel pairs)
    365     HPT = [tuple(hp) for hp in np.column_stack((eqH,HT))]
    366     HPT.sort()
    367     UPT = HPT[::dup]
    368     UPT.reverse()
    369     if Inv and not Friedel:
    370         UPT = UPT[:len(UPT)/2]           
    371     Uniq = np.split(UPT,[3,4],axis=1)
    372     Phs = Uniq[1].T[0]
    373     Uniq = Uniq[0]
    374     # determine if extinct
    375     HPT = np.split(HPT,[3,4],axis=1)[1].T[0]
    376     HPT = np.reshape(HPT,(mulp,-1)).T
    377     HPT = np.around(HPT-HPT[0],4)%1.
    378     iabsnt |= np.any(HPT)
    379     return iabsnt,mulp,Uniq
    380                                    
    381 #def GenHKLs(H,SGData,Friedel=False):
    382 #    '''
    383 #    Generate set of equivalent reflections
    384 #    input:
    385 #        H - np.array of [h,k,l] assumed to exclude lattice centering extinctions
    386 #        SGData - space group data obtained from SpcGroup
    387 #        Friedel = True to retain Friedel pairs for centrosymmetric case
    388 #    returns:
    389 #        iabsnt = True is reflection is forbidden by symmetry
    390 #        mulp = reflection multiplicity including Fridel pairs
    391 #        Uniq = numpy array of equivalent hkl in descending order of h,k,l
    392 #        Phs = numpy array of corresponding phase factors (multiples of 2pi)
    393 #    this is not any faster than GenHKL above - oh well
    394 #    '''
    395 #    H = np.array(H,dtype=float)
    396 #    Inv = SGData['SGInv']
    397 #    Ops = SGData['SGOps']
    398 #    opM = np.array([op[0] for op in Ops])
    399 #    opT = np.array([op[1] for op in Ops])
    400 #    if Inv:
    401 #        opM = np.concatenate((opM,-opM))
    402 #        opT = np.concatenate((opT,-opT))
    403 #    # generate full hkl table and phase factors; find duplicates for multiplicity
    404 #    eqH = np.swapaxes(np.swapaxes(np.inner(opM,H),0,2),1,2)
    405 #    HeqH = zip(H,eqH)
    406 #    dup = np.array([sum([ma.allclose(h,hkl,atol=0.001) for hkl in HKL]) for h,HKL in HeqH])
    407 #    mulp = len(eqH[0])/dup
    408 #    # generate unique reflection set (with or without Friedel pairs)
    409 #    HT = np.dot(H,opT.T)%1.
    410 #    HPT = zip(mulp,dup,[[tuple(hp) for hp in HP] for HP in np.dstack((eqH,HT))])
    411 #    Uniq = []
    412 #    Phs = []
    413 #    iabsnt = []
    414 #    for mp,dp,hpt in HPT:
    415 #        hpt.sort()
    416 #        upt = hpt[::dp]
    417 #        upt.reverse()
    418 #        if Inv and not Friedel:
    419 #            upt = upt[:len(upt)/2]           
    420 #        uniq = np.split(upt,[3,4],axis=1)
    421 #        Phs.append(uniq[1].T[0])
    422 #        Uniq.append(uniq[0])
    423 #        # determine if extinct
    424 #        hpt = np.split(hpt,[3,4],axis=1)[1].T[0]
    425 #        hpt = np.reshape(hpt,(mp,-1)).T
    426 #        hpt = np.around(hpt-hpt[0],4)%1.
    427 #        iabsnt.append(np.any(hpt))
    428 #    return iabsnt,mulp,Uniq,Phs
    429                                    
     327    return iabsnt,2*mulp,Uniq,phi       #include Friedel pairs in powder mulp
     328                                 
    430329def GetOprPtrName(key):           
    431330    OprPtrName = {
  • trunk/GSASIIstruct.py

    r379 r380  
    1111import os.path as ospath
    1212import numpy as np
     13import numpy.linalg as nl
    1314import cPickle
    1415import time
     
    2627tand = lambda x: np.tan(x*np.pi/180.)
    2728asind = lambda x: 180.*np.arcsin(x)/np.pi
     29atan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
    2830
    2931
     
    287289        FFs = G2el.GetFormFactorCoeff(El.split('+')[0].split('-')[0])
    288290        for item in FFs:
    289             if item['Symbol'] == El:
     291            if item['Symbol'] == El.upper():
    290292                FFtable[El] = item
    291293    return FFtable
     
    343345    elif SGData['SGLaue'] in ['m3m','m3']:
    344346        return [pfx+'A0']
     347       
    345348           
    346349def GetPhaseData(PhaseData,Print=True):
    347350           
     351    def PrintFFtable(FFtable):
     352        print '\n Scattering factors:'
     353        print '   Symbol     fa                                      fb                                      fc'
     354        print 99*'-'
     355        for Ename in FFtable:
     356            ffdata = FFtable[Ename]
     357            fa = ffdata['fa']
     358            fb = ffdata['fb']
     359            print ' %8s %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f' %  \
     360                (Ename.ljust(8),fa[0],fa[1],fa[2],fa[3],fb[0],fb[1],fb[2],fb[3],ffdata['fc'])
     361               
     362    def PrintAtoms(General,Atoms):
     363        print '\n Atoms:'
     364        line = '   name    type  refine?   x         y         z    '+ \
     365            '  frac site sym  mult I/A   Uiso     U11     U22     U33     U12     U13     U23'
     366        if General['Type'] == 'magnetic':
     367            line += '   Mx     My     Mz'
     368        elif General['Type'] == 'macromolecular':
     369            line = ' res no  residue  chain '+line
     370        print line
     371        if General['Type'] == 'nuclear':
     372            print 135*'-'
     373            for i,at in enumerate(Atoms):
     374                line = '%7s'%(at[0])+'%7s'%(at[1])+'%7s'%(at[2])+'%10.5f'%(at[3])+'%10.5f'%(at[4])+ \
     375                    '%10.5f'%(at[5])+'%8.3f'%(at[6])+'%7s'%(at[7])+'%5d'%(at[8])+'%5s'%(at[9])
     376                if at[9] == 'I':
     377                    line += '%8.4f'%(at[10])+48*' '
     378                else:
     379                    line += 8*' '
     380                    for j in range(6):
     381                        line += '%8.4f'%(at[11+j])
     382                print line
     383       
     384       
    348385    if Print: print ' Phases:'
    349386    phaseVary = []
     
    351388    phaseConstr = {}
    352389    pawleyLookup = {}
     390    FFtables = {}
     391    Natoms = {}
     392    AtMults = {}
     393    AtIA = {}
    353394    for name in PhaseData:
    354395        General = PhaseData[name]['General']
    355396        pId = PhaseData[name]['pId']
    356397        pfx = str(pId)+'::'
     398        FFtable = GetFFtable(General)
     399        FFtables.update(FFtable)
    357400        Atoms = PhaseData[name]['Atoms']
    358401        try:
     
    373416            ' alpha =','%.3f'%(cell[4]),' beta =','%.3f'%(cell[5]),' gamma =', \
    374417            '%.3f'%(cell[6]),' volume =','%.3f'%(cell[7]),' Refine?',cell[0]
     418        if Print and FFtable:
     419            PrintFFtable(FFtable)
     420        Natoms[pfx] = 0
    375421        if Atoms:
    376             if Print: print '\n Atoms:'
    377             line = '   name    type  refine?   x         y         z    '+ \
    378                 '  frac site sym  mult I/A   Uiso     U11     U22     U33     U12     U13     U23'
    379             if General['Type'] == 'magnetic':
    380                 line += '   Mx     My     Mz'
    381             elif General['Type'] == 'macromolecular':
    382                 line = ' res no  residue  chain '+line
    383             if Print: print line
    384422            if General['Type'] == 'nuclear':
    385                 if Print: print 135*'-'
    386                 for at in Atoms:
    387                     line = '%7s'%(at[0])+'%7s'%(at[1])+'%7s'%(at[2])+'%10.5f'%(at[3])+'%10.5f'%(at[4])+ \
    388                         '%10.5f'%(at[5])+'%8.3f'%(at[6])+'%7s'%(at[7])+'%5d'%(at[8])+'%5s'%(at[9])
     423                Natoms[pfx] = len(Atoms)
     424                for i,at in enumerate(Atoms):
     425                    phaseDict.update({pfx+'Atype:'+str(i):at[1],pfx+'Afrac:'+str(i):at[6],pfx+'Amul:'+str(i):at[8],
     426                        pfx+'Ax:'+str(i):at[3],pfx+'Ay:'+str(i):at[4],pfx+'Az:'+str(i):at[5],
     427                        pfx+'dAx:'+str(i):0.,pfx+'dAy:'+str(i):0.,pfx+'dAz:'+str(i):0.,         #refined shifts for x,y,z
     428                        pfx+'AI/A:'+str(i):at[9],})
    389429                    if at[9] == 'I':
    390                         line += '%8.4f'%(at[10])+48*' '
     430                        phaseDict[pfx+'AUiso:'+str(i)] = at[10]
    391431                    else:
    392                         line += 8*' '
    393                         for i in range(6):
    394                             line += '%8.4f'%(at[11+i])
    395                     if Print: print line
     432                        phaseDict.update({pfx+'AU11:'+str(i):at[11],pfx+'AU22:'+str(i):at[12],pfx+'AU33:'+str(i):at[13],
     433                            pfx+'AU12:'+str(i):at[14],pfx+'AU13:'+str(i):at[15],pfx+'AU23:'+str(i):at[16]})
     434                    if 'F' in at[2]:
     435                        phaseVary.append(pfx+'Afrac:'+str(i))
    396436                    if 'X' in at[2]:
    397437                        xId,xCoef = G2spc.GetCSxinel(at[7])
     438                        delnames = [pfx+'dAx:'+str(i),pfx+'dAy:'+str(i),pfx+'dAz:'+str(i)]
     439                        for j in range(3):
     440                            if xId[j] > 0:                               
     441                                phaseVary.append(delnames[j])
     442                                for k in range(j):
     443                                    if xId[j] == xId[k]:
     444                                        G2mv.StoreEquivalence(delnames[k],((delnames[j],xCoef[j]),))
    398445                    if 'U' in at[2]:
    399                         uId,uCoef = G2spc.GetCSuinel(at[7])[:2]
     446                        if at[9] == 'I':
     447                            phaseVary.append(pfx+'AUiso:'+str(i))
     448                        else:
     449                            uId,uCoef = G2spc.GetCSuinel(at[7])[:2]
     450                            names = [pfx+'AU11:'+str(i),pfx+'AU22:'+str(i),pfx+'AU33:'+str(i),
     451                                pfx+'AU12:'+str(i),pfx+'AU13:'+str(i),pfx+'AU23:'+str(i)]
     452                            for j in range(6):
     453                                if uId[j] > 0:                               
     454                                    phaseVary.append(names[j])
     455                                    for k in range(j):
     456                                        if uId[j] == uId[k]:
     457                                            G2mv.StoreEquivalence(names[k],((names[j],uCoef[j]),))
     458            if Print:
     459                PrintAtoms(General,Atoms)
    400460#        elif General['Type'] == 'magnetic':
    401461#        elif General['Type'] == 'macromolecular':
    402 #       PWDR: harmonics texture, unit cell, atom parameters
     462#       PWDR: harmonics texture
    403463        elif PawleyRef:
    404464            pawleyVary = []
     
    411471            phaseVary += pawleyVary
    412472               
    413     return phaseVary,phaseDict,pawleyLookup
     473    return Natoms,phaseVary,phaseDict,pawleyLookup,FFtables
    414474   
    415475def SetPhaseData(parmDict,sigDict,Phases):
     
    455515            sigA = [sigDict[pfx+'A0'],0,0,0,0,0]
    456516        return A,sigA
     517       
     518    def PrintAtomsAndSig(General,Atoms,atomsSig):
     519        print '\n Atoms:'
     520        line = '   name      x         y         z      frac   Uiso     U11     U22     U33     U12     U13     U23'
     521        if General['Type'] == 'magnetic':
     522            line += '   Mx     My     Mz'
     523        elif General['Type'] == 'macromolecular':
     524            line = ' res no  residue  chain '+line
     525        print line
     526        if General['Type'] == 'nuclear':
     527            print 135*'-'
     528            fmt = {0:'%7s',1:'%7s',3:'%10.5f',4:'%10.5f',5:'%10.5f',6:'%8.3f',10:'%8.5f',
     529                11:'%8.5f',12:'%8.5f',13:'%8.5f',14:'%8.5f',15:'%8.5f',16:'%8.5f'}
     530            noFXsig = {3:[10*' ','%10s'],4:[10*' ','%10s'],5:[10*' ','%10s'],6:[8*' ','%8s']}
     531            for i,at in enumerate(Atoms):
     532                name = fmt[0]%(at[0])+fmt[1]%(at[1])+':'
     533                valstr = ' values:'
     534                sigstr = ' sig   :'
     535                for ind in [3,4,5,6]:
     536                    sigind = str(i)+':'+str(ind)
     537                    valstr += fmt[ind]%(at[ind])                   
     538                    if sigind in atomsSig:
     539                        sigstr += fmt[ind]%(atomsSig[sigind])
     540                    else:
     541                        sigstr += noFXsig[ind][1]%(noFXsig[ind][0])
     542                if at[9] == 'I':
     543                    valstr += fmt[10]%(at[10])
     544                    if str(i)+':10' in atomsSig:
     545                        sigstr += fmt[10]%(atomsSig[str(i)+':10'])
     546                    else:
     547                        sigstr += 8*' '
     548                else:
     549                    valstr += 8*' '
     550                    sigstr += 8*' '
     551                    for ind in [11,12,13,14,15,16]:
     552                        sigind = str(i)+':'+str(ind)
     553                        valstr += fmt[ind]%(at[ind])
     554                        if sigind in atomsSig:                       
     555                            sigstr += fmt[ind]%(atomsSig[sigind])
     556                        else:
     557                            sigstr += 8*' '
     558                print name
     559                print valstr
     560                print sigstr
    457561           
    458562    print '\n Phases:'
     
    462566        General = Phase['General']
    463567        SGData = General['SGData']
     568        Atoms = Phase['Atoms']
    464569        cell = General['Cell']
    465570        pId = Phase['pId']
     
    497602                else:
    498603                    refl[7] = 0
     604        else:
     605            atomsSig = {}
     606            if General['Type'] == 'nuclear':
     607                for i,at in enumerate(Atoms):
     608                    names = {3:pfx+'Ax:'+str(i),4:pfx+'Ay:'+str(i),5:pfx+'Az:'+str(i),6:pfx+'Afrac:'+str(i),
     609                        10:pfx+'AUiso:'+str(i),11:pfx+'AU11:'+str(i),12:pfx+'AU22:'+str(i),13:pfx+'AU33:'+str(i),
     610                        14:pfx+'AU12:'+str(i),15:pfx+'AU13:'+str(i),16:pfx+'AU23:'+str(i)}
     611                    for ind in [3,4,5,6]:
     612                        at[ind] = parmDict[names[ind]]
     613                        if ind in [3,4,5]:
     614                            name = names[ind].replace('A','dA')
     615                        else:
     616                            name = names[ind]
     617                        if name in sigDict:
     618                            atomsSig[str(i)+':'+str(ind)] = sigDict[name]
     619                    if at[9] == 'I':
     620                        at[10] = parmDict[names[10]]
     621                        if names[10] in sigDict:
     622                            atomsSig[str(i)+':10'] = sigDict[names[10]]
     623                    else:
     624                        for ind in [11,12,13,14,15,16]:
     625                            at[ind] = parmDict[names[ind]]
     626                            if names[ind] in sigDict:
     627                                atomsSig[str(i)+':'+str(ind)] = sigDict[names[ind]]
     628            PrintAtomsAndSig(General,Atoms,atomsSig)   
    499629
    500630def GetHistogramPhaseData(Phases,Histograms,Print=True):
     
    659789            refList = []
    660790            for h,k,l,d in HKLd:
    661                 ext,mul,Uniq = G2spc.GenHKLf([h,k,l],SGData)
     791                ext,mul,Uniq,phi = G2spc.GenHKLf([h,k,l],SGData)
    662792                if ext:
    663793                    continue
     
    665795                    pos = 2.0*asind(wave/(2.0*d))
    666796                    if limits[0] < pos < limits[1]:
    667                         refList.append([h,k,l,mul,d,pos,0.0,0.0,0.0,0.0,Uniq])
     797                        refList.append([h,k,l,mul,d,pos,0.0,0.0,0.0,0.0,0.0,Uniq,phi])
    668798                else:
    669799                    raise ValueError
     
    9381068        Type,instDict,insVary = GetInstParms(hId,Inst)
    9391069        controlDict[pfx+'histType'] = Type
     1070        if pfx+'Lam1' in instDict:
     1071            controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam1']
     1072        else:
     1073            controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam']           
    9401074        histDict.update(instDict)
    9411075        histVary += insVary
     
    9911125        elif 'Debye' in Sample['Type']:        #Debye-Scherrer
    9921126            sampSig = [0 for i in range(4)]
    993             for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
     1127            for i,item in enumerate(['Scale','Absorption','DisplaceX','DisplaceY']):
    9941128                Sample[item][0] = parmDict[pfx+item]
    9951129                if pfx+item in sigDict:
     
    10361170                ptlbls += '%14s'%(item)
    10371171                ptstr += '%14.4f'%(Sample[item][0])
    1038                 sigstr += '%14.4f'%(sampleSig[i])
     1172                if sampleSig[i]:
     1173                    sigstr += '%14.4f'%(sampleSig[i])
     1174                else:
     1175                    sigstr += 14*' '
    10391176           
    10401177        elif 'Debye' in Sample['Type']:        #Debye-Scherrer
     
    10421179                ptlbls += '%14s'%(item)
    10431180                ptstr += '%14.4f'%(Sample[item][0])
    1044                 sigstr += '%14.4f'%(sampleSig[i])
     1181                if sampleSig[i]:
     1182                    sigstr += '%14.4f'%(sampleSig[i])
     1183                else:
     1184                    sigstr += 14*' '
    10451185
    10461186        print ptlbls
     
    10691209            PrintBackgroundSig(Background,backSig)
    10701210
    1071 
    1072 def SetupSFcalc(General,Atoms):
    1073     ''' setup data for use in StructureFactor; mostly rearranging arrays
     1211def GetAtomFXU(pfx,FFtables,calcControls,parmDict):
     1212    Natoms = calcControls['Natoms'][pfx]
     1213    Tdata = Natoms*[' ',]
     1214    Mdata = np.zeros(Natoms)
     1215    IAdata = Natoms*[' ',]
     1216    Fdata = np.zeros(Natoms)
     1217    FFdata = []
     1218    Xdata = np.zeros((3,Natoms))
     1219    dXdata = np.zeros((3,Natoms))
     1220    Uisodata = np.zeros(Natoms)
     1221    Uijdata = np.zeros((6,Natoms))
     1222    keys = {'Atype:':Tdata,'Amul:':Mdata,'Afrac:':Fdata,'AI/A:':IAdata,
     1223        'dAx:':dXdata[0],'dAy:':dXdata[1],'dAz:':dXdata[2],
     1224        'Ax:':Xdata[0],'Ay:':Xdata[1],'Az:':Xdata[2],'AUiso:':Uisodata,
     1225        'AU11:':Uijdata[0],'AU22:':Uijdata[1],'AU33:':Uijdata[2],
     1226        'AU12:':Uijdata[3],'AU13:':Uijdata[4],'AU23:':Uijdata[5]}
     1227    for iatm in range(Natoms):
     1228        for key in keys:
     1229            parm = pfx+key+str(iatm)
     1230            if parm in parmDict:
     1231                keys[key][iatm] = parmDict[parm]
     1232        FFdata.append(FFtables[Tdata[iatm]])
     1233    return FFdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata
     1234   
     1235def StructureFactor(refList,G,hfx,pfx,SGData,calcControls,parmDict):
     1236    ''' Compute structure factors for all h,k,l for phase
    10741237    input:
    1075         General = dictionary of general phase info.
    1076         Atoms = list of atom parameters
    1077     returns:
    1078         G = reciprocal metric tensor
    1079         StrData = list of arrays; one entry per atom:
    1080             T = atom types
    1081             R = refinement flags, e.g. 'FXU'
    1082             F = site fractions
    1083             X = atom coordinates as numpy array
    1084             M = atom multiplicities
    1085             IA = isotropic/anisotropic thermal flags
    1086             Uiso = isotropic thermal parameters if IA = 'I'; else = 0
    1087             Uij = numpy array of 6 anisotropic thermal parameters if IA='A'; else = zeros
    1088     '''           
    1089     SGData = General['SGData']
    1090     Cell = General['Cell']
    1091     G,g = G2lat.cell2Gmat(Cell[1:7])        #skip refine & volume; get recip & real metric tensors
    1092     cx,ct,cs,cia = General['AtomPtrs']
    1093     X = [];F = [];T = [];IA = [];Uiso = [];Uij = [];R = [];M = []
    1094     for atom in Atoms:
    1095         T.append(atom[ct])
    1096         R.append(atom[ct+1])
    1097         F.append(atom[cx+3])
    1098         X.append(np.array(atom[cx:cx+3]))
    1099         M.append(atom[cia-1])
    1100         IA.append(atom[cia])
    1101         Uiso.append(atom[cia+1])
    1102         Uij.append(np.array(atom[cia+2:cia+8]))
    1103     return G,[T,R,np.array(F),np.array(X),np.array(M),IA,np.array(Uiso),np.array(Uij)]
    1104            
    1105 def StructureFactor(H,G,SGData,StrData,FFtable):
    1106     ''' Compute structure factor for a single h,k,l
    1107     H = np.array(h,k,l)
    1108     G = reciprocal metric tensor
    1109     SGData = space group info. dictionary output from SpcGroup
    1110     StrData = [
    1111         [atom types],
    1112         [refinement flags],
    1113         [site fractions],
    1114         np.array(coordinates),
    1115         [multiplicities],
    1116         [I/A flag],
    1117         [Uiso],
    1118         np.array(Uij)]
    1119     FFtable = dictionary of form factor coeff. for atom types used in StrData
    1120     '''
    1121     twopi = 2.0*math.pi
    1122     twopisq = 2.0*math.pi**2
    1123     SQ = G2lat.calc_rDsq2(H,G)/4.0          # SQ = (sin(theta)/lambda)**2
    1124     SQfactor = 8.0*SQ*math.pi**2
    1125     print 'SQ',SQfactor
    1126     FF = {}
    1127     for type in FFtable.keys():
    1128         FF[type] = G2el.ScatFac(FFtable[type],SQ)
    1129     print 'FF',FF
    1130     iabsnt,mulp,Uniq,Phs = G2spc.GenHKL(H,SGData,False)
    1131     fa = [0,0,0]        #real
    1132     fb = [0,0,0]        #imaginary
    1133     if not iabsnt:
    1134         phase = twopi*np.inner(Uniq,StrData[3])       #2pi[hx+ky+lz] for each atom in each equiv. hkl
     1238        refList: [ref] where each ref = h,k,l,m,d,...,[equiv h,k,l],phase[equiv]
     1239        G:      reciprocal metric tensor
     1240        pfx:    phase id string
     1241        SGData: space group info. dictionary output from SpcGroup
     1242        calcControls:
     1243        ParmDict:
     1244    puts result F^2 in each ref[8] in refList
     1245    '''       
     1246    twopi = 2.0*np.pi
     1247    twopisq = 2.0*np.pi**2
     1248    ast = np.sqrt(np.diag(G))
     1249    Mast = twopisq*np.multiply.outer(ast,ast)
     1250    FFtables = calcControls['FFtables']
     1251    FFdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,FFtables,calcControls,parmDict)
     1252    FP = np.array([El[hfx+'FP'] for El in FFdata])
     1253    FPP = np.array([El[hfx+'FPP'] for El in FFdata])
     1254    maxPos = len(SGData['SGOps'])
     1255    Uij = np.array(G2lat.U6toUij(Uijdata))
     1256    bij = Mast*Uij.T
     1257    for refl in refList:
     1258        fb = [0,0]
     1259        H = refl[:3]
     1260        SQ = 1./(2.*refl[4])**2
     1261        FF = np.array([G2el.ScatFac(El,SQ)[0] for El in FFdata])
     1262        SQfactor = 4.0*SQ*twopisq
     1263        Uniq = refl[11]
     1264        phi = refl[12]
     1265        phase = twopi*(np.inner(Uniq,(dXdata.T+Xdata.T))+phi[:,np.newaxis])
    11351266        sinp = np.sin(phase)
    11361267        cosp = np.cos(phase)
    1137         occ = StrData[2]*StrData[4]/mulp
    1138         Tiso = np.multiply(StrData[6],SQfactor)
    1139         Tuij = np.multiply(-SQfactor,np.inner(H,np.inner(G2spc.Uij2U(StrData[7]),H)))
    1140         print 'sinp',sinp
    1141         print 'cosp',cosp
    1142         print 'occ',occ
    1143         print 'Tiso',Tiso
    1144         print 'Tuij',Tuij
    1145     else:
    1146         print 'Absent'
     1268        occ = Mdata*Fdata/len(Uniq)
     1269        biso = -SQfactor*Uisodata
     1270        Tiso = np.where(biso<1.,np.exp(biso),1.0)
     1271        HbH = np.array([-np.inner(h,np.inner(bij,h)) for h in Uniq])
     1272        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
     1273        Tcorr = Tiso*Tuij
     1274        fa = np.array([(FF+FP)*occ*cosp*Tcorr,-FPP*occ*sinp*Tcorr])
     1275        fa = np.sum(np.sum(fa,axis=1),axis=1)        #real
     1276        if not SGData['SGInv']:
     1277            fb = np.array([(FF+FP)*occ*sinp*Tcorr,FPP*occ*cosp*Tcorr])
     1278            fb = np.sum(np.sum(fb,axis=1),axis=1)        #imaginary
     1279        refl[9] = fa[0]**2+fb[1]**2+fb[0]+fa[1]**2
     1280        refl[10] = atan2d(fb[0],fa[0])
     1281    return refList
     1282   
     1283def StructureFactorDerv(refList,G,hfx,pfx,SGData,calcControls,parmDict):
     1284    twopi = 2.0*np.pi
     1285    twopisq = 2.0*np.pi**2
     1286    ast = np.sqrt(np.diag(G))
     1287    Mast = twopisq*np.multiply.outer(ast,ast)
     1288    FFtables = calcControls['FFtables']
     1289    FFdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,FFtables,calcControls,parmDict)
     1290    FP = np.array([El[hfx+'FP'] for El in FFdata])
     1291    FPP = np.array([El[hfx+'FPP'] for El in FFdata])
     1292    maxPos = len(SGData['SGOps'])       
     1293    Uij = np.array(G2lat.U6toUij(Uijdata))
     1294    bij = Mast*Uij.T
     1295    dFdvDict = {}
     1296    dFdfr = np.zeros((len(refList),len(Mdata)))
     1297    dFdx = np.zeros((len(refList),len(Mdata),3))
     1298    dFdui = np.zeros((len(refList),len(Mdata)))
     1299    dFdua = np.zeros((len(refList),len(Mdata),6))
     1300    for iref,refl in enumerate(refList):
     1301        H = np.array(refl[:3])
     1302        SQ = 1./(2.*refl[4])**2          # or (sin(theta)/lambda)**2
     1303        FF = np.array([G2el.ScatFac(El,SQ)[0] for El in FFdata])
     1304        SQfactor = 8.0*SQ*np.pi**2
     1305        Uniq = refl[11]
     1306        phi = refl[12]
     1307        phase = twopi*(np.inner((dXdata.T+Xdata.T),Uniq)+phi[np.newaxis,:])
     1308        sinp = np.sin(phase)
     1309        cosp = np.cos(phase)
     1310        occ = Mdata*Fdata/len(Uniq)
     1311        biso = -SQfactor*Uisodata
     1312        Tiso = np.where(biso<1.,np.exp(biso),1.0)
     1313#        HbH = np.array([-np.inner(h,np.inner(bij,h)) for h in Uniq])
     1314        HbH = -np.inner(H,np.inner(bij,H))
     1315        Hij = np.array([Mast*np.multiply.outer(U,U) for U in Uniq])
     1316        Hij = np.array([G2lat.UijtoU6(Uij) for Uij in Hij])
     1317        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
     1318        Tcorr = Tiso*Tuij
     1319        fot = (FF+FP)*occ*Tcorr
     1320        fotp = FPP*occ*Tcorr
     1321        fa = np.array([fot[:,np.newaxis]*cosp,fotp[:,np.newaxis]*cosp])       #non positions
     1322        fb = np.array([fot[:,np.newaxis]*sinp,-fotp[:,np.newaxis]*sinp])
     1323       
     1324        fas = np.sum(np.sum(fa,axis=1),axis=1)
     1325        fbs = np.sum(np.sum(fb,axis=1),axis=1)
     1326        fax = np.array([-fot[:,np.newaxis]*sinp,-fotp[:,np.newaxis]*sinp])   #positions
     1327        fbx = np.array([fot[:,np.newaxis]*cosp,-fot[:,np.newaxis]*cosp])
     1328        #sum below is over Uniq
     1329        dfadfr = np.sum(fa/occ[:,np.newaxis],axis=2)
     1330        dfadx = np.sum(twopi*Uniq*fax[:,:,:,np.newaxis],axis=2)
     1331        dfadui = np.sum(-SQfactor*fa,axis=2)
     1332        dfadua = np.sum(-Hij*fa[:,:,:,np.newaxis],axis=2)
     1333        #NB: the above have been checked against PA(1:10,1:2) in strfctr.for     
     1334        dFdfr[iref] = 2.*(fas[0]*dfadfr[0]+fas[1]*dfadfr[1])*Mdata/len(Uniq)
     1335        dFdx[iref] = 2.*(fas[0]*dfadx[0]+fas[1]*dfadx[1])
     1336        dFdui[iref] = 2.*(fas[0]*dfadui[0]+fas[1]*dfadui[1])
     1337        dFdua[iref] = 2.*(fas[0]*dfadua[0]+fas[1]*dfadua[1])
     1338        if not SGData['SGInv']:
     1339            dfbdfr = np.sum(fb/occ[:,np.newaxis],axis=2)
     1340            dfbdx = np.sum(twopi*Uniq*fbx[:,:,:,np.newaxis],axis=2)         
     1341            dfbdui = np.sum(-SQfactor*fb,axis=2)
     1342            dfbdua = np.sum(-Hij*fb[:,:,:,np.newaxis],axis=2)
     1343            dFdfr[iref] += 2.*(fbs[0]*dfbdfr[0]-fbs[1]*dfbdfr[1])*Mdata/len(Uniq)
     1344            dFdx[iref] += 2.*(fbs[0]*dfbdx[0]+fbs[1]*dfbdx[1])
     1345            dFdui[iref] += 2.*(fbs[0]*dfbdui[0]-fbs[1]*dfbdui[1])
     1346            dFdua[iref] += 2.*(fbs[0]*dfbdua[0]+fbs[1]*dfbdua[1])
     1347        #loop over atoms - each dict entry is list of derivatives for all the reflections
     1348    for i in range(len(Mdata)):     
     1349        dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i]
     1350        dFdvDict[pfx+'dAx:'+str(i)] = dFdx.T[0][i]
     1351        dFdvDict[pfx+'dAy:'+str(i)] = dFdx.T[1][i]
     1352        dFdvDict[pfx+'dAz:'+str(i)] = dFdx.T[2][i]
     1353        dFdvDict[pfx+'AUiso:'+str(i)] = dFdui.T[i]
     1354        dFdvDict[pfx+'AU11:'+str(i)] = dFdua.T[0][i]
     1355        dFdvDict[pfx+'AU22:'+str(i)] = dFdua.T[1][i]
     1356        dFdvDict[pfx+'AU33:'+str(i)] = dFdua.T[2][i]
     1357        dFdvDict[pfx+'AU12:'+str(i)] = 2.*dFdua.T[3][i]
     1358        dFdvDict[pfx+'AU13:'+str(i)] = 2.*dFdua.T[4][i]
     1359        dFdvDict[pfx+'AU23:'+str(i)] = 2.*dFdua.T[5][i]
     1360    return dFdvDict
    11471361       
    11481362def Dict2Values(parmdict, varylist):
     
    11551369    values corresponding to keys in varylist'''
    11561370    parmdict.update(zip(varylist,values))
     1371   
     1372def ApplyXYZshifts(parmDict,varyList):
     1373    ''' takes atom x,y,z shift and applies it to corresponding atom x,y,z value
     1374    '''
     1375    for vary in varyList:
     1376        if 'dA' in vary:
     1377            parm = ''.join(vary.split('d'))
     1378            parmDict[parm] += parmDict[vary]
    11571379   
    11581380def SHPOcal(refl,g,phfx,hfx,SGData,calcControls,parmDict):
     
    12041426        MDAxis = calcControls[phfx+'MDAxis']
    12051427        sumMD = 0
    1206         for H in refl[10]:           
     1428        for H in refl[11]:           
    12071429            cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G)
    12081430            A = 1.0/np.sqrt((MD*cosP)**2+sinP**2/MD)
    12091431            sumMD += A**3
    1210         POcorr = sumMD/len(refl[10])
     1432        POcorr = sumMD/len(refl[11])
    12111433    else:   #spherical harmonics
    12121434        POcorr = SHPOcal(refl,g,phfx,hfx,SGData,calcControls,parmDict)
     
    12201442        sumMD = 0
    12211443        sumdMD = 0
    1222         for H in refl[10]:           
     1444        for H in refl[11]:           
    12231445            cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G)
    12241446            A = 1.0/np.sqrt((MD*cosP)**2+sinP**2/MD)
    12251447            sumMD += A**3
    12261448            sumdMD -= (1.5*A**5)*(2.0*MD*cosP**2-(sinP/MD)**2)
    1227         POcorr = sumMD/len(refl[10])
    1228         POderv[phfx+'MD'] = sumdMD/len(refl[10])
     1449        POcorr = sumMD/len(refl[11])
     1450        POderv[phfx+'MD'] = sumdMD/len(refl[11])
    12291451    else:   #spherical harmonics
    12301452        POcorr,POderv = SHPOcalDerv(refl,g,phfx,hfx,SGData,calcControls,parmDict)
     
    14131635        dDijDict[item] *= refl[4]**2*tand(refl[5]/2.0)
    14141636    return dDijDict
     1637   
     1638def GetFprime(controlDict,Histograms):
     1639    FFtables = controlDict['FFtables']
     1640    if not FFtables:
     1641        return
     1642    for histogram in Histograms:
     1643        if 'PWDR' in histogram[:4]:
     1644            Histogram = Histograms[histogram]
     1645            hId = Histogram['hId']
     1646            hfx = ':%d:'%(hId)
     1647            keV = controlDict[hfx+'keV']
     1648            for El in FFtables:
     1649                Orbs = G2el.GetXsectionCoeff(El.split('+')[0].split('-')[0])
     1650                FP,FPP,Mu = G2el.FPcalc(Orbs, keV)
     1651                FFtables[El][hfx+'FP'] = FP
     1652                FFtables[El][hfx+'FPP'] = FPP               
    14151653           
    14161654def getPowderProfile(parmDict,x,varylist,Histogram,Phases,calcControls,pawleyLookup):
     
    14261664        sig = max(0.001,sig)
    14271665        gam = X/cosd(refl[5]/2.0)+Y*tanPos+GetSampleGam(refl,wave,G,phfx,calcControls,parmDict,sizeEllipse) #save peak gamma
    1428         gam = max(0.01,gam)
     1666        gam = max(0.001,gam)
    14291667        return sig,gam
    14301668               
     
    14361674       
    14371675    if 'C' in calcControls[hfx+'histType']:   
    1438         shl = max(parmDict[hfx+'SH/L'],0.0005)
     1676        shl = max(parmDict[hfx+'SH/L'],0.002)
    14391677        Ka2 = False
    14401678        if hfx+'Lam1' in parmDict.keys():
     
    14541692        pfx = '%d::'%(pId)
    14551693        phfx = '%d:%d:'%(pId,hId)
     1694        hfx = ':%d:'%(hId)
    14561695        SGData = Phase['General']['SGData']
    14571696        A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
    14581697        G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
     1698        Vst = np.sqrt(nl.det(G))    #V*
     1699        if 'Pawley' not in Phase['General']['Type']:
     1700            refList = StructureFactor(refList,G,hfx,pfx,SGData,calcControls,parmDict)
    14591701        sizeEllipse = []
    14601702        if calcControls[phfx+'SizeType'] == 'ellipsoidal':
     
    14641706                h,k,l = refl[:3]
    14651707                refl[5] = GetReflPos(refl,wave,G,hfx,calcControls,parmDict)         #corrected reflection position
     1708                Lorenz = 1./(2.*sind(refl[5]/2.)**2*cosd(refl[5]/2.))           #Lorentz correction
    14661709                refl[5] += GetHStrainShift(refl,SGData,phfx,parmDict)               #apply hydrostatic strain shift
    14671710                refl[6:8] = GetReflSIgGam(refl,wave,G,hfx,phfx,calcControls,parmDict,sizeEllipse)    #peak sig & gam
    1468                 Icorr = GetIntensityCorr(refl,G,g,phfx,hfx,SGData,calcControls,parmDict)
     1711                Icorr = GetIntensityCorr(refl,G,g,phfx,hfx,SGData,calcControls,parmDict)*Vst*Lorenz
    14691712                if 'Pawley' in Phase['General']['Type']:
    14701713                    try:
    1471                         refl[8] = abs(parmDict[pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])])
     1714                        refl[9] = abs(parmDict[pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])])
    14721715                    except KeyError:
    14731716#                        print ' ***Error %d,%d,%d missing from Pawley reflection list ***'%(h,k,l)
    14741717                        continue
    1475                 else:
    1476                     raise ValueError       #wants strctrfacr calc here
    14771718                Wd,fmin,fmax = G2pwd.getWidths(refl[5],refl[6],refl[7],shl)
    14781719                iBeg = np.searchsorted(x,refl[5]-fmin)
     
    14821723                elif not iBeg-iFin:     #peak above high limit - done
    14831724                    break
    1484                 yc[iBeg:iFin] += Icorr*refl[8]*G2pwd.getFCJVoigt3(refl[5],refl[6],refl[7],shl,x[iBeg:iFin])    #>90% of time spent here
     1725                yc[iBeg:iFin] += Icorr*refl[9]*G2pwd.getFCJVoigt3(refl[5],refl[6],refl[7],shl,x[iBeg:iFin])    #>90% of time spent here
    14851726                if Ka2:
    14861727                    pos2 = refl[5]+lamRatio*tand(refl[5]/2.0)       # + 360/pi * Dlam/lam * tan(th)
     
    14921733                    elif not iBeg-iFin:     #peak above high limit - done
    14931734                        return yc,yb
    1494                     yc[iBeg:iFin] += Icorr*refl[8]*kRatio*G2pwd.getFCJVoigt3(pos2,refl[6],refl[7],shl,x[iBeg:iFin])        #and here
     1735                    yc[iBeg:iFin] += Icorr*refl[9]*kRatio*G2pwd.getFCJVoigt3(pos2,refl[6],refl[7],shl,x[iBeg:iFin])        #and here
    14951736            elif 'T' in calcControls[hfx+'histType']:
    14961737                print 'TOF Undefined at present'
     
    15551796        A = [parmDict[pfx+'A%d'%(i)] for i in range(6)]
    15561797        G,g = G2lat.A2Gmat(A)       #recip & real metric tensors
     1798        Vst = np.sqrt(nl.det(G))    #V*
     1799        if 'Pawley' not in Phase['General']['Type']:
     1800            dFdvDict = StructureFactorDerv(refList,G,hfx,pfx,SGData,calcControls,parmDict)
    15571801        sizeEllipse = []
    15581802        if calcControls[phfx+'SizeType'] == 'ellipsoidal':
    15591803            sizeEllipse = G2lat.U6toUij([parmDIct[phfx+'Size:%d'%(i)] for i in range(6)])
    1560         for refl in refList:
     1804        for iref,refl in enumerate(refList):
    15611805            if 'C' in calcControls[hfx+'histType']:        #CW powder
    15621806                h,k,l = refl[:3]
    1563                 iref = pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)]
     1807                Lorenz = 1./(2.*sind(refl[5]/2.)**2*cosd(refl[5]/2.))           #Lorentz correction
    15641808                Icorr,dIdsh,dIdsp,dIdpola,dIdPO = GetIntensityDerv(refl,G,g,phfx,hfx,SGData,calcControls,parmDict)
     1809                Icorr *= Vst*Lorenz
    15651810                if 'Pawley' in Phase['General']['Type']:
    15661811                    try:
    1567                         refl[8] = abs(parmDict[pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])])
     1812                        refl[9] = abs(parmDict[pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])])
    15681813                    except KeyError:
    15691814#                        print ' ***Error %d,%d,%d missing from Pawley reflection list ***'%(h,k,l)
    15701815                        continue
    1571                 else:
    1572                     raise ValueError       #wants strctrfacr deriv calc here
    15731816                Wd,fmin,fmax = G2pwd.getWidths(refl[5],refl[6],refl[7],shl)
    15741817                iBeg = np.searchsorted(x,refl[5]-fmin)
     
    15841827                dMdipk = G2pwd.getdFCJVoigt3(refl[5],refl[6],refl[7],shl,x[iBeg:iFin])
    15851828                for i in range(1,5):
    1586                     dMdpk[i][iBeg:iFin] += 100.*dx*Icorr*refl[8]*dMdipk[i]
    1587                 dMdpk[0][iBeg:iFin] += 100.*dx*Icorr*refl[8]*dMdipk[0]
     1829                    dMdpk[i][iBeg:iFin] += 100.*dx*Icorr*refl[9]*dMdipk[i]
     1830                dMdpk[0][iBeg:iFin] += 100.*dx*Icorr*refl[9]*dMdipk[0]
    15881831                dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4]}
    15891832                if Ka2:
     
    15951838                        dMdipk2 = G2pwd.getdFCJVoigt3(pos2,refl[6],refl[7],shl,x[iBeg2:iFin2])
    15961839                        for i in range(1,5):
    1597                             dMdpk[i][iBeg2:iFin2] += 100.*dx*Icorr*refl[8]*kRatio*dMdipk2[i]
    1598                         dMdpk[0][iBeg2:iFin2] += 100.*dx*Icorr*refl[8]*kRatio*dMdipk2[0]
     1840                            dMdpk[i][iBeg2:iFin2] += 100.*dx*Icorr*refl[9]*kRatio*dMdipk2[i]
     1841                        dMdpk[0][iBeg2:iFin2] += 100.*dx*Icorr*refl[9]*kRatio*dMdipk2[0]
    15991842                        dMdpk[5][iBeg2:iFin2] += 100.*dx*Icorr*dMdipk2[0]
    1600                         dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':dMdpk[5]*refl[8]}
    1601                 try:
    1602                     idx = varylist.index(pfx+'PWLref:'+str(iref))
    1603                     dMdv[idx] = dervDict['int']/refl[8]
    1604                 except ValueError:
    1605                     pass
     1843                        dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':dMdpk[5]*refl[9]}
     1844                if 'Pawley' in Phase['General']['Type']:
     1845                    try:
     1846                        idx = varylist.index(pfx+'PWLref:'+str(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)]))
     1847                        dMdv[idx] = dervDict['int']/refl[9]
     1848                    except ValueError:
     1849                        pass
    16061850                dpdA,dpdw,dpdZ,dpdSh,dpdTr,dpdX,dpdY = GetReflPosDerv(refl,wave,A,hfx,calcControls,parmDict)
    16071851                names = {hfx+'Scale':[dIdsh,'int'],hfx+'Polariz.':[dIdpola,'int'],phfx+'Scale':[dIdsp,'int'],
     
    16291873                for name in gamDict:
    16301874                    if name in varylist:
    1631                         dMdv[varylist.index(name)] += gamDict[name]*dervDict['gam']                       
     1875                        dMdv[varylist.index(name)] += gamDict[name]*dervDict['gam']
     1876                                               
    16321877            elif 'T' in calcControls[hfx+'histType']:
    16331878                print 'TOF Undefined at present'
    1634                 raise ValueError    #no TOF yet
     1879                raise Exception    #no TOF yet
     1880#do atom derivatives -  for F,X & U so far             
     1881            atNames = [pfx+'Afrac',pfx+'dAx',pfx+'dAy',pfx+'dAz',pfx+'AUiso',
     1882                pfx+'AU11',pfx+'AU22',pfx+'AU33',pfx+'AU12',pfx+'AU13',pfx+'AU23']
     1883            for atname in atNames:
     1884                for name in varylist:
     1885                    if atname in name:
     1886                        dMdv[varylist.index(name)] += dFdvDict[name][iref]*dervDict['int']/refl[9]
    16351887           
    16361888    return dMdv   
     
    17041956    varyList = []
    17051957    parmDict = {}
    1706     calcControls = {}   
     1958    calcControls = {}
     1959    G2mv.InitVars()   
    17071960    Controls = GetControls(GPXfile)
    17081961    ShowControls(Controls)           
     
    17161969        print ' *** Refine aborted ***'
    17171970        raise Exception
    1718     phaseVary,phaseDict,pawleyLookup = GetPhaseData(Phases)
     1971    Natoms,phaseVary,phaseDict,pawleyLookup,FFtables = GetPhaseData(Phases)
     1972    calcControls['Natoms'] = Natoms
     1973    calcControls['FFtables'] = FFtables
    17191974    hapVary,hapDict,controlDict = GetHistogramPhaseData(Phases,Histograms)
    17201975    calcControls.update(controlDict)
     
    17251980    parmDict.update(hapDict)
    17261981    parmDict.update(histDict)
     1982    GetFprime(calcControls,Histograms)
    17271983    constrDict,constrFlag,fixedList = G2mv.InputParse([])        #constraints go here?
    17281984    groups,parmlist = G2mv.GroupConstraints(constrDict)
     
    17442000                args=([Histograms,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg))
    17452001            ncyc = int(result[2]['nfev']/len(varyList))
    1746 #        table = dict(zip(varyList,zip(values,result[0],(result[0]-values))))
    1747 #        for item in table: print item,table[item]               #useful debug - are things shifting?
     2002        table = dict(zip(varyList,zip(values,result[0],(result[0]-values))))
     2003        for item in table: print item,table[item]               #useful debug - are things shifting?
    17482004        runtime = time.time()-begin
    17492005        chisq = np.sum(result[2]['fvec']**2)
    17502006        Values2Dict(parmDict, varyList, result[0])
     2007        ApplyXYZshifts(parmDict,varyList)
    17512008        G2mv.Dict2Map(parmDict)
     2009       
    17522010        Rwp = np.sqrt(chisq/Histograms['sumwYo'])*100.      #to %
    17532011        GOF = chisq/(Histograms['Nobs']-len(varyList))
Note: See TracChangeset for help on using the changeset viewer.