From e432c470a1ff8f5dc3c8234d7c9bdd816c9cc55e Mon Sep 17 00:00:00 2001
From: Emmanuel Giner <giner.emmanuel@gmail.com>
Date: Sun, 17 Apr 2016 23:30:04 +0200
Subject: [PATCH] begin DFT in qp

---
 plugins/DFT_Utils/EZFIO.cfg                |   4 +
 plugins/DFT_Utils/NEEDED_CHILDREN_MODULES  |   1 +
 plugins/DFT_Utils/grid_density.irp.f       | 104 ++++++++++
 plugins/DFT_Utils/integration_3d.irp.f     |  64 ++++++
 plugins/DFT_Utils/integration_radial.irp.f |  48 +++++
 plugins/DFT_Utils/routines_roland.irp.f    | 219 +++++++++++++++++++++
 src/Determinants/two_body_dm_map.irp.f     |  29 ++-
 7 files changed, 467 insertions(+), 2 deletions(-)
 create mode 100644 plugins/DFT_Utils/EZFIO.cfg
 create mode 100644 plugins/DFT_Utils/NEEDED_CHILDREN_MODULES
 create mode 100644 plugins/DFT_Utils/grid_density.irp.f
 create mode 100644 plugins/DFT_Utils/integration_3d.irp.f
 create mode 100644 plugins/DFT_Utils/integration_radial.irp.f
 create mode 100644 plugins/DFT_Utils/routines_roland.irp.f

