CheMPS2
TensorX.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 "TensorX.h"
24 #include "Lapack.h"
25 #include "Wigner.h"
26 
27 CheMPS2::TensorX::TensorX(const int boundary_index, const bool moving_right, const SyBookkeeper * denBK, const Problem * Prob) :
28 TensorOperator(boundary_index,
29  0, //two_j
30  0, //n_elec
31  0, //n_irrep
32  moving_right,
33  true, //prime_last (doesn't matter for spin-0 tensors)
34  false, //jw_phase (four 2nd quantized operators)
35  denBK,
36  denBK){
37 
38  this->Prob = Prob;
39 
40 }
41 
42 CheMPS2::TensorX::~TensorX(){ }
43 
44 void CheMPS2::TensorX::update(TensorT * denT){
45 
46  if (moving_right){
47  //PARALLEL
48  #pragma omp parallel for schedule(dynamic)
49  for (int ikappa=0; ikappa<nKappa; ikappa++){ makenewRight(ikappa, denT); }
50  } else {
51  //PARALLEL
52  #pragma omp parallel for schedule(dynamic)
53  for (int ikappa=0; ikappa<nKappa; ikappa++){ makenewLeft(ikappa, denT); }
54  }
55 
56 }
57 
58 void CheMPS2::TensorX::update(TensorT * denT, TensorL ** Ltensors, TensorX * Xtensor, TensorQ * Qtensor, TensorOperator * Atensor, TensorOperator * Ctensor, TensorOperator * Dtensor){
59 
60  if (moving_right){
61  //PARALLEL
62  #pragma omp parallel
63  {
64 
65  const bool doOtherThings = (index>1) ? true : false ;
66  const int dimL = (doOtherThings) ? bk_up->gMaxDimAtBound(index-1) : 0 ;
67  const int dimR = (doOtherThings) ? bk_up->gMaxDimAtBound(index) : 0 ;
68  double * workmemLL = (doOtherThings) ? new double[dimL*dimL] : NULL ;
69  double * workmemLR = (doOtherThings) ? new double[dimL*dimR] : NULL ;
70  double * workmemRR = (doOtherThings) ? new double[dimR*dimR] : NULL ;
71 
72  #pragma omp for schedule(dynamic)
73  for (int ikappa=0; ikappa<nKappa; ikappa++){
74  makenewRight(ikappa, denT);
75  if (doOtherThings){
76  update_moving_right(ikappa, Xtensor, denT, denT, workmemLR);
77  addTermQLRight(ikappa, denT, Ltensors, Qtensor, workmemRR, workmemLR, workmemLL);
78  addTermARight(ikappa, denT, Atensor, workmemRR, workmemLR);
79  addTermCRight(ikappa, denT, Ctensor, workmemLR);
80  addTermDRight(ikappa, denT, Dtensor, workmemLR);
81  }
82  }
83 
84  if (doOtherThings){
85  delete [] workmemLL;
86  delete [] workmemLR;
87  delete [] workmemRR;
88  }
89 
90  }
91  } else {
92  //PARALLEL
93  #pragma omp parallel
94  {
95 
96  const bool doOtherThings = (index<Prob->gL()-1) ? true : false ;
97  const int dimL = (doOtherThings) ? bk_up->gMaxDimAtBound(index) : 0 ;
98  const int dimR = (doOtherThings) ? bk_up->gMaxDimAtBound(index+1) : 0 ;
99  double * workmemLL = (doOtherThings) ? new double[dimL*dimL] : NULL ;
100  double * workmemLR = (doOtherThings) ? new double[dimL*dimR] : NULL ;
101  double * workmemRR = (doOtherThings) ? new double[dimR*dimR] : NULL ;
102 
103  #pragma omp for schedule(dynamic)
104  for (int ikappa=0; ikappa<nKappa; ikappa++){
105  makenewLeft(ikappa, denT);
106  if (doOtherThings){
107  update_moving_left(ikappa, Xtensor, denT, denT, workmemLR);
108  addTermQLLeft(ikappa, denT, Ltensors, Qtensor, workmemLL, workmemLR, workmemRR);
109  addTermALeft(ikappa, denT, Atensor, workmemLR, workmemLL);
110  addTermCLeft(ikappa, denT, Ctensor, workmemLR);
111  addTermDLeft(ikappa, denT, Dtensor, workmemLR);
112  }
113  }
114 
115  if (doOtherThings){
116  delete [] workmemLL;
117  delete [] workmemLR;
118  delete [] workmemRR;
119  }
120 
121  }
122  }
123 
124 }
125 
126 void CheMPS2::TensorX::makenewRight(const int ikappa, TensorT * denT){
127 
128  int dimR = bk_up->gCurrentDim(index, sector_nelec_up[ikappa], sector_spin_up[ikappa], sector_irrep_up[ikappa]);
129  int dimL = bk_up->gCurrentDim(index-1, sector_nelec_up[ikappa]-2, sector_spin_up[ikappa], sector_irrep_up[ikappa]);
130  double alpha = Prob->gMxElement(index-1,index-1,index-1,index-1);
131 
132  if ((dimL>0) && (fabs(alpha)>0.0)){
133 
134  double * BlockT = 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]);
135  char trans = 'T';
136  char notr = 'N';
137  double beta = 0.0; //because there's only 1 term contributing per kappa, we might as well set it i.o. adding
138  dgemm_(&trans,&notr,&dimR,&dimR,&dimL,&alpha,BlockT,&dimL,BlockT,&dimL,&beta,storage+kappa2index[ikappa],&dimR);
139 
140  } else {
141  for (int cnt=kappa2index[ikappa]; cnt<kappa2index[ikappa+1]; cnt++){ storage[cnt] = 0.0; }
142  }
143 
144 }
145 
146 void CheMPS2::TensorX::makenewLeft(const int ikappa, TensorT * denT){
147 
148  int dimL = bk_up->gCurrentDim(index, sector_nelec_up[ikappa], sector_spin_up[ikappa], sector_irrep_up[ikappa]);
149  int dimR = bk_up->gCurrentDim(index+1, sector_nelec_up[ikappa]+2, sector_spin_up[ikappa], sector_irrep_up[ikappa]);
150  double alpha = Prob->gMxElement(index,index,index,index);
151 
152  if ((dimR>0) && (fabs(alpha)>0.0)){
153 
154  double * BlockT = 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]);
155  char trans = 'T';
156  char notr = 'N';
157  double beta = 0.0; //set, not add (only 1 term)
158  dgemm_(&notr,&trans,&dimL,&dimL,&dimR,&alpha,BlockT,&dimL,BlockT,&dimL,&beta,storage+kappa2index[ikappa],&dimL);
159 
160  } else {
161  for (int cnt=kappa2index[ikappa]; cnt<kappa2index[ikappa+1]; cnt++){ storage[cnt] = 0.0; }
162  }
163 
164 }
165 
166 
167 void CheMPS2::TensorX::addTermQLRight(const int ikappa, TensorT * denT, TensorL ** Lprev, TensorQ * Qprev, double * workmemRR, double * workmemLR, double * workmemLL){
168 
169  int dimR = bk_up->gCurrentDim(index, sector_nelec_up[ikappa], sector_spin_up[ikappa], sector_irrep_up[ikappa]);
170  int dimTot = dimR * dimR;
171  for (int cnt=0; cnt<dimTot; cnt++){ workmemRR[cnt] = 0.0; }
172 
173  for (int geval=0; geval<4; geval++){
174  int NLup,TwoSLup,ILup,NLdown,TwoSLdown,ILdown;
175  switch(geval){
176  case 0:
177  NLup = sector_nelec_up[ikappa]-1;
178  TwoSLup = sector_spin_up[ikappa]-1;
179  ILup = Irreps::directProd( sector_irrep_up[ikappa] , bk_up->gIrrep(index-1) );
180  NLdown = sector_nelec_up[ikappa];
181  TwoSLdown = sector_spin_up[ikappa];
182  ILdown = sector_irrep_up[ikappa];
183  break;
184  case 1:
185  NLup = sector_nelec_up[ikappa]-1;
186  TwoSLup = sector_spin_up[ikappa]+1;
187  ILup = Irreps::directProd( sector_irrep_up[ikappa] , bk_up->gIrrep(index-1) );
188  NLdown = sector_nelec_up[ikappa];
189  TwoSLdown = sector_spin_up[ikappa];
190  ILdown = sector_irrep_up[ikappa];
191  break;
192  case 2:
193  NLup = sector_nelec_up[ikappa]-2;
194  TwoSLup = sector_spin_up[ikappa];
195  ILup = sector_irrep_up[ikappa];
196  NLdown = sector_nelec_up[ikappa]-1;
197  TwoSLdown = sector_spin_up[ikappa]-1;
198  ILdown = Irreps::directProd( sector_irrep_up[ikappa] , bk_up->gIrrep(index-1) );
199  break;
200  case 3:
201  NLup = sector_nelec_up[ikappa]-2;
202  TwoSLup = sector_spin_up[ikappa];
203  ILup = sector_irrep_up[ikappa];
204  NLdown = sector_nelec_up[ikappa]-1;
205  TwoSLdown = sector_spin_up[ikappa]+1;
206  ILdown = Irreps::directProd( sector_irrep_up[ikappa] , bk_up->gIrrep(index-1) );
207  break;
208  }
209  int dimLup = bk_up->gCurrentDim(index-1, NLup, TwoSLup, ILup);
210  int dimLdown = bk_up->gCurrentDim(index-1, NLdown, TwoSLdown, ILdown);
211 
212  if ((dimLup>0) && (dimLdown>0)){
213  double * BlockTup = denT->gStorage(NLup, TwoSLup, ILup, sector_nelec_up[ikappa], sector_spin_up[ikappa], sector_irrep_up[ikappa]);
214  double * BlockTdown = denT->gStorage(NLdown, TwoSLdown, ILdown, sector_nelec_up[ikappa], sector_spin_up[ikappa], sector_irrep_up[ikappa]);
215  double * BlockQ = Qprev->gStorage(NLup, TwoSLup, ILup, NLdown, TwoSLdown, ILdown);
216 
217  double factor;
218  double * ptr;
219  if (geval<2){
220 
221  factor = 1.0;
222  ptr = BlockQ;
223 
224  } else {
225 
226  int fase = ((((sector_spin_up[ikappa]+1-TwoSLdown)/2)%2)!=0)?-1:1;
227  factor = fase * sqrt((TwoSLdown + 1.0)/(sector_spin_up[ikappa] + 1.0));
228 
229  int dimLupdown = dimLup * dimLdown;
230  int inc = 1;
231  ptr = workmemLL;
232  dcopy_(&dimLupdown,BlockQ,&inc,ptr,&inc);
233 
234  for (int loca=0; loca<index-1; loca++){
235  if (bk_up->gIrrep(index-1) == bk_up->gIrrep(loca)){
236  double alpha = Prob->gMxElement(loca, index-1, index-1, index-1);
237  double * BlockL = Lprev[index-2-loca]->gStorage(NLup, TwoSLup, ILup, NLdown, TwoSLdown, ILdown);
238  daxpy_(&dimLupdown,&alpha,BlockL,&inc,ptr,&inc);
239  }
240  }
241 
242  }
243 
244  //factor * Tup^T * L --> mem2 //set
245  char trans = 'T';
246  char notr = 'N';
247  double beta = 0.0;
248  dgemm_(&trans,&notr,&dimR,&dimLdown,&dimLup,&factor,BlockTup,&dimLup,ptr,&dimLup,&beta,workmemLR,&dimR);
249 
250  //mem2 * Tdown --> mem //add
251  factor = 1.0;
252  beta = 1.0;
253  dgemm_(&notr,&notr,&dimR,&dimR,&dimLdown,&factor,workmemLR,&dimR,BlockTdown,&dimLdown,&beta,workmemRR,&dimR);
254 
255  }
256  }
257  //mem + mem^T --> storage
258  for (int irow = 0; irow<dimR; irow++){
259  for (int icol = irow; icol<dimR; icol++){
260  workmemRR[irow + dimR * icol] += workmemRR[icol + dimR * irow];
261  workmemRR[icol + dimR * irow] = workmemRR[irow + dimR * icol];
262  }
263  }
264  int inc = 1;
265  double alpha = 1.0;
266  daxpy_(&dimTot,&alpha,workmemRR,&inc,storage + kappa2index[ikappa],&inc);
267 
268 }
269 
270 void CheMPS2::TensorX::addTermQLLeft(const int ikappa, TensorT * denT, TensorL ** Lprev, TensorQ * Qprev, double * workmemLL, double * workmemLR, double * workmemRR){
271 
272  int dimL = bk_up->gCurrentDim(index, sector_nelec_up[ikappa], sector_spin_up[ikappa], sector_irrep_up[ikappa]);
273  int dimTot = dimL * dimL;
274  for (int cnt=0; cnt<dimTot; cnt++){ workmemLL[cnt] = 0.0; }
275 
276  for (int geval=0; geval<4; geval++){
277  int NRup,TwoSRup,IRup,NRdown,TwoSRdown,IRdown;
278  switch(geval){
279  case 0:
280  NRup = sector_nelec_up[ikappa];
281  TwoSRup = sector_spin_up[ikappa];
282  IRup = sector_irrep_up[ikappa];
283  NRdown = sector_nelec_up[ikappa]+1;
284  TwoSRdown = sector_spin_up[ikappa]-1;
285  IRdown = Irreps::directProd( sector_irrep_up[ikappa] , bk_up->gIrrep(index) );
286  break;
287  case 1:
288  NRup = sector_nelec_up[ikappa];
289  TwoSRup = sector_spin_up[ikappa];
290  IRup = sector_irrep_up[ikappa];
291  NRdown = sector_nelec_up[ikappa]+1;
292  TwoSRdown = sector_spin_up[ikappa]+1;
293  IRdown = Irreps::directProd( sector_irrep_up[ikappa] , bk_up->gIrrep(index) );
294  break;
295  case 2:
296  NRup = sector_nelec_up[ikappa]+1;
297  TwoSRup = sector_spin_up[ikappa]-1;
298  IRup = Irreps::directProd( sector_irrep_up[ikappa] , bk_up->gIrrep(index) );
299  NRdown = sector_nelec_up[ikappa]+2;
300  TwoSRdown = sector_spin_up[ikappa];
301  IRdown = sector_irrep_up[ikappa];
302  break;
303  case 3:
304  NRup = sector_nelec_up[ikappa]+1;
305  TwoSRup = sector_spin_up[ikappa]+1;
306  IRup = Irreps::directProd( sector_irrep_up[ikappa] , bk_up->gIrrep(index) );
307  NRdown = sector_nelec_up[ikappa]+2;
308  TwoSRdown = sector_spin_up[ikappa];
309  IRdown = sector_irrep_up[ikappa];
310  break;
311  }
312  int dimRup = bk_up->gCurrentDim(index+1, NRup, TwoSRup, IRup);
313  int dimRdown = bk_up->gCurrentDim(index+1, NRdown, TwoSRdown, IRdown);
314 
315  if ((dimRup>0) && (dimRdown>0)){
316  double * BlockTup = denT->gStorage(sector_nelec_up[ikappa], sector_spin_up[ikappa], sector_irrep_up[ikappa], NRup, TwoSRup, IRup);
317  double * BlockTdown = denT->gStorage(sector_nelec_up[ikappa], sector_spin_up[ikappa], sector_irrep_up[ikappa], NRdown, TwoSRdown, IRdown);
318  double * BlockQ = Qprev->gStorage(NRup, TwoSRup, IRup, NRdown, TwoSRdown, IRdown);
319 
320  double factor;
321  double * ptr;
322  if (geval<2){
323 
324  factor = (TwoSRdown + 1.0)/(sector_spin_up[ikappa] + 1.0);
325  ptr = BlockQ;
326 
327  } else {
328 
329  int fase = ((((sector_spin_up[ikappa]+1-TwoSRup)/2)%2)!=0)?-1:1;
330  factor = fase * sqrt((TwoSRup + 1.0)/(sector_spin_up[ikappa] + 1.0));
331 
332  int dimRupdown = dimRup * dimRdown;
333  ptr = workmemRR;
334  int inc = 1;
335  dcopy_(&dimRupdown,BlockQ,&inc,ptr,&inc);
336 
337  for (int loca=index+1; loca<Prob->gL(); loca++){
338  if (bk_up->gIrrep(index) == bk_up->gIrrep(loca)){
339  double alpha = Prob->gMxElement(index,index,index,loca);
340  double * BlockL = Lprev[loca-index-1]->gStorage(NRup, TwoSRup, IRup, NRdown, TwoSRdown, IRdown);
341  daxpy_(&dimRupdown,&alpha,BlockL,&inc,ptr,&inc);
342  }
343  }
344 
345  }
346 
347  //factor * Tup * L --> mem2 //set
348  char notr = 'N';
349  double beta = 0.0;//set
350  dgemm_(&notr,&notr,&dimL,&dimRdown,&dimRup,&factor,BlockTup,&dimL,ptr,&dimRup,&beta,workmemLR,&dimL);
351 
352  //mem2 * Tdown^T --> mem //add
353  char trans = 'T';
354  factor = 1.0;
355  beta = 1.0;
356  dgemm_(&notr,&trans,&dimL,&dimL,&dimRdown,&factor,workmemLR,&dimL,BlockTdown,&dimL,&beta,workmemLL,&dimL);
357 
358  }
359  }
360  //mem + mem^T --> storage
361  for (int irow = 0; irow<dimL; irow++){
362  for (int icol = irow; icol<dimL; icol++){
363  workmemLL[irow + dimL * icol] += workmemLL[icol + dimL * irow];
364  workmemLL[icol + dimL * irow] = workmemLL[irow + dimL * icol];
365  }
366  }
367  int inc = 1;
368  double alpha = 1.0;
369  daxpy_(&dimTot,&alpha,workmemLL,&inc,storage + kappa2index[ikappa],&inc);
370 
371 }
372 
373 void CheMPS2::TensorX::addTermARight(const int ikappa, TensorT * denT, TensorOperator * Aprev, double * workmemRR, double * workmemLR){
374 
375  int dimR = bk_up->gCurrentDim(index, sector_nelec_up[ikappa], sector_spin_up[ikappa], sector_irrep_up[ikappa]);
376  int dimLup = bk_up->gCurrentDim(index-1, sector_nelec_up[ikappa]-2, sector_spin_up[ikappa], sector_irrep_up[ikappa]);
377  int dimLdown = bk_up->gCurrentDim(index-1, sector_nelec_up[ikappa], sector_spin_up[ikappa], sector_irrep_up[ikappa]);
378 
379  if ((dimLup>0) && (dimLdown>0)){
380 
381  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]);
382  double * BlockTdown = 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]);
383  double * BlockA = Aprev->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]);
384 
385  //factor * Tup^T * A --> mem2 //set
386  char trans = 'T';
387  char notr = 'N';
388  double factor = sqrt(2.0);
389  double beta = 0.0; //set
390  dgemm_(&trans,&notr,&dimR,&dimLdown,&dimLup,&factor,BlockTup,&dimLup,BlockA,&dimLup,&beta,workmemLR,&dimR);
391 
392  //mem2 * Tdown --> mem //set
393  factor = 1.0;
394  dgemm_(&notr,&notr,&dimR,&dimR,&dimLdown,&factor,workmemLR,&dimR,BlockTdown,&dimLdown,&beta,workmemRR,&dimR);
395 
396  //mem + mem^T --> storage
397  for (int irow = 0; irow<dimR; irow++){
398  for (int icol = irow; icol<dimR; icol++){
399  workmemRR[irow + dimR * icol] += workmemRR[icol + dimR * irow];
400  workmemRR[icol + dimR * irow] = workmemRR[irow + dimR * icol];
401  }
402  }
403  int dimTot = dimR * dimR;
404  int inc = 1;
405  double alpha = 1.0;
406  daxpy_(&dimTot,&alpha,workmemRR,&inc,storage + kappa2index[ikappa],&inc);
407 
408  }
409 
410 }
411 
412 void CheMPS2::TensorX::addTermALeft(const int ikappa, TensorT * denT, TensorOperator * Aprev, double * workmemLR, double * workmemLL){
413 
414  int dimL = bk_up->gCurrentDim(index, sector_nelec_up[ikappa], sector_spin_up[ikappa], sector_irrep_up[ikappa]);
415  int dimRup = bk_up->gCurrentDim(index+1, sector_nelec_up[ikappa], sector_spin_up[ikappa], sector_irrep_up[ikappa]);
416  int dimRdown = bk_up->gCurrentDim(index+1, sector_nelec_up[ikappa]+2, sector_spin_up[ikappa], sector_irrep_up[ikappa]);
417 
418  if ((dimRup>0) && (dimRdown>0)){
419 
420  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]);
421  double * BlockTdown = 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]);
422  double * BlockA = Aprev->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]);
423 
424  //factor * Tup * A --> mem2 //set
425  char notr = 'N';
426  double factor = sqrt(2.0);
427  double beta = 0.0; //set
428  dgemm_(&notr,&notr,&dimL,&dimRdown,&dimRup,&factor,BlockTup,&dimL,BlockA,&dimRup,&beta,workmemLR,&dimL);
429 
430  //mem2 * Tdown^T --> mem //set
431  char trans = 'T';
432  factor = 1.0;
433  dgemm_(&notr,&trans,&dimL,&dimL,&dimRdown,&factor,workmemLR,&dimL,BlockTdown,&dimL,&beta,workmemLL,&dimL);
434 
435  //mem + mem^T --> storage
436  for (int irow = 0; irow<dimL; irow++){
437  for (int icol = irow; icol<dimL; icol++){
438  workmemLL[irow + dimL * icol] += workmemLL[icol + dimL * irow];
439  workmemLL[icol + dimL * irow] = workmemLL[irow + dimL * icol];
440  }
441  }
442  int dimTot = dimL * dimL;
443  int inc = 1;
444  double alpha = 1.0;
445  daxpy_(&dimTot,&alpha,workmemLL,&inc,storage + kappa2index[ikappa],&inc);
446 
447  }
448 
449 }
450 
451 void CheMPS2::TensorX::addTermCRight(const int ikappa, TensorT * denT, TensorOperator * denC, double * workmemLR){
452 
453  int dimR = bk_up->gCurrentDim(index, sector_nelec_up[ikappa], sector_spin_up[ikappa], sector_irrep_up[ikappa]);
454  for (int geval=0; geval<3; geval++){
455  int NL, TwoSL, IL;
456  switch(geval){
457  case 0:
458  NL = sector_nelec_up[ikappa]-1;
459  TwoSL = sector_spin_up[ikappa]-1;
460  IL = Irreps::directProd( sector_irrep_up[ikappa], bk_up->gIrrep(index-1) );
461  break;
462  case 1:
463  NL = sector_nelec_up[ikappa]-1;
464  TwoSL = sector_spin_up[ikappa]+1;
465  IL = Irreps::directProd( sector_irrep_up[ikappa], bk_up->gIrrep(index-1) );
466  break;
467  case 2:
468  NL = sector_nelec_up[ikappa]-2;
469  TwoSL = sector_spin_up[ikappa];
470  IL = sector_irrep_up[ikappa];
471  break;
472  }
473  int dimL = bk_up->gCurrentDim(index-1,NL,TwoSL,IL);
474  if (dimL>0){
475 
476  double * BlockC = denC->gStorage(NL,TwoSL,IL,NL,TwoSL,IL);
477  double * BlockT = denT->gStorage(NL,TwoSL,IL,sector_nelec_up[ikappa],sector_spin_up[ikappa],sector_irrep_up[ikappa]);
478 
479  double factor = (geval<2)?sqrt(0.5):sqrt(2.0);
480  double beta = 0.0; //set
481  char totrans = 'T';
482  dgemm_(&totrans, &totrans, &dimR, &dimL, &dimL, &factor, BlockT, &dimL, BlockC, &dimL, &beta, workmemLR, &dimR);
483 
484  totrans = 'N';
485  factor = 1.0;
486  beta = 1.0; //add
487  dgemm_(&totrans, &totrans, &dimR, &dimR, &dimL, &factor, workmemLR, &dimR, BlockT, &dimL, &beta, storage+kappa2index[ikappa], &dimR);
488 
489  }
490  }
491 
492 }
493 
494 void CheMPS2::TensorX::addTermCLeft(const int ikappa, TensorT * denT, TensorOperator * denC, double * workmemLR){
495 
496  int dimL = bk_up->gCurrentDim(index, sector_nelec_up[ikappa], sector_spin_up[ikappa], sector_irrep_up[ikappa]);
497  for (int geval=0; geval<3; geval++){
498  int NR, TwoSR, IR;
499  switch(geval){
500  case 0:
501  NR = sector_nelec_up[ikappa]+1;
502  TwoSR = sector_spin_up[ikappa]-1;
503  IR = Irreps::directProd( sector_irrep_up[ikappa], bk_up->gIrrep(index) );
504  break;
505  case 1:
506  NR = sector_nelec_up[ikappa]+1;
507  TwoSR = sector_spin_up[ikappa]+1;
508  IR = Irreps::directProd( sector_irrep_up[ikappa], bk_up->gIrrep(index) );
509  break;
510  case 2:
511  NR = sector_nelec_up[ikappa]+2;
512  TwoSR = sector_spin_up[ikappa];
513  IR = sector_irrep_up[ikappa];
514  break;
515  }
516  int dimR = bk_up->gCurrentDim(index+1,NR,TwoSR,IR);
517  if (dimR>0){
518 
519  double * BlockC = denC->gStorage(NR,TwoSR,IR,NR,TwoSR,IR);
520  double * BlockT = denT->gStorage(sector_nelec_up[ikappa],sector_spin_up[ikappa],sector_irrep_up[ikappa],NR,TwoSR,IR);
521 
522  double factor = (geval<2)?(sqrt(0.5)*(TwoSR+1.0)/(sector_spin_up[ikappa]+1.0)):sqrt(2.0);
523  double beta = 0.0; //set
524  char trans = 'T';
525  char notr = 'N';
526  dgemm_(&notr, &trans, &dimL, &dimR, &dimR, &factor, BlockT, &dimL, BlockC, &dimR, &beta, workmemLR, &dimL);
527 
528  factor = 1.0;
529  beta = 1.0; //add
530  dgemm_(&notr, &trans, &dimL, &dimL, &dimR, &factor, workmemLR, &dimL, BlockT, &dimL, &beta, storage+kappa2index[ikappa], &dimL);
531 
532  }
533  }
534 
535 }
536 
537 void CheMPS2::TensorX::addTermDRight(const int ikappa, TensorT * denT, TensorOperator * denD, double * workmemLR){
538 
539  int dimR = bk_up->gCurrentDim(index, sector_nelec_up[ikappa], sector_spin_up[ikappa], sector_irrep_up[ikappa]);
540 
541  const int IL = Irreps::directProd( sector_irrep_up[ikappa], bk_up->gIrrep(index-1) );
542  const int NL = sector_nelec_up[ikappa]-1;
543 
544  for (int geval=0; geval<4; geval++){
545  int TwoSLup, TwoSLdown;
546  switch(geval){
547  case 0:
548  TwoSLup = sector_spin_up[ikappa]-1;
549  TwoSLdown = sector_spin_up[ikappa]-1;
550  break;
551  case 1:
552  TwoSLup = sector_spin_up[ikappa]+1;
553  TwoSLdown = sector_spin_up[ikappa]-1;
554  break;
555  case 2:
556  TwoSLup = sector_spin_up[ikappa]-1;
557  TwoSLdown = sector_spin_up[ikappa]+1;
558  break;
559  case 3:
560  TwoSLup = sector_spin_up[ikappa]+1;
561  TwoSLdown = sector_spin_up[ikappa]+1;
562  break;
563  }
564 
565  int dimLup = bk_up->gCurrentDim(index-1,NL,TwoSLup, IL);
566  int dimLdown = bk_up->gCurrentDim(index-1,NL,TwoSLdown,IL);
567 
568  if ((dimLup>0) && (dimLdown>0)){
569 
570  double * BlockD = denD->gStorage(NL,TwoSLdown,IL,NL,TwoSLup,IL);
571  double * BlockTup = denT->gStorage(NL,TwoSLup, IL,sector_nelec_up[ikappa],sector_spin_up[ikappa],sector_irrep_up[ikappa]);
572  double * BlockTdown = (TwoSLup==TwoSLdown)? BlockTup : denT->gStorage(NL,TwoSLdown,IL,sector_nelec_up[ikappa],sector_spin_up[ikappa],sector_irrep_up[ikappa]);
573 
574  int fase = ((((TwoSLdown + sector_spin_up[ikappa] + 1)/2)%2)!=0)?-1:1;
575  double factor = fase * sqrt(3.0 * (TwoSLup+1))
576  * Wigner::wigner6j( 1, 1, 2, TwoSLup, TwoSLdown, sector_spin_up[ikappa] );
577  double beta = 0.0; //set
578  char totrans = 'T';
579  dgemm_(&totrans, &totrans, &dimR, &dimLdown, &dimLup, &factor, BlockTup, &dimLup, BlockD, &dimLdown, &beta, workmemLR, &dimR);
580 
581  totrans = 'N';
582  factor = 1.0;
583  beta = 1.0; //add
584  dgemm_(&totrans, &totrans, &dimR, &dimR, &dimLdown, &factor, workmemLR, &dimR, BlockTdown, &dimLdown, &beta, storage+kappa2index[ikappa], &dimR);
585 
586  }
587  }
588 
589 }
590 
591 void CheMPS2::TensorX::addTermDLeft(const int ikappa, TensorT * denT, TensorOperator * denD, double * workmemLR){
592 
593  int dimL = bk_up->gCurrentDim(index, sector_nelec_up[ikappa], sector_spin_up[ikappa], sector_irrep_up[ikappa]);
594 
595  const int NR = sector_nelec_up[ikappa]+1;
596  const int IR = Irreps::directProd( sector_irrep_up[ikappa], bk_up->gIrrep(index) );
597 
598  for (int geval=0; geval<4; geval++){
599  int TwoSRup, TwoSRdown;
600  switch(geval){
601  case 0:
602  TwoSRup = sector_spin_up[ikappa] - 1;
603  TwoSRdown = sector_spin_up[ikappa] - 1;
604  break;
605  case 1:
606  TwoSRup = sector_spin_up[ikappa] + 1;
607  TwoSRdown = sector_spin_up[ikappa] - 1;
608  break;
609  case 2:
610  TwoSRup = sector_spin_up[ikappa] - 1;
611  TwoSRdown = sector_spin_up[ikappa] + 1;
612  break;
613  case 3:
614  TwoSRup = sector_spin_up[ikappa] + 1;
615  TwoSRdown = sector_spin_up[ikappa] + 1;
616  break;
617  }
618 
619  int dimRup = bk_up->gCurrentDim(index+1,NR,TwoSRup, IR);
620  int dimRdown = bk_up->gCurrentDim(index+1,NR,TwoSRdown,IR);
621 
622  if ((dimRup>0) && (dimRdown>0)){
623 
624  double * BlockD = denD->gStorage(NR,TwoSRdown,IR,NR,TwoSRup,IR);
625  double * BlockTup = denT->gStorage(sector_nelec_up[ikappa],sector_spin_up[ikappa],sector_irrep_up[ikappa],NR,TwoSRup,IR);
626  double * BlockTdown = (TwoSRup == TwoSRdown)? BlockTup : denT->gStorage(sector_nelec_up[ikappa],sector_spin_up[ikappa],sector_irrep_up[ikappa],NR,TwoSRdown,IR);
627 
628  int fase = ((((sector_spin_up[ikappa] + TwoSRdown + 3)/2)%2)!=0)?-1:1;
629  double factor = fase*sqrt(3.0 *(TwoSRup+1))*((TwoSRdown + 1.0)/(sector_spin_up[ikappa]+1.0))
630  * Wigner::wigner6j( 1, 1, 2, TwoSRup, TwoSRdown, sector_spin_up[ikappa] );
631  double beta = 0.0; //set
632  char trans = 'T';
633  char notr = 'N';
634  dgemm_(&notr, &trans, &dimL, &dimRdown, &dimRup, &factor, BlockTup, &dimL, BlockD, &dimRdown, &beta, workmemLR, &dimL);
635 
636  factor = 1.0;
637  beta = 1.0; //add
638  dgemm_(&notr, &trans, &dimL, &dimL, &dimRdown, &factor, workmemLR, &dimL, BlockTdown, &dimL, &beta, storage+kappa2index[ikappa], &dimL);
639 
640  }
641  }
642 
643 }
644 
645 
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::TensorX::update
void update(TensorT *denT, TensorL **Ltensors, TensorX *Xtensor, TensorQ *Qtensor, TensorOperator *Atensor, TensorOperator *Ctensor, TensorOperator *Dtensor)
Clear and add the relevant terms to the TensorX.
Definition: TensorX.cpp:58
CheMPS2::TensorOperator::sector_nelec_up
int * sector_nelec_up
The up particle number sector.
Definition: TensorOperator.h:157
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::TensorX::TensorX
TensorX(const int boundary_index, const bool moving_right, const SyBookkeeper *denBK, const Problem *Prob)
Constructor.
Definition: TensorX.cpp:27
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::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::TensorOperator::update_moving_left
void update_moving_left(const int ikappa, TensorOperator *previous, TensorT *mps_tensor_up, TensorT *mps_tensor_down, double *workmem)
Update moving left.
Definition: TensorOperator.cpp:293
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::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::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::SyBookkeeper::gMaxDimAtBound
int gMaxDimAtBound(const int boundary) const
Get the maximum virtual dimension at a certain boundary.
Definition: SyBookkeeper.cpp:282
CheMPS2::Tensor::nKappa
int nKappa
Number of Tensor blocks.
Definition: Tensor.h:81
CheMPS2::TensorX::~TensorX
virtual ~TensorX()
Destructor.
Definition: TensorX.cpp:42
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
CheMPS2::TensorOperator::update_moving_right
void update_moving_right(const int ikappa, TensorOperator *previous, TensorT *mps_tensor_up, TensorT *mps_tensor_down, double *workmem)
Update moving right.
Definition: TensorOperator.cpp:179