PTEROSOR/QuantumPackage
10
0
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2025-01-08 20:33:20 +01:00
Code Releases Activity

Merge branch 'good-dev-tc' of https://github.com/QuantumPackage/qp2 into good-dev-tc

Browse Source
This commit is contained in:
eginer 2022-11-07 22:21:33 +01:00
commit cfcb02745b
9 changed files with 408 additions and 59 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
src/becke_numerical_grid/grid_becke_per_atom.irp.f
View File

@ -46,8 +46,4 @@ END_PROVIDER
enddo enddo
enddo enddo
enddo enddo
END_PROVIDER END_PROVIDER

27
src/dft_utils_in_r/mo_in_r.irp.f
View File

@ -181,3 +181,30 @@
END_PROVIDER END_PROVIDER
BEGIN_PROVIDER [double precision, mo_abs_dist_per_nuclei, (mo_num,nucl_num)]
implicit none
BEGIN_DOC
! mo_abs_dist_per_nuclei(j,i_nucl) = <phi_j| |r-R_nucl| |phi_j>
END_DOC
integer :: ipoint,i_nucl,m,j
double precision :: weight, r(3),r_nucl(3),dist
mo_abs_dist_per_nuclei = 0.d0
do ipoint = 1, n_points_final_grid
r(1) = final_grid_points(1,ipoint)
r(2) = final_grid_points(2,ipoint)
r(3) = final_grid_points(3,ipoint)
weight = final_weight_at_r_vector(ipoint)
do i_nucl = 1, nucl_num
dist = 0.d0
do m = 1, 3
r_nucl(m) = r(m) - nucl_coord_transp(m,i_nucl)
r_nucl(m) *= r_nucl(m)
dist += r_nucl(m)
enddo
dist = dsqrt(dist)
do j = 1, mo_num
mo_abs_dist_per_nuclei(j,i_nucl) += weight * mos_in_r_array(j,ipoint)*mos_in_r_array(j,ipoint) * dist
enddo
enddo
enddo
END_PROVIDER

6
src/tc_bi_ortho/slater_tc_3e.irp.f
View File

@ -169,8 +169,6 @@ subroutine single_htilde_three_body_ints_bi_ort(Nint, key_j, key_i, hthree)
! is == ispin in ::: s1 is is s1 is is s1 is is s1 is is ! is == ispin in ::: s1 is is s1 is is s1 is is s1 is is
! < h1 j i | p1 j i > - < h1 j i | p1 i j > ! < h1 j i | p1 j i > - < h1 j i | p1 i j >
! !
! direct_int = three_body_ints_bi_ort(jj,ii,p1,jj,ii,h1) ! USES THE 6-IDX tensor
! exchange_int = three_body_ints_bi_ort(jj,ii,p1,ii,jj,h1) ! USES THE 6-IDX tensor
direct_int = three_e_4_idx_direct_bi_ort(jj,ii,p1,h1) direct_int = three_e_4_idx_direct_bi_ort(jj,ii,p1,h1)
exchange_int = three_e_4_idx_exch23_bi_ort(jj,ii,p1,h1) exchange_int = three_e_4_idx_exch23_bi_ort(jj,ii,p1,h1)
hthree += direct_int - exchange_int hthree += direct_int - exchange_int
@ -182,14 +180,10 @@ subroutine single_htilde_three_body_ints_bi_ort(Nint, key_j, key_i, hthree)
ii = occ(i,ispin) ! other spin ii = occ(i,ispin) ! other spin
do j = 1, Ne(s1) ! same spin do j = 1, Ne(s1) ! same spin
jj = occ(j,s1) ! same spin jj = occ(j,s1) ! same spin
! direct_int = three_body_ints_bi_ort(jj,ii,p1,jj,ii,h1) ! USES THE 6-IDX tensor
! exchange_int = three_body_ints_bi_ort(jj,ii,p1,h1,ii,jj) ! exchange the spin s1 :: 6-IDX tensor
direct_int = three_e_4_idx_direct_bi_ort(jj,ii,p1,h1) direct_int = three_e_4_idx_direct_bi_ort(jj,ii,p1,h1)
exchange_int = three_e_4_idx_exch13_bi_ort(jj,ii,p1,h1) exchange_int = three_e_4_idx_exch13_bi_ort(jj,ii,p1,h1)
! < h1 j i | p1 j i > - < h1 j i | j p1 i > ! < h1 j i | p1 j i > - < h1 j i | j p1 i >
hthree += direct_int - exchange_int hthree += direct_int - exchange_int
! print*,'h1,p1,ii,jj = ',h1,p1,ii,jj
! print*,direct_int, exchange_int
enddo enddo
enddo enddo
! !

