Added cut_off_observable

Author: haggaila

lindbladmpo/LindbladMPOSolver.py

@@ -955,7 +955,7 @@ def verify_parameters(
                         "Error 410: " + key + " must be a non-negative integer\n"
                     )
                     continue
-            elif (key == "cut_off") or (key == "cut_off_rho"):
+            elif (key == "cut_off") or (key == "cut_off_rho") or (key == "cut_off_observable"):
                 if not LindbladMPOSolver.is_float(parameters[key]):
                     check_msg += "Error 420: " + key + " is not a float\n"
                     continue

lindbladmpo/examples/simulation_building/LindbladMatrixSolver.py

@@ -394,6 +394,7 @@ def solve(self):
 
             if parameters.get("b_save_final_state", False):  # Write final state
                 np.save(self.s_output_path + ".state", sol.y[-1])
+            cut_off_observable = parameters.get("cut_off_observable", .0)
 
             file_1q = open(self.s_output_path + ".obs-1q.dat", "w")
             file_2q = open(self.s_output_path + ".obs-2q.dat", "w")
@@ -419,6 +420,8 @@ def solve(self):
                             f"Warning: imaginary component {val.imag} in observable "
                             f"{key[0].upper()}, qubit {key[1] + 1}.\n"
                         )
+                    if cut_off_observable and (abs(val) < cut_off_observable):
+                        val = .0
                     file_1q.write(f"{t}\t{key[0].upper()}\t{key[1] + 1}\t{val.real}\n")
                 file_1q.write("\n")
                 file_1q.flush()
@@ -431,6 +434,8 @@ def solve(self):
                             f"Warning: imaginary component {val.imag} in observable "
                             f"{key[0].upper()}, qubits ({key[1] + 1}, {key[2] + 1}).\n"
                         )
+                    if cut_off_observable and (abs(val) < cut_off_observable):
+                        val = .0
                     file_2q.write(
                         f"{t}\t{key[0].upper()}\t{key[1] + 1}\t{key[2] + 1}\t{val.real}\n"
                     )
@@ -446,6 +451,8 @@ def solve(self):
                             f"{key[0].upper()}, qubits ({key[1] + 1}, {key[2] + 1}, "
                             f"{key[3] + 1}).\n"
                         )
+                    if cut_off_observable and (abs(val) < cut_off_observable):
+                        val = .0
                     file_3q.write(
                         f"{t}\t{key[0].upper()}\t{key[1] + 1}\t{key[2] + 1}\t{key[3] + 1}\t{val.real}\n"
                     )
@@ -459,6 +466,8 @@ def solve(self):
                             f"Warning: imaginary component {val.imag} in observable "
                             f"{obs_cu_name}.\n"
                         )
+                    if cut_off_observable and (abs(val) < cut_off_observable):
+                        val = .0
                     file_cu.write(f"{t}\t{obs_cu_name}\t{val.real}\n")
                     # print(f"{t}\t{obs_cu_name}\t{val.real}")
                 file_cu.write("\n")

src/SimulationParameters.h

@@ -24,8 +24,11 @@ class SimulationParameters : public Parameters
         operator[]("trotter_order") = "4"; // Possible choices are 2, 3, 4. 3 or 4 are recommended.
         operator[]("max_dim_rho") = "400"; // Maximum bond dimension for density matrices
         operator[]("cut_off_rho") =
-            "1e-16"; // Maximum truncation error for density matrices. The actual
-                     // truncation is done using the most severe condition between cut_off_rho and max_dim_rho.
+            "1e-16"; // Maximum truncation error for density matrices. The actual truncation is
+                     // done using the most severe condition between cut_off_rho and max_dim_rho.
+        operator[]("cut_off_observable") = "0"; // If a nonzero value is given, it is the maximum
+                     // (absolute) value of an observable to be written to output files -
+                     // smaller values are truncated to 0.
 
         operator[]("b_force_rho_trace") = "1";         // Whether to force the density matrix trace to 1,
                                                        // by substituting rho /= trace{rho} at every time step and every
@@ -41,7 +44,7 @@ class SimulationParameters : public Parameters
             "0"; // If nonzero, after reading rho from
                  // a saved file, perform a re-gauging/compression using ITensor's method orthogonalize().
         operator[]("b_quiet") = "0"; // If nonzero, after initialization of the simulation
-                                     // avoid the console output at every time step (but write it the log file).
+                 // avoid the console output at every time step (but write it to the log file).
 
         operator[]("init_product_state") =
             ""; // Initialize a product state. Can be specified

src/lindbladmpo.cc

@@ -130,14 +130,14 @@ int main(int argc, char *argv[])
     }
 
     C.ConstructIdentity(); // Construct the density matrix corresponding to infinite temperature
-
+    const double cut_off_rho = param.val("cut_off_rho");
     // Note on the option below: if we  "Normalize=true" (default for fitapplyMPO)
     // we would normalize rho such that Tr[rho^2]=1, (norm of the MPS)
     // which is of course not appropriate for a density matrix.
     auto argsRho = Args();
     argsRho.add("Normalize", false);
     argsRho.add("MaxDim", param.longval("max_dim_rho"));
