feat: generalize hopping matrix efficiently!

Duumbo · Jun 13, 2024 · 33e9efd · 33e9efd
1 parent 11565a3
commit 33e9efd
Show file tree

Hide file tree

Showing 3 changed files with 116 additions and 40 deletions.
diff --git a/src/bin/dvmc.rs b/src/bin/dvmc.rs
@@ -2,13 +2,13 @@ use rand::Rng;
 use std::ptr::addr_of;
 
 use impurity::{FockState, RandomStateGeneration, VarParams};
-use impurity::{FIJ, GI, VIJ};
+use impurity::{FIJ, GI, VIJ, SIZE};
 use impurity::density::compute_internal_product;
 use impurity::hamiltonian::{potential, kinetic};
+use impurity::HOP_BITMASKS;
 
-const NELEC: usize = 4;
-const SIZE: usize = 4;
-const NMCSAMP: usize = 100;
+const NELEC: usize = 8;
+const NMCSAMP: usize = 1000;
 
 
 fn propose_hopping<R: Rng + ?Sized>(state: &FockState<u8>, rng: &mut R, params: &VarParams) -> (f64, FockState<u8>) {
@@ -34,8 +34,8 @@ fn propose_hopping<R: Rng + ?Sized>(state: &FockState<u8>, rng: &mut R, params:
 //}
 
 fn compute_hamiltonian(state: FockState<u8>, ip: f64, params: &VarParams) -> f64 {
-    let kin = <f64>::ln(kinetic(state, params));
-    <f64>::exp(kin - ip) + potential(state)
+    let kin = kinetic(state, params);
+    (kin / <f64>::exp(ip)) + potential(state)
 }
 
 fn main() {
@@ -71,4 +71,5 @@ fn main() {
     }
     energy = energy / NMCSAMP as f64;
     println!("Energy: {}", energy);
+    println!("Hop bitmasks: {:08b}, {:08b}", HOP_BITMASKS[0], HOP_BITMASKS[1]);
 }
diff --git a/src/hamiltonian.rs b/src/hamiltonian.rs
@@ -1,4 +1,4 @@
-use crate::{density::compute_internal_product, BitOps, FockState, CONS_U, Hopper};
+use crate::{density::compute_internal_product, BitOps, FockState, CONS_U, Hopper, HOPPINGS, SIZE};
 use crate::{VarParams, CONS_T};
 
 /// Computes the potential term of the Hamiltonian.
@@ -49,7 +49,7 @@ where
             f_state.spin_up.set(hop.0);
             f_state.spin_up.set(hop.1);
         }
-        kin += compute_internal_product(f_state, params)*CONS_T;
+        kin += compute_internal_product(f_state, params)*CONS_T*HOPPINGS[hop.0 + hop.1*SIZE];
     }
 
     kin

diff --git a/src/params.rs b/src/params.rs
@@ -1,7 +1,8 @@
-use macros::generate_bitmasks_from_hoppings;
 /// Size of the system.
-pub const SIZE: usize = 4;
+pub const SIZE: usize = 8;
 pub const ARRAY_SIZE: usize = (SIZE + 7) / 8;
+pub type BitStruct = u8;
+pub const NBITS: usize = 8;
 
 /// Hubbard's model $U$ parameter
 pub static CONS_U: f64 = 0.0;
