Trouble computing Duschinsky matrix

Photonic Software Strawberry Fields
1

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,)
2

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:

3

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

4

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

5

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:

6

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

7

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.

8

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

9

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.

10

Thanks for all the help so far.

I am following the same approach but with different data. I have been able to plot the Duschinsky matrix but run into trouble when trying to plot the spectrum. I am not sure if it is an issue with the data or with the way I am generating the samples.

Thanks for any help :slight_smile:

import numpy as np
import strawberryfields as sf
from strawberryfields import ops
from strawberryfields.apps import qchem
from strawberryfields.apps import plot
from strawberryfields.apps.qchem.vibronic import energies
from strawberryfields.apps.plot import spectrum
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\newv.log")
rf, mf, wf, Lf = qchem.read_gamess(r"C:\Users\Public\gamess-64\wf.log")

# Temperature
T = 0

# Masses
m = np.array([12.     ,  1.00782,  1.00782,  1.00782,  1.00782])
# Triplate masses to match w and r vectors
m = np.repeat(m, 3)

# Coordinates for each atom
Ci = np.array([-0.0000000004, -0.0000000003, -0.0000000008])
H1i = np.array([0.6247332933, 0.6247332928, 0.6247332934])
H2i = np.array([-0.6247332934, -0.6247332932, 0.6247332943])
H3i = np.array([-0.6247332940, 0.6247332942, -0.6247332935])
H4i = np.array([0.6247332945, -0.6247332936, -0.6247332933])

# Combine all atoms into one array
ri = np.vstack([Ci, H1i, H2i, H3i, H4i])

# Flatten the array 
ri = ri.flatten()

# Coordinates for each atom
Cf = np.array([0.0000000005, 0.0000000002, 0.0000000005])
H1f = np.array([0.6252486296, 0.6252486308, 0.6252486301])
H2f = np.array([-0.6252486296, -0.6252486303, 0.6252486294])
H3f = np.array([-0.6252486293, 0.6252486297, -0.6252486298])
H4f = np.array([0.6252486288, -0.6252486304, -0.6252486301])

# Combine all atoms into one array
rf = np.vstack([Cf, H1f, H2f, H3f, H4f])

# Flatten the array 
rf = rf.flatten()

# Remove the first 6 arrays
Li = Li[6:]

# Reshape the array into a (15, 9) array
Li = np.reshape(Li, (15, 9))

# Remove the first 6 arrays
Lf = Lf[6:]

# Reshape the array into a (15, 9) array
Lf = np.reshape(Lf, (15, 9))

# Remove the first 6 arrays
wi = wi[6:]

# Remove the first 6 arrays
wf = wf[6:]

# Compute Duschinsky matrix and displacement vector
Ud, delta = qchem.duschinsky(Li, Lf, ri, rf, wf, m)

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

# map this information to GBS parameters
t, U1, r, U2, alpha = qchem.vibronic.gbs_params(wi, wf, Ud, delta, T)

# Number of modes
n_modes = len(wi)

# Initialize the engine and program objects
eng = sf.Engine("gaussian")
gbs = sf.Program(n_modes)

# Prepare the GBS state
with gbs.context as q:
    # Apply the interferometer U1
    ops.Interferometer(U1) | q

    # Apply the squeezing
    for i in range(n_modes):
        ops.Sgate(r[i]) | q[i]

    # Apply the interferometer U2
    ops.Interferometer(U2) | q

    # Displace the modes
    for i in range(n_modes):
        ops.Dgate(alpha[i]) | q[i]

    # Measure in the Fock basis
    ops.MeasureFock() | q

# Run the engine
results = eng.run(gbs, shots=400)

# Get the samples
samples = results.samples
energies = qchem.vibronic.energies(samples, wi, wf)
plot.spectrum(energies, xmin=-6000, xmax=6000)

newv.txt (48.4 KB)
wf.txt (42.5 KB)

This is the error I am receiving,

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[8], line 108
    105     ops.MeasureFock() | q
    107 # Run the engine
--> 108 results = eng.run(gbs, shots=400)
    110 # Get the samples
    111 samples = results.samples

File ~\anaconda3\Lib\site-packages\strawberryfields\engine.py:570, in LocalEngine.run(self, program, args, compile_options, **kwargs)
    565         if c.op.measurement_deps and eng_run_options["shots"] > 1:
    566             raise NotImplementedError(
    567                 "Feed-forwarding of measurements cannot be used together with multiple shots."
    568             )
--> 570 return super()._run(
    571     program_lst, args=args, compile_options=compile_options, **eng_run_options
    572 )

File ~\anaconda3\Lib\site-packages\strawberryfields\engine.py:306, in BaseEngine._run(self, program, args, compile_options, **kwargs)
    303 p.bind_params(args)
    304 p.lock()
--> 306 _, self.samples, self.samples_dict = self._run_program(p, **kwargs)
    307 self.run_progs.append(p)
    309 if isinstance(p, TDMProgram) and received_rolled:

