Changeset 1309 for trunk/GSASIIsasd.py


Ignore:
Timestamp:
Apr 30, 2014 2:05:57 PM (9 years ago)
Author:
vondreele
Message:

Implement structure factors for non-dilute small angle systems

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/GSASIIsasd.py

    r1302 r1309  
    435435    return []
    436436   
     437################################################################################
     438#### Structure factors for condensed systems
     439################################################################################
     440
     441def DiluteSF(Q,args=[]):
     442    ''' Default: no structure factor correction for dilute system
     443    '''
     444    return np.ones_like(Q)  #or return 1.0
     445
     446def HardSpheresSF(Q,args):
     447    ''' Computes structure factor for not dilute monodisperse hard spheres
     448    Refs.: PERCUS,YEVICK PHYS. REV. 110 1 (1958),THIELE J. CHEM PHYS. 39 474 (1968),
     449    WERTHEIM  PHYS. REV. LETT. 47 1462 (1981)
     450   
     451    param float Q: Q value array (A-1)
     452    param array args: [float R, float VolFrac]: interparticle distance & volume fraction
     453    returns numpy array S(Q)
     454    '''
     455   
     456    R,VolFr = args
     457    denom = (1.-VolFr)**4
     458    num = (1.+2.*VolFr)**2
     459    alp = num/denom
     460    bet = -6.*VolFr*(1.+VolFr/2.)**2/denom
     461    gamm = 0.5*VolFr*num/denom       
     462    A = 2.*Q*R
     463    A2 = A**2
     464    A3 = A**3
     465    A4 = A**4
     466    Rca = np.cos(A)
     467    Rsa = np.sin(A)
     468    calp = alp*(Rsa/A2-Rca/A)
     469    cbet = bet*(2.*Rsa/A2-(A2-2.)*Rca/A3-2./A3)
     470    cgam = gamm*(-Rca/A+(4./A)*((3.*A2-6.)*Rca/A4+(A2-6.)*Rsa/A3+6./A4))
     471    pfac = -24.*VolFr/A
     472    C = pfac*(calp+cbet+cgam)
     473    return 1./(1.-C)
     474       
     475def SquareWellSF(Q,args):
     476    ''' Computes structure factor for not dilute monodisperse hard sphere with a
     477    square well potential interaction.
     478    Refs.: SHARMA,SHARMA, PHYSICA 89A,(1977),212
     479   
     480    param float Q: Q value array (A-1)
     481    param array args: [float R, float VolFrac, float depth, float width]:
     482        interparticle distance, volume fraction (<0.08), well depth (e/kT<1.5kT),
     483        well width
     484    returns numpy array S(Q)
     485    well depth > 0 attractive & values outside above limits nonphysical cf.
     486    Monte Carlo simulations
     487    '''
     488    R,VolFr,Depth,Width = args
     489    eta = VolFr
     490    eta2 = eta*eta
     491    eta3 = eta*eta2
     492    eta4 = eta*eta3       
     493    etam1 = 1. - eta
     494    etam14 = etam1**4
     495    alp = (  ( (1. + 2.*eta)**2 ) + eta3*( eta-4.0 )  )/etam14
     496    bet = -(eta/3.0) * ( 18. + 20.*eta - 12.*eta2 + eta4 )/etam14
     497    gam = 0.5*eta*( (1. + 2.*eta)**2 + eta3*(eta-4.) )/etam14
     498
     499    SK = 2.*Q*R
     500    SK2 = SK*SK
     501    SK3 = SK*SK2
     502    SK4 = SK3*SK
     503    T1 = alp * SK3 * ( np.sin(SK) - SK * np.cos(SK) )
     504    T2 = bet * SK2 * ( 2.*SK*np.sin(SK) - (SK2-2.)*np.cos(SK) - 2.0 )
     505    T3 =   ( 4.0*SK3 - 24.*SK ) * np.sin(SK) 
     506    T3 = T3 - ( SK4 - 12.0*SK2 + 24.0 )*np.cos(SK) + 24.0   
     507    T3 = gam*T3
     508    T4 = -Depth*SK3*(np.sin(Width*SK) - Width*SK*np.cos(Width*SK)+ SK*np.cos(SK) - np.sin(SK) )
     509    CK =  -24.0*eta*( T1 + T2 + T3 + T4 )/SK3/SK3
     510    return 1./(1.-CK)
     511   
     512def StickyHardSpheresSF(Q,args):
     513    ''' Computes structure factor for not dilute monodisperse hard spheres
     514    Refs.: PERCUS,YEVICK PHYS. REV. 110 1 (1958),THIELE J. CHEM PHYS. 39 474 (1968),
     515    WERTHEIM  PHYS. REV. LETT. 47 1462 (1981)
     516   
     517    param float Q: Q value array (A-1)
     518    param array args: [float R, float VolFrac]: sphere radius & volume fraction
     519    returns numpy array S(Q)
     520    '''
     521    R,VolFr,epis,tau = args
     522#//      Input (fitting) variables are:
     523#       //radius = w[0]
     524#       //volume fraction = w[1]
     525#       //epsilon (perurbation param) = w[2]
     526#       //tau (stickiness) = w[3]
     527#       Variable rad,phi,eps,tau,eta
     528#       Variable sig,aa,etam1,qa,qb,qc,radic
     529#       Variable lam,lam2,test,mu,alpha,BetaVar
     530#       Variable qv,kk,k2,k3,ds,dc,aq1,aq2,aq3,aq,bq1,bq2,bq3,bq,sq
     531#       rad = w[0]
     532#       phi = w[1]
     533#       eps = w[2]
     534#       tau = w[3]
     535#       
     536#       eta = phi/(1.0-eps)/(1.0-eps)/(1.0-eps)
     537#       
     538#       sig = 2.0 * rad
     539#       aa = sig/(1.0 - eps)
     540#       etam1 = 1.0 - eta
     541#//C
     542#//C  SOLVE QUADRATIC FOR LAMBDA
     543#//C
     544#       qa = eta/12.0
     545#       qb = -1.0*(tau + eta/(etam1))
     546#       qc = (1.0 + eta/2.0)/(etam1*etam1)
     547#       radic = qb*qb - 4.0*qa*qc
     548#       if(radic<0)
     549#               if(x>0.01 && x<0.015)
     550#                       Print "Lambda unphysical - both roots imaginary"
     551#               endif
     552#               return(-1)
     553#       endif
     554#//C   KEEP THE SMALLER ROOT, THE LARGER ONE IS UNPHYSICAL
     555#       lam = (-1.0*qb-sqrt(radic))/(2.0*qa)
     556#       lam2 = (-1.0*qb+sqrt(radic))/(2.0*qa)
     557#       if(lam2<lam)
     558#               lam = lam2
     559#       endif
     560#       test = 1.0 + 2.0*eta
     561#       mu = lam*eta*etam1
     562#       if(mu>test)
     563#               if(x>0.01 && x<0.015)
     564#                Print "Lambda unphysical mu>test"
     565#               endif
     566#               return(-1)
     567#       endif
     568#       alpha = (1.0 + 2.0*eta - mu)/(etam1*etam1)
     569#       BetaVar = (mu - 3.0*eta)/(2.0*etam1*etam1)
     570#       
     571#//C
     572#//C   CALCULATE THE STRUCTURE FACTOR
     573#//C
     574#
     575#       qv = x
     576#       kk = qv*aa
     577#       k2 = kk*kk
     578#       k3 = kk*k2
     579#       ds = sin(kk)
     580#       dc = cos(kk)
     581#       aq1 = ((ds - kk*dc)*alpha)/k3
     582#       aq2 = (BetaVar*(1.0-dc))/k2
     583#       aq3 = (lam*ds)/(12.0*kk)
     584#       aq = 1.0 + 12.0*eta*(aq1+aq2-aq3)
     585#//
     586#       bq1 = alpha*(0.5/kk - ds/k2 + (1.0 - dc)/k3)
     587#       bq2 = BetaVar*(1.0/kk - ds/k2)
     588#       bq3 = (lam/12.0)*((1.0 - dc)/kk)
     589#       bq = 12.0*eta*(bq1+bq2-bq3)
     590#//
     591#       sq = 1.0/(aq*aq +bq*bq)
     592#
     593#       Return (sq)
     594    pass
     595   
     596#def HayterPenfoldSF(Q,args): #big & ugly function - do later (if ever)
     597#    pass
     598   
     599def InterPrecipitateSF(Q,args):
     600    ''' Computes structure factor for precipitates in a matrix
     601    Refs.: E-Wen Huang, Peter K. Liaw, Lionel Porcar, Yun Liu, Yee-Lang Liu,
     602    Ji-Jung Kai, and Wei-Ren Chen,APPLIED PHYSICS LETTERS 93, 161904 (2008)
     603    R. Giordano, A. Grasso, and J. Teixeira, Phys. Rev. A 43, 6894    1991   
     604    param float Q: Q value array (A-1)
     605    param array args: [float R, float VolFr]: "radius" & volume fraction
     606    returns numpy array S(Q)
     607    '''
     608    R,VolFr = args
     609    QV2 = Q**2*VolFr**2
     610    top = 1 - np.exp(-QV2/4)*np.cos(2.*Q*R)
     611    bot = 1-2*np.exp(-QV2/4)*np.cos(2.*Q*R) + np.exp(-QV2/2)
     612    return 2*(top/bot) - 1
    437613         
    438614################################################################################
    … …  
    9211097        'Unified disk':[UniDiskFF,UniDiskVol],'Sphere':[SphereFF,SphereVol],
    9221098        'Unified tube':[UniTubeFF,UniTubeVol],'Cylinder diam':[CylinderDFF,CylinderDVol]}
     1099           
     1100    sfxns = {'Dilute':DiluteSF,'Hard sphere':HardSpheresSF,'Square well':SquareWellSF,
     1101            'Sticky hard sphere':StickyHardSpheresSF,'InterPrecipitate':InterPrecipitateSF,}
     1102               
     1103    parmOrder = ['Volume','Radius','Mean','StdDev','MinSize','G','Rg','B','P','Cutoff',
     1104        'PkInt','PkPos','PkSig','PkGam',]
     1105       
     1106    FFparmOrder = ['Aspect ratio','Length','Diameter','Thickness']
     1107   
     1108    SFparmOrder = ['Dist','VolFr','epis','Sticky','Depth','Width']
    9231109
    9241110    def GetModelParms():
    925         parmOrder = ['Volume','Radius','Mean','StdDev','MinSize','G','Rg','B','P','Cutoff',
    926             'PkInt','PkPos','PkSig','PkGam',]
    9271111        parmDict = {'Scale':Sample['Scale'][0]}
    9281112        varyList = []
    … …  
    9471131                parmDict[cid+'FormFact'] = shapes[controls['FormFact']][0]
    9481132                parmDict[cid+'FFVolume'] = shapes[controls['FormFact']][1]
     1133                parmDict[cid+'StrFact'] = sfxns[controls['StrFact']]
    9491134                parmDict[cid+'XAnom density'] = Substances['Substances'][controls['Material']].get('XAnom density',0.0)
    950                 for item in ['Aspect ratio','Length','Thickness','Diameter',]:
     1135                for item in FFparmOrder:
    9511136                    if item in controls['FFargs']:
    9521137                        parmDict[cid+item] = controls['FFargs'][item][0]
    … …  
    9541139                            varyList.append(cid+item)
    9551140                            values.append(controls['FFargs'][item][0])
     1141                for item in SFparmOrder:
     1142                    if item in controls.get('SFargs',{}):
     1143                        parmDict[cid+item] = controls['SFargs'][item][0]
     1144                        if controls['SFargs'][item][1]:
     1145                            varyList.append(cid+item)
     1146                            values.append(controls['SFargs'][item][0])
    9561147            distDict = controls['DistType']
    9571148            for item in parmOrder:
    … …  
    9701161        partData = Model['Particle']
    9711162        for i,level in enumerate(partData['Levels']):
    972             Type = level['Controls']['DistType']
    973             print ' Component %d: Type: %s:'%(i,Type)
    974             cid = str(i)+':'
    9751163            controls = level['Controls']
    9761164            Type = controls['DistType']
    9771165            if Type in ['LogNormal','Gaussian','LSW','Schulz-Zimm','Monodisperse']:
    978                 for item in ['Aspect ratio','Length','Thickness','Diameter',]:
     1166                print ' Component %d: Type: %s: Structure Factor: %s'%(i,Type,controls['StrFact'])               
     1167            else:
     1168                print ' Component %d: Type: %s: '%(i,Type,)
     1169            cid = str(i)+':'
     1170            if Type in ['LogNormal','Gaussian','LSW','Schulz-Zimm','Monodisperse']:
     1171                for item in FFparmOrder:
    9791172                    if cid+item in varyList:
    9801173                        controls['FFargs'][item][0] = parmDict[cid+item]
     1174                        print ' %15s: %15.4g esd: %15.4g'%(cid+item,parmDict[cid+item],sigDict[cid+item])
     1175                for item in SFparmOrder:
     1176                    if cid+item in varyList:
     1177                        controls['SFargs'][item][0] = parmDict[cid+item]
    9811178                        print ' %15s: %15.4g esd: %15.4g'%(cid+item,parmDict[cid+item],sigDict[cid+item])
    9821179            distDict = controls['DistType']
    … …  
    9991196                FFfxn = parmDict[cid+'FormFact']
    10001197                Volfxn = parmDict[cid+'FFVolume']
     1198                SFfxn = parmDict[cid+'StrFact']
    10011199                FFargs = []
    1002                 for item in [cid+'Aspect ratio',cid+'Length',cid+'Thickness',cid+'Diameter',]:
     1200                SFargs = []
     1201                for item in [cid+'Aspect ratio',cid+'Length',cid+'Thickness',cid+'Diameter']:
    10031202                    if item in parmDict:
    10041203                        FFargs.append(parmDict[item])
     1204                for item in [cid+'Dist',cid+'VolFr',cid+'epis',cid+'Sticky',cid+'Depth',cid+'Width']:
     1205                    if item in parmDict:
     1206                        SFargs.append(parmDict[item])
    10051207                distDict = {}
    10061208                for item in [cid+'Volume',cid+'Mean',cid+'StdDev',cid+'MinSize',]:
    … …  
    10121214                Gmat = G_matrix(Q,rBins,contrast,FFfxn,Volfxn,FFargs).T
    10131215                dist *= parmDict[cid+'Volume']
    1014                 Ic += np.dot(Gmat,dist)
     1216                Ic += np.dot(Gmat,dist)*SFfxn(Q,args=SFargs)
    10151217            elif 'Unified' in Type:
    10161218                Rg,G,B,P,Rgco = parmDict[cid+'Rg'],parmDict[cid+'G'],parmDict[cid+'B'], \
    … …  
    10271229                FFfxn = parmDict[cid+'FormFact']
    10281230                Volfxn = parmDict[cid+'FFVolume']
     1231                SFfxn = parmDict[cid+'StrFact']
    10291232                FFargs = []
     1233                SFargs = []
    10301234                for item in [cid+'Aspect ratio',cid+'Length',cid+'Thickness',cid+'Diameter',]:
    10311235                    if item in parmDict:
    … …  
    10351239                R = parmDict[cid+'Radius']
    10361240                Gmat = G_matrix(Q,R,contrast,FFfxn,Volfxn,FFargs)             
    1037                 Ic += Gmat[0]*parmDict[cid+'Volume']
     1241                Ic += Gmat[0]*parmDict[cid+'Volume']*SFfxn(Q,args=SFargs)
    10381242            elif 'Bragg' in distFxn:
    10391243                Ic += parmDict[cid+'PkInt']*G2pwd.getPsVoigt(parmDict[cid+'PkPos'],
    … …  
    10841288   
    10851289def ModelFxn(Profile,ProfDict,Limits,Substances,Sample,sasdData):
     1290   
    10861291    shapes = {'Spheroid':[SpheroidFF,SpheroidVol],'Cylinder':[CylinderDFF,CylinderDVol],
    10871292        'Cylinder AR':[CylinderARFF,CylinderARVol],'Unified sphere':[UniSphereFF,UniSphereVol],
    … …  
    10891294        'Unified disk':[UniDiskFF,UniDiskVol],'Sphere':[SphereFF,SphereVol],
    10901295        'Unified tube':[UniTubeFF,UniTubeVol],'Cylinder diam':[CylinderDFF,CylinderDVol]}
     1296    sfxns = {'Dilute':DiluteSF,'Hard sphere':HardSpheresSF,'Square well':SquareWellSF,
     1297            'Sticky hard sphere':StickyHardSpheresSF,'InterPrecipitate':InterPrecipitateSF,}
     1298
    10911299#    pdb.set_trace()
    10921300    partData = sasdData['Particle']
    … …  
    11091317        distFxn = controls['DistType']
    11101318        if distFxn in ['LogNormal','Gaussian','LSW','Schulz-Zimm']:
     1319            parmDict = level[controls['DistType']]
    11111320            FFfxn = shapes[controls['FormFact']][0]
    11121321            Volfxn = shapes[controls['FormFact']][1]
     1322            SFfxn = sfxns[controls['StrFact']]
    11131323            FFargs = []
     1324            SFargs = []
     1325            for item in ['Dist','VolFr','epis','Sticky','Depth','Width',]:
     1326                if item in controls.get('SFargs',{}):
     1327                    SFargs.append(controls['SFargs'][item][0])
    11141328            for item in ['Aspect ratio','Length','Thickness','Diameter',]:
    11151329                if item in controls['FFargs']:
    … …  
    11171331            rho = Substances['Substances'][level['Controls']['Material']].get('XAnom density',0.0)
    11181332            contrast = rho**2-rhoMat**2
    1119             parmDict = level[controls['DistType']]
    11201333            distDict = {}
    11211334            for item in parmDict:
    … …  
    11241337            Gmat = G_matrix(Q[Ibeg:Ifin],rBins,contrast,FFfxn,Volfxn,FFargs).T
    11251338            dist *= level[distFxn]['Volume'][0]
    1126             Ic[Ibeg:Ifin] += np.dot(Gmat,dist)
     1339            Ic[Ibeg:Ifin] += np.dot(Gmat,dist)*SFfxn(Q[Ibeg:Ifin],args=SFargs)
    11271340            Rbins.append(rBins)
    11281341            Dist.append(dist/(4.*dBins))
    … …  
    11451358            Dist.append([])
    11461359        elif 'Mono' in distFxn:
     1360            parmDict = level[controls['DistType']]
     1361            R = level[controls['DistType']]['Radius'][0]
    11471362            FFfxn = shapes[controls['FormFact']][0]
    11481363            Volfxn = shapes[controls['FormFact']][1]
     1364            SFfxn = sfxns[controls['StrFact']]
    11491365            FFargs = []
     1366            SFargs = [parmDict['Volume'][0],R]
     1367            for item in ['epis','Sticky','Depth','Width',]:
     1368                if item in controls.get('SFargs',{}):
     1369                    SFargs.append(controls['SFargs'][item][0])
    11501370            for item in ['Aspect ratio','Length','Thickness','Diameter',]:
    11511371                if item in controls['FFargs']:
    … …  
    11531373            rho = Substances['Substances'][level['Controls']['Material']].get('XAnom density',0.0)
    11541374            contrast = rho**2-rhoMat**2
    1155             R = level[controls['DistType']]['Radius'][0]
    11561375            Gmat = G_matrix(Q[Ibeg:Ifin],R,contrast,FFfxn,Volfxn,FFargs)             
    1157             Ic[Ibeg:Ifin] += Gmat[0]*level[distFxn]['Volume'][0]
     1376            Ic[Ibeg:Ifin] += Gmat[0]*level[distFxn]['Volume'][0]*SFfxn(Q[Ibeg:Ifin],args=SFargs)
    11581377            Rbins.append([])
    11591378            Dist.append([])
Note: See TracChangeset for help on using the changeset viewer.