LCPQ/qp2
9
1
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-12-21 11:03:29 +01:00
Code Releases Activity

Merge branch 'dev-stable-tc-scf' of https://github.com/AbdAmmar/qp2 into dev-stable-tc-scf

Browse Source
This commit is contained in:
Abdallah Ammar 2023-06-02 20:55:54 +02:00
commit 381c98c3ca
10 changed files with 247 additions and 211 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

15
scripts/qp_import_trexio.py
View File

@ -17,6 +17,7 @@ import numpy as np
from functools import reduce from functools import reduce
from ezfio import ezfio from ezfio import ezfio
from docopt import docopt from docopt import docopt
import qp_bitmasks
try: try:
import trexio import trexio
@ -453,6 +454,20 @@ def write_ezfio(trexio_filename, filename):
else: else:
print("None") print("None")
print("Determinant\t\t...\t", end=' ')
alpha = [ i for i in range(num_alpha) ]
beta = [ i for i in range(num_beta) ]
if trexio.has_mo_spin(trexio_file):
spin = trexio.read_mo_spin(trexio_file)
beta = [ i for i in range(mo_num) if spin[i] == 1 ]
beta = [ beta[i] for i in range(num_beta) ]
alpha = qp_bitmasks.BitMask(alpha)
beta = qp_bitmasks.BitMask(beta )
print(alpha)
print(beta)
print("OK")

6
scripts/utility/qp_bitmasks.py
View File

@ -22,7 +22,7 @@ def int_to_string(s):
assert s>=0 assert s>=0
AssertionError AssertionError
""" """
assert type(s) in (int, long) assert type(s) == int
assert s>=0 assert s>=0
return '{s:0b}'.format(s=s) return '{s:0b}'.format(s=s)
@ -62,7 +62,7 @@ def int_to_bitmask(s,bit_kind_size=BIT_KIND_SIZE):
['1111111111111111111111111111111111111111111111111111111111110110'] ['1111111111111111111111111111111111111111111111111111111111110110']
>>> >>>
""" """
assert type(s) in (int, long) assert type(s) == int
if s < 0: if s < 0:
s = s + (1 << bit_kind_size) s = s + (1 << bit_kind_size)
return ['{s:0{width}b}'.format(s=s,width=bit_kind_size)] return ['{s:0{width}b}'.format(s=s,width=bit_kind_size)]
@ -104,7 +104,7 @@ class BitMask(object):
return self._data_int[i] return self._data_int[i]
def __setitem__(self,i,value): def __setitem__(self,i,value):
if type(value) in (int,long): if type(value) == int :
self._data_int[i] = value self._data_int[i] = value
elif type(value) == str: elif type(value) == str:
s = string_to_bitmask(value,bit_kind_size=self.bit_kind_size)[0] s = string_to_bitmask(value,bit_kind_size=self.bit_kind_size)[0]

3
src/ao_basis/cosgtos.irp.f
View File

@ -6,13 +6,14 @@ BEGIN_PROVIDER [ logical, use_cosgtos ]
logical :: has logical :: has
PROVIDE ezfio_filename PROVIDE ezfio_filename
use_cosgtos = .False.
if (mpi_master) then if (mpi_master) then
call ezfio_has_ao_basis_use_cosgtos(has) call ezfio_has_ao_basis_use_cosgtos(has)
if (has) then if (has) then
! write(6,'(A)') '.. >>>>> [ IO READ: use_cosgtos ] <<<<< ..' ! write(6,'(A)') '.. >>>>> [ IO READ: use_cosgtos ] <<<<< ..'
call ezfio_get_ao_basis_use_cosgtos(use_cosgtos) call ezfio_get_ao_basis_use_cosgtos(use_cosgtos)
else else
use_cosgtos = .False. call ezfio_set_ao_basis_use_cosgtos(use_cosgtos)
endif endif
endif endif
IRP_IF MPI_DEBUG IRP_IF MPI_DEBUG

115
src/bi_ort_ints/bi_ort_ints.irp.f
View File

@ -55,6 +55,7 @@ subroutine test_5idx
implicit none implicit none
integer :: i,k,j,l,m,n,ipoint integer :: i,k,j,l,m,n,ipoint
double precision :: accu, contrib,new,ref double precision :: accu, contrib,new,ref
double precision, external :: three_e_5_idx_exch12_bi_ort
i = 1 i = 1
k = 1 k = 1
n = 0 n = 0
@ -64,18 +65,21 @@ subroutine test_5idx
do j = 1, mo_num do j = 1, mo_num
do l = 1, mo_num do l = 1, mo_num
do m = 1, mo_num do m = 1, mo_num
! if (dabs(three_e_5_idx_direct_bi_ort(m,l,j,k,i) - three_e_5_idx_exch12_bi_ort(m,l,i,k,j)) > 1.d-10) then
! stop
! endif
new = three_e_5_idx_direct_bi_ort(m,l,j,k,i) ! new = three_e_5_idx_direct_bi_ort(m,l,j,k,i)
ref = three_e_5_idx_direct_bi_ort_old(m,l,j,k,i) ! ref = three_e_5_idx_direct_bi_ort_old(m,l,j,k,i)
contrib = dabs(new - ref) ! contrib = dabs(new - ref)
accu += contrib ! accu += contrib
if(contrib .gt. 1.d-10)then ! if(contrib .gt. 1.d-10)then
print*,'direct' ! print*,'direct'
print*,i,k,j,l,m ! print*,i,k,j,l,m
print*,ref,new,contrib ! print*,ref,new,contrib
stop ! stop
endif ! endif
!
new = three_e_5_idx_exch12_bi_ort(m,l,j,k,i) new = three_e_5_idx_exch12_bi_ort(m,l,j,k,i)
ref = three_e_5_idx_exch12_bi_ort_old(m,l,j,k,i) ref = three_e_5_idx_exch12_bi_ort_old(m,l,j,k,i)
contrib = dabs(new - ref) contrib = dabs(new - ref)
@ -86,51 +90,52 @@ subroutine test_5idx
print*,ref,new,contrib print*,ref,new,contrib
stop stop
endif endif
!
! new = three_e_5_idx_cycle_1_bi_ort(m,l,j,k,i)
! ref = three_e_5_idx_cycle_1_bi_ort_old(m,l,j,k,i)
! contrib = dabs(new - ref)
! accu += contrib
! if(contrib .gt. 1.d-10)then
! print*,'cycle1'
! print*,i,k,j,l,m
! print*,ref,new,contrib
! stop
! endif
!
! new = three_e_5_idx_cycle_2_bi_ort(m,l,j,k,i)
! ref = three_e_5_idx_cycle_2_bi_ort_old(m,l,j,k,i)
! contrib = dabs(new - ref)
! accu += contrib
! if(contrib .gt. 1.d-10)then
! print*,'cycle2'
! print*,i,k,j,l,m
! print*,ref,new,contrib
! stop
! endif
!
! new = three_e_5_idx_exch23_bi_ort(m,l,j,k,i)
! ref = three_e_5_idx_exch23_bi_ort_old(m,l,j,k,i)
! contrib = dabs(new - ref)
! accu += contrib
! if(contrib .gt. 1.d-10)then
! print*,'exch23'
! print*,i,k,j,l,m
! print*,ref,new,contrib
! stop
! endif
!
! new = three_e_5_idx_exch13_bi_ort(m,l,j,k,i)
! ref = three_e_5_idx_exch13_bi_ort_old(m,l,j,k,i)
! contrib = dabs(new - ref)
! accu += contrib
! if(contrib .gt. 1.d-10)then
! print*,'exch13'
! print*,i,k,j,l,m
! print*,ref,new,contrib
! stop
! endif
! !
new = three_e_5_idx_cycle_1_bi_ort(m,l,j,k,i)
ref = three_e_5_idx_cycle_1_bi_ort_old(m,l,j,k,i)
contrib = dabs(new - ref)
accu += contrib
if(contrib .gt. 1.d-10)then
print*,'cycle1'
print*,i,k,j,l,m
print*,ref,new,contrib
stop
endif
new = three_e_5_idx_cycle_2_bi_ort(m,l,j,k,i)
ref = three_e_5_idx_cycle_2_bi_ort_old(m,l,j,k,i)
contrib = dabs(new - ref)
accu += contrib
if(contrib .gt. 1.d-10)then
print*,'cycle2'
print*,i,k,j,l,m
print*,ref,new,contrib
stop
endif
new = three_e_5_idx_exch23_bi_ort(m,l,j,k,i)
ref = three_e_5_idx_exch23_bi_ort_old(m,l,j,k,i)
contrib = dabs(new - ref)
accu += contrib
if(contrib .gt. 1.d-10)then
print*,'exch23'
print*,i,k,j,l,m
print*,ref,new,contrib
stop
endif
new = three_e_5_idx_exch13_bi_ort(m,l,j,k,i)
ref = three_e_5_idx_exch13_bi_ort_old(m,l,j,k,i)
contrib = dabs(new - ref)
accu += contrib
if(contrib .gt. 1.d-10)then
print*,'exch13'
print*,i,k,j,l,m
print*,ref,new,contrib
stop
endif
enddo enddo
enddo enddo
enddo enddo

