Split pseudo in local/non-local

src/FCIdump/README.rst

@@ -10,9 +10,6 @@ Documentation
 .. Do not edit this section. It was auto-generated from the
 .. NEEDED_MODULES file.
 
-`fcidump <http://github.com/LCPQ/quantum_package/tree/master/src/FCIdump/fcidump.irp.f#L1>`_
-  Undocumented
-
 
 
 Needed Modules

src/Pseudo_integrals/README.rst

@@ -9,6 +9,12 @@ Documentation
 .. NEEDED_MODULES file.
 
 `ao_nucl_elec_integral_pseudo <http://github.com/LCPQ/quantum_package/tree/master/src/Pseudo_integrals/pot_ao_ints_pseudo.irp.f#L1>`_
+  Undocumented
+
+`ao_nucl_elec_integral_pseudo_local <http://github.com/LCPQ/quantum_package/tree/master/src/Pseudo_integrals/pot_ao_ints_pseudo.irp.f#L10>`_
+  Local component of the pseudopotential
+
+`ao_nucl_elec_integral_pseudo_non_local <http://github.com/LCPQ/quantum_package/tree/master/src/Pseudo_integrals/pot_ao_ints_pseudo.irp.f#L124>`_
   interaction nuclear electron
 
 

src/Pseudo_integrals/int.f90

@@ -1003,11 +1003,11 @@ end
 double precision function crochet(n,g)
 implicit none
 integer n
-double precision g,pi,dble_fact,expo
+double precision g,dble_fact,expo
-pi=dacos(-1.d0)
+double precision, parameter :: sq_pi_ov_2=dsqrt(dacos(-1.d0)*0.5d0)
 expo=0.5d0*dfloat(n+1)
 crochet=dble_fact(n-1)/(2.d0*g)**expo
-if(mod(n,2).eq.0)crochet=crochet*dsqrt(pi/2.d0)
+if(mod(n,2).eq.0)crochet=crochet*sq_pi_ov_2
 end
 
 !!
@@ -2043,10 +2043,10 @@ double precision function int_prod_bessel_loc(l,gam,n,a)
   int=0
   
   ! Int f_0
-  coef_nk=1.d0/dble_fact(2*n+1)
+  coef_nk=1.d0/dble_fact( n+n+1 )
   expo=0.5d0*dfloat(n+l+1)
   crochet=dble_fact(n+l-1)/(2.d0*gam)**expo
-  if(mod(n+l,2).eq.0)crochet=crochet*dsqrt(pi/2.d0)
+  if(mod(n+l,2).eq.0)crochet=crochet*dsqrt(0.5d0*pi)
 
   f_0 = coef_nk*a**n*crochet
 

src/Pseudo_integrals/pot_ao_ints_pseudo.irp.f

@@ -1,8 +1,17 @@
- BEGIN_PROVIDER [ double precision, ao_nucl_elec_integral_pseudo, (ao_num_align,ao_num)]
+BEGIN_PROVIDER [ double precision, ao_nucl_elec_integral_pseudo, (ao_num_align,ao_num) ]
- BEGIN_DOC
-! interaction nuclear electron
- END_DOC
  implicit none
+! Nucleus-pseudopotential interaction
+ END_DOC
+ write(output_monoints,*) 'Providing the pseudopotential integrals '
+ ao_nucl_elec_integral_pseudo = ao_nucl_elec_integral_pseudo_local + ao_nucl_elec_integral_pseudo_non_local
+END_PROVIDER
+
+
+BEGIN_PROVIDER [ double precision, ao_nucl_elec_integral_pseudo_local, (ao_num_align,ao_num)]
+ implicit none
+ BEGIN_DOC
+! Local component of the pseudopotential
+ END_DOC
  double precision  :: alpha, beta, gama, delta
  integer           :: num_A,num_B
  double precision  :: A_center(3),B_center(3),C_center(3)
@@ -13,13 +22,8 @@
  double precision  :: cpu_1, cpu_2, wall_1, wall_2, wall_0
  integer           :: thread_num
 