72
src/tc_bi_ortho/test_tc_fock.irp.f
View File

@ -120,44 +120,54 @@ subroutine routine_3()
implicit none implicit none
integer :: i, a, i_ok, s1 integer :: i, a, i_ok, s1
double precision :: hmono, htwoe, hthree, htilde_ij double precision :: hmono, htwoe, hthree, htilde_ij
double precision :: err_ai, err_tot double precision :: err_ai, err_tot, ref, new
integer(bit_kind), allocatable :: det_i(:,:) integer(bit_kind), allocatable :: det_i(:,:)
allocate(det_i(N_int,2)) allocate(det_i(N_int,2))
err_tot = 0.d0 err_tot = 0.d0
s1 = 1 do s1 = 1, 2
det_i = ref_bitmask det_i = ref_bitmask
call debug_det(det_i, N_int) call debug_det(det_i, N_int)
print*, ' HF det' print*, ' HF det'
call debug_det(det_i, N_int) call debug_det(det_i, N_int)
do i = 1, elec_alpha_num ! occupied do i = 1, elec_num_tab(s1)
do a = elec_alpha_num+1, mo_num ! virtual do a = elec_num_tab(s1)+1, mo_num ! virtual
det_i = ref_bitmask det_i = ref_bitmask
call do_single_excitation(det_i, i, a, s1, i_ok) call do_single_excitation(det_i, i, a, s1, i_ok)
if(i_ok == -1) then if(i_ok == -1) then
print*, 'PB !!' print*, 'PB !!'
print*, i, a print*, i, a
stop stop
endif endif
!print*, ' excited det' !print*, ' excited det'
!call debug_det(det_i, N_int) !call debug_det(det_i, N_int)
call htilde_mu_mat_bi_ortho(det_i, ref_bitmask, N_int, hmono, htwoe, hthree, htilde_ij) call htilde_mu_mat_bi_ortho(det_i, ref_bitmask, N_int, hmono, htwoe, hthree, htilde_ij)
err_ai = dabs(htilde_ij) if(dabs(hthree).lt.1.d-10)cycle
if(err_ai .gt. 1d-7) then ref = hthree
print*, ' warning on', i, a if(s1 == 1)then
print*, hmono, htwoe, htilde_ij new = fock_a_tot_3e_bi_orth(a,i)
endif else if(s1 == 2)then
err_tot += err_ai new = fock_b_tot_3e_bi_orth(a,i)
endif
write(22, *) htilde_ij err_ai = dabs(dabs(ref) - dabs(new))
enddo if(err_ai .gt. 1d-7) then
print*,'s1 = ',s1
print*, ' warning on', i, a
print*, ref,new,err_ai
endif
print*, ref,new,err_ai
err_tot += err_ai
write(22, *) htilde_ij
enddo
enddo
enddo enddo
print *, ' err_tot = ', err_tot print *, ' err_tot = ', err_tot

42
src/tc_scf/fock_tc.irp.f
View File

