/***************************************************************************
                          freeMQM-1.cpp  -  description
                             -------------------
    begin                : Thu Sept 17 2003
    copyright            : (C) 2003 by Brian J. Keay
    email                : ulisse2k@yahoo.com
 ***************************************************************************/

/***************************************************************************
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 ***************************************************************************/
#include <string>
#include <cctype>
#include <iostream>
#include <fstream>
#include <iomanip>
#include <ctime>
//#include <vector>
//#include <algorithm>
#include <math.h>
#include <gsl/gsl_sf.h>
#include <gsl/gsl_cblas.h> 
#include <gsl/gsl_sf_erf.h>
#include "freeMQM-1.h"
#include "../routines/Lebedev-Laikov.h"

const double PI = 4.0*atan(1.0) ;

using namespace std;


int main(void){
  
   vector<double> Alpha,d;
   vector<int> l,m,n,Ncontract;
   int functional;
ReadBasisSets(Alpha,d,l,m,n,Ncontract,functional);
int time0 = time(0);
Vec_DP C(N);
  double Energy, OldEner;
  const double Precision = 1.0e-11;
// Initialise the Coefficients
  for(int I = 0 ; I <= N-1 ; I++) C[I] = 1.0;

//////////////  Calculate all the overlap integrals and form the S matrix ///
   Mat_DP SMatr(Nbasis,Nbasis);
     for(int I = 0 ; I <= Nbasis-1 ; I++){
cout << "Alpha["<<I<<"]  =   "<< Alpha[I] << "  d["<<I<<"]  =   "<< d[I]<< "  l["<<I<<"]  =   "<< l[I]<< "  m["<<I<<"]  =   "<< m[I]<< "  n["<<I<<"]  =   "<< n[I]<<   endl;
}
   Smatrix(SMatr,Alpha,l,m,n);
//   Form the contracted S matrix   
   Mat_DP SMatrix(MSize,MSize);
   S_contract(SMatrix,SMatr,Alpha,d,l,m,n,Ncontract);
//   Calculate the V matrix that diagonalizes the S matrix 
   Mat_DP VMatrix(MSize,MSize);
   CalcVMatrix(VMatrix,SMatrix,MSize);
//   Calculate the kinetic matrix in the noncontracted basis set.
  Mat_DP Kinet2(Nbasis,Nbasis);
   KE(Kinet2,Alpha,l,m,n);
//   Calculate the electron-nucleous coulomb attraction matrix in the noncontracted basis set.    
   Mat_DP Coul2(Nbasis,Nbasis);
   e_nucl(Coul2,Alpha,l,m,n);
//   Contract the basis set for the electron-nucleous coulomb attraction matrix and 
//    the kinetic matrix and add the result to form the "core" Hamiltonian matrix.   

   Mat_DP HMatrix(MSize,MSize),KE_Contract(MSize,MSize),e_nucl_Contract(MSize,MSize);
   CalcHMatrix(HMatrix,KE_Contract,e_nucl_Contract,Kinet2,Coul2,Alpha,d,l,m,n,Ncontract);
   
   Mat4D_DP QMatrix2(Nbasis,Nbasis,Nbasis,Nbasis);
   ERI(QMatrix2,Alpha,l,m,n);
   Mat4D_DP QMatrix(Nbasis,Nbasis,Nbasis,Nbasis);
  ERI_Contract(QMatrix,QMatrix2,Alpha,d,l,m,n,Ncontract);  // 02 June 2003
      OldEner = -1.0;
      Energy = 0.0;

        Mat_DP GMatrix(MSize,MSize),JMatrix(MSize,MSize),KMatrix(MSize,MSize),
        FMatrix(MSize,MSize),DensMat(MSize,MSize);

        vector<double> Diag(MSize);
 //       double energyXC=0.0;
        double energyX=0.0;
        double energyC=0.0;
                             
//C Self-consistency loop:
       while (abs(OldEner-Energy) > Precision) {
      
        OldEner = Energy;

       if ((method==2)||(method==3)) BuildG(GMatrix,DensMat,QMatrix);
       if (method==1)  BuildG4(GMatrix,JMatrix,KMatrix,C,QMatrix);
        CalcFockMatrix(FMatrix, HMatrix, GMatrix);
         cout << " time =  " << time(0) - time0 <<" seconds " <<  (time(0) - time0 )/60.0 <<" minutes " << endl;
        if (method==2) R_Integrate(energyX, energyC,FMatrix, C, Alpha, d, l, m, n, Ncontract,functional);
        if (method==3)  xyz_Integrate(energyX, energyC,FMatrix, C, Alpha, d, l, m, n, Ncontract,functional);
        DiagFock( C, Diag,VMatrix,  FMatrix, SMatrix, MSize);
        Energy = Diag[0];
cout << "DiagFock Energy ********  "<< Energy << endl;
         BuildDensMat(DensMat, C, SMatrix);

//     CheckVMatrix();
DisplayEnergies(Energy, energyX,energyC,Diag, DensMat, HMatrix,
 GMatrix, FMatrix, JMatrix, KMatrix, C, KE_Contract, e_nucl_Contract);
} // End while

 cout << " time =  " << time(0) - time0 <<" seconds " <<  (time(0) - time0 )/60.0 <<" minutes " << endl;
 return 0;
 }

void ReadBasisSets(vector<double> &Alpha,vector<double> &d,vector<int> &l,vector<int> &m,vector<int> &n,vector<int> &Ncontract,int &functional){
int dum,basis;
 int count = 0;
// int count2 = 0;
 string junk,basisfile,element,text;   

//------Begin ReadBasisSets----------


 cout << "Do you want to model the Helium atom or the Hydrogen molecule?" << endl
  << "Type He for Helium atom or H for the Hydrogen molecule" <<endl;
  cin >> element;

  cout << "What type of calculation do you want to do?" << endl
  << "1) For the Hartree-Fock calculation type the number 1."<< endl
  << "2) For a Density Functional calculation type the number 2." <<endl;
    string Method;
getline(cin,Method);
getline(cin,Method);
  if (Method=="")  method=2;
  else method =atoi(Method.c_str());
  cout << " method = " << method << endl;

  if (method==2){
  cout << "Which functional do you want to use?" << endl
  << "1) For the BLYP exchange-correlation functional just press the return/enter key."<< endl
  << "2) For the Slater-LDA exchange functional type the number 2."<< endl
  << "3) For the Becke B88 exchange functional type the number 3."<< endl
  << "4) For the B88/VWNRPA exchange-correlation functional type the number 4."<< endl
//  << "4) For the Becke B88 exchange functional type the number 4."<< endl
//  << "5) For the LYP correlation functional type the number 5."<< endl
//  << "6) For the  B88-VWNRPA (BVWNRPA) exchange-correlation functional type the number 6." <<endl
//  << "4) To use a hybrid  B88+0.81*LYP+0.19*VWNRPA exchange-correlation functional type the number 4." <<endl
  << "5) To use the Becke (B88) exchange functional with the..." << endl
  << "...Perdew-Zunger (PZ81) correlation functional type the number 5." <<endl
  << "6) To use the Perdew-Burke-Ernzerhof (PBE) exchange functional with the..." << endl
  << "...Perdew-Zunger (PZ81) correlation functional type the number 6." <<endl;
    string Functional;
//getline(cin,Functional);
getline(cin,Functional);
  if (Functional=="")  functional=1;
  else functional =atoi(Functional.c_str());
  cout << " functional = " << functional << endl;

  }
  
  cout << "What type of basis set do you want to use?" << endl
  << "1) To use an noncontracted basis set type the number 1."<< endl
  << "2) To use 6-31G* type the number 2." <<endl
  << "3) To use 6-31G** type the number 3 (and see Revision Note 0.112)." <<endl;
  cin >> basis;

  element[0]= toupper(element[0]);
  element[1]= tolower(element[1]);

  if (element=="H") {atoms = 2;
                     Ratom = 1.0;
          if (method==2) method=3;

     cout << "The default internuclear distance for the hydrogen molecule is 1.40 a.u.." << endl
  << " If you would you like to change the distance type it in now?"<< endl
  << "Otherwise just press the return key." <<endl;
  string dist;
getline(cin,dist);
getline(cin,dist);
  if (dist=="")  Dist = 1.400;//61;
  else Dist =atof(dist.c_str());
dx = 0;
dy = 0;
dz = Dist;
  cout << "Dist = " <<Dist << endl;
  }else if (element=="He"){ atoms = 1;
                       Ratom = 0.5882;
                       if ((method==2)&&(basis==3)) method=3;}
                       
  cout << "atoms = " << atoms << "  Ratom = " << Ratom <<endl;
//----------------------------------
  if (basis==1) basisfile = "-noncontracted.txt";
  else if (basis==2) basisfile = "-6-31Gx.txt";
  else if (basis==3) basisfile = "-6-31Gxx.txt";
  else if (basis==4) basisfile = "-test.txt";
  basisfile = "../Basis_sets/"+element+basisfile;

  cout <<" The basis file that is being read is "<<basisfile<< endl;

 ifstream fp(basisfile.c_str());

 if(!fp) cout << "Cannot open output file.\n";

   

   fp >>Z[0]>> cbasis[0]>>ubasis[0]>>NS>>N2S>>NP>>NL6>>NL3>>ND;
  cout << "Z[0] = "<< Z[0]<< " cbasis[0] = "<< cbasis[0] << " ubasis[0] = "<< ubasis[0] << " NS = "<< NS<< " N2S = "<< N2S<< " NP  = "<< NP<< " NL6  =  "<< NL6<< " NL3  =  "<< NL3<<   endl;
  MSize=N;  // For atoms, must be changed for molecules

Alpha.resize(ubasis[0],0);
d.resize(ubasis[0],0);
l.resize(ubasis[0],0);
m.resize(ubasis[0],0);
    n.resize(ubasis[0],0);
    Ncontract.resize(cbasis[0],1);
    if ( NP != 0 ){

    l[NS+N2S]    = 1;
    m[NS+N2S +1] = 1;
    n[NS+N2S +2] = 1;
  }


cout << "Alpha.size() = "<< Alpha.size()<< " Alpha[0] = "<< Alpha[0] << endl;
cout << "l[ NS+N2S] = ";
cout << l[ NS+N2S] <<endl;
int i=0;
element[1]= toupper(element[1]);
int pos;
fp.seekg(0);
NS=0;NP=0;NL=0;ND=0;
cout << "cbasis[0] = "<< cbasis[0] << " ubasis[0] = "<< ubasis[0] << " NS = "<< NS<< " NP  = "<< NP<< " NL  =  "<< NL<<  " ND  =  "<< ND<<   endl;
while (!fp.eof()){

 pos = fp.tellg();
getline(fp,text);

      if (text.length()<2) continue;
     if ((text.substr(1,1)=="S") && (text.length() < 7 )){
     fp.seekg(pos);
     fp >> junk >> Ncontract[i];
       
        NS+= Ncontract[i];
     
     } else if ((text.substr(1,1)=="P") && (text.length() < 7 )){
       fp.seekg(pos);
       fp >> junk >> Ncontract[i];
       NP+= Ncontract[i];
     } else if ((text.substr(1,1)=="L") && (text.length() < 7 )){
       fp.seekg(pos);
       fp >> junk >> Ncontract[i];
       NL+= Ncontract[i];

     } else if ((text.substr(1,1)=="D") && (text.length() < 7 )){
       fp.seekg(pos);
       fp >> junk >> Ncontract[i];
       ND+= Ncontract[i];
      } else continue;

  for(int I = 0 ; I <= Ncontract[i]-1 ; I++){

    fp >> dum>>Alpha[count] >> d[count];
    cout << "Alpha["<<count<<"] = "<< Alpha[count] << " d["<<count<<"] = "<< d[count]
     << " l["<<count<<"] = "<< l[count] << " m["<<count<<"] = "<< m[count]
     << " n["<<count<<"] = "<< n[count]<<   endl;
    count++;

}
   i++;
   if (count >= ubasis[0]) break;
   if (i >= cbasis[0]) break;
}
cout << " ubasis[0] = " << NS  +3*NP << " cbasis[0] = " << i+2*NP << endl;
cout << "cbasis[0] = "<< cbasis[0] << " ubasis[0] = "<< ubasis[0] << " NS = "<< NS<< " N2S = "<< N2S<< " NP  = "<< NP<< " NL6  =  "<< NL6<< " NL3  =  "<< NL3<< " NL1  =  "<< NL1<< " ND  =  "<< ND<<   endl;
for(int I = NS ; I <= ubasis[0]-1 ; I++){
 Alpha[I]=Alpha[NS];
 d[I]=d[NS];
     cout << "Alpha["<<I<<"] = "<< Alpha[I] << " d["<<I<<"] = "<< d[I]
     << " l["<<I<<"] = "<< l[I] << " m["<<I<<"] = "<< m[I]
     << " n["<<I<<"] = "<< n[I]<<   endl;

  }
// Only considering diatomic molecules at the moment so just double everything.  
  if (atoms==2){

  Alpha.insert (Alpha.end(),Alpha.begin(),Alpha.end());
  d.insert (d.end(),d.begin(),d.end());
  l.insert (l.end(),l.begin(),l.end());
  m.insert (m.end(),m.begin(),m.end());
  n.insert (n.end(),n.begin(),n.end());
  Ncontract.insert (Ncontract.end(),Ncontract.begin(),Ncontract.end());
  Z.insert (Z.end(),Z.begin(),Z.end());
  ubasis.insert (ubasis.end(),ubasis.begin(),ubasis.end());
  cbasis.insert (cbasis.end(),cbasis.begin(),cbasis.end());
  Nbasis=ubasis[0]+ubasis[1];
    MSize=cbasis[0]+cbasis[1];
   
  for(int I = 0 ; I <= Nbasis-1 ; I++){

     cout << "Alpha["<<I<<"] = "<< Alpha[I] << " d["<<I<<"] = "<< d[I]
     << " l["<<I<<"] = "<< l[I] << " m["<<I<<"] = "<< m[I]
     << " n["<<I<<"] = "<< n[I]<<   endl;
  }

  }else if (atoms==1){
     Nbasis=ubasis[0];
    MSize=cbasis[0];
    
   }
   N=MSize;
 fp.close();
//------END ReadBasisSets----------
}
 
///////////       Smatrix      /////////////////
// 23 May 2003    This version returns "s11" and "s12" for the
//  Gaussians centered on different nulcei.
//  The values of l1,m1,n1,l2,m2,n2 must be specified by the calling function
// and aren't automatically identified by the indices "I" and "J".
void Smatrix(Mat_O_DP &SMatr,vector<double> &Alpha,vector<int> &l,vector<int> &m,vector<int> &n){
  double A,B,Factor;
  double rAx,rAy,rAz,rBx,rBy,rBz,BAx,BAy,BAz,AB2;
  int l1,m1,n1,l2,m2,n2;
  for(int I = 0 ; I <= Nbasis-1 ; I++){
     for(int J = 0 ; J <= I; J++){
  A = Alpha[I];
  B = Alpha[J];
   rAx =  (2*(I/ubasis[0])-1)*dx/2;
   rAy =  (2*(I/ubasis[0])-1)*dy/2;
   rAz =  (2*(I/ubasis[0])-1)*dz/2;
   rBx =  (2*(J/ubasis[0])-1)*dx/2;
   rBy =  (2*(J/ubasis[0])-1)*dy/2;
   rBz =  (2*(J/ubasis[0])-1)*dz/2;
   BAx = rBx -  rAx;
   BAy = rBy -  rAy;
   BAz = rBz -  rAz;
   AB2= BAx*BAx + BAy*BAy + BAz*BAz;
  double gamma = Alpha[I] +  Alpha[J];
  double PAx =  Alpha[J]*(BAx)/gamma;
  double PBx =  -Alpha[I]*(BAx)/gamma;
  double PAy =  Alpha[J]*(BAy)/gamma;
  double PBy =  -Alpha[I]*(BAy)/gamma;
  double PAz =  Alpha[J]*(BAz)/gamma;
  double PBz =  -Alpha[I]*(BAz)/gamma;
   l1 = l[I];
   m1 = m[I];
   n1 = n[I];
   l2 = l[J];
   m2 = m[J];
   n2 = n[J];

  double Sigmax12 = 0.0;
  double Sigmay12 = 0.0;
  double Sigmaz12 = 0.0;

  Factor = PI/(A+B);

  Factor = Factor*sqrt(Factor);
  for(int i = 0 ; i <= floor((l1 + l2)/2.0) ; i++){
     if (2*i >=1 ){

     Sigmax12 = Sigmax12 + f_(2*i,l1,l2,PAx,PBx)
                        *gsl_sf_doublefact(2*i -1)/pow(2.0*gamma,i);
       }else{
               Sigmax12 = Sigmax12 + f_(2*i,l1,l2,PAx,PBx);}

}

    for(int i = 0 ; i <= floor((m1 + m2)/2.0) ; i++){
      if (2*i >=1 ){

     Sigmay12 = Sigmay12 + f_(2*i,m1,m2,PAy,PBy)
                        *gsl_sf_doublefact(2*i -1)/pow(2.0*gamma,i);
        }else{
               Sigmay12 = Sigmay12 + f_(2*i,m1,m2,PAy,PBy);}

}

    for(int i = 0 ; i <= floor((n1 + n2)/2.0) ; i++){
      if (2*i >=1 ){
     Sigmaz12 = Sigmaz12 + f_(2*i,n1,n2,PAz,PBz)
                        *gsl_sf_doublefact(2*i -1)/pow(2.0*gamma,i);
       }else{
               Sigmaz12 = Sigmaz12 + f_(2*i,n1,n2,PAz,PBz);}

}

          SMatr[I][J] = Factor*exp(-A*B/(A+B)*AB2)*Sigmax12*Sigmay12*Sigmaz12;
        SMatr[J][I] = SMatr[I][J];

 }}}
////////////////////////////////////////////



double bico(int n, int k){
  return floor(0.5 + exp(gsl_sf_lnfact(n) -gsl_sf_lnfact(k) - gsl_sf_lnfact(n-k)));

}

double f_(int k,int l1,int l2,double PA,double PB){
  double f=0;
    for(int i = 0 ; i <= l1 ; i++){
      for(int j = 0 ; j <= l2 ; j++){
       if (i+j==k)  f =f + pow(PA,l1-i)*bico(l1, i)*pow(PB,l2-j)*bico(l2, j);
}}
  return f;
}


      double F(int m, double X){


        double f=0.0;
     const double FFac = 0.5*sqrt(PI);


      if (X == 0.0) f = 1.0/(2.0*m + 1.0 );
      else if (m==0) f = FFac/sqrt(X)*gsl_sf_erf(sqrt(X));
      else if (m==1){      f = FFac/sqrt(X)*gsl_sf_erf(sqrt(X));
                           f = ( f - exp(-X) )/2.0/X;
     }else if (m==2){   f = FFac/sqrt(X)*gsl_sf_erf(sqrt(X));
                        f = ( f - exp(-X) )/2.0/X;
                        f = ( 3.0*f  - exp(-X) )/2.0/X;
     }else if (m==3){   f = exp(-X)*(-30.0*sqrt(X) -20.0*pow(X,1.5) -8.0*pow(X,2.5) +15.0*exp(X)*sqrt(PI)*gsl_sf_erf(sqrt(X)) )/16.0/pow(X,3.5);
    }else if (m==4){   f = exp(-X)*(-210.0*sqrt(X) -140.0*pow(X,1.5) -56.0*pow(X,2.5) -16.0*pow(X,3.5) +105.0*exp(X)*sqrt(PI)*gsl_sf_erf(sqrt(X)) )/32.0/pow(X,4.5);
      }else cout << " SOMETHINGS ROTTEN IN DENMARK!!! m = " << m << endl;
      return f;
 }

/////////////////////////////////////////////////////////////

void S_contract(Mat_O_DP &SMatrix,Mat_I_DP &SMatr,vector<double> &Alpha,vector<double> &d,vector<int> &l,vector<int> &m,vector<int> &n,vector<int> &Ncontract){
//--------------------------------------------
       int ilow,iup,jlow,jup;
       int icontract=0;
       int jcontract;
      for(int I = 0 ; I <= MSize-1 ; I++){
        ilow=icontract;
        icontract+= Ncontract[I];
         iup=icontract-1;
        jcontract=0;
        for(int J = 0 ; J <= MSize-1 ; J++){
        jlow=jcontract;
        jcontract+= Ncontract[J];
         jup=jcontract-1;
       
        SMatrix[I][J] =0;

      for(int i = ilow ; i <= iup; i++){
        for(int j = jlow ; j <= jup; j++){

double factor = 1;
if  (2*l[i] >= 1) factor = gsl_sf_doublefact(2*l[i] -1);
else if (2*m[i] >= 1) factor = factor*gsl_sf_doublefact(2*m[i] -1);
else if (2*n[i] >= 1) factor = factor*gsl_sf_doublefact(2*n[i] -1);
else if (2*l[j] >= 1) factor = factor*gsl_sf_doublefact(2*l[j] -1);
else if (2*m[j] >= 1) factor = factor*gsl_sf_doublefact(2*m[j] -1);
else if (2*n[j] >= 1) factor = factor*gsl_sf_doublefact(2*n[j] -1);


SMatrix[I][J] = SMatrix[I][J] +   SMatr[i][j]
*d[i]*d[j]*pow(2.0/PI,1.50)*pow(Alpha[i]*Alpha[j],.750)
*pow(2.0,l[i]+m[i]+n[i]+l[j]+m[j]+n[j] )*pow(Alpha[i],0.50*(l[i]+m[i]+n[i]))*pow(Alpha[j],0.50*(l[j]+m[j]+n[j]))/pow(factor, 0.50);

}}

}}
//------------------------------------------------------
}

////////////////////


void KE(Mat_O_DP &Kinet2,vector<double> &Alpha,vector<int> &l,vector<int> &m,vector<int> &n){
  double Ix11,Iy11,Iz11;//Ix12,Iy12,Iz12,s11,s12;


//      double Kinet12 = 0;

      for(int I = 0 ; I <= Nbasis-1; I++){
        for(int J = 0 ; J <= I; J++){

  Ix11 =l[I]*l[J]*dSMatrix(Alpha, I,  J, l[I]-1,m[I],n[I],l[J]-1,m[J],n[J]);
  Ix11= Ix11 + 4.0*Alpha[I]*Alpha[J]*dSMatrix(Alpha, I,  J, l[I]+1,m[I],n[I],l[J]+1,m[J],n[J]);
 Ix11= Ix11  -2.0*Alpha[I]*l[J]*dSMatrix(Alpha, I,  J, l[I]+1,m[I],n[I],l[J]-1,m[J],n[J]);
  Ix11= 0.5*(Ix11 - 2.0*l[I]*Alpha[J]*dSMatrix(Alpha, I,  J, l[I]-1,m[I],n[I],l[J]+1,m[J],n[J]));


////////////////////////////////////////////////////////

  Iy11 =m[I]*m[J]*dSMatrix(Alpha, I,  J, l[I],m[I]-1,n[I],l[J],m[J]-1,n[J]);
  Iy11= Iy11 + 4.0*Alpha[I]*Alpha[J]*dSMatrix(Alpha, I,  J, l[I],m[I]+1,n[I],l[J],m[J]+1,n[J]);
 Iy11= Iy11  -2.0*Alpha[I]*m[J]*dSMatrix(Alpha, I,  J, l[I],m[I]+1,n[I],l[J],m[J]-1,n[J]);
  Iy11= 0.5*(Iy11 - 2.0*m[I]*Alpha[J]*dSMatrix(Alpha, I,  J, l[I],m[I]-1,n[I],l[J],m[J]+1,n[J]));


//////////////////////////////////////////////////////////////

  Iz11 =n[I]*n[J]*dSMatrix(Alpha, I,  J, l[I],m[I],n[I]-1,l[J],m[J],n[J]-1);
  Iz11= Iz11 + 4.0*Alpha[I]*Alpha[J]*dSMatrix(Alpha, I,  J, l[I],m[I],n[I]+1,l[J],m[J],n[J]+1);
 Iz11= Iz11  -2.0*Alpha[I]*n[J]*dSMatrix(Alpha, I,  J, l[I],m[I],n[I]+1,l[J],m[J],n[J]-1);
  Iz11= 0.5*(Iz11 - 2.0*n[I]*Alpha[J]*dSMatrix(Alpha, I,  J, l[I],m[I],n[I]-1,l[J],m[J],n[J]+1));


  Kinet2[I][J] = (Ix11 + Iy11 + Iz11);
  Kinet2[J][I] = Kinet2[I][J];

}}}

//-------------------------------------------


/////////////////////////////////////////////////////////////////////////

// 23 May 2003    This version returns "s11" and "s12" for the
//  Gaussians centered on different nulcei.
//  The values of l1,m1,n1,l2,m2,n2 must be specified by the calling function
// and aren't automatically identified by the indices "I" and "J".
double dSMatrix(vector<double> &Alpha,int I, int J, int l1,int m1,int n1,int l2,int m2,int n2){
  double A,B,Factor;
  double rAx,rAy,rAz,rBx,rBy,rBz,BAx,BAy,BAz,AB2;

  A = Alpha[I];
  B = Alpha[J];
   rAx =  (2*(I/ubasis[0])-1)*dx/2;
   rAy =  (2*(I/ubasis[0])-1)*dy/2;
   rAz =  (2*(I/ubasis[0])-1)*dz/2;
   rBx =  (2*(J/ubasis[0])-1)*dx/2;
   rBy =  (2*(J/ubasis[0])-1)*dy/2;
   rBz =  (2*(J/ubasis[0])-1)*dz/2;
   BAx = rBx -  rAx;
   BAy = rBy -  rAy;
   BAz = rBz -  rAz;
   AB2= BAx*BAx + BAy*BAy + BAz*BAz;
  double gamma = Alpha[I] +  Alpha[J];
  double PAx =  Alpha[J]*(BAx)/gamma;
  double PBx =  -Alpha[I]*(BAx)/gamma;
  double PAy =  Alpha[J]*(BAy)/gamma;
  double PBy =  -Alpha[I]*(BAy)/gamma;
  double PAz =  Alpha[J]*(BAz)/gamma;
  double PBz =  -Alpha[I]*(BAz)/gamma;

  double Sigmax12 = 0.0;
  double Sigmay12 = 0.0;
  double Sigmaz12 = 0.0;

  Factor = PI/(A+B);
  Factor = Factor*sqrt(Factor);
  for(int i = 0 ; i <= floor((l1 + l2)/2.0) ; i++){
     if (2*i >=1 ){

     Sigmax12 = Sigmax12 + f_(2*i,l1,l2,PAx,PBx)
                        *gsl_sf_doublefact(2*i -1)/pow(2.0*gamma,i);
       }else{
               Sigmax12 = Sigmax12 + f_(2*i,l1,l2,PAx,PBx);}
}

    for(int i = 0 ; i <= floor((m1 + m2)/2.0) ; i++){
      if (2*i >=1 ){
     Sigmay12 = Sigmay12 + f_(2*i,m1,m2,PAy,PBy)
                        *gsl_sf_doublefact(2*i -1)/pow(2.0*gamma,i);
        }else{
               Sigmay12 = Sigmay12 + f_(2*i,m1,m2,PAy,PBy);}
}

    for(int i = 0 ; i <= floor((n1 + n2)/2.0) ; i++){
      if (2*i >=1 ){
     Sigmaz12 = Sigmaz12 + f_(2*i,n1,n2,PAz,PBz)
                        *gsl_sf_doublefact(2*i -1)/pow(2.0*gamma,i);
       }else{
               Sigmaz12 = Sigmaz12 + f_(2*i,n1,n2,PAz,PBz);}
}

return Factor*exp(-A*B/(A+B)*AB2)*Sigmax12*Sigmay12*Sigmaz12;

}

//-------------------------------------------------


      void e_nucl(Mat_O_DP &Coul2,vector<double> &Alpha,vector<int> &l,vector<int> &m,vector<int> &n){
/* C Coulomb matrix element.
*/
 //     double K0, K1, t1, t2, t;
  double Factor;
  double Px,Py,Pz;//CP11x,CP11y,CP11z,CP211;

  double rAx,rAy,rAz,rBx,rBy,rBz,rCx,rCy,rCz;//rDx,rDy,rDz;
//  double QCx,QDx,QCy,QDy,QCz,QDz;
//  double AP2,AB2,BP2,CD2;
  double A,B,BAx,BAy,BAz;

//      int ilow,iup,jlow,jup;

      for(int I = 0 ; I <= Nbasis-1; I++){
        for(int J = 0 ; J <= I; J++){
         Coul2[I][J] = 0;
  A = Alpha[I];
  B = Alpha[J];

   rAx =  (2*(I/ubasis[0])-1)*dx/2;            
   rAy =  (2*(I/ubasis[0])-1)*dy/2;             
   rAz =  (2*(I/ubasis[0])-1)*dz/2;               
   rBx =  (2*(J/ubasis[0])-1)*dx/2;            
   rBy =  (2*(J/ubasis[0])-1)*dy/2;             
   rBz =  (2*(J/ubasis[0])-1)*dz/2;             
   BAx =rBx -  rAx;                       
   BAy =rBy -  rAy;                          
   BAz =rBz -  rAz;                           
  Px = (rAx*A +  rBx*B)/(A+B);              
  Py = (rAy*A +  rBy*B)/(A+B);              
  Pz = (rAz*A +  rBz*B)/(A+B);               


  double gamma = Alpha[I] +  Alpha[J];
  double eps = 0.25/gamma;
  double PAx =  Alpha[J]*(BAx)/gamma;
  double PBx =  -Alpha[I]*(BAx)/gamma;
  double PAy =  Alpha[J]*(BAy)/gamma;
  double PBy =  -Alpha[I]*(BAy)/gamma;
  double PAz =  Alpha[J]*(BAz)/gamma;
  double PBz =  -Alpha[I]*(BAz)/gamma;

  double CPx,CPy,CPz,CP2;
//  double AP2= PAx*PAx + PAy*PAy + PAz*PAz;
  double AB2= BAx*BAx + BAy*BAy + BAz*BAz;
//  double BP2= PBx*PBx + PBy*PBy + PBz*PBz;

//  double Aiur11 = 0.0;
//  double Ajsv11 = 0.0;
//  double Aktw11 = 0.0;
//  double Nucl11 = 0.0;
  double Aiur12 = 0.0;
  double Ajsv12 = 0.0;
  double Aktw12 = 0.0;
  double Nucl12 = 0.0;
//  double s;
  int l1 = l[I];
  int l2 = l[J];
  int m1 = m[I];
  int m2 = m[J];
  int n1 = n[I];
  int n2 = n[J];
  Factor = PI/(A+B);
  Factor   = Factor*exp(-A*B/(A+B)*AB2);


for(int nuc = 0 ; nuc <= atoms-1 ; nuc++){

   rCx =  (2.0*nuc-1.0)*dx/2;            
   rCy =  (2.0*nuc-1.0)*dy/2;              
   rCz =  (2.0*nuc-1.0)*dz/2;
   CPx =  rCx - Px ;
   CPy =  rCy - Py ;
   CPz =  rCz - Pz ;
   CP2 =  CPx*CPx  +   CPy*CPy  +    CPz*CPz;

    for(int i = 0 ; i <= l1 + l2 ; i++){
      for(int r = 0 ; r <= floor(i/2.0) ; r++){
         for(int u = 0 ; u <= floor((i - 2*r)/2.0) ; u++){

    Aiur12 =  pow(-1.0,u)*gsl_sf_fact(i)*f_(i,l1,l2,PAx,PBx)*pow(CPx,i-2*(r +u))
    *pow(eps,r+u)/gsl_sf_fact(r)/gsl_sf_fact(u)/gsl_sf_fact(i-2*r-2*u);

    for(int j = 0 ; j <= m1 + m2 ; j++){
      for(int s = 0 ; s <= floor(j/2.0) ; s++){
         for(int v = 0 ; v <= floor((j - 2*s)/2.0) ; v++){

    Ajsv12 =  pow(-1.0,v)*gsl_sf_fact(j)*f_(j,m1,m2,PAy,PBy)*pow(-CPy,j-2*(s +v))
    *pow(eps,s+v)/gsl_sf_fact(s)/gsl_sf_fact(v)/gsl_sf_fact(j-2*s-2*v);

    for(int k = 0 ; k <= n1 + n2 ; k++){
      for(int t = 0 ; t <= floor(k/2.0) ; t++){
         for(int w = 0 ; w <= floor((k - 2*t)/2.0) ; w++){

    Aktw12 =  pow(-1.0,w)*gsl_sf_fact(k)*f_(k,n1,n2,PAz,PBz)*pow(CPz,k-2*(t +w))
    *pow(eps,t+w)/gsl_sf_fact(t)/gsl_sf_fact(w)/gsl_sf_fact(k-2*t-2*w);


    int mu = i+j+k -2*(r + s + t) - u - v - w;

    Nucl12 = Nucl12 + Aiur12*Ajsv12*Aktw12*Z[nuc]*F(mu, gamma*CP2);

}}} }}} }}}
}

    Coul2[I][J] = -2.0*Factor*Nucl12;
    Coul2[J][I] =Coul2[I][J];
}} // I,J  Loop

}//      END

