Massively Parallel Trotter-Suzuki Solver  1.6.1
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 1
#define PARTICLES_NUM 1700000
#define KERNEL_TYPE "cpu"
#define SNAPSHOTS 1
#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, coupling_ab, coupling_b, 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;
}