Changeset 458


Ignore:
Timestamp:
Jan 27, 2012 4:43:30 PM (10 years ago)
Author:
vondreele
Message:

finish dist angle calculations

Location:
trunk
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIImath.py

    r457 r458  
    1616import math
    1717import GSASIIpath
     18import GSASIIspc as G2spc
    1819import scipy.optimize as so
    1920import scipy.linalg as sl
     
    2324tand = lambda x: np.tan(x*np.pi/180.)
    2425asind = lambda x: 180.*np.arcsin(x)/np.pi
     26acosd = lambda x: 180.*np.arccos(x)/np.pi
    2527atan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
    2628
     
    143145    return vcov
    144146   
     147def 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,(TxT-Ox))**2))
     153       
     154    inv = 1
     155    if Top < 0:
     156        inv = -1
     157    cent = abs(Top)/100
     158    op = abs(Top)%100-1
     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   
     176def 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,(TxT-Ox))
     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 = VecB-VecA
     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)%100-1
     204    opB = abs(TopB)%100-1
     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   
    145243def ValEsd(value,esd=0,nTZ=False):                  #NOT complete - don't use
    146244    # returns value(esd) string; nTZ=True for no trailing zeros
  • trunk/GSASIIphsGUI.py

    r457 r458  
    192192            pass
    193193        Obj.SetValue("%.3f"%(self.data[item[0]][item[1]]))          #reset in case of error
     194       
     195    def getData(self):
     196        return self.Data
    194197       
    195198    def OnOk(self,event):
     
    10591062            cn = colLabels.index('Name')
    10601063            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])
    10621065                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])
    10641067            DisAglData['OrigAtoms'] = Oxyz
    10651068            DisAglData['TargAtoms'] = xyz
  • trunk/GSASIIstruct.py

    r457 r458  
    2626tand = lambda x: np.tan(x*np.pi/180.)
    2727asind = lambda x: 180.*np.arcsin(x)/np.pi
     28acosd = lambda x: 180.*np.arccos(x)/np.pi
    2829atan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
    2930
     
    25692570   
    25702571    Amat,Bmat = G2lat.cell2AB(Cell[:6])
     2572    covData = {}
    25712573    if 'covData' in DisAglData:   
    25722574        covData = DisAglData['covData']
     2575        covMatrix = covData['covMatrix']
     2576        varyList = covData['varyList']
    25732577        pfx = str(DisAglData['pId'])+'::'
    25742578        A = G2lat.cell2A(Cell[:6])
     
    25802584            line += name+vals 
    25812585        print line
    2582     else:
     2586    else: 
    25832587        print '\n Unit cell: a = ','%.5f'%(Cell[0]),' b = ','%.5f'%(Cell[1]),' c = ','%.5f'%(Cell[2]), \
    25842588            ' alpha = ','%.3f'%(Cell[3]),' beta = ','%.3f'%(Cell[4]),' gamma = ', \
     
    25902594        [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],
    25912595        [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']
    25932597    targAtoms = DisAglData['TargAtoms']
    25942598    for Oatom in origAtoms:
     2599        OxyzNames = ''
     2600        IndBlist = []
    25952601        Dist = []
    25962602        Vect = []
     2603        VectA = []
     2604        angles = []
    25972605        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]
    26012615            for Txyz,Top,Tunit in result:
    2602                 Dx = (Txyz-np.array(Oatom[1][2:5]))+Units
     2616                Dx = (Txyz-np.array(Oatom[3:6]))+Units
    26032617                dx = np.inner(Amat,Dx)
    26042618                dist = ma.masked_less(np.sqrt(np.sum(dx**2,axis=0)),0.5)
    26052619                IndB = ma.nonzero(ma.masked_greater(dist-BsumR,0.))
    2606                 IndA = ma.nonzero(ma.masked_greater(dist-AsumR,0.))
    26072620                if np.any(IndB):
    26082621                    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       
    26142676           
    2615         for vect in Vect:
    2616             print vect
    2617            
    2618                  
    2619 
    2620 #    def FindBonds():                    #uses numpy & masks - very fast even for proteins!
    2621 #        import numpy.ma as ma
    2622 #        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.5
    2630 #        except:
    2631 #            pass           
    2632 #        for atom in atomData:
    2633 #            atom[-2] = []               #clear out old bond/angle tables
    2634 #            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 = Atoms-atomA[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.5A
    2653 #                sumR = atomA[3]+Radii
    2654 #                IndB = ma.nonzero(ma.masked_greater(dist-radiusFactor*sumR,0.))                 #get indices of bonded atoms
    2655 #                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 
    26622677def main():
    26632678    arg = sys.argv
Note: See TracChangeset for help on using the changeset viewer.