jscaux
/
ABACUS
2
0
Fork 0
Code Issues 1 Pull Requests Releases Wiki Activity
ABACUS/src/SCAN/Descendents.cc

1189 lines
52 KiB
C++
Raw Permalink Blame History

/**********************************************************
This software is part of J.-S. Caux's ABACUS library.
Copyright (c) J.-S. Caux.
-----------------------------------------------------------
File: Descendents.cc
Purpose: New implementation of state scanning:
functions to descend down hierarchy of states.
***********************************************************/
#include "ABACUS.h"
using namespace std;
using namespace ABACUS;
namespace ABACUS {
Vect<string> Descendent_States_with_iK_Stepped_Up
(string ScanIx2_label, const Vect<Vect<int> >& OriginIx2, const Vect<Vect<int> >& BaseScanIx2,
const Vect<int>& Ix2_min, const Vect<int>& Ix2_max, bool disperse_only_current_exc, bool preserve_nexc)
{
// Given an OriginIx2 and a BaseScanIx2 (by which we mean a set of Ix2 providing a basis for scanning),
// this function returns a vector of the labels of the states obtained by each allowable
// one-step increase of the quantum numbers.
// Rules for moving quantum numbers:
// - quantum numbers preserve their ordering (an Umklapp is thus made of many one-step movements)
// - if a quantum number has moved as compared to BaseScanIx2, then it can only move further in the same direction
// The ordering rule for allowable increases is that we start from the topmost, rightmost quantum number.
// The increases can thus be applied to quantum numbers below the bottom-most, left-most right-displaced quantum nr.
// The preserve_nexc boolean flags whether we maintain or increase the nr of ph pairs.
// This is only checked if disperse_only_current_exc == false (otherwise it's irrelevant).
Vect<Vect<int> > ScanIx2 = Return_Ix2_from_Label (ScanIx2_label, OriginIx2);
// Determine the level and index of the bottom-most left-most right-moving quantum number sits:
int exclevel = -1;
int excindex = 0;
bool excfound = false;
do {
exclevel++;
if (exclevel == ScanIx2.size()) { // there isn't a single right-moving quantum number in ScanIx2
break;
}
for (int alpha = 0; alpha < ScanIx2[exclevel].size(); ++alpha)
if (ScanIx2[exclevel][alpha] > BaseScanIx2[exclevel][alpha]) {
excindex = alpha;
excfound = true;
break;
}
} while (!excfound);
// If we haven't found an excitation, then exclevel == ScanIx2.size() and excindex = 0;
// The quantum numbers which we can move right are thus those
// with (j < exclevel) and (j == exclelev and alpha <= excindex)
int ndesc_possible = 1;
if (!disperse_only_current_exc) {
ndesc_possible = 0;
for (int j = 0; j <= ABACUS::min(exclevel, ScanIx2.size() -1); ++j) ndesc_possible += ScanIx2[j].size();
}
// Now construct the actual descendents:
Vect<string> desclabelfound (ndesc_possible);
int ndesc_found = 0;
if (!disperse_only_current_exc && !preserve_nexc)
// List descendents with a new excitation at a lower level:
// this can only create an additional ph pair
for (int j = 0; j < exclevel; ++j)
for (int alpha = 0; alpha < ScanIx2[j].size(); ++alpha) {
if (ScanIx2[j][alpha] >= BaseScanIx2[j][alpha] // candidate either undisplaced or right-moving
&& !ScanIx2[j].includes(ScanIx2[j][alpha] + 2) // there is a hole to its right
&& ScanIx2[j][alpha] + 2 <= Ix2_max[j]) { // this new quantum number fits within the limits
// We've found a descendent
ScanIx2[j][alpha] += 2;
desclabelfound[ndesc_found] = Return_State_Label (ScanIx2, OriginIx2);
ScanIx2[j][alpha] -= 2;
ndesc_found++;
}
}
// List descendents with a new excitation at exclevel:
if (exclevel < ScanIx2.size()) { // excfound == true, excindex is now guaranteed < ScanIx2[exclevel].size()
int alphamin = (disperse_only_current_exc ? excindex : 0);
int alphamax = (disperse_only_current_exc ? excindex : excindex - 1);
//for (int alpha = 0; alpha <= excindex; ++alpha) {
for (int alpha = alphamin; alpha <= alphamax; ++alpha) {
if (ScanIx2[exclevel][alpha] >= BaseScanIx2[exclevel][alpha]
&& !ScanIx2[exclevel].includes(ScanIx2[exclevel][alpha] + 2)
&& ScanIx2[exclevel][alpha] + 2 <= Ix2_max[exclevel]) {
// We've found a descendent
ScanIx2[exclevel][alpha] += 2;
desclabelfound[ndesc_found] = Return_State_Label (ScanIx2, OriginIx2);
ScanIx2[exclevel][alpha] -= 2;
// If we're dispersing a subleading Ix2, check whether we match the preserve_nexc condition:
// The descendent is acceptable if disperse_only_current_exc == true,
// or if preserve_nexc == true and nexc labels match,
// or if preserve_nexc == false and nexc labels don't match:
if (disperse_only_current_exc
//|| (preserve_nexc == BaseScanIx2[exclevel].includes(ScanIx2[exclevel][alpha] + 2)))
|| (preserve_nexc == (Extract_nexc_Label(desclabelfound[ndesc_found]
).compare(Extract_nexc_Label(ScanIx2_label)) == 0)))
ndesc_found++;
}
}
} //if (exclevel < ScanIx2.size())
// Resize desc:
Vect<string> desclabelfound_resized(ndesc_found);
for (int idesc = 0; idesc < ndesc_found; ++idesc)
desclabelfound_resized[idesc] = desclabelfound[idesc];
return(desclabelfound_resized);
}
Vect<string> Descendent_States_with_iK_Stepped_Down
(string ScanIx2_label, const Vect<Vect<int> >& OriginIx2, const Vect<Vect<int> >& BaseScanIx2,
const Vect<int>& Ix2_min, const Vect<int>& Ix2_max, bool disperse_only_current_exc, bool preserve_nexc)
{
// Given an OriginIx2 and a BaseScanIx2 (by which we mean a set of Ix2 providing a basis for scanning),
// this function returns a vector of the labels of the states obtained by each allowable
// one-step decrease of the quantum numbers.
// Rules for moving quantum numbers:
// - quantum numbers preserve their ordering (an Umklapp is thus made of many one-step movements)
// - if a quantum number has moved as compared to BaseScanIx2, then it can only move further in the same direction
// The ordering rule for allowable decreases is that we start from the topmost, leftmost quantum number.
// The decreases can thus be applied to quantum numbers below the bottom-most, right-most left-displaced quantum nr.
// The preserve_nexc boolean flags whether we maintain or increase the nr of ph pairs.
// This is only checked if disperse_only_current_exc == false (otherwise it's irrelevant).
Vect<Vect<int> > ScanIx2 = Return_Ix2_from_Label (ScanIx2_label, OriginIx2);
// Determine the level and index of the bottom-most right-most left-moving quantum number sits:
int exclevel = -1;
int excindex = 0;
bool excfound = false;
do {
exclevel++;
if (exclevel == ScanIx2.size()) { // there isn't a single left-moving quantum number in ScanIx2
break;
}
for (int alpha = ScanIx2[exclevel].size() - 1; alpha >= 0; --alpha) {
if (ScanIx2[exclevel][alpha] < BaseScanIx2[exclevel][alpha]) {
excindex = alpha;
excfound = true;
break;
}
}
} while (!excfound);
// If we haven't found an excitation, then exclevel == ScanIx2.size() and excindex = 0;
// The quantum numbers which we can move left are thus those
// with (j < exclevel) and (j == exclelev and alpha >= excindex)
int ndesc_possible = 1;
if (!disperse_only_current_exc) {
ndesc_possible = 0;
for (int j = 0; j <= ABACUS::min(exclevel, ScanIx2.size() - 1); ++j) ndesc_possible += ScanIx2[j].size();
}
// Now construct the actual descendents:
Vect<string> desclabelfound (ndesc_possible);
int ndesc_found = 0;
if (!disperse_only_current_exc && !preserve_nexc)
// List descendents with a new excitation at a lower level:
for (int j = 0; j < exclevel; ++j)
for (int alpha = 0; alpha < ScanIx2[j].size(); ++alpha) {
if (ScanIx2[j][alpha] <= BaseScanIx2[j][alpha] // candidate either undisplaced or left-moving
&& !ScanIx2[j].includes(ScanIx2[j][alpha] - 2) // there is a hole to its left
&& ScanIx2[j][alpha] - 2 >= Ix2_min[j]) { // this new quantum number fits within the limits
// We've found a descendent
ScanIx2[j][alpha] -= 2;
desclabelfound[ndesc_found] = Return_State_Label (ScanIx2, OriginIx2);
ScanIx2[j][alpha] += 2;
ndesc_found++;
}
}
// List descendents with a new excitation at exclevel:
if (exclevel < ScanIx2.size()) { // excfound == true, excindex is now guaranteed < ScanIx2[exclevel].size()
int alphamin = (disperse_only_current_exc ? excindex : excindex + 1);
int alphamax = (disperse_only_current_exc ? excindex : ScanIx2[exclevel].size() - 1);
for (int alpha = alphamax; alpha >= alphamin; --alpha) {
if (ScanIx2[exclevel][alpha] <= BaseScanIx2[exclevel][alpha]
&& !ScanIx2[exclevel].includes(ScanIx2[exclevel][alpha] - 2)
&& ScanIx2[exclevel][alpha] - 2 >= Ix2_min[exclevel]) {
ScanIx2[exclevel][alpha] -= 2;
desclabelfound[ndesc_found] = Return_State_Label (ScanIx2, OriginIx2);
ScanIx2[exclevel][alpha] += 2;
if (disperse_only_current_exc
|| (preserve_nexc == (Extract_nexc_Label(desclabelfound[ndesc_found]
).compare(Extract_nexc_Label(ScanIx2_label)) == 0)))
ndesc_found++;
}
}
} // if (exclevel < ScanIx2.size())
// Resize desc:
Vect<string> desclabelfound_resized(ndesc_found);
for (int idesc = 0; idesc < ndesc_found; ++idesc)
desclabelfound_resized[idesc] = desclabelfound[idesc];
return(desclabelfound_resized);
}
Vect<string> Descendent_States_with_iK_Preserved
(string ScanIx2_label, const Vect<Vect<int> >& OriginIx2, const Vect<Vect<int> >& BaseScanIx2,
const Vect<int>& Ix2_min, const Vect<int>& Ix2_max,
bool disperse_only_current_exc_up, bool preserve_nexc_up,
bool disperse_only_current_exc_down, bool preserve_nexc_down)
{
// Returns the labels of all states which are obtained from ScanIx2 by stepping p one step, and down one step.
Vect<string> labels_up = Descendent_States_with_iK_Stepped_Up
(ScanIx2_label, OriginIx2, BaseScanIx2, Ix2_min, Ix2_max, disperse_only_current_exc_up, preserve_nexc_up);
Vect<string> labels_found(0);
for (int i = 0; i < labels_up.size(); ++i) {
labels_found.Append (Descendent_States_with_iK_Stepped_Down
(labels_up[i], OriginIx2, BaseScanIx2, Ix2_min, Ix2_max,
disperse_only_current_exc_down, preserve_nexc_down));
}
return(labels_found);
}
// For symmetric state descending:
Vect<string> Descendent_States_with_iK_Stepped_Up_rightIx2only
(string ScanIx2_label, const Vect<Vect<int> >& OriginIx2, const Vect<Vect<int> >& BaseScanIx2,
const Vect<int>& Ix2_min, const Vect<int>& Ix2_max, bool disperse_only_current_exc, bool preserve_nexc)
{
// This function raises an Ix2 in the right half of ScanIx2, and does the symmetric lowering in the lower half.
// ASSUMPTIONS: OriginIx2 is a symmetric state.
Vect<Vect<int> > ScanIx2 = Return_Ix2_from_Label (ScanIx2_label, OriginIx2);
// Determine the level and index of the bottom-most left-most right-moving quantum number sits:
int exclevel = -1;
int excindex = 0;
bool excfound = false;
do {
exclevel++;
if (exclevel == ScanIx2.size()) { // there isn't a single right-moving quantum number in ScanIx2
break;
}
for (int alpha = (ScanIx2[exclevel].size() + 1)/2; alpha < ScanIx2[exclevel].size(); ++alpha)
if (ScanIx2[exclevel][alpha] > BaseScanIx2[exclevel][alpha]) {
excindex = alpha;
excfound = true;
break;
}
} while (!excfound);
// If we haven't found an excitation, then exclevel == ScanIx2.size() and excindex = 0;
// The quantum numbers which we can move right are thus those
// with (j < exclevel) and (j == exclelev and alpha <= excindex)
int ndesc_possible = 1;
if (!disperse_only_current_exc) {
ndesc_possible = 0;
for (int j = 0; j <= ABACUS::min(exclevel, ScanIx2.size() -1); ++j) ndesc_possible += ScanIx2[j].size();
}
// Now construct the actual descendents:
Vect<string> desclabelfound (ndesc_possible);
int ndesc_found = 0;
if (!disperse_only_current_exc && !preserve_nexc)
// List descendents with a new excitation at a lower level:
// this can only create an additional ph pair
for (int j = 0; j < exclevel; ++j)
for (int alpha = (ScanIx2[j].size() + 1)/2; alpha < ScanIx2[j].size(); ++alpha) {
if (ScanIx2[j][alpha] >= BaseScanIx2[j][alpha] // candidate either undisplaced or right-moving
&& !ScanIx2[j].includes(ScanIx2[j][alpha] + 2) // there is a hole to its right
&& ScanIx2[j][alpha] + 2 <= Ix2_max[j]) { // this new quantum number fits within the limits
// We've found a descendent
ScanIx2[j][alpha] += 2;
ScanIx2[j][ScanIx2[j].size() - 1 - alpha] -= 2;
desclabelfound[ndesc_found] = Return_State_Label (ScanIx2, OriginIx2);
ScanIx2[j][alpha] -= 2;
ScanIx2[j][ScanIx2[j].size() - 1 - alpha] += 2;
ndesc_found++;
}
}
// List descendents with a new excitation at exclevel:
if (exclevel < ScanIx2.size()) { // excfound == true, excindex is now guaranteed < ScanIx2[exclevel].size()
int alphamin = (disperse_only_current_exc ? excindex : 0);
alphamin = ABACUS::max((ScanIx2[exclevel].size() + 1)/2, alphamin);
int alphamax = (disperse_only_current_exc ? excindex : excindex - 1);
for (int alpha = alphamin; alpha <= alphamax; ++alpha) {
if (ScanIx2[exclevel][alpha] >= BaseScanIx2[exclevel][alpha]
&& !ScanIx2[exclevel].includes(ScanIx2[exclevel][alpha] + 2)
&& ScanIx2[exclevel][alpha] + 2 <= Ix2_max[exclevel]) {
// We've found a descendent
ScanIx2[exclevel][alpha] += 2;
ScanIx2[exclevel][ScanIx2[exclevel].size() - 1 - alpha] -= 2;
desclabelfound[ndesc_found] = Return_State_Label (ScanIx2, OriginIx2);
ScanIx2[exclevel][alpha] -= 2;
ScanIx2[exclevel][ScanIx2[exclevel].size() - 1 - alpha] += 2;
// If we're dispersing a subleading Ix2, check whether we match the preserve_nexc condition:
// The descendent is acceptable if disperse_only_current_exc == true,
// or if preserve_nexc == true and nexc labels match,
// or if preserve_nexc == false and nexc labels don't match:
if (disperse_only_current_exc
|| (preserve_nexc == (Extract_nexc_Label(desclabelfound[ndesc_found]
).compare(Extract_nexc_Label(ScanIx2_label)) == 0)))
ndesc_found++;
}
}
} //if (exclevel < ScanIx2.size())
// Resize desc:
Vect<string> desclabelfound_resized(ndesc_found);
for (int idesc = 0; idesc < ndesc_found; ++idesc)
desclabelfound_resized[idesc] = desclabelfound[idesc];
return(desclabelfound_resized);
}
Vect<string> Descendent_States_with_iK_Stepped_Down_rightIx2only
(string ScanIx2_label, const Vect<Vect<int> >& OriginIx2, const Vect<Vect<int> >& BaseScanIx2,
const Vect<int>& Ix2_min, const Vect<int>& Ix2_max, bool disperse_only_current_exc, bool preserve_nexc)
{
Vect<Vect<int> > ScanIx2 = Return_Ix2_from_Label (ScanIx2_label, OriginIx2);
// Determine the level and index of the bottom-most right-most left-moving quantum number sits:
int exclevel = -1;
int excindex = 0;
bool excfound = false;
do {
exclevel++;
if (exclevel == ScanIx2.size()) { // there isn't a single left-moving quantum number in ScanIx2
break;
}
for (int alpha = ScanIx2[exclevel].size() - 1; alpha >= (ScanIx2[exclevel].size() + 1)/2 ; --alpha) {
if (ScanIx2[exclevel][alpha] < BaseScanIx2[exclevel][alpha]) {
excindex = alpha;
excfound = true;
break;
}
}
} while (!excfound);
// If we haven't found an excitation, then exclevel == ScanIx2.size() and excindex = 0;
// The quantum numbers which we can move left are thus those
// with (j < exclevel) and (j == exclevel and alpha >= excindex)
int ndesc_possible = 1;
if (!disperse_only_current_exc) {
ndesc_possible = 0;
for (int j = 0; j <= ABACUS::min(exclevel, ScanIx2.size() - 1); ++j) ndesc_possible += ScanIx2[j].size();
}
// Now construct the actual descendents:
Vect<string> desclabelfound (ndesc_possible);
int ndesc_found = 0;
if (!disperse_only_current_exc && !preserve_nexc)
// List descendents with a new excitation at a lower level:
for (int j = 0; j < exclevel; ++j)
for (int alpha = (ScanIx2[j].size() + 1)/2; alpha < ScanIx2[j].size(); ++alpha) {
if (ScanIx2[j][alpha] <= BaseScanIx2[j][alpha] // candidate either undisplaced or left-moving
&& !ScanIx2[j].includes(ScanIx2[j][alpha] - 2) // there is a hole to its left
&& ScanIx2[j][alpha] - 2 >= Ix2_min[j]) { // this new quantum number fits within the limits
// We've found a descendent
ScanIx2[j][alpha] -= 2;
ScanIx2[j][ScanIx2[j].size() - 1 - alpha] += 2;
desclabelfound[ndesc_found] = Return_State_Label (ScanIx2, OriginIx2);
ScanIx2[j][alpha] += 2;
ScanIx2[j][ScanIx2[j].size() - 1 - alpha] -= 2;
ndesc_found++;
}
}
// List descendents with a new excitation at exclevel:
if (exclevel < ScanIx2.size()) { // excfound == true, excindex is now guaranteed < ScanIx2[exclevel].size()
int alphamin = (disperse_only_current_exc ? excindex : excindex + 1);
alphamin = ABACUS::max((ScanIx2[exclevel].size() + 1)/2, alphamin);
int alphamax = (disperse_only_current_exc ? excindex : ScanIx2[exclevel].size() - 1);
for (int alpha = alphamax; alpha >= alphamin; --alpha) {
if (ScanIx2[exclevel][alpha] <= BaseScanIx2[exclevel][alpha]
&& !ScanIx2[exclevel].includes(ScanIx2[exclevel][alpha] - 2)
&& ScanIx2[exclevel][alpha] - 2 >= Ix2_min[exclevel]) {
ScanIx2[exclevel][alpha] -= 2;
ScanIx2[exclevel][ScanIx2[exclevel].size() - 1 - alpha] += 2;
desclabelfound[ndesc_found] = Return_State_Label (ScanIx2, OriginIx2);
ScanIx2[exclevel][alpha] += 2;
ScanIx2[exclevel][ScanIx2[exclevel].size() - 1 - alpha] -= 2;
if (disperse_only_current_exc
|| (preserve_nexc == (Extract_nexc_Label(desclabelfound[ndesc_found]
).compare(Extract_nexc_Label(ScanIx2_label)) == 0)))
ndesc_found++;
}
}
} // if (exclevel < ScanIx2.size())
// Resize desc:
Vect<string> desclabelfound_resized(ndesc_found);
for (int idesc = 0; idesc < ndesc_found; ++idesc)
desclabelfound_resized[idesc] = desclabelfound[idesc];
return(desclabelfound_resized);
}
// Specialization for Lieb-Liniger case:
Vect<string> Descendent_States_with_iK_Stepped_Up (string ScanIx2_label, const LiebLin_Bethe_State& OriginState,
bool disperse_only_current_exc, bool preserve_nexc)
{
Vect<Vect<int> > OriginIx2here(1);
OriginIx2here[0] = OriginState.Ix2;
Vect<Vect<int> > BaseScanIx2here(1);
BaseScanIx2here[0] = OriginState.Ix2;
Vect<int> Ix2_min(1);
Ix2_min[0] = LIEBLIN_Ix2_MIN;
Vect<int> Ix2_max(1);
Ix2_max[0] = LIEBLIN_Ix2_MAX;
return (Descendent_States_with_iK_Stepped_Up (ScanIx2_label, OriginIx2here, BaseScanIx2here, Ix2_min, Ix2_max,
disperse_only_current_exc, preserve_nexc));
}
// Specialization for Lieb-Liniger case:
Vect<string> Descendent_States_with_iK_Stepped_Down (string ScanIx2_label, const LiebLin_Bethe_State& OriginState,
bool disperse_only_current_exc, bool preserve_nexc)
{
Vect<Vect<int> > OriginIx2here(1);
OriginIx2here[0] = OriginState.Ix2;
Vect<Vect<int> > BaseScanIx2here(1);
BaseScanIx2here[0] = OriginState.Ix2;
Vect<int> Ix2_min(1);
Ix2_min[0] = LIEBLIN_Ix2_MIN;
Vect<int> Ix2_max(1);
Ix2_max[0] = LIEBLIN_Ix2_MAX;
return (Descendent_States_with_iK_Stepped_Down (ScanIx2_label, OriginIx2here, BaseScanIx2here, Ix2_min, Ix2_max,
disperse_only_current_exc, preserve_nexc));
}
// Specialization for Lieb-Liniger case:
Vect<string> Descendent_States_with_iK_Preserved (string ScanIx2_label, const LiebLin_Bethe_State& OriginState,
bool disperse_only_current_exc_up, bool preserve_nexc_up,
bool disperse_only_current_exc_down, bool preserve_nexc_down)
{
Vect<Vect<int> > OriginIx2here(1);
OriginIx2here[0] = OriginState.Ix2;
Vect<Vect<int> > BaseScanIx2here(1);
BaseScanIx2here[0] = OriginState.Ix2;
Vect<int> Ix2_min(1);
Ix2_min[0] = LIEBLIN_Ix2_MIN;
Vect<int> Ix2_max(1);
Ix2_max[0] = LIEBLIN_Ix2_MAX;
return (Descendent_States_with_iK_Preserved (ScanIx2_label, OriginIx2here, BaseScanIx2here, Ix2_min, Ix2_max,
disperse_only_current_exc_up, preserve_nexc_up,
disperse_only_current_exc_down, preserve_nexc_down));
}
// Specialization for Lieb-Liniger case:
Vect<string> Descendent_States_with_iK_Stepped_Up_rightIx2only
(string ScanIx2_label, const LiebLin_Bethe_State& OriginState, bool disperse_only_current_exc, bool preserve_nexc)
{
Vect<Vect<int> > OriginIx2here(1);
OriginIx2here[0] = OriginState.Ix2;
Vect<Vect<int> > BaseScanIx2here(1);
BaseScanIx2here[0] = OriginState.Ix2;
Vect<int> Ix2_min(1);
Ix2_min[0] = LIEBLIN_Ix2_MIN;
Vect<int> Ix2_max(1);
Ix2_max[0] = LIEBLIN_Ix2_MAX;
return(Descendent_States_with_iK_Stepped_Up_rightIx2only (ScanIx2_label, OriginIx2here, BaseScanIx2here,
Ix2_min, Ix2_max, disperse_only_current_exc, preserve_nexc));
}
Vect<string> Descendent_States_with_iK_Stepped_Down_rightIx2only
(string ScanIx2_label, const LiebLin_Bethe_State& OriginState, bool disperse_only_current_exc, bool preserve_nexc)
{
Vect<Vect<int> > OriginIx2here(1);
OriginIx2here[0] = OriginState.Ix2;
Vect<Vect<int> > BaseScanIx2here(1);
BaseScanIx2here[0] = OriginState.Ix2;
Vect<int> Ix2_min(1);
Ix2_min[0] = LIEBLIN_Ix2_MIN;
Vect<int> Ix2_max(1);
Ix2_max[0] = LIEBLIN_Ix2_MAX;
return(Descendent_States_with_iK_Stepped_Down_rightIx2only (ScanIx2_label, OriginIx2here, BaseScanIx2here,
Ix2_min, Ix2_max, disperse_only_current_exc, preserve_nexc));
}
// Specializations for Heis states:
Vect<string> Descendent_States_with_iK_Stepped_Up
(string ScanIx2_label, const Heis_Bethe_State& OriginState, bool disperse_only_current_exc, bool preserve_nexc)
{
return(Descendent_States_with_iK_Stepped_Up (ScanIx2_label, OriginState.Ix2, OriginState.Ix2, OriginState.base.Ix2_min,
OriginState.base.Ix2_max, disperse_only_current_exc, preserve_nexc));
}
Vect<string> Descendent_States_with_iK_Stepped_Down
(string ScanIx2_label, const Heis_Bethe_State& OriginState, bool disperse_only_current_exc, bool preserve_nexc)
{
return(Descendent_States_with_iK_Stepped_Down (ScanIx2_label, OriginState.Ix2, OriginState.Ix2, OriginState.base.Ix2_min,
OriginState.base.Ix2_max, disperse_only_current_exc, preserve_nexc));
}
Vect<string> Descendent_States_with_iK_Preserved
(string ScanIx2_label, const Heis_Bethe_State& OriginState, bool disperse_only_current_exc_up, bool preserve_nexc_up,
bool disperse_only_current_exc_down, bool preserve_nexc_down)
{
return(Descendent_States_with_iK_Preserved
(ScanIx2_label, OriginState.Ix2, OriginState.Ix2,
OriginState.base.Ix2_min, OriginState.base.Ix2_max, disperse_only_current_exc_up, preserve_nexc_up,
disperse_only_current_exc_down, preserve_nexc_down));
}
Vect<string> Descendent_States_with_iK_Stepped_Up_rightIx2only
(string ScanIx2_label, const Heis_Bethe_State& OriginState, bool disperse_only_current_exc, bool preserve_nexc)
{
return(Descendent_States_with_iK_Stepped_Up_rightIx2only
(ScanIx2_label, OriginState.Ix2, OriginState.Ix2, OriginState.base.Ix2_min, OriginState.base.Ix2_max,
disperse_only_current_exc, preserve_nexc));
}
Vect<string> Descendent_States_with_iK_Stepped_Down_rightIx2only
(string ScanIx2_label, const Heis_Bethe_State& OriginState, bool disperse_only_current_exc, bool preserve_nexc)
{
return(Descendent_States_with_iK_Stepped_Down_rightIx2only
(ScanIx2_label, OriginState.Ix2, OriginState.Ix2, OriginState.base.Ix2_min, OriginState.base.Ix2_max,
disperse_only_current_exc, preserve_nexc));
}
Vect<bool> Is_Good_New_Hole_Position (const Vect<Vect<int> >& OriginIx2, State_Label_Data currentdata, int exclevel_newph)
{
// Given a state, returns the acceptable new hole positions.
// Define the objects for the newstatedata:
Vect<int> type_new = currentdata.type;
Vect<int> M_new = currentdata.M;
Vect<int> nexc_new = currentdata.nexc;
nexc_new[exclevel_newph] += 1; // we drill one more particle-hole pair at this level
int index_new = (currentdata.nexc[exclevel_newph] + 1)/2; // we put the new p-h pair at index index_new.
int ntypespresent = currentdata.type.size();
Vect<Vect<int> > Ix2old_new(ntypespresent);
Vect<Vect<int> > Ix2exc_new(ntypespresent);
for (int it = 0; it < ntypespresent; ++it) Ix2old_new[it] = Vect<int>(ABACUS::max(nexc_new[it],1));
for (int it = 0; it < ntypespresent; ++it) Ix2exc_new[it] = Vect<int>(ABACUS::max(nexc_new[it],1));
// Copy earlier data in:
for (int it = 0; it < ntypespresent; ++it) {
for (int i = 0; i < currentdata.nexc[it]; ++i) {
Ix2old_new[it][i + (it == exclevel_newph && i >= index_new)] = currentdata.Ix2old[it][i];
Ix2exc_new[it][i + (it == exclevel_newph && i >= index_new)] = currentdata.Ix2exc[it][i];
}
}
State_Label_Data descdatanewph (type_new, M_new, nexc_new, Ix2old_new, Ix2exc_new);
// We now look for all possible hole positions,
// the allowable ones being either at the edge of a block, or next to an existing hole excitation,
// under the condition of obeying the `towards the block center' rule.
Vect<bool> isgoodnewholepos(false, OriginIx2[descdatanewph.type[exclevel_newph] ].size());
for (int ih = 0; ih < OriginIx2[descdatanewph.type[exclevel_newph] ].size(); ++ih) {
int candidateIx2old = OriginIx2[descdatanewph.type[exclevel_newph] ][ih];
// candidateIx2old is an acceptable position for the new hole provided the following conditions are fulfilled:
// A- it is in OriginIx2
// B- it follows the ordering rule, i.e. it sits in the middle of previous particle excitations,
// namely between Ix2old[index_new - 1] (if this exists) and Ix2old[index_new] (if this exists)
// C- it is next to an OriginIx2 vacancy or immediately right of Ix2old[index_new - 1] (if this exists)
// or immediately left of Ix2old[index_new] (if this exists)
// D- it does not break the `towards the center' rule
// (it will break the rule at this point if it is created away from an OriginIx2 boundary
// (and thus next to a preexisting excitation),
// and if this excitation and it are not in the same sideblock
// (in other words: if there is a sideblock boundary between them)
// A
if (!OriginIx2[descdatanewph.type[exclevel_newph] ].includes (candidateIx2old)) {
// is not contained in OriginIx2
continue;
}
// B1
if (currentdata.nexc[exclevel_newph] > 0
&& candidateIx2old <= currentdata.Ix2old[exclevel_newph][index_new - 1]) {
// there is at least one hole exc to the left, and the candidate position isn't right of Ix2old[index_new - 1]
continue;
}
// B2
if (currentdata.nexc[exclevel_newph] > 1
&& candidateIx2old >= currentdata.Ix2old[exclevel_newph][index_new]) {
// there is at least one hole exc to the right, and the candidate position isn't left of Ix2old[index_new]
continue;
}
// C- it is next to an OriginIx2 vacancy or immediately right of Ix2old[index_new - 1] (if this exists)
// or immediately left of Ix2old[index_new] (if this exists)
if (!(!OriginIx2[descdatanewph.type[exclevel_newph] ].includes (candidateIx2old + 2)
|| !OriginIx2[descdatanewph.type[exclevel_newph] ].includes (candidateIx2old - 2))
// doesn't sit next to an OriginIx2 vacancy
&& (currentdata.nexc[exclevel_newph] == 0
|| candidateIx2old != currentdata.Ix2old[exclevel_newph][index_new - 1] + 2)
// doesn't sit immediately right of first hole excitation to the left
&& (currentdata.nexc[exclevel_newph] <= 1
|| candidateIx2old != currentdata.Ix2old[exclevel_newph][index_new] - 2)
// doesn't sit immediately left of first hole excitation to the right
) {
continue;
}
// D- it does not break the `towards the center' rule
// In other words, if created away from a block boundary but next to a preexisting hole,
// must be in same sideblock as this particle:
// Determine the size of the block of OriginIx2 in which this hole sits:
int nroccupiedtoleft = 0;
while (OriginIx2[descdatanewph.type[exclevel_newph] ].includes (candidateIx2old -2*(nroccupiedtoleft + 1)))
nroccupiedtoleft++;
int nroccupiedtoright = 0;
while (OriginIx2[descdatanewph.type[exclevel_newph] ].includes (candidateIx2old +2*(nroccupiedtoright + 1)))
nroccupiedtoright++;
// We can determine whether the new hole would be left- or right-moving
bool hole_candidate_is_left_moving = nroccupiedtoleft >= nroccupiedtoright;
bool hole_candidate_is_right_moving = nroccupiedtoleft < nroccupiedtoright;
if (hole_candidate_is_left_moving
&& (currentdata.nexc[exclevel_newph] > 0 && candidateIx2old == currentdata.Ix2old[exclevel_newph][index_new - 1] + 2)
// it is created to the right of a preexisting hole excitation
&& OriginIx2[descdatanewph.type[exclevel_newph] ].includes (candidateIx2old + 2)
// and is not sitting at the boundary
&& (currentdata.nexc[exclevel_newph] <= 1 || candidateIx2old != currentdata.Ix2old[exclevel_newph][index_new] - 2)
// and is not sitting just left of another preexisting hole excitation which is also left moving
)
{
continue;
}
if (hole_candidate_is_right_moving
&& (currentdata.nexc[exclevel_newph] > 1 && candidateIx2old == currentdata.Ix2old[exclevel_newph][index_new] - 2)
// it is created to the left of a preexisting hole excitation
&& OriginIx2[descdatanewph.type[exclevel_newph] ].includes (candidateIx2old - 2)
// and is not sitting at the boundary
&& (currentdata.nexc[exclevel_newph] == 0
|| candidateIx2old != currentdata.Ix2old[exclevel_newph][index_new - 1] + 2)
// and is not sitting just right of another preexisting hole excitation which is also right moving
)
{
continue;
}
// If we have reached this point, candidateIx2old is acceptable.
isgoodnewholepos[ih] = true;
} // for (int ih
return(isgoodnewholepos);
} // Vect<bool> Is_Good_New_Hole_Position
Vect<string> Descendents (string ScanStateLabel, const Vect<Vect<int> >& ScanIx2, const Vect<Vect<int> >& OriginIx2,
const Vect<int>& Ix2_min, const Vect<int>& Ix2_max, int type_required)
{
// This generates descendents for a fixed base.
// As compared to the OriginState, ScanState thus has only (in-level) particle-hole excitations
// ASSUMPTIONS:
// ScanStateLabel is consistent with ScanIx2 and OriginIx2.
// type_required == 0: inner (right)most hole moved
// type_required == 1: inner (right)most particle moved
// type_required == 2: new particle-hole pair added, keeping only closest p < h and p > h cases
// type_required == 3: at level 0, move leftmost particle one step left, rightmost one step right
// (fixed iK logic using skeletons).
// type_required == 4: generalized Umklapp: increase distance between latest-added p-h pair,
// staying on boundary of blocks
bool cinbreaks = false;
// Number of descendents:
int ndesc_possible = 1;
if (type_required == 1 || type_required == 2) {
for (int i = 0; i < ScanIx2.size(); ++i)
ndesc_possible = ABACUS::max (ndesc_possible, 2* ScanIx2[i].size() * 2 * ScanIx2[i].size());
// inexact, should be refined
}
Vect<string> desclabelfound (ndesc_possible);
Vect<int> desctypefound (ndesc_possible);
int ndesc_found = 0;
// First read the current state data:
State_Label_Data currentdata = Read_State_Label (ScanStateLabel, OriginIx2);
// Determine the level at which the highest current excitation sits
int exclevel = currentdata.type.size() - 1;
while (exclevel > 0 && currentdata.nexc[exclevel] == 0) exclevel--;
if (type_required == 0) {
// We move the inner(right)most hole of the highest excited level one more step towards
// the center of the block of Ix2 vacancies in which it sits.
// Is there already an excitation at level 0? If so, move inner(right)most particle one step further.
if (currentdata.nexc[0] > 0) {
// Produce a new state_label_data:
State_Label_Data descdata = Read_State_Label (ScanStateLabel, OriginIx2);
// Identify the inner(right)most excitation:
int innerindex = currentdata.nexc[exclevel]/2;
int newholepos = currentdata.Ix2old[exclevel][innerindex];
// Determine the size of the block of OriginState.Ix2 in which this hole sits:
int nroccupiedtoleft = 0;
while (OriginIx2[descdata.type[exclevel] ].includes (currentdata.Ix2old[exclevel][innerindex]
-2*(nroccupiedtoleft + 1)))
nroccupiedtoleft++;
int nroccupiedtoright = 0;
while (OriginIx2[descdata.type[exclevel] ].includes (currentdata.Ix2old[exclevel][innerindex]
+2*(nroccupiedtoright + 1)))
nroccupiedtoright++;
if (nroccupiedtoleft - nroccupiedtoright + 2 < 0) newholepos += 2;
else if (nroccupiedtoleft - nroccupiedtoright - 2 >= 0) newholepos -= 2;
if (newholepos != currentdata.Ix2old[exclevel][innerindex] // we have successfully moved the hole
&& !currentdata.Ix2old[exclevel].includes(newholepos) // new hole position is not already taken
)
{ // we have found a descendent
descdata.Ix2old[exclevel][innerindex] = newholepos;
desclabelfound[ndesc_found] = Return_State_Label (descdata, OriginIx2);
desctypefound[ndesc_found] = 0;
ndesc_found++;
if (cinbreaks) { char a; cin >> a;}
}
} // if (currentdata.nexc[exclevel] > 0)
} // if (type_required == 0)
if (type_required == 1) {
// We move the inner(right)most particle of the highest excited level one more step towards
// the center of the block of Ix2 vacancies in which it sits.
// Is there already an excitation at level 0? If so, move inner(right)most particle one step further.
if (currentdata.nexc[exclevel] > 0) {
// Produce a new state_label_data:
State_Label_Data descdata = Read_State_Label (ScanStateLabel, OriginIx2);
// Identify the inner(right)most excitation:
int innerindex = currentdata.nexc[exclevel]/2;
int partpos = currentdata.Ix2exc[exclevel][innerindex];
// Determine the size of the block of OriginState.Ix2 vacancies in which this particle sits:
bool partpos_is_left_of_all_OriginIx2 = (partpos <= OriginIx2[descdata.type[exclevel] ].min());
bool partpos_is_right_of_all_OriginIx2 = (partpos >= OriginIx2[descdata.type[exclevel] ].max());
int nremptytoleft = 0;
if (partpos_is_left_of_all_OriginIx2) nremptytoleft = (partpos - Ix2_min[descdata.type[exclevel] ])/2;
else while (!OriginIx2[descdata.type[exclevel] ].includes (partpos -2*(nremptytoleft + 1)))
nremptytoleft++;
int nremptytoright = 0;
if (partpos_is_right_of_all_OriginIx2) nremptytoright = (Ix2_max[descdata.type[exclevel] ] - partpos)/2;
else while (!OriginIx2[descdata.type[exclevel] ].includes (partpos +2*(nremptytoright + 1)))
nremptytoright++;
if (!partpos_is_left_of_all_OriginIx2 && (partpos_is_right_of_all_OriginIx2
|| nremptytoleft - nremptytoright + 2 < 0))
partpos += 2;
else if (!partpos_is_right_of_all_OriginIx2 && (partpos_is_left_of_all_OriginIx2
|| nremptytoleft - nremptytoright - 2 >= 0))
partpos -= 2;
if (partpos != currentdata.Ix2exc[exclevel][innerindex] // we have successfully moved the particle
&& !OriginIx2[descdata.type[exclevel] ].includes(partpos) // it's actually a new particle position
&& !currentdata.Ix2exc[exclevel].includes(partpos) // new particle position is not already taken
&& partpos >= Ix2_min[descdata.type[exclevel] ]
&& partpos <= Ix2_max[descdata.type[exclevel] ]
)
{ // we have found a descendent
descdata.Ix2exc[exclevel][innerindex] = partpos;
desclabelfound[ndesc_found] = Return_State_Label (descdata, OriginIx2);
desctypefound[ndesc_found] = 1;
ndesc_found++;
if (cinbreaks) { char a; cin >> a;}
}
} // if (currentdata.nexc[exclevel] > 0)
} // if (type_required == 1)
if (type_required == 2) {
// Now add a new p-h pair at the inner(right)most position, at each level from exclevel upwards,
// putting the particle and hole to each available positions at the edge of vacancy/occupancy blocks.
for (int exclevel_newph = exclevel; exclevel_newph < currentdata.nexc.size(); ++exclevel_newph) {
if (ScanIx2[currentdata.type[exclevel_newph] ].size() <= currentdata.nexc[exclevel_newph])
continue; // no space for another p-h
// Define the objects for the newstatedata:
Vect<int> type_new = currentdata.type;
Vect<int> M_new = currentdata.M;
Vect<int> nexc_new = currentdata.nexc;
nexc_new[exclevel_newph] += 1; // we drill one more particle-hole pair at this level
int index_new = (currentdata.nexc[exclevel_newph] + 1)/2; // we put the new p-h pair at index index_new.
//int ntypespresent = ScanIx2.size();
int ntypespresent = currentdata.type.size();
Vect<Vect<int> > Ix2old_new(ntypespresent);
Vect<Vect<int> > Ix2exc_new(ntypespresent);
for (int it = 0; it < ntypespresent; ++it) Ix2old_new[it] = Vect<int>(ABACUS::max(nexc_new[it],1));
for (int it = 0; it < ntypespresent; ++it) Ix2exc_new[it] = Vect<int>(ABACUS::max(nexc_new[it],1));
// Copy earlier data in:
for (int it = 0; it < ntypespresent; ++it) {
for (int i = 0; i < currentdata.nexc[it]; ++i) {
Ix2old_new[it][i + (it == exclevel_newph && i >= index_new)] = currentdata.Ix2old[it][i];
Ix2exc_new[it][i + (it == exclevel_newph && i >= index_new)] = currentdata.Ix2exc[it][i];
}
}
State_Label_Data descdatanewph (type_new, M_new, nexc_new, Ix2old_new, Ix2exc_new);
// We now look for all possible hole positions,
// the allowable ones being either at the edge of a block, or next to an existing hole excitation,
// under the condition of obeying the `towards the block center' rule.
Vect<bool> isgoodnewholepos = Is_Good_New_Hole_Position (OriginIx2, currentdata, exclevel_newph);
// We now look for all possible particle positions,
// the allowable ones being either at the edge of a block, or next to an existing particle excitation,
// under the condition of obeying the `towards the block center' rule.
// Determine range of possible Ix2exc:
int Ix2excmin = ScanIx2[descdatanewph.type[exclevel_newph] ].min() - 2;
if (currentdata.nexc[exclevel_newph] > 0)
Ix2excmin = ABACUS::min (Ix2excmin, currentdata.Ix2exc[exclevel_newph][index_new - 1] + 2);
Ix2excmin = ABACUS::max (Ix2excmin, Ix2_min[descdatanewph.type[exclevel_newph] ]);
int Ix2excmax = ScanIx2[descdatanewph.type[exclevel_newph] ].max() + 2;
if (currentdata.nexc[exclevel_newph] > 1)
Ix2excmax = ABACUS::max (Ix2excmax, currentdata.Ix2exc[exclevel_newph][index_new] - 2);
Ix2excmax = ABACUS::min (Ix2excmax, Ix2_max[descdatanewph.type[exclevel_newph] ]);
Vect<bool> isgoodnewpartpos(false, (Ix2excmax - Ix2excmin)/2 + 1);
for (int candidateIx2exc = Ix2excmin; candidateIx2exc <= Ix2excmax; candidateIx2exc += 2) {
// candidateIx2exc is an acceptable position for the new particle provided the following conditions are fulfilled:
// A- it is not in OriginIx2
// B- it follows the ordering rule, i.e. it sits in the middle of previous particle excitations,
// namely between Ix2exc[index_new - 1] (if this exists) and Ix2exc[index_new] (if this exists)
// C- it is next to an OriginIx2 or immediately right of Ix2exc[index_new - 1] (if this exists)
// or immediately left of Ix2exc[index_new] (if this exists)
// D- it does not break the `towards the center' rule
// (it will break the rule at this point if it is created away from an OriginIx2 boundary
// (and thus next to a preexisting excitation),
// and if this excitation and it are not in the same sideblock
// (in other words: if there is a sideblock boundary between them)
// A
if (OriginIx2[descdatanewph.type[exclevel_newph] ].includes (candidateIx2exc)) {
// is contained in OriginIx2
continue;
}
// B1
if (currentdata.nexc[exclevel_newph] > 0
&& candidateIx2exc <= currentdata.Ix2exc[exclevel_newph][index_new - 1]) {
// there is at least one particle exc to the left, and the candidate position isn't right of Ix2exc[index_new - 1]
continue;
}
// B2
if (currentdata.nexc[exclevel_newph] > 1
&& candidateIx2exc >= currentdata.Ix2exc[exclevel_newph][index_new]) {
// there is at least one particle exc to the right, and the candidate position isn't left of Ix2exc[index_new]
continue;
}
// C- it is next to an OriginIx2 or immediately right of Ix2exc[index_new - 1] (if this exists)
// or immediately left of Ix2exc[index_new] (if this exists)
if (!(OriginIx2[descdatanewph.type[exclevel_newph] ].includes (candidateIx2exc + 2)
|| OriginIx2[descdatanewph.type[exclevel_newph] ].includes (candidateIx2exc - 2))
// doesn't sit next to an OriginIx2
&& (currentdata.nexc[exclevel_newph] == 0
|| candidateIx2exc != currentdata.Ix2exc[exclevel_newph][index_new - 1] + 2)
// doesn't sit immediately right of first particle excitation to the left
&& (currentdata.nexc[exclevel_newph] <= 1
|| candidateIx2exc != currentdata.Ix2exc[exclevel_newph][index_new] - 2)
// doesn't sit immediately left of first particle excitation to the right
) {
continue;
}
// D- it does not break the `towards the center' rule
// In other words, if created away from a block boundary but next to a preexisting particle,
// must be in same sideblock as this particle:
// Determine the size of the block of OriginIx2 vacancies in which this particle sits:
bool candidate_is_left_of_all_OriginIx2 =
(candidateIx2exc < OriginIx2[descdatanewph.type[exclevel_newph] ].min()); // this makes it left-moving
bool candidate_is_right_of_all_OriginIx2 =
(candidateIx2exc > OriginIx2[descdatanewph.type[exclevel_newph] ].max()); // this makes it right-moving
int nremptytoleft = 0;
if (candidate_is_left_of_all_OriginIx2)
nremptytoleft = (candidateIx2exc - Ix2_min[descdatanewph.type[exclevel_newph] ])/2;
else while (!OriginIx2[descdatanewph.type[exclevel_newph] ].includes (candidateIx2exc -2*(nremptytoleft + 1)))
nremptytoleft++;
int nremptytoright = 0;
if (candidate_is_right_of_all_OriginIx2)
nremptytoright = (Ix2_max[descdatanewph.type[exclevel_newph] ] - candidateIx2exc)/2;
else while (!OriginIx2[descdatanewph.type[exclevel_newph] ].includes (candidateIx2exc +2*(nremptytoright + 1)))
nremptytoright++;
// We can determine whether the new particle would be left- or right-moving
bool candidate_is_left_moving = candidate_is_left_of_all_OriginIx2
|| (!candidate_is_right_of_all_OriginIx2 && nremptytoleft >= nremptytoright);
bool candidate_is_right_moving = candidate_is_right_of_all_OriginIx2
|| (!candidate_is_left_of_all_OriginIx2 && nremptytoleft < nremptytoright);
// Consistency checks:
if (candidate_is_left_moving && candidate_is_right_moving)
ABACUSerror("New particle moving left and right at same time");
if (!candidate_is_left_moving && !candidate_is_right_moving)
ABACUSerror("New particle not moving either left or right");
if (candidate_is_left_moving
&& (currentdata.nexc[exclevel_newph] > 0
&& candidateIx2exc == currentdata.Ix2exc[exclevel_newph][index_new - 1] + 2)
// it is created to the right of a preexisting particle excitation
&& !OriginIx2[descdatanewph.type[exclevel_newph] ].includes (candidateIx2exc + 2)
// and is not sitting at the boundary
&& (currentdata.nexc[exclevel_newph] <= 1
|| candidateIx2exc != currentdata.Ix2exc[exclevel_newph][index_new] - 2)
// and is not sitting just left of another preexisting particle excitation which is also left moving
)
{
continue;
}
if (candidate_is_right_moving
&& (currentdata.nexc[exclevel_newph] > 1
&& candidateIx2exc == currentdata.Ix2exc[exclevel_newph][index_new] - 2)
// it is created to the left of a preexisting particle excitation
&& !OriginIx2[descdatanewph.type[exclevel_newph] ].includes (candidateIx2exc - 2)
// and is not sitting at the boundary
&& (currentdata.nexc[exclevel_newph] == 0
|| candidateIx2exc != currentdata.Ix2exc[exclevel_newph][index_new - 1] + 2)
// and is not sitting just right of another preexisting particle excitation which is also right moving
)
{
continue;
}
// If we have reached this point, candidateIx2exc is acceptable.
// Immediately construct all descendents with this part position:
// We now select the hole position closest to left or right:
int ihclosestleft = -1;
int ihclosestright = -1;
for (int ih = 0; ih < OriginIx2[descdatanewph.type[exclevel_newph] ].size(); ++ih)
if (isgoodnewholepos[ih])
{
if (OriginIx2[descdatanewph.type[exclevel_newph] ][ih] - candidateIx2exc < 0)
// new hole is left of new particle
if (ihclosestleft == -1 // we hadn't found a hole previously
|| ihclosestleft >= 0 // we had aleady found a new hole to the left
&& OriginIx2[descdatanewph.type[exclevel_newph] ][ih] >
OriginIx2[descdatanewph.type[exclevel_newph] ][ihclosestleft])
ihclosestleft = ih;
if (OriginIx2[descdatanewph.type[exclevel_newph] ][ih] - candidateIx2exc > 0)
// new hole is right of new particle
if (ihclosestright == -1 // we hadn't found a hole previously
|| ihclosestright >= 0 // we had aleady found a new hole to the right
&& OriginIx2[descdatanewph.type[exclevel_newph] ][ih] <
OriginIx2[descdatanewph.type[exclevel_newph] ][ihclosestright])
ihclosestright = ih;
} // if (isgoodnewholepos[ih])
// for ih
if (ihclosestleft >= 0) { // Found a descendent with new hole to left of new particle
descdatanewph.Ix2old[exclevel_newph][index_new] = OriginIx2[descdatanewph.type[exclevel_newph] ][ihclosestleft];
descdatanewph.Ix2exc[exclevel_newph][index_new] = candidateIx2exc;
desclabelfound[ndesc_found] = Return_State_Label (descdatanewph, OriginIx2);
desctypefound[ndesc_found] = 2;
ndesc_found++;
}
if (ihclosestright >= 0) { // Found a descendent with new hole to right of new particle
descdatanewph.Ix2old[exclevel_newph][index_new] = OriginIx2[descdatanewph.type[exclevel_newph] ][ihclosestright];
descdatanewph.Ix2exc[exclevel_newph][index_new] = candidateIx2exc;
desclabelfound[ndesc_found] = Return_State_Label (descdatanewph, OriginIx2);
desctypefound[ndesc_found] = 2;
ndesc_found++;
}
} // for (int candidateIx2exc)
} // for (int exclevel_newph)
} // if (type_required == 2)
if (type_required == 3) {
// We here only move the leftmost Ix2[0] left, and rightmost Ix2[0] right (this preserves iK)
if (ScanIx2[0][0] > Ix2_min[0] && ScanIx2[0][ScanIx2[0].size() - 1] < Ix2_max[0]) {
Vect<Vect<int> > descIx2 = ScanIx2;
descIx2[0][0] -= 2;
descIx2[0][descIx2[0].size() - 1] += 2;
desclabelfound[ndesc_found] = Return_State_Label (descIx2, OriginIx2);
desctypefound[ndesc_found] = 3;
ndesc_found++;
}
} // if (type_required == 3)
if (type_required == 4) {
// We take the latest particle-hole pair, and move the hole to a further generalized Umklapp position
State_Label_Data descdata = currentdata;
// Start by finding the possible hole positions of the state without the innermost p-h pair:
Vect<int> nexc_new = currentdata.nexc;
if (currentdata.nexc[exclevel] < 1)
ABACUSerror("Should not call type 3 descendent if there are no ph exc at exclevel");
nexc_new[exclevel] -= 1; // we remove the innermost ph pair
int innerindex = currentdata.nexc[exclevel]/2; // index of ph pair we remove
int ntypespresent = currentdata.type.size();
Vect<Vect<int> > Ix2old_new(ntypespresent);
Vect<Vect<int> > Ix2exc_new(ntypespresent);
for (int it = 0; it < ntypespresent; ++it) Ix2old_new[it] = Vect<int>(ABACUS::max(nexc_new[it],1));
for (int it = 0; it < ntypespresent; ++it) Ix2exc_new[it] = Vect<int>(ABACUS::max(nexc_new[it],1));
// Copy earlier data in:
for (int it = 0; it < ntypespresent; ++it) {
for (int i = 0; i < nexc_new[it]; ++i) {
Ix2old_new[it][i] = currentdata.Ix2old[it][i + (it == exclevel && i >= innerindex)];
Ix2exc_new[it][i] = currentdata.Ix2exc[it][i + (it == exclevel && i >= innerindex)];
}
}
State_Label_Data data_latest_ph_removed (currentdata.type, currentdata.M, nexc_new, Ix2old_new, Ix2exc_new);
Vect<bool> isgoodnewholepos = Is_Good_New_Hole_Position (OriginIx2, data_latest_ph_removed, exclevel);
// Move hole further away from its corresponding particle, preserving the direction
bool holemovingleft = currentdata.Ix2old[exclevel][innerindex] < currentdata.Ix2exc[exclevel][innerindex];
int inewhole;
if (holemovingleft) {
inewhole = 0;
while (OriginIx2[currentdata.type[exclevel] ][inewhole] != currentdata.Ix2old[exclevel][innerindex]) inewhole++;
do {
inewhole--;
} while (inewhole >= 0 && !isgoodnewholepos[inewhole]);
}
if (!holemovingleft) {
inewhole = OriginIx2[currentdata.type[exclevel] ].size() - 1;
while (OriginIx2[currentdata.type[exclevel] ][inewhole] != currentdata.Ix2old[exclevel][innerindex]) inewhole--;
do {
inewhole++;
} while (inewhole <= OriginIx2[currentdata.type[exclevel] ].size() - 1 && !isgoodnewholepos[inewhole]);
}
if (inewhole >= 0 && inewhole <= OriginIx2[currentdata.type[exclevel] ].size() - 1) { // Found a descendent
descdata.Ix2old[exclevel][innerindex] = OriginIx2[currentdata.type[exclevel] ][inewhole];
desclabelfound[ndesc_found] = Return_State_Label (descdata, OriginIx2);
desctypefound[ndesc_found] = 4;
ndesc_found++;
}
} // if (type_required == 4)
Vect<string> desclabelfoundresized (ndesc_found);
for (int i = 0; i < ndesc_found; ++i) {
desclabelfoundresized[i] = desclabelfound[i];
}
if (cinbreaks) { char a; cin >> a;}
return(desclabelfoundresized);
} // Vect<string> Descendents
// Specialization for LiebLin gas:
Vect<string> Descendents (const LiebLin_Bethe_State& ScanState, const LiebLin_Bethe_State& OriginState, int type_required)
{
Vect<Vect<int> > ScanIx2here(1);
ScanIx2here[0] = ScanState.Ix2;
Vect<Vect<int> > OriginIx2here(1);
OriginIx2here[0] = OriginState.Ix2;
Vect<int> Ix2_min(1);
Ix2_min[0] = LIEBLIN_Ix2_MIN;
Vect<int> Ix2_max(1);
Ix2_max[0] = LIEBLIN_Ix2_MAX;
Vect<string> desc_found;
// If the state is an outer skeleton, we use a slightly modified OriginIx2 (in which the outermost
// rapidities at level zero are put to the outer skeleton position convention) for computing the descendents.
if (Is_Outer_Skeleton(ScanState)) {
// Construct and descend a state with 2 less particles:
if (ScanState.N < 3) ABACUSerror("Skeleton descendent logic at fixed iK not implemented for N < 3.");
LiebLin_Bethe_State ReducedScanState (ScanState.c_int, ScanState.L, ScanState.N - 2);
LiebLin_Bethe_State ReducedOriginState (ScanState.c_int, ScanState.L, ScanState.N - 2);
for (int i = 0; i < ScanState.N - 2; ++i) {
ReducedScanState.Ix2[i] = ScanState.Ix2[i+1];
ReducedOriginState.Ix2[i] = OriginState.Ix2[i+1];
}
ReducedScanState.Set_Label_from_Ix2 (ReducedOriginState.Ix2);
Vect<string> desc_mod = Descendents (ReducedScanState, ReducedOriginState, type_required);
// Now translate results back to ScanState: we return outer skeleton states
LiebLin_Bethe_State LabellingState = ScanState;
for (int idesc = 0; idesc < desc_mod.size(); ++idesc) {
ReducedScanState.Set_to_Label (desc_mod[idesc], ReducedOriginState.Ix2);
// Now set all quantum numbers of LabellingState:
LabellingState.Ix2[0] = LIEBLIN_Ix2_MIN + (ScanState.N % 2) + 1;
for (int i = 0; i < ScanState.N - 2; ++i) LabellingState.Ix2[i+1] = ReducedScanState.Ix2[i];
LabellingState.Ix2[ScanState.N - 1] = LIEBLIN_Ix2_MAX - (ScanState.N % 2) - 1;
LabellingState.Set_Label_from_Ix2 (OriginState.Ix2); // Now reset label on OriginIx2 from the correctly set Ix2
desc_mod[idesc] = LabellingState.label;
}
desc_found = desc_mod;
}
// If not outer skeleton, just return straight descendents
else desc_found = Descendents (ScanState.label, ScanIx2here, OriginIx2here, Ix2_min, Ix2_max, type_required);
return(desc_found);
}
// Specialization for Heisenberg:
Vect<string> Descendents (const Heis_Bethe_State& ScanState, const Heis_Bethe_State& OriginState, int type_required)
{
return(Descendents (ScanState.label, ScanState.Ix2, OriginState.Ix2,
ScanState.base.Ix2_min, ScanState.base.Ix2_max, type_required));
}
} // namespace ABACUS
Reference in New Issue View Git Blame Copy Permalink