diff --git a/projects/braidingTightBinding/bin/Frequencies.pdf b/projects/braidingTightBinding/bin/Frequencies.pdf
new file mode 100644
index 0000000000000000000000000000000000000000..e6be3f66ea472caa0f961a1bbeeac0a1869f54bd
Binary files /dev/null and b/projects/braidingTightBinding/bin/Frequencies.pdf differ
diff --git a/projects/braidingTightBinding/bin/Frequencies_zoomed.pdf b/projects/braidingTightBinding/bin/Frequencies_zoomed.pdf
new file mode 100644
index 0000000000000000000000000000000000000000..35dd75c002d035c662af1186e58c32ba65336d74
Binary files /dev/null and b/projects/braidingTightBinding/bin/Frequencies_zoomed.pdf differ
diff --git a/projects/braidingTightBinding/bin/Makefile b/projects/braidingTightBinding/bin/Makefile
new file mode 100644
index 0000000000000000000000000000000000000000..745727be5f5ea92e0afd5864043af03f3a2973c2
--- /dev/null
+++ b/projects/braidingTightBinding/bin/Makefile
@@ -0,0 +1,15 @@
+INCLUDEDIR = /Users/pascal/rbcomb_simulation/projects/braidingTightBinding/include
+LIBDIR = /Users/pascal/rbcomb_simulation/projects/braidingTightBinding/lib
+
+CXX = /usr/local/bin/g++-9
+CXXFLAGS = -std=c++2a -Wall -Wpedantic -I$(INCLUDEDIR) -I$(LIBDIR) -O5 -mtune=native -march=native -framework Accelerate -lblas
+
+run: all
+
+all: main_braiding
+
+main_braiding: main_braiding.cpp ${LIBDIR}/diagonalizer.hpp ${LIBDIR}/drum.hpp ${LIBDIR}/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
+ $(CXX) $< $(CXXFLAGS) -o $@
+
+clean:
+ rm -rf main_braiding
diff --git a/projects/braidingTightBinding/bin/ev.txt b/projects/braidingTightBinding/bin/ev.txt
new file mode 100644
index 0000000000000000000000000000000000000000..e403deecd75150bddb31082d47099534b4664771
--- /dev/null
+++ b/projects/braidingTightBinding/bin/ev.txt
@@ -0,0 +1,1860 @@
+0 21653.4
+1 21662
+2 21665.9
+3 21673
+4 21675.3
+5 21682.6
+6 21686.2
+7 21690.1
+8 21691.2
+9 21700.9
+10 21704.9
+11 21707.5
+12 21710.2
+13 21713.6
+14 21720.3
+15 21724.7
+16 21728.3
+17 21734.5
+18 21736
+19 21737.9
+20 21742.7
+21 21747.6
+22 21752.1
+23 21755.7
+24 21761.4
+25 21764.1
+26 21766.5
+27 21770
+28 21772.7
+29 21776.6
+30 21782.5
+31 21788.5
+32 21790.9
+33 21792.4
+34 21796.3
+35 21799.6
+36 21803.9
+37 21807.2
+38 21809.3
+39 21814.4
+40 21817.4
+41 21823.9
+42 21826.3
+43 21827.3
+44 21831
+45 21835.3
+46 21840.2
+47 21842.4
+48 21846.5
+49 21849.2
+50 21850.2
+51 21852.1
+52 21858.1
+53 21861.9
+54 21865.3
+55 21869.8
+56 21874.5
+57 21875.3
+58 21878.5
+59 21879.7
+60 21883.5
+61 21888
+62 21889.5
+63 21892.7
+64 21894.4
+65 21900.8
+66 21903.2
+67 21906.7
+68 21910
+69 21911.5
+70 21914
+71 21917
+72 21919
+73 21921.6
+74 21926
+75 21929.4
+76 21933.1
+77 21937.5
+78 21938.8
+79 21940.2
+80 21943.3
+81 21944.2
+82 21947.5
+83 21948.8
+84 21951.2
+85 21955.2
+86 21958
+87 21961.5
+88 21966.1
+89 21966.6
+90 21968.3
+91 21972.6
+92 21976
+93 21977.9
+94 21979
+95 21981.2
+96 21982.5
+97 21984
+98 21985.8
+99 21989.3
+100 21991.8
+101 21996.5
+102 21997
+103 21998.3
+104 22001.1
+105 22004.7
+106 22005.5
+107 22009.4
+108 22011.3
+109 22012.5
+110 22014.3
+111 22016.9
+112 22017.4
+113 22019.3
+114 22023
+115 22024.2
+116 22024.9
+117 22025.9
+118 22026.7
+119 22027.9
+120 22032.3
+121 22034.4
+122 22035.3
+123 22037.6
+124 22039.1
+125 22041.3
+126 22042.4
+127 22044.4
+128 22045.1
+129 22048.2
+130 22048.7
+131 22052
+132 22054.2
+133 22055.5
+134 22057.1
+135 22058.8
+136 22059.9
+137 22061.2
+138 22062
+139 22065.2
+140 22067.2
+141 22069.6
+142 22071.5
+143 22073.9
+144 22077.6
+145 22079.5
+146 22081.7
+147 22084.6
+148 22087.4
+149 22088.2
+150 22090.8
+151 22092.1
+152 22095.1
+153 22098.4
+154 22101.8
+155 22102.4
+156 22104.2
+157 22107.4
+158 22108.4
+159 22115.5
+160 22117
+161 22119.7
+162 22120.8
+163 22123.5
+164 22124.9
+165 22126
+166 22127.8
+167 22131.1
+168 22131.9
+169 22138.3
+170 22139.6
+171 22143.3
+172 22146.7
+173 22147.7
+174 22155
+175 22155.3
+176 22160.2
+177 22163.8
+178 22165.9
+179 22168.7
+180 22170.3
+181 22171.9
+182 22177.8
+183 22178.6
+184 22181.6
+185 22184.9
+186 22190.4
+187 22191.2
+188 22193.4
+189 22196.1
+190 22199.3
+191 22206.3
+192 22208.6
+193 22211.1
+194 22216.8
+195 22221
+196 22224.1
+197 22224.8
+198 22227.9
+199 22231.6
+200 22237.1
+201 22242.3
+202 22248.6
+203 22251.5
+204 22259.9
+205 22261.7
+206 22266.1
+207 22269.3
+208 22276.9
+209 22278.8
+210 22281.1
+211 22282.2
+212 22284.1
+213 22286.6
+214 22290.1
+215 22294.5
+216 22298.9
+217 22301.2
+218 22304
+219 22310.2
+220 22316.4
+221 22326
+222 22328.3
+223 22333.7
+224 22339.4
+225 22344.9
+226 22348.7
+227 22350.8
+228 22353.3
+229 22355.2
+230 22359.8
+231 22361.2
+232 22366
+233 22370.9
+234 22378.2
+235 22392.7
+236 22396.9
+237 22399.1
+238 22400.3
+239 22406.6
+240 22410
+241 22412.1
+242 22417.8
+243 22429.8
+244 22432.9
+245 22439.5
+246 22443.7
+247 22455.1
+248 22458.3
+249 22467.6
+250 22469.9
+251 22474
+252 22485.5
+253 22487
+254 22497.7
+255 22504.8
+256 22520.2
+257 22527.7
+258 22529.7
+259 22535.8
+260 22540.2
+261 22546.4
+262 22554.3
+263 22582.9
+264 22586
+265 22603.4
+266 22606.6
+267 22610.5
+268 22612.4
+269 22619.8
+270 22625.6
+271 22627.9
+272 22631.2
+273 22635.3
+274 22657.6
+275 22658.7
+276 22670.9
+277 22672.6
+278 22681.6
+279 22692.7
+280 22698.2
+281 22711.7
+282 22717.5
+283 22734.1
+284 22737.3
+285 22745.9
+286 22748.2
+287 22755.3
+288 22763.2
+289 22782.9
+290 22787
+291 22799.6
+292 22809.2
+293 22820.2
+294 22822.2
+295 22836.9
+296 22860.1
+297 22862.3
+298 22881
+299 22883.5
+300 22899.3
+301 22905.3
+302 22915
+303 22930.6
+304 22968.7
+305 22983.9
+306 22999.1
+307 23014.8
+308 23056.8
+309 23082
+310 23091
+311 23105.1
+312 23127.8
+313 23153.4
+314 23165.6
+315 23190.6
+316 23214.4
+317 23227.3
+318 23250.8
+319 23260.1
+320 23263.5
+321 23270
+322 23278.1
+323 23286.3
+324 23292.7
+325 23300.1
+326 23300.6
+327 23308.8
+328 23311.9
+329 23319.2
+330 23326.6
+331 23329.8
+332 23333.1
+333 23365.6
+334 23375.2
+335 23380.1
+336 23413.8
+337 23415.2
+338 23433.6
+339 23440.5
+340 23445.7
+341 23450.1
+342 23460.6
+343 23465.5
+344 23467.5
+345 23477.1
+346 23486.5
+347 23488.7
+348 23498.6
+349 23499.8
+350 23513.1
+351 23520.5
+352 23523.1
+353 23530.8
+354 23535.3
+355 23539.5
+356 23552.8
+357 23560.3
+358 23569.3
+359 23573.3
+360 23582.5
+361 23587
+362 23597.7
+363 23603.9
+364 23606.3
+365 23607.5
+366 23621.4
+367 23626.4
+368 23630.2
+369 23634.1
+370 23642.6
+371 23644
+372 23649.2
+373 23652.4
+374 23656.4
+375 23660.2
+376 23670
+377 23676.6
+378 23680.1
+379 23683.9
+380 23692.9
+381 23703.9
+382 23707.8
+383 23715.6
+384 23718.5
+385 23719.6
+386 23722.4
+387 23729.5
+388 23732.5
+389 23737.7
+390 23739.8
+391 23742.4
+392 23747.6
+393 23750.6
+394 23756.9
+395 23761.8
+396 23763.2
+397 23769.1
+398 23773.4
+399 23781.3
+400 23789.6
+401 23792.4
+402 23793.7
+403 23795.9
+404 23799.3
+405 23805.6
+406 23807.7
+407 23808.6
+408 23811.1
+409 23814.6
+410 23820.1
+411 23824.6
+412 23828.2
+413 23834.8
+414 23841.6
+415 23846.6
+416 23847.2
+417 23849.6
+418 23855.5
+419 23857.1
+420 23857.8
+421 23864.6
+422 23865.1
+423 23869.6
+424 23870.9
+425 23875.8
+426 23879
+427 23879.9
+428 23882.6
+429 23885.5
+430 23886
+431 23889.2
+432 23895.5
+433 23899.2
+434 23900.5
+435 23902.8
+436 23904.1
+437 23907.7
+438 23908.7
+439 23913.8
+440 23917.5
+441 23920.6
+442 23925.1
+443 23926.2
+444 23930.8
+445 23932.9
+446 23934
+447 23937.9
+448 23940.3
+449 23941.6
+450 23942.6
+451 23947.9
+452 23952.4
+453 23953.7
+454 23957.2
+455 23959.7
+456 23961.1
+457 23963.1
+458 23965.9
+459 23967.3
+460 23968.6
+461 23970.4
+462 23971.9
+463 23973.2
+464 23976
+465 23978.8
+466 23980
+467 23981
+468 23982.6
+469 23983.7
+470 23985.1
+471 23986.1
+472 23986.9
+473 23987.7
+474 23989.1
+475 23990.9
+476 23991.8
+477 23993.2
+478 23993.9
+479 23997.6
+480 23999.4
+481 24001.9
+482 24004
+483 24004.4
+484 24007.1
+485 24008.8
+486 24012.6
+487 24013.6
+488 24015
+489 24016.3
+490 24016.6
+491 24017.1
+492 24019.6
+493 24020
+494 24024.3
+495 24026.1
+496 24027.4
+497 24029.3
+498 24030.9
+499 24031.7
+500 24032.2
+501 24033.6
+502 24035.9
+503 24038.5
+504 24039
+505 24039.5
+506 24040.2
+507 24046
+508 24046.4
+509 24047.3
+510 24049.5
+511 24050.1
+512 24052.1
+513 24052.4
+514 24055.4
+515 24057.1
+516 24057.3
+517 24058.7
+518 24059
+519 24061.3
+520 24062.7
+521 24063.7
+522 24066.3
+523 24067.7
+524 24069.1
+525 24070.2
+526 24072.8
+527 24073.8
+528 24075.2
+529 24075.9
+530 24076.5
+531 24078
+532 24078.7
+533 24079.8
+534 24082.9
+535 24083.6
+536 24085.3
+537 24085.5
+538 24086.9
+539 24088.3
+540 24090.1
+541 24091.1
+542 24093
+543 24093.8
+544 24095.1
+545 24097.5
+546 24098.4
+547 24100.2
+548 24101.7
+549 24102.8
+550 24103.1
+551 24104.8
+552 24105.3
+553 24109
+554 24110.9
+555 24111.5
+556 24111.8
+557 24113.3
+558 24113.6
+559 24116
+560 24116.4
+561 24118.2
+562 24119.4
+563 24121.5
+564 24122
+565 24122.5
+566 24124.4
+567 24126.5
+568 24126.7
+569 24127.4
+570 24128.9
+571 24129.7
+572 24132.4
+573 24132.8
+574 24133.6
+575 24135.5
+576 24136.1
+577 24137.7
+578 24139.5
+579 24141.9
+580 24142.6
+581 24143.5
+582 24144.4
+583 24145.6
+584 24145.9
+585 24147.6
+586 24149.1
+587 24150.7
+588 24152.1
+589 24152.7
+590 24153.4
+591 24155.5
+592 24156.3
+593 24157.2
+594 24158.9
+595 24160.6
+596 24162.4
+597 24162.5
+598 24163.1
+599 24163.6
+600 24164.1
+601 24164.8
+602 24167.3
+603 24167.9
+604 24168.9
+605 24170.6
+606 24171.4
+607 24172.7
+608 24173.2
+609 24174.9
+610 24175.9
+611 24177
+612 24178.1
+613 24178.8
+614 24180
+615 24181.3
+616 24181.7
+617 24184.2
+618 24185.1
+619 24186.3
+620 24187.1
+621 24188.3
+622 24188.7
+623 24189.3
+624 24190.3
+625 24191
+626 24191.2
+627 24193
+628 24193.9
+629 24194.2
+630 24194.9
+631 24196.8
+632 24198
+633 24198.8
+634 24199.9
+635 24200.5
+636 24201
+637 24201.8
+638 24203.7
+639 24204.3
+640 24205.8
+641 24208.1
+642 24208.4
+643 24209.4
+644 24209.8
+645 24210.6
+646 24211.2
+647 24212
+648 24212.7
+649 24213.4
+650 24216.3
+651 24217
+652 24217.9
+653 24218.4
+654 24219.8
+655 24220.3
+656 24220.6
+657 24221.8
+658 24224.4
+659 24224.8
+660 24226
+661 24226.9
+662 24227.2
+663 24229
+664 24230.3
+665 24231.9
+666 24232
+667 24232.2
+668 24232.7
+669 24235
+670 24235.5
+671 24236
+672 24236.5
+673 24238.1
+674 24239.3
+675 24240.4
+676 24241.1
+677 24242
+678 24242.9
+679 24243.4
+680 24244.2
+681 24245.5
+682 24245.6
+683 24247.3
+684 24247.9
+685 24248.9
+686 24249.9
+687 24250.1
+688 24251.6
+689 24252.1
+690 24252.8
+691 24254.2
+692 24255.1
+693 24256.6
+694 24257.1
+695 24258.2
+696 24259.8
+697 24260.1
+698 24261.2
+699 24261.7
+700 24262.9
+701 24263.6
+702 24264.1
+703 24265.5
+704 24268.1
+705 24268.5
+706 24269.6
+707 24270.8
+708 24271.1
+709 24272
+710 24272.8
+711 24274.2
+712 24275.5
+713 24275.7
+714 24277.3
+715 24278.2
+716 24279.2
+717 24280
+718 24280.9
+719 24281.9
+720 24282.1
+721 24283.2
+722 24284.8
+723 24285.9
+724 24286.4
+725 24287.4
+726 24287.8
+727 24289.4
+728 24290.2
+729 24291.7
+730 24292.1
+731 24293.2
+732 24294.2
+733 24294.8
+734 24295.9
+735 24296.8
+736 24298.5
+737 24299.4
+738 24299.9
+739 24301.7
+740 24302.4
+741 24303.4
+742 24303.6
+743 24304.6
+744 24305.2
+745 24306.2
+746 24308
+747 24308.7
+748 24309
+749 24311.8
+750 24312.2
+751 24313.4
+752 24314.6
+753 24316.3
+754 24316.8
+755 24316.9
+756 24318.2
+757 24319
+758 24320.5
+759 24321.8
+760 24322
+761 24322.7
+762 24323.5
+763 24324.7
+764 24325.6
+765 24327.6
+766 24328.4
+767 24329.4
+768 24330.4
+769 24331.4
+770 24332.6
+771 24333
+772 24334.4
+773 24335.6
+774 24336.9
+775 24337.3
+776 24338.7
+777 24339.7
+778 24340.4
+779 24341.3
+780 24342.9
+781 24343.6
+782 24344.3
+783 24345.1
+784 24346.7
+785 24347.4
+786 24348.1
+787 24349.7
+788 24351.4
+789 24352.3
+790 24353.7
+791 24354.3
+792 24354.7
+793 24357.2
+794 24358.5
+795 24359.1
+796 24360.2
+797 24362.3
+798 24362.6
+799 24363.3
+800 24364.1
+801 24365.5
+802 24366.8
+803 24367.8
+804 24369.5
+805 24371.1
+806 24371.4
+807 24372.1
+808 24373.5
+809 24375.1
+810 24376.9
+811 24378.1
+812 24378.8
+813 24381.4
+814 24382.5
+815 24383
+816 24383.2
+817 24385.1
+818 24387.1
+819 24388.8
+820 24389.2
+821 24390.2
+822 24393.9
+823 24394.9
+824 24395.6
+825 24397.5
+826 24398.4
+827 24399.1
+828 24400.7
+829 24401.8
+830 24402.4
+831 24403.7
+832 24406.2
+833 24406.2
+834 24409.6
+835 24411.3
+836 24411.9
+837 24412.4
+838 24413.2
+839 24414.8
+840 24416.7
+841 24418.8
+842 24419.9
+843 24421.7
+844 24423.3
+845 24425.7
+846 24427.4
+847 24428.2
+848 24430.4
+849 24430.9
+850 24435.3
+851 24436.1
+852 24438.4
+853 24438.9
+854 24440.2
+855 24441.8
+856 24444.7
+857 24450.2
+858 24451.9
+859 24457.7
+860 24458.4
+861 24460.3
+862 24461.4
+863 24463.1
+864 24465.8
+865 24466.4
+866 24470.3
+867 24472.3
+868 24475.4
+869 24479.6
+870 24485.2
+871 24486.7
+872 24489.3
+873 24491.5
+874 24493.4
+875 24497.6
+876 24502.1
+877 24504
+878 24508.1
+879 24509.6
+880 24511.9
+881 24513.9
+882 24516.8
+883 24519.8
+884 24525.2
+885 24529
+886 24539
+887 24542.2
+888 24545.4
+889 24558.8
+890 24565
+891 24570.9
+892 24572.1
+893 24588.9
+894 24593.2
+895 24606.4
+896 24630.6
+897 24641.5
+898 24641.9
+899 24646.3
+900 24655
+901 24670.7
+902 24726.3
+903 24730.2
+904 24736.5
+905 24813.4
+906 24831.7
+907 24862.3
+908 24864.3
+909 24903.1
+910 24903.2
+911 25005.2
+912 25138.5
+913 25211.9
+914 25212.2
+915 26339.1
+916 26339.1
+917 26339.1
+918 26339.1
+919 26339.1
+920 26339.1
+921 26339.1
+922 26339.1
+923 26339.1
+924 26339.1
+925 26339.1
+926 26339.1
+927 26339.1
+928 26339.1
+929 26339.1
+930 26339.1
+931 26339.1
+932 26339.1
+933 26339.1
+934 26339.1
+935 26339.1
+936 26339.1
+937 26339.1
+938 26339.1
+939 26339.1
+940 26339.1
+941 26339.1
+942 26339.1
+943 26339.1
+944 26339.1
+945 27419.8
+946 27420.1
+947 27487.3
+948 27608.7
+949 27700.7
+950 27700.8
+951 27735.7
+952 27737.5
+953 27764.8
+954 27781.2
+955 27849.7
+956 27855.3
+957 27858.7
+958 27908
+959 27921.9
+960 27929.5
+961 27933.5
+962 27933.8
+963 27943.4
+964 27964.7
+965 27976.3
+966 27980.1
+967 27994.8
+968 27995.9
+969 28001.1
+970 28006.5
+971 28018.3
+972 28021
+973 28023.8
+974 28032.7
+975 28035.9
+976 28040.7
+977 28043.3
+978 28045.8
+979 28047.6
+980 28049.6
+981 28050.9
+982 28054.5
+983 28056.2
+984 28060.1
+985 28063.8
+986 28065.4
+987 28067.3
+988 28069.6
+989 28070.9
+990 28075.8
+991 28079.5
+992 28082.1
+993 28083.9
+994 28087.3
+995 28087.8
+996 28090.1
+997 28091.6
+998 28092.6
+999 28094.3
+1000 28094.8
+1001 28099.9
+1002 28101.4
+1003 28106.2
+1004 28108.7
+1005 28110.1
+1006 28111.2
+1007 28111.6
+1008 28113.6
+1009 28114.3
+1010 28118.2
+1011 28118.6
+1012 28120.5
+1013 28121.2
+1014 28122.7
+1015 28124.7
+1016 28126.2
+1017 28127.7
+1018 28128.7
+1019 28130.5
+1020 28132.1
+1021 28133.5
+1022 28134.2
+1023 28134.7
+1024 28135.2
+1025 28136.7
+1026 28139.6
+1027 28139.6
+1028 28141.7
+1029 28142.9
+1030 28143.5
+1031 28144.4
+1032 28145.8
+1033 28146.4
+1034 28147.1
+1035 28148.8
+1036 28149.4
+1037 28150.3
+1038 28153.5
+1039 28154.3
+1040 28154.7
+1041 28156.1
+1042 28157.9
+1043 28159.5
+1044 28159.7
+1045 28160.2
+1046 28161.1
+1047 28163.3
+1048 28164
+1049 28165
+1050 28166.5
+1051 28168
+1052 28169.1
+1053 28169.7
+1054 28170
+1055 28171.4
+1056 28172.9
+1057 28173.7
+1058 28174.9
+1059 28176.1
+1060 28176.8
+1061 28177.3
+1062 28177.6
+1063 28179.5
+1064 28180.4
+1065 28180.9
+1066 28182
+1067 28184.2
+1068 28184.5
+1069 28185
+1070 28186.2
+1071 28187
+1072 28188.5
+1073 28189.9
+1074 28190.5
+1075 28191.1
+1076 28192.5
+1077 28193.2
+1078 28193.8
+1079 28194.4
+1080 28195.8
+1081 28196.5
+1082 28197.2
+1083 28198
+1084 28199.2
+1085 28199.6
+1086 28200.7
+1087 28201.7
+1088 28202.9
+1089 28203.3
+1090 28204.3
+1091 28205.1
+1092 28206
+1093 28206.9
+1094 28207.6
+1095 28209.3
+1096 28210.1
+1097 28211.1
+1098 28211.8
+1099 28212.4
+1100 28212.6
+1101 28213.7
+1102 28215
+1103 28215.7
+1104 28216.8
+1105 28216.9
+1106 28217.4
+1107 28218.8
+1108 28219.9
+1109 28220.9
+1110 28221.2
+1111 28223.6
+1112 28223.9
+1113 28224.5
+1114 28226
+1115 28226.9
+1116 28227.4
+1117 28228.2
+1118 28228.4
+1119 28229.3
+1120 28229.9
+1121 28231.5
+1122 28231.9
+1123 28232.7
+1124 28234.1
+1125 28234.9
+1126 28235.8
+1127 28236.3
+1128 28237.2
+1129 28238.1
+1130 28238.5
+1131 28239.8
+1132 28240.5
+1133 28241.9
+1134 28242.2
+1135 28243.1
+1136 28243.5
+1137 28244.5
+1138 28245.8
+1139 28246.8
+1140 28247
+1141 28247.8
+1142 28248.5
+1143 28249.2
+1144 28250.1
+1145 28250.9
+1146 28252.3
+1147 28252.4
+1148 28253.6
+1149 28254.8
+1150 28255.4
+1151 28256.2
+1152 28256.5
+1153 28257.5
+1154 28258.5
+1155 28258.8
+1156 28261
+1157 28262.2
+1158 28262.7
+1159 28263.2
+1160 28264.3
+1161 28264.8
+1162 28265.7
+1163 28265.9
+1164 28267.3
+1165 28268.2
+1166 28268.6
+1167 28270
+1168 28270.7
+1169 28272
+1170 28272.5
+1171 28273
+1172 28274.3
+1173 28274.4
+1174 28275.3
+1175 28276.1
+1176 28276.6
+1177 28278.1
+1178 28278.2
+1179 28279.3
+1180 28280
+1181 28280.4
+1182 28281.2
+1183 28282
+1184 28282.6
+1185 28283.5
+1186 28284.5
+1187 28285.9
+1188 28286.3
+1189 28286.8
+1190 28287.2
+1191 28289.2
+1192 28289.6
+1193 28289.7
+1194 28289.8
+1195 28291.2
+1196 28292.3
+1197 28293.9
+1198 28294.1
+1199 28294.9
+1200 28295.9
+1201 28296.3
+1202 28298.5
+1203 28299.5
+1204 28299.7
+1205 28300.2
+1206 28301.4
+1207 28301.8
+1208 28302.6
+1209 28303.2
+1210 28305.7
+1211 28306.2
+1212 28306.9
+1213 28307.5
+1214 28308
+1215 28308.8
+1216 28309.1
+1217 28309.9
+1218 28310.2
+1219 28312.1
+1220 28313.5
+1221 28314
+1222 28315.6
+1223 28316.2
+1224 28316.7
+1225 28317.2
+1226 28318.2
+1227 28318.9
+1228 28319.9
+1229 28321.5
+1230 28322.1
+1231 28322.3
+1232 28323.1
+1233 28324.7
+1234 28324.8
+1235 28325.4
+1236 28326.2
+1237 28326.7
+1238 28327.1
+1239 28328.2
+1240 28328.9
+1241 28329.9
+1242 28330.6
+1243 28332.7
+1244 28333.1
+1245 28334.2
+1246 28335.2
+1247 28335.8
+1248 28336.8
+1249 28337.7
+1250 28338.6
+1251 28340
+1252 28340.4
+1253 28341.6
+1254 28342.3
+1255 28343.6
+1256 28344.5
+1257 28345
+1258 28347.2
+1259 28347.7
+1260 28348.2
+1261 28348.7
+1262 28349.1
+1263 28349.3
+1264 28350.7
+1265 28352.2
+1266 28353.7
+1267 28354.4
+1268 28355.1
+1269 28356.9
+1270 28357.4
+1271 28358
+1272 28359.2
+1273 28360.5
+1274 28361.8
+1275 28363.2
+1276 28363.6
+1277 28364.6
+1278 28365.3
+1279 28366.1
+1280 28366.7
+1281 28368.7
+1282 28370.3
+1283 28371.6
+1284 28372.1
+1285 28373.8
+1286 28374.4
+1287 28374.7
+1288 28377.1
+1289 28377.8
+1290 28379
+1291 28379.6
+1292 28379.8
+1293 28381.5
+1294 28383.2
+1295 28383.6
+1296 28384
+1297 28385.8
+1298 28386.8
+1299 28388.3
+1300 28388.7
+1301 28390.8
+1302 28391
+1303 28392.2
+1304 28392.6
+1305 28393
+1306 28394.7
+1307 28397.8
+1308 28398.2
+1309 28399.6
+1310 28399.9
+1311 28400.8
+1312 28402.1
+1313 28403.7
+1314 28404.4
+1315 28406.5
+1316 28407.6
+1317 28408.2
+1318 28409.9
+1319 28410.7
+1320 28412.2
+1321 28413.4
+1322 28414.6
+1323 28414.8
+1324 28416.2
+1325 28416.8
+1326 28419.4
+1327 28420.3
+1328 28421
+1329 28422.2
+1330 28422.7
+1331 28423.4
+1332 28424.5
+1333 28425.4
+1334 28427.5
+1335 28428.4
+1336 28429.6
+1337 28430.8
+1338 28433
+1339 28433.9
+1340 28435.1
+1341 28437.1
+1342 28437.3
+1343 28438.4
+1344 28438.7
+1345 28440
+1346 28442.6
+1347 28442.9
+1348 28444.6
+1349 28445.1
+1350 28446.9
+1351 28447.7
+1352 28448
+1353 28452.9
+1354 28453.6
+1355 28454
+1356 28454.3
+1357 28456.6
+1358 28458.5
+1359 28459.7
+1360 28460.1
+1361 28460.8
+1362 28462.1
+1363 28463.8
+1364 28464.8
+1365 28466.3
+1366 28469.9
+1367 28470.3
+1368 28472.4
+1369 28472.8
+1370 28473.1
+1371 28474.2
+1372 28475.3
+1373 28476.2
+1374 28479.4
+1375 28480.9
+1376 28483.1
+1377 28483.5
+1378 28485.2
+1379 28487.3
+1380 28488.9
+1381 28492
+1382 28492.6
+1383 28493.8
+1384 28494.5
+1385 28496
+1386 28497.2
+1387 28497.9
+1388 28498.6
+1389 28499.4
+1390 28500.5
+1391 28501.5
+1392 28502.8
+1393 28503.7
+1394 28504.7
+1395 28507
+1396 28509.4
+1397 28510.5
+1398 28511.8
+1399 28513.2
+1400 28514.3
+1401 28515.5
+1402 28517.9
+1403 28519.6
+1404 28520.8
+1405 28522.8
+1406 28525.8
+1407 28526.9
+1408 28530.7
+1409 28535.1
+1410 28535.9
+1411 28537
+1412 28539.1
+1413 28542.4
+1414 28543.2
+1415 28545
+1416 28548.8
+1417 28549.8
+1418 28553.6
+1419 28556.1
+1420 28559.3
+1421 28563.5
+1422 28564.3
+1423 28567.3
+1424 28568.4
+1425 28570.4
+1426 28571.5
+1427 28574.6
+1428 28579.8
+1429 28582.5
+1430 28582.9
+1431 28585.3
+1432 28587.6
+1433 28588.3
+1434 28591
+1435 28595.1
+1436 28596.2
+1437 28599.9
+1438 28600.3
+1439 28606
+1440 28606.6
+1441 28607.9
+1442 28612.8
+1443 28614.9
+1444 28615.4
+1445 28619.5
+1446 28625.2
+1447 28630.7
+1448 28633.7
+1449 28637.5
+1450 28642
+1451 28644.9
+1452 28647
+1453 28647.7
+1454 28649.5
+1455 28654.7
+1456 28657.5
+1457 28659.4
+1458 28660.4
+1459 28662.8
+1460 28669.6
+1461 28676.2
+1462 28679.8
+1463 28684.6
+1464 28685.8
+1465 28689.9
+1466 28695.1
+1467 28697.6
+1468 28701.9
+1469 28704.1
+1470 28705.8
+1471 28710.1
+1472 28712.5
+1473 28718.4
+1474 28720.7
+1475 28721.6
+1476 28724.1
+1477 28730.5
+1478 28733.7
+1479 28742.8
+1480 28750.2
+1481 28753.3
+1482 28756.2
+1483 28761.6
+1484 28769.7
+1485 28772.8
+1486 28776.1
+1487 28778.8
+1488 28783
+1489 28784.2
+1490 28791.1
+1491 28794.3
+1492 28797.5
+1493 28801.6
+1494 28813
+1495 28813.9
+1496 28815.9
+1497 28821
+1498 28829.7
+1499 28833.4
+1500 28841
+1501 28844.2
+1502 28851.6
+1503 28857.7
+1504 28868.5
+1505 28871.9
+1506 28875.6
+1507 28881.9
+1508 28884
+1509 28890.1
+1510 28900.8
+1511 28901.9
+1512 28909.9
+1513 28911.7
+1514 28919.3
+1515 28927.1
+1516 28928.7
+1517 28932.7
+1518 28941.2
+1519 28944.8
+1520 28949
+1521 28954.6
+1522 28969.4
+1523 28970.6
+1524 28997.8
+1525 29001.7
+1526 29009.5
+1527 29035.6
+1528 29038.2
+1529 29040.9
+1530 29046.8
+1531 29052.6
+1532 29055.1
+1533 29061.7
+1534 29062.1
+1535 29068
+1536 29073.2
+1537 29079.7
+1538 29086.2
+1539 29091.4
+1540 29094.1
+1541 29101.5
+1542 29120.3
+1543 29130.6
+1544 29149.6
+1545 29169.4
+1546 29179.1
+1547 29199.4
+1548 29217.4
+1549 29228.5
+1550 29235.6
+1551 29255.5
+1552 29288.6
+1553 29300.9
+1554 29312.8
+1555 29324.7
+1556 29354.5
+1557 29366.7
+1558 29374.3
+1559 29379
+1560 29391.2
+1561 29393.2
+1562 29407.7
+1563 29409.4
+1564 29427.5
+1565 29438.8
+1566 29440.5
+1567 29449
+1568 29456.4
+1569 29466.1
+1570 29469.3
+1571 29484.5
+1572 29490.6
+1573 29496.1
+1574 29497.9
+1575 29504.5
+1576 29506.9
+1577 29519.8
+1578 29524.2
+1579 29534.6
+1580 29538.8
+1581 29547.4
+1582 29554.2
+1583 29555.5
+1584 29564.9
+1585 29565.7
+1586 29582.8
+1587 29585.9
+1588 29588.5
+1589 29590.2
+1590 29594.7
+1591 29600.3
+1592 29601.8
+1593 29604.7
+1594 29607.2
+1595 29620.5
+1596 29622.9
+1597 29644.6
+1598 29650.6
+1599 29655.3
+1600 29658.7
+1601 29663.3
+1602 29664.8
+1603 29670.5
+1604 29682.3
+1605 29687.6
+1606 29695.7
+1607 29696.8
+1608 29705.5
+1609 29708.6
+1610 29710.4
+1611 29717.4
+1612 29719.8
+1613 29728.4
+1614 29731.6
+1615 29736.6
+1616 29738.9
+1617 29747.9
+1618 29752.3
+1619 29753.9
+1620 29756.4
+1621 29761.1
+1622 29762
+1623 29763.8
+1624 29766.9
+1625 29777.8
+1626 29783.3
+1627 29786.9
+1628 29790.5
+1629 29791.6
+1630 29795
+1631 29796.5
+1632 29798.4
+1633 29799.9
+1634 29802.8
+1635 29806.9
+1636 29811.1
+1637 29815.2
+1638 29816.9
+1639 29824.1
+1640 29828.8
+1641 29833.4
+1642 29835.5
+1643 29837.2
+1644 29840.5
+1645 29843.8
+1646 29846.4
+1647 29848.3
+1648 29849.7
+1649 29850.5
+1650 29852.3
+1651 29853.7
+1652 29859.3
+1653 29861.7
+1654 29865
+1655 29866.3
+1656 29872.5
+1657 29874.7
+1658 29879.4
+1659 29883.3
+1660 29887.4
+1661 29890.1
+1662 29892.5
+1663 29892.9
+1664 29895.3
+1665 29898.4
+1666 29902.7
+1667 29904.5
+1668 29906.2
+1669 29911.4
+1670 29913.8
+1671 29915.8
+1672 29917.4
+1673 29918
+1674 29922
+1675 29924.5
+1676 29926.8
+1677 29927.3
+1678 29931.7
+1679 29932.9
+1680 29934.1
+1681 29936.2
+1682 29937.7
+1683 29940.4
+1684 29944
+1685 29944.2
+1686 29949.6
+1687 29950.3
+1688 29952.9
+1689 29955.6
+1690 29956.6
+1691 29961.3
+1692 29961.9
+1693 29964.3
+1694 29965.6
+1695 29966.5
+1696 29967.5
+1697 29969.5
+1698 29970.3
+1699 29972.3
+1700 29973.4
+1701 29978.6
+1702 29979.4
+1703 29981.7
+1704 29983.1
+1705 29983.5
+1706 29986
+1707 29988.4
+1708 29990.6
+1709 29991.6
+1710 29993.6
+1711 29994.1
+1712 29996.1
+1713 29998.3
+1714 29999.9
+1715 30001.3
+1716 30004
+1717 30005.8
+1718 30007.2
+1719 30009
+1720 30010.4
+1721 30012.8
+1722 30013.4
+1723 30014.4
+1724 30015.1
+1725 30016.4
+1726 30017.6
+1727 30018.6
+1728 30020.1
+1729 30022.6
+1730 30023
+1731 30025.2
+1732 30025.7
+1733 30027.2
+1734 30028
+1735 30029.6
+1736 30030.7
+1737 30032.4
+1738 30033
+1739 30034.6
+1740 30037.8
+1741 30038.7
+1742 30039.3
+1743 30040
+1744 30040.6
+1745 30041.4
+1746 30044.1
+1747 30045.6
+1748 30045.9
+1749 30047.8
+1750 30049.2
+1751 30050
+1752 30051.4
+1753 30054.3
+1754 30054.8
+1755 30057.5
+1756 30059.5
+1757 30060.4
+1758 30060.8
+1759 30064.3
+1760 30066.1
+1761 30068.7
+1762 30070
+1763 30071.1
+1764 30072
+1765 30073.6
+1766 30074.4
+1767 30075.8
+1768 30078.3
+1769 30081.4
+1770 30082.7
+1771 30083
+1772 30086.4
+1773 30088.9
+1774 30091
+1775 30094
+1776 30095.7
+1777 30096.6
+1778 30099
+1779 30099.7
+1780 30102
+1781 30103
+1782 30103.9
+1783 30107.1
+1784 30109.8
+1785 30112.3
+1786 30115.5
+1787 30117.4
+1788 30118.9
+1789 30121.1
+1790 30122.8
+1791 30124
+1792 30126.3
+1793 30128.9
+1794 30130.6
+1795 30135.3
+1796 30136.6
+1797 30138.8
+1798 30140
+1799 30143.2
+1800 30146
+1801 30146.9
+1802 30149.2
+1803 30149.7
+1804 30153.1
+1805 30156.4
+1806 30158.9
+1807 30161.6
+1808 30166
+1809 30167.4
+1810 30168.1
+1811 30170
+1812 30173
+1813 30174.6
+1814 30178.2
+1815 30181.2
+1816 30183.9
+1817 30184.7
+1818 30186.4
+1819 30191
+1820 30193.3
+1821 30196.9
+1822 30198.4
+1823 30200.8
+1824 30204
+1825 30206.3
+1826 30209.1
+1827 30210.2
+1828 30212
+1829 30216.2
+1830 30220.5
+1831 30223.4
+1832 30225.3
+1833 30227.8
+1834 30229.5
+1835 30231.5
+1836 30235.5
+1837 30238.2
+1838 30241.4
+1839 30244.9
+1840 30248.4
+1841 30249.7
+1842 30250.8
+1843 30255.3
+1844 30257.9
+1845 30261
+1846 30265.8
+1847 30268.2
+1848 30270.2
+1849 30272
+1850 30274.9
+1851 30281.9
+1852 30282.7
+1853 30285.5
+1854 30288.1
+1855 30293.3
+1856 30294.9
+1857 30300
+1858 30302.8
+1859 30308.9
diff --git a/projects/braidingTightBinding/bin/main_braiding b/projects/braidingTightBinding/bin/main_braiding
new file mode 100755
index 0000000000000000000000000000000000000000..0c7b9140900c8d3e643fccd497a1ec1bfb96f37c
Binary files /dev/null and b/projects/braidingTightBinding/bin/main_braiding differ
diff --git a/projects/braidingTightBinding/bin/main_braiding.cpp b/projects/braidingTightBinding/bin/main_braiding.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..357332608e333530e4f5ebc00c0b8bdb7876d111
--- /dev/null
+++ b/projects/braidingTightBinding/bin/main_braiding.cpp
@@ -0,0 +1,86 @@
+#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_braid.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>
+
+int main(){
+ //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 = 27388152213.022964;
+ value_t k1 = 493480220.0544679;
+ value_t k2 = 2467401100.2723393;
+ value_t c = 5.235987755982988;
+
+ vortex_t vortex (vec2_t(17.,25.), 5., 1., 1., 1.); //as Eliska did it
+ std::vector<vortex_t> vortices ({vortex});
+
+ params_t drum_parameters (k0, k1, k2, c, vec2_t (0., 0.), 'A');
+ generator_t lattice_generator (30,15);
+ auto lat_pair = lattice_generator(drum_parameters);
+
+ //fix B site
+ for(auto d: lat_pair.first)
+ if(d.get_parameters().sublattice == 'B')
+ d.get_parameters().k0 = 21959869792.42382;
+
+ 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 = "PARAMS.TXT";
+ std::string adjacency_file = "ADJACENCY.TXT";
+ std::string dynamics_file = "DYNAMICS.TXT";
+ grabber_t grabber (save_every, parameters_file, adjacency_file, dynamics_file);
+
+ system_t system (1., 0.00001, lat_pair.first, stepper, force, system_parameters, grabber);
+
+ //set up matrix element calculator
+ matelecalc_t mec;
+
+ //set up diagonalizer
+ Diagonalizer diagonalizer;
+ //calculate eigenvalues
+ std::cout << "Getting matrix\n";
+ std::vector<value_t> matrix = system.get_matrix(mec);
+ std::cout << "Diaging\n";
+ std::vector<value_t> evals = diagonalizer.ev(system.get_matrix(mec), lat_pair.first.size()-1);
+
+ for(size_t i = 0; i < evals.size(); ++i)
+ std::cout << i << " \t" << std::sqrt(evals[i])/(2.*M_PI) << std::endl;
+
+
+}
diff --git a/projects/braidingTightBinding/bin/plot_ev.py b/projects/braidingTightBinding/bin/plot_ev.py
new file mode 100644
index 0000000000000000000000000000000000000000..55b694a6539ed47f4f1f5e07774b2f01d2a3dddb
--- /dev/null
+++ b/projects/braidingTightBinding/bin/plot_ev.py
@@ -0,0 +1,17 @@
+import matplotlib.pyplot as plt
+import numpy as np
+
+data = np.loadtxt('ev.txt')
+
+plt.plot(data[900:960,0],data[900:960,1],'o', ms=1)
+plt.title("Frequency levels under vortex Kekule")
+plt.xlabel("Level index")
+plt.ylabel("Frequency [Hz]")
+plt.savefig("Frequencies_zoomed.pdf", dpi = 500)
+
+plt.clf()
+plt.plot(data[:,0],data[:,1],'o', ms=1)
+plt.title("Frequency levels under vortex Kekule")
+plt.xlabel("Level index")
+plt.ylabel("Frequency [Hz]")
+plt.savefig("Frequencies.pdf", dpi = 500)
diff --git a/projects/braidingTightBinding/include/coupler_braid.hpp b/projects/braidingTightBinding/include/coupler_braid.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..9408009433277fc0481303b0fc4b4b52702d0bca
--- /dev/null
+++ b/projects/braidingTightBinding/include/coupler_braid.hpp
@@ -0,0 +1,57 @@
+#ifndef COUPLER_BRAID_HPP_INCLUDED
+#define COUPLER_BRAID_HPP_INCLUDED
+#include <vector>
+#include <vortex.hpp>
+
+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;
+
+ value_t ret_val = 1.; //vortices are multiplicative
+ for(auto tex: vortices_){
+ ret_val *= tex(v, s);
+ }
+
+ return 1. + 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
diff --git a/projects/braidingTightBinding/include/driver_braid.hpp b/projects/braidingTightBinding/include/driver_braid.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..e36a2adf4367a7224cb64f14767b5e9f33a87c3d
--- /dev/null
+++ b/projects/braidingTightBinding/include/driver_braid.hpp
@@ -0,0 +1,27 @@
+#ifndef DRIVER_BRAID_HPP_INCLUDED
+#define DRIVER_BRAID_HPP_INCLUDED
+#include <vector>
+#include <drum.hpp>
+
+template <typename value_t, typename drum_t>
+class DriverBraid: public Driver<value_t, drum_t>{
+ 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 step(value_t dt) noexcept final
+ {
+
+ }
+ value_t operator()(const size_t drum_index) const noexcept final override
+ {
+ return 0;
+ }
+};
+
+#endif
diff --git a/projects/braidingTightBinding/include/drum_parameters_braid.hpp b/projects/braidingTightBinding/include/drum_parameters_braid.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..7211663064acbbe15b6d87708e93c596bbedc989
--- /dev/null
+++ b/projects/braidingTightBinding/include/drum_parameters_braid.hpp
@@ -0,0 +1,65 @@
+#ifndef DRUM_PARAMETERS_BRAID_HPP_INCLUDED
+#define DRUM_PARAMETERS_BRAID_HPP_INCLUDED
+#include <vec2.hpp>
+#include <cmath>
+#include <iostream>
+
+//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;
+
+ 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
+
+};
+
+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;
+
+}
+
+#endif
diff --git a/projects/braidingTightBinding/include/drum_variables_braid.hpp b/projects/braidingTightBinding/include/drum_variables_braid.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..d8a7170ef82ca03abeed3fa3d8bbfc8032146e73
--- /dev/null
+++ b/projects/braidingTightBinding/include/drum_variables_braid.hpp
@@ -0,0 +1,23 @@
+#ifndef DRUM_VARIABLES_BRAID_HPP_INCLUDED
+#define DRUM_VARIABLES_BRAID_HPP_INCLUDED
+
+/*
+ * Stores all dynamic variables of a drum
+ */
+
+//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;
+
+ 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
diff --git a/projects/braidingTightBinding/include/force_braid.hpp b/projects/braidingTightBinding/include/force_braid.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..62748139d9021ddc31906ff8fcc7273a3e4c9133
--- /dev/null
+++ b/projects/braidingTightBinding/include/force_braid.hpp
@@ -0,0 +1,44 @@
+#ifndef FORCE_BRAID_HPP_INCLUDED
+#define FORCE_BRAID_HPP_INDLUDED
+#include <force.hpp>
+#include <drum_parameters_braid.hpp>
+#include <drum_variables_braid.hpp>
+#include <drum.hpp>
+
+//This force is to be used with params_t = DrumParametersBraid, vars_t = DrumVariablesBraid.
+
+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;
+
+ 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;
+
+ return part1 + part2 + part3;
+ }
+};
+
+
+#endif
diff --git a/projects/braidingTightBinding/include/matrix_element_calculator_braid.hpp b/projects/braidingTightBinding/include/matrix_element_calculator_braid.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..e09455181531bc4510ea629624333421502f8c20
--- /dev/null
+++ b/projects/braidingTightBinding/include/matrix_element_calculator_braid.hpp
@@ -0,0 +1,35 @@
+#ifndef MATRIX_ELEMENT_CALCULATOR_BRAID_HPP_INCLUDED
+#define MATRIX_ELEMENT_CALCULATOR_BRAID_HPP_INCLUDED
+#include <matrix_element_calculator.hpp>
+
+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;
+
+ //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;
+ }
+};
+
+#endif
diff --git a/projects/braidingTightBinding/include/rbcomb_generator_braid.hpp b/projects/braidingTightBinding/include/rbcomb_generator_braid.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..761ec39d7e2ea76414f666d170d0ace8f7159661
--- /dev/null
+++ b/projects/braidingTightBinding/include/rbcomb_generator_braid.hpp
@@ -0,0 +1,301 @@
+#ifndef RBCOMB_GENERATOR_BRAID_HPP_INCLUDED
+#define RBCOMB_GENERATOR_BRAID_HPP_INCLUDED
+#include <lattice_generator.hpp>
+#include <vec2.hpp>
+#include <utility>
+#include <vector>
+
+template <typename value_t, typename params_t, typename vars_t, typename sbuffer_t>
+class RbcombGeneratorBraid: public LatticeGenerator<value_t, params_t, vars_t, sbuffer_t> {
+ public:
+ RbcombGeneratorBraid(size_t x_extent, size_t y_extent) noexcept : x_extent_(x_extent), y_extent_(y_extent) {}
+ RbcombGeneratorBraid() = delete;
+ ~RbcombGeneratorBraid() = default;
+
+ //the params_t argument shall contain the position of the first drum that is placed.
+ //the last drum (return.first.back()) signifies inexistent neighbours, it is not part of the lattice.
+ //It shall be of sublattice 'A'.
+ std::pair<std::vector<Drum<value_t, params_t, vars_t, sbuffer_t> >, std::vector<std::vector<int> > >
+ operator()(const params_t& drum_params) noexcept final override
+ {
+ using drum_t = Drum<value_t, params_t, vars_t, sbuffer_t>;
+ std::vector<drum_t> current_line, next_line;
+ current_line.reserve(x_extent_);
+ next_line.reserve(x_extent_);
+
+ //calculate vectors to next drums
+ Vec2<value_t> vec_right = drum_params.s2 - drum_params.s1; //drum to the right in same row
+ Vec2<value_t> vec_down = drum_params.s0; //drum below
+ Vec2<value_t> vec_upleft = drum_params.s1; //drum at top left
+ Vec2<value_t> vec_upright = drum_params.s2; //drum at top right
+
+ //prepare memory
+ drums_.reserve(x_extent_);
+ adjacency_vector_.reserve(x_extent_);
+ //push first drum
+ drums_.push_back(drum_t(drum_params));
+ adjacency_vector_.push_back(std::vector<int>(3,-1)); //initialized with 'no neighbour here'
+
+ //initialize the first line to bootstrap the algorithm
+ for(size_t i = 1; i < x_extent_; ++i){
+ params_t new_params (drum_params);
+ new_params.position = drums_.back().get_parameters().position + vec_right;
+ drums_.push_back(drum_t(new_params));
+ adjacency_vector_.push_back(std::vector<int>(3,-1));
+ }
+
+ //add y_extent_ inequivalent hexagon rows
+ for(size_t i = 0; i < y_extent_; ++i){
+ generate_layer(vec_down, vec_upleft, vec_upright);
+ }
+
+ //push a drum to the end that has no neighbours and signifies the inexistent neighbour
+ drums_.push_back(drum_t(drum_params));
+ adjacency_vector_.push_back(std::vector<int>(3,-1));
+ //make no-neighbours point here
+ for(size_t i = 0; i < adjacency_vector_.size()-1; ++i){
+ for(size_t j = 0; j < 3; ++j){
+ if(adjacency_vector_[i][j] == -1){
+ adjacency_vector_[i][j] = drums_.size()-1;
+ }
+ }
+ }
+
+ /* Deploy without consistency check
+ if(check_consistency() == false)
+ std::cout << "INCONSISTENT LATTICE DETECTED!\n";
+ */
+
+ return std::pair<std::vector<drum_t>, std::vector<std::vector<int> > > (drums_, adjacency_vector_);
+ }
+
+ private:
+ //check if there exists a drum in drum_vec that is close to drum_pos.
+ //returns the index of the first close drum or -1.
+ int find_close(
+ const Vec2<value_t>& drum_pos,
+ const std::vector<Drum<value_t, params_t, vars_t, sbuffer_t> >& drum_vec) const noexcept
+ {
+ for(int i = 0; i < drum_vec.size(); ++i){
+ auto test_drum = drum_vec[i];
+ if(test_drum.get_parameters().position.r_wrt(drum_pos) < 1.e-8)
+ return i;
+ }
+ return -1; //no close drum found
+ }
+
+ //checks sanity of generated lattice
+ //returns true if lattice passed the test
+ bool check_consistency() noexcept
+ {
+ //check for duplicate drums
+ for(size_t i = 0; i < drums_.size(); ++i){
+ if(i != find_close(drums_[i].get_parameters().position, drums_))
+ return false;
+ }
+ return true;
+ }
+
+ //arguments are the neighbour vectors down, upleft, upright.
+ //adds a new layer of drums to drums_, filling in the adjency_vector_ as well.
+ //
+ // 4: | x_extent_ drums
+ // 3: \/ x_extent_ + 1 drums
+ // 2: | x_extent_ + 1 drums
+ // 1: \/ x_extent_ drums
+ void generate_layer(const Vec2<value_t>& v_d, const Vec2<value_t>& v_ul, const Vec2<value_t>& v_ur) noexcept
+ {
+ //shorten syntax
+ using drum_t = Drum<value_t, params_t, vars_t, sbuffer_t>;
+
+ //prepare vector to hold new row
+ std::vector<drum_t> new_row;
+ new_row.reserve(x_extent_);
+
+ //prepare parameter templates for both sublattices
+ params_t params_A(drums_[0].get_parameters()); //first drum will always be of sublattice 'A'
+ params_t params_B(params_A.k0, params_A.k1, params_A.k2, params_A.c, params_A.position, 'B'); //generate sublattice 'B' type of parameters
+
+ //row 1
+ //prepare iterator
+ auto drum_it = drums_.end();
+ drum_it -= x_extent_;
+
+ while(drum_it != drums_.end()){
+ Vec2<value_t> new_pos1 = drum_it->get_parameters().position + v_ul;
+ Vec2<value_t> new_pos2 = drum_it->get_parameters().position + v_ur;
+
+ //see if these drums already exist and add them to the list
+ int index_1 = find_close(new_pos1, new_row);
+ if(index_1 == -1){ //not in there yet, so push it
+ params_t new_params (params_B);
+ new_params.position = new_pos1;
+ new_row.push_back(drum_t(new_params));
+ index_1 = drums_.size()+new_row.size()-1;
+ //make adjacency vector longer
+ adjacency_vector_.push_back(std::vector<int>(3,-1));
+ }
+ else
+ index_1 += drums_.size();
+ int index_2 = find_close(new_pos2, new_row);
+ if(index_2 == -1){ //not in there yet, so push it
+ params_t new_params (params_B);
+ new_params.position = new_pos2;
+ new_row.push_back(drum_t(new_params));
+ index_2 = drums_.size()+new_row.size()-1;
+ //make adjacency vector longer
+ adjacency_vector_.push_back(std::vector<int>(3,-1));
+ }
+ else
+ index_2 += drums_.size();
+
+ //add new bonds to the adjacency vector
+ int old_index = std::distance(drums_.begin(), drum_it);
+ adjacency_vector_[old_index][1] = index_1;
+ adjacency_vector_[old_index][2] = index_2;
+ adjacency_vector_[index_1][1] = old_index;
+ adjacency_vector_[index_2][2] = old_index;
+
+ ++drum_it;
+ }
+
+ //push new drums
+ for(auto d: new_row)
+ drums_.push_back(d);
+
+ //clear new_row
+ new_row.clear();
+
+ //row 2
+ //prepare iterator
+ drum_it = drums_.end();
+ drum_it -= static_cast<int>(x_extent_ + 1); //this row has one drum more
+
+ while(drum_it != drums_.end()){
+ //These drums are guaranteed to be unique, hence inexistent.
+ params_t new_params (params_A);
+ new_params.position = drum_it->get_parameters().position - v_d;
+ new_row.push_back(drum_t(new_params));
+ int index = drums_.size() + new_row.size() - 1;
+ adjacency_vector_.push_back(std::vector<int>(3,-1));
+
+ int old_index = std::distance(drums_.begin(), drum_it);
+ adjacency_vector_[old_index][0] = index;
+ adjacency_vector_[index][0] = old_index;
+
+ ++drum_it;
+ }
+
+ //push new drums
+ for(auto d: new_row)
+ drums_.push_back(d);
+
+ //clear new_row
+ new_row.clear();
+
+ //row 3
+ drum_it = drums_.end();
+ drum_it -= static_cast<int>(x_extent_ + 1); //this row has one drum more
+
+ //first drum is special (only goes to ur)
+ params_t new_params (params_B);
+ new_params.position = drum_it->get_parameters().position + v_ur;
+ new_row.push_back(drum_t(new_params));
+ int indexx = drums_.size();
+ adjacency_vector_.push_back(std::vector<int>(3,-1));
+ int old_indexx = std::distance(drums_.begin(), drum_it); //avoid naming conflicts
+ adjacency_vector_[old_indexx][2] = indexx;
+ adjacency_vector_[indexx][2] = old_indexx;
+ ++drum_it;
+
+
+ //make an iterator that points to the last element
+ auto it_last_element = drums_.end();
+ --it_last_element;
+
+ //now the unspecial ones
+ while(drum_it != it_last_element){
+ Vec2<value_t> new_pos1 = drum_it->get_parameters().position + v_ul;
+ Vec2<value_t> new_pos2 = drum_it->get_parameters().position + v_ur;
+
+ //see if these drums already exist and add them to the list
+ int index_1 = find_close(new_pos1, new_row);
+ if(index_1 == -1){ //not in there yet, so push it
+ params_t new_params (params_B);
+ new_params.position = new_pos1;
+ new_row.push_back(drum_t(new_params));
+ index_1 = drums_.size()+new_row.size()-1;
+ //make adjacency vector longer
+ adjacency_vector_.push_back(std::vector<int>(3,-1));
+ }
+ else
+ index_1 += drums_.size();
+ int index_2 = find_close(new_pos2, new_row);
+ if(index_2 == -1){ //not in there yet, so push it
+ params_t new_params (params_B);
+ new_params.position = new_pos2;
+ new_row.push_back(drum_t(new_params));
+ index_2 = drums_.size()+new_row.size()-1;
+ //make adjacency vector longer
+ adjacency_vector_.push_back(std::vector<int>(3,-1));
+ }
+ else
+ index_2 += drums_.size();
+
+ //add new bonds to the adjacency vector
+ int old_index = std::distance(drums_.begin(), drum_it);
+ adjacency_vector_[old_index][1] = index_1;
+ adjacency_vector_[old_index][2] = index_2;
+ adjacency_vector_[index_1][1] = old_index;
+ adjacency_vector_[index_2][2] = old_index;
+
+ ++drum_it;
+ }
+
+ //the last drum is special too (only goes to ul, which already exists. only update adjacency.)
+ adjacency_vector_.back()[1] = drums_.size() - 1;
+ adjacency_vector_[drums_.size() - 1][1] = drums_.size() + new_row.size() - 1;
+
+ //push new drums
+ for(auto d: new_row)
+ drums_.push_back(d);
+
+ //clear new_row
+ new_row.clear();
+
+
+
+ //row 4
+ //prepare iterator
+ drum_it = drums_.end();
+ drum_it -= static_cast<int>(x_extent_);
+
+ while(drum_it != drums_.end()){
+ //These drums are guaranteed to be unique, hence inexistent.
+ params_t new_params (params_A);
+ new_params.position = drum_it->get_parameters().position - v_d;
+ new_row.push_back(drum_t(new_params));
+ int index = drums_.size() + new_row.size() - 1;
+ adjacency_vector_.push_back(std::vector<int>(3,-1));
+
+ int old_index = std::distance(drums_.begin(), drum_it);
+ adjacency_vector_[old_index][0] = index;
+ adjacency_vector_[index][0] = old_index;
+
+ ++drum_it;
+ }
+
+ //push new drums
+ for(auto d: new_row)
+ drums_.push_back(d);
+
+ //clear new_row
+ new_row.clear();
+ }
+
+ size_t x_extent_, y_extent_; //Number of equivalent hexagon rows in x and y direction
+ std::vector<Drum<value_t, params_t, vars_t, sbuffer_t> > drums_;
+ std::vector<std::vector<int> > adjacency_vector_;
+};
+
+#endif
diff --git a/projects/braidingTightBinding/include/vortex.hpp b/projects/braidingTightBinding/include/vortex.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..196e2c1c158c65fdde0488fa54f1adee01955c26
--- /dev/null
+++ b/projects/braidingTightBinding/include/vortex.hpp
@@ -0,0 +1,46 @@
+#ifndef VORTEX_HPP_INCLUDED
+#define VORTEX_HPP_INCLUDED
+#include <vec2.hpp>
+#include <cmath>
+
+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;
+
+ //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; }
+
+ //functional
+ //there is no divide_by_zero check for efficiency reasons. l0_ can not be zero upon call.
+ value_t operator()(const Vec2<value_t>& v, const Vec2<value_t>& s)
+ {
+ value_t prefactor = delta_*std::tanh(v.r_wrt(position_)/l0_);
+ value_t cosine = std::cos(alpha_ - v.phi_wrt(position_) + K_*(s + 2.*(v - position_)));
+ return prefactor*cosine;
+ }
+ private:
+ Vec2<value_t> position_;
+ const Vec2<value_t> K_;
+ value_t l0_, alpha_, delta_;
+};
+
+#endif
diff --git a/projects/braidingTightBinding/lib/coupler.hpp b/projects/braidingTightBinding/lib/coupler.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..ed35420a649e80033fd3202682ac4c36dc97156e
--- /dev/null
+++ b/projects/braidingTightBinding/lib/coupler.hpp
@@ -0,0 +1,17 @@
+#ifndef COUPLER_HPP_INCLUDED
+#define COUPLER_HPP_INCLUDED
+#include <vector>
+
+template <typename value_t, typename drum_t>
+class Coupler{
+ 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;
+};
+
+#endif
diff --git a/projects/braidingTightBinding/lib/diagonalizer.hpp b/projects/braidingTightBinding/lib/diagonalizer.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..5d5a8569b92628518d2d060ea40268715dcc0512
--- /dev/null
+++ b/projects/braidingTightBinding/lib/diagonalizer.hpp
@@ -0,0 +1,90 @@
+#ifndef DIAGONALIZER_HPP_INCLUDED
+#define DIAGONALIZER_HPP_INCLUDED
+#include <vector>
+#include <iostream>
+
+extern "C" void dsyev_(
+ char const & JOBZ, // 'N': Only eigenvalues, 'V': Eigenvalues and -vectors
+ char const & UPLO, // 'U': Upper triangular, 'V': Lower triangular
+ int const & N, // Matrix order
+ double * A, // Matrix and Eigenvector output
+ int const & LDA, // Matrix order (for my purposes)
+ double * W, // Eigenvalue output
+ double * WORK, // Workspace
+ int const & LWORK, // Size of workspace
+ int & INFO // Info
+ );
+
+//TODO: Rule of 7 conformity (nontrivial heap resource)
+class Diagonalizer{
+ public:
+ Diagonalizer() = default;
+ ~Diagonalizer() { delete[] work_; }
+
+ 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 if not done yet
+ if(N != N_ || dwork_ == 0.){
+ N_ = N;
+ prepare_workspace('N');
+ }
+
+ dsyev_('N', 'L', N_, matrix_.data(), N_, eigenvalues_.data(), work_, lwork_, info_);
+
+ if(info_ != 0)
+ throw("Diagonalization failed!");
+
+ return eigenvalues_;
+ }
+
+ std::pair<std::vector<double>, std::vector<double> > evv(const std::vector<double>& matrix, const size_t N){
+ //prepare memory
+ matrix_.clear();
+ matrix_.reserve(N*N - matrix_.capacity());
+ 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 if not done yet
+ if(dwork_ == 0.){
+ N_ = N;
+ prepare_workspace('V');
+ }
+
+ dsyev_('V', 'L', N_, matrix_.data(), N_, eigenvalues_.data(), work_, lwork_, info_);
+
+ if(info_ != 0)
+ throw("Diagonalization failed!");
+
+ 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_);
+ std::cout << "Allocating " << lwork_ << " doubles\n";
+ 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
diff --git a/projects/braidingTightBinding/lib/driver.hpp b/projects/braidingTightBinding/lib/driver.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..b1d4217d7f0257a3b189e187e2427c346690fa82
--- /dev/null
+++ b/projects/braidingTightBinding/lib/driver.hpp
@@ -0,0 +1,18 @@
+#ifndef DRIVER_HPP_INCLUDED
+#define DRIVER_HPP_INCLUDED
+#include <vector>
+#include <drum.hpp>
+
+template <typename value_t, typename drum_t>
+class Driver{
+ public:
+ //constructors
+ Driver() = default;
+ ~Driver() = 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 noexcept = 0;
+};
+
+#endif
diff --git a/projects/braidingTightBinding/lib/drum.hpp b/projects/braidingTightBinding/lib/drum.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..ee5f3f7252e624f2c9849bbfac659a284f458ca1
--- /dev/null
+++ b/projects/braidingTightBinding/lib/drum.hpp
@@ -0,0 +1,76 @@
+#ifndef DRUM_HPP_INCLUDED
+#define DRUM_HPP_INCLUDED
+#include <cmath>
+
+/* Neighbour enumeration:
+ * For drums in sublattice 'A', neighbour 0 is straight down, then clockwise.
+ * For drums in sublattice 'B', neighbour 0 is straight up, then clockwise.
+ * Hence neighbouring drums see each other as the same number neighbour.
+ */
+
+template <typename value_t, typename params_t, typename vars_t, typename sbuffer_t>
+class Drum{
+ public:
+ //'structors
+ Drum(const params_t& dp) noexcept: params_(dp) {}
+ Drum() = delete; //default construction not allowed
+ Drum(const Drum&) = default;
+ Drum& operator=(const Drum&) = default;
+ ~Drum() = default;
+
+ //access
+ params_t& get_parameters() noexcept
+ {
+ return params_;
+ }
+
+ const params_t& get_parameters() const noexcept
+ {
+ return params_;
+ }
+
+ vars_t& get_variables() noexcept
+ {
+ return vars_;
+ }
+
+ const vars_t& get_variables() const noexcept
+ {
+ return vars_;
+ }
+
+ sbuffer_t& get_sbuffer() noexcept
+ {
+ return sbuffer_;
+ }
+
+ const sbuffer_t& get_sbuffer() const noexcept
+ {
+ return sbuffer_;
+ }
+
+ //modifiers
+ void set_coupling_0(value_t t0)
+ {
+ vars_.t0 = t0;
+ }
+ void set_coupling_1(value_t t1)
+ {
+ vars_.t1 = t1;
+ }
+ void set_coupling_2(value_t t2)
+ {
+ vars_.t2 = t2;
+ }
+ void set_drive(value_t V)
+ {
+ vars_.V = V;
+ }
+
+ private:
+ params_t params_; // constant drum parameters
+ vars_t vars_; // drum variables
+ sbuffer_t sbuffer_; //stepper buffer
+};
+
+#endif
diff --git a/projects/braidingTightBinding/lib/force.hpp b/projects/braidingTightBinding/lib/force.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..094fffcfa4cbfb930090cd77860839ea80089668
--- /dev/null
+++ b/projects/braidingTightBinding/lib/force.hpp
@@ -0,0 +1,22 @@
+#ifndef FORCE_HPP_INCLUDED
+#define FORCE_HPP_INCLUDED
+#include <drum.hpp>
+
+//Force evaluation interface
+template <typename value_t, typename params_t, typename vars_t, typename buffer_t>
+class Force{
+ public:
+ //constructors
+ Force() = default;
+ ~Force() = default;
+
+ //virtual functional
+ virtual 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 = 0;
+};
+
+#endif
diff --git a/projects/braidingTightBinding/lib/grabber.hpp b/projects/braidingTightBinding/lib/grabber.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..911bd486690e09c742e831677b39ec696e65264a
--- /dev/null
+++ b/projects/braidingTightBinding/lib/grabber.hpp
@@ -0,0 +1,143 @@
+#ifndef GRABBER_HPP_INCLUDED
+#define GRABBER_HPP_INCLUDED
+#include <vector>
+#include <list>
+#include <fstream>
+#include <string>
+
+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;
+
+ 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;
+
+ //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;
+ }
+
+ //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 << "# a \t c \t f \t m \t Q \t gap \t coulomb_pref \t pos_x \t pos_y \t sublattice\n";
+ for(auto d: drums_){
+ par << d.get_parameters().a << " \t";
+ par << d.get_parameters().c << " \t";
+ par << d.get_parameters().f << " \t";
+ par << d.get_parameters().m << " \t";
+ par << d.get_parameters().Q << " \t";
+ par << d.get_parameters().gap << " \t";
+ par << d.get_parameters().coulomb_prefactor << " \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();
+
+ 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;
+ }
+
+ 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
diff --git a/projects/braidingTightBinding/lib/lattice_generator.hpp b/projects/braidingTightBinding/lib/lattice_generator.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..20819df1be678f9073606d752c15345dc8ec47be
--- /dev/null
+++ b/projects/braidingTightBinding/lib/lattice_generator.hpp
@@ -0,0 +1,20 @@
+#ifndef LATTICE_GENERATOR_HPP_INCLUDED
+#define LATTICE_GENERATOR_HPP_INCLUDED
+#include <utility>
+#include <vector>
+#include <drum.hpp>
+
+template <typename value_t, typename params_t, typename vars_t, typename sbuffer_t>
+class LatticeGenerator{
+ public:
+ LatticeGenerator() = default;
+ ~LatticeGenerator() = default;
+
+ //returns a std::pair of a std::vector of drums in the system and an adjacency vector of size_t.
+ //the last drum (.back()) is the drum seen when no neighbour is present.
+ //TODO: later, overload with a vector of drum parameters to set each drum's properties individually
+ virtual std::pair<std::vector<Drum<value_t, params_t, vars_t, sbuffer_t> >, std::vector<std::vector<int> > >
+ operator()(const params_t&) noexcept = 0;
+};
+
+#endif
diff --git a/projects/braidingTightBinding/lib/matrix_element_calculator.hpp b/projects/braidingTightBinding/lib/matrix_element_calculator.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..05aaa4bcafbeb4874275c853de2893fbf4aadbb9
--- /dev/null
+++ b/projects/braidingTightBinding/lib/matrix_element_calculator.hpp
@@ -0,0 +1,21 @@
+#ifndef MATRIX_ELEMENT_CALCULATOR_HPP_INCLUDED
+#define MATRIX_ELEMENT_CALCULATOR_HPP_INCLUDED
+
+template <typename value_t, typename params_t, typename vars_t, typename drum_t>
+class MatrixElementCalculator{
+ public:
+ MatrixElementCalculator() = default;
+ ~MatrixElementCalculator() = default;
+
+ //diagonal element at (index, index)
+ virtual value_t operator()(const size_t index, const std::vector<drum_t>& drums) const noexcept = 0;
+ //coupling elements at (index1, index2)
+ //for neighbour 0
+ virtual value_t operator()(const size_t index1, const size_t index2, const std::vector<drum_t>& drums) const noexcept = 0;
+ //for neighbour 1
+ virtual value_t operator()(const size_t index1, const size_t index2, const std::vector<drum_t>& drums, const int) const noexcept = 0;
+ //for neighbour 2
+ virtual value_t operator()(const size_t index1, const size_t index2, const std::vector<drum_t>& drums, const int, const int) const noexcept = 0;
+};
+
+#endif
diff --git a/projects/braidingTightBinding/lib/rk4_buffer.hpp b/projects/braidingTightBinding/lib/rk4_buffer.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..5e466227f195f00f9c191801332b5fa115a93b11
--- /dev/null
+++ b/projects/braidingTightBinding/lib/rk4_buffer.hpp
@@ -0,0 +1,16 @@
+#ifndef RK4_BUFFER_HPP_INCLUDED
+#define RK4_BUFFER_HPP_INCLUDED
+
+template <typename value_t>
+class RK4Buffer{
+ public:
+ RK4Buffer() = default;
+ ~RK4Buffer() = default;
+
+ public:
+ value_t k1, k2, k3, k4; //for velocity
+ value_t l1, l2, l3, l4; //for position
+
+};
+
+#endif
diff --git a/projects/braidingTightBinding/lib/rk4_stepper.hpp b/projects/braidingTightBinding/lib/rk4_stepper.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..1951a33ea67a7c0e880961fcd0ce56696c86a31d
--- /dev/null
+++ b/projects/braidingTightBinding/lib/rk4_stepper.hpp
@@ -0,0 +1,95 @@
+#ifndef RK4_STEPPER_HPP_INCLUDED
+#define RK4_STEPPER_HPP_INCLUDED
+
+template <typename value_t, typename params_t, typename vars_t, typename buffer_t,
+ typename force_t>
+class Rk4Stepper{
+ public:
+ Rk4Stepper() = default;
+ ~Rk4Stepper() = default;
+
+ //time signifies the state of the system, i.e. the time when the complete step was started.
+ //operations performed at t
+ void step_1(const force_t& force,
+ std::vector<Drum<value_t, params_t, vars_t, buffer_t> >& drums,
+ const std::vector<std::vector<int> >& adj_vec,
+ const value_t dt,
+ const value_t time) const noexcept
+ {
+ //step 1: calculate k1 and l1 for all drums
+ for(size_t i = 0; i < drums.size()-1; ++i){
+ //TODO: Empty neighbours are now signified by adj_vec[i][j] == -1. We could also push a 0-drum to the end and point to that.
+ drums[i].get_sbuffer().k1 = dt*(force(drums[i], drums[adj_vec[i][0]], drums[adj_vec[i][1]], drums[adj_vec[i][2]], time));
+ drums[i].get_sbuffer().l1 = dt*drums[i].get_variables().xdot;
+ }
+ //update x_temp, xdot_temp
+ for(size_t i = 0; i < drums.size()-1; ++i){
+ drums[i].get_variables().xdot_temp = drums[i].get_variables().xdot + drums[i].get_sbuffer().k1/2.;
+ drums[i].get_variables().x_temp = drums[i].get_variables().x + drums[i].get_sbuffer().l1/2.;
+ }
+ //now it's time for the owner to update the coupler and driver to t+dt/2.
+
+
+ }
+
+
+ //operations performed at t+dt/2
+ void step_2(const force_t& force,
+ std::vector<Drum<value_t, params_t, vars_t, buffer_t> >& drums,
+ const std::vector<std::vector<int> >& adj_vec,
+ const value_t dt,
+ const value_t time) const noexcept
+ {
+ //step 2: calculate k2 and l2 for all drums
+ for(size_t i = 0; i < drums.size()-1; ++i){
+ //TODO: Empty neighbours are now signified by adj_vec[i][j] == -1. We could also push a 0-drum to the end and point to that.
+ drums[i].get_sbuffer().k2 = dt*(force(drums[i], drums[adj_vec[i][0]], drums[adj_vec[i][1]], drums[adj_vec[i][2]], time));
+ drums[i].get_sbuffer().l2 = dt*drums[i].get_variables().xdot_temp;
+ }
+ //update x_temp, xdot_temp
+ for(size_t i = 0; i < drums.size()-1; ++i){
+ drums[i].get_variables().xdot_temp = drums[i].get_variables().xdot + drums[i].get_sbuffer().k2/2.;
+ drums[i].get_variables().x_temp = drums[i].get_variables().x + drums[i].get_sbuffer().l2/2.;
+ }
+
+
+
+ //step 3: calculate k3 and l3 for all drums
+ for(size_t i = 0; i < drums.size()-1; ++i){
+ //TODO: Empty neighbours are now signified by adj_vec[i][j] == -1. We could also push a 0-drum to the end and point to that.
+ drums[i].get_sbuffer().k3 = dt*(force(drums[i], drums[adj_vec[i][0]], drums[adj_vec[i][1]], drums[adj_vec[i][2]], time));
+ drums[i].get_sbuffer().l3 = dt*drums[i].get_variables().xdot_temp;
+ }
+ //update x_temp, xdot_temp
+ for(size_t i = 0; i < drums.size()-1; ++i){
+ drums[i].get_variables().xdot_temp = drums[i].get_variables().xdot + drums[i].get_sbuffer().k3;
+ drums[i].get_variables().x_temp = drums[i].get_variables().x + drums[i].get_sbuffer().l3;
+ }
+ //now it's time for the owner to update the coupler and driver to t+dt/2.
+ }
+
+
+ //operations performed at t+dt
+ void step_3(const force_t& force,
+ std::vector<Drum<value_t, params_t, vars_t, buffer_t> >& drums,
+ const std::vector<std::vector<int> >& adj_vec,
+ const value_t dt,
+ const value_t time) const noexcept
+ {
+ //step 4: calculate k4 and l4 for all drums
+ for(size_t i = 0; i < drums.size()-1; ++i){
+ //TODO: Empty neighbours are now signified by adj_vec[i][j] == -1. We could also push a 0-drum to the end and point to that.
+ drums[i].get_sbuffer().k4 = dt*(force(drums[i], drums[adj_vec[i][0]], drums[adj_vec[i][1]], drums[adj_vec[i][2]], time));
+ drums[i].get_sbuffer().l4 = dt*drums[i].get_variables().xdot_temp;
+ }
+ //update x, xdot, x_temp, xdot_temp to t+dt
+ for(size_t i = 0; i < drums.size()-1; ++i){
+ drums[i].get_variables().xdot_temp = drums[i].get_variables().xdot + drums[i].get_sbuffer().k4;
+ drums[i].get_variables().x_temp = drums[i].get_variables().x + drums[i].get_sbuffer().l4;
+ drums[i].get_variables().xdot = drums[i].get_variables().xdot + drums[i].get_sbuffer().k4;
+ drums[i].get_variables().x = drums[i].get_variables().x + drums[i].get_sbuffer().l4;
+ }
+ }
+};
+
+#endif
diff --git a/projects/braidingTightBinding/lib/system.hpp b/projects/braidingTightBinding/lib/system.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..37855c12fa61b561f5846b77ecca5f3b1a7f452b
--- /dev/null
+++ b/projects/braidingTightBinding/lib/system.hpp
@@ -0,0 +1,115 @@
+#ifndef SYSTEM_HPP_INCLUDED
+#define SYSTEM_HPP_INCLUDED
+#include <vector>
+#include <iostream>
+
+template<typename value_t, typename drum_t, typename grabber_t, typename sysparams_t, typename force_t, typename coupler_t, typename driver_t, typename stepper_t, typename matelecalc_t>
+class System{
+ public:
+ System(const value_t t_end, const value_t dt, const std::vector<drum_t>& drums, const stepper_t stepper, const force_t force, const sysparams_t sysparams, const grabber_t grabber)
+ : drums_(drums), stepper_(stepper), force_(force), sysparams_(sysparams), grabber_(grabber), t_end_(t_end), dt_(dt), time_(0.)
+ {
+ sysparams_.coupler.precompute(t_end, dt, drums);
+ sysparams_.driver.precompute(t_end, dt, drums);
+ grabber_.init(t_end, dt, drums, sysparams.adjacency_vector);
+
+ push_dc(); //push the initial values
+
+ }
+ ~System() = default;
+
+ void simulate(){
+ while(time_ <= t_end_){
+ grabber_.grab(drums_, time_);
+ step();
+ }
+ }
+
+ void step(){
+ push_dc();
+ stepper_.step_1(force_, drums_, sysparams_.adjacency_vector, dt_, time_);
+ step_dc(dt_/2.);
+ push_dc();
+ stepper_.step_2(force_, drums_, sysparams_.adjacency_vector, dt_, time_);
+ step_dc(dt_/2.);
+ push_dc();
+ stepper_.step_3(force_, drums_, sysparams_.adjacency_vector, dt_, time_);
+
+ time_ += dt_;
+ }
+
+ void reset_time(){
+ time_ = value_t(0.);
+ }
+
+ void set_step(value_t dt){
+ dt_ = dt;
+ }
+
+ bool save(){
+ return grabber_.save();
+ }
+
+ std::vector<value_t> get_matrix(const matelecalc_t& mec){
+ calculate_matrix(mec);
+ return matrix_;
+ }
+
+ private:
+ void step_dc(const value_t dt){
+ sysparams_.coupler.step(dt);
+ sysparams_.driver.step(dt);
+ }
+
+ void push_dc(){
+ for(size_t i = 0; i < drums_.size()-1; ++i){
+ drums_[i].get_variables().V = sysparams_.driver(i);
+ drums_[i].get_variables().t0 = sysparams_.coupler(i,0);
+ drums_[i].get_variables().t1 = sysparams_.coupler(i,1);
+ drums_[i].get_variables().t2 = sysparams_.coupler(i,2);
+ }
+ }
+
+ //Calculates dynamical matrix for the currently pushed values
+ void calculate_matrix(const matelecalc_t& mec){
+ //for convenience
+ size_t N = drums_.size()-1;
+ //clear the matrix
+ matrix_.clear();
+ //prepare memory
+ matrix_.reserve((drums_.size()-1)*(drums_.size()-1) - matrix_.capacity());
+ //initialize
+ for(size_t i = 0; i < (drums_.size()-1) * (drums_.size()-1); ++i)
+ matrix_.push_back(value_t(0.));
+
+ for(size_t i = 0; i < drums_.size()-1; ++i){
+ size_t n0_index = sysparams_.adjacency_vector[i][0];
+ size_t n1_index = sysparams_.adjacency_vector[i][1];
+ size_t n2_index = sysparams_.adjacency_vector[i][2];
+ //diagonal
+ matrix_[N*i + i] = mec(i, drums_);
+ //couplings
+ //neighbour 0
+ if(n0_index != drums_.size()-1) [[likely]] //actually a neighbour
+ matrix_[N*i + n0_index] = mec(i, n0_index, drums_);
+ //neighbour 1
+ if(n1_index != drums_.size()-1) [[likely]] //actually a neighbour
+ matrix_[N*i + n1_index] = mec(i, n1_index, drums_, 1);
+ //neighbour 2
+ if(n2_index != drums_.size()-1) [[likely]] //actually a neighbour
+ matrix_[N*i + n2_index] = mec(i, n2_index, drums_, 2, 2);
+ }
+ }
+
+ std::vector<drum_t> drums_; //The last drum is NOT TO BE TOUCHED!
+ sysparams_t sysparams_;
+ grabber_t grabber_;
+ stepper_t stepper_;
+ force_t force_;
+ value_t dt_, t_end_;
+ value_t time_;
+
+ std::vector<value_t> matrix_; //Matrix of this problem
+};
+
+#endif
diff --git a/projects/braidingTightBinding/lib/system_parameters.hpp b/projects/braidingTightBinding/lib/system_parameters.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..637d7792be595c26bc2c0eeb8a89652b5dbc19be
--- /dev/null
+++ b/projects/braidingTightBinding/lib/system_parameters.hpp
@@ -0,0 +1,18 @@
+#ifndef SYSTEM_PARAMETERS_HPP_INCLUDED
+#define SYSTEM_PARAMETERS_HPP_INCLUDED
+#include <vector>
+
+template <typename coupler_t, typename driver_t>
+class SystemParameters{
+ public:
+ //constructors
+ SystemParameters(coupler_t coupler, driver_t driver, std::vector< std::vector<int> > adjacency_vector) noexcept: coupler(coupler), driver(driver), adjacency_vector(adjacency_vector) {}
+
+
+ public:
+ coupler_t coupler;
+ driver_t driver;
+ std::vector< std::vector<int> > adjacency_vector;
+};
+
+#endif
diff --git a/projects/braidingTightBinding/lib/vec2.hpp b/projects/braidingTightBinding/lib/vec2.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..a2bb4c3aabc8ba0cace7bc557ab5d14c4afd1287
--- /dev/null
+++ b/projects/braidingTightBinding/lib/vec2.hpp
@@ -0,0 +1,150 @@
+#ifndef VEC2_HPP_INCLUDED
+#define VEC2_HPP_INCLUDED
+#include <cmath>
+#include <iostream>
+
+template <typename value_t>
+class Vec2{
+ public:
+ //constructors
+ Vec2(const value_t x, const value_t y) noexcept: x_(x), y_(y) {}
+ Vec2(const Vec2<value_t>& other) = default;
+ Vec2& operator=(const Vec2&) = default;
+ Vec2() = default;
+ ~Vec2() = default;
+
+
+ //access
+ value_t x() const noexcept { return x_; }
+ value_t y() const noexcept { return y_; }
+ Vec2 normalized() const { return *this/this->norm(); }
+
+ //TODO: Precalculate these and flag when recomputation is needed
+ value_t r() const noexcept { return this->norm(); }
+ value_t phi() const
+ {
+ value_t angle (0);
+ angle = std::atan2(y_, x_);
+ if(std::isnan(angle)){
+ throw("ATAN2 NAN OCCURRED");
+ return 0.;
+ }
+ return angle;
+ }
+
+ //function members
+ value_t r_wrt(const Vec2& origin) const noexcept { return (*this - origin).r(); }
+ value_t phi_wrt(const Vec2& origin) const { return (*this - origin).phi(); }
+
+ value_t norm() const noexcept { return std::sqrt(x_*x_ + y_*y_); }
+ value_t norm_sq() const noexcept { return x_*x_ + y_* y_; }
+
+ //modifiers
+ Vec2& rotate(const Vec2& center, const value_t degrees)
+ {
+ *this -= center;
+ value_t temp = x_;
+ const value_t radians = degrees*3.141592653589793/180.;
+
+ x_ = std::cos(radians)*x_ - std::sin(radians)*y_;
+ y_ = std::sin(radians)*temp + std::cos(radians)*y_;
+
+ *this += center;
+
+ return *this;
+ }
+
+ Vec2& normalize()
+ {
+ value_t length = norm();
+ if(length == value_t(0))
+ throw("NORMALIZE ZERO VEC2 ATTEMPTED");
+ x_ /= length;
+ y_ /= length;
+ return *this;
+ }
+
+ //overloading
+ Vec2& operator+=(const Vec2& rhs) noexcept
+ {
+ x_ += rhs.x_;
+ y_ += rhs.y_;
+ return *this;
+ }
+ Vec2& operator-=(const Vec2& rhs) noexcept
+ {
+ x_ -= rhs.x_;
+ y_ -= rhs.y_;
+ return *this;
+ }
+ Vec2& operator*=(const value_t s) noexcept
+ {
+ x_ *= s;
+ y_ *= s;
+ return *this;
+ }
+ Vec2& operator/=(const value_t s)
+ {
+ if(x_ == value_t(0) and y_ == value_t(0))
+ return *this;
+ else if(s == value_t(0))
+ throw("DIVISION BY ZERO ATTEMPTED");
+ x_ /= s;
+ y_ /= s;
+ return *this;
+ }
+
+ value_t& operator[](std::size_t i) noexcept
+ {
+ if(i == 0)
+ return x_;
+ else
+ return y_;
+ }
+ value_t operator[](std::size_t i) const noexcept
+ {
+ if(i == 0)
+ return x_;
+ else
+ return y_;
+ }
+
+ //non-ADL hidden friends
+ friend Vec2 operator+(const Vec2& lhs, Vec2 rhs) noexcept
+ {
+ rhs += lhs;
+ return rhs;
+ }
+ friend Vec2 operator-(Vec2 lhs, const Vec2& rhs) noexcept
+ {
+ lhs -= rhs;
+ return lhs;
+ }
+ friend Vec2 operator*(const value_t s, Vec2 v) noexcept
+ {
+ v *= s;
+ return v;
+ }
+ friend Vec2 operator*(Vec2 v, const value_t s) noexcept
+ {
+ v *= s;
+ return v;
+ }
+ friend Vec2 operator/(Vec2 v, const value_t s)
+ {
+ v /= s;
+ return v;
+ }
+ friend value_t operator*(const Vec2& v, const Vec2& w) noexcept { return v.x_*w.x_ + v.y_*w.y_; }//dot product
+
+ private:
+ value_t x_, y_;
+};
+
+template <typename value_t>
+std::ostream& operator<<(std::ostream& os, const Vec2<value_t>& v)
+{
+ return os << "(" << v.x() << ", " << v.y() << ")";
+}
+
+#endif