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 <assert.h>
22 #include <math.h>
23
24 #include "Tensor3RDM.h"
25 #include "Special.h"
26 #include "Lapack.h"
27 #include "Wigner.h"
28
Tensor3RDM(const int boundary,const int two_j1_in,const int two_j2,const int nelec,const int irrep,const bool prime_last,const SyBookkeeper * book)29 CheMPS2::Tensor3RDM::Tensor3RDM(const int boundary, const int two_j1_in, const int two_j2, const int nelec, const int irrep, const bool prime_last, const SyBookkeeper * book):
30 TensorOperator(boundary,
31 two_j2,
32 nelec,
33 irrep,
34 true, // moving_right
35 prime_last,
36 true, // jw_phase
37 book,
38 book){
39
40 two_j1 = two_j1_in;
41
42 }
43
~Tensor3RDM()44 CheMPS2::Tensor3RDM::~Tensor3RDM(){ }
45
get_two_j1() const46 int CheMPS2::Tensor3RDM::get_two_j1() const{ return two_j1; }
47
get_two_j2() const48 int CheMPS2::Tensor3RDM::get_two_j2() const{ return get_2j(); }
49
get_prime_last() const50 bool CheMPS2::Tensor3RDM::get_prime_last() const{ return prime_last; }
51
a1(TensorOperator * Sigma,TensorT * denT,double * workmem)52 void CheMPS2::Tensor3RDM::a1(TensorOperator * Sigma, TensorT * denT, double * workmem){
53
54 clear();
55 assert( two_j1 == Sigma->get_2j() );
56 assert( n_elec == 3 );
57 assert( n_irrep == Irreps::directProd( Sigma->get_irrep(), bk_up->gIrrep( index-1 ) ) );
58 const int two_j2 = two_j;
59
60 for ( int ikappa = 0; ikappa < nKappa; ikappa++ ){
61
62 const int two_jr_up = sector_spin_up[ ikappa ];
63 const int nr_up = sector_nelec_up[ ikappa ];
64 const int ir_up = sector_irrep_up[ ikappa ];
65 const int two_jr_down = sector_spin_down[ ikappa ];
66 const int ir_down = Irreps::directProd( ir_up, n_irrep );
67
68 int dimRup = bk_up->gCurrentDim( index, nr_up, two_jr_up, ir_up );
69 int dimRdown = bk_up->gCurrentDim( index, nr_up+3, two_jr_down, ir_down );
70
71 { // Contribution 1
72 const int il_down = Irreps::directProd( ir_down, bk_up->gIrrep( index-1 ) );
73 for ( int two_jl_down = two_jr_down-1; two_jl_down <= two_jr_down+1; two_jl_down+=2 ){
74
75 int dimLup = bk_up->gCurrentDim( index-1, nr_up, two_jr_up, ir_up );
76 int dimLdown = bk_up->gCurrentDim( index-1, nr_up+2, two_jl_down, il_down );
77
78 if (( dimLup > 0 ) && ( dimLdown > 0 ) && ( abs( two_jr_up - two_jl_down ) <= two_j1 )){
79
80 double * Sblock = Sigma->gStorage( nr_up, two_jr_up, ir_up, nr_up+2, two_jl_down, il_down );
81 double * Tup = denT->gStorage( nr_up, two_jr_up, ir_up, nr_up, two_jr_up, ir_up );
82 double * Tdown = denT->gStorage( nr_up+2, two_jl_down, il_down, nr_up+3, two_jr_down, ir_down );
83
84 double alpha = sqrt( 1.0 * ( two_j2 + 1 ) * ( two_jl_down + 1 ) )
85 * Wigner::wigner6j( 1, two_j1, two_j2, two_jr_up, two_jr_down, two_jl_down )
86 * Special::phase( two_jr_up + two_jr_down + two_j1 + 1 );
87 double beta = 0.0; //set
88 char trans = 'T';
89 char notrans = 'N';
90 dgemm_( ¬rans, ¬rans, &dimLup, &dimRdown, &dimLdown, &alpha, Sblock, &dimLup, Tdown, &dimLdown, &beta, workmem, &dimLup );
91 alpha = 1.0;
92 beta = 1.0; //add
93 dgemm_( &trans, ¬rans, &dimRup, &dimRdown, &dimLup, &alpha, Tup, &dimLup, workmem, &dimLup, &beta, storage + kappa2index[ikappa], &dimRup );
94
95 }
96 }
97 }
98 { // Contribution 2
99 const int il_up = Irreps::directProd( ir_up, bk_up->gIrrep( index-1 ) );
100 for ( int two_jl_up = two_jr_up-1; two_jl_up <= two_jr_up+1; two_jl_up+=2 ){
101
102 int dimLup = bk_up->gCurrentDim( index-1, nr_up-1, two_jl_up, il_up );
103 int dimLdown = bk_up->gCurrentDim( index-1, nr_up+1, two_jr_down, ir_down );
104
105 if (( dimLup > 0 ) && ( dimLdown > 0 ) && ( abs( two_jr_down - two_jl_up ) <= two_j1 )){
106
107 double * Sblock = Sigma->gStorage( nr_up-1, two_jl_up, il_up, nr_up+1, two_jr_down, ir_down );
108 double * Tup = denT->gStorage( nr_up-1, two_jl_up, il_up, nr_up, two_jr_up, ir_up );
109 double * Tdown = denT->gStorage( nr_up+1, two_jr_down, ir_down, nr_up+3, two_jr_down, ir_down );
110
111 double alpha = sqrt( 1.0 * ( two_j2 + 1 ) * ( two_jr_up + 1 ) )
112 * Wigner::wigner6j( 1, two_j1, two_j2, two_jr_down, two_jr_up, two_jl_up )
113 * Special::phase( two_jl_up + two_jr_down + two_j2 + 1 );
114 double beta = 0.0; //set
115 char trans = 'T';
116 char notrans = 'N';
117 dgemm_( ¬rans, ¬rans, &dimLup, &dimRdown, &dimLdown, &alpha, Sblock, &dimLup, Tdown, &dimLdown, &beta, workmem, &dimLup );
118 alpha = 1.0;
119 beta = 1.0; //add
120 dgemm_( &trans, ¬rans, &dimRup, &dimRdown, &dimLup, &alpha, Tup, &dimLup, workmem, &dimLup, &beta, storage + kappa2index[ikappa], &dimRup );
121
122 }
123 }
124 }
125 }
126 }
127
b1(TensorOperator * Sigma,TensorT * denT,double * workmem)128 void CheMPS2::Tensor3RDM::b1(TensorOperator * Sigma, TensorT * denT, double * workmem){
129
130 clear();
131 assert( two_j1 == Sigma->get_2j() );
132 assert( n_elec == 1 );
133 assert( n_irrep == Irreps::directProd( Sigma->get_irrep(), bk_up->gIrrep( index-1 ) ) );
134 const int two_j2 = two_j;
135
136 for ( int ikappa = 0; ikappa < nKappa; ikappa++ ){
137
138 const int two_jr_up = sector_spin_up[ ikappa ];
139 const int nr_up = sector_nelec_up[ ikappa ];
140 const int ir_up = sector_irrep_up[ ikappa ];
141 const int two_jr_down = sector_spin_down[ ikappa ];
142 const int ir_down = Irreps::directProd( ir_up, n_irrep );
143
144 int dimRup = bk_up->gCurrentDim( index, nr_up, two_jr_up, ir_up );
145 int dimRdown = bk_up->gCurrentDim( index, nr_up+1, two_jr_down, ir_down );
146
147 { // Contribution 1
148 const int il_up = Irreps::directProd( ir_up, bk_up->gIrrep( index-1 ) );
149 for ( int two_jl_up = two_jr_up-1; two_jl_up <= two_jr_up+1; two_jl_up+=2 ){
150
151 int dimLup = bk_up->gCurrentDim( index-1, nr_up-1, two_jl_up, il_up );
152 int dimLdown = bk_up->gCurrentDim( index-1, nr_up+1, two_jr_down, ir_down );
153
154 if (( dimLup > 0 ) && ( dimLdown > 0 ) && ( abs( two_jr_down - two_jl_up ) <= two_j1 )){
155
156 double * Sblock = Sigma->gStorage( nr_up-1, two_jl_up, il_up, nr_up+1, two_jr_down, ir_down );
157 double * Tup = denT->gStorage( nr_up-1, two_jl_up, il_up, nr_up, two_jr_up, ir_up );
158 double * Tdown = denT->gStorage( nr_up+1, two_jr_down, ir_down, nr_up+1, two_jr_down, ir_down );
159
160 double alpha = sqrt( 1.0 * ( two_j2 + 1 ) * ( two_jr_up + 1 ) )
161 * Wigner::wigner6j( 1, two_j1, two_j2, two_jr_down, two_jr_up, two_jl_up )
162 * Special::phase( two_jl_up + two_jr_down + two_j2 + 3 );
163 double beta = 0.0; //set
164 char trans = 'T';
165 char notrans = 'N';
166 dgemm_( ¬rans, ¬rans, &dimLup, &dimRdown, &dimLdown, &alpha, Sblock, &dimLup, Tdown, &dimLdown, &beta, workmem, &dimLup );
167 alpha = 1.0;
168 beta = 1.0; //add
169 dgemm_( &trans, ¬rans, &dimRup, &dimRdown, &dimLup, &alpha, Tup, &dimLup, workmem, &dimLup, &beta, storage + kappa2index[ikappa], &dimRup );
170
171 }
172 }
173 }
174 { // Contribution 2
175 const int il_down = Irreps::directProd( ir_down, bk_up->gIrrep( index-1 ) );
176 for ( int two_jl_down = two_jr_down-1; two_jl_down <= two_jr_down+1; two_jl_down+=2 ){
177
178 int dimLup = bk_up->gCurrentDim( index-1, nr_up-2, two_jr_up, ir_up );
179 int dimLdown = bk_up->gCurrentDim( index-1, nr_up, two_jl_down, il_down );
180
181 if (( dimLup > 0 ) && ( dimLdown > 0 ) && ( abs( two_jr_up - two_jl_down ) <= two_j1 )){
182
183 double * Sblock = Sigma->gStorage( nr_up-2, two_jr_up, ir_up, nr_up, two_jl_down, il_down );
184 double * Tup = denT->gStorage( nr_up-2, two_jr_up, ir_up, nr_up, two_jr_up, ir_up );
185 double * Tdown = denT->gStorage( nr_up, two_jl_down, il_down, nr_up+1, two_jr_down, ir_down );
186
187 double alpha = sqrt( 1.0 * ( two_j2 + 1 ) * ( two_jl_down + 1 ) )
188 * Wigner::wigner6j( 1, two_j1, two_j2, two_jr_up, two_jr_down, two_jl_down )
189 * Special::phase( two_jr_up + two_jr_down + two_j1 + 1 );
190 double beta = 0.0; //set
191 char trans = 'T';
192 char notrans = 'N';
193 dgemm_( ¬rans, ¬rans, &dimLup, &dimRdown, &dimLdown, &alpha, Sblock, &dimLup, Tdown, &dimLdown, &beta, workmem, &dimLup );
194 alpha = 1.0;
195 beta = 1.0; //add
196 dgemm_( &trans, ¬rans, &dimRup, &dimRdown, &dimLup, &alpha, Tup, &dimLup, workmem, &dimLup, &beta, storage + kappa2index[ikappa], &dimRup );
197
198 }
199 }
200 }
201 }
202 }
203
c1(TensorOperator * denF,TensorT * denT,double * workmem)204 void CheMPS2::Tensor3RDM::c1(TensorOperator * denF, TensorT * denT, double * workmem){
205
206 clear();
207 assert( two_j1 == denF->get_2j() );
208 assert( n_elec == 1 );
209 assert( n_irrep == Irreps::directProd( denF->get_irrep(), bk_up->gIrrep( index-1 ) ) );
210 const int two_j2 = two_j;
211
212 for ( int ikappa = 0; ikappa < nKappa; ikappa++ ){
213
214 const int two_jr_up = sector_spin_up[ ikappa ];
215 const int nr_up = sector_nelec_up[ ikappa ];
216 const int ir_up = sector_irrep_up[ ikappa ];
217 const int two_jr_down = sector_spin_down[ ikappa ];
218 const int ir_down = Irreps::directProd( ir_up, n_irrep );
219
220 int dimRup = bk_up->gCurrentDim( index, nr_up, two_jr_up, ir_up );
221 int dimRdown = bk_up->gCurrentDim( index, nr_up+1, two_jr_down, ir_down );
222
223 { // Contribution 1
224 const int il_down = Irreps::directProd( ir_down, bk_up->gIrrep( index-1 ) );
225 for ( int two_jl_down = two_jr_down-1; two_jl_down <= two_jr_down+1; two_jl_down+=2 ){
226
227 int dimLup = bk_up->gCurrentDim( index-1, nr_up, two_jr_up, ir_up );
228 int dimLdown = bk_up->gCurrentDim( index-1, nr_up, two_jl_down, il_down );
229
230 if (( dimLup > 0 ) && ( dimLdown > 0 ) && ( abs( two_jr_up - two_jl_down ) <= two_j1 )){
231
232 double * Fblock = denF->gStorage( nr_up, two_jr_up, ir_up, nr_up, two_jl_down, il_down );
233 double * Tup = denT->gStorage( nr_up, two_jr_up, ir_up, nr_up, two_jr_up, ir_up );
234 double * Tdown = denT->gStorage( nr_up, two_jl_down, il_down, nr_up+1, two_jr_down, ir_down );
235
236 double alpha = sqrt( 1.0 * ( two_j2 + 1 ) * ( two_jl_down + 1 ) )
237 * Wigner::wigner6j( 1, two_j1, two_j2, two_jr_up, two_jr_down, two_jl_down )
238 * Special::phase( two_jr_up + two_jr_down + two_j1 + 1 );
239 double beta = 0.0; //set
240 char trans = 'T';
241 char notrans = 'N';
242 dgemm_( ¬rans, ¬rans, &dimLup, &dimRdown, &dimLdown, &alpha, Fblock, &dimLup, Tdown, &dimLdown, &beta, workmem, &dimLup );
243 alpha = 1.0;
244 beta = 1.0; //add
245 dgemm_( &trans, ¬rans, &dimRup, &dimRdown, &dimLup, &alpha, Tup, &dimLup, workmem, &dimLup, &beta, storage + kappa2index[ikappa], &dimRup );
246
247 }
248 }
249 }
250 { // Contribution 2
251 const int il_up = Irreps::directProd( ir_up, bk_up->gIrrep( index-1 ) );
252 for ( int two_jl_up = two_jr_up-1; two_jl_up <= two_jr_up+1; two_jl_up+=2 ){
253
254 int dimLup = bk_up->gCurrentDim( index-1, nr_up-1, two_jl_up, il_up );
255 int dimLdown = bk_up->gCurrentDim( index-1, nr_up-1, two_jr_down, ir_down );
256
257 if (( dimLup > 0 ) && ( dimLdown > 0 ) && ( abs( two_jr_down - two_jl_up ) <= two_j1 )){
258
259 double * Fblock = denF->gStorage( nr_up-1, two_jl_up, il_up, nr_up-1, two_jr_down, ir_down );
260 double * Tup = denT->gStorage( nr_up-1, two_jl_up, il_up, nr_up, two_jr_up, ir_up );
261 double * Tdown = denT->gStorage( nr_up-1, two_jr_down, ir_down, nr_up+1, two_jr_down, ir_down );
262
263 double alpha = sqrt( 1.0 * ( two_j2 + 1 ) * ( two_jr_up + 1 ) )
264 * Wigner::wigner6j( 1, two_j1, two_j2, two_jr_down, two_jr_up, two_jl_up )
265 * Special::phase( two_jl_up + two_jr_down + two_j2 + 1 );
266 double beta = 0.0; //set
267 char trans = 'T';
268 char notrans = 'N';
269 dgemm_( ¬rans, ¬rans, &dimLup, &dimRdown, &dimLdown, &alpha, Fblock, &dimLup, Tdown, &dimLdown, &beta, workmem, &dimLup );
270 alpha = 1.0;
271 beta = 1.0; //add
272 dgemm_( &trans, ¬rans, &dimRup, &dimRdown, &dimLup, &alpha, Tup, &dimLup, workmem, &dimLup, &beta, storage + kappa2index[ikappa], &dimRup );
273
274 }
275 }
276 }
277 }
278 }
279
d1(TensorOperator * denF,TensorT * denT,double * workmem)280 void CheMPS2::Tensor3RDM::d1(TensorOperator * denF, TensorT * denT, double * workmem){
281
282 clear();
283 assert( two_j1 == denF->get_2j() );
284 assert( n_elec == 1 );
285 assert( n_irrep == Irreps::directProd( denF->get_irrep(), bk_up->gIrrep( index-1 ) ) );
286 const int two_j2 = two_j;
287
288 for ( int ikappa = 0; ikappa < nKappa; ikappa++ ){
289
290 const int two_jr_up = sector_spin_up[ ikappa ];
291 const int nr_up = sector_nelec_up[ ikappa ];
292 const int ir_up = sector_irrep_up[ ikappa ];
293 const int two_jr_down = sector_spin_down[ ikappa ];
294 const int ir_down = Irreps::directProd( ir_up, n_irrep );
295
296 int dimRup = bk_up->gCurrentDim( index, nr_up, two_jr_up, ir_up );
297 int dimRdown = bk_up->gCurrentDim( index, nr_up+1, two_jr_down, ir_down );
298
299 { // Contribution 1
300 const int il_down = Irreps::directProd( ir_down, bk_up->gIrrep( index-1 ) );
301 for ( int two_jl_down = two_jr_down-1; two_jl_down <= two_jr_down+1; two_jl_down+=2 ){
302
303 int dimLup = bk_up->gCurrentDim( index-1, nr_up, two_jr_up, ir_up );
304 int dimLdown = bk_up->gCurrentDim( index-1, nr_up, two_jl_down, il_down );
305
306 if (( dimLup > 0 ) && ( dimLdown > 0 ) && ( abs( two_jr_up - two_jl_down ) <= two_j1 )){
307
308 double * Fblock = denF->gStorage( nr_up, two_jl_down, il_down, nr_up, two_jr_up, ir_up );
309 double * Tup = denT->gStorage( nr_up, two_jr_up, ir_up, nr_up, two_jr_up, ir_up );
310 double * Tdown = denT->gStorage( nr_up, two_jl_down, il_down, nr_up+1, two_jr_down, ir_down );
311
312 double alpha = sqrt( 1.0 * ( two_j2 + 1 ) * ( two_jr_down + 1 ) )
313 * Wigner::wigner6j( 1, two_j1, two_j2, two_jr_up, two_jr_down, two_jl_down )
314 * Special::phase( two_jr_up + two_jl_down + two_j2 + 3 );
315 double beta = 0.0; //set
316 char trans = 'T';
317 char notrans = 'N';
318 dgemm_( &trans, ¬rans, &dimLup, &dimRdown, &dimLdown, &alpha, Fblock, &dimLdown, Tdown, &dimLdown, &beta, workmem, &dimLup );
319 alpha = 1.0;
320 beta = 1.0; //add
321 dgemm_( &trans, ¬rans, &dimRup, &dimRdown, &dimLup, &alpha, Tup, &dimLup, workmem, &dimLup, &beta, storage + kappa2index[ikappa], &dimRup );
322
323 }
324 }
325 }
326 { // Contribution 2
327 const int il_up = Irreps::directProd( ir_up, bk_up->gIrrep( index-1 ) );
328 for ( int two_jl_up = two_jr_up-1; two_jl_up <= two_jr_up+1; two_jl_up+=2 ){
329
330 int dimLup = bk_up->gCurrentDim( index-1, nr_up-1, two_jl_up, il_up );
331 int dimLdown = bk_up->gCurrentDim( index-1, nr_up-1, two_jr_down, ir_down );
332
333 if (( dimLup > 0 ) && ( dimLdown > 0 ) && ( abs( two_jr_down - two_jl_up ) <= two_j1 )){
334
335 double * Fblock = denF->gStorage( nr_up-1, two_jr_down, ir_down, nr_up-1, two_jl_up, il_up );
336 double * Tup = denT->gStorage( nr_up-1, two_jl_up, il_up, nr_up, two_jr_up, ir_up );
337 double * Tdown = denT->gStorage( nr_up-1, two_jr_down, ir_down, nr_up+1, two_jr_down, ir_down );
338
339 double alpha = sqrt( 1.0 * ( two_j2 + 1 ) * ( two_jl_up + 1 ) )
340 * Wigner::wigner6j( 1, two_j1, two_j2, two_jr_down, two_jr_up, two_jl_up )
341 * Special::phase( two_jr_up + two_jr_down + two_j1 + 1 );
342 double beta = 0.0; //set
343 char trans = 'T';
344 char notrans = 'N';
345 dgemm_( &trans, ¬rans, &dimLup, &dimRdown, &dimLdown, &alpha, Fblock, &dimLdown, Tdown, &dimLdown, &beta, workmem, &dimLup );
346 alpha = 1.0;
347 beta = 1.0; //add
348 dgemm_( &trans, ¬rans, &dimRup, &dimRdown, &dimLup, &alpha, Tup, &dimLup, workmem, &dimLup, &beta, storage + kappa2index[ikappa], &dimRup );
349
350 }
351 }
352 }
353 }
354 }
355
extra1(TensorT * denT)356 void CheMPS2::Tensor3RDM::extra1(TensorT * denT){
357
358 clear();
359 assert( n_elec == 1 );
360 assert( n_irrep == bk_up->gIrrep( index-1 ) );
361 const int two_j2 = two_j;
362 assert( two_j2 == 1 );
363
364 for ( int ikappa = 0; ikappa < nKappa; ikappa++ ){
365
366 const int two_jr_up = sector_spin_up[ ikappa ];
367 const int nr_up = sector_nelec_up[ ikappa ];
368 const int ir_up = sector_irrep_up[ ikappa ];
369 const int two_jr_down = sector_spin_down[ ikappa ];
370 const int ir_down = Irreps::directProd( ir_up, n_irrep );
371
372 int dimRup = bk_up->gCurrentDim( index, nr_up, two_jr_up, ir_up );
373 int dimRdown = bk_up->gCurrentDim( index, nr_up+1, two_jr_down, ir_down );
374 int dimL = bk_up->gCurrentDim( index-1, nr_up-1, two_jr_down, ir_down );
375
376 if ( dimL > 0 ){
377
378 double * Tup = denT->gStorage( nr_up-1, two_jr_down, ir_down, nr_up, two_jr_up, ir_up );
379 double * Tdown = denT->gStorage( nr_up-1, two_jr_down, ir_down, nr_up+1, two_jr_down, ir_down );
380
381 double alpha = sqrt( 0.5 * ( two_j1 + 1 ) ) * Special::phase( two_j1 );
382 double beta = 0.0; //set --> only contribution to this symmetry sector
383 char trans = 'T';
384 char notrans = 'N';
385 dgemm_( &trans, ¬rans, &dimRup, &dimRdown, &dimL, &alpha, Tup, &dimL, Tdown, &dimL, &beta, storage + kappa2index[ikappa], &dimRup );
386
387 }
388 }
389 }
390
extra2(TensorL * denL,TensorT * denT,double * workmem)391 void CheMPS2::Tensor3RDM::extra2(TensorL * denL, TensorT * denT, double * workmem){
392
393 clear();
394 assert( n_elec == 3 );
395 assert( n_irrep == denL->get_irrep() );
396 const int two_j2 = two_j;
397 assert( two_j2 == 1 );
398
399 for ( int ikappa = 0; ikappa < nKappa; ikappa++ ){
400
401 const int two_jr_up = sector_spin_up[ ikappa ];
402 const int nr_up = sector_nelec_up[ ikappa ];
403 const int ir_up = sector_irrep_up[ ikappa ];
404 const int two_jr_down = sector_spin_down[ ikappa ];
405 const int ir_down = Irreps::directProd( ir_up, n_irrep );
406
407 int dimRup = bk_up->gCurrentDim( index, nr_up, two_jr_up, ir_up );
408 int dimRdown = bk_up->gCurrentDim( index, nr_up+3, two_jr_down, ir_down );
409 int dimLup = bk_up->gCurrentDim( index-1, nr_up, two_jr_up, ir_up );
410 int dimLdown = bk_up->gCurrentDim( index-1, nr_up+1, two_jr_down, ir_down );
411
412 if (( dimLup > 0 ) && ( dimLdown > 0 )){
413
414 double * Tup = denT->gStorage( nr_up, two_jr_up, ir_up, nr_up, two_jr_up, ir_up );
415 double * Tdown = denT->gStorage( nr_up+1, two_jr_down, ir_down, nr_up+3, two_jr_down, ir_down );
416 double * Lblock = denL->gStorage( nr_up, two_jr_up, ir_up, nr_up+1, two_jr_down, ir_down );
417
418 double alpha = sqrt( 0.5 * ( two_j1 + 1 ) ) * Special::phase( two_j1 + 2 );
419 double beta = 0.0; //set
420 char trans = 'T';
421 char notrans = 'N';
422 dgemm_( ¬rans, ¬rans, &dimLup, &dimRdown, &dimLdown, &alpha, Lblock, &dimLup, Tdown, &dimLdown, &beta, workmem, &dimLup );
423 alpha = 1.0;
424 dgemm_( &trans, ¬rans, &dimRup, &dimRdown, &dimLup, &alpha, Tup, &dimLup, workmem, &dimLup, &beta, storage + kappa2index[ikappa], &dimRup );
425
426 }
427 }
428 }
429
extra3(TensorL * denL,TensorT * denT,double * workmem)430 void CheMPS2::Tensor3RDM::extra3(TensorL * denL, TensorT * denT, double * workmem){
431
432 clear();
433 assert( n_elec == 1 );
434 assert( n_irrep == denL->get_irrep() );
435 const int two_j2 = two_j;
436 assert( two_j2 == 1 );
437
438 for ( int ikappa = 0; ikappa < nKappa; ikappa++ ){
439
440 const int two_jr_up = sector_spin_up[ ikappa ];
441 const int nr_up = sector_nelec_up[ ikappa ];
442 const int ir_up = sector_irrep_up[ ikappa ];
443 const int two_jr_down = sector_spin_down[ ikappa ];
444 const int ir_down = Irreps::directProd( ir_up, n_irrep );
445
446 int dimRup = bk_up->gCurrentDim( index, nr_up, two_jr_up, ir_up );
447 int dimRdown = bk_up->gCurrentDim( index, nr_up+1, two_jr_down, ir_down );
448 int dimLup = bk_up->gCurrentDim( index-1, nr_up, two_jr_up, ir_up );
449 int dimLdown = bk_up->gCurrentDim( index-1, nr_up-1, two_jr_down, ir_down );
450
451 if (( dimLup > 0 ) && ( dimLdown > 0 )){
452
453 double * Tup = denT->gStorage( nr_up, two_jr_up, ir_up, nr_up, two_jr_up, ir_up );
454 double * Tdown = denT->gStorage( nr_up-1, two_jr_down, ir_down, nr_up+1, two_jr_down, ir_down );
455 double * Lblock = denL->gStorage( nr_up-1, two_jr_down, ir_down, nr_up, two_jr_up, ir_up );
456
457 double alpha = sqrt( 0.5 * ( two_j1 + 1 ) ) * Special::phase( two_j1 );
458 double beta = 0.0; //set
459 char trans = 'T';
460 char notrans = 'N';
461 dgemm_( &trans, ¬rans, &dimLup, &dimRdown, &dimLdown, &alpha, Lblock, &dimLdown, Tdown, &dimLdown, &beta, workmem, &dimLup );
462 alpha = 1.0;
463 dgemm_( &trans, ¬rans, &dimRup, &dimRdown, &dimLup, &alpha, Tup, &dimLup, workmem, &dimLup, &beta, storage + kappa2index[ikappa], &dimRup );
464
465 }
466 }
467 }
468
extra4(TensorL * denL,TensorT * denT,double * workmem)469 void CheMPS2::Tensor3RDM::extra4(TensorL * denL, TensorT * denT, double * workmem){
470
471 clear();
472 assert( n_elec == 1 );
473 assert( n_irrep == denL->get_irrep() );
474 const int two_j2 = two_j;
475
476 for ( int ikappa = 0; ikappa < nKappa; ikappa++ ){
477
478 const int two_jr_up = sector_spin_up[ ikappa ];
479 const int nr_up = sector_nelec_up[ ikappa ];
480 const int ir_up = sector_irrep_up[ ikappa ];
481 const int two_jr_down = sector_spin_down[ ikappa ];
482 const int ir_down = Irreps::directProd( ir_up, n_irrep );
483
484 int dimRup = bk_up->gCurrentDim( index, nr_up, two_jr_up, ir_up );
485 int dimRdown = bk_up->gCurrentDim( index, nr_up+1, two_jr_down, ir_down );
486
487 if ( two_j2 == 1 ){ // Contribution 2
488
489 int dimLup = bk_up->gCurrentDim( index-1, nr_up-2, two_jr_up, ir_up );
490 int dimLdown = bk_up->gCurrentDim( index-1, nr_up-1, two_jr_down, ir_down );
491
492 if (( dimLup > 0 ) && ( dimLdown > 0 )){
493
494 double * Tup = denT->gStorage( nr_up-2, two_jr_up, ir_up, nr_up, two_jr_up, ir_up );
495 double * Tdown = denT->gStorage( nr_up-1, two_jr_down, ir_down, nr_up+1, two_jr_down, ir_down );
496 double * Lblock = denL->gStorage( nr_up-2, two_jr_up, ir_up, nr_up-1, two_jr_down, ir_down );
497
498 double alpha = sqrt( 0.5 * ( two_j1 + 1 ) ) * Special::phase( two_j1 + 2 );
499 double beta = 0.0; //set
500 char trans = 'T';
501 char notrans = 'N';
502 dgemm_( ¬rans, ¬rans, &dimLup, &dimRdown, &dimLdown, &alpha, Lblock, &dimLup, Tdown, &dimLdown, &beta, workmem, &dimLup );
503 alpha = 1.0;
504 beta = 1.0; // add
505 dgemm_( &trans, ¬rans, &dimRup, &dimRdown, &dimLup, &alpha, Tup, &dimLup, workmem, &dimLup, &beta, storage + kappa2index[ikappa], &dimRup );
506
507 }
508 }
509 { // Contribution 1
510 const int il_up = Irreps::directProd( ir_up, bk_up->gIrrep( index-1 ) );
511 const int il_down = Irreps::directProd( ir_down, bk_up->gIrrep( index-1 ) );
512 for ( int two_jl_up = two_jr_up-1; two_jl_up <= two_jr_up+1; two_jl_up+=2 ){
513 for ( int two_jl_down = two_jr_down-1; two_jl_down <= two_jr_down+1; two_jl_down+=2 ){
514
515 int dimLup = bk_up->gCurrentDim( index-1, nr_up-1, two_jl_up, il_up );
516 int dimLdown = bk_up->gCurrentDim( index-1, nr_up, two_jl_down, il_down );
517
518 if (( dimLup > 0 ) && ( dimLdown > 0 ) && ( abs( two_jl_up - two_jl_down ) <= 1 )){
519
520 double * Tup = denT->gStorage( nr_up-1, two_jl_up, il_up, nr_up, two_jr_up, ir_up );
521 double * Tdown = denT->gStorage( nr_up, two_jl_down, il_down, nr_up+1, two_jr_down, ir_down );
522 double * Lblock = denL->gStorage( nr_up-1, two_jl_up, il_up, nr_up, two_jl_down, il_down );
523
524 double alpha = sqrt( 1.0 * ( two_j1 + 1 ) * ( two_j2 + 1 ) * ( two_jr_up + 1 ) * ( two_jl_down + 1 ) )
525 * Special::phase( two_jr_up + two_jr_down - two_jl_up - two_jl_down )
526 * Wigner::wigner6j( 1, 1, two_j1, two_jl_down, two_jr_up, two_jl_up )
527 * Wigner::wigner6j( 1, two_j1, two_j2, two_jr_up, two_jr_down, two_jl_down );
528 double beta = 0.0; //set
529 char trans = 'T';
530 char notrans = 'N';
531 dgemm_( ¬rans, ¬rans, &dimLup, &dimRdown, &dimLdown, &alpha, Lblock, &dimLup, Tdown, &dimLdown, &beta, workmem, &dimLup );
532 alpha = 1.0;
533 beta = 1.0; // add
534 dgemm_( &trans, ¬rans, &dimRup, &dimRdown, &dimLup, &alpha, Tup, &dimLup, workmem, &dimLup, &beta, storage + kappa2index[ikappa], &dimRup );
535
536 }
537 }
538 }
539 }
540 }
541 }
542
contract(Tensor3RDM * buddy) const543 double CheMPS2::Tensor3RDM::contract( Tensor3RDM * buddy ) const{
544
545 if ( buddy == NULL ){ return 0.0; }
546
547 assert( get_two_j2() == buddy->get_two_j2() );
548 assert( n_elec == buddy->get_nelec() );
549 assert( n_irrep == buddy->get_irrep() );
550
551 double value = 0.0;
552
553 if ( buddy->get_prime_last() ){ // Tensor abc
554
555 int length = kappa2index[ nKappa ];
556 int inc = 1;
557 value = ddot_( &length, storage, &inc, buddy->gStorage(), &inc );
558 return value;
559
560 } else { // Tensor d
561
562 for ( int ikappa = 0; ikappa < nKappa; ikappa++ ){
563
564 int offset = kappa2index[ ikappa ];
565 int length = kappa2index[ ikappa + 1 ] - offset;
566 int inc = 1;
567 double prefactor = sqrt( ( sector_spin_up[ ikappa ] + 1.0 ) / ( sector_spin_down[ ikappa ] + 1 ) )
568 * Special::phase( sector_spin_up[ ikappa ] + 1 - sector_spin_down[ ikappa ] );
569 value += prefactor * ddot_( &length, storage + offset, &inc, buddy->gStorage() + offset, &inc );
570
571 }
572
573 return value;
574
575 }
576
577 return value;
578
579 }
580
581
582