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

Corrected PT2 for MRCC

Browse Source
This commit is contained in:
Anthony Scemama 2016-09-17 23:33:06 +02:00
parent bd9d322068
commit cc947afc89
5 changed files with 142 additions and 59 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
plugins/MRCC_Utils/H_apply.irp.f
View File

@ -32,7 +32,7 @@ s.unset_openmp()
print s
s = H_apply_zmq("mrcepa_PT2")
s.energy = "psi_ref_energy_diagonalized"
s.energy = "psi_energy"
s.set_perturbation("epstein_nesbet_2x2")
s.unset_openmp()
print s

2
plugins/MRCC_Utils/NEEDED_CHILDREN_MODULES
View File

@ -1 +1 @@
Perturbation Selectors_full Generators_full Psiref_Utils Psiref_CAS
Perturbation Selectors_full Generators_full Psiref_Utils Psiref_CAS

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

@ -7,7 +7,7 @@ END_PROVIDER
BEGIN_PROVIDER [ double precision, lambda_mrcc, (N_states, N_det_non_ref) ]
&BEGIN_PROVIDER [ integer, lambda_mrcc_pt2, (0:psi_det_size) ]
&BEGIN_PROVIDER [ integer, lambda_mrcc_pt3, (0:psi_det_size) ]
&BEGIN_PROVIDER [ integer, lambda_mrcc_kept, (0:psi_det_size) ]
implicit none
BEGIN_DOC
! cm/<Psi_0|H|D_m> or perturbative 1/Delta_E(m)
@ -23,7 +23,7 @@ END_PROVIDER
N_lambda_mrcc_pt2 = 0
N_lambda_mrcc_pt3 = 0
lambda_mrcc_pt2(0) = 0
lambda_mrcc_pt3(0) = 0
lambda_mrcc_kept(0) = 0
do i=1,N_det_non_ref
call i_h_psi(psi_non_ref(1,1,i), psi_ref, psi_ref_coef, N_int, N_det_ref,&
@ -36,6 +36,7 @@ END_PROVIDER
lambda_mrcc(k,i) = min(-1.d-32,psi_non_ref_coef(i,k)/ihpsi_current(k) )
lambda_pert = 1.d0 / (psi_ref_energy_diagonalized(k)-hii)
if (lambda_pert / lambda_mrcc(k,i) < 0.5d0) then
! Ignore lamdba
i_pert_count += 1
lambda_mrcc(k,i) = 0.d0
if (lambda_mrcc_pt2(N_lambda_mrcc_pt2) /= i) then
@ -43,15 +44,16 @@ END_PROVIDER
lambda_mrcc_pt2(N_lambda_mrcc_pt2) = i
endif
else
if (lambda_mrcc_pt3(N_lambda_mrcc_pt3) /= i) then
! Keep lamdba
if (lambda_mrcc_kept(N_lambda_mrcc_pt3) /= i) then
N_lambda_mrcc_pt3 += 1
lambda_mrcc_pt3(N_lambda_mrcc_pt3) = i
lambda_mrcc_kept(N_lambda_mrcc_pt3) = i
endif
endif
enddo
enddo
lambda_mrcc_pt2(0) = N_lambda_mrcc_pt2
lambda_mrcc_pt3(0) = N_lambda_mrcc_pt3
lambda_mrcc_kept(0) = N_lambda_mrcc_pt3
print*,'N_det_non_ref = ',N_det_non_ref
print*,'psi_coef_ref_ratio = ',psi_ref_coef(2,1)/psi_ref_coef(1,1)
print*,'lambda max = ',maxval(dabs(lambda_mrcc))

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

@ -59,49 +59,89 @@ subroutine run(N_st,energy)
end
subroutine run_pt2(N_st,energy)
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
subroutine run_pt2_old(N_st,energy)
implicit none
integer :: i,j,k
double precision, allocatable :: pt2(:), norm_pert(:), H_pert_diag(:)
integer, intent(in) :: N_st
double precision, intent(in) :: energy(N_st)
double precision :: pt3(N_st)
allocate (pt2(N_st), norm_pert(N_st),H_pert_diag(N_st))
pt2 = 0.d0
pt3 = 0d0
double precision :: pt2_redundant(N_st), pt2(N_st)
double precision :: norm_pert(N_st),H_pert_diag(N_st)
pt2_redundant = 0.d0
pt2 = 0d0
!if(lambda_mrcc_pt2(0) == 0) return
print*,'Last iteration only to compute the PT2'
threshold_selectors = 1.d0
threshold_generators = 0.999d0
N_det_generators = lambda_mrcc_pt3(0)
N_det_selectors = lambda_mrcc_pt3(0)
! psi_det_generators(:,:,:N_det_ref) = psi_ref(:,:,:N_det_ref)
! psi_selectors(:,:,:N_det_ref) = psi_ref(:,:,:N_det_ref)
! psi_coef_generators(:N_det_ref,:) = psi_ref_coef(:N_det_ref,:)
! psi_selectors_coef(:N_det_ref,:) = psi_ref_coef(:N_det_ref,:)
do i=1,N_det_generators
j = lambda_mrcc_pt3(i)
do k=1,N_int
psi_det_generators(k,1,i) = psi_non_ref(k,1,j)
psi_det_generators(k,2,i) = psi_non_ref(k,2,j)
psi_selectors(k,1,i) = psi_non_ref(k,1,j)
psi_selectors(k,2,i) = psi_non_ref(k,2,j)
print * ,'Computing the redundant PT2 contribution'
if (mrmode == 1) then
N_det_generators = lambda_mrcc_kept(0)
N_det_selectors = lambda_mrcc_kept(0)
do i=1,N_det_generators
j = lambda_mrcc_kept(i)
do k=1,N_int
psi_det_generators(k,1,i) = psi_non_ref(k,1,j)
psi_det_generators(k,2,i) = psi_non_ref(k,2,j)
psi_selectors(k,1,i) = psi_non_ref(k,1,j)
psi_selectors(k,2,i) = psi_non_ref(k,2,j)
enddo
do k=1,N_st
psi_coef_generators(i,k) = psi_non_ref_coef(j,k)
psi_selectors_coef(i,k) = psi_non_ref_coef(j,k)
enddo
enddo
do k=1,N_st
psi_coef_generators(i,k) = psi_non_ref_coef(j,k)
psi_selectors_coef(i,k) = psi_non_ref_coef(j,k)
else
N_det_generators = N_det_non_ref
N_det_selectors = N_det_non_ref
do i=1,N_det_generators
j = i
do k=1,N_int
psi_det_generators(k,1,i) = psi_non_ref(k,1,j)
psi_det_generators(k,2,i) = psi_non_ref(k,2,j)
psi_selectors(k,1,i) = psi_non_ref(k,1,j)
psi_selectors(k,2,i) = psi_non_ref(k,2,j)
enddo
do k=1,N_st
psi_coef_generators(i,k) = psi_non_ref_coef(j,k)
psi_selectors_coef(i,k) = psi_non_ref_coef(j,k)
enddo
enddo
enddo
endif
SOFT_TOUCH N_det_selectors psi_selectors_coef psi_selectors N_det_generators psi_det_generators psi_coef_generators ci_eigenvectors_dressed ci_eigenvectors_s2_dressed ci_electronic_energy_dressed
SOFT_TOUCH psi_ref_coef_diagonalized psi_ref_energy_diagonalized! psi_coef_energy_diagonalized
call H_apply_mrcepa_PT2(pt2, norm_pert, H_pert_diag, N_st)
SOFT_TOUCH psi_ref_coef_diagonalized psi_ref_energy_diagonalized
call H_apply_mrcepa_PT2(pt2_redundant, norm_pert, H_pert_diag, N_st)
print * ,'Computing the remaining contribution'
threshold_selectors = 1.d0
threshold_generators = 0.999d0
N_det_generators = N_det_non_ref + N_det_ref
N_det_selectors = N_det_non_ref + N_det_ref
@ -125,17 +165,15 @@ subroutine run_pt2(N_st,energy)
enddo
SOFT_TOUCH N_det_selectors psi_selectors_coef psi_selectors N_det_generators psi_det_generators psi_coef_generators ci_eigenvectors_dressed ci_eigenvectors_s2_dressed ci_electronic_energy_dressed
SOFT_TOUCH psi_ref_coef_diagonalized psi_ref_energy_diagonalized! psi_coef_energy_diagonalized
call H_apply_mrcepa_PT2(pt3, norm_pert, H_pert_diag, N_st)
SOFT_TOUCH psi_ref_coef_diagonalized psi_ref_energy_diagonalized
call H_apply_mrcepa_PT2(pt2, norm_pert, H_pert_diag, N_st)
!!!!!!!!!!!!!!!!
print *, "3-2 :",pt3, pt2
print *, lambda_mrcc_pt3(0), N_det, N_det_ref, psi_coef(1,1), psi_ref_coef(1,1)
pt2 = pt3 - pt2
print *, "Redundant PT2 :",pt2_redundant
print *, "Full PT2 :",pt2
print *, lambda_mrcc_kept(0), N_det, N_det_ref, psi_coef(1,1), psi_ref_coef(1,1)
pt2 = pt2 - pt2_redundant
print *, 'Final step'
print *, 'N_det = ', N_det
@ -146,24 +184,59 @@ subroutine run_pt2(N_st,energy)
print *, '-----'
call ezfio_set_full_ci_energy_pt2(energy+pt2)
deallocate(pt2,norm_pert)
! call ezfio_set_full_ci_energy_pt2(energy+pt2)
end
subroutine run_pt2(N_st,energy)
implicit none
integer :: i,j,k
integer, intent(in) :: N_st
double precision, intent(in) :: energy(N_st)
double precision :: pt2(N_st)
double precision :: norm_pert(N_st),H_pert_diag(N_st)
pt2 = 0d0
!if(lambda_mrcc_pt2(0) == 0) return
print*,'Last iteration only to compute the PT2'
N_det_generators = N_det_cas
N_det_selectors = N_det_non_ref
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)
do i=1,N_det_generators
do k=1,N_int
psi_det_generators(k,1,i) = psi_ref(k,1,i)
psi_det_generators(k,2,i) = psi_ref(k,2,i)
enddo
do k=1,N_st
psi_coef_generators(i,k) = psi_ref_coef(i,k)
enddo
enddo
do i=1,N_det
do k=1,N_int
psi_selectors(k,1,i) = psi_det_sorted(k,1,i)
psi_selectors(k,2,i) = psi_det_sorted(k,2,i)
enddo
do k=1,N_st
psi_selectors_coef(i,k) = psi_coef_sorted(i,k)
enddo
enddo
call write_double(6,ci_energy(1),"Initial CI energy")
end
SOFT_TOUCH N_det_selectors psi_selectors_coef psi_selectors N_det_generators psi_det_generators psi_coef_generators ci_eigenvectors_dressed ci_eigenvectors_s2_dressed ci_electronic_energy_dressed
SOFT_TOUCH psi_ref_coef_diagonalized psi_ref_energy_diagonalized
call H_apply_mrcepa_PT2(pt2, norm_pert, H_pert_diag, N_st)
! call ezfio_set_full_ci_energy_pt2(energy+pt2)
print *, 'Final step'
print *, 'N_det = ', N_det
print *, 'N_states = ', N_states
print *, 'PT2 = ', pt2
print *, 'E = ', energy
print *, 'E+PT2 = ', energy+pt2
print *, '-----'
end

8
src/Determinants/slater_rules.irp.f
View File

@ -1935,3 +1935,11 @@ subroutine get_phase(key1,key2,phase,Nint)
!DIR$ FORCEINLINE
call get_excitation(key1, key2, exc, degree, phase, Nint)
end
BEGIN_PROVIDER [ double precision, psi_energy, (N_states) ]
implicit none
BEGIN_DOC
! Energy of the current wave function
END_DOC
call u_0_H_u_0_nstates(psi_energy,psi_coef,N_det,psi_det,N_int,N_states,psi_det_size)
END_PROVIDER