CheMPS2
TensorF0.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 <math.h>
22 
23 #include "TensorF0.h"
24 #include "Lapack.h"
25 
26 CheMPS2::TensorF0::TensorF0( const int boundary_index, const int Idiff, const bool moving_right, const SyBookkeeper * denBK ) :
27 TensorOperator(boundary_index,
28  0, // two_j
29  0, // 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 
37 CheMPS2::TensorF0::~TensorF0(){ }
38 
39 void CheMPS2::TensorF0::makenew(TensorT * denT){
40 
41  if (moving_right){ makenewRight(denT); }
42  else{ makenewLeft( denT); }
43 
44 }
45 
46 void CheMPS2::TensorF0::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 
53 void CheMPS2::TensorF0::makenewRight(TensorT * denT){ //Idiff = Itrivial
54 
55  clear();
56 
57  for (int ikappa=0; ikappa<nKappa; ikappa++){
58  int dimR = bk_up->gCurrentDim(index, sector_nelec_up[ikappa], sector_spin_up[ikappa], sector_irrep_up[ikappa]);
59  const double sqrt_of_2 = sqrt(2.0);
60  for (int geval=0; geval<3; geval++){
61  int TwoSL, NL, IL;
62  switch(geval){
63  case 0:
64  NL = sector_nelec_up[ikappa]-2;
65  TwoSL = sector_spin_up[ikappa];
66  IL = sector_irrep_up[ikappa];
67  break;
68  case 1:
69  NL = sector_nelec_up[ikappa]-1;
70  TwoSL = sector_spin_up[ikappa]-1;
71  IL = Irreps::directProd( sector_irrep_up[ikappa] , bk_up->gIrrep(index-1) );
72  break;
73  case 2:
74  NL = sector_nelec_up[ikappa]-1;
75  TwoSL = sector_spin_up[ikappa]+1;
76  IL = Irreps::directProd( sector_irrep_up[ikappa] , bk_up->gIrrep(index-1) );
77  break;
78  }
79  int dimL = bk_up->gCurrentDim(index-1, NL, TwoSL, IL);
80  if (dimL>0){
81 
82  double * BlockT = denT->gStorage(NL, TwoSL, IL, sector_nelec_up[ikappa], sector_spin_up[ikappa], sector_irrep_up[ikappa]);
83 
84  char trans = 'T';
85  char notrans = 'N';
86  double alpha = (geval==0)?sqrt_of_2:(0.5*sqrt_of_2);
87  double beta = 1.0; //add
88  dgemm_(&trans,&notrans,&dimR,&dimR,&dimL,&alpha,BlockT,&dimL,BlockT,&dimL,&beta,storage+kappa2index[ikappa],&dimR);
89 
90  }
91  }
92  }
93 
94 }
95 
96 void CheMPS2::TensorF0::makenewLeft(TensorT * denT){ //Idiff = Itrivial
97 
98  clear();
99 
100  for (int ikappa=0; ikappa<nKappa; ikappa++){
101  int dimL = bk_up->gCurrentDim(index, sector_nelec_up[ikappa], sector_spin_up[ikappa], sector_irrep_up[ikappa]);
102  const double sqrt_of_2 = sqrt(2.0);
103  for (int geval=0; geval<3; geval++){
104  int TwoSR, NR, IR;
105  switch(geval){
106  case 0:
107  NR = sector_nelec_up[ikappa]+2;
108  TwoSR = sector_spin_up[ikappa];
109  IR = sector_irrep_up[ikappa];
110  break;
111  case 1:
112  NR = sector_nelec_up[ikappa]+1;
113  TwoSR = sector_spin_up[ikappa]-1;
114  IR = Irreps::directProd( sector_irrep_up[ikappa] , bk_up->gIrrep(index) );
115  break;
116  case 2:
117  NR = sector_nelec_up[ikappa]+1;
118  TwoSR = sector_spin_up[ikappa]+1;
119  IR = Irreps::directProd( sector_irrep_up[ikappa] , bk_up->gIrrep(index) );
120  break;
121  }
122  int dimR = bk_up->gCurrentDim(index+1, NR, TwoSR, IR);
123  if (dimR>0){
124 
125  double * BlockT = denT->gStorage(sector_nelec_up[ikappa], sector_spin_up[ikappa], sector_irrep_up[ikappa], NR, TwoSR, IR);
126 
127  char trans = 'T';
128  char notrans = 'N';
129  double alpha = sqrt_of_2;
130  if (geval>=1){ alpha *= 0.5 * (TwoSR + 1.0) / (sector_spin_up[ikappa] + 1.0); }
131  double beta = 1.0; //add
132  dgemm_(&notrans,&trans,&dimL,&dimL,&dimR,&alpha,BlockT,&dimL,BlockT,&dimL,&beta,storage+kappa2index[ikappa],&dimL);
133 
134  }
135  }
136  }
137 
138 }
139 
140 void CheMPS2::TensorF0::makenewRight(TensorL * denL, TensorT * denT, double * workmem){
141 
142  clear();
143 
144  for (int ikappa=0; ikappa<nKappa; ikappa++){
145  const int IDR = Irreps::directProd(n_irrep,sector_irrep_up[ikappa]);
146  int dimUR = bk_up->gCurrentDim(index, sector_nelec_up[ikappa], sector_spin_up[ikappa], sector_irrep_up[ikappa]);
147  int dimDR = bk_up->gCurrentDim(index, sector_nelec_up[ikappa], sector_spin_up[ikappa], IDR );
148 
149  for (int geval=0; geval<4; geval++){
150  int NLU,TwoSLU,ILU,TwoSLD,ILD; //NLD = NLU+1
151  switch(geval){
152  case 0:
153  NLU = sector_nelec_up[ikappa]-1;
154  TwoSLU = sector_spin_up[ikappa]-1;
155  ILU = Irreps::directProd( sector_irrep_up[ikappa], bk_up->gIrrep(index-1) );
156  TwoSLD = sector_spin_up[ikappa];
157  ILD = IDR;
158  break;
159  case 1:
160  NLU = sector_nelec_up[ikappa]-1;
161  TwoSLU = sector_spin_up[ikappa]+1;
162  ILU = Irreps::directProd( sector_irrep_up[ikappa], bk_up->gIrrep(index-1) );
163  TwoSLD = sector_spin_up[ikappa];
164  ILD = IDR;
165  break;
166  case 2:
167  NLU = sector_nelec_up[ikappa]-2;
168  TwoSLU = sector_spin_up[ikappa];
169  ILU = sector_irrep_up[ikappa];
170  TwoSLD = sector_spin_up[ikappa]-1;
171  ILD = Irreps::directProd( ILU, denL->get_irrep() );
172  break;
173  case 3:
174  NLU = sector_nelec_up[ikappa]-2;
175  TwoSLU = sector_spin_up[ikappa];
176  ILU = sector_irrep_up[ikappa];
177  TwoSLD = sector_spin_up[ikappa]+1;
178  ILD = Irreps::directProd( ILU, denL->get_irrep() );
179  break;
180  }
181  int dimLU = bk_up->gCurrentDim(index-1, NLU, TwoSLU, ILU);
182  int dimLD = bk_up->gCurrentDim(index-1, NLU+1, TwoSLD, ILD);
183  if ((dimLU>0) && (dimLD>0)){
184 
185  double * BlockTup = denT->gStorage(NLU, TwoSLU, ILU, sector_nelec_up[ikappa], sector_spin_up[ikappa], sector_irrep_up[ikappa]);
186  double * BlockTdown = denT->gStorage(NLU+1, TwoSLD, ILD, sector_nelec_up[ikappa], sector_spin_up[ikappa], IDR);
187  double * BlockL = denL->gStorage(NLU, TwoSLU, ILU, NLU+1, TwoSLD, ILD);
188 
189  //factor * Tup^T * L -> mem
190  char trans = 'T';
191  char notrans = 'N';
192  double alpha = 1.0;
193  if (geval<=1){
194  alpha = sqrt(0.5);
195  } else {
196  int fase = ((((sector_spin_up[ikappa] + 1 - TwoSLD)/2)%2)!=0)?-1:1;
197  alpha = fase * sqrt( 0.5 * ( TwoSLD + 1.0 ) / ( sector_spin_up[ikappa] + 1.0 ) );
198  }
199  double beta = 0.0; //set
200  dgemm_(&trans,&notrans,&dimUR,&dimLD,&dimLU,&alpha,BlockTup,&dimLU,BlockL,&dimLU,&beta,workmem,&dimUR);
201 
202  //mem * Tdown -> storage
203  alpha = 1.0;
204  beta = 1.0; // add
205  dgemm_(&notrans,&notrans,&dimUR,&dimDR,&dimLD,&alpha,workmem,&dimUR,BlockTdown,&dimLD,&beta,storage+kappa2index[ikappa],&dimUR);
206 
207  }
208  }
209  }
210 
211 }
212 
213 void CheMPS2::TensorF0::makenewLeft(TensorL * denL, TensorT * denT, double * workmem){
214 
215  clear();
216 
217  for (int ikappa=0; ikappa<nKappa; ikappa++){
218  const int IDL = Irreps::directProd(n_irrep,sector_irrep_up[ikappa]);
219  int dimUL = bk_up->gCurrentDim(index, sector_nelec_up[ikappa], sector_spin_up[ikappa], sector_irrep_up[ikappa]);
220  int dimDL = bk_up->gCurrentDim(index, sector_nelec_up[ikappa], sector_spin_up[ikappa], IDL );
221 
222  for (int geval=0; geval<4; geval++){
223  int NRU,TwoSRU,IRU,TwoSRD,IRD; //NRD = NRU+1
224  switch(geval){
225  case 0:
226  NRU = sector_nelec_up[ikappa];
227  TwoSRU = sector_spin_up[ikappa];
228  IRU = sector_irrep_up[ikappa];
229  TwoSRD = sector_spin_up[ikappa] - 1;
230  IRD = Irreps::directProd( IRU , denL->get_irrep() );
231  break;
232  case 1:
233  NRU = sector_nelec_up[ikappa];
234  TwoSRU = sector_spin_up[ikappa];
235  IRU = sector_irrep_up[ikappa];
236  TwoSRD = sector_spin_up[ikappa] + 1;
237  IRD = Irreps::directProd( IRU , denL->get_irrep() );
238  break;
239  case 2:
240  NRU = sector_nelec_up[ikappa] + 1;
241  TwoSRU = sector_spin_up[ikappa] - 1;
242  IRU = Irreps::directProd( sector_irrep_up[ikappa] , bk_up->gIrrep(index) );
243  TwoSRD = sector_spin_up[ikappa];
244  IRD = IDL;
245  break;
246  case 3:
247  NRU = sector_nelec_up[ikappa] + 1;
248  TwoSRU = sector_spin_up[ikappa] + 1;
249  IRU = Irreps::directProd( sector_irrep_up[ikappa] , bk_up->gIrrep(index) );
250  TwoSRD = sector_spin_up[ikappa];
251  IRD = IDL;
252  break;
253  }
254  int dimRU = bk_up->gCurrentDim(index+1, NRU, TwoSRU, IRU);
255  int dimRD = bk_up->gCurrentDim(index+1, NRU+1, TwoSRD, IRD);
256  if ((dimRU>0) && (dimRD>0)){
257 
258  double * BlockTup = denT->gStorage(sector_nelec_up[ikappa], sector_spin_up[ikappa], sector_irrep_up[ikappa], NRU, TwoSRU, IRU);
259  double * BlockTdown = denT->gStorage(sector_nelec_up[ikappa], sector_spin_up[ikappa], IDL, NRU+1, TwoSRD, IRD);
260  double * BlockL = denL->gStorage(NRU, TwoSRU, IRU, NRU+1, TwoSRD, IRD);
261 
262  //factor * Tup * L -> mem
263  char notrans = 'N';
264  double alpha = 1.0;
265  if (geval<=1){
266  alpha = sqrt(0.5) * (TwoSRD + 1.0) / (sector_spin_up[ikappa] + 1.0);
267  } else {
268  int fase = ((((sector_spin_up[ikappa] - TwoSRU + 1)/2)%2)!=0)?-1:1;
269  alpha = fase * sqrt( 0.5 * (TwoSRU + 1.0) / (sector_spin_up[ikappa] + 1.0) ) ;
270  }
271  double beta = 0.0; //set
272  dgemm_(&notrans,&notrans,&dimUL,&dimRD,&dimRU,&alpha,BlockTup,&dimUL,BlockL,&dimRU,&beta,workmem,&dimUL);
273 
274  //mem * Tdown^T -> storage
275  char trans = 'T';
276  alpha = 1.0;
277  beta = 1.0; // add
278  dgemm_(&notrans,&trans,&dimUL,&dimDL,&dimRD,&alpha,workmem,&dimUL,BlockTdown,&dimDL,&beta,storage+kappa2index[ikappa],&dimUL);
279 
280  }
281  }
282  }
283 
284 }
285 
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::TensorF0::~TensorF0
virtual ~TensorF0()
Destructor.
Definition: TensorF0.cpp:37
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::gStorage
double * gStorage()
Get the pointer to the storage.
Definition: TensorOperator.cpp:117
CheMPS2::TensorOperator::moving_right
bool moving_right
Whether or not moving right.
Definition: TensorOperator.h:154
CheMPS2::TensorF0::TensorF0
TensorF0(const int boundary_index, const int Idiff, const bool moving_right, const SyBookkeeper *denBK)
Constructor.
Definition: TensorF0.cpp:26
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::TensorOperator::bk_up
const SyBookkeeper * bk_up
The bookkeeper of the upper MPS.
Definition: TensorOperator.h:139
CheMPS2::TensorF0::makenew
void makenew(TensorT *denT)
Definition: TensorF0.cpp:39
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::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::Tensor::storage
double * storage
The actual variables. Tensor block kappa begins at storage+kappa2index[kappa] and ends at storage+kap...
Definition: Tensor.h:78