Trouble computing Duschinsky matrix

Hi, I am having trouble computing the Duschinsky matrix for methane. I think the issue is to do with the modes being given as a 15x15 array. Shouldn’t there only be 9 modes? The rest of the data from GAMESS looks okay.

Excited data
Efreq.txt (86.2 KB)
Ground data
Gfreq.txt (52.3 KB)

Thanks for any help

import numpy as np
import strawberryfields as sf
from strawberryfields.apps import qchem
import matplotlib.pyplot as plt

# r = coords, m = masses, w = freq, L = modes

ri, m, wi, L1 = qchem.read_gamess(r"C:\Users\Public\gamess-64\Methane\Gfreq.log")
rf, mf, wf, L2 = qchem.read_gamess(r'C:\Users\Public\gamess-64\Methane\Efreq.log')

Ud, delta = qchem.duschinsky(Li, Lf, ri, rf, wf, m)

plt.imshow(abs(Ud), cmap="Greens")
plt.colorbar()
plt.xlabel("Mode index")
plt.ylabel("Mode index")
plt.tight_layout()
plt.show()
# operands could not be broadcast together with shapes (15,15) (5,) 

Hello @ReeceOB ! Welcome to the forum!

I believe there are some small typos in your code, I took the liberty of fixing the relevant part. You passing as arguments for the variable Lf and Li instead of L1 and L2 as you defined in previous lines.

import numpy as np
from strawberryfields.apps import qchem
import matplotlib.pyplot as plt

# r = coords, m = masses, w = freq, L = modes

ri, m, wi, Li = qchem.read_gamess(r"C:\Users\Public\gamess-64\Methane\Gfreq.log")
rf, mf, wf, Lf = qchem.read_gamess(r'C:\Users\Public\gamess-64\Methane\Efreq.log')

Ud, delta = qchem.duschinsky(Li, Lf, ri, rf, wf, m)

I agree with you that the problem is related to the modes array shape. It seems that the shape of the array m is (5,0), so you have to match the arrays to the proper dimensions, such that Lf.T * m**0.5 @ (ri-rf).

Does it help? :slight_smile:

Yes, thank you very much.

I have reshaped m so that it is a (5,1) array.

m = np.reshape(m, (m.shape[0], 1)) 

However, Lf should be a (1,9) array, as there are 3N-6 modes for methane. Yet the read_gamess function is outputting a (15,15) array. Can you advise further on how I would compute the correct amount of modes and format them into the correct shape.

Many thanks

I have realised that 6 of the elements are rotational and translational modes which can be disregarded, leaving 9 vibrational modes.

I’m still encountering some difficulties in completing this computation. I’m assuming that the normal modes for the molecule must be represented in an array with the shape (15,9), as there are 5 molecules with 3 Cartesian coordinates and 9 vibronic modes. However, this doesn’t seem to make sense when multiplying this array with m (5,1), as the dimensions don’t match.

Referencing this tutorial using pyrrole and examining the normal modes and mass arrays, the modes appear as expected (3N, 24 modes); however, the masses are given as (30,). It seems that the elements have been triplicated for each atom

If someone could clarify this and explain where I am mistaken that would be very much appreciated :blush:

Sorry for the late reply! I am checking with the dev team, let’s see what is going on :thinking:

Hi @ReeceOB.

Yes the mass vector has been extended to match the size of the frequency and atomic coordinates vector. You can simply triplicate the masses for each atom.

Regarding the normal modes, the read_gammes function is an utility function provided to assist the users and it simply reads the output and returns all the normal modes. But 6 of these modes belong to the transnational and rotational degrees of freedom and should be excluded such that you get 3N - 6 vibrational modes. You should selects the modes manually for your molecule but they are typically the last 3N-6 modes. You can visualise the modes with a software such as Jmol that can read your gamess files to be sure. So, please make sure that you properly reshape the normal mode matrices. This post discusses a similar subject and might be useful.

Please let me know if you have any questions.

Thankyou @sjahangiri

