Tensor Network Theory Library  Beta release 1.2.1 A library of routines for performing TNT-based operations
tntMps1sDmrg.c
1 /*
2  Authors: Sarah Al-Assam, Stephen Clark and Dieter Jaksch
3  Date: $LastChangedDate$
4  (c) University of Oxford 2013
5  */
6
11 /* See the paper PRB 91, 155115(2015) Hubig, McCulloch, Schollwoeck and Wolf (page 6) for a full exaplanation of these terms. */
12 double tntMpsVarMinMpo1sConvergenceFactor = 0.3;
13 double tntMpsVarMinMpo1sMixingFactor = 0.1;
14
15 /* Include the internal header for the TNT library */
16 #include "tntMpsInternal.h"
17 #include <math.h>
18
43 double tntMpsVarMinMpo1sSweep(tntNetwork mps,
44  tntLegLabel dir,
45  int chi,
46  tntNetwork mpo,
47  tntNodeArray *betas,
48  tntNodeArray *gammas,
49  double *err,
50  int start_site,
51  int end_site)
52 {
53  tntNode A; /* The current MPS node */
54  tntNode O; /* The current MPO node */
55  tntNode S, VT; /* The nodes returned from an SVD */
56  int L, j; /* number of sites and site index. */
57  double Ei, Emin, Ef; /* Initial energy, energy after optimisation, energy after truncation */
58  double DE_0, DE_T, ratio; /* Energy difference after optimisation, and after expansion and truncation */
59
60  /* Find the length of the MPS */
61  L = tntMpsLength(mps);
62
63  if ('R' == dir[0]) {
64  /* Set the current tensor equal to the start tensor */
65  A = tntNodeFindFirst(mps);
66  O = tntNodeFindFirst(mpo);
67  for (j = 0; j < start_site; j++) {
68  A = tntNodeFindConn(A,"R");
69  O = tntNodeFindConn(O,"R");
70  }
71
72  Ei = tntMpsVarMinMpo1sFullContract(betas->vals[j], gammas->vals[j], A, O);
73
74  for (j = start_site; j < end_site; j++) {
75  /* Perform local minimisation */
76  A = tntDmrg1sLocalOptimise(A, O, betas->vals[j], gammas->vals[j], &Emin);
77
78  /* Expand local subspace and truncate */
80  betas->vals[j], gammas->vals[j+1], "R", chi, err);
81
82  /* Update the node for beta with the new value of A */
83  tntMpsVarMinMpoUpdate(A, O, "R", j, betas);
84
85  /* Find other nodes that are the result of the SVD */
86  S = tntNodeFindConn(A,"R");
87  VT = tntNodeFindConn(S,"R");
88
89  /* Contract S and VT with the node to the right of them:
90  * this tensor will contain the orthogonality twist and be the new site to optimise */
91  A = tntNodeListContract(NULL, S, VT, tntNodeFindConn(VT,"R"));
92
93  /* Move operator along one site right */
94  O = tntNodeFindConn(O,"R");
95
96  /* Calculate the new energy */
97  Ef = tntMpsVarMinMpo1sFullContract(betas->vals[j+1], gammas->vals[j+1], A, O);
98
99  DE_0 = Ei-Emin;
100  DE_T = Emin-Ef;
101  ratio = -DE_T/(tntMpsVarMinMpo1sConvergenceFactor*DE_0);
102
103  // tntPrintf("j = %d, alpha = %2g, DE_0 = %7g, DE_T = %7g, -DE_T/DE_0 = %7g\n", j, tntMpsVarMinMpo1sMixingFactor, DE_0, DE_T, ratio);
104
105  if (DE_0 > 0) {
106  if (ratio > 1) {
107  tntMpsVarMinMpo1sMixingFactor = MAX(tntMpsVarMinMpo1sMixingFactor/pow(2,ratio-1),1e-12);
108  } else {
109  tntMpsVarMinMpo1sMixingFactor = MIN(tntMpsVarMinMpo1sMixingFactor*pow(2,1-ratio),1e3);
110  }
111  }
112
113  Ei = Ef;
114  }
115
116  } else {
117  /* Set the current tensor equal to the start tensor */
118  A = tntNodeFindLast(mps);
119  O = tntNodeFindLast(mpo);
120  for (j = L-1; j > start_site; j--) {
121  A = tntNodeFindConn(A,"L");
122  O = tntNodeFindConn(O,"L");
123  }
124
125  Ei = tntMpsVarMinMpo1sFullContract(betas->vals[j], gammas->vals[j], A, O);
126
127  for (j = start_site; j > end_site; j--) {
128
129  /* Perform local minimisation */
130  A = tntDmrg1sLocalOptimise(A, O, betas->vals[j], gammas->vals[j], &Emin);
131
132  A = tntDmrg1sSubspaceExpand(tntNodeFindConn(A,"L"), A, tntNodeFindConn(O,"L"), O,
133  betas->vals[j-1], gammas->vals[j], "L", chi, err);
134
135  /* Update the node for gamma with the new value of A */
136  tntMpsVarMinMpoUpdate(A, O, "L", j, gammas);
137
138  /* Find other nodes that are the result of the SVD */
139  S = tntNodeFindConn(A,"L");
140  VT = tntNodeFindConn(S,"L");
141
142  /* Contract S and VT with the node to the left of them:
143  * this tensor will contain the orthogonality twist and be the new site to optimise */
144  A = tntNodeListContract(NULL, S, VT, tntNodeFindConn(VT,"L"));
145
146  /* Move operator along one site left */
147  O = tntNodeFindConn(O,"L");
148
149  /* Calculate the new energy */
150  Ef = tntMpsVarMinMpo1sFullContract(betas->vals[j-1], gammas->vals[j-1], A, O);
151
152  DE_0 = Ei-Emin;
153  DE_T = Emin-Ef;
154  ratio = -DE_T/(tntMpsVarMinMpo1sConvergenceFactor*DE_0);
155
156  // tntPrintf("j = %d, alpha = %2g, DE_0 = %7g, DE_T = %7g, -DE_T/DE_0 = %7g\n", j, tntMpsVarMinMpo1sMixingFactor, DE_0, DE_T, ratio);
157
158  if (DE_0 > 0) {
159  if (ratio > 1) {
160  tntMpsVarMinMpo1sMixingFactor = MAX(tntMpsVarMinMpo1sMixingFactor/pow(2,ratio-1),1e-12);
161  } else {
162  tntMpsVarMinMpo1sMixingFactor = MIN(tntMpsVarMinMpo1sMixingFactor*pow(2,1-ratio),1e3);
163  }
164  }
165
166  Ei = Ef;
167  }
168  }
169
170  return Ef;
171 }
172
173
181 tntNode tntDmrg1sLocalOptimise(tntNode A,
182  tntNode O,
183  tntNode beta,
184  tntNode gamma,
185  double *eigvp)
186 {
187
188  tntNetwork Oeff; /* The network representing the effective operator for the current site */
189  tntNode A_eig; /* Optimisied MPS node */
190  tntComplex E_eig; /* The eigenvalue */
191
192  /* Build the network representing the effective operator for the two sites to be updated */
193  Oeff = tntMpsVarMinMpo1sBuild(beta,gamma,A,O);
194
195  /* The eigenvector of the Operator Oeff with the smallest eigenvalue is found */
196  A_eig = tntNetworkMinSite(A, &Oeff, 1, &tntMpsVarMinMpo1sContract, NULL, &E_eig);
197  *eigvp = E_eig.re;
198
199  /* Replace A_eig in the network */
200  tntNodeReplace(A, A_eig);
201  tntNodeFree(&A);
202
203  /* Free the networks and nodes no longer needed */
204  tntNetworkFree(&Oeff);
205
206  return A_eig;
207 }
208
216 tntNode tntDmrg1sSubspaceExpand(tntNode Al,
217  tntNode Ar,
218  tntNode Ol,
219  tntNode Or,
220  tntNode beta,
221  tntNode gamma,
222  tntLegLabel dir,
223  int chi,
224  double *errp)
225 {
226
227  tntNode P, Q; /* Nodes that will be used to expand the subspace to the left and right respectively */
228  tntNode Alc, Arc; /* Copies of nodes */
229  tntNode Ale, Are; /* Expanded versions of Al and Ar */
230  double errc; /* Error from the SVD */
231
232  /* Make copies of the original nodes */
233  Alc = tntNodeCopy(Al);
234  Ol = tntNodeCopy(Ol);
235  Arc = tntNodeCopy(Ar);
236  Or = tntNodeCopy(Or);
237  beta = tntNodeCopy(beta);
238  gamma = tntNodeCopy(gamma);
239
240  /* Get rid of the singleton legs on beta and gamma */
241  tntNodeSqueeze(beta,"LMN");
242  tntNodeSqueeze(gamma,"RST");
243
244  /* Join the nodes needed to make the left expander, then contract */
245  tntNodeJoin(beta,"R", Alc, "L");
246  tntNodeJoin(beta,"S", Ol, "L");
247  tntNodeJoin(Alc,"D", Ol, "U");
248
249  P = tntNodeListContract("LRDS", beta, Alc, Ol);
250
251  /* Now fuse the two right facing legs together */
252  tntNodeFuse(P,"RS","R", NULL);
253
254  /* Join the nodes needed to make the right expander, then contract */
255  tntNodeJoin(gamma,"L", Arc, "R");
256  tntNodeJoin(gamma,"M", Or, "R");
257  tntNodeJoin(Arc,"D", Or, "U");
258  Q = tntNodeListContract("RLDM", gamma, Arc, Or);
259
260  /* Now fuse the two right facing legs together */
261  tntNodeFuse(Q,"LM","L", NULL);
262
263  /* Scale the expanders appropriately depending on the direction */
264  if ('R' == dir[0]) {
265  tntNodeScaleReal(P, tntMpsVarMinMpo1sMixingFactor);
266  tntNodeScaleReal(Q, 0);
267  } else {
268  tntNodeScaleReal(Q, tntMpsVarMinMpo1sMixingFactor);
269  tntNodeScaleReal(P, 0);
270  }
271
272  /* Create the expanded nodes */
273  Ale = tntNodeDirectSum(P, Al, "R");
274  Are = tntNodeDirectSum(Q, Ar, "L");
275
276  tntNodeFree(&P);
277  tntNodeFree(&Q);
278
279  /* Put these nodes in the network in the place of the original left and right nodes.
280  * Before doing this the connections between them must first be removed */
281  tntNodeSplit(Al, Ar);
282  tntNodeReplace(Al, Ale);
283  tntNodeReplace(Ar, Are);
284  tntNodeJoin(Ale, "R", Are, "L");
285
286  tntNodeFree(&Al);
287  tntNodeFree(&Ar);
288
289  /* Take an SVD to truncate and also return A to orthonormal form */
290  if ('R' == dir[0]) {
291  Al = tntNodeSVD(Ale, "LD", "R", "L", "R", "L", chi, &errc, NULL);
292  *errp += errc;
293  return Al;
294  } else {
295  Ar = tntNodeSVD(Are, "RD", "L", "R", "L", "R", chi, &errc, NULL);
296  *errp += errc;
297  return Ar;
298  }
299 }
300
317 #ifdef MAKINGPUBLICDOCS
318 double tntMpsVarMinMpo1sStep(tntNetwork mps,
319 #else
320 double tntMpsVarMinMpo1sStep_(tntNetwork mps,
321 #endif
322  int chi,
323  tntNetwork mpo,
324  tntNodeArray *HeffL,
325  tntNodeArray *HeffR,
326  double *err,
327  unsigned start_site)
328 {
329
330  double E, errstep;
331  unsigned L = tntMpsLength(mps);
332
333  if (L-1 != start_site) E = tntMpsVarMinMpo1sSweep(mps, "R", chi, mpo, HeffL, HeffR, err, start_site, L-1);
334  E = tntMpsVarMinMpo1sSweep(mps, "L", chi, mpo, HeffL, HeffR, &errstep, L-1, 0);
335  if (0 != start_site) E = tntMpsVarMinMpo1sSweep(mps, "R", chi, mpo, HeffL, HeffR, err, 0, start_site);
336  *err += errstep;
337
338  return E;
339
340 }
void tntNodeFuse(tntNode A, tntLegLabel fuseLegs, tntLegLabel newLabelA, tntLegLabel newLabelB)
Definition: tntNodeConn.c:90
tntNode tntDmrg1sSubspaceExpand(tntNode Al, tntNode Ar, tntNode Ol, tntNode Or, tntNode beta, tntNode gamma, tntLegLabel dir, int chi, double *errp)
Definition: tntMps1sDmrg.c:216
tntNode tntNodeCopy(tntNode A)
Definition: tntNodeUtil.c:304
tntNode tntNodeDirectSum(tntNode A, tntNode B, tntLegLabel expandedLegs)
void tntNetworkFree(tntNetwork *nwp)
Definition: tntNetwork.c:118
void tntNodeJoin(tntNode A, tntLegLabel legA, tntNode B, tntLegLabel legB)
Definition: tntNodeConn.c:52
tntNode tntNodeFindLast(tntNetwork nw)
Definition: tntNodeInfo.c:68
double tntMpsVarMinMpo1sSweep(tntNetwork mps, tntLegLabel dir, int chi, tntNetwork mpo, tntNodeArray *betas, tntNodeArray *gammas, double *err, int start_site, int end_site)
Definition: tntMps1sDmrg.c:43
double tntMpsVarMinMpo1sStep(tntNetwork mps, int chi, tntNetwork mpo, tntNodeArray *HeffL, tntNodeArray *HeffR, double *err, unsigned start_site)
Definition: tntMps1sDmrg.c:318
void tntNodeSqueeze(tntNode A, tntLegLabel squeezeLegs)
Definition: tntNodeConn.c:292
tntNetwork tntMpsVarMinMpo1sBuild(tntNode beta, tntNode gamma, tntNode A, tntNode tno)
void tntNodeReplace(tntNode A, tntNode B)
Definition: tntNodeConn.c:145
tntNode tntNodeSVD(tntNode A, tntLegLabel Ulegs, tntLegLabel legUS, tntLegLabel legSU, tntLegLabel legSV, tntLegLabel legVS, int chi, double *err, tntLegLabel legmap)
tntNode tntDmrg1sLocalOptimise(tntNode A, tntNode O, tntNode beta, tntNode gamma, double *eigvp)
Definition: tntMps1sDmrg.c:181
void tntNodeFree(tntNode *A)
Definition: tntNodeUtil.c:275
unsigned tntMpsLength(tntNetwork wf)
Definition: tntMpsUtil.c:17
tntNode tntMpsVarMinMpo1sContract(tntNode A, tntNetwork nwMV)
Definition: tnt.h:65
tntNode tntNetworkMinSite(tntNode V, tntNetwork *nwMV, unsigned NM, tntNode(*eig_contract)(tntNode, tntNetwork), void(*eig_prep)(tntNode, tntNetwork), tntComplex *eigval)
Definition: tntNetwork.c:153
tntNode tntNodeFindConn(tntNode A, tntLegLabel legA)
Definition: tntNodeInfo.c:29
double re
Definition: tnt.h:66
tntNode tntNodeFindFirst(tntNetwork nw)
Definition: tntNodeInfo.c:50
void tntNodeScaleReal(tntNode A, double val)
void tntMpsVarMinMpoUpdate(tntNode tnA, tntNode tno, tntLegLabel dir, unsigned j, tntNodeArray *Heff)
double tntMpsVarMinMpo1sFullContract(tntNode beta, tntNode gamma, tntNode A, tntNode O)
void tntNodeSplit(tntNode A, tntNode B)
Definition: tntNodeConn.c:161