1 //
2 // btest.cc --- test program
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 #include <scconfig.h>
29 #ifdef HAVE_SSTREAM
30 # include <sstream>
31 #else
32 # include <strstream.h>
33 #endif
34
35 #include <util/misc/formio.h>
36 #include <util/keyval/keyval.h>
37 #include <util/state/stateio.h>
38 #include <util/state/state_text.h>
39 #include <util/state/state_bin.h>
40 #include <chemistry/qc/basis/basis.h>
41 #include <chemistry/qc/basis/files.h>
42 #include <chemistry/qc/basis/petite.h>
43 #include <chemistry/qc/basis/gpetite.h>
44 #include <chemistry/qc/basis/symmint.h>
45 #include <chemistry/qc/intv3/intv3.h>
46 #include <chemistry/qc/basis/sobasis.h>
47 #include <chemistry/qc/basis/sointegral.h>
48 #include <chemistry/qc/basis/extent.h>
49 #include <chemistry/qc/basis/orthog.h>
50
51 using namespace std;
52 using namespace sc;
53
54 static void
do_so_shell_test(const Ref<SOBasis> & sobas,const Ref<TwoBodySOInt> & soer,int i,int j,int k,int l)55 do_so_shell_test(const Ref<SOBasis>& sobas, const Ref<TwoBodySOInt> &soer,
56 int i, int j, int k, int l)
57 {
58 if (i>=soer->basis1()->nshell()
59 ||j>=soer->basis2()->nshell()
60 ||k>=soer->basis3()->nshell()
61 ||l>=soer->basis4()->nshell()) return;
62
63 int p,q,r,s;
64 soer->compute_shell(i,j,k,l);
65 const double *buf = soer->buffer();
66 int np = sobas->nfunction(i);
67 int nq = sobas->nfunction(j);
68 int nr = sobas->nfunction(k);
69 int ns = sobas->nfunction(l);
70 int off = 0;
71 cout << "SHELL ("<<i<<j<<"|"<<k<<l<<"):" << endl;
72 for (p=0; p<np; p++) {
73 for (q=0; q<nq; q++) {
74 for (r=0; r<nr; r++) {
75 cout << " ("<<p<<q<<"|"<<r<<"*) =";
76 for (s=0; s<ns; s++, off++) {
77 cout << scprintf(" % 10.6f",buf[off]);
78 }
79 cout << endl;
80 }
81 }
82 }
83 }
84
85 static void
do_so_shell_test(const Ref<SOBasis> & sobas,const Ref<OneBodySOInt> & soov,int i,int j)86 do_so_shell_test(const Ref<SOBasis>& sobas, const Ref<OneBodySOInt> &soov,
87 int i, int j)
88 {
89 if (i>=soov->basis1()->nshell()
90 ||j>=soov->basis2()->nshell()) return;
91
92 int p,q;
93 soov->compute_shell(i,j);
94 const double *buf = soov->buffer();
95 int np = sobas->nfunction(i);
96 int nq = sobas->nfunction(j);
97 int off = 0;
98 cout << "SHELL ("<<i<<"|"<<j<<"):" << endl;
99 for (p=0; p<np; p++) {
100 cout << " ("<<p<<"|"<<"*) =";
101 for (q=0; q<nq; q++, off++) {
102 cout << scprintf(" % 10.6f",buf[off]);
103 }
104 cout << endl;
105 }
106 }
107
108 static void
do_so_test(const Ref<KeyVal> & keyval,const Ref<Integral> & intgrl,const Ref<GaussianBasisSet> & gbs)109 do_so_test(const Ref<KeyVal> &keyval,
110 const Ref<Integral>& intgrl, const Ref<GaussianBasisSet> &gbs)
111 {
112 intgrl->set_basis(gbs);
113
114 Ref<SOBasis> sobas = new SOBasis(gbs, intgrl);
115 sobas->print(cout);
116
117 Ref<TwoBodyInt> aoer = intgrl->electron_repulsion();
118 Ref<TwoBodySOInt> soer = new TwoBodySOInt(aoer);
119
120 Ref<OneBodyInt> aoov = intgrl->overlap();
121 Ref<OneBodySOInt> soov = new OneBodySOInt(aoov);
122
123 sobas = soer->basis();
124 sobas->print(cout);
125
126 if (keyval->exists(":shell")) {
127 do_so_shell_test(sobas, soer,
128 keyval->intvalue(":shell",0),
129 keyval->intvalue(":shell",1),
130 keyval->intvalue(":shell",2),
131 keyval->intvalue(":shell",3));
132 do_so_shell_test(sobas, soov,
133 keyval->intvalue(":shell",0),
134 keyval->intvalue(":shell",1));
135 }
136 else {
137 int i,j,k,l;
138 cout << "SO Electron Repulsion:" << endl;
139 for (i=0; i<sobas->nshell(); i++) {
140 for (j=0; j<sobas->nshell(); j++) {
141 for (k=0; k<sobas->nshell(); k++) {
142 for (l=0; l<sobas->nshell(); l++) {
143 do_so_shell_test(sobas, soer, i, j, k, l);
144 }
145 }
146 }
147 }
148 cout << "SO Overlap:" << endl;
149 for (i=0; i<sobas->nshell(); i++) {
150 for (j=0; j<sobas->nshell(); j++) {
151 do_so_shell_test(sobas, soov, i, j);
152 }
153 }
154 }
155 }
156
157 static void
test_overlap(const Ref<GaussianBasisSet> & gbs,const Ref<GaussianBasisSet> & gbs2,const Ref<Integral> & intgrl)158 test_overlap(const Ref<GaussianBasisSet>& gbs, const Ref<GaussianBasisSet>& gbs2,
159 const Ref<Integral>& intgrl)
160 {
161 intgrl->set_basis(gbs);
162
163 // first form AO basis overlap
164 RefSymmSCMatrix s(gbs->basisdim(), gbs->matrixkit());
165 Ref<SCElementOp> ov
166 = new OneBodyIntOp(new OneBodyIntIter(intgrl->overlap()));
167 s.assign(0.0);
168 s.element_op(ov);
169 ov=0;
170 s.print("overlap");
171
172 // now transform s to SO basis
173 Ref<PetiteList> pl = intgrl->petite_list();
174 RefSymmSCMatrix sb = pl->to_SO_basis(s);
175 sb.print("blocked s");
176
177 // and back to AO basis
178 s = pl->to_AO_basis(sb);
179 s.print("reconstituted s");
180
181 // form skeleton overlap
182 ov = new OneBodyIntOp(new SymmOneBodyIntIter(intgrl->overlap(),pl));
183 s.assign(0.0);
184 s.element_op(ov);
185 ov=0;
186 s.print("overlap");
187
188 // and symmetrize to get blocked overlap again
189 sb.assign(0.0);
190 pl->symmetrize(s,sb);
191 sb.print("blocked again");
192
193 s=0; sb=0;
194
195 // now try overlap between two basis sets
196 RefSCMatrix ssq(gbs2->basisdim(),gbs->basisdim(),gbs2->matrixkit());
197 intgrl->set_basis(gbs2,gbs);
198
199 ov = new OneBodyIntOp(new OneBodyIntIter(intgrl->overlap()));
200 ssq.assign(0.0);
201 ssq.element_op(ov);
202 ssq.print("overlap sq");
203 ov=0;
204
205 Ref<PetiteList> pl2 = intgrl->petite_list(gbs2);
206 RefSCMatrix ssqb(pl2->AO_basisdim(), pl->AO_basisdim(), gbs->so_matrixkit());
207 ssqb->convert(ssq);
208
209 RefSCMatrix syms2 = pl2->aotoso().t() * ssqb * pl->aotoso();
210 syms2.print("symm S2");
211 }
212
213 static void
test_aoorthog(const Ref<GaussianBasisSet> & gbs,const Ref<Integral> & intgrl)214 test_aoorthog(const Ref<GaussianBasisSet>& gbs,
215 const Ref<Integral>& intgrl)
216 {
217 cout << "Beginning AO Orthog test" << endl;
218
219 intgrl->set_basis(gbs);
220
221 // first form AO basis overlap
222 RefSymmSCMatrix s(gbs->basisdim(), gbs->matrixkit());
223 Ref<SCElementOp> ov
224 = new OneBodyIntOp(new OneBodyIntIter(intgrl->overlap()));
225 s.assign(0.0);
226 s.element_op(ov);
227 ov=0;
228
229 Ref<OverlapOrthog> orthog;
230 orthog = new OverlapOrthog(OverlapOrthog::Symmetric,
231 s,
232 s.kit(),
233 0.0001,
234 1);
235 orthog->basis_to_orthog_basis().print("basis to orthog basis");
236
237 s = 0;
238 orthog = 0;
239
240 cout << "Ending AO Orthog test" << endl;
241 }
242
243 static void
test_eigvals(const Ref<GaussianBasisSet> & gbs,const Ref<Integral> & intgrl)244 test_eigvals(const Ref<GaussianBasisSet>& gbs, const Ref<Integral>& intgrl)
245 {
246 intgrl->set_basis(gbs);
247 Ref<PetiteList> pl = intgrl->petite_list();
248
249 // form AO Hcore and evecs
250 RefSymmSCMatrix hcore_ao(gbs->basisdim(), gbs->matrixkit());
251 RefSCMatrix ao_evecs(gbs->basisdim(), gbs->basisdim(), gbs->matrixkit());
252 RefDiagSCMatrix ao_evals(gbs->basisdim(), gbs->matrixkit());
253
254 hcore_ao.assign(0.0);
255
256 Ref<SCElementOp> op
257 = new OneBodyIntOp(new OneBodyIntIter(intgrl->kinetic()));
258 hcore_ao.element_op(op);
259 op=0;
260
261 Ref<OneBodyInt> nuc = intgrl->nuclear();
262 nuc->reinitialize();
263 op = new OneBodyIntOp(nuc);
264 hcore_ao.element_op(op);
265 op=0;
266
267 hcore_ao.print("Hcore (AO)");
268
269 hcore_ao.diagonalize(ao_evals, ao_evecs);
270 ao_evecs.print("AO Evecs");
271 ao_evals.print("AO Evals");
272
273 // form SO Hcore and evecs
274 RefSymmSCMatrix hcore_so(pl->SO_basisdim(), gbs->so_matrixkit());
275 RefSCMatrix so_evecs(pl->SO_basisdim(), pl->SO_basisdim(),
276 gbs->so_matrixkit());
277 RefDiagSCMatrix so_evals(pl->SO_basisdim(), gbs->so_matrixkit());
278
279 // reuse hcore_ao to get skeleton Hcore
280 hcore_ao.assign(0.0);
281
282 op = new OneBodyIntOp(new SymmOneBodyIntIter(intgrl->kinetic(),pl));
283 hcore_ao.element_op(op);
284 op=0;
285
286 nuc = intgrl->nuclear();
287 nuc->reinitialize();
288 op = new OneBodyIntOp(new SymmOneBodyIntIter(nuc,pl));
289 hcore_ao.element_op(op);
290 op=0;
291
292 pl->symmetrize(hcore_ao, hcore_so);
293
294 hcore_so.print("Hcore (SO)");
295
296 hcore_so.diagonalize(so_evals, so_evecs);
297 so_evecs.print("SO Evecs");
298 so_evals.print("SO Evals");
299
300 RefSCMatrix new_ao_evecs = pl->evecs_to_AO_basis(so_evecs);
301 new_ao_evecs.print("AO Evecs again");
302
303 //RefSCMatrix new_so_evecs = pl->evecs_to_SO_basis(ao_evecs);
304 //new_so_evecs.print("SO Evecs again");
305
306 pl->to_AO_basis(hcore_so).print("Hcore (AO) again");
307 }
308
309 void
checkerror(const char * name,int shell,int func,double numerical,double check)310 checkerror(const char *name, int shell, int func,
311 double numerical, double check)
312 {
313 double mag = fabs(check);
314 double err = fabs(numerical - check);
315 cout << scprintf("%2s %2d %2d %12.8f %12.8f er = %6.4f",
316 name, shell, func,
317 numerical, check, err/mag) << endl;
318 if (mag > 0.001) {
319 if (err/mag > 0.05) {
320 cout << scprintf("ERROR %2s %2d %2d %12.8f %12.8f er = %6.4f",
321 name, shell, func,
322 numerical, check, err/mag) << endl;
323 }
324 }
325 else if (err > 0.02) {
326 cout << scprintf("ERROR %2s %2d %2d %12.8f %12.8f ea = %16.14f",
327 name, shell, func, numerical, check, err) << endl;
328 }
329 }
330
331 void
test_func_values(const Ref<GaussianBasisSet> & gbs,const Ref<Integral> & integral)332 test_func_values(const Ref<GaussianBasisSet> &gbs,
333 const Ref<Integral> &integral)
334 {
335 cout << "testing basis function value gradient and hessian numerically"
336 << endl;
337
338 int nbasis = gbs->nbasis();
339 double *b_val = new double[nbasis];
340 double *b_val_plsx = new double[nbasis];
341 double *b_val_mnsx = new double[nbasis];
342 double *b_val_plsy = new double[nbasis];
343 double *b_val_mnsy = new double[nbasis];
344 double *b_val_plsz = new double[nbasis];
345 double *b_val_mnsz = new double[nbasis];
346 double *b_val_plsyx = new double[nbasis];
347 double *b_val_mnsyx = new double[nbasis];
348 double *b_val_plszy = new double[nbasis];
349 double *b_val_mnszy = new double[nbasis];
350 double *b_val_plszx = new double[nbasis];
351 double *b_val_mnszx = new double[nbasis];
352 double *g_val = new double[3*nbasis];
353 double *h_val = new double[6*nbasis];
354
355 const int x_ = 0;
356 const int y_ = 1;
357 const int z_ = 2;
358 const int xx_ = 0;
359 const int yx_ = 1;
360 const int yy_ = 2;
361 const int zx_ = 3;
362 const int zy_ = 4;
363 const int zz_ = 5;
364
365 SCVector3 r;
366 SCVector3 d;
367 double delta = 0.001;
368 SCVector3 dx(delta, 0., 0.);
369 SCVector3 dy(0., delta, 0.);
370 SCVector3 dz(0., 0., delta);
371 SCVector3 dxy(delta, delta, 0.);
372 SCVector3 dxz(delta, 0., delta);
373 SCVector3 dyz(0., delta, delta);
374 double deltax = 0.1;
375 GaussianBasisSet::ValueData vdat(gbs, integral);
376 for (r.x()=0.0; r.x() < 1.0; r.x() += deltax) {
377 deltax *= 2.;
378 double deltay = 0.1;
379 for (r.y()=0.0; r.y() < 1.0; r.y() += deltay) {
380 deltay *= 2.;
381 double deltaz = 0.1;
382 for (r.z()=0.0; r.z() < 1.0; r.z() += deltaz) {
383 deltaz *= 2.;
384 cout << "R = " << r << endl;
385 gbs->hessian_values(r, &vdat, h_val, g_val, b_val);
386 gbs->values(r + dx, &vdat, b_val_plsx);
387 gbs->values(r - dx, &vdat, b_val_mnsx);
388 gbs->values(r + dy, &vdat, b_val_plsy);
389 gbs->values(r - dy, &vdat, b_val_mnsy);
390 gbs->values(r + dz, &vdat, b_val_plsz);
391 gbs->values(r - dz, &vdat, b_val_mnsz);
392 gbs->values(r + dxy, &vdat, b_val_plsyx);
393 gbs->values(r - dxy, &vdat, b_val_mnsyx);
394 gbs->values(r + dyz, &vdat, b_val_plszy);
395 gbs->values(r - dyz, &vdat, b_val_mnszy);
396 gbs->values(r + dxz, &vdat, b_val_plszx);
397 gbs->values(r - dxz, &vdat, b_val_mnszx);
398 for (int i=0; i<nbasis; i++) {
399 int shell = gbs->function_to_shell(i);
400 int func = i - gbs->shell_to_function(shell);
401 double g_val_test[3];
402 double h_val_test[6];
403 int x = i*3+x_;
404 int y = i*3+y_;
405 int z = i*3+z_;
406 g_val_test[x_] = 0.5*(b_val_plsx[i] - b_val_mnsx[i])/delta;
407 g_val_test[y_] = 0.5*(b_val_plsy[i] - b_val_mnsy[i])/delta;
408 g_val_test[z_] = 0.5*(b_val_plsz[i] - b_val_mnsz[i])/delta;
409 int xx = i*6+xx_;
410 int yx = i*6+yx_;
411 int yy = i*6+yy_;
412 int zx = i*6+zx_;
413 int zy = i*6+zy_;
414 int zz = i*6+zz_;
415 h_val_test[xx_]
416 = (b_val_plsx[i] + b_val_mnsx[i] - 2. * b_val[i])
417 * 1./(delta*delta);
418 h_val_test[yy_]
419 = (b_val_plsy[i] + b_val_mnsy[i] - 2. * b_val[i])
420 * 1./(delta*delta);
421 h_val_test[zz_]
422 = (b_val_plsz[i] + b_val_mnsz[i] - 2. * b_val[i])
423 * 1./(delta*delta);
424 h_val_test[yx_]
425 = 0.5 * ((b_val_plsyx[i]+b_val_mnsyx[i]-2.*b_val[i])
426 * 1. / (delta * delta)
427 - h_val_test[xx_] - h_val_test[yy_]);
428 h_val_test[zx_]
429 = 0.5 * ((b_val_plszx[i]+b_val_mnszx[i]-2.*b_val[i])
430 * 1. / (delta * delta)
431 - h_val_test[xx_] - h_val_test[zz_]);
432 h_val_test[zy_]
433 = 0.5 * ((b_val_plszy[i]+b_val_mnszy[i]-2.*b_val[i])
434 * 1. / (delta * delta)
435 - h_val_test[zz_] - h_val_test[yy_]);
436 checkerror("x", shell, func, g_val_test[x_], g_val[x]);
437 checkerror("y", shell, func, g_val_test[y_], g_val[y]);
438 checkerror("z", shell, func, g_val_test[z_], g_val[z]);
439 checkerror("xx", shell, func, h_val_test[xx_], h_val[xx]);
440 checkerror("yy", shell, func, h_val_test[yy_], h_val[yy]);
441 checkerror("zz", shell, func, h_val_test[zz_], h_val[zz]);
442 checkerror("yx", shell, func, h_val_test[yx_], h_val[yx]);
443 checkerror("zx", shell, func, h_val_test[zx_], h_val[zx]);
444 checkerror("zy", shell, func, h_val_test[zy_], h_val[zy]);
445 }
446 }
447 }
448 }
449 delete[] b_val;
450 delete[] b_val_plsx;
451 delete[] b_val_mnsx;
452 delete[] b_val_plsy;
453 delete[] b_val_mnsy;
454 delete[] b_val_plsz;
455 delete[] b_val_mnsz;
456 delete[] b_val_plsyx;
457 delete[] b_val_mnsyx;
458 delete[] b_val_plszy;
459 delete[] b_val_mnszy;
460 delete[] b_val_plszx;
461 delete[] b_val_mnszx;
462 delete[] g_val;
463 delete[] h_val;
464 }
465
466 void
do_extent_test(const Ref<GaussianBasisSet> & gbs)467 do_extent_test(const Ref<GaussianBasisSet> &gbs)
468 {
469 int i, j;
470 for (i=0; i<gbs->nshell(); i++) {
471 gbs->shell(i).print();
472 for (j=0; j<10; j++) {
473 cout << " " << gbs->shell(i).monobound(0.1*j);
474 }
475 cout << endl;
476 for (j=0; j<10; j++) {
477 double threshold = pow(10.0, -j);
478 //cout << " threshold = " << threshold << endl;
479 cout << " " << gbs->shell(i).extent(threshold);
480 }
481 cout << endl;
482 }
483
484 Ref<ShellExtent> extent = new ShellExtent;
485 extent->init(gbs);
486 extent->print();
487 }
488
489 void
do_gpetite_test(const Ref<GaussianBasisSet> & b1,const Ref<GaussianBasisSet> & b2)490 do_gpetite_test(const Ref<GaussianBasisSet> &b1,
491 const Ref<GaussianBasisSet> &b2)
492 {
493 canonical_aabc c4(b1,b1,b1,b2);
494 Ref<GPetite4<canonical_aabc> > p4
495 = new GPetite4<canonical_aabc>(b1,b1,b1,b2,c4);
496 cout << "tesing GPetite4<canonical_aabc>" << endl;
497 for (int i=0; i<b1->nshell(); i++) {
498 for (int j=0; j<=i; j++) {
499 for (int k=0; k<b1->nshell(); k++) {
500 for (int l=0; l<b2->nshell(); l++) {
501 cout << " " << i
502 << " " << j
503 << " " << k
504 << " " << l
505 << " in p4: " << p4->in_p4(i,j,k,l)
506 << endl;
507 }
508 }
509 }
510 }
511 }
512
513 int
main(int,char * argv[])514 main(int, char *argv[])
515 {
516 int i, j;
517
518 char o[10000];
519 #ifdef HAVE_SSTREAM
520 ostringstream perlout(o);
521 #else
522 ostrstream perlout(o,sizeof(o));
523 #endif
524
525 const char *filename = (argv[1]) ? argv[1] : SRCDIR "/btest.kv";
526
527 Ref<KeyVal> keyval = new ParsedKeyVal(filename);
528
529 Ref<Integral> intgrl = new IntegralV3;
530
531 int doconcat = keyval->booleanvalue("concat");
532 int dooverlap = keyval->booleanvalue("overlap");
533 int doeigvals = keyval->booleanvalue("eigvals");
534 int dostate = keyval->booleanvalue("state");
535 int doso = keyval->booleanvalue("so");
536 int doatoms = keyval->booleanvalue("atoms");
537 int dopetite = keyval->booleanvalue("petite");
538 int dogpetite = keyval->booleanvalue("gpetite");
539 int dovalues = keyval->booleanvalue("values");
540 int doextent = keyval->booleanvalue("extent");
541 int doaoorthog = keyval->booleanvalue("aoorthog");
542
543 if (doconcat) {
544 Ref<GaussianBasisSet> b1, b2;
545 b1 << keyval->describedclassvalue("concat1");
546 b2 << keyval->describedclassvalue("concat2");
547 Ref<GaussianBasisSet> b12 = b1->operator+(b2);
548 Ref<GaussianBasisSet> b121 = b12->operator+(b1);
549 b1->print();
550 b2->print();
551 b12->print();
552 b121->print();
553 }
554
555 for (i=0; i<keyval->count("test"); i++) {
556 Ref<GaussianBasisSet> gbs;
557 gbs << keyval->describedclassvalue("test", i);
558 Ref<GaussianBasisSet> gbs2;
559 gbs2 << keyval->describedclassvalue("test2", i);
560
561 if (dooverlap) test_overlap(gbs,gbs2,intgrl);
562
563 if (doaoorthog) test_aoorthog(gbs,intgrl);
564
565 if (doeigvals) test_eigvals(gbs,intgrl);
566
567 if (dostate) {
568 StateOutBin out("btest.out");
569 SavableState::save_state(gbs.pointer(),out);
570 out.close();
571 StateInBin in("btest.out");
572 gbs << SavableState::restore_state(in);
573 gbs->print();
574 }
575
576 if (dopetite) {
577 intgrl->set_basis(gbs);
578 intgrl->petite_list()->print();
579 }
580
581 if (dogpetite) {
582 do_gpetite_test(gbs,gbs2);
583 }
584
585 if (doso) {
586 do_so_test(keyval, intgrl, gbs);
587 }
588
589 if (dovalues) {
590 intgrl->set_basis(gbs);
591 test_func_values(gbs,intgrl);
592 }
593
594 if (doextent) {
595 do_extent_test(gbs);
596 }
597 }
598
599 if (doatoms) {
600 const int nelem = 37;
601
602 // Make H, C, and P molecules
603 Ref<Molecule> hmol = new Molecule();
604 hmol->add_atom(hmol->atominfo()->string_to_Z("H"),0,0,0);
605 Ref<Molecule> cmol = new Molecule();
606 cmol->add_atom(cmol->atominfo()->string_to_Z("C"),0,0,0);
607 Ref<Molecule> pmol = new Molecule();
608 pmol->add_atom(pmol->atominfo()->string_to_Z("P"),0,0,0);
609
610 perlout << "%basissets = (" << endl;
611 int nbasis = keyval->count("basislist");
612 Ref<KeyVal> nullkv = new AssignedKeyVal();
613 for (i=0; i<nbasis; i++) {
614 int first_element = 1;
615 char *basisname = keyval->pcharvalue("basislist",i);
616 perlout << " \"" << basisname << "\" => (";
617 BasisFileSet bfs(nullkv);
618 Ref<KeyVal> basiskv = bfs.keyval(nullkv, basisname);
619 char elemstr[512];
620 elemstr[0] = '\0';
621 int last_elem_exists = 0;
622 int n0 = 0;
623 int n1 = 0;
624 int n2 = 0;
625 for (j=0; j<nelem; j++) {
626 Ref<AssignedKeyVal> atombaskv_a(new AssignedKeyVal());
627 Ref<KeyVal> atombaskv(atombaskv_a);
628 char keyword[256];
629 strcpy(keyword,":basis:");
630 strcat(keyword,hmol->atominfo()->name(j+1).c_str());
631 strcat(keyword,":");
632 strcat(keyword,basisname);
633 if (basiskv->exists(keyword)) {
634 if (!first_element) {
635 perlout << ",";
636 }
637 else {
638 first_element = 0;
639 }
640 // This will print the symbol:
641 //perlout << "\"" << AtomInfo::symbol(j+1) << "\"";
642 // This will print the atomic number:
643 perlout << j+1;
644 if (!last_elem_exists) {
645 if (elemstr[0] != '\0') strcat(elemstr,", ");
646 strcat(elemstr,hmol->atominfo()->symbol(j+1).c_str());
647 }
648 else if (last_elem_exists == 2) {
649 strcat(elemstr,"-");
650 }
651 last_elem_exists++;
652 if (j+1 == 1) {
653 atombaskv_a->assign("name", basisname);
654 atombaskv_a->assign("molecule", hmol.pointer());
655 Ref<GaussianBasisSet> gbs=new GaussianBasisSet(atombaskv);
656 n0 = gbs->nbasis();
657 }
658 if (j+1 == 6) {
659 atombaskv_a->assign("name", basisname);
660 atombaskv_a->assign("molecule", cmol.pointer());
661 Ref<GaussianBasisSet> gbs=new GaussianBasisSet(atombaskv);
662 n1 = gbs->nbasis();
663 }
664 if (j+1 == 15) {
665 atombaskv_a->assign("name", basisname);
666 atombaskv_a->assign("molecule", pmol.pointer());
667 Ref<GaussianBasisSet> gbs=new GaussianBasisSet(atombaskv);
668 n2 = gbs->nbasis();
669 }
670 }
671 else {
672 if (last_elem_exists > 1) {
673 if (last_elem_exists == 2) strcat(elemstr,", ");
674 strcat(elemstr, hmol->atominfo()->symbol(j).c_str());
675 }
676 last_elem_exists = 0;
677 }
678 }
679 perlout << ")";
680 if (i != nbasis-1) perlout << "," << endl;
681 perlout << endl;
682 cout << "<tr><td><tt>" << basisname
683 << "</tt><td>" << elemstr << "<td>";
684 if (n0>0) cout << n0;
685 cout << "<td>";
686 if (n1>0) cout << n1;
687 cout << "<td>";
688 if (n2>0) cout << n2;
689 cout << endl;
690 delete[] basisname;
691 }
692
693 perlout << ")" << endl;
694
695 #ifdef HAVE_SSTREAM
696 const char *perlout_s = perlout.str().c_str();
697 #else
698 perlout << ")" << ends;
699 char *perlout_s = perlout.str();
700 #endif
701 cout << perlout_s;
702 }
703
704 return 0;
705 }
706
707 /////////////////////////////////////////////////////////////////////////////
708
709 // Local Variables:
710 // mode: c++
711 // c-file-style: "CLJ"
712 // End:
713