Merge pull request #189 from QuantumPackage/cleaning_dft

frozen core unset
2024-06-02 03:15:34 +02:00 · 2022-03-07 16:41:35 +01:00 · 2022-03-07 16:41:35 +01:00 · 3d63d39df8
commit 3d63d39df8
parent 95638079d7 980b48d906
8 changed files with 64 additions and 47 deletions
--- a/bin/qp_set_frozen_core
+++ b/bin/qp_set_frozen_core
@ -7,12 +7,13 @@ setting all MOs as Active, except the n/2 first ones which are set as Core.
 If pseudo-potentials are used, all the MOs are set as Active.

 Usage:
-      qp_set_frozen_core [-q|--query] [(-l|-s|--large|--small)] EZFIO_DIR
+      qp_set_frozen_core [-q|--query] [(-l|-s|-u|--large|--small|--unset)] EZFIO_DIR

 Options:
    -q --query   Prints in the standard output the number of frozen MOs
    -l --large   Use a small core
    -s --small   Use a large core
+    -u --unset   Unset frozen core


 Default numbers of frozen electrons:
@ -88,7 +89,9 @@ def main(arguments):
                elif charge <=  54: n_frozen += 9
                elif charge <=  86: n_frozen += 18
                elif charge <= 118: n_frozen += 27
+        elif arguments["--unset"]:

+              n_frozen = 0
        else:  # default                                                                      
            for charge in ezfio.nuclei_nucl_charge:
                if   charge <=   4: pass
--- a/src/csf/configurations.irp.f
+++ b/src/csf/configurations.irp.f
@ -779,6 +779,7 @@ subroutine binary_search_cfg(cfgInp,addcfg)
 end subroutine

 BEGIN_PROVIDER [ integer, psi_configuration_to_psi_det, (2,N_configuration) ]
+&BEGIN_PROVIDER [ integer, psi_configuration_n_det, (N_configuration) ]
 &BEGIN_PROVIDER [ integer, psi_configuration_to_psi_det_data, (N_det) ]

 implicit none
@ -867,6 +868,29 @@ end subroutine
 enddo

 deallocate(dets, old_order)
+ integer :: ndet_conf
+ do i = 1, N_configuration
+  ndet_conf = psi_configuration_to_psi_det(2,i) - psi_configuration_to_psi_det(1,i) + 1
+  psi_configuration_n_det(i) = ndet_conf
+ enddo

 END_PROVIDER

+
+BEGIN_PROVIDER [ integer, n_elec_alpha_for_psi_configuration, (N_configuration)]
+ implicit none
+ integer :: i,j,k,l
+ integer(bit_kind) :: det_tmp(N_int,2),det_alpha(N_int)
+ n_elec_alpha_for_psi_configuration = 0
+ do i = 1, N_configuration
+  j = psi_configuration_to_psi_det(2,i) 
+  det_tmp(:,:) = psi_det(:,:,j)
+  k = 0
+  do l = 1, N_int
+   det_alpha(N_int) = iand(det_tmp(l,1),psi_configuration(l,1,i))
+   k += popcnt(det_alpha(l))
+  enddo
+  n_elec_alpha_for_psi_configuration(i) = k
+ enddo
+
+END_PROVIDER 
--- a/src/dav_general_mat/dav_diag_dressed_ext_rout.irp.f
+++ b/src/dav_general_mat/dav_diag_dressed_ext_rout.irp.f
@ -1,5 +1,5 @@

-subroutine davidson_general_ext_rout(u_in,H_jj,Dress_jj,energies,sze,N_st,N_st_diag_in,converged,hcalc)
+subroutine davidson_general_ext_rout_diag_dressed(u_in,H_jj,Dress_jj,energies,sze,N_st,N_st_diag_in,converged,hcalc)
  use mmap_module
  implicit none
  BEGIN_DOC
@ -412,36 +412,6 @@ subroutine davidson_general_ext_rout(u_in,H_jj,Dress_jj,energies,sze,N_st,N_st_d
  FREE nthreads_davidson
 end

-subroutine hcalc_template(v,u,N_st,sze)
-  use bitmasks
-  implicit none
-  BEGIN_DOC
-  ! Template of routine for the application of H
-  !
-  ! Here, it is done with the Hamiltonian matrix 
-  !
-  ! on the set of determinants of psi_det 
-  !
-  ! Computes $v = H | u \rangle$ 
-  !
-  END_DOC
-  integer, intent(in)              :: N_st,sze
-  double precision, intent(in)     :: u(sze,N_st)
-  double precision, intent(inout)  :: v(sze,N_st)
-  integer :: i,j,istate
-  v = 0.d0
-  do istate = 1, N_st
-   do i = 1, sze
-    do j = 1, sze
-      v(i,istate) += H_matrix_all_dets(j,i) * u(j,istate)
-    enddo
-   enddo
-   do i = 1, sze
-    v(i,istate) += u(i,istate) * nuclear_repulsion
-   enddo
-  enddo
-end
-
 subroutine dressing_diag_uv(v,u,dress_diag,N_st,sze)
  implicit none
  BEGIN_DOC
--- a/src/determinants/s2.irp.f
+++ b/src/determinants/s2.irp.f
@ -103,13 +103,17 @@ BEGIN_PROVIDER [ double precision, expected_s2]

 END_PROVIDER

-BEGIN_PROVIDER [ double precision, s2_values, (N_states) ]
+ BEGIN_PROVIDER [ double precision, s2_values, (N_states) ]
+&BEGIN_PROVIDER [ double precision, s_values, (N_states) ]
 implicit none
 BEGIN_DOC
 ! array of the averaged values of the S^2 operator on the various states
 END_DOC
 integer :: i
 call u_0_S2_u_0(s2_values,psi_coef,n_det,psi_det,N_int,N_states,psi_det_size)
+ do i = 1, N_states
+  s_values(i) = 0.5d0 *(-1.d0 + dsqrt(1.d0 + 4 * s2_values(i)))
+ enddo

 END_PROVIDER

--- a/src/dft_utils_in_r/ao_in_r.irp.f
+++ b/src/dft_utils_in_r/ao_in_r.irp.f
@ -91,7 +91,19 @@
 enddo
 END_PROVIDER

- BEGIN_PROVIDER[double precision, aos_lapl_in_r_array, (ao_num,n_points_final_grid,3)]
+ BEGIN_PROVIDER [double precision, aos_lapl_in_r_array_transp, (ao_num, n_points_final_grid,3)]
+ implicit none
+ integer :: i,j,m
+ do i = 1, n_points_final_grid
+  do j = 1, ao_num
+   do m = 1, 3
+    aos_lapl_in_r_array_transp(j,i,m) =  aos_lapl_in_r_array(m,j,i)
+   enddo
+  enddo
+ enddo
+ END_PROVIDER 
+
+ BEGIN_PROVIDER [double precision, aos_lapl_in_r_array, (3,ao_num,n_points_final_grid)]
 implicit none
 BEGIN_DOC
 ! aos_lapl_in_r_array(i,j,k)   = value of the kth component of the laplacian of jth ao on the ith grid point
@ -100,20 +112,20 @@
 END_DOC
 integer :: i,j,m
 double precision :: aos_array(ao_num), r(3)
- double precision :: aos_grad_array(ao_num,3)
- double precision :: aos_lapl_array(ao_num,3)
+ double precision :: aos_grad_array(3,ao_num)
+ double precision :: aos_lapl_array(3,ao_num)
 !$OMP PARALLEL DO &
 !$OMP DEFAULT (NONE)  &
 !$OMP PRIVATE (i,r,aos_array,aos_grad_array,aos_lapl_array,j,m) & 
 !$OMP SHARED(aos_lapl_in_r_array,n_points_final_grid,ao_num,final_grid_points)
- do m = 1, 3
-  do i = 1, n_points_final_grid
-   r(1) = final_grid_points(1,i)
-   r(2) = final_grid_points(2,i)
-   r(3) = final_grid_points(3,i)
-   call give_all_aos_and_grad_and_lapl_at_r(r,aos_array,aos_grad_array,aos_lapl_array)
-   do j = 1, ao_num
-    aos_lapl_in_r_array(j,i,m) = aos_lapl_array(j,m)
+ do i = 1, n_points_final_grid
+  r(1) = final_grid_points(1,i)
+  r(2) = final_grid_points(2,i)
+  r(3) = final_grid_points(3,i)
+  call give_all_aos_and_grad_and_lapl_at_r(r,aos_array,aos_grad_array,aos_lapl_array)
+  do j = 1, ao_num
+   do m = 1, 3
+    aos_lapl_in_r_array(m,j,i) = aos_lapl_array(m,j)
   enddo
  enddo
 enddo
--- a/src/dft_utils_in_r/mo_in_r.irp.f
+++ b/src/dft_utils_in_r/mo_in_r.irp.f
@ -138,7 +138,7 @@
 integer :: m
 mos_lapl_in_r_array = 0.d0
 do m=1,3
-  call dgemm('N','N',mo_num,n_points_final_grid,ao_num,1.d0,mo_coef_transp,mo_num,aos_lapl_in_r_array(1,1,m),ao_num,0.d0,mos_lapl_in_r_array(1,1,m),mo_num)
+  call dgemm('N','N',mo_num,n_points_final_grid,ao_num,1.d0,mo_coef_transp,mo_num,aos_lapl_in_r_array_transp(1,1,m),ao_num,0.d0,mos_lapl_in_r_array(1,1,m),mo_num)
 enddo
 END_PROVIDER

--- a/src/mo_guess/h_core_guess_routine.irp.f
+++ b/src/mo_guess/h_core_guess_routine.irp.f
@ -7,7 +7,7 @@ subroutine hcore_guess
  label = 'Guess'
  call mo_as_eigvectors_of_mo_matrix(mo_one_e_integrals,          &
                                     size(mo_one_e_integrals,1),  &
-                                     size(mo_one_e_integrals,2),label,1,.false.)
+                                     size(mo_one_e_integrals,2),label,1,.true.)
  call nullify_small_elements(ao_num, mo_num, mo_coef, size(mo_coef,1), 1.d-12 )
  call save_mos
  TOUCH mo_coef mo_label
