From ed0336a85b7cfc5272130d87691febe8e4193a19 Mon Sep 17 00:00:00 2001
From: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Thu, 8 Jun 2017 10:24:24 +0200
Subject: [PATCH] SlaterDressed

---
 .../LinearSystem.irp.f                        | 49 ++++++++++++++
 .../Hartree_Fock_SlaterDressed/at_nucl.irp.f  | 65 +++++++++++++++++++
 2 files changed, 114 insertions(+)
 create mode 100644 plugins/Hartree_Fock_SlaterDressed/LinearSystem.irp.f
 create mode 100644 plugins/Hartree_Fock_SlaterDressed/at_nucl.irp.f

diff --git a/plugins/Hartree_Fock_SlaterDressed/LinearSystem.irp.f b/plugins/Hartree_Fock_SlaterDressed/LinearSystem.irp.f
new file mode 100644
index 00000000..579dee57
--- /dev/null
+++ b/plugins/Hartree_Fock_SlaterDressed/LinearSystem.irp.f
@@ -0,0 +1,49 @@
+BEGIN_PROVIDER [ double precision, cusp_A, (nucl_num, nucl_num) ]
+   implicit none
+   BEGIN_DOC
+   ! Equations to solve : A.X = B
+   END_DOC
+   
+   integer                        :: mu, A, B
+
+   do B=1,nucl_num
+     do A=1,nucl_num
+        cusp_A(A,B) = 0.d0
+        if (A/=B) then
+          cusp_A(A,B) -= slater_value_at_nucl(A,B)
+        endif
+        do mu=1,ao_num
+          cusp_A(A,B) += slater_overlap(mu,B) * ao_value_at_nucl(mu,A)
+        enddo
+     enddo
+   enddo
+
+END_PROVIDER
+
+BEGIN_PROVIDER [ double precision, cusp_C, (nucl_num, mo_tot_num) ]
+   implicit none
+   BEGIN_DOC
+   ! Equations to solve : A.C = B
+   END_DOC
+   
+   integer                        :: i, A, info
+   
+   do i=1,mo_tot_num
+     do A=1,nucl_num
+       cusp_C(A,i) = mo_value_at_nucl(i,A)
+     enddo
+   enddo
+   
+   integer, allocatable           :: ipiv(:)
+   allocate ( ipiv(nucl_num) )
+   call dgegv(nucl_num, mo_tot_num, cusp_A, size(cusp_A,1),          &
+       ipiv, cusp_C, size(cusp_C,1), info)
+   deallocate (ipiv)
+   
+   if (info /= 0) then
+     print *, 'Cusp : linear solve failed'
+     stop -1
+   endif
+     
+END_PROVIDER
+
diff --git a/plugins/Hartree_Fock_SlaterDressed/at_nucl.irp.f b/plugins/Hartree_Fock_SlaterDressed/at_nucl.irp.f
new file mode 100644
index 00000000..48179628
--- /dev/null
+++ b/plugins/Hartree_Fock_SlaterDressed/at_nucl.irp.f
@@ -0,0 +1,65 @@
+BEGIN_PROVIDER [ double precision , ao_value_at_nucl, (ao_num,nucl_num) ]
+  implicit none
+  BEGIN_DOC
+! Values of the atomic orbitals at the nucleus 
+  END_DOC
+  integer :: i,j,k
+  double precision :: x,y,z,expo,poly, r2
+
+  do k=1,nucl_num
+    do i=1,ao_num
+      x = nucl_coord(ao_nucl(i),1) - nucl_coord(k,1)
+      y = nucl_coord(ao_nucl(i),2) - nucl_coord(k,2)
+      z = nucl_coord(ao_nucl(i),3) - nucl_coord(k,3)
+      poly = x**(ao_power(i,1)) * y**(ao_power(i,2)) * z**(ao_power(i,3)) 
+      if (poly == 0.d0) cycle
+
+      r2 = (x*x) + (y*y) + (z*z)
+      ao_value_at_nucl(i,k) = 0.d0
+      do j=1,ao_prim_num(i)
+        expo = ao_expo_ordered_transp(j,i)*r2
+        if (expo > 40.d0) cycle
+        ao_value_at_nucl(i,k) = ao_value_at_nucl(i,k) +             &
+            ao_coef_normalized_ordered_transp(j,i) *                &
+            dexp(-expo)
+      enddo
+      ao_value_at_nucl(i,k) *= poly
+    enddo
+  enddo
+END_PROVIDER
+
+BEGIN_PROVIDER [ double precision, mo_value_at_nucl, (mo_tot_num,nucl_num) ]
+  implicit none
+  BEGIN_DOC
+! Values of the molecular orbitals at the nucleus 
+  END_DOC
+
+  call dgemm('N','N',mo_tot_num,nucl_num,ao_num,1.d0,                &
+      mo_coef_transp, size(mo_coef_transp,1),                        &
+      ao_value_at_nucl, size(ao_value_at_nucl,1),                    &
+      0.d0, mo_value_at_nucl, size(mo_value_at_nucl,1))
+END_PROVIDER
+
+
+BEGIN_PROVIDER [ double precision , slater_value_at_nucl, (nucl_num,nucl_num) ]
+  implicit none
+  BEGIN_DOC
+! Values of the Slater orbitals (1) at the nucleus (2)
+  END_DOC
+  integer :: i,j,k
+  double precision :: x,y,z,expo,poly, r
+
+  do k=1,nucl_num
+    do i=1,nucl_num
+      x = nucl_coord(ao_nucl(i),1) - nucl_coord(k,1)
+      y = nucl_coord(ao_nucl(i),2) - nucl_coord(k,2)
+      z = nucl_coord(ao_nucl(i),3) - nucl_coord(k,3)
+
+      expo = slater_expo(i)*slater_expo(i)*((x*x) + (y*y) + (z*z))
+      if (expo > 160.d0) cycle
+      expo = dsqrt(expo)
+      slater_value_at_nucl(i,k) = dexp(-expo)
+    enddo
+  enddo
+END_PROVIDER
+