This helps a lot. I have reshaped the matrices to fit. The modes were given as (9,15), but I have taken the transpose of this. However I am still getting the error message

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[11], line 124
    116 wi = np.array([1.24800e+01, 1.21200e+01, 7.51000e+00, 0.00000e+00, 3.00000e-02,
    117        4.00000e-02, 1.67567e+03, 1.67576e+03, 1.67577e+03, 1.90370e+03,
    118        1.90373e+03, 3.52600e+03, 3.78663e+03, 3.78675e+03, 3.78678e+03])
    120 wf = np.array([9.58000e+00, 5.22000e+00, 4.88000e+00, 4.70000e+00, 2.00000e-02,
    121        1.30000e-01, 2.10000e-01, 1.90700e+01, 2.37100e+01, 8.65110e+02,
    122        1.71559e+03, 1.71581e+03, 3.55823e+03, 3.82900e+03, 3.82936e+03])
--> 124 Ud, delta = qchem.duschinsky(Li, Lf, ri, rf, wf, m)
    126 plt.imshow(abs(Ud), cmap="Greens")
    127 plt.colorbar()

File ~\anaconda3\Lib\site-packages\strawberryfields\apps\qchem\utils.py:93, in duschinsky(Li, Lf, ri, rf, wf, m)
     28 r"""Generate the Duschinsky rotation matrix :math:`U` and displacement vector :math:`\delta`.
     29 
     30 The Duschinsky transformation relates the normal coordinates of the initial and
   (...)
     89     :math:`\delta`
     90 """
     91 U = Lf.T @ Li
---> 93 d = Lf.T * m**0.5 @ (ri - rf)
     95 l0_inv = np.diag((h / (wf * 100.0 * c)) ** (-0.5) * 2.0 * pi) / 1.0e10 * m_u**0.5
     97 delta = np.array(d @ l0_inv)
ValueError: operands could not be broadcast together with shapes (9,15) (15,1)

But I do not understand, as I have reshaped all the arrays, see below;
image

Here is my code. I have inserted the data manually rather than using the read_gamess function.

import numpy as np
import strawberryfields as sf
from strawberryfields.apps import qchem
import matplotlib.pyplot as plt

# r = coords, m = masses, w = freq, L = modes

m = np.array([12.     ,  1.00782,  1.00782,  1.00782,  1.00782,12.     ,  1.00782,  1.00782,  1.00782,  1.00782,12.     ,  1.00782,  1.00782,  1.00782,  1.00782])

m = np.reshape(m, (m.shape[0], 1)) 

ri = np.array([-1.34848393e-01,  1.45102586e+00, -7.78600000e-07,
        9.48151120e-01,  1.45094531e+00, -2.18030000e-05,
       -4.95844774e-01,  2.12550214e+00, -7.66575787e-01,
       -4.95881634e-01,  1.77756817e+00,  9.67422229e-01,
       -4.95826320e-01,  4.49888516e-01, -2.00823860e-01])

rf = np.array ([0.18911291,  2.34712294,  0.18027572,
        1.12702634,  1.82103898,  0.07528684,
       -0.37809845,  2.51594969, -0.72380278,
       -0.3765942 ,  2.157796  ,  1.08122229,
       -1.2356966 , -1.58697761, -0.61298207])

