Subject: Re: CCL:G:OVERLAP INTEGRAL AND SLATER FUNNC

From: Herbert Homeier t4720 <[email-masked]>
Date: Mon, 3 Nov 97 17:59:30 +0100

Hello,

> Date: Sun, 2 Nov 1997 16:54:11 -0500 (EST)
> From: "E. Lewars" <[email-masked]>
> Message-Id: <[email-masked]>
> To: [email-masked]
> Subject: CCL:G:OVERLAP INTEGRAL AND SLATER FUNNC
> 
> 1997 Nov 2
> 
> SLATER OVERLAP INTEGRALS
> >From E. Lewars
> To CCL
> 
> Hello,
> 
> There is some info' about overlap integrals and Slater functions that I am
> finding it hard to dig out of the literature. I should state at the outset
> that _I know_:
> 
> that there are programs for calculating overlap integrals
> 
> that in most molecular quantum chem Gaussian functions, not Slater
> functions,
> are used
> 
> that Mulliken et al published a long paper (J Chem Phys, 17 (1949)
> 1248-1267))
> with recipes for calculating overlap integrals

that contains quite a few errors ....

> 
> that all kinds of tricks exist for "simplifying" the calculation of overlap
> integrals
> ---------
> OK, but what I want is simply this:
> 
> (1) In the molecule H-He+ (yes, protonated helium) if the H nucleus has
> Cartesians H1(x_1, y_1, z_1) and the He nucleus He2(x_2, y_2, z_2),
> where the coordinates are in Angstroms (_not_ atomic units or Bohrs)
> 
> and if we write
> 
> for H: phi_1 = a_1 exp{b_1[(x-x_1)^2 + (y-y_1)^2 + (z-z_1)^2]^1/2}
> and for He: phi_2 = a_2 exp{b_2[(x-x_2)^2 + (y-y_2)^2 + (z-z_2)^2]^1/2}

The best source of information on such questions is still the
"classical" paper
E. Clementi and C. Roetti, At. Data & Nucl. Data Tables 14 (1974) 177
[maybe in atomic units, though 8^)]

