Changeset 458
 Timestamp:
 Jan 27, 2012 4:43:30 PM (10 years ago)
 Location:
 trunk
 Files:

 3 edited
Legend:
 Unmodified
 Added
 Removed

trunk/GSASIImath.py
r457 r458 16 16 import math 17 17 import GSASIIpath 18 import GSASIIspc as G2spc 18 19 import scipy.optimize as so 19 20 import scipy.linalg as sl … … 23 24 tand = lambda x: np.tan(x*np.pi/180.) 24 25 asind = lambda x: 180.*np.arcsin(x)/np.pi 26 acosd = lambda x: 180.*np.arccos(x)/np.pi 25 27 atan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi 26 28 … … 143 145 return vcov 144 146 147 def getDistDerv(Oxyz,Txyz,Amat,Tunit,Top,SGData): 148 149 def calcDist(Ox,Tx,U,inv,C,M,T,Amat): 150 TxT = inv*(np.inner(M,Tx)+T)+C 151 TxT = G2spc.MoveToUnitCell(TxT)+U 152 return np.sqrt(np.sum(np.inner(Amat,(TxTOx))**2)) 153 154 inv = 1 155 if Top < 0: 156 inv = 1 157 cent = abs(Top)/100 158 op = abs(Top)%1001 159 M,T = SGData['SGOps'][op] 160 C = SGData['SGCen'][cent] 161 dx = .00001 162 deriv = np.zeros(6) 163 for i in [0,1,2]: 164 Oxyz[i] += dx 165 d0 = calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat) 166 Oxyz[i] = 2*dx 167 deriv[i] = (calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat)d0)/(2.*dx) 168 Oxyz[i] += dx 169 Txyz[i] += dx 170 d0 = calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat) 171 Txyz[i] = 2*dx 172 deriv[i+3] = (calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat)d0)/(2.*dx) 173 Txyz[i] += dx 174 return deriv 175 176 def getAngSig(VA,VB,Amat,SGData,covData={}): 177 178 def calcVec(Ox,Tx,U,inv,C,M,T,Amat): 179 TxT = inv*(np.inner(M,Tx)+T)+C 180 TxT = G2spc.MoveToUnitCell(TxT)+U 181 return np.inner(Amat,(TxTOx)) 182 183 def calcAngle(Ox,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat): 184 VecA = calcVec(Ox,TxA,unitA,invA,CA,MA,TA,Amat) 185 VecA /= np.sqrt(np.sum(VecA**2)) 186 VecB = calcVec(Ox,TxB,unitB,invB,CB,MB,TB,Amat) 187 VecB /= np.sqrt(np.sum(VecB**2)) 188 edge = VecBVecA 189 edge = np.sum(edge**2) 190 angle = (2.edge)/2. 191 angle = max(angle,1.) 192 return acosd(angle) 193 194 OxAN,OxA,TxAN,TxA,unitA,TopA = VA 195 OxBN,OxB,TxBN,TxB,unitB,TopB = VB 196 invA = invB = 1 197 if TopA < 0: 198 invA = 1 199 if TopB < 0: 200 invB = 1 201 centA = abs(TopA)/100 202 centB = abs(TopB)/100 203 opA = abs(TopA)%1001 204 opB = abs(TopB)%1001 205 MA,TA = SGData['SGOps'][opA] 206 MB,TB = SGData['SGOps'][opB] 207 CA = SGData['SGCen'][centA] 208 CB = SGData['SGCen'][centB] 209 if 'covMatrix' in covData: 210 covMatrix = covData['covMatrix'] 211 varyList = covData['varyList'] 212 AngVcov = getVCov(OxAN+TxAN+TxBN,varyList,covMatrix) 213 dx = .00001 214 dadx = np.zeros(9) 215 Ang = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat) 216 for i in [0,1,2]: 217 OxA[i] += dx 218 a0 = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat) 219 OxA[i] = 2*dx 220 dadx[i] = (calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)a0)/dx 221 OxA[i] += dx 222 223 TxA[i] += dx 224 a0 = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat) 225 TxA[i] = 2*dx 226 dadx[i+3] = (calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)a0)/dx 227 TxA[i] += dx 228 229 TxB[i] += dx 230 a0 = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat) 231 TxB[i] = 2*dx 232 dadx[i+6] = (calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)a0)/dx 233 TxB[i] += dx 234 235 sigAng = np.sqrt(np.inner(dadx,np.inner(AngVcov,dadx))) 236 if sigAng < 0.01: 237 sigAng = 0.0 238 return Ang,sigAng 239 else: 240 return calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat),0.0 241 242 145 243 def ValEsd(value,esd=0,nTZ=False): #NOT complete  don't use 146 244 # returns value(esd) string; nTZ=True for no trailing zeros 
trunk/GSASIIphsGUI.py
r457 r458 192 192 pass 193 193 Obj.SetValue("%.3f"%(self.data[item[0]][item[1]])) #reset in case of error 194 195 def getData(self): 196 return self.Data 194 197 195 198 def OnOk(self,event): … … 1059 1062 cn = colLabels.index('Name') 1060 1063 for i,atom in enumerate(atomData): 1061 xyz.append( atom[cn:cn+2]+atom[cx:cx+3])1064 xyz.append([i,]+atom[cn:cn+2]+atom[cx:cx+3]) 1062 1065 if i in indx: 1063 Oxyz.append( atom[cn:cn+2]+atom[cx:cx+3])1066 Oxyz.append([i,]+atom[cn:cn+2]+atom[cx:cx+3]) 1064 1067 DisAglData['OrigAtoms'] = Oxyz 1065 1068 DisAglData['TargAtoms'] = xyz 
trunk/GSASIIstruct.py
r457 r458 26 26 tand = lambda x: np.tan(x*np.pi/180.) 27 27 asind = lambda x: 180.*np.arcsin(x)/np.pi 28 acosd = lambda x: 180.*np.arccos(x)/np.pi 28 29 atan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi 29 30 … … 2569 2570 2570 2571 Amat,Bmat = G2lat.cell2AB(Cell[:6]) 2572 covData = {} 2571 2573 if 'covData' in DisAglData: 2572 2574 covData = DisAglData['covData'] 2575 covMatrix = covData['covMatrix'] 2576 varyList = covData['varyList'] 2573 2577 pfx = str(DisAglData['pId'])+'::' 2574 2578 A = G2lat.cell2A(Cell[:6]) … … 2580 2584 line += name+vals 2581 2585 print line 2582 else: 2586 else: 2583 2587 print '\n Unit cell: a = ','%.5f'%(Cell[0]),' b = ','%.5f'%(Cell[1]),' c = ','%.5f'%(Cell[2]), \ 2584 2588 ' alpha = ','%.3f'%(Cell[3]),' beta = ','%.3f'%(Cell[4]),' gamma = ', \ … … 2590 2594 [0,1,1],[0,1,0],[0,1,1],[0,0,1],[0,0,0],[0,0,1],[0,1,1],[0,1,0],[0,1,1], 2591 2595 [1,1,1],[1,1,0],[1,1,1],[1,0,1],[1,0,0],[1,0,1],[1,1,1],[1,1,0],[1,1,1]]) 2592 origAtoms = zip(DisAglData['OrigIndx'],DisAglData['OrigAtoms'])2596 origAtoms = DisAglData['OrigAtoms'] 2593 2597 targAtoms = DisAglData['TargAtoms'] 2594 2598 for Oatom in origAtoms: 2599 OxyzNames = '' 2600 IndBlist = [] 2595 2601 Dist = [] 2596 2602 Vect = [] 2603 VectA = [] 2604 angles = [] 2597 2605 for Tatom in targAtoms: 2598 result = G2spc.GenAtom(Tatom[2:5],SGData,False) 2599 BsumR = (Radii[Oatom[1][1]][0]+Radii[Tatom[1]][0])*Factor[0] 2600 AsumR = (Radii[Oatom[1][1]][1]+Radii[Tatom[1]][1])*Factor[1] 2606 Xvcov = [] 2607 TxyzNames = '' 2608 if 'covData' in DisAglData: 2609 OxyzNames = [pfx+'dAx:%d'%(Oatom[0]),pfx+'dAy:%d'%(Oatom[0]),pfx+'dAz:%d'%(Oatom[0])] 2610 TxyzNames = [pfx+'dAx:%d'%(Tatom[0]),pfx+'dAy:%d'%(Tatom[0]),pfx+'dAz:%d'%(Tatom[0])] 2611 Xvcov = G2mth.getVCov(OxyzNames+TxyzNames,varyList,covMatrix) 2612 result = G2spc.GenAtom(Tatom[3:6],SGData,False) 2613 BsumR = (Radii[Oatom[2]][0]+Radii[Tatom[2]][0])*Factor[0] 2614 AsumR = (Radii[Oatom[2]][1]+Radii[Tatom[2]][1])*Factor[1] 2601 2615 for Txyz,Top,Tunit in result: 2602 Dx = (Txyznp.array(Oatom[ 1][2:5]))+Units2616 Dx = (Txyznp.array(Oatom[3:6]))+Units 2603 2617 dx = np.inner(Amat,Dx) 2604 2618 dist = ma.masked_less(np.sqrt(np.sum(dx**2,axis=0)),0.5) 2605 2619 IndB = ma.nonzero(ma.masked_greater(distBsumR,0.)) 2606 IndA = ma.nonzero(ma.masked_greater(distAsumR,0.))2607 2620 if np.any(IndB): 2608 2621 for indb in IndB: 2609 Dist.append([Oatom[1][0],Tatom[0],Tunit,Top,ma.getdata(dist[indb])]) 2610 Vect += dx.T[indb] 2611 2612 for dist in Dist: 2613 print dist 2622 for i in range(len(indb)): 2623 if str(dx.T[indb][i]) not in IndBlist: 2624 IndBlist.append(str(dx.T[indb][i])) 2625 unit = Units[indb][i] 2626 tunit = '[%2d%2d%2d]'%(unit[0],unit[1],unit[2]) 2627 pdpx = G2mth.getDistDerv(Oatom[3:6],Tatom[3:6],Amat,unit,Top,SGData) 2628 sig = 0.0 2629 if len(Xvcov): 2630 sig = np.sqrt(np.inner(pdpx,np.inner(Xvcov,pdpx))) 2631 Dist.append([Oatom[1],Tatom[1],tunit,Top,ma.getdata(dist[indb])[i],sig]) 2632 if (Dist[1][1]AsumR) <= 0.: 2633 Vect.append(dx.T[indb][i]/Dist[1][2]) 2634 VectA.append([OxyzNames,np.array(Oatom[3:6]),TxyzNames,np.array(Tatom[3:6]),unit,Top]) 2635 else: 2636 Vect.append([0.,0.,0.]) 2637 VectA.append([]) 2638 Vect = np.array(Vect) 2639 angles = np.zeros((len(Vect),len(Vect))) 2640 angsig = np.zeros((len(Vect),len(Vect))) 2641 for i,veca in enumerate(Vect): 2642 if np.any(veca): 2643 for j,vecb in enumerate(Vect): 2644 if np.any(vecb): 2645 angles[i][j],angsig[i][j] = G2mth.getAngSig(VectA[i],VectA[j],Amat,SGData,covData) 2646 line = '' 2647 for i,x in enumerate(Oatom[3:6]): 2648 if len(Xvcov): 2649 line += '%12s'%(G2mth.ValEsd(x,np.sqrt(Xvcov[i][i]),True)) 2650 else: 2651 line += '%12.5f'%(x) 2652 print '\n Distances & angles for ',Oatom[1],' at ',line 2653 print 80*'*' 2654 line = '' 2655 for dist in Dist[:1]: 2656 line += '%12s'%(dist[1].center(12)) 2657 print ' To cell +(sym. op.) dist. ',line 2658 for i,dist in enumerate(Dist): 2659 line = '' 2660 for j,angle in enumerate(angles[i][0:i]): 2661 sig = angsig[i][j] 2662 if angle: 2663 if sig: 2664 line += '%12s'%(G2mth.ValEsd(angle,sig,True).center(12)) 2665 else: 2666 val = '%.3f'%(angle) 2667 line += '%12s'%(val.center(12)) 2668 else: 2669 line += 12*' ' 2670 if dist[5]: #sig exists! 2671 val = G2mth.ValEsd(dist[4],dist[5]) 2672 else: 2673 val = '%8.4f'%(dist[4]) 2674 print ' %8s%10s+(%4d) %12s'%(dist[1].ljust(8),dist[2].ljust(10),dist[3],val.center(12)),line 2675 2614 2676 2615 for vect in Vect:2616 print vect2617 2618 2619 2620 # def FindBonds(): #uses numpy & masks  very fast even for proteins!2621 # import numpy.ma as ma2622 # atomData = data['Atoms']2623 # generalData = data['General']2624 # Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])2625 # radii = generalData['BondRadii']2626 # atomTypes = generalData['AtomTypes']2627 # try:2628 # indH = atomTypes.index('H')2629 # radii[indH] = 0.52630 # except:2631 # pass2632 # for atom in atomData:2633 # atom[2] = [] #clear out old bond/angle tables2634 # atom[1] = []2635 # Indx = range(len(atomData))2636 # Atoms = []2637 # Radii = []2638 # for atom in atomData:2639 # Atoms.append(np.array(atom[cx:cx+3]))2640 # try:2641 # if not hydro and atom[ct] == 'H':2642 # Radii.append(0.0)2643 # else:2644 # Radii.append(radii[atomTypes.index(atom[ct])])2645 # except ValueError: #changed atom type!2646 # Radii.append(0.20)2647 # Atoms = np.array(Atoms)2648 # Radii = np.array(Radii)2649 # IASR = zip(Indx,Atoms,Radii)2650 # for atomA in IASR:2651 # Dx = AtomsatomA[1]2652 # dist = ma.masked_less(np.sqrt(np.sum(np.inner(Amat,Dx)**2,axis=0)),0.5) #gets rid of self & disorder "bonds" < 0.5A2653 # sumR = atomA[3]+Radii2654 # IndB = ma.nonzero(ma.masked_greater(distradiusFactor*sumR,0.)) #get indices of bonded atoms2655 # i = atomA[0]2656 # for j in IndB[0]:2657 # if j > i:2658 # atomData[i][2].append(Dx[j]*Radii[i]/sumR[j])2659 # atomData[j][2].append(Dx[j]*Radii[j]/sumR[j])2660 #2661 2662 2677 def main(): 2663 2678 arg = sys.argv
Note: See TracChangeset
for help on using the changeset viewer.