Annotation of /home/cljanss/src/SC/src/lib/chemistry/qc/scf/scf.cc for ./mpqc.vmon.0018
1 //
2 // scf.cc --- implementation of the SCF abstract base class
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 #ifdef __GNUC__
29 #pragma implementation
30 #endif
31
32 #include <math.h>
33 #include <limits.h>
34 #include <sys/stat.h>
35
36 #include <util/misc/formio.h>
37 #include <util/state/stateio.h>
38 #include <util/group/mstate.h>
39
40 #include <math/scmat/local.h>
41 #include <math/scmat/repl.h>
42 #include <math/scmat/offset.h>
43 #include <math/scmat/blocked.h>
44
45 #include <math/optimize/diis.h>
46
47 #include <chemistry/qc/basis/petite.h>
48 #include <chemistry/qc/scf/scf.h>
49
50 using namespace std;
51
52 ///////////////////////////////////////////////////////////////////////////
53 // SCF
54
55 static ClassDesc SCF_cd(
56 typeid(SCF),"SCF",2,"public OneBodyWavefunction",
57 0, 0, 0);
58
59 SCF::SCF(StateIn& s) :
60 SavableState(s),
61 OneBodyWavefunction(s)
62 {
63 need_vec_ = 1;
64 compute_guess_ = 0;
65
66 s.get(maxiter_,"maxiter");
67 s.get(dens_reset_freq_);
68 s.get(reset_occ_);
69 s.get(local_dens_);
70 s.get(storage_);
71 if (s.version(::class_desc<SCF>()) >= 2) {
72 s.get(print_all_evals_);
73 s.get(print_occ_evals_);
74 }
75 s.get(level_shift_);
76
77 extrap_ << SavableState::restore_state(s);
78 accumdih_ << SavableState::restore_state(s);
79 accumddh_ << SavableState::restore_state(s);
80
81 scf_grp_ = basis()->matrixkit()->messagegrp();
82 threadgrp_ = ThreadGrp::get_default_threadgrp();
83 }
84
85 SCF::SCF(const Ref<KeyVal>& keyval) :
86 OneBodyWavefunction(keyval),
87 need_vec_(1),
88 compute_guess_(0),
89 maxiter_(40),
90 dens_reset_freq_(10),
91 reset_occ_(0),
92 local_dens_(1),
93 storage_(0),
94 level_shift_(0)
95 {
96 if (keyval->exists("maxiter"))
97 maxiter_ = keyval->intvalue("maxiter");
98
99 if (keyval->exists("density_reset_frequency"))
100 dens_reset_freq_ = keyval->intvalue("density_reset_frequency");
101
102 if (keyval->exists("reset_occupations"))
103 reset_occ_ = keyval->booleanvalue("reset_occupations");
104
105 if (keyval->exists("level_shift"))
106 level_shift_ = keyval->doublevalue("level_shift");
107
108 extrap_ << keyval->describedclassvalue("extrap");
109 if (extrap_.null())
110 extrap_ = new DIIS;
111
112 accumdih_ << keyval->describedclassvalue("accumdih");
113 if (accumdih_.null())
114 accumdih_ = new AccumHNull;
115
116 accumddh_ << keyval->describedclassvalue("accumddh");
117 if (accumddh_.null())
118 accumddh_ = new AccumHNull;
119
120 KeyValValueint defaultmem(8000000);
121 storage_ = keyval->intvalue("memory",defaultmem);
122
123 if (keyval->exists("local_density"))
124 local_dens_ = keyval->booleanvalue("local_density");
125
126 print_all_evals_ = keyval->booleanvalue("print_evals");
127 print_occ_evals_ = keyval->booleanvalue("print_occupied_evals");
128
129 scf_grp_ = basis()->matrixkit()->messagegrp();
130 threadgrp_ = ThreadGrp::get_default_threadgrp();
131
132 // first see if guess_wavefunction is a wavefunction, then check to
133 // see if it's a string.
134 if (keyval->exists("guess_wavefunction")) {
135 ExEnv::out() << incindent << incindent;
136 guess_wfn_ << keyval->describedclassvalue("guess_wavefunction");
137 compute_guess_=1;
138 if (guess_wfn_.null()) {
139 compute_guess_=0;
140 char *path = keyval->pcharvalue("guess_wavefunction");
141 struct stat sb;
142 if (path && stat(path, &sb)==0 && sb.st_size) {
143 BcastStateInBin s(scf_grp_, path);
144
145 // reset the default matrixkit so that the matrices in the guess
146 // wavefunction will match those in this wavefunction
147 Ref<SCMatrixKit> oldkit = SCMatrixKit::default_matrixkit();
148 SCMatrixKit::set_default_matrixkit(basis()->matrixkit());
149
150 guess_wfn_ << SavableState::restore_state(s);
151
152 // go back to the original default matrixkit
153 SCMatrixKit::set_default_matrixkit(oldkit);
154 delete[] path;
155 }
156 }
157 ExEnv::out() << decindent << decindent;
158 }
159 }
160
161 SCF::~SCF()
162 {
163 }
164
165 void
166 SCF::save_data_state(StateOut& s)
167 {
168 OneBodyWavefunction::save_data_state(s);
169 s.put(maxiter_);
170 s.put(dens_reset_freq_);
171 s.put(reset_occ_);
172 s.put(local_dens_);
173 s.put(storage_);
174 s.put(print_all_evals_);
175 s.put(print_occ_evals_);
176 s.put(level_shift_);
177 SavableState::save_state(extrap_.pointer(),s);
178 SavableState::save_state(accumdih_.pointer(),s);
179 SavableState::save_state(accumddh_.pointer(),s);
180 }
181
182 RefSCMatrix
183 SCF::oso_eigenvectors()
184 {
185 return oso_eigenvectors_.result();
186 }
187
188 RefDiagSCMatrix
189 SCF::eigenvalues()
190 {
191 return eigenvalues_.result();
192 }
193
194 int
195 SCF::spin_unrestricted()
196 {
197 return 0;
198 }
199
200 void
201 SCF::print(ostream&o) const
202 {
203 OneBodyWavefunction::print(o);
204 o << node0 << indent << "SCF Parameters:\n" << incindent
205 << indent << "maxiter = " << maxiter_ << endl
206 << indent << "density_reset_frequency = " << dens_reset_freq_ << endl
207 << indent << scprintf("level_shift = %f\n",level_shift_)
208 << decindent << endl;
209 }
210
211 //////////////////////////////////////////////////////////////////////////////
212
213 void
214 SCF::compute()
215 {
216 local_ = (dynamic_cast<LocalSCMatrixKit*>(basis()->matrixkit().pointer()) ||
217 dynamic_cast<ReplSCMatrixKit*>(basis()->matrixkit().pointer())) ? 1:0;
218
219 const double hess_to_grad_acc = 1.0/100.0;
220 if (hessian_needed())
221 set_desired_gradient_accuracy(desired_hessian_accuracy()*hess_to_grad_acc);
222
223 const double grad_to_val_acc = 1.0/100.0;
224 if (gradient_needed())
225 set_desired_value_accuracy(desired_gradient_accuracy()*grad_to_val_acc);
226
227 double delta;
228 if (value_needed()) {
229 ExEnv::out() << node0 << endl << indent
230 << scprintf("SCF::compute: energy accuracy = %10.7e\n",
231 desired_value_accuracy())
232 << endl;
233
234 double eelec;
235 delta = compute_vector(eelec);
236
237 // this will be done elsewhere eventually
238 double nucrep = molecule()->nuclear_repulsion_energy();
239 double eother = 0.0;
240 if (accumddh_.nonnull()) eother = accumddh_->e();
241 ExEnv::out() << node0 << endl << indent
242 << scprintf("total scf energy = %15.10f", eelec+eother+nucrep)
243 << endl;
244
245 set_energy(eelec+eother+nucrep);
246 set_actual_value_accuracy(delta);
247 }
248 else {
249 delta = actual_value_accuracy();
250 }
251
252 if (gradient_needed()) {
253 RefSCVector gradient = matrixkit()->vector(moldim());
254
255 ExEnv::out() << node0 << endl << indent
256 << scprintf("SCF::compute: gradient accuracy = %10.7e\n",
257 desired_gradient_accuracy())
258 << endl;
259
260 compute_gradient(gradient);
261 print_natom_3(gradient,"Total Gradient:");
262 set_gradient(gradient);
263
264 set_actual_gradient_accuracy(delta/grad_to_val_acc);
265 }
266
267 if (hessian_needed()) {
268 RefSymmSCMatrix hessian = matrixkit()->symmmatrix(moldim());
269
270 ExEnv::out() << node0 << endl << indent
271 << scprintf("SCF::compute: hessian accuracy = %10.7e\n",
272 desired_hessian_accuracy())
273 << endl;
274
275 compute_hessian(hessian);
276 set_hessian(hessian);
277
278 set_actual_hessian_accuracy(delta/grad_to_val_acc/hess_to_grad_acc);
279 }
280 }
281
282 //////////////////////////////////////////////////////////////////////////////
283
284 signed char *
285 SCF::init_pmax(double *pmat_data)
286 {
287 double l2inv = 1.0/log(2.0);
288 double tol = pow(2.0,-126.0);
289
290 GaussianBasisSet& gbs = *basis().pointer();
291
292 signed char * pmax = new signed char[i_offset(gbs.nshell())];
293
294 int ish, jsh, ij;
295 for (ish=ij=0; ish < gbs.nshell(); ish++) {
296 int istart = gbs.shell_to_function(ish);
297 int iend = istart + gbs(ish).nfunction();
298
299 1 0 0 0 1 0 0 0 0 0 for (jsh=0; jsh <= ish; jsh++,ij++) {
300 int jstart = gbs.shell_to_function(jsh);
301 int jend = jstart + gbs(jsh).nfunction();
302
303 double maxp=0, tmp;
304
305 for (int i=istart; i < iend; i++) {
306 3 0 0 0 0 0 0 0 0 0 int ijoff = i_offset(i) + jstart;
307 2 0 0 0 1 0 0 0 0 0 for (int j=jstart; j < ((ish==jsh) ? i+1 : jend); j++,ijoff++)
308 if ((tmp=fabs(pmat_data[ijoff])) > maxp)
309 2 0 0 0 3 1 0 0 0 0 maxp=tmp;
310 17 0 9 2 9 3 0 0 0 0 }
311
312 2 0 2 0 1 0 0 0 0 0 if (maxp <= tol)
313 1 0 0 0 0 0 0 0 0 0 maxp=tol;
314
315 long power = long(ceil(log(maxp)*l2inv));
316 1 0 0 0 0 0 0 0 0 0 if (power < SCHAR_MIN) pmax[ij] = SCHAR_MIN;
317 else if (power > SCHAR_MAX) pmax[ij] = SCHAR_MAX;
318 else pmax[ij] = (signed char) power;
319 1 0 0 0 1 0 0 0 0 0 }
320 }
321
322 return pmax;
323 }
324
325 //////////////////////////////////////////////////////////////////////////////
326
327 RefSymmSCMatrix
328 SCF::get_local_data(const RefSymmSCMatrix& m, double*& p, Access access)
329 1 0 1 0 0 0 0 0 0 0 {
330 RefSymmSCMatrix l = m;
331
332 if (!dynamic_cast<LocalSymmSCMatrix*>(l.pointer())
333 && !dynamic_cast<ReplSymmSCMatrix*>(l.pointer())) {
334 Ref<SCMatrixKit> k = new ReplSCMatrixKit;
335 l = k->symmmatrix(m.dim());
336 l->convert(m);
337
338 if (access == Accum)
339 l->assign(0.0);
340 } else if (scf_grp_->n() > 1 && access==Accum) {
341 l = m.clone();
342 l.assign(0.0);
343 }
344
345 if (dynamic_cast<ReplSymmSCMatrix*>(l.pointer()))
346 p = dynamic_cast<ReplSymmSCMatrix*>(l.pointer())->get_data();
347 else
348 p = dynamic_cast<LocalSymmSCMatrix*>(l.pointer())->get_data();
349
350 return l;
351 }
352
353 //////////////////////////////////////////////////////////////////////////////
354
355 void
356 SCF::initial_vector(int needv)
357 {
358 1 0 0 0 0 0 0 0 0 0 if (need_vec_) {
359 if (oso_eigenvectors_.result_noupdate().null()) {
360 // if guess_wfn_ is non-null then try to get a guess vector from it.
361 // First check that the same basis is used...if not, then project the
362 // guess vector into the present basis.
363 if (guess_wfn_.nonnull()) {
364 if (basis()->equiv(guess_wfn_->basis())
365 &&orthog_method() == guess_wfn_->orthog_method()) {
366 ExEnv::out() << node0 << indent
367 << "Using guess wavefunction as starting vector" << endl;
368
369 // indent output of eigenvectors() call if there is any
370 ExEnv::out() << incindent << incindent;
371 SCF *g = dynamic_cast<SCF*>(guess_wfn_.pointer());
372 if (!g || compute_guess_) {
373 oso_eigenvectors_ = guess_wfn_->oso_eigenvectors().copy();
374 eigenvalues_ = guess_wfn_->eigenvalues().copy();
375 } else {
376 oso_eigenvectors_ = g->oso_eigenvectors_.result_noupdate().copy();
377 eigenvalues_ = g->eigenvalues_.result_noupdate().copy();
378 }
379 ExEnv::out() << decindent << decindent;
380 } else {
381 ExEnv::out() << node0 << indent
382 << "Projecting guess wavefunction into the present basis set"
383 << endl;
384
385 // indent output of projected_eigenvectors() call if there is any
386 ExEnv::out() << incindent << incindent;
387 oso_eigenvectors_ = projected_eigenvectors(guess_wfn_);
388 eigenvalues_ = projected_eigenvalues(guess_wfn_);
389 ExEnv::out() << decindent << decindent;
390 }
391
392 // we should only have to do this once, so free up memory used
393 // for the old wavefunction
394 guess_wfn_=0;
395
396 ExEnv::out() << node0 << endl;
397
398 } else {
399 ExEnv::out() << node0 << indent << "Starting from core Hamiltonian guess\n"
400 << endl;
401 oso_eigenvectors_ = hcore_guess(eigenvalues_.result_noupdate());
402 }
403 } else {
404 // this is just an old vector
405 }
406 }
407
408 need_vec_=needv;
409 }
410
411 //////////////////////////////////////////////////////////////////////////////
412
413 void
414 SCF::init_mem(int nm)
415 {
416 // if local_den_ is already 0, then that means it was set to zero by
417 // the user.
418 if (!local_dens_) {
419 integral()->set_storage(storage_);
420 return;
421 }
422
423 int nmem = i_offset(basis()->nbasis())*nm*sizeof(double);
424
425 // if we're actually using local matrices, then there's no choice
426 if (dynamic_cast<LocalSCMatrixKit*>(basis()->matrixkit().pointer())
427 ||dynamic_cast<ReplSCMatrixKit*>(basis()->matrixkit().pointer())) {
428 if (nmem > storage_)
429 return;
430 } else {
431 if (nmem > storage_) {
432 local_dens_=0;
433 integral()->set_storage(storage_);
434 return;
435 }
436 }
437
438 integral()->set_storage(storage_-nmem);
439 }
440
441 /////////////////////////////////////////////////////////////////////////////
442
443 void
444 SCF::so_density(const RefSymmSCMatrix& d, double occ, int alp)
445 {
446 int i,j,k;
447 int me=scf_grp_->me();
448 int nproc=scf_grp_->n();
449 int uhf = spin_unrestricted();
450
451 d->assign(0.0);
452
453 RefSCMatrix oso_vector;
454 if (alp || !uhf) {
455 if (oso_scf_vector_.nonnull())
456 oso_vector = oso_scf_vector_;
457 }
458 else {
459 if (oso_scf_vector_beta_.nonnull())
460 oso_vector = oso_scf_vector_beta_;
461 }
462
463 if (oso_vector.null()) {
464 if (uhf) {
465 if (alp)
466 oso_vector = oso_alpha_eigenvectors();
467 else
468 oso_vector = oso_beta_eigenvectors();
469 } else
470 oso_vector = oso_eigenvectors();
471 }
472
473 1 0 0 0 0 0 0 0 0 0 if (debug_ > 1) oso_vector.print("ortho SO vector");
474
475 RefSCMatrix vector = so_to_orthog_so().t() * oso_vector;
476 oso_vector = 0;
477
478 if (debug_ > 1) vector.print("SO vector");
479
480 BlockedSCMatrix *bvec = require_dynamic_cast<BlockedSCMatrix*>(
481 vector, "SCF::so_density: blocked vector");
482
483 BlockedSymmSCMatrix *bd = require_dynamic_cast<BlockedSymmSCMatrix*>(
484 d, "SCF::so_density: blocked density");
485
486 for (int ir=0; ir < oso_dimension()->blocks()->nblock(); ir++) {
487 RefSCMatrix vir = bvec->block(ir);
488 RefSymmSCMatrix dir = bd->block(ir);
489
490 if (vir.null() || vir.ncol()==0)
491 continue;
492
493 int n_orthoSO = oso_dimension()->blocks()->size(ir);
494 int n_SO = so_dimension()->blocks()->size(ir);
495
496 // figure out which columns of the scf vector we'll need
497 int col0 = -1, coln = -1;
498 for (i=0; i < n_orthoSO; i++) {
499 double occi;
500 1 0 0 0 0 0 0 0 0 0 if (!uhf)
501 0 0 1 0 0 0 0 0 0 0 occi = occupation(ir, i);
502 else if (alp)
503 occi = alpha_occupation(ir, i);
504 else
505 occi = beta_occupation(ir, i);
506
507 if (fabs(occi-occ) < 1.0e-8) {
508 if (col0 == -1)
509 col0 = i;
510 continue;
511 0 0 1 0 0 0 0 0 0 0 } else if (col0 != -1) {
512 coln = i-1;
513 break;
514 }
515 }
516
517 if (col0 == -1)
518 continue;
519
520 1 0 0 0 0 0 0 0 0 0 if (coln == -1)
521 coln = n_orthoSO-1;
522
523 1 0 0 0 0 0 0 0 0 0 if (local_ || local_dens_) {
524 RefSymmSCMatrix ldir = dir;
525
526 RefSCMatrix occbits; // holds the occupied bits of the scf vector
527
528 // get local copies of vector and density matrix
529 if (!local_) {
530 Ref<SCMatrixKit> rk = new ReplSCMatrixKit;
531 RefSCMatrix lvir = rk->matrix(vir.rowdim(), vir.coldim());
532 lvir->convert(vir);
533 occbits = lvir->get_subblock(0, n_SO-1, col0, coln);
534 lvir = 0;
535
536 ldir = rk->symmmatrix(dir.dim());
537 ldir->convert(dir);
538
539 } else {
540 occbits = vir->get_subblock(0, n_SO-1, col0, coln);
541 }
542
543 double **c;
544 double *dens;
545
546 if (dynamic_cast<LocalSCMatrix*>(occbits.pointer()))
547 c = dynamic_cast<LocalSCMatrix*>(occbits.pointer())->get_rows();
548 else if (dynamic_cast<ReplSCMatrix*>(occbits.pointer()))
549 c = dynamic_cast<ReplSCMatrix*>(occbits.pointer())->get_rows();
550 else
551 abort();
552
553 if (dynamic_cast<LocalSymmSCMatrix*>(ldir.pointer()))
554 dens = dynamic_cast<LocalSymmSCMatrix*>(ldir.pointer())->get_data();
555 else if (dynamic_cast<ReplSymmSCMatrix*>(ldir.pointer()))
556 dens = dynamic_cast<ReplSymmSCMatrix*>(ldir.pointer())->get_data();
557 else
558 abort();
559
560 int ij=0;
561 for (i=0; i < n_SO; i++) {
562 for (j=0; j <= i; j++, ij++) {
563 17 2 11 8 9 0 0 0 0 0 if (ij%nproc != me)
564 continue;
565
566 double dv = 0;
567
568 0 0 0 1 1 0 0 0 0 0 int kk=0;
569 1 0 3 0 0 0 0 0 1 0 for (k=col0; k <= coln; k++, kk++)
570 11 1 10 6 6 1 0 0 0 0 dv += c[i][kk]*c[j][kk];
571
572 1 0 0 0 0 0 0 0 0 0 dens[ij] = dv;
573 3 0 0 1 0 0 0 0 0 0 }
574 }
575
576 if (nproc > 1)
577 scf_grp_->sum(dens, i_offset(n_SO));
578
579 if (!local_) {
580 dir->convert(ldir);
581 }
582 }
583
584 // for now quit
585 else {
586 ExEnv::err() << node0 << indent
587 << "Cannot yet use anything but Local matrices"
588 << endl;
589 abort();
590 }
591 }
592
593 if (debug_ > 0) {
594 ExEnv::out() << node0 << indent
595 << "Nelectron = " << 2.0 * (d * overlap()).trace() << endl;
596 }
597
598 int use_alternate_density = 0;
599 if (use_alternate_density || debug_ > 2) {
600 // double check the density with this simpler, slower way to compute
601 // the density matrix
602 RefSymmSCMatrix occ(oso_dimension(), basis_matrixkit());
603 occ.assign(0.0);
604 for (i=0; i<oso_dimension()->n(); i++) {
605 occ(i,i) = occupation(i);
606 }
607 occ.scale(0.5);
608 RefSymmSCMatrix d2(so_dimension(), basis_matrixkit());
609 d2.assign(0.0);
610 d2.accumulate_transform(vector, occ);
611 if (debug_ > 2) {
612 d2.print("d2 density");
613 ExEnv::out() << node0 << indent << "d2 Nelectron = "
614 << 2.0 * (d2 * overlap()).trace() << endl;
615 }
616 if (use_alternate_density) {
617 d.assign(d2);
618 ExEnv::out() << node0 << indent
619 << "using alternate density algorithm" << endl;
620 }
621 }
622
623 1 0 0 0 0 0 0 0 0 0 if (debug_ > 1) {
624 d.print("SO Density");
625 RefSCMatrix rd(d.dim(), d.dim(), basis_matrixkit());
626 rd.assign(0.0);
627 rd.accumulate(d);
628 (d*overlap()*d-rd).print("SO Density idempotency error");
629 }
630 0 0 0 1 0 0 0 0 0 0 }
631
632 double
633 SCF::one_body_energy()
634 {
635 RefSymmSCMatrix dens = ao_density().copy();
636 RefSymmSCMatrix hcore = dens->clone();
637 hcore.assign(0.0);
638 Ref<SCElementOp> hcore_op = new OneBodyIntOp(integral()->hcore());
639 hcore.element_op(hcore_op);
640
641 dens->scale_diagonal(0.5);
642 SCElementScalarProduct *prod = new SCElementScalarProduct;
643 prod->reference();
644 Ref<SCElementOp2> op = prod;
645 hcore->element_op(prod, dens);
646 double e = prod->result();
647 op = 0;
648 prod->dereference();
649 delete prod;
650 return 2.0 * e;
651 }
652
653 void
654 SCF::two_body_energy(double &ec, double &ex)
655 {
656 ExEnv::err() << class_name() << ": two_body_energy not implemented" << endl;
657 }
658
659 /////////////////////////////////////////////////////////////////////////////
660
661 void
662 SCF::init_threads()
663 1 0 0 0 0 0 0 0 0 0 {
664 int nthread = threadgrp_->nthread();
665 int int_store = integral()->storage_unused()/nthread;
666
667 // initialize the two electron integral classes
668 tbis_ = new Ref<TwoBodyInt>[nthread];
669 for (int i=0; i < nthread; i++) {
670 tbis_[i] = integral()->electron_repulsion();
671 tbis_[i]->set_integral_storage(int_store);
672 }
673
674 }
675
676 void
677 SCF::done_threads()
678 {
679 for (int i=0; i < threadgrp_->nthread(); i++) tbis_[i] = 0;
680 delete[] tbis_;
681 tbis_ = 0;
682 }
683
684 int *
685 SCF::read_occ(const Ref<KeyVal> &keyval, const char *name, int nirrep)
686 {
687 int *occ = 0;
688 if (keyval->exists(name)) {
689 if (keyval->count(name) != nirrep) {
690 ExEnv::err() << node0 << indent
691 << "ERROR: SCF: have " << nirrep << " irreps but "
692 << name << " vector is length " << keyval->count(name)
693 << endl;
694 abort();
695 }
696 occ = new int[nirrep];
697 for (int i=0; i<nirrep; i++) {
698 occ[i] = keyval->intvalue(name,i);
699 }
700 }
701 return occ;
702 }
703
704 /////////////////////////////////////////////////////////////////////////////
705
706 // Local Variables:
707 // mode: c++
708 // c-file-style: "ETS"
709 // End:
710