Ignore:
Timestamp:
May 18, 2021 1:22:50 PM (5 months ago)
Author:
vondreele
Message:

use eigenvalue/vector math to find molecule orientation in onSetPlane

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIIconstrGUI.py

    r4906 r4907  
    22792279
    22802280            def onSetPlane(event):
    2281                 '''Compute least-squares plane for selected atoms;
    2282                 move atoms so that LS plane aligned with x-y plane,
    2283                 with minimum required change to x
     2281                '''Finds eigen vector/matrix for best "ellipsoid" about atoms;
     2282                rotate atoms so that smallest axis is along choice.
    22842283                '''
    22852284                grid.completeEdits()
    22862285                selList = [i==1 for i in rd.Phase['RBselection']]
    22872286                XYZ = rd.Phase['RBcoords'][selList]
    2288                 if len(XYZ) < 3:
    2289                     G2G.G2MessageBox(G2frame,'A plane requires three or more atoms',
    2290                                      'Need more atoms')
     2287                Natoms = len(XYZ)
     2288                if Natoms < 3:
     2289                    G2G.G2MessageBox(G2frame,'A plane requires three or more atoms','Need more atoms')
    22912290                    return
    2292                 # fit 3 ways (in case of singularity) and take result with lowest residual
    2293                 X,Y,Z = [XYZ[:,i] for i in (0,1,2)]
    2294                 XZ = copy.copy(XYZ)
    2295                 XZ[:,1] = 1
    2296                 (a,d,b), resd, rank, sing = nl.lstsq(XZ, -Y,rcond=-1)
    2297                 resid_min = resd
    2298                 normal = a,1,b
    2299                 YZ = copy.copy(XYZ)
    2300                 YZ[:,0] = 1
    2301                 (d,a,b), resd, rank, sing = nl.lstsq(YZ, -X,rcond=-1)
    2302                 if resid_min > resd:
    2303                     resid_min = resd
    2304                     normal = 1,a,b
    2305                 XY = copy.copy(XYZ)
    2306                 XY[:,2] = 1
    2307                 (a,b,d), resd, rank, sing = nl.lstsq(XY, -Z,rcond=-1)
    2308                 if resid_min > resd:
    2309                     resid_min = resd
    2310                     normal = a,b,1
    2311                 # solve for  ax + bx + z + c = 0 or equivalently ax + bx + c = -z
    2312                 # try:
    2313                 # except:
    2314                 #     G2G.G2MessageBox(G2frame,
    2315                 #             'Error computing plane; are atoms in a line?',
    2316                 #             'Computation error')
    2317                 #     return
    2318                 if bntOpts['plane'] == 'xy':
    2319                     # new coordinate system is
    2320                     #   ZP, z' normal to plane
    2321                     #   YP, y' = z' cross x (= [0,1,-b])
    2322                     #   XP, x' = (z' cross x) cross z'
    2323                     # this puts XP as close as possible to X with XP & YP in plane
    2324                     ZP = np.array(normal)
    2325                     ZP /= np.sqrt(np.sum(ZP**2))
    2326                     YP = np.cross(ZP,[1,0,0])
    2327                     if sum(YP*YP) < .1: # pathological condition: z' along x
    2328                         YP = np.cross(ZP,[0,1,0])
    2329                     YP /= np.sqrt(np.sum(YP**2))
    2330                     XP = np.cross(YP,ZP)
    2331                 elif bntOpts['plane'] == 'yz':
    2332                     # new coordinate system is
    2333                     #   XP, x' normal to plane
    2334                     #   ZP, z' = x' cross y
    2335                     #   YP, y' = (x' cross y) cross x'
    2336                     # this puts y' as close as possible to y with z' & y' in plane
    2337                     XP = np.array(normal)
    2338                     XP /= np.sqrt(np.sum(XP**2))
    2339                     ZP = np.cross(XP,[0,1,0])
    2340                     if sum(ZP*ZP) < .1: # pathological condition: x' along y
    2341                         ZP = np.cross(XP,(0,0,1))
    2342                     ZP /= np.sqrt(np.sum(ZP**2))
    2343                     YP = np.cross(ZP,XP)
    2344                 elif bntOpts['plane'] == 'xz':
    2345                     # new coordinate system is
    2346                     #   YP, y' normal to plane
    2347                     #   ZP, z' = x cross y'
    2348                     #   XP, y' = (x cross y') cross z'
    2349                     # this puts XP as close as possible to X with XP & YP in plane
    2350                     YP = np.array(normal)
    2351                     YP /= np.sqrt(np.sum(YP**2))
    2352                     ZP = np.cross([1,0,0],YP)
    2353                     if sum(ZP*ZP) < .1: # pathological condition: y' along x
    2354                         ZP = np.cross([0,1,0],YP)
    2355                     ZP /= np.sqrt(np.sum(ZP**2))
    2356                     XP = np.cross(YP,ZP)
     2291                Zmat = np.zeros((3,3))
     2292                for xyz in XYZ:
     2293                    Zmat += np.outer(xyz.T,xyz)
     2294                Evec,Emat = nl.eig(Zmat)
     2295                Order = np.argsort(np.nan_to_num(Evec))     #short-long order
     2296                if bntOpts['plane'] == 'xy':        #short along z
     2297                    trans = np.array([Emat[Order[2]],Emat[Order[1]],Emat[Order[0]]])
     2298                elif bntOpts['plane'] == 'yz':      #short along x
     2299                    trans = np.array([Emat[Order[0]],Emat[Order[2]],Emat[Order[1]]])
     2300                elif bntOpts['plane'] == 'xz':      #short along y
     2301                    trans = np.array([Emat[Order[1]],Emat[Order[0]],Emat[Order[2]]])
    23572302                else:
    23582303                    print('unexpected plane',bntOpts['plane'])
    23592304                    return
    2360                 trans = np.array((XP,YP,ZP))
    23612305                # update atoms in place
    23622306                rd.Phase['RBcoords'][:] = np.inner(trans,rd.Phase['RBcoords']).T
Note: See TracChangeset for help on using the changeset viewer.