-    argsRho.add("Cutoff", param.val("cut_off_rho"));
+    argsRho.add("Cutoff", cut_off_rho);
 
     vector<long> init_graph_state = param.longvec("init_graph_state");
     vector<long> init_cz_gates = param.longvec("init_cz_gates");
@@ -680,7 +680,7 @@ int main(int argc, char *argv[])
 
         cout2 << "C.rho.orthogonalize...";
         cout2.flush();
-        C.rho.orthogonalize(Args("Cutoff", param.val("cut_off_rho"), "MaxDim", param.longval("max_dim_rho")));
+        C.rho.orthogonalize(Args("Cutoff", cut_off_rho, "MaxDim", param.longval("max_dim_rho")));
         cout2 << "done.\n";
         cout2 << "New max bond dimension of rho:" << maxLinkDim(C.rho) << "\n";
         cout2.flush();
@@ -912,7 +912,7 @@ int main(int argc, char *argv[])
                 if (std::abs(z) < _2_N)
                     cout2 << "\t\tNote: " << "Tr{rho} = " << z << " encountered during collapse projectors, "
                           << "this is smaller than 2^(-N)!";
-                // C.rho /= z;
+                 C.rho /= z;
             }
             if (i_coll == 0)
                 rho_c = MPS(C.rho);
@@ -928,16 +928,17 @@ int main(int argc, char *argv[])
             cout2 << "\t\tNote: this is smaller than 2^(-N)!"
                   << "\n";
         // TODO: This is a somewhat arbitrary threshold for the warning.
-        // C.rho /= z;
+         C.rho /= z;
         auto t_collapse_end = steady_clock::now();
         duration_ms = duration_cast<milliseconds>(t_collapse_end - t_collapse_start);
         cout2 << "Collapse evaluation terminated. Duration: " << duration_ms.count() / 1000. << "s\n";
     }
 
     char buf[100];
+    const bool b_quiet = param.boolval("b_quiet");
+    const double cut_off_observable = param.val("cut_off_observable");
     const long force_rho_hermitian_step = param.longval("force_rho_hermitian_step");
     const long force_rho_hermitian_gates = param.longval("force_rho_hermitian_gates");
-    const bool b_quiet = param.boolval("b_quiet");
     cout2.quiet(b_quiet);
     double t = t_0;
     for (int n = 0; n <= n_steps; n++)
@@ -1093,6 +1094,8 @@ int main(int argc, char *argv[])
                                 cout2 << "\tWarning: <S^" << s << "(" << i << ")> = " << expectation_value
                                       << "; it should be real, but has an imaginary part > " << IMAGINARY_THRESHOLD
                                       << ".\n";
+                            if (cut_off_observable && (abs(expectation_value) < cut_off_observable))
+                                expectation_value = .0;
                             file_1q << t << "\t" << char(toupper(s[0])) << "\t" << i << "\t" << expectation_value.real()
                                     << endl;
                             count++;
@@ -1129,6 +1132,8 @@ int main(int argc, char *argv[])
                                       << ")> = " << expectation_value
                                       << "; it should be real, but has an imaginary part > " << IMAGINARY_THRESHOLD
                                       << ".\n";
+                            if (cut_off_observable && (abs(expectation_value) < cut_off_observable))
+                                expectation_value = .0;
                             file_2q << t << "\t" << char(toupper(s[0])) << char(toupper(s[1])) << "\t" << i << "\t" << j
                                     << "\t" << expectation_value.real() << endl;
                             count++;
@@ -1165,6 +1170,8 @@ int main(int argc, char *argv[])
                                       << k << ")" << expectation_value
                                       << "; it should be real, but has an imaginary part > " << IMAGINARY_THRESHOLD
                                       << ".\n";
+                            if (cut_off_observable && (abs(expectation_value) < cut_off_observable))
+                                expectation_value = .0;
                             file_3q << t << "\t" << char(toupper(s[0])) << char(toupper(s[1])) << char(toupper(s[2]))
                                     << "\t" << i << "\t" << j << "\t" << k << "\t" << expectation_value.real() << endl;
                             count++;
@@ -1192,6 +1199,8 @@ int main(int argc, char *argv[])
                         for (MPS &proj : ProjectorList)
                         {
                             Cplx op_val = innerC(proj, C.rho);
+                            if (cut_off_observable && (abs(op_val) < cut_off_observable))
+                                op_val = .0;
                             file_custom << t << " \t" << ProjectorNames[c] << "\t" << op_val.real() << endl;
                             c++;
                         }
@@ -1203,6 +1212,8 @@ int main(int argc, char *argv[])
                         for (string &s_op_obs : OperatorObsNames)
                         {
                             Cplx op_val = C.Expect(OperatorObs[c], OperatorObsQubits[c]);
+                            if (cut_off_observable && (abs(op_val) < cut_off_observable))
+                                op_val = .0;
                             // cout2 << "\nCalculating: " << OperatorObs[c] << " " << OperatorObsQubits[c];
                             file_custom << t << " \t" << s_op_obs << "\t" << op_val.real() << endl;
                             c++;