Merge pull request #91 from issp-center-dev/coeff

enable to specify coeff of observable

docs/sphinx/en/file_specification/observable_section.rst

@@ -18,6 +18,8 @@ Define one-body operators that indicate physical quantities defined at each site
    ``sites``,    "Site number",         Integer or a list of integer
    ``dim``,      "Dimension of an operator",       Integer
    ``elements``, "Non-zero elements of an operator", String
+   ``coeff``,    "Coefficient of operator (real part)", Float
+   ``coeff_im``,  "Coefficient of operator (imaginary part)", Float
 
 ``name``  specifies an operator name.
 
@@ -35,6 +37,9 @@ One element is specified by one line consisting of two integers and two floating
 - The first two integers are the state numbers before and after the act of the operator, respectively.
 - The latter two floats indicate the real and imaginary parts of the elements of the operator, respectively.
 
+``coeff`` and ``coeff_im`` are real and imaginary parts of the coefficient of the operator, respectively.
+If omitted, they are set to 1.0 and 0.0, respectively.
+
 Example
 .......
 
@@ -93,6 +98,8 @@ Define two-body operators that indicate physical quantities defined on two sites
    ``dim``,      "Dimension of an operator",           Integer
    ``elements``, "Non-zero elements of an operator",   String
    ``ops``,      "Index of onesite operators",         A list of integer
+   ``coeff``,    "Coefficient of operator (real part)", Float
+   ``coeff_im``,  "Coefficient of operator (imaginary part)", Float
 
 
 ``name``  specifies an operator name.
