CheMPS2
Tensor3RDM.cpp
1 /*
2  CheMPS2: a spin-adapted implementation of DMRG for ab initio quantum chemistry
3  Copyright (C) 2013-2016 Sebastian Wouters
4 
5  This program is free software; you can redistribute it and/or modify
6  it under the terms of the GNU General Public License as published by
7  the Free Software Foundation; either version 2 of the License, or
8  (at your option) any later version.
9 
10  This program is distributed in the hope that it will be useful,
11  but WITHOUT ANY WARRANTY; without even the implied warranty of
12  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13  GNU General Public License for more details.
14 
15  You should have received a copy of the GNU General Public License along
16  with this program; if not, write to the Free Software Foundation, Inc.,
17  51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
18 */
19 
20 #include <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 
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 
44 CheMPS2::Tensor3RDM::~Tensor3RDM(){ }
45 
46 int CheMPS2::Tensor3RDM::get_two_j1() const{ return two_j1; }
47 
48 int CheMPS2::Tensor3RDM::get_two_j2() const{ return get_2j(); }
49 
50 bool CheMPS2::Tensor3RDM::get_prime_last() const{ return prime_last; }
51 
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 
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 
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 
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 
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 
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 
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 
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 
543 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 
CheMPS2::Wigner::wigner6j
static double wigner6j(const int two_ja, const int two_jb, const int two_jc, const int two_jd, const int two_je, const int two_jf)
Wigner-6j symbol (gsl api)
Definition: Wigner.cpp:294
CheMPS2::Tensor::index
int index
Index of the Tensor object. For TensorT: a site index; for other tensors: a boundary index...
Definition: Tensor.h:75
CheMPS2::TensorOperator::sector_nelec_up
int * sector_nelec_up
The up particle number sector.
Definition: TensorOperator.h:157
CheMPS2::TensorOperator::get_2j
int get_2j() const
Get twice the spin of the tensor operator.
Definition: TensorOperator.cpp:151
CheMPS2::Tensor3RDM::extra1
void extra1(TensorT *denT)
Make diagram extra1.
Definition: Tensor3RDM.cpp:356
CheMPS2::TensorOperator::sector_spin_down
int * sector_spin_down
The down spin symmetry sector (pointer points to sectorTwoS1 if two_j == 0)
Definition: TensorOperator.h:166
CheMPS2::SyBookkeeper::gCurrentDim
int gCurrentDim(const int boundary, const int N, const int TwoS, const int irrep) const
Get the current virtual dimensions.
Definition: SyBookkeeper.cpp:165
CheMPS2::TensorOperator::get_nelec
int get_nelec() const
Get how many electrons there are more in the symmetry sector of the lower leg compared to the upper l...
Definition: TensorOperator.cpp:153
CheMPS2::Tensor3RDM::c1
void c1(TensorOperator *denF, TensorT *denT, double *workmem)
Make diagram c1.
Definition: Tensor3RDM.cpp:204
CheMPS2::TensorOperator::prime_last
bool prime_last
Convention in which the tensor operator is stored (see class information)
Definition: TensorOperator.h:185
CheMPS2::Special::phase
static int phase(const int power_times_two)
Phase function.
Definition: Special.h:36
CheMPS2::Tensor3RDM::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)
Constructor.
Definition: Tensor3RDM.cpp:29
CheMPS2::Tensor3RDM::get_prime_last
bool get_prime_last() const
Get whether the tensor convention is prime last.
Definition: Tensor3RDM.cpp:50
CheMPS2::Tensor3RDM::contract
double contract(Tensor3RDM *buddy) const
Make the in-product of two Tensor3RDMs.
Definition: Tensor3RDM.cpp:543
CheMPS2::Tensor3RDM::extra3
void extra3(TensorL *denL, TensorT *denT, double *workmem)
Make diagram extra3.
Definition: Tensor3RDM.cpp:430
CheMPS2::Tensor3RDM::a1
void a1(TensorOperator *Sigma, TensorT *denT, double *workmem)
Make diagram a1.
Definition: Tensor3RDM.cpp:52
CheMPS2::Tensor3RDM::extra2
void extra2(TensorL *denL, TensorT *denT, double *workmem)
Make diagram extra2.
Definition: Tensor3RDM.cpp:391
CheMPS2::TensorOperator::gStorage
double * gStorage()
Get the pointer to the storage.
Definition: TensorOperator.cpp:117
CheMPS2::Tensor3RDM::get_two_j1
int get_two_j1() const
Get the intermediary spin coupling value of the Tensor3RDM.
Definition: Tensor3RDM.cpp:46
CheMPS2::Tensor3RDM::d1
void d1(TensorOperator *denF, TensorT *denT, double *workmem)
Make diagram d1.
Definition: Tensor3RDM.cpp:280
CheMPS2::Tensor3RDM::~Tensor3RDM
virtual ~Tensor3RDM()
Destructor.
Definition: Tensor3RDM.cpp:44
CheMPS2::TensorOperator::sector_spin_up
int * sector_spin_up
The up spin symmetry sector.
Definition: TensorOperator.h:163
CheMPS2::TensorT::gStorage
double * gStorage()
Get the pointer to the storage.
Definition: TensorT.cpp:134
CheMPS2::TensorOperator::n_irrep
int n_irrep
The (real-valued abelian) point group irrep difference between the symmetry sectors of the lower and ...
Definition: TensorOperator.h:151
CheMPS2::Irreps::directProd
static int directProd(const int Irrep1, const int Irrep2)
Get the direct product of the irreps with numbers Irrep1 and Irrep2: a bitwise XOR for psi4&#39;s convent...
Definition: Irreps.h:123
CheMPS2::TensorOperator::sector_irrep_up
int * sector_irrep_up
The up spin symmetry sector.
Definition: TensorOperator.h:160
CheMPS2::Tensor3RDM::extra4
void extra4(TensorL *denL, TensorT *denT, double *workmem)
Make diagram extra4.
Definition: Tensor3RDM.cpp:469
CheMPS2::TensorOperator::bk_up
const SyBookkeeper * bk_up
The bookkeeper of the upper MPS.
Definition: TensorOperator.h:139
CheMPS2::Tensor::kappa2index
int * kappa2index
kappa2index[kappa] indicates the start of tensor block kappa in storage. kappa2index[nKappa] gives th...
Definition: Tensor.h:84
CheMPS2::Tensor3RDM::b1
void b1(TensorOperator *Sigma, TensorT *denT, double *workmem)
Make diagram b1.
Definition: Tensor3RDM.cpp:128
CheMPS2::Tensor3RDM::get_two_j2
int get_two_j2() const
Get the spin value of the Tensor3RDM.
Definition: Tensor3RDM.cpp:48
CheMPS2::TensorOperator::get_irrep
int get_irrep() const
Get the (real-valued abelian) point group irrep difference between the symmetry sectors of the lower ...
Definition: TensorOperator.cpp:155
CheMPS2::TensorOperator::clear
void clear()
Set all storage variables to 0.0.
Definition: TensorOperator.cpp:157
CheMPS2::Tensor::nKappa
int nKappa
Number of Tensor blocks.
Definition: Tensor.h:81
CheMPS2::SyBookkeeper::gIrrep
int gIrrep(const int orbital) const
Get an orbital irrep.
Definition: SyBookkeeper.cpp:145
CheMPS2::TensorOperator::n_elec
int n_elec
How many electrons there are more in the symmetry sector of the lower leg compared to the upper leg...
Definition: TensorOperator.h:148
CheMPS2::TensorOperator::two_j
int two_j
Twice the spin of the tensor operator.
Definition: TensorOperator.h:145
CheMPS2::Tensor::storage
double * storage
The actual variables. Tensor block kappa begins at storage+kappa2index[kappa] and ends at storage+kap...
Definition: Tensor.h:78