Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found

Target

Select target project
  • engelerp/rbcomb-simulation
1 result
Show changes
Commits on Source (8)
Expand all Hide whitespace changes
Inline Side-by-side
Showing
with 1156 additions and 855 deletions
parameter_calculation/Untitled.ipynb
View file @ 5f273894
......@@ -328,7 +328,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.7.4"
"version": "3.9.4"
}
},
"nbformat": 4,
......
projects/braidingTightBinding/bin/Makefile
View file @ 5f273894
......@@ -2,7 +2,7 @@ INCLUDEDIR = /Users/pascal/rbcomb_simulation/projects/braidingTightBinding/inclu
LIBDIR = /Users/pascal/rbcomb_simulation/projects/braidingTightBinding/lib
HDF5DIR = /Users/pascal/rbcomb_simulation/projects/braidingTightBinding/lib/src
CXX = /usr/local/bin/g++-9
CXX = /usr/local/bin/g++-10
CXXFLAGS = -std=c++2a -Wall -Wpedantic -I$(INCLUDEDIR) -I$(LIBDIR) -I${HDF5DIR} -O5 -mtune=native -march=native -framework Accelerate -lblas
run: all
......@@ -11,10 +11,22 @@ run: all
./main_braiding_zerostate ../results/data/zerostate/
python3 plot_zerostate.py
all: main_braiding_zerostate
run_nocoupling: main_nocoupling
mkdir -p ../results/data/nocoupling/
mkdir -p ../results/plots/nocoupling/
./main_nocoupling ../results/data/nocoupling/
python3 plot.py ../results/data/nocoupling/ ../results/plots/nocoupling/ "No Coupling" 0 10
main_braiding_zerostate: main_braiding_zerostate.cpp ${LIBDIR}/diagonalizer.hpp ${LIBDIR}/drum.hpp ${INCLUDEDIR}/grabber.hpp ${LIBDIR}/rk4_buffer.hpp ${LIBDIR}/rk4_stepper.hpp ${LIBDIR}/system.hpp ${LIBDIR}/system_parameters.hpp ${LIBDIR}/vec2.hpp ${INCLUDEDIR}/coupler_braid.hpp ${INCLUDEDIR}/driver_braid.hpp ${INCLUDEDIR}/drum_parameters_braid.hpp ${INCLUDEDIR}/drum_variables_braid.hpp ${INCLUDEDIR}/force_braid.hpp ${INCLUDEDIR}/matrix_element_calculator_braid.hpp ${INCLUDEDIR}/rbcomb_generator_braid.hpp ${INCLUDEDIR}/vortex.hpp
all: main_braiding_zerostate main_nocoupling
main_braiding_zerostate: main_braiding_zerostate.cpp ${LIBDIR}/diagonalizer.hpp ${LIBDIR}/drum.hpp ${INCLUDEDIR}/grabber.hpp ${LIBDIR}/rk4_buffer.hpp ${LIBDIR}/rk4_stepper.hpp ${LIBDIR}/system.hpp ${LIBDIR}/system_parameters.hpp ${LIBDIR}/vec2.hpp ${INCLUDEDIR}/coupler_braid.hpp ${INCLUDEDIR}/coupler_const.hpp ${INCLUDEDIR}/driver_braid.hpp ${INCLUDEDIR}/drum_parameters_braid.hpp ${INCLUDEDIR}/drum_variables_braid.hpp ${INCLUDEDIR}/force_braid.hpp ${INCLUDEDIR}/matrix_element_calculator_braid.hpp ${INCLUDEDIR}/rbcomb_generator_braid.hpp ${INCLUDEDIR}/vortex.hpp
$(CXX) $< $(CXXFLAGS) -o $@
main_nocoupling: main_nocoupling.cpp ${LIBDIR}/diagonalizer.hpp ${LIBDIR}/drum.hpp ${INCLUDEDIR}/grabber.hpp ${LIBDIR}/rk4_buffer.hpp ${LIBDIR}/rk4_stepper.hpp ${LIBDIR}/system.hpp ${LIBDIR}/system_parameters.hpp ${LIBDIR}/vec2.hpp ${INCLUDEDIR}/coupler_const.hpp ${INCLUDEDIR}/driver_braid.hpp ${INCLUDEDIR}/drum_parameters_braid.hpp ${INCLUDEDIR}/drum_variables_braid.hpp ${INCLUDEDIR}/force_braid.hpp ${INCLUDEDIR}/matrix_element_calculator_braid.hpp ${INCLUDEDIR}/rbcomb_generator_braid.hpp
$(CXX) $< $(CXXFLAGS) -o $@
plot:
python3 plot_zerostate.py
clean:
rm -rf main_braiding_zerostate ../results/data/* ../results/plots/*
projects/braidingTightBinding/bin/main_braiding_zerostate deleted 100755 → 0
View file @ 92fc03a4
File deleted
projects/braidingTightBinding/bin/main_braiding_zerostate.cpp
View file @ 5f273894
......@@ -25,111 +25,121 @@
template <typename iterator_t>
void writeToFile(const std::string filename, const std::string header, int line_length, iterator_t begin, iterator_t end) noexcept
{
//no line breaks if line_length == -1
if(line_length == -1)
line_length = std::distance(begin, end);
std::fstream file (filename, file.out);
file << header;
file << std::endl;
for(auto it = begin; it != end; ++it){
if(std::distance(begin, it) % line_length == 0)
file << std::endl;
file << *it << " ";
}
file << std::flush;
file.close();
//no line breaks if line_length == -1
if(line_length == -1){
line_length = std::distance(begin, end);
}
std::fstream file (filename, file.out);
file << header;
file << std::endl;
for(auto it = begin; it != end; ++it){
if(std::distance(begin, it) % line_length == 0){
file << std::endl;
}
file << *it << " ";
}
file << std::flush;
file.close();
}
int main(int argc, char** argv){
assert(argc == 2);
//shorthands
using value_t = double;
using params_t = DrumParametersBraid<value_t>;
using vars_t = DrumVariablesBraid<value_t>;
using buffer_t = RK4Buffer<value_t>;
using drum_t = Drum<value_t, params_t, vars_t, buffer_t>;
using coupler_t = CouplerBraid<value_t, drum_t, params_t>;
using driver_t = DriverBraid<value_t, drum_t>;
using sysparams_t = SystemParameters<coupler_t, driver_t>;
using force_t = ForceBraid<value_t, params_t, vars_t, buffer_t>;
using stepper_t = Rk4Stepper<value_t, params_t, vars_t, buffer_t, force_t>;
using generator_t = RbcombGeneratorBraid<value_t, params_t, vars_t, buffer_t>;
using grabber_t = Grabber<value_t, drum_t>;
using matelecalc_t = MatrixElementCalculatorBraid<value_t, params_t, vars_t, drum_t>;
using system_t = System<value_t, drum_t, grabber_t, sysparams_t, force_t, coupler_t, driver_t, stepper_t, matelecalc_t>;
using vortex_t = Vortex<value_t>;
using vec2_t = Vec2<value_t>;
value_t k0 =24674011002.723392;
value_t k1 = 493480220.0544679;
value_t k2 = 2467401100.2723393;
value_t c = 5.235987755982988;
vortex_t vortex0 (vec2_t(10.5,24.), 5., 1., 1., 1.); //as Eliska did it
vortex_t vortex1 (vec2_t(31.5,24.), 5., 1., 1., 1.); //as Eliska did it
vortex_t vortex2 (vec2_t(52.5,24.), 5., 1., 1., 1.); //as Eliska did it
//std::vector<vortex_t> vortices ({vortex0, vortex1, vortex2});
std::vector<vortex_t> vortices ({vortex0, vortex1, vortex2});
params_t drum_parameters (k0, k1, k2, c, vec2_t (0., 0.), 'A');
generator_t lattice_generator (37,15);
auto lat_pair = lattice_generator(drum_parameters);
//change edge potentials to push edge modes around
assert(argc == 2);
//shorthands
using value_t = double;
using params_t = DrumParametersBraid<value_t>;
using vars_t = DrumVariablesBraid<value_t>;
using buffer_t = RK4Buffer<value_t>;
using drum_t = Drum<value_t, params_t, vars_t, buffer_t>;
using coupler_t = CouplerBraid<value_t, drum_t, params_t>;
using driver_t = DriverBraid<value_t, drum_t>;
using sysparams_t = SystemParameters<coupler_t, driver_t>;
using force_t = ForceBraid<value_t, params_t, vars_t, buffer_t>;
using stepper_t = Rk4Stepper<value_t, params_t, vars_t, buffer_t, force_t>;
using generator_t = RbcombGeneratorBraid<value_t, params_t, vars_t, buffer_t>;
using grabber_t = Grabber<value_t, drum_t>;
using matelecalc_t = MatrixElementCalculatorBraid<value_t, params_t, vars_t, drum_t>;
using system_t = System<value_t, drum_t, grabber_t, sysparams_t, force_t, coupler_t, driver_t, stepper_t, matelecalc_t>;
using vortex_t = Vortex<value_t>;
using vec2_t = Vec2<value_t>;
value_t k0 = 24674011002.723392;
value_t k1 = 493480220.0544679;
value_t k2 = 2467401100.2723393;
value_t c = 5.235987755982988;
vortex_t vortex0 (vec2_t(10.5,24.), 5., 1., 1., 1.); //as Eliska did it
vortex_t vortex1 (vec2_t(31.5,24.), 5., 1., 1., 1.); //as Eliska did it
vortex_t vortex2 (vec2_t(52.5,24.), 5., 1., 1., 1.); //as Eliska did it
//std::vector<vortex_t> vortices ({vortex0, vortex1, vortex2});
std::vector<vortex_t> vortices ({vortex1});
params_t drum_parameters (k0, k1, k2, c, vec2_t (0., 0.), 'A');
generator_t lattice_generator (37,15);
auto lat_pair = lattice_generator(drum_parameters);
//change edge potentials to push edge modes around
value_t onsite_A = 27388152213.022964;
value_t onsite_B = 21959869792.42382;
for(size_t i = 0; i < lat_pair.first.size()-1; ++i){
for(size_t i = 0; i < lat_pair.first.size()-1; ++i){
if(lat_pair.second[i][0] == lat_pair.first.size() || lat_pair.second[i][1] == lat_pair.first.size() || lat_pair.second[i][2] == lat_pair.first.size()){
if(lat_pair.first[i].get_parameters().sublattice == 'A'){ //'A' site
lat_pair.first[i].get_parameters().k0 = onsite_A;
for(int j: lat_pair.second[i])
if(j != lat_pair.first.size())
for(int j: lat_pair.second[i]){
if(j != lat_pair.first.size()){
lat_pair.first[j].get_parameters().k0 = onsite_B;
}
}
}
else{ //'B' site
lat_pair.first[i].get_parameters().k0 = onsite_B;
for(int j: lat_pair.second[i])
if(j != lat_pair.first.size())
for(int j: lat_pair.second[i]){
if(j != lat_pair.first.size()){
lat_pair.first[j].get_parameters().k0 = onsite_A;
}
}
}
}
}
//set initial conditions
double x_0 = 0.;
double v_0 = 0.;
for(auto drum: lat_pair.first){
drum.get_variables().x = x_0;
drum.get_variables().xdot = v_0;
}
sysparams_t system_parameters (coupler_t(vortices), driver_t(), lat_pair.second);
force_t force;
stepper_t stepper;
//generate grabber
size_t save_every = 10;
std::string parameters_file = std::string(argv[1]) + "PARAMS.TXT";
std::string adjacency_file = std::string(argv[1]) + "ADJACENCY.TXT";
std::string dynamics_file = std::string(argv[1]) + "DYNAMICS.TXT";
grabber_t grabber (save_every, parameters_file, adjacency_file, dynamics_file);
sysparams_t system_parameters (coupler_t(vortices), driver_t(), lat_pair.second);
force_t force;
stepper_t stepper;
system_t system (1., 0.00001, lat_pair.first, stepper, force, system_parameters, grabber);
system.save();
//generate grabber
size_t save_every = 10;
std::string parameters_file = std::string(argv[1]) + "PARAMS.TXT";
std::string adjacency_file = std::string(argv[1]) + "ADJACENCY.TXT";
std::string dynamics_file = std::string(argv[1]) + "DYNAMICS.TXT";
grabber_t grabber (save_every, parameters_file, adjacency_file, dynamics_file);
//set up matrix element calculator
matelecalc_t mec;
system_t system (1., 0.00001, lat_pair.first, stepper, force, system_parameters, grabber);
system.save();
//set up diagonalizer
Diagonalizer diagonalizer;
//calculate eigenvalues
std::vector<value_t> matrix = system.get_matrix(mec);
//std::vector<value_t> evals = diagonalizer.ev(system.get_matrix(mec), lat_pair.first.size()-1);
std::pair<std::vector<value_t>, std::vector<value_t> > diag_pair = diagonalizer.evv(system.get_matrix(mec), lat_pair.first.size()-1);
//set up matrix element calculator
matelecalc_t mec;
//make eigenvalues to frequencies
std::transform(diag_pair.first.begin(), diag_pair.first.end(), diag_pair.first.begin(), [](value_t x) -> value_t {return std::sqrt(x)/(2.*M_PI);});
//set up diagonalizer
Diagonalizer diagonalizer;
//calculate eigenvalues
std::vector<value_t> matrix = system.get_matrix(mec);
//std::vector<value_t> evals = diagonalizer.ev(system.get_matrix(mec), lat_pair.first.size()-1);
std::pair<std::vector<value_t>, std::vector<value_t> > diag_pair = diagonalizer.evv(system.get_matrix(mec), lat_pair.first.size()-1);
//find the interesting eigenvector range
size_t ev_begin = 913 * diag_pair.first.size();
size_t ev_end = 943 * diag_pair.first.size();
//make eigenvalues to frequencies
std::transform(diag_pair.first.begin(), diag_pair.first.end(), diag_pair.first.begin(), [](value_t x) -> value_t {return std::sqrt(x)/(2.*M_PI);});
//write to file
writeToFile(std::string(argv[1]) + "eigenvalues.txt", "", 1, diag_pair.first.begin(), diag_pair.first.end());
writeToFile(std::string(argv[1]) + "eigenvectors.txt", "", diag_pair.first.size(), diag_pair.second.begin(), diag_pair.second.end());
//write to file
writeToFile(std::string(argv[1]) + "eigenvalues.txt", "", 1, diag_pair.first.begin(), diag_pair.first.end());
writeToFile(std::string(argv[1]) + "eigenvectors.txt", "", diag_pair.first.size(), diag_pair.second.begin(), diag_pair.second.end());
}
projects/braidingTightBinding/bin/main_nocoupling.cpp 0 → 100644
View file @ 5f273894
#include <coupler.hpp>
#include <diagonalizer.hpp>
#include <driver.hpp>
#include <drum.hpp>
#include <force.hpp>
#include <grabber.hpp>
#include <lattice_generator.hpp>
#include <rk4_buffer.hpp>
#include <rk4_stepper.hpp>
#include <system.hpp>
#include <system_parameters.hpp>
#include <vec2.hpp>
#include <coupler_const.hpp>
#include <driver_braid.hpp>
#include <drum_parameters_braid.hpp>
#include <drum_variables_braid.hpp>
#include <force_braid.hpp>
#include <matrix_element_calculator_braid.hpp>
#include <rbcomb_generator_braid.hpp>
#include <vortex.hpp>
#include <string>
#include <fstream>
#include <cassert>
template <typename iterator_t>
void writeToFile(const std::string filename, const std::string header, int line_length, iterator_t begin, iterator_t end) noexcept
{
//no line breaks if line_length == -1
if(line_length == -1){
line_length = std::distance(begin, end);
}
std::fstream file (filename, file.out);
file << header;
file << std::endl;
for(auto it = begin; it != end; ++it){
if(std::distance(begin, it) % line_length == 0){
file << std::endl;
}
file << *it << " ";
}
file << std::flush;
file.close();
}
int main(int argc, char** argv){
assert(argc == 2);
//shorthands
using value_t = double;
using params_t = DrumParametersBraid<value_t>;
using vars_t = DrumVariablesBraid<value_t>;
using buffer_t = RK4Buffer<value_t>;
using drum_t = Drum<value_t, params_t, vars_t, buffer_t>;
using coupler_t = CouplerConst<value_t, drum_t, params_t>;
using driver_t = DriverBraid<value_t, drum_t>;
using sysparams_t = SystemParameters<coupler_t, driver_t>;
using force_t = ForceBraid<value_t, params_t, vars_t, buffer_t>;
using stepper_t = Rk4Stepper<value_t, params_t, vars_t, buffer_t, force_t>;
using generator_t = RbcombGeneratorBraid<value_t, params_t, vars_t, buffer_t>;
using grabber_t = Grabber<value_t, drum_t>;
using matelecalc_t = MatrixElementCalculatorBraid<value_t, params_t, vars_t, drum_t>;
using system_t = System<value_t, drum_t, grabber_t, sysparams_t, force_t, coupler_t, driver_t, stepper_t, matelecalc_t>;
using vortex_t = Vortex<value_t>;
using vec2_t = Vec2<value_t>;
value_t k0 = 24674011002.723392;
value_t k1 = 0.;
value_t k2 = 2467401100.2723393;
value_t c = 5.235987755982988;
params_t drum_parameters (k0, k1, k2, c, vec2_t (0., 0.), 'A');
generator_t lattice_generator (37,15);
auto lat_pair = lattice_generator(drum_parameters);
//set initial conditions
double x_0 = 1.;
double xdot_0 = 0.;
for(drum_t& drum: lat_pair.first){
drum.get_variables().x = x_0;
drum.get_variables().xdot = xdot_0;
}
sysparams_t system_parameters (coupler_t(0.), driver_t(), lat_pair.second);
force_t force;
stepper_t stepper;
//generate grabber
double period = 0.00004;
double points_per_period = 100.;
size_t save_every = 10;
std::string parameters_file = std::string(argv[1]) + "PARAMS.TXT";
std::string adjacency_file = std::string(argv[1]) + "ADJACENCY.TXT";
std::string dynamics_file = std::string(argv[1]) + "DYNAMICS.TXT";
grabber_t grabber (save_every, parameters_file, adjacency_file, dynamics_file);
system_t system (10000*period, period/points_per_period, lat_pair.first, stepper, force, system_parameters, grabber);
system.simulate();
system.save();
//set up matrix element calculator
matelecalc_t mec;
//set up diagonalizer
Diagonalizer diagonalizer;
//calculate eigenvalues
std::vector<value_t> matrix = system.get_matrix(mec);
//std::vector<value_t> evals = diagonalizer.ev(system.get_matrix(mec), lat_pair.first.size()-1);
std::pair<std::vector<value_t>, std::vector<value_t> > diag_pair = diagonalizer.evv(system.get_matrix(mec), lat_pair.first.size()-1);
//make eigenvalues to frequencies
std::transform(diag_pair.first.begin(), diag_pair.first.end(), diag_pair.first.begin(), [](value_t x) -> value_t {return std::sqrt(x)/(2.*M_PI);});
//write to file
writeToFile(std::string(argv[1]) + "eigenvalues.txt", "", 1, diag_pair.first.begin(), diag_pair.first.end());
writeToFile(std::string(argv[1]) + "eigenvectors.txt", "", diag_pair.first.size(), diag_pair.second.begin(), diag_pair.second.end());
}
projects/braidingTightBinding/bin/plot.py 0 → 100644
View file @ 5f273894
import matplotlib.pyplot as plt
import numpy as np
import sys
if len(sys.argv) < 4:
print(f"Usage: python3.9 {argv[0]} ../path/to/data/ ../path/to/plots/ name [modes_range_start modes_range_end]")
sys.exit()
readpath = sys.argv[1]
writepath = sys.argv[2]
name = sys.argv[3]
#eigenvalues
ev = np.loadtxt(readpath+"eigenvalues.txt")
starting_mode = 1110
end_mode = 1180
if len(sys.argv) == 6:
starting_mode = int(sys.argv[4])
end_mode = int(sys.argv[5])
if end_mode - starting_mode > 70:
print(f"Maximum number of modes: 70 (entered: {end_mode - starting_mode})")
end_mode = starting_mode + 70
print(f"Using Modes {starting_mode}-{end_mode}")
plt.plot(np.array(list(range(starting_mode,end_mode))),ev[starting_mode:end_mode],'o', ms=1)
plt.title(f"Frequency levels under {name}")
plt.xlabel("Level index")
plt.ylabel("Frequency [Hz]")
plt.savefig(writepath + "Frequencies_zoomed.pdf", dpi = 500)
plt.clf()
plt.plot(ev[:],'o', ms=1)
plt.title(f"Frequency levels under {name}")
plt.xlabel("Level index")
plt.ylabel("Frequency [Hz]")
plt.savefig(writepath + "Frequencies.pdf", dpi = 500)
#eigenvectors
data = np.loadtxt(readpath+"eigenvectors.txt")
pos_data = np.loadtxt(readpath+"PARAMS.txt", dtype=str) #positions are rows 4 and 5
desired_length = len(data[0,:])
pos_x = pos_data[0:desired_length,4].astype('float')
pos_y = pos_data[:desired_length,5].astype('float')
plt.figure(figsize=(8,5*(end_mode - starting_mode + 1)))
for i in range(starting_mode,end_mode):
plt.subplot(int(end_mode-starting_mode),1,i-starting_mode+1)
cmax = np.amax(data[i,:])
cmin = np.amin(data[i,:])
if(-cmin > cmax):
cmax = -cmin
elif(-cmin <= cmax):
cmin = -cmax
plt.scatter(pos_x, pos_y, c=data[i,:], edgecolors='k', linewidth=0.1, cmap='seismic', s=16, vmin=cmin, vmax=cmax)
plt.colorbar()
plt.title(f"Mode {i+1}, ({ev[i]} Hz)")
plt.xlabel("x")
plt.ylabel("y")
plt.savefig(writepath + f"Modes.pdf", dpi = 40, bbox_inches='tight', pad_inches=0)
projects/braidingTightBinding/bin/plot_states.py 0 → 100644
View file @ 5f273894
import matplotlib.pyplot as plt
import numpy as np
import sys
readpath = "../results/data/zerostate/"
writepath = "../results/plots/zerostate/"
if(len(sys.argv) < 2):
print(f"Usage: python3.9 {sys.argv[0]} mode1 mode2 mode3 ...")
else:
num_modes = int(len(sys.argv)-1)
#eigenvalues
ev = np.loadtxt(readpath+"eigenvalues.txt")
#eigenvectors
data = np.loadtxt(readpath+"eigenvectors.txt")
pos_data = np.loadtxt(readpath+"PARAMS.txt", dtype=str) #positions are rows 4 and 5
desired_length = len(data[0,:])
pos_x = pos_data[0:desired_length,4].astype('float')
pos_y = pos_data[:desired_length,5].astype('float')
plt.figure(figsize=(8,5*num_modes))
count = 1
for str_i in sys.argv[1:]:
i = int(str_i)-1
plt.subplot(num_modes, 1, count)
cmax = np.amax(data[i,:])
cmin = np.amin(data[i,:])
if(-cmin > cmax):
cmax = -cmin
elif(-cmin <= cmax):
cmin = -cmax
plt.scatter(pos_x, pos_y, c=data[i,:], edgecolors='k', linewidth=0.1, cmap='seismic', s=16, vmin=cmin, vmax=cmax)
plt.colorbar()
plt.title(f"Mode {i+1}, ({ev[i]} Hz)")
plt.xlabel("x")
plt.ylabel("y")
count += 1
name = "Modes"
for str_i in sys.argv[1:]:
name += f"_{str_i}"
plt.savefig(writepath + name, dpi=500, bbox_inches='tight', pad_inches=0)
projects/braidingTightBinding/bin/plot_zerostate.py
View file @ 5f273894
......@@ -7,8 +7,8 @@ writepath = "../results/plots/zerostate/"
#eigenvalues
ev = np.loadtxt(readpath+"eigenvalues.txt")
starting_mode = 1115
end_mode = 1174
starting_mode = 1110
end_mode = 1180
plt.plot(np.array(list(range(starting_mode,end_mode))),ev[starting_mode:end_mode],'o', ms=1)
plt.title("Frequency levels under vortex Kekule")
......@@ -35,7 +35,7 @@ pos_y = pos_data[:desired_length,5].astype('float')
plt.figure(figsize=(8,300))
for i in range(starting_mode,end_mode):
plt.subplot(60,1,i-starting_mode+1)
plt.subplot(int(end_mode-starting_mode),1,i-starting_mode+1)
cmax = np.amax(data[i,:])
cmin = np.amin(data[i,:])
if(-cmin > cmax):
......@@ -48,4 +48,4 @@ for i in range(starting_mode,end_mode):
plt.xlabel("x")
plt.ylabel("y")
plt.savefig(writepath + f"Modes.pdf", dpi = 50, bbox_inches='tight', pad_inches=0)
plt.savefig(writepath + f"Modes.pdf", dpi = 40, bbox_inches='tight', pad_inches=0)
projects/braidingTightBinding/include/coupler_braid.hpp
View file @ 5f273894
......@@ -6,55 +6,57 @@
template <typename value_t, typename drum_t, typename params_t>
class CouplerBraid: public Coupler<value_t, drum_t>{
public:
//constructors
CouplerBraid(const std::vector<Vortex<value_t> >& vortices): vortices_(vortices) {}
~CouplerBraid() = default;
void precompute(const value_t t_end, const value_t dt, const std::vector<drum_t>& drum_vec) noexcept final override
{
//we need the drum parameters later on, so extract it
params_.reserve(drum_vec.size());
for(auto d: drum_vec)
params_.push_back(d.get_parameters());
//we also need the neighbour vectors nicely indexed
for(auto d: drum_vec){
if(d.get_parameters().sublattice == 'A'){
n_vecs_.push_back(d.get_parameters().s0);
n_vecs_.push_back(d.get_parameters().s1);
n_vecs_.push_back(d.get_parameters().s2);
break;
}
}
}
void step(value_t dt) noexcept final
{
//move vortices or alpha here
}
value_t operator()(const size_t drum_index, const size_t neighbour_index) const noexcept final override
{
Vec2<value_t> v = params_[drum_index].position;
Vec2<value_t> s = n_vecs_[neighbour_index];
//the formula is only valid on the 'A' sublattice
if(params_[drum_index].sublattice == 'B')
v = v - s;
std::complex<value_t> ret_val = 1.; //vortices are multiplicative
for(auto tex: vortices_){
ret_val *= tex.distortion(v);
}
//add kekule
ret_val *= vortices_.begin()->kekule(v, s);
return 1. + std::real<value_t>(ret_val); //return with added offset, as in eliska's code
}
private:
std::vector<Vortex<value_t> > vortices_;
std::vector<params_t> params_;
std::vector<Vec2<value_t> > n_vecs_;
public:
//constructors
CouplerBraid(const std::vector<Vortex<value_t> >& vortices): vortices_(vortices) {}
~CouplerBraid() = default;
void precompute(const value_t t_end, const value_t dt, const std::vector<drum_t>& drum_vec) noexcept final override
{
//we need the drum parameters later on, so extract it
params_.reserve(drum_vec.size());
for(auto d: drum_vec){
params_.push_back(d.get_parameters());
}
//we also need the neighbour vectors nicely indexed
for(auto d: drum_vec){
if(d.get_parameters().sublattice == 'A'){
n_vecs_.push_back(d.get_parameters().s0);
n_vecs_.push_back(d.get_parameters().s1);
n_vecs_.push_back(d.get_parameters().s2);
break;
}
}
}
void step(value_t dt) noexcept final
{
//move vortices or alpha here
}
value_t operator()(const size_t drum_index, const size_t neighbour_index) const noexcept final override
{
Vec2<value_t> v = params_[drum_index].position;
Vec2<value_t> s = n_vecs_[neighbour_index];
//the formula is only valid on the 'A' sublattice
if(params_[drum_index].sublattice == 'B'){
v = v - s;
}
std::complex<value_t> ret_val = 1.; //vortices are multiplicative
for(auto tex: vortices_){
ret_val *= tex.distortion(v);
}
//add kekule
ret_val *= vortices_.begin()->kekule(v, s);
return 1. + std::real<value_t>(ret_val); //return with added offset, as in eliska's code
}
private:
std::vector<Vortex<value_t> > vortices_;
std::vector<params_t> params_;
std::vector<Vec2<value_t> > n_vecs_;
};
#endif
projects/braidingTightBinding/include/coupler_const.hpp 0 → 100644
View file @ 5f273894
#ifndef COUPLER_CONST_HPP_INCLUDED
#define COUPLER_CONST_HPP_INCLUDED
#include <vector>
#include <coupler.hpp>
template <typename value_t, typename drum_t, typename params_t>
class CouplerConst: public Coupler<value_t, drum_t>{
public:
//constructors
CouplerConst(value_t value): value_(value) {}
~CouplerConst() = default;
void precompute(const value_t t_end, const value_t dt, const std::vector<drum_t>& drum_vec) noexcept final override
{
//nothing to do
}
void step(value_t dt) noexcept final
{
//nothing to do
}
value_t operator()(const size_t drum_index, const size_t neighbour_index) const noexcept final override
{
return value_;
}
private:
value_t value_;
};
#endif
projects/braidingTightBinding/include/driver_braid.hpp
View file @ 5f273894
......@@ -5,23 +5,23 @@
template <typename value_t, typename drum_t>
class DriverBraid: public Driver<value_t, drum_t>{
public:
//constructors
DriverBraid() = default;
~DriverBraid() = default;
public:
//constructors
DriverBraid() = default;
~DriverBraid() = default;
void precompute(const value_t t_end, const value_t dt, const std::vector<drum_t>& drum_vec) noexcept final override
{
void precompute(const value_t t_end, const value_t dt, const std::vector<drum_t>& drum_vec) noexcept final override
{
}
void step(value_t dt) noexcept final
{
}
void step(value_t dt) noexcept final
{
}
value_t operator()(const size_t drum_index) const noexcept final override
{
return 0;
}
}
value_t operator()(const size_t drum_index) const noexcept final override
{
return 0;
}
};
#endif
projects/braidingTightBinding/include/drum_parameters_braid.hpp
View file @ 5f273894
......@@ -7,58 +7,59 @@
//TODO: Create an interface. It would be a trivial one, though.
template <typename value_t>
class DrumParametersBraid{
public:
/* Arguments:
* k0: diagonal matrix element
* k1: coupling matrix element prefactor 1
* k2: coupling matrix element prefactor 2
* c: damping coefficient
* position: position of the drum in real space
* sublattice: 'A'/'a' or 'B'/'b', identifies sublattice drum is part of
*/
DrumParametersBraid(const value_t k0, const value_t k1, const value_t k2, const value_t c, const Vec2<value_t>& position, char sublattice) noexcept
: k0(k0), k1(k1), k2(k2), c(c), position(position), sublattice(sublattice)
{
generate_sublattice(sublattice, value_t(1.));
}
DrumParametersBraid() = delete; //no default initialization allowed
DrumParametersBraid(const DrumParametersBraid&) = default;
DrumParametersBraid& operator=(const DrumParametersBraid&) = default;
~DrumParametersBraid() = default;
public:
/* Arguments:
* k0: diagonal matrix element
* k1: coupling matrix element prefactor 1
* k2: coupling matrix element prefactor 2
* c: damping coefficient
* position: position of the drum in real space
* sublattice: 'A'/'a' or 'B'/'b', identifies sublattice drum is part of
*/
DrumParametersBraid(const value_t k0, const value_t k1, const value_t k2, const value_t c, const Vec2<value_t>& position, char sublattice) noexcept
: k0(k0), k1(k1), k2(k2), c(c), position(position), sublattice(sublattice)
{
generate_sublattice(sublattice, value_t(1.));
}
DrumParametersBraid() = delete; //no default initialization allowed
DrumParametersBraid(const DrumParametersBraid&) = default;
DrumParametersBraid& operator=(const DrumParametersBraid&) = default;
~DrumParametersBraid() = default;
private:
int generate_sublattice(char sublattice, value_t lc) noexcept{
if(sublattice == 'A' or sublattice == 'a'){
s0 = Vec2<value_t> (0, -lc);
s1 = Vec2<value_t> (-std::sqrt(3)*lc/2., lc/2.);
s2 = Vec2<value_t> (std::sqrt(3)*lc/2., lc/2.);
return 0;
}
else if(sublattice == 'B' or sublattice == 'b'){
s0 = -1*(Vec2<value_t> (0, -lc));
s1 = -1*(Vec2<value_t> (-std::sqrt(3)*lc/2., lc/2.));
s2 = -1*(Vec2<value_t> (std::sqrt(3)*lc/2., lc/2.));
return 0;
}
else
return -1; //invalid sublattice identifier
}
private:
int generate_sublattice(char sublattice, value_t lc) noexcept{
if(sublattice == 'A' or sublattice == 'a'){
s0 = Vec2<value_t> (0, -lc);
s1 = Vec2<value_t> (-std::sqrt(3)*lc/2., lc/2.);
s2 = Vec2<value_t> (std::sqrt(3)*lc/2., lc/2.);
return 0;
}
else if(sublattice == 'B' or sublattice == 'b'){
s0 = -1*(Vec2<value_t> (0, -lc));
s1 = -1*(Vec2<value_t> (-std::sqrt(3)*lc/2., lc/2.));
s2 = -1*(Vec2<value_t> (std::sqrt(3)*lc/2., lc/2.));
return 0;
}
else{
return -1; //invalid sublattice identifier
}
}
public:
value_t k0; //onsite prefactor
value_t k1; //coupling prefactor 1
value_t k2; //coupling prefactor 2
value_t c; //damping coefficient
char sublattice; //sublattice the drum is part of (governs electrode geometry)
Vec2<value_t> s0, s1, s2; //vectors connecting to neighbours, ordered according to drum.hpp
Vec2<value_t> position; //position of the drum in real space
public:
value_t k0; //onsite prefactor
value_t k1; //coupling prefactor 1
value_t k2; //coupling prefactor 2
value_t c; //damping coefficient
char sublattice; //sublattice the drum is part of (governs electrode geometry)
Vec2<value_t> s0, s1, s2; //vectors connecting to neighbours, ordered according to drum.hpp
Vec2<value_t> position; //position of the drum in real space
};
template<typename value_t>
std::ostream& operator<<(std::ostream& os, const DrumParametersBraid<value_t>& dp)
{
return os << "Sublattice: " << dp.sublattice << ", k0: " << dp.k0 << ", k1: " << dp.k1 << ", k2: " << dp.k2 << std::endl << "s0: " << dp.s0 << ", s1: " << dp.s1 << ", s2: " << dp.s2;
return os << "Sublattice: " << dp.sublattice << ", k0: " << dp.k0 << ", k1: " << dp.k1 << ", k2: " << dp.k2 << std::endl << "s0: " << dp.s0 << ", s1: " << dp.s1 << ", s2: " << dp.s2;
}
......
projects/braidingTightBinding/include/drum_variables_braid.hpp
View file @ 5f273894
......@@ -8,16 +8,16 @@
//TODO: Generate an interface. It would be trivial, though.
template <typename value_t>
class DrumVariablesBraid{
public:
//constructors
DrumVariablesBraid() = default;
DrumVariablesBraid& operator=(const DrumVariablesBraid&) = default;
~DrumVariablesBraid() = default;
public:
//constructors
DrumVariablesBraid() = default;
DrumVariablesBraid& operator=(const DrumVariablesBraid&) = default;
~DrumVariablesBraid() = default;
value_t t0, t1, t2; //t + dt, i.e. baseline plus vortices along bond 0,1,2
value_t V; //Coupling to central electrode
value_t x, xdot; //Current Elongation (position) and velocity
value_t x_temp, xdot_temp; //used as calculation arguments in rk steps
value_t t0, t1, t2; //t + dt, i.e. baseline plus vortices along bond 0,1,2
value_t V; //Coupling to central electrode
value_t x, xdot; //Current Elongation (position) and velocity
value_t x_temp, xdot_temp; //used as calculation arguments in rk steps
};
#endif
projects/braidingTightBinding/include/force_braid.hpp
View file @ 5f273894
......@@ -9,35 +9,35 @@
template <typename value_t, typename params_t, typename vars_t, typename buffer_t>
class ForceBraid: public Force<value_t, params_t, vars_t, buffer_t>{
public:
ForceBraid() = default;
~ForceBraid() = default;
public:
ForceBraid() = default;
~ForceBraid() = default;
value_t operator()(
const Drum<value_t, params_t, vars_t, buffer_t>& drum, //Drum we're calculating the force on
const Drum<value_t, params_t, vars_t, buffer_t>& neighbour0, //Neighbour 0
const Drum<value_t, params_t, vars_t, buffer_t>& neighbour1, //Neighbour 1
const Drum<value_t, params_t, vars_t, buffer_t>& neighbour2, //Neighbour 2
const value_t time) const noexcept final override
{
//fetch data
//note that no copying is done here, these are all references
//TODO: implement drive
const params_t& drum_params = drum.get_parameters(); //parameters of drum
const vars_t& drum_vars = drum.get_variables(); //variables of drum
const vars_t& n0_vars = neighbour0.get_variables(); //variables of neighbour 0
const vars_t& n1_vars = neighbour1.get_variables(); //variables of neighbour 1
const vars_t& n2_vars = neighbour2.get_variables(); //variables of neighbour 2
value_t operator()(
const Drum<value_t, params_t, vars_t, buffer_t>& drum, //Drum we're calculating the force on
const Drum<value_t, params_t, vars_t, buffer_t>& neighbour0, //Neighbour 0
const Drum<value_t, params_t, vars_t, buffer_t>& neighbour1, //Neighbour 1
const Drum<value_t, params_t, vars_t, buffer_t>& neighbour2, //Neighbour 2
const value_t time) const noexcept final override
{
//fetch data
//note that no copying is done here, these are all references
//TODO: implement drive
const params_t& drum_params = drum.get_parameters(); //parameters of drum
const vars_t& drum_vars = drum.get_variables(); //variables of drum
const vars_t& n0_vars = neighbour0.get_variables(); //variables of neighbour 0
const vars_t& n1_vars = neighbour1.get_variables(); //variables of neighbour 1
const vars_t& n2_vars = neighbour2.get_variables(); //variables of neighbour 2
value_t part1 = - drum_params.k0 * drum_vars.x_temp;
value_t part2 = - drum_params.c * drum_vars.xdot_temp;
value_t part3 = 0.;
part3 -= drum_params.k1 + drum_params.k2 * n0_vars.t0 * n0_vars.x_temp;
part3 -= drum_params.k1 + drum_params.k2 * n1_vars.t1 * n1_vars.x_temp;
part3 -= drum_params.k1 + drum_params.k2 * n2_vars.t2 * n2_vars.x_temp;
value_t part1 = - drum_params.k0 * drum_vars.x_temp;
value_t part2 = - drum_params.c * drum_vars.xdot_temp;
value_t part3 = 0.;
part3 -= drum_params.k1 + drum_params.k2 * n0_vars.t0 * n0_vars.x_temp;
part3 -= drum_params.k1 + drum_params.k2 * n1_vars.t1 * n1_vars.x_temp;
part3 -= drum_params.k1 + drum_params.k2 * n2_vars.t2 * n2_vars.x_temp;
return part1 + part2 + part3;
}
return part1 + part2 + part3;
}
};
......
projects/braidingTightBinding/include/grabber.hpp
View file @ 5f273894
......@@ -7,134 +7,137 @@
template<typename value_t>
struct DataStorage{
DataStorage(const size_t num_drums)
{
x.reserve(num_drums);
xdot.reserve(num_drums);
t0.reserve(num_drums);
t1.reserve(num_drums);
t2.reserve(num_drums);
drive.reserve(num_drums);
}
~DataStorage() = default;
DataStorage(const size_t num_drums)
{
x.reserve(num_drums);
xdot.reserve(num_drums);
t0.reserve(num_drums);
t1.reserve(num_drums);
t2.reserve(num_drums);
drive.reserve(num_drums);
}
~DataStorage() = default;
value_t time;
std::vector<value_t> x;
std::vector<value_t> xdot;
std::vector<value_t> t0;
std::vector<value_t> t1;
std::vector<value_t> t2;
std::vector<value_t> drive;
value_t time;
std::vector<value_t> x;
std::vector<value_t> xdot;
std::vector<value_t> t0;
std::vector<value_t> t1;
std::vector<value_t> t2;
std::vector<value_t> drive;
};
template <typename value_t, typename drum_t>
class Grabber{
public:
Grabber(const size_t grab_every, const std::string params_file, const std::string adjacency_file, const std::string dynamic_file)
: grab_every_(grab_every), grab_cnt_(0), par_f_(params_file), adj_f_(adjacency_file), dyn_f_(dynamic_file) {}
Grabber() = delete; //no default initialization
Grabber(const Grabber&) = default;
~Grabber() = default;
public:
Grabber(const size_t grab_every, const std::string params_file, const std::string adjacency_file, const std::string dynamic_file)
: grab_every_(grab_every), grab_cnt_(0), par_f_(params_file), adj_f_(adjacency_file), dyn_f_(dynamic_file) {}
Grabber() = delete; //no default initialization
Grabber(const Grabber&) = default;
~Grabber() = default;
//call in system constructor
void init(const value_t t_end, const value_t dt, const std::vector<drum_t>& drums, const std::vector<std::vector<int> >& adjacency) noexcept
{
drums_ = drums;
adjacency_ = adjacency;
t_end_ = t_end;
dt_ = dt;
}
//call in system constructor
void init(const value_t t_end, const value_t dt, const std::vector<drum_t>& drums, const std::vector<std::vector<int> >& adjacency) noexcept
{
drums_ = drums;
adjacency_ = adjacency;
t_end_ = t_end;
dt_ = dt;
}
//call after each simulation step, checks itself if it should save
bool grab(const std::vector<drum_t>& drums, const value_t time)
{
if(grab_cnt_++ % grab_every_ == 0){
data_.push_back(DataStorage<value_t>(drums_.size()));
data_.back().time = time;
for(auto d: drums){
data_.back().x.push_back(d.get_variables().x);
data_.back().xdot.push_back(d.get_variables().xdot);
data_.back().t0.push_back(d.get_variables().t0);
data_.back().t1.push_back(d.get_variables().t1);
data_.back().t2.push_back(d.get_variables().t2);
data_.back().drive.push_back(d.get_variables().V);
}
return true;
}
else
return false;
}
//call after each simulation step, checks itself if it should save
bool grab(const std::vector<drum_t>& drums, const value_t time)
{
if(grab_cnt_++ % grab_every_ == 0){
data_.push_back(DataStorage<value_t>(drums_.size()));
data_.back().time = time;
for(auto d: drums){
data_.back().x.push_back(d.get_variables().x);
data_.back().xdot.push_back(d.get_variables().xdot);
data_.back().t0.push_back(d.get_variables().t0);
data_.back().t1.push_back(d.get_variables().t1);
data_.back().t2.push_back(d.get_variables().t2);
data_.back().drive.push_back(d.get_variables().V);
}
return true;
}
else{
return false;
}
}
//prints the data out. call at end of it all.
bool save(){
//TODO: Check for file open
//TODO: Overload operator<< for DrumParameters and DataStorage.
//print parameters
std::fstream par (par_f_, par.out);
par << "# The last drum is to be ignored\n";
par << "# k0 \t k1 \t k2 \t c \t pos_x \t pos_y \t sublattice\n";
for(auto d: drums_){
par << d.get_parameters().k0 << " \t";
par << d.get_parameters().k1 << " \t";
par << d.get_parameters().k2 << " \t";
par << d.get_parameters().c << " \t";
par << d.get_parameters().position.x() << " \t";
par << d.get_parameters().position.y() << " \t";
par << d.get_parameters().sublattice << "\n";
}
par.close();
//prints the data out. call at end of it all.
bool save(){
//TODO: Check for file open
//TODO: Overload operator<< for DrumParameters and DataStorage.
//print parameters
std::fstream par (par_f_, par.out);
par << "# The last drum is to be ignored\n";
par << "# k0 \t k1 \t k2 \t c \t pos_x \t pos_y \t sublattice\n";
for(auto d: drums_){
par << d.get_parameters().k0 << " \t";
par << d.get_parameters().k1 << " \t";
par << d.get_parameters().k2 << " \t";
par << d.get_parameters().c << " \t";
par << d.get_parameters().position.x() << " \t";
par << d.get_parameters().position.y() << " \t";
par << d.get_parameters().sublattice << "\n";
}
par.close();
//print adjacency vectors
std::fstream adj (adj_f_, adj.out);
adj << "# <-1> means <no neighbour here>\n";
adj << "# The last drum is to be ignored\n";
adj << "# n0 \t n1 \t n2\n";
for(auto v: adjacency_){
for(auto i: v){
if(i == drums_.size()-1)
adj << -1 << " \t"; //no neighbour present
else
adj << i << " \t"; //neighbour
}
adj << "\n";
}
adj.close();
//print adjacency vectors
std::fstream adj (adj_f_, adj.out);
adj << "# <-1> means <no neighbour here>\n";
adj << "# The last drum is to be ignored\n";
adj << "# n0 \t n1 \t n2\n";
for(auto v: adjacency_){
for(auto i: v){
if(i == drums_.size()-1){
adj << -1 << " \t"; //no neighbour present
}
else{
adj << i << " \t"; //neighbour
}
}
adj << "\n";
}
adj.close();
std::fstream dyn (dyn_f_, dyn.out);
dyn << "# If adjacency of a bond is -1 in " << adj_f_;
dyn << ", the coupling is to be ignored\n";
dyn << "# The last drum is to be ignored\n";
dyn << "# Number of drums : " << drums_.size() << "\n";
dyn << "# t \t drum_index \t x \t xdot \t t0 \t t1 \t t2 \t drive\n";
for(auto entry: data_){
for(size_t i = 0; i < entry.x.size(); ++i){
dyn << entry.time << " \t";
dyn << i << " \t";
dyn << entry.x[i] << " \t";
dyn << entry.xdot[i] << " \t";
dyn << entry.t0[i] << " \t";
dyn << entry.t1[i] << " \t";
dyn << entry.t2[i] << " \t";
dyn << entry.drive[i] << "\n";
}
}
dyn.close();
std::fstream dyn (dyn_f_, dyn.out);
dyn << "# If adjacency of a bond is -1 in " << adj_f_;
dyn << ", the coupling is to be ignored\n";
dyn << "# The last drum is to be ignored\n";
dyn << "# Number of drums : " << drums_.size() << "\n";
dyn << "# t \t drum_index \t x \t xdot \t t0 \t t1 \t t2 \t drive\n";
for(auto entry: data_){
for(size_t i = 0; i < entry.x.size(); ++i){
dyn << entry.time << " \t";
dyn << i << " \t";
dyn << entry.x[i] << " \t";
dyn << entry.xdot[i] << " \t";
dyn << entry.t0[i] << " \t";
dyn << entry.t1[i] << " \t";
dyn << entry.t2[i] << " \t";
dyn << entry.drive[i] << "\n";
}
}
dyn.close();
return true;
}
return true;
}
private:
size_t grab_every_; //grab data every grab_every_ steps
size_t grab_cnt_; //grab count
value_t t_end_, dt_; //maximal simulation time and timestep
std::vector<std::vector<int> > adjacency_; //adjacency vector
std::vector<drum_t> drums_; //drums in initial configuration, grab params from here
std::list<DataStorage<value_t> > data_; //here goes the dynamical data
//filenames
std::string par_f_; //parameters file
std::string adj_f_; //adjacency file
std::string dyn_f_; //dynamics file (this is the gold)
private:
size_t grab_every_; //grab data every grab_every_ steps
size_t grab_cnt_; //grab count
value_t t_end_, dt_; //maximal simulation time and timestep
std::vector<std::vector<int> > adjacency_; //adjacency vector
std::vector<drum_t> drums_; //drums in initial configuration, grab params from here
std::list<DataStorage<value_t> > data_; //here goes the dynamical data
//filenames
std::string par_f_; //parameters file
std::string adj_f_; //adjacency file
std::string dyn_f_; //dynamics file (this is the gold)
};
#endif
projects/braidingTightBinding/include/matrix_element_calculator_braid.hpp
View file @ 5f273894
......@@ -4,32 +4,32 @@
template <typename value_t, typename params_t, typename vars_t, typename drum_t>
class MatrixElementCalculatorBraid: public MatrixElementCalculator<value_t, params_t, vars_t, drum_t>{
public:
MatrixElementCalculatorBraid() = default;
~MatrixElementCalculatorBraid() = default;
public:
MatrixElementCalculatorBraid() = default;
~MatrixElementCalculatorBraid() = default;
//diagonal element at (index, index)
value_t operator()(const size_t index, const std::vector<drum_t>& drums) const noexcept final override
{
return drums[index].get_parameters().k0;
}
//diagonal element at (index, index)
value_t operator()(const size_t index, const std::vector<drum_t>& drums) const noexcept final override
{
return drums[index].get_parameters().k0;
}
//coupling elements at (index1, index2)
//for neighbour 0
value_t operator()(const size_t index1, const size_t index2, const std::vector<drum_t>& drums) const noexcept final override
{
return drums[index1].get_parameters().k1 + drums[index1].get_parameters().k2 * drums[index2].get_variables().t0;
}
//for neighbour 1
value_t operator()(const size_t index1, const size_t index2, const std::vector<drum_t>& drums, const int) const noexcept final override
{
return drums[index1].get_parameters().k1 + drums[index1].get_parameters().k2 * drums[index2].get_variables().t1;
}
//for neighbour 2
value_t operator()(const size_t index1, const size_t index2, const std::vector<drum_t>& drums, const int, const int) const noexcept final override
{
return drums[index1].get_parameters().k1 + drums[index1].get_parameters().k2 * drums[index2].get_variables().t2;
}
//coupling elements at (index1, index2)
//for neighbour 0
value_t operator()(const size_t index1, const size_t index2, const std::vector<drum_t>& drums) const noexcept final override
{
return drums[index1].get_parameters().k1 + drums[index1].get_parameters().k2 * drums[index2].get_variables().t0;
}
//for neighbour 1
value_t operator()(const size_t index1, const size_t index2, const std::vector<drum_t>& drums, const int) const noexcept final override
{
return drums[index1].get_parameters().k1 + drums[index1].get_parameters().k2 * drums[index2].get_variables().t1;
}
//for neighbour 2
value_t operator()(const size_t index1, const size_t index2, const std::vector<drum_t>& drums, const int, const int) const noexcept final override
{
return drums[index1].get_parameters().k1 + drums[index1].get_parameters().k2 * drums[index2].get_variables().t2;
}
};
#endif
projects/braidingTightBinding/include/vortex.hpp
View file @ 5f273894
......@@ -6,43 +6,43 @@
template <typename value_t>
class Vortex{
public:
//constructors
Vortex(const Vec2<value_t> position, const value_t l0, const value_t alpha, const value_t delta, const value_t a) noexcept: position_(position), l0_(l0), alpha_(alpha), delta_(delta), K_(value_t(4.*3.141592653589793/(3.*std::sqrt(3)*a)), value_t(0.)) {}
Vortex(const Vec2<value_t> position, const Vortex& o) noexcept: position_(position), l0_(o.l0_), alpha_(o.alpha_), delta_(o.delta_), K_(o.K_) {}
Vortex() = default;
Vortex(const Vortex&) = default;
Vortex& operator=(const Vortex&) = default;
~Vortex() = default;
public:
//constructors
Vortex(const Vec2<value_t> position, const value_t l0, const value_t alpha, const value_t delta, const value_t a) noexcept: position_(position), l0_(l0), alpha_(alpha), delta_(delta), K_(value_t(4.*3.141592653589793/(3.*std::sqrt(3)*a)), value_t(0.)) {}
Vortex(const Vec2<value_t> position, const Vortex& o) noexcept: position_(position), l0_(o.l0_), alpha_(o.alpha_), delta_(o.delta_), K_(o.K_) {}
Vortex() = default;
Vortex(const Vortex&) = default;
Vortex& operator=(const Vortex&) = default;
~Vortex() = default;
//access
value_t l0() const noexcept { return l0_; }
value_t alpha() const noexcept { return alpha_; }
value_t delta() const noexcept { return delta_; }
Vec2<value_t>& position() noexcept { return position_; }
const Vec2<value_t>& position() const noexcept { return position_; }
//access
value_t l0() const noexcept { return l0_; }
value_t alpha() const noexcept { return alpha_; }
value_t delta() const noexcept { return delta_; }
Vec2<value_t>& position() noexcept { return position_; }
const Vec2<value_t>& position() const noexcept { return position_; }
//modifiers
void set_l0(value_t l0) noexcept { l0_ = l0; }
void set_alpha(value_t alpha) noexcept { alpha_ = alpha; }
void set_delta(value_t delta) noexcept { delta_ = delta; }
void set_position(const Vec2<value_t>& position) noexcept { position_ = position; }
void move_by(const Vec2<value_t>& translation) noexcept { position_ += translation; }
void move_to(const Vec2<value_t>& position) noexcept { position_ = position; }
//modifiers
void set_l0(value_t l0) noexcept { l0_ = l0; }
void set_alpha(value_t alpha) noexcept { alpha_ = alpha; }
void set_delta(value_t delta) noexcept { delta_ = delta; }
void set_position(const Vec2<value_t>& position) noexcept { position_ = position; }
void move_by(const Vec2<value_t>& translation) noexcept { position_ += translation; }
void move_to(const Vec2<value_t>& position) noexcept { position_ = position; }
//functional
//there is no divide_by_zero check for efficiency reasons. l0_ can not be zero upon call.
std::complex<value_t> distortion(const Vec2<value_t>& v) const
{
return std::polar<value_t>(delta_*std::tanh(v.r_wrt(position_)/l0_), alpha_ - v.phi_wrt(position_));
}
std::complex<value_t> kekule(const Vec2<value_t>& v, const Vec2<value_t>& s) const{
return std::polar<value_t>(1., K_*(s + 2.*(v - position_)));
}
private:
Vec2<value_t> position_;
const Vec2<value_t> K_;
value_t l0_, alpha_, delta_;
//functional
//there is no divide_by_zero check for efficiency reasons. l0_ can not be zero upon call.
std::complex<value_t> distortion(const Vec2<value_t>& v) const
{
return std::polar<value_t>(delta_*std::tanh(v.r_wrt(position_)/l0_), alpha_ - v.phi_wrt(position_));
}
std::complex<value_t> kekule(const Vec2<value_t>& v, const Vec2<value_t>& s) const{
return std::polar<value_t>(1., K_*(s + 2.*(v - position_)));
}
private:
Vec2<value_t> position_;
const Vec2<value_t> K_;
value_t l0_, alpha_, delta_;
};
#endif
projects/braidingTightBinding/lib/coupler.hpp
View file @ 5f273894
......@@ -4,14 +4,14 @@
template <typename value_t, typename drum_t>
class Coupler{
public:
//constructors
Coupler() = default;
~Coupler() = default;
public:
//constructors
Coupler() = default;
~Coupler() = default;
virtual void precompute(const value_t t_end, const value_t dt, const std::vector<drum_t>& drum_vec) noexcept = 0;
virtual void step(value_t dt) noexcept = 0;
virtual value_t operator()(const size_t drum_index, const size_t neighbour_index) const noexcept = 0;
virtual void precompute(const value_t t_end, const value_t dt, const std::vector<drum_t>& drum_vec) noexcept = 0;
virtual void step(value_t dt) noexcept = 0;
virtual value_t operator()(const size_t drum_index, const size_t neighbour_index) const noexcept = 0;
};
#endif
projects/braidingTightBinding/lib/diagonalizer.hpp
View file @ 5f273894
......@@ -18,74 +18,78 @@ extern "C" void dsyev_(
//TODO: Rule of 7 conformity (nontrivial heap resource)
class Diagonalizer{
public:
Diagonalizer() = default;
~Diagonalizer() = default;
std::vector<double> ev(const std::vector<double>& matrix, const size_t N){
//prepare memory
matrix_ = matrix;
if(eigenvalues_.size() != N){
eigenvalues_.clear();
eigenvalues_.reserve(N - eigenvalues_.capacity());
}
for(size_t i = 0; i < N; ++i)
eigenvalues_.push_back(0.);
//prepare workspace
N_ = N;
prepare_workspace('N');
dsyev_('N', 'L', N_, matrix_.data(), N_, eigenvalues_.data(), work_, lwork_, info_);
if(info_ != 0)
throw("Diagonalization failed!");
//clean up
delete[] work_;
return eigenvalues_;
}
std::pair<std::vector<double>, std::vector<double> > evv(const std::vector<double>& matrix, const size_t N){
//prepare memory
matrix_ = matrix;
if(eigenvalues_.size() != N){
eigenvalues_.clear();
eigenvalues_.reserve(N - eigenvalues_.capacity());
}
for(size_t i = 0; i < N; ++i)
eigenvalues_.push_back(0.);
//prepare workspace
N_ = N;
prepare_workspace('V');
dsyev_('V', 'L', N_, matrix_.data(), N_, eigenvalues_.data(), work_, lwork_, info_);
if(info_ != 0)
throw("Diagonalization failed!");
//clean up
delete[] work_;
return std::pair<std::vector<double>, std::vector<double> > (eigenvalues_, matrix_);
}
private:
void prepare_workspace(char type){
dsyev_(type, 'L', N_, matrix_.data(), N_, eigenvalues_.data(), &dwork_, -1, info_);
lwork_ = static_cast<int>(dwork_);
work_ = new double[lwork_];
}
int info_;
double dwork_;
int lwork_;
int N_; //linear matrix dimension
std::vector<double> matrix_; //the matrix
std::vector<double> eigenvalues_; //eigenvalues
double* work_;
public:
Diagonalizer() = default;
~Diagonalizer() = default;
std::vector<double> ev(const std::vector<double>& matrix, const size_t N){
//prepare memory
matrix_ = matrix;
if(eigenvalues_.size() != N){
eigenvalues_.clear();
eigenvalues_.reserve(N - eigenvalues_.capacity());
}
for(size_t i = 0; i < N; ++i){
eigenvalues_.push_back(0.);
}
//prepare workspace
N_ = N;
prepare_workspace('N');
dsyev_('N', 'L', N_, matrix_.data(), N_, eigenvalues_.data(), work_, lwork_, info_);
if(info_ != 0){
throw("Diagonalization failed!");
}
//clean up
delete[] work_;
return eigenvalues_;
}
std::pair<std::vector<double>, std::vector<double> > evv(const std::vector<double>& matrix, const size_t N){
//prepare memory
matrix_ = matrix;
if(eigenvalues_.size() != N){
eigenvalues_.clear();
eigenvalues_.reserve(N - eigenvalues_.capacity());
}
for(size_t i = 0; i < N; ++i){
eigenvalues_.push_back(0.);
}
//prepare workspace
N_ = N;
prepare_workspace('V');
dsyev_('V', 'L', N_, matrix_.data(), N_, eigenvalues_.data(), work_, lwork_, info_);
if(info_ != 0){
throw("Diagonalization failed!");
}
//clean up
delete[] work_;
return std::pair<std::vector<double>, std::vector<double> > (eigenvalues_, matrix_);
}
private:
void prepare_workspace(char type){
dsyev_(type, 'L', N_, matrix_.data(), N_, eigenvalues_.data(), &dwork_, -1, info_);
lwork_ = static_cast<int>(dwork_);
work_ = new double[lwork_];
}
int info_;
double dwork_;
int lwork_;
int N_; //linear matrix dimension
std::vector<double> matrix_; //the matrix
std::vector<double> eigenvalues_; //eigenvalues
double* work_;
};
#endif