From af60a02919b9923772398abf877e79dd0eed160d Mon Sep 17 00:00:00 2001 From: eginer Date: Tue, 24 Jan 2023 13:49:04 +0100 Subject: [PATCH 1/2] it works with same left and right coefficients for 2 dets in WF --- src/cipsi_tc_bi_ortho/get_d.irp.f | 612 +++++++++++------- .../pt2_stoch_routines.irp.f | 2 +- .../run_selection_slave.irp.f | 2 +- src/cipsi_tc_bi_ortho/selection.irp.f | 97 ++- src/cipsi_tc_bi_ortho/stochastic_cipsi.irp.f | 1 + src/determinants/fock_diag.irp.f | 34 +- src/fci_tc_bi/selectors.irp.f | 10 +- src/tc_bi_ortho/slater_tc_opt.irp.f | 42 ++ src/tc_bi_ortho/slater_tc_opt_diag.irp.f | 194 ++++++ src/tc_bi_ortho/slater_tc_opt_double.irp.f | 55 ++ src/tc_bi_ortho/slater_tc_opt_single.irp.f | 112 ++++ 11 files changed, 890 insertions(+), 271 deletions(-) diff --git a/src/cipsi_tc_bi_ortho/get_d.irp.f b/src/cipsi_tc_bi_ortho/get_d.irp.f index c642f420..58b1972a 100644 --- a/src/cipsi_tc_bi_ortho/get_d.irp.f +++ b/src/cipsi_tc_bi_ortho/get_d.irp.f @@ -194,7 +194,7 @@ end subroutine get_d3_h ! --- -subroutine get_d2(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs) +subroutine get_d2(gen, phasemask, bannedOrb, banned, mat_p, mat_m, mask, h, p, sp, coefs) use bitmasks implicit none @@ -202,15 +202,15 @@ subroutine get_d2(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs) integer(bit_kind), intent(in) :: mask(N_int, 2), gen(N_int, 2) integer(bit_kind), intent(in) :: phasemask(N_int,2) logical, intent(in) :: bannedOrb(mo_num, 2), banned(mo_num, mo_num,2) - double precision, intent(in) :: coefs(N_states) - double precision, intent(inout) :: mat(N_states, mo_num, mo_num) + double precision, intent(in) :: coefs(N_states,2) + double precision, intent(inout) :: mat_p(N_states, mo_num, mo_num), mat_m(N_states, mo_num, mo_num) integer, intent(in) :: h(0:2,2), p(0:4,2), sp - double precision, external :: get_phase_bi, mo_two_e_integral + double precision, external :: get_phase_bi integer :: i, j, k, tip, ma, mi, puti, putj integer :: h1, h2, p1, p2, i1, i2 - double precision :: hij, phase + double precision :: hij, hji, phase integer, parameter:: turn2d(2,3,4) = reshape((/0,0, 0,0, 0,0, 3,4, 0,0, 0,0, 2,4, 1,4, 0,0, 2,3, 1,3, 1,2 /), (/2,3,4/)) integer, parameter :: turn2(2) = (/2, 1/) @@ -222,11 +222,13 @@ subroutine get_d2(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs) tip = p(0,1) * p(0,2) ma = sp + print*,'in get d2' if(p(0,1) > p(0,2)) ma = 1 if(p(0,1) < p(0,2)) ma = 2 mi = mod(ma, 2) + 1 if(sp == 3) then + print*,'in sp == 3' if(ma == 2) bant = 2 if(tip == 3) then @@ -247,20 +249,24 @@ subroutine get_d2(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs) ! ! - ! < p2 p1 | H^tilde^dag| h1 h2 > = < h1 h2 | w_ee^h + t^nh | p1 p2 > - hij = mo_two_e_integral(p1, p2, h1, h2) - mo_two_e_integral(p2, p1, h1, h2) + hji = mo_bi_ortho_tc_two_e(p1, p2, h1, h2) - mo_bi_ortho_tc_two_e(p2, p1, h1, h2) + hij = mo_bi_ortho_tc_two_e(h1, h2, p1, p2) - mo_bi_ortho_tc_two_e(h2, h1, p1, p2) if (hij == 0.d0) cycle hij = hij * get_phase_bi(phasemask, ma, ma, h1, p1, h2, p2, N_int) + hji = hji * get_phase_bi(phasemask, ma, ma, h1, p1, h2, p2, N_int) if(ma == 1) then !DIR$ LOOP COUNT AVG(4) do k=1,N_states - mat(k, putj, puti) = mat(k, putj, puti) + coefs(k) * hij + mat_p(k, putj, puti) = mat_p(k, putj, puti) + coefs(k,1) * hij + mat_m(k, putj, puti) = mat_m(k, putj, puti) + coefs(k,2) * hji enddo else !DIR$ LOOP COUNT AVG(4) do k=1,N_states - mat(k, puti, putj) = mat(k, puti, putj) + coefs(k) * hij + mat_p(k, puti, putj) = mat_p(k, puti, putj) + coefs(k,1) * hij + mat_m(k, puti, putj) = mat_m(k, puti, putj) + coefs(k,2) * hji enddo end if end do @@ -277,12 +283,15 @@ subroutine get_d2(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs) if(banned(puti,putj,bant) .or. bannedOrb(puti,1)) cycle p1 = p(turn2(i), 1) - hij = mo_two_e_integral(p1, p2, h1, h2) + hji = mo_bi_ortho_tc_two_e(p1, p2, h1, h2) + hij = mo_bi_ortho_tc_two_e(h1, h2, p1, p2) if (hij /= 0.d0) then hij = hij * get_phase_bi(phasemask, 1, 2, h1, p1, h2, p2, N_int) + hji = hji * get_phase_bi(phasemask, 1, 2, h1, p1, h2, p2, N_int) !DIR$ LOOP COUNT AVG(4) do k=1,N_states - mat(k, puti, putj) = mat(k, puti, putj) + coefs(k) * hij + mat_p(k, puti, putj) = mat_p(k, puti, putj) + coefs(k,1) * hij + mat_m(k, puti, putj) = mat_m(k, puti, putj) + coefs(k,2) * hji enddo endif end do @@ -290,6 +299,7 @@ subroutine get_d2(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs) end if else + print*,'NOT in sp == 3' if(tip == 0) then h1 = h(1, ma) h2 = h(2, ma) @@ -305,13 +315,16 @@ subroutine get_d2(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs) i2 = turn2d(2, i, j) p1 = p(i1, ma) p2 = p(i2, ma) - hij = mo_two_e_integral(p1, p2, h1, h2) - mo_two_e_integral(p2,p1, h1, h2) + hij = mo_bi_ortho_tc_two_e(p1, p2, h1, h2) - mo_bi_ortho_tc_two_e(p2,p1, h1, h2) + hji = mo_bi_ortho_tc_two_e(h1, h2, p1, p2) - mo_bi_ortho_tc_two_e(h2,h1, p1, p2) if (hij == 0.d0) cycle hij = hij * get_phase_bi(phasemask, ma, ma, h1, p1, h2, p2, N_int) + hji = hji * get_phase_bi(phasemask, ma, ma, h1, p1, h2, p2, N_int) !DIR$ LOOP COUNT AVG(4) do k=1,N_states - mat(k, puti, putj) = mat(k, puti, putj) +coefs(k) * hij + mat_p(k, puti, putj) = mat_p(k, puti, putj) +coefs(k,1) * hij + mat_m(k, puti, putj) = mat_m(k, puti, putj) +coefs(k,2) * hji enddo end do end do @@ -327,19 +340,23 @@ subroutine get_d2(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs) if(banned(puti,putj,1)) cycle p2 = p(i, ma) - hij = mo_two_e_integral(p1, p2, h1, h2) + hij = mo_bi_ortho_tc_two_e(p1, p2, h1, h2) + hji = mo_bi_ortho_tc_two_e(h1, h2, p1, p2) if (hij == 0.d0) cycle hij = hij * get_phase_bi(phasemask, mi, ma, h1, p1, h2, p2, N_int) + hji = hji * get_phase_bi(phasemask, mi, ma, h1, p1, h2, p2, N_int) if (puti < putj) then !DIR$ LOOP COUNT AVG(4) do k=1,N_states - mat(k, puti, putj) = mat(k, puti, putj) + coefs(k) * hij + mat_p(k, puti, putj) = mat_p(k, puti, putj) + coefs(k,1) * hij + mat_m(k, puti, putj) = mat_m(k, puti, putj) + coefs(k,2) * hji enddo else !DIR$ LOOP COUNT AVG(4) do k=1,N_states - mat(k, putj, puti) = mat(k, putj, puti) + coefs(k) * hij + mat_p(k, putj, puti) = mat_p(k, putj, puti) + coefs(k,1) * hij + mat_m(k, putj, puti) = mat_m(k, putj, puti) + coefs(k,2) * hji enddo endif end do @@ -351,12 +368,15 @@ subroutine get_d2(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs) p2 = p(2, mi) h1 = h(1, mi) h2 = h(2, mi) - hij = (mo_two_e_integral(p1, p2, h1, h2) - mo_two_e_integral(p2,p1, h1, h2)) + hij = (mo_bi_ortho_tc_two_e(p1, p2, h1, h2) - mo_bi_ortho_tc_two_e(p2,p1, h1, h2)) + hji = (mo_bi_ortho_tc_two_e(h1, h2, p1, p2) - mo_bi_ortho_tc_two_e(h2,h1, p1, p2)) if (hij /= 0.d0) then hij = hij * get_phase_bi(phasemask, mi, mi, h1, p1, h2, p2, N_int) + hji = hji * get_phase_bi(phasemask, mi, mi, h1, p1, h2, p2, N_int) !DIR$ LOOP COUNT AVG(4) do k=1,N_states - mat(k, puti, putj) = mat(k, puti, putj) + coefs(k) * hij + mat_p(k, puti, putj) = mat_p(k, puti, putj) + coefs(k,1) * hij + mat_m(k, puti, putj) = mat_m(k, puti, putj) + coefs(k,2) * hji enddo end if end if @@ -367,7 +387,7 @@ end subroutine get_d2 ! --- -subroutine get_d1(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs) +subroutine get_d1(gen, phasemask, bannedOrb, banned, mat_p, mat_m, mask, h, p, sp, coefs) use bitmasks implicit none @@ -376,28 +396,34 @@ subroutine get_d1(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs) integer(bit_kind), intent(in) :: phasemask(N_int,2) logical, intent(in) :: bannedOrb(mo_num, 2), banned(mo_num, mo_num,2) integer(bit_kind) :: det(N_int, 2) - double precision, intent(in) :: coefs(N_states) - double precision, intent(inout) :: mat(N_states, mo_num, mo_num) + double precision, intent(in) :: coefs(N_states,2) + double precision, intent(inout) :: mat_p(N_states, mo_num, mo_num), mat_m(N_states, mo_num, mo_num) integer, intent(in) :: h(0:2,2), p(0:4,2), sp - double precision, external :: get_phase_bi, mo_two_e_integral + double precision, external :: get_phase_bi logical :: ok logical, allocatable :: lbanned(:,:) integer :: puti, putj, ma, mi, s1, s2, i, i1, i2, j integer :: hfix, pfix, h1, h2, p1, p2, ib, k, l + integer :: mm integer, parameter :: turn2(2) = (/2,1/) integer, parameter :: turn3(2,3) = reshape((/2,3, 1,3, 1,2/), (/2,3/)) integer :: bant - double precision, allocatable :: hij_cache(:,:) - double precision :: hij, tmp_row(N_states, mo_num), tmp_row2(N_states, mo_num) + double precision, allocatable :: hij_cache(:,:) + double precision :: hij, tmp_row_ij(N_states, mo_num), tmp_row_ij2(N_states, mo_num) + double precision, allocatable :: hji_cache(:,:) + double precision :: hji, tmp_row_ji(N_states, mo_num), tmp_row_ji2(N_states, mo_num) PROVIDE mo_integrals_map N_int allocate (lbanned(mo_num, 2)) allocate (hij_cache(mo_num,2)) + allocate (hji_cache(mo_num,2)) lbanned = bannedOrb + print*,'in get d1' + call debug_det(gen, N_int) do i=1, p(0,1) lbanned(p(i,1), 1) = .true. @@ -413,16 +439,26 @@ subroutine get_d1(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs) bant = 1 if(sp == 3) then + print*,'in sp == 3' !move MA if(ma == 2) bant = 2 puti = p(1,mi) hfix = h(1,ma) p1 = p(1,ma) p2 = p(2,ma) + print*,puti, hfix,p1,p2 if(.not. bannedOrb(puti, mi)) then - call get_mo_two_e_integrals(hfix,p1,p2,mo_num,hij_cache(1,1),mo_integrals_map) - call get_mo_two_e_integrals(hfix,p2,p1,mo_num,hij_cache(1,2),mo_integrals_map) - tmp_row = 0d0 +! print*,'not banned' + do mm = 1, mo_num + hji_cache(mm,1) = mo_bi_ortho_tc_two_e(p1,p2,mm,hfix) + hji_cache(mm,2) = mo_bi_ortho_tc_two_e(p2,p1,mm,hfix) + hij_cache(mm,1) = mo_bi_ortho_tc_two_e(mm,hfix,p1,p2) + hij_cache(mm,2) = mo_bi_ortho_tc_two_e(mm,hfix,p2,p1) + enddo +! call get_mo_bi_ortho_tc_two_es(hfix,p1,p2,mo_num,hij_cache(1,1),mo_integrals_map) +! call get_mo_bi_ortho_tc_two_es(hfix,p2,p1,mo_num,hij_cache(1,2),mo_integrals_map) + tmp_row_ij = 0d0 + tmp_row_ji = 0d0 do putj=1, hfix-1 if(lbanned(putj, ma)) cycle if(banned(putj, puti,bant)) cycle @@ -431,7 +467,15 @@ subroutine get_d1(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs) hij = hij * get_phase_bi(phasemask, ma, ma, putj, p1, hfix, p2, N_int) !DIR$ LOOP COUNT AVG(4) do k=1,N_states - tmp_row(k,putj) = tmp_row(k,putj) + hij * coefs(k) + tmp_row_ij(k,putj) = tmp_row_ij(k,putj) + hij * coefs(k,2) + enddo + endif + hji = hji_cache(putj,1) - hji_cache(putj,2) + if (hji /= 0.d0) then + hji = hji * get_phase_bi(phasemask, ma, ma, putj, p1, hfix, p2, N_int) + !DIR$ LOOP COUNT AVG(4) + do k=1,N_states + tmp_row_ji(k,putj) = tmp_row_ji(k,putj) + hji * coefs(k,1) enddo endif end do @@ -443,18 +487,28 @@ subroutine get_d1(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs) hij = hij * get_phase_bi(phasemask, ma, ma, hfix, p1, putj, p2, N_int) !DIR$ LOOP COUNT AVG(4) do k=1,N_states - tmp_row(k,putj) = tmp_row(k,putj) + hij * coefs(k) + tmp_row_ij(k,putj) = tmp_row_ij(k,putj) + hij * coefs(k,2) + enddo + endif + hji = hji_cache(putj,2) - hji_cache(putj,1) + if (hji /= 0.d0) then + hji = hji * get_phase_bi(phasemask, ma, ma, hfix, p1, putj, p2, N_int) + !DIR$ LOOP COUNT AVG(4) + do k=1,N_states + tmp_row_ji(k,putj) = tmp_row_ji(k,putj) + hji * coefs(k,1) enddo endif end do if(ma == 1) then - mat(1:N_states,1:mo_num,puti) = mat(1:N_states,1:mo_num,puti) + tmp_row(1:N_states,1:mo_num) + mat_p(1:N_states,1:mo_num,puti) = mat_p(1:N_states,1:mo_num,puti) + tmp_row_ij(1:N_states,1:mo_num) + mat_m(1:N_states,1:mo_num,puti) = mat_m(1:N_states,1:mo_num,puti) + tmp_row_ji(1:N_states,1:mo_num) else do l=1,mo_num !DIR$ LOOP COUNT AVG(4) do k=1,N_states - mat(k,puti,l) = mat(k,puti,l) + tmp_row(k,l) + mat_p(k,puti,l) = mat_p(k,puti,l) + tmp_row_ij(k,l) + mat_m(k,puti,l) = mat_m(k,puti,l) + tmp_row_ji(k,l) enddo enddo end if @@ -462,10 +516,18 @@ subroutine get_d1(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs) !MOVE MI pfix = p(1,mi) - tmp_row = 0d0 - tmp_row2 = 0d0 - call get_mo_two_e_integrals(hfix,pfix,p1,mo_num,hij_cache(1,1),mo_integrals_map) - call get_mo_two_e_integrals(hfix,pfix,p2,mo_num,hij_cache(1,2),mo_integrals_map) + tmp_row_ij = 0d0 + tmp_row_ij2 = 0d0 + tmp_row_ji = 0d0 + tmp_row_ji2 = 0d0 +! call get_mo_bi_ortho_tc_two_es(hfix,pfix,p1,mo_num,hij_cache(1,1),mo_integrals_map) +! call get_mo_bi_ortho_tc_two_es(hfix,pfix,p2,mo_num,hij_cache(1,2),mo_integrals_map) + do mm = 1, mo_num + hji_cache(mm,1) = mo_bi_ortho_tc_two_e(pfix,p1,mm,hfix) + hji_cache(mm,2) = mo_bi_ortho_tc_two_e(pfix,p2,mm,hfix) + hij_cache(mm,1) = mo_bi_ortho_tc_two_e(mm,hfix,pfix,p1) + hij_cache(mm,2) = mo_bi_ortho_tc_two_e(mm,hfix,pfix,p2) + enddo putj = p1 do puti = 1, mo_num !HOT @@ -478,7 +540,15 @@ subroutine get_d1(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs) hij = hij * get_phase_bi(phasemask, ma, mi, hfix, p2, puti, pfix, N_int) !DIR$ LOOP COUNT AVG(4) do k=1,N_states - tmp_row(k,puti) = tmp_row(k,puti) + hij * coefs(k) + tmp_row_ij(k,puti) = tmp_row_ij(k,puti) + hij * coefs(k,2) + enddo + endif + hji = hji_cache(puti,2) + if (hji /= 0.d0) then + hji = hji * get_phase_bi(phasemask, ma, mi, hfix, p2, puti, pfix, N_int) + !DIR$ LOOP COUNT AVG(4) + do k=1,N_states + tmp_row_ji(k,puti) = tmp_row_ji(k,puti) + hji * coefs(k,1) enddo endif endif @@ -489,7 +559,14 @@ subroutine get_d1(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs) if (hij /= 0.d0) then hij = hij * get_phase_bi(phasemask, ma, mi, hfix, p1, puti, pfix, N_int) do k=1,N_states - tmp_row2(k,puti) = tmp_row2(k,puti) + hij * coefs(k) + tmp_row_ij2(k,puti) = tmp_row_ij2(k,puti) + hij * coefs(k,2) + enddo + endif + hji = hji_cache(puti,1) + if (hji /= 0.d0) then + hji = hji * get_phase_bi(phasemask, ma, mi, hfix, p1, puti, pfix, N_int) + do k=1,N_states + tmp_row_ji2(k,puti) = tmp_row_ji2(k,puti) + hji * coefs(k,1) enddo endif endif @@ -497,19 +574,24 @@ subroutine get_d1(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs) enddo if(mi == 1) then - mat(:,:,p1) = mat(:,:,p1) + tmp_row(:,:) - mat(:,:,p2) = mat(:,:,p2) + tmp_row2(:,:) + mat_p(:,:,p1) = mat_p(:,:,p1) + tmp_row_ij(:,:) + mat_p(:,:,p2) = mat_p(:,:,p2) + tmp_row_ij2(:,:) + mat_m(:,:,p1) = mat_m(:,:,p1) + tmp_row_ji(:,:) + mat_m(:,:,p2) = mat_m(:,:,p2) + tmp_row_ji2(:,:) else do l=1,mo_num !DIR$ LOOP COUNT AVG(4) do k=1,N_states - mat(k,p1,l) = mat(k,p1,l) + tmp_row(k,l) - mat(k,p2,l) = mat(k,p2,l) + tmp_row2(k,l) + mat_p(k,p1,l) = mat_p(k,p1,l) + tmp_row_ij(k,l) + mat_p(k,p2,l) = mat_p(k,p2,l) + tmp_row_ij2(k,l) + mat_m(k,p1,l) = mat_m(k,p1,l) + tmp_row_ji(k,l) + mat_m(k,p2,l) = mat_m(k,p2,l) + tmp_row_ji2(k,l) enddo enddo end if else ! sp /= 3 + print*,'not in sp == 3' if(p(0,ma) == 3) then do i=1,3 @@ -517,16 +599,28 @@ subroutine get_d1(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs) puti = p(i, ma) p1 = p(turn3(1,i), ma) p2 = p(turn3(2,i), ma) - call get_mo_two_e_integrals(hfix,p1,p2,mo_num,hij_cache(1,1),mo_integrals_map) - call get_mo_two_e_integrals(hfix,p2,p1,mo_num,hij_cache(1,2),mo_integrals_map) - tmp_row = 0d0 +! call get_mo_bi_ortho_tc_two_es(hfix,p1,p2,mo_num,hij_cache(1,1),mo_integrals_map) +! call get_mo_bi_ortho_tc_two_es(hfix,p2,p1,mo_num,hij_cache(1,2),mo_integrals_map) + do mm = 1, mo_num + hji_cache(mm,1) = mo_bi_ortho_tc_two_e(p1,p2,mm,hfix) + hji_cache(mm,2) = mo_bi_ortho_tc_two_e(p2,p1,mm,hfix) + hij_cache(mm,1) = mo_bi_ortho_tc_two_e(mm,hfix,p1,p2) + hij_cache(mm,2) = mo_bi_ortho_tc_two_e(mm,hfix,p2,p1) + enddo + tmp_row_ij = 0d0 + tmp_row_ji = 0d0 do putj=1,hfix-1 if(banned(putj,puti,1)) cycle if(lbanned(putj,ma)) cycle hij = hij_cache(putj,1) - hij_cache(putj,2) if (hij /= 0.d0) then hij = hij * get_phase_bi(phasemask, ma, ma, putj, p1, hfix, p2, N_int) - tmp_row(:,putj) = tmp_row(:,putj) + hij * coefs(:) + tmp_row_ij(:,putj) = tmp_row_ij(:,putj) + hij * coefs(:,1) + endif + hji = hji_cache(putj,1) - hji_cache(putj,2) + if (hji /= 0.d0) then + hji = hji * get_phase_bi(phasemask, ma, ma, putj, p1, hfix, p2, N_int) + tmp_row_ji(:,putj) = tmp_row_ji(:,putj) + hji * coefs(:,2) endif end do do putj=hfix+1,mo_num @@ -535,15 +629,22 @@ subroutine get_d1(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs) hij = hij_cache(putj,2) - hij_cache(putj,1) if (hij /= 0.d0) then hij = hij * get_phase_bi(phasemask, ma, ma, hfix, p1, putj, p2, N_int) - tmp_row(:,putj) = tmp_row(:,putj) + hij * coefs(:) + tmp_row_ij(:,putj) = tmp_row_ij(:,putj) + hij * coefs(:,1) + endif + hji = hji_cache(putj,2) - hji_cache(putj,1) + if (hji /= 0.d0) then + hji = hji * get_phase_bi(phasemask, ma, ma, hfix, p1, putj, p2, N_int) + tmp_row_ji(:,putj) = tmp_row_ji(:,putj) + hji * coefs(:,2) endif end do - mat(:, :puti-1, puti) = mat(:, :puti-1, puti) + tmp_row(:,:puti-1) + mat_p(:, :puti-1, puti) = mat_p(:, :puti-1, puti) + tmp_row_ij(:,:puti-1) + mat_m(:, :puti-1, puti) = mat_m(:, :puti-1, puti) + tmp_row_ji(:,:puti-1) do l=puti,mo_num !DIR$ LOOP COUNT AVG(4) do k=1,N_states - mat(k, puti, l) = mat(k, puti,l) + tmp_row(k,l) + mat_p(k, puti, l) = mat_p(k, puti,l) + tmp_row_ij(k,l) + mat_m(k, puti, l) = mat_m(k, puti,l) + tmp_row_ji(k,l) enddo enddo end do @@ -552,10 +653,18 @@ subroutine get_d1(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs) pfix = p(1,mi) p1 = p(1,ma) p2 = p(2,ma) - tmp_row = 0d0 - tmp_row2 = 0d0 - call get_mo_two_e_integrals(hfix,p1,pfix,mo_num,hij_cache(1,1),mo_integrals_map) - call get_mo_two_e_integrals(hfix,p2,pfix,mo_num,hij_cache(1,2),mo_integrals_map) + tmp_row_ij = 0d0 + tmp_row_ij2 = 0d0 + tmp_row_ji = 0d0 + tmp_row_ji2 = 0d0 +! call get_mo_bi_ortho_tc_two_es(hfix,p1,pfix,mo_num,hij_cache(1,1),mo_integrals_map) +! call get_mo_bi_ortho_tc_two_es(hfix,p2,pfix,mo_num,hij_cache(1,2),mo_integrals_map) + do mm = 1, mo_num + hji_cache(mm,1) = mo_bi_ortho_tc_two_e(p1,pfix,mm,hfix) + hji_cache(mm,2) = mo_bi_ortho_tc_two_e(p2,pfix,mm,hfix) + hij_cache(mm,1) = mo_bi_ortho_tc_two_e(mm,hfix,p1,pfix) + hij_cache(mm,2) = mo_bi_ortho_tc_two_e(mm,hfix,p2,pfix) + enddo putj = p2 do puti=1,mo_num if(lbanned(puti,ma)) cycle @@ -566,7 +675,15 @@ subroutine get_d1(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs) hij = hij * get_phase_bi(phasemask, mi, ma, hfix, pfix, puti, p1, N_int) !DIR$ LOOP COUNT AVG(4) do k=1,N_states - tmp_row(k,puti) = tmp_row(k,puti) + hij * coefs(k) + tmp_row_ij(k,puti) = tmp_row_ij(k,puti) + hij * coefs(k,1) + enddo + endif + hji = hji_cache(puti,1) + if (hji /= 0.d0) then + hji = hji * get_phase_bi(phasemask, mi, ma, hfix, pfix, puti, p1, N_int) + !DIR$ LOOP COUNT AVG(4) + do k=1,N_states + tmp_row_ji(k,puti) = tmp_row_ji(k,puti) + hji * coefs(k,2) enddo endif end if @@ -577,23 +694,34 @@ subroutine get_d1(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs) if (hij /= 0.d0) then hij = hij * get_phase_bi(phasemask, mi, ma, hfix, pfix, puti, p2, N_int) do k=1,N_states - tmp_row2(k,puti) = tmp_row2(k,puti) + hij * coefs(k) + tmp_row_ij2(k,puti) = tmp_row_ij2(k,puti) + hij * coefs(k,1) + enddo + endif + hji = hji_cache(puti,2) + if (hji /= 0.d0) then + hji = hji * get_phase_bi(phasemask, mi, ma, hfix, pfix, puti, p2, N_int) + do k=1,N_states + tmp_row_ji2(k,puti) = tmp_row_ji2(k,puti) + hji * coefs(k,2) enddo endif end if end do - mat(:,:p2-1,p2) = mat(:,:p2-1,p2) + tmp_row(:,:p2-1) + mat_p(:,:p2-1,p2) = mat_p(:,:p2-1,p2) + tmp_row_ij(:,:p2-1) + mat_m(:,:p2-1,p2) = mat_m(:,:p2-1,p2) + tmp_row_ji(:,:p2-1) do l=p2,mo_num !DIR$ LOOP COUNT AVG(4) do k=1,N_states - mat(k,p2,l) = mat(k,p2,l) + tmp_row(k,l) + mat_p(k,p2,l) = mat_p(k,p2,l) + tmp_row_ij(k,l) + mat_m(k,p2,l) = mat_m(k,p2,l) + tmp_row_ji(k,l) enddo enddo - mat(:,:p1-1,p1) = mat(:,:p1-1,p1) + tmp_row2(:,:p1-1) + mat_p(:,:p1-1,p1) = mat_p(:,:p1-1,p1) + tmp_row_ij2(:,:p1-1) + mat_m(:,:p1-1,p1) = mat_m(:,:p1-1,p1) + tmp_row_ji2(:,:p1-1) do l=p1,mo_num !DIR$ LOOP COUNT AVG(4) do k=1,N_states - mat(k,p1,l) = mat(k,p1,l) + tmp_row2(k,l) + mat_p(k,p1,l) = mat_p(k,p1,l) + tmp_row_ij2(k,l) + mat_m(k,p1,l) = mat_m(k,p1,l) + tmp_row_ji2(k,l) enddo enddo end if @@ -617,10 +745,17 @@ subroutine get_d1(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs) p2 = p(i2,s2) if(bannedOrb(p1, s1) .or. bannedOrb(p2, s2) .or. banned(p1, p2, 1)) cycle call apply_particles(mask, s1, p1, s2, p2, det, ok, N_int) - call i_h_j(gen, det, N_int, hij) +! call i_h_j(gen, det, N_int, hij) + !!!! GUESS ON THE ORDER OF DETS + print*,'compute hij' +! hij = 0.d0 +! hji = 0.d0 + call htilde_mu_mat_opt_bi_ortho_no_3e(gen,det,N_int, hji) + call htilde_mu_mat_opt_bi_ortho_no_3e(det,gen,N_int, hij) !DIR$ LOOP COUNT AVG(4) do k = 1, N_states - mat(k, p1, p2) = mat(k, p1, p2) + coefs(k) * hij + mat_p(k, p1, p2) = mat_p(k, p1, p2) + coefs(k,1) * hij + mat_m(k, p1, p2) = mat_m(k, p1, p2) + coefs(k,2) * hji enddo enddo enddo @@ -629,7 +764,7 @@ end subroutine get_d1 ! --- -subroutine get_d0(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs) +subroutine get_d0(gen, phasemask, bannedOrb, banned, mat_p, mat_m, mask, h, p, sp, coefs) use bitmasks implicit none @@ -638,72 +773,103 @@ subroutine get_d0(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs) integer(bit_kind), intent(in) :: phasemask(N_int,2) logical, intent(in) :: bannedOrb(mo_num, 2), banned(mo_num, mo_num,2) integer(bit_kind) :: det(N_int, 2) - double precision, intent(in) :: coefs(N_states) - double precision, intent(inout) :: mat(N_states, mo_num, mo_num) + double precision, intent(in) :: coefs(N_states,2) + double precision, intent(inout) :: mat_m(N_states, mo_num, mo_num) + double precision, intent(inout) :: mat_p(N_states, mo_num, mo_num) integer, intent(in) :: h(0:2,2), p(0:4,2), sp - integer :: i, j, k, s, h1, h2, p1, p2, puti, putj - double precision :: hij, phase - double precision, external :: get_phase_bi, mo_two_e_integral + integer :: i, j, k, s, h1, h2, p1, p2, puti, putj, mm + double precision :: hij, phase, hji + double precision, external :: get_phase_bi logical :: ok integer, parameter :: bant=1 double precision, allocatable :: hij_cache1(:), hij_cache2(:) allocate (hij_cache1(mo_num),hij_cache2(mo_num)) + double precision, allocatable :: hji_cache1(:), hji_cache2(:) + allocate (hji_cache1(mo_num),hji_cache2(mo_num)) + print*,'in get d0' +! call debug_det(gen, N_int) if(sp == 3) then ! AB h1 = p(1,1) h2 = p(1,2) +! print*,'in AB' do p1=1, mo_num if(bannedOrb(p1, 1)) cycle - call get_mo_two_e_integrals(p1,h2,h1,mo_num,hij_cache1,mo_integrals_map) +! call get_mo_bi_ortho_tc_two_es(p1,h2,h1,mo_num,hij_cache1,mo_integrals_map) + do mm =1, mo_num + hij_cache1(mm) = mo_bi_ortho_tc_two_e(mm,p1,h2,h1) + hji_cache1(mm) = mo_bi_ortho_tc_two_e(h2,h1,mm,p1) + enddo do p2=1, mo_num if(bannedOrb(p2,2)) cycle if(banned(p1, p2, bant)) cycle ! rentable? if(p1 == h1 .or. p2 == h2) then +! print*,'in p1 == h1 or p2 == h2' call apply_particles(mask, 1,p1,2,p2, det, ok, N_int) - call i_h_j(gen, det, N_int, hij) +! call i_h_j(gen, det, N_int, hij) + !!! GUESS ON THE ORDER + call htilde_mu_mat_opt_bi_ortho_no_3e(det,gen,N_int, hij) + call htilde_mu_mat_opt_bi_ortho_no_3e(gen,det,N_int, hji) else +! print*,'ELSE ' phase = get_phase_bi(phasemask, 1, 2, h1, p1, h2, p2, N_int) hij = hij_cache1(p2) * phase + hji = hji_cache1(p2) * phase end if if (hij == 0.d0) cycle !DIR$ LOOP COUNT AVG(4) do k=1,N_states - mat(k, p1, p2) = mat(k, p1, p2) + coefs(k) * hij ! HOTSPOT + mat_p(k, p1, p2) = mat_p(k, p1, p2) + coefs(k,1) * hij ! HOTSPOT + mat_m(k, p1, p2) = mat_m(k, p1, p2) + coefs(k,2) * hji ! HOTSPOT enddo end do end do else ! AA BB +! print*, 'in AA BB' p1 = p(1,sp) p2 = p(2,sp) do puti=1, mo_num if(bannedOrb(puti, sp)) cycle - call get_mo_two_e_integrals(puti,p2,p1,mo_num,hij_cache1,mo_integrals_map) - call get_mo_two_e_integrals(puti,p1,p2,mo_num,hij_cache2,mo_integrals_map) + do mm = 1, mo_num + hij_cache1(mm) = mo_bi_ortho_tc_two_e(p2,p1,mm,puti) + hij_cache2(mm) = mo_bi_ortho_tc_two_e(p1,p2,mm,puti) + hji_cache1(mm) = mo_bi_ortho_tc_two_e(mm,puti,p2,p1) + hji_cache2(mm) = mo_bi_ortho_tc_two_e(mm,puti,p1,p2) + enddo +! call get_mo_bi_ortho_tc_two_es(puti,p2,p1,mo_num,hij_cache1,mo_integrals_map) +! call get_mo_bi_ortho_tc_two_es(puti,p1,p2,mo_num,hij_cache2,mo_integrals_map) do putj=puti+1, mo_num if(bannedOrb(putj, sp)) cycle if(banned(puti, putj, bant)) cycle ! rentable? if(puti == p1 .or. putj == p2 .or. puti == p2 .or. putj == p1) then call apply_particles(mask, sp,puti,sp,putj, det, ok, N_int) - call i_h_j(gen, det, N_int, hij) - if (hij == 0.d0) cycle +! call i_h_j(gen, det, N_int, hij) + !!! GUESS + call htilde_mu_mat_opt_bi_ortho_no_3e(gen,det,N_int, hij) + call htilde_mu_mat_opt_bi_ortho_no_3e(det,gen,N_int, hji) + if (hij == 0.d0.or.hji == 0.d0) cycle else - hij = (mo_two_e_integral(p1, p2, puti, putj) - mo_two_e_integral(p2, p1, puti, putj)) - if (hij == 0.d0) cycle + hji = (mo_bi_ortho_tc_two_e(p1, p2, puti, putj) - mo_bi_ortho_tc_two_e(p2, p1, puti, putj)) + hij = (mo_bi_ortho_tc_two_e(puti, putj, p1, p2) - mo_bi_ortho_tc_two_e(puti, putj, p2, p1)) + if (hij == 0.d0.or.hji==0.d0) cycle hij = hij * get_phase_bi(phasemask, sp, sp, puti, p1 , putj, p2, N_int) + hji = hji * get_phase_bi(phasemask, sp, sp, puti, p1 , putj, p2, N_int) end if !DIR$ LOOP COUNT AVG(4) do k=1,N_states - mat(k, puti, putj) = mat(k, puti, putj) + coefs(k) * hij + mat_p(k, puti, putj) = mat_p(k, puti, putj) + coefs(k,1) * hij + mat_m(k, puti, putj) = mat_m(k, puti, putj) + coefs(k,2) * hji enddo end do end do end if - deallocate(hij_cache1,hij_cache2) +! deallocate(hij_cache1,hij_cache2) +! deallocate(hji_cache1,hji_cache2) end subroutine get_d0 @@ -734,7 +900,7 @@ end subroutine get_d0 ! double precision :: hij_p, hij_m, phase ! ! double precision, external :: get_phase_bi -! double precision, external :: get_mo_two_e_integral_tc_int, get_mo_two_e_integral_tcdag_int +! double precision, external :: get_mo_bi_ortho_tc_two_e_tc_int, get_mo_bi_ortho_tc_two_e_tcdag_int ! ! PROVIDE mo_integrals_tc_int_map mo_integrals_tcdag_int_map ! @@ -763,10 +929,10 @@ end subroutine get_d0 ! p1 = p(i1, ma) ! p2 = p(i2, ma) ! -! hij_p = get_mo_two_e_integral_tc_int (p1, p2, h1, h2, mo_integrals_tc_int_map ) & -! - get_mo_two_e_integral_tc_int (p2, p1, h1, h2, mo_integrals_tc_int_map ) -! hij_m = get_mo_two_e_integral_tcdag_int(p1, p2, h1, h2, mo_integrals_tcdag_int_map) & -! - get_mo_two_e_integral_tcdag_int(p2, p1, h1, h2, mo_integrals_tcdag_int_map) +! hij_p = get_mo_bi_ortho_tc_two_e_tc_int (p1, p2, h1, h2, mo_integrals_tc_int_map ) & +! - get_mo_bi_ortho_tc_two_e_tc_int (p2, p1, h1, h2, mo_integrals_tc_int_map ) +! hij_m = get_mo_bi_ortho_tc_two_e_tcdag_int(p1, p2, h1, h2, mo_integrals_tcdag_int_map) & +! - get_mo_bi_ortho_tc_two_e_tcdag_int(p2, p1, h1, h2, mo_integrals_tcdag_int_map) ! ! if( (hij_p.eq.0.d0) .and. (hij_m.eq.0.d0) ) cycle ! @@ -802,8 +968,8 @@ end subroutine get_d0 ! if(banned(puti,putj,bant) .or. bannedOrb(puti,1)) cycle ! p1 = p(turn2(i), 1) ! -! hij_p = get_mo_two_e_integral_tc_int (p1, p2, h1, h2, mo_integrals_tc_int_map ) -! hij_m = get_mo_two_e_integral_tcdag_int(p1, p2, h1, h2, mo_integrals_tcdag_int_map) +! hij_p = get_mo_bi_ortho_tc_two_e_tc_int (p1, p2, h1, h2, mo_integrals_tc_int_map ) +! hij_m = get_mo_bi_ortho_tc_two_e_tcdag_int(p1, p2, h1, h2, mo_integrals_tcdag_int_map) ! ! if( (hij_p.ne.0.d0) .and. (hij_m.ne.0.d0) ) then ! hij_p = hij_p * get_phase_bi(phasemask, 1, 2, h1, p1, h2, p2, N_int) @@ -835,10 +1001,10 @@ end subroutine get_d0 ! p1 = p(i1, ma) ! p2 = p(i2, ma) ! -! hij_p = get_mo_two_e_integral_tc_int (p1, p2, h1, h2, mo_integrals_tc_int_map ) & -! - get_mo_two_e_integral_tc_int (p2, p1, h1, h2, mo_integrals_tc_int_map ) -! hij_m = get_mo_two_e_integral_tcdag_int(p1, p2, h1, h2, mo_integrals_tcdag_int_map) & -! - get_mo_two_e_integral_tcdag_int(p2, p1, h1, h2, mo_integrals_tcdag_int_map) +! hij_p = get_mo_bi_ortho_tc_two_e_tc_int (p1, p2, h1, h2, mo_integrals_tc_int_map ) & +! - get_mo_bi_ortho_tc_two_e_tc_int (p2, p1, h1, h2, mo_integrals_tc_int_map ) +! hij_m = get_mo_bi_ortho_tc_two_e_tcdag_int(p1, p2, h1, h2, mo_integrals_tcdag_int_map) & +! - get_mo_bi_ortho_tc_two_e_tcdag_int(p2, p1, h1, h2, mo_integrals_tcdag_int_map) ! ! if( (hij_p.eq.0.d0) .and. (hij_m.eq.0.d0) ) cycle ! @@ -865,8 +1031,8 @@ end subroutine get_d0 ! if(banned(puti,putj,1)) cycle ! p2 = p(i, ma) ! -! hij_p = get_mo_two_e_integral_tc_int (p1, p2, h1, h2, mo_integrals_tc_int_map ) -! hij_m = get_mo_two_e_integral_tcdag_int(p1, p2, h1, h2, mo_integrals_tcdag_int_map) +! hij_p = get_mo_bi_ortho_tc_two_e_tc_int (p1, p2, h1, h2, mo_integrals_tc_int_map ) +! hij_m = get_mo_bi_ortho_tc_two_e_tcdag_int(p1, p2, h1, h2, mo_integrals_tcdag_int_map) ! ! if( (hij_p.eq.0.d0) .and. (hij_m.eq.0.d0) ) cycle ! @@ -895,10 +1061,10 @@ end subroutine get_d0 ! h1 = h(1, mi) ! h2 = h(2, mi) ! -! hij_p = get_mo_two_e_integral_tc_int (p1, p2, h1, h2, mo_integrals_tc_int_map ) & -! - get_mo_two_e_integral_tc_int (p2, p1, h1, h2, mo_integrals_tc_int_map ) -! hij_m = get_mo_two_e_integral_tcdag_int(p1, p2, h1, h2, mo_integrals_tcdag_int_map) & -! - get_mo_two_e_integral_tcdag_int(p2, p1, h1, h2, mo_integrals_tcdag_int_map) +! hij_p = get_mo_bi_ortho_tc_two_e_tc_int (p1, p2, h1, h2, mo_integrals_tc_int_map ) & +! - get_mo_bi_ortho_tc_two_e_tc_int (p2, p1, h1, h2, mo_integrals_tc_int_map ) +! hij_m = get_mo_bi_ortho_tc_two_e_tcdag_int(p1, p2, h1, h2, mo_integrals_tcdag_int_map) & +! - get_mo_bi_ortho_tc_two_e_tcdag_int(p2, p1, h1, h2, mo_integrals_tcdag_int_map) ! ! if( (hij_p.ne.0.d0) .and. (hij_m.ne.0.d0) ) then ! hij_p = hij_p * get_phase_bi(phasemask, mi, mi, h1, p1, h2, p2, N_int) @@ -937,15 +1103,15 @@ end subroutine get_d0 ! double precision, intent(inout) :: mat_p(N_states, mo_num, mo_num), mat_m(N_states, mo_num, mo_num) ! ! double precision, external :: get_phase_bi -! double precision, external :: get_mo_two_e_integral_tc_int, get_mo_two_e_integral_tcdag_int +! double precision, external :: get_mo_bi_ortho_tc_two_e_tc_int, get_mo_bi_ortho_tc_two_e_tcdag_int ! ! logical :: ok ! logical, allocatable :: lbanned(:,:) ! integer :: bant ! integer :: puti, putj, ma, mi, s1, s2, i, i1, i2, j ! integer :: hfix, pfix, h1, h2, p1, p2, ib, k, l -! double precision :: tmp_row_p (N_states, mo_num), tmp_row_m (N_states, mo_num) -! double precision :: hij_p, hij_m, tmp_row2_p(N_states, mo_num), tmp_row2_m(N_states, mo_num) +! double precision :: tmp_row_ij_p (N_states, mo_num), tmp_row_ij_m (N_states, mo_num) +! double precision :: hij_p, hij_m, tmp_row_ij2_p(N_states, mo_num), tmp_row_ij2_m(N_states, mo_num) ! double precision, allocatable :: hijp_cache(:,:), hijm_cache(:,:) ! ! integer, parameter :: turn2(2) = (/2,1/) @@ -979,13 +1145,13 @@ end subroutine get_d0 ! p2 = p(2,ma) ! if(.not. bannedOrb(puti, mi)) then ! -! call get_mo_two_e_integrals_tc_int (hfix, p1, p2, mo_num, hijp_cache(1,1), mo_integrals_tc_int_map ) -! call get_mo_two_e_integrals_tc_int (hfix, p2, p1, mo_num, hijp_cache(1,2), mo_integrals_tc_int_map ) -! call get_mo_two_e_integrals_tcdag_int(hfix, p1, p2, mo_num, hijm_cache(1,1), mo_integrals_tcdag_int_map) -! call get_mo_two_e_integrals_tcdag_int(hfix, p2, p1, mo_num, hijm_cache(1,2), mo_integrals_tcdag_int_map) +! call get_mo_bi_ortho_tc_two_es_tc_int (hfix, p1, p2, mo_num, hijp_cache(1,1), mo_integrals_tc_int_map ) +! call get_mo_bi_ortho_tc_two_es_tc_int (hfix, p2, p1, mo_num, hijp_cache(1,2), mo_integrals_tc_int_map ) +! call get_mo_bi_ortho_tc_two_es_tcdag_int(hfix, p1, p2, mo_num, hijm_cache(1,1), mo_integrals_tcdag_int_map) +! call get_mo_bi_ortho_tc_two_es_tcdag_int(hfix, p2, p1, mo_num, hijm_cache(1,2), mo_integrals_tcdag_int_map) ! -! tmp_row_p = 0d0 -! tmp_row_m = 0d0 +! tmp_row_ij_p = 0d0 +! tmp_row_ij_m = 0d0 ! do putj=1, hfix-1 ! if(lbanned(putj, ma)) cycle ! if(banned(putj, puti,bant)) cycle @@ -998,8 +1164,8 @@ end subroutine get_d0 ! hij_m = hij_m * get_phase_bi(phasemask, ma, ma, putj, p1, hfix, p2, N_int) ! !DIR$ LOOP COUNT AVG(4) ! do k=1,N_states -! tmp_row_p(k,putj) = tmp_row_p(k,putj) + hij_p * coefs(k) -! tmp_row_m(k,putj) = tmp_row_m(k,putj) + hij_m * coefs(k) +! tmp_row_ij_p(k,putj) = tmp_row_ij_p(k,putj) + hij_p * coefs(k) +! tmp_row_ij_m(k,putj) = tmp_row_ij_m(k,putj) + hij_m * coefs(k) ! enddo ! endif ! end do @@ -1015,21 +1181,21 @@ end subroutine get_d0 ! hij_m = hij_m * get_phase_bi(phasemask, ma, ma, hfix, p1, putj, p2, N_int) ! !DIR$ LOOP COUNT AVG(4) ! do k=1,N_states -! tmp_row_p(k,putj) = tmp_row_p(k,putj) + hij_p * coefs(k) -! tmp_row_m(k,putj) = tmp_row_m(k,putj) + hij_m * coefs(k) +! tmp_row_ij_p(k,putj) = tmp_row_ij_p(k,putj) + hij_p * coefs(k) +! tmp_row_ij_m(k,putj) = tmp_row_ij_m(k,putj) + hij_m * coefs(k) ! enddo ! endif ! end do ! ! if(ma == 1) then -! mat_p(1:N_states,1:mo_num,puti) = mat_p(1:N_states,1:mo_num,puti) + tmp_row_p(1:N_states,1:mo_num) -! mat_m(1:N_states,1:mo_num,puti) = mat_m(1:N_states,1:mo_num,puti) + tmp_row_m(1:N_states,1:mo_num) +! mat_p(1:N_states,1:mo_num,puti) = mat_p(1:N_states,1:mo_num,puti) + tmp_row_ij_p(1:N_states,1:mo_num) +! mat_m(1:N_states,1:mo_num,puti) = mat_m(1:N_states,1:mo_num,puti) + tmp_row_ij_m(1:N_states,1:mo_num) ! else ! do l=1,mo_num ! !DIR$ LOOP COUNT AVG(4) ! do k=1,N_states -! mat_p(k,puti,l) = mat_p(k,puti,l) + tmp_row_p(k,l) -! mat_m(k,puti,l) = mat_m(k,puti,l) + tmp_row_m(k,l) +! mat_p(k,puti,l) = mat_p(k,puti,l) + tmp_row_ij_p(k,l) +! mat_m(k,puti,l) = mat_m(k,puti,l) + tmp_row_ij_m(k,l) ! enddo ! enddo ! end if @@ -1037,15 +1203,15 @@ end subroutine get_d0 ! ! !MOVE MI ! pfix = p(1,mi) -! tmp_row_p = 0d0 -! tmp_row_m = 0d0 -! tmp_row2_p = 0d0 -! tmp_row2_m = 0d0 +! tmp_row_ij_p = 0d0 +! tmp_row_ij_m = 0d0 +! tmp_row_ij2_p = 0d0 +! tmp_row_ij2_m = 0d0 ! -! call get_mo_two_e_integrals_tc_int (hfix, pfix, p1, mo_num, hijp_cache(1,1), mo_integrals_tc_int_map ) -! call get_mo_two_e_integrals_tc_int (hfix, pfix, p2, mo_num, hijp_cache(1,2), mo_integrals_tc_int_map ) -! call get_mo_two_e_integrals_tcdag_int(hfix, pfix, p1, mo_num, hijm_cache(1,1), mo_integrals_tcdag_int_map) -! call get_mo_two_e_integrals_tcdag_int(hfix, pfix, p2, mo_num, hijm_cache(1,2), mo_integrals_tcdag_int_map) +! call get_mo_bi_ortho_tc_two_es_tc_int (hfix, pfix, p1, mo_num, hijp_cache(1,1), mo_integrals_tc_int_map ) +! call get_mo_bi_ortho_tc_two_es_tc_int (hfix, pfix, p2, mo_num, hijp_cache(1,2), mo_integrals_tc_int_map ) +! call get_mo_bi_ortho_tc_two_es_tcdag_int(hfix, pfix, p1, mo_num, hijm_cache(1,1), mo_integrals_tcdag_int_map) +! call get_mo_bi_ortho_tc_two_es_tcdag_int(hfix, pfix, p2, mo_num, hijm_cache(1,2), mo_integrals_tcdag_int_map) ! ! putj = p1 ! do puti=1,mo_num !HOT @@ -1062,8 +1228,8 @@ end subroutine get_d0 ! hij_m = hij_m * get_phase_bi(phasemask, ma, mi, hfix, p2, puti, pfix, N_int) ! !DIR$ LOOP COUNT AVG(4) ! do k=1,N_states -! tmp_row_p(k,puti) = tmp_row_p(k,puti) + hij_p * coefs(k) -! tmp_row_m(k,puti) = tmp_row_m(k,puti) + hij_m * coefs(k) +! tmp_row_ij_p(k,puti) = tmp_row_ij_p(k,puti) + hij_p * coefs(k) +! tmp_row_ij_m(k,puti) = tmp_row_ij_m(k,puti) + hij_m * coefs(k) ! enddo ! endif ! end if @@ -1078,26 +1244,26 @@ end subroutine get_d0 ! hij_p = hij_p * get_phase_bi(phasemask, ma, mi, hfix, p1, puti, pfix, N_int) ! hij_m = hij_m * get_phase_bi(phasemask, ma, mi, hfix, p1, puti, pfix, N_int) ! do k=1,N_states -! tmp_row2_p(k,puti) = tmp_row2_p(k,puti) + hij_p * coefs(k) -! tmp_row2_m(k,puti) = tmp_row2_m(k,puti) + hij_m * coefs(k) +! tmp_row_ij2_p(k,puti) = tmp_row_ij2_p(k,puti) + hij_p * coefs(k) +! tmp_row_ij2_m(k,puti) = tmp_row_ij2_m(k,puti) + hij_m * coefs(k) ! enddo ! endif ! end if ! end do ! ! if(mi == 1) then -! mat_p(:,:,p1) = mat_p(:,:,p1) + tmp_row_p (:,:) -! mat_p(:,:,p2) = mat_p(:,:,p2) + tmp_row2_p(:,:) -! mat_m(:,:,p1) = mat_m(:,:,p1) + tmp_row_m (:,:) -! mat_m(:,:,p2) = mat_m(:,:,p2) + tmp_row2_m(:,:) +! mat_p(:,:,p1) = mat_p(:,:,p1) + tmp_row_ij_p (:,:) +! mat_p(:,:,p2) = mat_p(:,:,p2) + tmp_row_ij2_p(:,:) +! mat_m(:,:,p1) = mat_m(:,:,p1) + tmp_row_ij_m (:,:) +! mat_m(:,:,p2) = mat_m(:,:,p2) + tmp_row_ij2_m(:,:) ! else ! do l=1,mo_num ! !DIR$ LOOP COUNT AVG(4) ! do k=1,N_states -! mat_p(k,p1,l) = mat_p(k,p1,l) + tmp_row_p (k,l) -! mat_p(k,p2,l) = mat_p(k,p2,l) + tmp_row2_p(k,l) -! mat_m(k,p1,l) = mat_m(k,p1,l) + tmp_row_m (k,l) -! mat_m(k,p2,l) = mat_m(k,p2,l) + tmp_row2_m(k,l) +! mat_p(k,p1,l) = mat_p(k,p1,l) + tmp_row_ij_p (k,l) +! mat_p(k,p2,l) = mat_p(k,p2,l) + tmp_row_ij2_p(k,l) +! mat_m(k,p1,l) = mat_m(k,p1,l) + tmp_row_ij_m (k,l) +! mat_m(k,p2,l) = mat_m(k,p2,l) + tmp_row_ij2_m(k,l) ! enddo ! enddo ! end if @@ -1111,13 +1277,13 @@ end subroutine get_d0 ! p1 = p(turn3(1,i), ma) ! p2 = p(turn3(2,i), ma) ! -! call get_mo_two_e_integrals_tc_int (hfix, p1, p2, mo_num, hijp_cache(1,1), mo_integrals_tc_int_map ) -! call get_mo_two_e_integrals_tc_int (hfix, p2, p1, mo_num, hijp_cache(1,2), mo_integrals_tc_int_map ) -! call get_mo_two_e_integrals_tcdag_int(hfix, p1, p2, mo_num, hijm_cache(1,1), mo_integrals_tcdag_int_map) -! call get_mo_two_e_integrals_tcdag_int(hfix, p2, p1, mo_num, hijm_cache(1,2), mo_integrals_tcdag_int_map) +! call get_mo_bi_ortho_tc_two_es_tc_int (hfix, p1, p2, mo_num, hijp_cache(1,1), mo_integrals_tc_int_map ) +! call get_mo_bi_ortho_tc_two_es_tc_int (hfix, p2, p1, mo_num, hijp_cache(1,2), mo_integrals_tc_int_map ) +! call get_mo_bi_ortho_tc_two_es_tcdag_int(hfix, p1, p2, mo_num, hijm_cache(1,1), mo_integrals_tcdag_int_map) +! call get_mo_bi_ortho_tc_two_es_tcdag_int(hfix, p2, p1, mo_num, hijm_cache(1,2), mo_integrals_tcdag_int_map) ! -! tmp_row_p = 0d0 -! tmp_row_m = 0d0 +! tmp_row_ij_p = 0d0 +! tmp_row_ij_m = 0d0 ! do putj=1,hfix-1 ! if(banned(putj,puti,1)) cycle ! if(lbanned(putj,ma)) cycle @@ -1128,8 +1294,8 @@ end subroutine get_d0 ! if( (hij_p.ne.0.d0) .and. (hij_m.ne.0.d0) ) then ! hij_p = hij_p * get_phase_bi(phasemask, ma, ma, putj, p1, hfix, p2, N_int) ! hij_m = hij_m * get_phase_bi(phasemask, ma, ma, putj, p1, hfix, p2, N_int) -! tmp_row_p(:,putj) = tmp_row_p(:,putj) + hij_p * coefs(:) -! tmp_row_m(:,putj) = tmp_row_m(:,putj) + hij_m * coefs(:) +! tmp_row_ij_p(:,putj) = tmp_row_ij_p(:,putj) + hij_p * coefs(:) +! tmp_row_ij_m(:,putj) = tmp_row_ij_m(:,putj) + hij_m * coefs(:) ! endif ! end do ! do putj=hfix+1,mo_num @@ -1142,18 +1308,18 @@ end subroutine get_d0 ! if( (hij_p.ne.0.d0) .and. (hij_m.ne.0.d0) ) then ! hij_p = hij_p * get_phase_bi(phasemask, ma, ma, hfix, p1, putj, p2, N_int) ! hij_m = hij_m * get_phase_bi(phasemask, ma, ma, hfix, p1, putj, p2, N_int) -! tmp_row_p(:,putj) = tmp_row_p(:,putj) + hij_p * coefs(:) -! tmp_row_m(:,putj) = tmp_row_m(:,putj) + hij_m * coefs(:) +! tmp_row_ij_p(:,putj) = tmp_row_ij_p(:,putj) + hij_p * coefs(:) +! tmp_row_ij_m(:,putj) = tmp_row_ij_m(:,putj) + hij_m * coefs(:) ! endif ! end do ! -! mat_p(:, :puti-1, puti) = mat_p(:, :puti-1, puti) + tmp_row_p(:,:puti-1) -! mat_m(:, :puti-1, puti) = mat_m(:, :puti-1, puti) + tmp_row_m(:,:puti-1) +! mat_p(:, :puti-1, puti) = mat_p(:, :puti-1, puti) + tmp_row_ij_p(:,:puti-1) +! mat_m(:, :puti-1, puti) = mat_m(:, :puti-1, puti) + tmp_row_ij_m(:,:puti-1) ! do l=puti,mo_num ! !DIR$ LOOP COUNT AVG(4) ! do k=1,N_states -! mat_p(k, puti, l) = mat_p(k, puti,l) + tmp_row_p(k,l) -! mat_m(k, puti, l) = mat_m(k, puti,l) + tmp_row_m(k,l) +! mat_p(k, puti, l) = mat_p(k, puti,l) + tmp_row_ij_p(k,l) +! mat_m(k, puti, l) = mat_m(k, puti,l) + tmp_row_ij_m(k,l) ! enddo ! enddo ! end do @@ -1162,15 +1328,15 @@ end subroutine get_d0 ! pfix = p(1,mi) ! p1 = p(1,ma) ! p2 = p(2,ma) -! tmp_row_p = 0d0 -! tmp_row_m = 0d0 -! tmp_row2_p = 0d0 -! tmp_row2_m = 0d0 +! tmp_row_ij_p = 0d0 +! tmp_row_ij_m = 0d0 +! tmp_row_ij2_p = 0d0 +! tmp_row_ij2_m = 0d0 ! -! call get_mo_two_e_integrals_tc_int (hfix, p1, pfix, mo_num, hijp_cache(1,1), mo_integrals_tc_int_map ) -! call get_mo_two_e_integrals_tc_int (hfix, p2, pfix, mo_num, hijp_cache(1,2), mo_integrals_tc_int_map ) -! call get_mo_two_e_integrals_tcdag_int(hfix, p1, pfix, mo_num, hijp_cache(1,1), mo_integrals_tcdag_int_map) -! call get_mo_two_e_integrals_tcdag_int(hfix, p2, pfix, mo_num, hijp_cache(1,2), mo_integrals_tcdag_int_map) +! call get_mo_bi_ortho_tc_two_es_tc_int (hfix, p1, pfix, mo_num, hijp_cache(1,1), mo_integrals_tc_int_map ) +! call get_mo_bi_ortho_tc_two_es_tc_int (hfix, p2, pfix, mo_num, hijp_cache(1,2), mo_integrals_tc_int_map ) +! call get_mo_bi_ortho_tc_two_es_tcdag_int(hfix, p1, pfix, mo_num, hijp_cache(1,1), mo_integrals_tcdag_int_map) +! call get_mo_bi_ortho_tc_two_es_tcdag_int(hfix, p2, pfix, mo_num, hijp_cache(1,2), mo_integrals_tcdag_int_map) ! ! putj = p2 ! do puti=1,mo_num @@ -1186,8 +1352,8 @@ end subroutine get_d0 ! hij_m = hij_m * get_phase_bi(phasemask, mi, ma, hfix, pfix, puti, p1, N_int) ! !DIR$ LOOP COUNT AVG(4) ! do k=1,N_states -! tmp_row_p(k,puti) = tmp_row_p(k,puti) + hij_p * coefs(k) -! tmp_row_m(k,puti) = tmp_row_m(k,puti) + hij_m * coefs(k) +! tmp_row_ij_p(k,puti) = tmp_row_ij_p(k,puti) + hij_p * coefs(k) +! tmp_row_ij_m(k,puti) = tmp_row_ij_m(k,puti) + hij_m * coefs(k) ! enddo ! endif ! end if @@ -1200,28 +1366,28 @@ end subroutine get_d0 ! hij_p = hij_p * get_phase_bi(phasemask, mi, ma, hfix, pfix, puti, p2, N_int) ! hij_m = hij_m * get_phase_bi(phasemask, mi, ma, hfix, pfix, puti, p2, N_int) ! do k=1,N_states -! tmp_row2_p(k,puti) = tmp_row2_p(k,puti) + hij_p * coefs(k) -! tmp_row2_m(k,puti) = tmp_row2_m(k,puti) + hij_m * coefs(k) +! tmp_row_ij2_p(k,puti) = tmp_row_ij2_p(k,puti) + hij_p * coefs(k) +! tmp_row_ij2_m(k,puti) = tmp_row_ij2_m(k,puti) + hij_m * coefs(k) ! enddo ! endif ! end if ! end do -! mat_p(:,:p2-1,p2) = mat_p(:,:p2-1,p2) + tmp_row_p(:,:p2-1) -! mat_m(:,:p2-1,p2) = mat_m(:,:p2-1,p2) + tmp_row_m(:,:p2-1) +! mat_p(:,:p2-1,p2) = mat_p(:,:p2-1,p2) + tmp_row_ij_p(:,:p2-1) +! mat_m(:,:p2-1,p2) = mat_m(:,:p2-1,p2) + tmp_row_ij_m(:,:p2-1) ! do l=p2,mo_num ! !DIR$ LOOP COUNT AVG(4) ! do k=1,N_states -! mat_p(k,p2,l) = mat_p(k,p2,l) + tmp_row_p(k,l) -! mat_m(k,p2,l) = mat_m(k,p2,l) + tmp_row_m(k,l) +! mat_p(k,p2,l) = mat_p(k,p2,l) + tmp_row_ij_p(k,l) +! mat_m(k,p2,l) = mat_m(k,p2,l) + tmp_row_ij_m(k,l) ! enddo ! enddo -! mat_p(:,:p1-1,p1) = mat_p(:,:p1-1,p1) + tmp_row2_p(:,:p1-1) -! mat_m(:,:p1-1,p1) = mat_m(:,:p1-1,p1) + tmp_row2_m(:,:p1-1) +! mat_p(:,:p1-1,p1) = mat_p(:,:p1-1,p1) + tmp_row_ij2_p(:,:p1-1) +! mat_m(:,:p1-1,p1) = mat_m(:,:p1-1,p1) + tmp_row_ij2_m(:,:p1-1) ! do l=p1,mo_num ! !DIR$ LOOP COUNT AVG(4) ! do k=1,N_states -! mat_p(k,p1,l) = mat_p(k,p1,l) + tmp_row2_p(k,l) -! mat_m(k,p1,l) = mat_m(k,p1,l) + tmp_row2_m(k,l) +! mat_p(k,p1,l) = mat_p(k,p1,l) + tmp_row_ij2_p(k,l) +! mat_m(k,p1,l) = mat_m(k,p1,l) + tmp_row_ij2_m(k,l) ! enddo ! enddo ! end if @@ -1280,8 +1446,8 @@ end subroutine get_d0 ! double precision, intent(in) :: coefs(N_states) ! double precision, intent(inout) :: mat_p(N_states, mo_num, mo_num), mat_m(N_states, mo_num, mo_num) ! -! double precision, external :: get_phase_bi, mo_two_e_integral -! double precision, external :: get_mo_two_e_integral_tc_int, get_mo_two_e_integral_tcdag_int +! double precision, external :: get_phase_bi +! double precision, external :: get_mo_bi_ortho_tc_two_e_tc_int, get_mo_bi_ortho_tc_two_e_tcdag_int ! integer, parameter :: bant=1 ! integer :: i, j, k, s, h1, h2, p1, p2, puti, putj ! logical :: ok @@ -1299,8 +1465,8 @@ end subroutine get_d0 ! do p1=1, mo_num ! if(bannedOrb(p1, 1)) cycle ! -! call get_mo_two_e_integrals_tc_int (p1, h2, h1, mo_num, hijp_cache1, mo_integrals_tc_int_map ) -! call get_mo_two_e_integrals_tcdag_int(p1, h2, h1, mo_num, hijm_cache1, mo_integrals_tcdag_int_map) +! call get_mo_bi_ortho_tc_two_es_tc_int (p1, h2, h1, mo_num, hijp_cache1, mo_integrals_tc_int_map ) +! call get_mo_bi_ortho_tc_two_es_tcdag_int(p1, h2, h1, mo_num, hijm_cache1, mo_integrals_tcdag_int_map) ! ! do p2 = 1, mo_num ! if(bannedOrb(p2,2)) cycle @@ -1329,10 +1495,10 @@ end subroutine get_d0 ! do puti=1, mo_num ! if(bannedOrb(puti, sp)) cycle ! -! call get_mo_two_e_integrals_tc_int (puti, p2, p1, mo_num, hijp_cache1, mo_integrals_tc_int_map ) -! call get_mo_two_e_integrals_tc_int (puti, p1, p2, mo_num, hijp_cache2, mo_integrals_tc_int_map ) -! call get_mo_two_e_integrals_tcdag_int(puti, p2, p1, mo_num, hijm_cache1, mo_integrals_tcdag_int_map) -! call get_mo_two_e_integrals_tcdag_int(puti, p1, p2, mo_num, hijm_cache2, mo_integrals_tcdag_int_map) +! call get_mo_bi_ortho_tc_two_es_tc_int (puti, p2, p1, mo_num, hijp_cache1, mo_integrals_tc_int_map ) +! call get_mo_bi_ortho_tc_two_es_tc_int (puti, p1, p2, mo_num, hijp_cache2, mo_integrals_tc_int_map ) +! call get_mo_bi_ortho_tc_two_es_tcdag_int(puti, p2, p1, mo_num, hijm_cache1, mo_integrals_tcdag_int_map) +! call get_mo_bi_ortho_tc_two_es_tcdag_int(puti, p1, p2, mo_num, hijm_cache2, mo_integrals_tcdag_int_map) ! ! do putj=puti+1, mo_num ! if(bannedOrb(putj, sp)) cycle @@ -1344,10 +1510,10 @@ end subroutine get_d0 ! if( (hij_p.eq.0.d0).and.(hij_m.eq.0.d0) ) cycle ! else ! -! hij_p = get_mo_two_e_integral_tc_int (p1, p2, puti, putj, mo_integrals_tc_int_map ) & -! - get_mo_two_e_integral_tc_int (p2, p1, puti, putj, mo_integrals_tc_int_map ) -! hij_m = get_mo_two_e_integral_tcdag_int(p1, p2, puti, putj, mo_integrals_tcdag_int_map) & -! - get_mo_two_e_integral_tcdag_int(p2, p1, puti, putj, mo_integrals_tcdag_int_map) +! hij_p = get_mo_bi_ortho_tc_two_e_tc_int (p1, p2, puti, putj, mo_integrals_tc_int_map ) & +! - get_mo_bi_ortho_tc_two_e_tc_int (p2, p1, puti, putj, mo_integrals_tc_int_map ) +! hij_m = get_mo_bi_ortho_tc_two_e_tcdag_int(p1, p2, puti, putj, mo_integrals_tcdag_int_map) & +! - get_mo_bi_ortho_tc_two_e_tcdag_int(p2, p1, puti, putj, mo_integrals_tcdag_int_map) ! ! if( (hij_p.eq.0.d0).and.(hij_m.eq.0.d0) ) cycle ! @@ -1391,7 +1557,7 @@ subroutine get_d0_reference(gen, phasemask, bannedOrb, banned, mat, mask, h, p, integer :: i, j, s, h1, h2, p1, p2, puti, putj double precision :: hij, phase - double precision, external :: get_phase_bi, mo_two_e_integral + double precision, external :: get_phase_bi logical :: ok integer :: bant @@ -1411,7 +1577,7 @@ subroutine get_d0_reference(gen, phasemask, bannedOrb, banned, mat, mask, h, p, call i_h_j(gen, det, N_int, hij) else phase = get_phase_bi(phasemask, 1, 2, h1, p1, h2, p2, N_int) - hij = mo_two_e_integral(p1, p2, h1, h2) * phase + hij = mo_bi_ortho_tc_two_e(p1, p2, h1, h2) * phase end if mat(:, p1, p2) = mat(:, p1, p2) + coefs(:) * hij end do @@ -1428,7 +1594,7 @@ subroutine get_d0_reference(gen, phasemask, bannedOrb, banned, mat, mask, h, p, call apply_particles(mask, sp,puti,sp,putj, det, ok, N_int) call i_h_j(gen, det, N_int, hij) else - hij = (mo_two_e_integral(p1, p2, puti, putj) - mo_two_e_integral(p2, p1, puti, putj))* get_phase_bi(phasemask, sp, sp, puti, p1 , putj, p2, N_int) + hij = (mo_bi_ortho_tc_two_e(p1, p2, puti, putj) - mo_bi_ortho_tc_two_e(p2, p1, puti, putj))* get_phase_bi(phasemask, sp, sp, puti, p1 , putj, p2, N_int) end if mat(:, puti, putj) = mat(:, puti, putj) + coefs(:) * hij end do @@ -1451,8 +1617,8 @@ subroutine get_d1_reference(gen, phasemask, bannedOrb, banned, mat, mask, h, p, double precision, intent(in) :: coefs(N_states) double precision, intent(inout) :: mat(N_states, mo_num, mo_num) integer, intent(in) :: h(0:2,2), p(0:4,2), sp - double precision :: hij, tmp_row(N_states, mo_num), tmp_row2(N_states, mo_num) - double precision, external :: get_phase_bi, mo_two_e_integral + double precision :: hij, tmp_row_ij(N_states, mo_num), tmp_row_ij2(N_states, mo_num), hji + double precision, external :: get_phase_bi logical :: ok logical, allocatable :: lbanned(:,:) @@ -1489,51 +1655,51 @@ subroutine get_d1_reference(gen, phasemask, bannedOrb, banned, mat, mask, h, p, p1 = p(1,ma) p2 = p(2,ma) if(.not. bannedOrb(puti, mi)) then - tmp_row = 0d0 + tmp_row_ij = 0d0 do putj=1, hfix-1 if(lbanned(putj, ma) .or. banned(putj, puti,bant)) cycle - hij = (mo_two_e_integral(p1, p2, putj, hfix)-mo_two_e_integral(p2,p1,putj,hfix)) * get_phase_bi(phasemask, ma, ma, putj, p1, hfix, p2, N_int) - tmp_row(1:N_states,putj) = tmp_row(1:N_states,putj) + hij * coefs(1:N_states) + hij = (mo_bi_ortho_tc_two_e(p1, p2, putj, hfix)-mo_bi_ortho_tc_two_e(p2,p1,putj,hfix)) * get_phase_bi(phasemask, ma, ma, putj, p1, hfix, p2, N_int) + tmp_row_ij(1:N_states,putj) = tmp_row_ij(1:N_states,putj) + hij * coefs(1:N_states) end do do putj=hfix+1, mo_num if(lbanned(putj, ma) .or. banned(putj, puti,bant)) cycle - hij = (mo_two_e_integral(p1, p2, hfix, putj)-mo_two_e_integral(p2,p1,hfix,putj)) * get_phase_bi(phasemask, ma, ma, hfix, p1, putj, p2, N_int) - tmp_row(1:N_states,putj) = tmp_row(1:N_states,putj) + hij * coefs(1:N_states) + hij = (mo_bi_ortho_tc_two_e(p1, p2, hfix, putj)-mo_bi_ortho_tc_two_e(p2,p1,hfix,putj)) * get_phase_bi(phasemask, ma, ma, hfix, p1, putj, p2, N_int) + tmp_row_ij(1:N_states,putj) = tmp_row_ij(1:N_states,putj) + hij * coefs(1:N_states) end do if(ma == 1) then - mat(1:N_states,1:mo_num,puti) = mat(1:N_states,1:mo_num,puti) + tmp_row(1:N_states,1:mo_num) + mat(1:N_states,1:mo_num,puti) = mat(1:N_states,1:mo_num,puti) + tmp_row_ij(1:N_states,1:mo_num) else - mat(1:N_states,puti,1:mo_num) = mat(1:N_states,puti,1:mo_num) + tmp_row(1:N_states,1:mo_num) + mat(1:N_states,puti,1:mo_num) = mat(1:N_states,puti,1:mo_num) + tmp_row_ij(1:N_states,1:mo_num) end if end if !MOVE MI pfix = p(1,mi) - tmp_row = 0d0 - tmp_row2 = 0d0 + tmp_row_ij = 0d0 + tmp_row_ij2 = 0d0 do puti=1,mo_num if(lbanned(puti,mi)) cycle !p1 fixed putj = p1 if(.not. banned(putj,puti,bant)) then - hij = mo_two_e_integral(p2,pfix,hfix,puti) * get_phase_bi(phasemask, ma, mi, hfix, p2, puti, pfix, N_int) - tmp_row(:,puti) = tmp_row(:,puti) + hij * coefs(:) + hij = mo_bi_ortho_tc_two_e(p2,pfix,hfix,puti) * get_phase_bi(phasemask, ma, mi, hfix, p2, puti, pfix, N_int) + tmp_row_ij(:,puti) = tmp_row_ij(:,puti) + hij * coefs(:) end if putj = p2 if(.not. banned(putj,puti,bant)) then - hij = mo_two_e_integral(p1,pfix,hfix,puti) * get_phase_bi(phasemask, ma, mi, hfix, p1, puti, pfix, N_int) - tmp_row2(:,puti) = tmp_row2(:,puti) + hij * coefs(:) + hij = mo_bi_ortho_tc_two_e(p1,pfix,hfix,puti) * get_phase_bi(phasemask, ma, mi, hfix, p1, puti, pfix, N_int) + tmp_row_ij2(:,puti) = tmp_row_ij2(:,puti) + hij * coefs(:) end if end do if(mi == 1) then - mat(:,:,p1) = mat(:,:,p1) + tmp_row(:,:) - mat(:,:,p2) = mat(:,:,p2) + tmp_row2(:,:) + mat(:,:,p1) = mat(:,:,p1) + tmp_row_ij(:,:) + mat(:,:,p2) = mat(:,:,p2) + tmp_row_ij2(:,:) else - mat(:,p1,:) = mat(:,p1,:) + tmp_row(:,:) - mat(:,p2,:) = mat(:,p2,:) + tmp_row2(:,:) + mat(:,p1,:) = mat(:,p1,:) + tmp_row_ij(:,:) + mat(:,p2,:) = mat(:,p2,:) + tmp_row_ij2(:,:) end if else if(p(0,ma) == 3) then @@ -1542,46 +1708,46 @@ subroutine get_d1_reference(gen, phasemask, bannedOrb, banned, mat, mask, h, p, puti = p(i, ma) p1 = p(turn3(1,i), ma) p2 = p(turn3(2,i), ma) - tmp_row = 0d0 + tmp_row_ij = 0d0 do putj=1,hfix-1 if(lbanned(putj,ma) .or. banned(puti,putj,1)) cycle - hij = (mo_two_e_integral(p1, p2, putj, hfix)-mo_two_e_integral(p2,p1,putj,hfix)) * get_phase_bi(phasemask, ma, ma, putj, p1, hfix, p2, N_int) - tmp_row(:,putj) = tmp_row(:,putj) + hij * coefs(:) + hij = (mo_bi_ortho_tc_two_e(p1, p2, putj, hfix)-mo_bi_ortho_tc_two_e(p2,p1,putj,hfix)) * get_phase_bi(phasemask, ma, ma, putj, p1, hfix, p2, N_int) + tmp_row_ij(:,putj) = tmp_row_ij(:,putj) + hij * coefs(:) end do do putj=hfix+1,mo_num if(lbanned(putj,ma) .or. banned(puti,putj,1)) cycle - hij = (mo_two_e_integral(p1, p2, hfix, putj)-mo_two_e_integral(p2,p1,hfix,putj)) * get_phase_bi(phasemask, ma, ma, hfix, p1, putj, p2, N_int) - tmp_row(:,putj) = tmp_row(:,putj) + hij * coefs(:) + hij = (mo_bi_ortho_tc_two_e(p1, p2, hfix, putj)-mo_bi_ortho_tc_two_e(p2,p1,hfix,putj)) * get_phase_bi(phasemask, ma, ma, hfix, p1, putj, p2, N_int) + tmp_row_ij(:,putj) = tmp_row_ij(:,putj) + hij * coefs(:) end do - mat(:, :puti-1, puti) = mat(:, :puti-1, puti) + tmp_row(:,:puti-1) - mat(:, puti, puti:) = mat(:, puti, puti:) + tmp_row(:,puti:) + mat(:, :puti-1, puti) = mat(:, :puti-1, puti) + tmp_row_ij(:,:puti-1) + mat(:, puti, puti:) = mat(:, puti, puti:) + tmp_row_ij(:,puti:) end do else hfix = h(1,mi) pfix = p(1,mi) p1 = p(1,ma) p2 = p(2,ma) - tmp_row = 0d0 - tmp_row2 = 0d0 + tmp_row_ij = 0d0 + tmp_row_ij2 = 0d0 do puti=1,mo_num if(lbanned(puti,ma)) cycle putj = p2 if(.not. banned(puti,putj,1)) then - hij = mo_two_e_integral(pfix, p1, hfix, puti) * get_phase_bi(phasemask, mi, ma, hfix, pfix, puti, p1, N_int) - tmp_row(:,puti) = tmp_row(:,puti) + hij * coefs(:) + hij = mo_bi_ortho_tc_two_e(pfix, p1, hfix, puti) * get_phase_bi(phasemask, mi, ma, hfix, pfix, puti, p1, N_int) + tmp_row_ij(:,puti) = tmp_row_ij(:,puti) + hij * coefs(:) end if putj = p1 if(.not. banned(puti,putj,1)) then - hij = mo_two_e_integral(pfix, p2, hfix, puti) * get_phase_bi(phasemask, mi, ma, hfix, pfix, puti, p2, N_int) - tmp_row2(:,puti) = tmp_row2(:,puti) + hij * coefs(:) + hij = mo_bi_ortho_tc_two_e(pfix, p2, hfix, puti) * get_phase_bi(phasemask, mi, ma, hfix, pfix, puti, p2, N_int) + tmp_row_ij2(:,puti) = tmp_row_ij2(:,puti) + hij * coefs(:) end if end do - mat(:,:p2-1,p2) = mat(:,:p2-1,p2) + tmp_row(:,:p2-1) - mat(:,p2,p2:) = mat(:,p2,p2:) + tmp_row(:,p2:) - mat(:,:p1-1,p1) = mat(:,:p1-1,p1) + tmp_row2(:,:p1-1) - mat(:,p1,p1:) = mat(:,p1,p1:) + tmp_row2(:,p1:) + mat(:,:p2-1,p2) = mat(:,:p2-1,p2) + tmp_row_ij(:,:p2-1) + mat(:,p2,p2:) = mat(:,p2,p2:) + tmp_row_ij(:,p2:) + mat(:,:p1-1,p1) = mat(:,:p1-1,p1) + tmp_row_ij2(:,:p1-1) + mat(:,p1,p1:) = mat(:,p1,p1:) + tmp_row_ij2(:,p1:) end if end if deallocate(lbanned) @@ -1624,11 +1790,11 @@ subroutine get_d2_reference(gen, phasemask, bannedOrb, banned, mat, mask, h, p, double precision, intent(inout) :: mat(N_states, mo_num, mo_num) integer, intent(in) :: h(0:2,2), p(0:4,2), sp - double precision, external :: get_phase_bi, mo_two_e_integral + double precision, external :: get_phase_bi integer :: i, j, tip, ma, mi, puti, putj - integer :: h1, h2, p1, p2, i1, i2 - double precision :: hij, phase + integer :: h1, h2, p1, p2, i1, i2, mm + double precision :: hij, phase, hji integer, parameter:: turn2d(2,3,4) = reshape((/0,0, 0,0, 0,0, 3,4, 0,0, 0,0, 2,4, 1,4, 0,0, 2,3, 1,3, 1,2 /), (/2,3,4/)) integer, parameter :: turn2(2) = (/2, 1/) @@ -1659,7 +1825,7 @@ subroutine get_d2_reference(gen, phasemask, bannedOrb, banned, mat, mask, h, p, h1 = h(1, ma) h2 = h(2, ma) - hij = (mo_two_e_integral(p1, p2, h1, h2) - mo_two_e_integral(p2,p1, h1, h2)) * get_phase_bi(phasemask, ma, ma, h1, p1, h2, p2, N_int) + hij = (mo_bi_ortho_tc_two_e(p1, p2, h1, h2) - mo_bi_ortho_tc_two_e(p2,p1, h1, h2)) * get_phase_bi(phasemask, ma, ma, h1, p1, h2, p2, N_int) if(ma == 1) then mat(:, putj, puti) = mat(:, putj, puti) + coefs(:) * hij else @@ -1678,7 +1844,7 @@ subroutine get_d2_reference(gen, phasemask, bannedOrb, banned, mat, mask, h, p, if(banned(puti,putj,bant)) cycle p1 = p(turn2(i), 1) - hij = mo_two_e_integral(p1, p2, h1, h2) * get_phase_bi(phasemask, 1, 2, h1, p1, h2, p2,N_int) + hij = mo_bi_ortho_tc_two_e(p1, p2, h1, h2) * get_phase_bi(phasemask, 1, 2, h1, p1, h2, p2,N_int) mat(:, puti, putj) = mat(:, puti, putj) + coefs(:) * hij end do end do @@ -1698,7 +1864,7 @@ subroutine get_d2_reference(gen, phasemask, bannedOrb, banned, mat, mask, h, p, i2 = turn2d(2, i, j) p1 = p(i1, ma) p2 = p(i2, ma) - hij = (mo_two_e_integral(p1, p2, h1, h2) - mo_two_e_integral(p2,p1, h1, h2)) * get_phase_bi(phasemask, ma, ma, h1, p1, h2, p2,N_int) + hij = (mo_bi_ortho_tc_two_e(p1, p2, h1, h2) - mo_bi_ortho_tc_two_e(p2,p1, h1, h2)) * get_phase_bi(phasemask, ma, ma, h1, p1, h2, p2,N_int) mat(:, puti, putj) = mat(:, puti, putj) + coefs(:) * hij end do end do @@ -1712,7 +1878,7 @@ subroutine get_d2_reference(gen, phasemask, bannedOrb, banned, mat, mask, h, p, if(banned(puti,putj,1)) cycle p2 = p(i, ma) - hij = mo_two_e_integral(p1, p2, h1, h2) * get_phase_bi(phasemask, mi, ma, h1, p1, h2, p2,N_int) + hij = mo_bi_ortho_tc_two_e(p1, p2, h1, h2) * get_phase_bi(phasemask, mi, ma, h1, p1, h2, p2,N_int) mat(:, min(puti, putj), max(puti, putj)) = mat(:, min(puti, putj), max(puti, putj)) + coefs(:) * hij end do else ! tip == 4 @@ -1723,7 +1889,7 @@ subroutine get_d2_reference(gen, phasemask, bannedOrb, banned, mat, mask, h, p, p2 = p(2, mi) h1 = h(1, mi) h2 = h(2, mi) - hij = (mo_two_e_integral(p1, p2, h1, h2) - mo_two_e_integral(p2,p1, h1, h2)) * get_phase_bi(phasemask, mi, mi, h1, p1, h2, p2,N_int) + hij = (mo_bi_ortho_tc_two_e(p1, p2, h1, h2) - mo_bi_ortho_tc_two_e(p2,p1, h1, h2)) * get_phase_bi(phasemask, mi, mi, h1, p1, h2, p2,N_int) mat(:, puti, putj) = mat(:, puti, putj) + coefs(:) * hij end if end if diff --git a/src/cipsi_tc_bi_ortho/pt2_stoch_routines.irp.f b/src/cipsi_tc_bi_ortho/pt2_stoch_routines.irp.f index 56e6bd14..e146efb1 100644 --- a/src/cipsi_tc_bi_ortho/pt2_stoch_routines.irp.f +++ b/src/cipsi_tc_bi_ortho/pt2_stoch_routines.irp.f @@ -130,7 +130,7 @@ subroutine ZMQ_pt2(E, pt2_data, pt2_data_err, relative_error, N_in) PROVIDE psi_bilinear_matrix_columns_loc psi_det_alpha_unique psi_det_beta_unique PROVIDE psi_bilinear_matrix_rows psi_det_sorted_tc_order psi_bilinear_matrix_order PROVIDE psi_bilinear_matrix_transp_rows_loc psi_bilinear_matrix_transp_columns - PROVIDE psi_bilinear_matrix_transp_order psi_selectors_coef_transp psi_det_sorted_tc + PROVIDE psi_bilinear_matrix_transp_order psi_selectors_coef_transp_tc psi_det_sorted_tc PROVIDE psi_det_hii selection_weight pseudo_sym PROVIDE n_act_orb n_inact_orb n_core_orb n_virt_orb n_del_orb seniority_max PROVIDE excitation_beta_max excitation_alpha_max excitation_max diff --git a/src/cipsi_tc_bi_ortho/run_selection_slave.irp.f b/src/cipsi_tc_bi_ortho/run_selection_slave.irp.f index e6b016fa..d351cc79 100644 --- a/src/cipsi_tc_bi_ortho/run_selection_slave.irp.f +++ b/src/cipsi_tc_bi_ortho/run_selection_slave.irp.f @@ -23,7 +23,7 @@ subroutine run_selection_slave(thread, iproc, energy) PROVIDE psi_bilinear_matrix_rows psi_det_sorted_tc_order psi_bilinear_matrix_order PROVIDE psi_bilinear_matrix_transp_rows_loc psi_bilinear_matrix_transp_columns PROVIDE psi_bilinear_matrix_transp_order N_int pt2_F pseudo_sym - PROVIDE psi_selectors_coef_transp psi_det_sorted_tc weight_selection + PROVIDE psi_selectors_coef_transp_tc psi_det_sorted_tc weight_selection call pt2_alloc(pt2_data,N_states) diff --git a/src/cipsi_tc_bi_ortho/selection.irp.f b/src/cipsi_tc_bi_ortho/selection.irp.f index 6b93f663..9c695ba8 100644 --- a/src/cipsi_tc_bi_ortho/selection.irp.f +++ b/src/cipsi_tc_bi_ortho/selection.irp.f @@ -81,7 +81,7 @@ subroutine select_singles_and_doubles(i_generator, hole_mask,particle_mask, fock PROVIDE psi_bilinear_matrix_columns_loc psi_det_alpha_unique psi_det_beta_unique PROVIDE psi_bilinear_matrix_rows psi_det_sorted_tc_order psi_bilinear_matrix_order PROVIDE psi_bilinear_matrix_transp_rows_loc psi_bilinear_matrix_transp_columns - PROVIDE psi_bilinear_matrix_transp_order psi_selectors_coef_transp + PROVIDE psi_bilinear_matrix_transp_order psi_selectors_coef_transp_tc PROVIDE psi_selectors_rcoef_bi_orth_transp psi_selectors_lcoef_bi_orth_transp PROVIDE banned_excitation @@ -511,7 +511,7 @@ subroutine splash_pq(mask, sp, det, i_gen, N_sel, bannedOrb, banned, mat, intere integer(bit_kind) :: phasemask(N_int,2) - PROVIDE psi_selectors_coef_transp psi_det_sorted_tc + PROVIDE psi_selectors_coef_transp_tc psi_det_sorted_tc PROVIDE psi_selectors_rcoef_bi_orth_transp psi_selectors_lcoef_bi_orth_transp @@ -564,29 +564,30 @@ subroutine splash_pq(mask, sp, det, i_gen, N_sel, bannedOrb, banned, mat, intere call bitstring_to_list_in_selection(mobMask(1,1), p(1,1), p(0,1), N_int) call bitstring_to_list_in_selection(mobMask(1,2), p(1,2), p(0,2), N_int) - call get_d3_h ( det(1,1,i), bannedOrb, banned, mat , mask, p, sp, psi_selectors_coef_transp (1, interesting(i)) ) - call get_d3_htc( det(1,1,i), bannedOrb, banned, mat_m, mat_p, mask, p, sp, psi_selectors_rcoef_bi_orth_transp(1, interesting(i)) & - , psi_selectors_lcoef_bi_orth_transp(1, interesting(i)) ) + perMask(1,1) = iand(mask(1,1), not(det(1,1,i))) + perMask(1,2) = iand(mask(1,2), not(det(1,2,i))) + do j=2,N_int + perMask(j,1) = iand(mask(j,1), not(det(j,1,i))) + perMask(j,2) = iand(mask(j,2), not(det(j,2,i))) + end do +! call get_d3_h ( det(1,1,i), bannedOrb, banned, mat , mask, p, sp, psi_selectors_coef_transp_tc (1, interesting(i)) ) +! call get_d3_htc( det(1,1,i), bannedOrb, banned, mat_m, mat_p, mask, p, sp, psi_selectors_rcoef_bi_orth_transp(1, interesting(i)) & +! , psi_selectors_lcoef_bi_orth_transp(1, interesting(i)) ) - !perMask(1,1) = iand(mask(1,1), not(det(1,1,i))) - !perMask(1,2) = iand(mask(1,2), not(det(1,2,i))) - !do j=2,N_int - ! perMask(j,1) = iand(mask(j,1), not(det(j,1,i))) - ! perMask(j,2) = iand(mask(j,2), not(det(j,2,i))) - !end do - !call bitstring_to_list_in_selection(perMask(1,1), h(1,1), h(0,1), N_int) - !call bitstring_to_list_in_selection(perMask(1,2), h(1,2), h(0,2), N_int) - !call get_mask_phase(psi_det_sorted_tc(1,1,interesting(i)), phasemask,N_int) - !if(nt == 4) then - ! call get_d2 (det(1,1,i), phasemask, bannedOrb, banned, mat, mask, h, p, sp, psi_selectors_coef_transp(1, interesting(i))) - ! call get_pm2(det(1,1,i), phasemask, bannedOrb, banned, mat_p, mat_m, mask, h, p, sp, psi_selectors_coef_transp(1, interesting(i))) - !elseif(nt == 3) then - ! call get_d1 (det(1,1,i), phasemask, bannedOrb, banned, mat , mask, h, p, sp, psi_selectors_coef_transp(1, interesting(i))) - ! call get_pm1(det(1,1,i), phasemask, bannedOrb, banned, mat_p, mat_m, mask, h, p, sp, psi_selectors_coef_transp(1, interesting(i))) - !else - ! call get_d0 (det(1,1,i), phasemask, bannedOrb, banned, mat , mask, h, p, sp, psi_selectors_coef_transp(1, interesting(i))) - ! call get_pm0(det(1,1,i), phasemask, bannedOrb, banned, mat_p, mat_m, mask, h, p, sp, psi_selectors_coef_transp(1, interesting(i))) - !endif + call bitstring_to_list_in_selection(perMask(1,1), h(1,1), h(0,1), N_int) + call bitstring_to_list_in_selection(perMask(1,2), h(1,2), h(0,2), N_int) + + call get_mask_phase(psi_det_sorted_tc(1,1,interesting(i)), phasemask,N_int) + if(nt == 4) then + call get_d2 (det(1,1,i), phasemask, bannedOrb, banned, mat_p, mat_m, mask, h, p, sp, psi_selectors_coef_transp_tc(1, 1, interesting(i))) +! call get_pm2(det(1,1,i), phasemask, bannedOrb, banned, mat_p, mat_m, mask, h, p, sp, psi_selectors_coef_transp_tc(1, interesting(i))) + elseif(nt == 3) then + call get_d1 (det(1,1,i), phasemask, bannedOrb, banned, mat_p, mat_m, mask, h, p, sp, psi_selectors_coef_transp_tc(1, 1, interesting(i))) +! call get_pm1(det(1,1,i), phasemask, bannedOrb, banned, mat_p, mat_m, mask, h, p, sp, psi_selectors_coef_transp_tc(1, interesting(i))) + else + call get_d0 (det(1,1,i), phasemask, bannedOrb, banned, mat_p, mat_m, mask, h, p, sp, psi_selectors_coef_transp_tc(1, 1, interesting(i))) +! call get_pm0(det(1,1,i), phasemask, bannedOrb, banned, mat_p, mat_m, mask, h, p, sp, psi_selectors_coef_transp_tc(1, interesting(i))) + endif elseif(nt == 4) then call bitstring_to_list_in_selection(mobMask(1,1), p(1,1), p(0,1), N_int) call bitstring_to_list_in_selection(mobMask(1,2), p(1,2), p(0,2), N_int) @@ -775,17 +776,57 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d ! call get_excitation_degree( HF_bitmask, det, degree, N_int) -! psi_h_alpha = mat_m(istate, p1, p2) -! alpha_h_psi = mat_p(istate, p1, p2) + double precision :: alpha_h_psi_tmp, psi_h_alpha_tmp + psi_h_alpha_tmp = mat_m(istate, p1, p2) + alpha_h_psi_tmp = mat_p(istate, p1, p2) ! psi_h_alpha = 0.d0 alpha_h_psi = 0.d0 do iii = 1, N_det call htilde_mu_mat_bi_ortho_tot(psi_det(1,1,iii), det, N_int, i_h_alpha) call htilde_mu_mat_bi_ortho_tot(det, psi_det(1,1,iii), N_int, alpha_h_i) - psi_h_alpha += i_h_alpha * leigvec_tc_bi_orth(iii,1) - alpha_h_psi += alpha_h_i * reigvec_tc_bi_orth(iii,1) +! psi_h_alpha += i_h_alpha * leigvec_tc_bi_orth(iii,1) +! alpha_h_psi += alpha_h_i * reigvec_tc_bi_orth(iii,1) + psi_h_alpha += i_h_alpha * 1.d0 + alpha_h_psi += alpha_h_i * 1.d0 enddo +! print*,'---',p1,p2 +! call debug_det(det,N_int) +! print*,psi_h_alpha *alpha_h_psi, psi_h_alpha, alpha_h_psi +! print*,psi_h_alpha_tmp*alpha_h_psi_tmp,psi_h_alpha_tmp,alpha_h_psi_tmp +! if(dabs(psi_h_alpha - psi_h_alpha_tmp).gt.1.d-10 .or. dabs(alpha_h_psi - alpha_h_psi_tmp).gt.1.d-10)then +! if(dabs(psi_h_alpha_tmp*alpha_h_psi_tmp).gt.1.d+10)then + if(dabs(psi_h_alpha*alpha_h_psi - psi_h_alpha_tmp*alpha_h_psi_tmp).gt.1.d-10)then +! print*,'---' +! print*,psi_h_alpha *alpha_h_psi, psi_h_alpha, alpha_h_psi +! print*,psi_h_alpha_tmp*alpha_h_psi_tmp,psi_h_alpha_tmp,alpha_h_psi_tmp + call debug_det(det,N_int) + print*,dabs(psi_h_alpha*alpha_h_psi - psi_h_alpha_tmp*alpha_h_psi_tmp),psi_h_alpha *alpha_h_psi,psi_h_alpha_tmp*alpha_h_psi_tmp + print*,'-- Good ' + print*, psi_h_alpha, alpha_h_psi + print*,'-- bad ' + print*,psi_h_alpha_tmp,alpha_h_psi_tmp + print*,'-- details good' + do iii = 1, N_det + call get_excitation_degree( psi_det(1,1,iii), det, degree, N_int) + call htilde_mu_mat_bi_ortho_tot(psi_det(1,1,iii), det, N_int, i_h_alpha) + call htilde_mu_mat_bi_ortho_tot(det, psi_det(1,1,iii), N_int, alpha_h_i) + print*,iii,degree,i_h_alpha,alpha_h_i + + enddo +! if(dabs(psi_h_alpha*alpha_h_psi).gt.1.d-10)then +! print*,p1,p2 +! print*,det(1,1), det(1,2) +! call debug_det(det,N_int) +! print*,psi_h_alpha *alpha_h_psi, psi_h_alpha, alpha_h_psi +! print*,psi_h_alpha_tmp*alpha_h_psi_tmp,psi_h_alpha_tmp,alpha_h_psi_tmp +! print*, dabs(psi_h_alpha*alpha_h_psi - psi_h_alpha_tmp*alpha_h_psi_tmp),& +! psi_h_alpha *alpha_h_psi,psi_h_alpha_tmp*alpha_h_psi_tmp + stop + endif +! endif +! stop +! endif !if(alpha_h_psi*psi_h_alpha/delta_E.gt.1.d-10)then ! print*, 'E0,Hii,E_shift' diff --git a/src/cipsi_tc_bi_ortho/stochastic_cipsi.irp.f b/src/cipsi_tc_bi_ortho/stochastic_cipsi.irp.f index 33fe23fc..01cb57dd 100644 --- a/src/cipsi_tc_bi_ortho/stochastic_cipsi.irp.f +++ b/src/cipsi_tc_bi_ortho/stochastic_cipsi.irp.f @@ -89,6 +89,7 @@ subroutine run_stochastic_cipsi call pt2_alloc(pt2_data, N_states) call pt2_alloc(pt2_data_err, N_states) call ZMQ_pt2(E_denom, pt2_data, pt2_data_err, relative_error,to_select) ! Stochastic PT2 and selection + stop N_iter += 1 diff --git a/src/determinants/fock_diag.irp.f b/src/determinants/fock_diag.irp.f index a8ce33b8..c7c951b3 100644 --- a/src/determinants/fock_diag.irp.f +++ b/src/determinants/fock_diag.irp.f @@ -33,59 +33,59 @@ subroutine build_fock_tmp(fock_diag_tmp,det_ref,Nint) ! Occupied MOs do ii=1,elec_alpha_num i = occ(ii,1) - fock_diag_tmp(1,i) = fock_diag_tmp(1,i) + mo_one_e_integrals(i,i) - E0 = E0 + mo_one_e_integrals(i,i) + fock_diag_tmp(1,i) = fock_diag_tmp(1,i) + mo_bi_ortho_tc_one_e(i,i) + E0 = E0 + mo_bi_ortho_tc_one_e(i,i) do jj=1,elec_alpha_num j = occ(jj,1) if (i==j) cycle - fock_diag_tmp(1,i) = fock_diag_tmp(1,i) + mo_two_e_integrals_jj_anti(i,j) - E0 = E0 + 0.5d0*mo_two_e_integrals_jj_anti(i,j) + fock_diag_tmp(1,i) = fock_diag_tmp(1,i) + mo_bi_ortho_tc_two_e_jj_anti(i,j) + E0 = E0 + 0.5d0*mo_bi_ortho_tc_two_e_jj_anti(i,j) enddo do jj=1,elec_beta_num j = occ(jj,2) - fock_diag_tmp(1,i) = fock_diag_tmp(1,i) + mo_two_e_integrals_jj(i,j) - E0 = E0 + mo_two_e_integrals_jj(i,j) + fock_diag_tmp(1,i) = fock_diag_tmp(1,i) + mo_bi_ortho_tc_two_e_jj(i,j) + E0 = E0 + mo_bi_ortho_tc_two_e_jj(i,j) enddo enddo do ii=1,elec_beta_num i = occ(ii,2) - fock_diag_tmp(2,i) = fock_diag_tmp(2,i) + mo_one_e_integrals(i,i) - E0 = E0 + mo_one_e_integrals(i,i) + fock_diag_tmp(2,i) = fock_diag_tmp(2,i) + mo_bi_ortho_tc_one_e(i,i) + E0 = E0 + mo_bi_ortho_tc_one_e(i,i) do jj=1,elec_beta_num j = occ(jj,2) if (i==j) cycle - fock_diag_tmp(2,i) = fock_diag_tmp(2,i) + mo_two_e_integrals_jj_anti(i,j) - E0 = E0 + 0.5d0*mo_two_e_integrals_jj_anti(i,j) + fock_diag_tmp(2,i) = fock_diag_tmp(2,i) + mo_bi_ortho_tc_two_e_jj_anti(i,j) + E0 = E0 + 0.5d0*mo_bi_ortho_tc_two_e_jj_anti(i,j) enddo do jj=1,elec_alpha_num j = occ(jj,1) - fock_diag_tmp(2,i) = fock_diag_tmp(2,i) + mo_two_e_integrals_jj(i,j) + fock_diag_tmp(2,i) = fock_diag_tmp(2,i) + mo_bi_ortho_tc_two_e_jj(i,j) enddo enddo ! Virtual MOs do i=1,mo_num if (fock_diag_tmp(1,i) /= 0.d0) cycle - fock_diag_tmp(1,i) = fock_diag_tmp(1,i) + mo_one_e_integrals(i,i) + fock_diag_tmp(1,i) = fock_diag_tmp(1,i) + mo_bi_ortho_tc_one_e(i,i) do jj=1,elec_alpha_num j = occ(jj,1) - fock_diag_tmp(1,i) = fock_diag_tmp(1,i) + mo_two_e_integrals_jj_anti(i,j) + fock_diag_tmp(1,i) = fock_diag_tmp(1,i) + mo_bi_ortho_tc_two_e_jj_anti(i,j) enddo do jj=1,elec_beta_num j = occ(jj,2) - fock_diag_tmp(1,i) = fock_diag_tmp(1,i) + mo_two_e_integrals_jj(i,j) + fock_diag_tmp(1,i) = fock_diag_tmp(1,i) + mo_bi_ortho_tc_two_e_jj(i,j) enddo enddo do i=1,mo_num if (fock_diag_tmp(2,i) /= 0.d0) cycle - fock_diag_tmp(2,i) = fock_diag_tmp(2,i) + mo_one_e_integrals(i,i) + fock_diag_tmp(2,i) = fock_diag_tmp(2,i) + mo_bi_ortho_tc_one_e(i,i) do jj=1,elec_beta_num j = occ(jj,2) - fock_diag_tmp(2,i) = fock_diag_tmp(2,i) + mo_two_e_integrals_jj_anti(i,j) + fock_diag_tmp(2,i) = fock_diag_tmp(2,i) + mo_bi_ortho_tc_two_e_jj_anti(i,j) enddo do jj=1,elec_alpha_num j = occ(jj,1) - fock_diag_tmp(2,i) = fock_diag_tmp(2,i) + mo_two_e_integrals_jj(i,j) + fock_diag_tmp(2,i) = fock_diag_tmp(2,i) + mo_bi_ortho_tc_two_e_jj(i,j) enddo enddo diff --git a/src/fci_tc_bi/selectors.irp.f b/src/fci_tc_bi/selectors.irp.f index 734c8ed0..502c2b7d 100644 --- a/src/fci_tc_bi/selectors.irp.f +++ b/src/fci_tc_bi/selectors.irp.f @@ -32,6 +32,7 @@ END_PROVIDER BEGIN_PROVIDER [ integer(bit_kind), psi_selectors, (N_int,2,psi_selectors_size) ] &BEGIN_PROVIDER [ double precision, psi_selectors_coef, (psi_selectors_size,N_states) ] +&BEGIN_PROVIDER [ double precision, psi_selectors_coef_tc, (psi_selectors_size,2,N_states) ] implicit none BEGIN_DOC ! Determinants on which we apply for perturbation. @@ -47,12 +48,17 @@ END_PROVIDER do k=1,N_states do i=1,N_det_selectors psi_selectors_coef(i,k) = psi_coef_sorted_tc_gen(i,k) +! psi_selectors_coef_tc(i,1,k) = psi_r_coef_sorted_bi_ortho(i,k) +! psi_selectors_coef_tc(i,2,k) = psi_l_coef_sorted_bi_ortho(i,k) + psi_selectors_coef_tc(i,1,k) = 1.d0 + psi_selectors_coef_tc(i,2,k) = 1.d0 enddo enddo END_PROVIDER -BEGIN_PROVIDER [ double precision, psi_selectors_coef_transp, (N_states,psi_selectors_size) ] + BEGIN_PROVIDER [ double precision, psi_selectors_coef_transp, (N_states,psi_selectors_size) ] +&BEGIN_PROVIDER [ double precision, psi_selectors_coef_transp_tc, (N_states,2,psi_selectors_size) ] implicit none BEGIN_DOC ! Transposed psi_selectors @@ -62,6 +68,8 @@ BEGIN_PROVIDER [ double precision, psi_selectors_coef_transp, (N_states,psi_sele do i=1,N_det_selectors do k=1,N_states psi_selectors_coef_transp(k,i) = psi_selectors_coef(i,k) + psi_selectors_coef_transp_tc(k,1,i) = psi_selectors_coef_tc(i,1,k) + psi_selectors_coef_transp_tc(k,2,i) = psi_selectors_coef_tc(i,2,k) enddo enddo END_PROVIDER diff --git a/src/tc_bi_ortho/slater_tc_opt.irp.f b/src/tc_bi_ortho/slater_tc_opt.irp.f index 8ab3388c..1d7e7a5f 100644 --- a/src/tc_bi_ortho/slater_tc_opt.irp.f +++ b/src/tc_bi_ortho/slater_tc_opt.irp.f @@ -42,3 +42,45 @@ subroutine htilde_mu_mat_opt_bi_ortho(key_j, key_i, Nint, hmono, htwoe, hthree, end ! --- + +subroutine htilde_mu_mat_opt_bi_ortho_no_3e(key_j, key_i, Nint, htot) + + BEGIN_DOC + ! + ! where |key_j> is developed on the LEFT basis and |key_i> is developed on the RIGHT basis + !! + ! Returns the detail of the matrix element WITHOUT ANY CONTRIBUTION FROM THE THREE ELECTRON TERMS + !! WARNING !! + ! + ! Non hermitian !! + ! + END_DOC + + use bitmasks + + implicit none + integer, intent(in) :: Nint + integer(bit_kind), intent(in) :: key_i(Nint,2), key_j(Nint,2) + double precision, intent(out) :: htot + integer :: degree + + htot = 0.d0 + + call get_excitation_degree(key_i, key_j, degree, Nint) + if(degree.gt.2) return + + if(degree == 0)then + call diag_htilde_mu_mat_fock_bi_ortho_no_3e(Nint, key_i,htot) + else if (degree == 1)then + call single_htilde_mu_mat_fock_bi_ortho_no_3e(Nint,key_j, key_i , htot) + else if(degree == 2)then + call double_htilde_mu_mat_fock_bi_ortho_no_3e(Nint, key_j, key_i, htot) + endif + + if(degree==0) then + htot += nuclear_repulsion + endif + +end + +! --- diff --git a/src/tc_bi_ortho/slater_tc_opt_diag.irp.f b/src/tc_bi_ortho/slater_tc_opt_diag.irp.f index c0b59969..68f647dd 100644 --- a/src/tc_bi_ortho/slater_tc_opt_diag.irp.f +++ b/src/tc_bi_ortho/slater_tc_opt_diag.irp.f @@ -277,3 +277,197 @@ subroutine a_tc_operator(iorb,ispin,key,hmono,htwoe,hthree,Nint,na,nb) end + +subroutine diag_htilde_mu_mat_fock_bi_ortho_no_3e(Nint, det_in,htot) + implicit none + BEGIN_DOC + ! Computes $\langle i|H|i \rangle$. WITHOUT ANY CONTRIBUTIONS FROM 3E TERMS + END_DOC + integer,intent(in) :: Nint + integer(bit_kind),intent(in) :: det_in(Nint,2) + double precision, intent(out) :: htot + double precision :: hmono,htwoe + + integer(bit_kind) :: hole(Nint,2) + integer(bit_kind) :: particle(Nint,2) + integer :: i, nexc(2), ispin + integer :: occ_particle(Nint*bit_kind_size,2) + integer :: occ_hole(Nint*bit_kind_size,2) + integer(bit_kind) :: det_tmp(Nint,2) + integer :: na, nb + + ASSERT (Nint > 0) + ASSERT (sum(popcnt(det_in(:,1))) == elec_alpha_num) + ASSERT (sum(popcnt(det_in(:,2))) == elec_beta_num) + + + nexc(1) = 0 + nexc(2) = 0 + do i=1,Nint + hole(i,1) = xor(det_in(i,1),ref_bitmask(i,1)) + hole(i,2) = xor(det_in(i,2),ref_bitmask(i,2)) + particle(i,1) = iand(hole(i,1),det_in(i,1)) + particle(i,2) = iand(hole(i,2),det_in(i,2)) + hole(i,1) = iand(hole(i,1),ref_bitmask(i,1)) + hole(i,2) = iand(hole(i,2),ref_bitmask(i,2)) + nexc(1) = nexc(1) + popcnt(hole(i,1)) + nexc(2) = nexc(2) + popcnt(hole(i,2)) + enddo + + if (nexc(1)+nexc(2) == 0) then + hmono = ref_tc_energy_1e + htwoe = ref_tc_energy_2e + htot = ref_tc_energy_tot + return + endif + + !call debug_det(det_in,Nint) + integer :: tmp(2) + !DIR$ FORCEINLINE + call bitstring_to_list_ab(particle, occ_particle, tmp, Nint) + ASSERT (tmp(1) == nexc(1)) ! Number of particles alpha + ASSERT (tmp(2) == nexc(2)) ! Number of particle beta + !DIR$ FORCEINLINE + call bitstring_to_list_ab(hole, occ_hole, tmp, Nint) + ASSERT (tmp(1) == nexc(1)) ! Number of holes alpha + ASSERT (tmp(2) == nexc(2)) ! Number of holes beta + + + det_tmp = ref_bitmask + hmono = ref_tc_energy_1e + htwoe = ref_tc_energy_2e + do ispin=1,2 + na = elec_num_tab(ispin) + nb = elec_num_tab(iand(ispin,1)+1) + do i=1,nexc(ispin) + !DIR$ FORCEINLINE + call ac_tc_operator_no_3e( occ_particle(i,ispin), ispin, det_tmp, hmono,htwoe, Nint,na,nb) + !DIR$ FORCEINLINE + call a_tc_operator_no_3e ( occ_hole (i,ispin), ispin, det_tmp, hmono,htwoe, Nint,na,nb) + enddo + enddo + htot = hmono+htwoe +end + +subroutine ac_tc_operator_no_3e(iorb,ispin,key,hmono,htwoe,Nint,na,nb) + use bitmasks + implicit none + BEGIN_DOC + ! Routine that computes one- and two-body energy corresponding + ! + ! to the ADDITION of an electron in an orbital 'iorb' of spin 'ispin' + ! + ! onto a determinant 'key'. + ! + ! in output, the determinant key is changed by the ADDITION of that electron + ! + ! and the quantities hmono,htwoe are INCREMENTED + END_DOC + integer, intent(in) :: iorb, ispin, Nint + integer, intent(inout) :: na, nb + integer(bit_kind), intent(inout) :: key(Nint,2) + double precision, intent(inout) :: hmono,htwoe + + integer :: occ(Nint*bit_kind_size,2) + integer :: other_spin + integer :: k,l,i,jj,mm,j,m + double precision :: direct_int, exchange_int + + + if (iorb < 1) then + print *, irp_here, ': iorb < 1' + print *, iorb, mo_num + stop -1 + endif + if (iorb > mo_num) then + print *, irp_here, ': iorb > mo_num' + print *, iorb, mo_num + stop -1 + endif + + ASSERT (ispin > 0) + ASSERT (ispin < 3) + ASSERT (Nint > 0) + + integer :: tmp(2) + !DIR$ FORCEINLINE + call bitstring_to_list_ab(key, occ, tmp, Nint) + ASSERT (tmp(1) == elec_alpha_num) + ASSERT (tmp(2) == elec_beta_num) + + k = shiftr(iorb-1,bit_kind_shift)+1 + ASSERT (k >0) + l = iorb - shiftl(k-1,bit_kind_shift)-1 + ASSERT (l >= 0) + key(k,ispin) = ibset(key(k,ispin),l) + other_spin = iand(ispin,1)+1 + + hmono = hmono + mo_bi_ortho_tc_one_e(iorb,iorb) + + ! Same spin + do i=1,na + htwoe = htwoe + mo_bi_ortho_tc_two_e_jj_anti(occ(i,ispin),iorb) + enddo + + ! Opposite spin + do i=1,nb + htwoe = htwoe + mo_bi_ortho_tc_two_e_jj(occ(i,other_spin),iorb) + enddo + + na = na+1 +end + +subroutine a_tc_operator_no_3e(iorb,ispin,key,hmono,htwoe,Nint,na,nb) + use bitmasks + implicit none + BEGIN_DOC + ! Routine that computes one- and two-body energy corresponding + ! + ! to the REMOVAL of an electron in an orbital 'iorb' of spin 'ispin' + ! + ! onto a determinant 'key'. + ! + ! in output, the determinant key is changed by the REMOVAL of that electron + ! + ! and the quantities hmono,htwoe are INCREMENTED + END_DOC + integer, intent(in) :: iorb, ispin, Nint + integer, intent(inout) :: na, nb + integer(bit_kind), intent(inout) :: key(Nint,2) + double precision, intent(inout) :: hmono,htwoe + + double precision :: direct_int, exchange_int + integer :: occ(Nint*bit_kind_size,2) + integer :: other_spin + integer :: k,l,i,jj,mm,j,m + integer :: tmp(2) + + ASSERT (iorb > 0) + ASSERT (ispin > 0) + ASSERT (ispin < 3) + ASSERT (Nint > 0) + + k = shiftr(iorb-1,bit_kind_shift)+1 + ASSERT (k>0) + l = iorb - shiftl(k-1,bit_kind_shift)-1 + key(k,ispin) = ibclr(key(k,ispin),l) + other_spin = iand(ispin,1)+1 + + !DIR$ FORCEINLINE + call bitstring_to_list_ab(key, occ, tmp, Nint) + na = na-1 + + hmono = hmono - mo_bi_ortho_tc_one_e(iorb,iorb) + + ! Same spin + do i=1,na + htwoe= htwoe- mo_bi_ortho_tc_two_e_jj_anti(occ(i,ispin),iorb) + enddo + + ! Opposite spin + do i=1,nb + htwoe= htwoe- mo_bi_ortho_tc_two_e_jj(occ(i,other_spin),iorb) + enddo + +end + diff --git a/src/tc_bi_ortho/slater_tc_opt_double.irp.f b/src/tc_bi_ortho/slater_tc_opt_double.irp.f index 9d33523b..d094d76e 100644 --- a/src/tc_bi_ortho/slater_tc_opt_double.irp.f +++ b/src/tc_bi_ortho/slater_tc_opt_double.irp.f @@ -419,3 +419,58 @@ subroutine give_contrib_for_bbbb(h1,h2,p1,p2,occ,Ne,contrib) enddo end + +subroutine double_htilde_mu_mat_fock_bi_ortho_no_3e(Nint, key_j, key_i, htot) + + BEGIN_DOC + ! for double excitation ONLY FOR ONE- AND TWO-BODY TERMS + !! + !! WARNING !! + ! + ! Non hermitian !! + END_DOC + + use bitmasks + + implicit none + integer, intent(in) :: Nint + integer(bit_kind), intent(in) :: key_j(Nint,2), key_i(Nint,2) + double precision, intent(out) :: htot + double precision :: hmono, htwoe + integer :: occ(Nint*bit_kind_size,2) + integer :: Ne(2), i, j, ii, jj, ispin, jspin, k, kk + integer :: degree,exc(0:2,2,2) + integer :: h1, p1, h2, p2, s1, s2 + double precision :: get_mo_two_e_integral_tc_int,phase + + + call get_excitation_degree(key_i, key_j, degree, Nint) + + hmono = 0.d0 + htwoe = 0.d0 + htot = 0.d0 + + if(degree.ne.2)then + return + endif + integer :: degree_i,degree_j + call get_excitation_degree(ref_bitmask,key_i,degree_i,N_int) + call get_excitation_degree(ref_bitmask,key_j,degree_j,N_int) + call get_double_excitation(key_i, key_j, exc, phase, Nint) + call decode_exc(exc, 2, h1, p1, h2, p2, s1, s2) + + if(s1.ne.s2)then + ! opposite spin two-body + htwoe = mo_bi_ortho_tc_two_e(p2,p1,h2,h1) + else + ! same spin two-body + ! direct terms + htwoe = mo_bi_ortho_tc_two_e(p2,p1,h2,h1) + ! exchange terms + htwoe -= mo_bi_ortho_tc_two_e(p1,p2,h2,h1) + endif + htwoe *= phase + htot = htwoe + +end + diff --git a/src/tc_bi_ortho/slater_tc_opt_single.irp.f b/src/tc_bi_ortho/slater_tc_opt_single.irp.f index ae41591a..7cff3c73 100644 --- a/src/tc_bi_ortho/slater_tc_opt_single.irp.f +++ b/src/tc_bi_ortho/slater_tc_opt_single.irp.f @@ -458,3 +458,115 @@ BEGIN_PROVIDER [double precision, fock_op_2_e_tc_closed_shell, (mo_num, mo_num) END_PROVIDER + +subroutine single_htilde_mu_mat_fock_bi_ortho_no_3e(Nint, key_j, key_i, htot) + BEGIN_DOC + ! for single excitation ONLY FOR ONE- AND TWO-BODY TERMS + !! + !! WARNING !! + ! + ! Non hermitian !! + END_DOC + + use bitmasks + + implicit none + integer, intent(in) :: Nint + integer(bit_kind), intent(in) :: key_j(Nint,2), key_i(Nint,2) + double precision, intent(out) :: htot + double precision :: hmono, htwoe + integer :: occ(Nint*bit_kind_size,2) + integer :: Ne(2), i, j, ii, jj, ispin, jspin, k, kk + integer :: degree,exc(0:2,2,2) + integer :: h1, p1, h2, p2, s1, s2 + double precision :: get_mo_two_e_integral_tc_int, phase + double precision :: direct_int, exchange_int_12, exchange_int_23, exchange_int_13 + integer :: other_spin(2) + integer(bit_kind) :: key_j_core(Nint,2), key_i_core(Nint,2) + + other_spin(1) = 2 + other_spin(2) = 1 + + hmono = 0.d0 + htwoe = 0.d0 + htot = 0.d0 + call get_excitation_degree(key_i, key_j, degree, Nint) + if(degree.ne.1)then + return + endif + call bitstring_to_list_ab(key_i, occ, Ne, Nint) + + call get_single_excitation(key_i, key_j, exc, phase, Nint) + call decode_exc(exc,1,h1,p1,h2,p2,s1,s2) + call get_single_excitation_from_fock_tc_no_3e(key_i,key_j,h1,p1,s1,phase,hmono,htwoe,htot) +end + + +subroutine get_single_excitation_from_fock_tc_no_3e(key_i,key_j,h,p,spin,phase,hmono,htwoe,htot) + use bitmasks + implicit none + integer,intent(in) :: h,p,spin + double precision, intent(in) :: phase + integer(bit_kind), intent(in) :: key_i(N_int,2), key_j(N_int,2) + double precision, intent(out) :: hmono,htwoe,htot + integer(bit_kind) :: differences(N_int,2) + integer(bit_kind) :: hole(N_int,2) + integer(bit_kind) :: partcl(N_int,2) + integer :: occ_hole(N_int*bit_kind_size,2) + integer :: occ_partcl(N_int*bit_kind_size,2) + integer :: n_occ_ab_hole(2),n_occ_ab_partcl(2) + integer :: i0,i + double precision :: buffer_c(mo_num),buffer_x(mo_num) + do i=1, mo_num + buffer_c(i) = tc_2e_3idx_coulomb_integrals(i,p,h) + buffer_x(i) = tc_2e_3idx_exchange_integrals(i,p,h) + enddo + do i = 1, N_int + differences(i,1) = xor(key_i(i,1),ref_closed_shell_bitmask(i,1)) + differences(i,2) = xor(key_i(i,2),ref_closed_shell_bitmask(i,2)) + hole(i,1) = iand(differences(i,1),ref_closed_shell_bitmask(i,1)) + hole(i,2) = iand(differences(i,2),ref_closed_shell_bitmask(i,2)) + partcl(i,1) = iand(differences(i,1),key_i(i,1)) + partcl(i,2) = iand(differences(i,2),key_i(i,2)) + enddo + call bitstring_to_list_ab(hole, occ_hole, n_occ_ab_hole, N_int) + call bitstring_to_list_ab(partcl, occ_partcl, n_occ_ab_partcl, N_int) + hmono = mo_bi_ortho_tc_one_e(p,h) + htwoe = fock_op_2_e_tc_closed_shell(p,h) + ! holes :: direct terms + do i0 = 1, n_occ_ab_hole(1) + i = occ_hole(i0,1) + htwoe -= buffer_c(i) + enddo + do i0 = 1, n_occ_ab_hole(2) + i = occ_hole(i0,2) + htwoe -= buffer_c(i) + enddo + + ! holes :: exchange terms + do i0 = 1, n_occ_ab_hole(spin) + i = occ_hole(i0,spin) + htwoe += buffer_x(i) + enddo + + ! particles :: direct terms + do i0 = 1, n_occ_ab_partcl(1) + i = occ_partcl(i0,1) + htwoe += buffer_c(i) + enddo + do i0 = 1, n_occ_ab_partcl(2) + i = occ_partcl(i0,2) + htwoe += buffer_c(i) + enddo + + ! particles :: exchange terms + do i0 = 1, n_occ_ab_partcl(spin) + i = occ_partcl(i0,spin) + htwoe -= buffer_x(i) + enddo + htwoe = htwoe * phase + hmono = hmono * phase + htot = htwoe + hmono + +end + From c0f60753c23d73a360409c137ad909ea99017e9a Mon Sep 17 00:00:00 2001 From: eginer Date: Thu, 26 Jan 2023 11:29:36 +0100 Subject: [PATCH 2/2] moved the point charges definition in nuclei and added the interaction between the charges and nuclei --- src/ao_one_e_ints/EZFIO.cfg | 23 -- src/ao_one_e_ints/point_charges.irp.f | 272 ------------------ src/ao_one_e_ints/pot_pt_charges.irp.f | 108 +++++++ src/cas_based_on_top/two_body_dens_rout.irp.f | 2 +- src/cipsi_tc_bi_ortho/get_d.irp.f | 10 +- src/cipsi_tc_bi_ortho/selection.irp.f | 6 + src/cipsi_tc_bi_ortho/stochastic_cipsi.irp.f | 7 +- src/nuclei/EZFIO.cfg | 24 ++ src/nuclei/nuclei.irp.f | 3 + src/nuclei/point_charges.irp.f | 209 ++++++++++++++ .../write_pt_charges.py | 4 +- 11 files changed, 362 insertions(+), 306 deletions(-) delete mode 100644 src/ao_one_e_ints/point_charges.irp.f create mode 100644 src/ao_one_e_ints/pot_pt_charges.irp.f create mode 100644 src/nuclei/point_charges.irp.f rename src/{ao_one_e_ints => nuclei}/write_pt_charges.py (94%) diff --git a/src/ao_one_e_ints/EZFIO.cfg b/src/ao_one_e_ints/EZFIO.cfg index 262301e0..8d4fff57 100644 --- a/src/ao_one_e_ints/EZFIO.cfg +++ b/src/ao_one_e_ints/EZFIO.cfg @@ -106,26 +106,3 @@ interface: ezfio,provider,ocaml default: 1.e-15 ezfio_name: threshold_ao -[n_pts_charge] -type: integer -doc: Number of point charges to be added to the potential -interface: ezfio -default: 0 - -[pts_charge_z] -type: double precision -doc: Charge associated to each point charge -interface: ezfio -size: (ao_one_e_ints.n_pts_charge) - -[pts_charge_coord] -type: double precision -doc: Coordinate of each point charge. -interface: ezfio -size: (ao_one_e_ints.n_pts_charge,3) - -[point_charges] -type: logical -doc: If |true|, point charges (see ao_one_e_ints/write_pt_charges.py) are added to the one-electron potential -interface: ezfio,provider,ocaml -default: False diff --git a/src/ao_one_e_ints/point_charges.irp.f b/src/ao_one_e_ints/point_charges.irp.f deleted file mode 100644 index c038458d..00000000 --- a/src/ao_one_e_ints/point_charges.irp.f +++ /dev/null @@ -1,272 +0,0 @@ - -! --- - - -BEGIN_PROVIDER [ integer, n_pts_charge ] - implicit none - BEGIN_DOC -! Number of point charges to be added to the potential - END_DOC - - logical :: has - PROVIDE ezfio_filename - if (mpi_master) then - - call ezfio_has_ao_one_e_ints_n_pts_charge(has) - if (has) then - write(6,'(A)') '.. >>>>> [ IO READ: n_pts_charge ] <<<<< ..' - call ezfio_get_ao_one_e_ints_n_pts_charge(n_pts_charge) - else - print *, 'ao_one_e_ints/n_pts_charge not found in EZFIO file' - stop 1 - endif - endif - IRP_IF MPI_DEBUG - print *, irp_here, mpi_rank - call MPI_BARRIER(MPI_COMM_WORLD, ierr) - IRP_ENDIF - IRP_IF MPI - include 'mpif.h' - integer :: ierr - call MPI_BCAST( n_pts_charge, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) - if (ierr /= MPI_SUCCESS) then - stop 'Unable to read n_pts_charge with MPI' - endif - IRP_ENDIF - - call write_time(6) - -END_PROVIDER - -BEGIN_PROVIDER [ double precision, pts_charge_z, (n_pts_charge) ] - - BEGIN_DOC - ! Charge associated to each point charge. - END_DOC - - implicit none - logical :: exists - - PROVIDE ezfio_filename - - if (mpi_master) then - call ezfio_has_ao_one_e_ints_pts_charge_z(exists) - endif - - IRP_IF MPI_DEBUG - print *, irp_here, mpi_rank - call MPI_BARRIER(MPI_COMM_WORLD, ierr) - IRP_ENDIF - - IRP_IF MPI - include 'mpif.h' - integer :: ierr - call MPI_BCAST(pts_charge_z, (n_pts_charge), MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) - if (ierr /= MPI_SUCCESS) then - stop 'Unable to read pts_charge_z with MPI' - endif - IRP_ENDIF - - if (exists) then - - if (mpi_master) then - write(6,'(A)') '.. >>>>> [ IO READ: pts_charge_z ] <<<<< ..' - call ezfio_get_ao_one_e_ints_pts_charge_z(pts_charge_z) - IRP_IF MPI - call MPI_BCAST(pts_charge_z, (n_pts_charge), MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) - if (ierr /= MPI_SUCCESS) then - stop 'Unable to read pts_charge_z with MPI' - endif - IRP_ENDIF - endif - - else - - integer :: i - do i = 1, n_pts_charge - pts_charge_z(i) = 0.d0 - enddo - - endif - print*,'Point charges ' - do i = 1, n_pts_charge - print*,'i,pts_charge_z(i)',i,pts_charge_z(i) - enddo - -END_PROVIDER - - -BEGIN_PROVIDER [ double precision, pts_charge_coord, (n_pts_charge,3) ] - - BEGIN_DOC - ! Coordinates of each point charge. - END_DOC - - implicit none - logical :: exists - - PROVIDE ezfio_filename - - if (mpi_master) then - call ezfio_has_ao_one_e_ints_pts_charge_coord(exists) - endif - - IRP_IF MPI_DEBUG - print *, irp_here, mpi_rank - call MPI_BARRIER(MPI_COMM_WORLD, ierr) - IRP_ENDIF - - IRP_IF MPI - include 'mpif.h' - integer :: ierr - call MPI_BCAST(pts_charge_coord, (n_pts_charge), MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) - if (ierr /= MPI_SUCCESS) then - stop 'Unable to read pts_charge_coord with MPI' - endif - IRP_ENDIF - - if (exists) then - - if (mpi_master) then - double precision, allocatable :: buffer(:,:) - allocate (buffer(n_pts_charge,3)) - write(6,'(A)') '.. >>>>> [ IO READ: pts_charge_coord ] <<<<< ..' - call ezfio_get_ao_one_e_ints_pts_charge_coord(buffer) - integer :: i,j - do i=1,3 - do j=1,n_pts_charge - pts_charge_coord(j,i) = buffer(j,i) - enddo - enddo - deallocate(buffer) - IRP_IF MPI - call MPI_BCAST(pts_charge_coord, (n_pts_charge), MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) - if (ierr /= MPI_SUCCESS) then - stop 'Unable to read pts_charge_coord with MPI' - endif - IRP_ENDIF - endif - - else - - do i = 1, n_pts_charge - pts_charge_coord(i,:) = 0.d0 - enddo - - endif - print*,'Coordinates for the point charges ' - do i = 1, n_pts_charge - write(*,'(I3,X,3(F16.8,X))'),i,pts_charge_coord(i,1:3) - enddo - -END_PROVIDER - -! --- - -BEGIN_PROVIDER [ double precision, ao_integrals_pt_chrg, (ao_num,ao_num)] - - BEGIN_DOC - ! Point charge-electron interaction, in the |AO| basis set. - ! - ! :math:`\langle \chi_i | -\sum_A \frac{1}{|r-R_A|} | \chi_j \rangle` - ! - ! These integrals also contain the pseudopotential integrals. - END_DOC - - implicit none - integer :: num_A, num_B, power_A(3), power_B(3) - integer :: i, j, k, l, n_pt_in, m - double precision :: alpha, beta - double precision :: A_center(3),B_center(3),C_center(3) - double precision :: overlap_x,overlap_y,overlap_z,overlap,dx,NAI_pol_mult - - ao_integrals_pt_chrg = 0.d0 - -! if (read_ao_integrals_pt_chrg) then -! -! call ezfio_get_ao_one_e_ints_ao_integrals_pt_chrg(ao_integrals_pt_chrg) -! print *, 'AO N-e integrals read from disk' -! -! else - -! if(use_cosgtos) then -! !print *, " use_cosgtos for ao_integrals_pt_chrg ?", use_cosgtos -! -! do j = 1, ao_num -! do i = 1, ao_num -! ao_integrals_pt_chrg(i,j) = ao_integrals_pt_chrg_cosgtos(i,j) -! enddo -! enddo -! -! else - - !$OMP PARALLEL & - !$OMP DEFAULT (NONE) & - !$OMP PRIVATE (i,j,k,l,m,alpha,beta,A_center,B_center,C_center,power_A,power_B,& - !$OMP num_A,num_B,Z,c,c1,n_pt_in) & - !$OMP SHARED (ao_num,ao_prim_num,ao_expo_ordered_transp,ao_power,ao_nucl,pts_charge_coord,ao_coef_normalized_ordered_transp,nucl_coord,& - !$OMP n_pt_max_integrals,ao_integrals_pt_chrg,n_pts_charge,pts_charge_z) - - n_pt_in = n_pt_max_integrals - - !$OMP DO SCHEDULE (dynamic) - - do j = 1, ao_num - num_A = ao_nucl(j) - power_A(1:3)= ao_power(j,1:3) - A_center(1:3) = nucl_coord(num_A,1:3) - - do i = 1, ao_num - - num_B = ao_nucl(i) - power_B(1:3)= ao_power(i,1:3) - B_center(1:3) = nucl_coord(num_B,1:3) - - do l=1,ao_prim_num(j) - alpha = ao_expo_ordered_transp(l,j) - - do m=1,ao_prim_num(i) - beta = ao_expo_ordered_transp(m,i) - - double precision :: c, c1 - c = 0.d0 - - do k = 1, n_pts_charge - double precision :: Z - Z = pts_charge_z(k) - - C_center(1:3) = pts_charge_coord(k,1:3) - - c1 = NAI_pol_mult( A_center, B_center, power_A, power_B & - , alpha, beta, C_center, n_pt_in ) - - c = c + Z * c1 - - enddo - ao_integrals_pt_chrg(i,j) = ao_integrals_pt_chrg(i,j) & - + ao_coef_normalized_ordered_transp(l,j) & - * ao_coef_normalized_ordered_transp(m,i) * c - enddo - enddo - enddo - enddo - - !$OMP END DO - !$OMP END PARALLEL - -! endif - - -! IF(do_pseudo) THEN -! ao_integrals_pt_chrg += ao_pseudo_integrals -! ENDIF - -! endif - - -! if (write_ao_integrals_pt_chrg) then -! call ezfio_set_ao_one_e_ints_ao_integrals_pt_chrg(ao_integrals_pt_chrg) -! print *, 'AO N-e integrals written to disk' -! endif - -END_PROVIDER diff --git a/src/ao_one_e_ints/pot_pt_charges.irp.f b/src/ao_one_e_ints/pot_pt_charges.irp.f new file mode 100644 index 00000000..93f1acff --- /dev/null +++ b/src/ao_one_e_ints/pot_pt_charges.irp.f @@ -0,0 +1,108 @@ + +BEGIN_PROVIDER [ double precision, ao_integrals_pt_chrg, (ao_num,ao_num)] + + BEGIN_DOC + ! Point charge-electron interaction, in the |AO| basis set. + ! + ! :math:`\langle \chi_i | -\sum_charge charge * \frac{1}{|r-R_charge|} | \chi_j \rangle` + ! + ! Notice the minus sign convention as it is supposed to be for electrons. + END_DOC + + implicit none + integer :: num_A, num_B, power_A(3), power_B(3) + integer :: i, j, k, l, n_pt_in, m + double precision :: alpha, beta + double precision :: A_center(3),B_center(3),C_center(3) + double precision :: overlap_x,overlap_y,overlap_z,overlap,dx,NAI_pol_mult + + ao_integrals_pt_chrg = 0.d0 + +! if (read_ao_integrals_pt_chrg) then +! +! call ezfio_get_ao_one_e_ints_ao_integrals_pt_chrg(ao_integrals_pt_chrg) +! print *, 'AO N-e integrals read from disk' +! +! else + +! if(use_cosgtos) then +! !print *, " use_cosgtos for ao_integrals_pt_chrg ?", use_cosgtos +! +! do j = 1, ao_num +! do i = 1, ao_num +! ao_integrals_pt_chrg(i,j) = ao_integrals_pt_chrg_cosgtos(i,j) +! enddo +! enddo +! +! else + + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (i,j,k,l,m,alpha,beta,A_center,B_center,C_center,power_A,power_B,& + !$OMP num_A,num_B,Z,c,c1,n_pt_in) & + !$OMP SHARED (ao_num,ao_prim_num,ao_expo_ordered_transp,ao_power,ao_nucl,pts_charge_coord,ao_coef_normalized_ordered_transp,nucl_coord,& + !$OMP n_pt_max_integrals,ao_integrals_pt_chrg,n_pts_charge,pts_charge_z) + + n_pt_in = n_pt_max_integrals + + !$OMP DO SCHEDULE (dynamic) + + do j = 1, ao_num + num_A = ao_nucl(j) + power_A(1:3)= ao_power(j,1:3) + A_center(1:3) = nucl_coord(num_A,1:3) + + do i = 1, ao_num + + num_B = ao_nucl(i) + power_B(1:3)= ao_power(i,1:3) + B_center(1:3) = nucl_coord(num_B,1:3) + + do l=1,ao_prim_num(j) + alpha = ao_expo_ordered_transp(l,j) + + do m=1,ao_prim_num(i) + beta = ao_expo_ordered_transp(m,i) + + double precision :: c, c1 + c = 0.d0 + + do k = 1, n_pts_charge + double precision :: Z + Z = pts_charge_z(k) + + C_center(1:3) = pts_charge_coord(k,1:3) + + c1 = NAI_pol_mult( A_center, B_center, power_A, power_B & + , alpha, beta, C_center, n_pt_in ) + + c = c - Z * c1 + + enddo + ao_integrals_pt_chrg(i,j) = ao_integrals_pt_chrg(i,j) & + + ao_coef_normalized_ordered_transp(l,j) & + * ao_coef_normalized_ordered_transp(m,i) * c + enddo + enddo + enddo + enddo + + !$OMP END DO + !$OMP END PARALLEL + +! endif + + +! IF(do_pseudo) THEN +! ao_integrals_pt_chrg += ao_pseudo_integrals +! ENDIF + +! endif + + +! if (write_ao_integrals_pt_chrg) then +! call ezfio_set_ao_one_e_ints_ao_integrals_pt_chrg(ao_integrals_pt_chrg) +! print *, 'AO N-e integrals written to disk' +! endif + +END_PROVIDER diff --git a/src/cas_based_on_top/two_body_dens_rout.irp.f b/src/cas_based_on_top/two_body_dens_rout.irp.f index 4a57a868..5d066831 100644 --- a/src/cas_based_on_top/two_body_dens_rout.irp.f +++ b/src/cas_based_on_top/two_body_dens_rout.irp.f @@ -132,7 +132,7 @@ end subroutine give_n2_cas(r1,r2,istate,n2_psi) implicit none BEGIN_DOC -! returns mu(r), f_psi, n2_psi for a general cas wave function +! returns n2_psi for a general cas wave function END_DOC integer, intent(in) :: istate double precision, intent(in) :: r1(3),r2(3) diff --git a/src/cipsi_tc_bi_ortho/get_d.irp.f b/src/cipsi_tc_bi_ortho/get_d.irp.f index 58b1972a..7fdc5e12 100644 --- a/src/cipsi_tc_bi_ortho/get_d.irp.f +++ b/src/cipsi_tc_bi_ortho/get_d.irp.f @@ -523,10 +523,10 @@ subroutine get_d1(gen, phasemask, bannedOrb, banned, mat_p, mat_m, mask, h, p, s ! call get_mo_bi_ortho_tc_two_es(hfix,pfix,p1,mo_num,hij_cache(1,1),mo_integrals_map) ! call get_mo_bi_ortho_tc_two_es(hfix,pfix,p2,mo_num,hij_cache(1,2),mo_integrals_map) do mm = 1, mo_num - hji_cache(mm,1) = mo_bi_ortho_tc_two_e(pfix,p1,mm,hfix) - hji_cache(mm,2) = mo_bi_ortho_tc_two_e(pfix,p2,mm,hfix) hij_cache(mm,1) = mo_bi_ortho_tc_two_e(mm,hfix,pfix,p1) hij_cache(mm,2) = mo_bi_ortho_tc_two_e(mm,hfix,pfix,p2) + hji_cache(mm,1) = mo_bi_ortho_tc_two_e(pfix,p1,mm,hfix) + hji_cache(mm,2) = mo_bi_ortho_tc_two_e(pfix,p2,mm,hfix) enddo putj = p1 do puti = 1, mo_num !HOT @@ -800,7 +800,7 @@ subroutine get_d0(gen, phasemask, bannedOrb, banned, mat_p, mat_m, mask, h, p, s if(bannedOrb(p1, 1)) cycle ! call get_mo_bi_ortho_tc_two_es(p1,h2,h1,mo_num,hij_cache1,mo_integrals_map) do mm =1, mo_num - hij_cache1(mm) = mo_bi_ortho_tc_two_e(mm,p1,h2,h1) + hji_cache1(mm) = mo_bi_ortho_tc_two_e(mm,p1,h2,h1) hji_cache1(mm) = mo_bi_ortho_tc_two_e(h2,h1,mm,p1) enddo do p2=1, mo_num @@ -811,8 +811,8 @@ subroutine get_d0(gen, phasemask, bannedOrb, banned, mat_p, mat_m, mask, h, p, s call apply_particles(mask, 1,p1,2,p2, det, ok, N_int) ! call i_h_j(gen, det, N_int, hij) !!! GUESS ON THE ORDER - call htilde_mu_mat_opt_bi_ortho_no_3e(det,gen,N_int, hij) - call htilde_mu_mat_opt_bi_ortho_no_3e(gen,det,N_int, hji) + call htilde_mu_mat_opt_bi_ortho_no_3e(det,gen,N_int, hji) + call htilde_mu_mat_opt_bi_ortho_no_3e(gen,det,N_int, hij) else ! print*,'ELSE ' phase = get_phase_bi(phasemask, 1, 2, h1, p1, h2, p2, N_int) diff --git a/src/cipsi_tc_bi_ortho/selection.irp.f b/src/cipsi_tc_bi_ortho/selection.irp.f index 9c695ba8..659e50a8 100644 --- a/src/cipsi_tc_bi_ortho/selection.irp.f +++ b/src/cipsi_tc_bi_ortho/selection.irp.f @@ -807,11 +807,17 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d print*,'-- bad ' print*,psi_h_alpha_tmp,alpha_h_psi_tmp print*,'-- details good' + double precision :: accu_1, accu_2 + accu_1 = 0.d0 + accu_2 = 0.d0 do iii = 1, N_det call get_excitation_degree( psi_det(1,1,iii), det, degree, N_int) call htilde_mu_mat_bi_ortho_tot(psi_det(1,1,iii), det, N_int, i_h_alpha) call htilde_mu_mat_bi_ortho_tot(det, psi_det(1,1,iii), N_int, alpha_h_i) print*,iii,degree,i_h_alpha,alpha_h_i + accu_1 += i_h_alpha + accu_2 += alpha_h_i + print*,accu_1,accu_2 enddo ! if(dabs(psi_h_alpha*alpha_h_psi).gt.1.d-10)then diff --git a/src/cipsi_tc_bi_ortho/stochastic_cipsi.irp.f b/src/cipsi_tc_bi_ortho/stochastic_cipsi.irp.f index 01cb57dd..dc8a4c07 100644 --- a/src/cipsi_tc_bi_ortho/stochastic_cipsi.irp.f +++ b/src/cipsi_tc_bi_ortho/stochastic_cipsi.irp.f @@ -44,9 +44,10 @@ subroutine run_stochastic_cipsi pt2_data % overlap= 0.d0 pt2_data % variance = huge(1.e0) - if (s2_eig) then - call make_s2_eigenfunction - endif + !!!! WARNING !!!! SEEMS TO BE PROBLEM WTH make_s2_eigenfunction !!!! THE DETERMINANTS CAN APPEAR TWICE IN THE WFT DURING SELECTION +! if (s2_eig) then +! call make_s2_eigenfunction +! endif print_pt2 = .False. call diagonalize_CI_tc_bi_ortho(ndet, E_tc,norm,pt2_data,print_pt2) ! call routine_save_right diff --git a/src/nuclei/EZFIO.cfg b/src/nuclei/EZFIO.cfg index 34c27c46..bd25e38a 100644 --- a/src/nuclei/EZFIO.cfg +++ b/src/nuclei/EZFIO.cfg @@ -37,3 +37,27 @@ type: logical doc: If true, the calculation uses periodic boundary conditions interface: ezfio, provider, ocaml default: false + +[n_pts_charge] +type: integer +doc: Number of point charges to be added to the potential +interface: ezfio +default: 0 + +[pts_charge_z] +type: double precision +doc: Charge associated to each point charge +interface: ezfio +size: (nuclei.n_pts_charge) + +[pts_charge_coord] +type: double precision +doc: Coordinate of each point charge. +interface: ezfio +size: (nuclei.n_pts_charge,3) + +[point_charges] +type: logical +doc: If |true|, point charges (see nuclei/write_pt_charges.py) are added to the one-electron potential +interface: ezfio,provider,ocaml +default: False diff --git a/src/nuclei/nuclei.irp.f b/src/nuclei/nuclei.irp.f index c1b5f52f..f765e107 100644 --- a/src/nuclei/nuclei.irp.f +++ b/src/nuclei/nuclei.irp.f @@ -205,6 +205,9 @@ BEGIN_PROVIDER [ double precision, nuclear_repulsion ] enddo enddo nuclear_repulsion *= 0.5d0 + if(point_charges)then + nuclear_repulsion += pt_chrg_nuclei_repulsion + pt_chrg_repulsion + endif end if call write_time(6) diff --git a/src/nuclei/point_charges.irp.f b/src/nuclei/point_charges.irp.f new file mode 100644 index 00000000..86038742 --- /dev/null +++ b/src/nuclei/point_charges.irp.f @@ -0,0 +1,209 @@ +! --- + + +BEGIN_PROVIDER [ integer, n_pts_charge ] + implicit none + BEGIN_DOC +! Number of point charges to be added to the potential + END_DOC + + logical :: has + PROVIDE ezfio_filename + if (mpi_master) then + + call ezfio_has_nuclei_n_pts_charge(has) + if (has) then + write(6,'(A)') '.. >>>>> [ IO READ: n_pts_charge ] <<<<< ..' + call ezfio_get_nuclei_n_pts_charge(n_pts_charge) + else + print *, 'nuclei/n_pts_charge not found in EZFIO file' + stop 1 + endif + endif + IRP_IF MPI_DEBUG + print *, irp_here, mpi_rank + call MPI_BARRIER(MPI_COMM_WORLD, ierr) + IRP_ENDIF + IRP_IF MPI + include 'mpif.h' + integer :: ierr + call MPI_BCAST( n_pts_charge, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) + if (ierr /= MPI_SUCCESS) then + stop 'Unable to read n_pts_charge with MPI' + endif + IRP_ENDIF + + call write_time(6) + +END_PROVIDER + +BEGIN_PROVIDER [ double precision, pts_charge_z, (n_pts_charge) ] + + BEGIN_DOC + ! Charge associated to each point charge. + END_DOC + + implicit none + logical :: exists + + PROVIDE ezfio_filename + + if (mpi_master) then + call ezfio_has_nuclei_pts_charge_z(exists) + endif + + IRP_IF MPI_DEBUG + print *, irp_here, mpi_rank + call MPI_BARRIER(MPI_COMM_WORLD, ierr) + IRP_ENDIF + + IRP_IF MPI + include 'mpif.h' + integer :: ierr + call MPI_BCAST(pts_charge_z, (n_pts_charge), MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) + if (ierr /= MPI_SUCCESS) then + stop 'Unable to read pts_charge_z with MPI' + endif + IRP_ENDIF + + if (exists) then + + if (mpi_master) then + write(6,'(A)') '.. >>>>> [ IO READ: pts_charge_z ] <<<<< ..' + call ezfio_get_nuclei_pts_charge_z(pts_charge_z) + IRP_IF MPI + call MPI_BCAST(pts_charge_z, (n_pts_charge), MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) + if (ierr /= MPI_SUCCESS) then + stop 'Unable to read pts_charge_z with MPI' + endif + IRP_ENDIF + endif + + else + + integer :: i + do i = 1, n_pts_charge + pts_charge_z(i) = 0.d0 + enddo + + endif + print*,'Point charges ' + do i = 1, n_pts_charge + print*,'i,pts_charge_z(i)',i,pts_charge_z(i) + enddo + +END_PROVIDER + + +BEGIN_PROVIDER [ double precision, pts_charge_coord, (n_pts_charge,3) ] + + BEGIN_DOC + ! Coordinates of each point charge. + END_DOC + + implicit none + logical :: exists + + PROVIDE ezfio_filename + + if (mpi_master) then + call ezfio_has_nuclei_pts_charge_coord(exists) + endif + + IRP_IF MPI_DEBUG + print *, irp_here, mpi_rank + call MPI_BARRIER(MPI_COMM_WORLD, ierr) + IRP_ENDIF + + IRP_IF MPI + include 'mpif.h' + integer :: ierr + call MPI_BCAST(pts_charge_coord, (n_pts_charge), MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) + if (ierr /= MPI_SUCCESS) then + stop 'Unable to read pts_charge_coord with MPI' + endif + IRP_ENDIF + + if (exists) then + + if (mpi_master) then + double precision, allocatable :: buffer(:,:) + allocate (buffer(n_pts_charge,3)) + write(6,'(A)') '.. >>>>> [ IO READ: pts_charge_coord ] <<<<< ..' + call ezfio_get_nuclei_pts_charge_coord(buffer) + integer :: i,j + do i=1,3 + do j=1,n_pts_charge + pts_charge_coord(j,i) = buffer(j,i) + enddo + enddo + deallocate(buffer) + IRP_IF MPI + call MPI_BCAST(pts_charge_coord, (n_pts_charge), MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) + if (ierr /= MPI_SUCCESS) then + stop 'Unable to read pts_charge_coord with MPI' + endif + IRP_ENDIF + endif + + else + + do i = 1, n_pts_charge + pts_charge_coord(i,:) = 0.d0 + enddo + + endif + print*,'Coordinates for the point charges ' + do i = 1, n_pts_charge + write(*,'(I3,X,3(F16.8,X))'),i,pts_charge_coord(i,1:3) + enddo + +END_PROVIDER + +! --- +BEGIN_PROVIDER [ double precision, pt_chrg_repulsion] + implicit none + BEGIN_DOC + ! repulsion between the point charges + END_DOC + integer :: i,j + double precision :: Z_A, z_B,A_center(3), B_center(3), dist + pt_chrg_repulsion = 0.d0 + do i = 1, n_pts_charge + Z_A = pts_charge_z(i) + A_center(1:3) = pts_charge_coord(i,1:3) + do j = i+1, n_pts_charge + Z_B = pts_charge_z(j) + B_center(1:3) = pts_charge_coord(j,1:3) + dist = (A_center(1)-B_center(1))**2 + (A_center(2)-B_center(2))**2 + (A_center(3)-B_center(3))**2 + dist = dsqrt(dist) + pt_chrg_repulsion += Z_A*Z_B/dist + enddo + enddo + print*,'Repulsion of point charges ' + print*,'pt_chrg_repulsion = ',pt_chrg_repulsion +END_PROVIDER + +BEGIN_PROVIDER [ double precision, pt_chrg_nuclei_repulsion] + implicit none + BEGIN_DOC + ! repulsion between the point charges and the nuclei + END_DOC + integer :: i,j + double precision :: Z_A, z_B,A_center(3), B_center(3), dist + pt_chrg_nuclei_repulsion = 0.d0 + do i = 1, n_pts_charge + Z_A = pts_charge_z(i) + A_center(1:3) = pts_charge_coord(i,1:3) + do j = 1, nucl_num + Z_B = nucl_charge(j) + B_center(1:3) = nucl_coord(j,1:3) + dist = (A_center(1)-B_center(1))**2 + (A_center(2)-B_center(2))**2 + (A_center(3)-B_center(3))**2 + dist = dsqrt(dist) + pt_chrg_nuclei_repulsion += Z_A*Z_B/dist + enddo + enddo + print*,'Repulsion between point charges and nuclei' + print*,'pt_chrg_nuclei_repulsion = ',pt_chrg_nuclei_repulsion +END_PROVIDER + diff --git a/src/ao_one_e_ints/write_pt_charges.py b/src/nuclei/write_pt_charges.py similarity index 94% rename from src/ao_one_e_ints/write_pt_charges.py rename to src/nuclei/write_pt_charges.py index d4b6d251..d722faa8 100755 --- a/src/ao_one_e_ints/write_pt_charges.py +++ b/src/nuclei/write_pt_charges.py @@ -13,11 +13,11 @@ def zip_in_ezfio(ezfio,tmp): cmdzip="gzip -c "+tmp+" > "+tmpzip os.system(cmdzip) os.system("rm "+tmp) - cmdmv="mv "+tmpzip+" "+EZFIO+"/ao_one_e_ints/"+tmpzip + cmdmv="mv "+tmpzip+" "+EZFIO+"/nuclei/"+tmpzip os.system(cmdmv) def mv_in_ezfio(ezfio,tmp): - cmdmv="mv "+tmp+" "+EZFIO+"/ao_one_e_ints/"+tmp + cmdmv="mv "+tmp+" "+EZFIO+"/nuclei/"+tmp os.system(cmdmv)