Changeset 4595

Timestamp:

Oct 15, 2020 5:43:24 PM (15 months ago)

Author:

toby

Message:

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

Location:

trunk

Files:

• GSASIImath.py (modified) (3 diffs)
• GSASIIphsGUI.py (modified) (14 diffs)
• GSASIIplot.py (modified) (3 diffs)
• GSASIIstrMath.py (modified) (5 diffs)

trunk/GSASIImath.py

@@ -979 +979 @@
     return np.array([X,Mat3,Mat2])        
                 
-def RotateRBXYZ(Bmat,Cart,oriQ):
+def RotateRBXYZ(Bmat,Cart,oriQ,symAxis=None):
     '''rotate & transform cartesian coordinates to crystallographic ones
     no translation applied. To be used for numerical derivatives 
@@ -986 +986 @@
     :param array Cart: 2D array of coordinates
     :param array Q: quaternion as an np.array
-    
+    :param tuple symAxis: if not None (default), specifies the symmetry
+      axis of the rigid body, which will be aligned to the quaternion vector.
     :returns: 2D array of fractional coordinates, without translation to origin
     '''
+    if symAxis is None:
+        Q = oriQ
+    else:
+        a,v = Q2AV(oriQ)
+        symaxis = np.array(symAxis)
+        xformAng = np.arccos(np.dot(v,symaxis))
+        xformVec = np.cross(symaxis,v)
+        Q = prodQQ(oriQ,AV2Q(xformAng,xformVec))
     XYZ = np.zeros_like(Cart)
     for i,xyz in enumerate(Cart):
-        XYZ[i] = np.inner(Bmat,prodQVQ(oriQ,xyz))
+        XYZ[i] = np.inner(Bmat,prodQVQ(Q,xyz))
     return XYZ
 
 def UpdateRBXYZ(Bmat,RBObj,RBData,RBType):
-    '''returns crystal coordinates for atoms described by RBObj
+    '''returns crystal coordinates for atoms described by RBObj. 
+    Note that RBObj['symAxis'], if present, determines the symmetry
+    axis of the rigid body, which will be aligned to the 
+    quaternion direction.
     
     :param np.array Bmat: see :func:`GSASIIlattice.cell2AB`
@@ -1017 +1029 @@
             QuatA = AVdeg2Q(tor[0],Cart[seq[0]]-Cart[seq[1]])
             Cart[seq[3]] = prodQVQ(QuatA,(Cart[seq[3]]-Cart[seq[1]]))+Cart[seq[1]]
+    # if symmetry axis is defined, place symmetry axis along quaternion
+    if RBObj.get('symAxis') is None:
+        Q = RBObj['Orient'][0]
+    else:
+        a,v = Q2AV(RBObj['Orient'][0])
+        symaxis = np.array(RBObj.get('symAxis'))
+        xformAng = np.arccos(np.dot(v,symaxis))
+        xformVec = np.cross(symaxis,v)
+        Q = prodQQ(RBObj['Orient'][0],AV2Q(xformAng,xformVec))
     XYZ = np.zeros_like(Cart)
     for i,xyz in enumerate(Cart):
-        XYZ[i] = np.inner(Bmat,prodQVQ(RBObj['Orient'][0],xyz))+RBObj['Orig'][0]
+        XYZ[i] = np.inner(Bmat,prodQVQ(Q,xyz))+RBObj['Orig'][0]
     return XYZ,Cart
 

trunk/GSASIIphsGUI.py

@@ -1565 +1565 @@
     '''Update all orientation text on the Add RB panel
     '''
-    for i,sizer in enumerate(G2frame.testRBObjSizers['Osizers']):
-        sizer.SetLabel('%8.5f'%(testRBObj['rbObj']['Orient'][0][i]))
     A,V = G2mth.Q2AVdeg(testRBObj['rbObj']['Orient'][0])
     testRBObj['rbObj']['OrientVec'][0] = A
@@ -9934 +9932 @@
                 G2plt.PlotStructure(G2frame,data)
                 
-            def OnOrien(event):
-                event.Skip()
-                Obj = event.GetEventObject()
-                item = Indx[Obj.GetId()]
-                A,V = G2mth.Q2AVdeg(RBObj['Orient'][0])
-                V = np.inner(Bmat,V)
+            def OnOrien(*args, **kwargs):
+                '''Called when the orientation info is changed (vector 
+                or azimuth)
+                '''
                 try:
