Port to Pomerol 2.0

.github/workflows/CI.yml

@@ -56,17 +56,32 @@ jobs:
  python3-mako
  libnfft3-dev
 
+ - name: Install libcommute
+ env:
+ CXX: ${{ matrix.compiler.cxx }}
+ run: |
+ git clone https://github.com/krivenko/libcommute libcommute
+ mkdir libcommute/build && pushd libcommute/build
+ cmake .. \
+ -DCMAKE_INSTALL_PREFIX=$HOME/install \
+ -DTESTS=OFF \
+ -DEXAMPLES=OFF
+ make install
+ popd
+
  - name: Build & install Pomerol
  env:
  CC: ${{ matrix.cc }}
  CXX: ${{ matrix.cxx }}
  run: |
  git clone https://github.com/aeantipov/pomerol
  mkdir pomerol/build && pushd pomerol/build
- cmake .. \
- -DCMAKE_BUILD_TYPE=Debug \
- -DCMAKE_INSTALL_PREFIX=$HOME/install \
- -DTesting=OFF
+ cmake .. \
+ -DCMAKE_BUILD_TYPE=Debug \
+ -DCMAKE_INSTALL_PREFIX=$HOME/install \
+ -Dlibcommute_DIR=$HOME/install/lib/cmake \
+ -DTesting=OFF \
+ -DDocumentation=OFF
  make -j2 install VERBOSE=1
  popd
 

CMakeLists.txt

@@ -1,5 +1,5 @@
 # Version number of the application
-set(POMEROL2TRIQS_VERSION "0.5")
+set(POMEROL2TRIQS_VERSION "0.6")
 
 # Start configuration
 cmake_minimum_required(VERSION 3.8 FATAL_ERROR)
@@ -40,7 +40,8 @@ if(PythonSupport AND NOT TRIQS_WITH_PYTHON_SUPPORT)
 endif()
 
 # Load Pomerol
-find_package(pomerol REQUIRED)
+find_package(pomerol 2.0 CONFIG REQUIRED)
+message(STATUS "Found pomerol ${pomerol_VERSION}")
 
 # Find MPI
 find_package(MPI)

README.md

@@ -24,7 +24,8 @@ Notation for the two-particle Green's functions is adopted from the
 Installation
 ------------
 
