|
| 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 heisenberg_model_mpi.cpp -o a.out && |
| 12 | +// mpiexec -np <N> ./a.out |
| 13 | +// ``` |
| 14 | + |
| 15 | +#include "cudaq/algorithms/evolve.h" |
| 16 | +#include "cudaq/algorithms/integrator.h" |
| 17 | +#include "cudaq/operators.h" |
| 18 | +#include "export_csv_helper.h" |
| 19 | +#include <cudaq.h> |
| 20 | + |
| 21 | +int main() { |
| 22 | + cudaq::mpi::initialize(); |
| 23 | + std::cout << "Number of ranks = " << cudaq::mpi::num_ranks() << "\n"; |
| 24 | + // Set up a 15-spin chain, where each spin is a two-level system. |
| 25 | + const int num_spins = 15; |
| 26 | + cudaq::dimension_map dimensions; |
| 27 | + for (int i = 0; i < num_spins; i++) { |
| 28 | + dimensions[i] = 2; // Each spin (site) has dimension 2. |
| 29 | + } |
| 30 | + |
| 31 | + // Initial state |
| 32 | + // Prepare an initial state where the spins are arranged in a staggered |
| 33 | + // configuration. Even indices get the value '0' and odd indices get '1'. For |
| 34 | + // example, for 15 spins: spins: 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 |
| 35 | + std::string spin_state; |
| 36 | + for (int i = 0; i < num_spins; i++) { |
| 37 | + spin_state.push_back((i % 2 == 0) ? '0' : '1'); |
| 38 | + } |
| 39 | + |
| 40 | + // Convert the binary string to an integer index |
| 41 | + // In the Hilbert space of 15 spins (size 2^15 = 32768), this index |
| 42 | + // corresponds to the state |0 1 0 1 0 1 0 1 0 1 0 1 0 1 0> |
| 43 | + int initial_state_index = std::stoi(spin_state, nullptr, 2); |
| 44 | + |
| 45 | + // Build the staggered magnetization operator |
| 46 | + // The staggered magnetization operator is used to measure antiferromagnetic |
| 47 | + // order. It is defined as a sum over all spins of the Z operator, alternating |
| 48 | + // in sign. For even sites, we add `sz`; for odd sites, we subtract `sz`. |
| 49 | + auto staggered_magnetization_t = cudaq::spin_op::empty(); |
| 50 | + for (int i = 0; i < num_spins; i++) { |
| 51 | + auto sz = cudaq::spin_op::z(i); |
| 52 | + if (i % 2 == 0) { |
| 53 | + staggered_magnetization_t += sz; |
| 54 | + } else { |
| 55 | + staggered_magnetization_t -= sz; |
| 56 | + } |
| 57 | + } |
| 58 | + |
| 59 | + // Normalize the number of spins so that the observable is intensive. |
| 60 | + auto staggered_magnetization_op = |
| 61 | + (1 / static_cast<double>(num_spins)) * staggered_magnetization_t; |
| 62 | + |
| 63 | + // Each entry will associate a value of g (the `anisotropy` in the Z coupling) |
| 64 | + // with its corresponding time-series of expectation values of the staggered |
| 65 | + // magnetization. |
| 66 | + std::vector<std::pair<double, std::vector<double>>> observe_results; |
| 67 | + |
| 68 | + // Simulate the dynamics over 100 time steps spanning from time 0 to 5. |
| 69 | + const int num_steps = 100; |
| 70 | + std::vector<double> steps = cudaq::linspace(0.0, 5.0, num_steps); |
| 71 | + |
| 72 | + // For three different values of g, which sets the strength of the Z-Z |
| 73 | + // interaction: g = 0.0 (isotropic in the `XY` plane), 0.25, and 4.0 (strongly |
| 74 | + // `anisotropy`). |
| 75 | + std::vector<double> g_values = {0.0, 0.25, 4.0}; |
| 76 | + |
| 77 | + for (auto g : g_values) { |
| 78 | + // Set the coupling strengths: |
| 79 | + // `Jx` and `Jy` are set to 1.0 (coupling along X and Y axes), while `Jz` is |
| 80 | + // set to the current g value (coupling along the Z axis). |
| 81 | + double Jx = 1.0, Jy = 1.0, Jz = g; |
| 82 | + |
| 83 | + // The Hamiltonian is built from the nearest-neighbor interactions: |
| 84 | + // H = H + `Jx` * `Sx`_i * `Sx`_{i+1} |
| 85 | + // H = H + `Jy` * `Sy`_i * `Sy`_{i+1} |
| 86 | + // H = H + `Jz` * `Sz`_i * `Sz`_{i+1} |
| 87 | + // This is a form of the `anisotropic` Heisenberg (or `XYZ`) model. |
| 88 | + auto hamiltonian = cudaq::spin_op::empty(); |
| 89 | + for (int i = 0; i < num_spins - 1; i++) { |
| 90 | + hamiltonian = |
| 91 | + hamiltonian + Jx * cudaq::spin_op::x(i) * cudaq::spin_op::x(i + 1); |
| 92 | + hamiltonian = |
| 93 | + hamiltonian + Jy * cudaq::spin_op::y(i) * cudaq::spin_op::y(i + 1); |
| 94 | + hamiltonian = |
| 95 | + hamiltonian + Jz * cudaq::spin_op::z(i) * cudaq::spin_op::z(i + 1); |
| 96 | + } |
| 97 | + |
| 98 | + // Initial state vector |
| 99 | + // For a 9-spin system, the Hilbert space dimension is 2^9 = 512. |
| 100 | + // Initialize the state as a vector with all zeros except for a 1 at the |
| 101 | + // index corresponding to our staggered state. |
| 102 | + const int state_size = 1 << num_spins; |
| 103 | + std::vector<std::complex<double>> psi0_data(state_size, {0.0, 0.0}); |
| 104 | + psi0_data[initial_state_index] = {1.0, 0.0}; |
| 105 | + auto psi0 = cudaq::state::from_data(psi0_data); |
| 106 | + |
| 107 | + // The schedule is built using the time steps array. |
| 108 | + cudaq::schedule schedule(steps); |
| 109 | + |
| 110 | + // Use a Runge-`Kutta` integrator (4`th` order) with a small time step `dt` |
| 111 | + // = 0.001. |
| 112 | + cudaq::integrators::runge_kutta integrator(4, 0.001); |
| 113 | + |
| 114 | + // Evolve the initial state psi0 under the Hamiltonian, using the specified |
| 115 | + // schedule and integrator. No collapse operators are included (closed |
| 116 | + // system evolution). Measure the expectation value of the staggered |
| 117 | + // magnetization operator at each time step. |
| 118 | + auto evolve_result = |
| 119 | + cudaq::evolve(hamiltonian, dimensions, schedule, psi0, integrator, {}, |
| 120 | + {staggered_magnetization_op}, true); |
| 121 | + |
| 122 | + // Lambda to extract expectation values for a given observable index |
| 123 | + auto get_expectation = [](int idx, auto &result) -> std::vector<double> { |
| 124 | + std::vector<double> expectations; |
| 125 | + |
| 126 | + auto all_exps = result.expectation_values.value(); |
| 127 | + for (auto exp_vals : all_exps) { |
| 128 | + expectations.push_back((double)exp_vals[idx]); |
| 129 | + } |
| 130 | + return expectations; |
| 131 | + }; |
| 132 | + |
| 133 | + observe_results.push_back({g, get_expectation(0, evolve_result)}); |
| 134 | + } |
| 135 | + |
| 136 | + if (observe_results.size() != 3) { |
| 137 | + std::cerr << "Unexpected number of g values" << std::endl; |
| 138 | + return 1; |
| 139 | + } |
| 140 | + |
| 141 | + if (cudaq::mpi::rank() == 0) { |
| 142 | + // Only save on the first rank to prevent race condition. |
| 143 | + // The `CSV` file "`heisenberg`_model.`csv`" will contain column with: |
| 144 | + // - The time steps |
| 145 | + // - The expectation values of the staggered magnetization for each g |
| 146 | + // value (labeled g_0, g_0.25, g_4). |
| 147 | + export_csv("heisenberg_model_mpi_result.csv", "time", steps, "g_0", |
| 148 | + observe_results[0].second, "g_0.25", observe_results[1].second, |
| 149 | + "g_4", observe_results[2].second); |
| 150 | + |
| 151 | + std::cout << "Simulation complete. The results are saved in " |
| 152 | + "heisenberg_model_mpi_result.csv file." |
| 153 | + << std::endl; |
| 154 | + } |
| 155 | + cudaq::mpi::finalize(); |
| 156 | + return 0; |
| 157 | +} |
0 commit comments