3
#include <libint/libint.h>
11
/*------------------------------------------------------
12
Use Newton-Raphson method to find maximum r for which
13
the radial parts of basis functions drop below
15
------------------------------------------------------*/
17
void init_rad_extent(double thresh)
21
int shell, prim, first_prim, last_prim, am;
22
const double r0 = 4.0; /* Start at 4.0 bohr ---*/
23
double func, dfuncdr, sum, dsumdr;
26
BasisSet.thresh = thresh;
28
for (shell=0;shell<BasisSet.num_shells;shell++) {
30
first_prim = BasisSet.shells[shell].fprim - 1;
31
last_prim = first_prim + BasisSet.shells[shell].n_prims - 1;
32
am = BasisSet.shells[shell].am-1;
37
/*--------------------------------------
38
evaluate the radial part of the shell
39
and its first derivative
40
--------------------------------------*/
43
for(prim=first_prim;prim<=last_prim;prim++) {
44
tmp = BasisSet.cgtos[prim].ccoeff[am]*exp(-BasisSet.cgtos[prim].exp*(r*r));
48
dfuncdr += -2.0*r*BasisSet.cgtos[prim].exp*tmp;
52
dfuncdr += (1.0-2.0*r*r*BasisSet.cgtos[prim].exp)*tmp;
58
dfuncdr += (1.0-2.0*r*r*BasisSet.cgtos[prim].exp)*tmp;
62
/*--- We actually want to look at the
63
absolute value of the basis function ---*/
64
if (func+thresh < 0.0) {
65
func = -(func+thresh)-thresh;
69
/*--- If the function is too tight - reduce r ---*/
70
if (func+thresh == 0.0 && dfuncdr == 0.0) {
75
/*--- If before or at the maximum (for p-functions and higher)
76
step hopefully past it ---*/
82
r_new = r - func/dfuncdr;
88
punt("Too many iterations while computing radial extents");
90
} while (fabs(func/thresh) >= SOFT_ZERO);
92
BasisSet.shells[shell].rad_extent = r_new;
93
if (UserOptions.print_lvl > PRINT_DEBUG)
94
fprintf(outfile,"Shell# = %d Radial extent = %lf\n",shell,r_new);