Implemented generator of spin Hamiltonians; Added new ctor for TensorOperator

Signed-off-by: Dmitry I. Lyakh <quant4me@gmail.com>

src/exatn/quantum.cpp

 /** ExaTN: Quantum computing related
-REVISION: 2021/09/25
+REVISION: 2021/10/26
 
 Copyright (C) 2018-2021 Dmitry I. Lyakh (Liakh)
 Copyright (C) 2018-2021 Oak Ridge National Laboratory (UT-Battelle) **/
@@ -179,6 +179,41 @@ std::shared_ptr<exatn::numerics::TensorOperator> readSpinHamiltonian(const std::
  return tens_oper;
 }
 
+
+std::shared_ptr<exatn::numerics::TensorOperator> generateSpinHamiltonian(const std::string & operator_name,
+                                                  std::function<PauliProduct ()> hamiltonian_generator,
+                                                  TensorElementType precision)
+{
+ auto hamiltonian = exatn::makeSharedTensorOperator(operator_name);
+ while(true){
+  const auto pauli_prod = hamiltonian_generator();
+  if(pauli_prod.product.size() == 0 && pauli_prod.coefficient == std::complex<double>{0.0,0.0}) break;
+  std::string paulis = "[";
+  bool not_first = false;
+  for(const auto & pauli: pauli_prod.product){
+   if(not_first) paulis += " ";
+   if(pauli.pauli_gate == Gate::gate_I){
+    paulis += "I";
+   }else if(pauli.pauli_gate == Gate::gate_X){
+    paulis += "X";
+   }else if(pauli.pauli_gate == Gate::gate_Y){
+    paulis += "Y";
+   }else if(pauli.pauli_gate == Gate::gate_Z){
+    paulis += "Z";
+   }else{
+    std::cout << "#ERROR(exatn::quantum::generateSpinHamiltonian): Invalid gate returned by the generator: "
+              << static_cast<int>(pauli.pauli_gate) << std::endl;
+    assert(false);
+   }
+   paulis += std::to_string(pauli.qubit);
+   not_first = true;
+  }
+  paulis += "]";
+  auto success = appendPauliComponent(*hamiltonian,paulis,pauli_prod.coefficient,precision); assert(success);
+ }
+ return hamiltonian;
+}
+
 } //namespace quantum
 
 } //namespace exatn

src/exatn/quantum.hpp

 /** ExaTN: Quantum computing related
-REVISION: 2021/10/01
+REVISION: 2021/10/26
 
 Copyright (C) 2018-2021 Dmitry I. Lyakh (Liakh)
 Copyright (C) 2018-2021 Oak Ridge National Laboratory (UT-Battelle) **/
 
 /** Rationale:
  a) Provides utilities related to quantum circuit simulations, like quantum gates,
-    Pauli matrix based Hamiltonian reading, etc.
+    Pauli matrix based Hamiltonian reading or generation, etc.
  b) Normal gate action translates to the following tensor notation:
      q(j0) * G(j0|i0) --> v(i0),
      q(j1,j0) * G(j1,j0|i1,i0) --> v(i1,i0), etc,
@@ -15,9 +15,10 @@ Copyright (C) 2018-2021 Oak Ridge National Laboratory (UT-Battelle) **/
     storage of matrix G(j1,j0|i1,i0), to match the textbook
     definitions of quantum gates. Note that if G(j1,j0|i1,i0)
     is a controlled 2-body gate, the control (senior) indices
-    are i0 and j0. A convenient way to remember this is to use
-    the standard bra-ket convention such that we will have:
-     <v(i1,i0)| = <q(j1,j0)| * CX(j1,j0|i1,i0)
+    are i0 and j0. A convenient way to remember this is to
+    adhere to the bra convention such that we will have:
+     <v(i1,i0)| = <q(j1,j0)| * CX(j1,j0|i1,i0), where
+    the CX gate is applied to a 2-qubit register q.
 **/
 
 #ifndef EXATN_QUANTUM_HPP_
@@ -26,7 +27,9 @@ Copyright (C) 2018-2021 Oak Ridge National Laboratory (UT-Battelle) **/
 #include "exatn_numerics.hpp"
 #include "tensor_operator.hpp"
 #include "tensor_symbol.hpp"
+#include "tensor_range.hpp"
 
+#include <functional>
 #include <vector>
 #include <complex>
 #include <string>
@@ -53,8 +56,22 @@ enum class Gate{
  gate_ISWAP
 };
 
