1 /*
2 CheMPS2: a spin-adapted implementation of DMRG for ab initio quantum chemistry
3 Copyright (C) 2013-2018 Sebastian Wouters
4
5 This program is free software; you can redistribute it and/or modify
6 it under the terms of the GNU General Public License as published by
7 the Free Software Foundation; either version 2 of the License, or
8 (at your option) any later version.
9
10 This program is distributed in the hope that it will be useful,
11 but WITHOUT ANY WARRANTY; without even the implied warranty of
12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 GNU General Public License for more details.
14
15 You should have received a copy of the GNU General Public License along
16 with this program; if not, write to the Free Software Foundation, Inc.,
17 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
18 */
19
20 #include <stdlib.h>
21 #include <math.h>
22
23 #include "TensorS0.h"
24 #include "Lapack.h"
25
TensorS0(const int boundary_index,const int Idiff,const bool moving_right,const SyBookkeeper * denBK)26 CheMPS2::TensorS0::TensorS0(const int boundary_index, const int Idiff, const bool moving_right, const SyBookkeeper * denBK) :
27 TensorOperator(boundary_index,
28 0, // two_j
29 2, // n_elec
30 Idiff,
31 moving_right,
32 true, // prime_last (doesn't matter for spin-0)
33 false, // jw_phase (two 2nd quantized operators)
34 denBK,
35 denBK){ }
36
~TensorS0()37 CheMPS2::TensorS0::~TensorS0(){ }
38
makenew(TensorT * denT)39 void CheMPS2::TensorS0::makenew(TensorT * denT){
40
41 if (moving_right){ makenewRight(denT); }
42 else{ makenewLeft( denT); }
43
44 }
45
makenew(TensorL * denL,TensorT * denT,double * workmem)46 void CheMPS2::TensorS0::makenew(TensorL * denL, TensorT * denT, double * workmem){
47
48 if (moving_right){ makenewRight(denL, denT, workmem); }
49 else{ makenewLeft( denL, denT, workmem); }
50
51 }
52
makenewRight(TensorT * denT)53 void CheMPS2::TensorS0::makenewRight(TensorT * denT){ //Idiff = Itrivial
54
55 clear();
56
57 for (int ikappa=0; ikappa<nKappa; ikappa++){
58 int dimUR = bk_up->gCurrentDim(index, sector_nelec_up[ikappa], sector_spin_up[ikappa], sector_irrep_up[ikappa]);
59 int dimDR = bk_up->gCurrentDim(index, sector_nelec_up[ikappa]+2, sector_spin_up[ikappa], sector_irrep_up[ikappa]);
60 int dimL = bk_up->gCurrentDim(index-1, sector_nelec_up[ikappa], sector_spin_up[ikappa], sector_irrep_up[ikappa]);
61 if (dimL>0){
62
63 double * BlockTup = denT->gStorage(sector_nelec_up[ikappa], sector_spin_up[ikappa], sector_irrep_up[ikappa], sector_nelec_up[ikappa], sector_spin_up[ikappa], sector_irrep_up[ikappa]);
64 double * BlockTdown = denT->gStorage(sector_nelec_up[ikappa], sector_spin_up[ikappa], sector_irrep_up[ikappa], sector_nelec_up[ikappa]+2, sector_spin_up[ikappa], sector_irrep_up[ikappa]);
65
66 char trans = 'T';
67 char notrans = 'N';
68 double alpha = sqrt(2.0);
69 double beta = 1.0; //add
70 dgemm_(&trans,¬rans,&dimUR,&dimDR,&dimL,&alpha,BlockTup,&dimL,BlockTdown,&dimL,&beta,storage+kappa2index[ikappa],&dimUR);
71
72 }
73 }
74
75 }
76
makenewLeft(TensorT * denT)77 void CheMPS2::TensorS0::makenewLeft(TensorT * denT){ //Idiff = Itrivial
78
79 clear();
80
81 for (int ikappa=0; ikappa<nKappa; ikappa++){
82 int dimUL = bk_up->gCurrentDim(index, sector_nelec_up[ikappa], sector_spin_up[ikappa], sector_irrep_up[ikappa]);
83 int dimDL = bk_up->gCurrentDim(index, sector_nelec_up[ikappa]+2, sector_spin_up[ikappa], sector_irrep_up[ikappa]);
84 int dimR = bk_up->gCurrentDim(index+1, sector_nelec_up[ikappa]+2, sector_spin_up[ikappa], sector_irrep_up[ikappa]);
85 if (dimR>0){
86
87 double * BlockTup = denT->gStorage(sector_nelec_up[ikappa], sector_spin_up[ikappa], sector_irrep_up[ikappa], sector_nelec_up[ikappa]+2, sector_spin_up[ikappa], sector_irrep_up[ikappa]);
88 double * BlockTdown = denT->gStorage(sector_nelec_up[ikappa]+2, sector_spin_up[ikappa], sector_irrep_up[ikappa], sector_nelec_up[ikappa]+2, sector_spin_up[ikappa], sector_irrep_up[ikappa]);
89
90 char trans = 'T';
91 char notrans = 'N';
92 double alpha = sqrt(2.0);
93 double beta = 1.0; //add
94 dgemm_(¬rans,&trans,&dimUL,&dimDL,&dimR,&alpha,BlockTup,&dimUL,BlockTdown,&dimDL,&beta,storage+kappa2index[ikappa],&dimUL);
95
96 }
97 }
98
99 }
100
makenewRight(TensorL * denL,TensorT * denT,double * workmem)101 void CheMPS2::TensorS0::makenewRight(TensorL * denL, TensorT * denT, double * workmem){
102
103 clear();
104
105 for (int ikappa=0; ikappa<nKappa; ikappa++){
106 const int IDR = Irreps::directProd(n_irrep,sector_irrep_up[ikappa]);
107 int dimUR = bk_up->gCurrentDim(index, sector_nelec_up[ikappa], sector_spin_up[ikappa], sector_irrep_up[ikappa]);
108 int dimDR = bk_up->gCurrentDim(index, sector_nelec_up[ikappa]+2, sector_spin_up[ikappa], IDR );
109
110 for (int geval=0; geval<4; geval++){
111 int NLU,TwoSLU,ILU,TwoSLD,ILD; //NLD = NLU+1
112 switch(geval){
113 case 0:
114 NLU = sector_nelec_up[ikappa];
115 TwoSLU = sector_spin_up[ikappa];
116 ILU = sector_irrep_up[ikappa];
117 TwoSLD = sector_spin_up[ikappa]-1;
118 ILD = Irreps::directProd( ILU, denL->get_irrep() );
119 break;
120 case 1:
121 NLU = sector_nelec_up[ikappa];
122 TwoSLU = sector_spin_up[ikappa];
123 ILU = sector_irrep_up[ikappa];
124 TwoSLD = sector_spin_up[ikappa]+1;
125 ILD = Irreps::directProd( ILU, denL->get_irrep() );
126 break;
127 case 2:
128 NLU = sector_nelec_up[ikappa]-1;
129 TwoSLU = sector_spin_up[ikappa]-1;
130 ILU = Irreps::directProd( sector_irrep_up[ikappa] , bk_up->gIrrep(index-1) );
131 TwoSLD = sector_spin_up[ikappa];
132 ILD = IDR;
133 break;
134 case 3:
135 NLU = sector_nelec_up[ikappa]-1;
136 TwoSLU = sector_spin_up[ikappa]+1;
137 ILU = Irreps::directProd( sector_irrep_up[ikappa] , bk_up->gIrrep(index-1) );
138 TwoSLD = sector_spin_up[ikappa];
139 ILD = IDR;
140 break;
141 }
142 int dimLU = bk_up->gCurrentDim(index-1, NLU, TwoSLU, ILU);
143 int dimLD = bk_up->gCurrentDim(index-1, NLU+1, TwoSLD, ILD);
144 if ((dimLU>0) && (dimLD>0)){
145
146 double * BlockTup = denT->gStorage(NLU, TwoSLU, ILU, sector_nelec_up[ikappa], sector_spin_up[ikappa], sector_irrep_up[ikappa]);
147 double * BlockTdown = denT->gStorage(NLU+1, TwoSLD, ILD, sector_nelec_up[ikappa]+2, sector_spin_up[ikappa], IDR);
148 double * BlockL = denL->gStorage(NLU, TwoSLU, ILU, NLU+1, TwoSLD, ILD);
149
150 //factor * Tup^T * L -> mem
151 char trans = 'T';
152 char notrans = 'N';
153 double alpha;
154 if (geval<=1){
155 int fase = ((((sector_spin_up[ikappa] - TwoSLD + 1)/2)%2)!=0)?-1:1;
156 alpha = fase * sqrt(0.5 * (TwoSLD+1.0) / (sector_spin_up[ikappa]+1.0) );
157 } else {
158 alpha = - sqrt(0.5);
159 }
160 double beta = 0.0; //set
161 dgemm_(&trans,¬rans,&dimUR,&dimLD,&dimLU,&alpha,BlockTup,&dimLU,BlockL,&dimLU,&beta,workmem,&dimUR);
162
163 //mem * Tdown -> storage
164 alpha = 1.0;
165 beta = 1.0; // add
166 dgemm_(¬rans,¬rans,&dimUR,&dimDR,&dimLD,&alpha,workmem,&dimUR,BlockTdown,&dimLD,&beta,storage+kappa2index[ikappa],&dimUR);
167
168 }
169 }
170 }
171
172 }
173
makenewLeft(TensorL * denL,TensorT * denT,double * workmem)174 void CheMPS2::TensorS0::makenewLeft(TensorL * denL, TensorT * denT, double * workmem){
175
176 clear();
177
178 for (int ikappa=0; ikappa<nKappa; ikappa++){
179 const int IDL = Irreps::directProd(n_irrep,sector_irrep_up[ikappa]);
180 int dimUL = bk_up->gCurrentDim(index, sector_nelec_up[ikappa], sector_spin_up[ikappa], sector_irrep_up[ikappa]);
181 int dimDL = bk_up->gCurrentDim(index, sector_nelec_up[ikappa]+2, sector_spin_up[ikappa], IDL );
182
183 for (int geval=0; geval<4; geval++){
184 int NRU,TwoSRU,IRU,TwoSRD,IRD; //NRD = NRU+1
185 switch(geval){
186 case 0:
187 NRU = sector_nelec_up[ikappa]+1;
188 TwoSRU = sector_spin_up[ikappa]-1;
189 IRU = Irreps::directProd( sector_irrep_up[ikappa] , bk_up->gIrrep(index) );
190 TwoSRD = sector_spin_up[ikappa];
191 IRD = IDL;
192 break;
193 case 1:
194 NRU = sector_nelec_up[ikappa]+1;
195 TwoSRU = sector_spin_up[ikappa]+1;
196 IRU = Irreps::directProd( sector_irrep_up[ikappa] , bk_up->gIrrep(index) );
197 TwoSRD = sector_spin_up[ikappa];
198 IRD = IDL;
199 break;
200 case 2:
201 NRU = sector_nelec_up[ikappa]+2;
202 TwoSRU = sector_spin_up[ikappa];
203 IRU = sector_irrep_up[ikappa];
204 TwoSRD = sector_spin_up[ikappa]-1;
205 IRD = Irreps::directProd( sector_irrep_up[ikappa] , denL->get_irrep() );
206 break;
207 case 3:
208 NRU = sector_nelec_up[ikappa]+2;
209 TwoSRU = sector_spin_up[ikappa];
210 IRU = sector_irrep_up[ikappa];
211 TwoSRD = sector_spin_up[ikappa]+1;
212 IRD = Irreps::directProd( sector_irrep_up[ikappa] , denL->get_irrep() );
213 break;
214 }
215 int dimRU = bk_up->gCurrentDim(index+1, NRU, TwoSRU, IRU);
216 int dimRD = bk_up->gCurrentDim(index+1, NRU+1, TwoSRD, IRD);
217 if ((dimRU>0) && (dimRD>0)){
218
219 double * BlockTup = denT->gStorage(sector_nelec_up[ikappa], sector_spin_up[ikappa], sector_irrep_up[ikappa], NRU, TwoSRU, IRU);
220 double * BlockTdown = denT->gStorage(sector_nelec_up[ikappa]+2, sector_spin_up[ikappa], IDL, NRU+1, TwoSRD, IRD);
221 double * BlockL = denL->gStorage(NRU, TwoSRU, IRU, NRU+1, TwoSRD, IRD);
222
223 //factor * Tup * L -> mem
224 char notrans = 'N';
225 double alpha = 1.0;
226 if (geval<=1){
227 int fase = ((((sector_spin_up[ikappa] - TwoSRU + 1)/2)%2)!=0)?-1:1;
228 alpha = fase * sqrt(0.5 * (TwoSRU+1.0) / (sector_spin_up[ikappa]+1.0) );
229 } else {
230 alpha = - sqrt(0.5) * (TwoSRD+1.0) / (sector_spin_up[ikappa]+1.0);
231 }
232 double beta = 0.0; //set
233 dgemm_(¬rans,¬rans,&dimUL,&dimRD,&dimRU,&alpha,BlockTup,&dimUL,BlockL,&dimRU,&beta,workmem,&dimUL);
234
235 //mem * Tdown^T -> storage
236 char trans = 'T';
237 alpha = 1.0;
238 beta = 1.0; // add
239 dgemm_(¬rans,&trans,&dimUL,&dimDL,&dimRD,&alpha,workmem,&dimUL,BlockTdown,&dimDL,&beta,storage+kappa2index[ikappa],&dimUL);
240
241 }
242 }
243 }
244
245 }
246
247