-                    val = float(Obj.GetValue())
-                    if item:
-                        V[item-1] = val
-                    else:
-                        A = val
-                    Obj.SetValue('%8.5f'%(val))
-                    V = np.inner(Amat,V)
+                    orient = [float(Indx['Orien'][i].GetValue()) for i in range(4)]
+                    A = orient[0]
+                    V = np.inner(Amat,orient[1:]) # normalized in AVdeg2Q
                     Q = G2mth.AVdeg2Q(A,V)
                     if not any(Q):
@@ -9965 +9957 @@
                 if GSASIIpath.GetConfigValue('debug'): print('patching origin!')
                 RBObj['Orig'][0] = list(RBObj['Orig'][0])
-            Orien,OrienV = G2mth.Q2AVdeg(RBObj['Orient'][0])
-            Orien = [Orien,]
-            Orien.extend(OrienV/nl.norm(OrienV))
-            topSizer.Add(wx.StaticText(RigidBodies,-1,'Origin x,y,z:'),0,WACV)
+            topSizer.Add(wx.StaticText(RigidBodies,-1,'Origin x,y,z (frac)'),0,WACV)
             for ix in range(3):
                 origX = G2G.ValidatedTxtCtrl(RigidBodies,RBObj['Orig'][0],ix,nDig=(8,5),
-                                    typeHint=float,OnLeave=OnOrigX)
+                                    typeHint=float,OnLeave=OnOrigX,
+                                    xmin=-1,xmax=1.,size=(70,-1))
                 topSizer.Add(origX,0,WACV)
             topSizer.Add((5,0),)
@@ -9978 +9968 @@
             Ocheck.SetValue(RBObj['Orig'][1])
             topSizer.Add(Ocheck,0,WACV)
-            topSizer.Add(wx.StaticText(RigidBodies,-1,'Rotation angle, vector:'),0,WACV)
+            topSizer.Add(wx.StaticText(RigidBodies,-1,
+                    'Rotation angle (deg)\n&& Orient. vector (frac)'),0,WACV)
+            Indx['Orien'] = {}
+            Orien,OrienV = G2mth.Q2AVdeg(RBObj['Orient'][0])
+            Orien = [Orien,]
+            Orien.extend(np.inner(Bmat,OrienV)) # fractional coords
+            dp,xmin,xmax = 2,-180.,360.
             for ix,x in enumerate(Orien):
