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

Corrected diagonalize_CI

Browse Source
This commit is contained in:
Anthony Scemama 2016-09-30 15:12:17 +02:00
parent 7219dda33b
commit cbdc8f68d7
8 changed files with 180 additions and 72 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

8
plugins/CAS_SD/H_apply.irp.f
View File

@ -3,6 +3,7 @@ BEGIN_SHELL [ /usr/bin/env python ]
from generate_h_apply import *
s = H_apply("CAS_SD")
s.unset_skip()
print s
s = H_apply("CAS_SD_selected_no_skip")
@ -12,6 +13,7 @@ print s
s = H_apply("CAS_SD_selected")
s.set_selection_pt2("epstein_nesbet_2x2")
s.unset_skip()
print s
s = H_apply("CAS_SD_PT2")
@ -22,13 +24,9 @@ print s
s = H_apply("CAS_S",do_double_exc=False)
print s
s = H_apply("CAS_S_selected_no_skip",do_double_exc=False)
s.set_selection_pt2("epstein_nesbet_2x2")
s.unset_skip()
print s
s = H_apply("CAS_S_selected",do_double_exc=False)
s.set_selection_pt2("epstein_nesbet_2x2")
s.unset_skip()
print s
s = H_apply("CAS_S_PT2",do_double_exc=False)

73
plugins/CAS_SD/cas_s_selected.irp.f
View File

@ -12,6 +12,7 @@ program full_ci
pt2 = 1.d0
diag_algorithm = "Lapack"
if (N_det > N_det_max) then
call diagonalize_CI
call save_wavefunction
@ -28,49 +29,84 @@ program full_ci
print *, 'E+PT2 = ', CI_energy+pt2
print *, '-----'
endif
double precision :: i_H_psi_array(N_states),diag_H_mat_elem,h,i_O1_psi_array(N_states)
double precision :: E_CI_before(N_states)
if(read_wf)then
call i_H_psi(psi_det(1,1,N_det),psi_det,psi_coef,N_int,N_det,psi_det_size,N_states,i_H_psi_array)
h = diag_H_mat_elem(psi_det(1,1,N_det),N_int)
selection_criterion = dabs(psi_coef(N_det,1) * (i_H_psi_array(1) - h * psi_coef(N_det,1))) * 0.1d0
soft_touch selection_criterion
endif
integer :: n_det_before
print*,'Beginning the selection ...'
E_CI_before(1:N_states) = CI_energy(1:N_states)
do while (N_det < N_det_max.and.maxval(abs(pt2(1:N_st))) > pt2_max)
call H_apply_CAS_S_selected_no_skip(pt2, norm_pert, H_pert_diag, N_st)
n_det_before = N_det
call H_apply_CAS_SD_selected(pt2, norm_pert, H_pert_diag, N_st)
PROVIDE psi_coef
PROVIDE psi_det
PROVIDE psi_det_sorted
if (N_det > N_det_max) then
psi_det = psi_det_sorted
psi_coef = psi_coef_sorted
N_det = N_det_max
soft_touch N_det psi_det psi_coef
endif
call diagonalize_CI
if (N_det > N_det_max) then
N_det = N_det_max
psi_det = psi_det_sorted
psi_coef = psi_coef_sorted
touch N_det psi_det psi_coef psi_det_sorted psi_coef_sorted psi_average_norm_contrib_sorted
endif
call save_wavefunction
print *, 'N_det = ', N_det
print *, 'N_states = ', N_states
print *, 'PT2 = ', pt2
print *, 'E = ', CI_energy
print *, 'E+PT2 = ', CI_energy+pt2
if(n_det_before == N_det)then
selection_criterion = selection_criterion * 0.5d0
endif
print *, 'N_det = ', N_det
print *, 'N_states = ', N_states
do k = 1, N_states
print*,'State ',k
print *, 'PT2 = ', pt2(k)
print *, 'E = ', CI_energy(k)
print *, 'E(before)+PT2 = ', E_CI_before(k)+pt2(k)
enddo
print *, '-----'
if(N_states.gt.1)then
print*,'Variational Energy difference'
do i = 2, N_states
print*,'Delta E = ',CI_energy(i) - CI_energy(1)
enddo
endif
if(N_states.gt.1)then
print*,'Variational + perturbative Energy difference'
do i = 2, N_states
print*,'Delta E = ',E_CI_before(i)+ pt2(i) - (E_CI_before(1) + pt2(1))
enddo
endif
E_CI_before(1:N_states) = CI_energy(1:N_states)
call ezfio_set_cas_sd_energy(CI_energy(1))
enddo
call diagonalize_CI
N_det = min(N_det_max,N_det)
touch N_det psi_det psi_coef
call diagonalize_CI
if(do_pt2_end)then
print*,'Last iteration only to compute the PT2'
threshold_selectors = 1.d0
threshold_generators = 0.999d0
call H_apply_CAS_S_PT2(pt2, norm_pert, H_pert_diag, N_st)
call H_apply_CAS_SD_PT2(pt2, norm_pert, H_pert_diag, N_st)
print *, 'Final step'
print *, 'N_det = ', N_det
print *, 'N_states = ', N_states
print *, 'PT2 = ', pt2
print *, 'E = ', CI_energy
print *, 'E+PT2 = ', CI_energy+pt2
print *, 'E = ', CI_energy(1:N_states)
print *, 'E+PT2 = ', CI_energy(1:N_states)+pt2(1:N_states)
print *, '-----'
call ezfio_set_cas_sd_energy_pt2(CI_energy(1)+pt2(1))
endif
integer :: exc_max, degree_min
exc_max = 0
print *, 'CAS determinants : ', N_det_cas
@ -79,6 +115,7 @@ program full_ci
call get_excitation_degree(psi_cas(1,1,k),psi_cas(1,1,i),degree,N_int)
exc_max = max(exc_max,degree)
enddo
print *, psi_coef_cas_diagonalized(i,:)
call debug_det(psi_cas(1,1,i),N_int)
print *, ''
enddo

