10
0
mirror of https://github.com/LCPQ/quantum_package synced 2024-12-27 06:43:48 +01:00
quantum_package/plugins/QMC/truncate_wf_spin.irp.f

118 lines
3.5 KiB
Fortran
Raw Normal View History

2017-07-14 16:15:34 +02:00
program truncate
read_wf = .True.
SOFT_TOUCH read_wf
call run
end
subroutine run
2017-07-12 23:44:37 +02:00
use bitmasks
implicit none
integer :: i,j,k, kk, nab, m, l
double precision :: norm, E, hij, num, ci, cj
integer, allocatable :: iorder(:)
double precision , allocatable :: norm_sort(:)
double precision :: e_0(N_states)
PROVIDE mo_bielec_integrals_in_map H_apply_buffer_allocated
nab = n_det_alpha_unique+n_det_beta_unique
allocate ( norm_sort(0:nab), iorder(0:nab) )
integer(bit_kind), allocatable :: det_i(:,:), det_j(:,:)
double precision, allocatable :: u_t(:,:), v_t(:,:), s_t(:,:)
double precision, allocatable :: u_0(:,:), v_0(:,:)
allocate(u_t(N_states,N_det),v_t(N_states,N_det),s_t(N_states,N_det))
2017-07-14 03:12:46 +02:00
allocate(u_0(N_det,N_states),v_0(N_det,N_states))
2017-07-12 23:44:37 +02:00
norm_sort(0) = 0.d0
iorder(0) = 0
do i=1,n_det_alpha_unique
norm_sort(i) = det_alpha_norm(i)
iorder(i) = i
enddo
do i=1,n_det_beta_unique
norm_sort(i+n_det_alpha_unique) = det_beta_norm(i)
iorder(i+n_det_alpha_unique) = -i
enddo
call dsort(norm_sort(1),iorder(1),nab)
2017-10-09 14:14:27 +02:00
PROVIDE psi_bilinear_matrix_values psi_bilinear_matrix_rows psi_bilinear_matrix_columns
PROVIDE nuclear_repulsion
2017-07-12 23:44:37 +02:00
print *, ''
do j=0,nab
i = iorder(j)
if (i<0) then
2017-07-14 03:24:51 +02:00
!$OMP PARALLEL DO PRIVATE(k)
2017-07-12 23:44:37 +02:00
do k=1,n_det
if (psi_bilinear_matrix_columns(k) == -i) then
2017-10-09 14:14:27 +02:00
do l=1,N_states
psi_bilinear_matrix_values(k,l) = 0.d0
enddo
2017-07-12 23:44:37 +02:00
endif
enddo
2017-07-14 03:24:51 +02:00
!$OMP END PARALLEL DO
2017-07-12 23:44:37 +02:00
else
2017-07-14 03:24:51 +02:00
!$OMP PARALLEL DO PRIVATE(k)
2017-07-12 23:44:37 +02:00
do k=1,n_det
if (psi_bilinear_matrix_rows(k) == i) then
2017-10-09 14:14:27 +02:00
do l=1,N_states
psi_bilinear_matrix_values(k,l) = 0.d0
enddo
2017-07-12 23:44:37 +02:00
endif
enddo
2017-07-14 03:24:51 +02:00
!$OMP END PARALLEL DO
2017-07-12 23:44:37 +02:00
endif
2017-07-14 03:24:51 +02:00
if (ci_threshold > norm_sort(j)) then
2017-07-12 23:44:37 +02:00
cycle
endif
2017-10-09 14:14:27 +02:00
u_0(1:N_det,1:N_states) = psi_bilinear_matrix_values(1:N_det,1:N_states)
v_0(1:N_det,1:N_states) = 0.d0
u_t(1:N_states,1:N_det) = 0.d0
v_t(1:N_states,1:N_det) = 0.d0
s_t(1:N_states,1:N_det) = 0.d0
2017-07-12 23:44:37 +02:00
call dtranspose( &
u_0, &
size(u_0, 1), &
u_t, &
size(u_t, 1), &
N_det, N_states)
2017-07-14 03:24:51 +02:00
print *, 'Computing H|Psi> ...'
2017-07-12 23:44:37 +02:00
call H_S2_u_0_nstates_openmp_work(v_t,s_t,u_t,N_states,N_det,1,N_det,0,1)
2017-07-14 03:24:51 +02:00
print *, 'Done'
2017-07-12 23:44:37 +02:00
call dtranspose( &
v_t, &
size(v_t, 1), &
v_0, &
size(v_0, 1), &
N_states, N_det)
double precision, external :: u_dot_u, u_dot_v
do i=1,N_states
2017-10-09 14:14:27 +02:00
e_0(i) = u_dot_v(u_0(1,i),v_0(1,i),N_det)/u_dot_u(u_0(1,i),N_det)
print *, 'E = ', e_0(i) + nuclear_repulsion
2017-07-12 23:44:37 +02:00
enddo
m = 0
do k=1,n_det
2017-10-06 15:41:44 +02:00
if (sum(psi_bilinear_matrix_values(k,1:N_states)) /= 0.d0) then
2017-07-12 23:44:37 +02:00
m = m+1
endif
enddo
2017-10-06 15:41:44 +02:00
do k=1,N_states
E = E_0(k) + nuclear_repulsion
enddo
2017-07-12 23:44:37 +02:00
print *, 'Number of determinants:', m
exit
enddo
2017-07-14 21:31:57 +02:00
call wf_of_psi_bilinear_matrix(.True.)
2017-07-12 23:44:37 +02:00
call save_wavefunction
deallocate (iorder, norm_sort)
end