-- Install the latest version of [Pomerol](http://aeantipov.github.io/pomerol/) exact diagonalization library (`master` branch).
+- Install the [Pomerol](http://aeantipov.github.io/pomerol/) exact diagonalization library, version 2.0 or newer
+ **(Pomerol 1.x is supported only by pomerol2triqs 0.5 and older)**.
 - Install the [TRIQS](http://triqs.github.io/triqs/3.0.x/install.html) library version 3.0.x.
 - `source <path_to_triqs_install_dir>/share/triqsvars.sh`
 - `git clone https://github.com/krivenko/pomerol2triqs.git pomerol2triqs.git`
@@ -37,7 +38,7 @@ Installation
 License
 -------
 
-Copyright (C) 2017-2020 Igor Krivenko <igor.s.krivenko @ gmail.com>
+Copyright (C) 2017-2021 Igor Krivenko <igor.s.krivenko @ gmail.com>
 
 With contributions from Hugo U.R. Strand and Nils Wentzell.
 

c++/CMakeLists.txt

@@ -1,10 +1,6 @@
 add_library(pomerol2triqs_c pomerol_ed.cpp g.cpp g2.cpp chi.cpp)
 
-find_package(Boost COMPONENTS mpi serialization)
-
-target_link_libraries(pomerol2triqs_c PUBLIC triqs ${pomerol_LIBRARIES})
-
-target_include_directories(pomerol2triqs_c PUBLIC ${pomerol_INCLUDE_DIRS})
+target_link_libraries(pomerol2triqs_c PUBLIC triqs pomerol::pomerol)
 target_compile_options(pomerol2triqs_c PRIVATE
  -DPOMEROL2TRIQS_GIT_HASH=${POMEROL2TRIQS_GIT_HASH}
  -DPOMEROL2TRIQS_GIT_HASH=${POMEROL2TRIQS_GIT_HASH}

c++/chi.cpp

@@ -1,7 +1,7 @@
 /**
  * pomerol2triqs
  *
- * Copyright (C) 2017-2020 Igor Krivenko <igor.s.krivenko @ gmail.com>
+ * Copyright (C) 2017-2021 Igor Krivenko <igor.s.krivenko @ gmail.com>
  *
  * This program is free software: you can redistribute it and/or modify
  * it under the terms of the GNU General Public License as published by
@@ -31,23 +31,23 @@ namespace pomerol2triqs {
 
  Pomerol::ParticleIndex pom_i = lookup_pomerol_index(i);
  if (pom_i == -1)
- TRIQS_RUNTIME_ERROR << "ensemble_average: unexpected index i = " << i;
+ TRIQS_RUNTIME_ERROR << "ensemble_average: Unexpected index i = " << i;
  Pomerol::ParticleIndex pom_j = lookup_pomerol_index(j);
  if (pom_j == -1)
- TRIQS_RUNTIME_ERROR << "ensemble_average: unexpected index j = " << j;
+ TRIQS_RUNTIME_ERROR << "ensemble_average: Unexpected index j = " << j;
 
  if (!matrix_h)
- TRIQS_RUNTIME_ERROR << "ensemble_average: no Hamiltonian has been diagonalized";
+ TRIQS_RUNTIME_ERROR << "ensemble_average: No Hamiltonian has been diagonalized";
  compute_rho(beta);
 
- Pomerol::QuadraticOperator op(index_info, *states_class, *matrix_h, pom_i, pom_j);
- op.prepare();
+ Pomerol::QuadraticOperator op(index_info, *hs, *states_class, *matrix_h, pom_i, pom_j);
+ op.prepare(*hs);
  op.compute();
 
- Pomerol::EnsembleAverage EA(*states_class, *matrix_h, op, *rho);
- EA.prepare();
+ Pomerol::EnsembleAverage EA(op, *rho);
+ EA.compute();
 
- return EA.getResult();
+ return EA();
  }
 
  template <typename Mesh, typename Filler>
@@ -59,12 +59,12 @@ namespace pomerol2triqs {
  gf_mesh<Mesh> const &mesh,
  Filler filler) const {
  if (!states_class || !matrix_h || !rho)
- TRIQS_RUNTIME_ERROR << "compute_chi: internal error!";
+ TRIQS_RUNTIME_ERROR << "compute_chi: Internal error!";
 
  auto checked_lookup = [&](indices_t const &i) {
  Pomerol::ParticleIndex pom_i = lookup_pomerol_index(i);
  if (pom_i == -1)
- TRIQS_RUNTIME_ERROR << "compute_chi: unexpected index " << i;
+ TRIQS_RUNTIME_ERROR << "compute_chi: Unexpected index " << i;
  return pom_i;
  };
 
@@ -73,11 +73,11 @@ namespace pomerol2triqs {
  Pomerol::ParticleIndex pom_i2 = checked_lookup(i2);
  Pomerol::ParticleIndex pom_j2 = checked_lookup(j2);
 
- Pomerol::QuadraticOperator A(index_info, *states_class, *matrix_h, pom_i1, pom_j1);
- Pomerol::QuadraticOperator B(index_info, *states_class, *matrix_h, pom_i2, pom_j2);
+ Pomerol::QuadraticOperator A(index_info, *hs, *states_class, *matrix_h, pom_i1, pom_j1);
+ Pomerol::QuadraticOperator B(index_info, *hs, *states_class, *matrix_h, pom_i2, pom_j2);
 
- A.prepare(); A.compute();
- B.prepare(); B.compute();
+ A.prepare(*hs); A.compute();
+ B.prepare(*hs); B.compute();
 
  Pomerol::Susceptibility pom_chi(*states_class, *matrix_h, A, B, *rho);
  pom_chi.prepare();
@@ -95,7 +95,7 @@ namespace pomerol2triqs {
  indices_t const &i2,
  indices_t const &j2,
  double beta, int n_tau, bool connected) {
- if (!matrix_h) TRIQS_RUNTIME_ERROR << "chi_tau: no Hamiltonian has been diagonalized";
+ if (!matrix_h) TRIQS_RUNTIME_ERROR << "chi_tau: No Hamiltonian has been diagonalized";
  compute_rho(beta);
 
  auto filler = [](gf_view<imtime, scalar_valued> chi, Pomerol::Susceptibility const &pom_chi) {
@@ -109,7 +109,7 @@ namespace pomerol2triqs {
  indices_t const &i2,
  indices_t const &j2,
  double beta, int n_inu, bool connected) {
- if (!matrix_h) TRIQS_RUNTIME_ERROR << "chi_inu: no Hamiltonian has been diagonalized";
+ if (!matrix_h) TRIQS_RUNTIME_ERROR << "chi_inu: No Hamiltonian has been diagonalized";
  compute_rho(beta);
 
  auto filler = [](gf_view<imfreq, scalar_valued> chi, Pomerol::Susceptibility const &pom_chi) {

c++/g.cpp

@@ -1,7 +1,7 @@
 /**
  * pomerol2triqs
  *
- * Copyright (C) 2017-2020 Igor Krivenko <igor.s.krivenko @ gmail.com>
+ * Copyright (C) 2017-2021 Igor Krivenko <igor.s.krivenko @ gmail.com>
  *
  * This program is free software: you can redistribute it and/or modify
  * it under the terms of the GNU General Public License as published by
@@ -28,7 +28,7 @@ namespace pomerol2triqs {
  template <typename Mesh, typename Filler>
  block_gf<Mesh> pomerol_ed::compute_gf(gf_struct_t const &gf_struct, gf_mesh<Mesh> const &mesh, Filler filler) const {
 
- if (!states_class || !matrix_h || !rho || !ops_container) TRIQS_RUNTIME_ERROR << "compute_gf: internal error!";
+ if (!states_class || !matrix_h || !rho || !ops_container) TRIQS_RUNTIME_ERROR << "compute_gf: Internal error!";
 
  struct index_visitor {
  std::vector<std::string> indices;
@@ -56,7 +56,7 @@ namespace pomerol2triqs {
  for (int i2 : range(n)) {
  Pomerol::ParticleIndex pom_i2 = lookup_pomerol_index({bl.first, bl.second[i2]});
 
- if (verbose && !comm.rank())
+ if (verbose && !pMPI::rank(comm))
  std::cout << "fill_gf: Filling GF component (" << bl.first << "," << bl.second[i1] << ")(" << bl.first << "," << bl.second[i2] << ")"
  << std::endl;
  auto g_el = slice_target_to_scalar(g, i1, i2);
@@ -74,7 +74,7 @@ namespace pomerol2triqs {
  }
 
  block_gf<imfreq> pomerol_ed::G_iw(gf_struct_t const &gf_struct, double beta, int n_iw) {
- if (!matrix_h) TRIQS_RUNTIME_ERROR << "G_iw: no Hamiltonian has been diagonalized";
+ if (!matrix_h) TRIQS_RUNTIME_ERROR << "G_iw: No Hamiltonian has been diagonalized";
  compute_rho(beta);
  compute_field_operators(gf_struct);
 
@@ -85,7 +85,7 @@ namespace pomerol2triqs {
  }
 
  block_gf<imtime> pomerol_ed::G_tau(gf_struct_t const &gf_struct, double beta, int n_tau) {
- if (!matrix_h) TRIQS_RUNTIME_ERROR << "G_tau: no Hamiltonian has been diagonalized";
+ if (!matrix_h) TRIQS_RUNTIME_ERROR << "G_tau: No Hamiltonian has been diagonalized";
  compute_rho(beta);
  compute_field_operators(gf_struct);
 
@@ -97,7 +97,7 @@ namespace pomerol2triqs {
 
  block_gf<refreq> pomerol_ed::G_w(gf_struct_t const &gf_struct, double beta, std::pair<double, double> const &energy_window, int n_w,
  double im_shift) {
- if (!matrix_h) TRIQS_RUNTIME_ERROR << "G_w: no Hamiltonian has been diagonalized";
+ if (!matrix_h) TRIQS_RUNTIME_ERROR << "G_w: No Hamiltonian has been diagonalized";
  compute_rho(beta);
  compute_field_operators(gf_struct);
 

c++/g2.cpp

@@ -1,7 +1,7 @@
 /**
  * pomerol2triqs
  *
- * Copyright (C) 2017-2020 Igor Krivenko <igor.s.krivenko @ gmail.com>
+ * Copyright (C) 2017-2021 Igor Krivenko <igor.s.krivenko @ gmail.com>
  *
  * This program is free software: you can redistribute it and/or modify
  * it under the terms of the GNU General Public License as published by
@@ -50,7 +50,7 @@ namespace pomerol2triqs {
  auto pomerol_ed::compute_g2(gf_struct_t const &gf_struct, gf_mesh<Mesh> const &mesh, block_order_t block_order, g2_blocks_t const &g2_blocks,
  Filler filler) const -> block2_gf<Mesh, tensor_valued<4>> {
 
- if (!states_class || !matrix_h || !rho || !ops_container) TRIQS_RUNTIME_ERROR << "compute_g2: internal error!";
+ if (!states_class || !matrix_h || !rho || !ops_container) TRIQS_RUNTIME_ERROR << "compute_g2: Internal error!";
 
  bool compute_all_blocks = g2_blocks.empty();
 
@@ -85,8 +85,8 @@ namespace pomerol2triqs {
  for (int c : range(B_size))
  for (int d : range(B_size)) {
 
- if (verbose && !comm.rank()) {
- std::cout << "compute_g2: filling G^2 element ";
+ if (verbose && !pMPI::rank(comm)) {
+ std::cout << "compute_g2: Filling G^2 element ";
  if (block_order == AABB) {
  std::cout << "(" << A << "," << a << ")";
  std::cout << "(" << A << "," << b << ")";
@@ -125,16 +125,16 @@ namespace pomerol2triqs {
  }
 
  auto pomerol_ed::G2_iw_inu_inup(g2_iw_inu_inup_params_t const &p) -> block2_gf<w_nu_nup_t, tensor_valued<4>> {
- if (!matrix_h) TRIQS_RUNTIME_ERROR << "G2_iw_inu_inup: no Hamiltonian has been diagonalized";
+ if (!matrix_h) TRIQS_RUNTIME_ERROR << "G2_iw_inu_inup: No Hamiltonian has been diagonalized";
  compute_rho(p.beta);
  compute_field_operators(p.gf_struct);
 
- if (verbose && !comm.rank()) std::cout << "G2_iw_inu_inup: filling output container" << std::endl;
+ if (verbose && !pMPI::rank(comm)) std::cout << "G2_iw_inu_inup: Filling output container" << std::endl;
 
  auto filler = [&p, this](gf_view<w_nu_nup_t, scalar_valued> g2_el, auto const &pom_g2) {
  long mesh_index = 0;
  for (auto w_nu_nup : g2_el.mesh()) {
- if ((mesh_index++) % comm.size() != comm.rank()) continue;
+ if ((mesh_index++) % pMPI::size(comm) != pMPI::rank(comm)) continue;
 
  if (p.channel == AllFermionic) {
 
@@ -183,13 +183,13 @@ namespace pomerol2triqs {
  }
 
  auto pomerol_ed::G2_iw_l_lp(g2_iw_l_lp_params_t const &p) -> block2_gf<w_l_lp_t, tensor_valued<4>> {
- if (!matrix_h) TRIQS_RUNTIME_ERROR << "G2_iw_l_lp: no Hamiltonian has been diagonalized";
+ if (!matrix_h) TRIQS_RUNTIME_ERROR << "G2_iw_l_lp: No Hamiltonian has been diagonalized";
  compute_rho(p.beta);
  compute_field_operators(p.gf_struct);
 
  gf_mesh<w_l_lp_t> mesh{{p.beta, Boson, p.n_iw}, {p.beta, Fermion, static_cast<size_t>(p.n_l)}, {p.beta, Fermion, static_cast<size_t>(p.n_l)}};
 
- if (verbose && !comm.rank()) std::cout << "G2_iw_l_lp: filling output container" << std::endl;
+ if (verbose && !pMPI::rank(comm)) std::cout << "G2_iw_l_lp: Filling output container" << std::endl;
 
  auto filler = [&p, this](gf_view<w_l_lp_t, scalar_valued> g2_el, auto const &pom_g2) {
 
@@ -207,7 +207,7 @@ namespace pomerol2triqs {
 
  long mesh_index = 0;
  for (auto iw : std::get<0>(g2_el.mesh())) {
- if((mesh_index++) % comm.size() != comm.rank()) continue;
+ if((mesh_index++) % pMPI::size(comm) != pMPI::rank(comm)) continue;
 
  int w_m = iw.index();