First modified version that installs with pip

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}")
+        message("${Green}Target optimization level: ${OptimizationLevel}${ColourReset}")
-        SET( CMAKE_CXX_FLAGS_RELEASE "-Ofast")
+        SET( CMAKE_CXX_FLAGS_RELEASE "${OptimizationLevel}")
-        SET( CMAKE_C_FLAGS_RELEASE "-Ofast")
+        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()

MANIFEST.in

@@ -0,0 +1,5 @@
+recursive-include spintrum/OpenBLAS_install *
+recursive-include src *
+recursive-include examples *
+include CMakeLists.txt
+

run_build.sh

@@ -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

run_install.sh

@@ -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

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.")
+class DependenciesBuilder(_build):
-    os._exit(1)
+    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']},
+      include_package_data=True, # enable including files with MANIFEST.in
-      zip_safe=False
+      package_data={'': [lib_file]},
+      zip_safe=False,
+      cmdclass=dict(build=DependenciesBuilder)
       )

spintrum/meta.py

@@ -1 +1 @@
-__version__ = "0.1.6"
+__version__ = "0.2.0"

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)

src/common/SpinInfo.cpp

@@ -0,0 +1,79 @@
+#include "SpinInfo.h"
+#include <algorithm>
+
+
+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<std::string> jCouplingsLines, std::vector<std::string> gammaLines)
+{
+    SpinInfo spinInfo;
+    SamArray<char> parseBy;
+    parseBy.push_back('\t');
+    parseBy.push_back(' ');
+    SamArray<double> buffer;
+    buffer = ParseByCharacters<double>(jCouplingsLines[0],parseBy,0);
+    long len = static_cast<long>(buffer.size());
+    if(static_cast<long>(jCouplingsLines.size()) < len)
+    {
+        std::cout<<"Error: J-couplings provided do not represent a square matrix"<<std::endl;
+        std::exit(2);
+    }
+    spinInfo.jCouplings.resize(len,len);
+    spinInfo.gammas.resize(len,1);
+    spinInfo.jCouplings. ZEROS;
+    spinInfo.gammas. ZEROS;
+    for(long i = 0; i < len; i++)
+    {
+        //function to delete newline character from windows files
+        jCouplingsLines[i].erase(std::remove(jCouplingsLines[i].begin(), jCouplingsLines[i].end(), '\r'), jCouplingsLines[i].end());
+
+        if(stringCheck(jCouplingsLines[i]) < 0)
+        {
+            throw std::runtime_error("Error in J-Coupling file: only 1,2,3,4,5,6,7,8,9,+,-,. and tabs are allowed as symbols in the file");
+        }
+
+        buffer = ParseByCharacters<double>(jCouplingsLines[i],parseBy,0);
+        if(buffer.size() < len)
+        {
+            std::cout<<"Error: J-couplings provided do not represent a square matrix"<<std::endl;
+            std::exit(3);
+        }
+        for(long j = 0; j < len; j++)
+        {
+
+            spinInfo.jCouplings.AT(i,j) = buffer[j];
+        }
+    }
+
+
+    if(static_cast<long>(gammaLines.size()) < len)
+    {
+        std::cout<<"Error: gammas provided do not match the side length of the J-couplings matrix"<<std::endl;
+        std::exit(3);
+    }
+    for(long i = 0; i < len; i++)
+    {
+        //function to delete newline character from windows files
+        gammaLines[i].erase(std::remove(gammaLines[i].begin(), gammaLines[i].end(), '\r'), gammaLines[i].end());
+
+        if(stringCheck(gammaLines[i]) < 0)
+        {
+            throw std::runtime_error("Error in Gamma file: only 1,2,3,4,5,6,7,8,9,+,-,. and tabs are allowed as symbols in the file");
+        }
+
+        spinInfo.gammas[i] = FromString<double>(gammaLines[i]);
+    }
+    return spinInfo;
+}

src/common/SpinInfo.h

