Changeset 882


Ignore:
Timestamp:
Apr 5, 2013 2:36:03 PM (9 years ago)
Author:
vondreele
Message:

more RB progress

Location:
trunk
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIImath.py

    r881 r882  
    257257            QuatA = AVdeg2Q(tor[0],Cart[seq[0]]-Cart[seq[1]])
    258258            for ride in seq[3]:
    259                 XYZ[ride] = prodQVQ(QuatA,Cart[ride]-Cart[seq[1]])+Cart[seq[1]]
     259                Cart[ride] = prodQVQ(QuatA,Cart[ride]-Cart[seq[1]])+Cart[seq[1]]
    260260    XYZ = np.zeros_like(Cart)
    261261    for i,xyz in enumerate(Cart):
    262         xyz = prodQVQ(RBObj['Orient'][0],xyz)
    263         XYZ[i] = np.inner(Bmat,xyz)+RBObj['Orig'][0]
     262        X = prodQVQ(RBObj['Orient'][0],xyz)
     263        XYZ[i] = np.inner(Bmat,X)+RBObj['Orig'][0]
    264264    return XYZ,Cart
    265265   
     
    268268    '''
    269269    TLStype,TLS = RBObj['ThermalMotion'][:2]
    270 #    T = TLS[:6]
    271 #    L = TLS[6:12]
    272 #    S = TLS[12:]
    273 #    Uout = []
    274 #    for X in Cart:
    275 #        U = [0,0,0,0,0,0]
    276 #        U[0] = T[0]+L[1]*X[2]**2+L[2]*X[1]**2-2.0*L[5]*X[1]*X[2]+2.0*(S[2]*X[2]-S[4]*X[1])
    277 #        U[1] = T[1]+L[0]*X[2]**2+L[2]*X[0]**2-2.0*L[4]*X[0]*X[2]+2.0*(S[5]*X[0]-S[0]*X[2])
    278 #        U[2] = T[2]+L[1]*X[0]**2+L[0]*X[1]**2-2.0*L[3]*X[1]*X[0]+2.0*(S[1]*X[1]-S[3]*X[0])
    279 #        U[3] = T[3]+L[4]*X[1]*X[2]+L[5]*X[0]*X[2]-L[3]*X[2]**2-L[2]*X[0]*X[1]+
    280 #            S[4]*X[0]-S[5]*X[1]-(S[6]+S[7])*X[2]
    281 #        U[4] = T[4]+L[3]*X[1]*X[2]+L[5]*X[0]*X[1]-L[4]*X[1]**2-L[1]*X[0]*X[2]+
    282 #            S[3]*X[2]-S[2]*X[0]+S[6]*X[1]
    283 #        U[5] = T[5]+L[3]*X[0]*X[2]+L[4]*X[0]*X[1]-L[5]*X[0]**2-L[0]*X[2]*X[1]+
    284 #            S[0]*X[1]-S[1]*X[2]+S[7]*X[0]
    285 #        UIJ = G2lat.U6toUij(U)
    286 #        Uout.append(
    287    
    288270    Tmat = np.zeros((3,3))
    289271    Lmat = np.zeros((3,3))
    290272    Smat = np.zeros((3,3))
    291273    if 'T' in TLStype:
    292         Tmat = G2lat.U6toUij(TLS[:6])
     274        T = TLS[:6]
     275#        Tmat = G2lat.U6toUij(T)
    293276    if 'L' in TLStype:
    294         Lmat = np.array(G2lat.U6toUij(TLS[6:12]))*(np.pi/180.)**2
     277        L = np.array(TLS[6:12])*(np.pi/180.)**2
     278#        Lmat = np.array(G2lat.U6toUij(L))
    295279    if 'S' in TLStype:
    296         Smat = np.array([[TLS[18],TLS[12],TLS[13]],[TLS[14],TLS[19],TLS[15]],[TLS[16],TLS[17],0]])*(np.pi/180.)
     280        S = np.array(TLS[12:])*(np.pi/180.)
     281#        Smat = np.array([[S[6],S[0],S[1]],[S[2],S[7],S[3]],[S[4],S[5],0]])
    297282    g = np.inner(Bmat,Bmat.T)
    298283    gvec = 1./np.sqrt(np.array([g[0][0]**2,g[1][1]**2,g[2][2]**2,
    299284        g[0][0]*g[1][1],g[0][0]*g[2][2],g[1][1]*g[2][2]]))
    300285    Uout = []
     286    Q = RBObj['Orient'][0]
    301287    for X in Cart:
    302         print 'X: ',X
    303         Axyz = np.array([[0,X[2],-X[1]], [-X[2],0,X[0]], [X[1],-X[0],0]])
     288        X = prodQVQ(Q,X)
     289#        Axyz = np.array([[0,X[2],-X[1]], [-X[2],0,X[0]], [X[1],-X[0],0]])
    304290        if 'U' in TLStype:
    305291            Uout.append(['I',TLS[0],0,0,0,0,0,0])
    306         else:
    307             Umat = Tmat+np.inner(Axyz,Smat)+np.inner(Smat.T,Axyz.T)+np.inner(np.inner(Axyz,Lmat),Axyz.T)
    308             print 'Umat: ',Umat
    309             beta = np.inner(np.inner(g,Umat),g)
     292        elif not 'N' in TLStype:
     293            U = [0,0,0,0,0,0]
     294            U[0] = T[0]+L[1]*X[2]**2+L[2]*X[1]**2-2.0*L[5]*X[1]*X[2]+2.0*(S[2]*X[2]-S[4]*X[1])
     295            U[1] = T[1]+L[0]*X[2]**2+L[2]*X[0]**2-2.0*L[4]*X[0]*X[2]+2.0*(S[5]*X[0]-S[0]*X[2])
     296            U[2] = T[2]+L[1]*X[0]**2+L[0]*X[1]**2-2.0*L[3]*X[1]*X[0]+2.0*(S[1]*X[1]-S[3]*X[0])
     297            U[3] = T[3]+L[4]*X[1]*X[2]+L[5]*X[0]*X[2]-L[3]*X[2]**2-L[2]*X[0]*X[1]+  \
     298                S[4]*X[0]-S[5]*X[1]-(S[6]+S[7])*X[2]
     299            U[4] = T[4]+L[3]*X[1]*X[2]+L[5]*X[0]*X[1]-L[4]*X[1]**2-L[1]*X[0]*X[2]+  \
     300                S[3]*X[2]-S[2]*X[0]+S[6]*X[1]
     301            U[5] = T[5]+L[3]*X[0]*X[2]+L[4]*X[0]*X[1]-L[5]*X[0]**2-L[0]*X[2]*X[1]+  \
     302                S[0]*X[1]-S[1]*X[2]+S[7]*X[0]
     303            Umat = G2lat.U6toUij(U)
     304#            print 'Umat: ',Umat
     305#wrong?      Umat1 = Tmat+np.inner(Axyz,Smat)+np.inner(Smat.T,Axyz.T)+np.inner(np.inner(Axyz.T,Lmat),Axyz)
     306#            print 'Umat1: ',Umat1
     307            beta = np.inner(np.inner(Bmat,Umat),Bmat.T)
    310308            Uout.append(['A',0.0,]+list(G2lat.UijtoU6(beta)*gvec))
    311309    return Uout
  • trunk/GSASIIphsGUI.py

    r881 r882  
    32973297            thermSel.Bind(wx.EVT_COMBOBOX,OnThermSel)
    32983298            thermSizer.Add(thermSel,0,wx.ALIGN_CENTER_VERTICAL)
     3299            thermSizer.Add(wx.StaticText(RigidBodies,-1,' Units: T A^2, L deg^2, S deg-A'),0,wx.ALIGN_CENTER_VERTICAL)
    32993300            resrbSizer.Add(thermSizer)
    33003301            if RBObj['ThermalMotion'][0] != 'None':
     
    33363337            thermSel.Bind(wx.EVT_COMBOBOX,OnThermSel)
    33373338            thermSizer.Add(thermSel,0,wx.ALIGN_CENTER_VERTICAL)
     3339            thermSizer.Add(wx.StaticText(RigidBodies,-1,' Units: T A^2, L deg^2, S deg-A'),0,wx.ALIGN_CENTER_VERTICAL)
    33383340            vecrbSizer.Add(thermSizer)
    33393341            if RBObj['ThermalMotion'][0] != 'None':
  • trunk/GSASIIstruct.py

    r881 r882  
    596596    current RB values in parmDict & modifies atom contents (xyz & Uij) of parmDict
    597597    '''
     598    atxIds = ['Ax:','Ay:','Az:']
     599    atuIds = ['AU11:','AU22:','AU33:','AU12:','AU13:','AU23:']
    598600    RBIds = rigidbodyDict['RBIds']
    599601    if not RBIds['Vector'] and not RBIds['Residue']:
    600602        return
     603    VRBIds = RBIds['Vector']
     604    RRBIds = RBIds['Residue']
    601605    RBData = copy.deepcopy(rigidbodyDict)     # don't mess with original!
    602606    if RBIds['Vector']:                       # first update the vector magnitudes
    603         VRBIds = RBIds['Vector']
    604607        VRBData = RBData['Vector']
    605608        for i,rbId in enumerate(VRBIds):
     
    613616        Amat,Bmat = G2lat.cell2AB(cell)
    614617        AtLookup = G2mth.FillAtomLookUp(Phase['Atoms'])
    615         pId = Phase['pId']
     618        pfx = str(Phase['pId'])+'::'
    616619        RBModels =  copy.deepcopy(Phase['RBModels']) # again don't mess with original!
    617         for irb,RBObj in enumerate(RBModels['Vector']):
    618             print irb,RBObj
    619            
     620        for irb,rbId in enumerate(VRBIds):
     621            RBObj = RBModels['Vector'][irb]
    620622            XYZ,Cart = G2mth.UpdateRBXYZ(Bmat,RBObj,RBData,'Vector')
    621             for x in XYZ: print x
    622623            UIJ = G2mth.UpdateRBUIJ(Bmat,Cart,RBObj)
    623             for u in UIJ: print u
    624            
    625         for irb,RBObj in enumerate(RBModels['Residue']):
    626             print irb,RBObj
    627            
    628            
    629        
    630             XYZ = G2mth.UpdateRBXYZ(Bmat,RBObj,RBData,'Residue')
    631             UIJ = G2mth.UpdateRBUIJ(Amat,XYZ,RBObj)
    632            
    633     raise Exception
    634 
     624            for i,x in enumerate(XYZ):
     625                atId = RBObj['Ids'][i]
     626                for j in [0,1,2]:
     627                    parmDict[pfx+atxIds[j]+str(AtLookup[atId])] = x[j]
     628                if UIJ[0] == 'A':
     629                    for j in range(6):
     630                        parmDict[pfx+atuIds[j]+str(AtLookup[atId])] = Uij[j+2]
     631                elif UIJ[0] == 'I':
     632                    parmDict[pfx+'AUiso:'+str(AtLookup[atId])] = Uij[j+1]
     633        for irb,rbId in enumerate(RRBIds):
     634            RBObj = RBModels['Residue'][irb]           
     635            XYZ,Cart = G2mth.UpdateRBXYZ(Bmat,RBObj,RBData,'Residue')
     636            UIJ = G2mth.UpdateRBUIJ(Bmat,Cart,RBObj)
     637            for i,x in enumerate(XYZ):
     638                atId = RBObj['Ids'][i]
     639                for j in [0,1,2]:
     640                    parmDict[pfx+atxIds[j]+str(AtLookup[atId])] = x[j]
     641                if UIJ[0] == 'A':
     642                    for j in range(6):
     643                        parmDict[pfx+atuIds[j]+str(AtLookup[atId])] = Uij[j+2]
     644                elif UIJ[0] == 'I':
     645                    parmDict[pfx+'AUiso:'+str(AtLookup[atId])] = Uij[j+1]
    635646       
    636647################################################################################
     
    33143325        GA,GB = G2lat.Gmat2AB(G)    #Orthogonalization matricies
    33153326        Vst = np.sqrt(nl.det(G))    #V*
    3316 #apply RB parms to atom parms in parmDict?
    33173327        if not Phase['General'].get('doPawley'):
    33183328            refList = StructureFactor(refList,G,hfx,pfx,SGData,calcControls,parmDict)
Note: See TracChangeset for help on using the changeset viewer.