81
plugins/CAS_SD/cas_sd.irp.f
View File

@ -1,7 +1,6 @@
program full_ci
implicit none
integer :: i,k
integer :: N_det_old
double precision, allocatable :: pt2(:), norm_pert(:), H_pert_diag(:)
@ -11,9 +10,9 @@ program full_ci
character*(64) :: perturbation
PROVIDE N_det_cas
N_det_old = 0
pt2 = 1.d0
diag_algorithm = "Lapack"
if (N_det > N_det_max) then
call diagonalize_CI
call save_wavefunction
@ -30,36 +29,68 @@ program full_ci
print *, 'E+PT2 = ', CI_energy+pt2
print *, '-----'
endif
double precision :: i_H_psi_array(N_states),diag_H_mat_elem,h,i_O1_psi_array(N_states)
double precision :: E_CI_before(N_states)
if(read_wf)then
call i_H_psi(psi_det(1,1,N_det),psi_det,psi_coef,N_int,N_det,psi_det_size,N_states,i_H_psi_array)
h = diag_H_mat_elem(psi_det(1,1,N_det),N_int)
selection_criterion = dabs(psi_coef(N_det,1) * (i_H_psi_array(1) - h * psi_coef(N_det,1))) * 0.1d0
soft_touch selection_criterion
endif
integer :: n_det_before
print*,'Beginning the selection ...'
E_CI_before(1:N_states) = CI_energy(1:N_states)
do while (N_det < N_det_max.and.maxval(abs(pt2(1:N_st))) > pt2_max)
N_det_old = N_det
n_det_before = N_det
call H_apply_CAS_SD(pt2, norm_pert, H_pert_diag, N_st)
PROVIDE psi_coef
PROVIDE psi_det
PROVIDE psi_det_sorted
if (N_det > N_det_max) then
psi_det = psi_det_sorted
psi_coef = psi_coef_sorted
N_det = N_det_max
soft_touch N_det psi_det psi_coef
endif
call diagonalize_CI
call save_wavefunction
print *, 'N_det = ', N_det
print *, 'N_states = ', N_states
print *, 'PT2 = ', pt2
print *, 'E = ', CI_energy
print *, 'E+PT2 = ', CI_energy+pt2
print *, '-----'
call ezfio_set_cas_sd_energy(CI_energy(1))
if (N_det == N_det_old) then
exit
endif
enddo
call diagonalize_CI
if (N_det > N_det_max) then
N_det = N_det_max
psi_det = psi_det_sorted
psi_coef = psi_coef_sorted
touch N_det psi_det psi_coef psi_det_sorted psi_coef_sorted psi_average_norm_contrib_sorted
endif
call save_wavefunction
if(n_det_before == N_det)then
selection_criterion = selection_criterion * 0.5d0
endif
print *, 'N_det = ', N_det
print *, 'N_states = ', N_states
do k = 1, N_states
print*,'State ',k
print *, 'PT2 = ', pt2(k)
print *, 'E = ', CI_energy(k)
print *, 'E(before)+PT2 = ', E_CI_before(k)+pt2(k)
enddo
print *, '-----'
if(N_states.gt.1)then
print*,'Variational Energy difference'
do i = 2, N_states
print*,'Delta E = ',CI_energy(i) - CI_energy(1)
enddo
endif
if(N_states.gt.1)then
print*,'Variational + perturbative Energy difference'
do i = 2, N_states
print*,'Delta E = ',E_CI_before(i)+ pt2(i) - (E_CI_before(1) + pt2(1))
enddo
endif
E_CI_before(1:N_states) = CI_energy(1:N_states)
call ezfio_set_cas_sd_energy(CI_energy(1))
enddo
N_det = min(N_det_max,N_det)
touch N_det psi_det psi_coef
call diagonalize_CI
if(do_pt2_end)then
print*,'Last iteration only to compute the PT2'
threshold_selectors = 1.d0
@ -70,13 +101,12 @@ program full_ci
print *, 'N_det = ', N_det
print *, 'N_states = ', N_states
print *, 'PT2 = ', pt2
print *, 'E = ', CI_energy
print *, 'E+PT2 = ', CI_energy+pt2
print *, 'E = ', CI_energy(1:N_states)
print *, 'E+PT2 = ', CI_energy(1:N_states)+pt2(1:N_states)
print *, '-----'
call ezfio_set_cas_sd_energy_pt2(CI_energy(1)+pt2(1))
endif
integer :: exc_max, degree_min
exc_max = 0
print *, 'CAS determinants : ', N_det_cas
@ -85,6 +115,7 @@ program full_ci
call get_excitation_degree(psi_cas(1,1,k),psi_cas(1,1,i),degree,N_int)
exc_max = max(exc_max,degree)
enddo
print *, psi_coef_cas_diagonalized(i,:)
call debug_det(psi_cas(1,1,i),N_int)
print *, ''
enddo

