config/ifort_avx.cfg

@@ -32,7 +32,7 @@ OPENMP  : 1          ; Append OpenMP flags
 #
 [OPT]
 FC       : -traceback
-FCFLAGS  : -march=corei7-avx -O2 -ip -ftz -g 
+FCFLAGS  : -xAVX -O2 -ip -ftz -g 
 
 # Profiling flags
 #################

config/ifort_avx_mpi.cfg

@@ -31,14 +31,14 @@ OPENMP  : 1          ; Append OpenMP flags
 # -ftz                       : Flushes denormal results to zero
 #
 [OPT]
-FCFLAGS  : -march=corei7-avx -O2 -ip -ftz -g -traceback  
+FCFLAGS  : -xAVX -O2 -ip -ftz -g -traceback  
 
 # Profiling flags
 #################
 #
 [PROFILE]
 FC       : -p -g
-FCFLAGS  : -march=corei7 -O2 -ip -ftz 
+FCFLAGS  : -xSSE4.2 -O2 -ip -ftz 
 
 
 # Debugging flags

src/ao_two_e_erf_ints/two_e_integrals_erf.irp.f

@@ -16,8 +16,6 @@ double precision function ao_two_e_integral_erf(i,j,k,l)
   integer                        :: iorder_p(3), iorder_q(3)
   double precision               :: ao_two_e_integral_schwartz_accel_erf
 
-  provide mu_erf
-
    if (ao_prim_num(i) * ao_prim_num(j) * ao_prim_num(k) * ao_prim_num(l) > 1024 ) then
      ao_two_e_integral_erf = ao_two_e_integral_schwartz_accel_erf(i,j,k,l)
      return

src/ao_two_e_ints/map_integrals.irp.f

