Merge pull request #71 from issp-center-dev/fix_skew

NEWS.md

@@ -1,43 +1,8 @@
 # TeNeS v1.3.x Release Notes
 
-## Changes between v1.3.2 and v1.3.1
+## Changes between v1.3.3 and v1.3.2
 
 ### Bugfixes
 
-- `tenes_simple`
-    - Antiferromagnetic initial state for S>1/2 is wrong (#61 #62)
-    - Bilinear biquadratic term $B S \cdot S$ is wrong (#63)
-    - `tenes_simple` fails to treat longer interaction and field simultaneously unless using `--use-site-hamiltonian` option (#66)
-
-## Changes between v1.3.1 and v1.3.0
-
-- Update bundled doctest.h
-    - New version supports Intel compiler 19
-
-## Changes between v1.3.0 and v1.2
-
-### New Features
-
-- `tenes_simple`
-    - Enable to generate site hamiltonian (#56)
-- `tenes_std`
-    - Introduce site hamiltonian (#56)
-    - Copy `hamiltonian` as `observable` with `group=0` when `observable` with `group=0` is not defined (#53)
-- `tenes`
-    - Introduce site imaginary time evolutionary operator (#56)
-
-### Fixes
-
-- `tenes_simple`
-    - Fix a bug that an initial state always becomes "antiferro"  in the case of honeycomb lattice (#59)
 - `tenes`
-    - Undefined observables (skipped groups) will no longer be written in `density.dat` (#60)
-
-### Documents and tutorials
-
-- Update (#52, #55, #57, #58)
-
-### For developers
-
-- Fix bugs that compilation fails in some environments (#47, #51)
-- CMake can pass `MPIEXEC_PREFLAGS` and `MPIEXEC_POSTFLAGS` options to `mpiexec` commands in test (#48)
+    - Fixed a bug in the CTMRG method for the skewed boundary condition (#71)

src/SquareLattice.cpp

@@ -16,18 +16,37 @@
 
 #include "SquareLattice.hpp"
 
+#include <cassert>
 #include <fstream>  // for operator<<, basic_ostream, char_traits, endl
 #include <tuple>    // for tuple, make_tuple, tie
 
 #include "exception.hpp"
 
+namespace {
+int gcd(int a, int b) {
+  assert(a >= 0);
+  assert(b >= 0);
+  while (b != 0) {
+    int r = a % b;
+    a = b;
+    b = r;
+  }
+  return a;
+}
+
+int lcm(int a, int b) {
+  return (a / gcd(a, b)) * b;
+}
+
+}
+
 namespace tenes {
 
 SquareLattice::SquareLattice(int X, int Y, int skew)
     : LX(X),
       LY(Y),
       N_UNIT(LX * LY),
-      skew(skew),
+      skew(skew%LX),
       physical_dims(N_UNIT, -1),
       virtual_dims(N_UNIT, std::array<int, 4>{-1, -1, -1, -1}),
       initial_dirs(N_UNIT, std::vector<double>(1)),
@@ -38,6 +57,13 @@ SquareLattice::SquareLattice(int X, int Y, int skew)
   if (Y <= 0) {
     throw tenes::input_error("Lattice.Y should be positive");
   }
+  LX_noskew = LX;
+  if(skew == 0){
+    LY_noskew = LY;
+  }else{
+    LY_noskew = LY * (::lcm(LX, skew) / skew);
+  }
+  N_UNIT_noskew = LX_noskew * LY_noskew;
   calc_neighbors();
 }
 

src/SquareLattice.hpp

@@ -53,9 +53,9 @@ namespace tenes {
  */
 class SquareLattice {
  public:
-  int LX;
-  int LY;
-  int N_UNIT;
+  int LX, LX_noskew;
+  int LY, LY_noskew;
+  int N_UNIT, N_UNIT_noskew;
 
   /*!
    * @brief Skew boundary condition

src/iTPS/core/ctm.cpp

@@ -487,11 +487,11 @@ void Left_move(std::vector<tensor> &C1, const std::vector<tensor> &C2,
      part of C1, C4, eTl will be modified */
 
   std::vector<tensor> PUs, PLs;
-  PUs.resize(lattice.LY);
-  PLs.resize(lattice.LY);
+  PUs.resize(lattice.LY_noskew);
+  PLs.resize(lattice.LY_noskew);
   int i, j, k, l;
-  for (int iy = 0; iy < lattice.LY; ++iy) {
-    i = lattice.index_fast(ix, iy);
+  for (int iy = 0; iy < lattice.LY_noskew; ++iy) {
+    i = lattice.index(ix, iy);
     j = lattice.right(i);
     k = lattice.bottom(j);
     l = lattice.left(k);
@@ -516,13 +516,13 @@ void Left_move(std::vector<tensor> &C1, const std::vector<tensor> &C2,
     eTl_bak[num] = eTl[num];
   }
   int iy_up, iy_down;
-  for (int iy = 0; iy < lattice.LY; ++iy) {
-    i = lattice.index_fast(ix, iy);
+  for (int iy = 0; iy < lattice.LY_noskew; ++iy) {
+    i = lattice.index(ix, iy);
     j = lattice.right(i);
     k = lattice.bottom(j);
     l = lattice.left(k);
-    iy_up = (iy + 1) % lattice.LY;
-    iy_down = (iy - 1 + lattice.LY) % lattice.LY;
+    iy_up = (iy + 1) % lattice.LY_noskew;
+    iy_down = (iy - 1 + lattice.LY_noskew) % lattice.LY_noskew;
 
     Calc_Next_CTM(C1_bak[i], C4_bak[l], eTt[i], eTb[l], PUs[iy_up],
                   PLs[iy_down], C1[j], C4[k]);
@@ -544,11 +544,11 @@ void Right_move(const std::vector<tensor> &C1, std::vector<tensor> &C2,
     part of C2, C3, eTr will be modified
   */
   std::vector<tensor> PUs, PLs;
-  PUs.resize(lattice.LY);
-  PLs.resize(lattice.LY);
+  PUs.resize(lattice.LY_noskew);
+  PLs.resize(lattice.LY_noskew);
   int i, j, k, l;
-  for (int iy = 0; iy < lattice.LY; ++iy) {
-    k = lattice.index_fast(ix, iy);
+  for (int iy = 0; iy < lattice.LY_noskew; ++iy) {
+    k = lattice.index(ix, iy);
     l = lattice.left(k);
     i = lattice.top(l);
     j = lattice.right(i);
@@ -577,14 +577,14 @@ void Right_move(const std::vector<tensor> &C1, std::vector<tensor> &C2,
     eTr_bak[num] = eTr[num];
   }
   int iy_up, iy_down;
-  for (int iy = 0; iy < lattice.LY; ++iy) {
-    k = lattice.index_fast(ix, iy);
+  for (int iy = 0; iy < lattice.LY_noskew; ++iy) {
+    k = lattice.index(ix, iy);
     l = lattice.left(k);
     i = lattice.top(l);
     j = lattice.right(i);
 
-    iy_up = (iy + 1) % lattice.LY;
-    iy_down = (iy - 1 + lattice.LY) % lattice.LY;
+    iy_up = (iy + 1) % lattice.LY_noskew;
+    iy_down = (iy - 1 + lattice.LY_noskew) % lattice.LY_noskew;
 
     Calc_Next_CTM(C3_bak[k], C2_bak[j], eTb[k], eTt[j], PUs[iy_down],
                   PLs[iy_up], C3[l], C2[i]);
@@ -609,11 +609,11 @@ void Top_move(std::vector<tensor> &C1, std::vector<tensor> &C2,
     ## part of C1, C2, eTt will be modified
   */
   std::vector<tensor> PUs, PLs;
-  PUs.resize(lattice.LX);
-  PLs.resize(lattice.LX);
+  PUs.resize(lattice.LX_noskew);
+  PLs.resize(lattice.LX_noskew);
   int i, j, k, l;
-  for (int ix = 0; ix < lattice.LX; ++ix) {
-    j = lattice.index_fast(ix, iy);
+  for (int ix = 0; ix < lattice.LX_noskew; ++ix) {
+    j = lattice.index(ix, iy);
     k = lattice.bottom(j);
     l = lattice.left(k);
     i = lattice.top(l);
@@ -642,14 +642,14 @@ void Top_move(std::vector<tensor> &C1, std::vector<tensor> &C2,
     eTt_bak[num] = eTt[num];
   }
   int ix_right, ix_left;
-  for (int ix = 0; ix < lattice.LX; ++ix) {
-    j = lattice.index_fast(ix, iy);
+  for (int ix = 0; ix < lattice.LX_noskew; ++ix) {
+    j = lattice.index(ix, iy);
     k = lattice.bottom(j);
     l = lattice.left(k);
     i = lattice.top(l);
 
-    ix_right = (ix + 1) % lattice.LX;
-    ix_left = (ix - 1 + lattice.LX) % lattice.LX;
+    ix_right = (ix + 1) % lattice.LX_noskew;
+    ix_left = (ix - 1 + lattice.LX_noskew) % lattice.LX_noskew;
 
     Calc_Next_CTM(C2_bak[j], C1_bak[i], eTr[j], eTl[i], PUs[ix_right],
                   PLs[ix_left], C2[k], C1[l]);
@@ -675,11 +675,11 @@ void Bottom_move(const std::vector<tensor> &C1, const std::vector<tensor> &C2,
   */
 
   std::vector<tensor> PUs, PLs;
-  PUs.resize(lattice.LX);
-  PLs.resize(lattice.LX);
+  PUs.resize(lattice.LX_noskew);
+  PLs.resize(lattice.LX_noskew);
   int i, j, k, l;
-  for (int ix = 0; ix < lattice.LX; ++ix) {
-    l = lattice.index_fast(ix, iy);
+  for (int ix = 0; ix < lattice.LX_noskew; ++ix) {
+    l = lattice.index(ix, iy);
     i = lattice.top(l);
     j = lattice.right(i);
     k = lattice.bottom(j);
@@ -709,14 +709,14 @@ void Bottom_move(const std::vector<tensor> &C1, const std::vector<tensor> &C2,
     eTb_bak[num] = eTb[num];
   }
   int ix_left, ix_right;
-  for (int ix = 0; ix < lattice.LX; ++ix) {
-    l = lattice.index_fast(ix, iy);
+  for (int ix = 0; ix < lattice.LX_noskew; ++ix) {
+    l = lattice.index(ix, iy);
     i = lattice.top(l);
     j = lattice.right(i);
     k = lattice.bottom(j);
 
-    ix_right = (ix + 1) % lattice.LX;
-    ix_left = (ix - 1 + lattice.LX) % lattice.LX;
+    ix_right = (ix + 1) % lattice.LX_noskew;
+    ix_left = (ix - 1 + lattice.LX_noskew) % lattice.LX_noskew;
 
     Calc_Next_CTM(C4_bak[l], C3_bak[k], eTl[l], eTr[k], PUs[ix_left],
                   PLs[ix_right], C4[i], C3[j]);
@@ -1119,26 +1119,26 @@ int Calc_CTM_Environment(std::vector<tensor> &C1, std::vector<tensor> &C2,
   double sig_max = 0.0;
   while ((!convergence) && (count < peps_parameters.Max_CTM_Iteration)) {
     // left move
-    for (int ix = 0; ix < lattice.LX; ++ix) {
+    for (int ix = 0; ix < lattice.LX_noskew; ++ix) {
       Left_move(C1, C2, C3, C4, eTt, eTr, eTb, eTl, Tn, ix, peps_parameters,
                 lattice);
     }
 
     // right move
-    for (int ix = 0; ix > -lattice.LX; --ix) {
+    for (int ix = 0; ix > -lattice.LX_noskew; --ix) {
       Right_move(C1, C2, C3, C4, eTt, eTr, eTb, eTl, Tn,
-                 (ix + 1 + lattice.LX) % lattice.LX, peps_parameters, lattice);
+                 (ix + 1 + lattice.LX_noskew) % lattice.LX_noskew, peps_parameters, lattice);
     }
 
     // top move
-    for (int iy = 0; iy > -lattice.LY; --iy) {
+    for (int iy = 0; iy > -lattice.LY_noskew; --iy) {
       Top_move(C1, C2, C3, C4, eTt, eTr, eTb, eTl, Tn,
-               (iy + 1 + lattice.LY) % lattice.LY, peps_parameters, lattice);
+               (iy + 1 + lattice.LY_noskew) % lattice.LY_noskew, peps_parameters, lattice);
     }
 
     // bottom move
 
-    for (int iy = 0; iy < lattice.LY; ++iy) {
+    for (int iy = 0; iy < lattice.LY_noskew; ++iy) {
       Bottom_move(C1, C2, C3, C4, eTt, eTr, eTb, eTl, Tn, iy, peps_parameters,
                   lattice);
     }