Examples

The examples below show how to use the C++ API. If you are interested in the Python version, refer to Read the Docs.

Examples of using the API are included in the source tree. The respective files are in the examples folder.

Time evolution of a particle in a box with an exponential initial state with periodic boundary conditions

#include <iostream>
#include <fstream>
#include <iomanip>
#include <sys/stat.h>
#ifdef HAVE_MPI
#include <mpi.h>
#endif
#include "trottersuzuki.h"

#define DIM 640
#define ITERATIONS 1000
#define KERNEL_TYPE "cpu"
#define SNAPSHOTS 10

int main(int argc, char** argv) {
    double length = double(DIM);
    const double particle_mass = 1.;
    bool imag_time = false;
    double delta_t = 5.e-4;

#ifdef HAVE_MPI
    MPI_Init(&argc, &argv);
#endif
    //set lattice
    Lattice2D *grid = new Lattice2D(DIM, length, true, true);
    //set initial state
    State *state = new ExponentialState(grid, 1, 0);
    //set hamiltonian
    Hamiltonian *hamiltonian = new Hamiltonian(grid, NULL, particle_mass);
    //set evolution
    Solver *solver = new Solver(grid, state, hamiltonian, delta_t, KERNEL_TYPE);

    if(grid->mpi_rank == 0) {
        cout << "\n* This source provides an example of the trotter-suzuki program.\n";
        cout << "* It calculates the time-evolution of a particle in a box\n";
        cout << "* with periodic boundary conditions, where the initial\n";
        cout << "* state is the following:\n";
        cout << "* \texp(i2M_PI / L (x + y))\n\n";
    }

    //set file output directory
    stringstream dirname, fileprefix, file_info;
    dirname << "D" << DIM << "_I" << ITERATIONS << "_S" << SNAPSHOTS << "";
    mkdir(dirname.str().c_str(), S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
    file_info  << dirname.str() << "/file_info.txt";
    ofstream out(file_info.str().c_str());

    double mean_X, var_X;
    double mean_Y, var_Y;
    double mean_Px, var_Px;
    double mean_Py, var_Py;

    //get expected values
    double norm2 = solver->get_squared_norm();
    double tot_energy = solver->get_total_energy();
    double kin_energy = solver->get_kinetic_energy();

    mean_X = state->get_mean_x();
    var_X = state->get_mean_xx() - state->get_mean_x() * state->get_mean_x();
    mean_Y = state->get_mean_y();
    var_Y = state->get_mean_yy() - state->get_mean_y() * state->get_mean_y();
    mean_Px = state->get_mean_px();
    var_Px = state->get_mean_pxpx() - state->get_mean_px() * state->get_mean_px();
    mean_Py = state->get_mean_py();
    var_Py = state->get_mean_pypy() - state->get_mean_py() * state->get_mean_py();

    if (grid->mpi_rank == 0) {
        out << std::setw(11) << "time" << std::setw(14) << "squared norm" << std::setw(14) << "tot energy" << std::setw(14) << "kin energy" << std::setw(14) << "pot energy"  << std::setw(14) << "rot energy"
            << std::setw(14) << "<X>" << std::setw(14) << "<(X-<X>)^2>" << std::setw(14) << "<Y>" << std::setw(14) << "<(Y-<Y>)^2>"
            << std::setw(14) << "<Px>" << std::setw(14) << "<(Px-<Px>)^2>" << std::setw(14) << "<Py>" << std::setw(14) << "<(Py-<Py>)^2>\n";
        out << std::setw(11) << "0" << std::setw(14) << norm2 << std::setw(14) << std::setw(14) << tot_energy << std::setw(14) << kin_energy << std::setw(14) << solver->get_potential_energy() << std::setw(14) << solver->get_rotational_energy() << std::setw(14)
            << mean_X << std::setw(14) << var_X << std::setw(14) << mean_Y << std::setw(14) << var_Y << std::setw(14)
            << mean_Px << std::setw(14) << var_Px << std::setw(14) << mean_Py << std::setw(14) << var_Py << endl;
    }

    //evolve and stamp the state
    for(int count_snap = 0; count_snap < SNAPSHOTS; count_snap++) {
        solver->evolve(ITERATIONS, imag_time);
        fileprefix.str("");
        fileprefix << dirname.str() << "/" << ITERATIONS * count_snap;
        state->write_to_file(fileprefix.str());

        //get expected values
        norm2 = solver->get_squared_norm();
        tot_energy = solver->get_total_energy();
        kin_energy = solver->get_kinetic_energy();
        mean_X = state->get_mean_x();
        var_X = state->get_mean_xx() - state->get_mean_x() * state->get_mean_x();
        mean_Y = state->get_mean_y();
        var_Y = state->get_mean_yy() - state->get_mean_y() * state->get_mean_y();
        mean_Px = state->get_mean_px();
        var_Px = state->get_mean_pxpx() - state->get_mean_px() * state->get_mean_px();
        mean_Py = state->get_mean_py();
        var_Py = state->get_mean_pypy() - state->get_mean_py() * state->get_mean_py();

        if(grid->mpi_rank == 0) {
            out << std::setw(11) << (count_snap + 1) * ITERATIONS * delta_t << std::setw(14) << norm2 << std::setw(14) << tot_energy << std::setw(14) << kin_energy << std::setw(14) << solver->get_potential_energy() << std::setw(14) << solver->get_rotational_energy() << std::setw(14) <<
                mean_X << std::setw(14) << var_X << std::setw(14) << mean_Y << std::setw(14) << var_Y << std::setw(14) <<
                mean_Px << std::setw(14) << var_Px << std::setw(14) << mean_Py << std::setw(14) << var_Py << endl;
        }
    }
    out.close();

    if (grid->mpi_rank == 0) {
        cout << "TROTTER " << DIM << "x" << DIM << " kernel:" << KERNEL_TYPE << " np:" << grid->mpi_procs << endl;
    }
    delete solver;
    delete hamiltonian;
    delete state;
    delete grid;
#ifdef HAVE_MPI
    MPI_Finalize();
#endif
    return 0;
}

Imaginary time evolution of an exponential initial state with periodic boundary conditions

#include <iostream>
#include <sys/stat.h>
#ifdef HAVE_MPI
#include <mpi.h>
#endif
#include "trottersuzuki.h"

#define DIM 640
#define ITERATIONS 1000
#define KERNEL_TYPE "cpu"
#define SNAPSHOTS 10

complex<double> super_position_two_exp_state(double x, double y) {
    double L_x = double(DIM);

    return exp(complex<double>(0. , 2.*M_PI / L_x * x)) +
           exp(complex<double>(0. , 10.*2.*M_PI / L_x * x));
}

int main(int argc, char** argv) {
    bool verbose = true;
    double length = double(DIM);
    double delta_t = 0.08;
    const double particle_mass = 1.;
    bool imag_time = true;

#ifdef HAVE_MPI
    MPI_Init(&argc, &argv);
#endif

    //set lattice
    Lattice2D *grid = new Lattice2D(DIM, length, true, true);
    //set initial state
    State *state = new State(grid);
    state->init_state(super_position_two_exp_state);
    //set hamiltonian
    Hamiltonian *hamiltonian = new Hamiltonian(grid, NULL, particle_mass);
    //set evolution
    Solver *solver = new Solver(grid, state, hamiltonian, delta_t, KERNEL_TYPE);

    if(grid->mpi_rank == 0) {
        cout << "\n* This source provides an example of the trotter-suzuki program.\n";
        cout << "* It calculates the imaginary time-evolution of a free particle in a box\n";
        cout << "* with periodic boundary conditions, where the initial\n";
        cout << "* state is the following:\n";
        cout << "* \texp(i2M_PI / L * x) + exp(i20M_PI / L * x)\n\n";
        cout << "* The state will reach the eigenfunction of the Hamiltonian with the lowest\n";
        cout << "* eigenvalue:   exp(i2M_PI / L * x)\n\n";
    }

    //set file output directory
    stringstream dirname, fileprefix;
    dirname << "D" << DIM << "_I" << ITERATIONS << "_S" << SNAPSHOTS << "";
    mkdir(dirname.str().c_str(), S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
    //evolve and stamp the state
    for(int count_snap = 0; count_snap < SNAPSHOTS; count_snap++) {
        solver->evolve(ITERATIONS, imag_time);
        fileprefix.str("");
        fileprefix << dirname.str() << "/" << 1 << "-" << ITERATIONS * count_snap;
        state->write_to_file(fileprefix.str());
    }
    if (grid->mpi_rank == 0 && verbose == true) {
        cout << "TROTTER " << DIM << "x" << DIM << " kernel:" << KERNEL_TYPE << " np:" << grid->mpi_procs << endl;
    }
    delete solver;
    delete hamiltonian;
    delete state;
    delete grid;
#ifdef HAVE_MPI
    MPI_Finalize();
#endif
    return 0;
}

Imaginary time evolution of a Bose-Einstein Condensate trapped in a harmonic potential

#include <iostream>
#include <sys/stat.h>
#ifdef HAVE_MPI
#include <mpi.h>
#endif
#include "trottersuzuki.h"

#define LENGTH 30
#define DIM 640
#define ITERATIONS 1000
#define PARTICLES_NUM 1700000
#define KERNEL_TYPE "cpu"
#define SNAPSHOTS 10
#define SCATTER_LENGTH_2D 5.662739242e-5

int main(int argc, char** argv) {
    bool verbose = true;
    bool imag_time = true;
    double delta_t = 1.e-5;
    const double particle_mass = 1.;
    double coupling_a = 4. * M_PI * double(SCATTER_LENGTH_2D);

#ifdef HAVE_MPI
    MPI_Init(&argc, &argv);
#endif

    //set lattice
    Lattice2D *grid = new Lattice2D(DIM, (double)LENGTH);
    //set initial state
    State *state = new GaussianState(grid, 1., 1., 0., 0., PARTICLES_NUM);
    //set hamiltonian
    Potential *potential = new HarmonicPotential(grid, 1., sqrt(2));
    Hamiltonian *hamiltonian = new Hamiltonian(grid, potential, particle_mass, coupling_a);
    //set solver
    Solver *solver = new Solver(grid, state, hamiltonian, delta_t, KERNEL_TYPE);

    if(grid->mpi_rank == 0) {
        cout << "\n* This source provides an example of the trotter-suzuki program.\n";
        cout << "* It calculates the ground state of a BEC trapped in a harmonic potential.\n";
    }

    //set file output directory
    stringstream dirname, fileprefix;
    dirname << "D" << DIM << "_I" << ITERATIONS << "_S" << SNAPSHOTS << "";
    mkdir(dirname.str().c_str(), S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);

    fileprefix.str("");
    fileprefix << dirname.str() << "/" << 0;
    state->write_to_file(fileprefix.str());
    for(int count_snap = 0; count_snap < SNAPSHOTS; count_snap++) {
        solver->evolve(ITERATIONS, imag_time);
        fileprefix.str("");
        fileprefix << dirname.str() << "/" << ITERATIONS * (count_snap + 1);
        state->write_to_file(fileprefix.str());
    }
    if (grid->mpi_rank == 0 && verbose == true) {
        cout << "TROTTER " << DIM << "x" << DIM << " kernel:" << KERNEL_TYPE << " np:" << grid->mpi_procs << endl;
    }
    delete solver;
    delete hamiltonian;
    delete potential;
    delete state;
    delete grid;
#ifdef HAVE_MPI
    MPI_Finalize();
#endif
    return 0;
}

Time evolution of a two-component Bose-Einstein Condensate trapped in a harmonic potential

#include <iostream>
#include <fstream>
#include <iomanip>
#include <sys/stat.h>

#include "trottersuzuki.h"
#ifdef HAVE_MPI
#include <mpi.h>
#endif

#define LENGTH 10
#define DIM 50
#define ITERATIONS 1000
#define PARTICLES_NUM 1700000
#define KERNEL_TYPE "cpu"
#define SNAPSHOTS 5
#define SNAP_PER_STAMP 5

int main(int argc, char** argv) {
    double particle_mass_a = 1., particle_mass_b = 1.;
    bool imag_time = false;
    double delta_t = 1.e-3;
    double length = double(LENGTH);
    double coupling_a = 0;//7.116007999594e-4;
    double coupling_b = 0;//7.116007999594e-4;
    double coupling_ab = 0;
    double omega_i = 0.;
    double omega_r = 2.*M_PI / 20.;
#ifdef HAVE_MPI
    MPI_Init(&argc, &argv);
#endif
    //set lattice
    Lattice2D *grid = new Lattice2D(DIM, length);
    //set initial state
    State *state1 = new GaussianState(grid, 1, 1, 0., 0., PARTICLES_NUM);
    State *state2 = new GaussianState(grid, 1, 1, 0., 0., PARTICLES_NUM, M_PI / 2.);

    //set hamiltonian
    Potential *potential = new HarmonicPotential(grid, 1., 1.);
    Hamiltonian2Component *hamiltonian = new Hamiltonian2Component(grid, potential, potential, particle_mass_a, particle_mass_b, coupling_a, 0., coupling_ab, coupling_b, 0., omega_r, omega_i);

    //set evolution
    Solver *solver = new Solver(grid, state1, state2, hamiltonian, delta_t, KERNEL_TYPE);

    //set file output directory
    stringstream fileprefix;
    string dirname = "coupledGPE";
    mkdir(dirname.c_str(), S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);

    fileprefix << dirname << "/file_info.txt";
    ofstream out(fileprefix.str().c_str());

    double norm2[2];
    norm2[0] = state1->get_squared_norm();
    norm2[1] = state2->get_squared_norm();
    //get expected values
    double tot_energy = solver->get_total_energy();
    double kin_energy = solver->get_kinetic_energy();
    double rabi_energy = solver->get_rabi_energy();
    double pot_energy = solver->get_potential_energy();
    double intra_energy = solver->get_intra_species_energy();
    double inter_energy = solver->get_inter_species_energy();


    if (grid->mpi_rank == 0) {
        out << std::setw(11) << "time" << std::setw(14) << "squared norm" << std::setw(14) << "sq norm1" << std::setw(14) << "sq norm2" << std::setw(14) << "tot energy" << std::setw(14)
            << "kin energy" << std::setw(14) << "pot energy" << std::setw(14) << "intra energy" << std::setw(14) << "inter energy"  << std::setw(14) << "rabi energy\n";
        out << std::setw(11) << "0" << std::setw(14) << norm2[0] + norm2[1] << std::setw(14) << norm2[0] << std::setw(14) << norm2[1] << std::setw(14) << tot_energy << std::setw(14)
            << kin_energy << std::setw(14) << pot_energy << std::setw(14) << intra_energy << std::setw(14) << inter_energy << std::setw(14) << rabi_energy << endl;
    }

    //write phase and density
    fileprefix.str("");
    fileprefix << dirname << "/1-" << 0;
    state1->write_phase(fileprefix.str());
    state1->write_particle_density(fileprefix.str());
    fileprefix.str("");
    fileprefix << dirname << "/2-" << 0;
    state2->write_phase(fileprefix.str());
    state2->write_particle_density(fileprefix.str());

    for (int count_snap = 0; count_snap < SNAPSHOTS; count_snap++) {
        solver->evolve(ITERATIONS, imag_time);
        //norm calculation
        norm2[0] = state1->get_squared_norm();
        norm2[1] = state2->get_squared_norm();
        //get expected values
        tot_energy = solver->get_total_energy();
        kin_energy = solver->get_kinetic_energy();
        rabi_energy = solver->get_rabi_energy();
        pot_energy = solver->get_potential_energy();
        intra_energy = solver->get_intra_species_energy();
        inter_energy = solver->get_inter_species_energy();

        if(grid->mpi_rank == 0) {
            out << std::setw(11) << (count_snap + 1) * ITERATIONS * delta_t << std::setw(14) << norm2[0] + norm2[1] << std::setw(14) << norm2[0] << std::setw(14) << norm2[1] << std::setw(14) << tot_energy << std::setw(14)
                << kin_energy << std::setw(14) << pot_energy << std::setw(14) << intra_energy << std::setw(14) << inter_energy << std::setw(14) << rabi_energy << endl;
        }

        //stamp phase and particles density
        if(count_snap % SNAP_PER_STAMP == 0.) {
            //write phase and density
            fileprefix.str("");
            fileprefix << dirname << "/1-" << ITERATIONS * (count_snap + 1);
            state1->write_phase(fileprefix.str());
            state1->write_particle_density(fileprefix.str());
            fileprefix.str("");
            fileprefix << dirname << "/2-" << ITERATIONS * (count_snap + 1);
            state2->write_phase(fileprefix.str());
            state2->write_particle_density(fileprefix.str());
        }
    }

    out.close();
    fileprefix.str("");
    fileprefix << dirname << "/" << 1 << "-" << ITERATIONS * SNAPSHOTS;
    state1->write_to_file(fileprefix.str());
    delete solver;
    delete hamiltonian;
    delete potential;
    delete state1;
    delete state2;
    delete grid;
#ifdef HAVE_MPI
    MPI_Finalize();
#endif
    return 0;
}