From 7cb17e0a486efbc9aa20a5bb1a5a6a5ca5d7c865 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 13 Mar 2017 12:26:16 +0100 Subject: [PATCH] Super fast density matrix --- src/Determinants/density_matrix.irp.f | 87 +++++++++++++++++++------ src/Determinants/spindeterminants.irp.f | 48 +++++++++++++- 2 files changed, 113 insertions(+), 22 deletions(-) diff --git a/src/Determinants/density_matrix.irp.f b/src/Determinants/density_matrix.irp.f index ed2f49bd..e0ec6a6b 100644 --- a/src/Determinants/density_matrix.irp.f +++ b/src/Determinants/density_matrix.irp.f @@ -27,25 +27,33 @@ END_PROVIDER double precision :: ck, cl, ckl double precision :: phase integer :: h1,h2,p1,p2,s1,s2, degree + integer(bit_kind) :: tmp_det(N_int,2), tmp_det2(N_int,2) integer :: exc(0:2,2,2),n_occ(2) double precision, allocatable :: tmp_a(:,:,:), tmp_b(:,:,:) + integer :: krow, kcol, lrow, lcol one_body_dm_mo_alpha = 0.d0 one_body_dm_mo_beta = 0.d0 !$OMP PARALLEL DEFAULT(NONE) & !$OMP PRIVATE(j,k,l,m,occ,ck, cl, ckl,phase,h1,h2,p1,p2,s1,s2, degree,exc, & - !$OMP tmp_a, tmp_b, n_occ)& + !$OMP tmp_a, tmp_b, n_occ, krow, kcol, lrow, lcol, tmp_det, tmp_det2)& !$OMP SHARED(psi_det,psi_coef,N_int,N_states,elec_alpha_num,& !$OMP elec_beta_num,one_body_dm_mo_alpha,one_body_dm_mo_beta,N_det,mo_tot_num_align,& - !$OMP mo_tot_num) + !$OMP mo_tot_num,psi_bilinear_matrix_rows,psi_bilinear_matrix_columns, & + !$OMP psi_bilinear_matrix_transp_rows, psi_bilinear_matrix_transp_columns, & + !$OMP psi_bilinear_matrix_values, psi_bilinear_matrix_transp_values) allocate(tmp_a(mo_tot_num_align,mo_tot_num,N_states), tmp_b(mo_tot_num_align,mo_tot_num,N_states) ) tmp_a = 0.d0 tmp_b = 0.d0 - !$OMP DO SCHEDULE(dynamic) + !$OMP DO SCHEDULE(guided) do k=1,N_det - call bitstring_to_list_ab(psi_det(1,1,k), occ, n_occ, N_int) + krow = psi_bilinear_matrix_rows(k) + kcol = psi_bilinear_matrix_columns(k) + tmp_det(:,1) = psi_det(:,1, krow) + tmp_det(:,2) = psi_det(:,2, kcol) + call bitstring_to_list_ab(tmp_det, occ, n_occ, N_int) do m=1,N_states - ck = psi_coef(k,m)*psi_coef(k,m) + ck = psi_bilinear_matrix_values(k,m)*psi_bilinear_matrix_values(k,m) do l=1,elec_alpha_num j = occ(l,1) tmp_a(j,j,m) += ck @@ -55,24 +63,61 @@ END_PROVIDER tmp_b(j,j,m) += ck enddo enddo - do l=1,k-1 - call get_excitation_degree(psi_det(1,1,k),psi_det(1,1,l),degree,N_int) - if (degree /= 1) then - cycle + + l = k+1 + lrow = psi_bilinear_matrix_rows(l) + lcol = psi_bilinear_matrix_columns(l) + do while ( lcol == kcol ) + tmp_det2(:,1) = psi_det(:,1, lrow) + tmp_det2(:,2) = psi_det(:,2, lcol) + call get_excitation_degree(tmp_det,tmp_det2,degree,N_int) + if (degree == 1) then + call get_mono_excitation(psi_det(1,1,k),psi_det(1,1,l),exc,phase,N_int) + call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) + do m=1,N_states + ckl = psi_bilinear_matrix_values(k,m)*psi_bilinear_matrix_values(l,m) * phase + if (s1==1) then + tmp_a(h1,p1,m) += ckl + tmp_a(p1,h1,m) += ckl + else + tmp_b(h1,p1,m) += ckl + tmp_b(p1,h1,m) += ckl + endif + enddo endif - call get_mono_excitation(psi_det(1,1,k),psi_det(1,1,l),exc,phase,N_int) - call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) - do m=1,N_states - ckl = psi_coef(k,m) * psi_coef(l,m) * phase - if (s1==1) then - tmp_a(h1,p1,m) += ckl - tmp_a(p1,h1,m) += ckl - else - tmp_b(h1,p1,m) += ckl - tmp_b(p1,h1,m) += ckl - endif - enddo + l = l+1 + if (l>N_det) exit + lrow = psi_bilinear_matrix_rows(l) + lcol = psi_bilinear_matrix_columns(l) enddo + + l = k+1 + lrow = psi_bilinear_matrix_transp_rows(l) + lcol = psi_bilinear_matrix_transp_columns(l) + do while ( lrow == krow ) + tmp_det2(:,1) = psi_det(:,1, lrow) + tmp_det2(:,2) = psi_det(:,2, lcol) + call get_excitation_degree(tmp_det,tmp_det2,degree,N_int) + if (degree == 1) then + call get_mono_excitation(psi_det(1,1,k),psi_det(1,1,l),exc,phase,N_int) + call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) + do m=1,N_states + ckl = psi_bilinear_matrix_values(k,m)*psi_bilinear_matrix_transp_values(l,m) * phase + if (s1==1) then + tmp_a(h1,p1,m) += ckl + tmp_a(p1,h1,m) += ckl + else + tmp_b(h1,p1,m) += ckl + tmp_b(p1,h1,m) += ckl + endif + enddo + endif + l = l+1 + if (l>N_det) exit + lrow = psi_bilinear_matrix_transp_rows(l) + lcol = psi_bilinear_matrix_transp_columns(l) + enddo + enddo !$OMP END DO NOWAIT !$OMP CRITICAL diff --git a/src/Determinants/spindeterminants.irp.f b/src/Determinants/spindeterminants.irp.f index 2eec0dea..b30f84cb 100644 --- a/src/Determinants/spindeterminants.irp.f +++ b/src/Determinants/spindeterminants.irp.f @@ -393,6 +393,8 @@ BEGIN_PROVIDER [ double precision, psi_bilinear_matrix_values, (N_det,N_states) BEGIN_DOC ! Sparse coefficient matrix if the wave function is expressed in a bilinear form : ! D_a^t C D_b +! +! Rows are alpha determinants and columns are beta. END_DOC integer :: i,j,k, l integer(bit_kind) :: tmp_det(N_int,2) @@ -421,10 +423,54 @@ BEGIN_PROVIDER [ double precision, psi_bilinear_matrix_values, (N_det,N_states) call isort(to_sort, iorder, N_det) call iset_order(psi_bilinear_matrix_rows,iorder,N_det) call iset_order(psi_bilinear_matrix_columns,iorder,N_det) - call dset_order(psi_bilinear_matrix_values,iorder,N_det) + do l=1,N_states + call dset_order(psi_bilinear_matrix_values(1,l),iorder,N_det) + enddo deallocate(iorder,to_sort) END_PROVIDER + +BEGIN_PROVIDER [ double precision, psi_bilinear_matrix_transp_values, (N_det,N_states) ] +&BEGIN_PROVIDER [ integer, psi_bilinear_matrix_transp_rows, (N_det) ] +&BEGIN_PROVIDER [ integer, psi_bilinear_matrix_transp_columns, (N_det) ] + use bitmasks + implicit none + BEGIN_DOC +! Sparse coefficient matrix if the wave function is expressed in a bilinear form : +! D_a^t C D_b +! +! Rows are Beta determinants and columns are alpha + END_DOC + integer :: i,j,k,l + + + PROVIDE psi_coef_sorted_bit + + integer, allocatable :: iorder(:), to_sort(:) + allocate(iorder(N_det), to_sort(N_det)) + do l=1,N_states + do k=1,N_det + psi_bilinear_matrix_transp_values (k,l) = psi_bilinear_matrix_values (k,l) + enddo + enddo + do k=1,N_det + psi_bilinear_matrix_transp_columns(k) = psi_bilinear_matrix_columns(k) + psi_bilinear_matrix_transp_rows (k) = psi_bilinear_matrix_rows (k) + i = psi_bilinear_matrix_transp_columns(k) + j = psi_bilinear_matrix_transp_rows (k) + to_sort(k) = N_det_beta_unique * (j-1) + i + iorder(k) = k + enddo + call isort(to_sort, iorder, N_det) + call iset_order(psi_bilinear_matrix_transp_rows,iorder,N_det) + call iset_order(psi_bilinear_matrix_transp_columns,iorder,N_det) + do l=1,N_states + call dset_order(psi_bilinear_matrix_values(1,l),iorder,N_det) + enddo + deallocate(iorder,to_sort) +END_PROVIDER + + BEGIN_PROVIDER [ double precision, psi_bilinear_matrix, (N_det_alpha_unique,N_det_beta_unique,N_states) ] implicit none BEGIN_DOC