From 3ae3ad98e1ca8c0532cb9585a4d5f51ff702a60c Mon Sep 17 00:00:00 2001 From: Samer Afach Date: Thu, 12 Jan 2017 18:21:29 +0100 Subject: [PATCH] First modified version that installs with pip --- CMakeLists.txt | 71 ++-- MANIFEST.in | 5 + run_build.sh | 32 -- run_install.sh | 20 - setup.py | 231 +++++++++- spintrum/meta.py | 2 +- spintrum/spintrum.py | 29 +- src/common/SpinInfo.cpp | 79 ++++ src/common/SpinInfo.h | 28 ++ src/common/SpinSystem.cpp | 1 + src/common/SpinSystem.h | 866 ++++++++++++++++++++++++++++++++++++++ src/common/definitions.h | 71 ++++ 12 files changed, 1335 insertions(+), 100 deletions(-) create mode 100644 MANIFEST.in delete mode 100644 run_build.sh delete mode 100644 run_install.sh create mode 100644 src/common/SpinInfo.cpp create mode 100644 src/common/SpinInfo.h create mode 100644 src/common/SpinSystem.cpp create mode 100644 src/common/SpinSystem.h create mode 100644 src/common/definitions.h diff --git a/CMakeLists.txt b/CMakeLists.txt index ecdce20..de471bb 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,20 +1,56 @@ cmake_minimum_required(VERSION 3.0) PROJECT(Spintrum) +if(NOT WIN32) + string(ASCII 27 Esc) + set(ColourReset "${Esc}[m") + set(ColourBold "${Esc}[1m") + set(Red "${Esc}[31m") + set(Green "${Esc}[32m") + set(Yellow "${Esc}[33m") + set(Blue "${Esc}[34m") + set(Magenta "${Esc}[35m") + set(Cyan "${Esc}[36m") + set(White "${Esc}[37m") + set(BoldRed "${Esc}[1;31m") + set(BoldGreen "${Esc}[1;32m") + set(BoldYellow "${Esc}[1;33m") + set(BoldBlue "${Esc}[1;34m") + set(BoldMagenta "${Esc}[1;35m") + set(BoldCyan "${Esc}[1;36m") + set(BoldWhite "${Esc}[1;37m") +endif() + +#This allows in-source build in case it's a pip installation. +#This is OK because the directory is temporary and will be removed eventually +if(NOT DEFINED PythonInstallation) + if(EXISTS "setup.py") + message(FATAL_ERROR "${BoldRed}Error: You cannot make an in-source build. Please create a sub-dir and cmake in it.${ColourReset}") + endif() +endif() + +#Get rid of warning about unused variable +set(PythonInstallation "OFF") + + if(NOT DEFINED SpintrumPath) - set(SpintrumPath ../) + set(SpintrumPath src/) endif() if(NOT DEFINED PolymathPath) - set(PolymathPath ../Polymath/) + set(PolymathPath 3rdparty/Polymath/) endif() if(NOT DEFINED OpenBLASPath) - set(OpenBLASPath ../OpenBLAS_install/) + set(OpenBLASPath spintrum/OpenBLAS_install/) +endif() +if(NOT DEFINED OptimizationLevel) + set(OptimizationLevel -Ofast) endif() set(CMAKE_BUILD_TYPE Release) #add_definitions(-DPOLYMATH_DEBUG) -set(VERSION 0.1.5) +set(VERSION 0.1.6) +execute_process(COMMAND "python3 -c 'from spintrum import meta; print(meta.__version__)'" OUTPUT_VARIABLE VERSION) add_definitions(-DSPINTRUM_VERSION="${VERSION}") INCLUDE_DIRECTORIES(${PolymathPath}/include) @@ -40,32 +76,13 @@ ADD_LIBRARY(spintrum SHARED ${PolymathSrc} ${includesList} ) -if(NOT WIN32) - string(ASCII 27 Esc) - set(ColourReset "${Esc}[m") - set(ColourBold "${Esc}[1m") - set(Red "${Esc}[31m") - set(Green "${Esc}[32m") - set(Yellow "${Esc}[33m") - set(Blue "${Esc}[34m") - set(Magenta "${Esc}[35m") - set(Cyan "${Esc}[36m") - set(White "${Esc}[37m") - set(BoldRed "${Esc}[1;31m") - set(BoldGreen "${Esc}[1;32m") - set(BoldYellow "${Esc}[1;33m") - set(BoldBlue "${Esc}[1;34m") - set(BoldMagenta "${Esc}[1;35m") - set(BoldCyan "${Esc}[1;36m") - set(BoldWhite "${Esc}[1;37m") -endif() if (UNIX) IF(CMAKE_BUILD_TYPE MATCHES Release) message("${Blue}Building in Release mode${ColourReset}") - message("${Red}Be aware that the optimization Ofast is used. If results don't make sense, change that to O3${ColourReset}") - SET( CMAKE_CXX_FLAGS_RELEASE "-Ofast") - SET( CMAKE_C_FLAGS_RELEASE "-Ofast") + message("${Green}Target optimization level: ${OptimizationLevel}${ColourReset}") + SET( CMAKE_CXX_FLAGS_RELEASE "${OptimizationLevel}") + SET( CMAKE_C_FLAGS_RELEASE "${OptimizationLevel}") ELSE() message("${BoldRed}Building in Debug mode${ColourReset}") SET( CMAKE_CXX_FLAGS_RELEASE "-g") @@ -94,7 +111,7 @@ if (UNIX) TARGET_LINK_LIBRARIES(spintrum -L"${OpenBLASPath}/lib") TARGET_LINK_LIBRARIES(spintrum -L"../${OpenBLASPath}/lib") else() - message("${BoldRed}WARNING: OpenBLAS shared library was not found in ${OpenBLASPath}. This leads to a failure in compilation. Consider pulling and compiling OpenBLAS, or alternatively change CMakeLists.txt to fit your configuration.${ColourReset}") + message("${BoldRed}WARNING: OpenBLAS shared library was not found in ${OpenBLASPath}. This may lead to a failure in compilation, unless OpenBLAS is present in the operating system.${ColourReset}") message("${BoldGreen}I will try to fix this problem by considering that your system may have OpenBLAS installed.${ColourReset}") add_definitions(-DOPENBLAS_FROM_SYSTEM) endif() diff --git a/MANIFEST.in b/MANIFEST.in new file mode 100644 index 0000000..a4df41b --- /dev/null +++ b/MANIFEST.in @@ -0,0 +1,5 @@ +recursive-include spintrum/OpenBLAS_install * +recursive-include src * +recursive-include examples * +include CMakeLists.txt + diff --git a/run_build.sh b/run_build.sh deleted file mode 100644 index fd78ef6..0000000 --- a/run_build.sh +++ /dev/null @@ -1,32 +0,0 @@ -#!/bin/bash - -#go to the directory of the script -reldir=`dirname $0` -cd $reldir -directory=`pwd` - -rm -rf dist -rm -rf build -rm -rf spintrum.egg-info/ - -mkdir build -cd build -cmake -DOpenBLASPath=../OpenBLAS_install/ .. - -if [ $? -ne 0 ]; then - echo 'CMake failed to run. Are you sure that CMake is installed?' - exit -fi - -make -if [ $? -ne 0 ]; then - echo 'make failed to run.' - exit -fi - -cd .. -python3 setup.py --verbose sdist -if [ $? -ne 0 ]; then - echo 'Python package creator installation script failed.' - exit -fi diff --git a/run_install.sh b/run_install.sh deleted file mode 100644 index 27c60ae..0000000 --- a/run_install.sh +++ /dev/null @@ -1,20 +0,0 @@ -#!/bin/bash - -#go to the directory of the script -reldir=`dirname $0` -cd $reldir -directory=`pwd` - -RED='\033[1;31m' -NC='\033[0m' # No Color - -if [ "$(id -u)" != "0" ]; then - echo -e "${RED}ERROR: This script must be run as root${NC}" 1>&2 - exit 1 -fi - -pip3 install dist/spintrum* --upgrade -if [ $? -ne 0 ]; then - echo 'Package installation failed. Did you build the package? Do you have pip3 installed?' - exit -fi diff --git a/setup.py b/setup.py index 92685a4..75a9483 100644 --- a/setup.py +++ b/setup.py @@ -1,28 +1,245 @@ from setuptools import setup, find_packages +from distutils.command.build import build as _build +from distutils.command.build_ext import build_ext as _build_ext +from setuptools import Extension from spintrum import meta import os import sys +import shutil +import errno +import tempfile +import subprocess +import glob +import fnmatch + +try: + from setuptools import setup, Extension +except ImportError: + from distutils.core import setup + from distutils.extension import Extension del os.link #solve hardlinking problem when using NTFS drive +script_dir = os.path.dirname(os.path.realpath(__file__)) +start_working_dir = os.getcwd() +temp_dir = tempfile.gettempdir() + + def _print_error(msg): sys.stderr.write("\n\nERROR: " + msg + "\n\n\n") -dirname = "spintrum" + +def which(program): + """ + Checks if a program exists, and returns its path + :param program: global program name + :return: program path if it exists, otherwise None + """ + import os + def is_exe(fpath): + return os.path.isfile(fpath) and os.access(fpath, os.X_OK) + + fpath, fname = os.path.split(program) + if fpath: + if is_exe(program): + return program + else: + for path in os.environ["PATH"].split(os.pathsep): + path = path.strip('"') + exe_file = os.path.join(path, program) + if is_exe(exe_file): + return exe_file + + return None + + +def check_programs_availability(list_of_programs): + for prog in list_of_programs: + if which(prog) is None: + _print_error("The required program " + prog + " was not found in the global scope (in PATH). Please install it. Exiting...") + sys.exit(1) + print("All required programs were found to be installed. Proceeding with building and installation.") + + +def mkdir_p(path): + try: + os.makedirs(path) + except OSError as exc: + if exc.errno == errno.EEXIST and os.path.isdir(path): + pass + else: + raise + + +def deal_with_error_code(err_code, operation_name): + if err_code != 0: + sys.stderr.write("A non-zero error code was returned at operation: " + + str(operation_name) + ". Code returned is: " + + str(err_code) + ". Exiting!\n") + sys.exit(1) + + +def inplace_change(filename, old_string, new_string): + # Safely read the input filename using 'with' + with open(filename) as f: + s = f.read() + if old_string not in s: + print('"{old_string}" not found in {filename}.'.format(**locals())) + return + + # Safely write the changed content, if found in the file + with open(filename, 'w') as f: + print('Changing "{old_string}" to "{new_string}" in {filename}'.format(**locals())) + s = s.replace(old_string, new_string) + f.write(s) + + +def clone_git_repository(dir, repos_link): + try: + shutil.rmtree(dir) + except FileNotFoundError: + pass + mkdir_p(dir) + print("Cloning repository: " + repos_link + "...") + err_code = subprocess.call(["git clone " + repos_link + " " + dir], shell=True) + deal_with_error_code(err_code, "Cloning: " + repos_link) + print("Done cloning repository: " + repos_link) + + +def ask_openblas_optimization_level(): + print("\n") + print("Which optimization option would you like to use to compile OpenBLAS?") + print("For more information on what the options mean, please visit:") + print("https://gcc.gnu.org/onlinedocs/gcc/Optimize-Options.html") + print("1) -O2") + print("2) -O3") + print("3) -Ofast") + print("4) Keep the default") + print("5) Skip compiling OpenBLAS and assume it exists in the system") + option_chosen = input("Please enter the option number: ") + if option_chosen.replace(" ", "") == '1': + return "-O2" + elif option_chosen.replace(" ", "") == '2': + return "-O3" + elif option_chosen.replace(" ", "") == '3': + return "-Ofast" + elif option_chosen.replace(" ", "") == '4': + return "" + elif option_chosen.replace(" ", "") == '5': + return "SKIPOPENBLAS" + else: + raise IndexError("Error: An unknown option was entered") + + +def ask_spintrum_optimization_level(): + print("\n") + print("Which optimization option would you like to use to compile Spintrum?") + print("For more information on what the options mean, please visit:") + print("https://gcc.gnu.org/onlinedocs/gcc/Optimize-Options.html") + print("1) -g") + print("2) -O1") + print("3) -O2") + print("4) -O3") + print("5) -Ofast") + option_chosen = input("Please enter the option number: ") + if option_chosen.replace(" ", "") == '1': + return "-g" + elif option_chosen.replace(" ", "") == '2': + return "-O1" + elif option_chosen.replace(" ", "") == '3': + return "-O2" + elif option_chosen.replace(" ", "") == '4': + return "-O3" + elif option_chosen.replace(" ", "") == '5': + return "-Ofast" + else: + raise IndexError("Error: An unknown option was entered") + + +def recursive_glob(rootdir, pattern='*'): + lst = [os.path.join(looproot, filename) + for looproot, _, filenames in os.walk(rootdir) + for filename in filenames + if fnmatch.fnmatch(filename, pattern)] + return lst + +dir_name = "spintrum" lib_file = "lib/libspintrum.so" -if not os.path.isfile(os.path.join(dirname,lib_file)): - _print_error("Unable to find the compiled library header. Please compile Spintrum before installing it.") - os._exit(1) + +class DependenciesBuilder(_build): + def run(self): + if sys.platform.startswith('win'): + _print_error("You cannot build Spintrum on Windows. It is not supported.") + + python_version = float(str(sys.version_info[0]) + "." + str(sys.version_info[1])) + if python_version < 3.3: + _print_error("You must have python version >= 3.3 to use Spintrum") + + check_programs_availability(["cmake","make","g++","gcc","gfortran","git","python3"]) + + print("Building Spintrum C/C++ extension prerequisites") + openblas_dir = "3rdparty/OpenBLAS" + openblas_dir_install = "spintrum/OpenBLAS_install" + openblas_full_path = os.path.join(start_working_dir, openblas_dir) + openblas_full_path_install = os.path.join(start_working_dir, openblas_dir_install) + openblas_git = "https://github.com/xianyi/OpenBLAS.git" + + openblas_target_optimization = ask_openblas_optimization_level() + if openblas_target_optimization != "SKIPOPENBLAS": + clone_git_repository(openblas_full_path, openblas_git) + if openblas_target_optimization != "": + if not os.path.isdir(openblas_dir): + sys.stderr.write("Cannot open expected OpenBLAS directory " + openblas_dir + "\n") + sys.exit(1) + + makefiles_to_glob = os.path.join(openblas_dir, "Makefile*") + makefiles = glob.glob(makefiles_to_glob) + + print("Number of files found: " + str(len(makefiles))) + + for f in makefiles: + print(f) + inplace_change(f, "-O2", openblas_target_optimization) + + # make openblas + os.chdir(openblas_dir) + make_openblas_err_code = subprocess.call(["make"], shell=True) + deal_with_error_code(make_openblas_err_code, "Making OpenBLAS") + + # install openblas + install_openblas_err_code = subprocess.call( + ["make PREFIX=" + openblas_full_path_install + + " install"], shell=True) + deal_with_error_code(install_openblas_err_code, "Installing OpenBLAS") + os.chdir(start_working_dir) #restore working dir + + ############## get Polymath + polymath_git = "https://git.afach.de/samerafach/Polymath" + polymath_dir = "3rdparty/Polymath" + polymath_full_path = os.path.join(start_working_dir, polymath_dir) + clone_git_repository(polymath_full_path, polymath_git) + + ############## build spintrum + spintrum_optimization_target = ask_spintrum_optimization_level() + print("Building Spintrum") + err_code = subprocess.call(["cmake . -DPythonInstallation=1 -DOptimizationLevel=" + spintrum_optimization_target], shell=True) + deal_with_error_code(err_code, "CMaking spintrum") + err_code = subprocess.call(["make"], shell=True) + deal_with_error_code(err_code, "Making spintrum") + print("Done building spintrum.") + setup(name='spintrum', version=meta.__version__, description='Software for spin systems simulation', - url='http://www.afach.de', + url='http://www.afach.de/', author='Samer Afach', author_email='samer@afach.de', license='MPL', packages=['spintrum'], - package_data={'': ['lib/libspintrum.so']}, - zip_safe=False + include_package_data=True, # enable including files with MANIFEST.in + package_data={'': [lib_file]}, + zip_safe=False, + cmdclass=dict(build=DependenciesBuilder) ) diff --git a/spintrum/meta.py b/spintrum/meta.py index 0a8da88..d3ec452 100644 --- a/spintrum/meta.py +++ b/spintrum/meta.py @@ -1 +1 @@ -__version__ = "0.1.6" +__version__ = "0.2.0" diff --git a/spintrum/spintrum.py b/spintrum/spintrum.py index c2edea1..58d73b3 100644 --- a/spintrum/spintrum.py +++ b/spintrum/spintrum.py @@ -6,7 +6,23 @@ import numpy as np import os import mpmath +def _print_error(msg_str): + sys.stderr.write(msg_str + "\n") + + +try: + _OpenBLASLibPath = os.path.join(os.path.dirname(os.path.abspath(__file__)),"OpenBLAS_install/lib/libopenblas.so.0") + c.CDLL(_OpenBLASLibPath) +except OSError: + _print_error("Warning: Failed to load OpenBLAS in its default place. Excepting the library to be available from the system.") + _SpintrumLibPath = os.path.join(os.path.dirname(os.path.abspath(__file__)),"lib/libspintrum.so") +if not os.path.exists(_SpintrumLibPath): + msg = "WARNING: Could not load the spintrum shared library in path " \ + + _SpintrumLibPath + ". This means that the installation is not valid " \ + "and will most likely not work." + _print_error(msg) + print(msg) mpmath.mp.dps = 50 @@ -41,19 +57,6 @@ def who_am_i(): return func_name - -def _print_error(msg_str): - sys.stderr.write(msg_str + "\n") - - -if not os.path.exists(_SpintrumLibPath): - msg = "WARNING: Could not load the spintrum shared library in path " \ - + _SpintrumLibPath + ". This means that the installation is not valid " \ - "and will most likely not work." - _print_error(msg) - print(msg) - - def _raise_error(msg_str): sys.stderr.write(msg_str + "\n") raise Exception(msg_str) diff --git a/src/common/SpinInfo.cpp b/src/common/SpinInfo.cpp new file mode 100644 index 0000000..7c3920f --- /dev/null +++ b/src/common/SpinInfo.cpp @@ -0,0 +1,79 @@ +#include "SpinInfo.h" +#include + + +int stringCheck(const std::string &stringToCheck) +{ + bool res = std::any_of(stringToCheck.begin(),stringToCheck.end(), + [&](char c) -> bool { + return (isdigit(c)|| + c=='+'|| + c=='-'|| + c=='.'|| + c=='\t'|| + c=='e' + );}); + return (res ? 0 : -1); +} + +SpinInfo SpinInfo::ParseJCouplingsAndGammas(std::vector jCouplingsLines, std::vector gammaLines) +{ + SpinInfo spinInfo; + SamArray parseBy; + parseBy.push_back('\t'); + parseBy.push_back(' '); + SamArray buffer; + buffer = ParseByCharacters(jCouplingsLines[0],parseBy,0); + long len = static_cast(buffer.size()); + if(static_cast(jCouplingsLines.size()) < len) + { + std::cout<<"Error: J-couplings provided do not represent a square matrix"<(jCouplingsLines[i],parseBy,0); + if(buffer.size() < len) + { + std::cout<<"Error: J-couplings provided do not represent a square matrix"<(gammaLines.size()) < len) + { + std::cout<<"Error: gammas provided do not match the side length of the J-couplings matrix"<(gammaLines[i]); + } + return spinInfo; +} diff --git a/src/common/SpinInfo.h b/src/common/SpinInfo.h new file mode 100644 index 0000000..d3e2775 --- /dev/null +++ b/src/common/SpinInfo.h @@ -0,0 +1,28 @@ +#ifndef SPININFO_H +#define SPININFO_H + +#include "definitions.h" +#include "SamArray.h" +#include "StringManipulation.h" +#include +#include +#ifdef USE_ARMA +#include +#endif +#ifdef USE_POLYMATH +#include +#endif + +int stringCheck(const std::string &stringToCheck); + +struct SpinInfo +{ + CX_MAT jCouplings; + CX_VEC gammas; + static SpinInfo ParseJCouplingsAndGammas(std::vector jCouplingsLines, std::vector gammaLines); +}; + + + + +#endif // SPININFO_H diff --git a/src/common/SpinSystem.cpp b/src/common/SpinSystem.cpp new file mode 100644 index 0000000..518b547 --- /dev/null +++ b/src/common/SpinSystem.cpp @@ -0,0 +1 @@ +#include "SpinSystem.h" diff --git a/src/common/SpinSystem.h b/src/common/SpinSystem.h new file mode 100644 index 0000000..ac3ce31 --- /dev/null +++ b/src/common/SpinSystem.h @@ -0,0 +1,866 @@ +#ifndef SpinSystem_H +#define SpinSystem_H + +#define PRINT(INFO) if(print) std::cout << INFO << std::endl +#include "definitions.h" +#include +#include +#include +#include +#include +#include +#include +#ifndef ISWIN32 +#ifdef OPENBLAS_FROM_SYSTEM +#include +#else +#include "include/cblas.h" +#endif +#else +#include +#endif +#ifdef USE_ARMA +#include +#endif +#ifdef USE_POLYMATH +#include +#endif + +#ifndef M_PI +#define M_PI 3.1415926535897932384626433832795028841971693993751 +#endif + +#define TemplateHeader template +#define TemplateTail + +TemplateHeader +class SpinSystem +{ +private: + long NumSpins; + long SideLength; + Real sampleRate; + ComplexMatrix JCouplings; + ComplexMatrix rho; + ComplexVector gamma; + G_MAT multiplicities; + ComplexMatrix evOp; + ComplexMatrix TotalTimeIndependentHamiltonian; + ComplexMatrix magnetizationMatrix; + ComplexMatrix evMat; + ComplexMatrix eigenVectors; + RealVector eigenValues; + Real magnetization; + std::vector internalTraceArray; //this stores all evolution data of this system + //for the time evolution method we use here, the diagonal is constant + //therefore, for optimized calculation of the time evolution we sum the diagonals once only + Real diagonalsSum; + + //time dependent evolution can be done only before time independent evolution for efficeincy reasons. + //This variable ensures that this is fulfilled + bool densityMatrixLock; + bool hamiltonianCalculated; + + void _constructJointSpinOperators(int direction, std::vector &output); + void _constructEffectiveSpinOperators(const ComplexVector& gamma, std::vector &output); + void _constructMagneticFieldHamiltonian(const std::vector >& SpinOperators, const ComplexVector& gamma, const Real& Bx, const Real& By, const Real& Bz, ComplexMatrix& output); + void _constructJCouplingsHamiltonian(const ComplexMatrix &JCouplings, const std::vector >& SpinOperators, ComplexMatrix &output); + void _calculate3DirsSpinOperators(std::vector >& output); + void _fillSingleMultiplicityJMatrices(int multiplicity, std::vector& matrices); + void _fillJMatrices(); + ComplexMatrix JCouplingHamiltonian; + ComplexMatrix MagneticFieldHamiltonian; + std::vector > SpinOperators; + std::vector effectiveSpinOperators; + std::vector > JMatrices; + +// void _initializeDensityMatrix(const ComplexVector& diagonalElements); + void _initializeDensityMatrix(const ComplexMatrix& densityMatrix); + + Real getOperatorExpectationValue(const ComplexMatrix& OperatorMatrix); + + void _calculateEigensystem(const ComplexMatrix& Hamiltonian, RealVector& eigenValues, ComplexMatrix& eigenVectors); + const ComplexMatrix& _calculateMagnetizationMatrix(int direction); + const ComplexMatrix& _calculateElementWiseEvolutionMatrix(const Real& sampleRate); + template + std::vector _vectorizeMatrix_upperTriangular_columnMajor(const G_MAT& mat); + int directionToIndex(char dir); + +public: + SpinSystem(); + SpinSystem(const ComplexVector& Gammas, const ComplexMatrix& JCouplings); + SpinSystem(const ComplexVector& Gammas, const ComplexMatrix& JCouplings, const G_MAT& SpinMultiplicities); + void init(const ComplexVector& Gammas, const ComplexMatrix& JCouplings); + void init(const ComplexVector& Gammas, const ComplexMatrix& JCouplings, const G_MAT& SpinMultiplicities); + void tipSpins(int direction, Real FieldVsTimeArea); + const long& getNumOfSpins(); + const long& getSideLength(); + ComplexMatrix& getDensityMatrix(); + Real getSpinExpectation(int direction); + const ComplexMatrix& getSpinEffectiveOperator(int direction); +// const std::vector >& getEvolutionOperator(Real samplingRate); + void calculateHamiltonians(const Real& Bx, const Real& By, const Real& Bz); + void initTimeIndependentEvolution(const Real& sampleRate, int MeasurementDirection); + void setJCouplings(const ComplexMatrix& JCoupings); + void setSpinMultiplicities(const G_MAT& spinMultiplicities); + void setGammas(const ComplexVector& gammas); + std::vector evolveDensityMatrix(const Real& sampleRate, + const std::vector& variableNames, + const CX_MAT& variableData, + int threads = 0, + const std::set& store_trace_directions = std::set()); + std::vector solveMagnetization(const long numOfPoints); + ~SpinSystem(); + bool print = 1; + const Real& getMagnetization(); + std::vector solveMagnetization_parallel(const long numOfPoints, int threads); + void OptimizeEvMatAndMagnetizationMat(); + void populate_thermal(const Real& Bx, const Real& By, const Real& Bz, const Real& Temperature = 293.778); + const std::vector &getTraceArray(); +}; + +TemplateHeader +SpinSystem TemplateTail::SpinSystem() +{ +} + +TemplateHeader +SpinSystem TemplateTail::SpinSystem(const ComplexVector& Gammas, const ComplexMatrix& JCouplings) +{ + this->init(Gammas,JCouplings); +} + +TemplateHeader +SpinSystem TemplateTail::SpinSystem(const ComplexVector& Gammas, const ComplexMatrix& JCouplings, const G_MAT& SpinMultiplicities) +{ + this->init(Gammas,JCouplings,SpinMultiplicities); +} + +TemplateHeader +SpinSystem TemplateTail::~SpinSystem() +{ + +} + +TemplateHeader +const Real& +SpinSystem TemplateTail::getMagnetization() +{ + magnetization = std::real(std::accumulate(magnetizationMatrix.begin(),magnetizationMatrix.end(),Complex(0),std::plus())); + return magnetization; +} + +TemplateHeader +const ComplexMatrix& +SpinSystem TemplateTail::_calculateMagnetizationMatrix(int direction) +{ + ComplexMatrix transposed; + CONJTRANSPOSE(eigenVectors,transposed); + ComplexMatrix tmp; + ELEMENTWISEPRODUCT((transposed*rho*eigenVectors),(transposed*effectiveSpinOperators[direction]*eigenVectors),tmp); + +// magnetizationMatrix = arma::conv_to< std::vector > >::from(vectorise(tmp)); + //above line sums all matrix elements + //The below code sums the upper triangular and considers the diagonals constant to be summed only once for every point + Complex tr; + MATTRACE(tmp,tr); + diagonalsSum = std::real(tr); + HALFVECTORIZE(tmp,magnetizationMatrix); + return magnetizationMatrix; +} + +TemplateHeader +const ComplexMatrix& +SpinSystem TemplateTail::_calculateElementWiseEvolutionMatrix(const Real &sampleRate) +{ + ComplexMatrix evMatTemp; + evMatTemp.resize(SideLength,SideLength); + for(long i = 0; i < SideLength; i++) + { + evMatTemp.at(i,i) = 1; + for(long j = i+1; j < SideLength; j++) + { + evMatTemp.at(i,j) = std::exp(-Complex(0,1)*(eigenValues[i]-eigenValues[j])/sampleRate); + evMatTemp.at(j,i) = std::conj(evMatTemp.at(i,j)); + } + } +// evMat = arma::conv_to< std::vector > >::from(vectorise(evMatTemp)); + //above line sums all matrix elements + //The below code sums the upper triangular and considers the diagonals constant to be summed only once for every point + HALFVECTORIZE(evMatTemp,evMat); +// std::cout< SpinMultiplicities(Gammas.rows(),1,2); + this->init(Gammas,JCouplings,SpinMultiplicities); +} + +TemplateHeader +void +SpinSystem TemplateTail::init(const ComplexVector& Gammas, + const ComplexMatrix& JCouplings, + const G_MAT& SpinMultiplicities) +{ + densityMatrixLock = 0; + hamiltonianCalculated = 0; + NumSpins = Gammas.size(); + setGammas(Gammas); + setSpinMultiplicities(SpinMultiplicities); + setJCouplings(JCouplings); + //sidelength is the product of all matrix sizes, hence all multiplicities + SideLength = std::accumulate(multiplicities.begin(),multiplicities.end(),1,std::multiplies()); + if(print) PRINT("Constructing Pauli matrices"); + _fillJMatrices(); + if(print) PRINT("Constructing spin operators"); + _calculate3DirsSpinOperators(SpinOperators); + _constructEffectiveSpinOperators(gamma,effectiveSpinOperators); + if(print) PRINT("Done initalizing the system."); +} + +//TemplateHeader +//void SpinSystem TemplateTail::_initializeDensityMatrix(const ComplexVector &diagonalElements) +//{ +// for(long i = 0; i < SideLength; i++) +// { +// rho.at(i,i) = diagonalElements[i]; +// } +//} + +TemplateHeader +void +SpinSystem TemplateTail::_initializeDensityMatrix(const ComplexMatrix &densityMatrix) +{ + rho = densityMatrix; +} + +TemplateHeader +void +SpinSystem TemplateTail::_calculateEigensystem(const ComplexMatrix &Hamiltonian, + RealVector& eigenValues, + ComplexMatrix& eigenVectors) +{ +// eig_sym(eigenValues, eigenVectors, Hamiltonian, "std"); + EIGSYM(eigenValues,eigenVectors,Hamiltonian); +} + +TemplateHeader +void +SpinSystem TemplateTail::calculateHamiltonians(const Real &Bx, + const Real &By, + const Real &Bz) +{ + if(print) PRINT("Constructing J-couplings Hamiltonian"); + _constructJCouplingsHamiltonian(JCouplings,SpinOperators,JCouplingHamiltonian); + if(print) PRINT("Constructing the magnetic field Hamiltonian"); + _constructMagneticFieldHamiltonian(SpinOperators,gamma,Bx,By,Bz,MagneticFieldHamiltonian); + TotalTimeIndependentHamiltonian = MagneticFieldHamiltonian; + if(JCouplings.rows() != 0 && JCouplings.columns() != 0) + TotalTimeIndependentHamiltonian += JCouplingHamiltonian; + hamiltonianCalculated = 1; + if(print) PRINT("Done preparing the Hamiltonian"); +} + +TemplateHeader +void +SpinSystem TemplateTail::initTimeIndependentEvolution(const Real &SampleRate, + int MeasurementDirection) +{ + if(!hamiltonianCalculated) + { + throw std::logic_error("The time-independent Hamiltonian must be set before initializing the time-independent evolution."); + } + this->sampleRate = SampleRate; + if(print) PRINT("Solving eigen problem"); + _calculateEigensystem(TotalTimeIndependentHamiltonian,eigenValues,eigenVectors); + if(print) PRINT("Calculating evolution and magnetization matrices"); + if(MeasurementDirection < 0 || MeasurementDirection > 2) { + throw std::logic_error("Measurement direction can only be 0 for x, 1 for y, 2 for z"); + } + _calculateMagnetizationMatrix(MeasurementDirection); + _calculateElementWiseEvolutionMatrix(sampleRate); + densityMatrixLock = 1; + if(print) PRINT("Done initializing time independent evolution."); +} + +TemplateHeader +void +SpinSystem TemplateTail::setJCouplings(const ComplexMatrix& J_Couplings) +{ + this->JCouplings = J_Couplings; +// std::cout<<"Parameters in SpinSystem: "<< std::setprecision(15)< &spinMultiplicities) +{ + this->multiplicities = spinMultiplicities; + if(multiplicities.rows() != gamma.rows()) + throw std::length_error("Multiplicities and gyromagnetic ratios must have the same length."); +} + +TemplateHeader +void +SpinSystem TemplateTail::setGammas(const ComplexVector &gammas) +{ + this->gamma = gammas; +} + +TemplateHeader +int +SpinSystem TemplateTail::directionToIndex(char dir) +{ + switch (dir) { + case 'x': + return 0; + break; + case 'y': + return 1; + break; + case 'z': + return 2; + break; + default: + throw std::logic_error("Directions, for which trace has to be stored must be given as a string; eg x or xyz or xz, etc."); + break; + } +} + +TemplateHeader +std::vector +SpinSystem TemplateTail::evolveDensityMatrix( + const Real &sampleRate, + const std::vector &variableNames, + const CX_MAT &variableData, + int threads, + const std::set& store_trace_directions) +{ + if(densityMatrixLock) { + throw std::logic_error("You cannot use time dependent evolution matrix after initializing the time independent evolution."); + } + if(!hamiltonianCalculated) { + throw std::logic_error("The time-independent Hamiltonian must be set before evolving your density matrix."); + } + if(variableData.columns() != variableNames.size()) { + throw std::logic_error("Number of data columns should be equal to variable names."); + } + if(std::any_of(store_trace_directions.begin(),store_trace_directions.end(),[](char c) {return (c != 'x' && c != 'y' && c != 'z');})) { + throw std::logic_error("Directions, for which trace has to be stored must be given as a string; eg x or xyz or xz, etc."); + } + if(threads > 0) + { + openblas_set_num_threads(threads); + } + else + { + int n = static_cast(std::thread::hardware_concurrency()); + if(n > 0) + openblas_set_num_threads(n); + else + throw std::runtime_error("Unable to automatically detect the number of threads. Please specify that manually."); + } + ComplexMatrix evOp(SideLength,SideLength); + ComplexMatrix evOpT(SideLength,SideLength); + ComplexMatrix resMat(SideLength,SideLength); + ComplexMatrix hamiltonian(SideLength,SideLength); + Real dt = static_cast(1)/sampleRate; + long startPoint = internalTraceArray.size(); + std::vector TraceArray; + if(!store_trace_directions.empty()) { + //current size + number_of_directions_to_measure*number_of_points + this->internalTraceArray.resize(internalTraceArray.size() + variableData.rows()*store_trace_directions.size()); + TraceArray.resize(variableData.rows()*store_trace_directions.size()); + } + + for(long p = 0; p < variableData.rows(); p++) { + hamiltonian.ZEROS; + try { + for(long i = 0; i < (long)variableNames.size(); i++) { + if(variableNames[i] == std::string("Bx")) { + hamiltonian -= ((variableData.at(p,i))*effectiveSpinOperators[0]); + } + else if(variableNames[i] == std::string("By")) { + hamiltonian -= ((variableData.at(p,i))*effectiveSpinOperators[1]); + } + else if(variableNames[i] == std::string("Bz")) { + hamiltonian -= ((variableData.at(p,i))*effectiveSpinOperators[2]); + } + else { + throw std::runtime_error("Unknown variable to be evolved: " + variableNames[i]); + } + } + } + catch(std::exception& ex) { + throw std::runtime_error("Exception thrown while performing time dependent evolution. It says: " + std::string(ex.what())); + } +// hamiltonian *= (2.*M_PI); //not necessary as effectiveSpinOperators[] already contain the 2*M_PI factor. This was tested. + hamiltonian += TotalTimeIndependentHamiltonian; + hamiltonian *= Complex(0,dt); + EXPMAT_AH(hamiltonian,evOp); + CONJTRANSPOSE(evOp,evOpT); + Poly::MultiplyMatrices(rho,evOp,resMat,Poly::MATRIX_HERMITIAN,Poly::MATRIX_GENERAL); + Poly::MultiplyMatrices(evOpT,resMat,rho,Poly::MATRIX_GENERAL,Poly::MATRIX_GENERAL); +// rho = evOpT*rho*evOp; //this like is equivalent to the 2 lines above + if(!store_trace_directions.empty()) + { + Complex tr = 0; + int axis = 0; + for(auto it = store_trace_directions.begin(); it != store_trace_directions.end(); ++it) + { + Poly::MultiplyMatrices(rho,effectiveSpinOperators[directionToIndex(*it)],resMat,Poly::MATRIX_HERMITIAN,Poly::MATRIX_HERMITIAN); + MATTRACE(resMat,tr); + TraceArray[p*store_trace_directions.size()+axis] = std::real(tr); + axis++; + } + } + } + std::copy(TraceArray.begin(),TraceArray.end(),internalTraceArray.begin()+startPoint); + return TraceArray; +} + +TemplateHeader +template +std::vector SpinSystem TemplateTail::_vectorizeMatrix_upperTriangular_columnMajor(const G_MAT& mat) +{ + std::vector output((COLS(mat)*(COLS(mat)-1))/2); + typename G_MAT::const_iterator it = mat.begin() + ROWS(mat); //iterator at rows to skip in every step, starts at second column + long toSkipInVec = 0; + for(int i = 1; i < COLS(mat); i++) //Starts with 1 to skip the diagonal + { + std::copy(it, it + i, output.begin() + toSkipInVec); + toSkipInVec += i; + it += ROWS(mat); + } + return output; +} + +TemplateHeader +std::vector +SpinSystem TemplateTail::solveMagnetization(const long numOfPoints) +{ + std::vector TraceArray(numOfPoints); + if(print) + { + std::cout<<"Evolving time until t="< const c) -> Real& { return acc += c.real(); }); + + TraceArray[j] = magnetization; + std::transform( magnetizationMatrix.begin(), magnetizationMatrix.end(), + evMat.begin(), magnetizationMatrix.begin(), + std::multiplies >()); +// magnetization = real(accu(magnetizationMatrix)); +// traceArray[j] = magnetization; +// ELEMENTWISEPRODUCT(magnetizationMatrix,evMat,magnetizationMatrix); + } + if(print) + { + std::cout.precision(10); + std::cout<internalTraceArray.resize(internalTraceArray.size() + numOfPoints); + return TraceArray; +} + +TemplateHeader +std::vector +SpinSystem TemplateTail::solveMagnetization_parallel(const long numOfPoints, + int threads) +{ +// Real dt = static_cast(1)/sampleRate; + std::vector TraceArray(numOfPoints); + std::vector > > par_evMat(threads); + std::vector > > par_magnetizationMatrix(threads); + for(long i = 0; i < threads; i++) + { + par_evMat[i].resize (evMat.size() / threads + (i < static_cast(evMat.size() % threads) ? 1 : 0)); + par_magnetizationMatrix[i].resize(evMat.size() / threads + (i < static_cast(evMat.size() % threads) ? 1 : 0)); + } + std::vector sizes(threads); + for(long i = 0; i < threads; i++) + { + sizes[i] = par_evMat[i].size(); + } + for(long i = 1; i < (long)sizes.size(); i++) + { + sizes[i] += sizes[i-1]; + } + int index = 0; + for(long i = 0; i < static_cast(evMat.size()); i++) + { + if(i >= sizes[index]) + { + ++index; + } + par_evMat[index][i - (index > 0 ? sizes[index - 1] : 0)] = evMat[i]; + par_magnetizationMatrix[index][i - (index > 0 ? sizes[index - 1] : 0)] = magnetizationMatrix[i]; + } + std::vector > > par_traceArray(numOfPoints); + for(long i = 0; i < (long)par_traceArray.size(); i++) + { + par_traceArray[i].resize(threads); + } + if(print) + { + std::cout<<"Evolving time until t="<(std::ceil(std::log10(sampleRate)))); + std::cout<(numOfPoints)/sampleRate; + std::cout.flush(); + } + } +// par_traceArray[i] = accu(par_magnetizationMatrix[i]); +// par_magnetizationMatrix[i] = par_magnetizationMatrix[i] % par_evMat[i]; + par_traceArray[j][i] = std::accumulate(par_magnetizationMatrix[i].begin(), par_magnetizationMatrix[i].end(), std::complex(0,0)); +// par_traceArray[j][i] = std::accumulate(par_magnetizationMatrix[i].cbegin(), par_magnetizationMatrix[i].cend(), Real{}, [](Real const acc, std::complex const& c) { return acc + c.real(); }); + + std::transform( par_magnetizationMatrix[i].begin(), par_magnetizationMatrix[i].end(), + par_evMat[i].begin(), par_magnetizationMatrix[i].begin(), + std::multiplies >()); + } + } +#pragma omp parallel for schedule(dynamic,1) num_threads(threads) + for(long i = 0; i < (long)par_traceArray.size(); i++) + { + TraceArray[i] = diagonalsSum + 2*real(std::accumulate(par_traceArray[i].begin(), par_traceArray[i].end(), std::complex(0,0))); +// traceArray[i] = 2*std::accumulate(par_traceArray[i].cbegin(), par_traceArray[i].cend(), Real{}, [](Real const acc, std::complex const& c) { return acc + c.real(); }); + } + if(print) + { + std::cout<<"\r"; + std::cout<<"t = "<< numOfPoints/sampleRate << ":" << static_cast(numOfPoints)/sampleRate; + std::cout.flush(); + + std::cout.precision(10); + std::cout<internalTraceArray.resize(internalTraceArray.size() + numOfPoints); + return TraceArray; +} + +TemplateHeader +void +SpinSystem TemplateTail::OptimizeEvMatAndMagnetizationMat() +{ + for(long i = 0; i < magnetizationMatrix.size(); i++) + { + if((std::abs(magnetizationMatrix.at(i))) <= 1.e-15) + { + evMat.erase(evMat.begin() + i); + magnetizationMatrix.erase(evMat.begin() + i); + } + } +} + +TemplateHeader +void +SpinSystem TemplateTail::populate_thermal(const Real &Bx, + const Real &By, + const Real &Bz, + const Real &Temperature) +{ + if (print) PRINT("Populating the density matrix in thermal equilibrium"); + ComplexMatrix Z(SideLength,SideLength); + Real B[3] = {Bx, By, Bz}; + Z.ZEROS; + const Real k = (1.38064852e-23); + const Real hbar = (1.054571800e-34); + for(long dir = 0; dir < 3; dir++) + { + for(long i = 0; i < (long)SpinOperators[dir].size(); i++) + { + Z += (2*M_PI*hbar*B[dir]*gamma[i]/(k*Temperature))*SpinOperators[dir][i]; + } + } + ComplexMatrix Zexp; + EXPMAT_H(Z,Zexp); + Complex traceVal; + MATTRACE(Zexp,traceVal); + rho = Zexp*static_cast(1)/traceVal; + densityMatrixLock = 0; +// rho.ZEROS; +// rho[0] = 1; +// std::cout< &SpinSystem TemplateTail::getTraceArray() +{ + return internalTraceArray; +} + +TemplateHeader +ComplexMatrix& +SpinSystem TemplateTail::getDensityMatrix() +{ + if(densityMatrixLock) { + throw std::logic_error("You cannot get the density matrix after initializing time independent evolution."); + } + return rho; +} + +TemplateHeader +Real +SpinSystem TemplateTail::getSpinExpectation(int direction) +{ + if(densityMatrixLock) { + throw std::logic_error("You cannot modify the density matrix after initializing time independent evolution."); + } + return real(trace(rho*effectiveSpinOperators[direction])); +} + +TemplateHeader +const ComplexMatrix& +SpinSystem TemplateTail::getSpinEffectiveOperator(int direction) +{ + if(densityMatrixLock) { + throw std::logic_error("You cannot modify the density matrix after initializing time independent evolution."); + } + return effectiveSpinOperators[direction]; +} + +//TemplateHeader +//const std::vector > &SpinSystem TemplateTail::getEvolutionOperator(Real samplingRate) +//{ +// evOp = arma::conv_to > >::from(vectorise(expmat(std::complex(0,-1)*TotalTimeIndependentHamiltonian*1/samplingRate))); +// return evOp; +//} + +/** + * @brief SpinSystem::tipSpins + * @param direction + * @param FieldVsTimeArea + * + * Tips the spins instantaneously around the axis (direction) chose. Direction can be x=0, y=1, z=2. The amount of tipping is determined + * basic magnetic resonance, where the area under magnetic field vs Time determines the tipping angle. + */ + +TemplateHeader +void +SpinSystem TemplateTail::tipSpins(int direction, + Real FieldVsTimeArea) +{ + if(densityMatrixLock) { + throw std::logic_error("You cannot modify the density matrix after initializing time independent evolution (error from tipSpins())."); + } + ComplexMatrix initialStateEvolution; + EXPMAT_AH(Complex(0,1)*FieldVsTimeArea*effectiveSpinOperators[direction],initialStateEvolution); + ComplexMatrix inv; + CONJTRANSPOSE(initialStateEvolution,inv); + rho = inv*rho*initialStateEvolution; +} + +TemplateHeader +const long& +SpinSystem TemplateTail::getNumOfSpins() +{ + return NumSpins; +} + +TemplateHeader +const long& +SpinSystem TemplateTail::getSideLength() +{ + return SideLength; +} + +TemplateHeader +void +SpinSystem TemplateTail::_fillSingleMultiplicityJMatrices(int multiplicity, + std::vector& matrices) +{ + matrices.resize(3); + const Real SpinQ = (static_cast(multiplicity)-static_cast(1))/static_cast(2); + MAT Jp(multiplicity,multiplicity), Jm(multiplicity,multiplicity); //J+ and J- operators + for(long i = 0; i < 3; i++) + { + matrices[i].resize(multiplicity,multiplicity,Complex(0)); + } + Real SpinValue_forLoop = SpinQ; + for(long i = 0; i < multiplicity; i++) + { + matrices[2].at(i,i) = SpinValue_forLoop; + SpinValue_forLoop -= Real(1); + } + //Build J+ using off-diagonal elements being sqrt((J - m) (J + m + 1)) = sqrt(j(j+1)-m(m+1)) + SpinValue_forLoop = SpinQ-1; + for(long m = 0; m < multiplicity-1; m++) + { + Jp.at(m,m+1) = std::sqrt((SpinQ-SpinValue_forLoop)*(SpinQ+SpinValue_forLoop+1)); + SpinValue_forLoop -= Real(1); + } + //Build J- using off-diagonal elements being sqrt((J + m) (J - m + 1)) = sqrt(j(j+1)-m(m-1)) + SpinValue_forLoop = SpinQ-1; + for(long m = 1; m < multiplicity; m++) + { + Jm.at(m,m-1) = std::sqrt((SpinQ-SpinValue_forLoop)*(SpinQ+SpinValue_forLoop+1)); + SpinValue_forLoop -= Real(1); + } + matrices[0] = Complex(0.5)*(Jp+Jm); + matrices[1] = Complex(0,-0.5)*(Jp-Jm); +} + +TemplateHeader +void +SpinSystem TemplateTail::_fillJMatrices() +{ + std::set uniqueMultiplicities(multiplicities.begin(),multiplicities.end()); + int maxMultiplicity = *std::max_element(uniqueMultiplicities.begin(),uniqueMultiplicities.end()); + //JMatrices are stored by multiplicity, so multiplicity 2 is element number 0; 3 is element number 1, etc... + JMatrices.resize(maxMultiplicity-1); + for(auto it = uniqueMultiplicities.begin(); it != uniqueMultiplicities.end(); ++it) + { + _fillSingleMultiplicityJMatrices(*it,JMatrices[*it - 2]); + } +} + +/** + A function that creates spin operators for each spin in the system in the joint Hilbert space using kronecker product of + unity matrices and Pauli matrices + + direction is the spin operator direction, x=0, y=1, z=2 + */ +TemplateHeader +void +SpinSystem TemplateTail::_constructJointSpinOperators(int direction, + std::vector &output) +{ + output.resize(NumSpins); + ComplexMatrix mat; + mat. ZEROS; +// ComplexMatrix identity(ROWS(JMatrices[direction]),COLS(JMatrices[direction])); + ComplexMatrix initMat(1,1); + initMat.at(0,0) = 1; +// for(long i = 0; i < (long)ROWS(identity); i++) +// { +// identity.at(i,i) = 1; +// } + for(long i = NumSpins - 1; i >= 0; i--) + { + ComplexMatrix res = initMat; + for(long j = NumSpins - 1; j >= 0; j--) + { + if(i == j) + { + //J-Matrix of multiplicity S is stored in JMatrices[S-2] + mat = JMatrices[multiplicities[j]-2][direction]; + } + else + { + mat.resize(multiplicities[j],multiplicities[j]); + mat.ZEROS; + for(long k = 0; k < mat.rows(); k++) + mat.at(k,k) = 1; + } + KRONPROD(res,mat,res); + } + output[i] = res; + } +} + +TemplateHeader +void +SpinSystem TemplateTail::_constructEffectiveSpinOperators(const ComplexVector &gamma, + std::vector &output) +{ + output.resize(3); + for(long i = 0; i < 3; i++) + { + output[i].resize(SideLength,SideLength); + for(long j = 0; j < NumSpins; j++) //2 denotes the [x,y,z] component of combined spin operators + { + output[i] += (2*M_PI)*gamma[j]*SpinOperators[i][j]; //according to Ledbetter et al. 2009 Eq. 2 + } + } +} + +TemplateHeader +void +SpinSystem TemplateTail::_constructMagneticFieldHamiltonian(const std::vector >& SpinOperators, + const ComplexVector& gamma, + const Real& Bx, + const Real& By, + const Real& Bz, + ComplexMatrix& output) +{ + output.resize(SideLength,SideLength); + output.ZEROS; + //removed hbar to save computation as it will be removed also from the evolution operator + output = ((-Bx)*effectiveSpinOperators[0]+(-By)*effectiveSpinOperators[1]+(-Bz)*effectiveSpinOperators[2]); //-mu.B +} + +TemplateHeader +void +SpinSystem TemplateTail::_constructJCouplingsHamiltonian(const ComplexMatrix &JCouplings, + const std::vector >& SpinOperators, + ComplexMatrix &output) +{ + if(JCouplings.rows() == 0 || JCouplings.columns() == 0) + return; + output.resize(SideLength,SideLength); + output.ZEROS; + for(long i = 0; i < NumSpins; i++) + { + for(long j = i+1; j < NumSpins; j++) + { + //removed hbar to save computation as it will be removed also from the evolution operator + output += ((2*M_PI)*JCouplings.AT(i,j))*(SpinOperators[0][i]*SpinOperators[0][j] + SpinOperators[1][i]*SpinOperators[1][j] + SpinOperators[2][i]*SpinOperators[2][j]); + } + } +} + +TemplateHeader +void +SpinSystem TemplateTail::_calculate3DirsSpinOperators(std::vector > &output) +{ + output.resize(3); //combined spin operators for 3 directions + for(long i = 0; i < 3; i++) + { + _constructJointSpinOperators(i,output[i]); + } +} + +#undef TemplateHeader +#undef TemplateTail + +#endif // SpinSystem_H diff --git a/src/common/definitions.h b/src/common/definitions.h new file mode 100644 index 0000000..14f070a --- /dev/null +++ b/src/common/definitions.h @@ -0,0 +1,71 @@ +#ifndef DEFINITIONS_H +#define DEFINITIONS_H + +//#define USE_ARMA +#define USE_POLYMATH + +//#define ARMA_NO_DEBUG + +#include + +#ifndef REPLACEMENTS +#define REPLACEMENTS + +#define REALTYPE_DOUBLE +//#define REALTYPE_LONGDOUBLE + + +#ifdef REALTYPE_DOUBLE + typedef double _Real; +#endif +#ifdef REALTYPE_LONGDOUBLE + typedef long double _Real; +#endif + +#ifdef USE_ARMA +#define VEC arma::vec +#define CX_VEC arma::cx_vec +#define MAT arma::Mat<_Real> +#define CX_MAT arma::Mat> +#define G_MAT arma::Mat +#define ZEROS zeros() +#define AT(iIndex,jIndex) at(iIndex,jIndex) +#define EXPMAT_H(in,out) out = expmat(in) +#define EXPMAT_AH(in,out) out = expmat(in) +#define MATINVERSE(in,out) out = in.i() +#define ELEMENTWISEPRODUCT(in1,in2,out) out = in1 % in2 +#define EIGSYM(vals,vecs,mat) eig_sym(vals, vecs, mat, "std"); +#define MATTRACE(in,out) out = trace(in) +#define CONJTRANSPOSE(in,out) out = in.t() +#define TRANSPOSE(in,out) out = in.st() +#define HALFVECTORIZE(in,out) out = _vectorizeMatrix_upperTriangular_columnMajor(in) +#define KRONPROD(in1,in2,out) out = kron(in1,in2) +#define ROWS(mat) mat.n_rows +#define COLS(mat) mat.n_cols +#endif + +#ifdef USE_POLYMATH +#define VEC Poly::Matrix<_Real> +#define CX_VEC Poly::Matrix> +#define MAT Poly::Matrix<_Real> +#define CX_MAT Poly::Matrix> +#define G_MAT Poly::Matrix +#define ZEROS zeros() +#define AT(iIndex,jIndex) at(iIndex,jIndex) +#define EXPMAT_H(in,out) out = MatrixExp(in,Poly::MATRIX_HERMITIAN) +#define EXPMAT_AH(in,out) out = MatrixExp(in,Poly::MATRIX_ANTIHERMITIAN) +#define MATINVERSE(in,out) out = in.inverse() +#define EIGSYM(vals,vecs,mat) EigenVecsVals(vals, vecs, mat, Poly::MATRIX_HERMITIAN, Poly::MATRIX_HERMITIAN); +#define ELEMENTWISEPRODUCT(in1,in2,out) out = in1.elementWiseProduct(in2) +#define MATTRACE(in,out) out = Trace(in) +#define CONJTRANSPOSE(in,out) out = Poly::ConjugateTranspose(in) +#define TRANSPOSE(in,out) out = Poly::Transpose(in) +#define HALFVECTORIZE(in,out) out = in.vectorize(Poly::VECMODE_UPPER_TRIANGULAR_NO_DIAGONAL) +#define KRONPROD(in1,in2,out) out = KroneckerProduct(in1,in2) +#define ROWS(mat) mat.rows() +#define COLS(mat) mat.columns() +#endif + +#endif +#endif // DEFINITIONS_H +