Merge pull request #107 from QuantumPackage/cleaning_dft

Cleaning dft

src/becke_numerical_grid/EZFIO.cfg

@@ -14,3 +14,22 @@ type: double precision
 doc: threshold on the weight of a given grid point
 interface: ezfio,provider,ocaml
 default: 1.e-20
+
+[my_grid_becke]
+type: logical
+doc: if True, the number of angular and radial grid points are read from EZFIO
+interface: ezfio,provider,ocaml
+default: False
+
+[my_n_pt_r_grid]
+type: integer
+doc: Number of radial grid points given from input
+interface: ezfio,provider,ocaml
+default: 300
+
+[my_n_pt_a_grid]
+type: integer
+doc: Number of angular grid points given from input. Warning, this number cannot be any integer. See file list_angular_grid
+interface: ezfio,provider,ocaml
+default: 1202
+

src/becke_numerical_grid/grid_becke.irp.f

@@ -8,6 +8,7 @@
  !
  ! These numbers are automatically set by setting the grid_type_sgn parameter
  END_DOC
+if(.not.my_grid_becke)then
  select case (grid_type_sgn)
     case(0)
      n_points_radial_grid = 23
@@ -25,6 +26,10 @@ select case (grid_type_sgn)
       write(*,*) '!!! Quadrature grid not available !!!'
       stop
   end select
+else
+ n_points_radial_grid = my_n_pt_r_grid
+ n_points_integration_angular = my_n_pt_a_grid
+endif
 END_PROVIDER
 
 BEGIN_PROVIDER [integer, n_points_grid_per_atom]

src/becke_numerical_grid/list_angular_grid

@@ -0,0 +1,32 @@
+0006
+0014
+0026
+0038
+0050
+0074
+0086
+0110
+0146
+0170
+0194
+0230
+0266
+0302
+0350
+0434
+0590
+0770
+0974
+1202
+1454
+1730
+2030
+2354
+2702
+3074
+3470
+3890
+4334
+4802
+5294
+5810

src/cas_based_on_top/two_body_dens_rout.irp.f

