← Back to branch summary

~ubuntu-branches/ubuntu/oneiric/mpqc/oneiric

« back to all changes in this revision

Viewing changes to src/lib/chemistry/qc/mbptr12/ri_basis.cc

  • Committer: Bazaar Package Importer
  • Author(s): Michael Banck
  • Date: 2005-11-27 11:41:49 UTC
  • mfrom: (1.1.2 upstream)
  • Revision ID: james.westby@ubuntu.com-20051127114149-zgz9r3gk50w8ww2q
Tags: 2.3.0-1
http://bugs.debian.org/339232
* New upstream release.
* debian/rules (SONAME): Activate awk snippet for automatic so-name
  detection again, resulting in a bump to `7' and making a `c2a' for
  the C++ allocator change unnecessary; closes: #339232.
* debian/patches/00list (08_gcc-4.0_fixes): Removed, no longer needed.
* debian/rules (test): Remove workarounds, do not abort build if tests
  fail.
* debian/ref: Removed.
* debian/control.in (libsc): Added Conflict against libsc6c2.

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/lib/chemistry/qc/mbptr12/ri_basis.cc
 
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
 
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
 
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
 
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
 
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
 
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
 
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
 
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
 
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
 
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
 
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>
 
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>
 
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>
 
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: