Changeset 477


Ignore:
Timestamp:
Feb 9, 2012 1:23:01 PM (10 years ago)
Author:
vondreele
Message:

final fixes to dist-angle-torsion calc.

Location:
trunk
Files:
5 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIImath.py

    r470 r477  
    148148   
    149149    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
     150        TxT = inv*(np.inner(M,Tx)+T)+C+U
    152151        return np.sqrt(np.sum(np.inner(Amat,(TxT-Ox))**2))
    153152       
     
    236235        return calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat),0.0
    237236       
    238 def GetTorsionSig(Atoms,Amat,SGData,covData={}):
    239     XYZ = []
    240     for atom in Atoms:   
    241         XYZ.append(np.array(atom[3:6]))
    242     XYZ = np.inner(Amat,XYZ).T
    243     V1 = XYZ[1]-XYZ[0]
    244     V2 = XYZ[2]-XYZ[1]
    245     V3 = XYZ[3]-XYZ[2]
    246     V1 /= np.sqrt(np.sum(V1**2))
    247     V2 /= np.sqrt(np.sum(V2**2))
    248     V3 /= np.sqrt(np.sum(V3**2))
    249     M = np.array([V1,V2,V3])
    250     D = nl.det(M)
    251     Ang = 1.0
    252     P12 = np.dot(V1,V2)
    253     P13 = np.dot(V1,V3)
    254     P23 = np.dot(V2,V3)
    255     Tors = acosd((P12*P23-P13)/(np.sqrt(1.-P12**2)*np.sqrt(1.-P23**2)))*D/abs(D)
    256    
    257     sig = 0.0
     237def GetTorsionSig(Oatoms,Atoms,Amat,SGData,covData={}):
     238   
     239    def calcTorsion(Atoms,SyOps,Amat):
     240       
     241        XYZ = []
     242        for i,atom in enumerate(Atoms):
     243            Inv,M,T,C,U = SyOps[i]
     244            XYZ.append(np.array(atom[1:4]))
     245            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
     246            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
     247        V1 = XYZ[1]-XYZ[0]
     248        V2 = XYZ[2]-XYZ[1]
     249        V3 = XYZ[3]-XYZ[2]
     250        V1 /= np.sqrt(np.sum(V1**2))
     251        V2 /= np.sqrt(np.sum(V2**2))
     252        V3 /= np.sqrt(np.sum(V3**2))
     253        M = np.array([V1,V2,V3])
     254        D = nl.det(M)
     255        Ang = 1.0
     256        P12 = np.dot(V1,V2)
     257        P13 = np.dot(V1,V3)
     258        P23 = np.dot(V2,V3)
     259        Tors = acosd((P12*P23-P13)/(np.sqrt(1.-P12**2)*np.sqrt(1.-P23**2)))*D/abs(D)
     260        return Tors
     261           
     262    Inv = []
     263    SyOps = []
     264    names = []
     265    for i,atom in enumerate(Oatoms):
     266        names += atom[-1]
     267        Op,unit = Atoms[i][-1]
     268        inv = Op/abs(Op)
     269        m,t = SGData['SGOps'][abs(Op)%100-1]
     270        c = SGData['SGCen'][abs(Op)/100]
     271        SyOps.append([inv,m,t,c,unit])
     272    Tors = calcTorsion(Oatoms,SyOps,Amat)
     273   
     274    sig = -0.01
    258275    if 'covMatrix' in covData:
    259         for atom in Atoms:
    260             xyz = atom[3:6]
    261             Op,Unit = atom[-1]
    262             inv = Op/abs(Op)
    263             M,T = SGData['SGOps'][abs(Op)%100-1]
    264             C = SGData['SGCen'][abs(Op)/100]
    265             #reverse inv*(np.inner(M,Tx)+T)+C
    266             XYZ = np.inner(nl.inv(M),((xyz-C)*inv-T))
    267             print Op,Unit,xyz,XYZ
     276        parmNames = []
     277        dx = .00001
     278        dadx = np.zeros(12)
     279        ang = calcTorsion(Oatoms,SyOps,Amat)
     280        for i in range(12):
     281            ia = i/3
     282            ix = i%3
     283            Oatoms[ia][ix+1] += dx
     284            a0 = calcTorsion(Oatoms,SyOps,Amat)
     285            Oatoms[ia][ix+1] -= 2*dx
     286            dadx[i] = (calcTorsion(Oatoms,SyOps,Amat)-a0)/(2.*dx)
    268287        covMatrix = covData['covMatrix']
    269288        varyList = covData['varyList']
     289        TorVcov = getVCov(names,varyList,covMatrix)
     290        sig = np.sqrt(np.inner(dadx,np.inner(TorVcov,dadx)))
     291        if sig < 0.01:
     292            sig = -0.01
    270293   
    271294    return Tors,sig
     
    286309        return str(fmt%(value,int(round(esd*10**(mdec(esd)))))).strip('"')
    287310    elif esd < 0:
    288          return str(round(value,mdec(esd)))
     311         return str(round(value,mdec(esd)-1))
    289312    else:
    290313        text = str("%f"%(value))
  • trunk/GSASIIphsGUI.py

    r476 r477  
    16481648                if atom[cuia] == 'A':
    16491649                    Uij = atom[cuij:cuij+6]
    1650                     result = G2spc.GenAtom(XYZ,SGData,False,Uij)
     1650                    result = G2spc.GenAtom(XYZ,SGData,False,Uij,False)
    16511651                    for item in result:
    16521652                        atom = copy.copy(atomData[ind])
     
    16651665                                atomData.append(atom[:])
    16661666                else:
    1667                     result = G2spc.GenAtom(XYZ,SGData,False)
     1667                    result = G2spc.GenAtom(XYZ,SGData,False,Move=False)
    16681668                    for item in result:
    16691669                        atom = copy.copy(atomData[ind])
     
    18021802        event.StopPropagation()
    18031803       
    1804     def FindAtomIndexByIDs(atomData,IDs):
     1804    def FindAtomIndexByIDs(atomData,IDs,Draw=True):
    18051805        indx = []
    18061806        for i,atom in enumerate(atomData):
    1807             if atom[-2] in IDs:
     1807            if Draw and atom[-3] in IDs:
     1808                indx.append(i)
     1809            elif atom[-1] in IDs:
    18081810                indx.append(i)
    18091811        return indx
     
    18351837            return
    18361838        TorsionData = {}
     1839        ocx,oct,ocs,cia = data['General']['AtomPtrs']
    18371840        drawingData = data['Drawing']
    18381841        atomData = data['Atoms']
     1842#        for atom in atomData: print atom
    18391843        atomDData = drawingData['Atoms']
     1844#        for atom in atomDData: print atom
    18401845        colLabels = [drawAtoms.GetColLabelValue(c) for c in range(drawAtoms.GetNumberCols())]
    18411846        cx = colLabels.index('x')
    18421847        cn = colLabels.index('Name')
     1848        cid = colLabels.index('I/A')+8
    18431849        xyz = []
     1850        Oxyz = []
    18441851        for i,atom in enumerate(atomDData):
    18451852            if i in indx:
    18461853                xyz.append([i,]+atom[cn:cn+2]+atom[cx:cx+4]) #also gets Sym Op
     1854                id = FindAtomIndexByIDs(atomData,[atom[cid],],False)[0]
     1855                Oxyz.append([id,]+atomData[id][cx+1:cx+4])
    18471856        TorsionData['Datoms'] = xyz
     1857        TorsionData['Oatoms'] = Oxyz
    18481858        generalData = data['General']
    18491859        TorsionData['Name'] = generalData['Name']
  • trunk/GSASIIplot.py

    r468 r477  
    22662266            return G2frame.dataDisplay.GetSelection()
    22672267        except AttributeError:
    2268             print G2frame.dataDisplay.GetLabel()
    22692268            G2frame.G2plotNB.status.SetStatusText('Select this from Phase data window!')
    22702269            return 0
     
    25722571                           
    25732572    def Draw():
    2574         import numpy.linalg as nl
    25752573        Ind = GetSelectedAtoms()
    25762574        VS = np.array(Page.canvas.GetSize())
  • trunk/GSASIIspc.py

    r459 r477  
    252252    return new
    253253       
    254 def GenAtom(XYZ,SGData,All=False,Uij=[]):
     254def GenAtom(XYZ,SGData,All=False,Uij=[],Move=True):
    255255    '''
    256256    Generates the equivalent positions for a specified coordinate and space group
     
    273273    Cell = []
    274274    X = np.array(XYZ)
    275     X = MoveToUnitCell(X)
     275    if Move:
     276        X = MoveToUnitCell(X)
    276277    for ic,cen in enumerate(SGData['SGCen']):
    277278        C = np.array(cen)
  • trunk/GSASIIstruct.py

    r474 r477  
    26102610                TxyzNames = [pfx+'dAx:%d'%(Tatom[0]),pfx+'dAy:%d'%(Tatom[0]),pfx+'dAz:%d'%(Tatom[0])]
    26112611                Xvcov = G2mth.getVCov(OxyzNames+TxyzNames,varyList,covMatrix)
    2612             result = G2spc.GenAtom(Tatom[3:6],SGData,False)
     2612            result = G2spc.GenAtom(Tatom[3:6],SGData,False,Move=False)
    26132613            BsumR = (Radii[Oatom[2]][0]+Radii[Tatom[2]][0])*Factor[0]
    26142614            AsumR = (Radii[Oatom[2]][1]+Radii[Tatom[2]][1])*Factor[1]
     
    26792679        print 80*'*'
    26802680        print '   Torsion angle for phase '+name
    2681         print 80*'*','\n'
     2681        print 80*'*'
    26822682
    26832683    ShowBanner(TorsionData['Name'])
    26842684    SGData = TorsionData['SGData']
    2685     SGtext = G2spc.SGPrint(SGData)
    2686     for line in SGtext: print line
    26872685    Cell = TorsionData['Cell']
    26882686   
    26892687    Amat,Bmat = G2lat.cell2AB(Cell[:6])
    26902688    covData = {}
     2689    pfx = ''
    26912690    if 'covData' in TorsionData:   
    26922691        covData = TorsionData['covData']
     
    26942693        varyList = covData['varyList']
    26952694        pfx = str(TorsionData['pId'])+'::'
    2696         A = G2lat.cell2A(Cell[:6])
    2697         cellSig = getCellEsd(pfx,SGData,A,covData)
    2698         names = [' a = ',' b = ',' c = ',' alpha = ',' beta = ',' gamma = ',' Volume = ']
    2699         valEsd = [G2mth.ValEsd(Cell[i],cellSig[i],True) for i in range(7)]
    2700         line = '\n Unit cell:'
    2701         for name,vals in zip(names,valEsd):
    2702             line += name+vals 
    2703         print line
    2704     else:
    2705         print '\n Unit cell: a = ','%.5f'%(Cell[0]),' b = ','%.5f'%(Cell[1]),' c = ','%.5f'%(Cell[2]), \
    2706             ' alpha = ','%.3f'%(Cell[3]),' beta = ','%.3f'%(Cell[4]),' gamma = ', \
    2707             '%.3f'%(Cell[5]),' volume = ','%.3f'%(Cell[6])
    27082695    #find one end of 4 atom string - involved in longest distance
    27092696    dist = {}
     
    27172704    dist = {}
    27182705    X1 = TorsionData['Datoms'][end]
    2719     for i,X2 in enumerate(TorsionData['Datoms']):               
     2706    for i,X2 in enumerate(TorsionData['Datoms']):
    27202707        dist[np.sqrt(np.sum(np.inner(Amat,np.array(X2[3:6])-np.array(X1[3:6]))**2))] = i
    27212708    sortdist = dist.keys()
    27222709    sortdist.sort()
    27232710    Datoms = []
     2711    Oatoms = []
    27242712    for d in sortdist:
    27252713        atom = TorsionData['Datoms'][dist[d]]
     
    27302718        symop[1] = eval(symop[1])
    27312719        atom[-1] = symop
    2732         print atom
    27332720        Datoms.append(atom)
    2734     Tors,sig = G2mth.GetTorsionSig(Datoms,Amat,SGData,covData={})
    2735     print ' Torsion: ',G2mth.ValEsd(Tors,sig)
    2736        
    2737        
     2721        oatom = TorsionData['Oatoms'][dist[d]]
     2722        names = ['','','']
     2723        if pfx:
     2724            names = [pfx+'dAx:'+str(oatom[0]),pfx+'dAy:'+str(oatom[0]),pfx+'dAz:'+str(oatom[0])]
     2725        oatom += [names,]
     2726        Oatoms.append(oatom)
     2727    Tors,sig = G2mth.GetTorsionSig(Oatoms,Datoms,Amat,SGData,covData)
     2728    print ' Torsion angle for atom sequence: ',[Datoms[i][1] for i in range(4)],'=',G2mth.ValEsd(Tors,sig)
     2729    print ' NB: Atom sequence determined by interatomic distances'
     2730               
    27382731def BestPlane(PlaneData):
    27392732
Note: See TracChangeset for help on using the changeset viewer.