Merge pull request #34 from QuantumPackage/dev-lct

Dev lct

src/becke_numerical_grid/grid_becke.irp.f

@@ -195,6 +195,20 @@ BEGIN_PROVIDER [double precision, weight_at_r, (n_points_integration_angular,n_p
         enddo
         accu = 1.d0/accu
         weight_at_r(l,k,j) = tmp_array(j) * accu
+        if(isnan(weight_at_r(l,k,j)))then
+         print*,'isnan(weight_at_r(l,k,j))'
+         print*,l,k,j
+         accu = 0.d0
+         do i = 1, nucl_num
+           ! function defined for each atom "i" by equation (13) and (21) with k == 3
+           tmp_array(i) = cell_function_becke(r,i) ! P_n(r)
+           print*,i,tmp_array(i)
+           ! Then you compute the summ the P_n(r) function for each of the "r" points
+           accu += tmp_array(i)
+         enddo
+         write(*,'(100(F16.10,X))')tmp_array(j) , accu
+          stop
+        endif
       enddo
     enddo
   enddo
@@ -221,6 +235,12 @@ BEGIN_PROVIDER [double precision, final_weight_at_r, (n_points_integration_angul
         contrib_integration = derivative_knowles_function(alpha_knowles(int(nucl_charge(j))),m_knowles,x)&
             *knowles_function(alpha_knowles(int(nucl_charge(j))),m_knowles,x)**2
         final_weight_at_r(k,i,j) = weights_angular_points(k)  * weight_at_r(k,i,j) * contrib_integration * dr_radial_integral
+        if(isnan(final_weight_at_r(k,i,j)))then
+         print*,'isnan(final_weight_at_r(k,i,j))' 
+         print*,k,i,j
+         write(*,'(100(F16.10,X))')weights_angular_points(k)  , weight_at_r(k,i,j) , contrib_integration , dr_radial_integral
+         stop 
+        endif
       enddo
     enddo
   enddo

src/becke_numerical_grid/step_function_becke.irp.f

@@ -31,6 +31,10 @@ double precision function cell_function_becke(r,atom_number)
   double precision               :: mu_ij,nu_ij
   double precision               :: distance_i,distance_j,step_function_becke
   integer                        :: j
+  if(int(nucl_charge(atom_number))==0)then
+   cell_function_becke = 0.d0
+   return
+  endif
   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))
@@ -38,6 +42,7 @@ double precision function cell_function_becke(r,atom_number)
   cell_function_becke = 1.d0
   do j = 1, nucl_num
     if(j==atom_number)cycle
+    if(int(nucl_charge(j))==0)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))

src/density_for_dft/EZFIO.cfg

@@ -16,3 +16,10 @@ doc: Type of density
 doc: if [no_core_dm] then all elements of the density matrix involving at least one orbital set as core are set to zero
 interface: ezfio, provider, ocaml
 default: full_density
+
+[normalize_dm]
+type: logical
+doc: Type of density
+doc: if .True., then you normalize the no_core_dm to elec_alpha_num - n_core_orb  and elec_beta_num - n_core_orb
+interface: ezfio, provider, ocaml
+default: True

src/density_for_dft/density_for_dft.irp.f

@@ -29,6 +29,20 @@ BEGIN_PROVIDER [double precision, one_e_dm_mo_alpha_for_dft, (mo_num,mo_num, N_s
     one_e_dm_mo_alpha_for_dft(i,j,:) = 0.d0
    enddo
   enddo
+  if(normalize_dm)then
+   double precision :: elec_alpha_frozen_num, elec_alpha_valence(N_states)
+   elec_alpha_frozen_num = elec_alpha_num - n_core_orb 
+   elec_alpha_valence = 0.d0
+   integer :: istate 
+   do istate = 1, N_states
+    do i = 1, mo_num
+     elec_alpha_valence(istate) += one_e_dm_mo_alpha_for_dft(i,i,istate)
+    enddo
+    elec_alpha_valence(istate) = elec_alpha_frozen_num/elec_alpha_valence(istate)
+    one_e_dm_mo_alpha_for_dft(:,:,istate) = one_e_dm_mo_alpha_for_dft(:,:,istate) * elec_alpha_valence(istate)
+   enddo 
+   
+  endif
  endif
 
 END_PROVIDER
@@ -64,6 +78,19 @@ BEGIN_PROVIDER [double precision, one_e_dm_mo_beta_for_dft, (mo_num,mo_num, N_st
     one_e_dm_mo_beta_for_dft(i,j,:) = 0.d0
    enddo
   enddo
+  if(normalize_dm)then
+   double precision :: elec_beta_valence(N_states),elec_beta_frozen_num
+   elec_beta_frozen_num = elec_beta_num - n_core_orb 
+   elec_beta_valence = 0.d0
+   integer :: istate 
+   do istate = 1, N_states
+    do i = 1, mo_num
+     elec_beta_valence(istate) += one_e_dm_mo_beta_for_dft(i,i,istate)
+    enddo
+    elec_beta_valence(istate) = elec_beta_frozen_num/elec_beta_valence(istate)
+    one_e_dm_mo_beta_for_dft(:,:,istate) = one_e_dm_mo_beta_for_dft(:,:,istate) * elec_beta_valence(istate)
+   enddo 
+  endif
  endif
 END_PROVIDER
 

src/dft_utils_in_r/dm_in_r.irp.f

@@ -227,6 +227,8 @@ END_PROVIDER
 
  BEGIN_PROVIDER [double precision, one_e_dm_alpha_at_r, (n_points_final_grid,N_states) ]
 &BEGIN_PROVIDER [double precision, one_e_dm_beta_at_r, (n_points_final_grid,N_states) ]
+&BEGIN_PROVIDER [double precision, elec_beta_num_grid_becke , (N_states) ]
+&BEGIN_PROVIDER [double precision, elec_alpha_num_grid_becke , (N_states) ]
  implicit none
  BEGIN_DOC
 ! one_e_dm_alpha_at_r(i,istate) = n_alpha(r_i,istate)

src/nuclei/atomic_radii.irp.f

@@ -1,4 +1,4 @@
-BEGIN_PROVIDER [ double precision, slater_bragg_radii, (100)]
+BEGIN_PROVIDER [ double precision, slater_bragg_radii, (0:100)]
  implicit none
  BEGIN_DOC
  ! atomic radii in Angstrom defined in table I of JCP 41, 3199 (1964) Slater
@@ -54,10 +54,10 @@ BEGIN_PROVIDER [ double precision, slater_bragg_radii, (100)]
 
 END_PROVIDER
 
-BEGIN_PROVIDER [double precision, slater_bragg_radii_ua, (100)]
+BEGIN_PROVIDER [double precision, slater_bragg_radii_ua, (0:100)]
  implicit none
  integer :: i
- do i = 1, 100
+ do i = 0, 100
   slater_bragg_radii_ua(i) = slater_bragg_radii(i) * 1.889725989d0
  enddo
 END_PROVIDER