65
plugins/CAS_SD/cas_sd_selected.irp.f
View File

@ -12,6 +12,7 @@ program full_ci
pt2 = 1.d0
diag_algorithm = "Lapack"
if (N_det > N_det_max) then
call diagonalize_CI
call save_wavefunction
@ -28,32 +29,68 @@ program full_ci
print *, 'E+PT2 = ', CI_energy+pt2
print *, '-----'
endif
double precision :: i_H_psi_array(N_states),diag_H_mat_elem,h,i_O1_psi_array(N_states)
double precision :: E_CI_before(N_states)
if(read_wf)then
call i_H_psi(psi_det(1,1,N_det),psi_det,psi_coef,N_int,N_det,psi_det_size,N_states,i_H_psi_array)
h = diag_H_mat_elem(psi_det(1,1,N_det),N_int)
selection_criterion = dabs(psi_coef(N_det,1) * (i_H_psi_array(1) - h * psi_coef(N_det,1))) * 0.1d0
soft_touch selection_criterion
endif
integer :: n_det_before
print*,'Beginning the selection ...'
E_CI_before(1:N_states) = CI_energy(1:N_states)
do while (N_det < N_det_max.and.maxval(abs(pt2(1:N_st))) > pt2_max)
n_det_before = N_det
call H_apply_CAS_SD_selected(pt2, norm_pert, H_pert_diag, N_st)
PROVIDE psi_coef
PROVIDE psi_det
PROVIDE psi_det_sorted
if (N_det > N_det_max) then
psi_det = psi_det_sorted
psi_coef = psi_coef_sorted
N_det = N_det_max
soft_touch N_det psi_det psi_coef
endif
call diagonalize_CI
if (N_det > N_det_max) then
N_det = N_det_max
psi_det = psi_det_sorted
psi_coef = psi_coef_sorted
touch N_det psi_det psi_coef psi_det_sorted psi_coef_sorted psi_average_norm_contrib_sorted
endif
call save_wavefunction
print *, 'N_det = ', N_det
print *, 'N_states = ', N_states
print *, 'PT2 = ', pt2
print *, 'E = ', CI_energy(1:N_states)
print *, 'E+PT2 = ', CI_energy(1:N_states)+pt2(1:N_states)
if(n_det_before == N_det)then
selection_criterion = selection_criterion * 0.5d0
endif
print *, 'N_det = ', N_det
print *, 'N_states = ', N_states
do k = 1, N_states
print*,'State ',k
print *, 'PT2 = ', pt2(k)
print *, 'E = ', CI_energy(k)
print *, 'E(before)+PT2 = ', E_CI_before(k)+pt2(k)
enddo
print *, '-----'
if(N_states.gt.1)then
print*,'Variational Energy difference'
do i = 2, N_states
print*,'Delta E = ',CI_energy(i) - CI_energy(1)
enddo
endif
if(N_states.gt.1)then
print*,'Variational + perturbative Energy difference'
do i = 2, N_states
print*,'Delta E = ',E_CI_before(i)+ pt2(i) - (E_CI_before(1) + pt2(1))
enddo
endif
E_CI_before(1:N_states) = CI_energy(1:N_states)
call ezfio_set_cas_sd_energy(CI_energy(1))
enddo
call diagonalize_CI
N_det = min(N_det_max,N_det)
touch N_det psi_det psi_coef
call diagonalize_CI
if(do_pt2_end)then
print*,'Last iteration only to compute the PT2'
threshold_selectors = 1.d0
@ -70,7 +107,6 @@ program full_ci
call ezfio_set_cas_sd_energy_pt2(CI_energy(1)+pt2(1))
endif
integer :: exc_max, degree_min
exc_max = 0
print *, 'CAS determinants : ', N_det_cas
@ -79,6 +115,7 @@ program full_ci
call get_excitation_degree(psi_cas(1,1,k),psi_cas(1,1,i),degree,N_int)
exc_max = max(exc_max,degree)
enddo
print *, psi_coef_cas_diagonalized(i,:)
call debug_det(psi_cas(1,1,i),N_int)
print *, ''
enddo

