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

Barycentric

Browse Source
This commit is contained in:
Anthony Scemama 2018-05-07 15:21:28 +02:00
parent b91c9ebad1
commit 4e55e36260
8 changed files with 48 additions and 61 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
plugins/DavidsonDressed/diagonalize_CI.irp.f
View File

@ -63,8 +63,8 @@ END_PROVIDER
call davidson_diag_HS2(psi_det,CI_eigenvectors_dressed, CI_eigenvectors_s2_dressed,& call davidson_diag_HS2(psi_det,CI_eigenvectors_dressed, CI_eigenvectors_s2_dressed,&
size(CI_eigenvectors_dressed,1), CI_electronic_energy_dressed,& size(CI_eigenvectors_dressed,1), CI_electronic_energy_dressed,&
N_det,min(N_det,N_states),min(N_det,N_states_diag),N_int,1) N_det,min(N_det,N_states),min(N_det,N_states_diag),N_int,1)
call u_0_S2_u_0(CI_eigenvectors_s2_dressed,CI_eigenvectors_dressed,N_det,psi_det,N_int,& ! call u_0_S2_u_0(CI_eigenvectors_s2_dressed,CI_eigenvectors_dressed,N_det,psi_det,N_int,&
N_states_diag,size(CI_eigenvectors_dressed,1)) ! N_states_diag,size(CI_eigenvectors_dressed,1))
else if (diag_algorithm == "Lapack") then else if (diag_algorithm == "Lapack") then

11
plugins/dress_zmq/dress_general.irp.f
View File

@ -29,27 +29,26 @@ subroutine run_dressing(N_st,energy)
E_new = 0.d0 E_new = 0.d0
delta_E = 1.d0 delta_E = 1.d0
iteration = 0 iteration = 0
do while (delta_E > thresh_dress) do iteration=1,n_it_dress_max
iteration += 1
print *, '===============================================' print *, '==============================================='
print *, 'Iteration', iteration, '/', n_it_dress_max print *, 'Iteration', iteration, '/', n_it_dress_max
print *, '===============================================' print *, '==============================================='
print *, '' print *, ''
E_old = dress_e0_denominator(1) !sum(ci_energy_dressed(1:N_states)) E_old = sum(psi_energy(:))
do i=1,N_st do i=1,N_st
call write_double(6,ci_energy_dressed(i),"Energy") call write_double(6,ci_energy_dressed(i),"Energy")
enddo enddo
call diagonalize_ci_dressed call diagonalize_ci_dressed
E_new = dress_e0_denominator(1) !sum(ci_energy_dressed(1:N_states)) E_new = sum(psi_energy(:))
delta_E = (E_new - E_old)/dble(N_states) delta_E = (E_new - E_old)/dble(N_states)
print *, '' print *, ''
call write_double(6,thresh_dress,"thresh_dress") call write_double(6,thresh_dress,"thresh_dress")
call write_double(6,delta_E,"delta_E") call write_double(6,delta_E,"delta_E (undressed)")
delta_E = dabs(delta_E) delta_E = dabs(delta_E)
call save_wavefunction call save_wavefunction
! call ezfio_set_dress_zmq_energy(ci_energy_dressed(1)) ! call ezfio_set_dress_zmq_energy(ci_energy_dressed(1))
if (iteration >= n_it_dress_max) then if (delta_E < thresh_dress) then
exit exit
endif endif
enddo enddo

2
plugins/dress_zmq/dressing.irp.f
View File

@ -78,7 +78,7 @@ BEGIN_PROVIDER [ double precision, delta_ij, (N_states,N_det,2) ]
delta_ij = 0d0 delta_ij = 0d0
E_CI_before(:) = dress_E0_denominator(:) + nuclear_repulsion E_CI_before(:) = psi_energy(:) + nuclear_repulsion
threshold_selectors = 1.d0 threshold_selectors = 1.d0
threshold_generators = 1d0 threshold_generators = 1d0
! if(errr /= 0d0) then ! if(errr /= 0d0) then

6
plugins/dress_zmq/dressing_vector.irp.f
View File

