PTEROSOR/QuantumPackage
10
0
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-07-03 09:56:04 +02:00
Code Releases Activity

Merge pull request #37 from QuantumPackage/dev-stable

Browse Source 
Dev stable
This commit is contained in:
AbdAmmar 2024-06-26 08:27:56 +02:00 committed by GitHub
commit 0c9245d1b3
No known key found for this signature in database
GPG Key ID: B5690EEEBB952194
116 changed files with 3773 additions and 3084 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

23
bin/zcat
View File

@ -1,23 +0,0 @@
#!/bin/bash
# On Darwin: try gzcat if available, otherwise use Python
if [[ $(uname -s) = Darwin ]] ; then
which gzcat &> /dev/null
if [[ $? -eq 0 ]] ; then
exec gzcat $@
else
exec python3 << EOF
import sys
import gzip
with gzip.open("$1", "rt") as f:
print(f.read())
EOF
fi
else
SCRIPTPATH="$( cd -- "$(dirname "$0")" >/dev/null 2>&1 ; pwd -P )"
command=$(which -a zcat | grep -v "$SCRIPTPATH/" | head -1)
exec $command $@
fi

63
config/gfortran_debug_mkl.cfg Normal file
View File

@ -0,0 +1,63 @@
# Common flags
##############
#
# -ffree-line-length-none : Needed for IRPF90 which produces long lines
# -lblas -llapack : Link with libblas and liblapack libraries provided by the system
# -I . : Include the curent directory (Mandatory)
#
# --ninja : Allow the utilisation of ninja. (Mandatory)
# --align=32 : Align all provided arrays on a 32-byte boundary
#
#
[COMMON]
FC : gfortran -g -ffree-line-length-none -I . -fPIC -std=legacy
LAPACK_LIB : -I${MKLROOT}/include -L${MKLROOT}/lib/intel64 -Wl,--no-as-needed -lmkl_gf_lp64 -lmkl_core -lpthread -lm -ldl -lmkl_gnu_thread -lgomp -fopenmp
IRPF90 : irpf90
IRPF90_FLAGS : --ninja --align=32 --assert -DSET_NESTED
# Global options
################
#
# 1 : Activate
# 0 : Deactivate
#
[OPTION]
MODE : DEBUG ; [ OPT | PROFILE | DEBUG ] : Chooses the section below
CACHE : 0 ; Enable cache_compile.py
OPENMP : 1 ; Append OpenMP flags
# Optimization flags
####################
#
# -Ofast : Disregard strict standards compliance. Enables all -O3 optimizations.
# It also enables optimizations that are not valid
# for all standard-compliant programs. It turns on
# -ffast-math and the Fortran-specific
# -fno-protect-parens and -fstack-arrays.
[OPT]
FCFLAGS : -Ofast
# Profiling flags
#################
#
[PROFILE]
FC : -p -g
FCFLAGS : -Ofast
# Debugging flags
#################
#
# -fcheck=all : Checks uninitialized variables, array subscripts, etc...
# -g : Extra debugging information
#
[DEBUG]
#FCFLAGS : -g -msse4.2 -fcheck=all -Waliasing -Wampersand -Wconversion -Wsurprising -Wintrinsics-std -Wno-tabs -Wintrinsic-shadow -Wline-truncation -Wreal-q-constant -Wuninitialized -fbacktrace -ffpe-trap=zero,overflow,underflow -finit-real=nan
FCFLAGS : -g -mavx -fcheck=all -Waliasing -Wampersand -Wconversion -Wsurprising -Wintrinsics-std -Wno-tabs -Wintrinsic-shadow -Wline-truncation -Wreal-q-constant -Wuninitialized -fbacktrace -ffpe-trap=zero,overflow -finit-real=nan
# OpenMP flags
#################
#
[OPENMP]
FC : -fopenmp
IRPF90_FLAGS : --openmp

9
etc/paths.rc
View File

@ -28,6 +28,15 @@ function qp_prepend_export () {
fi fi
} }
function qp_append_export () {
eval "value_1="\${$1}""
if [[ -z $value_1 ]] ; then
echo "${2}:"
else
echo "${value_1}:${2}"
fi
}
export PYTHONPATH=$(qp_prepend_export "PYTHONPATH" "${QP_EZFIO}/Python":"${QP_PYTHON}") export PYTHONPATH=$(qp_prepend_export "PYTHONPATH" "${QP_EZFIO}/Python":"${QP_PYTHON}")
export PATH=$(qp_prepend_export "PATH" "${QP_PYTHON}":"${QP_ROOT}"/bin:"${QP_ROOT}"/ocaml) export PATH=$(qp_prepend_export "PATH" "${QP_PYTHON}":"${QP_ROOT}"/bin:"${QP_ROOT}"/ocaml)

8
plugins/local/basis_correction/51.basis_c.bats
View File