-#            Zstep = G2G.ValidatedTxtCtrl(drawOptions,drawingData,'Zstep',nDig=(10,2),xmin=0.01,xmax=4.0)
-                orien = wx.TextCtrl(RigidBodies,-1,value='%8.4f'%(x),style=wx.TE_PROCESS_ENTER)
-                orien.Bind(wx.EVT_TEXT_ENTER,OnOrien)
-                orien.Bind(wx.EVT_KILL_FOCUS,OnOrien)
-                Indx[orien.GetId()] = ix
+                orien = G2G.ValidatedTxtCtrl(RigidBodies,Orien,ix,nDig=(8,dp),
+                                    typeHint=float,OnLeave=OnOrien,
+                                    xmin=xmin,xmax=xmax,size=(70,-1))
+                dp, xmin,xmax = 4,-1.,1.
+                Indx['Orien'][ix] = orien
                 topSizer.Add(orien,0,WACV)
             Qcheck = wx.ComboBox(RigidBodies,-1,value='',choices=[' ','A','AV'],
@@ -9990 +9986 @@
             Qcheck.Bind(wx.EVT_COMBOBOX,OnOrienRef)
             Qcheck.SetValue(RBObj['Orient'][1])
-            topSizer.Add(Qcheck)
+            topSizer.Add(Qcheck,0,WACV)
             return topSizer
                          
@@ -10030 +10026 @@
             Indx[delRB.GetId()] = rbId
             topLine.Add(delRB,0,WACV)
+            symAxis = RBObj.get('symAxis')
+            if symAxis:
+                if symAxis[0]:
+                    lbl = 'x'
+                elif symAxis[1]:
+                    lbl = 'y'
+                else:
+                    lbl = 'z'                    
+                topLine.Add(wx.StaticText(RigidBodies,-1,
+                    '   Rigid body {} axis is aligned along oriention vector'
+                        .format(lbl)),0,WACV)
             resrbSizer.Add(topLine)
             resrbSizer.Add(LocationSizer(RBObj,'Residue'))
-            resrbSizer.Add(wx.StaticText(RigidBodies,-1,'Torsions:'),0)
+            if len(RBObj['Torsions']):
+                resrbSizer.Add(wx.StaticText(RigidBodies,-1,'Torsions:'),0)
             torSizer = wx.FlexGridSizer(0,6,5,5)
             for itors,tors in enumerate(RBObj['Torsions']):
@@ -10059 +10067 @@
             
         def VecrbSizer(RBObj):
-            G2frame.GetStatusBar().SetStatusText('NB: Rotation vector is in crystallographic space',1)
-                   
             def OnDelVecRB(event):
                 Obj = event.GetEventObject()
@@ -10167 +10173 @@
             G2frame.bottomSizer.Add(ResrbSizer(rbObj))
             mainSizer.Add(G2frame.bottomSizer)
-            G2frame.GetStatusBar().SetStatusText('NB: Rotation vector is in crystallographic space',1)
             G2plt.PlotStructure(G2frame,data)
+            G2plt.PlotStructure(G2frame,data) # draw twice initially for mac
         if 'Vector' in data['RBModels'] and len(data['RBModels']['Vector']):
             nobody = False
@@ -10418 +10424 @@
                 A,V = G2mth.Q2AVdeg(Q)
                 rbObj['OrientVec'][1:] = np.inner(Bmat,V)
-                for i,sizer in enumerate(G2frame.testRBObjSizers['Osizers']):
-                    sizer.SetLabel('%8.5f'%(rbObj['Orient'][0][i]))
                 G2plt.PlotStructure(G2frame,data,False,UpdateTable)
                 UpdateTable()
@@ -10430 +10434 @@
                     np.inner(Amat,rbObj['OrientVec'][1:]))
                 rbObj['Orient'][0] = Q
-                for i,sizer in enumerate(G2frame.testRBObjSizers['Osizers']):
-                    sizer.SetLabel('%8.5f'%(rbObj['Orient'][0][i]))
                 G2frame.testRBObjSizers['OrientVecSiz'][4].SetValue(
                     int(10*rbObj['OrientVec'][0]))                
@@ -10581 +10583 @@
                 G2plt.PlotStructure(G2frame,data,False,UpdateTable)
                 data['testRBObj']['showSelect'].SetLabelText(cryatom)
+            def OnSymRadioSet(event):
+                '''Set the symmetry axis for the body as 
+                data['testRBObj']['rbObj']['symAxis']. This may never be 
+                set, so use data['testRBObj']['rbObj'].get('symAxis') to 
+                access this so the default value is None. 
+                '''
+                data['testRBObj']['rbObj']['symAxis'] = (None,[1,0,0],[0,1,0],[0,0,1])[event.GetEventObject().GetSelection()]
+                UpdateTablePlot()                
             showAtom = [None]
             def showCryAtom(*args,**kwargs):
@@ -10601 +10611 @@
             mainSizer = wx.BoxSizer(wx.VERTICAL)
             mainSizer.Add((5,5),0)
-            Osizers = []
             rbObj = data['testRBObj']['rbObj']
             rbName = rbObj['RBname']
@@ -10650 +10659 @@
                 choice.sort()
                 G2frame.testRBObjSizers.update({'Xsizers':Xsizers})
-            OriSizer.Add(wx.StaticText(RigidBodies,wx.ID_ANY,'Orientation quaternion: '),0,WACV)
-            for ix,x in enumerate(Orien):
-                orien = wx.StaticText(RigidBodies,wx.ID_ANY,'%10.4f'%(x))
-                OriSizer.Add(orien,0,WACV)
-                Osizers.append(orien)
-            G2frame.testRBObjSizers.update({'Osizers':Osizers})
             mainSizer.Add(OriSizer)
             mainSizer.Add((5,5),0)
@@ -10683 +10686 @@
             OriSizer.Add(wx.StaticText(RigidBodies,wx.ID_ANY,' (frac coords)'),0,WACV)
             G2frame.testRBObjSizers.update({'OrientVecSiz':OrientVecSiz})
