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 <stdlib.h>
21 #include <math.h>
22 
23 #include "TensorQ.h"
24 #include "Lapack.h"
25 #include "Wigner.h"
26 
TensorQ(const int boundary_index,const int Idiff,const bool moving_right,const SyBookkeeper * denBK,const Problem * Prob,const int site)27 CheMPS2::TensorQ::TensorQ(const int boundary_index, const int Idiff, const bool moving_right, const SyBookkeeper * denBK, const Problem * Prob, const int site) :
28 TensorOperator(boundary_index,
29                1, //two_j
30                1, //n_elec
31                Idiff,
32                moving_right,
33                true, //prime_last
34                true, //jw_phase (three 2nd quantized operators)
35                denBK,
36                denBK){
37 
38    this->Prob = Prob;
39    this->site = site;
40 
41 }
42 
~TensorQ()43 CheMPS2::TensorQ::~TensorQ(){ }
44 
AddTermSimple(TensorT * denT)45 void CheMPS2::TensorQ::AddTermSimple(TensorT * denT){
46 
47    if (( moving_right) && (bk_up->gIrrep(denT->gIndex()) == n_irrep)){ AddTermSimpleRight(denT); }
48    if ((!moving_right) && (bk_up->gIrrep(denT->gIndex()) == n_irrep)){ AddTermSimpleLeft( denT); }
49 
50 }
51 
AddTermSimpleRight(TensorT * denT)52 void CheMPS2::TensorQ::AddTermSimpleRight(TensorT * denT){
53 
54    const double mxElement = Prob->gMxElement(index-1,index-1,index-1,site);
55 
56    for (int ikappa=0; ikappa<nKappa; ikappa++){
57       const int ID = Irreps::directProd( n_irrep, sector_irrep_up[ikappa] );
58       int dimRU = bk_up->gCurrentDim(index,   sector_nelec_up[ikappa],   sector_spin_up[ikappa],    sector_irrep_up[ikappa]);
59       int dimRD = bk_up->gCurrentDim(index,   sector_nelec_up[ikappa]+1, sector_spin_down[ikappa], ID);
60       int dimL  = bk_up->gCurrentDim(index-1, sector_nelec_up[ikappa]-1, sector_spin_down[ikappa], ID);
61       if (dimL>0){
62 
63          double * BlockTup = denT->gStorage(sector_nelec_up[ikappa]-1, sector_spin_down[ikappa], ID, sector_nelec_up[ikappa],   sector_spin_up[ikappa],    sector_irrep_up[ikappa]);
64          double * BlockTdo = denT->gStorage(sector_nelec_up[ikappa]-1, sector_spin_down[ikappa], ID, sector_nelec_up[ikappa]+1, sector_spin_down[ikappa], ID);
65 
66          int fase = ((((sector_spin_down[ikappa]+1-sector_spin_up[ikappa])/2)%2)!=0)?-1:1;
67          double alpha = fase * mxElement * sqrt((sector_spin_up[ikappa]+1.0)/(sector_spin_down[ikappa]+1.0));
68          double beta = 1.0; //add
69          char trans = 'T';
70          char notr = 'N';
71 
72          dgemm_(&trans,&notr,&dimRU,&dimRD,&dimL,&alpha,BlockTup,&dimL,BlockTdo,&dimL,&beta,storage+kappa2index[ikappa],&dimRU);
73 
74       }
75    }
76 
77 }
78 
AddTermSimpleLeft(TensorT * denT)79 void CheMPS2::TensorQ::AddTermSimpleLeft(TensorT * denT){
80 
81    const double mxElement = Prob->gMxElement(site,index,index,index);
82 
83    for (int ikappa=0; ikappa<nKappa; ikappa++){
84       const int ID = Irreps::directProd( n_irrep, sector_irrep_up[ikappa] );
85       int dimLU = bk_up->gCurrentDim(index,   sector_nelec_up[ikappa],   sector_spin_up[ikappa],    sector_irrep_up[ikappa]);
86       int dimLD = bk_up->gCurrentDim(index,   sector_nelec_up[ikappa]+1, sector_spin_down[ikappa], ID);
87       int dimR  = bk_up->gCurrentDim(index+1, sector_nelec_up[ikappa]+2, sector_spin_up[ikappa],    sector_irrep_up[ikappa]);
88       if (dimR>0){
89 
90          double * BlockTup = denT->gStorage(sector_nelec_up[ikappa],   sector_spin_up[ikappa],    sector_irrep_up[ikappa], sector_nelec_up[ikappa]+2, sector_spin_up[ikappa], sector_irrep_up[ikappa]);
91          double * BlockTdo = denT->gStorage(sector_nelec_up[ikappa]+1, sector_spin_down[ikappa], ID,               sector_nelec_up[ikappa]+2, sector_spin_up[ikappa], sector_irrep_up[ikappa]);
92 
93          int fase = ((((sector_spin_up[ikappa]+1-sector_spin_down[ikappa])/2)%2)!=0)?-1:1;
94          double alpha = fase * mxElement * sqrt((sector_spin_up[ikappa]+1.0)/(sector_spin_down[ikappa]+1.0));
95          double beta = 1.0; //add
96          char trans = 'T';
97          char notr = 'N';
98 
99          dgemm_(&notr,&trans,&dimLU,&dimLD,&dimR,&alpha,BlockTup,&dimLU,BlockTdo,&dimLD,&beta,storage+kappa2index[ikappa],&dimLU);
100 
101       }
102    }
103 
104 }
105 
AddTermsL(TensorL ** Ltensors,TensorT * denT,double * workmem,double * workmem2)106 void CheMPS2::TensorQ::AddTermsL(TensorL ** Ltensors, TensorT * denT, double * workmem, double * workmem2){
107 
108    if (moving_right){ AddTermsLRight(Ltensors, denT, workmem, workmem2); }
109    else{ AddTermsLLeft( Ltensors, denT, workmem, workmem2); }
110 
111 }
112 
AddTermsLRight(TensorL ** Ltensors,TensorT * denT,double * workmem,double * workmem2)113 void CheMPS2::TensorQ::AddTermsLRight(TensorL ** Ltensors, TensorT * denT, double * workmem, double * workmem2){
114 
115    bool OneToAdd = false;
116    for (int loca=0; loca<index-1; loca++){
117       if (Ltensors[index-2-loca]->get_irrep() == n_irrep){ OneToAdd = true; }
118    }
119 
120    if (OneToAdd){
121       for (int ikappa=0; ikappa<nKappa; ikappa++){
122 
123          const int ID = Irreps::directProd( sector_irrep_up[ikappa], n_irrep );
124          int dimRU = bk_up->gCurrentDim(index,   sector_nelec_up[ikappa],   sector_spin_up[ikappa],    sector_irrep_up[ikappa]);
125          int dimRD = bk_up->gCurrentDim(index,   sector_nelec_up[ikappa]+1, sector_spin_down[ikappa], ID);
126 
127          //case 1
128          int dimLU = bk_up->gCurrentDim(index-1, sector_nelec_up[ikappa],   sector_spin_up[ikappa],    sector_irrep_up[ikappa]);
129          int dimLD = bk_up->gCurrentDim(index-1, sector_nelec_up[ikappa]-1, sector_spin_down[ikappa], ID);
130 
131          if ((dimLU>0) && (dimLD>0)){
132 
133             int dimLUxLD = dimLU * dimLD;
134             for (int cnt=0; cnt<dimLUxLD; cnt++){ workmem[cnt] = 0.0; }
135 
136             for (int loca=0; loca<index-1; loca++){
137                if (Ltensors[index-2-loca]->get_irrep() == n_irrep){
138                   double * BlockL = Ltensors[index-2-loca]->gStorage(sector_nelec_up[ikappa]-1,sector_spin_down[ikappa],ID,sector_nelec_up[ikappa],sector_spin_up[ikappa],sector_irrep_up[ikappa]);
139                   double alpha = Prob->gMxElement(loca,site,index-1,index-1);
140                   int inc = 1;
141                   daxpy_(&dimLUxLD, &alpha, BlockL, &inc, workmem, &inc);
142                }
143             }
144 
145             int fase = ((((sector_spin_down[ikappa]+1-sector_spin_up[ikappa])/2)%2)!=0)?-1:1;
146             double alpha = fase * sqrt((sector_spin_up[ikappa]+1.0)/(sector_spin_down[ikappa]+1.0));
147             double beta = 0.0; //set
148 
149             double * BlockTup = denT->gStorage(sector_nelec_up[ikappa], sector_spin_up[ikappa], sector_irrep_up[ikappa], sector_nelec_up[ikappa], sector_spin_up[ikappa], sector_irrep_up[ikappa]);
150             double * BlockTdo = denT->gStorage(sector_nelec_up[ikappa]-1, sector_spin_down[ikappa], ID, sector_nelec_up[ikappa]+1, sector_spin_down[ikappa], ID);
151 
152             char totrans = 'T';
153             // factor * Tup^T * L^T --> mem2
154             dgemm_(&totrans,&totrans,&dimRU,&dimLD,&dimLU,&alpha,BlockTup,&dimLU,workmem,&dimLD,&beta,workmem2,&dimRU);
155 
156             alpha = 1.0;
157             beta = 1.0; //add
158             totrans = 'N';
159             // mem2 * Tdo --> storage
160             dgemm_(&totrans,&totrans,&dimRU,&dimRD,&dimLD,&alpha,workmem2,&dimRU, BlockTdo, &dimLD, &beta, storage+kappa2index[ikappa], &dimRU);
161 
162          }
163 
164          //case 2
165          dimLU = bk_up->gCurrentDim(index-1, sector_nelec_up[ikappa]-2, sector_spin_up[ikappa], sector_irrep_up[ikappa]);
166          //dimLD same as case1
167          if ((dimLU>0) && (dimLD>0)){
168 
169             int dimLUxLD = dimLU * dimLD;
170             for (int cnt=0; cnt<dimLUxLD; cnt++){ workmem[cnt] = 0.0; }
171 
172             for (int loca=0; loca<index-1; loca++){
173                if (Ltensors[index-2-loca]->get_irrep() == n_irrep){
174                   double * BlockL = Ltensors[index-2-loca]->gStorage(sector_nelec_up[ikappa]-2,sector_spin_up[ikappa],sector_irrep_up[ikappa],sector_nelec_up[ikappa]-1,sector_spin_down[ikappa],ID);
175                   double alpha = 2*Prob->gMxElement(loca,index-1,site,index-1) - Prob->gMxElement(loca,index-1,index-1,site);
176                   int inc = 1;
177                   daxpy_(&dimLUxLD, &alpha, BlockL, &inc, workmem, &inc);
178                }
179             }
180 
181             double alpha = 1.0; //factor = 1 in this case
182             double beta = 0.0; //set
183 
184             double * BlockTup = denT->gStorage(sector_nelec_up[ikappa]-2, sector_spin_up[ikappa],    sector_irrep_up[ikappa], sector_nelec_up[ikappa],   sector_spin_up[ikappa],    sector_irrep_up[ikappa]);
185             double * BlockTdo = denT->gStorage(sector_nelec_up[ikappa]-1, sector_spin_down[ikappa], ID,               sector_nelec_up[ikappa]+1, sector_spin_down[ikappa], ID);
186 
187             char trans = 'T';
188             char notr = 'N';
189             // factor * Tup^T * L --> mem2
190             dgemm_(&trans,&notr,&dimRU,&dimLD,&dimLU,&alpha,BlockTup,&dimLU,workmem,&dimLU,&beta,workmem2,&dimRU);
191 
192             beta = 1.0; //add
193             // mem2 * Tdo --> storage
194             dgemm_(&notr,&notr,&dimRU,&dimRD,&dimLD,&alpha,workmem2,&dimRU, BlockTdo, &dimLD, &beta, storage+kappa2index[ikappa], &dimRU);
195 
196          }
197 
198          //case 3
199          for (int TwoSLU=sector_spin_up[ikappa]-1; TwoSLU<=sector_spin_up[ikappa]+1; TwoSLU+=2){
200             for (int TwoSLD=sector_spin_down[ikappa]-1; TwoSLD<=sector_spin_down[ikappa]+1; TwoSLD+=2){
201                if ((TwoSLD>=0) && (TwoSLU>=0) && (abs(TwoSLD-TwoSLU)<2)){
202                   const int ILU = Irreps::directProd(sector_irrep_up[ikappa],bk_up->gIrrep(index-1));
203                   const int ILD = Irreps::directProd(ID,              bk_up->gIrrep(index-1));
204                   dimLU = bk_up->gCurrentDim(index-1, sector_nelec_up[ikappa]-1, TwoSLU, ILU);
205                   dimLD = bk_up->gCurrentDim(index-1, sector_nelec_up[ikappa]  , TwoSLD, ILD);
206                   if ((dimLU>0) && (dimLD>0)){
207                      int fase = ((((sector_spin_up[ikappa]+TwoSLD)/2)%2)!=0)?-1:1;
208                      double factor = fase * sqrt((TwoSLD+1)*(sector_spin_up[ikappa]+1.0))
209                                    * Wigner::wigner6j( sector_spin_up[ikappa], sector_spin_down[ikappa], 1, TwoSLD, TwoSLU, 1 );
210 
211                      int dimLUxLD = dimLU * dimLD;
212                      for (int cnt=0; cnt<dimLUxLD; cnt++){ workmem[cnt] = 0.0; }
213 
214                      for (int loca=0; loca<index-1; loca++){
215                         if (Ltensors[index-2-loca]->get_irrep() == n_irrep){
216                            double * BlockL = Ltensors[index-2-loca]->gStorage(sector_nelec_up[ikappa]-1, TwoSLU, ILU, sector_nelec_up[ikappa], TwoSLD, ILD);
217                            double alpha = factor * Prob->gMxElement(loca,index-1,site,index-1);
218                            if (TwoSLD==sector_spin_up[ikappa]){ alpha += Prob->gMxElement(loca,index-1,index-1,site); }
219                            int inc = 1;
220                            daxpy_(&dimLUxLD, &alpha, BlockL, &inc, workmem, &inc);
221                         }
222                      }
223 
224                      double alpha = 1.0;
225                      double beta = 0.0; //set
226 
227                      double * BlockTup = denT->gStorage(sector_nelec_up[ikappa]-1, TwoSLU, ILU, sector_nelec_up[ikappa],   sector_spin_up[ikappa],    sector_irrep_up[ikappa]);
228                      double * BlockTdo = denT->gStorage(sector_nelec_up[ikappa]  , TwoSLD, ILD, sector_nelec_up[ikappa]+1, sector_spin_down[ikappa], ID);
229 
230                      char trans = 'T';
231                      char notr = 'N';
232                      // Tup^T * mem --> mem2
233                      dgemm_(&trans,&notr,&dimRU,&dimLD,&dimLU,&alpha,BlockTup,&dimLU,workmem,&dimLU,&beta,workmem2,&dimRU);
234 
235                      beta = 1.0; //add
236                      // mem2 * Tdo --> storage
237                      dgemm_(&notr,&notr,&dimRU,&dimRD,&dimLD,&alpha,workmem2,&dimRU, BlockTdo, &dimLD, &beta, storage+kappa2index[ikappa], &dimRU);
238 
239                   }
240                }
241             }
242          }
243 
244       }
245    }
246 
247 }
248 
AddTermsLLeft(TensorL ** Ltensors,TensorT * denT,double * workmem,double * workmem2)249 void CheMPS2::TensorQ::AddTermsLLeft(TensorL ** Ltensors, TensorT * denT, double * workmem, double * workmem2){
250 
251    bool OneToAdd = false;
252    for (int loca=index+1; loca<Prob->gL(); loca++){
253       if (Ltensors[loca-index-1]->get_irrep() == n_irrep){ OneToAdd = true; }
254    }
255 
256    if (OneToAdd){
257       for (int ikappa=0; ikappa<nKappa; ikappa++){
258 
259          const int ID = Irreps::directProd( sector_irrep_up[ikappa], n_irrep );
260          int dimLU = bk_up->gCurrentDim(index,   sector_nelec_up[ikappa],   sector_spin_up[ikappa],    sector_irrep_up[ikappa]);
261          int dimLD = bk_up->gCurrentDim(index,   sector_nelec_up[ikappa]+1, sector_spin_down[ikappa], ID);
262 
263          //case 1
264          int dimRU = bk_up->gCurrentDim(index+1, sector_nelec_up[ikappa]+2, sector_spin_up[ikappa],    sector_irrep_up[ikappa]);
265          int dimRD = bk_up->gCurrentDim(index+1, sector_nelec_up[ikappa]+1, sector_spin_down[ikappa], ID);
266 
267          if ((dimRU>0) && (dimRD>0)){
268 
269             int dimRUxRD = dimRU * dimRD;
270             for (int cnt=0; cnt<dimRUxRD; cnt++){ workmem[cnt] = 0.0; }
271 
272             for (int loca=index+1; loca<Prob->gL(); loca++){
273                if (Ltensors[loca-index-1]->get_irrep() == n_irrep){
274                   double * BlockL = Ltensors[loca-index-1]->gStorage(sector_nelec_up[ikappa]+1,sector_spin_down[ikappa],ID,sector_nelec_up[ikappa]+2,sector_spin_up[ikappa],sector_irrep_up[ikappa]);
275                   double alpha = Prob->gMxElement(site,loca,index,index);
276                   int inc = 1;
277                   daxpy_(&dimRUxRD, &alpha, BlockL, &inc, workmem, &inc);
278                }
279             }
280 
281             int fase = ((((sector_spin_up[ikappa]+1-sector_spin_down[ikappa])/2)%2)!=0)?-1:1;
282             double alpha = fase * sqrt((sector_spin_up[ikappa]+1.0)/(sector_spin_down[ikappa]+1.0));
283             double beta = 0.0; //set
284 
285             double * BlockTup = denT->gStorage(sector_nelec_up[ikappa],   sector_spin_up[ikappa],    sector_irrep_up[ikappa], sector_nelec_up[ikappa]+2, sector_spin_up[ikappa],    sector_irrep_up[ikappa]);
286             double * BlockTdo = denT->gStorage(sector_nelec_up[ikappa]+1, sector_spin_down[ikappa], ID,               sector_nelec_up[ikappa]+1, sector_spin_down[ikappa], ID);
287 
288             char trans = 'T';
289             char notr = 'N';
290             // factor * Tup * L^T --> mem2
291             dgemm_(&notr,&trans,&dimLU,&dimRD,&dimRU,&alpha,BlockTup,&dimLU,workmem,&dimRD,&beta,workmem2,&dimLU);
292 
293             alpha = 1.0;
294             beta = 1.0; //add
295             // mem2 * Tdo^T --> storage
296             dgemm_(&notr,&trans,&dimLU,&dimLD,&dimRD,&alpha,workmem2,&dimLU, BlockTdo, &dimLD, &beta, storage+kappa2index[ikappa], &dimLU);
297 
298          }
299 
300          //case 2
301          dimRD = bk_up->gCurrentDim(index+1, sector_nelec_up[ikappa]+3, sector_spin_down[ikappa], ID);
302          //dimRU same as case1
303          if ((dimRU>0) && (dimRD>0)){
304 
305             int dimRUxRD = dimRU * dimRD;
306             for (int cnt=0; cnt<dimRUxRD; cnt++){ workmem[cnt] = 0.0; }
307 
308             for (int loca=index+1; loca<Prob->gL(); loca++){
309                if (Ltensors[loca-index-1]->get_irrep() == n_irrep){
310                   double * BlockL = Ltensors[loca-index-1]->gStorage(sector_nelec_up[ikappa]+2,sector_spin_up[ikappa],sector_irrep_up[ikappa],sector_nelec_up[ikappa]+3,sector_spin_down[ikappa],ID);
311                   double alpha = 2*Prob->gMxElement(site,index,loca,index) - Prob->gMxElement(site,index,index,loca);
312                   int inc = 1;
313                   daxpy_(&dimRUxRD, &alpha, BlockL, &inc, workmem, &inc);
314                }
315             }
316 
317             double alpha = 1.0; //factor = 1 in this case
318             double beta = 0.0; //set
319 
320             double * BlockTup = denT->gStorage(sector_nelec_up[ikappa],   sector_spin_up[ikappa],    sector_irrep_up[ikappa], sector_nelec_up[ikappa]+2, sector_spin_up[ikappa],    sector_irrep_up[ikappa]);
321             double * BlockTdo = denT->gStorage(sector_nelec_up[ikappa]+1, sector_spin_down[ikappa], ID,               sector_nelec_up[ikappa]+3, sector_spin_down[ikappa], ID);
322 
323             char notr = 'N';
324             // factor * Tup * L --> mem2
325             dgemm_(&notr,&notr,&dimLU,&dimRD,&dimRU,&alpha,BlockTup,&dimLU,workmem,&dimRU,&beta,workmem2,&dimLU);
326 
327             beta = 1.0; //add
328             // mem2 * Tdo^T --> storage
329             char trans = 'T';
330             dgemm_(&notr,&trans,&dimLU,&dimLD,&dimRD,&alpha,workmem2,&dimLU, BlockTdo, &dimLD, &beta, storage+kappa2index[ikappa], &dimLU);
331 
332          }
333 
334          //case 3
335          for (int TwoSRU=sector_spin_up[ikappa]-1; TwoSRU<=sector_spin_up[ikappa]+1; TwoSRU+=2){
336             for (int TwoSRD=sector_spin_down[ikappa]-1; TwoSRD<=sector_spin_down[ikappa]+1; TwoSRD+=2){
337                if ((TwoSRD>=0) && (TwoSRU>=0) && (abs(TwoSRD-TwoSRU)<2)){
338                   const int IRU = Irreps::directProd(sector_irrep_up[ikappa],bk_up->gIrrep(index));
339                   const int IRD = Irreps::directProd(ID,              bk_up->gIrrep(index));
340                   dimRU = bk_up->gCurrentDim(index+1, sector_nelec_up[ikappa]+1, TwoSRU, IRU);
341                   dimRD = bk_up->gCurrentDim(index+1, sector_nelec_up[ikappa]+2, TwoSRD, IRD);
342                   if ((dimRU>0) && (dimRD>0)){
343                      int fase = ((((sector_spin_down[ikappa]+TwoSRU)/2)%2)!=0)?-1:1;
344                      double factor1 = fase * sqrt((TwoSRU+1.0)/(sector_spin_down[ikappa]+1.0)) * (TwoSRD+1)
345                                     * Wigner::wigner6j( sector_spin_up[ikappa], sector_spin_down[ikappa], 1, TwoSRD, TwoSRU, 1 );
346                      double factor2 = (TwoSRD+1.0)/(sector_spin_down[ikappa]+1.0);
347 
348                      int dimRUxRD = dimRU * dimRD;
349                      for (int cnt=0; cnt<dimRUxRD; cnt++){ workmem[cnt] = 0.0; }
350 
351                      for (int loca=index+1; loca<Prob->gL(); loca++){
352                         if (Ltensors[loca-index-1]->get_irrep() == n_irrep){
353                            double * BlockL = Ltensors[loca-index-1]->gStorage(sector_nelec_up[ikappa]+1, TwoSRU, IRU, sector_nelec_up[ikappa]+2, TwoSRD, IRD);
354                            double alpha = factor1 * Prob->gMxElement(site,index,loca,index);
355                            if (TwoSRU==sector_spin_down[ikappa]){ alpha += factor2 * Prob->gMxElement(site,index,index,loca); }
356                            int inc = 1;
357                            daxpy_(&dimRUxRD, &alpha, BlockL, &inc, workmem, &inc);
358                         }
359                      }
360 
361                      double alpha = 1.0;
362                      double beta = 0.0; //set
363 
364                      double * BlockTup = denT->gStorage(sector_nelec_up[ikappa],   sector_spin_up[ikappa],    sector_irrep_up[ikappa], sector_nelec_up[ikappa]+1, TwoSRU, IRU);
365                      double * BlockTdo = denT->gStorage(sector_nelec_up[ikappa]+1, sector_spin_down[ikappa], ID,               sector_nelec_up[ikappa]+2, TwoSRD, IRD);
366 
367                      char notr = 'N';
368                      // Tup * mem --> mem2
369                      dgemm_(&notr,&notr,&dimLU,&dimRD,&dimRU,&alpha,BlockTup,&dimLU,workmem,&dimRU,&beta,workmem2,&dimLU);
370 
371                      beta = 1.0; //add
372                      // mem2 * Tdo^T --> storage
373                      char trans = 'T';
374                      dgemm_(&notr,&trans,&dimLU,&dimLD,&dimRD,&alpha,workmem2,&dimLU, BlockTdo, &dimLD, &beta, storage+kappa2index[ikappa], &dimLU);
375 
376                   }
377                }
378             }
379          }
380 
381       }
382    }
383 
384 }
385 
AddTermsAB(TensorOperator * denA,TensorOperator * denB,TensorT * denT,double * workmem,double * workmem2)386 void CheMPS2::TensorQ::AddTermsAB(TensorOperator * denA, TensorOperator * denB, TensorT * denT, double * workmem, double * workmem2){
387 
388    if (moving_right){ AddTermsABRight(denA, denB, denT, workmem, workmem2); }
389    else{ AddTermsABLeft( denA, denB, denT, workmem, workmem2); }
390 
391 }
392 
AddTermsABRight(TensorOperator * denA,TensorOperator * denB,TensorT * denT,double * workmem,double * workmem2)393 void CheMPS2::TensorQ::AddTermsABRight(TensorOperator * denA, TensorOperator * denB, TensorT * denT, double * workmem, double * workmem2){
394 
395    for (int ikappa=0; ikappa<nKappa; ikappa++){
396       const int ID = Irreps::directProd( sector_irrep_up[ikappa], n_irrep );
397       int dimRU = bk_up->gCurrentDim(index, sector_nelec_up[ikappa],   sector_spin_up[ikappa],    sector_irrep_up[ikappa]);
398       int dimRD = bk_up->gCurrentDim(index, sector_nelec_up[ikappa]+1, sector_spin_down[ikappa], ID);
399 
400       //case 1
401       const int ILU = Irreps::directProd(sector_irrep_up[ikappa],bk_up->gIrrep(index-1));
402       for (int TwoSLU=sector_spin_up[ikappa]-1; TwoSLU<=sector_spin_up[ikappa]+1; TwoSLU+=2){
403          int dimLU = bk_up->gCurrentDim(index-1, sector_nelec_up[ikappa]-1, TwoSLU,                 ILU);
404          int dimLD = bk_up->gCurrentDim(index-1, sector_nelec_up[ikappa]+1, sector_spin_down[ikappa], ID);
405          if ((dimLU>0) && (dimLD>0)){
406 
407             int fase = ((((TwoSLU + sector_spin_down[ikappa] + 2)/2)%2)!=0)?-1:1;
408             double factorB = fase * sqrt(3.0*(sector_spin_up[ikappa]+1))
409                            * Wigner::wigner6j( 1, 2, 1, sector_spin_down[ikappa], sector_spin_up[ikappa], TwoSLU );
410 
411             double alpha;
412             double * mem;
413 
414             if (TwoSLU == sector_spin_down[ikappa]){
415 
416                fase = ((((sector_spin_down[ikappa]+1-sector_spin_up[ikappa])/2)%2)!=0)?-1:1;
417                double factorA = fase * sqrt( 0.5 * (sector_spin_up[ikappa]+1.0) / (sector_spin_down[ikappa]+1.0) );
418 
419                double * BlockA = denA->gStorage( sector_nelec_up[ikappa]-1,TwoSLU,ILU,sector_nelec_up[ikappa]+1,sector_spin_down[ikappa],ID );
420                double * BlockB = denB->gStorage( sector_nelec_up[ikappa]-1,TwoSLU,ILU,sector_nelec_up[ikappa]+1,sector_spin_down[ikappa],ID );
421 
422                mem = workmem;
423                for (int cnt=0; cnt<dimLU*dimLD; cnt++){ mem[cnt] = factorA * BlockA[cnt] + factorB * BlockB[cnt]; }
424                alpha = 1.0;
425 
426             } else {
427 
428                alpha = factorB;
429                mem = denB->gStorage(sector_nelec_up[ikappa]-1,TwoSLU,ILU,sector_nelec_up[ikappa]+1,sector_spin_down[ikappa],ID);
430 
431             }
432 
433             double * BlockTup = denT->gStorage(sector_nelec_up[ikappa]-1, TwoSLU,                 ILU, sector_nelec_up[ikappa],   sector_spin_up[ikappa],    sector_irrep_up[ikappa]);
434             double * BlockTdo = denT->gStorage(sector_nelec_up[ikappa]+1, sector_spin_down[ikappa], ID,  sector_nelec_up[ikappa]+1, sector_spin_down[ikappa], ID);
435 
436             char trans = 'T';
437             char notr = 'N';
438             double beta = 0.0; //set
439             dgemm_(&trans,&notr,&dimRU,&dimLD,&dimLU,&alpha,BlockTup,&dimLU,mem,&dimLU,&beta,workmem2,&dimRU);
440 
441             alpha = 1.0;
442             beta = 1.0; //add
443             dgemm_(&notr,&notr,&dimRU,&dimRD,&dimLD,&alpha,workmem2,&dimRU,BlockTdo,&dimLD,&beta,storage+kappa2index[ikappa],&dimRU);
444 
445          }
446       }
447 
448       //case 2
449       const int ILD = Irreps::directProd(ID,bk_up->gIrrep(index-1));
450       for (int TwoSLD=sector_spin_down[ikappa]-1; TwoSLD<=sector_spin_down[ikappa]+1; TwoSLD+=2){
451          int dimLU = bk_up->gCurrentDim(index-1, sector_nelec_up[ikappa]-2, sector_spin_up[ikappa], sector_irrep_up[ikappa]);
452          int dimLD = bk_up->gCurrentDim(index-1, sector_nelec_up[ikappa],   TwoSLD,              ILD);
453          if ((dimLU>0) && (dimLD>0)){
454 
455             int fase = ((((sector_spin_up[ikappa] + sector_spin_down[ikappa] + 1)/2)%2)!=0)?-1:1;
456             double factorB = fase * sqrt(3.0*(TwoSLD+1)) * Wigner::wigner6j( 1, 2, 1, sector_spin_up[ikappa], sector_spin_down[ikappa], TwoSLD );
457 
458             double alpha;
459             double * mem;
460 
461             if (TwoSLD == sector_spin_up[ikappa]){
462 
463                double factorA = - sqrt(0.5);
464 
465                double * BlockA = denA->gStorage( sector_nelec_up[ikappa]-2, sector_spin_up[ikappa], sector_irrep_up[ikappa], sector_nelec_up[ikappa], TwoSLD, ILD );
466                double * BlockB = denB->gStorage( sector_nelec_up[ikappa]-2, sector_spin_up[ikappa], sector_irrep_up[ikappa], sector_nelec_up[ikappa], TwoSLD, ILD );
467 
468                mem = workmem;
469                for (int cnt=0; cnt<dimLU*dimLD; cnt++){ mem[cnt] = factorA * BlockA[cnt] + factorB * BlockB[cnt]; }
470                alpha = 1.0;
471 
472             } else {
473 
474                alpha = factorB;
475                mem = denB->gStorage( sector_nelec_up[ikappa]-2, sector_spin_up[ikappa], sector_irrep_up[ikappa], sector_nelec_up[ikappa], TwoSLD, ILD );
476 
477             }
478 
479             double * BlockTup = denT->gStorage(sector_nelec_up[ikappa]-2,sector_spin_up[ikappa],sector_irrep_up[ikappa],sector_nelec_up[ikappa],   sector_spin_up[ikappa],    sector_irrep_up[ikappa]);
480             double * BlockTdo = denT->gStorage(sector_nelec_up[ikappa]  ,TwoSLD,             ILD,             sector_nelec_up[ikappa]+1, sector_spin_down[ikappa], ID);
481 
482             char trans = 'T';
483             char notr = 'N';
484             double beta = 0.0; //set
485             dgemm_(&trans,&notr,&dimRU,&dimLD,&dimLU,&alpha,BlockTup,&dimLU,mem,&dimLU,&beta,workmem2,&dimRU);
486 
487             alpha = 1.0;
488             beta = 1.0; //add
489             dgemm_(&notr,&notr,&dimRU,&dimRD,&dimLD,&alpha,workmem2,&dimRU,BlockTdo,&dimLD,&beta,storage+kappa2index[ikappa],&dimRU);
490 
491          }
492       }
493    }
494 
495 }
496 
AddTermsABLeft(TensorOperator * denA,TensorOperator * denB,TensorT * denT,double * workmem,double * workmem2)497 void CheMPS2::TensorQ::AddTermsABLeft(TensorOperator * denA, TensorOperator * denB, TensorT * denT, double * workmem, double * workmem2){
498 
499    for (int ikappa=0; ikappa<nKappa; ikappa++){
500       const int ID  = Irreps::directProd( sector_irrep_up[ikappa], n_irrep );
501       int dimLU = bk_up->gCurrentDim(index, sector_nelec_up[ikappa],   sector_spin_up[ikappa],    sector_irrep_up[ikappa]);
502       int dimLD = bk_up->gCurrentDim(index, sector_nelec_up[ikappa]+1, sector_spin_down[ikappa], ID);
503 
504       //case 1
505       const int IRD = Irreps::directProd(ID, bk_up->gIrrep(index));
506       for (int TwoSRD=sector_spin_down[ikappa]-1; TwoSRD<=sector_spin_down[ikappa]+1; TwoSRD+=2){
507          int dimRU = bk_up->gCurrentDim(index+1, sector_nelec_up[ikappa],   sector_spin_up[ikappa], sector_irrep_up[ikappa]);
508          int dimRD = bk_up->gCurrentDim(index+1, sector_nelec_up[ikappa]+2, TwoSRD,              IRD);
509          if ((dimRU>0) && (dimRD>0)){
510 
511             int fase = ((((TwoSRD + sector_spin_up[ikappa] + 2)/2)%2)!=0)?-1:1;
512             const double factorB = fase * sqrt(3.0/(sector_spin_down[ikappa]+1.0)) * (TwoSRD+1)
513                                  * Wigner::wigner6j( 1, 1, 2, sector_spin_up[ikappa], TwoSRD, sector_spin_down[ikappa] );
514 
515             double alpha;
516             double * mem;
517 
518             if (TwoSRD == sector_spin_up[ikappa]){
519 
520                fase = ((((sector_spin_up[ikappa]+1-sector_spin_down[ikappa])/2)%2)!=0)?-1:1;
521                double factorA = fase * sqrt( 0.5 * (sector_spin_up[ikappa]+1.0) / (sector_spin_down[ikappa]+1.0) );
522 
523                double * BlockA = denA->gStorage( sector_nelec_up[ikappa], sector_spin_up[ikappa], sector_irrep_up[ikappa], sector_nelec_up[ikappa]+2, TwoSRD, IRD );
524                double * BlockB = denB->gStorage( sector_nelec_up[ikappa], sector_spin_up[ikappa], sector_irrep_up[ikappa], sector_nelec_up[ikappa]+2, TwoSRD, IRD );
525 
526                mem = workmem;
527                for (int cnt=0; cnt<dimRU*dimRD; cnt++){ mem[cnt] = factorA * BlockA[cnt] + factorB * BlockB[cnt]; }
528                alpha = 1.0;
529 
530             } else {
531 
532                alpha = factorB;
533                mem = denB->gStorage( sector_nelec_up[ikappa], sector_spin_up[ikappa], sector_irrep_up[ikappa], sector_nelec_up[ikappa]+2, TwoSRD, IRD );
534 
535             }
536 
537             double * BlockTup = denT->gStorage(sector_nelec_up[ikappa],  sector_spin_up[ikappa],   sector_irrep_up[ikappa],sector_nelec_up[ikappa],  sector_spin_up[ikappa], sector_irrep_up[ikappa]);
538             double * BlockTdo = denT->gStorage(sector_nelec_up[ikappa]+1,sector_spin_down[ikappa],ID,              sector_nelec_up[ikappa]+2,TwoSRD,              IRD);
539 
540             char notr = 'N';
541             double beta = 0.0; //set
542             dgemm_(&notr,&notr,&dimLU,&dimRD,&dimRU,&alpha,BlockTup,&dimLU,mem,&dimRU,&beta,workmem2,&dimLU);
543 
544             alpha = 1.0;
545             beta = 1.0; //add
546             char trans = 'T';
547             dgemm_(&notr,&trans,&dimLU,&dimLD,&dimRD,&alpha,workmem2,&dimLU,BlockTdo,&dimLD,&beta,storage+kappa2index[ikappa],&dimLU);
548 
549          }
550       }
551 
552       //case 2
553       const int IRU = Irreps::directProd(sector_irrep_up[ikappa],bk_up->gIrrep(index));
554       for (int TwoSRU=sector_spin_up[ikappa]-1; TwoSRU<=sector_spin_up[ikappa]+1; TwoSRU+=2){
555          int dimRU = bk_up->gCurrentDim(index+1, sector_nelec_up[ikappa]+1, TwoSRU,                 IRU);
556          int dimRD = bk_up->gCurrentDim(index+1, sector_nelec_up[ikappa]+3, sector_spin_down[ikappa], ID);
557          if ((dimRU>0) && (dimRD>0)){
558 
559             int fase = ((((sector_spin_up[ikappa] + sector_spin_down[ikappa] + 1)/2)%2)!=0)?-1:1;
560             double factorB = fase * sqrt(3.0*(TwoSRU+1)) * Wigner::wigner6j( 1, 1, 2, TwoSRU, sector_spin_down[ikappa], sector_spin_up[ikappa] );
561 
562             double alpha;
563             double * mem;
564 
565             if (TwoSRU == sector_spin_down[ikappa]){
566 
567                double factorA = - sqrt(0.5);
568 
569                double * BlockA = denA->gStorage( sector_nelec_up[ikappa]+1, TwoSRU, IRU, sector_nelec_up[ikappa]+3, sector_spin_down[ikappa], ID );
570                double * BlockB = denB->gStorage( sector_nelec_up[ikappa]+1, TwoSRU, IRU, sector_nelec_up[ikappa]+3, sector_spin_down[ikappa], ID );
571 
572                mem = workmem;
573                for (int cnt=0; cnt<dimRU*dimRD; cnt++){ mem[cnt] = factorA * BlockA[cnt] + factorB * BlockB[cnt]; }
574                alpha = 1.0;
575 
576             } else {
577 
578                alpha = factorB;
579                mem = denB->gStorage( sector_nelec_up[ikappa]+1, TwoSRU, IRU, sector_nelec_up[ikappa]+3, sector_spin_down[ikappa], ID );
580 
581             }
582 
583             double * BlockTup = denT->gStorage(sector_nelec_up[ikappa],   sector_spin_up[ikappa],    sector_irrep_up[ikappa], sector_nelec_up[ikappa]+1, TwoSRU,                 IRU);
584             double * BlockTdo = denT->gStorage(sector_nelec_up[ikappa]+1, sector_spin_down[ikappa], ID,               sector_nelec_up[ikappa]+3, sector_spin_down[ikappa], ID);
585 
586             char notr = 'N';
587             double beta = 0.0; //set
588             dgemm_(&notr,&notr,&dimLU,&dimRD,&dimRU,&alpha,BlockTup,&dimLU,mem,&dimRU,&beta,workmem2,&dimLU);
589 
590             alpha = 1.0;
591             beta = 1.0; //add
592             char trans = 'T';
593             dgemm_(&notr,&trans,&dimLU,&dimLD,&dimRD,&alpha,workmem2,&dimLU,BlockTdo,&dimLD,&beta,storage+kappa2index[ikappa],&dimLU);
594 
595          }
596       }
597    }
598 
599 }
600 
AddTermsCD(TensorOperator * denC,TensorOperator * denD,TensorT * denT,double * workmem,double * workmem2)601 void CheMPS2::TensorQ::AddTermsCD(TensorOperator * denC, TensorOperator * denD, TensorT * denT, double * workmem, double * workmem2){
602 
603    if (moving_right){ AddTermsCDRight(denC, denD, denT, workmem, workmem2); }
604    else{ AddTermsCDLeft(denC, denD, denT, workmem, workmem2); }
605 
606 }
607 
AddTermsCDRight(TensorOperator * denC,TensorOperator * denD,TensorT * denT,double * workmem,double * workmem2)608 void CheMPS2::TensorQ::AddTermsCDRight(TensorOperator * denC, TensorOperator * denD, TensorT * denT, double * workmem, double * workmem2){
609 
610    for (int ikappa=0; ikappa<nKappa; ikappa++){
611       const int IRD  = Irreps::directProd( sector_irrep_up[ikappa], n_irrep );
612       int dimRU = bk_up->gCurrentDim(index, sector_nelec_up[ikappa],   sector_spin_up[ikappa],    sector_irrep_up[ikappa]);
613       int dimRD = bk_up->gCurrentDim(index, sector_nelec_up[ikappa]+1, sector_spin_down[ikappa], IRD);
614 
615       //case 1
616       const int ILD = Irreps::directProd(IRD,bk_up->gIrrep(index-1));
617       for (int TwoSLD=sector_spin_down[ikappa]-1; TwoSLD<=sector_spin_down[ikappa]+1; TwoSLD+=2){
618          int dimLU = bk_up->gCurrentDim(index-1, sector_nelec_up[ikappa], sector_spin_up[ikappa], sector_irrep_up[ikappa]);
619          int dimLD = bk_up->gCurrentDim(index-1, sector_nelec_up[ikappa], TwoSLD,              ILD);
620 
621          if ((dimLU>0) && (dimLD>0)){
622 
623             int dimLUxLD = dimLU * dimLD;
624 
625             //first set to D
626             int fase = ((((sector_spin_up[ikappa]+sector_spin_down[ikappa]+1)/2)%2)!=0)?-1:1;
627             double factor = fase * sqrt(3.0*(TwoSLD+1)) * Wigner::wigner6j( 1, 2, 1, sector_spin_up[ikappa], sector_spin_down[ikappa], TwoSLD );
628             double * block = denD->gStorage( sector_nelec_up[ikappa], sector_spin_up[ikappa], sector_irrep_up[ikappa], sector_nelec_up[ikappa], TwoSLD, ILD );
629             for (int cnt=0; cnt<dimLUxLD; cnt++){ workmem[cnt] = factor * block[cnt]; }
630 
631             //add C
632             if (TwoSLD==sector_spin_up[ikappa]){
633                factor = sqrt(0.5);
634                block = denC->gStorage( sector_nelec_up[ikappa], sector_spin_up[ikappa], sector_irrep_up[ikappa], sector_nelec_up[ikappa], TwoSLD, ILD );
635                int inc = 1;
636                daxpy_(&dimLUxLD, &factor, block, &inc, workmem, &inc);
637             }
638 
639             double * BlockTup = denT->gStorage( sector_nelec_up[ikappa], sector_spin_up[ikappa], sector_irrep_up[ikappa], sector_nelec_up[ikappa], sector_spin_up[ikappa], sector_irrep_up[ikappa] );
640             double * BlockTdo = denT->gStorage( sector_nelec_up[ikappa], TwoSLD, ILD, sector_nelec_up[ikappa]+1, sector_spin_down[ikappa], IRD );
641 
642             char trans = 'T';
643             char notr = 'N';
644             double alpha = 1.0;
645             double beta = 0.0;
646             dgemm_(&trans,&notr,&dimRU,&dimLD,&dimLU,&alpha,BlockTup,&dimLU,workmem,&dimLU,&beta,workmem2,&dimRU);
647             beta = 1.0;
648             dgemm_(&notr,&notr,&dimRU,&dimRD,&dimLD,&alpha,workmem2,&dimRU,BlockTdo,&dimLD,&beta,storage+kappa2index[ikappa],&dimRU);
649 
650          }
651       }
652 
653       //case 2
654       const int ILU = Irreps::directProd(sector_irrep_up[ikappa],bk_up->gIrrep(index-1));
655       for (int TwoSLU=sector_spin_up[ikappa]-1; TwoSLU<=sector_spin_up[ikappa]+1; TwoSLU+=2){
656          int dimLU = bk_up->gCurrentDim(index-1, sector_nelec_up[ikappa]-1, TwoSLU,                 ILU);
657          int dimLD = bk_up->gCurrentDim(index-1, sector_nelec_up[ikappa]-1, sector_spin_down[ikappa], IRD);
658 
659          if ((dimLU>0) && (dimLD>0)){
660 
661             int dimLUxLD = dimLU * dimLD;
662 
663             //first set to D
664             int fase = ((((TwoSLU + sector_spin_down[ikappa])/2)%2)!=0)?-1:1;
665             double factor = fase * sqrt(3.0*(sector_spin_up[ikappa]+1)) * Wigner::wigner6j( 1, 2, 1, sector_spin_down[ikappa], sector_spin_up[ikappa], TwoSLU );
666             double * block = denD->gStorage( sector_nelec_up[ikappa]-1, TwoSLU, ILU, sector_nelec_up[ikappa]-1, sector_spin_down[ikappa], IRD );
667             for (int cnt=0; cnt<dimLUxLD; cnt++){ workmem[cnt] = factor * block[cnt]; }
668 
669             //add C
670             if (TwoSLU==sector_spin_down[ikappa]){
671                fase = ((((sector_spin_down[ikappa]+1-sector_spin_up[ikappa])/2)%2)!=0)?-1:1;
672                factor = fase * sqrt(0.5 * ( sector_spin_up[ikappa]+1.0 ) / ( sector_spin_down[ikappa] + 1.0 ) );
673                block = denC->gStorage( sector_nelec_up[ikappa]-1, TwoSLU, ILU, sector_nelec_up[ikappa]-1, sector_spin_down[ikappa], IRD );
674                int inc = 1;
675                daxpy_(&dimLUxLD, &factor, block, &inc, workmem, &inc);
676             }
677 
678             double * BlockTup = denT->gStorage( sector_nelec_up[ikappa]-1, TwoSLU,                 ILU, sector_nelec_up[ikappa],   sector_spin_up[ikappa],    sector_irrep_up[ikappa] );
679             double * BlockTdo = denT->gStorage( sector_nelec_up[ikappa]-1, sector_spin_down[ikappa], IRD, sector_nelec_up[ikappa]+1, sector_spin_down[ikappa], IRD );
680 
681             char trans = 'T';
682             char notr = 'N';
683             double alpha = 1.0;
684             double beta = 0.0;
685             dgemm_(&trans,&notr,&dimRU,&dimLD,&dimLU,&alpha,BlockTup,&dimLU,workmem,&dimLU,&beta,workmem2,&dimRU);
686             beta = 1.0;
687             dgemm_(&notr,&notr,&dimRU,&dimRD,&dimLD,&alpha,workmem2,&dimRU,BlockTdo,&dimLD,&beta,storage+kappa2index[ikappa],&dimRU);
688 
689          }
690       }
691    }
692 
693 }
694 
AddTermsCDLeft(TensorOperator * denC,TensorOperator * denD,TensorT * denT,double * workmem,double * workmem2)695 void CheMPS2::TensorQ::AddTermsCDLeft(TensorOperator * denC, TensorOperator * denD, TensorT * denT, double * workmem, double * workmem2){
696 
697    for (int ikappa=0; ikappa<nKappa; ikappa++){
698       const int ILD  = Irreps::directProd( sector_irrep_up[ikappa], n_irrep );
699       int dimLU = bk_up->gCurrentDim(index, sector_nelec_up[ikappa],   sector_spin_up[ikappa],    sector_irrep_up[ikappa]);
700       int dimLD = bk_up->gCurrentDim(index, sector_nelec_up[ikappa]+1, sector_spin_down[ikappa], ILD);
701 
702       //case 1
703       const int IRU = Irreps::directProd(sector_irrep_up[ikappa],bk_up->gIrrep(index));
704       for (int TwoSRU=sector_spin_up[ikappa]-1; TwoSRU<=sector_spin_up[ikappa]+1; TwoSRU+=2){
705          int dimRU = bk_up->gCurrentDim(index+1, sector_nelec_up[ikappa]+1, TwoSRU,                 IRU);
706          int dimRD = bk_up->gCurrentDim(index+1, sector_nelec_up[ikappa]+1, sector_spin_down[ikappa], ILD);
707 
708          if ((dimRU>0) && (dimRD>0)){
709 
710             int dimRUxRD = dimRU * dimRD;
711 
712             //first set to D
713             int fase = ((((sector_spin_up[ikappa]+TwoSRU+3)/2)%2)!=0)?-1:1;
714             double factor = fase * sqrt(3.0/(sector_spin_down[ikappa]+1.0)) * ( TwoSRU + 1 )
715                           * Wigner::wigner6j( 1, 1, 2, TwoSRU, sector_spin_down[ikappa], sector_spin_up[ikappa] );
716             double * block = denD->gStorage( sector_nelec_up[ikappa]+1, TwoSRU, IRU, sector_nelec_up[ikappa]+1, sector_spin_down[ikappa], ILD );
717             for (int cnt=0; cnt<dimRUxRD; cnt++){ workmem[cnt] = factor * block[cnt]; }
718 
719             //add C
720             if (TwoSRU==sector_spin_down[ikappa]){
721                factor = sqrt(0.5);
722                block = denC->gStorage( sector_nelec_up[ikappa]+1, TwoSRU, IRU, sector_nelec_up[ikappa]+1, sector_spin_down[ikappa], ILD );
723                int inc = 1;
724                daxpy_(&dimRUxRD, &factor, block, &inc, workmem, &inc);
725             }
726 
727             double * BlockTup = denT->gStorage( sector_nelec_up[ikappa],   sector_spin_up[ikappa],    sector_irrep_up[ikappa], sector_nelec_up[ikappa]+1, TwoSRU,                 IRU );
728             double * BlockTdo = denT->gStorage( sector_nelec_up[ikappa]+1, sector_spin_down[ikappa], ILD,              sector_nelec_up[ikappa]+1, sector_spin_down[ikappa], ILD );
729 
730             char notr = 'N';
731             double alpha = 1.0;
732             double beta = 0.0;
733             dgemm_(&notr,&notr,&dimLU,&dimRD,&dimRU,&alpha,BlockTup,&dimLU,workmem,&dimRU,&beta,workmem2,&dimLU);
734             beta = 1.0;
735             char trans = 'T';
736             dgemm_(&notr,&trans,&dimLU,&dimLD,&dimRD,&alpha,workmem2,&dimLU,BlockTdo,&dimLD,&beta,storage+kappa2index[ikappa],&dimLU);
737 
738          }
739       }
740 
741       //case 2
742       const int IRD = Irreps::directProd(ILD,bk_up->gIrrep(index));
743       for (int TwoSRD=sector_spin_down[ikappa]-1; TwoSRD<=sector_spin_down[ikappa]+1; TwoSRD+=2){
744          int dimRU = bk_up->gCurrentDim(index+1, sector_nelec_up[ikappa]+2, sector_spin_up[ikappa], sector_irrep_up[ikappa]);
745          int dimRD = bk_up->gCurrentDim(index+1, sector_nelec_up[ikappa]+2, TwoSRD,              IRD);
746 
747          if ((dimRU>0) && (dimRD>0)){
748 
749             int dimRUxRD = dimRU * dimRD;
750 
751             //first set to D
752             int fase = (((TwoSRD+1)%2)!=0)?-1:1;
753             double factor = fase * sqrt(3.0*(TwoSRD+1.0)*(sector_spin_up[ikappa]+1.0)/(sector_spin_down[ikappa]+1.0))
754                           * Wigner::wigner6j( 1, 1, 2, sector_spin_up[ikappa], TwoSRD, sector_spin_down[ikappa] );
755             double * block = denD->gStorage( sector_nelec_up[ikappa]+2, sector_spin_up[ikappa], sector_irrep_up[ikappa], sector_nelec_up[ikappa]+2, TwoSRD, IRD );
756             for (int cnt=0; cnt<dimRUxRD; cnt++){ workmem[cnt] = factor * block[cnt]; }
757 
758             //add C
759             if (TwoSRD==sector_spin_up[ikappa]){
760                fase = ((((sector_spin_up[ikappa]+1-sector_spin_down[ikappa])/2)%2)!=0)?-1:1;
761                factor = fase * sqrt(0.5 * ( sector_spin_up[ikappa]+1.0 ) / ( sector_spin_down[ikappa] + 1.0 ) );
762                block = denC->gStorage( sector_nelec_up[ikappa]+2, sector_spin_up[ikappa], sector_irrep_up[ikappa], sector_nelec_up[ikappa]+2, TwoSRD, IRD );
763                int inc = 1;
764                daxpy_(&dimRUxRD, &factor, block, &inc, workmem, &inc);
765             }
766 
767             double * BlockTup = denT->gStorage( sector_nelec_up[ikappa],  sector_spin_up[ikappa],   sector_irrep_up[ikappa],sector_nelec_up[ikappa]+2,sector_spin_up[ikappa],sector_irrep_up[ikappa] );
768             double * BlockTdo = denT->gStorage( sector_nelec_up[ikappa]+1,sector_spin_down[ikappa],ILD,             sector_nelec_up[ikappa]+2,TwoSRD,             IRD);
769 
770             char notr = 'N';
771             double alpha = 1.0;
772             double beta = 0.0;
773             dgemm_(&notr,&notr,&dimLU,&dimRD,&dimRU,&alpha,BlockTup,&dimLU,workmem,&dimRU,&beta,workmem2,&dimLU);
774             beta = 1.0;
775             char trans = 'T';
776             dgemm_(&notr,&trans,&dimLU,&dimLD,&dimRD,&alpha,workmem2,&dimLU,BlockTdo,&dimLD,&beta,storage+kappa2index[ikappa],&dimLU);
777 
778          }
779       }
780    }
781 
782 }
783 
784 
785 
786