add mode to make LGS state

src/iTPS/CMakeLists.txt

@@ -21,6 +21,7 @@ add_library(
   optimize.cpp
   time_evolution.cpp
   finite_temperature.cpp
+  make_LGS.cpp
   correlation_function.cpp
   correlation_length.cpp
   transfer_matrix.cpp)

src/iTPS/PEPS_Parameters.cpp

@@ -287,6 +287,8 @@ void PEPS_Parameters::save(const char *filename, bool append) {
       ofs << "time evolution" << std::endl;
     case PEPS_Parameters::CalculationMode::finite_temperature:
       ofs << "finite temperature" << std::endl;
+    case PEPS_Parameters::CalculationMode::make_lgs:
+      ofs << "Make LGS" << std::endl;
     default:
       break;
   }

src/iTPS/PEPS_Parameters.hpp

@@ -79,6 +79,7 @@ struct PEPS_Parameters {
       ground_state,
       time_evolution,
       finite_temperature,
+      make_lgs,
   };
   CalculationMode calcmode;
 

src/iTPS/iTPS.cpp

@@ -244,8 +244,9 @@ iTPS<tensor>::iTPS(MPI_Comm comm_, PEPS_Parameters peps_parameters_,
     }
   }
 
-  if (peps_parameters.calcmode !=
-      PEPS_Parameters::CalculationMode::ground_state) {
+  if (peps_parameters.calcmode ==
+      PEPS_Parameters::CalculationMode::finite_temperature || peps_parameters.calcmode ==
+      PEPS_Parameters::CalculationMode::time_evolution) {
     if (peps_parameters.num_simple_step[0] > 0 &&
         peps_parameters.num_full_step[0]) {
       std::string msg =

src/iTPS/iTPS.hpp

@@ -141,6 +141,7 @@ class iTPS {
 
   void finite_temperature();
 
+  void make_LGS();
   //! measure expectation value of onesite observables
   std::vector<std::vector<tensor_type>> measure_onesite();
   std::vector<std::vector<tensor_type>> measure_onesite_density();

src/iTPS/load_toml.cpp

@@ -341,7 +341,9 @@ PEPS_Parameters gen_param(decltype(cpptoml::parse_file("")) param) {
       pparam.calcmode = PEPS_Parameters::CalculationMode::time_evolution;
     } else if (util::startswith(mode_str, "finite")) {
       pparam.calcmode = PEPS_Parameters::CalculationMode::finite_temperature;
-    } else {
+    } else if (util::startswith(mode_str, "lgs")) {
+      pparam.calcmode = PEPS_Parameters::CalculationMode::make_lgs;
+    }else {
       throw input_error("Invalid mode: " + mode_str);
     }
     load_if(pparam.measure_interval, general, "measure_interval");

src/iTPS/main.cpp

@@ -195,6 +195,24 @@ int run_finitetemperature(MPI_Comm comm, PEPS_Parameters peps_parameters,
   return 0;
 }
 
+template <class tensor>
+int run_make_LGS(MPI_Comm comm, PEPS_Parameters peps_parameters,
+			  SquareLattice lattice,
+			  EvolutionOperators<tensor> simple_updates,
+			  EvolutionOperators<tensor> full_updates,
+			  Operators<tensor> onesite_operators,
+			  Operators<tensor> twosite_operators,
+			  Operators<tensor> multisite_operators,
+			  CorrelationParameter corparam,
+			  TransferMatrix_Parameters clength_param) {
+  iTPS<tensor> tns(comm, peps_parameters, lattice, simple_updates, full_updates,
+                   onesite_operators, twosite_operators, multisite_operators,
+                   corparam, clength_param);
+  tns.make_LGS();
+  tns.save_tensors();
+  tns.summary();
+  return 0;
+}
 
 
 int itps_main(const char* input_filename, MPI_Comm comm,
@@ -341,6 +359,11 @@ int itps_main(std::string input_filename, MPI_Comm comm,
                              full_updates, onesite_obs, twosite_obs,
                              multisite_obs, corparam, clength_param);
     }
+  } else if (peps_parameters.calcmode ==
+             PEPS_Parameters::CalculationMode::make_lgs) {
+    return run_make_LGS(comm, peps_parameters, lattice, simple_updates,
+                             full_updates, onesite_obs, twosite_obs,
+                             multisite_obs, corparam, clength_param);
   } else {
 
     return 1;

src/iTPS/make_LGS.cpp

@@ -0,0 +1,264 @@
+/* TeNeS - Massively parallel tensor network solver /
+/ Copyright (C) 2019- The University of Tokyo */
+
+/* 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 /
+/ 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 General Public License for more details. */
+
+/* You should have received a copy of the GNU General Public License /
+/ along with this program. If not, see http://www.gnu.org/licenses/. */
+
+#include "iTPS.hpp"
+#include <cmath>
+
+namespace tenes {
+namespace itps {
+template <class tensor>
+void iTPS<tensor>::make_LGS() {
+  using mptensor::Shape, mptensor::Axes, mptensor::Index;
+
+  // initialize iTPS for LGS //  
+  std::vector<tensor> LGS;
+  LGS.clear(); 
+  for (int i = 0; i < N_UNIT; ++i) {
+    const auto pdim = lattice.physical_dims[i];
+    if (pdim > 3){
+      std::stringstream ss;
+      ss << "ERROR: LGS support only S=1/2 and S=1" << std::endl;
+      throw std::runtime_error(ss.str());
+    }
+    LGS.push_back(tensor(comm, Shape(1, 1, 1, 1, pdim)));
+  }
+
+  // assuming all pdim s are identical
+  int pdim;
+  pdim = lattice.physical_dims[0];
+  std::vector<double> real_part,img_part;
+  double factor;
+  factor = (std::sqrt(3.0) - 1.0) * 0.5;
+  Index index;
+
+  if (pdim == 2){
+    // S = 1/2
+    real_part.push_back(1.0);
+    real_part.push_back(factor);
+
+    img_part.push_back(0.0);
+    img_part.push_back(factor);
+
+    for (int i = 0; i < lattice.N_UNIT; ++i){
+      for (int n = 0; n < LGS[i].local_size(); ++n){
+	index = LGS[i].global_index(n);
+	// in the present case, local_size == pdim
+	auto v = std::complex<double>(real_part[index[4]], img_part[index[4]]);
+	LGS[i].set_value(index, to_tensor_type(v));
+      }
+    }
+
+    // projectors
+    tensor Q = tensor(comm, Shape(2, 2, 2, pdim, pdim));
+
+    // Identity
+    Q.set_value(Index(0 ,0 ,0, 0, 0), to_tensor_type(1.0));
+    Q.set_value(Index(0 ,0 ,0, 1, 1), to_tensor_type(1.0));
+
+    // sigma_x
+    Q.set_value(Index(0 ,1 ,1, 0, 1), to_tensor_type(1.0));
+    Q.set_value(Index(0 ,1 ,1, 1, 0), to_tensor_type(1.0));
+
+    // sigma_y
+    Q.set_value(Index(1 ,0 ,1, 0, 1), to_tensor_type(std::complex<double>(0.0, -1.0)));
+    Q.set_value(Index(1 ,0 ,1, 1, 0), to_tensor_type(std::complex<double>(0.0, 1.0)));
+
+    // sigma_z
+    Q.set_value(Index(1 ,1 ,0, 0, 0), to_tensor_type(1.0));
+    Q.set_value(Index(1 ,1 ,0, 1, 1), to_tensor_type(-1.0));
+
+    if (lattice.initial_dirs[0][0] == lattice.initial_dirs[1][0]){
+      // ferro LGS
+      std::cout<< "Creating Ferro LGS for S = 1/2"<<std::endl;
+      for (int i = 0; i < lattice.N_UNIT; ++i){
+	int ix = i % LX;
+	int iy = i / LX;
+	if ((ix + iy) % 2 == 0){
+	  // sublattice 1
+	  LGS[i] = mptensor::reshape(mptensor::tensordot(Q, LGS[i], Axes(4), Axes(4)).transpose(Axes(1, 0, 2, 3, 4, 5, 6, 7)), Shape(2, 1, 2, 2, 2));
+	} else{
+	  // sublattice 2
+	  LGS[i] = mptensor::reshape(mptensor::tensordot(Q, LGS[i], Axes(4), Axes(4)).transpose(Axes(0, 2, 1, 3, 4, 5, 6, 7)), Shape(2, 2, 2, 1, 2));
+	}
+      }
+    } else {
+      // antiferro LGS
+      std::cout<< "Creating Antiferro LGS for S = 1/2"<<std::endl;
+      tensor v = tensor(comm, Shape(2, 2));
+      tensor vd = tensor(comm, Shape(2, 2));
+
+      v.set_value(Index(0,1), to_tensor_type(std::complex<double>(0.0, 1.0)));
+      v.set_value(Index(1,0), to_tensor_type(std::complex<double>(1.0, 0.0)));
+
+      vd.set_value(Index(0,1), to_tensor_type(std::complex<double>(1.0, 0.0)));
+      vd.set_value(Index(1,0), to_tensor_type(std::complex<double>(0.0, -1.0)));
+
+      for (int i = 0; i < lattice.N_UNIT; ++i){
+	int ix = i % LX;
+	int iy = i / LX;
+	if ((ix -  iy + LX) % 4 == 0){
+	  // sublattice 1
+	  LGS[i] = mptensor::reshape(mptensor::tensordot(
+							 v, mptensor::tensordot(
+									     vd, mptensor::tensordot(
+												 v, mptensor::tensordot(
+														     Q, LGS[i], Axes(4), Axes(4)), 
+												 Axes(1), Axes(2)),
+									     Axes(0), Axes(1)),
+							 Axes(0), Axes(2)),
+				     Shape(2, 1, 2, 2, 2));
+
+	} else if (((ix - iy + LX) % 4) % 2 == 1){
+	  // sublattice 2 or 4
+	  LGS[i] = mptensor::reshape(mptensor::tensordot(Q, LGS[i], Axes(4), Axes(4)).transpose(Axes(0, 2, 1, 3, 4, 5, 6, 7)), Shape(2, 2, 2, 1, 2));
+	} else if ((ix - iy + LX) % 4 == 2){
+	  // sublattice 3
+	  LGS[i] = mptensor::reshape(mptensor::tensordot(
+							 vd, mptensor::tensordot(
+									     v, mptensor::tensordot(
+												 vd, mptensor::tensordot(
+														     Q, LGS[i], Axes(4), Axes(4)), 
+												 Axes(1), Axes(2)),
+									     Axes(1), Axes(1)),
+							 Axes(0), Axes(2)),
+				     Shape(2, 1, 2, 2, 2));
+	}
+      }
+    }
+
+  } else {
+    // S = 1
+    real_part.push_back(1.0);
+    real_part.push_back(std::sqrt(2.0) * factor);
+    real_part.push_back(0.0);
+
+    img_part.push_back(0.0);
+    img_part.push_back(std::sqrt(2.0) * factor);
+    img_part.push_back(2.0 * factor * factor);
+
+    for (int i = 0; i < lattice.N_UNIT; ++i){
+      for (int n = 0; n < LGS[i].local_size(); ++n){
+	index = LGS[i].global_index(n);
+	// in the present case, local_size == pdim
+	auto v = std::complex<double>(real_part[index[4]], img_part[index[4]]);
+	LGS[i].set_value(index, to_tensor_type(v));
+      }
+    }
+    // projectors
+    tensor Q = tensor(comm, Shape(2, 2, 2, pdim, pdim));
+
+    // Identity
+    Q.set_value(Index(0 ,0 ,0, 0, 0), to_tensor_type(1.0));
+    Q.set_value(Index(0 ,0 ,0, 1, 1), to_tensor_type(1.0));
+    Q.set_value(Index(0 ,0 ,0, 2, 2), to_tensor_type(1.0));
+
+    // sigma_x
+    Q.set_value(Index(0 ,1 ,1, 0, 2), to_tensor_type(-1.0));
+    Q.set_value(Index(0 ,1 ,1, 1, 1), to_tensor_type(-1.0));
+    Q.set_value(Index(0 ,1 ,1, 2, 0), to_tensor_type(-1.0));
+
+    // sigma_y
+    Q.set_value(Index(1 ,0 ,1, 0, 2), to_tensor_type(1.0));
+    Q.set_value(Index(1 ,0 ,1, 1, 1), to_tensor_type(-1.0));
+    Q.set_value(Index(1 ,0 ,1, 2, 0), to_tensor_type(1.0));
+
+    // sigma_z
+    Q.set_value(Index(1 ,1 ,0, 0, 0), to_tensor_type(-1.0));
+    Q.set_value(Index(1 ,1 ,0, 1, 1), to_tensor_type(1.0));
+    Q.set_value(Index(1 ,1 ,0, 2, 2), to_tensor_type(-1.0));
+
+
+    if (lattice.initial_dirs[0][0] == lattice.initial_dirs[1][0]){
+      // ferro LGS
+      std::cout<< "Creating Ferro LGS for S = 1"<<std::endl;
+
+      for (int i = 0; i < lattice.N_UNIT; ++i){
+	int ix = i % LX;
+	int iy = i / LX;
+	if ((ix + iy) % 2 == 0){
+	  // sublattice 1
+	  LGS[i] = mptensor::reshape(mptensor::tensordot(Q, LGS[i], Axes(4), Axes(4)).transpose(Axes(1, 0, 2, 3, 4, 5, 6, 7)), Shape(2, 1, 2, 2, 3));
+	} else{
+	  // sublattice 2
+	  LGS[i] = mptensor::reshape(mptensor::tensordot(Q, LGS[i], Axes(4), Axes(4)).transpose(Axes(0, 2, 1, 3, 4, 5, 6, 7)), Shape(2, 2, 2, 1, 3));
+	}
+      }
+    } else {
+      // antiferro LGS
+      std::cout<< "Creating Antiferro LGS for S = 1"<<std::endl;
+      tensor x = tensor(comm, Shape(2, 2));
+
+      x.set_value(Index(0,1), to_tensor_type(1.0));
+      x.set_value(Index(1,0), to_tensor_type(1.0));
+
+      for (int i = 0; i < lattice.N_UNIT; ++i){
+	int ix = i % LX;
+	int iy = i / LX;
+	if ((ix + iy) % 2 == 0){
+	  // sublattice 1
+	  LGS[i] = mptensor::reshape(mptensor::tensordot(
+							 x, mptensor::tensordot(
+									     x, mptensor::tensordot(
+												 x, mptensor::tensordot(
+														     Q, LGS[i], Axes(4), Axes(4)), 
+												 Axes(1), Axes(2)),
+									     Axes(1), Axes(1)),
+							 Axes(1), Axes(2)),
+				     Shape(2, 1, 2, 2, 3));
+	} else{
+	  // sublattice 2
+	  LGS[i] = mptensor::reshape(mptensor::tensordot(Q, LGS[i], Axes(4), Axes(4)).transpose(Axes(0, 2, 1, 3, 4, 5, 6, 7)), Shape(2, 2, 2, 1, 3));
+	}
+      }
+    }
+  }
+
+  // copy to Tn
+
+  Tn.clear();
+  for (int i = 0; i < N_UNIT; ++i) {
+    const auto pdim = lattice.physical_dims[i];
+    const auto vdim = lattice.virtual_dims[i];
+
+    Tn.push_back(tensor(comm, Shape(vdim[0], vdim[1], vdim[2], vdim[3], pdim)));
+  }
+
+  for (int i = 0; i < lattice.N_UNIT; ++i) {
+    const auto pdim = lattice.physical_dims[i];
+    const auto vdim = lattice.virtual_dims[i];
+
+
+    for (int n = 0; n < Tn[i].local_size(); ++n) {
+      index = Tn[i].global_index(n);
+      if (index[0] < 2 && index[1] < 2 && index[2] < 2 && index[3] < 2) {
+	typename tensor::value_type v;
+	bool test;
+	test = LGS[i].get_value(index, v);
+	Tn[i].set_value(index, to_tensor_type(v));
+      } else {
+	Tn[i].set_value(index, to_tensor_type(0.0));
+      }
+    }
+  }
+  // output tensors
+
+}
+// template specialization
+template class iTPS<real_tensor>;
+template class iTPS<complex_tensor>;
+
+}  // namespace itps
+}  // namespace tenes