From 366afb2933baba919db1ad85b7eee965ea56d0c6 Mon Sep 17 00:00:00 2001 From: eginer Date: Mon, 6 May 2024 18:53:20 +0200 Subject: [PATCH] compiling after some cleaning --- plugins/local/old_delta_tc_qmc/NEED | 1 + plugins/local/old_delta_tc_qmc/README.rst | 4 + .../compute_deltamu_right.irp.f | 0 .../dressing_vectors_lr.irp.f | 0 .../old_delta_tc_qmc/old_delta_tc_qmc.irp.f | 7 + plugins/local/slater_tc/h_mat_triple.irp.f | 198 ++++++++++++++++++ .../local/slater_tc_no_opt/h_mat_triple.irp.f | 193 ----------------- .../test_tc_bi_ortho.irp.f | 0 plugins/local/tc_bi_ortho/pt2_tc_cisd.irp.f | 129 ------------ plugins/local/tc_bi_ortho/tc_cisd_sc2.irp.f | 36 ---- .../local/tc_bi_ortho/tc_cisd_sc2_utils.irp.f | 145 ------------- plugins/local/tc_bi_ortho/test_s2_tc.irp.f | 170 --------------- 12 files changed, 210 insertions(+), 673 deletions(-) create mode 100644 plugins/local/old_delta_tc_qmc/NEED create mode 100644 plugins/local/old_delta_tc_qmc/README.rst rename plugins/local/{tc_bi_ortho => old_delta_tc_qmc}/compute_deltamu_right.irp.f (100%) rename plugins/local/{tc_bi_ortho => old_delta_tc_qmc}/dressing_vectors_lr.irp.f (100%) create mode 100644 plugins/local/old_delta_tc_qmc/old_delta_tc_qmc.irp.f create mode 100644 plugins/local/slater_tc/h_mat_triple.irp.f delete mode 100644 plugins/local/slater_tc_no_opt/h_mat_triple.irp.f rename plugins/local/{tc_bi_ortho => slater_tc_no_opt}/test_tc_bi_ortho.irp.f (100%) delete mode 100644 plugins/local/tc_bi_ortho/pt2_tc_cisd.irp.f delete mode 100644 plugins/local/tc_bi_ortho/tc_cisd_sc2.irp.f delete mode 100644 plugins/local/tc_bi_ortho/tc_cisd_sc2_utils.irp.f delete mode 100644 plugins/local/tc_bi_ortho/test_s2_tc.irp.f diff --git a/plugins/local/old_delta_tc_qmc/NEED b/plugins/local/old_delta_tc_qmc/NEED new file mode 100644 index 00000000..8b137891 --- /dev/null +++ b/plugins/local/old_delta_tc_qmc/NEED @@ -0,0 +1 @@ + diff --git a/plugins/local/old_delta_tc_qmc/README.rst b/plugins/local/old_delta_tc_qmc/README.rst new file mode 100644 index 00000000..1d56f96c --- /dev/null +++ b/plugins/local/old_delta_tc_qmc/README.rst @@ -0,0 +1,4 @@ +================ +old_delta_tc_qmc +================ + diff --git a/plugins/local/tc_bi_ortho/compute_deltamu_right.irp.f b/plugins/local/old_delta_tc_qmc/compute_deltamu_right.irp.f similarity index 100% rename from plugins/local/tc_bi_ortho/compute_deltamu_right.irp.f rename to plugins/local/old_delta_tc_qmc/compute_deltamu_right.irp.f diff --git a/plugins/local/tc_bi_ortho/dressing_vectors_lr.irp.f b/plugins/local/old_delta_tc_qmc/dressing_vectors_lr.irp.f similarity index 100% rename from plugins/local/tc_bi_ortho/dressing_vectors_lr.irp.f rename to plugins/local/old_delta_tc_qmc/dressing_vectors_lr.irp.f diff --git a/plugins/local/old_delta_tc_qmc/old_delta_tc_qmc.irp.f b/plugins/local/old_delta_tc_qmc/old_delta_tc_qmc.irp.f new file mode 100644 index 00000000..5ff08bd6 --- /dev/null +++ b/plugins/local/old_delta_tc_qmc/old_delta_tc_qmc.irp.f @@ -0,0 +1,7 @@ +program old_delta_tc_qmc + implicit none + BEGIN_DOC +! TODO : Put the documentation of the program here + END_DOC + print *, 'Hello world' +end diff --git a/plugins/local/slater_tc/h_mat_triple.irp.f b/plugins/local/slater_tc/h_mat_triple.irp.f new file mode 100644 index 00000000..9cb4b60a --- /dev/null +++ b/plugins/local/slater_tc/h_mat_triple.irp.f @@ -0,0 +1,198 @@ +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_no_opt/h_mat_triple.irp.f b/plugins/local/slater_tc_no_opt/h_mat_triple.irp.f deleted file mode 100644 index e2c8f982..00000000 --- a/plugins/local/slater_tc_no_opt/h_mat_triple.irp.f +++ /dev/null @@ -1,193 +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 - diff --git a/plugins/local/tc_bi_ortho/test_tc_bi_ortho.irp.f b/plugins/local/slater_tc_no_opt/test_tc_bi_ortho.irp.f similarity index 100% rename from plugins/local/tc_bi_ortho/test_tc_bi_ortho.irp.f rename to plugins/local/slater_tc_no_opt/test_tc_bi_ortho.irp.f diff --git a/plugins/local/tc_bi_ortho/pt2_tc_cisd.irp.f b/plugins/local/tc_bi_ortho/pt2_tc_cisd.irp.f deleted file mode 100644 index 8940a4f6..00000000 --- a/plugins/local/tc_bi_ortho/pt2_tc_cisd.irp.f +++ /dev/null @@ -1,129 +0,0 @@ -program pt2_tc_cisd - - BEGIN_DOC - ! - ! TODO : Reads psi_det in the EZFIO folder and prints out the left- and right-eigenvectors together - ! with the energy. Saves the left-right wave functions at the end. - ! - END_DOC - - implicit none - - my_grid_becke = .True. - PROVIDE tc_grid1_a tc_grid1_r - my_n_pt_r_grid = tc_grid1_r - my_n_pt_a_grid = tc_grid1_a - touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid - - read_wf = .True. - touch read_wf - - print*, ' nb of states = ', N_states - print*, ' nb of det = ', N_det - call routine_diag() - - call routine -end - -subroutine routine - implicit none - integer :: i,h1,p1,h2,p2,s1,s2,degree - double precision :: h0i,hi0,e00,ei,delta_e - double precision :: norm,e_corr,coef,e_corr_pos,e_corr_neg,e_corr_abs - - integer :: exc(0:2,2,2) - double precision :: phase - double precision :: eh1,ep1,eh2,ep2 - - norm = 0.d0 - e_corr = 0.d0 - e_corr_abs = 0.d0 - e_corr_pos = 0.d0 - e_corr_neg = 0.d0 - call htilde_mu_mat_bi_ortho_tot_slow(psi_det(1,1,1), psi_det(1,1,1), N_int, e00) - do i = 2, N_det - call htilde_mu_mat_bi_ortho_tot_slow(psi_det(1,1,i), psi_det(1,1,1), N_int, hi0) - call htilde_mu_mat_bi_ortho_tot_slow(psi_det(1,1,1), psi_det(1,1,i), N_int, h0i) - call htilde_mu_mat_bi_ortho_tot_slow(psi_det(1,1,i), psi_det(1,1,i), N_int, ei) - call get_excitation_degree(psi_det(1,1,1), psi_det(1,1,i),degree,N_int) - call get_excitation(psi_det(1,1,1), psi_det(1,1,i),exc,degree,phase,N_int) - call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) - eh1 = Fock_matrix_tc_diag_mo_tot(h1) - ep1 = Fock_matrix_tc_diag_mo_tot(p1) - delta_e = eh1 - ep1 - if (degree==2)then - eh2 = Fock_matrix_tc_diag_mo_tot(h2) - ep2 = Fock_matrix_tc_diag_mo_tot(p2) - delta_e += eh2 - ep2 - endif -! delta_e = e00 - ei - coef = hi0/delta_e - norm += coef*coef - e_corr = coef* h0i - if(e_corr.lt.0.d0)then - e_corr_neg += e_corr - elseif(e_corr.gt.0.d0)then - e_corr_pos += e_corr - endif - e_corr_abs += dabs(e_corr) - enddo - print*,'e_corr_abs = ',e_corr_abs - print*,'e_corr_pos = ',e_corr_pos - print*,'e_corr_neg = ',e_corr_neg - print*,'norm = ',dsqrt(norm) - -end - -subroutine routine_diag() - - implicit none - integer :: i, j, k - double precision :: dE - - ! provide eigval_right_tc_bi_orth - ! provide overlap_bi_ortho - ! provide htilde_matrix_elmt_bi_ortho - - if(N_states .eq. 1) then - - print*,'eigval_right_tc_bi_orth = ',eigval_right_tc_bi_orth(1) - print*,'e_tc_left_right = ',e_tc_left_right - print*,'e_tilde_bi_orth_00 = ',e_tilde_bi_orth_00 - print*,'e_pt2_tc_bi_orth = ',e_pt2_tc_bi_orth - print*,'e_pt2_tc_bi_orth_single = ',e_pt2_tc_bi_orth_single - print*,'e_pt2_tc_bi_orth_double = ',e_pt2_tc_bi_orth_double - print*,'***' - print*,'e_corr_bi_orth = ',e_corr_bi_orth - print*,'e_corr_bi_orth_proj = ',e_corr_bi_orth_proj - print*,'e_corr_bi_orth_proj_abs = ',e_corr_bi_orth_proj_abs - print*,'e_corr_single_bi_orth = ',e_corr_single_bi_orth - print*,'e_corr_double_bi_orth = ',e_corr_double_bi_orth - print*,'e_corr_single_bi_orth_abs = ',e_corr_single_bi_orth_abs - print*,'e_corr_double_bi_orth_abs = ',e_corr_double_bi_orth_abs - print*,'Left/right eigenvectors' - do i = 1,N_det - write(*,'(I5,X,(100(F12.7,X)))')i,leigvec_tc_bi_orth(i,1),reigvec_tc_bi_orth(i,1),leigvec_tc_bi_orth(i,1)*reigvec_tc_bi_orth(i,1) - enddo - - else - - print*,'eigval_right_tc_bi_orth : ' - do i = 1, N_states - print*, i, eigval_right_tc_bi_orth(i) - enddo - - print*,'' - print*,'******************************************************' - print*,'TC Excitation energies (au) (eV)' - do i = 2, N_states - dE = eigval_right_tc_bi_orth(i) - eigval_right_tc_bi_orth(1) - print*, i, dE, dE/0.0367502d0 - enddo - print*,'' - - endif - -end - - - diff --git a/plugins/local/tc_bi_ortho/tc_cisd_sc2.irp.f b/plugins/local/tc_bi_ortho/tc_cisd_sc2.irp.f deleted file mode 100644 index d4c8c55d..00000000 --- a/plugins/local/tc_bi_ortho/tc_cisd_sc2.irp.f +++ /dev/null @@ -1,36 +0,0 @@ - -! --- - -program tc_cisd_sc2 - - BEGIN_DOC - ! TODO : Put the documentation of the program here - END_DOC - - implicit none - - print *, 'Hello world' - - my_grid_becke = .True. - PROVIDE tc_grid1_a tc_grid1_r - my_n_pt_r_grid = tc_grid1_r - my_n_pt_a_grid = tc_grid1_a - touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid - - read_wf = .True. - touch read_wf - - call test - -end - -! --- - -subroutine test() - implicit none -! double precision, allocatable :: dressing_dets(:),e_corr_dets(:) -! allocate(dressing_dets(N_det),e_corr_dets(N_det)) -! e_corr_dets = 0.d0 -! call get_cisd_sc2_dressing(psi_det,e_corr_dets,N_det,dressing_dets) - provide eigval_tc_cisd_sc2_bi_ortho -end diff --git a/plugins/local/tc_bi_ortho/tc_cisd_sc2_utils.irp.f b/plugins/local/tc_bi_ortho/tc_cisd_sc2_utils.irp.f deleted file mode 100644 index 5cbf26d2..00000000 --- a/plugins/local/tc_bi_ortho/tc_cisd_sc2_utils.irp.f +++ /dev/null @@ -1,145 +0,0 @@ - BEGIN_PROVIDER [ double precision, reigvec_tc_cisd_sc2_bi_ortho, (N_det,N_states)] -&BEGIN_PROVIDER [ double precision, leigvec_tc_cisd_sc2_bi_ortho, (N_det,N_states)] -&BEGIN_PROVIDER [ double precision, eigval_tc_cisd_sc2_bi_ortho, (N_states)] - implicit none - integer :: it,n_real,degree,i,istate - double precision :: e_before, e_current,thr, hmono,htwoe,hthree,accu - double precision, allocatable :: e_corr_dets(:),h0j(:), h_sc2(:,:), dressing_dets(:) - double precision, allocatable :: leigvec_tc_bi_orth_tmp(:,:),reigvec_tc_bi_orth_tmp(:,:),eigval_right_tmp(:) - allocate(leigvec_tc_bi_orth_tmp(N_det,N_det),reigvec_tc_bi_orth_tmp(N_det,N_det),eigval_right_tmp(N_det)) - allocate(e_corr_dets(N_det),h0j(N_det),h_sc2(N_det,N_det),dressing_dets(N_det)) - allocate(H_jj(N_det),vec_tmp(N_det,n_states_diag),eigval_tmp(N_states)) - dressing_dets = 0.d0 - do i = 1, N_det - call htilde_mu_mat_bi_ortho_tot_slow(psi_det(1,1,i), psi_det(1,1,i), N_int, H_jj(i)) - call get_excitation_degree(HF_bitmask,psi_det(1,1,i),degree,N_int) - if(degree == 1 .or. degree == 2)then - call htilde_mu_mat_opt_bi_ortho(HF_bitmask,psi_det(1,1,i),N_int,hmono,htwoe,hthree,h0j(i)) - endif - enddo - reigvec_tc_bi_orth_tmp = 0.d0 - do i = 1, N_det - reigvec_tc_bi_orth_tmp(i,1) = psi_r_coef_bi_ortho(i,1) - enddo - vec_tmp = 0.d0 - do istate = 1, N_states - vec_tmp(:,istate) = reigvec_tc_bi_orth_tmp(:,istate) - enddo - do istate = N_states+1, n_states_diag - vec_tmp(istate,istate) = 1.d0 - enddo - print*,'Diagonalizing the TC CISD ' - call davidson_general_diag_dressed_ext_rout_nonsym_b1space(vec_tmp, H_jj, dressing_dets,eigval_tmp, N_det, n_states, n_states_diag, converged, htc_bi_ortho_calc_tdav_slow) - do i = 1, N_det - e_corr_dets(i) = reigvec_tc_bi_orth_tmp(i,1) * h0j(i)/reigvec_tc_bi_orth_tmp(1,1) - enddo - E_before = eigval_tmp(1) - print*,'Starting from ',E_before - - e_current = 10.d0 - thr = 1.d-5 - it = 0 - dressing_dets = 0.d0 - double precision, allocatable :: H_jj(:),vec_tmp(:,:),eigval_tmp(:) - external htc_bi_ortho_calc_tdav_slow - external htcdag_bi_ortho_calc_tdav_slow - logical :: converged - do while (dabs(E_before-E_current).gt.thr) - it += 1 - E_before = E_current -! h_sc2 = htilde_matrix_elmt_bi_ortho - call get_cisd_sc2_dressing(psi_det,e_corr_dets,N_det,dressing_dets) - do i = 1, N_det -! print*,'dressing_dets(i) = ',dressing_dets(i) - h_sc2(i,i) += dressing_dets(i) - enddo - print*,'********************' - print*,'iteration ',it -! call non_hrmt_real_diag(N_det,h_sc2,& -! leigvec_tc_bi_orth_tmp,reigvec_tc_bi_orth_tmp,& -! n_real,eigval_right_tmp) -! print*,'eigval_right_tmp(1)',eigval_right_tmp(1) - vec_tmp = 0.d0 - do istate = 1, N_states - vec_tmp(:,istate) = reigvec_tc_bi_orth_tmp(:,istate) - enddo - do istate = N_states+1, n_states_diag - vec_tmp(istate,istate) = 1.d0 - enddo - call davidson_general_diag_dressed_ext_rout_nonsym_b1space(vec_tmp, H_jj, dressing_dets,eigval_tmp, N_det, n_states, n_states_diag, converged, htc_bi_ortho_calc_tdav_slow) - print*,'outside Davidson' - print*,'eigval_tmp(1) = ',eigval_tmp(1) - do i = 1, N_det - reigvec_tc_bi_orth_tmp(i,1) = vec_tmp(i,1) - e_corr_dets(i) = reigvec_tc_bi_orth_tmp(i,1) * h0j(i)/reigvec_tc_bi_orth_tmp(1,1) - enddo -! E_current = eigval_right_tmp(1) - E_current = eigval_tmp(1) - print*,'it, E(SC)^2 = ',it,E_current - enddo - eigval_tc_cisd_sc2_bi_ortho(1:N_states) = eigval_right_tmp(1:N_states) - reigvec_tc_cisd_sc2_bi_ortho(1:N_det,1:N_states) = reigvec_tc_bi_orth_tmp(1:N_det,1:N_states) - leigvec_tc_cisd_sc2_bi_ortho(1:N_det,1:N_states) = leigvec_tc_bi_orth_tmp(1:N_det,1:N_states) - -END_PROVIDER - -subroutine get_cisd_sc2_dressing(dets,e_corr_dets,ndet,dressing_dets) - implicit none - use bitmasks - integer, intent(in) :: ndet - integer(bit_kind), intent(in) :: dets(N_int,2,ndet) - double precision, intent(in) :: e_corr_dets(ndet) - double precision, intent(out) :: dressing_dets(ndet) - integer, allocatable :: degree(:),hole(:,:),part(:,:),spin(:,:) - integer(bit_kind), allocatable :: hole_part(:,:,:) - integer :: i,j,k, exc(0:2,2,2),h1,p1,h2,p2,s1,s2 - integer(bit_kind) :: xorvec(2,N_int) - - double precision :: phase - dressing_dets = 0.d0 - allocate(degree(ndet),hole(2,ndet),part(2,ndet), spin(2,ndet),hole_part(N_int,2,ndet)) - do i = 2, ndet - call get_excitation_degree(HF_bitmask,dets(1,1,i),degree(i),N_int) - do j = 1, N_int - hole_part(j,1,i) = xor( HF_bitmask(j,1), dets(j,1,i)) - hole_part(j,2,i) = xor( HF_bitmask(j,2), dets(j,2,i)) - enddo - if(degree(i) == 1)then - call get_single_excitation(HF_bitmask,psi_det(1,1,i),exc,phase,N_int) - else if(degree(i) == 2)then - call get_double_excitation(HF_bitmask,psi_det(1,1,i),exc,phase,N_int) - endif - call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) - hole(1,i) = h1 - hole(2,i) = h2 - part(1,i) = p1 - part(2,i) = p2 - spin(1,i) = s1 - spin(2,i) = s2 - enddo - - integer :: same - if(elec_alpha_num+elec_beta_num<3)return - do i = 2, ndet - do j = i+1, ndet - same = 0 - if(degree(i) == degree(j) .and. degree(i)==1)cycle - do k = 1, N_int - xorvec(k,1) = iand(hole_part(k,1,i),hole_part(k,1,j)) - xorvec(k,2) = iand(hole_part(k,2,i),hole_part(k,2,j)) - same += popcnt(xorvec(k,1)) + popcnt(xorvec(k,2)) - enddo -! print*,'i,j',i,j -! call debug_det(dets(1,1,i),N_int) -! call debug_det(hole_part(1,1,i),N_int) -! call debug_det(dets(1,1,j),N_int) -! call debug_det(hole_part(1,1,j),N_int) -! print*,'same = ',same - if(same.eq.0)then - dressing_dets(i) += e_corr_dets(j) - dressing_dets(j) += e_corr_dets(i) - endif - enddo - enddo - -end diff --git a/plugins/local/tc_bi_ortho/test_s2_tc.irp.f b/plugins/local/tc_bi_ortho/test_s2_tc.irp.f deleted file mode 100644 index 7c70b119..00000000 --- a/plugins/local/tc_bi_ortho/test_s2_tc.irp.f +++ /dev/null @@ -1,170 +0,0 @@ - -! --- - -program test_tc - - implicit none - - my_grid_becke = .True. - PROVIDE tc_grid1_a tc_grid1_r - my_n_pt_r_grid = tc_grid1_r - my_n_pt_a_grid = tc_grid1_a - touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid - - read_wf = .True. - touch read_wf - - call provide_all_three_ints_bi_ortho() - call routine_h_triple_left - call routine_h_triple_right -! call routine_test_s2_davidson - -end - -subroutine routine_h_triple_right - implicit none - logical :: do_right - integer :: sze ,i, N_st, j - double precision :: sij, accu_e, accu_s, accu_e_0, accu_s_0 - double precision, allocatable :: v_0_ref(:,:),u_0(:,:),s_0_ref(:,:) - double precision, allocatable :: v_0_new(:,:),s_0_new(:,:) - sze = N_det - N_st = 1 - allocate(v_0_ref(N_det,1),u_0(N_det,1),s_0_ref(N_det,1),s_0_new(N_det,1),v_0_new(N_det,1)) - print*,'Checking first the Right ' - do i = 1, sze - u_0(i,1) = psi_r_coef_bi_ortho(i,1) - enddo - double precision :: wall0,wall1 - call wall_time(wall0) - call H_tc_s2_u_0_with_pure_three_omp(v_0_ref,s_0_ref, u_0,N_st,sze) - call wall_time(wall1) - print*,'time for omp',wall1 - wall0 - call wall_time(wall0) - call H_tc_s2_u_0_with_pure_three(v_0_new, s_0_new, u_0, N_st, sze) - call wall_time(wall1) - print*,'time serial ',wall1 - wall0 - accu_e = 0.d0 - accu_s = 0.d0 - do i = 1, sze - accu_e += dabs(v_0_ref(i,1) - v_0_new(i,1)) - accu_s += dabs(s_0_ref(i,1) - s_0_new(i,1)) - enddo - print*,'accu_e = ',accu_e - print*,'accu_s = ',accu_s - -end - -subroutine routine_h_triple_left - implicit none - logical :: do_right - integer :: sze ,i, N_st, j - double precision :: sij, accu_e, accu_s, accu_e_0, accu_s_0 - double precision, allocatable :: v_0_ref(:,:),u_0(:,:),s_0_ref(:,:) - double precision, allocatable :: v_0_new(:,:),s_0_new(:,:) - sze = N_det - N_st = 1 - allocate(v_0_ref(N_det,1),u_0(N_det,1),s_0_ref(N_det,1),s_0_new(N_det,1),v_0_new(N_det,1)) - print*,'Checking the Left ' - do i = 1, sze - u_0(i,1) = psi_l_coef_bi_ortho(i,1) - enddo - double precision :: wall0,wall1 - call wall_time(wall0) - call H_tc_s2_dagger_u_0_with_pure_three_omp(v_0_ref,s_0_ref, u_0,N_st,sze) - call wall_time(wall1) - print*,'time for omp',wall1 - wall0 - call wall_time(wall0) - call H_tc_s2_dagger_u_0_with_pure_three(v_0_new, s_0_new, u_0, N_st, sze) - call wall_time(wall1) - print*,'time serial ',wall1 - wall0 - accu_e = 0.d0 - accu_s = 0.d0 - do i = 1, sze - accu_e += dabs(v_0_ref(i,1) - v_0_new(i,1)) - accu_s += dabs(s_0_ref(i,1) - s_0_new(i,1)) - enddo - print*,'accu_e = ',accu_e - print*,'accu_s = ',accu_s - -end - - -subroutine routine_test_s2_davidson - implicit none - double precision, allocatable :: H_jj(:),vec_tmp(:,:), energies(:) , s2(:) - integer :: i,istate - logical :: converged - external H_tc_s2_dagger_u_0_opt - external H_tc_s2_u_0_opt - allocate(H_jj(N_det),vec_tmp(N_det,n_states_diag),energies(n_states_diag), s2(n_states_diag)) - do i = 1, N_det - call htilde_mu_mat_bi_ortho_tot_slow(psi_det(1,1,i), psi_det(1,1,i), N_int, H_jj(i)) - enddo - ! Preparing the left-eigenvector - print*,'Computing the left-eigenvector ' - vec_tmp = 0.d0 - do istate = 1, N_states - vec_tmp(1:N_det,istate) = psi_l_coef_bi_ortho(1:N_det,istate) - enddo - do istate = N_states+1, n_states_diag - vec_tmp(istate,istate) = 1.d0 - enddo - do istate = 1, N_states - leigvec_tc_bi_orth(1:N_det,istate) = vec_tmp(1:N_det,istate) - enddo - integer :: n_it_max - n_it_max = 1 - call davidson_hs2_nonsym_b1space(vec_tmp, H_jj, s2, energies, N_det, n_states, n_states_diag, n_it_max, converged, H_tc_s2_dagger_u_0_opt) - double precision, allocatable :: v_0_new(:,:),s_0_new(:,:) - integer :: sze,N_st - logical :: do_right - sze = N_det - N_st = 1 - do_right = .False. - allocate(s_0_new(N_det,1),v_0_new(N_det,1)) - call H_tc_s2_u_0_nstates_openmp(v_0_new,s_0_new,vec_tmp,N_st,sze, do_right) - double precision :: accu_e_0, accu_s_0 - accu_e_0 = 0.d0 - accu_s_0 = 0.d0 - do i = 1, sze - accu_e_0 += v_0_new(i,1) * vec_tmp(i,1) - accu_s_0 += s_0_new(i,1) * vec_tmp(i,1) - enddo - print*,'energies = ',energies - print*,'s2 = ',s2 - print*,'accu_e_0',accu_e_0 - print*,'accu_s_0',accu_s_0 - - ! Preparing the right-eigenvector - print*,'Computing the right-eigenvector ' - vec_tmp = 0.d0 - do istate = 1, N_states - vec_tmp(1:N_det,istate) = psi_r_coef_bi_ortho(1:N_det,istate) - enddo - do istate = N_states+1, n_states_diag - vec_tmp(istate,istate) = 1.d0 - enddo - do istate = 1, N_states - leigvec_tc_bi_orth(1:N_det,istate) = vec_tmp(1:N_det,istate) - enddo - n_it_max = 1 - call davidson_hs2_nonsym_b1space(vec_tmp, H_jj, s2, energies, N_det, n_states, n_states_diag, n_it_max, converged, H_tc_s2_u_0_opt) - sze = N_det - N_st = 1 - do_right = .True. - v_0_new = 0.d0 - s_0_new = 0.d0 - call H_tc_s2_u_0_nstates_openmp(v_0_new,s_0_new,vec_tmp,N_st,sze, do_right) - accu_e_0 = 0.d0 - accu_s_0 = 0.d0 - do i = 1, sze - accu_e_0 += v_0_new(i,1) * vec_tmp(i,1) - accu_s_0 += s_0_new(i,1) * vec_tmp(i,1) - enddo - print*,'energies = ',energies - print*,'s2 = ',s2 - print*,'accu_e_0',accu_e_0 - print*,'accu_s_0',accu_s_0 - -end