Merge pull request #171 from QuantumPackage/cleaning_dft

Cleaning dft

etc/openmp.rc

@@ -0,0 +1 @@
+export OMP_NESTED=True

src/ao_two_e_ints/map_integrals.irp.f

@@ -327,6 +327,8 @@ double precision function get_ao_two_e_integral(i,j,k,l,map) result(result)
   implicit none
   BEGIN_DOC
   ! Gets one AO bi-electronic integral from the AO map
+  !
+  ! i,j,k,l in physicist notation <ij|kl>
   END_DOC
   integer, intent(in)            :: i,j,k,l
   integer(key_kind)              :: idx

src/basis_correction/print_routine.irp.f

@@ -38,7 +38,7 @@ subroutine print_basis_correction
     write(*, '(A29,X,I3,X,A3,X,F16.10)') '  ECMD PBE-UEG       , state ',istate,' = ',ecmd_pbe_ueg_mu_of_r(istate)
    enddo
 
-  else if(mu_of_r_potential.EQ."cas_ful")then
+  else if(mu_of_r_potential.EQ."cas_ful".or.mu_of_r_potential.EQ."cas_truncated".or.mu_of_r_potential.EQ."pure_act")then
    print*, ''
    print*,'Using a CAS-like two-body density to define mu(r)'
    print*,'This assumes that the CAS is a qualitative representation of the wave function ' 

src/determinants/density_matrix.irp.f

@@ -262,17 +262,48 @@ subroutine set_natural_mos
      iorb = list_virt(i)
      do j = 1, n_core_inact_act_orb
       jorb = list_core_inact_act(j)
-      if(one_e_dm_mo(iorb,jorb).ne. 0.d0)then
-        print*,'AHAHAH'
-        print*,iorb,jorb,one_e_dm_mo(iorb,jorb)
-        stop
-      endif
      enddo
     enddo
    call mo_as_svd_vectors_of_mo_matrix_eig(one_e_dm_mo,size(one_e_dm_mo,1),mo_num,mo_num,mo_occ,label)
    soft_touch mo_occ
-
 end
+
+subroutine set_natorb_no_ov_rot
+   implicit none
+   BEGIN_DOC
+   ! Set natural orbitals, obtained by diagonalization of the one-body density matrix
+   ! in the |MO| basis
+   END_DOC
+   character*(64)                 :: label
+   double precision, allocatable  :: tmp(:,:)
+   allocate(tmp(mo_num, mo_num))
+   label = "Natural"
+   tmp = one_e_dm_mo
+    integer :: i,j,iorb,jorb
+    do i = 1, n_virt_orb
+     iorb = list_virt(i)
+     do j = 1, n_core_inact_act_orb
+      jorb = list_core_inact_act(j)
+      tmp(iorb, jorb) = 0.d0
+      tmp(jorb, iorb) = 0.d0
+     enddo
+    enddo
+   call mo_as_svd_vectors_of_mo_matrix_eig(tmp,size(tmp,1),mo_num,mo_num,mo_occ,label)
+   soft_touch mo_occ
+end
+
+subroutine save_natural_mos_no_ov_rot
+   implicit none
+   BEGIN_DOC
+   ! Save natural orbitals, obtained by diagonalization of the one-body density matrix in
+   ! the |MO| basis
+   END_DOC
+   call set_natorb_no_ov_rot
+   call nullify_small_elements(ao_num,mo_num,mo_coef,size(mo_coef,1),1.d-10)
+   call orthonormalize_mos
+   call save_mos
+end
+
 subroutine save_natural_mos
    implicit none
    BEGIN_DOC
@@ -384,6 +415,14 @@ END_PROVIDER
 
 END_PROVIDER
 
