Cusp dressing almost OK

plugins/Hartree_Fock/EZFIO.cfg

@@ -26,13 +26,13 @@ default: 1.e-12
 type: Strictly_positive_int
 doc: Maximum number of SCF iterations
 interface: ezfio,provider,ocaml
-default: 128
+default: 500
 
 [level_shift]
 type: Positive_float
 doc: Energy shift on the virtual MOs to improve SCF convergence
 interface: ezfio,provider,ocaml
-default: 0.0
+default: 0.3
 
 [scf_algorithm]
 type: character*(32)

plugins/Hartree_Fock/Roothaan_Hall_SCF.irp.f

@@ -71,7 +71,12 @@ END_DOC
 
       Fock_matrix_AO_alpha = Fock_matrix_AO*0.5d0
       Fock_matrix_AO_beta  = Fock_matrix_AO*0.5d0
+      level_shift_save = level_shift
+      if (max_error_DIIS > 0.1d0) then
+        level_shift += 0.5d0
+      endif
       TOUCH Fock_matrix_AO_alpha Fock_matrix_AO_beta
+      level_shift = level_shift_save
 
     endif
 
@@ -97,6 +102,9 @@ END_DOC
 
     double precision :: level_shift_save
     level_shift_save = level_shift
+    if (max_error_DIIS > 0.5d0) then
+      level_shift += 0.5d0
+    endif
     do while (Delta_Energy_SCF > 0.d0)
       level_shift = level_shift + 0.1d0
       MO_coef = eigenvectors_Fock_matrix_MO

plugins/Hartree_Fock_SlaterDressed/LinearSystem.irp.f

@@ -11,6 +11,7 @@ BEGIN_PROVIDER [ double precision, cusp_A, (nucl_num, nucl_num) ]
     cusp_A(A,A) = slater_expo(A)/nucl_charge(A) * slater_value_at_nucl(A,A) 
     do B=1,nucl_num
       cusp_A(A,B) -= slater_value_at_nucl(B,A)
+      ! Projector
       do mu=1,mo_tot_num
         cusp_A(A,B) += MOSlaOverlap_matrix(mu,B) * mo_value_at_nucl(mu,A)
       enddo

plugins/Hartree_Fock_SlaterDressed/SCF_dressed.irp.f

@@ -76,14 +76,21 @@ subroutine run
   double precision               :: EHF
   integer                        :: i_it, i, j, k
    
-  EHF = HF_energy 
-
   mo_label = "CuspDressed"
 
+
+  print *,  HF_energy
+
   call ezfio_set_Hartree_Fock_SlaterDressed_slater_coef_ezfio(cusp_C)
-! Choose SCF algorithm
 
+  do i=1,ao_num
+    print *,  mo_coef(i,1), cusp_corrected_mos(i,1)
+  enddo
+  mo_coef(1:ao_num,1:mo_tot_num) = cusp_corrected_mos(1:ao_num,1:mo_tot_num)
+  TOUCH mo_coef
 
+  EHF = HF_energy 
+  print *,  HF_energy
 !  call Roothaan_Hall_SCF
   
 end

plugins/Hartree_Fock_SlaterDressed/integrals.irp.f

@@ -24,44 +24,93 @@ subroutine GauSlaOverlap(expGau,cGau,aGau,expSla,cSla,result)
   
 ! calculate the length AB between the two centers and other usful quantities
 
-  AB = (cGau(1)-cSla(1))**2d0 + (cGau(2)-cSla(2))**2d0 + (cGau(3)-cSla(3))**2d0
+  AB = (cGau(1)-cSla(1))**2 + (cGau(2)-cSla(2))**2 + (cGau(3)-cSla(3))**2
-  AB = sqrt(AB)
+  AB = dsqrt(AB)
 
   AxBx = (cGau(1)-cSla(1))/2d0
   AyBy = (cGau(2)-cSla(2))/2d0
   AzBz = (cGau(3)-cSla(3))/2d0
+  ds = 0.d0
 
 ! intermediate variables
 
-  t = expSla*sqrt(0.25d0/expGau)
+  t = expSla*dsqrt(0.25d0/expGau)
-  u = sqrt(expGau)*AB
+  u = dsqrt(expGau)*AB
 
+  double precision :: d, et2
   if(AB > 0d0) then 
 
 !   (s|s) 
-    ss = (t+u)*erfc(t+u)*exp(2d0*t*(t+u)) - (t-u)*erfc(t-u)*exp(2d0*t*(t-u))
+    ss = 0.d0
+
+    d = derfc(t+u)
+    if (dabs(d) > 1.d-30) then
+      ss  = (t+u)*d*dexp(2d0*t*(t+u)) 
+    endif
+
+    d = derfc(t-u)
+    if (dabs(d) > 1.d-30) then
+      ss -= (t-u)*d*dexp(2d0*t*(t-u))
+    endif
 
 !   (p|s)
