mps-0.0.0
MatrixProductStatesinC++
 All Data Structures Namespaces Functions Variables Enumerations Enumerator Groups Pages
itebd_itime.hpp
1 // -*- mode: c++; fill-column: 80; c-basic-offset: 2; indent-tabs-mode: nil -*-
2 /*
3  Copyright (c) 2012 Juan Jose Garcia Ripoll
4 
5  Tensor is free software; you can redistribute it and/or modify it
6  under the terms of the GNU Library General Public License as published
7  by the Free Software Foundation; either version 2 of the License, or
8  (at your option) any later version.
9 
10  This program is distributed in the hope that it will be useful,
11  but WITHOUT ANY WARRANTY; without even the implied warranty of
12  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13  GNU Library General Public License for more details.
14 
15  You should have received a copy of the GNU General Public License along
16  with this program; if not, write to the Free Software Foundation, Inc.,
17  51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
18 */
19 
20 #include <numeric>
21 #include <tensor/linalg.h>
22 #include <tensor/io.h>
23 #include <mps/itebd.h>
24 
25 namespace mps {
26 
27  template<class Tensor>
28  const iTEBD<Tensor>
29  evolve_itime(iTEBD<Tensor> psi, const Tensor &H12,
30  double dt, tensor::index nsteps,
31  double tolerance, tensor::index max_dim,
32  tensor::index deltan, int method,
33  std::vector<double> *energies,
34  std::vector<double> *entropies)
35  {
36  static const double FR_param[5] =
37  {0.67560359597983, 1.35120719195966, -0.17560359597983, -1.70241438391932};
38 
39  Tensor eH12[4];
40  //int method = 2;
41  switch (method) {
42  case 1:
43  /* Second order Trotter expansion */
44  eH12[1] = linalg::expm((-dt/2) * H12);
45  case 0:
46  /* First order Trotter expansion */
47  eH12[0] = linalg::expm((-dt) * H12);
48  break;
49  default:
50  /* Fourth order Trotter expansion */
51  for (int i = 0; i < 4; i++) {
52  eH12[i] = linalg::expm((-dt*FR_param[i]) * H12);
53  }
54  }
55  Tensor Id = Tensor::eye(H12.rows());
56  double time = 0;
57  psi = psi.canonical_form();
58  double E = energy(psi, H12), S = psi.entropy();
59  if (energies)
60  energies->push_back(E);
61  if (entropies)
62  entropies->push_back(S);
63 
64  if (!deltan) {
65  deltan = 1;
66  }
67  std::cout.precision(5);
68  std::cout << nsteps << ", " << dt << " x " << deltan
69  << " = " << dt * deltan << std::endl;
70  std::cout << "t=" << time
71  << ";\tE=" << E << "; dE=" << 0.0
72  << ";\tS=" << S << "; dS=" << 0.0
73  << ";\tl=" << std::max(psi.left_dimension(0),
74  psi.right_dimension(0))
75  << std::endl
76  << "l = " << matrix_form(real(psi.left_vector(0)))
77  << std::endl;
78  for (size_t phases = (nsteps + deltan - 1) / deltan; phases; phases--) {
79  for (size_t i = 0; (i < deltan); i++) {
80  switch (method) {
81  case 0:
82  psi = psi.apply_operator(eH12[0], 0, tolerance, max_dim);
83  psi = psi.apply_operator(eH12[0], 1, tolerance, max_dim);
84  break;
85  case 1:
86  psi = psi.apply_operator(eH12[1], 0, tolerance, max_dim);
87  psi = psi.apply_operator(eH12[0], 1, tolerance, max_dim);
88  psi = psi.apply_operator(eH12[1], 0, tolerance, max_dim);
89  break;
90  default:
91  psi = psi.apply_operator(eH12[0], 0, tolerance, max_dim);
92  psi = psi.apply_operator(eH12[1], 1, tolerance, max_dim);
93  psi = psi.apply_operator(eH12[2], 0, tolerance, max_dim);
94  psi = psi.apply_operator(eH12[3], 1, tolerance, max_dim);
95  psi = psi.apply_operator(eH12[2], 0, tolerance, max_dim);
96  psi = psi.apply_operator(eH12[1], 1, tolerance, max_dim);
97  psi = psi.apply_operator(eH12[0], 0, tolerance, max_dim);
98  }
99  time += dt;
100  }
101  psi = psi.canonical_form();
102  double newE = energy(psi, H12);
103  double newS = psi.entropy();
104  if (energies)
105  energies->push_back(newE);
106  if (entropies)
107  entropies->push_back(newS);
108  std::cout << "t=" << time
109  << ";\tE=" << newE << "; dE=" << newE-E
110  << ";\tS=" << newS << "; dS=" << newS-S
111  << ";\tl=" << std::max(psi.left_dimension(0),
112  psi.right_dimension(0))
113  << std::endl
114  << "l = " << matrix_form(real(psi.left_vector(0)))
115  << std::endl;
116  E = newE;
117  S = newS;
118  }
119  return psi;
120  }
121 
122 }
tensor::index left_dimension(int site) const
Return dimension of the MPS to the left of this site.
Definition: itebd.h:113
double energy(const RiTEBD &psi, const RTensor &Op12)
Expected value of the two-site operator Op12 acting on 'site' and 'site+1'.
tensor::index right_dimension(int site) const
Return dimension of the MPS to the right of this site.
Definition: itebd.h:118
const Tensor & left_vector(int site) const
Return the vector to the left of this site.
Definition: itebd.h:98
const iTEBD< Tensor > canonical_form(double tolerance=-1, tensor::index max_dim=0) const
Return a new state which is in canonical form.
const iTEBD< Tensor > apply_operator(const Tensor &U, int odd=0, double tolerance=-1, tensor::index max_dim=0) const
Construct a new state after acting on 'odd' or 'even' pair of sites with the two-site operator U...
An infinite Matrix Product State with translational invariance but using two tensors: one for odd and...
Definition: itebd.h:41
double entropy(int site) const
Estimate the entanglement entropy associated to splitting the state around 'site'.
Definition: itebd.hpp:83
const iTEBD< Tensor > evolve_itime(iTEBD< Tensor > psi, const Tensor &H12, double dt, tensor::index nsteps, double tolerance=-1, tensor::index max_dim=0, tensor::index deltan=1, int method=1, std::vector< double > *energies=0, std::vector< double > *entropies=0)
Evolve an iTEBD in imaginary time, using the local Hamiltonian H12 on state psi.
Definition: itebd_itime.hpp:29