@@ -0,0 +1,152 @@
+
+subroutine give_n2_ii_val_ab(r1,r2,two_bod_dens)
+ implicit none
+ BEGIN_DOC
+! contribution from purely inactive orbitals to n2_{\Psi^B}(r_1,r_2) for a CAS wave function 
+ END_DOC
+ double precision, intent(in) :: r1(3),r2(3)
+ double precision, intent(out):: two_bod_dens
+ integer :: i,j,m,n,i_m,i_n
+ integer :: i_i,i_j
+ double precision, allocatable :: mos_array_inact_r1(:),mos_array_inact_r2(:)
+ double precision, allocatable :: mos_array_r1(:) , mos_array_r2(:) 
+ ! You get all orbitals in r_1 and r_2
+ allocate(mos_array_r1(mo_num) , mos_array_r2(mo_num) )
+ call give_all_mos_at_r(r1,mos_array_r1) 
+ call give_all_mos_at_r(r2,mos_array_r2) 
+ ! You extract the inactive orbitals 
+ allocate(mos_array_inact_r1(n_inact_orb) , mos_array_inact_r2(n_inact_orb) )
+ do i_m = 1, n_inact_orb
+  mos_array_inact_r1(i_m) = mos_array_r1(list_inact(i_m))
+ enddo
+ do i_m = 1, n_inact_orb
+  mos_array_inact_r2(i_m) = mos_array_r2(list_inact(i_m))
+ enddo
+
+ two_bod_dens = 0.d0
+ ! You browse all OCCUPIED ALPHA electrons in the \mathcal{A} space
+ do m = 1, n_inact_orb ! electron 1 
+ ! You browse all OCCUPIED BETA  electrons in the \mathcal{A} space
+  do n = 1, n_inact_orb ! electron 2 
+   ! two_bod_dens(r_1,r_2) = n_alpha(r_1) * n_beta(r_2)
+   two_bod_dens += mos_array_inact_r1(m) * mos_array_inact_r1(m) * mos_array_inact_r2(n) * mos_array_inact_r2(n) 
+  enddo
+ enddo
+end
+
+
+subroutine give_n2_ia_val_ab(r1,r2,two_bod_dens,istate)
+ BEGIN_DOC
+! contribution from inactive and active orbitals to n2_{\Psi^B}(r_1,r_2) for the "istate" state of a CAS wave function 
+ END_DOC
+ implicit none
+ integer, intent(in) :: istate
+ double precision, intent(in) :: r1(3),r2(3)
+ double precision, intent(out):: two_bod_dens
+ integer :: i,orb_i,a,orb_a,n,m,b
+ double precision :: rho
+ double precision, allocatable :: mos_array_r1(:) , mos_array_r2(:) 
+ double precision, allocatable :: mos_array_inact_r1(:),mos_array_inact_r2(:)
+ double precision, allocatable :: mos_array_act_r1(:),mos_array_act_r2(:)
+
+ two_bod_dens = 0.d0
+ ! You get all orbitals in r_1 and r_2
+ allocate(mos_array_r1(mo_num) , mos_array_r2(mo_num) )
+ call give_all_mos_at_r(r1,mos_array_r1) 
+ call give_all_mos_at_r(r2,mos_array_r2) 
+
+ ! You extract the inactive orbitals 
+ allocate( mos_array_inact_r1(n_inact_orb) , mos_array_inact_r2(n_inact_orb) )
+ do i = 1, n_inact_orb
+  mos_array_inact_r1(i) = mos_array_r1(list_inact(i))
+ enddo
+ do i= 1, n_inact_orb
+  mos_array_inact_r2(i) = mos_array_r2(list_inact(i))
+ enddo
+
+ ! You extract the active orbitals 
+ allocate( mos_array_act_r1(n_act_orb) , mos_array_act_r2(n_act_orb) )
+ do i= 1, n_act_orb
+  mos_array_act_r1(i) = mos_array_r1(list_act(i))
+ enddo
+ do i= 1, n_act_orb
+  mos_array_act_r2(i) = mos_array_r2(list_act(i))
+ enddo
+
+ ! Contracted density : intermediate quantity
+ two_bod_dens = 0.d0
+ do a = 1, n_act_orb 
+  do i = 1, n_inact_orb
+   do b = 1, n_act_orb 
+    rho = one_e_act_dm_beta_mo_for_dft(b,a,istate) + one_e_act_dm_alpha_mo_for_dft(b,a,istate) 
+    two_bod_dens += mos_array_inact_r1(i) * mos_array_inact_r1(i) * mos_array_act_r2(a) * mos_array_act_r2(b) * rho
+   enddo
+  enddo
+ enddo
+end
+
+
+subroutine give_n2_aa_val_ab(r1,r2,two_bod_dens,istate)
+ BEGIN_DOC
+! contribution from purely active orbitals to n2_{\Psi^B}(r_1,r_2) for the "istate" state of a CAS wave function 
+ END_DOC
+ implicit none
+ integer, intent(in) :: istate
+ double precision, intent(in) :: r1(3),r2(3)
+ double precision, intent(out):: two_bod_dens
+ integer :: i,orb_i,a,orb_a,n,m,b,c,d
+ double precision :: rho
+ double precision, allocatable :: mos_array_r1(:) , mos_array_r2(:) 
+ double precision, allocatable :: mos_array_act_r1(:),mos_array_act_r2(:)
+
+ two_bod_dens = 0.d0
+ ! You get all orbitals in r_1 and r_2
+ allocate(mos_array_r1(mo_num) , mos_array_r2(mo_num) )
+ call give_all_mos_at_r(r1,mos_array_r1) 
+ call give_all_mos_at_r(r2,mos_array_r2) 
+
+ ! You extract the active orbitals 
+ allocate( mos_array_act_r1(n_act_orb) , mos_array_act_r2(n_act_orb) )
+ do i= 1, n_act_orb
+  mos_array_act_r1(i) = mos_array_r1(list_act(i))
+ enddo
+ do i= 1, n_act_orb
+  mos_array_act_r2(i) = mos_array_r2(list_act(i))
+ enddo
+
+ ! Contracted density : intermediate quantity
+ two_bod_dens = 0.d0
+ do a = 1, n_act_orb  ! 1
+  do b = 1, n_act_orb  ! 2
+   do c = 1, n_act_orb  ! 1
+    do d = 1, n_act_orb  ! 2
+     rho = mos_array_act_r1(c) * mos_array_act_r2(d) * act_2_rdm_ab_mo(d,c,b,a,istate)
+     two_bod_dens += rho * mos_array_act_r1(a) * mos_array_act_r2(b) 
+    enddo
+   enddo
+  enddo
+ enddo
+
+end
+
+subroutine give_n2_cas(r1,r2,istate,n2_psi)
+ implicit none
+ BEGIN_DOC
+! returns mu(r), f_psi, n2_psi for a general cas wave function
+ END_DOC
+ integer, intent(in) :: istate
+ double precision, intent(in)  :: r1(3),r2(3)
+ double precision, intent(out) :: n2_psi
+ double precision :: two_bod_dens_ii
+ double precision :: two_bod_dens_ia
+ double precision :: two_bod_dens_aa
+ ! inactive-inactive part of n2_psi(r1,r2)
+ call give_n2_ii_val_ab(r1,r2,two_bod_dens_ii)
+ ! inactive-active part of n2_psi(r1,r2)
+ call give_n2_ia_val_ab(r1,r2,two_bod_dens_ia,istate)
+ ! active-active part of n2_psi(r1,r2)
+ call give_n2_aa_val_ab(r1,r2,two_bod_dens_aa,istate)
+
+ n2_psi = two_bod_dens_ii + two_bod_dens_ia + two_bod_dens_aa
+ 
+end

src/casscf/casscf.irp.f

@@ -5,6 +5,7 @@ program casscf
   END_DOC
   call reorder_orbitals_for_casscf
   no_vvvv_integrals = .True.