+//Pauli gate acting on a specific qubit:
+struct PauliMap {
+ Gate pauli_gate;   //Pauli gate (I,X,Y,Z)
+ std::size_t qubit; //flattened qubit id
+};
+
+//Product of Pauli gates:
+struct PauliProduct {
+ std::vector<PauliMap> product;              //product of Pauli maps
+ std::complex<double> coefficient {0.0,0.0}; //linear combination coefficient
+};
+
+
 /** Returns the data initialization vector for a specific quantum gate
-    that can subsequently be used for initializing its tensor. **/
+    that can subsequently be used for initializing its tensor.
+**/
 std::vector<std::complex<double>> getGateData(const Gate gate_name,
                                               std::initializer_list<double> angles = {});
 
@@ -67,6 +84,23 @@ std::shared_ptr<exatn::numerics::TensorOperator> readSpinHamiltonian(const std::
                                                                      const std::string & filename,
                                                                      TensorElementType precision = TensorElementType::COMPLEX64,
                                                                      const std::string & format = "OpenFermion");
+
+/** Generates a grid-based spin Hamiltonian of the form:
+     H = Sum{i1,...,iK} [H(i1,...,iK) * P(i1) * ... * P(iK)],
+    where P(i) is a Pauli matrix acting on spin site i. In general,
+    the index i enumerating the spin sites is a flattened D-dimensional
+    multi-index, where D is the dimensionality of the spin grid.
+    In general, the Sum can contain Pauli matrix products of varying
+    length, K. The constituting individual Pauli matrix products are
+    generated by the provided Lambda generator <hamiltonian_generator>
+    which is supposed to return an individual PauliProduct upon each
+    invocation (in their respective order). The end of the generated
+    sequence is signalled by an empty PauliProduct with zero coefficient.
+**/
+std::shared_ptr<exatn::numerics::TensorOperator> generateSpinHamiltonian(const std::string & operator_name,
+                                                  std::function<PauliProduct ()> hamiltonian_generator,
+                                                  TensorElementType precision = TensorElementType::COMPLEX64);
+
 } //namespace quantum
 
 } //namespace exatn

src/exatn/tests/NumServerTester.cpp

 #include <gtest/gtest.h>
 
 #include "exatn.hpp"
-#include "talshxx.hpp"
 #include "quantum.hpp"
+#include "talshxx.hpp"
 
 #ifdef MPI_ENABLED
 #include "mpi.h"
@@ -18,7 +18,7 @@
 #include "errors.hpp"
 
 //Test activation:
-#define EXATN_TEST0
+/*#define EXATN_TEST0
 #define EXATN_TEST1
 #define EXATN_TEST2
 #define EXATN_TEST3
@@ -47,8 +47,9 @@
 #define EXATN_TEST26
 //#define EXATN_TEST27 //requires input file from source
 //#define EXATN_TEST28 //requires input file from source
-#define EXATN_TEST29
+#define EXATN_TEST29*/
 #define EXATN_TEST30
+//#define EXATN_TEST31
 
 
 #ifdef EXATN_TEST0
@@ -1549,7 +1550,7 @@ TEST(NumServerTester, IsingTNO)
 
  bool success = true;
 
- exatn::resetLoggingLevel(2,2); //debug
+ //exatn::resetLoggingLevel(2,2); //debug
 
  //Define Ising Hamiltonian constants:
  constexpr std::complex<double> ZERO{ 0.0, 0.0};