+            mainSizer.Add(OriSizer)
+            mainSizer.Add((5,5),0)
+            OriSizer = wx.BoxSizer(wx.HORIZONTAL)
+            OriSizer.Add(wx.StaticText(RigidBodies,wx.ID_ANY,
+                                    'Rigid body symmetry axis: '),0, WACV)
+            choices = ['None','x','y','z']
+            symRadioSet = wx.RadioBox(RigidBodies,wx.ID_ANY,choices=choices)
+            symRadioSet.Bind(wx.EVT_RADIOBOX, OnSymRadioSet)
+            OriSizer.Add(symRadioSet)
             mainSizer.Add(OriSizer)
             mainSizer.Add((5,5),0)

trunk/GSASIIplot.py

@@ -8727 +8727 @@
         GL.glLightfv(GL.GL_LIGHT0,GL.GL_AMBIENT,[.2,.2,.2,1])
 
-    def RenderRBtriplet(orig,Q,Bmat):
+    def RenderRBtriplet(orig,Q,Bmat,symAxis=None):
         '''draw an axes triplet located at the origin of a rigid body
         and with the x, y & z axes drawn as red, green and blue. 
@@ -8740 +8740 @@
         GL.glTranslate(*orig)
         GL.glBegin(GL.GL_LINES)
-        lines = G2mth.RotateRBXYZ(Bmat,np.eye(3),Q)
+        lines = G2mth.RotateRBXYZ(Bmat,np.eye(3),Q,symAxis)
         colors = [Rd,Gr,Bl]
         # lines along axial directions
@@ -9363 +9363 @@
             RenderRBtriplet(testRBObj['rbObj']['Orig'][0],
                             testRBObj['rbObj']['Orient'][0],
-                            Bmat)
+                            Bmat,testRBObj['rbObj'].get('symAxis'))
         if len(mcsaModels) > 1 and pageName == 'MC/SA':             #skip the default MD entry
             for ind,[x,y,z] in enumerate(mcsaXYZ):

trunk/GSASIIstrMath.py

@@ -151 +151 @@
                     
 def ApplyRBModelDervs(dFdvDict,parmDict,rigidbodyDict,Phase):
-    'Needs a doc string'
+    'Computes rigid body derivatives'
     atxIds = ['dAx:','dAy:','dAz:']
     atuIds = ['AU11:','AU22:','AU33:','AU12:','AU13:','AU23:']
@@ -177 +177 @@
     RBModels =  Phase['RBModels']
     for irb,RBObj in enumerate(RBModels.get('Vector',[])):
+        symAxis = RBObj.get('symAxis')
         VModel = RBData['Vector'][RBObj['RBId']]
         Q = RBObj['Orient'][0]
@@ -198 +199 @@
             for iv in range(4):
                 Q[iv] -= dx
-                XYZ1 = G2mth.RotateRBXYZ(Bmat,Cart,G2mth.normQ(Q))
+                XYZ1 = G2mth.RotateRBXYZ(Bmat,Cart,G2mth.normQ(Q),symAxis)
                 Q[iv] += 2.*dx
-                XYZ2 = G2mth.RotateRBXYZ(Bmat,Cart,G2mth.normQ(Q))
+                XYZ2 = G2mth.RotateRBXYZ(Bmat,Cart,G2mth.normQ(Q),symAxis)
                 Q[iv] -= dx
                 dXdO = (XYZ2[ia]-XYZ1[ia])/(2.*dx)
@@ -238 +239 @@
 
     for irb,RBObj in enumerate(RBModels.get('Residue',[])):
+        symAxis = RBObj.get('symAxis')
         Q = RBObj['Orient'][0]
         jrb = RRBIds.index(RBObj['RBId'])
@@ -263 +265 @@
             for iv in range(4):
                 Q[iv] -= dx
-                XYZ1 = G2mth.RotateRBXYZ(Bmat,Cart,G2mth.normQ(Q))
+                XYZ1 = G2mth.RotateRBXYZ(Bmat,Cart,G2mth.normQ(Q),symAxis)
                 Q[iv] += 2.*dx
-                XYZ2 = G2mth.RotateRBXYZ(Bmat,Cart,G2mth.normQ(Q))
+                XYZ2 = G2mth.RotateRBXYZ(Bmat,Cart,G2mth.normQ(Q),symAxis)
                 Q[iv] -= dx
                 dXdO = (XYZ2[ia]-XYZ1[ia])/(2.*dx)