source: ThirdParty/mpqc_open/src/lib/chemistry/qc/basis/orthog.cc@ 1ca493a

Candidate_v1.7.0 stable
Last change on this file since 1ca493a was 860145, checked in by Frederik Heber <heber@…>, 9 years ago

Merge commit '0b990dfaa8c6007a996d030163a25f7f5fc8a7e7' as 'ThirdParty/mpqc_open'
  • Property mode set to 100644
File size: 15.4 KB
Line 
1//
2// orthog.cc
3//
4// Copyright (C) 1996 Limit Point Systems, Inc.
5//
6// Author: Curtis Janssen <cljanss@ca.sandia.gov>
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#ifdef __GNUC__
29#pragma implementation
30#endif
31
32#include <stdlib.h>
33#include <math.h>
34
35#include <iostream>
36#include <stdexcept>
37
38#include <util/state/stateio.h>
39#include <chemistry/qc/basis/orthog.h>
40#include <math/scmat/elemop.h>
41#include <math/scmat/blocked.h>
42
43using namespace std;
44using namespace sc;
45
46static ClassDesc OverlapOrthog_cd(typeid(OverlapOrthog),"OverlapOrthog",1,
47 "virtual public SavableState",
48 0, 0, create<OverlapOrthog>);
49
50OverlapOrthog::OverlapOrthog(
51 OrthogMethod method,
52 const RefSymmSCMatrix &overlap,
53 const Ref<SCMatrixKit> &result_kit,
54 double lindep_tolerance,
55 int debug) : nlindep_(0), min_orthog_res_(0.0), max_orthog_res_(0.0)
56{
57 reinit(method,overlap,result_kit,lindep_tolerance,debug);
58}
59
60OverlapOrthog::OverlapOrthog(StateIn& si):
61 SavableState(si)
62{
63 Ref<SCMatrixKit> kit;
64 kit << si.override()->describedclassvalue("matrixkit");
65
66 if (kit.null()) {
67 throw std::runtime_error("OverlapOrthog::OverlapOrthog(StateIn& si): requires that a matrixkit be set up in the override info");
68 }
69
70 si.get(debug_);
71 dim_ << SavableState::restore_state(si);
72 orthog_dim_ << SavableState::restore_state(si);
73 si.get(lindep_tol_);
74 int i_orthog_method;
75 si.get(i_orthog_method);
76 orthog_method_ = OrthogMethod(i_orthog_method);
77
78 orthog_trans_ = kit->matrix(orthog_dim_, dim_);
79 orthog_trans_.restore(si);
80 orthog_trans_inverse_ = kit->matrix(dim_, orthog_dim_);
81 orthog_trans_inverse_.restore(si);
82 si.get(min_orthog_res_);
83 si.get(max_orthog_res_);
84 si.get(nlindep_);
85}
86
87OverlapOrthog::~OverlapOrthog()
88{
89}
90
91void
92OverlapOrthog::save_data_state(StateOut& so)
93{
94 so.put(debug_);
95 SavableState::save_state(dim_.pointer(), so);
96 SavableState::save_state(orthog_dim_.pointer(), so);
97 so.put(lindep_tol_);
98 so.put(int(orthog_method_));
99
100 orthog_trans_.save(so);
101 orthog_trans_inverse_.save(so);
102 so.put(min_orthog_res_);
103 so.put(max_orthog_res_);
104 so.put(nlindep_);
105 // The overlap_ member is not saved since should not be needed.
106 // The result_kit_ member is not saved since it depends on the
107 // runtime environment. It is given to the StateIn CTOR by
108 // an overriding KeyVal.
109}
110
111void
112OverlapOrthog::reinit(
113 OrthogMethod method,
114 const RefSymmSCMatrix &overlap,
115 const Ref<SCMatrixKit> &result_kit,
116 double lindep_tolerance,
117 int debug)
118{
119 orthog_method_ = method;
120 overlap_ = overlap;
121 lindep_tol_ = lindep_tolerance;
122 debug_ = debug;
123 dim_ = overlap_.dim();
124 result_kit_ = result_kit;
125}
126
127Ref<OverlapOrthog>
128OverlapOrthog::copy() const
129{
130 Ref<OverlapOrthog> orthog
131 = new OverlapOrthog(orthog_method_,
132 overlap_,
133 result_kit_,
134 lindep_tol_,
135 debug_);
136
137 orthog->orthog_trans_ = orthog_trans_.copy();
138 orthog->orthog_trans_inverse_ = orthog_trans_inverse_.copy();
139 orthog->orthog_dim_ = orthog_dim_;
140 orthog->min_orthog_res_ = min_orthog_res_;
141 orthog->max_orthog_res_ = max_orthog_res_;
142 orthog->nlindep_ = nlindep_;
143 return orthog;
144}
145
146// computes intermediates needed to form orthogonalization matrices
147// and their inverses.
148void
149OverlapOrthog::compute_overlap_eig(RefSCMatrix& overlap_eigvec,
150 RefDiagSCMatrix& overlap_isqrt_eigval,
151 RefDiagSCMatrix& overlap_sqrt_eigval)
152{
153 // first calculate S
154 RefSymmSCMatrix M = overlap_;
155
156 // Diagonalize M to get m and U
157 RefSCMatrix U(M.dim(), M.dim(), M.kit());
158 RefDiagSCMatrix m(M.dim(), M.kit());
159 M.diagonalize(m,U);
160 M = 0;
161
162 Ref<SCElementMaxAbs> maxabsop = new SCElementMaxAbs;
163 m.element_op(maxabsop.pointer());
164 double maxabs = maxabsop->result();
165 double s_tol = lindep_tol_ * maxabs;
166
167 double minabs = maxabs;
168 BlockedDiagSCMatrix *bm = dynamic_cast<BlockedDiagSCMatrix*>(m.pointer());
169 bool blocked;
170 int nblocks;
171 if (bm == 0) {
172 blocked = false;
173 nblocks = 1;
174 }
175 else {
176 blocked = true;
177 nblocks = bm->nblocks();
178 }
179 int i, j;
180 double *pm_sqrt = new double[m->dim()->n()];
181 double *pm_isqrt = new double[m->dim()->n()];
182 int *pm_index = new int[m->dim()->n()];
183 int *nfunc = new int[nblocks];
184 int nfunctot = 0;
185 nlindep_ = 0;
186 for (i=0; i<nblocks; i++) {
187 nfunc[i] = 0;
188 if (blocked && bm->block(i).null()) continue;
189 int n;
190 if (blocked) n = bm->block(i)->dim()->n();
191 else n = m->dim()->n();
192 double *pm = new double[n];
193 if (blocked) bm->block(i)->convert(pm);
194 else m->convert(pm);
195 for (j=0; j<n; j++) {
196 if (pm[j] > s_tol) {
197 if (pm[j] < minabs) { minabs = pm[j]; }
198 pm_sqrt[nfunctot] = sqrt(pm[j]);
199 pm_isqrt[nfunctot] = 1.0/pm_sqrt[nfunctot];
200 pm_index[nfunctot] = j;
201 nfunc[i]++;
202 nfunctot++;
203 }
204 else if (orthog_method_ == Symmetric) {
205 pm_sqrt[nfunctot] = 0.0;
206 pm_isqrt[nfunctot] = 0.0;
207 pm_index[nfunctot] = j;
208 nfunc[i]++;
209 nfunctot++;
210 nlindep_++;
211 }
212 else {
213 nlindep_++;
214 }
215 }
216 delete[] pm;
217 }
218
219 if (nlindep_ > 0 && orthog_method_ == Symmetric) {
220 ExEnv::out0() << indent
221 << "WARNING: " << nlindep_
222 << " basis function"
223 << (dim_.n()-orthog_dim_.n()>1?"s":"")
224 << " ignored in symmetric orthogonalization."
225 << endl;
226 }
227
228 // make sure all nodes end up with exactly the same data
229 MessageGrp::get_default_messagegrp()->bcast(nfunctot);
230 MessageGrp::get_default_messagegrp()->bcast(nfunc, nblocks);
231 MessageGrp::get_default_messagegrp()->bcast(pm_sqrt,nfunctot);
232 MessageGrp::get_default_messagegrp()->bcast(pm_isqrt,nfunctot);
233 MessageGrp::get_default_messagegrp()->bcast(pm_index,nfunctot);
234
235 if (orthog_method_ == Symmetric) {
236 orthog_dim_ = new SCDimension(m->dim()->blocks(),
237 "ortho basis (symmetric)");
238 }
239 else {
240 orthog_dim_ = new SCDimension(nfunctot, nblocks,
241 nfunc, "ortho basis (canonical)");
242 for (i=0; i<nblocks; i++) {
243 orthog_dim_->blocks()->set_subdim(i, new SCDimension(nfunc[i]));
244 }
245 }
246
247 overlap_eigvec = result_kit_->matrix(dim_, orthog_dim_);
248 if (orthog_method_ == Symmetric) {
249 overlap_eigvec.assign(U);
250 }
251 else {
252 BlockedSCMatrix *bev
253 = dynamic_cast<BlockedSCMatrix*>(overlap_eigvec.pointer());
254 BlockedSCMatrix *bU
255 = dynamic_cast<BlockedSCMatrix*>(U.pointer());
256 int ifunc = 0;
257 for (i=0; i<bev->nblocks(); i++) {
258 if (bev->block(i).null()) continue;
259 for (j=0; j<nfunc[i]; j++) {
260 RefSCVector col = bU->block(i)->get_column(pm_index[ifunc]);
261 bev->block(i)->assign_column(col,j);
262 col = 0;
263 ifunc++;
264 }
265 }
266 }
267
268 overlap_sqrt_eigval = result_kit_->diagmatrix(orthog_dim_);
269 overlap_sqrt_eigval->assign(pm_sqrt);
270 overlap_isqrt_eigval = result_kit_->diagmatrix(orthog_dim_);
271 overlap_isqrt_eigval->assign(pm_isqrt);
272
273 delete[] nfunc;
274 delete[] pm_sqrt;
275 delete[] pm_isqrt;
276 delete[] pm_index;
277
278 max_orthog_res_ = maxabs;
279 min_orthog_res_ = minabs;
280
281 if (debug_ > 1) {
282 overlap_.print("S");
283 overlap_eigvec.print("S eigvec");
284 overlap_isqrt_eigval.print("s^(-1/2) eigval");
285 overlap_sqrt_eigval.print("s^(1/2) eigval");
286 }
287}
288
289void
290OverlapOrthog::compute_symmetric_orthog()
291{
292 RefSCMatrix overlap_eigvec;
293 RefDiagSCMatrix overlap_isqrt_eigval;
294 RefDiagSCMatrix overlap_sqrt_eigval;
295 compute_overlap_eig(overlap_eigvec,
296 overlap_isqrt_eigval,
297 overlap_sqrt_eigval);
298
299 orthog_trans_ = overlap_eigvec
300 * overlap_isqrt_eigval
301 * overlap_eigvec.t();
302 orthog_trans_inverse_ = overlap_eigvec
303 * overlap_sqrt_eigval
304 * overlap_eigvec.t();
305}
306
307void
308OverlapOrthog::compute_canonical_orthog()
309{
310 RefSCMatrix overlap_eigvec;
311 RefDiagSCMatrix overlap_isqrt_eigval;
312 RefDiagSCMatrix overlap_sqrt_eigval;
313 compute_overlap_eig(overlap_eigvec,
314 overlap_isqrt_eigval,
315 overlap_sqrt_eigval);
316
317 orthog_trans_ = overlap_isqrt_eigval * overlap_eigvec.t();
318 orthog_trans_inverse_ = overlap_eigvec * overlap_sqrt_eigval;
319}
320
321void
322OverlapOrthog::compute_gs_orthog()
323{
324 // Orthogonalize each subblock of the overlap.
325 max_orthog_res_ = 1.0;
326 min_orthog_res_ = 1.0;
327 nlindep_ = 0;
328 BlockedSymmSCMatrix *S
329 = dynamic_cast<BlockedSymmSCMatrix *>(overlap_.pointer());
330 int nblock = S->nblocks();
331 Ref<BlockedSCMatrixKit> kit
332 = dynamic_cast<BlockedSCMatrixKit*>(S->kit().pointer());
333 Ref<SCMatrixKit> subkit = kit->subkit();
334 RefSCMatrix *blockorthogs = new RefSCMatrix[nblock];
335 int *nblockorthogs = new int[nblock];
336 int northog = 0;
337 for (int i=0; i<nblock; i++) {
338 RefSymmSCMatrix Sblock = S->block(i);
339 if (Sblock.null()) {
340 blockorthogs[i] = 0;
341 nblockorthogs[i] = 0;
342 continue;
343 }
344 RefSCDimension dim = Sblock->dim();
345 RefSCMatrix blockorthog(dim,dim,subkit);
346 blockorthog->unit();
347 double res;
348 int nblockorthog = blockorthog->schmidt_orthog_tol(Sblock, lindep_tol_,
349 &res);
350 if (res < min_orthog_res_) min_orthog_res_ = res;
351 blockorthogs[i] = blockorthog;
352 nblockorthogs[i] = nblockorthog;
353 northog += nblockorthog;
354 nlindep_ += dim.n() - nblockorthog;
355 }
356
357 // Construct the orthog basis SCDimension object.
358 Ref<SCBlockInfo> blockinfo
359 = new SCBlockInfo(northog, nblock, nblockorthogs);
360 for (int i=0; i<nblock; i++) {
361 blockinfo->set_subdim(i, new SCDimension(nblockorthogs[i]));
362 }
363 orthog_dim_ = new SCDimension(blockinfo, "ortho (Gram-Schmidt)");
364
365 // Replace each blockorthog by a matrix with only linear independent columns
366 for (int i=0; i<nblock; i++) {
367 if (nblockorthogs[i] == 0) continue;
368 RefSCMatrix old_blockorthog = blockorthogs[i];
369 blockorthogs[i] = subkit->matrix(dim_->blocks()->subdim(i),
370 orthog_dim_->blocks()->subdim(i));
371 blockorthogs[i].assign_subblock(old_blockorthog,
372 0, dim_->blocks()->subdim(i).n()-1,
373 0, orthog_dim_->blocks()->subdim(i).n()-1);
374 }
375
376 // Compute the inverse of each orthogonalization block.
377 RefSCMatrix *inverse_blockorthogs = new RefSCMatrix[nblock];
378 for (int i=0; i<nblock; i++) {
379 if (nblockorthogs[i] == 0) {
380 inverse_blockorthogs[i] = 0;
381 }
382 else {
383 inverse_blockorthogs[i] = blockorthogs[i].gi();
384 }
385 }
386
387 // Construct the complete transformation matrices
388 orthog_trans_ = result_kit_->matrix(dim_, orthog_dim_);
389 orthog_trans_inverse_ = result_kit_->matrix(orthog_dim_, dim_);
390 orthog_trans_.assign(0.0);
391 orthog_trans_inverse_.assign(0.0);
392 BlockedSCMatrix *X
393 = dynamic_cast<BlockedSCMatrix*>(orthog_trans_.pointer());
394 BlockedSCMatrix *Xi
395 = dynamic_cast<BlockedSCMatrix*>(orthog_trans_inverse_.pointer());
396 for (int i=0; i<nblock; i++) {
397 if (nblockorthogs[i] == 0) continue;
398 int nrow = blockorthogs[i].rowdim().n();
399 int ncol = blockorthogs[i].coldim().n();
400 X->block(i).assign_subblock(blockorthogs[i],
401 0, nrow-1, 0, ncol-1,
402 0, 0);
403 Xi->block(i).assign_subblock(inverse_blockorthogs[i],
404 0, ncol-1, 0, nrow-1,
405 0, 0);
406 }
407 orthog_trans_ = orthog_trans_.t();
408 orthog_trans_inverse_ = orthog_trans_inverse_.t();
409
410 delete[] blockorthogs;
411 delete[] inverse_blockorthogs;
412 delete[] nblockorthogs;
413}
414
415void
416OverlapOrthog::compute_orthog_trans()
417{
418 switch(orthog_method_) {
419 case GramSchmidt:
420 ExEnv::out0() << indent
421 << "Using Gram-Schmidt orthogonalization."
422 << endl;
423 compute_gs_orthog();
424 break;
425 case Symmetric:
426 compute_symmetric_orthog();
427 ExEnv::out0() << indent
428 << "Using symmetric orthogonalization."
429 << endl;
430 break;
431 case Canonical:
432 compute_canonical_orthog();
433 ExEnv::out0() << indent
434 << "Using canonical orthogonalization."
435 << endl;
436 break;
437 default:
438 ExEnv::outn() << "OverlapOrthog::compute_orthog_trans(): bad orthog method"
439 << endl;
440 abort();
441 }
442
443 ExEnv::out0() << indent
444 << "n(basis): ";
445 for (int i=0; i<dim_->blocks()->nblock(); i++) {
446 ExEnv::out0() << scprintf(" %5d", dim_->blocks()->size(i));
447 }
448 ExEnv::out0() << endl;
449
450 if (dim_.n() != orthog_dim_.n()) {
451 ExEnv::out0() << indent
452 << "n(orthog basis): ";
453 for (int i=0; i<orthog_dim_->blocks()->nblock(); i++) {
454 ExEnv::out0() << scprintf(" %5d", orthog_dim_->blocks()->size(i));
455 }
456 ExEnv::out0() << endl;
457
458 ExEnv::out0() << indent
459 << "WARNING: " << dim_.n() - orthog_dim_.n()
460 << " basis function"
461 << (dim_.n()-orthog_dim_.n()>1?"s":"")
462 << " discarded."
463 << endl;
464 }
465 ExEnv::out0() << indent
466 << "Maximum orthogonalization residual = "
467 << max_orthog_res_ << endl
468 << indent
469 << "Minimum orthogonalization residual = "
470 << min_orthog_res_ << endl;
471
472 if (debug_ > 0) {
473 dim_.print();
474 orthog_dim_.print();
475 if (debug_ > 1) {
476 orthog_trans_.print("basis to orthog basis");
477 orthog_trans_inverse_.print("basis to orthog basis inverse");
478 (orthog_trans_*overlap_
479 *orthog_trans_.t()).print("X*S*X'",ExEnv::out0(),14);
480 (orthog_trans_inverse_.t()*overlap_.gi()
481 *orthog_trans_inverse_).print("X'^(-1)*S^(-1)*X^(-1)",
482 ExEnv::out0(),14);
483 (orthog_trans_
484 *orthog_trans_inverse_).print("X*X^(-1)",ExEnv::out0(),14);
485 }
486 }
487}
488
489// returns the orthogonalization matrix
490RefSCMatrix
491OverlapOrthog::basis_to_orthog_basis()
492{
493 if (orthog_trans_.null()) {
494 compute_orthog_trans();
495 }
496 return orthog_trans_;
497}
498
499RefSCMatrix
500OverlapOrthog::basis_to_orthog_basis_inverse()
501{
502 if (orthog_trans_inverse_.null()) {
503 compute_orthog_trans();
504 }
505 return orthog_trans_inverse_;
506}
507
508RefSCDimension
509OverlapOrthog::dim()
510{
511 return dim_;
512}
513
514RefSCDimension
515OverlapOrthog::orthog_dim()
516{
517 if (orthog_dim_.null()) compute_orthog_trans();
518 return orthog_dim_;
519}
520
521int
522OverlapOrthog::nlindep()
523{
524 if (orthog_dim_.null()) compute_orthog_trans();
525 return nlindep_;
526}
527
528/////////////////////////////////////////////////////////////////////////////
529
530// Local Variables:
531// mode: c++
532// c-file-style: "CLJ-CONDENSED"
533// End:
Note: See TracBrowser for help on using the repository browser.