LCPQ/quantum_package
10
0
Fork 0
mirror of https://github.com/LCPQ/quantum_package synced 2025-01-05 19:08:52 +01:00
Code Issues Releases Wiki Activity

ready for integral driven version - dirty copypastas

Browse Source
This commit is contained in:
Yann Garniron 2016-07-25 14:10:28 +02:00
parent 98ecbb8a25
commit a117ac1a6b
2 changed files with 531 additions and 414 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

25
plugins/Full_CI_ZMQ/fci_zmq.irp.f
View File

@ -7,8 +7,8 @@ program fci_zmq
double precision, allocatable :: pt2(:), norm_pert(:), H_pert_diag(:) double precision, allocatable :: pt2(:), norm_pert(:), H_pert_diag(:)
integer :: N_st, degree integer :: N_st, degree
integer :: it, mit(0:5) integer :: it, mit(0:6)
mit = (/1, 268, 1517, 10018, 45096, 100000/) mit = (/1, 246, 1600, 17528, 112067, 519459, 2685970/)
it = 0 it = 0
N_st = N_states N_st = N_states
allocate (pt2(N_st), norm_pert(N_st),H_pert_diag(N_st)) allocate (pt2(N_st), norm_pert(N_st),H_pert_diag(N_st))
@ -41,18 +41,18 @@ program fci_zmq
E_CI_before = CI_energy E_CI_before = CI_energy
do while (N_det < N_det_max.and.maxval(abs(pt2(1:N_st))) > pt2_max) do while (N_det < N_det_max.and.maxval(abs(pt2(1:N_st))) > pt2_max)
n_det_before = N_det n_det_before = N_det
! call H_apply_FCI(pt2, norm_pert, H_pert_diag, N_st) ! call H_apply_FCI(pt2, norm_pert, H_pert_diag, N_st)
it += 1 it += 1
if(it > 5) stop if(it > 6) stop
call ZMQ_selection(mit(it) - mit(it-1), pt2) !(max(N_det*3, 1000-N_det)) call ZMQ_selection(mit(it) - mit(it-1), pt2) ! max(1000-N_det, N_det), pt2)
do i=1, N_det !do i=1, N_det
!if(popcnt(psi_det(1,1,i)) + popcnt(psi_det(2,1,i)) /= 23) stop "ZZ1" -2099.2504682049275 !if(popcnt(psi_det(1,1,i)) + popcnt(psi_det(2,1,i)) /= 23) stop "ZZ1" -2099.2504682049275
!if(popcnt(psi_det(1,2,i)) + popcnt(psi_det(2,2,i)) /= 23) stop "ZZ2" !if(popcnt(psi_det(1,2,i)) + popcnt(psi_det(2,2,i)) /= 23) stop "ZZ2"
do k=1,i-1 ! do k=1,i-1
if(detEq(psi_det(1,1,i), psi_det(1,1,k), N_int)) stop "ATRRGRZER" ! if(detEq(psi_det(1,1,i), psi_det(1,1,k), N_int)) stop "ATRRGRZER"
end do ! end do
end do !end do
PROVIDE psi_coef PROVIDE psi_coef
PROVIDE psi_det PROVIDE psi_det
PROVIDE psi_det_sorted PROVIDE psi_det_sorted
@ -128,11 +128,10 @@ subroutine ZMQ_selection(N, pt2)
integer :: i integer :: i
integer, external :: omp_get_thread_num integer, external :: omp_get_thread_num
double precision, intent(out) :: pt2(N_states) double precision, intent(out) :: pt2(N_states)
call flip_generators() !call flip_generators()
call new_parallel_job(zmq_to_qp_run_socket,'selection') call new_parallel_job(zmq_to_qp_run_socket,'selection')
call create_selection_buffer(N, N*2, b) call create_selection_buffer(N, N*2, b)
do i= N_det_generators, 1, -1 do i= N_det_generators, 1, -1
write(task,*) i, N write(task,*) i, N
call add_task_to_taskserver(zmq_to_qp_run_socket,task) call add_task_to_taskserver(zmq_to_qp_run_socket,task)
@ -149,7 +148,7 @@ subroutine ZMQ_selection(N, pt2)
endif endif
!$OMP END PARALLEL !$OMP END PARALLEL
call end_parallel_job(zmq_to_qp_run_socket, 'selection') call end_parallel_job(zmq_to_qp_run_socket, 'selection')
call flip_generators() !call flip_generators()
call fill_H_apply_buffer_no_selection(b%cur,b%det,N_int,0) !!! PAS DE ROBIN call fill_H_apply_buffer_no_selection(b%cur,b%det,N_int,0) !!! PAS DE ROBIN
call copy_H_apply_buffer_to_wf() call copy_H_apply_buffer_to_wf()
end subroutine end subroutine

580
plugins/Full_CI_ZMQ/selection.irp.f
View File