@@ -0,0 +1,28 @@
+#ifndef SPININFO_H
+#define SPININFO_H
+
+#include "definitions.h"
+#include "SamArray.h"
+#include "StringManipulation.h"
+#include <string>
+#include <vector>
+#ifdef USE_ARMA
+#include <armadillo>
+#endif
+#ifdef USE_POLYMATH
+#include <Polymath.h>
+#endif
+
+int stringCheck(const std::string &stringToCheck);
+
+struct SpinInfo
+{
+    CX_MAT jCouplings;
+    CX_VEC gammas;
+    static SpinInfo ParseJCouplingsAndGammas(std::vector<std::string> jCouplingsLines, std::vector<std::string> gammaLines);
+};
+
+
+
+
+#endif // SPININFO_H

src/common/SpinSystem.cpp

@@ -0,0 +1 @@
+#include "SpinSystem.h"

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 <complex>
+#include <algorithm>
+#include <iomanip>
+#include <set>
+#include <vector>
+#include <thread>
+#include <numeric>
+#ifndef ISWIN32
+#ifdef OPENBLAS_FROM_SYSTEM
+#include <openblas/cblas.h>
+#else
+#include "include/cblas.h"
+#endif
+#else
+#include <cblas.h>
+#endif
+#ifdef USE_ARMA
+#include <armadillo>
+#endif
+#ifdef USE_POLYMATH
+#include <Polymath.h>
+#endif
+
+#ifndef M_PI
+#define M_PI 3.1415926535897932384626433832795028841971693993751
+#endif
+
+#define TemplateHeader template <typename Real, typename Complex, typename RealVector, typename ComplexVector, typename RealMatrix, typename ComplexMatrix>
+#define TemplateTail <Real, Complex, RealVector, ComplexVector, RealMatrix, ComplexMatrix>
+
+TemplateHeader
+class SpinSystem
+{
+private:
+    long NumSpins;
+    long SideLength;
+    Real sampleRate;
+    ComplexMatrix JCouplings;
+    ComplexMatrix rho;
+    ComplexVector gamma;
+    G_MAT<int> multiplicities;
+    ComplexMatrix evOp;
+    ComplexMatrix TotalTimeIndependentHamiltonian;
+    ComplexMatrix magnetizationMatrix;
+    ComplexMatrix evMat;
+    ComplexMatrix eigenVectors;
+    RealVector eigenValues;
+    Real magnetization;
+    std::vector<Real> 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<ComplexMatrix> &output);
+    void _constructEffectiveSpinOperators(const ComplexVector& gamma, std::vector<ComplexMatrix> &output);
+    void _constructMagneticFieldHamiltonian(const std::vector<std::vector<ComplexMatrix> >& SpinOperators, const ComplexVector& gamma, const Real& Bx, const Real& By, const Real& Bz, ComplexMatrix& output);
+    void _constructJCouplingsHamiltonian(const ComplexMatrix &JCouplings, const std::vector<std::vector<ComplexMatrix> >& SpinOperators, ComplexMatrix &output);
+    void _calculate3DirsSpinOperators(std::vector<std::vector<ComplexMatrix> >& output);
+    void _fillSingleMultiplicityJMatrices(int multiplicity, std::vector<ComplexMatrix>& matrices);
+    void _fillJMatrices();
+    ComplexMatrix JCouplingHamiltonian;
+    ComplexMatrix MagneticFieldHamiltonian;
+    std::vector<std::vector<ComplexMatrix> > SpinOperators;
+    std::vector<ComplexMatrix> effectiveSpinOperators;
+    std::vector<std::vector<ComplexMatrix> > 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 <typename T>
+    std::vector<T> _vectorizeMatrix_upperTriangular_columnMajor(const G_MAT<T>& mat);
+    int directionToIndex(char dir);
+
+public:
+    SpinSystem();
+    SpinSystem(const ComplexVector& Gammas, const ComplexMatrix& JCouplings);
+    SpinSystem(const ComplexVector& Gammas, const ComplexMatrix& JCouplings, const G_MAT<int>& SpinMultiplicities);
+    void init(const ComplexVector& Gammas, const ComplexMatrix& JCouplings);
+    void init(const ComplexVector& Gammas, const ComplexMatrix& JCouplings, const G_MAT<int>& 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<std::complex<Real> >& 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<int>& spinMultiplicities);
+    void setGammas(const ComplexVector& gammas);
+    std::vector<Real> evolveDensityMatrix(const Real& sampleRate,
+                             const std::vector<std::string>& variableNames,
+                             const CX_MAT& variableData,
+                             int threads = 0,
+                             const std::set<char>& store_trace_directions = std::set<char>());
+    std::vector<Real> solveMagnetization(const long numOfPoints);
+    ~SpinSystem();
+    bool print = 1;
+    const Real& getMagnetization();
+    std::vector<Real> 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<Real> &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<int>& 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<Complex>()));
+    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<std::complex<Real> > >::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<std::complex<Real> > >::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<<evMat.size()<<std::endl;
+    return evMat;
+}
+
+TemplateHeader
+void
+SpinSystem TemplateTail::init(const ComplexVector& Gammas,
+                              const ComplexMatrix& JCouplings)
+{
+    G_MAT<int> SpinMultiplicities(Gammas.rows(),1,2);
+    this->init(Gammas,JCouplings,SpinMultiplicities);
+}
+
+TemplateHeader
+void
+SpinSystem TemplateTail::init(const ComplexVector& Gammas,
+                              const ComplexMatrix& JCouplings,
+                              const G_MAT<int>& 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<int>());
+    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)<<J_Couplings<<std::endl;
+
+}
+
+TemplateHeader
+void
+SpinSystem TemplateTail::setSpinMultiplicities(const Poly::Matrix<int> &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<Real>
+SpinSystem TemplateTail::evolveDensityMatrix(
+        const Real &sampleRate,
+        const std::vector<std::string> &variableNames,
+        const CX_MAT &variableData,
+        int threads,
+        const std::set<char>& 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<int>(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<Real>(1)/sampleRate;
+    long startPoint = internalTraceArray.size();
+    std::vector<Real> 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 <typename T>
+std::vector<T> SpinSystem TemplateTail::_vectorizeMatrix_upperTriangular_columnMajor(const G_MAT<T>& mat)
+{
+    std::vector<T> output((COLS(mat)*(COLS(mat)-1))/2);
+    typename G_MAT<T>::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<Real>
+SpinSystem TemplateTail::solveMagnetization(const long numOfPoints)
+{
+    std::vector<Real> TraceArray(numOfPoints);
+    if(print)
+    {
+        std::cout<<"Evolving time until t="<<numOfPoints/sampleRate<<":"<<std::endl;
+
+        std::cout.precision(std::ceil(std::log10(sampleRate)));
+        std::cout<<std::fixed;
+    }
+    for(long j = 0; j < numOfPoints; j++)
+    {
+        if(print)
+        {
+            std::cout<<"\r";
+            std::cout<<"t = "<< j/sampleRate << ":" <<  numOfPoints/sampleRate;
+            std::cout.flush();
+        }
+        magnetization = diagonalsSum + 2*std::real(std::accumulate(magnetizationMatrix.begin(), magnetizationMatrix.end(), Complex(0,0)));
+//        magnetization = 2*std::accumulate(magnetizationMatrix.cbegin(), magnetizationMatrix.cend(), Real{}, [](Real& acc, std::complex<Real> const c) -> Real& { return acc += c.real(); });
+
+        TraceArray[j] = magnetization;
+        std::transform( magnetizationMatrix.begin(), magnetizationMatrix.end(),
+                        evMat.begin(), magnetizationMatrix.begin(),
+                        std::multiplies<std::complex<Real> >());
+//        magnetization = real(accu(magnetizationMatrix));
+//        traceArray[j] = magnetization;
+//        ELEMENTWISEPRODUCT(magnetizationMatrix,evMat,magnetizationMatrix);
+    }
+    if(print)
+    {
+        std::cout.precision(10);
+        std::cout<<std::fixed;
+        std::cout<<"\r";
+        std::cout.flush();
+        std::cout<<std::endl;
+    }
+    this->internalTraceArray.resize(internalTraceArray.size() + numOfPoints);
+    return TraceArray;
+}
+
+TemplateHeader
+std::vector<Real>
+SpinSystem TemplateTail::solveMagnetization_parallel(const long numOfPoints,
+                                                     int threads)
+{
+//    Real dt = static_cast<Real>(1)/sampleRate;
+    std::vector<Real> TraceArray(numOfPoints);
+    std::vector<std::vector<std::complex<Real> > > par_evMat(threads);
+    std::vector<std::vector<std::complex<Real> > > par_magnetizationMatrix(threads);
+    for(long i = 0; i < threads; i++)
+    {
+        par_evMat[i].resize              (evMat.size() / threads + (i < static_cast<long>(evMat.size() % threads) ? 1 : 0));
+        par_magnetizationMatrix[i].resize(evMat.size() / threads + (i < static_cast<long>(evMat.size() % threads) ? 1 : 0));
+    }
+    std::vector<long> 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<long>(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<std::vector< std::complex<Real> > > 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="<<numOfPoints/sampleRate<<":"<<std::endl;
+        std::cout.precision(static_cast<std::streamsize>(std::ceil(std::log10(sampleRate))));
+        std::cout<<std::fixed;
+    }
+
+#pragma omp parallel for schedule(dynamic,1) num_threads(threads)
+    for(long i = 0; i < threads; i++)
+    {
+        for(long j = 0; j < numOfPoints; j++)
+        {
+            if(i == 0) //print only in the first thread
+            if(print)
+            {
+                if(j % 200 == 0)
+                {
+                    std::cout<<"\r";
+                    std::cout<<"t = "<< static_cast<Real>(j)/sampleRate << ":" <<  static_cast<Real>(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<Real>(0,0));
+//            par_traceArray[j][i] = std::accumulate(par_magnetizationMatrix[i].cbegin(), par_magnetizationMatrix[i].cend(), Real{}, [](Real const acc, std::complex<Real> 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<std::complex<Real> >());
+        }
+    }
+#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<Real>(0,0)));
+//        traceArray[i] = 2*std::accumulate(par_traceArray[i].cbegin(), par_traceArray[i].cend(), Real{}, [](Real const acc, std::complex<Real> const& c) { return acc + c.real(); });
+    }
+    if(print)
+    {
+        std::cout<<"\r";
+        std::cout<<"t = "<< numOfPoints/sampleRate << ":" <<  static_cast<Real>(numOfPoints)/sampleRate;
+        std::cout.flush();
+
+        std::cout.precision(10);
+        std::cout<<std::fixed;
+        std::cout<<"\r";
+        std::cout.flush();
+        std::cout<<std::endl;
+    }
+    this->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<Real>(1)/traceVal;
+    densityMatrixLock = 0;
+//    rho.ZEROS;
+//    rho[0] = 1;
+//    std::cout<<rho<<std::endl;
+    PRINT("Done populating thermally.");
+}
+
+TemplateHeader
+const
+std::vector<Real> &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<std::complex<Real> > &SpinSystem TemplateTail::getEvolutionOperator(Real samplingRate)
+//{
+//    evOp = arma::conv_to<std::vector<std::complex<Real> > >::from(vectorise(expmat(std::complex<Real>(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<ComplexMatrix>& matrices)
+{
+    matrices.resize(3);
+    const Real SpinQ = (static_cast<Real>(multiplicity)-static_cast<Real>(1))/static_cast<Real>(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<int> 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<ComplexMatrix> &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<ComplexMatrix> &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<std::vector<ComplexMatrix> >& 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<std::vector<ComplexMatrix> >& 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<std::vector<ComplexMatrix> > &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

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 <complex>
+
+#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<std::complex<_Real>>
+#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<std::complex<_Real>>
+#define MAT    Poly::Matrix<_Real>
+#define CX_MAT Poly::Matrix<std::complex<_Real>>
+#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
+