359 lines
14 KiB
C++
359 lines
14 KiB
C++
#ifdef _MSC_VER //If the compiler is the Visual Studio Compiler
|
|
#define NOMINMAX //solves a problem with std::max on Windows
|
|
#pragma warning(disable: 4297) //remove throw warnings
|
|
#endif
|
|
|
|
#include "common/definitions.h"
|
|
#include <iostream>
|
|
#include <vector>
|
|
|
|
#include "common/SpinSystem.h"
|
|
|
|
const bool debug = false;
|
|
|
|
const char* get_version()
|
|
{
|
|
return SPINTRUM_VERSION;
|
|
}
|
|
|
|
void print_software_info()
|
|
{
|
|
std::cout<<"Spintrum version "<<get_version()<<". ";
|
|
std::cout<<"This software is provided under a license that should have been delivered with its package in the file LICENSE.txt."
|
|
"For more information, contact Samer Afach: samer [[at]] afach.de."<<std::endl<<std::endl;
|
|
}
|
|
|
|
/**
|
|
* A quantum mechanics operation (populate spins, tip spins, etc...) is a struct that contains a specific, predefined number of parameters.
|
|
* Each operation has an Identifier, which is simply an integer. From this identifier, the code knows what parameters to seek
|
|
* and how many parameters should be expected.
|
|
* There are integer parameters, real parameters, and array of real parameters.
|
|
*/
|
|
struct SpinOperation
|
|
{
|
|
long OperationIdentifier;
|
|
long Params_int_count; //number of integer parameters
|
|
long* Params_int; //integer parameters
|
|
long Params_real_count; //number of real parameters
|
|
double* Params_real; //real parameters
|
|
long Params_array_count; //number of arrays
|
|
long* Params_array_lengths; //lengths of suppled arrays
|
|
char** Params_array_names; //names of the parameters of the arrays
|
|
double** Params_array; //array data
|
|
};
|
|
const static int OPERATION__THERMAL_POPULATE = 0;
|
|
const static int OPERATION__TIP_SPINS = 1;
|
|
const static int OPERATION__SET_HAMILTONIAN = 2;
|
|
const static int OPERATION__EVOLVE_TIME_INDEPENDENT = 3;
|
|
const static int OPERATION__INIT_TIME_INDEPENDENT_EVOLUTION = 4;
|
|
const static int OPERATION__EVOLVE_TIME_DEPENDENT = 5;
|
|
|
|
typedef double Real;
|
|
|
|
typedef SpinSystem<Real,std::complex<Real>,VEC,CX_VEC,MAT,CX_MAT> SpinSystemType;
|
|
|
|
struct SpinInfo
|
|
{
|
|
CX_MAT jCouplings;
|
|
CX_VEC gammas;
|
|
G_MAT<int> multiplicities;
|
|
};
|
|
|
|
void devectorizeJCouplings(const long numOfSpins, const double* const jCouplingsInput, CX_MAT& jCouplingsOutput)
|
|
{
|
|
if(jCouplingsInput == 0)
|
|
{
|
|
jCouplingsOutput.resize(0,0);
|
|
return;
|
|
}
|
|
int i = 0; //element of given array
|
|
int j = 0; //element of current row
|
|
int row = 0;
|
|
while (i < numOfSpins*(numOfSpins-1)/2)
|
|
{
|
|
if(j <= row)
|
|
{
|
|
j++;
|
|
continue;
|
|
}
|
|
if(j >= numOfSpins)
|
|
{
|
|
j = 0;
|
|
row++;
|
|
continue;
|
|
}
|
|
jCouplingsOutput.at(row,j) = jCouplingsInput[i];
|
|
j++;
|
|
i++;
|
|
}
|
|
}
|
|
extern "C" int InitializeSimulator(double* gammas, double* jCouplings, int* spinMultiplicities, int numOfSpins, bool doPrintOutput, void** SpinSystemObject)
|
|
{
|
|
try
|
|
{
|
|
print_software_info();
|
|
// std::cout<<"Number of spins: "<<numOfSpins<<std::endl;
|
|
SpinInfo spinInfo;
|
|
spinInfo.gammas.resize(numOfSpins,1);
|
|
spinInfo.multiplicities.resize(numOfSpins,1);
|
|
spinInfo.jCouplings.resize(numOfSpins,numOfSpins);
|
|
for(int i = 0; i < numOfSpins; i++)
|
|
{
|
|
// std::cout<<gammas[i]<<"\t";
|
|
spinInfo.gammas[i] = gammas[i];
|
|
}
|
|
// std::cout<<spinInfo.gammas<<std::endl;
|
|
if(spinMultiplicities == 0)
|
|
{
|
|
for(int i = 0; i < numOfSpins; i++)
|
|
spinInfo.multiplicities[i] = 2;
|
|
}
|
|
else
|
|
{
|
|
for(int i = 0; i < numOfSpins; i++)
|
|
spinInfo.multiplicities[i] = spinMultiplicities[i];
|
|
}
|
|
|
|
devectorizeJCouplings(numOfSpins, jCouplings, spinInfo.jCouplings);
|
|
// std::cout<<spinInfo.jCouplings<<std::endl;
|
|
|
|
SpinSystemType *Spins = new SpinSystem<Real,std::complex<Real>,VEC,CX_VEC,MAT,CX_MAT>();
|
|
Spins->print = doPrintOutput;
|
|
Spins->init(spinInfo.gammas,spinInfo.jCouplings,spinInfo.multiplicities);
|
|
*SpinSystemObject = static_cast<void*>(Spins);
|
|
return 0;
|
|
}
|
|
catch(std::exception &ex)
|
|
{
|
|
std::cout<<"An exception was thrown while trying to initialize the spin system. Exception message is: " << ex.what() << std::endl;
|
|
return -1;
|
|
}
|
|
catch(...)
|
|
{
|
|
std::cout<<"An unhandled exception was thrown while trying to initialize the spin system."<< std::endl;
|
|
return -2;
|
|
}
|
|
}
|
|
|
|
extern "C" int DenitializeSimulator(void* SpinSystemObject)
|
|
{
|
|
try
|
|
{
|
|
delete (static_cast<SpinSystemType*>(SpinSystemObject));
|
|
return 0;
|
|
}
|
|
catch(std::exception &ex)
|
|
{
|
|
std::cout<<"An exception was thrown while trying to delete the spin system from memory. Exception message is: " << ex.what() << std::endl;
|
|
return -1;
|
|
}
|
|
catch(...)
|
|
{
|
|
std::cout<<"An unhandled exception was thrown while trying to delete the spin system from memory."<< std::endl;
|
|
return -2;
|
|
}
|
|
}
|
|
|
|
|
|
extern "C" int SetGammasInSimulator(double* gammas, void* SpinSystemObject)
|
|
{
|
|
try
|
|
{
|
|
CX_VEC gammasVec((static_cast<SpinSystemType*>(SpinSystemObject))->getNumOfSpins(),1);
|
|
for(int i = 0; i < gammasVec.size(); i++)
|
|
{
|
|
gammasVec[i] = gammas[i];
|
|
}
|
|
(static_cast<SpinSystemType*>(SpinSystemObject))->setGammas(gammasVec);
|
|
return 0;
|
|
}
|
|
catch(std::exception &ex)
|
|
{
|
|
std::cout<<"An exception was thrown while trying to delete the spin system from memory. Exception message is: " << ex.what() << std::endl;
|
|
return -1;
|
|
}
|
|
catch(...)
|
|
{
|
|
std::cout<<"An unhandled exception was thrown while trying to delete the spin system from memory."<< std::endl;
|
|
return -2;
|
|
}
|
|
}
|
|
|
|
extern "C" int SetJCouplingsInSimulator(double* jCouplings, void* SpinSystemObject)
|
|
{
|
|
try
|
|
{
|
|
CX_MAT jCouplingsMat((static_cast<SpinSystemType*>(SpinSystemObject))->getNumOfSpins(),(static_cast<SpinSystemType*>(SpinSystemObject))->getNumOfSpins());
|
|
devectorizeJCouplings((static_cast<SpinSystemType*>(SpinSystemObject))->getNumOfSpins(), jCouplings, jCouplingsMat);
|
|
(static_cast<SpinSystemType*>(SpinSystemObject))->setJCouplings(jCouplingsMat);
|
|
return 0;
|
|
}
|
|
catch(std::exception &ex)
|
|
{
|
|
std::cout<<"An exception was thrown while trying to delete the spin system from memory. Exception message is: " << ex.what() << std::endl;
|
|
return -1;
|
|
}
|
|
catch(...)
|
|
{
|
|
std::cout<<"An unhandled exception was thrown while trying to delete the spin system from memory."<< std::endl;
|
|
return -2;
|
|
}
|
|
}
|
|
void _do_operation_evolve_time_dependent(SpinOperation* op, SpinSystemType* SpinsPtr, double* outputArray, long ¤tDataPoint)
|
|
{
|
|
if(debug) std::cout<<"Printing received data:"<<std::endl;
|
|
if(op->Params_array_count == 0)
|
|
{
|
|
throw std::runtime_error("There's no data provided for the time-dependent evolution.");
|
|
}
|
|
else
|
|
{
|
|
//make sure that all lengths are equal
|
|
long initialLength = op->Params_array_lengths[0];
|
|
for(long i = 1; i < op->Params_array_count; i++)
|
|
{
|
|
if(op->Params_array_lengths[i] != initialLength)
|
|
{
|
|
throw std::length_error("Data arrays in time-dependent evolution are expected all to be of the same length.");
|
|
}
|
|
}
|
|
}
|
|
std::vector<std::string> varsNames(static_cast<std::vector<std::string>::size_type>(op->Params_array_count));
|
|
std::copy(op->Params_array_names,op->Params_array_names+op->Params_array_count,varsNames.begin()); //copy names to vector<string>
|
|
|
|
CX_MAT evData(op->Params_array_lengths[0],op->Params_array_count);
|
|
for(long i = 0; i < op->Params_array_count; i++)
|
|
{
|
|
std::copy(&op->Params_array[i][0],&op->Params_array[i][op->Params_array_lengths[i]],evData.begin()+evData.rows()*i);
|
|
}
|
|
std::set<char> dirString;
|
|
if(op->Params_int[1]) dirString.insert('x');
|
|
if(op->Params_int[2]) dirString.insert('y');
|
|
if(op->Params_int[3]) dirString.insert('z');
|
|
std::vector<Real> tr = SpinsPtr->evolveDensityMatrix(op->Params_real[0],varsNames,evData,static_cast<int>(op->Params_int[0]),dirString);
|
|
std::copy(std::make_move_iterator(tr.begin()), std::make_move_iterator(tr.end()), outputArray + currentDataPoint);
|
|
currentDataPoint += op->Params_array_lengths[0]*static_cast<long>(dirString.size());
|
|
}
|
|
|
|
long ExecuteOperations(void* SpinSystemObject, void** ops, int numOfOperations, double* output, bool debug)
|
|
{
|
|
if(debug)
|
|
{
|
|
std::cout<<"Number of operations: "<<numOfOperations<<std::endl;
|
|
}
|
|
|
|
SpinSystemType* SpinsPtr = static_cast<SpinSystemType*>(SpinSystemObject);
|
|
long currentDataPoint = 0; //for a simulation with many operations that consecutively fill the output, this keeps track of the current point
|
|
for(int i = 0; i < numOfOperations; i++) {
|
|
void** currOperationBytePosition = (static_cast<void**>(ops));
|
|
long OperationIdentifier = *(static_cast<long*>(currOperationBytePosition[i]));
|
|
if(debug) std::cout<<"Identifier: "<<OperationIdentifier<<std::endl;
|
|
if(OperationIdentifier == OPERATION__THERMAL_POPULATE) {
|
|
SpinOperation* op;
|
|
op = static_cast<SpinOperation*>(currOperationBytePosition[i]);
|
|
SpinsPtr->populate_thermal(op->Params_real[0],op->Params_real[1],op->Params_real[2],op->Params_real[3]);
|
|
}
|
|
else if (OperationIdentifier == OPERATION__TIP_SPINS) {
|
|
SpinOperation* op;
|
|
op = static_cast<SpinOperation*>(currOperationBytePosition[i]);
|
|
SpinsPtr->tipSpins(static_cast<int>(op->Params_int[0]),op->Params_real[0]);
|
|
}
|
|
else if (OperationIdentifier == OPERATION__SET_HAMILTONIAN) {
|
|
SpinOperation* op;
|
|
op = static_cast<SpinOperation*>(currOperationBytePosition[i]);
|
|
SpinsPtr->calculateHamiltonians(op->Params_real[0],op->Params_real[1],op->Params_real[2]);
|
|
}
|
|
else if (OperationIdentifier == OPERATION__EVOLVE_TIME_INDEPENDENT) {
|
|
SpinOperation* op;
|
|
op = static_cast<SpinOperation*>(currOperationBytePosition[i]);
|
|
std::vector<Real> trVals = SpinsPtr->solveMagnetization_parallel(op->Params_int[0],static_cast<int>(op->Params_int[1]));
|
|
std::copy(trVals.begin(), trVals.end(), output + currentDataPoint);
|
|
currentDataPoint += op->Params_int[0];
|
|
}
|
|
else if (OperationIdentifier == OPERATION__INIT_TIME_INDEPENDENT_EVOLUTION) {
|
|
SpinOperation* op;
|
|
op = static_cast<SpinOperation*>(currOperationBytePosition[i]);
|
|
SpinsPtr->initTimeIndependentEvolution(op->Params_real[0],static_cast<int>(op->Params_int[0]));
|
|
}
|
|
else if (OperationIdentifier == OPERATION__EVOLVE_TIME_DEPENDENT) {
|
|
_do_operation_evolve_time_dependent(static_cast<SpinOperation*>(currOperationBytePosition[i]), SpinsPtr, output, currentDataPoint);
|
|
}
|
|
else {
|
|
std::cout << "Error: An unknown operation was requested... Exiting." << std::endl;
|
|
exit(1);
|
|
}
|
|
}
|
|
return currentDataPoint;
|
|
}
|
|
|
|
extern "C" long CreateSignalFromInitializedSimulator(void* SpinSystemObject, void** ops, int numOfOperations, long numOfPoints, double* output)
|
|
{
|
|
try {
|
|
if(debug)
|
|
{
|
|
std::cout<<"Number of operations: "<<numOfOperations<<std::endl;
|
|
std::cout<<"Number of points: "<<numOfPoints<<std::endl;
|
|
}
|
|
|
|
long processedNumOfPoints = ExecuteOperations(SpinSystemObject,ops,numOfOperations,output,debug);
|
|
if (debug) std::cout<<processedNumOfPoints<<std::endl;
|
|
return processedNumOfPoints;
|
|
}
|
|
catch(std::exception& ex) {
|
|
std::cerr << "An exception was thrown by Spintrum. It says: " << ex.what() << std::endl;
|
|
return 0;
|
|
}
|
|
catch(...) {
|
|
std::cerr << "An unknown exception was thrown by Spintrum. " << std::endl;
|
|
return 0;
|
|
}
|
|
}
|
|
void _fill_spin_info(int numOfSpins, double* gammas, double* jCouplings, int* spinMultiplicities, SpinInfo& spinInfo)
|
|
{
|
|
spinInfo.gammas.resize(numOfSpins,1);
|
|
spinInfo.multiplicities.resize(numOfSpins,1);
|
|
spinInfo.jCouplings.resize(numOfSpins,numOfSpins);
|
|
for(int i = 0; i < numOfSpins; i++) {
|
|
spinInfo.gammas[i] = gammas[i];
|
|
}
|
|
if(spinMultiplicities == 0) {
|
|
for(int i = 0; i < numOfSpins; i++)
|
|
spinInfo.multiplicities[i] = 2;
|
|
}
|
|
else {
|
|
for(int i = 0; i < numOfSpins; i++)
|
|
spinInfo.multiplicities[i] = spinMultiplicities[i];
|
|
}
|
|
if(debug) std::cout<<"Gammas:"<<std::endl;
|
|
if(debug) std::cout<<spinInfo.gammas<<std::endl;
|
|
devectorizeJCouplings(numOfSpins, jCouplings, spinInfo.jCouplings);
|
|
if(debug) std::cout<<"J-couplings:"<<std::endl;
|
|
if(debug) std::cout<<spinInfo.jCouplings<<std::endl;
|
|
}
|
|
|
|
extern "C" long SimulateSignal(void** ops, double* gammas, double* jCouplings, int* spinMultiplicities, int numOfSpins, int numOfOperations, long numOfPoints, double* output)
|
|
{
|
|
try {
|
|
print_software_info();
|
|
if(debug) {
|
|
std::cout<<"Number of operations: "<<numOfOperations<<std::endl;
|
|
std::cout<<"Number of points: "<<numOfPoints<<std::endl;
|
|
std::cout<<"Number of spins: "<<numOfSpins<<std::endl;
|
|
}
|
|
SpinInfo spinInfo;
|
|
_fill_spin_info(numOfSpins, gammas, jCouplings, spinMultiplicities, spinInfo);
|
|
SpinSystem<Real,std::complex<Real>,VEC,CX_VEC,MAT,CX_MAT > Spins(spinInfo.gammas,spinInfo.jCouplings,spinInfo.multiplicities);
|
|
long processedNumOfPoints = ExecuteOperations(static_cast<void*>(&Spins),ops,numOfOperations,output,debug);
|
|
if (debug) std::cout<<processedNumOfPoints<<std::endl;
|
|
return processedNumOfPoints;
|
|
}
|
|
catch(std::exception& ex) {
|
|
std::cerr << "An exception was thrown by Spintrum. It says: " << ex.what() << std::endl;
|
|
return 0;
|
|
}
|
|
catch(...) {
|
|
std::cerr << "An unknown exception was thrown by Spintrum. " << std::endl;
|
|
return 0;
|
|
}
|
|
}
|