feat: Compute internal product and random state generation

data/orbitale.csv

@@ -27,3 +27,11 @@ i,j,n
 6,5,26
 7,5,27
 7,6,28
+0,0,29
+1,1,30
+2,2,31
+3,3,32
+4,4,33
+5,5,34
+6,6,35
+7,7,36

data/orbitale.yaml

@@ -0,0 +1,30 @@
+---
+fij:
+  - [1,0,1]
+  - [2,0,2]
+  - [3,0,3]
+  - [4,0,4]
+  - [5,0,5]
+  - [6,0,6]
+  - [7,0,7]
+  - [2,1,8]
+  - [3,1,9]
+  - [4,1,10]
+  - [5,1,11]
+  - [6,1,12]
+  - [7,1,13]
+  - [3,2,14]
+  - [4,2,15]
+  - [5,2,16]
+  - [6,2,17]
+  - [7,2,18]
+  - [4,3,19]
+  - [5,3,20]
+  - [6,3,21]
+  - [7,3,22]
+  - [5,4,23]
+  - [6,4,24]
+  - [7,4,25]
+  - [6,5,26]
+  - [7,5,27]
+  - [7,6,28]

src/fock_state.rs

@@ -362,6 +362,7 @@ impl<T: BitOps> Distribution<FockState<T>> for Standard where Standard: Distribu
 
 pub trait RandomStateGeneration {
     fn generate_from_nelec<R: Rng + ?Sized>(rng: &mut R, nelec: usize, max_size: usize) -> Self;
+    fn generate_hopping<R: Rng + ?Sized>(self: &Self, rng: &mut R, max_size: u32) -> Self;
 }
 
 impl<T: BitOps> RandomStateGeneration for FockState<T> where Standard: Distribution<T> {
@@ -386,6 +387,50 @@ impl<T: BitOps> RandomStateGeneration for FockState<T> where Standard: Distribut
         }
         state
     }
+
+    fn generate_hopping<R: Rng + ?Sized>(self: &FockState<T>, rng: &mut R, max_size: u32) -> FockState<T> {
+        // Test for empty state
+        if (self.spin_up.count_ones() == 0) && (self.spin_down.count_ones() == 0) {
+            return self.clone();
+        }
+        // Test for full state
+        if (self.spin_up.count_ones() == max_size) && (self.spin_down.count_ones() == max_size) {
+            return self.clone();
+        }
+        // Choose up or down
+        let mut spin = rng.gen::<bool>();
+        if spin && (self.spin_up.count_ones() == max_size) { spin = ! spin;}
+        if !spin && (self.spin_down.count_ones() == max_size) { spin = ! spin;}
+        let mut random_to = <u32>::MAX;
+        let random_from;
+        let mut sup = self.spin_up.clone();
+        let mut sdown = self.spin_down.clone();
+        if spin {
+            random_from = rng.gen::<u32>() % self.spin_up.count_ones();
+            while random_to == <u32>::MAX {
+                let index = rng.gen::<u32>() % max_size;
+                if !self.spin_up.check(index as usize) {random_to = index;}
+            }
+            sup.set(random_to as usize);
+            sup.set(random_from as usize);
+        }
+        else {
+            random_from = rng.gen::<u32>() % self.spin_down.count_ones();
+            while random_to == <u32>::MAX {
+                let index = rng.gen::<u32>() % max_size;
+                if !self.spin_down.check(index as usize) {random_to = index;}
+            }
+            sdown.set(random_to as usize);
+            sdown.set(random_from as usize);
+        }
+
+        FockState{
+            n_sites: max_size as usize,
+            spin_up: sup,
+            spin_down: sdown,
+        }
+
+    }
 }
 
 

src/hamiltonian.rs

@@ -1,4 +1,4 @@
-use crate::{BitOps, FockState, CONS_U, SIZE};
+use crate::{BitOps, FockState, CONS_U};
 
 /// Computes the potential term of the Hamiltonian.
 /// # Arguments
@@ -31,7 +31,7 @@ where
 /// $$
 /// H_T=-t\sum_{<i,j>,\sigma}c^\dagger_{i\sigma}c_{j\sigma}+c^\dagger_{j\sigma}c_{i\sigma}
 /// $$
-pub fn kinetic<T>(spin_up: T, spin_down: T) -> Vec<FockState<T>>
+pub fn kinetic<T>(spin_up: T, spin_down: T, size: usize) -> Vec<FockState<T>>
 where
     T: BitOps + From<u8>,
 {
@@ -50,65 +50,65 @@ where
     let can_gol_spin_down = spin_down_shr1 ^ spin_down;
 
     out.append(
-        &mut tm(spin_up, can_gor_spin_up, 2)
+        &mut tm(spin_up, can_gor_spin_up, 2, size)
             .into_iter()
             .map(|s| FockState {
                 spin_up: s,
                 spin_down,
-                n_sites: SIZE,
+                n_sites: size,
             })
             .collect::<Vec<FockState<T>>>(),
     );
     out.append(
-        &mut tm(spin_down, can_gor_spin_down, 2)
+        &mut tm(spin_down, can_gor_spin_down, 2, size)
             .into_iter()
             .map(|s| FockState {
                 spin_up,
                 spin_down: s,
-                n_sites: SIZE,
+                n_sites: size,
             })
             .collect::<Vec<FockState<T>>>(),
     );
     out.append(
-        &mut tm(spin_up, can_gol_spin_up, 1)
+        &mut tm(spin_up, can_gol_spin_up, 1, size)
             .into_iter()
             .map(|s| FockState {
                 spin_up: s,
                 spin_down,
-                n_sites: SIZE,
+                n_sites: size,
             })
             .collect::<Vec<FockState<T>>>(),
     );
     out.append(
-        &mut tm(spin_down, can_gol_spin_down, 1)
+        &mut tm(spin_down, can_gol_spin_down, 1, size)
             .into_iter()
             .map(|s| FockState {
                 spin_up,
                 spin_down: s,
-                n_sites: SIZE,
+                n_sites: size,
             })
             .collect::<Vec<FockState<T>>>(),
     );
 
     out
 }
 