@ -76,13 +76,13 @@ BEGIN_PROVIDER [ double precision, Fock_matrix_tc_ao_beta, (ao_num, ao_num)]
END_PROVIDER END_PROVIDER
! --- ! ---
BEGIN_PROVIDER [ double precision, Fock_matrix_tc_ao_tot, (ao_num, ao_num) ] !BEGIN_PROVIDER [ double precision, Fock_matrix_tc_ao_tot, (ao_num, ao_num) ]
implicit none ! implicit none
BEGIN_DOC ! BEGIN_DOC
! Total alpha+beta TC Fock matrix : h_c + Two-e^TC terms on the AO basis ! ! Total alpha+beta TC Fock matrix : h_c + Two-e^TC terms on the AO basis
END_DOC ! END_DOC
Fock_matrix_tc_ao_tot = 0.5d0 * (Fock_matrix_tc_ao_alpha + Fock_matrix_tc_ao_beta) ! Fock_matrix_tc_ao_tot = 0.5d0 * (Fock_matrix_tc_ao_alpha + Fock_matrix_tc_ao_beta)
END_PROVIDER !END_PROVIDER
! --- ! ---
@ -94,6 +94,9 @@ BEGIN_PROVIDER [ double precision, Fock_matrix_tc_mo_alpha, (mo_num, mo_num) ]
if(bi_ortho)then if(bi_ortho)then
call ao_to_mo_bi_ortho( Fock_matrix_tc_ao_alpha, size(Fock_matrix_tc_ao_alpha, 1) & call ao_to_mo_bi_ortho( Fock_matrix_tc_ao_alpha, size(Fock_matrix_tc_ao_alpha, 1) &
, Fock_matrix_tc_mo_alpha, size(Fock_matrix_tc_mo_alpha, 1) ) , Fock_matrix_tc_mo_alpha, size(Fock_matrix_tc_mo_alpha, 1) )
if(three_body_h_tc)then
Fock_matrix_tc_mo_alpha += fock_a_tot_3e_bi_orth
endif
else else
call ao_to_mo( Fock_matrix_tc_ao_alpha, size(Fock_matrix_tc_ao_alpha, 1) & call ao_to_mo( Fock_matrix_tc_ao_alpha, size(Fock_matrix_tc_ao_alpha, 1) &
, Fock_matrix_tc_mo_alpha, size(Fock_matrix_tc_mo_alpha, 1) ) , Fock_matrix_tc_mo_alpha, size(Fock_matrix_tc_mo_alpha, 1) )
@ -110,6 +113,9 @@ BEGIN_PROVIDER [ double precision, Fock_matrix_tc_mo_beta, (mo_num,mo_num) ]
if(bi_ortho)then if(bi_ortho)then
call ao_to_mo_bi_ortho( Fock_matrix_tc_ao_beta, size(Fock_matrix_tc_ao_beta, 1) & call ao_to_mo_bi_ortho( Fock_matrix_tc_ao_beta, size(Fock_matrix_tc_ao_beta, 1) &
, Fock_matrix_tc_mo_beta, size(Fock_matrix_tc_mo_beta, 1) ) , Fock_matrix_tc_mo_beta, size(Fock_matrix_tc_mo_beta, 1) )
if(three_body_h_tc)then
Fock_matrix_tc_mo_beta += fock_b_tot_3e_bi_orth
endif
else else
call ao_to_mo( Fock_matrix_tc_ao_beta, size(Fock_matrix_tc_ao_beta, 1) & call ao_to_mo( Fock_matrix_tc_ao_beta, size(Fock_matrix_tc_ao_beta, 1) &
, Fock_matrix_tc_mo_beta, size(Fock_matrix_tc_mo_beta, 1) ) , Fock_matrix_tc_mo_beta, size(Fock_matrix_tc_mo_beta, 1) )
@ -117,17 +123,17 @@ BEGIN_PROVIDER [ double precision, Fock_matrix_tc_mo_beta, (mo_num,mo_num) ]
END_PROVIDER END_PROVIDER
BEGIN_PROVIDER [ double precision, Fock_matrix_tc_mo_tot, (mo_num, mo_num)] !BEGIN_PROVIDER [ double precision, Fock_matrix_tc_mo_tot, (mo_num, mo_num)]
implicit none ! implicit none
BEGIN_DOC ! BEGIN_DOC
! Total alpha+beta TC Fock matrix : h_c + Two-e^TC terms on the MO basis ! ! Total alpha+beta TC Fock matrix : h_c + Two-e^TC terms on the MO basis
END_DOC ! END_DOC
Fock_matrix_tc_mo_tot = 0.5d0 * (Fock_matrix_tc_mo_alpha + Fock_matrix_tc_mo_beta) ! Fock_matrix_tc_mo_tot = 0.5d0 * (Fock_matrix_tc_mo_alpha + Fock_matrix_tc_mo_beta)
if(three_body_h_tc) then ! if(three_body_h_tc) then
Fock_matrix_tc_mo_tot += fock_3_mat ! Fock_matrix_tc_mo_tot += fock_3_mat
endif ! endif
!call restore_symmetry(mo_num, mo_num, Fock_matrix_tc_mo_tot, mo_num, 1.d-10) ! !call restore_symmetry(mo_num, mo_num, Fock_matrix_tc_mo_tot, mo_num, 1.d-10)
END_PROVIDER !END_PROVIDER
! --- ! ---