@@ -20,44 +21,118 @@ pub const NVIJ: usize = SIZE*SIZE;
 pub const NGI: usize = SIZE;
 
 pub static mut FIJ: [f64; NFIJ] = [
-    0.65583263, 0.2579242 , 0.78342288, 0.86818364, 0.99423766,
-    0.38035571, 0.47056418, 0.56113084, 0.1475486 , 0.5369956 ,
-    0.26664678, 0.42718298, 0.23611604, 0.43433408, 0.53575446,
-    0.07166779, 0.15895308, 0.71804482, 0.2714126 , 0.71162796,
-    0.77433389, 0.28221235, 0.14508085, 0.49846785, 0.66801568,
-    0.32978509, 0.6571238 , 0.03039801, 0.3468509 , 0.15403364,
-    0.37306052, 0.90458229, 0.91791243, 0.43872965, 0.56242938,
-    0.23684173, 0.46798428, 0.8399697 , 0.20428968, 0.46700014,
-    0.08779668, 0.4836716 , 0.82304382, 0.46428091, 0.35912961,
-    0.52023224, 0.12116029, 0.2411605 , 0.21063809, 0.19513569,
-    0.5357807 , 0.59983385, 0.11489229, 0.79493793, 0.63347575,
-    0.45736727, 0.85445574, 0.00426934, 0.41213109, 0.41722305,
-    0.91049378, 0.26227231, 0.67512998, 0.89924604
+    0.73454851, 0.79079811, 0.02890209, 0.96743271, 0.55712194,
+       0.53366545, 0.90578359, 0.2536    , 0.14239049, 0.39043369,
+       0.27090031, 0.95716023, 0.89890399, 0.90824588, 0.75315743,
+       0.90976429, 0.9039633 , 0.00135397, 0.31248781, 0.87812967,
+       0.31596079, 0.04218799, 0.07100452, 0.46693654, 0.8694935 ,
+       0.2659505 , 0.32777126, 0.21211653, 0.10695527, 0.89737699,
+       0.77005737, 0.65921572, 0.10743515, 0.42278062, 0.176886  ,
+       0.3833599 , 0.84193656, 0.57276593, 0.60844303, 0.27030017,
+       0.90519023, 0.90095001, 0.29416643, 0.30291901, 0.43374941,
+       0.75226538, 0.15869649, 0.14401949, 0.34860453, 0.94639767,
+       0.21215991, 0.42571408, 0.96594105, 0.32025106, 0.11556635,
+       0.7382251 , 0.01082574, 0.3741212 , 0.79539587, 0.3839405 ,
+       0.31366624, 0.99487885, 0.12603209, 0.63806619, 0.54221907,
+       0.30298915, 0.29624267, 0.93190217, 0.07538637, 0.12080278,
+       0.45480728, 0.39749362, 0.53047916, 0.54561027, 0.97786655,
+       0.1031861 , 0.31242597, 0.29231074, 0.89732391, 0.81015248,
+       0.16532583, 0.79491769, 0.20926697, 0.34544669, 0.79878245,
+       0.37496588, 0.35622598, 0.79120737, 0.51771781, 0.98350729,
+       0.92669342, 0.53872625, 0.01591008, 0.28435145, 0.43073094,
+       0.84013006, 0.69663584, 0.67117945, 0.21636238, 0.43222978,
+       0.91599643, 0.73676706, 0.90278944, 0.80216004, 0.31563004,
+       0.81211759, 0.80539688, 0.6039476 , 0.49381156, 0.55864535,
+       0.94775156, 0.71903457, 0.73560788, 0.37442695, 0.16568087,
+       0.06257877, 0.17409916, 0.49573729, 0.36348517, 0.49636986,
+       0.16317322, 0.5147186 , 0.30106074, 0.1724178 , 0.08779814,
+       0.28648642, 0.05373693, 0.10446211, 0.12991774, 0.99196729,
+       0.3663934 , 0.56518385, 0.96664491, 0.79773777, 0.77428983,
+       0.31597117, 0.11301903, 0.59287381, 0.60232572, 0.36598577,
+       0.39785319, 0.02049467, 0.74161357, 0.7027563 , 0.31616185,
+       0.0948578 , 0.67187836, 0.59165484, 0.33743617, 0.31508259,
+       0.68012331, 0.6193884 , 0.50465435, 0.02363093, 0.99257409,
+       0.21343555, 0.05923286, 0.06809502, 0.89038897, 0.28383551,
+       0.66777553, 0.36549757, 0.31072894, 0.23127301, 0.3547174 ,
+       0.97725038, 0.34559845, 0.35087415, 0.90526068, 0.25512045,
+       0.31396858, 0.22995697, 0.94415834, 0.52877933, 0.37381948,
+       0.9492286 , 0.81304415, 0.56846197, 0.42051353, 0.36764267,
+       0.20839462, 0.9916981 , 0.66552448, 0.93551247, 0.1033957 ,
+       0.56389839, 0.83568559, 0.93820606, 0.51842815, 0.15110758,
+       0.50414457, 0.44204805, 0.15785273, 0.19707327, 0.91630047,
+       0.89832097, 0.083598  , 0.50487228, 0.03308039, 0.29779481,
+       0.05836903, 0.63053828, 0.8276419 , 0.77818981, 0.54681855,
+       0.46639428, 0.47894859, 0.05500165, 0.07338289, 0.69701156,
+       0.29120426, 0.423073  , 0.63817691, 0.19338117, 0.63491751,
+       0.78108888, 0.59436187, 0.83914825, 0.48350799, 0.04073812,
+       0.83191958, 0.36907401, 0.92450273, 0.35130647, 0.54268754,
+       0.37796119, 0.34816465, 0.8030052 , 0.15185223, 0.09852173,
+       0.91883544, 0.94385193, 0.12352969, 0.83164358, 0.5486104 ,
+       0.67315083, 0.78069208, 0.77094661, 0.7384394 , 0.75453324,
+       0.84711911, 0.53715719, 0.51538202, 0.28472542, 0.89775889,
+       0.71417763, 0.62104784, 0.75099927, 0.90738298, 0.43109713,
+       0.79622524, 0.50342189, 0.05948419, 0.66338672, 0.39778654,
+       0.69899412
 ];
 
 pub static mut VIJ: [f64; NVIJ] = [
-    0.38200108, 0.39699789, 0.99773244, 0.06024556, 0.89587208,
-    0.72621538, 0.75943934, 0.47029373, 0.99266814, 0.09317458,
-    0.4747971 , 0.00407951, 0.61192636, 0.98595942, 0.24153617,
-    0.07105668
+    0.16441764, 0.57378353, 0.08483113, 0.85996846, 0.44735341,
+       0.28084081, 0.12064007, 0.06950027, 0.72987928, 0.73925121,
+       0.02788582, 0.22349883, 0.80942185, 0.22165408, 0.45886613,
+       0.29862472, 0.15722897, 0.6747521 , 0.95284056, 0.45173082,
+       0.9726399 , 0.88055554, 0.94868543, 0.56493827, 0.77562641,
+       0.21161496, 0.4228778 , 0.402329  , 0.61999778, 0.47456151,
+       0.069231  , 0.64847089, 0.91819044, 0.35393685, 0.60485968,
+       0.73834106, 0.97126728, 0.21180482, 0.66585592, 0.07816157,
+       0.45303348, 0.35036158, 0.19635927, 0.00101243, 0.95448494,
+       0.58267052, 0.09669098, 0.89369989, 0.44613375, 0.32236403,
+       0.8426723 , 0.98566725, 0.9600568 , 0.56439038, 0.10817061,
+       0.07101563, 0.10087598, 0.76554334, 0.81188881, 0.26653177,
+       0.70735203, 0.79009915, 0.61491369, 0.39787751
 ];
 
 pub static mut GI: [f64; NGI] = [
-    0.41160929, 0.37086116, 0.828815, 0.78192488
+    0.91300853, 0.85431752, 0.40335153, 0.80450021, 0.61176878,
+       0.34797097, 0.52353966, 0.89905992
 ];
 
 pub const HOPPINGS: [f64; SIZE*SIZE] = [
-    0.0, 1.0, 0.0, 1.0,
-    1.0, 0.0, 1.0, 0.0,
-    0.0, 1.0, 0.0, 1.0,
-    1.0, 0.0, 1.0, 0.0
+    0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0,
+    1.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0,
+    0.0, 1.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0,
+    0.0, 0.0, 1.0, 0.0, 1.0, 0.0, 0.0, 0.0,
+    0.0, 0.0, 0.0, 1.0, 0.0, 1.0, 0.0, 0.0,
+    0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 1.0, 0.0,
+    0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 1.0,
+    1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0,
 ];
 
