My own (not S.F. tutorial) Formic acid vibronic spectra sample calculation not computing/stuck

Photonic Software Strawberry Fields
1

Hi there,

I am trying to compute/reproduce my own vibronic spectra of formic acid using my own simulated molecular parameters via Guassian 16. However, when I run my code (identical to tutorial code but with my own data to input for the molecular parameters) the code seems to be stuck and never computes, but also gives no errors. I have tried going stepwise through the working tutorial S.F. code and my own via PyCharm debugging to try and find where my calculation is getting “stuck”, but it has not been apparent and I have not been able to identify the issue.
I have run this calculation on local hardware and on a supercomputer cluster affiliated with my university and it does not appear to be an expense issue, as I have also tried this calcualtion for very few samples. In comparison, if I run the tutorial code and data provided by Strawberry Fields/Xanadu, calculation of 1000+ samples happens very quickly, even on local hardware. However, with my data, a calculation of even 2 samples has not been successfully completed on local hardware or the supercomputing hardware.

I have been trying to solve this for a long time now, I would love any sort of help or assistance. Please let me know what else I can provide.

My data is as such:

Ground state frequencies (w_array):
[ 614.6721 687.0858 1065.1329 1110.3491 1315.0635 1424.1593 1775.7993 3193.522 3644.1215]

Excited state frequencies (wp_array):
[ 328.1703 374.58 845.5295 1028.8746 1199.1784 1223.0602 1368.7532
3223.2968 3666.4988]

Shift vector/displacement vector (deltavector_array):
[ 18.1633 -47.1063 41.5326 6.01971 -17.8149 0.132368
-24.1867 9.39005 -9.9606 ]

Duschinsky matrix (Ud):
[[ 4.15543e-01 -8.63575e-01 -1.62455e-01 -1.87696e-02 2.61311e-02
-2.15749e-01 4.24389e-03 -3.40365e-03 -2.13763e-02]
[ 7.27906e-01 3.03558e-01 -4.15271e-02 -1.90889e-01 1.71165e-02
3.11590e-01 -1.63186e-01 -1.43962e-02 -4.54687e-01]
[ 1.22947e-01 -6.60346e-02 8.93439e-01 1.29866e-02 -7.58260e-02
-1.98558e-01 -4.05041e-02 -3.52725e-01 -6.45187e-02]
[ 2.36294e-01 1.86639e-01 -1.01143e-01 8.97067e-01 -9.70724e-02
-2.73998e-01 -2.82340e-02 3.32640e-04 -3.70688e-02]
[ 1.64662e-01 2.94501e-01 -1.21936e-01 -2.81116e-01 3.95358e-01
-7.20733e-01 3.28468e-01 -9.51115e-03 -2.49967e-02]
[-1.11157e-01 -4.90874e-02 3.10860e-02 1.10414e-01 7.32418e-01
1.09592e-02 -6.56250e-01 -3.44750e-02 9.25105e-03]
[ 5.82274e-02 -1.01300e-01 1.59049e-01 2.33492e-01 5.33609e-01
4.45137e-01 6.40852e-01 3.05726e-02 -7.62547e-03]
[ 4.20004e-02 -8.88412e-03 3.45030e-01 -1.19081e-02 -3.23134e-06
-9.94103e-02 -4.65681e-02 9.17859e-01 -1.36203e-02]
[ 4.24395e-01 1.34269e-01 2.95371e-02 -5.91448e-02 -8.34568e-03
1.16467e-01 -8.15665e-02 1.76883e-02 8.67626e-01]]

My code is as such using these parameters:
t, U1, r, U2, alpha = qchem.vibronic.gbs_params(w_array, wp_array, Ud, deltavector_Array, T)
nr_samples_test = 2
s = qchem.vibronic.sample(t, U1, r, U2, alpha, nr_samples_test)
e = qchem.vibronic.energies(s, w_array, wp_array)
print(np.around(e[:5], 4))

2

Hey @Tony_Hampton! Welcome to the forum :raised_hands:

The code example you sent me is incomplete (I don’t have T). That said, I plugged in a dummy value of 0.1. Here is the code I’m running:

import strawberryfields as sf
from strawberryfields.apps import qchem

import numpy as np

w_array = np.array(
    [
        614.6721,
        687.0858,
        1065.1329,
        1110.3491,
        1315.0635,
        1424.1593,
        1775.7993,
        3193.522,
        3644.1215,
    ]
)

wp_array = np.array(
    [
        328.1703,
        374.58,
        845.5295,
        1028.8746,
        1199.1784,
        1223.0602,
        1368.7532,
        3223.2968,
        3666.4988,
    ]
)

deltavector_array = np.array(
    [
        18.1633,
        -47.1063,
        41.5326,
        6.01971,
        -17.8149,
        0.132368,
        -24.1867,
        9.39005,
        -9.9606,
    ]
)