--- a/src/tools/print_wf.irp.f
+++ b/src/tools/print_wf.irp.f
@ -32,8 +32,9 @@ subroutine routine
 double precision :: norm_mono_a,norm_mono_b
 double precision :: norm_mono_a_2,norm_mono_b_2
 double precision :: norm_mono_a_pert_2,norm_mono_b_pert_2
- double precision :: norm_mono_a_pert,norm_mono_b_pert
+ double precision :: norm_mono_a_pert,norm_mono_b_pert,norm_double_1
 double precision :: delta_e,coef_2_2
+
 norm_mono_a = 0.d0
 norm_mono_b = 0.d0
 norm_mono_a_2 = 0.d0
@ -42,6 +43,7 @@ subroutine routine
 norm_mono_b_pert = 0.d0
 norm_mono_a_pert_2 = 0.d0
 norm_mono_b_pert_2 = 0.d0
+ norm_double_1 = 0.d0
 do i = 1, min(N_det_print_wf,N_det)
  print*,''
  print*,'i = ',i
@ -93,6 +95,7 @@ subroutine routine
    print*,'h1,p1 = ',h1,p1
    print*,'s2',s2
    print*,'h2,p2 = ',h2,p2
+    norm_double_1 += dabs(psi_coef_sorted(i,1)/psi_coef_sorted(1,1))
   endif

   print*,'<Ref| H |D_I> = ',hij
@ -109,6 +112,7 @@ subroutine routine
 print*,''
 print*,'L1 norm of mono alpha = ',norm_mono_a
 print*,'L1 norm of mono beta  = ',norm_mono_b
+ print*,'L1 norm of double exc = ',norm_double_1
 print*, '---'
 print*,'L2 norm of mono alpha = ',norm_mono_a_2
 print*,'L2 norm of mono beta  = ',norm_mono_b_2