-    ps = (exp(t**2d0-u**2d0)*(-4d0*t+sqrt(pi)*(exp((t-u)**2d0)*(1d0+2d0*t*(t-u))*erfc(t-u) & 
+    ps = 0.d0
-       + exp((t+u)**2d0)*(1d0+2d0*t*(t+u))*erfc(t+u))))/sqrt(pi)
+    if (t*t-u*u > 300.d0) then
+      et2 = huge(1.0)
+    else
+      et2 = dexp(t*t-u*u)
+    endif
+    if (et2 /= 0.d0) then
+      d = derfc(t-u) 
+      if (d /= 0.d0) then
+        ps += dexp((t-u)**2)*(1d0+2d0*t*(t-u))*d
+      endif
+      d = derfc(t+u)
+      if (d /= 0.d0) then
+        ps += dexp((t+u)**2)*(1d0+2d0*t*(t+u))*d
+      endif
+      ps *= dsqrt(pi)
+      ps -= 4d0*t 
+      ps *= et2/dsqrt(pi)
+    endif
 
 !   (d|s)
-    ds = 4d0*exp(2d0*t*(t-u))*t*(-((1d0+t**2d0-t*u)*erfc(t-u))+exp(4d0*t*u)*(1d0+t*(t+u))*erfc(t+u))
+!    ds = 4d0*dexp(2d0*t*(t-u))*t*(-((1d0+t**2-t*u)*derfc(t-u))+dexp(4d0*t*u)*(1d0+t*(t+u))*derfc(t+u))
+    ds = 0.d0
+    d = derfc(t+u)
+    if (d /= 0.d0) then
+      ds  = dexp(4d0*t*u)*(1d0+t*(t+u))*d
+    endif
+    d = derfc(t-u)
+    if (d /= 0.d0) then
+      ds -= (1d0+t*t-t*u)*d
+    endif
+
+    if ( dabs(ds) > 1.d-100) then
+      ds *= 4d0*dexp(2d0*t*(t-u))*t
+    endif
 
 !   backward scaling
     ds = 3d0*ss/u**5d0 - 3d0*ps/u**4d0 + ds/u**3d0
-    ps = ps/u**2d0-ss/u**3d0
+    ps = ps/u**2-ss/u**3d0
     ss = ss/u
 
   else
 
 !   concentric case
-    ss = 2d0*exp(t**2d0)*((-2d0*t)/sqrt(pi)+exp(t**2d0)*(1d0+2d0*t**2d0)*erfc(t))
+    d = derfc(t)
-    ps = (8d0*exp(t**2d0)*t*(-2d0*(1d0+t**2d0)+exp(t**2d0)*sqrt(pi)*t*(3d0+2d0*t**2d0)*erfc(t)))/(3d0*sqrt(pi))
+    if (d /= 0.d0) then
+      et2 = dexp(t*t)
+      ss = 2d0*et2*((-2d0*t)/dsqrt(pi)+et2*(1d0+2d0*t*t)*d)
+      ps = (8d0*et2*t*(-2d0*(1d0+t*t)+et2*dsqrt(pi)*t*(3d0+2d0*t*t)*d))/(3d0*dsqrt(pi))
+    else
+      ss = 0.d0
+      ps = 0.d0
+    endif
 
   endif
 
-  k = t**3d0*exp(-t**2d0)*4d0*pi/expSla**(3d0/2d0)
+  k = t**3d0*dexp(-t*t)*4d0*pi/expSla**(3d0/2d0)
 
 ! (s|s) 
   ss = k*ss
@@ -76,9 +125,134 @@ subroutine GauSlaOverlap(expGau,cGau,aGau,expSla,cSla,result)
 ! (d|s) 
   ds = k*ds
 
