1 /*
2    CheMPS2: a spin-adapted implementation of DMRG for ab initio quantum chemistry
3    Copyright (C) 2013-2018 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 
addDiagram3Aand3D(const int ikappa,double * memS,double * memHeff,const Sobject * denS,TensorQ * Qleft,TensorL ** Lleft,double * temp) const28 void CheMPS2::Heff::addDiagram3Aand3D(const int ikappa, double * memS, double * memHeff, const Sobject * denS, TensorQ * Qleft, TensorL ** Lleft, double * temp) const{
29 
30    int NL = denS->gNL(ikappa);
31    int TwoSL = denS->gTwoSL(ikappa);
32    int IL = denS->gIL(ikappa);
33    int N1 = denS->gN1(ikappa);
34    int N2 = denS->gN2(ikappa);
35    int TwoJ = denS->gTwoJ(ikappa);
36    int NR = denS->gNR(ikappa);
37    int TwoSR = denS->gTwoSR(ikappa);
38    int IR = denS->gIR(ikappa);
39 
40    int theindex = denS->gIndex();
41    int ILdown = Irreps::directProd(IL,denBK->gIrrep(theindex));
42    int TwoS2 = (N2==1)?1:0;
43 
44    int dimR = denBK->gCurrentDim(theindex+2,NR,TwoSR,IR);
45    int dimLup = denBK->gCurrentDim(theindex,NL,TwoSL,IL);
46 
47    if (N1==2){ //3A1A and 3D1
48       for (int TwoSLdown=TwoSL-1; TwoSLdown<=TwoSL+1; TwoSLdown+=2){
49 
50          int dimLdown = denBK->gCurrentDim(theindex, NL+1, TwoSLdown, ILdown);
51          if (dimLdown>0){
52 
53             int TwoJstart = ((TwoSR!=TwoSLdown) || (TwoS2==0)) ? 1 + TwoS2 : 0;
54             for (int TwoJdown=TwoJstart; TwoJdown<=1+TwoS2; TwoJdown+=2){
55                if (abs(TwoSLdown-TwoSR)<=TwoJdown){
56 
57                   int memSkappa = denS->gKappa(NL+1,TwoSLdown,ILdown,1,N2,TwoJdown,NR,TwoSR,IR);
58                   if (memSkappa!=-1){
59                      int fase = phase(TwoSL+TwoSR+2+TwoS2);
60                      double factor = sqrt((TwoJdown+1)*(TwoSLdown+1.0))*fase*Wigner::wigner6j(TwoJdown,TwoS2,1,TwoSL,TwoSLdown,TwoSR);
61                      double beta = 1.0; //add
62                      char notr = 'N';
63 
64                      double * BlockQ = Qleft->gStorage(NL,TwoSL,IL,NL+1,TwoSLdown,ILdown);
65                      int inc = 1;
66                      int size = dimLup * dimLdown;
67                      dcopy_(&size, BlockQ, &inc, temp, &inc);
68 
69                      for (int l_index=0; l_index<theindex; l_index++){
70                         if (denBK->gIrrep(l_index) == denBK->gIrrep(theindex)){
71                            double alpha = Prob->gMxElement(l_index,theindex,theindex,theindex);
72                            double * BlockL = Lleft[theindex-1-l_index]->gStorage(NL,TwoSL,IL,NL+1,TwoSLdown,ILdown);
73                            daxpy_(&size, &alpha, BlockL, &inc, temp, &inc);
74                         }
75                      }
76 
77                      dgemm_(&notr,&notr,&dimLup,&dimR,&dimLdown,&factor,temp,&dimLup,memS+denS->gKappa2index(memSkappa),&dimLdown,&beta,memHeff+denS->gKappa2index(ikappa),&dimLup);
78                   }
79                }
80             }
81          }
82       }
83    }
84 
85    if (N1==1){ //3A1B
86       for (int TwoSLdown=TwoSL-1; TwoSLdown<=TwoSL+1; TwoSLdown+=2){
87          if ((abs(TwoSLdown-TwoSR)<=TwoS2) && (TwoSLdown>=0)){
88             int dimLdown = denBK->gCurrentDim(theindex, NL+1, TwoSLdown, ILdown);
89             int memSkappa = denS->gKappa(NL+1,TwoSLdown,ILdown,0,N2,TwoS2,NR,TwoSR,IR);
90             if (memSkappa!=-1){
91                int fase = phase(TwoSL+TwoSR+1+TwoS2);
92                double factor = sqrt((TwoSLdown+1)*(TwoJ+1.0))*fase*Wigner::wigner6j(TwoS2,TwoJ,1,TwoSL,TwoSLdown,TwoSR);
93                double beta = 1.0;
94                char notr = 'N';
95                double * BlockQ = Qleft->gStorage(NL,TwoSL,IL,NL+1,TwoSLdown,ILdown);
96                dgemm_(&notr,&notr,&dimLup,&dimR,&dimLdown,&factor,BlockQ,&dimLup,memS+denS->gKappa2index(memSkappa),&dimLdown,&beta,memHeff+denS->gKappa2index(ikappa),&dimLup);
97             }
98          }
99       }
100    }
101 
102    if (N1==0){ //3A2A
103       for (int TwoSLdown=TwoSL-1;TwoSLdown<=TwoSL+1;TwoSLdown+=2){
104 
105          int dimLdown = denBK->gCurrentDim(theindex, NL-1, TwoSLdown, ILdown);
106          if (dimLdown>0){
107 
108             int TwoJstart = ((TwoSR!=TwoSLdown) || (TwoS2==0)) ? 1 + TwoS2 : 0;
109             for (int TwoJdown=TwoJstart; TwoJdown<=1+TwoS2; TwoJdown+=2){
110                if (abs(TwoSLdown-TwoSR)<=TwoJdown){
111 
112                   int memSkappa = denS->gKappa(NL-1,TwoSLdown,ILdown,1,N2,TwoJdown,NR,TwoSR,IR);
113                   if (memSkappa!=-1){
114                      int fase = phase(TwoSLdown+TwoSR+1+TwoS2);
115                      double factor = fase*sqrt((TwoSL+1)*(TwoJdown+1.0))*Wigner::wigner6j(TwoJdown,TwoS2,1,TwoSL,TwoSLdown,TwoSR);
116                      double beta = 1.0;
117                      char notr = 'N';
118                      char trans = 'T';
119                      double * BlockQ = Qleft->gStorage(NL-1,TwoSLdown,ILdown,NL,TwoSL,IL);
120                      dgemm_(&trans,&notr,&dimLup,&dimR,&dimLdown,&factor,BlockQ,&dimLdown,memS+denS->gKappa2index(memSkappa),&dimLdown,&beta,memHeff+denS->gKappa2index(ikappa),&dimLup);
121                   }
122                }
123             }
124          }
125       }
126    }
127 
128    if (N1==1){ //3A2B ans 3D2
129       for (int TwoSLdown=TwoSL-1;TwoSLdown<=TwoSL+1;TwoSLdown+=2){
130          if ((abs(TwoSLdown-TwoSR)<=TwoS2) && (TwoSLdown>=0)){
131             int dimLdown = denBK->gCurrentDim(theindex, NL-1, TwoSLdown, ILdown);
132             int memSkappa = denS->gKappa(NL-1,TwoSLdown,ILdown,2,N2,TwoS2,NR,TwoSR,IR);
133             if (memSkappa!=-1){
134                int fase = phase(TwoSLdown+TwoSR+2+TwoS2);
135                double factor = fase*sqrt((TwoSL+1)*(TwoJ+1.0))*Wigner::wigner6j(TwoS2,TwoJ,1,TwoSL,TwoSLdown,TwoSR);
136                double beta = 1.0;
137                char notr = 'N';
138                char trans = 'T';
139 
140                double * BlockQ = Qleft->gStorage(NL-1,TwoSLdown,ILdown,NL,TwoSL,IL);
141                int inc = 1;
142                int size = dimLup * dimLdown;
143                dcopy_(&size, BlockQ, &inc, temp, &inc);
144 
145                for (int l_index=0; l_index<theindex; l_index++){
146                   if (denBK->gIrrep(l_index) == denBK->gIrrep(theindex)){
147                      double alpha = Prob->gMxElement(l_index,theindex,theindex,theindex);
148                      double * BlockL = Lleft[theindex-1-l_index]->gStorage(NL-1,TwoSLdown,ILdown,NL,TwoSL,IL);
149                      daxpy_(&size, &alpha, BlockL, &inc, temp, &inc);
150                   }
151                }
152 
153                dgemm_(&trans,&notr,&dimLup,&dimR,&dimLdown,&factor,temp,&dimLdown,memS+denS->gKappa2index(memSkappa),&dimLdown,&beta,memHeff+denS->gKappa2index(ikappa),&dimLup);
154             }
155          }
156       }
157    }
158 
159 }
160 
addDiagram3Band3I(const int ikappa,double * memS,double * memHeff,const Sobject * denS,TensorQ * Qleft,TensorL ** Lleft,double * temp) const161 void CheMPS2::Heff::addDiagram3Band3I(const int ikappa, double * memS, double * memHeff, const Sobject * denS, TensorQ * Qleft, TensorL ** Lleft, double * temp) const{
162 
163    int NL = denS->gNL(ikappa);
164    int TwoSL = denS->gTwoSL(ikappa);
165    int IL = denS->gIL(ikappa);
166    int N1 = denS->gN1(ikappa);
167    int N2 = denS->gN2(ikappa);
168    int TwoJ = denS->gTwoJ(ikappa);
169    int NR = denS->gNR(ikappa);
170    int TwoSR = denS->gTwoSR(ikappa);
171    int IR = denS->gIR(ikappa);
172 
173    int theindex = denS->gIndex();
174    int ILdown = Irreps::directProd(IL,denBK->gIrrep(theindex+1));
175    int TwoS1 = (N1==1)?1:0;
176 
177    int dimR = denBK->gCurrentDim(theindex+2,NR,TwoSR,IR);
178    int dimLup = denBK->gCurrentDim(theindex,NL,TwoSL,IL);
179 
180    if (N2==2){ //3B1A and 3I2
181       for (int TwoSLdown=TwoSL-1; TwoSLdown<=TwoSL+1; TwoSLdown+=2){
182 
183          int dimLdown = denBK->gCurrentDim(theindex, NL+1, TwoSLdown, ILdown);
184          if (dimLdown>0){
185 
186             int TwoJstart = ((TwoSR!=TwoSLdown) || (TwoS1==0)) ? 1 + TwoS1 : 0;
187             for (int TwoJdown=TwoJstart; TwoJdown<=1+TwoS1; TwoJdown+=2){
188                if (abs(TwoSLdown-TwoSR)<=TwoJdown){
189 
190                   int memSkappa = denS->gKappa(NL+1,TwoSLdown,ILdown,N1,1,TwoJdown,NR,TwoSR,IR);
191                   if (memSkappa!=-1){
192                      int fase = phase(TwoSL+TwoSR+3-TwoJdown);
193                      double factor = sqrt((TwoJdown+1)*(TwoSLdown+1.0))*fase*Wigner::wigner6j(TwoJdown,TwoS1,1,TwoSL,TwoSLdown,TwoSR);
194                      double beta = 1.0; //add
195                      char notr = 'N';
196 
197                      double * BlockQ = Qleft->gStorage(NL,TwoSL,IL,NL+1,TwoSLdown,ILdown);
198                      int inc = 1;
199                      int size = dimLup * dimLdown;
200                      dcopy_(&size, BlockQ, &inc, temp, &inc);
201 
202                      for (int l_index=0; l_index<theindex; l_index++){
203                         if (denBK->gIrrep(l_index) == denBK->gIrrep(theindex+1)){
204                            double alpha = Prob->gMxElement(l_index,theindex+1,theindex+1,theindex+1);
205                            double * BlockL = Lleft[theindex-1-l_index]->gStorage(NL,TwoSL,IL,NL+1,TwoSLdown,ILdown);
206                            daxpy_(&size, &alpha, BlockL, &inc, temp, &inc);
207                         }
208                      }
209 
210                      dgemm_(&notr,&notr,&dimLup,&dimR,&dimLdown,&factor,temp,&dimLup,memS+denS->gKappa2index(memSkappa),&dimLdown,&beta,memHeff+denS->gKappa2index(ikappa),&dimLup);
211                   }
212                }
213             }
214          }
215       }
216    }
217 
218    if (N2==1){ //3B1B
219       for (int TwoSLdown=TwoSL-1; TwoSLdown<=TwoSL+1; TwoSLdown+=2){
220          if ((abs(TwoSLdown-TwoSR)<=TwoS1) && (TwoSLdown>=0)){
221             int dimLdown = denBK->gCurrentDim(theindex, NL+1, TwoSLdown, ILdown);
222             int memSkappa = denS->gKappa(NL+1,TwoSLdown,ILdown,N1,0,TwoS1,NR,TwoSR,IR);
223             if (memSkappa!=-1){
224                int fase = phase(TwoSL+TwoSR+2-TwoJ);
225                double factor = sqrt((TwoSLdown+1)*(TwoJ+1.0))*fase*Wigner::wigner6j(TwoS1,TwoJ,1,TwoSL,TwoSLdown,TwoSR);
226                double beta = 1.0;
227                char notr = 'N';
228                double * BlockQ = Qleft->gStorage(NL,TwoSL,IL,NL+1,TwoSLdown,ILdown);
229                dgemm_(&notr,&notr,&dimLup,&dimR,&dimLdown,&factor,BlockQ,&dimLup,memS+denS->gKappa2index(memSkappa),&dimLdown,&beta,memHeff+denS->gKappa2index(ikappa),&dimLup);
230             }
231          }
232       }
233    }
234 
235    if (N2==0){ //3B2A
236       for (int TwoSLdown=TwoSL-1;TwoSLdown<=TwoSL+1;TwoSLdown+=2){
237 
238          int dimLdown = denBK->gCurrentDim(theindex, NL-1, TwoSLdown, ILdown);
239          if (dimLdown>0){
240 
241             int TwoJstart = ((TwoSR!=TwoSLdown) || (TwoS1==0)) ? 1 + TwoS1 : 0;
242             for (int TwoJdown=TwoJstart; TwoJdown<=1+TwoS1; TwoJdown+=2){
243                if (abs(TwoSLdown-TwoSR)<=TwoJdown){
244 
245                   int memSkappa = denS->gKappa(NL-1,TwoSLdown,ILdown,N1,1,TwoJdown,NR,TwoSR,IR);
246                   if (memSkappa!=-1){
247                      int fase = phase(TwoSLdown+TwoSR+2-TwoJdown);
248                      double factor = fase*sqrt((TwoSL+1)*(TwoJdown+1.0))*Wigner::wigner6j(TwoJdown,TwoS1,1,TwoSL,TwoSLdown,TwoSR);
249                      double beta = 1.0;
250                      char notr = 'N';
251                      char trans = 'T';
252                      double * BlockQ = Qleft->gStorage(NL-1,TwoSLdown,ILdown,NL,TwoSL,IL);
253                      dgemm_(&trans,&notr,&dimLup,&dimR,&dimLdown,&factor,BlockQ,&dimLdown,memS+denS->gKappa2index(memSkappa),&dimLdown,&beta,memHeff+denS->gKappa2index(ikappa),&dimLup);
254                   }
255                }
256             }
257          }
258       }
259    }
260 
261    if (N2==1){ //3B2B and 3I1
262       for (int TwoSLdown=TwoSL-1;TwoSLdown<=TwoSL+1;TwoSLdown+=2){
263          if ((abs(TwoSLdown-TwoSR)<=TwoS1) && (TwoSLdown>=0)){
264             int dimLdown = denBK->gCurrentDim(theindex, NL-1, TwoSLdown, ILdown);
265             int memSkappa = denS->gKappa(NL-1,TwoSLdown,ILdown,N1,2,TwoS1,NR,TwoSR,IR);
266             if (memSkappa!=-1){
267                int fase = phase(TwoSLdown+TwoSR+3-TwoJ);
268                double factor = fase*sqrt((TwoSL+1)*(TwoJ+1.0))*Wigner::wigner6j(TwoS1,TwoJ,1,TwoSL,TwoSLdown,TwoSR);
269                double beta = 1.0;
270                char notr = 'N';
271                char trans = 'T';
272 
273                double * BlockQ = Qleft->gStorage(NL-1,TwoSLdown,ILdown,NL,TwoSL,IL);
274                int inc = 1;
275                int size = dimLup * dimLdown;
276                dcopy_(&size, BlockQ, &inc, temp, &inc);
277 
278                for (int l_index=0; l_index<theindex; l_index++){
279                   if (denBK->gIrrep(l_index) == denBK->gIrrep(theindex+1)){
280                      double alpha = Prob->gMxElement(l_index,theindex+1,theindex+1,theindex+1);
281                      double * BlockL = Lleft[theindex-1-l_index]->gStorage(NL-1,TwoSLdown,ILdown,NL,TwoSL,IL);
282                      daxpy_(&size, &alpha, BlockL, &inc, temp, &inc);
283                   }
284                }
285 
286                dgemm_(&trans,&notr,&dimLup,&dimR,&dimLdown,&factor,temp,&dimLdown,memS+denS->gKappa2index(memSkappa),&dimLdown,&beta,memHeff+denS->gKappa2index(ikappa),&dimLup);
287             }
288          }
289       }
290    }
291 
292 }
293 
addDiagram3C(const int ikappa,double * memS,double * memHeff,const Sobject * denS,TensorQ ** Qleft,TensorL ** Lright,double * temp) const294 void CheMPS2::Heff::addDiagram3C(const int ikappa, double * memS, double * memHeff, const Sobject * denS, TensorQ ** Qleft, TensorL ** Lright, double * temp) const{
295 
296    #ifdef CHEMPS2_MPI_COMPILATION
297    const int MPIRANK = MPIchemps2::mpi_rank();
298    #endif
299 
300    int NL = denS->gNL(ikappa);
301    int TwoSL = denS->gTwoSL(ikappa);
302    int IL = denS->gIL(ikappa);
303    int N1 = denS->gN1(ikappa);
304    int N2 = denS->gN2(ikappa);
305    int TwoJ = denS->gTwoJ(ikappa);
306    int NR = denS->gNR(ikappa);
307    int TwoSR = denS->gTwoSR(ikappa);
308    int IR = denS->gIR(ikappa);
309 
310    int theindex = denS->gIndex();
311 
312    int dimRup = denBK->gCurrentDim(theindex+2,NR,TwoSR,IR);
313    int dimLup = denBK->gCurrentDim(theindex,  NL,TwoSL,IL);
314 
315    //First do 3C1
316    for (int TwoSLdown=TwoSL-1; TwoSLdown<=TwoSL+1; TwoSLdown+=2){
317       for (int TwoSRdown=TwoSR-1; TwoSRdown<=TwoSR+1; TwoSRdown+=2){
318          if ((abs(TwoSLdown-TwoSRdown)<=TwoJ) && (TwoSLdown>=0) && (TwoSRdown>=0)){
319 
320             int fase = phase(TwoSLdown+TwoSR+TwoJ+1 + ((N1==1)?2:0) + ((N2==1)?2:0) );
321             const double factor = fase * sqrt((TwoSLdown+1)*(TwoSRdown+1.0)) * Wigner::wigner6j(TwoSL,TwoSR,TwoJ,TwoSRdown,TwoSLdown,1);
322 
323             for (int l_index=theindex+2; l_index<Prob->gL(); l_index++){
324 
325                #ifdef CHEMPS2_MPI_COMPILATION
326                if ( MPIchemps2::owner_q( Prob->gL(), l_index ) == MPIRANK )
327                #endif
328                {
329                   int ILdown = Irreps::directProd(IL,denBK->gIrrep(l_index));
330                   int IRdown = Irreps::directProd(IR,denBK->gIrrep(l_index));
331                   int memSkappa = denS->gKappa(NL+1, TwoSLdown, ILdown, N1, N2, TwoJ, NR+1, TwoSRdown, IRdown);
332                   if (memSkappa!=-1){
333                      int dimRdown = denBK->gCurrentDim(theindex+2, NR+1, TwoSRdown, IRdown);
334                      int dimLdown = denBK->gCurrentDim(theindex,   NL+1, TwoSLdown, ILdown);
335 
336                      double * Qblock = Qleft[ l_index-theindex  ]->gStorage(NL,TwoSL,IL,NL+1,TwoSLdown,ILdown);
337                      double * Lblock = Lright[l_index-theindex-2]->gStorage(NR,TwoSR,IR,NR+1,TwoSRdown,IRdown);
338 
339                      char trans = 'T';
340                      char notra = 'N';
341                      double beta = 0.0; //set
342                      double alpha = factor;
343                      dgemm_(&notra,&notra,&dimLup,&dimRdown,&dimLdown,&alpha,Qblock,&dimLup,memS+denS->gKappa2index(memSkappa),&dimLdown,&beta,temp,&dimLup);
344 
345                      beta = 1.0; //add
346                      alpha = 1.0;
347                      dgemm_(&notra,&trans,&dimLup,&dimRup,&dimRdown,&alpha,temp,&dimLup,Lblock,&dimRup,&beta,memHeff+denS->gKappa2index(ikappa),&dimLup);
348 
349                   }
350                }
351             }
352          }
353       }
354    }
355 
356    //Then do 3C2
357    for (int TwoSLdown=TwoSL-1; TwoSLdown<=TwoSL+1; TwoSLdown+=2){
358       for (int TwoSRdown=TwoSR-1; TwoSRdown<=TwoSR+1; TwoSRdown+=2){
359          if ((abs(TwoSLdown-TwoSRdown)<=TwoJ) && (TwoSLdown>=0) && (TwoSRdown>=0)){
360 
361             int fase = phase(TwoSL+TwoSRdown+TwoJ+1 + ((N1==1)?2:0) + ((N2==1)?2:0) );
362             const double factor = fase * sqrt((TwoSL+1)*(TwoSR+1.0)) * Wigner::wigner6j(TwoSL,TwoSR,TwoJ,TwoSRdown,TwoSLdown,1);
363 
364             for (int l_index=theindex+2; l_index<Prob->gL(); l_index++){
365 
366                #ifdef CHEMPS2_MPI_COMPILATION
367                if ( MPIchemps2::owner_q( Prob->gL(), l_index ) == MPIRANK )
368                #endif
369                {
370                   int ILdown = Irreps::directProd(IL,denBK->gIrrep(l_index));
371                   int IRdown = Irreps::directProd(IR,denBK->gIrrep(l_index));
372                   int memSkappa = denS->gKappa(NL-1, TwoSLdown, ILdown, N1, N2, TwoJ, NR-1, TwoSRdown, IRdown);
373                   if (memSkappa!=-1){
374                      int dimRdown = denBK->gCurrentDim(theindex+2, NR-1, TwoSRdown, IRdown);
375                      int dimLdown = denBK->gCurrentDim(theindex,   NL-1, TwoSLdown, ILdown);
376 
377                      double * Qblock = Qleft[ l_index-theindex  ]->gStorage(NL-1,TwoSLdown,ILdown,NL,TwoSL,IL);
378                      double * Lblock = Lright[l_index-theindex-2]->gStorage(NR-1,TwoSRdown,IRdown,NR,TwoSR,IR);
379 
380                      char trans = 'T';
381                      char notra = 'N';
382                      double beta = 0.0; //set
383                      double alpha = factor;
384                      dgemm_(&trans,&notra,&dimLup,&dimRdown,&dimLdown,&alpha,Qblock,&dimLdown,memS+denS->gKappa2index(memSkappa),&dimLdown,&beta,temp,&dimLup);
385 
386                      beta = 1.0; //add
387                      alpha = 1.0;
388                      dgemm_(&notra,&notra,&dimLup,&dimRup,&dimRdown,&alpha,temp,&dimLup,Lblock,&dimRdown,&beta,memHeff+denS->gKappa2index(ikappa),&dimLup);
389 
390                   }
391                }
392             }
393          }
394       }
395    }
396 
397 }
398 
addDiagram3Eand3H(const int ikappa,double * memS,double * memHeff,const Sobject * denS) const399 void CheMPS2::Heff::addDiagram3Eand3H(const int ikappa, double * memS, double * memHeff, const Sobject * denS) const{ //TwoJ = TwoJdown
400 
401    int theindex = denS->gIndex();
402 
403    if (denBK->gIrrep(theindex) != denBK->gIrrep(theindex+1)){ return; }
404 
405    int NL = denS->gNL(ikappa);
406    int TwoSL = denS->gTwoSL(ikappa);
407    int IL = denS->gIL(ikappa);
408    int N1 = denS->gN1(ikappa);
409    int N2 = denS->gN2(ikappa);
410    int TwoJ = denS->gTwoJ(ikappa);
411    int NR = denS->gNR(ikappa);
412    int TwoSR = denS->gTwoSR(ikappa);
413    int IR = denS->gIR(ikappa);
414 
415    int size = ( denBK->gCurrentDim(theindex,NL,TwoSL,IL) ) * ( denBK->gCurrentDim(theindex+2,NR,TwoSR,IR) );
416    int inc = 1;
417 
418    if ((N1==2) && (N2==0)){ //3E1A
419 
420       int memSkappa = denS->gKappa(NL,TwoSL,IL,1,1,0,NR,TwoSR,IR);
421       if (memSkappa!=-1){
422          double alpha = sqrt(2.0) * Prob->gMxElement(theindex,theindex,theindex,theindex+1);
423          daxpy_(&size,&alpha,memS+denS->gKappa2index(memSkappa),&inc,memHeff+denS->gKappa2index(ikappa),&inc);
424       }
425 
426    }
427 
428    if ((N1==2) && (N2==1)){ //3E1B and 3H1B
429 
430       int memSkappa = denS->gKappa(NL,TwoSL,IL,1,2,1,NR,TwoSR,IR);
431       if (memSkappa!=-1){
432          double alpha = - ( Prob->gMxElement(theindex,theindex,theindex,theindex+1) + Prob->gMxElement(theindex,theindex+1,theindex+1,theindex+1) );
433          daxpy_(&size,&alpha,memS+denS->gKappa2index(memSkappa),&inc,memHeff+denS->gKappa2index(ikappa),&inc);
434       }
435 
436    }
437 
438    if ((N1==1) && (N2==1) && (TwoJ==0)){ //3E2A and 3H1A
439 
440       int memSkappa = denS->gKappa(NL,TwoSL,IL,2,0,0,NR,TwoSR,IR);
441       if (memSkappa!=-1){
442          double alpha = sqrt(2.0) * Prob->gMxElement(theindex,theindex,theindex,theindex+1);
443          daxpy_(&size,&alpha,memS+denS->gKappa2index(memSkappa),&inc,memHeff+denS->gKappa2index(ikappa),&inc);
444       }
445 
446       memSkappa = denS->gKappa(NL,TwoSL,IL,0,2,0,NR,TwoSR,IR);
447       if (memSkappa!=-1){
448          double alpha = sqrt(2.0) * Prob->gMxElement(theindex,theindex+1,theindex+1,theindex+1);
449          daxpy_(&size,&alpha,memS+denS->gKappa2index(memSkappa),&inc,memHeff+denS->gKappa2index(ikappa),&inc);
450       }
451 
452    }
453 
454    if ((N1==1) && (N2==2)){ //3E2B and 3H2B
455 
456       int memSkappa = denS->gKappa(NL,TwoSL,IL,2,1,1,NR,TwoSR,IR);
457       if (memSkappa!=-1){
458          double alpha = - ( Prob->gMxElement(theindex,theindex,theindex,theindex+1) + Prob->gMxElement(theindex,theindex+1,theindex+1,theindex+1) );
459          daxpy_(&size,&alpha,memS+denS->gKappa2index(memSkappa),&inc,memHeff+denS->gKappa2index(ikappa),&inc);
460       }
461 
462    }
463 
464    if ((N1==0) && (N2==2)){ //3H2A
465 
466       int memSkappa = denS->gKappa(NL,TwoSL,IL,1,1,0,NR,TwoSR,IR);
467       if (memSkappa!=-1){
468          double alpha = sqrt(2.0) * Prob->gMxElement(theindex,theindex+1,theindex+1,theindex+1);
469          daxpy_(&size,&alpha,memS+denS->gKappa2index(memSkappa),&inc,memHeff+denS->gKappa2index(ikappa),&inc);
470       }
471 
472    }
473 
474 }
475 
addDiagram3Kand3F(const int ikappa,double * memS,double * memHeff,const Sobject * denS,TensorQ * Qright,TensorL ** Lright,double * temp) const476 void CheMPS2::Heff::addDiagram3Kand3F(const int ikappa, double * memS, double * memHeff, const Sobject * denS, TensorQ * Qright, TensorL ** Lright, double * temp) const{
477 
478    int NL = denS->gNL(ikappa);
479    int TwoSL = denS->gTwoSL(ikappa);
480    int IL = denS->gIL(ikappa);
481    int N1 = denS->gN1(ikappa);
482    int N2 = denS->gN2(ikappa);
483    int TwoJ = denS->gTwoJ(ikappa);
484    int NR = denS->gNR(ikappa);
485    int TwoSR = denS->gTwoSR(ikappa);
486    int IR = denS->gIR(ikappa);
487 
488    int theindex = denS->gIndex();
489    int IRdown = Irreps::directProd(IR,denBK->gIrrep(theindex));
490    int TwoS2 = (N2==1)?1:0;
491 
492    int dimL   = denBK->gCurrentDim(theindex,  NL,TwoSL,IL);
493    int dimRup = denBK->gCurrentDim(theindex+2,NR,TwoSR,IR);
494 
495    if (N1==1){ //3K1A
496       for (int TwoSRdown=TwoSR-1; TwoSRdown<=TwoSR+1; TwoSRdown+=2){
497          if ((abs(TwoSL-TwoSRdown)<=TwoS2) && (TwoSRdown>=0)){
498             int dimRdown = denBK->gCurrentDim(theindex+2, NR-1, TwoSRdown, IRdown);
499             int memSkappa = denS->gKappa(NL,TwoSL,IL,0,N2,TwoS2,NR-1,TwoSRdown,IRdown);
500             if (memSkappa!=-1){
501                int fase = phase(TwoSL+TwoSR+TwoJ+2*TwoS2);
502                double factor = sqrt((TwoJ+1)*(TwoSR+1.0)) * fase * Wigner::wigner6j(TwoS2,TwoJ,1,TwoSR,TwoSRdown,TwoSL);
503                double beta = 1.0; //add
504                char notr = 'N';
505                double * BlockQ = Qright->gStorage(NR-1,TwoSRdown,IRdown,NR,TwoSR,IR);
506                dgemm_(&notr,&notr,&dimL,&dimRup,&dimRdown,&factor,memS+denS->gKappa2index(memSkappa),&dimL,BlockQ,&dimRdown,&beta,memHeff+denS->gKappa2index(ikappa),&dimL);
507             }
508          }
509       }
510    }
511 
512    if (N1==2){ //3K1B and 3F1
513       for (int TwoSRdown=TwoSR-1; TwoSRdown<=TwoSR+1; TwoSRdown+=2){
514 
515          int dimRdown = denBK->gCurrentDim(theindex+2, NR-1, TwoSRdown, IRdown);
516          if (dimRdown>0){
517 
518             int TwoJstart = ((TwoSRdown!=TwoSL) || (TwoS2==0)) ? 1 + TwoS2 : 0;
519             for (int TwoJdown=TwoJstart; TwoJdown<=1+TwoS2; TwoJdown+=2){
520                if (abs(TwoSL-TwoSRdown)<=TwoJdown){
521 
522                   int memSkappa = denS->gKappa(NL,TwoSL,IL,1,N2,TwoJdown,NR-1,TwoSRdown,IRdown);
523                   if (memSkappa!=-1){
524                      int fase = phase(TwoSL+TwoSR+TwoJdown+1+2*TwoS2);
525                      double factor = sqrt((TwoJdown+1)*(TwoSR+1.0)) * fase * Wigner::wigner6j( TwoJdown, TwoS2, 1, TwoSR, TwoSRdown, TwoSL );
526                      double beta = 1.0; //add
527                      char notr = 'N';
528 
529                      double * BlockQ = Qright->gStorage(NR-1,TwoSRdown,IRdown,NR,TwoSR,IR);
530                      int inc = 1;
531                      int size = dimRup * dimRdown;
532                      dcopy_(&size,BlockQ,&inc,temp,&inc);
533 
534                      for (int l_index=theindex+2; l_index<Prob->gL(); l_index++){
535                         if (denBK->gIrrep(l_index) == denBK->gIrrep(theindex)){
536                            double alpha = Prob->gMxElement(theindex,theindex,theindex,l_index);
537                            double * BlockL = Lright[l_index-theindex-2]->gStorage(NR-1,TwoSRdown,IRdown,NR,TwoSR,IR);
538                            daxpy_(&size, &alpha, BlockL, &inc, temp, &inc);
539                         }
540                      }
541 
542                      dgemm_(&notr,&notr,&dimL,&dimRup,&dimRdown,&factor,memS+denS->gKappa2index(memSkappa),&dimL,temp,&dimRdown,&beta,memHeff+denS->gKappa2index(ikappa),&dimL);
543                   }
544                }
545             }
546          }
547       }
548    }
549 
550    if (N1==0){ //3K2A
551       for (int TwoSRdown=TwoSR-1; TwoSRdown<=TwoSR+1; TwoSRdown+=2){
552 
553          int dimRdown = denBK->gCurrentDim(theindex+2, NR+1, TwoSRdown, IRdown);
554          if (dimRdown>0){
555 
556             int TwoJstart = ((TwoSRdown!=TwoSL) || (TwoS2==0)) ? 1 + TwoS2 : 0;
557             for (int TwoJdown=TwoJstart; TwoJdown<=1+TwoS2; TwoJdown+=2){
558                if (abs(TwoSL-TwoSRdown)<=TwoJdown){
559 
560                   int memSkappa = denS->gKappa(NL,TwoSL,IL,1,N2,TwoJdown,NR+1,TwoSRdown,IRdown);
561                   if (memSkappa!=-1){
562                      int fase = phase(TwoSL+TwoSRdown+TwoJdown+2*TwoS2);
563                      double factor = sqrt((TwoJdown+1)*(TwoSRdown+1.0)) * fase * Wigner::wigner6j(TwoJdown,TwoS2,1,TwoSR,TwoSRdown,TwoSL);
564                      double beta = 1.0; //add
565                      char notr = 'N';
566                      char tran = 'T';
567                      double * BlockQ = Qright->gStorage(NR,TwoSR,IR,NR+1,TwoSRdown,IRdown);
568                      dgemm_(&notr,&tran,&dimL,&dimRup,&dimRdown,&factor,memS+denS->gKappa2index(memSkappa),&dimL,BlockQ,&dimRup,&beta,memHeff+denS->gKappa2index(ikappa),&dimL);
569                   }
570                }
571             }
572          }
573       }
574    }
575 
576    if (N1==1){ //3K2B and 3F2
577       for (int TwoSRdown=TwoSR-1; TwoSRdown<=TwoSR+1; TwoSRdown+=2){
578          if ((abs(TwoSL-TwoSRdown)<=TwoS2) && (TwoSRdown>=0)){
579             int dimRdown = denBK->gCurrentDim(theindex+2, NR+1, TwoSRdown, IRdown);
580             int memSkappa = denS->gKappa(NL,TwoSL,IL,2,N2,TwoS2,NR+1,TwoSRdown,IRdown);
581             if (memSkappa!=-1){
582                int fase = phase(TwoSL+TwoSRdown+TwoJ+1+2*TwoS2);
583                double factor = sqrt((TwoJ+1)*(TwoSRdown+1.0)) * fase * Wigner::wigner6j(TwoS2,TwoJ,1,TwoSR,TwoSRdown,TwoSL);
584                double beta = 1.0; //add
585                char notr = 'N';
586                char tran = 'T';
587 
588                double * BlockQ = Qright->gStorage(NR,TwoSR,IR,NR+1,TwoSRdown,IRdown);
589                int inc = 1;
590                int size = dimRup * dimRdown;
591                dcopy_(&size,BlockQ,&inc,temp,&inc);
592 
593                for (int l_index=theindex+2; l_index<Prob->gL(); l_index++){
594                   if (denBK->gIrrep(l_index) == denBK->gIrrep(theindex)){
595                      double alpha = Prob->gMxElement(theindex,theindex,theindex,l_index);
596                      double * BlockL = Lright[l_index-theindex-2]->gStorage(NR,TwoSR,IR,NR+1,TwoSRdown,IRdown);
597                      daxpy_(&size, &alpha, BlockL, &inc, temp, &inc);
598                   }
599                }
600 
601                dgemm_(&notr,&tran,&dimL,&dimRup,&dimRdown,&factor,memS+denS->gKappa2index(memSkappa),&dimL,temp,&dimRup,&beta,memHeff+denS->gKappa2index(ikappa),&dimL);
602             }
603          }
604       }
605    }
606 
607 }
608 
addDiagram3Land3G(const int ikappa,double * memS,double * memHeff,const Sobject * denS,TensorQ * Qright,TensorL ** Lright,double * temp) const609 void CheMPS2::Heff::addDiagram3Land3G(const int ikappa, double * memS, double * memHeff, const Sobject * denS, TensorQ * Qright, TensorL ** Lright, double * temp) const{
610 
611    int NL = denS->gNL(ikappa);
612    int TwoSL = denS->gTwoSL(ikappa);
613    int IL = denS->gIL(ikappa);
614    int N1 = denS->gN1(ikappa);
615    int N2 = denS->gN2(ikappa);
616    int TwoJ = denS->gTwoJ(ikappa);
617    int NR = denS->gNR(ikappa);
618    int TwoSR = denS->gTwoSR(ikappa);
619    int IR = denS->gIR(ikappa);
620 
621    int theindex = denS->gIndex();
622    int IRdown = Irreps::directProd(IR,denBK->gIrrep(theindex+1));
623    int TwoS1 = (N1==1)?1:0;
624 
625    int dimL   = denBK->gCurrentDim(theindex,  NL,TwoSL,IL);
626    int dimRup = denBK->gCurrentDim(theindex+2,NR,TwoSR,IR);
627 
628    if (N2==1){ //3L1A
629       for (int TwoSRdown=TwoSR-1; TwoSRdown<=TwoSR+1; TwoSRdown+=2){
630          if ((abs(TwoSL-TwoSRdown)<=TwoS1) && (TwoSRdown>=0)){
631             int dimRdown = denBK->gCurrentDim(theindex+2, NR-1, TwoSRdown, IRdown);
632             int memSkappa = denS->gKappa(NL,TwoSL,IL,N1,0,TwoS1,NR-1,TwoSRdown,IRdown);
633             if (memSkappa!=-1){
634                int fase = phase(TwoSL+TwoSR+TwoS1+1);
635                double factor = sqrt((TwoJ+1)*(TwoSR+1.0)) * fase * Wigner::wigner6j(TwoS1,TwoJ,1,TwoSR,TwoSRdown,TwoSL);
636                double beta = 1.0; //add
637                char notr = 'N';
638                double * BlockQ = Qright->gStorage(NR-1,TwoSRdown,IRdown,NR,TwoSR,IR);
639                dgemm_(&notr,&notr,&dimL,&dimRup,&dimRdown,&factor,memS+denS->gKappa2index(memSkappa),&dimL,BlockQ,&dimRdown,&beta,memHeff+denS->gKappa2index(ikappa),&dimL);
640             }
641          }
642       }
643    }
644 
645    if (N2==2){ //3L1B and 3G1
646       for (int TwoSRdown=TwoSR-1; TwoSRdown<=TwoSR+1; TwoSRdown+=2){
647 
648          int dimRdown = denBK->gCurrentDim(theindex+2, NR-1, TwoSRdown, IRdown);
649          if (dimRdown>0){
650 
651             int TwoJstart = ((TwoSRdown!=TwoSL) || (TwoS1==0)) ? 1 + TwoS1 : 0;
652             for (int TwoJdown=TwoJstart; TwoJdown<=1+TwoS1; TwoJdown+=2){
653                if (abs(TwoSL-TwoSRdown)<=TwoJdown){
654 
655                   int memSkappa = denS->gKappa(NL,TwoSL,IL,N1,1,TwoJdown,NR-1,TwoSRdown,IRdown);
656                   if (memSkappa!=-1){
657                      int fase = phase(TwoSL+TwoSR+TwoS1+2);
658                      double factor = sqrt((TwoJdown+1)*(TwoSR+1.0)) * fase * Wigner::wigner6j(TwoJdown,TwoS1,1,TwoSR,TwoSRdown,TwoSL);
659                      double beta = 1.0; //add
660                      char notr = 'N';
661 
662                      double * BlockQ = Qright->gStorage(NR-1,TwoSRdown,IRdown,NR,TwoSR,IR);
663                      int inc = 1;
664                      int size = dimRup * dimRdown;
665                      dcopy_(&size,BlockQ,&inc,temp,&inc);
666 
667                      for (int l_index=theindex+2; l_index<Prob->gL(); l_index++){
668                         if (denBK->gIrrep(l_index) == denBK->gIrrep(theindex+1)){
669                            double alpha = Prob->gMxElement(theindex+1,theindex+1,theindex+1,l_index);
670                            double * BlockL = Lright[l_index-theindex-2]->gStorage(NR-1,TwoSRdown,IRdown,NR,TwoSR,IR);
671                            daxpy_(&size, &alpha, BlockL, &inc, temp, &inc);
672                         }
673                      }
674 
675                      dgemm_(&notr,&notr,&dimL,&dimRup,&dimRdown,&factor,memS+denS->gKappa2index(memSkappa),&dimL,temp,&dimRdown,&beta,memHeff+denS->gKappa2index(ikappa),&dimL);
676                   }
677                }
678             }
679          }
680       }
681    }
682 
683    if (N2==0){ //3L2A
684       for (int TwoSRdown=TwoSR-1; TwoSRdown<=TwoSR+1; TwoSRdown+=2){
685 
686          int dimRdown = denBK->gCurrentDim(theindex+2, NR+1, TwoSRdown, IRdown);
687          if (dimRdown>0){
688 
689             int TwoJstart = ((TwoSRdown!=TwoSL) || (TwoS1==0)) ? 1 + TwoS1 : 0;
690             for (int TwoJdown=TwoJstart; TwoJdown<=1+TwoS1; TwoJdown+=2){
691                if (abs(TwoSL-TwoSRdown)<=TwoJdown){
692 
693                   int memSkappa = denS->gKappa(NL,TwoSL,IL,N1,1,TwoJdown,NR+1,TwoSRdown,IRdown);
694                   if (memSkappa!=-1){
695                      int fase = phase(TwoSL+TwoSRdown+TwoS1+1);
696                      double factor = sqrt((TwoJdown+1)*(TwoSRdown+1.0)) * fase * Wigner::wigner6j(TwoJdown,TwoS1,1,TwoSR,TwoSRdown,TwoSL);
697                      double beta = 1.0; //add
698                      char notr = 'N';
699                      char tran = 'T';
700                      double * BlockQ = Qright->gStorage(NR,TwoSR,IR,NR+1,TwoSRdown,IRdown);
701                      dgemm_(&notr,&tran,&dimL,&dimRup,&dimRdown,&factor,memS+denS->gKappa2index(memSkappa),&dimL,BlockQ,&dimRup,&beta,memHeff+denS->gKappa2index(ikappa),&dimL);
702                   }
703                }
704             }
705          }
706       }
707    }
708 
709    if (N2==1){ //3L2B and 3G2
710       for (int TwoSRdown=TwoSR-1; TwoSRdown<=TwoSR+1; TwoSRdown+=2){
711          if ((abs(TwoSL-TwoSRdown)<=TwoS1) && (TwoSRdown>=0)){
712             int dimRdown = denBK->gCurrentDim(theindex+2, NR+1, TwoSRdown, IRdown);
713             int memSkappa = denS->gKappa(NL,TwoSL,IL,N1,2,TwoS1,NR+1,TwoSRdown,IRdown);
714             if (memSkappa!=-1){
715                int fase = phase(TwoSL+TwoSRdown+TwoS1+2);
716                double factor = sqrt((TwoJ+1)*(TwoSRdown+1.0)) * fase * Wigner::wigner6j(TwoS1,TwoJ,1,TwoSR,TwoSRdown,TwoSL);
717                double beta = 1.0; //add
718                char notr = 'N';
719                char tran = 'T';
720 
721                double * BlockQ = Qright->gStorage(NR,TwoSR,IR,NR+1,TwoSRdown,IRdown);
722                int inc = 1;
723                int size = dimRup * dimRdown;
724                dcopy_(&size,BlockQ,&inc,temp,&inc);
725 
726                for (int l_index=theindex+2; l_index<Prob->gL(); l_index++){
727                   if (denBK->gIrrep(l_index) == denBK->gIrrep(theindex+1)){
728                      double alpha = Prob->gMxElement(theindex+1,theindex+1,theindex+1,l_index);
729                      double * BlockL = Lright[l_index-theindex-2]->gStorage(NR,TwoSR,IR,NR+1,TwoSRdown,IRdown);
730                      daxpy_(&size, &alpha, BlockL, &inc, temp, &inc);
731                   }
732                }
733 
734                dgemm_(&notr,&tran,&dimL,&dimRup,&dimRdown,&factor,memS+denS->gKappa2index(memSkappa),&dimL,temp,&dimRup,&beta,memHeff+denS->gKappa2index(ikappa),&dimL);
735             }
736          }
737       }
738    }
739 
740 }
741 
addDiagram3J(const int ikappa,double * memS,double * memHeff,const Sobject * denS,TensorQ ** Qright,TensorL ** Lleft,double * temp) const742 void CheMPS2::Heff::addDiagram3J(const int ikappa, double * memS, double * memHeff, const Sobject * denS, TensorQ ** Qright, TensorL ** Lleft, double * temp) const{
743 
744    #ifdef CHEMPS2_MPI_COMPILATION
745    const int MPIRANK = MPIchemps2::mpi_rank();
746    #endif
747 
748    int NL = denS->gNL(ikappa);
749    int TwoSL = denS->gTwoSL(ikappa);
750    int IL = denS->gIL(ikappa);
751    int N1 = denS->gN1(ikappa);
752    int N2 = denS->gN2(ikappa);
753    int TwoJ = denS->gTwoJ(ikappa);
754    int NR = denS->gNR(ikappa);
755    int TwoSR = denS->gTwoSR(ikappa);
756    int IR = denS->gIR(ikappa);
757 
758    int theindex = denS->gIndex();
759 
760    int dimRup = denBK->gCurrentDim(theindex+2,NR,TwoSR,IR);
761    int dimLup = denBK->gCurrentDim(theindex,  NL,TwoSL,IL);
762 
763    //First do 3J2
764    for (int TwoSLdown=TwoSL-1; TwoSLdown<=TwoSL+1; TwoSLdown+=2){
765       for (int TwoSRdown=TwoSR-1; TwoSRdown<=TwoSR+1; TwoSRdown+=2){
766          if ((abs(TwoSLdown-TwoSRdown)<=TwoJ) && (TwoSLdown>=0) && (TwoSRdown>=0)){
767 
768             int fase = phase(TwoSLdown+TwoSR+TwoJ+1 + ((N1==1)?2:0) + ((N2==1)?2:0) );
769             const double factor = fase * sqrt((TwoSLdown+1)*(TwoSRdown+1.0)) * Wigner::wigner6j(TwoSL,TwoSR,TwoJ,TwoSRdown,TwoSLdown,1);
770 
771             for (int l_index=0; l_index<theindex; l_index++){
772 
773                #ifdef CHEMPS2_MPI_COMPILATION
774                if ( MPIchemps2::owner_q( Prob->gL(), l_index ) == MPIRANK )
775                #endif
776                {
777                   int ILdown = Irreps::directProd(IL,denBK->gIrrep(l_index));
778                   int IRdown = Irreps::directProd(IR,denBK->gIrrep(l_index));
779                   int memSkappa = denS->gKappa(NL+1, TwoSLdown, ILdown, N1, N2, TwoJ, NR+1, TwoSRdown, IRdown);
780                   if (memSkappa!=-1){
781 
782                      int dimRdown = denBK->gCurrentDim(theindex+2, NR+1, TwoSRdown, IRdown);
783                      int dimLdown = denBK->gCurrentDim(theindex,   NL+1, TwoSLdown, ILdown);
784 
785                      double * Lblock = Lleft[ theindex-1-l_index]->gStorage(NL,TwoSL,IL,NL+1,TwoSLdown,ILdown);
786                      double * Qblock = Qright[theindex+1-l_index]->gStorage(NR,TwoSR,IR,NR+1,TwoSRdown,IRdown);
787 
788                      char trans = 'T';
789                      char notra = 'N';
790                      double beta = 0.0; //set
791                      double alpha = factor;
792                      dgemm_(&notra,&notra,&dimLup,&dimRdown,&dimLdown,&alpha,Lblock,&dimLup,memS+denS->gKappa2index(memSkappa),&dimLdown,&beta,temp,&dimLup);
793 
794                      beta = 1.0; //add
795                      alpha = 1.0;
796                      dgemm_(&notra,&trans,&dimLup,&dimRup,&dimRdown,&alpha,temp,&dimLup,Qblock,&dimRup,&beta,memHeff+denS->gKappa2index(ikappa),&dimLup);
797 
798                   }
799                }
800             }
801          }
802       }
803    }
804 
805    //Then do 3J1
806    for (int TwoSLdown=TwoSL-1; TwoSLdown<=TwoSL+1; TwoSLdown+=2){
807       for (int TwoSRdown=TwoSR-1; TwoSRdown<=TwoSR+1; TwoSRdown+=2){
808          if ((abs(TwoSLdown-TwoSRdown)<=TwoJ) && (TwoSLdown>=0) && (TwoSRdown>=0)){
809 
810             int fase = phase(TwoSL+TwoSRdown+TwoJ+1 + ((N1==1)?2:0) + ((N2==1)?2:0) );
811             const double factor = fase * sqrt((TwoSL+1)*(TwoSR+1.0)) * Wigner::wigner6j(TwoSL,TwoSR,TwoJ,TwoSRdown,TwoSLdown,1);
812 
813             for (int l_index=0; l_index<theindex; l_index++){
814 
815                #ifdef CHEMPS2_MPI_COMPILATION
816                if ( MPIchemps2::owner_q( Prob->gL(), l_index ) == MPIRANK )
817                #endif
818                {
819                   int ILdown = Irreps::directProd(IL,denBK->gIrrep(l_index));
820                   int IRdown = Irreps::directProd(IR,denBK->gIrrep(l_index));
821                   int memSkappa = denS->gKappa(NL-1, TwoSLdown, ILdown, N1, N2, TwoJ, NR-1, TwoSRdown, IRdown);
822                   if (memSkappa!=-1){
823                      int dimRdown = denBK->gCurrentDim(theindex+2, NR-1, TwoSRdown, IRdown);
824                      int dimLdown = denBK->gCurrentDim(theindex,   NL-1, TwoSLdown, ILdown);
825 
826                      double * Lblock = Lleft[ theindex-1-l_index]->gStorage(NL-1,TwoSLdown,ILdown,NL,TwoSL,IL);
827                      double * Qblock = Qright[theindex+1-l_index]->gStorage(NR-1,TwoSRdown,IRdown,NR,TwoSR,IR);
828 
829                      char trans = 'T';
830                      char notra = 'N';
831                      double beta = 0.0; //set
832                      double alpha = factor;
833                      dgemm_(&trans,&notra,&dimLup,&dimRdown,&dimLdown,&alpha,Lblock,&dimLdown,memS+denS->gKappa2index(memSkappa),&dimLdown,&beta,temp,&dimLup);
834 
835                      beta = 1.0; //add
836                      alpha = 1.0;
837                      dgemm_(&notra,&notra,&dimLup,&dimRup,&dimRdown,&alpha,temp,&dimLup,Qblock,&dimRdown,&beta,memHeff+denS->gKappa2index(ikappa),&dimLup);
838 
839                   }
840                }
841             }
842          }
843       }
844    }
845 
846 }
847 
848 
849