src/heg_solver.cc
// Calculate hamiltonian between two determinants.
Calculate hamiltonian between two determinants.
Calculate hamiltonian between two determinants .
double HEGSolver::get_hamiltonian_elem(
    const Det& det_pq,
    const Det& det_rs,
    const int n_up,
    const int n_dn) const {
  auto& k_vectors = heg.k_vectors;
  const double k_unit = heg.k_unit;
  const double one_over_pi_l = 1.0 / (M_PI * heg.cell_length);
  double H = 0.0;
  if (det_pq == det_rs) {
    const auto& occ_pq_up = det_pq.up.get_elec_orbs(n_up);
    const auto& occ_pq_dn = det_pq.dn.get_elec_orbs(n_dn);
    const int n_occ_up = occ_pq_up.size();
    const int n_occ_dn = occ_pq_dn.size();
    for (const int p : occ_pq_up) {
      H += sum(square(k_vectors[p] * k_unit) * 0.5);
    }
    for (const int p : occ_pq_dn) {
      H += sum(square(k_vectors[p] * k_unit) * 0.5);
    }
    for (int i = 0; i < n_occ_up; i++) {
      const int p = occ_pq_up[i];
      for (int j = i + 1; j < n_occ_up; j++) {
        const int q = occ_pq_up[j];
        H -= one_over_pi_l / sum(square(k_vectors[p] - k_vectors[q]));
      }
    }
    for (int i = 0; i < n_occ_dn; i++) {
      const int p = occ_pq_dn[i];
      for (int j = i + 1; j < n_occ_dn; j++) {
        const int q = occ_pq_dn[j];
        H -= one_over_pi_l / sum(square(k_vectors[p] - k_vectors[q]));
      }
    }
  } else {
    Det det_eor;
    det_eor.from_eor(det_pq, det_rs);
    const int n_eor_up = det_eor.up.get_n_elecs();
    const int n_eor_dn = det_eor.dn.get_n_elecs();
    if (n_eor_up + n_eor_dn != 4) return 0.0;
    const auto& eor_up_set_bits = det_eor.up.get_elec_orbs(n_eor_up);
    const auto& eor_dn_set_bits = det_eor.dn.get_elec_orbs(n_eor_dn);
    bool k_p_set = false, k_r_set = false;
    int orb_p = 0, orb_r = 0, orb_s = 0;
    Int3 k_change;
    k_change.fill(0);
    for (const auto& orb_i : eor_up_set_bits) {
      if (det_pq.up.get_orb(orb_i)) {
        k_change -= k_vectors[orb_i];
        if (!k_p_set) {
          orb_p = orb_i;
          k_p_set = true;
        }
      } else {
        k_change += k_vectors[orb_i];
        if (!k_r_set) {
          orb_r = orb_i;
          k_r_set = true;
        } else {
          orb_s = orb_i;
        }
      }
    }
    for (const auto& orb_i : eor_dn_set_bits) {
      if (det_pq.dn.get_orb(orb_i)) {
        k_change -= k_vectors[orb_i];
        if (!k_p_set) {
          orb_p = orb_i;
          k_p_set = true;
        }
      } else {
        k_change += k_vectors[orb_i];
        if (!k_r_set) {
          orb_r = orb_i;
          k_r_set = true;
        } else {
          orb_s = orb_i;
        }
      }
    }
    if (k_change != 0) return 0.0;
    H = one_over_pi_l / sum(square(k_vectors[orb_p] - k_vectors[orb_r]));
    if (n_eor_up != 2) {
      H -= one_over_pi_l / sum(square(k_vectors[orb_p] - k_vectors[orb_s]));
    }
    int gamma_exp = get_gamma_exp(det_pq.up, n_up, eor_up_set_bits) +
                    get_gamma_exp(det_pq.dn, n_dn, eor_dn_set_bits) +
                    get_gamma_exp(det_rs.up, n_up, eor_up_set_bits) +
                    get_gamma_exp(det_rs.dn, n_dn, eor_dn_set_bits);
    if ((gamma_exp & 1) == 1) H = -H;
  }
  return H;
}
double HEGSolver :: get_hamiltonian_elem ( const Det & det_pq , const Det & det_rs , const int n_up , const int n_dn ) const { auto & k_vectors = heg . k_vectors ; const double k_unit = heg . k_unit ; const double one_over_pi_l = 1.0 / ( M_PI * heg . cell_length ) ; double H = 0.0 ; if ( det_pq == det_rs ) { const auto & occ_pq_up = det_pq . up . get_elec_orbs ( n_up ) ; const auto & occ_pq_dn = det_pq . dn . get_elec_orbs ( n_dn ) ; const int n_occ_up = occ_pq_up . size ( ) ; const int n_occ_dn = occ_pq_dn . size ( ) ; for ( const int p : occ_pq_up ) { H += sum ( square ( k_vectors [ p ] * k_unit ) * 0.5 ) ; } for ( const int p : occ_pq_dn ) { H += sum ( square ( k_vectors [ p ] * k_unit ) * 0.5 ) ; } for ( int i = 0 ; i < n_occ_up ; i ++ ) { const int p = occ_pq_up [ i ] ; for ( int j = i + 1 ; j < n_occ_up ; j ++ ) { const int q = occ_pq_up [ j ] ; H -= one_over_pi_l / sum ( square ( k_vectors [ p ] - k_vectors [ q ] ) ) ; } } for ( int i = 0 ; i < n_occ_dn ; i ++ ) { const int p = occ_pq_dn [ i ] ; for ( int j = i + 1 ; j < n_occ_dn ; j ++ ) { const int q = occ_pq_dn [ j ] ; H -= one_over_pi_l / sum ( square ( k_vectors [ p ] - k_vectors [ q ] ) ) ; } } } else { Det det_eor ; det_eor . from_eor ( det_pq , det_rs ) ; const int n_eor_up = det_eor . up . get_n_elecs ( ) ; const int n_eor_dn = det_eor . dn . get_n_elecs ( ) ; if ( n_eor_up + n_eor_dn != 4 ) return 0.0 ; const auto & eor_up_set_bits = det_eor . up . get_elec_orbs ( n_eor_up ) ; const auto & eor_dn_set_bits = det_eor . dn . get_elec_orbs ( n_eor_dn ) ; bool k_p_set = false , k_r_set = false ; int orb_p = 0 , orb_r = 0 , orb_s = 0 ; Int3 k_change ; k_change . fill ( 0 ) ; for ( const auto & orb_i : eor_up_set_bits ) { if ( det_pq . up . get_orb ( orb_i ) ) { k_change -= k_vectors [ orb_i ] ; if ( ! k_p_set ) { orb_p = orb_i ; k_p_set = true ; } } else { k_change += k_vectors [ orb_i ] ; if ( ! k_r_set ) { orb_r = orb_i ; k_r_set = true ; } else { orb_s = orb_i ; } } } for ( const auto & orb_i : eor_dn_set_bits ) { if ( det_pq . dn . get_orb ( orb_i ) ) { k_change -= k_vectors [ orb_i ] ; if ( ! k_p_set ) { orb_p = orb_i ; k_p_set = true ; } } else { k_change += k_vectors [ orb_i ] ; if ( ! k_r_set ) { orb_r = orb_i ; k_r_set = true ; } else { orb_s = orb_i ; } } } if ( k_change != 0 ) return 0.0 ; H = one_over_pi_l / sum ( square ( k_vectors [ orb_p ] - k_vectors [ orb_r ] ) ) ; if ( n_eor_up != 2 ) { H -= one_over_pi_l / sum ( square ( k_vectors [ orb_p ] - k_vectors [ orb_s ] ) ) ; } int gamma_exp = get_gamma_exp ( det_pq . up , n_up , eor_up_set_bits ) + get_gamma_exp ( det_pq . dn , n_dn , eor_dn_set_bits ) + get_gamma_exp ( det_rs . up , n_up , eor_up_set_bits ) + get_gamma_exp ( det_rs . dn , n_dn , eor_dn_set_bits ) ; if ( ( gamma_exp & 1 ) == 1 ) H = - H ; } return H ; }
Calculate hamiltonian between two determinants.
Calculate hamiltonian between two determinants .
// n_up n_dn are optional inputs for speed up get_elec_orbs. // Diagonal elements. // One electron operator. // Two electrons operator. // Off-diagonal elements. // Obtain p, q, s. // Check for momentum conservation.