Li = np.array([[-5.7550730e-02, -9.6378750e-02,  1.0274670e-02,  6.7000000e-06,
        -1.8810000e-05,  7.1000000e-07, -4.2092350e-02, -7.8580440e-02,
        -1.7383430e-02],
       [-8.4813860e-02,  5.5856730e-02,  4.8922830e-02, -1.1420000e-05,
         7.2300000e-06,  1.9520000e-05, -7.5951850e-02,  3.2305320e-02,
         3.7897600e-02],
       [-4.6922800e-02,  1.7246500e-02, -1.0103026e-01,  8.6500000e-06,
         1.1950000e-05, -3.3000000e-07, -2.6602530e-02,  3.2102340e-02,
        -8.0689500e-02],
       [-5.4276260e-02, -9.0969330e-02,  9.6846400e-03,  1.6900000e-05,
        -3.9380000e-05,  4.9805072e-01,  3.7957733e-01,  7.0841836e-01,
         1.5669270e-01],
       [ 4.1870525e-01, -2.7563321e-01, -2.4155155e-01,  1.7597873e-01,
        -4.6598604e-01, -2.5080000e-05, -3.2751400e-03,  1.3428900e-03,
         1.6125100e-03],
       [ 2.3164028e-01, -8.5139940e-02,  4.9885581e-01,  4.6590732e-01,
         1.7586280e-01, -2.2100000e-06, -1.1425900e-03,  1.3639200e-03,
        -3.4493000e-03],
       [ 2.8494248e-01,  3.6864020e-01, -2.4397831e-01, -2.2023407e-01,
        -4.1462052e-01, -1.6599249e-01, -4.5435940e-02,  6.7879840e-02,
         2.6047695e-01],
       [ 4.1715417e-01, -7.5456950e-02,  1.1938824e-01, -3.8167052e-01,
         7.8489720e-02,  3.1013301e-01,  7.8275000e-02, -1.3172194e-01,
        -4.8643895e-01],
       [ 2.3341572e-01, -3.1292895e-01,  8.8542380e-02, -2.3208301e-01,
         2.6431321e-01, -3.5247987e-01, -9.3794910e-02,  1.5264852e-01,
         5.5124794e-01],
       [ 3.7879001e-01,  3.4966671e-01,  1.0403589e-01,  4.6923181e-01,
         1.6677270e-02, -1.6601344e-01, -1.0038778e-01,  1.9175435e-01,
        -1.6636720e-01],
       [ 3.3304536e-01, -1.6169095e-01, -3.8143984e-01,  3.6423050e-02,
         4.7347115e-01,  1.5017193e-01,  8.5921640e-02, -1.7510752e-01,
         1.5143653e-01],
       [-2.2094120e-02,  2.5307597e-01,  8.4205880e-02,  1.6281085e-01,
        -1.5358352e-01,  4.4487181e-01,  2.6303948e-01, -5.2146573e-01,
         4.4037443e-01],
       [ 7.5790670e-02,  5.2022804e-01,  7.9190600e-03, -2.4909445e-01,
         3.9820662e-01, -1.6605323e-01,  2.6743282e-01, -3.2408650e-02,
        -4.3820950e-02],
       [-1.5904029e-01, -1.5229553e-01, -7.8912780e-02,  1.6940475e-01,
        -8.6060900e-02, -4.6051208e-01,  7.4342420e-01, -7.9167300e-02,
        -1.1785036e-01],
       [ 1.1574001e-01, -6.0358260e-02,  5.3134631e-01, -3.9673811e-01,
        -2.8673472e-01, -9.2385800e-02,  1.4864984e-01, -1.4783800e-02,
        -2.7417010e-02]])

Lf = np.array([[ 2.7984740e-02,  1.2029700e-02, -2.1681160e-02, -4.2138570e-02,
         2.9407500e-02, -7.2854000e-02,  7.1841200e-03,  2.4478200e-03,
        -8.9376960e-02],
       [-5.0101500e-03, -9.7562700e-03, -5.9721690e-02, -1.1824614e-01,
         4.9165700e-03,  3.1043620e-02,  1.9920510e-02,  1.7565100e-02,
         3.1184390e-02],
       [-2.3819040e-02,  2.3526110e-02, -1.2121640e-02, -2.3433470e-02,
        -7.7712070e-02, -2.5602210e-02,  3.9579700e-03, -9.3046850e-02,
         3.5288400e-03],
       [-3.1015851e-01, -1.6636281e-01, -2.0597010e-02,  2.5743664e-01,
         5.1614910e-02, -1.2619182e-01,  5.0655275e-01, -1.8757240e-02,
         6.6736628e-01],
       [-6.2183044e-01, -3.1179174e-01, -5.7887230e-02,  4.3796489e-01,
        -1.0740452e-01, -1.3004910e-01, -2.6243299e-01,  1.1210910e-02,
        -3.7352556e-01],
       [ 4.5869200e-02, -5.6777600e-02, -1.1619380e-02,  8.6733060e-02,
         6.9910964e-01,  2.6575130e-01, -5.2421660e-02, -1.4818600e-03,
        -7.4667720e-02],
       [ 1.0512270e-02,  2.9484893e-01, -2.2763060e-02,  1.2204827e-01,
        -4.9578676e-01,  3.7721102e-01, -2.9649757e-01,  3.4480243e-01,
         2.0855672e-01],
       [ 7.5235200e-03,  7.3556686e-01, -6.3105700e-02,  5.0102207e-01,
         2.0031011e-01,  2.9899650e-02,  1.0824659e-01, -1.0311496e-01,
        -6.2620710e-02],
       [-1.0530860e-02, -1.4710520e-02, -1.2073830e-02,  1.4872750e-02,
         2.9128660e-01, -3.0428475e-01, -4.7932048e-01,  5.4627526e-01,
         3.3828667e-01],
       [ 3.2361560e-02, -2.7793500e-01, -2.0555880e-02,  1.2225306e-01,
         9.4021680e-02,  6.1644112e-01, -2.9559583e-01, -3.5519094e-01,
         1.8827311e-01],
       [ 6.7806630e-01, -3.1486974e-01, -5.8859340e-02,  4.6893594e-01,
        -1.5144619e-01, -2.6948154e-01, -8.3005570e-02, -1.1724050e-01,
         6.4839130e-02],
       [ 1.2246510e-01, -2.2264460e-01, -1.1234940e-02,  1.7741348e-01,
        -6.5091490e-02,  3.4337489e-01,  4.8461518e-01,  5.6309954e-01,
        -3.0563622e-01],
       [-9.8430410e-02,  1.0625540e-02,  3.2209196e-01,  3.7000000e-07,
         6.0000000e-08, -1.6000000e-07,  4.7000000e-07, -0.0000000e+00,
         8.0000000e-08],
       [ 2.3008060e-02, -3.6648300e-03,  8.9105393e-01,  5.7000000e-07,
         2.0000000e-08,  8.0000000e-08,  1.3300000e-06, -3.0000000e-08,
        -4.0000000e-08],
       [ 6.7576590e-02,  1.9393890e-02,  1.7925217e-01,  2.1000000e-07,
        -1.9000000e-07, -7.0000000e-08,  2.6000000e-07,  9.0000000e-08,
        -0.0000000e+00]])