@ -28,6 +28,7 @@ subroutine selection_slave(thread,iproc)
buf%N = 0 buf%N = 0
ctask = 1 ctask = 1
pt2 = 0d0 pt2 = 0d0
do do
call get_task_from_taskserver(zmq_to_qp_run_socket,worker_id, task_id(ctask), task) call get_task_from_taskserver(zmq_to_qp_run_socket,worker_id, task_id(ctask), task)
done = task_id(ctask) == 0 done = task_id(ctask) == 0
@ -45,26 +46,26 @@ subroutine selection_slave(thread,iproc)
if(done) ctask = ctask - 1 if(done) ctask = ctask - 1
if(done .or. ctask == size(task_id)) then if(done .or. ctask == size(task_id)) then
if(ctask > 0 .and. buf%N /= 0) then if(buf%N == 0 .and. ctask > 0) stop "uninitialized selection_buffer"
do i=1, ctask
call task_done_to_taskserver(zmq_to_qp_run_socket,worker_id,task_id(i))
end do
if(ctask > 0) then
call push_selection_results(zmq_socket_push, pt2, buf, task_id(1), ctask) call push_selection_results(zmq_socket_push, pt2, buf, task_id(1), ctask)
pt2 = 0d0 pt2 = 0d0
buf%cur = 0 buf%cur = 0
end if end if
do i=1, ctask
call task_done_to_taskserver(zmq_to_qp_run_socket,worker_id,task_id(i))
end do
ctask = 0 ctask = 0
end if end if
if(done) exit if(done) exit
ctask = ctask + 1 ctask = ctask + 1
end do end do
call disconnect_from_taskserver(zmq_to_qp_run_socket,zmq_socket_push,worker_id) call disconnect_from_taskserver(zmq_to_qp_run_socket,zmq_socket_push,worker_id)
call end_zmq_to_qp_run_socket(zmq_to_qp_run_socket) call end_zmq_to_qp_run_socket(zmq_to_qp_run_socket)
call end_zmq_push_socket(zmq_socket_push,thread) call end_zmq_push_socket(zmq_socket_push,thread)
end subroutine end subroutine
@ -82,16 +83,21 @@ subroutine push_selection_results(zmq_socket_push, pt2, b, task_id, ntask)
call sort_selection_buffer(b) call sort_selection_buffer(b)
rc = f77_zmq_send( zmq_socket_push, b%cur, 4, ZMQ_SNDMORE) rc = f77_zmq_send( zmq_socket_push, b%cur, 4, ZMQ_SNDMORE)
if(rc /= 4) stop "push"
rc = f77_zmq_send( zmq_socket_push, pt2, 8*N_states, ZMQ_SNDMORE) rc = f77_zmq_send( zmq_socket_push, pt2, 8*N_states, ZMQ_SNDMORE)
if(rc /= 8*N_states) stop "push"
rc = f77_zmq_send( zmq_socket_push, b%val(1), 8*b%cur, ZMQ_SNDMORE) rc = f77_zmq_send( zmq_socket_push, b%val(1), 8*b%cur, ZMQ_SNDMORE)
if(rc /= 8*b%cur) stop "push"
rc = f77_zmq_send( zmq_socket_push, b%det(1,1,1), bit_kind*N_int*2*b%cur, ZMQ_SNDMORE) rc = f77_zmq_send( zmq_socket_push, b%det(1,1,1), bit_kind*N_int*2*b%cur, ZMQ_SNDMORE)
if(rc /= bit_kind*N_int*2*b%cur) stop "push"
rc = f77_zmq_send( zmq_socket_push, ntask, 4, ZMQ_SNDMORE) rc = f77_zmq_send( zmq_socket_push, ntask, 4, ZMQ_SNDMORE)
if(rc /= 4) stop "push"
rc = f77_zmq_send( zmq_socket_push, task_id(1), ntask*4, 0) rc = f77_zmq_send( zmq_socket_push, task_id(1), ntask*4, 0)
if(rc /= 4*ntask) stop "push"
end subroutine end subroutine
@ -108,16 +114,22 @@ subroutine pull_selection_results(zmq_socket_pull, pt2, val, det, N, task_id, nt
integer :: rc, rn, i integer :: rc, rn, i
rc = f77_zmq_recv( zmq_socket_pull, N, 4, ZMQ_SNDMORE) rc = f77_zmq_recv( zmq_socket_pull, N, 4, ZMQ_SNDMORE)
if(rc /= 4) stop "pull"
rc = f77_zmq_recv( zmq_socket_pull, pt2, N_states*8, ZMQ_SNDMORE) rc = f77_zmq_recv( zmq_socket_pull, pt2, N_states*8, ZMQ_SNDMORE)
if(rc /= 8*N_states) stop "pull"
rc = f77_zmq_recv( zmq_socket_pull, val(1), 8*N, ZMQ_SNDMORE) rc = f77_zmq_recv( zmq_socket_pull, val(1), 8*N, ZMQ_SNDMORE)
if(rc /= 8*N) stop "pull"
rc = f77_zmq_recv( zmq_socket_pull, det(1,1,1), bit_kind*N_int*2*N, ZMQ_SNDMORE) rc = f77_zmq_recv( zmq_socket_pull, det(1,1,1), bit_kind*N_int*2*N, ZMQ_SNDMORE)
if(rc /= bit_kind*N_int*2*N) stop "pull"
rc = f77_zmq_recv( zmq_socket_pull, ntask, 4, ZMQ_SNDMORE) rc = f77_zmq_recv( zmq_socket_pull, ntask, 4, ZMQ_SNDMORE)
if(rc /= 4) stop "pull"
rc = f77_zmq_recv( zmq_socket_pull, task_id(1), ntask*4, 0) rc = f77_zmq_recv( zmq_socket_pull, task_id(1), ntask*4, 0)
if(rc /= 4*ntask) stop "pull"
end subroutine end subroutine
@ -247,12 +259,15 @@ subroutine selection_collector(b, pt2)
double precision, allocatable :: val(:) double precision, allocatable :: val(:)
integer(bit_kind), allocatable :: det(:,:,:) integer(bit_kind), allocatable :: det(:,:,:)
integer, allocatable :: task_id(:) integer, allocatable :: task_id(:)
integer :: done
real :: time, time0
zmq_to_qp_run_socket = new_zmq_to_qp_run_socket() zmq_to_qp_run_socket = new_zmq_to_qp_run_socket()
zmq_socket_pull = new_zmq_pull_socket() zmq_socket_pull = new_zmq_pull_socket()
allocate(val(b%N), det(N_int, 2, b%N), task_id(N_det)) allocate(val(b%N), det(N_int, 2, b%N), task_id(N_det))
pt2 = 0d0 done = 0
more = 1 more = 1
pt2(:) = 0d0
call CPU_TIME(time0)
do while (more == 1) do while (more == 1)
call pull_selection_results(zmq_socket_pull, pt2_mwen, val(1), det(1,1,1), N, task_id, ntask) call pull_selection_results(zmq_socket_pull, pt2_mwen, val(1), det(1,1,1), N, task_id, ntask)
pt2 += pt2_mwen pt2 += pt2_mwen
@ -261,12 +276,15 @@ subroutine selection_collector(b, pt2)
end do end do
do i=1, ntask do i=1, ntask
if (task_id(i) /= 0) then if(task_id(i) == 0) stop "collector"
call zmq_delete_task(zmq_to_qp_run_socket,zmq_socket_pull,task_id(i),more) call zmq_delete_task(zmq_to_qp_run_socket,zmq_socket_pull,task_id(i),more)
endif
end do end do
done += ntask
call CPU_TIME(time)
print *, "DONE" , done, time - time0
end do end do
call end_zmq_to_qp_run_socket(zmq_to_qp_run_socket) call end_zmq_to_qp_run_socket(zmq_to_qp_run_socket)
call end_zmq_pull_socket(zmq_socket_pull) call end_zmq_pull_socket(zmq_socket_pull)
call sort_selection_buffer(b) call sort_selection_buffer(b)
@ -328,11 +346,6 @@ subroutine select_singles(i_generator,hole_mask,particle_mask,fock_diag_tmp,E0,p
enddo enddo
! Create the mini wave function where <i|H|psi_mini> = <i|H|psi>
! --------------------------------------------------------------
! integer(bit_kind) :: psi_det_connected(N_int,2,psi_selectors_size)
! double precision :: psi_coef_connected(psi_selectors_size,N_states)
integer :: ptr_microlist(0:mo_tot_num * 2 + 1), N_microlist(0:mo_tot_num * 2) integer :: ptr_microlist(0:mo_tot_num * 2 + 1), N_microlist(0:mo_tot_num * 2)
integer, allocatable :: idx_microlist(:) integer, allocatable :: idx_microlist(:)
@ -342,9 +355,6 @@ subroutine select_singles(i_generator,hole_mask,particle_mask,fock_diag_tmp,E0,p
allocate(microlist(N_int, 2, N_det_selectors * 3), psi_coef_microlist(psi_selectors_size * 3, N_states), idx_microlist(N_det_selectors * 3)) allocate(microlist(N_int, 2, N_det_selectors * 3), psi_coef_microlist(psi_selectors_size * 3, N_states), idx_microlist(N_det_selectors * 3))
do ispin=1,2 do ispin=1,2
! do k=1,N_int
! ion_det(k,ispin) = psi_det_generators(k,ispin,i_generator)
! enddo
do i=1, N_holes(ispin) do i=1, N_holes(ispin)
@ -356,7 +366,6 @@ subroutine select_singles(i_generator,hole_mask,particle_mask,fock_diag_tmp,E0,p
integer :: j_hole, k_hole integer :: j_hole, k_hole
k_hole = ishft(i_hole-1,-bit_kind_shift)+1 ! N_int k_hole = ishft(i_hole-1,-bit_kind_shift)+1 ! N_int
j_hole = i_hole-ishft(k_hole-1,bit_kind_shift)-1 ! bit index j_hole = i_hole-ishft(k_hole-1,bit_kind_shift)-1 ! bit index
! ion_det(k_hole,ispin) = ibclr(psi_det_generators(k_hole,ispin,i_generator),j_hole)
ion_det(k_hole,ispin) = ibclr(ion_det(k_hole,ispin),j_hole) ion_det(k_hole,ispin) = ibclr(ion_det(k_hole,ispin),j_hole)
@ -370,49 +379,36 @@ subroutine select_singles(i_generator,hole_mask,particle_mask,fock_diag_tmp,E0,p
cycle cycle
endif endif
! Create particles
! ----------------
do j=1,N_particles(ispin) do j=1,N_particles(ispin)
! exc_det(k_hole,ispin) = ion_det(k_hole,ispin)
exc_det(:,:) = ion_det(:,:) exc_det(:,:) = ion_det(:,:)
integer :: i_particle integer :: i_particle
i_particle = particle_list(j,ispin) i_particle = particle_list(j,ispin)
! Apply the particle
integer :: j_particle, k_particle integer :: j_particle, k_particle
k_particle = ishft(i_particle-1,-bit_kind_shift)+1 ! N_int k_particle = ishft(i_particle-1,-bit_kind_shift)+1 ! N_int
j_particle = i_particle-ishft(k_particle-1,bit_kind_shift)-1 ! bit index j_particle = i_particle-ishft(k_particle-1,bit_kind_shift)-1 ! bit index
! exc_det(k_particle,ispin) = ibset(ion_det(k_particle,ispin),j_particle)
exc_det(k_particle,ispin) = ibset(exc_det(k_particle,ispin),j_particle) exc_det(k_particle,ispin) = ibset(exc_det(k_particle,ispin),j_particle)
! TODO
logical, external :: is_in_wavefunction logical, external :: is_in_wavefunction
logical :: nok logical :: nok
! TODO : Check connected to ref
if (.not. is_in_wavefunction(exc_det,N_int)) then if (.not. is_in_wavefunction(exc_det,N_int)) then
! Compute perturbative contribution and select determinant
double precision :: i_H_psi_value(N_states), i_H_psi_value2(N_states) double precision :: i_H_psi_value(N_states), i_H_psi_value2(N_states)
double precision :: i_H_full(N_states)
i_H_psi_value = 0d0 i_H_psi_value = 0d0
i_H_psi_value2 = 0d0 i_H_psi_value2 = 0d0
i_H_full = 0d0
integer :: sporb integer :: sporb
! call i_H_psi(exc_det,psi_selectors,psi_selectors_coef,N_int,N_det_selectors,psi_selectors_size,N_states,i_H_full)
!
nok = .false. nok = .false.
sporb = i_particle + (ispin - 1) * mo_tot_num sporb = i_particle + (ispin - 1) * mo_tot_num
! ! ! subroutine check_past(det, list, idx, N, cur, ok, Nint)
if(N_microlist(sporb) > 0) call check_past(exc_det, microlist(1,1,ptr_microlist(sporb)), idx_microlist(ptr_microlist(sporb)), N_microlist(sporb), i_generator,nok, N_int) if(N_microlist(sporb) > 0) call check_past(exc_det, microlist(1,1,ptr_microlist(sporb)), idx_microlist(ptr_microlist(sporb)), N_microlist(sporb), i_generator, nok, N_int)
if(nok) cycle if(nok) cycle
!
if(N_microlist(0) > 0) call i_H_psi(exc_det,microlist,psi_coef_microlist,N_int,N_microlist(0),psi_selectors_size*3,N_states,i_H_psi_value) if(N_microlist(0) > 0) call i_H_psi(exc_det,microlist,psi_coef_microlist,N_int,N_microlist(0),psi_selectors_size*3,N_states,i_H_psi_value)
if(N_microlist(sporb) > 0) call i_H_psi(exc_det,microlist(1,1,ptr_microlist(sporb)),psi_coef_microlist(ptr_microlist(sporb), 1),N_int,N_microlist(sporb),psi_selectors_size*3,N_states,i_H_psi_value2) if(N_microlist(sporb) > 0) call i_H_psi(exc_det,microlist(1,1,ptr_microlist(sporb)),psi_coef_microlist(ptr_microlist(sporb), 1),N_int,N_microlist(sporb),psi_selectors_size*3,N_states,i_H_psi_value2)
i_H_psi_value(:) = i_H_psi_value(:) + i_H_psi_value2(:) i_H_psi_value(:) = i_H_psi_value(:) + i_H_psi_value2(:)
double precision :: Hii, diag_H_mat_elem_fock double precision :: Hii, diag_H_mat_elem_fock
Hii = diag_H_mat_elem_fock(psi_det_generators(1,1,i_generator),exc_det,fock_diag_tmp,N_int) Hii = diag_H_mat_elem_fock(psi_det_generators(1,1,i_generator),exc_det,fock_diag_tmp,N_int)
@ -422,9 +418,6 @@ subroutine select_singles(i_generator,hole_mask,particle_mask,fock_diag_tmp,E0,p
e_pertm = 0d0 e_pertm = 0d0
do k=1,N_states do k=1,N_states
! if(dabs(1d0 - i_H_psi_value(k)/i_H_full(k)) > 1d-6) then
! stop "PAS BON, PAS BOOOOON!! (single)"
! endif
if (i_H_psi_value(k) == 0.d0) cycle if (i_H_psi_value(k) == 0.d0) cycle
delta_E = E0(k) - Hii delta_E = E0(k) - Hii
if (delta_E < 0.d0) then if (delta_E < 0.d0) then
@ -466,10 +459,6 @@ subroutine select_doubles(i_generator,hole_mask,particle_mask,fock_diag_tmp,E0,p
type(selection_buffer), intent(inout) :: buf type(selection_buffer), intent(inout) :: buf
logical :: isinwf(mo_tot_num*2, mo_tot_num*2) logical :: isinwf(mo_tot_num*2, mo_tot_num*2)
double precision :: d0s(mo_tot_num, mo_tot_num, N_states) double precision :: d0s(mo_tot_num, mo_tot_num, N_states)
d0s = 0d0
! double precision, save :: d0 = 0d0
! double precision, save :: d1 = 0d0
! double precision, save :: d2 = 0d0
integer :: i,j,k,l,j1,j2,i1,i2,ib,jb integer :: i,j,k,l,j1,j2,i1,i2,ib,jb
@ -536,7 +525,7 @@ subroutine select_doubles(i_generator,hole_mask,particle_mask,fock_diag_tmp,E0,p
if(ispin1 == ispin2) ib = i1+1 if(ispin1 == ispin2) ib = i1+1
do i2=ib, N_holes(ispin2) do i2=ib, N_holes(ispin2)
ion_det(:,:) = psi_det_generators(:,:,i_generator) ion_det(:,:) = psi_det_generators(:,:,i_generator)
! call set_hole(ion_det, hole_list(i1,ispin1), ispin1, hole_list(i1,ispin1), ispin1, Nint)
i_hole1 = hole_list(i1,ispin1) i_hole1 = hole_list(i1,ispin1)
k_hole = ishft(i_hole1-1,-bit_kind_shift)+1 ! N_int k_hole = ishft(i_hole1-1,-bit_kind_shift)+1 ! N_int
j_hole = i_hole1-ishft(k_hole-1,bit_kind_shift)-1 ! bit index j_hole = i_hole1-ishft(k_hole-1,bit_kind_shift)-1 ! bit index
@ -551,7 +540,7 @@ subroutine select_doubles(i_generator,hole_mask,particle_mask,fock_diag_tmp,E0,p
microlist, idx_microlist, N_microlist, ptr_microlist, & microlist, idx_microlist, N_microlist, ptr_microlist, &
tmicrolist, idx_tmicrolist, N_tmicrolist, ptr_tmicrolist, & tmicrolist, idx_tmicrolist, N_tmicrolist, ptr_tmicrolist, &
isinwf, d0s, N_int) isinwf, d0s, N_int)
if(N_microlist(0) > 0 .and. idx_microlist(1) > i_generator) stop "wtf..."
if(ptr_microlist(mo_tot_num * 2 + 1) == 1 .and. ptr_tmicrolist(mo_tot_num * 2 + 1) == 1) cycle if(ptr_microlist(mo_tot_num * 2 + 1) == 1 .and. ptr_tmicrolist(mo_tot_num * 2 + 1) == 1) cycle
call finish_isinwf(ion_det, psi_det_sorted(1,1,N_det_selectors+1), N_det - N_det_selectors, isinwf) call finish_isinwf(ion_det, psi_det_sorted(1,1,N_det_selectors+1), N_det - N_det_selectors, isinwf)
@ -577,7 +566,7 @@ subroutine select_doubles(i_generator,hole_mask,particle_mask,fock_diag_tmp,E0,p
do j1=1,N_particles(ispin1) do j1=1,N_particles(ispin1)
i_particle1 = particle_list(j1, ispin1) i_particle1 = particle_list(j1, ispin1)
p1 = i_particle1 + (ispin1 - 1) * mo_tot_num p1 = i_particle1 + (ispin1 - 1) * mo_tot_num
if(N_tmicrolist(p1) > 0 .and. idx_tmicrolist(ptr_tmicrolist(p1+1)-1) > i_generator) cycle if(N_tmicrolist(p1) > 0 .and. idx_tmicrolist(ptr_tmicrolist(p1)) < i_generator) cycle
jb = 1 jb = 1
if(ispin1 == ispin2) jb = j1+1 if(ispin1 == ispin2) jb = j1+1
do j2=jb,N_particles(ispin2) do j2=jb,N_particles(ispin2)
@ -587,7 +576,7 @@ subroutine select_doubles(i_generator,hole_mask,particle_mask,fock_diag_tmp,E0,p
p2 = i_particle2 + (ispin2 - 1) * mo_tot_num p2 = i_particle2 + (ispin2 - 1) * mo_tot_num
if(N_tmicrolist(p2) > 0 .and. idx_tmicrolist(ptr_tmicrolist(p2+1)-1) > i_generator) cycle if(N_tmicrolist(p2) > 0 .and. idx_tmicrolist(ptr_tmicrolist(p2)) < i_generator) cycle
if(isinwf(p1, p2)) cycle if(isinwf(p1, p2)) cycle
exc_det = ion_det exc_det = ion_det
@ -611,63 +600,28 @@ subroutine select_doubles(i_generator,hole_mask,particle_mask,fock_diag_tmp,E0,p
j_particle = i_particle1-ishft(k_particle-1,bit_kind_shift)-1 ! bit index j_particle = i_particle1-ishft(k_particle-1,bit_kind_shift)-1 ! bit index
exc_det(k_particle,ispin1) = ibset(exc_det(k_particle,ispin1),j_particle) exc_det(k_particle,ispin1) = ibset(exc_det(k_particle,ispin1),j_particle)
! if(.false. .or. (is_in_wavefunction(exc_det,N_int) .and. .not. isinwf(p1,p2))) then
! print *, p1, p2
! call debug_det(ion_det, N_int)
! call debug_det(exc_det, N_int)
! ! do i=1,mo_tot_num*2
! ! print *, isinwf(:, i)
! ! end do
! print *, isinwf(p1, p2)
! stop "isw"
! end if
! TODO
logical, external :: is_in_wavefunction logical, external :: is_in_wavefunction
logical :: nok logical :: nok
! TODO : Check connected to ref ! Compute perturbative contribution and select determinant
! if (.not. is_in_wavefunction(exc_det,N_int)) then double precision :: i_H_psi_value(N_states), i_H_psi_value2(N_states)
i_H_psi_value = 0d0
i_H_psi_value2 = 0d0
nok = .false.
call check_past_s(exc_det, microlist(1,1,ptr_microlist(sporb)), N_microlist(sporb) - N_futur_microlist(sporb), nok, N_int)
if(nok) cycle
!DET DRIVEN
if(N_futur_microlist(0) > 0) then
call i_H_psi(exc_det,microlist(1,1,ptr_futur_microlist(0)),psi_coef_microlist(ptr_futur_microlist(0), 1),N_int,N_futur_microlist(0),psi_selectors_size*4,N_states,i_H_psi_value)
end if
!INTEGRAL DRIVEN
! i_H_psi_value = d0s(mod(p1-1, mo_tot_num)+1, mod(p2-1, mo_tot_num)+1, :)
! Compute perturbative contribution and select determinant if(N_futur_microlist(sporb) > 0) then
double precision :: i_H_psi_value(N_states), i_H_psi_value2(N_states) !!! COMPUTE INTERSECTION
double precision :: i_H_full(N_states) !!!!!!!!!!!!!
i_H_psi_value = 0d0 ! if(dfloat(N_futur_microlist(lorb)) / dfloat(N_futur_microlist(sporb)) < 2d0) then
i_H_psi_value2 = 0d0
i_H_full = 0d0
! call i_H_psi(exc_det,psi_selectors,psi_selectors_coef,N_int,N_det_selectors,psi_selectors_size,N_states,i_H_full)
! call check_past(exc_det, microlist, idx_microlist, N_microlist(0), i_generator, nok, N_int)
! if(nok) cycle
nok = .false.
call check_futur(exc_det, microlist(1,1,ptr_futur_microlist(sporb)), N_futur_microlist(sporb), nok, N_int)
if(nok) cycle
if(N_microlist(0)-N_futur_microlist(0) > 0) then
call i_H_psi(exc_det,microlist(1,1,ptr_microlist(0)),psi_coef_microlist(ptr_microlist(0), 1),N_int,N_microlist(0)-N_futur_microlist(0),psi_selectors_size*4,N_states,i_H_psi_value)
! if(i_H_psi_value(1) /= d0s(p1, p2, 1) .and. d0s(p1, p2, 1) /= 0d0) then
! print *, d0s(p1, p2, 1), i_H_psi_value(1)
! print *, d0s(:3, :3, 1)
! stop "SKSL"
! end if
end if
! if(N_futur_microlist(sporb) > 0) call i_H_psi(exc_det,microlist(1,1,ptr_futur_microlist(sporb)),psi_coef_microlist(ptr_futur_microlist(sporb), 1),N_int,N_futur_microlist(sporb),psi_selectors_size*4,N_states,i_H_psi_value2)
! !$OMP ATOMIC
! d0 += dabs(i_H_psi_value(1))
! d2 += N_futur_microlist(sporb)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
if(N_microlist(sporb)-N_futur_microlist(sporb) > 0) then
! ! if(dfloat(N_futur_microlist(lorb)) / dfloat(N_futur_microlist(sporb)) < 2d0) then
! c1 = ptr_futur_microlist(p1) ! c1 = ptr_futur_microlist(p1)
! c2 = ptr_futur_microlist(p2) ! c2 = ptr_futur_microlist(p2)
! do while(c1 < ptr_microlist(p1+1) .and. c2 < ptr_microlist(p2+1)) ! do while(c1 < ptr_microlist(p1+1) .and. c2 < ptr_microlist(p2+1))
@ -685,112 +639,72 @@ subroutine select_doubles(i_generator,hole_mask,particle_mask,fock_diag_tmp,E0,p
! endif ! endif
! end do ! end do
! else ! else
call i_H_psi(exc_det,microlist(1,1,ptr_microlist(sporb)),psi_coef_microlist(ptr_microlist(sporb), 1),N_int,N_microlist(sporb)-N_futur_microlist(sporb),psi_selectors_size*4,N_states,i_H_psi_value2) call i_H_psi(exc_det,microlist(1,1,ptr_futur_microlist(sporb)),psi_coef_microlist(ptr_futur_microlist(sporb), 1),N_int,N_futur_microlist(sporb),psi_selectors_size*4,N_states,i_H_psi_value2)
i_H_psi_value = i_H_psi_value + i_H_psi_value2
end if
integer :: c1, c2
double precision :: hij
c1 = ptr_futur_tmicrolist(p1)
c2 = ptr_futur_tmicrolist(p2)
do while(.true.)
if(c1 >= ptr_tmicrolist(p1+1) .or. c2 >= ptr_tmicrolist(p2+1)) then
if(ptr_tmicrolist(p1+1) /= c1) then
call i_H_psi(exc_det,tmicrolist(1,1,c1),psi_coef_tmicrolist(c1, 1),N_int, ptr_tmicrolist(p1+1)-c1 ,psi_selectors_size*3,N_states,i_H_psi_value2)
i_H_psi_value = i_H_psi_value + i_H_psi_value2 i_H_psi_value = i_H_psi_value + i_H_psi_value2
! !$OMP ATOMIC
! d2 += dabs(i_H_psi_value2(1))
! d2 += N_futur_microlist(sporb)
! ! end if
end if
! enddo 2099.3283623955813
!!!!!!!!!!!!!!!!!!
! if(N_microlist(0) > 0) call i_H_psi(exc_det,microlist,psi_coef_microlist(ptr_microlist(0), 1),N_int,N_microlist(0),psi_selectors_size*4,N_states,i_H_psi_value)
! if(N_microlist(sporb) > 0) call i_H_psi(exc_det,microlist(1,1,ptr_microlist(sporb)),psi_coef_microlist(ptr_microlist(sporb), 1),N_int,N_microlist(sporb),psi_selectors_size*4,N_states,i_H_psi_value2)
! i_H_psi_value = i_H_psi_value + i_H_psi_value2
integer :: c1, c2
double precision :: hij
c1 = ptr_tmicrolist(p1)
c2 = ptr_tmicrolist(p2)
do while(.true.)
if(c1 >= ptr_futur_tmicrolist(p1) .or. c2 >= ptr_futur_tmicrolist(p2)) then
if(ptr_futur_tmicrolist(p1) /= c1) then
call i_H_psi(exc_det,tmicrolist(1,1,c1),psi_coef_tmicrolist(c1, 1),N_int, ptr_futur_tmicrolist(p1)-c1 ,psi_selectors_size*3,N_states,i_H_psi_value2)
i_H_psi_value = i_H_psi_value + i_H_psi_value2
! ! !$OMP ATOMIC
! d1 += dabs(i_H_psi_value2(1))
end if
if(ptr_futur_tmicrolist(p2) /= c2) then
call i_H_psi(exc_det,tmicrolist(1,1,c2),psi_coef_tmicrolist(c2, 1),N_int, ptr_futur_tmicrolist(p2)-c2 ,psi_selectors_size*3,N_states,i_H_psi_value2)
i_H_psi_value = i_H_psi_value + i_H_psi_value2
! !$OMP ATOMIC
! d1 += dabs(i_H_psi_value2(1))
endif
exit
endif
if(idx_tmicrolist(c1) < idx_tmicrolist(c2)) then
call i_H_j(exc_det,tmicrolist(1,1,c1),N_int,hij)
do j = 1, N_states
i_H_psi_value(j) = i_H_psi_value(j) + psi_coef_tmicrolist(c1,j)*hij
enddo
! !$OMP ATOMIC
! d1 += dabs(psi_coef_tmicrolist(c1,1)*hij)
c1 += 1
else
call i_H_j(exc_det,tmicrolist(1,1,c2),N_int,hij)
do j = 1, N_states
i_H_psi_value(j) = i_H_psi_value(j) + psi_coef_tmicrolist(c2,j)*hij
enddo
if(idx_tmicrolist(c1) == idx_tmicrolist(c2)) c1 = c1 + 1
! !$OMP ATOMIC
! d1 += dabs(psi_coef_tmicrolist(c2,1)*hij)
c2 += 1
end if end if
enddo
double precision :: Hii, diag_H_mat_elem_fock if(ptr_tmicrolist(p2+1) /= c2) then
Hii = diag_H_mat_elem_fock(psi_det_generators(1,1,i_generator),exc_det,fock_diag_tmp,N_int) call i_H_psi(exc_det,tmicrolist(1,1,c2),psi_coef_tmicrolist(c2, 1),N_int, ptr_tmicrolist(p2+1)-c2 ,psi_selectors_size*3,N_states,i_H_psi_value2)
double precision :: delta_E, e_pert(N_states), e_pertm i_H_psi_value = i_H_psi_value + i_H_psi_value2
e_pert(:) = 0d0
e_pertm = 0d0
do k=1,N_states
! if(dabs(1d0 - i_H_psi_value(k)/i_H_full(k)) > 1d-6) then
! print *, i_H_psi_value(k), i_H_full(k), i_H_psi_value(k)/i_H_full(k)
! stop "PAS BON, PAS BOOON (double)"
!
! endif
if (i_H_psi_value(k) == 0.d0) cycle
delta_E = E0(k) - Hii
if (delta_E < 0.d0) then
e_pert(k) = 0.5d0 * (-dsqrt(delta_E * delta_E + 4.d0 * i_H_psi_value(k) * i_H_psi_value(k)) - delta_E)
else
e_pert(k) = 0.5d0 * ( dsqrt(delta_E * delta_E + 4.d0 * i_H_psi_value(k) * i_H_psi_value(k)) - delta_E)
endif endif
if(dabs(e_pert(k)) > dabs(e_pertm)) e_pertm = e_pert(k)
pt2(k) += e_pert(k) exit
enddo endif
if(dabs(e_pertm) > dabs(buf%mini)) then
call add_to_selection_buffer(buf, exc_det, e_pertm) if(idx_tmicrolist(c1) < idx_tmicrolist(c2)) then
call i_H_j(exc_det,tmicrolist(1,1,c1),N_int,hij)
do j = 1, N_states
i_H_psi_value(j) = i_H_psi_value(j) + psi_coef_tmicrolist(c1,j)*hij
enddo
c1 += 1
else
call i_H_j(exc_det,tmicrolist(1,1,c2),N_int,hij)
do j = 1, N_states
i_H_psi_value(j) = i_H_psi_value(j) + psi_coef_tmicrolist(c2,j)*hij
enddo
if(idx_tmicrolist(c1) == idx_tmicrolist(c2)) c1 = c1 + 1
c2 += 1
end if end if
! endif ! iwf enddo
double precision :: Hii, diag_H_mat_elem_fock
Hii = diag_H_mat_elem_fock(psi_det_generators(1,1,i_generator),exc_det,fock_diag_tmp,N_int)
double precision :: delta_E, e_pert(N_states), e_pertm
e_pert(:) = 0d0
e_pertm = 0d0
do k=1,N_states
if (i_H_psi_value(k) == 0.d0) cycle
delta_E = E0(k) - Hii
if (delta_E < 0.d0) then
e_pert(k) = 0.5d0 * (-dsqrt(delta_E * delta_E + 4.d0 * i_H_psi_value(k) * i_H_psi_value(k)) - delta_E)
! Reset exc_det else
! exc_det(k_particle,ispin) = psi_det_generators(k_particle,ispin,i_generator) e_pert(k) = 0.5d0 * ( dsqrt(delta_E * delta_E + 4.d0 * i_H_psi_value(k) * i_H_psi_value(k)) - delta_E)
enddo ! j endif
if(dabs(e_pert(k)) > dabs(e_pertm)) e_pertm = e_pert(k)
pt2(k) += e_pert(k)
enddo
if(dabs(e_pertm) > dabs(buf%mini)) then
call add_to_selection_buffer(buf, exc_det, e_pertm)
end if
enddo
enddo enddo
! Reset ion_det
! ion_det(k_hole,ispin) = psi_det_generators(k_hole,ispin,i_generator)
enddo ! i
enddo enddo
enddo ! ispin enddo
enddo
enddo enddo
!print *, "D ::: ", d0/1000000, d1/1000000, d2/1000000
end end
@ -804,7 +718,7 @@ subroutine create_futur_ptr(ptr_microlist, idx_microlist, ptr_futur_microlist, N
do i=0,mo_tot_num*2 do i=0,mo_tot_num*2
ptr_futur_microlist(i) = ptr_microlist(i+1) ptr_futur_microlist(i) = ptr_microlist(i+1)
do j=ptr_microlist(i), ptr_microlist(i+1) - 1 do j=ptr_microlist(i), ptr_microlist(i+1) - 1
if(idx_microlist(j) > i_generator) then if(idx_microlist(j) >= i_generator) then
ptr_futur_microlist(i) = j ptr_futur_microlist(i) = j
N_futur_microlist(i) = ptr_microlist(i+1) - j N_futur_microlist(i) = ptr_microlist(i+1) - j
exit exit
@ -849,7 +763,7 @@ subroutine create_microlist_single(minilist, i_cur, N_minilist, key_mask, microl
if(nt > 3) then !! TOO MANY DIFFERENCES if(nt > 3) then !! TOO MANY DIFFERENCES
continue continue
else if(nt < 3) then else if(nt < 3) then
if(i > i_cur) then if(i < i_cur) then !!!!!!!!!!!!!!!!!!!!! DESACTIVADO
N_microlist(:) = 0 !!!! PAST LINKED TO EVERYBODY! N_microlist(:) = 0 !!!! PAST LINKED TO EVERYBODY!
ptr_microlist(:) = 1 ptr_microlist(:) = 1
return return
@ -862,8 +776,8 @@ subroutine create_microlist_single(minilist, i_cur, N_minilist, key_mask, microl
do s=1,2 do s=1,2
do j=1,n_element(s) do j=1,n_element(s)
nt = list(j,s) + mo_tot_num * (s-1) nt = list(j,s) + mo_tot_num * (s-1)
N_microlist(nt) = N_microlist(nt) + 1 N_microlist(nt) = N_microlist(nt) + 1
end do end do
end do end do
endif endif
@ -970,13 +884,17 @@ subroutine create_microlist_double(minilist, i_cur, N_minilist, key_mask, microl
integer(bit_kind), intent(out) :: tmicrolist(Nint,2,N_minilist*4) integer(bit_kind), intent(out) :: tmicrolist(Nint,2,N_minilist*4)
integer :: i,j,k,s,nt,nt2,n_element(2,N_minilist), idx(0:N_minilist) integer :: i,j,k,s,nt,nt2
integer :: list(4,2,N_minilist), cur_microlist(0:mo_tot_num*2+1), cur_tmicrolist(0:mo_tot_num*2+1) integer, allocatable :: n_element(:,:), idx(:), list(:,:,:)
integer :: cur_microlist(0:mo_tot_num*2+1), cur_tmicrolist(0:mo_tot_num*2+1)
integer(bit_kind) :: key_mask_neg(Nint,2), mobileMask(Nint,2) integer(bit_kind) :: key_mask_neg(Nint,2), mobileMask(Nint,2)
integer :: mo_tot_num_2 integer :: mo_tot_num_2
logical,intent(out) :: isinwf(mo_tot_num*2, mo_tot_num*2) logical,intent(out) :: isinwf(mo_tot_num*2, mo_tot_num*2)
double precision, intent(out) :: d0s(mo_tot_num, mo_tot_num, N_states) double precision, intent(out) :: d0s(mo_tot_num, mo_tot_num, N_states)
double precision :: integ(mo_tot_num, mo_tot_num) double precision :: integ(mo_tot_num, mo_tot_num)
allocate(list(4,2,N_minilist), n_element(2,N_minilist), idx(0:N_minilist))
isinwf = .false. isinwf = .false.
integ = 0d0 integ = 0d0
d0s = 0d0 d0s = 0d0
@ -1002,7 +920,6 @@ subroutine create_microlist_double(minilist, i_cur, N_minilist, key_mask, microl
end do end do
if(nt > 4) cycle !! TOO MANY DIFFERENCES if(nt > 4) cycle !! TOO MANY DIFFERENCES
idx(0) += 1 idx(0) += 1
idx(idx(0)) = i idx(idx(0)) = i
@ -1010,34 +927,15 @@ subroutine create_microlist_double(minilist, i_cur, N_minilist, key_mask, microl
call bitstring_to_list(mobileMask(1,2), list(1,2,idx(0)), n_element(2, idx(0)), Nint) call bitstring_to_list(mobileMask(1,2), list(1,2,idx(0)), n_element(2, idx(0)), Nint)
if(nt <= 2) then if(nt == 2) then
if(i > i_cur) then if(i < i_cur) then
N_microlist = 0 N_microlist(:) = 0
ptr_microlist = 1 ptr_microlist = 1
N_tmicrolist = 0 N_tmicrolist = 0
ptr_tmicrolist = 1 ptr_tmicrolist = 1
return return
else else
!n_element(:, idx(0)) = (/2, 0/)
N_microlist(0) = N_microlist(0) + 1 N_microlist(0) = N_microlist(0) + 1
if(n_element(1,idx(0)) >= 1) nt = list(1,1,idx(0))
if(n_element(1,idx(0)) == 2) nt2 = list(2,1,idx(0))
if(n_element(2,idx(0)) == 2) nt = list(2,2,idx(0)) + mo_tot_num
if(n_element(2,idx(0)) >= 1) nt2 = list(1,2,idx(0)) + mo_tot_num
isinwf(nt, nt2) = .true.
isinwf(nt2, nt) = .true.
double precision, external :: get_mo_bielec_integral
nt = mod(nt, mo_tot_num)
nt2 = mod(nt2, mo_tot_num)
! call get_mo_bielec_integrals_ij(nt, nt2 ,mo_tot_num,integ,mo_integrals_map)
! do j=1, N_states
! call i_h_j
! d0s(:,:,j) += integ(:,:) * psi_selectors_coef(i,j) !!!!!!!!!!!!!!!!!!! MOOOOOCHE !!!!! suppose que minilist = psi_selectors .....
! end do
! print *, "TO ", integ(mod(nt, mo_tot_num), mod(nt2, mo_tot_num))
endif endif
else else
do s=1,2 do s=1,2
@ -1064,7 +962,6 @@ subroutine create_microlist_double(minilist, i_cur, N_minilist, key_mask, microl
do i=1, idx(0) do i=1, idx(0)
if(n_element(1, i) + n_element(2, i) > 4) stop "wired"
if(n_element(1, i) + n_element(2, i) <= 2) then if(n_element(1, i) + n_element(2, i) <= 2) then
idx_microlist(cur_microlist(0)) = idx(i) idx_microlist(cur_microlist(0)) = idx(i)
do k=1,Nint do k=1,Nint
@ -1072,6 +969,26 @@ subroutine create_microlist_double(minilist, i_cur, N_minilist, key_mask, microl
microlist(k,2,cur_microlist(0)) = minilist(k,2,idx(i)) microlist(k,2,cur_microlist(0)) = minilist(k,2,idx(i))
enddo enddo
cur_microlist(0) = cur_microlist(0) + 1 cur_microlist(0) = cur_microlist(0) + 1
if(n_element(1,i) >= 1) nt = list(1,1,i)
if(n_element(1,i) == 2) nt2 = list(2,1,i)
if(n_element(2,i) == 2) nt = list(2,2,i) + mo_tot_num
if(n_element(2,i) >= 1) nt2 = list(1,2,i) + mo_tot_num
isinwf(nt, nt2) = .true.
isinwf(nt2, nt) = .true.
!!!! INTEGRAL DRIVEN
!!!!!!!!!!!!!!!!!!!!
! call get_d0(minilist(1,1,idx(i)), integ, key_mask, 1+(nt2-1)/mo_tot_num, 1+(nt-1)/mo_tot_num, &
! mod(nt2-1, mo_tot_num)+1, mod(nt-1, mo_tot_num)+1)
!
! do j=1, N_states
! do nt=1, mo_tot_num
! do nt2=1, mo_tot_num
! d0s(nt,nt2,j) = d0s(nt,nt2,j) + (integ(nt,nt2) * psi_selectors_coef(idx(i), j)) !!! SUPPOSE MINILIST = SELECTORS !!!!
! end do
! end do
! end do
else else
do s = 1, 2 do s = 1, 2
do j=1,n_element(s,i) do j=1,n_element(s,i)
@ -1109,8 +1026,8 @@ subroutine check_past(det, list, idx, N, cur, ok, Nint)
integer :: i,s,ni integer :: i,s,ni
ok = .false. ok = .false.
do i=N,1,-1 do i=1,N
if(idx(i) <= cur) exit if(idx(i) >= cur) exit
s = 0 s = 0
do ni=1,Nint do ni=1,Nint
s += popcnt(xor(det(ni,1), list(ni,1,i))) + popcnt(xor(det(ni,2), list(ni,2,i))) s += popcnt(xor(det(ni,1), list(ni,1,i))) + popcnt(xor(det(ni,2), list(ni,2,i)))
@ -1123,7 +1040,7 @@ subroutine check_past(det, list, idx, N, cur, ok, Nint)
end subroutine end subroutine
subroutine check_futur(det, list, N, ok, Nint) subroutine check_past_s(det, list, N, ok, Nint)
implicit none implicit none
use bitmasks use bitmasks
@ -1145,3 +1062,204 @@ subroutine check_futur(det, list, N, ok, Nint)
end do end do
end subroutine end subroutine
subroutine get_d0(gen, mat, mask, s1, s2, h1, h2)
use bitmasks
implicit none
double precision, intent(out) :: mat(mo_tot_num, mo_tot_num)
double precision :: mat_mwen(mo_tot_num, mo_tot_num)
integer, intent(in) :: h1, h2, s1, s2
integer(bit_kind), intent(in) :: mask(N_int, 2), gen(N_int, 2)
integer(bit_kind) :: det1(N_int, 2), det2(N_int, 2)
logical :: ok, mono
double precision :: phase, phase2, inv
integer :: p1, p2, hmi, hma
logical, external :: detEq
integer :: exc(0:2, 2, 2), exc2(0:2,2,2)
exc = 0
call get_mo_bielec_integrals_ij(h1, h2 ,mo_tot_num,mat_mwen,mo_integrals_map)
mat = 0d0
if(s1 == s2) then
hmi = min(h1, h2)
hma = max(h1, h2)
inv = 1d0
if(h1 > h2) inv = -1d0
exc(0, :, s1) = 2
exc(1, 1, s1) = hmi
exc(2, 1, s1) = hma
do p1=1,mo_tot_num
do p2=1,mo_tot_num
if(p1 == p2) cycle
call apply_particle(mask, (/s1,p1,s2,p2/), det2, ok, N_int)
if(.not. ok) cycle
mono = (hmi == p1 .or. hma == p2 .or. hmi == p2 .or. hma == p1)
if(mono) then
call i_h_j(gen, det2, N_int, mat(p1, p2))
!mat(p1, p2) = phase
else
exc(1, 2, s1) = min(p1, p2)
exc(2, 2, s2) = max(p2, p1)
call get_double_excitation_phase(gen, det2, exc, phase, N_int)
mat(p1, p2) = inv * (mat_mwen(p1, p2) - mat_mwen(p2, p1)) * phase
end if
end do
end do
else
exc(0, :, 1) = 1
exc(0, :, 2) = 1
if(s1 /= 2) stop "alpha beta inversified"
exc(1, 1, 1) = h2
exc(1, 1, 2) = h1
do p1=1, mo_tot_num
do p2=1, mo_tot_num
call apply_particle(mask, (/s1,p1,s2,p2/), det2, ok, N_int)
if(.not. ok) cycle
mono = (h1 == p1 .or. h2 == p2)
if(mono) then
call i_h_j(gen, det2, N_int, phase)
mat(p1, p2) = phase
else
exc(1, 2, s1) = p1
exc(1, 2, s2) = p2
call get_double_excitation_phase(gen, det2, exc, phase, N_int)
mat(p1, p2) = mat_mwen(p1, p2) * phase
end if
end do
end do
end if
end subroutine
subroutine apply_particle(det, exc, res, ok, Nint)
use bitmasks
implicit none
integer, intent(in) :: Nint
integer, intent(in) :: exc(4)
integer :: s1, s2, p1, p2
integer(bit_kind),intent(in) :: det(Nint, 2)
integer(bit_kind),intent(out) :: res(Nint, 2)
logical, intent(out) :: ok
integer :: ii, pos
ok = .false.
s1 = exc(1)
p1 = exc(2)
s2 = exc(3)
p2 = exc(4)
res = det
if(p1 /= 0) then
ii = (p1-1)/bit_kind_size + 1
pos = mod(p1-1, 64)!iand(p1-1,bit_kind_size-1)
if(iand(det(ii, s1), ishft(1_bit_kind, pos)) /= 0_8) return
res(ii, s1) = ibset(res(ii, s1), pos)
end if
ii = (p2-1)/bit_kind_size + 1
pos = mod(p2-1, 64)!iand(p2-1,bit_kind_size-1)
if(iand(det(ii, s2), ishft(1_bit_kind, pos)) /= 0_8) return
res(ii, s2) = ibset(res(ii, s2), pos)
ok = .true.
end subroutine
subroutine get_double_excitation_phase(det1,det2,exc,phase,Nint)
use bitmasks
implicit none
BEGIN_DOC
! Returns the two excitation operators between two doubly excited determinants and the phase
END_DOC
integer, intent(in) :: Nint
integer(bit_kind), intent(in) :: det1(Nint,2)
integer(bit_kind), intent(in) :: det2(Nint,2)
integer, intent(in) :: exc(0:2,2,2)
double precision, intent(out) :: phase
integer :: tz
integer :: l, ispin, idx_hole, idx_particle, ishift
integer :: nperm
integer :: i,j,k,m,n
integer :: high, low
integer :: a,b,c,d
integer(bit_kind) :: hole, particle, tmp
double precision, parameter :: phase_dble(0:1) = (/ 1.d0, -1.d0 /)
ASSERT (Nint > 0)
nperm = 0
do ispin = 1,2
select case (exc(0,1,ispin))
case(0)
cycle
case(1)
low = min(exc(1,1,ispin), exc(1,2,ispin))
high = max(exc(1,1,ispin), exc(1,2,ispin))
ASSERT (low > 0)
j = ishft(low-1,-bit_kind_shift)+1 ! Find integer in array(Nint)
n = iand(low-1,bit_kind_size-1)+1 ! mod(low,bit_kind_size)
ASSERT (high > 0)
k = ishft(high-1,-bit_kind_shift)+1
m = iand(high-1,bit_kind_size-1)+1
if (j==k) then
nperm = nperm + popcnt(iand(det1(j,ispin), &
iand( ibset(0_bit_kind,m-1)-1_bit_kind, &
ibclr(-1_bit_kind,n)+1_bit_kind ) ))
else
nperm = nperm + popcnt(iand(det1(k,ispin), &
ibset(0_bit_kind,m-1)-1_bit_kind))
if (n < bit_kind_size) then
nperm = nperm + popcnt(iand(det1(j,ispin), ibclr(-1_bit_kind,n) +1_bit_kind))
endif
do i=j+1,k-1
nperm = nperm + popcnt(det1(i,ispin))
end do
endif
case (2)
do i=1,2
low = min(exc(i,1,ispin), exc(i,2,ispin))
high = max(exc(i,1,ispin), exc(i,2,ispin))
ASSERT (low > 0)
j = ishft(low-1,-bit_kind_shift)+1 ! Find integer in array(Nint)
n = iand(low-1,bit_kind_size-1)+1 ! mod(low,bit_kind_size)
ASSERT (high > 0)
k = ishft(high-1,-bit_kind_shift)+1
m = iand(high-1,bit_kind_size-1)+1
if (j==k) then
nperm = nperm + popcnt(iand(det1(j,ispin), &
iand( ibset(0_bit_kind,m-1)-1_bit_kind, &
ibclr(-1_bit_kind,n)+1_bit_kind ) ))
else
nperm = nperm + popcnt(iand(det1(k,ispin), &
ibset(0_bit_kind,m-1)-1_bit_kind))
if (n < bit_kind_size) then
nperm = nperm + popcnt(iand(det1(j,ispin), ibclr(-1_bit_kind,n) +1_bit_kind))
endif
do l=j+1,k-1
nperm = nperm + popcnt(det1(l,ispin))
end do
endif
enddo
a = min(exc(1,1,ispin), exc(1,2,ispin))
b = max(exc(1,1,ispin), exc(1,2,ispin))
c = min(exc(2,1,ispin), exc(2,2,ispin))
d = max(exc(2,1,ispin), exc(2,2,ispin))
if (c>a .and. c<b .and. d>b) then
nperm = nperm + 1
endif
exit
end select
enddo
phase = phase_dble(iand(nperm,1))
end