2
plugins/Full_CI/full_ci.irp.f
View File

@ -11,7 +11,7 @@ program full_ci
pt2 = 1.d0
diag_algorithm = "Lapack"
if (N_det > N_det_max) then
call diagonalize_CI
call save_wavefunction

2
src/Davidson/EZFIO.cfg
View File

@ -7,6 +7,6 @@ default: 1.e-12
[n_states_diag]
type: States_number
doc: n_states_diag
default: 1
default: 10
interface: ezfio,provider,ocaml

17
src/Davidson/diagonalization_hs2.irp.f
View File

@ -71,7 +71,7 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sze,N_s
!
! N_st : Number of eigenstates
!
! N_st_diag : Number of states in which H is diagonalized
! N_st_diag : Number of states in which H is diagonalized. Assumed > sze
!
! iunit : Unit for the I/O
!
@ -107,6 +107,9 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sze,N_s
include 'constants.include.F'
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: U, W, R, S, y, h, lambda
if (N_st_diag > sze) then
stop 'error in Davidson : N_st_diag > sze'
endif
PROVIDE nuclear_repulsion
@ -174,16 +177,13 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sze,N_s
converged = .False.
do k=1,N_st_diag
if (k > N_st) then
do k=N_st+1,N_st_diag
do i=1,sze
double precision :: r1, r2
call random_number(r1)
call random_number(r2)
u_in(i,k) = dsqrt(-2.d0*dlog(r1))*dcos(dtwo_pi*r2)
enddo
endif
! Gram-Schmidt
! ------------
@ -247,12 +247,14 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sze,N_s
S(i,k,iter+1) = 0.d0
enddo
enddo
! do k=1,N_st_diag
! do iter2=1,iter
! do l=1,N_st_diag
! do i=1,sze
! U(i,k,iter+1) = U(i,k,iter+1) + U(i,l,iter2)*y(l,iter2,k,1)
! W(i,k,iter+1) = W(i,k,iter+1) + W(i,l,iter2)*y(l,iter2,k,1)
! S(i,k,iter+1) = W(i,k,iter+1) + S(i,l,iter2)*y(l,iter2,k,1)
! enddo
! enddo
! enddo
@ -276,13 +278,16 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sze,N_s
do k=1,N_st_diag
do i=1,sze
R(i,k) = (lambda(k) * U(i,k,iter+1) - W(i,k,iter+1) ) &
* (1.d0 + s2(k) * U(i,k,iter+1) - S(i,k,iter+1) )
* (1.d0 + s2(k) * U(i,k,iter+1) - S(i,k,iter+1) - S_z2_Sz)
enddo
if (k <= N_st) then
residual_norm(k) = u_dot_u(R(1,k),sze)
to_print(1,k) = lambda(k) + nuclear_repulsion
to_print(2,k) = s2(k)
to_print(3,k) = residual_norm(k)
if (residual_norm(k) > 1.e9) then
stop 'Davidson failed'
endif
endif
enddo

4
src/Davidson/diagonalize_CI.irp.f
View File

@ -57,10 +57,10 @@ END_PROVIDER
call davidson_diag_HS2(psi_det,CI_eigenvectors, &
size(CI_eigenvectors,1),CI_electronic_energy, &
N_det,N_states,N_states_diag,N_int,output_determinants)
N_det,min(N_det,N_states),min(N_det,N_states_diag),N_int,output_determinants)
call u_0_S2_u_0(CI_eigenvectors_s2,CI_eigenvectors,N_det,psi_det,N_int,&
N_states_diag,size(CI_eigenvectors,1))
min(N_det,N_states_diag),size(CI_eigenvectors,1))
else if (diag_algorithm == "Lapack") then