1 //
2 // ri_basis.cc
3 //
4 // Copyright (C) 2004 Edward Valeev
5 //
6 // Author: Edward Valeev <edward.valeev@chemistry.gatech.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 #ifdef __GNUC__
29 #pragma implementation
30 #endif
31
32 #include <stdexcept>
33 #include <sstream>
34
35 #include <util/misc/formio.h>
36 #include <util/misc/exenv.h>
37 #include <chemistry/qc/basis/basis.h>
38 #include <chemistry/qc/basis/symmint.h>
39 #include <chemistry/qc/mbptr12/linearr12.h>
40 #include <chemistry/qc/mbptr12/vxb_eval_info.h>
41 #include <chemistry/qc/mbptr12/svd.h>
42
43 using namespace sc;
44 using namespace std;
45
46 void
construct_ri_basis_(bool safe)47 R12IntEvalInfo::construct_ri_basis_(bool safe)
48 {
49 if (bs_aux_->equiv(bs_)) {
50 bs_ri_ = bs_;
51 if (abs_method_ == LinearR12::ABS_CABS ||
52 abs_method_ == LinearR12::ABS_CABSPlus)
53 throw std::runtime_error("R12IntEvalInfo::construct_ri_basis_ -- ABS methods CABS and CABS+ can only be used when ABS != OBS");
54 }
55 else {
56 switch(abs_method_) {
57 case LinearR12::ABS_ABS:
58 construct_ri_basis_ks_(safe);
59 break;
60 case LinearR12::ABS_ABSPlus:
61 construct_ri_basis_ksplus_(safe);
62 break;
63 case LinearR12::ABS_CABS:
64 construct_ri_basis_ev_(safe);
65 break;
66 case LinearR12::ABS_CABSPlus:
67 construct_ri_basis_evplus_(safe);
68 break;
69 default:
70 throw std::runtime_error("R12IntEvalInfo::construct_ri_basis_ -- invalid ABS method");
71 }
72 }
73 }
74
75 void
construct_ri_basis_ks_(bool safe)76 R12IntEvalInfo::construct_ri_basis_ks_(bool safe)
77 {
78 bs_ri_ = bs_aux_;
79 if (!abs_spans_obs_()) {
80 ExEnv::out0() << endl << indent << "WARNING: the auxiliary basis is not safe to use with the given orbital basis" << endl << endl;
81 if (safe)
82 throw std::runtime_error("R12IntEvalInfo::construct_ri_basis_ks_ -- auxiliary basis is not safe to use with the given orbital basis");
83 }
84 }
85
86 void
construct_ri_basis_ksplus_(bool safe)87 R12IntEvalInfo::construct_ri_basis_ksplus_(bool safe)
88 {
89 GaussianBasisSet& abs = *(bs_aux_.pointer());
90 bs_ri_ = abs + bs_;
91 construct_orthog_ri_();
92 }
93
94 void
construct_ri_basis_ev_(bool safe)95 R12IntEvalInfo::construct_ri_basis_ev_(bool safe)
96 {
97 bs_ri_ = bs_aux_;
98 if (!abs_spans_obs_()) {
99 ExEnv::out0() << endl << indent << "WARNING: the auxiliary basis is not safe to use with the given orbital basis" << endl << endl;
100 if (safe)
101 throw std::runtime_error("R12IntEvalInfo::construct_ri_basis_ev_ -- auxiliary basis is not safe to use with the given orbital basis");
102 }
103 construct_ortho_comp_svd_();
104 }
105
106 void
construct_ri_basis_evplus_(bool safe)107 R12IntEvalInfo::construct_ri_basis_evplus_(bool safe)
108 {
109 GaussianBasisSet& abs = *(bs_aux_.pointer());
110 bs_ri_ = abs + bs_;
111 construct_ortho_comp_svd_();
112 }
113
114 void
construct_orthog_aux_()115 R12IntEvalInfo::construct_orthog_aux_()
116 {
117 if (abs_space_.nonnull())
118 return;
119
120 abs_space_ = orthogonalize("ABS", bs_aux_, integral(), ref_->orthog_method(), ref_->lindep_tol(), nlindep_aux_);
121
122 if (bs_aux_ == bs_ri_)
123 ribs_space_ = abs_space_;
124 }
125
126 void
construct_orthog_vir_()127 R12IntEvalInfo::construct_orthog_vir_()
128 {
129 if (vir_space_.nonnull())
130 return;
131
132 if (bs_ == bs_vir_) {
133 // If virtuals are from the same space as occupieds, then everything is easy
134 vir_space_ = new MOIndexSpace("unoccupied MOs sorted by energy", mo_space_->coefs(),
135 mo_space_->basis(), mo_space_->integral(),
136 mo_space_->evals(), nocc_, 0);
137 // If virtuals are from the same space as occupieds, then everything is easy
138 vir_space_symblk_ = new MOIndexSpace("unoccupied MOs symmetry-blocked", mo_space_->coefs(),
139 mo_space_->basis(), mo_space_->integral(),
140 mo_space_->evals(), nocc_, 0, MOIndexSpace::symmetry);
141
142 if (nfzv_ == 0)
143 act_vir_space_ = vir_space_;
144 else
145 act_vir_space_ = new MOIndexSpace("active unoccupied MOs sorted by energy", mo_space_->coefs(),
146 mo_space_->basis(), mo_space_->integral(),
147 mo_space_->evals(), nocc_, nfzv_);
148 nlindep_vir_ = 0;
149 }
150 else {
151 // This is a set of orthonormal functions that span VBS
152 Ref<MOIndexSpace> vir_space = orthogonalize("VBS", bs_vir_, integral(), ref_->orthog_method(), ref_->lindep_tol(), nlindep_vir_);
153 // Now project out occupied MOs
154 vir_space_symblk_ = orthog_comp(occ_space_symblk_, vir_space, "VBS", ref_->lindep_tol());
155
156 // Design flaw!!! Need to compute Fock matrix right here but can't since Fock is built into R12IntEval
157 // Need to move all relevant code outside of MBPT2-R12 code
158 if (nfzv_ != 0)
159 throw std::runtime_error("R12IntEvalInfo::construct_orthog_vir_() -- nfzv_ != 0 is not allowed yet");
160 vir_space_ = vir_space_symblk_;
161 act_vir_space_ = vir_space_symblk_;
162 }
163 }
164
165 void
construct_orthog_ri_()166 R12IntEvalInfo::construct_orthog_ri_()
167 {
168 if (bs_ri_.null())
169 throw std::runtime_error("R12IntEvalInfo::construct_orthog_ri_ -- RI basis has not been set yet");
170 if (bs_aux_ == bs_ri_)
171 construct_orthog_aux_();
172 if (ribs_space_.nonnull())
173 return;
174
175 ribs_space_ = orthogonalize("RI-BS", bs_ri_, integral(), ref_->orthog_method(), ref_->lindep_tol(), nlindep_ri_);
176 }
177
178 bool
abs_spans_obs_()179 R12IntEvalInfo::abs_spans_obs_()
180 {
181 construct_orthog_aux_();
182
183 // Compute the bumber of linear dependencies in OBS+ABS
184 GaussianBasisSet& abs = *(bs_aux_.pointer());
185 Ref<GaussianBasisSet> ri_basis = abs + bs_;
186 int nlindep_ri = 0;
187 if (bs_ri_.nonnull() && ri_basis->equiv(bs_ri_)) {
188 construct_orthog_ri_();
189 nlindep_ri = nlindep_ri_;
190 }
191 else {
192 Ref<MOIndexSpace> ribs_space = orthogonalize("OBS+ABS", ri_basis, integral(), ref_->orthog_method(), ref_->lindep_tol(), nlindep_ri);
193 }
194
195 if (nlindep_ri - nlindep_aux_ - mo_space_->rank() == 0)
196 return true;
197 else
198 return false;
199 }
200
201 /////////////////////////////////////////////////////////////////////////////
202
203 void
construct_ortho_comp_svd_()204 R12IntEvalInfo::construct_ortho_comp_svd_()
205 {
206 construct_orthog_aux_();
207 construct_orthog_vir_();
208 construct_orthog_ri_();
209
210 if (debug_ > 1) {
211 occ_space_symblk_->coefs().print("Occupied MOs (symblocked)");
212 vir_space_symblk_->coefs().print("Virtual MOs (symblocked)");
213 obs_space_->coefs().print("All MOs");
214 act_occ_space_->coefs().print("Active occupied MOs");
215 act_vir_space_->coefs().print("Active virtual MOs");
216 ribs_space_->coefs().print("Orthogonal RI-BS");
217 }
218
219 ribs_space_ = orthog_comp(occ_space_symblk_, ribs_space_, "RI-BS", ref_->lindep_tol());
220 ribs_space_ = orthog_comp(vir_space_symblk_, ribs_space_, "RI-BS", ref_->lindep_tol());
221 }
222
223 Ref<MOIndexSpace>
orthogonalize(const std::string & name,const Ref<GaussianBasisSet> & bs,const Ref<Integral> & ints,OverlapOrthog::OrthogMethod orthog_method,double lindep_tol,int & nlindep)224 R12IntEvalInfo::orthogonalize(const std::string& name, const Ref<GaussianBasisSet>& bs,
225 const Ref<Integral>& ints,
226 OverlapOrthog::OrthogMethod orthog_method, double lindep_tol,
227 int& nlindep)
228 {
229 // Make an Integral and initialize with bs_aux
230 Ref<Integral> integral = ints->clone();
231 integral->set_basis(bs);
232 Ref<PetiteList> plist = integral->petite_list();
233 Ref<OneBodyInt> ov_engine = integral->overlap();
234
235 // form skeleton s matrix
236 Ref<SCMatrixKit> matrixkit = bs->matrixkit();
237 RefSymmSCMatrix s(bs->basisdim(), matrixkit);
238 Ref<SCElementOp> ov =
239 new OneBodyIntOp(new SymmOneBodyIntIter(ov_engine, plist));
240
241 s.assign(0.0);
242 s.element_op(ov);
243 ov=0;
244 //if (debug_ > 1) {
245 // std::string s_label = "AO skeleton overlap (" + name + "/" + name + ")";
246 // s.print(s_label.c_str());
247 //}
248
249 // then symmetrize it
250 RefSCDimension sodim = plist->SO_basisdim();
251 Ref<SCMatrixKit> so_matrixkit = bs->so_matrixkit();
252 RefSymmSCMatrix overlap(sodim, so_matrixkit);
253 plist->symmetrize(s,overlap);
254
255 // and clean up a bit
256 ov_engine = 0;
257 s = 0;
258
259 //
260 // Compute orthogonalizer for bs
261 //
262 ExEnv::out0() << indent << "Orthogonalizing basis for space " << name << ":" << endl << incindent;
263 OverlapOrthog orthog(orthog_method,
264 overlap,
265 so_matrixkit,
266 lindep_tol,
267 0);
268 RefSCMatrix orthog_so = orthog.basis_to_orthog_basis();
269 orthog_so = orthog_so.t();
270 RefSCMatrix orthog_ao = plist->evecs_to_AO_basis(orthog_so);
271 orthog_so = 0;
272 ExEnv::out0() << decindent;
273
274 nlindep = orthog.nlindep();
275 Ref<MOIndexSpace> space = new MOIndexSpace(name,orthog_ao,bs,integral);
276
277 return space;
278 }
279
280
281 Ref<MOIndexSpace>
orthog_comp(const Ref<MOIndexSpace> & space1,const Ref<MOIndexSpace> & space2,const std::string & name,double lindep_tol)282 R12IntEvalInfo::orthog_comp(const Ref<MOIndexSpace>& space1, const Ref<MOIndexSpace>& space2,
283 const std::string& name, double lindep_tol)
284 {
285 if (!space1->integral()->equiv(space2->integral()))
286 throw ProgrammingError("Two MOIndexSpaces use incompatible Integral factories");
287 // Both spaces must be ordered in the same way
288 if (space1->moorder() != space2->moorder())
289 throw std::runtime_error("R12IntEvalInfo::orthog_comp() -- space1 and space2 are ordered differently ");
290
291 ExEnv::out0() << indent
292 << "SVD-projecting out " << space1->name() << " out of " << space2->name()
293 << " to obtain space " << name << endl << incindent;
294
295 // C12 = C1 * S12 * C2
296 RefSCMatrix C12;
297 compute_overlap_ints(space1,space2,C12);
298
299 //
300 // SVDecompose C12 = U Sigma V and throw out columns of V
301 //
302 Ref<SCMatrixKit> ao_matrixkit = space1->basis()->matrixkit();
303 Ref<SCMatrixKit> so_matrixkit = space1->basis()->so_matrixkit();
304 int nblocks = C12.nblock();
305 const double toler = lindep_tol;
306 double min_sigma = 1.0;
307 double max_sigma = 0.0;
308 int* nvec_per_block = new int[nblocks];
309 // basis for orthogonal complement is a vector of nvecs by nbasis2
310 // we don't know nvecs yet, so use rank2
311 RefSCMatrix orthog2 = space2->coefs();
312 int rank2 = orthog2.coldim().n();
313 int nbasis2 = orthog2.rowdim().n();
314 double* vecs = new double[rank2 * nbasis2];
315 int nlindep = 0;
316
317 int v_offset = 0;
318 for(int b=0; b<nblocks; b++) {
319
320 RefSCDimension rowd = C12.rowdim()->blocks()->subdim(b);
321 RefSCDimension cold = C12.coldim()->blocks()->subdim(b);
322 int nrow = rowd.n();
323 int ncol = cold.n();
324 if (nrow && ncol) {
325
326 RefSCMatrix C12_b = C12.block(b);
327 RefSCDimension sigd = nrow < ncol ? rowd : cold;
328 int nsigmas = sigd.n();
329
330 RefSCMatrix U(rowd, rowd, ao_matrixkit);
331 RefSCMatrix V(cold, cold, ao_matrixkit);
332 RefDiagSCMatrix Sigma(sigd, ao_matrixkit);
333
334 // C12_b.svd(U,Sigma,V);
335 exp::lapack_svd(C12_b,U,Sigma,V);
336
337 // Transform V into AO basis. Vectors are in rows
338 RefSCMatrix orthog2_b = orthog2.block(b);
339 V = V * orthog2_b.t();
340
341 // Figure out how many sigmas are too small, i.e. how many vectors from space2 overlap
342 // only weakly with space1.
343 // NOTE: Sigma values returned by svd() are in descending order
344 int nzeros = 0;
345 for(int s=0; s<nsigmas; s++) {
346 double sigma = Sigma(s);
347 if (sigma < toler)
348 nzeros++;
349 if (sigma < min_sigma)
350 min_sigma = sigma;
351 if (sigma > max_sigma)
352 max_sigma = sigma;
353 }
354
355 // number of vectors that span the orthogonal space
356 nvec_per_block[b] = nzeros + ncol - nsigmas;
357 nlindep += nsigmas - nzeros;
358
359 if (nvec_per_block[b]) {
360 int v_first = nsigmas - nzeros;
361 int v_last = ncol - 1;
362 double* v_ptr = vecs + v_offset*nbasis2;
363 RefSCMatrix vtmp = V.get_subblock(v_first,v_last,0,nbasis2-1);
364 vtmp.convert(v_ptr);
365 }
366 }
367 else {
368 nvec_per_block[b] = ncol;
369
370 if (nvec_per_block[b]) {
371 RefSCMatrix orthog2_b = orthog2.block(b);
372 orthog2_b = orthog2_b.t();
373 double* v_ptr = vecs + v_offset*nbasis2;
374 orthog2_b.convert(v_ptr);
375 }
376 }
377
378 v_offset += nvec_per_block[b];
379 }
380
381 // Modify error message
382 if (v_offset == 0) {
383 const std::string errmsg = "R12IntEvalInfo::orthog_comp() -- " + space2->name()
384 + " has null projection on orthogonal complement to " + space2->name()
385 + "Modify/increase basis for " + space2->name() + ".";
386 throw std::runtime_error(errmsg.c_str());
387 }
388
389 // convert vecs into orthog2
390 // modify for the dimension
391 RefSCDimension orthog_dim = new SCDimension(v_offset, nblocks, nvec_per_block, "");
392 for(int b=0; b<nblocks; b++)
393 orthog_dim->blocks()->set_subdim(b, new SCDimension(nvec_per_block[b]));
394 RefSCMatrix orthog_vecs(orthog_dim,orthog2.rowdim(),so_matrixkit);
395 orthog_vecs.assign(vecs);
396 orthog2 = orthog_vecs.t();
397
398 ExEnv::out0() << indent
399 << nlindep << " basis function"
400 << (nlindep>1?"s":"")
401 << " projected out of " << space2->name() << "."
402 << endl;
403 ExEnv::out0() << indent
404 << "n(basis): ";
405 for (int i=0; i<orthog_dim->blocks()->nblock(); i++) {
406 ExEnv::out0() << scprintf(" %5d", orthog_dim->blocks()->size(i));
407 }
408 ExEnv::out0() << endl;
409 ExEnv::out0() << indent
410 << "Maximum singular value = "
411 << max_sigma << endl
412 << indent
413 << "Minimum singular value = "
414 << min_sigma << endl;
415 ExEnv::out0() << decindent;
416
417 delete[] vecs;
418 delete[] nvec_per_block;
419
420 Ref<MOIndexSpace> orthog_comp_space = new MOIndexSpace(name,orthog2,space2->basis(),space2->integral());
421
422 return orthog_comp_space;
423 }
424
425
426 Ref<MOIndexSpace>
gen_project(const Ref<MOIndexSpace> & space1,const Ref<MOIndexSpace> & space2,const std::string & name,double lindep_tol)427 R12IntEvalInfo::gen_project(const Ref<MOIndexSpace>& space1, const Ref<MOIndexSpace>& space2,
428 const std::string& name, double lindep_tol)
429 {
430 //
431 // Projection works as follows:
432 // 1) Compute overlap matrix between orthonormal spaces 1 and 2: C12 = C1 * S12 * C2
433 // 2) SVDecompose C12 = U Sigma V^t, throw out (near)singular triplets
434 // 3) Projected vectors (in AO basis) are X2 = C2 * V * Sigma^{-1} * U^t, where Sigma^{-1} is the generalized inverse
435 //
436
437 // Check integral factories
438 if (!space1->integral()->equiv(space2->integral()))
439 throw ProgrammingError("Two MOIndexSpaces use incompatible Integral factories");
440 // Both spaces must be ordered in the same way
441 if (space1->moorder() != space2->moorder())
442 throw std::runtime_error("R12IntEvalInfo::orthog_comp() -- space1 and space2 are ordered differently ");
443
444 ExEnv::out0() << indent
445 << "Projecting " << space1->name() << " onto " << space2->name()
446 << " exactly to obtain space " << name << endl << incindent;
447
448 // C12 = C1 * S12 * C2
449 RefSCMatrix C12;
450 compute_overlap_ints(space1,space2,C12);
451 C12.print("C12 matrix");
452
453 // Check dimensions of C12 to make sure that projection makes sense
454
455
456 Ref<SCMatrixKit> ao_matrixkit = space1->basis()->matrixkit();
457 Ref<SCMatrixKit> so_matrixkit = space1->basis()->so_matrixkit();
458 int nblocks = C12.nblock();
459 const double toler = lindep_tol;
460 double min_sigma = 1.0;
461 double max_sigma = 0.0;
462 int* nvec_per_block = new int[nblocks];
463
464 // projected vectors are a matrix of nvecs by nbasis2
465 // we don't know nvecs yet, so use rank1
466 RefSCMatrix C1 = space1->coefs();
467 RefSCMatrix C2 = space2->coefs();
468 int rank1 = space1->coefs()->ncol();
469 int nbasis2 = C2->nrow();
470 double* vecs = new double[rank1 * nbasis2];
471 int nweakovlp = 0;
472
473 int v_offset = 0;
474 for(int b=0; b<nblocks; b++) {
475
476 RefSCDimension rowd = C12.rowdim()->blocks()->subdim(b);
477 RefSCDimension cold = C12.coldim()->blocks()->subdim(b);
478 int nrow = rowd.n();
479 int ncol = cold.n();
480
481 // Cannot project if rank of the target space is smaller than the rank of the source space
482 if (nrow > ncol)
483 throw std::runtime_error("R12IntEvalInfo::svd_project() -- rank of the target space is smaller than the rank of the source space");
484
485 if (nrow && ncol) {
486
487 RefSCMatrix C12_b = C12.block(b);
488 RefSCDimension sigd = rowd;
489 int nsigmas = sigd.n();
490
491 RefSCMatrix U(rowd, rowd, ao_matrixkit);
492 RefSCMatrix V(cold, cold, ao_matrixkit);
493 RefDiagSCMatrix Sigma(sigd, ao_matrixkit);
494
495 //
496 // Compute C12 = U * Sigma * V
497 //
498 /* C12_b.svd(U,Sigma,V); */
499 exp::lapack_svd(C12_b,U,Sigma,V);
500
501 // Figure out how many sigmas are too small, i.e. how many vectors from space2 overlap
502 // only weakly with space1.
503 // NOTE: Sigma values returned by svd() are in descending order
504 int nzeros = 0;
505 for(int s=0; s<nsigmas; s++) {
506 double sigma = Sigma(s);
507 if (sigma < toler)
508 nzeros++;
509 if (sigma < min_sigma)
510 min_sigma = sigma;
511 if (sigma > max_sigma)
512 max_sigma = sigma;
513 }
514
515 // number of vectors that span the projected space
516 nvec_per_block[b] = nsigmas - nzeros;
517 if (nvec_per_block[b] < nrow)
518 throw std::runtime_error("R12IntEvalInfo::gen_project() -- space 1 is not fully spanned by space 2");
519 nweakovlp += nzeros + ncol - nrow;
520
521 if (nvec_per_block[b]) {
522 int s_first = 0;
523 int s_last = nvec_per_block[b]-1;
524 RefSCMatrix vtmp = V.get_subblock(s_first,s_last,0,ncol-1);
525 RefSCDimension rowdim = vtmp.rowdim();
526 RefDiagSCMatrix stmp = vtmp.kit()->diagmatrix(rowdim);
527 for(int i=0; i<nvec_per_block[b]; i++)
528 stmp(i) = 1.0/(Sigma(i));
529 RefSCMatrix utmp = U.get_subblock(0,nrow-1,s_first,s_last);
530 RefSCMatrix C12_inv_t = (utmp * stmp) * vtmp;
531
532 (C12_b * C12_inv_t.t()).print("C12 * C12^{-1}");
533 (C12_inv_t * C12_b.t()).print("C12^{-1} * C12");
534
535 // Transform V into AO basis and transpose so that vectors are in rows
536 RefSCMatrix C2_b = C2.block(b);
537 RefSCMatrix X2_t = C12_inv_t * C2_b.t();
538 double* x2t_ptr = vecs + v_offset*nbasis2;
539 X2_t.convert(x2t_ptr);
540 }
541 }
542 else {
543 nvec_per_block[b] = 0;
544 }
545
546
547 v_offset += nvec_per_block[b];
548 }
549
550 // convert vecs into proj
551 RefSCMatrix proj(C1.coldim(),C2.rowdim(),so_matrixkit);
552 proj.assign(vecs);
553 proj = proj.t();
554
555 ExEnv::out0() << indent
556 << nweakovlp << " basis function"
557 << (nweakovlp>1?"s":"")
558 << " in " << space2->name() << " did not overlap significantly with "
559 << space1->name() << "." << endl;
560 ExEnv::out0() << indent
561 << "n(basis): ";
562 for (int i=0; i<proj.coldim()->blocks()->nblock(); i++) {
563 ExEnv::out0() << scprintf(" %5d", proj.coldim()->blocks()->size(i));
564 }
565 ExEnv::out0() << endl;
566 ExEnv::out0() << indent
567 << "Maximum singular value = "
568 << max_sigma << endl
569 << indent
570 << "Minimum singular value = "
571 << min_sigma << endl;
572 ExEnv::out0() << decindent;
573
574 delete[] vecs;
575 delete[] nvec_per_block;
576
577 Ref<MOIndexSpace> proj_space = new MOIndexSpace(name,proj,space2->basis(),space2->integral());
578
579 RefSCMatrix S12; compute_overlap_ints(space1,proj_space,S12);
580 S12.print("Check: overlap between space1 and projected space");
581
582 return proj_space;
583 }
584
585 /////////////////////////////////////////////////////////////////////////////
586
587 // Local Variables:
588 // mode: c++
589 // c-file-style: "CLJ-CONDENSED"
590 // End:
591