//------------------------------
//------------------------------       

      void CalcHMatrix(Mat_O_DP &HMatrix,Mat_O_DP &KE_Contract,Mat_O_DP &e_nucl_Contract,Mat_I_DP &Kinet2,Mat_I_DP &Coul2,vector<double> &Alpha,vector<double> &d,vector<int> &l,vector<int> &m,vector<int> &n,vector<int> &Ncontract){
/* C Coulomb matrix element. The subscript 11 indicates that the
C basis functions are centred on the same nucleus, 12 means that
C the kinetic matrix element is calculated for two basis functions
*/

        double Coul12;
        double Kinet12;


       int ilow,iup,jlow,jup;
       int icontract=0;
       int jcontract;
      for(int R = 0 ; R <= MSize-1 ; R++){
        ilow=icontract;
        icontract+= Ncontract[R];
         iup=icontract-1;
        jcontract=0;
         for(int S = 0 ; S <= R; S++){
        jlow=jcontract;
        jcontract+= Ncontract[S];
         jup=jcontract-1;

         Coul12 =  0;
         Kinet12 = 0;

      for(int I = ilow ; I <= iup; I++){
        for(int J = jlow ; J <= jup; J++){


double factor = 1;
if  (2*l[I] >= 1) factor = gsl_sf_doublefact(2*l[I] -1);
else if (2*m[I] >= 1) factor = factor*gsl_sf_doublefact(2*m[I] -1);
else if (2*n[I] >= 1) factor = factor*gsl_sf_doublefact(2*n[I] -1);
else if (2*l[J] >= 1) factor = factor*gsl_sf_doublefact(2*l[J] -1);
else if (2*m[J] >= 1) factor = factor*gsl_sf_doublefact(2*m[J] -1);
else if (2*n[J] >= 1) factor = factor*gsl_sf_doublefact(2*n[J] -1);

      Coul12  = Coul12  +  Coul2[I][J]* d[I]*d[J]*pow(2.0/PI,1.50)*pow(Alpha[I]*Alpha[J],.750)
*pow(2.0,l[I]+m[I]+n[I]+l[J]+m[J]+n[J] )*pow(Alpha[I],0.50*(l[I]+m[I]+n[I]))*pow(Alpha[J],0.50*(l[J]+m[J]+n[J]))/pow(factor, 0.50);

  Kinet12 = Kinet12 + Kinet2[I][J]*
                    d[I]*d[J]*pow(2.0/PI,1.50)*pow(Alpha[I]*Alpha[J],.750)
*pow(2.0,l[I]+m[I]+n[I]+l[J]+m[J]+n[J] )*pow(Alpha[I],0.50*(l[I]+m[I]+n[I]))*pow(Alpha[J],0.50*(l[J]+m[J]+n[J]))/pow(factor, 0.50);
           
    
    }}
    KE_Contract[R][S]= Kinet12;
    e_nucl_Contract[R][S]= Coul12;
    HMatrix[R][S] = Kinet12 +  Coul12;
    }}
   
}

// -----------------------------------------

// -----------------------------------------

     void ERI(Mat4D_O_DP &QMatrix2,vector<double> &Alpha,vector<int> &l,vector<int> &m,vector<int> &n){

  double Factor;//CPx[2],CPy[2],CPz[2],CP2[2],Cx[2],Cy[2],Cz[2];
  double Px,Py,Pz,Qx,Qy,Qz;//CP11x,CP11y,CP11z,CP211;
  double PQx,PQy,PQz,PQ2,P,Q;    
  double DCx,DCy,DCz;
  double BAx,BAy,BAz;
  double rAx,rAy,rAz,rBx,rBy,rBz,rCx,rCy,rCz,rDx,rDy,rDz;
  double PAx,PBx,PAy,PBy,PAz,PBz,QCx,QDx,QCy,QDy,QCz,QDz;
  double AP2,AB2,BP2,CD2;
  double A,B,CC,D;
  int MaxU,L1,L2,M1,M2,N1,N2;//R, S, T, U, 
  int l1,m1,n1,l2,m2,n2,l3,m3,n3,l4,m4,n4;
      for(int i = 0 ; i <= Nbasis -1; i++){
        for(int j = 0 ; j <= i; j++){
        for(int k = 0 ; k <= i; k++){
            if (k < i) MaxU = k;
            else MaxU = j;

      for(int ll = 0 ; ll <= MaxU; ll++){

QMatrix2[i][j][k][ll]= 0;


   
   A = Alpha[i];
   l1 = l[i];
   m1 = m[i];
   n1 = n[i];

   B = Alpha[j];
   l2 = l[j];
   m2 = m[j];
   n2 = n[j];

     CC = Alpha[k];
     l3 = l[k];
     m3 = m[k];
     n3 = n[k];

     D = Alpha[ll];
     l4 = l[ll];
     m4 = m[ll];
     n4 = n[ll];
  

  double gamma1 = A +  B;
  double gamma2 = CC +  D;
  P = ((2*(i/ubasis[0])-1)*A+(2*(j/ubasis[0])-1)*B)/(A+B)*Dist/2;
  Q = ((2*(k/ubasis[0])-1)*CC+(2*(ll/ubasis[0])-1)*D)/(CC+D)*Dist/2;

  PQ2  = (P-Q)*(P-Q);
  double eps1 = 0.25/gamma1;
  double eps2 = 0.25/gamma2;
  double delta = eps1 + eps2;

   rAx =  (2*(i/ubasis[0])-1)*dx/2;
   rAy =  (2*(i/ubasis[0])-1)*dy/2;
   rAz =  (2*(i/ubasis[0])-1)*dz/2;
   rBx =  (2*(j/ubasis[0])-1)*dx/2;
   rBy =  (2*(j/ubasis[0])-1)*dy/2;
   rBz =  (2*(j/ubasis[0])-1)*dz/2;
   rCx =  (2*(k/ubasis[0])-1)*dx/2;
   rCy =  (2*(k/ubasis[0])-1)*dy/2;
   rCz =  (2*(k/ubasis[0])-1)*dz/2;
   rDx =  (2*(ll/ubasis[0])-1)*dx/2;
   rDy =  (2*(ll/ubasis[0])-1)*dy/2;
   rDz =  (2*(ll/ubasis[0])-1)*dz/2;
   BAx =rBx -  rAx;
   BAy =rBy -  rAy;
   BAz =rBz -  rAz;
   DCx = rDx - rCx;
   DCy = rDy - rCy;
   DCz = rDz - rCz;
  Px = (rAx*A +  rBx*B)/(A+B);
  Qx = (rCx*CC + rDx*D)/(CC+D);
  Py = (rAy*A +  rBy*B)/(A+B);
  Qy = (rCy*CC + rDy*D)/(CC+D);
  Pz = (rAz*A +  rBz*B)/(A+B);
  Qz = (rCz*CC + rDz*D)/(CC+D);

  PQx = Px-Qx;
  PQy = Py-Qy;
  PQz = Pz-Qz;

   PAx =  Alpha[j]*(BAx)/gamma1;
   PBx =  -Alpha[i]*(BAx)/gamma1;
   PAy =  Alpha[j]*(BAy)/gamma1;
   PBy =  -Alpha[i]*(BAy)/gamma1;
   PAz =  Alpha[j]*(BAz)/gamma1;
   PBz =  -Alpha[i]*(BAz)/gamma1;

   QCx =  Alpha[ll]*(DCx)/gamma2;
   QDx =  -Alpha[k]*(DCx)/gamma2;
   QCy =  Alpha[ll]*(DCy)/gamma2;
   QDy =  -Alpha[k]*(DCy)/gamma2;
   QCz =  Alpha[ll]*(DCz)/gamma2;
   QDz =  -Alpha[k]*(DCz)/gamma2;

   AP2= PAx*PAx + PAy*PAy + PAz*PAz;
   AB2= BAx*BAx + BAy*BAy + BAz*BAz;
   CD2 = DCx*DCx + DCy*DCy + DCz*DCz;

   BP2= PBx*PBx + PBy*PBy + PBz*PBz;


  double Biur = 0.0;
  double Bjsv = 0.0;
  double Bktw = 0.0;
  double Nucl11 = 0.0;


    Factor = PI/(gamma1+gamma2);
    Factor = 2.0*PI*PI*sqrt(Factor)/gamma1/gamma2*exp(-A*B/(A+B)*AB2 - CC*D/(CC+D)*CD2);

    for(int i1 = 0 ; i1 <= l1 + l2 ; i1++){
      for(int i2 = 0 ; i2 <= l3 + l4 ; i2++){
      for(int r1 = 0 ; r1 <= floor(i1/2.0) ; r1++){
      for(int r2 = 0 ; r2 <= floor(i2/2.0) ; r2++){
     L1 = i1 - 2*r1;
     L2 = i2 - 2*r2;

         for(int u = 0 ; u <= floor((L1 + L2)/2.0) ; u++){

    Biur =  pow(-1.0,i1+u)*gsl_sf_fact(i1)*gsl_sf_fact(i2)*gsl_sf_fact(L1 + L2)*f_(i1,l1,l2,PAx,PBx)*f_(i2,l3,l4,QCx,QDx)*
    pow(PQx,L1 + L2 - 2*u)*pow(eps1,i1 - r1)*pow(eps2,i2 - r2)*pow(delta,u - L1 - L2 )
    /gsl_sf_fact(r1)/gsl_sf_fact(r2)/gsl_sf_fact(L1)/gsl_sf_fact(L2)/gsl_sf_fact(u)/gsl_sf_fact(L1 + L2 - 2*u);


    for(int j1 = 0 ; j1 <= m1 + m2 ; j1++){
      for(int j2 = 0 ; j2 <= m3 + m4 ; j2++){
      for(int s1 = 0 ; s1 <= floor(j1/2.0) ; s1++){
        for(int s2 = 0 ; s2 <= floor(j2/2.0) ; s2++){
         for(int v = 0 ; v <= floor((j1 + j2)/2.0 - s1 -s2); v++){

     M1 = j1 - 2*s1;
     M2 = j2 - 2*s2;
    Bjsv =  pow(-1.0,j1+v)*gsl_sf_fact(j1)*gsl_sf_fact(j2)*gsl_sf_fact(M1 + M2)*f_(j1,m1,m2,PAy,PBy)*f_(j2,m3,m4,QCy,QDy)*
    pow(PQy,M1 + M2 - 2*v)*pow(eps1,j1 - s1)*pow(eps2,j2 - s2)*pow(delta,v - M1 - M2 )
    /gsl_sf_fact(s1)/gsl_sf_fact(s2)/gsl_sf_fact(M1)/gsl_sf_fact(M2)/gsl_sf_fact(v)/gsl_sf_fact(M1 + M2 - 2*v);


    for(int k1 = 0 ; k1 <= n1 + n2 ; k1++){
      for(int k2 = 0 ; k2 <= n3 + n4 ; k2++){
      for(int t1 = 0 ; t1 <= floor(k1/2.0) ; t1++){
        for(int t2 = 0 ; t2 <= floor(k2/2.0) ; t2++){
         for(int w = 0 ; w <= floor((k1 + k2)/2.0 - t1 -t2); w++){

     N1 = k1 - 2*t1;
     N2 = k2 - 2*t2;

    Bktw =  pow(-1.0,k1+w)*gsl_sf_fact(k1)*gsl_sf_fact(k2)*gsl_sf_fact(N1 + N2)*f_(k1,n1,n2,PAz,PBz)*f_(k2,n3,n4,QCz,QDz)*
    pow(PQz,N1 + N2 - 2*w)*pow(eps1,k1 - t1)*pow(eps2,k2 - t2)*pow(delta,w - N1 - N2 )
    /gsl_sf_fact(t1)/gsl_sf_fact(t2)/gsl_sf_fact(N1)/gsl_sf_fact(N2)/gsl_sf_fact(w)/gsl_sf_fact(N1 + N2 - 2*w);


    int mu = i1 + i2 + j1 + j2 + k1 + k2  -2*(r1 + r2 + s1 + s2 + t1 + t2) - u - v - w;

    Nucl11 = Nucl11 + Biur*Bjsv*Bktw*F(mu, 0.25*PQ2/delta);

}}}}}  }}}}}  }}}}}

QMatrix2[i][j][k][ll] = Factor*Nucl11;


      A = QMatrix2[i][j][k][ll];
      QMatrix2[j][i][k][ll] = A;
      QMatrix2[i][j][ll][k] = A;
      QMatrix2[j][i][ll][k] = A;
      QMatrix2[k][ll][i][j] = A;
      QMatrix2[ll][k][i][j] = A;
      QMatrix2[k][ll][j][i] = A;
      QMatrix2[ll][k][j][i] = A;


}}}} // END R,S,T,U Loops
}

