TRIQS / triqs

a Toolbox for Research on Interacting Quantum Systems
https://triqs.github.io
GNU General Public License v3.0
139 stars 72 forks source link

Clef assignment with repeated indices assigns erroneously and surprisingly #475

Closed HugoStrand closed 5 years ago

HugoStrand commented 7 years ago

I am trying to build a tensor two-particle Green's function from a sampled block2 Green's function on AABB block form. To do so I want to first populate the AABB blocks of the tensor Green's function using a simple clef expression

g2_tensor(n1, n2, n3)(A,A,B,B) << block_gf2(A,B)(n1, n2, n3)(0, 0, 0, 0);

The problem is that this expression (surprisingly) does not take into account that A and B are repeated on the left hand side. Hence, it also assigns to, e.g. g2_tensor(n1, n2, n3)(0, 0, 0, 1), i.e. not following the (A, A, B, B) index constraint given in the above clef expression.

This is a serious bug. Is it possible to fix it?

Seemingly we need more test coverage of the clef machinery, debugging it while trying to solve the BSE is a bit painful. The current workaround is plain stupid

  g2_tensor(n1, n2, n3)(A,B,C,D) << kronecker(A,B)*kronecker(C,D)*block_gf2(A,C)(n1, n2, n3)(0, 0, 0, 0);

Best, Hugo

Below is a test that currently fails for unstable to be put in ./test/triqs/gfs/multivar/clef_block_tensor.cpp

#include <triqs/test_tools/gfs.hpp>

using namespace triqs::clef;

placeholder<7> n1;
placeholder<8> n2;
placeholder<9> n3;

placeholder<10> A;
placeholder<11> B;

TEST(Gf, ClefBlockTensor) {

  double beta = 1.234;

  int N = 3;
  auto m = gf_mesh<imfreq>{beta, Fermion, N};
  auto mp = gf_mesh<cartesian_product<imfreq, imfreq, imfreq>>{m, m, m};

  auto g2_tensor = gf<cartesian_product<imfreq, imfreq, imfreq>, tensor_valued<4>>
    {{m, m, m}, {2, 2, 2, 2}};

  auto g2 = gf<cartesian_product<imfreq, imfreq, imfreq>, tensor_valued<4>>
    {{m, m, m}, {1, 1, 1, 1}};
  auto block_gf2 = make_block2_gf(2, 2, g2);

  // -- setup some values
  block_gf2(A, B)(n1, n2, n3) << 1. + 10. * A + 100. * B + (n1 + 10*n2 + 100*n3);
  g2_tensor *= 0.0;

  // -- Assign with clef using repeated indices
  // -- setting only blocks with structure AABB

  g2_tensor(n1, n2, n3)(A,A,B,B) << block_gf2(A,B)(n1, n2, n3)(0, 0, 0, 0);

  /*
  // Debug print
  for( auto i1 : range(N) ) {
  for( auto i2 : range(N) ) {
  for( auto i3 : range(N) ) {
    std::cout << i1 << ", " << i2 << ", " << i3 << " : " << g2_tensor(i1, i2, i3)(0, 0, 0, 1) << "\n";
  }}}
  */

  // AAAB blocks not set, should be zero
  EXPECT_CLOSE(g2_tensor(0, 0, 0)(0, 0, 0, 1), 0.0); 
  EXPECT_CLOSE(g2_tensor(0, 0, 0)(0, 0, 1, 0), 0.0); 
  EXPECT_CLOSE(g2_tensor(0, 0, 0)(0, 1, 0, 0), 0.0); 
  EXPECT_CLOSE(g2_tensor(0, 0, 0)(1, 0, 0, 0), 0.0); 

  // ABBA blocks not set, should be zero
  EXPECT_CLOSE(g2_tensor(0, 0, 0)(0, 1, 1, 0), 0.0); 
  EXPECT_CLOSE(g2_tensor(0, 0, 0)(1, 0, 0, 1), 0.0); 

}

MAKE_MAIN;
parcollet commented 7 years ago

One was not supposed to use the same placeholder twice in the LHS, the lib should reject this. Anyway, it should be fixed, or not compile. @Wentzell : or use 17 for loop + structured bindings.

parcollet commented 7 years ago

a(i,j) << anything (i,j) is rewritten triqs_auto_assign( a, [](auto i, auto j) { return anything(i,j);}) It will iterate on the whole domain, so the placeholders must be different.

parcollet commented 5 years ago

Fixed in unstable