LCPQ/quantum_package
10
0
Fork 0
mirror of https://github.com/LCPQ/quantum_package synced 2024-12-22 12:23:48 +01:00
Code Issues Releases Wiki Activity

minor modifs in FOBOCI

Browse Source
This commit is contained in:
Emmanuel Giner 2017-03-20 16:21:00 +01:00
parent cf55bf8093
commit 3b8239976a
12 changed files with 121 additions and 160 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
plugins/DDCI_selected/NEEDED_CHILDREN_MODULES
View File

@ -1 +1 @@
Perturbation Selectors_full Generators_CAS Davidson
Perturbation Selectors_full Generators_CAS Davidson Psiref_CAS

2
plugins/DDCI_selected/ddci.irp.f
View File

@ -5,7 +5,7 @@ program ddci
double precision, allocatable :: pt2(:), norm_pert(:), H_pert_diag(:),E_before(:)
integer :: N_st, degree
N_st = N_states_diag
N_st = N_states
allocate (pt2(N_st), norm_pert(N_st),H_pert_diag(N_st),E_before(N_st))
character*(64) :: perturbation

46
plugins/FOBOCI/density_matrix.irp.f
View File

@ -1,12 +1,12 @@
BEGIN_PROVIDER [ double precision, one_body_dm_mo_alpha_generators_restart, (mo_tot_num_align,mo_tot_num) ]
&BEGIN_PROVIDER [ double precision, one_body_dm_mo_beta_generators_restart, (mo_tot_num_align,mo_tot_num) ]
&BEGIN_PROVIDER [ double precision, norm_generators_restart]
&BEGIN_PROVIDER [ double precision, norm_generators_restart, (N_states)]
implicit none
BEGIN_DOC
! Alpha and beta one-body density matrix for the generators restart
END_DOC
integer :: j,k,l,m
integer :: j,k,l,m,istate
integer :: occ(N_int*bit_kind_size,2)
double precision :: ck, cl, ckl
double precision :: phase
@ -14,29 +14,37 @@
integer :: exc(0:2,2,2),n_occ_alpha
double precision, allocatable :: tmp_a(:,:), tmp_b(:,:)
integer :: degree_respect_to_HF_k
integer :: degree_respect_to_HF_l,index_ref_generators_restart
double precision :: inv_coef_ref_generators_restart
integer :: degree_respect_to_HF_l,index_ref_generators_restart(N_states)
double precision :: inv_coef_ref_generators_restart(N_states)
integer :: i
print*, 'providing the one_body_dm_mo_alpha_generators_restart'
do i = 1, N_det_generators_restart
! Find the reference determinant for intermediate normalization
call get_excitation_degree(ref_generators_restart,psi_det_generators_restart(1,1,i),degree,N_int)
if(degree == 0)then
index_ref_generators_restart = i
inv_coef_ref_generators_restart = 1.d0/psi_coef_generators_restart(i,1)
exit
endif
do istate = 1, N_states
do i = 1, N_det_generators_restart
! Find the reference determinant for intermediate normalization
call get_excitation_degree(ref_generators_restart(1,1,istate),psi_det_generators_restart(1,1,i),degree,N_int)
if(degree == 0)then
index_ref_generators_restart(istate) = i
inv_coef_ref_generators_restart(istate) = 1.d0/psi_coef_generators_restart(i,istate)
exit
endif
enddo
enddo
norm_generators_restart = 0.d0
do i = 1, N_det_generators_restart
psi_coef_generators_restart(i,1) = psi_coef_generators_restart(i,1) * inv_coef_ref_generators_restart
norm_generators_restart += psi_coef_generators_restart(i,1)**2
do istate = 1, N_states
do i = 1, N_det_generators_restart
psi_coef_generators_restart(i,istate) = psi_coef_generators_restart(i,istate) * inv_coef_ref_generators_restart(istate)
norm_generators_restart(istate) += psi_coef_generators_restart(i,istate)**2
enddo
enddo
double precision :: inv_norm
inv_norm = 1.d0/dsqrt(norm_generators_restart)
do i = 1, N_det_generators_restart
psi_coef_generators_restart(i,1) = psi_coef_generators_restart(i,1) * inv_norm
double precision :: inv_norm(N_States)
do istate = 1, N_states
inv_norm(istate) = 1.d0/dsqrt(norm_generators_restart(istate))
enddo
do istate = 1, N_states
do i = 1, N_det_generators_restart
psi_coef_generators_restart(i,istate) = psi_coef_generators_restart(i,istate) * inv_norm(istate)
enddo
enddo