@ -18,12 +18,6 @@
dressing_column_h(j,k) = delta_ij(k,j,1) dressing_column_h(j,k) = delta_ij(k,j,1)
dressing_column_s(j,k) = delta_ij(k,j,2) dressing_column_s(j,k) = delta_ij(k,j,2)
enddo enddo
! tmp = u_dot_v(dressing_column_h(1,k), psi_coef(1,k), N_det) &
! - dressing_column_h(l,k) * psi_coef(l,k)
! dressing_column_h(l,k) -= tmp * f
! tmp = u_dot_v(dressing_column_s(1,k), psi_coef(1,k), N_det) &
! - dressing_column_s(l,k) * psi_coef(l,k)
! dressing_column_s(l,k) -= tmp * f
enddo enddo
END_PROVIDER END_PROVIDER

19
plugins/dress_zmq/energy.irp.f
View File

@ -13,13 +13,24 @@ BEGIN_PROVIDER [ double precision, dress_E0_denominator, (N_states) ]
END_DOC END_DOC
integer :: i integer :: i
if (initialize_dress_E0_denominator) then if (initialize_dress_E0_denominator) then
call u_0_H_u_0(dress_E0_denominator,psi_coef,N_det,psi_det,N_int,N_states,size(psi_coef,1)) if (h0_type == "EN") then
do i=N_det+1,N_states dress_E0_denominator(1:N_states) = psi_energy(1:N_states)
dress_E0_denominator(i) = 0.d0 else if (h0_type == "Barycentric") then
enddo ! dress_E0_denominator(1:N_states) = barycentric_electronic_energy(1:N_states)
dress_E0_denominator(1:N_states) = minval(diagonal_H_matrix_on_psi_det(1:N_det))
else
print *, h0_type, ' not implemented'
stop
endif
! call u_0_H_u_0(dress_E0_denominator,psi_coef,N_det,psi_det,N_int,N_states,size(psi_coef,1))
! do i=N_det+1,N_states
! dress_E0_denominator(i) = 0.d0
! enddo
call write_double(6,dress_E0_denominator(1)+nuclear_repulsion, 'dress Energy denominator') call write_double(6,dress_E0_denominator(1)+nuclear_repulsion, 'dress Energy denominator')
else else
dress_E0_denominator = -huge(1.d0) dress_E0_denominator = -huge(1.d0)
endif endif
END_PROVIDER END_PROVIDER

22
plugins/shiftedbk/shifted_bk_routines.irp.f
View File

