Changeset 520
 Timestamp:
 Mar 14, 2012 3:17:17 PM (12 years ago)
 Location:
 trunk
 Files:

 3 edited
Legend:
 Unmodified
 Added
 Removed

trunk/GSASIIlattice.py
r432 r520 590 590 HKL.append([h,k,l,rdsq2d(rdsq,6),1]) 591 591 return sortHKLd(HKL,True,False) 592 593 def getHKLmax(dmin,SGData,A): 594 #finds maximum allowed hkl for given A within dmin 595 SGLaue = SGData['SGLaue'] 596 if SGLaue in ['3R','3mR']: #Rhombohedral axes 597 Hmax = [0,0,0] 598 cell = A2cell(A) 599 aHx = cell[0]*math.sqrt(2.0*(1.0cosd(cell[3]))) 600 cHx = cell[0]*math.sqrt(3.0*(1.0+2.0*cosd(cell[3]))) 601 Hmax[0] = Hmax[1] = int(round(aHx/dmin)) 602 Hmax[2] = int(round(cHx/dmin)) 603 #print Hmax,aHx,cHx 604 else: # all others 605 Hmax = MaxIndex(dmin,A) 606 return Hmax 592 607 593 608 def GenHLaue(dmin,SGData,A): … … 614 629 SGUniq = SGData['SGUniq'] 615 630 #finds maximum allowed hkl for given A within dmin 616 if SGLaue in ['3R','3mR']: #Rhombohedral axes 617 Hmax = [0,0,0] 618 cell = A2cell(A) 619 aHx = cell[0]*math.sqrt(2.0*(1.0cosd(cell[3]))) 620 cHx = cell[0]*math.sqrt(3.0*(1.0+2.0*cosd(cell[3]))) 621 Hmax[0] = Hmax[1] = int(round(aHx/dmin)) 622 Hmax[2] = int(round(cHx/dmin)) 623 #print Hmax,aHx,cHx 624 else: # all others 625 Hmax = MaxIndex(dmin,A) 631 Hmax = getHKLmax(dmin,SGData,A) 626 632 627 633 dminsq = 1./(dmin**2) 
trunk/GSASIIphsGUI.py
r519 r520 252 252 generalData['POhkl'] = [0,0,1] 253 253 if 'Map' not in generalData: 254 generalData['Map'] = {'MapType':'','RefList':'','Resolution': 4.0}254 generalData['Map'] = {'MapType':'','RefList':'','Resolution':1.0,'rhoMax':100.,'rho':[],'rhoMax':0.} 255 255 # if 'SH Texture' not in generalData: 256 256 # generalData['SH Texture'] = data['SH Texture'] … … 274 274 Info = G2elem.GetAtomInfo(atom[ct]) 275 275 generalData['AtomTypes'].append(atom[ct]) 276 generalData['Z'] = Info['Z'] 276 277 generalData['Isotopes'][atom[ct]] = Info['Isotopes'] 277 278 generalData['BondRadii'].append(Info['Drad']) … … 285 286 generalData['NoAtoms'][atom[ct]] = atom[cs1]*float(atom[cs+1]) 286 287 generalData['Color'].append(Info['Color']) 288 F000X = 0. 289 F000N = 0. 290 for i,elem in enumerate(generalData['AtomTypes']): 291 F000X += generalData['NoAtoms'][elem]*generalData['Z'] 292 isotope = generalData['Isotope'][elem] 293 F000N += generalData['NoAtoms'][elem]*generalData['Isotopes'][elem][isotope][1] 294 generalData['F000X'] = F000X 295 generalData['F000N'] = F000N 296 287 297 288 298 ################################################################################ … … 308 318 SetupGeneral() 309 319 generalData = data['General'] 310 Map = generalData['Map'] # {'MapType':'','RefList':'','Resolution': 4.0}320 Map = generalData['Map'] # {'MapType':'','RefList':'','Resolution':1.0,'rho':[],'rhoMax':0.} 311 321 312 322 def NameSizer(): … … 614 624 except ValueError: 615 625 pass 616 mapRes.SetValue("%. 1f"%(Map['Resolution'])) #reset in case of error626 mapRes.SetValue("%.2f"%(Map['Resolution'])) #reset in case of error 617 627 618 628 mapTypes = ['Fobs','Fcalc','deltF','2*FoFc','Patterson'] … … 633 643 mapSizer.Add(refList,0,wx.ALIGN_CENTER_VERTICAL) 634 644 mapSizer.Add(wx.StaticText(dataDisplay,label=' Resolution: '),0,wx.ALIGN_CENTER_VERTICAL) 635 mapRes = wx.TextCtrl(dataDisplay,value='%. 1f'%(Map['Resolution']),style=wx.TE_PROCESS_ENTER)645 mapRes = wx.TextCtrl(dataDisplay,value='%.2f'%(Map['Resolution']),style=wx.TE_PROCESS_ENTER) 636 646 mapRes.Bind(wx.EVT_TEXT_ENTER,OnResVal) 637 647 mapRes.Bind(wx.EVT_KILL_FOCUS,OnResVal) … … 1190 1200 'L','K','M','F','P','S','T','W','Y','V','M',' ',' ',' '] 1191 1201 defaultDrawing = {'viewPoint':[[0.5,0.5,0.5],[]],'showHydrogen':True,'backColor':[0,0,0],'depthFog':False, 1192 'Zclip':50.0,'cameraPos':50.,'radiusFactor':0.85, 1202 'Zclip':50.0,'cameraPos':50.,'radiusFactor':0.85,'contourLevel':1., 1193 1203 'bondRadius':0.1,'ballScale':0.33,'vdwScale':0.67,'ellipseProb':50,'sizeH':0.50, 1194 1204 'unitCellBox':False,'showABC':True,'selectedAtoms':[], … … 1202 1212 drawingData = copy.copy(defaultDrawing) 1203 1213 drawingData['Atoms'] = [] 1214 if 'contourLevel' not in drawingData: 1215 drawingData['contourLevel'] = 1. 1204 1216 cx,ct,cs = [0,0,0] 1205 1217 if generalData['Type'] == 'nuclear': … … 2007 2019 bondRadiusTxt.SetLabel('Bond radius, A: '+'%.2f'%(drawingData['bondRadius'])) 2008 2020 G2plt.PlotStructure(G2frame,data) 2021 2022 def OnContourLevel(event): 2023 drawingData['contourLevel'] = contourLevel.GetValue()/100. 2024 contourLevelTxt.SetLabel('Contour level: '+'%.2f'%(drawingData['contourLevel']*generalData['Map']['rhoMax'])) 2025 G2plt.PlotStructure(G2frame,data) 2009 2026 2010 2027 slopSizer = wx.BoxSizer(wx.HORIZONTAL) 2011 slideSizer = wx.FlexGridSizer( 6,2)2028 slideSizer = wx.FlexGridSizer(7,2) 2012 2029 slideSizer.AddGrowableCol(1,1) 2013 2030 … … 2052 2069 bondRadius.Bind(wx.EVT_SLIDER, OnBondRadius) 2053 2070 slideSizer.Add(bondRadius,1,wx.EXPANDwx.RIGHT) 2071 2072 if generalData['Map']['rhoMax']: 2073 contourLevelTxt = wx.StaticText(dataDisplay,1,' Contour level: '+'%.2f'%(drawingData['contourLevel']*generalData['Map']['rhoMax'])) 2074 slideSizer.Add(contourLevelTxt,0,wx.ALIGN_CENTER_VERTICAL) 2075 contourLevel = wx.Slider(dataDisplay,style=wx.SL_HORIZONTAL,value=int(100*drawingData['contourLevel'])) 2076 contourLevel.SetRange(1,100) 2077 contourLevel.Bind(wx.EVT_SLIDER, OnContourLevel) 2078 slideSizer.Add(contourLevel,1,wx.EXPANDwx.RIGHT) 2054 2079 2055 2080 slopSizer.Add(slideSizer,1,wx.EXPANDwx.RIGHT) … … 3314 3339 3315 3340 def OnFourierMaps(event): 3341 import scipy.fftpack as fft 3316 3342 generalData = data['General'] 3343 # for item in generalData: 3344 # print item,generalData[item] 3317 3345 if not generalData['Map']['MapType']: 3318 3346 print '**** ERROR  Fourier map not defined' 3319 3347 return 3320 print 'Calculate Fourier maps' 3321 print generalData['Map'] 3348 mapData = generalData['Map'] 3349 dmin = mapData['Resolution'] 3350 phaseName = generalData['Name'] 3351 SGData = generalData['SGData'] 3352 cell = generalData['Cell'][1:8] 3353 A = G2lat.cell2A(cell[:6]) 3354 Hmax = np.asarray(G2lat.getHKLmax(dmin,SGData,A),dtype='i')+1 3355 Fhkl = np.zeros(shape=2*Hmax,dtype='c16') 3356 reflName = mapData['RefList'] 3357 if 'PWDR' in reflName: 3358 PatternId = G2gd.GetPatternTreeItemId(G2frame,G2frame.root, reflName) 3359 reflSets = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,PatternId,'Reflection Lists')) 3360 reflData = reflSets[phaseName] 3361 elif 'HKLF' in reflName: 3362 print 'single crystal reflections' 3363 return 3364 # Fhkl[0,0,0] = generalData['F000X'] 3365 time0 = time.time() 3366 for ref in reflData: 3367 if ref[4] >= dmin: 3368 for i,hkl in enumerate(ref[11]): 3369 hkl = np.asarray(hkl,dtype='i') 3370 Fosq,Fcsq,ph = ref[8:11] 3371 dp = 360.*ref[12][i] 3372 a = cosd(ph+dp) 3373 b = sind(ph+dp) 3374 phasep = complex(a,b) 3375 phasem = complex(a,b) 3376 if 'Fobs' in mapData['MapType']: 3377 F = np.sqrt(Fosq) 3378 h,k,l = hkl+Hmax 3379 Fhkl[h,k,l] = F*phasep 3380 h,k,l = hkl+Hmax 3381 Fhkl[h,k,l] = F*phasem 3382 elif 'Fcalc' in mapData['MapType']: 3383 F = np.sqrt(Fcsq) 3384 h,k,l = hkl+Hmax 3385 Fhkl[h,k,l] = F*phasep 3386 h,k,l = hkl+Hmax 3387 Fhkl[h,k,l] = F*phasem 3388 elif 'deltF' in mapData['MapType']: 3389 dF = np.sqrt(Fosq)np.sqrt(Fcsq) 3390 h,k,l = hkl+Hmax 3391 Fhkl[h,k,l] = dF*phasep 3392 h,k,l = hkl+Hmax 3393 Fhkl[h,k,l] = dF*phasem 3394 elif '2*FoFc' in mapData['MapType']: 3395 F = 2.*np.sqrt(Fosq)np.sqrt(Fcsq) 3396 h,k,l = hkl+Hmax 3397 Fhkl[h,k,l] = F*phasep 3398 h,k,l = hkl+Hmax 3399 Fhkl[h,k,l] = F*phasem 3400 elif 'Patterson' in mapData['MapType']: 3401 h,k,l = hkl+Hmax 3402 Fhkl[h,k,l] = complex(Fosq,0.) 3403 h,k,l = hkl+Hmax 3404 Fhkl[h,k,l] = complex(Fosq,0.) 3405 Fhkl = fft.fftshift(Fhkl) 3406 rho = fft.fftn(Fhkl)/cell[6] 3407 print 'Fourier map time: %.4f'%(time.time()time0),'no. elements: %d'%(Fhkl.size) 3408 mapData['rho'] = np.real(rho) 3409 mapData['rhoMax'] = np.max(np.real(rho)) 3410 data['Drawing']['contourLevel'] = 1. 3411 print mapData['MapType']+' computed: rhomax = %.3f rhomin = %.3f'%(mapData['rhoMax'],np.min(np.real(rho))) 3412 ## map printing for testing purposes 3413 # ix,jy,kz = mapData['rho'].shape 3414 # for k in range(kz): 3415 # print 'k = ',k 3416 # for j in range(jy): 3417 # line = '' 3418 # if SGData['SGLaue'] in ['3','3m1','31m','6/m','6/mmm']: 3419 # line += (jyj)*' ' 3420 # for i in range(ix): 3421 # r = int(100*mapData['rho'][i,j,k]/mapData['rhoMax']) 3422 # line += '%4d'%(r) 3423 # print line+'\n' 3424 ### keep this 3425 3322 3426 3323 3427 def OnTextureRefine(event): 
trunk/GSASIIplot.py
r515 r520 12 12 import os.path 13 13 import numpy as np 14 import numpy.ma as ma 14 15 import numpy.linalg as nl 15 16 import wx … … 435 436 if G2frame.PickId: 436 437 found = [] 437 if G2frame.PatternTree.GetItemText(G2frame.PickId) in ['Index Peak List','Unit Cells List','Reflection Lists']: 438 if G2frame.PatternTree.GetItemText(G2frame.PickId) in ['Index Peak List','Unit Cells List','Reflection Lists'] or \ 439 'PWDR' in G2frame.PatternTree.GetItemText(PickId): 438 440 if len(HKL): 439 441 view = Page.toolbar._views.forward()[0][:2] … … 499 501 else: #picked a limit line 500 502 G2frame.itemPicked = pick 501 elif G2frame.PatternTree.GetItemText(PickId) == 'Reflection Lists': 503 elif G2frame.PatternTree.GetItemText(PickId) == 'Reflection Lists' or \ 504 'PWDR' in G2frame.PatternTree.GetItemText(PickId): 502 505 G2frame.itemPicked = pick 503 506 pick = str(pick) … … 540 543 G2frame.PatternTree.SetItemPyData(PeakId,data) 541 544 G2pdG.UpdatePeakGrid(G2frame,data) 542 elif G2frame.PatternTree.GetItemText(PickId) == 'Reflection Lists' and xpos: 545 elif (G2frame.PatternTree.GetItemText(PickId) == 'Reflection Lists' or \ 546 'PWDR' in G2frame.PatternTree.GetItemText(PickId)) and xpos: 543 547 Phases = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,PatternId,'Reflection Lists')) 544 548 pick = str(G2frame.itemPicked).split('(')[1].strip(')') … … 738 742 else: 739 743 Plot.axvline(hkl[5],color='r',dashes=(5,5)) 740 elif G2frame.PatternTree.GetItemText(PickId) in ['Reflection Lists']: 744 elif G2frame.PatternTree.GetItemText(PickId) in ['Reflection Lists'] or \ 745 'PWDR' in G2frame.PatternTree.GetItemText(PickId): 741 746 Phases = G2frame.PatternTree.GetItemPyData(G2gd.GetPatternTreeItemId(G2frame,PatternId,'Reflection Lists')) 742 747 for pId,phase in enumerate(Phases): … … 2271 2276 drawingData = data['Drawing'] 2272 2277 drawAtoms = drawingData['Atoms'] 2278 mapData = {} 2279 rhoXYZ = [] 2280 if 'Map' in generalData: 2281 mapData = generalData['Map'] 2282 contLevel = drawingData['contourLevel']*mapData['rhoMax'] 2283 if 'deltF' in mapData['MapType']: 2284 rho = ma.array(mapData['rho'],mask=(np.abs(mapData['rho'])<contLevel)) 2285 else: 2286 rho = ma.array(mapData['rho'],mask=(mapData['rho']<contLevel)) 2287 indx = np.array(ma.nonzero(rho)).T 2288 steps = 1./np.array(rho.shape) 2289 rhoXYZ = indx*steps 2273 2290 cx,ct,cs = drawingData['atomPtrs'] 2274 Wt = [255,255,255]2275 Rd = [255,0,0]2276 Gr = [0,255,0]2277 Bl = [0,0,255]2291 Wt = np.array([255,255,255]) 2292 Rd = np.array([255,0,0]) 2293 Gr = np.array([0,255,0]) 2294 Bl = np.array([0,0,255]) 2278 2295 uBox = np.array([[0,0,0],[1,0,0],[1,1,0],[0,1,0],[0,0,1],[1,0,1],[1,1,1],[0,1,1]]) 2279 2296 uEdges = np.array([ … … 2607 2624 glPopMatrix() 2608 2625 2626 def RenderSmallSphere(x,y,z,radius,color): 2627 glMaterialfv(GL_FRONT_AND_BACK,GL_DIFFUSE,color) 2628 glPushMatrix() 2629 glTranslate(x,y,z) 2630 glMultMatrixf(B4mat.T) 2631 q = gluNewQuadric() 2632 gluSphere(q,radius,4,2) 2633 glPopMatrix() 2634 2609 2635 def RenderEllipsoid(x,y,z,ellipseProb,E,R4,color): 2610 2636 s1,s2,s3 = E … … 2691 2717 glEnable(GL_LIGHTING) 2692 2718 glPopMatrix() 2719 2720 def RenderMap(rhoXYZ,indx,rho,cLevel): 2721 for i,xyz in enumerate(rhoXYZ): 2722 x,y,z = xyz 2723 I,J,K = indx[i] 2724 alpha = 1.0 2725 if cLevel < 1.: 2726 alpha = (abs(rho[I,J,K])/mapData['rhoMax']cLevel)/(1.cLevel) 2727 if rho[I,J,K] < 0.: 2728 RenderSmallSphere(x,y,z,0.1*alpha,Rd) 2729 else: 2730 RenderSmallSphere(x,y,z,0.1*alpha,Bl) 2693 2731 2694 2732 def Draw(): … … 2739 2777 BackboneColor = [] 2740 2778 time0 = time.time() 2779 # glEnable(GL_BLEND) 2780 # glBlendFunc(GL_SRC_ALPHA,GL_ONE_MINUS_SRC_ALPHA) 2741 2781 for iat,atom in enumerate(drawingData['Atoms']): 2742 2782 x,y,z = atom[cx:cx+3] … … 2751 2791 if iat in Ind: 2752 2792 color = np.array(Gr)/255. 2793 # color += [.25,] 2753 2794 radius = 0.5 2754 2795 if atom[cs] != '': 2755 glLoadName(atom[3]) 2796 try: 2797 glLoadName(atom[3]) 2798 except: #problem with old files  missing code 2799 pass 2756 2800 if 'balls' in atom[cs]: 2757 2801 vdwScale = drawingData['vdwScale'] … … 2820 2864 elif atom[cs+1] == 'chain' and atom[ct1] == 'CA': 2821 2865 RenderLabel(x,y,z,atom[ct2],radius) 2866 # glDisable(GL_BLEND) 2867 if len(rhoXYZ): 2868 RenderMap(rhoXYZ,indx,rho,drawingData['contourLevel']) 2822 2869 if Backbone: 2823 2870 RenderBackbone(Backbone,BackboneColor,bondR)
Note: See TracChangeset
for help on using the changeset viewer.