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 #include <assert.h>
23
24 #include "TensorL.h"
25 #include "Lapack.h"
26 #include "Special.h"
27
TensorL(const int boundary_index,const int Idiff,const bool moving_right,const SyBookkeeper * book_up,const SyBookkeeper * book_down)28 CheMPS2::TensorL::TensorL( const int boundary_index, const int Idiff, const bool moving_right, const SyBookkeeper * book_up, const SyBookkeeper * book_down ) :
29 TensorOperator( boundary_index,
30 1, //two_j
31 1, //n_elec
32 Idiff,
33 moving_right,
34 true, //prime_last
35 true, //jw_phase (one 2nd quantized operator)
36 book_up,
37 book_down ){ }
38
~TensorL()39 CheMPS2::TensorL::~TensorL(){ }
40
create(TensorT * mps_tensor)41 void CheMPS2::TensorL::create( TensorT * mps_tensor ){
42
43 clear();
44 assert( bk_up == bk_down );
45
46 if ( moving_right ){
47 for ( int ikappa = 0; ikappa < nKappa; ikappa++ ){ create_right( ikappa, mps_tensor, mps_tensor, NULL, NULL ); }
48 } else {
49 for ( int ikappa = 0; ikappa < nKappa; ikappa++ ){ create_left( ikappa, mps_tensor, mps_tensor, NULL, NULL ); }
50 }
51
52 }
53
create(TensorT * mps_tensor_up,TensorT * mps_tensor_down,TensorO * previous,double * workmem)54 void CheMPS2::TensorL::create( TensorT * mps_tensor_up, TensorT * mps_tensor_down, TensorO * previous, double * workmem ){
55
56 clear();
57
58 if ( moving_right ){
59 for ( int ikappa = 0; ikappa < nKappa; ikappa++ ){ create_right( ikappa, mps_tensor_up, mps_tensor_down, previous, workmem ); }
60 } else {
61 for ( int ikappa = 0; ikappa < nKappa; ikappa++ ){ create_left( ikappa, mps_tensor_up, mps_tensor_down, previous, workmem ); }
62 }
63
64 }
65
create_right(const int ikappa,TensorT * mps_tensor_up,TensorT * mps_tensor_down,TensorO * previous,double * workmem)66 void CheMPS2::TensorL::create_right( const int ikappa, TensorT * mps_tensor_up, TensorT * mps_tensor_down, TensorO * previous, double * workmem ){
67
68 const int NRup = sector_nelec_up[ ikappa ];
69 const int NRdown = NRup + 1;
70 const int IRup = sector_irrep_up[ ikappa ];
71 const int IRdown = Irreps::directProd( IRup, n_irrep );
72 const int TwoSRup = sector_spin_up[ ikappa ];
73 const int TwoSRdown = sector_spin_down[ ikappa ];
74
75 int dimRup = bk_up ->gCurrentDim( index, NRup, TwoSRup, IRup );
76 int dimRdown = bk_down->gCurrentDim( index, NRdown, TwoSRdown, IRdown );
77
78 for ( int geval = 0; geval < 2; geval++ ){
79 int NL, TwoSL, IL;
80 switch ( geval ){
81 case 0:
82 NL = NRup;
83 TwoSL = TwoSRup;
84 IL = IRup;
85 break;
86 case 1:
87 NL = NRup - 1;
88 TwoSL = TwoSRdown;
89 IL = IRdown;
90 break;
91 }
92
93 int dimLup = bk_up ->gCurrentDim( index - 1, NL, TwoSL, IL );
94 int dimLdown = bk_down->gCurrentDim( index - 1, NL, TwoSL, IL );
95
96 if ( previous == NULL ){
97 assert( dimLup == dimLdown );
98 if ( dimLup > 0 ){
99
100 double * Tup = mps_tensor_up ->gStorage( NL, TwoSL, IL, NRup, TwoSRup, IRup );
101 double * Tdown = mps_tensor_down->gStorage( NL, TwoSL, IL, NRdown, TwoSRdown, IRdown );
102
103 char trans = 'T';
104 char notrans = 'N';
105 double alpha = 1.0;
106 if ( geval == 1 ){
107 alpha = Special::phase( TwoSRdown - TwoSRup + 1 ) * sqrt( ( TwoSRup + 1.0 ) / ( TwoSRdown + 1 ) );
108 }
109 double add = 1.0;
110 dgemm_( &trans, ¬rans, &dimRup, &dimRdown, &dimLup, &alpha, Tup, &dimLup, Tdown, &dimLup, &add, storage + kappa2index[ ikappa ], &dimRup );
111
112 }
113 } else {
114 if (( dimLup > 0 ) && ( dimLdown > 0 )){
115
116 double * Tup = mps_tensor_up ->gStorage( NL, TwoSL, IL, NRup, TwoSRup, IRup );
117 double * Tdown = mps_tensor_down->gStorage( NL, TwoSL, IL, NRdown, TwoSRdown, IRdown );
118 double * Opart = previous->gStorage( NL, TwoSL, IL, NL, TwoSL, IL );
119
120 char trans = 'T';
121 char notrans = 'N';
122 double alpha = 1.0;
123 if ( geval == 1 ){
124 alpha = Special::phase( TwoSRdown - TwoSRup + 1 ) * sqrt( ( TwoSRup + 1.0 ) / ( TwoSRdown + 1 ) );
125 }
126 double set = 0.0;
127 dgemm_( &trans, ¬rans, &dimRup, &dimLdown, &dimLup, &alpha, Tup, &dimLup, Opart, &dimLup, &set, workmem, &dimRup );
128 double one = 1.0;
129 dgemm_( ¬rans, ¬rans, &dimRup, &dimRdown, &dimLdown, &one, workmem, &dimRup, Tdown, &dimLdown, &one, storage + kappa2index[ ikappa ], &dimRup );
130
131 }
132 }
133 }
134
135 }
136
create_left(const int ikappa,TensorT * mps_tensor_up,TensorT * mps_tensor_down,TensorO * previous,double * workmem)137 void CheMPS2::TensorL::create_left( const int ikappa, TensorT * mps_tensor_up, TensorT * mps_tensor_down, TensorO * previous, double * workmem ){
138
139 const int NLup = sector_nelec_up[ ikappa ];
140 const int NLdown = NLup + 1;
141 const int ILup = sector_irrep_up[ ikappa ];
142 const int ILdown = Irreps::directProd( ILup, n_irrep );
143 const int TwoSLup = sector_spin_up[ ikappa ];
144 const int TwoSLdown = sector_spin_down[ ikappa ];
145
146 int dimLup = bk_up ->gCurrentDim( index, NLup, TwoSLup, ILup );
147 int dimLdown = bk_down->gCurrentDim( index, NLdown, TwoSLdown, ILdown );
148
149 for ( int geval = 0; geval < 2; geval++ ){
150 int NR, TwoSR, IR;
151 switch ( geval ){
152 case 0:
153 NR = NLdown;
154 TwoSR = TwoSLdown;
155 IR = ILdown;
156 break;
157 case 1:
158 NR = NLup + 2;
159 TwoSR = TwoSLup;
160 IR = ILup;
161 break;
162 }
163
164 int dimRup = bk_up ->gCurrentDim( index + 1, NR, TwoSR, IR );
165 int dimRdown = bk_down->gCurrentDim( index + 1, NR, TwoSR, IR );
166
167 if ( previous == NULL ){
168 assert( dimRup == dimRdown );
169 if ( dimRup > 0 ){
170
171 double * Tup = mps_tensor_up ->gStorage( NLup, TwoSLup, ILup, NR, TwoSR, IR );
172 double * Tdown = mps_tensor_down->gStorage( NLdown, TwoSLdown, ILdown, NR, TwoSR, IR );
173
174 char trans = 'T';
175 char notrans = 'N';
176 double alpha = 1.0;
177 if ( geval == 1 ){
178 alpha = Special::phase( TwoSLup - TwoSLdown + 1 ) * sqrt( ( TwoSLup + 1.0 ) / ( TwoSLdown + 1 ) );
179 }
180 double add = 1.0;
181 dgemm_( ¬rans, &trans, &dimLup, &dimLdown, &dimRup, &alpha, Tup, &dimLup, Tdown, &dimLdown, &add, storage + kappa2index[ ikappa ], &dimLup );
182
183 }
184 } else {
185 if (( dimRup > 0 ) && ( dimRdown > 0 )){
186
187 double * Tup = mps_tensor_up ->gStorage( NLup, TwoSLup, ILup, NR, TwoSR, IR );
188 double * Tdown = mps_tensor_down->gStorage( NLdown, TwoSLdown, ILdown, NR, TwoSR, IR );
189 double * Opart = previous->gStorage( NR, TwoSR, IR, NR, TwoSR, IR );
190
191 char trans = 'T';
192 char notrans = 'N';
193 double alpha = 1.0;
194 if ( geval == 1 ){
195 alpha = Special::phase( TwoSLup - TwoSLdown + 1 ) * sqrt( ( TwoSLup + 1.0 ) / ( TwoSLdown + 1 ) );
196 }
197 double set = 0.0;
198 dgemm_( ¬rans, ¬rans, &dimLup, &dimRdown, &dimRup, &alpha, Tup, &dimLup, Opart, &dimRup, &set, workmem, &dimLup );
199 double one = 1.0;
200 dgemm_( ¬rans, &trans, &dimLup, &dimLdown, &dimRdown, &one, workmem, &dimLup, Tdown, &dimLdown, &one, storage + kappa2index[ ikappa ], &dimLup );
201
202 }
203 }
204 }
205
206 }
207
208
209