@@ -3387,6 +3388,117 @@ TEST(NumServerTester, TensorOperatorReconstruction) {
 #endif
 
 #ifdef EXATN_TEST30
+TEST(NumServerTester, SpinHamiltonians) {
+ using exatn::TensorShape;
+ using exatn::TensorSignature;
+ using exatn::Tensor;
+ using exatn::TensorComposite;
+ using exatn::TensorNetwork;
+ using exatn::TensorExpansion;
+ using exatn::TensorOperator;
+ using exatn::TensorElementType;
+ using exatn::TensorRange;
+ using exatn::quantum::Gate;
+ using exatn::quantum::PauliMap;
+ using exatn::quantum::PauliProduct;
+
+ const auto TENS_ELEM_TYPE = TensorElementType::COMPLEX64;
+
+ const std::complex<double> j_param {-1.0,0.0};
+ const std::complex<double> h_param {-0.1,0.0};
+ const int num_spin_sites = 4;
+ const int bond_dim_lim = 4;
+ const int max_bond_dim = std::min(static_cast<int>(std::pow(2,num_spin_sites/2)),bond_dim_lim);
+ const int arity = 2;
+ const std::string tn_type = "TTN"; //MPS or TTN
+
+ //exatn::resetLoggingLevel(2,2); //debug
+
+ bool success = true;
+
+ //Define the 1D-Ising Hamiltonian generator:
+ TensorRange spin_sites({num_spin_sites});
+ auto ising_generator = [j_param,
+                         spin_sites,
+                         num_sites = spin_sites.localVolume(),
+                         finished = false] () mutable -> PauliProduct {
+  assert(num_sites > 1);
+  PauliProduct pauli_product;
+  if(!finished){
+   const auto spin_site = spin_sites.localOffset();
+   pauli_product.product.emplace_back(PauliMap{Gate::gate_Z,spin_site});
+   pauli_product.product.emplace_back(PauliMap{Gate::gate_Z,spin_site+1});
+   pauli_product.coefficient = j_param;
+   if(spin_site < (num_sites - 2)){
+    spin_sites.next();
+   }else{
+    finished = true;
+   }
+  }
+  return pauli_product;
+ };
+
+ //Construct the 1D-Ising Hamiltonian using the generator:
+ auto ising_hamiltonian = exatn::quantum::generateSpinHamiltonian("IsingHamiltonian",
+                                                                  ising_generator,
+                                                                  TENS_ELEM_TYPE);
+ //ising_hamiltonian->printIt(); //debug
+
+ //Configure the tensor network builder:
+ auto tn_builder = exatn::getTensorNetworkBuilder(tn_type);
+ if(tn_type == "MPS"){
+  success = tn_builder->setParameter("max_bond_dim",max_bond_dim); assert(success);
+ }else if(tn_type == "TTN"){
+  success = tn_builder->setParameter("max_bond_dim",max_bond_dim); assert(success);
+  success = tn_builder->setParameter("arity",arity); assert(success);
+ }else{
+  assert(false);
+ }
+
+ //Build tensor network vectors:
+ auto ket_tensor = exatn::makeSharedTensor("TensorSpace",std::vector<int>(num_spin_sites,2));
+ auto vec_net = exatn::makeSharedTensorNetwork("VectorNet",ket_tensor,*tn_builder,false);
+ vec_net->markOptimizableAllTensors();
+ //vec_net->printIt(); //debug
+ auto vec_tns = exatn::makeSharedTensorExpansion("VectorTNS",vec_net,std::complex<double>{1.0,0.0});
+ auto rhs_net = exatn::makeSharedTensorNetwork("RightHandSideNet",ket_tensor,*tn_builder,false);
+ auto rhs_tns = exatn::makeSharedTensorExpansion("RightHandSideTNS",rhs_net,std::complex<double>{1.0,0.0});
+
+ //Numerical processing:
+ {
+  //Create and initialize tensor network vector tensors:
+  std::cout << "Creating and initializing tensor network vector tensors ... ";
+  success = exatn::createTensorsSync(*vec_net,TENS_ELEM_TYPE); assert(success);
+  success = exatn::initTensorsRndSync(*vec_net); assert(success);
+  success = exatn::createTensorsSync(*rhs_net,TENS_ELEM_TYPE); assert(success);
+  success = exatn::initTensorsRndSync(*rhs_net); assert(success);
+  std::cout << "Ok" << std::endl;
+
+  //Ground state search for the original Hamiltonian:
+  std::cout << "Ground state search for the original Hamiltonian:" << std::endl;
+  exatn::TensorNetworkOptimizer::resetDebugLevel(1,0);
+  exatn::TensorNetworkOptimizer optimizer(ising_hamiltonian,vec_tns,1e-5);
+  success = exatn::sync(); assert(success);
+  bool converged = optimizer.optimize();
+  success = exatn::sync(); assert(success);
+  if(converged){
+   std::cout << "Search succeeded: ";
+  }else{
+   std::cout << "Search failed!" << std::endl;
+   assert(false);
+  }
+  const auto expect_val = optimizer.getExpectationValue();
+  std::cout << "Expectation value = " << expect_val << std::endl;
+ }
+
+ //Synchronize:
+ success = exatn::sync(); assert(success);
+ exatn::resetLoggingLevel(0,0);
+ //Grab a beer!
+}
+#endif
+
+#ifdef EXATN_TEST31
 TEST(NumServerTester, TensorComposite) {
  using exatn::TensorShape;
  using exatn::TensorSignature;

src/numerics/tensor_operator.cpp

 /** ExaTN::Numerics: Tensor operator
-REVISION: 2021/10/22
+REVISION: 2021/10/26
 
 Copyright (C) 2018-2021 Dmitry I. Lyakh (Liakh)
 Copyright (C) 2018-2021 Oak Ridge National Laboratory (UT-Battelle) **/
@@ -11,6 +11,37 @@ namespace exatn{
 
 namespace numerics{
 
+TensorOperator::TensorOperator(const std::string & name,
+                               std::shared_ptr<TensorNetwork> network,
+                               const std::vector<std::pair<unsigned int, unsigned int>> & ket_pairing,
+                               const std::vector<std::pair<unsigned int, unsigned int>> & bra_pairing,
+                               const std::complex<double> coefficient):
+ name_(name)
+{
+ auto success = appendComponent(network,ket_pairing,bra_pairing,coefficient);
+ assert(success);
+}
+
+
+TensorOperator::TensorOperator(const std::string & name,
+                               std::shared_ptr<TensorNetwork> ket_network,
+                               std::shared_ptr<TensorNetwork> bra_network,
+                               const std::vector<std::pair<unsigned int, unsigned int>> & ket_pairing,
+                               const std::vector<std::pair<unsigned int, unsigned int>> & bra_pairing,
+                               const std::complex<double> coefficient):
+ name_(name)
+{
+ auto shifted_bra_pairing = bra_pairing;
+ const auto shift = ket_network->getRank();
+ for(auto & pairing: shifted_bra_pairing) pairing.second += shift;
+ auto combined_network = makeSharedTensorNetwork(*ket_network,true,ket_network->getName());
+ auto success = combined_network->appendTensorNetwork(TensorNetwork(*bra_network,true,bra_network->getName()),{});
+ assert(success);
+ success = appendComponent(combined_network,ket_pairing,shifted_bra_pairing,coefficient);
+ assert(success);
+}
+
+
 bool TensorOperator::appendComponent(std::shared_ptr<TensorNetwork> network, //in: tensor network (or single tensor as a tensor network)
      const std::vector<std::pair<unsigned int, unsigned int>> & ket_pairing, //in: ket pairing: Global tensor mode id <-- Output tensor leg
      const std::vector<std::pair<unsigned int, unsigned int>> & bra_pairing, //in: bra pairing: Global tensor mode id <-- Output tensor leg

src/numerics/tensor_operator.hpp

 /** ExaTN::Numerics: Tensor operator
-REVISION: 2021/10/22
+REVISION: 2021/10/26
 
 Copyright (C) 2018-2021 Dmitry I. Lyakh (Liakh)
 Copyright (C) 2018-2021 Oak Ridge National Laboratory (UT-Battelle) **/
@@ -60,7 +60,26 @@ public:
  using Iterator = typename std::vector<OperatorComponent>::iterator;
  using ConstIterator = typename std::vector<OperatorComponent>::const_iterator;
 
- TensorOperator(const std::string & name): name_(name) {}
+ /** Creates an empty named tensor network operator. **/
+ TensorOperator(const std::string & name): name_(name) {} //in: tensor operator name
+
+ /** Creates a named tensor network operator from a single tensor network whose
+     legs are distributed among the ket and bra dimensions of the operator space map. **/
+ TensorOperator(const std::string & name,                                               //in: tensor operator name
+                std::shared_ptr<TensorNetwork> network,                                 //in: tensor network (or single tensor as a tensor network)
+                const std::vector<std::pair<unsigned int, unsigned int>> & ket_pairing, //in: ket pairing: Global tensor mode id <-- Output tensor leg
+                const std::vector<std::pair<unsigned int, unsigned int>> & bra_pairing, //in: bra pairing: Global tensor mode id <-- Output tensor leg
+                const std::complex<double> coefficient);                                //in: expansion coefficient
+
+ /** Creates a named tensor network operator from an outer product of two
+     tensor networks: |BraNetwork><KetNetwork|, where the KetNetwork connects
+     to the ket space and BraNetwork connects to the bra space. **/
+ TensorOperator(const std::string & name,                                               //in: tensor operator name
+                std::shared_ptr<TensorNetwork> ket_network,                             //in: ket tensor network (or single tensor as a tensor network)
+                std::shared_ptr<TensorNetwork> bra_network,                             //in: bra tensor network (or single tensor as a tensor network)
+                const std::vector<std::pair<unsigned int, unsigned int>> & ket_pairing, //in: ket pairing: Global tensor mode id <-- Output tensor leg (from the ket network)
+                const std::vector<std::pair<unsigned int, unsigned int>> & bra_pairing, //in: bra pairing: Global tensor mode id <-- Output tensor leg (from the bra network)
+                const std::complex<double> coefficient);                                //in: expansion coefficient
 
  TensorOperator(const TensorOperator &) = default;
  TensorOperator & operator=(const TensorOperator &) = default;