10
src/bi_ort_ints/three_body_ijmkl.irp.f
View File

@ -1,7 +1,11 @@
! --- ! ---
double precision function three_e_5_idx_exch12_bi_ort(m,l,i,k,j) result(integral)
implicit none
integer, intent(in) :: m,l,j,k,i
integral = three_e_5_idx_direct_bi_ort(m,l,j,k,i)
end
BEGIN_PROVIDER [ double precision, three_e_5_idx_direct_bi_ort , (mo_num, mo_num, mo_num, mo_num, mo_num)] BEGIN_PROVIDER [ double precision, three_e_5_idx_direct_bi_ort , (mo_num, mo_num, mo_num, mo_num, mo_num)]
&BEGIN_PROVIDER [ double precision, three_e_5_idx_exch12_bi_ort , (mo_num, mo_num, mo_num, mo_num, mo_num)]
&BEGIN_PROVIDER [ double precision, three_e_5_idx_exch23_bi_ort , (mo_num, mo_num, mo_num, mo_num, mo_num)] &BEGIN_PROVIDER [ double precision, three_e_5_idx_exch23_bi_ort , (mo_num, mo_num, mo_num, mo_num, mo_num)]
&BEGIN_PROVIDER [ double precision, three_e_5_idx_exch13_bi_ort , (mo_num, mo_num, mo_num, mo_num, mo_num)] &BEGIN_PROVIDER [ double precision, three_e_5_idx_exch13_bi_ort , (mo_num, mo_num, mo_num, mo_num, mo_num)]
&BEGIN_PROVIDER [ double precision, three_e_5_idx_cycle_1_bi_ort, (mo_num, mo_num, mo_num, mo_num, mo_num)] &BEGIN_PROVIDER [ double precision, three_e_5_idx_cycle_1_bi_ort, (mo_num, mo_num, mo_num, mo_num, mo_num)]
@ -14,6 +18,8 @@
! three_e_5_idx_direct_bi_ort(m,l,j,k,i) = <mlk|-L|mji> ::: notice that i is the RIGHT MO and k is the LEFT MO ! three_e_5_idx_direct_bi_ort(m,l,j,k,i) = <mlk|-L|mji> ::: notice that i is the RIGHT MO and k is the LEFT MO
! !
! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign ! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign
!
END_DOC END_DOC
implicit none implicit none
@ -74,7 +80,6 @@
do j = 1, mo_num do j = 1, mo_num
do l = 1, mo_num do l = 1, mo_num
three_e_5_idx_direct_bi_ort(m,l,j,k,i) = - tmp_mat(l,j,k,i) - tmp_mat(k,i,l,j) three_e_5_idx_direct_bi_ort(m,l,j,k,i) = - tmp_mat(l,j,k,i) - tmp_mat(k,i,l,j)
three_e_5_idx_exch12_bi_ort(m,l,j,k,i) = - tmp_mat(l,i,k,j) - tmp_mat(k,j,l,i)
enddo enddo
enddo enddo
enddo enddo
@ -128,7 +133,6 @@
do j = 1, mo_num do j = 1, mo_num
do l = 1, mo_num do l = 1, mo_num
three_e_5_idx_direct_bi_ort(m,l,j,k,i) = three_e_5_idx_direct_bi_ort(m,l,j,k,i) - tmp_mat(l,j,k,i) three_e_5_idx_direct_bi_ort(m,l,j,k,i) = three_e_5_idx_direct_bi_ort(m,l,j,k,i) - tmp_mat(l,j,k,i)
three_e_5_idx_exch12_bi_ort(m,l,j,k,i) = three_e_5_idx_exch12_bi_ort(m,l,j,k,i) - tmp_mat(l,i,k,j)
enddo enddo
enddo enddo
enddo enddo

34
src/tc_bi_ortho/31.tc_bi_ortho.bats
View File

@ -4,46 +4,50 @@ source $QP_ROOT/tests/bats/common.bats.sh
source $QP_ROOT/quantum_package.rc source $QP_ROOT/quantum_package.rc
function get_e() {
grep "eigval_right_tc_bi_orth" $1 | cut -d '=' -f 2 | xargs
}
function run_Ne() { function run_Ne() {
qp set_file Ne_tc_scf qp set_file Ne_tc_scf
qp run cisd qp run cisd
qp run tc_bi_ortho | tee Ne_tc_scf.cisd_tc_bi_ortho.out qp run tc_bi_ortho | tee Ne_tc_scf.cisd_tc_bi_ortho.out
eref=-128.77020441279302 eref=-128.77020441279302
energy="$(grep "eigval_right_tc_bi_orth =" Ne_tc_scf.cisd_tc_bi_ortho.out)" energy=$(get_e Ne_tc_scf.cisd_tc_bi_ortho.out)
eq $energy $eref 1e-6 eq $energy $eref 1e-6
} }
@test "Ne" { @test "Ne" {
run_Ne run_Ne
} }
function run_C() { function run_C() {
qp set_file C_tc_scf qp set_file C_tc_scf
qp run cisd qp run cisd
qp run tc_bi_ortho | tee C_tc_scf.cisd_tc_bi_ortho.out qp run tc_bi_ortho | tee C_tc_scf.cisd_tc_bi_ortho.out
eref=-37.757536149952514 eref=-37.757536149952514
energy="$(grep "eigval_right_tc_bi_orth =" C_tc_scf.cisd_tc_bi_ortho.out)" energy=$(get_e C_tc_scf.cisd_tc_bi_ortho.out)
eq $energy $eref 1e-6 eq $energy $eref 1e-6
} }
@test "C" { @test "C" {
run_C run_C
} }
function run_O() { function run_O() {
qp set_file C_tc_scf qp set_file C_tc_scf
qp run cisd qp run cisd
qp run tc_bi_ortho | tee O_tc_scf.cisd_tc_bi_ortho.out qp run tc_bi_ortho | tee O_tc_scf.cisd_tc_bi_ortho.out
eref=-74.908518517716161 eref=-74.908518517716161
energy="$(grep "eigval_right_tc_bi_orth =" O_tc_scf.cisd_tc_bi_ortho.out)" energy=$(get_e O_tc_scf.cisd_tc_bi_ortho.out)
eq $energy $eref 1e-6 eq $energy $eref 1e-6
} }
@test "O" { @test "O" {
run_O run_O
} }

89
src/tc_bi_ortho/slater_tc_3e_slow.irp.f
View File