-  dxxs = (2d0*ss+ps)/(4d0*expGau) + AxBx**2d0*ds
+  dxxs = (2d0*ss+ps)/(4d0*expGau) + AxBx**2*ds
-  dyys = (2d0*ss+ps)/(4d0*expGau) + AyBy**2d0*ds
+  dyys = (2d0*ss+ps)/(4d0*expGau) + AyBy**2*ds
-  dzzs = (2d0*ss+ps)/(4d0*expGau) + AzBz**2d0*ds
+  dzzs = (2d0*ss+ps)/(4d0*expGau) + AzBz**2*ds
+
+  dxys = AxBx*AyBy*ds
+  dxzs = AxBx*AzBz*ds
+  dyzs = AyBy*AzBz*ds
+
+  select case (sum(aGau))
+    case (0)
+      result = ss
+
+    case (1)
+      if (aGau(1) == 1) then
+        result = pxs
+      else if (aGau(2) == 1) then
+        result = pys
+      else if (aGau(3) == 1) then
+        result = pzs
+      endif
+
+    case (2)
+      if (aGau(1) == 2) then
+        result = dxxs
+      else if (aGau(2) == 2) then
+        result = dyys
+      else if (aGau(3) == 2) then
+        result = dzzs
+      else if (aGau(1)+aGau(2) == 2) then
+        result = dxys
+      else if (aGau(1)+aGau(3) == 2) then
+        result = dxzs
+      else if (aGau(2)+aGau(3) == 2) then
+        result = dyzs
+      endif
+
+    case default
+      stop 'GauSlaOverlap not implemented'
+
+  end select
+
+end
+!*****************************************************************************
+
+!*****************************************************************************
+subroutine GauSlaKinetic(expGau,cGau,aGau,expSla,cSla,result)
+
+  implicit none
+
+  BEGIN_DOC
+  ! Compute the kinetic energy integral between a Gaussian function
+  ! with arbitrary angular momemtum and a s-type Slater function
+  END_DOC
+
+! Input variables 
+  double precision,intent(in)   :: expGau,expSla
+  double precision,intent(in)   :: cGau(3),cSla(3)
+  integer,intent(in)            :: aGau(3)
+  double precision,intent(out)  :: result
+
+! Final value of the integrals
+  double precision              :: ss,ps,ds
+  double precision              :: pxs,pys,pzs
+  double precision              :: dxxs,dyys,dzzs,dxys,dxzs,dyzs
+
+  double precision              :: pi,E,AB,AxBx,AyBy,AzBz,t,u,k
+
+  pi = 4d0*atan(1d0)
+  
+! calculate the length AB between the two centers
+
+  AB = (cGau(1)-cSla(1))**2 + (cGau(2)-cSla(2))**2 + (cGau(3)-cSla(3))**2
+  AB = dsqrt(AB)
+
+  AxBx = (cGau(1)-cSla(1))/2d0
+  AyBy = (cGau(2)-cSla(2))/2d0
+  AzBz = (cGau(3)-cSla(3))/2d0
+
+! intermediate variables
+
+  t = expSla*dsqrt(0.25d0/expGau)
+  u = dsqrt(expGau)*AB
+
+  if(AB > 0d0) then 
+
+!   (s|s) 
+    ss = (1d0+t*(t-u))*derfc(t-u)*dexp(2d0*t*(t-u)) - (1d0+t*(t+u))*derfc(t+u)*dexp(2d0*t*(t+u))
+
+!   (p|s)
+    ps = (dexp(t**2-2d0*t*u-u**2)*(4d0*dexp(2d0*t*u)*(1d0+t**2)    & 
+       + dsqrt(pi)*t*(-(dexp(t**2+u**2)*(3d0+2d0*t*(t-u))*derfc(t-u)) & 
+       - dexp(2d0*t*u+(t+u)**2)*(3d0+2d0*t*(t+u))*derfc(t+u))))/dsqrt(pi)
+
+!   (d|s)
+    ds = (-8d0*dexp(t**2-u**2)*u+4d0*dexp(2d0*t*(t-u))*dsqrt(pi)*t**2*((2d0+t**2-t*u)*derfc(t-u) &
+       - dexp(4d0*t*u)*(2d0+t*(t+u))*derfc(t+u)))/dsqrt(pi)
+
+!   backward scaling
+    ds = 3d0*ss/u**5d0 - 3d0*ps/u**4d0 + ds/u**3d0
+    ps = ps/u**2-ss/u**3d0
+    ss = ss/u
+
+  else
+
+!   concentric case
+    ss = (4d0*dexp(t**2)*(1d0+t**2))/dsqrt(pi)-2d0*dexp(2d0*t**2)*t*(3d0+2d0*t**2)*derfc(t)
+    ps = (8d0*dexp(t**2)*(-1d0+4d0*t**2+2d0*t**4d0-dexp(t**2)*dsqrt(pi)*t**3d0*(5d0+2d0*t**2)*derfc(t)))/(3d0*dsqrt(pi))
+
+  endif
+
+  k = expSla*dsqrt(expGau)*t**3d0*dexp(-t*t)*4d0*pi/expSla**(3d0/2d0)
+
+! (s|s) 
+  ss = k*ss
+
+! (p|s) 
+  ps = k*ps
+
+  pxs = AxBx*ps
+  pys = AyBy*ps
+  pzs = AzBz*ps
+
+! (d|s) 
+  ds = k*ds
+
+  dxxs = (2d0*ss+ps)/(4d0*expGau) + AxBx**2*ds
+  dyys = (2d0*ss+ps)/(4d0*expGau) + AyBy**2*ds
+  dzzs = (2d0*ss+ps)/(4d0*expGau) + AzBz**2*ds
 
   dxys = AxBx*AyBy*ds
   dxzs = AxBx*AzBz*ds
@@ -121,108 +295,9 @@ end
 !*****************************************************************************
 
 
