2016-06-13 10:29:47 +02:00
|
|
|
program dressed_dmc
|
|
|
|
|
implicit none
|
|
|
|
|
double precision :: E0, hij
|
|
|
|
|
double precision, allocatable :: H_jj(:), energies(:), delta_jj(:), cj(:), hj(:)
|
|
|
|
|
integer :: i
|
|
|
|
|
double precision, external :: diag_h_mat_elem
|
|
|
|
|
|
|
|
|
|
if (.not.read_wf) then
|
|
|
|
|
stop 'read_wf should be true'
|
|
|
|
|
endif
|
|
|
|
|
|
|
|
|
|
PROVIDE mo_bielec_integrals_in_map
|
|
|
|
|
allocate ( H_jj(N_det), delta_jj(N_det), hj(N_det), cj(N_det), energies(N_states) )
|
|
|
|
|
|
|
|
|
|
! Read <i|\Phi_0>
|
|
|
|
|
! -=-=-=-==-=-=-=
|
|
|
|
|
|
|
|
|
|
character*(32) :: w, w2
|
|
|
|
|
integer :: k
|
|
|
|
|
do while (.True.)
|
|
|
|
|
read(*,*) w
|
|
|
|
|
if ( trim(w) == 'Ci_h_psidet' ) then
|
|
|
|
|
exit
|
|
|
|
|
endif
|
|
|
|
|
enddo
|
|
|
|
|
do i=1,N_det
|
|
|
|
|
read(*,*) k, w, hj(i)
|
|
|
|
|
enddo
|
|
|
|
|
|
|
|
|
|
do while (.True.)
|
|
|
|
|
read(*,*) w
|
|
|
|
|
if ( trim(w) == 'Ci_overlap_psidet' ) then
|
|
|
|
|
exit
|
|
|
|
|
endif
|
|
|
|
|
enddo
|
|
|
|
|
do i=1,N_det
|
|
|
|
|
read(*,*) k, w, cj(i)
|
|
|
|
|
enddo
|
|
|
|
|
|
|
|
|
|
read(*,*)
|
|
|
|
|
read(*,*) w, w2, E0
|
|
|
|
|
print *, 'E0=', E0
|
|
|
|
|
print *, 'Ndet = ', N_det
|
|
|
|
|
|
|
|
|
|
! Compute delta_ii
|
|
|
|
|
! -=-=-=-==-=-=-=-
|
|
|
|
|
|
|
|
|
|
do i=1,N_det
|
|
|
|
|
call i_H_psi(psi_det(1,1,i),psi_det,cj,N_int,N_det,size(psi_coef,1),N_states,energies)
|
|
|
|
|
if (dabs(cj(i)) < 1.d-6) then
|
|
|
|
|
delta_jj(i) = 0.d0
|
|
|
|
|
else
|
|
|
|
|
delta_jj(i) = (hj(i) - energies(1))/cj(i)
|
|
|
|
|
endif
|
|
|
|
|
H_jj(i) = diag_h_mat_elem(psi_det(1,1,i),N_int) + delta_jj(i)
|
|
|
|
|
print *, 'Delta_jj(',i,') = ', Delta_jj(i), H_jj(i)
|
|
|
|
|
enddo
|
|
|
|
|
|
|
|
|
|
|
2016-10-07 19:21:04 +02:00
|
|
|
call davidson_diag_hjj(psi_det,psi_coef,H_jj,energies,size(psi_coef,1),N_det,N_states,N_states_diag,N_int,6)
|
2016-06-13 10:29:47 +02:00
|
|
|
|
|
|
|
|
call save_wavefunction
|
|
|
|
|
call write_spindeterminants
|
|
|
|
|
|
|
|
|
|
E0 = 0.d0
|
|
|
|
|
do i=1,N_det
|
|
|
|
|
call i_H_psi(psi_det(1,1,i),psi_det,psi_coef(1,1),N_int,N_det,size(psi_coef,1),N_states,energies)
|
|
|
|
|
E0 += psi_coef(i,1) * energies(1)
|
|
|
|
|
enddo
|
|
|
|
|
print *, 'Trial energy: ', E0 + nuclear_repulsion
|
|
|
|
|
|
|
|
|
|
deallocate (H_jj, delta_jj, energies, cj)
|
|
|
|
|
end
|
