diff --git a/plugins/MRCC_Utils/mrcc_utils.irp.f b/plugins/MRCC_Utils/mrcc_utils.irp.f index 1e2f974d..866a6e3b 100644 --- a/plugins/MRCC_Utils/mrcc_utils.irp.f +++ b/plugins/MRCC_Utils/mrcc_utils.irp.f @@ -110,6 +110,7 @@ END_PROVIDER BEGIN_PROVIDER [ double precision, delta_ij, (N_det_ref,N_det_non_ref,N_states) ] &BEGIN_PROVIDER [ double precision, delta_ii, (N_det_ref,N_states) ] +&BEGIN_PROVIDER [ double precision, delta_cas, (N_det_ref,N_det_ref,N_states) ] implicit none BEGIN_DOC ! Dressing matrix in N_det basis @@ -117,6 +118,7 @@ END_PROVIDER integer :: i,j,m delta_ij = 0.d0 delta_ii = 0.d0 + call H_apply_mrcc(delta_ij,delta_ii,N_det_ref,N_det_non_ref) double precision :: max_delta double precision :: accu @@ -157,6 +159,17 @@ BEGIN_PROVIDER [ double precision, h_matrix_dressed, (N_det,N_det,N_states) ] h_matrix_dressed(i,j,istate) = h_matrix_all_dets(i,j) enddo enddo + + !!!!!!!!!! + do ii = 1, N_det_ref + do jj = 1, N_det_ref + i = idx_ref(ii) + j = idx_ref(jj) + h_matrix_dressed(i,j,istate) += delta_cas(ii,jj,istate) + h_matrix_dressed(j,i,istate) += delta_cas(ii,jj,istate) + end do + end do + !!!!!!!!!!!!! do ii = 1, N_det_ref i =idx_ref(ii) h_matrix_dressed(i,i,istate) += delta_ii(ii,istate) diff --git a/plugins/mrcepa0/dressing.irp.f b/plugins/mrcepa0/dressing.irp.f new file mode 100644 index 00000000..53a9417b --- /dev/null +++ b/plugins/mrcepa0/dressing.irp.f @@ -0,0 +1,234 @@ +use bitmasks + + + BEGIN_PROVIDER [ integer, cepa0_shortcut, (0:N_det_non_ref+1) ] +&BEGIN_PROVIDER [ integer, det_cepa0_idx, (N_det_non_ref) ] +&BEGIN_PROVIDER [ integer(bit_kind), det_cepa0_active, (N_det_non_ref) ] +&BEGIN_PROVIDER [ integer(bit_kind), det_ref_active, (N_det_ref) ] + use bitmasks + implicit none + + integer(bit_kind) :: det_noactive(N_int, 2, N_det_non_ref), active_sorb(2) + integer i, II, j, k + logical, external :: detEq + + print *, "provide cepa0" + active_sorb = (/b'001100000000', b'001100000000'/) + do i=1, N_det_non_ref + det_noactive(1,1,i) = iand(psi_non_ref(1,1,i), not(active_sorb(1))) + det_noactive(1,2,i) = iand(psi_non_ref(1,2,i), not(active_sorb(2))) + end do + + call sort_dets_ab(det_noactive, det_cepa0_idx, cepa0_shortcut, N_det_non_ref, N_int) + do i=1,N_det_ref + det_ref_active(i) = iand(psi_ref(1,1,i), active_sorb(1)) + det_ref_active(i) = det_ref_active(i) + iand(psi_ref(1,2,i), active_sorb(2)) * 2_8**32_8 + end do + + cepa0_shortcut(0) = 1 + cepa0_shortcut(1) = 1 + det_cepa0_active(1) = iand(psi_non_ref(1,1,det_cepa0_idx(1)), active_sorb(1)) + det_cepa0_active(1) = det_cepa0_active(1) + iand(psi_non_ref(1,2,det_cepa0_idx(1)), active_sorb(2)) * 2_8**32_8 + + do i=2,N_det_non_ref + det_cepa0_active(i) = iand(psi_non_ref(1,1,det_cepa0_idx(i)), active_sorb(1)) + det_cepa0_active(i) = det_cepa0_active(i) + iand(psi_non_ref(1,2,det_cepa0_idx(i)), active_sorb(2)) * 2_8**32_8 + + if(.not. detEq(det_noactive(1,1,i), det_noactive(1,1,i-1), N_int)) then + cepa0_shortcut(0) += 1 + cepa0_shortcut(cepa0_shortcut(0)) = i + end if + end do + cepa0_shortcut(0) += 1 + cepa0_shortcut(cepa0_shortcut(0)) = N_det_non_ref+1 + +END_PROVIDER + + + + BEGIN_PROVIDER [ double precision, delta_cas, (N_det_ref, N_det_ref, N_states) ] + use bitmasks + implicit none + integer :: i,j,k + double precision :: Hjk, Hki + integer i_state + + i_state = 1 + + do i=1,N_det_ref + do j=1,i + delta_cas(i,j,i_state) = 0d0 + do k=1,N_det_non_ref + call i_h_j(psi_ref(1,1,j), psi_non_ref(1,1,k),N_int,Hjk) + call i_h_j(psi_ref(1,1,i), psi_non_ref(1,1,k),N_int,Hki) + delta_cas(i,j,i_state) += Hjk * Hki * lambda_mrcc(i_state, k) + end do + delta_cas(j,i,i_state) = delta_cas(i,j,i_state) + end do + end do + + print *, "mrcepa0_cas_dressing", delta_cas(:,:,1) + END_PROVIDER + +logical function detEq(a,b,Nint) + use bitmasks + implicit none + integer, intent(in) :: Nint + integer(bit_kind), intent(in) :: a(Nint,2), b(Nint,2) + integer :: ni, i + + detEq = .false. + do i=1,2 + do ni=1,Nint + if(a(ni,i) /= b(ni,i)) return + end do + end do + detEq = .true. +end function + + + + BEGIN_PROVIDER [ double precision, delta_ij, (N_det_ref,N_det_non_ref,N_states) ] +&BEGIN_PROVIDER [ double precision, delta_ii, (N_det_ref,N_states) ] + use bitmasks + implicit none + + integer :: i_state, i, i_I, J, k, degree, degree2, m, l, deg, ni + integer :: p1,p2,h1,h2,s1,s2, p1_,p2_,h1_,h2_,s1_,s2_ + logical :: ok + double precision :: phase_iI, phase_Ik, phase_Jl, phase_IJ, phase_al, diI, hIi, hJi, delta_JI, dkI(1), HkI, ci_inv(1), dia_hla(1) + double precision :: contrib + integer, dimension(0:2,2,2) :: exc_iI, exc_Ik, exc_IJ + integer(bit_kind) :: det_tmp(N_int, 2), made_hole, made_particle + integer, allocatable :: idx_sorted_bit(:) + integer, external :: get_index_in_psi_det_sorted_bit + logical, external :: is_in_wavefunction + + integer :: II, blok + + provide det_cepa0_active + + if(N_int /= 1) then + print *, "mrcepa0 experimental N_int==1" + stop + end if + + i_state = 1 + delta_ii(:,:) = 0 + delta_ij(:,:,:) = 0 + +! do i=1,N_det_ref +! delta_ii(i,i_state) = delta_cas(i,i) +! end do + + provide mo_bielec_integrals_in_map + allocate(idx_sorted_bit(N_det)) + + idx_sorted_bit(:) = -1 + do i=1,N_det_non_ref + idx_sorted_bit(get_index_in_psi_det_sorted_bit(psi_non_ref(1,1,i), N_int)) = i + enddo + + !sd $OMP PARALLEL DO schedule(dynamic,10) default(firstprivate) shared(delta_ij, delta_ii) + do blok=1,cepa0_shortcut(0) + do i=cepa0_shortcut(blok), cepa0_shortcut(blok+1)-1 + do II=1,N_det_ref + + made_hole = iand(det_cepa0_active(i), xor(det_cepa0_active(i), det_ref_active(II))) + made_particle = iand(det_ref_active(II), xor(det_cepa0_active(i), det_ref_active(II))) + + if(popcnt(made_hole) + popcnt(made_particle) > 2) cycle + + do k=cepa0_shortcut(blok), cepa0_shortcut(blok+1)-1 + do J=1,N_det_ref + if(iand(made_hole, det_ref_active(J)) /= made_hole) cycle + if(iand(made_particle, det_ref_active(J)) /= 0) cycle + call i_h_j(psi_ref(1,1,II), psi_non_ref(1,1,det_cepa0_idx(k)),N_int,Hki) + contrib = Hki * lambda_mrcc(i_state, det_cepa0_idx(k)) * delta_cas(II,J,i_state) + delta_ij(II, det_cepa0_idx(i), i_state) += contrib + +! if(dabs(psi_ref_coef(II,i_state)).ge.5.d-5) then +! !qs$OMP CRITICAL +! delta_ii(II,i_state) -= contrib / psi_ref_coef(II, i_state) * psi_non_ref_coef(det_cepa0_idx(k),i_state) +! !dsd$OMP END CRITICAL +! endif + + end do + end do + end do + end do + end do + !qsd $OMP END PARALLEL DO + deallocate(idx_sorted_bit) +END_PROVIDER + + +subroutine set_det_bit(det, p, s) + use bitmasks + implicit none + integer(bit_kind),intent(inout) :: det(N_int, 2) + integer, intent(in) :: p, s + integer :: ni, pos + + ni = (p-1)/bit_kind_size + 1 + pos = mod(p-1, bit_kind_size) + det(ni,s) = ibset(det(ni,s), pos) +end subroutine + + + +subroutine apply_excitation(det, exc, res, ok, Nint) + use bitmasks + implicit none + + integer, intent(in) :: Nint + integer, intent(in) :: exc(0:2,2,2) + integer(bit_kind),intent(in) :: det(Nint, 2) + integer(bit_kind),intent(out) :: res(Nint, 2) + logical, intent(out) :: ok + integer :: h1,p1,h2,p2,s1,s2,degree + integer :: ii, pos + + + ok = .false. + degree = exc(0,1,1) + exc(0,1,2) + if(.not. (degree > 0 .and. degree <= 2)) then + print *, "apply ex" + STOP + endif + + call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) + res = det + + ii = (h1-1)/bit_kind_size + 1 + pos = mod(h1-1, 64)!iand(h1-1,bit_kind_size-1) ! mod 64 + if(iand(det(ii, s1), ishft(1_bit_kind, pos)) == 0_8) return + res(ii, s1) = ibclr(res(ii, s1), pos) + + ii = (p1-1)/bit_kind_size + 1 + pos = mod(p1-1, 64)!iand(p1-1,bit_kind_size-1) + if(iand(det(ii, s1), ishft(1_bit_kind, pos)) /= 0_8) return + res(ii, s1) = ibset(res(ii, s1), pos) + + if(degree == 2) then + ii = (h2-1)/bit_kind_size + 1 + pos = mod(h2-1, 64)!iand(h2-1,bit_kind_size-1) + if(iand(det(ii, s2), ishft(1_bit_kind, pos)) == 0_8) return + res(ii, s2) = ibclr(res(ii, s2), pos) + + ii = (p2-1)/bit_kind_size + 1 + pos = mod(p2-1, 64)!iand(p2-1,bit_kind_size-1) + if(iand(det(ii, s2), ishft(1_bit_kind, pos)) /= 0_8) return + res(ii, s2) = ibset(res(ii, s2), pos) + endif + + ok = .true. +end subroutine + + + + + + + + diff --git a/plugins/mrcepa0/mrcepa0.irp.f b/plugins/mrcepa0/mrcepa0.irp.f new file mode 100644 index 00000000..b9eb2fc5 --- /dev/null +++ b/plugins/mrcepa0/mrcepa0.irp.f @@ -0,0 +1,24 @@ +program mrcepa0 + implicit none + if (.not.read_wf) then + print *, 'read_wf has to be true.' + stop 1 + endif + call print_cas_coefs + call run_mrcepa0 +end + +subroutine print_cas_coefs + implicit none + + integer :: i,j + print *, 'CAS' + print *, '===' + do i=1,N_det_cas + print *, psi_cas_coef(i,:) + call debug_det(psi_cas(1,1,i),N_int) + enddo + + call write_double(6,ci_energy(1),"Initial CI energy") +end + diff --git a/plugins/mrcepa0/mrcepa0_general.irp.f b/plugins/mrcepa0/mrcepa0_general.irp.f new file mode 100644 index 00000000..3f892887 --- /dev/null +++ b/plugins/mrcepa0/mrcepa0_general.irp.f @@ -0,0 +1,96 @@ +subroutine run_mrcepa0 + implicit none + call set_generators_bitmasks_as_holes_and_particles + call mrcepa0_iterations +end + +subroutine mrcepa0_iterations + implicit none + + integer :: i,j + + double precision :: E_new, E_old, delta_e + integer :: iteration,i_oscillations + double precision :: E_past(4) + E_new = 0.d0 + delta_E = 1.d0 + iteration = 0 + j = 1 + i_oscillations = 0 + do while (delta_E > 1.d-7) + iteration += 1 + print *, '===========================' + print *, 'MRCEPA0 Iteration', iteration + print *, '===========================' + print *, '' + E_old = sum(ci_energy_dressed) + call write_double(6,ci_energy_dressed(1),"MRCEPA0 energy") + call diagonalize_ci_dressed + E_new = sum(ci_energy_dressed) + delta_E = dabs(E_new - E_old) + E_past(j) = E_new + j +=1 + if(j>4)then + j=1 + endif + if(iteration > 4) then + if(delta_E > 1.d-10)then + if(dabs(E_past(1) - E_past(3)) .le. delta_E .and. dabs(E_past(2) - E_past(4)).le. delta_E)then + print*,'OSCILLATIONS !!!' +! oscillations = .True. + i_oscillations +=1 + lambda_mrcc_tmp = lambda_mrcc + endif + endif + endif + call save_wavefunction +! if (i_oscillations > 5) then +! exit +! endif + if (iteration > 100) then + exit + endif + print*,'------------' + print*,'VECTOR' + do i = 1, N_det_ref + print*,'' + print*,'psi_ref_coef(i,1) = ',psi_ref_coef(i,1) + print*,'delta_ii(i,1) = ',delta_ii(i,1) + enddo + print*,'------------' + enddo + call write_double(6,ci_energy_dressed(1),"Final MRCEPA0 energy") + call ezfio_set_mrcc_cassd_energy(ci_energy_dressed(1)) + call save_wavefunction + +end + +subroutine set_generators_bitmasks_as_holes_and_particles + implicit none + integer :: i,k + do k = 1, N_generators_bitmask + do i = 1, N_int + ! Pure single part + generators_bitmask(i,1,1,k) = holes_operators(i,1) ! holes for pure single exc alpha + generators_bitmask(i,1,2,k) = particles_operators(i,1) ! particles for pure single exc alpha + generators_bitmask(i,2,1,k) = holes_operators(i,2) ! holes for pure single exc beta + generators_bitmask(i,2,2,k) = particles_operators(i,2) ! particles for pure single exc beta + + ! Double excitation + generators_bitmask(i,1,3,k) = holes_operators(i,1) ! holes for first single exc alpha + generators_bitmask(i,1,4,k) = particles_operators(i,1) ! particles for first single exc alpha + generators_bitmask(i,2,3,k) = holes_operators(i,2) ! holes for first single exc beta + generators_bitmask(i,2,4,k) = particles_operators(i,2) ! particles for first single exc beta + + generators_bitmask(i,1,5,k) = holes_operators(i,1) ! holes for second single exc alpha + generators_bitmask(i,1,6,k) = particles_operators(i,1) ! particles for second single exc alpha + generators_bitmask(i,2,5,k) = holes_operators(i,2) ! holes for second single exc beta + generators_bitmask(i,2,6,k) = particles_operators(i,2) ! particles for second single exc beta + + enddo + enddo + touch generators_bitmask + + + +end