1 //
2 // gaussshval.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 #include <stdlib.h>
29 #include <math.h>
30
31 #include <util/misc/formio.h>
32 #include <util/keyval/keyval.h>
33
34 #include <chemistry/qc/basis/gaussshell.h>
35 #include <chemistry/qc/basis/integral.h>
36 #include <chemistry/qc/basis/cartiter.h>
37 #include <chemistry/qc/basis/transform.h>
38
39 using namespace std;
40 using namespace sc;
41
42 #define MAX_NPRIM 20
43 #define MAX_NCON 10
44 #define MAX_AM 8
45
46 int
values(CartesianIter ** civec,SphericalTransformIter ** sivec,const SCVector3 & r,double * basis_values)47 GaussianShell::values(CartesianIter **civec, SphericalTransformIter **sivec,
48 const SCVector3& r, double* basis_values)
49 {
50 return hessian_values(civec, sivec, r, 0, 0, basis_values);
51 }
52
53 int
grad_values(CartesianIter ** civec,SphericalTransformIter ** sivec,const SCVector3 & r,double * g_values,double * basis_values) const54 GaussianShell::grad_values(CartesianIter **civec,
55 SphericalTransformIter **sivec,
56 const SCVector3& r,
57 double* g_values,
58 double* basis_values) const
59 {
60 return hessian_values(civec, sivec, r, 0, g_values, basis_values);
61 }
62
63 int
hessian_values(CartesianIter ** civec,SphericalTransformIter ** sivec,const SCVector3 & r,double * h_values,double * g_values,double * basis_values) const64 GaussianShell::hessian_values(CartesianIter **civec,
65 SphericalTransformIter **sivec,
66 const SCVector3& r,
67 double* h_values,
68 double* g_values,
69 double* basis_values) const
70 {
71
72 // compute the maximum angular momentum component of the shell
73 int maxam = max_am();
74 if (g_values || h_values) maxam++;
75 if (h_values) maxam++;
76
77 // check limitations
78 if (nprim > MAX_NPRIM || ncon > MAX_NCON || maxam >= MAX_AM) {
79 ExEnv::err0() << indent
80 << "GaussianShell::grad_values: limit exceeded:\n"
81 << indent
82 << scprintf(
83 "ncon = %d (%d max) nprim = %d (%d max) maxam = %d (%d max)\n",
84 ncon,MAX_NCON,nprim,MAX_NPRIM,maxam,MAX_AM-1);
85 abort();
86 }
87
88 // loop variables
89 int i,j;
90
91 // precompute powers of x, y, and z
92 double xs[MAX_AM];
93 double ys[MAX_AM];
94 double zs[MAX_AM];
95 xs[0] = ys[0] = zs[0] = 1.0;
96 if (maxam>0) {
97 xs[1] = r[0];
98 ys[1] = r[1];
99 zs[1] = r[2];
100 }
101 for (i=2; i<=maxam; i++) {
102 xs[i] = xs[i-1]*r[0];
103 ys[i] = ys[i-1]*r[1];
104 zs[i] = zs[i-1]*r[2];
105 }
106
107 // precompute r*r
108 double r2;
109 if (maxam<2) {
110 r2 = 0.0;
111 for (i=0; i<3; i++) {
112 r2+=r[i]*r[i];
113 }
114 }
115 else {
116 r2 = xs[2] + ys[2] + zs[2];
117 }
118
119 // precompute exponentials
120 double exps[MAX_NPRIM];
121 for (i=0; i<nprim; i++) {
122 exps[i]=::exp(-r2*exp[i]);
123 }
124
125 // precompute contractions over exponentials
126 double precon[MAX_NCON];
127 for (i=0; i<ncon; i++) {
128 precon[i] = 0.0;
129 for (j=0; j<nprim; j++) {
130 precon[i] += coef[i][j] * exps[j];
131 }
132 }
133
134 // precompute contractions over exponentials with exponent weighting
135 double precon_g[MAX_NCON];
136 if (g_values || h_values) {
137 for (i=0; i<ncon; i++) {
138 precon_g[i] = 0.0;
139 for (j=0; j<nprim; j++) {
140 precon_g[i] += exp[j] * coef[i][j] * exps[j];
141 }
142 precon_g[i] *= 2.0;
143 }
144 }
145
146 // precompute contractions over exponentials with exponent^2 weighting
147 double precon_h[MAX_NCON];
148 if (h_values) {
149 for (i=0; i<ncon; i++) {
150 precon_h[i] = 0.0;
151 for (j=0; j<nprim; j++) {
152 precon_h[i] += exp[j] * exp[j] * coef[i][j] * exps[j];
153 }
154 precon_h[i] *= 4.0;
155 }
156 }
157
158 // compute the shell values
159 int i_basis=0; // Basis function counter
160 if (basis_values) {
161 for (i=0; i<ncon; i++) {
162 // handle s functions with a special case to speed things up
163 if (l[i] == 0) {
164 basis_values[i_basis] = precon[i];
165 i_basis++;
166 }
167 else if (!puream[i]) {
168 CartesianIter *jp = civec[l[i]];
169 CartesianIter& j = *jp;
170 for (j.start(); j; j.next()) {
171 basis_values[i_basis] = xs[j.a()]*ys[j.b()]*zs[j.c()]
172 *precon[i];
173 i_basis++;
174 }
175 }
176 else {
177 double cart_basis_values[((MAX_AM+1)*(MAX_AM+2))/2];
178 CartesianIter *jp = civec[l[i]];
179 CartesianIter& j = *jp;
180 int i_cart = 0;
181 for (j.start(); j; j.next()) {
182 cart_basis_values[i_cart] = xs[j.a()]*ys[j.b()]*zs[j.c()]
183 *precon[i];
184 i_cart++;
185 }
186 SphericalTransformIter *ti = sivec[l[i]];
187 int n = ti->n();
188 memset(&basis_values[i_basis], 0, sizeof(double)*n);
189 for (ti->start(); ti->ready(); ti->next()) {
190 basis_values[i_basis + ti->pureindex()]
191 += ti->coef() * cart_basis_values[ti->cartindex()];
192 }
193 i_basis += n;
194 }
195 }
196 }
197
198 // compute the gradient of the shell values
199 if (g_values) {
200 int i_grad=0; // Basis function counter
201 for (i=0; i<ncon; i++) {
202 // handle s functions with a special case to speed things up
203 if (l[i] == 0) {
204 double norm_precon_g = precon_g[i];
205 g_values[i_grad] = -xs[1]*norm_precon_g;
206 i_grad++;
207 g_values[i_grad] = -ys[1]*norm_precon_g;
208 i_grad++;
209 g_values[i_grad] = -zs[1]*norm_precon_g;
210 i_grad++;
211 }
212 else if (!puream[i]) {
213 CartesianIter *jp = civec[l[i]];
214 CartesianIter& j = *jp;
215 for (j.start(); j; j.next()) {
216 double norm_precon = precon[i];
217 double norm_precon_g = precon_g[i];
218 g_values[i_grad] = - norm_precon_g
219 * xs[j.a()+1] * ys[j.b()] * zs[j.c()];
220 if (j.a()) g_values[i_grad] += j.a() * norm_precon
221 * xs[j.a()-1] * ys[j.b()] * zs[j.c()];
222 i_grad++;
223
224 g_values[i_grad] = - norm_precon_g
225 * xs[j.a()] * ys[j.b()+1] * zs[j.c()];
226 if (j.b()) g_values[i_grad] += j.b() * norm_precon
227 * xs[j.a()] * ys[j.b()-1] * zs[j.c()];
228 i_grad++;
229
230 g_values[i_grad] = - norm_precon_g
231 * xs[j.a()] * ys[j.b()] * zs[j.c()+1];
232 if (j.c()) g_values[i_grad] += j.c() * norm_precon
233 * xs[j.a()] * ys[j.b()] * zs[j.c()-1];
234 i_grad++;
235 }
236 }
237 else {
238 double cart_g_values[3*((MAX_AM+1)*(MAX_AM+2))/2];
239 CartesianIter *jp = civec[l[i]];
240 CartesianIter& j = *jp;
241 int i_cart = 0;
242 for (j.start(); j; j.next()) {
243 double norm_precon = precon[i];
244 double norm_precon_g = precon_g[i];
245 cart_g_values[i_cart] = - norm_precon_g
246 * xs[j.a()+1] * ys[j.b()] * zs[j.c()];
247 if (j.a()) cart_g_values[i_cart] += j.a() * norm_precon
248 * xs[j.a()-1] * ys[j.b()] * zs[j.c()];
249 i_cart++;
250
251 cart_g_values[i_cart] = - norm_precon_g
252 * xs[j.a()] * ys[j.b()+1] * zs[j.c()];
253 if (j.b()) cart_g_values[i_cart] += j.b() * norm_precon
254 * xs[j.a()] * ys[j.b()-1] * zs[j.c()];
255 i_cart++;
256
257 cart_g_values[i_cart] = - norm_precon_g
258 * xs[j.a()] * ys[j.b()] * zs[j.c()+1];
259 if (j.c()) cart_g_values[i_cart] += j.c() * norm_precon
260 * xs[j.a()] * ys[j.b()] * zs[j.c()-1];
261 i_cart++;
262 }
263 SphericalTransformIter *ti = sivec[l[i]];
264 int n = ti->n();
265 memset(&g_values[i_grad], 0, sizeof(double)*n*3);
266 for (ti->start(); ti->ready(); ti->next()) {
267 double coef = ti->coef();
268 int pi = ti->pureindex();
269 int ci = ti->cartindex();
270 for (int xyz=0; xyz<3; xyz++) {
271 g_values[i_grad + pi*3 + xyz]
272 += coef * cart_g_values[ci*3 + xyz];
273 }
274 }
275 i_grad += 3*n;
276 }
277 }
278 }
279
280 // compute the hessian of the shell values
281 if (h_values) {
282 int i_hess=0; // Basis function counter
283 for (i=0; i<ncon; i++) {
284 // handle s functions with a special case to speed things up
285 if (l[i] == 0) {
286 double norm_precon_g = precon_g[i];
287 double norm_precon_h = precon_h[i];
288 // xx
289 h_values[i_hess] = norm_precon_h*xs[2] - norm_precon_g;
290 i_hess++;
291 // yx
292 h_values[i_hess] = norm_precon_h*xs[1]*ys[1];
293 i_hess++;
294 // yy
295 h_values[i_hess] = norm_precon_h*ys[2] - norm_precon_g;
296 i_hess++;
297 // zx
298 h_values[i_hess] = norm_precon_h*zs[1]*xs[1];
299 i_hess++;
300 // zy
301 h_values[i_hess] = norm_precon_h*zs[1]*ys[1];
302 i_hess++;
303 // zz
304 h_values[i_hess] = norm_precon_h*zs[2] - norm_precon_g;
305 i_hess++;
306 }
307 else {
308 double *cart_h;
309 double tmp_cart_h[6*((MAX_AM+1)*(MAX_AM+2))/2];
310 if (!puream[i]) {
311 cart_h = &h_values[i_hess];
312 }
313 else {
314 cart_h = tmp_cart_h;
315 }
316 CartesianIter *jp = civec[l[i]];
317 CartesianIter& j = *jp;
318 int i_cart = 0;
319 for (j.start(); j; j.next()) {
320 double pre = precon[i];
321 double pre_g = - precon_g[i];
322 double pre_h = precon_h[i];
323 int a = j.a();
324 int b = j.b();
325 int c = j.c();
326 // xx
327 cart_h[i_cart] = pre_h * xs[a+2]*ys[b]*zs[c]
328 + pre_g * (a+1) * xs[a]*ys[b]*zs[c];
329 if (a>0) {
330 cart_h[i_cart] += pre_g * a*xs[a]*ys[b]*zs[c];
331 if (a>1) cart_h[i_cart] += pre * a*(a-1)
332 * xs[a-2]*ys[b]*zs[c];
333 }
334 i_cart++;
335
336 // yx
337 cart_h[i_cart] = pre_h * xs[a+1]*ys[b+1]*zs[c];
338 if (a>0)
339 cart_h[i_cart] += pre_g * a * xs[a-1]*ys[b+1]*zs[c];
340 if (b>0)
341 cart_h[i_cart] += pre_g * b * xs[a+1]*ys[b-1]*zs[c];
342 if (a>0 && b>0)
343 cart_h[i_cart] += pre * a*b * xs[a-1]*ys[b-1]*zs[c];
344 i_cart++;
345
346 // yy
347 cart_h[i_cart] = pre_h * xs[a]*ys[b+2]*zs[c]
348 + pre_g * (b+1) * xs[a]*ys[b]*zs[c];
349 if (b>0) {
350 cart_h[i_cart] += pre_g * b*xs[a]*ys[b]*zs[c];
351 if (b>1) cart_h[i_cart] += pre * b*(b-1)
352 * xs[a]*ys[b-2]*zs[c];
353 }
354 i_cart++;
355
356 // zx
357 cart_h[i_cart] = pre_h * xs[a+1]*ys[b]*zs[c+1];
358 if (a>0)
359 cart_h[i_cart] += pre_g * a * xs[a-1]*ys[b]*zs[c+1];
360 if (c>0)
361 cart_h[i_cart] += pre_g * c * xs[a+1]*ys[b]*zs[c-1];
362 if (a>0 && c>0)
363 cart_h[i_cart] += pre * a*c * xs[a-1]*ys[b]*zs[c-1];
364 i_cart++;
365
366 // zy
367 cart_h[i_cart] = pre_h * xs[a]*ys[b+1]*zs[c+1];
368 if (c>0)
369 cart_h[i_cart] += pre_g * c * xs[a]*ys[b+1]*zs[c-1];
370 if (b>0)
371 cart_h[i_cart] += pre_g * b * xs[a]*ys[b-1]*zs[c+1];
372 if (c>0 && b>0)
373 cart_h[i_cart] += pre * c*b * xs[a]*ys[b-1]*zs[c-1];
374 i_cart++;
375
376 // zz
377 cart_h[i_cart] = pre_h * xs[a]*ys[b]*zs[c+2]
378 + pre_g * (c+1) * xs[a]*ys[b]*zs[c];
379 if (c>0) {
380 cart_h[i_cart] += pre_g * c*xs[a]*ys[b]*zs[c];
381 if (c>1) cart_h[i_cart] += pre * c*(c-1)
382 * xs[a]*ys[b]*zs[c-2];
383 }
384 i_cart++;
385 }
386 if (puream[i]) {
387 SphericalTransformIter *ti = sivec[l[i]];
388 int n = ti->n();
389 memset(&h_values[i_hess], 0, sizeof(double)*n*6);
390 for (ti->start(); ti->ready(); ti->next()) {
391 double coef = ti->coef();
392 int pi = ti->pureindex();
393 int ci = ti->cartindex();
394 for (int xyz2=0; xyz2<6; xyz2++) {
395 h_values[i_hess + pi*6 + xyz2]
396 += coef * cart_h[ci*6 + xyz2];
397 }
398 }
399 i_hess += 6*n;
400 }
401 else {
402 i_hess += 3*(l[i]+1)*(l[i]+2);
403 }
404 }
405 }
406 }
407
408 return i_basis;
409 }
410
411 int
test_monobound(double & r,double & bound) const412 GaussianShell::test_monobound(double &r, double &bound) const
413 {
414 // compute the maximum angular momentum component of the shell
415 // add one since derivatives will be needed
416 int maxam = max_am() + 1;
417
418 // check limitations
419 if (nprim > MAX_NPRIM || ncon > MAX_NCON || maxam >= MAX_AM) {
420 ExEnv::err0() << indent
421 << "GaussianShell::gaussshval: limit exceeded:\n"
422 << indent
423 << scprintf(
424 "ncon = %d (%d max) nprim = %d (%d max) maxam = %d (%d max)\n",
425 ncon,MAX_NCON,nprim,MAX_NPRIM,maxam,MAX_AM-1);
426 abort();
427 }
428
429 // loop variables
430 int i,j;
431
432 // precompute powers of r
433 double rs[MAX_AM+1];
434 rs[0] = 1.0;
435 if (maxam>0) {
436 rs[1] = r;
437 }
438 for (i=2; i<=maxam; i++) {
439 rs[i] = rs[i-1]*r;
440 }
441
442 // precompute r*r
443 double r2 = r*r;
444
445 // precompute exponentials
446 double exps[MAX_NPRIM];
447 for (i=0; i<nprim; i++) {
448 exps[i]=::exp(-r2*exp[i]);
449 }
450
451 // precompute contractions over exponentials
452 double precon[MAX_NCON];
453 for (i=0; i<ncon; i++) {
454 precon[i] = 0.0;
455 for (j=0; j<nprim; j++) {
456 // using fabs since we want a monotonically decreasing bound
457 precon[i] += fabs(coef[i][j]) * exps[j];
458 }
459 }
460
461 // precompute contractions over exponentials with exponent weighting
462 double precon_w[MAX_NCON];
463 for (i=0; i<ncon; i++) {
464 precon_w[i] = 0.0;
465 for (j=0; j<nprim; j++) {
466 precon_w[i] += exp[j] * fabs(coef[i][j]) * exps[j];
467 }
468 }
469
470 double max_bound = 0.0;
471 bound = 0.0;
472 for (i=0; i<ncon; i++) {
473 // using r^l since r^l >= x^a y^b z^c
474 double component_bound = rs[l[i]]*precon[i];
475 if (l[i] > 0) {
476 double d1 = -2.0*rs[l[i]+1]*precon_w[i];
477 double d2 = l[i]*rs[l[i]-1]*precon[i];
478 if (d1+d2 > 0) {
479 // This bound is no good since the contraction is increasing
480 // at this position. Move r out and return to let the driver
481 // call again.
482 double rold = r;
483 r = sqrt(l[i]*precon[i]/(2.0*precon_w[i]));
484 if (r<rold+0.01) r = rold+0.01;
485 //ExEnv::outn() << "rejected at " << rold << " trying again at "
486 // << r << endl;
487 return 1;
488 }
489 }
490 if (component_bound > max_bound) {
491 max_bound = component_bound;
492 }
493 }
494
495 bound = max_bound;
496 return 0;
497 }
498
499 double
monobound(double r) const500 GaussianShell::monobound(double r) const
501 {
502 // doesn't work at r <= zero
503 if (r<=0.001) r = 0.001;
504 double b;
505 while (test_monobound(r, b));
506 return b;
507 }
508
509 /////////////////////////////////////////////////////////////////////////////
510
511 // Local Variables:
512 // mode: c++
513 // c-file-style: "CLJ"
514 // End:
515