File ~\anaconda3\Lib\site-packages\strawberryfields\engine.py:430, in LocalEngine._run_program(self, prog, **kwargs)
    427 for cmd in prog.circuit:
    428     try:
    429         # try to apply it to the backend and, if op is a measurement, store it in values
--> 430         val = cmd.op.apply(cmd.reg, self.backend, **kwargs)
    431         if val is not None:
    432             for i, r in enumerate(cmd.reg):

File ~\anaconda3\Lib\site-packages\strawberryfields\ops.py:325, in Measurement.apply(self, reg, backend, **kwargs)
    322 if _shots is None:
    323     return None
--> 325 values = super().apply(reg, backend, **kwargs)
    327 # store the results in the register reference objects
    328 for v, r in zip(np.transpose(values), reg):

File ~\anaconda3\Lib\site-packages\strawberryfields\ops.py:228, in Operation.apply(self, reg, backend, **kwargs)
    226 temp = [rr.ind for rr in reg]
    227 # call the child class specialized _apply method
--> 228 return self._apply(temp, backend, **kwargs)

File ~\anaconda3\Lib\site-packages\strawberryfields\ops.py:1179, in MeasureFock._apply(self, reg, backend, shots, **kwargs)
   1178 def _apply(self, reg, backend, shots=1, **kwargs):
-> 1179     samples = backend.measure_fock(reg, shots=shots, select=self.select, **kwargs)
   1181     if isinstance(samples, list):
   1182         samples = np.array(samples)

File ~\anaconda3\Lib\site-packages\strawberryfields\backends\gaussianbackend\backend.py:235, in GaussianBackend.measure_fock(self, modes, shots, select, **kwargs)
    233     samples = hafnian_sample_state(reduced_cov, shots)
    234 else:
--> 235     samples = hafnian_sample_state(reduced_cov, shots, mean=reduced_mean)
    237 return samples

File ~\anaconda3\Lib\site-packages\thewalrus\samples.py:361, in hafnian_sample_state(cov, samples, mean, hbar, cutoff, max_photons, parallel)
    358     return np.vstack(results)
    360 params = [cov, samples, mean, hbar, cutoff, max_photons]
--> 361 return _hafnian_sample(params)

File ~\anaconda3\Lib\site-packages\thewalrus\samples.py:311, in _hafnian_sample(args)
    308 j = 0
    310 while j < samples:
--> 311     result = generate_hafnian_sample(
    312         cov, mean=mean, hbar=hbar, cutoff=cutoff, max_photons=max_photons
    313     )
    315     if result != -1:
    316         # if result == -1, then you never get anything beyond cutoff
    317         samples_array.append(result)

File ~\anaconda3\Lib\site-packages\thewalrus\samples.py:250, in generate_hafnian_sample(cov, mean, hbar, cutoff, max_photons)
    247     norm_probs = probs.sum()
    248     probs /= norm_probs
--> 250     det_outcome_i = np.random.choice(det_outcomes, p=probs)
    251     det_pattern[mode] = det_outcome_i
    253 if det_pattern[order_inv][-1] == cutoff:

File mtrand.pyx:954, in numpy.random.mtrand.RandomState.choice()

ValueError: probabilities contain NaN

Thanks

11

Hi @ReeceOB. Could you please try simulating one of the molecules in the demos here and see if your code works for them? As you pointed out, the issue could stem from the electronic structure calculations, which are a bit tricky for excited states.

12

Hi,

I have inserted the data for formic acid into the script. The samples are generated fine, but the results are slightly different to that shown in the demo. Could you explain why this is? Could I improve the circuit so the samples generated are more accurate?

Also, referring to this demo, how is the U_vibronic array calculated, and using which parameters? I would like to create a circuit to run on the X8 to be able to compare the results between the two engines, but I understand there are some limitations with only 8 modes and limited gates. Using data such as the equilibrium geometry, normal modes and frequencies of both states, is it possible to create a unitary matrix suitable to be applied to the interferometer of the X8.

Thanks for the help :slight_smile:

import numpy as np
import strawberryfields as sf
from strawberryfields import ops
from strawberryfields.apps import qchem, data, plot
from strawberryfields.apps.qchem.vibronic import energies
from strawberryfields.apps.plot import spectrum

formic = data.Formic()
wi = formic.w  # ground state frequencies
wf = formic.wp  # excited state frequencies
Ud = formic.Ud  # Duschinsky matrix
delta = formic.delta  # displacement vector
T = 0  # temperature

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


# map this information to GBS parameters
t, U1, r, U2, alpha = qchem.vibronic.gbs_params(wi, wf, Ud, delta, T)

# Number of modes
n_modes = len(wi)

# Initialize the engine and program objects
eng = sf.Engine("gaussian")
gbs = sf.Program(n_modes)