wi = np.array([1.24800e+01, 1.21200e+01, 7.51000e+00, 0.00000e+00, 3.00000e-02,
       4.00000e-02, 1.67567e+03, 1.67576e+03, 1.67577e+03, 1.90370e+03,
       1.90373e+03, 3.52600e+03, 3.78663e+03, 3.78675e+03, 3.78678e+03])

wf = np.array([9.58000e+00, 5.22000e+00, 4.88000e+00, 4.70000e+00, 2.00000e-02,
       1.30000e-01, 2.10000e-01, 1.90700e+01, 2.37100e+01, 8.65110e+02,
       1.71559e+03, 1.71581e+03, 3.55823e+03, 3.82900e+03, 3.82936e+03])

Ud, delta = qchem.duschinsky(Li, Lf, ri, rf, wf, m)

plt.imshow(abs(Ud), cmap="Greens")
plt.colorbar()
plt.xlabel("Mode index")
plt.ylabel("Mode index")
plt.tight_layout()
plt.show()

Thanks again for any help

Hi @ReeceOB and thanks for the code. There are two small issue that you can easily fix:

  1. The m vector should be generated as follows. Please note that the mass for each atom is repeat 3 times to make sure it matches the atomic coordinates, where you have x, y, z for each atom.
m = np.array([12.     ,  1.00782,  1.00782,  1.00782,  1.00782])
m = np.repeat(m, 3)
  1. The frequencies still have 6 components which should be discarded. They belong to the transitional and rotational degrees of freedom. I just discarded the first 6 but please make sure that you visualise your output with Jmol to select the correct modes.
wi = np.array([1.24800e+01, 1.21200e+01, 7.51000e+00, 0.00000e+00, 3.00000e-02,
       4.00000e-02, 1.67567e+03, 1.67576e+03, 1.67577e+03, 1.90370e+03,
       1.90373e+03, 3.52600e+03, 3.78663e+03, 3.78675e+03, 3.78678e+03])[6:]

wf = np.array([9.58000e+00, 5.22000e+00, 4.88000e+00, 4.70000e+00, 2.00000e-02,
       1.30000e-01, 2.10000e-01, 1.90700e+01, 2.37100e+01, 8.65110e+02,
       1.71559e+03, 1.71581e+03, 3.55823e+03, 3.82900e+03, 3.82936e+03])[6:]

With these two corrections, I could generate the matrix. Please let me know if you have any other questions.