121
src/tc_scf/fock_tc_mo_tot.irp.f Normal file
View File

@ -0,0 +1,121 @@
BEGIN_PROVIDER [ double precision, Fock_matrix_tc_mo_tot, (mo_num,mo_num) ]
&BEGIN_PROVIDER [ double precision, Fock_matrix_tc_diag_mo_tot, (mo_num)]
implicit none
BEGIN_DOC
! Fock matrix on the MO basis.
! For open shells, the ROHF Fock Matrix is ::
!
! | F-K | F + K/2 | F |
! |---------------------------------|
! | F + K/2 | F | F - K/2 |
! |---------------------------------|
! | F | F - K/2 | F + K |
!
!
! F = 1/2 (Fa + Fb)
!
! K = Fb - Fa
!
END_DOC
integer :: i,j,n
if (elec_alpha_num == elec_beta_num) then
Fock_matrix_tc_mo_tot = Fock_matrix_tc_mo_alpha
else
do j=1,elec_beta_num
! F-K
do i=1,elec_beta_num !CC
Fock_matrix_tc_mo_tot(i,j) = 0.5d0*(Fock_matrix_tc_mo_alpha(i,j)+Fock_matrix_tc_mo_beta(i,j))&
- (Fock_matrix_tc_mo_beta(i,j) - Fock_matrix_tc_mo_alpha(i,j))
enddo
! F+K/2
do i=elec_beta_num+1,elec_alpha_num !CA
Fock_matrix_tc_mo_tot(i,j) = 0.5d0*(Fock_matrix_tc_mo_alpha(i,j)+Fock_matrix_tc_mo_beta(i,j))&
+ 0.5d0*(Fock_matrix_tc_mo_beta(i,j) - Fock_matrix_tc_mo_alpha(i,j))
enddo
! F
do i=elec_alpha_num+1, mo_num !CV
Fock_matrix_tc_mo_tot(i,j) = 0.5d0*(Fock_matrix_tc_mo_alpha(i,j)+Fock_matrix_tc_mo_beta(i,j))
enddo
enddo
do j=elec_beta_num+1,elec_alpha_num
! F+K/2
do i=1,elec_beta_num !AC
Fock_matrix_tc_mo_tot(i,j) = 0.5d0*(Fock_matrix_tc_mo_alpha(i,j)+Fock_matrix_tc_mo_beta(i,j))&
+ 0.5d0*(Fock_matrix_tc_mo_beta(i,j) - Fock_matrix_tc_mo_alpha(i,j))
enddo
! F
do i=elec_beta_num+1,elec_alpha_num !AA
Fock_matrix_tc_mo_tot(i,j) = 0.5d0*(Fock_matrix_tc_mo_alpha(i,j)+Fock_matrix_tc_mo_beta(i,j))
enddo
! F-K/2
do i=elec_alpha_num+1, mo_num !AV
Fock_matrix_tc_mo_tot(i,j) = 0.5d0*(Fock_matrix_tc_mo_alpha(i,j)+Fock_matrix_tc_mo_beta(i,j))&
- 0.5d0*(Fock_matrix_tc_mo_beta(i,j) - Fock_matrix_tc_mo_alpha(i,j))
enddo
enddo
do j=elec_alpha_num+1, mo_num
! F
do i=1,elec_beta_num !VC
Fock_matrix_tc_mo_tot(i,j) = 0.5d0*(Fock_matrix_tc_mo_alpha(i,j)+Fock_matrix_tc_mo_beta(i,j))
enddo
! F-K/2
do i=elec_beta_num+1,elec_alpha_num !VA
Fock_matrix_tc_mo_tot(i,j) = 0.5d0*(Fock_matrix_tc_mo_alpha(i,j)+Fock_matrix_tc_mo_beta(i,j))&
- 0.5d0*(Fock_matrix_tc_mo_beta(i,j) - Fock_matrix_tc_mo_alpha(i,j))
enddo
! F+K
do i=elec_alpha_num+1,mo_num !VV
Fock_matrix_tc_mo_tot(i,j) = 0.5d0*(Fock_matrix_tc_mo_alpha(i,j)+Fock_matrix_tc_mo_beta(i,j)) &
+ (Fock_matrix_tc_mo_beta(i,j) - Fock_matrix_tc_mo_alpha(i,j))
enddo
enddo
endif
do i = 1, mo_num
Fock_matrix_tc_diag_mo_tot(i) = Fock_matrix_tc_mo_tot(i,i)
enddo
if(frozen_orb_scf)then
integer :: iorb,jorb
do i = 1, n_core_orb
iorb = list_core(i)
do j = 1, n_act_orb
jorb = list_act(j)
Fock_matrix_tc_mo_tot(iorb,jorb) = 0.d0
Fock_matrix_tc_mo_tot(jorb,iorb) = 0.d0
enddo
enddo
endif
if(no_oa_or_av_opt)then
do i = 1, n_act_orb
iorb = list_act(i)
do j = 1, n_inact_orb
jorb = list_inact(j)
Fock_matrix_tc_mo_tot(iorb,jorb) = 0.d0
Fock_matrix_tc_mo_tot(jorb,iorb) = 0.d0
enddo
do j = 1, n_virt_orb
jorb = list_virt(j)
Fock_matrix_tc_mo_tot(iorb,jorb) = 0.d0
Fock_matrix_tc_mo_tot(jorb,iorb) = 0.d0
enddo
do j = 1, n_core_orb
jorb = list_core(j)
Fock_matrix_tc_mo_tot(iorb,jorb) = 0.d0
Fock_matrix_tc_mo_tot(jorb,iorb) = 0.d0
enddo
enddo
endif
if(.not.bi_ortho .and. three_body_h_tc)then
Fock_matrix_tc_mo_tot += fock_3_mat
endif
END_PROVIDER

