LCPQ/quantum_package
10
0
Fork 0
mirror of https://github.com/LCPQ/quantum_package synced 2024-06-26 07:02:14 +02:00
Code Issues Releases Wiki Activity

different binaries for mrcepa0/mrsc2/mrsc2sub + corrected reversed index in mrcc_utils

Browse Source
This commit is contained in:
Yann Garniron 2016-04-01 12:00:03 +02:00
parent 993a034f07
commit 48cb3b3ddc
6 changed files with 275 additions and 201 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
plugins/MRCC_Utils/davidson.irp.f
View File

@ -486,8 +486,8 @@ subroutine H_u_0_mrcc(v_0,u_0,H_jj,n,keys_tmp,Nint,istate)
i = idx_ref(ii)
do jj = 1, n_det_non_ref
j = idx_non_ref(jj)
vt (i) = vt (i) + delta_ij(ii,jj,istate)*u_0(j)
vt (j) = vt (j) + delta_ij(ii,jj,istate)*u_0(i)
vt (i) = vt (i) + delta_ij(istate,jj,ii)*u_0(j)
vt (j) = vt (j) + delta_ij(istate,jj,ii)*u_0(i)
enddo
enddo
!$OMP END DO

4
plugins/MRCC_Utils/mrcc_general.irp.f
View File

@ -35,8 +35,8 @@ subroutine mrcc_iterations
! lambda = min(1.d0, lambda * 1.1d0)
! endif
! print *, 'energy lambda ', lambda
E_past(j) = E_new
j +=1
! E_past(j) = E_new
! j +=1
call save_wavefunction
if (iteration > 200) then
exit

1
plugins/MRCC_Utils/mrcc_utils.irp.f
View File

@ -191,4 +191,3 @@ subroutine diagonalize_CI_dressed(lambda)
SOFT_TOUCH psi_coef
end

426
plugins/mrcepa0/dressing.irp.f
View File

