Changeset 477

Timestamp:

Feb 9, 2012 1:23:01 PM (12 years ago)

Author:

vondreele

Message:

final fixes to dist-angle-torsion calc.

Location:

trunk

Files:

• GSASIImath.py (modified) (3 diffs)
• GSASIIphsGUI.py (modified) (4 diffs)
• GSASIIplot.py (modified) (2 diffs)
• GSASIIspc.py (modified) (2 diffs)
• GSASIIstruct.py (modified) (5 diffs)

trunk/GSASIImath.py

@@ -148 +148 @@
     
     def calcDist(Ox,Tx,U,inv,C,M,T,Amat):
-        TxT = inv*(np.inner(M,Tx)+T)+C
-        TxT = G2spc.MoveToUnitCell(TxT)+U
+        TxT = inv*(np.inner(M,Tx)+T)+C+U
         return np.sqrt(np.sum(np.inner(Amat,(TxT-Ox))**2))
         
@@ -236 +235 @@
         return calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat),0.0
         
-def GetTorsionSig(Atoms,Amat,SGData,covData={}):
-    XYZ = []
-    for atom in Atoms:    
-        XYZ.append(np.array(atom[3:6]))
-    XYZ = np.inner(Amat,XYZ).T
-    V1 = XYZ[1]-XYZ[0]
-    V2 = XYZ[2]-XYZ[1]
-    V3 = XYZ[3]-XYZ[2]
-    V1 /= np.sqrt(np.sum(V1**2))
-    V2 /= np.sqrt(np.sum(V2**2))
-    V3 /= np.sqrt(np.sum(V3**2))
-    M = np.array([V1,V2,V3])
-    D = nl.det(M)
-    Ang = 1.0
-    P12 = np.dot(V1,V2)
-    P13 = np.dot(V1,V3)
-    P23 = np.dot(V2,V3)
-    Tors = acosd((P12*P23-P13)/(np.sqrt(1.-P12**2)*np.sqrt(1.-P23**2)))*D/abs(D)
-    
-    sig = 0.0
+def GetTorsionSig(Oatoms,Atoms,Amat,SGData,covData={}):
+    
+    def calcTorsion(Atoms,SyOps,Amat):
+        
+        XYZ = []
+        for i,atom in enumerate(Atoms):
+            Inv,M,T,C,U = SyOps[i]
+            XYZ.append(np.array(atom[1:4]))
+            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
+            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
+        V1 = XYZ[1]-XYZ[0]
+        V2 = XYZ[2]-XYZ[1]
+        V3 = XYZ[3]-XYZ[2]
+        V1 /= np.sqrt(np.sum(V1**2))
+        V2 /= np.sqrt(np.sum(V2**2))
+        V3 /= np.sqrt(np.sum(V3**2))
+        M = np.array([V1,V2,V3])
+        D = nl.det(M)
+        Ang = 1.0
+        P12 = np.dot(V1,V2)
+        P13 = np.dot(V1,V3)
+        P23 = np.dot(V2,V3)
+        Tors = acosd((P12*P23-P13)/(np.sqrt(1.-P12**2)*np.sqrt(1.-P23**2)))*D/abs(D)
+        return Tors
+            
+    Inv = []
+    SyOps = []
+    names = []
+    for i,atom in enumerate(Oatoms):
+        names += atom[-1]
+        Op,unit = Atoms[i][-1]
+        inv = Op/abs(Op)
+        m,t = SGData['SGOps'][abs(Op)%100-1]
+        c = SGData['SGCen'][abs(Op)/100]
+        SyOps.append([inv,m,t,c,unit])
+    Tors = calcTorsion(Oatoms,SyOps,Amat)
+    
+    sig = -0.01
     if 'covMatrix' in covData:
-        for atom in Atoms:
-            xyz = atom[3:6]
-            Op,Unit = atom[-1]
-            inv = Op/abs(Op)
-            M,T = SGData['SGOps'][abs(Op)%100-1]
-            C = SGData['SGCen'][abs(Op)/100]
-            #reverse inv*(np.inner(M,Tx)+T)+C
-            XYZ = np.inner(nl.inv(M),((xyz-C)*inv-T))
-            print Op,Unit,xyz,XYZ
+        parmNames = []
+        dx = .00001
+        dadx = np.zeros(12)
+        ang = calcTorsion(Oatoms,SyOps,Amat)
+        for i in range(12):
+            ia = i/3
+            ix = i%3
+            Oatoms[ia][ix+1] += dx
+            a0 = calcTorsion(Oatoms,SyOps,Amat)
+            Oatoms[ia][ix+1] -= 2*dx
+            dadx[i] = (calcTorsion(Oatoms,SyOps,Amat)-a0)/(2.*dx)
         covMatrix = covData['covMatrix']
         varyList = covData['varyList']
+        TorVcov = getVCov(names,varyList,covMatrix)
+        sig = np.sqrt(np.inner(dadx,np.inner(TorVcov,dadx)))
+        if sig < 0.01:
+            sig = -0.01
     
     return Tors,sig
@@ -286 +309 @@
         return str(fmt%(value,int(round(esd*10**(mdec(esd)))))).strip('"')
     elif esd < 0:
-         return str(round(value,mdec(esd)))
+         return str(round(value,mdec(esd)-1))
     else:
         text = str("%f"%(value))

trunk/GSASIIphsGUI.py

@@ -1648 +1648 @@
                 if atom[cuia] == 'A':
                     Uij = atom[cuij:cuij+6]
-                    result = G2spc.GenAtom(XYZ,SGData,False,Uij)
+                    result = G2spc.GenAtom(XYZ,SGData,False,Uij,False)
                     for item in result:
                         atom = copy.copy(atomData[ind])
@@ -1665 +1665 @@
                                 atomData.append(atom[:])
                 else:
-                    result = G2spc.GenAtom(XYZ,SGData,False)
+                    result = G2spc.GenAtom(XYZ,SGData,False,Move=False)
                     for item in result:
                         atom = copy.copy(atomData[ind])
@@ -1802 +1802 @@
         event.StopPropagation()
         
-    def FindAtomIndexByIDs(atomData,IDs):
+    def FindAtomIndexByIDs(atomData,IDs,Draw=True):
         indx = []
         for i,atom in enumerate(atomData):
-            if atom[-2] in IDs:
+            if Draw and atom[-3] in IDs:
+                indx.append(i)
+            elif atom[-1] in IDs:
                 indx.append(i)
         return indx
@@ -1835 +1837 @@
             return
         TorsionData = {}
+        ocx,oct,ocs,cia = data['General']['AtomPtrs']
         drawingData = data['Drawing']
         atomData = data['Atoms']
+#        for atom in atomData: print atom
         atomDData = drawingData['Atoms']
+#        for atom in atomDData: print atom
         colLabels = [drawAtoms.GetColLabelValue(c) for c in range(drawAtoms.GetNumberCols())]
         cx = colLabels.index('x')
         cn = colLabels.index('Name')
+        cid = colLabels.index('I/A')+8
         xyz = []
+        Oxyz = []
         for i,atom in enumerate(atomDData):
             if i in indx:
                 xyz.append([i,]+atom[cn:cn+2]+atom[cx:cx+4]) #also gets Sym Op
+                id = FindAtomIndexByIDs(atomData,[atom[cid],],False)[0]
+                Oxyz.append([id,]+atomData[id][cx+1:cx+4])
         TorsionData['Datoms'] = xyz
+        TorsionData['Oatoms'] = Oxyz
         generalData = data['General']
         TorsionData['Name'] = generalData['Name']

trunk/GSASIIplot.py

@@ -2266 +2266 @@
             return G2frame.dataDisplay.GetSelection()
         except AttributeError:
-            print G2frame.dataDisplay.GetLabel()
             G2frame.G2plotNB.status.SetStatusText('Select this from Phase data window!')
             return 0
@@ -2572 +2571 @@
                             
     def Draw():
-        import numpy.linalg as nl
         Ind = GetSelectedAtoms()
         VS = np.array(Page.canvas.GetSize())

trunk/GSASIIspc.py

@@ -252 +252 @@
     return new
         
