From ca49adfa896bf2cd41b77b6b181246766da80e3a Mon Sep 17 00:00:00 2001 From: TApplencourt Date: Fri, 4 Mar 2016 18:29:59 +0100 Subject: [PATCH] First try of libint --- plugins/Hartree_Fock/debug_libinit.irp.f | 118 +++++++++++++++++++++++ src/AO_Basis/aos.irp.f | 16 +-- src/Integrals_Bielec/libint_module.f90 | 51 ++++++++++ 3 files changed, 174 insertions(+), 11 deletions(-) create mode 100644 plugins/Hartree_Fock/debug_libinit.irp.f create mode 100644 src/Integrals_Bielec/libint_module.f90 diff --git a/plugins/Hartree_Fock/debug_libinit.irp.f b/plugins/Hartree_Fock/debug_libinit.irp.f new file mode 100644 index 00000000..44b6b8f7 --- /dev/null +++ b/plugins/Hartree_Fock/debug_libinit.irp.f @@ -0,0 +1,118 @@ + program debug_libint + use libint_module + + implicit none + double precision :: ao_bielec_integral + + integer, allocatable :: s2bf(:) + double precision, allocatable :: buffer_int(:) + + call init_libint(trim(ezfio_filename)//char(0)) + + integer :: nb_shell_f + nb_shell_f = get_nb_shell() + + allocate(s2bf(2*nb_shell_f)) + call map_shell_to_basis_function_interval(2*nb_shell_f,s2bf) + + integer :: s1, s2,s3,s4 + integer :: bf1,bf2,bf3,bf4 + integer :: bf1_begin,bf2_begin,bf3_begin,bf4_begin + integer :: bf1_end,bf2_end,bf3_end,bf4_end + integer :: n1,n2,n3,n4 + integer :: f1,f2,f3,f4,f1234 + + ! =================== ! + ! Loop over the shell ! + ! =================== ! + + do s1 = 1, nb_shell_f + + print*, s1, "/", nb_shell_f + + bf1_begin = s2bf(2*s1-1) + bf1_end = s2bf(2*s1) + n1 = 1 + bf1_end - bf1_begin + + do s2 = 1, nb_shell_f + + bf2_begin = s2bf(2*s2-1) + bf2_end = s2bf(2*s2) + n2 = 1 + bf2_end - bf2_begin + + do s3 = 1, nb_shell_f + + bf3_begin = s2bf(2*s3-1) + bf3_end = s2bf(2*s3) + n3 = 1 + bf3_end - bf3_begin + + do s4 = 1, nb_shell_f + + bf4_begin = s2bf(2*s4-1) + bf4_end = s2bf(2*s4) + n4 = 1 + bf4_end - bf4_begin + + ! ========================== ! + ! Compute the shell integral ! + ! ========================== ! + integer :: sze + sze = n1*n2*n3*n4 + allocate(buffer_int(sze)) + call compute_ao_bielec_integrals_shell(s1,s2,s3,s4,sze,buffer_int) + + ! ============================ ! + ! Loop over the basis function ! + ! ============================ ! + + do bf1 = bf1_begin, bf1_end + do bf2 = bf2_begin, bf2_end + do bf3 = bf3_begin, bf3_end + do bf4 = bf4_begin, bf4_end + + f1 = bf1 - bf1_begin + f2 = bf2 - bf2_begin + f3 = bf3 - bf3_begin + f4 = bf4 - bf4_begin + + !Get the integral from the buffer + f1234 = f1*n2*n3*n4+f2*n3*n4+f3*n4+f4 + 1; + + !Compute the norm + double precision:: coef1, coef2, coef3, coef4, norm + + coef1 = ao_coef_normalization_factor(bf1) + coef2 = ao_coef_normalization_factor(bf2) + coef3 = ao_coef_normalization_factor(bf3) + coef4 = ao_coef_normalization_factor(bf4) + norm = coef1*coef2*coef3*coef4 + + double precision:: libint, ref + + !Value of itegral bf1,bf2,bf3, bf4 + libint = buffer_int(f1234) * norm + + !Verify with the manu's one +! ref = ao_bielec_integral(bf1,bf2,bf3,bf4) +! +! if ( (ABS(ABS(ref) - ABS(libint)).GE.1e-6) ) THEN +! print*, bf1,bf2,bf3,bf4 +! print*, ref +! print*, libint +! call exit(1) +! end if + + enddo + enddo + enddo + enddo + + !Deallocate the buffer_intergral for the shell + deallocate(buffer_int) + + enddo + enddo + enddo + enddo + + call finalize_libint() + end debug_libint diff --git a/src/AO_Basis/aos.irp.f b/src/AO_Basis/aos.irp.f index 04c90ca7..793a7e2c 100644 --- a/src/AO_Basis/aos.irp.f +++ b/src/AO_Basis/aos.irp.f @@ -34,9 +34,9 @@ END_PROVIDER C_A(3) = 0.d0 ao_coef_normalized = 0.d0 do i=1,ao_num - powA(1) = ao_power(i,1) - powA(2) = ao_power(i,2) - powA(3) = ao_power(i,3) + powA(1) = ao_l(i) + powA(2) = 0 + powA(3) = 0 do j=1,ao_prim_num(i) call overlap_gaussian_xyz(C_A,C_A,ao_expo(i,j),ao_expo(i,j),powA,powA,overlap_x,overlap_y,overlap_z,norm,nz) ao_coef_normalized(i,j) = ao_coef(i,j)/sqrt(norm) @@ -55,12 +55,6 @@ END_PROVIDER enddo enddo -! do i=1,ao_num -! do j=1,ao_prim_num(i) -! ao_coef_normalized(i,j) = ao_coef(i,j) -! enddo -! enddo - END_PROVIDER BEGIN_PROVIDER [ double precision, ao_coef_normalized_ordered, (ao_num_align,ao_prim_num_max) ] @@ -79,8 +73,8 @@ END_PROVIDER d(j,1) = ao_expo(i,j) d(j,2) = ao_coef_normalized(i,j) enddo -! call dsort(d(1,1),iorder,ao_prim_num(i)) -! call dset_order(d(1,2),iorder,ao_prim_num(i)) + call dsort(d(1,1),iorder,ao_prim_num(i)) + call dset_order(d(1,2),iorder,ao_prim_num(i)) do j=1,ao_prim_num(i) ao_expo_ordered(i,j) = d(j,1) ao_coef_normalized_ordered(i,j) = d(j,2) diff --git a/src/Integrals_Bielec/libint_module.f90 b/src/Integrals_Bielec/libint_module.f90 new file mode 100644 index 00000000..2cedf283 --- /dev/null +++ b/src/Integrals_Bielec/libint_module.f90 @@ -0,0 +1,51 @@ +module libint_module + use iso_c_binding + + implicit none + interface + subroutine init_libint(str) bind(c,name='init_libint') + import :: c_char + character(len=1,kind=C_char), dimension(*), intent(in) :: str + end subroutine + + integer(c_int) function get_nb_shell() bind(c,name='nb_shell') + import :: c_int + end function + + subroutine finalize_libint() bind(c,name='finalize_libint') + end subroutine + + subroutine map_shell_to_basis_function_interval(sze, out_val) bind(c,name='map_shell_to_basis_function_interval') + import :: c_ptr + import :: c_int + + integer(c_int), INTENT(IN), value :: sze + integer(c_int), INTENT(OUT) :: out_val(sze) + end subroutine + + real(c_double) function ao_bielec_integral_libint(i,j,k,l) bind(c,name='ao_bielec_integral') + import :: c_int + import :: c_double + + integer(c_int), value :: i + integer(c_int), value :: j + integer(c_int), value :: k + integer(c_int), value :: l + end function + + subroutine compute_ao_bielec_integrals_shell(i,j,k,l,sze,values) bind(c,name='compute_ao_bielec_integrals_shell') + import :: c_ptr + import :: c_int + import :: c_double + + integer(c_int), value :: i + integer(c_int), value :: j + integer(c_int), value :: k + integer(c_int), value :: l + integer(c_int), INTENT(IN), value :: sze + real(c_double), INTENT(OUT) :: values(sze) + end subroutine + + + end interface +end module libint_module