@@ -123,6 +130,9 @@ For example, if :math:`S^z` is defined as ``group = 0`` in ``observable.onesite`
 
 If both ``elements`` and ``ops`` are defined, the process will end in error.
 
+``coeff`` and ``coeff_im`` are real and imaginary parts of the coefficient of the operator, respectively.
+If omitted, they are set to 1.0 and 0.0, respectively.
+
 Example
 .......
 
@@ -193,6 +203,8 @@ It is defined as a direct product of one-body operators defined in ``observable.
    ``group``,      "Identification number of operators", Integer
    ``multisites``, "Sites",                              String
    ``ops``,        "Index of onesite operators",         List of integers
+   ``coeff``,      "Coefficient of operator (real part)", Float
+   ``coeff_im``,    "Coefficient of operator (imaginary part)", Float
 
 ``name``  specifies an operator name.
 
@@ -210,3 +222,5 @@ One line consisting of integers means a set sites.
 Using ``ops``, a multi-body operator can be defined as a direct product of the one-body operators defined in ``observable.onesite``.
 For example, if :math:`S^z` is defined as ``group = 0`` in ``observable.onesite``,  :math:`S^z_i S^z_j S^z_k` can be expressed as ``ops = [0,0,0]``.
 
+``coeff`` and ``coeff_im`` are real and imaginary parts of the coefficient of the operator, respectively.
+If omitted, they are set to 1.0 and 0.0, respectively.

docs/sphinx/ja/file_specification/observable_section.rst

@@ -14,11 +14,13 @@
    :header: "名前", "説明", "型"
    :widths: 15, 30, 20
 
-   ``name``,     "演算子の名前",       文字列
-   ``group``,    "演算子の識別番号",   整数
-   ``sites``,    "サイト番号",         整数 or 整数のリスト
-   ``dim``,      "演算子の次元",       整数
-   ``elements``, "演算子の非ゼロ要素", 文字列
+   ``name``,     "演算子の名前",               文字列
+   ``group``,    "演算子の識別番号",           整数
+   ``sites``,    "サイト番号",                 整数 or 整数のリスト
+   ``dim``,      "演算子の次元",               整数
+   ``elements``, "演算子の非ゼロ要素",         文字列
+   ``coeff``,    "演算子にかかる係数（実部）", 実数
+   ``coeff_im``, "演算子にかかる係数（虚部）", 実数
 
 ``name`` は演算子の名前です。
 
@@ -36,6 +38,9 @@
 - 最初の2つはそれぞれ演算子が作用する前と後の状態番号を示します。
 - あとの2つはそれぞれ演算子の要素の実部と虚部を示します。
 
+``coeff``, ``coeff_im`` は演算子にかかる係数の実部と虚部です。
+省略した場合はそれぞれ 1.0, 0.0 になります。
+
 例
 ....
 
@@ -95,6 +100,8 @@ S=1/2 のSz 演算子
    ``dim``,      "演算子の次元",             整数のリスト
    ``elements``, "演算子の非ゼロ要素",       文字列
    ``ops``,      "onesite 演算子の識別番号", 整数のリスト
+   ``coeff``,    "演算子にかかる係数（実部）", 実数
+   ``coeff_im``, "演算子にかかる係数（虚部）", 実数
 
 ``name`` は演算子の名前です。
 
@@ -124,6 +131,9 @@ S=1/2 のSz 演算子
 
 ``elements`` と ``ops`` を同時に定義した場合にはエラー終了します。
 
+``coeff``, ``coeff_im`` は演算子にかかる係数の実部と虚部です。
+省略した場合はそれぞれ 1.0, 0.0 になります。
+
 例
 ....
 ここでは具体例として、``Lsub=[2,2]`` の正方格子 S=1/2 ハイゼンベルグ模型のボンドハミルトニアンのエネルギーを求めるため、
@@ -197,6 +207,8 @@ S=1/2 のSz 演算子
    ``group``,      "演算子の識別番号",         整数
    ``multisites``, "サイトの組み合わせ",       文字列
    ``ops``,        "onesite 演算子の識別番号", 整数のリスト
+   ``coeff``,    "演算子にかかる係数（実部）", 実数
+   ``coeff_im``, "演算子にかかる係数（虚部）", 実数
 
 ``name`` は演算子の名前です。
 
@@ -215,3 +227,5 @@ S=1/2 のSz 演算子
 例えば ``observable.onesite`` の ``group=0`` として :math:`S^z` を定義していた場合には、
 ``ops = [0,0,0]`` として :math:`S^z_iS^z_jS^z_k` を表現できます。
 
+``coeff``, ``coeff_im`` は演算子にかかる係数の実部と虚部です。
+省略した場合はそれぞれ 1.0, 0.0 になります。

src/iTPS/load_toml.cpp

@@ -78,6 +78,13 @@ struct bittype<long long int> {
   typedef int64_t type;
 };
 
+template <class T>
+T reduce(double re, double im);
+template <>
+std::complex<double> reduce<std::complex<double>>(double re, double im) { return std::complex<double>(re, im); }
+template <>
+double reduce<double>(double re, double im) { return re; }
+
 }  // end of namespace detail
 
 template <class T>
@@ -516,9 +523,17 @@ Operators<tensor> load_operator(decltype(cpptoml::parse_file("")) param,
 
   auto group = find<int>(param, "group");
   auto name = find_or(param, "name", std::string(""));
+  auto coeff_re = find_or(param, "coeff", 1.0);
+  auto coeff_im = find_or(param, "coeff_im", 0.0);
+  typename tensor::value_type coeff = detail::reduce<typename tensor::value_type>(coeff_re, coeff_im);
+  auto is_real = std::is_same<typename tensor::value_type, double>::value;
+  if(is_real && coeff_im != 0.0){
+    std::stringstream ss;
+    ss << "parameter.general.is_real is true but coeff_im is not zero in a section " << tablename;
+    throw tenes::input_error(ss.str());
+  }
 
   std::vector<Operator<tensor>> ret;
-
   if (nbody == 1) {
     auto site_int = param->get_as<int>("sites");
     auto site_arr = param->get_array_of<int64_t>("sites");
@@ -536,18 +551,18 @@ Operators<tensor> load_operator(decltype(cpptoml::parse_file("")) param,
       throw input_error(detail::msg_cannot_find("sites", tablename));
     }
     for (int s : sites) {
-      ret.emplace_back(name, group, s, A);
+      ret.emplace_back(name, group, s, A, coeff);
     }
   } else if (nbody == 2) {
     auto bonds_str = find<std::string>(param, "bonds");
     auto bonds = read_bonds(bonds_str);
     for (auto bond : bonds) {
       if (elements) {
         ret.emplace_back(name, group, std::get<0>(bond), std::get<1>(bond),
-                         std::get<2>(bond), A);
+                         std::get<2>(bond), A, coeff);
       } else {
         ret.emplace_back(name, group, std::get<0>(bond), std::get<1>(bond),
-                         std::get<2>(bond), op_ind);
+                         std::get<2>(bond), op_ind, coeff);
       }
     }
   } else {
@@ -556,10 +571,10 @@ Operators<tensor> load_operator(decltype(cpptoml::parse_file("")) param,
     for (auto ms : mss) {
       if (elements) {
         ret.emplace_back(name, group, std::get<0>(ms), std::get<1>(ms),
-                         std::get<2>(ms), A);
+                         std::get<2>(ms), A, coeff);
       } else {
         ret.emplace_back(name, group, std::get<0>(ms), std::get<1>(ms),
-                         std::get<2>(ms), op_ind);
+                         std::get<2>(ms), op_ind, coeff);
       }
     }
   }

src/iTPS/main.cpp

@@ -48,23 +48,24 @@ tenes::Operators<mptensor::Tensor<tenes::mptensor_matrix_type, double>> to_real(
   }
   MPI_Comm comm = ops[0].op.get_comm();
   for (auto const& op : ops) {
+    double coeff = std::real(op.coeff);
     if (op.is_onesite()) {
       mptensor::Tensor<tenes::mptensor_matrix_type, double> A(comm,
                                                               op.op.shape());
       for (size_t lindex = 0; lindex < op.op.local_size(); ++lindex) {
         A[lindex] = op.op[lindex].real();
       }
-      ret.emplace_back(op.name, op.group, op.source_site, A);
+      ret.emplace_back(op.name, op.group, op.source_site, A, coeff);
     } else if (op.ops_indices.empty()) {
       mptensor::Tensor<tenes::mptensor_matrix_type, double> A(comm,
                                                               op.op.shape());
       for (size_t lindex = 0; lindex < op.op.local_size(); ++lindex) {
         A[lindex] = op.op[lindex].real();
       }
-      ret.emplace_back(op.name, op.group, op.source_site, op.dx, op.dy, A);
+      ret.emplace_back(op.name, op.group, op.source_site, op.dx, op.dy, A, coeff);
     } else {
       ret.emplace_back(op.name, op.group, op.source_site, op.dx, op.dy,
-                       op.ops_indices);
+                       op.ops_indices, coeff);
     }
   }
   return ret;

src/iTPS/multisite_obs.cpp

@@ -202,7 +202,7 @@ auto iTPS<ptensor>::measure_multisite()
           core::Contract(C_, eTt_, eTr_, eTb_, eTl_, Tn_, op_, is_density, is_meanfield);
       value += localvalue;
     }
-    ret[op.group][{op.source_site, op.dx, op.dy}] = value / norm;
+    ret[op.group][{op.source_site, op.dx, op.dy}] = op.coeff * value / norm;
   }
   // ret.push_back(norms);
 

src/iTPS/onesite_obs.cpp

@@ -66,7 +66,7 @@ auto iTPS<tensor>::measure_onesite()
     for (auto const &op : onesite_operators) {
       const int i = op.source_site;
       const auto val = core::Contract_one_site_iTPS_MF(Tn_[i], op.op);
-      local_obs[op.group][i] = val / norm[i];
+      local_obs[op.group][i] = op.coeff * val / norm[i];
     }
   } else {  // CTM
     if (is_density) {
@@ -80,7 +80,7 @@ auto iTPS<tensor>::measure_onesite()
         const auto val = core::Contract_one_site_density_CTM(
             C1[i], C2[i], C3[i], C4[i], eTt[i], eTr[i], eTb[i], eTl[i], Tn[i],
             op.op);
-        local_obs[op.group][i] = val / norm[i];
+        local_obs[op.group][i] = op.coeff * val / norm[i];
       }
     } else {
       for (int i = 0; i < N_UNIT; ++i) {
@@ -93,7 +93,7 @@ auto iTPS<tensor>::measure_onesite()
         const auto val = core::Contract_one_site_iTPS_CTM(
             C1[i], C2[i], C3[i], C4[i], eTt[i], eTr[i], eTb[i], eTl[i], Tn[i],
             op.op);
-        local_obs[op.group][i] = val / norm[i];
+        local_obs[op.group][i] = op.coeff * val / norm[i];
       }
     }
   }
@@ -240,7 +240,7 @@ auto iTPS<tensor>::measure_onesite_density()
     const auto val = core::Contract_one_site_density_CTM(
         C1[i], C2[i], C3[i], C4[i], eTt[i], eTr[i], eTb[i], eTl[i], Tn[i],
         op.op);
-    local_obs[op.group][i] = val / norm[i];
+    local_obs[op.group][i] = op.coeff * val / norm[i];
   }
   //  }
   double norm_real_min = 1e100;

src/iTPS/twosite_obs.cpp

@@ -258,7 +258,7 @@ auto iTPS<ptensor>::measure_twosite()
       //         : core::Contract_CTM(C_, eTt_, eTr_, eTb_, eTl_, Tn_, op_);
       value += localvalue;
     }
-    ret[op.group][{op.source_site, op.dx[0], op.dy[0]}] = value / norm;
+    ret[op.group][{op.source_site, op.dx[0], op.dy[0]}] = op.coeff * value / norm;
   }
   ret.push_back(norms);
 
@@ -583,7 +583,7 @@ auto iTPS<ptensor>::measure_twosite_density()
         }
       }
     }
-    ret[op.group][{op.source_site, op.dx[0], op.dy[0]}] = value / norm;
+    ret[op.group][{op.source_site, op.dx[0], op.dy[0]}] = op.coeff * value / norm;
   }
   ret.push_back(norms);
 

src/operator.hpp

@@ -19,18 +19,22 @@
 
 #include <vector>
 #include <string>
+#include <iostream>
 
 namespace tenes {
 
 template <class tensor>
 struct Operator {
+  using value_type = typename tensor::value_type;
+
   std::string name;
   int group;
   int source_site;
   std::vector<int> dx;
   std::vector<int> dy;
   tensor op;
   std::vector<int> ops_indices;
+  typename tensor::value_type coeff;
 
   /*!
    * @brief Constructor for a one-site operator.
@@ -39,9 +43,10 @@ struct Operator {
    * @param[in] group Group of the operator.
    * @param[in] site Site of the operator.
    * @param[in] op Operator tensor.
+   * @param[in] coeff coefficient of operator (default: 1.0).
    */
-  Operator(std::string const &name, int group, int site, tensor const &op)
-      : name(name), group(group), source_site(site), dx(0), dy(0), op(op) {
+  Operator(std::string const &name, int group, int site, tensor const &op, value_type coeff=static_cast<value_type>(1.0))
+      : name(name), group(group), source_site(site), dx(0), dy(0), op(op), coeff(coeff) {
     if (op.rank() != 2) {
       throw std::runtime_error("Operator tensor must be rank 2.");
     }
@@ -56,15 +61,17 @@ struct Operator {
    * @param[in] dx X displacement of the other site.
    * @param[in] dy Y displacement of the other site.
    * @param[in] op Operator tensor.
+   * @param[in] coeff coefficient of operator (default: 1.0).
    */
   Operator(std::string const &name, int group, int source_site, int dx, int dy,
-           tensor const &op)
+           tensor const &op, value_type coeff=static_cast<value_type>(1.0))
       : name(name),
         group(group),
         source_site(source_site),
         dx(1, dx),
         dy(1, dy),
-        op(op) {
+        op(op),
+        coeff(coeff) {
     if (op.rank() != 4) {
       throw std::runtime_error("Operator tensor must be rank 4.");
     }
@@ -80,15 +87,17 @@ struct Operator {
    * @param[in] dx X displacement of the other site.
    * @param[in] dy Y displacement of the other site.
    * @param[in] ops_indices Onesite operator indices.
+   * @param[in] coeff coefficient of operator (default: 1.0).
    */
   Operator(std::string const &name, int group, int source_site, int dx, int dy,
-           std::vector<int> const &ops_indices)
+           std::vector<int> const &ops_indices, value_type coeff=static_cast<value_type>(1.0))
       : name(name),
         group(group),
         source_site(source_site),
         dx(1, dx),
         dy(1, dy),
-        ops_indices(ops_indices) {
+        ops_indices(ops_indices),
+        coeff(coeff){
     if (ops_indices.size() != 2) {
       throw std::runtime_error(
           "Operator must be a product of two one-site operators.");
@@ -104,16 +113,18 @@ struct Operator {
    * @param[in] dx X displacement of the other sites.
    * @param[in] dy Y displacement of the other sites.
    * @param[in] op Operator tensor.
+   * @param[in] coeff coefficient of operator (default: 1.0).
    */
   Operator(std::string const &name, int group, int source_site,
            std::vector<int> const &dx, std::vector<int> const &dy,
-           tensor const &op)
+           tensor const &op, value_type coeff=static_cast<value_type>(1.0))
       : name(name),
         group(group),
         source_site(source_site),
         dx(dx),
         dy(dy),
-        op(op) {
+        op(op),
+        coeff(coeff) {
     if (dx.size() != dy.size()) {
       throw std::runtime_error("dx and dy must have the same size.");
     }
@@ -132,16 +143,18 @@ struct Operator {
    * @param[in] dx X displacement of the other sites.
    * @param[in] dy Y displacement of the other sites.
    * @param[in] ops_indices Onesite operator indices.
+   * @param[in] coeff coefficient of operator (default: 1.0).
    */
   Operator(std::string const &name, int group, int source_site,
            std::vector<int> const &dx, std::vector<int> const &dy,
-           std::vector<int> const &ops_indices)
+           std::vector<int> const &ops_indices, value_type coeff=static_cast<value_type>(1.0))
       : name(name),
         group(group),
         source_site(source_site),
         dx(dx),
         dy(dy),
-        ops_indices(ops_indices) {
+        ops_indices(ops_indices),
+        coeff(coeff) {
     if (dx.size() != dy.size()) {
       throw std::runtime_error("dx and dy must have the same size.");
     }