@ -32,28 +32,28 @@ subroutine diag_htilde_three_body_ints_bi_ort_slow(Nint, key_i, hthree)
if(Ne(1)+Ne(2).ge.3)then if(Ne(1)+Ne(2).ge.3)then
!! ! alpha/alpha/beta three-body !! ! alpha/alpha/beta three-body
do i = 1, Ne(1) do i = 1, Ne(1)
ii = occ(i,1) ii = occ(i,1)
do j = i+1, Ne(1) do j = i+1, Ne(1)
jj = occ(j,1) jj = occ(j,1)
do m = 1, Ne(2) do m = 1, Ne(2)
mm = occ(m,2) mm = occ(m,2)
! direct_int = three_body_ints_bi_ort(mm,jj,ii,mm,jj,ii) USES THE 6-IDX TENSOR ! direct_int = three_body_ints_bi_ort(mm,jj,ii,mm,jj,ii) USES THE 6-IDX TENSOR
! exchange_int = three_body_ints_bi_ort(mm,jj,ii,mm,ii,jj) USES THE 6-IDX TENSOR ! exchange_int = three_body_ints_bi_ort(mm,jj,ii,mm,ii,jj) USES THE 6-IDX TENSOR
direct_int = three_e_3_idx_direct_bi_ort(mm,jj,ii) ! USES 3-IDX TENSOR direct_int = three_e_3_idx_direct_bi_ort(mm,jj,ii) ! USES 3-IDX TENSOR
exchange_int = three_e_3_idx_exch12_bi_ort(mm,jj,ii) ! USES 3-IDX TENSOR exchange_int = three_e_3_idx_exch12_bi_ort(mm,jj,ii) ! USES 3-IDX TENSOR
hthree += direct_int - exchange_int hthree += direct_int - exchange_int
enddo enddo
enddo enddo
enddo enddo
! beta/beta/alpha three-body ! beta/beta/alpha three-body
do i = 1, Ne(2) do i = 1, Ne(2)
ii = occ(i,2) ii = occ(i,2)
do j = i+1, Ne(2) do j = i+1, Ne(2)
jj = occ(j,2) jj = occ(j,2)
do m = 1, Ne(1) do m = 1, Ne(1)
mm = occ(m,1) mm = occ(m,1)
direct_int = three_e_3_idx_direct_bi_ort(mm,jj,ii) direct_int = three_e_3_idx_direct_bi_ort(mm,jj,ii)
exchange_int = three_e_3_idx_exch12_bi_ort(mm,jj,ii) exchange_int = three_e_3_idx_exch12_bi_ort(mm,jj,ii)
hthree += direct_int - exchange_int hthree += direct_int - exchange_int
enddo enddo
@ -64,10 +64,10 @@ subroutine diag_htilde_three_body_ints_bi_ort_slow(Nint, key_i, hthree)
do i = 1, Ne(1) do i = 1, Ne(1)
ii = occ(i,1) ! 1 ii = occ(i,1) ! 1
do j = i+1, Ne(1) do j = i+1, Ne(1)
jj = occ(j,1) ! 2 jj = occ(j,1) ! 2
do m = j+1, Ne(1) do m = j+1, Ne(1)
mm = occ(m,1) ! 3 mm = occ(m,1) ! 3
! ref = sym_3_e_int_from_6_idx_tensor(mm,jj,ii,mm,jj,ii) USES THE 6 IDX TENSOR ! ref = sym_3_e_int_from_6_idx_tensor(mm,jj,ii,mm,jj,ii) USES THE 6 IDX TENSOR
hthree += three_e_diag_parrallel_spin(mm,jj,ii) ! USES ONLY 3-IDX TENSORS hthree += three_e_diag_parrallel_spin(mm,jj,ii) ! USES ONLY 3-IDX TENSORS
enddo enddo
enddo enddo
@ -80,7 +80,7 @@ subroutine diag_htilde_three_body_ints_bi_ort_slow(Nint, key_i, hthree)
jj = occ(j,2) ! 2 jj = occ(j,2) ! 2
do m = j+1, Ne(2) do m = j+1, Ne(2)
mm = occ(m,2) ! 3 mm = occ(m,2) ! 3
! ref = sym_3_e_int_from_6_idx_tensor(mm,jj,ii,mm,jj,ii) USES THE 6 IDX TENSOR ! ref = sym_3_e_int_from_6_idx_tensor(mm,jj,ii,mm,jj,ii) USES THE 6 IDX TENSOR
hthree += three_e_diag_parrallel_spin(mm,jj,ii) ! USES ONLY 3-IDX TENSORS hthree += three_e_diag_parrallel_spin(mm,jj,ii) ! USES ONLY 3-IDX TENSORS
enddo enddo
enddo enddo
@ -96,7 +96,7 @@ subroutine single_htilde_three_body_ints_bi_ort_slow(Nint, key_j, key_i, hthree)
! <key_j | H_tilde | key_i> for single excitation ONLY FOR THREE-BODY TERMS WITH BI ORTHONORMAL ORBITALS ! <key_j | H_tilde | key_i> for single excitation ONLY FOR THREE-BODY TERMS WITH BI ORTHONORMAL ORBITALS
!! !!
!! WARNING !! !! WARNING !!
! !
! Non hermitian !! ! Non hermitian !!
END_DOC END_DOC
@ -110,7 +110,7 @@ subroutine single_htilde_three_body_ints_bi_ort_slow(Nint, key_j, key_i, hthree)
integer :: Ne(2),i,j,ii,jj,ispin,jspin,k,kk integer :: Ne(2),i,j,ii,jj,ispin,jspin,k,kk
integer :: degree,exc(0:2,2,2) integer :: degree,exc(0:2,2,2)
integer :: h1, p1, h2, p2, s1, s2 integer :: h1, p1, h2, p2, s1, s2
double precision :: direct_int,phase,exchange_int,three_e_single_parrallel_spin double precision :: direct_int,phase,exchange_int,three_e_single_parrallel_spin
double precision :: sym_3_e_int_from_6_idx_tensor double precision :: sym_3_e_int_from_6_idx_tensor
integer :: other_spin(2) integer :: other_spin(2)
integer(bit_kind) :: key_j_core(Nint,2),key_i_core(Nint,2) integer(bit_kind) :: key_j_core(Nint,2),key_i_core(Nint,2)
@ -142,26 +142,26 @@ subroutine single_htilde_three_body_ints_bi_ort_slow(Nint, key_j, key_i, hthree)
! alpha/alpha/beta three-body ! alpha/alpha/beta three-body
! print*,'IN SLAT RULES' ! print*,'IN SLAT RULES'
if(Ne(1)+Ne(2).ge.3)then if(Ne(1)+Ne(2).ge.3)then
! hole of spin s1 :: contribution from purely other spin ! hole of spin s1 :: contribution from purely other spin
ispin = other_spin(s1) ! ispin is the other spin than s1 ispin = other_spin(s1) ! ispin is the other spin than s1
do i = 1, Ne(ispin) ! i is the orbitals of the other spin than s1 do i = 1, Ne(ispin) ! i is the orbitals of the other spin than s1
ii = occ(i,ispin) ii = occ(i,ispin)
do j = i+1, Ne(ispin) ! j has the same spin than s1 do j = i+1, Ne(ispin) ! j has the same spin than s1
jj = occ(j,ispin) jj = occ(j,ispin)
! 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_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
enddo enddo
enddo enddo
! hole of spin s1 :: contribution from mixed other spin / same spin ! hole of spin s1 :: contribution from mixed other spin / same spin
do i = 1, Ne(ispin) ! other spin do i = 1, Ne(ispin) ! other spin
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_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 >
@ -174,8 +174,8 @@ subroutine single_htilde_three_body_ints_bi_ort_slow(Nint, key_j, key_i, hthree)
ii = occ(i,s1) ii = occ(i,s1)
do j = i+1, Ne(s1) do j = i+1, Ne(s1)
jj = occ(j,s1) jj = occ(j,s1)
! ref = sym_3_e_int_from_6_idx_tensor(jj,ii,p1,jj,ii,h1) ! ref = sym_3_e_int_from_6_idx_tensor(jj,ii,p1,jj,ii,h1)
hthree += three_e_single_parrallel_spin(jj,ii,p1,h1) ! USES THE 4-IDX TENSOR hthree += three_e_single_parrallel_spin(jj,ii,p1,h1) ! USES THE 4-IDX TENSOR
enddo enddo
enddo enddo
endif endif
@ -191,7 +191,7 @@ subroutine double_htilde_three_body_ints_bi_ort_slow(Nint, key_j, key_i, hthree)
! <key_j | H_tilde | key_i> for double excitation ONLY FOR THREE-BODY TERMS WITH BI ORTHONORMAL ORBITALS ! <key_j | H_tilde | key_i> for double excitation ONLY FOR THREE-BODY TERMS WITH BI ORTHONORMAL ORBITALS
!! !!
!! WARNING !! !! WARNING !!
! !
! Non hermitian !! ! Non hermitian !!
END_DOC END_DOC
@ -235,29 +235,30 @@ subroutine double_htilde_three_body_ints_bi_ort_slow(Nint, key_j, key_i, hthree)
call get_double_excitation(key_i, key_j, exc, phase, Nint) call get_double_excitation(key_i, key_j, exc, phase, Nint)
call decode_exc(exc, 2, h1, p1, h2, p2, s1, s2) call decode_exc(exc, 2, h1, p1, h2, p2, s1, s2)
if(Ne(1)+Ne(2).ge.3)then if(Ne(1)+Ne(2).ge.3)then
if(s1==s2)then ! same spin excitation if(s1==s2)then ! same spin excitation
ispin = other_spin(s1) ispin = other_spin(s1)
do m = 1, Ne(ispin) ! direct(other_spin) - exchange(s1) do m = 1, Ne(ispin) ! direct(other_spin) - exchange(s1)
mm = occ(m,ispin) mm = occ(m,ispin)
direct_int = three_e_5_idx_direct_bi_ort(mm,p2,h2,p1,h1) direct_int = three_e_5_idx_direct_bi_ort(mm,p2,h2,p1,h1)
exchange_int = three_e_5_idx_exch12_bi_ort(mm,p2,h2,p1,h1) ! exchange_int = three_e_5_idx_exch12_bi_ort(mm,p2,h2,p1,h1)
exchange_int = three_e_5_idx_direct_bi_ort(mm,p2,h1,p1,h2)
hthree += direct_int - exchange_int hthree += direct_int - exchange_int
enddo enddo
do m = 1, Ne(s1) ! pure contribution from s1 do m = 1, Ne(s1) ! pure contribution from s1
mm = occ(m,s1) mm = occ(m,s1)
hthree += three_e_double_parrallel_spin(mm,p2,h2,p1,h1) hthree += three_e_double_parrallel_spin(mm,p2,h2,p1,h1)
enddo enddo
else ! different spin excitation else ! different spin excitation
do m = 1, Ne(s1) do m = 1, Ne(s1)
mm = occ(m,s1) ! mm = occ(m,s1) !
direct_int = three_e_5_idx_direct_bi_ort(mm,p2,h2,p1,h1) direct_int = three_e_5_idx_direct_bi_ort(mm,p2,h2,p1,h1)
exchange_int = three_e_5_idx_exch13_bi_ort(mm,p2,h2,p1,h1) exchange_int = three_e_5_idx_exch13_bi_ort(mm,p2,h2,p1,h1)
hthree += direct_int - exchange_int hthree += direct_int - exchange_int
enddo enddo
do m = 1, Ne(s2) do m = 1, Ne(s2)
mm = occ(m,s2) ! mm = occ(m,s2) !
direct_int = three_e_5_idx_direct_bi_ort(mm,p2,h2,p1,h1) direct_int = three_e_5_idx_direct_bi_ort(mm,p2,h2,p1,h1)
exchange_int = three_e_5_idx_exch23_bi_ort(mm,p2,h2,p1,h1) exchange_int = three_e_5_idx_exch23_bi_ort(mm,p2,h2,p1,h1)
hthree += direct_int - exchange_int hthree += direct_int - exchange_int