> 
> (usually a is expressed in terms of zeta and pi, b = -zeta, and
> the variable is a radius vector r)
> QUESTION: what are a_1 and b_1, as numbers, not Greek letters?
> what are a_2,and b_2 as numbers, not Greek letters?
> (i.e a_1 = 0.5647 or whatever, etc etc (4 decimals).
> 
> S_12 = the triple integral from 0--> infinity of phi_1.phi_2 dxdydz
> must = the correct overlap, if the six Cartesian coordinates are in
> Angstroms. And of course integral of S_11 = S_22 = 1.
> 
> (2) If (x_1, y_1, z_1) = (0, 0, 0) and (x_2, y_2, z_2) = (0, 0, 0.800),
> (0.800 Angstroms) what is S_12 (as a number, e.g 0.4273 or whatever)?
> 
> (3) For (x_1, y_1, z_1) = (0, 0, 0) and (x_2, y_2, z_2) = (0, 0, d),
> is there a simple function S_12 = f(d) ? (for H, He? for any two
> 1s orbitals)? Of course f will depend parametrically on the a's and b's.
> -------------
> Thanks E. Lewars
> ===============
>

The answer to your final question is given below, however, in somewhat
different notation:

alpha : -b_1
beta  : -b_2
R     : d
my S  : Overlap of exp(-alpha*r_1) and exp(-beta*r_2), ie., for a_1=a_2=1,
        where r_1=sqrt(x**2+y**2+z**2) and
              r_2=sqrt(x**2+y**2+(z-R)**2)

The result of the fortran program listed below contains also the
value where a_1 and a_2 are chosen so that the 1s orbitals are
normalized. See the output.

There is an WWW interface to this:

http://www.chemie.uni-regensburg.de/~hoh05008/form-Soverlap.html

For 1s functions, of course, the result for the overlap integral
does only depends on a_1,a_2,b_1,b_2 and d, not on the particular
geometry of the centers (only their distance d matters).

As regards units:

If you enter R in Angstroem, and alpha and beta in 1/Angstroem,
everything should be fine in the case of the overlap integral of
normalized 1s functions.

You may be interested also in the following papers:

@ARTICLE{HomeierSteinborn92ote,
        author={H. H. H. Homeier and E. O. Steinborn},
        title={On the
Evaluation of Overlap Integrals with Exponential--type Basis
Functions},
        journal="Int. J. Quantum Chem.",
        volume=42,
        pages="761--778",
        year="1992"
}

@ARTICLE{HomeierWenigerSteinborn92pft,
        author={H. H. H. Homeier and E. J. Weniger and E. O. Steinborn},
        title={Programs for the Evaluation of Overlap Integrals with
{$B$} Functions},
        journal="Comput. Phys. Commun.",
        volume=72,
        pages="269--287",
        year="1992"
}

The latter describes a fortran program package S_INT for more general orbitals.
The output generated by the programs given below has been checked against
the output of this package.

One should be careful (!), however, with the closed form expression
coded in the programs below in the case that alpha and beta are not
equal, but very close, and in the case that R is small but not exactly
zero. The reason is that in these cases, cancelling singularities turn
up. For instance, one gets a result somewhat greater than one for the
normalized functions for alpha=1.000 beta=1.001 R=0.01 and this, of
course is wrong (overlaps of normalized functions have to be smaller
than 1.

The package S_INT generates accurate output even in the case described above.

Best regards

Herbert Homeier

---------------- C program -------------------------
double S(alpha,beta,R)
double alpha;
double beta;
double R;
{
  if (R == 0)
    return(8*0.3141592653589793E1/pow(alpha+beta,3));
  else 
    if (alpha != beta)
      return(8*0.3141592653589793E1*(-4*alpha*beta+4*exp(-(alpha-beta)*R)*alpha
*beta-exp(-(alpha-beta)*R)*R*beta*beta*beta+R*alpha*alpha*alpha-R*alpha*beta*
beta+exp(-(alpha-beta)*R)*R*alpha*alpha*beta)*exp(-R*beta)/R/pow(alpha+beta,3)/
pow(alpha-beta,3));
    else 
      return((R*R*alpha*alpha/3+1+R*alpha)*0.3141592653589793E1*exp(-R*alpha)/(
alpha*alpha*alpha));
}

--------------  maple procedure
S:=
proc(alpha, beta, R)
    if R = 0 then 8*Pi/(alpha + beta)^3
    elif alpha <> beta then 8*Pi*(-4*alpha*beta + 4*exp(- (alpha - beta)*R)*alpha*beta
         - exp(- (alpha - beta)*R)*R*beta^3 + R*alpha^3 - R*alpha*beta^2
         + exp(- (alpha - beta)*R)*R*alpha^2*beta)*exp(- R*beta)
           /(R*(alpha + beta)^3*(alpha - beta)^3)
    else (1/3*R^2*alpha^2 + 1 + R*alpha)*Pi*exp(- R*alpha)/alpha^3
    fi
end;

--------------- f77 program including driver and overlap integral for normalized orbitals
      doubleprecision function S(alpha,beta,R)
      doubleprecision alpha
      doubleprecision beta
      doubleprecision R
 
        if (R .eq. 0.D0) then
          S = 8.D0*0.3141592653589793D1/(alpha+beta)**3.D0
        else 
          if (alpha .ne. beta) then
            S = 8.D0*0.3141592653589793D1*(-4.D0*alpha*beta+4.D0*exp(-(a
     #lpha-beta)*R)*alpha*beta-exp(-(alpha-beta)*R)*R*beta**3+R*alpha**3
     #-R*alpha*beta**2+exp(-(alpha-beta)*R)*R*alpha**2*beta)*exp(-R*beta
     #)/R/(alpha+beta)**3.D0/(alpha-beta)**3.D0
          else 
            S = (R**2*alpha**2/3.D0+1.D0+R*alpha)*0.3141592653589793D1*e
     #xp(-R*alpha)/alpha**3
          endif
        endif
      return
      end
C
      program mainS
      double precision a,b,r,S,SS,SSS
      external S
      read(*,*,ERR=99999,END=99999) a,b,r
      if(a.le.0.d0) goto 99999
      if(b.le.0.d0) goto 99999
      if(R.lt.0.d0) goto 99999
      SS=S(a,b,r)
      SSS=SS/0.3141592653589793D1*DSQRT(a*b)**3
      write(*,*) 'int(exp(-',a,'*sqrt(x**2+y**2+z**2))'
      write(*,*) '   *exp(-',b,'*sqrt(x**2+y**2+(z-',r,')**2))'
      write(*,*) '    dx dy dz ) = '
      write(*,*) ' = ',SS
      write(*,*)
      write(*,*) 'int(exp(-',a,'*sqrt(x**2+y**2+z**2))'
      write(*,*) '   *exp(-',b,'*sqrt(x**2+y**2+(z-',r,')**2))'
      write(*,*) '    dx dy dz ) /Pi *(',a,'*',b,')**(3/2) = '
      write(*,*) ' = ',SSS
      stop
99999 write(*,*) 'Input ends or errs'
      stop 
      end

--------------------------------------------------------------
Priv.-Doz. Dr. Herbert H. H. Homeier
Institut fuer Physikalische und Theoretische Chemie
Universitaet Regensburg, D-93040 Regensburg, Germany
Phone: +49-941-943 4720  FAX: +49-941-943 4719/+49-941-943 2305
email: [email-masked]
WWW: http://www.chemie.uni-regensburg.de/~hoh05008