CheMPS2
HeffDiagonal.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 <math.h>
21 
22 #include "Heff.h"
23 #include "Lapack.h"
24 #include "MPIchemps2.h"
25 #include "Wigner.h"
26 
27 void CheMPS2::Heff::addDiagonal1A(const int ikappa, double * memHeffDiag, const Sobject * denS, TensorX * Xleft) const{
28  int dimL = denBK->gCurrentDim(denS->gIndex(), denS->gNL(ikappa), denS->gTwoSL(ikappa), denS->gIL(ikappa));
29  int dimR = denBK->gCurrentDim(denS->gIndex()+2, denS->gNR(ikappa), denS->gTwoSR(ikappa), denS->gIR(ikappa));
30  double * BlockX = Xleft->gStorage( denS->gNL(ikappa), denS->gTwoSL(ikappa), denS->gIL(ikappa), denS->gNL(ikappa), denS->gTwoSL(ikappa), denS->gIL(ikappa) );
31  int ptr = denS->gKappa2index(ikappa);
32 
33  for (int cnt=0; cnt<dimL; cnt++){
34  for (int cnt2=0; cnt2<dimR; cnt2++){
35  memHeffDiag[ptr + cnt + dimL*cnt2] += BlockX[cnt*(dimL+1)];
36  }
37  }
38 
39 }
40 
41 void CheMPS2::Heff::addDiagonal1B(const int ikappa, double * memHeffDiag, const Sobject * denS, TensorX * Xright) const{
42  int dimL = denBK->gCurrentDim(denS->gIndex(), denS->gNL(ikappa), denS->gTwoSL(ikappa), denS->gIL(ikappa));
43  int dimR = denBK->gCurrentDim(denS->gIndex()+2, denS->gNR(ikappa), denS->gTwoSR(ikappa), denS->gIR(ikappa));
44  double * BlockX = Xright->gStorage( denS->gNR(ikappa), denS->gTwoSR(ikappa), denS->gIR(ikappa), denS->gNR(ikappa), denS->gTwoSR(ikappa), denS->gIR(ikappa) );
45  int ptr = denS->gKappa2index(ikappa);
46 
47  for (int cnt=0; cnt<dimL; cnt++){
48  for (int cnt2=0; cnt2<dimR; cnt2++){
49  memHeffDiag[ptr + cnt + dimL*cnt2] += BlockX[cnt2*(dimR+1)];
50  }
51  }
52 
53 }
54 
55 void CheMPS2::Heff::addDiagonal1C(const int ikappa, double * memHeffDiag, const Sobject * denS, const double Helem_links) const{
56  if (denS->gN1(ikappa)==2){
57  int ptr = denS->gKappa2index(ikappa);
58  int dim = denS->gKappa2index(ikappa+1) - ptr;
59  for (int cnt=0; cnt<dim; cnt++){ memHeffDiag[ptr + cnt] += Helem_links; }
60  }
61 }
62 
63 void CheMPS2::Heff::addDiagonal1D(const int ikappa, double * memHeffDiag, const Sobject * denS, const double Helem_rechts) const{
64  if (denS->gN2(ikappa)==2){
65  int ptr = denS->gKappa2index(ikappa);
66  int dim = denS->gKappa2index(ikappa+1) - ptr;
67  for (int cnt=0; cnt<dim; cnt++){ memHeffDiag[ptr + cnt] += Helem_rechts; }
68  }
69 }
70 
71 void CheMPS2::Heff::addDiagonal2d3all(const int ikappa, double * memHeffDiag, const Sobject * denS) const{
72 
73  if ((denS->gN1(ikappa)==2)&&(denS->gN2(ikappa)==2)){ //2d3a
74 
75  const int theindex = denS->gIndex();
76  const int ptr = denS->gKappa2index(ikappa);
77  const int dim = denS->gKappa2index(ikappa+1) - ptr;
78  const double factor = 4 * Prob->gMxElement(theindex,theindex+1,theindex,theindex+1)
79  - 2 * Prob->gMxElement(theindex,theindex+1,theindex+1,theindex);
80 
81  for (int cnt=0; cnt<dim; cnt++){ memHeffDiag[ptr + cnt] += factor; }
82 
83  }
84 
85  if ((denS->gN1(ikappa)==1)&&(denS->gN2(ikappa)==1)){ //2d3b
86 
87  const int theindex = denS->gIndex();
88  const int ptr = denS->gKappa2index(ikappa);
89  const int dim = denS->gKappa2index(ikappa+1) - ptr;
90  const int fase = (denS->gTwoJ(ikappa) == 2)? -1: 1;
91  const double factor = Prob->gMxElement(theindex,theindex+1,theindex,theindex+1)
92  + fase * Prob->gMxElement(theindex,theindex+1,theindex+1,theindex);
93 
94  for (int cnt=0; cnt<dim; cnt++){ memHeffDiag[ptr + cnt] += factor; }
95 
96  }
97 
98  if ((denS->gN1(ikappa)==2)&&(denS->gN2(ikappa)==1)){ //2d3c
99 
100  const int theindex = denS->gIndex();
101  const int ptr = denS->gKappa2index(ikappa);
102  const int dim = denS->gKappa2index(ikappa+1) - ptr;
103  const double factor = 2 * Prob->gMxElement(theindex,theindex+1,theindex,theindex+1)
104  - Prob->gMxElement(theindex,theindex+1,theindex+1,theindex);
105 
106  for (int cnt=0; cnt<dim; cnt++){ memHeffDiag[ptr + cnt] += factor; }
107 
108  }
109 
110  if ((denS->gN1(ikappa)==1)&&(denS->gN2(ikappa)==2)){ //2d3d
111 
112  const int theindex = denS->gIndex();
113  const int ptr = denS->gKappa2index(ikappa);
114  const int dim = denS->gKappa2index(ikappa+1) - ptr;
115  const double factor = 2 * Prob->gMxElement(theindex,theindex+1,theindex,theindex+1)
116  - Prob->gMxElement(theindex,theindex+1,theindex+1,theindex);
117 
118  for (int cnt=0; cnt<dim; cnt++){ memHeffDiag[ptr + cnt] += factor; }
119 
120  }
121 
122 }
123 
124 void CheMPS2::Heff::addDiagonal2b3spin0(const int ikappa, double * memHeffDiag, const Sobject * denS, TensorOperator * Ctensor) const{
125 
126  int N1 = denS->gN1(ikappa);
127 
128  if (N1!=0){
129 
130  int theindex = denS->gIndex();
131  int ptr = denS->gKappa2index(ikappa);
132 
133  double sqrt0p5 = sqrt(0.5);
134 
135  int NL = denS->gNL(ikappa);
136  int TwoSL = denS->gTwoSL(ikappa);
137  int IL = denS->gIL(ikappa);
138 
139  int dimL = denBK->gCurrentDim(theindex, NL, TwoSL, IL);
140  int dimR = denBK->gCurrentDim(theindex+2,denS->gNR(ikappa),denS->gTwoSR(ikappa),denS->gIR(ikappa));
141 
142  double * Cblock = Ctensor->gStorage(NL,TwoSL,IL,NL,TwoSL,IL);
143  for (int cntR=0; cntR<dimR; cntR++){
144  for (int cntL=0; cntL<dimL; cntL++){
145  memHeffDiag[ptr + cntL + dimL*cntR] += N1*sqrt0p5*Cblock[(dimL+1)*cntL];
146  }
147  }
148 
149  }
150 
151 }
152 
153 void CheMPS2::Heff::addDiagonal2c3spin0(const int ikappa, double * memHeffDiag, const Sobject * denS, TensorOperator * Ctensor) const{
154 
155  int N2 = denS->gN2(ikappa);
156 
157  if (N2!=0){
158 
159  int theindex = denS->gIndex();
160  int ptr = denS->gKappa2index(ikappa);
161 
162  double sqrt0p5 = sqrt(0.5);
163 
164  int NL = denS->gNL(ikappa);
165  int TwoSL = denS->gTwoSL(ikappa);
166  int IL = denS->gIL(ikappa);
167 
168  int dimL = denBK->gCurrentDim(theindex, NL, TwoSL, IL);
169  int dimR = denBK->gCurrentDim(theindex+2,denS->gNR(ikappa),denS->gTwoSR(ikappa),denS->gIR(ikappa));
170 
171  double * Cblock = Ctensor->gStorage(NL,TwoSL,IL,NL,TwoSL,IL);
172  for (int cntR=0; cntR<dimR; cntR++){
173  for (int cntL=0; cntL<dimL; cntL++){
174  memHeffDiag[ptr + cntL + dimL*cntR] += N2*sqrt0p5*Cblock[(dimL+1)*cntL];
175  }
176  }
177 
178  }
179 
180 }
181 
182 void CheMPS2::Heff::addDiagonal2e3spin0(const int ikappa, double * memHeffDiag, const Sobject * denS, TensorOperator * Ctensor) const{
183 
184  int N1 = denS->gN1(ikappa);
185 
186  if (N1!=0){
187 
188  int theindex = denS->gIndex();
189  int ptr = denS->gKappa2index(ikappa);
190 
191  double sqrt0p5 = sqrt(0.5);
192 
193  int NR = denS->gNR(ikappa);
194  int TwoSR = denS->gTwoSR(ikappa);
195  int IR = denS->gIR(ikappa);
196 
197  int dimR = denBK->gCurrentDim(theindex+2,NR, TwoSR, IR);
198  int dimL = denBK->gCurrentDim(theindex ,denS->gNL(ikappa),denS->gTwoSL(ikappa),denS->gIL(ikappa));
199 
200  double * Cblock = Ctensor->gStorage(NR,TwoSR,IR,NR,TwoSR,IR);
201  for (int cntR=0; cntR<dimR; cntR++){
202  for (int cntL=0; cntL<dimL; cntL++){
203  memHeffDiag[ptr + cntL + dimL*cntR] += N1*sqrt0p5*Cblock[(dimR+1)*cntR];
204  }
205  }
206 
207  }
208 
209 }
210 
211 void CheMPS2::Heff::addDiagonal2f3spin0(const int ikappa, double * memHeffDiag, const Sobject * denS, TensorOperator * Ctensor) const{
212 
213  int N2 = denS->gN2(ikappa);
214 
215  if (N2!=0){
216 
217  int theindex = denS->gIndex();
218  int ptr = denS->gKappa2index(ikappa);
219 
220  double sqrt0p5 = sqrt(0.5);
221 
222  int NR = denS->gNR(ikappa);
223  int TwoSR = denS->gTwoSR(ikappa);
224  int IR = denS->gIR(ikappa);
225 
226  int dimR = denBK->gCurrentDim(theindex+2,NR, TwoSR, IR);
227  int dimL = denBK->gCurrentDim(theindex ,denS->gNL(ikappa),denS->gTwoSL(ikappa),denS->gIL(ikappa));
228 
229  double * Cblock = Ctensor->gStorage(NR,TwoSR,IR,NR,TwoSR,IR);
230  for (int cntR=0; cntR<dimR; cntR++){
231  for (int cntL=0; cntL<dimL; cntL++){
232  memHeffDiag[ptr + cntL + dimL*cntR] += N2*sqrt0p5*Cblock[(dimR+1)*cntR];
233  }
234  }
235 
236  }
237 
238 }
239 
240 void CheMPS2::Heff::addDiagonal2b3spin1(const int ikappa, double * memHeffDiag, const Sobject * denS, TensorOperator * Dtensor) const{
241 
242  int N1 = denS->gN1(ikappa);
243 
244  if (N1==1){
245 
246  int theindex = denS->gIndex();
247  int ptr = denS->gKappa2index(ikappa);
248 
249  int NL = denS->gNL(ikappa);
250  int TwoSL = denS->gTwoSL(ikappa);
251  int IL = denS->gIL(ikappa);
252  int NR = denS->gNR(ikappa);
253  int TwoSR = denS->gTwoSR(ikappa);
254  int IR = denS->gIR(ikappa);
255  int TwoJ = denS->gTwoJ(ikappa);
256  int N2 = denS->gN2(ikappa);
257 
258  int dimL = denBK->gCurrentDim(theindex, NL,TwoSL,IL);
259  int dimR = denBK->gCurrentDim(theindex+2,NR,TwoSR,IR);
260 
261  int fase = phase(TwoSL + TwoSR + 2*TwoJ + ((N2==1)?1:0) - 1);
262  const double alpha = fase * (TwoJ+1) * sqrt(3.0*(TwoSL+1)) * Wigner::wigner6j(TwoJ,TwoJ,2,1,1,((N2==1)?1:0)) * Wigner::wigner6j(TwoJ,TwoJ,2,TwoSL,TwoSL,TwoSR);
263 
264  double * Dblock = Dtensor->gStorage(NL,TwoSL,IL,NL,TwoSL,IL);
265  for (int cntR=0; cntR<dimR; cntR++){
266  for (int cntL=0; cntL<dimL; cntL++){
267  memHeffDiag[ptr + cntL + dimL*cntR] += alpha * Dblock[(dimL+1)*cntL];
268  }
269  }
270 
271  }
272 
273 }
274 
275 void CheMPS2::Heff::addDiagonal2c3spin1(const int ikappa, double * memHeffDiag, const Sobject * denS, TensorOperator * Dtensor) const{
276 
277  int N2 = denS->gN2(ikappa);
278 
279  if (N2==1){
280 
281  int theindex = denS->gIndex();
282  int ptr = denS->gKappa2index(ikappa);
283 
284  int NL = denS->gNL(ikappa);
285  int TwoSL = denS->gTwoSL(ikappa);
286  int IL = denS->gIL(ikappa);
287  int NR = denS->gNR(ikappa);
288  int TwoSR = denS->gTwoSR(ikappa);
289  int IR = denS->gIR(ikappa);
290  int TwoJ = denS->gTwoJ(ikappa);
291  int N1 = denS->gN1(ikappa);
292 
293  int dimL = denBK->gCurrentDim(theindex, NL,TwoSL,IL);
294  int dimR = denBK->gCurrentDim(theindex+2,NR,TwoSR,IR);
295 
296  int fase = phase(TwoSL + TwoSR + 2*TwoJ + ((N1==1)?1:0) - 1);
297  const double alpha = fase * (TwoJ+1) * sqrt(3.0*(TwoSL+1)) * Wigner::wigner6j(TwoJ,TwoJ,2,1,1,((N1==1)?1:0)) * Wigner::wigner6j(TwoJ,TwoJ,2,TwoSL,TwoSL,TwoSR);
298 
299  double * Dblock = Dtensor->gStorage(NL,TwoSL,IL,NL,TwoSL,IL);
300  for (int cntR=0; cntR<dimR; cntR++){
301  for (int cntL=0; cntL<dimL; cntL++){
302  memHeffDiag[ptr + cntL + dimL*cntR] += alpha * Dblock[(dimL+1)*cntL];
303  }
304  }
305 
306  }
307 
308 }
309 
310 void CheMPS2::Heff::addDiagonal2e3spin1(const int ikappa, double * memHeffDiag, const Sobject * denS, TensorOperator * Dtensor) const{
311 
312  int N1 = denS->gN1(ikappa);
313 
314  if (N1==1){
315 
316  int theindex = denS->gIndex();
317  int ptr = denS->gKappa2index(ikappa);
318 
319  int NL = denS->gNL(ikappa);
320  int TwoSL = denS->gTwoSL(ikappa);
321  int IL = denS->gIL(ikappa);
322  int NR = denS->gNR(ikappa);
323  int TwoSR = denS->gTwoSR(ikappa);
324  int IR = denS->gIR(ikappa);
325  int TwoJ = denS->gTwoJ(ikappa);
326  int N2 = denS->gN2(ikappa);
327 
328  int dimL = denBK->gCurrentDim(theindex, NL,TwoSL,IL);
329  int dimR = denBK->gCurrentDim(theindex+2,NR,TwoSR,IR);
330 
331  int fase = phase(TwoSR + TwoSL + 2*TwoJ + ((N2==1)?1:0) + 1);
332  const double alpha = fase * (TwoJ+1) * sqrt(3.0*(TwoSR+1)) * Wigner::wigner6j(TwoJ,TwoJ,2,1,1,((N2==1)?1:0)) * Wigner::wigner6j(TwoJ,TwoJ,2,TwoSR,TwoSR,TwoSL);
333 
334  double * Dblock = Dtensor->gStorage(NR,TwoSR,IR,NR,TwoSR,IR);
335  for (int cntR=0; cntR<dimR; cntR++){
336  for (int cntL=0; cntL<dimL; cntL++){
337  memHeffDiag[ptr + cntL + dimL*cntR] += alpha * Dblock[(dimR+1)*cntR];
338  }
339  }
340 
341  }
342 
343 }
344 
345 void CheMPS2::Heff::addDiagonal2f3spin1(const int ikappa, double * memHeffDiag, const Sobject * denS, TensorOperator * Dtensor) const{
346 
347  int N2 = denS->gN2(ikappa);
348 
349  if (N2==1){
350 
351  int theindex = denS->gIndex();
352  int ptr = denS->gKappa2index(ikappa);
353 
354  int NL = denS->gNL(ikappa);
355  int TwoSL = denS->gTwoSL(ikappa);
356  int IL = denS->gIL(ikappa);
357  int NR = denS->gNR(ikappa);
358  int TwoSR = denS->gTwoSR(ikappa);
359  int IR = denS->gIR(ikappa);
360  int TwoJ = denS->gTwoJ(ikappa);
361  int N1 = denS->gN1(ikappa);
362 
363  int dimL = denBK->gCurrentDim(theindex, NL,TwoSL,IL);
364  int dimR = denBK->gCurrentDim(theindex+2,NR,TwoSR,IR);
365 
366  int fase = phase(TwoSR + TwoSL + 2*TwoJ + ((N1==1)?1:0) + 1);
367  const double alpha = fase * (TwoJ+1) * sqrt(3.0*(TwoSR+1)) * Wigner::wigner6j(TwoJ,TwoJ,2,1,1,((N1==1)?1:0)) * Wigner::wigner6j(TwoJ,TwoJ,2,TwoSR,TwoSR,TwoSL);
368 
369  double * Dblock = Dtensor->gStorage(NR,TwoSR,IR,NR,TwoSR,IR);
370  for (int cntR=0; cntR<dimR; cntR++){
371  for (int cntL=0; cntL<dimL; cntL++){
372  memHeffDiag[ptr + cntL + dimL*cntR] += alpha * Dblock[(dimR+1)*cntR];
373  }
374  }
375 
376  }
377 
378 }
379 
380 void CheMPS2::Heff::addDiagonal2a3spin0(const int ikappa, double * memHeffDiag, const Sobject * denS, TensorOperator **** Ctensors, TensorF0 **** F0tensors) const{
381 
382  #ifdef CHEMPS2_MPI_COMPILATION
383  const int MPIRANK = MPIchemps2::mpi_rank();
384  #endif
385 
386  int NL = denS->gNL(ikappa);
387  int TwoSL = denS->gTwoSL(ikappa);
388  int IL = denS->gIL(ikappa);
389 
390  int NR = denS->gNR(ikappa);
391  int TwoSR = denS->gTwoSR(ikappa);
392  int IR = denS->gIR(ikappa);
393 
394  int theindex = denS->gIndex();
395  int ptr = denS->gKappa2index(ikappa);
396 
397  int dimL = denBK->gCurrentDim(theindex ,NL,TwoSL,IL);
398  int dimR = denBK->gCurrentDim(theindex+2,NR,TwoSR,IR);
399 
400  bool leftSum = ( theindex < Prob->gL()*0.5 )?true:false;
401 
402  if (leftSum){
403 
404  for (int l_gamma=0; l_gamma<theindex; l_gamma++){
405  for (int l_alpha=l_gamma+1; l_alpha<theindex; l_alpha++){
406 
407  #ifdef CHEMPS2_MPI_COMPILATION
408  if ( MPIchemps2::owner_cdf( Prob->gL(), l_gamma, l_alpha ) == MPIRANK )
409  #endif
410  {
411  if (denBK->gIrrep(l_alpha) == denBK->gIrrep(l_gamma)){
412 
413  double * Cblock = Ctensors[theindex+1][l_alpha-l_gamma][theindex+1-l_alpha]->gStorage(NR,TwoSR,IR,NR,TwoSR,IR);
414  double * BlockF0 = F0tensors[theindex-1][l_alpha-l_gamma][theindex-1-l_alpha]->gStorage(NL,TwoSL,IL,NL,TwoSL,IL);
415 
416  for (int cntL=0; cntL<dimL; cntL++){
417  for (int cntR=0; cntR<dimR; cntR++){
418  memHeffDiag[ptr + cntL + dimL*cntR] += BlockF0[(dimL+1)*cntL] * Cblock[(dimR+1)*cntR];
419  }
420  }
421  }
422  }
423  }
424  }
425 
426  for (int l_alpha=0; l_alpha<theindex; l_alpha++){
427  for (int l_gamma=l_alpha; l_gamma<theindex; l_gamma++){
428 
429  #ifdef CHEMPS2_MPI_COMPILATION
430  if ( MPIchemps2::owner_cdf( Prob->gL(), l_alpha, l_gamma ) == MPIRANK )
431  #endif
432  {
433  if (denBK->gIrrep(l_alpha) == denBK->gIrrep(l_gamma)){
434 
435  double * Cblock = Ctensors[theindex+1][l_gamma-l_alpha][theindex+1-l_gamma]->gStorage(NR,TwoSR,IR,NR,TwoSR,IR);
436  double * BlockF0 = F0tensors[theindex-1][l_gamma-l_alpha][theindex-1-l_gamma]->gStorage(NL,TwoSL,IL,NL,TwoSL,IL);
437 
438  for (int cntL=0; cntL<dimL; cntL++){
439  for (int cntR=0; cntR<dimR; cntR++){
440  memHeffDiag[ptr + cntL + dimL*cntR] += BlockF0[(dimL+1)*cntL] * Cblock[(dimR+1)*cntR];
441  }
442  }
443  }
444  }
445  }
446  }
447 
448  } else {
449 
450  for (int l_delta=theindex+2; l_delta<Prob->gL(); l_delta++){
451  for (int l_beta=l_delta+1; l_beta<Prob->gL(); l_beta++){
452 
453  #ifdef CHEMPS2_MPI_COMPILATION
454  if ( MPIchemps2::owner_cdf( Prob->gL(), l_delta, l_beta ) == MPIRANK )
455  #endif
456  {
457  if (denBK->gIrrep(l_delta) == denBK->gIrrep(l_beta)){
458 
459  double * Cblock = Ctensors[theindex-1][l_beta-l_delta][l_delta-theindex]->gStorage(NL,TwoSL,IL,NL,TwoSL,IL);
460  double * BlockF0 = F0tensors[theindex+1][l_beta-l_delta][l_delta-theindex-2]->gStorage(NR,TwoSR,IR,NR,TwoSR,IR);
461 
462  for (int cntL=0; cntL<dimL; cntL++){
463  for (int cntR=0; cntR<dimR; cntR++){
464  memHeffDiag[ptr + cntL + dimL*cntR] += Cblock[(dimL+1)*cntL] * BlockF0[(dimR+1)*cntR];
465  }
466  }
467  }
468  }
469  }
470  }
471 
472  for (int l_beta=theindex+2; l_beta<Prob->gL(); l_beta++){
473  for (int l_delta=l_beta; l_delta<Prob->gL(); l_delta++){
474 
475  #ifdef CHEMPS2_MPI_COMPILATION
476  if ( MPIchemps2::owner_cdf( Prob->gL(), l_beta, l_delta ) == MPIRANK )
477  #endif
478  {
479  if (denBK->gIrrep(l_delta) == denBK->gIrrep(l_beta)){
480 
481  double * Cblock = Ctensors[theindex-1][l_delta-l_beta][l_beta-theindex]->gStorage(NL,TwoSL,IL,NL,TwoSL,IL);
482  double * BlockF0 = F0tensors[theindex+1][l_delta-l_beta][l_beta-theindex-2]->gStorage(NR,TwoSR,IR,NR,TwoSR,IR);
483 
484  for (int cntL=0; cntL<dimL; cntL++){
485  for (int cntR=0; cntR<dimR; cntR++){
486  memHeffDiag[ptr + cntL + dimL*cntR] += Cblock[(dimL+1)*cntL] * BlockF0[(dimR+1)*cntR];
487  }
488  }
489  }
490  }
491  }
492  }
493 
494  }
495 
496 }
497 
498 void CheMPS2::Heff::addDiagonal2a3spin1(const int ikappa, double * memHeffDiag, const Sobject * denS, TensorOperator **** Dtensors, TensorF1 **** F1tensors) const{
499 
500  #ifdef CHEMPS2_MPI_COMPILATION
501  const int MPIRANK = MPIchemps2::mpi_rank();
502  #endif
503 
504  int NL = denS->gNL(ikappa);
505  int TwoSL = denS->gTwoSL(ikappa);
506  int IL = denS->gIL(ikappa);
507 
508  int NR = denS->gNR(ikappa);
509  int TwoSR = denS->gTwoSR(ikappa);
510  int IR = denS->gIR(ikappa);
511 
512  int TwoJ = denS->gTwoJ(ikappa);
513  const int fase = phase(TwoSL+TwoSR+TwoJ+2);
514  const double alpha = fase * sqrt((TwoSR + 1)*(TwoSL + 1.0)) * Wigner::wigner6j(TwoSL,TwoSR,TwoJ,TwoSR,TwoSL,2);
515 
516  int theindex = denS->gIndex();
517  int ptr = denS->gKappa2index(ikappa);
518 
519  int dimL = denBK->gCurrentDim(theindex ,NL,TwoSL,IL);
520  int dimR = denBK->gCurrentDim(theindex+2,NR,TwoSR,IR);
521 
522  bool leftSum = ( theindex < Prob->gL()*0.5 )?true:false;
523 
524  if (leftSum){
525 
526  for (int l_gamma=0; l_gamma<theindex; l_gamma++){
527  for (int l_alpha=l_gamma+1; l_alpha<theindex; l_alpha++){
528 
529  #ifdef CHEMPS2_MPI_COMPILATION
530  if ( MPIchemps2::owner_cdf( Prob->gL(), l_gamma, l_alpha ) == MPIRANK )
531  #endif
532  {
533  if (denBK->gIrrep(l_alpha) == denBK->gIrrep(l_gamma)){
534 
535  double * Dblock = Dtensors[theindex+1][l_alpha-l_gamma][theindex+1-l_alpha]->gStorage(NR,TwoSR,IR,NR,TwoSR,IR);
536  double * BlockF1 = F1tensors[theindex-1][l_alpha-l_gamma][theindex-1-l_alpha]->gStorage(NL,TwoSL,IL,NL,TwoSL,IL);
537 
538  for (int cntL=0; cntL<dimL; cntL++){
539  for (int cntR=0; cntR<dimR; cntR++){
540  memHeffDiag[ptr + cntL + dimL*cntR] += alpha * BlockF1[(dimL+1)*cntL] * Dblock[(dimR+1)*cntR];
541  }
542  }
543  }
544  }
545  }
546  }
547 
548  for (int l_alpha=0; l_alpha<theindex; l_alpha++){
549  for (int l_gamma=l_alpha; l_gamma<theindex; l_gamma++){
550 
551  #ifdef CHEMPS2_MPI_COMPILATION
552  if ( MPIchemps2::owner_cdf( Prob->gL(), l_alpha, l_gamma ) == MPIRANK )
553  #endif
554  {
555 
556  if (denBK->gIrrep(l_alpha) == denBK->gIrrep(l_gamma)){
557 
558  double * Dblock = Dtensors[theindex+1][l_gamma-l_alpha][theindex+1-l_gamma]->gStorage(NR,TwoSR,IR,NR,TwoSR,IR);
559  double * BlockF1 = F1tensors[theindex-1][l_gamma-l_alpha][theindex-1-l_gamma]->gStorage(NL,TwoSL,IL,NL,TwoSL,IL);
560 
561  for (int cntL=0; cntL<dimL; cntL++){
562  for (int cntR=0; cntR<dimR; cntR++){
563  memHeffDiag[ptr + cntL + dimL*cntR] += alpha * BlockF1[(dimL+1)*cntL] * Dblock[(dimR+1)*cntR];
564  }
565  }
566  }
567  }
568  }
569  }
570 
571  } else {
572 
573  for (int l_delta=theindex+2; l_delta<Prob->gL(); l_delta++){
574  for (int l_beta=l_delta+1; l_beta<Prob->gL(); l_beta++){
575 
576  #ifdef CHEMPS2_MPI_COMPILATION
577  if ( MPIchemps2::owner_cdf( Prob->gL(), l_delta, l_beta ) == MPIRANK )
578  #endif
579  {
580  if (denBK->gIrrep(l_delta) == denBK->gIrrep(l_beta)){
581 
582  double * Dblock = Dtensors[theindex-1][l_beta-l_delta][l_delta-theindex]->gStorage(NL,TwoSL,IL,NL,TwoSL,IL);
583  double * BlockF1 = F1tensors[theindex+1][l_beta-l_delta][l_delta-theindex-2]->gStorage(NR,TwoSR,IR,NR,TwoSR,IR);
584 
585  for (int cntL=0; cntL<dimL; cntL++){
586  for (int cntR=0; cntR<dimR; cntR++){
587  memHeffDiag[ptr + cntL + dimL*cntR] += alpha * Dblock[(dimL+1)*cntL] * BlockF1[(dimR+1)*cntR];
588  }
589  }
590  }
591  }
592  }
593  }
594 
595  for (int l_beta=theindex+2; l_beta<Prob->gL(); l_beta++){
596  for (int l_delta=l_beta; l_delta<Prob->gL(); l_delta++){
597 
598  #ifdef CHEMPS2_MPI_COMPILATION
599  if ( MPIchemps2::owner_cdf( Prob->gL(), l_beta, l_delta ) == MPIRANK )
600  #endif
601  {
602  if (denBK->gIrrep(l_delta) == denBK->gIrrep(l_beta)){
603 
604  double * Dblock = Dtensors[theindex-1][l_delta-l_beta][l_beta-theindex]->gStorage(NL,TwoSL,IL,NL,TwoSL,IL);
605  double * BlockF1 = F1tensors[theindex+1][l_delta-l_beta][l_beta-theindex-2]->gStorage(NR,TwoSR,IR,NR,TwoSR,IR);
606 
607  for (int cntL=0; cntL<dimL; cntL++){
608  for (int cntR=0; cntR<dimR; cntR++){
609  memHeffDiag[ptr + cntL + dimL*cntR] += alpha * Dblock[(dimL+1)*cntL] * BlockF1[(dimR+1)*cntR];
610  }
611  }
612  }
613  }
614  }
615  }
616 
617  }
618 
619 }
620 
621 void CheMPS2::Heff::addDiagonalExcitations(const int ikappa, double * memHeffDiag, const Sobject * denS, int nLower, double ** VeffTilde) const{
622 
623  const int loc = denS->gKappa2index(ikappa);
624  const int dim = denS->gKappa2index(ikappa+1) - loc;
625  #ifdef CHEMPS2_MPI_COMPILATION
626  const int MPIRANK = MPIchemps2::mpi_rank();
627  #endif
628 
629  for (int state=0; state<nLower; state++){
630  #ifdef CHEMPS2_MPI_COMPILATION
631  if ( MPIchemps2::owner_specific_excitation( Prob->gL(), state ) == MPIRANK )
632  #endif
633  {
634  for (int cnt=0; cnt<dim; cnt++){
635  memHeffDiag[loc + cnt] += VeffTilde[state][loc+cnt] * VeffTilde[state][loc+cnt];
636  }
637  }
638  }
639 
640 }
641 
642 
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::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::MPIchemps2::mpi_rank
static int mpi_rank()
Get the rank of this MPI process.
Definition: MPIchemps2.h:131
CheMPS2::Problem::gMxElement
double gMxElement(const int alpha, const int beta, const int gamma, const int delta) const
Get a specific interaction matrix element.
Definition: Problem.cpp:350
CheMPS2::Problem::gL
int gL() const
Get the number of orbitals.
Definition: Problem.h:51
CheMPS2::Heff::phase
static int phase(const int TwoTimesPower)
Phase function.
Definition: Heff.h:75
CheMPS2::SyBookkeeper::gIrrep
int gIrrep(const int orbital) const
Get an orbital irrep.
Definition: SyBookkeeper.cpp:145