262 lines
11 KiB
C++
262 lines
11 KiB
C++
/**********************************************************
|
|
|
|
This software is part of J.-S. Caux's ABACUS library.
|
|
|
|
Copyright (c) J.-S. Caux.
|
|
|
|
-----------------------------------------------------------
|
|
|
|
File: NRG_DFF_Matrix_Block_builder.cc
|
|
|
|
Purpose: calculate matrix elements of density operator between
|
|
selected states in NRG for a selected block of the whole matrix.
|
|
For collaboration with Robert Konik.
|
|
|
|
***********************************************************/
|
|
|
|
#include "ABACUS.h"
|
|
|
|
using namespace std;
|
|
using namespace ABACUS;
|
|
|
|
namespace ABACUS {
|
|
|
|
void Build_DME_Matrix_Block_for_NRG (DP c_int, DP L, int N, int iKmin, int iKmax,
|
|
int Nstates_required, bool symmetric_states, int iKmod,
|
|
int weighing_option, int nrglabel_left_begin, int nrglabel_left_end,
|
|
int nrglabel_right_begin, int nrglabel_right_end,
|
|
int block_option, DP* DME_block_1, DP* DME_block_2, Vect_DP Kweight)
|
|
{
|
|
// Given a list of states produced by Select_States_for_NRG, this function
|
|
// computes all the density form factors < lstate| \rho | rstate >
|
|
// for all id_left_begin <= id_left <= id_left_end (similarly for right states).
|
|
|
|
// The block_option flag determines what kind of data is written in the 2 DME blocks.
|
|
// If !symmetric_states, then only block 1 is filled.
|
|
// If symmetric_states, then
|
|
// if block_option == 1, only block 1 is filled with the symmetric states ME.
|
|
// if block_option == 2, both blocks 1 and 2 are filled with ME and ME (after parity) respectively.
|
|
|
|
// ASSUMPTIONS:
|
|
// We assume the DME_blocks are already reserved in memory.
|
|
// If !symmetric states, DME_block_2 doesn't need to be allocated.
|
|
|
|
if (nrglabel_left_begin < 0 || nrglabel_right_begin < 0)
|
|
ABACUSerror("beginning nrglabels negative in Build_DME_Matrix_Block_for_NRG");
|
|
if (nrglabel_left_end >= Nstates_required)
|
|
ABACUSerror("nrglabel_left_end too large in Build_DME_Matric_Block_for_NRG");
|
|
if (nrglabel_right_end >= Nstates_required)
|
|
ABACUSerror("nrglabel_right_end too large in Build_DME_Matric_Block_for_NRG");
|
|
|
|
if (nrglabel_left_begin > nrglabel_left_end)
|
|
ABACUSerror("nrglabels of left states improperly chosen in DME block builder.");
|
|
if (nrglabel_right_begin > nrglabel_right_end)
|
|
ABACUSerror("nrglabels of right states improperly chosen in DME block builder.");
|
|
|
|
// DME_block is a pointer of size row_l * col_l, where
|
|
int block_row_length = nrglabel_right_end - nrglabel_right_begin + 1;
|
|
int block_col_length = nrglabel_left_end - nrglabel_left_begin + 1;
|
|
|
|
// iKmod is an integer specifying that we're only interested in momenta multiples of
|
|
// a base unit. In other words, if the perturbation potential repeats itself iKmod times
|
|
// on the periodic cycle on which the theory is defined, then we only need states such
|
|
// that iKl - iKr = (integer) * iKmod. We impose this constraint at the same time as
|
|
// the check on the weight of the state.
|
|
|
|
// We start by reading off the list of selected states:
|
|
|
|
stringstream NRG_stringstream;
|
|
string NRG_string;
|
|
|
|
NRG_stringstream << "States_";
|
|
NRG_stringstream << "c_" << c_int << "_L_" << L << "_N_" << N
|
|
<< "_iKmin_" << iKmin << "_iKmax_" << iKmax << "_Nstates_" << Nstates_required
|
|
<< "_Sym_" << symmetric_states << "_iKmod_" << iKmod << "_wopt_" << weighing_option << ".nrg";
|
|
|
|
NRG_string = NRG_stringstream.str();
|
|
const char* NRG_Cstr = NRG_string.c_str();
|
|
|
|
ifstream infile;
|
|
infile.open(NRG_Cstr);
|
|
|
|
if (infile.fail()) {
|
|
cout << NRG_Cstr << endl;
|
|
ABACUSerror("The input file was not opened successfully in Build_DME_Matrix_between_Selected_States_for_NRG. ");
|
|
}
|
|
|
|
// Read the whole data file:
|
|
const int MAXDATA = ABACUS::max(nrglabel_left_end, nrglabel_right_end) + 10; // 10 for safety...
|
|
|
|
int* nrglabel = new int[MAXDATA];
|
|
DP* omega = new DP[MAXDATA];
|
|
int* iK = new int[MAXDATA];
|
|
string* label = new string[MAXDATA];
|
|
bool* sym = new bool[MAXDATA];
|
|
DP* weight = new DP[MAXDATA];
|
|
|
|
int Ndata = 0;
|
|
while (((infile.peek()) != EOF) && (Ndata < MAXDATA)) {
|
|
|
|
infile >> nrglabel[Ndata];
|
|
if (nrglabel[Ndata] != Ndata) ABACUSerror("reading states nrglabels wrong in NRG_DME_Matrix_Block_builder");
|
|
infile >> omega[Ndata];
|
|
infile >> iK[Ndata];
|
|
infile >> label[Ndata];
|
|
infile >> sym[Ndata];
|
|
infile >> weight[Ndata];
|
|
|
|
Ndata++;
|
|
}
|
|
infile.close();
|
|
|
|
// We first reconstruct all the needed eigenstates:
|
|
|
|
// Vector containing all the kept states:
|
|
Vect<LiebLin_Bethe_State> Kept_States_left(block_col_length);
|
|
Vect<LiebLin_Bethe_State> Kept_States_right(block_row_length);
|
|
|
|
// Dummy state for calculations
|
|
LiebLin_Bethe_State Scanstate (c_int, L, N);
|
|
|
|
// State list containing all the types of states
|
|
Scan_State_List<LiebLin_Bethe_State> List_of_Basic_States ('Z', Scanstate);
|
|
|
|
for (int nrglabel_left = nrglabel_left_begin; nrglabel_left <= nrglabel_left_end; ++nrglabel_left) {
|
|
|
|
Scanstate = List_of_Basic_States.Return_State (Extract_Base_Label(label[nrglabel_left]));
|
|
Scanstate.Set_to_Label (label[nrglabel_left]);
|
|
|
|
Scanstate.Compute_All(true);
|
|
|
|
Kept_States_left[nrglabel_left - nrglabel_left_begin] = Scanstate;
|
|
|
|
if (!Scanstate.conv)
|
|
cout << "State of label " << label[nrglabel_left] << " did not converge after " << Scanstate.iter_Newton
|
|
<< " Newton step; diffsq = " << Scanstate.diffsq << endl;
|
|
} // for nrglabel_left
|
|
|
|
for (int nrglabel_right = nrglabel_right_begin; nrglabel_right <= nrglabel_right_end; ++nrglabel_right) {
|
|
|
|
Scanstate = List_of_Basic_States.Return_State (Extract_Base_Label(label[nrglabel_right]));
|
|
Scanstate.Set_to_Label (label[nrglabel_right]);
|
|
|
|
Scanstate.Compute_All(true);
|
|
|
|
Kept_States_right[nrglabel_right - nrglabel_right_begin] = Scanstate;
|
|
|
|
if (!Scanstate.conv)
|
|
cout << "State of label " << label[nrglabel_right] << " did not converge after " << Scanstate.iter_Newton
|
|
<< " Newton step; diffsq = " << Scanstate.diffsq << endl;
|
|
} // for nrglabel_left
|
|
|
|
|
|
// Now that we have all the states, we can do the full calculation of all cross form factors:
|
|
|
|
for (int ll = 0; ll < block_col_length; ++ll) {
|
|
for (int lr = 0; lr < block_row_length; ++lr) {
|
|
|
|
if (Kept_States_left[ll].conv && Kept_States_right[lr].conv
|
|
&& abs(Kept_States_left[ll].iK) % iKmod == 0
|
|
&& abs(Kept_States_right[lr].iK) % iKmod == 0) {
|
|
|
|
if (!symmetric_states) {
|
|
DME_block_1[ll* block_row_length + lr] = // By convention, put the lowest energy state to the left
|
|
Kweight[abs(Kept_States_left[ll].iK - Kept_States_right[lr].iK)]
|
|
* (Kept_States_left[ll].E <= Kept_States_right[lr].E
|
|
? real(exp(ln_Density_ME (Kept_States_left[ll], Kept_States_right[lr])))
|
|
: real(exp(ln_Density_ME (Kept_States_right[lr], Kept_States_left[ll]))));
|
|
// We don't do anything with block 2, we don't even assume it's been allocated.
|
|
}
|
|
|
|
else if (symmetric_states) {
|
|
// Check parity of both states:
|
|
bool Lstate_parity_inv = Kept_States_left[ll].Check_Symmetry();
|
|
bool Rstate_parity_inv = Kept_States_right[lr].Check_Symmetry();
|
|
|
|
if (Lstate_parity_inv && Rstate_parity_inv) {
|
|
DME_block_1[ll* block_row_length + lr] = Kweight[abs(Kept_States_left[ll].iK - Kept_States_right[lr].iK)]
|
|
* real(exp(ln_Density_ME (Kept_States_left[ll], Kept_States_right[lr])));
|
|
if (block_option == 2) DME_block_2[ll* block_row_length + lr] = 0.0;
|
|
}
|
|
else if (!Lstate_parity_inv && !Rstate_parity_inv) {
|
|
LiebLin_Bethe_State PRstate = Kept_States_right[lr];
|
|
PRstate.Parity_Flip();
|
|
|
|
if (block_option == 1) {
|
|
DME_block_1[ll* block_row_length + lr] =
|
|
(Kweight[abs(Kept_States_left[ll].iK - Kept_States_right[lr].iK)]
|
|
* real(exp(ln_Density_ME (Kept_States_left[ll], Kept_States_right[lr])))
|
|
+ Kweight[abs(Kept_States_left[ll].iK - PRstate.iK)]
|
|
* real(exp(ln_Density_ME (Kept_States_left[ll], PRstate))));
|
|
// We don't do anything with block 2, we don't even assume it's been allocated.
|
|
}
|
|
else if (block_option == 2) {
|
|
DME_block_1[ll* block_row_length + lr] =
|
|
Kweight[abs(Kept_States_left[ll].iK - Kept_States_right[lr].iK)]
|
|
* real(exp(ln_Density_ME (Kept_States_left[ll], Kept_States_right[lr])));
|
|
DME_block_2[ll* block_row_length + lr] =
|
|
Kweight[abs(Kept_States_left[ll].iK - PRstate.iK)]
|
|
* real(exp(ln_Density_ME (Kept_States_left[ll], PRstate)));
|
|
}
|
|
}
|
|
|
|
else if (Lstate_parity_inv) {
|
|
LiebLin_Bethe_State PRstate = Kept_States_right[lr];
|
|
PRstate.Parity_Flip();
|
|
|
|
if (block_option == 1) {
|
|
DME_block_1[ll* block_row_length + lr] = sqrt(0.5) *
|
|
(Kweight[abs(Kept_States_left[ll].iK - Kept_States_right[lr].iK)]
|
|
* real(exp(ln_Density_ME (Kept_States_left[ll], Kept_States_right[lr])))
|
|
+ Kweight[abs(Kept_States_left[ll].iK - PRstate.iK)]
|
|
* real(exp(ln_Density_ME (Kept_States_left[ll], PRstate))));
|
|
// We don't do anything with block 2, we don't even assume it's been allocated.
|
|
}
|
|
else if (block_option == 2) {
|
|
DME_block_1[ll* block_row_length + lr] = sqrt(0.5) *
|
|
Kweight[abs(Kept_States_left[ll].iK - Kept_States_right[lr].iK)]
|
|
* real(exp(ln_Density_ME (Kept_States_left[ll], Kept_States_right[lr])));
|
|
DME_block_2[ll* block_row_length + lr] = sqrt(0.5) *
|
|
Kweight[abs(Kept_States_left[ll].iK - PRstate.iK)]
|
|
* real(exp(ln_Density_ME (Kept_States_left[ll], PRstate)));
|
|
}
|
|
}
|
|
|
|
else { // only R state is parity invariant, flip left:
|
|
LiebLin_Bethe_State PLstate = Kept_States_left[ll];
|
|
PLstate.Parity_Flip();
|
|
|
|
if (block_option == 1) {
|
|
DME_block_1[ll* block_row_length + lr] = sqrt(0.5) *
|
|
(Kweight[abs(Kept_States_left[ll].iK - Kept_States_right[lr].iK)]
|
|
* real(exp(ln_Density_ME (Kept_States_left[ll], Kept_States_right[lr])))
|
|
+ Kweight[abs(PLstate.iK - Kept_States_right[lr].iK)]
|
|
* real(exp(ln_Density_ME (PLstate, Kept_States_right[lr]))));
|
|
// We don't do anything with block 2, we don't even assume it's been allocated.
|
|
}
|
|
else if (block_option == 2) {
|
|
DME_block_1[ll* block_row_length + lr] = sqrt(0.5) *
|
|
Kweight[abs(Kept_States_left[ll].iK - Kept_States_right[lr].iK)]
|
|
* real(exp(ln_Density_ME (Kept_States_left[ll], Kept_States_right[lr])));
|
|
DME_block_2[ll* block_row_length + lr] = sqrt(0.5) *
|
|
Kweight[abs(PLstate.iK - Kept_States_right[lr].iK)]
|
|
* real(exp(ln_Density_ME (PLstate, Kept_States_right[lr])));
|
|
}
|
|
}
|
|
|
|
} // if (symmetric_states)
|
|
} // if conv & iKmod condition fulfilled
|
|
|
|
else {
|
|
DME_block_1[ll* block_row_length + lr] = 0.0;
|
|
// condition, to prevent segfault if DME_block_2 not allocated.
|
|
if (symmetric_states && block_option == 2) DME_block_2[ll* block_row_length + lr] = 0.0;
|
|
}
|
|
} // for lr
|
|
} // for ll
|
|
|
|
return;
|
|
}
|
|
|
|
} // namespace ABACUS
|