@ -4,6 +4,7 @@
&BEGIN_PROVIDER [ double precision, a_s2_i, (N_det, Nproc) ] &BEGIN_PROVIDER [ double precision, a_s2_i, (N_det, Nproc) ]
implicit none implicit none
current_generator_(:) = 0 current_generator_(:) = 0
fock_diag_tmp_(:,:,:) = 0.d0
a_h_i = 0d0 a_h_i = 0d0
a_s2_i = 0d0 a_s2_i = 0d0
END_PROVIDER END_PROVIDER
@ -52,7 +53,7 @@ subroutine dress_with_alpha_buffer(Nstates,Ndet,Nint,delta_ij_loc, i_gen, minili
do i=1,Nstates do i=1,Nstates
de = E0_denominator(i) - haa de = dress_E0_denominator(i) - haa
if(DABS(de) < 1D-5) cycle if(DABS(de) < 1D-5) cycle
c_alpha = a_h_psi(i) / de c_alpha = a_h_psi(i) / de
@ -76,23 +77,4 @@ BEGIN_PROVIDER [ logical, initialize_E0_denominator ]
END_PROVIDER END_PROVIDER
BEGIN_PROVIDER [ double precision, E0_denominator, (N_states) ]
implicit none
BEGIN_DOC
! E0 in the denominator of the PT2
END_DOC
if (initialize_E0_denominator) then
if (h0_type == "EN") then
E0_denominator(1:N_states) = psi_energy(1:N_states)
else if (h0_type == "Barycentric") then
E0_denominator(1:N_states) = barycentric_electronic_energy(1:N_states)
else
print *, h0_type, ' not implemented'
stop
endif
else
E0_denominator = -huge(1.d0)
endif
END_PROVIDER

26
src/Davidson/diagonalization_hs2_dressed.irp.f
View File

@ -35,7 +35,7 @@ subroutine davidson_diag_hs2(dets_in,u_in,s2_out,dim_in,energies,sze,N_st,N_st_d
double precision, intent(inout) :: u_in(dim_in,N_st_diag) double precision, intent(inout) :: u_in(dim_in,N_st_diag)
double precision, intent(out) :: energies(N_st_diag), s2_out(N_st_diag) double precision, intent(out) :: energies(N_st_diag), s2_out(N_st_diag)
integer, intent(in) :: dressing_state integer, intent(in) :: dressing_state
double precision, allocatable :: H_jj(:), S2_jj(:) double precision, allocatable :: H_jj(:)
double precision, external :: diag_H_mat_elem, diag_S_mat_elem double precision, external :: diag_H_mat_elem, diag_S_mat_elem
integer :: i,k integer :: i,k
@ -44,7 +44,7 @@ subroutine davidson_diag_hs2(dets_in,u_in,s2_out,dim_in,energies,sze,N_st,N_st_d
ASSERT (Nint > 0) ASSERT (Nint > 0)
ASSERT (Nint == N_int) ASSERT (Nint == N_int)
PROVIDE mo_bielec_integrals_in_map PROVIDE mo_bielec_integrals_in_map
allocate(H_jj(sze),S2_jj(sze)) allocate(H_jj(sze))
H_jj(1) = diag_h_mat_elem(dets_in(1,1,1),Nint) H_jj(1) = diag_h_mat_elem(dets_in(1,1,1),Nint)
!$OMP PARALLEL DEFAULT(NONE) & !$OMP PARALLEL DEFAULT(NONE) &
@ -66,7 +66,7 @@ subroutine davidson_diag_hs2(dets_in,u_in,s2_out,dim_in,energies,sze,N_st,N_st_d
endif endif
call davidson_diag_hjj_sjj(dets_in,u_in,H_jj,S2_out,energies,dim_in,sze,N_st,N_st_diag,Nint,dressing_state) call davidson_diag_hjj_sjj(dets_in,u_in,H_jj,S2_out,energies,dim_in,sze,N_st,N_st_diag,Nint,dressing_state)
deallocate (H_jj,S2_jj) deallocate (H_jj)
end end
@ -254,9 +254,9 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_
dressing_column_h, size(dressing_column_h,1), s_tmp, size(s_tmp,1), & dressing_column_h, size(dressing_column_h,1), s_tmp, size(s_tmp,1), &
1.d0, W(1,shift+1), size(W,1)) 1.d0, W(1,shift+1), size(W,1))
call dgemm('N','N', sze, N_st_diag, N_st, 0.5d0, & ! call dgemm('N','N', sze, N_st_diag, N_st, 0.5d0, &
dressing_column_s, size(dressing_column_s,1), s_tmp, size(s_tmp,1), & ! dressing_column_s, size(dressing_column_s,1), s_tmp, size(s_tmp,1), &
1.d0, S(1,shift+1), size(S,1)) ! 1.d0, S(1,shift+1), size(S,1))
call dgemm('T','N', N_st, N_st_diag, sze, 1.d0, & call dgemm('T','N', N_st, N_st_diag, sze, 1.d0, &
@ -267,13 +267,13 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_
psi_coef, size(psi_coef,1), s_tmp, size(s_tmp,1), & psi_coef, size(psi_coef,1), s_tmp, size(s_tmp,1), &
1.d0, W(1,shift+1), size(W,1)) 1.d0, W(1,shift+1), size(W,1))
call dgemm('T','N', N_st, N_st_diag, sze, 1.d0, & ! call dgemm('T','N', N_st, N_st_diag, sze, 1.d0, &
dressing_column_s, size(dressing_column_s,1), & ! dressing_column_s, size(dressing_column_s,1), &
U(1,shift+1), size(U,1), 0.d0, s_tmp, size(s_tmp,1)) ! U(1,shift+1), size(U,1), 0.d0, s_tmp, size(s_tmp,1))
!
call dgemm('N','N', sze, N_st_diag, N_st, 0.5d0, & ! call dgemm('N','N', sze, N_st_diag, N_st, 0.5d0, &
psi_coef, size(psi_coef,1), s_tmp, size(s_tmp,1), & ! psi_coef, size(psi_coef,1), s_tmp, size(s_tmp,1), &
1.d0, S(1,shift+1), size(S,1)) ! 1.d0, S(1,shift+1), size(S,1))
endif endif

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

@ -509,11 +509,12 @@ subroutine i_H_j_s2(key_i,key_j,Nint,hij,s2)
select case (degree) select case (degree)
case (2) case (2)
call get_double_excitation(key_i,key_j,exc,phase,Nint) call get_double_excitation(key_i,key_j,exc,phase,Nint)
if (exc(0,1,1) == 1) then
! Mono alpha, mono beta ! Mono alpha, mono beta
if (exc(0,1,1) == 1) then
if ( (exc(1,1,1) == exc(1,2,2)).and.(exc(1,1,2) == exc(1,2,1)) ) then
s2 = -phase
endif
if(exc(1,1,1) == exc(1,2,2) )then if(exc(1,1,1) == exc(1,2,2) )then
if(exc(1,1,2) == exc(1,2,1)) s2 = -phase !!!!!
hij = phase * big_array_exchange_integrals(exc(1,1,1),exc(1,1,2),exc(1,2,1)) hij = phase * big_array_exchange_integrals(exc(1,1,1),exc(1,1,2),exc(1,2,1))
else if (exc(1,2,1) ==exc(1,1,2))then else if (exc(1,2,1) ==exc(1,1,2))then
hij = phase * big_array_exchange_integrals(exc(1,2,1),exc(1,1,1),exc(1,2,2)) hij = phase * big_array_exchange_integrals(exc(1,2,1),exc(1,1,1),exc(1,2,2))
@ -524,8 +525,8 @@ subroutine i_H_j_s2(key_i,key_j,Nint,hij,s2)
exc(1,2,1), & exc(1,2,1), &
exc(1,2,2) ,mo_integrals_map) exc(1,2,2) ,mo_integrals_map)
endif endif
else if (exc(0,1,1) == 2) then
! Double alpha ! Double alpha
else if (exc(0,1,1) == 2) then
hij = phase*(get_mo_bielec_integral( & hij = phase*(get_mo_bielec_integral( &
exc(1,1,1), & exc(1,1,1), &
exc(2,1,1), & exc(2,1,1), &
@ -536,8 +537,8 @@ subroutine i_H_j_s2(key_i,key_j,Nint,hij,s2)
exc(2,1,1), & exc(2,1,1), &
exc(2,2,1), & exc(2,2,1), &
exc(1,2,1) ,mo_integrals_map) ) exc(1,2,1) ,mo_integrals_map) )
else if (exc(0,1,2) == 2) then
! Double beta ! Double beta
else if (exc(0,1,2) == 2) then
hij = phase*(get_mo_bielec_integral( & hij = phase*(get_mo_bielec_integral( &
exc(1,1,2), & exc(1,1,2), &
exc(2,1,2), & exc(2,1,2), &
@ -553,13 +554,13 @@ subroutine i_H_j_s2(key_i,key_j,Nint,hij,s2)
call get_mono_excitation(key_i,key_j,exc,phase,Nint) call get_mono_excitation(key_i,key_j,exc,phase,Nint)
!DIR$ FORCEINLINE !DIR$ FORCEINLINE
call bitstring_to_list_ab(key_i, occ, n_occ_ab, Nint) call bitstring_to_list_ab(key_i, occ, n_occ_ab, Nint)
if (exc(0,1,1) == 1) then
! Mono alpha ! Mono alpha
if (exc(0,1,1) == 1) then
m = exc(1,1,1) m = exc(1,1,1)
p = exc(1,2,1) p = exc(1,2,1)
spin = 1 spin = 1
else
! Mono beta ! Mono beta
else
m = exc(1,1,2) m = exc(1,1,2)
p = exc(1,2,2) p = exc(1,2,2)
spin = 2 spin = 2
@ -567,7 +568,7 @@ subroutine i_H_j_s2(key_i,key_j,Nint,hij,s2)
call get_mono_excitation_from_fock(key_i,key_j,p,m,spin,phase,hij) call get_mono_excitation_from_fock(key_i,key_j,p,m,spin,phase,hij)
case (0) case (0)
print *," ZERO" print *,irp_here,": ZERO"
double precision, external :: diag_S_mat_elem double precision, external :: diag_S_mat_elem
s2 = diag_S_mat_elem(key_i,Nint) s2 = diag_S_mat_elem(key_i,Nint)
hij = diag_H_mat_elem(key_i,Nint) hij = diag_H_mat_elem(key_i,Nint)