-!*****************************************************************************
-subroutine GauSlaKinetic(expGau,cGau,aGau,expSla,cSla)
-
-  implicit none
-
-  BEGIN_DOC
-  ! Compute the kinetic energy integral between a Gaussian function
-  ! with arbitrary angular momemtum and a s-type Slater function
-  END_DOC
-
-! Input variables 
-  double precision,intent(in)   :: expGau,expSla
-  double precision,intent(in)   :: cGau(3),cSla(3)
-  integer,intent(in)            :: aGau(3)
-
-! Final value of the integrals
-  double precision              :: ss,ps,ds
-  double precision              :: pxs,pys,pzs
-  double precision              :: dxxs,dyys,dzzs,dxys,dxzs,dyzs
-
-  double precision              :: pi,E,AB,AxBx,AyBy,AzBz,t,u,k
-
-  pi = 4d0*atan(1d0)
-  
-! calculate the length AB between the two centers
-
-  AB = (cGau(1)-cSla(1))**2d0 + (cGau(2)-cSla(2))**2d0 + (cGau(3)-cSla(3))**2d0
-  AB = sqrt(AB)
-
-  AxBx = (cGau(1)-cSla(1))/2d0
-  AyBy = (cGau(2)-cSla(2))/2d0
-  AzBz = (cGau(3)-cSla(3))/2d0
-
-! intermediate variables
-
-  t = expSla*sqrt(0.25d0/expGau)
-  u = sqrt(expGau)*AB
-
-  if(AB > 0d0) then 
-
-!   (s|s) 
-    ss = (1d0+t*(t-u))*erfc(t-u)*exp(2d0*t*(t-u)) - (1d0+t*(t+u))*erfc(t+u)*exp(2d0*t*(t+u))
-
-!   (p|s)
-    ps = (exp(t**2d0-2d0*t*u-u**2d0)*(4d0*exp(2d0*t*u)*(1d0+t**2d0)    & 
-       + sqrt(pi)*t*(-(exp(t**2d0+u**2d0)*(3d0+2d0*t*(t-u))*erfc(t-u)) & 
-       - exp(2d0*t*u+(t+u)**2d0)*(3d0+2d0*t*(t+u))*erfc(t+u))))/sqrt(pi)
-
-!   (d|s)
-    ds = (-8d0*exp(t**2d0-u**2d0)*u+4d0*exp(2d0*t*(t-u))*sqrt(pi)*t**2d0*((2d0+t**2d0-t*u)*erfc(t-u) &
-       - exp(4d0*t*u)*(2d0+t*(t+u))*erfc(t+u)))/sqrt(pi)
-
-!   backward scaling
-    ds = 3d0*ss/u**5d0 - 3d0*ps/u**4d0 + ds/u**3d0
-    ps = ps/u**2d0-ss/u**3d0
-    ss = ss/u
-
-  else
-
-!   concentric case
-    ss = (4d0*exp(t**2d0)*(1d0+t**2d0))/sqrt(pi)-2d0*exp(2d0*t**2d0)*t*(3d0+2d0*t**2d0)*erfc(t)
-    ps = (8d0*exp(t**2d0)*(-1d0+4d0*t**2d0+2d0*t**4d0-exp(t**2)*sqrt(pi)*t**3d0*(5d0+2d0*t**2d0)*erfc(t)))/(3d0*sqrt(pi))
-
-  endif
-
-  k = expSla*sqrt(expGau)*t**3d0*exp(-t**2d0)*4d0*pi/expSla**(3d0/2d0)
-
-! (s|s) 
-  ss = k*ss
-
-! (p|s) 
-  ps = k*ps
-
-  pxs = AxBx*ps
-  pys = AyBy*ps
-  pzs = AzBz*ps
-
-! (d|s) 
-  ds = k*ds
-
-  dxxs = (2d0*ss+ps)/(4d0*expGau) + AxBx**2d0*ds
-  dyys = (2d0*ss+ps)/(4d0*expGau) + AyBy**2d0*ds
-  dzzs = (2d0*ss+ps)/(4d0*expGau) + AzBz**2d0*ds
-
-  dxys = AxBx*AyBy*ds
-  dxzs = AxBx*AzBz*ds
-  dyzs = AyBy*AzBz*ds
-
-! Print result
-  write(*,'(A12,F16.10)') & 
-    '(s|T|s) = ',ss
-  write(*,'(A12,F16.10,3X,A12,F16.10,3X,A12,F16.10)') & 
-    '(px|T|s) = ',pxs,'(py|T|s) = ',pys,'(pz|T|s) = ',pzs
-  write(*,'(A12,F16.10,3X,A12,F16.10,3X,A12,F16.10,3X,A12,F16.10,3X,A12,F16.10,3X,A12,F16.10)') & 
-    '(dx2|T|s) = ',dxxs,'(dy2|T|s) = ',dyys,'(dz2|T|s) = ',dzzs,'(dxy|T|s) = ',dxys,'(dxz|T|s) = ',dxzs,'(dyz|T|s) = ',dyzs
-
-end
-!*****************************************************************************
-
 
 !*****************************************************************************
-subroutine GauSlaNuclear(expGau,cGau,aGau,expSla,cSla,ZNuc,cNuc)
+subroutine GauSlaNuclear(expGau,cGau,aGau,expSla,cSla,ZNuc,cNuc,result)
 
   implicit none
 
@@ -237,6 +312,7 @@ subroutine GauSlaNuclear(expGau,cGau,aGau,expSla,cSla,ZNuc,cNuc)
   integer,intent(in)            :: aGau(3)
   double precision,intent(in)   :: cNuc(3)
   double precision,intent(in)   :: ZNuc
+  double precision,intent(out)  :: result
 
 ! Final value of the overlap integral
   double precision              :: ss,ps,ds,fs
@@ -249,29 +325,31 @@ subroutine GauSlaNuclear(expGau,cGau,aGau,expSla,cSla,ZNuc,cNuc)
   
 ! calculate the length AB between the two centers
 
-  AB = (cGau(1)-cSla(1))**2d0 + (cGau(2)-cSla(2))**2d0 + (cGau(3)-cSla(3))**2d0
+  AB = (cGau(1)-cSla(1))**2 + (cGau(2)-cSla(2))**2 + (cGau(3)-cSla(3))**2
-  AB = sqrt(AB)
+  AB = dsqrt(AB)
 
 ! intermediate variables
 
-  x = sqrt(expSla**2d0/(4d0*expGau))
+  x = dsqrt(expSla**2/(4d0*expGau))
-  y = sqrt(expGau)*AB
+  y = dsqrt(expGau)*AB
 
   if(AB > 0d0) then 
