CheMPS2
HeffDiagrams2.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 <math.h>
21 #include <stdlib.h>
22 
23 #include "Heff.h"
24 #include "Lapack.h"
25 #include "MPIchemps2.h"
26 #include "Wigner.h"
27 
28 void CheMPS2::Heff::addDiagram2a1spin0(const int ikappa, double * memS, double * memHeff, const Sobject * denS, TensorOperator **** Atensors, TensorS0 **** S0tensors, double * workspace) const{
29 
30  #ifdef CHEMPS2_MPI_COMPILATION
31  const int MPIRANK = MPIchemps2::mpi_rank();
32  #endif
33 
34  int NL = denS->gNL(ikappa);
35  int TwoSL = denS->gTwoSL(ikappa);
36  int IL = denS->gIL(ikappa);
37 
38  int NR = denS->gNR(ikappa);
39  int TwoSR = denS->gTwoSR(ikappa);
40  int IR = denS->gIR(ikappa);
41 
42  int N1 = denS->gN1(ikappa);
43  int N2 = denS->gN2(ikappa);
44  int TwoJ = denS->gTwoJ(ikappa);
45 
46  int theindex = denS->gIndex();
47  int dimL = denBK->gCurrentDim(theindex ,NL,TwoSL,IL);
48  int dimR = denBK->gCurrentDim(theindex+2,NR,TwoSR,IR);
49 
50  const bool leftSum = ( theindex < Prob->gL()*0.5 )?true:false;
51 
52  if (leftSum){
53 
54  for (int l_alpha=0; l_alpha<theindex; l_alpha++){
55  for (int l_beta=l_alpha; l_beta<theindex; l_beta++){
56 
57  #ifdef CHEMPS2_MPI_COMPILATION
58  if ( MPIchemps2::owner_absigma( l_alpha, l_beta ) == MPIRANK )
59  #endif
60  {
61  int ILdown = Irreps::directProd(IL,S0tensors[theindex-1][l_beta-l_alpha][theindex-1-l_beta]->get_irrep());
62  int IRdown = Irreps::directProd(IR,Atensors[theindex+1][l_beta-l_alpha][theindex+1-l_beta]->get_irrep());
63  int memSkappa = denS->gKappa(NL-2,TwoSL,ILdown,N1,N2,TwoJ,NR-2,TwoSR,IRdown);
64 
65  if (memSkappa!=-1){
66 
67  int dimLdown = denBK->gCurrentDim(theindex ,NL-2,TwoSL,ILdown);
68  int dimRdown = denBK->gCurrentDim(theindex+2,NR-2,TwoSR,IRdown);
69 
70  double * BlockS0 = S0tensors[theindex-1][l_beta-l_alpha][theindex-1-l_beta]->gStorage(NL-2,TwoSL,ILdown,NL,TwoSL,IL);
71  double * BlockA = Atensors[theindex+1][l_beta-l_alpha][theindex+1-l_beta]->gStorage(NR-2,TwoSR,IRdown,NR,TwoSR,IR);
72 
73  char trans = 'T';
74  char notrans = 'N';
75  double alpha = 1.0;
76  double beta = 0.0;
77 
78  dgemm_(&trans,&notrans,&dimL,&dimRdown,&dimLdown,&alpha,BlockS0,&dimLdown,memS+denS->gKappa2index(memSkappa),&dimLdown,&beta,workspace,&dimL);
79 
80  beta = 1.0;
81 
82  dgemm_(&notrans,&notrans,&dimL,&dimR,&dimRdown,&alpha,workspace,&dimL,BlockA,&dimRdown,&beta,memHeff+denS->gKappa2index(ikappa),&dimL);
83 
84  }
85  }
86  }
87  }
88 
89  } else {
90 
91  for (int l_gamma=theindex+2; l_gamma<Prob->gL(); l_gamma++){
92  for (int l_delta=l_gamma; l_delta<Prob->gL(); l_delta++){
93 
94  #ifdef CHEMPS2_MPI_COMPILATION
95  if ( MPIchemps2::owner_absigma( l_gamma, l_delta ) == MPIRANK )
96  #endif
97  {
98  int ILdown = Irreps::directProd(IL,Atensors[theindex-1][l_delta-l_gamma][l_gamma-theindex]->get_irrep());
99  int IRdown = Irreps::directProd(IR,S0tensors[theindex+1][l_delta-l_gamma][l_gamma-theindex-2]->get_irrep());
100  int memSkappa = denS->gKappa(NL-2,TwoSL,ILdown,N1,N2,TwoJ,NR-2,TwoSR,IRdown);
101 
102  if (memSkappa!=-1){
103 
104  int dimLdown = denBK->gCurrentDim(theindex ,NL-2,TwoSL,ILdown);
105  int dimRdown = denBK->gCurrentDim(theindex+2,NR-2,TwoSR,IRdown);
106 
107  double * BlockA = Atensors[theindex-1][l_delta-l_gamma][l_gamma-theindex]->gStorage(NL-2,TwoSL,ILdown,NL,TwoSL,IL);
108  double * BlockS0 = S0tensors[theindex+1][l_delta-l_gamma][l_gamma-theindex-2]->gStorage(NR-2,TwoSR,IRdown,NR,TwoSR,IR);
109 
110  char trans = 'T';
111  char notrans = 'N';
112  double alpha = 1.0;
113  double beta = 0.0;
114 
115  dgemm_(&trans,&notrans,&dimL,&dimRdown,&dimLdown,&alpha,BlockA,&dimLdown,memS+denS->gKappa2index(memSkappa),&dimLdown,&beta,workspace,&dimL);
116 
117  beta = 1.0;
118 
119  dgemm_(&notrans,&notrans,&dimL,&dimR,&dimRdown,&alpha,workspace,&dimL,BlockS0,&dimRdown,&beta,memHeff+denS->gKappa2index(ikappa),&dimL);
120 
121  }
122  }
123  }
124  }
125  }
126 
127 }
128 
129 void CheMPS2::Heff::addDiagram2a2spin0(const int ikappa, double * memS, double * memHeff, const Sobject * denS, TensorOperator **** Atensors, TensorS0 **** S0tensors, double * workspace) const{
130 
131  #ifdef CHEMPS2_MPI_COMPILATION
132  const int MPIRANK = MPIchemps2::mpi_rank();
133  #endif
134 
135  int NL = denS->gNL(ikappa);
136  int TwoSL = denS->gTwoSL(ikappa);
137  int IL = denS->gIL(ikappa);
138 
139  int NR = denS->gNR(ikappa);
140  int TwoSR = denS->gTwoSR(ikappa);
141  int IR = denS->gIR(ikappa);
142 
143  int N1 = denS->gN1(ikappa);
144  int N2 = denS->gN2(ikappa);
145  int TwoJ = denS->gTwoJ(ikappa);
146 
147  int theindex = denS->gIndex();
148  int dimL = denBK->gCurrentDim(theindex ,NL,TwoSL,IL);
149  int dimR = denBK->gCurrentDim(theindex+2,NR,TwoSR,IR);
150 
151  const bool leftSum = ( theindex < Prob->gL()*0.5 )?true:false;
152 
153  if (leftSum){
154 
155  for (int l_alpha=0; l_alpha<theindex; l_alpha++){
156  for (int l_beta=l_alpha; l_beta<theindex; l_beta++){
157 
158  #ifdef CHEMPS2_MPI_COMPILATION
159  if ( MPIchemps2::owner_absigma( l_alpha, l_beta ) == MPIRANK )
160  #endif
161  {
162  int ILdown = Irreps::directProd(IL,S0tensors[theindex-1][l_beta-l_alpha][theindex-1-l_beta]->get_irrep());
163  int IRdown = Irreps::directProd(IR,Atensors[theindex+1][l_beta-l_alpha][theindex+1-l_beta]->get_irrep());
164  int memSkappa = denS->gKappa(NL+2,TwoSL,ILdown,N1,N2,TwoJ,NR+2,TwoSR,IRdown);
165 
166  if (memSkappa!=-1){
167 
168  int dimLdown = denBK->gCurrentDim(theindex ,NL+2,TwoSL,ILdown);
169  int dimRdown = denBK->gCurrentDim(theindex+2,NR+2,TwoSR,IRdown);
170 
171  double * BlockS0 = S0tensors[theindex-1][l_beta-l_alpha][theindex-1-l_beta]->gStorage(NL,TwoSL,IL,NL+2,TwoSL,ILdown);
172  double * BlockA = Atensors[theindex+1][l_beta-l_alpha][theindex+1-l_beta]->gStorage(NR,TwoSR,IR,NR+2,TwoSR,IRdown);
173 
174  char trans = 'T';
175  char notrans = 'N';
176  double alpha = 1.0;
177  double beta = 0.0;
178 
179  dgemm_(&notrans,&notrans,&dimL,&dimRdown,&dimLdown,&alpha,BlockS0,&dimL,memS+denS->gKappa2index(memSkappa),&dimLdown,&beta,workspace,&dimL);
180 
181  beta = 1.0;
182 
183  dgemm_(&notrans,&trans,&dimL,&dimR,&dimRdown,&alpha,workspace,&dimL,BlockA,&dimR,&beta,memHeff+denS->gKappa2index(ikappa),&dimL);
184 
185  }
186  }
187  }
188  }
189 
190  } else {
191 
192  for (int l_gamma=theindex+2; l_gamma<Prob->gL(); l_gamma++){
193  for (int l_delta=l_gamma; l_delta<Prob->gL(); l_delta++){
194 
195  #ifdef CHEMPS2_MPI_COMPILATION
196  if ( MPIchemps2::owner_absigma( l_gamma, l_delta ) == MPIRANK )
197  #endif
198  {
199  int ILdown = Irreps::directProd(IL,Atensors[theindex-1][l_delta-l_gamma][l_gamma-theindex]->get_irrep());
200  int IRdown = Irreps::directProd(IR,S0tensors[theindex+1][l_delta-l_gamma][l_gamma-theindex-2]->get_irrep());
201  int memSkappa = denS->gKappa(NL+2,TwoSL,ILdown,N1,N2,TwoJ,NR+2,TwoSR,IRdown);
202 
203  if (memSkappa!=-1){
204 
205  int dimLdown = denBK->gCurrentDim(theindex ,NL+2,TwoSL,ILdown);
206  int dimRdown = denBK->gCurrentDim(theindex+2,NR+2,TwoSR,IRdown);
207 
208  double * BlockA = Atensors[theindex-1][l_delta-l_gamma][l_gamma-theindex]->gStorage(NL,TwoSL,IL,NL+2,TwoSL,ILdown);
209  double * BlockS0 = S0tensors[theindex+1][l_delta-l_gamma][l_gamma-theindex-2]->gStorage(NR,TwoSR,IR,NR+2,TwoSR,IRdown);
210 
211  char trans = 'T';
212  char notrans = 'N';
213  double alpha = 1.0;
214  double beta = 0.0;
215 
216  dgemm_(&notrans,&notrans,&dimL,&dimRdown,&dimLdown,&alpha,BlockA,&dimL,memS+denS->gKappa2index(memSkappa),&dimLdown,&beta,workspace,&dimL);
217 
218  beta = 1.0;
219 
220  dgemm_(&notrans,&trans,&dimL,&dimR,&dimRdown,&alpha,workspace,&dimL,BlockS0,&dimR,&beta,memHeff+denS->gKappa2index(ikappa),&dimL);
221 
222  }
223  }
224  }
225  }
226  }
227 
228 }
229 
230 void CheMPS2::Heff::addDiagram2a1spin1(const int ikappa, double * memS, double * memHeff, const Sobject * denS, TensorOperator **** Btensors, TensorS1 **** S1tensors, double * workspace) const{
231 
232  #ifdef CHEMPS2_MPI_COMPILATION
233  const int MPIRANK = MPIchemps2::mpi_rank();
234  #endif
235 
236  int NL = denS->gNL(ikappa);
237  int TwoSL = denS->gTwoSL(ikappa);
238  int IL = denS->gIL(ikappa);
239 
240  int NR = denS->gNR(ikappa);
241  int TwoSR = denS->gTwoSR(ikappa);
242  int IR = denS->gIR(ikappa);
243 
244  int N1 = denS->gN1(ikappa);
245  int N2 = denS->gN2(ikappa);
246  int TwoJ = denS->gTwoJ(ikappa);
247 
248  int theindex = denS->gIndex();
249  int dimL = denBK->gCurrentDim(theindex ,NL,TwoSL,IL);
250  int dimR = denBK->gCurrentDim(theindex+2,NR,TwoSR,IR);
251 
252  const bool leftSum = ( theindex < Prob->gL()*0.5 )?true:false;
253 
254  if (leftSum){
255 
256  for (int TwoSLdown=TwoSL-2; TwoSLdown<=TwoSL+2; TwoSLdown+=2){
257  for (int TwoSRdown=TwoSR-2; TwoSRdown<=TwoSR+2; TwoSRdown+=2){
258 
259  if ((TwoSLdown>=0) && (TwoSRdown>=0) && (abs(TwoSLdown-TwoSRdown)<=TwoJ)){
260 
261  int fase = phase(TwoSRdown+TwoSL+TwoJ+2);
262  const double thefactor = fase * sqrt((TwoSR + 1)*(TwoSL + 1.0)) * Wigner::wigner6j( TwoSLdown, TwoSRdown, TwoJ, TwoSR, TwoSL, 2 );
263 
264  for (int l_alpha=0; l_alpha<theindex; l_alpha++){
265  for (int l_beta=l_alpha+1; l_beta<theindex; l_beta++){
266 
267  #ifdef CHEMPS2_MPI_COMPILATION
268  if ( MPIchemps2::owner_absigma( l_alpha, l_beta ) == MPIRANK )
269  #endif
270  {
271  int ILdown = Irreps::directProd(IL,S1tensors[theindex-1][l_beta-l_alpha][theindex-1-l_beta]->get_irrep());
272  int IRdown = Irreps::directProd(IR,Btensors[theindex+1][l_beta-l_alpha][theindex+1-l_beta]->get_irrep());
273  int memSkappa = denS->gKappa(NL-2,TwoSLdown,ILdown,N1,N2,TwoJ,NR-2,TwoSRdown,IRdown);
274 
275  if (memSkappa!=-1){
276 
277  int dimLdown = denBK->gCurrentDim(theindex ,NL-2,TwoSLdown,ILdown);
278  int dimRdown = denBK->gCurrentDim(theindex+2,NR-2,TwoSRdown,IRdown);
279 
280  double * BlockS1 = S1tensors[theindex-1][l_beta-l_alpha][theindex-1-l_beta]->gStorage(NL-2,TwoSLdown,ILdown,NL,TwoSL,IL);
281  double * BlockB = Btensors[theindex+1][l_beta-l_alpha][theindex+1-l_beta]->gStorage(NR-2,TwoSRdown,IRdown,NR,TwoSR,IR);
282 
283  char trans = 'T';
284  char notrans = 'N';
285  double alpha = thefactor;
286  double beta = 0.0;
287 
288  dgemm_(&trans,&notrans,&dimL,&dimRdown,&dimLdown,&alpha,BlockS1,&dimLdown,memS+denS->gKappa2index(memSkappa),&dimLdown,&beta,workspace,&dimL);
289 
290  alpha = beta = 1.0;
291 
292  dgemm_(&notrans,&notrans,&dimL,&dimR,&dimRdown,&alpha,workspace,&dimL,BlockB,&dimRdown,&beta,memHeff+denS->gKappa2index(ikappa),&dimL);
293 
294  }
295  }
296  }
297  }
298  }
299  }
300  }
301 
302  } else {
303 
304  for (int TwoSLdown=TwoSL-2; TwoSLdown<=TwoSL+2; TwoSLdown+=2){
305  for (int TwoSRdown=TwoSR-2; TwoSRdown<=TwoSR+2; TwoSRdown+=2){
306 
307  if ((TwoSLdown>=0) && (TwoSRdown>=0) && (abs(TwoSLdown-TwoSRdown)<=TwoJ)){
308 
309  int fase = phase(TwoSRdown+TwoSL+TwoJ+2);
310  const double thefactor = fase * sqrt((TwoSR + 1)*(TwoSL + 1.0)) * Wigner::wigner6j(TwoSLdown,TwoSRdown,TwoJ,TwoSR,TwoSL,2);
311 
312  for (int l_gamma=theindex+2; l_gamma<Prob->gL(); l_gamma++){
313  for (int l_delta=l_gamma+1; l_delta<Prob->gL(); l_delta++){
314 
315  #ifdef CHEMPS2_MPI_COMPILATION
316  if ( MPIchemps2::owner_absigma( l_gamma, l_delta ) == MPIRANK )
317  #endif
318  {
319  int ILdown = Irreps::directProd(IL,Btensors[theindex-1][l_delta-l_gamma][l_gamma-theindex]->get_irrep());
320  int IRdown = Irreps::directProd(IR,S1tensors[theindex+1][l_delta-l_gamma][l_gamma-theindex-2]->get_irrep());
321  int memSkappa = denS->gKappa(NL-2,TwoSLdown,ILdown,N1,N2,TwoJ,NR-2,TwoSRdown,IRdown);
322 
323  if (memSkappa!=-1){
324 
325  int dimLdown = denBK->gCurrentDim(theindex ,NL-2,TwoSLdown,ILdown);
326  int dimRdown = denBK->gCurrentDim(theindex+2,NR-2,TwoSRdown,IRdown);
327 
328  double * BlockB = Btensors[theindex-1][l_delta-l_gamma][l_gamma-theindex]->gStorage(NL-2,TwoSLdown,ILdown,NL,TwoSL,IL);
329  double * BlockS1 = S1tensors[theindex+1][l_delta-l_gamma][l_gamma-theindex-2]->gStorage(NR-2,TwoSRdown,IRdown,NR,TwoSR,IR);
330 
331  char trans = 'T';
332  char notrans = 'N';
333  double alpha = thefactor;
334  double beta = 0.0;
335 
336  dgemm_(&trans,&notrans,&dimL,&dimRdown,&dimLdown,&alpha,BlockB,&dimLdown,memS+denS->gKappa2index(memSkappa),&dimLdown,&beta,workspace,&dimL);
337 
338  alpha = beta = 1.0;
339 
340  dgemm_(&notrans,&notrans,&dimL,&dimR,&dimRdown,&alpha,workspace,&dimL,BlockS1,&dimRdown,&beta,memHeff+denS->gKappa2index(ikappa),&dimL);
341 
342  }
343  }
344  }
345  }
346  }
347  }
348  }
349  }
350 
351 }
352 
353 void CheMPS2::Heff::addDiagram2a2spin1(const int ikappa, double * memS, double * memHeff, const Sobject * denS, TensorOperator **** Btensors, TensorS1 **** S1tensors, double * workspace) const{
354 
355  #ifdef CHEMPS2_MPI_COMPILATION
356  const int MPIRANK = MPIchemps2::mpi_rank();
357  #endif
358 
359  int NL = denS->gNL(ikappa);
360  int TwoSL = denS->gTwoSL(ikappa);
361  int IL = denS->gIL(ikappa);
362 
363  int NR = denS->gNR(ikappa);
364  int TwoSR = denS->gTwoSR(ikappa);
365  int IR = denS->gIR(ikappa);
366 
367  int N1 = denS->gN1(ikappa);
368  int N2 = denS->gN2(ikappa);
369  int TwoJ = denS->gTwoJ(ikappa);
370 
371  int theindex = denS->gIndex();
372  int dimL = denBK->gCurrentDim(theindex ,NL,TwoSL,IL);
373  int dimR = denBK->gCurrentDim(theindex+2,NR,TwoSR,IR);
374 
375  const bool leftSum = ( theindex < Prob->gL()*0.5 )?true:false;
376 
377  if (leftSum){
378 
379  for (int TwoSLdown=TwoSL-2; TwoSLdown<=TwoSL+2; TwoSLdown+=2){
380  for (int TwoSRdown=TwoSR-2; TwoSRdown<=TwoSR+2; TwoSRdown+=2){
381 
382  if ((TwoSLdown>=0) && (TwoSRdown>=0) && (abs(TwoSLdown-TwoSRdown)<=TwoJ)){
383 
384  int fase = phase(TwoSLdown+TwoSR+TwoJ+2);
385  const double thefactor = fase * sqrt((TwoSRdown + 1)*(TwoSLdown + 1.0)) * Wigner::wigner6j(TwoSLdown,TwoSRdown,TwoJ,TwoSR,TwoSL,2);
386 
387  for (int l_alpha=0; l_alpha<theindex; l_alpha++){
388  for (int l_beta=l_alpha+1; l_beta<theindex; l_beta++){
389 
390  #ifdef CHEMPS2_MPI_COMPILATION
391  if ( MPIchemps2::owner_absigma( l_alpha, l_beta ) == MPIRANK )
392  #endif
393  {
394  int ILdown = Irreps::directProd(IL,S1tensors[theindex-1][l_beta-l_alpha][theindex-1-l_beta]->get_irrep());
395  int IRdown = Irreps::directProd(IR,Btensors[theindex+1][l_beta-l_alpha][theindex+1-l_beta]->get_irrep());
396  int memSkappa = denS->gKappa(NL+2,TwoSLdown,ILdown,N1,N2,TwoJ,NR+2,TwoSRdown,IRdown);
397 
398  if (memSkappa!=-1){
399 
400  int dimLdown = denBK->gCurrentDim(theindex ,NL+2,TwoSLdown,ILdown);
401  int dimRdown = denBK->gCurrentDim(theindex+2,NR+2,TwoSRdown,IRdown);
402 
403  double * BlockS1 = S1tensors[theindex-1][l_beta-l_alpha][theindex-1-l_beta]->gStorage(NL,TwoSL,IL,NL+2,TwoSLdown,ILdown);
404  double * BlockB = Btensors[ theindex+1][l_beta-l_alpha][theindex+1-l_beta]->gStorage(NR,TwoSR,IR,NR+2,TwoSRdown,IRdown);
405 
406  char trans = 'T';
407  char notr = 'N';
408  double alpha = thefactor;
409  double beta = 0.0;
410 
411  dgemm_(&notr,&notr,&dimL,&dimRdown,&dimLdown,&alpha,BlockS1,&dimL,memS+denS->gKappa2index(memSkappa),&dimLdown,&beta,workspace,&dimL);
412 
413  alpha = beta = 1.0;
414 
415  dgemm_(&notr,&trans,&dimL,&dimR,&dimRdown,&alpha,workspace,&dimL,BlockB,&dimR,&beta,memHeff+denS->gKappa2index(ikappa),&dimL);
416 
417  }
418  }
419  }
420  }
421  }
422  }
423  }
424 
425  } else {
426 
427  for (int TwoSLdown=TwoSL-2; TwoSLdown<=TwoSL+2; TwoSLdown+=2){
428  for (int TwoSRdown=TwoSR-2; TwoSRdown<=TwoSR+2; TwoSRdown+=2){
429 
430  if ((TwoSLdown>=0) && (TwoSRdown>=0) && (abs(TwoSLdown-TwoSRdown)<=TwoJ)){
431 
432  int fase = phase(TwoSLdown+TwoSR+TwoJ+2);
433  const double thefactor = fase * sqrt((TwoSRdown + 1)*(TwoSLdown + 1.0)) * Wigner::wigner6j(TwoSLdown,TwoSRdown,TwoJ,TwoSR,TwoSL,2);
434 
435  for (int l_gamma=theindex+2; l_gamma<Prob->gL(); l_gamma++){
436  for (int l_delta=l_gamma+1; l_delta<Prob->gL(); l_delta++){
437 
438  #ifdef CHEMPS2_MPI_COMPILATION
439  if ( MPIchemps2::owner_absigma( l_gamma, l_delta ) == MPIRANK )
440  #endif
441  {
442  int ILdown = Irreps::directProd(IL,Btensors[theindex-1][l_delta-l_gamma][l_gamma-theindex]->get_irrep());
443  int IRdown = Irreps::directProd(IR,S1tensors[theindex+1][l_delta-l_gamma][l_gamma-theindex-2]->get_irrep());
444  int memSkappa = denS->gKappa(NL+2,TwoSLdown,ILdown,N1,N2,TwoJ,NR+2,TwoSRdown,IRdown);
445 
446  if (memSkappa!=-1){
447 
448  int dimLdown = denBK->gCurrentDim(theindex ,NL+2,TwoSLdown,ILdown);
449  int dimRdown = denBK->gCurrentDim(theindex+2,NR+2,TwoSRdown,IRdown);
450 
451  double * BlockB = Btensors[theindex-1][l_delta-l_gamma][l_gamma-theindex]->gStorage(NL,TwoSL,IL,NL+2,TwoSLdown,ILdown);
452  double * BlockS1 = S1tensors[theindex+1][l_delta-l_gamma][l_gamma-theindex-2]->gStorage(NR,TwoSR,IR,NR+2,TwoSRdown,IRdown);
453 
454  char trans = 'T';
455  char notrans = 'N';
456  double alpha = thefactor;
457  double beta = 0.0;
458 
459  dgemm_(&notrans,&notrans,&dimL,&dimRdown,&dimLdown,&alpha,BlockB,&dimL,memS+denS->gKappa2index(memSkappa),&dimLdown,&beta,workspace,&dimL);
460 
461  alpha = beta = 1.0;
462 
463  dgemm_(&notrans,&trans,&dimL,&dimR,&dimRdown,&alpha,workspace,&dimL,BlockS1,&dimR,&beta,memHeff+denS->gKappa2index(ikappa),&dimL);
464 
465  }
466  }
467  }
468  }
469  }
470  }
471  }
472  }
473 
474 }
475 
476 void CheMPS2::Heff::addDiagram2a3spin0(const int ikappa, double * memS, double * memHeff, const Sobject * denS, TensorOperator **** Ctensors, TensorF0 **** F0tensors, double * workspace) const{
477 
478  #ifdef CHEMPS2_MPI_COMPILATION
479  const int MPIRANK = MPIchemps2::mpi_rank();
480  #endif
481 
482  int NL = denS->gNL(ikappa);
483  int TwoSL = denS->gTwoSL(ikappa);
484  int IL = denS->gIL(ikappa);
485 
486  int NR = denS->gNR(ikappa);
487  int TwoSR = denS->gTwoSR(ikappa);
488  int IR = denS->gIR(ikappa);
489 
490  int N1 = denS->gN1(ikappa);
491  int N2 = denS->gN2(ikappa);
492  int TwoJ = denS->gTwoJ(ikappa);
493 
494  int theindex = denS->gIndex();
495  int dimL = denBK->gCurrentDim(theindex ,NL,TwoSL,IL);
496  int dimR = denBK->gCurrentDim(theindex+2,NR,TwoSR,IR);
497 
498  const bool leftSum = ( theindex < Prob->gL()*0.5 )?true:false;
499 
500  if (leftSum){
501 
502  for (int l_gamma=0; l_gamma<theindex; l_gamma++){
503  for (int l_alpha=l_gamma+1; l_alpha<theindex; l_alpha++){
504 
505  #ifdef CHEMPS2_MPI_COMPILATION
506  if ( MPIchemps2::owner_cdf( Prob->gL(), l_gamma, l_alpha ) == MPIRANK )
507  #endif
508  {
509  int ILdown = Irreps::directProd(IL,F0tensors[theindex-1][l_alpha-l_gamma][theindex-1-l_alpha]->get_irrep());
510  int IRdown = Irreps::directProd(IR,Ctensors[theindex+1][l_alpha-l_gamma][theindex+1-l_alpha]->get_irrep());
511  int memSkappa = denS->gKappa(NL,TwoSL,ILdown,N1,N2,TwoJ,NR,TwoSR,IRdown);
512 
513  if (memSkappa!=-1){
514 
515  int dimLdown = denBK->gCurrentDim(theindex ,NL,TwoSL,ILdown);
516  int dimRdown = denBK->gCurrentDim(theindex+2,NR,TwoSR,IRdown);
517 
518  //no transpose
519  double * ptr = Ctensors[theindex+1][l_alpha-l_gamma][theindex+1-l_alpha]->gStorage(NR,TwoSR,IR,NR,TwoSR,IRdown);
520  double * BlockF0 = F0tensors[theindex-1][l_alpha-l_gamma][theindex-1-l_alpha]->gStorage(NL,TwoSL,IL,NL,TwoSL,ILdown);
521 
522  char trans = 'T';
523  char notrans = 'N';
524  double alpha = 1.0;
525  double beta = 0.0;
526 
527  dgemm_(&notrans,&notrans,&dimL,&dimRdown,&dimLdown,&alpha,BlockF0,&dimL,memS+denS->gKappa2index(memSkappa),&dimLdown,&beta,workspace,&dimL);
528 
529  beta = 1.0;
530 
531  dgemm_(&notrans,&trans,&dimL,&dimR,&dimRdown,&alpha,workspace,&dimL,ptr,&dimR,&beta,memHeff+denS->gKappa2index(ikappa),&dimL);
532 
533  }
534  }
535  }
536  }
537 
538  for (int l_alpha=0; l_alpha<theindex; l_alpha++){
539  for (int l_gamma=l_alpha; l_gamma<theindex; l_gamma++){
540 
541  #ifdef CHEMPS2_MPI_COMPILATION
542  if ( MPIchemps2::owner_cdf( Prob->gL(), l_alpha, l_gamma ) == MPIRANK )
543  #endif
544  {
545  int ILdown = Irreps::directProd(IL,F0tensors[theindex-1][l_gamma-l_alpha][theindex-1-l_gamma]->get_irrep());
546  int IRdown = Irreps::directProd(IR,Ctensors[theindex+1][l_gamma-l_alpha][theindex+1-l_gamma]->get_irrep());
547  int memSkappa = denS->gKappa(NL,TwoSL,ILdown,N1,N2,TwoJ,NR,TwoSR,IRdown);
548 
549  if (memSkappa!=-1){
550 
551  int dimLdown = denBK->gCurrentDim(theindex ,NL,TwoSL,ILdown);
552  int dimRdown = denBK->gCurrentDim(theindex+2,NR,TwoSR,IRdown);
553 
554  //transpose
555  double * ptr = Ctensors[theindex+1][l_gamma-l_alpha][theindex+1-l_gamma]->gStorage(NR,TwoSR,IRdown,NR,TwoSR,IR);
556  double * BlockF0 = F0tensors[theindex-1][l_gamma-l_alpha][theindex-1-l_gamma]->gStorage(NL,TwoSL,ILdown,NL,TwoSL,IL);
557 
558  char trans = 'T';
559  char notrans = 'N';
560  double alpha = 1.0;
561  double beta = 0.0;
562 
563  dgemm_(&trans,&notrans,&dimL,&dimRdown,&dimLdown,&alpha,BlockF0,&dimLdown,memS+denS->gKappa2index(memSkappa),&dimLdown,&beta,workspace,&dimL);
564 
565  beta = 1.0;
566 
567  dgemm_(&notrans,&notrans,&dimL,&dimR,&dimRdown,&alpha,workspace,&dimL,ptr,&dimRdown,&beta,memHeff+denS->gKappa2index(ikappa),&dimL);
568 
569  }
570  }
571  }
572  }
573 
574  } else {
575 
576  for (int l_delta=theindex+2; l_delta<Prob->gL(); l_delta++){
577  for (int l_beta=l_delta+1; l_beta<Prob->gL(); l_beta++){
578 
579  #ifdef CHEMPS2_MPI_COMPILATION
580  if ( MPIchemps2::owner_cdf( Prob->gL(), l_delta, l_beta ) == MPIRANK )
581  #endif
582  {
583  int ILdown = Irreps::directProd(IL,Ctensors[theindex-1][l_beta-l_delta][l_delta-theindex]->get_irrep());
584  int IRdown = Irreps::directProd(IR,F0tensors[theindex+1][l_beta-l_delta][l_delta-theindex-2]->get_irrep());
585  int memSkappa = denS->gKappa(NL,TwoSL,ILdown,N1,N2,TwoJ,NR,TwoSR,IRdown);
586 
587  if (memSkappa!=-1){
588 
589  int dimLdown = denBK->gCurrentDim(theindex ,NL,TwoSL,ILdown);
590  int dimRdown = denBK->gCurrentDim(theindex+2,NR,TwoSR,IRdown);
591 
592  //no transpose
593  double * ptr = Ctensors[theindex-1][l_beta-l_delta][l_delta-theindex]->gStorage(NL,TwoSL,IL,NL,TwoSL,ILdown);
594  double * BlockF0 = F0tensors[theindex+1][l_beta-l_delta][l_delta-theindex-2]->gStorage(NR,TwoSR,IR,NR,TwoSR,IRdown);
595 
596  char trans = 'T';
597  char notrans = 'N';
598  double alpha = 1.0;
599  double beta = 0.0;
600 
601  dgemm_(&notrans,&notrans,&dimL,&dimRdown,&dimLdown,&alpha,ptr,&dimL,memS+denS->gKappa2index(memSkappa),&dimLdown,&beta,workspace,&dimL);
602 
603  beta = 1.0;
604 
605  dgemm_(&notrans,&trans,&dimL,&dimR,&dimRdown,&alpha,workspace,&dimL,BlockF0,&dimR,&beta,memHeff+denS->gKappa2index(ikappa),&dimL);
606 
607  }
608  }
609  }
610  }
611 
612  for (int l_beta=theindex+2; l_beta<Prob->gL(); l_beta++){
613  for (int l_delta=l_beta; l_delta<Prob->gL(); l_delta++){
614 
615  #ifdef CHEMPS2_MPI_COMPILATION
616  if ( MPIchemps2::owner_cdf( Prob->gL(), l_beta, l_delta ) == MPIRANK )
617  #endif
618  {
619  int ILdown = Irreps::directProd(IL,Ctensors[theindex-1][l_delta-l_beta][l_beta-theindex]->get_irrep());
620  int IRdown = Irreps::directProd(IR,F0tensors[theindex+1][l_delta-l_beta][l_beta-theindex-2]->get_irrep());
621  int memSkappa = denS->gKappa(NL,TwoSL,ILdown,N1,N2,TwoJ,NR,TwoSR,IRdown);
622 
623  if (memSkappa!=-1){
624 
625  int dimLdown = denBK->gCurrentDim(theindex ,NL,TwoSL,ILdown);
626  int dimRdown = denBK->gCurrentDim(theindex+2,NR,TwoSR,IRdown);
627 
628  //transpose
629  double * ptr = Ctensors[theindex-1][l_delta-l_beta][l_beta-theindex]->gStorage(NL,TwoSL,ILdown,NL,TwoSL,IL);
630  double * BlockF0 = F0tensors[theindex+1][l_delta-l_beta][l_beta-theindex-2]->gStorage(NR,TwoSR,IRdown,NR,TwoSR,IR);
631 
632  char trans = 'T';
633  char notrans = 'N';
634  double alpha = 1.0;
635  double beta = 0.0;
636 
637  dgemm_(&trans,&notrans,&dimL,&dimRdown,&dimLdown,&alpha,ptr,&dimLdown,memS+denS->gKappa2index(memSkappa),&dimLdown,&beta,workspace,&dimL);
638 
639  beta = 1.0;
640 
641  dgemm_(&notrans,&notrans,&dimL,&dimR,&dimRdown,&alpha,workspace,&dimL,BlockF0,&dimRdown,&beta,memHeff+denS->gKappa2index(ikappa),&dimL);
642 
643  }
644  }
645  }
646  }
647  }
648 
649 }
650 
651 void CheMPS2::Heff::addDiagram2a3spin1(const int ikappa, double * memS, double * memHeff, const Sobject * denS, TensorOperator **** Dtensors, TensorF1 **** F1tensors, double * workspace) const{
652 
653  #ifdef CHEMPS2_MPI_COMPILATION
654  const int MPIRANK = MPIchemps2::mpi_rank();
655  #endif
656 
657  int NL = denS->gNL(ikappa);
658  int TwoSL = denS->gTwoSL(ikappa);
659  int IL = denS->gIL(ikappa);
660 
661  int NR = denS->gNR(ikappa);
662  int TwoSR = denS->gTwoSR(ikappa);
663  int IR = denS->gIR(ikappa);
664 
665  int N1 = denS->gN1(ikappa);
666  int N2 = denS->gN2(ikappa);
667  int TwoJ = denS->gTwoJ(ikappa);
668 
669  int theindex = denS->gIndex();
670  int dimL = denBK->gCurrentDim(theindex ,NL,TwoSL,IL);
671  int dimR = denBK->gCurrentDim(theindex+2,NR,TwoSR,IR);
672 
673  const bool leftSum = ( theindex < Prob->gL()*0.5 )?true:false;
674 
675  if (leftSum){
676 
677  for (int TwoSLdown=TwoSL-2; TwoSLdown<=TwoSL+2; TwoSLdown+=2){
678  for (int TwoSRdown=TwoSR-2; TwoSRdown<=TwoSR+2; TwoSRdown+=2){
679 
680  if ((TwoSLdown>=0) && (TwoSRdown>=0) && (abs(TwoSLdown-TwoSRdown)<=TwoJ)){
681 
682  int fase = phase(TwoSLdown+TwoSRdown+TwoJ+2);
683  double prefactor = fase * sqrt((TwoSR + 1)*(TwoSLdown + 1.0)) * Wigner::wigner6j(TwoSLdown,TwoSRdown,TwoJ,TwoSR,TwoSL,2);
684 
685  for (int l_gamma=0; l_gamma<theindex; l_gamma++){
686  for (int l_alpha=l_gamma+1; l_alpha<theindex; l_alpha++){
687 
688  #ifdef CHEMPS2_MPI_COMPILATION
689  if ( MPIchemps2::owner_cdf( Prob->gL(), l_gamma, l_alpha ) == MPIRANK )
690  #endif
691  {
692  int ILdown = Irreps::directProd(IL,F1tensors[theindex-1][l_alpha-l_gamma][theindex-1-l_alpha]->get_irrep());
693  int IRdown = Irreps::directProd(IR,Dtensors[theindex+1][l_alpha-l_gamma][theindex+1-l_alpha]->get_irrep());
694  int memSkappa = denS->gKappa(NL,TwoSLdown,ILdown,N1,N2,TwoJ,NR,TwoSRdown,IRdown);
695 
696  if (memSkappa!=-1){
697 
698  int dimLdown = denBK->gCurrentDim(theindex ,NL,TwoSLdown,ILdown);
699  int dimRdown = denBK->gCurrentDim(theindex+2,NR,TwoSRdown,IRdown);
700 
701  //no transpose
702  double * ptr = Dtensors[theindex+1][l_alpha-l_gamma][theindex+1-l_alpha]->gStorage(NR,TwoSR,IR,NR,TwoSRdown,IRdown);
703  double * BlockF1 = F1tensors[theindex-1][l_alpha-l_gamma][theindex-1-l_alpha]->gStorage(NL,TwoSL,IL,NL,TwoSLdown,ILdown);
704 
705  char trans = 'T';
706  char notr = 'N';
707  double beta = 0.0;
708 
709  dgemm_(&notr,&notr,&dimL,&dimRdown,&dimLdown,&prefactor,BlockF1,&dimL,memS+denS->gKappa2index(memSkappa),&dimLdown,&beta,workspace,&dimL);
710 
711  beta = 1.0;
712 
713  dgemm_(&notr,&trans,&dimL,&dimR,&dimRdown,&beta,workspace,&dimL,ptr,&dimR,&beta,memHeff+denS->gKappa2index(ikappa),&dimL);
714 
715  }
716  }
717  }
718  }
719 
720  fase = phase(TwoSL+TwoSR+TwoJ+2);
721  prefactor = fase * sqrt((TwoSRdown + 1)*(TwoSL + 1.0)) * Wigner::wigner6j(TwoSLdown,TwoSRdown,TwoJ,TwoSR,TwoSL,2);
722 
723  for (int l_alpha=0; l_alpha<theindex; l_alpha++){
724  for (int l_gamma=l_alpha; l_gamma<theindex; l_gamma++){
725 
726  #ifdef CHEMPS2_MPI_COMPILATION
727  if ( MPIchemps2::owner_cdf( Prob->gL(), l_alpha, l_gamma ) == MPIRANK )
728  #endif
729  {
730  int ILdown = Irreps::directProd(IL,F1tensors[theindex-1][l_gamma-l_alpha][theindex-1-l_gamma]->get_irrep());
731  int IRdown = Irreps::directProd(IR,Dtensors[theindex+1][l_gamma-l_alpha][theindex+1-l_gamma]->get_irrep());
732  int memSkappa = denS->gKappa(NL,TwoSLdown,ILdown,N1,N2,TwoJ,NR,TwoSRdown,IRdown);
733 
734  if (memSkappa!=-1){
735 
736  int dimLdown = denBK->gCurrentDim(theindex ,NL,TwoSLdown,ILdown);
737  int dimRdown = denBK->gCurrentDim(theindex+2,NR,TwoSRdown,IRdown);
738 
739  //transpose
740  double * ptr = Dtensors[theindex+1][l_gamma-l_alpha][theindex+1-l_gamma]->gStorage(NR,TwoSRdown,IRdown,NR,TwoSR,IR);
741  double * BlockF1 = F1tensors[theindex-1][l_gamma-l_alpha][theindex-1-l_gamma]->gStorage(NL,TwoSLdown,ILdown,NL,TwoSL,IL);
742 
743  char trans = 'T';
744  char notr = 'N';
745  double beta = 0.0;
746 
747  dgemm_(&trans,&notr,&dimL,&dimRdown,&dimLdown,&prefactor,BlockF1,&dimLdown,memS+denS->gKappa2index(memSkappa),&dimLdown,&beta,workspace,&dimL);
748 
749  beta = 1.0;
750 
751  dgemm_(&notr,&notr,&dimL,&dimR,&dimRdown,&beta,workspace,&dimL,ptr,&dimRdown,&beta,memHeff+denS->gKappa2index(ikappa),&dimL);
752 
753  }
754  }
755  }
756  }
757  }
758  }
759  }
760 
761  } else {
762 
763  for (int TwoSLdown=TwoSL-2; TwoSLdown<=TwoSL+2; TwoSLdown+=2){
764  for (int TwoSRdown=TwoSR-2; TwoSRdown<=TwoSR+2; TwoSRdown+=2){
765 
766  if ((TwoSLdown>=0) && (TwoSRdown>=0) && (abs(TwoSLdown-TwoSRdown)<=TwoJ)){
767 
768  int fase = phase(TwoSLdown+TwoSRdown+TwoJ+2);
769  double prefactor = fase * sqrt((TwoSR + 1)*(TwoSLdown + 1.0)) * Wigner::wigner6j(TwoSLdown,TwoSRdown,TwoJ,TwoSR,TwoSL,2);
770 
771  for (int l_delta=theindex+2; l_delta<Prob->gL(); l_delta++){
772  for (int l_beta=l_delta+1; l_beta<Prob->gL(); l_beta++){
773 
774  #ifdef CHEMPS2_MPI_COMPILATION
775  if ( MPIchemps2::owner_cdf( Prob->gL(), l_delta, l_beta ) == MPIRANK )
776  #endif
777  {
778  int ILdown = Irreps::directProd(IL,Dtensors[theindex-1][l_beta-l_delta][l_delta-theindex]->get_irrep());
779  int IRdown = Irreps::directProd(IR,F1tensors[theindex+1][l_beta-l_delta][l_delta-theindex-2]->get_irrep());
780  int memSkappa = denS->gKappa(NL,TwoSLdown,ILdown,N1,N2,TwoJ,NR,TwoSRdown,IRdown);
781 
782  if (memSkappa!=-1){
783 
784  int dimLdown = denBK->gCurrentDim(theindex ,NL,TwoSLdown,ILdown);
785  int dimRdown = denBK->gCurrentDim(theindex+2,NR,TwoSRdown,IRdown);
786 
787  //no transpose
788  double * ptr = Dtensors[theindex-1][l_beta-l_delta][l_delta-theindex]->gStorage(NL,TwoSL,IL,NL,TwoSLdown,ILdown);
789  double * BlockF1 = F1tensors[theindex+1][l_beta-l_delta][l_delta-theindex-2]->gStorage(NR,TwoSR,IR,NR,TwoSRdown,IRdown);
790 
791  char trans = 'T';
792  char notr = 'N';
793  double beta = 0.0;
794 
795  dgemm_(&notr,&notr,&dimL,&dimRdown,&dimLdown,&prefactor,ptr,&dimL,memS+denS->gKappa2index(memSkappa),&dimLdown,&beta,workspace,&dimL);
796 
797  beta = 1.0;
798 
799  dgemm_(&notr,&trans,&dimL,&dimR,&dimRdown,&beta,workspace,&dimL,BlockF1,&dimR,&beta,memHeff+denS->gKappa2index(ikappa),&dimL);
800 
801  }
802  }
803  }
804  }
805 
806  fase = phase(TwoSL+TwoSR+TwoJ+2);
807  prefactor = fase * sqrt((TwoSRdown + 1)*(TwoSL + 1.0)) * Wigner::wigner6j(TwoSLdown,TwoSRdown,TwoJ,TwoSR,TwoSL,2);
808 
809  for (int l_beta=theindex+2; l_beta<Prob->gL(); l_beta++){
810  for (int l_delta=l_beta; l_delta<Prob->gL(); l_delta++){
811 
812  #ifdef CHEMPS2_MPI_COMPILATION
813  if ( MPIchemps2::owner_cdf( Prob->gL(), l_beta, l_delta ) == MPIRANK )
814  #endif
815  {
816  int ILdown = Irreps::directProd(IL,Dtensors[theindex-1][l_delta-l_beta][l_beta-theindex]->get_irrep());
817  int IRdown = Irreps::directProd(IR,F1tensors[theindex+1][l_delta-l_beta][l_beta-theindex-2]->get_irrep());
818  int memSkappa = denS->gKappa(NL,TwoSLdown,ILdown,N1,N2,TwoJ,NR,TwoSRdown,IRdown);
819 
820  if (memSkappa!=-1){
821 
822  int dimLdown = denBK->gCurrentDim(theindex ,NL,TwoSLdown,ILdown);
823  int dimRdown = denBK->gCurrentDim(theindex+2,NR,TwoSRdown,IRdown);
824 
825  //transpose
826  double * ptr = Dtensors[theindex-1][l_delta-l_beta][l_beta-theindex]->gStorage(NL,TwoSLdown,ILdown,NL,TwoSL,IL);
827  double * BlockF1 = F1tensors[theindex+1][l_delta-l_beta][l_beta-theindex-2]->gStorage(NR,TwoSRdown,IRdown,NR,TwoSR,IR);
828 
829  char trans = 'T';
830  char notr = 'N';
831  double beta = 0.0;
832 
833  dgemm_(&trans,&notr,&dimL,&dimRdown,&dimLdown,&prefactor,ptr,&dimLdown,memS+denS->gKappa2index(memSkappa),&dimLdown,&beta,workspace,&dimL);
834 
835  beta = 1.0;
836 
837  dgemm_(&notr,&notr,&dimL,&dimR,&dimRdown,&beta,workspace,&dimL,BlockF1,&dimRdown,&beta,memHeff+denS->gKappa2index(ikappa),&dimL);
838 
839  }
840  }
841  }
842  }
843  }
844  }
845  }
846  }
847 
848 }
849 
850 void CheMPS2::Heff::addDiagram2b1and2b2(const int ikappa, double * memS, double * memHeff, const Sobject * denS, TensorOperator * Atensor) const{
851 
852  int N1 = denS->gN1(ikappa);
853 
854  if (N1==0){
855 
856  int theindex = denS->gIndex();
857  int NL = denS->gNL(ikappa);
858  int TwoSL = denS->gTwoSL(ikappa);
859  int IL = denS->gIL(ikappa);
860 
861  int dimLdown = denBK->gCurrentDim(theindex,NL-2,TwoSL,IL);
862 
863  if (dimLdown>0){
864 
865  int NR = denS->gNR(ikappa);
866  int TwoSR = denS->gTwoSR(ikappa);
867  int IR = denS->gIR(ikappa);
868 
869  int dimLup = denBK->gCurrentDim(theindex ,NL,TwoSL,IL);
870  int dimR = denBK->gCurrentDim(theindex+2,NR,TwoSR,IR);
871 
872  int memSkappa = denS->gKappa(NL-2,TwoSL,IL,2, denS->gN2(ikappa), denS->gTwoJ(ikappa), NR,TwoSR,IR);
873 
874  if (memSkappa!=-1){
875 
876  double * BlockA = Atensor->gStorage(NL-2,TwoSL,IL,NL,TwoSL,IL);
877  char trans = 'T';
878  char notrans = 'N';
879  double alpha = sqrt(2.0);
880  double beta = 1.0;
881 
882  dgemm_(&trans,&notrans,&dimLup,&dimR,&dimLdown,&alpha,BlockA,&dimLdown,memS+denS->gKappa2index(memSkappa),&dimLdown,&beta,memHeff+denS->gKappa2index(ikappa),&dimLup);
883 
884  }
885  }
886  }
887 
888  if (N1==2){
889 
890  int theindex = denS->gIndex();
891  int NL = denS->gNL(ikappa);
892  int TwoSL = denS->gTwoSL(ikappa);
893  int IL = denS->gIL(ikappa);
894 
895  int dimLdown = denBK->gCurrentDim(theindex,NL+2,TwoSL,IL);
896 
897  if (dimLdown>0){
898 
899  int NR = denS->gNR(ikappa);
900  int TwoSR = denS->gTwoSR(ikappa);
901  int IR = denS->gIR(ikappa);
902 
903  int dimLup = denBK->gCurrentDim(theindex ,NL,TwoSL,IL);
904  int dimR = denBK->gCurrentDim(theindex+2,NR,TwoSR,IR);
905 
906  int memSkappa = denS->gKappa(NL+2,TwoSL,IL,0, denS->gN2(ikappa), denS->gTwoJ(ikappa), NR,TwoSR,IR);
907 
908  if (memSkappa!=-1){
909 
910  double * BlockA = Atensor->gStorage(NL,TwoSL,IL,NL+2,TwoSL,IL);
911  char notrans = 'N';
912  double alpha = sqrt(2.0);
913  double beta = 1.0;
914 
915  dgemm_(&notrans,&notrans,&dimLup,&dimR,&dimLdown,&alpha,BlockA,&dimLup,memS+denS->gKappa2index(memSkappa),&dimLdown,&beta,memHeff+denS->gKappa2index(ikappa),&dimLup);
916 
917  }
918  }
919  }
920 
921 }
922 
923 void CheMPS2::Heff::addDiagram2c1and2c2(const int ikappa, double * memS, double * memHeff, const Sobject * denS, TensorOperator * Atensor) const{
924 
925  int N2 = denS->gN2(ikappa);
926 
927  if (N2==0){
928 
929  int theindex = denS->gIndex();
930  int NL = denS->gNL(ikappa);
931  int TwoSL = denS->gTwoSL(ikappa);
932  int IL = denS->gIL(ikappa);
933 
934  int dimLdown = denBK->gCurrentDim(theindex,NL-2,TwoSL,IL);
935 
936  if (dimLdown>0){
937 
938  int NR = denS->gNR(ikappa);
939  int TwoSR = denS->gTwoSR(ikappa);
940  int IR = denS->gIR(ikappa);
941 
942  int dimLup = denBK->gCurrentDim(theindex ,NL,TwoSL,IL);
943  int dimR = denBK->gCurrentDim(theindex+2,NR,TwoSR,IR);
944 
945  int memSkappa = denS->gKappa(NL-2,TwoSL,IL, denS->gN1(ikappa), 2, denS->gTwoJ(ikappa), NR,TwoSR,IR);
946 
947  if (memSkappa!=-1){
948 
949  double * BlockA = Atensor->gStorage(NL-2,TwoSL,IL,NL,TwoSL,IL);
950  char trans = 'T';
951  char notrans = 'N';
952  double alpha = sqrt(2.0);
953  double beta = 1.0;
954 
955  dgemm_(&trans,&notrans,&dimLup,&dimR,&dimLdown,&alpha,BlockA,&dimLdown,memS+denS->gKappa2index(memSkappa),&dimLdown,&beta,memHeff+denS->gKappa2index(ikappa),&dimLup);
956 
957  }
958  }
959  }
960 
961  if (N2==2){
962 
963  int theindex = denS->gIndex();
964  int NL = denS->gNL(ikappa);
965  int TwoSL = denS->gTwoSL(ikappa);
966  int IL = denS->gIL(ikappa);
967 
968  int dimLdown = denBK->gCurrentDim(theindex,NL+2,TwoSL,IL);
969 
970  if (dimLdown>0){
971 
972  int NR = denS->gNR(ikappa);
973  int TwoSR = denS->gTwoSR(ikappa);
974  int IR = denS->gIR(ikappa);
975 
976  int dimLup = denBK->gCurrentDim(theindex ,NL,TwoSL,IL);
977  int dimR = denBK->gCurrentDim(theindex+2,NR,TwoSR,IR);
978 
979  int memSkappa = denS->gKappa(NL+2,TwoSL,IL, denS->gN1(ikappa), 0, denS->gTwoJ(ikappa), NR,TwoSR,IR);
980 
981  if (memSkappa!=-1){
982 
983  double * BlockA = Atensor->gStorage(NL,TwoSL,IL,NL+2,TwoSL,IL);
984  char notrans = 'N';
985  double alpha = sqrt(2.0);
986  double beta = 1.0;
987 
988  dgemm_(&notrans,&notrans,&dimLup,&dimR,&dimLdown,&alpha,BlockA,&dimLup,memS+denS->gKappa2index(memSkappa),&dimLdown,&beta,memHeff+denS->gKappa2index(ikappa),&dimLup);
989 
990  }
991  }
992  }
993 
994 }
995 
996 void CheMPS2::Heff::addDiagram2dall(const int ikappa, double * memS, double * memHeff, const Sobject * denS) const{
997 
998  const int N1 = denS->gN1(ikappa);
999  const int N2 = denS->gN2(ikappa);
1000  const int theindex = denS->gIndex();
1001  int size = denBK->gCurrentDim(theindex, denS->gNL(ikappa),denS->gTwoSL(ikappa),denS->gIL(ikappa))
1002  * denBK->gCurrentDim(theindex+2,denS->gNR(ikappa),denS->gTwoSR(ikappa),denS->gIR(ikappa));
1003  int inc = 1;
1004 
1005  if ((N1==2) && (N2==0)){ //2d1
1006 
1007  int memSkappa = denS->gKappa(denS->gNL(ikappa),denS->gTwoSL(ikappa),denS->gIL(ikappa), 0, 2, 0, denS->gNR(ikappa),denS->gTwoSR(ikappa),denS->gIR(ikappa));
1008 
1009  if (memSkappa!=-1){
1010  double factor = Prob->gMxElement(theindex, theindex, theindex+1, theindex+1);
1011  daxpy_(&size,&factor,memS+denS->gKappa2index(memSkappa),&inc,memHeff+denS->gKappa2index(ikappa),&inc);
1012  }
1013  }
1014 
1015  if ((N1==0) && (N2==2)){ //2d2
1016 
1017  int memSkappa = denS->gKappa(denS->gNL(ikappa),denS->gTwoSL(ikappa),denS->gIL(ikappa), 2, 0, 0, denS->gNR(ikappa),denS->gTwoSR(ikappa),denS->gIR(ikappa));
1018 
1019  if (memSkappa!=-1){
1020  double factor = Prob->gMxElement(theindex, theindex, theindex+1, theindex+1);
1021  daxpy_(&size,&factor,memS+denS->gKappa2index(memSkappa),&inc,memHeff+denS->gKappa2index(ikappa),&inc);
1022  }
1023  }
1024 
1025  if ((N1==2) && (N2==2)){ //2d3a
1026 
1027  double factor = 4 * Prob->gMxElement(theindex, theindex+1, theindex, theindex+1)
1028  - 2 * Prob->gMxElement(theindex, theindex+1, theindex+1, theindex);
1029  daxpy_(&size,&factor,memS+denS->gKappa2index(ikappa),&inc,memHeff+denS->gKappa2index(ikappa),&inc);
1030 
1031  }
1032 
1033  if ((N1==1) && (N2==1)){ //2d3b
1034 
1035  int fase = (denS->gTwoJ(ikappa) == 0)? 1: -1;
1036  double factor = Prob->gMxElement(theindex, theindex+1, theindex, theindex+1)
1037  + fase * Prob->gMxElement(theindex, theindex+1, theindex+1, theindex);
1038  daxpy_(&size,&factor,memS+denS->gKappa2index(ikappa),&inc,memHeff+denS->gKappa2index(ikappa),&inc);
1039 
1040  }
1041 
1042  if ((N1==2) && (N2==1)){ //2d3c
1043 
1044  double factor = 2 * Prob->gMxElement(theindex, theindex+1, theindex, theindex+1)
1045  - Prob->gMxElement(theindex, theindex+1, theindex+1, theindex);
1046  daxpy_(&size,&factor,memS+denS->gKappa2index(ikappa),&inc,memHeff+denS->gKappa2index(ikappa),&inc);
1047 
1048  }
1049 
1050  if ((N1==1) && (N2==2)){ //2d3d
1051 
1052  double factor = 2 * Prob->gMxElement(theindex, theindex+1, theindex, theindex+1)
1053  - Prob->gMxElement(theindex, theindex+1, theindex+1, theindex);
1054  daxpy_(&size,&factor,memS+denS->gKappa2index(ikappa),&inc,memHeff+denS->gKappa2index(ikappa),&inc);
1055 
1056  }
1057 
1058 }
1059 
1060 void CheMPS2::Heff::addDiagram2e1and2e2(const int ikappa, double * memS, double * memHeff, const Sobject * denS, TensorOperator * Atensor) const{
1061 
1062  int N1 = denS->gN1(ikappa);
1063 
1064  if (N1==2){
1065 
1066  int theindex = denS->gIndex();
1067  int NL = denS->gNL(ikappa);
1068  int TwoSL = denS->gTwoSL(ikappa);
1069  int IL = denS->gIL(ikappa);
1070  int NR = denS->gNR(ikappa);
1071  int TwoSR = denS->gTwoSR(ikappa);
1072  int IR = denS->gIR(ikappa);
1073 
1074  int dimL = denBK->gCurrentDim(theindex, NL, TwoSL,IL);
1075  int dimRdown = denBK->gCurrentDim(theindex+2,NR-2,TwoSR,IR);
1076  int dimRup = denBK->gCurrentDim(theindex+2,NR, TwoSR,IR);
1077 
1078  int memSkappa = denS->gKappa(NL,TwoSL,IL, 0, denS->gN2(ikappa), denS->gTwoJ(ikappa), NR-2,TwoSR,IR);
1079 
1080  if (memSkappa!=-1){
1081 
1082  double * BlockA = Atensor->gStorage(NR-2,TwoSR,IR,NR,TwoSR,IR);
1083  char notrans = 'N';
1084  double alpha = sqrt(2.0);
1085  double beta = 1.0;
1086 
1087  dgemm_(&notrans,&notrans,&dimL,&dimRup,&dimRdown,&alpha,memS+denS->gKappa2index(memSkappa),&dimL,BlockA,&dimRdown,&beta,memHeff+denS->gKappa2index(ikappa),&dimL);
1088 
1089  }
1090  }
1091 
1092  if (N1==0){
1093 
1094  int theindex = denS->gIndex();
1095  int NL = denS->gNL(ikappa);
1096  int TwoSL = denS->gTwoSL(ikappa);
1097  int IL = denS->gIL(ikappa);
1098  int NR = denS->gNR(ikappa);
1099  int TwoSR = denS->gTwoSR(ikappa);
1100  int IR = denS->gIR(ikappa);
1101 
1102  int dimL = denBK->gCurrentDim(theindex, NL, TwoSL,IL);
1103  int dimRdown = denBK->gCurrentDim(theindex+2,NR+2,TwoSR,IR);
1104  int dimRup = denBK->gCurrentDim(theindex+2,NR, TwoSR,IR);
1105 
1106  int memSkappa = denS->gKappa(NL,TwoSL,IL, 2, denS->gN2(ikappa), denS->gTwoJ(ikappa), NR+2,TwoSR,IR);
1107 
1108  if (memSkappa!=-1){
1109 
1110  double * BlockA = Atensor->gStorage(NR,TwoSR,IR,NR+2,TwoSR,IR);
1111  char notrans = 'N';
1112  char trans = 'T';
1113  double alpha = sqrt(2.0);
1114  double beta = 1.0;
1115 
1116  dgemm_(&notrans,&trans,&dimL,&dimRup,&dimRdown,&alpha,memS+denS->gKappa2index(memSkappa),&dimL,BlockA,&dimRup,&beta,memHeff+denS->gKappa2index(ikappa),&dimL);
1117 
1118  }
1119  }
1120 
1121 }
1122 
1123 void CheMPS2::Heff::addDiagram2f1and2f2(const int ikappa, double * memS, double * memHeff, const Sobject * denS, TensorOperator * Atensor) const{
1124 
1125  int N2 = denS->gN2(ikappa);
1126 
1127  if (N2==2){
1128 
1129  int theindex = denS->gIndex();
1130  int NL = denS->gNL(ikappa);
1131  int TwoSL = denS->gTwoSL(ikappa);
1132  int IL = denS->gIL(ikappa);
1133  int NR = denS->gNR(ikappa);
1134  int TwoSR = denS->gTwoSR(ikappa);
1135  int IR = denS->gIR(ikappa);
1136 
1137  int dimL = denBK->gCurrentDim(theindex, NL, TwoSL,IL);
1138  int dimRdown = denBK->gCurrentDim(theindex+2,NR-2,TwoSR,IR);
1139  int dimRup = denBK->gCurrentDim(theindex+2,NR, TwoSR,IR);
1140 
1141  int memSkappa = denS->gKappa(NL,TwoSL,IL, denS->gN1(ikappa), 0, denS->gTwoJ(ikappa), NR-2,TwoSR,IR);
1142 
1143  if (memSkappa!=-1){
1144 
1145  double * BlockA = Atensor->gStorage(NR-2,TwoSR,IR,NR,TwoSR,IR);
1146  char notrans = 'N';
1147  double alpha = sqrt(2.0);
1148  double beta = 1.0;
1149 
1150  dgemm_(&notrans,&notrans,&dimL,&dimRup,&dimRdown,&alpha,memS+denS->gKappa2index(memSkappa),&dimL,BlockA,&dimRdown,&beta,memHeff+denS->gKappa2index(ikappa),&dimL);
1151 
1152  }
1153  }
1154 
1155  if (N2==0){
1156 
1157  int theindex = denS->gIndex();
1158  int NL = denS->gNL(ikappa);
1159  int TwoSL = denS->gTwoSL(ikappa);
1160  int IL = denS->gIL(ikappa);
1161  int NR = denS->gNR(ikappa);
1162  int TwoSR = denS->gTwoSR(ikappa);
1163  int IR = denS->gIR(ikappa);
1164 
1165  int dimL = denBK->gCurrentDim(theindex, NL, TwoSL,IL);
1166  int dimRdown = denBK->gCurrentDim(theindex+2,NR+2,TwoSR,IR);
1167  int dimRup = denBK->gCurrentDim(theindex+2,NR, TwoSR,IR);
1168 
1169  int memSkappa = denS->gKappa(NL,TwoSL,IL, denS->gN1(ikappa), 2, denS->gTwoJ(ikappa), NR+2,TwoSR,IR);
1170 
1171  if (memSkappa!=-1){
1172 
1173  double * BlockA = Atensor->gStorage(NR,TwoSR,IR,NR+2,TwoSR,IR);
1174  char notrans = 'N';
1175  char trans = 'T';
1176  double alpha = sqrt(2.0);
1177  double beta = 1.0;
1178 
1179  dgemm_(&notrans,&trans,&dimL,&dimRup,&dimRdown,&alpha,memS+denS->gKappa2index(memSkappa),&dimL,BlockA,&dimRup,&beta,memHeff+denS->gKappa2index(ikappa),&dimL);
1180 
1181  }
1182  }
1183 
1184 }
1185 
1186 void CheMPS2::Heff::addDiagram2b3spin0(const int ikappa, double * memS, double * memHeff, const Sobject * denS, TensorOperator * Ctensor) const{
1187 
1188  int N1 = denS->gN1(ikappa);
1189 
1190  if (N1!=0){
1191 
1192  int theindex = denS->gIndex();
1193 
1194  int NL = denS->gNL(ikappa);
1195  int TwoSL = denS->gTwoSL(ikappa);
1196  int IL = denS->gIL(ikappa);
1197 
1198  int dimL = denBK->gCurrentDim(theindex, NL, TwoSL, IL);
1199  int dimR = denBK->gCurrentDim(theindex+2,denS->gNR(ikappa),denS->gTwoSR(ikappa),denS->gIR(ikappa));
1200 
1201  double * Cblock = Ctensor->gStorage(NL,TwoSL,IL,NL,TwoSL,IL);
1202 
1203  char trans = 'T';
1204  char notrans = 'N';
1205  double alpha = ((N1==2)?1.0:0.5)*sqrt(2.0);
1206  double beta = 1.0;
1207 
1208  dgemm_(&trans,&notrans,&dimL,&dimR,&dimL,&alpha,Cblock,&dimL,memS+denS->gKappa2index(ikappa),&dimL,&beta,memHeff+denS->gKappa2index(ikappa),&dimL);
1209 
1210  }
1211 
1212 }
1213 
1214 void CheMPS2::Heff::addDiagram2c3spin0(const int ikappa, double * memS, double * memHeff, const Sobject * denS, TensorOperator * Ctensor) const{
1215 
1216  int N2 = denS->gN2(ikappa);
1217 
1218  if (N2!=0){
1219 
1220  int theindex = denS->gIndex();
1221 
1222  int NL = denS->gNL(ikappa);
1223  int TwoSL = denS->gTwoSL(ikappa);
1224  int IL = denS->gIL(ikappa);
1225 
1226  int dimL = denBK->gCurrentDim(theindex, NL, TwoSL, IL);
1227  int dimR = denBK->gCurrentDim(theindex+2,denS->gNR(ikappa),denS->gTwoSR(ikappa),denS->gIR(ikappa));
1228 
1229  double * Cblock = Ctensor->gStorage(NL,TwoSL,IL,NL,TwoSL,IL);
1230 
1231  char trans = 'T';
1232  char notrans = 'N';
1233  double alpha = ((N2==2)?1.0:0.5)*sqrt(2.0);
1234  double beta = 1.0;
1235 
1236  dgemm_(&trans,&notrans,&dimL,&dimR,&dimL,&alpha,Cblock,&dimL,memS+denS->gKappa2index(ikappa),&dimL,&beta,memHeff+denS->gKappa2index(ikappa),&dimL);
1237 
1238  }
1239 
1240 }
1241 
1242 void CheMPS2::Heff::addDiagram2e3spin0(const int ikappa, double * memS, double * memHeff, const Sobject * denS, TensorOperator * Ctensor) const{
1243 
1244  int N1 = denS->gN1(ikappa);
1245 
1246  if (N1!=0){
1247 
1248  int theindex = denS->gIndex();
1249 
1250  int NR = denS->gNR(ikappa);
1251  int TwoSR = denS->gTwoSR(ikappa);
1252  int IR = denS->gIR(ikappa);
1253 
1254  int dimR = denBK->gCurrentDim(theindex+2,NR, TwoSR, IR);
1255  int dimL = denBK->gCurrentDim(theindex ,denS->gNL(ikappa),denS->gTwoSL(ikappa),denS->gIL(ikappa));
1256 
1257  double * Cblock = Ctensor->gStorage(NR,TwoSR,IR,NR,TwoSR,IR);
1258 
1259  char notrans = 'N';
1260  double alpha = ((N1==2)?1.0:0.5)*sqrt(2.0);
1261  double beta = 1.0;
1262 
1263  dgemm_(&notrans,&notrans,&dimL,&dimR,&dimR,&alpha,memS+denS->gKappa2index(ikappa),&dimL,Cblock,&dimR,&beta,memHeff+denS->gKappa2index(ikappa),&dimL);
1264 
1265  }
1266 
1267 }
1268 
1269 void CheMPS2::Heff::addDiagram2f3spin0(const int ikappa, double * memS, double * memHeff, const Sobject * denS, TensorOperator * Ctensor) const{
1270 
1271  int N2 = denS->gN2(ikappa);
1272 
1273  if (N2!=0){
1274 
1275  int theindex = denS->gIndex();
1276 
1277  int NR = denS->gNR(ikappa);
1278  int TwoSR = denS->gTwoSR(ikappa);
1279  int IR = denS->gIR(ikappa);
1280 
1281  int dimR = denBK->gCurrentDim(theindex+2,NR, TwoSR, IR);
1282  int dimL = denBK->gCurrentDim(theindex ,denS->gNL(ikappa),denS->gTwoSL(ikappa),denS->gIL(ikappa));
1283 
1284  double * Cblock = Ctensor->gStorage(NR,TwoSR,IR,NR,TwoSR,IR);
1285 
1286  char notrans = 'N';
1287  double alpha = ((N2==2)?1.0:0.5)*sqrt(2.0);
1288  double beta = 1.0;
1289 
1290  dgemm_(&notrans,&notrans,&dimL,&dimR,&dimR,&alpha,memS+denS->gKappa2index(ikappa),&dimL,Cblock,&dimR,&beta,memHeff+denS->gKappa2index(ikappa),&dimL);
1291 
1292  }
1293 
1294 }
1295 
1296 void CheMPS2::Heff::addDiagram2b3spin1(const int ikappa, double * memS, double * memHeff, const Sobject * denS, TensorOperator * Dtensor) const{
1297 
1298  int N1 = denS->gN1(ikappa);
1299 
1300  if (N1==1){
1301 
1302  int theindex = denS->gIndex();
1303 
1304  int NL = denS->gNL(ikappa);
1305  int TwoSL = denS->gTwoSL(ikappa);
1306  int IL = denS->gIL(ikappa);
1307  int NR = denS->gNR(ikappa);
1308  int TwoSR = denS->gTwoSR(ikappa);
1309  int IR = denS->gIR(ikappa);
1310  int TwoJ = denS->gTwoJ(ikappa);
1311  int N2 = denS->gN2(ikappa);
1312 
1313  int dimLup = denBK->gCurrentDim(theindex, NL,TwoSL,IL);
1314  int dimR = denBK->gCurrentDim(theindex+2,NR,TwoSR,IR);
1315 
1316  for (int TwoSLdown=TwoSL-2; TwoSLdown<=TwoSL+2; TwoSLdown+=2){
1317 
1318  int dimLdown = denBK->gCurrentDim(theindex, NL, TwoSLdown, IL);
1319 
1320  if (dimLdown>0){
1321 
1322  double * Dblock = Dtensor->gStorage(NL,TwoSLdown,IL,NL,TwoSL,IL);
1323 
1324  int TwoS2 = (N2==1)?1:0;
1325  int TwoJstart = ((TwoSR!=TwoSLdown) || (TwoS2==0)) ? 1 + TwoS2 : 0;
1326  for (int TwoJdown=TwoJstart; TwoJdown<=1+TwoS2; TwoJdown+=2){
1327  if (abs(TwoSLdown-TwoSR)<=TwoJdown){
1328 
1329  int memSkappa = denS->gKappa(NL, TwoSLdown, IL, N1, N2, TwoJdown, NR, TwoSR, IR);
1330 
1331  if (memSkappa!=-1){
1332 
1333  int fase = phase(TwoSLdown + TwoSR + TwoJ + TwoS2 + TwoJdown - 1);
1334  double alpha = fase * sqrt(3.0*(TwoJ+1)*(TwoJdown+1)*(TwoSL+1)) * Wigner::wigner6j(TwoJdown,TwoJ,2,1,1,TwoS2) * Wigner::wigner6j(TwoJdown,TwoJ,2,TwoSL,TwoSLdown,TwoSR);
1335  char trans = 'T';
1336  char notra = 'N';
1337  double beta = 1.0;
1338 
1339  dgemm_(&trans,&notra,&dimLup,&dimR,&dimLdown,&alpha,Dblock,&dimLdown,memS+denS->gKappa2index(memSkappa),&dimLdown,&beta,memHeff+denS->gKappa2index(ikappa),&dimLup);
1340 
1341  }
1342  }
1343  }
1344  }
1345  }
1346  }
1347 
1348 }
1349 
1350 void CheMPS2::Heff::addDiagram2c3spin1(const int ikappa, double * memS, double * memHeff, const Sobject * denS, TensorOperator * Dtensor) const{
1351 
1352  int N2 = denS->gN2(ikappa);
1353 
1354  if (N2==1){
1355 
1356  int theindex = denS->gIndex();
1357 
1358  int NL = denS->gNL(ikappa);
1359  int TwoSL = denS->gTwoSL(ikappa);
1360  int IL = denS->gIL(ikappa);
1361  int NR = denS->gNR(ikappa);
1362  int TwoSR = denS->gTwoSR(ikappa);
1363  int IR = denS->gIR(ikappa);
1364  int TwoJ = denS->gTwoJ(ikappa);
1365  int N1 = denS->gN1(ikappa);
1366 
1367  int dimLup = denBK->gCurrentDim(theindex, NL,TwoSL,IL);
1368  int dimR = denBK->gCurrentDim(theindex+2,NR,TwoSR,IR);
1369 
1370  for (int TwoSLdown=TwoSL-2; TwoSLdown<=TwoSL+2; TwoSLdown+=2){
1371 
1372  int dimLdown = denBK->gCurrentDim(theindex, NL, TwoSLdown, IL);
1373 
1374  if (dimLdown>0){
1375 
1376  double * Dblock = Dtensor->gStorage(NL,TwoSLdown,IL,NL,TwoSL,IL);
1377 
1378  int TwoS1 = (N1==1)?1:0;
1379  int TwoJstart = ((TwoSR!=TwoSLdown) || (TwoS1==0)) ? 1 + TwoS1 : 0;
1380  for (int TwoJdown=TwoJstart; TwoJdown<=1+TwoS1; TwoJdown+=2){
1381  if (abs(TwoSLdown-TwoSR)<=TwoJdown){
1382 
1383  int memSkappa = denS->gKappa(NL, TwoSLdown, IL, N1, N2, TwoJdown, NR, TwoSR, IR);
1384 
1385  if (memSkappa!=-1){
1386 
1387  int fase = phase(TwoSLdown + TwoSR + 2*TwoJ + TwoS1 - 1);
1388  double alpha = fase * sqrt(3.0*(TwoJ+1)*(TwoJdown+1)*(TwoSL+1)) * Wigner::wigner6j(TwoJdown,TwoJ,2,1,1,TwoS1) * Wigner::wigner6j(TwoJdown,TwoJ,2,TwoSL,TwoSLdown,TwoSR);
1389  char trans = 'T';
1390  char notra = 'N';
1391  double beta = 1.0;
1392 
1393  dgemm_(&trans,&notra,&dimLup,&dimR,&dimLdown,&alpha,Dblock,&dimLdown,memS+denS->gKappa2index(memSkappa),&dimLdown,&beta,memHeff+denS->gKappa2index(ikappa),&dimLup);
1394 
1395  }
1396  }
1397  }
1398  }
1399  }
1400  }
1401 
1402 }
1403 
1404 void CheMPS2::Heff::addDiagram2e3spin1(const int ikappa, double * memS, double * memHeff, const Sobject * denS, TensorOperator * Dtensor) const{
1405 
1406  int N1 = denS->gN1(ikappa);
1407 
1408  if (N1==1){
1409 
1410  int theindex = denS->gIndex();
1411 
1412  int NL = denS->gNL(ikappa);
1413  int TwoSL = denS->gTwoSL(ikappa);
1414  int IL = denS->gIL(ikappa);
1415  int NR = denS->gNR(ikappa);
1416  int TwoSR = denS->gTwoSR(ikappa);
1417  int IR = denS->gIR(ikappa);
1418  int TwoJ = denS->gTwoJ(ikappa);
1419  int N2 = denS->gN2(ikappa);
1420 
1421  int dimL = denBK->gCurrentDim(theindex, NL,TwoSL,IL);
1422  int dimRup = denBK->gCurrentDim(theindex+2,NR,TwoSR,IR);
1423 
1424  for (int TwoSRdown=TwoSR-2; TwoSRdown<=TwoSR+2; TwoSRdown+=2){
1425 
1426  int dimRdown = denBK->gCurrentDim(theindex+2, NR, TwoSRdown, IR);
1427 
1428  if (dimRdown>0){
1429 
1430  double * Dblock = Dtensor->gStorage(NR,TwoSRdown,IR,NR,TwoSR,IR);
1431 
1432  int TwoS2 = (N2==1)?1:0;
1433  int TwoJstart = ((TwoSRdown!=TwoSL) || (TwoS2==0)) ? 1 + TwoS2 : 0;
1434  for (int TwoJdown=TwoJstart; TwoJdown<=1+TwoS2; TwoJdown+=2){
1435  if (abs(TwoSL-TwoSRdown)<=TwoJdown){
1436 
1437  int memSkappa = denS->gKappa(NL, TwoSL, IL, N1, N2, TwoJdown, NR, TwoSRdown, IR);
1438 
1439  if (memSkappa!=-1){
1440 
1441  int fase = phase(TwoSRdown + TwoSL + 2*TwoJ + TwoS2 + 1);
1442  double alpha = fase * sqrt(3.0*(TwoJ+1)*(TwoJdown+1)*(TwoSRdown+1)) * Wigner::wigner6j(TwoJdown,TwoJ,2,1,1,TwoS2) * Wigner::wigner6j(TwoJdown,TwoJ,2,TwoSR,TwoSRdown,TwoSL);
1443  char notr = 'N';
1444  double beta = 1.0;
1445 
1446  dgemm_(&notr,&notr,&dimL,&dimRup,&dimRdown,&alpha,memS+denS->gKappa2index(memSkappa),&dimL,Dblock,&dimRdown,&beta,memHeff+denS->gKappa2index(ikappa),&dimL);
1447 
1448  }
1449  }
1450  }
1451  }
1452  }
1453  }
1454 
1455 }
1456 
1457 void CheMPS2::Heff::addDiagram2f3spin1(const int ikappa, double * memS, double * memHeff, const Sobject * denS, TensorOperator * Dtensor) const{
1458 
1459  int N2 = denS->gN2(ikappa);
1460 
1461  if (N2==1){
1462 
1463  int theindex = denS->gIndex();
1464 
1465  int NL = denS->gNL(ikappa);
1466  int TwoSL = denS->gTwoSL(ikappa);
1467  int IL = denS->gIL(ikappa);
1468  int NR = denS->gNR(ikappa);
1469  int TwoSR = denS->gTwoSR(ikappa);
1470  int IR = denS->gIR(ikappa);
1471  int TwoJ = denS->gTwoJ(ikappa);
1472  int N1 = denS->gN1(ikappa);
1473 
1474  int dimL = denBK->gCurrentDim(theindex, NL,TwoSL,IL);
1475  int dimRup = denBK->gCurrentDim(theindex+2,NR,TwoSR,IR);
1476 
1477  for (int TwoSRdown=TwoSR-2; TwoSRdown<=TwoSR+2; TwoSRdown+=2){
1478 
1479  int dimRdown = denBK->gCurrentDim(theindex+2, NR, TwoSRdown, IR);
1480 
1481  if (dimRdown>0){
1482 
1483  double * Dblock = Dtensor->gStorage(NR,TwoSRdown,IR,NR,TwoSR,IR);
1484 
1485  int TwoS1 = (N1==1)?1:0;
1486  int TwoJstart = ((TwoSRdown!=TwoSL) || (TwoS1==0)) ? 1 + TwoS1 : 0;
1487  for (int TwoJdown=TwoJstart; TwoJdown<=1+TwoS1; TwoJdown+=2){
1488  if (abs(TwoSL-TwoSRdown)<=TwoJdown){
1489 
1490  int memSkappa = denS->gKappa(NL, TwoSL, IL, N1, N2, TwoJdown, NR, TwoSRdown, IR);
1491 
1492  if (memSkappa!=-1){
1493 
1494  int fase = phase(TwoSRdown + TwoSL + TwoJ + TwoS1 + TwoJdown + 1);
1495  double alpha = fase * sqrt(3.0*(TwoJ+1)*(TwoJdown+1)*(TwoSRdown+1)) * Wigner::wigner6j(TwoJdown,TwoJ,2,1,1,TwoS1) * Wigner::wigner6j(TwoJdown,TwoJ,2,TwoSR,TwoSRdown,TwoSL);
1496  char notr = 'N';
1497  double beta = 1.0;
1498 
1499  dgemm_(&notr,&notr,&dimL,&dimRup,&dimRdown,&alpha,memS+denS->gKappa2index(memSkappa),&dimL,Dblock,&dimRdown,&beta,memHeff+denS->gKappa2index(ikappa),&dimL);
1500 
1501  }
1502  }
1503  }
1504  }
1505  }
1506  }
1507 
1508 }
1509 
1510 
1511 
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::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::MPIchemps2::mpi_rank
static int mpi_rank()
Get the rank of this MPI process.
Definition: MPIchemps2.h:131
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::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::Heff::phase
static int phase(const int TwoTimesPower)
Phase function.
Definition: Heff.h:75