From jkl@ccl.net  Fri Aug 18 12:48:59 1995
Received: from bedrock.ccl.net  for jkl@ccl.net
	by www.ccl.net (8.6.10/950810.1506) id MAA23949; Fri, 18 Aug 1995 12:35:45 -0400
Received: from hrz-sun1.hrz.uni-kassel.de  for gdanitz@hrz.uni-kassel.de
	by bedrock.ccl.net (8.6.10/930601.1506) id MAA16007; Fri, 18 Aug 1995 12:34:54 -0400
Received: from hrz-serv1.hrz.uni-kassel.de by hrz-sun1.hrz.uni-kassel.de (4.1/SMI-4.1)
	id AA19695; Fri, 18 Aug 95 18:34:28 +0200
Received: by hrz-serv1.hrz.uni-kassel.de (AIX 3.2/UCB 5.64/HRZ-GhK/HRZ-SERV1/pm)
          id AA21649; Fri, 18 Aug 1995 18:34:46 +0200
From: gdanitz@hrz.uni-kassel.de (Robert Gdanitz)
Message-Id: <9508181634.AA21649@hrz-serv1.hrz.uni-kassel.de>
Subject: Num. stability of 4-index transf. in CI
To: chemistry@ccl.net (Computational Chemistry List)
Date: Fri, 18 Aug 1995 18:34:45 +0200 (MES)


+-----+  Robert J. Gdanitz                      email: gdanitz@hrz.uni-kassel.de
| GhK |  Gesamthochschule Kassel                Tel.: +(49) 561-804-4120
|     |  Fachbereich 18 (Physik)                Fax:  +(49) 561-804-4006
+-----+  34109 Kassel 
X-Mailer: ELM [version 2.4 PL24 PGP3 *ALPHA*]
MIME-Version: 1.0
Content-Type: text/plain; charset=ISO-8859-1
Content-Transfer-Encoding: 8bit
Content-Length: 2381      

Dear Netters:

In CI-calculations one faces the problem of transforming 2e-integrals over
(primitive) basis functions (ij|g|kl) to an atomic or a molecular orbital basis
yielding (IJ|g|KL) via the (in)famous 4-index transformation

(IJ|g|KL) =  Sum i,j,k,l  U(i,I) * U(j,J) * U(k,K) * U(l,L) * (ij|g|kl)

Clearly, numerical stability problems are to be expected if the orbital expan-
sion coefficients U "oscillate" with an amplitude that reaches the value of the
4-th root of the reciprocal floating point accuracy, i.e. approx. 1000 in 
REAL*8. Unfortunately this value (corresponding to a condition number of 1.0E8
of the overlap matrix in the primitive basis set) is quickly reached in the case
of a) rather large basis sets, b) diffuse functions for molecules or c) a low or
even absent symmetry of the system under investigation and the results may be 
spurious states that mix with the one to be optimized or even states with almost
arbitrary negative values in energy.

A few possible solutions that spontaneously came to my mind:

* Promoting everything to REAL*16 might be accurate enough to yield numerical 
  stability, however, REAL*16 is not well supported by computer manufacturers 
  and its usage may easily cause a degradation in performance of orders of 
  magnitude (on our SG Power-Challenge, a test program computing dot-products of
  two vectors was 350 times slower than in REAL*8!)
* Only accumulating the partial sums in REAL*16 might help, however, this is,
  where almost all numerical work is performed, so it should not be much faster
  than doing everything in REAL*16. In terms of stability, this option should be
  worse compared to the first, however, it is implemented much easier.
* Accumulating the partial sums in the order of descending absolute magnitude
  increases CPU-time by a constant factor, depending on the algorithm used for
  the ranking/sorting. The benefits and drawbacks seem to be similar to option
  two.

Unfortunately it seems to be impossible to create a significant (and small) 
test-case such that a primitive program could be used to find the best 
solution. 

I would like to encourage anyone who has comments or suggestions, especially
those, who already have experience with performing integral transformations in 
relatively large basis sets, to write me! If there is enough interest, I will 
summarize.

R. Gdanitz