ABACUS/src/LIEBLIN/LiebLin_Matrix_Element_Contrib.cc

/**********************************************************

This software is part of J.-S. Caux's ABACUS library.

Copyright (c) J.-S. Caux.

-----------------------------------------------------------

File:  src/LIEBLIN/LiebLin_Matrix_Element_Contrib.cc

Purpose:  handles the generic call for a Lieb-Liniger matrix element contribution.

***********************************************************/

#include "ABACUS.h"

using namespace std;
using namespace ABACUS;

namespace ABACUS {

  DP Compute_Matrix_Element_Contrib (char whichDSF, int iKmin, int iKmax, LiebLin_Bethe_State& LeftState,
				     LiebLin_Bethe_State& RefState, DP Chem_Pot, stringstream& DAT_outfile)
  {
    // This function prints the matrix element information to the fstream,
    // and returns a weighed `data_value' to be multiplied by sumrule_factor,
    // to determine if scanning along this thread should be continued.

    bool fixed_iK = (iKmin == iKmax);

    // Identify which matrix element is needed from the number of particles in Left and Right states:
    DP ME = 0.0;
    complex<DP> ME_CX = 0.0;

    if (whichDSF == 'Z')
      ME = LeftState.E - RefState.E;
    else if (whichDSF == 'd' || whichDSF == '1')
      ME = real(exp(ln_Density_ME (LeftState, RefState)));
    else if (whichDSF == 'o')
      ME = real(exp(ln_Psi_ME (LeftState, RefState)));
    else if (whichDSF == 'g')
      ME = real(exp(ln_Psi_ME (RefState, LeftState)));
    else if (whichDSF == 'q') // geometrical quench
      ME_CX = LiebLin_ln_Overlap(LeftState.lambdaoc, LeftState.lnnorm, RefState);
    else if (whichDSF == 'B') { // BEC to finite c quench: g2(x=0)
      // overlap part, relative to saddle-point state
      ME_CX = exp(ln_Overlap_with_BEC (LeftState) - ln_Overlap_with_BEC (RefState));
      // matrix element part
      ME = real(exp(ln_g2_ME (RefState, LeftState)));
      // The product is ME_CX * ME = e^{-\delta S_e} \langle \rho_{sp} | g2 (x=0) | \rho_{sp} + e \rangle
      // and is the first half of the t-dep expectation value formula coming from the Quench Action formalism
    }
    else if (whichDSF == 'C') { // BEC to finite c quench: overlap
      ME_CX = exp(2.0* ln_Overlap_with_BEC (LeftState)); // overlap sq part
      ME = 0.0;
    }
    else ABACUSerror("Wrong whichDSF in Compute_Matrix_Element_Contrib.");

    if (is_nan(ME)) ME = (whichDSF == 'Z') ?  1.0e+200 : 0.0;
    if (is_nan(norm(ME_CX))) ME_CX = -100.0;


    DAT_outfile << setprecision(16);
    // Print information to fstream:
    if (LeftState.iK - RefState.iK >= iKmin && LeftState.iK - RefState.iK <= iKmax) {
      if (whichDSF == 'Z') {
	DAT_outfile << endl << LeftState.E - RefState.E - (LeftState.N - RefState.N) * Chem_Pot << "\t"
		    << LeftState.iK - RefState.iK
		    << 0 << "\t" // This is the deviation, here always 0
		    << "\t" << LeftState.label;
      }
      else if (whichDSF == 'q') {
	DAT_outfile << endl << LeftState.E - RefState.E - (LeftState.N - RefState.N) * Chem_Pot << "\t"
		    << LeftState.iK - RefState.iK << "\t"
		    << real(ME_CX) << "\t" << imag(ME_CX) - twoPI * int(imag(ME_CX)/twoPI + 1.0e-10) << "\t"
		    << 0 << "\t" // This is the deviation, here always 0
		    << LeftState.label;
      }
      else if (whichDSF == '1') {
	DAT_outfile << endl << LeftState.E - RefState.E - (LeftState.N - RefState.N) * Chem_Pot << "\t"
		    << LeftState.iK - RefState.iK << "\t"
		    << ME << "\t"
		    << 0 << "\t" // This is the deviation, here always 0
		    << LeftState.label;
	DAT_outfile << "\t" << Momentum_Right_Excitations(LeftState) << "\t" << Momentum_Left_Excitations(LeftState);
      }
      else if (whichDSF == 'B') { // BEC to finite c > 0 quench; g2 (x=0)
	if (fabs(real(ME_CX) * ME) > 1.0e-100)
	  DAT_outfile << endl << LeftState.E - RefState.E << "\t"
		      << LeftState.iK - RefState.iK << "\t"
		      << real(ME_CX) << "\t" // the overlap is always real
		      << ME << "\t"
		      << LeftState.label;
      }
      else if (whichDSF == 'C') { // BEC to finite c, overlap sq
	if (fabs(real(ME_CX)) > 1.0e-100)
	  DAT_outfile << endl << LeftState.E - RefState.E << "\t"
		      << LeftState.iK - RefState.iK << "\t"
		      << real(ME_CX) << "\t" // the overlap is always real
		      << ME << "\t"
		      << LeftState.label;
      }
      else {
	DAT_outfile << endl << LeftState.E - RefState.E - (LeftState.N - RefState.N) * Chem_Pot << "\t"
		    << LeftState.iK - RefState.iK << "\t"
		    << ME << "\t"
		    << 0 << "\t" // This is the deviation, here always 0
		    << LeftState.label;
      }
    } // if iKmin <= iK <= iKmax

    // Calculate and return the data_value:
    DP data_value = ME * ME;
    if (whichDSF == 'Z') // use 1/(1 + omega)
      data_value = 1.0/(1.0 + LeftState.E - RefState.E - (LeftState.N - RefState.N) * Chem_Pot);
    else if (whichDSF == 'd' || whichDSF == '1') {
      if (fixed_iK)
	data_value = (LeftState.E - RefState.E - (LeftState.N - RefState.N) * Chem_Pot) * ME * ME;
      else if (!fixed_iK)
	data_value = (LeftState.E - RefState.E - (LeftState.N - RefState.N) * Chem_Pot) * ME * ME
	  /ABACUS::max(1, (LeftState.iK - RefState.iK) * (LeftState.iK - RefState.iK));
    }
    else if (whichDSF == 'g' || whichDSF == 'o') {
      if (fixed_iK)
	data_value = (LeftState.E - RefState.E - (LeftState.N - RefState.N) * Chem_Pot) * ME * ME;
      else if (!fixed_iK) {
	data_value = ME * ME;
      }
    }
    else if (whichDSF == 'q')
      data_value = exp(2.0 * real(ME_CX));
    else if (whichDSF == 'B')
      data_value = abs(ME_CX * ME)/(1.0 + sqrt(fabs(LeftState.E - RefState.E)));
    else if (whichDSF == 'C')
      data_value = abs(ME_CX);

    if (whichDSF == '1') // hard cutoff for nr, nl when calculating exponents
      if (Momentum_Right_Excitations(LeftState) >= 30 || Momentum_Left_Excitations(LeftState) >= 30)
	data_value = 0.0;

    return(data_value);
  }

} // namespace ABACUS