Ud = np.array(
    [
        [
            4.15543e-01,
            -8.63575e-01,
            -1.62455e-01,
            -1.87696e-02,
            2.61311e-02,
            -2.15749e-01,
            4.24389e-03,
            -3.40365e-03,
            -2.13763e-02,
        ],
        [
            7.27906e-01,
            3.03558e-01,
            -4.15271e-02,
            -1.90889e-01,
            1.71165e-02,
            3.11590e-01,
            -1.63186e-01,
            -1.43962e-02,
            -4.54687e-01,
        ],
        [
            1.22947e-01,
            -6.60346e-02,
            8.93439e-01,
            1.29866e-02,
            -7.58260e-02,
            -1.98558e-01,
            -4.05041e-02,
            -3.52725e-01,
            -6.45187e-02,
        ],
        [
            2.36294e-01,
            1.86639e-01,
            -1.01143e-01,
            8.97067e-01,
            -9.70724e-02,
            -2.73998e-01,
            -2.82340e-02,
            3.32640e-04,
            -3.70688e-02,
        ],
        [
            1.64662e-01,
            2.94501e-01,
            -1.21936e-01,
            -2.81116e-01,
            3.95358e-01,
            -7.20733e-01,
            3.28468e-01,
            -9.51115e-03,
            -2.49967e-02,
        ],
        [
            -1.11157e-01,
            -4.90874e-02,
            3.10860e-02,
            1.10414e-01,
            7.32418e-01,
            1.09592e-02,
            -6.56250e-01,
            -3.44750e-02,
            9.25105e-03,
        ],
        [
            5.82274e-02,
            -1.01300e-01,
            1.59049e-01,
            2.33492e-01,
            5.33609e-01,
            4.45137e-01,
            6.40852e-01,
            3.05726e-02,
            -7.62547e-03,
        ],
        [
            4.20004e-02,
            -8.88412e-03,
            3.45030e-01,
            -1.19081e-02,
            -3.23134e-06,
            -9.94103e-02,
            -4.65681e-02,
            9.17859e-01,
            -1.36203e-02,
        ],
        [
            4.24395e-01,
            1.34269e-01,
            2.95371e-02,
            -5.91448e-02,
            -8.34568e-03,
            1.16467e-01,
            -8.15665e-02,
            1.76883e-02,
            8.67626e-01,
        ],
    ]
)

t, U1, r, U2, alpha = qchem.vibronic.gbs_params(w_array, wp_array, Ud, deltavector_array, 0.1)
nr_samples_test = 2
s = qchem.vibronic.sample(t, U1, r, U2, alpha, nr_samples_test)
e = qchem.vibronic.energies(s, w_array, wp_array)
print(np.around(e[:5], 4))

Indeed, it is taking a while compared to this:

import strawberryfields as sf
from strawberryfields.apps import qchem
from strawberryfields.apps import data

formic = data.Formic()
w = formic.w
wp = formic.wp
Ud = formic.Ud
delta = formic.delta
T = 0
t, U1, r, U2, alpha = qchem.vibronic.gbs_params(w, wp, Ud, delta, T)

nr_samples_test = 2
s = qchem.vibronic.sample(t, U1, r, U2, alpha, nr_samples_test)

I’ll have to consult some colleagues about this. Sit tight!

3

Just consulted someone internally, and it looks like, given your parameters, the mean photon number on all modes is above 900, which is quite computationally intensive! Embedding your molecule into a GBS device with a lower photon number is necessary here, and I’m unsure of a method for doing so :sweat_smile:

4

Hi Isaac,

Thank you for your response! I have a few follow up questions:

Can you provide insight into how photons are assigned to modes? As in, what would be the source of the photons assigned to modes and how would that result in me obtaining over 900 mean photon count per mode? Is this process a part of the Strawberry Fields framework or is this arising from the simulation side of things, i.e., my choice for basis set and level of simulation in Gaussian?

Additionally, would computational issues here also arise from my inclusion of all 9 vibrational modes of formic acid? The tutorial and Nature Phtonics paper it is based off of only uses the 7 vibrational modes that contribute to the totally symmetric portion of the Duschinsky matrix. Would me including the additional 2 vibrational modes that make up the non-totally symmetric portion of the Duschinsky matrix be poses some of these issues?

Lastly, how is it that you determined the mean number of photons assigned to each mode? Is this something that is easily done in the Python code?

Sorry for the many questions, but thank you very much in advance for your help.

Tony

5

Hi @Tony_Hampton,

I cannot provide an answer to your first two questions but can happily describe how I computed the mean photon number: to do so I basically took the implementation of strawberryfields.apps.qchem.vibronic.sample but instead of sampling I compute the covariance and means of the gaussian circuit before Fock detection; using this you can then compute the mean photon number per mode. Here is the code

import strawberryfields as sf
from strawberryfields.apps.qchem.vibronic import VibronicTransition
from thewalrus.quantum import photon_number_mean

def vibronic_cov_and_means(
    t: np.ndarray,
    U1: np.ndarray,
    r: np.ndarray,
    U2: np.ndarray,
    alpha: np.ndarray,
    n_samples: int,
    loss: float = 0.0,
) -> list:
    if n_samples < 1:
        raise ValueError("Number of samples must be at least one")
    if not 0 <= loss <= 1:
        raise ValueError("Loss parameter must take a value between zero and one")

    n_modes = len(t)

    eng = sf.LocalEngine(backend="gaussian")

    if np.any(t != 0):
        gbs = sf.Program(n_modes * 2)
    else:
        gbs = sf.Program(n_modes)

    with gbs.context as q:

        if np.any(t != 0):
            for i in range(n_modes):
                sf.ops.S2gate(t[i]) | (q[i], q[i + n_modes])

        VibronicTransition(U1, r, U2, alpha) | q[:n_modes]

        if loss:
            for _q in q:
                sf.ops.LossChannel(1 - loss) | _q

    result = eng.run(gbs)

    return result.state.cov(), result.state.means()


cov, means = vibronic_cov_and_means(t, U1, r, U2, alpha, nr_samples_test)
photon_means_per_mode = [photon_number_mean(means, cov, mode) for mode in range(U1.shape[0])]
print(f"mean photon number per mode: {photon_means_per_mode}")
1 Like