Changeset 1309 for trunk/GSASIIsasd.py
 Timestamp:
 Apr 30, 2014 2:05:57 PM (9 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/GSASIIsasd.py
r1302 r1309 435 435 return [] 436 436 437 ################################################################################ 438 #### Structure factors for condensed systems 439 ################################################################################ 440 441 def DiluteSF(Q,args=[]): 442 ''' Default: no structure factor correction for dilute system 443 ''' 444 return np.ones_like(Q) #or return 1.0 445 446 def 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 (A1) 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/A2Rca/A) 469 cbet = bet*(2.*Rsa/A2(A22.)*Rca/A32./A3) 470 cgam = gamm*(Rca/A+(4./A)*((3.*A26.)*Rca/A4+(A26.)*Rsa/A3+6./A4)) 471 pfac = 24.*VolFr/A 472 C = pfac*(calp+cbet+cgam) 473 return 1./(1.C) 474 475 def 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 (A1) 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*( eta4.0 ) )/etam14 496 bet = (eta/3.0) * ( 18. + 20.*eta  12.*eta2 + eta4 )/etam14 497 gam = 0.5*eta*( (1. + 2.*eta)**2 + eta3*(eta4.) )/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)  (SK22.)*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 512 def 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 (A1) 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.0eps)/(1.0eps)/(1.0eps) 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*qbsqrt(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.0dc))/k2 583 # aq3 = (lam*ds)/(12.0*kk) 584 # aq = 1.0 + 12.0*eta*(aq1+aq2aq3) 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+bq2bq3) 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 599 def InterPrecipitateSF(Q,args): 600 ''' Computes structure factor for precipitates in a matrix 601 Refs.: EWen Huang, Peter K. Liaw, Lionel Porcar, Yun Liu, YeeLang Liu, 602 JiJung Kai, and WeiRen 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 (A1) 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 = 12*np.exp(QV2/4)*np.cos(2.*Q*R) + np.exp(QV2/2) 612 return 2*(top/bot)  1 437 613 438 614 ################################################################################ … … 921 1097 'Unified disk':[UniDiskFF,UniDiskVol],'Sphere':[SphereFF,SphereVol], 922 1098 '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'] 923 1109 924 1110 def GetModelParms(): 925 parmOrder = ['Volume','Radius','Mean','StdDev','MinSize','G','Rg','B','P','Cutoff',926 'PkInt','PkPos','PkSig','PkGam',]927 1111 parmDict = {'Scale':Sample['Scale'][0]} 928 1112 varyList = [] … … 947 1131 parmDict[cid+'FormFact'] = shapes[controls['FormFact']][0] 948 1132 parmDict[cid+'FFVolume'] = shapes[controls['FormFact']][1] 1133 parmDict[cid+'StrFact'] = sfxns[controls['StrFact']] 949 1134 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: 951 1136 if item in controls['FFargs']: 952 1137 parmDict[cid+item] = controls['FFargs'][item][0] … … 954 1139 varyList.append(cid+item) 955 1140 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]) 956 1147 distDict = controls['DistType'] 957 1148 for item in parmOrder: … … 970 1161 partData = Model['Particle'] 971 1162 for i,level in enumerate(partData['Levels']): 972 Type = level['Controls']['DistType']973 print ' Component %d: Type: %s:'%(i,Type)974 cid = str(i)+':'975 1163 controls = level['Controls'] 976 1164 Type = controls['DistType'] 977 1165 if Type in ['LogNormal','Gaussian','LSW','SchulzZimm','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','SchulzZimm','Monodisperse']: 1171 for item in FFparmOrder: 979 1172 if cid+item in varyList: 980 1173 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] 981 1178 print ' %15s: %15.4g esd: %15.4g'%(cid+item,parmDict[cid+item],sigDict[cid+item]) 982 1179 distDict = controls['DistType'] … … 999 1196 FFfxn = parmDict[cid+'FormFact'] 1000 1197 Volfxn = parmDict[cid+'FFVolume'] 1198 SFfxn = parmDict[cid+'StrFact'] 1001 1199 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']: 1003 1202 if item in parmDict: 1004 1203 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]) 1005 1207 distDict = {} 1006 1208 for item in [cid+'Volume',cid+'Mean',cid+'StdDev',cid+'MinSize',]: … … 1012 1214 Gmat = G_matrix(Q,rBins,contrast,FFfxn,Volfxn,FFargs).T 1013 1215 dist *= parmDict[cid+'Volume'] 1014 Ic += np.dot(Gmat,dist) 1216 Ic += np.dot(Gmat,dist)*SFfxn(Q,args=SFargs) 1015 1217 elif 'Unified' in Type: 1016 1218 Rg,G,B,P,Rgco = parmDict[cid+'Rg'],parmDict[cid+'G'],parmDict[cid+'B'], \ … … 1027 1229 FFfxn = parmDict[cid+'FormFact'] 1028 1230 Volfxn = parmDict[cid+'FFVolume'] 1231 SFfxn = parmDict[cid+'StrFact'] 1029 1232 FFargs = [] 1233 SFargs = [] 1030 1234 for item in [cid+'Aspect ratio',cid+'Length',cid+'Thickness',cid+'Diameter',]: 1031 1235 if item in parmDict: … … 1035 1239 R = parmDict[cid+'Radius'] 1036 1240 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) 1038 1242 elif 'Bragg' in distFxn: 1039 1243 Ic += parmDict[cid+'PkInt']*G2pwd.getPsVoigt(parmDict[cid+'PkPos'], … … 1084 1288 1085 1289 def ModelFxn(Profile,ProfDict,Limits,Substances,Sample,sasdData): 1290 1086 1291 shapes = {'Spheroid':[SpheroidFF,SpheroidVol],'Cylinder':[CylinderDFF,CylinderDVol], 1087 1292 'Cylinder AR':[CylinderARFF,CylinderARVol],'Unified sphere':[UniSphereFF,UniSphereVol], … … 1089 1294 'Unified disk':[UniDiskFF,UniDiskVol],'Sphere':[SphereFF,SphereVol], 1090 1295 '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 1091 1299 # pdb.set_trace() 1092 1300 partData = sasdData['Particle'] … … 1109 1317 distFxn = controls['DistType'] 1110 1318 if distFxn in ['LogNormal','Gaussian','LSW','SchulzZimm']: 1319 parmDict = level[controls['DistType']] 1111 1320 FFfxn = shapes[controls['FormFact']][0] 1112 1321 Volfxn = shapes[controls['FormFact']][1] 1322 SFfxn = sfxns[controls['StrFact']] 1113 1323 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]) 1114 1328 for item in ['Aspect ratio','Length','Thickness','Diameter',]: 1115 1329 if item in controls['FFargs']: … … 1117 1331 rho = Substances['Substances'][level['Controls']['Material']].get('XAnom density',0.0) 1118 1332 contrast = rho**2rhoMat**2 1119 parmDict = level[controls['DistType']]1120 1333 distDict = {} 1121 1334 for item in parmDict: … … 1124 1337 Gmat = G_matrix(Q[Ibeg:Ifin],rBins,contrast,FFfxn,Volfxn,FFargs).T 1125 1338 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) 1127 1340 Rbins.append(rBins) 1128 1341 Dist.append(dist/(4.*dBins)) … … 1145 1358 Dist.append([]) 1146 1359 elif 'Mono' in distFxn: 1360 parmDict = level[controls['DistType']] 1361 R = level[controls['DistType']]['Radius'][0] 1147 1362 FFfxn = shapes[controls['FormFact']][0] 1148 1363 Volfxn = shapes[controls['FormFact']][1] 1364 SFfxn = sfxns[controls['StrFact']] 1149 1365 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]) 1150 1370 for item in ['Aspect ratio','Length','Thickness','Diameter',]: 1151 1371 if item in controls['FFargs']: … … 1153 1373 rho = Substances['Substances'][level['Controls']['Material']].get('XAnom density',0.0) 1154 1374 contrast = rho**2rhoMat**2 1155 R = level[controls['DistType']]['Radius'][0]1156 1375 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) 1158 1377 Rbins.append([]) 1159 1378 Dist.append([])
Note: See TracChangeset
for help on using the changeset viewer.