3
src/tc_bi_ortho/slater_tc_opt.irp.f
View File

@ -13,8 +13,7 @@ subroutine provide_all_three_ints_bi_ortho
PROVIDE three_e_4_idx_exch23_bi_ort three_e_4_idx_exch13_bi_ort three_e_4_idx_exch12_bi_ort PROVIDE three_e_4_idx_exch23_bi_ort three_e_4_idx_exch13_bi_ort three_e_4_idx_exch12_bi_ort
endif endif
if(.not.double_normal_ord.and.three_e_5_idx_term)then if(.not.double_normal_ord.and.three_e_5_idx_term)then
PROVIDE three_e_5_idx_direct_bi_ort three_e_5_idx_cycle_1_bi_ort three_e_5_idx_cycle_2_bi_ort PROVIDE three_e_5_idx_direct_bi_ort
PROVIDE three_e_5_idx_exch23_bi_ort three_e_5_idx_exch13_bi_ort three_e_5_idx_exch12_bi_ort
elseif (double_normal_ord .and. (.not. three_e_5_idx_term))then elseif (double_normal_ord .and. (.not. three_e_5_idx_term))then
PROVIDE normal_two_body_bi_orth PROVIDE normal_two_body_bi_orth
endif endif

180
src/tc_bi_ortho/slater_tc_opt_double.irp.f
View File

