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_( &notrans, &notrans, &dimLup, &dimRdown, &dimLdown, &alpha, Sblock, &dimLup, Tdown, &dimLdown, &beta, workmem, &dimLup );
91                alpha = 1.0;
92                beta  = 1.0; //add
93                dgemm_( &trans, &notrans, &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_( &notrans, &notrans, &dimLup, &dimRdown, &dimLdown, &alpha, Sblock, &dimLup, Tdown, &dimLdown, &beta, workmem, &dimLup );
118                alpha = 1.0;
119                beta  = 1.0; //add
120                dgemm_( &trans, &notrans, &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_( &notrans, &notrans, &dimLup, &dimRdown, &dimLdown, &alpha, Sblock, &dimLup, Tdown, &dimLdown, &beta, workmem, &dimLup );
167                alpha = 1.0;
168                beta  = 1.0; //add
169                dgemm_( &trans, &notrans, &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_( &notrans, &notrans, &dimLup, &dimRdown, &dimLdown, &alpha, Sblock, &dimLup, Tdown, &dimLdown, &beta, workmem, &dimLup );
194                alpha = 1.0;
195                beta  = 1.0; //add
196                dgemm_( &trans, &notrans, &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_( &notrans, &notrans, &dimLup, &dimRdown, &dimLdown, &alpha, Fblock, &dimLup, Tdown, &dimLdown, &beta, workmem, &dimLup );
243                alpha = 1.0;
244                beta  = 1.0; //add
245                dgemm_( &trans, &notrans, &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_( &notrans, &notrans, &dimLup, &dimRdown, &dimLdown, &alpha, Fblock, &dimLup, Tdown, &dimLdown, &beta, workmem, &dimLup );
270                alpha = 1.0;
271                beta  = 1.0; //add
272                dgemm_( &trans, &notrans, &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, &notrans, &dimLup, &dimRdown, &dimLdown, &alpha, Fblock, &dimLdown, Tdown, &dimLdown, &beta, workmem, &dimLup );
319                alpha = 1.0;
320                beta  = 1.0; //add
321                dgemm_( &trans, &notrans, &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, &notrans, &dimLup, &dimRdown, &dimLdown, &alpha, Fblock, &dimLdown, Tdown, &dimLdown, &beta, workmem, &dimLup );
346                alpha = 1.0;
347                beta  = 1.0; //add
348                dgemm_( &trans, &notrans, &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, &notrans, &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_( &notrans, &notrans, &dimLup, &dimRdown, &dimLdown, &alpha, Lblock, &dimLup, Tdown, &dimLdown, &beta, workmem, &dimLup );
423          alpha = 1.0;
424          dgemm_( &trans, &notrans, &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, &notrans, &dimLup, &dimRdown, &dimLdown, &alpha, Lblock, &dimLdown, Tdown, &dimLdown, &beta, workmem, &dimLup );
462          alpha = 1.0;
463          dgemm_( &trans, &notrans, &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_( &notrans, &notrans, &dimLup, &dimRdown, &dimLdown, &alpha, Lblock, &dimLup, Tdown, &dimLdown, &beta, workmem, &dimLup );
503             alpha = 1.0;
504             beta  = 1.0; // add
505             dgemm_( &trans, &notrans, &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_( &notrans, &notrans, &dimLup, &dimRdown, &dimLdown, &alpha, Lblock, &dimLup, Tdown, &dimLdown, &beta, workmem, &dimLup );
532                   alpha = 1.0;
533                   beta  = 1.0; // add
534                   dgemm_( &trans, &notrans, &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