From 3979677a82b1008a30938d3b61e09e9e06b22f21 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Wed, 11 Jan 2017 15:04:00 +0100 Subject: [PATCH] MRSC2 no amplitudes --- plugins/Psiref_CAS/psi_ref.irp.f | 33 +++++++ plugins/mrsc2_no_amp/NEEDED_CHILDREN_MODULES | 1 + plugins/mrsc2_no_amp/README.rst | 12 +++ plugins/mrsc2_no_amp/mrsc2_no_amp.irp.f | 78 ++++++++++++++++ plugins/mrsc2_no_amp/sc2_no_amp.irp.f | 9 ++ src/Determinants/filter_connected.irp.f | 98 ++++++++++++++++++++ 6 files changed, 231 insertions(+) create mode 100644 plugins/mrsc2_no_amp/NEEDED_CHILDREN_MODULES create mode 100644 plugins/mrsc2_no_amp/README.rst create mode 100644 plugins/mrsc2_no_amp/mrsc2_no_amp.irp.f create mode 100644 plugins/mrsc2_no_amp/sc2_no_amp.irp.f diff --git a/plugins/Psiref_CAS/psi_ref.irp.f b/plugins/Psiref_CAS/psi_ref.irp.f index d3b6c28f..ab9e6943 100644 --- a/plugins/Psiref_CAS/psi_ref.irp.f +++ b/plugins/Psiref_CAS/psi_ref.irp.f @@ -67,3 +67,36 @@ END_PROVIDER END_PROVIDER + BEGIN_PROVIDER [double precision, norm_psi_ref, (N_states)] +&BEGIN_PROVIDER [double precision, inv_norm_psi_ref, (N_states)] + implicit none + integer :: i,j + norm_psi_ref = 0.d0 + do j = 1, N_states + do i = 1, N_det_ref + norm_psi_ref(j) += psi_ref_coef(i,j) * psi_ref_coef(i,j) + enddo + inv_norm_psi_ref(j) = 1.d0/(dsqrt(norm_psi_Ref(j))) + enddo + + END_PROVIDER + + BEGIN_PROVIDER [double precision, psi_ref_coef_interm_norm, (N_det_ref,N_states)] + implicit none + integer :: i,j + do j = 1, N_states + do i = 1, N_det_ref + psi_ref_coef_interm_norm(i,j) = inv_norm_psi_ref(j) * psi_ref_coef(i,j) + enddo + enddo + END_PROVIDER + + BEGIN_PROVIDER [double precision, psi_non_ref_coef_interm_norm, (N_det_non_ref,N_states)] + implicit none + integer :: i,j + do j = 1, N_states + do i = 1, N_det_non_ref + psi_non_ref_coef_interm_norm(i,j) = psi_non_ref_coef(i,j) * inv_norm_psi_ref(j) + enddo + enddo + END_PROVIDER diff --git a/plugins/mrsc2_no_amp/NEEDED_CHILDREN_MODULES b/plugins/mrsc2_no_amp/NEEDED_CHILDREN_MODULES new file mode 100644 index 00000000..f04fe3b0 --- /dev/null +++ b/plugins/mrsc2_no_amp/NEEDED_CHILDREN_MODULES @@ -0,0 +1 @@ +Psiref_CAS Determinants Davidson diff --git a/plugins/mrsc2_no_amp/README.rst b/plugins/mrsc2_no_amp/README.rst new file mode 100644 index 00000000..b24848f7 --- /dev/null +++ b/plugins/mrsc2_no_amp/README.rst @@ -0,0 +1,12 @@ +============ +mrsc2_no_amp +============ + +Needed Modules +============== +.. Do not edit this section It was auto-generated +.. by the `update_README.py` script. +Documentation +============= +.. Do not edit this section It was auto-generated +.. by the `update_README.py` script. diff --git a/plugins/mrsc2_no_amp/mrsc2_no_amp.irp.f b/plugins/mrsc2_no_amp/mrsc2_no_amp.irp.f new file mode 100644 index 00000000..b8b021e8 --- /dev/null +++ b/plugins/mrsc2_no_amp/mrsc2_no_amp.irp.f @@ -0,0 +1,78 @@ + BEGIN_PROVIDER [double precision, CI_eigenvectors_sc2_no_amp, (N_det,N_states_diag)] +&BEGIN_PROVIDER [double precision, CI_eigenvectors_s2_sc2_no_amp, (N_states_diag)] +&BEGIN_PROVIDER [double precision, CI_electronic_energy_sc2_no_amp, (N_states_diag)] + implicit none + integer :: i,j,k,l + integer, allocatable :: idx(:) + double precision, allocatable :: e_corr(:,:) + double precision, allocatable :: accu(:) + double precision, allocatable :: ihpsi_current(:) + double precision, allocatable :: H_jj(:),H_jj_total(:),S2_jj(:) + allocate(e_corr(N_det_non_ref,N_states),ihpsi_current(N_states),accu(N_states),H_jj(N_det_non_ref),idx(0:N_det_non_ref)) + allocate(H_jj_total(N_det),S2_jj(N_det)) + accu = 0.d0 + do i = 1, N_det_non_ref + call i_h_psi(psi_non_ref(1,1,i), psi_ref, psi_ref_coef_interm_norm, N_int, N_det_ref,& + size(psi_ref_coef_interm_norm,1), N_states,ihpsi_current) + do j = 1, N_states + e_corr(i,j) = psi_non_ref_coef_interm_norm(i,j) * ihpsi_current(j) + accu(j) += e_corr(i,j) + enddo + enddo + double precision :: hjj,diag_h_mat_elem + do i = 1, N_det_non_ref + call filter_not_connected(psi_non_ref,psi_non_ref(1,1,i),N_int,N_det_non_ref,idx) + H_jj(i) = 0.d0 + do j = 1, idx(0) + H_jj(i) += e_corr(idx(j),1) + enddo + enddo + do i=1,N_Det + H_jj_total(i) = diag_h_mat_elem(psi_det(1,1,i),N_int) + call get_s2(psi_det(1,1,i),psi_det(1,1,i),N_int,S2_jj(i)) + enddo + do i=1, N_det_non_ref + H_jj_total(idx_non_ref(i)) += H_jj(i) + enddo + + + call davidson_diag_hjj_sjj(psi_det,CI_eigenvectors_sc2_no_amp,H_jj_total,S2_jj,CI_electronic_energy_sc2_no_amp,size(CI_eigenvectors_sc2_no_amp,1),N_Det,N_states,N_states_diag,N_int,6) + do i=1,N_states_diag + CI_eigenvectors_s2_sc2_no_amp(i) = S2_jj(i) + enddo + + deallocate(e_corr,ihpsi_current,accu,H_jj,idx,H_jj_total,s2_jj) +END_PROVIDER + +BEGIN_PROVIDER [ double precision, CI_energy_sc2_no_amp, (N_states_diag) ] + implicit none + BEGIN_DOC + ! N_states lowest eigenvalues of the CI matrix + END_DOC + + integer :: j + character*(8) :: st + call write_time(output_determinants) + do j=1,min(N_det,N_states_diag) + CI_energy_sc2_no_amp(j) = CI_electronic_energy_sc2_no_amp(j) + nuclear_repulsion + enddo + do j=1,min(N_det,N_states) + write(st,'(I4)') j + call write_double(output_determinants,CI_energy_sc2_no_amp(j),'Energy of state '//trim(st)) + call write_double(output_determinants,CI_eigenvectors_s2_sc2_no_amp(j),'S^2 of state '//trim(st)) + enddo + +END_PROVIDER + +subroutine diagonalize_CI_sc2_no_amp + implicit none + integer :: i,j + do j=1,N_states + do i=1,N_det + psi_coef(i,j) = CI_eigenvectors_sc2_no_amp(i,j) + enddo + enddo + SOFT_TOUCH ci_eigenvectors_s2_sc2_no_amp ci_eigenvectors_sc2_no_amp ci_electronic_energy_sc2_no_amp ci_energy_sc2_no_amp psi_coef + +end + diff --git a/plugins/mrsc2_no_amp/sc2_no_amp.irp.f b/plugins/mrsc2_no_amp/sc2_no_amp.irp.f new file mode 100644 index 00000000..622d7236 --- /dev/null +++ b/plugins/mrsc2_no_amp/sc2_no_amp.irp.f @@ -0,0 +1,9 @@ +program pouet + implicit none + integer :: i + do i = 1, 10 + call diagonalize_CI_sc2_no_amp + TOUCH psi_coef + enddo + +end diff --git a/src/Determinants/filter_connected.irp.f b/src/Determinants/filter_connected.irp.f index da333b1e..b76540f7 100644 --- a/src/Determinants/filter_connected.irp.f +++ b/src/Determinants/filter_connected.irp.f @@ -1,4 +1,102 @@ +subroutine filter_not_connected(key1,key2,Nint,sze,idx) + use bitmasks + implicit none + BEGIN_DOC + ! Returns the array idx which contains the index of the + ! + ! determinants in the array key1 that DO NOT interact + ! + ! via the H operator with key2. + ! + ! idx(0) is the number of determinants that DO NOT interact with key1 + END_DOC + integer, intent(in) :: Nint, sze + integer(bit_kind), intent(in) :: key1(Nint,2,sze) + integer(bit_kind), intent(in) :: key2(Nint,2) + integer, intent(out) :: idx(0:sze) + + integer :: i,j,l + integer :: degree_x2 + + + ASSERT (Nint > 0) + ASSERT (sze >= 0) + + l=1 + + if (Nint==1) then + + !DIR$ LOOP COUNT (1000) + do i=1,sze + degree_x2 = popcnt( xor( key1(1,1,i), key2(1,1))) & + + popcnt( xor( key1(1,2,i), key2(1,2))) + if (degree_x2 > 4) then + idx(l) = i + l = l+1 + else + cycle + endif + enddo + + else if (Nint==2) then + + !DIR$ LOOP COUNT (1000) + do i=1,sze + degree_x2 = popcnt(xor( key1(1,1,i), key2(1,1))) + & + popcnt(xor( key1(2,1,i), key2(2,1))) + & + popcnt(xor( key1(1,2,i), key2(1,2))) + & + popcnt(xor( key1(2,2,i), key2(2,2))) + if (degree_x2 > 4) then + idx(l) = i + l = l+1 + else + cycle + endif + enddo + + else if (Nint==3) then + + !DIR$ LOOP COUNT (1000) + do i=1,sze + degree_x2 = popcnt(xor( key1(1,1,i), key2(1,1))) + & + popcnt(xor( key1(1,2,i), key2(1,2))) + & + popcnt(xor( key1(2,1,i), key2(2,1))) + & + popcnt(xor( key1(2,2,i), key2(2,2))) + & + popcnt(xor( key1(3,1,i), key2(3,1))) + & + popcnt(xor( key1(3,2,i), key2(3,2))) + if (degree_x2 > 4) then + idx(l) = i + l = l+1 + else + cycle + endif + enddo + + else + + !DIR$ LOOP COUNT (1000) + do i=1,sze + degree_x2 = 0 + !DEC$ LOOP COUNT MIN(4) + do j=1,Nint + degree_x2 = degree_x2+ popcnt(xor( key1(j,1,i), key2(j,1))) +& + popcnt(xor( key1(j,2,i), key2(j,2))) + if (degree_x2 > 4) then + idx(l) = i + l = l+1 + endif + enddo + if (degree_x2 <= 5) then + exit + endif + enddo + + endif + idx(0) = l-1 +end + + subroutine filter_connected(key1,key2,Nint,sze,idx) use bitmasks implicit none