diff --git a/plugins/local/slater_tc/h_mat_triple.irp.f b/plugins/local/slater_tc/h_mat_triple.irp.f deleted file mode 100644 index 6f5697a2..00000000 --- a/plugins/local/slater_tc/h_mat_triple.irp.f +++ /dev/null @@ -1,391 +0,0 @@ -subroutine get_excitation_general(key_i,key_j, Nint,degree_array,holes_array, particles_array,phase) - use bitmasks - BEGIN_DOC -! returns the array, for each spin, of holes/particles between key_i and key_j -! -! with the following convention: a^+_{particle} a_{hole}|key_i> = |key_j> - END_DOC - include 'utils/constants.include.F' - implicit none - integer, intent(in) :: Nint - integer(bit_kind), intent(in) :: key_j(Nint,2),key_i(Nint,2) - integer, intent(out) :: holes_array(100,2),particles_array(100,2),degree_array(2) - double precision, intent(out) :: phase - integer :: ispin,k,i,pos - integer(bit_kind) :: key_hole, key_particle - integer(bit_kind) :: xorvec(N_int_max,2) - holes_array = -1 - particles_array = -1 - degree_array = 0 - do i = 1, N_int - xorvec(i,1) = xor( key_i(i,1), key_j(i,1)) - xorvec(i,2) = xor( key_i(i,2), key_j(i,2)) - degree_array(1) += popcnt(xorvec(i,1)) - degree_array(2) += popcnt(xorvec(i,2)) - enddo - degree_array(1) = shiftr(degree_array(1),1) - degree_array(2) = shiftr(degree_array(2),1) - - do ispin = 1, 2 - k = 1 - !!! GETTING THE HOLES - do i = 1, N_int - key_hole = iand(xorvec(i,ispin),key_i(i,ispin)) - do while(key_hole .ne.0_bit_kind) - pos = trailz(key_hole) - holes_array(k,ispin) = 1+ bit_kind_size * (i-1) + pos - key_hole = ibclr(key_hole,pos) - k += 1 - if(k .gt.100)then - print*,'WARNING in get_excitation_general' - print*,'More than a 100-th excitation for spin ',ispin - print*,'stoping ...' - stop - endif - enddo - enddo - enddo - do ispin = 1, 2 - k = 1 - !!! GETTING THE PARTICLES - do i = 1, N_int - key_particle = iand(xor(key_i(i,ispin),key_j(i,ispin)),key_j(i,ispin)) - do while(key_particle .ne.0_bit_kind) - pos = trailz(key_particle) - particles_array(k,ispin) = 1+ bit_kind_size * (i-1) + pos - key_particle = ibclr(key_particle,pos) - k += 1 - if(k .gt.100)then - print*,'WARNING in get_excitation_general ' - print*,'More than a 100-th excitation for spin ',ispin - print*,'stoping ...' - stop - endif - enddo - enddo - enddo - integer :: h,p, i_ok - integer(bit_kind), allocatable :: det_i(:,:),det_ip(:,:) - integer :: exc(0:2,2,2) - double precision :: phase_tmp - allocate(det_i(Nint,2),det_ip(N_int,2)) - det_i = key_i - phase = 1.d0 - do ispin = 1, 2 - do i = 1, degree_array(ispin) - h = holes_array(i,ispin) - p = particles_array(i,ispin) - det_ip = det_i - call do_single_excitation(det_ip,h,p,ispin,i_ok) - if(i_ok == -1)then - print*,'excitation was not possible ' - stop - endif - call get_single_excitation(det_i,det_ip,exc,phase_tmp,Nint) - phase *= phase_tmp - det_i = det_ip - enddo - enddo - -end - -subroutine get_holes_general(key_i, key_j,Nint, holes_array) - use bitmasks - BEGIN_DOC -! returns the array, per spin, of holes between key_i and key_j -! -! with the following convention: a_{hole}|key_i> --> |key_j> - END_DOC - implicit none - integer, intent(in) :: Nint - integer(bit_kind), intent(in) :: key_j(Nint,2),key_i(Nint,2) - integer, intent(out) :: holes_array(100,2) - integer(bit_kind) :: key_hole - integer :: ispin,k,i,pos - holes_array = -1 - do ispin = 1, 2 - k = 1 - do i = 1, N_int - key_hole = iand(xor(key_i(i,ispin),key_j(i,ispin)),key_i(i,ispin)) - do while(key_hole .ne.0_bit_kind) - pos = trailz(key_hole) - holes_array(k,ispin) = 1+ bit_kind_size * (i-1) + pos - key_hole = ibclr(key_hole,pos) - k += 1 - if(k .gt.100)then - print*,'WARNING in get_holes_general' - print*,'More than a 100-th excitation for spin ',ispin - print*,'stoping ...' - stop - endif - enddo - enddo - enddo -end - -subroutine get_particles_general(key_i, key_j,Nint,particles_array) - use bitmasks - BEGIN_DOC -! returns the array, per spin, of particles between key_i and key_j -! -! with the following convention: a^dagger_{particle}|key_i> --> |key_j> - END_DOC - implicit none - integer, intent(in) :: Nint - integer(bit_kind), intent(in) :: key_j(Nint,2),key_i(Nint,2) - integer, intent(out) :: particles_array(100,2) - integer(bit_kind) :: key_particle - integer :: ispin,k,i,pos - particles_array = -1 - do ispin = 1, 2 - k = 1 - do i = 1, N_int - key_particle = iand(xor(key_i(i,ispin),key_j(i,ispin)),key_j(i,ispin)) - do while(key_particle .ne.0_bit_kind) - pos = trailz(key_particle) - particles_array(k,ispin) = 1+ bit_kind_size * (i-1) + pos - key_particle = ibclr(key_particle,pos) - k += 1 - if(k .gt.100)then - print*,'WARNING in get_holes_general' - print*,'More than a 100-th excitation for spin ',ispin - print*,'Those are the two determinants' - call debug_det(key_i, N_int) - call debug_det(key_j, N_int) - print*,'stoping ...' - stop - endif - enddo - enddo - enddo -end - -subroutine get_phase_general(key_i,Nint,degree, holes_array, particles_array,phase) - implicit none - integer, intent(in) :: degree(2), Nint - integer(bit_kind), intent(in) :: key_i(Nint,2) - integer, intent(in) :: holes_array(100,2),particles_array(100,2) - double precision, intent(out) :: phase - integer :: i,ispin,h,p, i_ok - integer(bit_kind), allocatable :: det_i(:,:),det_ip(:,:) - integer :: exc(0:2,2,2) - double precision :: phase_tmp - allocate(det_i(Nint,2),det_ip(N_int,2)) - det_i = key_i - phase = 1.d0 - do ispin = 1, 2 - do i = 1, degree(ispin) - h = holes_array(i,ispin) - p = particles_array(i,ispin) - det_ip = det_i - call do_single_excitation(det_ip,h,p,ispin,i_ok) - if(i_ok == -1)then - print*,'excitation was not possible ' - stop - endif - call get_single_excitation(det_i,det_ip,exc,phase_tmp,Nint) - phase *= phase_tmp - det_i = det_ip - enddo - enddo - -end - -subroutine H_tc_s2_u_0_with_pure_three(v_0, s_0, u_0, N_st, sze) - BEGIN_DOC - ! Computes $v_0 = H^TC | u_0\rangle$ WITH PURE TRIPLE EXCITATION TERMS - ! - ! Assumes that the determinants are in psi_det - ! - ! istart, iend, ishift, istep are used in ZMQ parallelization. - END_DOC - - use bitmasks - implicit none - - integer, intent(in) :: N_st,sze - double precision, intent(in) :: u_0(sze,N_st) - double precision, intent(out) :: v_0(sze,N_st), s_0(sze,N_st) - call H_tc_s2_u_0_opt(v_0, s_0, u_0, N_st, sze) - integer :: i,j,degree,ist - double precision :: hmono, htwoe, hthree, htot - do i = 1, N_det - do j = 1, N_det - call get_excitation_degree(psi_det(1,1,i),psi_det(1,1,j),degree,N_int) - if(degree .ne. 3)cycle - call triple_htilde_mu_mat_fock_bi_ortho(N_int, psi_det(1,1,i), psi_det(1,1,j), hmono, htwoe, hthree, htot) - do ist = 1, N_st - v_0(i,ist) += htot * u_0(j,ist) - enddo - enddo - enddo -end - -subroutine H_tc_s2_u_0_with_pure_three_omp(v_0, s_0, u_0, N_st, sze) - BEGIN_DOC - ! Computes $v_0 = H^TC | u_0\rangle$ WITH PURE TRIPLE EXCITATION TERMS - ! - ! Assumes that the determinants are in psi_det - ! - ! istart, iend, ishift, istep are used in ZMQ parallelization. - END_DOC - - use bitmasks - implicit none - - integer, intent(in) :: N_st,sze - double precision, intent(in) :: u_0(sze,N_st) - double precision, intent(out) :: v_0(sze,N_st), s_0(sze,N_st) - call H_tc_s2_u_0_opt(v_0, s_0, u_0, N_st, sze) - integer :: i,j,degree,ist - double precision :: hmono, htwoe, hthree, htot - !$OMP PARALLEL DO DEFAULT(NONE) SCHEDULE(dynamic,8) & - !$OMP SHARED(N_st, N_det, N_int, psi_det, u_0, v_0) & - !$OMP PRIVATE(ist, i, j, degree, hmono, htwoe, hthree,htot) - do i = 1, N_det - do j = 1, N_det - call get_excitation_degree(psi_det(1,1,i),psi_det(1,1,j),degree,N_int) - if(degree .ne. 3)cycle - call triple_htilde_mu_mat_fock_bi_ortho(N_int, psi_det(1,1,i), psi_det(1,1,j), hmono, htwoe, hthree, htot) - do ist = 1, N_st - v_0(i,ist) += htot * u_0(j,ist) - enddo - enddo - enddo - !$OMP END PARALLEL DO -end - -! --- - -subroutine H_tc_s2_dagger_u_0_with_pure_three(v_0, s_0, u_0, N_st, sze) - BEGIN_DOC - ! Computes $v_0 = (H^TC)^dagger | u_0\rangle$ WITH PURE TRIPLE EXCITATION TERMS - ! - ! Assumes that the determinants are in psi_det - ! - ! istart, iend, ishift, istep are used in ZMQ parallelization. - END_DOC - - use bitmasks - implicit none - - integer, intent(in) :: N_st,sze - double precision, intent(in) :: u_0(sze,N_st) - double precision, intent(out) :: v_0(sze,N_st), s_0(sze,N_st) - call H_tc_s2_dagger_u_0_opt(v_0, s_0, u_0, N_st, sze) - integer :: i,j,degree,ist - double precision :: hmono, htwoe, hthree, htot - do i = 1, N_det - do j = 1, N_det - call get_excitation_degree(psi_det(1,1,i),psi_det(1,1,j),degree,N_int) - if(degree .ne. 3)cycle - call triple_htilde_mu_mat_fock_bi_ortho(N_int, psi_det(1,1,j), psi_det(1,1,i), hmono, htwoe, hthree, htot) - do ist = 1, N_st - v_0(i,ist) += htot * u_0(j,ist) - enddo - enddo - enddo -end - -subroutine H_tc_s2_dagger_u_0_with_pure_three_omp(v_0, s_0, u_0, N_st, sze) - BEGIN_DOC - ! Computes $v_0 = (H^TC)^dagger | u_0\rangle$ WITH PURE TRIPLE EXCITATION TERMS - ! - ! Assumes that the determinants are in psi_det - ! - ! istart, iend, ishift, istep are used in ZMQ parallelization. - END_DOC - - use bitmasks - implicit none - - integer, intent(in) :: N_st,sze - double precision, intent(in) :: u_0(sze,N_st) - double precision, intent(out) :: v_0(sze,N_st), s_0(sze,N_st) - call H_tc_s2_dagger_u_0_opt(v_0, s_0, u_0, N_st, sze) - integer :: i,j,degree,ist - double precision :: hmono, htwoe, hthree, htot - !$OMP PARALLEL DO DEFAULT(NONE) SCHEDULE(dynamic,8) & - !$OMP SHARED(N_st, N_det, N_int, psi_det, u_0, v_0) & - !$OMP PRIVATE(ist, i, j, degree, hmono, htwoe, hthree,htot) - do i = 1, N_det - do j = 1, N_det - call get_excitation_degree(psi_det(1,1,i),psi_det(1,1,j),degree,N_int) - if(degree .ne. 3)cycle - call triple_htilde_mu_mat_fock_bi_ortho(N_int, psi_det(1,1,j), psi_det(1,1,i), hmono, htwoe, hthree, htot) - do ist = 1, N_st - v_0(i,ist) += htot * u_0(j,ist) - enddo - enddo - enddo - !$OMP END PARALLEL DO -end - -! --- -subroutine triple_htilde_mu_mat_fock_bi_ortho(Nint, key_j, key_i, hmono, htwoe, hthree, htot) - use bitmasks - BEGIN_DOC -! for triple excitation -!! -!! WARNING !! -! -! Genuine triple excitations of the same spin are not yet implemented - END_DOC - implicit none - integer(bit_kind), intent(in) :: key_j(N_int,2),key_i(N_int,2) - integer, intent(in) :: Nint - double precision, intent(out) :: hmono, htwoe, hthree, htot - integer :: degree - integer :: h1, p1, h2, p2, s1, s2, h3, p3, s3 - integer :: holes_array(100,2),particles_array(100,2),degree_array(2) - double precision :: phase,sym_3_e_int_from_6_idx_tensor - - hmono = 0.d0 - htwoe = 0.d0 - hthree = 0.d0 - htot = 0.d0 - call get_excitation_general(key_j, key_i, Nint,degree_array,holes_array, particles_array,phase) - degree = degree_array(1) + degree_array(2) - if(degree .ne. 3)return - if(degree_array(1)==3.or.degree_array(2)==3)then - if(degree_array(1) == 3)then - h1 = holes_array(1,1) - h2 = holes_array(2,1) - h3 = holes_array(3,1) - p1 = particles_array(1,1) - p2 = particles_array(2,1) - p3 = particles_array(3,1) - else - h1 = holes_array(1,2) - h2 = holes_array(2,2) - h3 = holes_array(3,2) - p1 = particles_array(1,2) - p2 = particles_array(2,2) - p3 = particles_array(3,2) - endif - hthree = sym_3_e_int_from_6_idx_tensor(p3, p2, p1, h3, h2, h1) - else - if(degree_array(1) == 2.and.degree_array(2) == 1)then ! double alpha + single beta - h1 = holes_array(1,1) - h2 = holes_array(2,1) - h3 = holes_array(1,2) - p1 = particles_array(1,1) - p2 = particles_array(2,1) - p3 = particles_array(1,2) - else if(degree_array(2) == 2 .and. degree_array(1) == 1)then ! double beta + single alpha - h1 = holes_array(1,2) - h2 = holes_array(2,2) - h3 = holes_array(1,1) - p1 = particles_array(1,2) - p2 = particles_array(2,2) - p3 = particles_array(1,1) - else - print*,'PB !!' - stop - endif - hthree = three_body_ints_bi_ort(p3,p2,p1,h3,h2,h1) - three_body_ints_bi_ort(p3,p2,p1,h3,h1,h2) - endif - hthree *= phase - htot = hthree - end - diff --git a/plugins/local/slater_tc/slater_tc_opt_diag.irp.f b/plugins/local/slater_tc/slater_tc_opt_diag.irp.f index 78f9dc66..3c5a5d12 100644 --- a/plugins/local/slater_tc/slater_tc_opt_diag.irp.f +++ b/plugins/local/slater_tc/slater_tc_opt_diag.irp.f @@ -19,13 +19,13 @@ PROVIDE HF_bitmask PROVIDE mo_l_coef mo_r_coef - call diag_htilde_mu_mat_bi_ortho_slow(N_int, HF_bitmask, hmono, htwoe, htot) + call diag_htc_bi_orth_2e_brute(N_int, HF_bitmask, hmono, htwoe, htot) ref_tc_energy_1e = hmono ref_tc_energy_2e = htwoe if(three_body_h_tc) then - call diag_htilde_three_body_ints_bi_ort_slow(N_int, HF_bitmask, hthree) + call diag_htc_bi_orth_3e_brute(N_int, HF_bitmask, hthree) ref_tc_energy_3e = hthree else ref_tc_energy_3e = 0.d0 @@ -524,3 +524,310 @@ end ! --- +subroutine diag_htc_bi_orth_2e_brute(Nint, key_i, hmono, htwoe, htot) + + BEGIN_DOC + ! + ! diagonal element of htilde ONLY FOR ONE- AND TWO-BODY TERMS + ! + END_DOC + + use bitmasks + + implicit none + integer, intent(in) :: Nint + integer(bit_kind), intent(in) :: key_i(Nint,2) + double precision, intent(out) :: hmono,htwoe,htot + integer :: occ(Nint*bit_kind_size,2) + integer :: Ne(2), i, j, ii, jj, ispin, jspin, k, kk + double precision :: get_mo_two_e_integral_tc_int + integer(bit_kind) :: key_i_core(Nint,2) + + PROVIDE mo_bi_ortho_tc_two_e + + hmono = 0.d0 + htwoe = 0.d0 + htot = 0.d0 + + call bitstring_to_list_ab(key_i, occ, Ne, Nint) + + do ispin = 1, 2 + do i = 1, Ne(ispin) + ii = occ(i,ispin) + hmono += mo_bi_ortho_tc_one_e(ii,ii) + enddo + enddo + + ! alpha/beta two-body + ispin = 1 + jspin = 2 + do i = 1, Ne(ispin) ! electron 1 (so it can be associated to mu(r1)) + ii = occ(i,ispin) + do j = 1, Ne(jspin) ! electron 2 + jj = occ(j,jspin) + htwoe += mo_bi_ortho_tc_two_e(jj,ii,jj,ii) + enddo + enddo + + ! alpha/alpha two-body + do i = 1, Ne(ispin) + ii = occ(i,ispin) + do j = i+1, Ne(ispin) + jj = occ(j,ispin) + htwoe += mo_bi_ortho_tc_two_e(ii,jj,ii,jj) - mo_bi_ortho_tc_two_e(ii,jj,jj,ii) + enddo + enddo + + ! beta/beta two-body + do i = 1, Ne(jspin) + ii = occ(i,jspin) + do j = i+1, Ne(jspin) + jj = occ(j,jspin) + htwoe += mo_bi_ortho_tc_two_e(ii,jj,ii,jj) - mo_bi_ortho_tc_two_e(ii,jj,jj,ii) + enddo + enddo + + htot = hmono + htwoe + +end + +! --- + +subroutine diag_htc_bi_orth_3e_brute(Nint, key_i, hthree) + + BEGIN_DOC + ! diagonal element of htilde ONLY FOR THREE-BODY TERMS WITH BI ORTHONORMAL ORBITALS + END_DOC + + use bitmasks + + implicit none + integer, intent(in) :: Nint + integer(bit_kind), intent(in) :: key_i(Nint,2) + double precision, intent(out) :: hthree + integer :: occ(Nint*bit_kind_size,2) + integer :: Ne(2),i,j,ii,jj,ispin,jspin,m,mm + integer(bit_kind) :: key_i_core(Nint,2) + double precision :: direct_int, exchange_int, ref + double precision, external :: sym_3_e_int_from_6_idx_tensor + double precision, external :: three_e_diag_parrallel_spin + + PROVIDE mo_l_coef mo_r_coef + + if(core_tc_op) then + do i = 1, Nint + key_i_core(i,1) = xor(key_i(i,1), core_bitmask(i,1)) + key_i_core(i,2) = xor(key_i(i,2), core_bitmask(i,2)) + enddo + call bitstring_to_list_ab(key_i_core, occ, Ne, Nint) + else + call bitstring_to_list_ab(key_i, occ, Ne, Nint) + endif + + hthree = 0.d0 + + if((Ne(1)+Ne(2)) .ge. 3) then + + ! alpha/alpha/beta three-body + do i = 1, Ne(1) + ii = occ(i,1) + do j = i+1, Ne(1) + jj = occ(j,1) + do m = 1, Ne(2) + mm = occ(m,2) + !direct_int = three_body_ints_bi_ort(mm,jj,ii,mm,jj,ii) !uses the 6-idx tensor + !exchange_int = three_body_ints_bi_ort(mm,jj,ii,mm,ii,jj) !uses the 6-idx tensor + direct_int = three_e_3_idx_direct_bi_ort(mm,jj,ii) !uses 3-idx tensor + exchange_int = three_e_3_idx_exch12_bi_ort(mm,jj,ii) !uses 3-idx tensor + hthree += direct_int - exchange_int + enddo + enddo + enddo + + ! beta/beta/alpha three-body + do i = 1, Ne(2) + ii = occ(i,2) + do j = i+1, Ne(2) + jj = occ(j,2) + do m = 1, Ne(1) + mm = occ(m,1) + !direct_int = three_body_ints_bi_ort(mm,jj,ii,mm,jj,ii) !uses the 6-idx tensor + !exchange_int = three_body_ints_bi_ort(mm,jj,ii,mm,ii,jj) !uses the 6-idx tensor + direct_int = three_e_3_idx_direct_bi_ort(mm,jj,ii) + exchange_int = three_e_3_idx_exch12_bi_ort(mm,jj,ii) + hthree += direct_int - exchange_int + enddo + enddo + enddo + + ! alpha/alpha/alpha three-body + do i = 1, Ne(1) + ii = occ(i,1) ! 1 + do j = i+1, Ne(1) + jj = occ(j,1) ! 2 + do m = j+1, Ne(1) + mm = occ(m,1) ! 3 + !hthree += sym_3_e_int_from_6_idx_tensor(mm,jj,ii,mm,jj,ii) !uses the 6 idx tensor + hthree += three_e_diag_parrallel_spin(mm,jj,ii) !uses only 3-idx tensors + enddo + enddo + enddo + + ! beta/beta/beta three-body + do i = 1, Ne(2) + ii = occ(i,2) ! 1 + do j = i+1, Ne(2) + jj = occ(j,2) ! 2 + do m = j+1, Ne(2) + mm = occ(m,2) ! 3 + !hthree += sym_3_e_int_from_6_idx_tensor(mm,jj,ii,mm,jj,ii) !uses the 6 idx tensor + hthree += three_e_diag_parrallel_spin(mm,jj,ii) !uses only 3-idx tensors + enddo + enddo + enddo + + endif + +end + + + +BEGIN_PROVIDER [ double precision, three_e_diag_parrallel_spin_prov, (mo_num, mo_num, mo_num)] + + BEGIN_DOC + ! + ! matrix element of the -L three-body operator ON A BI ORTHONORMAL BASIS + ! + ! three_e_diag_parrallel_spin_prov(m,j,i) = All combinations of the form for same spin matrix elements + ! + ! notice the -1 sign: in this way three_e_diag_parrallel_spin_prov can be directly used to compute Slater rules with a + sign + ! + END_DOC + + implicit none + integer :: i, j, m + double precision :: integral, wall1, wall0, three_e_diag_parrallel_spin + + three_e_diag_parrallel_spin_prov = 0.d0 + print *, ' Providing the three_e_diag_parrallel_spin_prov ...' + + integral = three_e_diag_parrallel_spin(1,1,1) ! to provide all stuffs + call wall_time(wall0) + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (i,j,m,integral) & + !$OMP SHARED (mo_num,three_e_diag_parrallel_spin_prov) + !$OMP DO SCHEDULE (dynamic) + do i = 1, mo_num + do j = 1, mo_num + do m = j, mo_num + three_e_diag_parrallel_spin_prov(m,j,i) = three_e_diag_parrallel_spin(m,j,i) + enddo + enddo + enddo + !$OMP END DO + !$OMP END PARALLEL + + do i = 1, mo_num + do j = 1, mo_num + do m = 1, j + three_e_diag_parrallel_spin_prov(m,j,i) = three_e_diag_parrallel_spin_prov(j,m,i) + enddo + enddo + enddo + + call wall_time(wall1) + print *, ' wall time for three_e_diag_parrallel_spin_prov', wall1 - wall0 + +END_PROVIDER + +BEGIN_PROVIDER [ double precision, three_e_single_parrallel_spin_prov, (mo_num, mo_num, mo_num, mo_num)] + + BEGIN_DOC + ! + ! matrix element of the -L three-body operator FOR THE DIRECT TERMS OF SINGLE EXCITATIONS AND BI ORTHO MOs + ! + ! three_e_single_parrallel_spin_prov(m,j,k,i) = All combination of for same spin matrix elements + ! + ! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign + ! + END_DOC + + implicit none + integer :: i, j, k, m + double precision :: integral, wall1, wall0, three_e_single_parrallel_spin + + three_e_single_parrallel_spin_prov = 0.d0 + print *, ' Providing the three_e_single_parrallel_spin_prov ...' + + integral = three_e_single_parrallel_spin(1,1,1,1) + call wall_time(wall0) + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (i,j,k,m,integral) & + !$OMP SHARED (mo_num,three_e_single_parrallel_spin_prov) + !$OMP DO SCHEDULE (dynamic) + do i = 1, mo_num + do k = 1, mo_num + do j = 1, mo_num + do m = 1, mo_num + three_e_single_parrallel_spin_prov(m,j,k,i) = three_e_single_parrallel_spin(m,j,k,i) + enddo + enddo + enddo + enddo + !$OMP END DO + !$OMP END PARALLEL + + call wall_time(wall1) + print *, ' wall time for three_e_single_parrallel_spin_prov', wall1 - wall0 + +END_PROVIDER + + +! --- + +BEGIN_PROVIDER [ double precision, three_e_double_parrallel_spin_prov, (mo_num, mo_num, mo_num, mo_num, mo_num)] + + BEGIN_DOC + ! + ! matrix element of the -L three-body operator FOR THE DIRECT TERMS OF DOUBLE EXCITATIONS AND BI ORTHO MOs + ! + ! three_e_double_parrallel_spin_prov(m,l,j,k,i) = ::: notice that i is the RIGHT MO and k is the LEFT MO + ! + ! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign + END_DOC + + implicit none + integer :: i, j, k, m, l + double precision :: integral, wall1, wall0, three_e_double_parrallel_spin + + three_e_double_parrallel_spin_prov = 0.d0 + print *, ' Providing the three_e_double_parrallel_spin_prov ...' + call wall_time(wall0) + + integral = three_e_double_parrallel_spin(1,1,1,1,1) + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (i,j,k,m,l,integral) & + !$OMP SHARED (mo_num,three_e_double_parrallel_spin_prov) + !$OMP DO SCHEDULE (dynamic) + do i = 1, mo_num + do k = 1, mo_num + do j = 1, mo_num + do l = 1, mo_num + do m = 1, mo_num + three_e_double_parrallel_spin_prov(m,l,j,k,i) = three_e_double_parrallel_spin(m,l,j,k,i) + enddo + enddo + enddo + enddo + enddo + !$OMP END DO + !$OMP END PARALLEL + + call wall_time(wall1) + print *, ' wall time for three_e_double_parrallel_spin_prov', wall1 - wall0 + +END_PROVIDER + diff --git a/plugins/local/slater_tc/symmetrized_3_e_int_prov.irp.f b/plugins/local/slater_tc/symmetrized_3_e_int_prov.irp.f deleted file mode 100644 index e8277a74..00000000 --- a/plugins/local/slater_tc/symmetrized_3_e_int_prov.irp.f +++ /dev/null @@ -1,140 +0,0 @@ - -BEGIN_PROVIDER [ double precision, three_e_diag_parrallel_spin_prov, (mo_num, mo_num, mo_num)] - - BEGIN_DOC - ! - ! matrix element of the -L three-body operator ON A BI ORTHONORMAL BASIS - ! - ! three_e_diag_parrallel_spin_prov(m,j,i) = All combinations of the form for same spin matrix elements - ! - ! notice the -1 sign: in this way three_e_diag_parrallel_spin_prov can be directly used to compute Slater rules with a + sign - ! - END_DOC - - implicit none - integer :: i, j, m - double precision :: integral, wall1, wall0, three_e_diag_parrallel_spin - - three_e_diag_parrallel_spin_prov = 0.d0 - print *, ' Providing the three_e_diag_parrallel_spin_prov ...' - - integral = three_e_diag_parrallel_spin(1,1,1) ! to provide all stuffs - call wall_time(wall0) - !$OMP PARALLEL & - !$OMP DEFAULT (NONE) & - !$OMP PRIVATE (i,j,m,integral) & - !$OMP SHARED (mo_num,three_e_diag_parrallel_spin_prov) - !$OMP DO SCHEDULE (dynamic) - do i = 1, mo_num - do j = 1, mo_num - do m = j, mo_num - three_e_diag_parrallel_spin_prov(m,j,i) = three_e_diag_parrallel_spin(m,j,i) - enddo - enddo - enddo - !$OMP END DO - !$OMP END PARALLEL - - do i = 1, mo_num - do j = 1, mo_num - do m = 1, j - three_e_diag_parrallel_spin_prov(m,j,i) = three_e_diag_parrallel_spin_prov(j,m,i) - enddo - enddo - enddo - - call wall_time(wall1) - print *, ' wall time for three_e_diag_parrallel_spin_prov', wall1 - wall0 - -END_PROVIDER - -BEGIN_PROVIDER [ double precision, three_e_single_parrallel_spin_prov, (mo_num, mo_num, mo_num, mo_num)] - - BEGIN_DOC - ! - ! matrix element of the -L three-body operator FOR THE DIRECT TERMS OF SINGLE EXCITATIONS AND BI ORTHO MOs - ! - ! three_e_single_parrallel_spin_prov(m,j,k,i) = All combination of for same spin matrix elements - ! - ! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign - ! - END_DOC - - implicit none - integer :: i, j, k, m - double precision :: integral, wall1, wall0, three_e_single_parrallel_spin - - three_e_single_parrallel_spin_prov = 0.d0 - print *, ' Providing the three_e_single_parrallel_spin_prov ...' - - integral = three_e_single_parrallel_spin(1,1,1,1) - call wall_time(wall0) - !$OMP PARALLEL & - !$OMP DEFAULT (NONE) & - !$OMP PRIVATE (i,j,k,m,integral) & - !$OMP SHARED (mo_num,three_e_single_parrallel_spin_prov) - !$OMP DO SCHEDULE (dynamic) - do i = 1, mo_num - do k = 1, mo_num - do j = 1, mo_num - do m = 1, mo_num - three_e_single_parrallel_spin_prov(m,j,k,i) = three_e_single_parrallel_spin(m,j,k,i) - enddo - enddo - enddo - enddo - !$OMP END DO - !$OMP END PARALLEL - - call wall_time(wall1) - print *, ' wall time for three_e_single_parrallel_spin_prov', wall1 - wall0 - -END_PROVIDER - - -! --- - -BEGIN_PROVIDER [ double precision, three_e_double_parrallel_spin_prov, (mo_num, mo_num, mo_num, mo_num, mo_num)] - - BEGIN_DOC - ! - ! matrix element of the -L three-body operator FOR THE DIRECT TERMS OF DOUBLE EXCITATIONS AND BI ORTHO MOs - ! - ! three_e_double_parrallel_spin_prov(m,l,j,k,i) = ::: notice that i is the RIGHT MO and k is the LEFT MO - ! - ! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign - END_DOC - - implicit none - integer :: i, j, k, m, l - double precision :: integral, wall1, wall0, three_e_double_parrallel_spin - - three_e_double_parrallel_spin_prov = 0.d0 - print *, ' Providing the three_e_double_parrallel_spin_prov ...' - call wall_time(wall0) - - integral = three_e_double_parrallel_spin(1,1,1,1,1) - !$OMP PARALLEL & - !$OMP DEFAULT (NONE) & - !$OMP PRIVATE (i,j,k,m,l,integral) & - !$OMP SHARED (mo_num,three_e_double_parrallel_spin_prov) - !$OMP DO SCHEDULE (dynamic) - do i = 1, mo_num - do k = 1, mo_num - do j = 1, mo_num - do l = 1, mo_num - do m = 1, mo_num - three_e_double_parrallel_spin_prov(m,l,j,k,i) = three_e_double_parrallel_spin(m,l,j,k,i) - enddo - enddo - enddo - enddo - enddo - !$OMP END DO - !$OMP END PARALLEL - - call wall_time(wall1) - print *, ' wall time for three_e_double_parrallel_spin_prov', wall1 - wall0 - -END_PROVIDER - diff --git a/plugins/local/slater_tc_no_opt/.gitignore b/plugins/local/slater_tc_no_opt/.gitignore new file mode 100644 index 00000000..1561915b --- /dev/null +++ b/plugins/local/slater_tc_no_opt/.gitignore @@ -0,0 +1,59 @@ +IRPF90_temp/ +IRPF90_man/ +build.ninja +irpf90.make +ezfio_interface.irp.f +irpf90_entities +tags +Makefile +ao_basis +ao_one_e_ints +ao_two_e_erf_ints +ao_two_e_ints +aux_quantities +becke_numerical_grid +bitmask +cis +cisd +cipsi +davidson +davidson_dressed +davidson_undressed +density_for_dft +determinants +dft_keywords +dft_utils_in_r +dft_utils_one_e +dft_utils_two_body +dressing +dummy +electrons +ezfio_files +fci +generators_cas +generators_full +hartree_fock +iterations +kohn_sham +kohn_sham_rs +mo_basis +mo_guess +mo_one_e_ints +mo_two_e_erf_ints +mo_two_e_ints +mpi +mrpt_utils +nuclei +perturbation +pseudo +psiref_cas +psiref_utils +scf_utils +selectors_cassd +selectors_full +selectors_utils +single_ref_method +slave +tools +utils +zmq diff --git a/plugins/local/slater_tc_no_opt/NEED b/plugins/local/slater_tc_no_opt/NEED new file mode 100644 index 00000000..a8669866 --- /dev/null +++ b/plugins/local/slater_tc_no_opt/NEED @@ -0,0 +1,8 @@ +determinants +normal_order_old +bi_ort_ints +bi_ortho_mos +tc_keywords +non_hermit_dav +dav_general_mat +tc_scf diff --git a/plugins/local/slater_tc_no_opt/README.rst b/plugins/local/slater_tc_no_opt/README.rst new file mode 100644 index 00000000..90679e4c --- /dev/null +++ b/plugins/local/slater_tc_no_opt/README.rst @@ -0,0 +1,4 @@ +================ +slater_tc_no_opt +================ + diff --git a/plugins/local/slater_tc/h_biortho.irp.f b/plugins/local/slater_tc_no_opt/h_biortho.irp.f similarity index 100% rename from plugins/local/slater_tc/h_biortho.irp.f rename to plugins/local/slater_tc_no_opt/h_biortho.irp.f diff --git a/plugins/local/slater_tc_no_opt/h_mat_triple.irp.f b/plugins/local/slater_tc_no_opt/h_mat_triple.irp.f new file mode 100644 index 00000000..e2c8f982 --- /dev/null +++ b/plugins/local/slater_tc_no_opt/h_mat_triple.irp.f @@ -0,0 +1,193 @@ +subroutine get_excitation_general(key_i,key_j, Nint,degree_array,holes_array, particles_array,phase) + use bitmasks + BEGIN_DOC +! returns the array, for each spin, of holes/particles between key_i and key_j +! +! with the following convention: a^+_{particle} a_{hole}|key_i> = |key_j> + END_DOC + include 'utils/constants.include.F' + implicit none + integer, intent(in) :: Nint + integer(bit_kind), intent(in) :: key_j(Nint,2),key_i(Nint,2) + integer, intent(out) :: holes_array(100,2),particles_array(100,2),degree_array(2) + double precision, intent(out) :: phase + integer :: ispin,k,i,pos + integer(bit_kind) :: key_hole, key_particle + integer(bit_kind) :: xorvec(N_int_max,2) + holes_array = -1 + particles_array = -1 + degree_array = 0 + do i = 1, N_int + xorvec(i,1) = xor( key_i(i,1), key_j(i,1)) + xorvec(i,2) = xor( key_i(i,2), key_j(i,2)) + degree_array(1) += popcnt(xorvec(i,1)) + degree_array(2) += popcnt(xorvec(i,2)) + enddo + degree_array(1) = shiftr(degree_array(1),1) + degree_array(2) = shiftr(degree_array(2),1) + + do ispin = 1, 2 + k = 1 + !!! GETTING THE HOLES + do i = 1, N_int + key_hole = iand(xorvec(i,ispin),key_i(i,ispin)) + do while(key_hole .ne.0_bit_kind) + pos = trailz(key_hole) + holes_array(k,ispin) = 1+ bit_kind_size * (i-1) + pos + key_hole = ibclr(key_hole,pos) + k += 1 + if(k .gt.100)then + print*,'WARNING in get_excitation_general' + print*,'More than a 100-th excitation for spin ',ispin + print*,'stoping ...' + stop + endif + enddo + enddo + enddo + do ispin = 1, 2 + k = 1 + !!! GETTING THE PARTICLES + do i = 1, N_int + key_particle = iand(xor(key_i(i,ispin),key_j(i,ispin)),key_j(i,ispin)) + do while(key_particle .ne.0_bit_kind) + pos = trailz(key_particle) + particles_array(k,ispin) = 1+ bit_kind_size * (i-1) + pos + key_particle = ibclr(key_particle,pos) + k += 1 + if(k .gt.100)then + print*,'WARNING in get_excitation_general ' + print*,'More than a 100-th excitation for spin ',ispin + print*,'stoping ...' + stop + endif + enddo + enddo + enddo + integer :: h,p, i_ok + integer(bit_kind), allocatable :: det_i(:,:),det_ip(:,:) + integer :: exc(0:2,2,2) + double precision :: phase_tmp + allocate(det_i(Nint,2),det_ip(N_int,2)) + det_i = key_i + phase = 1.d0 + do ispin = 1, 2 + do i = 1, degree_array(ispin) + h = holes_array(i,ispin) + p = particles_array(i,ispin) + det_ip = det_i + call do_single_excitation(det_ip,h,p,ispin,i_ok) + if(i_ok == -1)then + print*,'excitation was not possible ' + stop + endif + call get_single_excitation(det_i,det_ip,exc,phase_tmp,Nint) + phase *= phase_tmp + det_i = det_ip + enddo + enddo + +end + +subroutine get_holes_general(key_i, key_j,Nint, holes_array) + use bitmasks + BEGIN_DOC +! returns the array, per spin, of holes between key_i and key_j +! +! with the following convention: a_{hole}|key_i> --> |key_j> + END_DOC + implicit none + integer, intent(in) :: Nint + integer(bit_kind), intent(in) :: key_j(Nint,2),key_i(Nint,2) + integer, intent(out) :: holes_array(100,2) + integer(bit_kind) :: key_hole + integer :: ispin,k,i,pos + holes_array = -1 + do ispin = 1, 2 + k = 1 + do i = 1, N_int + key_hole = iand(xor(key_i(i,ispin),key_j(i,ispin)),key_i(i,ispin)) + do while(key_hole .ne.0_bit_kind) + pos = trailz(key_hole) + holes_array(k,ispin) = 1+ bit_kind_size * (i-1) + pos + key_hole = ibclr(key_hole,pos) + k += 1 + if(k .gt.100)then + print*,'WARNING in get_holes_general' + print*,'More than a 100-th excitation for spin ',ispin + print*,'stoping ...' + stop + endif + enddo + enddo + enddo +end + +subroutine get_particles_general(key_i, key_j,Nint,particles_array) + use bitmasks + BEGIN_DOC +! returns the array, per spin, of particles between key_i and key_j +! +! with the following convention: a^dagger_{particle}|key_i> --> |key_j> + END_DOC + implicit none + integer, intent(in) :: Nint + integer(bit_kind), intent(in) :: key_j(Nint,2),key_i(Nint,2) + integer, intent(out) :: particles_array(100,2) + integer(bit_kind) :: key_particle + integer :: ispin,k,i,pos + particles_array = -1 + do ispin = 1, 2 + k = 1 + do i = 1, N_int + key_particle = iand(xor(key_i(i,ispin),key_j(i,ispin)),key_j(i,ispin)) + do while(key_particle .ne.0_bit_kind) + pos = trailz(key_particle) + particles_array(k,ispin) = 1+ bit_kind_size * (i-1) + pos + key_particle = ibclr(key_particle,pos) + k += 1 + if(k .gt.100)then + print*,'WARNING in get_holes_general' + print*,'More than a 100-th excitation for spin ',ispin + print*,'Those are the two determinants' + call debug_det(key_i, N_int) + call debug_det(key_j, N_int) + print*,'stoping ...' + stop + endif + enddo + enddo + enddo +end + +subroutine get_phase_general(key_i,Nint,degree, holes_array, particles_array,phase) + implicit none + integer, intent(in) :: degree(2), Nint + integer(bit_kind), intent(in) :: key_i(Nint,2) + integer, intent(in) :: holes_array(100,2),particles_array(100,2) + double precision, intent(out) :: phase + integer :: i,ispin,h,p, i_ok + integer(bit_kind), allocatable :: det_i(:,:),det_ip(:,:) + integer :: exc(0:2,2,2) + double precision :: phase_tmp + allocate(det_i(Nint,2),det_ip(N_int,2)) + det_i = key_i + phase = 1.d0 + do ispin = 1, 2 + do i = 1, degree(ispin) + h = holes_array(i,ispin) + p = particles_array(i,ispin) + det_ip = det_i + call do_single_excitation(det_ip,h,p,ispin,i_ok) + if(i_ok == -1)then + print*,'excitation was not possible ' + stop + endif + call get_single_excitation(det_i,det_ip,exc,phase_tmp,Nint) + phase *= phase_tmp + det_i = det_ip + enddo + enddo + +end + diff --git a/plugins/local/slater_tc/h_tc_bi_ortho_psi.irp.f b/plugins/local/slater_tc_no_opt/h_tc_bi_ortho_psi.irp.f similarity index 100% rename from plugins/local/slater_tc/h_tc_bi_ortho_psi.irp.f rename to plugins/local/slater_tc_no_opt/h_tc_bi_ortho_psi.irp.f diff --git a/plugins/local/slater_tc/slater_tc_3e_slow.irp.f b/plugins/local/slater_tc_no_opt/slater_tc_3e_slow.irp.f similarity index 99% rename from plugins/local/slater_tc/slater_tc_3e_slow.irp.f rename to plugins/local/slater_tc_no_opt/slater_tc_3e_slow.irp.f index cb33d343..f7919653 100644 --- a/plugins/local/slater_tc/slater_tc_3e_slow.irp.f +++ b/plugins/local/slater_tc_no_opt/slater_tc_3e_slow.irp.f @@ -1,7 +1,7 @@ ! --- -subroutine diag_htilde_three_body_ints_bi_ort_slow(Nint, key_i, hthree) +subroutine diag_htc_bi_orth_3e_brute(Nint, key_i, hthree) BEGIN_DOC ! diagonal element of htilde ONLY FOR THREE-BODY TERMS WITH BI ORTHONORMAL ORBITALS diff --git a/plugins/local/slater_tc/slater_tc.irp.f b/plugins/local/slater_tc_no_opt/slater_tc_no_opt.irp.f similarity index 82% rename from plugins/local/slater_tc/slater_tc.irp.f rename to plugins/local/slater_tc_no_opt/slater_tc_no_opt.irp.f index 27ab47c5..0fcc587f 100644 --- a/plugins/local/slater_tc/slater_tc.irp.f +++ b/plugins/local/slater_tc_no_opt/slater_tc_no_opt.irp.f @@ -1,4 +1,4 @@ -program slater_tc +program slater_tc_no_opt implicit none BEGIN_DOC ! TODO : Put the documentation of the program here diff --git a/plugins/local/slater_tc/slater_tc_slow.irp.f b/plugins/local/slater_tc_no_opt/slater_tc_slow.irp.f similarity index 80% rename from plugins/local/slater_tc/slater_tc_slow.irp.f rename to plugins/local/slater_tc_no_opt/slater_tc_slow.irp.f index caf7d665..b06fd12f 100644 --- a/plugins/local/slater_tc/slater_tc_slow.irp.f +++ b/plugins/local/slater_tc_no_opt/slater_tc_slow.irp.f @@ -61,7 +61,7 @@ subroutine htilde_mu_mat_bi_ortho_slow(key_j, key_i, Nint, hmono, htwoe, hthree, if(degree.gt.2) return if(degree == 0) then - call diag_htilde_mu_mat_bi_ortho_slow(Nint, key_i, hmono, htwoe, htot) + call diag_htc_bi_orth_2e_brute(Nint, key_i, hmono, htwoe, htot) else if (degree == 1) then call single_htilde_mu_mat_bi_ortho_slow(Nint, key_j, key_i, hmono, htwoe, htot) else if(degree == 2) then @@ -76,7 +76,7 @@ subroutine htilde_mu_mat_bi_ortho_slow(key_j, key_i, Nint, hmono, htwoe, hthree, else if((degree == 1) .and. (elec_num .gt. 2) .and. three_e_4_idx_term) then call single_htilde_three_body_ints_bi_ort_slow(Nint, key_j, key_i, hthree) else if((degree == 0) .and. (elec_num .gt. 2) .and. three_e_3_idx_term) then - call diag_htilde_three_body_ints_bi_ort_slow(Nint, key_i, hthree) + call diag_htc_bi_orth_3e_brute(Nint, key_i, hthree) endif endif @@ -95,75 +95,6 @@ end ! --- -subroutine diag_htilde_mu_mat_bi_ortho_slow(Nint, key_i, hmono, htwoe, htot) - - BEGIN_DOC - ! - ! diagonal element of htilde ONLY FOR ONE- AND TWO-BODY TERMS - ! - END_DOC - - use bitmasks - - implicit none - integer, intent(in) :: Nint - integer(bit_kind), intent(in) :: key_i(Nint,2) - double precision, intent(out) :: hmono,htwoe,htot - integer :: occ(Nint*bit_kind_size,2) - integer :: Ne(2), i, j, ii, jj, ispin, jspin, k, kk - double precision :: get_mo_two_e_integral_tc_int - integer(bit_kind) :: key_i_core(Nint,2) - - PROVIDE mo_bi_ortho_tc_two_e - - hmono = 0.d0 - htwoe = 0.d0 - htot = 0.d0 - - call bitstring_to_list_ab(key_i, occ, Ne, Nint) - - do ispin = 1, 2 - do i = 1, Ne(ispin) - ii = occ(i,ispin) - hmono += mo_bi_ortho_tc_one_e(ii,ii) - enddo - enddo - - ! alpha/beta two-body - ispin = 1 - jspin = 2 - do i = 1, Ne(ispin) ! electron 1 (so it can be associated to mu(r1)) - ii = occ(i,ispin) - do j = 1, Ne(jspin) ! electron 2 - jj = occ(j,jspin) - htwoe += mo_bi_ortho_tc_two_e(jj,ii,jj,ii) - enddo - enddo - - ! alpha/alpha two-body - do i = 1, Ne(ispin) - ii = occ(i,ispin) - do j = i+1, Ne(ispin) - jj = occ(j,ispin) - htwoe += mo_bi_ortho_tc_two_e(ii,jj,ii,jj) - mo_bi_ortho_tc_two_e(ii,jj,jj,ii) - enddo - enddo - - ! beta/beta two-body - do i = 1, Ne(jspin) - ii = occ(i,jspin) - do j = i+1, Ne(jspin) - jj = occ(j,jspin) - htwoe += mo_bi_ortho_tc_two_e(ii,jj,ii,jj) - mo_bi_ortho_tc_two_e(ii,jj,jj,ii) - enddo - enddo - - htot = hmono + htwoe - -end - -! --- - subroutine double_htilde_mu_mat_bi_ortho_slow(Nint, key_j, key_i, hmono, htwoe, htot) BEGIN_DOC diff --git a/src/determinants/slater_rules_general.irp.f b/src/determinants/slater_rules_general.irp.f new file mode 100644 index 00000000..e987c846 --- /dev/null +++ b/src/determinants/slater_rules_general.irp.f @@ -0,0 +1,192 @@ +subroutine get_excitation_general(key_i,key_j, Nint,degree_array,holes_array, particles_array,phase) + use bitmasks + BEGIN_DOC +! returns the array, for each spin, of holes/particles between key_i and key_j +! +! with the following convention: a^+_{particle} a_{hole}|key_i> = |key_j> + END_DOC + include 'utils/constants.include.F' + implicit none + integer, intent(in) :: Nint + integer(bit_kind), intent(in) :: key_j(Nint,2),key_i(Nint,2) + integer, intent(out) :: holes_array(100,2),particles_array(100,2),degree_array(2) + double precision, intent(out) :: phase + integer :: ispin,k,i,pos + integer(bit_kind) :: key_hole, key_particle + integer(bit_kind) :: xorvec(N_int_max,2) + holes_array = -1 + particles_array = -1 + degree_array = 0 + do i = 1, N_int + xorvec(i,1) = xor( key_i(i,1), key_j(i,1)) + xorvec(i,2) = xor( key_i(i,2), key_j(i,2)) + degree_array(1) += popcnt(xorvec(i,1)) + degree_array(2) += popcnt(xorvec(i,2)) + enddo + degree_array(1) = shiftr(degree_array(1),1) + degree_array(2) = shiftr(degree_array(2),1) + + do ispin = 1, 2 + k = 1 + !!! GETTING THE HOLES + do i = 1, N_int + key_hole = iand(xorvec(i,ispin),key_i(i,ispin)) + do while(key_hole .ne.0_bit_kind) + pos = trailz(key_hole) + holes_array(k,ispin) = 1+ bit_kind_size * (i-1) + pos + key_hole = ibclr(key_hole,pos) + k += 1 + if(k .gt.100)then + print*,'WARNING in get_excitation_general' + print*,'More than a 100-th excitation for spin ',ispin + print*,'stoping ...' + stop + endif + enddo + enddo + enddo + do ispin = 1, 2 + k = 1 + !!! GETTING THE PARTICLES + do i = 1, N_int + key_particle = iand(xor(key_i(i,ispin),key_j(i,ispin)),key_j(i,ispin)) + do while(key_particle .ne.0_bit_kind) + pos = trailz(key_particle) + particles_array(k,ispin) = 1+ bit_kind_size * (i-1) + pos + key_particle = ibclr(key_particle,pos) + k += 1 + if(k .gt.100)then + print*,'WARNING in get_excitation_general ' + print*,'More than a 100-th excitation for spin ',ispin + print*,'stoping ...' + stop + endif + enddo + enddo + enddo + integer :: h,p, i_ok + integer(bit_kind), allocatable :: det_i(:,:),det_ip(:,:) + integer :: exc(0:2,2,2) + double precision :: phase_tmp + allocate(det_i(Nint,2),det_ip(N_int,2)) + det_i = key_i + phase = 1.d0 + do ispin = 1, 2 + do i = 1, degree_array(ispin) + h = holes_array(i,ispin) + p = particles_array(i,ispin) + det_ip = det_i + call do_single_excitation(det_ip,h,p,ispin,i_ok) + if(i_ok == -1)then + print*,'excitation was not possible ' + stop + endif + call get_single_excitation(det_i,det_ip,exc,phase_tmp,Nint) + phase *= phase_tmp + det_i = det_ip + enddo + enddo + +end + +subroutine get_holes_general(key_i, key_j,Nint, holes_array) + use bitmasks + BEGIN_DOC +! returns the array, per spin, of holes between key_i and key_j +! +! with the following convention: a_{hole}|key_i> --> |key_j> + END_DOC + implicit none + integer, intent(in) :: Nint + integer(bit_kind), intent(in) :: key_j(Nint,2),key_i(Nint,2) + integer, intent(out) :: holes_array(100,2) + integer(bit_kind) :: key_hole + integer :: ispin,k,i,pos + holes_array = -1 + do ispin = 1, 2 + k = 1 + do i = 1, N_int + key_hole = iand(xor(key_i(i,ispin),key_j(i,ispin)),key_i(i,ispin)) + do while(key_hole .ne.0_bit_kind) + pos = trailz(key_hole) + holes_array(k,ispin) = 1+ bit_kind_size * (i-1) + pos + key_hole = ibclr(key_hole,pos) + k += 1 + if(k .gt.100)then + print*,'WARNING in get_holes_general' + print*,'More than a 100-th excitation for spin ',ispin + print*,'stoping ...' + stop + endif + enddo + enddo + enddo +end + +subroutine get_particles_general(key_i, key_j,Nint,particles_array) + use bitmasks + BEGIN_DOC +! returns the array, per spin, of particles between key_i and key_j +! +! with the following convention: a^dagger_{particle}|key_i> --> |key_j> + END_DOC + implicit none + integer, intent(in) :: Nint + integer(bit_kind), intent(in) :: key_j(Nint,2),key_i(Nint,2) + integer, intent(out) :: particles_array(100,2) + integer(bit_kind) :: key_particle + integer :: ispin,k,i,pos + particles_array = -1 + do ispin = 1, 2 + k = 1 + do i = 1, N_int + key_particle = iand(xor(key_i(i,ispin),key_j(i,ispin)),key_j(i,ispin)) + do while(key_particle .ne.0_bit_kind) + pos = trailz(key_particle) + particles_array(k,ispin) = 1+ bit_kind_size * (i-1) + pos + key_particle = ibclr(key_particle,pos) + k += 1 + if(k .gt.100)then + print*,'WARNING in get_holes_general' + print*,'More than a 100-th excitation for spin ',ispin + print*,'Those are the two determinants' + call debug_det(key_i, N_int) + call debug_det(key_j, N_int) + print*,'stoping ...' + stop + endif + enddo + enddo + enddo +end + +subroutine get_phase_general(key_i,Nint,degree, holes_array, particles_array,phase) + implicit none + integer, intent(in) :: degree(2), Nint + integer(bit_kind), intent(in) :: key_i(Nint,2) + integer, intent(in) :: holes_array(100,2),particles_array(100,2) + double precision, intent(out) :: phase + integer :: i,ispin,h,p, i_ok + integer(bit_kind), allocatable :: det_i(:,:),det_ip(:,:) + integer :: exc(0:2,2,2) + double precision :: phase_tmp + allocate(det_i(Nint,2),det_ip(N_int,2)) + det_i = key_i + phase = 1.d0 + do ispin = 1, 2 + do i = 1, degree(ispin) + h = holes_array(i,ispin) + p = particles_array(i,ispin) + det_ip = det_i + call do_single_excitation(det_ip,h,p,ispin,i_ok) + if(i_ok == -1)then + print*,'excitation was not possible ' + stop + endif + call get_single_excitation(det_i,det_ip,exc,phase_tmp,Nint) + phase *= phase_tmp + det_i = det_ip + enddo + enddo + +end