1 //
2 // aotoso.cc --- more symmetry stuff
3 //
4 // Copyright (C) 1996 Limit Point Systems, Inc.
5 //
6 // Author: Edward Seidl <seidl@janed.com>
7 // Maintainer: LPS
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 <float.h>
29
30 #include <util/misc/formio.h>
31
32 #include <chemistry/qc/basis/basis.h>
33 #include <chemistry/qc/basis/integral.h>
34 #include <chemistry/qc/basis/shellrot.h>
35 #include <chemistry/qc/basis/petite.h>
36 #include <chemistry/qc/basis/f77sym.h>
37
38 using namespace std;
39 using namespace sc;
40
41 extern "C" {
42 void
43 F77_DGESVD(const char * JOBU, const char *JOBVT,
44 int *M, int *N, double *A, int *LDA,
45 double *S, double *U, int *LDU, double *VT, int *LDVT,
46 double *WORK, int *LWORK, int *INFO );
47 }
48
49 ////////////////////////////////////////////////////////////////////////////
50
contribution()51 contribution::contribution()
52 {
53 }
54
~contribution()55 contribution::~contribution()
56 {
57 }
58
contribution(int b,double c)59 contribution::contribution(int b, double c) : bfn(b), coef(c)
60 {
61 }
62
63 ////////////////////////////////////////////////////////////////////////////
64
SO()65 SO::SO() : len(0), length(0), cont(0)
66 {
67 }
68
SO(int l)69 SO::SO(int l) : len(0), length(0), cont(0)
70 {
71 set_length(l);
72 }
73
~SO()74 SO::~SO()
75 {
76 set_length(0);
77 }
78
79 SO&
operator =(const SO & so)80 SO::operator=(const SO& so)
81 {
82 set_length(so.length);
83 length = so.length;
84 for (int i=0; i < length; i++)
85 cont[i] = so.cont[i];
86 return *this;
87 }
88
89 void
set_length(int l)90 SO::set_length(int l)
91 {
92 len=l;
93 length=l;
94 if (cont) {
95 delete[] cont;
96 cont=0;
97 }
98
99 if (l)
100 cont = new contribution[l];
101 }
102
103 void
reset_length(int l)104 SO::reset_length(int l)
105 {
106 length=l;
107
108 if (l <= len)
109 return;
110
111 l=l+10;
112
113 contribution *newcont = new contribution[l];
114
115 if (cont) {
116 for (int i=0; i < len; i++)
117 newcont[i] = cont[i];
118
119 delete[] cont;
120 }
121
122 cont=newcont;
123 len=l;
124 }
125
126 int
equiv(const SO & so)127 SO::equiv(const SO& so)
128 {
129 int i;
130
131 if (so.length != length)
132 return 0;
133
134 double c=0;
135 for (i=0; i < length; i++) {
136 if (cont[i].bfn != so.cont[i].bfn)
137 return 0;
138 c += cont[i].coef*so.cont[i].coef;
139 }
140
141 // if the overlap == 1.0, they're equal (SO's should have been
142 // normalized by now)
143 if (fabs(fabs(c)-1.0) < 1.0e-3)
144 return 1;
145
146 return 0;
147 }
148
149 ////////////////////////////////////////////////////////////////////////////
150
SO_block()151 SO_block::SO_block() : len(0), so(0)
152 {
153 }
154
SO_block(int l)155 SO_block::SO_block(int l) : len(0), so(0)
156 {
157 set_length(l);
158 }
159
~SO_block()160 SO_block::~SO_block()
161 {
162 set_length(0);
163 }
164
165 void
set_length(int l)166 SO_block::set_length(int l)
167 {
168 len=l;
169 if (so) {
170 delete[] so;
171 so=0;
172 }
173
174 if (l)
175 so = new SO[l];
176 }
177
178 void
reset_length(int l)179 SO_block::reset_length(int l)
180 {
181 if (len == l) return;
182
183 SO *newso = new SO[l];
184
185 if (so) {
186 for (int i=0; i < len; i++)
187 newso[i] = so[i];
188
189 delete[] so;
190 }
191
192 so=newso;
193 len=l;
194 }
195
196 int
add(SO & s,int i)197 SO_block::add(SO& s, int i)
198 {
199 // first check to see if s is already here
200 for (int j=0; j < ((i < len) ? i : len); j++)
201 if (so[j].equiv(s))
202 return 0;
203
204 if (i >= len)
205 reset_length(i+1);
206 so[i] = s;
207
208 return 1;
209 }
210
211 void
print(const char * title)212 SO_block::print(const char *title)
213 {
214 int i,j;
215 ExEnv::out0() << indent << "SO block " << title << endl;
216 for (i=0; i < len; i++) {
217 ExEnv::out0() << indent << "SO " << i+1 << endl << indent;
218 for (j=0; j < so[i].length; j++)
219 ExEnv::out0() << scprintf(" %10d",so[i].cont[j].bfn);
220 ExEnv::out0() << endl << indent;
221 for (j=0; j < so[i].length; j++)
222 ExEnv::out0() << scprintf(" %10.7f",so[i].cont[j].coef);
223 ExEnv::out0() << endl;
224 }
225 }
226
227 ////////////////////////////////////////////////////////////////////////////
228
229 struct lin_comb {
230 int ns;
231 int f0;
232 int mapf0;
233 double **c;
234
lin_comblin_comb235 lin_comb(int ins, int if0, int imf0) : ns(ins), f0(if0), mapf0(imf0) {
236 int i;
237
238 c = new double*[ns];
239 for (i=0; i < ns; i++) {
240 c[i] = new double[ns];
241 memset(c[i],0,sizeof(double)*ns);
242 }
243 }
244
~lin_comblin_comb245 ~lin_comb() {
246 if (c) {
247 for (int i=0; i < ns; i++)
248 if (c[i])
249 delete[] c[i];
250 delete[] c;
251 c=0;
252 }
253 }
254
printlin_comb255 void print() const {
256 int i;
257 ExEnv::out0() << indent;
258 for (i=0; i < ns; i++)
259 ExEnv::out0() << scprintf(" %10d",mapf0+i);
260 ExEnv::out0() << endl;
261
262 for (i=0; i < ns; i++) {
263 ExEnv::out0() << indent << scprintf("%2d",f0+i);
264 for (int j=0; j < ns; j++)
265 ExEnv::out0() << scprintf(" %10.7f",c[i][j]);
266 ExEnv::out0() << endl;
267 }
268 }
269 };
270
271 ////////////////////////////////////////////////////////////////////////////
272
273 SO_block *
aotoso_info()274 PetiteList::aotoso_info()
275 {
276 int iuniq, i, j, d, ii, jj, g, s, c, ir;
277
278 GaussianBasisSet& gbs = *gbs_.pointer();
279 Molecule& mol = *gbs.molecule().pointer();
280
281 // create the character table for the point group
282 CharacterTable ct = mol.point_group()->char_table();
283 SymmetryOperation so;
284
285 if (c1_) {
286 SO_block *SOs = new SO_block[1];
287 SOs[0].set_length(gbs.nbasis());
288 for (i=0; i < gbs.nbasis(); i++) {
289 SOs[0].so[i].set_length(1);
290 SOs[0].so[i].cont[0].bfn=i;
291 SOs[0].so[i].cont[0].coef=1.0;
292 }
293 return SOs;
294 }
295
296 // ncomp is the number of symmetry blocks we have. for point groups with
297 // complex E representations, this will be cut in two eventually
298 int ncomp=0;
299 for (i=0; i < nirrep_; i++)
300 ncomp += ct.gamma(i).degeneracy();
301
302 // saoelem is the current SO in a block
303 int *saoelem = new int[ncomp];
304 memset(saoelem,0,sizeof(int)*ncomp);
305
306 int *whichir = new int[ncomp];
307 int *whichcmp = new int[ncomp];
308 for (i=ii=0; i < nirrep_; i++) {
309 for (int j=0; j < ct.gamma(i).degeneracy(); j++,ii++) {
310 whichir[ii] = i;
311 whichcmp[ii] = j;
312 }
313 }
314
315 // SOs is an array of SO_blocks which holds the redundant SO's
316 SO_block *SOs = new SO_block[ncomp];
317
318 for (i=0; i < ncomp; i++) {
319 ir = whichir[i];
320 int len = (ct.gamma(ir).complex()) ? nbf_in_ir_[ir]/2 : nbf_in_ir_[ir];
321 SOs[i].set_length(len);
322 }
323
324 // loop over all unique shells
325 for (iuniq=0; iuniq < gbs.molecule()->nunique(); iuniq++) {
326 int nequiv = gbs.molecule()->nequivalent(iuniq);
327 i = gbs.molecule()->unique(iuniq);
328 for (s=0; s < gbs.nshell_on_center(i); s++) {
329 int shell_i = gbs.shell_on_center(i,s);
330
331 // test to see if there are any high am cartesian functions in this
332 // shell. for now don't allow symmetry with cartesian functions...I
333 // just can't seem to get them working.
334 for (c=0; c < gbs(i,s).ncontraction(); c++) {
335 if (gbs(i,s).am(c) > 1 && gbs(i,s).is_cartesian(c)) {
336 if (ng_ != nirrep_) {
337 ExEnv::err0() << indent
338 << "PetiteList::aotoso: cannot yet handle"
339 << " symmetry for angular momentum >= 2\n";
340 abort();
341 }
342 }
343 }
344
345 // the functions do not mix between contractions
346 // so the contraction loop can be done outside the symmetry
347 // operation loop
348 int bfn_offset_in_shell = 0;
349 for (c=0; c < gbs(i,s).ncontraction(); c++) {
350 int nfuncuniq = gbs(i,s).nfunction(c);
351 int nfuncall = nfuncuniq * nequiv;
352
353 // allocate an array to store linear combinations of orbitals
354 // this is destroyed by the SVD routine
355 double **linorb = new double*[nfuncuniq];
356 linorb[0] = new double[nfuncuniq*nfuncall];
357 for (j=1; j<nfuncuniq; j++) {
358 linorb[j] = &linorb[j-1][nfuncall];
359 }
360
361 // a copy of linorb
362 double **linorbcop = new double*[nfuncuniq];
363 linorbcop[0] = new double[nfuncuniq*nfuncall];
364 for (j=1; j<nfuncuniq; j++) {
365 linorbcop[j] = &linorbcop[j-1][nfuncall];
366 }
367
368 // allocate an array for the SVD routine
369 double **u = new double*[nfuncuniq];
370 u[0] = new double[nfuncuniq*nfuncuniq];
371 for (j=1; j<nfuncuniq; j++) {
372 u[j] = &u[j-1][nfuncuniq];
373 }
374
375 // loop through each irrep to form the linear combination
376 // of orbitals of that symmetry
377 int irnum = 0;
378 for (ir=0; ir<ct.nirrep(); ir++) {
379 int cmplx = (ct.complex() && ct.gamma(ir).complex());
380 for (int comp=0; comp < ct.gamma(ir).degeneracy(); comp++) {
381 memset(linorb[0], 0, nfuncuniq*nfuncall*sizeof(double));
382 for (d=0; d < ct.gamma(ir).degeneracy(); d++) {
383 // if this is a point group with a complex E representation,
384 // then only do the first set of projections for E
385 if (d && cmplx) break;
386
387 // operate on each function in this contraction with each
388 // symmetry operation to form symmetrized linear combinations
389 // of orbitals
390
391 for (g=0; g<ng_; g++) {
392 // the character
393 double ccdg = ct.gamma(ir).p(comp,d,g);
394
395 so = ct.symm_operation(g);
396 int equivatom = atom_map_[i][g];
397 int atomoffset
398 = gbs.molecule()->atom_to_unique_offset(equivatom);
399
400 ShellRotation rr
401 = ints_->shell_rotation(gbs(i,s).am(c),
402 so,gbs(i,s).is_pure(c));
403
404 for (ii=0; ii < rr.dim(); ii++) {
405 for (jj=0; jj < rr.dim(); jj++) {
406 linorb[ii][atomoffset*nfuncuniq+jj] += ccdg * rr(ii,jj);
407 }
408 }
409 }
410 }
411 // find the linearly independent SO's for this irrep/component
412 memcpy(linorbcop[0], linorb[0], nfuncuniq*nfuncall*sizeof(double));
413 double *singval = new double[nfuncuniq];
414 double djunk=0; int ijunk=1;
415 int lwork = 5*nfuncall;
416 double *work = new double[lwork];
417 int info;
418 // solves At = V SIGMA Ut (since FORTRAN array layout is used)
419 F77_DGESVD("N","A",&nfuncall,&nfuncuniq,linorb[0],&nfuncall,
420 singval, &djunk, &ijunk, u[0], &nfuncuniq,
421 work, &lwork, &info);
422 // put the lin indep symm orbs into the so array
423 for (j=0; j<nfuncuniq; j++) {
424 if (singval[j] > 1.0e-6) {
425 SO tso;
426 tso.set_length(nfuncall);
427 int ll = 0, llnonzero = 0;
428 for (int k=0; k<nequiv; k++) {
429 for (int l=0; l<nfuncuniq; l++,ll++) {
430 double tmp = 0.0;
431 for (int m=0; m<nfuncuniq; m++) {
432 tmp += u[m][j] * linorbcop[m][ll] / singval[j];
433 }
434 if (fabs(tmp) > DBL_EPSILON) {
435 int equivatom = gbs.molecule()->equivalent(iuniq,k);
436 tso.cont[llnonzero].bfn
437 = l
438 + bfn_offset_in_shell
439 + gbs.shell_to_function(gbs.shell_on_center(equivatom,
440 s));
441 tso.cont[llnonzero].coef = tmp;
442 llnonzero++;
443 }
444 }
445 }
446 tso.reset_length(llnonzero);
447 if (llnonzero == 0) {
448 ExEnv::errn() << "aotoso: internal error: no bfns in SO"
449 << endl;
450 abort();
451 }
452 if (SOs[irnum+comp].add(tso,saoelem[irnum+comp])) {
453 saoelem[irnum+comp]++;
454 }
455 else {
456 ExEnv::errn() << "aotoso: internal error: "
457 << "impossible duplicate SO"
458 << endl;
459 abort();
460 }
461 }
462 }
463 delete[] singval;
464 delete[] work;
465 }
466 irnum += ct.gamma(ir).degeneracy();
467 }
468 bfn_offset_in_shell += gbs(i,s).nfunction(c);
469
470 delete[] linorb[0];
471 delete[] linorb;
472 delete[] linorbcop[0];
473 delete[] linorbcop;
474 delete[] u[0];
475 delete[] u;
476 }
477 }
478 }
479
480 // Make sure all the nodes agree on what the symmetry orbitals are.
481 // (All the above work for me > 0 is ignored.)
482 Ref<MessageGrp> grp = MessageGrp::get_default_messagegrp();
483 for (i=0; i<ncomp; i++) {
484 int len = SOs[i].len;
485 grp->bcast(len);
486 SOs[i].reset_length(len);
487 for (j=0; j<len; j++) {
488 int solen = SOs[i].so[j].length;
489 grp->bcast(solen);
490 SOs[i].so[j].reset_length(solen);
491 for (int k=0; k<solen; k++) {
492 grp->bcast(SOs[i].so[j].cont[k].bfn);
493 grp->bcast(SOs[i].so[j].cont[k].coef);
494 }
495 }
496 }
497
498 for (i=0; i < ncomp; i++) {
499 ir = whichir[i];
500 int scal = ct.gamma(ir).complex() ? 2 : 1;
501
502 if (saoelem[i] < nbf_in_ir_[ir]/scal) {
503 // if we found too few, there are big problems
504
505 ExEnv::err0() << indent
506 << scprintf("trouble making SO's for irrep %s\n",
507 ct.gamma(ir).symbol());
508 ExEnv::err0() << indent
509 << scprintf(" only found %d out of %d SO's\n",
510 saoelem[i], nbf_in_ir_[ir]/scal);
511 SOs[i].print("");
512 abort();
513
514 } else if (saoelem[i] > nbf_in_ir_[ir]/scal) {
515 // there are some redundant so's left...need to do something to get
516 // the elements we want
517
518 ExEnv::err0() << indent
519 << scprintf("trouble making SO's for irrep %s\n",
520 ct.gamma(ir).symbol());
521 ExEnv::err0() << indent
522 << scprintf(" found %d SO's, but there should only be %d\n",
523 saoelem[i], nbf_in_ir_[ir]/scal);
524 SOs[i].print("");
525 abort();
526 }
527 }
528
529 if (ct.complex()) {
530 SO_block *nSOs = new SO_block[nblocks_];
531
532 int in=0;
533
534 for (i=ir=0; ir < nirrep_; ir++) {
535 if (ct.gamma(ir).complex()) {
536 nSOs[in].set_length(nbf_in_ir_[ir]);
537 int k;
538 for (k=0; k < SOs[i].len; k++)
539 nSOs[in].add(SOs[i].so[k],k);
540 i++;
541
542 for (j=0; j < SOs[i].len; j++,k++)
543 nSOs[in].add(SOs[i].so[j],k);
544
545 i++;
546 in++;
547 } else {
548 for (j=0; j < ct.gamma(ir).degeneracy(); j++,i++,in++) {
549 nSOs[in].set_length(nbf_in_ir_[ir]);
550 for (int k=0; k < SOs[i].len; k++)
551 nSOs[in].add(SOs[i].so[k],k);
552 }
553 }
554 }
555
556 SO_block *tmp= SOs;
557 SOs = nSOs;
558 delete[] tmp;
559 }
560
561 delete[] saoelem;
562 delete[] whichir;
563 delete[] whichcmp;
564
565 return SOs;
566 }
567
568 RefSCMatrix
aotoso()569 PetiteList::aotoso()
570 {
571 RefSCMatrix aoso(AO_basisdim(), SO_basisdim(), gbs_->so_matrixkit());
572 aoso.assign(0.0);
573
574 if (c1_) {
575 aoso->unit();
576 return aoso;
577 }
578
579 SO_block *sos = aotoso_info();
580
581 BlockedSCMatrix *aosop = dynamic_cast<BlockedSCMatrix*>(aoso.pointer());
582
583 for (int b=0; b < aosop->nblocks(); b++) {
584 RefSCMatrix aosb = aosop->block(b);
585
586 if (aosb.null())
587 continue;
588
589 SO_block& sob = sos[b];
590
591 Ref<SCMatrixSubblockIter> iter =
592 aosb->local_blocks(SCMatrixSubblockIter::Write);
593
594 for (iter->begin(); iter->ready(); iter->next()) {
595 if (dynamic_cast<SCMatrixRectBlock*>(iter->block())) {
596 SCMatrixRectBlock *blk = dynamic_cast<SCMatrixRectBlock*>(iter->block());
597
598 int jlen = blk->jend-blk->jstart;
599
600 for (int j=0; j < sob.len; j++) {
601 if (j < blk->jstart || j >= blk->jend)
602 continue;
603
604 SO& soj = sob.so[j];
605
606 for (int i=0; i < soj.len; i++) {
607 int ii=soj.cont[i].bfn;
608
609 if (ii < blk->istart || ii >= blk->iend)
610 continue;
611
612 blk->data[(ii-blk->istart)*jlen+(j-blk->jstart)] =
613 soj.cont[i].coef;
614 }
615 }
616 } else {
617 SCMatrixRectSubBlock *blk =
618 dynamic_cast<SCMatrixRectSubBlock*>(iter->block());
619
620 for (int j=0; j < sob.len; j++) {
621 if (j < blk->jstart || j >= blk->jend)
622 continue;
623
624 SO& soj = sob.so[j];
625
626 for (int i=0; i < soj.len; i++) {
627 int ii=soj.cont[i].bfn;
628
629 if (ii < blk->istart || ii >= blk->iend)
630 continue;
631
632 blk->data[ii*blk->istride+j] = soj.cont[i].coef;
633 }
634 }
635 }
636 }
637 }
638
639 delete[] sos;
640
641 return aoso;
642 }
643
644 RefSCMatrix
sotoao()645 PetiteList::sotoao()
646 {
647 if (c1_)
648 return aotoso();
649 else if (nirrep_ == ng_) // subgroup of d2h
650 return aotoso().t();
651 else
652 return aotoso().i();
653 }
654
655 RefSymmSCMatrix
to_SO_basis(const RefSymmSCMatrix & a)656 PetiteList::to_SO_basis(const RefSymmSCMatrix& a)
657 {
658 // SO basis is always blocked, so first make sure a is blocked
659 RefSymmSCMatrix aomatrix = dynamic_cast<BlockedSymmSCMatrix*>(a.pointer());
660 if (aomatrix.null()) {
661 aomatrix = gbs_->so_matrixkit()->symmmatrix(AO_basisdim());
662 aomatrix->convert(a);
663 }
664
665 // if C1, then do nothing
666 if (c1_)
667 return aomatrix.copy();
668
669 RefSymmSCMatrix somatrix(SO_basisdim(), gbs_->so_matrixkit());
670 somatrix.assign(0.0);
671 somatrix->accumulate_transform(aotoso(), aomatrix,
672 SCMatrix::TransposeTransform);
673
674 return somatrix;
675 }
676
677 RefSymmSCMatrix
to_AO_basis(const RefSymmSCMatrix & somatrix)678 PetiteList::to_AO_basis(const RefSymmSCMatrix& somatrix)
679 {
680 // if C1, then do nothing
681 if (c1_)
682 return somatrix.copy();
683
684 RefSymmSCMatrix aomatrix(AO_basisdim(), gbs_->so_matrixkit());
685 aomatrix.assign(0.0);
686
687 if (nirrep_ == ng_) // subgroup of d2h
688 aomatrix->accumulate_transform(aotoso(), somatrix);
689 else
690 aomatrix->accumulate_transform(aotoso().i(), somatrix,
691 SCMatrix::TransposeTransform);
692
693 return aomatrix;
694 }
695
696 RefSCMatrix
evecs_to_SO_basis(const RefSCMatrix & aoev)697 PetiteList::evecs_to_SO_basis(const RefSCMatrix& aoev)
698 {
699 ExEnv::err0() << indent
700 << "PetiteList::evecs_to_SO_basis: don't work yet\n";
701 abort();
702
703 RefSCMatrix aoevecs = dynamic_cast<BlockedSCMatrix*>(aoev.pointer());
704 if (aoevecs.null()) {
705 aoevecs = gbs_->so_matrixkit()->matrix(AO_basisdim(), AO_basisdim());
706 aoevecs->convert(aoev);
707 }
708
709 RefSCMatrix soev = aotoso().t() * aoevecs;
710 soev.print("soev");
711
712 RefSCMatrix soevecs(SO_basisdim(), SO_basisdim(), gbs_->so_matrixkit());
713 soevecs->convert(soev);
714
715 return soevecs;
716 }
717
718 RefSCMatrix
evecs_to_AO_basis(const RefSCMatrix & soevecs)719 PetiteList::evecs_to_AO_basis(const RefSCMatrix& soevecs)
720 {
721 // if C1, then do nothing
722 if (c1_)
723 return soevecs.copy();
724
725 RefSCMatrix aoev = aotoso() * soevecs;
726
727 return aoev;
728 }
729
730 /////////////////////////////////////////////////////////////////////////////
731
732 void
symmetrize(const RefSymmSCMatrix & skel,const RefSymmSCMatrix & sym)733 PetiteList::symmetrize(const RefSymmSCMatrix& skel,
734 const RefSymmSCMatrix& sym)
735 {
736 GaussianBasisSet& gbs = *gbs_.pointer();
737
738 // SO basis is always blocked, so first make sure skel is blocked
739 RefSymmSCMatrix bskel = dynamic_cast<BlockedSymmSCMatrix*>(skel.pointer());
740 if (bskel.null()) {
741 bskel = gbs.so_matrixkit()->symmmatrix(AO_basisdim());
742 bskel->convert(skel);
743 }
744
745 // if C1, then do nothing
746 if (c1_) {
747 sym.assign(bskel);
748 return;
749 }
750
751 int b,c;
752
753 CharacterTable ct = gbs.molecule()->point_group()->char_table();
754
755 RefSCMatrix aoso = aotoso();
756 BlockedSCMatrix *lu = dynamic_cast<BlockedSCMatrix*>(aoso.pointer());
757
758 for (b=0; b < lu->nblocks(); b++) {
759 if (lu->block(b).null())
760 continue;
761
762 int ir = ct.which_irrep(b);
763
764 double skal = (double)ct.order()/(double)ct.gamma(ir).degeneracy();
765 skal = sqrt(skal);
766 lu->block(b).scale(skal);
767 }
768
769 sym.assign(0.0);
770 sym.accumulate_transform(aoso,bskel,SCMatrix::TransposeTransform);
771 aoso=0;
772
773 BlockedSymmSCMatrix *la = dynamic_cast<BlockedSymmSCMatrix*>(sym.pointer());
774
775 // loop through blocks and finish symmetrizing degenerate blocks
776 for (b=0; b < la->nblocks(); b++) {
777 if (la->block(b).null())
778 continue;
779
780 int ir=ct.which_irrep(b);
781
782 if (ct.gamma(ir).degeneracy()==1)
783 continue;
784
785 if (ct.gamma(ir).complex()) {
786 int nbf = nbf_in_ir_[ir]/2;
787
788 RefSymmSCMatrix irrep = la->block(b).get_subblock(0, nbf-1);
789 irrep.accumulate(la->block(b).get_subblock(nbf, 2*nbf-1));
790
791 RefSCMatrix sub = la->block(b).get_subblock(nbf, 2*nbf-1, 0, nbf-1);
792 RefSCMatrix subt = sub.t();
793 subt.scale(-1.0);
794 sub.accumulate(subt);
795 subt=0;
796
797 la->block(b).assign_subblock(irrep, 0, nbf-1);
798 la->block(b).assign_subblock(irrep,nbf, 2*nbf-1);
799 la->block(b).assign_subblock(sub, nbf, 2*nbf-1, 0, nbf-1);
800
801 } else {
802 RefSymmSCMatrix irrep = la->block(b).copy();
803 for (c=1; c < ct.gamma(ir).degeneracy(); c++)
804 irrep.accumulate(la->block(b+c));
805
806 for (c=0; c < ct.gamma(ir).degeneracy(); c++)
807 la->block(b+c).assign(irrep);
808
809 b += ct.gamma(ir).degeneracy()-1;
810 }
811 }
812 }
813
814 /////////////////////////////////////////////////////////////////////////////
815
816 // Local Variables:
817 // mode: c++
818 // c-file-style: "ETS"
819 // End:
820