src/heg_solver.cc
// Find determinants connected to a det passed in.
// Return as a Wavefunction object with coefs all equal zeroes.
Find determinants connected to a det passed in.
Return as a Wavefunction object with coefs all equal zeroes.
Find determinants connected to a det passed in . Return as a Wavefunction object with coefs all equal zeroes .
std::list<Det> HEGSolver::find_connected_dets(
    const Det& det, const double eps) {
  std::list<Det> connected_dets;
  connected_dets.push_back(det);
  if (max_abs_H < eps) return connected_dets;
  const auto& k_vectors = heg.k_vectors;
  std::forward_list<IntPair> pq_pairs;
  const auto& occ_up = det.up.get_elec_orbs(n_up);
  const auto& occ_dn = det.dn.get_elec_orbs(n_dn);
  for (int i = 0; i < n_up; i++) {
    for (int j = i + 1; j < n_up; j++) {
      pq_pairs.push_front(IntPair(occ_up[i], occ_up[j]));
    }
  }
  for (int i = 0; i < n_dn; i++) {
    for (int j = i + 1; j < n_dn; j++) {
      pq_pairs.push_front(IntPair(occ_dn[i] + n_orbs, occ_dn[j] + n_orbs));
    }
  }
  for (int i = 0; i < n_up; i++) {
    for (int j = 0; j < n_dn; j++) {
      pq_pairs.push_front(IntPair(occ_up[i], occ_dn[j] + n_orbs));
    }
  }
  for (const auto& pq_pair : pq_pairs) {
    const int p = pq_pair.first;
    const int q = pq_pair.second;
    int pp = p, qq = q;
    if (p >= n_orbs && q >= n_orbs) {
      pp -= n_orbs;
      qq -= n_orbs;
    } else if (p < n_orbs && q >= n_orbs && p > q - n_orbs) {
      pp = q - n_orbs;
      qq = p + n_orbs;
    }
    bool same_spin = false;
    std::vector<Int3Double>* items_ptr;
    if (pp < n_orbs && qq < n_orbs) {
      same_spin = true;
      const auto& diff_pq = k_vectors[qq] - k_vectors[pp];
      items_ptr = &(heg.same_spin_queue[diff_pq]);
    } else {
      items_ptr = &(heg.opposite_spin_queue);
    }
    const auto& items = *items_ptr;
    int qs_offset = 0;
    if (!same_spin) qs_offset = n_orbs;
    for (const auto& item : items) {
      if (item.second < eps) break;
      const auto& diff_pr = item.first;
      int r = find_orb_id(diff_pr + k_vectors[pp]);
      if (r < 0) continue;
      int s =
          find_orb_id(k_vectors[pp] + k_vectors[qq - qs_offset] - k_vectors[r]);
      if (s < 0) continue;
      if (same_spin && s < r) continue;
      s += qs_offset;
      if (p >= n_orbs && q >= n_orbs) {
        r += n_orbs;
        s += n_orbs;
      } else if (p < n_orbs && q >= n_orbs && p > q - n_orbs) {
        const int tmp = s;
        s = r + n_orbs;
        r = tmp - n_orbs;
      }
      if (det.get_orb(r, n_orbs) || det.get_orb(s, n_orbs)) continue;
      connected_dets.push_back(det);
      Det& new_det = connected_dets.back();
      new_det.set_orb(p, n_orbs, false).set_orb(q, n_orbs, false);
      new_det.set_orb(r, n_orbs).set_orb(s, n_orbs);
    }
  }
  return connected_dets;
}
std :: list < Det > HEGSolver :: find_connected_dets ( const Det & det , const double eps ) { std :: list < Det > connected_dets ; connected_dets . push_back ( det ) ; if ( max_abs_H < eps ) return connected_dets ; const auto & k_vectors = heg . k_vectors ; std :: forward_list < IntPair > pq_pairs ; const auto & occ_up = det . up . get_elec_orbs ( n_up ) ; const auto & occ_dn = det . dn . get_elec_orbs ( n_dn ) ; for ( int i = 0 ; i < n_up ; i ++ ) { for ( int j = i + 1 ; j < n_up ; j ++ ) { pq_pairs . push_front ( IntPair ( occ_up [ i ] , occ_up [ j ] ) ) ; } } for ( int i = 0 ; i < n_dn ; i ++ ) { for ( int j = i + 1 ; j < n_dn ; j ++ ) { pq_pairs . push_front ( IntPair ( occ_dn [ i ] + n_orbs , occ_dn [ j ] + n_orbs ) ) ; } } for ( int i = 0 ; i < n_up ; i ++ ) { for ( int j = 0 ; j < n_dn ; j ++ ) { pq_pairs . push_front ( IntPair ( occ_up [ i ] , occ_dn [ j ] + n_orbs ) ) ; } } for ( const auto & pq_pair : pq_pairs ) { const int p = pq_pair . first ; const int q = pq_pair . second ; int pp = p , qq = q ; if ( p >= n_orbs && q >= n_orbs ) { pp -= n_orbs ; qq -= n_orbs ; } else if ( p < n_orbs && q >= n_orbs && p > q - n_orbs ) { pp = q - n_orbs ; qq = p + n_orbs ; } bool same_spin = false ; std :: vector < Int3Double > * items_ptr ; if ( pp < n_orbs && qq < n_orbs ) { same_spin = true ; const auto & diff_pq = k_vectors [ qq ] - k_vectors [ pp ] ; items_ptr = & ( heg . same_spin_queue [ diff_pq ] ) ; } else { items_ptr = & ( heg . opposite_spin_queue ) ; } const auto & items = * items_ptr ; int qs_offset = 0 ; if ( ! same_spin ) qs_offset = n_orbs ; for ( const auto & item : items ) { if ( item . second < eps ) break ; const auto & diff_pr = item . first ; int r = find_orb_id ( diff_pr + k_vectors [ pp ] ) ; if ( r < 0 ) continue ; int s = find_orb_id ( k_vectors [ pp ] + k_vectors [ qq - qs_offset ] - k_vectors [ r ] ) ; if ( s < 0 ) continue ; if ( same_spin && s < r ) continue ; s += qs_offset ; if ( p >= n_orbs && q >= n_orbs ) { r += n_orbs ; s += n_orbs ; } else if ( p < n_orbs && q >= n_orbs && p > q - n_orbs ) { const int tmp = s ; s = r + n_orbs ; r = tmp - n_orbs ; } if ( det . get_orb ( r , n_orbs ) || det . get_orb ( s , n_orbs ) ) continue ; connected_dets . push_back ( det ) ; Det & new_det = connected_dets . back ( ) ; new_det . set_orb ( p , n_orbs , false ) . set_orb ( q , n_orbs , false ) ; new_det . set_orb ( r , n_orbs ) . set_orb ( s , n_orbs ) ; } } return connected_dets ; }
Find determinants connected to a det passed in.
Find determinants connected to a det passed in .
// Get pq pairs. // Get rs pairs. // Test whether pqrs is a valid excitation for det.