+  touch no_vvvv_integrals
   pt2_max = 0.02
   SOFT_TOUCH no_vvvv_integrals pt2_max
   call run_stochastic_cipsi

src/mo_two_e_ints/EZFIO.cfg

@@ -11,3 +11,9 @@ interface: ezfio,provider,ocaml
 default: 1.e-15
 ezfio_name: threshold_mo
 
+[no_vvvv_integrals]
+type: logical
+doc: If `True`, computes all integrals except for the integrals having 3 or 4 virtual indices
+interface: ezfio,provider,ocaml
+default: false
+

src/mo_two_e_ints/four_idx_novvvv.irp.f

@@ -1,11 +1,11 @@
-BEGIN_PROVIDER [ logical, no_vvvv_integrals  ]
+!BEGIN_PROVIDER [ logical, no_vvvv_integrals  ]
-  implicit none
+!  implicit none
-  BEGIN_DOC
+!  BEGIN_DOC
 ! If `True`, computes all integrals except for the integrals having 3 or 4 virtual indices
-  END_DOC
+!  END_DOC
-
+!
-  no_vvvv_integrals = .False.
+!  no_vvvv_integrals = .False.
-END_PROVIDER
+!END_PROVIDER
 
 BEGIN_PROVIDER [ double precision, mo_coef_novirt, (ao_num,n_core_inact_act_orb) ]
  implicit none
@@ -56,6 +56,8 @@ subroutine four_idx_novvvv
   BEGIN_DOC
   ! Retransform MO integrals for next CAS-SCF step
   END_DOC
+  print*,'Using partial transformation'
+  print*,'It will not transform all integrals with at least 3 indices within the virtuals'
   integer                        :: i,j,k,l,n_integrals
   double precision, allocatable  :: f(:,:,:), f2(:,:,:), d(:,:), T(:,:,:,:), T2(:,:,:,:)
   double precision, external     :: get_ao_two_e_integral

src/two_body_rdm/act_2_rdm.irp.f

@@ -2,6 +2,8 @@
  BEGIN_PROVIDER [double precision, act_2_rdm_ab_mo, (n_act_orb,n_act_orb,n_act_orb,n_act_orb,N_states)]
  implicit none
  BEGIN_DOC
+!                             12 12
+!                 1 2 1 2 == <ij|kl>
 ! act_2_rdm_ab_mo(i,j,k,l,istate) =  STATE SPECIFIC physicist notation for 2RDM of alpha/beta electrons 
 ! 
 ! <Psi_{istate}| a^{\dagger}_{i \alpha} a^{\dagger}_{j \beta} a_{l \beta} a_{k \alpha} |Psi_{istate}>

src/utils/integration.irp.f

@@ -75,7 +75,6 @@ subroutine give_explicit_poly_and_gaussian(P_new,P_center,p,fact_k,iorder,alpha,
   P_new(0,1) = 0.d0
   P_new(0,2) = 0.d0
   P_new(0,3) = 0.d0
-
   !DIR$ FORCEINLINE
   call gaussian_product(alpha,A_center,beta,B_center,fact_k,p,P_center)
   if (fact_k < thresh) then

src/utils/shank.irp.f

@@ -0,0 +1,54 @@
+double precision function shank_general(array,n,nmax)
+ implicit none                                                                                                                                                         
+ integer, intent(in) :: n,nmax
+ double precision, intent(in) :: array(0:nmax) ! array of the partial sums
+ integer :: ntmp,i
+ double precision :: sum(0:nmax),shank1(0:nmax)
+ if(n.lt.3)then
+  print*,'You asked to Shank a sum but the order is smaller than 3 ...'
+  print*,'n = ',n
+  print*,'stopping ....'
+  stop
+ endif
+ ntmp = n
+ sum = array
+ i = 0
+ do while(ntmp.ge.2)
+  i += 1
+!  print*,'i = ',i
+  call shank(sum,ntmp,nmax,shank1)
+  ntmp = ntmp - 2
+  sum = shank1
+  shank_general = shank1(ntmp)
+ enddo
+end
+
+
+subroutine shank(array,n,nmax,shank1)
+ implicit none
+ integer, intent(in) :: n,nmax
+ double precision, intent(in)  :: array(0:nmax)
+ double precision, intent(out) :: shank1(0:nmax)
+ integer :: i,j
+ double precision :: shank_function
+ do i = 1, n-1
+  shank1(i-1) = shank_function(array,i,nmax)
+ enddo
+end
+
+double precision function shank_function(array,i,n)
+ implicit none
+ integer, intent(in) :: i,n
+ double precision, intent(in) :: array(0:n)
+ double precision :: b_n, b_n1
+ b_n = array(i) - array(i-1)
+ b_n1 = array(i+1) - array(i)
+ if(dabs(b_n1-b_n).lt.1.d-12)then
+  shank_function = array(i+1)
+ else
+  shank_function = array(i+1) - b_n1*b_n1/(b_n1-b_n)
+ endif
+
+end
+
+