/////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////


void ERI_Contract(Mat4D_O_DP &QMatrix,Mat4D_I_DP &QMatrix2,vector<double> &Alpha,vector<double> &d,vector<int> &l,vector<int> &m,vector<int> &n,vector<int> &Ncontract){

      double A;
       int ilow,iup,jlow,jup,klow,kup,llow,lup,MaxU;
       int icontract=0;
       int jcontract,kcontract,lcontract;
      for(int R = 0 ; R <= MSize-1 ; R++){
        ilow=icontract;
        icontract+= Ncontract[R];
         iup=icontract-1;
        jcontract=0;
         for(int S = 0 ; S <= R; S++){
          jlow=jcontract;
          jcontract+= Ncontract[S];
          jup=jcontract-1;
          kcontract=0;
        for(int T = 0 ; T <= R /* !!! */; T++){
          klow=kcontract;
          kcontract+= Ncontract[T];
          kup=kcontract-1;

            if (T < R) MaxU = T;
            else MaxU = S;
            lcontract=0;
      for(int U = 0 ; U <= MaxU; U++){
          llow=lcontract;
          lcontract+= Ncontract[U];
          lup=lcontract-1;

QMatrix[R][S][T][U]= 0;
      for(int i = ilow ; i <= iup; i++){
        for(int j = jlow ; j <= jup; j++){
          for(int k = klow ; k <= kup; k++){
            for(int ll = llow ; ll <= lup; ll++){


double factor = 1;
if  (2*l[i] >= 1) factor = gsl_sf_doublefact(2*l[i] -1);
else if (2*m[i] >= 1) factor = factor*gsl_sf_doublefact(2*m[i] -1);
else if (2*n[i] >= 1) factor = factor*gsl_sf_doublefact(2*n[i] -1);
else if (2*l[j] >= 1) factor = factor*gsl_sf_doublefact(2*l[j] -1);
else if (2*m[j] >= 1) factor = factor*gsl_sf_doublefact(2*m[j] -1);
else if (2*n[j] >= 1) factor = factor*gsl_sf_doublefact(2*n[j] -1);
else if (2*l[k] >= 1) factor = factor*gsl_sf_doublefact(2*l[k] -1);
else if (2*m[k] >= 1) factor = factor*gsl_sf_doublefact(2*m[k] -1);
else if (2*n[k] >= 1) factor = factor*gsl_sf_doublefact(2*n[k] -1);
else if (2*l[ll] >= 1) factor = factor*gsl_sf_doublefact(2*l[ll] -1);
else if (2*m[ll] >= 1) factor = factor*gsl_sf_doublefact(2*m[ll] -1);
else if (2*n[ll] >= 1) factor = factor*gsl_sf_doublefact(2*n[ll] -1);

QMatrix[R][S][T][U] = QMatrix[R][S][T][U] + QMatrix2[i][j][k][ll]
              *d[i]*d[j]*d[k]*d[ll]*pow(2.0/PI,3)*pow(Alpha[i]*Alpha[j]*Alpha[k]*Alpha[ll],.750)
*pow(2.0,l[i]+m[i]+n[i]+l[j]+m[j]+n[j]+l[k]+m[k]+n[k]+l[ll]+m[ll]+n[ll] )
*pow(Alpha[i],0.50*(l[i]+m[i]+n[i]))*pow(Alpha[j],0.50*(l[j]+m[j]+n[j]))
*pow(Alpha[k],0.50*(l[k]+m[k]+n[k]))*pow(Alpha[ll],0.50*(l[ll]+m[ll]+n[ll]))/pow(factor, 0.50);

}}}}


      A = QMatrix[R][S][T][U];
      QMatrix[S][R][T][U] = A;
      QMatrix[R][S][U][T] = A;
      QMatrix[S][R][U][T] = A;
      QMatrix[T][U][R][S] = A;
      QMatrix[U][T][R][S] = A;
      QMatrix[T][U][S][R] = A;
      QMatrix[U][T][S][R] = A;



}}}} // END R,S,T,U Loops

}
        int DetermineType(int R, int S, int T, int U){
//      INTEGER FUNCTION DetermineType(R,S,T,U)
//C Determines relative values of coeffs rstu. See the
//C discussion on page 75 of "Computational Physics".
// Reference: Jos M. Thijssen, "Computational Physics", (Cambridge University Press, 1999).

      int Type;

      Type = 0;
      if ((R == S) && (T == U) && (R==T))  Type = 1;
      else if ((R == S) && (T == U))    Type = 2;
      else if ((R == S) && (T > U)) Type = 3;
      else if ((R > S) && (T == U)) Type = 4;
      else if ((R > S) && (T > U) &&
             (R == T) && (S == U)) Type = 5;
      else Type = 6;

    return Type;
 }


      void BuildG(Mat_O_DP &GMatrix,Mat_I_DP &DensMat,Mat4D_I_DP &QMatrix){
//C Builds the matrix G which is a contraction of the two-electron
//C matrix elements with the density matrix, see Eq. (4.86)
//C Full use is made of symmetries
//Reference: Jos M. Thijssen, "Computational Physics", (Cambridge University Press, 1999).


      double SuperElem;
      int MaxU, Count;//R, S, T, U, 
      cout <<  "     BuildG      " << endl;
      for(int R = 0 ; R <= N-1 ; R++){
        for(int S = 0 ; S <= R ; S++){
            GMatrix[R][S] = 0.0;
}}

      Count = 0;

      for(int R = 0 ; R <= N -1; R++){
        for(int S = 0 ; S <= R ; S++){
            for(int T = 0 ; T <= R ; T++){

            if (T < R)  MaxU = T;
            else MaxU = S;
            for(int U = 0 ; U <= MaxU; U++){

              Count = Count + 1;
              SuperElem = 2.0*QMatrix[R][S][T][U]; // factor of 2 


switch(DetermineType(R,S,T,U))
  {
    case 1:
      GMatrix[R][R] = GMatrix[R][R]+0.5*SuperElem*DensMat[R][R];

      break;
    case 2:
      GMatrix[R][R] = GMatrix[R][R]+0.50*SuperElem*DensMat[T][T];
      GMatrix[T][T] = GMatrix[T][T]+0.50*SuperElem*DensMat[R][R];

      break;
    case 3:
       GMatrix[R][R] = GMatrix[R][R]+SuperElem*DensMat[T][U];
       GMatrix[T][U] = GMatrix[T][U]+0.50*SuperElem*DensMat[R][R];

      break;
    case 4:
       GMatrix[T][T] = GMatrix[T][T]+SuperElem*DensMat[R][S];
       GMatrix[R][S] = GMatrix[R][S]+0.50*SuperElem*DensMat[T][T];

      break;
    case 5:
      GMatrix[R][S] = GMatrix[R][S]+SuperElem*DensMat[R][S];

     break;
    default:
      GMatrix[R][S] = GMatrix[R][S]+SuperElem*DensMat[T][U];
      GMatrix[T][U] = GMatrix[T][U]+SuperElem*DensMat[R][S];


      break;
  }

}}}}

}


      void BuildG4(Mat_O_DP &GMatrix,Mat_O_DP &JMatrix,Mat_O_DP &KMatrix,Vec_I_DP &C,Mat4D_I_DP &QMatrix){

//      int R, S, T, U, MaxU, Count;
      double G[N][N];//error,
      cout <<  "     BuildG4      " << endl;
      for(int R = 0 ; R <= N-1 ; R++){
        for(int S = 0 ; S <= N-1 ; S++){
           GMatrix[R][S] = 0.0;
            JMatrix[R][S] = 0.0;
            KMatrix[R][S] = 0.0;
            G[R][S] = 0.0;
            for(int T = 0 ; T <= N-1 ; T++){
            for(int U = 0 ; U <= N-1; U++){

 //             GMatrix[R][S] = GMatrix[R][S]+2.0*C[T]*C[U]*QMatrix[R][S][T][U];
              JMatrix[R][S] = JMatrix[R][S]+2.0*C[T]*C[U]*QMatrix[R][S][T][U];
              KMatrix[R][S] = KMatrix[R][S]+2.0*C[T]*C[U]*QMatrix[R][U][T][S];
              GMatrix[R][S] = GMatrix[R][S]+C[T]*C[U]*(2.0*QMatrix[R][S][T][U]-QMatrix[R][U][T][S]);
             } // end do U
//              GMatrix[R][S] = GMatrix[R][S] + 0.50*DensMat[T][T]*QMatrix[R][S][T][T];
             }}}  
 }


      void BuildDensMat(Mat_O_DP &DensMat,Vec_I_DP &C,Mat_I_DP &SMatrix){
/*
C Calculates the normalies density matrix D=2 C C^T/(C^T C),
C where C is the eigenvector corresponding to the ground state.

*/
cout << "////// BUILDING THE DENSITY MATRIX WITH THE \"BuildDensMat\" FUNCTION //////"<< endl;
      int R, S;
      double Norm;

      Norm = 0.0;
        for(R = 0 ; R <= N-1 ; R++){
          for(S = 0 ; S <= R-1; S++){
          DensMat[R][S] = 2.0*C[R]*C[S];

Norm = Norm + 2.0*C[R]*C[S]*SMatrix[R][S];
}
        DensMat[R][R] = 2*C[R]*C[R];  
//        Norm = Norm + 0.50*DensMat[R][R]*SMatrix[R][R]; // commented out 23 July 2003
Norm = Norm + C[R]*C[R]*SMatrix[R][R];  // added 23 July 2003
}
cout << "BuildDensMat Norm =  " <<  Norm<< endl;

}


      void CalcFockMatrix(Mat_O_DP &FMatrix,Mat_I_DP &HMatrix,Mat_I_DP &GMatrix){
// Calculate Fock matrix: F = H + G,  F(r,s); r>s
cout << "////// CALCULATING THE FOCK MATRIX WITH THE \"CalcFockMatrix\" FUNCTION //////"<< endl;
//        double Energy;
        for(int R = 0 ; R <= N-1 ; R++){
          for(int S = 0 ; S <= R-1; S++){
          FMatrix[R][S] = HMatrix[R][S] + GMatrix[R][S];
}
          FMatrix[R][R] = HMatrix[R][R] + GMatrix[R][R];
}}

      void DiagFock(Vec_O_DP &C,vector<double> &Diag,Mat_I_DP &VMatrix, Mat_I_DP &FMatrix,Mat_I_DP &SMatrix,int MSize){
cout << "////// DIAGONALIZING THE FOCK MATRIX WITH THE \"DiagFock\" FUNCTION //////"<< endl;
        
      Mat_DP Eigen(MSize,MSize);
      
      LowdinDiag(VMatrix, FMatrix, Eigen, Diag,  MSize);
      for(int I = 0 ; I <= N-1 ; I++) C[I] = Eigen[0][I];

}

