MPQC  3.0.0-alpha
utils.impl.h
1 //
2 // utils.impl.h
3 //
4 // Copyright (C) 2005 Edward Valeev
5 //
6 // Author: Edward Valeev <evaleev@vt.edu>
7 // Maintainer: EV
8 //
9 // This file is part of the SC Toolkit.
10 //
11 // The SC Toolkit is free software; you can redistribute it and/or modify
12 // it under the terms of the GNU Library General Public License as published by
13 // the Free Software Foundation; either version 2, or (at your option)
14 // any later version.
15 //
16 // The SC Toolkit is distributed in the hope that it will be useful,
17 // but WITHOUT ANY WARRANTY; without even the implied warranty of
18 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19 // GNU Library General Public License for more details.
20 //
21 // You should have received a copy of the GNU Library General Public License
22 // along with the SC Toolkit; see the file COPYING.LIB. If not, write to
23 // the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
24 //
25 // The U.S. Government is granted a limited license as per AL 91-7.
26 //
27 
28 #include <util/misc/scexception.h>
29 #include <chemistry/qc/wfn/orbitalspace.h>
30 #include <math/mmisc/pairiter.h>
31 
32 #ifndef _chemistry_qc_lcao_utilsimpl_h
33 #define _chemistry_qc_lcao_utilsimpl_h
34 
35 namespace sc {
36 
37  template <PureSpinCase2 spin>
38  RefSCMatrix spinadapt(const RefSCMatrix &A,
39  const Ref<OrbitalSpace> &bra,
40  const Ref<OrbitalSpace> &ket){
41  SpatialMOPairIter_eq ij_iter(bra->rank());
42  SpatialMOPairIter_eq kl_iter(ket->rank());
43  const unsigned int brablock_size_ab = ij_iter.nij_ab();
44  const unsigned int ketblock_size_ab = kl_iter.nij_ab();
45  if (A.rowdim().n()%brablock_size_ab)
46  throw ProgrammingError("sc::spinadapt() -- row dimension is not integer multiple of bra-space rank",__FILE__,__LINE__);
47  if (A.coldim().n()%ketblock_size_ab)
48  throw ProgrammingError("sc::spinadapt() -- col dimension is not integer multiple of ket-space rank",__FILE__,__LINE__);
49  const unsigned int nbra_blocks = A.rowdim().n() / brablock_size_ab;
50  const unsigned int nket_blocks = A.coldim().n() / ketblock_size_ab;
51  const unsigned int brablock_size_aa = ij_iter.nij_aa();
52  const unsigned int ketblock_size_aa = kl_iter.nij_aa();
53  RefSCMatrix Aspinadapted(A.rowdim(),A.coldim(),A->kit());
54  const double one_divby_sqrt_two=1.0/sqrt(2.0);
55  unsigned int bra_offset_ab = 0;
56  unsigned int bra_offset_aa = 0;
57  for(int brablock=0; brablock<nbra_blocks; brablock++, bra_offset_ab += brablock_size_ab, bra_offset_aa += brablock_size_aa) {
58  for(ij_iter.start();int(ij_iter);ij_iter.next()) {
59 
60  const int ij_ab = ij_iter.ij_ab();
61 
62  unsigned int ket_offset_ab = 0;
63  unsigned int ket_offset_aa = 0;
64  for(int ketblock=0; ketblock<nket_blocks; ketblock++, ket_offset_ab += ketblock_size_ab, ket_offset_aa += ketblock_size_aa) {
65  for(kl_iter.start();int(kl_iter);kl_iter.next()) {
66 
67  const int kl_ab = kl_iter.ij_ab();
68  const int lk_ab = kl_iter.ij_ba();
69  double Aspinadapted_element;
70 
71  if(spin==Singlet){
72  double prefactor=1.0;
73  if(ij_iter.i()==ij_iter.j()) {
74  prefactor*=one_divby_sqrt_two;
75  }
76  if(kl_iter.i()==kl_iter.j()) {
77  prefactor*=one_divby_sqrt_two;
78  }
79 
80  Aspinadapted_element=prefactor*(A.get_element(ij_ab+bra_offset_ab,kl_ab+ket_offset_ab)+A.get_element(ij_ab+bra_offset_ab,lk_ab+ket_offset_ab));
81  }
82  else if(spin==Triplet){
83  Aspinadapted_element=A.get_element(ij_ab+bra_offset_ab,kl_ab+ket_offset_ab)-A.get_element(ij_ab+bra_offset_ab,lk_ab+ket_offset_ab);
84  }
85  else {
86  throw ProgrammingError("sc::spinadapt() -- PureSpinCase2 spin is not adeqate",__FILE__,__LINE__);
87  }
88 
89  Aspinadapted.set_element(ij_ab+bra_offset_ab,kl_ab+ket_offset_ab,Aspinadapted_element);
90  }
91  }
92  }
93  }
94 
95  return(Aspinadapted);
96  }
97 
98  template <bool accumulate>
99  void
100  antisymmetrize(RefSCMatrix& Aanti, const RefSCMatrix& A,
101  const Ref<OrbitalSpace>& bra1,
102  const Ref<OrbitalSpace>& bra2,
103  const Ref<OrbitalSpace>& ket1,
104  const Ref<OrbitalSpace>& ket2)
105  {
106  const bool bra1_eq_bra2 = (bra1 == bra2);
107  const bool ket1_eq_ket2 = (ket1 == ket2);
108  if (!bra1_eq_bra2 && !ket1_eq_ket2)
109  throw ProgrammingError("sc::antisymmetrize() -- the operation does not make sense: bra1!=bra2, ket1!=ket2",__FILE__,__LINE__);
110 
111  SpatialMOPairIter* ij_iter;
112  SpatialMOPairIter* kl_iter;
113  if (!bra1_eq_bra2)
114  ij_iter = new SpatialMOPairIter_neq(bra1->rank(),bra2->rank());
115  else
116  ij_iter = new SpatialMOPairIter_eq(bra1->rank());
117  if (!ket1_eq_ket2)
118  kl_iter = new SpatialMOPairIter_neq(ket1->rank(),ket2->rank());
119  else
120  kl_iter = new SpatialMOPairIter_eq(ket1->rank());
121 
122  const unsigned int brablock_size_ab = ij_iter->nij_ab();
123  const unsigned int ketblock_size_ab = kl_iter->nij_ab();
124  const unsigned int brablock_size_aa = ij_iter->nij_aa();
125  const unsigned int ketblock_size_aa = kl_iter->nij_aa();
126  if (brablock_size_ab==0 || ketblock_size_ab==0)
127  return;
128  if (A.rowdim().n()%brablock_size_ab)
129  throw ProgrammingError("sc::antisymmetrize() -- row dimension of Source is not integer multiple of bra-space rank",__FILE__,__LINE__);
130  if (A.coldim().n()%ketblock_size_ab)
131  throw ProgrammingError("sc::antisymmetrize() -- col dimension of Source is not integer multiple of ket-space rank",__FILE__,__LINE__);
132  if (Aanti.rowdim().n()%brablock_size_aa)
133  throw ProgrammingError("sc::antisymmetrize() -- row dimension of Result is not integer multiple of bra-space rank",__FILE__,__LINE__);
134  if (Aanti.coldim().n()%ketblock_size_aa)
135  throw ProgrammingError("sc::antisymmetrize() -- col dimension of Result is not integer multiple of ket-space rank",__FILE__,__LINE__);
136  const unsigned int nbra_blocks = A.rowdim().n() / brablock_size_ab;
137  const unsigned int nket_blocks = A.coldim().n() / ketblock_size_ab;
138  if (Aanti.rowdim().n() / brablock_size_aa != nbra_blocks)
139  throw ProgrammingError("sc::antisymmetrize() -- bra dimensions of Source and Result do not match",__FILE__,__LINE__);
140  if (Aanti.coldim().n() / ketblock_size_aa != nket_blocks)
141  throw ProgrammingError("sc::antisymmetrize() -- ket dimensions of Source and Result do not match",__FILE__,__LINE__);
142 
143  unsigned int bra_offset_ab = 0;
144  unsigned int bra_offset_aa = 0;
145  for(int brablock=0; brablock<nbra_blocks; brablock++, bra_offset_ab += brablock_size_ab, bra_offset_aa += brablock_size_aa) {
146  for(ij_iter->start();int(*ij_iter);ij_iter->next()) {
147 
148  const int ij_aa = ij_iter->ij_aa();
149  if (ij_aa == -1)
150  continue;
151  const int ij_ab = ij_iter->ij_ab();
152  const int ji_ab = ij_iter->ij_ba();
153 
154  unsigned int ket_offset_ab = 0;
155  unsigned int ket_offset_aa = 0;
156  for(int ketblock=0; ketblock<nket_blocks; ketblock++, ket_offset_ab += ketblock_size_ab, ket_offset_aa += ketblock_size_aa) {
157  for(kl_iter->start();int(*kl_iter);kl_iter->next()) {
158 
159  const int kl_aa = kl_iter->ij_aa();
160  if (kl_aa == -1)
161  continue;
162  const int kl_ab = kl_iter->ij_ab();
163  const int lk_ab = kl_iter->ij_ba();
164 
165  double Aanti_ijkl = (ket1_eq_ket2) ?
166  A.get_element(ij_ab+bra_offset_ab,kl_ab+ket_offset_ab) -
167  A.get_element(ij_ab+bra_offset_ab,lk_ab+ket_offset_ab)
168  :
169  A.get_element(ij_ab+bra_offset_ab,kl_ab+ket_offset_ab) -
170  A.get_element(ji_ab+bra_offset_ab,kl_ab+ket_offset_ab);
171  if (accumulate)
172  Aanti.accumulate_element(ij_aa+bra_offset_aa,kl_aa+ket_offset_aa,Aanti_ijkl);
173  else
174  Aanti.set_element(ij_aa+bra_offset_aa,kl_aa+ket_offset_aa,Aanti_ijkl);
175  }
176  }
177  }
178  }
179 
180  delete ij_iter;
181  delete kl_iter;
182  }
183 
184  template <bool accumulate>
185  void
186  antisymmetrize(RefSymmSCMatrix& Aanti, const RefSymmSCMatrix& A,
187  const Ref<OrbitalSpace>& bra1)
188  {
189  SpatialMOPairIter* ij_iter = new SpatialMOPairIter_eq(bra1->rank());
190  SpatialMOPairIter* kl_iter = new SpatialMOPairIter_eq(bra1->rank());
191 
192  const unsigned int block_size_ab = ij_iter->nij_ab();
193  const unsigned int block_size_aa = ij_iter->nij_aa();
194  if (block_size_ab==0)
195  return;
196  if (A.dim().n()%block_size_ab)
197  throw ProgrammingError("sc::antisymmetrize() -- dimension of Source is not integer multiple of (bra1->rank)^2",__FILE__,__LINE__);
198  if (Aanti.dim().n()%block_size_aa)
199  throw ProgrammingError("sc::antisymmetrize() -- dimension of Result is not integer multiple of (bra1->rank * (bra1->rank-1)/2)",__FILE__,__LINE__);
200  const unsigned int nblocks = A.dim().n() / block_size_ab;
201  if (Aanti.dim().n() / block_size_aa != nblocks)
202  throw ProgrammingError("sc::antisymmetrize() -- dimensions of Source and Result do not match",__FILE__,__LINE__);
203 
204  unsigned int bra_offset_ab = 0;
205  unsigned int bra_offset_aa = 0;
206  for(unsigned int brablock=0; brablock<nblocks; brablock++, bra_offset_ab += block_size_ab, bra_offset_aa += block_size_aa) {
207  for(ij_iter->start();int(*ij_iter);ij_iter->next()) {
208 
209  const int ij_aa = ij_iter->ij_aa();
210  if (ij_aa == -1)
211  continue;
212  const int ij_ab = ij_iter->ij_ab();
213  const int ji_ab = ij_iter->ij_ba();
214 
215  unsigned int ket_offset_ab = 0;
216  unsigned int ket_offset_aa = 0;
217  for(unsigned int ketblock=0; ketblock<nblocks; ketblock++, ket_offset_ab += block_size_ab, ket_offset_aa += block_size_aa) {
218  for(kl_iter->start();int(*kl_iter);kl_iter->next()) {
219 
220  const int kl_aa = kl_iter->ij_aa();
221  if (kl_aa == -1)
222  continue;
223  const int kl_ab = kl_iter->ij_ab();
224  const int lk_ab = kl_iter->ij_ba();
225 
226  double Aanti_ijkl = A.get_element(ij_ab+bra_offset_ab,kl_ab+ket_offset_ab) -
227  A.get_element(ij_ab+bra_offset_ab,lk_ab+ket_offset_ab);
228  if (accumulate)
229  Aanti.accumulate_element(ij_aa+bra_offset_aa,kl_aa+ket_offset_aa,Aanti_ijkl);
230  else
231  Aanti.set_element(ij_aa+bra_offset_aa,kl_aa+ket_offset_aa,Aanti_ijkl);
232  }
233  }
234  }
235  }
236 
237  delete ij_iter;
238  delete kl_iter;
239  }
240 
241  template <bool accumulate>
242  void antisymmetrize(double* Aanti, const double* A,
243  const int n) {
244  double * result_ptr = Aanti;
245  for(int r=0; r<n; ++r) {
246  int RC_offset = r*n;
247  int CR = r;
248  for(int c=0; c<r; ++c, ++result_ptr, CR+=n) {
249  const double v = A[RC_offset + c] - A[CR];
250  if (accumulate)
251  *result_ptr += v;
252  else
253  *result_ptr = v;
254  }
255  }
256  }
257 
258 
259  template <bool Accumulate>
260  void symmetrize(RefSCMatrix& Asymm, const RefSCMatrix& A,
261  const Ref<OrbitalSpace>& bra,
262  const Ref<OrbitalSpace>& ket)
263  {
264  if (A.rowdim().n() != Asymm.rowdim().n())
265  throw ProgrammingError("sc::symmetrize() -- source and target matrices have different row dimensions",__FILE__,__LINE__);
266  if (A.coldim().n() != Asymm.coldim().n())
267  throw ProgrammingError("sc::symmetrize() -- source and target matrices have different column dimensions",__FILE__,__LINE__);
268  SpatialMOPairIter_eq ij_iter(bra->rank());
269  SpatialMOPairIter_eq kl_iter(ket->rank());
270  const unsigned int brablock_size_ab = ij_iter.nij_ab();
271  const unsigned int ketblock_size_ab = kl_iter.nij_ab();
272  if (A.rowdim().n()%brablock_size_ab)
273  throw ProgrammingError("sc::symmetrize() -- row dimension is not integer multiple of bra-space rank",__FILE__,__LINE__);
274  if (A.coldim().n()%ketblock_size_ab)
275  throw ProgrammingError("sc::symmetrize() -- col dimension is not integer multiple of ket-space rank",__FILE__,__LINE__);
276  const unsigned int nbra_blocks = A.rowdim().n() / brablock_size_ab;
277  const unsigned int nket_blocks = A.coldim().n() / ketblock_size_ab;
278 
279  unsigned int bra_offset_ab = 0;
280  for(int brablock=0; brablock<nbra_blocks; brablock++, bra_offset_ab += brablock_size_ab) {
281  for(ij_iter.start();int(ij_iter);ij_iter.next()) {
282 
283  const unsigned int IJ_ab = ij_iter.ij_ab() + bra_offset_ab;
284  const unsigned int JI_ab = ij_iter.ij_ba() + bra_offset_ab;
285 
286  unsigned int ket_offset_ab = 0;
287  for(int ketblock=0; ketblock<nket_blocks; ketblock++, ket_offset_ab += ketblock_size_ab) {
288  for(kl_iter.start();int(kl_iter);kl_iter.next()) {
289 
290  const unsigned int KL_ab = kl_iter.ij_ab() + ket_offset_ab;
291  const unsigned int LK_ab = kl_iter.ij_ba() + ket_offset_ab;
292 
293  const double A_ijkl = A.get_element(IJ_ab,KL_ab);
294  const double A_ijlk = A.get_element(IJ_ab,LK_ab);
295  const double A_jikl = A.get_element(JI_ab,KL_ab);
296  const double A_jilk = A.get_element(JI_ab,LK_ab);
297 
298  {
299  const double Asymm_ijkl = 0.5 * (A_ijkl + A_jilk);
300  if (Accumulate) {
301  Asymm.accumulate_element(IJ_ab,KL_ab,Asymm_ijkl);
302  Asymm.accumulate_element(JI_ab,LK_ab,Asymm_ijkl);
303  }
304  else {
305  Asymm.set_element(IJ_ab,KL_ab,Asymm_ijkl);
306  Asymm.set_element(JI_ab,LK_ab,Asymm_ijkl);
307  }
308  }
309  {
310  const double Asymm_ijlk = 0.5 * (A_ijlk + A_jikl);
311  if (Accumulate) {
312  Asymm.accumulate_element(IJ_ab,LK_ab,Asymm_ijlk);
313  Asymm.accumulate_element(JI_ab,KL_ab,Asymm_ijlk);
314  }
315  else {
316  Asymm.set_element(IJ_ab,LK_ab,Asymm_ijlk);
317  Asymm.set_element(JI_ab,KL_ab,Asymm_ijlk);
318  }
319  }
320 
321  } // end of kl
322  } // endo f ket blocks
323  } // end of ij
324  } // end of bra blocks
325 
326  }
327 
328  template <bool Accumulate, sc::fastpairiter::PairSymm BraSymm, sc::fastpairiter::PairSymm KetSymm>
329  void symmetrize12(RefSCMatrix& Asymm, const RefSCMatrix& A,
330  const Ref<OrbitalSpace>& bra1,
331  const Ref<OrbitalSpace>& bra2,
332  const Ref<OrbitalSpace>& ket1,
333  const Ref<OrbitalSpace>& ket2)
334  {
335  // doesn't make sense if Asymm == A
336  MPQC_ASSERT(&Asymm != &A);
337  using namespace sc::fastpairiter;
338  using sc::fastpairiter::MOPairIter;
339 
340  // Detect inputs that violate semantics of this function
341  if ( (BraSymm == AntiSymm && KetSymm == AntiSymm) ||
342  (BraSymm == Symm && KetSymm == Symm) )
343  throw ProgrammingError("sc::symmetrize12() -- the matrix is already symmetric",__FILE__,__LINE__);
344  if ( (BraSymm == AntiSymm || BraSymm == Symm) && bra1 != bra2)
345  throw ProgrammingError("sc::symmetrize12() -- bra dimension is anti/symmetrized, but bra1!=bra2",__FILE__,__LINE__);
346  if ( (KetSymm == AntiSymm || KetSymm == Symm) && ket1 != ket2)
347  throw ProgrammingError("sc::symmetrize12() -- ket dimension is anti/symmetrized, but ket1!=ket2",__FILE__,__LINE__);
348 
349  if (A.rowdim().n() != Asymm.rowdim().n())
350  throw ProgrammingError("sc::symmetrize() -- source and target matrices have different row dimensions",__FILE__,__LINE__);
351  if (A.coldim().n() != Asymm.coldim().n())
352  throw ProgrammingError("sc::symmetrize() -- source and target matrices have different column dimensions",__FILE__,__LINE__);
353  MOPairIter<BraSymm> ij_iter(bra1->rank(),bra2->rank());
354  MOPairIter<KetSymm> kl_iter(ket1->rank(),ket2->rank());
355  const unsigned int brablock_size = ij_iter.nij();
356  const unsigned int ketblock_size = kl_iter.nij();
357  if (A.rowdim().n()%brablock_size)
358  throw ProgrammingError("sc::symmetrize() -- row dimension is not integer multiple of bra-space rank",__FILE__,__LINE__);
359  if (A.coldim().n()%ketblock_size)
360  throw ProgrammingError("sc::symmetrize() -- col dimension is not integer multiple of ket-space rank",__FILE__,__LINE__);
361  const unsigned int nbra_blocks = A.rowdim().n() / brablock_size;
362  const unsigned int nket_blocks = A.coldim().n() / ketblock_size;
363 
364  unsigned int bra_offset = 0;
365  for(int brablock=0; brablock<nbra_blocks; brablock++, bra_offset += brablock_size) {
366  for(ij_iter.start();int(ij_iter);ij_iter.next()) {
367 
368  const unsigned int IJ = ij_iter.ij() + bra_offset;
369  const unsigned int JI = (BraSymm == ASymm) ? ij_iter.ij(ij_iter.j(),ij_iter.i()) + bra_offset
370  : IJ;
371  const double ij_permfac = (BraSymm == AntiSymm) ? -1.0 : 1.0;
372 
373  unsigned int ket_offset = 0;
374  for(int ketblock=0; ketblock<nket_blocks; ketblock++, ket_offset += ketblock_size) {
375  for(kl_iter.start();int(kl_iter);kl_iter.next()) {
376 
377  const unsigned int KL = kl_iter.ij() + ket_offset;
378  const unsigned int LK = (KetSymm == ASymm) ? kl_iter.ij(kl_iter.j(),kl_iter.i()) + ket_offset
379  : KL;
380  const double kl_permfac = (KetSymm == AntiSymm) ? -1.0 : 1.0;
381 
382  const double A_ijkl = A.get_element(IJ,KL);
383  const double A_jilk = A.get_element(JI,LK);
384  const double Asymm_ijkl = 0.5 * (A_ijkl + ij_permfac*kl_permfac*A_jilk);
385 
386  if (Accumulate)
387  Asymm.accumulate_element(IJ,KL,Asymm_ijkl);
388  else
389  Asymm.set_element(IJ,KL,Asymm_ijkl);
390 
391  //ExEnv::out0() << IJ << " " << KL << " " << JI << " " << LK << " " << A_ijkl << " " << A_jilk << " " << Asymm_ijkl << std::endl;
392 
393  } // end of kl
394  } // end of ket blocks
395  } // end of ij
396  } // end of bra blocks
397 
398  }
399 
400  template <bool Accumulate,
401  sc::fastpairiter::PairSymm SrcBraSymm, sc::fastpairiter::PairSymm SrcKetSymm,
402  sc::fastpairiter::PairSymm DstBraSymm, sc::fastpairiter::PairSymm DstKetSymm
403  >
404  void symmetrize(RefSCMatrix& Asymm, const RefSCMatrix& A,
405  const Ref<OrbitalSpace>& bra1,
406  const Ref<OrbitalSpace>& bra2,
407  const Ref<OrbitalSpace>& ket1,
408  const Ref<OrbitalSpace>& ket2)
409  {
410  using namespace sc::fastpairiter;
411  using sc::fastpairiter::MOPairIter;
412 
413  // Detect inputs that violate semantics of this function
414  if (SrcBraSymm == DstBraSymm && SrcKetSymm == DstKetSymm)
415  throw ProgrammingError("sc::symmetrize() -- nothing to be done",__FILE__,__LINE__);
416  if ( (SrcBraSymm != ASymm && SrcKetSymm != ASymm) )
417  throw ProgrammingError("sc::symmetrize() -- either bra of ket must be asymmetric",__FILE__,__LINE__);
418  if ( (SrcBraSymm != ASymm && SrcBraSymm != DstBraSymm) )
419  throw ProgrammingError("sc::symmetrize() -- can only change bra that is asymmetric",__FILE__,__LINE__);
420  if ( (SrcKetSymm != ASymm && SrcKetSymm != DstKetSymm) )
421  throw ProgrammingError("sc::symmetrize() -- can only change ket that is asymmetric",__FILE__,__LINE__);
422  if (DstBraSymm!=ASymm && bra1 != bra2)
423  throw ProgrammingError("sc::symmetrize12() -- bra is to be anti/symmetrized, but bra1!=bra2",__FILE__,__LINE__);
424  if (DstKetSymm!=ASymm && ket1 != ket2)
425  throw ProgrammingError("sc::symmetrize12() -- ket is to be anti/symmetrized, but ket1!=ket2",__FILE__,__LINE__);
426  const bool transform_bra = DstBraSymm != ASymm && DstBraSymm != SrcBraSymm;
427  const bool transform_ket = DstKetSymm != ASymm && DstKetSymm != SrcKetSymm;
428 
429  MOPairIter<SrcBraSymm> ij_srciter(bra1->rank(),bra2->rank());
430  MOPairIter<SrcKetSymm> kl_srciter(ket1->rank(),ket2->rank());
431  MOPairIter<DstBraSymm> ij_dstiter(bra1->rank(),bra2->rank());
432  MOPairIter<DstKetSymm> kl_dstiter(ket1->rank(),ket2->rank());
433  const unsigned int brablock_size_src = ij_srciter.nij();
434  const unsigned int ketblock_size_src = kl_srciter.nij();
435  const unsigned int brablock_size_dst = ij_dstiter.nij();
436  const unsigned int ketblock_size_dst = kl_dstiter.nij();
437  if (A.rowdim().n()%brablock_size_src)
438  throw ProgrammingError("sc::symmetrize() -- row dimension of Source is not integer multiple of bra-space rank",__FILE__,__LINE__);
439  if (A.coldim().n()%ketblock_size_src)
440  throw ProgrammingError("sc::symmetrize() -- col dimension of Source is not integer multiple of ket-space rank",__FILE__,__LINE__);
441  if (Asymm.rowdim().n()%brablock_size_dst)
442  throw ProgrammingError("sc::symmetrize() -- row dimension of Result is not integer multiple of bra-space rank",__FILE__,__LINE__);
443  if (Asymm.coldim().n()%ketblock_size_dst)
444  throw ProgrammingError("sc::symmetrize() -- col dimension of Result is not integer multiple of ket-space rank",__FILE__,__LINE__);
445  const unsigned int nbra_blocks = A.rowdim().n() / brablock_size_src;
446  const unsigned int nket_blocks = A.coldim().n() / ketblock_size_src;
447  if (nbra_blocks != Asymm.rowdim().n() / brablock_size_dst)
448  throw ProgrammingError("sc::symmetrize() -- # of bra blocks in Source and Result do not match",__FILE__,__LINE__);
449  if (nket_blocks != Asymm.coldim().n() / ketblock_size_dst)
450  throw ProgrammingError("sc::symmetrize() -- # of bra blocks in Source and Result do not match",__FILE__,__LINE__);
451 
452  unsigned int bra_offset_src = 0;
453  unsigned int bra_offset_dst = 0;
454  for(int brablock=0; brablock<nbra_blocks;
455  brablock++, bra_offset_src += brablock_size_src, bra_offset_dst += brablock_size_dst) {
456  for(ij_dstiter.start();int(ij_dstiter);ij_dstiter.next()) {
457 
458  const unsigned int IJ_dst = ij_dstiter.ij() + bra_offset_dst;
459 
460  const unsigned int I = ij_srciter.i();
461  const unsigned int J = ij_srciter.j();
462  const unsigned IJ = ij_srciter.ij(I,J) + bra_offset_src;
463  unsigned int JI;
464  double ij_permfac;
465  if (transform_bra) {
466  JI = ij_srciter.ij(J,I) + bra_offset_src;
467  ij_permfac = (DstBraSymm == AntiSymm) ? -1.0 : 1.0;
468  }
469 
470  unsigned int ket_offset_src = 0;
471  unsigned int ket_offset_dst = 0;
472  for(int ketblock=0; ketblock<nket_blocks;
473  ketblock++, ket_offset_src += ketblock_size_src, ket_offset_dst += ketblock_size_dst) {
474  for(kl_dstiter.start();int(kl_dstiter);kl_dstiter.next()) {
475 
476  const unsigned int KL_dst = kl_dstiter.ij() + ket_offset_dst;
477 
478  const unsigned int K = kl_srciter.i();
479  const unsigned int L = kl_srciter.j();
480  const unsigned int KL = kl_srciter.ij(K,L) + ket_offset_src;
481  unsigned int LK;
482  double kl_permfac;
483  if (transform_ket) {
484  LK = kl_srciter.ij(L,K) + ket_offset_src;
485  kl_permfac = (DstKetSymm == AntiSymm) ? -1.0 : 1.0;
486  }
487 
488  double Asymm_ijkl;
489  if (transform_bra) {
490  const double A_ijkl = A.get_element(IJ,KL);
491  const double A_jikl = A.get_element(JI,KL);
492  Asymm_ijkl = A_ijkl + ij_permfac*A_jikl;
493  }
494  if (transform_ket) {
495  const double A_ijkl = A.get_element(IJ,KL);
496  const double A_ijlk = A.get_element(IJ,LK);
497  Asymm_ijkl = A_ijkl + kl_permfac*A_ijlk;
498  }
499  if (transform_bra && transform_ket) {
500  const double A_ijkl = A.get_element(IJ,KL);
501  const double A_jikl = A.get_element(JI,KL);
502  const double A_ijlk = A.get_element(IJ,LK);
503  const double A_jilk = A.get_element(JI,LK);
504  Asymm_ijkl = A_ijkl + ij_permfac*A_jikl + kl_permfac*A_ijlk + ij_permfac*kl_permfac*A_jilk;
505  }
506 
507  if (Accumulate)
508  Asymm.accumulate_element(IJ_dst,KL_dst,Asymm_ijkl);
509  else
510  Asymm.set_element(IJ_dst,KL_dst,Asymm_ijkl);
511 
512  } // end of kl
513  } // end of ket blocks
514  } // end of ij
515  } // end of bra blocks
516 
517  }
518 
519 }
520 
521 #endif
522 
sc::fastpairiter::MOPairIter
SpinMOPairIter iterates over pairs of spinorbitals of spin case Spin12 This class differs from other ...
Definition: pairiter.h:329
sc::MOPairIter
MOPairIter gives the ordering of orbital pairs.
Definition: pairiter.h:38
sc::RefSymmSCMatrix
The RefSymmSCMatrix class is a smart pointer to an SCSymmSCMatrix specialization.
Definition: matrix.h:265
sc::RefSCMatrix
The RefSCMatrix class is a smart pointer to an SCMatrix specialization.
Definition: matrix.h:135
sc::Ref
A template class that maintains references counts.
Definition: ref.h:361
sc::SpatialMOPairIter_neq::ij_aa
int ij_aa() const
Returns compound index ij for alpha-alpha case.
Definition: pairiter.h:266
sc::symmetrize
void symmetrize(const Ref< GPetiteList2 > &plist2, const Ref< Integral > &integral, const RefSymmSCMatrix &skel, const RefSymmSCMatrix &sym)
Uses plist2 to convert the "skeleton" matrix into the full matrix. Only applicable when the two basis...
sc::symmetrize12
void symmetrize12(RefSCMatrix &Asymm, const RefSCMatrix &A, const Ref< OrbitalSpace > &bra1, const Ref< OrbitalSpace > &bra2, const Ref< OrbitalSpace > &ket1, const Ref< OrbitalSpace > &ket2)
Generalization of the above.
Definition: utils.impl.h:329
sc::SpatialMOPairIter_neq
SpatialMOPairIter_neq gives the ordering of pairs of spatial orbitals from different spaces.
Definition: pairiter.h:204
sc::SpatialMOPairIter_neq::nij_ab
int nij_ab() const
Returns the number of functions in alpha-beta space.
Definition: pairiter.h:264
sc::SpatialMOPairIter::ij_aa
virtual int ij_aa() const =0
Returns compound index ij for alpha-alpha case.
sc::MOPairIter::start
virtual void start(const int first_ij=0)=0
Start the iteration.
sc::OrbitalSpace::rank
unsigned int rank() const
Returns the rank of the space.
sc::SpatialMOPairIter_neq::nij_aa
int nij_aa() const
Returns the number of functions in alpha-alpha space.
Definition: pairiter.h:262
sc::spinadapt
RefSCMatrix spinadapt(const RefSCMatrix &A, const Ref< OrbitalSpace > &bra, const Ref< OrbitalSpace > &ket)
Takes the 4-index quantity <ij|A|kl> and returns, depending on the value of the PureSpinCase2 spin,...
Definition: utils.impl.h:38
sc::SpatialMOPairIter::ij_ab
virtual int ij_ab() const =0
Returns compound index ij for alpha-beta case.
sc::MOPairIter::next
virtual void next()=0
Move to the next pair.
sc::SpatialMOPairIter_eq
SpatialMOPairIter_eq gives the ordering of same-spin and different-spin orbital pairs if both orbital...
Definition: pairiter.h:109
sc::SpatialMOPairIter_neq::ij_ab
int ij_ab() const
Returns compound index ij for alpha-beta case.
Definition: pairiter.h:268
sc::MOPairIter::nij
int nij() const
Returns the number of pair combinations for this iterator.
Definition: pairiter.h:77
sc::SpatialMOPairIter
SpatialMOPairIter gives the ordering of pairs of spatial orbitals.
Definition: pairiter.h:85
sc::ProgrammingError
This is thrown when a situations arises that should be impossible.
Definition: scexception.h:92
sc::SpatialMOPairIter::nij_aa
virtual int nij_aa() const =0
Returns the number of functions in alpha-alpha space.
sc::SpatialMOPairIter::ij_ba
virtual int ij_ba() const =0
Returns compound index ij for beta-alpha case.
sc::SpatialMOPairIter_neq::ij_ba
int ij_ba() const
Returns compound index ij for beta-alpha case.
Definition: pairiter.h:270
sc::antisymmetrize
void antisymmetrize(RefSCMatrix &Aanti, const RefSCMatrix &A, const Ref< OrbitalSpace > &bra, const Ref< OrbitalSpace > &ket, bool accumulate=false)
Antisymmetrizes 4-index quantity <ij|A|kl> -> <ij|A|kl> - <ij|A|lk> and saves to Aanti.
sc::SpatialMOPairIter::nij_ab
virtual int nij_ab() const =0
Returns the number of functions in alpha-beta space.
sc::SpatialMOPairIter_neq::start
void start(const int ij_offset=0)
Initialize the iterator assuming that iteration will start with pair ij_offset.
Definition: pairiter.h:248
sc
Contains all MPQC code up to version 3.
Definition: mpqcin.h:14
sc::SpatialMOPairIter_eq::nij_ab
int nij_ab() const
Returns the number of functions in alpha-beta space.
Definition: pairiter.h:191
sc::SpatialMOPairIter_neq::next
void next()
Move to the next pair.
Definition: pairiter.h:255
Generated at Sun Jan 26 2020 23:23:58 for MPQC 3.0.0-alpha using the documentation package Doxygen 1.8.16.