+BEGIN_PROVIDER [ double precision, one_e_dm_ao, (ao_num, ao_num)]
+ implicit none
+   BEGIN_DOC
+   !  one_e_dm_ao = one_e_dm_ao_alpha + one_e_dm_ao_beta 
+   END_DOC
+    one_e_dm_ao = one_e_dm_ao_alpha + one_e_dm_ao_beta 
+END_PROVIDER 
+
 
 subroutine get_occupation_from_dets(istate,occupation)
   implicit none

src/mu_of_r/EZFIO.cfg

@@ -6,7 +6,7 @@ size: (becke_numerical_grid.n_points_final_grid,determinants.n_states)
 
 [mu_of_r_potential]
 type: character*(32)
-doc: type of potential for the mu(r) interaction: can be [ hf| cas_ful | cas_truncated]
+doc: type of potential for the mu(r) interaction: can be [ hf| cas_ful | cas_truncated | pure_act]
 interface: ezfio, provider, ocaml
 default: hf
 

src/mu_of_r/basis_def.irp.f

@@ -76,7 +76,11 @@ BEGIN_PROVIDER [integer, n_basis_orb]
  !
  ! It corresponds to all MOs except those defined as "deleted" 
  END_DOC
+ if(mu_of_r_potential == "pure_act")then
+  n_basis_orb = n_act_orb
+ else
   n_basis_orb = n_all_but_del_orb
+ endif
 END_PROVIDER 
 
 BEGIN_PROVIDER [integer, list_basis, (n_basis_orb)]
@@ -89,9 +93,15 @@ BEGIN_PROVIDER [integer, list_basis, (n_basis_orb)]
  ! It corresponds to all MOs except those defined as "deleted" 
  END_DOC
  integer :: i
+ if(mu_of_r_potential == "pure_act")then
+  do i = 1, n_act_orb
+   list_basis(i) = list_act(i)
+  enddo
+ else
   do i = 1, n_all_but_del_orb
     list_basis(i) = list_all_but_del_orb(i)
   enddo
+ endif
 END_PROVIDER 
 
 BEGIN_PROVIDER [double precision, basis_mos_in_r_array, (n_basis_orb,n_points_final_grid)]

src/mu_of_r/mu_of_r_conditions.irp.f

@@ -26,7 +26,7 @@
   do ipoint = 1, n_points_final_grid
    if(mu_of_r_potential.EQ."hf")then
     mu_of_r_prov(ipoint,istate) =  mu_of_r_hf(ipoint)
-   else if(mu_of_r_potential.EQ."cas_ful".or.mu_of_r_potential.EQ."cas_truncated")then
+   else if(mu_of_r_potential.EQ."cas_ful".or.mu_of_r_potential.EQ."cas_truncated".or.mu_of_r_potential.EQ."pure_act")then
     mu_of_r_prov(ipoint,istate) =  mu_of_r_psi_cas(ipoint,istate)
    else 
     print*,'you requested the following mu_of_r_potential'

src/tools/save_natorb_no_ov_rot.irp.f

@@ -0,0 +1,25 @@
+program save_natorb
+  implicit none
+  BEGIN_DOC
+! Save natural |MOs| into the |EZFIO|.
+!
+! This program reads the wave function stored in the |EZFIO| directory,
+! extracts the corresponding natural orbitals and setd them as the new
+! |MOs|.
+!
+! If this is a multi-state calculation, the density matrix that produces
+! the natural orbitals is obtained from an average of the density
+! matrices of each state with the corresponding
+! :option:`determinants state_average_weight`
+  END_DOC
+  read_wf = .True.
+  touch read_wf
+  call save_natural_mos_no_ov_rot
+  call save_ref_determinant
+  call ezfio_set_mo_two_e_ints_io_mo_two_e_integrals('None')
+  call ezfio_set_mo_one_e_ints_io_mo_one_e_integrals('None')
+  call ezfio_set_mo_one_e_ints_io_mo_integrals_kinetic('None')
+  call ezfio_set_mo_one_e_ints_io_mo_integrals_n_e('None')
+  call ezfio_set_mo_one_e_ints_io_mo_integrals_pseudo('None')
+end
+