7
$\begingroup$

My question is somewhat related to molecular rotation. I have calculated the polarizability tensor of $\ce{HCHO}$ molecule in PSI4. The output tensor is this,

0     20.834919148327081     1.691888747924788    -0.000000000000510
1      1.691888747924788    10.794718480348456    -0.000000000000070
2     -0.000000000000510    -0.000000000000070    15.501236458190954

However, keeping the bond parameters same, if I rotate the molecule and again calculate the same, the output is as follows,

               0                     1                     2        

0     20.899104179948107    -0.873319044510326    -0.000000000000008
1     -0.873319044510326    15.698712047619225     0.000000000000006
2     -0.000000000000008     0.000000000000006    10.526399871832178

Looking at the diagonal elements, it is clear that depending upon the rotation, the matrix elements changes.

The problem is, this calculation takes very long with larger basis (Q-Zeta). So, I was wondering if there is any mechanism (probably using euler rotation) where I need only to calculate the tensor once and rotate the tensor as the molecule is rotated in the real system.

I would be grateful for any suggestion/help. If more clarity is/are needed, please comment/edit.


Following @Tyberius and @Susi's suggestion. I have gathered a small code.

def kabsch_umeyama(A, B):
    assert A.shape == B.shape
    n, m = A.shape

    EA = np.mean(A, axis=0)
    EB = np.mean(B, axis=0)
    VarA = np.mean(np.linalg.norm(A - EA, axis=1) ** 2) #Variation of reference CoM and points

    H = ((A - EA).T @ (B - EB)) / n                     #Calculation of covariance matrix
    U, D, VT = np.linalg.svd(H)                         #Singular value decomposition
    d = np.sign(np.linalg.det(U) * np.linalg.det(VT))
    S = np.diag([1] * (m - 1) + [d])

    R = U @ S @ VT                                      #R= optimal rotation matrix
    
    c = VarA / np.trace(np.diag(D) @ S)                 #c=Scale factor
    t = EA - c * R @ EB

    return R, c, t

A and B refers to P and Q, two set of points (coordinates of formaldehyde), as P=

[[ 0.          0.6012398  -1.35383976]
 [ 0.         -0.59301814 -1.55209021]
 [ 0.9354225   1.17427624 -1.26515132]
 [-0.9354225   1.17427624 -1.26515132]]

and Q=

      [[ 0.        , -0.60200476,  1.55228866],
       [ 0.        ,  0.59238638,  1.35511328],
       [ 0.        , -1.00937982,  2.57524635],
       [ 0.        , -1.32002906,  0.71694997]]

The rotation matrix R=

      [[ 0.00000112,  0.16405875,  0.98645057],
       [-0.16029666, -0.97369463,  0.16193746],
       [ 0.98706888, -0.15812491,  0.02629698]]

applying scaling and Rotation on Q, as,

print(np.array([t + c * R @ q for q in Q]))

[[-0.0004728   0.6063693  -1.35626887]
 [ 0.0009733  -0.58816112 -1.55025658]
 [ 0.94149768  1.16850837 -1.26498045]
 [-0.94199817  1.17005758 -1.26472672]]

Which suggest the rotation has been performed very well since, rotation on Q produces P very nearly.

Now I have calculated polarizability tensor of P and Q. Which are as follows, P=

      [[22.1081711 ,  1.51237211,  0.        ],
       [ 1.51237211, 13.15927185,  0.        ],
       [ 0.        ,  0.        , 17.44240983]]

and Q =

      [[22.17406111, -0.76701797, -0.        ],
       [-0.76701797, 17.61697616, -0.        ],
       [-0.        , -0.        , 12.91504335]]

My thinking is that, application of Rotation on Polarizability tensor of Q (Q_pol) should give Polarizability tensor of P (P_pol). So, $P_{pol}=RQ_{pol}$ or P_pol=R@Q_pol (in python3) However, I am getting R@Q_pol as,

