Jean-Sébastien Caux
197 lines
6.8 KiB
C++
197 lines
6.8 KiB
C++
/*! Type II quasisolitons in Lieb-Liniger: amplitudes_IDeltax
|
|
|
|
Purpose:
|
|
Produce the .amplitudes file
|
|
containing the (complex) amplitudes of each basis state
|
|
for amplitudes optimizing depletion density integral I(Δx)
|
|
for (multi-)type II hole wavepackets in Lieb-Liniger.
|
|
|
|
This executable requires Abacus version 2.
|
|
|
|
See README for compilation instructions.
|
|
|
|
Copyright © Jean-Sébastien Caux, Anahita Sarvi and Cesare Vianello.
|
|
|
|
This program is free software: you can redistribute it and/or modify it under the terms of the GNU Affero General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.
|
|
|
|
This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Affero General Public License for more details.
|
|
|
|
You should have received a copy of the GNU Affero General Public License along with this program. If not, see <https://www.gnu.org/licenses/>.
|
|
*/
|
|
|
|
|
|
import std;
|
|
|
|
import abacus;
|
|
|
|
|
|
int main(int argc, char* argv[])
|
|
{
|
|
using namespace std::complex_literals;
|
|
|
|
if (argc != 8) {
|
|
std::cout << "Executable amplitudes_IDeltax\n"
|
|
<< " part of the Type II quasisolitons in Lieb-Liniger software suite\n"
|
|
<< " copyright © Jean-Sébastien Caux, Anahita Sarvi and Cesare Vianello.\n";
|
|
std::cout << "\nPurpose:\n"
|
|
<< " Produce the .amplitudes file\n"
|
|
<< " containing the (complex) amplitudes of each basis state\n"
|
|
<< " for amplitudes optimizing depletion density integral I(Δx)\n"
|
|
<< " for (multi-)type II hole wavepackets in Lieb-Liniger.\n";
|
|
|
|
std::cout << "\nPrerequirements:\n----------------\n"
|
|
<< " - .states and .rho files produced by executable rho\n";
|
|
|
|
std::cout << "\nUsage:\n------\n";
|
|
std::cout << "holeIx2 <c> <L> <N> <nr holes> <width> <offset>\n\n";
|
|
|
|
int warg { 16 }, wtype { 10 }, wcons { 26 };
|
|
std::cout << std::left << std::setw(warg) << "Argument" << std::setw(wtype) << "Type"
|
|
<< std::setw(wcons) << "Constraints" << "Description\n";
|
|
std::cout << std::left << std::setw(warg) << "--------" << std::setw(wtype) << "----"
|
|
<< std::setw(wcons) << "-----------" << "-----------\n";
|
|
std::cout << std::left << std::setw(warg) << "c" << std::setw(wtype) << "Real"
|
|
<< std::setw(wcons) << "> 0"
|
|
<< "Value of the interaction parameter\n";
|
|
std::cout << std::left << std::setw(warg) << "L" << std::setw(wtype) << "Real"
|
|
<< std::setw(wcons) << "> 0"
|
|
<< "System size\n";
|
|
std::cout << std::left << std::setw(warg) << "N" << std::setw(wtype) << "int"
|
|
<< std::setw(wcons) << "> 0"
|
|
<< "Number of particles\n";
|
|
std::cout << std::left << std::setw(warg) << "nr holes" << std::setw(wtype) << "int"
|
|
<< std::setw(wcons) << "1 <= nr holes <= N"
|
|
<< "Number of holes (Type II modes)\n";
|
|
std::cout << std::left << std::setw(warg) << "width" << std::setw(wtype) << "int"
|
|
<< std::setw(wcons) << "nr holes < width <= N"
|
|
<< "Width of the hole window\n";
|
|
std::cout << std::left << std::setw(warg) << "offset" << std::setw(wtype) << "int"
|
|
<< std::setw(wcons) << "0 <= offset <= N-width"
|
|
<< "Offset of the hole window w/r to the right Fermi edge\n";
|
|
std::cout << std::left << std::setw(warg) << "Δx" << std::setw(wtype) << "Real"
|
|
<< std::setw(wcons) << "0 < Δx <= L"
|
|
<< "Spatial width of depletion density to optimize\n";
|
|
|
|
return 0;
|
|
}
|
|
|
|
std::cout << std::setprecision(std::numeric_limits<Real>::digits10 + 1);
|
|
|
|
Real c { std::stold(argv[1]) };
|
|
Real L { std::stold(argv[2]) };
|
|
int N { std::stoi(argv[3]) };
|
|
int nholes { std::stoi(argv[4]) };
|
|
int width { std::stoi(argv[5]) };
|
|
int offset { std::stoi(argv[6]) };
|
|
Real Δx {std::stold(argv[7]) };
|
|
|
|
|
|
// Define the input files
|
|
|
|
std::stringstream filename_base;
|
|
filename_base << "c_" << c << "_N_" << N << "_L_" << L << "_nholes_" << nholes;
|
|
filename_base << "_width_" << width;
|
|
filename_base << "_offset_" << offset;
|
|
|
|
std::stringstream states_filename;
|
|
states_filename << filename_base.str() << ".states";
|
|
std::ifstream states_file;
|
|
states_file.open(states_filename.str());
|
|
states_file >> std::setprecision(std::numeric_limits<Real>::digits10 + 1);
|
|
|
|
std::stringstream rho_ME_filename;
|
|
rho_ME_filename << filename_base.str() << ".rho";
|
|
std::ifstream rho_ME_file;
|
|
rho_ME_file.open(rho_ME_filename.str());
|
|
rho_ME_file >> std::setprecision(std::numeric_limits<Real>::digits10 + 1);
|
|
|
|
|
|
// Input the states info
|
|
|
|
std::vector<std::string> label;
|
|
std::vector<int> iK;
|
|
std::vector<Real> E;
|
|
std::string tmp_label;
|
|
int tmp_iK;
|
|
Real tmp_E;
|
|
states_file >> tmp_label;
|
|
do {
|
|
states_file >> tmp_iK >> tmp_E;
|
|
label.push_back(tmp_label);
|
|
iK.push_back(tmp_iK);
|
|
E.push_back(tmp_E);
|
|
} while (states_file >> tmp_label);
|
|
states_file.close();
|
|
|
|
int nr_states { int(label.size()) };
|
|
|
|
|
|
// Input the density operator matrix elements
|
|
|
|
std::map<std::string, std::map<std::string, std::complex<Real>>> rho_ME;
|
|
|
|
for (int ibra { 0 }; ibra < nr_states; ++ibra) {
|
|
for (int iket { 0 }; iket <= ibra; ++iket) {
|
|
rho_ME_file >> rho_ME[label[ibra]][label[iket]];
|
|
}
|
|
}
|
|
rho_ME_file.close();
|
|
|
|
|
|
|
|
// Define the M matrix for the given Δx
|
|
|
|
Matrix<Real> M (nr_states);
|
|
Real K;
|
|
M.setZero();
|
|
|
|
for (int ibra { 0 }; ibra < nr_states; ++ibra) {
|
|
for (int iket { 0 }; iket < ibra; ++iket) {
|
|
K = pi_r * (iK[ibra] - iK[iket]) * Δx/L; // K_{ab}*Δx/2
|
|
// To balance the matrix numerically, it is multiplied by sqrt(N)
|
|
M(ibra, iket) = std::sqrt(N) *
|
|
std::real(rho_ME[label[ibra]][label[iket]]) * (K != Real(0) ? std::sin(K)/K : 1);
|
|
M(iket, ibra) = M(ibra, iket); // matrix is symmetric
|
|
}
|
|
}
|
|
|
|
|
|
// Householder reduction and QL factorization
|
|
|
|
std::vector<Real> d(nr_states);
|
|
std::vector<Real> e(nr_states);
|
|
tred2(M, d, e);
|
|
tqli(d, e, M);
|
|
// std::cout << "Eigenvalues according to QL: " << d << std::endl;
|
|
|
|
// Find minimal eigenvalue; the optimal wavefunction is then the corresponding eigenvector
|
|
int min_red_ev_index { 0 };
|
|
Real min_red_ev { d[0] };
|
|
for (int ibra { 0 }; ibra < nr_states; ++ibra) {
|
|
if (d[ibra] < min_red_ev) {
|
|
min_red_ev_index = ibra;
|
|
min_red_ev = d[ibra];
|
|
}
|
|
}
|
|
|
|
|
|
// Define the output files
|
|
|
|
std::stringstream amplitudes_filename;
|
|
amplitudes_filename << filename_base.str() << "_IDeltax_" << Δx << ".amplitudes";
|
|
std::ofstream amplitudes_file;
|
|
amplitudes_file.open(amplitudes_filename.str(), std::ios::out | std::ios::trunc);
|
|
amplitudes_file << std::setprecision(std::numeric_limits<Real>::digits10 + 1);
|
|
|
|
|
|
// Output the state amplitudes
|
|
|
|
for (int ibra { 0 }; ibra < nr_states; ++ibra) {
|
|
amplitudes_file << M(ibra, min_red_ev_index) << "\t";
|
|
}
|
|
amplitudes_file.close();
|
|
|
|
|
|
return 0;
|
|
}
|