Changeset 458 for trunk/GSASIIstruct.py


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

finish dist angle calculations

File:
1 edited

Legend:

Unmodified
Added
Removed
  • 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.