1
src/tc_scf/fock_three.irp.f
View File

@ -16,6 +16,7 @@ BEGIN_PROVIDER [ double precision, fock_3_mat, (mo_num, mo_num)]
fock_3_mat(j,i) = -contrib fock_3_mat(j,i) = -contrib
enddo enddo
enddo enddo
else if(bi_ortho.and.three_body_h_tc)then
!! !$OMP END DO !! !$OMP END DO
!! !$OMP END PARALLEL !! !$OMP END PARALLEL
!! do i = 1, mo_num !! do i = 1, mo_num

177
src/tc_scf/fock_three_bi_ortho_new.irp.f Normal file
View File

@ -0,0 +1,177 @@
BEGIN_PROVIDER [ double precision, fock_a_abb_3e_bi_orth, (mo_num, mo_num)]
implicit none
BEGIN_DOC
! fock_a_abb_3e_bi_orth(a,i) = bi-ortho 3-e Fock matrix for alpha electrons from alpha,beta,beta contribution
END_DOC
fock_a_abb_3e_bi_orth = 0.d0
integer :: i,a,j,k
double precision :: direct_int, exch_23_int
do i = 1, mo_num
do a = 1, mo_num
do j = 1, elec_beta_num
do k = j+1, elec_beta_num
call give_integrals_3_body_bi_ort(a, k, j, i, k, j, direct_int) ! < a k j | i k j >
call give_integrals_3_body_bi_ort(a, k, j, i, j, k, exch_23_int)! < a k j | i j k > : E_23
fock_a_abb_3e_bi_orth(a,i) += direct_int - exch_23_int
enddo
enddo
enddo
enddo
fock_a_abb_3e_bi_orth = - fock_a_abb_3e_bi_orth
END_PROVIDER
BEGIN_PROVIDER [ double precision, fock_a_aba_3e_bi_orth, (mo_num, mo_num)]
implicit none
BEGIN_DOC
! fock_a_aba_3e_bi_orth(a,i) = bi-ortho 3-e Fock matrix for alpha electrons from alpha,alpha,beta contribution
END_DOC
fock_a_aba_3e_bi_orth = 0.d0
integer :: i,a,j,k
double precision :: direct_int, exch_13_int
do i = 1, mo_num
do a = 1, mo_num
do j = 1, elec_alpha_num ! a
do k = 1, elec_beta_num ! b
! a b a a b a
call give_integrals_3_body_bi_ort(a, k, j, i, k, j, direct_int )! < a k j | i k j >
call give_integrals_3_body_bi_ort(a, k, j, j, k, i, exch_13_int)! < a k j | j k i > : E_13
fock_a_aba_3e_bi_orth(a,i) += direct_int - exch_13_int
enddo
enddo
enddo
enddo
fock_a_aba_3e_bi_orth = - fock_a_aba_3e_bi_orth
END_PROVIDER
BEGIN_PROVIDER [ double precision, fock_a_aaa_3e_bi_orth, (mo_num, mo_num)]
implicit none
BEGIN_DOC
! fock_a_aaa_3e_bi_orth(a,i) = bi-ortho 3-e Fock matrix for alpha electrons from alpha,alpha,alpha contribution
END_DOC
fock_a_aaa_3e_bi_orth = 0.d0
integer :: i,a,j,k
double precision :: direct_int, exch_13_int, exch_23_int, exch_12_int, c_3_int, c_minus_3_int
do i = 1, mo_num
do a = 1, mo_num
do j = 1, elec_alpha_num
do k = j+1, elec_alpha_num
! positive terms :: cycle contrib
call give_integrals_3_body_bi_ort(a, k, j, i, k, j, direct_int )!!! < a k j | i k j >
call give_integrals_3_body_bi_ort(a, k, j, j, i, k, c_3_int) ! < a k j | j i k >
call give_integrals_3_body_bi_ort(a, k, j, k, j, i, c_minus_3_int)! < a k j | k j i >
fock_a_aaa_3e_bi_orth(a,i) += direct_int + c_3_int + c_minus_3_int
! negative terms :: exchange contrib
call give_integrals_3_body_bi_ort(a, k, j, j, k, i, exch_13_int)!!! < a k j | j k i > : E_13
call give_integrals_3_body_bi_ort(a, k, j, i, j, k, exch_23_int)!!! < a k j | i j k > : E_23
call give_integrals_3_body_bi_ort(a, k, j, k, i, j, exch_12_int)!!! < a k j | k i j > : E_12
fock_a_aaa_3e_bi_orth(a,i) += - exch_13_int - exch_23_int - exch_12_int
enddo
enddo
enddo
enddo
fock_a_aaa_3e_bi_orth = - fock_a_aaa_3e_bi_orth
END_PROVIDER
BEGIN_PROVIDER [double precision, fock_a_tot_3e_bi_orth, (mo_num, mo_num)]
implicit none
BEGIN_DOC
! fock_a_tot_3e_bi_orth = bi-ortho 3-e Fock matrix for alpha electrons from all possible spin contributions
END_DOC
fock_a_tot_3e_bi_orth = fock_a_abb_3e_bi_orth + fock_a_aba_3e_bi_orth + fock_a_aaa_3e_bi_orth
END_PROVIDER
BEGIN_PROVIDER [ double precision, fock_b_baa_3e_bi_orth, (mo_num, mo_num)]
implicit none
BEGIN_DOC
! fock_b_baa_3e_bi_orth(a,i) = bi-ortho 3-e Fock matrix for beta electrons from beta,alpha,alpha contribution
END_DOC
fock_b_baa_3e_bi_orth = 0.d0
integer :: i,a,j,k
double precision :: direct_int, exch_23_int
do i = 1, mo_num
do a = 1, mo_num
do j = 1, elec_alpha_num
do k = j+1, elec_alpha_num
call give_integrals_3_body_bi_ort(a, k, j, i, k, j, direct_int) ! < a k j | i k j >
call give_integrals_3_body_bi_ort(a, k, j, i, j, k, exch_23_int)! < a k j | i j k > : E_23
fock_b_baa_3e_bi_orth(a,i) += direct_int - exch_23_int
enddo
enddo
enddo
enddo
fock_b_baa_3e_bi_orth = - fock_b_baa_3e_bi_orth
END_PROVIDER
BEGIN_PROVIDER [ double precision, fock_b_bab_3e_bi_orth, (mo_num, mo_num)]
implicit none
BEGIN_DOC
! fock_b_bab_3e_bi_orth(a,i) = bi-ortho 3-e Fock matrix for beta electrons from beta,alpha,beta contribution
END_DOC
fock_b_bab_3e_bi_orth = 0.d0
integer :: i,a,j,k
double precision :: direct_int, exch_13_int
do i = 1, mo_num
do a = 1, mo_num
do j = 1, elec_beta_num
do k = 1, elec_alpha_num
! b a b b a b
call give_integrals_3_body_bi_ort(a, k, j, i, k, j, direct_int) ! < a k j | i k j >
call give_integrals_3_body_bi_ort(a, k, j, j, k, i, exch_13_int)! < a k j | j k i > : E_13
fock_b_bab_3e_bi_orth(a,i) += direct_int - exch_13_int
enddo
enddo
enddo
enddo
fock_b_bab_3e_bi_orth = - fock_b_bab_3e_bi_orth
END_PROVIDER
BEGIN_PROVIDER [ double precision, fock_b_bbb_3e_bi_orth, (mo_num, mo_num)]
implicit none
BEGIN_DOC
! fock_b_bbb_3e_bi_orth(a,i) = bi-ortho 3-e Fock matrix for alpha electrons from alpha,alpha,alpha contribution
END_DOC
fock_b_bbb_3e_bi_orth = 0.d0
integer :: i,a,j,k
double precision :: direct_int, exch_13_int, exch_23_int, exch_12_int, c_3_int, c_minus_3_int
do i = 1, mo_num
do a = 1, mo_num
do j = 1, elec_beta_num
do k = j+1, elec_beta_num
! positive terms :: cycle contrib
call give_integrals_3_body_bi_ort(a, k, j, i, k, j, direct_int )!!! < a k j | i k j >
call give_integrals_3_body_bi_ort(a, k, j, j, i, k, c_3_int) ! < a k j | j i k >
call give_integrals_3_body_bi_ort(a, k, j, k, j, i, c_minus_3_int)! < a k j | k j i >
fock_b_bbb_3e_bi_orth(a,i) += direct_int + c_3_int + c_minus_3_int
! negative terms :: exchange contrib
call give_integrals_3_body_bi_ort(a, k, j, j, k, i, exch_13_int)!!! < a k j | j k i > : E_13
call give_integrals_3_body_bi_ort(a, k, j, i, j, k, exch_23_int)!!! < a k j | i j k > : E_23
call give_integrals_3_body_bi_ort(a, k, j, k, i, j, exch_12_int)!!! < a k j | k i j > : E_12
fock_b_bbb_3e_bi_orth(a,i) += - exch_13_int - exch_23_int - exch_12_int
enddo
enddo
enddo
enddo
fock_b_bbb_3e_bi_orth = - fock_b_bbb_3e_bi_orth
END_PROVIDER
BEGIN_PROVIDER [ double precision, fock_b_tot_3e_bi_orth, (mo_num, mo_num)]
implicit none
BEGIN_DOC
! fock_b_tot_3e_bi_orth = bi-ortho 3-e Fock matrix for alpha electrons from all possible spin contributions
END_DOC
fock_b_tot_3e_bi_orth = fock_b_bbb_3e_bi_orth + fock_b_bab_3e_bi_orth + fock_b_baa_3e_bi_orth
END_PROVIDER

17
src/tools/print_wf.irp.f
View File

@ -17,6 +17,7 @@ program print_wf
! psi_coef_sorted are the wave function stored in the |EZFIO| directory. ! psi_coef_sorted are the wave function stored in the |EZFIO| directory.
read_wf = .True. read_wf = .True.
touch read_wf touch read_wf
call write_wf
call routine call routine
end end
@ -120,3 +121,19 @@ subroutine routine
print*,'L2 norm of pert beta = ',norm_mono_b_pert_2 print*,'L2 norm of pert beta = ',norm_mono_b_pert_2
end end
subroutine write_wf
implicit none
character*(128) :: output
integer :: i_unit_output,getUnitAndOpen
output=trim(ezfio_filename)//'.wf'
i_unit_output = getUnitAndOpen(output,'w')
integer :: i
print*,'Writing the sorted wf'
do i = 1, N_det
write(i_unit_output,*)i,psi_coef_sorted(i,1)/psi_coef_sorted(1,1)
enddo
end