From 4cc8dae42010e062f82ace4373e2d5927e9074b0 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 2 Jun 2023 20:32:31 +0200 Subject: [PATCH] Improve 5idx --- scripts/qp_import_trexio.py | 15 ++ scripts/utility/qp_bitmasks.py | 6 +- src/ao_basis/cosgtos.irp.f | 3 +- src/bi_ort_ints/bi_ort_ints.irp.f | 115 ++++++------- src/bi_ort_ints/three_body_ijmkl.irp.f | 10 +- src/tc_bi_ortho/31.tc_bi_ortho.bats | 34 ++-- src/tc_bi_ortho/slater_tc_3e_slow.irp.f | 89 +++++----- src/tc_bi_ortho/slater_tc_opt.irp.f | 3 +- src/tc_bi_ortho/slater_tc_opt_double.irp.f | 180 +++++++++++---------- src/tc_bi_ortho/symmetrized_3_e_int.irp.f | 3 +- 10 files changed, 247 insertions(+), 211 deletions(-) diff --git a/scripts/qp_import_trexio.py b/scripts/qp_import_trexio.py index 89096387..2c829f5c 100755 --- a/scripts/qp_import_trexio.py +++ b/scripts/qp_import_trexio.py @@ -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") + diff --git a/scripts/utility/qp_bitmasks.py b/scripts/utility/qp_bitmasks.py index 38aa48d7..11965b72 100644 --- a/scripts/utility/qp_bitmasks.py +++ b/scripts/utility/qp_bitmasks.py @@ -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] diff --git a/src/ao_basis/cosgtos.irp.f b/src/ao_basis/cosgtos.irp.f index 721a3e57..dfa7d6b9 100644 --- a/src/ao_basis/cosgtos.irp.f +++ b/src/ao_basis/cosgtos.irp.f @@ -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 diff --git a/src/bi_ort_ints/bi_ort_ints.irp.f b/src/bi_ort_ints/bi_ort_ints.irp.f index f7a42f37..42bbe315 100644 --- a/src/bi_ort_ints/bi_ort_ints.irp.f +++ b/src/bi_ort_ints/bi_ort_ints.irp.f @@ -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 diff --git a/src/bi_ort_ints/three_body_ijmkl.irp.f b/src/bi_ort_ints/three_body_ijmkl.irp.f index bd669163..7b39235b 100644 --- a/src/bi_ort_ints/three_body_ijmkl.irp.f +++ b/src/bi_ort_ints/three_body_ijmkl.irp.f @@ -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) = ::: 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 diff --git a/src/tc_bi_ortho/31.tc_bi_ortho.bats b/src/tc_bi_ortho/31.tc_bi_ortho.bats index f5b9d8c0..93bed2ab 100644 --- a/src/tc_bi_ortho/31.tc_bi_ortho.bats +++ b/src/tc_bi_ortho/31.tc_bi_ortho.bats @@ -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 } diff --git a/src/tc_bi_ortho/slater_tc_3e_slow.irp.f b/src/tc_bi_ortho/slater_tc_3e_slow.irp.f index 6abb6b78..49977f37 100644 --- a/src/tc_bi_ortho/slater_tc_3e_slow.irp.f +++ b/src/tc_bi_ortho/slater_tc_3e_slow.irp.f @@ -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) ! 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) ! 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 diff --git a/src/tc_bi_ortho/slater_tc_opt.irp.f b/src/tc_bi_ortho/slater_tc_opt.irp.f index 3fd2576a..882470ed 100644 --- a/src/tc_bi_ortho/slater_tc_opt.irp.f +++ b/src/tc_bi_ortho/slater_tc_opt.irp.f @@ -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 diff --git a/src/tc_bi_ortho/slater_tc_opt_double.irp.f b/src/tc_bi_ortho/slater_tc_opt_double.irp.f index 2d6bfb27..12bbbec0 100644 --- a/src/tc_bi_ortho/slater_tc_opt_double.irp.f +++ b/src/tc_bi_ortho/slater_tc_opt_double.irp.f @@ -2,17 +2,17 @@ subroutine double_htilde_mu_mat_fock_bi_ortho(Nint, key_j, key_i, hmono, htwoe, hthree, htot) BEGIN_DOC - ! for double excitation ONLY FOR ONE- AND TWO-BODY TERMS + ! 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 - ! for double excitation ONLY FOR ONE- AND TWO-BODY TERMS + ! 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 diff --git a/src/tc_bi_ortho/symmetrized_3_e_int.irp.f b/src/tc_bi_ortho/symmetrized_3_e_int.irp.f index e4f7ca93..e725d8e5 100644 --- a/src/tc_bi_ortho/symmetrized_3_e_int.irp.f +++ b/src/tc_bi_ortho/symmetrized_3_e_int.irp.f @@ -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