diff --git a/plugins/DFT_Utils/EZFIO.cfg b/plugins/DFT_Utils/EZFIO.cfg
new file mode 100644
index 00000000..21cc5b98
--- /dev/null
+++ b/plugins/DFT_Utils/EZFIO.cfg
@@ -0,0 +1,4 @@
+[energy]
+type: double precision
+doc: Calculated energy
+interface: ezfio
diff --git a/plugins/DFT_Utils/NEEDED_CHILDREN_MODULES b/plugins/DFT_Utils/NEEDED_CHILDREN_MODULES
new file mode 100644
index 00000000..bff2467f
--- /dev/null
+++ b/plugins/DFT_Utils/NEEDED_CHILDREN_MODULES
@@ -0,0 +1 @@
+Determinants 
diff --git a/plugins/DFT_Utils/grid_density.irp.f b/plugins/DFT_Utils/grid_density.irp.f
new file mode 100644
index 00000000..7cb93a91
--- /dev/null
+++ b/plugins/DFT_Utils/grid_density.irp.f
@@ -0,0 +1,104 @@
+BEGIN_PROVIDER [integer, n_points_angular_grid]
+ implicit none
+ n_points_angular_grid = 18
+END_PROVIDER 
+
+BEGIN_PROVIDER [integer, n_points_radial_grid]
+ implicit none
+ n_points_radial_grid = 10
+END_PROVIDER 
+
+
+ BEGIN_PROVIDER [double precision, angular_quadrature_points, (n_points_angular_grid,3) ]
+&BEGIN_PROVIDER [double precision, weights_angular_points, (n_points_angular_grid)]
+ implicit none
+ BEGIN_DOC
+! weights and grid points for the integration on the angular variables on 
+! the unit sphere centered on (0,0,0)
+ END_DOC
+ call cal_quad(n_points_aangular_grid, angular_quadrature_points,weights_angular_points)
+
+END_PROVIDER 
+
+BEGIN_PROVIDER [double precision, grid_points_per_atom, (3,n_points_angular_grid,n_points_radial_grid,nucl_num)]
+ BEGIN_DOC
+! points for integration over space
+ END_DOC
+ implicit none
+ integer :: i,j,k
+ double precision :: dr,x_ref,y_ref,z_ref
+ dr = 1.d0/dble(n_points_radial_grid-1)
+ do i = 1, nucl_num
+  x_ref = nucl_coord(i,1)
+  y_ref = nucl_coord(i,2)
+  z_ref = nucl_coord(i,3)
+  do j = 1, n_points_radial_grid
+   do k = 1, n_points_angular_grid
+    grid_points_per_atom(1,k,j,i) = x_ref + angular_quadrature_points(k,1) * dr
+    grid_points_per_atom(2,k,j,i) = y_ref + angular_quadrature_points(k,2) * dr
+    grid_points_per_atom(3,k,j,i) = z_ref + angular_quadrature_points(k,3) * dr
+   enddo
+  enddo
+ enddo
+END_PROVIDER 
+
+BEGIN_PROVIDER [double precision, weight_functions_at_grid_points, (nucl_num,n_points_angular_grid,n_points_radial_grid) ]
+ BEGIN_DOC 
+! Weight function at grid points : w_n(r) according to the equation (22) of Becke original paper (JCP, 88, 1988)
+! the "n" discrete variable represents the nucleis (j=1,nucl_num) 
+ END_DOC
+ implicit none
+ integer :: i,j,k,l,m
+ double precision :: r(3)
+ double precision :: accu,cell_function_becke
+ double precision :: tmp_array(nucl_num)
+  do j = 1, nucl_num
+   do k = 1, n_points_radial_grid
+    do l = 1, n_points_angular_grid
+     r(1) = grid_points_per_atom(1,j,k,l)
+     r(2) = grid_points_per_atom(2,j,k,l)
+     r(3) = grid_points_per_atom(3,j,k,l)
+     accu = 0.d0
+     do i = 1, nucl_num
+      tmp_array(i) = cell_function_becke(r,i)
+      accu += tmp_array(i)
+     enddo
+     accu = 1.d0/accu
+     do i = 1, nucl_num
+      weight_functions_at_grid_points(i,j,k,l) = tmp_array(i)*accu
+     enddo
+    enddo
+   enddo
+  enddo
+
+
+END_PROVIDER 
+
+ BEGIN_PROVIDER [double precision, one_body_dm_mo_alpha_at_grid_points, (n_points_angular_grid,n_points_radial_grid,nucl_num) ]
+&BEGIN_PROVIDER [double precision, one_body_dm_mo_beta_at_grid_points, (n_points_angular_grid,n_points_radial_grid,nucl_num) ]
+ implicit none
+ integer :: i,j,k,l,m
+ double precision :: contrib
+ double precision :: r(3)
+ double precision :: aos_array(ao_num)
+  do j = 1, nucl_num
+   do k = 1, n_points_radial_grid
+    do l = 1, n_points_angular_grid
+     r(1) = grid_points_per_atom(1,j,k,l)
+     r(2) = grid_points_per_atom(2,j,k,l)
+     r(3) = grid_points_per_atom(3,j,k,l)
+     call give_all_aos_at_r(r,aos_array)
+     one_body_dm_mo_alpha_at_grid_points(j,k,l) = 0.d0
+     do i = 1, ao_num
+       do m = 1, ao_num
+        contrib = aos_array(i) * aos_array(m)
+        one_body_dm_mo_alpha_at_grid_points(j,k,l) +=  one_body_dm_ao_alpha(i,m) * contrib
+        one_body_dm_mo_beta_at_grid_points(j,k,l) +=  one_body_dm_ao_beta(i,m)  * contrib
+       enddo
+     enddo
+    enddo
+   enddo
+  enddo
+
+END_PROVIDER 
+
diff --git a/plugins/DFT_Utils/integration_3d.irp.f b/plugins/DFT_Utils/integration_3d.irp.f
new file mode 100644
index 00000000..f4088302
--- /dev/null
+++ b/plugins/DFT_Utils/integration_3d.irp.f
@@ -0,0 +1,64 @@
+double precision function step_function_becke(x)
+ implicit none
+ double precision, intent(in) :: x
+ double precision :: f_function_becke
+ integer :: i,n_max_becke
+
+ step_function_becke = f_function_becke(x)
+ n_max_becke = 2
+ do i = 1, n_max_becke
+  step_function_becke = f_function_becke(step_function_becke)
+ enddo
+ step_function_becke = 0.5d0*(1.d0 - step_function_becke)
+end
+
+double precision function f_function_becke(x)
+ implicit none
+ double precision, intent(in) :: x
+ f_function_becke = 1.5d0 * x - 0.5d0 * x*x*x
+end
+
+double precision function cell_function_becke(r,atom_number)
+ implicit none
+ double precision, intent(in) :: r(3)
+ integer, intent(in) :: atom_number
+ BEGIN_DOC
+ ! atom_number :: atom on which the cell function of Becke (1988, JCP,88(4))
+ ! r(1:3) :: x,y,z coordinantes of the current point 
+ END_DOC
+ double precision :: mu_ij,nu_ij
+ double precision :: distance_i,distance_j,step_function_becke
+ integer :: j
+ distance_i = (r(1) - nucl_coord_transp(1,atom_number) ) * (r(1) - nucl_coord_transp(1,atom_number))  
+ distance_i += (r(2) - nucl_coord_transp(2,atom_number) ) * (r(2) - nucl_coord_transp(2,atom_number))  
+ distance_i += (r(3) - nucl_coord_transp(3,atom_number) ) * (r(3) - nucl_coord_transp(3,atom_number))
+ distance_i = dsqrt(distance_i)
+ cell_function_becke = 1.d0
+ do j = 1, nucl_num
+  if(j==atom_number)cycle
+  distance_j = (r(1) - nucl_coord_transp(1,j) ) * (r(1) - nucl_coord_transp(1,j))  
+  distance_j+= (r(2) - nucl_coord_transp(2,j) ) * (r(2) - nucl_coord_transp(2,j))  
+  distance_j+= (r(3) - nucl_coord_transp(3,j) ) * (r(3) - nucl_coord_transp(3,j))
+  distance_j = dsqrt(distance_j)
+  mu_ij = (distance_i - distance_j)/nucl_dist(atom_number,j)
+  nu_ij = mu_ij + slater_bragg_type_inter_distance_ua(atom_number,j) * (1.d0 - mu_ij*mu_ij)
+  cell_function_becke *= step_function_becke(nu_ij)
+ enddo
+end
+
+double precision function weight_function_becke(r,atom_number)
+ implicit none
+ double precision, intent(in) :: r(3)
+ integer, intent(in) :: atom_number
+ BEGIN_DOC
+ ! atom_number :: atom on which the weight function of Becke (1988, JCP,88(4))
+ ! r(1:3) :: x,y,z coordinantes of the current point 
+ END_DOC
+ double precision :: cell_function_becke,accu
+ integer :: j
+ accu = 0.d0
+ do j = 1, nucl_num
+  accu += cell_function_becke(r,j)
+ enddo
+ weight_function_becke = cell_function_becke(r,atom_number)/accu
+end
diff --git a/plugins/DFT_Utils/integration_radial.irp.f b/plugins/DFT_Utils/integration_radial.irp.f
new file mode 100644
index 00000000..59874b8b
--- /dev/null
+++ b/plugins/DFT_Utils/integration_radial.irp.f
@@ -0,0 +1,48 @@
+ BEGIN_PROVIDER [ double precision, integral_density_alpha_knowles_becke_per_atom, (nucl_num)]
+&BEGIN_PROVIDER [ double precision, integral_density_beta_knowles_becke_per_atom, (nucl_num)]
+ implicit none
+ double precision :: accu
+ integer :: i,j,k,l
+ integer :: m_param_knowles
+ double precision :: dx,x
+ integer :: n_pt_int_radial
+ double precision :: integrand(n_points_angular_grid), weights(n_points_angular_grid)
+ double precision :: f_average_angular_alpha,f_average_angular_beta
+ double precision :: derivative_knowles_function,knowles_function
+ n_pt_int_radial = 10
+ dx = 1.d0/dble(n_pt_int_radial-1)
+ x = 0.d0
+ m_param_knowles = 3
+ do j = 1, nucl_num
+  integral_density_alpha_knowles_becke_per_atom(j) = 0.d0
+  do i = 1, n_points_radial_grid
+   ! Angular integration 
+   f_average_angular_alpha = 0.d0
+   f_average_angular_beta = 0.d0
+   do k = 1, n_points_angular_grid
+    f_average_angular_alpha += weights_angular_points(k) * one_body_dm_mo_alpha_at_grid_points(k,i,j) * weight_functions_at_grid_points(k,i,j)
+    f_average_angular_beta += weights_angular_points(k) * one_body_dm_mo_beta_at_grid_points(k,i,j) * weight_functions_at_grid_points(k,i,j)
+   enddo
+   integral_density_alpha_knowles_becke_per_atom(j) += derivative_knowles_function(alpha,m_param_knowles,x) & 
+                                                      *knowles_function(alpha,m_param_knowles,x)**2         &
+                                                      *f_average_angular_alpha
+   x += dx
+  enddo
+  integral_density_alpha_knowles_becke_per_atom(j) *= dx
+ enddo
+ 
+END_PROVIDER 
+
+ double precision function knowles_function(alpha,m,x)
+ implicit none
+ double precision, intent(in) :: alpha,x
+ integer, intent(in) :: m
+ knowles_function = -alpha * dlog(1.d0-x**m)
+ end
+
+ double precision function derivative_knowles_function(alpha,m,x)
+ implicit none
+ double precision, intent(in) :: alpha,x
+ integer, intent(in) :: m
+ derivative_knowles_function = m x**(m-1) / (alpha * dlog(1.d0-x**m))
+ end
diff --git a/plugins/DFT_Utils/routines_roland.irp.f b/plugins/DFT_Utils/routines_roland.irp.f
new file mode 100644
index 00000000..0f555902
--- /dev/null
+++ b/plugins/DFT_Utils/routines_roland.irp.f
@@ -0,0 +1,219 @@
+
+  subroutine cal_quad(n_quad, quad, weight)
+! --------------------------------------------------------------------------------
+! 
+! Arguments  : subroutine cal_quad
+! Description: evaluates quadrature points an weights
+! 
+! Authors    : B. Lévy, P. Pernot
+! Date       : 15 Nov 2000
+! --------------------------------------------------------------------------------
+  implicit none
+  integer, intent(in) :: n_quad
+  double precision, intent(out) :: weight(n_quad)
+  double precision, intent(out) :: quad(n_quad,3)
+  
+! local:
+    double precision, parameter :: zero=0.d0, one= 1.d0
+
+    double precision, parameter :: p=0.707106781186547462d0
+    double precision, parameter :: q=0.577350269189625842d0
+    double precision, parameter :: r=0.301511344577763629d0
+    double precision, parameter :: s=0.904534033733290888d0
+
+    double precision, parameter :: fourpi= 12.5663706143591725d0
+
+    double precision, parameter :: a6=0.166666666666666657d0
+    double precision, parameter :: a18=0.333333333333333329d-01
+    double precision, parameter :: b18=0.666666666666666657d-01
+    double precision, parameter :: a26=0.476190476190476164d-01
+    double precision, parameter :: b26=0.380952380952380987d-01
+    double precision, parameter :: c26=0.321428571428571397d-01
+    double precision, parameter :: a50=0.126984126984126984d-01
+    double precision, parameter :: b50=0.225749559082892431d-01
+    double precision, parameter :: c50=0.210937500000000014d-01
+    double precision, parameter :: d50=0.201733355379188697d-01
+
+    double precision            :: apt(3,6),bpt(3,12),cpt(3,8),dpt(3,24)
+    double precision            :: awght,bwght,cwght,dwght
+    double precision            :: s1, s2, s3
+    integer             :: idim, ipt, i1, i2, i3, is1, is2, is3
+    integer             :: iquad
+
+! begin:
+! l_here ='cal_quad'
+! call enter (l_here,3)
+
+! verifications:
+!   message = 'in '//trim(l_here)//', number of dimensions='//&
+!             trim(encode(dimensions_nb))//', must be 3'
+!   call ensure(message, dimensions_nb .eq. 3 )
+
+!   message = 'in '//trim(l_here)//', invalid number of quadrature points ='&
+!             //trim(encode(n_quad))
+!   call ensure(message,(n_quad-2)*(n_quad-6)*(n_quad-18)*(n_quad-26)*(n_quad-50) .eq. 0)
+
+! initialize weights
+    awght = zero
+    bwght = zero
+    cwght = zero
+    dwght = zero
+
+! type A points : (+/-1,0,0)
+    awght=a6*fourpi
+    ipt= 1
+    apt=0.
+    do idim = 1, 3
+      apt(idim,ipt)=one
+      ipt=ipt+1
+      apt(idim,ipt)=-one
+      ipt=ipt+1
+    enddo     
+
+! type B points : (+/-p,+/-p,0) with p= 1/sqrt(2)
+    if(n_quad.gt.6) then
+
+      awght=a18*fourpi
+      bwght=b18*fourpi
+
+      s1=p
+      s2=p
+      ipt= 1
+      bpt=0.
+      do idim = 1, 3
+        i1=idim+1
+        if(i1.gt.3) i1=i1-3
+        i2=idim+2
+        if(i2.gt.3) i2=i2-3
+        do is1= 1,2
+          do is2= 1,2
+            bpt(i1,ipt)=s1
+            bpt(i2,ipt)=s2
+            s2=-s2
+            ipt=ipt+1
+          enddo      
+          s1=-s1
+        enddo       
+      enddo    
+    endif
+
+! type C points : (+/-q,+/-q,+/-q) with q= 1/sqrt(3)
+   if(n_quad.gt.18) then
+
+      awght=a26*fourpi
+      bwght=b26*fourpi
+      cwght=c26*fourpi
+
+      s1=q
+      s2=q
+      s3=q
+      ipt= 1
+      cpt=0.
+      do is1= 1,2
+        do is2= 1,2
+          do is3= 1,2
+            cpt(1,ipt)=s1
+            cpt(2,ipt)=s2
+            cpt(3,ipt)=s3
+            s3=-s3
+            ipt=ipt+1
+          enddo       
+          s2=-s2
+        enddo       
+        s1=-s1
+      enddo    
+    endif
+
+! type D points : (+/-r,+/-r,+/-s)
+    if(n_quad.gt.26) then
+
+      awght=a50*fourpi
+      bwght=b50*fourpi
+      cwght=c50*fourpi
+      dwght=d50*fourpi
+
+      ipt= 1
+      dpt=0.
+      do i1= 1, 3
+        s1=s
+        s2=r
+        s3=r
+        i2=i1+1
+        if(i2.gt.3) i2=i2-3
+        i3=i1+2
+        if(i3.gt.3) i3=i3-3
+        do is1= 1,2
+          do is2= 1,2
+            do is3= 1,2
+              dpt(i1,ipt)=s1
+              dpt(i2,ipt)=s2
+              dpt(i3,ipt)=s3
+              s3=-s3
+              ipt=ipt+1
+            enddo         
+            s2=-s2
+          enddo        
+          s1=-s1
+        enddo         
+      enddo     
+    endif
+
+! fill the points and weights tables
+    iquad= 1
+    do ipt= 1, 6
+      do idim = 1, 3
+        quad(iquad,idim)=apt(idim,ipt)
+      enddo
+      weight(iquad)=awght
+      iquad=iquad+1
+    enddo     
+
+    if(n_quad.gt.6) then
+      do ipt= 1,12
+        do idim = 1, 3
+          quad(iquad,idim)=bpt(idim,ipt)
+        enddo
+        weight(iquad)=bwght
+        iquad=iquad+1
+      enddo            
+    endif
+
+    if(n_quad.gt.18) then
+      do ipt= 1,8
+        do idim = 1, 3
+          quad(iquad,idim)=cpt(idim,ipt)
+        enddo
+        weight(iquad)=cwght
+        iquad=iquad+1
+      enddo            
+    endif
+
+    if(n_quad.gt.26) then
+      do ipt= 1,24
+        do idim = 1, 3
+          quad(iquad,idim)=dpt(idim,ipt)
+        enddo
+        weight(iquad)=dwght
+        iquad=iquad+1
+      enddo            
+    endif
+
+!   if (debug) then
+!     write(6,*)
+!     write(6,'(1X,a)') trim(l_here)//'-d : '//&
+!       '------------------------------------------------------'
+!     write(6,'(1X,a)') trim(l_here)//'-d : '//'  I    Weight      Quad_points'
+!     write(6,'(1X,a)') trim(l_here)//'-d : '//&
+!       '-----  ----------  -----------------------------------'
+!     do iquad= 1, n_quad
+!        write(6,'(1X,A,i5,4e12.3)') trim(l_here)//'-d : ',&
+!          iquad,weight(iquad),quad(iquad,1:3)
+!     enddo
+!     write(6,'(1X,a)') trim(l_here)//'-d : '//&
+!       '------------------------------------------------------'
+!     write(6,*)
+!   endif
+
+! call exit (l_here,3)
+
+  end subroutine cal_quad
diff --git a/src/Determinants/two_body_dm_map.irp.f b/src/Determinants/two_body_dm_map.irp.f
index f570e2bf..38467040 100644
--- a/src/Determinants/two_body_dm_map.irp.f
+++ b/src/Determinants/two_body_dm_map.irp.f
@@ -268,8 +268,6 @@ BEGIN_PROVIDER [double precision, two_body_dm_ab_big_array, (n_act_orb,n_act_orb
  enddo
  print*,'Big array for density matrix provided !'
 
-
-
 END_PROVIDER 
 
 subroutine insert_into_two_body_dm_big_array(big_array,dim1,dim2,dim3,dim4,contrib,h1,p1,h2,p2)
@@ -291,3 +289,30 @@ subroutine insert_into_two_body_dm_big_array(big_array,dim1,dim2,dim3,dim4,contr
 !endif
  
 end
+
+double precision function compute_two_body_dm_ab(r1,r2)
+ implicit none
+ double precision :: r1(3),r2(3)
+ integer :: i,j,k,l
+ double precision :: mos_array_r1(mo_tot_num),mos_array_r2(mo_tot_num)
+ double precision :: contrib
+ compute_two_body_dm_ab = 0.d0
+ call give_all_mos_at_r(r1,mos_array_r1)
+ call give_all_mos_at_r(r2,mos_array_r2)
+ do l = 1, n_act_orb  ! p2 
+  contrib = mos_array_r2(l)
+  if(dabs(contrib).lt.1.d-6)cycle
+  do k = 1, n_act_orb  ! h2 
+    contrib *= mos_array_r2(k)
+   if(dabs(contrib).lt.1.d-6)cycle
+   do j = 1, n_act_orb  ! p1 
+    contrib *= mos_array_r1(j)
+    if(dabs(contrib).lt.1.d-6)cycle
+    do i = 1,n_act_orb   ! h1 
+     compute_two_body_dm_ab += two_body_dm_ab_big_array(i,j,k,l) * mos_array_r1(i) * contrib
+    enddo
+   enddo
+  enddo
+ enddo
+
+end