void DisplayEnergies(double &Energy,double &energyX,double &energyC,vector<double> &Diag,Mat_I_DP &DensMat,Mat_I_DP &HMatrix,
Mat_I_DP &GMatrix,Mat_I_DP &FMatrix,Mat_I_DP &JMatrix,Mat_I_DP &KMatrix,Vec_I_DP &C,Mat_I_DP &KE_Contract,Mat_I_DP &e_nucl_Contract){
cout << "//////// CALLING THE \"DisplayEnergies\" FUNCTION ////////////"<< endl;


      double Ener;
      double Ecore = 0.0;
      double GEnergy =0.0;
      double JEnergy =0.0;
      double KEnergy =0.0;
      double kinetic =  0.0;
      double coulomb = 0.0;
      double E2electron = 0.0;
      double energyXC = 0.0;
      int I, J;
       for(I = 0 ; I <= MSize-1 ; I++){
cout << " C["<<I<<"]" << " Diag["<<I<<"]  "<< C[I] << "  "  << "  "<< Diag[I] << endl;

}


Energy = 0.0;
Ener = 0.0;
      for(I = 0 ; I <= N-1 ; I++){
        for(J = 0 ; J <= I-1 ; J++){
         Ecore   = Ecore +  2.0*DensMat[I][J]*HMatrix[I][J];
         GEnergy = GEnergy  + DensMat[I][J]*GMatrix[I][J];
         Energy = Energy + DensMat[I][J]*(2.0*HMatrix[I][J]+ GMatrix[I][J]);
         Ener = Ener+2.0*C[I]*C[J]*(FMatrix[I][J]+ HMatrix[I][J]);
         kinetic =  kinetic + 2.0*DensMat[I][J]*KE_Contract[I][J];
         coulomb =  coulomb + 2.0*DensMat[I][J]*e_nucl_Contract[I][J];
         if (method==1){
          E2electron= E2electron + 2.0*C[I]*C[J]*(JMatrix[I][J]-0.5*KMatrix[I][J]);
          JEnergy= JEnergy + 2.0*C[I]*C[J]*JMatrix[I][J];
          KEnergy= KEnergy + C[I]*C[J]*KMatrix[I][J]; }
}//         ENDDO

// 16 July 2003           Energy = Energy  + 0.50*DensMat[I][I]*(FMatrix[I][I] + HMatrix[I][I]);
         Energy = Energy  + 0.50*DensMat[I][I]*(2.0*HMatrix[I][I]+ GMatrix[I][I]);
         Ener =   Ener+C[I]*C[I]*(FMatrix[I][I]+HMatrix[I][I]);
         Ecore =  Ecore + DensMat[I][I]*HMatrix[I][I];
         kinetic =  kinetic + DensMat[I][I]*KE_Contract[I][I];
         coulomb =  coulomb + DensMat[I][I]*e_nucl_Contract[I][I];
         GEnergy = GEnergy  + 0.50*DensMat[I][I]*GMatrix[I][I];
         if (method==1){         
          E2electron= E2electron + C[I]*C[J]*(JMatrix[I][I]-0.5*KMatrix[I][I]);
          JEnergy= JEnergy +     C[I]*C[J]*JMatrix[I][I];
          KEnergy= KEnergy + 0.5*C[I]*C[J]*KMatrix[I][I];  }

}//      ENDDO

      cout <<  "------------------------------------------   " <<   endl;
      cout <<  "The Core Energy is .......   " << Ecore <<  endl;
      
      cout <<  "Energy = Hcore + GMatrix =  " << Ecore + GEnergy << endl;
      cout <<  "Diag[0]+0.5*Ecore  =  " << Diag[0]+0.5*Ecore <<  endl;
      cout <<  "------------------------------------------   " <<   endl;
      cout <<  "ener   DensMat[I][J]*(FMatrix[I][J]+ HMatrix[I][J]) = " << Ener <<  endl;

      cout <<  "Energy verify  = Diag[0]+0.5*Ecore...   " << Diag[0]+0.5*Ecore <<  endl;
      cout <<  "Energy verify2  = 2*Diag[0]-JEnergy...   " << 2.0*Diag[0]-GEnergy <<  endl;
      cout <<  "The Hartree-Fock orbital energy Diag[0] is .......   " << Diag[0] <<  endl;
      cout <<  "The two electron G Energy .......   " << GEnergy <<  endl;
      cout <<  "kinetic = " << kinetic <<  "   ||  electron-nucleous attracton energy = " << coulomb<<  endl;

if (method==1){
      cout <<  "Coulomb repulsion energy = " << JEnergy <<  "  ||    Exchange energy = " << KEnergy<<  endl;
      cout <<  "E2electron= two-electron repulsion energy- exchange energy =...." << E2electron<<  endl;
      cout <<  "Ecore + JEnergy - KEnergy =   " << Ecore + GEnergy <<  endl;
}
      cout <<  "DensMat[I][J]*(2.0*HMatrix[I][J]+ GMatrix[I][J]) = Energy = " << Energy << endl;
if ((method==2)||(method==3)){
      energyXC = energyX  + energyC;
      cout <<  "Diag[0]....................... " << Diag[0]  << endl;
      cout <<  "kinetic....................... " << kinetic  << endl;      
      cout <<  "electron-nucleous............. " << coulomb  << endl;
      cout <<  "Hartree....................... " << GEnergy  << endl;           
      cout <<  "energyX....................... " << energyX  << endl;
      cout <<  "energyC....................... " << energyC  << endl;
      cout <<  "energyXC...................... " << energyXC  << endl;
      cout  << endl;      
      cout <<  "Diag[1]....................... " << Diag[1]  << endl;
      cout <<  "Diag[1]-Diag[0]............... " << Diag[1]-Diag[0]  << endl;
      cout <<  "electron-nucleous + kinetic + Hartree....................... " << coulomb + kinetic + GEnergy<< endl;
      Energy = Energy  + energyXC; // 16 July 2003
      cout <<  "Energy  + energyXC =  " << Energy << endl;
      }
if (atoms==2){      
      cout <<  "The bond distance you input was.................... " <<Dist<<endl
      << "The Nuclear Repulsion Energy = +Z[0]*Z[1]/Dist is .......   " << Z[0]*Z[1]/Dist <<  endl;
      cout <<  "************  THE TOTAL ENERGY  *************   " <<endl;
      cout <<  setprecision(9);
      cout <<  "The TOTAL energy is .......................  " << Energy +Z[0]*Z[1]/Dist << endl;
//      cout <<  " TOTAL energy + 1.3...............  " << Energy +Z[0]*Z[1]/Dist + 1.3<< endl;
      }else{
      cout <<  "************  THE TOTAL ENERGY  *************   " <<endl;
      cout <<  setprecision(9);      
      cout <<  "The TOTAL Energy   =  " << Energy << endl;
 
//       cout << "The TOTAL Energy  +3.0 =  " << Energy +3.0 << endl;
      }
      cout <<  setprecision(6) ;
}