@@ -279,100 +279,6 @@ subroutine get_ao_two_e_integrals_non_zero(j,k,l,sze,out_val,out_val_index,non_z
 end
 
 
-subroutine get_ao_two_e_integrals_non_zero_jl(j,l,thresh,sze_max,sze,out_val,out_val_index,non_zero_int)
-  use map_module
-  implicit none
-  BEGIN_DOC
-  ! Gets multiple AO bi-electronic integral from the AO map .
-  ! All non-zero i are retrieved for j,k,l fixed.
-  END_DOC
-  double precision, intent(in)   :: thresh
-  integer, intent(in)            :: j,l, sze,sze_max
-  real(integral_kind), intent(out) :: out_val(sze_max)
-  integer, intent(out)           :: out_val_index(2,sze_max),non_zero_int
-
-  integer                        :: i,k
-  integer(key_kind)              :: hash
-  double precision               :: tmp
-
-  PROVIDE ao_two_e_integrals_in_map
-  non_zero_int = 0
-  if (ao_overlap_abs(j,l) < thresh) then
-    out_val = 0.d0
-    return
-  endif
-
-  non_zero_int = 0
-  do k = 1, sze
-   do i = 1, sze
-     integer, external :: ao_l4
-     double precision, external :: ao_two_e_integral
-     !DIR$ FORCEINLINE
-     if (ao_two_e_integral_schwartz(i,k)*ao_two_e_integral_schwartz(j,l) < thresh) then
-       cycle
-     endif
-     call two_e_integrals_index(i,j,k,l,hash)
-     call map_get(ao_integrals_map, hash,tmp)
-     if (dabs(tmp) < thresh ) cycle
-     non_zero_int = non_zero_int+1
-     out_val_index(1,non_zero_int) = i
-     out_val_index(2,non_zero_int) = k
-     out_val(non_zero_int) = tmp
-   enddo
-  enddo
-
-end
-
-
-subroutine get_ao_two_e_integrals_non_zero_jl_from_list(j,l,thresh,list,n_list,sze_max,out_val,out_val_index,non_zero_int)
-  use map_module
-  implicit none
-  BEGIN_DOC
-  ! Gets multiple AO two-electron integrals from the AO map .
-  ! All non-zero i are retrieved for j,k,l fixed.
-  END_DOC
-  double precision, intent(in)   :: thresh
-  integer, intent(in)            :: sze_max
-  integer, intent(in)            :: j,l, n_list,list(2,sze_max)
-  real(integral_kind), intent(out) :: out_val(sze_max)
-  integer, intent(out)           :: out_val_index(2,sze_max),non_zero_int
-
-  integer                        :: i,k
-  integer(key_kind)              :: hash
-  double precision               :: tmp
-
-  PROVIDE ao_two_e_integrals_in_map
-  non_zero_int = 0
-  if (ao_overlap_abs(j,l) < thresh) then
-    out_val = 0.d0
-    return
-  endif
-
-  non_zero_int = 0
- integer :: kk
-  do kk = 1, n_list
-   k = list(1,kk)
-   i = list(2,kk)
-   integer, external :: ao_l4
-   double precision, external :: ao_two_e_integral
-   !DIR$ FORCEINLINE
-   if (ao_two_e_integral_schwartz(i,k)*ao_two_e_integral_schwartz(j,l) < thresh) then
-     cycle
-   endif
-   call two_e_integrals_index(i,j,k,l,hash)
-   call map_get(ao_integrals_map, hash,tmp)
-   if (dabs(tmp) < thresh ) cycle
-   non_zero_int = non_zero_int+1
-   out_val_index(1,non_zero_int) = i
-   out_val_index(2,non_zero_int) = k
-   out_val(non_zero_int) = tmp
-  enddo
-
-end
-
-
-
-
 function get_ao_map_size()
   implicit none
   integer (map_size_kind) :: get_ao_map_size

src/becke_numerical_grid/EZFIO.cfg

@@ -8,9 +8,3 @@ default: 2
 type: integer
 doc: Total number of grid points 
 interface: ezfio
-
-[thresh_grid]
-type: double precision
-doc: threshold on the weight of a given grid point
-interface: ezfio,provider,ocaml
-default: 1.e-20

src/becke_numerical_grid/atomic_number.irp.f

@@ -1,9 +0,0 @@
-BEGIN_PROVIDER [ integer, grid_atomic_number, (nucl_num) ]
- implicit none
- BEGIN_DOC
- ! Atomic number used to adjust the grid
- END_DOC
- grid_atomic_number(:) = max(1,int(nucl_charge(:)))
-
-END_PROVIDER
-

src/becke_numerical_grid/grid_becke.irp.f

@@ -146,7 +146,7 @@ BEGIN_PROVIDER [double precision, grid_points_per_atom, (3,n_points_integration_
       x = grid_points_radial(j)
 
       ! value of the radial coordinate for the integration
-      r = knowles_function(alpha_knowles(grid_atomic_number(i)),m_knowles,x)
+      r = knowles_function(alpha_knowles(int(nucl_charge(i))),m_knowles,x)
 
       ! explicit values of the grid points centered around each atom
       do k = 1, n_points_integration_angular
@@ -232,8 +232,8 @@ BEGIN_PROVIDER [double precision, final_weight_at_r, (n_points_integration_angul
     do i = 1, n_points_radial_grid -1 !for each radial grid attached to the "jth" atom
       x = grid_points_radial(i) ! x value for the mapping of the [0, +\infty] to [0,1]
       do k = 1, n_points_integration_angular  ! for each angular point attached to the "jth" atom
-        contrib_integration = derivative_knowles_function(alpha_knowles(grid_atomic_number(j)),m_knowles,x)&
+        contrib_integration = derivative_knowles_function(alpha_knowles(int(nucl_charge(j))),m_knowles,x)&
-            *knowles_function(alpha_knowles(grid_atomic_number(j)),m_knowles,x)**2
+            *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))' 

src/becke_numerical_grid/grid_becke_per_atom.irp.f

@@ -1,53 +0,0 @@
-
-
- BEGIN_PROVIDER [integer, n_pts_per_atom, (nucl_num)]
-&BEGIN_PROVIDER [integer, n_pts_max_per_atom]
-  BEGIN_DOC
-  ! Number of points which are non zero
-  END_DOC
-  integer                        :: i,j,k,l
-  n_pts_per_atom = 0
-  do j = 1, nucl_num
-    do i = 1, n_points_radial_grid -1
-      do k = 1, n_points_integration_angular
-        if(dabs(final_weight_at_r(k,i,j)) < thresh_grid)then
-          cycle
-        endif
-        n_pts_per_atom(j) += 1
-      enddo
-    enddo
-  enddo
-  n_pts_max_per_atom = maxval(n_pts_per_atom)
-END_PROVIDER
-
- BEGIN_PROVIDER [double precision, final_grid_points_per_atom, (3,n_pts_max_per_atom,nucl_num)]
-&BEGIN_PROVIDER [double precision, final_weight_at_r_vector_per_atom, (n_pts_max_per_atom,nucl_num) ]
-&BEGIN_PROVIDER [integer, index_final_points_per_atom, (3,n_pts_max_per_atom,nucl_num) ]
-&BEGIN_PROVIDER [integer, index_final_points_per_atom_reverse, (n_points_integration_angular,n_points_radial_grid,nucl_num) ]
- implicit none
- integer                        :: i,j,k,l,i_count(nucl_num)
- double precision               :: r(3)
- i_count = 0
-  do j = 1, nucl_num
-    do i = 1, n_points_radial_grid -1
-      do k = 1, n_points_integration_angular
-        if(dabs(final_weight_at_r(k,i,j)) < thresh_grid)then
-          cycle
-        endif
-        i_count(j) += 1
-        final_grid_points_per_atom(1,i_count(j),j) = grid_points_per_atom(1,k,i,j)
-        final_grid_points_per_atom(2,i_count(j),j) = grid_points_per_atom(2,k,i,j)
-        final_grid_points_per_atom(3,i_count(j),j) = grid_points_per_atom(3,k,i,j)
-        final_weight_at_r_vector_per_atom(i_count(j),j) = final_weight_at_r(k,i,j)
-        index_final_points_per_atom(1,i_count(j),j) = k
-        index_final_points_per_atom(2,i_count(j),j) = i
-        index_final_points_per_atom(3,i_count(j),j) = j
-        index_final_points_per_atom_reverse(k,i,j) = i_count(j)
-      enddo
-    enddo
-  enddo
-
-
-
-
-END_PROVIDER 

src/becke_numerical_grid/grid_becke_vector.irp.f

@@ -1,6 +1,5 @@
 
 BEGIN_PROVIDER [integer, n_points_final_grid]
-  implicit none
   BEGIN_DOC
   ! Number of points which are non zero
   END_DOC
@@ -9,9 +8,9 @@ BEGIN_PROVIDER [integer, n_points_final_grid]
   do j = 1, nucl_num
     do i = 1, n_points_radial_grid -1
       do k = 1, n_points_integration_angular
-        if(dabs(final_weight_at_r(k,i,j)) < thresh_grid)then
+!       if(dabs(final_weight_at_r(k,i,j)) < 1.d-30)then
-          cycle
+!         cycle
-        endif
+!       endif
         n_points_final_grid += 1
       enddo
     enddo
@@ -40,9 +39,9 @@ END_PROVIDER
   do j = 1, nucl_num
     do i = 1, n_points_radial_grid -1
       do k = 1, n_points_integration_angular
-        if(dabs(final_weight_at_r(k,i,j)) < thresh_grid)then
+       !if(dabs(final_weight_at_r(k,i,j)) < 1.d-30)then
-          cycle
+       !  cycle
-        endif
+       !endif
         i_count += 1
         final_grid_points(1,i_count) = grid_points_per_atom(1,k,i,j)
         final_grid_points(2,i_count) = grid_points_per_atom(2,k,i,j)

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/dft_utils_in_r/ao_in_r.irp.f

@@ -121,26 +121,3 @@
  enddo
  END_PROVIDER
 
-
- BEGIN_PROVIDER[double precision, aos_in_r_array_per_atom, (ao_num,n_pts_max_per_atom,nucl_num)]
-&BEGIN_PROVIDER[double precision, aos_in_r_array_per_atom_transp, (n_pts_max_per_atom,ao_num,nucl_num)]
- implicit none
- BEGIN_DOC
- ! aos_in_r_array_per_atom(i,j,k)        = value of the ith ao on the jth grid point attached on the kth atom 
- END_DOC
- integer :: i,j,k
- double precision :: aos_array(ao_num), r(3)
- do k = 1, nucl_num
-  do i = 1, n_pts_per_atom(k)
-   r(1) = final_grid_points_per_atom(1,i,k)
-   r(2) = final_grid_points_per_atom(2,i,k)
-   r(3) = final_grid_points_per_atom(3,i,k)
-   call give_all_aos_at_r(r,aos_array)
-   do j = 1, ao_num
-    aos_in_r_array_per_atom(j,i,k) = aos_array(j)
-    aos_in_r_array_per_atom_transp(i,j,k) = aos_array(j)
-   enddo
-  enddo
- enddo
- END_PROVIDER
-

src/nuclei/atomic_radii.irp.f

@@ -66,7 +66,7 @@ BEGIN_PROVIDER [double precision, slater_bragg_radii_per_atom, (nucl_num)]
  implicit none
  integer :: i
  do i = 1, nucl_num
-  slater_bragg_radii_per_atom(i) = slater_bragg_radii(max(1,int(nucl_charge(i))))
+  slater_bragg_radii_per_atom(i) = slater_bragg_radii(int(nucl_charge(i)))
  enddo
 END_PROVIDER
 
@@ -74,7 +74,7 @@ BEGIN_PROVIDER [double precision, slater_bragg_radii_per_atom_ua, (nucl_num)]
  implicit none
  integer :: i
  do i = 1, nucl_num
-  slater_bragg_radii_per_atom_ua(i) = slater_bragg_radii_ua(max(1,int(nucl_charge(i))))
+  slater_bragg_radii_per_atom_ua(i) = slater_bragg_radii_ua(int(nucl_charge(i)))
  enddo
 END_PROVIDER