Changeset 5460
- Timestamp:
- Dec 31, 2022 12:58:26 PM (11 months ago)
- Location:
- trunk
- Files:
-
- 5 edited
-
GSASIIElem.py (modified) (1 diff)
-
GSASIIlattice.py (modified) (3 diffs)
-
GSASIIphsGUI.py (modified) (1 diff)
-
GSASIIplot.py (modified) (4 diffs)
-
GSASIIstrIO.py (modified) (4 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/GSASIIElem.py
r5454 r5460 728 728 if len(landeg) < len(generalData['AtomTypes']): 729 729 landeg.append(2.0) 730 if 'Q' in atom[ct]: #RBData must exist #TODO: need code here to add to atom count for atom type in spinrb730 if 'Q' in atom[ct]: 731 731 for Srb in RBModels.get('Spin',[]): 732 if Srb['Ids'][0] != atom[cia+8]: 733 continue 732 734 nSh = len(Srb['RBId']) 733 735 for iSh in range(nSh): 734 Info = G2elem.GetAtomInfo(Srb['atType'][iSh]) 735 if Info['Symbol'] not in generalData['AtomTypes']: 736 generalData['AtomTypes'].append(Info['Symbol']) 737 generalData['Z'] = Info['Z'] 738 generalData['Isotopes'][Info['Symbol']] = Info['Isotopes'] 739 generalData['BondRadii'].append(Info['Drad']) 740 generalData['AngleRadii'].append(Info['Arad']) 741 generalData['vdWRadii'].append(Info['Vdrad']) 742 if Info['Symbol'] in generalData['Isotope']: 743 if generalData['Isotope'][Info['Symbol']] not in generalData['Isotopes'][Info['Symbol']]: 744 isotope = list(generalData['Isotopes'][Info['Symbol']].keys())[-1] 745 generalData['Isotope'][Info['Symbol']] = isotope 746 generalData['AtomMass'].append(Info['Isotopes'][generalData['Isotope'][Info['Symbol']]]['Mass']) 747 else: 748 generalData['Isotope'][Info['Symbol']] = 'Nat. Abund.' 749 if 'Nat. Abund.' not in generalData['Isotopes'][Info['Symbol']]: 750 isotope = list(generalData['Isotopes'][Info['Symbol']].keys())[-1] 751 generalData['Isotope'][Info['Symbol']] = isotope 752 generalData['AtomMass'].append(Info['Mass']) 753 generalData['NoAtoms'][Info['Symbol']] = atom[cx+3]*atom[cs+1]*Srb['Natoms'][iSh] 754 generalData['Color'].append(Info['Color']) 736 Info = G2elem.GetAtomInfo(Srb['atType'][iSh]) 737 if Info['Symbol'] not in generalData['AtomTypes']: 738 generalData['AtomTypes'].append(Info['Symbol']) 739 generalData['Z'] = Info['Z'] 740 generalData['Isotopes'][Info['Symbol']] = Info['Isotopes'] 741 generalData['BondRadii'].append(Info['Drad']) 742 generalData['AngleRadii'].append(Info['Arad']) 743 generalData['vdWRadii'].append(Info['Vdrad']) 744 if Info['Symbol'] in generalData['Isotope']: 745 if generalData['Isotope'][Info['Symbol']] not in generalData['Isotopes'][Info['Symbol']]: 746 isotope = list(generalData['Isotopes'][Info['Symbol']].keys())[-1] 747 generalData['Isotope'][Info['Symbol']] = isotope 748 generalData['AtomMass'].append(Info['Isotopes'][generalData['Isotope'][Info['Symbol']]]['Mass']) 755 749 else: 756 generalData['NoAtoms'][Info['Symbol']] += atom[cx+3]*atom[cs+1]*Srb['Natoms'][iSh] 757 # break 758 # else: 759 # break 750 generalData['Isotope'][Info['Symbol']] = 'Nat. Abund.' 751 if 'Nat. Abund.' not in generalData['Isotopes'][Info['Symbol']]: 752 isotope = list(generalData['Isotopes'][Info['Symbol']].keys())[-1] 753 generalData['Isotope'][Info['Symbol']] = isotope 754 generalData['AtomMass'].append(Info['Mass']) 755 generalData['NoAtoms'][Info['Symbol']] = atom[cx+3]*atom[cs+1]*Srb['Natoms'][iSh] 756 generalData['Color'].append(Info['Color']) 757 else: 758 generalData['NoAtoms'][Info['Symbol']] += atom[cx+3]*atom[cs+1]*Srb['Natoms'][iSh] 760 759 761 760 if generalData['Type'] == 'magnetic': -
trunk/GSASIIlattice.py
r5453 r5460 2614 2614 return Th,Ph 2615 2615 2616 def SHarmcal(SytSym,SHFln,psi,gam): 2617 '''Perform a surface spherical harmonics computation. 2618 Note that the the number of gam values must either be 1 or must match psi 2619 2620 :param str SytSym: sit symmetry - only looking for cubics 2621 :param dict SHFln: spherical harmonics coefficients; key has L & M 2622 :param float/array psi: Azimuthal coordinate 0 <= Th <= 360 2623 :param float/array gam: Polar coordinate 0<= Ph <= 180 2624 2625 :returns array SHVal: spherical harmonics array for psi,gam values 2626 ''' 2627 SHVal = np.ones_like(psi) 2628 for term in SHFln: 2629 if 'C(' in term[:2]: 2630 l,m = eval(term.strip('C')) 2631 if SytSym in ['m3m','m3','43m','432','23']: 2632 Ksl = CubicSHarm(l,m,psi,gam) 2633 else: 2634 p = SHFln[term][2] 2635 Ksl = SphHarmAng(l,m,p,psi,gam) 2636 SHVal += SHFln[term][0]*Ksl 2637 return SHVal 2638 2616 2639 def SphHarmAng(L,M,P,Th,Ph): 2617 2640 ''' Compute spherical harmonics values using scipy.special.sph_harm … … 2629 2652 2630 2653 if M == 0: 2631 return np.real(ylmp)/ 4.2654 return np.real(ylmp)/SQ2 2632 2655 if P>0: 2633 2656 return np.real(ylmp) … … 2818 2841 return PolVal 2819 2842 2820 def SHarmcal(SytSym,SHFln,psi,gam):2821 '''Perform a surface spherical harmonics computation.2822 Note that the the number of gam values must either be 1 or must match psi2823 2824 :param str SytSym: sit symmetry - only looking for cubics2825 :param dict SHFln: spherical harmonics coefficients; key has L & M2826 :param float/array psi: Azimuthal coordinate 0 <= Th <= 3602827 :param float/array gam: Polar coordinate 0<= Ph <= 1802828 2829 :returns array SHVal: spherical harmonics array for psi,gam values2830 '''2831 SHVal = np.ones_like(gam)2832 for term in SHFln:2833 if 'C(' in term[:2]:2834 l,m = eval(term.strip('C'))2835 if SytSym in ['m3m','m3','43m','432','23']:2836 Ksl = CubicSHarm(l,m,psi,gam)2837 else:2838 p = SHFln[term][2]2839 Ksl = SphHarmAng(l,m,p,psi,gam)2840 SHVal += SHFln[term][0]*Ksl2841 return SHVal2842 2843 2843 def invpolfcal(ODFln,SGData,phi,beta): 2844 2844 'needs doc string' -
trunk/GSASIIphsGUI.py
r5454 r5460 11554 11554 data['RBModels']['Spin'][-1][name].append(rbData[name]) 11555 11555 wx.CallAfter(FillRigidBodyGrid,True,spnId=rbId) 11556 11557 11556 11558 11557 def SHsizer(): -
trunk/GSASIIplot.py
r5453 r5460 9782 9782 9783 9783 def RenderTextureSphere(x,y,z,radius,Qmat,color,shape=[20,10],Fade=None): 9784 global spID9785 9784 SpFade = np.zeros(list(Fade.shape)+[4,],dtype=np.dtype('B')) 9786 9785 SpFade[:,:,:3] = Fade[:,:,nxs]*list(color) 9787 9786 SpFade[:,:,3] = 128 9788 newSP = False9789 9787 spID = GL.glGenTextures(1) 9790 # if not spID:9791 # spID = GL.glGenTextures(1)9792 # newSP = True9793 9788 GL.glPixelStorei(GL.GL_UNPACK_ALIGNMENT, 1) 9794 9789 GL.glEnable(GL.GL_BLEND) 9790 GL.glFrontFace(GL.GL_CCW) #shows outside 9791 GL.glEnable(GL.GL_CULL_FACE) #removes striping 9795 9792 GL.glBlendFunc(GL.GL_SRC_ALPHA,GL.GL_ONE_MINUS_SRC_ALPHA) 9796 9793 GL.glEnable(GL.GL_TEXTURE_2D) 9797 9794 GL.glBindTexture(GL.GL_TEXTURE_2D, spID) 9798 GL.glFrontFace(GL.GL_CCW) #shows outside9799 GL.glEnable(GL.GL_CULL_FACE) #removes striping9800 9795 GL.glTexEnvf(GL.GL_TEXTURE_ENV, GL.GL_TEXTURE_ENV_MODE, GL.GL_REPLACE) 9801 9796 GL.glTexEnvf(GL.GL_TEXTURE_ENV, GL.GL_ALPHA_SCALE, 1.0) 9802 GL.glTexParameterf(GL.GL_TEXTURE_2D, GL.GL_TEXTURE_WRAP_S, GL.GL_CLAMP)9803 GL.glTexParameterf(GL.GL_TEXTURE_2D, GL.GL_TEXTURE_WRAP_T, GL.GL_CLAMP)9804 GL.glTexParameterf(GL.GL_TEXTURE_2D, GL.GL_TEXTURE_WRAP_S, GL.GL_REPEAT)9805 GL.glTexParameterf(GL.GL_TEXTURE_2D, GL.GL_TEXTURE_WRAP_T, GL.GL_REPEAT)9806 9797 GL.glTexParameteri(GL.GL_TEXTURE_2D, GL.GL_TEXTURE_BASE_LEVEL, 0) 9807 9798 GL.glTexParameteri(GL.GL_TEXTURE_2D, GL.GL_TEXTURE_MAX_LEVEL, 0) … … 9809 9800 GL.glTexParameter(GL.GL_TEXTURE_2D,GL.GL_TEXTURE_MAG_FILTER,GL.GL_LINEAR) 9810 9801 GL.glTexImage2D(GL.GL_TEXTURE_2D,0,GL.GL_RGBA,shape[0], shape[1],0,GL.GL_RGBA,GL.GL_UNSIGNED_BYTE,SpFade) 9811 # if newSP:9812 # GL.glTexImage2D(GL.GL_TEXTURE_2D,0,GL.GL_RGBA,shape[0], shape[1],0,GL.GL_RGBA,GL.GL_UNSIGNED_BYTE,SpFade)9813 # else:9814 # GL.glTexSubImage2D(GL.GL_TEXTURE_2D,0,0,0,shape[0], shape[1],GL.GL_RGBA,GL.GL_UNSIGNED_BYTE,SpFade)9815 9802 q = GLU.gluNewQuadric() 9816 9803 GLU.gluQuadricDrawStyle(q,GLU.GLU_FILL) … … 10309 10296 QR,R = G2mth.make2Quat(V,np.array([0.,0.,1.0])) 10310 10297 QA = G2mth.AVdeg2Q(A,np.array([0.,0.,1.0])) 10311 Q2 = G2mth.prodQQ(QA,QR) 10298 Q2 = G2mth.prodQQ(QA,QR) #correct - rotates about V axis 10312 10299 Qmat = G2mth.Q2Mat(Q2) 10313 Q4mat = np.concatenate((np.concatenate((Qmat,[[0],[0],[0]]),axis=1),[[0,0,0,1],]),axis=0) 10300 Q4mat = np.eye(4) 10301 Q4mat[:3,:3] = Qmat 10302 Npsi,Ngam = 120,60 10303 PSI,GAM = np.mgrid[0:Npsi,0:Ngam] #[azm,pol] 10304 PSI = PSI.flatten()*360./Npsi #azimuth 0-360 ncl 10305 GAM = GAM.flatten()*180./Ngam #polar 0-180 incl 10314 10306 for ish,nSH in enumerate(SpnData['nSH']): 10315 10307 if nSH > 0: 10316 10308 SHC = SpnData['SHC'][ish] 10317 PSI,GAM = np.mgrid[0:120,0:60] #[azm,pol] 10318 PSI = PSI.flatten()*3. #azimuth 0-360 ncl 10319 GAM = GAM.flatten()*3. #polar 0-180 incl 10320 P = G2lat.SHarmcal(SytSym,SHC,GAM,PSI).reshape((120,60)) 10309 P = G2lat.SHarmcal(SytSym,SHC,PSI,GAM).reshape((Npsi,Ngam)) 10321 10310 if np.min(P) < np.max(P): 10322 10311 P = (P-np.min(P))/(np.max(P)-np.min(P)) 10323 RenderTextureSphere(x,y,z,radius[ish][0],Q4mat,atColor[ish],shape=[ 120,60],Fade=P)10312 RenderTextureSphere(x,y,z,radius[ish][0],Q4mat,atColor[ish],shape=[Npsi,Ngam],Fade=P.T) 10324 10313 else: 10325 10314 RenderSphere(x,y,z,radius[ish][0],atColor[ish],fade,shape=[60,30]) … … 10535 10524 10536 10525 #### PlotStructure starts here 10537 global mcsaXYZ,mcsaTypes,mcsaBonds,txID,contourSet ,spID10526 global mcsaXYZ,mcsaTypes,mcsaBonds,txID,contourSet 10538 10527 global cell, Vol, Amat, Bmat, A4mat, B4mat, BondRadii 10539 10528 txID = -1 10540 spID = 010541 10529 ForthirdPI = 4.0*math.pi/3.0 10542 10530 RBId = G2gd.GetGPXtreeItemId(G2frame, G2frame.root, 'Rigid bodies') -
trunk/GSASIIstrIO.py
r5443 r5460 1297 1297 Angle = 2.0*acosd(Qrijk[0]) 1298 1298 pFile.write('\nRBObject %s at %10.4f %10.4f %10.4f Refine? %s\n'% 1299 (RB['RBname'] ,Oxyz[0],Oxyz[1],Oxyz[2],RB['Orig'][1]))1299 (RB['RBname'][0],Oxyz[0],Oxyz[1],Oxyz[2],RB['Orig'][1])) 1300 1300 pFile.write('Orientation angle,vector: %10.3f %10.4f %10.4f %10.4f Refine? %s\n'% 1301 1301 (Angle,Qrijk[1],Qrijk[2],Qrijk[3],RB['Orient'][1])) 1302 pFile.write('Bessel/Spherical Harmonics coefficients; symmetry required sign shown') 1303 SHC = RB['SHC'] 1304 if len(SHC): 1305 SHkeys = list(SHC.keys()) 1306 nCoeff = len(SHC) 1307 nBlock = nCoeff//6+1 1308 iBeg = 0 1309 iFin = min(iBeg+6,nCoeff) 1310 for block in range(nBlock): 1311 ptlbls = ' names :' 1312 ptstr = ' values:' 1313 ptref = ' refine:' 1314 for item in SHkeys[iBeg:iFin]: 1315 ptlbls += '%12s'%(item) 1316 ptstr += '%9.4f*%2.0f'%(SHC[item][0],SHC[item][1]) 1317 ptref += '%12s'%(SHC[item][2]) 1318 pFile.write(ptlbls+'\n') 1319 pFile.write(ptstr+'\n') 1320 pFile.write(ptref+'\n') 1321 iBeg += 6 1302 pFile.write('Bessel/Spherical Harmonics coefficients; symmetry required sign shown\n') 1303 for ish,SHC in enumerate(RB['SHC']): 1304 if len(SHC): 1305 pFile.write('Spin shell (%s) no: %d\n'%(RB['RBname'][ish],ish)) 1306 SHkeys = list(SHC.keys()) 1307 nCoeff = len(SHC) 1308 nBlock = nCoeff//6+1 1309 iBeg = 0 1322 1310 iFin = min(iBeg+6,nCoeff) 1311 for block in range(nBlock): 1312 ptlbls = ' names :' 1313 ptstr = ' values:' 1314 ptref = ' refine:' 1315 for item in SHkeys[iBeg:iFin]: 1316 ptlbls += '%12s'%(item) 1317 ptstr += '%9.4f*%2.0f'%(SHC[item][0],SHC[item][1]) 1318 ptref += '%12s'%(SHC[item][2]) 1319 pFile.write(ptlbls+'\n') 1320 pFile.write(ptstr+'\n') 1321 pFile.write(ptref+'\n') 1322 iBeg += 6 1323 iFin = min(iBeg+6,nCoeff) 1323 1324 1324 1325 def PrintAtoms(General,Atoms): … … 1440 1441 if 'AtomFrac' not in RB and rbKey != 'S': raise Exception('out of date RB: edit in RB Models') 1441 1442 # end patch 1442 rbid = str(rbids.index(RB['RBId'])) 1443 if rbKey == 'S': 1444 rbid = str(rbids.index(RB['RBId'][0])) #for spin RBs 1445 else: 1446 rbid = str(rbids.index(RB['RBId'])) #others 1443 1447 pfxRB = pfx+'RB'+rbKey+'P' 1444 1448 pstr = ['x','y','z'] … … 1485 1489 name = pfx+'RB'+rbKey+'AtNo:'+str(iRB)+':'+rbid 1486 1490 phaseDict[name] = atomIndx[RB['Ids'][0]][1] 1487 name = pfx+'RB'+rbKey+'AtType:'+str(iRB)+':'+rbid1488 phaseDict[name] = RB['atType']1489 1491 name = pfx+'RB'+rbKey+'SytSym:'+str(iRB)+':'+rbid 1490 1492 phaseDict[name] = RB['SytSym'] 1493 for ish,atype in enumerate(RB['atType']): 1494 name = '%sRBSAtType;%d:%d:%s'%(pfx,ish,iRB,rbid) 1495 phaseDict[name] = RB['atType'][ish] 1491 1496 1492 1497 def MakeRBThermals(rbKey,phaseVary,phaseDict): … … 1531 1536 1532 1537 def MakeRBSphHarm(rbKey,phaseVary,phaseDict): 1533 rbid = str(rbids.index(RB['RBId'])) 1534 pfxRB = pfx+'RB'+rbKey+'Sh;' 1535 for i,shcof in enumerate(RB['SHC']): 1536 SHcof = RB['SHC'][shcof] 1537 name = pfxRB+shcof+':'+str(iRB)+':'+rbid 1538 phaseDict[name] = SHcof[0]*SHcof[1] 1539 if SHcof[2]: 1540 phaseVary += [name,] 1538 rbid = str(rbids.index(RB['RBId'][0])) 1539 for ish,Shcof in enumerate(RB['SHC']): 1540 pfxRB = '%sRBSSh;%d;'%(pfx,ish) 1541 for i,shcof in enumerate(Shcof): 1542 SHcof = Shcof[shcof] 1543 name = pfxRB+shcof+':'+str(iRB)+':'+rbid 1544 phaseDict[name] = SHcof[0]*SHcof[1] 1545 if SHcof[2]: 1546 phaseVary += [name,] 1541 1547 1542 1548 if Print and pFile is None: raise Exception("specify pFile or Print=False")
Note: See TracChangeset
for help on using the changeset viewer.