////////////////////////////////////////////////////////////

       void R_Integrate(double &energyX,double &energyC,Mat_IO_DP &FMatrix,Vec_I_DP &C,vector<double> &Alpha,vector<double> &d,vector<int> &l,vector<int> &m,vector<int> &n,vector<int> &Ncontract,int functional){
cout << "//////// STARTING THE DFT 1-D INTEGRATION ROUTINE R_Integrate ////////////"<< endl;

  ofstream fout("Density_Output.txt");  // create output file
  if(!fout) {
    cout << "Cannot open output file.\n";
//    return 1;
  }

       double  ri,w,rint,D; //Ratom,,theta
       double dphi,phi[N],delphix[N],delphiy[N],delphiz[N];
       double fock;//delrhox,delrhoy,delrhoz;//,lyp_wt,vwnrpa_wt,pz81_wt,LDA_wt,b88_wt,pbe_wt;
       double delrhox =0.0;
       double delrhoy =0.0;
       double delrhoz =0.0;
       double energyXC =0.0;
       energyX =0.0;       // Thursday 13 December 2003
       energyC =0.0;       // Thursday 13 December 2003

       double f,dfdra,dfdgaa,dfdgab;   //dfdrb, ,dfdgbb Thursday 10 July 2003
       double fx,dfxdra,dfxdgaa,dfxdgab,fc,dfcdra,dfcdgaa,dfcdgab;   // dfcdgbb dfxdrb dfxdgbb dfcdrb Thursday 13 December 2003
       double int_Vxc_rho =0.0;   // 18 July 2003

//       double dpsi=0.0;   // 18 June 2003 another test (r2,Dist2,sum2)
       double psi=0.0;   // 18 June 2003 another test (r2,Dist2,sum2)
       double rho=0.0;   // Thursday 10 July 2003

//       double ddelpsix=0.0;   // 18 June 2003
       double delpsix=0.0;    // 18 June 2003
//       double ddelpsiy=0.0;   // 18 June 2003
       double delpsiy=0.0;    // 18 June 2003           int_gamma
//       double ddelpsiz=0.0;   // 18 June 2003
       double delpsiz=0.0;    // 18 June 2003

       double int_gamma=0.0;   // Thursday 10 July 2003
       double gamma=0.0;   // Thursday 10 July 2003
       double Electrons=0.0;   // Thursday 10 July 2003
       double slater=0.0;   // Wednesday 16 July 2003

//       int Ns =195;
//       double x[Ns],y[Ns],z[Ns],W[Ns];
       Mat_DP Fxc(0.0,MSize,MSize);
        Vec_INT l_1(0,MSize),m_1(0,MSize),n_1(0,MSize);

       int nr;
       nr = 50;
//       Ratom = 1.0;    // Hydrogen
//       Ratom = 0.5882;    // Helium
cout << " R_Integrate functional = " << functional << endl;
       for(double i = 1 ; i <= nr; i++){
          w=2*pow(Ratom,3)*(nr +1)*pow(i,5)/pow(nr+1-i,7);
          ri =Ratom*pow(i/(nr +1-i),2  );
        
for(int iatom = 0 ; iatom <= atoms-1 ; iatom++){
            D=0.0;
          psi=0.0;
          delpsix =0.0;
          delpsiy =0.0;
          delpsiz =0.0;
       int ilow,iup;//,jlow,jup;
       int icontract=0;
 //      int jcontract;
      for(int I = 0 ; I <= MSize-1 ; I++){
        ilow=icontract;
        icontract+= Ncontract[I];
         iup=icontract-1;

//---------
          rint =ri;
          phi[I] = 0.0;
          delphix[I] = 0.0;
          delphiy[I] = 0.0;
          delphiz[I] = 0.0;

            for(int k = ilow ; k <= iup; k++){

double factor = 1;
if  (2*l[k] >= 1) factor = gsl_sf_doublefact(2*l[k] -1);
else if (2*m[k] >= 1) factor = factor*gsl_sf_doublefact(2*m[k] -1);
else if (2*n[k] >= 1) factor = factor*gsl_sf_doublefact(2*n[k] -1);


        dphi = d[k]*pow(2.0*Alpha[k]/PI,.750)
*pow(2.0,l[k]+m[k]+n[k] )*pow(Alpha[k],0.50*(l[k]+m[k]+n[k]))/pow(factor, 0.50)
      *exp(-pow(rint,2)*Alpha[k]);


phi[I] = phi[I]+dphi;
delphix[I] = delphix[I] - 2.0*Alpha[k]*ri*dphi;
delphiy[I] = 0.0;
delphiz[I] = 0.0;


} //////////////// END DO k =  Contraction over Nbasis ///////////////////

          psi = psi +  C[I]*phi[I];
          delpsix = delpsix + C[I]*delphix[I];
//***           delpsiy = delpsiy + C[I]*delphiy[I];
//***           delpsiz = delpsiz + C[I]*delphiz[I];
} ///////// END DO I =  Sum over basis states N   //////////////////

          D = pow(2*psi*psi,3.0/3.0);
          rho = D;
          delrhox = 4*psi*delpsix;
//*** Comment these out on 28 July 2003
//***          delrhoy = 4*psi*delpsiy;
//***          delrhoz = 4*psi*delpsiz;

          Electrons = Electrons  + w*rho;
//          slater = slater -   W[j]*wA*w*(4*PI*0.930525736349100013861335177933319*pow(rho,4.0/3.0) );
//          slater = slater -   W[j]*wA*w*4*PI*0.75*pow(3.0/PI,1.0/3.0)*pow(rho,4.0/3.0);
          slater = slater +  w*4*PI*0.25*pow(3.0/PI,1.0/3.0)*pow(rho,4.0/3.0);
          gamma = rho*8.0*(delpsix*delpsix + delpsiy*delpsiy + delpsiz*delpsiz);
          int_gamma = int_gamma  + w*gamma;
          f=0.0;
          dfdra=0.0;
          dfdgaa=0.0;
          dfdgab=0.0;
          fx=0.0;
          dfxdra=0.0;
          dfxdgaa=0.0;
          dfxdgab=0.0;
          fc=0.0;
          dfcdra=0.0;
          dfcdgaa=0.0;
          dfcdgab=0.0;
          
          if (functional==1){
          x_rks_b88(rho,gamma,fx,dfxdra,dfxdgaa);
          c_rks_lyp(rho,gamma,fc,dfcdra,dfcdgaa,dfdgab);

          
          }else if (functional==2){
          x_rks_slater(rho,fx,dfxdra);

          }else if (functional==3){
          x_rks_b88(rho,gamma,fx,dfxdra,dfxdgaa);        
          
          }else if (functional==4){
          x_rks_b88(rho,gamma,fx,dfxdra,dfxdgaa);  
          c_rks_vwnrpa(rho,fc,dfcdra);

/*          }else if (functional==4){
          lyp_wt   =  0.81;
          vwnrpa_wt=  0.19;            
          x_rks_b88(rho,gamma,fx,dfxdra,dfxdgaa);
          
          c_rks_lyp(rho,gamma,f,dfdra,dfdgaa,dfdgab);
          fc =lyp_wt*f;
          dfcdra =lyp_wt*dfdra;
          dfcdgaa =lyp_wt*dfdgaa;
          dfcdgab =lyp_wt*dfdgab;
                             
          c_rks_vwnrpa(rho,f,dfdra);                    
         fc+=vwnrpa_wt*f;
          dfcdra+=vwnrpa_wt*dfdra;
          dfcdgaa+=vwnrpa_wt*dfdgaa;
          dfcdgab+=vwnrpa_wt*dfdgab;                     
*/
          }else if (functional==5){
           x_rks_b88(rho,gamma,fx,dfxdra,dfxdgaa);
           c_rks_pz81(rho,fc,dfcdra); 

          
          }else if (functional==6){
            x_rks_pbe(rho,gamma,fx,dfxdra,dfxdgaa);
            c_rks_pz81(rho,fc,dfcdra);

          }

          
          if ((rho==0.0)){
          f=0.0;
          dfdra=0.0;
          dfdgaa=0.0;
          dfdgab=0.0;
          fx=0.0;
          dfxdra=0.0;
          dfxdgaa=0.0;
          dfxdgab=0.0;
          fc=0.0;
          dfcdra=0.0;
          dfcdgaa=0.0;
          dfcdgab=0.0;
}
          energyX+= w*fx;
          energyC+= w*fc;
          energyXC= energyX + energyC;
          f = fx + fc;
          dfdra = dfxdra + dfcdra;
          dfdgaa = dfxdgaa + dfcdgaa;
          dfdgab = dfxdgab + dfcdgab;

          int_Vxc_rho =  int_Vxc_rho + w*rho*dfdra;

          for(int R = 0 ; R <= N-1 ; R++){
             for(int S = 0 ; S <= R /* !!! */; S++){
               fock = dfdra*phi[R]*phi[S]  +  ( dfdgaa + 0.5*dfdgab )*
               (delrhox*(phi[S]*delphix[R] + phi[R]*delphix[S]) +
                delrhoy*(phi[S]*delphiy[R] + phi[R]*delphiy[S]) +
                delrhoz*(phi[S]*delphiz[R] + phi[R]*delphiz[S]) );

               Fxc[R][S] +=  w*fock;

          }}

} ///////////// END DO for Atomic Loop here   //////////////
} // END DO i =  radial integration

// 21 July 2003 } // END DO j =   Angular integration ////////////////////
       cout <<" Ratom =  "<< Ratom <<"  nr =  "<< nr  <<"  ri =  "<< ri << endl;
       cout << "  # of Electrons =  "<<4*PI*Electrons<< endl;
       cout <<"  int_gamma =  " << 4*PI*int_gamma  << endl;

       for(int I = 0 ; I <= N-1 ; I++){
cout << "  C["<<I<<"]   " <<  C[I] << "  Alpha["<<I<<"]   " <<  Alpha[I]  << "  d["<<I<<"]   " <<  d[I]<<  endl;
}

      for(int I = 0 ; I <= N-1 ; I++){
        for(int J = 0 ; J <= I ; J++){
FMatrix[I][J] = FMatrix[I][J] + 4.0*PI*Fxc[I][J];

}}  /**/

      energyXC = 4.0*PI*energyXC;
      energyX = 4.0*PI*energyX;
      energyC = 4.0*PI*energyC;
cout << " energyXC =   "<< energyXC << endl;
cout << " energyX =   "<< energyX << " energyC =   "<< energyC << endl;
       }


////////////////////////////////////////////////////////////