@ -37,14 +37,6 @@ function run_sd() {
eq $energy1 $1 $thresh eq $energy1 $1 $thresh
} }
@test "O2 CAS" {
qp set_file o2_cas.gms.ezfio
qp set_mo_class -c "[1-2]" -a "[3-10]" -d "[11-46]"
run -149.72435425 3.e-4 10000
qp set_mo_class -c "[1-2]" -a "[3-10]" -v "[11-46]"
run_md -0.1160222327 1.e-6
}
@test "LiF RHF" { @test "LiF RHF" {
qp set_file lif.ezfio qp set_file lif.ezfio

4
plugins/local/basis_correction/basis_correction.irp.f
View File

@ -7,10 +7,6 @@ program basis_correction
touch read_wf touch read_wf
no_core_density = .True. no_core_density = .True.
touch no_core_density touch no_core_density
if(io_mo_two_e_integrals .ne. "Read")then
provide ao_two_e_integrals_in_map
endif
provide mo_two_e_integrals_in_map
call print_basis_correction call print_basis_correction
end end

2
plugins/local/basis_correction/print_routine.irp.f
View File

@ -22,7 +22,7 @@ subroutine print_basis_correction
print*, '****************************************' print*, '****************************************'
print*, '****************************************' print*, '****************************************'
print*, 'mu_of_r_potential = ',mu_of_r_potential print*, 'mu_of_r_potential = ',mu_of_r_potential
if(mu_of_r_potential.EQ."hf")then if(mu_of_r_potential.EQ."hf".or.mu_of_r_potential.EQ."hf_old".or.mu_of_r_potential.EQ."hf_sparse")then
print*, '' print*, ''
print*,'Using a HF-like two-body density to define mu(r)' print*,'Using a HF-like two-body density to define mu(r)'
print*,'This assumes that HF is a qualitative representation of the wave function ' print*,'This assumes that HF is a qualitative representation of the wave function '

18
plugins/local/basis_correction/test_chol_bas.irp.f Normal file
View File

@ -0,0 +1,18 @@
program pouet
implicit none
call test
end
subroutine test
implicit none
! provide mos_times_cholesky_r1
! provide mos_times_cholesky_r2
integer :: ipoint
double precision :: accu,weight
accu = 0.d0
do ipoint = 1, n_points_final_grid
weight = final_weight_at_r_vector(ipoint)
! accu += dabs(mu_of_r_hf(ipoint) - mu_of_r_hf_old(ipoint)) * weight
accu += dabs(f_hf_cholesky_sparse(ipoint) - f_hf_cholesky(ipoint)) * weight
enddo
print*,'accu = ',accu
end

28
plugins/local/bi_ort_ints/total_twoe_pot.irp.f
View File

@ -259,15 +259,21 @@ BEGIN_PROVIDER [ double precision, mo_bi_ortho_tc_two_e_transp, (mo_num, mo_num,
END_DOC END_DOC
integer :: i,j,k,l integer :: i,j,k,l
print*,'Providing mo_bi_ortho_tc_two_e_transp'
double precision :: t0,t1
call wall_time(t0)
do i = 1, mo_num do i = 1, mo_num
do j = 1, mo_num do j = 1, mo_num
do k = 1, mo_num do k = 1, mo_num
do l = 1, mo_num do l = 1, mo_num
mo_bi_ortho_tc_two_e_transp(i,j,k,l) = mo_bi_ortho_tc_two_e_transp(k,l,i,j) mo_bi_ortho_tc_two_e_transp(i,j,k,l) = mo_bi_ortho_tc_two_e(k,l,i,j)
enddo enddo
enddo enddo
enddo enddo
enddo enddo
call wall_time(t1)
print *, ' WALL TIME for PROVIDING mo_bi_ortho_tc_two_e_transp (min)', (t1-t0)/60.d0
END_PROVIDER END_PROVIDER
! --- ! ---
@ -326,3 +332,23 @@ END_PROVIDER
! --- ! ---
BEGIN_PROVIDER [double precision, tc_2e_3idx_coulomb_integrals_transp , (mo_num,mo_num,mo_num)]
&BEGIN_PROVIDER [double precision, tc_2e_3idx_exchange_integrals_transp, (mo_num,mo_num,mo_num)]
BEGIN_DOC
! tc_2e_3idx_coulomb_integrals_transp (j,k,i) = <jk|ji>
! tc_2e_3idx_exchange_integrals_transp(j,k,i) = <kj|ji>
END_DOC
implicit none
integer :: i, j, k
do i = 1, mo_num
do k = 1, mo_num
do j = 1, mo_num
tc_2e_3idx_coulomb_integrals_transp(j, k,i) = mo_bi_ortho_tc_two_e_transp(j ,k ,j ,i )
tc_2e_3idx_exchange_integrals_transp(j,k,i) = mo_bi_ortho_tc_two_e_transp(k ,j ,j ,i )
enddo
enddo
enddo
END_PROVIDER

108
plugins/local/cipsi_tc_bi_ortho/get_d0_transp.irp.f Normal file
View File

@ -0,0 +1,108 @@
subroutine get_d0_transp(gen, phasemask, bannedOrb, banned, mat_l, mat_r, mask, h, p, sp, coefs)
!todo: indices/conjg should be okay for complex
use bitmasks
implicit none
integer(bit_kind), intent(in) :: gen(N_int, 2), mask(N_int, 2)
integer(bit_kind), intent(in) :: phasemask(N_int,2)
logical, intent(in) :: bannedOrb(mo_num, 2), banned(mo_num, mo_num,2)
integer(bit_kind) :: det(N_int, 2)
double precision, intent(in) :: coefs(N_states,2)
double precision, intent(inout) :: mat_l(N_states, mo_num, mo_num)
double precision, intent(inout) :: mat_r(N_states, mo_num, mo_num)
integer, intent(in) :: h(0:2,2), p(0:4,2), sp
integer :: i, j, k, s, h1, h2, p1, p2, puti, putj, mm
double precision :: phase
double precision :: hij,hji
double precision, external :: get_phase_bi
logical :: ok
integer, parameter :: bant=1
double precision, allocatable :: hij_cache1(:), hij_cache2(:)
allocate (hij_cache1(mo_num),hij_cache2(mo_num))
double precision, allocatable :: hji_cache1(:), hji_cache2(:)
allocate (hji_cache1(mo_num),hji_cache2(mo_num))
! print*,'in get_d0_new'
! call debug_det(gen,N_int)
! print*,'coefs',coefs(1,:)
if(sp == 3) then ! AB
h1 = p(1,1)
h2 = p(1,2)
do p1=1, mo_num
if(bannedOrb(p1, 1)) cycle
! call get_mo_two_e_integrals_complex(p1,h2,h1,mo_num,hij_cache1,mo_integrals_map)
do mm = 1, mo_num
hij_cache1(mm) = mo_bi_ortho_tc_two_e(mm,p1,h2,h1)
hji_cache1(mm) = mo_bi_ortho_tc_two_e_transp(mm,p1,h2,h1)
enddo
!!!!!!!!!! <alpha|H|psi>
do p2=1, mo_num
if(bannedOrb(p2,2)) cycle
if(banned(p1, p2, bant)) cycle ! rentable?
if(p1 == h1 .or. p2 == h2) then
call apply_particles(mask, 1,p1,2,p2, det, ok, N_int)
! call i_h_j_complex(gen, det, N_int, hij) ! need to take conjugate of this
! call i_h_j_complex(det, gen, N_int, hij)
call htilde_mu_mat_opt_bi_ortho_no_3e_both(det,gen,N_int, hij,hji)
else
phase = get_phase_bi(phasemask, 1, 2, h1, p1, h2, p2, N_int)
hij = hij_cache1(p2) * phase
hji = hji_cache1(p2) * phase
end if
if (hij == 0.d0.or.hji == 0.d0) cycle
!DIR$ LOOP COUNT AVG(4)
do k=1,N_states
mat_r(k, p1, p2) = mat_r(k, p1, p2) + coefs(k,2) * hij ! HOTSPOT
mat_l(k, p1, p2) = mat_l(k, p1, p2) + coefs(k,1) * hji ! HOTSPOT
enddo
end do
end do
else ! AA BB
p1 = p(1,sp)
p2 = p(2,sp)
do puti=1, mo_num
if(bannedOrb(puti, sp)) cycle
! call get_mo_two_e_integrals_complex(puti,p2,p1,mo_num,hij_cache1,mo_integrals_map,mo_integrals_map_2)
! call get_mo_two_e_integrals_complex(puti,p1,p2,mo_num,hij_cache2,mo_integrals_map,mo_integrals_map_2)
do mm = 1, mo_num
hij_cache1(mm) = mo_bi_ortho_tc_two_e(mm,puti,p2,p1)
hij_cache2(mm) = mo_bi_ortho_tc_two_e(mm,puti,p1,p2)
hji_cache1(mm) = mo_bi_ortho_tc_two_e_transp(mm,puti,p2,p1)
hji_cache2(mm) = mo_bi_ortho_tc_two_e_transp(mm,puti,p1,p2)
enddo
!!!!!!!!!! <alpha|H|psi>
do putj=puti+1, mo_num
if(bannedOrb(putj, sp)) cycle
if(banned(puti, putj, bant)) cycle ! rentable?
if(puti == p1 .or. putj == p2 .or. puti == p2 .or. putj == p1) then
call apply_particles(mask, sp,puti,sp,putj, det, ok, N_int)
!call i_h_j_complex(gen, det, N_int, hij) ! need to take conjugate of this
! call i_h_j_complex(det, gen, N_int, hij)
call htilde_mu_mat_opt_bi_ortho_no_3e_both(det,gen,N_int, hij,hji)
if (hij == 0.d0.or.hji == 0.d0) cycle
else
! hij = (mo_two_e_integral_complex(p1, p2, puti, putj) - mo_two_e_integral_complex(p2, p1, puti, putj))
! hij = (mo_bi_ortho_tc_two_e(p1, p2, puti, putj) - mo_bi_ortho_tc_two_e(p2, p1, puti, putj))
hij = (mo_bi_ortho_tc_two_e(puti, putj, p1, p2) - mo_bi_ortho_tc_two_e(puti, putj, p2, p1))
hji = (mo_bi_ortho_tc_two_e_transp(puti, putj, p1, p2) - mo_bi_ortho_tc_two_e_transp(puti, putj, p2, p1))
if (hij == 0.d0.or.hji == 0.d0) cycle
phase = get_phase_bi(phasemask, sp, sp, puti, p1 , putj, p2, N_int)
hij = (hij) * phase
hji = (hji) * phase
end if
!DIR$ LOOP COUNT AVG(4)
do k=1,N_states
mat_r(k, puti, putj) = mat_r(k, puti, putj) + coefs(k,2) * hij
mat_l(k, puti, putj) = mat_l(k, puti, putj) + coefs(k,1) * hji
enddo
end do
end do
end if
deallocate(hij_cache1,hij_cache2)
end

358
plugins/local/cipsi_tc_bi_ortho/get_d1_transp.irp.f Normal file
View File

@ -0,0 +1,358 @@
subroutine get_d1_transp(gen, phasemask, bannedOrb, banned, mat_l, mat_r, mask, h, p, sp, coefs)
!todo: indices should be okay for complex?
use bitmasks
implicit none
integer(bit_kind), intent(in) :: mask(N_int, 2), gen(N_int, 2)
integer(bit_kind), intent(in) :: phasemask(N_int,2)
logical, intent(in) :: bannedOrb(mo_num, 2), banned(mo_num, mo_num,2)
integer(bit_kind) :: det(N_int, 2)
double precision, intent(in) :: coefs(N_states,2)
double precision, intent(inout) :: mat_l(N_states, mo_num, mo_num)
double precision, intent(inout) :: mat_r(N_states, mo_num, mo_num)
integer, intent(in) :: h(0:2,2), p(0:4,2), sp
double precision, external :: get_phase_bi
double precision, external :: mo_two_e_integral_complex
logical :: ok
logical, allocatable :: lbanned(:,:)
integer :: puti, putj, ma, mi, s1, s2, i, i1, i2, j, istate
integer :: hfix, pfix, h1, h2, p1, p2, ib, k, l, mm
integer, parameter :: turn2(2) = (/2,1/)
integer, parameter :: turn3(2,3) = reshape((/2,3, 1,3, 1,2/), (/2,3/))
integer :: bant
double precision, allocatable :: hij_cache(:,:)
double precision :: hij, tmp_rowij(N_states, mo_num), tmp_rowij2(N_states, mo_num),phase
double precision, allocatable :: hji_cache(:,:)
double precision :: hji, tmp_rowji(N_states, mo_num), tmp_rowji2(N_states, mo_num)
! PROVIDE mo_integrals_map N_int
! print*,'in get_d1_new'
! call debug_det(gen,N_int)
! print*,'coefs',coefs(1,:)
allocate (lbanned(mo_num, 2))
allocate (hij_cache(mo_num,2))
allocate (hji_cache(mo_num,2))
lbanned = bannedOrb
do i=1, p(0,1)
lbanned(p(i,1), 1) = .true.
end do
do i=1, p(0,2)
lbanned(p(i,2), 2) = .true.
end do
ma = 1
if(p(0,2) >= 2) ma = 2
mi = turn2(ma)
bant = 1
if(sp == 3) then
!move MA
if(ma == 2) bant = 2
puti = p(1,mi)
hfix = h(1,ma)
p1 = p(1,ma)
p2 = p(2,ma)
if(.not. bannedOrb(puti, mi)) then
! call get_mo_two_e_integrals_complex(hfix,p1,p2,mo_num,hij_cache(1,1),mo_integrals_map,mo_integrals_map_2)
! call get_mo_two_e_integrals_complex(hfix,p2,p1,mo_num,hij_cache(1,2),mo_integrals_map,mo_integrals_map_2)
do mm = 1, mo_num
hij_cache(mm,1) = mo_bi_ortho_tc_two_e(mm,hfix,p1,p2)
hij_cache(mm,2) = mo_bi_ortho_tc_two_e(mm,hfix,p2,p1)
hji_cache(mm,1) = mo_bi_ortho_tc_two_e_transp(mm,hfix,p1,p2)
hji_cache(mm,2) = mo_bi_ortho_tc_two_e_transp(mm,hfix,p2,p1)
do istate = 1,N_states
tmp_rowij(istate,mm) = 0.d0
tmp_rowji(istate,mm) = 0.d0
enddo
enddo
!! <alpha|H|psi>
do putj=1, hfix-1
if(lbanned(putj, ma)) cycle
if(banned(putj, puti,bant)) cycle
hij = hij_cache(putj,1) - hij_cache(putj,2)
hji = hji_cache(putj,1) - hji_cache(putj,2)
if (hij /= 0.d0.and.hji/=0.d0) then
phase = get_phase_bi(phasemask, ma, ma, putj, p1, hfix, p2, N_int)
hij = hij * phase
hji = hji * phase
!DIR$ LOOP COUNT AVG(4)
do k=1,N_states
tmp_rowij(k,putj) = tmp_rowij(k,putj) + hij * coefs(k,2)
tmp_rowji(k,putj) = tmp_rowji(k,putj) + hji * coefs(k,1)
enddo
endif
end do
do putj=hfix+1, mo_num
if(lbanned(putj, ma)) cycle
if(banned(putj, puti,bant)) cycle
hij = hij_cache(putj,2) - hij_cache(putj,1)
hji = hji_cache(putj,2) - hji_cache(putj,1)
if (hij /= 0.d0.and.hji/=0.d0) then
phase = get_phase_bi(phasemask, ma, ma, hfix, p1, putj, p2, N_int)
hij = hij * phase
hji = hji * phase
!DIR$ LOOP COUNT AVG(4)
do k=1,N_states
tmp_rowij(k,putj) = tmp_rowij(k,putj) + hij * coefs(k,2)
tmp_rowji(k,putj) = tmp_rowji(k,putj) + hji * coefs(k,1)
enddo
endif
end do
if(ma == 1) then
mat_r(1:N_states,1:mo_num,puti) = mat_r(1:N_states,1:mo_num,puti) + tmp_rowij(1:N_states,1:mo_num)
mat_l(1:N_states,1:mo_num,puti) = mat_l(1:N_states,1:mo_num,puti) + tmp_rowji(1:N_states,1:mo_num)
else
do l=1,mo_num
!DIR$ LOOP COUNT AVG(4)
do k=1,N_states
mat_r(k,puti,l) = mat_r(k,puti,l) + tmp_rowij(k,l)
mat_l(k,puti,l) = mat_l(k,puti,l) + tmp_rowji(k,l)
enddo
enddo
end if
end if
!MOVE MI
pfix = p(1,mi)
! call get_mo_two_e_integrals_complex(hfix,pfix,p1,mo_num,hij_cache(1,1),mo_integrals_map,mo_integrals_map_2)
! call get_mo_two_e_integrals_complex(hfix,pfix,p2,mo_num,hij_cache(1,2),mo_integrals_map,mo_integrals_map_2)
do mm = 1, mo_num
do istate = 1,N_states
tmp_rowij(istate,mm) = 0.d0
tmp_rowij2(istate,mm) = 0.d0
tmp_rowji(istate,mm) = 0.d0
tmp_rowji2(istate,mm) = 0.d0
enddo
hij_cache(mm,1) = mo_bi_ortho_tc_two_e(mm,hfix,pfix,p1)
hij_cache(mm,2) = mo_bi_ortho_tc_two_e(mm,hfix,pfix,p2)
hji_cache(mm,1) = mo_bi_ortho_tc_two_e_transp(mm,hfix,pfix,p1)
hji_cache(mm,2) = mo_bi_ortho_tc_two_e_transp(mm,hfix,pfix,p2)
enddo
putj = p1
!! <alpha|H|psi>
do puti=1,mo_num !HOT
if(lbanned(puti,mi)) cycle
!p1 fixed
putj = p1
if(.not. banned(putj,puti,bant)) then
hij = hij_cache(puti,2)
hji = hji_cache(puti,2)
if (hij /= 0.d0.and.hji/=0.d0) then
phase = get_phase_bi(phasemask, ma, mi, hfix, p2, puti, pfix, N_int)
hij = hij * phase
hji = hji * phase
!DIR$ LOOP COUNT AVG(4)
do k=1,N_states
tmp_rowij(k,puti) = tmp_rowij(k,puti) + hij * coefs(k,2)
tmp_rowji(k,puti) = tmp_rowji(k,puti) + hji * coefs(k,1)
enddo
endif
end if
!
putj = p2
if(.not. banned(putj,puti,bant)) then
hij = hij_cache(puti,1)
hji = hji_cache(puti,1)
if (hij /= 0.d0.and.hji/=0.d0) then
phase = get_phase_bi(phasemask, ma, mi, hfix, p1, puti, pfix, N_int)
hij = hij * phase
hji = hji * phase
do k=1,N_states
tmp_rowij2(k,puti) = tmp_rowij2(k,puti) + hij * coefs(k,2)
tmp_rowji2(k,puti) = tmp_rowji2(k,puti) + hji * coefs(k,1)
enddo
endif
end if
end do
if(mi == 1) then
mat_r(:,:,p1) = mat_r(:,:,p1) + tmp_rowij(:,:)
mat_r(:,:,p2) = mat_r(:,:,p2) + tmp_rowij2(:,:)
mat_l(:,:,p1) = mat_l(:,:,p1) + tmp_rowji(:,:)
mat_l(:,:,p2) = mat_l(:,:,p2) + tmp_rowji2(:,:)
else
do l=1,mo_num
!DIR$ LOOP COUNT AVG(4)
do k=1,N_states
mat_r(k,p1,l) = mat_r(k,p1,l) + tmp_rowij(k,l)
mat_r(k,p2,l) = mat_r(k,p2,l) + tmp_rowij2(k,l)
mat_l(k,p1,l) = mat_l(k,p1,l) + tmp_rowji(k,l)
mat_l(k,p2,l) = mat_l(k,p2,l) + tmp_rowji2(k,l)
enddo
enddo
end if
else ! sp /= 3
if(p(0,ma) == 3) then
do i=1,3
hfix = h(1,ma)
puti = p(i, ma)
p1 = p(turn3(1,i), ma)
p2 = p(turn3(2,i), ma)
! call get_mo_two_e_integrals_complex(hfix,p1,p2,mo_num,hij_cache(1,1),mo_integrals_map,mo_integrals_map_2)
! call get_mo_two_e_integrals_complex(hfix,p2,p1,mo_num,hij_cache(1,2),mo_integrals_map,mo_integrals_map_2)
do mm = 1, mo_num
hij_cache(mm,1) = mo_bi_ortho_tc_two_e(mm,hfix,p1,p2)
hij_cache(mm,2) = mo_bi_ortho_tc_two_e(mm,hfix,p2,p1)
hji_cache(mm,1) = mo_bi_ortho_tc_two_e_transp(mm,hfix,p1,p2)
hji_cache(mm,2) = mo_bi_ortho_tc_two_e_transp(mm,hfix,p2,p1)
do istate = 1, N_states
tmp_rowij(istate,mm) = 0.d0
tmp_rowji(istate,mm) = 0.d0
enddo
enddo
!! <alpha|H|psi>
do putj=1,hfix-1
if(banned(putj,puti,1)) cycle
if(lbanned(putj,ma)) cycle
hij = hij_cache(putj,1) - hij_cache(putj,2)
hji = hji_cache(putj,1) - hji_cache(putj,2)
if (hij /= 0.d0.and.hji/=0.d0) then
phase = get_phase_bi(phasemask, ma, ma, putj, p1, hfix, p2, N_int)
hij = hij * phase
hji = hji * phase
tmp_rowij(:,putj) = tmp_rowij(:,putj) + hij * coefs(:,2)
tmp_rowji(:,putj) = tmp_rowji(:,putj) + hji * coefs(:,1)
endif
end do
do putj=hfix+1,mo_num
if(banned(putj,puti,1)) cycle
if(lbanned(putj,ma)) cycle
hij = hij_cache(putj,2) - hij_cache(putj,1)
hji = hji_cache(putj,2) - hji_cache(putj,1)
if (hij /= 0.d0.and.hji/=0.d0) then
phase = get_phase_bi(phasemask, ma, ma, hfix, p1, putj, p2, N_int)
hij = hij * phase
hji = hji * phase
tmp_rowij(:,putj) = tmp_rowij(:,putj) + hij * coefs(:,2)
tmp_rowji(:,putj) = tmp_rowji(:,putj) + hji * coefs(:,1)
endif
end do
mat_r(:, :puti-1, puti) = mat_r(:, :puti-1, puti) + tmp_rowij(:,:puti-1)
mat_l(:, :puti-1, puti) = mat_l(:, :puti-1, puti) + tmp_rowji(:,:puti-1)
do l=puti,mo_num
!DIR$ LOOP COUNT AVG(4)
do k=1,N_states
mat_r(k, puti, l) = mat_r(k, puti,l) + tmp_rowij(k,l)
mat_l(k, puti, l) = mat_l(k, puti,l) + tmp_rowji(k,l)
enddo
enddo
end do
else
hfix = h(1,mi)
pfix = p(1,mi)
p1 = p(1,ma)
p2 = p(2,ma)
! call get_mo_two_e_integrals_complex(hfix,p1,pfix,mo_num,hij_cache(1,1),mo_integrals_map,mo_integrals_map_2)
! call get_mo_two_e_integrals_complex(hfix,p2,pfix,mo_num,hij_cache(1,2),mo_integrals_map,mo_integrals_map_2)
do mm = 1, mo_num
hij_cache(mm,1) = mo_bi_ortho_tc_two_e(mm,hfix,p1,pfix)
hij_cache(mm,2) = mo_bi_ortho_tc_two_e(mm,hfix,p2,pfix)
hji_cache(mm,1) = mo_bi_ortho_tc_two_e_transp(mm,hfix,p1,pfix)
hji_cache(mm,2) = mo_bi_ortho_tc_two_e_transp(mm,hfix,p2,pfix)
do istate = 1,N_states
tmp_rowij (istate,mm) = 0.d0
tmp_rowij2(istate,mm) = 0.d0
tmp_rowji (istate,mm) = 0.d0
tmp_rowji2(istate,mm) = 0.d0
enddo
enddo
putj = p2
!! <alpha|H|psi>
do puti=1,mo_num
if(lbanned(puti,ma)) cycle
putj = p2
if(.not. banned(puti,putj,1)) then
hij = hij_cache(puti,1)
hji = hji_cache(puti,1)
if (hij /= 0.d0.and.hji/=0.d0) then
phase = get_phase_bi(phasemask, mi, ma, hfix, pfix, puti, p1, N_int)
hij = hij * phase
hji = hji * phase
!DIR$ LOOP COUNT AVG(4)
do k=1,N_states
tmp_rowij(k,puti) = tmp_rowij(k,puti) + hij * coefs(k,2)
tmp_rowji(k,puti) = tmp_rowji(k,puti) + hji * coefs(k,1)
enddo
endif
end if
putj = p1
if(.not. banned(puti,putj,1)) then
hij = hij_cache(puti,2)
hji = hji_cache(puti,2)
if (hij /= 0.d0.and.hji/=0.d0) then
phase = get_phase_bi(phasemask, mi, ma, hfix, pfix, puti, p2, N_int)
hij = hij * phase
hji = hji * phase
do k=1,N_states
tmp_rowij2(k,puti) = tmp_rowij2(k,puti) + hij * coefs(k,2)
tmp_rowji2(k,puti) = tmp_rowji2(k,puti) + hji * coefs(k,1)
enddo
endif
end if
end do
mat_r(:,:p2-1,p2) = mat_r(:,:p2-1,p2) + tmp_rowij(:,:p2-1)
mat_l(:,:p2-1,p2) = mat_l(:,:p2-1,p2) + tmp_rowji(:,:p2-1)
do l=p2,mo_num
!DIR$ LOOP COUNT AVG(4)
do k=1,N_states
mat_r(k,p2,l) = mat_r(k,p2,l) + tmp_rowij(k,l)
mat_l(k,p2,l) = mat_l(k,p2,l) + tmp_rowji(k,l)
enddo
enddo
mat_r(:,:p1-1,p1) = mat_r(:,:p1-1,p1) + tmp_rowij2(:,:p1-1)
mat_l(:,:p1-1,p1) = mat_l(:,:p1-1,p1) + tmp_rowji2(:,:p1-1)
do l=p1,mo_num
!DIR$ LOOP COUNT AVG(4)
do k=1,N_states
mat_r(k,p1,l) = mat_r(k,p1,l) + tmp_rowij2(k,l)
mat_l(k,p1,l) = mat_l(k,p1,l) + tmp_rowji2(k,l)
enddo
enddo
end if
end if
deallocate(lbanned,hij_cache, hji_cache)
!! MONO
if(sp == 3) then
s1 = 1
s2 = 2
else
s1 = sp
s2 = sp
end if
do i1=1,p(0,s1)
ib = 1
if(s1 == s2) ib = i1+1
do i2=ib,p(0,s2)
p1 = p(i1,s1)
p2 = p(i2,s2)
if(bannedOrb(p1, s1) .or. bannedOrb(p2, s2) .or. banned(p1, p2, 1)) cycle
call apply_particles(mask, s1, p1, s2, p2, det, ok, N_int)
! gen is a selector; mask is ionized generator; det is alpha
! hij is contribution to <psi|H|alpha>
! call i_h_j_complex(gen, det, N_int, hij)
call htilde_mu_mat_opt_bi_ortho_no_3e_both(det, gen, N_int, hij,hji)
! call htilde_mu_mat_opt_bi_ortho_no_3e(gen, det, N_int, hji)
!DIR$ LOOP COUNT AVG(4)
do k=1,N_states
! take conjugate to get contribution to <alpha|H|psi> instead of <psi|H|alpha>
! mat_r(k, p1, p2) = mat_r(k, p1, p2) + coefs(k,1) * dconjg(hij)
mat_r(k, p1, p2) = mat_r(k, p1, p2) + coefs(k,2) * hij
mat_l(k, p1, p2) = mat_l(k, p1, p2) + coefs(k,1) * hji
enddo
end do
end do
end

3
plugins/local/cipsi_tc_bi_ortho/get_d2_good.irp.f
View File

@ -25,9 +25,6 @@ subroutine get_d2_new(gen, phasemask, bannedOrb, banned, mat_l, mat_r, mask, h,
integer :: bant integer :: bant
bant = 1 bant = 1
! print*, 'in get_d2_new'
! call debug_det(gen,N_int)
! print*,'coefs',coefs(1,:)
tip = p(0,1) * p(0,2) ! number of alpha particles times number of beta particles tip = p(0,1) * p(0,2) ! number of alpha particles times number of beta particles

235
plugins/local/cipsi_tc_bi_ortho/get_d2_transp.irp.f Normal file
View File

@ -0,0 +1,235 @@
subroutine get_d2_new_transp(gen, phasemask, bannedOrb, banned, mat_l, mat_r, mask, h, p, sp, coefs)
!todo: indices/conjg should be correct for complex
use bitmasks
implicit none
integer(bit_kind), intent(in) :: mask(N_int, 2), gen(N_int, 2)
integer(bit_kind), intent(in) :: phasemask(N_int,2)
logical, intent(in) :: bannedOrb(mo_num, 2), banned(mo_num, mo_num,2)
double precision, intent(in) :: coefs(N_states,2)
double precision, intent(inout) :: mat_r(N_states, mo_num, mo_num)
double precision, intent(inout) :: mat_l(N_states, mo_num, mo_num)
integer, intent(in) :: h(0:2,2), p(0:4,2), sp
double precision, external :: get_phase_bi
integer :: i, j, k, tip, ma, mi, puti, putj
integer :: h1, h2, p1, p2, i1, i2
double precision :: phase
double precision :: hij,hji
integer, parameter:: turn2d(2,3,4) = reshape((/0,0, 0,0, 0,0, 3,4, 0,0, 0,0, 2,4, 1,4, 0,0, 2,3, 1,3, 1,2 /), (/2,3,4/))
integer, parameter :: turn2(2) = (/2, 1/)
integer, parameter :: turn3(2,3) = reshape((/2,3, 1,3, 1,2/), (/2,3/))
integer :: bant
bant = 1
tip = p(0,1) * p(0,2) ! number of alpha particles times number of beta particles
ma = sp !1:(alpha,alpha); 2:(b,b); 3:(a,b)
if(p(0,1) > p(0,2)) ma = 1 ! more alpha particles than beta particles
if(p(0,1) < p(0,2)) ma = 2 ! fewer alpha particles than beta particles
mi = mod(ma, 2) + 1
if(sp == 3) then ! if one alpha and one beta xhole
!(where xholes refer to the ionizations from the generator, not the holes occupied in the ionized generator)
if(ma == 2) bant = 2 ! if more beta particles than alpha particles
if(tip == 3) then ! if 3 of one particle spin and 1 of the other particle spin
puti = p(1, mi)
if(bannedOrb(puti, mi)) return
h1 = h(1, ma)
h2 = h(2, ma)
!! <alpha|H|psi>
do i = 1, 3 ! loop over all 3 combinations of 2 particles with spin ma
putj = p(i, ma)
if(banned(putj,puti,bant)) cycle
i1 = turn3(1,i)
i2 = turn3(2,i)
p1 = p(i1, ma)
p2 = p(i2, ma)
! |G> = |psi_{gen,i}>
! |G'> = a_{x1} a_{x2} |G>
! |alpha> = a_{puti}^{\dagger} a_{putj}^{\dagger} |G'>
! |alpha> = t_{x1,x2}^{puti,putj} |G>
! hij = <psi_{selectors,i}|H|alpha>
! |alpha> = t_{p1,p2}^{h1,h2}|psi_{selectors,i}>
!todo: <i|H|j> = (<h1,h2|p1,p2> - <h1,h2|p2,p1>) * phase
! <psi|H|j> += dconjg(c_i) * <i|H|j>
! <j|H|i> = (<p1,p2|h1,h2> - <p2,p1|h1,h2>) * phase
! <j|H|psi> += <j|H|i> * c_i
!!!!!!!!!!!!! WARNING !!!!!!!!!!!!!!!!
! take the transpose of what's written above because later use the complex conjugate
! hij = mo_bi_ortho_tc_two_e(h1, h2, p1, p2) - mo_bi_ortho_tc_two_e( h1, h2, p2, p1)
! hji = mo_bi_ortho_tc_two_e_transp(h1, h2, p1, p2) - mo_bi_ortho_tc_two_e_transp( h1, h2, p2, p1)
hij = mo_bi_ortho_tc_two_e_transp(p1, p2,h1, h2) - mo_bi_ortho_tc_two_e_transp( p1, p2, h2, h1)
hji = mo_bi_ortho_tc_two_e(p1, p2, h1, h2) - mo_bi_ortho_tc_two_e( p1, p2, h2, h1)
if (hij == 0.d0.or.hji==0.d0) cycle
! take conjugate to get contribution to <alpha|H|psi> instead of <psi|H|alpha>
! hij = dconjg(hij) * get_phase_bi(phasemask, ma, ma, h1, p1, h2, p2, N_int)
phase = get_phase_bi(phasemask, ma, ma, h1, p1, h2, p2, N_int)
hij = hij * phase
hji = hji * phase
if(ma == 1) then ! if particle spins are (alpha,alpha,alpha,beta), then puti is beta and putj is alpha
!DIR$ LOOP COUNT AVG(4)
do k=1,N_states
mat_r(k, putj, puti) = mat_r(k, putj, puti) + coefs(k,2) * hij
mat_l(k, putj, puti) = mat_l(k, putj, puti) + coefs(k,1) * hji
enddo
else ! if particle spins are (beta,beta,beta,alpha), then puti is alpha and putj is beta
!DIR$ LOOP COUNT AVG(4)
do k=1,N_states
mat_r(k, puti, putj) = mat_r(k, puti, putj) + coefs(k,2) * hij
mat_l(k, puti, putj) = mat_l(k, puti, putj) + coefs(k,1) * hji
enddo
end if
end do
else ! if 2 alpha and 2 beta particles
h1 = h(1,1)
h2 = h(1,2)
!! <alpha|H|psi>
do j = 1,2 ! loop over all 4 combinations of one alpha and one beta particle
putj = p(j, 2)
if(bannedOrb(putj, 2)) cycle
p2 = p(turn2(j), 2)
do i = 1,2
puti = p(i, 1)
if(banned(puti,putj,bant) .or. bannedOrb(puti,1)) cycle
p1 = p(turn2(i), 1)
! hij = <psi_{selectors,i}|H|alpha>
! hij = mo_bi_ortho_tc_two_e(p1, p2, h1, h2)
!!!!!!!!!!!!! WARNING !!!!!!!!!!!!!!!!
! take the transpose of what's written above because later use the complex conjugate
! hij = mo_bi_ortho_tc_two_e(h1, h2, p1, p2 )
! hji = mo_bi_ortho_tc_two_e_transp(h1, h2, p1, p2 )
hij = mo_bi_ortho_tc_two_e_transp(p1, p2 ,h1, h2 )
hji = mo_bi_ortho_tc_two_e( p1, p2, h1, h2)
if (hij /= 0.d0.or.hji==0.d0) then
! take conjugate to get contribution to <alpha|H|psi> instead of <psi|H|alpha>
! hij = dconjg(hij) * get_phase_bi(phasemask, 1, 2, h1, p1, h2, p2, N_int)
phase = get_phase_bi(phasemask, 1, 2, h1, p1, h2, p2, N_int)
hij = hij * phase
hji = hji * phase
!DIR$ LOOP COUNT AVG(4)
do k=1,N_states
mat_r(k, puti, putj) = mat_r(k, puti, putj) + coefs(k,2) * hij
mat_l(k, puti, putj) = mat_l(k, puti, putj) + coefs(k,1) * hji
enddo
endif
end do
end do
end if
else ! if holes are (a,a) or (b,b)
if(tip == 0) then ! if particles are (a,a,a,a) or (b,b,b,b)
h1 = h(1, ma)
h2 = h(2, ma)
!! <alpha|H|psi>
do i=1,3
puti = p(i, ma)
if(bannedOrb(puti,ma)) cycle
do j=i+1,4
putj = p(j, ma)
if(bannedOrb(putj,ma)) cycle
if(banned(puti,putj,1)) cycle
i1 = turn2d(1, i, j)
i2 = turn2d(2, i, j)
p1 = p(i1, ma)
p2 = p(i2, ma)
! hij = mo_bi_ortho_tc_two_e(p1, p2, h1, h2) - mo_bi_ortho_tc_two_e(p2,p1, h1, h2)
!!!!!!!!!!!!! WARNING !!!!!!!!!!!!!!!!
! take the transpose of what's written above because later use the complex conjugate
hij = mo_bi_ortho_tc_two_e_transp(p1, p2, h1, h2) - mo_bi_ortho_tc_two_e_transp(p1, p2, h2,h1 )
hji = mo_bi_ortho_tc_two_e(p1, p2, h1, h2) - mo_bi_ortho_tc_two_e(p1, p2, h2,h1 )
if (hij == 0.d0.or.hji == 0.d0) cycle
! take conjugate to get contribution to <alpha|H|psi> instead of <psi|H|alpha>
! hij = dconjg(hij) * get_phase_bi(phasemask, ma, ma, h1, p1, h2, p2, N_int)
phase = get_phase_bi(phasemask, ma, ma, h1, p1, h2, p2, N_int)
hij = hij * phase
hji = hji * phase
!DIR$ LOOP COUNT AVG(4)
do k=1,N_states
mat_r(k, puti, putj) = mat_r(k, puti, putj) +coefs(k,2) * hij
mat_l(k, puti, putj) = mat_l(k, puti, putj) +coefs(k,1) * hji
enddo
end do
end do
else if(tip == 3) then ! if particles are (a,a,a,b) (ma=1,mi=2) or (a,b,b,b) (ma=2,mi=1)
h1 = h(1, mi)
h2 = h(1, ma)
p1 = p(1, mi)
!! <alpha|H|psi>
do i=1,3
puti = p(turn3(1,i), ma)
if(bannedOrb(puti,ma)) cycle
putj = p(turn3(2,i), ma)
if(bannedOrb(putj,ma)) cycle
if(banned(puti,putj,1)) cycle
p2 = p(i, ma)
! hij = mo_bi_ortho_tc_two_e(p1, p2, h1, h2)
!!!!!!!!!!!!! WARNING !!!!!!!!!!!!!!!!
! take the transpose of what's written above because later use the complex conjugate
hij = mo_bi_ortho_tc_two_e_transp(p1, p2 ,h1, h2)
hji = mo_bi_ortho_tc_two_e(p1, p2,h1, h2 )
if (hij == 0.d0) cycle
! take conjugate to get contribution to <alpha|H|psi> instead of <psi|H|alpha>
! hij = dconjg(hij) * get_phase_bi(phasemask, mi, ma, h1, p1, h2, p2, N_int)
phase = get_phase_bi(phasemask, mi, ma, h1, p1, h2, p2, N_int)
hij = hij * phase
hji = hji * phase
if (puti < putj) then
!DIR$ LOOP COUNT AVG(4)
do k=1,N_states
mat_r(k, puti, putj) = mat_r(k, puti, putj) + coefs(k,2) * hij
mat_l(k, puti, putj) = mat_l(k, puti, putj) + coefs(k,1) * hji
enddo
else
!DIR$ LOOP COUNT AVG(4)
do k=1,N_states
mat_r(k, putj, puti) = mat_r(k, putj, puti) + coefs(k,2) * hij
mat_l(k, putj, puti) = mat_l(k, putj, puti) + coefs(k,1) * hji
enddo
endif
end do
else ! tip == 4 (a,a,b,b)
puti = p(1, sp)
putj = p(2, sp)
if(.not. banned(puti,putj,1)) then
p1 = p(1, mi)
p2 = p(2, mi)
h1 = h(1, mi)
h2 = h(2, mi)
!! <alpha|H|psi>
! hij = (mo_bi_ortho_tc_two_e(p1, p2, h1, h2) - mo_bi_ortho_tc_two_e(p2,p1, h1, h2))
!!!!!!!!!!!!! WARNING !!!!!!!!!!!!!!!!
! take the transpose of what's written above because later use the complex conjugate
hij = (mo_bi_ortho_tc_two_e_transp(p1, p2,h1, h2) - mo_bi_ortho_tc_two_e_transp(p2,p1,h1, h2))
hji = (mo_bi_ortho_tc_two_e(p1, p2,h1, h2) - mo_bi_ortho_tc_two_e(p2,p1,h1, h2))
if (hij /= 0.d0.or.hji==0.d0) then
! take conjugate to get contribution to <alpha|H|psi> instead of <psi|H|alpha>
! hij = dconjg(hij) * get_phase_bi(phasemask, mi, mi, h1, p1, h2, p2, N_int)
phase = get_phase_bi(phasemask, mi, mi, h1, p1, h2, p2, N_int)
hij = hij * phase
hji = hji* phase
!DIR$ LOOP COUNT AVG(4)
do k=1,N_states
mat_r(k, puti, putj) = mat_r(k, puti, putj) + coefs(k,2) * hij
mat_l(k, puti, putj) = mat_l(k, puti, putj) + coefs(k,1) * hji
enddo
end if
end if
end if
end if
end

4
plugins/local/cipsi_tc_bi_ortho/pt2.irp.f
View File

@ -65,8 +65,12 @@ subroutine tc_pt2
call pt2_dealloc(pt2_data_err) call pt2_dealloc(pt2_data_err)
call pt2_alloc(pt2_data, N_states) call pt2_alloc(pt2_data, N_states)
call pt2_alloc(pt2_data_err, N_states) call pt2_alloc(pt2_data_err, N_states)
if(transpose_two_e_int)then
provide mo_bi_ortho_tc_two_e_transp tc_2e_3idx_coulomb_integrals_transp
endif
call ZMQ_pt2(E_tc, pt2_data, pt2_data_err, relative_error,0) ! Stochastic PT2 and selection call ZMQ_pt2(E_tc, pt2_data, pt2_data_err, relative_error,0) ! Stochastic PT2 and selection
call diagonalize_CI_tc_bi_ortho(ndet, E_tc,norm,pt2_data,print_pt2) call diagonalize_CI_tc_bi_ortho(ndet, E_tc,norm,pt2_data,print_pt2)
call print_summary_tc(psi_energy_with_nucl_rep, pt2_data, pt2_data_err, N_det, N_configuration, N_states, psi_s2)
end end

96
plugins/local/cipsi_tc_bi_ortho/selection.irp.f
View File

@ -636,10 +636,7 @@ subroutine splash_pq(mask, sp, det, i_gen, N_sel, bannedOrb, banned, mat, intere
negMask(i,2) = not(mask(i,2)) negMask(i,2) = not(mask(i,2))
end do end do
! print*,'in selection '
do i = 1, N_sel do i = 1, N_sel
! call debug_det(det(1,1,i),N_int)
! print*,i,dabs(psi_selectors_coef_transp_tc(1,2,i) * psi_selectors_coef_transp_tc(1,1,i))
if(interesting(i) < 0) then if(interesting(i) < 0) then
stop 'prefetch interesting(i) and det(i)' stop 'prefetch interesting(i) and det(i)'
endif endif
@ -691,12 +688,24 @@ subroutine splash_pq(mask, sp, det, i_gen, N_sel, bannedOrb, banned, mat, intere
call get_mask_phase(psi_det_sorted_tc(1,1,interesting(i)), phasemask,N_int) call get_mask_phase(psi_det_sorted_tc(1,1,interesting(i)), phasemask,N_int)
if(nt == 4) then if(nt == 4) then
call get_d2_new(det(1,1,i), phasemask, bannedOrb, banned, mat_l, mat_r, mask, h, p, sp, psi_selectors_coef_transp_tc(1, 1, interesting(i))) if(transpose_two_e_int)then
call get_d2_new_transp(det(1,1,i), phasemask, bannedOrb, banned, mat_l, mat_r, mask, h, p, sp, psi_selectors_coef_transp_tc(1, 1, interesting(i)))
else
call get_d2_new (det(1,1,i), phasemask, bannedOrb, banned, mat_l, mat_r, mask, h, p, sp, psi_selectors_coef_transp_tc(1, 1, interesting(i)))
endif
elseif(nt == 3) then elseif(nt == 3) then
call get_d1_new(det(1,1,i), phasemask, bannedOrb, banned, mat_l, mat_r, mask, h, p, sp, psi_selectors_coef_transp_tc(1, 1, interesting(i))) if(transpose_two_e_int)then
call get_d1_transp(det(1,1,i), phasemask, bannedOrb, banned, mat_l, mat_r, mask, h, p, sp, psi_selectors_coef_transp_tc(1, 1, interesting(i)))
else
call get_d1_new (det(1,1,i), phasemask, bannedOrb, banned, mat_l, mat_r, mask, h, p, sp, psi_selectors_coef_transp_tc(1, 1, interesting(i)))
endif
else
if(transpose_two_e_int)then
call get_d0_transp (det(1,1,i), phasemask, bannedOrb, banned, mat_l, mat_r, mask, h, p, sp, psi_selectors_coef_transp_tc(1, 1, interesting(i)))
else else
call get_d0_new (det(1,1,i), phasemask, bannedOrb, banned, mat_l, mat_r, mask, h, p, sp, psi_selectors_coef_transp_tc(1, 1, interesting(i))) call get_d0_new (det(1,1,i), phasemask, bannedOrb, banned, mat_l, mat_r, mask, h, p, sp, psi_selectors_coef_transp_tc(1, 1, interesting(i)))
endif endif
endif
elseif(nt == 4) then elseif(nt == 4) then
call bitstring_to_list_in_selection(mobMask(1,1), p(1,1), p(0,1), N_int) call bitstring_to_list_in_selection(mobMask(1,1), p(1,1), p(0,1), N_int)
call bitstring_to_list_in_selection(mobMask(1,2), p(1,2), p(0,2), N_int) call bitstring_to_list_in_selection(mobMask(1,2), p(1,2), p(0,2), N_int)
@ -887,79 +896,11 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d
call diag_htilde_mu_mat_fock_bi_ortho(N_int, det, hmono, htwoe, hthree, hii) call diag_htilde_mu_mat_fock_bi_ortho(N_int, det, hmono, htwoe, hthree, hii)
do istate = 1,N_states do istate = 1,N_states
delta_E = E0(istate) - Hii + E_shift delta_E = E0(istate) - Hii + E_shift
double precision :: alpha_h_psi_tmp, psi_h_alpha_tmp, error
if(debug_tc_pt2 == 1)then !! Using the old version
psi_h_alpha = 0.d0
alpha_h_psi = 0.d0
do iii = 1, N_det_selectors
call htilde_mu_mat_bi_ortho_tot_slow(psi_selectors(1,1,iii), det, N_int, i_h_alpha)
call htilde_mu_mat_bi_ortho_tot_slow(det, psi_selectors(1,1,iii), N_int, alpha_h_i)
call get_excitation_degree(psi_selectors(1,1,iii), det,degree,N_int)
if(degree == 0)then
print*,'problem !!!'
print*,'a determinant is already in the wave function !!'
print*,'it corresponds to the selector number ',iii
call debug_det(det,N_int)
stop
endif
! call htilde_mu_mat_opt_bi_ortho_no_3e(psi_selectors(1,1,iii), det, N_int, i_h_alpha)
! call htilde_mu_mat_opt_bi_ortho_no_3e(det, psi_selectors(1,1,iii), N_int, alpha_h_i)
psi_h_alpha += i_h_alpha * psi_selectors_coef_tc(iii,2,1) ! left function
alpha_h_psi += alpha_h_i * psi_selectors_coef_tc(iii,1,1) ! right function
enddo
else if(debug_tc_pt2 == 2)then !! debugging the new version
! psi_h_alpha_tmp = 0.d0
! alpha_h_psi_tmp = 0.d0
! do iii = 1, N_det_selectors ! old version
! call htilde_mu_mat_opt_bi_ortho_no_3e(psi_selectors(1,1,iii), det, N_int, i_h_alpha)
! call htilde_mu_mat_opt_bi_ortho_no_3e(det, psi_selectors(1,1,iii), N_int, alpha_h_i)
! psi_h_alpha_tmp += i_h_alpha * psi_selectors_coef_tc(iii,1,1) ! left function
! alpha_h_psi_tmp += alpha_h_i * psi_selectors_coef_tc(iii,2,1) ! right function
! enddo
psi_h_alpha_tmp = mat_l(istate, p1, p2) ! new version
alpha_h_psi_tmp = mat_r(istate, p1, p2) ! new version
psi_h_alpha = 0.d0
alpha_h_psi = 0.d0
do iii = 1, N_det ! old version
call htilde_mu_mat_opt_bi_ortho_no_3e(psi_det(1,1,iii), det, N_int, i_h_alpha)
call htilde_mu_mat_opt_bi_ortho_no_3e(det, psi_det(1,1,iii), N_int, alpha_h_i)
psi_h_alpha += i_h_alpha * psi_l_coef_bi_ortho(iii,1) ! left function
alpha_h_psi += alpha_h_i * psi_r_coef_bi_ortho(iii,1) ! right function
enddo
if(dabs(psi_h_alpha*alpha_h_psi/delta_E).gt.1.d-10)then
error = dabs(psi_h_alpha * alpha_h_psi - psi_h_alpha_tmp * alpha_h_psi_tmp)/dabs(psi_h_alpha * alpha_h_psi)
if(error.gt.1.d-2)then
call debug_det(det, N_int)
print*,'error =',error,psi_h_alpha * alpha_h_psi/delta_E,psi_h_alpha_tmp * alpha_h_psi_tmp/delta_E
print*,psi_h_alpha , alpha_h_psi
print*,psi_h_alpha_tmp , alpha_h_psi_tmp
print*,'selectors '
do iii = 1, N_det_selectors ! old version
print*,'iii',iii,psi_selectors_coef_tc(iii,1,1),psi_selectors_coef_tc(iii,2,1)
call htilde_mu_mat_opt_bi_ortho_no_3e(psi_selectors(1,1,iii), det, N_int, i_h_alpha)
call htilde_mu_mat_opt_bi_ortho_no_3e(det, psi_selectors(1,1,iii), N_int, alpha_h_i)
print*,i_h_alpha,alpha_h_i
call debug_det(psi_selectors(1,1,iii),N_int)
enddo
! print*,'psi_det '
! do iii = 1, N_det! old version
! print*,'iii',iii,psi_l_coef_bi_ortho(iii,1),psi_r_coef_bi_ortho(iii,1)
! call debug_det(psi_det(1,1,iii),N_int)
! enddo
stop
endif
endif
else
psi_h_alpha = mat_l(istate, p1, p2) psi_h_alpha = mat_l(istate, p1, p2)
alpha_h_psi = mat_r(istate, p1, p2) alpha_h_psi = mat_r(istate, p1, p2)
endif
val = 4.d0 * psi_h_alpha * alpha_h_psi val = 4.d0 * psi_h_alpha * alpha_h_psi
tmp = dsqrt(delta_E * delta_E + val) tmp = dsqrt(delta_E * delta_E + val)
! if (delta_E < 0.d0) then
! tmp = -tmp
! endif
e_pert(istate) = 0.25 * val / delta_E e_pert(istate) = 0.25 * val / delta_E
! e_pert(istate) = 0.5d0 * (tmp - delta_E)
if(dsqrt(tmp).gt.1.d-4.and.dabs(psi_h_alpha).gt.1.d-4)then if(dsqrt(tmp).gt.1.d-4.and.dabs(psi_h_alpha).gt.1.d-4)then
coef(istate) = e_pert(istate) / psi_h_alpha coef(istate) = e_pert(istate) / psi_h_alpha
else else
@ -976,15 +917,6 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d
if(e_pert(istate).gt.0.d0)e_pert(istate)=0.d0 if(e_pert(istate).gt.0.d0)e_pert(istate)=0.d0
endif endif
! if(selection_tc == 1 )then
! if(e_pert(istate).lt.0.d0)then
! e_pert(istate) = 0.d0
! endif
! else if(selection_tc == -1)then
! if(e_pert(istate).gt.0.d0)then
! e_pert(istate) = 0.d0
! endif
! endif
enddo enddo

3
plugins/local/cipsi_tc_bi_ortho/stochastic_cipsi.irp.f
View File

@ -88,6 +88,9 @@ subroutine run_stochastic_cipsi
call pt2_dealloc(pt2_data_err) call pt2_dealloc(pt2_data_err)
call pt2_alloc(pt2_data, N_states) call pt2_alloc(pt2_data, N_states)
call pt2_alloc(pt2_data_err, N_states) call pt2_alloc(pt2_data_err, N_states)
if(transpose_two_e_int)then
provide mo_bi_ortho_tc_two_e_transp tc_2e_3idx_coulomb_integrals_transp
endif
call ZMQ_pt2(E_tc, pt2_data, pt2_data_err, relative_error,to_select) ! Stochastic PT2 and selection call ZMQ_pt2(E_tc, pt2_data, pt2_data_err, relative_error,to_select) ! Stochastic PT2 and selection
! stop ! stop

2
plugins/local/fci_tc_bi/pt2_tc.irp.f
View File

@ -13,6 +13,8 @@ program tc_pt2_prog
pruning = -1.d0 pruning = -1.d0
touch pruning touch pruning
read_wf = .True.
touch read_wf
! pt2_relative_error = 0.01d0 ! pt2_relative_error = 0.01d0
! touch pt2_relative_error ! touch pt2_relative_error

101
plugins/local/non_h_ints_mu/deb_mos.irp.f Normal file
View File

@ -0,0 +1,101 @@
! ---
program deb_mos
implicit none
my_grid_becke = .True.
PROVIDE tc_grid1_a tc_grid1_r
my_n_pt_r_grid = tc_grid1_r
my_n_pt_a_grid = tc_grid1_a
touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
if(tc_integ_type .eq. "numeric") then
my_extra_grid_becke = .True.
PROVIDE tc_grid2_a tc_grid2_r
my_n_pt_r_extra_grid = tc_grid2_r
my_n_pt_a_extra_grid = tc_grid2_a
touch my_extra_grid_becke my_n_pt_r_extra_grid my_n_pt_a_extra_grid
endif
call print_mos()
end
! ---
subroutine print_mos()
implicit none
integer :: i, ipoint
double precision :: r(3)
double precision :: mo_val, mo_der(3), mo_lap
PROVIDE final_grid_points mos_in_r_array mos_grad_in_r_array mos_lapl_in_r_array
! do ipoint = 1, n_points_final_grid
! r(:) = final_grid_points(:,ipoint)
! print*, r
! enddo
double precision :: accu_vgl(5)
double precision :: accu_vgl_nrm(5)
do ipoint = 1, n_points_final_grid
r(1) = final_grid_points(1,i)
r(2) = final_grid_points(2,i)
r(3) = final_grid_points(3,i)
write(1111, '(5(f15.7, 3X))') r
do i = 1, mo_num
mo_val = mos_in_r_array (i,ipoint)
mo_der(:) = mos_grad_in_r_array(i,ipoint,:)
mo_lap = mos_lapl_in_r_array(i,ipoint,1) + mos_lapl_in_r_array(i,ipoint,2) + mos_lapl_in_r_array(i,ipoint,3)
write(1111, '(5(f15.7, 3X))') mo_val, mo_der, mo_lap
enddo
enddo
do ipoint = 1, n_points_final_grid
r(1) = final_grid_points(1,i)
r(2) = final_grid_points(2,i)
r(3) = final_grid_points(3,i)
write(2222, '(5(f15.7, 3X))') r
do i = 1, mo_num
mo_val = mos_in_r_array_qmckl (i,ipoint)
mo_der(:) = mos_grad_in_r_array_qmckl(i,ipoint,:)
mo_lap = mos_lapl_in_r_array_qmckl(i,ipoint)
write(2222, '(5(f15.7, 3X))') mo_val, mo_der, mo_lap
enddo
enddo
accu_vgl = 0.d0
accu_vgl_nrm = 0.d0
do ipoint = 1, n_points_final_grid
do i = 1, mo_num
mo_val = mos_in_r_array (i,ipoint)
mo_der(:) = mos_grad_in_r_array(i,ipoint,:)
mo_lap = mos_lapl_in_r_array(i,ipoint,1) + mos_lapl_in_r_array(i,ipoint,2) + mos_lapl_in_r_array(i,ipoint,3)
accu_vgl_nrm(1) += dabs(mo_val)
accu_vgl_nrm(2) += dabs(mo_der(1))
accu_vgl_nrm(3) += dabs(mo_der(2))
accu_vgl_nrm(4) += dabs(mo_der(3))
accu_vgl_nrm(5) += dabs(mo_lap)
mo_val -= mos_in_r_array_qmckl (i,ipoint)
mo_der(:) -= mos_grad_in_r_array_qmckl(i,ipoint,:)
mo_lap -= mos_lapl_in_r_array_qmckl(i,ipoint)
accu_vgl(1) += dabs(mo_val)
accu_vgl(2) += dabs(mo_der(1))
accu_vgl(3) += dabs(mo_der(2))
accu_vgl(4) += dabs(mo_der(3))
accu_vgl(5) += dabs(mo_lap)
enddo
enddo
accu_vgl(:) *= 1.d0 / accu_vgl_nrm(:)
print *, accu_vgl
return
end
! ---

26
plugins/local/non_h_ints_mu/jast_deriv_utils_vect.irp.f
View File

@ -340,8 +340,8 @@ subroutine grad1_j12_r1_seq(r1, n_grid2, gradx, grady, gradz)
endif endif
tmp1 = double_p(mpA) * f1A_power(mpA-1) * f2A_power(npA) + double_p(npA) * f1A_power(npA-1) * f2A_power(mpA) tmp1 = double_p(mpA) * f1A_power(mpA-1) * f2A_power(npA) + double_p(npA) * f1A_power(npA-1) * f2A_power(mpA)
tmp1 = tmp1 * g12_power(opA) tmp1 = tmp1 * g12_power(opA) * tmp
tmp2 = double_p(opA) * g12_power(opA-1) * (f1A_power(mpA) * f2A_power(npA) + f1A_power(npA) * f2A_power(mpA)) tmp2 = double_p(opA) * g12_power(opA-1) * (f1A_power(mpA) * f2A_power(npA) + f1A_power(npA) * f2A_power(mpA)) * tmp
!tmp1 = 0.d0 !tmp1 = 0.d0
!if(mpA .gt. 0) then !if(mpA .gt. 0) then
@ -356,9 +356,12 @@ subroutine grad1_j12_r1_seq(r1, n_grid2, gradx, grady, gradz)
! tmp2 = tmp2 + dble(opA) * g12**dble(opA-1) * (f1A**dble(mpA) * f2A**dble(npA) + f1A**dble(npA) * f2A**dble(mpA)) ! tmp2 = tmp2 + dble(opA) * g12**dble(opA-1) * (f1A**dble(mpA) * f2A**dble(npA) + f1A**dble(npA) * f2A**dble(mpA))
!endif !endif
gradx(jpoint) = gradx(jpoint) + tmp * (tmp1 * grad1_f1A(1) + tmp2 * grad1_g12(1)) ! gradx(jpoint) = gradx(jpoint) + tmp * (tmp1 * grad1_f1A(1) + tmp2 * grad1_g12(1))
grady(jpoint) = grady(jpoint) + tmp * (tmp1 * grad1_f1A(2) + tmp2 * grad1_g12(2)) ! grady(jpoint) = grady(jpoint) + tmp * (tmp1 * grad1_f1A(2) + tmp2 * grad1_g12(2))
gradz(jpoint) = gradz(jpoint) + tmp * (tmp1 * grad1_f1A(3) + tmp2 * grad1_g12(3)) ! gradz(jpoint) = gradz(jpoint) + tmp * (tmp1 * grad1_f1A(3) + tmp2 * grad1_g12(3))
gradx(jpoint) = gradx(jpoint) + tmp1 * grad1_f1A(1) + tmp2 * grad1_g12(1)
grady(jpoint) = grady(jpoint) + tmp1 * grad1_f1A(2) + tmp2 * grad1_g12(2)
gradz(jpoint) = gradz(jpoint) + tmp1 * grad1_f1A(3) + tmp2 * grad1_g12(3)
enddo ! p enddo ! p
enddo ! i_nucl enddo ! i_nucl
enddo ! jpoint enddo ! jpoint
@ -864,19 +867,20 @@ subroutine jBH_elem_fct_grad(alpha, r1, r2, fct, grad1_fct)
+ (r1(2) - r2(2)) * (r1(2) - r2(2)) & + (r1(2) - r2(2)) * (r1(2) - r2(2)) &
+ (r1(3) - r2(3)) * (r1(3) - r2(3)) ) + (r1(3) - r2(3)) * (r1(3) - r2(3)) )
if(dist .ge. 1d-10) then
tmp1 = 1.d0 / (1.d0 + alpha * dist) tmp1 = 1.d0 / (1.d0 + alpha * dist)
fct = alpha * dist * tmp1 fct = alpha * dist * tmp1
if(dist .lt. 1d-10) then
grad1_fct(1) = 0.d0
grad1_fct(2) = 0.d0
grad1_fct(3) = 0.d0
else
tmp2 = alpha * tmp1 * tmp1 / dist tmp2 = alpha * tmp1 * tmp1 / dist
grad1_fct(1) = tmp2 * (r1(1) - r2(1)) grad1_fct(1) = tmp2 * (r1(1) - r2(1))
grad1_fct(2) = tmp2 * (r1(2) - r2(2)) grad1_fct(2) = tmp2 * (r1(2) - r2(2))
grad1_fct(3) = tmp2 * (r1(3) - r2(3)) grad1_fct(3) = tmp2 * (r1(3) - r2(3))
else
grad1_fct(1) = 0.d0
grad1_fct(2) = 0.d0
grad1_fct(3) = 0.d0
fct = 0.d0
endif endif
return return

2
plugins/local/non_h_ints_mu/qmckl.irp.f
View File

@ -158,7 +158,7 @@ END_PROVIDER
double precision, allocatable :: vgl(:,:,:) double precision, allocatable :: vgl(:,:,:)
allocate( vgl(mo_num,5,n_points_final_grid)) allocate( vgl(mo_num,5,n_points_final_grid))
rc = qmckl_get_mo_basis_mo_vgl_inplace(qmckl_ctx, vgl, n_points_final_grid*mo_num*5_8) rc = qmckl_get_mo_basis_mo_vgl(qmckl_ctx, vgl, n_points_final_grid*mo_num*5_8)
if (rc /= QMCKL_SUCCESS) then if (rc /= QMCKL_SUCCESS) then
print *, irp_here, 'qmckl error in get_mo_vgl' print *, irp_here, 'qmckl error in get_mo_vgl'
rc = qmckl_check(qmckl_ctx, rc) rc = qmckl_check(qmckl_ctx, rc)

2
plugins/local/normal_order_old/NEED
View File

@ -1 +1 @@
tc_scf

1
plugins/local/old_delta_tc_qmc/NEED Normal file
View File

@ -0,0 +1 @@

4
plugins/local/old_delta_tc_qmc/README.rst Normal file
View File

@ -0,0 +1,4 @@
================
old_delta_tc_qmc
================

0
plugins/local/tc_bi_ortho/compute_deltamu_right.irp.f → plugins/local/old_delta_tc_qmc/compute_deltamu_right.irp.f
View File

8
plugins/local/tc_bi_ortho/dressing_vectors_lr.irp.f → plugins/local/old_delta_tc_qmc/dressing_vectors_lr.irp.f
View File

@ -27,7 +27,7 @@ subroutine get_delta_bitc_right(psidet, psicoef, ndet, Nint, delta)
i = 1 i = 1
j = 1 j = 1
call htilde_mu_mat_bi_ortho_slow(psidet(1,1,i), psidet(1,1,j), Nint, htc_mono, htc_twoe, htc_three, htc_tot) call htilde_mu_mat_opt_bi_ortho(psidet(1,1,i), psidet(1,1,j), Nint, htc_mono, htc_twoe, htc_three, htc_tot)
call hmat_bi_ortho (psidet(1,1,i), psidet(1,1,j), Nint, h_mono, h_twoe, h_tot) call hmat_bi_ortho (psidet(1,1,i), psidet(1,1,j), Nint, h_mono, h_twoe, h_tot)
delta = 0.d0 delta = 0.d0
@ -39,7 +39,7 @@ subroutine get_delta_bitc_right(psidet, psicoef, ndet, Nint, delta)
do j = 1, ndet do j = 1, ndet
! < I |Htilde | J > ! < I |Htilde | J >
call htilde_mu_mat_bi_ortho_slow(psidet(1,1,i), psidet(1,1,j), Nint, htc_mono, htc_twoe, htc_three, htc_tot) call htilde_mu_mat_opt_bi_ortho(psidet(1,1,i), psidet(1,1,j), Nint, htc_mono, htc_twoe, htc_three, htc_tot)
! < I |H | J > ! < I |H | J >
call hmat_bi_ortho(psidet(1,1,i), psidet(1,1,j), Nint, h_mono, h_twoe, h_tot) call hmat_bi_ortho(psidet(1,1,i), psidet(1,1,j), Nint, h_mono, h_twoe, h_tot)
@ -78,7 +78,7 @@ subroutine get_htc_bitc_right(psidet, psicoef, ndet, Nint, delta)
i = 1 i = 1
j = 1 j = 1
call htilde_mu_mat_bi_ortho_slow(psidet(1,1,i), psidet(1,1,j), Nint, htc_mono, htc_twoe, htc_three, htc_tot) call htilde_mu_mat_opt_bi_ortho(psidet(1,1,i), psidet(1,1,j), Nint, htc_mono, htc_twoe, htc_three, htc_tot)
delta = 0.d0 delta = 0.d0
!$OMP PARALLEL DO DEFAULT(NONE) SCHEDULE(dynamic,8) & !$OMP PARALLEL DO DEFAULT(NONE) SCHEDULE(dynamic,8) &
@ -88,7 +88,7 @@ subroutine get_htc_bitc_right(psidet, psicoef, ndet, Nint, delta)
do j = 1, ndet do j = 1, ndet
! < I |Htilde | J > ! < I |Htilde | J >
call htilde_mu_mat_bi_ortho_slow(psidet(1,1,i), psidet(1,1,j), Nint, htc_mono, htc_twoe, htc_three, htc_tot) call htilde_mu_mat_opt_bi_ortho(psidet(1,1,i), psidet(1,1,j), Nint, htc_mono, htc_twoe, htc_three, htc_tot)
delta(i) = delta(i) + psicoef(j) * htc_tot delta(i) = delta(i) + psicoef(j) * htc_tot
enddo enddo

2
plugins/local/slater_tc/slater_tc.irp.f → plugins/local/old_delta_tc_qmc/old_delta_tc_qmc.irp.f
View File

@ -1,4 +1,4 @@
program slater_tc program old_delta_tc_qmc
implicit none implicit none
BEGIN_DOC BEGIN_DOC
! TODO : Put the documentation of the program here ! TODO : Put the documentation of the program here

193
plugins/local/slater_tc/h_mat_triple.irp.f
View File

@ -1,196 +1,3 @@
subroutine get_excitation_general(key_i,key_j, Nint,degree_array,holes_array, particles_array,phase)
use bitmasks
BEGIN_DOC
! returns the array, for each spin, of holes/particles between key_i and key_j
!
! with the following convention: a^+_{particle} a_{hole}|key_i> = |key_j>
END_DOC
include 'utils/constants.include.F'
implicit none
integer, intent(in) :: Nint
integer(bit_kind), intent(in) :: key_j(Nint,2),key_i(Nint,2)
integer, intent(out) :: holes_array(100,2),particles_array(100,2),degree_array(2)
double precision, intent(out) :: phase
integer :: ispin,k,i,pos
integer(bit_kind) :: key_hole, key_particle
integer(bit_kind) :: xorvec(N_int_max,2)
holes_array = -1
particles_array = -1
degree_array = 0
do i = 1, N_int
xorvec(i,1) = xor( key_i(i,1), key_j(i,1))
xorvec(i,2) = xor( key_i(i,2), key_j(i,2))
degree_array(1) += popcnt(xorvec(i,1))
degree_array(2) += popcnt(xorvec(i,2))
enddo
degree_array(1) = shiftr(degree_array(1),1)
degree_array(2) = shiftr(degree_array(2),1)
do ispin = 1, 2
k = 1
!!! GETTING THE HOLES
do i = 1, N_int
key_hole = iand(xorvec(i,ispin),key_i(i,ispin))
do while(key_hole .ne.0_bit_kind)
pos = trailz(key_hole)
holes_array(k,ispin) = 1+ bit_kind_size * (i-1) + pos
key_hole = ibclr(key_hole,pos)
k += 1
if(k .gt.100)then
print*,'WARNING in get_excitation_general'
print*,'More than a 100-th excitation for spin ',ispin
print*,'stoping ...'
stop
endif
enddo
enddo
enddo
do ispin = 1, 2
k = 1
!!! GETTING THE PARTICLES
do i = 1, N_int
key_particle = iand(xor(key_i(i,ispin),key_j(i,ispin)),key_j(i,ispin))
do while(key_particle .ne.0_bit_kind)
pos = trailz(key_particle)
particles_array(k,ispin) = 1+ bit_kind_size * (i-1) + pos
key_particle = ibclr(key_particle,pos)
k += 1
if(k .gt.100)then
print*,'WARNING in get_excitation_general '
print*,'More than a 100-th excitation for spin ',ispin
print*,'stoping ...'
stop
endif
enddo
enddo
enddo
integer :: h,p, i_ok
integer(bit_kind), allocatable :: det_i(:,:),det_ip(:,:)
integer :: exc(0:2,2,2)
double precision :: phase_tmp
allocate(det_i(Nint,2),det_ip(N_int,2))
det_i = key_i
phase = 1.d0
do ispin = 1, 2
do i = 1, degree_array(ispin)
h = holes_array(i,ispin)
p = particles_array(i,ispin)
det_ip = det_i
call do_single_excitation(det_ip,h,p,ispin,i_ok)
if(i_ok == -1)then
print*,'excitation was not possible '
stop
endif
call get_single_excitation(det_i,det_ip,exc,phase_tmp,Nint)
phase *= phase_tmp
det_i = det_ip
enddo
enddo
end
subroutine get_holes_general(key_i, key_j,Nint, holes_array)
use bitmasks
BEGIN_DOC
! returns the array, per spin, of holes between key_i and key_j
!
! with the following convention: a_{hole}|key_i> --> |key_j>
END_DOC
implicit none
integer, intent(in) :: Nint
integer(bit_kind), intent(in) :: key_j(Nint,2),key_i(Nint,2)
integer, intent(out) :: holes_array(100,2)
integer(bit_kind) :: key_hole
integer :: ispin,k,i,pos
holes_array = -1
do ispin = 1, 2
k = 1
do i = 1, N_int
key_hole = iand(xor(key_i(i,ispin),key_j(i,ispin)),key_i(i,ispin))
do while(key_hole .ne.0_bit_kind)
pos = trailz(key_hole)
holes_array(k,ispin) = 1+ bit_kind_size * (i-1) + pos
key_hole = ibclr(key_hole,pos)
k += 1
if(k .gt.100)then
print*,'WARNING in get_holes_general'
print*,'More than a 100-th excitation for spin ',ispin
print*,'stoping ...'
stop
endif
enddo
enddo
enddo
end
subroutine get_particles_general(key_i, key_j,Nint,particles_array)
use bitmasks
BEGIN_DOC
! returns the array, per spin, of particles between key_i and key_j
!
! with the following convention: a^dagger_{particle}|key_i> --> |key_j>
END_DOC
implicit none
integer, intent(in) :: Nint
integer(bit_kind), intent(in) :: key_j(Nint,2),key_i(Nint,2)
integer, intent(out) :: particles_array(100,2)
integer(bit_kind) :: key_particle
integer :: ispin,k,i,pos
particles_array = -1
do ispin = 1, 2
k = 1
do i = 1, N_int
key_particle = iand(xor(key_i(i,ispin),key_j(i,ispin)),key_j(i,ispin))
do while(key_particle .ne.0_bit_kind)
pos = trailz(key_particle)
particles_array(k,ispin) = 1+ bit_kind_size * (i-1) + pos
key_particle = ibclr(key_particle,pos)
k += 1
if(k .gt.100)then
print*,'WARNING in get_holes_general'
print*,'More than a 100-th excitation for spin ',ispin
print*,'Those are the two determinants'
call debug_det(key_i, N_int)
call debug_det(key_j, N_int)
print*,'stoping ...'
stop
endif
enddo
enddo
enddo
end
subroutine get_phase_general(key_i,Nint,degree, holes_array, particles_array,phase)
implicit none
integer, intent(in) :: degree(2), Nint
integer(bit_kind), intent(in) :: key_i(Nint,2)
integer, intent(in) :: holes_array(100,2),particles_array(100,2)
double precision, intent(out) :: phase
integer :: i,ispin,h,p, i_ok
integer(bit_kind), allocatable :: det_i(:,:),det_ip(:,:)
integer :: exc(0:2,2,2)
double precision :: phase_tmp
allocate(det_i(Nint,2),det_ip(N_int,2))
det_i = key_i
phase = 1.d0
do ispin = 1, 2
do i = 1, degree(ispin)
h = holes_array(i,ispin)
p = particles_array(i,ispin)
det_ip = det_i
call do_single_excitation(det_ip,h,p,ispin,i_ok)
if(i_ok == -1)then
print*,'excitation was not possible '
stop
endif
call get_single_excitation(det_i,det_ip,exc,phase_tmp,Nint)
phase *= phase_tmp
det_i = det_ip
enddo
enddo
end
subroutine H_tc_s2_u_0_with_pure_three(v_0, s_0, u_0, N_st, sze) subroutine H_tc_s2_u_0_with_pure_three(v_0, s_0, u_0, N_st, sze)
BEGIN_DOC BEGIN_DOC
! Computes $v_0 = H^TC | u_0\rangle$ WITH PURE TRIPLE EXCITATION TERMS ! Computes $v_0 = H^TC | u_0\rangle$ WITH PURE TRIPLE EXCITATION TERMS

45
plugins/local/slater_tc/slater_tc_opt.irp.f
View File

@ -181,3 +181,48 @@ end
! --- ! ---
subroutine htilde_mu_mat_opt_bi_ortho_no_3e_both(key_j, key_i, Nint, hji,hij)
BEGIN_DOC
!
! <key_j |H_tilde | key_i> where |key_j> is developed on the LEFT basis and |key_i> is developed on the RIGHT basis
!!
! Returns the detail of the matrix element WITHOUT ANY CONTRIBUTION FROM THE THREE ELECTRON TERMS
!! WARNING !!
!
! Non hermitian !!
!
END_DOC
use bitmasks
implicit none
integer, intent(in) :: Nint
integer(bit_kind), intent(in) :: key_i(Nint,2), key_j(Nint,2)
double precision, intent(out) :: hji,hij
integer :: degree
hji = 0.d0
hij = 0.d0
call get_excitation_degree(key_i, key_j, degree, Nint)
if(degree.gt.2) return
if(degree == 0) then
call diag_htilde_mu_mat_fock_bi_ortho_no_3e(Nint, key_i,hji)
hij = hji
else if (degree == 1) then
call single_htilde_mu_mat_fock_bi_ortho_no_3e_both(Nint,key_j, key_i , hji,hij)
else if(degree == 2) then
call double_htilde_mu_mat_fock_bi_ortho_no_3e_both(Nint, key_j, key_i, hji,hij)
endif
if(degree==0) then
hji += nuclear_repulsion
hij += nuclear_repulsion
endif
end
! ---

311
plugins/local/slater_tc/slater_tc_opt_diag.irp.f
View File

@ -19,13 +19,13 @@
PROVIDE HF_bitmask PROVIDE HF_bitmask
PROVIDE mo_l_coef mo_r_coef PROVIDE mo_l_coef mo_r_coef
call diag_htilde_mu_mat_bi_ortho_slow(N_int, HF_bitmask, hmono, htwoe, htot) call diag_htc_bi_orth_2e_brute(N_int, HF_bitmask, hmono, htwoe, htot)
ref_tc_energy_1e = hmono ref_tc_energy_1e = hmono
ref_tc_energy_2e = htwoe ref_tc_energy_2e = htwoe
if(three_body_h_tc) then if(three_body_h_tc) then
call diag_htilde_three_body_ints_bi_ort_slow(N_int, HF_bitmask, hthree) call diag_htc_bi_orth_3e_brute(N_int, HF_bitmask, hthree)
ref_tc_energy_3e = hthree ref_tc_energy_3e = hthree
else else
ref_tc_energy_3e = 0.d0 ref_tc_energy_3e = 0.d0
@ -524,3 +524,310 @@ end
! --- ! ---
subroutine diag_htc_bi_orth_2e_brute(Nint, key_i, hmono, htwoe, htot)
BEGIN_DOC
!
! diagonal element of htilde ONLY FOR ONE- AND TWO-BODY TERMS
!
END_DOC
use bitmasks
implicit none
integer, intent(in) :: Nint
integer(bit_kind), intent(in) :: key_i(Nint,2)
double precision, intent(out) :: hmono,htwoe,htot
integer :: occ(Nint*bit_kind_size,2)
integer :: Ne(2), i, j, ii, jj, ispin, jspin, k, kk
double precision :: get_mo_two_e_integral_tc_int
integer(bit_kind) :: key_i_core(Nint,2)
PROVIDE mo_bi_ortho_tc_two_e
hmono = 0.d0
htwoe = 0.d0
htot = 0.d0
call bitstring_to_list_ab(key_i, occ, Ne, Nint)
do ispin = 1, 2
do i = 1, Ne(ispin)
ii = occ(i,ispin)
hmono += mo_bi_ortho_tc_one_e(ii,ii)
enddo
enddo
! alpha/beta two-body
ispin = 1
jspin = 2
do i = 1, Ne(ispin) ! electron 1 (so it can be associated to mu(r1))
ii = occ(i,ispin)
do j = 1, Ne(jspin) ! electron 2
jj = occ(j,jspin)
htwoe += mo_bi_ortho_tc_two_e(jj,ii,jj,ii)
enddo
enddo
! alpha/alpha two-body
do i = 1, Ne(ispin)
ii = occ(i,ispin)
do j = i+1, Ne(ispin)
jj = occ(j,ispin)
htwoe += mo_bi_ortho_tc_two_e(ii,jj,ii,jj) - mo_bi_ortho_tc_two_e(ii,jj,jj,ii)
enddo
enddo
! beta/beta two-body
do i = 1, Ne(jspin)
ii = occ(i,jspin)
do j = i+1, Ne(jspin)
jj = occ(j,jspin)
htwoe += mo_bi_ortho_tc_two_e(ii,jj,ii,jj) - mo_bi_ortho_tc_two_e(ii,jj,jj,ii)
enddo
enddo
htot = hmono + htwoe
end
! ---
subroutine diag_htc_bi_orth_3e_brute(Nint, key_i, hthree)
BEGIN_DOC
! diagonal element of htilde ONLY FOR THREE-BODY TERMS WITH BI ORTHONORMAL ORBITALS
END_DOC
use bitmasks
implicit none
integer, intent(in) :: Nint
integer(bit_kind), intent(in) :: key_i(Nint,2)
double precision, intent(out) :: hthree
integer :: occ(Nint*bit_kind_size,2)
integer :: Ne(2),i,j,ii,jj,ispin,jspin,m,mm
integer(bit_kind) :: key_i_core(Nint,2)
double precision :: direct_int, exchange_int, ref
double precision, external :: sym_3_e_int_from_6_idx_tensor
double precision, external :: three_e_diag_parrallel_spin
PROVIDE mo_l_coef mo_r_coef
if(core_tc_op) then
do i = 1, Nint
key_i_core(i,1) = xor(key_i(i,1), core_bitmask(i,1))
key_i_core(i,2) = xor(key_i(i,2), core_bitmask(i,2))
enddo
call bitstring_to_list_ab(key_i_core, occ, Ne, Nint)
else
call bitstring_to_list_ab(key_i, occ, Ne, Nint)
endif
hthree = 0.d0
if((Ne(1)+Ne(2)) .ge. 3) then
! alpha/alpha/beta three-body
do i = 1, Ne(1)
ii = occ(i,1)
do j = i+1, Ne(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
hthree += direct_int - exchange_int
enddo
enddo
enddo
! beta/beta/alpha three-body
do i = 1, Ne(2)
ii = occ(i,2)
do j = i+1, Ne(2)
jj = occ(j,2)
do m = 1, Ne(1)
mm = occ(m,1)
!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)
exchange_int = three_e_3_idx_exch12_bi_ort(mm,jj,ii)
hthree += direct_int - exchange_int
enddo
enddo
enddo
! alpha/alpha/alpha three-body
do i = 1, Ne(1)
ii = occ(i,1) ! 1
do j = i+1, Ne(1)
jj = occ(j,1) ! 2
do m = j+1, Ne(1)
mm = occ(m,1) ! 3
!hthree += 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
enddo
! beta/beta/beta three-body
do i = 1, Ne(2)
ii = occ(i,2) ! 1
do j = i+1, Ne(2)
jj = occ(j,2) ! 2
do m = j+1, Ne(2)
mm = occ(m,2) ! 3
!hthree += 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
enddo
endif
end
BEGIN_PROVIDER [ double precision, three_e_diag_parrallel_spin_prov, (mo_num, mo_num, mo_num)]
BEGIN_DOC
!
! matrix element of the -L three-body operator ON A BI ORTHONORMAL BASIS
!
! three_e_diag_parrallel_spin_prov(m,j,i) = All combinations of the form <mji|-L|mji> for same spin matrix elements
!
! notice the -1 sign: in this way three_e_diag_parrallel_spin_prov can be directly used to compute Slater rules with a + sign
!
END_DOC
implicit none
integer :: i, j, m
double precision :: integral, wall1, wall0, three_e_diag_parrallel_spin
three_e_diag_parrallel_spin_prov = 0.d0
print *, ' Providing the three_e_diag_parrallel_spin_prov ...'
integral = three_e_diag_parrallel_spin(1,1,1) ! to provide all stuffs
call wall_time(wall0)
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,j,m,integral) &
!$OMP SHARED (mo_num,three_e_diag_parrallel_spin_prov)
!$OMP DO SCHEDULE (dynamic)
do i = 1, mo_num
do j = 1, mo_num
do m = j, mo_num
three_e_diag_parrallel_spin_prov(m,j,i) = three_e_diag_parrallel_spin(m,j,i)
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
do i = 1, mo_num
do j = 1, mo_num
do m = 1, j
three_e_diag_parrallel_spin_prov(m,j,i) = three_e_diag_parrallel_spin_prov(j,m,i)
enddo
enddo
enddo
call wall_time(wall1)
print *, ' wall time for three_e_diag_parrallel_spin_prov', wall1 - wall0
END_PROVIDER
BEGIN_PROVIDER [ double precision, three_e_single_parrallel_spin_prov, (mo_num, mo_num, mo_num, mo_num)]
BEGIN_DOC
!
! matrix element of the -L three-body operator FOR THE DIRECT TERMS OF SINGLE EXCITATIONS AND BI ORTHO MOs
!
! three_e_single_parrallel_spin_prov(m,j,k,i) = All combination of <mjk|-L|mji> for same spin matrix elements
!
! 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
integer :: i, j, k, m
double precision :: integral, wall1, wall0, three_e_single_parrallel_spin
three_e_single_parrallel_spin_prov = 0.d0
print *, ' Providing the three_e_single_parrallel_spin_prov ...'
integral = three_e_single_parrallel_spin(1,1,1,1)
call wall_time(wall0)
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,j,k,m,integral) &
!$OMP SHARED (mo_num,three_e_single_parrallel_spin_prov)
!$OMP DO SCHEDULE (dynamic)
do i = 1, mo_num
do k = 1, mo_num
do j = 1, mo_num
do m = 1, mo_num
three_e_single_parrallel_spin_prov(m,j,k,i) = three_e_single_parrallel_spin(m,j,k,i)
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
call wall_time(wall1)
print *, ' wall time for three_e_single_parrallel_spin_prov', wall1 - wall0
END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, three_e_double_parrallel_spin_prov, (mo_num, mo_num, mo_num, mo_num, mo_num)]
BEGIN_DOC
!
! matrix element of the -L three-body operator FOR THE DIRECT TERMS OF DOUBLE EXCITATIONS AND BI ORTHO MOs
!
! three_e_double_parrallel_spin_prov(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
integer :: i, j, k, m, l
double precision :: integral, wall1, wall0, three_e_double_parrallel_spin
three_e_double_parrallel_spin_prov = 0.d0
print *, ' Providing the three_e_double_parrallel_spin_prov ...'
call wall_time(wall0)
integral = three_e_double_parrallel_spin(1,1,1,1,1)
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,j,k,m,l,integral) &
!$OMP SHARED (mo_num,three_e_double_parrallel_spin_prov)
!$OMP DO SCHEDULE (dynamic)
do i = 1, mo_num
do k = 1, mo_num
do j = 1, mo_num
do l = 1, mo_num
do m = 1, mo_num
three_e_double_parrallel_spin_prov(m,l,j,k,i) = three_e_double_parrallel_spin(m,l,j,k,i)
enddo
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
call wall_time(wall1)
print *, ' wall time for three_e_double_parrallel_spin_prov', wall1 - wall0
END_PROVIDER

60
plugins/local/slater_tc/slater_tc_opt_double.irp.f
View File

@ -505,3 +505,63 @@ subroutine double_htilde_mu_mat_fock_bi_ortho_no_3e(Nint, key_j, key_i, htot)
end end
subroutine double_htilde_mu_mat_fock_bi_ortho_no_3e_both(Nint, key_j, key_i, hji,hij)
BEGIN_DOC
! <key_j |H_tilde | key_i> and <key_i |H_tilde | key_j> for double excitation ONLY FOR ONE- AND TWO-BODY TERMS
!!
!! WARNING !!
!
! Non hermitian !!
END_DOC
use bitmasks
implicit none
integer, intent(in) :: Nint
integer(bit_kind), intent(in) :: key_j(Nint,2), key_i(Nint,2)
double precision, intent(out) :: hji,hij
double precision :: hmono, htwoe_ji, htwoe_ij
integer :: occ(Nint*bit_kind_size,2)
integer :: Ne(2), i, j, ii, jj, ispin, jspin, k, kk
integer :: degree,exc(0:2,2,2)
integer :: h1, p1, h2, p2, s1, s2
double precision :: get_mo_two_e_integral_tc_int,phase
call get_excitation_degree(key_i, key_j, degree, Nint)
hmono = 0.d0
htwoe_ji = 0.d0
htwoe_ij = 0.d0
hji = 0.d0
hij = 0.d0
if(degree.ne.2)then
return
endif
integer :: degree_i,degree_j
call get_excitation_degree(ref_bitmask,key_i,degree_i,N_int)
call get_excitation_degree(ref_bitmask,key_j,degree_j,N_int)
call get_double_excitation(key_i, key_j, exc, phase, Nint)
call decode_exc(exc, 2, h1, p1, h2, p2, s1, s2)
if(s1.ne.s2)then
! opposite spin two-body
htwoe_ji = mo_bi_ortho_tc_two_e(p2,p1,h2,h1)
htwoe_ij = mo_bi_ortho_tc_two_e_transp(p2,p1,h2,h1)
else
! same spin two-body
! direct terms
htwoe_ji = mo_bi_ortho_tc_two_e(p2,p1,h2,h1)
htwoe_ij = mo_bi_ortho_tc_two_e_transp(p2,p1,h2,h1)
! exchange terms
htwoe_ji -= mo_bi_ortho_tc_two_e(p1,p2,h2,h1)
htwoe_ij -= mo_bi_ortho_tc_two_e_transp(p1,p2,h2,h1)
endif
htwoe_ji *= phase
hji = htwoe_ji
htwoe_ij *= phase
hij = htwoe_ij
end

142
plugins/local/slater_tc/slater_tc_opt_single.irp.f
View File

@ -618,3 +618,145 @@ subroutine get_single_excitation_from_fock_tc_no_3e(Nint, key_i, key_j, h, p, sp
end end
subroutine single_htilde_mu_mat_fock_bi_ortho_no_3e_both(Nint, key_j, key_i, hji,hij)
BEGIN_DOC
! <key_j |H_tilde | key_i> and <key_i |H_tilde | key_j> for single excitation ONLY FOR ONE- AND TWO-BODY TERMS
!!
!! WARNING !!
!
! Non hermitian !!
END_DOC
use bitmasks
implicit none
integer, intent(in) :: Nint
integer(bit_kind), intent(in) :: key_j(Nint,2), key_i(Nint,2)
double precision, intent(out) :: hji,hij
double precision :: hmono, htwoe
integer :: occ(Nint*bit_kind_size,2)
integer :: Ne(2), i, j, ii, jj, ispin, jspin, k, kk
integer :: degree,exc(0:2,2,2)
integer :: h1, p1, h2, p2, s1, s2
double precision :: get_mo_two_e_integral_tc_int, phase
double precision :: direct_int, exchange_int_12, exchange_int_23, exchange_int_13
integer :: other_spin(2)
integer(bit_kind) :: key_j_core(Nint,2), key_i_core(Nint,2)
other_spin(1) = 2
other_spin(2) = 1
hmono = 0.d0
htwoe = 0.d0
hji = 0.d0
hij = 0.d0
call get_excitation_degree(key_i, key_j, degree, Nint)
if(degree.ne.1)then
return
endif
call bitstring_to_list_ab(key_i, occ, Ne, Nint)
call get_single_excitation(key_i, key_j, exc, phase, Nint)
call decode_exc(exc,1,h1,p1,h2,p2,s1,s2)
call get_single_excitation_from_fock_tc_no_3e_both(Nint, key_i, key_j, h1, p1, s1, phase, hji,hij)
end
! ---
subroutine get_single_excitation_from_fock_tc_no_3e_both(Nint, key_i, key_j, h, p, spin, phase, hji,hij)
use bitmasks
implicit none
integer, intent(in) :: Nint
integer, intent(in) :: h, p, spin
double precision, intent(in) :: phase
integer(bit_kind), intent(in) :: key_i(Nint,2), key_j(Nint,2)
double precision, intent(out) :: hji,hij
double precision :: hmono_ji,htwoe_ji
double precision :: hmono_ij,htwoe_ij
integer(bit_kind) :: differences(Nint,2)
integer(bit_kind) :: hole(Nint,2)
integer(bit_kind) :: partcl(Nint,2)
integer :: occ_hole(Nint*bit_kind_size,2)
integer :: occ_partcl(Nint*bit_kind_size,2)
integer :: n_occ_ab_hole(2),n_occ_ab_partcl(2)
integer :: i0,i
double precision :: buffer_c_ji(mo_num), buffer_x_ji(mo_num)
double precision :: buffer_c_ij(mo_num), buffer_x_ij(mo_num)
do i = 1, mo_num
buffer_c_ji(i) = tc_2e_3idx_coulomb_integrals(i,p,h)
buffer_x_ji(i) = tc_2e_3idx_exchange_integrals(i,p,h)
buffer_c_ij(i) = tc_2e_3idx_coulomb_integrals_transp(i,p,h)
buffer_x_ij(i) = tc_2e_3idx_exchange_integrals_transp(i,p,h)
enddo
do i = 1, Nint
differences(i,1) = xor(key_i(i,1),ref_closed_shell_bitmask(i,1))
differences(i,2) = xor(key_i(i,2),ref_closed_shell_bitmask(i,2))
hole(i,1) = iand(differences(i,1),ref_closed_shell_bitmask(i,1))
hole(i,2) = iand(differences(i,2),ref_closed_shell_bitmask(i,2))
partcl(i,1) = iand(differences(i,1),key_i(i,1))
partcl(i,2) = iand(differences(i,2),key_i(i,2))
enddo
call bitstring_to_list_ab(hole, occ_hole, n_occ_ab_hole, Nint)
call bitstring_to_list_ab(partcl, occ_partcl, n_occ_ab_partcl, Nint)
hmono_ji = mo_bi_ortho_tc_one_e(p,h)
htwoe_ji = fock_op_2_e_tc_closed_shell(p,h)
hmono_ij = mo_bi_ortho_tc_one_e(h,p)
htwoe_ij = fock_op_2_e_tc_closed_shell(h,p)
! holes :: direct terms
do i0 = 1, n_occ_ab_hole(1)
i = occ_hole(i0,1)
htwoe_ji -= buffer_c_ji(i)
htwoe_ij -= buffer_c_ij(i)
enddo
do i0 = 1, n_occ_ab_hole(2)
i = occ_hole(i0,2)
htwoe_ji -= buffer_c_ji(i)
htwoe_ij -= buffer_c_ij(i)
enddo
! holes :: exchange terms
do i0 = 1, n_occ_ab_hole(spin)
i = occ_hole(i0,spin)
htwoe_ji += buffer_x_ji(i)
htwoe_ij += buffer_x_ij(i)
enddo
! particles :: direct terms
do i0 = 1, n_occ_ab_partcl(1)
i = occ_partcl(i0,1)
htwoe_ji += buffer_c_ji(i)
htwoe_ij += buffer_c_ij(i)
enddo
do i0 = 1, n_occ_ab_partcl(2)
i = occ_partcl(i0,2)
htwoe_ji += buffer_c_ji(i)
htwoe_ij += buffer_c_ij(i)
enddo
! particles :: exchange terms
do i0 = 1, n_occ_ab_partcl(spin)
i = occ_partcl(i0,spin)
htwoe_ji -= buffer_x_ji(i)
htwoe_ij -= buffer_x_ij(i)
enddo
htwoe_ji = htwoe_ji * phase
hmono_ji = hmono_ji * phase
hji = htwoe_ji + hmono_ji
htwoe_ij = htwoe_ij * phase
hmono_ij = hmono_ij * phase
hij = htwoe_ij + hmono_ij
end

140
plugins/local/slater_tc/symmetrized_3_e_int_prov.irp.f
View File

@ -1,140 +0,0 @@
BEGIN_PROVIDER [ double precision, three_e_diag_parrallel_spin_prov, (mo_num, mo_num, mo_num)]
BEGIN_DOC
!
! matrix element of the -L three-body operator ON A BI ORTHONORMAL BASIS
!
! three_e_diag_parrallel_spin_prov(m,j,i) = All combinations of the form <mji|-L|mji> for same spin matrix elements
!
! notice the -1 sign: in this way three_e_diag_parrallel_spin_prov can be directly used to compute Slater rules with a + sign
!
END_DOC
implicit none
integer :: i, j, m
double precision :: integral, wall1, wall0, three_e_diag_parrallel_spin
three_e_diag_parrallel_spin_prov = 0.d0
print *, ' Providing the three_e_diag_parrallel_spin_prov ...'
integral = three_e_diag_parrallel_spin(1,1,1) ! to provide all stuffs
call wall_time(wall0)
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,j,m,integral) &
!$OMP SHARED (mo_num,three_e_diag_parrallel_spin_prov)
!$OMP DO SCHEDULE (dynamic)
do i = 1, mo_num
do j = 1, mo_num
do m = j, mo_num
three_e_diag_parrallel_spin_prov(m,j,i) = three_e_diag_parrallel_spin(m,j,i)
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
do i = 1, mo_num
do j = 1, mo_num
do m = 1, j
three_e_diag_parrallel_spin_prov(m,j,i) = three_e_diag_parrallel_spin_prov(j,m,i)
enddo
enddo
enddo
call wall_time(wall1)
print *, ' wall time for three_e_diag_parrallel_spin_prov', wall1 - wall0
END_PROVIDER
BEGIN_PROVIDER [ double precision, three_e_single_parrallel_spin_prov, (mo_num, mo_num, mo_num, mo_num)]
BEGIN_DOC
!
! matrix element of the -L three-body operator FOR THE DIRECT TERMS OF SINGLE EXCITATIONS AND BI ORTHO MOs
!
! three_e_single_parrallel_spin_prov(m,j,k,i) = All combination of <mjk|-L|mji> for same spin matrix elements
!
! 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
integer :: i, j, k, m
double precision :: integral, wall1, wall0, three_e_single_parrallel_spin
three_e_single_parrallel_spin_prov = 0.d0
print *, ' Providing the three_e_single_parrallel_spin_prov ...'
integral = three_e_single_parrallel_spin(1,1,1,1)
call wall_time(wall0)
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,j,k,m,integral) &
!$OMP SHARED (mo_num,three_e_single_parrallel_spin_prov)
!$OMP DO SCHEDULE (dynamic)
do i = 1, mo_num
do k = 1, mo_num
do j = 1, mo_num
do m = 1, mo_num
three_e_single_parrallel_spin_prov(m,j,k,i) = three_e_single_parrallel_spin(m,j,k,i)
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
call wall_time(wall1)
print *, ' wall time for three_e_single_parrallel_spin_prov', wall1 - wall0
END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, three_e_double_parrallel_spin_prov, (mo_num, mo_num, mo_num, mo_num, mo_num)]
BEGIN_DOC
!
! matrix element of the -L three-body operator FOR THE DIRECT TERMS OF DOUBLE EXCITATIONS AND BI ORTHO MOs
!
! three_e_double_parrallel_spin_prov(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
integer :: i, j, k, m, l
double precision :: integral, wall1, wall0, three_e_double_parrallel_spin
three_e_double_parrallel_spin_prov = 0.d0
print *, ' Providing the three_e_double_parrallel_spin_prov ...'
call wall_time(wall0)
integral = three_e_double_parrallel_spin(1,1,1,1,1)
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,j,k,m,l,integral) &
!$OMP SHARED (mo_num,three_e_double_parrallel_spin_prov)
!$OMP DO SCHEDULE (dynamic)
do i = 1, mo_num
do k = 1, mo_num
do j = 1, mo_num
do l = 1, mo_num
do m = 1, mo_num
three_e_double_parrallel_spin_prov(m,l,j,k,i) = three_e_double_parrallel_spin(m,l,j,k,i)
enddo
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
call wall_time(wall1)
print *, ' wall time for three_e_double_parrallel_spin_prov', wall1 - wall0
END_PROVIDER

