-
Notifications
You must be signed in to change notification settings - Fork 2
/
target_pt2_qmc.irp.f
121 lines (108 loc) · 3.79 KB
/
target_pt2_qmc.irp.f
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
program e_curve
use bitmasks
implicit none
integer :: i,j,k, nab, m, l, n_up, n_dn, n
double precision :: norm, E, hij, num, ci, cj
integer, allocatable :: iorder(:)
double precision , allocatable :: norm_sort(:), psi_bilinear_matrix_values_save(:)
nab = n_det_alpha_unique+n_det_beta_unique
allocate ( norm_sort(0:nab), iorder(0:nab), psi_bilinear_matrix_values_save(N_det) )
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)
if (.not.read_wf) then
stop 'Please set read_wf to true'
endif
psi_bilinear_matrix_values_save = psi_bilinear_matrix_values(:,1)
print *, '=========================================================='
print '(A8,2X,A8,2X,A12,2X,A10,2X,A12)', 'Thresh.', 'Ndet', 'Cost', 'Norm', 'E'
print *, '=========================================================='
integer(bit_kind), allocatable :: det_i(:,:), det_j(:,:)
double precision :: thresh, E_min, E_max, E_prev
thresh = 0.d0
call compute_energy(psi_bilinear_matrix_values_save,E_max,m,norm)
call i_h_j(psi_det_sorted(1,1,1), psi_det_sorted(1,1,1), N_int, E_min)
print *, E_min, E_max
n_up = nab
n_dn = 0
do while (n_up > n_dn)
n = n_dn + (n_up-n_dn)/2
psi_bilinear_matrix_values_save = psi_bilinear_matrix_values(:,1)
do j=1,n
i = iorder(j)
if (i<0) then
do k=1,n_det
if (psi_bilinear_matrix_columns(k) == -i) then
psi_bilinear_matrix_values_save(k) = 0.d0
endif
enddo
else
do k=1,n_det
if (psi_bilinear_matrix_rows(k) == i) then
psi_bilinear_matrix_values_save(k) = 0.d0
endif
enddo
endif
enddo
call compute_energy(psi_bilinear_matrix_values_save,E,m,norm)
print '(E9.1,2X,I8,2X,F10.2,2X,F10.6,2X,F12.6)', norm_sort(n), m, &
dble( elec_alpha_num**3 + elec_alpha_num**2 * m ) / &
dble( elec_alpha_num**3 + elec_alpha_num**2 * n ), norm, E
if (E < target_energy) then
n_dn = n+1
else
n_up = n
endif
enddo
print *, '=========================================================='
print *, norm_sort(n), target_energy
deallocate (iorder, norm_sort, psi_bilinear_matrix_values_save)
end
subroutine compute_energy(psi_bilinear_matrix_values_save, E, m, norm)
implicit none
BEGIN_DOC
! Compute an energy when a threshold is applied
END_DOC
double precision, intent(in) :: psi_bilinear_matrix_values_save(n_det)
integer(bit_kind), allocatable :: det_i(:,:), det_j(:,:)
integer :: i,j, k, l, m
double precision :: num, norm, ci, cj, hij, E
num = 0.d0
norm = 0.d0
m = 0
!$OMP PARALLEL DEFAULT(SHARED) PRIVATE(k,l,det_i,det_j,ci,cj,hij) REDUCTION(+:norm,m,num)
allocate( det_i(N_int,2), det_j(N_int,2))
!$OMP DO schedule(guided)
do k=1,n_det
if (psi_bilinear_matrix_values_save(k) == 0.d0) then
cycle
endif
ci = psi_bilinear_matrix_values_save(k)
det_i(:,1) = psi_det_alpha_unique(:,psi_bilinear_matrix_rows(k))
det_i(:,2) = psi_det_beta_unique(:,psi_bilinear_matrix_columns(k))
do l=1,n_det
if (psi_bilinear_matrix_values_save(l) == 0.d0) then
cycle
endif
cj = psi_bilinear_matrix_values_save(l)
det_j(:,1) = psi_det_alpha_unique(:,psi_bilinear_matrix_rows(l))
det_j(:,2) = psi_det_beta_unique(:,psi_bilinear_matrix_columns(l))
call i_h_j(det_i, det_j, N_int, hij)
num = num + ci*cj*hij
enddo
norm = norm + ci*ci
m = m+1
enddo
!$OMP END DO
deallocate (det_i,det_j)
!$OMP END PARALLEL
E = num / norm + nuclear_repulsion
end