CheMPS2
TensorQ.cpp
1 /*
2  CheMPS2: a spin-adapted implementation of DMRG for ab initio quantum chemistry
3  Copyright (C) 2013-2016 Sebastian Wouters
4 
5  This program is free software; you can redistribute it and/or modify
6  it under the terms of the GNU General Public License as published by
7  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 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 <stdlib.h>
21 #include <math.h>
22 
23 #include "TensorQ.h"
24 #include "Lapack.h"
25 #include "Wigner.h"
26 
27 CheMPS2::TensorQ::TensorQ(const int boundary_index, const int Idiff, const bool moving_right, const SyBookkeeper * denBK, const Problem * Prob, const int site) :
28 TensorOperator(boundary_index,
29  1, //two_j
30  1, //n_elec
31  Idiff,
32  moving_right,
33  true, //prime_last
34  true, //jw_phase (three 2nd quantized operators)
35  denBK,
36  denBK){
37 
38  this->Prob = Prob;
39  this->site = site;
40 
41 }
42 
43 CheMPS2::TensorQ::~TensorQ(){ }
44 
45 void CheMPS2::TensorQ::AddTermSimple(TensorT * denT){
46 
47  if (( moving_right) && (bk_up->gIrrep(denT->gIndex()) == n_irrep)){ AddTermSimpleRight(denT); }
48  if ((!moving_right) && (bk_up->gIrrep(denT->gIndex()) == n_irrep)){ AddTermSimpleLeft( denT); }
49 
50 }
51 
52 void CheMPS2::TensorQ::AddTermSimpleRight(TensorT * denT){
53 
54  const double mxElement = Prob->gMxElement(index-1,index-1,index-1,site);
55 
56  for (int ikappa=0; ikappa<nKappa; ikappa++){
57  const int ID = Irreps::directProd( n_irrep, sector_irrep_up[ikappa] );
58  int dimRU = bk_up->gCurrentDim(index, sector_nelec_up[ikappa], sector_spin_up[ikappa], sector_irrep_up[ikappa]);
59  int dimRD = bk_up->gCurrentDim(index, sector_nelec_up[ikappa]+1, sector_spin_down[ikappa], ID);
60  int dimL = bk_up->gCurrentDim(index-1, sector_nelec_up[ikappa]-1, sector_spin_down[ikappa], ID);
61  if (dimL>0){
62 
63  double * BlockTup = denT->gStorage(sector_nelec_up[ikappa]-1, sector_spin_down[ikappa], ID, sector_nelec_up[ikappa], sector_spin_up[ikappa], sector_irrep_up[ikappa]);
64  double * BlockTdo = denT->gStorage(sector_nelec_up[ikappa]-1, sector_spin_down[ikappa], ID, sector_nelec_up[ikappa]+1, sector_spin_down[ikappa], ID);
65 
66  int fase = ((((sector_spin_down[ikappa]+1-sector_spin_up[ikappa])/2)%2)!=0)?-1:1;
67  double alpha = fase * mxElement * sqrt((sector_spin_up[ikappa]+1.0)/(sector_spin_down[ikappa]+1.0));
68  double beta = 1.0; //add
69  char trans = 'T';
70  char notr = 'N';
71 
72  dgemm_(&trans,&notr,&dimRU,&dimRD,&dimL,&alpha,BlockTup,&dimL,BlockTdo,&dimL,&beta,storage+kappa2index[ikappa],&dimRU);
73 
74  }
75  }
76 
77 }
78 
79 void CheMPS2::TensorQ::AddTermSimpleLeft(TensorT * denT){
80 
81  const double mxElement = Prob->gMxElement(site,index,index,index);
82 
83  for (int ikappa=0; ikappa<nKappa; ikappa++){
84  const int ID = Irreps::directProd( n_irrep, sector_irrep_up[ikappa] );
85  int dimLU = bk_up->gCurrentDim(index, sector_nelec_up[ikappa], sector_spin_up[ikappa], sector_irrep_up[ikappa]);
86  int dimLD = bk_up->gCurrentDim(index, sector_nelec_up[ikappa]+1, sector_spin_down[ikappa], ID);
87  int dimR = bk_up->gCurrentDim(index+1, sector_nelec_up[ikappa]+2, sector_spin_up[ikappa], sector_irrep_up[ikappa]);
88  if (dimR>0){
89 
90  double * BlockTup = denT->gStorage(sector_nelec_up[ikappa], sector_spin_up[ikappa], sector_irrep_up[ikappa], sector_nelec_up[ikappa]+2, sector_spin_up[ikappa], sector_irrep_up[ikappa]);
91  double * BlockTdo = denT->gStorage(sector_nelec_up[ikappa]+1, sector_spin_down[ikappa], ID, sector_nelec_up[ikappa]+2, sector_spin_up[ikappa], sector_irrep_up[ikappa]);
92 
93  int fase = ((((sector_spin_up[ikappa]+1-sector_spin_down[ikappa])/2)%2)!=0)?-1:1;
94  double alpha = fase * mxElement * sqrt((sector_spin_up[ikappa]+1.0)/(sector_spin_down[ikappa]+1.0));
95  double beta = 1.0; //add
96  char trans = 'T';
97  char notr = 'N';
98 
99  dgemm_(&notr,&trans,&dimLU,&dimLD,&dimR,&alpha,BlockTup,&dimLU,BlockTdo,&dimLD,&beta,storage+kappa2index[ikappa],&dimLU);
100 
101  }
102  }
103 
104 }
105 
106 void CheMPS2::TensorQ::AddTermsL(TensorL ** Ltensors, TensorT * denT, double * workmem, double * workmem2){
107 
108  if (moving_right){ AddTermsLRight(Ltensors, denT, workmem, workmem2); }
109  else{ AddTermsLLeft( Ltensors, denT, workmem, workmem2); }
110 
111 }
112 
113 void CheMPS2::TensorQ::AddTermsLRight(TensorL ** Ltensors, TensorT * denT, double * workmem, double * workmem2){
114 
115  bool OneToAdd = false;
116  for (int loca=0; loca<index-1; loca++){
117  if (Ltensors[index-2-loca]->get_irrep() == n_irrep){ OneToAdd = true; }
118  }
119 
120  if (OneToAdd){
121  for (int ikappa=0; ikappa<nKappa; ikappa++){
122 
123  const int ID = Irreps::directProd( sector_irrep_up[ikappa], n_irrep );
124  int dimRU = bk_up->gCurrentDim(index, sector_nelec_up[ikappa], sector_spin_up[ikappa], sector_irrep_up[ikappa]);
125  int dimRD = bk_up->gCurrentDim(index, sector_nelec_up[ikappa]+1, sector_spin_down[ikappa], ID);
126 
127  //case 1
128  int dimLU = bk_up->gCurrentDim(index-1, sector_nelec_up[ikappa], sector_spin_up[ikappa], sector_irrep_up[ikappa]);
129  int dimLD = bk_up->gCurrentDim(index-1, sector_nelec_up[ikappa]-1, sector_spin_down[ikappa], ID);
130 
131  if ((dimLU>0) && (dimLD>0)){
132 
133  int dimLUxLD = dimLU * dimLD;
134  for (int cnt=0; cnt<dimLUxLD; cnt++){ workmem[cnt] = 0.0; }
135 
136  for (int loca=0; loca<index-1; loca++){
137  if (Ltensors[index-2-loca]->get_irrep() == n_irrep){
138  double * BlockL = Ltensors[index-2-loca]->gStorage(sector_nelec_up[ikappa]-1,sector_spin_down[ikappa],ID,sector_nelec_up[ikappa],sector_spin_up[ikappa],sector_irrep_up[ikappa]);
139  double alpha = Prob->gMxElement(loca,site,index-1,index-1);
140  int inc = 1;
141  daxpy_(&dimLUxLD, &alpha, BlockL, &inc, workmem, &inc);
142  }
143  }
144 
145  int fase = ((((sector_spin_down[ikappa]+1-sector_spin_up[ikappa])/2)%2)!=0)?-1:1;
146  double alpha = fase * sqrt((sector_spin_up[ikappa]+1.0)/(sector_spin_down[ikappa]+1.0));
147  double beta = 0.0; //set
148 
149  double * BlockTup = denT->gStorage(sector_nelec_up[ikappa], sector_spin_up[ikappa], sector_irrep_up[ikappa], sector_nelec_up[ikappa], sector_spin_up[ikappa], sector_irrep_up[ikappa]);
150  double * BlockTdo = denT->gStorage(sector_nelec_up[ikappa]-1, sector_spin_down[ikappa], ID, sector_nelec_up[ikappa]+1, sector_spin_down[ikappa], ID);
151 
152  char totrans = 'T';
153  // factor * Tup^T * L^T --> mem2
154  dgemm_(&totrans,&totrans,&dimRU,&dimLD,&dimLU,&alpha,BlockTup,&dimLU,workmem,&dimLD,&beta,workmem2,&dimRU);
155 
156  alpha = 1.0;
157  beta = 1.0; //add
158  totrans = 'N';
159  // mem2 * Tdo --> storage
160  dgemm_(&totrans,&totrans,&dimRU,&dimRD,&dimLD,&alpha,workmem2,&dimRU, BlockTdo, &dimLD, &beta, storage+kappa2index[ikappa], &dimRU);
161 
162  }
163 
164  //case 2
165  dimLU = bk_up->gCurrentDim(index-1, sector_nelec_up[ikappa]-2, sector_spin_up[ikappa], sector_irrep_up[ikappa]);
166  //dimLD same as case1
167  if ((dimLU>0) && (dimLD>0)){
168 
169  int dimLUxLD = dimLU * dimLD;
170  for (int cnt=0; cnt<dimLUxLD; cnt++){ workmem[cnt] = 0.0; }
171 
172  for (int loca=0; loca<index-1; loca++){
173  if (Ltensors[index-2-loca]->get_irrep() == n_irrep){
174  double * BlockL = Ltensors[index-2-loca]->gStorage(sector_nelec_up[ikappa]-2,sector_spin_up[ikappa],sector_irrep_up[ikappa],sector_nelec_up[ikappa]-1,sector_spin_down[ikappa],ID);
175  double alpha = 2*Prob->gMxElement(loca,index-1,site,index-1) - Prob->gMxElement(loca,index-1,index-1,site);
176  int inc = 1;
177  daxpy_(&dimLUxLD, &alpha, BlockL, &inc, workmem, &inc);
178  }
179  }
180 
181  double alpha = 1.0; //factor = 1 in this case
182  double beta = 0.0; //set
183 
184  double * BlockTup = denT->gStorage(sector_nelec_up[ikappa]-2, sector_spin_up[ikappa], sector_irrep_up[ikappa], sector_nelec_up[ikappa], sector_spin_up[ikappa], sector_irrep_up[ikappa]);
185  double * BlockTdo = denT->gStorage(sector_nelec_up[ikappa]-1, sector_spin_down[ikappa], ID, sector_nelec_up[ikappa]+1, sector_spin_down[ikappa], ID);
186 
187  char trans = 'T';
188  char notr = 'N';
189  // factor * Tup^T * L --> mem2
190  dgemm_(&trans,&notr,&dimRU,&dimLD,&dimLU,&alpha,BlockTup,&dimLU,workmem,&dimLU,&beta,workmem2,&dimRU);
191 
192  beta = 1.0; //add
193  // mem2 * Tdo --> storage
194  dgemm_(&notr,&notr,&dimRU,&dimRD,&dimLD,&alpha,workmem2,&dimRU, BlockTdo, &dimLD, &beta, storage+kappa2index[ikappa], &dimRU);
195 
196  }
197 
198  //case 3
199  for (int TwoSLU=sector_spin_up[ikappa]-1; TwoSLU<=sector_spin_up[ikappa]+1; TwoSLU+=2){
200  for (int TwoSLD=sector_spin_down[ikappa]-1; TwoSLD<=sector_spin_down[ikappa]+1; TwoSLD+=2){
201  if ((TwoSLD>=0) && (TwoSLU>=0) && (abs(TwoSLD-TwoSLU)<2)){
202  const int ILU = Irreps::directProd(sector_irrep_up[ikappa],bk_up->gIrrep(index-1));
203  const int ILD = Irreps::directProd(ID, bk_up->gIrrep(index-1));
204  dimLU = bk_up->gCurrentDim(index-1, sector_nelec_up[ikappa]-1, TwoSLU, ILU);
205  dimLD = bk_up->gCurrentDim(index-1, sector_nelec_up[ikappa] , TwoSLD, ILD);
206  if ((dimLU>0) && (dimLD>0)){
207  int fase = ((((sector_spin_up[ikappa]+TwoSLD)/2)%2)!=0)?-1:1;
208  double factor = fase * sqrt((TwoSLD+1)*(sector_spin_up[ikappa]+1.0))
209  * Wigner::wigner6j( sector_spin_up[ikappa], sector_spin_down[ikappa], 1, TwoSLD, TwoSLU, 1 );
210 
211  int dimLUxLD = dimLU * dimLD;
212  for (int cnt=0; cnt<dimLUxLD; cnt++){ workmem[cnt] = 0.0; }
213 
214  for (int loca=0; loca<index-1; loca++){
215  if (Ltensors[index-2-loca]->get_irrep() == n_irrep){
216  double * BlockL = Ltensors[index-2-loca]->gStorage(sector_nelec_up[ikappa]-1, TwoSLU, ILU, sector_nelec_up[ikappa], TwoSLD, ILD);
217  double alpha = factor * Prob->gMxElement(loca,index-1,site,index-1);
218  if (TwoSLD==sector_spin_up[ikappa]){ alpha += Prob->gMxElement(loca,index-1,index-1,site); }
219  int inc = 1;
220  daxpy_(&dimLUxLD, &alpha, BlockL, &inc, workmem, &inc);
221  }
222  }
223 
224  double alpha = 1.0;
225  double beta = 0.0; //set
226 
227  double * BlockTup = denT->gStorage(sector_nelec_up[ikappa]-1, TwoSLU, ILU, sector_nelec_up[ikappa], sector_spin_up[ikappa], sector_irrep_up[ikappa]);
228  double * BlockTdo = denT->gStorage(sector_nelec_up[ikappa] , TwoSLD, ILD, sector_nelec_up[ikappa]+1, sector_spin_down[ikappa], ID);
229 
230  char trans = 'T';
231  char notr = 'N';
232  // Tup^T * mem --> mem2
233  dgemm_(&trans,&notr,&dimRU,&dimLD,&dimLU,&alpha,BlockTup,&dimLU,workmem,&dimLU,&beta,workmem2,&dimRU);
234 
235  beta = 1.0; //add
236  // mem2 * Tdo --> storage
237  dgemm_(&notr,&notr,&dimRU,&dimRD,&dimLD,&alpha,workmem2,&dimRU, BlockTdo, &dimLD, &beta, storage+kappa2index[ikappa], &dimRU);
238 
239  }
240  }
241  }
242  }
243 
244  }
245  }
246 
247 }
248 
249 void CheMPS2::TensorQ::AddTermsLLeft(TensorL ** Ltensors, TensorT * denT, double * workmem, double * workmem2){
250 
251  bool OneToAdd = false;
252  for (int loca=index+1; loca<Prob->gL(); loca++){
253  if (Ltensors[loca-index-1]->get_irrep() == n_irrep){ OneToAdd = true; }
254  }
255 
256  if (OneToAdd){
257  for (int ikappa=0; ikappa<nKappa; ikappa++){
258 
259  const int ID = Irreps::directProd( sector_irrep_up[ikappa], n_irrep );
260  int dimLU = bk_up->gCurrentDim(index, sector_nelec_up[ikappa], sector_spin_up[ikappa], sector_irrep_up[ikappa]);
261  int dimLD = bk_up->gCurrentDim(index, sector_nelec_up[ikappa]+1, sector_spin_down[ikappa], ID);
262 
263  //case 1
264  int dimRU = bk_up->gCurrentDim(index+1, sector_nelec_up[ikappa]+2, sector_spin_up[ikappa], sector_irrep_up[ikappa]);
265  int dimRD = bk_up->gCurrentDim(index+1, sector_nelec_up[ikappa]+1, sector_spin_down[ikappa], ID);
266 
267  if ((dimRU>0) && (dimRD>0)){
268 
269  int dimRUxRD = dimRU * dimRD;
270  for (int cnt=0; cnt<dimRUxRD; cnt++){ workmem[cnt] = 0.0; }
271 
272  for (int loca=index+1; loca<Prob->gL(); loca++){
273  if (Ltensors[loca-index-1]->get_irrep() == n_irrep){
274  double * BlockL = Ltensors[loca-index-1]->gStorage(sector_nelec_up[ikappa]+1,sector_spin_down[ikappa],ID,sector_nelec_up[ikappa]+2,sector_spin_up[ikappa],sector_irrep_up[ikappa]);
275  double alpha = Prob->gMxElement(site,loca,index,index);
276  int inc = 1;
277  daxpy_(&dimRUxRD, &alpha, BlockL, &inc, workmem, &inc);
278  }
279  }
280 
281  int fase = ((((sector_spin_up[ikappa]+1-sector_spin_down[ikappa])/2)%2)!=0)?-1:1;
282  double alpha = fase * sqrt((sector_spin_up[ikappa]+1.0)/(sector_spin_down[ikappa]+1.0));
283  double beta = 0.0; //set
284 
285  double * BlockTup = denT->gStorage(sector_nelec_up[ikappa], sector_spin_up[ikappa], sector_irrep_up[ikappa], sector_nelec_up[ikappa]+2, sector_spin_up[ikappa], sector_irrep_up[ikappa]);
286  double * BlockTdo = denT->gStorage(sector_nelec_up[ikappa]+1, sector_spin_down[ikappa], ID, sector_nelec_up[ikappa]+1, sector_spin_down[ikappa], ID);
287 
288  char trans = 'T';
289  char notr = 'N';
290  // factor * Tup * L^T --> mem2
291  dgemm_(&notr,&trans,&dimLU,&dimRD,&dimRU,&alpha,BlockTup,&dimLU,workmem,&dimRD,&beta,workmem2,&dimLU);
292 
293  alpha = 1.0;
294  beta = 1.0; //add
295  // mem2 * Tdo^T --> storage
296  dgemm_(&notr,&trans,&dimLU,&dimLD,&dimRD,&alpha,workmem2,&dimLU, BlockTdo, &dimLD, &beta, storage+kappa2index[ikappa], &dimLU);
297 
298  }
299 
300  //case 2
301  dimRD = bk_up->gCurrentDim(index+1, sector_nelec_up[ikappa]+3, sector_spin_down[ikappa], ID);
302  //dimRU same as case1
303  if ((dimRU>0) && (dimRD>0)){
304 
305  int dimRUxRD = dimRU * dimRD;
306  for (int cnt=0; cnt<dimRUxRD; cnt++){ workmem[cnt] = 0.0; }
307 
308  for (int loca=index+1; loca<Prob->gL(); loca++){
309  if (Ltensors[loca-index-1]->get_irrep() == n_irrep){
310  double * BlockL = Ltensors[loca-index-1]->gStorage(sector_nelec_up[ikappa]+2,sector_spin_up[ikappa],sector_irrep_up[ikappa],sector_nelec_up[ikappa]+3,sector_spin_down[ikappa],ID);
311  double alpha = 2*Prob->gMxElement(site,index,loca,index) - Prob->gMxElement(site,index,index,loca);
312  int inc = 1;
313  daxpy_(&dimRUxRD, &alpha, BlockL, &inc, workmem, &inc);
314  }
315  }
316 
317  double alpha = 1.0; //factor = 1 in this case
318  double beta = 0.0; //set
319 
320  double * BlockTup = denT->gStorage(sector_nelec_up[ikappa], sector_spin_up[ikappa], sector_irrep_up[ikappa], sector_nelec_up[ikappa]+2, sector_spin_up[ikappa], sector_irrep_up[ikappa]);
321  double * BlockTdo = denT->gStorage(sector_nelec_up[ikappa]+1, sector_spin_down[ikappa], ID, sector_nelec_up[ikappa]+3, sector_spin_down[ikappa], ID);
322 
323  char notr = 'N';
324  // factor * Tup * L --> mem2
325  dgemm_(&notr,&notr,&dimLU,&dimRD,&dimRU,&alpha,BlockTup,&dimLU,workmem,&dimRU,&beta,workmem2,&dimLU);
326 
327  beta = 1.0; //add
328  // mem2 * Tdo^T --> storage
329  char trans = 'T';
330  dgemm_(&notr,&trans,&dimLU,&dimLD,&dimRD,&alpha,workmem2,&dimLU, BlockTdo, &dimLD, &beta, storage+kappa2index[ikappa], &dimLU);
331 
332  }
333 
334  //case 3
335  for (int TwoSRU=sector_spin_up[ikappa]-1; TwoSRU<=sector_spin_up[ikappa]+1; TwoSRU+=2){
336  for (int TwoSRD=sector_spin_down[ikappa]-1; TwoSRD<=sector_spin_down[ikappa]+1; TwoSRD+=2){
337  if ((TwoSRD>=0) && (TwoSRU>=0) && (abs(TwoSRD-TwoSRU)<2)){
338  const int IRU = Irreps::directProd(sector_irrep_up[ikappa],bk_up->gIrrep(index));
339  const int IRD = Irreps::directProd(ID, bk_up->gIrrep(index));
340  dimRU = bk_up->gCurrentDim(index+1, sector_nelec_up[ikappa]+1, TwoSRU, IRU);
341  dimRD = bk_up->gCurrentDim(index+1, sector_nelec_up[ikappa]+2, TwoSRD, IRD);
342  if ((dimRU>0) && (dimRD>0)){
343  int fase = ((((sector_spin_down[ikappa]+TwoSRU)/2)%2)!=0)?-1:1;
344  double factor1 = fase * sqrt((TwoSRU+1.0)/(sector_spin_down[ikappa]+1.0)) * (TwoSRD+1)
345  * Wigner::wigner6j( sector_spin_up[ikappa], sector_spin_down[ikappa], 1, TwoSRD, TwoSRU, 1 );
346  double factor2 = (TwoSRD+1.0)/(sector_spin_down[ikappa]+1.0);
347 
348  int dimRUxRD = dimRU * dimRD;
349  for (int cnt=0; cnt<dimRUxRD; cnt++){ workmem[cnt] = 0.0; }
350 
351  for (int loca=index+1; loca<Prob->gL(); loca++){
352  if (Ltensors[loca-index-1]->get_irrep() == n_irrep){
353  double * BlockL = Ltensors[loca-index-1]->gStorage(sector_nelec_up[ikappa]+1, TwoSRU, IRU, sector_nelec_up[ikappa]+2, TwoSRD, IRD);
354  double alpha = factor1 * Prob->gMxElement(site,index,loca,index);
355  if (TwoSRU==sector_spin_down[ikappa]){ alpha += factor2 * Prob->gMxElement(site,index,index,loca); }
356  int inc = 1;
357  daxpy_(&dimRUxRD, &alpha, BlockL, &inc, workmem, &inc);
358  }
359  }
360 
361  double alpha = 1.0;
362  double beta = 0.0; //set
363 
364  double * BlockTup = denT->gStorage(sector_nelec_up[ikappa], sector_spin_up[ikappa], sector_irrep_up[ikappa], sector_nelec_up[ikappa]+1, TwoSRU, IRU);
365  double * BlockTdo = denT->gStorage(sector_nelec_up[ikappa]+1, sector_spin_down[ikappa], ID, sector_nelec_up[ikappa]+2, TwoSRD, IRD);
366 
367  char notr = 'N';
368  // Tup * mem --> mem2
369  dgemm_(&notr,&notr,&dimLU,&dimRD,&dimRU,&alpha,BlockTup,&dimLU,workmem,&dimRU,&beta,workmem2,&dimLU);
370 
371  beta = 1.0; //add
372  // mem2 * Tdo^T --> storage
373  char trans = 'T';
374  dgemm_(&notr,&trans,&dimLU,&dimLD,&dimRD,&alpha,workmem2,&dimLU, BlockTdo, &dimLD, &beta, storage+kappa2index[ikappa], &dimLU);
375 
376  }
377  }
378  }
379  }
380 
381  }
382  }
383 
384 }
385 
386 void CheMPS2::TensorQ::AddTermsAB(TensorOperator * denA, TensorOperator * denB, TensorT * denT, double * workmem, double * workmem2){
387 
388  if (moving_right){ AddTermsABRight(denA, denB, denT, workmem, workmem2); }
389  else{ AddTermsABLeft( denA, denB, denT, workmem, workmem2); }
390 
391 }
392 
393 void CheMPS2::TensorQ::AddTermsABRight(TensorOperator * denA, TensorOperator * denB, TensorT * denT, double * workmem, double * workmem2){
394 
395  for (int ikappa=0; ikappa<nKappa; ikappa++){
396  const int ID = Irreps::directProd( sector_irrep_up[ikappa], n_irrep );
397  int dimRU = bk_up->gCurrentDim(index, sector_nelec_up[ikappa], sector_spin_up[ikappa], sector_irrep_up[ikappa]);
398  int dimRD = bk_up->gCurrentDim(index, sector_nelec_up[ikappa]+1, sector_spin_down[ikappa], ID);
399 
400  //case 1
401  const int ILU = Irreps::directProd(sector_irrep_up[ikappa],bk_up->gIrrep(index-1));
402  for (int TwoSLU=sector_spin_up[ikappa]-1; TwoSLU<=sector_spin_up[ikappa]+1; TwoSLU+=2){
403  int dimLU = bk_up->gCurrentDim(index-1, sector_nelec_up[ikappa]-1, TwoSLU, ILU);
404  int dimLD = bk_up->gCurrentDim(index-1, sector_nelec_up[ikappa]+1, sector_spin_down[ikappa], ID);
405  if ((dimLU>0) && (dimLD>0)){
406 
407  int fase = ((((TwoSLU + sector_spin_down[ikappa] + 2)/2)%2)!=0)?-1:1;
408  double factorB = fase * sqrt(3.0*(sector_spin_up[ikappa]+1))
409  * Wigner::wigner6j( 1, 2, 1, sector_spin_down[ikappa], sector_spin_up[ikappa], TwoSLU );
410 
411  double alpha;
412  double * mem;
413 
414  if (TwoSLU == sector_spin_down[ikappa]){
415 
416  fase = ((((sector_spin_down[ikappa]+1-sector_spin_up[ikappa])/2)%2)!=0)?-1:1;
417  double factorA = fase * sqrt( 0.5 * (sector_spin_up[ikappa]+1.0) / (sector_spin_down[ikappa]+1.0) );
418 
419  double * BlockA = denA->gStorage( sector_nelec_up[ikappa]-1,TwoSLU,ILU,sector_nelec_up[ikappa]+1,sector_spin_down[ikappa],ID );
420  double * BlockB = denB->gStorage( sector_nelec_up[ikappa]-1,TwoSLU,ILU,sector_nelec_up[ikappa]+1,sector_spin_down[ikappa],ID );
421 
422  mem = workmem;
423  for (int cnt=0; cnt<dimLU*dimLD; cnt++){ mem[cnt] = factorA * BlockA[cnt] + factorB * BlockB[cnt]; }
424  alpha = 1.0;
425 
426  } else {
427 
428  alpha = factorB;
429  mem = denB->gStorage(sector_nelec_up[ikappa]-1,TwoSLU,ILU,sector_nelec_up[ikappa]+1,sector_spin_down[ikappa],ID);
430 
431  }
432 
433  double * BlockTup = denT->gStorage(sector_nelec_up[ikappa]-1, TwoSLU, ILU, sector_nelec_up[ikappa], sector_spin_up[ikappa], sector_irrep_up[ikappa]);
434  double * BlockTdo = denT->gStorage(sector_nelec_up[ikappa]+1, sector_spin_down[ikappa], ID, sector_nelec_up[ikappa]+1, sector_spin_down[ikappa], ID);
435 
436  char trans = 'T';
437  char notr = 'N';
438  double beta = 0.0; //set
439  dgemm_(&trans,&notr,&dimRU,&dimLD,&dimLU,&alpha,BlockTup,&dimLU,mem,&dimLU,&beta,workmem2,&dimRU);
440 
441  alpha = 1.0;
442  beta = 1.0; //add
443  dgemm_(&notr,&notr,&dimRU,&dimRD,&dimLD,&alpha,workmem2,&dimRU,BlockTdo,&dimLD,&beta,storage+kappa2index[ikappa],&dimRU);
444 
445  }
446  }
447 
448  //case 2
449  const int ILD = Irreps::directProd(ID,bk_up->gIrrep(index-1));
450  for (int TwoSLD=sector_spin_down[ikappa]-1; TwoSLD<=sector_spin_down[ikappa]+1; TwoSLD+=2){
451  int dimLU = bk_up->gCurrentDim(index-1, sector_nelec_up[ikappa]-2, sector_spin_up[ikappa], sector_irrep_up[ikappa]);
452  int dimLD = bk_up->gCurrentDim(index-1, sector_nelec_up[ikappa], TwoSLD, ILD);
453  if ((dimLU>0) && (dimLD>0)){
454 
455  int fase = ((((sector_spin_up[ikappa] + sector_spin_down[ikappa] + 1)/2)%2)!=0)?-1:1;
456  double factorB = fase * sqrt(3.0*(TwoSLD+1)) * Wigner::wigner6j( 1, 2, 1, sector_spin_up[ikappa], sector_spin_down[ikappa], TwoSLD );
457 
458  double alpha;
459  double * mem;
460 
461  if (TwoSLD == sector_spin_up[ikappa]){
462 
463  double factorA = - sqrt(0.5);
464 
465  double * BlockA = denA->gStorage( sector_nelec_up[ikappa]-2, sector_spin_up[ikappa], sector_irrep_up[ikappa], sector_nelec_up[ikappa], TwoSLD, ILD );
466  double * BlockB = denB->gStorage( sector_nelec_up[ikappa]-2, sector_spin_up[ikappa], sector_irrep_up[ikappa], sector_nelec_up[ikappa], TwoSLD, ILD );
467 
468  mem = workmem;
469  for (int cnt=0; cnt<dimLU*dimLD; cnt++){ mem[cnt] = factorA * BlockA[cnt] + factorB * BlockB[cnt]; }
470  alpha = 1.0;
471 
472  } else {
473 
474  alpha = factorB;
475  mem = denB->gStorage( sector_nelec_up[ikappa]-2, sector_spin_up[ikappa], sector_irrep_up[ikappa], sector_nelec_up[ikappa], TwoSLD, ILD );
476 
477  }
478 
479  double * BlockTup = denT->gStorage(sector_nelec_up[ikappa]-2,sector_spin_up[ikappa],sector_irrep_up[ikappa],sector_nelec_up[ikappa], sector_spin_up[ikappa], sector_irrep_up[ikappa]);
480  double * BlockTdo = denT->gStorage(sector_nelec_up[ikappa] ,TwoSLD, ILD, sector_nelec_up[ikappa]+1, sector_spin_down[ikappa], ID);
481 
482  char trans = 'T';
483  char notr = 'N';
484  double beta = 0.0; //set
485  dgemm_(&trans,&notr,&dimRU,&dimLD,&dimLU,&alpha,BlockTup,&dimLU,mem,&dimLU,&beta,workmem2,&dimRU);
486 
487  alpha = 1.0;
488  beta = 1.0; //add
489  dgemm_(&notr,&notr,&dimRU,&dimRD,&dimLD,&alpha,workmem2,&dimRU,BlockTdo,&dimLD,&beta,storage+kappa2index[ikappa],&dimRU);
490 
491  }
492  }
493  }
494 
495 }
496 
497 void CheMPS2::TensorQ::AddTermsABLeft(TensorOperator * denA, TensorOperator * denB, TensorT * denT, double * workmem, double * workmem2){
498 
499  for (int ikappa=0; ikappa<nKappa; ikappa++){
500  const int ID = Irreps::directProd( sector_irrep_up[ikappa], n_irrep );
501  int dimLU = bk_up->gCurrentDim(index, sector_nelec_up[ikappa], sector_spin_up[ikappa], sector_irrep_up[ikappa]);
502  int dimLD = bk_up->gCurrentDim(index, sector_nelec_up[ikappa]+1, sector_spin_down[ikappa], ID);
503 
504  //case 1
505  const int IRD = Irreps::directProd(ID, bk_up->gIrrep(index));
506  for (int TwoSRD=sector_spin_down[ikappa]-1; TwoSRD<=sector_spin_down[ikappa]+1; TwoSRD+=2){
507  int dimRU = bk_up->gCurrentDim(index+1, sector_nelec_up[ikappa], sector_spin_up[ikappa], sector_irrep_up[ikappa]);
508  int dimRD = bk_up->gCurrentDim(index+1, sector_nelec_up[ikappa]+2, TwoSRD, IRD);
509  if ((dimRU>0) && (dimRD>0)){
510 
511  int fase = ((((TwoSRD + sector_spin_up[ikappa] + 2)/2)%2)!=0)?-1:1;
512  const double factorB = fase * sqrt(3.0/(sector_spin_down[ikappa]+1.0)) * (TwoSRD+1)
513  * Wigner::wigner6j( 1, 1, 2, sector_spin_up[ikappa], TwoSRD, sector_spin_down[ikappa] );
514 
515  double alpha;
516  double * mem;
517 
518  if (TwoSRD == sector_spin_up[ikappa]){
519 
520  fase = ((((sector_spin_up[ikappa]+1-sector_spin_down[ikappa])/2)%2)!=0)?-1:1;
521  double factorA = fase * sqrt( 0.5 * (sector_spin_up[ikappa]+1.0) / (sector_spin_down[ikappa]+1.0) );
522 
523  double * BlockA = denA->gStorage( sector_nelec_up[ikappa], sector_spin_up[ikappa], sector_irrep_up[ikappa], sector_nelec_up[ikappa]+2, TwoSRD, IRD );
524  double * BlockB = denB->gStorage( sector_nelec_up[ikappa], sector_spin_up[ikappa], sector_irrep_up[ikappa], sector_nelec_up[ikappa]+2, TwoSRD, IRD );
525 
526  mem = workmem;
527  for (int cnt=0; cnt<dimRU*dimRD; cnt++){ mem[cnt] = factorA * BlockA[cnt] + factorB * BlockB[cnt]; }
528  alpha = 1.0;
529 
530  } else {
531 
532  alpha = factorB;
533  mem = denB->gStorage( sector_nelec_up[ikappa], sector_spin_up[ikappa], sector_irrep_up[ikappa], sector_nelec_up[ikappa]+2, TwoSRD, IRD );
534 
535  }
536 
537  double * BlockTup = denT->gStorage(sector_nelec_up[ikappa], sector_spin_up[ikappa], sector_irrep_up[ikappa],sector_nelec_up[ikappa], sector_spin_up[ikappa], sector_irrep_up[ikappa]);
538  double * BlockTdo = denT->gStorage(sector_nelec_up[ikappa]+1,sector_spin_down[ikappa],ID, sector_nelec_up[ikappa]+2,TwoSRD, IRD);
539 
540  char notr = 'N';
541  double beta = 0.0; //set
542  dgemm_(&notr,&notr,&dimLU,&dimRD,&dimRU,&alpha,BlockTup,&dimLU,mem,&dimRU,&beta,workmem2,&dimLU);
543 
544  alpha = 1.0;
545  beta = 1.0; //add
546  char trans = 'T';
547  dgemm_(&notr,&trans,&dimLU,&dimLD,&dimRD,&alpha,workmem2,&dimLU,BlockTdo,&dimLD,&beta,storage+kappa2index[ikappa],&dimLU);
548 
549  }
550  }
551 
552  //case 2
553  const int IRU = Irreps::directProd(sector_irrep_up[ikappa],bk_up->gIrrep(index));
554  for (int TwoSRU=sector_spin_up[ikappa]-1; TwoSRU<=sector_spin_up[ikappa]+1; TwoSRU+=2){
555  int dimRU = bk_up->gCurrentDim(index+1, sector_nelec_up[ikappa]+1, TwoSRU, IRU);
556  int dimRD = bk_up->gCurrentDim(index+1, sector_nelec_up[ikappa]+3, sector_spin_down[ikappa], ID);
557  if ((dimRU>0) && (dimRD>0)){
558 
559  int fase = ((((sector_spin_up[ikappa] + sector_spin_down[ikappa] + 1)/2)%2)!=0)?-1:1;
560  double factorB = fase * sqrt(3.0*(TwoSRU+1)) * Wigner::wigner6j( 1, 1, 2, TwoSRU, sector_spin_down[ikappa], sector_spin_up[ikappa] );
561 
562  double alpha;
563  double * mem;
564 
565  if (TwoSRU == sector_spin_down[ikappa]){
566 
567  double factorA = - sqrt(0.5);
568 
569  double * BlockA = denA->gStorage( sector_nelec_up[ikappa]+1, TwoSRU, IRU, sector_nelec_up[ikappa]+3, sector_spin_down[ikappa], ID );
570  double * BlockB = denB->gStorage( sector_nelec_up[ikappa]+1, TwoSRU, IRU, sector_nelec_up[ikappa]+3, sector_spin_down[ikappa], ID );
571 
572  mem = workmem;
573  for (int cnt=0; cnt<dimRU*dimRD; cnt++){ mem[cnt] = factorA * BlockA[cnt] + factorB * BlockB[cnt]; }
574  alpha = 1.0;
575 
576  } else {
577 
578  alpha = factorB;
579  mem = denB->gStorage( sector_nelec_up[ikappa]+1, TwoSRU, IRU, sector_nelec_up[ikappa]+3, sector_spin_down[ikappa], ID );
580 
581  }
582 
583  double * BlockTup = denT->gStorage(sector_nelec_up[ikappa], sector_spin_up[ikappa], sector_irrep_up[ikappa], sector_nelec_up[ikappa]+1, TwoSRU, IRU);
584  double * BlockTdo = denT->gStorage(sector_nelec_up[ikappa]+1, sector_spin_down[ikappa], ID, sector_nelec_up[ikappa]+3, sector_spin_down[ikappa], ID);
585 
586  char notr = 'N';
587  double beta = 0.0; //set
588  dgemm_(&notr,&notr,&dimLU,&dimRD,&dimRU,&alpha,BlockTup,&dimLU,mem,&dimRU,&beta,workmem2,&dimLU);
589 
590  alpha = 1.0;
591  beta = 1.0; //add
592  char trans = 'T';
593  dgemm_(&notr,&trans,&dimLU,&dimLD,&dimRD,&alpha,workmem2,&dimLU,BlockTdo,&dimLD,&beta,storage+kappa2index[ikappa],&dimLU);
594 
595  }
596  }
597  }
598 
599 }
600 
601 void CheMPS2::TensorQ::AddTermsCD(TensorOperator * denC, TensorOperator * denD, TensorT * denT, double * workmem, double * workmem2){
602 
603  if (moving_right){ AddTermsCDRight(denC, denD, denT, workmem, workmem2); }
604  else{ AddTermsCDLeft(denC, denD, denT, workmem, workmem2); }
605 
606 }
607 
608 void CheMPS2::TensorQ::AddTermsCDRight(TensorOperator * denC, TensorOperator * denD, TensorT * denT, double * workmem, double * workmem2){
609 
610  for (int ikappa=0; ikappa<nKappa; ikappa++){
611  const int IRD = Irreps::directProd( sector_irrep_up[ikappa], n_irrep );
612  int dimRU = bk_up->gCurrentDim(index, sector_nelec_up[ikappa], sector_spin_up[ikappa], sector_irrep_up[ikappa]);
613  int dimRD = bk_up->gCurrentDim(index, sector_nelec_up[ikappa]+1, sector_spin_down[ikappa], IRD);
614 
615  //case 1
616  const int ILD = Irreps::directProd(IRD,bk_up->gIrrep(index-1));
617  for (int TwoSLD=sector_spin_down[ikappa]-1; TwoSLD<=sector_spin_down[ikappa]+1; TwoSLD+=2){
618  int dimLU = bk_up->gCurrentDim(index-1, sector_nelec_up[ikappa], sector_spin_up[ikappa], sector_irrep_up[ikappa]);
619  int dimLD = bk_up->gCurrentDim(index-1, sector_nelec_up[ikappa], TwoSLD, ILD);
620 
621  if ((dimLU>0) && (dimLD>0)){
622 
623  int dimLUxLD = dimLU * dimLD;
624 
625  //first set to D
626  int fase = ((((sector_spin_up[ikappa]+sector_spin_down[ikappa]+1)/2)%2)!=0)?-1:1;
627  double factor = fase * sqrt(3.0*(TwoSLD+1)) * Wigner::wigner6j( 1, 2, 1, sector_spin_up[ikappa], sector_spin_down[ikappa], TwoSLD );
628  double * block = denD->gStorage( sector_nelec_up[ikappa], sector_spin_up[ikappa], sector_irrep_up[ikappa], sector_nelec_up[ikappa], TwoSLD, ILD );
629  for (int cnt=0; cnt<dimLUxLD; cnt++){ workmem[cnt] = factor * block[cnt]; }
630 
631  //add C
632  if (TwoSLD==sector_spin_up[ikappa]){
633  factor = sqrt(0.5);
634  block = denC->gStorage( sector_nelec_up[ikappa], sector_spin_up[ikappa], sector_irrep_up[ikappa], sector_nelec_up[ikappa], TwoSLD, ILD );
635  int inc = 1;
636  daxpy_(&dimLUxLD, &factor, block, &inc, workmem, &inc);
637  }
638 
639  double * BlockTup = denT->gStorage( sector_nelec_up[ikappa], sector_spin_up[ikappa], sector_irrep_up[ikappa], sector_nelec_up[ikappa], sector_spin_up[ikappa], sector_irrep_up[ikappa] );
640  double * BlockTdo = denT->gStorage( sector_nelec_up[ikappa], TwoSLD, ILD, sector_nelec_up[ikappa]+1, sector_spin_down[ikappa], IRD );
641 
642  char trans = 'T';
643  char notr = 'N';
644  double alpha = 1.0;
645  double beta = 0.0;
646  dgemm_(&trans,&notr,&dimRU,&dimLD,&dimLU,&alpha,BlockTup,&dimLU,workmem,&dimLU,&beta,workmem2,&dimRU);
647  beta = 1.0;
648  dgemm_(&notr,&notr,&dimRU,&dimRD,&dimLD,&alpha,workmem2,&dimRU,BlockTdo,&dimLD,&beta,storage+kappa2index[ikappa],&dimRU);
649 
650  }
651  }
652 
653  //case 2
654  const int ILU = Irreps::directProd(sector_irrep_up[ikappa],bk_up->gIrrep(index-1));
655  for (int TwoSLU=sector_spin_up[ikappa]-1; TwoSLU<=sector_spin_up[ikappa]+1; TwoSLU+=2){
656  int dimLU = bk_up->gCurrentDim(index-1, sector_nelec_up[ikappa]-1, TwoSLU, ILU);
657  int dimLD = bk_up->gCurrentDim(index-1, sector_nelec_up[ikappa]-1, sector_spin_down[ikappa], IRD);
658 
659  if ((dimLU>0) && (dimLD>0)){
660 
661  int dimLUxLD = dimLU * dimLD;
662 
663  //first set to D
664  int fase = ((((TwoSLU + sector_spin_down[ikappa])/2)%2)!=0)?-1:1;
665  double factor = fase * sqrt(3.0*(sector_spin_up[ikappa]+1)) * Wigner::wigner6j( 1, 2, 1, sector_spin_down[ikappa], sector_spin_up[ikappa], TwoSLU );
666  double * block = denD->gStorage( sector_nelec_up[ikappa]-1, TwoSLU, ILU, sector_nelec_up[ikappa]-1, sector_spin_down[ikappa], IRD );
667  for (int cnt=0; cnt<dimLUxLD; cnt++){ workmem[cnt] = factor * block[cnt]; }
668 
669  //add C
670  if (TwoSLU==sector_spin_down[ikappa]){
671  fase = ((((sector_spin_down[ikappa]+1-sector_spin_up[ikappa])/2)%2)!=0)?-1:1;
672  factor = fase * sqrt(0.5 * ( sector_spin_up[ikappa]+1.0 ) / ( sector_spin_down[ikappa] + 1.0 ) );
673  block = denC->gStorage( sector_nelec_up[ikappa]-1, TwoSLU, ILU, sector_nelec_up[ikappa]-1, sector_spin_down[ikappa], IRD );
674  int inc = 1;
675  daxpy_(&dimLUxLD, &factor, block, &inc, workmem, &inc);
676  }
677 
678  double * BlockTup = denT->gStorage( sector_nelec_up[ikappa]-1, TwoSLU, ILU, sector_nelec_up[ikappa], sector_spin_up[ikappa], sector_irrep_up[ikappa] );
679  double * BlockTdo = denT->gStorage( sector_nelec_up[ikappa]-1, sector_spin_down[ikappa], IRD, sector_nelec_up[ikappa]+1, sector_spin_down[ikappa], IRD );
680 
681  char trans = 'T';
682  char notr = 'N';
683  double alpha = 1.0;
684  double beta = 0.0;
685  dgemm_(&trans,&notr,&dimRU,&dimLD,&dimLU,&alpha,BlockTup,&dimLU,workmem,&dimLU,&beta,workmem2,&dimRU);
686  beta = 1.0;
687  dgemm_(&notr,&notr,&dimRU,&dimRD,&dimLD,&alpha,workmem2,&dimRU,BlockTdo,&dimLD,&beta,storage+kappa2index[ikappa],&dimRU);
688 
689  }
690  }
691  }
692 
693 }
694 
695 void CheMPS2::TensorQ::AddTermsCDLeft(TensorOperator * denC, TensorOperator * denD, TensorT * denT, double * workmem, double * workmem2){
696 
697  for (int ikappa=0; ikappa<nKappa; ikappa++){
698  const int ILD = Irreps::directProd( sector_irrep_up[ikappa], n_irrep );
699  int dimLU = bk_up->gCurrentDim(index, sector_nelec_up[ikappa], sector_spin_up[ikappa], sector_irrep_up[ikappa]);
700  int dimLD = bk_up->gCurrentDim(index, sector_nelec_up[ikappa]+1, sector_spin_down[ikappa], ILD);
701 
702  //case 1
703  const int IRU = Irreps::directProd(sector_irrep_up[ikappa],bk_up->gIrrep(index));
704  for (int TwoSRU=sector_spin_up[ikappa]-1; TwoSRU<=sector_spin_up[ikappa]+1; TwoSRU+=2){
705  int dimRU = bk_up->gCurrentDim(index+1, sector_nelec_up[ikappa]+1, TwoSRU, IRU);
706  int dimRD = bk_up->gCurrentDim(index+1, sector_nelec_up[ikappa]+1, sector_spin_down[ikappa], ILD);
707 
708  if ((dimRU>0) && (dimRD>0)){
709 
710  int dimRUxRD = dimRU * dimRD;
711 
712  //first set to D
713  int fase = ((((sector_spin_up[ikappa]+TwoSRU+3)/2)%2)!=0)?-1:1;
714  double factor = fase * sqrt(3.0/(sector_spin_down[ikappa]+1.0)) * ( TwoSRU + 1 )
715  * Wigner::wigner6j( 1, 1, 2, TwoSRU, sector_spin_down[ikappa], sector_spin_up[ikappa] );
716  double * block = denD->gStorage( sector_nelec_up[ikappa]+1, TwoSRU, IRU, sector_nelec_up[ikappa]+1, sector_spin_down[ikappa], ILD );
717  for (int cnt=0; cnt<dimRUxRD; cnt++){ workmem[cnt] = factor * block[cnt]; }
718 
719  //add C
720  if (TwoSRU==sector_spin_down[ikappa]){
721  factor = sqrt(0.5);
722  block = denC->gStorage( sector_nelec_up[ikappa]+1, TwoSRU, IRU, sector_nelec_up[ikappa]+1, sector_spin_down[ikappa], ILD );
723  int inc = 1;
724  daxpy_(&dimRUxRD, &factor, block, &inc, workmem, &inc);
725  }
726 
727  double * BlockTup = denT->gStorage( sector_nelec_up[ikappa], sector_spin_up[ikappa], sector_irrep_up[ikappa], sector_nelec_up[ikappa]+1, TwoSRU, IRU );
728  double * BlockTdo = denT->gStorage( sector_nelec_up[ikappa]+1, sector_spin_down[ikappa], ILD, sector_nelec_up[ikappa]+1, sector_spin_down[ikappa], ILD );
729 
730  char notr = 'N';
731  double alpha = 1.0;
732  double beta = 0.0;
733  dgemm_(&notr,&notr,&dimLU,&dimRD,&dimRU,&alpha,BlockTup,&dimLU,workmem,&dimRU,&beta,workmem2,&dimLU);
734  beta = 1.0;
735  char trans = 'T';
736  dgemm_(&notr,&trans,&dimLU,&dimLD,&dimRD,&alpha,workmem2,&dimLU,BlockTdo,&dimLD,&beta,storage+kappa2index[ikappa],&dimLU);
737 
738  }
739  }
740 
741  //case 2
742  const int IRD = Irreps::directProd(ILD,bk_up->gIrrep(index));
743  for (int TwoSRD=sector_spin_down[ikappa]-1; TwoSRD<=sector_spin_down[ikappa]+1; TwoSRD+=2){
744  int dimRU = bk_up->gCurrentDim(index+1, sector_nelec_up[ikappa]+2, sector_spin_up[ikappa], sector_irrep_up[ikappa]);
745  int dimRD = bk_up->gCurrentDim(index+1, sector_nelec_up[ikappa]+2, TwoSRD, IRD);
746 
747  if ((dimRU>0) && (dimRD>0)){
748 
749  int dimRUxRD = dimRU * dimRD;
750 
751  //first set to D
752  int fase = (((TwoSRD+1)%2)!=0)?-1:1;
753  double factor = fase * sqrt(3.0*(TwoSRD+1.0)*(sector_spin_up[ikappa]+1.0)/(sector_spin_down[ikappa]+1.0))
754  * Wigner::wigner6j( 1, 1, 2, sector_spin_up[ikappa], TwoSRD, sector_spin_down[ikappa] );
755  double * block = denD->gStorage( sector_nelec_up[ikappa]+2, sector_spin_up[ikappa], sector_irrep_up[ikappa], sector_nelec_up[ikappa]+2, TwoSRD, IRD );
756  for (int cnt=0; cnt<dimRUxRD; cnt++){ workmem[cnt] = factor * block[cnt]; }
757 
758  //add C
759  if (TwoSRD==sector_spin_up[ikappa]){
760  fase = ((((sector_spin_up[ikappa]+1-sector_spin_down[ikappa])/2)%2)!=0)?-1:1;
761  factor = fase * sqrt(0.5 * ( sector_spin_up[ikappa]+1.0 ) / ( sector_spin_down[ikappa] + 1.0 ) );
762  block = denC->gStorage( sector_nelec_up[ikappa]+2, sector_spin_up[ikappa], sector_irrep_up[ikappa], sector_nelec_up[ikappa]+2, TwoSRD, IRD );
763  int inc = 1;
764  daxpy_(&dimRUxRD, &factor, block, &inc, workmem, &inc);
765  }
766 
767  double * BlockTup = denT->gStorage( sector_nelec_up[ikappa], sector_spin_up[ikappa], sector_irrep_up[ikappa],sector_nelec_up[ikappa]+2,sector_spin_up[ikappa],sector_irrep_up[ikappa] );
768  double * BlockTdo = denT->gStorage( sector_nelec_up[ikappa]+1,sector_spin_down[ikappa],ILD, sector_nelec_up[ikappa]+2,TwoSRD, IRD);
769 
770  char notr = 'N';
771  double alpha = 1.0;
772  double beta = 0.0;
773  dgemm_(&notr,&notr,&dimLU,&dimRD,&dimRU,&alpha,BlockTup,&dimLU,workmem,&dimRU,&beta,workmem2,&dimLU);
774  beta = 1.0;
775  char trans = 'T';
776  dgemm_(&notr,&trans,&dimLU,&dimLD,&dimRD,&alpha,workmem2,&dimLU,BlockTdo,&dimLD,&beta,storage+kappa2index[ikappa],&dimLU);
777 
778  }
779  }
780  }
781 
782 }
783 
784 
785 
CheMPS2::Wigner::wigner6j
static double wigner6j(const int two_ja, const int two_jb, const int two_jc, const int two_jd, const int two_je, const int two_jf)
Wigner-6j symbol (gsl api)
Definition: Wigner.cpp:294
CheMPS2::Tensor::index
int index
Index of the Tensor object. For TensorT: a site index; for other tensors: a boundary index...
Definition: Tensor.h:75
CheMPS2::TensorOperator::sector_nelec_up
int * sector_nelec_up
The up particle number sector.
Definition: TensorOperator.h:157
CheMPS2::TensorQ::AddTermsL
void AddTermsL(TensorL **Ltensors, TensorT *denT, double *workmem, double *workmem2)
Add terms after update/clear with previous TensorL&#39;s.
Definition: TensorQ.cpp:106
CheMPS2::TensorOperator::sector_spin_down
int * sector_spin_down
The down spin symmetry sector (pointer points to sectorTwoS1 if two_j == 0)
Definition: TensorOperator.h:166
CheMPS2::TensorT::gIndex
int gIndex() const
Get the location index.
Definition: TensorT.cpp:161
CheMPS2::SyBookkeeper::gCurrentDim
int gCurrentDim(const int boundary, const int N, const int TwoS, const int irrep) const
Get the current virtual dimensions.
Definition: SyBookkeeper.cpp:165
CheMPS2::TensorOperator::gStorage
double * gStorage()
Get the pointer to the storage.
Definition: TensorOperator.cpp:117
CheMPS2::Problem::gMxElement
double gMxElement(const int alpha, const int beta, const int gamma, const int delta) const
Get a specific interaction matrix element.
Definition: Problem.cpp:350
CheMPS2::TensorOperator::moving_right
bool moving_right
Whether or not moving right.
Definition: TensorOperator.h:154
CheMPS2::TensorQ::AddTermsCD
void AddTermsCD(TensorOperator *denC, TensorOperator *denD, TensorT *denT, double *workmem, double *workmem2)
Add terms after update/clear with previous C-tensors and D-tensors.
Definition: TensorQ.cpp:601
CheMPS2::TensorOperator::sector_spin_up
int * sector_spin_up
The up spin symmetry sector.
Definition: TensorOperator.h:163
CheMPS2::TensorT::gStorage
double * gStorage()
Get the pointer to the storage.
Definition: TensorT.cpp:134
CheMPS2::TensorQ::~TensorQ
virtual ~TensorQ()
Destructor.
Definition: TensorQ.cpp:43
CheMPS2::TensorOperator::n_irrep
int n_irrep
The (real-valued abelian) point group irrep difference between the symmetry sectors of the lower and ...
Definition: TensorOperator.h:151
CheMPS2::Irreps::directProd
static int directProd(const int Irrep1, const int Irrep2)
Get the direct product of the irreps with numbers Irrep1 and Irrep2: a bitwise XOR for psi4&#39;s convent...
Definition: Irreps.h:123
CheMPS2::TensorQ::AddTermSimple
void AddTermSimple(TensorT *denT)
Add terms after update/clear without previous tensors.
Definition: TensorQ.cpp:45
CheMPS2::Problem::gL
int gL() const
Get the number of orbitals.
Definition: Problem.h:51
CheMPS2::TensorOperator::sector_irrep_up
int * sector_irrep_up
The up spin symmetry sector.
Definition: TensorOperator.h:160
CheMPS2::TensorQ::TensorQ
TensorQ(const int boundary_index, const int Idiff, const bool moving_right, const SyBookkeeper *denBK, const Problem *Prob, const int site)
Constructor.
Definition: TensorQ.cpp:27
CheMPS2::TensorOperator::bk_up
const SyBookkeeper * bk_up
The bookkeeper of the upper MPS.
Definition: TensorOperator.h:139
CheMPS2::Tensor::kappa2index
int * kappa2index
kappa2index[kappa] indicates the start of tensor block kappa in storage. kappa2index[nKappa] gives th...
Definition: Tensor.h:84
CheMPS2::TensorQ::AddTermsAB
void AddTermsAB(TensorOperator *denA, TensorOperator *denB, TensorT *denT, double *workmem, double *workmem2)
Add terms after update/clear with previous A-tensors and B-tensors.
Definition: TensorQ.cpp:386
CheMPS2::TensorOperator::get_irrep
int get_irrep() const
Get the (real-valued abelian) point group irrep difference between the symmetry sectors of the lower ...
Definition: TensorOperator.cpp:155
CheMPS2::Tensor::nKappa
int nKappa
Number of Tensor blocks.
Definition: Tensor.h:81
CheMPS2::SyBookkeeper::gIrrep
int gIrrep(const int orbital) const
Get an orbital irrep.
Definition: SyBookkeeper.cpp:145
CheMPS2::Tensor::storage
double * storage
The actual variables. Tensor block kappa begins at storage+kappa2index[kappa] and ends at storage+kap...
Definition: Tensor.h:78