Changeset 4595
- Timestamp:
- Oct 15, 2020 5:43:24 PM (2 years ago)
- Location:
- trunk
- Files:
-
- 4 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/GSASIImath.py
r4583 r4595 979 979 return np.array([X,Mat3,Mat2]) 980 980 981 def RotateRBXYZ(Bmat,Cart,oriQ ):981 def RotateRBXYZ(Bmat,Cart,oriQ,symAxis=None): 982 982 '''rotate & transform cartesian coordinates to crystallographic ones 983 983 no translation applied. To be used for numerical derivatives … … 986 986 :param array Cart: 2D array of coordinates 987 987 :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. 989 990 :returns: 2D array of fractional coordinates, without translation to origin 990 991 ''' 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)) 991 1000 XYZ = np.zeros_like(Cart) 992 1001 for i,xyz in enumerate(Cart): 993 XYZ[i] = np.inner(Bmat,prodQVQ( oriQ,xyz))1002 XYZ[i] = np.inner(Bmat,prodQVQ(Q,xyz)) 994 1003 return XYZ 995 1004 996 1005 def 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. 998 1010 999 1011 :param np.array Bmat: see :func:`GSASIIlattice.cell2AB` … … 1017 1029 QuatA = AVdeg2Q(tor[0],Cart[seq[0]]-Cart[seq[1]]) 1018 1030 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)) 1019 1040 XYZ = np.zeros_like(Cart) 1020 1041 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] 1022 1043 return XYZ,Cart 1023 1044 -
trunk/GSASIIphsGUI.py
r4594 r4595 1565 1565 '''Update all orientation text on the Add RB panel 1566 1566 ''' 1567 for i,sizer in enumerate(G2frame.testRBObjSizers['Osizers']):1568 sizer.SetLabel('%8.5f'%(testRBObj['rbObj']['Orient'][0][i]))1569 1567 A,V = G2mth.Q2AVdeg(testRBObj['rbObj']['Orient'][0]) 1570 1568 testRBObj['rbObj']['OrientVec'][0] = A … … 9934 9932 G2plt.PlotStructure(G2frame,data) 9935 9933 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 ''' 9942 9938 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 9950 9942 Q = G2mth.AVdeg2Q(A,V) 9951 9943 if not any(Q): … … 9965 9957 if GSASIIpath.GetConfigValue('debug'): print('patching origin!') 9966 9958 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) 9971 9960 for ix in range(3): 9972 9961 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)) 9974 9964 topSizer.Add(origX,0,WACV) 9975 9965 topSizer.Add((5,0),) … … 9978 9968 Ocheck.SetValue(RBObj['Orig'][1]) 9979 9969 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. 9981 9977 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()] = ix9978 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 9987 9983 topSizer.Add(orien,0,WACV) 9988 9984 Qcheck = wx.ComboBox(RigidBodies,-1,value='',choices=[' ','A','AV'], … … 9990 9986 Qcheck.Bind(wx.EVT_COMBOBOX,OnOrienRef) 9991 9987 Qcheck.SetValue(RBObj['Orient'][1]) 9992 topSizer.Add(Qcheck )9988 topSizer.Add(Qcheck,0,WACV) 9993 9989 return topSizer 9994 9990 … … 10030 10026 Indx[delRB.GetId()] = rbId 10031 10027 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) 10032 10039 resrbSizer.Add(topLine) 10033 10040 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) 10035 10043 torSizer = wx.FlexGridSizer(0,6,5,5) 10036 10044 for itors,tors in enumerate(RBObj['Torsions']): … … 10059 10067 10060 10068 def VecrbSizer(RBObj): 10061 G2frame.GetStatusBar().SetStatusText('NB: Rotation vector is in crystallographic space',1)10062 10063 10069 def OnDelVecRB(event): 10064 10070 Obj = event.GetEventObject() … … 10167 10173 G2frame.bottomSizer.Add(ResrbSizer(rbObj)) 10168 10174 mainSizer.Add(G2frame.bottomSizer) 10169 G2frame.GetStatusBar().SetStatusText('NB: Rotation vector is in crystallographic space',1)10170 10175 G2plt.PlotStructure(G2frame,data) 10176 G2plt.PlotStructure(G2frame,data) # draw twice initially for mac 10171 10177 if 'Vector' in data['RBModels'] and len(data['RBModels']['Vector']): 10172 10178 nobody = False … … 10418 10424 A,V = G2mth.Q2AVdeg(Q) 10419 10425 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]))10422 10426 G2plt.PlotStructure(G2frame,data,False,UpdateTable) 10423 10427 UpdateTable() … … 10430 10434 np.inner(Amat,rbObj['OrientVec'][1:])) 10431 10435 rbObj['Orient'][0] = Q 10432 for i,sizer in enumerate(G2frame.testRBObjSizers['Osizers']):10433 sizer.SetLabel('%8.5f'%(rbObj['Orient'][0][i]))10434 10436 G2frame.testRBObjSizers['OrientVecSiz'][4].SetValue( 10435 10437 int(10*rbObj['OrientVec'][0])) … … 10581 10583 G2plt.PlotStructure(G2frame,data,False,UpdateTable) 10582 10584 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() 10583 10593 showAtom = [None] 10584 10594 def showCryAtom(*args,**kwargs): … … 10601 10611 mainSizer = wx.BoxSizer(wx.VERTICAL) 10602 10612 mainSizer.Add((5,5),0) 10603 Osizers = []10604 10613 rbObj = data['testRBObj']['rbObj'] 10605 10614 rbName = rbObj['RBname'] … … 10650 10659 choice.sort() 10651 10660 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})10658 10661 mainSizer.Add(OriSizer) 10659 10662 mainSizer.Add((5,5),0) … … 10683 10686 OriSizer.Add(wx.StaticText(RigidBodies,wx.ID_ANY,' (frac coords)'),0,WACV) 10684 10687 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) 10685 10697 mainSizer.Add(OriSizer) 10686 10698 mainSizer.Add((5,5),0) -
trunk/GSASIIplot.py
r4593 r4595 8727 8727 GL.glLightfv(GL.GL_LIGHT0,GL.GL_AMBIENT,[.2,.2,.2,1]) 8728 8728 8729 def RenderRBtriplet(orig,Q,Bmat ):8729 def RenderRBtriplet(orig,Q,Bmat,symAxis=None): 8730 8730 '''draw an axes triplet located at the origin of a rigid body 8731 8731 and with the x, y & z axes drawn as red, green and blue. … … 8740 8740 GL.glTranslate(*orig) 8741 8741 GL.glBegin(GL.GL_LINES) 8742 lines = G2mth.RotateRBXYZ(Bmat,np.eye(3),Q )8742 lines = G2mth.RotateRBXYZ(Bmat,np.eye(3),Q,symAxis) 8743 8743 colors = [Rd,Gr,Bl] 8744 8744 # lines along axial directions … … 9363 9363 RenderRBtriplet(testRBObj['rbObj']['Orig'][0], 9364 9364 testRBObj['rbObj']['Orient'][0], 9365 Bmat )9365 Bmat,testRBObj['rbObj'].get('symAxis')) 9366 9366 if len(mcsaModels) > 1 and pageName == 'MC/SA': #skip the default MD entry 9367 9367 for ind,[x,y,z] in enumerate(mcsaXYZ): -
trunk/GSASIIstrMath.py
r4557 r4595 151 151 152 152 def ApplyRBModelDervs(dFdvDict,parmDict,rigidbodyDict,Phase): 153 ' Needs a doc string'153 'Computes rigid body derivatives' 154 154 atxIds = ['dAx:','dAy:','dAz:'] 155 155 atuIds = ['AU11:','AU22:','AU33:','AU12:','AU13:','AU23:'] … … 177 177 RBModels = Phase['RBModels'] 178 178 for irb,RBObj in enumerate(RBModels.get('Vector',[])): 179 symAxis = RBObj.get('symAxis') 179 180 VModel = RBData['Vector'][RBObj['RBId']] 180 181 Q = RBObj['Orient'][0] … … 198 199 for iv in range(4): 199 200 Q[iv] -= dx 200 XYZ1 = G2mth.RotateRBXYZ(Bmat,Cart,G2mth.normQ(Q) )201 XYZ1 = G2mth.RotateRBXYZ(Bmat,Cart,G2mth.normQ(Q),symAxis) 201 202 Q[iv] += 2.*dx 202 XYZ2 = G2mth.RotateRBXYZ(Bmat,Cart,G2mth.normQ(Q) )203 XYZ2 = G2mth.RotateRBXYZ(Bmat,Cart,G2mth.normQ(Q),symAxis) 203 204 Q[iv] -= dx 204 205 dXdO = (XYZ2[ia]-XYZ1[ia])/(2.*dx) … … 238 239 239 240 for irb,RBObj in enumerate(RBModels.get('Residue',[])): 241 symAxis = RBObj.get('symAxis') 240 242 Q = RBObj['Orient'][0] 241 243 jrb = RRBIds.index(RBObj['RBId']) … … 263 265 for iv in range(4): 264 266 Q[iv] -= dx 265 XYZ1 = G2mth.RotateRBXYZ(Bmat,Cart,G2mth.normQ(Q) )267 XYZ1 = G2mth.RotateRBXYZ(Bmat,Cart,G2mth.normQ(Q),symAxis) 266 268 Q[iv] += 2.*dx 267 XYZ2 = G2mth.RotateRBXYZ(Bmat,Cart,G2mth.normQ(Q) )269 XYZ2 = G2mth.RotateRBXYZ(Bmat,Cart,G2mth.normQ(Q),symAxis) 268 270 Q[iv] -= dx 269 271 dXdO = (XYZ2[ia]-XYZ1[ia])/(2.*dx)
Note: See TracChangeset
for help on using the changeset viewer.