@ -2,19 +2,47 @@ use bitmasks
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_ij, (N_states,N_det_non_ref,N_det_ref) ]
&BEGIN_PROVIDER [ double precision, delta_ii, (N_states, N_det_ref) ]
use bitmasks
implicit none
integer :: i, j, i_state
! delta_ij(:,:,:) = delta_ij_old(:,:,:)
! delta_ii(:,:) = delta_ii_old(:,:)
delta_ij(:,:,:) = delta_mrcepa0_ij(:,:,:)! - delta_sub_ij(:,:,:)
delta_ii(:,:)= delta_mrcepa0_ii(:,:)! - delta_sub_ii(:,:)
!mrmode : 1=mrcepa0, 2=mrsc2 add, 3=mrsc2 sub
do i_state = 1, N_states
if(mrmode == 3) then
do i = 1, N_det_ref
delta_ii(i_state,i)= delta_mrcepa0_ii(i,i_state) - delta_sub_ii(i,i_state)
do j = 1, N_det_non_ref
delta_ij(i_state,j,i) = delta_mrcepa0_ij(i,j,i_state) - delta_sub_ij(i,j,i_state)
end do
end do
else if(mrmode == 2) then
do i = 1, N_det_ref
delta_ii(i_state,i)= delta_ii_old(i,i_state)
do j = 1, N_det_non_ref
delta_ij(i_state,j,i) = delta_ij_old(i,j,i_state)
end do
end do
else if(mrmode == 1) then
do i = 1, N_det_ref
delta_ii(i_state,i)= delta_mrcepa0_ii(i,i_state)
do j = 1, N_det_non_ref
delta_ij(i_state,j,i) = delta_mrcepa0_ij(i,j,i_state)
end do
end do
else
stop "invalid mrmode"
end if
end do
END_PROVIDER
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) ]
@ -27,7 +55,6 @@ END_PROVIDER
integer i, II, j, k
logical, external :: detEq
print *, "provide cepa0"
active_sorb(:) = 0_8
nonactive_sorb(:) = not(0_8)
@ -49,10 +76,6 @@ END_PROVIDER
call sort_dets_ab(det_noactive, det_cepa0_idx, cepa0_shortcut, N_det_non_ref, N_int)
! do i=1, N_det_non_ref
! print "(B30,B30)", det_noactive(1,1,i), det_noactive(1,2,i)
! end do
! stop
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
@ -83,17 +106,25 @@ END_PROVIDER
use bitmasks
implicit none
integer :: i,j,k
double precision :: Hjk, Hki, Hij, mat(2,2)
integer i_state
double precision :: Hjk, Hki, Hij
integer i_state, degree
provide lambda_mrcc
i_state = 1
do i=1,N_det_ref
do j=1,i
call get_excitation_degree(psi_ref(1,1,i), psi_ref(1,1,j), degree, N_int)
if(degree /= 2 .and. degree /= 0) cycle
delta_cas(i,j,i_state) = 0d0
do k=1,N_det_non_ref
! call get_excitation_degree(psi_ref(1,1,j), psi_non_ref(1,1,k), degree, N_int)
! if(degree /= 2) cycle
! call get_excitation_degree(psi_ref(1,1,i), psi_non_ref(1,1,k), degree, N_int)
! if(degree /= 2) cycle
call i_h_j(psi_ref(1,1,j), psi_non_ref(1,1,k),N_int,Hjk)
call i_h_j(psi_non_ref(1,1,k),psi_ref(1,1,i), 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)
@ -101,6 +132,9 @@ END_PROVIDER
end do
END_PROVIDER
!-199.0906497310625
!-199.0913388716010
logical function detEq(a,b,Nint)
use bitmasks
implicit none
@ -120,8 +154,8 @@ end function
BEGIN_PROVIDER [ double precision, delta_mrcepa0_ij, (N_det_ref,N_det_non_ref,N_states) ]
&BEGIN_PROVIDER [ double precision, delta_mrcepa0_ii, (N_det_ref,N_states) ]
BEGIN_PROVIDER [ double precision, delta_mrcepa0_ij_old, (N_det_ref,N_det_non_ref,N_states) ]
&BEGIN_PROVIDER [ double precision, delta_mrcepa0_ii_old, (N_det_ref,N_states) ]
use bitmasks
implicit none
@ -135,11 +169,14 @@ end function
integer, allocatable :: idx_sorted_bit(:)
integer, external :: get_index_in_psi_det_sorted_bit
logical, external :: is_in_wavefunction
double precision, allocatable :: hab(:,:)
integer :: II, blok
provide det_cepa0_active delta_cas lambda_mrcc
if(N_int /= 1) then
print *, "mrcepa0 experimental N_int==1"
stop
@ -148,7 +185,8 @@ end function
i_state = 1
delta_mrcepa0_ii(:,:) = 0d0
delta_mrcepa0_ij(:,:,:) = 0d0
!allocate(hab(N_det_non_ref, N_det_non_ref))
!hab(:,:) = 0d0
! do i=1,N_det_ref
! delta_ii(i,i_state) = delta_cas(i,i,i_state)
! end do
@ -168,7 +206,7 @@ end function
made_hole = iand(det_ref_active(II), xor(det_cepa0_active(i), det_ref_active(II)))
made_particle = iand(det_cepa0_active(i), xor(det_cepa0_active(i), det_ref_active(II)))
call get_excitation_degree(psi_ref(1,1,II),psi_non_ref(1,1,det_cepa0_idx(i)),degree,N_int)
if (degree > 2 .or. popcnt(made_hole) * popcnt(made_particle) /= degree*2) cycle
if (degree > 2 .or. popcnt(made_hole) + popcnt(made_particle) == 7650) cycle
do k=cepa0_shortcut(blok), cepa0_shortcut(blok+1)-1
if(iand(not(active_sorb(1)), xor(psi_non_ref(1,1,det_cepa0_idx(k)), psi_non_ref(1,1,det_cepa0_idx(i)))) /= 0) stop "STOOOP"
@ -185,19 +223,24 @@ end function
!!!!!
call get_excitation_degree(psi_ref(1,1,J),psi_non_ref(1,1,det_cepa0_idx(k)),degree,N_int)
if(degree > 2) stop "BBBB"
!!!!!!!!!
! if(i/=k .and. popcnt(made_hole) /= popcnt(made_particle)) then
! print *, "=================", made_hole, made_particle
! call debug_det(psi_ref(1,1,II),N_int)
! call debug_det(psi_non_ref(1,1,det_cepa0_idx(i)),N_int)
! call debug_det(psi_ref(1,1,J),N_int)
! call debug_det(psi_non_ref(1,1,det_cepa0_idx(k)),N_int)
! print *, "================="
! end if
call i_h_j(psi_non_ref(1,1,det_cepa0_idx(k)),psi_ref(1,1,II),N_int,Hki)
!call i_h_j(psi_non_ref(1,1,det_cepa0_idx(k)),psi_ref(1,1,J),N_int,Hki)
contrib = Hki * lambda_mrcc(i_state, det_cepa0_idx(k)) * delta_cas(II,J,i_state)
!contrib = Hki * lambda_mrcc(i_state, det_cepa0_idx(k)) * delta_cas(II,J,i_state)! * psi_ref_coef(II, I_state)*psi_ref_coef(J, I_state)/(psi_ref_coef(1, I_state)**2 + psi_ref_coef(2, I_state)**2)
contrib = Hki * lambda_mrcc(i_state, det_cepa0_idx(k)) * delta_cas(II,J,i_state)! * psi_ref_coef(II, I_state)*psi_ref_coef(J, I_state)/(psi_ref_coef(1, I_state)**2 + psi_ref_coef(2, I_state)**2)
! (psi_ref_coef(II, I_state) * psi_ref_coef(J, I_state)) / (psi_ref_coef(1, I_state)**2 + psi_ref_coef(2, I_state)**2)
! if(hab(det_cepa0_idx(k), det_cepa0_idx(i)) /= 0) then
! print *, "HAB ", contrib, hab(det_cepa0_idx(k), det_cepa0_idx(i))
! !contrib = 0d0
! !stop
! else
! hab(det_cepa0_idx(k), det_cepa0_idx(i)) = contrib
! hab(det_cepa0_idx(i), det_cepa0_idx(k)) = contrib
! end if
delta_mrcepa0_ij(II, det_cepa0_idx(i), i_state) += contrib
!
if(dabs(psi_ref_coef(II,i_state)).ge.5.d-5) then
@ -222,6 +265,152 @@ END_PROVIDER
BEGIN_PROVIDER [ double precision, delta_mrcepa0_ij, (N_det_ref,N_det_non_ref,N_states) ]
&BEGIN_PROVIDER [ double precision, delta_mrcepa0_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, HIIi, HJk
integer, dimension(0:2,2,2) :: exc_iI, exc_Ik, exc_IJ
integer(bit_kind) :: det_tmp(N_int, 2), made_hole, made_particle, myActive
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 delta_cas lambda_mrcc
if(N_int /= 1) then
print *, "mrcepa0 experimental N_int==1"
stop
end if
i_state = 1
delta_mrcepa0_ii(:,:) = 0d0
delta_mrcepa0_ij(:,:,:) = 0d0
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
!- qsd $OMP PARALLEL DO schedule(dynamic,10) default(firstprivate) shared(delta_mrcepa0_ij, delta_mrcepa0_ii)
do blok=1,cepa0_shortcut(0)
do i=cepa0_shortcut(blok), cepa0_shortcut(blok+1)-1
do II=1,N_det_ref
call get_excitation_degree(psi_ref(1,1,II),psi_non_ref(1,1,det_cepa0_idx(i)),degree,N_int)
if (degree > 2 ) cycle
made_hole = iand(det_ref_active(II), xor(det_cepa0_active(i), det_ref_active(II)))
made_particle = iand(det_cepa0_active(i), xor(det_cepa0_active(i), det_ref_active(II)))
do k=cepa0_shortcut(blok), cepa0_shortcut(blok+1)-1
!if(i==k) cycle
if(iand(not(active_sorb(1)), xor(psi_non_ref(1,1,det_cepa0_idx(k)), psi_non_ref(1,1,det_cepa0_idx(i)))) /= 0) stop "STOOOP"
if(iand(made_hole, det_cepa0_active(k)) /= 0) cycle
if(iand(made_particle, det_cepa0_active(k)) /= made_particle) cycle
myActive = xor(det_cepa0_active(k), made_hole)
myActive = xor(myActive, made_particle)
if(i==k .and. myActive /= det_ref_active(II)) stop "AAAA"
do J=1,N_det_ref
if(det_ref_active(J) /= myActive) cycle
call i_h_j(psi_non_ref(1,1,det_cepa0_idx(k)),psi_ref(1,1,J),N_int,HJk)
contrib = delta_cas(II, J, i_state) * HJk * lambda_mrcc(i_state, det_cepa0_idx(k))
delta_mrcepa0_ij(J, det_cepa0_idx(i), i_state) += contrib
if(dabs(psi_ref_coef(J,i_state)).ge.5.d-5) then
delta_mrcepa0_ii(J,i_state) -= contrib / psi_ref_coef(J, i_state) * psi_non_ref_coef(det_cepa0_idx(k),i_state)
end if
end do
end do
end do
end do
end do
deallocate(idx_sorted_bit)
END_PROVIDER
BEGIN_PROVIDER [ double precision, delta_mrcepa0_ij_exp, (N_det_ref,N_det_non_ref,N_states) ]
&BEGIN_PROVIDER [ double precision, delta_mrcepa0_ii_exp, (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, HIIi, HJk
integer, dimension(0:2,2,2) :: exc_iI, exc_Ik, exc_IJ
integer(bit_kind) :: det_tmp(N_int, 2), made_hole, made_particle, myActive
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 delta_cas lambda_mrcc
print *, "mrcepa0 experimental"
if(N_int /= 1) then
print *, "mrcepa0 experimental N_int==1"
stop
end if
i_state = 1
delta_mrcepa0_ii_exp(:,:) = 0d0
delta_mrcepa0_ij_exp(:,:,:) = 0d0
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
!- qsd $OMP PARALLEL DO schedule(dynamic,10) default(firstprivate) shared(delta_mrcepa0_ii_exp, delta_mrcepa0_ij_exp)
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_ref_active(II), xor(det_cepa0_active(i), det_ref_active(II)))
made_particle = iand(det_cepa0_active(i), xor(det_cepa0_active(i), det_ref_active(II)))
do k=cepa0_shortcut(blok), cepa0_shortcut(blok+1)-1
do J=1,N_det_ref
if(made_hole /= iand(det_ref_active(J), xor(det_cepa0_active(k), det_ref_active(J)))) cycle
if(made_particle /= iand(det_cepa0_active(k), xor(det_cepa0_active(k), det_ref_active(J)))) cycle
call i_h_j(psi_non_ref(1,1,det_cepa0_idx(k)),psi_ref(1,1,J),N_int,HJk)
contrib = delta_cas(II, J, i_state) * HJk * lambda_mrcc(i_state, det_cepa0_idx(k))
delta_mrcepa0_ij_exp(J, det_cepa0_idx(i), i_state) += contrib
if(dabs(psi_ref_coef(J,i_state)).ge.5.d-5) then
delta_mrcepa0_ii_exp(J,i_state) -= contrib / psi_ref_coef(J, i_state) * psi_non_ref_coef(det_cepa0_idx(k),i_state)
end if
end do
end do
end do
end do
end do
deallocate(idx_sorted_bit)
END_PROVIDER
BEGIN_PROVIDER [ double precision, delta_sub_ij, (N_det_ref,N_det_non_ref,N_states) ]
&BEGIN_PROVIDER [ double precision, delta_sub_ii, (N_det_ref, N_states) ]
use bitmasks
@ -260,7 +449,7 @@ END_PROVIDER
!$OMP PARALLEL DO schedule(dynamic,10) default(firstprivate) shared(delta_sub_ij, delta_sub_ii)
do i=1,N_det_non_ref
if(mod(i,1000) == 0) print "(A,I3,A)", "♫ sloubi", i/1000, " ♪"
if(mod(i,1000) == 0) print "(A,I3,A)", "♫ sloubi", i/1000, " ♪ (sub)"
do J=1,N_det_ref
call get_excitation(psi_ref(1,1,J),psi_non_ref(1,1,i),exc_Ji,degree,phase_Ji,N_int)
if(degree == -1) cycle
@ -279,6 +468,8 @@ END_PROVIDER
call i_h_j(psi_ref(1,1,II), det_tmp, N_int, HIl)
do k=1,N_det_non_ref
! if(lambda_mrcc(i_state, k) == 0d0) cycle
if(lambda_mrcc(i_state, k) == 0d0) cycle
call get_excitation(psi_ref(1,1,II),psi_non_ref(1,1,k),exc_Ik,degree2,phase_Ik,N_int)
det_tmp(:,:) = 0_bit_kind
@ -295,42 +486,16 @@ END_PROVIDER
if(ok) cycle
! call decode_exc(exc_Ii,degree,h1_,p1_,h2_,p2_,s1_,s2_)
! call decode_exc(exc_Ik,degree2,h1,p1,h2,p2,s1,s2)
!
!
! det_tmp(:,:) = 0_bit_kind
! call set_det_bit(det_tmp, p1, s1)
! call set_det_bit(det_tmp, h1, s1)
! call set_det_bit(det_tmp, p1_, s1_)
! call set_det_bit(det_tmp, h1_, s1_)
! if(degree == 2) then
! call set_det_bit(det_tmp, p2_, s2_)
! call set_det_bit(det_tmp, h2_, s2_)
! end if
! if(degree2 == 2) then
! call set_det_bit(det_tmp, p2, s2)
! call set_det_bit(det_tmp, h2, s2)
! end if
! deg = 0
! do ni = 1, N_int
! deg += popcnt(det_tmp(ni,1)) + popcnt(det_tmp(ni,2))
! end do
! if(deg == 2*degree2 + 2*degree) cycle
! if(degree == -1) cycle
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,II), psi_non_ref(1,1,k), N_int, HIk)
if(HJk == 0) cycle
!assert HIk == 0
!delta_IJk = HJk * HIk * lambda_mrcc(i_state, k)
delta_IJk = HJk * HIk * lambda_mrcc(i_state, k)
call apply_excitation(psi_non_ref(1,1,i),exc_Ik,det_tmp,ok,N_int)
if(ok) cycle
contrib = delta_IJk * HIl * lambda_mrcc(i_state, l)
! contrib = delta_IJk * HIl * lambda_mrcc(i_state, l)
contrib = delta_IJk * HIl * lambda_mrcc(i_state,l)
!$OMP CRITICAL
delta_sub_ij(II, i, i_state) += contrib
if(dabs(psi_ref_coef(II,i_state)).ge.5.d-5) then
@ -376,6 +541,7 @@ implicit none
delta_ii_old(:,:) = 0
delta_ij_old(:,:,:) = 0
i_state = 1
provide mo_bielec_integrals_in_map
allocate(idx_sorted_bit(N_det))
@ -387,7 +553,8 @@ implicit none
!$OMP PARALLEL DO schedule(dynamic,10) default(firstprivate) shared(delta_ij_old, delta_ii_old)
do i = 1 , N_det_non_ref
if(mod(i,1000) == 0) print *, i, N_det_non_ref
if(lambda_mrcc(i_state, i) == 0d0) cycle
if(mod(i,1000) == 0) print "(A,I3,A)", "♫ sloubi", i/1000, " ♪ (old)"
do i_I = 1 , N_det_ref
call get_excitation(psi_ref(1,1,i_I),psi_non_ref(1,1,i),exc_iI,degree2,phase_iI,N_int)
if(degree2 == -1) cycle
@ -396,12 +563,14 @@ implicit none
call decode_exc(exc_iI,degree2,h1,p1,h2,p2,s1,s2)
call i_h_j(psi_non_ref(1,1,i), psi_ref(1,1,i_I),N_int,hIi)
diI = hIi * lambda_mrcc(i_state,i)
diI = hIi * lambda_mrcc(i_state, i)
do J = 1 , N_det_ref !!!
call get_excitation(psi_ref(1,1,i_I),psi_ref(1,1,J),exc_IJ,degree,phase_IJ,N_int)
call i_h_j(psi_non_ref(1,1,i), psi_ref(1,1,J),N_int,hJi)
delta_JI = hJi * diI
do k = 1 , N_det_non_ref
if(lambda_mrcc(i_state, k) == 0d0) cycle
call get_excitation(psi_ref(1,1,i_I),psi_non_ref(1,1,k),exc_Ik,degree,phase_Ik,N_int)
if(degree == -1) cycle
@ -411,11 +580,6 @@ implicit none
det_tmp(:,:) = 0_bit_kind
det_tmp2(:,:) = 0_bit_kind
!!!!!!!!!!!!!!!
det_tmp(1,1) = iand(xor(HF_bitmask(1,1), psi_non_ref(1,1,k)), not(active_sorb(1)))
det_tmp(1,2) = iand(xor(HF_bitmask(1,1), psi_non_ref(1,1,i)), not(active_sorb(1)))
@ -427,60 +591,6 @@ implicit none
if(.not. ok) cycle
!if(ok) cycle
!!!!!!!!!!!!!!
! call set_det_bit(det_tmp, p1, s1)
!
! call set_det_bit(det_tmp, p1_, s1_)
!
! if(degree == 2) then
! call set_det_bit(det_tmp, p2_, s2_)
!
! end if
! if(degree2 == 2) then
! call set_det_bit(det_tmp, p2, s2)
! end if
!
! x(:) = 0
! do ni=1,N_int
! x(1) += popcnt(iand(det_tmp(ni, 1), cas_bitmask(ni, 1, 1)))
! x(2) += popcnt(iand(det_tmp(ni, 2), cas_bitmask(ni, 2, 1)))
! end do
!
!
! !det_tmp(:,:) = 0_bit_kind
!
! call set_det_bit(det_tmp, h1, s1)
! call set_det_bit(det_tmp, h1_, s1_)
! if(degree == 2) then
! call set_det_bit(det_tmp, h2_, s2_)
! end if
! if(degree2 == 2) then
! call set_det_bit(det_tmp, h2, s2)
! end if
!
! y(1) = -x(1)
! y(2) = -x(2)
! do ni=1,N_int
! y(1) += popcnt(iand(det_tmp(ni, 1), cas_bitmask(ni, 1, 1)))
! y(2) += popcnt(iand(det_tmp(ni, 2), cas_bitmask(ni, 2, 1)))
! end do
!
! ! print *, x, y
!
! if(x(1) * y(1) /= 0) cycle
! if(x(2) * y(2) /= 0) cycle
!
!
!
! deg = 0
! do ni = 1, N_int
! deg += popcnt(det_tmp(ni,1)) + popcnt(det_tmp(ni,2))
! end do
! if(deg /= 2*degree2 + 2*degree) cycle
call apply_excitation(psi_non_ref(1,1,i),exc_Ik,det_tmp,ok,N_int)
@ -500,17 +610,22 @@ implicit none
l = get_index_in_psi_det_sorted_bit(det_tmp, N_int)
if(l == 0) cycle
l = idx_sorted_bit(get_index_in_psi_det_sorted_bit(det_tmp, N_int))
if(l ==-1) cycle
call i_h_j(psi_non_ref(1,1,k), psi_ref(1,1,i_I),N_int,HkI)
dkI(i_state) = HkI * lambda_mrcc(i_state,k) * phase_Jl * phase_Ik
!dkI(i_state) = HkI * lambda_mrcc(i_state,k) * phase_Jl * phase_Ik * Xref(I_i)
dkI(i_state) = HkI * lambda_mrcc(i_state, k) * phase_Jl * phase_Ik
!!!!
call get_excitation(psi_ref(1,1,J),psi_non_ref(1,1,k),exc_Ik,degree,phase_Ik,N_int)
if(degree /= 2 .and. degree /= 0) cycle
!!!!!!
!$OMP CRITICAL
contrib = dkI(i_state) * delta_JI
!erro += abs(dkI(i_state) - psi_non_ref_coef(k,i_state) / psi_ref_coef(1,i_state))
delta_ij_old(i_I,l,i_state) += contrib
if(dabs(psi_ref_coef(i_I,i_state)).ge.5.d-5) then
delta_ii_old(i_I,i_state) -= contrib * ci_inv(i_state) * psi_non_ref_coef(k,i_state)
@ -521,61 +636,24 @@ implicit none
enddo
enddo
!$OMP END PARALLEL DO
! double precision :: error, acc
! integer :: II
! error = 0d0
! do i=1, N_det_non_ref
! acc = 0d0
! do II=1, N_det_ref
! call i_h_j(psi_ref(1,1,II), psi_non_ref(1,1,i), N_int, HIi)
! acc += HIi * lambda_mrcc(i_state, i) * Xref(II) * psi_ref_coef(II, i_state)
! end do
! error += (psi_non_ref_coef(i, i_state) - acc)**2
! end do
! print *, "QUALITY ", error
deallocate(idx_sorted_bit)
END_PROVIDER
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