-    ss = (1d0+x*(x+y))*erfc(x+y)*exp(2d0*x*(x+y)) - (1d0+x*(x-y))*erfc(x-y)*exp(2d0*x*(x-y))
+    ss = (1d0+x*(x+y))*derfc(x+y)*dexp(2d0*x*(x+y)) - (1d0+x*(x-y))*derfc(x-y)*dexp(2d0*x*(x-y))
     ss = ss/y
   else
-    ss = (4d0*E**x**2d0*(1d0+x**2d0))/sqrt(Pi)-2d0*E**(2d0*x**2d0)*x*(3d0+2d0*x**2d0)*Erfc(x)
+    ss = (4d0*E**x**2*(1d0+x**2))/dsqrt(Pi)-2d0*E**(2d0*x**2)*x*(3d0+2d0*x**2)*dErfc(x)
   endif
 
-  k = expSla*sqrt(expGau)*x**3d0*exp(-x**2)*4d0*pi/expSla**(3d0/2d0)
+  k = expSla*dsqrt(expGau)*x**3d0*dexp(-x*x)*4d0*pi/expSla**(3d0/2d0)
   ss = k*ss
 
 ! Print result
   write(*,*) ss
+  result = 0.d0
 
 end
 !*****************************************************************************
+
 double precision function BoysF0(t)
   implicit none
   double precision, intent(in)  :: t
@@ -280,7 +358,7 @@ double precision function BoysF0(t)
   pi = 4d0*atan(1d0)
 
   if(t > 0d0) then
-    BoysF0 = 0.5d0*sqrt(pi/t)*erf(sqrt(t))
+    BoysF0 = 0.5d0*dsqrt(pi/t)*derf(dsqrt(t))
   else
     BoysF0 = 1d0
   endif
@@ -288,8 +366,81 @@ double precision function BoysF0(t)
 end
 !*****************************************************************************
 
+!TODO
+subroutine GauSlaOverlap_write(expGau,cGau,aGau,expSla,cSla,result,iunit)
+  implicit none
+  double precision,intent(in)   :: expGau,expSla
+  double precision,intent(in)   :: cGau(3),cSla(3)
+  integer,intent(in)            :: aGau(3)
+  integer,intent(in)            :: iunit
+  double precision,intent(out)  :: result
+  write(iunit, *)  &
+  'SDrSla[ {',expGau,',{',cGau(1),',',cGau(2),',',cGau(3),'},{',aGau(1),',',aGau(2),',',aGau(3),'} },{', expSla,', {',cSla(1),',',cSla(2),',',cSla(3),'} } ],'
+  result = 0.d0
+end
 
-BEGIN_PROVIDER [ double precision, GauSlaOverlap_matrix, (ao_num, nucl_num) ]
+subroutine GauSlaOverlap_read(expGau,cGau,aGau,expSla,cSla,result,iunit)
+  implicit none
+  double precision,intent(in)   :: expGau,expSla
+  double precision,intent(in)   :: cGau(3),cSla(3)
+  integer,intent(in)            :: aGau(3)
+  integer,intent(in)            :: iunit
+  double precision,intent(out)  :: result
+  read(iunit, *)  result
+end
+
+subroutine GauSlaKinetic_write(expGau,cGau,aGau,expSla,cSla,result,iunit)
+  implicit none
+  double precision,intent(in)   :: expGau,expSla
+  double precision,intent(in)   :: cGau(3),cSla(3)
+  integer,intent(in)            :: aGau(3)
+  integer,intent(in)            :: iunit
+  double precision,intent(out)  :: result
+  write(iunit, *)  &
+  'TDrSla[ {',expGau,',{',cGau(1),',',cGau(2),',',cGau(3),'},{',aGau(1),',',aGau(2),',',aGau(3),'} },{', expSla,',{',cSla(1),',',cSla(2),',',cSla(3),'} } ],'
+  result = 0.d0
+end
+
+subroutine GauSlaKinetic_read(expGau,cGau,aGau,expSla,cSla,result,iunit)
+  implicit none
+  double precision,intent(in)   :: expGau,expSla
+  double precision,intent(in)   :: cGau(3),cSla(3)
+  integer,intent(in)            :: aGau(3)
+  integer,intent(in)            :: iunit
+  double precision,intent(out)  :: result
+  read(iunit, *)  result
+end
+
+subroutine GauSlaNuclear_write(expGau,cGau,aGau,expSla,cSla,ZNuc,cNuc,result,iunit)
+  implicit none
+  double precision,intent(in)   :: expGau,expSla
+  double precision,intent(in)   :: cGau(3),cSla(3)
+  integer,intent(in)            :: aGau(3)
+  double precision,intent(in)   :: cNuc(3)
+  double precision,intent(in)   :: ZNuc
+  integer,intent(in)            :: iunit
+  double precision,intent(out)  :: result
+  write(iunit, *)  &
+  'VDrSla[ {',expGau,',{',cGau(1),',',cGau(2),',',cGau(3),'},{',aGau(1),',',aGau(2),',',aGau(3),'} },{ ', expSla,',{',cSla(1),',',cSla(2),',',cSla(3),'} }, {', ZNuc, ',{', cNuc(1),',', cNuc(2),',', cNuc(3), '} } ],'
+  result = 0.d0
+end
+
+subroutine GauSlaNuclear_read(expGau,cGau,aGau,expSla,cSla,ZNuc,cNuc,result,iunit)
+  implicit none
+  double precision,intent(in)   :: expGau,expSla
+  double precision,intent(in)   :: cGau(3),cSla(3)
+  integer,intent(in)            :: aGau(3)
+  double precision,intent(in)   :: cNuc(3)
+  double precision,intent(in)   :: ZNuc
+  integer,intent(in)            :: iunit
+  double precision,intent(out)  :: result
+  read(iunit, *)  result
+end
+!TODO
+
+BEGIN_TEMPLATE
+
+BEGIN_PROVIDER [ double precision, GauSla$X_matrix, (ao_num, nucl_num) ]
   implicit none
   BEGIN_DOC
   ! <Gaussian | Slater> overlap matrix