[[ -0.12581111   2.89021823  12.74005187]
 [ -2.80758665 -17.03060472   2.09142938]
 [ 22.00861039  -3.5427824    0.33962659]]

So, as the problem was, main diagonal elements are on the other direction and it doesn't conserve the trace.

However, $RR^\dagger=I$ (or R@np.transpose(R)=I) as, I =

[[1. 0. 0.]
 [0. 1. 0.]
 [0. 0. 1.]]

Any help?

$\endgroup$
1
  • 1
    $\begingroup$ I added an answer that I think resolves the issue. I left it as a community wiki since the core issue was actually already addressed in Susi's answer (e.g. you were using the wrong formula for rotation of a matrix). $\endgroup$
    – Tyberius
    Feb 8 at 15:15

2 Answers 2

4
$\begingroup$

As Susi already mentioned, rotating a matrix is different from rotating a vector. While a vector can rotated using $Rv$ , a matrix requires you to do $RMR^\dagger$.

With the Kabsch algorithm, its also important to consider which points you are passing in as $A$ and $B$, since it solves for the optimal rotation of $A$ into $B$ and since it allows for scaling and translating the two transformations the $B$ into $A$ transformation will not necessarily be the reverse transformation.

You also have a slight deviation from the typical notation in your code. You write R = U @ S @ VT, but the optimal rotation is usually written R = V @ S @ U.T, which is the transpose of what you have. This isn't an issue as long as you are consistent, but if you are using formulas related to this matrix from the literature, you will have to account for this transpose to make them agree.

Modified code

import numpy as np    
    
def kabsch_umeyama(A, B):    
    assert A.shape == B.shape    
    n, m = A.shape    
    
    EA = np.mean(A, axis=0)    
    EB = np.mean(B, axis=0)    
    VarA = np.mean(np.linalg.norm(A - EA, axis=1) ** 2) #Variation of reference CoM and points    
    
    H = ((A - EA).T @ (B - EB))                      #Calculation of covariance matrix    
    U, D, VT = np.linalg.svd(H)                         #Singular value decomposition    
    d = np.sign(np.linalg.det(U.T) * np.linalg.det(VT.T))    
    S = np.diag([1] * (m - 1) + [d])    
        
    R = VT.T @ S @ U                                      #R= optimal rotation matrix    
    
    c = VarA / np.trace(np.diag(D) @ S)                 #c=Scale factor    
    t = EA - c * R @ EB    
    
    return R, c, t    
    
P=np.array([[ 0.     ,     0.6012398 , -1.35383976],    
 [ 0.       ,  -0.59301814, -1.55209021],    
 [ 0.9354225 ,  1.17427624 ,-1.26515132],    
 [-0.9354225 ,  1.17427624 ,-1.26515132]])    
    
Q=np.array([[ 0.        , -0.60200476,  1.55228866],    
       [ 0.        ,  0.59238638,  1.35511328],    
       [ 0.        , -1.00937982,  2.57524635],    
       [ 0.        , -1.32002906,  0.71694997]])    

R,c,t=kabsch_umeyama(Q,P)
print(R)

Q_pol=np.array([[22.17406111, -0.76701797, -0.        ],
       [-0.76701797, 17.61697616, -0.        ],
       [-0.        , -0.        , 12.91504335]])

print(R@Q_pol@R.T)
print(c,t)

Result

R=[[-0.98709715  0.16007156  0.00403816]
 [ 0.02172461  0.15886771 -0.98706083]
 [-0.1586419  -0.97423721 -0.16029536]]
R@Q_pol@R.T=[[ 2.22995315e+01  3.86328068e-02 -1.47789384e-03]
 [ 3.86328068e-02  1.30327906e+01 -7.24086379e-01]
 [-1.47789384e-03 -7.24086379e-01  1.73737585e+01]]
