CheMPS2
HeffDiagrams5.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::addDiagram5A(const int ikappa, double * memS, double * memHeff, const Sobject * denS, TensorL ** Lleft, TensorL ** Lright, double * temp, double * temp2) 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 dimLup = denBK->gCurrentDim(theindex ,NL,TwoSL,IL);
48  int dimRup = denBK->gCurrentDim(theindex+2,NR,TwoSR,IR);
49 
50  char trans = 'T';
51  char notrans = 'N';
52  int inc = 1;
53 
54  int IprodMID = Irreps::directProd(denBK->gIrrep(theindex),denBK->gIrrep(theindex+1));
55 
56  //5A1
57  #ifdef CHEMPS2_MPI_COMPILATION
58  if (( MPIchemps2::owner_specific_diagram( Prob->gL(), MPI_CHEMPS2_5A1 ) == MPIRANK ) && (N1==0) && (N2==0)){
59  #else
60  if ((N1==0) && (N2==0)){
61  #endif
62 
63  for (int TwoSLdown=TwoSL-1; TwoSLdown<=TwoSL+1; TwoSLdown+=2){
64  for (int TwoSRdown=TwoSR-1; TwoSRdown<=TwoSR+1; TwoSRdown+=2){
65  for (int TwoJdown=0; TwoJdown<=2; TwoJdown+=2){
66  if ((abs(TwoSLdown-TwoSRdown)<=TwoJdown) && (TwoSLdown>=0) && (TwoSRdown>=0)){
67 
68  int fase = phase(TwoSLdown+TwoSRdown+2);
69  const double factor = fase * sqrt((TwoJdown+1)*(TwoSRdown+1.0)) * Wigner::wigner6j( 1, 1, TwoJdown, TwoSLdown, TwoSRdown, TwoSR );
70 
71  for (int Irrep=0; Irrep<(denBK->getNumberOfIrreps()); Irrep++){
72 
73  int IrrepTimesMid = Irreps::directProd(Irrep,IprodMID);
74 
75  bool leftOK = false;
76  for (int l_alpha=0; l_alpha<theindex; l_alpha++){
77  if (denBK->gIrrep(l_alpha) == Irrep){ leftOK = true; }
78  }
79  bool rightOK = false;
80  for (int l_beta=theindex+2; l_beta<Prob->gL(); l_beta++){
81  if (denBK->gIrrep(l_beta) == IrrepTimesMid){ rightOK = true; }
82  }
83 
84  if ((leftOK) && (rightOK)){
85 
86  for (int l_alpha=0; l_alpha<theindex; l_alpha++){
87  if (denBK->gIrrep(l_alpha) == Irrep){
88 
89  int ILdown = Irreps::directProd(IL, Irrep);
90  int IRdown = Irreps::directProd(IR, IrrepTimesMid);
91 
92  int dimLdown = denBK->gCurrentDim(theindex, NL-1, TwoSLdown, ILdown);
93  int dimRdown = denBK->gCurrentDim(theindex+2, NR+1, TwoSRdown, IRdown);
94 
95  if ((dimLdown>0) && (dimRdown>0)){
96 
97  int size = dimRup * dimRdown;
98  for (int cnt=0; cnt<size; cnt++){ temp[cnt] = 0.0; }
99 
100  for (int l_beta=theindex+2; l_beta<Prob->gL(); l_beta++){
101  if (denBK->gIrrep(l_beta) == IrrepTimesMid){
102  double prefac = factor * ( Prob->gMxElement(l_alpha, l_beta, theindex, theindex+1)
103  + ((TwoJdown==0)?1:-1) * Prob->gMxElement(l_alpha, l_beta, theindex+1, theindex) );
104  double * LblockR = Lright[l_beta-theindex-2]->gStorage(NR,TwoSR,IR,NR+1,TwoSRdown,IRdown);
105  daxpy_(&size,&prefac,LblockR,&inc,temp,&inc);
106  }
107  }
108 
109  double alpha = 1.0;
110  double beta = 0.0; //set
111  int memSkappa = denS->gKappa(NL-1, TwoSLdown, ILdown, 1, 1, TwoJdown, NR+1, TwoSRdown, IRdown);
112 
113  dgemm_(&notrans,&trans,&dimLdown,&dimRup,&dimRdown,&alpha,memS+denS->gKappa2index(memSkappa),&dimLdown,temp,&dimRup,&beta,temp2,&dimLdown);
114 
115  beta = 1.0; //add
116  double * LblockL = Lleft[theindex-1-l_alpha]->gStorage(NL-1,TwoSLdown,ILdown,NL,TwoSL,IL);
117  dgemm_(&trans,&notrans,&dimLup,&dimRup,&dimLdown,&alpha,LblockL,&dimLdown,temp2,&dimLdown,&beta,memHeff+denS->gKappa2index(ikappa),&dimLup);
118 
119  }
120  }
121  }
122  }
123  }
124  }
125  }
126  }
127  }
128  }
129 
130  //5A2
131  #ifdef CHEMPS2_MPI_COMPILATION
132  if (( MPIchemps2::owner_specific_diagram( Prob->gL(), MPI_CHEMPS2_5A2 ) == MPIRANK ) && (N1==1) && (N2==0)){
133  #else
134  if ((N1==1) && (N2==0)){
135  #endif
136 
137  for (int TwoSLdown=TwoSL-1; TwoSLdown<=TwoSL+1; TwoSLdown+=2){
138  for (int TwoSRdown=TwoSR-1; TwoSRdown<=TwoSR+1; TwoSRdown+=2){
139 
140  if (((TwoSL == TwoSRdown) || (TwoSLdown == TwoSR)) && (abs(TwoSLdown-TwoSRdown)<=1) && (TwoSLdown>=0) && (TwoSRdown>=0)){
141 
142  double factor1 = 0.0;
143  double factor2 = 0.0;
144 
145  if (TwoSL == TwoSRdown){ factor2 = phase(TwoSRdown + 1 - TwoSR); }
146  if (TwoSLdown == TwoSR){
147  int fase = phase(TwoSLdown + 1 - TwoSL);
148  factor1 = fase * sqrt(((TwoSRdown+1)*(TwoSL+1.0))/((TwoSR+1)*(TwoSLdown+1.0)));
149  }
150 
151  for (int Irrep=0; Irrep<(denBK->getNumberOfIrreps()); Irrep++){
152 
153  int IrrepTimesMid = Irreps::directProd(Irrep,IprodMID);
154 
155  bool leftOK = false;
156  for (int l_alpha=0; l_alpha<theindex; l_alpha++){
157  if (denBK->gIrrep(l_alpha) == Irrep){ leftOK = true; }
158  }
159  bool rightOK = false;
160  for (int l_beta=theindex+2; l_beta<Prob->gL(); l_beta++){
161  if (denBK->gIrrep(l_beta) == IrrepTimesMid){ rightOK = true; }
162  }
163 
164  if ((leftOK) && (rightOK)){
165 
166  for (int l_alpha=0; l_alpha<theindex; l_alpha++){
167  if (denBK->gIrrep(l_alpha) == Irrep){
168 
169  int ILdown = Irreps::directProd(IL, Irrep);
170  int IRdown = Irreps::directProd(IR, IrrepTimesMid);
171 
172  int dimLdown = denBK->gCurrentDim(theindex, NL-1, TwoSLdown, ILdown);
173  int dimRdown = denBK->gCurrentDim(theindex+2, NR+1, TwoSRdown, IRdown);
174 
175  if ((dimLdown>0) && (dimRdown>0)){
176 
177  int size = dimRup * dimRdown;
178  for (int cnt=0; cnt<size; cnt++){ temp[cnt] = 0.0; }
179 
180  for (int l_beta=theindex+2; l_beta<Prob->gL(); l_beta++){
181  if (denBK->gIrrep(l_beta) == IrrepTimesMid){
182  double prefac = factor1 * Prob->gMxElement(l_alpha, l_beta, theindex, theindex+1)
183  + factor2 * Prob->gMxElement(l_alpha, l_beta, theindex+1, theindex);
184  double * LblockR = Lright[l_beta-theindex-2]->gStorage(NR,TwoSR,IR,NR+1,TwoSRdown,IRdown);
185  daxpy_(&size,&prefac,LblockR,&inc,temp,&inc);
186  }
187  }
188 
189  double alpha = 1.0;
190  double beta = 0.0; //set
191  int memSkappa = denS->gKappa(NL-1, TwoSLdown, ILdown, 2, 1, 1, NR+1, TwoSRdown, IRdown);
192 
193  dgemm_(&notrans,&trans,&dimLdown,&dimRup,&dimRdown,&alpha,memS+denS->gKappa2index(memSkappa),&dimLdown,temp,&dimRup,&beta,temp2,&dimLdown);
194 
195  beta = 1.0; //add
196  double * LblockL = Lleft[theindex-1-l_alpha]->gStorage(NL-1,TwoSLdown,ILdown,NL,TwoSL,IL);
197  dgemm_(&trans,&notrans,&dimLup,&dimRup,&dimLdown,&alpha,LblockL,&dimLdown,temp2,&dimLdown,&beta,memHeff+denS->gKappa2index(ikappa),&dimLup);
198 
199  }
200  }
201  }
202  }
203  }
204  }
205  }
206  }
207  }
208 
209  //5A3
210  #ifdef CHEMPS2_MPI_COMPILATION
211  if (( MPIchemps2::owner_specific_diagram( Prob->gL(), MPI_CHEMPS2_5A3 ) == MPIRANK ) && (N1==0) && (N2==1)){
212  #else
213  if ((N1==0) && (N2==1)){
214  #endif
215 
216  for (int TwoSLdown=TwoSL-1; TwoSLdown<=TwoSL+1; TwoSLdown+=2){
217  for (int TwoSRdown=TwoSR-1; TwoSRdown<=TwoSR+1; TwoSRdown+=2){
218 
219  if (((TwoSL == TwoSRdown) || (TwoSLdown == TwoSR)) && (abs(TwoSLdown-TwoSRdown)<=1) && (TwoSLdown>=0) && (TwoSRdown>=0)){
220 
221  double factor1 = 0.0;
222  double factor2 = 0.0;
223 
224  if (TwoSL == TwoSRdown){ factor1 = phase(TwoSRdown + 1 - TwoSR); }
225  if (TwoSLdown == TwoSR){
226  int fase = phase(TwoSLdown + 1 - TwoSL);
227  factor2 = fase * sqrt(((TwoSRdown+1)*(TwoSL+1.0))/((TwoSR+1)*(TwoSLdown+1.0)));
228  }
229 
230  for (int Irrep=0; Irrep<(denBK->getNumberOfIrreps()); Irrep++){
231 
232  int IrrepTimesMid = Irreps::directProd(Irrep,IprodMID);
233 
234  bool leftOK = false;
235  for (int l_alpha=0; l_alpha<theindex; l_alpha++){
236  if (denBK->gIrrep(l_alpha) == Irrep){ leftOK = true; }
237  }
238  bool rightOK = false;
239  for (int l_beta=theindex+2; l_beta<Prob->gL(); l_beta++){
240  if (denBK->gIrrep(l_beta) == IrrepTimesMid){ rightOK = true; }
241  }
242 
243  if ((leftOK) && (rightOK)){
244 
245  for (int l_alpha=0; l_alpha<theindex; l_alpha++){
246  if (denBK->gIrrep(l_alpha) == Irrep){
247 
248  int ILdown = Irreps::directProd(IL, Irrep);
249  int IRdown = Irreps::directProd(IR, IrrepTimesMid);
250 
251  int dimLdown = denBK->gCurrentDim(theindex, NL-1, TwoSLdown, ILdown);
252  int dimRdown = denBK->gCurrentDim(theindex+2, NR+1, TwoSRdown, IRdown);
253 
254  if ((dimLdown>0) && (dimRdown>0)){
255 
256  int size = dimRup * dimRdown;
257  for (int cnt=0; cnt<size; cnt++){ temp[cnt] = 0.0; }
258 
259  for (int l_beta=theindex+2; l_beta<Prob->gL(); l_beta++){
260  if (denBK->gIrrep(l_beta) == IrrepTimesMid){
261  double prefac = factor1 * Prob->gMxElement(l_alpha, l_beta, theindex, theindex+1)
262  + factor2 * Prob->gMxElement(l_alpha, l_beta, theindex+1, theindex);
263  double * LblockR = Lright[l_beta-theindex-2]->gStorage(NR,TwoSR,IR,NR+1,TwoSRdown,IRdown);
264  daxpy_(&size,&prefac,LblockR,&inc,temp,&inc);
265  }
266  }
267 
268  double alpha = 1.0;
269  double beta = 0.0; //set
270  int memSkappa = denS->gKappa(NL-1, TwoSLdown, ILdown, 1, 2, 1, NR+1, TwoSRdown, IRdown);
271 
272  dgemm_(&notrans,&trans,&dimLdown,&dimRup,&dimRdown,&alpha,memS+denS->gKappa2index(memSkappa),&dimLdown,temp,&dimRup,&beta,temp2,&dimLdown);
273 
274  beta = 1.0; //add
275  double * LblockL = Lleft[theindex-1-l_alpha]->gStorage(NL-1,TwoSLdown,ILdown,NL,TwoSL,IL);
276  dgemm_(&trans,&notrans,&dimLup,&dimRup,&dimLdown,&alpha,LblockL,&dimLdown,temp2,&dimLdown,&beta,memHeff+denS->gKappa2index(ikappa),&dimLup);
277 
278  }
279  }
280  }
281  }
282  }
283  }
284  }
285  }
286  }
287 
288  //5A4
289  #ifdef CHEMPS2_MPI_COMPILATION
290  if (( MPIchemps2::owner_specific_diagram( Prob->gL(), MPI_CHEMPS2_5A4 ) == MPIRANK ) && (N1==1) && (N2==1)){
291  #else
292  if ((N1==1) && (N2==1)){
293  #endif
294 
295  for (int TwoSLdown=TwoSL-1; TwoSLdown<=TwoSL+1; TwoSLdown+=2){
296  if ((TwoSLdown>=0) && (abs(TwoSR-TwoSLdown)<=1)){
297  int TwoSRdown = TwoSLdown;
298 
299  int fase = (((TwoSL+1)%2)!=0)?-1:1;
300  const double factor = fase * sqrt((TwoJ+1)*(TwoSL+1.0)) * Wigner::wigner6j( 1, 1, TwoJ, TwoSL, TwoSR, TwoSRdown );
301 
302  for (int Irrep=0; Irrep<(denBK->getNumberOfIrreps()); Irrep++){
303 
304  int IrrepTimesMid = Irreps::directProd(Irrep,IprodMID);
305 
306  bool leftOK = false;
307  for (int l_alpha=0; l_alpha<theindex; l_alpha++){
308  if (denBK->gIrrep(l_alpha) == Irrep){ leftOK = true; }
309  }
310  bool rightOK = false;
311  for (int l_beta=theindex+2; l_beta<Prob->gL(); l_beta++){
312  if (denBK->gIrrep(l_beta) == IrrepTimesMid){ rightOK = true; }
313  }
314 
315  if ((leftOK) && (rightOK)){
316 
317  for (int l_alpha=0; l_alpha<theindex; l_alpha++){
318  if (denBK->gIrrep(l_alpha) == Irrep){
319 
320  int ILdown = Irreps::directProd(IL, Irrep);
321  int IRdown = Irreps::directProd(IR, IrrepTimesMid);
322 
323  int dimLdown = denBK->gCurrentDim(theindex, NL-1, TwoSLdown, ILdown);
324  int dimRdown = denBK->gCurrentDim(theindex+2, NR+1, TwoSRdown, IRdown);
325 
326  if ((dimLdown>0) && (dimRdown>0)){
327 
328  int size = dimRup * dimRdown;
329  for (int cnt=0; cnt<size; cnt++){ temp[cnt] = 0.0; }
330 
331  for (int l_beta=theindex+2; l_beta<Prob->gL(); l_beta++){
332  if (denBK->gIrrep(l_beta) == IrrepTimesMid){
333  double prefac = factor * ( Prob->gMxElement(l_alpha, l_beta, theindex, theindex+1)
334  + ((TwoJ==0)?1:-1) * Prob->gMxElement(l_alpha, l_beta, theindex+1, theindex) );
335  double * LblockR = Lright[l_beta-theindex-2]->gStorage(NR,TwoSR,IR,NR+1,TwoSRdown,IRdown);
336  daxpy_(&size,&prefac,LblockR,&inc,temp,&inc);
337  }
338  }
339 
340  double alpha = 1.0;
341  double beta = 0.0; //set
342  int memSkappa = denS->gKappa(NL-1, TwoSLdown, ILdown, 2, 2, 0, NR+1, TwoSRdown, IRdown);
343 
344  dgemm_(&notrans,&trans,&dimLdown,&dimRup,&dimRdown,&alpha,memS+denS->gKappa2index(memSkappa),&dimLdown,temp,&dimRup,&beta,temp2,&dimLdown);
345 
346  beta = 1.0; //add
347  double * LblockL = Lleft[theindex-1-l_alpha]->gStorage(NL-1,TwoSLdown,ILdown,NL,TwoSL,IL);
348  dgemm_(&trans,&notrans,&dimLup,&dimRup,&dimLdown,&alpha,LblockL,&dimLdown,temp2,&dimLdown,&beta,memHeff+denS->gKappa2index(ikappa),&dimLup);
349  }
350  }
351  }
352  }
353  }
354  }
355  }
356  }
357 
358 }
359 
360 void CheMPS2::Heff::addDiagram5B(const int ikappa, double * memS, double * memHeff, const Sobject * denS, TensorL ** Lleft, TensorL ** Lright, double * temp, double * temp2) const{
361 
362  #ifdef CHEMPS2_MPI_COMPILATION
363  const int MPIRANK = MPIchemps2::mpi_rank();
364  #endif
365 
366  int NL = denS->gNL(ikappa);
367  int TwoSL = denS->gTwoSL(ikappa);
368  int IL = denS->gIL(ikappa);
369 
370  int NR = denS->gNR(ikappa);
371  int TwoSR = denS->gTwoSR(ikappa);
372  int IR = denS->gIR(ikappa);
373 
374  int N1 = denS->gN1(ikappa);
375  int N2 = denS->gN2(ikappa);
376  int TwoJ = denS->gTwoJ(ikappa);
377 
378  int theindex = denS->gIndex();
379  int dimLup = denBK->gCurrentDim(theindex ,NL,TwoSL,IL);
380  int dimRup = denBK->gCurrentDim(theindex+2,NR,TwoSR,IR);
381 
382  char notrans = 'N';
383  int inc = 1;
384 
385  int IprodMID = Irreps::directProd(denBK->gIrrep(theindex),denBK->gIrrep(theindex+1));
386 
387  //5B1
388  #ifdef CHEMPS2_MPI_COMPILATION
389  if (( MPIchemps2::owner_specific_diagram( Prob->gL(), MPI_CHEMPS2_5B1 ) == MPIRANK ) && (N1==1) && (N2==1)){
390  #else
391  if ((N1==1) && (N2==1)){
392  #endif
393 
394  for (int TwoSLdown=TwoSL-1; TwoSLdown<=TwoSL+1; TwoSLdown+=2){
395  if ((TwoSLdown>=0) && (abs(TwoSR-TwoSLdown)<=1)){
396  int TwoSRdown = TwoSLdown;
397 
398  int fase = phase(TwoSL+TwoSR+2);
399  const double factor = fase * sqrt((TwoJ+1)*(TwoSR+1.0)) * Wigner::wigner6j( 1, 1, TwoJ, TwoSL, TwoSR, TwoSRdown );
400 
401  for (int Irrep=0; Irrep<(denBK->getNumberOfIrreps()); Irrep++){
402 
403  int IrrepTimesMid = Irreps::directProd(Irrep,IprodMID);
404 
405  bool leftOK = false;
406  for (int l_alpha=0; l_alpha<theindex; l_alpha++){
407  if (denBK->gIrrep(l_alpha) == Irrep){ leftOK = true; }
408  }
409  bool rightOK = false;
410  for (int l_beta=theindex+2; l_beta<Prob->gL(); l_beta++){
411  if (denBK->gIrrep(l_beta) == IrrepTimesMid){ rightOK = true; }
412  }
413 
414  if ((leftOK) && (rightOK)){
415 
416  for (int l_alpha=0; l_alpha<theindex; l_alpha++){
417  if (denBK->gIrrep(l_alpha) == Irrep){
418 
419  int ILdown = Irreps::directProd(IL, Irrep);
420  int IRdown = Irreps::directProd(IR, IrrepTimesMid);
421 
422  int dimLdown = denBK->gCurrentDim(theindex, NL+1, TwoSLdown, ILdown);
423  int dimRdown = denBK->gCurrentDim(theindex+2, NR-1, TwoSRdown, IRdown);
424 
425  if ((dimLdown>0) && (dimRdown>0)){
426 
427  int size = dimRup * dimRdown;
428  for (int cnt=0; cnt<size; cnt++){ temp[cnt] = 0.0; }
429 
430  for (int l_beta=theindex+2; l_beta<Prob->gL(); l_beta++){
431  if (denBK->gIrrep(l_beta) == IrrepTimesMid){
432  double prefac = factor * ( Prob->gMxElement(l_alpha, l_beta, theindex, theindex+1)
433  + ((TwoJ==0)?1:-1) * Prob->gMxElement(l_alpha, l_beta, theindex+1, theindex) );
434  double * LblockR = Lright[l_beta-theindex-2]->gStorage(NR-1,TwoSRdown,IRdown,NR,TwoSR,IR);
435  daxpy_(&size,&prefac,LblockR,&inc,temp,&inc);
436  }
437  }
438 
439  double alpha = 1.0;
440  double beta = 0.0; //set
441  int memSkappa = denS->gKappa(NL+1, TwoSLdown, ILdown, 0, 0, 0, NR-1, TwoSRdown, IRdown);
442 
443  dgemm_(&notrans,&notrans,&dimLdown,&dimRup,&dimRdown,&alpha,memS+denS->gKappa2index(memSkappa),&dimLdown,temp,&dimRdown,&beta,temp2,&dimLdown);
444 
445  beta = 1.0; //add
446  double * LblockL = Lleft[theindex-1-l_alpha]->gStorage(NL,TwoSL,IL,NL+1,TwoSLdown,ILdown);
447  dgemm_(&notrans,&notrans,&dimLup,&dimRup,&dimLdown,&alpha,LblockL,&dimLup,temp2,&dimLdown,&beta,memHeff+denS->gKappa2index(ikappa),&dimLup);
448  }
449  }
450  }
451  }
452  }
453  }
454  }
455  }
456 
457  //5B2
458  #ifdef CHEMPS2_MPI_COMPILATION
459  if (( MPIchemps2::owner_specific_diagram( Prob->gL(), MPI_CHEMPS2_5B2 ) == MPIRANK ) && (N1==2) && (N2==1)){
460  #else
461  if ((N1==2) && (N2==1)){
462  #endif
463 
464  for (int TwoSLdown=TwoSL-1; TwoSLdown<=TwoSL+1; TwoSLdown+=2){
465  for (int TwoSRdown=TwoSR-1; TwoSRdown<=TwoSR+1; TwoSRdown+=2){
466 
467  if (((TwoSL == TwoSRdown) || (TwoSLdown == TwoSR)) && (abs(TwoSLdown-TwoSRdown)<=1) && (TwoSLdown>=0) && (TwoSRdown>=0)){
468 
469  double factor1 = 0.0;
470  double factor2 = 0.0;
471 
472  if (TwoSLdown == TwoSR){ factor2 = phase(TwoSR + 1 - TwoSRdown); }
473  if (TwoSL == TwoSRdown){
474  int fase = phase(TwoSL + 1 - TwoSLdown);
475  factor1 = fase * sqrt(((TwoSR+1)*(TwoSLdown+1.0))/((TwoSRdown+1)*(TwoSL+1.0)));
476  }
477 
478  for (int Irrep=0; Irrep<(denBK->getNumberOfIrreps()); Irrep++){
479 
480  int IrrepTimesMid = Irreps::directProd(Irrep,IprodMID);
481 
482  bool leftOK = false;
483  for (int l_alpha=0; l_alpha<theindex; l_alpha++){
484  if (denBK->gIrrep(l_alpha) == Irrep){ leftOK = true; }
485  }
486  bool rightOK = false;
487  for (int l_beta=theindex+2; l_beta<Prob->gL(); l_beta++){
488  if (denBK->gIrrep(l_beta) == IrrepTimesMid){ rightOK = true; }
489  }
490 
491  if ((leftOK) && (rightOK)){
492 
493  for (int l_alpha=0; l_alpha<theindex; l_alpha++){
494  if (denBK->gIrrep(l_alpha) == Irrep){
495 
496  int ILdown = Irreps::directProd(IL, Irrep);
497  int IRdown = Irreps::directProd(IR, IrrepTimesMid);
498 
499  int dimLdown = denBK->gCurrentDim(theindex, NL+1, TwoSLdown, ILdown);
500  int dimRdown = denBK->gCurrentDim(theindex+2, NR-1, TwoSRdown, IRdown);
501 
502  if ((dimLdown>0) && (dimRdown>0)){
503 
504  int size = dimRup * dimRdown;
505  for (int cnt=0; cnt<size; cnt++){ temp[cnt] = 0.0; }
506 
507  for (int l_beta=theindex+2; l_beta<Prob->gL(); l_beta++){
508  if (denBK->gIrrep(l_beta) == IrrepTimesMid){
509  double prefac = factor1 * Prob->gMxElement(l_alpha, l_beta, theindex, theindex+1)
510  + factor2 * Prob->gMxElement(l_alpha, l_beta, theindex+1, theindex);
511  double * LblockR = Lright[l_beta-theindex-2]->gStorage(NR-1,TwoSRdown,IRdown,NR,TwoSR,IR);
512  daxpy_(&size,&prefac,LblockR,&inc,temp,&inc);
513  }
514  }
515 
516  double alpha = 1.0;
517  double beta = 0.0; //set
518  int memSkappa = denS->gKappa(NL+1, TwoSLdown, ILdown, 1, 0, 1, NR-1, TwoSRdown, IRdown);
519 
520  dgemm_(&notrans,&notrans,&dimLdown,&dimRup,&dimRdown,&alpha,memS+denS->gKappa2index(memSkappa),&dimLdown,temp,&dimRdown,&beta,temp2,&dimLdown);
521 
522  beta = 1.0; //add
523  double * LblockL = Lleft[theindex-1-l_alpha]->gStorage(NL,TwoSL,IL,NL+1,TwoSLdown,ILdown);
524  dgemm_(&notrans,&notrans,&dimLup,&dimRup,&dimLdown,&alpha,LblockL,&dimLup,temp2,&dimLdown,&beta,memHeff+denS->gKappa2index(ikappa),&dimLup);
525 
526  }
527  }
528  }
529  }
530  }
531  }
532  }
533  }
534  }
535 
536  //5B3
537  #ifdef CHEMPS2_MPI_COMPILATION
538  if (( MPIchemps2::owner_specific_diagram( Prob->gL(), MPI_CHEMPS2_5B3 ) == MPIRANK ) && (N1==1) && (N2==2)){
539  #else
540  if ((N1==1) && (N2==2)){
541  #endif
542 
543  for (int TwoSLdown=TwoSL-1; TwoSLdown<=TwoSL+1; TwoSLdown+=2){
544  for (int TwoSRdown=TwoSR-1; TwoSRdown<=TwoSR+1; TwoSRdown+=2){
545 
546  if (((TwoSL == TwoSRdown) || (TwoSLdown == TwoSR)) && (abs(TwoSLdown-TwoSRdown)<=1) && (TwoSLdown>=0) && (TwoSRdown>=0)){
547 
548  double factor1 = 0.0;
549  double factor2 = 0.0;
550 
551  if (TwoSLdown == TwoSR){ factor1 = phase(TwoSR + 1 - TwoSRdown); }
552  if (TwoSL == TwoSRdown){
553  factor2 = phase(TwoSL + 1 - TwoSLdown) * sqrt(((TwoSR+1)*(TwoSLdown+1.0))/((TwoSRdown+1)*(TwoSL+1.0)));
554  }
555 
556  for (int Irrep=0; Irrep<(denBK->getNumberOfIrreps()); Irrep++){
557 
558  int IrrepTimesMid = Irreps::directProd(Irrep,IprodMID);
559 
560  bool leftOK = false;
561  for (int l_alpha=0; l_alpha<theindex; l_alpha++){
562  if (denBK->gIrrep(l_alpha) == Irrep){ leftOK = true; }
563  }
564  bool rightOK = false;
565  for (int l_beta=theindex+2; l_beta<Prob->gL(); l_beta++){
566  if (denBK->gIrrep(l_beta) == IrrepTimesMid){ rightOK = true; }
567  }
568 
569  if ((leftOK) && (rightOK)){
570 
571  for (int l_alpha=0; l_alpha<theindex; l_alpha++){
572  if (denBK->gIrrep(l_alpha) == Irrep){
573 
574  int ILdown = Irreps::directProd(IL, Irrep);
575  int IRdown = Irreps::directProd(IR, IrrepTimesMid);
576 
577  int dimLdown = denBK->gCurrentDim(theindex, NL+1, TwoSLdown, ILdown);
578  int dimRdown = denBK->gCurrentDim(theindex+2, NR-1, TwoSRdown, IRdown);
579 
580  if ((dimLdown>0) && (dimRdown>0)){
581 
582  int size = dimRup * dimRdown;
583  for (int cnt=0; cnt<size; cnt++){ temp[cnt] = 0.0; }
584 
585  for (int l_beta=theindex+2; l_beta<Prob->gL(); l_beta++){
586  if (denBK->gIrrep(l_beta) == IrrepTimesMid){
587  double prefac = factor1 * Prob->gMxElement(l_alpha, l_beta, theindex, theindex+1)
588  + factor2 * Prob->gMxElement(l_alpha, l_beta, theindex+1, theindex);
589  double * LblockR = Lright[l_beta-theindex-2]->gStorage(NR-1,TwoSRdown,IRdown,NR,TwoSR,IR);
590  daxpy_(&size,&prefac,LblockR,&inc,temp,&inc);
591  }
592  }
593 
594  double alpha = 1.0;
595  double beta = 0.0; //set
596  int memSkappa = denS->gKappa(NL+1, TwoSLdown, ILdown, 0, 1, 1, NR-1, TwoSRdown, IRdown);
597 
598  dgemm_(&notrans,&notrans,&dimLdown,&dimRup,&dimRdown,&alpha,memS+denS->gKappa2index(memSkappa),&dimLdown,temp,&dimRdown,&beta,temp2,&dimLdown);
599 
600  beta = 1.0; //add
601  double * LblockL = Lleft[theindex-1-l_alpha]->gStorage(NL,TwoSL,IL,NL+1,TwoSLdown,ILdown);
602  dgemm_(&notrans,&notrans,&dimLup,&dimRup,&dimLdown,&alpha,LblockL,&dimLup,temp2,&dimLdown,&beta,memHeff+denS->gKappa2index(ikappa),&dimLup);
603 
604  }
605  }
606  }
607  }
608  }
609  }
610  }
611  }
612  }
613 
614  //5B4
615  #ifdef CHEMPS2_MPI_COMPILATION
616  if (( MPIchemps2::owner_specific_diagram( Prob->gL(), MPI_CHEMPS2_5B4 ) == MPIRANK ) && (N1==2) && (N2==2)){
617  #else
618  if ((N1==2) && (N2==2)){
619  #endif
620 
621  for (int TwoSLdown=TwoSL-1; TwoSLdown<=TwoSL+1; TwoSLdown+=2){
622  for (int TwoSRdown=TwoSR-1; TwoSRdown<=TwoSR+1; TwoSRdown+=2){
623  for (int TwoJdown=0; TwoJdown<=2; TwoJdown+=2){
624  if ((abs(TwoSLdown-TwoSRdown)<=TwoJdown) && (TwoSLdown>=0) && (TwoSRdown>=0)){
625 
626  int fase = (((TwoSLdown+1)%2)!=0)?-1:1;
627  const double factor = fase * sqrt((TwoJdown+1)*(TwoSLdown+1.0)) * Wigner::wigner6j( 1, 1, TwoJdown, TwoSLdown, TwoSRdown, TwoSR );
628 
629  for (int Irrep=0; Irrep<(denBK->getNumberOfIrreps()); Irrep++){
630 
631  int IrrepTimesMid = Irreps::directProd(Irrep,IprodMID);
632 
633  bool leftOK = false;
634  for (int l_alpha=0; l_alpha<theindex; l_alpha++){
635  if (denBK->gIrrep(l_alpha) == Irrep){ leftOK = true; }
636  }
637  bool rightOK = false;
638  for (int l_beta=theindex+2; l_beta<Prob->gL(); l_beta++){
639  if (denBK->gIrrep(l_beta) == IrrepTimesMid){ rightOK = true; }
640  }
641 
642  if ((leftOK) && (rightOK)){
643 
644  for (int l_alpha=0; l_alpha<theindex; l_alpha++){
645  if (denBK->gIrrep(l_alpha) == Irrep){
646 
647  int ILdown = Irreps::directProd(IL, Irrep);
648  int IRdown = Irreps::directProd(IR, IrrepTimesMid);
649 
650  int dimLdown = denBK->gCurrentDim(theindex, NL+1, TwoSLdown, ILdown);
651  int dimRdown = denBK->gCurrentDim(theindex+2, NR-1, TwoSRdown, IRdown);
652 
653  if ((dimLdown>0) && (dimRdown>0)){
654 
655  int size = dimRup * dimRdown;
656  for (int cnt=0; cnt<size; cnt++){ temp[cnt] = 0.0; }
657 
658  for (int l_beta=theindex+2; l_beta<Prob->gL(); l_beta++){
659  if (denBK->gIrrep(l_beta) == IrrepTimesMid){
660  double prefac = factor * ( Prob->gMxElement(l_alpha, l_beta, theindex, theindex+1)
661  + ((TwoJdown==0)?1:-1) * Prob->gMxElement(l_alpha, l_beta, theindex+1, theindex) );
662  double * LblockR = Lright[l_beta-theindex-2]->gStorage(NR-1,TwoSRdown,IRdown,NR,TwoSR,IR);
663  daxpy_(&size,&prefac,LblockR,&inc,temp,&inc);
664  }
665  }
666 
667  double alpha = 1.0;
668  double beta = 0.0; //set
669  int memSkappa = denS->gKappa(NL+1, TwoSLdown, ILdown, 1, 1, TwoJdown, NR-1, TwoSRdown, IRdown);
670 
671  dgemm_(&notrans,&notrans,&dimLdown,&dimRup,&dimRdown,&alpha,memS+denS->gKappa2index(memSkappa),&dimLdown,temp,&dimRdown,&beta,temp2,&dimLdown);
672 
673  beta = 1.0; //add
674  double * LblockL = Lleft[theindex-1-l_alpha]->gStorage(NL,TwoSL,IL,NL+1,TwoSLdown,ILdown);
675  dgemm_(&notrans,&notrans,&dimLup,&dimRup,&dimLdown,&alpha,LblockL,&dimLup,temp2,&dimLdown,&beta,memHeff+denS->gKappa2index(ikappa),&dimLup);
676  }
677  }
678  }
679  }
680  }
681  }
682  }
683  }
684  }
685  }
686 
687 }
688 
689 void CheMPS2::Heff::addDiagram5C(const int ikappa, double * memS, double * memHeff, const Sobject * denS, TensorL ** Lleft, TensorL ** Lright, double * temp, double * temp2) const{
690 
691  #ifdef CHEMPS2_MPI_COMPILATION
692  const int MPIRANK = MPIchemps2::mpi_rank();
693  #endif
694 
695  int NL = denS->gNL(ikappa);
696  int TwoSL = denS->gTwoSL(ikappa);
697  int IL = denS->gIL(ikappa);
698 
699  int NR = denS->gNR(ikappa);
700  int TwoSR = denS->gTwoSR(ikappa);
701  int IR = denS->gIR(ikappa);
702 
703  int N1 = denS->gN1(ikappa);
704  int N2 = denS->gN2(ikappa);
705  int TwoJ = denS->gTwoJ(ikappa);
706 
707  int theindex = denS->gIndex();
708  int dimLup = denBK->gCurrentDim(theindex ,NL,TwoSL,IL);
709  int dimRup = denBK->gCurrentDim(theindex+2,NR,TwoSR,IR);
710 
711  char notrans = 'N';
712  char trans = 'T';
713  int inc = 1;
714 
715  int IprodMID = Irreps::directProd(denBK->gIrrep(theindex),denBK->gIrrep(theindex+1));
716 
717  //5C1
718  #ifdef CHEMPS2_MPI_COMPILATION
719  if (( MPIchemps2::owner_specific_diagram( Prob->gL(), MPI_CHEMPS2_5C1 ) == MPIRANK ) && (N1==1) && (N2==0)){
720  #else
721  if ((N1==1) && (N2==0)){
722  #endif
723 
724  for (int TwoSLdown=TwoSL-1; TwoSLdown<=TwoSL+1; TwoSLdown+=2){
725  for (int TwoSRdown=TwoSR-1; TwoSRdown<=TwoSR+1; TwoSRdown+=2){
726  if ((abs(TwoSLdown-TwoSRdown)<=1) && (TwoSLdown>=0) && (TwoSRdown>=0)){
727 
728  int fase = phase(TwoSL+TwoSRdown);
729  const double factor2 = fase * sqrt((TwoSL+1)*(TwoSR+1.0)) * Wigner::wigner6j( TwoSLdown, TwoSRdown, 1, TwoSR, TwoSL, 1 );
730  const double factor1 = (TwoSL==TwoSRdown)?sqrt((TwoSR+1.0)/(TwoSRdown+1.0)):0.0;
731 
732  for (int Irrep=0; Irrep<(denBK->getNumberOfIrreps()); Irrep++){
733 
734  int IrrepTimesMid = Irreps::directProd(Irrep,IprodMID);
735 
736  bool leftOK = false;
737  for (int l_alpha=0; l_alpha<theindex; l_alpha++){
738  if (denBK->gIrrep(l_alpha) == Irrep){ leftOK = true; }
739  }
740  bool rightOK = false;
741  for (int l_beta=theindex+2; l_beta<Prob->gL(); l_beta++){
742  if (denBK->gIrrep(l_beta) == IrrepTimesMid){ rightOK = true; }
743  }
744 
745  if ((leftOK) && (rightOK)){
746 
747  for (int l_alpha=0; l_alpha<theindex; l_alpha++){
748  if (denBK->gIrrep(l_alpha) == Irrep){
749 
750  int ILdown = Irreps::directProd(IL, Irrep);
751  int IRdown = Irreps::directProd(IR, IrrepTimesMid);
752 
753  int dimLdown = denBK->gCurrentDim(theindex, NL-1, TwoSLdown, ILdown);
754  int dimRdown = denBK->gCurrentDim(theindex+2, NR-1, TwoSRdown, IRdown);
755 
756  if ((dimLdown>0) && (dimRdown>0)){
757 
758  int size = dimRup * dimRdown;
759  for (int cnt=0; cnt<size; cnt++){ temp[cnt] = 0.0; }
760 
761  for (int l_beta=theindex+2; l_beta<Prob->gL(); l_beta++){
762  if (denBK->gIrrep(l_beta) == IrrepTimesMid){
763  double prefac = factor1 * Prob->gMxElement(l_alpha,theindex,theindex+1,l_beta)
764  + factor2 * Prob->gMxElement(l_alpha,theindex,l_beta,theindex+1);
765  double * LblockR = Lright[l_beta-theindex-2]->gStorage(NR-1,TwoSRdown,IRdown,NR,TwoSR,IR);
766  daxpy_(&size,&prefac,LblockR,&inc,temp,&inc);
767  }
768  }
769 
770  double alpha = 1.0;
771  double beta = 0.0; //set
772  int memSkappa = denS->gKappa(NL-1, TwoSLdown, ILdown, 0, 1, 1, NR-1, TwoSRdown, IRdown);
773 
774  dgemm_(&notrans,&notrans,&dimLdown,&dimRup,&dimRdown,&alpha,memS+denS->gKappa2index(memSkappa),&dimLdown,temp,&dimRdown,&beta,temp2,&dimLdown);
775 
776  beta = 1.0; //add
777  double * LblockL = Lleft[theindex-1-l_alpha]->gStorage(NL-1,TwoSLdown,ILdown,NL,TwoSL,IL);
778  dgemm_(&trans,&notrans,&dimLup,&dimRup,&dimLdown,&alpha,LblockL,&dimLdown,temp2,&dimLdown,&beta,memHeff+denS->gKappa2index(ikappa),&dimLup);
779  }
780  }
781  }
782  }
783  }
784  }
785  }
786  }
787  }
788 
789  //5C2
790  #ifdef CHEMPS2_MPI_COMPILATION
791  if (( MPIchemps2::owner_specific_diagram( Prob->gL(), MPI_CHEMPS2_5C2 ) == MPIRANK ) && (N1==2) && (N2==0)){
792  #else
793  if ((N1==2) && (N2==0)){
794  #endif
795 
796  for (int TwoSLdown=TwoSL-1; TwoSLdown<=TwoSL+1; TwoSLdown+=2){
797  for (int TwoSRdown=TwoSR-1; TwoSRdown<=TwoSR+1; TwoSRdown+=2){
798  for (int TwoJdown=0; TwoJdown<=2; TwoJdown+=2){
799  if ((abs(TwoSLdown-TwoSRdown)<=TwoJdown) && (TwoSLdown>=0) && (TwoSRdown>=0)){
800 
801  int fase = phase(TwoSR + TwoSLdown + 3 + TwoJdown);
802  const double factor1 = fase * sqrt((TwoSR+1)*(TwoJdown+1.0)) * Wigner::wigner6j( 1, 1, TwoJdown, TwoSLdown, TwoSRdown, TwoSR );
803  const double factor2 = (TwoJdown==0)?sqrt(2.0*(TwoSR+1.0)/(TwoSRdown+1.0)):0.0;
804 
805  for (int Irrep=0; Irrep<(denBK->getNumberOfIrreps()); Irrep++){
806 
807  int IrrepTimesMid = Irreps::directProd(Irrep,IprodMID);
808 
809  bool leftOK = false;
810  for (int l_alpha=0; l_alpha<theindex; l_alpha++){
811  if (denBK->gIrrep(l_alpha) == Irrep){ leftOK = true; }
812  }
813  bool rightOK = false;
814  for (int l_beta=theindex+2; l_beta<Prob->gL(); l_beta++){
815  if (denBK->gIrrep(l_beta) == IrrepTimesMid){ rightOK = true; }
816  }
817 
818  if ((leftOK) && (rightOK)){
819 
820  for (int l_alpha=0; l_alpha<theindex; l_alpha++){
821  if (denBK->gIrrep(l_alpha) == Irrep){
822 
823  int ILdown = Irreps::directProd(IL, Irrep);
824  int IRdown = Irreps::directProd(IR, IrrepTimesMid);
825 
826  int dimLdown = denBK->gCurrentDim(theindex, NL-1, TwoSLdown, ILdown);
827  int dimRdown = denBK->gCurrentDim(theindex+2, NR-1, TwoSRdown, IRdown);
828 
829  if ((dimLdown>0) && (dimRdown>0)){
830 
831  int size = dimRup * dimRdown;
832  for (int cnt=0; cnt<size; cnt++){ temp[cnt] = 0.0; }
833 
834  for (int l_beta=theindex+2; l_beta<Prob->gL(); l_beta++){
835  if (denBK->gIrrep(l_beta) == IrrepTimesMid){
836  double prefac = factor1 * Prob->gMxElement(l_alpha, theindex, theindex+1, l_beta)
837  + factor2 * Prob->gMxElement(l_alpha, theindex, l_beta, theindex+1);
838  double * LblockR = Lright[l_beta-theindex-2]->gStorage(NR-1,TwoSRdown,IRdown,NR,TwoSR,IR);
839  daxpy_(&size,&prefac,LblockR,&inc,temp,&inc);
840  }
841  }
842 
843  double alpha = 1.0;
844  double beta = 0.0; //set
845  int memSkappa = denS->gKappa(NL-1, TwoSLdown, ILdown, 1, 1, TwoJdown, NR-1, TwoSRdown, IRdown);
846 
847  dgemm_(&notrans,&notrans,&dimLdown,&dimRup,&dimRdown,&alpha,memS+denS->gKappa2index(memSkappa),&dimLdown,temp,&dimRdown,&beta,temp2,&dimLdown);
848 
849  beta = 1.0; //add
850  double * LblockL = Lleft[theindex-1-l_alpha]->gStorage(NL-1,TwoSLdown,ILdown,NL,TwoSL,IL);
851  dgemm_(&trans,&notrans,&dimLup,&dimRup,&dimLdown,&alpha,LblockL,&dimLdown,temp2,&dimLdown,&beta,memHeff+denS->gKappa2index(ikappa),&dimLup);
852  }
853  }
854  }
855  }
856  }
857  }
858  }
859  }
860  }
861  }
862 
863  //5C3
864  #ifdef CHEMPS2_MPI_COMPILATION
865  if (( MPIchemps2::owner_specific_diagram( Prob->gL(), MPI_CHEMPS2_5C3 ) == MPIRANK ) && (N1==1) && (N2==1)){
866  #else
867  if ((N1==1) && (N2==1)){
868  #endif
869 
870  for (int TwoSLdown=TwoSL-1; TwoSLdown<=TwoSL+1; TwoSLdown+=2){
871  if ((TwoSLdown>=0) && (abs(TwoSR-TwoSLdown)<=1)){
872  int TwoSRdown = TwoSLdown;
873 
874  int fase = phase(TwoSR + TwoSRdown + 3 + TwoJ);
875  double factor1 = fase * sqrt((TwoSL+1.0)*(TwoJ+1.0)*(TwoSR+1.0)/(TwoSRdown+1.0)) * Wigner::wigner6j( 1, 1, TwoJ, TwoSL, TwoSR, TwoSRdown );
876  double factor2 = (TwoJ==0)?sqrt(2.0*(TwoSR+1.0)/(TwoSRdown+1.0)):0.0;
877 
878  for (int Irrep=0; Irrep<(denBK->getNumberOfIrreps()); Irrep++){
879 
880  int IrrepTimesMid = Irreps::directProd(Irrep,IprodMID);
881 
882  bool leftOK = false;
883  for (int l_alpha=0; l_alpha<theindex; l_alpha++){
884  if (denBK->gIrrep(l_alpha) == Irrep){ leftOK = true; }
885  }
886  bool rightOK = false;
887  for (int l_beta=theindex+2; l_beta<Prob->gL(); l_beta++){
888  if (denBK->gIrrep(l_beta) == IrrepTimesMid){ rightOK = true; }
889  }
890 
891  if ((leftOK) && (rightOK)){
892 
893  for (int l_alpha=0; l_alpha<theindex; l_alpha++){
894  if (denBK->gIrrep(l_alpha) == Irrep){
895 
896  int ILdown = Irreps::directProd(IL, Irrep);
897  int IRdown = Irreps::directProd(IR, IrrepTimesMid);
898 
899  int dimLdown = denBK->gCurrentDim(theindex, NL-1, TwoSLdown, ILdown);
900  int dimRdown = denBK->gCurrentDim(theindex+2, NR-1, TwoSRdown, IRdown);
901 
902  if ((dimLdown>0) && (dimRdown>0)){
903 
904  int size = dimRup * dimRdown;
905  for (int cnt=0; cnt<size; cnt++){ temp[cnt] = 0.0; }
906 
907  for (int l_beta=theindex+2; l_beta<Prob->gL(); l_beta++){
908  if (denBK->gIrrep(l_beta) == IrrepTimesMid){
909  double prefac = factor1 * Prob->gMxElement(l_alpha, theindex, theindex+1, l_beta)
910  + factor2 * Prob->gMxElement(l_alpha, theindex, l_beta, theindex+1);
911  double * LblockR = Lright[l_beta-theindex-2]->gStorage(NR-1,TwoSRdown,IRdown,NR,TwoSR,IR);
912  daxpy_(&size,&prefac,LblockR,&inc,temp,&inc);
913  }
914  }
915 
916  double alpha = 1.0;
917  double beta = 0.0; //set
918  int memSkappa = denS->gKappa(NL-1, TwoSLdown, ILdown, 0, 2, 0, NR-1, TwoSRdown, IRdown);
919 
920  dgemm_(&notrans,&notrans,&dimLdown,&dimRup,&dimRdown,&alpha,memS+denS->gKappa2index(memSkappa),&dimLdown,temp,&dimRdown,&beta,temp2,&dimLdown);
921 
922  beta = 1.0; //add
923  double * LblockL = Lleft[theindex-1-l_alpha]->gStorage(NL-1,TwoSLdown,ILdown,NL,TwoSL,IL);
924  dgemm_(&trans,&notrans,&dimLup,&dimRup,&dimLdown,&alpha,LblockL,&dimLdown,temp2,&dimLdown,&beta,memHeff+denS->gKappa2index(ikappa),&dimLup);
925  }
926  }
927  }
928  }
929  }
930  }
931  }
932  }
933 
934  //5C4
935  #ifdef CHEMPS2_MPI_COMPILATION
936  if (( MPIchemps2::owner_specific_diagram( Prob->gL(), MPI_CHEMPS2_5C4 ) == MPIRANK ) && (N1==2) && (N2==1)){
937  #else
938  if ((N1==2) && (N2==1)){
939  #endif
940 
941  for (int TwoSLdown=TwoSL-1; TwoSLdown<=TwoSL+1; TwoSLdown+=2){
942  for (int TwoSRdown=TwoSR-1; TwoSRdown<=TwoSR+1; TwoSRdown+=2){
943  if ((abs(TwoSLdown-TwoSRdown)<=1) && (TwoSLdown>=0) && (TwoSRdown>=0)){
944 
945  const double factor1 = (TwoSLdown==TwoSR) ? phase(TwoSL-TwoSRdown) * sqrt((TwoSL+1.0)/(TwoSLdown+1.0)) : 0.0;
946  const double factor2 = phase(TwoSL+TwoSRdown+2) * sqrt((TwoSL+1)*(TwoSR+1.0)) * Wigner::wigner6j( TwoSLdown, TwoSRdown, 1, TwoSR, TwoSL, 1 );
947 
948  for (int Irrep=0; Irrep<(denBK->getNumberOfIrreps()); Irrep++){
949 
950  int IrrepTimesMid = Irreps::directProd(Irrep,IprodMID);
951 
952  bool leftOK = false;
953  for (int l_alpha=0; l_alpha<theindex; l_alpha++){
954  if (denBK->gIrrep(l_alpha) == Irrep){ leftOK = true; }
955  }
956  bool rightOK = false;
957  for (int l_beta=theindex+2; l_beta<Prob->gL(); l_beta++){
958  if (denBK->gIrrep(l_beta) == IrrepTimesMid){ rightOK = true; }
959  }
960 
961  if ((leftOK) && (rightOK)){
962 
963  for (int l_alpha=0; l_alpha<theindex; l_alpha++){
964  if (denBK->gIrrep(l_alpha) == Irrep){
965 
966  int ILdown = Irreps::directProd(IL, Irrep);
967  int IRdown = Irreps::directProd(IR, IrrepTimesMid);
968 
969  int dimLdown = denBK->gCurrentDim(theindex, NL-1, TwoSLdown, ILdown);
970  int dimRdown = denBK->gCurrentDim(theindex+2, NR-1, TwoSRdown, IRdown);
971 
972  if ((dimLdown>0) && (dimRdown>0)){
973 
974  int size = dimRup * dimRdown;
975  for (int cnt=0; cnt<size; cnt++){ temp[cnt] = 0.0; }
976 
977  for (int l_beta=theindex+2; l_beta<Prob->gL(); l_beta++){
978  if (denBK->gIrrep(l_beta) == IrrepTimesMid){
979  double prefac = factor1 * Prob->gMxElement(l_alpha, theindex, theindex+1, l_beta)
980  + factor2 * Prob->gMxElement(l_alpha, theindex, l_beta, theindex+1);
981  double * LblockR = Lright[l_beta-theindex-2]->gStorage(NR-1,TwoSRdown,IRdown,NR,TwoSR,IR);
982  daxpy_(&size,&prefac,LblockR,&inc,temp,&inc);
983  }
984  }
985 
986  double alpha = 1.0;
987  double beta = 0.0; //set
988  int memSkappa = denS->gKappa(NL-1, TwoSLdown, ILdown, 1, 2, 1, NR-1, TwoSRdown, IRdown);
989 
990  dgemm_(&notrans,&notrans,&dimLdown,&dimRup,&dimRdown,&alpha,memS+denS->gKappa2index(memSkappa),&dimLdown,temp,&dimRdown,&beta,temp2,&dimLdown);
991 
992  beta = 1.0; //add
993  double * LblockL = Lleft[theindex-1-l_alpha]->gStorage(NL-1,TwoSLdown,ILdown,NL,TwoSL,IL);
994  dgemm_(&trans,&notrans,&dimLup,&dimRup,&dimLdown,&alpha,LblockL,&dimLdown,temp2,&dimLdown,&beta,memHeff+denS->gKappa2index(ikappa),&dimLup);
995  }
996  }
997  }
998  }
999  }
1000  }
1001  }
1002  }
1003  }
1004 
1005 }
1006 
1007 void CheMPS2::Heff::addDiagram5D(const int ikappa, double * memS, double * memHeff, const Sobject * denS, TensorL ** Lleft, TensorL ** Lright, double * temp, double * temp2) const{
1008 
1009  #ifdef CHEMPS2_MPI_COMPILATION
1010  const int MPIRANK = MPIchemps2::mpi_rank();
1011  #endif
1012 
1013  int NL = denS->gNL(ikappa);
1014  int TwoSL = denS->gTwoSL(ikappa);
1015  int IL = denS->gIL(ikappa);
1016 
1017  int NR = denS->gNR(ikappa);
1018  int TwoSR = denS->gTwoSR(ikappa);
1019  int IR = denS->gIR(ikappa);
1020 
1021  int N1 = denS->gN1(ikappa);
1022  int N2 = denS->gN2(ikappa);
1023  int TwoJ = denS->gTwoJ(ikappa);
1024 
1025  int theindex = denS->gIndex();
1026  int dimLup = denBK->gCurrentDim(theindex ,NL,TwoSL,IL);
1027  int dimRup = denBK->gCurrentDim(theindex+2,NR,TwoSR,IR);
1028 
1029  char notrans = 'N';
1030  char trans = 'T';
1031  int inc = 1;
1032 
1033  int IprodMID = Irreps::directProd(denBK->gIrrep(theindex),denBK->gIrrep(theindex+1));
1034 
1035  //5D1
1036  #ifdef CHEMPS2_MPI_COMPILATION
1037  if (( MPIchemps2::owner_specific_diagram( Prob->gL(), MPI_CHEMPS2_5D1 ) == MPIRANK ) && (N1==0) && (N2==1)){
1038  #else
1039  if ((N1==0) && (N2==1)){
1040  #endif
1041 
1042  for (int TwoSLdown=TwoSL-1; TwoSLdown<=TwoSL+1; TwoSLdown+=2){
1043  for (int TwoSRdown=TwoSR-1; TwoSRdown<=TwoSR+1; TwoSRdown+=2){
1044  if ((abs(TwoSLdown-TwoSRdown)<=1) && (TwoSLdown>=0) && (TwoSRdown>=0)){
1045 
1046  int fase = phase(TwoSLdown+TwoSR);
1047  const double factor2 = fase * sqrt((TwoSLdown+1)*(TwoSRdown+1.0)) * Wigner::wigner6j( TwoSL, TwoSR, 1, TwoSRdown, TwoSLdown, 1 );
1048  const double factor1 = (TwoSLdown==TwoSR)?sqrt((TwoSRdown+1.0)/(TwoSR+1.0)):0.0;
1049 
1050  for (int Irrep=0; Irrep<(denBK->getNumberOfIrreps()); Irrep++){
1051 
1052  int IrrepTimesMid = Irreps::directProd(Irrep,IprodMID);
1053 
1054  bool leftOK = false;
1055  for (int l_alpha=0; l_alpha<theindex; l_alpha++){
1056  if (denBK->gIrrep(l_alpha) == Irrep){ leftOK = true; }
1057  }
1058  bool rightOK = false;
1059  for (int l_beta=theindex+2; l_beta<Prob->gL(); l_beta++){
1060  if (denBK->gIrrep(l_beta) == IrrepTimesMid){ rightOK = true; }
1061  }
1062 
1063  if ((leftOK) && (rightOK)){
1064 
1065  for (int l_alpha=0; l_alpha<theindex; l_alpha++){
1066  if (denBK->gIrrep(l_alpha) == Irrep){
1067 
1068  int ILdown = Irreps::directProd(IL, Irrep);
1069  int IRdown = Irreps::directProd(IR, IrrepTimesMid);
1070 
1071  int dimLdown = denBK->gCurrentDim(theindex, NL+1, TwoSLdown, ILdown);
1072  int dimRdown = denBK->gCurrentDim(theindex+2, NR+1, TwoSRdown, IRdown);
1073 
1074  if ((dimLdown>0) && (dimRdown>0)){
1075 
1076  int size = dimRup * dimRdown;
1077  for (int cnt=0; cnt<size; cnt++){ temp[cnt] = 0.0; }
1078 
1079  for (int l_beta=theindex+2; l_beta<Prob->gL(); l_beta++){
1080  if (denBK->gIrrep(l_beta) == IrrepTimesMid){
1081  double prefac = factor1 * Prob->gMxElement(l_alpha,theindex,theindex+1,l_beta)
1082  + factor2 * Prob->gMxElement(l_alpha,theindex,l_beta,theindex+1);
1083  double * LblockR = Lright[l_beta-theindex-2]->gStorage(NR,TwoSR,IR,NR+1,TwoSRdown,IRdown);
1084  daxpy_(&size,&prefac,LblockR,&inc,temp,&inc);
1085  }
1086  }
1087 
1088  double alpha = 1.0;
1089  double beta = 0.0; //set
1090  int memSkappa = denS->gKappa(NL+1, TwoSLdown, ILdown, 1, 0, 1, NR+1, TwoSRdown, IRdown);
1091 
1092  dgemm_(&notrans,&trans,&dimLdown,&dimRup,&dimRdown,&alpha,memS+denS->gKappa2index(memSkappa),&dimLdown,temp,&dimRup,&beta,temp2,&dimLdown);
1093 
1094  beta = 1.0; //add
1095  double * LblockL = Lleft[theindex-1-l_alpha]->gStorage(NL,TwoSL,IL,NL+1,TwoSLdown,ILdown);
1096  dgemm_(&notrans,&notrans,&dimLup,&dimRup,&dimLdown,&alpha,LblockL,&dimLup,temp2,&dimLdown,&beta,memHeff+denS->gKappa2index(ikappa),&dimLup);
1097 
1098  }
1099  }
1100  }
1101  }
1102  }
1103  }
1104  }
1105  }
1106  }
1107 
1108  //5D2
1109  #ifdef CHEMPS2_MPI_COMPILATION
1110  if (( MPIchemps2::owner_specific_diagram( Prob->gL(), MPI_CHEMPS2_5D2 ) == MPIRANK ) && (N1==1) && (N2==1)){
1111  #else
1112  if ((N1==1) && (N2==1)){
1113  #endif
1114 
1115  for (int TwoSLdown=TwoSL-1; TwoSLdown<=TwoSL+1; TwoSLdown+=2){
1116  if ((TwoSLdown>=0) && (abs(TwoSR-TwoSLdown)<=1)){
1117  int TwoSRdown = TwoSLdown;
1118 
1119  int fase = phase(TwoSRdown + TwoSL + 3 + TwoJ);
1120  const double factor1 = fase * sqrt((TwoSRdown+1)*(TwoJ+1.0)) * Wigner::wigner6j( 1, 1, TwoJ, TwoSL, TwoSR, TwoSRdown );
1121  const double factor2 = (TwoJ==0)?sqrt(2.0*(TwoSRdown+1.0)/(TwoSR+1.0)):0.0;
1122 
1123  for (int Irrep=0; Irrep<(denBK->getNumberOfIrreps()); Irrep++){
1124 
1125  int IrrepTimesMid = Irreps::directProd(Irrep,IprodMID);
1126 
1127  bool leftOK = false;
1128  for (int l_alpha=0; l_alpha<theindex; l_alpha++){
1129  if (denBK->gIrrep(l_alpha) == Irrep){ leftOK = true; }
1130  }
1131  bool rightOK = false;
1132  for (int l_beta=theindex+2; l_beta<Prob->gL(); l_beta++){
1133  if (denBK->gIrrep(l_beta) == IrrepTimesMid){ rightOK = true; }
1134  }
1135 
1136  if ((leftOK) && (rightOK)){
1137 
1138  for (int l_alpha=0; l_alpha<theindex; l_alpha++){
1139  if (denBK->gIrrep(l_alpha) == Irrep){
1140 
1141  int ILdown = Irreps::directProd(IL, Irrep);
1142  int IRdown = Irreps::directProd(IR, IrrepTimesMid);
1143 
1144  int dimLdown = denBK->gCurrentDim(theindex, NL+1, TwoSLdown, ILdown);
1145  int dimRdown = denBK->gCurrentDim(theindex+2, NR+1, TwoSRdown, IRdown);
1146 
1147  if ((dimLdown>0) && (dimRdown>0)){
1148 
1149  int size = dimRup * dimRdown;
1150  for (int cnt=0; cnt<size; cnt++){ temp[cnt] = 0.0; }
1151 
1152  for (int l_beta=theindex+2; l_beta<Prob->gL(); l_beta++){
1153  if (denBK->gIrrep(l_beta) == IrrepTimesMid){
1154  double prefac = factor1 * Prob->gMxElement(l_alpha, theindex, theindex+1, l_beta)
1155  + factor2 * Prob->gMxElement(l_alpha, theindex, l_beta, theindex+1);
1156  double * LblockR = Lright[l_beta-theindex-2]->gStorage(NR,TwoSR,IR,NR+1,TwoSRdown,IRdown);
1157  daxpy_(&size,&prefac,LblockR,&inc,temp,&inc);
1158  }
1159  }
1160 
1161  double alpha = 1.0;
1162  double beta = 0.0; //set
1163  int memSkappa = denS->gKappa(NL+1, TwoSLdown, ILdown, 2, 0, 0, NR+1, TwoSRdown, IRdown);
1164 
1165  dgemm_(&notrans,&trans,&dimLdown,&dimRup,&dimRdown,&alpha,memS+denS->gKappa2index(memSkappa),&dimLdown,temp,&dimRup,&beta,temp2,&dimLdown);
1166 
1167  beta = 1.0; //add
1168  double * LblockL = Lleft[theindex-1-l_alpha]->gStorage(NL,TwoSL,IL,NL+1,TwoSLdown,ILdown);
1169  dgemm_(&notrans,&notrans,&dimLup,&dimRup,&dimLdown,&alpha,LblockL,&dimLup,temp2,&dimLdown,&beta,memHeff+denS->gKappa2index(ikappa),&dimLup);
1170  }
1171  }
1172  }
1173  }
1174  }
1175  }
1176  }
1177  }
1178 
1179  //5D3
1180  #ifdef CHEMPS2_MPI_COMPILATION
1181  if (( MPIchemps2::owner_specific_diagram( Prob->gL(), MPI_CHEMPS2_5D3 ) == MPIRANK ) && (N1==0) && (N2==2)){
1182  #else
1183  if ((N1==0) && (N2==2)){
1184  #endif
1185 
1186  for (int TwoSLdown=TwoSL-1; TwoSLdown<=TwoSL+1; TwoSLdown+=2){
1187  for (int TwoSRdown=TwoSR-1; TwoSRdown<=TwoSR+1; TwoSRdown+=2){
1188  for (int TwoJdown=0; TwoJdown<=2; TwoJdown+=2){
1189  if ((abs(TwoSLdown-TwoSRdown)<=TwoJdown) && (TwoSLdown>=0) && (TwoSRdown>=0)){
1190 
1191  int fase = phase(TwoSR + TwoSRdown + 3 + TwoJdown);
1192  double factor1 = fase * sqrt((TwoSLdown+1.0)*(TwoJdown+1.0)*(TwoSRdown+1.0)/(TwoSR+1.0)) * Wigner::wigner6j( 1, 1, TwoJdown, TwoSLdown, TwoSRdown, TwoSR );
1193  double factor2 = (TwoJdown==0)?sqrt(2.0*(TwoSRdown+1.0)/(TwoSR+1.0)):0.0;
1194 
1195  for (int Irrep=0; Irrep<(denBK->getNumberOfIrreps()); Irrep++){
1196 
1197  int IrrepTimesMid = Irreps::directProd(Irrep,IprodMID);
1198 
1199  bool leftOK = false;
1200  for (int l_alpha=0; l_alpha<theindex; l_alpha++){
1201  if (denBK->gIrrep(l_alpha) == Irrep){ leftOK = true; }
1202  }
1203  bool rightOK = false;
1204  for (int l_beta=theindex+2; l_beta<Prob->gL(); l_beta++){
1205  if (denBK->gIrrep(l_beta) == IrrepTimesMid){ rightOK = true; }
1206  }
1207 
1208  if ((leftOK) && (rightOK)){
1209 
1210  for (int l_alpha=0; l_alpha<theindex; l_alpha++){
1211  if (denBK->gIrrep(l_alpha) == Irrep){
1212 
1213  int ILdown = Irreps::directProd(IL, Irrep);
1214  int IRdown = Irreps::directProd(IR, IrrepTimesMid);
1215 
1216  int dimLdown = denBK->gCurrentDim(theindex, NL+1, TwoSLdown, ILdown);
1217  int dimRdown = denBK->gCurrentDim(theindex+2, NR+1, TwoSRdown, IRdown);
1218 
1219  if ((dimLdown>0) && (dimRdown>0)){
1220 
1221  int size = dimRup * dimRdown;
1222  for (int cnt=0; cnt<size; cnt++){ temp[cnt] = 0.0; }
1223 
1224  for (int l_beta=theindex+2; l_beta<Prob->gL(); l_beta++){
1225  if (denBK->gIrrep(l_beta) == IrrepTimesMid){
1226  double prefac = factor1 * Prob->gMxElement(l_alpha, theindex, theindex+1, l_beta)
1227  + factor2 * Prob->gMxElement(l_alpha, theindex, l_beta, theindex+1);
1228  double * LblockR = Lright[l_beta-theindex-2]->gStorage(NR,TwoSR,IR,NR+1,TwoSRdown,IRdown);
1229  daxpy_(&size,&prefac,LblockR,&inc,temp,&inc);
1230  }
1231  }
1232 
1233  double alpha = 1.0;
1234  double beta = 0.0; //set
1235  int memSkappa = denS->gKappa(NL+1, TwoSLdown, ILdown, 1, 1, TwoJdown, NR+1, TwoSRdown, IRdown);
1236 
1237  dgemm_(&notrans,&trans,&dimLdown,&dimRup,&dimRdown,&alpha,memS+denS->gKappa2index(memSkappa),&dimLdown,temp,&dimRup,&beta,temp2,&dimLdown);
1238 
1239  beta = 1.0; //add
1240  double * LblockL = Lleft[theindex-1-l_alpha]->gStorage(NL,TwoSL,IL,NL+1,TwoSLdown,ILdown);
1241  dgemm_(&notrans,&notrans,&dimLup,&dimRup,&dimLdown,&alpha,LblockL,&dimLup,temp2,&dimLdown,&beta,memHeff+denS->gKappa2index(ikappa),&dimLup);
1242  }
1243  }
1244  }
1245  }
1246  }
1247  }
1248  }
1249  }
1250  }
1251  }
1252 
1253  //5D4
1254  #ifdef CHEMPS2_MPI_COMPILATION
1255  if (( MPIchemps2::owner_specific_diagram( Prob->gL(), MPI_CHEMPS2_5D4 ) == MPIRANK ) && (N1==1) && (N2==2)){
1256  #else
1257  if ((N1==1) && (N2==2)){
1258  #endif
1259 
1260  for (int TwoSLdown=TwoSL-1; TwoSLdown<=TwoSL+1; TwoSLdown+=2){
1261  for (int TwoSRdown=TwoSR-1; TwoSRdown<=TwoSR+1; TwoSRdown+=2){
1262  if ((abs(TwoSLdown-TwoSRdown)<=1) && (TwoSLdown>=0) && (TwoSRdown>=0)){
1263 
1264  const double factor1 = (TwoSL==TwoSRdown) ? phase(TwoSLdown-TwoSR) * sqrt((TwoSLdown+1.0)/(TwoSL+1.0)) : 0.0;
1265  const double factor2 = phase(TwoSLdown+TwoSR+2) * sqrt((TwoSLdown+1)*(TwoSRdown+1.0)) * Wigner::wigner6j( TwoSL, TwoSR, 1, TwoSRdown, TwoSLdown, 1 );
1266 
1267  for (int Irrep=0; Irrep<(denBK->getNumberOfIrreps()); Irrep++){
1268 
1269  int IrrepTimesMid = Irreps::directProd(Irrep,IprodMID);
1270 
1271  bool leftOK = false;
1272  for (int l_alpha=0; l_alpha<theindex; l_alpha++){
1273  if (denBK->gIrrep(l_alpha) == Irrep){ leftOK = true; }
1274  }
1275  bool rightOK = false;
1276  for (int l_beta=theindex+2; l_beta<Prob->gL(); l_beta++){
1277  if (denBK->gIrrep(l_beta) == IrrepTimesMid){ rightOK = true; }
1278  }
1279 
1280  if ((leftOK) && (rightOK)){
1281 
1282  for (int l_alpha=0; l_alpha<theindex; l_alpha++){
1283  if (denBK->gIrrep(l_alpha) == Irrep){
1284 
1285  int ILdown = Irreps::directProd(IL, Irrep);
1286  int IRdown = Irreps::directProd(IR, IrrepTimesMid);
1287 
1288  int dimLdown = denBK->gCurrentDim(theindex, NL+1, TwoSLdown, ILdown);
1289  int dimRdown = denBK->gCurrentDim(theindex+2, NR+1, TwoSRdown, IRdown);
1290 
1291  if ((dimLdown>0) && (dimRdown>0)){
1292 
1293  int size = dimRup * dimRdown;
1294  for (int cnt=0; cnt<size; cnt++){ temp[cnt] = 0.0; }
1295 
1296  for (int l_beta=theindex+2; l_beta<Prob->gL(); l_beta++){
1297  if (denBK->gIrrep(l_beta) == IrrepTimesMid){
1298  double prefac = factor1 * Prob->gMxElement(l_alpha, theindex, theindex+1, l_beta)
1299  + factor2 * Prob->gMxElement(l_alpha, theindex, l_beta, theindex+1);
1300  double * LblockR = Lright[l_beta-theindex-2]->gStorage(NR,TwoSR,IR,NR+1,TwoSRdown,IRdown);
1301  daxpy_(&size,&prefac,LblockR,&inc,temp,&inc);
1302  }
1303  }
1304 
1305  double alpha = 1.0;
1306  double beta = 0.0; //set
1307  int memSkappa = denS->gKappa(NL+1, TwoSLdown, ILdown, 2, 1, 1, NR+1, TwoSRdown, IRdown);
1308 
1309  dgemm_(&notrans,&trans,&dimLdown,&dimRup,&dimRdown,&alpha,memS+denS->gKappa2index(memSkappa),&dimLdown,temp,&dimRup,&beta,temp2,&dimLdown);
1310 
1311  beta = 1.0; //add
1312  double * LblockL = Lleft[theindex-1-l_alpha]->gStorage(NL,TwoSL,IL,NL+1,TwoSLdown,ILdown);
1313  dgemm_(&notrans,&notrans,&dimLup,&dimRup,&dimLdown,&alpha,LblockL,&dimLup,temp2,&dimLdown,&beta,memHeff+denS->gKappa2index(ikappa),&dimLup);
1314  }
1315  }
1316  }
1317  }
1318  }
1319  }
1320  }
1321  }
1322  }
1323 
1324 }
1325 
1326 void CheMPS2::Heff::addDiagram5E(const int ikappa, double * memS, double * memHeff, const Sobject * denS, TensorL ** Lleft, TensorL ** Lright, double * temp, double * temp2) const{
1327 
1328  #ifdef CHEMPS2_MPI_COMPILATION
1329  const int MPIRANK = MPIchemps2::mpi_rank();
1330  #endif
1331 
1332  int NL = denS->gNL(ikappa);
1333  int TwoSL = denS->gTwoSL(ikappa);
1334  int IL = denS->gIL(ikappa);
1335 
1336  int NR = denS->gNR(ikappa);
1337  int TwoSR = denS->gTwoSR(ikappa);
1338  int IR = denS->gIR(ikappa);
1339 
1340  int N1 = denS->gN1(ikappa);
1341  int N2 = denS->gN2(ikappa);
1342  int TwoJ = denS->gTwoJ(ikappa);
1343 
1344  int theindex = denS->gIndex();
1345  int dimLup = denBK->gCurrentDim(theindex ,NL,TwoSL,IL);
1346  int dimRup = denBK->gCurrentDim(theindex+2,NR,TwoSR,IR);
1347 
1348  char notrans = 'N';
1349  char trans = 'T';
1350  int inc = 1;
1351 
1352  int IprodMID = Irreps::directProd(denBK->gIrrep(theindex),denBK->gIrrep(theindex+1));
1353 
1354  //5E1
1355  #ifdef CHEMPS2_MPI_COMPILATION
1356  if (( MPIchemps2::owner_specific_diagram( Prob->gL(), MPI_CHEMPS2_5E1 ) == MPIRANK ) && (N1==0) && (N2==1)){
1357  #else
1358  if ((N1==0) && (N2==1)){
1359  #endif
1360 
1361  for (int TwoSLdown=TwoSL-1; TwoSLdown<=TwoSL+1; TwoSLdown+=2){
1362  for (int TwoSRdown=TwoSR-1; TwoSRdown<=TwoSR+1; TwoSRdown+=2){
1363  if ((abs(TwoSLdown-TwoSRdown)<=1) && (TwoSLdown>=0) && (TwoSRdown>=0)){
1364 
1365  const double factor2 = phase(TwoSL+TwoSRdown) * sqrt((TwoSL+1)*(TwoSR+1.0)) * Wigner::wigner6j( TwoSLdown, TwoSRdown, 1, TwoSR, TwoSL, 1 );
1366  const double factor1 = (TwoSL==TwoSRdown)?sqrt((TwoSR+1.0)/(TwoSRdown+1.0)):0.0;
1367 
1368  for (int Irrep=0; Irrep<(denBK->getNumberOfIrreps()); Irrep++){
1369 
1370  int IrrepTimesMid = Irreps::directProd(Irrep,IprodMID);
1371 
1372  bool leftOK = false;
1373  for (int l_alpha=0; l_alpha<theindex; l_alpha++){
1374  if (denBK->gIrrep(l_alpha) == Irrep){ leftOK = true; }
1375  }
1376  bool rightOK = false;
1377  for (int l_beta=theindex+2; l_beta<Prob->gL(); l_beta++){
1378  if (denBK->gIrrep(l_beta) == IrrepTimesMid){ rightOK = true; }
1379  }
1380 
1381  if ((leftOK) && (rightOK)){
1382 
1383  for (int l_alpha=0; l_alpha<theindex; l_alpha++){
1384  if (denBK->gIrrep(l_alpha) == Irrep){
1385 
1386  int ILdown = Irreps::directProd(IL, Irrep);
1387  int IRdown = Irreps::directProd(IR, IrrepTimesMid);
1388 
1389  int dimLdown = denBK->gCurrentDim(theindex, NL-1, TwoSLdown, ILdown);
1390  int dimRdown = denBK->gCurrentDim(theindex+2, NR-1, TwoSRdown, IRdown);
1391 
1392  if ((dimLdown>0) && (dimRdown>0)){
1393 
1394  int size = dimRup * dimRdown;
1395  for (int cnt=0; cnt<size; cnt++){ temp[cnt] = 0.0; }
1396 
1397  for (int l_beta=theindex+2; l_beta<Prob->gL(); l_beta++){
1398  if (denBK->gIrrep(l_beta) == IrrepTimesMid){
1399  double prefac = factor1 * Prob->gMxElement(l_alpha,theindex+1,theindex,l_beta)
1400  + factor2 * Prob->gMxElement(l_alpha,theindex+1,l_beta,theindex);
1401  double * LblockR = Lright[l_beta-theindex-2]->gStorage(NR-1,TwoSRdown,IRdown,NR,TwoSR,IR);
1402  daxpy_(&size,&prefac,LblockR,&inc,temp,&inc);
1403  }
1404  }
1405 
1406  double alpha = 1.0;
1407  double beta = 0.0; //set
1408  int memSkappa = denS->gKappa(NL-1, TwoSLdown, ILdown, 1, 0, 1, NR-1, TwoSRdown, IRdown);
1409 
1410  dgemm_(&notrans,&notrans,&dimLdown,&dimRup,&dimRdown,&alpha,memS+denS->gKappa2index(memSkappa),&dimLdown,temp,&dimRdown,&beta,temp2,&dimLdown);
1411 
1412  beta = 1.0; //add
1413  double * LblockL = Lleft[theindex-1-l_alpha]->gStorage(NL-1,TwoSLdown,ILdown,NL,TwoSL,IL);
1414  dgemm_(&trans,&notrans,&dimLup,&dimRup,&dimLdown,&alpha,LblockL,&dimLdown,temp2,&dimLdown,&beta,memHeff+denS->gKappa2index(ikappa),&dimLup);
1415  }
1416  }
1417  }
1418  }
1419  }
1420  }
1421  }
1422  }
1423  }
1424 
1425  //5E2
1426  #ifdef CHEMPS2_MPI_COMPILATION
1427  if (( MPIchemps2::owner_specific_diagram( Prob->gL(), MPI_CHEMPS2_5E2 ) == MPIRANK ) && (N1==0) && (N2==2)){
1428  #else
1429  if ((N1==0) && (N2==2)){
1430  #endif
1431 
1432  for (int TwoSLdown=TwoSL-1; TwoSLdown<=TwoSL+1; TwoSLdown+=2){
1433  for (int TwoSRdown=TwoSR-1; TwoSRdown<=TwoSR+1; TwoSRdown+=2){
1434  for (int TwoJdown=0; TwoJdown<=2; TwoJdown+=2){
1435  if ((abs(TwoSLdown-TwoSRdown)<=TwoJdown) && (TwoSLdown>=0) && (TwoSRdown>=0)){
1436 
1437  int fase = phase(TwoSR + TwoSLdown + 3);
1438  const double factor1 = fase * sqrt((TwoSR+1)*(TwoJdown+1.0)) * Wigner::wigner6j( 1, 1, TwoJdown, TwoSLdown, TwoSRdown, TwoSR );
1439  const double factor2 = (TwoJdown==0)?sqrt(2.0*(TwoSR+1.0)/(TwoSRdown+1.0)):0.0;
1440 
1441  for (int Irrep=0; Irrep<(denBK->getNumberOfIrreps()); Irrep++){
1442 
1443  int IrrepTimesMid = Irreps::directProd(Irrep,IprodMID);
1444 
1445  bool leftOK = false;
1446  for (int l_alpha=0; l_alpha<theindex; l_alpha++){
1447  if (denBK->gIrrep(l_alpha) == Irrep){ leftOK = true; }
1448  }
1449  bool rightOK = false;
1450  for (int l_beta=theindex+2; l_beta<Prob->gL(); l_beta++){
1451  if (denBK->gIrrep(l_beta) == IrrepTimesMid){ rightOK = true; }
1452  }
1453 
1454  if ((leftOK) && (rightOK)){
1455 
1456  for (int l_alpha=0; l_alpha<theindex; l_alpha++){
1457  if (denBK->gIrrep(l_alpha) == Irrep){
1458 
1459  int ILdown = Irreps::directProd(IL, Irrep);
1460  int IRdown = Irreps::directProd(IR, IrrepTimesMid);
1461 
1462  int dimLdown = denBK->gCurrentDim(theindex, NL-1, TwoSLdown, ILdown);
1463  int dimRdown = denBK->gCurrentDim(theindex+2, NR-1, TwoSRdown, IRdown);
1464 
1465  if ((dimLdown>0) && (dimRdown>0)){
1466 
1467  int size = dimRup * dimRdown;
1468  for (int cnt=0; cnt<size; cnt++){ temp[cnt] = 0.0; }
1469 
1470  for (int l_beta=theindex+2; l_beta<Prob->gL(); l_beta++){
1471  if (denBK->gIrrep(l_beta) == IrrepTimesMid){
1472  double prefac = factor1 * Prob->gMxElement(l_alpha, theindex+1, theindex, l_beta)
1473  + factor2 * Prob->gMxElement(l_alpha, theindex+1, l_beta, theindex);
1474  double * LblockR = Lright[l_beta-theindex-2]->gStorage(NR-1,TwoSRdown,IRdown,NR,TwoSR,IR);
1475  daxpy_(&size,&prefac,LblockR,&inc,temp,&inc);
1476  }
1477  }
1478 
1479  double alpha = 1.0;
1480  double beta = 0.0; //set
1481  int memSkappa = denS->gKappa(NL-1, TwoSLdown, ILdown, 1, 1, TwoJdown, NR-1, TwoSRdown, IRdown);
1482 
1483  dgemm_(&notrans,&notrans,&dimLdown,&dimRup,&dimRdown,&alpha,memS+denS->gKappa2index(memSkappa),&dimLdown,temp,&dimRdown,&beta,temp2,&dimLdown);
1484 
1485  beta = 1.0; //add
1486  double * LblockL = Lleft[theindex-1-l_alpha]->gStorage(NL-1,TwoSLdown,ILdown,NL,TwoSL,IL);
1487  dgemm_(&trans,&notrans,&dimLup,&dimRup,&dimLdown,&alpha,LblockL,&dimLdown,temp2,&dimLdown,&beta,memHeff+denS->gKappa2index(ikappa),&dimLup);
1488  }
1489  }
1490  }
1491  }
1492  }
1493  }
1494  }
1495  }
1496  }
1497  }
1498 
1499  //5E3
1500  #ifdef CHEMPS2_MPI_COMPILATION
1501  if (( MPIchemps2::owner_specific_diagram( Prob->gL(), MPI_CHEMPS2_5E3 ) == MPIRANK ) && (N1==1) && (N2==1)){
1502  #else
1503  if ((N1==1) && (N2==1)){
1504  #endif
1505 
1506  for (int TwoSLdown=TwoSL-1; TwoSLdown<=TwoSL+1; TwoSLdown+=2){
1507  if ((TwoSLdown>=0) && (abs(TwoSR-TwoSLdown)<=1)){
1508  int TwoSRdown = TwoSLdown;
1509 
1510  int fase = phase(TwoSR + TwoSRdown + 3);
1511  double factor1 = fase * sqrt((TwoSL+1.0)*(TwoJ+1.0)*(TwoSR+1.0)/(TwoSRdown+1.0)) * Wigner::wigner6j( 1, 1, TwoJ, TwoSL, TwoSR, TwoSRdown );
1512  double factor2 = (TwoJ==0)?sqrt(2.0*(TwoSR+1.0)/(TwoSRdown+1.0)):0.0;
1513 
1514  for (int Irrep=0; Irrep<(denBK->getNumberOfIrreps()); Irrep++){
1515 
1516  int IrrepTimesMid = Irreps::directProd(Irrep,IprodMID);
1517 
1518  bool leftOK = false;
1519  for (int l_alpha=0; l_alpha<theindex; l_alpha++){
1520  if (denBK->gIrrep(l_alpha) == Irrep){ leftOK = true; }
1521  }
1522  bool rightOK = false;
1523  for (int l_beta=theindex+2; l_beta<Prob->gL(); l_beta++){
1524  if (denBK->gIrrep(l_beta) == IrrepTimesMid){ rightOK = true; }
1525  }
1526 
1527  if ((leftOK) && (rightOK)){
1528 
1529  for (int l_alpha=0; l_alpha<theindex; l_alpha++){
1530  if (denBK->gIrrep(l_alpha) == Irrep){
1531 
1532  int ILdown = Irreps::directProd(IL, Irrep);
1533  int IRdown = Irreps::directProd(IR, IrrepTimesMid);
1534 
1535  int dimLdown = denBK->gCurrentDim(theindex, NL-1, TwoSLdown, ILdown);
1536  int dimRdown = denBK->gCurrentDim(theindex+2, NR-1, TwoSRdown, IRdown);
1537 
1538  if ((dimLdown>0) && (dimRdown>0)){
1539 
1540  int size = dimRup * dimRdown;
1541  for (int cnt=0; cnt<size; cnt++){ temp[cnt] = 0.0; }
1542 
1543  for (int l_beta=theindex+2; l_beta<Prob->gL(); l_beta++){
1544  if (denBK->gIrrep(l_beta) == IrrepTimesMid){
1545  double prefac = factor1 * Prob->gMxElement(l_alpha, theindex+1, theindex, l_beta)
1546  + factor2 * Prob->gMxElement(l_alpha, theindex+1, l_beta, theindex);
1547  double * LblockR = Lright[l_beta-theindex-2]->gStorage(NR-1,TwoSRdown,IRdown,NR,TwoSR,IR);
1548  daxpy_(&size,&prefac,LblockR,&inc,temp,&inc);
1549  }
1550  }
1551 
1552  double alpha = 1.0;
1553  double beta = 0.0; //set
1554  int memSkappa = denS->gKappa(NL-1, TwoSLdown, ILdown, 2, 0, 0, NR-1, TwoSRdown, IRdown);
1555 
1556  dgemm_(&notrans,&notrans,&dimLdown,&dimRup,&dimRdown,&alpha,memS+denS->gKappa2index(memSkappa),&dimLdown,temp,&dimRdown,&beta,temp2,&dimLdown);
1557 
1558  beta = 1.0; //add
1559  double * LblockL = Lleft[theindex-1-l_alpha]->gStorage(NL-1,TwoSLdown,ILdown,NL,TwoSL,IL);
1560  dgemm_(&trans,&notrans,&dimLup,&dimRup,&dimLdown,&alpha,LblockL,&dimLdown,temp2,&dimLdown,&beta,memHeff+denS->gKappa2index(ikappa),&dimLup);
1561  }
1562  }
1563  }
1564  }
1565  }
1566  }
1567  }
1568  }
1569 
1570  //5E4
1571  #ifdef CHEMPS2_MPI_COMPILATION
1572  if (( MPIchemps2::owner_specific_diagram( Prob->gL(), MPI_CHEMPS2_5E4 ) == MPIRANK ) && (N1==1) && (N2==2)){
1573  #else
1574  if ((N1==1) && (N2==2)){
1575  #endif
1576 
1577  for (int TwoSLdown=TwoSL-1; TwoSLdown<=TwoSL+1; TwoSLdown+=2){
1578  for (int TwoSRdown=TwoSR-1; TwoSRdown<=TwoSR+1; TwoSRdown+=2){
1579  if ((abs(TwoSLdown-TwoSRdown)<=1) && (TwoSLdown>=0) && (TwoSRdown>=0)){
1580 
1581  const double factor1 = (TwoSLdown==TwoSR) ? phase(TwoSL-TwoSRdown) * sqrt((TwoSL+1.0)/(TwoSLdown+1.0)) : 0.0;
1582  const double factor2 = phase(TwoSL+TwoSRdown+2) * sqrt((TwoSL+1)*(TwoSR+1.0)) * Wigner::wigner6j( TwoSLdown, TwoSRdown, 1, TwoSR, TwoSL, 1 );
1583 
1584  for (int Irrep=0; Irrep<(denBK->getNumberOfIrreps()); Irrep++){
1585 
1586  int IrrepTimesMid = Irreps::directProd(Irrep,IprodMID);
1587 
1588  bool leftOK = false;
1589  for (int l_alpha=0; l_alpha<theindex; l_alpha++){
1590  if (denBK->gIrrep(l_alpha) == Irrep){ leftOK = true; }
1591  }
1592  bool rightOK = false;
1593  for (int l_beta=theindex+2; l_beta<Prob->gL(); l_beta++){
1594  if (denBK->gIrrep(l_beta) == IrrepTimesMid){ rightOK = true; }
1595  }
1596 
1597  if ((leftOK) && (rightOK)){
1598 
1599  for (int l_alpha=0; l_alpha<theindex; l_alpha++){
1600  if (denBK->gIrrep(l_alpha) == Irrep){
1601 
1602  int ILdown = Irreps::directProd(IL, Irrep);
1603  int IRdown = Irreps::directProd(IR, IrrepTimesMid);
1604 
1605  int dimLdown = denBK->gCurrentDim(theindex, NL-1, TwoSLdown, ILdown);
1606  int dimRdown = denBK->gCurrentDim(theindex+2, NR-1, TwoSRdown, IRdown);
1607 
1608  if ((dimLdown>0) && (dimRdown>0)){
1609 
1610  int size = dimRup * dimRdown;
1611  for (int cnt=0; cnt<size; cnt++){ temp[cnt] = 0.0; }
1612 
1613  for (int l_beta=theindex+2; l_beta<Prob->gL(); l_beta++){
1614  if (denBK->gIrrep(l_beta) == IrrepTimesMid){
1615  double prefac = factor1 * Prob->gMxElement(l_alpha, theindex+1, theindex, l_beta)
1616  + factor2 * Prob->gMxElement(l_alpha, theindex+1, l_beta, theindex);
1617  double * LblockR = Lright[l_beta-theindex-2]->gStorage(NR-1,TwoSRdown,IRdown,NR,TwoSR,IR);
1618  daxpy_(&size,&prefac,LblockR,&inc,temp,&inc);
1619  }
1620  }
1621 
1622  double alpha = 1.0;
1623  double beta = 0.0; //set
1624  int memSkappa = denS->gKappa(NL-1, TwoSLdown, ILdown, 2, 1, 1, NR-1, TwoSRdown, IRdown);
1625 
1626  dgemm_(&notrans,&notrans,&dimLdown,&dimRup,&dimRdown,&alpha,memS+denS->gKappa2index(memSkappa),&dimLdown,temp,&dimRdown,&beta,temp2,&dimLdown);
1627 
1628  beta = 1.0; //add
1629  double * LblockL = Lleft[theindex-1-l_alpha]->gStorage(NL-1,TwoSLdown,ILdown,NL,TwoSL,IL);
1630  dgemm_(&trans,&notrans,&dimLup,&dimRup,&dimLdown,&alpha,LblockL,&dimLdown,temp2,&dimLdown,&beta,memHeff+denS->gKappa2index(ikappa),&dimLup);
1631  }
1632  }
1633  }
1634  }
1635  }
1636  }
1637  }
1638  }
1639  }
1640 
1641 }
1642 
1643 void CheMPS2::Heff::addDiagram5F(const int ikappa, double * memS, double * memHeff, const Sobject * denS, TensorL ** Lleft, TensorL ** Lright, double * temp, double * temp2) const{
1644 
1645  #ifdef CHEMPS2_MPI_COMPILATION
1646  const int MPIRANK = MPIchemps2::mpi_rank();
1647  #endif
1648 
1649  int NL = denS->gNL(ikappa);
1650  int TwoSL = denS->gTwoSL(ikappa);
1651  int IL = denS->gIL(ikappa);
1652 
1653  int NR = denS->gNR(ikappa);
1654  int TwoSR = denS->gTwoSR(ikappa);
1655  int IR = denS->gIR(ikappa);
1656 
1657  int N1 = denS->gN1(ikappa);
1658  int N2 = denS->gN2(ikappa);
1659  int TwoJ = denS->gTwoJ(ikappa);
1660 
1661  int theindex = denS->gIndex();
1662  int dimLup = denBK->gCurrentDim(theindex ,NL,TwoSL,IL);
1663  int dimRup = denBK->gCurrentDim(theindex+2,NR,TwoSR,IR);
1664 
1665  char notrans = 'N';
1666  char trans = 'T';
1667  int inc = 1;
1668 
1669  int IprodMID = Irreps::directProd(denBK->gIrrep(theindex),denBK->gIrrep(theindex+1));
1670 
1671  //5F1
1672  #ifdef CHEMPS2_MPI_COMPILATION
1673  if (( MPIchemps2::owner_specific_diagram( Prob->gL(), MPI_CHEMPS2_5F1 ) == MPIRANK ) && (N1==1) && (N2==0)){
1674  #else
1675  if ((N1==1) && (N2==0)){
1676  #endif
1677 
1678  for (int TwoSLdown=TwoSL-1; TwoSLdown<=TwoSL+1; TwoSLdown+=2){
1679  for (int TwoSRdown=TwoSR-1; TwoSRdown<=TwoSR+1; TwoSRdown+=2){
1680  if ((abs(TwoSLdown-TwoSRdown)<=1) && (TwoSLdown>=0) && (TwoSRdown>=0)){
1681 
1682  const double factor2 = phase(TwoSLdown+TwoSR) * sqrt((TwoSLdown+1)*(TwoSRdown+1.0)) * Wigner::wigner6j( TwoSL, TwoSR, 1, TwoSRdown, TwoSLdown, 1 );
1683  const double factor1 = (TwoSLdown==TwoSR)?sqrt((TwoSRdown+1.0)/(TwoSR+1.0)):0.0;
1684 
1685  for (int Irrep=0; Irrep<(denBK->getNumberOfIrreps()); Irrep++){
1686 
1687  int IrrepTimesMid = Irreps::directProd(Irrep,IprodMID);
1688 
1689  bool leftOK = false;
1690  for (int l_alpha=0; l_alpha<theindex; l_alpha++){
1691  if (denBK->gIrrep(l_alpha) == Irrep){ leftOK = true; }
1692  }
1693  bool rightOK = false;
1694  for (int l_beta=theindex+2; l_beta<Prob->gL(); l_beta++){
1695  if (denBK->gIrrep(l_beta) == IrrepTimesMid){ rightOK = true; }
1696  }
1697 
1698  if ((leftOK) && (rightOK)){
1699 
1700  for (int l_alpha=0; l_alpha<theindex; l_alpha++){
1701  if (denBK->gIrrep(l_alpha) == Irrep){
1702 
1703  int ILdown = Irreps::directProd(IL, Irrep);
1704  int IRdown = Irreps::directProd(IR, IrrepTimesMid);
1705 
1706  int dimLdown = denBK->gCurrentDim(theindex, NL+1, TwoSLdown, ILdown);
1707  int dimRdown = denBK->gCurrentDim(theindex+2, NR+1, TwoSRdown, IRdown);
1708 
1709  if ((dimLdown>0) && (dimRdown>0)){
1710 
1711  int size = dimRup * dimRdown;
1712  for (int cnt=0; cnt<size; cnt++){ temp[cnt] = 0.0; }
1713 
1714  for (int l_beta=theindex+2; l_beta<Prob->gL(); l_beta++){
1715  if (denBK->gIrrep(l_beta) == IrrepTimesMid){
1716  double prefac = factor1 * Prob->gMxElement(l_alpha,theindex+1,theindex,l_beta)
1717  + factor2 * Prob->gMxElement(l_alpha,theindex+1,l_beta,theindex);
1718  double * LblockR = Lright[l_beta-theindex-2]->gStorage(NR,TwoSR,IR,NR+1,TwoSRdown,IRdown);
1719  daxpy_(&size,&prefac,LblockR,&inc,temp,&inc);
1720  }
1721  }
1722 
1723  double alpha = 1.0;
1724  double beta = 0.0; //set
1725  int memSkappa = denS->gKappa(NL+1, TwoSLdown, ILdown, 0, 1, 1, NR+1, TwoSRdown, IRdown);
1726 
1727  dgemm_(&notrans,&trans,&dimLdown,&dimRup,&dimRdown,&alpha,memS+denS->gKappa2index(memSkappa),&dimLdown,temp,&dimRup,&beta,temp2,&dimLdown);
1728 
1729  beta = 1.0; //add
1730  double * LblockL = Lleft[theindex-1-l_alpha]->gStorage(NL,TwoSL,IL,NL+1,TwoSLdown,ILdown);
1731  dgemm_(&notrans,&notrans,&dimLup,&dimRup,&dimLdown,&alpha,LblockL,&dimLup,temp2,&dimLdown,&beta,memHeff+denS->gKappa2index(ikappa),&dimLup);
1732  }
1733  }
1734  }
1735  }
1736  }
1737  }
1738  }
1739  }
1740  }
1741 
1742  //5F2
1743  #ifdef CHEMPS2_MPI_COMPILATION
1744  if (( MPIchemps2::owner_specific_diagram( Prob->gL(), MPI_CHEMPS2_5F2 ) == MPIRANK ) && (N1==1) && (N2==1)){
1745  #else
1746  if ((N1==1) && (N2==1)){
1747  #endif
1748 
1749  for (int TwoSLdown=TwoSL-1; TwoSLdown<=TwoSL+1; TwoSLdown+=2){
1750  if ((TwoSLdown>=0) && (abs(TwoSR-TwoSLdown)<=1)){
1751  int TwoSRdown = TwoSLdown;
1752 
1753  const double factor1 = phase(TwoSRdown + TwoSL + 3) * sqrt((TwoSRdown+1)*(TwoJ+1.0)) * Wigner::wigner6j( 1, 1, TwoJ, TwoSL, TwoSR, TwoSRdown );
1754  const double factor2 = (TwoJ==0)?sqrt(2.0*(TwoSRdown+1.0)/(TwoSR+1.0)):0.0;
1755 
1756  for (int Irrep=0; Irrep<(denBK->getNumberOfIrreps()); Irrep++){
1757 
1758  int IrrepTimesMid = Irreps::directProd(Irrep,IprodMID);
1759 
1760  bool leftOK = false;
1761  for (int l_alpha=0; l_alpha<theindex; l_alpha++){
1762  if (denBK->gIrrep(l_alpha) == Irrep){ leftOK = true; }
1763  }
1764  bool rightOK = false;
1765  for (int l_beta=theindex+2; l_beta<Prob->gL(); l_beta++){
1766  if (denBK->gIrrep(l_beta) == IrrepTimesMid){ rightOK = true; }
1767  }
1768 
1769  if ((leftOK) && (rightOK)){
1770 
1771  for (int l_alpha=0; l_alpha<theindex; l_alpha++){
1772  if (denBK->gIrrep(l_alpha) == Irrep){
1773 
1774  int ILdown = Irreps::directProd(IL, Irrep);
1775  int IRdown = Irreps::directProd(IR, IrrepTimesMid);
1776 
1777  int dimLdown = denBK->gCurrentDim(theindex, NL+1, TwoSLdown, ILdown);
1778  int dimRdown = denBK->gCurrentDim(theindex+2, NR+1, TwoSRdown, IRdown);
1779 
1780  if ((dimLdown>0) && (dimRdown>0)){
1781 
1782  int size = dimRup * dimRdown;
1783  for (int cnt=0; cnt<size; cnt++){ temp[cnt] = 0.0; }
1784 
1785  for (int l_beta=theindex+2; l_beta<Prob->gL(); l_beta++){
1786  if (denBK->gIrrep(l_beta) == IrrepTimesMid){
1787  double prefac = factor1 * Prob->gMxElement(l_alpha, theindex+1, theindex, l_beta)
1788  + factor2 * Prob->gMxElement(l_alpha, theindex+1, l_beta, theindex);
1789  double * LblockR = Lright[l_beta-theindex-2]->gStorage(NR,TwoSR,IR,NR+1,TwoSRdown,IRdown);
1790  daxpy_(&size,&prefac,LblockR,&inc,temp,&inc);
1791  }
1792  }
1793 
1794  double alpha = 1.0;
1795  double beta = 0.0; //set
1796  int memSkappa = denS->gKappa(NL+1, TwoSLdown, ILdown, 0, 2, 0, NR+1, TwoSRdown, IRdown);
1797 
1798  dgemm_(&notrans,&trans,&dimLdown,&dimRup,&dimRdown,&alpha,memS+denS->gKappa2index(memSkappa),&dimLdown,temp,&dimRup,&beta,temp2,&dimLdown);
1799 
1800  beta = 1.0; //add
1801  double * LblockL = Lleft[theindex-1-l_alpha]->gStorage(NL,TwoSL,IL,NL+1,TwoSLdown,ILdown);
1802  dgemm_(&notrans,&notrans,&dimLup,&dimRup,&dimLdown,&alpha,LblockL,&dimLup,temp2,&dimLdown,&beta,memHeff+denS->gKappa2index(ikappa),&dimLup);
1803  }
1804  }
1805  }
1806  }
1807  }
1808  }
1809  }
1810  }
1811 
1812  //5F3
1813  #ifdef CHEMPS2_MPI_COMPILATION
1814  if (( MPIchemps2::owner_specific_diagram( Prob->gL(), MPI_CHEMPS2_5F3 ) == MPIRANK ) && (N1==2) && (N2==0)){
1815  #else
1816  if ((N1==2) && (N2==0)){
1817  #endif
1818 
1819  for (int TwoSLdown=TwoSL-1; TwoSLdown<=TwoSL+1; TwoSLdown+=2){
1820  for (int TwoSRdown=TwoSR-1; TwoSRdown<=TwoSR+1; TwoSRdown+=2){
1821  for (int TwoJdown=0; TwoJdown<=2; TwoJdown+=2){
1822  if ((abs(TwoSLdown-TwoSRdown)<=TwoJdown) && (TwoSLdown>=0) && (TwoSRdown>=0)){
1823 
1824  double factor1 = phase(TwoSR + TwoSRdown + 3) * sqrt((TwoSLdown+1.0)*(TwoJdown+1.0)*(TwoSRdown+1.0)/(TwoSR+1.0)) * Wigner::wigner6j( 1, 1, TwoJdown, TwoSLdown, TwoSRdown, TwoSR );
1825  double factor2 = (TwoJdown==0)?sqrt(2.0*(TwoSRdown+1.0)/(TwoSR+1.0)):0.0;
1826 
1827  for (int Irrep=0; Irrep<(denBK->getNumberOfIrreps()); Irrep++){
1828 
1829  int IrrepTimesMid = Irreps::directProd(Irrep,IprodMID);
1830 
1831  bool leftOK = false;
1832  for (int l_alpha=0; l_alpha<theindex; l_alpha++){
1833  if (denBK->gIrrep(l_alpha) == Irrep){ leftOK = true; }
1834  }
1835  bool rightOK = false;
1836  for (int l_beta=theindex+2; l_beta<Prob->gL(); l_beta++){
1837  if (denBK->gIrrep(l_beta) == IrrepTimesMid){ rightOK = true; }
1838  }
1839 
1840  if ((leftOK) && (rightOK)){
1841 
1842  for (int l_alpha=0; l_alpha<theindex; l_alpha++){
1843  if (denBK->gIrrep(l_alpha) == Irrep){
1844 
1845  int ILdown = Irreps::directProd(IL, Irrep);
1846  int IRdown = Irreps::directProd(IR, IrrepTimesMid);
1847 
1848  int dimLdown = denBK->gCurrentDim(theindex, NL+1, TwoSLdown, ILdown);
1849  int dimRdown = denBK->gCurrentDim(theindex+2, NR+1, TwoSRdown, IRdown);
1850 
1851  if ((dimLdown>0) && (dimRdown>0)){
1852 
1853  int size = dimRup * dimRdown;
1854  for (int cnt=0; cnt<size; cnt++){ temp[cnt] = 0.0; }
1855 
1856  for (int l_beta=theindex+2; l_beta<Prob->gL(); l_beta++){
1857  if (denBK->gIrrep(l_beta) == IrrepTimesMid){
1858  double prefac = factor1 * Prob->gMxElement(l_alpha, theindex+1, theindex, l_beta)
1859  + factor2 * Prob->gMxElement(l_alpha, theindex+1, l_beta, theindex);
1860  double * LblockR = Lright[l_beta-theindex-2]->gStorage(NR,TwoSR,IR,NR+1,TwoSRdown,IRdown);
1861  daxpy_(&size,&prefac,LblockR,&inc,temp,&inc);
1862  }
1863  }
1864 
1865  double alpha = 1.0;
1866  double beta = 0.0; //set
1867  int memSkappa = denS->gKappa(NL+1, TwoSLdown, ILdown, 1, 1, TwoJdown, NR+1, TwoSRdown, IRdown);
1868 
1869  dgemm_(&notrans,&trans,&dimLdown,&dimRup,&dimRdown,&alpha,memS+denS->gKappa2index(memSkappa),&dimLdown,temp,&dimRup,&beta,temp2,&dimLdown);
1870 
1871  beta = 1.0; //add
1872  double * LblockL = Lleft[theindex-1-l_alpha]->gStorage(NL,TwoSL,IL,NL+1,TwoSLdown,ILdown);
1873  dgemm_(&notrans,&notrans,&dimLup,&dimRup,&dimLdown,&alpha,LblockL,&dimLup,temp2,&dimLdown,&beta,memHeff+denS->gKappa2index(ikappa),&dimLup);
1874  }
1875  }
1876  }
1877  }
1878  }
1879  }
1880  }
1881  }
1882  }
1883  }
1884 
1885  //5F4
1886  #ifdef CHEMPS2_MPI_COMPILATION
1887  if (( MPIchemps2::owner_specific_diagram( Prob->gL(), MPI_CHEMPS2_5F4 ) == MPIRANK ) && (N1==2) && (N2==1)){
1888  #else
1889  if ((N1==2) && (N2==1)){
1890  #endif
1891 
1892  for (int TwoSLdown=TwoSL-1; TwoSLdown<=TwoSL+1; TwoSLdown+=2){
1893  for (int TwoSRdown=TwoSR-1; TwoSRdown<=TwoSR+1; TwoSRdown+=2){
1894  if ((abs(TwoSLdown-TwoSRdown)<=1) && (TwoSLdown>=0) && (TwoSRdown>=0)){
1895 
1896  const double factor1 = (TwoSL==TwoSRdown) ? phase(TwoSLdown-TwoSR) * sqrt((TwoSLdown+1.0)/(TwoSL+1.0)) : 0.0;
1897  const double factor2 = phase(TwoSLdown+TwoSR+2) * sqrt((TwoSLdown+1)*(TwoSRdown+1.0)) * Wigner::wigner6j( TwoSL, TwoSR, 1, TwoSRdown, TwoSLdown, 1 );
1898 
1899  for (int Irrep=0; Irrep<(denBK->getNumberOfIrreps()); Irrep++){
1900 
1901  int IrrepTimesMid = Irreps::directProd(Irrep,IprodMID);
1902 
1903  bool leftOK = false;
1904  for (int l_alpha=0; l_alpha<theindex; l_alpha++){
1905  if (denBK->gIrrep(l_alpha) == Irrep){ leftOK = true; }
1906  }
1907  bool rightOK = false;
1908  for (int l_beta=theindex+2; l_beta<Prob->gL(); l_beta++){
1909  if (denBK->gIrrep(l_beta) == IrrepTimesMid){ rightOK = true; }
1910  }
1911 
1912  if ((leftOK) && (rightOK)){
1913 
1914  for (int l_alpha=0; l_alpha<theindex; l_alpha++){
1915  if (denBK->gIrrep(l_alpha) == Irrep){
1916 
1917  int ILdown = Irreps::directProd(IL, Irrep);
1918  int IRdown = Irreps::directProd(IR, IrrepTimesMid);
1919 
1920  int dimLdown = denBK->gCurrentDim(theindex, NL+1, TwoSLdown, ILdown);
1921  int dimRdown = denBK->gCurrentDim(theindex+2, NR+1, TwoSRdown, IRdown);
1922 
1923  if ((dimLdown>0) && (dimRdown>0)){
1924 
1925  int size = dimRup * dimRdown;
1926  for (int cnt=0; cnt<size; cnt++){ temp[cnt] = 0.0; }
1927 
1928  for (int l_beta=theindex+2; l_beta<Prob->gL(); l_beta++){
1929  if (denBK->gIrrep(l_beta) == IrrepTimesMid){
1930  double prefac = factor1 * Prob->gMxElement(l_alpha, theindex+1, theindex, l_beta)
1931  + factor2 * Prob->gMxElement(l_alpha, theindex+1, l_beta, theindex);
1932  double * LblockR = Lright[l_beta-theindex-2]->gStorage(NR,TwoSR,IR,NR+1,TwoSRdown,IRdown);
1933  daxpy_(&size,&prefac,LblockR,&inc,temp,&inc);
1934  }
1935  }
1936 
1937  double alpha = 1.0;
1938  double beta = 0.0; //set
1939  int memSkappa = denS->gKappa(NL+1, TwoSLdown, ILdown, 1, 2, 1, NR+1, TwoSRdown, IRdown);
1940 
1941  dgemm_(&notrans,&trans,&dimLdown,&dimRup,&dimRdown,&alpha,memS+denS->gKappa2index(memSkappa),&dimLdown,temp,&dimRup,&beta,temp2,&dimLdown);
1942 
1943  beta = 1.0; //add
1944  double * LblockL = Lleft[theindex-1-l_alpha]->gStorage(NL,TwoSL,IL,NL+1,TwoSLdown,ILdown);
1945  dgemm_(&notrans,&notrans,&dimLup,&dimRup,&dimLdown,&alpha,LblockL,&dimLup,temp2,&dimLdown,&beta,memHeff+denS->gKappa2index(ikappa),&dimLup);
1946  }
1947  }
1948  }
1949  }
1950  }
1951  }
1952  }
1953  }
1954  }
1955 
1956 }
1957 
1958 
1959 
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::SyBookkeeper::getNumberOfIrreps
int getNumberOfIrreps() const
Get the total number of irreps.
Definition: SyBookkeeper.cpp:153
CheMPS2::Heff::phase
static int phase(const int TwoTimesPower)
Phase function.
Definition: Heff.h:75
CheMPS2::SyBookkeeper::gIrrep
int gIrrep(const int orbital) const
Get an orbital irrep.
Definition: SyBookkeeper.cpp:145