@ -2,17 +2,17 @@
subroutine double_htilde_mu_mat_fock_bi_ortho(Nint, key_j, key_i, hmono, htwoe, hthree, htot) subroutine double_htilde_mu_mat_fock_bi_ortho(Nint, key_j, key_i, hmono, htwoe, hthree, htot)
BEGIN_DOC BEGIN_DOC
! <key_j | H_tilde | key_i> for double excitation ONLY FOR ONE- AND TWO-BODY TERMS ! <key_j | H_tilde | key_i> for double excitation ONLY FOR ONE- AND TWO-BODY TERMS
!! !!
!! WARNING !! !! WARNING !!
! !
! Non hermitian !! ! Non hermitian !!
END_DOC END_DOC
use bitmasks use bitmasks
implicit none implicit none
integer, intent(in) :: Nint integer, intent(in) :: Nint
integer(bit_kind), intent(in) :: key_j(Nint,2), key_i(Nint,2) integer(bit_kind), intent(in) :: key_j(Nint,2), key_i(Nint,2)
double precision, intent(out) :: hmono, htwoe, hthree, htot double precision, intent(out) :: hmono, htwoe, hthree, htot
integer :: occ(Nint*bit_kind_size,2) integer :: occ(Nint*bit_kind_size,2)
@ -39,8 +39,8 @@ subroutine double_htilde_mu_mat_fock_bi_ortho(Nint, key_j, key_i, hmono, htwoe,
call decode_exc(exc, 2, h1, p1, h2, p2, s1, s2) call decode_exc(exc, 2, h1, p1, h2, p2, s1, s2)
if(s1.ne.s2)then if(s1.ne.s2)then
! opposite spin two-body ! opposite spin two-body
htwoe = mo_bi_ortho_tc_two_e(p2,p1,h2,h1) htwoe = mo_bi_ortho_tc_two_e(p2,p1,h2,h1)
if(three_body_h_tc.and.elec_num.gt.2)then if(three_body_h_tc.and.elec_num.gt.2)then
if(.not.double_normal_ord.and.three_e_5_idx_term)then if(.not.double_normal_ord.and.three_e_5_idx_term)then
if(degree_i>degree_j)then if(degree_i>degree_j)then
@ -53,11 +53,11 @@ subroutine double_htilde_mu_mat_fock_bi_ortho(Nint, key_j, key_i, hmono, htwoe,
endif endif
endif endif
else else
! same spin two-body ! same spin two-body
! direct terms ! direct terms
htwoe = mo_bi_ortho_tc_two_e(p2,p1,h2,h1) htwoe = mo_bi_ortho_tc_two_e(p2,p1,h2,h1)
! exchange terms ! exchange terms
htwoe -= mo_bi_ortho_tc_two_e(p1,p2,h2,h1) htwoe -= mo_bi_ortho_tc_two_e(p1,p2,h2,h1)
if(three_body_h_tc.and.elec_num.gt.2)then if(three_body_h_tc.and.elec_num.gt.2)then
if(.not.double_normal_ord.and.three_e_5_idx_term)then if(.not.double_normal_ord.and.three_e_5_idx_term)then
if(degree_i>degree_j)then if(degree_i>degree_j)then
@ -112,72 +112,76 @@ subroutine three_comp_two_e_elem(key_i,h1,h2,p1,p2,s1,s2,hthree)
!DIR$ FORCEINLINE !DIR$ FORCEINLINE
call bitstring_to_list_ab(particle, occ_particle, tmp, N_int) call bitstring_to_list_ab(particle, occ_particle, tmp, N_int)
ASSERT (tmp(1) == nexc(1)) ! Number of particles alpha ASSERT (tmp(1) == nexc(1)) ! Number of particles alpha
ASSERT (tmp(2) == nexc(2)) ! Number of particle beta ASSERT (tmp(2) == nexc(2)) ! Number of particle beta
!DIR$ FORCEINLINE !DIR$ FORCEINLINE
call bitstring_to_list_ab(hole, occ_hole, tmp, N_int) call bitstring_to_list_ab(hole, occ_hole, tmp, N_int)
ASSERT (tmp(1) == nexc(1)) ! Number of holes alpha ASSERT (tmp(1) == nexc(1)) ! Number of holes alpha
ASSERT (tmp(2) == nexc(2)) ! Number of holes beta ASSERT (tmp(2) == nexc(2)) ! Number of holes beta
if(s1==s2.and.s1==1)then if(s1==s2.and.s1==1)then
!!!!!!!!!!!!!!!!!!!!!!!!!! alpha/alpha double exc !!!!!!!!!!!!!!!!!!!!!!!!!! alpha/alpha double exc
hthree = eff_2_e_from_3_e_aa(p2,p1,h2,h1) hthree = eff_2_e_from_3_e_aa(p2,p1,h2,h1)
if(nexc(1)+nexc(2) ==0)return !! if you're on the reference determinant if(nexc(1)+nexc(2) ==0)return !! if you're on the reference determinant
!!!!!!!! the matrix element is already exact !!!!!!!! the matrix element is already exact
!!!!!!!! else you need to take care of holes and particles !!!!!!!! else you need to take care of holes and particles
!!!!!!!!!!!!! Holes and particles !!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!! Holes and particles !!!!!!!!!!!!!!!!!!!!!!!
ispin = 1 ! i==alpha ==> pure same spin terms ispin = 1 ! i==alpha ==> pure same spin terms
do i = 1, nexc(ispin) ! number of couple of holes/particles do i = 1, nexc(ispin) ! number of couple of holes/particles
ipart=occ_particle(i,ispin) ipart=occ_particle(i,ispin)
hthree += three_e_double_parrallel_spin_prov(ipart,p2,h2,p1,h1) hthree += three_e_double_parrallel_spin_prov(ipart,p2,h2,p1,h1)
ihole=occ_hole(i,ispin) ihole=occ_hole(i,ispin)
hthree -= three_e_double_parrallel_spin_prov(ihole,p2,h2,p1,h1) hthree -= three_e_double_parrallel_spin_prov(ihole,p2,h2,p1,h1)
enddo enddo
ispin = 2 ! i==beta ==> alpha/alpha/beta terms ispin = 2 ! i==beta ==> alpha/alpha/beta terms
do i = 1, nexc(ispin) ! number of couple of holes/particles do i = 1, nexc(ispin) ! number of couple of holes/particles
! exchange between (h1,p1) and (h2,p2) ! exchange between (h1,p1) and (h2,p2)
ipart=occ_particle(i,ispin) ipart=occ_particle(i,ispin)
direct_int = three_e_5_idx_direct_bi_ort(ipart,p2,h2,p1,h1) direct_int = three_e_5_idx_direct_bi_ort(ipart,p2,h2,p1,h1)
exchange_int = three_e_5_idx_exch12_bi_ort(ipart,p2,h2,p1,h1) ! exchange_int = three_e_5_idx_exch12_bi_ort(ipart,p2,h2,p1,h1)
exchange_int = three_e_5_idx_direct_bi_ort(ipart,p2,h1,p1,h2)
hthree += direct_int - exchange_int hthree += direct_int - exchange_int
ihole=occ_hole(i,ispin) ihole=occ_hole(i,ispin)
direct_int = three_e_5_idx_direct_bi_ort(ihole,p2,h2,p1,h1) direct_int = three_e_5_idx_direct_bi_ort(ihole,p2,h2,p1,h1)
exchange_int = three_e_5_idx_exch12_bi_ort(ihole,p2,h2,p1,h1) ! exchange_int = three_e_5_idx_exch12_bi_ort(ihole,p2,h2,p1,h1)
exchange_int = three_e_5_idx_direct_bi_ort(ihole,p2,h1,p1,h2)
hthree -= direct_int - exchange_int hthree -= direct_int - exchange_int
enddo enddo
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
elseif(s1==s2.and.s1==2)then elseif(s1==s2.and.s1==2)then
!!!!!!!!!!!!!!!!!!!!!!!!!! beta/beta double exc !!!!!!!!!!!!!!!!!!!!!!!!!! beta/beta double exc
hthree = eff_2_e_from_3_e_bb(p2,p1,h2,h1) hthree = eff_2_e_from_3_e_bb(p2,p1,h2,h1)
if(nexc(1)+nexc(2) ==0)return !! if you're on the reference determinant if(nexc(1)+nexc(2) ==0)return !! if you're on the reference determinant
!!!!!!!! the matrix element is already exact !!!!!!!! the matrix element is already exact
!!!!!!!! else you need to take care of holes and particles !!!!!!!! else you need to take care of holes and particles
!!!!!!!!!!!!! Holes and particles !!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!! Holes and particles !!!!!!!!!!!!!!!!!!!!!!!
ispin = 2 ! i==beta ==> pure same spin terms ispin = 2 ! i==beta ==> pure same spin terms
do i = 1, nexc(ispin) ! number of couple of holes/particles do i = 1, nexc(ispin) ! number of couple of holes/particles
ipart=occ_particle(i,ispin) ipart=occ_particle(i,ispin)
hthree += three_e_double_parrallel_spin_prov(ipart,p2,h2,p1,h1) hthree += three_e_double_parrallel_spin_prov(ipart,p2,h2,p1,h1)
ihole=occ_hole(i,ispin) ihole=occ_hole(i,ispin)
hthree -= three_e_double_parrallel_spin_prov(ihole,p2,h2,p1,h1) hthree -= three_e_double_parrallel_spin_prov(ihole,p2,h2,p1,h1)
enddo enddo
ispin = 1 ! i==alpha==> beta/beta/alpha terms ispin = 1 ! i==alpha==> beta/beta/alpha terms
do i = 1, nexc(ispin) ! number of couple of holes/particles do i = 1, nexc(ispin) ! number of couple of holes/particles
! exchange between (h1,p1) and (h2,p2) ! exchange between (h1,p1) and (h2,p2)
ipart=occ_particle(i,ispin) ipart=occ_particle(i,ispin)
direct_int = three_e_5_idx_direct_bi_ort(ipart,p2,h2,p1,h1) direct_int = three_e_5_idx_direct_bi_ort(ipart,p2,h2,p1,h1)
exchange_int = three_e_5_idx_exch12_bi_ort(ipart,p2,h2,p1,h1) ! exchange_int = three_e_5_idx_exch12_bi_ort(ipart,p2,h2,p1,h1)
exchange_int = three_e_5_idx_direct_bi_ort(ipart,p2,h1,p1,h2)
hthree += direct_int - exchange_int hthree += direct_int - exchange_int
ihole=occ_hole(i,ispin) ihole=occ_hole(i,ispin)
direct_int = three_e_5_idx_direct_bi_ort(ihole,p2,h2,p1,h1) direct_int = three_e_5_idx_direct_bi_ort(ihole,p2,h2,p1,h1)
exchange_int = three_e_5_idx_exch12_bi_ort(ihole,p2,h2,p1,h1) ! exchange_int = three_e_5_idx_exch12_bi_ort(ihole,p2,h2,p1,h1)
exchange_int = three_e_5_idx_direct_bi_ort(ihole,p2,h1,p1,h2)
hthree -= direct_int - exchange_int hthree -= direct_int - exchange_int
enddo enddo
else ! (h1,p1) == alpha/(h2,p2) == beta else ! (h1,p1) == alpha/(h2,p2) == beta
hthree = eff_2_e_from_3_e_ab(p2,p1,h2,h1) hthree = eff_2_e_from_3_e_ab(p2,p1,h2,h1)
if(nexc(1)+nexc(2) ==0)return !! if you're on the reference determinant if(nexc(1)+nexc(2) ==0)return !! if you're on the reference determinant
!!!!!!!! the matrix element is already exact !!!!!!!! the matrix element is already exact
!!!!!!!! else you need to take care of holes and particles !!!!!!!! else you need to take care of holes and particles
!!!!!!!!!!!!! Holes and particles !!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!! Holes and particles !!!!!!!!!!!!!!!!!!!!!!!
ispin = 1 ! i==alpha ==> alpha/beta/alpha terms ispin = 1 ! i==alpha ==> alpha/beta/alpha terms
do i = 1, nexc(ispin) ! number of couple of holes/particles do i = 1, nexc(ispin) ! number of couple of holes/particles
! exchange between (h1,p1) and i ! exchange between (h1,p1) and i
ipart=occ_particle(i,ispin) ipart=occ_particle(i,ispin)
direct_int = three_e_5_idx_direct_bi_ort(ipart,p2,h2,p1,h1) direct_int = three_e_5_idx_direct_bi_ort(ipart,p2,h2,p1,h1)
@ -188,8 +192,8 @@ subroutine three_comp_two_e_elem(key_i,h1,h2,p1,p2,s1,s2,hthree)
exchange_int = three_e_5_idx_exch13_bi_ort(ihole,p2,h2,p1,h1) exchange_int = three_e_5_idx_exch13_bi_ort(ihole,p2,h2,p1,h1)
hthree -= direct_int - exchange_int hthree -= direct_int - exchange_int
enddo enddo
ispin = 2 ! i==beta ==> alpha/beta/beta terms ispin = 2 ! i==beta ==> alpha/beta/beta terms
do i = 1, nexc(ispin) ! number of couple of holes/particles do i = 1, nexc(ispin) ! number of couple of holes/particles
! exchange between (h2,p2) and i ! exchange between (h2,p2) and i
ipart=occ_particle(i,ispin) ipart=occ_particle(i,ispin)
direct_int = three_e_5_idx_direct_bi_ort(ipart,p2,h2,p1,h1) direct_int = three_e_5_idx_direct_bi_ort(ipart,p2,h2,p1,h1)
@ -207,7 +211,7 @@ end
BEGIN_PROVIDER [ double precision, eff_2_e_from_3_e_ab, (mo_num, mo_num, mo_num, mo_num)] BEGIN_PROVIDER [ double precision, eff_2_e_from_3_e_ab, (mo_num, mo_num, mo_num, mo_num)]
implicit none implicit none
BEGIN_DOC BEGIN_DOC
! eff_2_e_from_3_e_ab(p2,p1,h2,h1) = Effective Two-electron operator for alpha/beta double excitations ! eff_2_e_from_3_e_ab(p2,p1,h2,h1) = Effective Two-electron operator for alpha/beta double excitations
! !
! from contraction with HF density = a^{dagger}_p1_alpha a^{dagger}_p2_beta a_h2_beta a_h1_alpha ! from contraction with HF density = a^{dagger}_p1_alpha a^{dagger}_p2_beta a_h2_beta a_h1_alpha
END_DOC END_DOC
@ -222,16 +226,16 @@ BEGIN_PROVIDER [ double precision, eff_2_e_from_3_e_ab, (mo_num, mo_num, mo_num,
eff_2_e_from_3_e_ab = 0.d0 eff_2_e_from_3_e_ab = 0.d0
!$OMP PARALLEL & !$OMP PARALLEL &
!$OMP DEFAULT (NONE) & !$OMP DEFAULT (NONE) &
!$OMP PRIVATE (hh1, h1, hh2, h2, pp1, p1, pp2, p2, contrib) & !$OMP PRIVATE (hh1, h1, hh2, h2, pp1, p1, pp2, p2, contrib) &
!$OMP SHARED (n_act_orb, list_act, Ne,occ, eff_2_e_from_3_e_ab) !$OMP SHARED (n_act_orb, list_act, Ne,occ, eff_2_e_from_3_e_ab)
!$OMP DO SCHEDULE (static) !$OMP DO SCHEDULE (static)
do hh1 = 1, n_act_orb !! alpha do hh1 = 1, n_act_orb !! alpha
h1 = list_act(hh1) h1 = list_act(hh1)
do hh2 = 1, n_act_orb !! beta do hh2 = 1, n_act_orb !! beta
h2 = list_act(hh2) h2 = list_act(hh2)
do pp1 = 1, n_act_orb !! alpha do pp1 = 1, n_act_orb !! alpha
p1 = list_act(pp1) p1 = list_act(pp1)
do pp2 = 1, n_act_orb !! beta do pp2 = 1, n_act_orb !! beta
p2 = list_act(pp2) p2 = list_act(pp2)
call give_contrib_for_abab(h1,h2,p1,p2,occ,Ne,contrib) call give_contrib_for_abab(h1,h2,p1,p2,occ,Ne,contrib)
eff_2_e_from_3_e_ab(p2,p1,h2,h1) = contrib eff_2_e_from_3_e_ab(p2,p1,h2,h1) = contrib
@ -242,25 +246,25 @@ BEGIN_PROVIDER [ double precision, eff_2_e_from_3_e_ab, (mo_num, mo_num, mo_num,
!$OMP END DO !$OMP END DO
!$OMP END PARALLEL !$OMP END PARALLEL
END_PROVIDER END_PROVIDER
subroutine give_contrib_for_abab(h1,h2,p1,p2,occ,Ne,contrib) subroutine give_contrib_for_abab(h1,h2,p1,p2,occ,Ne,contrib)
implicit none implicit none
BEGIN_DOC BEGIN_DOC
! gives the contribution for a double excitation (h1,p1)_alpha (h2,p2)_beta ! gives the contribution for a double excitation (h1,p1)_alpha (h2,p2)_beta
! !
! on top of a determinant whose occupied orbitals is in (occ, Ne) ! on top of a determinant whose occupied orbitals is in (occ, Ne)
END_DOC END_DOC
integer, intent(in) :: h1,h2,p1,p2,occ(N_int*bit_kind_size,2),Ne(2) integer, intent(in) :: h1,h2,p1,p2,occ(N_int*bit_kind_size,2),Ne(2)
double precision, intent(out) :: contrib double precision, intent(out) :: contrib
integer :: mm,m integer :: mm,m
double precision :: direct_int, exchange_int double precision :: direct_int, exchange_int
!! h1,p1 == alpha !! h1,p1 == alpha
!! h2,p2 == beta !! h2,p2 == beta
contrib = 0.d0 contrib = 0.d0
do mm = 1, Ne(1) !! alpha do mm = 1, Ne(1) !! alpha
m = occ(mm,1) m = occ(mm,1)
direct_int = three_e_5_idx_direct_bi_ort(mm,p2,h2,p1,h1) direct_int = three_e_5_idx_direct_bi_ort(mm,p2,h2,p1,h1)
! exchange between (h1,p1) and m ! exchange between (h1,p1) and m
exchange_int = three_e_5_idx_exch13_bi_ort(mm,p2,h2,p1,h1) exchange_int = three_e_5_idx_exch13_bi_ort(mm,p2,h2,p1,h1)
contrib += direct_int - exchange_int contrib += direct_int - exchange_int
@ -268,7 +272,7 @@ subroutine give_contrib_for_abab(h1,h2,p1,p2,occ,Ne,contrib)
do mm = 1, Ne(2) !! beta do mm = 1, Ne(2) !! beta
m = occ(mm,2) m = occ(mm,2)
direct_int = three_e_5_idx_direct_bi_ort(mm,p2,h2,p1,h1) direct_int = three_e_5_idx_direct_bi_ort(mm,p2,h2,p1,h1)
! exchange between (h2,p2) and m ! exchange between (h2,p2) and m
exchange_int = three_e_5_idx_exch23_bi_ort(mm,p2,h2,p1,h1) exchange_int = three_e_5_idx_exch23_bi_ort(mm,p2,h2,p1,h1)
contrib += direct_int - exchange_int contrib += direct_int - exchange_int
@ -278,11 +282,11 @@ end
BEGIN_PROVIDER [ double precision, eff_2_e_from_3_e_aa, (mo_num, mo_num, mo_num, mo_num)] BEGIN_PROVIDER [ double precision, eff_2_e_from_3_e_aa, (mo_num, mo_num, mo_num, mo_num)]
implicit none implicit none
BEGIN_DOC BEGIN_DOC
! eff_2_e_from_3_e_ab(p2,p1,h2,h1) = Effective Two-electron operator for alpha/alpha double excitations ! eff_2_e_from_3_e_ab(p2,p1,h2,h1) = Effective Two-electron operator for alpha/alpha double excitations
! !
! from contractionelec_alpha_num with HF density = a^{dagger}_p1_alpha a^{dagger}_p2_alpha a_h2_alpha a_h1_alpha ! from contractionelec_alpha_num with HF density = a^{dagger}_p1_alpha a^{dagger}_p2_alpha a_h2_alpha a_h1_alpha
! !
! WARNING :: to be coherent with the phase convention used in the Hamiltonian matrix elements, you must fulfill ! WARNING :: to be coherent with the phase convention used in the Hamiltonian matrix elements, you must fulfill
! !
! |||| h2>h1, p2>p1 |||| ! |||| h2>h1, p2>p1 ||||
END_DOC END_DOC
@ -297,13 +301,13 @@ BEGIN_PROVIDER [ double precision, eff_2_e_from_3_e_aa, (mo_num, mo_num, mo_num,
eff_2_e_from_3_e_aa = 100000000.d0 eff_2_e_from_3_e_aa = 100000000.d0
!$OMP PARALLEL & !$OMP PARALLEL &
!$OMP DEFAULT (NONE) & !$OMP DEFAULT (NONE) &
!$OMP PRIVATE (hh1, h1, hh2, h2, pp1, p1, pp2, p2, contrib) & !$OMP PRIVATE (hh1, h1, hh2, h2, pp1, p1, pp2, p2, contrib) &
!$OMP SHARED (n_act_orb, list_act, Ne,occ, eff_2_e_from_3_e_aa) !$OMP SHARED (n_act_orb, list_act, Ne,occ, eff_2_e_from_3_e_aa)
!$OMP DO SCHEDULE (static) !$OMP DO SCHEDULE (static)
do hh1 = 1, n_act_orb !! alpha do hh1 = 1, n_act_orb !! alpha
h1 = list_act(hh1) h1 = list_act(hh1)
do hh2 = hh1+1, n_act_orb !! alpha do hh2 = hh1+1, n_act_orb !! alpha
h2 = list_act(hh2) h2 = list_act(hh2)
do pp1 = 1, n_act_orb !! alpha do pp1 = 1, n_act_orb !! alpha
p1 = list_act(pp1) p1 = list_act(pp1)
do pp2 = pp1+1, n_act_orb !! alpha do pp2 = pp1+1, n_act_orb !! alpha
@ -317,20 +321,20 @@ BEGIN_PROVIDER [ double precision, eff_2_e_from_3_e_aa, (mo_num, mo_num, mo_num,
!$OMP END DO !$OMP END DO
!$OMP END PARALLEL !$OMP END PARALLEL
END_PROVIDER END_PROVIDER
subroutine give_contrib_for_aaaa(h1,h2,p1,p2,occ,Ne,contrib) subroutine give_contrib_for_aaaa(h1,h2,p1,p2,occ,Ne,contrib)
implicit none implicit none
BEGIN_DOC BEGIN_DOC
! gives the contribution for a double excitation (h1,p1)_alpha (h2,p2)_alpha ! gives the contribution for a double excitation (h1,p1)_alpha (h2,p2)_alpha
! !
! on top of a determinant whose occupied orbitals is in (occ, Ne) ! on top of a determinant whose occupied orbitals is in (occ, Ne)
END_DOC END_DOC
integer, intent(in) :: h1,h2,p1,p2,occ(N_int*bit_kind_size,2),Ne(2) integer, intent(in) :: h1,h2,p1,p2,occ(N_int*bit_kind_size,2),Ne(2)
double precision, intent(out) :: contrib double precision, intent(out) :: contrib
integer :: mm,m integer :: mm,m
double precision :: direct_int, exchange_int double precision :: direct_int, exchange_int
!! h1,p1 == alpha !! h1,p1 == alpha
!! h2,p2 == alpha !! h2,p2 == alpha
contrib = 0.d0 contrib = 0.d0
do mm = 1, Ne(1) !! alpha ==> pure parallele spin contribution do mm = 1, Ne(1) !! alpha ==> pure parallele spin contribution
@ -340,9 +344,10 @@ subroutine give_contrib_for_aaaa(h1,h2,p1,p2,occ,Ne,contrib)
do mm = 1, Ne(2) !! beta do mm = 1, Ne(2) !! beta
m = occ(mm,2) m = occ(mm,2)
direct_int = three_e_5_idx_direct_bi_ort(mm,p2,h2,p1,h1) direct_int = three_e_5_idx_direct_bi_ort(mm,p2,h2,p1,h1)
! exchange between (h1,p1) and (h2,p2) ! exchange between (h1,p1) and (h2,p2)
exchange_int = three_e_5_idx_exch12_bi_ort(mm,p2,h2,p1,h1) ! exchange_int = three_e_5_idx_exch12_bi_ort(mm,p2,h2,p1,h1)
exchange_int = three_e_5_idx_direct_bi_ort(mm,p2,h1,p1,h2)
contrib += direct_int - exchange_int contrib += direct_int - exchange_int
enddo enddo
end end
@ -351,11 +356,11 @@ end
BEGIN_PROVIDER [ double precision, eff_2_e_from_3_e_bb, (mo_num, mo_num, mo_num, mo_num)] BEGIN_PROVIDER [ double precision, eff_2_e_from_3_e_bb, (mo_num, mo_num, mo_num, mo_num)]
implicit none implicit none
BEGIN_DOC BEGIN_DOC
! eff_2_e_from_3_e_ab(p2,p1,h2,h1) = Effective Two-electron operator for beta/beta double excitations ! eff_2_e_from_3_e_ab(p2,p1,h2,h1) = Effective Two-electron operator for beta/beta double excitations
! !
! from contractionelec_beta_num with HF density = a^{dagger}_p1_beta a^{dagger}_p2_beta a_h2_beta a_h1_beta ! from contractionelec_beta_num with HF density = a^{dagger}_p1_beta a^{dagger}_p2_beta a_h2_beta a_h1_beta
! !
! WARNING :: to be coherent with the phase convention used in the Hamiltonian matrix elements, you must fulfill ! WARNING :: to be coherent with the phase convention used in the Hamiltonian matrix elements, you must fulfill
! !
! |||| h2>h1, p2>p1 |||| ! |||| h2>h1, p2>p1 ||||
END_DOC END_DOC
@ -370,13 +375,13 @@ BEGIN_PROVIDER [ double precision, eff_2_e_from_3_e_bb, (mo_num, mo_num, mo_num,
eff_2_e_from_3_e_bb = 100000000.d0 eff_2_e_from_3_e_bb = 100000000.d0
!$OMP PARALLEL & !$OMP PARALLEL &
!$OMP DEFAULT (NONE) & !$OMP DEFAULT (NONE) &
!$OMP PRIVATE (hh1, h1, hh2, h2, pp1, p1, pp2, p2, contrib) & !$OMP PRIVATE (hh1, h1, hh2, h2, pp1, p1, pp2, p2, contrib) &
!$OMP SHARED (n_act_orb, list_act, Ne,occ, eff_2_e_from_3_e_bb) !$OMP SHARED (n_act_orb, list_act, Ne,occ, eff_2_e_from_3_e_bb)
!$OMP DO SCHEDULE (static) !$OMP DO SCHEDULE (static)
do hh1 = 1, n_act_orb !! beta do hh1 = 1, n_act_orb !! beta
h1 = list_act(hh1) h1 = list_act(hh1)
do hh2 = hh1+1, n_act_orb !! beta do hh2 = hh1+1, n_act_orb !! beta
h2 = list_act(hh2) h2 = list_act(hh2)
do pp1 = 1, n_act_orb !! beta do pp1 = 1, n_act_orb !! beta
p1 = list_act(pp1) p1 = list_act(pp1)
do pp2 = pp1+1, n_act_orb !! beta do pp2 = pp1+1, n_act_orb !! beta
@ -390,18 +395,18 @@ BEGIN_PROVIDER [ double precision, eff_2_e_from_3_e_bb, (mo_num, mo_num, mo_num,
!$OMP END DO !$OMP END DO
!$OMP END PARALLEL !$OMP END PARALLEL
END_PROVIDER END_PROVIDER
subroutine give_contrib_for_bbbb(h1,h2,p1,p2,occ,Ne,contrib) subroutine give_contrib_for_bbbb(h1,h2,p1,p2,occ,Ne,contrib)
implicit none implicit none
BEGIN_DOC BEGIN_DOC
! gives the contribution for a double excitation (h1,p1)_beta (h2,p2)_beta ! gives the contribution for a double excitation (h1,p1)_beta (h2,p2)_beta
! !
! on top of a determinant whose occupied orbitals is in (occ, Ne) ! on top of a determinant whose occupied orbitals is in (occ, Ne)
END_DOC END_DOC
integer, intent(in) :: h1,h2,p1,p2,occ(N_int*bit_kind_size,2),Ne(2) integer, intent(in) :: h1,h2,p1,p2,occ(N_int*bit_kind_size,2),Ne(2)
double precision, intent(out) :: contrib double precision, intent(out) :: contrib
integer :: mm,m integer :: mm,m
double precision :: direct_int, exchange_int double precision :: direct_int, exchange_int
!! h1,p1 == beta !! h1,p1 == beta
!! h2,p2 == beta !! h2,p2 == beta
@ -413,9 +418,10 @@ subroutine give_contrib_for_bbbb(h1,h2,p1,p2,occ,Ne,contrib)
do mm = 1, Ne(1) !! alpha do mm = 1, Ne(1) !! alpha
m = occ(mm,1) m = occ(mm,1)
direct_int = three_e_5_idx_direct_bi_ort(mm,p2,h2,p1,h1) direct_int = three_e_5_idx_direct_bi_ort(mm,p2,h2,p1,h1)
! exchange between (h1,p1) and (h2,p2) ! exchange between (h1,p1) and (h2,p2)
exchange_int = three_e_5_idx_exch12_bi_ort(mm,p2,h2,p1,h1) ! exchange_int = three_e_5_idx_exch12_bi_ort(mm,p2,h2,p1,h1)
exchange_int = three_e_5_idx_direct_bi_ort(mm,p2,h1,p1,h2)
contrib += direct_int - exchange_int contrib += direct_int - exchange_int
enddo enddo
end end
@ -424,17 +430,17 @@ end
subroutine double_htilde_mu_mat_fock_bi_ortho_no_3e(Nint, key_j, key_i, htot) subroutine double_htilde_mu_mat_fock_bi_ortho_no_3e(Nint, key_j, key_i, htot)
BEGIN_DOC BEGIN_DOC
! <key_j | H_tilde | key_i> for double excitation ONLY FOR ONE- AND TWO-BODY TERMS ! <key_j | H_tilde | key_i> for double excitation ONLY FOR ONE- AND TWO-BODY TERMS
!! !!
!! WARNING !! !! WARNING !!
! !
! Non hermitian !! ! Non hermitian !!
END_DOC END_DOC
use bitmasks use bitmasks
implicit none implicit none
integer, intent(in) :: Nint integer, intent(in) :: Nint
integer(bit_kind), intent(in) :: key_j(Nint,2), key_i(Nint,2) integer(bit_kind), intent(in) :: key_j(Nint,2), key_i(Nint,2)
double precision, intent(out) :: htot double precision, intent(out) :: htot
double precision :: hmono, htwoe double precision :: hmono, htwoe
@ -461,17 +467,17 @@ subroutine double_htilde_mu_mat_fock_bi_ortho_no_3e(Nint, key_j, key_i, htot)
call decode_exc(exc, 2, h1, p1, h2, p2, s1, s2) call decode_exc(exc, 2, h1, p1, h2, p2, s1, s2)
if(s1.ne.s2)then if(s1.ne.s2)then
! opposite spin two-body ! opposite spin two-body
htwoe = mo_bi_ortho_tc_two_e(p2,p1,h2,h1) htwoe = mo_bi_ortho_tc_two_e(p2,p1,h2,h1)
else else
! same spin two-body ! same spin two-body
! direct terms ! direct terms
htwoe = mo_bi_ortho_tc_two_e(p2,p1,h2,h1) htwoe = mo_bi_ortho_tc_two_e(p2,p1,h2,h1)
! exchange terms ! exchange terms
htwoe -= mo_bi_ortho_tc_two_e(p1,p2,h2,h1) htwoe -= mo_bi_ortho_tc_two_e(p1,p2,h2,h1)
endif endif
htwoe *= phase htwoe *= phase
htot = htwoe htot = htwoe
end end

3
src/tc_bi_ortho/symmetrized_3_e_int.irp.f
View File

@ -107,5 +107,6 @@ double precision function three_e_double_parrallel_spin(m,l,j,k,i)
three_e_double_parrallel_spin = three_e_5_idx_direct_bi_ort(m,l,j,k,i) ! direct three_e_double_parrallel_spin = three_e_5_idx_direct_bi_ort(m,l,j,k,i) ! direct
three_e_double_parrallel_spin += three_e_5_idx_cycle_1_bi_ort(m,l,j,k,i) + three_e_5_idx_cycle_2_bi_ort(m,l,j,k,i) & ! two cyclic permutations three_e_double_parrallel_spin += three_e_5_idx_cycle_1_bi_ort(m,l,j,k,i) + three_e_5_idx_cycle_2_bi_ort(m,l,j,k,i) & ! two cyclic permutations
- three_e_5_idx_exch23_bi_ort(m,l,j,k,i) - three_e_5_idx_exch13_bi_ort(m,l,j,k,i) & ! two first exchange - three_e_5_idx_exch23_bi_ort(m,l,j,k,i) - three_e_5_idx_exch13_bi_ort(m,l,j,k,i) & ! two first exchange
- three_e_5_idx_exch12_bi_ort(m,l,j,k,i) ! last exchange ! - three_e_5_idx_exch12_bi_ort(m,l,j,k,i) ! last exchange
- three_e_5_idx_direct_bi_ort(m,l,i,k,j) ! last exchange
end end