62
plugins/FOBOCI/dress_simple.irp.f
View File

@ -107,7 +107,6 @@ subroutine is_a_good_candidate(threshold,is_ok,e_pt2,verbose,exit_loop,is_ok_per
!enddo
!soft_touch psi_selectors psi_selectors_coef
!if(do_it_perturbative)then
print*, 'is_ok_perturbative',is_ok_perturbative
if(is_ok.or.is_ok_perturbative)then
N_det = N_det_generators
do m = 1, N_states
@ -117,7 +116,6 @@ subroutine is_a_good_candidate(threshold,is_ok,e_pt2,verbose,exit_loop,is_ok_per
psi_det(l,2,k) = psi_det_generators_input(l,2,k)
enddo
psi_coef(k,m) = psi_coef_diagonalized_tmp(k,m)
print*, 'psi_coef(k,m)',psi_coef(k,m)
enddo
enddo
soft_touch psi_det psi_coef N_det
@ -150,7 +148,7 @@ subroutine dress_H_matrix_from_psi_det_input(psi_det_generators_input,Ndet_gener
double precision, intent(inout) :: dressed_H_matrix(Ndet_generators, Ndet_generators)
integer :: i,j,degree,index_ref_generators_restart,i_count,k,i_det_no_ref
integer :: i,j,degree,index_ref_generators_restart(N_states),i_count,k,i_det_no_ref
double precision :: eigvalues(Ndet_generators), eigvectors(Ndet_generators,Ndet_generators),hij
double precision :: psi_coef_ref(Ndet_generators,N_states),diag_h_mat_average,diag_h_mat_no_ref_average
logical :: is_a_ref_det(Ndet_generators)
@ -168,17 +166,19 @@ subroutine dress_H_matrix_from_psi_det_input(psi_det_generators_input,Ndet_gener
enddo
integer :: istate
do istate = 1, N_states
do i = 1, Ndet_generators
call get_excitation_degree(ref_generators_restart(1,1,istate),psi_det_generators_input(1,1,i),degree,N_int)
if(degree == 0)then
index_ref_generators_restart(istate) = i
exit
endif
enddo
enddo
do i = 1, Ndet_generators
call get_excitation_degree(ref_generators_restart,psi_det_generators_input(1,1,i),degree,N_int)
if(degree == 0)then
index_ref_generators_restart = i
endif
do j = 1, Ndet_generators
call i_h_j(psi_det_generators_input(1,1,j),psi_det_generators_input(1,1,i),N_int,hij) ! Fill the zeroth order H matrix
if(i==j)then
call debug_det(psi_det_generators_input(1,1,i),N_int)
print*, hij
endif
dressed_H_matrix(i,j) = hij
enddo
enddo
@ -189,15 +189,21 @@ subroutine dress_H_matrix_from_psi_det_input(psi_det_generators_input,Ndet_gener
i_det_no_ref +=1
diag_h_mat_average+=dressed_H_matrix(i,i)
enddo
double precision :: average_ref_h_mat
average_ref_h_mat = 0.d0
do istate = 1, N_states
average_ref_h_mat += dressed_H_matrix(index_ref_generators_restart(istate),index_ref_generators_restart(istate))
enddo
average_ref_h_mat = 1.d0/dble(N_states)
diag_h_mat_average = diag_h_mat_average/dble(i_det_no_ref)
print*,'diag_h_mat_average = ',diag_h_mat_average
print*,'ref h_mat = ',dressed_H_matrix(index_ref_generators_restart,index_ref_generators_restart)
print*,'ref h_mat average = ',average_ref_h_mat
integer :: number_of_particles, number_of_holes
! Filter the the MLCT that are higher than 27.2 eV in energy with respect to the reference determinant
do i = 1, Ndet_generators
if(is_a_ref_det(i))cycle
if(number_of_holes(psi_det_generators_input(1,1,i)).eq.0 .and. number_of_particles(psi_det_generators_input(1,1,i)).eq.1)then
if(diag_h_mat_average - dressed_H_matrix(index_ref_generators_restart,index_ref_generators_restart) .gt.2.d0)then
if(diag_h_mat_average - average_ref_h_mat .gt.2.d0)then
is_ok = .False.
exit_loop = .True.
return
@ -206,7 +212,7 @@ subroutine dress_H_matrix_from_psi_det_input(psi_det_generators_input,Ndet_gener
! Filter the the LMCT that are higher than 54.4 eV in energy with respect to the reference determinant
if(number_of_holes(psi_det_generators_input(1,1,i)).eq.1 .and. number_of_particles(psi_det_generators_input(1,1,i)).eq.0)then
if(diag_h_mat_average - dressed_H_matrix(index_ref_generators_restart,index_ref_generators_restart) .gt.2.d0)then
if(diag_h_mat_average - average_ref_h_mat .gt.1.d0)then
is_ok = .False.
return
endif
@ -214,7 +220,7 @@ subroutine dress_H_matrix_from_psi_det_input(psi_det_generators_input,Ndet_gener
exit
enddo
call lapack_diagd(eigvalues,eigvectors,dressed_H_matrix,Ndet_generators,Ndet_generators) ! Diagonalize the Dressed_H_matrix
call lapack_diagd(eigvalues,eigvectors,dressed_H_matrix,Ndet_generators,Ndet_generators) ! Diagonalize the naked matrix
double precision :: s2(N_det_generators),E_ref(N_states)
integer :: i_state(N_states)
@ -238,18 +244,12 @@ subroutine dress_H_matrix_from_psi_det_input(psi_det_generators_input,Ndet_gener
do i = 1, N_states
i_state(i) = i
E_ref(i) = eigvalues(i)
print*, 'E_ref(i)',E_ref(i)
enddo
endif
do i = 1,N_states
print*,'i_state = ',i_state(i)
enddo
do k = 1, N_states
print*,'state ',k
do i = 1, Ndet_generators
psi_coef_diagonalized_tmp(i,k) = eigvectors(i,i_state(k)) / eigvectors(index_ref_generators_restart,i_state(k))
psi_coef_diagonalized_tmp(i,k) = eigvectors(i,i_state(k)) / eigvectors(index_ref_generators_restart(k),i_state(k))
psi_coef_ref(i,k) = eigvectors(i,i_state(k))
print*,'psi_coef_ref(i) = ',psi_coef_ref(i,k)
enddo
enddo
if(verbose)then
@ -262,7 +262,7 @@ subroutine dress_H_matrix_from_psi_det_input(psi_det_generators_input,Ndet_gener
do k = 1, N_states
print*,'state ',k
do i = 1, Ndet_generators
print*,'coef, <I|H|I> = ',psi_coef_diagonalized_tmp(i,k),dressed_H_matrix(i,i)-dressed_H_matrix(index_ref_generators_restart,index_ref_generators_restart),is_a_ref_det(i)
print*,'coef, <I|H|I> = ',psi_coef_diagonalized_tmp(i,k),dressed_H_matrix(i,i)-dressed_H_matrix(index_ref_generators_restart(k),index_ref_generators_restart(k)),is_a_ref_det(i)
enddo
enddo
endif
@ -283,18 +283,20 @@ subroutine dress_H_matrix_from_psi_det_input(psi_det_generators_input,Ndet_gener
call lapack_diagd(eigvalues,eigvectors,dressed_H_matrix,Ndet_generators,Ndet_generators) ! Diagonalize the Dressed_H_matrix
integer :: i_good_state(0:N_states)
i_good_state(0) = 0
do i = 1, Ndet_generators
do k = 1, N_states
! print*,'state',k
do i = 1, Ndet_generators
! State following
do k = 1, N_states
accu = 0.d0
do j =1, Ndet_generators
print*,'',eigvectors(j,i) , psi_coef_ref(j,k)
accu += eigvectors(j,i) * psi_coef_ref(j,k)
enddo
print*,'accu = ',accu
! print*,i,accu
if(dabs(accu).ge.0.60d0)then
i_good_state(0) +=1
i_good_state(i_good_state(0)) = i
print*, 'state, ovrlap',k,i,accu
exit
endif
enddo
if(i_good_state(0)==N_states)then
@ -309,14 +311,14 @@ subroutine dress_H_matrix_from_psi_det_input(psi_det_generators_input,Ndet_gener
accu = 0.d0
do k = 1, N_states
do i = 1, Ndet_generators
psi_coef_diagonalized_tmp(i,k) = eigvectors(i,i_state(k)) / eigvectors(index_ref_generators_restart,i_state(k))
psi_coef_diagonalized_tmp(i,k) = eigvectors(i,i_state(k)) / eigvectors(index_ref_generators_restart(k),i_state(k))
enddo
enddo
if(verbose)then
do k = 1, N_states
print*,'state ',k
do i = 1, Ndet_generators
print*,'coef, <I|H+Delta H|I> = ',psi_coef_diagonalized_tmp(i,k),dressed_H_matrix(i,i)-dressed_H_matrix(index_ref_generators_restart,index_ref_generators_restart),is_a_ref_det(i)
print*,'coef, <I|H+Delta H|I> = ',psi_coef_diagonalized_tmp(i,k),dressed_H_matrix(i,i)-dressed_H_matrix(index_ref_generators_restart(k),index_ref_generators_restart(k)),is_a_ref_det(i)
enddo
enddo
endif
@ -338,7 +340,7 @@ subroutine dress_H_matrix_from_psi_det_input(psi_det_generators_input,Ndet_gener
do i = 1, Ndet_generators
if(is_a_ref_det(i))cycle
do k = 1, N_states
print*, psi_coef_diagonalized_tmp(i,k),threshold_perturbative
! print*, psi_coef_diagonalized_tmp(i,k),threshold_perturbative
if(dabs(psi_coef_diagonalized_tmp(i,k)) .gt.threshold_perturbative)then
is_ok_perturbative = .False.
exit

23
plugins/FOBOCI/foboci_lmct_mlct_threshold_old.irp.f
View File

@ -46,6 +46,7 @@ subroutine FOBOCI_lmct_mlct_old_thr(iter)
lmct = .True.
integer :: i_hole_osoci
i_hole_osoci = list_inact(i)
! if(i_hole_osoci.ne.26)cycle
print*,'--------------------------'
! First set the current generators to the one of restart
call check_symetry(i_hole_osoci,thr,test_sym)
@ -55,11 +56,6 @@ subroutine FOBOCI_lmct_mlct_old_thr(iter)
print*,'i_hole_osoci = ',i_hole_osoci
call create_restart_and_1h(i_hole_osoci)
call set_generators_to_psi_det
print*,'Passed set generators'
integer :: m
do m = 1, N_det_generators
call debug_det(psi_det_generators(1,1,m),N_int)
enddo
call set_bitmask_particl_as_input(reunion_of_bitmask)
call set_bitmask_hole_as_input(reunion_of_bitmask)
double precision :: e_pt2
@ -88,9 +84,9 @@ subroutine FOBOCI_lmct_mlct_old_thr(iter)
call set_bitmask_hole_as_input(reunion_of_bitmask)
call all_single(e_pt2)
! call make_s2_eigenfunction_first_order
threshold_davidson = 1.d-6
soft_touch threshold_davidson davidson_criterion
call diagonalize_ci
! threshold_davidson = 1.d-6
! soft_touch threshold_davidson davidson_criterion
! call diagonalize_ci
double precision :: hkl
call provide_matrix_dressing(dressing_matrix,n_det_generators,psi_det_generators)
hkl = dressing_matrix(1,1)
@ -123,6 +119,7 @@ subroutine FOBOCI_lmct_mlct_old_thr(iter)
do i = 1, n_virt_orb
integer :: i_particl_osoci
i_particl_osoci = list_virt(i)
! cycle
print*,'--------------------------'
! First set the current generators to the one of restart
@ -157,11 +154,11 @@ subroutine FOBOCI_lmct_mlct_old_thr(iter)
enddo
enddo
call all_single(e_pt2)
call make_s2_eigenfunction_first_order
threshold_davidson = 1.d-6
soft_touch threshold_davidson davidson_criterion
call diagonalize_ci
! call make_s2_eigenfunction_first_order
! threshold_davidson = 1.d-6
! soft_touch threshold_davidson davidson_criterion
!
! call diagonalize_ci
deallocate(dressing_matrix)
else
if(exit_loop)then

20
plugins/FOBOCI/generators_restart_save.irp.f
View File

@ -21,23 +21,19 @@ END_PROVIDER
BEGIN_PROVIDER [ integer(bit_kind), psi_det_generators_restart, (N_int,2,N_det_generators_restart) ]
&BEGIN_PROVIDER [ integer(bit_kind), ref_generators_restart, (N_int,2) ]
&BEGIN_PROVIDER [ integer(bit_kind), ref_generators_restart, (N_int,2,N_states) ]
&BEGIN_PROVIDER [ double precision, psi_coef_generators_restart, (N_det_generators_restart,N_states) ]
implicit none
BEGIN_DOC
! read wf
!
END_DOC
integer :: i, k
integer :: i, k,j
integer, save :: ifirst = 0
double precision, allocatable :: psi_coef_read(:,:)
print*, ' Providing psi_det_generators_restart'
if(ifirst == 0)then
call read_dets(psi_det_generators_restart,N_int,N_det_generators_restart)
do k = 1, N_int
ref_generators_restart(k,1) = psi_det_generators_restart(k,1,1)
ref_generators_restart(k,2) = psi_det_generators_restart(k,2,1)
enddo
allocate (psi_coef_read(N_det_generators_restart,N_states))
call ezfio_get_determinants_psi_coef(psi_coef_read)
do k = 1, N_states
@ -45,6 +41,18 @@ END_PROVIDER
psi_coef_generators_restart(i,k) = psi_coef_read(i,k)
enddo
enddo
do k = 1, N_states
do i = 1, N_det_generators_restart
if(dabs(psi_coef_generators_restart(i,k)).gt.0.5d0)then
do j = 1, N_int
ref_generators_restart(j,1,k) = psi_det_generators_restart(j,1,i)
ref_generators_restart(j,2,k) = psi_det_generators_restart(j,2,i)
enddo
exit
endif
enddo
call debug_det(ref_generators_restart(1,1,k),N_int)
enddo
ifirst = 1
deallocate(psi_coef_read)
else

48
plugins/FOBOCI/routines_foboci.irp.f
View File

@ -2,7 +2,7 @@ subroutine set_intermediate_normalization_lmct_old(norm,i_hole)
implicit none
integer, intent(in) :: i_hole
double precision, intent(out) :: norm(N_states)
integer :: i,j,degree,index_ref_generators_restart,k
integer :: i,j,degree,index_ref_generators_restart(N_states),k
integer:: number_of_holes,n_h, number_of_particles,n_p
integer, allocatable :: index_one_hole(:),index_one_hole_one_p(:),index_two_hole_one_p(:),index_two_hole(:)
integer, allocatable :: index_one_p(:)
@ -24,16 +24,18 @@ subroutine set_intermediate_normalization_lmct_old(norm,i_hole)
n_good_hole = 0
! Find the one holes and one hole one particle
is_a_ref_det = .False.
integer :: istate
do istate = 1, N_States
do i = 1, N_det
! Find the reference determinant for intermediate normalization
call get_excitation_degree(ref_generators_restart(1,1,istate),psi_det(1,1,i),degree,N_int)
if(degree == 0)then
index_ref_generators_restart(istate) = i
inv_coef_ref_generators_restart(istate) = 1.d0/psi_coef(i,istate)
endif
enddo
enddo
do i = 1, N_det
! Find the reference determinant for intermediate normalization
call get_excitation_degree(ref_generators_restart,psi_det(1,1,i),degree,N_int)
if(degree == 0)then
index_ref_generators_restart = i
do k = 1, N_states
inv_coef_ref_generators_restart(k) = 1.d0/psi_coef(i,k)
enddo
endif
! Find all the determinants present in the reference wave function
do j = 1, N_det_generators_restart
call get_excitation_degree(psi_det(1,1,i),psi_det_generators_restart(1,1,j),degree,N_int)
@ -67,7 +69,7 @@ subroutine set_intermediate_normalization_lmct_old(norm,i_hole)
do k = 1,N_states
print*,'state ',k
do i = 1, n_good_hole
print*,'psi_coef(index_good_hole) = ',psi_coef(index_good_hole(i),k)/psi_coef(index_ref_generators_restart,k)
print*,'psi_coef(index_good_hole) = ',psi_coef(index_good_hole(i),k)/psi_coef(index_ref_generators_restart(k),k)
enddo
print*,''
enddo
@ -110,7 +112,7 @@ subroutine set_intermediate_normalization_mlct_old(norm,i_particl)
implicit none
integer, intent(in) :: i_particl
double precision, intent(out) :: norm(N_states)
integer :: i,j,degree,index_ref_generators_restart,k
integer :: i,j,degree,index_ref_generators_restart(N_states),k
integer:: number_of_holes,n_h, number_of_particles,n_p
integer, allocatable :: index_one_hole(:),index_one_hole_one_p(:),index_two_hole_one_p(:),index_two_hole(:)
integer, allocatable :: index_one_p(:),index_one_hole_two_p(:)
@ -139,16 +141,18 @@ subroutine set_intermediate_normalization_mlct_old(norm,i_particl)
! Find the one holes and one hole one particle
i_count = 0
is_a_ref_det = .False.
do i = 1, N_det
call get_excitation_degree(ref_generators_restart,psi_det(1,1,i),degree,N_int)
if(degree == 0)then
index_ref_generators_restart = i
do k = 1, N_states
inv_coef_ref_generators_restart(k) = 1.d0/psi_coef(i,k)
enddo
! cycle
endif
integer :: istate
do istate = 1, N_states
do i = 1, N_det
call get_excitation_degree(ref_generators_restart(1,1,istate),psi_det(1,1,i),degree,N_int)
if(degree == 0)then
index_ref_generators_restart(istate) = i
inv_coef_ref_generators_restart(istate) = 1.d0/psi_coef(i,istate)
endif
enddo
enddo
do i = 1, N_det
! Find all the determinants present in the reference wave function
do j = 1, N_det_generators_restart
call get_excitation_degree(psi_det(1,1,i),psi_det_generators_restart(1,1,j),degree,N_int)
@ -184,7 +188,7 @@ subroutine set_intermediate_normalization_mlct_old(norm,i_particl)
do k = 1, N_states
print*,'state ',k
do i = 1, n_good_particl
print*,'psi_coef(index_good_particl,1) = ',psi_coef(index_good_particl(i),k)/psi_coef(index_ref_generators_restart,k)
print*,'psi_coef(index_good_particl,1) = ',psi_coef(index_good_particl(i),k)/psi_coef(index_ref_generators_restart(k),k)
enddo
print*,''
enddo

8
plugins/MRPT/MRPT_Utils.main.irp.f
View File

@ -32,13 +32,13 @@ subroutine routine_3
call save_wavefunction_general(N_det_ref,N_states_diag_heff,psi_ref,N_det_ref,CI_dressed_pt2_new_eigenvectors)
endif
if(N_states.gt.1)then
print*, 'Energy differences : E(0) - E(i)'
print*, 'Energy differences : E(i) - E(0)'
do i = 2, N_States
print*,'State',i
write(*,'(A12,X,I3,A3,XX,F20.16)') ' S^2 ', i,' = ', CI_dressed_pt2_new_eigenvectors_s2(i)
write(*,'(A12,X,I3,A3,XX,F20.16)') 'Variational ', i,' = ', psi_ref_average_value(1) - psi_ref_average_value(i)
write(*,'(A12,X,I3,A3,XX,F20.16)') 'Perturbative', i,' = ', psi_ref_average_value(1)+second_order_pt_new(1) - (psi_ref_average_value(i)+second_order_pt_new(i))
write(*,'(A12,X,I3,A3,XX,F20.16)') 'Dressed ', i,' = ', CI_dressed_pt2_new_energy(1) - CI_dressed_pt2_new_energy(i)
write(*,'(A12,X,I3,A3,XX,F20.16)') 'Variational ', i,' = ', -(psi_ref_average_value(1) - psi_ref_average_value(i))
write(*,'(A12,X,I3,A3,XX,F20.16)') 'Perturbative', i,' = ', -(psi_ref_average_value(1)+second_order_pt_new(1) - (psi_ref_average_value(i)+second_order_pt_new(i)))
write(*,'(A12,X,I3,A3,XX,F20.16)') 'Dressed ', i,' = ', -( CI_dressed_pt2_new_energy(1) - CI_dressed_pt2_new_energy(i) )
enddo
endif

32
plugins/MRPT_Utils/energies_cas.irp.f
View File

@ -688,10 +688,6 @@ END_PROVIDER
call get_mono_excitation(psi_in_out(1,1,i),psi_ref(1,1,i),exc,phase,N_int)
do j = 1, n_states
psi_in_out_coef(i,j) = psi_ref_coef(i,j)* hij * phase
! if(orb_i == 5 .and. orb_v == 20)then
if(orb_i == 2 .and. orb_v == 6 )then
print*, i, ispin
endif
norm(j,ispin) += psi_in_out_coef(i,j) * psi_in_out_coef(i,j)
enddo
enddo
@ -702,12 +698,7 @@ END_PROVIDER
one_anhil_one_creat_inact_virt_norm(iorb,vorb,j,ispin) = 0.d0
else
norm_no_inv(j,ispin) = norm(j,ispin)
! one_anhil_one_creat_inact_virt_norm(iorb,vorb,j,ispin) = 1.d0 / norm(j,ispin)
norm(j,ispin) = 1.d0/dsqrt(norm(j,ispin))
! if(orb_i == 5 .and. orb_v == 20)then
if(orb_i == 2 .and. orb_v == 6 )then
print*,ispin ,norm(j,ispin)
endif
endif
enddo
integer :: iorb_annil,hole_particle,spin_exc,orb
@ -721,12 +712,8 @@ END_PROVIDER
do i = 1, N_det_ref
do j = 1, N_int
! psi_in_out(j,1,i) = iand(psi_in_out(j,1,i),cas_bitmask(j,1,1))
! psi_in_out(j,2,i) = iand(psi_in_out(j,2,i),cas_bitmask(j,1,1))
psi_in_out(j,1,i) = psi_active(j,1,i)
psi_in_out(j,2,i) = psi_active(j,2,i)
! psi_in_out(j,1,i) = psi_ref(j,1,i)
! psi_in_out(j,2,i) = psi_ref(j,2,i)
enddo
enddo
do state_target = 1, N_states
@ -740,30 +727,11 @@ END_PROVIDER
do state_target = 1, N_states
if((norm_no_inv(state_target,1) + norm_no_inv(state_target,2)) .ne. 0.d0)then
one_anhil_one_creat_inact_virt(iorb,vorb,state_target) = energy_cas_dyall_no_exchange(state_target) - &
! 0.5d0 * (energies_alpha_beta(state_target,1) + energies_alpha_beta(state_target,2))
( energies_alpha_beta(state_target,1) + energies_alpha_beta(state_target,2) ) &
/( norm_bis(state_target,1) + norm_bis(state_target,2) )
else
one_anhil_one_creat_inact_virt(iorb,vorb,state_target) = 0.d0
endif
if(dabs(dabs(one_anhil_one_creat_inact_virt(iorb,vorb,state_target)) - 1.30584271462d0) < 1.d-11)then
! if(dabs(dabs(one_anhil_one_creat_inact_virt(iorb,vorb,state_target)) - 1.29269686324d0) < 1.d-11)then
print*, ''
print*, orb_i,orb_v
print*, energy_cas_dyall_no_exchange(1) - energies_alpha_beta(state_target,1) !/ norm_bis(state_target,1)
print*, energy_cas_dyall_no_exchange(1) - energies_alpha_beta(state_target,2) !/ norm_bis(state_target,2)
print*, fock_core_inactive_total_spin_trace(orb_i,1)
print*, fock_virt_total_spin_trace(orb_v,1)
print*, one_anhil_one_creat_inact_virt(iorb,vorb,state_target)
print*, ''
endif
if(dabs(one_anhil_one_creat_inact_virt(iorb,vorb,state_target)) .gt. 1.d-10)then
write(*,'(F11.8)'), one_anhil_one_creat_inact_virt(iorb,vorb,state_target)
! if(dabs(one_anhil_one_creat_inact_virt(iorb,vorb,state_target)) .lt. 1.d-2)then
! one_anhil_one_creat_inact_virt(iorb,vorb,state_target) = 0.d0
! print*, orb_i,orb_v
! endif
endif
enddo
enddo
enddo

9
plugins/MRPT_Utils/mrpt_utils.irp.f
View File

@ -305,8 +305,8 @@ END_PROVIDER
END_PROVIDER
BEGIN_PROVIDER [ double precision, CI_electronic_dressed_pt2_new_energy, (N_states_diag_heff) ]
&BEGIN_PROVIDER [ double precision, CI_dressed_pt2_new_eigenvectors, (N_det_ref,N_states_diag_heff) ]
&BEGIN_PROVIDER [ double precision, CI_dressed_pt2_new_eigenvectors_s2, (N_states_diag_heff) ]
&BEGIN_PROVIDER [ double precision, CI_dressed_pt2_new_eigenvectors, (N_det_ref,N_states) ]
&BEGIN_PROVIDER [ double precision, CI_dressed_pt2_new_eigenvectors_s2, (N_states) ]
BEGIN_DOC
! Eigenvectors/values of the CI matrix
END_DOC
@ -394,6 +394,7 @@ END_PROVIDER
psi_tmp(i) = eigenvectors(i,iorder(1))
enddo
CI_electronic_dressed_pt2_new_energy(i_state) = eigenvalues(iorder(1))
print*, 'CI_electronic_dressed_pt2_new_energy',CI_electronic_dressed_pt2_new_energy(i_state)
call u_0_S2_u_0(CI_dressed_pt2_new_eigenvectors_s2(i_state),psi_tmp,N_det_ref,psi_det,N_int,1,N_det_ref)
print*,'S^2 = ', CI_dressed_pt2_new_eigenvectors_s2(i_state)
enddo
@ -502,7 +503,7 @@ END_PROVIDER
END_PROVIDER
BEGIN_PROVIDER [ double precision, CI_dressed_pt2_new_energy, (N_states_diag_heff) ]
BEGIN_PROVIDER [ double precision, CI_dressed_pt2_new_energy, (N_states) ]
implicit none
BEGIN_DOC
! N_states lowest eigenvalues of the CI matrix
@ -511,7 +512,7 @@ BEGIN_PROVIDER [ double precision, CI_dressed_pt2_new_energy, (N_states_diag_hef
integer :: j
character*(8) :: st
call write_time(output_determinants)
do j=1,N_states_diag_heff
do j=1,N_states
CI_dressed_pt2_new_energy(j) = CI_electronic_dressed_pt2_new_energy(j) + nuclear_repulsion
write(st,'(I4)') j
call write_double(output_determinants,CI_dressed_pt2_new_energy(j),'Energy of state '//trim(st))

2
src/Davidson/diagonalize_restart_and_save_all_states.irp.f
View File

@ -9,7 +9,7 @@ subroutine routine
implicit none
call diagonalize_CI
print*,'N_det = ',N_det
call save_wavefunction_general(N_det,N_states_diag,psi_det_sorted,size(psi_coef_sorted,1),psi_coef_sorted)
call save_wavefunction_general(N_det,N_states,psi_det_sorted,size(psi_coef_sorted,1),psi_coef_sorted)

27
src/Determinants/diagonalize_restart_and_save_two_states.irp.f
View File

@ -1,27 +0,0 @@
program diag_and_save
implicit none
read_wf = .True.
touch read_wf
call routine
end
subroutine routine
implicit none
integer :: igood_state_1,igood_state_2
double precision, allocatable :: psi_coef_tmp(:,:)
integer :: i
print*,'N_det = ',N_det
!call diagonalize_CI
write(*,*)'Which couple of states would you like to save ?'
read(5,*)igood_state_1,igood_state_2
allocate(psi_coef_tmp(n_det,2))
do i = 1, N_det
psi_coef_tmp(i,1) = psi_coef(i,igood_state_1)
psi_coef_tmp(i,2) = psi_coef(i,igood_state_2)
enddo
call save_wavefunction_general(N_det,2,psi_det,n_det,psi_coef_tmp)
deallocate(psi_coef_tmp)
end