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