CheMPS2
TensorS0.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 "TensorS0.h"
24 #include "Lapack.h"
25 
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 
37 CheMPS2::TensorS0::~TensorS0(){ }
38 
39 void CheMPS2::TensorS0::makenew(TensorT * denT){
40 
41  if (moving_right){ makenewRight(denT); }
42  else{ makenewLeft( denT); }
43 
44 }
45 
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 
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,&notrans,&dimUR,&dimDR,&dimL,&alpha,BlockTup,&dimL,BlockTdown,&dimL,&beta,storage+kappa2index[ikappa],&dimUR);
71 
72  }
73  }
74 
75 }
76 
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_(&notrans,&trans,&dimUL,&dimDL,&dimR,&alpha,BlockTup,&dimUL,BlockTdown,&dimDL,&beta,storage+kappa2index[ikappa],&dimUL);
95 
96  }
97  }
98 
99 }
100 
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,&notrans,&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_(&notrans,&notrans,&dimUR,&dimDR,&dimLD,&alpha,workmem,&dimUR,BlockTdown,&dimLD,&beta,storage+kappa2index[ikappa],&dimUR);
167 
168  }
169  }
170  }
171 
172 }
173 
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_(&notrans,&notrans,&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_(&notrans,&trans,&dimUL,&dimDL,&dimRD,&alpha,workmem,&dimUL,BlockTdown,&dimDL,&beta,storage+kappa2index[ikappa],&dimUL);
240 
241  }
242  }
243  }
244 
245 }
246 
CheMPS2::TensorS0::makenew
void makenew(TensorT *denT)
Definition: TensorS0.cpp:39
CheMPS2::TensorS0::TensorS0
TensorS0(const int boundary_index, const int Idiff, const bool moving_right, const SyBookkeeper *denBK)
Constructor.
Definition: TensorS0.cpp:26
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::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::TensorS0::~TensorS0
virtual ~TensorS0()
Destructor.
Definition: TensorS0.cpp:37
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::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::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