Changeset 4595 for trunk


Ignore:
Timestamp:
Oct 15, 2020 5:43:24 PM (14 months ago)
Author:
toby
Message:

add RB sym axis; implement in RB add; minor RB GUI reorg

Location:
trunk
Files:
4 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIImath.py

    r4583 r4595  
    979979    return np.array([X,Mat3,Mat2])       
    980980               
    981 def RotateRBXYZ(Bmat,Cart,oriQ):
     981def RotateRBXYZ(Bmat,Cart,oriQ,symAxis=None):
    982982    '''rotate & transform cartesian coordinates to crystallographic ones
    983983    no translation applied. To be used for numerical derivatives
     
    986986    :param array Cart: 2D array of coordinates
    987987    :param array Q: quaternion as an np.array
    988    
     988    :param tuple symAxis: if not None (default), specifies the symmetry
     989      axis of the rigid body, which will be aligned to the quaternion vector.
    989990    :returns: 2D array of fractional coordinates, without translation to origin
    990991    '''
     992    if symAxis is None:
     993        Q = oriQ
     994    else:
     995        a,v = Q2AV(oriQ)
     996        symaxis = np.array(symAxis)
     997        xformAng = np.arccos(np.dot(v,symaxis))
     998        xformVec = np.cross(symaxis,v)
     999        Q = prodQQ(oriQ,AV2Q(xformAng,xformVec))
    9911000    XYZ = np.zeros_like(Cart)
    9921001    for i,xyz in enumerate(Cart):
    993         XYZ[i] = np.inner(Bmat,prodQVQ(oriQ,xyz))
     1002        XYZ[i] = np.inner(Bmat,prodQVQ(Q,xyz))
    9941003    return XYZ
    9951004
    9961005def UpdateRBXYZ(Bmat,RBObj,RBData,RBType):
    997     '''returns crystal coordinates for atoms described by RBObj
     1006    '''returns crystal coordinates for atoms described by RBObj.
     1007    Note that RBObj['symAxis'], if present, determines the symmetry
     1008    axis of the rigid body, which will be aligned to the
     1009    quaternion direction.
    9981010   
    9991011    :param np.array Bmat: see :func:`GSASIIlattice.cell2AB`
     
    10171029            QuatA = AVdeg2Q(tor[0],Cart[seq[0]]-Cart[seq[1]])
    10181030            Cart[seq[3]] = prodQVQ(QuatA,(Cart[seq[3]]-Cart[seq[1]]))+Cart[seq[1]]
     1031    # if symmetry axis is defined, place symmetry axis along quaternion
     1032    if RBObj.get('symAxis') is None:
     1033        Q = RBObj['Orient'][0]
     1034    else:
     1035        a,v = Q2AV(RBObj['Orient'][0])
     1036        symaxis = np.array(RBObj.get('symAxis'))
     1037        xformAng = np.arccos(np.dot(v,symaxis))
     1038        xformVec = np.cross(symaxis,v)
     1039        Q = prodQQ(RBObj['Orient'][0],AV2Q(xformAng,xformVec))
    10191040    XYZ = np.zeros_like(Cart)
    10201041    for i,xyz in enumerate(Cart):
    1021         XYZ[i] = np.inner(Bmat,prodQVQ(RBObj['Orient'][0],xyz))+RBObj['Orig'][0]
     1042        XYZ[i] = np.inner(Bmat,prodQVQ(Q,xyz))+RBObj['Orig'][0]
    10221043    return XYZ,Cart
    10231044
  • trunk/GSASIIphsGUI.py

    r4594 r4595  
    15651565    '''Update all orientation text on the Add RB panel
    15661566    '''
    1567     for i,sizer in enumerate(G2frame.testRBObjSizers['Osizers']):
    1568         sizer.SetLabel('%8.5f'%(testRBObj['rbObj']['Orient'][0][i]))
    15691567    A,V = G2mth.Q2AVdeg(testRBObj['rbObj']['Orient'][0])
    15701568    testRBObj['rbObj']['OrientVec'][0] = A
     
    99349932                G2plt.PlotStructure(G2frame,data)
    99359933               
    9936             def OnOrien(event):
    9937                 event.Skip()
    9938                 Obj = event.GetEventObject()
    9939                 item = Indx[Obj.GetId()]
    9940                 A,V = G2mth.Q2AVdeg(RBObj['Orient'][0])
    9941                 V = np.inner(Bmat,V)
     9934            def OnOrien(*args, **kwargs):
     9935                '''Called when the orientation info is changed (vector
     9936                or azimuth)
     9937                '''
    99429938                try:
    9943                     val = float(Obj.GetValue())
    9944                     if item:
    9945                         V[item-1] = val
    9946                     else:
    9947                         A = val
    9948                     Obj.SetValue('%8.5f'%(val))
    9949                     V = np.inner(Amat,V)
     9939                    orient = [float(Indx['Orien'][i].GetValue()) for i in range(4)]
     9940                    A = orient[0]
     9941                    V = np.inner(Amat,orient[1:]) # normalized in AVdeg2Q
    99509942                    Q = G2mth.AVdeg2Q(A,V)
    99519943                    if not any(Q):
     
    99659957                if GSASIIpath.GetConfigValue('debug'): print('patching origin!')
    99669958                RBObj['Orig'][0] = list(RBObj['Orig'][0])
    9967             Orien,OrienV = G2mth.Q2AVdeg(RBObj['Orient'][0])
    9968             Orien = [Orien,]
    9969             Orien.extend(OrienV/nl.norm(OrienV))
    9970             topSizer.Add(wx.StaticText(RigidBodies,-1,'Origin x,y,z:'),0,WACV)
     9959            topSizer.Add(wx.StaticText(RigidBodies,-1,'Origin x,y,z (frac)'),0,WACV)
    99719960            for ix in range(3):
    99729961                origX = G2G.ValidatedTxtCtrl(RigidBodies,RBObj['Orig'][0],ix,nDig=(8,5),
    9973                                     typeHint=float,OnLeave=OnOrigX)
     9962                                    typeHint=float,OnLeave=OnOrigX,
     9963                                    xmin=-1,xmax=1.,size=(70,-1))
    99749964                topSizer.Add(origX,0,WACV)
    99759965            topSizer.Add((5,0),)
     
    99789968            Ocheck.SetValue(RBObj['Orig'][1])
    99799969            topSizer.Add(Ocheck,0,WACV)
    9980             topSizer.Add(wx.StaticText(RigidBodies,-1,'Rotation angle, vector:'),0,WACV)
     9970            topSizer.Add(wx.StaticText(RigidBodies,-1,
     9971                    'Rotation angle (deg)\n&& Orient. vector (frac)'),0,WACV)
     9972            Indx['Orien'] = {}
     9973            Orien,OrienV = G2mth.Q2AVdeg(RBObj['Orient'][0])
     9974            Orien = [Orien,]
     9975            Orien.extend(np.inner(Bmat,OrienV)) # fractional coords
     9976            dp,xmin,xmax = 2,-180.,360.
    99819977            for ix,x in enumerate(Orien):
    9982 #            Zstep = G2G.ValidatedTxtCtrl(drawOptions,drawingData,'Zstep',nDig=(10,2),xmin=0.01,xmax=4.0)
    9983                 orien = wx.TextCtrl(RigidBodies,-1,value='%8.4f'%(x),style=wx.TE_PROCESS_ENTER)
    9984                 orien.Bind(wx.EVT_TEXT_ENTER,OnOrien)
    9985                 orien.Bind(wx.EVT_KILL_FOCUS,OnOrien)
    9986                 Indx[orien.GetId()] = ix
     9978                orien = G2G.ValidatedTxtCtrl(RigidBodies,Orien,ix,nDig=(8,dp),
     9979                                    typeHint=float,OnLeave=OnOrien,
     9980                                    xmin=xmin,xmax=xmax,size=(70,-1))
     9981                dp, xmin,xmax = 4,-1.,1.
     9982                Indx['Orien'][ix] = orien
    99879983                topSizer.Add(orien,0,WACV)
    99889984            Qcheck = wx.ComboBox(RigidBodies,-1,value='',choices=[' ','A','AV'],
     
    99909986            Qcheck.Bind(wx.EVT_COMBOBOX,OnOrienRef)
    99919987            Qcheck.SetValue(RBObj['Orient'][1])
    9992             topSizer.Add(Qcheck)
     9988            topSizer.Add(Qcheck,0,WACV)
    99939989            return topSizer
    99949990                         
     
    1003010026            Indx[delRB.GetId()] = rbId
    1003110027            topLine.Add(delRB,0,WACV)
     10028            symAxis = RBObj.get('symAxis')
     10029            if symAxis:
     10030                if symAxis[0]:
     10031                    lbl = 'x'
     10032                elif symAxis[1]:
     10033                    lbl = 'y'
     10034                else:
     10035                    lbl = 'z'                   
     10036                topLine.Add(wx.StaticText(RigidBodies,-1,
     10037                    '   Rigid body {} axis is aligned along oriention vector'
     10038                        .format(lbl)),0,WACV)
    1003210039            resrbSizer.Add(topLine)
    1003310040            resrbSizer.Add(LocationSizer(RBObj,'Residue'))
    10034             resrbSizer.Add(wx.StaticText(RigidBodies,-1,'Torsions:'),0)
     10041            if len(RBObj['Torsions']):
     10042                resrbSizer.Add(wx.StaticText(RigidBodies,-1,'Torsions:'),0)
    1003510043            torSizer = wx.FlexGridSizer(0,6,5,5)
    1003610044            for itors,tors in enumerate(RBObj['Torsions']):
     
    1005910067           
    1006010068        def VecrbSizer(RBObj):
    10061             G2frame.GetStatusBar().SetStatusText('NB: Rotation vector is in crystallographic space',1)
    10062                    
    1006310069            def OnDelVecRB(event):
    1006410070                Obj = event.GetEventObject()
     
    1016710173            G2frame.bottomSizer.Add(ResrbSizer(rbObj))
    1016810174            mainSizer.Add(G2frame.bottomSizer)
    10169             G2frame.GetStatusBar().SetStatusText('NB: Rotation vector is in crystallographic space',1)
    1017010175            G2plt.PlotStructure(G2frame,data)
     10176            G2plt.PlotStructure(G2frame,data) # draw twice initially for mac
    1017110177        if 'Vector' in data['RBModels'] and len(data['RBModels']['Vector']):
    1017210178            nobody = False
     
    1041810424                A,V = G2mth.Q2AVdeg(Q)
    1041910425                rbObj['OrientVec'][1:] = np.inner(Bmat,V)
    10420                 for i,sizer in enumerate(G2frame.testRBObjSizers['Osizers']):
    10421                     sizer.SetLabel('%8.5f'%(rbObj['Orient'][0][i]))
    1042210426                G2plt.PlotStructure(G2frame,data,False,UpdateTable)
    1042310427                UpdateTable()
     
    1043010434                    np.inner(Amat,rbObj['OrientVec'][1:]))
    1043110435                rbObj['Orient'][0] = Q
    10432                 for i,sizer in enumerate(G2frame.testRBObjSizers['Osizers']):
    10433                     sizer.SetLabel('%8.5f'%(rbObj['Orient'][0][i]))
    1043410436                G2frame.testRBObjSizers['OrientVecSiz'][4].SetValue(
    1043510437                    int(10*rbObj['OrientVec'][0]))               
     
    1058110583                G2plt.PlotStructure(G2frame,data,False,UpdateTable)
    1058210584                data['testRBObj']['showSelect'].SetLabelText(cryatom)
     10585            def OnSymRadioSet(event):
     10586                '''Set the symmetry axis for the body as
     10587                data['testRBObj']['rbObj']['symAxis']. This may never be
     10588                set, so use data['testRBObj']['rbObj'].get('symAxis') to
     10589                access this so the default value is None.
     10590                '''
     10591                data['testRBObj']['rbObj']['symAxis'] = (None,[1,0,0],[0,1,0],[0,0,1])[event.GetEventObject().GetSelection()]
     10592                UpdateTablePlot()               
    1058310593            showAtom = [None]
    1058410594            def showCryAtom(*args,**kwargs):
     
    1060110611            mainSizer = wx.BoxSizer(wx.VERTICAL)
    1060210612            mainSizer.Add((5,5),0)
    10603             Osizers = []
    1060410613            rbObj = data['testRBObj']['rbObj']
    1060510614            rbName = rbObj['RBname']
     
    1065010659                choice.sort()
    1065110660                G2frame.testRBObjSizers.update({'Xsizers':Xsizers})
    10652             OriSizer.Add(wx.StaticText(RigidBodies,wx.ID_ANY,'Orientation quaternion: '),0,WACV)
    10653             for ix,x in enumerate(Orien):
    10654                 orien = wx.StaticText(RigidBodies,wx.ID_ANY,'%10.4f'%(x))
    10655                 OriSizer.Add(orien,0,WACV)
    10656                 Osizers.append(orien)
    10657             G2frame.testRBObjSizers.update({'Osizers':Osizers})
    1065810661            mainSizer.Add(OriSizer)
    1065910662            mainSizer.Add((5,5),0)
     
    1068310686            OriSizer.Add(wx.StaticText(RigidBodies,wx.ID_ANY,' (frac coords)'),0,WACV)
    1068410687            G2frame.testRBObjSizers.update({'OrientVecSiz':OrientVecSiz})
     10688            mainSizer.Add(OriSizer)
     10689            mainSizer.Add((5,5),0)
     10690            OriSizer = wx.BoxSizer(wx.HORIZONTAL)
     10691            OriSizer.Add(wx.StaticText(RigidBodies,wx.ID_ANY,
     10692                                    'Rigid body symmetry axis: '),0, WACV)
     10693            choices = ['None','x','y','z']
     10694            symRadioSet = wx.RadioBox(RigidBodies,wx.ID_ANY,choices=choices)
     10695            symRadioSet.Bind(wx.EVT_RADIOBOX, OnSymRadioSet)
     10696            OriSizer.Add(symRadioSet)
    1068510697            mainSizer.Add(OriSizer)
    1068610698            mainSizer.Add((5,5),0)
  • trunk/GSASIIplot.py

    r4593 r4595  
    87278727        GL.glLightfv(GL.GL_LIGHT0,GL.GL_AMBIENT,[.2,.2,.2,1])
    87288728
    8729     def RenderRBtriplet(orig,Q,Bmat):
     8729    def RenderRBtriplet(orig,Q,Bmat,symAxis=None):
    87308730        '''draw an axes triplet located at the origin of a rigid body
    87318731        and with the x, y & z axes drawn as red, green and blue.
     
    87408740        GL.glTranslate(*orig)
    87418741        GL.glBegin(GL.GL_LINES)
    8742         lines = G2mth.RotateRBXYZ(Bmat,np.eye(3),Q)
     8742        lines = G2mth.RotateRBXYZ(Bmat,np.eye(3),Q,symAxis)
    87438743        colors = [Rd,Gr,Bl]
    87448744        # lines along axial directions
     
    93639363            RenderRBtriplet(testRBObj['rbObj']['Orig'][0],
    93649364                            testRBObj['rbObj']['Orient'][0],
    9365                             Bmat)
     9365                            Bmat,testRBObj['rbObj'].get('symAxis'))
    93669366        if len(mcsaModels) > 1 and pageName == 'MC/SA':             #skip the default MD entry
    93679367            for ind,[x,y,z] in enumerate(mcsaXYZ):
  • trunk/GSASIIstrMath.py

    r4557 r4595  
    151151                   
    152152def ApplyRBModelDervs(dFdvDict,parmDict,rigidbodyDict,Phase):
    153     'Needs a doc string'
     153    'Computes rigid body derivatives'
    154154    atxIds = ['dAx:','dAy:','dAz:']
    155155    atuIds = ['AU11:','AU22:','AU33:','AU12:','AU13:','AU23:']
     
    177177    RBModels =  Phase['RBModels']
    178178    for irb,RBObj in enumerate(RBModels.get('Vector',[])):
     179        symAxis = RBObj.get('symAxis')
    179180        VModel = RBData['Vector'][RBObj['RBId']]
    180181        Q = RBObj['Orient'][0]
     
    198199            for iv in range(4):
    199200                Q[iv] -= dx
    200                 XYZ1 = G2mth.RotateRBXYZ(Bmat,Cart,G2mth.normQ(Q))
     201                XYZ1 = G2mth.RotateRBXYZ(Bmat,Cart,G2mth.normQ(Q),symAxis)
    201202                Q[iv] += 2.*dx
    202                 XYZ2 = G2mth.RotateRBXYZ(Bmat,Cart,G2mth.normQ(Q))
     203                XYZ2 = G2mth.RotateRBXYZ(Bmat,Cart,G2mth.normQ(Q),symAxis)
    203204                Q[iv] -= dx
    204205                dXdO = (XYZ2[ia]-XYZ1[ia])/(2.*dx)
     
    238239
    239240    for irb,RBObj in enumerate(RBModels.get('Residue',[])):
     241        symAxis = RBObj.get('symAxis')
    240242        Q = RBObj['Orient'][0]
    241243        jrb = RRBIds.index(RBObj['RBId'])
     
    263265            for iv in range(4):
    264266                Q[iv] -= dx
    265                 XYZ1 = G2mth.RotateRBXYZ(Bmat,Cart,G2mth.normQ(Q))
     267                XYZ1 = G2mth.RotateRBXYZ(Bmat,Cart,G2mth.normQ(Q),symAxis)
    266268                Q[iv] += 2.*dx
    267                 XYZ2 = G2mth.RotateRBXYZ(Bmat,Cart,G2mth.normQ(Q))
     269                XYZ2 = G2mth.RotateRBXYZ(Bmat,Cart,G2mth.normQ(Q),symAxis)
    268270                Q[iv] -= dx
    269271                dXdO = (XYZ2[ia]-XYZ1[ia])/(2.*dx)
Note: See TracChangeset for help on using the changeset viewer.