From cd418ee484509df68e12b834bcd00ccd67a0bc33 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 17 Apr 2014 23:50:51 +0200 Subject: [PATCH] Added Bi-electronic integrals module --- src/BiInts/ASSUMPTIONS.rst | 0 src/BiInts/Makefile | 8 + src/BiInts/NEEDED_MODULES | 1 + src/BiInts/README.rst | 27 + src/BiInts/ao_bi_integrals.irp.f | 1016 ++++++++++++++++++++ src/BiInts/bi_integrals.ezfio_config | 8 + src/BiInts/gauss_legendre.irp.f | 57 ++ src/BiInts/map_integrals.irp.f | 334 +++++++ src/BiInts/mo_bi_integrals.irp.f | 356 +++++++ src/BiInts/options.irp.f | 109 +++ src/BiInts/tests/Makefile | 33 + src/BiInts/tests/test_ao.ref | 424 ++++++++ src/BiInts/tests/test_ao_map.irp.f | 35 + src/BiInts/tests/test_ao_map.ref | 635 ++++++++++++ src/BiInts/tests/test_mo.irp.f | 15 + src/BiInts/tests/test_mo.ref | 846 ++++++++++++++++ src/Bitmask/ASSUMPTIONS.rst | 4 +- src/Bitmask/README.rst | 4 +- src/MOGuess/ASSUMPTIONS.rst | 0 src/MOGuess/H_CORE_guess.irp.f | 11 + src/MOGuess/Makefile | 8 + src/MOGuess/README.rst | 4 + src/MOGuess/localize.irp.f | 114 +++ src/MOGuess/mo_ortho_lowdin.irp.f | 46 + src/MOGuess/pot_mo_ortho_lowdin_ints.irp.f | 25 + src/MOGuess/tests/Makefile | 33 + src/Utils/Makefile | 4 +- src/Utils/map_module.f90 | 640 ++++++++++++ 28 files changed, 4793 insertions(+), 4 deletions(-) create mode 100644 src/BiInts/ASSUMPTIONS.rst create mode 100644 src/BiInts/Makefile create mode 100644 src/BiInts/NEEDED_MODULES create mode 100644 src/BiInts/README.rst create mode 100644 src/BiInts/ao_bi_integrals.irp.f create mode 100644 src/BiInts/bi_integrals.ezfio_config create mode 100644 src/BiInts/gauss_legendre.irp.f create mode 100644 src/BiInts/map_integrals.irp.f create mode 100644 src/BiInts/mo_bi_integrals.irp.f create mode 100644 src/BiInts/options.irp.f create mode 100644 src/BiInts/tests/Makefile create mode 100644 src/BiInts/tests/test_ao.ref create mode 100644 src/BiInts/tests/test_ao_map.irp.f create mode 100644 src/BiInts/tests/test_ao_map.ref create mode 100644 src/BiInts/tests/test_mo.irp.f create mode 100644 src/BiInts/tests/test_mo.ref create mode 100644 src/MOGuess/ASSUMPTIONS.rst create mode 100644 src/MOGuess/H_CORE_guess.irp.f create mode 100644 src/MOGuess/Makefile create mode 100644 src/MOGuess/README.rst create mode 100644 src/MOGuess/localize.irp.f create mode 100644 src/MOGuess/mo_ortho_lowdin.irp.f create mode 100644 src/MOGuess/pot_mo_ortho_lowdin_ints.irp.f create mode 100644 src/MOGuess/tests/Makefile create mode 100644 src/Utils/map_module.f90 diff --git a/src/BiInts/ASSUMPTIONS.rst b/src/BiInts/ASSUMPTIONS.rst new file mode 100644 index 00000000..e69de29b diff --git a/src/BiInts/Makefile b/src/BiInts/Makefile new file mode 100644 index 00000000..b2ea1de1 --- /dev/null +++ b/src/BiInts/Makefile @@ -0,0 +1,8 @@ +default: all + +# Define here all new external source files and objects.Don't forget to prefix the +# object files with IRPF90_temp/ +SRC= +OBJ= + +include $(QPACKAGE_ROOT)/src/Makefile.common diff --git a/src/BiInts/NEEDED_MODULES b/src/BiInts/NEEDED_MODULES new file mode 100644 index 00000000..5f709ce4 --- /dev/null +++ b/src/BiInts/NEEDED_MODULES @@ -0,0 +1 @@ +AOs Bitmask Electrons Ezfio_files MOs Nuclei Output Utils MonoInts diff --git a/src/BiInts/README.rst b/src/BiInts/README.rst new file mode 100644 index 00000000..eb7cefaf --- /dev/null +++ b/src/BiInts/README.rst @@ -0,0 +1,27 @@ +============= +BiInts Module +============= + +Here, all bi-electronic integrals (:math:`1/r_{12}`) are computed. As they have +4 indices and many are zero, they are stored in a map. +To fetch an AO integral, use the ``get_ao_bielec_integral(i,j,k,l,ao_integrals_map)`` + +function, and to fetch and MO integral, use +. + +Needed Modules +============== + +.. Do not edit this section. It was auto-generated from the +.. NEEDED_MODULES file. + +* `AOs `_ +* `Bitmask `_ +* `Electrons `_ +* `Ezfio_files `_ +* `MOs `_ +* `Nuclei `_ +* `Output `_ +* `Utils `_ +* `MonoInts `_ + diff --git a/src/BiInts/ao_bi_integrals.irp.f b/src/BiInts/ao_bi_integrals.irp.f new file mode 100644 index 00000000..904122d6 --- /dev/null +++ b/src/BiInts/ao_bi_integrals.irp.f @@ -0,0 +1,1016 @@ +double precision function ao_bielec_integral(i,j,k,l) + implicit none + BEGIN_DOC + ! integral of the AO basis or (ij|kl) + ! i(r1) j(r1) 1/r12 k(r2) l(r2) + END_DOC + integer,intent(in) :: i,j,k,l + integer :: p,q,r,s + double precision :: I_center(3),J_center(3),K_center(3),L_center(3) + integer :: num_i,num_j,num_k,num_l,dim1,I_power(3),J_power(3),K_power(3),L_power(3) + double precision :: integral + include 'include/constants.F' + double precision :: P_new(0:max_dim,3),P_center(3),fact_p,pp + double precision :: Q_new(0:max_dim,3),Q_center(3),fact_q,qq + integer :: iorder_p(3), iorder_q(3) + logical :: do_p, do_q, do_r + + PROVIDE all_utils + + dim1 = n_pt_max_integrals + + num_i = ao_nucl(i) + num_j = ao_nucl(j) + num_k = ao_nucl(k) + num_l = ao_nucl(l) + ao_bielec_integral = 0.d0 + double precision :: thresh + thresh = ao_integrals_threshold + thresh = 0.d0 + + if (num_i /= num_j .or. num_k /= num_l .or. num_j /= num_k)then + do p = 1, 3 + I_power(p) = ao_power(i,p) + J_power(p) = ao_power(j,p) + K_power(p) = ao_power(k,p) + L_power(p) = ao_power(l,p) + I_center(p) = nucl_coord(num_i,p) + J_center(p) = nucl_coord(num_j,p) + K_center(p) = nucl_coord(num_k,p) + L_center(p) = nucl_coord(num_l,p) + enddo + + do_p = .True. + do p = 1, ao_prim_num(i) + double precision :: coef1 + if (.not.do_p) exit + do_p = .False. + do_q = .True. + coef1 = ao_coef_transp(p,i) + do q = 1, ao_prim_num(j) + double precision :: coef2 + if (.not.do_q) exit + do_q = .False. + do_r = .True. + coef2 = coef1*ao_coef_transp(q,j) + double precision :: p_inv,q_inv + call give_explicit_poly_and_gaussian(P_new,P_center,pp,fact_p,iorder_p,& + ao_expo_transp(p,i),ao_expo_transp(q,j), & + I_power,J_power,I_center,J_center,dim1) + if (abs(fact_p) < 1.d-8) then + exit + endif + p_inv = 1.d0/pp + do r = 1, ao_prim_num(k) + double precision :: coef3 + if (.not.do_r) exit + do_r = .False. + coef3 = coef2*ao_coef_transp(r,k) + do s = 1, ao_prim_num(l) + double precision :: coef4 + coef4 = coef3*ao_coef_transp(s,l) + double precision :: general_primitive_integral + call give_explicit_poly_and_gaussian(Q_new,Q_center,qq,fact_q,iorder_q,& + ao_expo_transp(r,k),ao_expo_transp(s,l), & + K_power,L_power,K_center,L_center,dim1) + if (abs(fact_p*fact_q) < 1.d-8) then + exit + endif + q_inv = 1.d0/qq + integral = general_primitive_integral(dim1, & + P_new,P_center,fact_p,pp,p_inv,iorder_p, & + Q_new,Q_center,fact_q,qq,q_inv,iorder_q) + ao_bielec_integral += coef4 * integral + if (abs(integral) < thresh) then + exit + endif + do_r = .True. + do_q = .True. + do_p = .True. + enddo ! s + enddo ! r + enddo ! q + enddo ! p + + else + + do p = 1, 3 + I_power(p) = ao_power(i,p) + J_power(p) = ao_power(j,p) + K_power(p) = ao_power(k,p) + L_power(p) = ao_power(l,p) + enddo + double precision :: ERI + do_p = .True. + do p = 1, ao_prim_num(i) + if (.not.do_p) exit + do_p = .False. + do_q = .True. + coef1 = ao_coef_transp(p,i) + do q = 1, ao_prim_num(j) + if (.not.do_q) exit + do_q = .False. + do_r = .True. + coef2 = coef1*ao_coef_transp(q,j) + do r = 1, ao_prim_num(k) + if (.not.do_r) exit + do_r = .False. + coef3 = coef2*ao_coef_transp(r,k) + do s = 1, ao_prim_num(l) + coef4 = coef3*ao_coef_transp(s,l) + integral = ERI( & + ao_expo_transp(p,i),ao_expo_transp(q,j),ao_expo_transp(r,k),ao_expo_transp(s,l),& + I_power(1),J_power(1),K_power(1),L_power(1), & + I_power(2),J_power(2),K_power(2),L_power(2), & + I_power(3),J_power(3),K_power(3),L_power(3)) + ao_bielec_integral += coef4 * integral + if (abs(integral) < thresh) then + exit + endif + do_r = .True. + do_q = .True. + do_p = .True. + enddo ! s + enddo ! r + enddo ! q + enddo ! p + + endif + +end + + + + + + + +subroutine compute_ao_bielec_integrals(j,k,l,sze,buffer_value) + implicit none + use map_module + + BEGIN_DOC + ! Compute AO 1/r12 integrals for all i and fixed j,k,l + END_DOC + + integer, intent(in) :: j,k,l,sze + real(integral_kind), intent(out) :: buffer_value(sze) + double precision :: ao_bielec_integral + double precision :: thresh + thresh = ao_integrals_threshold + + integer :: n_centers, i + integer*1 :: center_count(nucl_num) + + PROVIDE gauleg_t2 ao_nucl all_utils + + if (ao_overlap_abs(j,l) < thresh) then + buffer_value = 0. + return + endif + + center_count = 0 + + do i = 1, ao_num + if (ao_overlap_abs(i,k)*ao_overlap_abs(j,l) < thresh) then + buffer_value(i) = 0. + cycle + endif + !DIR$ FORCEINLINE + buffer_value(i) = ao_bielec_integral(i,k,j,l) + enddo + +end + + + + + + + + +BEGIN_PROVIDER [ logical, ao_bielec_integrals_in_map ] + implicit none + use map_module + BEGIN_DOC + ! Map of Atomic integrals + END_DOC + + integer :: i,j,k,l + double precision :: ao_bielec_integral,cpu_1,cpu_2, wall_1, wall_2 + double precision :: integral + double precision :: thresh + thresh = ao_integrals_threshold + + ! For integrals file + integer*8,allocatable :: buffer_i(:) + integer,parameter :: size_buffer = 1024*64 + real(integral_kind),allocatable :: buffer_value(:) + integer(omp_lock_kind) :: lock + + integer :: n_integrals, n_centers + integer :: jl_pairs(2,ao_num*(ao_num+1)/2), kk, m, j1, i1, lmax + integer*1 :: center_count(nucl_num) + + PROVIDE gauleg_t2 ao_integrals_map all_utils + integral = ao_bielec_integral(1,1,1,1) + + real :: map_mb + if (read_ao_integrals) then + integer :: load_ao_integrals + if (load_ao_integrals(trim(ezfio_filename)//'/work/ao_integrals.bin') == 0) then + write(output_BiInts,*) 'AO integrals provided' + ao_bielec_integrals_in_map = .True. + return + endif + endif + + kk=1 + do l = 1, ao_num ! r2 + do j = 1, l ! r2 + jl_pairs(1,kk) = j + jl_pairs(2,kk) = l + kk += 1 + enddo + enddo + + PROVIDE ao_nucl + call omp_init_lock(lock) + lmax = ao_num*(ao_num+1)/2 + write(output_BiInts,*) 'providing the AO integrals' + call wall_time(wall_1) + call cpu_time(cpu_1) + !$OMP PARALLEL PRIVATE(i,j,k,l,kk, & + !$OMP integral,buffer_i,buffer_value,n_integrals, & + !$OMP cpu_2,wall_2,i1,j1,center_count) & + !$OMP DEFAULT(NONE) & + !$OMP SHARED (ao_num, jl_pairs, ao_integrals_map,thresh, & + !$OMP cpu_1,wall_1,lock, lmax,n_centers,ao_nucl, & + !$OMP ao_overlap_abs,ao_overlap,output_BiInts) + + allocate(buffer_i(size_buffer)) + allocate(buffer_value(size_buffer)) + n_integrals = 0 + center_count = 0 + + !$OMP DO SCHEDULE(dynamic) + do kk=1,lmax + j = jl_pairs(1,kk) + l = jl_pairs(2,kk) + j1 = j+ishft(l*l-l,-1) + if (ao_overlap_abs(j,l) < thresh) then + cycle + endif + do k = 1, ao_num ! r1 + i1 = ishft(k*k-k,-1) + if (i1 > j1) then + exit + endif + do i = 1, k + i1 += 1 + if (i1 > j1) then + exit + endif + if (ao_overlap_abs(i,k)*ao_overlap_abs(j,l) < thresh) then + cycle + endif + !DIR$ FORCEINLINE + integral = ao_bielec_integral(i,k,j,l) + if (abs(integral) < thresh) then + cycle + endif + n_integrals += 1 + !DIR$ FORCEINLINE + call bielec_integrals_index(i,j,k,l,buffer_i(n_integrals)) + buffer_value(n_integrals) = integral + if (n_integrals > 1024 ) then + if (omp_test_lock(lock)) then + call insert_into_ao_integrals_map(n_integrals,buffer_i,buffer_value) + n_integrals = 0 + call omp_unset_lock(lock) + endif + endif + if (n_integrals == size(buffer_i)) then + call insert_into_ao_integrals_map(n_integrals,buffer_i,buffer_value) + n_integrals = 0 + endif + enddo + enddo + call wall_time(wall_2) + write(output_BiInts,*) 100.*float(kk)/float(lmax), '% in ', wall_2-wall_1, 's' + enddo + !$OMP END DO NOWAIT + call insert_into_ao_integrals_map(n_integrals,buffer_i,buffer_value) + deallocate(buffer_i) + deallocate(buffer_value) + !$OMP END PARALLEL + call omp_destroy_lock(lock) + write(output_BiInts,*) 'Sorting the map' + call map_sort(ao_integrals_map) + call cpu_time(cpu_2) + call wall_time(wall_2) + integer*8 :: get_ao_map_size, ao_map_size + ao_map_size = get_ao_map_size() + + write(output_BiInts,*) 'AO integrals provided:' + write(output_BiInts,*) ' Size of AO map : ', map_mb(ao_integrals_map) ,'MB' + write(output_BiInts,*) ' Number of AO integrals :', ao_map_size + write(output_BiInts,*) ' cpu time :',cpu_2 - cpu_1, 's' + write(output_BiInts,*) ' wall time :',wall_2 - wall_1, 's ( x ', (cpu_2-cpu_1)/(wall_2-wall_1+tiny(1.d0)), ' )' + + ao_bielec_integrals_in_map = .True. + if (write_ao_integrals) then + call dump_ao_integrals(trim(ezfio_filename)//'/work/ao_integrals.bin') + endif + +END_PROVIDER + + + +BEGIN_PROVIDER [ double precision, ao_bielec_integral_schwartz,(ao_num,ao_num) ] + implicit none + BEGIN_DOC + ! Needed to compuet Schwartz inequalities + END_DOC + + integer :: i,k + double precision :: ao_bielec_integral,cpu_1,cpu_2, wall_1, wall_2 + + PROVIDE gauleg_t2 ao_integrals_map + ao_bielec_integral_schwartz(1,1) = ao_bielec_integral(1,1,1,1) + !$OMP PARALLEL DO PRIVATE(i,k) & + !$OMP DEFAULT(NONE) & + !$OMP SHARED (ao_num,ao_bielec_integral_schwartz) & + !$OMP SCHEDULE(dynamic) + do i=1,ao_num + do k=1,i + ao_bielec_integral_schwartz(i,k) = dsqrt(ao_bielec_integral(i,k,i,k)) + ao_bielec_integral_schwartz(k,i) = ao_bielec_integral_schwartz(i,k) + enddo + enddo + !$OMP END PARALLEL DO + +END_PROVIDER + + +double precision function general_primitive_integral(dim, & + P_new,P_center,fact_p,p,p_inv,iorder_p, & + Q_new,Q_center,fact_q,q,q_inv,iorder_q) + implicit none + BEGIN_DOC + ! Computes the integral where p,q,r,s are Gaussian primitives + END_DOC + integer,intent(in) :: dim + include 'constants.F' + double precision, intent(in) :: P_new(0:max_dim,3),P_center(3),fact_p,p,p_inv + double precision, intent(in) :: Q_new(0:max_dim,3),Q_center(3),fact_q,q,q_inv + integer, intent(in) :: iorder_p(3) + integer, intent(in) :: iorder_q(3) + + double precision :: r_cut,gama_r_cut,rho,dist + double precision :: dx(0:max_dim),Ix_pol(0:max_dim),dy(0:max_dim),Iy_pol(0:max_dim),dz(0:max_dim),Iz_pol(0:max_dim) + integer :: n_Ix,n_Iy,n_Iz,nx,ny,nz + double precision :: bla + integer :: ix,iy,iz,jx,jy,jz,i + double precision :: a,b,c,d,e,f,accu,pq,const + double precision :: pq_inv, p10_1, p10_2, p01_1, p01_2,pq_inv_2 + integer :: n_pt_tmp,n_pt_out, iorder + double precision :: d1(0:max_dim),d_poly(0:max_dim),rint,d1_screened(0:max_dim) + + general_primitive_integral = 0.d0 + + !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: dx,Ix_pol,dy,Iy_pol,dz,Iz_pol + !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: d1, d_poly + + ! Gaussian Product + ! ---------------- + + pq = p_inv*0.5d0*q_inv + pq_inv = 0.5d0/(p+q) + p10_1 = q*pq ! 1/(2p) + p01_1 = p*pq ! 1/(2q) + pq_inv_2 = pq_inv+pq_inv + p10_2 = pq_inv_2 * p10_1*q !0.5d0*q/(pq + p*p) + p01_2 = pq_inv_2 * p01_1*p !0.5d0*p/(q*q + pq) + + + accu = 0.d0 + iorder = iorder_p(1)+iorder_q(1)+iorder_p(1)+iorder_q(1) + !DIR$ VECTOR ALIGNED + do ix=0,iorder + Ix_pol(ix) = 0.d0 + enddo + n_Ix = 0 + do ix = 0, iorder_p(1) + if (abs(P_new(ix,1)) < 1.d-8) cycle + a = P_new(ix,1) + do jx = 0, iorder_q(1) + d = a*Q_new(jx,1) + if (abs(d) < 1.d-8) cycle + !DEC$ FORCEINLINE + call give_polynom_mult_center_x(P_center(1),Q_center(1),ix,jx,p,q,iorder,pq_inv,pq_inv_2,p10_1,p01_1,p10_2,p01_2,dx,nx) + !DEC$ FORCEINLINE + call add_poly_multiply(dx,nx,d,Ix_pol,n_Ix) + enddo + enddo + if (n_Ix == -1) then + return + endif + iorder = iorder_p(2)+iorder_q(2)+iorder_p(2)+iorder_q(2) + !DIR$ VECTOR ALIGNED + do ix=0, iorder + Iy_pol(ix) = 0.d0 + enddo + n_Iy = 0 + do iy = 0, iorder_p(2) + if (abs(P_new(iy,2)) > 1.d-8) then + b = P_new(iy,2) + do jy = 0, iorder_q(2) + e = b*Q_new(jy,2) + if (abs(e) < 1.d-8) cycle + !DEC$ FORCEINLINE + call give_polynom_mult_center_x(P_center(2),Q_center(2),iy,jy,p,q,iorder,pq_inv,pq_inv_2,p10_1,p01_1,p10_2,p01_2,dy,ny) + !DEC$ FORCEINLINE + call add_poly_multiply(dy,ny,e,Iy_pol,n_Iy) + enddo + endif + enddo + if (n_Iy == -1) then + return + endif + + iorder = iorder_p(3)+iorder_q(3)+iorder_p(3)+iorder_q(3) + do ix=0,iorder + Iz_pol(ix) = 0.d0 + enddo + n_Iz = 0 + do iz = 0, iorder_p(3) + if (abs(P_new(iz,3)) > 1.d-8) then + c = P_new(iz,3) + do jz = 0, iorder_q(3) + f = c*Q_new(jz,3) + if (abs(f) < 1.d-8) cycle + !DEC$ FORCEINLINE + call give_polynom_mult_center_x(P_center(3),Q_center(3),iz,jz,p,q,iorder,pq_inv,pq_inv_2,p10_1,p01_1,p10_2,p01_2,dz,nz) + !DEC$ FORCEINLINE + call add_poly_multiply(dz,nz,f,Iz_pol,n_Iz) + enddo + endif + enddo + if (n_Iz == -1) then + return + endif + + rho = p*q *pq_inv_2 + dist = (P_center(1) - Q_center(1))*(P_center(1) - Q_center(1)) + & + (P_center(2) - Q_center(2))*(P_center(2) - Q_center(2)) + & + (P_center(3) - Q_center(3))*(P_center(3) - Q_center(3)) + const = dist*rho + + n_pt_tmp = n_Ix+n_Iy + do i=0,n_pt_tmp + d_poly(i)=0.d0 + enddo + + !DEC$ FORCEINLINE + call multiply_poly(Ix_pol,n_Ix,Iy_pol,n_Iy,d_poly,n_pt_tmp) + n_pt_out = n_pt_tmp+n_Iz + do i=0,n_pt_out + d1(i)=0.d0 + enddo + + !DEC$ FORCEINLINE + call multiply_poly(d_poly ,n_pt_tmp ,Iz_pol,n_Iz,d1,n_pt_out) + double precision :: rint_sum + accu = accu + rint_sum(n_pt_out,const,d1) + + general_primitive_integral = fact_p * fact_q * accu *pi_5_2*p_inv*q_inv/dsqrt(p+q) +end + + +double precision function ERI(alpha,beta,delta,gama,a_x,b_x,c_x,d_x,a_y,b_y,c_y,d_y,a_z,b_z,c_z,d_z) + implicit none + BEGIN_DOC + ! ATOMIC PRIMTIVE bielectronic integral between the 4 primitives :: + ! primitive_1 = x1**(a_x) y1**(a_y) z1**(a_z) exp(-alpha * r1**2) + ! primitive_2 = x1**(b_x) y1**(b_y) z1**(b_z) exp(- beta * r1**2) + ! primitive_3 = x2**(c_x) y2**(c_y) z2**(c_z) exp(-delta * r2**2) + ! primitive_4 = x2**(d_x) y2**(d_y) z2**(d_z) exp(- gama * r2**2) + END_DOC + double precision, intent(in) :: delta,gama,alpha,beta + integer, intent(in) :: a_x,b_x,c_x,d_x,a_y,b_y,c_y,d_y,a_z,b_z,c_z,d_z + integer :: a_x_2,b_x_2,c_x_2,d_x_2,a_y_2,b_y_2,c_y_2,d_y_2,a_z_2,b_z_2,c_z_2,d_z_2 + integer :: i,j,k,l,n_pt + integer :: n_pt_sup + double precision :: p,q,denom,coeff + double precision :: I_f + include 'constants.F' + if(iand(a_x+b_x+c_x+d_x,1).eq.1.or.iand(a_y+b_y+c_y+d_y,1).eq.1.or.iand(a_z+b_z+c_z+d_z,1).eq.1)then + ERI = 0.d0 + return + else + + ASSERT (alpha >= 0.d0) + ASSERT (beta >= 0.d0) + ASSERT (delta >= 0.d0) + ASSERT (gama >= 0.d0) + p = alpha + beta + q = delta + gama + ASSERT (p+q >= 0.d0) + n_pt = n_pt_sup(a_x,b_x,c_x,d_x,a_y,b_y,c_y,d_y,a_z,b_z,c_z,d_z) + coeff = pi_5_2 / (p * q * dsqrt(p+q)) + + if (n_pt == 0) then + ERI = coeff + return + endif + + call integrale_new(I_f,a_x,b_x,c_x,d_x,a_y,b_y,c_y,d_y,a_z,b_z,c_z,d_z,p,q,n_pt) + + ERI = I_f * coeff + endif +end + + +subroutine integrale_new(I_f,a_x,b_x,c_x,d_x,a_y,b_y,c_y,d_y,a_z,b_z,c_z,d_z,p,q,n_pt) + BEGIN_DOC + ! calculate the integral of the polynom :: + ! I_x1(a_x+b_x, c_x+d_x,p,q) * I_x1(a_y+b_y, c_y+d_y,p,q) * I_x1(a_z+b_z, c_z+d_z,p,q) + ! between ( 0 ; 1) + END_DOC + + + implicit none + double precision :: p,q + double precision :: I_x1_new, I_x2_new + integer :: a_x,b_x,c_x,d_x,a_y,b_y,c_y,d_y,a_z,b_z,c_z,d_z + integer :: i, n_iter, n_pt, j + double precision :: accu, I_f, pq_inv, p10_1, p10_2, p01_1, p01_2 + + j = ishft(n_pt,-1) + accu = 0.d0 + !if(iand(a_x+b_x+c_x+d_x,1).eq.1.or.iand(a_y+b_y+c_y+d_y,1).eq.1.or.iand(a_z+b_z+c_z+d_z,1).eq.1)then + ! I_f = 0.d0 + !else + ASSERT (n_pt > 1) + pq_inv = 0.5d0/(p+q) + pq_inv_2 = pq_inv + pq_inv + p10_1 = 0.5d0/p + p01_1 = 0.5d0/q + p10_2 = 0.5d0 * q /(p * q + p * p) + p01_2 = 0.5d0 * p /(q * q + q * p) + double precision :: B10, B01, B00, flush_buffer,flush_buffer_one,rho,pq_inv_2 + do i = 1,n_pt + B10 = p10_1 - gauleg_t2(i,j)* p10_2 + B01 = p01_1 - gauleg_t2(i,j)* p01_2 + B00 = gauleg_t2(i,j)*pq_inv + flush_buffer = 0.d0 + flush_buffer_one = 1.d0 + accu += gauleg_w(i,j)*I_x1_new(a_x+b_x,c_x+d_x,B10,B01,B00,flush_buffer)& + *I_x1_new(a_y+b_y,c_y+d_y,B10,B01,B00,flush_buffer_one) & + *I_x1_new(a_z+b_z,c_z+d_z,B10,B01,B00,flush_buffer_one) + enddo + + I_f= accu + + !endif + +end + +recursive double precision function I_x1_new(a,c,B_10,B_01,B_00,I_0000) result(res) + BEGIN_DOC + ! recursive function involved in the bielectronic integral + END_DOC + implicit none + integer, intent(in) :: a,c + double precision, intent(in) :: B_10,B_01,B_00 + double precision, intent(in) :: I_0000 + double precision :: I_x2_new + + + if(c<0)then + res = 0.d0 + else if (a==0) then + res = I_x2_new(c,B_10,B_01,B_00,I_0000) + else if (a==1) then + res = c * B_00 * I_x2_new(c-1,B_10,B_01,B_00,I_0000) + else + res = (a-1) * B_10 * I_x1_new(a-2,c,B_10,B_01,B_00,I_0000) & + &+ c * B_00 * I_x1_new(a-1,c-1,B_10,B_01,B_00,I_0000) + endif +end + +recursive double precision function I_x2_new(c,B_10,B_01,B_00,I_0000) result(res) + implicit none + BEGIN_DOC + ! recursive function involved in the bielectronic integral + END_DOC + integer, intent(in) :: c + double precision, intent(in) :: B_10,B_01,B_00,I_0000 + double precision :: I_x1_new + + if(c==1)then + res = 0.d0 + elseif(c==0) then + res = 1.d0 + else + res = (c-1) * B_01 * I_x1_new(0,c-2,B_10,B_01,B_00,I_0000) + endif +end + + +integer function n_pt_sup(a_x,b_x,c_x,d_x,a_y,b_y,c_y,d_y,a_z,b_z,c_z,d_z) + implicit none + BEGIN_DOC + ! Returns the upper boundary of the degree of the polynom involved in the + ! bielctronic integral : + ! Ix(a_x,b_x,c_x,d_x) * Iy(a_y,b_y,c_y,d_y) * Iz(a_z,b_z,c_z,d_z) + END_DOC + integer :: a_x,b_x,c_x,d_x,a_y,b_y,c_y,d_y,a_z,b_z,c_z,d_z + n_pt_sup = 2 * ( (a_x+b_x+c_x+d_x) + (a_y+b_y+c_y+d_y) + (a_z+b_z+c_z+d_z) ) +end + + + + +subroutine give_polynom_mult_center_x(P_center,Q_center,a_x,d_x,p,q,n_pt_in,pq_inv,pq_inv_2,p10_1,p01_1,p10_2,p01_2,d,n_pt_out) + implicit none + BEGIN_DOC + ! subroutine that returns the explicit polynom in term of the "t" + ! variable of the following polynomw : + ! I_x1(a_x, d_x,p,q) * I_x1(a_y, d_y,p,q) * I_x1(a_z, d_z,p,q) + END_DOC + integer, intent(in) :: n_pt_in + integer,intent(out) :: n_pt_out + integer, intent(in) :: a_x,d_x + double precision, intent(in) :: P_center, Q_center + double precision, intent(in) :: p,q,pq_inv,p10_1,p01_1,p10_2,p01_2,pq_inv_2 + include 'constants.F' + double precision,intent(out) :: d(0:max_dim) + double precision :: accu + accu = 0.d0 + ASSERT (n_pt_in >= 0) + ! pq_inv = 0.5d0/(p+q) + ! pq_inv_2 = 1.d0/(p+q) + ! p10_1 = 0.5d0/p + ! p01_1 = 0.5d0/q + ! p10_2 = 0.5d0 * q /(p * q + p * p) + ! p01_2 = 0.5d0 * p /(q * q + q * p) + double precision :: B10(0:2), B01(0:2), B00(0:2),C00(0:2),D00(0:2) + B10(0) = p10_1 + B10(1) = 0.d0 + B10(2) = - p10_2 + ! B10 = p01_1 - t**2 * p10_2 + B01(0) = p01_1 + B01(1) = 0.d0 + B01(2) = - p01_2 + ! B01 = p01_1- t**2 * pq_inv + B00(0) = 0.d0 + B00(1) = 0.d0 + B00(2) = pq_inv + ! B00 = t**2 * pq_inv + do i = 0,n_pt_in + d(i) = 0.d0 + enddo + integer :: n_pt1,dim,i + n_pt1 = n_pt_in + ! C00 = -q/(p+q)*(Px-Qx) * t^2 + C00(0) = 0.d0 + C00(1) = 0.d0 + C00(2) = -q*(P_center-Q_center) * pq_inv_2 + ! D00 = -p/(p+q)*(Px-Qx) * t^2 + D00(0) = 0.d0 + D00(1) = 0.d0 + D00(2) = -p*(Q_center-P_center) * pq_inv_2 + !D00(2) = -p*(Q_center(1)-P_center(1)) /(p+q) + !DIR$ FORCEINLINE + call I_x1_pol_mult(a_x,d_x,B10,B01,B00,C00,D00,d,n_pt1,n_pt_in) + n_pt_out = n_pt1 + if(n_pt1<0)then + n_pt_out = -1 + do i = 0,n_pt_in + d(i) = 0.d0 + enddo + return + endif + +end + +subroutine I_x1_pol_mult(a,c,B_10,B_01,B_00,C_00,D_00,d,nd,n_pt_in) + implicit none + BEGIN_DOC + ! recursive function involved in the bielectronic integral + END_DOC + integer , intent(in) :: n_pt_in + include 'constants.F' + double precision,intent(inout) :: d(0:max_dim) + integer,intent(inout) :: nd + integer, intent(in) :: a,c + double precision, intent(in) :: B_10(0:2),B_01(0:2),B_00(0:2),C_00(0:2),D_00(0:2) + if( (c>=0).and.(nd>=0) )then + + if (a==1) then + call I_x1_pol_mult_a1(c,B_10,B_01,B_00,C_00,D_00,d,nd,n_pt_in) + else if (a==2) then + call I_x1_pol_mult_a2(c,B_10,B_01,B_00,C_00,D_00,d,nd,n_pt_in) + else if (a>2) then + call I_x1_pol_mult_recurs(a,c,B_10,B_01,B_00,C_00,D_00,d,nd,n_pt_in) + else ! a == 0 + + if( c==0 )then + nd = 0 + d(0) = 1.d0 + return + endif + + call I_x2_pol_mult(c,B_10,B_01,B_00,C_00,D_00,d,nd,n_pt_in) + endif + else + nd = -1 + endif +end + +recursive subroutine I_x1_pol_mult_recurs(a,c,B_10,B_01,B_00,C_00,D_00,d,nd,n_pt_in) + implicit none + BEGIN_DOC + ! recursive function involved in the bielectronic integral + END_DOC + integer , intent(in) :: n_pt_in + include 'constants.F' + double precision,intent(inout) :: d(0:max_dim) + integer,intent(inout) :: nd + integer, intent(in) :: a,c + double precision, intent(in) :: B_10(0:2),B_01(0:2),B_00(0:2),C_00(0:2),D_00(0:2) + double precision :: X(0:max_dim) + double precision :: Y(0:max_dim) + !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: X,Y + integer :: nx, ix,iy,ny + + ASSERT (a>2) + !DIR$ VECTOR ALIGNED + !DIR$ LOOP COUNT(8) + do ix=0,n_pt_in + X(ix) = 0.d0 + enddo + nx = 0 + if (a==3) then + call I_x1_pol_mult_a1(c,B_10,B_01,B_00,C_00,D_00,X,nx,n_pt_in) + else if (a==4) then + call I_x1_pol_mult_a2(c,B_10,B_01,B_00,C_00,D_00,X,nx,n_pt_in) + else + ASSERT (a>=5) + call I_x1_pol_mult_recurs(a-2,c,B_10,B_01,B_00,C_00,D_00,X,nx,n_pt_in) + endif + + !DIR$ VECTOR ALIGNED + !DIR$ LOOP COUNT(8) + do ix=0,nx + X(ix) *= dble(a-1) + enddo + + !DEC$ FORCEINLINE + call multiply_poly(X,nx,B_10,2,d,nd) + + nx = nd + !DIR$ VECTOR ALIGNED + !DIR$ LOOP COUNT(8) + do ix=0,n_pt_in + X(ix) = 0.d0 + enddo + + if (c>0) then + if (a==3) then + call I_x1_pol_mult_a2(c-1,B_10,B_01,B_00,C_00,D_00,X,nx,n_pt_in) + else + ASSERT(a >= 4) + call I_x1_pol_mult_recurs(a-1,c-1,B_10,B_01,B_00,C_00,D_00,X,nx,n_pt_in) + endif + if (c>1) then + !DIR$ VECTOR ALIGNED + !DIR$ LOOP COUNT(8) + do ix=0,nx + X(ix) *= c + enddo + endif + !DEC$ FORCEINLINE + call multiply_poly(X,nx,B_00,2,d,nd) + endif + + ny=0 + + !DIR$ VECTOR ALIGNED + !DIR$ LOOP COUNT(8) + do ix=0,n_pt_in + Y(ix) = 0.d0 + enddo + ASSERT(a > 2) + if (a==3) then + call I_x1_pol_mult_a2(c,B_10,B_01,B_00,C_00,D_00,Y,ny,n_pt_in) + else + ASSERT(a >= 4) + call I_x1_pol_mult_recurs(a-1,c,B_10,B_01,B_00,C_00,D_00,Y,ny,n_pt_in) + endif + + !DEC$ FORCEINLINE + call multiply_poly(Y,ny,C_00,2,d,nd) + +end + +recursive subroutine I_x1_pol_mult_a1(c,B_10,B_01,B_00,C_00,D_00,d,nd,n_pt_in) + implicit none + BEGIN_DOC + ! recursive function involved in the bielectronic integral + END_DOC + integer , intent(in) :: n_pt_in + include 'constants.F' + double precision,intent(inout) :: d(0:max_dim) + integer,intent(inout) :: nd + integer, intent(in) :: c + double precision, intent(in) :: B_10(0:2),B_01(0:2),B_00(0:2),C_00(0:2),D_00(0:2) + double precision :: X(0:max_dim) + double precision :: Y(0:max_dim) + !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: X,Y + integer :: nx, ix,iy,ny + + if( (c<0).or.(nd<0) )then + nd = -1 + return + endif + + nx = nd + !DIR$ VECTOR ALIGNED + !DIR$ LOOP COUNT(8) + do ix=0,n_pt_in + X(ix) = 0.d0 + enddo + call I_x2_pol_mult(c-1,B_10,B_01,B_00,C_00,D_00,X,nx,n_pt_in) + + if (c>1) then + !DIR$ VECTOR ALIGNED + !DIR$ LOOP COUNT(8) + do ix=0,nx + X(ix) *= dble(c) + enddo + endif + + !DEC$ FORCEINLINE + call multiply_poly(X,nx,B_00,2,d,nd) + + ny=0 + + !DIR$ VECTOR ALIGNED + !DIR$ LOOP COUNT(8) + do ix=0,n_pt_in + Y(ix) = 0.d0 + enddo + call I_x2_pol_mult(c,B_10,B_01,B_00,C_00,D_00,Y,ny,n_pt_in) + + !DEC$ FORCEINLINE + call multiply_poly(Y,ny,C_00,2,d,nd) + +end + +recursive subroutine I_x1_pol_mult_a2(c,B_10,B_01,B_00,C_00,D_00,d,nd,n_pt_in) + implicit none + BEGIN_DOC + ! recursive function involved in the bielectronic integral + END_DOC + integer , intent(in) :: n_pt_in + include 'constants.F' + double precision,intent(inout) :: d(0:max_dim) + integer,intent(inout) :: nd + integer, intent(in) :: c + double precision, intent(in) :: B_10(0:2),B_01(0:2),B_00(0:2),C_00(0:2),D_00(0:2) + double precision :: X(0:max_dim) + double precision :: Y(0:max_dim) + !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: X,Y + integer :: nx, ix,iy,ny + + !DIR$ VECTOR ALIGNED + !DIR$ LOOP COUNT(8) + do ix=0,n_pt_in + X(ix) = 0.d0 + enddo + nx = 0 + call I_x2_pol_mult(c,B_10,B_01,B_00,C_00,D_00,X,nx,n_pt_in) + + !DEC$ FORCEINLINE + call multiply_poly(X,nx,B_10,2,d,nd) + + nx = nd + !DIR$ VECTOR ALIGNED + !DIR$ LOOP COUNT(8) + do ix=0,n_pt_in + X(ix) = 0.d0 + enddo + + !DIR$ FORCEINLINE + call I_x1_pol_mult_a1(c-1,B_10,B_01,B_00,C_00,D_00,X,nx,n_pt_in) + + if (c>1) then + !DIR$ VECTOR ALIGNED + !DIR$ LOOP COUNT(8) + do ix=0,nx + X(ix) *= dble(c) + enddo + endif + + !DEC$ FORCEINLINE + call multiply_poly(X,nx,B_00,2,d,nd) + + ny=0 + !DIR$ VECTOR ALIGNED + !DIR$ LOOP COUNT(8) + do ix=0,n_pt_in + Y(ix) = 0.d0 + enddo + !DEC$ FORCEINLINE + call I_x1_pol_mult_a1(c,B_10,B_01,B_00,C_00,D_00,Y,ny,n_pt_in) + + !DEC$ FORCEINLINE + call multiply_poly(Y,ny,C_00,2,d,nd) + +end + +recursive subroutine I_x2_pol_mult(c,B_10,B_01,B_00,C_00,D_00,d,nd,dim) + implicit none + BEGIN_DOC + ! recursive function involved in the bielectronic integral + END_DOC + integer , intent(in) :: dim + include 'constants.F' + double precision :: d(0:max_dim) + integer,intent(inout) :: nd + integer, intent(in) :: c + double precision, intent(in) :: B_10(0:2),B_01(0:2),B_00(0:2),C_00(0:2),D_00(0:2) + integer :: nx, ix,ny + double precision :: X(0:max_dim),Y(0:max_dim) + !DEC$ ATTRIBUTES ALIGN : $IRP_ALIGN :: X, Y + integer :: i + + select case (c) + case (0) + nd = 0 + d(0) = 1.d0 + return + + case (:-1) + nd = -1 + return + + case (1) + nd = 2 + d(0) = D_00(0) + d(1) = D_00(1) + d(2) = D_00(2) + return + + case (2) + nd = 2 + d(0) = B_01(0) + d(1) = B_01(1) + d(2) = B_01(2) + + ny = 2 + Y(0) = D_00(0) + Y(1) = D_00(1) + Y(2) = D_00(2) + + !DEC$ FORCEINLINE + call multiply_poly(Y,ny,D_00,2,d,nd) + return + + case default + + !DIR$ VECTOR ALIGNED + !DIR$ LOOP COUNT(6) + do ix=0,c+c + X(ix) = 0.d0 + enddo + nx = 0 + call I_x2_pol_mult(c-2,B_10,B_01,B_00,C_00,D_00,X,nx,dim) + + !DIR$ VECTOR ALIGNED + !DIR$ LOOP COUNT(6) + do ix=0,nx + X(ix) *= dble(c-1) + enddo + + !DEC$ FORCEINLINE + call multiply_poly(X,nx,B_01,2,d,nd) + + ny = 0 + !DIR$ VECTOR ALIGNED + !DIR$ LOOP COUNT(6) + do ix=0,c+c + Y(ix) = 0.d0 + enddo + call I_x2_pol_mult(c-1,B_10,B_01,B_00,C_00,D_00,Y,ny,dim) + + !DEC$ FORCEINLINE + call multiply_poly(Y,ny,D_00,2,d,nd) + + end select +end + + diff --git a/src/BiInts/bi_integrals.ezfio_config b/src/BiInts/bi_integrals.ezfio_config new file mode 100644 index 00000000..c0a3989a --- /dev/null +++ b/src/BiInts/bi_integrals.ezfio_config @@ -0,0 +1,8 @@ +bielec_integrals + read_ao_integrals logical + read_mo_integrals logical + write_ao_integrals logical + write_mo_integrals logical + threshold_ao double precision + threshold_mo double precision + diff --git a/src/BiInts/gauss_legendre.irp.f b/src/BiInts/gauss_legendre.irp.f new file mode 100644 index 00000000..e7545d58 --- /dev/null +++ b/src/BiInts/gauss_legendre.irp.f @@ -0,0 +1,57 @@ + BEGIN_PROVIDER [ double precision, gauleg_t2, (n_pt_max_integrals,n_pt_max_integrals/2) ] +&BEGIN_PROVIDER [ double precision, gauleg_w, (n_pt_max_integrals,n_pt_max_integrals/2) ] + implicit none + BEGIN_DOC + ! t_w(i,1,k) = w(i) + ! t_w(i,2,k) = t(i) + END_DOC + integer :: i,j,l + l=0 + do i = 2,n_pt_max_integrals,2 + l = l+1 + call gauleg(0.d0,1.d0,gauleg_t2(1,l),gauleg_w(1,l),i) + do j=1,i + gauleg_t2(j,l) *= gauleg_t2(j,l) + enddo + enddo + +END_PROVIDER + +subroutine gauleg(x1,x2,x,w,n) + implicit none + BEGIN_DOC + ! Gauss-Legendre + END_DOC + integer, intent(in) :: n + double precision, intent(in) :: x1, x2 + double precision, intent (out) :: x(n),w(n) + double precision, parameter :: eps=3.d-14 + + integer :: m,i,j + double precision :: xm, xl, z, z1, p1, p2, p3, pp, dn + m=(n+1)/2 + xm=0.5d0*(x2+x1) + xl=0.5d0*(x2-x1) + dn = dble(n) + do i=1,m + z=dcos(3.141592654d0*(dble(i)-.25d0)/(dble(n)+.5d0)) + z1 = z+1.d0 + do while (dabs(z-z1) > eps) + p1=1.d0 + p2=0.d0 + do j=1,n + p3=p2 + p2=p1 + p1=(dble(j+j-1)*z*p2-dble(j-1)*p3)/j + enddo + pp=dn*(z*p1-p2)/(z*z-1.d0) + z1=z + z=z1-p1/pp + end do + x(i)=xm-xl*z + x(n+1-i)=xm+xl*z + w(i)=(xl+xl)/((1.d0-z*z)*pp*pp) + w(n+1-i)=w(i) + enddo +end + diff --git a/src/BiInts/map_integrals.irp.f b/src/BiInts/map_integrals.irp.f new file mode 100644 index 00000000..5d639702 --- /dev/null +++ b/src/BiInts/map_integrals.irp.f @@ -0,0 +1,334 @@ +use map_module + +!! AO Map +!! ====== + +BEGIN_PROVIDER [ type(map_type), ao_integrals_map ] + implicit none + BEGIN_DOC + ! AO integrals + END_DOC + integer*8 :: sze + call bielec_integrals_index(ao_num,ao_num,ao_num,ao_num,sze) + call map_init(ao_integrals_map,sze) + write(output_BiInts,*) 'AO map initialized' +END_PROVIDER + +subroutine bielec_integrals_index(i,j,k,l,i1) + implicit none + integer, intent(in) :: i,j,k,l + integer*8, intent(out) :: i1 + integer*8 :: p,q,r,s,i2 + p = min(i,k) + r = max(i,k) + p = p+ishft(r*r-r,-1) + q = min(j,l) + s = max(j,l) + q = q+ishft(s*s-s,-1) + i1 = min(p,q) + i2 = max(p,q) + i1 = i1+ishft(i2*i2-i2,-1) +end + +double precision function get_ao_bielec_integral(i,j,k,l,map) + use map_module + implicit none + BEGIN_DOC + ! Gets one AO bi-electronic integral from the AO map + END_DOC + integer, intent(in) :: i,j,k,l + integer*8 :: idx + type(map_type), intent(inout) :: map + real(integral_kind) :: tmp + PROVIDE ao_bielec_integrals_in_map + !DIR$ FORCEINLINE + call bielec_integrals_index(i,j,k,l,idx) + call map_get(map,idx,tmp) + get_ao_bielec_integral = tmp +end + + +subroutine get_ao_bielec_integrals(j,k,l,sze,out_val) + use map_module + BEGIN_DOC + ! Gets multiple AO bi-electronic integral from the AO map . + ! All i are retrieved for j,k,l fixed. + END_DOC + implicit none + integer, intent(in) :: j,k,l, sze + real(integral_kind), intent(out) :: out_val(sze) + + integer :: i + integer*8 :: hash + double precision :: thresh + PROVIDE ao_bielec_integrals_in_map + thresh = ao_integrals_threshold + + if (ao_overlap_abs(j,l) < thresh) then + out_val = 0.d0 + return + endif + + do i=1,sze + if (ao_overlap_abs(i,k)*ao_overlap_abs(j,l) < thresh ) then + out_val(i) = 0.d0 + else + !DIR$ FORCEINLINE + call bielec_integrals_index(i,j,k,l,hash) + call map_get(ao_integrals_map, hash, out_val(i)) + endif + enddo + +end + +subroutine get_ao_bielec_integrals_non_zero(j,k,l,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 + integer, intent(in) :: j,k,l, sze + real(integral_kind), intent(out) :: out_val(sze) + integer, intent(out) :: out_val_index(sze),non_zero_int + + integer :: i + integer*8 :: hash + double precision :: thresh,tmp + PROVIDE ao_bielec_integrals_in_map + thresh = ao_integrals_threshold + + if (ao_overlap_abs(j,l) < thresh) then + out_val = 0.d0 + return + endif + + non_zero_int = 0 + do i=1,sze + !DIR$ FORCEINLINE + call bielec_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(non_zero_int) = i + out_val(non_zero_int) = tmp + enddo + +end + + +integer*8 function get_ao_map_size() + implicit none + BEGIN_DOC + ! Returns the number of elements in the AO map + END_DOC + get_ao_map_size = ao_integrals_map % n_elements +end + +subroutine clear_ao_map + implicit none + BEGIN_DOC + ! Frees the memory of the AO map + END_DOC + call map_deinit(ao_integrals_map) + FREE ao_integrals_map +end + + + +!! MO Map +!! ====== + +BEGIN_PROVIDER [ type(map_type), mo_integrals_map ] + implicit none + BEGIN_DOC + ! MO integrals + END_DOC + integer*8 :: sze + call bielec_integrals_index(mo_tot_num,mo_tot_num,mo_tot_num,mo_tot_num,sze) + call map_init(mo_integrals_map,sze) + write(output_BiInts,*) 'MO map initialized' +END_PROVIDER + +subroutine insert_into_ao_integrals_map(n_integrals, & + buffer_i, buffer_values) + use map_module + implicit none + BEGIN_DOC + ! Create new entry into AO map + END_DOC + + integer, intent(in) :: n_integrals + integer*8, intent(inout) :: buffer_i(n_integrals) + real(integral_kind), intent(inout) :: buffer_values(n_integrals) + + call map_append(ao_integrals_map, buffer_i, buffer_values, n_integrals) +end + +subroutine insert_into_mo_integrals_map(n_integrals, & + buffer_i, buffer_values, thr) + use map_module + implicit none + + BEGIN_DOC + ! Create new entry into MO map, or accumulate in an existing entry + END_DOC + + integer, intent(in) :: n_integrals + integer*8, intent(inout) :: buffer_i(n_integrals) + real(integral_kind), intent(inout) :: buffer_values(n_integrals) + real(integral_kind), intent(in) :: thr + call map_update(mo_integrals_map, buffer_i, buffer_values, n_integrals, thr) +end + +double precision function get_mo_bielec_integral(i,j,k,l,map) + use map_module + implicit none + BEGIN_DOC + ! Returns one integral in the MO basis + END_DOC + integer, intent(in) :: i,j,k,l + integer*8 :: idx + type(map_type), intent(inout) :: map + real(integral_kind) :: tmp + PROVIDE mo_bielec_integrals_in_map + !DIR$ FORCEINLINE + call bielec_integrals_index(i,j,k,l,idx) + call map_get(map,idx,tmp) + get_mo_bielec_integral = dble(tmp) +end + +double precision function mo_bielec_integral(i,j,k,l) + implicit none + BEGIN_DOC + ! Returns one integral in the MO basis + END_DOC + integer, intent(in) :: i,j,k,l + double precision :: get_mo_bielec_integral + PROVIDE mo_bielec_integrals_in_map + mo_bielec_integral = get_mo_bielec_integral(i,j,k,l,mo_integrals_map) + return +end + +subroutine get_mo_bielec_integrals(j,k,l,sze,out_val,map) + use map_module + implicit none + BEGIN_DOC + ! Returns multiple integrals in the MO basis, all + ! i for j,k,l fixed. + END_DOC + integer, intent(in) :: j,k,l, sze + real(integral_kind), intent(out) :: out_val(sze) + type(map_type), intent(inout) :: map + integer :: i + integer*8 :: hash(sze) + PROVIDE mo_bielec_integrals_in_map + + do i=1,sze + !DIR$ FORCEINLINE + call bielec_integrals_index(i,j,k,l,hash(i)) + enddo + + call map_get_many(map, hash, out_val, sze) +end + +integer*8 function get_mo_map_size() + implicit none + BEGIN_DOC + ! Return the number of elements in the MO map + END_DOC + get_mo_map_size = mo_integrals_map % n_elements +end + +subroutine clear_mo_map + implicit none + BEGIN_DOC + ! Frees the memory of the MO map + END_DOC + call map_deinit(mo_integrals_map) + FREE mo_integrals_map +end + +BEGIN_TEMPLATE + +subroutine dump_$ao_integrals(filename) + use map_module + implicit none + BEGIN_DOC + ! Save to disk the $ao integrals + END_DOC + character*(*), intent(in) :: filename + integer(cache_key_kind), pointer :: key(:) + real(integral_kind), pointer :: val(:) + integer*8 :: i,j, n + open(unit=66,file=filename,FORM='unformatted') + write(66) integral_kind, key_kind + write(66) $ao_integrals_map%sorted, $ao_integrals_map%map_size, & + $ao_integrals_map%n_elements + do i=0_8,$ao_integrals_map%map_size + write(66) $ao_integrals_map%map(i)%sorted, $ao_integrals_map%map(i)%map_size,& + $ao_integrals_map%map(i)%n_elements + enddo + do i=0_8,$ao_integrals_map%map_size + key => $ao_integrals_map%map(i)%key + val => $ao_integrals_map%map(i)%value + n = $ao_integrals_map%map(i)%n_elements + write(66) (key(j), j=1,n), (val(j), j=1,n) + enddo + close(66) + +end + + +integer function load_$ao_integrals(filename) + implicit none + BEGIN_DOC + ! Read from disk the $ao integrals + END_DOC + character*(*), intent(in) :: filename + integer*8 :: i + integer(cache_key_kind), pointer :: key(:) + real(integral_kind), pointer :: val(:) + integer :: iknd, kknd, n, j + load_$ao_integrals = 1 + open(unit=66,file=filename,FORM='unformatted',STATUS='UNKNOWN') + read(66,err=98,end=98) iknd, kknd + if (iknd /= integral_kind) then + print *, 'Wrong integrals kind in file :', iknd + stop 1 + endif + if (kknd /= key_kind) then + print *, 'Wrong key kind in file :', kknd + stop 1 + endif + read(66,err=98,end=98) $ao_integrals_map%sorted, $ao_integrals_map%map_size,& + $ao_integrals_map%n_elements + do i=0_8, $ao_integrals_map%map_size + read(66,err=99,end=99) $ao_integrals_map%map(i)%sorted, & + $ao_integrals_map%map(i)%map_size, $ao_integrals_map%map(i)%n_elements + call cache_map_reallocate($ao_integrals_map%map(i),$ao_integrals_map%map(i)%map_size) + enddo + do i=0_8, $ao_integrals_map%map_size + key => $ao_integrals_map%map(i)%key + val => $ao_integrals_map%map(i)%value + n = $ao_integrals_map%map(i)%n_elements + read(66,err=99,end=99) (key(j), j=1,n), (val(j), j=1,n) + enddo + call map_sort($ao_integrals_map) + load_$ao_integrals = 0 + return + 99 continue + call map_deinit($ao_integrals_map) + FREE $ao_integrals_map + if (.True.) then + PROVIDE $ao_integrals_map + endif + stop 'Problem reading $ao_integrals_map file in work/' + 98 continue + +end + +SUBST [ ao_integrals_map, ao_integrals, ao_num , get_ao_bielec_integral ] +ao_integrals_map ; ao_integrals ; ao_num ; get_ao_bielec_integral ;; +mo_integrals_map ; mo_integrals ; n_act ; get_mo_bielec_integral ;; +END_TEMPLATE diff --git a/src/BiInts/mo_bi_integrals.irp.f b/src/BiInts/mo_bi_integrals.irp.f new file mode 100644 index 00000000..a0c77fe1 --- /dev/null +++ b/src/BiInts/mo_bi_integrals.irp.f @@ -0,0 +1,356 @@ +subroutine mo_bielec_integrals_index(i,j,k,l,i1) + implicit none + BEGIN_DOC + ! Computes an unique index for i,j,k,l integrals + END_DOC + integer, intent(in) :: i,j,k,l + integer*8, intent(out) :: i1 + integer*8 :: p,q,r,s,i2 + p = min(i,k) + r = max(i,k) + p = p+ishft(r*r-r,-1) + q = min(j,l) + s = max(j,l) + q = q+ishft(s*s-s,-1) + i1 = min(p,q) + i2 = max(p,q) + i1 = i1+ishft(i2*i2-i2,-1) +end + + +BEGIN_PROVIDER [ logical, mo_bielec_integrals_in_map ] + implicit none + + BEGIN_DOC + ! If True, the map of MO bielectronic integrals is provided + END_DOC + + mo_bielec_integrals_in_map = .True. + if (read_mo_integrals) then + integer :: load_mo_integrals + if (load_mo_integrals(trim(ezfio_filename)//'/work/mo_integrals.bin') == 0) then + write(output_BiInts,*) 'MO integrals provided' + return + endif + endif + + call add_integrals_to_map(full_ijkl_bitmask) +END_PROVIDER + +subroutine add_integrals_to_map(mask_ijkl) + use bitmasks + implicit none + + BEGIN_DOC + ! Adds integrals to tha MO map according to some bitmask + END_DOC + + integer(bit_kind), intent(in) :: mask_ijkl(N_int,4) + + integer :: i,j,k,l + integer :: i0,j0,k0,l0 + double precision :: c, cpu_1, cpu_2, wall_1, wall_2 + + integer, allocatable :: list_ijkl(:,:) + integer :: n_i, n_j, n_k, n_l + integer, allocatable :: bielec_tmp_0_idx(:) + real(integral_kind), allocatable :: bielec_tmp_0(:,:) + double precision, allocatable :: bielec_tmp_1(:) + double precision, allocatable :: bielec_tmp_2(:,:) + double precision, allocatable :: bielec_tmp_3(:,:,:) + !DEC$ ATTRIBUTES ALIGN : 64 :: bielec_tmp_1, bielec_tmp_2, bielec_tmp_3 + + integer :: n_integrals + integer :: size_buffer + integer*8,allocatable :: buffer_i(:) + real(integral_kind),allocatable :: buffer_value(:) + real :: map_mb + + integer :: i1,j1,k1,l1, ii1, kmax, l1_global + integer :: i2,i3,i4 + double precision,parameter :: thr_coef = 0.d0 + + PROVIDE N_int ao_bielec_integrals_in_map ao_integrals_map mo_coef mo_coef_transp + + + !Get list of MOs for i,j,k and l + !------------------------------- + + allocate(list_ijkl(mo_tot_num,4)) + call bitstring_to_list( mask_ijkl(1,1), list_ijkl(1,1), n_i, N_int ) + call bitstring_to_list( mask_ijkl(1,2), list_ijkl(1,2), n_j, N_int ) + call bitstring_to_list( mask_ijkl(1,3), list_ijkl(1,3), n_k, N_int ) + call bitstring_to_list( mask_ijkl(1,4), list_ijkl(1,4), n_l, N_int ) + + + + + l1_global = 0 + size_buffer = min(ao_num*ao_num*ao_num,16000000) + write(output_BiInts,*) 'Providing the molecular integrals ' + write(output_BiInts,*) 'Buffers : ', 8.*(mo_tot_num_align*(n_j)*(n_k+1) + mo_tot_num_align +& + ao_num+ao_num*ao_num+ size_buffer*3)/(1024*1024), 'MB / core' + + call wall_time(wall_1) + call cpu_time(cpu_1) + mo_integrals_threshold = 0.d0 + + !$OMP PARALLEL PRIVATE(l1,k1,j1,i1,i2,i3,i4,i,j,k,l,c, ii1,kmax, & + !$OMP bielec_tmp_0_idx, bielec_tmp_0, bielec_tmp_1,bielec_tmp_2,bielec_tmp_3,& + !$OMP buffer_i,buffer_value,n_integrals,wall_2,i0,j0,k0,l0) & + !$OMP DEFAULT(NONE) & + !$OMP SHARED(size_buffer,ao_num,mo_tot_num,n_i,n_j,n_k,n_l,mo_tot_num_align,& + !$OMP mo_coef_transp,output_BiInts, & + !$OMP mo_coef_transp_is_built, list_ijkl, & + !$OMP mo_coef_is_built, wall_1, & + !$OMP mo_coef,mo_integrals_threshold,l1_global,ao_integrals_map,mo_integrals_map) + n_integrals = 0 + allocate(bielec_tmp_3(mo_tot_num_align, n_j, n_k), & + bielec_tmp_1(mo_tot_num_align), & + bielec_tmp_0(ao_num,ao_num), & + bielec_tmp_0_idx(ao_num), & + bielec_tmp_2(mo_tot_num_align, n_j), & + buffer_i(size_buffer), & + buffer_value(size_buffer) ) + + !$OMP DO SCHEDULE(guided) + do l1 = 1,ao_num + !DEC$ VECTOR ALIGNED + bielec_tmp_3 = 0.d0 + do k1 = 1,ao_num + !DEC$ VECTOR ALIGNED + bielec_tmp_2 = 0.d0 + do j1 = 1,ao_num + call get_ao_bielec_integrals(j1,k1,l1,ao_num,bielec_tmp_0(1,j1)) + ! call compute_ao_bielec_integrals(j1,k1,l1,ao_num,bielec_tmp_0(1,j1)) + enddo + do j1 = 1,ao_num + kmax = 0 + do i1 = 1,ao_num + c = bielec_tmp_0(i1,j1) + if (c == 0.d0) then + cycle + endif + kmax += 1 + bielec_tmp_0(kmax,j1) = c + bielec_tmp_0_idx(kmax) = i1 + enddo + + if (kmax==0) then + cycle + endif + + !DEC$ VECTOR ALIGNED + bielec_tmp_1 = 0.d0 + ii1=1 + do ii1 = 1,kmax-4,4 + i1 = bielec_tmp_0_idx(ii1) + i2 = bielec_tmp_0_idx(ii1+1) + i3 = bielec_tmp_0_idx(ii1+2) + i4 = bielec_tmp_0_idx(ii1+3) + do i = list_ijkl(1,1), list_ijkl(n_i,1) + bielec_tmp_1(i) = bielec_tmp_1(i) + & + mo_coef_transp(i,i1) * bielec_tmp_0(ii1,j1) + & + mo_coef_transp(i,i2) * bielec_tmp_0(ii1+1,j1) + & + mo_coef_transp(i,i3) * bielec_tmp_0(ii1+2,j1) + & + mo_coef_transp(i,i4) * bielec_tmp_0(ii1+3,j1) + enddo ! i + enddo ! ii1 + + i2 = ii1 + do ii1 = i2,kmax + i1 = bielec_tmp_0_idx(ii1) + do i = list_ijkl(1,1), list_ijkl(n_i,1) + bielec_tmp_1(i) = bielec_tmp_1(i) + mo_coef_transp(i,i1) * bielec_tmp_0(ii1,j1) + enddo ! i + enddo ! ii1 + c = 0.d0 + + do i = list_ijkl(1,1), list_ijkl(n_i,1) + c = max(c,abs(bielec_tmp_1(i))) + if (c>mo_integrals_threshold) exit + enddo + if ( c < mo_integrals_threshold ) then + cycle + endif + + do j0 = 1, n_j + j = list_ijkl(j0,2) + c = mo_coef_transp(j,j1) + if (abs(c) < thr_coef) then + cycle + endif + do i = list_ijkl(1,1), list_ijkl(n_i,1) + bielec_tmp_2(i,j0) = bielec_tmp_2(i,j0) + c * bielec_tmp_1(i) + enddo ! i + enddo ! j + enddo !j1 + if ( maxval(abs(bielec_tmp_2)) < mo_integrals_threshold ) then + cycle + endif + + + do k0 = 1, n_k + k = list_ijkl(k0,3) + c = mo_coef_transp(k,k1) + if (abs(c) < thr_coef) then + cycle + endif + + do j0 = 1, n_j + j = list_ijkl(j0,2) + do i = list_ijkl(1,1), k + bielec_tmp_3(i,j0,k0) = bielec_tmp_3(i,j0,k0) + c* bielec_tmp_2(i,j0) + enddo!i + enddo !j + + enddo !k + enddo !k1 + + + + do l0 = 1,n_l + l = list_ijkl(l0,4) + c = mo_coef_transp(l,l1) + if (abs(c) < thr_coef) then + cycle + endif + j1 = ishft((l*l-l),-1) + do j0 = 1, n_j + j = list_ijkl(j0,2) + if (j > l) then + exit + endif + j1 += 1 + do k0 = 1, n_k + k = list_ijkl(k0,3) + i1 = ishft((k*k-k),-1) + if (i1<=j1) then + continue + else + exit + endif + bielec_tmp_1 = 0.d0 + do i0 = 1, n_i + i = list_ijkl(i0,1) + if (i>k) then + exit + endif + bielec_tmp_1(i) = c*bielec_tmp_3(i,j0,k0) + enddo + + do i = 1, min(k,j1-i1) + if (abs(bielec_tmp_1(i)) < mo_integrals_threshold) then + cycle + endif + n_integrals += 1 + buffer_value(n_integrals) = bielec_tmp_1(i) + !DEC$ FORCEINLINE + call mo_bielec_integrals_index(i,j,k,l,buffer_i(n_integrals)) + if (n_integrals == size_buffer) then + call insert_into_mo_integrals_map(n_integrals,buffer_i,buffer_value,& + real(mo_integrals_threshold,integral_kind)) + n_integrals = 0 + endif + enddo + enddo + enddo + enddo + + !$OMP ATOMIC + l1_global +=1 + call wall_time(wall_2) + write(output_BiInts,*) 100.*float(l1_global)/float(ao_num), '% in ',& + wall_2-wall_1, 's', map_mb(mo_integrals_map) ,'MB' + enddo + !$OMP END DO NOWAIT + deallocate (bielec_tmp_1,bielec_tmp_2,bielec_tmp_3) + + call insert_into_mo_integrals_map(n_integrals,buffer_i,buffer_value,& + real(mo_integrals_threshold,integral_kind)) + deallocate(buffer_i, buffer_value) + !$OMP END PARALLEL + call map_unique(mo_integrals_map) + + call wall_time(wall_2) + call cpu_time(cpu_2) + integer*8 :: get_mo_map_size, mo_map_size + mo_map_size = get_mo_map_size() + + deallocate(list_ijkl) + + + write(output_BiInts,*)'Molecular integrals provided:' + write(output_BiInts,*)' Size of MO map ', map_mb(mo_integrals_map) ,'MB' + write(output_BiInts,*)' Number of MO integrals: ', mo_map_size + write(output_BiInts,*)' cpu time :',cpu_2 - cpu_1, 's' + write(output_BiInts,*)' wall time :',wall_2 - wall_1, 's ( x ', (cpu_2-cpu_1)/(wall_2-wall_1), ')' + + if (write_mo_integrals) then + call dump_mo_integrals(trim(ezfio_filename)//'/work/mo_integrals.bin') + endif + + +end + + +BEGIN_PROVIDER [ double precision, mo_bielec_integral_jj, (mo_tot_num_align,mo_tot_num)] + implicit none + BEGIN_DOC + ! Bi-electronic integrals + END_DOC + integer :: i,j + double precision :: get_mo_bielec_integral + PROVIDE mo_bielec_integrals_in_map + do j= 1, mo_tot_num + !DIR$ VECTOR ALIGNED + do i= 1, mo_tot_num + mo_bielec_integral_jj(i,j) = get_mo_bielec_integral(i,j,i,j,mo_integrals_map) + enddo + ! Padding + do i= mo_tot_num+1, mo_tot_num_align + mo_bielec_integral_jj(i,j) = 0.d0 + enddo + enddo + +END_PROVIDER + +BEGIN_PROVIDER [ double precision, mo_bielec_integral_jj_exchange, (mo_tot_num_align,mo_tot_num)] + implicit none + BEGIN_DOC + ! Bi-electronic integrals - + END_DOC + integer :: i,j + double precision :: get_mo_bielec_integral + PROVIDE mo_bielec_integrals_in_map + + do j = 1, mo_tot_num + !DIR$ VECTOR ALIGNED + do i = 1,mo_tot_num + mo_bielec_integral_jj_exchange(i,j) = get_mo_bielec_integral(i,j,j,i,mo_integrals_map) + enddo + ! Padding + do i= mo_tot_num+1, mo_tot_num_align + mo_bielec_integral_jj_exchange(i,j) = 0.d0 + enddo + enddo + +END_PROVIDER + +BEGIN_PROVIDER [ double precision, mo_bielec_integral_jj_anti, (mo_tot_num_align,mo_tot_num)] + implicit none + BEGIN_DOC + ! Bi-electronic integrals - + END_DOC + integer :: i,j + PROVIDE mo_bielec_integrals_in_map + + do j = 1, mo_tot_num + !DIR$ VECTOR ALIGNED + do i = 1,mo_tot_num_align + mo_bielec_integral_jj_anti(i,j) = mo_bielec_integral_jj(i,j) - mo_bielec_integral_jj_exchange(i,j) + enddo + enddo + +END_PROVIDER + diff --git a/src/BiInts/options.irp.f b/src/BiInts/options.irp.f new file mode 100644 index 00000000..c821661c --- /dev/null +++ b/src/BiInts/options.irp.f @@ -0,0 +1,109 @@ +BEGIN_PROVIDER [ logical, write_mo_integrals ] + implicit none + BEGIN_DOC + ! If true, write MO integrals in EZFIO + END_DOC + + logical :: has + PROVIDE ezfio_filename + call ezfio_has_bielec_integrals_write_mo_integrals(has) + if (has) then + call ezfio_get_bielec_integrals_write_mo_integrals(write_mo_integrals) + else + write_mo_integrals = .False. + call ezfio_set_bielec_integrals_write_mo_integrals(write_mo_integrals) + endif + +END_PROVIDER + +BEGIN_PROVIDER [ logical, write_ao_integrals ] + implicit none + BEGIN_DOC + ! If true, write AO integrals in EZFIO + END_DOC + + logical :: has + PROVIDE ezfio_filename + call ezfio_has_bielec_integrals_write_ao_integrals(has) + if (has) then + call ezfio_get_bielec_integrals_write_ao_integrals(write_ao_integrals) + else + write_ao_integrals = .False. + call ezfio_set_bielec_integrals_write_ao_integrals(write_ao_integrals) + endif + +END_PROVIDER + +BEGIN_PROVIDER [ logical, read_mo_integrals ] + implicit none + BEGIN_DOC + ! If true, read MO integrals in EZFIO + END_DOC + + logical :: has + PROVIDE ezfio_filename + call ezfio_has_bielec_integrals_read_mo_integrals(has) + if (has) then + call ezfio_get_bielec_integrals_read_mo_integrals(read_mo_integrals) + else + read_mo_integrals = .False. + call ezfio_set_bielec_integrals_read_mo_integrals(read_mo_integrals) + endif + +END_PROVIDER + +BEGIN_PROVIDER [ logical, read_ao_integrals ] + implicit none + BEGIN_DOC + ! If true, read AO integrals in EZFIO + END_DOC + + logical :: has + PROVIDE ezfio_filename + call ezfio_has_bielec_integrals_read_ao_integrals(has) + if (has) then + call ezfio_get_bielec_integrals_read_ao_integrals(read_ao_integrals) + else + read_ao_integrals = .False. + call ezfio_set_bielec_integrals_read_ao_integrals(read_ao_integrals) + endif + +END_PROVIDER + +BEGIN_PROVIDER [ double precision, ao_integrals_threshold ] + implicit none + BEGIN_DOC + ! If < ao_integrals_threshold, = 0 + END_DOC + + logical :: has + PROVIDE ezfio_filename + call ezfio_has_bielec_integrals_threshold_ao(has) + if (has) then + call ezfio_get_bielec_integrals_threshold_ao(ao_integrals_threshold) + else + ao_integrals_threshold = 1.d-15 + call ezfio_set_bielec_integrals_threshold_ao(ao_integrals_threshold) + endif + +END_PROVIDER + + +BEGIN_PROVIDER [ double precision, mo_integrals_threshold ] + implicit none + BEGIN_DOC + ! If < mo_integrals_threshold, = 0 + END_DOC + + logical :: has + PROVIDE ezfio_filename + call ezfio_has_bielec_integrals_threshold_mo(has) + if (has) then + call ezfio_get_bielec_integrals_threshold_mo(mo_integrals_threshold) + else + mo_integrals_threshold = 1.d-11 + call ezfio_set_bielec_integrals_threshold_mo(mo_integrals_threshold) + endif + +END_PROVIDER + diff --git a/src/BiInts/tests/Makefile b/src/BiInts/tests/Makefile new file mode 100644 index 00000000..77bd84ba --- /dev/null +++ b/src/BiInts/tests/Makefile @@ -0,0 +1,33 @@ +OPENMP =1 +PROFILE =0 +DEBUG = 0 + +IRPF90+= -I tests + +REF_FILES=$(subst %.irp.f, %.ref, $(wildcard *.irp.f)) + +.PHONY: clean executables serial_tests parallel_tests + +all: clean executables serial_tests parallel_tests + +parallel_tests: $(REF_FILES) + @echo ; echo " ---- Running parallel tests ----" ; echo + @OMP_NUM_THREADS=10 ${QPACKAGE_ROOT}/scripts/run_tests.py + +serial_tests: $(REF_FILES) + @echo ; echo " ---- Running serial tests ----" ; echo + @OMP_NUM_THREADS=1 ${QPACKAGE_ROOT}/scripts/run_tests.py + +executables: $(wildcard *.irp.f) veryclean + $(MAKE) -C .. + +%.ref: $(wildcard $(QPACKAGE_ROOT)/data/inputs/*.md5) executables + $(QPACKAGE_ROOT)/scripts/create_test_ref.sh $* + +clean: + $(MAKE) -C .. clean + +veryclean: + $(MAKE) -C .. veryclean + + diff --git a/src/BiInts/tests/test_ao.ref b/src/BiInts/tests/test_ao.ref new file mode 100644 index 00000000..d34ef769 --- /dev/null +++ b/src/BiInts/tests/test_ao.ref @@ -0,0 +1,424 @@ +data = { + 'AlCl.ezfio' : { + }, + 'Al.ezfio' : { + }, + 'Al+.ezfio' : { + }, + 'AlH2.ezfio' : { + }, + 'AlH3.ezfio' : { + }, + 'AlH.ezfio' : { + }, + 'BCl.ezfio' : { + }, + 'BeCl.ezfio' : { + }, + 'Be.ezfio' : { + }, + 'Be+.ezfio' : { + }, + 'BeF.ezfio' : { + }, + 'BeH2.ezfio' : { + }, + 'BeH.ezfio' : { + }, + 'BeO.ezfio' : { + }, + 'BeOH.ezfio' : { + }, + 'BeS.ezfio' : { + }, + 'B.ezfio' : { + }, + 'B+.ezfio' : { + }, + 'BH2.ezfio' : { + }, + 'BH3.ezfio' : { + }, + 'BH.ezfio' : { + }, + 'BO.ezfio' : { + }, + 'BS.ezfio' : { + }, + 'C2.ezfio' : { + }, + 'C2H2.ezfio' : { + }, + 'C2H2+.ezfio' : { + }, + 'C2H3.ezfio' : { + }, + 'C2H3+.ezfio' : { + }, + 'C2H4.ezfio' : { + }, + 'C2H4+.ezfio' : { + }, + 'C2H5.ezfio' : { + }, + 'C2H6.ezfio' : { + }, + 'C2H.ezfio' : { + }, + 'CCl.ezfio' : { + }, + 'C-.ezfio' : { + }, + 'C.ezfio' : { + }, + 'C+.ezfio' : { + }, + 'CF.ezfio' : { + }, + 'CH2_1A1.ezfio' : { + }, + 'CH2_3B1.ezfio' : { + }, + 'CH2-.ezfio' : { + }, + 'CH3-.ezfio' : { + }, + 'CH3.ezfio' : { + }, + 'CH4.ezfio' : { + }, + 'CH4+.ezfio' : { + }, + 'CH-.ezfio' : { + }, + 'CH.ezfio' : { + }, + 'Cl2-.ezfio' : { + }, + 'Cl2.ezfio' : { + }, + 'Cl2+.ezfio' : { + }, + 'Cl-.ezfio' : { + }, + 'Cl.ezfio' : { + }, + 'Cl+.ezfio' : { + }, + 'ClH2+.ezfio' : { + }, + 'ClH.ezfio' : { + }, + 'ClH+.ezfio' : { + }, + 'ClS.ezfio' : { + }, + 'ClSiH3.ezfio' : { + }, + 'CN-.ezfio' : { + }, + 'CN.ezfio' : { + }, + 'CO2.ezfio' : { + }, + 'CO.ezfio' : { + }, + 'CO+.ezfio' : { + }, + 'COS.ezfio' : { + }, + 'CP.ezfio' : { + }, + 'CS2.ezfio' : { + }, + 'CS.ezfio' : { + }, + 'CS+.ezfio' : { + }, + 'CSi.ezfio' : { + }, + 'F2.ezfio' : { + }, + 'FAl.ezfio' : { + }, + 'FCl.ezfio' : { + }, + 'FCl+.ezfio' : { + }, + 'F-.ezfio' : { + }, + 'F.ezfio' : { + }, + 'F+.ezfio' : { + }, + 'FH.ezfio' : { + }, + 'FH+.ezfio' : { + }, + 'FMg.ezfio' : { + }, + 'FNa.ezfio' : { + }, + 'FP.ezfio' : { + }, + 'FS.ezfio' : { + }, + 'FSi.ezfio' : { + }, + 'FSiH3.ezfio' : { + }, + 'H2CNH.ezfio' : { + }, + 'H2CO.ezfio' : { + }, + 'H2CPH.ezfio' : { + }, + 'H2CS.ezfio' : { + }, + 'H2.ezfio' : { + }, + 'H2NNH2.ezfio' : { + }, + 'H2PPH2.ezfio' : { + }, + 'H3CCl.ezfio' : { + }, + 'H3CF.ezfio' : { + }, + 'H3CNH2.ezfio' : { + }, + 'H3COH.ezfio' : { + }, + 'H3CSH.ezfio' : { + }, + 'H3SiSiH3.ezfio' : { + }, + 'HBO.ezfio' : { + }, + 'HBS.ezfio' : { + }, + 'HCF.ezfio' : { + }, + 'HCN.ezfio' : { + }, + 'HCO.ezfio' : { + }, + 'HCP.ezfio' : { + }, + 'H.ezfio' : { + }, + 'HNO.ezfio' : { + }, + 'HOCl.ezfio' : { + }, + 'HOF.ezfio' : { + }, + 'HOMg.ezfio' : { + }, + 'HONa.ezfio' : { + }, + 'HOO.ezfio' : { + }, + 'HOOH.ezfio' : { + }, + 'HSSH.ezfio' : { + }, + 'Li2.ezfio' : { + }, + 'LiCl.ezfio' : { + }, + 'Li.ezfio' : { + }, + 'Li+.ezfio' : { + }, + 'LiF.ezfio' : { + }, + 'LiH.ezfio' : { + }, + 'LiN.ezfio' : { + }, + 'LiO.ezfio' : { + }, + 'LiOH.ezfio' : { + }, + 'MgCl.ezfio' : { + }, + 'Mg.ezfio' : { + }, + 'Mg+.ezfio' : { + }, + 'MgH.ezfio' : { + }, + 'MgS.ezfio' : { + }, + 'N2.ezfio' : { + }, + 'N2+.ezfio' : { + }, + 'Na2.ezfio' : { + }, + 'NaCl.ezfio' : { + }, + 'Na.ezfio' : { + }, + 'Na+.ezfio' : { + }, + 'NaH.ezfio' : { + }, + 'N.ezfio' : { + }, + 'N+.ezfio' : { + }, + 'NF.ezfio' : { + }, + 'NH2-.ezfio' : { + }, + 'NH2.ezfio' : { + }, + 'NH3.ezfio' : { + }, + 'NH3+.ezfio' : { + }, + 'NH4+.ezfio' : { + }, + 'NH-.ezfio' : { + }, + 'NH.ezfio' : { + }, + 'NO-.ezfio' : { + }, + 'NO.ezfio' : { + }, + 'NP.ezfio' : { + }, + 'NS.ezfio' : { + }, + 'NSi.ezfio' : { + }, + 'O2Cl.ezfio' : { + }, + 'O2-.ezfio' : { + }, + 'O2.ezfio' : { + }, + 'O2+.ezfio' : { + }, + 'O2S.ezfio' : { + }, + 'O2Si.ezfio' : { + }, + 'O3.ezfio' : { + }, + 'OCl.ezfio' : { + }, + 'O-.ezfio' : { + }, + 'O.ezfio' : { + }, + 'O+.ezfio' : { + }, + 'OH2.ezfio' : { + }, + 'OH2+.ezfio' : { + }, + 'OH3+.ezfio' : { + }, + 'OH-.ezfio' : { + }, + 'OH.ezfio' : { + }, + 'OH+.ezfio' : { + }, + 'OMg.ezfio' : { + }, + 'ONa.ezfio' : { + }, + 'OP-.ezfio' : { + }, + 'OP.ezfio' : { + }, + 'OPH.ezfio' : { + }, + 'OS.ezfio' : { + }, + 'OSi.ezfio' : { + }, + 'P2.ezfio' : { + }, + 'P2+.ezfio' : { + }, + 'PCl.ezfio' : { + }, + 'P-.ezfio' : { + }, + 'P.ezfio' : { + }, + 'PH2-.ezfio' : { + }, + 'PH2.ezfio' : { + }, + 'PH2+.ezfio' : { + }, + 'PH3.ezfio' : { + }, + 'PH3+.ezfio' : { + }, + 'PH4+.ezfio' : { + }, + 'PH-.ezfio' : { + }, + 'PH.ezfio' : { + }, + 'PS.ezfio' : { + }, + 'S2-.ezfio' : { + }, + 'S2.ezfio' : { + }, + 'S-.ezfio' : { + }, + 'S.ezfio' : { + }, + 'S+.ezfio' : { + }, + 'SH2.ezfio' : { + }, + 'SH2+.ezfio' : { + }, + 'SH3+.ezfio' : { + }, + 'SH-.ezfio' : { + }, + 'SH.ezfio' : { + }, + 'SH+.ezfio' : { + }, + 'Si2.ezfio' : { + }, + 'SiCl.ezfio' : { + }, + 'Si-.ezfio' : { + }, + 'Si.ezfio' : { + }, + 'SiH2_1A1.ezfio' : { + }, + 'SiH2_3B1.ezfio' : { + }, + 'SiH2-.ezfio' : { + }, + 'SiH3-.ezfio' : { + }, + 'SiH3.ezfio' : { + }, + 'SiH4.ezfio' : { + }, + 'SiH4+.ezfio' : { + }, + 'SiH-.ezfio' : { + }, + 'SiH.ezfio' : { + }, + 'SiS.ezfio' : { + }, +} diff --git a/src/BiInts/tests/test_ao_map.irp.f b/src/BiInts/tests/test_ao_map.irp.f new file mode 100644 index 00000000..2d41df81 --- /dev/null +++ b/src/BiInts/tests/test_ao_map.irp.f @@ -0,0 +1,35 @@ +program test_ao_map + implicit none + integer :: i,j,k,l + double precision :: get_ao_bielec_integral, tmp + double precision, allocatable :: buffer_value(:,:) + double precision :: ao_bielec_integral,tmp2 + integer :: non_zero_int + double precision :: s + integer :: OK + + allocate(buffer_value(ao_num,2)) + + OK = 1 + do l=1,ao_num,7 + s = 0.d0 + do k=1, ao_num + do j=1,l + call compute_ao_bielec_integrals(j,k,l,ao_num,buffer_value(1,1)) + call get_ao_bielec_integrals(j,k,l,ao_num,buffer_value(1,2)) + do i=1,k + s += (buffer_value(i,1)-buffer_value(i,2))*(buffer_value(i,1)-buffer_value(i,2)) + enddo + enddo + enddo + if (s > 1.d-20) then + print *, l, ' : ', s + OK=0 + endif + enddo + + deallocate(buffer_value) + print *, ' OK : ', OK + + +end diff --git a/src/BiInts/tests/test_ao_map.ref b/src/BiInts/tests/test_ao_map.ref new file mode 100644 index 00000000..ca3e78dc --- /dev/null +++ b/src/BiInts/tests/test_ao_map.ref @@ -0,0 +1,635 @@ +data = { + 'AlCl.ezfio' : { + 'OK' : 1, + }, + 'Al.ezfio' : { + 'OK' : 1, + }, + 'Al+.ezfio' : { + 'OK' : 1, + }, + 'AlH2.ezfio' : { + 'OK' : 1, + }, + 'AlH3.ezfio' : { + 'OK' : 1, + }, + 'AlH.ezfio' : { + 'OK' : 1, + }, + 'BCl.ezfio' : { + 'OK' : 1, + }, + 'BeCl.ezfio' : { + 'OK' : 1, + }, + 'Be.ezfio' : { + 'OK' : 1, + }, + 'Be+.ezfio' : { + 'OK' : 1, + }, + 'BeF.ezfio' : { + 'OK' : 1, + }, + 'BeH2.ezfio' : { + 'OK' : 1, + }, + 'BeH.ezfio' : { + 'OK' : 1, + }, + 'BeO.ezfio' : { + 'OK' : 1, + }, + 'BeOH.ezfio' : { + 'OK' : 1, + }, + 'BeS.ezfio' : { + 'OK' : 1, + }, + 'B.ezfio' : { + 'OK' : 1, + }, + 'B+.ezfio' : { + 'OK' : 1, + }, + 'BH2.ezfio' : { + 'OK' : 1, + }, + 'BH3.ezfio' : { + 'OK' : 1, + }, + 'BH.ezfio' : { + 'OK' : 1, + }, + 'BO.ezfio' : { + 'OK' : 1, + }, + 'BS.ezfio' : { + 'OK' : 1, + }, + 'C2.ezfio' : { + 'OK' : 1, + }, + 'C2H2.ezfio' : { + 'OK' : 1, + }, + 'C2H2+.ezfio' : { + 'OK' : 1, + }, + 'C2H3.ezfio' : { + 'OK' : 1, + }, + 'C2H3+.ezfio' : { + 'OK' : 1, + }, + 'C2H4.ezfio' : { + 'OK' : 1, + }, + 'C2H4+.ezfio' : { + 'OK' : 1, + }, + 'C2H5.ezfio' : { + 'OK' : 1, + }, + 'C2H6.ezfio' : { + 'OK' : 1, + }, + 'C2H.ezfio' : { + 'OK' : 1, + }, + 'CCl.ezfio' : { + 'OK' : 1, + }, + 'C-.ezfio' : { + 'OK' : 1, + }, + 'C.ezfio' : { + 'OK' : 1, + }, + 'C+.ezfio' : { + 'OK' : 1, + }, + 'CF.ezfio' : { + 'OK' : 1, + }, + 'CH2_1A1.ezfio' : { + 'OK' : 1, + }, + 'CH2_3B1.ezfio' : { + 'OK' : 1, + }, + 'CH2-.ezfio' : { + 'OK' : 1, + }, + 'CH3-.ezfio' : { + 'OK' : 1, + }, + 'CH3.ezfio' : { + 'OK' : 1, + }, + 'CH4.ezfio' : { + 'OK' : 1, + }, + 'CH4+.ezfio' : { + 'OK' : 1, + }, + 'CH-.ezfio' : { + 'OK' : 1, + }, + 'CH.ezfio' : { + 'OK' : 1, + }, + 'Cl2-.ezfio' : { + 'OK' : 1, + }, + 'Cl2.ezfio' : { + 'OK' : 1, + }, + 'Cl2+.ezfio' : { + 'OK' : 1, + }, + 'Cl-.ezfio' : { + 'OK' : 1, + }, + 'Cl.ezfio' : { + 'OK' : 1, + }, + 'Cl+.ezfio' : { + 'OK' : 1, + }, + 'ClH2+.ezfio' : { + 'OK' : 1, + }, + 'ClH.ezfio' : { + 'OK' : 1, + }, + 'ClH+.ezfio' : { + 'OK' : 1, + }, + 'ClS.ezfio' : { + 'OK' : 1, + }, + 'ClSiH3.ezfio' : { + 'OK' : 1, + }, + 'CN-.ezfio' : { + 'OK' : 1, + }, + 'CN.ezfio' : { + 'OK' : 1, + }, + 'CO2.ezfio' : { + 'OK' : 1, + }, + 'CO.ezfio' : { + 'OK' : 1, + }, + 'CO+.ezfio' : { + 'OK' : 1, + }, + 'COS.ezfio' : { + 'OK' : 1, + }, + 'CP.ezfio' : { + 'OK' : 1, + }, + 'CS2.ezfio' : { + 'OK' : 1, + }, + 'CS.ezfio' : { + 'OK' : 1, + }, + 'CS+.ezfio' : { + 'OK' : 1, + }, + 'CSi.ezfio' : { + 'OK' : 1, + }, + 'F2.ezfio' : { + 'OK' : 1, + }, + 'FAl.ezfio' : { + 'OK' : 1, + }, + 'FCl.ezfio' : { + 'OK' : 1, + }, + 'FCl+.ezfio' : { + 'OK' : 1, + }, + 'F-.ezfio' : { + 'OK' : 1, + }, + 'F.ezfio' : { + 'OK' : 1, + }, + 'F+.ezfio' : { + 'OK' : 1, + }, + 'FH.ezfio' : { + 'OK' : 1, + }, + 'FH+.ezfio' : { + 'OK' : 1, + }, + 'FMg.ezfio' : { + 'OK' : 1, + }, + 'FNa.ezfio' : { + 'OK' : 1, + }, + 'FP.ezfio' : { + 'OK' : 1, + }, + 'FS.ezfio' : { + 'OK' : 1, + }, + 'FSi.ezfio' : { + 'OK' : 1, + }, + 'FSiH3.ezfio' : { + 'OK' : 1, + }, + 'H2CNH.ezfio' : { + 'OK' : 1, + }, + 'H2CO.ezfio' : { + 'OK' : 1, + }, + 'H2CPH.ezfio' : { + 'OK' : 1, + }, + 'H2CS.ezfio' : { + 'OK' : 1, + }, + 'H2.ezfio' : { + 'OK' : 1, + }, + 'H2NNH2.ezfio' : { + 'OK' : 1, + }, + 'H2PPH2.ezfio' : { + 'OK' : 1, + }, + 'H3CCl.ezfio' : { + 'OK' : 1, + }, + 'H3CF.ezfio' : { + 'OK' : 1, + }, + 'H3CNH2.ezfio' : { + 'OK' : 1, + }, + 'H3COH.ezfio' : { + 'OK' : 1, + }, + 'H3CSH.ezfio' : { + 'OK' : 1, + }, + 'H3SiSiH3.ezfio' : { + 'OK' : 1, + }, + 'HBO.ezfio' : { + 'OK' : 1, + }, + 'HBS.ezfio' : { + 'OK' : 1, + }, + 'HCF.ezfio' : { + 'OK' : 1, + }, + 'HCN.ezfio' : { + 'OK' : 1, + }, + 'HCO.ezfio' : { + 'OK' : 1, + }, + 'HCP.ezfio' : { + 'OK' : 1, + }, + 'H.ezfio' : { + 'OK' : 1, + }, + 'HNO.ezfio' : { + 'OK' : 1, + }, + 'HOCl.ezfio' : { + 'OK' : 1, + }, + 'HOF.ezfio' : { + 'OK' : 1, + }, + 'HOMg.ezfio' : { + 'OK' : 1, + }, + 'HONa.ezfio' : { + 'OK' : 1, + }, + 'HOO.ezfio' : { + 'OK' : 1, + }, + 'HOOH.ezfio' : { + 'OK' : 1, + }, + 'HSSH.ezfio' : { + 'OK' : 1, + }, + 'Li2.ezfio' : { + 'OK' : 1, + }, + 'LiCl.ezfio' : { + 'OK' : 1, + }, + 'Li.ezfio' : { + 'OK' : 1, + }, + 'Li+.ezfio' : { + 'OK' : 1, + }, + 'LiF.ezfio' : { + 'OK' : 1, + }, + 'LiH.ezfio' : { + 'OK' : 1, + }, + 'LiN.ezfio' : { + 'OK' : 1, + }, + 'LiO.ezfio' : { + 'OK' : 1, + }, + 'LiOH.ezfio' : { + 'OK' : 1, + }, + 'MgCl.ezfio' : { + 'OK' : 1, + }, + 'Mg.ezfio' : { + 'OK' : 1, + }, + 'Mg+.ezfio' : { + 'OK' : 1, + }, + 'MgH.ezfio' : { + 'OK' : 1, + }, + 'MgS.ezfio' : { + 'OK' : 1, + }, + 'N2.ezfio' : { + 'OK' : 1, + }, + 'N2+.ezfio' : { + 'OK' : 1, + }, + 'Na2.ezfio' : { + 'OK' : 1, + }, + 'NaCl.ezfio' : { + 'OK' : 1, + }, + 'Na.ezfio' : { + 'OK' : 1, + }, + 'Na+.ezfio' : { + 'OK' : 1, + }, + 'NaH.ezfio' : { + 'OK' : 1, + }, + 'N.ezfio' : { + 'OK' : 1, + }, + 'N+.ezfio' : { + 'OK' : 1, + }, + 'NF.ezfio' : { + 'OK' : 1, + }, + 'NH2-.ezfio' : { + 'OK' : 1, + }, + 'NH2.ezfio' : { + 'OK' : 1, + }, + 'NH3.ezfio' : { + 'OK' : 1, + }, + 'NH3+.ezfio' : { + 'OK' : 1, + }, + 'NH4+.ezfio' : { + 'OK' : 1, + }, + 'NH-.ezfio' : { + 'OK' : 1, + }, + 'NH.ezfio' : { + 'OK' : 1, + }, + 'NO-.ezfio' : { + 'OK' : 1, + }, + 'NO.ezfio' : { + 'OK' : 1, + }, + 'NP.ezfio' : { + 'OK' : 1, + }, + 'NS.ezfio' : { + 'OK' : 1, + }, + 'NSi.ezfio' : { + 'OK' : 1, + }, + 'O2Cl.ezfio' : { + 'OK' : 1, + }, + 'O2-.ezfio' : { + 'OK' : 1, + }, + 'O2.ezfio' : { + 'OK' : 1, + }, + 'O2+.ezfio' : { + 'OK' : 1, + }, + 'O2S.ezfio' : { + 'OK' : 1, + }, + 'O2Si.ezfio' : { + 'OK' : 1, + }, + 'O3.ezfio' : { + 'OK' : 1, + }, + 'OCl.ezfio' : { + 'OK' : 1, + }, + 'O-.ezfio' : { + 'OK' : 1, + }, + 'O.ezfio' : { + 'OK' : 1, + }, + 'O+.ezfio' : { + 'OK' : 1, + }, + 'OH2.ezfio' : { + 'OK' : 1, + }, + 'OH2+.ezfio' : { + 'OK' : 1, + }, + 'OH3+.ezfio' : { + 'OK' : 1, + }, + 'OH-.ezfio' : { + 'OK' : 1, + }, + 'OH.ezfio' : { + 'OK' : 1, + }, + 'OH+.ezfio' : { + 'OK' : 1, + }, + 'OMg.ezfio' : { + 'OK' : 1, + }, + 'ONa.ezfio' : { + 'OK' : 1, + }, + 'OP-.ezfio' : { + 'OK' : 1, + }, + 'OP.ezfio' : { + 'OK' : 1, + }, + 'OPH.ezfio' : { + 'OK' : 1, + }, + 'OS.ezfio' : { + 'OK' : 1, + }, + 'OSi.ezfio' : { + 'OK' : 1, + }, + 'P2.ezfio' : { + 'OK' : 1, + }, + 'P2+.ezfio' : { + 'OK' : 1, + }, + 'PCl.ezfio' : { + 'OK' : 1, + }, + 'P-.ezfio' : { + 'OK' : 1, + }, + 'P.ezfio' : { + 'OK' : 1, + }, + 'PH2-.ezfio' : { + 'OK' : 1, + }, + 'PH2.ezfio' : { + 'OK' : 1, + }, + 'PH2+.ezfio' : { + 'OK' : 1, + }, + 'PH3.ezfio' : { + 'OK' : 1, + }, + 'PH3+.ezfio' : { + 'OK' : 1, + }, + 'PH4+.ezfio' : { + 'OK' : 1, + }, + 'PH-.ezfio' : { + 'OK' : 1, + }, + 'PH.ezfio' : { + 'OK' : 1, + }, + 'PS.ezfio' : { + 'OK' : 1, + }, + 'S2-.ezfio' : { + 'OK' : 1, + }, + 'S2.ezfio' : { + 'OK' : 1, + }, + 'S-.ezfio' : { + 'OK' : 1, + }, + 'S.ezfio' : { + 'OK' : 1, + }, + 'S+.ezfio' : { + 'OK' : 1, + }, + 'SH2.ezfio' : { + 'OK' : 1, + }, + 'SH2+.ezfio' : { + 'OK' : 1, + }, + 'SH3+.ezfio' : { + 'OK' : 1, + }, + 'SH-.ezfio' : { + 'OK' : 1, + }, + 'SH.ezfio' : { + 'OK' : 1, + }, + 'SH+.ezfio' : { + 'OK' : 1, + }, + 'Si2.ezfio' : { + 'OK' : 1, + }, + 'SiCl.ezfio' : { + 'OK' : 1, + }, + 'Si-.ezfio' : { + 'OK' : 1, + }, + 'Si.ezfio' : { + 'OK' : 1, + }, + 'SiH2_1A1.ezfio' : { + 'OK' : 1, + }, + 'SiH2_3B1.ezfio' : { + 'OK' : 1, + }, + 'SiH2-.ezfio' : { + 'OK' : 1, + }, + 'SiH3-.ezfio' : { + 'OK' : 1, + }, + 'SiH3.ezfio' : { + 'OK' : 1, + }, + 'SiH4.ezfio' : { + 'OK' : 1, + }, + 'SiH4+.ezfio' : { + 'OK' : 1, + }, + 'SiH-.ezfio' : { + 'OK' : 1, + }, + 'SiH.ezfio' : { + 'OK' : 1, + }, + 'SiS.ezfio' : { + 'OK' : 1, + }, +} diff --git a/src/BiInts/tests/test_mo.irp.f b/src/BiInts/tests/test_mo.irp.f new file mode 100644 index 00000000..014652bf --- /dev/null +++ b/src/BiInts/tests/test_mo.irp.f @@ -0,0 +1,15 @@ +program test + implicit none + integer :: i,j,k,l + double precision :: accu,accu2,get_mo_bielec_integral,accu3 + accu = 0.d0 + accu3 = 0.d0 + do i = 1, elec_beta_num + do j = 1, elec_beta_num + accu += (2.d0 * get_mo_bielec_integral(i,j,i,j,mo_integrals_map) - get_mo_bielec_integral(i,j,j,i,mo_integrals_map) ) + accu3 += (2.d0 * mo_bielec_integral_jj(i,j) - mo_bielec_integral_jj_exchange(i,j) ) + enddo + enddo + print *,'closed EE 1 : ',accu + print *,'closed EE 2 : ',accu3 +end diff --git a/src/BiInts/tests/test_mo.ref b/src/BiInts/tests/test_mo.ref new file mode 100644 index 00000000..0cee6ebe --- /dev/null +++ b/src/BiInts/tests/test_mo.ref @@ -0,0 +1,846 @@ +data = { + 'AlCl.ezfio' : { + 'closed EE 1' : 324.604108318567, +'closed EE 2' : 324.604108318567, + }, + 'Al.ezfio' : { + 'closed EE 1' : 90.6007174995339, +'closed EE 2' : 90.6007174995339, + }, + 'Al+.ezfio' : { + 'closed EE 1' : 91.1277124469289, +'closed EE 2' : 91.1277124469289, + }, + 'AlH2.ezfio' : { + 'closed EE 1' : 98.3932867799503, +'closed EE 2' : 98.3932867799503, + }, + 'AlH3.ezfio' : { + 'closed EE 1' : 108.322607071765, +'closed EE 2' : 108.322607071765, + }, + 'AlH.ezfio' : { + 'closed EE 1' : 98.9361676259577, +'closed EE 2' : 98.9361676259577, + }, + 'BCl.ezfio' : { + 'closed EE 1' : 208.669206621826, +'closed EE 2' : 208.669206621826, + }, + 'BeCl.ezfio' : { + 'closed EE 1' : 195.478154262549, +'closed EE 2' : 195.478154262549, + }, + 'Be.ezfio' : { + 'closed EE 1' : 4.48960803917240, +'closed EE 2' : 4.48960803917240, + }, + 'Be+.ezfio' : { + 'closed EE 1' : 2.27296598314244, +'closed EE 2' : 2.27296598314244, + }, + 'BeF.ezfio' : { + 'closed EE 1' : 55.2667378076567, +'closed EE 2' : 55.2667378076567, + }, + 'BeH2.ezfio' : { + 'closed EE 1' : 7.91863331588170, +'closed EE 2' : 7.91863331588170, + }, + 'BeH.ezfio' : { + 'closed EE 1' : 4.60494282657902, +'closed EE 2' : 4.60494282657902, + }, + 'BeO.ezfio' : { + 'closed EE 1' : 47.0178059698967, +'closed EE 2' : 47.0178059698967, + }, + 'BeOH.ezfio' : { + 'closed EE 1' : 47.4095700462092, +'closed EE 2' : 47.4095700462092, + }, + 'BeS.ezfio' : { + 'closed EE 1' : 177.044582851287, +'closed EE 2' : 177.044582851287, + }, + 'B.ezfio' : { + 'closed EE 1' : 5.87051480040891, +'closed EE 2' : 5.87051480040891, + }, + 'B+.ezfio' : { + 'closed EE 1' : 6.08872058329241, +'closed EE 2' : 6.08872058329241, + }, + 'BH2.ezfio' : { + 'closed EE 1' : 9.82092714954701, +'closed EE 2' : 9.82092714954701, + }, + 'BH3.ezfio' : { + 'closed EE 1' : 15.1738039588984, +'closed EE 2' : 15.1738039588984, + }, + 'BH.ezfio' : { + 'closed EE 1' : 9.89523045878027, +'closed EE 2' : 9.89523045878027, + }, + 'BO.ezfio' : { + 'closed EE 1' : 49.5244841428460, +'closed EE 2' : 49.5244841428460, + }, + 'BS.ezfio' : { + 'closed EE 1' : 179.847811135927, +'closed EE 2' : 179.847811135927, + }, + 'C2.ezfio' : { + 'closed EE 1' : 40.1266484535661, +'closed EE 2' : 40.1266484535661, + }, + 'C2H2.ezfio' : { + 'closed EE 1' : 50.0072453359780, +'closed EE 2' : 50.0072453359780, + }, + 'C2H2+.ezfio' : { + 'closed EE 1' : 38.8813663363468, +'closed EE 2' : 38.8813663363468, + }, + 'C2H3.ezfio' : { + 'closed EE 1' : 47.7960364243417, +'closed EE 2' : 47.7960364243417, + }, + 'C2H3+.ezfio' : { + 'closed EE 1' : 49.4483058762631, +'closed EE 2' : 49.4483058762631, + }, + 'C2H4.ezfio' : { + 'closed EE 1' : 58.5930731008834, +'closed EE 2' : 58.5930731008834, + }, + 'C2H4+.ezfio' : { + 'closed EE 1' : 46.8050385693024, +'closed EE 2' : 46.8050385693024, + }, + 'C2H5.ezfio' : { + 'closed EE 1' : 55.8648362479812, +'closed EE 2' : 55.8648362479812, + }, + 'C2H6.ezfio' : { + 'closed EE 1' : 67.4231505770701, +'closed EE 2' : 67.4231505770701, + }, + 'C2H.ezfio' : { + 'closed EE 1' : 39.7941997946625, +'closed EE 2' : 39.7941997946625, + }, + 'CCl.ezfio' : { + 'closed EE 1' : 211.332787485635, +'closed EE 2' : 211.332787485635, + }, + 'C-.ezfio' : { + 'closed EE 1' : 7.05668474724391, +'closed EE 2' : 7.05668474724391, + }, + 'C.ezfio' : { + 'closed EE 1' : 7.22275234701424, +'closed EE 2' : 7.22275234701424, + }, + 'C+.ezfio' : { + 'closed EE 1' : 7.43316493990123, +'closed EE 2' : 7.43316493990123, + }, + 'CF.ezfio' : { + 'closed EE 1' : 67.5749607457798, +'closed EE 2' : 67.5749607457798, + }, + 'CH2_1A1.ezfio' : { + 'closed EE 1' : 18.6032101860148, +'closed EE 2' : 18.6032101860148, + }, + 'CH2_3B1.ezfio' : { + 'closed EE 1' : 11.7811748791077, +'closed EE 2' : 11.7811748791077, + }, + 'CH2-.ezfio' : { + 'closed EE 1' : 17.4224231187863, +'closed EE 2' : 17.4224231187863, + }, + 'CH3-.ezfio' : { + 'closed EE 1' : 25.0349104359763, +'closed EE 2' : 25.0349104359763, + }, + 'CH3.ezfio' : { + 'closed EE 1' : 18.1595297280646, +'closed EE 2' : 18.1595297280646, + }, + 'CH4.ezfio' : { + 'closed EE 1' : 26.0266830207216, +'closed EE 2' : 26.0266830207216, + }, + 'CH4+.ezfio' : { + 'closed EE 1' : 18.4285252097398, +'closed EE 2' : 18.4285252097398, + }, + 'CH-.ezfio' : { + 'closed EE 1' : 11.4343449775700, +'closed EE 2' : 11.4343449775700, + }, + 'CH.ezfio' : { + 'closed EE 1' : 12.0838764958443, +'closed EE 2' : 12.0838764958443, + }, + 'Cl2-.ezfio' : { + 'closed EE 1' : 402.671444472536, +'closed EE 2' : 402.671444472536, + }, + 'Cl2.ezfio' : { + 'closed EE 1' : 425.050531461454, +'closed EE 2' : 425.050531461454, + }, + 'Cl2+.ezfio' : { + 'closed EE 1' : 403.762363169362, +'closed EE 2' : 403.762363169362, + }, + 'Cl-.ezfio' : { + 'closed EE 1' : 182.327075272228, +'closed EE 2' : 182.327075272228, + }, + 'Cl.ezfio' : { + 'closed EE 1' : 165.731725346204, +'closed EE 2' : 165.731725346204, + }, + 'Cl+.ezfio' : { + 'closed EE 1' : 149.193201694481, +'closed EE 2' : 149.193201694481, + }, + 'ClH2+.ezfio' : { + 'closed EE 1' : 182.780324107490, +'closed EE 2' : 182.780324107490, + }, + 'ClH.ezfio' : { + 'closed EE 1' : 182.691147725717, +'closed EE 2' : 182.691147725717, + }, + 'ClH+.ezfio' : { + 'closed EE 1' : 165.325459132812, +'closed EE 2' : 165.325459132812, + }, + 'ClS.ezfio' : { + 'closed EE 1' : 385.435861712905, +'closed EE 2' : 385.435861712905, + }, + 'ClSiH3.ezfio' : { + 'closed EE 1' : 372.487466205289, +'closed EE 2' : 372.487466205289, + }, + 'CN-.ezfio' : { + 'closed EE 1' : 55.3766122391237, +'closed EE 2' : 55.3766122391237, + }, + 'CN.ezfio' : { + 'closed EE 1' : 45.1366654659749, +'closed EE 2' : 45.1366654659749, + }, + 'CO2.ezfio' : { + 'closed EE 1' : 126.003917383945, +'closed EE 2' : 126.003917383945, + }, + 'CO.ezfio' : { + 'closed EE 1' : 62.7345020629198, +'closed EE 2' : 62.7345020629198, + }, + 'CO+.ezfio' : { + 'closed EE 1' : 52.0257806106742, +'closed EE 2' : 52.0257806106742, + }, + 'COS.ezfio' : { + 'closed EE 1' : 270.518850152602, +'closed EE 2' : 270.518850152602, + }, + 'CP.ezfio' : { + 'closed EE 1' : 165.996004819480, +'closed EE 2' : 165.996004819480, + }, + 'CS2.ezfio' : { + 'closed EE 1' : 423.926547892897, +'closed EE 2' : 423.926547892897, + }, + 'CS.ezfio' : { + 'closed EE 1' : 196.955019712066, +'closed EE 2' : 196.955019712066, + }, + 'CS+.ezfio' : { + 'closed EE 1' : 183.101079356045, +'closed EE 2' : 183.101079356045, + }, + 'CSi.ezfio' : { + 'closed EE 1' : 136.339370158432, +'closed EE 2' : 136.339370158432, + }, + 'F2.ezfio' : { + 'closed EE 1' : 109.589586860957, +'closed EE 2' : 109.589586860957, + }, + 'FAl.ezfio' : { + 'closed EE 1' : 172.329074560241, +'closed EE 2' : 172.329074560241, + }, + 'FCl.ezfio' : { + 'closed EE 1' : 263.250821798940, +'closed EE 2' : 263.250821798940, + }, + 'FCl+.ezfio' : { + 'closed EE 1' : 243.398671131879, +'closed EE 2' : 243.398671131879, + }, + 'F-.ezfio' : { + 'closed EE 1' : 45.5127217885110, +'closed EE 2' : 45.5127217885110, + }, + 'F.ezfio' : { + 'closed EE 1' : 32.5100994559105, +'closed EE 2' : 32.5100994559105, + }, + 'F+.ezfio' : { + 'closed EE 1' : 21.0616333765295, +'closed EE 2' : 21.0616333765295, + }, + 'FH.ezfio' : { + 'closed EE 1' : 45.4165968585922, +'closed EE 2' : 45.4165968585922, + }, + 'FH+.ezfio' : { + 'closed EE 1' : 32.0966404042882, +'closed EE 2' : 32.0966404042882, + }, + 'FMg.ezfio' : { + 'closed EE 1' : 147.723081703366, +'closed EE 2' : 147.723081703366, + }, + 'FNa.ezfio' : { + 'closed EE 1' : 136.229809908309, +'closed EE 2' : 136.229809908309, + }, + 'FP.ezfio' : { + 'closed EE 1' : 195.769918393858, +'closed EE 2' : 195.769918393858, + }, + 'FS.ezfio' : { + 'closed EE 1' : 226.988766049930, +'closed EE 2' : 226.988766049930, + }, + 'FSi.ezfio' : { + 'closed EE 1' : 184.557018579786, +'closed EE 2' : 184.557018579786, + }, + 'FSiH3.ezfio' : { + 'closed EE 1' : 213.848410787899, +'closed EE 2' : 213.848410787899, + }, + 'H2CNH.ezfio' : { + 'closed EE 1' : 64.8118892917432, +'closed EE 2' : 64.8118892917432, + }, + 'H2CO.ezfio' : { + 'closed EE 1' : 72.1212876559760, +'closed EE 2' : 72.1212876559760, + }, + 'H2CPH.ezfio' : { + 'closed EE 1' : 191.727599087985, +'closed EE 2' : 191.727599087985, + }, + 'H2CS.ezfio' : { + 'closed EE 1' : 209.231104868214, +'closed EE 2' : 209.231104868214, + }, + 'H2.ezfio' : { + 'closed EE 1' : 0.658003361873990, +'closed EE 2' : 0.658003361873990, + }, + 'H2NNH2.ezfio' : { + 'closed EE 1' : 80.4177432036452, +'closed EE 2' : 80.4177432036452, + }, + 'H2PPH2.ezfio' : { + 'closed EE 1' : 348.965444263036, +'closed EE 2' : 348.965444263036, + }, + 'H3CCl.ezfio' : { + 'closed EE 1' : 239.272625318642, +'closed EE 2' : 239.272625318642, + }, + 'H3CF.ezfio' : { + 'closed EE 1' : 89.9110566746969, +'closed EE 2' : 89.9110566746969, + }, + 'H3CNH2.ezfio' : { + 'closed EE 1' : 74.0553942748957, +'closed EE 2' : 74.0553942748957, + }, + 'H3COH.ezfio' : { + 'closed EE 1' : 81.4464317182977, +'closed EE 2' : 81.4464317182977, + }, + 'H3CSH.ezfio' : { + 'closed EE 1' : 220.631290723648, +'closed EE 2' : 220.631290723648, + }, + 'H3SiSiH3.ezfio' : { + 'closed EE 1' : 313.726268238800, +'closed EE 2' : 313.726268238800, + }, + 'HBO.ezfio' : { + 'closed EE 1' : 58.0559325869888, +'closed EE 2' : 58.0559325869888, + }, + 'HBS.ezfio' : { + 'closed EE 1' : 191.322628710909, +'closed EE 2' : 191.322628710909, + }, + 'HCF.ezfio' : { + 'closed EE 1' : 79.0947411389797, +'closed EE 2' : 79.0947411389797, + }, + 'HCN.ezfio' : { + 'closed EE 1' : 55.7130916664010, +'closed EE 2' : 55.7130916664010, + }, + 'HCO.ezfio' : { + 'closed EE 1' : 60.6786870182187, +'closed EE 2' : 60.6786870182187, + }, + 'HCP.ezfio' : { + 'closed EE 1' : 180.296919079404, +'closed EE 2' : 180.296919079404, + }, + 'H.ezfio' : { + 'closed EE 1' : 0.000000000000000E+000, +'closed EE 2' : 0.000000000000000E+000, + }, + 'HNO.ezfio' : { + 'closed EE 1' : 77.5573550103421, +'closed EE 2' : 77.5573550103421, + }, + 'HOCl.ezfio' : { + 'closed EE 1' : 253.799452346657, +'closed EE 2' : 253.799452346657, + }, + 'HOF.ezfio' : { + 'closed EE 1' : 101.663005370608, +'closed EE 2' : 101.663005370608, + }, + 'HOMg.ezfio' : { + 'closed EE 1' : 138.443726972836, +'closed EE 2' : 138.443726972836, + }, + 'HONa.ezfio' : { + 'closed EE 1' : 126.747456382832, +'closed EE 2' : 126.747456382832, + }, + 'HOO.ezfio' : { + 'closed EE 1' : 79.4438064643797, +'closed EE 2' : 79.4438064643797, + }, + 'HOOH.ezfio' : { + 'closed EE 1' : 93.6527230199755, +'closed EE 2' : 93.6527230199755, + }, + 'HSSH.ezfio' : { + 'closed EE 1' : 387.218345293825, +'closed EE 2' : 387.218345293825, + }, + 'Li2.ezfio' : { + 'closed EE 1' : 6.46521052077037, +'closed EE 2' : 6.46521052077037, + }, + 'LiCl.ezfio' : { + 'closed EE 1' : 193.091759412769, +'closed EE 2' : 193.091759412769, + }, + 'Li.ezfio' : { + 'closed EE 1' : 1.64991166351969, +'closed EE 2' : 1.64991166351969, + }, + 'Li+.ezfio' : { + 'closed EE 1' : 1.65066035543863, +'closed EE 2' : 1.65066035543863, + }, + 'LiF.ezfio' : { + 'closed EE 1' : 53.4448357111712, +'closed EE 2' : 53.4448357111712, + }, + 'LiH.ezfio' : { + 'closed EE 1' : 3.45968704863106, +'closed EE 2' : 3.45968704863106, + }, + 'LiN.ezfio' : { + 'closed EE 1' : 19.4901024446162, +'closed EE 2' : 19.4901024446162, + }, + 'LiO.ezfio' : { + 'closed EE 1' : 33.6641656402794, +'closed EE 2' : 33.6641656402794, + }, + 'LiOH.ezfio' : { + 'closed EE 1' : 45.3987062615096, +'closed EE 2' : 45.3987062615096, + }, + 'MgCl.ezfio' : { + 'closed EE 1' : 297.457458920467, +'closed EE 2' : 297.457458920467, + }, + 'Mg.ezfio' : { + 'closed EE 1' : 79.8152268844060, +'closed EE 2' : 79.8152268844060, + }, + 'Mg+.ezfio' : { + 'closed EE 1' : 72.1282818471152, +'closed EE 2' : 72.1282818471152, + }, + 'MgH.ezfio' : { + 'closed EE 1' : 79.3190204144503, +'closed EE 2' : 79.3190204144503, + }, + 'MgS.ezfio' : { + 'closed EE 1' : 279.954838755823, +'closed EE 2' : 279.954838755823, + }, + 'N2.ezfio' : { + 'closed EE 1' : 61.5022395201583, +'closed EE 2' : 61.5022395201583, + }, + 'N2+.ezfio' : { + 'closed EE 1' : 48.5010781414573, +'closed EE 2' : 48.5010781414573, + }, + 'Na2.ezfio' : { + 'closed EE 1' : 153.459222056420, +'closed EE 2' : 153.459222056420, + }, + 'NaCl.ezfio' : { + 'closed EE 1' : 285.552878597068, +'closed EE 2' : 285.552878597068, + }, + 'Na.ezfio' : { + 'closed EE 1' : 63.1438465495254, +'closed EE 2' : 63.1438465495254, + }, + 'Na+.ezfio' : { + 'closed EE 1' : 63.1844537555369, +'closed EE 2' : 63.1844537555369, + }, + 'NaH.ezfio' : { + 'closed EE 1' : 69.3524887715318, +'closed EE 2' : 69.3524887715318, + }, + 'N.ezfio' : { + 'closed EE 1' : 8.56327997038596, +'closed EE 2' : 8.56327997038596, + }, + 'N+.ezfio' : { + 'closed EE 1' : 8.76869142211657, +'closed EE 2' : 8.76869142211657, + }, + 'NF.ezfio' : { + 'closed EE 1' : 69.0051111837754, +'closed EE 2' : 69.0051111837754, + }, + 'NH2-.ezfio' : { + 'closed EE 1' : 30.6309164550511, +'closed EE 2' : 30.6309164550511, + }, + 'NH2.ezfio' : { + 'closed EE 1' : 21.9260220960358, +'closed EE 2' : 21.9260220960358, + }, + 'NH3.ezfio' : { + 'closed EE 1' : 31.4406045429538, +'closed EE 2' : 31.4406045429538, + }, + 'NH3+.ezfio' : { + 'closed EE 1' : 22.1802907966178, +'closed EE 2' : 22.1802907966178, + }, + 'NH4+.ezfio' : { + 'closed EE 1' : 31.8659857288766, +'closed EE 2' : 31.8659857288766, + }, + 'NH-.ezfio' : { + 'closed EE 1' : 21.6689399304262, +'closed EE 2' : 21.6689399304262, + }, + 'NH.ezfio' : { + 'closed EE 1' : 14.2929080048555, +'closed EE 2' : 14.2929080048555, + }, + 'NO-.ezfio' : { + 'closed EE 1' : 61.7198631648835, +'closed EE 2' : 61.7198631648835, + }, + 'NO.ezfio' : { + 'closed EE 1' : 64.6696341812441, +'closed EE 2' : 64.6696341812441, + }, + 'NP.ezfio' : { + 'closed EE 1' : 187.100108926969, +'closed EE 2' : 187.100108926969, + }, + 'NS.ezfio' : { + 'closed EE 1' : 200.128710693974, +'closed EE 2' : 200.128710693974, + }, + 'NSi.ezfio' : { + 'closed EE 1' : 155.279906480313, +'closed EE 2' : 155.279906480313, + }, + 'O2Cl.ezfio' : { + 'closed EE 1' : 301.747915430568, +'closed EE 2' : 301.747915430568, + }, + 'O2-.ezfio' : { + 'closed EE 1' : 79.2328963139684, +'closed EE 2' : 79.2328963139684, + }, + 'O2.ezfio' : { + 'closed EE 1' : 66.1914541099256, +'closed EE 2' : 66.1914541099256, + }, + 'O2+.ezfio' : { + 'closed EE 1' : 69.4144048026520, +'closed EE 2' : 69.4144048026520, + }, + 'O2S.ezfio' : { + 'closed EE 1' : 312.243487598150, +'closed EE 2' : 312.243487598150, + }, + 'O2Si.ezfio' : { + 'closed EE 1' : 257.255582658461, +'closed EE 2' : 257.255582658461, + }, + 'O3.ezfio' : { + 'closed EE 1' : 152.241108922056, +'closed EE 2' : 152.241108922056, + }, + 'OCl.ezfio' : { + 'closed EE 1' : 236.412413075442, +'closed EE 2' : 236.412413075442, + }, + 'O-.ezfio' : { + 'closed EE 1' : 26.9389565391700, +'closed EE 2' : 26.9389565391700, + }, + 'O.ezfio' : { + 'closed EE 1' : 17.7360424805971, +'closed EE 2' : 17.7360424805971, + }, + 'O+.ezfio' : { + 'closed EE 1' : 10.1017140483289, +'closed EE 2' : 10.1017140483289, + }, + 'OH2.ezfio' : { + 'closed EE 1' : 37.8440317785205, +'closed EE 2' : 37.8440317785205, + }, + 'OH2+.ezfio' : { + 'closed EE 1' : 26.6519594351168, +'closed EE 2' : 26.6519594351168, + }, + 'OH3+.ezfio' : { + 'closed EE 1' : 38.1275235387040, +'closed EE 2' : 38.1275235387040, + }, + 'OH-.ezfio' : { + 'closed EE 1' : 37.4117044437012, +'closed EE 2' : 37.4117044437012, + }, + 'OH.ezfio' : { + 'closed EE 1' : 26.7707296794756, +'closed EE 2' : 26.7707296794756, + }, + 'OH+.ezfio' : { + 'closed EE 1' : 17.2889916223005, +'closed EE 2' : 17.2889916223005, + }, + 'OMg.ezfio' : { + 'closed EE 1' : 138.801601522175, +'closed EE 2' : 138.801601522175, + }, + 'ONa.ezfio' : { + 'closed EE 1' : 112.111906727576, +'closed EE 2' : 112.111906727576, + }, + 'OP-.ezfio' : { + 'closed EE 1' : 188.573528250670, +'closed EE 2' : 188.573528250670, + }, + 'OP.ezfio' : { + 'closed EE 1' : 191.910766225339, +'closed EE 2' : 191.910766225339, + }, + 'OPH.ezfio' : { + 'closed EE 1' : 208.214788904911, +'closed EE 2' : 208.214788904911, + }, + 'OS.ezfio' : { + 'closed EE 1' : 203.101155235861, +'closed EE 2' : 203.101155235861, + }, + 'OSi.ezfio' : { + 'closed EE 1' : 179.456084306272, +'closed EE 2' : 179.456084306272, + }, + 'P2.ezfio' : { + 'closed EE 1' : 323.782785949801, +'closed EE 2' : 323.782785949801, + }, + 'P2+.ezfio' : { + 'closed EE 1' : 301.882738442470, +'closed EE 2' : 301.882738442470, + }, + 'PCl.ezfio' : { + 'closed EE 1' : 349.890715378447, +'closed EE 2' : 349.890715378447, + }, + 'P-.ezfio' : { + 'closed EE 1' : 122.657805807297, +'closed EE 2' : 122.657805807297, + }, + 'P.ezfio' : { + 'closed EE 1' : 111.445647659427, +'closed EE 2' : 111.445647659427, + }, + 'PH2-.ezfio' : { + 'closed EE 1' : 147.249988541779, +'closed EE 2' : 147.249988541779, + }, + 'PH2.ezfio' : { + 'closed EE 1' : 134.714659285095, +'closed EE 2' : 134.714659285095, + }, + 'PH2+.ezfio' : { + 'closed EE 1' : 136.440910105336, +'closed EE 2' : 136.440910105336, + }, + 'PH3.ezfio' : { + 'closed EE 1' : 148.202871342844, +'closed EE 2' : 148.202871342844, + }, + 'PH3+.ezfio' : { + 'closed EE 1' : 134.975246396259, +'closed EE 2' : 134.975246396259, + }, + 'PH4+.ezfio' : { + 'closed EE 1' : 148.989370984817, +'closed EE 2' : 148.989370984817, + }, + 'PH-.ezfio' : { + 'closed EE 1' : 134.328964674983, +'closed EE 2' : 134.328964674983, + }, + 'PH.ezfio' : { + 'closed EE 1' : 122.455536217569, +'closed EE 2' : 122.455536217569, + }, + 'PS.ezfio' : { + 'closed EE 1' : 337.453574046369, +'closed EE 2' : 337.453574046369, + }, + 'S2-.ezfio' : { + 'closed EE 1' : 366.675175225190, +'closed EE 2' : 366.675175225190, + }, + 'S2.ezfio' : { + 'closed EE 1' : 349.448868216573, +'closed EE 2' : 349.448868216573, + }, + 'S-.ezfio' : { + 'closed EE 1' : 149.700149385399, +'closed EE 2' : 149.700149385399, + }, + 'S.ezfio' : { + 'closed EE 1' : 135.897852467643, +'closed EE 2' : 135.897852467643, + }, + 'S+.ezfio' : { + 'closed EE 1' : 122.164730812357, +'closed EE 2' : 122.164730812357, + }, + 'SH2.ezfio' : { + 'closed EE 1' : 164.906914972667, +'closed EE 2' : 164.906914972667, + }, + 'SH2+.ezfio' : { + 'closed EE 1' : 149.596440849195, +'closed EE 2' : 149.596440849195, + }, + 'SH3+.ezfio' : { + 'closed EE 1' : 165.246577302836, +'closed EE 2' : 165.246577302836, + }, + 'SH-.ezfio' : { + 'closed EE 1' : 164.245517624070, +'closed EE 2' : 164.245517624070, + }, + 'SH.ezfio' : { + 'closed EE 1' : 149.754586732578, +'closed EE 2' : 149.754586732578, + }, + 'SH+.ezfio' : { + 'closed EE 1' : 135.225667181454, +'closed EE 2' : 135.225667181454, + }, + 'Si2.ezfio' : { + 'closed EE 1' : 252.066473929603, +'closed EE 2' : 252.066473929603, + }, + 'SiCl.ezfio' : { + 'closed EE 1' : 337.640041152210, +'closed EE 2' : 337.640041152210, + }, + 'Si-.ezfio' : { + 'closed EE 1' : 100.662878826237, +'closed EE 2' : 100.662878826237, + }, + 'Si.ezfio' : { + 'closed EE 1' : 101.098004923652, +'closed EE 2' : 101.098004923652, + }, + 'SiH2_1A1.ezfio' : { + 'closed EE 1' : 121.587026002297, +'closed EE 2' : 121.587026002297, + }, + 'SiH2_3B1.ezfio' : { + 'closed EE 1' : 109.761390873831, +'closed EE 2' : 109.761390873831, + }, + 'SiH2-.ezfio' : { + 'closed EE 1' : 120.021930833955, +'closed EE 2' : 120.021930833955, + }, + 'SiH3-.ezfio' : { + 'closed EE 1' : 131.352511652808, +'closed EE 2' : 131.352511652808, + }, + 'SiH3.ezfio' : { + 'closed EE 1' : 120.791742854126, +'closed EE 2' : 120.791742854126, + }, + 'SiH4.ezfio' : { + 'closed EE 1' : 132.863104887314, +'closed EE 2' : 132.863104887314, + }, + 'SiH4+.ezfio' : { + 'closed EE 1' : 121.801965974605, +'closed EE 2' : 121.801965974605, + }, + 'SiH-.ezfio' : { + 'closed EE 1' : 109.783147473057, +'closed EE 2' : 109.783147473057, + }, + 'SiH.ezfio' : { + 'closed EE 1' : 110.775916738616, +'closed EE 2' : 110.775916738616, + }, + 'SiS.ezfio' : { + 'closed EE 1' : 324.574435129198, +'closed EE 2' : 324.574435129198, + }, +} diff --git a/src/Bitmask/ASSUMPTIONS.rst b/src/Bitmask/ASSUMPTIONS.rst index b7cd9966..ac489a54 100644 --- a/src/Bitmask/ASSUMPTIONS.rst +++ b/src/Bitmask/ASSUMPTIONS.rst @@ -1,4 +1,6 @@ -* ``bit_kind_shift``, ``bit_kind_size`` and ``bit_kind`` are coherent:: +``bit_kind_shift``, ``bit_kind_size`` and ``bit_kind`` are coherent: + +.. code_block:: fortran 2**bit_kind_shift = bit_kind_size bit_kind = bit_kind_size / 8 diff --git a/src/Bitmask/README.rst b/src/Bitmask/README.rst index 44200c39..2a53efd0 100644 --- a/src/Bitmask/README.rst +++ b/src/Bitmask/README.rst @@ -24,7 +24,9 @@ Assumptions .. Do not edit this section. It was auto-generated from the .. ASSUMPTIONS.rst file. -* ``bit_kind_shift``, ``bit_kind_size`` and ``bit_kind`` are coherent:: +``bit_kind_shift``, ``bit_kind_size`` and ``bit_kind`` are coherent: + +.. code_block:: fortran 2**bit_kind_shift = bit_kind_size bit_kind = bit_kind_size / 8 diff --git a/src/MOGuess/ASSUMPTIONS.rst b/src/MOGuess/ASSUMPTIONS.rst new file mode 100644 index 00000000..e69de29b diff --git a/src/MOGuess/H_CORE_guess.irp.f b/src/MOGuess/H_CORE_guess.irp.f new file mode 100644 index 00000000..2c823baf --- /dev/null +++ b/src/MOGuess/H_CORE_guess.irp.f @@ -0,0 +1,11 @@ +program H_CORE_guess + implicit none + character*(64) :: label + mo_coef = ao_ortho_lowdin_coef + TOUCH mo_coef + label = "H_CORE_GUESS" + call mo_as_eigvectors_of_mo_matrix(mo_mono_elec_integral,size(mo_mono_elec_integral,1),size(mo_mono_elec_integral,2),label) + call save_mos + + +end diff --git a/src/MOGuess/Makefile b/src/MOGuess/Makefile new file mode 100644 index 00000000..b2ea1de1 --- /dev/null +++ b/src/MOGuess/Makefile @@ -0,0 +1,8 @@ +default: all + +# Define here all new external source files and objects.Don't forget to prefix the +# object files with IRPF90_temp/ +SRC= +OBJ= + +include $(QPACKAGE_ROOT)/src/Makefile.common diff --git a/src/MOGuess/README.rst b/src/MOGuess/README.rst new file mode 100644 index 00000000..14ade90f --- /dev/null +++ b/src/MOGuess/README.rst @@ -0,0 +1,4 @@ +============== +MOGuess Module +============== + diff --git a/src/MOGuess/localize.irp.f b/src/MOGuess/localize.irp.f new file mode 100644 index 00000000..53790e02 --- /dev/null +++ b/src/MOGuess/localize.irp.f @@ -0,0 +1,114 @@ +!TODO Ecrire un cholesky avec bitmask + + +subroutine localize_mos(mask, nint) + implicit none + use bitmasks + integer, intent(in) :: nint + integer(bit_kind), intent(in) :: mask(nint) + integer :: i,j,k,l + double precision, allocatable :: DM(:,:) + double precision, allocatable :: mo_coef_new(:,:), R(:,:) + integer :: n + !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: DM, mo_coef_new, R + integer :: rank + integer, parameter :: n_core = 2 + + allocate(R(mo_tot_num,mo_tot_num)) + allocate(DM(ao_num_align,ao_num)) + allocate(mo_coef_new(ao_num_align,mo_tot_num)) + n = ao_num + mo_coef_new = mo_coef + +BEGIN_TEMPLATE + DM = 0.d0 + if ($START < $END) then + do k=$START, $END + do j=1,n + !DEC$ VECTOR ALIGNED + do i=1,n + DM(i,j) += mo_coef_new(i,k)*mo_coef_new(j,k) + enddo + enddo + enddo + call cholesky_mo(n,$END-$START+1,DM,mo_coef_new(1,$START),size(mo_coef_new,1),-1.d0,rank) + endif +SUBST [ START, END ] + 1 ; n_core ;; +END_TEMPLATE + + deallocate(DM) + call find_rotation(mo_coef,ao_num_align,mo_coef_new,ao_num,R,mo_tot_num) + mo_coef = mo_coef_new + deallocate(mo_coef_new) + + double precision,allocatable :: mo_energy_new(:) + integer, allocatable :: iorder(:) + allocate(mo_energy_new(mo_tot_num),iorder(mo_tot_num)) + + do i=1,mo_tot_num + iorder(i) = i + mo_energy_new(i) = 0.d0 + do k=1,mo_tot_num + mo_energy_new(i) += R(k,i)*R(k,i)*mo_energy(k) + enddo + enddo + mo_energy = mo_energy_new + call dsort(mo_energy(1),iorder(1),n_core) + allocate (mo_coef_new(ao_num_align,mo_tot_num)) + mo_coef_new = mo_coef + do j=1,mo_tot_num + do i=1,ao_num + mo_coef(i,j) = mo_coef_new(i,iorder(j)) + enddo + enddo + deallocate (mo_coef_new,R) + deallocate(mo_energy_new,iorder) + mo_label = 'localized' + + SOFT_TOUCH mo_coef mo_energy mo_label +end + + + + + +subroutine cholesky_mo(n,m,P,C,LDC,tol_in,rank) + implicit none + integer, intent(in) :: n,m, LDC + double precision, intent(in) :: P(LDC,n) + double precision, intent(out) :: C(LDC,m) + double precision, intent(in) :: tol_in + integer, intent(out) :: rank + + integer :: info + integer :: i,k + integer :: ipiv(n) + double precision:: tol + double precision, allocatable :: W(:,:), work(:) + !DEC$ ATTRIBUTES ALIGN: 32 :: W + !DEC$ ATTRIBUTES ALIGN: 32 :: work + !DEC$ ATTRIBUTES ALIGN: 32 :: ipiv + + allocate(W(LDC,n),work(2*n)) + tol=tol_in + + info = 0 + do i=1,n + do k=1,i + W(i,k) = P(i,k) + enddo + do k=i+1,n + W(i,k) = 0. + enddo + enddo + call DPSTRF('L', n, W, LDC, ipiv, rank, tol, work, info ) + do i=1,n + do k=1,min(m,rank) + C(ipiv(i),k) = W(i,k) + enddo + enddo + + deallocate(W,work) +end + diff --git a/src/MOGuess/mo_ortho_lowdin.irp.f b/src/MOGuess/mo_ortho_lowdin.irp.f new file mode 100644 index 00000000..99063da2 --- /dev/null +++ b/src/MOGuess/mo_ortho_lowdin.irp.f @@ -0,0 +1,46 @@ + +BEGIN_PROVIDER [double precision, ao_ortho_lowdin_coef, (ao_num_align,ao_num)] + BEGIN_DOC +! matrix of the coefficients of the mos generated by the +! orthonormalization by the S^{-1/2} canonical transformation of the aos +! ao_ortho_lowdin_coef(i,j) = coefficient of the ith ao on the jth ao_ortho_lowdin orbital + END_DOC + implicit none + integer :: i,j,k,l + double precision :: tmp_matrix(ao_num_align,ao_num),accu + do i = 1, ao_num + do j = 1, ao_num + tmp_matrix(i,j) = 0.d0 + enddo + enddo + do i = 1, ao_num + tmp_matrix(i,i) = 1.d0 + enddo + call ortho_lowdin(ao_overlap,ao_num_align,ao_num,tmp_matrix,ao_num_align,ao_num) + do i = 1, ao_num + do j = 1, ao_num + ao_ortho_lowdin_coef(j,i) = tmp_matrix(i,j) + enddo + enddo +END_PROVIDER +BEGIN_PROVIDER [double precision, ao_ortho_lowdin_overlap, (ao_num_align,ao_num)] + BEGIN_DOC +! overlap matrix of the ao_ortho_lowdin +! supposed to be the Identity + END_DOC + implicit none + integer :: i,j,k,l + double precision :: accu + do i = 1, ao_num + do j = 1, ao_num + accu = 0.d0 + do k = 1, ao_num + do l = 1, ao_num + accu = accu + ao_ortho_lowdin_coef(i,k) * ao_ortho_lowdin_coef(j,l) * ao_overlap(k,l) + enddo + enddo + ao_ortho_lowdin_overlap(i,j) = accu + enddo + enddo + +END_PROVIDER diff --git a/src/MOGuess/pot_mo_ortho_lowdin_ints.irp.f b/src/MOGuess/pot_mo_ortho_lowdin_ints.irp.f new file mode 100644 index 00000000..3b1875f0 --- /dev/null +++ b/src/MOGuess/pot_mo_ortho_lowdin_ints.irp.f @@ -0,0 +1,25 @@ +BEGIN_PROVIDER [double precision, ao_ortho_lowdin_nucl_elec_integral, (mo_tot_num_align,mo_tot_num)] + implicit none + integer :: i1,j1,i,j + double precision :: c_i1,c_j1 + + ao_ortho_lowdin_nucl_elec_integral = 0.d0 + !$OMP PARALLEL DO DEFAULT(none) & + !$OMP PRIVATE(i,j,i1,j1,c_j1,c_i1) & + !$OMP SHARED(mo_tot_num,ao_num,ao_ortho_lowdin_coef, & + !$OMP ao_ortho_lowdin_nucl_elec_integral, ao_nucl_elec_integral) + do i = 1, mo_tot_num + do j = 1, mo_tot_num + do i1 = 1,ao_num + c_i1 = ao_ortho_lowdin_coef(i1,i) + do j1 = 1,ao_num + c_j1 = c_i1*ao_ortho_lowdin_coef(j1,j) + ao_ortho_lowdin_nucl_elec_integral(j,i) = ao_ortho_lowdin_nucl_elec_integral(j,i) + & + c_j1 * ao_nucl_elec_integral(j1,i1) + enddo + enddo + enddo + enddo + !$OMP END PARALLEL DO +END_PROVIDER + diff --git a/src/MOGuess/tests/Makefile b/src/MOGuess/tests/Makefile new file mode 100644 index 00000000..77bd84ba --- /dev/null +++ b/src/MOGuess/tests/Makefile @@ -0,0 +1,33 @@ +OPENMP =1 +PROFILE =0 +DEBUG = 0 + +IRPF90+= -I tests + +REF_FILES=$(subst %.irp.f, %.ref, $(wildcard *.irp.f)) + +.PHONY: clean executables serial_tests parallel_tests + +all: clean executables serial_tests parallel_tests + +parallel_tests: $(REF_FILES) + @echo ; echo " ---- Running parallel tests ----" ; echo + @OMP_NUM_THREADS=10 ${QPACKAGE_ROOT}/scripts/run_tests.py + +serial_tests: $(REF_FILES) + @echo ; echo " ---- Running serial tests ----" ; echo + @OMP_NUM_THREADS=1 ${QPACKAGE_ROOT}/scripts/run_tests.py + +executables: $(wildcard *.irp.f) veryclean + $(MAKE) -C .. + +%.ref: $(wildcard $(QPACKAGE_ROOT)/data/inputs/*.md5) executables + $(QPACKAGE_ROOT)/scripts/create_test_ref.sh $* + +clean: + $(MAKE) -C .. clean + +veryclean: + $(MAKE) -C .. veryclean + + diff --git a/src/Utils/Makefile b/src/Utils/Makefile index b2ea1de1..69065f55 100644 --- a/src/Utils/Makefile +++ b/src/Utils/Makefile @@ -2,7 +2,7 @@ default: all # Define here all new external source files and objects.Don't forget to prefix the # object files with IRPF90_temp/ -SRC= -OBJ= +SRC=map_module.f90 +OBJ=IRPF90_temp/map_module.o include $(QPACKAGE_ROOT)/src/Makefile.common diff --git a/src/Utils/map_module.f90 b/src/Utils/map_module.f90 new file mode 100644 index 00000000..ff6813e4 --- /dev/null +++ b/src/Utils/map_module.f90 @@ -0,0 +1,640 @@ +module map_module + + use omp_lib + + integer, parameter :: integral_kind = 8 + + integer, parameter :: cache_key_kind = 2 + integer, parameter :: cache_map_size_kind = 4 + + integer, parameter :: key_kind = 8 + integer, parameter :: map_size_kind = 8 + + integer, parameter :: map_shift = -15 + integer*8, parameter :: map_mask = ibset(0_8,15)-1_8 + + type cache_map_type + integer(cache_key_kind), pointer :: key(:) + real(integral_kind), pointer :: value(:) + logical :: sorted + integer(cache_map_size_kind) :: map_size + integer(cache_map_size_kind) :: n_elements + integer(omp_lock_kind) :: lock + end type cache_map_type + + type map_type + type(cache_map_type), allocatable :: map(:) + integer(map_size_kind) :: map_size + integer(map_size_kind) :: n_elements + logical :: sorted + integer(omp_lock_kind) :: lock + end type map_type + +end module map_module + + +real function map_mb(map) + use map_module + use omp_lib + implicit none + type (map_type), intent(in) :: map + integer(map_size_kind) :: i + + map_mb = 8+map_size_kind+map_size_kind+omp_lock_kind+4 + do i=0,map%map_size + map_mb = map_mb + map%map(i)%map_size*(cache_key_kind+integral_kind) + & + 8+8+4+cache_map_size_kind+cache_map_size_kind+omp_lock_kind + enddo + map_mb = map_mb / (1024.d0*1024.d0) +end + +subroutine cache_map_init(map,sze) + use map_module + implicit none + type (cache_map_type), intent(inout) :: map + integer(cache_map_size_kind) :: sze + call omp_init_lock(map%lock) + call omp_set_lock(map%lock) + map%n_elements = 0_8 + map%map_size = 0_8 + map%sorted = .True. + NULLIFY(map%value, map%key) + call cache_map_reallocate(map,sze) + call omp_unset_lock(map%lock) + +end + +subroutine map_init(map,keymax) + use map_module + implicit none + integer*8, intent(in) :: keymax + type (map_type), intent(inout) :: map + integer(map_size_kind) :: i + integer(cache_map_size_kind) :: sze + integer :: err + + call omp_init_lock(map%lock) + call omp_set_lock(map%lock) + + map%n_elements = 0_8 + map%map_size = ishft(keymax,map_shift) + + allocate(map%map(0_8:map%map_size),stat=err) + if (err /= 0) then + print *, 'Unable to allocate map' + stop 5 + endif +!sze = max(sqrt(map%map_size/16.d0),2048.d0) + sze = 2 + do i=0_8,map%map_size + call cache_map_init(map%map(i),sze) + enddo + map%sorted = .True. + + call omp_unset_lock(map%lock) + +end + +subroutine cache_map_reallocate(map,sze) + use map_module + implicit none + integer(cache_map_size_kind), intent(in) :: sze + type (cache_map_type), intent(inout) :: map + + integer(cache_key_kind), pointer :: key_new(:) + real(integral_kind), pointer :: value_new(:) + integer(map_size_kind) :: i + integer :: err + !DIR$ ATTRIBUTES ALIGN : 64 :: key_new, value_new + + if (sze < map%n_elements) then + print *, 'Unable to resize map : map too large' + stop 3 + endif + + ! Resize keys + allocate( key_new(sze), stat=err ) + if (err /= 0) then + print *, 'Unable to allocate map', sze + stop 1 + endif + if (associated(map%key)) then + do i=1_8,min(size(map%key),map%n_elements) + key_new(i) = map%key(i) + enddo + deallocate(map%key) + endif + + ! Resize values + allocate( value_new(sze), stat=err ) + if (err /= 0) then + print *, 'Unable to allocate map', sze + stop 2 + endif + if (associated(map%value)) then + do i=1_8,min(size(map%key),map%n_elements) + value_new(i) = map%value(i) + enddo + deallocate(map%value) + endif + + ! Set new pointers + map%key => key_new + map%value => value_new + map%map_size = sze + +end + + +subroutine cache_map_deinit(map) + use map_module + implicit none + type (cache_map_type), intent(inout) :: map + + integer :: err + + if (associated( map % value )) then + deallocate( map % value, stat=err ) + if (err /= 0) then + print *, 'Unable to deallocate map' + stop 2 + endif + NULLIFY(map%value) + endif + + if (associated( map % key )) then + deallocate( map % key, stat=err ) + if (err /= 0) then + print *, 'Unable to deallocate map' + stop 4 + endif + NULLIFY(map%key) + endif + + map%n_elements = 0_8 + map%map_size = 0_8 + call omp_destroy_lock(map%lock) +end + +subroutine map_deinit(map) + use map_module + implicit none + type (map_type), intent(inout) :: map + integer :: err + integer(map_size_kind) :: i + + if (allocated( map % map )) then + do i=0_8, map%map_size + call cache_map_deinit(map%map(i)) + enddo + deallocate( map % map, stat=err ) + if (err /= 0) then + print *, 'Unable to deallocate map' + stop 6 + endif + endif + + map%n_elements = 0_8 + map%map_size = 0_8 + call omp_destroy_lock(map%lock) +end + +subroutine cache_map_sort(map) + use map_module + implicit none + type (cache_map_type), intent(inout) :: map + integer(cache_map_size_kind), allocatable :: iorder(:) + integer(cache_map_size_kind) :: i + !DIR$ ATTRIBUTES ALIGN : 64 :: iorder + + if (.not.map%sorted) then + allocate(iorder(map%n_elements)) + do i=1,map%n_elements + iorder(i) = i + enddo + if (cache_key_kind == 2) then + call i2radix_sort(map%key,iorder,map%n_elements,-1) + else if (cache_key_kind == 4) then + call iradix_sort(map%key,iorder,map%n_elements,-1) + else if (cache_key_kind == 8) then + call i8radix_sort(map%key,iorder,map%n_elements,-1) + endif + if (integral_kind == 4) then + call set_order(map%value,iorder,map%n_elements) + else if (integral_kind == 8) then + call dset_order(map%value,iorder,map%n_elements) + endif + deallocate(iorder) + map%sorted = .True. + endif + +end + +subroutine map_sort(map) + use map_module + implicit none + type (map_type), intent(inout) :: map + integer(map_size_kind) :: i + + if (.not.map%sorted) then + !$OMP PARALLEL DO SCHEDULE(static,1024) DEFAULT(SHARED) PRIVATE(i) + do i=0_8,map%map_size + call omp_set_lock(map%map(i)%lock) + call cache_map_sort(map%map(i)) + call omp_unset_lock(map%map(i)%lock) + enddo + !$OMP END PARALLEL DO + map%sorted = .True. + endif + +end + +subroutine cache_map_unique(map) + use map_module + implicit none + type (cache_map_type), intent(inout) :: map + integer(cache_key_kind) :: prev_key + integer(cache_map_size_kind) :: i, j + + call cache_map_sort(map) + prev_key = -1_8 + j=0 + do i=1,map%n_elements + if (map%key(i) /= prev_key) then + j = j+1 + map%value(j) = map%value(i) + map%key(j) = map%key(i) + prev_key = map%key(i) + else + map%value(j) = map%value(j)+map%value(i) + endif + enddo + map%n_elements = j + +end + +subroutine cache_map_shrink(map,thr) + use map_module + implicit none + type (cache_map_type), intent(inout) :: map + real(integral_kind) , intent(in) :: thr + integer(cache_map_size_kind) :: i,j + + j=0 + do i=1,map%n_elements + if (abs(map%value(i)) > thr) then + j = j+1 + map%value(j) = map%value(i) + map%key(j) = map%key(i) + endif + enddo + map%n_elements = j + +end + +subroutine map_unique(map) + use map_module + implicit none + type (map_type), intent(inout) :: map + integer(map_size_kind) :: i + integer(map_size_kind) :: icount + + icount = 0_8 + !$OMP PARALLEL DO SCHEDULE(dynamic,1000) DEFAULT(SHARED) PRIVATE(i) & + !$OMP REDUCTION(+:icount) + do i=0_8,map%map_size + call omp_set_lock(map%map(i)%lock) + call cache_map_unique(map%map(i)) + call omp_unset_lock(map%map(i)%lock) + icount = icount + map%map(i)%n_elements + enddo + !$OMP END PARALLEL DO + map%n_elements = icount + +end + +subroutine map_shrink(map,thr) + use map_module + implicit none + type (map_type), intent(inout) :: map + type (map_type), intent(in) :: thr + integer(map_size_kind) :: i + integer(map_size_kind) :: icount + + icount = 0_8 + !$OMP PARALLEL DO SCHEDULE(dynamic,1000) DEFAULT(SHARED) PRIVATE(i) & + !$OMP REDUCTION(+:icount) + do i=0_8,map%map_size + call omp_set_lock(map%map(i)%lock) + call cache_map_shrink(map%map(i),thr) + call omp_unset_lock(map%map(i)%lock) + icount = icount + map%map(i)%n_elements + enddo + !$OMP END PARALLEL DO + map%n_elements = icount + +end + +subroutine map_update(map, key, value, sze, thr) + use map_module + implicit none + type (map_type), intent(inout) :: map + integer, intent(in) :: sze + integer(key_kind), intent(inout) :: key(sze) + real(integral_kind), intent(inout):: value(sze) + real(integral_kind), intent(in) :: thr + + integer :: i + integer(map_size_kind) :: idx_cache, idx_cache_new + integer(cache_map_size_kind) :: idx + integer :: sze2 + integer(cache_key_kind) :: cache_key + integer(map_size_kind) :: n_elements_temp + type (cache_map_type) :: local_map + logical :: map_sorted + + sze2 = sze + map_sorted = .True. + + n_elements_temp = 0_8 + n_elements_temp = n_elements_temp + 1_8 + do while (sze2>0) + i=1 + do while (i<=sze) + if (key(i) /= 0_8) then + idx_cache = ishft(key(i),map_shift) + if (omp_test_lock(map%map(idx_cache)%lock)) then + local_map%key => map%map(idx_cache)%key + local_map%value => map%map(idx_cache)%value + local_map%sorted = map%map(idx_cache)%sorted + local_map%map_size = map%map(idx_cache)%map_size + local_map%n_elements = map%map(idx_cache)%n_elements + do + !DIR$ FORCEINLINE + call search_key_big_interval(key(i),local_map%key, local_map%n_elements, idx, 1_8, local_map%n_elements) + if (idx > 0_8) then + local_map%value(idx) = local_map%value(idx) + value(i) + else + ! Assert that the map has a proper size + if (local_map%n_elements == local_map%map_size) then + call cache_map_unique(local_map) + call cache_map_reallocate(local_map, local_map%n_elements + local_map%n_elements) + call cache_map_shrink(local_map,thr) + endif + cache_key = iand(key(i),map_mask) + local_map%n_elements = local_map%n_elements + 1_8 + local_map%value(local_map%n_elements) = value(i) + local_map%key(local_map%n_elements) = cache_key + local_map%sorted = .False. + n_elements_temp = n_elements_temp + 1_8 + endif ! idx > 0 + key(i) = 0_8 + i = i+1 + sze2 = sze2-1 + if (i>sze) then + i=1 + endif + if ( (ishft(key(i),map_shift) /= idx_cache).or.(key(i)==0_8)) then + exit + endif + enddo + map%map(idx_cache)%key => local_map%key + map%map(idx_cache)%value => local_map%value + map%map(idx_cache)%sorted = local_map%sorted + map%map(idx_cache)%n_elements = local_map%n_elements + map%map(idx_cache)%map_size = local_map%map_size + map_sorted = map_sorted .and. local_map%sorted + call omp_unset_lock(map%map(idx_cache)%lock) + endif ! omp_test_lock + else + i=i+1 + endif ! key = 0 + enddo ! i + enddo ! sze2 > 0 + call omp_set_lock(map%lock) + map%n_elements = map%n_elements + n_elements_temp + map%sorted = map%sorted .and. map_sorted + call omp_unset_lock(map%lock) + +end + +subroutine map_append(map, key, value, sze) + use map_module + implicit none + type (map_type), intent(inout) :: map + integer, intent(in) :: sze + integer(key_kind), intent(inout) :: key(sze) + real(integral_kind), intent(inout):: value(sze) + + integer :: i + integer(cache_map_size_kind) :: n_elements + integer(map_size_kind) :: idx_cache + integer(cache_key_kind) :: cache_key + + do i=1,sze + idx_cache = ishft(key(i),map_shift) + call omp_set_lock(map%map(idx_cache)%lock) + n_elements = map%map(idx_cache)%n_elements + 1 + ! Assert that the map has a proper size + if (n_elements == map%map(idx_cache)%map_size) then + call cache_map_reallocate(map%map(idx_cache), n_elements+ ishft(n_elements,-1)) + endif + cache_key = iand(key(i),map_mask) + map%map(idx_cache)%value(n_elements) = value(i) + map%map(idx_cache)%key(n_elements) = cache_key + map%map(idx_cache)%n_elements = n_elements + if (map%map(idx_cache)%sorted.and.n_elements > 1) then + map%map(idx_cache)%sorted = (map%map(idx_cache)%key(n_elements-1) <= cache_key) + map%sorted = map%sorted .and. map%map(idx_cache)%sorted + endif + call omp_unset_lock(map%map(idx_cache)%lock) + enddo + call omp_set_lock(map%lock) + map%n_elements = map%n_elements + sze + call omp_unset_lock(map%lock) + +end + +subroutine map_get(map, key, value) + use map_module + implicit none + type (map_type), intent(in) :: map + integer(key_kind), intent(in) :: key + real(integral_kind), intent(out) :: value + integer(map_size_kind) :: idx_cache + integer(cache_map_size_kind) :: idx + + idx_cache = ishft(key,map_shift) + !DIR$ FORCEINLINE + call cache_map_get_interval(map%map(idx_cache), key, value, 1, map%map(idx_cache)%n_elements,idx) +end + +subroutine cache_map_get_interval(map, key, value, ibegin, iend, idx) + use map_module + implicit none + type (cache_map_type), intent(inout) :: map + integer(key_kind), intent(in) :: key + integer(cache_map_size_kind), intent(in) :: ibegin, iend + real(integral_kind), intent(out) :: value + integer(cache_map_size_kind), intent(in) :: idx + + call search_key_big_interval(key,map%key, map%n_elements, idx, ibegin, iend) + if (idx > 0) then + value = map%value(idx) + else + value = 0. + endif +end + + +subroutine map_get_many(map, key, value, sze) + use map_module + implicit none + type (map_type), intent(inout) :: map + integer, intent(in) :: sze + integer(key_kind), intent(in) :: key(sze) + real(integral_kind), intent(out) :: value(sze) + integer :: i + integer(map_size_kind) :: idx_cache, idx_cache_prev + integer(cache_map_size_kind) :: ibegin, iend + integer(cache_map_size_kind), allocatable :: idx(:) + !DIR$ ATTRIBUTES ALIGN : 64 :: idx + + allocate(idx(sze)) + do i=1,sze + idx_cache = ishft(key(i),map_shift) + iend = map%map(idx_cache)%n_elements + !DIR$ FORCEINLINE + call search_key_big_interval(key(i),map%map(idx_cache)%key, iend, idx(i), 1, iend) + enddo + do i=1,sze + idx_cache = ishft(key(i),map_shift) + if (idx(i) > 0) then + value(i) = map%map(idx_cache)%value(idx(i)) + else + value(i) = 0. + endif + enddo + deallocate(idx) +end + + +subroutine search_key_big(key,X,sze,idx) + use map_module + implicit none + integer(cache_map_size_kind), intent(in) :: sze + integer(key_kind) , intent(in) :: key + integer(cache_key_kind) , intent(in) :: X(sze) + integer(cache_map_size_kind), intent(out) :: idx + + call search_key_big_interval(key,X,sze,idx,1,sze) +end + + +subroutine search_key_big_interval(key,X,sze,idx,ibegin_in,iend_in) + use map_module + implicit none + integer(cache_map_size_kind), intent(in) :: sze + integer(key_kind) , intent(in) :: key + integer(cache_key_kind) , intent(in) :: X(sze) + integer(cache_map_size_kind), intent(in) :: ibegin_in, iend_in + integer(cache_map_size_kind), intent(out) :: idx + + integer(cache_map_size_kind) :: istep, ibegin, iend, i + integer(cache_key_kind) :: cache_key + + if (sze /= 0) then + continue + else + idx = -1 + return + endif + cache_key = iand(key,map_mask) + ibegin = min(ibegin_in,sze) + iend = min(iend_in,sze) + if ((cache_key > X(ibegin)) .and. (cache_key < X(iend))) then + + istep = ishft(iend-ibegin,-1) + idx = ibegin + istep + do while (istep > 32) + idx = ibegin + istep + if (cache_key < X(idx)) then + iend = idx + istep = ishft(idx-ibegin,-1) + idx = ibegin + istep + if (cache_key < X(idx)) then + iend = idx + istep = ishft(idx-ibegin,-1) + cycle + else if (cache_key > X(idx)) then + ibegin = idx + istep = ishft(iend-idx,-1) + cycle + else + return + endif + else if (cache_key > X(idx)) then + ibegin = idx + istep = ishft(iend-idx,-1) + idx = idx + istep + if (cache_key < X(idx)) then + iend = idx + istep = ishft(idx-ibegin,-1) + cycle + else if (cache_key > X(idx)) then + ibegin = idx + istep = ishft(iend-idx,-1) + cycle + else + return + endif + else + return + endif + enddo + idx = ibegin + if (min(iend_in,sze) > ibegin+64) then + iend = ibegin+64 + !DIR$ VECTOR ALIGNED + do while (cache_key > X(idx)) + idx = idx+1 + end do + else + !DIR$ VECTOR ALIGNED + do while (cache_key > X(idx)) + idx = idx+1 + if (idx /= iend) then + cycle + else + exit + endif + end do + endif + if (cache_key /= X(idx)) then + idx = 1-idx + endif + return + + else + + if (cache_key < X(ibegin)) then + idx = -ibegin + return + endif + if (cache_key > X(iend)) then + idx = -iend + return + endif + if (cache_key == X(ibegin)) then + idx = ibegin + return + endif + if (cache_key == X(iend)) then + idx = iend + return + endif + endif + +end + +