2
plugins/mrcepa0/mrcepa0.irp.f
View File

@ -1,5 +1,7 @@
program mrcepa0
implicit none
!mrmode : 1=mrcepa0, 2=mrsc2 add, 3=mrsc2 sub
mrmode = 1
if (.not.read_wf) then
print *, 'read_wf has to be true.'
stop 1

39
plugins/mrcepa0/mrcepa0_general.irp.f
View File

@ -4,6 +4,10 @@ subroutine run_mrcepa0
call mrcepa0_iterations
end
BEGIN_PROVIDER [ integer, mrmode ]
END_PROVIDER
subroutine mrcepa0_iterations
implicit none
@ -11,12 +15,13 @@ subroutine mrcepa0_iterations
double precision :: E_new, E_old, delta_e
integer :: iteration,i_oscillations
double precision :: E_past(4)
double precision :: E_past(4), lambda
E_new = 0.d0
delta_E = 1.d0
iteration = 0
j = 1
i_oscillations = 0
lambda = 1.d0
do while (delta_E > 1.d-7)
iteration += 1
print *, '==========================='
@ -25,29 +30,19 @@ subroutine mrcepa0_iterations
print *, ''
E_old = sum(ci_energy_dressed)
call write_double(6,ci_energy_dressed(1),"MRCEPA0 energy")
call diagonalize_ci_dressed
call diagonalize_ci_dressed(lambda)
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
! if (E_new > E_old) then
! lambda = lambda * 0.7d0
! else
! lambda = min(1.d0, lambda * 1.1d0)
! endif
! print *, 'energy lambda ', lambda
! E_past(j) = E_new
! j +=1
call save_wavefunction
! if (i_oscillations > 5) then
! exit
! endif
if (iteration > 100) then
if (iteration > 200) then
exit
endif
print*,'------------'
@ -55,7 +50,7 @@ subroutine mrcepa0_iterations
do i = 1, N_det_ref
print*,''
print*,'psi_ref_coef(i,1) = ',psi_ref_coef(i,1)
print*,'delta_ii(i,1) = ',delta_cas(i,i,1)
print*,'delta_ii(i,1) = ',delta_ii(i,1)
enddo
print*,'------------'
enddo