c=0.2500901704406381 
t=[-0.02221427 -0.94365544  1.63897257]
$\endgroup$
3
  • $\begingroup$ P = t + c * R @ b for b in Q, This formula also changes in order to recover one xyz file from another. Right? $\endgroup$ Feb 13 at 11:13
  • $\begingroup$ @PrasantaBandyopadhyay changes under what circumstances? If you change the order in which P and Q are passed into your function, the formula will technically be the same, but you will get a different R,c, and t depending on whether P or Q is the first argument. $\endgroup$
    – Tyberius
    Feb 14 at 15:03
  • $\begingroup$ Got it @tyberius....... 👍 $\endgroup$ Feb 14 at 15:06
7
$\begingroup$

Tensors are extremely useful in physics, since they allow writing down the laws of physics in a manner that doesn't depend on the used coordinate system. The definition of a tensor is in fact "that it transforms like a tensor".

Let's assume your original coordinate system is defined by the vectors $\hat{e}_i$, and you want to find the tensor ${\bf T}$ in the axis system defined by $\hat{e}'_i = Q_{ji} \hat{e}_j$ (Einstein summation assumed). How do we do this?

A rank-1 tensor is a vector. Let's assume our vector ${\bf v}$ has components $v_i$ in the original coordinate system, that is, ${\bf v} = v_i \hat{e}_i$. The vector has to be the same in the two coordinate systems, ${\bf v}=v_i \hat{e}_i = v_i' \hat{e}'_i$. The components thus transform as $v_i' = {\bf v} \cdot \hat{e}_i' = v_k\hat{e}_k \cdot Q_{ji} \hat{e}_j = v_k Q_{ji} \delta_{kj} = v_j Q_{ji} = Q_{ij}^\text{T} v_j = Q_{ij}^{-1} v_j$; that is, the components transform in the inverse direction of the coordinate system.

A rank-2 tensor is given by ${\bf T} = T_{ij} \hat{e}_i \otimes \hat{e}_j$; its tensor elements form a matrix $T_{ij}$. Repeating the same analysis, one finds that the tensor components transform as $T_{ij}' = Q_{ik}^{-1} Q_{jl}^{-1} T_{kl}$. You didn't specify the rotation matrix for the molecule, but you should be able to reproduce the tensor transform with this equation.

Note that if you employ DFT in your calculations, the quadrature grid may result in breaking of rotational invariance; if that is the case then you need to increase the size of the grid.

$\endgroup$
8
  • $\begingroup$ I am looking at your answer now more closely, but I am using CCSD/def2-QZVPPD for my original research. I was trying to write a program that reads the original tensor and provides the resultant tensor based upon the rotation. Thanks to your insight, I think it will be much easier now. $\endgroup$ Dec 24, 2021 at 5:31
  • $\begingroup$ regarding this, I have calculated the rotation matrix from Kabsch-Umeyama algorithm. taking two molecules in different plane, this algorithm works pretty well. However, when I compare the polarizability tensor between these two systems, there are some significant differences. 1. Main diagonals are rotated to anti-diagonal 2. Minor changes in off-diagonal elements. Any hint on this? $\endgroup$ Jan 3 at 4:14
  • 1
    $\begingroup$ What is anti-diagonal? $\endgroup$ Jan 4 at 16:44
  • 2
    $\begingroup$ It could just be regional differences in terminology, but I have also usually seen "off-diagonal" used to refer to any elements that are not along the top-left to bottom-right diagonal, while "anti-diagonal" would refer to a matrix with only elements along the bottom-left to top-right diagonal. A specific example of anti-diagonal matrices are the so called exchange matrices which are essentially identity matrices reversed along the rows. $\endgroup$
    – Tyberius
    Feb 4 at 2:33
  • 3
    $\begingroup$ @PrasantaBandyopadhyay As Susi said, the fact that your diagonal elements are moving off the diagonal seems to suggest the rotations you are using aren't unitary. You should check that your rotation matrix satisfies $RR^\dagger=I$. If it doesn't, this matrix will distort the geometry/the polarizability tensor rather than just reorient it. $\endgroup$
    – Tyberius
    Feb 4 at 2:36

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.