1 //
2 // gaussbas.cc
3 //
4 // Copyright (C) 1996 Limit Point Systems, Inc.
5 //
6 // Author: Curtis Janssen <cljanss@limitpt.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 <stdio.h>
33
34 #include <scconfig.h>
35 #ifdef HAVE_SSTREAM
36 # include <sstream>
37 #else
38 # include <strstream.h>
39 #endif
40
41 #include <util/keyval/keyval.h>
42 #include <util/misc/formio.h>
43 #include <util/misc/newstring.h>
44 #include <util/state/stateio.h>
45 #include <math/scmat/blocked.h>
46 #include <chemistry/molecule/molecule.h>
47 #include <chemistry/qc/basis/gaussshell.h>
48 #include <chemistry/qc/basis/gaussbas.h>
49 #include <chemistry/qc/basis/files.h>
50 #include <chemistry/qc/basis/cartiter.h>
51 #include <chemistry/qc/basis/transform.h>
52 #include <chemistry/qc/basis/integral.h>
53
54 using namespace std;
55
56 static ClassDesc GaussianBasisSet_cd(
57 typeid(GaussianBasisSet),"GaussianBasisSet",2,"public SavableState",
58 0, create<GaussianBasisSet>, create<GaussianBasisSet>);
59
60 GaussianBasisSet::GaussianBasisSet(const Ref<KeyVal>&topkeyval)
61 1 0 0 0 0 0 0 0 0 0 {
62 molecule_ << topkeyval->describedclassvalue("molecule");
63 if (molecule_.null()) {
64 ExEnv::err() << node0 << indent << "GaussianBasisSet: no \"molecule\"\n";
65 abort();
66 }
67
68 // see if the user requests pure am or cartesian functions
69 int pure;
70 pure = topkeyval->booleanvalue("puream");
71 if (topkeyval->error() != KeyVal::OK) pure = -1;
72
73 // construct a keyval that contains the basis library
74 Ref<KeyVal> keyval;
75
76 if (topkeyval->exists("basisfiles")) {
77 Ref<MessageGrp> grp = MessageGrp::get_default_messagegrp();
78 Ref<ParsedKeyVal> parsedkv = new ParsedKeyVal();
79 char *in_char_array;
80 if (grp->me() == 0) {
81 #ifdef HAVE_SSTREAM
82 ostringstream ostrs;
83 #else
84 ostrstream ostrs;
85 #endif
86 // Look at the basisdir and basisfiles variables to find out what
87 // basis set files are to be read in. The files are read on node
88 // 0 only.
89 ParsedKeyVal::cat_files("basis",topkeyval,ostrs);
90 #ifdef HAVE_SSTREAM
91 int n = 1 + strlen(ostrs.str().c_str());
92 in_char_array = strcpy(new char[n],ostrs.str().c_str());
93 #else
94 ostrs << ends;
95 in_char_array = ostrs.str();
96 int n = ostrs.pcount();
97 #endif
98 grp->bcast(n);
99 grp->bcast(in_char_array, n);
100 }
101 else {
102 int n;
103 grp->bcast(n);
104 in_char_array = new char[n];
105 grp->bcast(in_char_array, n);
106 }
107 parsedkv->parse_string(in_char_array);
108 delete[] in_char_array;
109 Ref<KeyVal> libkeyval = parsedkv.pointer();
110 keyval = new AggregateKeyVal(topkeyval,libkeyval);
111 }
112 else {
113 keyval = topkeyval;
114 }
115
116 // if there isn't a matrixkit in the input, init2() will assign the
117 // default matrixkit
118 matrixkit_ << keyval->describedclassvalue("matrixkit");
119
120 // Bases keeps track of what basis set data bases have already
121 // been read in. It also handles the conversion of basis
122 // names to file names.
123 BasisFileSet bases(keyval);
124 init(molecule_,keyval,bases,1,pure);
125 }
126
127 GaussianBasisSet::GaussianBasisSet(const GaussianBasisSet& gbs) :
128 molecule_(gbs.molecule_),
129 matrixkit_(gbs.matrixkit_),
130 basisdim_(gbs.basisdim_),
131 ncenter_(gbs.ncenter_),
132 nshell_(gbs.nshell_)
133 {
134 int i,j;
135
136 name_ = new_string(gbs.name_);
137
138 center_to_nshell_.set_length(ncenter_);
139 for (i=0; i < ncenter_; i++) {
140 center_to_nshell_(i) = gbs.center_to_nshell_(i);
141 }
142
143 shell_ = new GaussianShell*[nshell_];
144 for (i=0; i<nshell_; i++) {
145 const GaussianShell& gsi = gbs(i);
146
147 int nc=gsi.ncontraction();
148 int np=gsi.nprimitive();
149
150 int *ams = new int[nc];
151 int *pure = new int[nc];
152 double *exps = new double[np];
153 double **coefs = new double*[nc];
154
155 for (j=0; j < nc; j++) {
156 ams[j] = gsi.am(j);
157 pure[j] = gsi.is_pure(j);
158 coefs[j] = new double[np];
159 for (int k=0; k < np; k++)
160 coefs[j][k] = gsi.coefficient_unnorm(j,k);
161 }
162
163 for (j=0; j < np; j++)
164 exps[j] = gsi.exponent(j);
165
166 shell_[i] = new GaussianShell(nc, np, exps, ams, pure, coefs,
167 GaussianShell::Unnormalized);
168 }
169
170 init2();
171 }
172
173 GaussianBasisSet::GaussianBasisSet(StateIn&s):
174 SavableState(s),
175 center_to_nshell_(s)
176 {
177 matrixkit_ = SCMatrixKit::default_matrixkit();
178
179 molecule_ << SavableState::restore_state(s);
180 basisdim_ << SavableState::restore_state(s);
181
182 ncenter_ = center_to_nshell_.size();
183 s.getstring(name_);
184
185 nshell_ = 0;
186 int i;
187 for (i=0; i<ncenter_; i++) {
188 nshell_ += center_to_nshell_(i);
189 }
190
191 shell_ = new GaussianShell*[nshell_];
192 for (i=0; i<nshell_; i++) {
193 shell_[i] = new GaussianShell(s);
194 }
195
196 init2();
197 }
198
199 void
200 GaussianBasisSet::save_data_state(StateOut&s)
201 {
202 center_to_nshell_.save_object_state(s);
203
204 SavableState::save_state(molecule_.pointer(),s);
205 SavableState::save_state(basisdim_.pointer(),s);
206
207 s.putstring(name_);
208 for (int i=0; i<nshell_; i++) {
209 shell_[i]->save_object_state(s);
210 }
211 }
212
213 void
214 GaussianBasisSet::init(Ref<Molecule>&molecule,
215 Ref<KeyVal>&keyval,
216 BasisFileSet& bases,
217 int have_userkeyval,
218 int pur)
219 {
220 int pure, havepure;
221 pure = pur;
222 if (pur == -1) {
223 havepure = 0;
224 }
225 else {
226 havepure = 1;
227 }
228
229 int skip_ghosts = keyval->booleanvalue("skip_ghosts");
230
231 name_ = keyval->pcharvalue("name");
232 int have_custom, nelement;
233
234 if (keyval->exists("basis")) {
235 have_custom = 1;
236 nelement = keyval->count("element");
237 }
238 else {
239 have_custom = 0;
240 nelement = 0;
241 if (!name_) {
242 ExEnv::err() << node0 << indent
243 << "GaussianBasisSet: No name given for basis set\n";
244 abort();
245 }
246 }
247
248 // construct prefixes for each atom: :basis:<atom>:<basisname>:<shell#>
249 // and read in the shell
250 nbasis_ = 0;
251 int ishell = 0;
252 ncenter_ = molecule->natom();
253 int iatom;
254 for (iatom=0; iatom < ncenter_; iatom++) {
255 if (skip_ghosts && molecule->charge(iatom) == 0.0) continue;
256 int Z = molecule->Z(iatom);
257 // see if there is a specific basisname for this atom
258 char* sbasisname = 0;
259 if (have_custom && !nelement) {
260 sbasisname = keyval->pcharvalue("basis",iatom);
261 }
262 else if (nelement) {
263 int i;
264 for (i=0; i<nelement; i++) {
265 char *tmpelementname = keyval->pcharvalue("element", i);
266 int tmpZ = AtomInfo::string_to_Z(tmpelementname);
267 if (tmpZ == Z) {
268 sbasisname = keyval->pcharvalue("basis", i);
269 break;
270 }
271 delete[] tmpelementname;
272 }
273 }
274 if (!sbasisname) {
275 if (!name_) {
276 ExEnv::err() << node0 << indent << "GaussianBasisSet: "
277 << scprintf("no basis name for atom %d (%s)\n",
278 iatom, AtomInfo::name(Z));
279 abort();
280 }
281 sbasisname = strcpy(new char[strlen(name_)+1],name_);
282 }
283 recursively_get_shell(ishell,keyval,
284 AtomInfo::name(Z),
285 sbasisname,bases,havepure,pure,0);
286 delete[] sbasisname;
287 }
288 nshell_ = ishell;
289 shell_ = new GaussianShell*[nshell_];
290 ishell = 0;
291 center_to_nshell_.set_length(ncenter_);
292 for (iatom=0; iatom<ncenter_; iatom++) {
293 if (skip_ghosts && molecule->charge(iatom) == 0.0) {
294 center_to_nshell_[iatom] = 0;
295 continue;
296 }
297 int Z = molecule->Z(iatom);
298 // see if there is a specific basisname for this atom
299 char* sbasisname = 0;
300 if (have_custom && !nelement) {
301 sbasisname = keyval->pcharvalue("basis",iatom);
302 }
303 else if (nelement) {
304 int i;
305 for (i=0; i<nelement; i++) {
306 char *tmpelementname = keyval->pcharvalue("element", i);
307 int tmpZ = AtomInfo::string_to_Z(tmpelementname);
308 if (tmpZ == Z) {
309 sbasisname = keyval->pcharvalue("basis", i);
310 break;
311 }
312 delete[] tmpelementname;
313 }
314 }
315 if (!sbasisname) {
316 if (!name_) {
317 ExEnv::err() << node0 << indent << "GaussianBasisSet: "
318 << scprintf("no basis name for atom %d (%s)\n",
319 iatom, AtomInfo::name(Z));
320 abort();
321 }
322 sbasisname = strcpy(new char[strlen(name_)+1],name_);
323 }
324
325 int ishell_old = ishell;
326 recursively_get_shell(ishell,keyval,
327 AtomInfo::name(Z),
328 sbasisname,bases,havepure,pure,1);
329
330 center_to_nshell_[iatom] = ishell - ishell_old;
331
332 delete[] sbasisname;
333 }
334
335 // delete the name_ if the basis set is customized
336 if (have_custom) {
337 delete[] name_;
338 name_ = 0;
339 }
340
341 // finish with the initialization steps that don't require any
342 // external information
343 init2(skip_ghosts);
344 }
345
346 double
347 GaussianBasisSet::r(int icenter, int xyz) const
348 2 0 1 1 1 0 0 0 0 0 {
349 return molecule_->r(icenter,xyz);
350 2 0 1 0 0 0 0 0 0 0 }
351
352 void
353 GaussianBasisSet::init2(int skip_ghosts)
354 1 0 0 0 0 0 0 0 0 0 {
355 // center_to_shell_ and shell_to_center_
356 shell_to_center_.set_length(nshell_);
357 center_to_shell_.set_length(ncenter_);
358 center_to_nbasis_.set_length(ncenter_);
359 int ishell = 0;
360 for (int icenter=0; icenter<ncenter_; icenter++) {
361 if (skip_ghosts && molecule()->charge(icenter) == 0.0) {
362 center_to_shell_[icenter] = -1;
363 center_to_nbasis_[icenter] = 0;
364 continue;
365 }
366 int j;
367 center_to_shell_[icenter] = ishell;
368 center_to_nbasis_[icenter] = 0;
369 for (j = 0; j<center_to_nshell_[icenter]; j++) {
370 center_to_nbasis_[icenter] += shell_[ishell+j]->nfunction();
371 }
372 ishell += center_to_nshell_[icenter];
373 for (j = center_to_shell_[icenter]; j<ishell; j++) {
374 shell_to_center_[j] = icenter;
375 }
376 }
377
378 // compute nbasis_ and shell_to_function_[]
379 shell_to_function_.set_length(nshell_);
380 nbasis_ = 0;
381 nprim_ = 0;
382 for (ishell=0; ishell<nshell_; ishell++) {
383 shell_to_function_[ishell] = nbasis_;
384 nbasis_ += shell_[ishell]->nfunction();
385 nprim_ += shell_[ishell]->nprimitive();
386 }
387
388 // would like to do this in function_to_shell(), but it is const
389 int n = nbasis();
390 int nsh = nshell();
391 function_to_shell_.set_length(n);
392 int ifunc = 0;
393 for (int i=0; i<nsh; i++) {
394 int nfun = operator()(i).nfunction();
395 for (int j=0; j<nfun; j++) {
396 function_to_shell_[ifunc] = i;
397 ifunc++;
398 }
399 }
400
401 if (matrixkit_.null())
402 matrixkit_ = SCMatrixKit::default_matrixkit();
403
404 so_matrixkit_ = new BlockedSCMatrixKit(matrixkit_);
405
406 if (basisdim_.null()) {
407 int nb = nshell();
408 int *bs = new int[nb];
409 for (int s=0; s < nb; s++)
410 bs[s] = shell(s).nfunction();
411 basisdim_ = new SCDimension(nbasis(), nb, bs, "basis set dimension");
412 delete[] bs;
413 }
414 }
415
416 void
417 GaussianBasisSet::set_matrixkit(const Ref<SCMatrixKit>& mk)
418 {
419 matrixkit_ = mk;
420 so_matrixkit_ = new BlockedSCMatrixKit(matrixkit_);
421 }
422
423 void
424 GaussianBasisSet::
425 recursively_get_shell(int&ishell,Ref<KeyVal>&keyval,
426 const char*element,
427 const char*basisname,
428 BasisFileSet &bases,
429 int havepure,int pure,
430 int get)
431 {
432 char keyword[KeyVal::MaxKeywordLength],prefix[KeyVal::MaxKeywordLength];
433
434 sprintf(keyword,":basis:%s:%s",
435 element,basisname);
436 int count = keyval->count(keyword);
437 if (keyval->error() != KeyVal::OK) {
438 keyval = bases.keyval(keyval, basisname);
439 }
440 count = keyval->count(keyword);
441 if (keyval->error() != KeyVal::OK) {
442 ExEnv::err() << node0 << indent
443 << scprintf("GaussianBasisSet:: couldn't find \"%s\":\n", keyword);
444 keyval->errortrace(ExEnv::err());
445 exit(1);
446 }
447 if (!count) return;
448 1 0 0 0 0 0 0 0 0 0 for (int j=0; j<count; j++) {
449 sprintf(prefix,":basis:%s:%s",
450 element,basisname);
451 Ref<KeyVal> prefixkeyval = new PrefixKeyVal(prefix,keyval,j);
452 if (prefixkeyval->exists("get")) {
453 char* newbasis = prefixkeyval->pcharvalue("get");
454 if (!newbasis) {
455 ExEnv::err() << node0 << indent << "GaussianBasisSet: "
456 << scprintf("error processing get for \"%s\"\n", prefix);
457 keyval->errortrace(ExEnv::err());
458 exit(1);
459 }
460 recursively_get_shell(ishell,keyval,element,newbasis,bases,
461 havepure,pure,get);
462 delete[] newbasis;
463 }
464 else {
465 2 0 0 0 0 0 0 0 0 0 if (get) {
466 if (havepure) shell_[ishell] = new GaussianShell(prefixkeyval,pure);
467 else shell_[ishell] = new GaussianShell(prefixkeyval);
468 }
469 ishell++;
470 }
471 }
472 }
473
474 GaussianBasisSet::~GaussianBasisSet()
475 {
476 delete[] name_;
477
478 int ii;
479 for (ii=0; ii<nshell_; ii++) {
480 delete shell_[ii];
481 }
482 delete[] shell_;
483 }
484
485 int
486 GaussianBasisSet::max_nfunction_in_shell() const
487 {
488 int i;
489 int max = 0;
490 for (i=0; i<nshell_; i++) {
491 if (max < shell_[i]->nfunction()) max = shell_[i]->nfunction();
492 }
493 return max;
494 }
495
496 int
497 GaussianBasisSet::max_ncontraction() const
498 {
499 int i;
500 int max = 0;
501 for (i=0; i<nshell_; i++) {
502 if (max < shell_[i]->ncontraction()) max = shell_[i]->ncontraction();
503 }
504 return max;
505 }
506
507 int
508 GaussianBasisSet::max_angular_momentum() const
509 {
510 int i;
511 int max = 0;
512 for (i=0; i<nshell_; i++) {
513 int maxshi = shell_[i]->max_angular_momentum();
514 if (max < maxshi) max = maxshi;
515 }
516 return max;
517 }
518
519 int
520 GaussianBasisSet::max_cartesian() const
521 {
522 int i;
523 int max = 0;
524 for (i=0; i<nshell_; i++) {
525 int maxshi = shell_[i]->max_cartesian();
526 if (max < maxshi) max = maxshi;
527 }
528 return max;
529 }
530
531 int
532 GaussianBasisSet::max_ncartesian_in_shell(int aminc) const
533 {
534 int i;
535 int max = 0;
536 for (i=0; i<nshell_; i++) {
537 int maxshi = shell_[i]->ncartesian_with_aminc(aminc);
538 if (max < maxshi) max = maxshi;
539 }
540 return max;
541 }
542
543 int
544 GaussianBasisSet::max_am_for_contraction(int con) const
545 {
546 int i;
547 int max = 0;
548 for (i=0; i<nshell_; i++) {
549 if (shell_[i]->ncontraction() <= con) continue;
550 int maxshi = shell_[i]->am(con);
551 if (max < maxshi) max = maxshi;
552 }
553 return max;
554 }
555
556 int
557 GaussianBasisSet::function_to_shell(int func) const
558 2 0 3 1 2 0 0 0 0 0 {
559 return function_to_shell_[func];
560 2 0 2 1 0 0 0 0 0 0 }
561
562 int
563 GaussianBasisSet::ncenter() const
564 1 0 1 0 2 0 0 0 0 0 {
565 1 0 0 0 0 0 0 0 0 0 return ncenter_;
566 }
567
568 int
569 GaussianBasisSet::nshell_on_center(int icenter) const
570 5 0 3 0 2 0 0 0 0 0 {
571 return center_to_nshell_[icenter];
572 1 0 0 0 0 0 0 0 0 0 }
573
574 int
575 GaussianBasisSet::nbasis_on_center(int icenter) const
576 {
577 return center_to_nbasis_[icenter];
578 }
579
580 int
581 GaussianBasisSet::shell_on_center(int icenter, int ishell) const
582 12 0 4 0 3 1 0 0 0 0 {
583 return center_to_shell_(icenter) + ishell;
584 5 0 1 0 3 0 0 0 0 0 }
585
586 const GaussianShell&
587 GaussianBasisSet::operator()(int icenter,int ishell) const
588 {
589 return *shell_[center_to_shell_(icenter) + ishell];
590 }
591
592 GaussianShell&
593 GaussianBasisSet::operator()(int icenter,int ishell)
594 29 0 17 0 6 4 0 0 1 0 {
595 return *shell_[center_to_shell_(icenter) + ishell];
596 15 0 5 0 1 1 0 0 0 0 }
597
598 int
599 GaussianBasisSet::equiv(const Ref<GaussianBasisSet> &b)
600 {
601 if (nshell() != b->nshell()) return 0;
602 for (int i=0; i<nshell(); i++) {
603 if (!shell_[i]->equiv(b->shell_[i])) return 0;
604 }
605 return 1;
606 }
607
608 void
609 GaussianBasisSet::print_brief(ostream& os) const
610 {
611 os << node0 << indent
612 << "GaussianBasisSet:" << endl << incindent
613 << indent << "nbasis = " << nbasis_ << endl
614 << indent << "nshell = " << nshell_ << endl
615 << indent << "nprim = " << nprim_ << endl;
616 if (name_) {
617 os << node0 << indent
618 << "name = \"" << name_ << "\"" << endl;
619 }
620 os << decindent;
621 }
622
623 void
624 GaussianBasisSet::print(ostream& os) const
625 {
626 print_brief(os);
627 if (!SCFormIO::getverbose(os)) return;
628
629 os << incindent;
630
631 // Loop over centers
632 int icenter = 0;
633 int ioshell = 0;
634 for (icenter=0; icenter < ncenter_; icenter++) {
635 os << node0 << endl << indent
636 << scprintf(
637 "center %d: %12.8f %12.8f %12.8f, nshell = %d, shellnum = %d\n",
638 icenter,
639 r(icenter,0),
640 r(icenter,1),
641 r(icenter,2),
642 center_to_nshell_[icenter],
643 center_to_shell_[icenter]);
644 for (int ishell=0; ishell < center_to_nshell_[icenter]; ishell++) {
645 os << node0 << indent
646 << scprintf("Shell %d: functionnum = %d\n",
647 ishell,shell_to_function_[ioshell]);
648 os << node0 << incindent;
649 operator()(icenter,ishell).print(os);
650 os << node0 << decindent;
651 ioshell++;
652 }
653 }
654
655 os << node0 << decindent;
656 }
657
658 /////////////////////////////////////////////////////////////////////////////
659 // GaussianBasisSet::ValueData
660
661 GaussianBasisSet::ValueData::ValueData(
662 const Ref<GaussianBasisSet> &basis,
663 const Ref<Integral> &integral)
664 {
665 maxam_ = basis->max_angular_momentum();
666
667 civec_ = new CartesianIter *[maxam_+1];
668 sivec_ = new SphericalTransformIter *[maxam_+1];
669 for (int i=0; i<=maxam_; i++) {
670 civec_[i] = integral->new_cartesian_iter(i);
671 if (i>1) sivec_[i] = integral->new_spherical_transform_iter(i);
672 else sivec_[i] = 0;
673 }
674 }
675
676 GaussianBasisSet::ValueData::~ValueData()
677 {
678 for (int i=0; i<=maxam_; i++) {
679 delete civec_[i];
680 delete sivec_[i];
681 }
682 delete[] civec_;
683 delete[] sivec_;
684 }
685
686 /////////////////////////////////////////////////////////////////////////////
687
688 // Local Variables:
689 // mode: c++
690 // c-file-style: "CLJ"
691 // End:
692