Changeset 1775


Ignore:
Timestamp:
Apr 2, 2015 3:43:31 PM (7 years ago)
Author:
vondreele
Message:

fix to SH texture stuff

Location:
trunk
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIIstrIO.py

    r1773 r1775  
    20712071                    controlDict[pfx+'SHord'] = hapData['Pref.Ori.'][4]
    20722072                    controlDict[pfx+'SHncof'] = len(hapData['Pref.Ori.'][5])
    2073                     controlDict[pfx+'SHhkl'] = [eval(a.replace(' ',',')) for a in hapData['Pref.Ori.'][6]]
     2073                    controlDict[pfx+'SHhkl'] = []
     2074                    if hapData['Pref.Ori.'][6][0] != '':
     2075                        controlDict[pfx+'SHhkl'] = [eval(a.replace(' ',',')) for a in hapData['Pref.Ori.'][6]]
    20742076                    controlDict[pfx+'SHtoler'] = hapData['Pref.Ori.'][7]
    20752077                    for item in hapData['Pref.Ori.'][5]:
  • trunk/GSASIIstrMath.py

    r1770 r1775  
    304304################################################################################
    305305
    306 def penaltyFxn(HistoPhases,parmDict,varyList):
     306def penaltyFxn(HistoPhases,calcControls,parmDict,varyList):
    307307    'Needs a doc string'
    308308    Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
     
    327327        names = [['Bond','Bonds'],['Angle','Angles'],['Plane','Planes'],
    328328            ['Chiral','Volumes'],['Torsion','Torsions'],['Rama','Ramas'],
    329             ['ChemComp','Sites'],['Texture','HKLs']]
     329            ['ChemComp','Sites'],['Texture','HKLs'],]
    330330        for name,rest in names:
    331331            pWsum[name] = 0.
     
    389389                            pVals.append(Z1[ind[0]][ind[1]])
    390390                            pWt.append(wt/esd1**2)
    391                             pWsum[name] += wt*((obs-calc)/esd)**2
     391                            pWsum[name] += wt*(-Z1[ind[0]][ind[1]]/esd)**2
    392392                        if ifesd2:
    393393                            Z2 = 1.-Z
     
    397397                                pWt.append(wt/esd2**2)
    398398                                pWsum[name] += wt*((obs-calc)/esd)**2
    399          
     399       
     400        name = 'SH-Pref.Ori.'
     401        pWsum[name] = 0.0
     402        for hist in Phases[phase]['Histograms']:
     403            if hist in Histograms:
     404                hId = Histograms[hist]['hId']
     405                phfx = '%d:%d:'%(pId,hId)
     406                if calcControls[phfx+'poType'] == 'SH':
     407                    toler = calcControls[phfx+'SHtoler']
     408                    wt = 1./toler**2
     409                    HKLs = calcControls[phfx+'SHhkl']
     410                    SHcofNames = Phases[phase]['Histograms'][hist]['Pref.Ori.'][5].keys()
     411                    SHcof = dict(zip(SHcofNames,[parmDict[phfx+cof] for cof in SHcofNames]))
     412                    for i,PH in enumerate(HKLs):
     413                        phi,beta = G2lat.CrsAng(PH,cell,SGData)
     414                        SH3Coef = {}
     415                        for item in SHcof:
     416                            L,N = eval(item.strip('C'))
     417                            SH3Coef['C%d,0,%d'%(L,N)] = SHcof[item]                       
     418                        ODFln = G2lat.Flnh(False,SH3Coef,phi,beta,SGData)
     419                        X = np.linspace(0,90.0,26)
     420                        Y = -ma.masked_greater(G2lat.polfcal(ODFln,'0',X,0.0),0.0)
     421                        IndY = ma.nonzero(Y)
     422                        for ind in IndY[0]:
     423                            pNames.append('%d:%d:%s:%d:%.2f'%(pId,hId,name,i,X[ind]))
     424                            pVals.append(Y[ind])
     425                            pWt.append(wt)
     426                            pWsum[name] += wt*(-Y[ind])**2
    400427    pWsum['PWLref'] = 0.
    401428    for item in varyList:
     
    411438    return pNames,pVals,pWt,pWsum
    412439   
    413 def penaltyDeriv(pNames,pVal,HistoPhases,parmDict,varyList):
     440def penaltyDeriv(pNames,pVal,HistoPhases,calcControls,parmDict,varyList):
    414441    'Needs a doc string'
    415442    Histograms,Phases,restraintDict,rigidbodyDict = HistoPhases
     
    443470                if 'PWL' in pName:
    444471                    pDerv[varyList.index(pName)][ip] += 1.
     472                    continue
     473                elif 'SH-' in pName:
    445474                    continue
    446475                id = int(pnames[2])
     
    504533                    except ValueError:
    505534                        pass
     535       
     536#        for ip,pName in enumerate(pNames):
     537#            pnames = pNames.split(':')
     538#            if 'SH-' in pName and pId == int(pnames[0]):
     539#                name = pnames[2]
     540#                L,N = eval(name.strip('C'))
     541#                hId = int(pnames[1])
     542#                phfx = '%d:%d:'%(pId,hId)
     543#                hklId = int(pnames[3])
     544#                gam = float(pnames[4])
     545#                HKLs = calcControls[phfx+'SHhkl']
     546#                phi,beta = G2lat.CrsAng(HKLs[hklId],cell,SGData)
     547#                SHcofNames = Phases[phase]['Histograms'][hist]['Pref.Ori.'][5].keys()
     548#                SHcof = dict(zip(SHcofNames,[parmDict[phfx+cof] for cof in SHcofNames]))
     549#           
     550#            raise Exception
     551
    506552    return pDerv
    507553
     
    19391985                refl[5+im] = GetReflPos(refl,im,wave,A,pfx,hfx,calcControls,parmDict)         #corrected reflection position
    19401986                Lorenz = 1./(2.*sind(refl[5+im]/2.)**2*cosd(refl[5+im]/2.))           #Lorentz correction
    1941 #                refl[5+im] += GetHStrainShift(refl,im,SGData,phfx,hfx,calcControls,parmDict)               #apply hydrostatic strain shift
    19421987                refl[6+im:8+im] = GetReflSigGamCW(refl,im,wave,G,GB,phfx,calcControls,parmDict)    #peak sig & gam
    19431988                refl[11+im:15+im] = GetIntensityCorr(refl,im,Uniq,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict)
     
    24592504            dMdv = np.sqrt(wtFactor)*dMdvh
    24602505           
    2461     pNames,pVals,pWt,pWsum = penaltyFxn(HistoPhases,parmDict,varylist)
     2506    pNames,pVals,pWt,pWsum = penaltyFxn(HistoPhases,calcControls,parmDict,varylist)
    24622507    if np.any(pVals):
    2463         dpdv = penaltyDeriv(pNames,pVals,HistoPhases,parmDict,varylist)
     2508        dpdv = penaltyDeriv(pNames,pVals,HistoPhases,calcControls,parmDict,varylist)
    24642509        dMdv = np.concatenate((dMdv.T,(np.sqrt(pWt)*dpdv).T)).T
    24652510       
     
    25342579        else:
    25352580            continue        #skip non-histogram entries
    2536     pNames,pVals,pWt,pWsum = penaltyFxn(HistoPhases,parmDict,varylist)
     2581    pNames,pVals,pWt,pWsum = penaltyFxn(HistoPhases,calcControls,parmDict,varylist)
    25372582    if np.any(pVals):
    2538         dpdv = penaltyDeriv(pNames,pVals,HistoPhases,parmDict,varylist)
     2583        dpdv = penaltyDeriv(pNames,pVals,HistoPhases,calcControls,parmDict,varylist)
    25392584        Vec += np.sum(dpdv*pWt*pVals,axis=1)
    25402585        Hess += np.inner(dpdv*pWt,dpdv)
     
    26862731            dlg.Destroy()
    26872732            raise Exception         #Abort!!
    2688     pDict,pVals,pWt,pWsum = penaltyFxn(HistoPhases,parmDict,varylist)
     2733    pDict,pVals,pWt,pWsum = penaltyFxn(HistoPhases,calcControls,parmDict,varylist)
    26892734    if len(pVals):
    26902735        pSum = np.sum(pWt*pVals**2)
    26912736        for name in pWsum:
    2692             if pWsum:
     2737            if pWsum[name]:
    26932738                print '  Penalty function for %8s = %12.5g'%(name,pWsum[name])
    26942739        print 'Total penalty function: %12.5g on %d terms'%(pSum,len(pVals))
Note: See TracChangeset for help on using the changeset viewer.