59
plugins/local/slater_tc_no_opt/.gitignore vendored Normal file
View File

@ -0,0 +1,59 @@
IRPF90_temp/
IRPF90_man/
build.ninja
irpf90.make
ezfio_interface.irp.f
irpf90_entities
tags
Makefile
ao_basis
ao_one_e_ints
ao_two_e_erf_ints
ao_two_e_ints
aux_quantities
becke_numerical_grid
bitmask
cis
cisd
cipsi
davidson
davidson_dressed
davidson_undressed
density_for_dft
determinants
dft_keywords
dft_utils_in_r
dft_utils_one_e
dft_utils_two_body
dressing
dummy
electrons
ezfio_files
fci
generators_cas
generators_full
hartree_fock
iterations
kohn_sham
kohn_sham_rs
mo_basis
mo_guess
mo_one_e_ints
mo_two_e_erf_ints
mo_two_e_ints
mpi
mrpt_utils
nuclei
perturbation
pseudo
psiref_cas
psiref_utils
scf_utils
selectors_cassd
selectors_full
selectors_utils
single_ref_method
slave
tools
utils
zmq

8
plugins/local/slater_tc_no_opt/NEED Normal file
View File

@ -0,0 +1,8 @@
determinants
normal_order_old
bi_ort_ints
bi_ortho_mos
tc_keywords
non_hermit_dav
dav_general_mat
tc_scf

4
plugins/local/slater_tc_no_opt/README.rst Normal file
View File

@ -0,0 +1,4 @@
================
slater_tc_no_opt
================

0
plugins/local/slater_tc/h_biortho.irp.f → plugins/local/slater_tc_no_opt/h_biortho.irp.f
View File

0
plugins/local/slater_tc/h_tc_bi_ortho_psi.irp.f → plugins/local/slater_tc_no_opt/h_tc_bi_ortho_psi.irp.f
View File

2
plugins/local/slater_tc/slater_tc_3e_slow.irp.f → plugins/local/slater_tc_no_opt/slater_tc_3e_slow.irp.f
View File

@ -1,7 +1,7 @@
! --- ! ---
subroutine diag_htilde_three_body_ints_bi_ort_slow(Nint, key_i, hthree) subroutine diag_htc_bi_orth_3e_brute(Nint, key_i, hthree)
BEGIN_DOC BEGIN_DOC
! diagonal element of htilde ONLY FOR THREE-BODY TERMS WITH BI ORTHONORMAL ORBITALS ! diagonal element of htilde ONLY FOR THREE-BODY TERMS WITH BI ORTHONORMAL ORBITALS

7
plugins/local/slater_tc_no_opt/slater_tc_no_opt.irp.f Normal file
View File

