Changeset 395 for trunk/GSASIIstruct.py


Ignore:
Timestamp:
Oct 20, 2011 9:25:32 AM (10 years ago)
Author:
vondreele
Message:

Add goniometer omega, chi & phi to sample data
put SH texture in General
fix phase delete to remove it from reflection lists as well
continue development of constraints/restraints GUI
fixes to texture computations, GUI & least squares refinement

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIIstruct.py

    r391 r395  
    382382                print line
    383383       
     384    def PrintTexture(textureData):
     385        print '\n Spherical harmonics texture: Order:' + \
     386            str(textureData['Order'])+' Refine? '+str(textureData['SH Coeff'][0])
     387        names = ['omega','chi','phi']
     388        line = '\n'
     389        for name in names:
     390            line += ' SH '+name+':'+'%12.4f'%(textureData['Sample '+name][1])+' Refine? '+str(textureData['Sample '+name][0])
     391        print line
     392        print '\n Texture coefficients:'
     393        ptlbls = ' names :'
     394        ptstr =  ' values:'
     395        SHcoeff = textureData['SH Coeff'][1]
     396        for item in SHcoeff:
     397            ptlbls += '%12s'%(item)
     398            ptstr += '%12.4f'%(SHcoeff[item])
     399        print ptlbls
     400        print ptstr   
    384401       
    385402    if Print: print ' Phases:'
     
    392409    AtMults = {}
    393410    AtIA = {}
     411    shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
     412    SamSym = dict(zip(shModels,['0','-1','2/m','mmm']))
    394413    for name in PhaseData:
    395414        General = PhaseData[name]['General']
     
    456475                                        if uId[j] == uId[k]:
    457476                                            G2mv.StoreEquivalence(names[k],((names[j],uCoef[j]),))
     477#            elif General['Type'] == 'magnetic':
     478#            elif General['Type'] == 'macromolecular':
     479
     480            if 'SH Texture' in General:
     481                textureData = General['SH Texture']
     482                phaseDict[pfx+'SHmodel'] = SamSym[textureData['Model']]
     483                phaseDict[pfx+'SHorder'] = textureData['Order']
     484                for name in ['omega','chi','phi']:
     485                    phaseDict[pfx+'SH '+name] = textureData['Sample '+name][1]
     486                    if textureData['Sample '+name][0]:
     487                        phaseVary.append(pfx+'SH '+name)
     488                for name in textureData['SH Coeff'][1]:
     489                    phaseDict[pfx+name] = textureData['SH Coeff'][1][name]
     490                    if textureData['SH Coeff'][0]:
     491                        phaseVary.append(pfx+name)
     492               
    458493            if Print:
    459494                PrintAtoms(General,Atoms)
    460 #        elif General['Type'] == 'magnetic':
    461 #        elif General['Type'] == 'macromolecular':
    462 #       PWDR: harmonics texture
     495                if 'SH Texture' in General:
     496                    PrintTexture(textureData)
     497                   
    463498        elif PawleyRef:
    464499            pawleyVary = []
     
    627662                print valstr
    628663                print sigstr
     664               
     665    def PrintSHtextureAndSig(textureData,SHtextureSig):
     666        print '\n Spherical harmonics texture: Order:' + str(textureData['Order'])
     667        names = ['omega','chi','phi']
     668        namstr = '  names :'
     669        ptstr =  '  values:'
     670        sigstr = '  esds  :'
     671        for name in names:
     672            namstr += '%12s'%(name)
     673            ptstr += '%12.3f'%(textureData['Sample '+name][1])
     674            if 'Sample '+name in SHtextureSig:
     675                sigstr += '%12.3f'%(SHtextureSig['Sample '+name])
     676            else:
     677                sigstr += 12*' '
     678        print namstr
     679        print ptstr
     680        print sigstr
     681        print '\n Texture coefficients:'
     682        namstr = '  names :'
     683        ptstr =  '  values:'
     684        sigstr = '  esds  :'
     685        SHcoeff = textureData['SH Coeff'][1]
     686        for name in SHcoeff:
     687            namstr += '%12s'%(name)
     688            ptstr += '%12.3f'%(SHcoeff[name])
     689            if name in SHtextureSig:
     690                sigstr += '%12.3f'%(SHtextureSig[name])
     691            else:
     692                sigstr += 12*' '
     693        print namstr
     694        print ptstr
     695        print sigstr
     696       
    629697           
    630698    print '\n Phases:'
     
    636704        Atoms = Phase['Atoms']
    637705        cell = General['Cell']
     706        textureData = General['SH Texture']
    638707        pId = Phase['pId']
    639708        pfx = str(pId)+'::'
     
    709778                            if names[ind] in sigDict:
    710779                                atomsSig[str(i)+':'+str(ind)] = sigDict[names[ind]]
    711             PrintAtomsAndSig(General,Atoms,atomsSig)   
     780            PrintAtomsAndSig(General,Atoms,atomsSig)
     781           
     782        if textureData['Order']:
     783            SHtextureSig = {}
     784            for name in ['omega','chi','phi']:
     785                aname = pfx+'SH '+name
     786                textureData['Sample '+name][1] = parmDict[aname]
     787                if aname in sigDict:
     788                    SHtextureSig['Sample '+name] = sigDict[aname]
     789            for name in textureData['SH Coeff'][1]:
     790                aname = pfx+name
     791                textureData['SH Coeff'][1][name] = parmDict[aname]
     792                if aname in sigDict:
     793                    SHtextureSig[name] = sigDict[aname]
     794            PrintSHtextureAndSig(textureData,SHtextureSig)
    712795
    713796def GetHistogramPhaseData(Phases,Histograms,Print=True):
     
    776859        for item in hapData[5]:
    777860            ptlbls += '%12s'%(item)
    778             ptstr += '%12.4f'%(hapData[5][item])
     861            ptstr += '%12.3f'%(hapData[5][item])
    779862        print ptlbls
    780863        print ptstr
     
    9671050        for item in hapData[5]:
    9681051            ptlbls += '%12s'%(item)
    969             ptstr += '%12.4f'%(hapData[5][item])
     1052            ptstr += '%12.3f'%(hapData[5][item])
    9701053            if item in POsig:
    971                 sigstr += '%12.4f'%(POsig[item])
     1054                sigstr += '%12.3f'%(POsig[item])
    9721055            else:
    9731056                sigstr += 12*' '
     
    10741157        sampVary = []
    10751158        hfx = ':'+str(hId)+':'       
    1076         sampDict = {hfx+'Gonio. radius':Sample['Gonio. radius']}
     1159        sampDict = {hfx+'Gonio. radius':Sample['Gonio. radius'],hfx+'Omega':Sample['Omega'],
     1160            hfx+'Chi':Sample['Chi'],hfx+'Phi':Sample['Phi']}
    10771161        Type = Sample['Type']
    10781162        if 'Bragg' in Type:             #Bragg-Brentano
     
    11161200    def PrintSampleParms(Sample):
    11171201        print '\n Sample Parameters:'
     1202        print ' Goniometer omega = %.2f, chi = %.2f, phi = %.2f'% \
     1203            (Sample['Omega'],Sample['Chi'],Sample['Phi'])
    11181204        ptlbls = ' name  :'
    11191205        ptstr =  ' value :'
     
    14671553            parmDict[parm] += parmDict[item]
    14681554   
     1555def SHTXcal(refl,g,pfx,hfx,SGData,calcControls,parmDict):
     1556    IFCoup = 'Bragg' in calcControls[hfx+'instType']
     1557    odfCor = 1.0
     1558    H = refl[:3]
     1559    cell = G2lat.Gmat2cell(g)
     1560    Sangls = [parmDict[pfx+'SH omega'],parmDict[pfx+'SH chi'],parmDict[pfx+'SH phi']]
     1561    Gangls = [parmDict[hfx+'Omega'],parmDict[hfx+'Chi'],parmDict[hfx+'Phi'],parmDict[hfx+'Azimuth']]
     1562    phi,beta = G2lat.CrsAng(H,cell,SGData)
     1563    psi,gam,x,x = G2lat.SamAng(refl[5]/2.,Gangls,Sangls,IFCoup) #ignore 2 sets of angle derivs.
     1564    SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],parmDict[pfx+'SHmodel'],parmDict[pfx+'SHorder'])
     1565    for item in SHnames:
     1566        L,M,N = eval(item.strip('C'))
     1567        Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
     1568        Ksl,x,x = G2lat.GetKsl(L,M,parmDict[pfx+'SHmodel'],psi,gam)
     1569        Lnorm = G2lat.Lnorm(L)
     1570        odfCor += parmDict[pfx+item]*Lnorm*Kcl*Ksl
     1571    return odfCor
     1572   
     1573def SHTXcalDerv(refl,g,pfx,hfx,SGData,calcControls,parmDict):
     1574    FORPI = 12.5663706143592
     1575    IFCoup = 'Bragg' in calcControls[hfx+'instType']
     1576    odfCor = 1.0
     1577    dFdODF = {}
     1578    dFdSA = [0,0,0]
     1579    H = refl[:3]
     1580    cell = G2lat.Gmat2cell(g)
     1581    Sangls = [parmDict[pfx+'SH omega'],parmDict[pfx+'SH chi'],parmDict[pfx+'SH phi']]
     1582    Gangls = [parmDict[hfx+'Omega'],parmDict[hfx+'Chi'],parmDict[hfx+'Phi'],parmDict[hfx+'Azimuth']]
     1583    phi,beta = G2lat.CrsAng(H,cell,SGData)
     1584    psi,gam,dPSdA,dGMdA = G2lat.SamAng(refl[5]/2.,Gangls,Sangls,IFCoup)
     1585    SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],parmDict[pfx+'SHmodel'],parmDict[pfx+'SHorder'])
     1586    for item in SHnames:
     1587        L,M,N = eval(item.strip('C'))
     1588        Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta)
     1589        Ksl,dKsdp,dKsdg = G2lat.GetKsl(L,M,parmDict[pfx+'SHmodel'],psi,gam)
     1590        Lnorm = G2lat.Lnorm(L)
     1591        odfCor += parmDict[pfx+item]*Lnorm*Kcl*Ksl
     1592        dFdODF[pfx+item] = Lnorm*Kcl*Ksl
     1593        for i in range(3):
     1594            dFdSA[i] += parmDict[pfx+item]*Lnorm*Kcl*(dKsdp*dPSdA[i]+dKsdg*dGMdA[i])
     1595    return odfCor,dFdODF,dFdSA
     1596   
    14691597def SHPOcal(refl,g,phfx,hfx,SGData,calcControls,parmDict):
    14701598    odfCor = 1.0
     
    14801608    phi,beta = G2lat.CrsAng(H,cell,SGData)
    14811609    psi,gam,x,x = G2lat.SamAng(refl[5]/2.,Gangls,Sangl,IFCoup) #ignore 2 sets of angle derivs.
    1482     SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],None,calcControls[phfx+'SHord'])
     1610    SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],'0',calcControls[phfx+'SHord'],False)
    14831611    for item in SHnames:
    14841612        L,N = eval(item.strip('C'))
     
    15021630    phi,beta = G2lat.CrsAng(H,cell,SGData)
    15031631    psi,gam,x,x = G2lat.SamAng(refl[5]/2.,Gangls,Sangl,IFCoup) #ignore 2 sets of angle derivs.
    1504     SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],None,calcControls[phfx+'SHord'])
     1632    SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],'0',calcControls[phfx+'SHord'],False)
    15051633    for item in SHnames:
    15061634        L,N = eval(item.strip('C'))
    15071635        Kcsl,Lnorm = G2lat.GetKclKsl(L,N,SGData['SGLaue'],psi,phi,beta)
    15081636        odfCor += parmDict[phfx+item]*Lnorm*Kcsl
    1509         dFdODF[phfx+item] = FORPI*Kcsl
     1637        dFdODF[phfx+item] = Kcsl*Lnorm
    15101638    return odfCor,dFdODF
    15111639   
     
    15421670    return POcorr,POderv
    15431671   
    1544 def GetIntensityCorr(refl,G,g,phfx,hfx,SGData,calcControls,parmDict):
     1672def GetIntensityCorr(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict):
    15451673    Icorr = parmDict[phfx+'Scale']*parmDict[hfx+'Scale']*refl[3]               #scale*multiplicity
    15461674    Icorr *= G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5],parmDict[hfx+'Azimuth'])[0]
    15471675    Icorr *= GetPrefOri(refl,G,g,phfx,hfx,SGData,calcControls,parmDict)
     1676    if pfx+'SHorder' in parmDict:
     1677        Icorr *= SHTXcal(refl,g,pfx,hfx,SGData,calcControls,parmDict)
    15481678    refl[13] = Icorr       
    15491679   
    1550 def GetIntensityDerv(refl,G,g,phfx,hfx,SGData,calcControls,parmDict):
     1680def GetIntensityDerv(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict):
     1681    dIdsh = 1./parmDict[hfx+'Scale']
     1682    dIdsp = 1./parmDict[phfx+'Scale']
    15511683    pola,dIdPola = G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5],parmDict[hfx+'Azimuth'])
     1684    dIdPola /= pola
    15521685    POcorr,dIdPO = GetPrefOriDerv(refl,G,g,phfx,hfx,SGData,calcControls,parmDict)
    1553     dIdPola /= pola
    15541686    for iPO in dIdPO:
    15551687        dIdPO[iPO] /= POcorr
    1556     dIdsh = 1./parmDict[hfx+'Scale']
    1557     dIdsp = 1./parmDict[phfx+'Scale']
    1558     return dIdsh,dIdsp,dIdPola,dIdPO
     1688    dFdODF = {}
     1689    dFdSA = [0,0,0]
     1690    if pfx+'SHorder' in parmDict:
     1691        odfCor,dFdODF,dFdSA = SHTXcalDerv(refl,g,pfx,hfx,SGData,calcControls,parmDict)
     1692        for iSH in dFdODF:
     1693            dFdODF[iSH] /= odfCor
     1694        for i in range(3):
     1695            dFdSA[i] /= odfCor
     1696    return dIdsh,dIdsp,dIdPola,dIdPO,dFdODF,dFdSA
    15591697       
    15601698def GetSampleGam(refl,wave,G,phfx,calcControls,parmDict,sizeEllipse):
     
    17971935                refl[5] += GetHStrainShift(refl,SGData,phfx,parmDict)               #apply hydrostatic strain shift
    17981936                refl[6:8] = GetReflSIgGam(refl,wave,G,hfx,phfx,calcControls,parmDict,sizeEllipse)    #peak sig & gam
    1799                 GetIntensityCorr(refl,G,g,phfx,hfx,SGData,calcControls,parmDict)    #puts corrections in refl[13]
     1937                GetIntensityCorr(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)    #puts corrections in refl[13]
    18001938                refl[13] *= Vst*Lorenz
    18011939                if 'Pawley' in Phase['General']['Type']:
     
    19512089            if 'C' in calcControls[hfx+'histType']:        #CW powder
    19522090                h,k,l = refl[:3]
    1953                 dIdsh,dIdsp,dIdpola,dIdPO = GetIntensityDerv(refl,G,g,phfx,hfx,SGData,calcControls,parmDict)
     2091                dIdsh,dIdsp,dIdpola,dIdPO,dFdODF,dFdSA = GetIntensityDerv(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
    19542092                if 'Pawley' in Phase['General']['Type']:
    19552093                    try:
     
    20062144                    if iPO in varylist:
    20072145                        dMdv[varylist.index(iPO)] += dIdPO[iPO]*dervDict['int']
     2146                for i,name in enumerate(['omega','chi','phi']):
     2147                    aname = pfx+'SH '+name
     2148                    if aname in varylist:
     2149                        dMdv[varylist.index(aname)] += dFdSA[i]*dervDict['int']
     2150                for iSH in dFdODF:
     2151                    if iSH in varylist:
     2152                        dMdv[varylist.index(iSH)] += dFdODF[iSH]*dervDict['int']
    20082153                cellDervNames = cellVaryDerv(pfx,SGData,dpdA)
    20092154                for name,dpdA in cellDervNames:
     
    20982243            GoOn = dlg.Update(Rwp,newmsg='%s%8.3f%s'%('Powder profile wRp =',Rwp,'%'))[0]
    20992244            if not GoOn:
    2100                 return -M           #abort!!
     2245                parmDict['saved values'] = values
     2246                raise Exception         #Abort!!
    21012247        return M
    21022248   
     
    21392285        Ftol = Controls['min dM/M']
    21402286        Factor = Controls['shift factor']
    2141         if Controls['deriv type'] == 'analytic':
    2142             result = so.leastsq(errRefine,values,Dfun=dervRefine,full_output=True,
    2143                 ftol=Ftol,col_deriv=True,factor=Factor,
    2144                 args=([Histograms,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg))
    2145             ncyc = int(result[2]['nfev']/2)               
    2146         else:           #'numeric'
    2147             result = so.leastsq(errRefine,values,full_output=True,ftol=Ftol,epsfcn=1.e-8,factor=Factor,
    2148                 args=([Histograms,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg))
    2149             ncyc = int(result[2]['nfev']/len(varyList))
    2150 #        table = dict(zip(varyList,zip(values,result[0],(result[0]-values))))
    2151 #        for item in table: print item,table[item]               #useful debug - are things shifting?
     2287        try:
     2288            if Controls['deriv type'] == 'analytic':
     2289                result = so.leastsq(errRefine,values,Dfun=dervRefine,full_output=True,
     2290                    ftol=Ftol,col_deriv=True,factor=Factor,
     2291                    args=([Histograms,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg))
     2292                ncyc = int(result[2]['nfev']/2)               
     2293            else:           #'numeric'
     2294                result = so.leastsq(errRefine,values,full_output=True,ftol=Ftol,epsfcn=1.e-8,factor=Factor,
     2295                    args=([Histograms,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg))
     2296                ncyc = int(result[2]['nfev']/len(varyList))
     2297#            table = dict(zip(varyList,zip(values,result[0],(result[0]-values))))
     2298#            for item in table: print item,table[item]               #useful debug - are things shifting?
     2299        except Exception:
     2300            result = [parmDict['saved values'],None]
    21522301        runtime = time.time()-begin
    21532302        chisq = np.sum(result[2]['fvec']**2)
Note: See TracChangeset for help on using the changeset viewer.