diff --git a/src/AOs/tests/Makefile b/src/AOs/tests/Makefile new file mode 100644 index 00000000..77bd84ba --- /dev/null +++ b/src/AOs/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/README.rst b/src/BiInts/README.rst index b2cdffc7..48fa022b 100644 --- a/src/BiInts/README.rst +++ b/src/BiInts/README.rst @@ -107,7 +107,8 @@ Documentation AO integrals `bielec_integrals_index `_ -None + Undocumented + `clear_ao_map `_ Frees the memory of the AO map diff --git a/src/CISD/cisd.irp.f b/src/CISD/cisd.irp.f index 0b68124a..d2f6de94 100644 --- a/src/CISD/cisd.irp.f +++ b/src/CISD/cisd.irp.f @@ -1,11 +1,16 @@ program cisd implicit none - integer :: i + integer :: i,k double precision, allocatable :: eigvalues(:),eigvectors(:,:) + PROVIDE ref_bitmask_energy call H_apply_cisd - allocate(eigvalues(n_det),eigvectors(n_det,n_det)) + allocate(eigvalues(n_states),eigvectors(n_det,n_states)) print *, 'N_det = ', N_det - psi_coef = psi_coef - 1.d-4 + print *, 'N_states = ', N_states + psi_coef = - 1.d-4 + do k=1,N_states + psi_coef(k,k) = 1.d0 + enddo call davidson_diag(psi_det,psi_coef,eigvalues,size(psi_coef,1),N_det,N_states,N_int) print *, '---' diff --git a/src/Dets/H_apply_template.f b/src/Dets/H_apply_template.f index 21c73377..3c13d31f 100644 --- a/src/Dets/H_apply_template.f +++ b/src/Dets/H_apply_template.f @@ -26,7 +26,7 @@ subroutine $subroutine_diexc(key_in, hole_1,particl_1, hole_2, particl_2 $parame integer, allocatable :: ia_ja_pairs(:,:,:) double precision :: diag_H_mat_elem, E_ref - PROVIDE mo_integrals_map + PROVIDE mo_integrals_map ref_bitmask_energy PROVIDE mo_bielec_integrals_in_map $set_i_H_j_threshold @@ -257,7 +257,7 @@ subroutine $subroutine_monoexc(key_in, hole_1,particl_1 $parameters ) integer, allocatable :: ia_ja_pairs(:,:,:) double precision :: diag_H_mat_elem, E_ref - PROVIDE mo_integrals_map + PROVIDE mo_integrals_map ref_bitmask_energy PROVIDE mo_bielec_integrals_in_map $set_i_H_j_threshold diff --git a/src/Dets/davidson.irp.f b/src/Dets/davidson.irp.f index 928bcb72..780d54b5 100644 --- a/src/Dets/davidson.irp.f +++ b/src/Dets/davidson.irp.f @@ -12,7 +12,7 @@ BEGIN_PROVIDER [ integer, davidson_sze_max] ! Max number of Davidson sizes END_DOC ASSERT (davidson_sze_max <= davidson_iter_max) - davidson_sze_max = 4 + davidson_sze_max = 8 END_PROVIDER subroutine davidson_diag(dets_in,u_in,energies,dim_in,sze,N_st,Nint) @@ -50,15 +50,17 @@ subroutine davidson_diag(dets_in,u_in,energies,dim_in,sze,N_st,Nint) integer :: k_pairs, kl integer :: iter2 - double precision, allocatable :: W(:,:), H_jj(:), U(:,:,:), R(:,:) + double precision, allocatable :: W(:,:,:), H_jj(:), U(:,:,:), R(:,:) double precision, allocatable :: y(:,:,:,:), h(:,:,:,:), lambda(:) double precision :: diag_h_mat_elem double precision :: residual_norm(N_st) + PROVIDE ref_bitmask_energy + allocate( & kl_pairs(2,N_st*(N_st+1)/2), & H_jj(sze), & - W(sze,N_st), & + W(sze,N_st,davidson_sze_max), & U(sze,N_st,davidson_sze_max), & R(sze,N_st), & h(N_st,davidson_sze_max,N_st,davidson_sze_max), & @@ -109,6 +111,7 @@ subroutine davidson_diag(dets_in,u_in,energies,dim_in,sze,N_st,Nint) enddo !$OMP END DO !$OMP END PARALLEL + call ortho_lowdin(overlap,size(overlap,1),N_st,U_in,size(U_in,1),sze) ! Davidson iterations @@ -148,33 +151,24 @@ subroutine davidson_diag(dets_in,u_in,energies,dim_in,sze,N_st,Nint) ! ---------------------- do k=1,N_st - call H_u_0(W(1,k),U(1,k,iter),H_jj,sze,dets_in,Nint) + call H_u_0(W(1,k,iter),U(1,k,iter),H_jj,sze,dets_in,Nint) enddo ! Compute h_kl = = ! ------------------------------------------- - !$OMP PARALLEL - !$OMP SINGLE do l=1,N_st - do iter2=1,iter-1 - do k=1,N_st - !$OMP TASK FIRSTPRIVATE(k,iter,l,iter2) - h(k,iter2,l,iter) = u_dot_v(U(1,k,iter2),W(1,l),sze) + do k=1,N_st + do iter2=1,iter-1 + h(k,iter2,l,iter) = u_dot_v(U(1,k,iter2),W(1,l,iter),sze) h(k,iter,l,iter2) = h(k,iter2,l,iter) - !$OMP END TASK enddo enddo do k=1,l - !$OMP TASK FIRSTPRIVATE(k,iter,l) - h(k,iter,l,iter) = u_dot_v(U(1,k,iter),W(1,l),sze) + h(k,iter,l,iter) = u_dot_v(U(1,k,iter),W(1,l,iter),sze) h(l,iter,k,iter) = h(k,iter,l,iter) - !$OMP END TASK enddo enddo - !$OMP END SINGLE NOWAIT - !$OMP TASKWAIT - !$OMP END PARALLEL ! Diagonalize h ! ------------- @@ -184,13 +178,17 @@ subroutine davidson_diag(dets_in,u_in,energies,dim_in,sze,N_st,Nint) ! Express eigenvectors of h in the determinant basis ! -------------------------------------------------- - !TODO dgemm + ! call dgemm ( 'N','N', sze, N_st*iter, N_st, & + ! 1.d0, U(1,1,1), size(U,1), y(1,1,1,1), size(y,1)*size(y,2), & + ! 0.d0, U(1,1,iter+1), size(U,1) ) do k=1,N_st do i=1,sze U(i,k,iter+1) = 0.d0 - do iter2=1,iter - do l=1,N_st - U(i,k,iter+1) += U(i,l,iter2)*y(l,iter2,k,1) + W(i,k,iter+1) = 0.d0 + do l=1,N_st + do iter2=1,iter + U(i,k,iter+1) = U(i,k,iter+1) + U(i,l,iter2)*y(l,iter2,k,1) + W(i,k,iter+1) = W(i,k,iter+1) + W(i,l,iter2)*y(l,iter2,k,1) enddo enddo enddo @@ -199,13 +197,9 @@ subroutine davidson_diag(dets_in,u_in,energies,dim_in,sze,N_st,Nint) ! Compute residual vector ! ----------------------- - do k=1,N_st - call H_u_0(W(1,k),U(1,k,iter+1),H_jj,sze,dets_in,Nint) - enddo - do k=1,N_st do i=1,sze - R(i,k) = lambda(k) * U(i,k,iter+1) - W(i,k) + R(i,k) = lambda(k) * U(i,k,iter+1) - W(i,k,iter+1) enddo residual_norm(k) = u_dot_u(R(1,k),sze) enddo @@ -215,7 +209,7 @@ subroutine davidson_diag(dets_in,u_in,energies,dim_in,sze,N_st,Nint) print *, residual_norm(1:N_st) print *, '' - converged = maxval(residual_norm) < 1.d-10 + converged = maxval(residual_norm) < 1.d-5 if (converged) then exit endif @@ -253,7 +247,7 @@ subroutine davidson_diag(dets_in,u_in,energies,dim_in,sze,N_st,Nint) enddo if (.not.converged) then - iter = davidson_sze_max + iter = davidson_sze_max-1 endif ! Re-contract to u_in diff --git a/src/Dets/slater_rules.irp.f b/src/Dets/slater_rules.irp.f index 4e45f586..59dc68ad 100644 --- a/src/Dets/slater_rules.irp.f +++ b/src/Dets/slater_rules.irp.f @@ -964,26 +964,38 @@ subroutine H_u_0(v_0,u_0,H_jj,n,keys_tmp,Nint) integer, allocatable :: idx(:) double precision :: hij integer :: i,j,k,l, jj + integer :: i0, j0 ASSERT (Nint > 0) ASSERT (Nint == N_int) - ASSERT (n>0) + ASSERT (n>0) + PROVIDE ref_bitmask_energy + integer, parameter :: block_size = 157 !$OMP PARALLEL DEFAULT(NONE) & !$OMP PRIVATE(i,hij,j,k,idx,jj) SHARED(n,H_jj,u_0,keys_tmp,Nint)& !$OMP SHARED(v_0) - allocate(idx(0:n)) - !$OMP DO SCHEDULE(guided) + allocate(idx(0:block_size)) + !$OMP DO SCHEDULE(static) do i=1,n v_0(i) = H_jj(i) * u_0(i) - call filter_connected(keys_tmp,keys_tmp(1,1,i),Nint,n,idx) - do jj=1,idx(0) - j = idx(jj) - if (j/=i) then - call i_H_j(keys_tmp(1,1,j),keys_tmp(1,1,i),Nint,hij) - v_0(i) = v_0(i) + hij*u_0(j) - endif - enddo enddo !$OMP END DO + !$OMP DO SCHEDULE(guided) + do i0=1,n,block_size + do j0=1,n,block_size + do i=i0,min(i0+block_size-1,n) + call filter_connected(keys_tmp(1,1,j0),keys_tmp(1,1,i),Nint,min(block_size,i-j0+1),idx) + do jj=1,idx(0) + j = idx(jj)+j0-1 + if ( (j 1.d-8)) then + call i_H_j(keys_tmp(1,1,j),keys_tmp(1,1,i),Nint,hij) + v_0(i) = v_0(i) + hij*u_0(j) + v_0(j) = v_0(j) + hij*u_0(i) + endif + enddo + enddo + enddo + enddo + !$OMP END DO NOWAIT deallocate(idx) !$OMP END PARALLEL end diff --git a/src/Hartree_Fock/README.rst b/src/Hartree_Fock/README.rst index 94c85fc1..bd52e526 100644 --- a/src/Hartree_Fock/README.rst +++ b/src/Hartree_Fock/README.rst @@ -82,7 +82,8 @@ Documentation Diagonal Fock matrix in the MO basis `scf_iteration `_ -None + Undocumented + `do_diis `_ If True, compute integrals on the fly diff --git a/src/MOs/README.rst b/src/MOs/README.rst index 3e4981df..659db5d5 100644 --- a/src/MOs/README.rst +++ b/src/MOs/README.rst @@ -54,8 +54,10 @@ Documentation Aligned variable for dimensioning of arrays `mo_as_eigvectors_of_mo_matrix `_ -None + Undocumented + `save_mos `_ -None + Undocumented + diff --git a/src/MonoInts/kin_ao_ints.irp.f b/src/MonoInts/kin_ao_ints.irp.f index 444c3880..0a5cccd3 100644 --- a/src/MonoInts/kin_ao_ints.irp.f +++ b/src/MonoInts/kin_ao_ints.irp.f @@ -11,6 +11,7 @@ double precision :: A_center(3), B_center(3) integer :: power_A(3), power_B(3) double precision :: d_a_2,d_2 + PROVIDE all_utils dim1=100 BEGIN_DOC ! second derivatives matrix elements in the ao basis diff --git a/src/Utils/map_module.f90 b/src/Utils/map_module.f90 index ff6813e4..ba36e323 100644 --- a/src/Utils/map_module.f90 +++ b/src/Utils/map_module.f90 @@ -482,7 +482,7 @@ subroutine cache_map_get_interval(map, key, value, ibegin, iend, idx) if (idx > 0) then value = map%value(idx) else - value = 0. + value = 0._integral_kind endif end diff --git a/src/Utils/util.irp.f b/src/Utils/util.irp.f index e7116334..e25148f1 100644 --- a/src/Utils/util.irp.f +++ b/src/Utils/util.irp.f @@ -253,12 +253,16 @@ double precision function u_dot_u(u,sze) t3 = t2+t2 t4 = t3+t2 u_dot_u = 0.d0 - do i=1,t2 - u_dot_u = u_dot_u + u(t1+i)*u(t1+i) + u(t2+i)*u(t2+i) + & - u(t3+i)*u(t3+i) + u(t4+i)*u(t4+i) - enddo - do i=t4+t2+1,sze - u_dot_u = u_dot_u+u(i)*u(i) +! do i=1,t2 +! u_dot_u = u_dot_u + u(t1+i)*u(t1+i) + u(t2+i)*u(t2+i) + & +! u(t3+i)*u(t3+i) + u(t4+i)*u(t4+i) +! enddo +! do i=t4+t2+1,sze +! u_dot_u = u_dot_u+u(i)*u(i) +! enddo + + do i=1,sze + u_dot_u = u_dot_u + u(i)*u(i) enddo end