-  ao_nucl_elec_integral_pseudo = 0.d0
+  PROVIDE output_monoints
-
+  ao_nucl_elec_integral_pseudo_local = 0.d0
-  !                 
-  ! |   _   _  _. | 
-  ! |_ (_) (_ (_| | 
-  !                 
-  !! Parameters of the local part of pseudo:
 
   integer klocmax
   integer, allocatable ::  n_k(:,:)
@@ -40,41 +44,6 @@
   allocate(n_k_dump(1:klocmax), v_k_dump(1:klocmax), dz_k_dump(1:klocmax)) 
 
 
-  !                               
-  ! |\ |  _  ._    |  _   _  _. | 
-  ! | \| (_) | |   | (_) (_ (_| | 
-  !                              
-  !! Parameters of non local part of pseudo:
-
-  integer :: kmax,lmax
-  integer, allocatable ::  n_kl(:,:,:)
-  double precision, allocatable ::  v_kl(:,:,:), dz_kl(:,:,:) 
-
-  call ezfio_get_pseudo_integrals_lmaxpo(lmax)
-  call ezfio_get_pseudo_integrals_kmax(kmax)
-  !lmax plus one -> lmax
-  lmax = lmax - 1 
- 
-  allocate(n_kl(nucl_num,kmax,0:lmax), v_kl(nucl_num,kmax,0:lmax), dz_kl(nucl_num,kmax,0:lmax)) 
-
-  call ezfio_get_pseudo_integrals_n_kl(n_kl)
-  call ezfio_get_pseudo_integrals_v_kl(v_kl)
-  call ezfio_get_pseudo_integrals_dz_kl(dz_kl)
-
-
-  !! Dump array 
-  integer, allocatable ::  n_kl_dump(:,:)
-  double precision, allocatable ::  v_kl_dump(:,:), dz_kl_dump(:,:) 
- 
-  allocate(n_kl_dump(kmax,0:lmax), v_kl_dump(kmax,0:lmax), dz_kl_dump(kmax,0:lmax)) 
-
-  !  _                
-  ! /   _. |  _     | 
-  ! \_ (_| | (_ |_| | 
-  !                   
-
-  write(output_monoints,*) 'Providing the nuclear electron pseudo integrals '
-
   call wall_time(wall_1)
   call cpu_time(cpu_1)
 
@@ -82,11 +51,11 @@
   !$OMP DEFAULT (NONE) &
   !$OMP PRIVATE (i,j,k,l,m,alpha,beta,A_center,B_center,C_center,power_A,power_B, &
   !$OMP          num_A,num_B,Z,c,n_pt_in, &
-  !$OMP          v_k_dump,n_k_dump, dz_k_dump, n_kl_dump, v_kl_dump, dz_kl_dump, &
+  !$OMP          v_k_dump,n_k_dump, dz_k_dump, &
   !$OMP          wall_0,wall_2,thread_num, output_monoints) & 
   !$OMP SHARED (ao_num,ao_prim_num,ao_expo_ordered_transp,ao_power,ao_nucl,nucl_coord,ao_coef_normalized_ordered_transp, &
-  !$OMP         ao_nucl_elec_integral_pseudo,nucl_num,nucl_charge, &
+  !$OMP         ao_nucl_elec_integral_pseudo_local,nucl_num,nucl_charge, &
-  !$OMP         klocmax,lmax,kmax,v_k,n_k, dz_k, n_kl, v_kl, dz_kl, &
+  !$OMP         klocmax,v_k,n_k, dz_k, &
   !$OMP         wall_1)
   
   !$OMP DO SCHEDULE (guided)
@@ -124,15 +93,8 @@
         c = c + Vloc(klocmax, v_k_dump,n_k_dump, dz_k_dump, &
                      A_center,power_A,alpha,B_center,power_B,beta,C_center)
   
-
-        n_kl_dump = n_kl(k,1:kmax,0:lmax)
-        v_kl_dump = v_kl(k,1:kmax,0:lmax)
-        dz_kl_dump = dz_kl(k,1:kmax,0:lmax)
-
-        c = c + Vpseudo(lmax,kmax,v_kl_dump,n_kl_dump,dz_kl_dump,A_center,power_A,alpha,B_center,power_B,beta,C_center)
-  
       enddo
-      ao_nucl_elec_integral_pseudo(i,j) = ao_nucl_elec_integral_pseudo(i,j) + &
+      ao_nucl_elec_integral_pseudo_local(i,j) = ao_nucl_elec_integral_pseudo_local(i,j) + &
                                    ao_coef_normalized_ordered_transp(l,j)*ao_coef_normalized_ordered_transp(m,i)*c
      enddo
      enddo
@@ -152,16 +114,134 @@
  !$OMP END PARALLEL
 
 
-!  _                                 
-! | \  _   _. | |  _   _  _. _|_  _  
-! |_/ (/_ (_| | | (_) (_ (_|  |_ (/_ 
-!                                    
- 
   deallocate(n_k,v_k, dz_k)
   deallocate(n_k_dump,v_k_dump, dz_k_dump)
 
+
+
+END_PROVIDER
+
+
+ BEGIN_PROVIDER [ double precision, ao_nucl_elec_integral_pseudo_non_local, (ao_num_align,ao_num)]
+ BEGIN_DOC
+! interaction nuclear electron
+ END_DOC
+ implicit none
+ double precision  :: alpha, beta, gama, delta
+ integer           :: num_A,num_B
+ double precision  :: A_center(3),B_center(3),C_center(3)
+ integer           :: power_A(3),power_B(3)
+ integer           :: i,j,k,l,n_pt_in,m
+ double precision  ::  Vloc, Vpseudo
+
+ double precision  :: cpu_1, cpu_2, wall_1, wall_2, wall_0
+ integer           :: thread_num
+
+  PROVIDE output_monoints
+  ao_nucl_elec_integral_pseudo_non_local = 0.d0
+
+  integer :: kmax,lmax
+  integer, allocatable ::  n_kl(:,:,:)
+  double precision, allocatable ::  v_kl(:,:,:), dz_kl(:,:,:) 
+
+  call ezfio_get_pseudo_integrals_lmaxpo(lmax)
+  call ezfio_get_pseudo_integrals_kmax(kmax)
+  !lmax plus one -> lmax
+  lmax = lmax - 1 
+ 
+  allocate(n_kl(nucl_num,kmax,0:lmax), v_kl(nucl_num,kmax,0:lmax), dz_kl(nucl_num,kmax,0:lmax)) 
+
+  call ezfio_get_pseudo_integrals_n_kl(n_kl)
+  call ezfio_get_pseudo_integrals_v_kl(v_kl)
+  call ezfio_get_pseudo_integrals_dz_kl(dz_kl)
+
+
+  !! Dump array 
+  integer, allocatable ::  n_kl_dump(:,:)
+  double precision, allocatable ::  v_kl_dump(:,:), dz_kl_dump(:,:) 
+ 
+  allocate(n_kl_dump(kmax,0:lmax), v_kl_dump(kmax,0:lmax), dz_kl_dump(kmax,0:lmax)) 
+
+  !  _                
+  ! /   _. |  _     | 
+  ! \_ (_| | (_ |_| | 
+  !                   
+
+  call wall_time(wall_1)
+  call cpu_time(cpu_1)
+
+  !$OMP PARALLEL &
+  !$OMP DEFAULT (NONE) &
+  !$OMP PRIVATE (i,j,k,l,m,alpha,beta,A_center,B_center,C_center,power_A,power_B, &
+  !$OMP          num_A,num_B,Z,c,n_pt_in, &
+  !$OMP          n_kl_dump, v_kl_dump, dz_kl_dump, &
+  !$OMP          wall_0,wall_2,thread_num, output_monoints) & 
+  !$OMP SHARED (ao_num,ao_prim_num,ao_expo_ordered_transp,ao_power,ao_nucl,nucl_coord,ao_coef_normalized_ordered_transp, &
+  !$OMP         ao_nucl_elec_integral_pseudo_non_local,nucl_num,nucl_charge, &
+  !$OMP         lmax,kmax,v_kl, dz_kl, n_kl, &
+  !$OMP         wall_1)
+  
+  !$OMP DO SCHEDULE (guided)
+
+  do j = 1, ao_num
+
+   num_A = ao_nucl(j)
+   power_A(1:3)= ao_power(j,1:3)
+   A_center(1:3) = nucl_coord(num_A,1:3)
+
+  do i = 1, ao_num
+
+   num_B = ao_nucl(i)
+   power_B(1:3)= ao_power(i,1:3)
+   B_center(1:3) = nucl_coord(num_B,1:3)
+
+    do l=1,ao_prim_num(j)
+     alpha = ao_expo_ordered_transp(l,j)
+  
+    do m=1,ao_prim_num(i)
+      beta = ao_expo_ordered_transp(m,i)
+      double precision :: c
+      c = 0.d0
+       
+      do  k = 1, nucl_num
+        double precision :: Z
+        Z = nucl_charge(k)
+  
+        C_center(1:3) = nucl_coord(k,1:3)
+  
+        n_kl_dump = n_kl(k,1:kmax,0:lmax)
+        v_kl_dump = v_kl(k,1:kmax,0:lmax)
+        dz_kl_dump = dz_kl(k,1:kmax,0:lmax)
+
+        c = c + Vpseudo(lmax,kmax,v_kl_dump,n_kl_dump,dz_kl_dump,A_center,power_A,alpha,B_center,power_B,beta,C_center)
+  
+      enddo
+      ao_nucl_elec_integral_pseudo_non_local(i,j) = ao_nucl_elec_integral_pseudo_non_local(i,j) + &
+                                   ao_coef_normalized_ordered_transp(l,j)*ao_coef_normalized_ordered_transp(m,i)*c
+     enddo
+     enddo
+  enddo
+
+    call wall_time(wall_2)
+    if (thread_num == 0) then
+      if (wall_2 - wall_0 > 1.d0) then
+        wall_0 = wall_2
+        write(output_monoints,*) 100.*float(j)/float(ao_num), '% in ',  &
+                                 wall_2-wall_1, 's'
+      endif
+    endif
+  enddo
+
+ !$OMP END DO
+ !$OMP END PARALLEL
+
+ 
   deallocate(n_kl,v_kl, dz_kl)
   deallocate(n_kl_dump,v_kl_dump, dz_kl_dump)
 
 
 END_PROVIDER
+
+
+
+

src/Utils/README.rst

@@ -242,7 +242,7 @@ Documentation
 `align_double <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/util.irp.f#L48>`_
   Compute 1st dimension such that it is aligned for vectorization.
 
-`approx_dble <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/util.irp.f#L341>`_
+`approx_dble <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/util.irp.f#L380>`_
   Undocumented
 
 `binom <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/util.irp.f#L31>`_
@@ -257,10 +257,16 @@ Documentation
 `binom_transp <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/util.irp.f#L32>`_
   Binomial coefficients
 
-`dble_fact <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/util.irp.f#L137>`_
+`dble_fact <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/util.irp.f#L138>`_
+  Undocumented
+
+`dble_fact_odd <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/util.irp.f#L176>`_
   n!!
 
-`dble_logfact <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/util.irp.f#L171>`_
+`dble_fact_peer <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/util.irp.f#L155>`_
+  n!!
+
+`dble_logfact <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/util.irp.f#L210>`_
   n!!
 
 `fact <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/util.irp.f#L63>`_
@@ -269,29 +275,29 @@ Documentation
 `fact_inv <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/util.irp.f#L125>`_
   1/n!
 
-`inv_int <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/util.irp.f#L218>`_
+`inv_int <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/util.irp.f#L257>`_
   1/i
 
 `logfact <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/util.irp.f#L93>`_
   n!
 
-`normalize <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/util.irp.f#L317>`_
+`normalize <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/util.irp.f#L356>`_
   Normalizes vector u
   u is expected to be aligned in memory.
 
-`nproc <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/util.irp.f#L244>`_
+`nproc <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/util.irp.f#L283>`_
   Number of current OpenMP threads
 
-`u_dot_u <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/util.irp.f#L286>`_
+`u_dot_u <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/util.irp.f#L325>`_
   Compute <u|u>
 
-`u_dot_v <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/util.irp.f#L260>`_
+`u_dot_v <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/util.irp.f#L299>`_
   Compute <u|v>
 
-`wall_time <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/util.irp.f#L229>`_
+`wall_time <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/util.irp.f#L268>`_
   The equivalent of cpu_time, but for the wall time.
 
-`write_git_log <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/util.irp.f#L204>`_
+`write_git_log <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/util.irp.f#L243>`_
   Write the last git commit in file iunit.
 
 

src/Utils/util.irp.f

@@ -138,21 +138,21 @@ END_PROVIDER
 double precision function dble_fact(n)
   implicit none
   integer :: n
-  double precision :: dble_fact_peer, dble_fact_odd
+  double precision :: dble_fact_even, dble_fact_odd
 
   dble_fact = 1.d0
 
   if(n.lt.0) return
 
   if(iand(n,1).eq.0)then
-    dble_fact = dble_fact_peer(n)
+    dble_fact = dble_fact_even(n)
   else
     dble_fact= dble_fact_odd(n)
   endif
 
 end function
 
-double precision function dble_fact_peer(n) result(fact2)
+double precision function dble_fact_even(n) result(fact2)
   implicit none
   BEGIN_DOC
   ! n!!