-def GenAtom(XYZ,SGData,All=False,Uij=[]):
+def GenAtom(XYZ,SGData,All=False,Uij=[],Move=True):
     '''
     Generates the equivalent positions for a specified coordinate and space group
@@ -273 +273 @@
     Cell = []
     X = np.array(XYZ)
-    X = MoveToUnitCell(X)
+    if Move:
+        X = MoveToUnitCell(X)
     for ic,cen in enumerate(SGData['SGCen']):
         C = np.array(cen)

trunk/GSASIIstruct.py

@@ -2610 +2610 @@
                 TxyzNames = [pfx+'dAx:%d'%(Tatom[0]),pfx+'dAy:%d'%(Tatom[0]),pfx+'dAz:%d'%(Tatom[0])]
                 Xvcov = G2mth.getVCov(OxyzNames+TxyzNames,varyList,covMatrix)
-            result = G2spc.GenAtom(Tatom[3:6],SGData,False)
+            result = G2spc.GenAtom(Tatom[3:6],SGData,False,Move=False)
             BsumR = (Radii[Oatom[2]][0]+Radii[Tatom[2]][0])*Factor[0]
             AsumR = (Radii[Oatom[2]][1]+Radii[Tatom[2]][1])*Factor[1]
@@ -2679 +2679 @@
         print 80*'*'
         print '   Torsion angle for phase '+name
-        print 80*'*','\n'
+        print 80*'*'
 
     ShowBanner(TorsionData['Name'])
     SGData = TorsionData['SGData']
-    SGtext = G2spc.SGPrint(SGData)
-    for line in SGtext: print line
     Cell = TorsionData['Cell']
     
     Amat,Bmat = G2lat.cell2AB(Cell[:6])
     covData = {}
+    pfx = ''
     if 'covData' in TorsionData:   
         covData = TorsionData['covData']
@@ -2694 +2693 @@
         varyList = covData['varyList']
         pfx = str(TorsionData['pId'])+'::'
-        A = G2lat.cell2A(Cell[:6])
-        cellSig = getCellEsd(pfx,SGData,A,covData)
-        names = [' a = ',' b = ',' c = ',' alpha = ',' beta = ',' gamma = ',' Volume = ']
-        valEsd = [G2mth.ValEsd(Cell[i],cellSig[i],True) for i in range(7)]
-        line = '\n Unit cell:'
-        for name,vals in zip(names,valEsd):
-            line += name+vals  
-        print line
-    else: 
-        print '\n Unit cell: a = ','%.5f'%(Cell[0]),' b = ','%.5f'%(Cell[1]),' c = ','%.5f'%(Cell[2]), \
-            ' alpha = ','%.3f'%(Cell[3]),' beta = ','%.3f'%(Cell[4]),' gamma = ', \
-            '%.3f'%(Cell[5]),' volume = ','%.3f'%(Cell[6])
     #find one end of 4 atom string - involved in longest distance
     dist = {}
@@ -2717 +2704 @@
     dist = {}
     X1 = TorsionData['Datoms'][end]
-    for i,X2 in enumerate(TorsionData['Datoms']):                
+    for i,X2 in enumerate(TorsionData['Datoms']):
         dist[np.sqrt(np.sum(np.inner(Amat,np.array(X2[3:6])-np.array(X1[3:6]))**2))] = i
     sortdist = dist.keys()
     sortdist.sort()
     Datoms = []
+    Oatoms = []
     for d in sortdist:
         atom = TorsionData['Datoms'][dist[d]]
@@ -2730 +2718 @@
         symop[1] = eval(symop[1])
         atom[-1] = symop
-        print atom
         Datoms.append(atom)
-    Tors,sig = G2mth.GetTorsionSig(Datoms,Amat,SGData,covData={})
-    print ' Torsion: ',G2mth.ValEsd(Tors,sig)
-        
-        
+        oatom = TorsionData['Oatoms'][dist[d]]
+        names = ['','','']
+        if pfx:
+            names = [pfx+'dAx:'+str(oatom[0]),pfx+'dAy:'+str(oatom[0]),pfx+'dAz:'+str(oatom[0])]
+        oatom += [names,]
+        Oatoms.append(oatom)
+    Tors,sig = G2mth.GetTorsionSig(Oatoms,Datoms,Amat,SGData,covData)
+    print ' Torsion angle for atom sequence: ',[Datoms[i][1] for i in range(4)],'=',G2mth.ValEsd(Tors,sig)
+    print ' NB: Atom sequence determined by interatomic distances'
+                
 def BestPlane(PlaneData):