@ -0,0 +1,7 @@
program slater_tc_no_opt
implicit none
BEGIN_DOC
! TODO : Put the documentation of the program here
END_DOC
print *, 'Hello world'
end

73
plugins/local/slater_tc/slater_tc_slow.irp.f → plugins/local/slater_tc_no_opt/slater_tc_slow.irp.f
View File

@ -61,7 +61,7 @@ subroutine htilde_mu_mat_bi_ortho_slow(key_j, key_i, Nint, hmono, htwoe, hthree,
if(degree.gt.2) return if(degree.gt.2) return
if(degree == 0) then if(degree == 0) then
call diag_htilde_mu_mat_bi_ortho_slow(Nint, key_i, hmono, htwoe, htot) call diag_htc_bi_orth_2e_brute(Nint, key_i, hmono, htwoe, htot)
else if (degree == 1) then else if (degree == 1) then
call single_htilde_mu_mat_bi_ortho_slow(Nint, key_j, key_i, hmono, htwoe, htot) call single_htilde_mu_mat_bi_ortho_slow(Nint, key_j, key_i, hmono, htwoe, htot)
else if(degree == 2) then else if(degree == 2) then
@ -76,7 +76,7 @@ subroutine htilde_mu_mat_bi_ortho_slow(key_j, key_i, Nint, hmono, htwoe, hthree,
else if((degree == 1) .and. (elec_num .gt. 2) .and. three_e_4_idx_term) then else if((degree == 1) .and. (elec_num .gt. 2) .and. three_e_4_idx_term) then
call single_htilde_three_body_ints_bi_ort_slow(Nint, key_j, key_i, hthree) call single_htilde_three_body_ints_bi_ort_slow(Nint, key_j, key_i, hthree)
else if((degree == 0) .and. (elec_num .gt. 2) .and. three_e_3_idx_term) then else if((degree == 0) .and. (elec_num .gt. 2) .and. three_e_3_idx_term) then
call diag_htilde_three_body_ints_bi_ort_slow(Nint, key_i, hthree) call diag_htc_bi_orth_3e_brute(Nint, key_i, hthree)
endif endif
endif endif
@ -95,75 +95,6 @@ end
! --- ! ---
subroutine diag_htilde_mu_mat_bi_ortho_slow(Nint, key_i, hmono, htwoe, htot)
BEGIN_DOC
!
! diagonal element of htilde ONLY FOR ONE- AND TWO-BODY TERMS
!
END_DOC
use bitmasks
implicit none
integer, intent(in) :: Nint
integer(bit_kind), intent(in) :: key_i(Nint,2)
double precision, intent(out) :: hmono,htwoe,htot
integer :: occ(Nint*bit_kind_size,2)
integer :: Ne(2), i, j, ii, jj, ispin, jspin, k, kk
double precision :: get_mo_two_e_integral_tc_int
integer(bit_kind) :: key_i_core(Nint,2)
PROVIDE mo_bi_ortho_tc_two_e
hmono = 0.d0
htwoe = 0.d0
htot = 0.d0
call bitstring_to_list_ab(key_i, occ, Ne, Nint)
do ispin = 1, 2
do i = 1, Ne(ispin)
ii = occ(i,ispin)
hmono += mo_bi_ortho_tc_one_e(ii,ii)
enddo
enddo
! alpha/beta two-body
ispin = 1
jspin = 2
do i = 1, Ne(ispin) ! electron 1 (so it can be associated to mu(r1))
ii = occ(i,ispin)
do j = 1, Ne(jspin) ! electron 2
jj = occ(j,jspin)
htwoe += mo_bi_ortho_tc_two_e(jj,ii,jj,ii)
enddo
enddo
! alpha/alpha two-body
do i = 1, Ne(ispin)
ii = occ(i,ispin)
do j = i+1, Ne(ispin)
jj = occ(j,ispin)
htwoe += mo_bi_ortho_tc_two_e(ii,jj,ii,jj) - mo_bi_ortho_tc_two_e(ii,jj,jj,ii)
enddo
enddo
! beta/beta two-body
do i = 1, Ne(jspin)
ii = occ(i,jspin)
do j = i+1, Ne(jspin)
jj = occ(j,jspin)
htwoe += mo_bi_ortho_tc_two_e(ii,jj,ii,jj) - mo_bi_ortho_tc_two_e(ii,jj,jj,ii)
enddo
enddo
htot = hmono + htwoe
end
! ---
subroutine double_htilde_mu_mat_bi_ortho_slow(Nint, key_j, key_i, hmono, htwoe, htot) subroutine double_htilde_mu_mat_bi_ortho_slow(Nint, key_j, key_i, hmono, htwoe, htot)
BEGIN_DOC BEGIN_DOC

10
plugins/local/tc_bi_ortho/test_tc_bi_ortho.irp.f → plugins/local/slater_tc_no_opt/test_tc_bi_ortho.irp.f
View File

@ -88,7 +88,7 @@ subroutine test_slater_tc_opt
i_count = 0.d0 i_count = 0.d0
do i = 1, N_det do i = 1, N_det
do j = 1,N_det do j = 1,N_det
call htilde_mu_mat_bi_ortho_slow(psi_det(1,1,j), psi_det(1,1,i), N_int, hmono, htwoe, hthree, htot) call htilde_mu_mat_opt_bi_ortho(psi_det(1,1,j), psi_det(1,1,i), N_int, hmono, htwoe, hthree, htot)
call htilde_mu_mat_opt_bi_ortho(psi_det(1,1,j), psi_det(1,1,i), N_int, hnewmono, hnewtwoe, hnewthree, hnewtot) call htilde_mu_mat_opt_bi_ortho(psi_det(1,1,j), psi_det(1,1,i), N_int, hnewmono, hnewtwoe, hnewthree, hnewtot)
if(dabs(htot).gt.1.d-15)then if(dabs(htot).gt.1.d-15)then
i_count += 1.D0 i_count += 1.D0
@ -124,7 +124,7 @@ subroutine timing_tot
do j = 1, N_det do j = 1, N_det
! call get_excitation_degree(psi_det(1,1,j), psi_det(1,1,i),degree,N_int) ! call get_excitation_degree(psi_det(1,1,j), psi_det(1,1,i),degree,N_int)
i_count += 1.d0 i_count += 1.d0
call htilde_mu_mat_bi_ortho_slow(psi_det(1,1,j), psi_det(1,1,i), N_int, hmono, htwoe, hthree, htot) call htilde_mu_mat_opt_bi_ortho(psi_det(1,1,j), psi_det(1,1,i), N_int, hmono, htwoe, hthree, htot)
enddo enddo
enddo enddo
call wall_time(wall1) call wall_time(wall1)
@ -171,7 +171,7 @@ subroutine timing_diag
do i = 1, N_det do i = 1, N_det
do j = i,i do j = i,i
i_count += 1.d0 i_count += 1.d0
call htilde_mu_mat_bi_ortho_slow(psi_det(1,1,j), psi_det(1,1,i), N_int, hmono, htwoe, hthree, htot) call htilde_mu_mat_opt_bi_ortho(psi_det(1,1,j), psi_det(1,1,i), N_int, hmono, htwoe, hthree, htot)
enddo enddo
enddo enddo
call wall_time(wall1) call wall_time(wall1)
@ -208,7 +208,7 @@ subroutine timing_single
if(degree.ne.1)cycle if(degree.ne.1)cycle
i_count += 1.d0 i_count += 1.d0
call wall_time(wall0) call wall_time(wall0)
call htilde_mu_mat_bi_ortho_slow(psi_det(1,1,j), psi_det(1,1,i), N_int, hmono, htwoe, hthree, htot) call htilde_mu_mat_opt_bi_ortho(psi_det(1,1,j), psi_det(1,1,i), N_int, hmono, htwoe, hthree, htot)
call wall_time(wall1) call wall_time(wall1)
accu += wall1 - wall0 accu += wall1 - wall0
enddo enddo
@ -250,7 +250,7 @@ subroutine timing_double
if(degree.ne.2)cycle if(degree.ne.2)cycle
i_count += 1.d0 i_count += 1.d0
call wall_time(wall0) call wall_time(wall0)
call htilde_mu_mat_bi_ortho_slow(psi_det(1,1,j), psi_det(1,1,i), N_int, hmono, htwoe, hthree, htot) call htilde_mu_mat_opt_bi_ortho(psi_det(1,1,j), psi_det(1,1,i), N_int, hmono, htwoe, hthree, htot)
call wall_time(wall1) call wall_time(wall1)
accu += wall1 - wall0 accu += wall1 - wall0
enddo enddo

59
plugins/local/spher_harm/.gitignore vendored Normal file
View File

@ -0,0 +1,59 @@
IRPF90_temp/
IRPF90_man/
build.ninja
irpf90.make
ezfio_interface.irp.f
irpf90_entities
tags
Makefile
ao_basis
ao_one_e_ints
ao_two_e_erf_ints
ao_two_e_ints
aux_quantities
becke_numerical_grid
bitmask
cis
cisd
cipsi
davidson
davidson_dressed
davidson_undressed
density_for_dft
determinants
dft_keywords
dft_utils_in_r
dft_utils_one_e
dft_utils_two_body
dressing
dummy
electrons
ezfio_files
fci
generators_cas
generators_full
hartree_fock
iterations
kohn_sham
kohn_sham_rs
mo_basis
mo_guess
mo_one_e_ints
mo_two_e_erf_ints
mo_two_e_ints
mpi
mrpt_utils
nuclei
perturbation
pseudo
psiref_cas
psiref_utils
scf_utils
selectors_cassd
selectors_full
selectors_utils
single_ref_method
slave
tools
utils
zmq

1
plugins/local/spher_harm/NEED Normal file
View File

@ -0,0 +1 @@
dft_utils_in_r

7
plugins/local/spher_harm/README.rst Normal file
View File

@ -0,0 +1,7 @@
==========
spher_harm
==========
Routines for spherical Harmonics evaluation in real space.
The main routine is "spher_harm_func_r3(r,l,m,re_ylm, im_ylm)".
The test routine is "test_spher_harm" where everything is explained in details.

50
plugins/local/spher_harm/assoc_gaus_pol.irp.f Normal file
View File

@ -0,0 +1,50 @@
double precision function plgndr(l,m,x)
integer, intent(in) :: l,m
double precision, intent(in) :: x
BEGIN_DOC
! associated Legenre polynom P_l,m(x). Used for the Y_lm(theta,phi)
! Taken from https://iate.oac.uncor.edu/~mario/materia/nr/numrec/f6-8.pdf
END_DOC
integer :: i,ll
double precision :: fact,pll,pmm,pmmp1,somx2
if(m.lt.0.or.m.gt.l.or.dabs(x).gt.1.d0)then
print*,'bad arguments in plgndr'
pause
endif
pmm=1.d0
if(m.gt.0) then
somx2=dsqrt((1.d0-x)*(1.d0+x))
fact=1.d0
do i=1,m
pmm=-pmm*fact*somx2
fact=fact+2.d0
enddo
endif ! m > 0
if(l.eq.m) then
plgndr=pmm
else
pmmp1=x*(2*m+1)*pmm ! Compute P_m+1^m
if(l.eq.m+1) then
plgndr=pmmp1
else ! Compute P_l^m, l> m+1
do ll=m+2,l
pll=(x*dble(2*ll-1)*pmmp1-dble(ll+m-1)*pmm)/(ll-m)
pmm=pmmp1
pmmp1=pll
enddo
plgndr=pll
endif ! l.eq.m+1
endif ! l.eq.m
return
end
double precision function ortho_assoc_gaus_pol(l1,m1,l2)
implicit none
integer, intent(in) :: l1,m1,l2
double precision :: fact
if(l1.ne.l2)then
ortho_assoc_gaus_pol= 0.d0
else
ortho_assoc_gaus_pol = 2.d0*fact(l1+m1) / (dble(2*l1+1)*fact(l1-m1))
endif
end

231
plugins/local/spher_harm/routines_test.irp.f Normal file
View File

@ -0,0 +1,231 @@
subroutine test_spher_harm
implicit none
BEGIN_DOC
! routine to test the generic spherical harmonics routine "spher_harm_func_r3" from R^3 --> C
!
! We test <Y_l1,m1|Y_l2,m2> = delta_m1,m2 delta_l1,l2
!
! The test is done through the integration on a sphere with the Lebedev grid.
END_DOC
include 'constants.include.F'
integer :: l1,m1,i,l2,m2,lmax
double precision :: r(3),weight,accu_re, accu_im,accu
double precision :: re_ylm_1, im_ylm_1,re_ylm_2, im_ylm_2
double precision :: theta,phi,r_abs
lmax = 5 ! Maximum angular momentum until which we are going to test orthogonality conditions
do l1 = 0,lmax
do m1 = -l1 ,l1
do l2 = 0,lmax
do m2 = -l2 ,l2
accu_re = 0.d0 ! accumulator for the REAL part of <Y_l1,m1|Y_l2,m2>
accu_im = 0.d0 ! accumulator for the IMAGINARY part of <Y_l1,m1|Y_l2,m2>
accu = 0.d0 ! accumulator for the weights ==> should be \int dOmega == 4 pi
! <l1,m1|l2,m2> = \int dOmega Y_l1,m1^* Y_l2,m2
! \approx \sum_i W_i Y_l1,m1^*(r_i) Y_l2,m2(r_i) WITH r_i being on the spher of radius 1
do i = 1, n_points_integration_angular
r(1:3) = angular_quadrature_points(i,1:3) ! ith Lebedev point (x,y,z) on the sphere of radius 1
weight = weights_angular_points(i) ! associated Lebdev weight not necessarily positive
!!!!!!!!!!! Test of the Cartesian --> Spherical coordinates
! theta MUST belong to [0,pi] and phi to [0,2pi]
! gets the cartesian to spherical change of coordinates
call cartesian_to_spherical(r,theta,phi,r_abs)
if(theta.gt.pi.or.theta.lt.0.d0)then
print*,'pb with theta, it should be in [0,pi]',theta
print*,r
endif
if(phi.gt.2.d0*pi.or.phi.lt.0.d0)then
print*,'pb with phi, it should be in [0,2 pi]',phi/pi
print*,r
endif
!!!!!!!!!!! Routines returning the Spherical harmonics on the grid point
call spher_harm_func_r3(r,l1,m1,re_ylm_1, im_ylm_1)
call spher_harm_func_r3(r,l2,m2,re_ylm_2, im_ylm_2)
!!!!!!!!!!! Integration of Y_l1,m1^*(r) Y_l2,m2(r)
! = \int dOmega (re_ylm_1 -i im_ylm_1) * (re_ylm_2 +i im_ylm_2)
! = \int dOmega (re_ylm_1*re_ylm_2 + im_ylm_1*im_ylm_2) +i (im_ylm_2*re_ylm_1 - im_ylm_1*re_ylm_2)
accu_re += weight * (re_ylm_1*re_ylm_2 + im_ylm_1*im_ylm_2)
accu_im += weight * (im_ylm_2*re_ylm_1 - im_ylm_1*re_ylm_2)
accu += weight
enddo
! Test that the sum of the weights is 4 pi
if(dabs(accu - dfour_pi).gt.1.d-6)then
print*,'Problem !! The sum of the Lebedev weight is not 4 pi ..'
print*,accu
stop
endif
! Test for the delta l1,l2 and delta m1,m2
!
! Test for the off-diagonal part of the Kronecker delta
if(l1.ne.l2.or.m1.ne.m2)then
if(dabs(accu_re).gt.1.d-6.or.dabs(accu_im).gt.1.d-6)then
print*,'pb OFF DIAG !!!!! '
print*,'l1,m1,l2,m2',l1,m1,l2,m2
print*,'accu_re = ',accu_re
print*,'accu_im = ',accu_im
endif
endif
! Test for the diagonal part of the Kronecker delta
if(l1==l2.and.m1==m2)then
if(dabs(accu_re-1.d0).gt.1.d-5.or.dabs(accu_im).gt.1.d-6)then
print*,'pb DIAG !!!!! '
print*,'l1,m1,l2,m2',l1,m1,l2,m2
print*,'accu_re = ',accu_re
print*,'accu_im = ',accu_im
endif
endif
enddo
enddo
enddo
enddo
end
subroutine test_cart
implicit none
BEGIN_DOC
! test for the cartesian --> spherical change of coordinates
!
! test the routine "cartesian_to_spherical" such that the polar angle theta ranges in [0,pi]
!
! and the asymuthal angle phi ranges in [0,2pi]
END_DOC
include 'constants.include.F'
double precision :: r(3),theta,phi,r_abs
print*,''
r = 0.d0
r(1) = 1.d0
r(2) = 1.d0
call cartesian_to_spherical(r,theta,phi,r_abs)
print*,r
print*,phi/pi
print*,''
r = 0.d0
r(1) =-1.d0
r(2) = 1.d0
call cartesian_to_spherical(r,theta,phi,r_abs)
print*,r
print*,phi/pi
print*,''
r = 0.d0
r(1) =-1.d0
r(2) =-1.d0
call cartesian_to_spherical(r,theta,phi,r_abs)
print*,r
print*,phi/pi
print*,''
r = 0.d0
r(1) = 1.d0
r(2) =-1.d0
call cartesian_to_spherical(r,theta,phi,r_abs)
print*,r
print*,phi/pi
end
subroutine test_brutal_spheric
implicit none
include 'constants.include.F'
BEGIN_DOC
! Test for the <Y_l1,m1|Y_l2,m2> = delta_m1,m2 delta_l1,l2 using the following two dimentional integration
!
! \int_0^2pi d Phi \int_-1^+1 d(cos(Theta)) Y_l1,m1^*(Theta,Phi) Y_l2,m2(Theta,Phi)
!
!= \int_0^2pi d Phi \int_0^pi dTheta sin(Theta) Y_l1,m1^*(Theta,Phi) Y_l2,m2(Theta,Phi)
!
! Allows to test for the general functions "spher_harm_func_m_pos" with "spher_harm_func_expl"
END_DOC
integer :: itheta, iphi,ntheta,nphi
double precision :: theta_min, theta_max, dtheta,theta
double precision :: phi_min, phi_max, dphi,phi
double precision :: accu_re, accu_im,weight
double precision :: re_ylm_1, im_ylm_1 ,re_ylm_2, im_ylm_2,accu
integer :: l1,m1,i,l2,m2,lmax
phi_min = 0.d0
phi_max = 2.D0 * pi
theta_min = 0.d0
theta_max = 1.D0 * pi
ntheta = 1000
nphi = 1000
dphi = (phi_max - phi_min)/dble(nphi)
dtheta = (theta_max - theta_min)/dble(ntheta)
lmax = 2
do l1 = 0,lmax
do m1 = 0 ,l1
do l2 = 0,lmax
do m2 = 0 ,l2
accu_re = 0.d0
accu_im = 0.d0
accu = 0.d0
theta = theta_min
do itheta = 1, ntheta
phi = phi_min
do iphi = 1, nphi
! call spher_harm_func_expl(l1,m1,theta,phi,re_ylm_1, im_ylm_1)
! call spher_harm_func_expl(l2,m2,theta,phi,re_ylm_2, im_ylm_2)
call spher_harm_func_m_pos(l1,m1,theta,phi,re_ylm_1, im_ylm_1)
call spher_harm_func_m_pos(l2,m2,theta,phi,re_ylm_2, im_ylm_2)
weight = dtheta * dphi * dsin(theta)
accu_re += weight * (re_ylm_1*re_ylm_2 + im_ylm_1*im_ylm_2)
accu_im += weight * (im_ylm_2*re_ylm_1 - im_ylm_1*re_ylm_2)
accu += weight
phi += dphi
enddo
theta += dtheta
enddo
print*,'l1,m1,l2,m2',l1,m1,l2,m2
print*,'accu_re = ',accu_re
print*,'accu_im = ',accu_im
print*,'accu = ',accu
if(l1.ne.l2.or.m1.ne.m2)then
if(dabs(accu_re).gt.1.d-6.or.dabs(accu_im).gt.1.d-6)then
print*,'pb OFF DIAG !!!!! '
endif
endif
if(l1==l2.and.m1==m2)then
if(dabs(accu_re-1.d0).gt.1.d-5.or.dabs(accu_im).gt.1.d-6)then
print*,'pb DIAG !!!!! '
endif
endif
enddo
enddo
enddo
enddo
end
subroutine test_assoc_leg_pol
implicit none
BEGIN_DOC
! Test for the associated Legendre Polynoms. The test is done through the orthogonality condition.
END_DOC
print *, 'Hello world'
integer :: l1,m1,ngrid,i,l2,m2
l1 = 0
m1 = 0
l2 = 2
m2 = 0
double precision :: x, dx,xmax,accu,xmin
double precision :: plgndr,func_1,func_2,ortho_assoc_gaus_pol
ngrid = 100000
xmax = 1.d0
xmin = -1.d0
dx = (xmax-xmin)/dble(ngrid)
do l2 = 0,10
x = xmin
accu = 0.d0
do i = 1, ngrid
func_1 = plgndr(l1,m1,x)
func_2 = plgndr(l2,m2,x)
write(33,*)x, func_1,func_2
accu += func_1 * func_2 * dx
x += dx
enddo
print*,'l2 = ',l2
print*,'accu = ',accu
print*,ortho_assoc_gaus_pol(l1,m1,l2)
enddo
end

7
plugins/local/spher_harm/spher_harm.irp.f Normal file
View File

@ -0,0 +1,7 @@
program spher_harm
implicit none
! call test_spher_harm
! call test_cart
call test_brutal_spheric
end

151
plugins/local/spher_harm/spher_harm_func.irp.f Normal file
View File

@ -0,0 +1,151 @@
subroutine spher_harm_func_r3(r,l,m,re_ylm, im_ylm)
implicit none
integer, intent(in) :: l,m
double precision, intent(in) :: r(3)
double precision, intent(out) :: re_ylm, im_ylm
double precision :: theta, phi,r_abs
call cartesian_to_spherical(r,theta,phi,r_abs)
call spher_harm_func(l,m,theta,phi,re_ylm, im_ylm)
end
subroutine spher_harm_func_m_pos(l,m,theta,phi,re_ylm, im_ylm)
include 'constants.include.F'
implicit none
BEGIN_DOC
! Y_lm(theta,phi) with m >0
!
END_DOC
double precision, intent(in) :: theta, phi
integer, intent(in) :: l,m
double precision, intent(out):: re_ylm,im_ylm
double precision :: prefact,fact,cos_theta,plgndr,p_lm
double precision :: tmp
prefact = dble(2*l+1)*fact(l-m)/(dfour_pi * fact(l+m))
prefact = dsqrt(prefact)
cos_theta = dcos(theta)
p_lm = plgndr(l,m,cos_theta)
tmp = prefact * p_lm
re_ylm = dcos(dble(m)*phi) * tmp
im_ylm = dsin(dble(m)*phi) * tmp
end
subroutine spher_harm_func(l,m,theta,phi,re_ylm, im_ylm)
implicit none
BEGIN_DOC
! Y_lm(theta,phi) with -l<m<+l
!
END_DOC
double precision, intent(in) :: theta, phi
integer, intent(in) :: l,m
double precision, intent(out):: re_ylm,im_ylm
double precision :: re_ylm_pos,im_ylm_pos,tmp
integer :: minus_m
if(abs(m).gt.l)then
print*,'|m| > l in spher_harm_func !! stopping ...'
stop
endif
if(m.ge.0)then
call spher_harm_func_m_pos(l,m,theta,phi,re_ylm_pos, im_ylm_pos)
re_ylm = re_ylm_pos
im_ylm = im_ylm_pos
else
minus_m = -m !> 0
call spher_harm_func_m_pos(l,minus_m,theta,phi,re_ylm_pos, im_ylm_pos)
tmp = (-1)**minus_m
re_ylm = tmp * re_ylm_pos
im_ylm = -tmp * im_ylm_pos ! complex conjugate
endif
end
subroutine cartesian_to_spherical(r,theta,phi,r_abs)
implicit none
double precision, intent(in) :: r(3)
double precision, intent(out):: theta, phi,r_abs
double precision :: r_2,x_2_y_2,tmp
include 'constants.include.F'
x_2_y_2 = r(1)*r(1) + r(2)*r(2)
r_2 = x_2_y_2 + r(3)*r(3)
r_abs = dsqrt(r_2)
if(r_abs.gt.1.d-20)then
theta = dacos(r(3)/r_abs)
else
theta = 0.d0
endif
if(.true.)then
if(dabs(r(1)).gt.0.d0)then
tmp = datan(r(2)/r(1))
! phi = datan2(r(2),r(1))
endif
! From Wikipedia on Spherical Harmonics
if(r(1).gt.0.d0)then
phi = tmp
else if(r(1).lt.0.d0.and.r(2).ge.0.d0)then
phi = tmp + pi
else if(r(1).lt.0.d0.and.r(2).lt.0.d0)then
phi = tmp - pi
else if(r(1)==0.d0.and.r(2).gt.0.d0)then
phi = 0.5d0*pi
else if(r(1)==0.d0.and.r(2).lt.0.d0)then
phi =-0.5d0*pi
else if(r(1)==0.d0.and.r(2)==0.d0)then
phi = 0.d0
endif
if(r(2).lt.0.d0.and.r(1).le.0.d0)then
tmp = pi - dabs(phi)
phi = pi + tmp
else if(r(2).lt.0.d0.and.r(1).gt.0.d0)then
phi = dtwo_pi + phi
endif
endif
if(.false.)then
x_2_y_2 = dsqrt(x_2_y_2)
if(dabs(x_2_y_2).gt.1.d-20.and.dabs(r(2)).gt.1.d-20)then
phi = dabs(r(2))/r(2) * dacos(r(1)/x_2_y_2)
else
phi = 0.d0
endif
endif
end
subroutine spher_harm_func_expl(l,m,theta,phi,re_ylm, im_ylm)
implicit none
BEGIN_DOC
! Y_lm(theta,phi) with -l<m<+l and 0<= l <=2
!
END_DOC
double precision, intent(in) :: theta, phi
integer, intent(in) :: l,m
double precision, intent(out):: re_ylm,im_ylm
double precision :: tmp
include 'constants.include.F'
if(l==0.and.m==0)then
re_ylm = 0.5d0 * inv_sq_pi
im_ylm = 0.d0
else if(l==1.and.m==1)then
tmp = - inv_sq_pi * dsqrt(3.d0/8.d0) * dsin(theta)
re_ylm = tmp * dcos(phi)
im_ylm = tmp * dsin(phi)
else if(l==1.and.m==0)then
tmp = inv_sq_pi * dsqrt(3.d0/4.d0) * dcos(theta)
re_ylm = tmp
im_ylm = 0.d0
else if(l==2.and.m==2)then
tmp = 0.25d0 * inv_sq_pi * dsqrt(0.5d0*15.d0) * dsin(theta)*dsin(theta)
re_ylm = tmp * dcos(2.d0*phi)
im_ylm = tmp * dsin(2.d0*phi)
else if(l==2.and.m==1)then
tmp = - inv_sq_pi * dsqrt(15.d0/8.d0) * dsin(theta) * dcos(theta)
re_ylm = tmp * dcos(phi)
im_ylm = tmp * dsin(phi)
else if(l==2.and.m==0)then
tmp = dsqrt(5.d0/4.d0) * inv_sq_pi* (1.5d0*dcos(theta)*dcos(theta)-0.5d0)
re_ylm = tmp
im_ylm = 0.d0
endif
end

34
plugins/local/tc_bi_ortho/tc_bi_ortho.irp.f → plugins/local/tc_bi_ortho/diagonalize_tc_h.irp.f
View File

@ -35,8 +35,8 @@ program tc_bi_ortho
print*, ' nb of det = ', N_det print*, ' nb of det = ', N_det
call routine_diag() call routine_diag()
call write_tc_energy() ! call write_tc_energy()
call save_tc_bi_ortho_wavefunction() ! call save_tc_bi_ortho_wavefunction()
end end
@ -76,28 +76,26 @@ subroutine routine_diag()
PROVIDE noL_2e PROVIDE noL_2e
endif endif
PROVIDE htilde_matrix_elmt_bi_ortho
return
if(N_states .eq. 1) then if(N_states .eq. 1) then
print*,'eigval_right_tc_bi_orth = ',eigval_right_tc_bi_orth(1) print*,'eigval_right_tc_bi_orth = ',eigval_right_tc_bi_orth(1)
print*,'e_tc_left_right = ',e_tc_left_right ! print*,'e_tc_left_right = ',e_tc_left_right
print*,'e_tilde_bi_orth_00 = ',e_tilde_bi_orth_00 ! print*,'e_tilde_bi_orth_00 = ',e_tilde_bi_orth_00
print*,'e_pt2_tc_bi_orth = ',e_pt2_tc_bi_orth ! print*,'e_pt2_tc_bi_orth = ',e_pt2_tc_bi_orth
print*,'e_pt2_tc_bi_orth_single = ',e_pt2_tc_bi_orth_single ! print*,'e_pt2_tc_bi_orth_single = ',e_pt2_tc_bi_orth_single
print*,'e_pt2_tc_bi_orth_double = ',e_pt2_tc_bi_orth_double ! print*,'e_pt2_tc_bi_orth_double = ',e_pt2_tc_bi_orth_double
print*,'***' ! print*,'***'
print*,'e_corr_bi_orth = ',e_corr_bi_orth ! print*,'e_corr_bi_orth = ',e_corr_bi_orth
print*,'e_corr_bi_orth_proj = ',e_corr_bi_orth_proj ! print*,'e_corr_bi_orth_proj = ',e_corr_bi_orth_proj
print*,'e_corr_bi_orth_proj_abs = ',e_corr_bi_orth_proj_abs ! print*,'e_corr_bi_orth_proj_abs = ',e_corr_bi_orth_proj_abs
print*,'e_corr_single_bi_orth = ',e_corr_single_bi_orth ! print*,'e_corr_single_bi_orth = ',e_corr_single_bi_orth
print*,'e_corr_double_bi_orth = ',e_corr_double_bi_orth ! print*,'e_corr_double_bi_orth = ',e_corr_double_bi_orth
print*,'e_corr_single_bi_orth_abs = ',e_corr_single_bi_orth_abs ! print*,'e_corr_single_bi_orth_abs = ',e_corr_single_bi_orth_abs
print*,'e_corr_double_bi_orth_abs = ',e_corr_double_bi_orth_abs ! print*,'e_corr_double_bi_orth_abs = ',e_corr_double_bi_orth_abs
print*,'Left/right eigenvectors' print*,'Left/right eigenvectors'
do i = 1,N_det do i = 1,N_det
write(*,'(I5,X,(100(F12.7,X)))')i,leigvec_tc_bi_orth(i,1),reigvec_tc_bi_orth(i,1),leigvec_tc_bi_orth(i,1)*reigvec_tc_bi_orth(i,1) write(*,'(I6,X,(100(F12.7,X)))')i,leigvec_tc_bi_orth(i,1),reigvec_tc_bi_orth(i,1),leigvec_tc_bi_orth(i,1)*reigvec_tc_bi_orth(i,1)
enddo enddo
else else

19
plugins/local/tc_bi_ortho/e_corr_bi_ortho.irp.f
View File

@ -2,7 +2,7 @@
BEGIN_PROVIDER [ double precision, e_tilde_00] BEGIN_PROVIDER [ double precision, e_tilde_00]
implicit none implicit none
double precision :: hmono,htwoe,hthree,htot double precision :: hmono,htwoe,hthree,htot
call htilde_mu_mat_bi_ortho_slow(HF_bitmask,HF_bitmask,N_int,hmono,htwoe,hthree,htot) call htilde_mu_mat_opt_bi_ortho(HF_bitmask,HF_bitmask,N_int,hmono,htwoe,hthree,htot)
e_tilde_00 = htot e_tilde_00 = htot
END_PROVIDER END_PROVIDER
@ -18,16 +18,15 @@
do i = 1, N_det do i = 1, N_det
call get_excitation_degree(HF_bitmask,psi_det(1,1,i),degree,N_int) call get_excitation_degree(HF_bitmask,psi_det(1,1,i),degree,N_int)
if(degree == 1 .or. degree == 2)then if(degree == 1 .or. degree == 2)then
call htilde_mu_mat_bi_ortho_slow(psi_det(1,1,i),HF_bitmask,N_int,hmono,htwoe,hthree,htilde_ij) call htilde_mu_mat_opt_bi_ortho(psi_det(1,1,i),HF_bitmask,N_int,hmono,htwoe,hthree,htilde_ij)
call htilde_mu_mat_bi_ortho_slow(psi_det(1,1,i),psi_det(1,1,i),N_int,hmono,htwoe,hthree,e_i0) call htilde_mu_mat_opt_bi_ortho(psi_det(1,1,i),psi_det(1,1,i),N_int,hmono,htwoe,hthree,e_i0)
delta_e = e_tilde_00 - e_i0 delta_e = e_tilde_00 - e_i0
coef_pt1 = htilde_ij / delta_e coef_pt1 = htilde_ij / delta_e
call htilde_mu_mat_bi_ortho_slow(HF_bitmask,psi_det(1,1,i),N_int,hmono,htwoe,hthree,htilde_ij) call htilde_mu_mat_opt_bi_ortho(HF_bitmask,psi_det(1,1,i),N_int,hmono,htwoe,hthree,htilde_ij)
e_pt2_tc_bi_orth += coef_pt1 * htilde_ij e_pt2_tc_bi_orth += coef_pt1 * htilde_ij
if(degree == 1)then if(degree == 1)then
e_pt2_tc_bi_orth_single += coef_pt1 * htilde_ij e_pt2_tc_bi_orth_single += coef_pt1 * htilde_ij
else else
! print*,'coef_pt1, e_pt2',coef_pt1,coef_pt1 * htilde_ij
e_pt2_tc_bi_orth_double += coef_pt1 * htilde_ij e_pt2_tc_bi_orth_double += coef_pt1 * htilde_ij
endif endif
endif endif
@ -37,7 +36,7 @@
BEGIN_PROVIDER [ double precision, e_tilde_bi_orth_00] BEGIN_PROVIDER [ double precision, e_tilde_bi_orth_00]
implicit none implicit none
double precision :: hmono,htwoe,hthree,htilde_ij double precision :: hmono,htwoe,hthree,htilde_ij
call htilde_mu_mat_bi_ortho_slow(HF_bitmask,HF_bitmask,N_int,hmono,htwoe,hthree,e_tilde_bi_orth_00) call htilde_mu_mat_opt_bi_ortho(HF_bitmask,HF_bitmask,N_int,hmono,htwoe,hthree,e_tilde_bi_orth_00)
e_tilde_bi_orth_00 += nuclear_repulsion e_tilde_bi_orth_00 += nuclear_repulsion
END_PROVIDER END_PROVIDER
@ -57,7 +56,7 @@
e_corr_double_bi_orth = 0.d0 e_corr_double_bi_orth = 0.d0
do i = 1, N_det do i = 1, N_det
call get_excitation_degree(HF_bitmask,psi_det(1,1,i),degree,N_int) call get_excitation_degree(HF_bitmask,psi_det(1,1,i),degree,N_int)
call htilde_mu_mat_bi_ortho_slow(HF_bitmask,psi_det(1,1,i),N_int,hmono,htwoe,hthree,htilde_ij) call htilde_mu_mat_opt_bi_ortho(HF_bitmask,psi_det(1,1,i),N_int,hmono,htwoe,hthree,htilde_ij)
if(degree == 1)then if(degree == 1)then
e_corr_single_bi_orth += reigvec_tc_bi_orth(i,1) * htilde_ij/reigvec_tc_bi_orth(1,1) e_corr_single_bi_orth += reigvec_tc_bi_orth(i,1) * htilde_ij/reigvec_tc_bi_orth(1,1)
e_corr_single_bi_orth_abs += dabs(reigvec_tc_bi_orth(i,1) * htilde_ij/reigvec_tc_bi_orth(1,1)) e_corr_single_bi_orth_abs += dabs(reigvec_tc_bi_orth(i,1) * htilde_ij/reigvec_tc_bi_orth(1,1))
@ -80,7 +79,7 @@
do i = 1, N_det do i = 1, N_det
accu += reigvec_tc_bi_orth(i,1) * leigvec_tc_bi_orth(i,1) accu += reigvec_tc_bi_orth(i,1) * leigvec_tc_bi_orth(i,1)
do j = 1, N_det do j = 1, N_det
call htilde_mu_mat_bi_ortho_slow(psi_det(1,1,j),psi_det(1,1,i),N_int,hmono,htwoe,hthree,htilde_ij) call htilde_mu_mat_opt_bi_ortho(psi_det(1,1,j),psi_det(1,1,i),N_int,hmono,htwoe,hthree,htilde_ij)
e_tc_left_right += htilde_ij * reigvec_tc_bi_orth(i,1) * leigvec_tc_bi_orth(j,1) e_tc_left_right += htilde_ij * reigvec_tc_bi_orth(i,1) * leigvec_tc_bi_orth(j,1)
enddo enddo
enddo enddo
@ -99,8 +98,8 @@ BEGIN_PROVIDER [ double precision, coef_pt1_bi_ortho, (N_det)]
if(degree==0)then if(degree==0)then
coef_pt1_bi_ortho(i) = 1.d0 coef_pt1_bi_ortho(i) = 1.d0
else else
call htilde_mu_mat_bi_ortho_slow(psi_det(1,1,i),HF_bitmask,N_int,hmono,htwoe,hthree,htilde_ij) call htilde_mu_mat_opt_bi_ortho(psi_det(1,1,i),HF_bitmask,N_int,hmono,htwoe,hthree,htilde_ij)
call htilde_mu_mat_bi_ortho_slow(psi_det(1,1,i),psi_det(1,1,i),N_int,hmono,htwoe,hthree,e_i0) call htilde_mu_mat_opt_bi_ortho(psi_det(1,1,i),psi_det(1,1,i),N_int,hmono,htwoe,hthree,e_i0)
delta_e = e_tilde_00 - e_i0 delta_e = e_tilde_00 - e_i0
coef_pt1 = htilde_ij / delta_e coef_pt1 = htilde_ij / delta_e
coef_pt1_bi_ortho(i)= coef_pt1 coef_pt1_bi_ortho(i)= coef_pt1

6
plugins/local/tc_bi_ortho/psi_r_l_prov.irp.f
View File

@ -1,6 +1,7 @@
use bitmasks use bitmasks
BEGIN_PROVIDER [ double precision, psi_l_coef_bi_ortho, (psi_det_size,N_states) ] !BEGIN_PROVIDER [ double precision, psi_l_coef_bi_ortho, (psi_det_size,N_states) ]
BEGIN_PROVIDER [ double precision, psi_l_coef_bi_ortho, (N_det,N_states) ]
implicit none implicit none
BEGIN_DOC BEGIN_DOC
! The wave function coefficients. Initialized with Hartree-Fock if the |EZFIO| file ! The wave function coefficients. Initialized with Hartree-Fock if the |EZFIO| file
@ -68,7 +69,8 @@ BEGIN_PROVIDER [ double precision, psi_l_coef_bi_ortho, (psi_det_size,N_states)
END_PROVIDER END_PROVIDER
BEGIN_PROVIDER [ double precision, psi_r_coef_bi_ortho, (psi_det_size,N_states) ] !BEGIN_PROVIDER [ double precision, psi_r_coef_bi_ortho, (psi_det_size,N_states) ]
BEGIN_PROVIDER [ double precision, psi_r_coef_bi_ortho, (N_det,N_states) ]
implicit none implicit none
BEGIN_DOC BEGIN_DOC
! The wave function coefficients. Initialized with Hartree-Fock if the |EZFIO| file ! The wave function coefficients. Initialized with Hartree-Fock if the |EZFIO| file

129
plugins/local/tc_bi_ortho/pt2_tc_cisd.irp.f
View File

@ -1,129 +0,0 @@
program pt2_tc_cisd
BEGIN_DOC
!
! TODO : Reads psi_det in the EZFIO folder and prints out the left- and right-eigenvectors together
! with the energy. Saves the left-right wave functions at the end.
!
END_DOC
implicit none
my_grid_becke = .True.
PROVIDE tc_grid1_a tc_grid1_r
my_n_pt_r_grid = tc_grid1_r
my_n_pt_a_grid = tc_grid1_a
touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
read_wf = .True.
touch read_wf
print*, ' nb of states = ', N_states
print*, ' nb of det = ', N_det
call routine_diag()
call routine
end
subroutine routine
implicit none
integer :: i,h1,p1,h2,p2,s1,s2,degree
double precision :: h0i,hi0,e00,ei,delta_e
double precision :: norm,e_corr,coef,e_corr_pos,e_corr_neg,e_corr_abs
integer :: exc(0:2,2,2)
double precision :: phase
double precision :: eh1,ep1,eh2,ep2
norm = 0.d0
e_corr = 0.d0
e_corr_abs = 0.d0
e_corr_pos = 0.d0
e_corr_neg = 0.d0
call htilde_mu_mat_bi_ortho_tot_slow(psi_det(1,1,1), psi_det(1,1,1), N_int, e00)
do i = 2, N_det
call htilde_mu_mat_bi_ortho_tot_slow(psi_det(1,1,i), psi_det(1,1,1), N_int, hi0)
call htilde_mu_mat_bi_ortho_tot_slow(psi_det(1,1,1), psi_det(1,1,i), N_int, h0i)
call htilde_mu_mat_bi_ortho_tot_slow(psi_det(1,1,i), psi_det(1,1,i), N_int, ei)
call get_excitation_degree(psi_det(1,1,1), psi_det(1,1,i),degree,N_int)
call get_excitation(psi_det(1,1,1), psi_det(1,1,i),exc,degree,phase,N_int)
call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
eh1 = Fock_matrix_tc_diag_mo_tot(h1)
ep1 = Fock_matrix_tc_diag_mo_tot(p1)
delta_e = eh1 - ep1
if (degree==2)then
eh2 = Fock_matrix_tc_diag_mo_tot(h2)
ep2 = Fock_matrix_tc_diag_mo_tot(p2)
delta_e += eh2 - ep2
endif
! delta_e = e00 - ei
coef = hi0/delta_e
norm += coef*coef
e_corr = coef* h0i
if(e_corr.lt.0.d0)then
e_corr_neg += e_corr
elseif(e_corr.gt.0.d0)then
e_corr_pos += e_corr
endif
e_corr_abs += dabs(e_corr)
enddo
print*,'e_corr_abs = ',e_corr_abs
print*,'e_corr_pos = ',e_corr_pos
print*,'e_corr_neg = ',e_corr_neg
print*,'norm = ',dsqrt(norm)
end
subroutine routine_diag()
implicit none
integer :: i, j, k
double precision :: dE
! provide eigval_right_tc_bi_orth
! provide overlap_bi_ortho
! provide htilde_matrix_elmt_bi_ortho
if(N_states .eq. 1) then
print*,'eigval_right_tc_bi_orth = ',eigval_right_tc_bi_orth(1)
print*,'e_tc_left_right = ',e_tc_left_right
print*,'e_tilde_bi_orth_00 = ',e_tilde_bi_orth_00
print*,'e_pt2_tc_bi_orth = ',e_pt2_tc_bi_orth
print*,'e_pt2_tc_bi_orth_single = ',e_pt2_tc_bi_orth_single
print*,'e_pt2_tc_bi_orth_double = ',e_pt2_tc_bi_orth_double
print*,'***'
print*,'e_corr_bi_orth = ',e_corr_bi_orth
print*,'e_corr_bi_orth_proj = ',e_corr_bi_orth_proj
print*,'e_corr_bi_orth_proj_abs = ',e_corr_bi_orth_proj_abs
print*,'e_corr_single_bi_orth = ',e_corr_single_bi_orth
print*,'e_corr_double_bi_orth = ',e_corr_double_bi_orth
print*,'e_corr_single_bi_orth_abs = ',e_corr_single_bi_orth_abs
print*,'e_corr_double_bi_orth_abs = ',e_corr_double_bi_orth_abs
print*,'Left/right eigenvectors'
do i = 1,N_det
write(*,'(I5,X,(100(F12.7,X)))')i,leigvec_tc_bi_orth(i,1),reigvec_tc_bi_orth(i,1),leigvec_tc_bi_orth(i,1)*reigvec_tc_bi_orth(i,1)
enddo
else
print*,'eigval_right_tc_bi_orth : '
do i = 1, N_states
print*, i, eigval_right_tc_bi_orth(i)
enddo
print*,''
print*,'******************************************************'
print*,'TC Excitation energies (au) (eV)'
do i = 2, N_states
dE = eigval_right_tc_bi_orth(i) - eigval_right_tc_bi_orth(1)
print*, i, dE, dE/0.0367502d0
enddo
print*,''
endif
end

36
plugins/local/tc_bi_ortho/tc_cisd_sc2.irp.f
View File

@ -1,36 +0,0 @@
! ---
program tc_cisd_sc2
BEGIN_DOC
! TODO : Put the documentation of the program here
END_DOC
implicit none
print *, 'Hello world'
my_grid_becke = .True.
PROVIDE tc_grid1_a tc_grid1_r
my_n_pt_r_grid = tc_grid1_r
my_n_pt_a_grid = tc_grid1_a
touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
read_wf = .True.
touch read_wf
call test
end
! ---
subroutine test()
implicit none
! double precision, allocatable :: dressing_dets(:),e_corr_dets(:)
! allocate(dressing_dets(N_det),e_corr_dets(N_det))
! e_corr_dets = 0.d0
! call get_cisd_sc2_dressing(psi_det,e_corr_dets,N_det,dressing_dets)
provide eigval_tc_cisd_sc2_bi_ortho
end

145
plugins/local/tc_bi_ortho/tc_cisd_sc2_utils.irp.f
View File

@ -1,145 +0,0 @@
BEGIN_PROVIDER [ double precision, reigvec_tc_cisd_sc2_bi_ortho, (N_det,N_states)]
&BEGIN_PROVIDER [ double precision, leigvec_tc_cisd_sc2_bi_ortho, (N_det,N_states)]
&BEGIN_PROVIDER [ double precision, eigval_tc_cisd_sc2_bi_ortho, (N_states)]
implicit none
integer :: it,n_real,degree,i,istate
double precision :: e_before, e_current,thr, hmono,htwoe,hthree,accu
double precision, allocatable :: e_corr_dets(:),h0j(:), h_sc2(:,:), dressing_dets(:)
double precision, allocatable :: leigvec_tc_bi_orth_tmp(:,:),reigvec_tc_bi_orth_tmp(:,:),eigval_right_tmp(:)
allocate(leigvec_tc_bi_orth_tmp(N_det,N_det),reigvec_tc_bi_orth_tmp(N_det,N_det),eigval_right_tmp(N_det))
allocate(e_corr_dets(N_det),h0j(N_det),h_sc2(N_det,N_det),dressing_dets(N_det))
allocate(H_jj(N_det),vec_tmp(N_det,n_states_diag),eigval_tmp(N_states))
dressing_dets = 0.d0
do i = 1, N_det
call htilde_mu_mat_bi_ortho_tot_slow(psi_det(1,1,i), psi_det(1,1,i), N_int, H_jj(i))
call get_excitation_degree(HF_bitmask,psi_det(1,1,i),degree,N_int)
if(degree == 1 .or. degree == 2)then
call htilde_mu_mat_bi_ortho_slow(HF_bitmask,psi_det(1,1,i),N_int,hmono,htwoe,hthree,h0j(i))
endif
enddo
reigvec_tc_bi_orth_tmp = 0.d0
do i = 1, N_det
reigvec_tc_bi_orth_tmp(i,1) = psi_r_coef_bi_ortho(i,1)
enddo
vec_tmp = 0.d0
do istate = 1, N_states
vec_tmp(:,istate) = reigvec_tc_bi_orth_tmp(:,istate)
enddo
do istate = N_states+1, n_states_diag
vec_tmp(istate,istate) = 1.d0
enddo
print*,'Diagonalizing the TC CISD '
call davidson_general_diag_dressed_ext_rout_nonsym_b1space(vec_tmp, H_jj, dressing_dets,eigval_tmp, N_det, n_states, n_states_diag, converged, htc_bi_ortho_calc_tdav_slow)
do i = 1, N_det
e_corr_dets(i) = reigvec_tc_bi_orth_tmp(i,1) * h0j(i)/reigvec_tc_bi_orth_tmp(1,1)
enddo
E_before = eigval_tmp(1)
print*,'Starting from ',E_before
e_current = 10.d0
thr = 1.d-5
it = 0
dressing_dets = 0.d0
double precision, allocatable :: H_jj(:),vec_tmp(:,:),eigval_tmp(:)
external htc_bi_ortho_calc_tdav_slow
external htcdag_bi_ortho_calc_tdav_slow
logical :: converged
do while (dabs(E_before-E_current).gt.thr)
it += 1
E_before = E_current
! h_sc2 = htilde_matrix_elmt_bi_ortho
call get_cisd_sc2_dressing(psi_det,e_corr_dets,N_det,dressing_dets)
do i = 1, N_det
! print*,'dressing_dets(i) = ',dressing_dets(i)
h_sc2(i,i) += dressing_dets(i)
enddo
print*,'********************'
print*,'iteration ',it
! call non_hrmt_real_diag(N_det,h_sc2,&
! leigvec_tc_bi_orth_tmp,reigvec_tc_bi_orth_tmp,&
! n_real,eigval_right_tmp)
! print*,'eigval_right_tmp(1)',eigval_right_tmp(1)
vec_tmp = 0.d0
do istate = 1, N_states
vec_tmp(:,istate) = reigvec_tc_bi_orth_tmp(:,istate)
enddo
do istate = N_states+1, n_states_diag
vec_tmp(istate,istate) = 1.d0
enddo
call davidson_general_diag_dressed_ext_rout_nonsym_b1space(vec_tmp, H_jj, dressing_dets,eigval_tmp, N_det, n_states, n_states_diag, converged, htc_bi_ortho_calc_tdav_slow)
print*,'outside Davidson'
print*,'eigval_tmp(1) = ',eigval_tmp(1)
do i = 1, N_det
reigvec_tc_bi_orth_tmp(i,1) = vec_tmp(i,1)
e_corr_dets(i) = reigvec_tc_bi_orth_tmp(i,1) * h0j(i)/reigvec_tc_bi_orth_tmp(1,1)
enddo
! E_current = eigval_right_tmp(1)
E_current = eigval_tmp(1)
print*,'it, E(SC)^2 = ',it,E_current
enddo
eigval_tc_cisd_sc2_bi_ortho(1:N_states) = eigval_right_tmp(1:N_states)
reigvec_tc_cisd_sc2_bi_ortho(1:N_det,1:N_states) = reigvec_tc_bi_orth_tmp(1:N_det,1:N_states)
leigvec_tc_cisd_sc2_bi_ortho(1:N_det,1:N_states) = leigvec_tc_bi_orth_tmp(1:N_det,1:N_states)
END_PROVIDER
subroutine get_cisd_sc2_dressing(dets,e_corr_dets,ndet,dressing_dets)
implicit none
use bitmasks
integer, intent(in) :: ndet
integer(bit_kind), intent(in) :: dets(N_int,2,ndet)
double precision, intent(in) :: e_corr_dets(ndet)
double precision, intent(out) :: dressing_dets(ndet)
integer, allocatable :: degree(:),hole(:,:),part(:,:),spin(:,:)
integer(bit_kind), allocatable :: hole_part(:,:,:)
integer :: i,j,k, exc(0:2,2,2),h1,p1,h2,p2,s1,s2
integer(bit_kind) :: xorvec(2,N_int)
double precision :: phase
dressing_dets = 0.d0
allocate(degree(ndet),hole(2,ndet),part(2,ndet), spin(2,ndet),hole_part(N_int,2,ndet))
do i = 2, ndet
call get_excitation_degree(HF_bitmask,dets(1,1,i),degree(i),N_int)
do j = 1, N_int
hole_part(j,1,i) = xor( HF_bitmask(j,1), dets(j,1,i))
hole_part(j,2,i) = xor( HF_bitmask(j,2), dets(j,2,i))
enddo
if(degree(i) == 1)then
call get_single_excitation(HF_bitmask,psi_det(1,1,i),exc,phase,N_int)
else if(degree(i) == 2)then
call get_double_excitation(HF_bitmask,psi_det(1,1,i),exc,phase,N_int)
endif
call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
hole(1,i) = h1
hole(2,i) = h2
part(1,i) = p1
part(2,i) = p2
spin(1,i) = s1
spin(2,i) = s2
enddo
integer :: same
if(elec_alpha_num+elec_beta_num<3)return
do i = 2, ndet
do j = i+1, ndet
same = 0
if(degree(i) == degree(j) .and. degree(i)==1)cycle
do k = 1, N_int
xorvec(k,1) = iand(hole_part(k,1,i),hole_part(k,1,j))
xorvec(k,2) = iand(hole_part(k,2,i),hole_part(k,2,j))
same += popcnt(xorvec(k,1)) + popcnt(xorvec(k,2))
enddo
! print*,'i,j',i,j
! call debug_det(dets(1,1,i),N_int)
! call debug_det(hole_part(1,1,i),N_int)
! call debug_det(dets(1,1,j),N_int)
! call debug_det(hole_part(1,1,j),N_int)
! print*,'same = ',same
if(same.eq.0)then
dressing_dets(i) += e_corr_dets(j)
dressing_dets(j) += e_corr_dets(i)
endif
enddo
enddo
end

8
plugins/local/tc_bi_ortho/tc_h_eigvectors.irp.f
View File

@ -326,7 +326,13 @@ end
enddo enddo
double precision, allocatable :: buffer(:,:) double precision, allocatable :: buffer(:,:)
allocate(buffer(N_det,N_states)) allocate(buffer(psi_det_size,N_states))
! print*,N_det,N_states
! print*,size(psi_l_coef_bi_ortho,1),size(psi_l_coef_bi_ortho,2)
! print*,size(leigvec_tc_bi_orth,1),size(leigvec_tc_bi_orth,2)
! print*,size(reigvec_tc_bi_orth,1),size(reigvec_tc_bi_orth,2)
! print*,size(psi_r_coef_bi_ortho,1),size(psi_r_coef_bi_ortho,2)
buffer = 0.d0
do k = 1, N_states do k = 1, N_states
do i = 1, N_det do i = 1, N_det
psi_l_coef_bi_ortho(i,k) = leigvec_tc_bi_orth(i,k) psi_l_coef_bi_ortho(i,k) = leigvec_tc_bi_orth(i,k)

10
plugins/local/tc_bi_ortho/tc_utils.irp.f
View File

@ -25,7 +25,7 @@ subroutine write_tc_energy()
E_2e_tmp(i) = 0.d0 E_2e_tmp(i) = 0.d0
E_3e_tmp(i) = 0.d0 E_3e_tmp(i) = 0.d0
do j = 1, N_det do j = 1, N_det
call htilde_mu_mat_bi_ortho_slow(psi_det(1,1,i), psi_det(1,1,j), N_int, hmono, htwoe, hthree, htot) call htilde_mu_mat_opt_bi_ortho(psi_det(1,1,i), psi_det(1,1,j), N_int, hmono, htwoe, hthree, htot)
E_TC_tmp(i) = E_TC_tmp(i) + psi_l_coef_bi_ortho(i,1) * psi_r_coef_bi_ortho(j,1) * htot E_TC_tmp(i) = E_TC_tmp(i) + psi_l_coef_bi_ortho(i,1) * psi_r_coef_bi_ortho(j,1) * htot
E_1e_tmp(i) = E_1e_tmp(i) + psi_l_coef_bi_ortho(i,1) * psi_r_coef_bi_ortho(j,1) * hmono E_1e_tmp(i) = E_1e_tmp(i) + psi_l_coef_bi_ortho(i,1) * psi_r_coef_bi_ortho(j,1) * hmono
E_2e_tmp(i) = E_2e_tmp(i) + psi_l_coef_bi_ortho(i,1) * psi_r_coef_bi_ortho(j,1) * htwoe E_2e_tmp(i) = E_2e_tmp(i) + psi_l_coef_bi_ortho(i,1) * psi_r_coef_bi_ortho(j,1) * htwoe
@ -70,7 +70,7 @@ subroutine write_tc_energy()
E_3e = 0.d0 E_3e = 0.d0
do i = 1, N_det do i = 1, N_det
do j = 1, N_det do j = 1, N_det
call htilde_mu_mat_bi_ortho_slow(psi_det(1,1,i), psi_det(1,1,j), N_int, hmono, htwoe, hthree, htot) call htilde_mu_mat_opt_bi_ortho(psi_det(1,1,i), psi_det(1,1,j), N_int, hmono, htwoe, hthree, htot)
E_TC = E_TC + psi_l_coef_bi_ortho(i,k) * psi_r_coef_bi_ortho(j,k) * htot E_TC = E_TC + psi_l_coef_bi_ortho(i,k) * psi_r_coef_bi_ortho(j,k) * htot
E_1e = E_1e + psi_l_coef_bi_ortho(i,k) * psi_r_coef_bi_ortho(j,k) * hmono E_1e = E_1e + psi_l_coef_bi_ortho(i,k) * psi_r_coef_bi_ortho(j,k) * hmono
E_2e = E_2e + psi_l_coef_bi_ortho(i,k) * psi_r_coef_bi_ortho(j,k) * htwoe E_2e = E_2e + psi_l_coef_bi_ortho(i,k) * psi_r_coef_bi_ortho(j,k) * htwoe
@ -109,8 +109,8 @@ subroutine write_tc_var()
SIGMA_TC = 0.d0 SIGMA_TC = 0.d0
do j = 2, N_det do j = 2, N_det
call htilde_mu_mat_bi_ortho_slow(psi_det(1,1,1), psi_det(1,1,j), N_int, hmono, htwoe, hthree, htot_1j) call htilde_mu_mat_opt_bi_ortho(psi_det(1,1,1), psi_det(1,1,j), N_int, hmono, htwoe, hthree, htot_1j)
call htilde_mu_mat_bi_ortho_slow(psi_det(1,1,j), psi_det(1,1,1), N_int, hmono, htwoe, hthree, htot_j1) call htilde_mu_mat_opt_bi_ortho(psi_det(1,1,j), psi_det(1,1,1), N_int, hmono, htwoe, hthree, htot_j1)
SIGMA_TC = SIGMA_TC + htot_1j * htot_j1 SIGMA_TC = SIGMA_TC + htot_1j * htot_j1
enddo enddo
@ -132,7 +132,7 @@ subroutine write_tc_gs_var_HF()
SIGMA_TC = 0.d0 SIGMA_TC = 0.d0
do j = 2, N_det do j = 2, N_det
call htilde_mu_mat_bi_ortho_slow(psi_det(1,1,j), psi_det(1,1,1), N_int, hmono, htwoe, hthree, htot) call htilde_mu_mat_opt_bi_ortho(psi_det(1,1,j), psi_det(1,1,1), N_int, hmono, htwoe, hthree, htot)
SIGMA_TC = SIGMA_TC + htot * htot SIGMA_TC = SIGMA_TC + htot * htot
enddo enddo

64
plugins/local/tc_bi_ortho/test_natorb.irp.f
View File

@ -1,64 +0,0 @@
! ---
program test_natorb
BEGIN_DOC
! TODO : Reads psi_det in the EZFIO folder and prints out the left- and right-eigenvectors together with the energy. Saves the left-right wave functions at the end.
END_DOC
implicit none
print *, 'Hello world'
my_grid_becke = .True.
PROVIDE tc_grid1_a tc_grid1_r
my_n_pt_r_grid = tc_grid1_r
my_n_pt_a_grid = tc_grid1_a
touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
read_wf = .True.
touch read_wf
call routine()
! call test()
end
! ---
subroutine routine()
implicit none
double precision, allocatable :: fock_diag(:),eigval(:),leigvec(:,:),reigvec(:,:),mat_ref(:,:)
allocate(eigval(mo_num),leigvec(mo_num,mo_num),reigvec(mo_num,mo_num),fock_diag(mo_num),mat_ref(mo_num, mo_num))
double precision, allocatable :: eigval_ref(:),leigvec_ref(:,:),reigvec_ref(:,:)
allocate(eigval_ref(mo_num),leigvec_ref(mo_num,mo_num),reigvec_ref(mo_num,mo_num))
double precision :: thr_deg
integer :: i,n_real,j
print*,'fock_matrix'
do i = 1, mo_num
fock_diag(i) = Fock_matrix_mo(i,i)
print*,i,fock_diag(i)
enddo
thr_deg = 1.d-6
mat_ref = -one_e_dm_mo
print*,'diagonalization by block'
call diag_mat_per_fock_degen(fock_diag,mat_ref,mo_num,thr_deg,leigvec,reigvec,eigval)
call non_hrmt_bieig( mo_num, mat_ref&
, leigvec_ref, reigvec_ref&
, n_real, eigval_ref)
print*,'TEST ***********************************'
double precision :: accu_l, accu_r
do i = 1, mo_num
accu_l = 0.d0
accu_r = 0.d0
do j = 1, mo_num
accu_r += reigvec_ref(j,i) * reigvec(j,i)
accu_l += leigvec_ref(j,i) * leigvec(j,i)
enddo
print*,i
write(*,'(I3,X,100(F16.10,X))')i,eigval(i),eigval_ref(i),accu_l,accu_r
enddo
end

173
plugins/local/tc_bi_ortho/test_normal_order.irp.f
View File

@ -1,173 +0,0 @@
! ---
program test_normal_order
BEGIN_DOC
! TODO : Put the documentation of the program here
END_DOC
implicit none
print *, 'Hello world'
my_grid_becke = .True.
PROVIDE tc_grid1_a tc_grid1_r
my_n_pt_r_grid = tc_grid1_r
my_n_pt_a_grid = tc_grid1_a
touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
read_wf = .True.
touch read_wf
call provide_all_three_ints_bi_ortho()
call test()
end
! ---
subroutine test
implicit none
use bitmasks ! you need to include the bitmasks_module.f90 features
integer :: h1,h2,p1,p2,s1,s2,i_ok,degree,Ne(2)
integer :: exc(0:2,2,2)
integer(bit_kind), allocatable :: det_i(:,:)
double precision :: hmono,htwoe,hthree,htilde_ij,accu,phase,normal,hthree_tmp
integer, allocatable :: occ(:,:)
allocate( occ(N_int*bit_kind_size,2) )
call bitstring_to_list_ab(ref_bitmask, occ, Ne, N_int)
allocate(det_i(N_int,2))
s1 = 1
s2 = 2
accu = 0.d0
do h1 = 1, elec_beta_num
do p1 = elec_alpha_num+1, mo_num
do h2 = 1, elec_beta_num
do p2 = elec_beta_num+1, mo_num
hthree = 0.d0
det_i = ref_bitmask
s1 = 1
s2 = 2
call do_single_excitation(det_i,h1,p1,s1,i_ok)
if(i_ok.ne.1)cycle
call do_single_excitation(det_i,h2,p2,s2,i_ok)
if(i_ok.ne.1)cycle
call htilde_mu_mat_bi_ortho_slow(det_i,HF_bitmask,N_int,hmono,htwoe,hthree_tmp,htilde_ij)
call get_excitation_degree(ref_bitmask,det_i,degree,N_int)
call get_excitation(ref_bitmask,det_i,exc,degree,phase,N_int)
hthree_tmp *= phase
hthree += 0.5d0 * hthree_tmp
det_i = ref_bitmask
s1 = 2
s2 = 1
call do_single_excitation(det_i,h1,p1,s1,i_ok)
if(i_ok.ne.1)cycle
call do_single_excitation(det_i,h2,p2,s2,i_ok)
if(i_ok.ne.1)cycle
call htilde_mu_mat_bi_ortho_slow(det_i,HF_bitmask,N_int,hmono,htwoe,hthree_tmp,htilde_ij)
call get_excitation_degree(ref_bitmask,det_i,degree,N_int)
call get_excitation(ref_bitmask,det_i,exc,degree,phase,N_int)
hthree_tmp *= phase
hthree += 0.5d0 * hthree_tmp
! normal = normal_two_body_bi_orth_ab(p2,h2,p1,h1)
call give_aba_contraction(N_int, h1, h2, p1, p2, Ne, occ, normal)
if(dabs(hthree).lt.1.d-10)cycle
if(dabs(hthree-normal).gt.1.d-10)then
! print*,pp2,pp1,hh2,hh1
print*,p2,p1,h2,h1
print*,hthree,normal,dabs(hthree-normal)
stop
endif
! call three_comp_two_e_elem(det_i,h1,h2,p1,p2,s1,s2,normal)
! normal = eff_2_e_from_3_e_ab(p2,p1,h2,h1)
accu += dabs(hthree-normal)
enddo
enddo
enddo
enddo
print*,'accu opposite spin = ',accu
stop
! p2=6
! p1=5
! h2=2
! h1=1
s1 = 1
s2 = 1
accu = 0.d0
do h1 = 1, elec_alpha_num
do p1 = elec_alpha_num+1, mo_num
do p2 = p1+1, mo_num
do h2 = h1+1, elec_alpha_num
det_i = ref_bitmask
call do_single_excitation(det_i,h1,p1,s1,i_ok)
if(i_ok.ne.1)cycle
call do_single_excitation(det_i,h2,p2,s2,i_ok)
if(i_ok.ne.1)cycle
call htilde_mu_mat_bi_ortho_slow(det_i,ref_bitmask,N_int,hmono,htwoe,hthree,htilde_ij)
call get_excitation_degree(ref_bitmask,det_i,degree,N_int)
call get_excitation(ref_bitmask,det_i,exc,degree,phase,N_int)
integer :: hh1, pp1, hh2, pp2, ss1, ss2
call decode_exc(exc, 2, hh1, pp1, hh2, pp2, ss1, ss2)
hthree *= phase
normal = normal_two_body_bi_orth_aa_bb(p2,h2,p1,h1)
! normal = eff_2_e_from_3_e_aa(p2,p1,h2,h1)
if(dabs(hthree).lt.1.d-10)cycle
if(dabs(hthree-normal).gt.1.d-10)then
print*,pp2,pp1,hh2,hh1
print*,p2,p1,h2,h1
print*,hthree,normal,dabs(hthree-normal)
stop
endif
! print*,hthree,normal,dabs(hthree-normal)
accu += dabs(hthree-normal)
enddo
enddo
enddo
enddo
print*,'accu same spin alpha = ',accu
s1 = 2
s2 = 2
accu = 0.d0
do h1 = 1, elec_beta_num
do p1 = elec_beta_num+1, mo_num
do p2 = p1+1, mo_num
do h2 = h1+1, elec_beta_num
det_i = ref_bitmask
call do_single_excitation(det_i,h1,p1,s1,i_ok)
if(i_ok.ne.1)cycle
call do_single_excitation(det_i,h2,p2,s2,i_ok)
if(i_ok.ne.1)cycle
call htilde_mu_mat_bi_ortho_slow(det_i,ref_bitmask,N_int,hmono,htwoe,hthree,htilde_ij)
call get_excitation_degree(ref_bitmask,det_i,degree,N_int)
call get_excitation(ref_bitmask,det_i,exc,degree,phase,N_int)
call decode_exc(exc, 2, hh1, pp1, hh2, pp2, ss1, ss2)
hthree *= phase
! normal = normal_two_body_bi_orth_aa_bb(p2,h2,p1,h1)
normal = eff_2_e_from_3_e_bb(p2,p1,h2,h1)
if(dabs(hthree).lt.1.d-10)cycle
if(dabs(hthree-normal).gt.1.d-10)then
print*,pp2,pp1,hh2,hh1
print*,p2,p1,h2,h1
print*,hthree,normal,dabs(hthree-normal)
stop
endif
! print*,hthree,normal,dabs(hthree-normal)
accu += dabs(hthree-normal)
enddo
enddo
enddo
enddo
print*,'accu same spin beta = ',accu
end

170
plugins/local/tc_bi_ortho/test_s2_tc.irp.f
View File

@ -1,170 +0,0 @@
! ---
program test_tc
implicit none
my_grid_becke = .True.
PROVIDE tc_grid1_a tc_grid1_r
my_n_pt_r_grid = tc_grid1_r
my_n_pt_a_grid = tc_grid1_a
touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
read_wf = .True.
touch read_wf
call provide_all_three_ints_bi_ortho()
call routine_h_triple_left
call routine_h_triple_right
! call routine_test_s2_davidson
end
subroutine routine_h_triple_right
implicit none
logical :: do_right
integer :: sze ,i, N_st, j
double precision :: sij, accu_e, accu_s, accu_e_0, accu_s_0
double precision, allocatable :: v_0_ref(:,:),u_0(:,:),s_0_ref(:,:)
double precision, allocatable :: v_0_new(:,:),s_0_new(:,:)
sze = N_det
N_st = 1
allocate(v_0_ref(N_det,1),u_0(N_det,1),s_0_ref(N_det,1),s_0_new(N_det,1),v_0_new(N_det,1))
print*,'Checking first the Right '
do i = 1, sze
u_0(i,1) = psi_r_coef_bi_ortho(i,1)
enddo
double precision :: wall0,wall1
call wall_time(wall0)
call H_tc_s2_u_0_with_pure_three_omp(v_0_ref,s_0_ref, u_0,N_st,sze)
call wall_time(wall1)
print*,'time for omp',wall1 - wall0
call wall_time(wall0)
call H_tc_s2_u_0_with_pure_three(v_0_new, s_0_new, u_0, N_st, sze)
call wall_time(wall1)
print*,'time serial ',wall1 - wall0
accu_e = 0.d0
accu_s = 0.d0
do i = 1, sze
accu_e += dabs(v_0_ref(i,1) - v_0_new(i,1))
accu_s += dabs(s_0_ref(i,1) - s_0_new(i,1))
enddo
print*,'accu_e = ',accu_e
print*,'accu_s = ',accu_s
end
subroutine routine_h_triple_left
implicit none
logical :: do_right
integer :: sze ,i, N_st, j
double precision :: sij, accu_e, accu_s, accu_e_0, accu_s_0
double precision, allocatable :: v_0_ref(:,:),u_0(:,:),s_0_ref(:,:)
double precision, allocatable :: v_0_new(:,:),s_0_new(:,:)
sze = N_det
N_st = 1
allocate(v_0_ref(N_det,1),u_0(N_det,1),s_0_ref(N_det,1),s_0_new(N_det,1),v_0_new(N_det,1))
print*,'Checking the Left '
do i = 1, sze
u_0(i,1) = psi_l_coef_bi_ortho(i,1)
enddo
double precision :: wall0,wall1
call wall_time(wall0)
call H_tc_s2_dagger_u_0_with_pure_three_omp(v_0_ref,s_0_ref, u_0,N_st,sze)
call wall_time(wall1)
print*,'time for omp',wall1 - wall0
call wall_time(wall0)
call H_tc_s2_dagger_u_0_with_pure_three(v_0_new, s_0_new, u_0, N_st, sze)
call wall_time(wall1)
print*,'time serial ',wall1 - wall0
accu_e = 0.d0
accu_s = 0.d0
do i = 1, sze
accu_e += dabs(v_0_ref(i,1) - v_0_new(i,1))
accu_s += dabs(s_0_ref(i,1) - s_0_new(i,1))
enddo
print*,'accu_e = ',accu_e
print*,'accu_s = ',accu_s
end
subroutine routine_test_s2_davidson
implicit none
double precision, allocatable :: H_jj(:),vec_tmp(:,:), energies(:) , s2(:)
integer :: i,istate
logical :: converged
external H_tc_s2_dagger_u_0_opt
external H_tc_s2_u_0_opt
allocate(H_jj(N_det),vec_tmp(N_det,n_states_diag),energies(n_states_diag), s2(n_states_diag))
do i = 1, N_det
call htilde_mu_mat_bi_ortho_tot_slow(psi_det(1,1,i), psi_det(1,1,i), N_int, H_jj(i))
enddo
! Preparing the left-eigenvector
print*,'Computing the left-eigenvector '
vec_tmp = 0.d0
do istate = 1, N_states
vec_tmp(1:N_det,istate) = psi_l_coef_bi_ortho(1:N_det,istate)
enddo
do istate = N_states+1, n_states_diag
vec_tmp(istate,istate) = 1.d0
enddo
do istate = 1, N_states
leigvec_tc_bi_orth(1:N_det,istate) = vec_tmp(1:N_det,istate)
enddo
integer :: n_it_max
n_it_max = 1
call davidson_hs2_nonsym_b1space(vec_tmp, H_jj, s2, energies, N_det, n_states, n_states_diag, n_it_max, converged, H_tc_s2_dagger_u_0_opt)
double precision, allocatable :: v_0_new(:,:),s_0_new(:,:)
integer :: sze,N_st
logical :: do_right
sze = N_det
N_st = 1
do_right = .False.
allocate(s_0_new(N_det,1),v_0_new(N_det,1))
call H_tc_s2_u_0_nstates_openmp(v_0_new,s_0_new,vec_tmp,N_st,sze, do_right)
double precision :: accu_e_0, accu_s_0
accu_e_0 = 0.d0
accu_s_0 = 0.d0
do i = 1, sze
accu_e_0 += v_0_new(i,1) * vec_tmp(i,1)
accu_s_0 += s_0_new(i,1) * vec_tmp(i,1)
enddo
print*,'energies = ',energies
print*,'s2 = ',s2
print*,'accu_e_0',accu_e_0
print*,'accu_s_0',accu_s_0
! Preparing the right-eigenvector
print*,'Computing the right-eigenvector '
vec_tmp = 0.d0
do istate = 1, N_states
vec_tmp(1:N_det,istate) = psi_r_coef_bi_ortho(1:N_det,istate)
enddo
do istate = N_states+1, n_states_diag
vec_tmp(istate,istate) = 1.d0
enddo
do istate = 1, N_states
leigvec_tc_bi_orth(1:N_det,istate) = vec_tmp(1:N_det,istate)
enddo
n_it_max = 1
call davidson_hs2_nonsym_b1space(vec_tmp, H_jj, s2, energies, N_det, n_states, n_states_diag, n_it_max, converged, H_tc_s2_u_0_opt)
sze = N_det
N_st = 1
do_right = .True.
v_0_new = 0.d0
s_0_new = 0.d0
call H_tc_s2_u_0_nstates_openmp(v_0_new,s_0_new,vec_tmp,N_st,sze, do_right)
accu_e_0 = 0.d0
accu_s_0 = 0.d0
do i = 1, sze
accu_e_0 += v_0_new(i,1) * vec_tmp(i,1)
accu_s_0 += s_0_new(i,1) * vec_tmp(i,1)
enddo
print*,'energies = ',energies
print*,'s2 = ',s2
print*,'accu_e_0',accu_e_0
print*,'accu_s_0',accu_s_0
end

171
plugins/local/tc_bi_ortho/test_tc_fock.irp.f
View File

@ -1,171 +0,0 @@
! ---
program test_tc_fock
BEGIN_DOC
! TODO : Put the documentation of the program here
END_DOC
implicit none
print *, 'Hello world'
my_grid_becke = .True.
PROVIDE tc_grid1_a tc_grid1_r
my_n_pt_r_grid = tc_grid1_r
my_n_pt_a_grid = tc_grid1_a
touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
read_wf = .True.
touch read_wf
!call routine_1
!call routine_2
! call routine_3()
call routine_tot
end
! ---
subroutine routine_3()
use bitmasks ! you need to include the bitmasks_module.f90 features
implicit none
integer :: i, a, i_ok, s1
double precision :: hmono, htwoe, hthree, htilde_ij
double precision :: err_ai, err_tot, ref, new
integer(bit_kind), allocatable :: det_i(:,:)
allocate(det_i(N_int,2))
err_tot = 0.d0
do s1 = 1, 2
det_i = ref_bitmask
call debug_det(det_i, N_int)
print*, ' HF det'
call debug_det(det_i, N_int)
do i = 1, elec_num_tab(s1)
do a = elec_num_tab(s1)+1, mo_num ! virtual
det_i = ref_bitmask
call do_single_excitation(det_i, i, a, s1, i_ok)
if(i_ok == -1) then
print*, 'PB !!'
print*, i, a
stop
endif
print*, ' excited det'
call debug_det(det_i, N_int)
call htilde_mu_mat_bi_ortho_slow(det_i, ref_bitmask, N_int, hmono, htwoe, hthree, htilde_ij)
if(dabs(hthree).lt.1.d-10)cycle
ref = hthree
if(s1 == 1)then
new = fock_a_tot_3e_bi_orth(a,i)
else if(s1 == 2)then
new = fock_b_tot_3e_bi_orth(a,i)
endif
err_ai = dabs(dabs(ref) - dabs(new))
if(err_ai .gt. 1d-7) then
print*,'s1 = ',s1
print*, ' warning on', i, a
print*, ref,new,err_ai
endif
print*, ref,new,err_ai
err_tot += err_ai
write(22, *) htilde_ij
enddo
enddo
enddo
print *, ' err_tot = ', err_tot
deallocate(det_i)
end subroutine routine_3
! ---
subroutine routine_tot()
use bitmasks ! you need to include the bitmasks_module.f90 features
implicit none
integer :: i, a, i_ok, s1,other_spin(2)
double precision :: hmono, htwoe, hthree, htilde_ij
double precision :: err_ai, err_tot, ref, new
integer(bit_kind), allocatable :: det_i(:,:)
allocate(det_i(N_int,2))
other_spin(1) = 2
other_spin(2) = 1
err_tot = 0.d0
! do s1 = 1, 2
s1 = 2
det_i = ref_bitmask
call debug_det(det_i, N_int)
print*, ' HF det'
call debug_det(det_i, N_int)
! do i = 1, elec_num_tab(s1)
! do a = elec_num_tab(s1)+1, mo_num ! virtual
do i = 1, elec_beta_num
do a = elec_beta_num+1, mo_num! virtual
print*,i,a
det_i = ref_bitmask
call do_single_excitation(det_i, i, a, s1, i_ok)
if(i_ok == -1) then
print*, 'PB !!'
print*, i, a
stop
endif
call htilde_mu_mat_bi_ortho_slow(det_i, ref_bitmask, N_int, hmono, htwoe, hthree, htilde_ij)
print*,htilde_ij
! if(dabs(htilde_ij).lt.1.d-10)cycle
print*, ' excited det'
call debug_det(det_i, N_int)
if(s1 == 1)then
new = Fock_matrix_tc_mo_alpha(a,i)
else
new = Fock_matrix_tc_mo_beta(a,i)
endif
ref = htilde_ij
! if(s1 == 1)then
! new = fock_a_tot_3e_bi_orth(a,i)
! else if(s1 == 2)then
! new = fock_b_tot_3e_bi_orth(a,i)
! endif
err_ai = dabs(dabs(ref) - dabs(new))
if(err_ai .gt. 1d-7) then
print*,'---------'
print*,'s1 = ',s1
print*, ' warning on', i, a
print*, ref,new,err_ai
print*,hmono, htwoe, hthree
print*,'---------'
endif
print*, ref,new,err_ai
err_tot += err_ai
write(22, *) htilde_ij
enddo
enddo
! enddo
print *, ' err_tot = ', err_tot
deallocate(det_i)
end subroutine routine_3

21
plugins/local/tc_keywords/EZFIO.cfg
View File

@ -14,7 +14,7 @@ default: False
type: logical type: logical
doc: If |true|, three-body terms are included doc: If |true|, three-body terms are included
interface: ezfio,provider,ocaml interface: ezfio,provider,ocaml
default: True default: False
[three_e_3_idx_term] [three_e_3_idx_term]
type: logical type: logical
@ -50,7 +50,7 @@ default: False
type: logical type: logical
doc: If |true|, standard normal-ordering for L (to be used with three_body_h_tc |false|) doc: If |true|, standard normal-ordering for L (to be used with three_body_h_tc |false|)
interface: ezfio,provider,ocaml interface: ezfio,provider,ocaml
default: False default: True
[core_tc_op] [core_tc_op]
type: logical type: logical
@ -184,12 +184,6 @@ doc: Read/Write normal_two_body_bi_orth from/to disk [ Write | Read | None ]
interface: ezfio,provider,ocaml interface: ezfio,provider,ocaml
default: None default: None
[debug_tc_pt2]
type: integer
doc: If :: 1 then you compute the TC-PT2 the old way, :: 2 then you check with the new version but without three-body
interface: ezfio,provider,ocaml
default: -1
[only_spin_tc_right] [only_spin_tc_right]
type: logical type: logical
doc: If |true|, only the right part of WF is used to compute spin dens doc: If |true|, only the right part of WF is used to compute spin dens
@ -244,6 +238,12 @@ doc: If |true|, you minimize the angle between the left and right vectors associ
interface: ezfio,provider,ocaml interface: ezfio,provider,ocaml
default: False default: False
[thresh_de_tc_angles]
type: Threshold
doc: Thresholds on delta E for changing angles between orbitals
interface: ezfio,provider,ocaml
default: 1.e-6
[ao_to_mo_tc_n3] [ao_to_mo_tc_n3]
type: logical type: logical
doc: If |true|, memory scale of TC ao -> mo: O(N3) doc: If |true|, memory scale of TC ao -> mo: O(N3)
@ -268,3 +268,8 @@ doc: Thresholds on the Imag part of TC energy
interface: ezfio,provider,ocaml interface: ezfio,provider,ocaml
default: 1.e-7 default: 1.e-7
[transpose_two_e_int]
type: logical
doc: If |true|, you duplicate the two-electron TC integrals with the transpose matrix. Acceleates the PT2.
interface: ezfio,provider,ocaml
default: False

1
plugins/local/tc_progs/NEED Normal file
View File

@ -0,0 +1 @@
tc_bi_ortho

0
plugins/local/tc_bi_ortho/print_he_tc_energy.irp.f → plugins/local/tc_progs/print_he_tc_energy.irp.f
View File

0
plugins/local/tc_bi_ortho/print_tc_dump.irp.f → plugins/local/tc_progs/print_tc_dump.irp.f
View File

0
plugins/local/tc_bi_ortho/print_tc_energy.irp.f → plugins/local/tc_progs/print_tc_energy.irp.f
View File

0
plugins/local/tc_bi_ortho/print_tc_spin_dens.irp.f → plugins/local/tc_progs/print_tc_spin_dens.irp.f
View File

0
plugins/local/tc_bi_ortho/print_tc_var.irp.f → plugins/local/tc_progs/print_tc_var.irp.f
View File

8
plugins/local/tc_bi_ortho/print_tc_wf.irp.f → plugins/local/tc_progs/print_tc_wf.irp.f
View File

@ -37,7 +37,7 @@ subroutine write_l_r_wf
integer :: i integer :: i
print*,'Writing the left-right wf' print*,'Writing the left-right wf'
do i = 1, N_det do i = 1, N_det
write(i_unit_output,*)i, psi_coef_sorted_tc(i,1)/psi_coef_sorted_tc(i,1) & write(i_unit_output,'(I8,X,10(F16.10,X))')i, psi_coef_sorted_tc(i,1),psi_coef_sorted_tc(i,1)/psi_coef_sorted_tc(1,1)&
, psi_l_coef_sorted_bi_ortho_left(i)/psi_l_coef_sorted_bi_ortho_left(1) & , psi_l_coef_sorted_bi_ortho_left(i)/psi_l_coef_sorted_bi_ortho_left(1) &
, psi_r_coef_sorted_bi_ortho_right(i)/psi_r_coef_sorted_bi_ortho_right(1) , psi_r_coef_sorted_bi_ortho_right(i)/psi_r_coef_sorted_bi_ortho_right(1)
enddo enddo
@ -61,12 +61,12 @@ subroutine routine
do i = 1, N_det do i = 1, N_det
call get_excitation_degree(HF_bitmask,psi_det(1,1,i),degree,N_int) call get_excitation_degree(HF_bitmask,psi_det(1,1,i),degree,N_int)
if(degree == 1 .or. degree == 2)then if(degree == 1 .or. degree == 2)then
call htilde_mu_mat_bi_ortho_slow(psi_det(1,1,i),HF_bitmask,N_int,hmono,htwoe,hthree,htilde_ij) call htilde_mu_mat_opt_bi_ortho(psi_det(1,1,i),HF_bitmask,N_int,hmono,htwoe,hthree,htilde_ij)
call htilde_mu_mat_bi_ortho_slow(psi_det(1,1,i),psi_det(1,1,i),N_int,hmono,htwoe,hthree,e_i0) call htilde_mu_mat_opt_bi_ortho(psi_det(1,1,i),psi_det(1,1,i),N_int,hmono,htwoe,hthree,e_i0)
delta_e = e_tilde_00 - e_i0 delta_e = e_tilde_00 - e_i0
coef_pt1 = htilde_ij / delta_e coef_pt1 = htilde_ij / delta_e
call htilde_mu_mat_bi_ortho_slow(HF_bitmask,psi_det(1,1,i),N_int,hmono,htwoe,hthree,htilde_ij) call htilde_mu_mat_opt_bi_ortho(HF_bitmask,psi_det(1,1,i),N_int,hmono,htwoe,hthree,htilde_ij)
contrib_pt = coef_pt1 * htilde_ij contrib_pt = coef_pt1 * htilde_ij
e_pt2 += contrib_pt e_pt2 += contrib_pt

0
plugins/local/tc_bi_ortho/save_bitcpsileft_for_qmcchem.irp.f → plugins/local/tc_progs/save_bitcpsileft_for_qmcchem.irp.f
View File

0
plugins/local/tc_bi_ortho/save_tc_bi_ortho_nat.irp.f → plugins/local/tc_progs/save_tc_bi_ortho_nat.irp.f
View File

0
plugins/local/tc_bi_ortho/select_dets_bi_ortho.irp.f → plugins/local/tc_progs/select_dets_bi_ortho.irp.f
View File

0
plugins/local/tc_bi_ortho/tc_bi_ortho_prop.irp.f → plugins/local/tc_progs/tc_bi_ortho_prop.irp.f
View File

4
plugins/local/tc_bi_ortho/tc_som.irp.f → plugins/local/tc_progs/tc_som.irp.f
View File

@ -49,8 +49,8 @@ subroutine main()
U_SOM = 0.d0 U_SOM = 0.d0
do i = 1, N_det do i = 1, N_det
if(i == i_HF) cycle if(i == i_HF) cycle
call htilde_mu_mat_bi_ortho_slow(psi_det(1,1,i_HF), psi_det(1,1,i), N_int, hmono_1, htwoe_1, hthree_1, htot_1) call htilde_mu_mat_opt_bi_ortho(psi_det(1,1,i_HF), psi_det(1,1,i), N_int, hmono_1, htwoe_1, hthree_1, htot_1)
call htilde_mu_mat_bi_ortho_slow(psi_det(1,1,i), psi_det(1,1,i_HF), N_int, hmono_2, htwoe_2, hthree_2, htot_2) call htilde_mu_mat_opt_bi_ortho(psi_det(1,1,i), psi_det(1,1,i_HF), N_int, hmono_2, htwoe_2, hthree_2, htot_2)
U_SOM += htot_1 * htot_2 U_SOM += htot_1 * htot_2
enddo enddo
U_SOM = 0.5d0 * U_SOM U_SOM = 0.5d0 * U_SOM

0
plugins/local/tc_bi_ortho/test_tc_two_rdm.irp.f → plugins/local/tc_progs/test_tc_two_rdm.irp.f
View File

1
plugins/local/tc_scf/routines_rotates.irp.f
View File

@ -304,6 +304,7 @@ subroutine routine_save_rotated_mos(thr_deg, good_angles)
! check if TC energy has changed ! check if TC energy has changed
E_new = TC_HF_energy E_new = TC_HF_energy
E_thr = thresh_de_tc_angles
if(dabs(E_new - E_old) .gt. E_thr) then if(dabs(E_new - E_old) .gt. E_thr) then
mo_r_coef = mo_r_coef_old mo_r_coef = mo_r_coef_old
mo_l_coef = mo_l_coef_old mo_l_coef = mo_l_coef_old

53
plugins/local/tuto_plugins/tuto_I/test_cholesky.irp.f Normal file
View File

@ -0,0 +1,53 @@
program my_program_to_print_stuffs
implicit none
BEGIN_DOC
! TODO : Put the documentation of the program here
END_DOC
integer :: i,j,k,l,m
double precision :: integral, accu, accu_tot, integral_cholesky
double precision :: get_ao_two_e_integral, get_two_e_integral ! declaration of the functions
print*,'AO integrals, physicist notations : <i j|k l>'
accu_tot = 0.D0
do i = 1, ao_num
do j = 1, ao_num
do k = 1, ao_num
do l = 1, ao_num
integral = get_ao_two_e_integral(i, j, k, l, ao_integrals_map)
integral_cholesky = 0.D0
do m = 1, cholesky_ao_num
integral_cholesky += cholesky_ao_transp(m,i,k) * cholesky_ao_transp(m,j,l)
enddo
accu = dabs(integral_cholesky-integral)
accu_tot += accu
if(accu.gt.1.d-10)then
print*,i,j,k,l
print*,accu, integral, integral_cholesky
endif
enddo
enddo
enddo
enddo
print*,'accu_tot',accu_tot
print*,'MO integrals, physicist notations : <i j|k l>'
do i = 1, mo_num
do j = 1, mo_num
do k = 1, mo_num
do l = 1, mo_num
integral = get_two_e_integral(i, j, k, l, mo_integrals_map)
accu = 0.D0
integral_cholesky = 0.D0
do m = 1, cholesky_mo_num
integral_cholesky += cholesky_mo_transp(m,i,k) * cholesky_mo_transp(m,j,l)
enddo
accu = dabs(integral_cholesky-integral)
accu_tot += accu
if(accu.gt.1.d-10)then
print*,i,j,k,l
print*,accu, integral, integral_cholesky
endif
enddo
enddo
enddo
enddo
end

8
src/ao_two_e_ints/EZFIO.cfg
View File

@ -6,7 +6,7 @@ default: None
[io_ao_cholesky] [io_ao_cholesky]
type: Disk_access type: Disk_access
doc: Read/Write |AO| integrals from/to disk [ Write | Read | None ] doc: Read/Write |AO| Cholesky integrals from/to disk [ Write | Read | None ]
interface: ezfio,provider,ocaml interface: ezfio,provider,ocaml
default: None default: None
@ -25,16 +25,16 @@ default: 1.e-12
[do_direct_integrals] [do_direct_integrals]
type: logical type: logical
doc: Compute integrals on the fly (very slow, only for debugging) doc: Compute integrals on the fly (Useful only for Cholesky decomposition)
interface: ezfio,provider,ocaml interface: ezfio,provider,ocaml
default: False default: True
ezfio_name: direct ezfio_name: direct
[do_ao_cholesky] [do_ao_cholesky]
type: logical type: logical
doc: Perform Cholesky decomposition of AO integrals doc: Perform Cholesky decomposition of AO integrals
interface: ezfio,provider,ocaml interface: ezfio,provider,ocaml
default: False default: True
[io_ao_two_e_integrals_erf] [io_ao_two_e_integrals_erf]
type: Disk_access type: Disk_access

384
src/ao_two_e_ints/cholesky.irp.f
View File

@ -6,7 +6,7 @@ BEGIN_PROVIDER [ double precision, cholesky_ao_transp, (cholesky_ao_num, ao_num,
integer :: i,j,k integer :: i,j,k
do j=1,ao_num do j=1,ao_num
do i=1,ao_num do i=1,ao_num
do k=1,ao_num do k=1,cholesky_ao_num
cholesky_ao_transp(k,i,j) = cholesky_ao(i,j,k) cholesky_ao_transp(k,i,j) = cholesky_ao(i,j,k)
enddo enddo
enddo enddo
@ -16,27 +16,32 @@ END_PROVIDER
BEGIN_PROVIDER [ integer, cholesky_ao_num ] BEGIN_PROVIDER [ integer, cholesky_ao_num ]
&BEGIN_PROVIDER [ double precision, cholesky_ao, (ao_num, ao_num, 1) ] &BEGIN_PROVIDER [ double precision, cholesky_ao, (ao_num, ao_num, 1) ]
use mmap_module
implicit none implicit none
BEGIN_DOC BEGIN_DOC
! Cholesky vectors in AO basis: (ik|a): ! Cholesky vectors in AO basis: (ik|a):
! <ij|kl> = (ik|jl) = sum_a (ik|a).(a|jl) ! <ij|kl> = (ik|jl) = sum_a (ik|a).(a|jl)
! !
! Last dimension of cholesky_ao is cholesky_ao_num ! Last dimension of cholesky_ao is cholesky_ao_num
!
! https://mogp-emulator.readthedocs.io/en/latest/methods/proc/ProcPivotedCholesky.html
! https://doi.org/10.1016/j.apnum.2011.10.001 : Page 4, Algorithm 1
END_DOC END_DOC
integer :: rank, ndim integer*8 :: ndim8
double precision :: tau integer :: rank
double precision, pointer :: L(:,:), L_old(:,:) double precision :: tau, tau2
double precision, pointer :: L(:,:)
double precision :: s double precision :: s
double precision, parameter :: dscale = 1.d0
double precision, allocatable :: D(:), Delta(:,:), Ltmp_p(:,:), Ltmp_q(:,:) double precision, allocatable :: D(:), Ltmp_p(:,:), Ltmp_q(:,:), D_sorted(:), Delta_col(:), Delta(:,:)
integer, allocatable :: Lset(:), Dset(:), addr(:,:) integer, allocatable :: addr1(:), addr2(:)
integer*8, allocatable :: Lset(:), Dset(:)
logical, allocatable :: computed(:) logical, allocatable :: computed(:)
integer :: i,j,k,m,p,q, qj, dj, p2, q2 integer :: i,j,k,m,p,q, dj, p2, q2, ii, jj
integer*8 :: i8, j8, p8, qj8, rank_max, np8
integer :: N, np, nq integer :: N, np, nq
double precision :: Dmax, Dmin, Qmax, f double precision :: Dmax, Dmin, Qmax, f
@ -44,19 +49,32 @@ END_PROVIDER
logical, external :: ao_two_e_integral_zero logical, external :: ao_two_e_integral_zero
double precision, external :: ao_two_e_integral double precision, external :: ao_two_e_integral
integer :: block_size, iblock, ierr integer :: block_size, iblock
double precision :: mem double precision :: mem, mem0
double precision, external :: memory_of_double, memory_of_int double precision, external :: memory_of_double, memory_of_int
double precision, external :: memory_of_double8, memory_of_int8
integer, external :: getUnitAndOpen integer, external :: getUnitAndOpen
integer :: iunit integer :: iunit, ierr
ndim = ao_num*ao_num ndim8 = ao_num*ao_num*1_8
double precision :: wall0,wall1
type(c_ptr) :: c_pointer(2)
integer :: fd(2)
PROVIDE nproc ao_cholesky_threshold do_direct_integrals qp_max_mem
PROVIDE nucl_coord ao_two_e_integral_schwartz
call set_multiple_levels_omp(.False.)
call wall_time(wall0)
! Will be reallocated at the end
deallocate(cholesky_ao) deallocate(cholesky_ao)
if (read_ao_cholesky) then if (read_ao_cholesky) then
print *, 'Reading Cholesky vectors from disk...' print *, 'Reading Cholesky AO vectors from disk...'
iunit = getUnitAndOpen(trim(ezfio_work_dir)//'cholesky_ao', 'R') iunit = getUnitAndOpen(trim(ezfio_work_dir)//'cholesky_ao', 'R')
read(iunit) rank read(iunit) rank
allocate(cholesky_ao(ao_num,ao_num,rank), stat=ierr) allocate(cholesky_ao(ao_num,ao_num,rank), stat=ierr)
@ -66,7 +84,6 @@ END_PROVIDER
else else
PROVIDE nucl_coord ao_two_e_integral_schwartz
call set_multiple_levels_omp(.False.) call set_multiple_levels_omp(.False.)
if (do_direct_integrals) then if (do_direct_integrals) then
@ -79,14 +96,17 @@ END_PROVIDER
endif endif
tau = ao_cholesky_threshold tau = ao_cholesky_threshold
tau2 = tau*tau
mem = 6.d0 * memory_of_double(ndim) + 6.d0 * memory_of_int(ndim) rank = 0
call check_mem(mem, irp_here)
allocate( D(ndim8), Lset(ndim8), Dset(ndim8), D_sorted(ndim8))
allocate( addr1(ndim8), addr2(ndim8), Delta_col(ndim8), computed(ndim8) )
call resident_memory(mem0)
call print_memory_usage() call print_memory_usage()
allocate(L(ndim,1))
print *, '' print *, ''
print *, 'Cholesky decomposition of AO integrals' print *, 'Cholesky decomposition of AO integrals'
print *, '======================================' print *, '======================================'
@ -96,49 +116,58 @@ END_PROVIDER
print *, '============ =============' print *, '============ ============='
rank = 0
allocate( D(ndim), Lset(ndim), Dset(ndim) )
allocate( addr(3,ndim) )
! 1. ! 1.
k=0 i8=0
do j=1,ao_num do j=1,ao_num
do i=1,ao_num do i=1,ao_num
k = k+1 i8 = i8+1
addr(1,k) = i addr1(i8) = i
addr(2,k) = j addr2(i8) = j
addr(3,k) = (i-1)*ao_num + j
enddo enddo
enddo enddo
if (do_direct_integrals) then if (do_direct_integrals) then
!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i) SCHEDULE(guided) !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i8) SCHEDULE(dynamic,21)
do i=1,ndim do i8=ndim8,1,-1
D(i) = ao_two_e_integral(addr(1,i), addr(2,i), & D(i8) = ao_two_e_integral(addr1(i8), addr2(i8), &
addr(1,i), addr(2,i)) addr1(i8), addr2(i8))
enddo enddo
!$OMP END PARALLEL DO !$OMP END PARALLEL DO
else else
!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i) SCHEDULE(guided) !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i8) SCHEDULE(dynamic,21)
do i=1,ndim do i8=ndim8,1,-1
D(i) = get_ao_two_e_integral(addr(1,i), addr(1,i), & D(i8) = get_ao_two_e_integral(addr1(i8), addr1(i8), &
addr(2,i), addr(2,i), & addr2(i8), addr2(i8), ao_integrals_map)
ao_integrals_map)
enddo enddo
!$OMP END PARALLEL DO !$OMP END PARALLEL DO
endif endif
Dmax = maxval(D) D_sorted(:) = -D(:)
call dsort_noidx_big(D_sorted,ndim8)
D_sorted(:) = -D_sorted(:)
Dmax = D_sorted(1)
! 2. ! 2.
np=0 np8=0_8
do p=1,ndim do p8=1,ndim8
if ( dscale*dscale*Dmax*D(p) > tau*tau ) then if ( Dmax*D(p8) >= tau2 ) then
np = np+1 np8 = np8+1_8
Lset(np) = p Lset(np8) = p8
endif endif
enddo enddo
np = np8
if (np <= 0) stop 'np<=0'
if (np > ndim8) stop 'np>ndim8'
rank_max = min(np,20*elec_num*elec_num)
call mmap(trim(ezfio_work_dir)//'cholesky_ao_tmp', (/ ndim8, rank_max /), 8, fd(1), .False., .True., c_pointer(1))
call c_f_pointer(c_pointer(1), L, (/ ndim8, rank_max /))
! Deleting the file while it is open makes the file invisible on the filesystem,
! and automatically deleted, even if the program crashes
iunit = getUnitAndOpen(trim(ezfio_work_dir)//'cholesky_ao_tmp', 'R')
close(iunit,status='delete')
! 3. ! 3.
N = 0 N = 0
@ -146,77 +175,67 @@ END_PROVIDER
! 4. ! 4.
i = 0 i = 0
mem = memory_of_double(np) & ! Delta(np,nq)
+ (np+1)*memory_of_double(block_size) ! Ltmp_p(np,block_size) + Ltmp_q(nq,block_size)
! call check_mem(mem)
! 5. ! 5.
do while ( (Dmax > tau).and.(rank < ndim) ) do while ( (Dmax > tau).and.(np > 0) )
! a. ! a.
i = i+1 i = i+1
s = 0.01d0
! Inrease s until the arrays fit in memory block_size = max(N,24)
! Determine nq so that Delta fits in memory
s = 0.1d0
Dmin = max(s*Dmax,tau)
do nq=2,np-1
if (D_sorted(nq) < Dmin) exit
enddo
do while (.True.) do while (.True.)
! b. mem = mem0 &
Dmin = max(s*Dmax,tau) + np*memory_of_double(nq) & ! Delta(np,nq)
+ (np+nq)*memory_of_double(block_size) ! Ltmp_p(np,block_size) + Ltmp_q(nq,block_size)
! c. if (mem > qp_max_mem*0.5d0) then
Dmin = D_sorted(nq/2)
do ii=nq/2,np-1
if (D_sorted(ii) < Dmin) then
nq = ii
exit
endif
enddo
else
exit
endif
enddo
!call print_memory_usage
!print *, 'np, nq, Predicted memory: ', np, nq, mem
if (nq <= 0) then
print *, nq
stop 'bug in cholesky: nq <= 0'
endif
Dmin = D_sorted(nq)
nq=0 nq=0
do p=1,np do p=1,np
if ( D(Lset(p)) > Dmin ) then if ( D(Lset(p)) >= Dmin ) then
nq = nq+1 nq = nq+1
Dset(nq) = Lset(p) Dset(nq) = Lset(p)
endif endif
enddo enddo
call total_memory(mem)
mem = mem &
+ np*memory_of_double(nq) &! Delta(np,nq)
+ (rank+nq)* memory_of_double(ndim) &! L(ndim,rank+nq)
+ (np+nq)*memory_of_double(block_size) ! Ltmp_p(np,block_size) + Ltmp_q(nq,block_size)
if (mem > qp_max_mem) then
s = s*2.d0
else
exit
endif
if ((s > 1.d0).or.(nq == 0)) then
call print_memory_usage()
print *, 'Not enough memory. Reduce cholesky threshold'
stop -1
endif
enddo
! d., e.
block_size = max(N,24)
L_old => L
allocate(L(ndim,rank+nq), stat=ierr)
if (ierr /= 0) then
call print_memory_usage()
print *, irp_here, ': allocation failed : (L(ndim,rank+nq))'
stop -1
endif
!$OMP PARALLEL DO PRIVATE(k,j)
do k=1,rank
do j=1,ndim
L(j,k) = L_old(j,k)
enddo
enddo
!$OMP END PARALLEL DO
deallocate(L_old)
allocate(Delta(np,nq), stat=ierr)
if (ierr /= 0) then
call print_memory_usage()
print *, irp_here, ': allocation failed : (Delta(np,nq))'
stop -1
endif
allocate(Delta(np,nq))
allocate(Ltmp_p(np,block_size), stat=ierr) allocate(Ltmp_p(np,block_size), stat=ierr)
if (ierr /= 0) then if (ierr /= 0) then
call print_memory_usage() call print_memory_usage()
print *, irp_here, ': allocation failed : (Ltmp_p(np,block_size))' print *, irp_here, ': allocation failed : (Ltmp_p(np,block_size))'
@ -224,6 +243,7 @@ END_PROVIDER
endif endif
allocate(Ltmp_q(nq,block_size), stat=ierr) allocate(Ltmp_q(nq,block_size), stat=ierr)
if (ierr /= 0) then if (ierr /= 0) then
call print_memory_usage() call print_memory_usage()
print *, irp_here, ': allocation failed : (Ltmp_q(nq,block_size))' print *, irp_here, ': allocation failed : (Ltmp_q(nq,block_size))'
@ -231,36 +251,39 @@ END_PROVIDER
endif endif
allocate(computed(nq)) computed(1:nq) = .False.
!$OMP PARALLEL DEFAULT(SHARED) PRIVATE(m,k,p,q,j)
!$OMP DO !$OMP PARALLEL DEFAULT(SHARED) PRIVATE(k,p,q)
do q=1,nq
do j=1,np
Delta(j,q) = 0.d0
enddo
computed(q) = .False.
enddo
!$OMP ENDDO NOWAIT
!$OMP DO
do k=1,N do k=1,N
!$OMP DO
do p=1,np do p=1,np
Ltmp_p(p,k) = L(Lset(p),k) Ltmp_p(p,k) = L(Lset(p),k)
enddo enddo
!$OMP END DO NOWAIT
!$OMP DO
do q=1,nq do q=1,nq
Ltmp_q(q,k) = L(Dset(q),k) Ltmp_q(q,k) = L(Dset(q),k)
enddo enddo
enddo
!$OMP END DO NOWAIT !$OMP END DO NOWAIT
enddo
!$OMP BARRIER !$OMP BARRIER
!$OMP END PARALLEL !$OMP END PARALLEL
if (N>0) then if (N>0) then
call dgemm('N','T', np, nq, N, -1.d0, &
Ltmp_p, np, Ltmp_q, nq, 1.d0, Delta, np) call dgemm('N', 'T', np, nq, N, -1.d0, &
Ltmp_p(1,1), np, Ltmp_q(1,1), nq, 0.d0, Delta, np)
else
!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(q,j)
do q=1,nq
Delta(:,q) = 0.d0
enddo
!$OMP END PARALLEL DO
endif endif
! f. ! f.
@ -274,51 +297,81 @@ END_PROVIDER
iblock = 0 iblock = 0
do j=1,nq do j=1,nq
if ( (Qmax <= Dmin).or.(N+j > ndim) ) exit if ( (Qmax <= Dmin).or.(N+j*1_8 > ndim8) ) exit
! i. ! i.
rank = N+j rank = N+j
if (rank == rank_max) then
print *, 'cholesky: rank_max reached'
exit
endif
if (iblock == block_size) then if (iblock == block_size) then
call dgemm('N','T',np,nq,block_size,-1.d0, & call dgemm('N','T',np,nq,block_size,-1.d0, &
Ltmp_p, np, Ltmp_q, nq, 1.d0, Delta, np) Ltmp_p, np, Ltmp_q, nq, 1.d0, Delta, np)
iblock = 0 iblock = 0
endif endif
! ii. ! ii.
do dj=1,nq do dj=1,nq
qj = Dset(dj) qj8 = Dset(dj)
if (D(qj) == Qmax) then if (D(qj8) == Qmax) then
exit exit
endif endif
enddo enddo
L(1:ndim, rank) = 0.d0 do i8=1,ndim8
L(i8, rank) = 0.d0
if (.not.computed(dj)) then
m = dj
!$OMP PARALLEL DO PRIVATE(k) SCHEDULE(guided)
do k=np,1,-1
if (.not.ao_two_e_integral_zero( addr(1,Lset(k)), addr(1,Dset(m)),&
addr(2,Lset(k)), addr(2,Dset(m)) ) ) then
if (do_direct_integrals) then
Delta(k,m) = Delta(k,m) + &
ao_two_e_integral(addr(1,Lset(k)), addr(2,Lset(k)),&
addr(1,Dset(m)), addr(2,Dset(m)))
else
Delta(k,m) = Delta(k,m) + &
get_ao_two_e_integral( addr(1,Lset(k)), addr(1,Dset(m)),&
addr(2,Lset(k)), addr(2,Dset(m)), ao_integrals_map)
endif
endif
enddo enddo
!$OMP END PARALLEL DO
computed(dj) = .True.
endif
iblock = iblock+1 iblock = iblock+1
!$OMP PARALLEL DO PRIVATE(p)
do p=1,np do p=1,np
Ltmp_p(p,iblock) = Delta(p,dj) Ltmp_p(p,iblock) = Delta(p,dj)
enddo enddo
!$OMP END PARALLEL DO
if (.not.computed(dj)) then
m = dj
if (do_direct_integrals) then
!$OMP PARALLEL DO PRIVATE(k) SCHEDULE(dynamic,21)
do k=1,np
Delta_col(k) = 0.d0
if (.not.ao_two_e_integral_zero( addr1(Lset(k)), addr1(Dset(m)),&
addr2(Lset(k)), addr2(Dset(m)) ) ) then
Delta_col(k) = &
ao_two_e_integral(addr1(Lset(k)), addr2(Lset(k)),&
addr1(Dset(m)), addr2(Dset(m)))
endif
enddo
!$OMP END PARALLEL DO
else
PROVIDE ao_integrals_map
!$OMP PARALLEL DO PRIVATE(k) SCHEDULE(dynamic,21)
do k=1,np
Delta_col(k) = 0.d0
if (.not.ao_two_e_integral_zero( addr1(Lset(k)), addr1(Dset(m)),&
addr2(Lset(k)), addr2(Dset(m)) ) ) then
Delta_col(k) = &
get_ao_two_e_integral( addr1(Lset(k)), addr1(Dset(m)),&
addr2(Lset(k)), addr2(Dset(m)), ao_integrals_map)
endif
enddo
!$OMP END PARALLEL DO
endif
!$OMP PARALLEL DO PRIVATE(p)
do p=1,np
Ltmp_p(p,iblock) = Ltmp_p(p,iblock) + Delta_col(p)
Delta(p,dj) = Ltmp_p(p,iblock)
enddo
!$OMP END PARALLEL DO
computed(dj) = .True.
endif
! iv. ! iv.
if (iblock > 1) then if (iblock > 1) then
@ -329,7 +382,7 @@ END_PROVIDER
! iii. ! iii.
f = 1.d0/dsqrt(Qmax) f = 1.d0/dsqrt(Qmax)
!$OMP PARALLEL PRIVATE(m,p,q,k) DEFAULT(shared) !$OMP PARALLEL PRIVATE(p,q) DEFAULT(shared)
!$OMP DO !$OMP DO
do p=1,np do p=1,np
Ltmp_p(p,iblock) = Ltmp_p(p,iblock) * f Ltmp_p(p,iblock) = Ltmp_p(p,iblock) * f
@ -343,7 +396,6 @@ END_PROVIDER
Ltmp_q(q,iblock) = L(Dset(q), rank) Ltmp_q(q,iblock) = L(Dset(q), rank)
enddo enddo
!$OMP END DO !$OMP END DO
!$OMP END PARALLEL !$OMP END PARALLEL
Qmax = D(Dset(1)) Qmax = D(Dset(1))
@ -355,49 +407,62 @@ END_PROVIDER
print '(I10, 4X, ES12.3)', rank, Qmax print '(I10, 4X, ES12.3)', rank, Qmax
deallocate(computed)
deallocate(Delta)
deallocate(Ltmp_p) deallocate(Ltmp_p)
deallocate(Ltmp_q) deallocate(Ltmp_q)
deallocate(Delta)
! i. ! i.
N = rank N = rank
! j. ! j.
Dmax = D(Lset(1)) D_sorted(:) = -D(:)
do p=1,np call dsort_noidx_big(D_sorted,ndim8)
Dmax = max(Dmax, D(Lset(p))) D_sorted(:) = -D_sorted(:)
enddo
np=0 Dmax = D_sorted(1)
do p=1,ndim
if ( dscale*dscale*Dmax*D(p) > tau*tau ) then np8=0_8
np = np+1 do p8=1,ndim8
Lset(np) = p if ( Dmax*D(p8) >= tau2 ) then
np8 = np8+1_8
Lset(np8) = p8
endif endif
enddo enddo
np = np8
enddo enddo
print *, '============ ============='
print *, ''
deallocate( D, Lset, Dset, D_sorted )
deallocate( addr1, addr2, Delta_col, computed )
allocate(cholesky_ao(ao_num,ao_num,rank), stat=ierr) allocate(cholesky_ao(ao_num,ao_num,rank), stat=ierr)
if (ierr /= 0) then if (ierr /= 0) then
call print_memory_usage() call print_memory_usage()
print *, irp_here, ': Allocation failed' print *, irp_here, ': Allocation failed'
stop -1 stop -1
endif endif
!$OMP PARALLEL DO PRIVATE(k)
!$OMP PARALLEL DO PRIVATE(k,j)
do k=1,rank do k=1,rank
call dcopy(ndim, L(1,k), 1, cholesky_ao(1,1,k), 1) do j=1,ao_num
cholesky_ao(1:ao_num,j,k) = L((j-1_8)*ao_num+1_8:1_8*j*ao_num,k)
enddo
enddo enddo
!$OMP END PARALLEL DO !$OMP END PARALLEL DO
deallocate(L)
call munmap( (/ ndim8, rank_max /), 8, fd(1), c_pointer(1) )
cholesky_ao_num = rank cholesky_ao_num = rank
print *, '============ ============='
print *, ''
if (write_ao_cholesky) then if (write_ao_cholesky) then
print *, 'Writing Cholesky vectors to disk...' print *, 'Writing Cholesky AO vectors to disk...'
iunit = getUnitAndOpen(trim(ezfio_work_dir)//'cholesky_ao', 'W') iunit = getUnitAndOpen(trim(ezfio_work_dir)//'cholesky_ao', 'W')
write(iunit) rank write(iunit) rank
write(iunit) cholesky_ao write(iunit) cholesky_ao
@ -409,6 +474,9 @@ END_PROVIDER
print *, 'Rank : ', cholesky_ao_num, '(', 100.d0*dble(cholesky_ao_num)/dble(ao_num*ao_num), ' %)' print *, 'Rank : ', cholesky_ao_num, '(', 100.d0*dble(cholesky_ao_num)/dble(ao_num*ao_num), ' %)'
print *, '' print *, ''
call wall_time(wall1)
print*,'Time to provide AO cholesky vectors = ',(wall1-wall0)/60.d0, ' min'
END_PROVIDER END_PROVIDER

4
src/ao_two_e_ints/two_e_integrals.irp.f
View File

@ -460,8 +460,8 @@ BEGIN_PROVIDER [ double precision, ao_two_e_integral_schwartz, (ao_num, ao_num)
!$OMP PARALLEL DO PRIVATE(i,k) & !$OMP PARALLEL DO PRIVATE(i,k) &
!$OMP DEFAULT(NONE) & !$OMP DEFAULT(NONE) &
!$OMP SHARED (ao_num,ao_two_e_integral_schwartz) & !$OMP SHARED (ao_num,ao_two_e_integral_schwartz) &
!$OMP SCHEDULE(guided) !$OMP SCHEDULE(dynamic)
do i=1,ao_num do i=ao_num,1,-1
do k=1,i do k=1,i
ao_two_e_integral_schwartz(i,k) = dsqrt(ao_two_e_integral(i,i,k,k)) ao_two_e_integral_schwartz(i,k) = dsqrt(ao_two_e_integral(i,i,k,k))
ao_two_e_integral_schwartz(k,i) = ao_two_e_integral_schwartz(i,k) ao_two_e_integral_schwartz(k,i) = ao_two_e_integral_schwartz(i,k)

16
src/ccsd/ccsd_space_orb_sub.irp.f
View File

@ -18,6 +18,8 @@ subroutine run_ccsd_space_orb
integer(bit_kind) :: det(N_int,2) integer(bit_kind) :: det(N_int,2)
integer :: nO, nV, nOa, nVa integer :: nO, nV, nOa, nVa
call set_multiple_levels_omp(.False.)
if (do_ao_cholesky) then if (do_ao_cholesky) then
PROVIDE cholesky_mo_transp PROVIDE cholesky_mo_transp
FREE cholesky_ao FREE cholesky_ao
@ -192,7 +194,7 @@ subroutine run_ccsd_space_orb
deallocate(H_vv,H_oo,H_vo,r1,r2,tau) deallocate(H_vv,H_oo,H_vo,r1,r2,tau)
! CCSD(T) ! CCSD(T)
double precision :: e_t double precision :: e_t, e_t_err
e_t = 0.d0 e_t = 0.d0
if (cc_par_t .and. elec_alpha_num + elec_beta_num > 2) then if (cc_par_t .and. elec_alpha_num + elec_beta_num > 2) then
@ -210,22 +212,24 @@ subroutine run_ccsd_space_orb
!print*,'' !print*,''
! New ! New
e_t = uncorr_energy + energy ! For print in (T) call
e_t_err = 0.d0
print*,'Computing (T) correction...' print*,'Computing (T) correction...'
call wall_time(ta) call wall_time(ta)
! call ccsd_par_t_space_v3(nO,nV,t1,t2,cc_space_f_o,cc_space_f_v & ! call ccsd_par_t_space_v3(nO,nV,t1,t2,cc_space_f_o,cc_space_f_v &
! ,cc_space_v_vvvo,cc_space_v_vvoo,cc_space_v_vooo,e_t) ! ,cc_space_v_vvvo,cc_space_v_vvoo,cc_space_v_vooo,e_t)
e_t = uncorr_energy + energy ! For print in next call
call ccsd_par_t_space_stoch(nO,nV,t1,t2,cc_space_f_o,cc_space_f_v & call ccsd_par_t_space_stoch(nO,nV,t1,t2,cc_space_f_o,cc_space_f_v &
,cc_space_v_vvvo,cc_space_v_vvoo,cc_space_v_vooo,e_t) ,cc_space_v_vvvo,cc_space_v_vvoo,cc_space_v_vooo,e_t, e_t_err)
call wall_time(tb) call wall_time(tb)
print*,'Time: ',tb-ta, ' s' print*,'Time: ',tb-ta, ' s'
print*,'' print*,''
write(*,'(A15,F18.12,A3)') ' E(CCSD(T)) = ', uncorr_energy + energy + e_t, ' Ha' write(*,'(A15,F18.12,A7,F18.12)') ' E(CCSD(T)) = ', uncorr_energy + energy + e_t, ' Ha +/- ', e_t_err
write(*,'(A15,F18.12,A3)') ' E(T) = ', e_t, ' Ha' write(*,'(A15,F18.12,A7,F18.12)') ' E(T) = ', e_t, ' Ha +/- ', e_t_err
write(*,'(A15,F18.12,A3)') ' Correlation = ', energy + e_t, ' Ha' write(*,'(A15,F18.12,A7,F18.12)') ' Correlation = ', energy + e_t, ' Ha +/- ', e_t_err
print*,'' print*,''
endif endif

16
src/ccsd/ccsd_t_space_orb_stoch.irp.f
View File

@ -1,5 +1,5 @@
! Main ! Main
subroutine ccsd_par_t_space_stoch(nO,nV,t1,t2,f_o,f_v,v_vvvo,v_vvoo,v_vooo,energy) subroutine ccsd_par_t_space_stoch(nO,nV,t1,t2,f_o,f_v,v_vvvo,v_vvoo,v_vooo,energy,t_error)
implicit none implicit none
@ -7,7 +7,7 @@ subroutine ccsd_par_t_space_stoch(nO,nV,t1,t2,f_o,f_v,v_vvvo,v_vvoo,v_vooo,energ
double precision, intent(in) :: t1(nO,nV), f_o(nO), f_v(nV) double precision, intent(in) :: t1(nO,nV), f_o(nO), f_v(nV)
double precision, intent(in) :: t2(nO,nO,nV,nV) double precision, intent(in) :: t2(nO,nO,nV,nV)
double precision, intent(in) :: v_vvvo(nV,nV,nV,nO), v_vvoo(nV,nV,nO,nO), v_vooo(nV,nO,nO,nO) double precision, intent(in) :: v_vvvo(nV,nV,nV,nO), v_vvoo(nV,nV,nO,nO), v_vooo(nV,nO,nO,nO)
double precision, intent(inout) :: energy double precision, intent(inout) :: energy, t_error
double precision, allocatable :: X_vovv(:,:,:,:), X_ooov(:,:,:,:), X_oovv(:,:,:,:) double precision, allocatable :: X_vovv(:,:,:,:), X_ooov(:,:,:,:), X_oovv(:,:,:,:)
double precision, allocatable :: T_voov(:,:,:,:), T_oovv(:,:,:,:) double precision, allocatable :: T_voov(:,:,:,:), T_oovv(:,:,:,:)
@ -220,6 +220,8 @@ subroutine ccsd_par_t_space_stoch(nO,nV,t1,t2,f_o,f_v,v_vvvo,v_vvoo,v_vooo,energ
call set_multiple_levels_omp(.False.) call set_multiple_levels_omp(.False.)
call wall_time(t00) call wall_time(t00)
imin = 1_8 imin = 1_8
t_error = huge(1.d0)
!$OMP PARALLEL & !$OMP PARALLEL &
!$OMP PRIVATE(ieta,eta,a,b,c,kiter,isample) & !$OMP PRIVATE(ieta,eta,a,b,c,kiter,isample) &
!$OMP DEFAULT(SHARED) NUM_THREADS(nthreads_pt2) !$OMP DEFAULT(SHARED) NUM_THREADS(nthreads_pt2)
@ -328,15 +330,23 @@ subroutine ccsd_par_t_space_stoch(nO,nV,t1,t2,f_o,f_v,v_vvvo,v_vvoo,v_vooo,energ
if (norm > 0.d0) then if (norm > 0.d0) then
energy_stoch = ET / norm energy_stoch = ET / norm
variance = ET2 / norm - energy_stoch*energy_stoch variance = ET2 / norm - energy_stoch*energy_stoch
if (norm > 1.d0) then
t_error = dsqrt(variance/(norm-1.d0))
else
t_error = dsqrt(variance)
endif
endif endif
energy = energy_det + energy_stoch energy = energy_det + energy_stoch
print '('' '',F20.8, '' '', ES12.4,'' '', F8.2,'' '')', eccsd+energy, dsqrt(variance/(norm-1.d0)), 100.*real(Ncomputed)/real(Nabc) print '('' '',F20.8, '' '', ES12.4,'' '', F8.2,'' '')', eccsd+energy, t_error, 100.*real(Ncomputed)/real(Nabc)
endif endif
!$OMP END MASTER !$OMP END MASTER
if (t_error < cc_par_t_stop) exit
if (imin > Nabc) exit if (imin > Nabc) exit
enddo enddo
!$OMP TASKWAIT
!$OMP END PARALLEL !$OMP END PARALLEL
print '(A)', ' ======================= ============== ========== ' print '(A)', ' ======================= ============== ========== '

154
src/cipsi/selection.irp.f
View File

@ -178,7 +178,7 @@ subroutine select_singles_and_doubles(i_generator, hole_mask, particle_mask, foc
integer(bit_kind), allocatable :: minilist(:, :, :), fullminilist(:, :, :) integer(bit_kind), allocatable :: minilist(:, :, :), fullminilist(:, :, :)
logical, allocatable :: banned(:,:,:), bannedOrb(:,:) logical, allocatable :: banned(:,:,:), bannedOrb(:,:)
double precision, allocatable :: coef_fullminilist_rev(:,:) double precision, allocatable :: coef_fullminilist_rev(:,:)
double precision, allocatable :: mat(:,:,:) double precision, allocatable :: mat(:,:,:), hij_cache(:,:,:)
PROVIDE psi_bilinear_matrix_columns_loc psi_det_alpha_unique psi_det_beta_unique PROVIDE psi_bilinear_matrix_columns_loc psi_det_alpha_unique psi_det_beta_unique
@ -205,7 +205,7 @@ subroutine select_singles_and_doubles(i_generator, hole_mask, particle_mask, foc
! Removed to avoid introducing determinants already presents in the wf ! Removed to avoid introducing determinants already presents in the wf
!double precision, parameter :: norm_thr = 1.d-16 !double precision, parameter :: norm_thr = 1.d-16
allocate (indices(N_det), & allocate (indices(N_det), hij_cache(mo_num,mo_num,2), &
exc_degree(max(N_det_alpha_unique,N_det_beta_unique))) exc_degree(max(N_det_alpha_unique,N_det_beta_unique)))
! Pre-compute excitation degrees wrt alpha determinants ! Pre-compute excitation degrees wrt alpha determinants
@ -511,11 +511,15 @@ subroutine select_singles_and_doubles(i_generator, hole_mask, particle_mask, foc
maskInd = maskInd + 1 maskInd = maskInd + 1
if(mod(maskInd, csubset) == (subset-1)) then if(mod(maskInd, csubset) == (subset-1)) then
call get_mo_two_e_integrals_ij(h2,h1,mo_num,hij_cache(1,1,1),mo_integrals_map)
if (sp /= 3) then ! AA or BB
call get_mo_two_e_integrals_ij(h1,h2,mo_num,hij_cache(1,1,2),mo_integrals_map)
endif
call spot_isinwf(mask, fullminilist, i_generator, fullinteresting(0), banned, fullMatch, fullinteresting) call spot_isinwf(mask, fullminilist, i_generator, fullinteresting(0), banned, fullMatch, fullinteresting)
if(fullMatch) cycle if(fullMatch) cycle
call splash_pq(mask, sp, minilist, i_generator, interesting(0), bannedOrb, banned, mat, interesting) call splash_pq(mask, sp, minilist, i_generator, interesting(0), bannedOrb, banned, mat, interesting, hij_cache)
call fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_diag_tmp, E0, pt2_data, mat, buf) call fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_diag_tmp, E0, pt2_data, mat, buf)
end if end if
@ -531,7 +535,7 @@ subroutine select_singles_and_doubles(i_generator, hole_mask, particle_mask, foc
enddo enddo
enddo enddo
deallocate(preinteresting, prefullinteresting, interesting, fullinteresting) deallocate(preinteresting, prefullinteresting, interesting, fullinteresting)
deallocate(banned, bannedOrb,mat) deallocate(banned, bannedOrb, mat, hij_cache)
end subroutine end subroutine
BEGIN_TEMPLATE BEGIN_TEMPLATE
@ -556,7 +560,7 @@ subroutine fill_buffer_$DOUBLE(i_generator, sp, h1, h2, bannedOrb, banned, fock_
double precision, external :: diag_H_mat_elem_fock double precision, external :: diag_H_mat_elem_fock
double precision :: E_shift double precision :: E_shift
double precision :: s_weight(N_states,N_states) double precision :: s_weight(N_states,N_states)
PROVIDE dominant_dets_of_cfgs N_dominant_dets_of_cfgs PROVIDE dominant_dets_of_cfgs N_dominant_dets_of_cfgs thresh_sym excitation_ref hf_bitmask elec_alpha_num
do jstate=1,N_states do jstate=1,N_states
do istate=1,N_states do istate=1,N_states
s_weight(istate,jstate) = dsqrt(selection_weight(istate)*selection_weight(jstate)) s_weight(istate,jstate) = dsqrt(selection_weight(istate)*selection_weight(jstate))
@ -742,7 +746,7 @@ subroutine fill_buffer_$DOUBLE(i_generator, sp, h1, h2, bannedOrb, banned, fock_
do istate=1,N_states do istate=1,N_states
delta_E = E0(istate) - Hii + E_shift delta_E = E0(istate) - Hii + E_shift
alpha_h_psi = mat(istate, p1, p2) alpha_h_psi = mat(istate, p1, p2)
if (alpha_h_psi == 0.d0) cycle if (dabs(alpha_h_psi) < mo_integrals_threshold) cycle
val = alpha_h_psi + alpha_h_psi val = alpha_h_psi + alpha_h_psi
tmp = dsqrt(delta_E * delta_E + val * val) tmp = dsqrt(delta_E * delta_E + val * val)
@ -914,7 +918,7 @@ single ; do p1=1,mo_num ; enddo ; p2=1 ; ; .False. ;;
END_TEMPLATE END_TEMPLATE
subroutine splash_pq(mask, sp, det, i_gen, N_sel, bannedOrb, banned, mat, interesting) subroutine splash_pq(mask, sp, det, i_gen, N_sel, bannedOrb, banned, mat, interesting, hij_cache)
use bitmasks use bitmasks
implicit none implicit none
BEGIN_DOC BEGIN_DOC
@ -926,6 +930,7 @@ subroutine splash_pq(mask, sp, det, i_gen, N_sel, bannedOrb, banned, mat, intere
integer, intent(in) :: sp, i_gen, N_sel integer, intent(in) :: sp, i_gen, N_sel
integer, intent(in) :: interesting(0:N_sel) integer, intent(in) :: interesting(0:N_sel)
integer(bit_kind),intent(in) :: mask(N_int, 2), det(N_int, 2, N_sel) integer(bit_kind),intent(in) :: mask(N_int, 2), det(N_int, 2, N_sel)
double precision, intent(in) :: hij_cache(mo_num, mo_num, 2)
logical, intent(inout) :: bannedOrb(mo_num, 2), banned(mo_num, mo_num, 2) logical, intent(inout) :: bannedOrb(mo_num, 2), banned(mo_num, mo_num, 2)
double precision, intent(inout) :: mat(N_states, mo_num, mo_num) double precision, intent(inout) :: mat(N_states, mo_num, mo_num)
@ -995,18 +1000,36 @@ subroutine splash_pq(mask, sp, det, i_gen, N_sel, bannedOrb, banned, mat, intere
if(nt == 4) then if(nt == 4) then
call get_d2(det(1,1,i), phasemask, bannedOrb, banned, mat, mask, h, p, sp, psi_selectors_coef_transp(1, interesting(i))) call get_d2(det(1,1,i), phasemask, bannedOrb, banned, mat, mask, h, p, sp, psi_selectors_coef_transp(1, interesting(i)))
else if(nt == 3) then else if(nt == 3) then
call get_d1(det(1,1,i), phasemask, bannedOrb, banned, mat, mask, h, p, sp, psi_selectors_coef_transp(1, interesting(i))) call get_d1(det(1,1,i), phasemask, bannedOrb, banned, mat, mask, h, p, sp, psi_selectors_coef_transp(1, interesting(i)), hij_cache)
else else
call get_d0(det(1,1,i), phasemask, bannedOrb, banned, mat, mask, h, p, sp, psi_selectors_coef_transp(1, interesting(i))) call get_d0(det(1,1,i), phasemask, bannedOrb, banned, mat, mask, h, p, sp, psi_selectors_coef_transp(1, interesting(i)), hij_cache)
end if end if
else if(nt == 4) then else if(nt == 4) then
call bitstring_to_list_in_selection(mobMask(1,1), p(1,1), p(0,1), N_int) call bitstring_to_list_in_selection(mobMask(1,1), p(1,1), p(0,1), N_int)
call bitstring_to_list_in_selection(mobMask(1,2), p(1,2), p(0,2), N_int) call bitstring_to_list_in_selection(mobMask(1,2), p(1,2), p(0,2), N_int)
call past_d2(banned, p, sp) if(sp == 3) then
do j=1,p(0,2)
do ii=1,p(0,1)
banned(p(ii,1), p(j,2),1) = .true.
end do
end do
else
do ii=1,p(0, sp)
do j=1,ii-1
banned(p(j,sp), p(ii,sp),1) = .true.
banned(p(ii,sp), p(j,sp),1) = .true.
end do
end do
end if
else if(nt == 3) then else if(nt == 3) then
call bitstring_to_list_in_selection(mobMask(1,1), p(1,1), p(0,1), N_int) call bitstring_to_list_in_selection(mobMask(1,1), p(1,1), p(0,1), N_int)
call bitstring_to_list_in_selection(mobMask(1,2), p(1,2), p(0,2), N_int) call bitstring_to_list_in_selection(mobMask(1,2), p(1,2), p(0,2), N_int)
call past_d1(bannedOrb, p) do ii = 1, p(0, 1)
bannedOrb(p(ii, 1), 1) = .true.
end do
do ii = 1, p(0, 2)
bannedOrb(p(ii, 2), 2) = .true.
end do
end if end if
end do end do
@ -1037,6 +1060,7 @@ subroutine get_d2(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs)
integer :: bant integer :: bant
bant = 1 bant = 1
PROVIDE mo_integrals_threshold
tip = p(0,1) * p(0,2) tip = p(0,1) * p(0,2)
ma = sp ma = sp
@ -1062,7 +1086,7 @@ subroutine get_d2(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs)
p2 = p(i2, ma) p2 = p(i2, ma)
hij = mo_two_e_integral(p1, p2, h1, h2) - mo_two_e_integral(p2, p1, h1, h2) hij = mo_two_e_integral(p1, p2, h1, h2) - mo_two_e_integral(p2, p1, h1, h2)
if (hij == 0.d0) cycle if (dabs(hij) < mo_integrals_threshold) cycle
hij = hij * get_phase_bi(phasemask, ma, ma, h1, p1, h2, p2, N_int) hij = hij * get_phase_bi(phasemask, ma, ma, h1, p1, h2, p2, N_int)
@ -1092,7 +1116,7 @@ subroutine get_d2(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs)
p1 = p(turn2(i), 1) p1 = p(turn2(i), 1)
hij = mo_two_e_integral(p1, p2, h1, h2) hij = mo_two_e_integral(p1, p2, h1, h2)
if (hij /= 0.d0) then if (dabs(hij) > mo_integrals_threshold) then
hij = hij * get_phase_bi(phasemask, 1, 2, h1, p1, h2, p2, N_int) hij = hij * get_phase_bi(phasemask, 1, 2, h1, p1, h2, p2, N_int)
!DIR$ LOOP COUNT AVG(4) !DIR$ LOOP COUNT AVG(4)
do k=1,N_states do k=1,N_states
@ -1120,7 +1144,7 @@ subroutine get_d2(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs)
p1 = p(i1, ma) p1 = p(i1, ma)
p2 = p(i2, ma) p2 = p(i2, ma)
hij = mo_two_e_integral(p1, p2, h1, h2) - mo_two_e_integral(p2,p1, h1, h2) hij = mo_two_e_integral(p1, p2, h1, h2) - mo_two_e_integral(p2,p1, h1, h2)
if (hij == 0.d0) cycle if (dabs(hij) < mo_integrals_threshold) cycle
hij = hij * get_phase_bi(phasemask, ma, ma, h1, p1, h2, p2, N_int) hij = hij * get_phase_bi(phasemask, ma, ma, h1, p1, h2, p2, N_int)
!DIR$ LOOP COUNT AVG(4) !DIR$ LOOP COUNT AVG(4)
@ -1142,7 +1166,7 @@ subroutine get_d2(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs)
p2 = p(i, ma) p2 = p(i, ma)
hij = mo_two_e_integral(p1, p2, h1, h2) hij = mo_two_e_integral(p1, p2, h1, h2)
if (hij == 0.d0) cycle if (dabs(hij) < mo_integrals_threshold) cycle
hij = hij * get_phase_bi(phasemask, mi, ma, h1, p1, h2, p2, N_int) hij = hij * get_phase_bi(phasemask, mi, ma, h1, p1, h2, p2, N_int)
if (puti < putj) then if (puti < putj) then
@ -1179,7 +1203,7 @@ subroutine get_d2(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs)
end end
subroutine get_d1(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs) subroutine get_d1(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs, hij_cache)
use bitmasks use bitmasks
implicit none implicit none
@ -1190,6 +1214,8 @@ subroutine get_d1(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs)
double precision, intent(in) :: coefs(N_states) double precision, intent(in) :: coefs(N_states)
double precision, intent(inout) :: mat(N_states, mo_num, mo_num) double precision, intent(inout) :: mat(N_states, mo_num, mo_num)
integer, intent(in) :: h(0:2,2), p(0:4,2), sp integer, intent(in) :: h(0:2,2), p(0:4,2), sp
double precision, intent(in) :: hij_cache(mo_num, mo_num, 2)
double precision, external :: get_phase_bi, mo_two_e_integral double precision, external :: get_phase_bi, mo_two_e_integral
logical :: ok logical :: ok
@ -1201,12 +1227,12 @@ subroutine get_d1(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs)
integer, parameter :: turn3(2,3) = reshape((/2,3, 1,3, 1,2/), (/2,3/)) integer, parameter :: turn3(2,3) = reshape((/2,3, 1,3, 1,2/), (/2,3/))
integer :: bant integer :: bant
double precision, allocatable :: hij_cache(:,:) double precision, allocatable :: hij_cache1(:,:)
double precision :: hij, tmp_row(N_states, mo_num), tmp_row2(N_states, mo_num) double precision :: hij, tmp_row(N_states, mo_num), tmp_row2(N_states, mo_num)
PROVIDE mo_integrals_map N_int PROVIDE mo_integrals_map N_int
allocate (lbanned(mo_num, 2)) allocate (lbanned(mo_num, 2))
allocate (hij_cache(mo_num,2)) allocate (hij_cache1(mo_num,2))
lbanned = bannedOrb lbanned = bannedOrb
do i=1, p(0,1) do i=1, p(0,1)
@ -1230,13 +1256,11 @@ subroutine get_d1(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs)
p1 = p(1,ma) p1 = p(1,ma)
p2 = p(2,ma) p2 = p(2,ma)
if(.not. bannedOrb(puti, mi)) then if(.not. bannedOrb(puti, mi)) then
call get_mo_two_e_integrals(hfix,p1,p2,mo_num,hij_cache(1,1),mo_integrals_map)
call get_mo_two_e_integrals(hfix,p2,p1,mo_num,hij_cache(1,2),mo_integrals_map)
tmp_row = 0d0 tmp_row = 0d0
do putj=1, hfix-1 do putj=1, hfix-1
if(lbanned(putj, ma)) cycle if(lbanned(putj, ma)) cycle
if(banned(putj, puti,bant)) cycle if(banned(putj, puti,bant)) cycle
hij = hij_cache(putj,1) - hij_cache(putj,2) hij = hij_cache(hfix,putj,1) - hij_cache(putj,hfix,1)
if (hij /= 0.d0) then if (hij /= 0.d0) then
hij = hij * get_phase_bi(phasemask, ma, ma, putj, p1, hfix, p2, N_int) hij = hij * get_phase_bi(phasemask, ma, ma, putj, p1, hfix, p2, N_int)
!DIR$ LOOP COUNT AVG(4) !DIR$ LOOP COUNT AVG(4)
@ -1248,7 +1272,7 @@ subroutine get_d1(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs)
do putj=hfix+1, mo_num do putj=hfix+1, mo_num
if(lbanned(putj, ma)) cycle if(lbanned(putj, ma)) cycle
if(banned(putj, puti,bant)) cycle if(banned(putj, puti,bant)) cycle
hij = hij_cache(putj,2) - hij_cache(putj,1) hij = hij_cache(putj,hfix,1) - hij_cache(hfix,putj,1)
if (hij /= 0.d0) then if (hij /= 0.d0) then
hij = hij * get_phase_bi(phasemask, ma, ma, hfix, p1, putj, p2, N_int) hij = hij * get_phase_bi(phasemask, ma, ma, hfix, p1, putj, p2, N_int)
!DIR$ LOOP COUNT AVG(4) !DIR$ LOOP COUNT AVG(4)
@ -1274,15 +1298,15 @@ subroutine get_d1(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs)
pfix = p(1,mi) pfix = p(1,mi)
tmp_row = 0d0 tmp_row = 0d0
tmp_row2 = 0d0 tmp_row2 = 0d0
call get_mo_two_e_integrals(hfix,pfix,p1,mo_num,hij_cache(1,1),mo_integrals_map) call get_mo_two_e_integrals(hfix,pfix,p1,mo_num,hij_cache1(1,1),mo_integrals_map)
call get_mo_two_e_integrals(hfix,pfix,p2,mo_num,hij_cache(1,2),mo_integrals_map) call get_mo_two_e_integrals(hfix,pfix,p2,mo_num,hij_cache1(1,2),mo_integrals_map)
putj = p1 putj = p1
do puti=1,mo_num !HOT do puti=1,mo_num !HOT
if(lbanned(puti,mi)) cycle if(lbanned(puti,mi)) cycle
!p1 fixed !p1 fixed
putj = p1 putj = p1
if(.not. banned(putj,puti,bant)) then if(.not. banned(putj,puti,bant)) then
hij = hij_cache(puti,2) hij = hij_cache1(puti,2)
if (hij /= 0.d0) then if (hij /= 0.d0) then
hij = hij * get_phase_bi(phasemask, ma, mi, hfix, p2, puti, pfix, N_int) hij = hij * get_phase_bi(phasemask, ma, mi, hfix, p2, puti, pfix, N_int)
!DIR$ LOOP COUNT AVG(4) !DIR$ LOOP COUNT AVG(4)
@ -1296,7 +1320,7 @@ subroutine get_d1(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs)
putj = p2 putj = p2
! do puti=1,mo_num !HOT ! do puti=1,mo_num !HOT
if(.not. banned(putj,puti,bant)) then if(.not. banned(putj,puti,bant)) then
hij = hij_cache(puti,1) hij = hij_cache1(puti,1)
if (hij /= 0.d0) then if (hij /= 0.d0) then
hij = hij * get_phase_bi(phasemask, ma, mi, hfix, p1, puti, pfix, N_int) hij = hij * get_phase_bi(phasemask, ma, mi, hfix, p1, puti, pfix, N_int)
do k=1,N_states do k=1,N_states
@ -1327,13 +1351,13 @@ subroutine get_d1(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs)
puti = p(i, ma) puti = p(i, ma)
p1 = p(turn3(1,i), ma) p1 = p(turn3(1,i), ma)
p2 = p(turn3(2,i), ma) p2 = p(turn3(2,i), ma)
call get_mo_two_e_integrals(hfix,p1,p2,mo_num,hij_cache(1,1),mo_integrals_map) call get_mo_two_e_integrals(hfix,p1,p2,mo_num,hij_cache1(1,1),mo_integrals_map)
call get_mo_two_e_integrals(hfix,p2,p1,mo_num,hij_cache(1,2),mo_integrals_map) call get_mo_two_e_integrals(hfix,p2,p1,mo_num,hij_cache1(1,2),mo_integrals_map)
tmp_row = 0d0 tmp_row = 0d0
do putj=1,hfix-1 do putj=1,hfix-1
if(banned(putj,puti,1)) cycle if(banned(putj,puti,1)) cycle
if(lbanned(putj,ma)) cycle if(lbanned(putj,ma)) cycle
hij = hij_cache(putj,1) - hij_cache(putj,2) hij = hij_cache1(putj,1) - hij_cache1(putj,2)
if (hij /= 0.d0) then if (hij /= 0.d0) then
hij = hij * get_phase_bi(phasemask, ma, ma, putj, p1, hfix, p2, N_int) hij = hij * get_phase_bi(phasemask, ma, ma, putj, p1, hfix, p2, N_int)
tmp_row(:,putj) = tmp_row(:,putj) + hij * coefs(:) tmp_row(:,putj) = tmp_row(:,putj) + hij * coefs(:)
@ -1342,7 +1366,7 @@ subroutine get_d1(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs)
do putj=hfix+1,mo_num do putj=hfix+1,mo_num
if(banned(putj,puti,1)) cycle if(banned(putj,puti,1)) cycle
if(lbanned(putj,ma)) cycle if(lbanned(putj,ma)) cycle
hij = hij_cache(putj,2) - hij_cache(putj,1) hij = hij_cache1(putj,2) - hij_cache1(putj,1)
if (hij /= 0.d0) then if (hij /= 0.d0) then
hij = hij * get_phase_bi(phasemask, ma, ma, hfix, p1, putj, p2, N_int) hij = hij * get_phase_bi(phasemask, ma, ma, hfix, p1, putj, p2, N_int)
tmp_row(:,putj) = tmp_row(:,putj) + hij * coefs(:) tmp_row(:,putj) = tmp_row(:,putj) + hij * coefs(:)
@ -1364,14 +1388,14 @@ subroutine get_d1(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs)
p2 = p(2,ma) p2 = p(2,ma)
tmp_row = 0d0 tmp_row = 0d0
tmp_row2 = 0d0 tmp_row2 = 0d0
call get_mo_two_e_integrals(hfix,p1,pfix,mo_num,hij_cache(1,1),mo_integrals_map) call get_mo_two_e_integrals(hfix,p1,pfix,mo_num,hij_cache1(1,1),mo_integrals_map)
call get_mo_two_e_integrals(hfix,p2,pfix,mo_num,hij_cache(1,2),mo_integrals_map) call get_mo_two_e_integrals(hfix,p2,pfix,mo_num,hij_cache1(1,2),mo_integrals_map)
putj = p2 putj = p2
do puti=1,mo_num do puti=1,mo_num
if(lbanned(puti,ma)) cycle if(lbanned(puti,ma)) cycle
putj = p2 putj = p2
if(.not. banned(puti,putj,1)) then if(.not. banned(puti,putj,1)) then
hij = hij_cache(puti,1) hij = hij_cache1(puti,1)
if (hij /= 0.d0) then if (hij /= 0.d0) then
hij = hij * get_phase_bi(phasemask, mi, ma, hfix, pfix, puti, p1, N_int) hij = hij * get_phase_bi(phasemask, mi, ma, hfix, pfix, puti, p1, N_int)
!DIR$ LOOP COUNT AVG(4) !DIR$ LOOP COUNT AVG(4)
@ -1383,7 +1407,7 @@ subroutine get_d1(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs)
putj = p1 putj = p1
if(.not. banned(puti,putj,1)) then if(.not. banned(puti,putj,1)) then
hij = hij_cache(puti,2) hij = hij_cache1(puti,2)
if (hij /= 0.d0) then if (hij /= 0.d0) then
hij = hij * get_phase_bi(phasemask, mi, ma, hfix, pfix, puti, p2, N_int) hij = hij * get_phase_bi(phasemask, mi, ma, hfix, pfix, puti, p2, N_int)
do k=1,N_states do k=1,N_states
@ -1408,7 +1432,7 @@ subroutine get_d1(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs)
enddo enddo
end if end if
end if end if
deallocate(lbanned,hij_cache) deallocate(lbanned,hij_cache1)
!! MONO !! MONO
if(sp == 3) then if(sp == 3) then
@ -1439,7 +1463,7 @@ end
subroutine get_d0(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs) subroutine get_d0(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs, hij_cache)
use bitmasks use bitmasks
implicit none implicit none
@ -1450,6 +1474,7 @@ subroutine get_d0(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs)
double precision, intent(in) :: coefs(N_states) double precision, intent(in) :: coefs(N_states)
double precision, intent(inout) :: mat(N_states, mo_num, mo_num) double precision, intent(inout) :: mat(N_states, mo_num, mo_num)
integer, intent(in) :: h(0:2,2), p(0:4,2), sp integer, intent(in) :: h(0:2,2), p(0:4,2), sp
double precision, intent(in) :: hij_cache(mo_num, mo_num, 2)
integer :: i, j, k, s, h1, h2, p1, p2, puti, putj integer :: i, j, k, s, h1, h2, p1, p2, puti, putj
double precision :: hij, phase double precision :: hij, phase
@ -1457,16 +1482,14 @@ subroutine get_d0(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs)
logical :: ok logical :: ok
integer, parameter :: bant=1 integer, parameter :: bant=1
double precision, allocatable :: hij_cache1(:), hij_cache2(:)
allocate (hij_cache1(mo_num),hij_cache2(mo_num))
PROVIDE mo_integrals_threshold
if(sp == 3) then ! AB if(sp == 3) then ! AB
h1 = p(1,1) h1 = p(1,1)
h2 = p(1,2) h2 = p(1,2)
do p1=1, mo_num do p1=1, mo_num
if(bannedOrb(p1, 1)) cycle if(bannedOrb(p1, 1)) cycle
call get_mo_two_e_integrals(p1,h2,h1,mo_num,hij_cache1,mo_integrals_map)
do p2=1, mo_num do p2=1, mo_num
if(bannedOrb(p2,2)) cycle if(bannedOrb(p2,2)) cycle
if(banned(p1, p2, bant)) cycle ! rentable? if(banned(p1, p2, bant)) cycle ! rentable?
@ -1475,9 +1498,9 @@ subroutine get_d0(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs)
call i_h_j(gen, det, N_int, hij) call i_h_j(gen, det, N_int, hij)
else else
phase = get_phase_bi(phasemask, 1, 2, h1, p1, h2, p2, N_int) phase = get_phase_bi(phasemask, 1, 2, h1, p1, h2, p2, N_int)
hij = hij_cache1(p2) * phase hij = hij_cache(p2,p1,1) * phase
end if end if
if (hij == 0.d0) cycle if (dabs(hij) < mo_integrals_threshold) cycle
!DIR$ LOOP COUNT AVG(4) !DIR$ LOOP COUNT AVG(4)
do k=1,N_states do k=1,N_states
mat(k, p1, p2) = mat(k, p1, p2) + coefs(k) * hij ! HOTSPOT mat(k, p1, p2) = mat(k, p1, p2) + coefs(k) * hij ! HOTSPOT
@ -1490,18 +1513,16 @@ subroutine get_d0(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs)
p2 = p(2,sp) p2 = p(2,sp)
do puti=1, mo_num do puti=1, mo_num
if (bannedOrb(puti, sp)) cycle if (bannedOrb(puti, sp)) cycle
call get_mo_two_e_integrals(puti,p2,p1,mo_num,hij_cache1,mo_integrals_map)
call get_mo_two_e_integrals(puti,p1,p2,mo_num,hij_cache2,mo_integrals_map)
do putj=puti+1, mo_num do putj=puti+1, mo_num
if(bannedOrb(putj, sp)) cycle if(bannedOrb(putj, sp)) cycle
if(banned(puti, putj, bant)) cycle ! rentable? if(banned(puti, putj, bant)) cycle ! rentable?
if(puti == p1 .or. putj == p2 .or. puti == p2 .or. putj == p1) then if(puti == p1 .or. putj == p2 .or. puti == p2 .or. putj == p1) then
call apply_particles(mask, sp,puti,sp,putj, det, ok, N_int) call apply_particles(mask, sp,puti,sp,putj, det, ok, N_int)
call i_h_j(gen, det, N_int, hij) call i_h_j(gen, det, N_int, hij)
if (hij == 0.d0) cycle if (dabs(hij) < mo_integrals_threshold) cycle
else else
hij = hij_cache1(putj) - hij_cache2(putj) hij = hij_cache(putj,puti,1) - hij_cache(putj,puti,2)
if (hij == 0.d0) cycle if (dabs(hij) < mo_integrals_threshold) cycle
hij = hij * get_phase_bi(phasemask, sp, sp, puti, p1 , putj, p2, N_int) hij = hij * get_phase_bi(phasemask, sp, sp, puti, p1 , putj, p2, N_int)
end if end if
!DIR$ LOOP COUNT AVG(4) !DIR$ LOOP COUNT AVG(4)
@ -1512,50 +1533,9 @@ subroutine get_d0(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs)
end do end do
end if end if
deallocate(hij_cache1,hij_cache2)
end end
subroutine past_d1(bannedOrb, p)
use bitmasks
implicit none
logical, intent(inout) :: bannedOrb(mo_num, 2)
integer, intent(in) :: p(0:4, 2)
integer :: i,s
do s = 1, 2
do i = 1, p(0, s)
bannedOrb(p(i, s), s) = .true.
end do
end do
end
subroutine past_d2(banned, p, sp)
use bitmasks
implicit none
logical, intent(inout) :: banned(mo_num, mo_num)
integer, intent(in) :: p(0:4, 2), sp
integer :: i,j
if(sp == 3) then
do j=1,p(0,2)
do i=1,p(0,1)
banned(p(i,1), p(j,2)) = .true.
end do
end do
else
do i=1,p(0, sp)
do j=1,i-1
banned(p(j,sp), p(i,sp)) = .true.
banned(p(i,sp), p(j,sp)) = .true.
end do
end do
end if
end
subroutine spot_isinwf(mask, det, i_gen, N, banned, fullMatch, interesting) subroutine spot_isinwf(mask, det, i_gen, N, banned, fullMatch, interesting)
use bitmasks use bitmasks
implicit none implicit none

1
src/dav_general_mat/dav_diag_dressed_ext_rout.irp.f
View File

@ -1,6 +1,5 @@
subroutine davidson_general_ext_rout_diag_dressed(u_in,H_jj,Dress_jj,energies,sze,N_st,N_st_diag_in,converged,hcalc) subroutine davidson_general_ext_rout_diag_dressed(u_in,H_jj,Dress_jj,energies,sze,N_st,N_st_diag_in,converged,hcalc)
use mmap_module
implicit none implicit none
BEGIN_DOC BEGIN_DOC
! Generic Davidson diagonalization with ONE DIAGONAL DRESSING OPERATOR ! Generic Davidson diagonalization with ONE DIAGONAL DRESSING OPERATOR

2
src/dav_general_mat/dav_diag_dressed_ext_rout_nonsym_B1space.irp.f
View File

@ -3,8 +3,6 @@
subroutine davidson_general_diag_dressed_ext_rout_nonsym_b1space(u_in, H_jj, Dress_jj,energies, sze, N_st, N_st_diag_in, converged, hcalc) subroutine davidson_general_diag_dressed_ext_rout_nonsym_b1space(u_in, H_jj, Dress_jj,energies, sze, N_st, N_st_diag_in, converged, hcalc)
use mmap_module
BEGIN_DOC BEGIN_DOC
! Generic modified-Davidson diagonalization ! Generic modified-Davidson diagonalization
! !

1
src/dav_general_mat/dav_double_dress_ext_rout.irp.f
View File

@ -1,5 +1,4 @@
subroutine dav_double_dressed(u_in,H_jj,Dress_jj,Dressing_vec,idx_dress,energies,sze,N_st,N_st_diag,converged,hcalc) subroutine dav_double_dressed(u_in,H_jj,Dress_jj,Dressing_vec,idx_dress,energies,sze,N_st,N_st_diag,converged,hcalc)
use mmap_module
BEGIN_DOC BEGIN_DOC
! Generic Davidson diagonalization with TWO DRESSING VECTORS ! Generic Davidson diagonalization with TWO DRESSING VECTORS
! !

1
src/dav_general_mat/dav_dressed_ext_rout.irp.f
View File

@ -1,5 +1,4 @@
subroutine davidson_general_ext_rout_dressed(u_in,H_jj,energies,sze,N_st,N_st_diag,dressing_state,dressing_vec,idress,converged,hcalc) subroutine davidson_general_ext_rout_dressed(u_in,H_jj,energies,sze,N_st,N_st_diag,dressing_state,dressing_vec,idress,converged,hcalc)
use mmap_module
implicit none implicit none
BEGIN_DOC BEGIN_DOC
! Davidson diagonalization. ! Davidson diagonalization.

1
src/dav_general_mat/dav_ext_rout.irp.f
View File

@ -1,6 +1,5 @@
subroutine davidson_general_ext_rout(u_in,H_jj,energies,sze,N_st,N_st_diag_in,converged,hcalc) subroutine davidson_general_ext_rout(u_in,H_jj,energies,sze,N_st,N_st_diag_in,converged,hcalc)
use mmap_module
implicit none implicit none
BEGIN_DOC BEGIN_DOC
! Generic Davidson diagonalization ! Generic Davidson diagonalization

2
src/dav_general_mat/dav_ext_rout_nonsym_B1space.irp.f
View File

@ -3,8 +3,6 @@
subroutine davidson_general_ext_rout_nonsym_b1space(u_in, H_jj, energies, sze, N_st, N_st_diag_in, converged, hcalc) subroutine davidson_general_ext_rout_nonsym_b1space(u_in, H_jj, energies, sze, N_st, N_st_diag_in, converged, hcalc)
use mmap_module
BEGIN_DOC BEGIN_DOC
! Generic modified-Davidson diagonalization ! Generic modified-Davidson diagonalization
! !

6
src/dav_general_mat/dav_general.irp.f
View File

@ -1,6 +1,6 @@
subroutine davidson_general(u_in,H_jj,energies,dim_in,sze,N_st,N_st_diag_in,converged,h_mat) subroutine davidson_general(u_in,H_jj,energies,dim_in,sze,N_st,N_st_diag_in,converged,h_mat)
use mmap_module ! use mmap_module
implicit none implicit none
BEGIN_DOC BEGIN_DOC
! Davidson diagonalization with specific diagonal elements of the H matrix ! Davidson diagonalization with specific diagonal elements of the H matrix
@ -160,9 +160,9 @@ subroutine davidson_general(u_in,H_jj,energies,dim_in,sze,N_st,N_st_diag_in,conv
! type(c_ptr) :: ptr_w, ptr_s ! type(c_ptr) :: ptr_w, ptr_s
! integer :: fd_s, fd_w ! integer :: fd_s, fd_w
! call mmap(trim(ezfio_work_dir)//'davidson_w', (/int(sze,8),int(N_st_diag*itermax,8)/),& ! call mmap(trim(ezfio_work_dir)//'davidson_w', (/int(sze,8),int(N_st_diag*itermax,8)/),&
! 8, fd_w, .False., ptr_w) ! 8, fd_w, .False., .True., ptr_w)
! call mmap(trim(ezfio_work_dir)//'davidson_s', (/int(sze,8),int(N_st_diag*itermax,8)/),& ! call mmap(trim(ezfio_work_dir)//'davidson_s', (/int(sze,8),int(N_st_diag*itermax,8)/),&
! 4, fd_s, .False., ptr_s) ! 4, fd_s, .False., .True., ptr_s)
! call c_f_pointer(ptr_w, w, (/sze,N_st_diag*itermax/)) ! call c_f_pointer(ptr_w, w, (/sze,N_st_diag*itermax/))
! call c_f_pointer(ptr_s, s, (/sze,N_st_diag*itermax/)) ! call c_f_pointer(ptr_s, s, (/sze,N_st_diag*itermax/))
! else ! else

2
src/davidson/diagonalization_h_dressed.irp.f
View File

@ -228,7 +228,7 @@ subroutine davidson_diag_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,N_st_dia
type(c_ptr) :: ptr_w, ptr_s type(c_ptr) :: ptr_w, ptr_s
integer :: fd_s, fd_w integer :: fd_s, fd_w
call mmap(trim(ezfio_work_dir)//'davidson_w', (/int(sze,8),int(N_st_diag*itermax,8)/),& call mmap(trim(ezfio_work_dir)//'davidson_w', (/int(sze,8),int(N_st_diag*itermax,8)/),&
8, fd_w, .False., ptr_w) 8, fd_w, .False., .True., ptr_w)
call c_f_pointer(ptr_w, w, (/sze,N_st_diag*itermax/)) call c_f_pointer(ptr_w, w, (/sze,N_st_diag*itermax/))
else else
allocate(W(sze,N_st_diag*itermax)) allocate(W(sze,N_st_diag*itermax))

2
src/davidson/diagonalization_hcsf_dressed.irp.f
View File

@ -229,7 +229,7 @@ subroutine davidson_diag_csf_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,sze_csf,N
type(c_ptr) :: ptr_w, ptr_s type(c_ptr) :: ptr_w, ptr_s
integer :: fd_s, fd_w integer :: fd_s, fd_w
call mmap(trim(ezfio_work_dir)//'davidson_w', (/int(sze,8),int(N_st_diag*itermax,8)/),& call mmap(trim(ezfio_work_dir)//'davidson_w', (/int(sze,8),int(N_st_diag*itermax,8)/),&
8, fd_w, .False., ptr_w) 8, fd_w, .False., .True., ptr_w)
call c_f_pointer(ptr_w, W_csf, (/sze_csf,N_st_diag*itermax/)) call c_f_pointer(ptr_w, W_csf, (/sze_csf,N_st_diag*itermax/))
else else
allocate(W(sze,N_st_diag),W_csf(sze_csf,N_st_diag*itermax)) allocate(W(sze,N_st_diag),W_csf(sze_csf,N_st_diag*itermax))

4
src/davidson/diagonalization_hs2_dressed.irp.f
View File

@ -270,9 +270,9 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_
type(c_ptr) :: ptr_w, ptr_s type(c_ptr) :: ptr_w, ptr_s
integer :: fd_s, fd_w integer :: fd_s, fd_w
call mmap(trim(ezfio_work_dir)//'davidson_w', (/int(sze,8),int(N_st_diag*itermax,8)/),& call mmap(trim(ezfio_work_dir)//'davidson_w', (/int(sze,8),int(N_st_diag*itermax,8)/),&
8, fd_w, .False., ptr_w) 8, fd_w, .False., .True., ptr_w)
call mmap(trim(ezfio_work_dir)//'davidson_s', (/int(sze,8),int(N_st_diag*itermax,8)/),& call mmap(trim(ezfio_work_dir)//'davidson_s', (/int(sze,8),int(N_st_diag*itermax,8)/),&
4, fd_s, .False., ptr_s) 4, fd_s, .False., .True., ptr_s)
call c_f_pointer(ptr_w, w, (/sze,N_st_diag*itermax/)) call c_f_pointer(ptr_w, w, (/sze,N_st_diag*itermax/))
call c_f_pointer(ptr_s, s, (/sze,N_st_diag*itermax/)) call c_f_pointer(ptr_s, s, (/sze,N_st_diag*itermax/))
else else

2
src/davidson/diagonalization_nonsym_h_dressed.irp.f
View File

@ -251,7 +251,7 @@ subroutine davidson_diag_nonsym_hjj(dets_in, u_in, H_jj, energies, dim_in, sze,
type(c_ptr) :: ptr_w, ptr_s type(c_ptr) :: ptr_w, ptr_s
integer :: fd_s, fd_w integer :: fd_s, fd_w
call mmap(trim(ezfio_work_dir)//'davidson_w', (/int(sze,8),int(N_st_diag*itermax,8)/),& call mmap(trim(ezfio_work_dir)//'davidson_w', (/int(sze,8),int(N_st_diag*itermax,8)/),&
8, fd_w, .False., ptr_w) 8, fd_w, .False., .True., ptr_w)
call c_f_pointer(ptr_w, w, (/sze,N_st_diag*itermax/)) call c_f_pointer(ptr_w, w, (/sze,N_st_diag*itermax/))
else else
allocate(W(sze,N_st_diag*itermax)) allocate(W(sze,N_st_diag*itermax))

192
src/determinants/slater_rules_general.irp.f Normal file
View File

@ -0,0 +1,192 @@
subroutine get_excitation_general(key_i,key_j, Nint,degree_array,holes_array, particles_array,phase)
use bitmasks
BEGIN_DOC
! returns the array, for each spin, of holes/particles between key_i and key_j
!
! with the following convention: a^+_{particle} a_{hole}|key_i> = |key_j>
END_DOC
include 'utils/constants.include.F'
implicit none
integer, intent(in) :: Nint
integer(bit_kind), intent(in) :: key_j(Nint,2),key_i(Nint,2)
integer, intent(out) :: holes_array(100,2),particles_array(100,2),degree_array(2)
double precision, intent(out) :: phase
integer :: ispin,k,i,pos
integer(bit_kind) :: key_hole, key_particle
integer(bit_kind) :: xorvec(N_int_max,2)
holes_array = -1
particles_array = -1
degree_array = 0
do i = 1, N_int
xorvec(i,1) = xor( key_i(i,1), key_j(i,1))
xorvec(i,2) = xor( key_i(i,2), key_j(i,2))
degree_array(1) += popcnt(xorvec(i,1))
degree_array(2) += popcnt(xorvec(i,2))
enddo
degree_array(1) = shiftr(degree_array(1),1)
degree_array(2) = shiftr(degree_array(2),1)
do ispin = 1, 2
k = 1
!!! GETTING THE HOLES
do i = 1, N_int
key_hole = iand(xorvec(i,ispin),key_i(i,ispin))
do while(key_hole .ne.0_bit_kind)
pos = trailz(key_hole)
holes_array(k,ispin) = 1+ bit_kind_size * (i-1) + pos
key_hole = ibclr(key_hole,pos)
k += 1
if(k .gt.100)then
print*,'WARNING in get_excitation_general'
print*,'More than a 100-th excitation for spin ',ispin
print*,'stoping ...'
stop
endif
enddo
enddo
enddo
do ispin = 1, 2
k = 1
!!! GETTING THE PARTICLES
do i = 1, N_int
key_particle = iand(xor(key_i(i,ispin),key_j(i,ispin)),key_j(i,ispin))
do while(key_particle .ne.0_bit_kind)
pos = trailz(key_particle)
particles_array(k,ispin) = 1+ bit_kind_size * (i-1) + pos
key_particle = ibclr(key_particle,pos)
k += 1
if(k .gt.100)then
print*,'WARNING in get_excitation_general '
print*,'More than a 100-th excitation for spin ',ispin
print*,'stoping ...'
stop
endif
enddo
enddo
enddo
integer :: h,p, i_ok
integer(bit_kind), allocatable :: det_i(:,:),det_ip(:,:)
integer :: exc(0:2,2,2)
double precision :: phase_tmp
allocate(det_i(Nint,2),det_ip(N_int,2))
det_i = key_i
phase = 1.d0
do ispin = 1, 2
do i = 1, degree_array(ispin)
h = holes_array(i,ispin)
p = particles_array(i,ispin)
det_ip = det_i
call do_single_excitation(det_ip,h,p,ispin,i_ok)
if(i_ok == -1)then
print*,'excitation was not possible '
stop
endif
call get_single_excitation(det_i,det_ip,exc,phase_tmp,Nint)
phase *= phase_tmp
det_i = det_ip
enddo
enddo
end
subroutine get_holes_general(key_i, key_j,Nint, holes_array)
use bitmasks
BEGIN_DOC
! returns the array, per spin, of holes between key_i and key_j
!
! with the following convention: a_{hole}|key_i> --> |key_j>
END_DOC
implicit none
integer, intent(in) :: Nint
integer(bit_kind), intent(in) :: key_j(Nint,2),key_i(Nint,2)
integer, intent(out) :: holes_array(100,2)
integer(bit_kind) :: key_hole
integer :: ispin,k,i,pos
holes_array = -1
do ispin = 1, 2
k = 1
do i = 1, N_int
key_hole = iand(xor(key_i(i,ispin),key_j(i,ispin)),key_i(i,ispin))
do while(key_hole .ne.0_bit_kind)
pos = trailz(key_hole)
holes_array(k,ispin) = 1+ bit_kind_size * (i-1) + pos
key_hole = ibclr(key_hole,pos)
k += 1
if(k .gt.100)then
print*,'WARNING in get_holes_general'
print*,'More than a 100-th excitation for spin ',ispin
print*,'stoping ...'
stop
endif
enddo
enddo
enddo
end
subroutine get_particles_general(key_i, key_j,Nint,particles_array)
use bitmasks
BEGIN_DOC
! returns the array, per spin, of particles between key_i and key_j
!
! with the following convention: a^dagger_{particle}|key_i> --> |key_j>
END_DOC
implicit none
integer, intent(in) :: Nint
integer(bit_kind), intent(in) :: key_j(Nint,2),key_i(Nint,2)
integer, intent(out) :: particles_array(100,2)
integer(bit_kind) :: key_particle
integer :: ispin,k,i,pos
particles_array = -1
do ispin = 1, 2
k = 1
do i = 1, N_int
key_particle = iand(xor(key_i(i,ispin),key_j(i,ispin)),key_j(i,ispin))
do while(key_particle .ne.0_bit_kind)
pos = trailz(key_particle)
particles_array(k,ispin) = 1+ bit_kind_size * (i-1) + pos
key_particle = ibclr(key_particle,pos)
k += 1
if(k .gt.100)then
print*,'WARNING in get_holes_general'
print*,'More than a 100-th excitation for spin ',ispin
print*,'Those are the two determinants'
call debug_det(key_i, N_int)
call debug_det(key_j, N_int)
print*,'stoping ...'
stop
endif
enddo
enddo
enddo
end
subroutine get_phase_general(key_i,Nint,degree, holes_array, particles_array,phase)
implicit none
integer, intent(in) :: degree(2), Nint
integer(bit_kind), intent(in) :: key_i(Nint,2)
integer, intent(in) :: holes_array(100,2),particles_array(100,2)
double precision, intent(out) :: phase
integer :: i,ispin,h,p, i_ok
integer(bit_kind), allocatable :: det_i(:,:),det_ip(:,:)
integer :: exc(0:2,2,2)
double precision :: phase_tmp
allocate(det_i(Nint,2),det_ip(N_int,2))
det_i = key_i
phase = 1.d0
do ispin = 1, 2
do i = 1, degree(ispin)
h = holes_array(i,ispin)
p = particles_array(i,ispin)
det_ip = det_i
call do_single_excitation(det_ip,h,p,ispin,i_ok)
if(i_ok == -1)then
print*,'excitation was not possible '
stop
endif
call get_single_excitation(det_i,det_ip,exc,phase_tmp,Nint)
phase *= phase_tmp
det_i = det_ip
enddo
enddo
end

1
src/dft_utils_func/on_top_from_ueg.irp.f
View File

@ -32,7 +32,6 @@ double precision function g0_UEG_mu_inf(rho_a,rho_b)
C = 0.08193d0 C = 0.08193d0
D = -0.01277d0 D = -0.01277d0
E = 0.001859d0 E = 0.001859d0
x = -d2*rs
if (dabs(rho) > 1.d-20) then if (dabs(rho) > 1.d-20) then
rs = (3d0 / (4d0*pi*rho))**(1d0/3d0) ! JT: serious bug fixed 20/03/19 rs = (3d0 / (4d0*pi*rho))**(1d0/3d0) ! JT: serious bug fixed 20/03/19
x = -d2*rs x = -d2*rs

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

@ -48,7 +48,7 @@
integer :: i,j integer :: i,j
do i = 1, n_points_final_grid do i = 1, n_points_final_grid
do j = 1, mo_num do j = 1, mo_num
mos_in_r_array_transp(i,j) = mos_in_r_array(j,i) mos_in_r_array_transp(i,j) = mos_in_r_array_omp(j,i)
enddo enddo
enddo enddo
END_PROVIDER END_PROVIDER

0
src/ezfio_files/01.convert.bats → src/ezfio_files/convert_bats_old
View File

8
src/ezfio_files/get_unit_and_open.irp.f
View File

@ -47,11 +47,13 @@ integer function getUnitAndOpen(f,mode)
endif endif
open(unit=getUnitAndOpen,file=f,status='OLD',action='READ',form='UNFORMATTED') open(unit=getUnitAndOpen,file=f,status='OLD',action='READ',form='UNFORMATTED')
else if (mode.eq.'W') then else if (mode.eq.'W') then
open(unit=getUnitAndOpen,file=new_f,status='UNKNOWN',action='WRITE',form='UNFORMATTED') open(unit=getUnitAndOpen,file=new_f,status='UNKNOWN',action='READWRITE',form='UNFORMATTED')
else if (mode.eq.'A') then
open(unit=getUnitAndOpen,file=new_f,status='UNKNOWN',action='READWRITE',position='APPEND',form='UNFORMATTED')
else if (mode.eq.'w') then else if (mode.eq.'w') then
open(unit=getUnitAndOpen,file=new_f,status='UNKNOWN',action='WRITE',form='FORMATTED') open(unit=getUnitAndOpen,file=new_f,status='UNKNOWN',action='READWRITE',form='FORMATTED')
else if (mode.eq.'a') then else if (mode.eq.'a') then
open(unit=getUnitAndOpen,file=new_f,status='UNKNOWN',action='WRITE',position='APPEND',form='FORMATTED') open(unit=getUnitAndOpen,file=new_f,status='UNKNOWN',action='READWRITE',position='APPEND',form='FORMATTED')
else if (mode.eq.'x') then else if (mode.eq.'x') then
open(unit=getUnitAndOpen,file=new_f,form='FORMATTED') open(unit=getUnitAndOpen,file=new_f,form='FORMATTED')
endif endif

13
src/hartree_fock/10.hf.bats
View File

@ -115,9 +115,6 @@ rm -rf $EZFIO
run hco.ezfio -113.1841002944744 run hco.ezfio -113.1841002944744
} }
@test "HBO" { # 0.805600 1.4543s
run hbo.ezfio -100.018582259096
}
@test "H2S" { # 1.655600 4.21402s @test "H2S" { # 1.655600 4.21402s
run h2s.ezfio -398.6944130421982 run h2s.ezfio -398.6944130421982
@ -127,9 +124,6 @@ rm -rf $EZFIO
run h3coh.ezfio -114.9865030596373 run h3coh.ezfio -114.9865030596373
} }
@test "H2O" { # 1.811100 1.84387s
run h2o.ezfio -0.760270218692179E+02
}
@test "H2O2" { # 2.217000 8.50267s @test "H2O2" { # 2.217000 8.50267s
run h2o2.ezfio -150.7806608469964 run h2o2.ezfio -150.7806608469964
@ -187,13 +181,6 @@ rm -rf $EZFIO
run oh.ezfio -75.42025413469165 run oh.ezfio -75.42025413469165
} }
@test "[Cu(NH3)4]2+" { # 59.610100 4.18766m
[[ -n $TRAVIS ]] && skip
qp set_file cu_nh3_4_2plus.ezfio
qp set scf_utils thresh_scf 1.e-10
run cu_nh3_4_2plus.ezfio -1862.97590358903
}
@test "SO2" { # 71.894900 3.22567m @test "SO2" { # 71.894900 3.22567m
[[ -n $TRAVIS ]] && skip [[ -n $TRAVIS ]] && skip
run so2.ezfio -41.55800401346361 run so2.ezfio -41.55800401346361

19
src/hartree_fock/fock_matrix_hf.irp.f
View File

@ -194,14 +194,25 @@ END_PROVIDER
endif endif
double precision :: rss double precision :: rss, mem0, mem
double precision :: memory_of_double double precision :: memory_of_double
integer :: iblock integer :: iblock
integer, parameter :: block_size = 32 integer :: block_size
call resident_memory(mem0)
block_size = 1024
rss = memory_of_double(2.d0*ao_num*ao_num)
do
mem = mem0 + block_size*rss
if ( (block_size < 2).or.(mem < qp_max_mem) ) exit
block_size = block_size/2
enddo
call check_mem(block_size*rss, irp_here)
rss = memory_of_double(ao_num*ao_num)
call check_mem(2.d0*block_size*rss, irp_here)
allocate(X2(ao_num,ao_num,block_size,2)) allocate(X2(ao_num,ao_num,block_size,2))
allocate(X3(ao_num,block_size,ao_num,2)) allocate(X3(ao_num,block_size,ao_num,2))

18
src/mo_two_e_ints/EZFIO.cfg
View File

@ -1,9 +1,21 @@
[io_mo_cholesky]
type: Disk_access
doc: Read/Write |MO| Cholesky integrals from/to disk [ Write | Read | None ]
interface: ezfio,provider,ocaml
default: None
[io_mo_two_e_integrals] [io_mo_two_e_integrals]
type: Disk_access type: Disk_access
doc: Read/Write |MO| integrals from/to disk [ Write | Read | None ] doc: Read/Write |MO| integrals from/to disk [ Write | Read | None ]
interface: ezfio,provider,ocaml interface: ezfio,provider,ocaml
default: None default: None
[mo_integrals_cache_shift]
type: integer
doc: Adjusts the size of the MO integrals cache. 2: 2KB, 3: 32KB, 4: 512KB, 5: 8MB, 6: 128MB, 7: 2GB, 8: 32GB, 9: 512GB
interface: ezfio, provider, ocaml
default: 7
[mo_integrals_threshold] [mo_integrals_threshold]
type: Threshold type: Threshold
doc: If | <ij|kl> | < `mo_integrals_threshold` then <ij|kl> is zero doc: If | <ij|kl> | < `mo_integrals_threshold` then <ij|kl> is zero
@ -11,12 +23,6 @@ interface: ezfio,provider,ocaml
default: 1.e-15 default: 1.e-15
ezfio_name: threshold_mo ezfio_name: threshold_mo
[no_vvvv_integrals]
type: logical
doc: If `True`, computes all integrals except for the integrals having 3 or 4 virtual indices
interface: ezfio,provider,ocaml
default: false
[io_mo_two_e_integrals_erf] [io_mo_two_e_integrals_erf]
type: Disk_access type: Disk_access
doc: Read/Write MO integrals with the long range interaction from/to disk [ Write | Read | None ] doc: Read/Write MO integrals with the long range interaction from/to disk [ Write | Read | None ]

Some files were not shown because too many files have changed in this diff Show More