@@ -300,6 +451,22 @@ BEGIN_PROVIDER [ double precision, GauSlaOverlap_matrix, (ao_num, nucl_num) ]
   double precision               :: expSla, res, expGau
   integer                        :: aGau(3)
 
+!TODO
+  logical :: read
+  integer :: iunit
+  integer :: getunitandopen
+
+  inquire(FILE=trim(ezfio_filename)//'/work/GauSla$X.dat',EXIST=read)
+  if (read) then
+    print *,  'READ $X'
+    iunit = getunitandopen(trim(ezfio_filename)//'/work/GauSla$X.dat','r')
+  else
+    print *,  'WRITE $X'
+    iunit = getunitandopen(trim(ezfio_filename)//'/work/GauSla$X.inp','w')
+    write(iunit,*) '{'
+  endif
+!TODO
+
   do k=1,nucl_num
     cSla(1:3) = nucl_coord_transp(1:3,k)
     expSla    = slater_expo(k)
@@ -307,28 +474,151 @@ BEGIN_PROVIDER [ double precision, GauSlaOverlap_matrix, (ao_num, nucl_num) ]
     do i=1,ao_num
       cGau(1:3) = nucl_coord_transp(1:3, ao_nucl(i))
       aGau(1:3) = ao_power(i,1:3)
-      GauSlaOverlap_matrix(i,k) = 0.d0
+      GauSla$X_matrix(i,k) = 0.d0
 
       do j=1,ao_prim_num(i)
         expGau = ao_expo_ordered_transp(j,i)
-        call GauSlaOverlap(expGau,cGau,aGau,expSla,cSla,res)
+!        call GauSla$X(expGau,cGau,aGau,expSla,cSla,res)
-        GauSlaOverlap_matrix(i,k) += ao_coef_normalized_ordered_transp(j,i) * res
+        if (read) then
+          call GauSla$X_read(expGau,cGau,aGau,expSla,cSla,res,iunit)
+        else
+          call GauSla$X_write(expGau,cGau,aGau,expSla,cSla,res,iunit)
+        endif
+        GauSla$X_matrix(i,k) += ao_coef_normalized_ordered_transp(j,i) * res
       enddo
 
     enddo
 
   enddo
+  if (.not.read) then
+    write(iunit,*) '0.}'
+  endif
+  close(iunit)
   
 END_PROVIDER
 
-BEGIN_PROVIDER [ double precision, MOSlaOverlap_matrix, (mo_tot_num, nucl_num) ]
+BEGIN_PROVIDER [ double precision, MOSla$X_matrix, (mo_tot_num, nucl_num) ]
   implicit none
   BEGIN_DOC
 ! <MO | Slater>
   END_DOC
   call dgemm('N','N',mo_tot_num,nucl_num,ao_num,1.d0,                &
       mo_coef_transp, size(mo_coef_transp,1),                        &
-      GauSlaOverlap_matrix, size(GauSlaOverlap_matrix,1),                    &
+      GauSla$X_matrix, size(GauSla$X_matrix,1),                    &
-      0.d0, MOSlaOverlap_matrix, size(MOSlaOverlap_matrix,1))
+      0.d0, MOSla$X_matrix, size(MOSla$X_matrix,1))
+END_PROVIDER
+
+BEGIN_PROVIDER [ double precision, AO_orthoSla$X_matrix, (ao_num, nucl_num) ]
+ implicit none
+ BEGIN_DOC
+ ! <AO_ortho | Slater>
+ END_DOC
+  call dgemm('T','N',ao_num,nucl_num,ao_num,1.d0,                &
+      ao_ortho_canonical_coef, size(ao_ortho_canonical_coef,1),      &
+      GauSla$X_matrix, size(GauSla$X_matrix,1),                      &
+      0.d0, AO_orthoSla$X_matrix, size(AO_orthoSla$X_matrix,1))
+ 
+END_PROVIDER
+
+SUBST [ X ]
+
+Overlap ;;
+Kinetic ;;
+
+END_TEMPLATE
+
+BEGIN_PROVIDER [ double precision, GauSlaNuclear_matrix, (ao_num, nucl_num) ]
+  implicit none
+  BEGIN_DOC
+  ! <Gaussian | Slater> overlap matrix
+  END_DOC
+  integer                        :: i,j,k,A
+  double precision               :: cGau(3)
+  double precision               :: cSla(3)
+  double precision               :: expSla, res, expGau, Znuc, cNuc(3)
+  integer                        :: aGau(3)
+
+!TODO
+  logical :: read
+  integer :: iunit
+  integer :: getunitandopen
+
+  inquire(FILE=trim(ezfio_filename)//'/work/GauSlaNuclear.dat',EXIST=read)
+  if (read) then
+    print *,  'READ Nuclear'
+    iunit = getunitandopen(trim(ezfio_filename)//'/work/GauSlaNuclear.dat','r')
+  else
+    print *,  'WRITE Nuclear'
+    iunit = getunitandopen(trim(ezfio_filename)//'/work/GauSlaNuclear.inp','w')
+    write(iunit,*)'{'
+  endif
+!TODO
+
+  do k=1,nucl_num
+    cSla(1:3) = nucl_coord_transp(1:3,k)
+    expSla    = slater_expo(k)
+
+    do i=1,ao_num
+      cGau(1:3) = nucl_coord_transp(1:3, ao_nucl(i))
+      aGau(1:3) = ao_power(i,1:3)
+      GauSlaNuclear_matrix(i,k) = 0.d0
+
+      do j=1,ao_prim_num(i)
+        expGau = ao_expo_ordered_transp(j,i)
+        do A=1,nucl_num
+          cNuc(1:3) = nucl_coord_transp(1:3,A)
+          Znuc = nucl_charge(A)
+!          call GauSlaNuclear(expGau,cGau,aGau,expSla,cSla,Znuc,cNuc,res)
+          if (read) then
+            call GauSlaNuclear_read(expGau,cGau,aGau,expSla,cSla,Znuc,cNuc,res,iunit)
+          else
+            call GauSlaNuclear_write(expGau,cGau,aGau,expSla,cSla,Znuc,cNuc,res,iunit)
+          endif
+          GauSlaNuclear_matrix(i,k) += ao_coef_normalized_ordered_transp(j,i) * res
+        enddo
+      enddo
+
+    enddo
+
+  enddo
+  if (.not.read) then
+    write(iunit,*) '0.}'
+  endif
+  close(iunit)
+  
+END_PROVIDER
+
+BEGIN_PROVIDER [ double precision, GauSlaH_matrix, (ao_num, nucl_num) ]
+ implicit none
+ BEGIN_DOC
+ ! Core hamiltonian in AO basis
+ END_DOC
+ GauSlaH_matrix(1:ao_num,1:nucl_num) =                               &
+     GauSlaKinetic_matrix(1:ao_num,1:nucl_num) +                     &
+     GauSlaNuclear_matrix(1:ao_num,1:nucl_num)
+
+END_PROVIDER
+
+
+BEGIN_PROVIDER [ double precision, MOSlaNuclear_matrix, (mo_tot_num, nucl_num) ]
+  implicit none
+  BEGIN_DOC
+! <MO | Slater>
+  END_DOC
+  call dgemm('N','N',mo_tot_num,nucl_num,ao_num,1.d0,                &
+      mo_coef_transp, size(mo_coef_transp,1),                        &
+      GauSlaNuclear_matrix, size(GauSlaNuclear_matrix,1),                    &
+      0.d0, MOSlaNuclear_matrix, size(MOSlaNuclear_matrix,1))
+END_PROVIDER
+
+BEGIN_PROVIDER [ double precision, AO_orthoSlaH_matrix, (ao_num, nucl_num) ]
+  implicit none
+  BEGIN_DOC
+! <AO ortho | Slater>
+  END_DOC
+  call dgemm('T','N',ao_num,nucl_num,ao_num,1.d0,                &
+      ao_ortho_canonical_coef, size(ao_ortho_canonical_coef,1),                        &
+      GauSlaH_matrix, size(GauSlaH_matrix,1),                        &
+      0.d0, AO_orthoSlaH_matrix, size(AO_orthoSlaH_matrix,1))
 END_PROVIDER
 

src/MOGuess/mo_ortho_lowdin.irp.f

@@ -6,7 +6,9 @@ BEGIN_PROVIDER [double precision, ao_ortho_lowdin_coef, (ao_num_align,ao_num)]
 ! ao_ortho_lowdin_coef(i,j) = coefficient of the ith ao on the jth ao_ortho_lowdin orbital
   END_DOC
   integer                        :: i,j,k,l
-  double precision               :: tmp_matrix(ao_num_align,ao_num),accu
+  double precision               :: accu
+  double precision, allocatable  :: tmp_matrix(:,:)
+  allocate (tmp_matrix(ao_num_align,ao_num))
   tmp_matrix(:,:) = 0.d0
   do j=1, ao_num
     tmp_matrix(j,j) = 1.d0
@@ -17,6 +19,7 @@ BEGIN_PROVIDER [double precision, ao_ortho_lowdin_coef, (ao_num_align,ao_num)]
       ao_ortho_lowdin_coef(j,i) = tmp_matrix(i,j)
     enddo
   enddo
+  deallocate(tmp_matrix)
 END_PROVIDER
 
 BEGIN_PROVIDER [double precision, ao_ortho_lowdin_overlap, (ao_num_align,ao_num)]

src/MO_Basis/ao_ortho_canonical.irp.f

@@ -82,6 +82,15 @@ END_PROVIDER
 
 
 
+BEGIN_PROVIDER [ double precision, ao_ortho_canonical_coef_inv, (ao_num,ao_num)]
+ implicit none
+ BEGIN_DOC
+! ao_ortho_canonical_coef^(-1)
+ END_DOC
+ call get_pseudo_inverse(ao_ortho_canonical_coef,ao_num,ao_num, &
+   ao_ortho_canonical_coef_inv, size(ao_ortho_canonical_coef,1))
+END_PROVIDER
+
  BEGIN_PROVIDER [ double precision, ao_ortho_canonical_coef, (ao_num_align,ao_num)]
 &BEGIN_PROVIDER [ integer, ao_ortho_canonical_num ]
   implicit none

src/MO_Basis/mos.irp.f

@@ -68,6 +68,18 @@ END_PROVIDER
   endif
 END_PROVIDER
 
+BEGIN_PROVIDER [ double precision, mo_coef_in_ao_ortho_basis, (ao_num_align, mo_tot_num) ]
+ implicit none
+ BEGIN_DOC
+ ! MO coefficients in orthogonalized AO basis
+ END_DOC
+ call dgemm('T','N',ao_num,mo_tot_num,ao_num,1.d0,                   &
+     ao_ortho_canonical_coef_inv, size(ao_ortho_canonical_coef_inv,1),&
+     mo_coef, size(mo_coef,1), 0.d0,                                 &
+     mo_coef_in_ao_ortho_basis, size(mo_coef_in_ao_ortho_basis,1))
+
+END_PROVIDER
+
 BEGIN_PROVIDER [ character*(64), mo_label ]
   implicit none
   BEGIN_DOC
@@ -257,3 +269,62 @@ subroutine mix_mo_jk(j,k)
 
 end
 
+subroutine ao_ortho_cano_to_ao(A_ao,LDA_ao,A,LDA)
+  implicit none
+  BEGIN_DOC
+  ! Transform A from the AO basis to the orthogonal AO basis
+  END_DOC
+  integer, intent(in)            :: LDA_ao,LDA
+  double precision, intent(in)   :: A_ao(LDA_ao,*)
+  double precision, intent(out)  :: A(LDA,*)
+  double precision, allocatable  :: T(:,:)
+  
+  allocate ( T(ao_num_align,ao_num) )
+  !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: T
+  
+  call dgemm('T','N', ao_num, ao_num, ao_num,                    &
+      1.d0,                                                          &
+      ao_ortho_canonical_coef_inv, size(ao_ortho_canonical_coef_inv,1),  &
+      A_ao,LDA_ao,                                                   &
+      0.d0, T, ao_num_align)
+  
+  call dgemm('N','N', ao_num, ao_num, ao_num, 1.d0,          &
+      T, ao_num_align,                                               &
+      ao_ortho_canonical_coef_inv,size(ao_ortho_canonical_coef_inv,1),&
+      0.d0, A, LDA)
+  
+  deallocate(T)
+end
+
+subroutine mo_to_ao_ortho_cano(A_mo,LDA_mo,A_ao,LDA_ao)
+  implicit none
+  BEGIN_DOC
+  ! Transform A from the AO orthogonal basis to the AO basis
+  END_DOC
+  integer, intent(in)            :: LDA_ao,LDA_mo
+  double precision, intent(in)   :: A_mo(LDA_mo)
+  double precision, intent(out)  :: A_ao(LDA_ao)
+  double precision, allocatable  :: T(:,:), SC(:,:)
+  
+  allocate ( SC(ao_num_align,mo_tot_num) )
+  allocate ( T(mo_tot_num_align,ao_num) )
+  !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: T
+  
+  call dgemm('N','N', ao_num, mo_tot_num, ao_num,                    &
+      1.d0, ao_overlap,size(ao_overlap,1),                           &
+      ao_ortho_canonical_coef, size(ao_ortho_canonical_coef,1),                                      &
+      0.d0, SC, ao_num_align)
+  
+  call dgemm('N','T', mo_tot_num, ao_num, mo_tot_num,                &
+      1.d0, A_mo,LDA_mo,                                             &
+      SC, size(SC,1),                                      &
+      0.d0, T, mo_tot_num_align)
+  
+  call dgemm('N','N', ao_num, ao_num, mo_tot_num,                    &
+      1.d0, SC,size(SC,1),                                 &
+      T, mo_tot_num_align,                                           &
+      0.d0, A_ao, LDA_ao)
+  
+  deallocate(T,SC)
+end
+

src/Utils/LinearAlgebra.irp.f

@@ -72,6 +72,8 @@ subroutine ortho_canonical(overlap,LDA,N,C,LDC,m)
   double precision, allocatable  :: S_half(:,:)
   !DEC$ ATTRIBUTES ALIGN : 64    :: U, Vt, D
   integer                        :: info, i, j
+call ortho_lowdin(overlap,LDA,N,C,LDC,m)
+return
   
   if (n < 2) then
     return
@@ -297,12 +299,12 @@ subroutine get_pseudo_inverse(A,m,n,C,LDA)
   allocate(work(lwork))
   call dgesvd('S','A', m, n, A_tmp, m,D,U,m,Vt,n,work,lwork,info)
   if (info /= 0) then
-    print *,  info, ': SVD failed'
+    print *,  info, ':: SVD failed'
     stop 1
   endif
   
   do i=1,n
-    if (abs(D(i)) > 1.d-6) then
+    if (dabs(D(i)) > 1.d-6) then
       D(i) = 1.d0/D(i)
     else
       D(i) = 0.d0