LCPQ/quantum_package
10
0
Fork 0
mirror of https://github.com/LCPQ/quantum_package synced 2024-12-23 04:43:50 +01:00
Code Issues Releases Wiki Activity

Acceleration of filter_connected and Davidson

Browse Source
This commit is contained in:
Anthony Scemama 2014-05-15 17:58:30 +02:00
parent df9efae335
commit 13ee3531ee
11 changed files with 108 additions and 55 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

33
src/AOs/tests/Makefile Normal file
View File

@ -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

3
src/BiInts/README.rst
View File

@ -107,7 +107,8 @@ Documentation
AO integrals
`bielec_integrals_index <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/map_integrals.irp.f#L17>`_
None
Undocumented
`clear_ao_map <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/map_integrals.irp.f#L128>`_
Frees the memory of the AO map

11
src/CISD/cisd.irp.f
View File

@ -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 *, '---'

4
src/Dets/H_apply_template.f
View File

@ -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

50
src/Dets/davidson.irp.f
View File

@ -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 = <u_k | W_l> = <u_k| H |u_l>
! -------------------------------------------
!$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

32
src/Dets/slater_rules.irp.f
View File

@ -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)
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<i).and.(dabs(u_0(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

3
src/Hartree_Fock/README.rst
View File

@ -82,7 +82,8 @@ Documentation
Diagonal Fock matrix in the MO basis
`scf_iteration <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/mo_SCF_iterations.irp.f#L1>`_
None
Undocumented
`do_diis <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/options.irp.f#L41>`_
If True, compute integrals on the fly

6
src/MOs/README.rst
View File

@ -54,8 +54,10 @@ Documentation
Aligned variable for dimensioning of arrays
`mo_as_eigvectors_of_mo_matrix <http://github.com/LCPQ/quantum_package/tree/master/src/MOs/utils.irp.f#L21>`_
None
Undocumented
`save_mos <http://github.com/LCPQ/quantum_package/tree/master/src/MOs/utils.irp.f#L1>`_
None
Undocumented

1
src/MonoInts/kin_ao_ints.irp.f
View File

@ -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

2
src/Utils/map_module.f90
View File

@ -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

16
src/Utils/util.irp.f
View File

@ -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