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 "TensorX.h"
24 #include "Lapack.h"
25 #include "Wigner.h"
26
TensorX(const int boundary_index,const bool moving_right,const SyBookkeeper * denBK,const Problem * Prob)27 CheMPS2::TensorX::TensorX(const int boundary_index, const bool moving_right, const SyBookkeeper * denBK, const Problem * Prob) :
28 TensorOperator(boundary_index,
29 0, //two_j
30 0, //n_elec
31 0, //n_irrep
32 moving_right,
33 true, //prime_last (doesn't matter for spin-0 tensors)
34 false, //jw_phase (four 2nd quantized operators)
35 denBK,
36 denBK){
37
38 this->Prob = Prob;
39
40 }
41
~TensorX()42 CheMPS2::TensorX::~TensorX(){ }
43
update(TensorT * denT)44 void CheMPS2::TensorX::update(TensorT * denT){
45
46 if (moving_right){
47 //PARALLEL
48 #pragma omp parallel for schedule(dynamic)
49 for (int ikappa=0; ikappa<nKappa; ikappa++){ makenewRight(ikappa, denT); }
50 } else {
51 //PARALLEL
52 #pragma omp parallel for schedule(dynamic)
53 for (int ikappa=0; ikappa<nKappa; ikappa++){ makenewLeft(ikappa, denT); }
54 }
55
56 }
57
update(TensorT * denT,TensorL ** Ltensors,TensorX * Xtensor,TensorQ * Qtensor,TensorOperator * Atensor,TensorOperator * Ctensor,TensorOperator * Dtensor)58 void CheMPS2::TensorX::update(TensorT * denT, TensorL ** Ltensors, TensorX * Xtensor, TensorQ * Qtensor, TensorOperator * Atensor, TensorOperator * Ctensor, TensorOperator * Dtensor){
59
60 if (moving_right){
61 //PARALLEL
62 #pragma omp parallel
63 {
64
65 const bool doOtherThings = (index>1) ? true : false ;
66 const int dimL = (doOtherThings) ? bk_up->gMaxDimAtBound(index-1) : 0 ;
67 const int dimR = (doOtherThings) ? bk_up->gMaxDimAtBound(index) : 0 ;
68 double * workmemLL = (doOtherThings) ? new double[dimL*dimL] : NULL ;
69 double * workmemLR = (doOtherThings) ? new double[dimL*dimR] : NULL ;
70 double * workmemRR = (doOtherThings) ? new double[dimR*dimR] : NULL ;
71
72 #pragma omp for schedule(dynamic)
73 for (int ikappa=0; ikappa<nKappa; ikappa++){
74 makenewRight(ikappa, denT);
75 if (doOtherThings){
76 update_moving_right(ikappa, Xtensor, denT, denT, workmemLR);
77 addTermQLRight(ikappa, denT, Ltensors, Qtensor, workmemRR, workmemLR, workmemLL);
78 addTermARight(ikappa, denT, Atensor, workmemRR, workmemLR);
79 addTermCRight(ikappa, denT, Ctensor, workmemLR);
80 addTermDRight(ikappa, denT, Dtensor, workmemLR);
81 }
82 }
83
84 if (doOtherThings){
85 delete [] workmemLL;
86 delete [] workmemLR;
87 delete [] workmemRR;
88 }
89
90 }
91 } else {
92 //PARALLEL
93 #pragma omp parallel
94 {
95
96 const bool doOtherThings = (index<Prob->gL()-1) ? true : false ;
97 const int dimL = (doOtherThings) ? bk_up->gMaxDimAtBound(index) : 0 ;
98 const int dimR = (doOtherThings) ? bk_up->gMaxDimAtBound(index+1) : 0 ;
99 double * workmemLL = (doOtherThings) ? new double[dimL*dimL] : NULL ;
100 double * workmemLR = (doOtherThings) ? new double[dimL*dimR] : NULL ;
101 double * workmemRR = (doOtherThings) ? new double[dimR*dimR] : NULL ;
102
103 #pragma omp for schedule(dynamic)
104 for (int ikappa=0; ikappa<nKappa; ikappa++){
105 makenewLeft(ikappa, denT);
106 if (doOtherThings){
107 update_moving_left(ikappa, Xtensor, denT, denT, workmemLR);
108 addTermQLLeft(ikappa, denT, Ltensors, Qtensor, workmemLL, workmemLR, workmemRR);
109 addTermALeft(ikappa, denT, Atensor, workmemLR, workmemLL);
110 addTermCLeft(ikappa, denT, Ctensor, workmemLR);
111 addTermDLeft(ikappa, denT, Dtensor, workmemLR);
112 }
113 }
114
115 if (doOtherThings){
116 delete [] workmemLL;
117 delete [] workmemLR;
118 delete [] workmemRR;
119 }
120
121 }
122 }
123
124 }
125
makenewRight(const int ikappa,TensorT * denT)126 void CheMPS2::TensorX::makenewRight(const int ikappa, TensorT * denT){
127
128 int dimR = bk_up->gCurrentDim(index, sector_nelec_up[ikappa], sector_spin_up[ikappa], sector_irrep_up[ikappa]);
129 int dimL = bk_up->gCurrentDim(index-1, sector_nelec_up[ikappa]-2, sector_spin_up[ikappa], sector_irrep_up[ikappa]);
130 double alpha = Prob->gMxElement(index-1,index-1,index-1,index-1);
131
132 if ((dimL>0) && (fabs(alpha)>0.0)){
133
134 double * BlockT = 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]);
135 char trans = 'T';
136 char notr = 'N';
137 double beta = 0.0; //because there's only 1 term contributing per kappa, we might as well set it i.o. adding
138 dgemm_(&trans,¬r,&dimR,&dimR,&dimL,&alpha,BlockT,&dimL,BlockT,&dimL,&beta,storage+kappa2index[ikappa],&dimR);
139
140 } else {
141 for (int cnt=kappa2index[ikappa]; cnt<kappa2index[ikappa+1]; cnt++){ storage[cnt] = 0.0; }
142 }
143
144 }
145
makenewLeft(const int ikappa,TensorT * denT)146 void CheMPS2::TensorX::makenewLeft(const int ikappa, TensorT * denT){
147
148 int dimL = bk_up->gCurrentDim(index, sector_nelec_up[ikappa], sector_spin_up[ikappa], sector_irrep_up[ikappa]);
149 int dimR = bk_up->gCurrentDim(index+1, sector_nelec_up[ikappa]+2, sector_spin_up[ikappa], sector_irrep_up[ikappa]);
150 double alpha = Prob->gMxElement(index,index,index,index);
151
152 if ((dimR>0) && (fabs(alpha)>0.0)){
153
154 double * BlockT = 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]);
155 char trans = 'T';
156 char notr = 'N';
157 double beta = 0.0; //set, not add (only 1 term)
158 dgemm_(¬r,&trans,&dimL,&dimL,&dimR,&alpha,BlockT,&dimL,BlockT,&dimL,&beta,storage+kappa2index[ikappa],&dimL);
159
160 } else {
161 for (int cnt=kappa2index[ikappa]; cnt<kappa2index[ikappa+1]; cnt++){ storage[cnt] = 0.0; }
162 }
163
164 }
165
166
addTermQLRight(const int ikappa,TensorT * denT,TensorL ** Lprev,TensorQ * Qprev,double * workmemRR,double * workmemLR,double * workmemLL)167 void CheMPS2::TensorX::addTermQLRight(const int ikappa, TensorT * denT, TensorL ** Lprev, TensorQ * Qprev, double * workmemRR, double * workmemLR, double * workmemLL){
168
169 int dimR = bk_up->gCurrentDim(index, sector_nelec_up[ikappa], sector_spin_up[ikappa], sector_irrep_up[ikappa]);
170 int dimTot = dimR * dimR;
171 for (int cnt=0; cnt<dimTot; cnt++){ workmemRR[cnt] = 0.0; }
172
173 for (int geval=0; geval<4; geval++){
174 int NLup,TwoSLup,ILup,NLdown,TwoSLdown,ILdown;
175 switch(geval){
176 case 0:
177 NLup = sector_nelec_up[ikappa]-1;
178 TwoSLup = sector_spin_up[ikappa]-1;
179 ILup = Irreps::directProd( sector_irrep_up[ikappa] , bk_up->gIrrep(index-1) );
180 NLdown = sector_nelec_up[ikappa];
181 TwoSLdown = sector_spin_up[ikappa];
182 ILdown = sector_irrep_up[ikappa];
183 break;
184 case 1:
185 NLup = sector_nelec_up[ikappa]-1;
186 TwoSLup = sector_spin_up[ikappa]+1;
187 ILup = Irreps::directProd( sector_irrep_up[ikappa] , bk_up->gIrrep(index-1) );
188 NLdown = sector_nelec_up[ikappa];
189 TwoSLdown = sector_spin_up[ikappa];
190 ILdown = sector_irrep_up[ikappa];
191 break;
192 case 2:
193 NLup = sector_nelec_up[ikappa]-2;
194 TwoSLup = sector_spin_up[ikappa];
195 ILup = sector_irrep_up[ikappa];
196 NLdown = sector_nelec_up[ikappa]-1;
197 TwoSLdown = sector_spin_up[ikappa]-1;
198 ILdown = Irreps::directProd( sector_irrep_up[ikappa] , bk_up->gIrrep(index-1) );
199 break;
200 case 3:
201 NLup = sector_nelec_up[ikappa]-2;
202 TwoSLup = sector_spin_up[ikappa];
203 ILup = sector_irrep_up[ikappa];
204 NLdown = sector_nelec_up[ikappa]-1;
205 TwoSLdown = sector_spin_up[ikappa]+1;
206 ILdown = Irreps::directProd( sector_irrep_up[ikappa] , bk_up->gIrrep(index-1) );
207 break;
208 }
209 int dimLup = bk_up->gCurrentDim(index-1, NLup, TwoSLup, ILup);
210 int dimLdown = bk_up->gCurrentDim(index-1, NLdown, TwoSLdown, ILdown);
211
212 if ((dimLup>0) && (dimLdown>0)){
213 double * BlockTup = denT->gStorage(NLup, TwoSLup, ILup, sector_nelec_up[ikappa], sector_spin_up[ikappa], sector_irrep_up[ikappa]);
214 double * BlockTdown = denT->gStorage(NLdown, TwoSLdown, ILdown, sector_nelec_up[ikappa], sector_spin_up[ikappa], sector_irrep_up[ikappa]);
215 double * BlockQ = Qprev->gStorage(NLup, TwoSLup, ILup, NLdown, TwoSLdown, ILdown);
216
217 double factor;
218 double * ptr;
219 if (geval<2){
220
221 factor = 1.0;
222 ptr = BlockQ;
223
224 } else {
225
226 int fase = ((((sector_spin_up[ikappa]+1-TwoSLdown)/2)%2)!=0)?-1:1;
227 factor = fase * sqrt((TwoSLdown + 1.0)/(sector_spin_up[ikappa] + 1.0));
228
229 int dimLupdown = dimLup * dimLdown;
230 int inc = 1;
231 ptr = workmemLL;
232 dcopy_(&dimLupdown,BlockQ,&inc,ptr,&inc);
233
234 for (int loca=0; loca<index-1; loca++){
235 if (bk_up->gIrrep(index-1) == bk_up->gIrrep(loca)){
236 double alpha = Prob->gMxElement(loca, index-1, index-1, index-1);
237 double * BlockL = Lprev[index-2-loca]->gStorage(NLup, TwoSLup, ILup, NLdown, TwoSLdown, ILdown);
238 daxpy_(&dimLupdown,&alpha,BlockL,&inc,ptr,&inc);
239 }
240 }
241
242 }
243
244 //factor * Tup^T * L --> mem2 //set
245 char trans = 'T';
246 char notr = 'N';
247 double beta = 0.0;
248 dgemm_(&trans,¬r,&dimR,&dimLdown,&dimLup,&factor,BlockTup,&dimLup,ptr,&dimLup,&beta,workmemLR,&dimR);
249
250 //mem2 * Tdown --> mem //add
251 factor = 1.0;
252 beta = 1.0;
253 dgemm_(¬r,¬r,&dimR,&dimR,&dimLdown,&factor,workmemLR,&dimR,BlockTdown,&dimLdown,&beta,workmemRR,&dimR);
254
255 }
256 }
257 //mem + mem^T --> storage
258 for (int irow = 0; irow<dimR; irow++){
259 for (int icol = irow; icol<dimR; icol++){
260 workmemRR[irow + dimR * icol] += workmemRR[icol + dimR * irow];
261 workmemRR[icol + dimR * irow] = workmemRR[irow + dimR * icol];
262 }
263 }
264 int inc = 1;
265 double alpha = 1.0;
266 daxpy_(&dimTot,&alpha,workmemRR,&inc,storage + kappa2index[ikappa],&inc);
267
268 }
269
addTermQLLeft(const int ikappa,TensorT * denT,TensorL ** Lprev,TensorQ * Qprev,double * workmemLL,double * workmemLR,double * workmemRR)270 void CheMPS2::TensorX::addTermQLLeft(const int ikappa, TensorT * denT, TensorL ** Lprev, TensorQ * Qprev, double * workmemLL, double * workmemLR, double * workmemRR){
271
272 int dimL = bk_up->gCurrentDim(index, sector_nelec_up[ikappa], sector_spin_up[ikappa], sector_irrep_up[ikappa]);
273 int dimTot = dimL * dimL;
274 for (int cnt=0; cnt<dimTot; cnt++){ workmemLL[cnt] = 0.0; }
275
276 for (int geval=0; geval<4; geval++){
277 int NRup,TwoSRup,IRup,NRdown,TwoSRdown,IRdown;
278 switch(geval){
279 case 0:
280 NRup = sector_nelec_up[ikappa];
281 TwoSRup = sector_spin_up[ikappa];
282 IRup = sector_irrep_up[ikappa];
283 NRdown = sector_nelec_up[ikappa]+1;
284 TwoSRdown = sector_spin_up[ikappa]-1;
285 IRdown = Irreps::directProd( sector_irrep_up[ikappa] , bk_up->gIrrep(index) );
286 break;
287 case 1:
288 NRup = sector_nelec_up[ikappa];
289 TwoSRup = sector_spin_up[ikappa];
290 IRup = sector_irrep_up[ikappa];
291 NRdown = sector_nelec_up[ikappa]+1;
292 TwoSRdown = sector_spin_up[ikappa]+1;
293 IRdown = Irreps::directProd( sector_irrep_up[ikappa] , bk_up->gIrrep(index) );
294 break;
295 case 2:
296 NRup = sector_nelec_up[ikappa]+1;
297 TwoSRup = sector_spin_up[ikappa]-1;
298 IRup = Irreps::directProd( sector_irrep_up[ikappa] , bk_up->gIrrep(index) );
299 NRdown = sector_nelec_up[ikappa]+2;
300 TwoSRdown = sector_spin_up[ikappa];
301 IRdown = sector_irrep_up[ikappa];
302 break;
303 case 3:
304 NRup = sector_nelec_up[ikappa]+1;
305 TwoSRup = sector_spin_up[ikappa]+1;
306 IRup = Irreps::directProd( sector_irrep_up[ikappa] , bk_up->gIrrep(index) );
307 NRdown = sector_nelec_up[ikappa]+2;
308 TwoSRdown = sector_spin_up[ikappa];
309 IRdown = sector_irrep_up[ikappa];
310 break;
311 }
312 int dimRup = bk_up->gCurrentDim(index+1, NRup, TwoSRup, IRup);
313 int dimRdown = bk_up->gCurrentDim(index+1, NRdown, TwoSRdown, IRdown);
314
315 if ((dimRup>0) && (dimRdown>0)){
316 double * BlockTup = denT->gStorage(sector_nelec_up[ikappa], sector_spin_up[ikappa], sector_irrep_up[ikappa], NRup, TwoSRup, IRup);
317 double * BlockTdown = denT->gStorage(sector_nelec_up[ikappa], sector_spin_up[ikappa], sector_irrep_up[ikappa], NRdown, TwoSRdown, IRdown);
318 double * BlockQ = Qprev->gStorage(NRup, TwoSRup, IRup, NRdown, TwoSRdown, IRdown);
319
320 double factor;
321 double * ptr;
322 if (geval<2){
323
324 factor = (TwoSRdown + 1.0)/(sector_spin_up[ikappa] + 1.0);
325 ptr = BlockQ;
326
327 } else {
328
329 int fase = ((((sector_spin_up[ikappa]+1-TwoSRup)/2)%2)!=0)?-1:1;
330 factor = fase * sqrt((TwoSRup + 1.0)/(sector_spin_up[ikappa] + 1.0));
331
332 int dimRupdown = dimRup * dimRdown;
333 ptr = workmemRR;
334 int inc = 1;
335 dcopy_(&dimRupdown,BlockQ,&inc,ptr,&inc);
336
337 for (int loca=index+1; loca<Prob->gL(); loca++){
338 if (bk_up->gIrrep(index) == bk_up->gIrrep(loca)){
339 double alpha = Prob->gMxElement(index,index,index,loca);
340 double * BlockL = Lprev[loca-index-1]->gStorage(NRup, TwoSRup, IRup, NRdown, TwoSRdown, IRdown);
341 daxpy_(&dimRupdown,&alpha,BlockL,&inc,ptr,&inc);
342 }
343 }
344
345 }
346
347 //factor * Tup * L --> mem2 //set
348 char notr = 'N';
349 double beta = 0.0;//set
350 dgemm_(¬r,¬r,&dimL,&dimRdown,&dimRup,&factor,BlockTup,&dimL,ptr,&dimRup,&beta,workmemLR,&dimL);
351
352 //mem2 * Tdown^T --> mem //add
353 char trans = 'T';
354 factor = 1.0;
355 beta = 1.0;
356 dgemm_(¬r,&trans,&dimL,&dimL,&dimRdown,&factor,workmemLR,&dimL,BlockTdown,&dimL,&beta,workmemLL,&dimL);
357
358 }
359 }
360 //mem + mem^T --> storage
361 for (int irow = 0; irow<dimL; irow++){
362 for (int icol = irow; icol<dimL; icol++){
363 workmemLL[irow + dimL * icol] += workmemLL[icol + dimL * irow];
364 workmemLL[icol + dimL * irow] = workmemLL[irow + dimL * icol];
365 }
366 }
367 int inc = 1;
368 double alpha = 1.0;
369 daxpy_(&dimTot,&alpha,workmemLL,&inc,storage + kappa2index[ikappa],&inc);
370
371 }
372
addTermARight(const int ikappa,TensorT * denT,TensorOperator * Aprev,double * workmemRR,double * workmemLR)373 void CheMPS2::TensorX::addTermARight(const int ikappa, TensorT * denT, TensorOperator * Aprev, double * workmemRR, double * workmemLR){
374
375 int dimR = bk_up->gCurrentDim(index, sector_nelec_up[ikappa], sector_spin_up[ikappa], sector_irrep_up[ikappa]);
376 int dimLup = bk_up->gCurrentDim(index-1, sector_nelec_up[ikappa]-2, sector_spin_up[ikappa], sector_irrep_up[ikappa]);
377 int dimLdown = bk_up->gCurrentDim(index-1, sector_nelec_up[ikappa], sector_spin_up[ikappa], sector_irrep_up[ikappa]);
378
379 if ((dimLup>0) && (dimLdown>0)){
380
381 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]);
382 double * BlockTdown = 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]);
383 double * BlockA = Aprev->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]);
384
385 //factor * Tup^T * A --> mem2 //set
386 char trans = 'T';
387 char notr = 'N';
388 double factor = sqrt(2.0);
389 double beta = 0.0; //set
390 dgemm_(&trans,¬r,&dimR,&dimLdown,&dimLup,&factor,BlockTup,&dimLup,BlockA,&dimLup,&beta,workmemLR,&dimR);
391
392 //mem2 * Tdown --> mem //set
393 factor = 1.0;
394 dgemm_(¬r,¬r,&dimR,&dimR,&dimLdown,&factor,workmemLR,&dimR,BlockTdown,&dimLdown,&beta,workmemRR,&dimR);
395
396 //mem + mem^T --> storage
397 for (int irow = 0; irow<dimR; irow++){
398 for (int icol = irow; icol<dimR; icol++){
399 workmemRR[irow + dimR * icol] += workmemRR[icol + dimR * irow];
400 workmemRR[icol + dimR * irow] = workmemRR[irow + dimR * icol];
401 }
402 }
403 int dimTot = dimR * dimR;
404 int inc = 1;
405 double alpha = 1.0;
406 daxpy_(&dimTot,&alpha,workmemRR,&inc,storage + kappa2index[ikappa],&inc);
407
408 }
409
410 }
411
addTermALeft(const int ikappa,TensorT * denT,TensorOperator * Aprev,double * workmemLR,double * workmemLL)412 void CheMPS2::TensorX::addTermALeft(const int ikappa, TensorT * denT, TensorOperator * Aprev, double * workmemLR, double * workmemLL){
413
414 int dimL = bk_up->gCurrentDim(index, sector_nelec_up[ikappa], sector_spin_up[ikappa], sector_irrep_up[ikappa]);
415 int dimRup = bk_up->gCurrentDim(index+1, sector_nelec_up[ikappa], sector_spin_up[ikappa], sector_irrep_up[ikappa]);
416 int dimRdown = bk_up->gCurrentDim(index+1, sector_nelec_up[ikappa]+2, sector_spin_up[ikappa], sector_irrep_up[ikappa]);
417
418 if ((dimRup>0) && (dimRdown>0)){
419
420 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]);
421 double * BlockTdown = 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]);
422 double * BlockA = Aprev->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]);
423
424 //factor * Tup * A --> mem2 //set
425 char notr = 'N';
426 double factor = sqrt(2.0);
427 double beta = 0.0; //set
428 dgemm_(¬r,¬r,&dimL,&dimRdown,&dimRup,&factor,BlockTup,&dimL,BlockA,&dimRup,&beta,workmemLR,&dimL);
429
430 //mem2 * Tdown^T --> mem //set
431 char trans = 'T';
432 factor = 1.0;
433 dgemm_(¬r,&trans,&dimL,&dimL,&dimRdown,&factor,workmemLR,&dimL,BlockTdown,&dimL,&beta,workmemLL,&dimL);
434
435 //mem + mem^T --> storage
436 for (int irow = 0; irow<dimL; irow++){
437 for (int icol = irow; icol<dimL; icol++){
438 workmemLL[irow + dimL * icol] += workmemLL[icol + dimL * irow];
439 workmemLL[icol + dimL * irow] = workmemLL[irow + dimL * icol];
440 }
441 }
442 int dimTot = dimL * dimL;
443 int inc = 1;
444 double alpha = 1.0;
445 daxpy_(&dimTot,&alpha,workmemLL,&inc,storage + kappa2index[ikappa],&inc);
446
447 }
448
449 }
450
addTermCRight(const int ikappa,TensorT * denT,TensorOperator * denC,double * workmemLR)451 void CheMPS2::TensorX::addTermCRight(const int ikappa, TensorT * denT, TensorOperator * denC, double * workmemLR){
452
453 int dimR = bk_up->gCurrentDim(index, sector_nelec_up[ikappa], sector_spin_up[ikappa], sector_irrep_up[ikappa]);
454 for (int geval=0; geval<3; geval++){
455 int NL, TwoSL, IL;
456 switch(geval){
457 case 0:
458 NL = sector_nelec_up[ikappa]-1;
459 TwoSL = sector_spin_up[ikappa]-1;
460 IL = Irreps::directProd( sector_irrep_up[ikappa], bk_up->gIrrep(index-1) );
461 break;
462 case 1:
463 NL = sector_nelec_up[ikappa]-1;
464 TwoSL = sector_spin_up[ikappa]+1;
465 IL = Irreps::directProd( sector_irrep_up[ikappa], bk_up->gIrrep(index-1) );
466 break;
467 case 2:
468 NL = sector_nelec_up[ikappa]-2;
469 TwoSL = sector_spin_up[ikappa];
470 IL = sector_irrep_up[ikappa];
471 break;
472 }
473 int dimL = bk_up->gCurrentDim(index-1,NL,TwoSL,IL);
474 if (dimL>0){
475
476 double * BlockC = denC->gStorage(NL,TwoSL,IL,NL,TwoSL,IL);
477 double * BlockT = denT->gStorage(NL,TwoSL,IL,sector_nelec_up[ikappa],sector_spin_up[ikappa],sector_irrep_up[ikappa]);
478
479 double factor = (geval<2)?sqrt(0.5):sqrt(2.0);
480 double beta = 0.0; //set
481 char totrans = 'T';
482 dgemm_(&totrans, &totrans, &dimR, &dimL, &dimL, &factor, BlockT, &dimL, BlockC, &dimL, &beta, workmemLR, &dimR);
483
484 totrans = 'N';
485 factor = 1.0;
486 beta = 1.0; //add
487 dgemm_(&totrans, &totrans, &dimR, &dimR, &dimL, &factor, workmemLR, &dimR, BlockT, &dimL, &beta, storage+kappa2index[ikappa], &dimR);
488
489 }
490 }
491
492 }
493
addTermCLeft(const int ikappa,TensorT * denT,TensorOperator * denC,double * workmemLR)494 void CheMPS2::TensorX::addTermCLeft(const int ikappa, TensorT * denT, TensorOperator * denC, double * workmemLR){
495
496 int dimL = bk_up->gCurrentDim(index, sector_nelec_up[ikappa], sector_spin_up[ikappa], sector_irrep_up[ikappa]);
497 for (int geval=0; geval<3; geval++){
498 int NR, TwoSR, IR;
499 switch(geval){
500 case 0:
501 NR = sector_nelec_up[ikappa]+1;
502 TwoSR = sector_spin_up[ikappa]-1;
503 IR = Irreps::directProd( sector_irrep_up[ikappa], bk_up->gIrrep(index) );
504 break;
505 case 1:
506 NR = sector_nelec_up[ikappa]+1;
507 TwoSR = sector_spin_up[ikappa]+1;
508 IR = Irreps::directProd( sector_irrep_up[ikappa], bk_up->gIrrep(index) );
509 break;
510 case 2:
511 NR = sector_nelec_up[ikappa]+2;
512 TwoSR = sector_spin_up[ikappa];
513 IR = sector_irrep_up[ikappa];
514 break;
515 }
516 int dimR = bk_up->gCurrentDim(index+1,NR,TwoSR,IR);
517 if (dimR>0){
518
519 double * BlockC = denC->gStorage(NR,TwoSR,IR,NR,TwoSR,IR);
520 double * BlockT = denT->gStorage(sector_nelec_up[ikappa],sector_spin_up[ikappa],sector_irrep_up[ikappa],NR,TwoSR,IR);
521
522 double factor = (geval<2)?(sqrt(0.5)*(TwoSR+1.0)/(sector_spin_up[ikappa]+1.0)):sqrt(2.0);
523 double beta = 0.0; //set
524 char trans = 'T';
525 char notr = 'N';
526 dgemm_(¬r, &trans, &dimL, &dimR, &dimR, &factor, BlockT, &dimL, BlockC, &dimR, &beta, workmemLR, &dimL);
527
528 factor = 1.0;
529 beta = 1.0; //add
530 dgemm_(¬r, &trans, &dimL, &dimL, &dimR, &factor, workmemLR, &dimL, BlockT, &dimL, &beta, storage+kappa2index[ikappa], &dimL);
531
532 }
533 }
534
535 }
536
addTermDRight(const int ikappa,TensorT * denT,TensorOperator * denD,double * workmemLR)537 void CheMPS2::TensorX::addTermDRight(const int ikappa, TensorT * denT, TensorOperator * denD, double * workmemLR){
538
539 int dimR = bk_up->gCurrentDim(index, sector_nelec_up[ikappa], sector_spin_up[ikappa], sector_irrep_up[ikappa]);
540
541 const int IL = Irreps::directProd( sector_irrep_up[ikappa], bk_up->gIrrep(index-1) );
542 const int NL = sector_nelec_up[ikappa]-1;
543
544 for (int geval=0; geval<4; geval++){
545 int TwoSLup, TwoSLdown;
546 switch(geval){
547 case 0:
548 TwoSLup = sector_spin_up[ikappa]-1;
549 TwoSLdown = sector_spin_up[ikappa]-1;
550 break;
551 case 1:
552 TwoSLup = sector_spin_up[ikappa]+1;
553 TwoSLdown = sector_spin_up[ikappa]-1;
554 break;
555 case 2:
556 TwoSLup = sector_spin_up[ikappa]-1;
557 TwoSLdown = sector_spin_up[ikappa]+1;
558 break;
559 case 3:
560 TwoSLup = sector_spin_up[ikappa]+1;
561 TwoSLdown = sector_spin_up[ikappa]+1;
562 break;
563 }
564
565 int dimLup = bk_up->gCurrentDim(index-1,NL,TwoSLup, IL);
566 int dimLdown = bk_up->gCurrentDim(index-1,NL,TwoSLdown,IL);
567
568 if ((dimLup>0) && (dimLdown>0)){
569
570 double * BlockD = denD->gStorage(NL,TwoSLdown,IL,NL,TwoSLup,IL);
571 double * BlockTup = denT->gStorage(NL,TwoSLup, IL,sector_nelec_up[ikappa],sector_spin_up[ikappa],sector_irrep_up[ikappa]);
572 double * BlockTdown = (TwoSLup==TwoSLdown)? BlockTup : denT->gStorage(NL,TwoSLdown,IL,sector_nelec_up[ikappa],sector_spin_up[ikappa],sector_irrep_up[ikappa]);
573
574 int fase = ((((TwoSLdown + sector_spin_up[ikappa] + 1)/2)%2)!=0)?-1:1;
575 double factor = fase * sqrt(3.0 * (TwoSLup+1))
576 * Wigner::wigner6j( 1, 1, 2, TwoSLup, TwoSLdown, sector_spin_up[ikappa] );
577 double beta = 0.0; //set
578 char totrans = 'T';
579 dgemm_(&totrans, &totrans, &dimR, &dimLdown, &dimLup, &factor, BlockTup, &dimLup, BlockD, &dimLdown, &beta, workmemLR, &dimR);
580
581 totrans = 'N';
582 factor = 1.0;
583 beta = 1.0; //add
584 dgemm_(&totrans, &totrans, &dimR, &dimR, &dimLdown, &factor, workmemLR, &dimR, BlockTdown, &dimLdown, &beta, storage+kappa2index[ikappa], &dimR);
585
586 }
587 }
588
589 }
590
addTermDLeft(const int ikappa,TensorT * denT,TensorOperator * denD,double * workmemLR)591 void CheMPS2::TensorX::addTermDLeft(const int ikappa, TensorT * denT, TensorOperator * denD, double * workmemLR){
592
593 int dimL = bk_up->gCurrentDim(index, sector_nelec_up[ikappa], sector_spin_up[ikappa], sector_irrep_up[ikappa]);
594
595 const int NR = sector_nelec_up[ikappa]+1;
596 const int IR = Irreps::directProd( sector_irrep_up[ikappa], bk_up->gIrrep(index) );
597
598 for (int geval=0; geval<4; geval++){
599 int TwoSRup, TwoSRdown;
600 switch(geval){
601 case 0:
602 TwoSRup = sector_spin_up[ikappa] - 1;
603 TwoSRdown = sector_spin_up[ikappa] - 1;
604 break;
605 case 1:
606 TwoSRup = sector_spin_up[ikappa] + 1;
607 TwoSRdown = sector_spin_up[ikappa] - 1;
608 break;
609 case 2:
610 TwoSRup = sector_spin_up[ikappa] - 1;
611 TwoSRdown = sector_spin_up[ikappa] + 1;
612 break;
613 case 3:
614 TwoSRup = sector_spin_up[ikappa] + 1;
615 TwoSRdown = sector_spin_up[ikappa] + 1;
616 break;
617 }
618
619 int dimRup = bk_up->gCurrentDim(index+1,NR,TwoSRup, IR);
620 int dimRdown = bk_up->gCurrentDim(index+1,NR,TwoSRdown,IR);
621
622 if ((dimRup>0) && (dimRdown>0)){
623
624 double * BlockD = denD->gStorage(NR,TwoSRdown,IR,NR,TwoSRup,IR);
625 double * BlockTup = denT->gStorage(sector_nelec_up[ikappa],sector_spin_up[ikappa],sector_irrep_up[ikappa],NR,TwoSRup,IR);
626 double * BlockTdown = (TwoSRup == TwoSRdown)? BlockTup : denT->gStorage(sector_nelec_up[ikappa],sector_spin_up[ikappa],sector_irrep_up[ikappa],NR,TwoSRdown,IR);
627
628 int fase = ((((sector_spin_up[ikappa] + TwoSRdown + 3)/2)%2)!=0)?-1:1;
629 double factor = fase*sqrt(3.0 *(TwoSRup+1))*((TwoSRdown + 1.0)/(sector_spin_up[ikappa]+1.0))
630 * Wigner::wigner6j( 1, 1, 2, TwoSRup, TwoSRdown, sector_spin_up[ikappa] );
631 double beta = 0.0; //set
632 char trans = 'T';
633 char notr = 'N';
634 dgemm_(¬r, &trans, &dimL, &dimRdown, &dimRup, &factor, BlockTup, &dimL, BlockD, &dimRdown, &beta, workmemLR, &dimL);
635
636 factor = 1.0;
637 beta = 1.0; //add
638 dgemm_(¬r, &trans, &dimL, &dimL, &dimRdown, &factor, workmemLR, &dimL, BlockTdown, &dimL, &beta, storage+kappa2index[ikappa], &dimL);
639
640 }
641 }
642
643 }
644
645
646