CheMPS2
HeffDiagrams3.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::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 
161 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 
294 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 
399 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 
476 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 
609 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 
742 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 
CheMPS2::Wigner::wigner6j
static double wigner6j(const int two_ja, const int two_jb, const int two_jc, const int two_jd, const int two_je, const int two_jf)
Wigner-6j symbol (gsl api)
Definition: Wigner.cpp:294
CheMPS2::SyBookkeeper::gCurrentDim
int gCurrentDim(const int boundary, const int N, const int TwoS, const int irrep) const
Get the current virtual dimensions.
Definition: SyBookkeeper.cpp:165
CheMPS2::MPIchemps2::mpi_rank
static int mpi_rank()
Get the rank of this MPI process.
Definition: MPIchemps2.h:131
CheMPS2::Problem::gMxElement
double gMxElement(const int alpha, const int beta, const int gamma, const int delta) const
Get a specific interaction matrix element.
Definition: Problem.cpp:350
CheMPS2::Irreps::directProd
static int directProd(const int Irrep1, const int Irrep2)
Get the direct product of the irreps with numbers Irrep1 and Irrep2: a bitwise XOR for psi4&#39;s convent...
Definition: Irreps.h:123
CheMPS2::Problem::gL
int gL() const
Get the number of orbitals.
Definition: Problem.h:51
CheMPS2::Heff::phase
static int phase(const int TwoTimesPower)
Phase function.
Definition: Heff.h:75
CheMPS2::SyBookkeeper::gIrrep
int gIrrep(const int orbital) const
Get an orbital irrep.
Definition: SyBookkeeper.cpp:145