-fn tm<T>(spin: T, mut truth: T, shl_qt: usize) -> Vec<T>
+fn tm<T>(spin: T, mut truth: T, shl_qt: usize, size: usize) -> Vec<T>
 where
     T: BitOps + From<u8>,
 {
     let mut i = truth.leading_zeros();
     let mut out_vec: Vec<T> = Vec::with_capacity(8);
-    while (i as usize) < SIZE {
-        let n = SIZE as i32 - shl_qt as i32 - i as i32;
+    while (i as usize) < size {
+        let n = size as i32 - shl_qt as i32 - i as i32;
         let mask;
         if n < 0 {
             mask = T::from(3).rotate_right(n.abs() as u32);
         } else {
             mask = T::from(3).rotate_left(n as u32);
         }
         out_vec.push(spin ^ mask);
-        truth ^= T::from(1 << (SIZE - 1 - i as usize));
+        truth ^= T::from(1 << (size - 1 - i as usize));
         i = truth.leading_zeros();
     }
     out_vec

src/jastrow.rs

@@ -409,6 +409,7 @@ mod test {
 
             // Copy the state
             let fock_state2 = fock_state1.clone();
+
             let mut jastrow_params: Vec<f64> = Vec::with_capacity(SIZE * SIZE);
             for _ in 0..SIZE * SIZE {
                 jastrow_params.push(rng.gen::<f64>());

src/main.rs

@@ -6,10 +6,10 @@ use impurity::hamiltonian::{kinetic, potential};
 use impurity::jastrow::compute_jastrow_exp;
 use impurity::parse::orbitale::parse_orbitale_def;
 use impurity::pfaffian::construct_matrix_a_from_state;
-use impurity::{FockState, CONS_T, SIZE};
+use impurity::{FockState, CONS_T, RandomStateGeneration};
 
-const UP: u8 = 5;
-const DOWN: u8 = 12;
+const NELEC: usize = 6;
+const SIZE: usize = 8;
 
 fn compute_internal_product(
     state: FockState<u8>,
@@ -20,8 +20,11 @@ fn compute_internal_product(
     let mut pfaffian_state = construct_matrix_a_from_state(fij, state);
     let pfaffian = pfaffian_state.pfaff;
     pfaffian_state.rebuild_matrix();
+    println!("pfaffian");
     let jastrow_exp = compute_jastrow_exp(state, &vij, 8);
+    println!("jastrow");
     let gutz_exp = compute_gutzwiller_exp(state, &gi, 8);
+    println!("gutzwiller");
     let scalar_prod = <f64>::exp(jastrow_exp + gutz_exp) * pfaffian;
     scalar_prod
 }
@@ -38,15 +41,11 @@ fn main() {
     for _ in 0..SIZE {
         gi.push(rng.gen())
     }
-    let state = FockState {
-        spin_up: UP,
-        spin_down: DOWN,
-        n_sites: SIZE,
-    };
+    let state = FockState::generate_from_nelec(&mut rng, NELEC, SIZE);
     let internal_product = compute_internal_product(state, fij.clone(), vij.clone(), gi.clone());
     let rho_x = internal_product * internal_product;
-    let pot = potential(UP, DOWN);
-    let cin = kinetic(UP, DOWN);
+    let pot = potential(state.spin_up, state.spin_down);
+    let cin = kinetic(state.spin_up, state.spin_down, state.n_sites);
     println!("Terme cin to compute: {:?}", cin);
     println!("Terme pot to compute: {:?} * the internal product", pot);
     println!("rho(x): {}", rho_x);

src/parse/orbitale.rs

@@ -11,7 +11,7 @@ use std::path::PathBuf;
 /// * __`s`__ - Size of the system.
 pub fn parse_orbitale_def(fp: &PathBuf, s: usize) -> Result<Vec<f64>> {
     // Parameter reference vector. column packed.
-    let mut fij: Vec<f64> = vec![0.0; s * s];
+    let mut fij: Vec<f64> = vec![0.0; 4 * s * s];
     // Parameter vector. Ordered by "n".
     let mut orbitale_params: Vec<f64> = Vec::new();
     let mut n_orbitale_params: usize = 0;
@@ -43,11 +43,17 @@ pub fn parse_orbitale_def(fp: &PathBuf, s: usize) -> Result<Vec<f64>> {
             Ok(param) => {
                 if n_orbitale_params < param {
                     n_orbitale_params += 1;
-                    orbitale_params.push(rng.gen());
+                    orbitale_params.push(rng.gen::<f64>() * 2.0 - 1.0);
                 }
                 let var_param: f64 = orbitale_params[param - 1];
-                fij[i + j * s] = var_param;
-                fij[j + i * s] = -var_param;
+                fij[i + j * s] = 0.0;
+                fij[j + i * s] = -0.0;
+                fij[i + j * s + s*s] = var_param;
+                fij[j + i * s + s*s] = -var_param;
+                fij[i + j * s + 2*s*s] = var_param;
+                fij[j + i * s + 2*s*s] = -var_param;
+                fij[i + j * s + 3*s*s] = 0.0;
+                fij[j + i * s + 3*s*s] = -0.0;
             }
             Err(error) => {
                 // Add error message.