LCPQ
/
qp2
mirror of https://github.com/QuantumPackage/qp2.git
9
1
Fork 0
Code Releases Activity

Improve 5idx
continuous-integration/drone/push Build is failing Details

This commit is contained in:
Anthony Scemama 2023-06-02 20:32:31 +02:00
parent 81b7751b00
commit 4cc8dae420
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 ezfio import ezfio
from docopt import docopt
import qp_bitmasks
try:
import trexio
@ -453,6 +454,20 @@ def write_ezfio(trexio_filename, filename):
else:
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
AssertionError
"""
assert type(s) in (int, long)
assert type(s) == int
assert s>=0
return '{s:0b}'.format(s=s)
@ -62,7 +62,7 @@ def int_to_bitmask(s,bit_kind_size=BIT_KIND_SIZE):
['1111111111111111111111111111111111111111111111111111111111110110']
>>>
"""
assert type(s) in (int, long)
assert type(s) == int
if s < 0:
s = s + (1 << 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]
def __setitem__(self,i,value):
if type(value) in (int,long):
if type(value) == int :
self._data_int[i] = value
elif type(value) == str:
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
PROVIDE ezfio_filename
use_cosgtos = .False.
if (mpi_master) then
call ezfio_has_ao_basis_use_cosgtos(has)
if (has) then
! write(6,'(A)') '.. >>>>> [ IO READ: use_cosgtos ] <<<<< ..'
call ezfio_get_ao_basis_use_cosgtos(use_cosgtos)
else
use_cosgtos = .False.
call ezfio_set_ao_basis_use_cosgtos(use_cosgtos)
endif
endif
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
integer :: i,k,j,l,m,n,ipoint
double precision :: accu, contrib,new,ref
double precision, external :: three_e_5_idx_exch12_bi_ort
i = 1
k = 1
n = 0
@ -64,18 +65,21 @@ subroutine test_5idx
do j = 1, mo_num
do l = 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)
ref = three_e_5_idx_direct_bi_ort_old(m,l,j,k,i)
contrib = dabs(new - ref)
accu += contrib
if(contrib .gt. 1.d-10)then
print*,'direct'
print*,i,k,j,l,m
print*,ref,new,contrib
stop
endif
! 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)
! contrib = dabs(new - ref)
! accu += contrib
! if(contrib .gt. 1.d-10)then
! print*,'direct'
! print*,i,k,j,l,m
! print*,ref,new,contrib
! stop
! endif
!
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)
contrib = dabs(new - ref)
@ -86,51 +90,52 @@ subroutine test_5idx
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
!
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

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_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_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)]
@ -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
!
! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign
!
END_DOC
implicit none
@ -72,7 +78,6 @@
do j = 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_exch12_bi_ort(m,l,j,k,i) = - tmp_mat(l,i,k,j) - tmp_mat(k,j,l,i)
enddo
enddo
enddo
@ -125,7 +130,6 @@
do j = 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_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

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
function get_e() {
grep "eigval_right_tc_bi_orth" $1 | cut -d '=' -f 2 | xargs
}
function run_Ne() {
qp set_file Ne_tc_scf
qp run cisd
qp run tc_bi_ortho | tee Ne_tc_scf.cisd_tc_bi_ortho.out
qp set_file Ne_tc_scf
qp run cisd
qp run tc_bi_ortho | tee Ne_tc_scf.cisd_tc_bi_ortho.out
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
}
@test "Ne" {
run_Ne
run_Ne
}
function run_C() {
qp set_file C_tc_scf
qp run cisd
qp run tc_bi_ortho | tee C_tc_scf.cisd_tc_bi_ortho.out
qp set_file C_tc_scf
qp run cisd
qp run tc_bi_ortho | tee C_tc_scf.cisd_tc_bi_ortho.out
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
}
@test "C" {
run_C
run_C
}
function run_O() {
qp set_file C_tc_scf
qp run cisd
qp run tc_bi_ortho | tee O_tc_scf.cisd_tc_bi_ortho.out
qp set_file C_tc_scf
qp run cisd
qp run tc_bi_ortho | tee O_tc_scf.cisd_tc_bi_ortho.out
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
}
@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
!! ! alpha/alpha/beta three-body
do i = 1, Ne(1)
ii = occ(i,1)
ii = occ(i,1)
do j = i+1, Ne(1)
jj = occ(j,1)
jj = occ(j,1)
do m = 1, Ne(2)
mm = occ(m,2)
! direct_int = three_body_ints_bi_ort(mm,jj,ii,mm,jj,ii) USES THE 6-IDX TENSOR
! exchange_int = three_body_ints_bi_ort(mm,jj,ii,mm,ii,jj) USES THE 6-IDX TENSOR
direct_int = three_e_3_idx_direct_bi_ort(mm,jj,ii) ! USES 3-IDX TENSOR
exchange_int = three_e_3_idx_exch12_bi_ort(mm,jj,ii) ! USES 3-IDX TENSOR
mm = occ(m,2)
! direct_int = three_body_ints_bi_ort(mm,jj,ii,mm,jj,ii) USES THE 6-IDX TENSOR
! exchange_int = three_body_ints_bi_ort(mm,jj,ii,mm,ii,jj) USES THE 6-IDX TENSOR
direct_int = three_e_3_idx_direct_bi_ort(mm,jj,ii) ! USES 3-IDX TENSOR
exchange_int = three_e_3_idx_exch12_bi_ort(mm,jj,ii) ! USES 3-IDX TENSOR
hthree += direct_int - exchange_int
enddo
enddo
enddo
! beta/beta/alpha three-body
do i = 1, Ne(2)
ii = occ(i,2)
ii = occ(i,2)
do j = i+1, Ne(2)
jj = occ(j,2)
jj = occ(j,2)
do m = 1, Ne(1)
mm = occ(m,1)
direct_int = three_e_3_idx_direct_bi_ort(mm,jj,ii)
mm = occ(m,1)
direct_int = three_e_3_idx_direct_bi_ort(mm,jj,ii)
exchange_int = three_e_3_idx_exch12_bi_ort(mm,jj,ii)
hthree += direct_int - exchange_int
enddo
@ -64,10 +64,10 @@ subroutine diag_htilde_three_body_ints_bi_ort_slow(Nint, key_i, hthree)
do i = 1, Ne(1)
ii = occ(i,1) ! 1
do j = i+1, Ne(1)
jj = occ(j,1) ! 2
jj = occ(j,1) ! 2
do m = j+1, Ne(1)
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
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
hthree += three_e_diag_parrallel_spin(mm,jj,ii) ! USES ONLY 3-IDX TENSORS
enddo
enddo
@ -80,7 +80,7 @@ subroutine diag_htilde_three_body_ints_bi_ort_slow(Nint, key_i, hthree)
jj = occ(j,2) ! 2
do m = j+1, Ne(2)
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
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
!!
!! WARNING !!
!
!
! Non hermitian !!
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 :: degree,exc(0:2,2,2)
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
integer :: other_spin(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
! print*,'IN SLAT RULES'
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
do i = 1, Ne(ispin) ! i is the orbitals of the other spin than s1
ii = occ(i,ispin)
do j = i+1, Ne(ispin) ! j has the same spin than s1
jj = occ(j,ispin)
do i = 1, Ne(ispin) ! i is the orbitals of the other spin than s1
ii = occ(i,ispin)
do j = i+1, Ne(ispin) ! j has the same spin than s1
jj = occ(j,ispin)
! 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 >
!
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)
!
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)
hthree += direct_int - exchange_int
enddo
enddo
! hole of spin s1 :: contribution from mixed other spin / same spin
do i = 1, Ne(ispin) ! other spin
ii = occ(i,ispin) ! other spin
do j = 1, Ne(s1) ! same spin
jj = occ(j,s1) ! same spin
do i = 1, Ne(ispin) ! other spin
ii = occ(i,ispin) ! other spin
do j = 1, Ne(s1) ! same spin
jj = occ(j,s1) ! same spin
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)
! < 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)
do j = i+1, Ne(s1)
jj = occ(j,s1)
! 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
! 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
enddo
enddo
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
!!
!! WARNING !!
!
!
! Non hermitian !!
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 decode_exc(exc, 2, h1, p1, h2, p2, s1, s2)
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)
do m = 1, Ne(ispin) ! direct(other_spin) - exchange(s1)
mm = occ(m,ispin)
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)
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_direct_bi_ort(mm,p2,h1,p1,h2)
hthree += direct_int - exchange_int
enddo
do m = 1, Ne(s1) ! pure contribution from s1
do m = 1, Ne(s1) ! pure contribution from s1
mm = occ(m,s1)
hthree += three_e_double_parrallel_spin(mm,p2,h2,p1,h1)
enddo
else ! different spin excitation
enddo
else ! different spin excitation
do m = 1, Ne(s1)
mm = occ(m,s1) !
direct_int = three_e_5_idx_direct_bi_ort(mm,p2,h2,p1,h1)
mm = occ(m,s1) !
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)
hthree += direct_int - exchange_int
enddo
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)
exchange_int = three_e_5_idx_exch23_bi_ort(mm,p2,h2,p1,h1)
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
endif
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_exch23_bi_ort three_e_5_idx_exch13_bi_ort three_e_5_idx_exch12_bi_ort
PROVIDE three_e_5_idx_direct_bi_ort
elseif (double_normal_ord .and. (.not. three_e_5_idx_term))then
PROVIDE normal_two_body_bi_orth
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)
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 !!
!
!
! Non hermitian !!
END_DOC
use bitmasks
implicit none
integer, intent(in) :: Nint
integer, intent(in) :: Nint
integer(bit_kind), intent(in) :: key_j(Nint,2), key_i(Nint,2)
double precision, intent(out) :: hmono, htwoe, hthree, htot
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)
if(s1.ne.s2)then
! opposite spin two-body
htwoe = mo_bi_ortho_tc_two_e(p2,p1,h2,h1)
! opposite spin two-body
htwoe = mo_bi_ortho_tc_two_e(p2,p1,h2,h1)
if(three_body_h_tc.and.elec_num.gt.2)then
if(.not.double_normal_ord.and.three_e_5_idx_term)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
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)
! 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)
if(three_body_h_tc.and.elec_num.gt.2)then
if(.not.double_normal_ord.and.three_e_5_idx_term)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
call bitstring_to_list_ab(particle, occ_particle, tmp, N_int)
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
call bitstring_to_list_ab(hole, occ_hole, tmp, N_int)
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
!!!!!!!!!!!!!!!!!!!!!!!!!! alpha/alpha double exc
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
!!!!!!!! the matrix element is already exact
!!!!!!!! else you need to take care of holes and particles
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
!!!!!!!! the matrix element is already exact
!!!!!!!! else you need to take care of holes and particles
!!!!!!!!!!!!! Holes and particles !!!!!!!!!!!!!!!!!!!!!!!
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)
hthree += three_e_double_parrallel_spin_prov(ipart,p2,h2,p1,h1)
ihole=occ_hole(i,ispin)
hthree -= three_e_double_parrallel_spin_prov(ihole,p2,h2,p1,h1)
enddo
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)
ipart=occ_particle(i,ispin)
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
ihole=occ_hole(i,ispin)
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
enddo
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
elseif(s1==s2.and.s1==2)then
elseif(s1==s2.and.s1==2)then
!!!!!!!!!!!!!!!!!!!!!!!!!! beta/beta double exc
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
!!!!!!!! the matrix element is already exact
!!!!!!!! else you need to take care of holes and particles
if(nexc(1)+nexc(2) ==0)return !! if you're on the reference determinant
!!!!!!!! the matrix element is already exact
!!!!!!!! else you need to take care of holes and particles
!!!!!!!!!!!!! Holes and particles !!!!!!!!!!!!!!!!!!!!!!!
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)
hthree += three_e_double_parrallel_spin_prov(ipart,p2,h2,p1,h1)
ihole=occ_hole(i,ispin)
hthree -= three_e_double_parrallel_spin_prov(ihole,p2,h2,p1,h1)
enddo
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)
ipart=occ_particle(i,ispin)
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
ihole=occ_hole(i,ispin)
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
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)
if(nexc(1)+nexc(2) ==0)return !! if you're on the reference determinant
!!!!!!!! the matrix element is already exact
!!!!!!!! else you need to take care of holes and particles
if(nexc(1)+nexc(2) ==0)return !! if you're on the reference determinant
!!!!!!!! the matrix element is already exact
!!!!!!!! else you need to take care of holes and particles
!!!!!!!!!!!!! Holes and particles !!!!!!!!!!!!!!!!!!!!!!!
ispin = 1 ! i==alpha ==> alpha/beta/alpha terms
do i = 1, nexc(ispin) ! number of couple of holes/particles
ispin = 1 ! i==alpha ==> alpha/beta/alpha terms
do i = 1, nexc(ispin) ! number of couple of holes/particles
! exchange between (h1,p1) and i
ipart=occ_particle(i,ispin)
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)
hthree -= direct_int - exchange_int
enddo
ispin = 2 ! i==beta ==> alpha/beta/beta terms
do i = 1, nexc(ispin) ! number of couple of holes/particles
ispin = 2 ! i==beta ==> alpha/beta/beta terms
do i = 1, nexc(ispin) ! number of couple of holes/particles
! exchange between (h2,p2) and i
ipart=occ_particle(i,ispin)
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)]
implicit none
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
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
!$OMP PARALLEL &
!$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 DO SCHEDULE (static)
do hh1 = 1, n_act_orb !! alpha
h1 = list_act(hh1)
do hh2 = 1, n_act_orb !! beta
h2 = list_act(hh2)
!$OMP DO SCHEDULE (static)
do hh1 = 1, n_act_orb !! alpha
h1 = list_act(hh1)
do hh2 = 1, n_act_orb !! beta
h2 = list_act(hh2)
do pp1 = 1, n_act_orb !! alpha
p1 = list_act(pp1)
do pp2 = 1, n_act_orb !! beta
do pp2 = 1, n_act_orb !! beta
p2 = list_act(pp2)
call give_contrib_for_abab(h1,h2,p1,p2,occ,Ne,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 PARALLEL
END_PROVIDER
END_PROVIDER
subroutine give_contrib_for_abab(h1,h2,p1,p2,occ,Ne,contrib)
implicit none
BEGIN_DOC
BEGIN_DOC
! 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)
END_DOC
integer, intent(in) :: h1,h2,p1,p2,occ(N_int*bit_kind_size,2),Ne(2)
double precision, intent(out) :: contrib
integer :: mm,m
integer :: mm,m
double precision :: direct_int, exchange_int
!! h1,p1 == alpha
!! h1,p1 == alpha
!! h2,p2 == beta
contrib = 0.d0
do mm = 1, Ne(1) !! alpha
do mm = 1, Ne(1) !! alpha
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_int = three_e_5_idx_exch13_bi_ort(mm,p2,h2,p1,h1)
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
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_int = three_e_5_idx_exch23_bi_ort(mm,p2,h2,p1,h1)
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)]
implicit none
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
!
! 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 ||||
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
!$OMP PARALLEL &
!$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 DO SCHEDULE (static)
do hh1 = 1, n_act_orb !! alpha
h1 = list_act(hh1)
!$OMP DO SCHEDULE (static)
do hh1 = 1, n_act_orb !! alpha
h1 = list_act(hh1)
do hh2 = hh1+1, n_act_orb !! alpha
h2 = list_act(hh2)
h2 = list_act(hh2)
do pp1 = 1, n_act_orb !! alpha
p1 = list_act(pp1)
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 PARALLEL
END_PROVIDER
END_PROVIDER
subroutine give_contrib_for_aaaa(h1,h2,p1,p2,occ,Ne,contrib)
implicit none
BEGIN_DOC
BEGIN_DOC
! 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)
END_DOC
integer, intent(in) :: h1,h2,p1,p2,occ(N_int*bit_kind_size,2),Ne(2)
double precision, intent(out) :: contrib
integer :: mm,m
integer :: mm,m
double precision :: direct_int, exchange_int
!! h1,p1 == alpha
!! h1,p1 == alpha
!! h2,p2 == alpha
contrib = 0.d0
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
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_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
enddo
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)]
implicit none
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
!
! 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 ||||
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
!$OMP PARALLEL &
!$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 DO SCHEDULE (static)
do hh1 = 1, n_act_orb !! beta
h1 = list_act(hh1)
!$OMP DO SCHEDULE (static)
do hh1 = 1, n_act_orb !! beta
h1 = list_act(hh1)
do hh2 = hh1+1, n_act_orb !! beta
h2 = list_act(hh2)
h2 = list_act(hh2)
do pp1 = 1, n_act_orb !! beta
p1 = list_act(pp1)
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 PARALLEL
END_PROVIDER
END_PROVIDER
subroutine give_contrib_for_bbbb(h1,h2,p1,p2,occ,Ne,contrib)
implicit none
BEGIN_DOC
BEGIN_DOC
! 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)
END_DOC
integer, intent(in) :: h1,h2,p1,p2,occ(N_int*bit_kind_size,2),Ne(2)
double precision, intent(out) :: contrib
integer :: mm,m
integer :: mm,m
double precision :: direct_int, exchange_int
!! h1,p1 == 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
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_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
enddo
end
@ -424,17 +430,17 @@ end
subroutine double_htilde_mu_mat_fock_bi_ortho_no_3e(Nint, key_j, key_i, htot)
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 !!
!
!
! Non hermitian !!
END_DOC
use bitmasks
implicit none
integer, intent(in) :: Nint
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
@ -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)
if(s1.ne.s2)then
! opposite spin two-body
htwoe = mo_bi_ortho_tc_two_e(p2,p1,h2,h1)
! 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)
! 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
htot = htwoe
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_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_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