-generate_bitmasks_from_hoppings!(
-[
-    0.0, 1.0, 0.0, 1.0,
-    1.0, 0.0, 1.0, 0.0,
-    0.0, 1.0, 0.0, 1.0,
-    1.0, 0.0, 1.0, 0.0
-]
-);
+pub const HOP_BITMASKS: [BitStruct; SIZE / 2] = {
+    let mut hop_tmp: [BitStruct; SIZE / 2] = [BitStruct::from_be(0); SIZE / 2];
+    // Index for array
+    let mut i: usize = 0;
+    while i < SIZE / 2 {
+        hop_tmp[i] = {
+            let mut mask = BitStruct::from_be(0);
+            let one = BitStruct::from_be(1);
+            let mut j: usize = 0;
+            while j < (SIZE - 1 - i) {
+                if HOPPINGS[j + i + 1 + SIZE * j] != 0.0 {
+                    if i == 0 {
+                        mask ^= one << (NBITS - j - 2);
+                    } else {
+                        mask ^= one << (NBITS - j - 1);
+                    }
+                }
+                j += 1;
+            }
+            if i == 0 && HOPPINGS[SIZE - 1] != 0.0 {
+                mask ^= one << (NBITS - 1);
+            }
+            // Index for  HOPPINGS
+            mask
+        };
+        i += 1;
+    }
+    hop_tmp
+};