void xyz_Integrate(double &energyX,double &energyC,Mat_IO_DP &FMatrix,Vec_I_DP &C,vector<double> &Alpha,vector<double> &d,vector<int> &l,vector<int> &m,vector<int> &n,vector<int> &Ncontract,int functional){
cout << "//////// STARTING THE DFT 3-D INTEGRATION ROUTINE xyz_Integrate ////////////"<< endl;
  ofstream fout("Density_Output.txt");  // create output file
  if(!fout) {
    cout << "Cannot open output file.\n";
//    return 1;
  }

       double  w,sum,s1,s2,wA,wB,mu,D; //Ratom,,theta
       double dphi;
       double fock,delrhox,delrhoy,delrhoz;//lyp_wt,vwnrpa_wt,pz81_wt,LDA_wt,b88_wt,pbe_wt;
//       double sum1 =0.0;
	double rint =0.0;
       double energyXC =0.0;
       energyX =0.0;       // Thursday 13 December 2003
       energyC =0.0;       // Thursday 13 December 2003

       double f,dfdra,dfdgaa,dfdgab;   // dfdrb,dfdgbb  Thursday 10 July 2003
       double fx,dfxdra,dfxdgaa,dfxdgab,fc,dfcdra,dfcdgaa,dfcdgab;   // Thursday 13 December 2003
//       double dfxdrb,dfxdgbb,dfcdrb,dfcdgbb;   // Thursday 20 January 2005       
       double int_Vxc_rho =0.0;   // 18 July 2003
//       double sum2 =0.0;  // 11 June 2003 another test (r2,Dist2,sum2)
       double Dist2 =Dist;  // 11 June 2003 another test (r2,Dist2,sum2)

       double ri=0.0;
       double rB=0.0;	      
       double rA=0.0;
       double r2=0.0;   // 11 June 2003 another test (r2,Dist2,sum2)
//       double dpsi=0.0;   // 18 June 2003 another test (r2,Dist2,sum2)
       double psi=0.0;   // 18 June 2003 another test (r2,Dist2,sum2)
       double rho=0.0;   // Thursday 10 July 2003
//       double phi[N]=0.0;   // Thursday 13 July 2003
//       double ddelpsix=0.0;   // 18 June 2003
       double delpsix=0.0;    // 18 June 2003
//       double ddelpsiy=0.0;   // 18 June 2003
       double delpsiy=0.0;    // 18 June 2003           int_gamma
//       double ddelpsiz=0.0;   // 18 June 2003
       double delpsiz=0.0;    // 18 June 2003
       double int_gamma=0.0;   // Thursday 10 July 2003
       double gamma=0.0;   // Thursday 10 July 2003
       double Electrons=0.0;   // Thursday 10 July 2003


       int Ns =195;
       Vec_DP x(0.0,Ns),y(0.0,Ns),z(0.0,Ns),W(0.0,Ns);
       Vec_DP phi(0.0,MSize),delphix(0.0,MSize),delphiy(0.0,MSize),delphiz(0.0,MSize);
       Mat_DP Fxc(0.0,MSize,MSize);
        Vec_INT l_1(0,Nbasis),m_1(0,Nbasis),n_1(0,Nbasis);

       int nr;//ilow,iup;
       nr = 50;
       sum = 0;

       LD0194(&x[0],&y[0],&z[0],&W[0]);
//       LD0230(x,y,z,W);
cout << "Ratom = "<< Ratom << " atoms = "<< atoms<< endl;
       for(int j = 1 ; j <= Ns-1 ; j++){
       for(double i = 1 ; i <= nr; i++){
          w=2*pow(Ratom,3)*(nr +1)*pow(i,5)/pow(nr+1-i,7);
          ri =Ratom*pow(i/(nr +1-i),2  );
// 15 July 2003          rB = sqrt(Dist*Dist+ ri*ri -2*ri*Dist*z[j] );
          r2 = sqrt(Dist2*Dist2+ ri*ri -2*ri*Dist2*z[j] );

for(int iatom = 0 ; iatom <= atoms-1 ; iatom++){
            D=0.0;
          psi=0.0;
          delpsix =0.0;
          delpsiy =0.0;
          delpsiz =0.0;
        rB = sqrt(Dist*Dist+ ri*ri -2*ri*pow(-1.0,iatom)*Dist*z[j] );
/////////////////
       int ilow,iup;//jlow,jup;
       int icontract=0;
//cout << "i= " <<  i << " j= " <<  j << " iatom= " <<  iatom << " z[j]= " <<  z[j] << " ri= " <<  ri << " rA= " <<  rA << "  ri*z[j]-rA= " <<  ri*z[j]-rA <<   endl;
      for(int I = 0 ; I <= MSize-1 ; I++){
        ilow=icontract;
        icontract+= Ncontract[I];
         iup=icontract-1;
/////////////////                 
if (I==0){
         rint = ri;
         rA=-Dist*iatom;
         rint = ri*(1-iatom) + rB*iatom;

}else if (I>=cbasis[iatom]){
           rA=Dist*(1-iatom);
           rint = rB*(1-iatom) + ri*iatom;

        }
/**/
          phi[I] = 0.0;
          delphix[I] = 0.0;
          delphiy[I] = 0.0;
          delphiz[I] = 0.0;

            for(int k = ilow ; k <= iup; k++){

double factor = 1;
if  (2*l[k] >= 1) factor = gsl_sf_doublefact(2*l[k] -1);
else if (2*m[k] >= 1) factor = factor*gsl_sf_doublefact(2*m[k] -1);
else if (2*n[k] >= 1) factor = factor*gsl_sf_doublefact(2*n[k] -1);

/* Maybe the Integration loop should be moved here. */
// 14 July 2003 might need to make phi and delphi's arrays
        dphi = d[k]*pow(2.0*Alpha[k]/PI,.750)
*pow(2.0,l[k]+m[k]+n[k] )*pow(Alpha[k],0.50*(l[k]+m[k]+n[k]))/pow(factor, 0.50)
      *exp(-pow(rint,2)*Alpha[k]);


//delphix[I] += ( l[k]*pow(ri*x[j],l_1[k]) -  2.0*Alpha[k]*ri*x[j]*pow(ri*x[j],l[k]) )*dphi*pow(ri*y[j],m[k])*pow(ri*z[j]-rA,n[k]);
//delphiy[I] += ( m[k]*pow(ri*y[j],m_1[k]) -     2.0*Alpha[k]* ri*y[j]*pow(ri*y[j],m[k]) )*dphi*pow(ri*x[j],l[k])*pow(ri*z[j]-rA,n[k]);
//delphiz[I] += ( n[k]*pow(ri*z[j]-rA,n_1[k]) -  2.0*Alpha[k]*(ri*z[j]-rA)*pow(ri*z[j]-rA,n[k]) )*dphi*pow(ri*x[j],l[k])*pow(ri*y[j],m[k]);

delphix[I] += ( l[k]*pow(ri*x[j],l_1[k]) -     2.0*Alpha[k]*ri*x[j]*pow(ri*x[j],l[k]) )*dphi*pow(ri*y[j],m[k])*pow(ri*z[j]-rA,n[k]);
delphiy[I] += ( m[k]*pow(ri*y[j],m_1[k]) -     2.0*Alpha[k]*ri*y[j]*pow(ri*y[j],m[k]) )*dphi*pow(ri*x[j],l[k])*pow(ri*z[j]-rA,n[k]);
delphiz[I] += ( n[k]*pow(ri*z[j]-rA,n_1[k]) -  2.0*Alpha[k]*(ri*z[j]-rA)*pow(ri*z[j]-rA,n[k]) )*dphi*pow(ri*x[j],l[k])*pow(ri*y[j],m[k]);

phi[I] += dphi*pow(ri*x[j],l[k])*pow(ri*y[j],m[k])*pow(ri*z[j]-rA,n[k]);

//cout << "l[k]= " <<  l[k] << " l_1[k]= "<<l_1[k]<< "  m[k]= " <<  m[k] << "  m_1[k]= "<<m_1[k]  << "  n[k]= " <<  n[k] << " n_1[k]= "<<n_1[k]<<  endl;

//cout << "  delphix[I]= " <<  delphix[I] << "  delphiy[I]= " <<  delphiy[I] << "  delphiz[I]= " <<  delphiz[I] <<   endl;

/***************  Maybe this should end the integration ************/
} //////////////// END DO k =  Contraction over Nbasis ///////////////////
//          dpsi =


          psi += C[I]*phi[I];

          delpsix +=  C[I]*delphix[I];
          delpsiy +=  C[I]*delphiy[I];
          delpsiz +=  C[I]*delphiz[I];
} ///////// END DO I =  Sum over basis states N   //////////////////

          D = pow(2*psi*psi,3.0/3.0);
          rho = D;
          delrhox = 4*psi*delpsix;
          delrhoy = 4*psi*delpsiy;
          delrhoz = 4*psi*delpsiz;
          //  START Multicentre three-dimensional integration -- weight function
          mu = (ri - rB)/Dist;

          if (atoms==2){
          s1 =  StepS(mu);
          s2 =  StepS(-mu);
          wA = s1/(  s1 + s2 );          
          wB = s2/(  s1 + s2 );
          }else wA = 1.0;

          Electrons +=  W[j]*wA*w*rho;
          gamma = rho*8.0*(delpsix*delpsix + delpsiy*delpsiy + delpsiz*delpsiz);
          int_gamma +=  W[j]*wA*w*gamma;

//          int_gamma = int_gamma  + w*gamma;
          f=0.0;
          dfdra=0.0;
          dfdgaa=0.0;
          dfdgab=0.0;
          fx=0.0;
          dfxdra=0.0;
          dfxdgaa=0.0;
          dfxdgab=0.0;
          fc=0.0;
          dfcdra=0.0;
          dfcdgaa=0.0;
          dfcdgab=0.0;

          if (functional==1){
          x_rks_b88(rho,gamma,fx,dfxdra,dfxdgaa);
          c_rks_lyp(rho,gamma,fc,dfcdra,dfcdgaa,dfdgab);


          }else if (functional==2){
          x_rks_slater(rho,fx,dfxdra);

          }else if (functional==3){
          x_rks_b88(rho,gamma,fx,dfxdra,dfxdgaa);
          
          }else if (functional==4){
          x_rks_b88(rho,gamma,fx,dfxdra,dfxdgaa);
          c_rks_vwnrpa(rho,fc,dfcdra);

/*          }else if (functional==4){
          lyp_wt   =  0.81;
          vwnrpa_wt=  0.19;
          x_rks_b88(rho,gamma,fx,dfxdra,dfxdgaa);

          c_rks_lyp(rho,gamma,f,dfdra,dfdgaa,dfdgab);
          fc =lyp_wt*f;
          dfcdra =lyp_wt*dfdra;
          dfcdgaa =lyp_wt*dfdgaa;
          dfcdgab =lyp_wt*dfdgab;

          c_rks_vwnrpa(rho,f,dfdra);
         fc+=vwnrpa_wt*f;
          dfcdra+=vwnrpa_wt*dfdra;
          dfcdgaa+=vwnrpa_wt*dfdgaa;
          dfcdgab+=vwnrpa_wt*dfdgab;
*/
          }else if (functional==5){
           x_rks_b88(rho,gamma,fx,dfxdra,dfxdgaa);
           c_rks_pz81(rho,fc,dfcdra);


          }else if (functional==6){
            x_rks_pbe(rho,gamma,fx,dfxdra,dfxdgaa);
            c_rks_pz81(rho,fc,dfcdra);

          }

          if ((rho==0.0)){
          f=0.0;
          dfdra=0.0;
          dfdgaa=0.0;
          dfdgab=0.0;
          fx=0.0;
          dfxdra=0.0;
          dfxdgaa=0.0;
          dfxdgab=0.0;
          fc=0.0;
          dfcdra=0.0;
          dfcdgaa=0.0;
          dfcdgab=0.0;
}
          energyX+= 4.0*PI*W[j]*wA*w*fx;
          energyC+= 4.0*PI*W[j]*wA*w*fc;
          energyXC= energyX + energyC;
          f = fx + fc;
          dfdra = dfxdra + dfcdra;
          dfdgaa = dfxdgaa + dfcdgaa;
          dfdgab = dfxdgab + dfcdgab;

//          int_Vxc_rho =  int_Vxc_rho + w*rho*dfdra;

//          energyXC+= 4.0*PI*W[j]*wA*w*f;
      

            
          int_Vxc_rho =  int_Vxc_rho + W[j]*wA*w*rho*dfdra;
//          energyXC =  energyXC + W[j]*wA*w*f;
//cout << "blyp   rho= " <<  rho << " gamma= "<<gamma<< " f= " <<  f << " dfdra= " <<  dfdra  << " dfdgaa= " <<  dfdgaa<< " dfdgab= " <<  dfdgab<<  endl;

          for(int R = 0 ; R <= MSize-1 ; R++){
             for(int S = 0 ; S <= R /* !!! */; S++){
               fock = dfdra*phi[R]*phi[S]  +  ( dfdgaa + 0.5*dfdgab )*
               (delrhox*(phi[S]*delphix[R] + phi[R]*delphix[S]) +
                delrhoy*(phi[S]*delphiy[R] + phi[R]*delphiy[S]) +
                delrhoz*(phi[S]*delphiz[R] + phi[R]*delphiz[S]) );
               Fxc[R][S] = Fxc[R][S] + W[j]*wA*w*fock;
          }
              // Fxc[R][R] = Fxc[R][R] + W[j]*wA*w*(phi[R]*phi[R]);
          }



} ///////////// END DO for Atomic Loop here   //////////////

} // END DO i =  radial integration


} // END DO j =   Angular integration
       cout <<" Ratom =  "<< Ratom <<"  nr =  "<< nr <<"  sum =  "<< 4*PI*sum <<"  ri =  "<< ri << endl;
       cout << "  # of Electrons +=  W[j]*wA*w*rho=  "<<4*PI*Electrons<<  endl;
	cout <<"  int_Vxc_rho =  " << 4*PI*int_Vxc_rho  << endl;
       cout <<"  int_gamma =  " << 4*PI*int_gamma  << endl;
       for(int I = 0 ; I <= MSize-1 ; I++){
cout << "  C["<<I<<"]   " <<  C[I] << "  Alpha["<<I<<"]   " <<  Alpha[I]  << "  d["<<I<<"]   " <<  d[I]<<  endl;
}


      for(int I = 0 ; I <= MSize-1 ; I++){
        for(int J = 0 ; J <= I ; J++){
//cout << " FMatrix["<<I<<"]["<<J<<"] =   "<< FMatrix[I][J] << "  Fxc["<<I<<"]["<<J<<"] = " <<4*PI*Fxc[I][J] << endl;
FMatrix[I][J] = FMatrix[I][J] + 4*PI*Fxc[I][J];
}}
//fout.close();
cout << " energyXC =   "<< energyXC << endl;
cout << " energyX =   "<< energyX << " energyC =   "<< energyC << endl;

       }
/////////////////////////////////////////////////////////////////////////

       double StepS(double mu){
       double  p,s;

       p =  3*mu/2  - pow(mu,3)/2;
       p =  3*p/2  - pow(p,3)/2;
       p =  3*p/2  - pow(p,3)/2;
       s =   (1- p)/2.0;
       return s;
}

/////////////////////////////////////////////////////////////////////////

       void CheckVMatrix(Mat_I_DP &VMatrix,Mat_I_DP &SMatrix){
cout << ".................. CheckVMatrix ......................" << endl;
      double norm;    //  Coeffs[NSize],

//      int R, S,i,j;
/*
       for(int R = 0 ; R <= NSize -1; R++){
        for(int S = 0 ; S <= R -1; S++){
           SMatrix[S][R] = SMatrix[R][S];
          }}
*/

     for(int i = 0 ; i <= MSize -1; i++){
     for(int j = 0 ; j <= MSize -1; j++){
       norm = 0;
      for(int R = 0 ; R <= MSize -1; R++){
        for(int S = 0 ; S <= MSize -1; S++){
//      DO R=1, NSize
//        DO S=1, NSize
         norm = norm + VMatrix[i][R]*SMatrix[R][S]*VMatrix[j][S];
//          cout << "Norm Coeffs[R]  " << Norm << "   " << Coeffs[R] << endl;
//        END DO
//      END DO
}}
//        if ((abs(norm) > 1.0e-14) || (abs(VMatrix[i][j]) > 1.0e4) ){
            cout << "i  j  " << i << "   " << j << " norm = " << norm<< "   SMatrix =  " << SMatrix[i][j]<< "  VMatrix = " << VMatrix[i][j]<< endl;
//}
}}

}