Skip to content

Commit ec9dfce

Browse files
sacpisbettinaheim1tnguyen
authored
Dynamics cpp integration (#2683)
Dynamics cpp integration co-authored by: @1tnguyen --------- Signed-off-by: Sachin Pisal <[email protected]> Signed-off-by: Thien Nguyen <[email protected]> Co-authored-by: Bettina Heim <[email protected]> Co-authored-by: Thien Nguyen <[email protected]>
1 parent 98e4f1a commit ec9dfce

Some content is hidden

Large Commits have some content hidden by default. Use the searchbox below for content that may be hidden.

57 files changed

+5272
-344
lines changed

CMakeLists.txt

+3
Original file line numberDiff line numberDiff line change
@@ -156,6 +156,9 @@ endif()
156156
if(NOT CUTENSORNET_ROOT)
157157
SET(CUTENSORNET_ROOT "$ENV{CUQUANTUM_INSTALL_PREFIX}")
158158
endif()
159+
if(NOT CUDENSITYMAT_ROOT)
160+
SET(CUDENSITYMAT_ROOT "$ENV{CUQUANTUM_INSTALL_PREFIX}")
161+
endif()
159162
if(NOT CUTENSOR_ROOT)
160163
SET(CUTENSOR_ROOT "$ENV{CUTENSOR_INSTALL_PREFIX}")
161164
endif()

docker/build/assets.Dockerfile

+1-1
Original file line numberDiff line numberDiff line change
@@ -150,7 +150,7 @@ RUN source /cuda-quantum/scripts/configure_build.sh && \
150150
libname="$(basename "$binary")" && \
151151
# Linking cublas dynamically necessarily adds a libgcc_s dependency to the GPU-based simulators.
152152
# The same holds for a range of CUDA libraries, whereas libcudart.so does not have any GCC dependencies.
153-
if [ "${libname#libnvqir-custatevec}" != "$libname" ] || [ "${libname#libnvqir-tensornet}" != "$libname" ]; then \
153+
if [ "${libname#libnvqir-custatevec}" != "$libname" ] || [ "${libname#libnvqir-tensornet}" != "$libname" ] || [ "${libname#libnvqir-dynamics}" != "$libname" ]; then \
154154
echo "Skipping validation of $libname."; \
155155
elif [ -n "$(ldd "${binary}" 2>/dev/null | grep gcc)" ]; then \
156156
has_gcc_dependencies=true && \
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,140 @@
1+
/*******************************************************************************
2+
* Copyright (c) 2022 - 2025 NVIDIA Corporation & Affiliates. *
3+
* All rights reserved. *
4+
* *
5+
* This source code and the accompanying materials are made available under *
6+
* the terms of the Apache License 2.0 which accompanies this distribution. *
7+
******************************************************************************/
8+
9+
// Compile and run with:
10+
// ```
11+
// nvq++ --target dynamics cavity_qed.cpp -o a.out && ./a.out
12+
// ```
13+
14+
#include "cudaq/algorithms/evolve.h"
15+
#include "cudaq/dynamics_integrators.h"
16+
#include "cudaq/operators.h"
17+
#include "export_csv_helper.h"
18+
#include <cudaq.h>
19+
20+
int main() {
21+
22+
// Dimension of our composite quantum system:
23+
// subsystem 0 (atom) has 2 levels (ground and excited states).
24+
// subsystem 1 (cavity) has 10 levels (Fock states, representing photon number
25+
// states).
26+
cudaq::dimension_map dimensions{{0, 2}, {1, 10}};
27+
28+
// For the cavity subsystem 1
29+
// We create the annihilation (a) and creation (a+) operators.
30+
// These operators lower and raise the photon number, respectively.
31+
auto a = cudaq::boson_operator::annihilate(1);
32+
auto a_dag = cudaq::boson_operator::create(1);
33+
34+
// For the atom subsystem 0
35+
// We create the annihilation (`sm`) and creation (`sm_dag`) operators.
36+
// These operators lower and raise the excitation number, respectively.
37+
auto sm = cudaq::boson_operator::annihilate(0);
38+
auto sm_dag = cudaq::boson_operator::create(0);
39+
40+
// Number operators
41+
// These operators count the number of excitations.
42+
// For the atom (`subsytem` 0) and the cavity (`subsystem` 1) they give the
43+
// population in each subsystem.
44+
auto atom_occ_op = cudaq::matrix_operator::number(0);
45+
auto cavity_occ_op = cudaq::matrix_operator::number(1);
46+
47+
// Hamiltonian
48+
// The `hamiltonian` models the dynamics of the atom-cavity (cavity QED)
49+
// system It has 3 parts:
50+
// 1. Cavity energy: 2 * pi * photon_number_operator -> energy proportional to
51+
// the number of photons.
52+
// 2. Atomic energy: 2 * pi * atom_number_operator -> energy proportional to
53+
// the atomic excitation.
54+
// 3. Atomic-cavity interaction: 2 * pi * 0.25 * (`sm` * a_dag + `sm_dag` * a)
55+
// -> represents the exchange of energy between the atom and the cavity. It is
56+
// analogous to the Jaynes-Cummings model in cavity QED.
57+
auto hamiltonian = (2 * M_PI * cavity_occ_op) + (2 * M_PI * atom_occ_op) +
58+
(2 * M_PI * 0.25 * (sm * a_dag + sm_dag * a));
59+
60+
// Build the initial state
61+
// Atom (sub-system 0) in ground state.
62+
// Cavity (sub-system 1) has 5 photons (Fock space).
63+
// The overall Hilbert space is 2 * 10 = 20.
64+
const int num_photons = 5;
65+
std::vector<std::complex<double>> initial_state_vec(20, 0.0);
66+
// The index is chosen such that the atom is in the ground state while the
67+
// cavity is in the Fock state with 5 photons.
68+
initial_state_vec[dimensions[0] * num_photons] = 1;
69+
70+
// Define a time evolution schedule
71+
// We define a time grid from 0 to 10 (in arbitrary time units) with 201 time
72+
// steps. This schedule is used by the integrator to simulate the dynamics.
73+
const int num_steps = 201;
74+
std::vector<double> steps = cudaq::linspace(0.0, 10.0, num_steps);
75+
cudaq::schedule schedule(steps);
76+
77+
// Create a CUDA quantum state
78+
// The initial state is converted into a quantum state object for evolution.
79+
auto rho0 = cudaq::state::from_data(initial_state_vec);
80+
81+
// Numerical integrator
82+
// Here we choose a Runge-`Kutta` method for time evolution.
83+
// `dt` defines the time step for the numerical integration, and order 4
84+
// indicates a 4`th` order method.
85+
cudaq::integrators::runge_kutta integrator(4, 0.01);
86+
87+
// Evolve without collapse operators
88+
// This evolution is ideal (closed system) dynamics governed solely by the
89+
// Hamiltonian. The expectation values of the observables (cavity photon
90+
// number and atom excitation probability) are recorded.
91+
cudaq::evolve_result evolve_result =
92+
cudaq::evolve(hamiltonian, dimensions, schedule, rho0, integrator, {},
93+
{cavity_occ_op, atom_occ_op}, true);
94+
95+
// Adding dissipation
96+
// To simulate a realistic scenario, we introduce decay (dissipation).
97+
// Here, the collapse operator represents photon loss from the cavity.
98+
constexpr double decay_rate = 0.1;
99+
auto collapse_operator = std::sqrt(decay_rate) * a;
100+
// Evolve with the collapse operator to incorporate the effect of decay.
101+
cudaq::evolve_result evolve_result_decay =
102+
cudaq::evolve(hamiltonian, dimensions, schedule, rho0, integrator,
103+
{collapse_operator}, {cavity_occ_op, atom_occ_op}, true);
104+
105+
// Lambda to extract expectation values for a given observable index
106+
// Here, index 0 corresponds to the cavity photon number and index 1
107+
// corresponds to the atom excitation probability.
108+
auto get_expectation = [](int idx, auto &result) -> std::vector<double> {
109+
std::vector<double> expectations;
110+
111+
auto all_exps = result.get_expectation_values().value();
112+
for (auto exp_vals : all_exps) {
113+
expectations.push_back((double)exp_vals[idx]);
114+
}
115+
return expectations;
116+
};
117+
118+
// Retrieve expectation values from both the ideal and decaying `evolutions`.
119+
auto ideal_result0 = get_expectation(0, evolve_result);
120+
auto ideal_result1 = get_expectation(1, evolve_result);
121+
auto decay_result0 = get_expectation(0, evolve_result_decay);
122+
auto decay_result1 = get_expectation(1, evolve_result_decay);
123+
124+
// Export the results to `CSV` files
125+
// "cavity_`qed`_ideal_result.`csv`" contains the ideal evolution results.
126+
// "cavity_`qed`_decay_result.`csv`" contains the evolution results with
127+
// cavity decay.
128+
export_csv("cavity_qed_ideal_result.csv", "time", steps,
129+
"cavity_photon_number", ideal_result0,
130+
"atom_excitation_probability", ideal_result1);
131+
export_csv("cavity_qed_decay_result.csv", "time", steps,
132+
"cavity_photon_number", decay_result0,
133+
"atom_excitation_probability", decay_result1);
134+
135+
std::cout << "Simulation complete. The results are saved in "
136+
"cavity_qed_ideal_result.csv "
137+
"and cavity_qed_decay_result.csv files."
138+
<< std::endl;
139+
return 0;
140+
}
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,164 @@
1+
/*******************************************************************************
2+
* Copyright (c) 2022 - 2025 NVIDIA Corporation & Affiliates. *
3+
* All rights reserved. *
4+
* *
5+
* This source code and the accompanying materials are made available under *
6+
* the terms of the Apache License 2.0 which accompanies this distribution. *
7+
******************************************************************************/
8+
9+
// Compile and run with:
10+
// ```
11+
// nvq++ --target dynamics cavity_qed.cpp -o a.out && ./a.out
12+
// ```
13+
14+
#include "cudaq/algorithms/evolve.h"
15+
#include "cudaq/dynamics_integrators.h"
16+
#include "cudaq/operators.h"
17+
#include "export_csv_helper.h"
18+
#include <cudaq.h>
19+
20+
int main() {
21+
22+
// `delta` represents the detuning between the two qubits.
23+
// In physical terms, detuning is the energy difference (or frequency offset)
24+
// between qubit levels. Detuning term (in angular frequency units).
25+
double delta = 100 * 2 * M_PI;
26+
// `J` is the static coupling strength between the two qubits.
27+
// This terms facilitates energy exchange between the qubits, effectively
28+
// coupling their dynamics.
29+
double J = 7 * 2 * M_PI;
30+
// `m_12` models spurious electromagnetic `crosstalk`.
31+
// `Crosstalk` is an unwanted interaction , here represented as a fraction of
32+
// the drive strength applied to the second qubit.
33+
double m_12 = 0.2;
34+
// `Omega` is the drive strength applied to the qubits.
35+
// A driving field can induce transitions between qubit states.
36+
double Omega = 20 * 2 * M_PI;
37+
38+
// For a spin-1/2 system, the raising operator S^+ and lowering operator S^-
39+
// are defined as: S^+ = 0.5 * (X + `iY`) and S^- = 0.5 * (X - `iY`) These
40+
// operators allow transitions between the spin states (|0> and |1>).
41+
auto spin_plus = [](int degree) {
42+
return 0.5 *
43+
(cudaq::spin_operator::x(degree) +
44+
std::complex<double>(0.0, 1.0) * cudaq::spin_operator::y(degree));
45+
};
46+
47+
auto spin_minus = [](int degree) {
48+
return 0.5 *
49+
(cudaq::spin_operator::x(degree) -
50+
std::complex<double>(0.0, 1.0) * cudaq::spin_operator::y(degree));
51+
};
52+
53+
// The Hamiltonian describes the energy and dynamics of our 2-qubit system.
54+
// It consist of several parts:
55+
// 1. Detuning term for qubit 0: (delta / 2) * Z. This sets the energy
56+
// splitting for qubit 0.
57+
// 2. Exchange interaction: J * (S^-_1 * S^+_0 + S^+_1 * S^-_0). This couples
58+
// the two qubits, enabling excitation transfer.
59+
// 3. Drive on qubit 0: Omega * X. A control field that drives transition in
60+
// qubit 0.
61+
// 4. `Crosstalk` drive on qubit 1: m_12 * Omega * X. A reduces drive on qubit
62+
// 1 due to electromagnetic `crosstalk`.
63+
auto hamiltonian =
64+
(delta / 2.0) * cudaq::spin_operator::z(0) +
65+
J * (spin_minus(1) * spin_plus(0) + spin_plus(1) * spin_minus(0)) +
66+
Omega * cudaq::spin_operator::x(0) +
67+
m_12 * Omega * cudaq::spin_operator::x(1);
68+
69+
// Each qubit is a 2-level system (dimension 2).
70+
// The composite system (two qubits) has a total Hilbert space dimension of 2
71+
// * 2 = 4.
72+
cudaq::dimension_map dimensions{{0, 2}, {1, 2}};
73+
74+
// Build the initial state
75+
// psi_00 represents the state |00> (both qubits in the ground state).
76+
// psi_10 represents the state |10> (first qubit excited, second qubit in the
77+
// ground state).
78+
std::vector<std::complex<double>> psi00_data = {
79+
{1.0, 0.0}, {0.0, 0.0}, {0.0, 0.0}, {0.0, 0.0}};
80+
std::vector<std::complex<double>> psi10_data = {
81+
{0.0, 0.0}, {1.0, 0.0}, {0.0, 0.0}, {0.0, 0.0}};
82+
83+
// Two initial state vectors for the 2-qubit system (dimension 4)
84+
auto psi_00 = cudaq::state::from_data(psi00_data);
85+
auto psi_10 = cudaq::state::from_data(psi10_data);
86+
87+
// Create a list of time steps for the simulation.
88+
// Here we use 1001 points linearly spaced between time 0 and 1.
89+
// This schedule will be used to integrate the time evolution of the system.
90+
const int num_steps = 1001;
91+
std::vector<double> steps = cudaq::linspace(0.0, 1.0, num_steps);
92+
cudaq::schedule schedule(steps);
93+
94+
// Use Runge-`Kutta` integrator (4`th` order) to solve the time-dependent
95+
// evolution. `dt` is the integration time step, and `order` sets the accuracy
96+
// of the integrator method.
97+
cudaq::integrators::runge_kutta integrator(4, 0.0001);
98+
99+
// The observables are the spin components along the x, y, and z directions
100+
// for both qubits. These observables will be measured during the evolution.
101+
auto observables = {cudaq::spin_operator::x(0), cudaq::spin_operator::y(0),
102+
cudaq::spin_operator::z(0), cudaq::spin_operator::x(1),
103+
cudaq::spin_operator::y(1), cudaq::spin_operator::z(1)};
104+
105+
// Evolution with 2 initial states
106+
// We evolve the system under the defined Hamiltonian for both initial states
107+
// simultaneously. No collapsed operators are provided (closed system
108+
// evolution). The evolution returns expectation values for all defined
109+
// observables at each time step.
110+
auto evolution_results =
111+
cudaq::evolve(hamiltonian, dimensions, schedule, {psi_00, psi_10},
112+
integrator, {}, observables, true);
113+
114+
// Retrieve the evolution result corresponding to each initial state.
115+
auto &evolution_result_00 = evolution_results[0];
116+
auto &evolution_result_10 = evolution_results[1];
117+
118+
// Lambda to extract expectation values for a given observable index
119+
auto get_expectation = [](int idx, auto &result) -> std::vector<double> {
120+
std::vector<double> expectations;
121+
122+
auto all_exps = result.get_expectation_values().value();
123+
for (auto exp_vals : all_exps) {
124+
expectations.push_back((double)exp_vals[idx]);
125+
}
126+
return expectations;
127+
};
128+
129+
// For the two `evolutions`, extract the six observable trajectories.
130+
// For the |00> initial state, we extract the expectation trajectories for
131+
// each observable.
132+
auto result_00_0 = get_expectation(0, evolution_result_00);
133+
auto result_00_1 = get_expectation(1, evolution_result_00);
134+
auto result_00_2 = get_expectation(2, evolution_result_00);
135+
auto result_00_3 = get_expectation(3, evolution_result_00);
136+
auto result_00_4 = get_expectation(4, evolution_result_00);
137+
auto result_00_5 = get_expectation(5, evolution_result_00);
138+
139+
// Similarly, for the |10> initial state:
140+
auto result_10_0 = get_expectation(0, evolution_result_10);
141+
auto result_10_1 = get_expectation(1, evolution_result_10);
142+
auto result_10_2 = get_expectation(2, evolution_result_10);
143+
auto result_10_3 = get_expectation(3, evolution_result_10);
144+
auto result_10_4 = get_expectation(4, evolution_result_10);
145+
auto result_10_5 = get_expectation(5, evolution_result_10);
146+
147+
// Export the results to a `CSV` file
148+
// Export the Z-component of qubit 1's expectation values for both initial
149+
// states. The `CSV` file "cross_resonance_z.`csv`" will have time versus (Z1)
150+
// data for both |00> and |10> initial conditions.
151+
export_csv("cross_resonance_z.csv", "time", steps, "<Z1>_00", result_00_5,
152+
"<Z1>_10", result_10_5);
153+
// Export the Y-component of qubit 1's expectation values for both initial
154+
// states. The `CSV` file "cross_resonance_y.`csv`" will have time versus (Y1)
155+
// data.
156+
export_csv("cross_resonance_y.csv", "time", steps, "<Y1>_00", result_00_4,
157+
"<Y1>_10", result_10_4);
158+
159+
std::cout
160+
<< "Simulation complete. The results are saved in cross_resonance_z.csv "
161+
"and cross_resonance_y.csv files."
162+
<< std::endl;
163+
return 0;
164+
}

0 commit comments

Comments
 (0)