# Prepare the GBS state
with gbs.context as q:
    # Apply the interferometer U1
    ops.Interferometer(U1) | q

    # Apply the squeezing
    for i in range(n_modes):
        ops.Sgate(r[i]) | q[i]

    # Apply the interferometer U2
    ops.Interferometer(U2) | q

    # Displace the modes
    for i in range(n_modes):
        ops.Dgate(alpha[i]) | q[i]

    # Measure in the Fock basis
    ops.MeasureFock() | q

# Run the engine
results = eng.run(gbs, shots=20000)

# Get the samples
samples = results.samples
energies = qchem.vibronic.energies(samples, wi, wf)
plot.spectrum(energies, xmin=-1000, xmax=8000)

image
image956×518 42.6 KB

13

Hey @ReeceOB! We’ll be getting to your questions shortly :pray:. Apologies for the delay!

1 Like
14

Hi @ReeceOB. This difference in the data could be due to using different number of samples, if you are generating less samples in your calculations. Please also see here for details on the circuit used in those calculations.

For the details of obtaining U_vibronic in the demo mentioned, please see page 20 this reference.

For using the X8 device with a unitary obtained from data such as the equilibrium geometry, normal modes and frequencies of both states, you make to have many simplifications in the problem, otherwise it will not be possible to compare hardware and simulator results. You can also find some details on this in the reference cited.

Hope these help!

15

Thank you for the reply.

I am generating 20,000 samples, the same as shown in the demo, but my peaks are much broader than when simply using qchem.vibronic.sample(). Is it possible to increase the resolution? Why is the data slightly different if the same operations are being applied in the same order?

Thanks again for all the help

16

Hey @ReeceOB,

I might be missing something, but the demo uses 400k samples. That might be the explanation here.

vibronic_results = eng.run(vibronic_prog, shots=400000)
17

Thanks for the reply Isaac.

I was referring to this demo where 20,000 pre-generated samples are loaded from the data module. When I use the following code to generate samples, I get the same graph.

nr_samples = 20000

# Generate samples and energies using qchem functions
s = qchem.vibronic.sample(t, U1, r, U2, alpha, nr_samples)
e = qchem.vibronic.energies(s, wi, wf)
plot.spectrum(e, xmin=-1000, xmax=8000)

However, when I apply the operations directly in the quantum circuit and generate 20,000 samples as so, the graph is different.

with gbs.context as q:
    # Apply the interferometer U1
    ops.Interferometer(U1) | q

    # Apply the squeezing
    for i in range(n_modes):
        ops.Sgate(r[i]) | q[i]

    # Apply the interferometer U2
    ops.Interferometer(U2) | q

    # Displace the modes
    for i in range(n_modes):
        ops.Dgate(alpha[i]) | q[i]

    # Measure in the Fock basis
    ops.MeasureFock() | q

# Run the engine
results = eng.run(gbs, shots=20000)

# Get the samples
samples = results.samples
energies = qchem.vibronic.energies(samples, wi, wf)
plot.spectrum(energies, xmin=-1000, xmax=8000)

The resulting graph is posted earlier in this thread. If you could explain why the samples have much broader peaks, or if I am generating them incorrectly that would be appreciated.

Thanks :slight_smile:

18

Hi @ReeceOB, did you copy the circuit directly from VibronicTransition? I see a small difference in the application Dgate:

sf.ops.Dgate(np.abs(alpha[i]), np.angle(alpha[i])) | q[i]

Could you please have a look and see if this caused the issue?

Thanks.

19

Hi @sjahangiri, I did not copy the circuit directly. I have made the change in my code but still get the same graph.

# map this information to GBS parameters
t, U1, r, U2, alpha = qchem.vibronic.gbs_params(wi, wf, Ud, delta, T)

# Number of modes
n_modes = len(wi)

# Initialize the engine and program objects
eng = sf.Engine("gaussian")
gbs = sf.Program(n_modes)

# Prepare the GBS state
with gbs.context as q:
    # Apply the interferometer U1
    ops.Interferometer(U1) | q

    # Apply the squeezing
    for i in range(n_modes):
        ops.Sgate(r[i]) | q[i]

    # Apply the interferometer U2
    ops.Interferometer(U2) | q

    # Displace the modes
    for i in range(n_modes):
        ops.Dgate(np.abs(alpha[i]), np.angle(alpha[i])) | q[i]

    # Measure in the Fock basis
    ops.MeasureFock() | q

# Run the engine
results = eng.run(gbs, shots=20000)

# Get the samples
samples = results.samples
energies = qchem.vibronic.energies(samples, wi, wf)
plot.spectrum(energies, xmin=-1000, xmax=8000)

Thanks

20

Hi @ReeceOB, could you please make sure that all gates are applied as listed in the sample function? It is very likely that a set of gates are missing in your manual implementation.

The docstring in the sample function might be also helpful. Please let us know if the issue is not resolved. Thanks!