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

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

Browse Source 
Conflicts:
	src/tc_bi_ortho/tc_utils.irp.f
This commit is contained in:
Abdallah Ammar 2023-06-02 20:31:34 +02:00
commit 6a7c33aa39
65 changed files with 3054 additions and 1201 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
configure vendored
View File

@ -215,7 +215,6 @@ EOF
cd trexio-${VERSION}
./configure --prefix=\${QP_ROOT} --without-hdf5
make -j 8 && make -j 8 check && make -j 8 install
cp ${QP_ROOT}/include/trexio_f.f90 ${QP_ROOT}/src/ezfio_files
tar -zxvf "\${QP_ROOT}"/external/qp2-dependencies/${ARCHITECTURE}/ninja.tar.gz
mv ninja "\${QP_ROOT}"/bin/
EOF
@ -229,7 +228,6 @@ EOF
cd trexio-${VERSION}
./configure --prefix=\${QP_ROOT}
make -j 8 && make -j 8 check && make -j 8 install
cp ${QP_ROOT}/include/trexio_f.f90 ${QP_ROOT}/src/ezfio_files
EOF

11
ocaml/Input_ao_basis.ml
View File

@ -44,8 +44,12 @@ end = struct
let get_default = Qpackage.get_ezfio_default "ao_basis";;
let read_ao_basis () =
let result =
Ezfio.get_ao_basis_ao_basis ()
|> AO_basis_name.of_string
in
if result <> "None" then
AO_basis_name.of_string result
else failwith "No basis"
;;
let read_ao_num () =
@ -267,7 +271,10 @@ end = struct
|> Ezfio.set_ao_basis_ao_md5 ;
Some result
with
| _ -> (Ezfio.set_ao_basis_ao_md5 "None" ; None)
| _ -> ( "None"
|> Digest.string
|> Digest.to_hex
|> Ezfio.set_ao_basis_ao_md5 ; None)
;;

5
ocaml/Input_mo_basis.ml
View File

@ -56,7 +56,10 @@ end = struct
let read_ao_md5 () =
let ao_md5 =
match (Input_ao_basis.Ao_basis.read ()) with
| None -> failwith "Unable to read AO basis"
| None -> ("None"
|> Digest.string
|> Digest.to_hex
|> MD5.of_string)
| Some result -> Input_ao_basis.Ao_basis.to_md5 result
in
let result =

93
src/trexio/qp_import_trexio.py → scripts/qp_import_trexio.py
View File

@ -13,12 +13,17 @@ Options:
import sys
import os
import trexio
import numpy as np
from functools import reduce
from ezfio import ezfio
from docopt import docopt
try:
import trexio
except ImportError:
print("Error: trexio python module is not found. Try python3 -m pip install trexio")
sys.exit(1)
try:
QP_ROOT = os.environ["QP_ROOT"]
@ -90,14 +95,15 @@ def write_ezfio(trexio_filename, filename):
p = re.compile(r'(\d*)$')
label = [p.sub("", x).capitalize() for x in label]
ezfio.set_nuclei_nucl_label(label)
print("OK")
else:
ezfio.set_nuclei_nucl_num(1)
ezfio.set_nuclei_nucl_charge([0.])
ezfio.set_nuclei_nucl_coord([0.,0.,0.])
ezfio.set_nuclei_nucl_label(["X"])
print("None")
print("OK")
print("Electrons\t...\t", end=' ')
@ -105,12 +111,12 @@ def write_ezfio(trexio_filename, filename):
try:
num_beta = trexio.read_electron_dn_num(trexio_file)
except:
num_beta = sum(charge)//2
num_beta = int(sum(charge))//2
try:
num_alpha = trexio.read_electron_up_num(trexio_file)
except:
num_alpha = sum(charge) - num_beta
num_alpha = int(sum(charge)) - num_beta
if num_alpha == 0:
print("\n\nError: There are zero electrons in the TREXIO file.\n\n")
@ -118,7 +124,7 @@ def write_ezfio(trexio_filename, filename):
ezfio.set_electrons_elec_alpha_num(num_alpha)
ezfio.set_electrons_elec_beta_num(num_beta)
print("OK")
print(f"{num_alpha} {num_beta}")
print("Basis\t\t...\t", end=' ')
@ -126,9 +132,7 @@ def write_ezfio(trexio_filename, filename):
try:
basis_type = trexio.read_basis_type(trexio_file)
if basis_type.lower() not in ["gaussian", "slater"]:
raise TypeError
if basis_type.lower() in ["gaussian", "slater"]:
shell_num = trexio.read_basis_shell_num(trexio_file)
prim_num = trexio.read_basis_prim_num(trexio_file)
ang_mom = trexio.read_basis_shell_ang_mom(trexio_file)
@ -139,6 +143,7 @@ def write_ezfio(trexio_filename, filename):
ao_shell = trexio.read_ao_shell(trexio_file)
ezfio.set_basis_basis("Read from TREXIO")
ezfio.set_ao_basis_ao_basis("Read from TREXIO")
ezfio.set_basis_shell_num(shell_num)
ezfio.set_basis_prim_num(prim_num)
ezfio.set_basis_shell_ang_mom(ang_mom)
@ -179,7 +184,61 @@ def write_ezfio(trexio_filename, filename):
shell_factor = trexio.read_basis_shell_factor(trexio_file)
prim_factor = trexio.read_basis_prim_factor(trexio_file)
print("OK")
elif basis_type.lower() == "numerical":
shell_num = trexio.read_basis_shell_num(trexio_file)
prim_num = shell_num
ang_mom = trexio.read_basis_shell_ang_mom(trexio_file)
nucl_index = trexio.read_basis_nucleus_index(trexio_file)
exponent = [1.]*prim_num
coefficient = [1.]*prim_num
shell_index = [i for i in range(shell_num)]
ao_shell = trexio.read_ao_shell(trexio_file)
ezfio.set_basis_basis("None")
ezfio.set_ao_basis_ao_basis("None")
ezfio.set_basis_shell_num(shell_num)
ezfio.set_basis_prim_num(prim_num)
ezfio.set_basis_shell_ang_mom(ang_mom)
ezfio.set_basis_basis_nucleus_index([ x+1 for x in nucl_index ])
ezfio.set_basis_prim_expo(exponent)
ezfio.set_basis_prim_coef(coefficient)
nucl_shell_num = []
prev = None
m = 0
for i in ao_shell:
if i != prev:
m += 1
if prev is None or nucl_index[i] != nucl_index[prev]:
nucl_shell_num.append(m)
m = 0
prev = i
assert (len(nucl_shell_num) == nucl_num)
shell_prim_num = []
prev = shell_index[0]
count = 0
for i in shell_index:
if i != prev:
shell_prim_num.append(count)
count = 0
count += 1
prev = i
shell_prim_num.append(count)
assert (len(shell_prim_num) == shell_num)
ezfio.set_basis_shell_prim_num(shell_prim_num)
ezfio.set_basis_shell_index([x+1 for x in shell_index])
ezfio.set_basis_nucleus_shell_num(nucl_shell_num)
shell_factor = trexio.read_basis_shell_factor(trexio_file)
prim_factor = [1.]*prim_num
else:
raise TypeError
print(basis_type)
except:
print("None")
ezfio.set_ao_basis_ao_cartesian(True)
@ -256,10 +315,12 @@ def write_ezfio(trexio_filename, filename):
# ezfio.set_ao_basis_ao_prim_num_max(prim_num_max)
ezfio.set_ao_basis_ao_coef(coef)
ezfio.set_ao_basis_ao_expo(expo)
ezfio.set_ao_basis_ao_basis("Read from TREXIO")
print("OK")
else:
print("None")
# _
# |\/| _ _ |_) _. _ o _
@ -279,6 +340,7 @@ def write_ezfio(trexio_filename, filename):
except:
label = "None"
ezfio.set_mo_basis_mo_label(label)
ezfio.set_determinants_mo_label(label)
try:
clss = trexio.read_mo_class(trexio_file)
@ -303,10 +365,10 @@ def write_ezfio(trexio_filename, filename):
for i in range(num_beta):
mo_occ[i] += 1.
ezfio.set_mo_basis_mo_occ(mo_occ)
except:
pass
print("OK")
except:
print("None")
print("Pseudos\t\t...\t", end=' ')
@ -386,10 +448,11 @@ def write_ezfio(trexio_filename, filename):
ezfio.set_pseudo_pseudo_n_kl(pseudo_n_kl)
ezfio.set_pseudo_pseudo_v_kl(pseudo_v_kl)
ezfio.set_pseudo_pseudo_dz_kl(pseudo_dz_kl)
print("OK")
else:
print("None")

12
src/ao_basis/EZFIO.cfg
View File

@ -67,3 +67,15 @@ doc: Use normalized primitive functions
interface: ezfio, provider
default: true
[ao_expoim_cosgtos]
type: double precision
doc: imag part for Exponents for each primitive of each cosGTOs |AO|
size: (ao_basis.ao_num,ao_basis.ao_prim_num_max)
interface: ezfio, provider
[use_cosgtos]
type: logical
doc: If true, use cosgtos for AO integrals
interface: ezfio
default: False

33
src/ao_basis/cosgtos.irp.f Normal file
View File

@ -0,0 +1,33 @@
BEGIN_PROVIDER [ logical, use_cosgtos ]
implicit none
BEGIN_DOC
! If true, use cosgtos for AO integrals
END_DOC
logical :: has
PROVIDE ezfio_filename
if (mpi_master) then
call ezfio_has_ao_basis_use_cosgtos(has)
if (has) then
! write(6,'(A)') '.. >>>>> [ IO READ: use_cosgtos ] <<<<< ..'
call ezfio_get_ao_basis_use_cosgtos(use_cosgtos)
else
use_cosgtos = .False.
endif
endif
IRP_IF MPI_DEBUG
print *, irp_here, mpi_rank
call MPI_BARRIER(MPI_COMM_WORLD, ierr)
IRP_ENDIF
IRP_IF MPI
include 'mpif.h'
integer :: ierr
call MPI_BCAST( use_cosgtos, 1, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr)
if (ierr /= MPI_SUCCESS) then
stop 'Unable to read use_cosgtos with MPI'
endif
IRP_ENDIF
! call write_time(6)
END_PROVIDER

1
src/ao_one_e_ints/NEED
View File

@ -1,3 +1,2 @@
ao_basis
pseudo
cosgtos_ao_int

0
src/cosgtos_ao_int/aos_cosgtos.irp.f → src/ao_one_e_ints/aos_cosgtos.irp.f
View File

0
src/cosgtos_ao_int/one_e_Coul_integrals.irp.f → src/ao_one_e_ints/one_e_Coul_integrals_cosgtos.irp.f
View File

0
src/cosgtos_ao_int/one_e_kin_integrals.irp.f → src/ao_one_e_ints/one_e_kin_integrals_cosgtos.irp.f
View File

21
src/ao_one_e_ints/pot_ao_ints.irp.f
View File

@ -455,10 +455,12 @@ recursive subroutine I_x1_pol_mult_one_e(a,c,R1x,R1xp,R2x,d,nd,n_pt_in)
do ix=0,nx
X(ix) *= dble(c)
enddo
call multiply_poly(X,nx,R2x,2,d,nd)
! call multiply_poly(X,nx,R2x,2,d,nd)
call multiply_poly_c2(X,nx,R2x,d,nd)
ny=0
call I_x2_pol_mult_one_e(c,R1x,R1xp,R2x,Y,ny,n_pt_in)
call multiply_poly(Y,ny,R1x,2,d,nd)
! call multiply_poly(Y,ny,R1x,2,d,nd)
call multiply_poly_c2(Y,ny,R1x,d,nd)
else
do ix=0,n_pt_in
X(ix) = 0.d0
@ -469,7 +471,8 @@ recursive subroutine I_x1_pol_mult_one_e(a,c,R1x,R1xp,R2x,d,nd,n_pt_in)
do ix=0,nx
X(ix) *= dble(a-1)
enddo
call multiply_poly(X,nx,R2x,2,d,nd)
! call multiply_poly(X,nx,R2x,2,d,nd)
call multiply_poly_c2(X,nx,R2x,d,nd)
nx = nd
do ix=0,n_pt_in
@ -479,10 +482,12 @@ recursive subroutine I_x1_pol_mult_one_e(a,c,R1x,R1xp,R2x,d,nd,n_pt_in)
do ix=0,nx
X(ix) *= dble(c)
enddo
call multiply_poly(X,nx,R2x,2,d,nd)
! call multiply_poly(X,nx,R2x,2,d,nd)
call multiply_poly_c2(X,nx,R2x,d,nd)
ny=0
call I_x1_pol_mult_one_e(a-1,c,R1x,R1xp,R2x,Y,ny,n_pt_in)
call multiply_poly(Y,ny,R1x,2,d,nd)
! call multiply_poly(Y,ny,R1x,2,d,nd)
call multiply_poly_c2(Y,ny,R1x,d,nd)
endif
end
@ -519,7 +524,8 @@ recursive subroutine I_x2_pol_mult_one_e(c,R1x,R1xp,R2x,d,nd,dim)
do ix=0,nx
X(ix) *= dble(c-1)
enddo
call multiply_poly(X,nx,R2x,2,d,nd)
! call multiply_poly(X,nx,R2x,2,d,nd)
call multiply_poly_c2(X,nx,R2x,d,nd)
ny = 0
do ix=0,dim
Y(ix) = 0.d0
@ -527,7 +533,8 @@ recursive subroutine I_x2_pol_mult_one_e(c,R1x,R1xp,R2x,d,nd,dim)
call I_x1_pol_mult_one_e(0,c-1,R1x,R1xp,R2x,Y,ny,dim)
if(ny.ge.0)then
call multiply_poly(Y,ny,R1xp,2,d,nd)
! call multiply_poly(Y,ny,R1xp,2,d,nd)
call multiply_poly_c2(Y,ny,R1xp,d,nd)
endif
endif
end

13
src/ao_two_e_ints/EZFIO.cfg
View File

@ -4,6 +4,19 @@ doc: Read/Write |AO| integrals from/to disk [ Write | Read | None ]
interface: ezfio,provider,ocaml
default: None
[ao_integrals_threshold]
type: Threshold
doc: If | (pq|rs) | < `ao_integrals_threshold` then (pq|rs) is zero
interface: ezfio,provider,ocaml
default: 1.e-15
ezfio_name: threshold_ao
[ao_cholesky_threshold]
type: Threshold
doc: If | (ii|jj) | < `ao_cholesky_threshold` then (ii|jj) is zero
interface: ezfio,provider,ocaml
default: 1.e-12
[do_direct_integrals]
type: logical
doc: Compute integrals on the fly (very slow, only for debugging)

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

@ -4,29 +4,7 @@ BEGIN_PROVIDER [ integer, cholesky_ao_num_guess ]
! Number of Cholesky vectors in AO basis
END_DOC
integer :: i,j,k,l
double precision :: xnorm0, x, integral
double precision, external :: ao_two_e_integral
cholesky_ao_num_guess = 0
xnorm0 = 0.d0
x = 0.d0
do j=1,ao_num
do i=1,ao_num
integral = ao_two_e_integral(i,i,j,j)
if (integral > ao_integrals_threshold) then
cholesky_ao_num_guess += 1
else
x += integral
endif
enddo
enddo
print *, 'Cholesky decomposition of AO integrals'
print *, '--------------------------------------'
print *, ''
print *, 'Estimated Error: ', x
print *, 'Guess size: ', cholesky_ao_num_guess, '(', 100.d0*dble(cholesky_ao_num_guess)/dble(ao_num*ao_num), ' %)'
cholesky_ao_num_guess = ao_num*ao_num / 2
END_PROVIDER
BEGIN_PROVIDER [ integer, cholesky_ao_num ]
@ -39,7 +17,7 @@ END_PROVIDER
END_DOC
type(c_ptr) :: ptr
integer :: fd, i,j,k,l, rank
integer :: fd, i,j,k,l,m,rank
double precision, pointer :: ao_integrals(:,:,:,:)
double precision, external :: ao_two_e_integral
@ -49,28 +27,90 @@ END_PROVIDER
8, fd, .False., ptr)
call c_f_pointer(ptr, ao_integrals, (/ao_num, ao_num, ao_num, ao_num/))
double precision :: integral
print*, 'Providing the AO integrals (Cholesky)'
call wall_time(wall_1)
call cpu_time(cpu_1)
ao_integrals = 0.d0
double precision :: integral, cpu_1, cpu_2, wall_1, wall_2
logical, external :: ao_two_e_integral_zero
!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i,j,k,l, integral) SCHEDULE(dynamic)
do l=1,ao_num
double precision, external :: get_ao_two_e_integral
if (read_ao_two_e_integrals) then
PROVIDE ao_two_e_integrals_in_map
!$OMP PARALLEL DEFAULT(SHARED) PRIVATE(i,j,k,l, integral, wall_2)
do m=0,9
do l=1+m,ao_num,10
!$OMP DO SCHEDULE(dynamic)
do j=1,l
do k=1,ao_num
do i=1,k
do i=1,min(k,j)
if (ao_two_e_integral_zero(i,j,k,l)) cycle
integral = get_ao_two_e_integral(i,j,k,l, ao_integrals_map)
ao_integrals(i,k,j,l) = integral
ao_integrals(k,i,j,l) = integral
ao_integrals(i,k,l,j) = integral
ao_integrals(k,i,l,j) = integral
ao_integrals(j,l,i,k) = integral
ao_integrals(j,l,k,i) = integral
ao_integrals(l,j,i,k) = integral
ao_integrals(l,j,k,i) = integral
enddo
enddo
enddo
!$OMP END DO NOWAIT
enddo
!$OMP MASTER
call wall_time(wall_2)
print '(I10,'' % in'', 4X, F10.2, '' s.'')', (m+1) * 10, wall_2-wall_1
!$OMP END MASTER
enddo
!$OMP END PARALLEL
else
!$OMP PARALLEL DEFAULT(SHARED) PRIVATE(i,j,k,l, integral, wall_2)
do m=0,9
do l=1+m,ao_num,10
!$OMP DO SCHEDULE(dynamic)
do j=1,l
do k=1,ao_num
do i=1,min(k,j)
if (ao_two_e_integral_zero(i,j,k,l)) cycle
integral = ao_two_e_integral(i,k,j,l)
ao_integrals(i,k,j,l) = integral
ao_integrals(k,i,j,l) = integral
ao_integrals(i,k,l,j) = integral
ao_integrals(k,i,l,j) = integral
ao_integrals(j,l,i,k) = integral
ao_integrals(j,l,k,i) = integral
ao_integrals(l,j,i,k) = integral
ao_integrals(l,j,k,i) = integral
enddo
enddo
enddo
!$OMP END DO NOWAIT
enddo
!$OMP END PARALLEL DO
!$OMP MASTER
call wall_time(wall_2)
print '(I10,'' % in'', 4X, F10.2, '' s.'')', (m+1) * 10, wall_2-wall_1
!$OMP END MASTER
enddo
!$OMP END PARALLEL
call wall_time(wall_2)
call cpu_time(cpu_2)
print*, 'AO integrals provided:'
print*, ' cpu time :',cpu_2 - cpu_1, 's'
print*, ' wall time :',wall_2 - wall_1, 's ( x ', (cpu_2-cpu_1)/(wall_2-wall_1+tiny(1.d0)), ' )'
endif
! Call Lapack
cholesky_ao_num = cholesky_ao_num_guess
call pivoted_cholesky(ao_integrals, cholesky_ao_num, ao_integrals_threshold, ao_num*ao_num, cholesky_ao)
call pivoted_cholesky(ao_integrals, cholesky_ao_num, ao_cholesky_threshold, ao_num*ao_num, cholesky_ao)
print *, 'Rank: ', cholesky_ao_num, '(', 100.d0*dble(cholesky_ao_num)/dble(ao_num*ao_num), ' %)'
! Remove mmap

0
src/cosgtos_ao_int/gauss_legendre.irp.f → src/ao_two_e_ints/gauss_legendre.irp.f
View File

0
src/cosgtos_ao_int/two_e_Coul_integrals.irp.f → src/ao_two_e_ints/two_e_Coul_integrals_cosgtos.irp.f
View File

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

@ -590,8 +590,20 @@ double precision function general_primitive_integral(dim, &
d_poly(i)=0.d0
enddo
!DIR$ FORCEINLINE
call multiply_poly(Ix_pol,n_Ix,Iy_pol,n_Iy,d_poly,n_pt_tmp)
! call multiply_poly(Ix_pol,n_Ix,Iy_pol,n_Iy,d_poly,n_pt_tmp)
integer :: ib, ic
if (ior(n_Ix,n_Iy) >= 0) then
do ib=0,n_Ix
do ic = 0,n_Iy
d_poly(ib+ic) = d_poly(ib+ic) + Iy_pol(ic) * Ix_pol(ib)
enddo
enddo
do n_pt_tmp = n_Ix+n_Iy, 0, -1
if (d_poly(n_pt_tmp) /= 0.d0) exit
enddo
endif
if (n_pt_tmp == -1) then
return
endif
@ -600,8 +612,21 @@ double precision function general_primitive_integral(dim, &
d1(i)=0.d0
enddo
!DIR$ FORCEINLINE
call multiply_poly(d_poly ,n_pt_tmp ,Iz_pol,n_Iz,d1,n_pt_out)
! call multiply_poly(d_poly ,n_pt_tmp ,Iz_pol,n_Iz,d1,n_pt_out)
if (ior(n_pt_tmp,n_Iz) >= 0) then
! Bottleneck here
do ib=0,n_pt_tmp
do ic = 0,n_Iz
d1(ib+ic) = d1(ib+ic) + Iz_pol(ic) * d_poly(ib)
enddo
enddo
do n_pt_out = n_pt_tmp+n_Iz, 0, -1
if (d1(n_pt_out) /= 0.d0) exit
enddo
endif
double precision :: rint_sum
accu = accu + rint_sum(n_pt_out,const,d1)
@ -948,8 +973,9 @@ recursive subroutine I_x1_pol_mult_recurs(a,c,B_10,B_01,B_00,C_00,D_00,d,nd,n_pt
X(ix) *= dble(a-1)
enddo
!DIR$ FORCEINLINE
call multiply_poly(X,nx,B_10,2,d,nd)
! !DIR$ FORCEINLINE
! call multiply_poly(X,nx,B_10,2,d,nd)
call multiply_poly_c2(X,nx,B_10,d,nd)
nx = nd
!DIR$ LOOP COUNT(8)
@ -970,8 +996,9 @@ recursive subroutine I_x1_pol_mult_recurs(a,c,B_10,B_01,B_00,C_00,D_00,d,nd,n_pt
X(ix) *= c
enddo
endif
!DIR$ FORCEINLINE
call multiply_poly(X,nx,B_00,2,d,nd)
! !DIR$ FORCEINLINE
! call multiply_poly(X,nx,B_00,2,d,nd)
call multiply_poly_c2(X,nx,B_00,d,nd)
endif
ny=0
@ -988,9 +1015,9 @@ recursive subroutine I_x1_pol_mult_recurs(a,c,B_10,B_01,B_00,C_00,D_00,d,nd,n_pt
call I_x1_pol_mult_recurs(a-1,c,B_10,B_01,B_00,C_00,D_00,Y,ny,n_pt_in)
endif
!DIR$ FORCEINLINE
call multiply_poly(Y,ny,C_00,2,d,nd)
! !DIR$ FORCEINLINE
! call multiply_poly(Y,ny,C_00,2,d,nd)
call multiply_poly_c2(Y,ny,C_00,d,nd)
end
recursive subroutine I_x1_pol_mult_a1(c,B_10,B_01,B_00,C_00,D_00,d,nd,n_pt_in)
@ -1028,8 +1055,9 @@ recursive subroutine I_x1_pol_mult_a1(c,B_10,B_01,B_00,C_00,D_00,d,nd,n_pt_in)
enddo
endif
!DIR$ FORCEINLINE
call multiply_poly(X,nx,B_00,2,d,nd)
! !DIR$ FORCEINLINE
! call multiply_poly(X,nx,B_00,2,d,nd)
call multiply_poly_c2(X,nx,B_00,d,nd)
ny=0
@ -1039,8 +1067,9 @@ recursive subroutine I_x1_pol_mult_a1(c,B_10,B_01,B_00,C_00,D_00,d,nd,n_pt_in)
enddo
call I_x2_pol_mult(c,B_10,B_01,B_00,C_00,D_00,Y,ny,n_pt_in)
!DIR$ FORCEINLINE
call multiply_poly(Y,ny,C_00,2,d,nd)
! !DIR$ FORCEINLINE
! call multiply_poly(Y,ny,C_00,2,d,nd)
call multiply_poly_c2(Y,ny,C_00,d,nd)
end
@ -1067,8 +1096,9 @@ recursive subroutine I_x1_pol_mult_a2(c,B_10,B_01,B_00,C_00,D_00,d,nd,n_pt_in)
nx = 0
call I_x2_pol_mult(c,B_10,B_01,B_00,C_00,D_00,X,nx,n_pt_in)
!DIR$ FORCEINLINE
call multiply_poly(X,nx,B_10,2,d,nd)
! !DIR$ FORCEINLINE
! call multiply_poly(X,nx,B_10,2,d,nd)
call multiply_poly_c2(X,nx,B_10,d,nd)
nx = nd
!DIR$ LOOP COUNT(8)
@ -1086,8 +1116,9 @@ recursive subroutine I_x1_pol_mult_a2(c,B_10,B_01,B_00,C_00,D_00,d,nd,n_pt_in)
enddo
endif
!DIR$ FORCEINLINE
call multiply_poly(X,nx,B_00,2,d,nd)
! !DIR$ FORCEINLINE
! call multiply_poly(X,nx,B_00,2,d,nd)
call multiply_poly_c2(X,nx,B_00,d,nd)
ny=0
!DIR$ LOOP COUNT(8)
@ -1097,9 +1128,9 @@ recursive subroutine I_x1_pol_mult_a2(c,B_10,B_01,B_00,C_00,D_00,d,nd,n_pt_in)
!DIR$ FORCEINLINE
call I_x1_pol_mult_a1(c,B_10,B_01,B_00,C_00,D_00,Y,ny,n_pt_in)
!DIR$ FORCEINLINE
call multiply_poly(Y,ny,C_00,2,d,nd)
! !DIR$ FORCEINLINE
! call multiply_poly(Y,ny,C_00,2,d,nd)
call multiply_poly_c2(Y,ny,C_00,d,nd)
end
recursive subroutine I_x2_pol_mult(c,B_10,B_01,B_00,C_00,D_00,d,nd,dim)
@ -1146,8 +1177,10 @@ recursive subroutine I_x2_pol_mult(c,B_10,B_01,B_00,C_00,D_00,d,nd,dim)
Y(1) = D_00(1)
Y(2) = D_00(2)
!DIR$ FORCEINLINE
call multiply_poly(Y,ny,D_00,2,d,nd)
! !DIR$ FORCEINLINE
! call multiply_poly(Y,ny,D_00,2,d,nd)
call multiply_poly_c2(Y,ny,D_00,d,nd)
return
case default
@ -1164,8 +1197,9 @@ recursive subroutine I_x2_pol_mult(c,B_10,B_01,B_00,C_00,D_00,d,nd,dim)
X(ix) *= dble(c-1)
enddo
!DIR$ FORCEINLINE
call multiply_poly(X,nx,B_01,2,d,nd)
! !DIR$ FORCEINLINE
! call multiply_poly(X,nx,B_01,2,d,nd)
call multiply_poly_c2(X,nx,B_01,d,nd)
ny = 0
!DIR$ LOOP COUNT(6)
@ -1174,8 +1208,9 @@ recursive subroutine I_x2_pol_mult(c,B_10,B_01,B_00,C_00,D_00,d,nd,dim)
enddo
call I_x2_pol_mult(c-1,B_10,B_01,B_00,C_00,D_00,Y,ny,dim)
!DIR$ FORCEINLINE
call multiply_poly(Y,ny,D_00,2,d,nd)
! !DIR$ FORCEINLINE
! call multiply_poly(Y,ny,D_00,2,d,nd)
call multiply_poly_c2(Y,ny,D_00,d,nd)
end select
end
@ -1233,3 +1268,34 @@ subroutine compute_ao_integrals_jl(j,l,n_integrals,buffer_i,buffer_value)
enddo
end
subroutine multiply_poly_local(b,nb,c,nc,d,nd)
implicit none
BEGIN_DOC
! Multiply two polynomials
! D(t) += B(t)*C(t)
END_DOC
integer, intent(in) :: nb, nc
integer, intent(out) :: nd
double precision, intent(in) :: b(0:nb), c(0:nc)
double precision, intent(inout) :: d(0:nb+nc)
integer :: ndtmp
integer :: ib, ic, id, k
if(ior(nc,nb) < 0) return !False if nc>=0 and nb>=0
do ib=0,nb
do ic = 0,nc
d(ib+ic) = d(ib+ic) + c(ic) * b(ib)
enddo
enddo
do nd = nb+nc,0,-1
if (d(nd) /= 0.d0) exit
enddo
end

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

@ -7,7 +7,13 @@ program bi_ort_ints
my_n_pt_r_grid = 10
my_n_pt_a_grid = 14
touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
call test_3e
! call test_3e
call test_5idx
! call test_5idx2
end
subroutine test_5idx2
PROVIDE three_e_5_idx_cycle_2_bi_ort
end
subroutine test_3e
@ -16,6 +22,7 @@ subroutine test_3e
double precision :: accu, contrib,new,ref
i = 1
k = 1
n = 0
accu = 0.d0
do i = 1, mo_num
do k = 1, mo_num
@ -31,6 +38,7 @@ subroutine test_3e
print*,'pb !!'
print*,i,k,j,l,m,n
print*,ref,new,contrib
stop
endif
enddo
enddo
@ -42,3 +50,93 @@ subroutine test_3e
end
subroutine test_5idx
implicit none
integer :: i,k,j,l,m,n,ipoint
double precision :: accu, contrib,new,ref
i = 1
k = 1
n = 0
accu = 0.d0
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
new = three_e_5_idx_direct_bi_ort(m,l,j,k,i)
ref = three_e_5_idx_direct_bi_ort_old(m,l,j,k,i)
contrib = dabs(new - ref)
accu += contrib
if(contrib .gt. 1.d-10)then
print*,'direct'
print*,i,k,j,l,m
print*,ref,new,contrib
stop
endif
new = three_e_5_idx_exch12_bi_ort(m,l,j,k,i)
ref = three_e_5_idx_exch12_bi_ort_old(m,l,j,k,i)
contrib = dabs(new - ref)
accu += contrib
if(contrib .gt. 1.d-10)then
print*,'exch12'
print*,i,k,j,l,m
print*,ref,new,contrib
stop
endif
!
new = three_e_5_idx_cycle_1_bi_ort(m,l,j,k,i)
ref = three_e_5_idx_cycle_1_bi_ort_old(m,l,j,k,i)
contrib = dabs(new - ref)
accu += contrib
if(contrib .gt. 1.d-10)then
print*,'cycle1'
print*,i,k,j,l,m
print*,ref,new,contrib
stop
endif
new = three_e_5_idx_cycle_2_bi_ort(m,l,j,k,i)
ref = three_e_5_idx_cycle_2_bi_ort_old(m,l,j,k,i)
contrib = dabs(new - ref)
accu += contrib
if(contrib .gt. 1.d-10)then
print*,'cycle2'
print*,i,k,j,l,m
print*,ref,new,contrib
stop
endif
new = three_e_5_idx_exch23_bi_ort(m,l,j,k,i)
ref = three_e_5_idx_exch23_bi_ort_old(m,l,j,k,i)
contrib = dabs(new - ref)
accu += contrib
if(contrib .gt. 1.d-10)then
print*,'exch23'
print*,i,k,j,l,m
print*,ref,new,contrib
stop
endif
new = three_e_5_idx_exch13_bi_ort(m,l,j,k,i)
ref = three_e_5_idx_exch13_bi_ort_old(m,l,j,k,i)
contrib = dabs(new - ref)
accu += contrib
if(contrib .gt. 1.d-10)then
print*,'exch13'
print*,i,k,j,l,m
print*,ref,new,contrib
stop
endif
enddo
enddo
enddo
enddo
enddo
print*,'accu = ',accu/dble(mo_num)**5
end

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

@ -1,7 +1,11 @@
! ---
BEGIN_PROVIDER [ double precision, three_e_5_idx_direct_bi_ort, (mo_num, mo_num, mo_num, mo_num, mo_num)]
BEGIN_PROVIDER [ double precision, three_e_5_idx_direct_bi_ort , (mo_num, mo_num, mo_num, mo_num, mo_num)]
&BEGIN_PROVIDER [ double precision, three_e_5_idx_exch12_bi_ort , (mo_num, mo_num, mo_num, mo_num, mo_num)]
&BEGIN_PROVIDER [ double precision, three_e_5_idx_exch23_bi_ort , (mo_num, mo_num, mo_num, mo_num, mo_num)]
&BEGIN_PROVIDER [ double precision, three_e_5_idx_exch13_bi_ort , (mo_num, mo_num, mo_num, mo_num, mo_num)]
&BEGIN_PROVIDER [ double precision, three_e_5_idx_cycle_1_bi_ort, (mo_num, mo_num, mo_num, mo_num, mo_num)]
&BEGIN_PROVIDER [ double precision, three_e_5_idx_cycle_2_bi_ort, (mo_num, mo_num, mo_num, mo_num, mo_num)]
BEGIN_DOC
!
@ -14,289 +18,221 @@ BEGIN_PROVIDER [ double precision, three_e_5_idx_direct_bi_ort, (mo_num, mo_num,
implicit none
integer :: i, j, k, m, l
double precision :: integral, wall1, wall0
three_e_5_idx_direct_bi_ort = 0.d0
print *, ' Providing the three_e_5_idx_direct_bi_ort ...'
call wall_time(wall0)
double precision :: wall1, wall0
integer :: ipoint
double precision, allocatable :: grad_mli(:,:,:), orb_mat(:,:,:)
double precision, allocatable :: lk_grad_mi(:,:,:,:), rk_grad_im(:,:,:,:)
double precision, allocatable :: lm_grad_ik(:,:,:,:), rm_grad_ik(:,:,:,:)
double precision, allocatable :: tmp_mat(:,:,:,:)
allocate(tmp_mat(mo_num,mo_num,mo_num,mo_num))
provide mos_r_in_r_array_transp mos_l_in_r_array_transp
PROVIDE mo_l_coef mo_r_coef int2_grad1_u12_bimo_t
print *, ' Providing the three_e_5_idx_bi_ort ...'
call wall_time(wall0)
do m = 1, mo_num
allocate(grad_mli(n_points_final_grid,mo_num,mo_num))
allocate(orb_mat(n_points_final_grid,mo_num,mo_num))
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,j,k,m,l,integral) &
!$OMP SHARED (mo_num,three_e_5_idx_direct_bi_ort)
!$OMP DO SCHEDULE (dynamic) COLLAPSE(2)
!$OMP PRIVATE (i,l,ipoint) &
!$OMP SHARED (m,mo_num,n_points_final_grid, &
!$OMP mos_l_in_r_array_transp, mos_r_in_r_array_transp, &
!$OMP int2_grad1_u12_bimo_t, final_weight_at_r_vector, &
!$OMP grad_mli, orb_mat)
!$OMP DO COLLAPSE(2)
do i=1,mo_num
do l=1,mo_num
do ipoint=1, n_points_final_grid
grad_mli(ipoint,l,i) = final_weight_at_r_vector(ipoint) * ( &
int2_grad1_u12_bimo_t(ipoint,1,m,m) * int2_grad1_u12_bimo_t(ipoint,1,l,i) + &
int2_grad1_u12_bimo_t(ipoint,2,m,m) * int2_grad1_u12_bimo_t(ipoint,2,l,i) + &
int2_grad1_u12_bimo_t(ipoint,3,m,m) * int2_grad1_u12_bimo_t(ipoint,3,l,i) )
orb_mat(ipoint,l,i) = mos_l_in_r_array_transp(ipoint,l) * mos_r_in_r_array_transp(ipoint,i)
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
call dgemm('T','N', mo_num*mo_num, mo_num*mo_num, n_points_final_grid, 1.d0, &
orb_mat, n_points_final_grid, &
grad_mli, n_points_final_grid, 0.d0, &
tmp_mat, mo_num*mo_num)
!$OMP PARALLEL DO PRIVATE(i,j,k,l)
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
call give_integrals_3_body_bi_ort(m, l, k, m, j, i, integral)
three_e_5_idx_direct_bi_ort(m,l,j,k,i) = -1.d0 * integral
three_e_5_idx_direct_bi_ort(m,l,j,k,i) = - tmp_mat(l,j,k,i) - tmp_mat(k,i,l,j)
three_e_5_idx_exch12_bi_ort(m,l,j,k,i) = - tmp_mat(l,i,k,j) - tmp_mat(k,j,l,i)
enddo
enddo
enddo
enddo
!$OMP END PARALLEL DO
deallocate(orb_mat,grad_mli)
allocate(lm_grad_ik(n_points_final_grid,3,mo_num,mo_num))
allocate(rm_grad_ik(n_points_final_grid,3,mo_num,mo_num))
allocate(rk_grad_im(n_points_final_grid,3,mo_num,mo_num))
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,l,ipoint) &
!$OMP SHARED (m,mo_num,n_points_final_grid, &
!$OMP mos_l_in_r_array_transp, mos_r_in_r_array_transp, &
!$OMP int2_grad1_u12_bimo_t, final_weight_at_r_vector, &
!$OMP rm_grad_ik, lm_grad_ik, rk_grad_im, lk_grad_mi)
!$OMP DO COLLAPSE(2)
do i=1,mo_num
do l=1,mo_num
do ipoint=1, n_points_final_grid
lm_grad_ik(ipoint,1,l,i) = mos_l_in_r_array_transp(ipoint,m) * int2_grad1_u12_bimo_t(ipoint,1,l,i) * final_weight_at_r_vector(ipoint)
lm_grad_ik(ipoint,2,l,i) = mos_l_in_r_array_transp(ipoint,m) * int2_grad1_u12_bimo_t(ipoint,2,l,i) * final_weight_at_r_vector(ipoint)
lm_grad_ik(ipoint,3,l,i) = mos_l_in_r_array_transp(ipoint,m) * int2_grad1_u12_bimo_t(ipoint,3,l,i) * final_weight_at_r_vector(ipoint)
rm_grad_ik(ipoint,1,l,i) = mos_r_in_r_array_transp(ipoint,m) * int2_grad1_u12_bimo_t(ipoint,1,l,i)
rm_grad_ik(ipoint,2,l,i) = mos_r_in_r_array_transp(ipoint,m) * int2_grad1_u12_bimo_t(ipoint,2,l,i)
rm_grad_ik(ipoint,3,l,i) = mos_r_in_r_array_transp(ipoint,m) * int2_grad1_u12_bimo_t(ipoint,3,l,i)
rk_grad_im(ipoint,1,l,i) = mos_r_in_r_array_transp(ipoint,l) * int2_grad1_u12_bimo_t(ipoint,1,i,m)
rk_grad_im(ipoint,2,l,i) = mos_r_in_r_array_transp(ipoint,l) * int2_grad1_u12_bimo_t(ipoint,2,i,m)
rk_grad_im(ipoint,3,l,i) = mos_r_in_r_array_transp(ipoint,l) * int2_grad1_u12_bimo_t(ipoint,3,i,m)
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
call dgemm('T','N', mo_num*mo_num, mo_num*mo_num, 3*n_points_final_grid, 1.d0, &
lm_grad_ik, 3*n_points_final_grid, &
rm_grad_ik, 3*n_points_final_grid, 0.d0, &
tmp_mat, mo_num*mo_num)
!$OMP PARALLEL DO PRIVATE(i,j,k,l)
do i = 1, mo_num
do k = 1, mo_num
do j = 1, mo_num
do l = 1, mo_num
three_e_5_idx_direct_bi_ort(m,l,j,k,i) = three_e_5_idx_direct_bi_ort(m,l,j,k,i) - tmp_mat(l,j,k,i)
three_e_5_idx_exch12_bi_ort(m,l,j,k,i) = three_e_5_idx_exch12_bi_ort(m,l,j,k,i) - tmp_mat(l,i,k,j)
enddo
enddo
enddo
enddo
!$OMP END PARALLEL DO
call dgemm('T','N', mo_num*mo_num, mo_num*mo_num, 3*n_points_final_grid, 1.d0, &
lm_grad_ik, 3*n_points_final_grid, &
rk_grad_im, 3*n_points_final_grid, 0.d0, &
tmp_mat, mo_num*mo_num)
!$OMP PARALLEL DO PRIVATE(i,j,k,l)
do i = 1, mo_num
do k = 1, mo_num
do j = 1, mo_num
do l = 1, mo_num
three_e_5_idx_cycle_1_bi_ort(m,l,j,k,i) = - tmp_mat(l,i,j,k)
three_e_5_idx_cycle_2_bi_ort(m,l,j,k,i) = - tmp_mat(k,j,i,l)
three_e_5_idx_exch23_bi_ort (m,l,j,k,i) = - tmp_mat(k,i,j,l)
three_e_5_idx_exch13_bi_ort (m,l,j,k,i) = - tmp_mat(l,j,i,k)
enddo
enddo
enddo
enddo
!$OMP END PARALLEL DO
deallocate(lm_grad_ik)
allocate(lk_grad_mi(n_points_final_grid,3,mo_num,mo_num))
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,l,ipoint) &
!$OMP SHARED (m,mo_num,n_points_final_grid, &
!$OMP mos_l_in_r_array_transp, mos_r_in_r_array_transp, &
!$OMP int2_grad1_u12_bimo_t, final_weight_at_r_vector, &
!$OMP lk_grad_mi)
!$OMP DO COLLAPSE(2)
do i=1,mo_num
do l=1,mo_num
do ipoint=1, n_points_final_grid
lk_grad_mi(ipoint,1,l,i) = mos_l_in_r_array_transp(ipoint,l) * int2_grad1_u12_bimo_t(ipoint,1,m,i) * final_weight_at_r_vector(ipoint)
lk_grad_mi(ipoint,2,l,i) = mos_l_in_r_array_transp(ipoint,l) * int2_grad1_u12_bimo_t(ipoint,2,m,i) * final_weight_at_r_vector(ipoint)
lk_grad_mi(ipoint,3,l,i) = mos_l_in_r_array_transp(ipoint,l) * int2_grad1_u12_bimo_t(ipoint,3,m,i) * final_weight_at_r_vector(ipoint)
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
call wall_time(wall1)
print *, ' wall time for three_e_5_idx_direct_bi_ort', wall1 - wall0
call print_memory_usage()
call dgemm('T','N', mo_num*mo_num, mo_num*mo_num, 3*n_points_final_grid, 1.d0, &
lk_grad_mi, 3*n_points_final_grid, &
rm_grad_ik, 3*n_points_final_grid, 0.d0, &
tmp_mat, mo_num*mo_num)
END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, three_e_5_idx_cycle_1_bi_ort, (mo_num, mo_num, mo_num, mo_num, mo_num)]
BEGIN_DOC
!
! matrix element of the -L three-body operator FOR THE FIRST CYCLIC PERMUTATION TERMS OF DOUBLE EXCITATIONS AND BI ORTHO MOs
!
! three_e_5_idx_cycle_1_bi_ort(m,l,j,k,i) = <mlk|-L|jim> ::: 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_5_idx_cycle_1_bi_ort = 0.d0
print *, ' Providing the three_e_5_idx_cycle_1_bi_ort ...'
call wall_time(wall0)
provide mos_r_in_r_array_transp mos_l_in_r_array_transp
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,j,k,m,l,integral) &
!$OMP SHARED (mo_num,three_e_5_idx_cycle_1_bi_ort)
!$OMP DO SCHEDULE (dynamic) COLLAPSE(2)
!$OMP PARALLEL DO PRIVATE(i,j,k,l)
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
call give_integrals_3_body_bi_ort(m, l, k, j, i, m, integral)
three_e_5_idx_cycle_1_bi_ort(m,l,j,k,i) = -1.d0 * integral
three_e_5_idx_cycle_1_bi_ort(m,l,j,k,i) = three_e_5_idx_cycle_1_bi_ort(m,l,j,k,i) - tmp_mat(k,j,l,i)
three_e_5_idx_cycle_2_bi_ort(m,l,j,k,i) = three_e_5_idx_cycle_2_bi_ort(m,l,j,k,i) - tmp_mat(l,i,k,j)
three_e_5_idx_exch23_bi_ort (m,l,j,k,i) = three_e_5_idx_exch23_bi_ort (m,l,j,k,i) - tmp_mat(l,j,k,i)
three_e_5_idx_exch13_bi_ort (m,l,j,k,i) = three_e_5_idx_exch13_bi_ort (m,l,j,k,i) - tmp_mat(k,i,l,j)
enddo
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
!$OMP END PARALLEL DO
call wall_time(wall1)
print *, ' wall time for three_e_5_idx_cycle_1_bi_ort', wall1 - wall0
call print_memory_usage()
END_PROVIDER
call dgemm('T','N', mo_num*mo_num, mo_num*mo_num, 3*n_points_final_grid, 1.d0, &
lk_grad_mi, 3*n_points_final_grid, &
rk_grad_im, 3*n_points_final_grid, 0.d0, &
tmp_mat, mo_num*mo_num)
! ---
BEGIN_PROVIDER [ double precision, three_e_5_idx_cycle_2_bi_ort, (mo_num, mo_num, mo_num, mo_num, mo_num)]
BEGIN_DOC
!
! matrix element of the -L three-body operator FOR THE FIRST CYCLIC PERMUTATION TERMS OF DOUBLE EXCITATIONS AND BI ORTHO MOs
!
! three_e_5_idx_cycle_2_bi_ort(m,l,j,k,i) = <mlk|-L|imj> ::: 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_5_idx_cycle_2_bi_ort = 0.d0
print *, ' Providing the three_e_5_idx_cycle_2_bi_ort ...'
call wall_time(wall0)
provide mos_r_in_r_array_transp mos_l_in_r_array_transp
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,j,k,m,l,integral) &
!$OMP SHARED (mo_num,three_e_5_idx_cycle_2_bi_ort)
!$OMP DO SCHEDULE (dynamic) COLLAPSE(2)
do i = 1, mo_num
do k = 1, mo_num
do j = 1, mo_num
do m = 1, mo_num
do l = 1, mo_num
call give_integrals_3_body_bi_ort(m, l, k, i, m, j, integral)
three_e_5_idx_cycle_2_bi_ort(m,l,j,k,i) = -1.d0 * integral
enddo
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
call wall_time(wall1)
print *, ' wall time for three_e_5_idx_cycle_2_bi_ort', wall1 - wall0
call print_memory_usage()
END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, three_e_5_idx_exch23_bi_ort, (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_5_idx_exch23_bi_ort(m,l,j,k,i) = <mlk|-L|jmi> ::: 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_5_idx_exch23_bi_ort = 0.d0
print *, ' Providing the three_e_5_idx_exch23_bi_ort ...'
call wall_time(wall0)
provide mos_r_in_r_array_transp mos_l_in_r_array_transp
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,j,k,m,l,integral) &
!$OMP SHARED (mo_num,three_e_5_idx_exch23_bi_ort)
!$OMP DO SCHEDULE (dynamic) COLLAPSE(2)
!$OMP PARALLEL DO PRIVATE(i,j,k,l)
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
call give_integrals_3_body_bi_ort(m, l, k, j, m, i, integral)
three_e_5_idx_exch23_bi_ort(m,l,j,k,i) = -1.d0 * integral
three_e_5_idx_cycle_1_bi_ort(m,l,j,k,i) = three_e_5_idx_cycle_1_bi_ort(m,l,j,k,i) - tmp_mat(l,j,i,k)
three_e_5_idx_cycle_2_bi_ort(m,l,j,k,i) = three_e_5_idx_cycle_2_bi_ort(m,l,j,k,i) - tmp_mat(k,i,j,l)
three_e_5_idx_exch23_bi_ort (m,l,j,k,i) = three_e_5_idx_exch23_bi_ort (m,l,j,k,i) - tmp_mat(k,j,i,l)
three_e_5_idx_exch13_bi_ort (m,l,j,k,i) = three_e_5_idx_exch13_bi_ort (m,l,j,k,i) - tmp_mat(l,i,j,k)
enddo
enddo
enddo
enddo
!$OMP END PARALLEL DO
deallocate(lk_grad_mi)
deallocate(rm_grad_ik)
deallocate(rk_grad_im)
enddo
!$OMP END DO
!$OMP END PARALLEL
call wall_time(wall1)
print *, ' wall time for three_e_5_idx_exch23_bi_ort', wall1 - wall0
print *, ' wall time for three_e_5_idx_bi_ort', wall1 - wall0
call print_memory_usage()
END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, three_e_5_idx_exch13_bi_ort, (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_5_idx_exch13_bi_ort(m,l,j,k,i) = <mlk|-L|ijm> ::: 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_5_idx_exch13_bi_ort = 0.d0
print *, ' Providing the three_e_5_idx_exch13_bi_ort ...'
call wall_time(wall0)
provide mos_r_in_r_array_transp mos_l_in_r_array_transp
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,j,k,m,l,integral) &
!$OMP SHARED (mo_num,three_e_5_idx_exch13_bi_ort)
!$OMP DO SCHEDULE (dynamic) COLLAPSE(2)
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
call give_integrals_3_body_bi_ort(m, l, k, i, j, m, integral)
three_e_5_idx_exch13_bi_ort(m,l,j,k,i) = -1.d0 * integral
enddo
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
call wall_time(wall1)
print *, ' wall time for three_e_5_idx_exch13_bi_ort', wall1 - wall0
call print_memory_usage()
END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, three_e_5_idx_exch12_bi_ort, (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_5_idx_exch12_bi_ort(m,l,j,k,i) = <mlk|-L|mij> ::: 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_5_idx_exch12_bi_ort = 0.d0
print *, ' Providing the three_e_5_idx_exch12_bi_ort ...'
call wall_time(wall0)
provide mos_r_in_r_array_transp mos_l_in_r_array_transp
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,j,k,m,l,integral) &
!$OMP SHARED (mo_num,three_e_5_idx_exch12_bi_ort)
!$OMP DO SCHEDULE (dynamic) COLLAPSE(2)
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
call give_integrals_3_body_bi_ort(m, l, k, m, i, j, integral)
three_e_5_idx_exch12_bi_ort(m,l,j,k,i) = -1.d0 * integral
enddo
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
call wall_time(wall1)
print *, ' wall time for three_e_5_idx_exch12_bi_ort', wall1 - wall0
call print_memory_usage()
END_PROVIDER
! ---

295
src/bi_ort_ints/three_body_ijmkl_old.irp.f Normal file
View File

@ -0,0 +1,295 @@
! ---
BEGIN_PROVIDER [ double precision, three_e_5_idx_direct_bi_ort_old, (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_5_idx_direct_bi_ort_old(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_5_idx_direct_bi_ort_old = 0.d0
print *, ' Providing the three_e_5_idx_direct_bi_ort_old ...'
call wall_time(wall0)
provide mos_r_in_r_array_transp mos_l_in_r_array_transp
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,j,k,m,l,integral) &
!$OMP SHARED (mo_num,three_e_5_idx_direct_bi_ort_old)
!$OMP DO SCHEDULE (dynamic) COLLAPSE(2)
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
call give_integrals_3_body_bi_ort(m, l, k, m, j, i, integral)
three_e_5_idx_direct_bi_ort_old(m,l,j,k,i) = -1.d0 * integral
enddo
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
call wall_time(wall1)
print *, ' wall time for three_e_5_idx_direct_bi_ort_old', wall1 - wall0
END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, three_e_5_idx_cycle_1_bi_ort_old, (mo_num, mo_num, mo_num, mo_num, mo_num)]
BEGIN_DOC
!
! matrix element of the -L three-body operator FOR THE FIRST CYCLIC PERMUTATION TERMS OF DOUBLE EXCITATIONS AND BI ORTHO MOs
!
! three_e_5_idx_cycle_1_bi_ort_old(m,l,j,k,i) = <mlk|-L|jim> ::: 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_5_idx_cycle_1_bi_ort_old = 0.d0
print *, ' Providing the three_e_5_idx_cycle_1_bi_ort_old ...'
call wall_time(wall0)
provide mos_r_in_r_array_transp mos_l_in_r_array_transp
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,j,k,m,l,integral) &
!$OMP SHARED (mo_num,three_e_5_idx_cycle_1_bi_ort_old)
!$OMP DO SCHEDULE (dynamic) COLLAPSE(2)
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
call give_integrals_3_body_bi_ort(m, l, k, j, i, m, integral)
three_e_5_idx_cycle_1_bi_ort_old(m,l,j,k,i) = -1.d0 * integral
enddo
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
call wall_time(wall1)
print *, ' wall time for three_e_5_idx_cycle_1_bi_ort_old', wall1 - wall0
END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, three_e_5_idx_cycle_2_bi_ort_old, (mo_num, mo_num, mo_num, mo_num, mo_num)]
BEGIN_DOC
!
! matrix element of the -L three-body operator FOR THE FIRST CYCLIC PERMUTATION TERMS OF DOUBLE EXCITATIONS AND BI ORTHO MOs
!
! three_e_5_idx_cycle_2_bi_ort_old(m,l,j,k,i) = <mlk|-L|imj> ::: 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_5_idx_cycle_2_bi_ort_old = 0.d0
print *, ' Providing the three_e_5_idx_cycle_2_bi_ort_old ...'
call wall_time(wall0)
provide mos_r_in_r_array_transp mos_l_in_r_array_transp
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,j,k,m,l,integral) &
!$OMP SHARED (mo_num,three_e_5_idx_cycle_2_bi_ort_old)
!$OMP DO SCHEDULE (dynamic) COLLAPSE(2)
do i = 1, mo_num
do k = 1, mo_num
do j = 1, mo_num
do m = 1, mo_num
do l = 1, mo_num
call give_integrals_3_body_bi_ort(m, l, k, i, m, j, integral)
three_e_5_idx_cycle_2_bi_ort_old(m,l,j,k,i) = -1.d0 * integral
enddo
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
call wall_time(wall1)
print *, ' wall time for three_e_5_idx_cycle_2_bi_ort_old', wall1 - wall0
END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, three_e_5_idx_exch23_bi_ort_old, (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_5_idx_exch23_bi_ort_old(m,l,j,k,i) = <mlk|-L|jmi> ::: 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_5_idx_exch23_bi_ort_old = 0.d0
print *, ' Providing the three_e_5_idx_exch23_bi_ort_old ...'
call wall_time(wall0)
provide mos_r_in_r_array_transp mos_l_in_r_array_transp
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,j,k,m,l,integral) &
!$OMP SHARED (mo_num,three_e_5_idx_exch23_bi_ort_old)
!$OMP DO SCHEDULE (dynamic) COLLAPSE(2)
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
call give_integrals_3_body_bi_ort(m, l, k, j, m, i, integral)
three_e_5_idx_exch23_bi_ort_old(m,l,j,k,i) = -1.d0 * integral
enddo
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
call wall_time(wall1)
print *, ' wall time for three_e_5_idx_exch23_bi_ort_old', wall1 - wall0
END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, three_e_5_idx_exch13_bi_ort_old, (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_5_idx_exch13_bi_ort_old(m,l,j,k,i) = <mlk|-L|ijm> ::: 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_5_idx_exch13_bi_ort_old = 0.d0
print *, ' Providing the three_e_5_idx_exch13_bi_ort_old ...'
call wall_time(wall0)
provide mos_r_in_r_array_transp mos_l_in_r_array_transp
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,j,k,m,l,integral) &
!$OMP SHARED (mo_num,three_e_5_idx_exch13_bi_ort_old)
!$OMP DO SCHEDULE (dynamic) COLLAPSE(2)
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
call give_integrals_3_body_bi_ort(m, l, k, i, j, m, integral)
three_e_5_idx_exch13_bi_ort_old(m,l,j,k,i) = -1.d0 * integral
enddo
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
call wall_time(wall1)
print *, ' wall time for three_e_5_idx_exch13_bi_ort_old', wall1 - wall0
END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, three_e_5_idx_exch12_bi_ort_old, (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_5_idx_exch12_bi_ort_old(m,l,j,k,i) = <mlk|-L|mij> ::: 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
provide mos_r_in_r_array_transp mos_l_in_r_array_transp
PROVIDE mo_l_coef mo_r_coef int2_grad1_u12_bimo_t
three_e_5_idx_exch12_bi_ort_old = 0.d0
print *, ' Providing the three_e_5_idx_exch12_bi_ort_old ...'
call wall_time(wall0)
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,j,k,m,l,integral) &
!$OMP SHARED (mo_num,three_e_5_idx_exch12_bi_ort_old)
!$OMP DO SCHEDULE (dynamic) COLLAPSE(2)
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
call give_integrals_3_body_bi_ort(m, l, k, m, i, j, integral)
three_e_5_idx_exch12_bi_ort_old(m,l,j,k,i) = -1.d0 * integral
enddo
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
call wall_time(wall1)
print *, ' wall time for three_e_5_idx_exch12_bi_ort_old', wall1 - wall0
END_PROVIDER

11
src/bi_ort_ints/three_body_ints_bi_ort.irp.f
View File

@ -86,22 +86,25 @@ subroutine give_integrals_3_body_bi_ort(n, l, k, m, j, i, integral)
PROVIDE int2_grad1_u12_bimo_t
integral = 0.d0
! (n, l, k, m, j, i)
do ipoint = 1, n_points_final_grid
weight = final_weight_at_r_vector(ipoint)
integral += weight * mos_l_in_r_array_transp(ipoint,k) * mos_r_in_r_array_transp(ipoint,i) &
tmp = mos_l_in_r_array_transp(ipoint,k) * mos_r_in_r_array_transp(ipoint,i) &
* ( int2_grad1_u12_bimo_t(ipoint,1,n,m) * int2_grad1_u12_bimo_t(ipoint,1,l,j) &
+ int2_grad1_u12_bimo_t(ipoint,2,n,m) * int2_grad1_u12_bimo_t(ipoint,2,l,j) &
+ int2_grad1_u12_bimo_t(ipoint,3,n,m) * int2_grad1_u12_bimo_t(ipoint,3,l,j) )
integral += weight * mos_l_in_r_array_transp(ipoint,l) * mos_r_in_r_array_transp(ipoint,j) &
tmp = tmp + mos_l_in_r_array_transp(ipoint,l) * mos_r_in_r_array_transp(ipoint,j) &
* ( int2_grad1_u12_bimo_t(ipoint,1,n,m) * int2_grad1_u12_bimo_t(ipoint,1,k,i) &
+ int2_grad1_u12_bimo_t(ipoint,2,n,m) * int2_grad1_u12_bimo_t(ipoint,2,k,i) &
+ int2_grad1_u12_bimo_t(ipoint,3,n,m) * int2_grad1_u12_bimo_t(ipoint,3,k,i) )
integral += weight * mos_l_in_r_array_transp(ipoint,n) * mos_r_in_r_array_transp(ipoint,m) &
tmp = tmp + mos_l_in_r_array_transp(ipoint,n) * mos_r_in_r_array_transp(ipoint,m) &
* ( int2_grad1_u12_bimo_t(ipoint,1,l,j) * int2_grad1_u12_bimo_t(ipoint,1,k,i) &
+ int2_grad1_u12_bimo_t(ipoint,2,l,j) * int2_grad1_u12_bimo_t(ipoint,2,k,i) &
+ int2_grad1_u12_bimo_t(ipoint,3,l,j) * int2_grad1_u12_bimo_t(ipoint,3,k,i) )
integral = integral + tmp * final_weight_at_r_vector(ipoint)
enddo
end subroutine give_integrals_3_body_bi_ort

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

@ -16,20 +16,16 @@ subroutine run_ccsd_space_orb
double precision, allocatable :: all_err(:,:), all_t(:,:)
integer, allocatable :: list_occ(:), list_vir(:)
integer(bit_kind) :: det(N_int,2)
integer :: nO, nV, nOa, nOb, nVa, nVb, n_spin(4)
integer :: nO, nV, nOa, nVa
PROVIDE mo_two_e_integrals_in_map
! PROVIDE mo_two_e_integrals_in_map
det = psi_det(:,:,cc_ref)
print*,'Reference determinant:'
call print_det(det,N_int)
! Extract number of occ/vir alpha/beta spin orbitals
!call extract_n_spin(det,n_spin)
nOa = cc_nOa !n_spin(1)
nOb = cc_nOb !n_spin(2)
nVa = cc_nVa !n_spin(3)
nVb = cc_nVb !n_spin(4)
nOa = cc_nOa
nVa = cc_nVa
! Check that the reference is a closed shell determinant
if (cc_ref_is_open_shell) then
@ -109,7 +105,7 @@ subroutine run_ccsd_space_orb
call update_t1(nO,nV,cc_space_f_o,cc_space_f_v,r1,t1)
call update_t2(nO,nV,cc_space_f_o,cc_space_f_v,r2,t2)
else
print*,'Unkonw cc_method_method: '//cc_update_method
print*,'Unkown cc_method_method: '//cc_update_method
endif
call update_tau_space(nO,nV,t1,t2,tau)
@ -169,8 +165,13 @@ subroutine run_ccsd_space_orb
! New
print*,'Computing (T) correction...'
call wall_time(ta)
call ccsd_par_t_space_v2(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)
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 &
,cc_space_v_vvvo,cc_space_v_vvoo,cc_space_v_vooo,e_t)
call wall_time(tb)
print*,'Time: ',tb-ta, ' s'
@ -211,8 +212,8 @@ subroutine ccsd_energy_space(nO,nV,tau,t1,energy)
!$omp default(none)
e = 0d0
!$omp do
do i = 1, nO
do a = 1, nV
do i = 1, nO
e = e + 2d0 * cc_space_f_vo(a,i) * t1(i,a)
enddo
enddo
@ -255,7 +256,7 @@ subroutine update_tau_space(nO,nV,t1,t2,tau)
!$OMP SHARED(nO,nV,tau,t2,t1) &
!$OMP PRIVATE(i,j,a,b) &
!$OMP DEFAULT(NONE)
!$OMP DO collapse(3)
!$OMP DO
do b = 1, nV
do a = 1, nV
do j = 1, nO
@ -373,7 +374,7 @@ subroutine compute_r1_space(nO,nV,t1,t2,tau,H_oo,H_vv,H_vo,r1,max_r1)
!$omp shared(nO,nV,X_voov,t2,t1) &
!$omp private(u,beta,i,a) &
!$omp default(none)
!$omp do collapse(3)
!$omp do
do beta = 1, nV
do u = 1, nO
do i = 1, nO
@ -412,7 +413,7 @@ subroutine compute_r1_space(nO,nV,t1,t2,tau,H_oo,H_vv,H_vo,r1,max_r1)
!$omp shared(nO,nV,cc_space_v_ovov,cc_space_v_voov,X_ovov) &
!$omp private(u,beta,i,a) &
!$omp default(none)
!$omp do collapse(3)
!$omp do
do beta = 1, nV
do u = 1, nO
do a = 1, nv
@ -452,7 +453,7 @@ subroutine compute_r1_space(nO,nV,t1,t2,tau,H_oo,H_vv,H_vo,r1,max_r1)
!$omp shared(nO,nV,cc_space_v_vvov,W_vvov,T_vvoo,tau) &
!$omp private(b,beta,i,a) &
!$omp default(none)
!$omp do collapse(3)
!$omp do
do beta = 1, nV
do i = 1, nO
do b = 1, nV
@ -464,11 +465,11 @@ subroutine compute_r1_space(nO,nV,t1,t2,tau,H_oo,H_vv,H_vo,r1,max_r1)
enddo
!$omp end do nowait
!$omp do collapse(3)
!$omp do
do u = 1, nO
do i = 1, nO
do b = 1, nV
do a = 1, nV
do u = 1, nO
T_vvoo(a,b,i,u) = tau(i,u,a,b)
enddo
enddo
@ -504,8 +505,8 @@ subroutine compute_r1_space(nO,nV,t1,t2,tau,H_oo,H_vv,H_vo,r1,max_r1)
!$omp shared(nO,nV,cc_space_v_vooo,W_oovo) &
!$omp private(u,a,i,j) &
!$omp default(none)
!$omp do collapse(3)
do u = 1, nO
!$omp do
do a = 1, nV
do j = 1, nO
do i = 1, nO
@ -513,8 +514,8 @@ subroutine compute_r1_space(nO,nV,t1,t2,tau,H_oo,H_vv,H_vo,r1,max_r1)
enddo
enddo
enddo
!$omp end do nowait
enddo
!$omp end do
!$omp end parallel
call dgemm('T','N', nO, nV, nO*nO*nV, &
@ -527,9 +528,7 @@ subroutine compute_r1_space(nO,nV,t1,t2,tau,H_oo,H_vv,H_vo,r1,max_r1)
max_r1 = 0d0
do a = 1, nV
do i = 1, nO
if (dabs(r1(i,a)) > max_r1) then
max_r1 = dabs(r1(i,a))
endif
max_r1 = max(dabs(r1(i,a)), max_r1)
enddo
enddo
@ -657,7 +656,7 @@ subroutine compute_H_vv(nO,nV,t1,t2,tau,H_vv)
! H_vv(a,beta) = H_vv(a,beta) - cc_space_w_vvoo(a,b,i,j) * tau(i,j,beta,b)
! H_vv(a,beta) = H_vv(a,beta) - cc_space_w_vvoo(a,b,i,j) * tmp_tau(b,i,j,beta)
!$omp do collapse(3)
!$omp do
do beta = 1, nV
do j = 1, nO
do i = 1, nO
@ -727,7 +726,7 @@ subroutine compute_H_vo(nO,nV,t1,t2,H_vo)
! H_vo(a,i) = H_vo(a,i) + cc_space_w_vvoo(a,b,i,j) * t1(j,b)
! H_vo(a,i) = H_vo(a,i) + w(a,i,j,b) * t1(j,b)
!$omp do collapse(3)
!$omp do
do b = 1, nV
do j = 1, nO
do i = 1, nO
@ -765,7 +764,7 @@ subroutine compute_r2_space(nO,nV,t1,t2,tau,H_oo,H_vv,H_vo,r2,max_r2)
! internal
double precision, allocatable :: g_occ(:,:), g_vir(:,:), J1(:,:,:,:), K1(:,:,:,:)
double precision, allocatable :: A1(:,:,:,:), B1(:,:,:,:)
double precision, allocatable :: A1(:,:,:,:), B1_gam(:,:,:)
integer :: u,v,i,j,beta,gam,a,b
allocate(g_occ(nO,nO), g_vir(nV,nV))
@ -787,7 +786,7 @@ subroutine compute_r2_space(nO,nV,t1,t2,tau,H_oo,H_vv,H_vo,r2,max_r2)
!$omp shared(nO,nV,r2,cc_space_v_oovv) &
!$omp private(u,v,gam,beta) &
!$omp default(none)
!$omp do collapse(3)
!$omp do
do gam = 1, nV
do beta = 1, nV
do v = 1, nO
@ -835,13 +834,18 @@ subroutine compute_r2_space(nO,nV,t1,t2,tau,H_oo,H_vv,H_vo,r2,max_r2)
! enddo
!enddo
allocate(B1(nV,nV,nV,nV))
call compute_B1(nO,nV,t1,t2,B1)
call dgemm('N','N',nO*nO,nV*nV,nV*nV, &
! allocate(B1(nV,nV,nV,nV))
! call compute_B1(nO,nV,t1,t2,B1)
allocate(B1_gam(nV,nV,nV))
do gam=1,nV
call compute_B1_gam(nO,nV,t1,t2,B1_gam,gam)
call dgemm('N','N',nO*nO,nV,nV*nV, &
1d0, tau, size(tau,1) * size(tau,2), &
B1 , size(B1,1) * size(B1,2), &
1d0, r2, size(r2,1) * size(r2,2))
deallocate(B1)
B1_gam , size(B1_gam,1) * size(B1_gam,2), &
1d0, r2(1,1,1,gam), size(r2,1) * size(r2,2))
enddo
deallocate(B1_gam)
!do gam = 1, nV
! do beta = 1, nV
@ -863,7 +867,7 @@ subroutine compute_r2_space(nO,nV,t1,t2,tau,H_oo,H_vv,H_vo,r2,max_r2)
!$omp shared(nO,nV,t2,X_oovv) &
!$omp private(u,v,gam,a) &
!$omp default(none)
!$omp do collapse(3)
!$omp do
do a = 1, nV
do gam = 1, nV
do v = 1, nO
@ -885,7 +889,7 @@ subroutine compute_r2_space(nO,nV,t1,t2,tau,H_oo,H_vv,H_vo,r2,max_r2)
!$omp shared(nO,nV,r2,Y_oovv) &
!$omp private(u,v,gam,beta) &
!$omp default(none)
!$omp do collapse(3)
!$omp do
do gam = 1, nV
do beta = 1, nV
do v = 1, nO
@ -921,7 +925,7 @@ subroutine compute_r2_space(nO,nV,t1,t2,tau,H_oo,H_vv,H_vo,r2,max_r2)
!$omp shared(nO,nV,r2,X_oovv) &
!$omp private(u,v,gam,beta) &
!$omp default(none)
!$omp do collapse(3)
!$omp do
do gam = 1, nV
do beta = 1, nV
do v = 1, nO
@ -957,7 +961,7 @@ subroutine compute_r2_space(nO,nV,t1,t2,tau,H_oo,H_vv,H_vo,r2,max_r2)
!$omp shared(nO,nV,X_vovv,cc_space_v_ovvv) &
!$omp private(u,a,gam,beta) &
!$omp default(none)
!$omp do collapse(3)
!$omp do
do gam = 1, nV
do beta = 1, nV
do u = 1, nO
@ -979,7 +983,7 @@ subroutine compute_r2_space(nO,nV,t1,t2,tau,H_oo,H_vv,H_vo,r2,max_r2)
!$omp shared(nO,nV,r2,Y_oovv) &
!$omp private(u,v,gam,beta) &
!$omp default(none)
!$omp do collapse(3)
!$omp do
do gam = 1, nV
do beta = 1, nV
do v = 1, nO
@ -1014,8 +1018,8 @@ subroutine compute_r2_space(nO,nV,t1,t2,tau,H_oo,H_vv,H_vo,r2,max_r2)
!$omp shared(nO,nV,X_vovo,cc_space_v_ovov) &
!$omp private(u,v,gam,i) &
!$omp default(none)
!$omp do collapse(3)
do i = 1, nO
!$omp do
do gam = 1, nV
do u = 1, nO
do a = 1, nV
@ -1023,8 +1027,8 @@ subroutine compute_r2_space(nO,nV,t1,t2,tau,H_oo,H_vv,H_vo,r2,max_r2)
enddo
enddo
enddo
!$omp end do nowait
enddo
!$omp end do
!$omp end parallel
call dgemm('N','N',nV*nO*nV,nV,nO, &
@ -1041,7 +1045,7 @@ subroutine compute_r2_space(nO,nV,t1,t2,tau,H_oo,H_vv,H_vo,r2,max_r2)
!$omp shared(nO,nV,r2,X_oovv) &
!$omp private(u,v,gam,beta) &
!$omp default(none)
!$omp do collapse(3)
!$omp do
do gam = 1, nV
do beta = 1, nV
do v = 1, nO
@ -1079,7 +1083,7 @@ subroutine compute_r2_space(nO,nV,t1,t2,tau,H_oo,H_vv,H_vo,r2,max_r2)
!$omp shared(nO,nV,r2,X_oovv) &
!$omp private(u,v,gam,beta) &
!$omp default(none)
!$omp do collapse(3)
!$omp do
do gam = 1, nV
do beta = 1, nV
do v = 1, nO
@ -1116,8 +1120,8 @@ subroutine compute_r2_space(nO,nV,t1,t2,tau,H_oo,H_vv,H_vo,r2,max_r2)
!$omp shared(nO,nV,X_vovo,cc_space_v_ovvo) &
!$omp private(a,v,gam,i) &
!$omp default(none)
!$omp do collapse(3)
do i = 1, nO
!$omp do
do gam = 1, nV
do v = 1, nO
do a = 1, nV
@ -1125,8 +1129,8 @@ subroutine compute_r2_space(nO,nV,t1,t2,tau,H_oo,H_vv,H_vo,r2,max_r2)
enddo
enddo
enddo
!$omp end do nowait
enddo
!$omp end do
!$omp end parallel
call dgemm('N','N',nO,nO*nV*nO,nV, &
@ -1143,7 +1147,7 @@ subroutine compute_r2_space(nO,nV,t1,t2,tau,H_oo,H_vv,H_vo,r2,max_r2)
!$omp shared(nO,nV,r2,X_oovv) &
!$omp private(u,v,gam,beta) &
!$omp default(none)
!$omp do collapse(3)
!$omp do
do gam = 1, nV
do beta = 1, nV
do v = 1, nO
@ -1182,19 +1186,19 @@ subroutine compute_r2_space(nO,nV,t1,t2,tau,H_oo,H_vv,H_vo,r2,max_r2)
!$omp shared(nO,nV,X_ovvo,Y_voov,K1,J1,t2) &
!$omp private(u,v,gam,beta,i,a) &
!$omp default(none)
!$omp do collapse(3)
do i = 1, nO
!$omp do
do a = 1, nV
do beta = 1, nV
do u = 1, nO
X_ovvo(u,beta,a,i) = 0.5d0 * (2d0 * J1(u,a,beta,i) - K1(u,a,i,beta))
enddo
X_ovvo(u,beta,a,i) = (J1(u,a,beta,i) - 0.5d0 * K1(u,a,i,beta))
enddo
enddo
enddo
!$omp end do nowait
enddo
!$omp do collapse(3)
!$omp do
do gam = 1, nV
do v = 1, nO
do i = 1, nO
@ -1216,7 +1220,7 @@ subroutine compute_r2_space(nO,nV,t1,t2,tau,H_oo,H_vv,H_vo,r2,max_r2)
!$omp shared(nO,nV,r2,Z_ovov) &
!$omp private(u,v,gam,beta) &
!$omp default(none)
!$omp do collapse(3)
!$omp do
do gam = 1, nV
do beta = 1, nV
do v = 1, nO
@ -1252,7 +1256,7 @@ subroutine compute_r2_space(nO,nV,t1,t2,tau,H_oo,H_vv,H_vo,r2,max_r2)
!$omp shared(nO,nV,r2,K1,X_ovov,Y_ovov,t2) &
!$omp private(u,a,i,beta,gam) &
!$omp default(none)
!$omp do collapse(3)
!$omp do
do beta = 1, nV
do u = 1, nO
do a = 1, nV
@ -1264,7 +1268,7 @@ subroutine compute_r2_space(nO,nV,t1,t2,tau,H_oo,H_vv,H_vo,r2,max_r2)
enddo
!$omp end do nowait
!$omp do collapse(3)
!$omp do
do gam = 1, nV
do v = 1, nO
do a = 1, nV
@ -1286,7 +1290,7 @@ subroutine compute_r2_space(nO,nV,t1,t2,tau,H_oo,H_vv,H_vo,r2,max_r2)
!$omp shared(nO,nV,r2,Z_ovov) &
!$omp private(u,v,gam,beta) &
!$omp default(none)
!$omp do collapse(3)
!$omp do
do gam = 1, nV
do beta = 1, nV
do v = 1, nO
@ -1319,7 +1323,7 @@ subroutine compute_r2_space(nO,nV,t1,t2,tau,H_oo,H_vv,H_vo,r2,max_r2)
!$omp shared(nO,nV,K1,X_ovov,Z_ovov,t2) &
!$omp private(u,v,gam,beta,i,a) &
!$omp default(none)
!$omp do collapse(3)
!$omp do
do a = 1, nV
do i = 1, nO
do gam = 1, nV
@ -1331,7 +1335,7 @@ subroutine compute_r2_space(nO,nV,t1,t2,tau,H_oo,H_vv,H_vo,r2,max_r2)
enddo
!$omp end do nowait
!$omp do collapse(3)
!$omp do
do beta = 1, nV
do v = 1, nO
do a = 1, nV
@ -1353,7 +1357,7 @@ subroutine compute_r2_space(nO,nV,t1,t2,tau,H_oo,H_vv,H_vo,r2,max_r2)
!$omp shared(nO,nV,r2,Z_ovov) &
!$omp private(u,v,gam,beta) &
!$omp default(none)
!$omp do collapse(3)
!$omp do
do gam = 1, nV
do beta = 1, nV
do v = 1, nO
@ -1373,7 +1377,7 @@ subroutine compute_r2_space(nO,nV,t1,t2,tau,H_oo,H_vv,H_vo,r2,max_r2)
!$omp shared(nO,nV,r2) &
!$omp private(i,j,a,b) &
!$omp default(none)
!$omp do collapse(3)
!$omp do
do b = 1, nV
do a = 1, nV
do j = 1, nO
@ -1391,9 +1395,7 @@ subroutine compute_r2_space(nO,nV,t1,t2,tau,H_oo,H_vv,H_vo,r2,max_r2)
do a = 1, nV
do j = 1, nO
do i = 1, nO
if (dabs(r2(i,j,a,b)) > max_r2) then
max_r2 = dabs(r2(i,j,a,b))
endif
max_r2 = max(r2(i,j,a,b), max_r2)
enddo
enddo
enddo
@ -1448,7 +1450,7 @@ subroutine compute_A1(nO,nV,t1,t2,tau,A1)
!$omp shared(nO,nV,A1,cc_space_v_oooo,cc_space_v_ovoo,X_vooo) &
!$omp private(u,v,i,j) &
!$omp default(none)
!$omp do collapse(3)
!$omp do collapse(2)
do j = 1, nO
do i = 1, nO
do v = 1, nO
@ -1462,7 +1464,7 @@ subroutine compute_A1(nO,nV,t1,t2,tau,A1)
! A1(u,v,i,j) += cc_space_v_ovoo(u,a,i,j) * t1(v,a) &
!$omp do collapse(3)
!$omp do collapse(2)
do j = 1, nO
do i = 1, nO
do u = 1, nO
@ -1484,7 +1486,7 @@ subroutine compute_A1(nO,nV,t1,t2,tau,A1)
!$omp shared(nO,nV,A1,Y_oooo) &
!$omp private(u,v,i,j) &
!$omp default(none)
!$omp do collapse(3)
!$omp do collapse(2)
do j = 1, nO
do i = 1, nO
do v = 1, nO
@ -1515,6 +1517,90 @@ end
! B1
subroutine compute_B1_gam(nO,nV,t1,t2,B1,gam)
implicit none
integer, intent(in) :: nO,nV,gam
double precision, intent(in) :: t1(nO, nV)
double precision, intent(in) :: t2(nO, nO, nV, nV)
double precision, intent(out) :: B1(nV, nV, nV)
integer :: a,tmp_a,b,k,l,c,d,tmp_c,tmp_d,i,j,u,v, beta
! do beta = 1, nV
! do b = 1, nV
! do a = 1, nV
! B1(a,b,beta) = cc_space_v_vvvv(a,b,beta,gam)
!
! do i = 1, nO
! B1(a,b,beta) = B1(a,b,beta) &
! - cc_space_v_vvvo(a,b,beta,i) * t1(i,gam) &
! - cc_space_v_vvov(a,b,i,gam) * t1(i,beta)
! enddo
!
! enddo
! enddo
! enddo
double precision, allocatable :: X_vvvo(:,:,:), Y_vvvv(:,:,:)
allocate(X_vvvo(nV,nV,nO), Y_vvvv(nV,nV,nV))
! ! B1(a,b,beta,gam) = cc_space_v_vvvv(a,b,beta,gam)
!$omp parallel &
!$omp shared(nO,nV,B1,cc_space_v_vvvv,cc_space_v_vvov,X_vvvo,gam) &
!$omp private(a,b,beta) &
!$omp default(none)
!$omp do
do beta = 1, nV
do b = 1, nV
do a = 1, nV
B1(a,b,beta) = cc_space_v_vvvv(a,b,beta,gam)
enddo
enddo
enddo
!$omp end do nowait
do i = 1, nO
!$omp do
do b = 1, nV
do a = 1, nV
X_vvvo(a,b,i) = cc_space_v_vvov(a,b,i,gam)
enddo
enddo
!$omp end do nowait
enddo
!$omp end parallel
! ! B1(a,b,beta) -= cc_space_v_vvvo(a,b,beta,i) * t1(i,gam) &
call dgemm('N','N', nV*nV*nV, 1, nO, &
-1d0, cc_space_v_vvvo, size(cc_space_v_vvvo,1) * size(cc_space_v_vvvo,2) * size(cc_space_v_vvvo,3), &
t1(1,gam), size(t1,1), &
1d0, B1 , size(B1,1) * size(B1,2) * size(B1,3))
! B1(a,b,beta,gam) -= cc_space_v_vvov(a,b,i,gam) * t1(i,beta)
call dgemm('N','N', nV*nV, nV, nO, &
-1d0, X_vvvo, size(X_vvvo,1) * size(X_vvvo,2), &
t1 , size(t1,1), &
0d0, Y_vvvv, size(Y_vvvv,1) * size(Y_vvvv,2))
!$omp parallel &
!$omp shared(nV,B1,Y_vvvv,gam) &
!$omp private(a,b,beta) &
!$omp default(none)
!$omp do
do beta = 1, nV
do b = 1, nV
do a = 1, nV
B1(a,b,beta) = B1(a,b,beta) + Y_vvvv(a,b,beta)
enddo
enddo
enddo
!$omp end do
!$omp end parallel
deallocate(X_vvvo,Y_vvvv)
end
subroutine compute_B1(nO,nV,t1,t2,B1)
implicit none
@ -1553,7 +1639,7 @@ subroutine compute_B1(nO,nV,t1,t2,B1)
!$omp shared(nO,nV,B1,cc_space_v_vvvv,cc_space_v_vvov,X_vvvo) &
!$omp private(a,b,beta,gam) &
!$omp default(none)
!$omp do collapse(3)
!$omp do
do gam = 1, nV
do beta = 1, nV
do b = 1, nV
@ -1564,8 +1650,8 @@ subroutine compute_B1(nO,nV,t1,t2,B1)
enddo
enddo
!$omp end do nowait
!$omp do collapse(3)
do i = 1, nO
!$omp do
do gam = 1, nV
do b = 1, nV
do a = 1, nV
@ -1573,8 +1659,8 @@ subroutine compute_B1(nO,nV,t1,t2,B1)
enddo
enddo
enddo
!$omp end do nowait
enddo
!$omp end do
!$omp end parallel
! B1(a,b,beta,gam) -= cc_space_v_vvvo(a,b,beta,i) * t1(i,gam) &
@ -1594,7 +1680,7 @@ subroutine compute_B1(nO,nV,t1,t2,B1)
!$omp shared(nV,B1,Y_vvvv) &
!$omp private(a,b,beta,gam) &
!$omp default(none)
!$omp do collapse(3)
!$omp do
do gam = 1, nV
do beta = 1, nV
do b = 1, nV
@ -1658,7 +1744,7 @@ subroutine compute_g_occ(nO,nV,t1,t2,H_oo,g_occ)
enddo
!$omp end do
!$omp do collapse(1)
!$omp do
do i = 1, nO
do j = 1, nO
do a = 1, nV
@ -1720,7 +1806,7 @@ subroutine compute_g_vir(nO,nV,t1,t2,H_vv,g_vir)
enddo
!$omp end do
!$omp do collapse(1)
!$omp do
do beta = 1, nV
do i = 1, nO
do b = 1, nV
@ -1788,8 +1874,8 @@ subroutine compute_J1(nO,nV,t1,t2,v_ovvo,v_ovoo,v_vvvo,v_vvoo,J1)
!$omp shared(nO,nV,J1,v_ovvo,v_ovoo,X_ovoo) &
!$omp private(i,j,a,u,beta) &
!$omp default(none)
!$omp do collapse(3)
do i = 1, nO
!$omp do
do beta = 1, nV
do a = 1, nV
do u = 1, nO
@ -1797,10 +1883,10 @@ subroutine compute_J1(nO,nV,t1,t2,v_ovvo,v_ovoo,v_vvvo,v_vvoo,J1)
enddo
enddo
enddo
enddo
!$omp end do nowait
enddo
!$omp do collapse(3)
!$omp do collapse(2)
do j = 1, nO
do i = 1, nO
do a = 1, nV
@ -1822,8 +1908,8 @@ subroutine compute_J1(nO,nV,t1,t2,v_ovvo,v_ovoo,v_vvvo,v_vvoo,J1)
!$omp shared(nO,nV,J1,Y_ovov) &
!$omp private(i,beta,a,u) &
!$omp default(none)
!$omp do collapse(3)
do i = 1, nO
!$omp do
do beta = 1, nV
do a = 1, nV
do u = 1, nO
@ -1831,8 +1917,8 @@ subroutine compute_J1(nO,nV,t1,t2,v_ovvo,v_ovoo,v_vvvo,v_vvoo,J1)
enddo
enddo
enddo
!$omp end do nowait
enddo
!$omp end do
!$omp end parallel
deallocate(X_ovoo)
@ -1849,7 +1935,7 @@ subroutine compute_J1(nO,nV,t1,t2,v_ovvo,v_ovoo,v_vvvo,v_vvoo,J1)
!$omp shared(nO,nV,t2,t1,Y_ovov,X_voov,v_vvoo) &
!$omp private(i,beta,a,u,b,j) &
!$omp default(none)
!$omp do collapse(3)
!$omp do
do b = 1, nV
do j = 1, nO
do beta = 1, nV
@ -1861,7 +1947,7 @@ subroutine compute_J1(nO,nV,t1,t2,v_ovvo,v_ovoo,v_vvvo,v_vvoo,J1)
enddo
!$omp end do nowait
!$omp do collapse(3)
!$omp do
do b = 1, nV
do j = 1, nO
do i = 1, nO
@ -1886,8 +1972,8 @@ subroutine compute_J1(nO,nV,t1,t2,v_ovvo,v_ovoo,v_vvvo,v_vvoo,J1)
!$omp shared(nO,nV,J1,Z_ovvo,t2,Y_vovo,v_vvoo,X_ovvo) &
!$omp private(i,beta,a,u,j,b) &
!$omp default(none)
!$omp do collapse(3)
do i = 1, nO
!$omp do
do beta = 1, nV
do a = 1, nV
do u = 1, nO
@ -1895,12 +1981,12 @@ subroutine compute_J1(nO,nV,t1,t2,v_ovvo,v_ovoo,v_vvvo,v_vvoo,J1)
enddo
enddo
enddo
enddo
!$omp end do nowait
enddo
!+ 0.5d0 * (2d0 * cc_space_v_vvoo(a,b,i,j) - cc_space_v_vvoo(b,a,i,j)) * t2(u,j,beta,b)
!$omp do collapse(3)
do j = 1, nO
!$omp do
do b = 1, nV
do i = 1, nO
do a = 1, nV
@ -1908,11 +1994,11 @@ subroutine compute_J1(nO,nV,t1,t2,v_ovvo,v_ovoo,v_vvvo,v_vvoo,J1)
enddo
enddo
enddo
enddo
!$omp end do nowait
enddo
!$omp do collapse(3)
do j = 1, nO
!$omp do
do b = 1, nV
do beta = 1, nV
do u = 1, nO
@ -1920,8 +2006,8 @@ subroutine compute_J1(nO,nV,t1,t2,v_ovvo,v_ovoo,v_vvvo,v_vvoo,J1)
enddo
enddo
enddo
!$omp end do nowait
enddo
!$omp end do
!$omp end parallel
call dgemm('N','T',nO*nV,nV*nO,nV*nO, &
@ -1933,8 +2019,8 @@ subroutine compute_J1(nO,nV,t1,t2,v_ovvo,v_ovoo,v_vvvo,v_vvoo,J1)
!$omp shared(nO,nV,J1,Z_ovvo) &
!$omp private(i,beta,a,u) &
!$omp default(none)
!$omp do collapse(3)
do i = 1, nO
!$omp do
do beta = 1, nV
do a = 1, nV
do u = 1, nO
@ -1942,8 +2028,8 @@ subroutine compute_J1(nO,nV,t1,t2,v_ovvo,v_ovoo,v_vvvo,v_vvoo,J1)
enddo
enddo
enddo
!$omp end do nowait
enddo
!$omp end do
!$omp end parallel
deallocate(X_ovvo,Z_ovvo,Y_ovov)
@ -2003,7 +2089,7 @@ subroutine compute_K1(nO,nV,t1,t2,v_ovoo,v_vvoo,v_ovov,v_vvov,K1)
!$omp shared(nO,nV,K1,X,Y,v_vvoo,v_ovov,t1,t2) &
!$omp private(i,beta,a,u,j,b) &
!$omp default(none)
!$omp do collapse(3)
!$omp do
do beta = 1, nV
do i = 1, nO
do a = 1, nV
@ -2015,8 +2101,8 @@ subroutine compute_K1(nO,nV,t1,t2,v_ovoo,v_vvoo,v_ovov,v_vvov,K1)
enddo
!$omp end do nowait
!$omp do collapse(3)
do i = 1, nO
!$omp do
do a = 1, nV
do j = 1, nO
do b = 1, nV
@ -2024,11 +2110,11 @@ subroutine compute_K1(nO,nV,t1,t2,v_ovoo,v_vvoo,v_ovov,v_vvov,K1)
enddo
enddo
enddo
enddo
!$omp end do nowait
enddo
!$omp do collapse(3)
do j = 1, nO
!$omp do
do b = 1, nV
do beta = 1, nV
do u = 1, nO
@ -2036,8 +2122,8 @@ subroutine compute_K1(nO,nV,t1,t2,v_ovoo,v_vvoo,v_ovov,v_vvov,K1)
enddo
enddo
enddo
enddo
!$omp end do
enddo
!$omp end parallel
call dgemm('N','N',nO*nV*nO,nV,nO, &
@ -2060,7 +2146,7 @@ subroutine compute_K1(nO,nV,t1,t2,v_ovoo,v_vvoo,v_ovov,v_vvov,K1)
!$omp shared(nO,nV,K1,Z) &
!$omp private(i,beta,a,u) &
!$omp default(none)
!$omp do collapse(3)
!$omp do
do beta = 1, nV
do i = 1, nO
do a = 1, nV

464
src/ccsd/ccsd_t_space_orb_abc.irp.f
View File

@ -10,51 +10,43 @@ subroutine ccsd_par_t_space_v3(nO,nV,t1,t2,f_o,f_v,v_vvvo,v_vvoo,v_vooo,energy)
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(out) :: energy
double precision, allocatable :: W(:,:,:,:,:,:)
double precision, allocatable :: V(:,:,:,:,:,:)
double precision, allocatable :: W_abc(:,:,:), V_abc(:,:,:)
double precision, allocatable :: W_cab(:,:,:), W_cba(:,:,:)
double precision, allocatable :: W_bca(:,:,:), V_cba(:,:,:)
double precision, allocatable :: X_vvvo(:,:,:,:), X_ovoo(:,:,:,:), X_vvoo(:,:,:,:)
double precision, allocatable :: T_vvoo(:,:,:,:), T_ovvo(:,:,:,:), T_vo(:,:)
double precision, allocatable :: X_vovv(:,:,:,:), X_ooov(:,:,:,:), X_oovv(:,:,:,:)
double precision, allocatable :: T_voov(:,:,:,:), T_oovv(:,:,:,:)
integer :: i,j,k,l,a,b,c,d
double precision :: e,ta,tb, delta, delta_abc
double precision :: e,ta,tb
!allocate(W(nV,nV,nV,nO,nO,nO))
!allocate(V(nV,nV,nV,nO,nO,nO))
allocate(W_abc(nO,nO,nO), V_abc(nO,nO,nO), W_cab(nO,nO,nO))
allocate(W_bca(nO,nO,nO), V_cba(nO,nO,nO), W_cba(nO,nO,nO))
allocate(X_vvvo(nV,nV,nV,nO), X_ovoo(nO,nV,nO,nO), X_vvoo(nV,nV,nO,nO))
allocate(T_vvoo(nV,nV,nO,nO), T_ovvo(nO,nV,nV,nO), T_vo(nV,nO))
call set_multiple_levels_omp(.False.)
allocate(X_vovv(nV,nO,nV,nV), X_ooov(nO,nO,nO,nV), X_oovv(nO,nO,nV,nV))
allocate(T_voov(nV,nO,nO,nV),T_oovv(nO,nO,nV,nV))
! Temporary arrays
!$OMP PARALLEL &
!$OMP SHARED(nO,nV,T_vvoo,T_ovvo,T_vo,X_vvvo,X_ovoo,X_vvoo, &
!$OMP SHARED(nO,nV,T_voov,T_oovv,X_vovv,X_ooov,X_oovv, &
!$OMP t1,t2,v_vvvo,v_vooo,v_vvoo) &
!$OMP PRIVATE(a,b,c,d,i,j,k,l) &
!$OMP DEFAULT(NONE)
!v_vvvo(b,a,d,i) * t2(k,j,c,d) &
!X_vvvo(d,b,a,i) * T_vvoo(d,c,k,j)
!X_vovv(d,i,b,a,i) * T_voov(d,j,c,k)
!$OMP DO collapse(3)
do i = 1, nO
!$OMP DO
do a = 1, nV
do b = 1, nV
do i = 1, nO
do d = 1, nV
X_vvvo(d,b,a,i) = v_vvvo(b,a,d,i)
X_vovv(d,i,b,a) = v_vvvo(b,a,d,i)
enddo
enddo
enddo
enddo
!$OMP END DO nowait
!$OMP DO collapse(3)
!$OMP DO
do c = 1, nV
do j = 1, nO
do k = 1, nO
do c = 1, nV
do d = 1, nV
T_vvoo(d,c,k,j) = t2(k,j,c,d)
T_voov(d,k,j,c) = t2(k,j,c,d)
enddo
enddo
enddo
@ -62,191 +54,399 @@ subroutine ccsd_par_t_space_v3(nO,nV,t1,t2,f_o,f_v,v_vvvo,v_vvoo,v_vooo,energy)
!$OMP END DO nowait
!v_vooo(c,j,k,l) * t2(i,l,a,b) &
!X_ovoo(l,c,j,k) * T_ovvo(l,a,b,i) &
!X_ooov(l,j,k,c) * T_oovv(l,i,a,b) &
!$OMP DO collapse(3)
!$OMP DO
do c = 1, nV
do k = 1, nO
do j = 1, nO
do c = 1, nV
do l = 1, nO
X_ovoo(l,c,j,k) = v_vooo(c,j,k,l)
X_ooov(l,j,k,c) = v_vooo(c,j,k,l)
enddo
enddo
enddo
enddo
!$OMP END DO nowait
!$OMP DO collapse(3)
do i = 1, nO
!$OMP DO
do b = 1, nV
do a = 1, nV
do i = 1, nO
do l = 1, nO
T_ovvo(l,a,b,i) = t2(i,l,a,b)
T_oovv(l,i,a,b) = t2(i,l,a,b)
enddo
enddo
enddo
enddo
!$OMP END DO nowait
!v_vvoo(b,c,j,k) * t1(i,a) &
!X_vvoo(b,c,k,j) * T1_vo(a,i) &
!X_oovv(j,k,b,c) * T1_vo(a,i) &
!$OMP DO collapse(3)
do j = 1, nO
!$OMP DO
do c = 1, nV
do b = 1, nV
do k = 1, nO
do c = 1, nV
do b = 1, nV
X_vvoo(b,c,k,j) = v_vvoo(b,c,j,k)
do j = 1, nO
X_oovv(j,k,b,c) = v_vvoo(b,c,j,k)
enddo
enddo
enddo
enddo
!$OMP END DO nowait
!$OMP DO collapse(1)
do i = 1, nO
do a = 1, nV
T_vo(a,i) = t1(i,a)
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
call wall_time(ta)
energy = 0d0
do c = 1, nV
do b = 1, nV
do a = 1, nV
delta_abc = f_v(a) + f_v(b) + f_v(c)
call form_w_abc(nO,nV,a,b,c,T_vvoo,T_ovvo,X_vvvo,X_ovoo,W_abc)
call form_w_abc(nO,nV,b,c,a,T_vvoo,T_ovvo,X_vvvo,X_ovoo,W_bca)
call form_w_abc(nO,nV,c,a,b,T_vvoo,T_ovvo,X_vvvo,X_ovoo,W_cab)
call form_w_abc(nO,nV,c,b,a,T_vvoo,T_ovvo,X_vvvo,X_ovoo,W_cba)
double precision, external :: ccsd_t_task_aba
double precision, external :: ccsd_t_task_abc
call form_v_abc(nO,nV,a,b,c,T_vo,X_vvoo,W_abc,V_abc)
call form_v_abc(nO,nV,c,b,a,T_vo,X_vvoo,W_cba,V_cba)
!$OMP PARALLEL &
!$OMP SHARED(energy,nO,a,b,c,W_abc,W_cab,W_bca,V_abc,V_cba,f_o,f_v,delta_abc)&
!$OMP PRIVATE(i,j,k,e,delta) &
!$OMP DEFAULT(NONE)
!$OMP PARALLEL PRIVATE(a,b,c,e) DEFAULT(SHARED)
e = 0d0
!$OMP DO
do i = 1, nO
do j = 1, nO
do k = 1, nO
delta = 1d0 / (f_o(i) + f_o(j) + f_o(k) - delta_abc)
!energy = energy + (4d0 * W(i,j,k,a,b,c) + W(i,j,k,b,c,a) + W(i,j,k,c,a,b)) * (V(i,j,k,a,b,c) - V(i,j,k,c,b,a)) / (cc_space_f_o(i) + cc_space_f_o(j) + cc_space_f_o(k) - cc_space_f_v(a) - cc_space_f_v(b) - cc_space_f_v(c)) !delta_ooovvv(i,j,k,a,b,c)
e = e + (4d0 * W_abc(i,j,k) + W_bca(i,j,k) + W_cab(i,j,k))&
* (V_abc(i,j,k) - V_cba(i,j,k)) * delta
!$OMP DO SCHEDULE(dynamic)
do a = 1, nV
do b = a+1, nV
do c = b+1, nV
e = e + ccsd_t_task_abc(a,b,c,nO,nV,t1,T_oovv,T_voov, &
X_ooov,X_oovv,X_vovv,f_o,f_v)
enddo
e = e + ccsd_t_task_aba(a,b,nO,nV,t1,T_oovv,T_voov, &
X_ooov,X_oovv,X_vovv,f_o,f_v)
e = e + ccsd_t_task_aba(b,a,nO,nV,t1,T_oovv,T_voov, &
X_ooov,X_oovv,X_vovv,f_o,f_v)
enddo
enddo
!$OMP END DO NOWAIT
!$OMP CRITICAL
energy = energy + e
!$OMP END CRITICAL
!$OMP END PARALLEL
enddo
enddo
call wall_time(tb)
write(*,'(F12.2,A5,F12.2,A2)') dble(i)/dble(nO)*100d0, '% in ', tb - ta, ' s'
enddo
energy = energy / 3d0
energy = energy / 3.d0
deallocate(W_abc,V_abc,W_cab,V_cba,W_bca,X_vvvo,X_ovoo,T_vvoo,T_ovvo,T_vo)
!deallocate(V,W)
deallocate(X_vovv,X_ooov,T_voov,T_oovv)
end
subroutine form_w_abc(nO,nV,a,b,c,T_vvoo,T_ovvo,X_vvvo,X_ovoo,W_abc)
double precision function ccsd_t_task_abc(a,b,c,nO,nV,t1,T_oovv,T_voov,&
X_ooov,X_oovv,X_vovv,f_o,f_v) result(e)
implicit none
integer, intent(in) :: nO,nV,a,b,c
double precision, intent(in) :: t1(nO,nV), f_o(nO), f_v(nV)
double precision, intent(in) :: X_oovv(nO,nO,nV,nV)
double precision, intent(in) :: T_voov(nV,nO,nO,nV), T_oovv(nO,nO,nV,nV)
double precision, intent(in) :: X_vovv(nV,nO,nV,nV), X_ooov(nO,nO,nO,nV)
double precision :: delta, delta_abc
integer :: i,j,k
double precision, allocatable :: W_abc(:,:,:), W_cab(:,:,:), W_bca(:,:,:)
double precision, allocatable :: W_bac(:,:,:), W_cba(:,:,:), W_acb(:,:,:)
double precision, allocatable :: V_abc(:,:,:), V_cab(:,:,:), V_bca(:,:,:)
double precision, allocatable :: V_bac(:,:,:), V_cba(:,:,:), V_acb(:,:,:)
allocate( W_abc(nO,nO,nO), W_cab(nO,nO,nO), W_bca(nO,nO,nO), &
W_bac(nO,nO,nO), W_cba(nO,nO,nO), W_acb(nO,nO,nO), &
V_abc(nO,nO,nO), V_cab(nO,nO,nO), V_bca(nO,nO,nO), &
V_bac(nO,nO,nO), V_cba(nO,nO,nO), V_acb(nO,nO,nO) )
call form_w_abc(nO,nV,a,b,c,T_voov,T_oovv,X_vovv,X_ooov,W_abc,W_cba,W_bca,W_cab,W_bac,W_acb)
call form_v_abc(nO,nV,a,b,c,t1,X_oovv,W_abc,V_abc,W_cba,V_cba,W_bca,V_bca,W_cab,V_cab,W_bac,V_bac,W_acb,V_acb)
delta_abc = f_v(a) + f_v(b) + f_v(c)
e = 0.d0
do k = 1, nO
do j = 1, nO
do i = 1, nO
delta = 1.d0 / (f_o(i) + f_o(j) + f_o(k) - delta_abc)
e = e + delta * ( &
(4d0 * (W_abc(i,j,k) - W_cba(i,j,k)) + &
W_bca(i,j,k) - W_bac(i,j,k) + &
W_cab(i,j,k) - W_acb(i,j,k) ) * (V_abc(i,j,k) - V_cba(i,j,k)) +&
(4d0 * (W_acb(i,j,k) - W_bca(i,j,k)) + &
W_cba(i,j,k) - W_cab(i,j,k) + &
W_bac(i,j,k) - W_abc(i,j,k) ) * (V_acb(i,j,k) - V_bca(i,j,k)) +&
(4d0 * (W_bac(i,j,k) - W_cab(i,j,k)) + &
W_acb(i,j,k) - W_abc(i,j,k) + &
W_cba(i,j,k) - W_bca(i,j,k) ) * (V_bac(i,j,k) - V_cab(i,j,k)) )
enddo
enddo
enddo
deallocate(W_abc, W_cab, W_bca, W_bac, W_cba, W_acb, &
V_abc, V_cab, V_bca, V_bac, V_cba, V_acb )
end
double precision function ccsd_t_task_aba(a,b,nO,nV,t1,T_oovv,T_voov,&
X_ooov,X_oovv,X_vovv,f_o,f_v) result(e)
implicit none
integer, intent(in) :: nO,nV,a,b
double precision, intent(in) :: t1(nO,nV), f_o(nO), f_v(nV)
double precision, intent(in) :: X_oovv(nO,nO,nV,nV)
double precision, intent(in) :: T_voov(nV,nO,nO,nV), T_oovv(nO,nO,nV,nV)
double precision, intent(in) :: X_vovv(nV,nO,nV,nV), X_ooov(nO,nO,nO,nV)
double precision :: delta, delta_abc
integer :: i,j,k
double precision, allocatable :: W_abc(:,:,:), W_cab(:,:,:), W_bca(:,:,:)
double precision, allocatable :: W_bac(:,:,:), W_cba(:,:,:), W_acb(:,:,:)
double precision, allocatable :: V_abc(:,:,:), V_cab(:,:,:), V_bca(:,:,:)
double precision, allocatable :: V_bac(:,:,:), V_cba(:,:,:), V_acb(:,:,:)
allocate( W_abc(nO,nO,nO), W_cab(nO,nO,nO), W_bca(nO,nO,nO), &
W_bac(nO,nO,nO), W_cba(nO,nO,nO), W_acb(nO,nO,nO), &
V_abc(nO,nO,nO), V_cab(nO,nO,nO), V_bca(nO,nO,nO), &
V_bac(nO,nO,nO), V_cba(nO,nO,nO), V_acb(nO,nO,nO) )
call form_w_abc(nO,nV,a,b,a,T_voov,T_oovv,X_vovv,X_ooov,W_abc,W_cba,W_bca,W_cab,W_bac,W_acb)
call form_v_abc(nO,nV,a,b,a,t1,X_oovv,W_abc,V_abc,W_cba,V_cba,W_bca,V_bca,W_cab,V_cab,W_bac,V_bac,W_acb,V_acb)
delta_abc = f_v(a) + f_v(b) + f_v(a)
e = 0.d0
do k = 1, nO
do j = 1, nO
do i = 1, nO
delta = 1.d0 / (f_o(i) + f_o(j) + f_o(k) - delta_abc)
e = e + delta * ( &
(4d0 * W_abc(i,j,k) + W_bca(i,j,k) + W_cab(i,j,k)) * (V_abc(i,j,k) - V_cba(i,j,k)) + &
(4d0 * W_acb(i,j,k) + W_cba(i,j,k) + W_bac(i,j,k)) * (V_acb(i,j,k) - V_bca(i,j,k)) + &
(4d0 * W_bac(i,j,k) + W_acb(i,j,k) + W_cba(i,j,k)) * (V_bac(i,j,k) - V_cab(i,j,k)) )
enddo
enddo
enddo
deallocate(W_abc, W_cab, W_bca, W_bac, W_cba, W_acb, &
V_abc, V_cab, V_bca, V_bac, V_cba, V_acb )
end
subroutine form_w_abc(nO,nV,a,b,c,T_voov,T_oovv,X_vovv,X_ooov,W_abc,W_cba,W_bca,W_cab,W_bac,W_acb)
implicit none
integer, intent(in) :: nO,nV,a,b,c
!double precision, intent(in) :: t2(nO,nO,nV,nV)
double precision, intent(in) :: T_vvoo(nV,nV,nO,nO), T_ovvo(nO,nV,nV,nO)
double precision, intent(in) :: X_vvvo(nV,nV,nV,nO), X_ovoo(nO,nV,nO,nO)
double precision, intent(in) :: T_voov(nV,nO,nO,nV), T_oovv(nO,nO,nV,nV)
double precision, intent(in) :: X_vovv(nV,nO,nV,nV), X_ooov(nO,nO,nO,nV)
double precision, intent(out) :: W_abc(nO,nO,nO)
double precision, intent(out) :: W_cba(nO,nO,nO)
double precision, intent(out) :: W_bca(nO,nO,nO)
double precision, intent(out) :: W_cab(nO,nO,nO)
double precision, intent(out) :: W_bac(nO,nO,nO)
double precision, intent(out) :: W_acb(nO,nO,nO)
integer :: l,i,j,k,d
double precision, allocatable, dimension(:,:,:,:) :: W_ikj
double precision, allocatable :: X(:,:,:,:)
allocate(W_ikj(nO,nO,nO,6))
allocate(X(nV,nO,nO,3))
!$OMP PARALLEL &
!$OMP SHARED(nO,nV,a,b,c,T_vvoo,T_ovvo,X_vvvo,X_ovoo,W_abc) &
!$OMP PRIVATE(i,j,k,d,l) &
!$OMP DEFAULT(NONE)
!$OMP DO collapse(3)
do k = 1, nO
do j = 1, nO
do i = 1, nO
W_abc(i,j,k) = 0.d0
do d = 1, nV
W_abc(i,j,k) = W_abc(i,j,k) &
+ X_vvvo(d,b,a,i) * T_vvoo(d,c,k,j) &
+ X_vvvo(d,c,a,i) * T_vvoo(d,b,j,k) &
+ X_vvvo(d,a,c,k) * T_vvoo(d,b,j,i) &
+ X_vvvo(d,b,c,k) * T_vvoo(d,a,i,j) &
+ X_vvvo(d,c,b,j) * T_vvoo(d,a,i,k) &
+ X_vvvo(d,a,b,j) * T_vvoo(d,c,k,i)
enddo
do l = 1, nO
W_abc(i,j,k) = W_abc(i,j,k) &
- T_ovvo(l,a,b,i) * X_ovoo(l,c,j,k) &
- T_ovvo(l,a,c,i) * X_ovoo(l,b,k,j) & ! bc kj
- T_ovvo(l,c,a,k) * X_ovoo(l,b,i,j) & ! prev ac ik
- T_ovvo(l,c,b,k) * X_ovoo(l,a,j,i) & ! prev ab ij
- T_ovvo(l,b,c,j) * X_ovoo(l,a,k,i) & ! prev bc kj
- T_ovvo(l,b,a,j) * X_ovoo(l,c,i,k) ! prev ac ik
enddo
do k=1,nO
do i=1,nO
do d=1,nV
X(d,i,k,1) = T_voov(d,k,i,a)
X(d,i,k,2) = T_voov(d,k,i,b)
X(d,i,k,3) = T_voov(d,k,i,c)
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
! X_vovv(d,i,c,a) * T_voov(d,j,k,b) : i jk
call dgemm('T','N', nO, nO*nO, nV, 1.d0, X_vovv(1,1,c,a), nV, T_voov(1,1,1,b), nV, 0.d0, W_abc, nO)
call dgemm('T','N', nO, nO*nO, nV, 1.d0, X_vovv(1,1,c,b), nV, T_voov(1,1,1,a), nV, 0.d0, W_bac, nO)
call dgemm('T','N', nO, nO*nO, nV, 1.d0, X_vovv(1,1,a,c), nV, T_voov(1,1,1,b), nV, 0.d0, W_cba, nO)
call dgemm('T','N', nO, nO*nO, nV, 1.d0, X_vovv(1,1,a,b), nV, T_voov(1,1,1,c), nV, 0.d0, W_bca, nO)
call dgemm('T','N', nO, nO*nO, nV, 1.d0, X_vovv(1,1,b,c), nV, T_voov(1,1,1,a), nV, 0.d0, W_cab, nO)
call dgemm('T','N', nO, nO*nO, nV, 1.d0, X_vovv(1,1,b,a), nV, T_voov(1,1,1,c), nV, 0.d0, W_acb, nO)
! T_voov(d,i,j,a) * X_vovv(d,k,b,c) : ij k
call dgemm('T','N', nO*nO, nO, nV, 1.d0, T_voov(1,1,1,a), nV, X_vovv(1,1,b,c), nV, 1.d0, W_abc, nO*nO)
call dgemm('T','N', nO*nO, nO, nV, 1.d0, T_voov(1,1,1,b), nV, X_vovv(1,1,a,c), nV, 1.d0, W_bac, nO*nO)
call dgemm('T','N', nO*nO, nO, nV, 1.d0, T_voov(1,1,1,c), nV, X_vovv(1,1,b,a), nV, 1.d0, W_cba, nO*nO)
call dgemm('T','N', nO*nO, nO, nV, 1.d0, T_voov(1,1,1,b), nV, X_vovv(1,1,c,a), nV, 1.d0, W_bca, nO*nO)
call dgemm('T','N', nO*nO, nO, nV, 1.d0, T_voov(1,1,1,c), nV, X_vovv(1,1,a,b), nV, 1.d0, W_cab, nO*nO)
call dgemm('T','N', nO*nO, nO, nV, 1.d0, T_voov(1,1,1,a), nV, X_vovv(1,1,c,b), nV, 1.d0, W_acb, nO*nO)
! X_vovv(d,k,a,c) * T_voov(d,j,i,b) : k ji
call dgemm('T','N', nO*nO, nO, nV, 1.d0, X(1,1,1,2), nV, X_vovv(1,1,a,c), nV, 1.d0, W_abc, nO*nO)
call dgemm('T','N', nO*nO, nO, nV, 1.d0, X(1,1,1,1), nV, X_vovv(1,1,b,c), nV, 1.d0, W_bac, nO*nO)
call dgemm('T','N', nO*nO, nO, nV, 1.d0, X(1,1,1,2), nV, X_vovv(1,1,c,a), nV, 1.d0, W_cba, nO*nO)
call dgemm('T','N', nO*nO, nO, nV, 1.d0, X(1,1,1,3), nV, X_vovv(1,1,b,a), nV, 1.d0, W_bca, nO*nO)
call dgemm('T','N', nO*nO, nO, nV, 1.d0, X(1,1,1,1), nV, X_vovv(1,1,c,b), nV, 1.d0, W_cab, nO*nO)
call dgemm('T','N', nO*nO, nO, nV, 1.d0, X(1,1,1,3), nV, X_vovv(1,1,a,b), nV, 1.d0, W_acb, nO*nO)
! X_vovv(d,i,b,a) * T_voov(d,k,j,c) : i kj
call dgemm('T','N', nO, nO*nO, nV, 1.d0, X_vovv(1,1,b,a), nV, X(1,1,1,3), nV, 1.d0, W_abc, nO)
call dgemm('T','N', nO, nO*nO, nV, 1.d0, X_vovv(1,1,a,b), nV, X(1,1,1,3), nV, 1.d0, W_bac, nO)
call dgemm('T','N', nO, nO*nO, nV, 1.d0, X_vovv(1,1,b,c), nV, X(1,1,1,1), nV, 1.d0, W_cba, nO)
call dgemm('T','N', nO, nO*nO, nV, 1.d0, X_vovv(1,1,c,b), nV, X(1,1,1,1), nV, 1.d0, W_bca, nO)
call dgemm('T','N', nO, nO*nO, nV, 1.d0, X_vovv(1,1,a,c), nV, X(1,1,1,2), nV, 1.d0, W_cab, nO)
call dgemm('T','N', nO, nO*nO, nV, 1.d0, X_vovv(1,1,c,a), nV, X(1,1,1,2), nV, 1.d0, W_acb, nO)
! T_voov(d,k,i,c) * X_vovv(d,j,a,b) : ki j
call dgemm('T','N', nO*nO, nO, nV, 1.d0, X(1,1,1,3), nV, X_vovv(1,1,a,b), nV, 0.d0, W_ikj(1,1,1,1), nO*nO)
call dgemm('T','N', nO*nO, nO, nV, 1.d0, X(1,1,1,3), nV, X_vovv(1,1,b,a), nV, 0.d0, W_ikj(1,1,1,2), nO*nO)
call dgemm('T','N', nO*nO, nO, nV, 1.d0, X(1,1,1,1), nV, X_vovv(1,1,c,b), nV, 0.d0, W_ikj(1,1,1,3), nO*nO)
call dgemm('T','N', nO*nO, nO, nV, 1.d0, X(1,1,1,1), nV, X_vovv(1,1,b,c), nV, 0.d0, W_ikj(1,1,1,4), nO*nO)
call dgemm('T','N', nO*nO, nO, nV, 1.d0, X(1,1,1,2), nV, X_vovv(1,1,c,a), nV, 0.d0, W_ikj(1,1,1,5), nO*nO)
call dgemm('T','N', nO*nO, nO, nV, 1.d0, X(1,1,1,2), nV, X_vovv(1,1,a,c), nV, 0.d0, W_ikj(1,1,1,6), nO*nO)
! T_voov(d,i,k,a) * X_vovv(d,j,c,b) : ik j
call dgemm('T','N', nO*nO, nO, nV, 1.d0, T_voov(1,1,1,a), nV, X_vovv(1,1,c,b), nV, 1.d0, W_ikj(1,1,1,1), nO*nO)
call dgemm('T','N', nO*nO, nO, nV, 1.d0, T_voov(1,1,1,b), nV, X_vovv(1,1,c,a), nV, 1.d0, W_ikj(1,1,1,2), nO*nO)
call dgemm('T','N', nO*nO, nO, nV, 1.d0, T_voov(1,1,1,c), nV, X_vovv(1,1,a,b), nV, 1.d0, W_ikj(1,1,1,3), nO*nO)
call dgemm('T','N', nO*nO, nO, nV, 1.d0, T_voov(1,1,1,b), nV, X_vovv(1,1,a,c), nV, 1.d0, W_ikj(1,1,1,4), nO*nO)
call dgemm('T','N', nO*nO, nO, nV, 1.d0, T_voov(1,1,1,c), nV, X_vovv(1,1,b,a), nV, 1.d0, W_ikj(1,1,1,5), nO*nO)
call dgemm('T','N', nO*nO, nO, nV, 1.d0, T_voov(1,1,1,a), nV, X_vovv(1,1,b,c), nV, 1.d0, W_ikj(1,1,1,6), nO*nO)
deallocate(X)
allocate(X(nO,nO,nO,3))
do k=1,nO
do j=1,nO
do l=1,nO
X(l,j,k,1) = X_ooov(l,k,j,a)
X(l,j,k,2) = X_ooov(l,k,j,b)
X(l,j,k,3) = X_ooov(l,k,j,c)
enddo
enddo
enddo
! - T_oovv(l,i,a,b) * X_ooov(l,j,k,c) : i jk
call dgemm('T','N', nO, nO*nO, nO, -1.d0, T_oovv(1,1,a,b), nO, X_ooov(1,1,1,c), nO, 1.d0, W_abc, nO)
call dgemm('T','N', nO, nO*nO, nO, -1.d0, T_oovv(1,1,b,a), nO, X_ooov(1,1,1,c), nO, 1.d0, W_bac, nO)
call dgemm('T','N', nO, nO*nO, nO, -1.d0, T_oovv(1,1,c,b), nO, X_ooov(1,1,1,a), nO, 1.d0, W_cba, nO)
call dgemm('T','N', nO, nO*nO, nO, -1.d0, T_oovv(1,1,b,c), nO, X_ooov(1,1,1,a), nO, 1.d0, W_bca, nO)
call dgemm('T','N', nO, nO*nO, nO, -1.d0, T_oovv(1,1,c,a), nO, X_ooov(1,1,1,b), nO, 1.d0, W_cab, nO)
call dgemm('T','N', nO, nO*nO, nO, -1.d0, T_oovv(1,1,a,c), nO, X_ooov(1,1,1,b), nO, 1.d0, W_acb, nO)
! - T_oovv(l,i,a,c) * X_ooov(l,k,j,b) : i kj
call dgemm('T','N', nO, nO*nO, nO, -1.d0, T_oovv(1,1,a,c), nO, X(1,1,1,2), nO, 1.d0, W_abc, nO)
call dgemm('T','N', nO, nO*nO, nO, -1.d0, T_oovv(1,1,b,c), nO, X(1,1,1,1), nO, 1.d0, W_bac, nO)
call dgemm('T','N', nO, nO*nO, nO, -1.d0, T_oovv(1,1,c,a), nO, X(1,1,1,2), nO, 1.d0, W_cba, nO)
call dgemm('T','N', nO, nO*nO, nO, -1.d0, T_oovv(1,1,b,a), nO, X(1,1,1,3), nO, 1.d0, W_bca, nO)
call dgemm('T','N', nO, nO*nO, nO, -1.d0, T_oovv(1,1,c,b), nO, X(1,1,1,1), nO, 1.d0, W_cab, nO)
call dgemm('T','N', nO, nO*nO, nO, -1.d0, T_oovv(1,1,a,b), nO, X(1,1,1,3), nO, 1.d0, W_acb, nO)
! - X_ooov(l,i,j,b) * T_oovv(l,k,c,a) : ij k
call dgemm('T','N', nO*nO, nO, nO, -1.d0, X_ooov(1,1,1,b), nO, T_oovv(1,1,c,a), nO, 1.d0, W_abc, nO*nO)
call dgemm('T','N', nO*nO, nO, nO, -1.d0, X_ooov(1,1,1,a), nO, T_oovv(1,1,c,b), nO, 1.d0, W_bac, nO*nO)
call dgemm('T','N', nO*nO, nO, nO, -1.d0, X_ooov(1,1,1,b), nO, T_oovv(1,1,a,c), nO, 1.d0, W_cba, nO*nO)
call dgemm('T','N', nO*nO, nO, nO, -1.d0, X_ooov(1,1,1,c), nO, T_oovv(1,1,a,b), nO, 1.d0, W_bca, nO*nO)
call dgemm('T','N', nO*nO, nO, nO, -1.d0, X_ooov(1,1,1,a), nO, T_oovv(1,1,b,c), nO, 1.d0, W_cab, nO*nO)
call dgemm('T','N', nO*nO, nO, nO, -1.d0, X_ooov(1,1,1,c), nO, T_oovv(1,1,b,a), nO, 1.d0, W_acb, nO*nO)
! - X_ooov(l,j,i,a) * T_oovv(l,k,c,b) : ji k
call dgemm('T','N', nO*nO, nO, nO, -1.d0, X(1,1,1,1), nO, T_oovv(1,1,c,b), nO, 1.d0, W_abc, nO*nO)
call dgemm('T','N', nO*nO, nO, nO, -1.d0, X(1,1,1,2), nO, T_oovv(1,1,c,a), nO, 1.d0, W_bac, nO*nO)
call dgemm('T','N', nO*nO, nO, nO, -1.d0, X(1,1,1,3), nO, T_oovv(1,1,a,b), nO, 1.d0, W_cba, nO*nO)
call dgemm('T','N', nO*nO, nO, nO, -1.d0, X(1,1,1,2), nO, T_oovv(1,1,a,c), nO, 1.d0, W_bca, nO*nO)
call dgemm('T','N', nO*nO, nO, nO, -1.d0, X(1,1,1,3), nO, T_oovv(1,1,b,a), nO, 1.d0, W_cab, nO*nO)
call dgemm('T','N', nO*nO, nO, nO, -1.d0, X(1,1,1,1), nO, T_oovv(1,1,b,c), nO, 1.d0, W_acb, nO*nO)
! - X_ooov(l,k,i,a) * T_oovv(l,j,b,c) : ki j
call dgemm('T','N', nO*nO, nO, nO, -1.d0, X(1,1,1,1), nO, T_oovv(1,1,b,c), nO, 1.d0, W_ikj(1,1,1,1), nO*nO)
call dgemm('T','N', nO*nO, nO, nO, -1.d0, X(1,1,1,2), nO, T_oovv(1,1,a,c), nO, 1.d0, W_ikj(1,1,1,2), nO*nO)
call dgemm('T','N', nO*nO, nO, nO, -1.d0, X(1,1,1,3), nO, T_oovv(1,1,b,a), nO, 1.d0, W_ikj(1,1,1,3), nO*nO)
call dgemm('T','N', nO*nO, nO, nO, -1.d0, X(1,1,1,2), nO, T_oovv(1,1,c,a), nO, 1.d0, W_ikj(1,1,1,4), nO*nO)
call dgemm('T','N', nO*nO, nO, nO, -1.d0, X(1,1,1,3), nO, T_oovv(1,1,a,b), nO, 1.d0, W_ikj(1,1,1,5), nO*nO)
call dgemm('T','N', nO*nO, nO, nO, -1.d0, X(1,1,1,1), nO, T_oovv(1,1,c,b), nO, 1.d0, W_ikj(1,1,1,6), nO*nO)
! - X_ooov(l,i,k,c) * T_oovv(l,j,b,a) : ik j
call dgemm('T','N', nO*nO, nO, nO, -1.d0, X_ooov(1,1,1,c), nO, T_oovv(1,1,b,a), nO, 1.d0, W_ikj(1,1,1,1), nO*nO)
call dgemm('T','N', nO*nO, nO, nO, -1.d0, X_ooov(1,1,1,c), nO, T_oovv(1,1,a,b), nO, 1.d0, W_ikj(1,1,1,2), nO*nO)
call dgemm('T','N', nO*nO, nO, nO, -1.d0, X_ooov(1,1,1,a), nO, T_oovv(1,1,b,c), nO, 1.d0, W_ikj(1,1,1,3), nO*nO)
call dgemm('T','N', nO*nO, nO, nO, -1.d0, X_ooov(1,1,1,a), nO, T_oovv(1,1,c,b), nO, 1.d0, W_ikj(1,1,1,4), nO*nO)
call dgemm('T','N', nO*nO, nO, nO, -1.d0, X_ooov(1,1,1,b), nO, T_oovv(1,1,a,c), nO, 1.d0, W_ikj(1,1,1,5), nO*nO)
call dgemm('T','N', nO*nO, nO, nO, -1.d0, X_ooov(1,1,1,b), nO, T_oovv(1,1,c,a), nO, 1.d0, W_ikj(1,1,1,6), nO*nO)
do k=1,nO
do j=1,nO
do i=1,nO
W_abc(i,j,k) = W_abc(i,j,k) + W_ikj(i,k,j,1)
W_bac(i,j,k) = W_bac(i,j,k) + W_ikj(i,k,j,2)
W_cba(i,j,k) = W_cba(i,j,k) + W_ikj(i,k,j,3)
W_bca(i,j,k) = W_bca(i,j,k) + W_ikj(i,k,j,4)
W_cab(i,j,k) = W_cab(i,j,k) + W_ikj(i,k,j,5)
W_acb(i,j,k) = W_acb(i,j,k) + W_ikj(i,k,j,6)
enddo
enddo
enddo
deallocate(X,W_ikj)
end
! V_abc
subroutine form_v_abc(nO,nV,a,b,c,T_vo,X_vvoo,W,V)
subroutine form_v_abc(nO,nV,a,b,c,T_ov,X_oovv,W_abc,V_abc,W_cba,V_cba,W_bca,V_bca,W_cab,V_cab,W_bac,V_bac,W_acb,V_acb)
implicit none
integer, intent(in) :: nO,nV,a,b,c
!double precision, intent(in) :: t1(nO,nV)
double precision, intent(in) :: T_vo(nV,nO)
double precision, intent(in) :: X_vvoo(nV,nV,nO,nO)
double precision, intent(in) :: W(nO,nO,nO)
double precision, intent(out) :: V(nO,nO,nO)
double precision, intent(in) :: T_ov(nO,nV)
double precision, intent(in) :: X_oovv(nO,nO,nV,nV)
double precision, intent(in) :: W_abc(nO,nO,nO), W_cab(nO,nO,nO), W_bca(nO,nO,nO)
double precision, intent(in) :: W_bac(nO,nO,nO), W_cba(nO,nO,nO), W_acb(nO,nO,nO)
double precision, intent(out) :: V_abc(nO,nO,nO), V_cab(nO,nO,nO), V_bca(nO,nO,nO)
double precision, intent(out) :: V_bac(nO,nO,nO), V_cba(nO,nO,nO), V_acb(nO,nO,nO)
integer :: i,j,k
!$OMP PARALLEL &
!$OMP SHARED(nO,nV,a,b,c,T_vo,X_vvoo,W,V) &
!$OMP PRIVATE(i,j,k) &
!$OMP DEFAULT(NONE)
!$OMP DO collapse(2)
do k = 1, nO
do j = 1, nO
do i = 1, nO
!V(i,j,k,a,b,c) = V(i,j,k,a,b,c) + W(i,j,k,a,b,c) &
V(i,j,k) = W(i,j,k) &
+ X_vvoo(b,c,k,j) * T_vo(a,i) &
+ X_vvoo(a,c,k,i) * T_vo(b,j) &
+ X_vvoo(a,b,j,i) * T_vo(c,k)
V_abc(i,j,k) = W_abc(i,j,k) &
+ X_oovv(j,k,b,c) * T_ov(i,a) &
+ X_oovv(i,k,a,c) * T_ov(j,b) &
+ X_oovv(i,j,a,b) * T_ov(k,c)
V_cba(i,j,k) = W_cba(i,j,k) &
+ X_oovv(j,k,b,a) * T_ov(i,c) &
+ X_oovv(i,k,c,a) * T_ov(j,b) &
+ X_oovv(i,j,c,b) * T_ov(k,a)
V_bca(i,j,k) = W_bca(i,j,k) &
+ X_oovv(j,k,c,a) * T_ov(i,b) &
+ X_oovv(i,k,b,a) * T_ov(j,c) &
+ X_oovv(i,j,b,c) * T_ov(k,a)
V_cab(i,j,k) = W_cab(i,j,k) &
+ X_oovv(j,k,a,b) * T_ov(i,c) &
+ X_oovv(i,k,c,b) * T_ov(j,a) &
+ X_oovv(i,j,c,a) * T_ov(k,b)
V_bac(i,j,k) = W_bac(i,j,k) &
+ X_oovv(j,k,a,c) * T_ov(i,b) &
+ X_oovv(i,k,b,c) * T_ov(j,a) &
+ X_oovv(i,j,b,a) * T_ov(k,c)
V_acb(i,j,k) = W_acb(i,j,k) &
+ X_oovv(j,k,c,b) * T_ov(i,a) &
+ X_oovv(i,k,a,b) * T_ov(j,c) &
+ X_oovv(i,j,a,c) * T_ov(k,b)
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
end

363
src/ccsd/ccsd_t_space_orb_stoch.irp.f Normal file
View File

@ -0,0 +1,363 @@
! Main
subroutine ccsd_par_t_space_stoch(nO,nV,t1,t2,f_o,f_v,v_vvvo,v_vvoo,v_vooo,energy)
implicit none
integer, intent(in) :: nO,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) :: 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, allocatable :: X_vovv(:,:,:,:), X_ooov(:,:,:,:), X_oovv(:,:,:,:)
double precision, allocatable :: T_voov(:,:,:,:), T_oovv(:,:,:,:)
integer :: i,j,k,l,a,b,c,d
double precision :: e,ta,tb,eccsd
eccsd = energy
call set_multiple_levels_omp(.False.)
allocate(X_vovv(nV,nO,nV,nV), X_ooov(nO,nO,nO,nV), X_oovv(nO,nO,nV,nV))
allocate(T_voov(nV,nO,nO,nV),T_oovv(nO,nO,nV,nV))
!$OMP PARALLEL &
!$OMP SHARED(nO,nV,T_voov,T_oovv,X_vovv,X_ooov,X_oovv, &
!$OMP t1,t2,v_vvvo,v_vooo,v_vvoo) &
!$OMP PRIVATE(a,b,c,d,i,j,k,l) &
!$OMP DEFAULT(NONE)
!v_vvvo(b,a,d,i) * t2(k,j,c,d) &
!X_vovv(d,i,b,a,i) * T_voov(d,j,c,k)
!$OMP DO
do a = 1, nV
do b = 1, nV
do i = 1, nO
do d = 1, nV
X_vovv(d,i,b,a) = v_vvvo(b,a,d,i)
enddo
enddo
enddo
enddo
!$OMP END DO nowait
!$OMP DO
do c = 1, nV
do j = 1, nO
do k = 1, nO
do d = 1, nV
T_voov(d,k,j,c) = t2(k,j,c,d)
enddo
enddo
enddo
enddo
!$OMP END DO nowait
!v_vooo(c,j,k,l) * t2(i,l,a,b) &
!X_ooov(l,j,k,c) * T_oovv(l,i,a,b) &
!$OMP DO
do c = 1, nV
do k = 1, nO
do j = 1, nO
do l = 1, nO
X_ooov(l,j,k,c) = v_vooo(c,j,k,l)
enddo
enddo
enddo
enddo
!$OMP END DO nowait
!$OMP DO
do b = 1, nV
do a = 1, nV
do i = 1, nO
do l = 1, nO
T_oovv(l,i,a,b) = t2(i,l,a,b)
enddo
enddo
enddo
enddo
!$OMP END DO nowait
!X_oovv(j,k,b,c) * T1_vo(a,i) &
!$OMP DO
do c = 1, nV
do b = 1, nV
do k = 1, nO
do j = 1, nO
X_oovv(j,k,b,c) = v_vvoo(b,c,j,k)
enddo
enddo
enddo
enddo
!$OMP END DO nowait
!$OMP END PARALLEL
double precision, external :: ccsd_t_task_aba
double precision, external :: ccsd_t_task_abc
! logical, external :: omp_test_lock
double precision, allocatable :: memo(:), Pabc(:), waccu(:)
integer*8, allocatable :: sampled(:)
! integer(omp_lock_kind), allocatable :: lock(:)
integer*2 , allocatable :: abc(:,:)
integer*8 :: Nabc, i8
integer*8, allocatable :: iorder(:)
double precision :: eocc
double precision :: norm
integer :: kiter, isample
! Prepare table of triplets (a,b,c)
Nabc = (int(nV,8) * int(nV+1,8) * int(nV+2,8))/6_8 - nV
allocate (memo(Nabc), sampled(Nabc), Pabc(Nabc), waccu(Nabc))
allocate (abc(4,Nabc), iorder(Nabc)) !, lock(Nabc))
! eocc = 3.d0/dble(nO) * sum(f_o(1:nO))
Nabc = 0_8
do a = 1, nV
do b = a+1, nV
do c = b+1, nV
Nabc = Nabc + 1_8
Pabc(Nabc) = -1.d0/(f_v(a) + f_v(b) + f_v(c))
abc(1,Nabc) = a
abc(2,Nabc) = b
abc(3,Nabc) = c
enddo
Nabc = Nabc + 1_8
abc(1,Nabc) = a
abc(2,Nabc) = b
abc(3,Nabc) = a
Pabc(Nabc) = -1.d0/(2.d0*f_v(a) + f_v(b))
Nabc = Nabc + 1_8
abc(1,Nabc) = b
abc(2,Nabc) = a
abc(3,Nabc) = b
Pabc(Nabc) = -1.d0/(f_v(a) + 2.d0*f_v(b))
enddo
enddo
do i8=1,Nabc
iorder(i8) = i8
enddo
! Sort triplets in decreasing Pabc
call dsort_big(Pabc, iorder, Nabc)
! Normalize
norm = 0.d0
do i8=Nabc,1,-1
norm = norm + Pabc(i8)
enddo
norm = 1.d0/norm
do i8=1,Nabc
Pabc(i8) = Pabc(i8) * norm
enddo
call i8set_order_big(abc, iorder, Nabc)
! Cumulative distribution for sampling
waccu(Nabc) = 0.d0
do i8=Nabc-1,1,-1
waccu(i8) = waccu(i8+1) - Pabc(i8+1)
enddo
waccu(:) = waccu(:) + 1.d0
logical :: converged, do_comp
double precision :: eta, variance, error, sample
double precision :: t00, t01
integer*8 :: ieta, Ncomputed
integer*8, external :: binary_search
integer :: nbuckets
nbuckets = 100
double precision, allocatable :: wsum(:)
allocate(wsum(nbuckets))
converged = .False.
Ncomputed = 0_8
energy = 0.d0
variance = 0.d0
memo(:) = 0.d0
sampled(:) = -1_8
integer*8 :: ileft, iright, imin
ileft = 1_8
iright = Nabc
integer*8, allocatable :: bounds(:,:)
allocate (bounds(2,nbuckets))
do isample=1,nbuckets
eta = 1.d0/dble(nbuckets) * dble(isample)
ieta = binary_search(waccu,eta,Nabc)
bounds(1,isample) = ileft
bounds(2,isample) = ieta
ileft = ieta+1
wsum(isample) = sum( Pabc(bounds(1,isample):bounds(2,isample) ) )
enddo
Pabc(:) = 1.d0/Pabc(:)
print '(A)', ''
print '(A)', ' +----------------------+--------------+----------+'
print '(A)', ' | E(CCSD(T)) | Error | % |'
print '(A)', ' +----------------------+--------------+----------+'
call wall_time(t00)
imin = 1_8
!$OMP PARALLEL &
!$OMP PRIVATE(ieta,eta,a,b,c,kiter,isample) &
!$OMP DEFAULT(SHARED)
do kiter=1,Nabc
!$OMP MASTER
do while ((imin <= Nabc).and.(sampled(imin)>-1_8))
imin = imin+1
enddo
! Deterministic part
if (imin < Nabc) then
ieta=imin
sampled(ieta) = 0_8
a = abc(1,ieta)
b = abc(2,ieta)
c = abc(3,ieta)
Ncomputed += 1_8
!$OMP TASK DEFAULT(SHARED) FIRSTPRIVATE(a,b,c,ieta)
if (a/=c) then
memo(ieta) = ccsd_t_task_abc(a,b,c,nO,nV,t1,T_oovv,T_voov, &
X_ooov,X_oovv,X_vovv,f_o,f_v) / 3.d0
else
memo(ieta) = ccsd_t_task_aba(a,b,nO,nV,t1,T_oovv,T_voov, &
X_ooov,X_oovv,X_vovv,f_o,f_v) / 3.d0
endif
!$OMP END TASK
endif
! Stochastic part
call random_number(eta)
do isample=1,nbuckets
if (imin >= bounds(2,isample)) then
cycle
endif
ieta = binary_search(waccu,(eta + dble(isample-1))/dble(nbuckets),Nabc)
if (sampled(ieta) == -1_8) then
sampled(ieta) = 0_8
a = abc(1,ieta)
b = abc(2,ieta)
c = abc(3,ieta)
Ncomputed += 1_8
!$OMP TASK DEFAULT(SHARED) FIRSTPRIVATE(a,b,c,ieta)
if (a/=c) then
memo(ieta) = ccsd_t_task_abc(a,b,c,nO,nV,t1,T_oovv,T_voov, &
X_ooov,X_oovv,X_vovv,f_o,f_v) / 3.d0
else
memo(ieta) = ccsd_t_task_aba(a,b,nO,nV,t1,T_oovv,T_voov, &
X_ooov,X_oovv,X_vovv,f_o,f_v) / 3.d0
endif
!$OMP END TASK
endif
sampled(ieta) = sampled(ieta)+1_8
enddo
call wall_time(t01)
if ((t01-t00 > 1.0d0).or.(imin >= Nabc)) then
t00 = t01
!$OMP TASKWAIT
double precision :: ET, ET2
double precision :: energy_stoch, energy_det
double precision :: scale
double precision :: w
double precision :: tmp
energy_stoch = 0.d0
energy_det = 0.d0
norm = 0.d0
scale = 1.d0
ET = 0.d0
ET2 = 0.d0
do isample=1,nbuckets
if (imin >= bounds(2,isample)) then
energy_det = energy_det + sum(memo(bounds(1,isample):bounds(2,isample)))
scale = scale - wsum(isample)
else
exit
endif
enddo
do ieta=bounds(1,isample), Nabc
w = dble(max(sampled(ieta),0_8))
tmp = w * memo(ieta) * Pabc(ieta)
ET = ET + tmp
ET2 = ET2 + tmp * memo(ieta) * Pabc(ieta)
norm = norm + w
enddo
norm = norm/scale
if (norm > 0.d0) then
energy_stoch = ET / norm
variance = ET2 / norm - energy_stoch*energy_stoch
endif
energy = energy_det + energy_stoch
print '('' | '',F20.8, '' | '', E12.4,'' | '', F8.2,'' |'')', eccsd+energy, dsqrt(variance/(norm-1.d0)), 100.*real(Ncomputed)/real(Nabc)
endif
!$OMP END MASTER
if (imin >= Nabc) exit
enddo
!$OMP END PARALLEL
print '(A)', ' +----------------------+--------------+----------+'
print '(A)', ''
deallocate(X_vovv,X_ooov,T_voov,T_oovv)
end
integer*8 function binary_search(arr, key, size)
implicit none
BEGIN_DOC
! Searches the key in array arr(1:size) between l_in and r_in, and returns its index
END_DOC
integer*8 :: size, i, j, mid, l_in, r_in
double precision, dimension(size) :: arr(1:size)
double precision :: key
i = 1_8
j = size
do while (j >= i)
mid = i + (j - i) / 2
if (arr(mid) >= key) then
if (mid > 1 .and. arr(mid - 1) < key) then
binary_search = mid
return
end if
j = mid - 1
else if (arr(mid) < key) then
i = mid + 1
else
binary_search = mid + 1
return
end if
end do
binary_search = i
end function binary_search

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

@ -76,6 +76,8 @@ subroutine select_connected(i_generator,E0,pt2_data,b,subset,csubset)
double precision, allocatable :: fock_diag_tmp(:,:)
if (csubset == 0) return
allocate(fock_diag_tmp(2,mo_num+1))
call build_fock_tmp(fock_diag_tmp,psi_det_generators(1,1,i_generator),N_int)
@ -177,6 +179,7 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d
monoAdo = .true.
monoBdo = .true.
if (csubset == 0) return
do k=1,N_int
hole (k,1) = iand(psi_det_generators(k,1,i_generator), hole_mask(k,1))

5
src/cipsi_tc_bi_ortho/selection.irp.f
View File

@ -868,7 +868,6 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d
! <det|H(j)|psi_0> and transpose
! -------------------------------------------
! call htilde_mu_mat_bi_ortho_tot(det, det, N_int, Hii)
double precision :: hmono, htwoe, hthree
call diag_htilde_mu_mat_fock_bi_ortho(N_int, det, hmono, htwoe, hthree, hii)
do istate = 1,N_states
@ -878,8 +877,8 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d
psi_h_alpha = 0.d0
alpha_h_psi = 0.d0
do iii = 1, N_det_selectors
call htilde_mu_mat_bi_ortho_tot(psi_selectors(1,1,iii), det, N_int, i_h_alpha)
call htilde_mu_mat_bi_ortho_tot(det, psi_selectors(1,1,iii), N_int, alpha_h_i)
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 !!!'

19
src/cosgtos_ao_int/EZFIO.cfg
View File

@ -1,19 +0,0 @@
[ao_expoim_cosgtos]
type: double precision
doc: imag part for Exponents for each primitive of each cosGTOs |AO|
size: (ao_basis.ao_num,ao_basis.ao_prim_num_max)
interface: ezfio, provider
[use_cosgtos]
type: logical
doc: If true, use cosgtos for AO integrals
interface: ezfio,provider,ocaml
default: False
[ao_integrals_threshold]
type: Threshold
doc: If | (pq|rs) | < `ao_integrals_threshold` then (pq|rs) is zero
interface: ezfio,provider,ocaml
default: 1.e-15
ezfio_name: threshold_ao

2
src/cosgtos_ao_int/NEED
View File

@ -1,2 +0,0 @@
ezfio_files
ao_basis

4
src/cosgtos_ao_int/README.rst
View File

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

7
src/cosgtos_ao_int/cosgtos_ao_int.irp.f
View File

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

2
src/fci_tc_bi/diagonalize_ci.irp.f
View File

@ -39,7 +39,7 @@ subroutine diagonalize_CI_tc_bi_ortho(ndet, E_tc,norm,pt2_data,print_pt2)
write(*,'(A28,X,I10,X,100(F16.8,X))')'Ndet,E,E+PT2,E+RPT2,|PT2|=',ndet,E_tc ,E_tc + pt2_tmp/norm,E_tc + rpt2_tmp/norm,abs_pt2
print*,'*****'
endif
psi_energy(1:N_states) = eigval_right_tc_bi_orth(1:N_states)
psi_energy(1:N_states) = eigval_right_tc_bi_orth(1:N_states) - nuclear_repulsion
psi_s2(1:N_states) = s2_eigvec_tc_bi_orth(1:N_states)
E_tc = eigval_right_tc_bi_orth(1)

31
src/mo_two_e_ints/cholesky.irp.f
View File

@ -6,11 +6,42 @@ BEGIN_PROVIDER [ double precision, cholesky_mo, (mo_num, mo_num, cholesky_ao_num
integer :: k
call set_multiple_levels_omp(.False.)
print *, 'AO->MO Transformation of Cholesky vectors'
!$OMP PARALLEL DO PRIVATE(k)
do k=1,cholesky_ao_num
call ao_to_mo(cholesky_ao(1,1,k),ao_num,cholesky_mo(1,1,k),mo_num)
enddo
!$OMP END PARALLEL DO
print *, ''
END_PROVIDER
BEGIN_PROVIDER [ double precision, cholesky_mo_transp, (cholesky_ao_num, mo_num, mo_num) ]
implicit none
BEGIN_DOC
! Cholesky vectors in MO basis
END_DOC
integer :: i,j,k
double precision, allocatable :: buffer(:,:)
print *, 'AO->MO Transformation of Cholesky vectors .'
!$OMP PARALLEL PRIVATE(i,j,k,buffer)
allocate(buffer(mo_num,mo_num))
!$OMP DO SCHEDULE(static)
do k=1,cholesky_ao_num
call ao_to_mo(cholesky_ao(1,1,k),ao_num,buffer,mo_num)
do j=1,mo_num
do i=1,mo_num
cholesky_mo_transp(k,i,j) = buffer(i,j)
enddo
enddo
enddo
!$OMP END DO
deallocate(buffer)
!$OMP END PARALLEL
print *, ''
END_PROVIDER

48
src/mo_two_e_ints/integrals_3_index.irp.f
View File

@ -4,12 +4,54 @@
BEGIN_DOC
! big_array_coulomb_integrals(j,i,k) = <ij|kj> = (ik|jj)
!
! big_array_exchange_integrals(i,j,k) = <ij|jk> = (ij|kj)
! big_array_exchange_integrals(j,i,k) = <ij|jk> = (ij|kj)
END_DOC
integer :: i,j,k,l
integer :: i,j,k,l,a
double precision :: get_two_e_integral
double precision :: integral
if (do_ao_cholesky) then
double precision, allocatable :: buffer_jj(:,:), buffer(:,:,:)
allocate(buffer_jj(cholesky_ao_num,mo_num), buffer(mo_num,mo_num,mo_num))
do j=1,mo_num
buffer_jj(:,j) = cholesky_mo_transp(:,j,j)
enddo
call dgemm('T','N', mo_num*mo_num,mo_num,cholesky_ao_num, 1.d0, &
cholesky_mo_transp, cholesky_ao_num, &
buffer_jj, cholesky_ao_num, 0.d0, &
buffer, mo_num*mo_num)
do k = 1, mo_num
do i = 1, mo_num
do j = 1, mo_num
big_array_coulomb_integrals(j,i,k) = buffer(i,k,j)
enddo
enddo
enddo
deallocate(buffer_jj)
allocate(buffer_jj(mo_num,mo_num))
do j = 1, mo_num
call dgemm('T','N',mo_num,mo_num,cholesky_ao_num, 1.d0, &
cholesky_mo_transp(1,1,j), cholesky_ao_num, &
cholesky_mo_transp(1,1,j), cholesky_ao_num, 0.d0, &
buffer_jj, mo_num)
do k=1,mo_num
do i=1,mo_num
big_array_exchange_integrals(j,i,k) = buffer_jj(i,k)
enddo
enddo
enddo
deallocate(buffer_jj)
else
do k = 1, mo_num
do i = 1, mo_num
do j = 1, mo_num
@ -23,5 +65,7 @@
enddo
enddo
endif
END_PROVIDER

21
src/mo_two_e_ints/mo_bi_integrals.irp.f
View File

@ -1353,14 +1353,29 @@ END_PROVIDER
integer :: i,j
double precision :: get_two_e_integral
PROVIDE mo_two_e_integrals_in_map
mo_two_e_integrals_jj = 0.d0
mo_two_e_integrals_jj_exchange = 0.d0
if (do_ao_cholesky) then
do j=1,mo_num
do i=1,mo_num
!TODO: use dgemm
mo_two_e_integrals_jj(i,j) = sum(cholesky_mo_transp(:,i,i)*cholesky_mo_transp(:,j,j))
mo_two_e_integrals_jj_exchange(i,j) = sum(cholesky_mo_transp(:,i,j)*cholesky_mo_transp(:,j,i))
enddo
enddo
else
do j=1,mo_num
do i=1,mo_num
mo_two_e_integrals_jj(i,j) = get_two_e_integral(i,j,i,j,mo_integrals_map)
mo_two_e_integrals_jj_exchange(i,j) = get_two_e_integral(i,j,j,i,mo_integrals_map)
enddo
enddo
endif
do j=1,mo_num
do i=1,mo_num
mo_two_e_integrals_jj_anti(i,j) = mo_two_e_integrals_jj(i,j) - mo_two_e_integrals_jj_exchange(i,j)
enddo
enddo

167
src/non_h_ints_mu/jast_deriv.irp.f
View File

@ -187,6 +187,19 @@ end function j12_mu
subroutine grad1_j12_mu(r1, r2, grad)
BEGIN_DOC
! gradient of j(mu(r1,r2),r12) form of jastrow.
!
! if mu(r1,r2) = cst ---> j1b_type < 200 and
!
! d/dx1 j(mu,r12) = 0.5 * (1 - erf(mu *r12))/r12 * (x1 - x2)
!
! if mu(r1,r2) /= cst ---> 200 < j1b_type < 300 and
!
! d/dx1 j(mu(r1,r2),r12) = exp(-(mu(r1,r2)*r12)**2) /(2 *sqrt(pi) * mu(r1,r2)**2 ) d/dx1 mu(r1,r2)
!
! + 0.5 * (1 - erf(mu(r1,r2) *r12))/r12 * (x1 - x2)
END_DOC
include 'constants.include.F'
implicit none
@ -515,6 +528,9 @@ subroutine mu_r_val_and_grad(r1, r2, mu_val, mu_der)
double precision :: r(3)
double precision :: dm_a(1), dm_b(1), grad_dm_a(3,1), grad_dm_b(3,1)
double precision :: dm_tot, tmp1, tmp2, tmp3
double precision :: rho1, grad_rho1(3),rho2,rho_tot,inv_rho_tot
double precision :: f_rho1, f_rho2, d_drho_f_rho1
double precision :: d_dx1_f_rho1(3),d_dx_rho_f_rho(3),nume
if(j1b_type .eq. 200) then
@ -578,8 +594,84 @@ subroutine mu_r_val_and_grad(r1, r2, mu_val, mu_der)
mu_der(2) = tmp3 * (grad_dm_a(2,1) + grad_dm_b(2,1))
mu_der(3) = tmp3 * (grad_dm_a(3,1) + grad_dm_b(3,1))
else
elseif(j1b_type .eq. 202) then
! mu(r1,r2) = {rho(r1) f[rho(r1)] + rho(r2) f[rho(r2)]} / RHO
!
! RHO = rho(r1) + rho(r2)
!
! f[rho] = alpha rho^beta + mu0 exp(-rho)
!
! d/dx1 mu(r1,r2) = 1/RHO^2 * {RHO * d/dx1 (rho(r1) f[rho(r1)])
! - d/dx1 rho(r1) * [rho(r1) f[rho(r1)] + rho(r2) f[rho(r2)]] }
!
! d/dx1 f[rho(r1)] = [0.5 alpha / sqrt(rho(r1)) - mu0 exp(-rho(r1))] (d rho(r1) / dx1)
!
! d/dx1 (rho(r1) f[rho(r1)] = rho(r1) * d/dx1 f[rho(r1)] + f[rho(r1)] * d/dx1 rho(r1)
!!!!!!!!! rho1,rho2,rho1+rho2
call get_all_rho_grad_rho(r1,r2,rho1,rho2,grad_rho1)
rho_tot = rho1 + rho2
if(rho_tot.lt.1.d-10)rho_tot = 1.d-10
inv_rho_tot = 1.d0/rho_tot
! f(rho) = mu_r_ct * rho**beta_rho_power + mu_erf * exp(-rho)
call get_all_f_rho(rho1,rho2,mu_r_ct,mu_erf,beta_rho_power,f_rho1,d_drho_f_rho1,f_rho2)
d_dx1_f_rho1(1:3) = d_drho_f_rho1 * grad_rho1(1:3)
d_dx_rho_f_rho(1:3) = rho1 * d_dx1_f_rho1(1:3) + f_rho1 * grad_rho1(1:3)
nume = rho1 * f_rho1 + rho2 * f_rho2
mu_val = nume * inv_rho_tot
mu_der(1:3) = inv_rho_tot*inv_rho_tot * (rho_tot * d_dx_rho_f_rho(1:3) - grad_rho1(1:3) * nume)
elseif(j1b_type .eq. 203) then
! mu(r1,r2) = {rho(r1) f[rho(r1)] + rho(r2) f[rho(r2)]} / RHO
!
! RHO = rho(r1) + rho(r2)
!
! f[rho] = alpha rho^beta + mu0
!
! d/dx1 mu(r1,r2) = 1/RHO^2 * {RHO * d/dx1 (rho(r1) f[rho(r1)])
! - d/dx1 rho(r1) * [rho(r1) f[rho(r1)] + rho(r2) f[rho(r2)]] }
!
! d/dx1 f[rho(r1)] = [0.5 alpha / sqrt(rho(r1)) ] (d rho(r1) / dx1)
!
! d/dx1 (rho(r1) f[rho(r1)] = rho(r1) * d/dx1 f[rho(r1)] + f[rho(r1)] * d/dx1 rho(r1)
!!!!!!!!! rho1,rho2,rho1+rho2
call get_all_rho_grad_rho(r1,r2,rho1,rho2,grad_rho1)
rho_tot = rho1 + rho2
if(rho_tot.lt.1.d-10)rho_tot = 1.d-10
inv_rho_tot = 1.d0/rho_tot
! f(rho) = mu_r_ct * rho**beta_rho_power + mu_erf
call get_all_f_rho_simple(rho1,rho2,mu_r_ct,mu_erf,beta_rho_power,f_rho1,d_drho_f_rho1,f_rho2)
d_dx1_f_rho1(1:3) = d_drho_f_rho1 * grad_rho1(1:3)
d_dx_rho_f_rho(1:3) = rho1 * d_dx1_f_rho1(1:3) + f_rho1 * grad_rho1(1:3)
nume = rho1 * f_rho1 + rho2 * f_rho2
mu_val = nume * inv_rho_tot
mu_der(1:3) = inv_rho_tot*inv_rho_tot * (rho_tot * d_dx_rho_f_rho(1:3) - grad_rho1(1:3) * nume)
elseif(j1b_type .eq. 204) then
! mu(r1,r2) = 1/2 * (f[rho(r1)] + f[rho(r2)]}
!
! f[rho] = alpha rho^beta + mu0
!
! d/dx1 mu(r1,r2) = 1/2 * d/dx1 (rho(r1) f[rho(r1)])
!
! d/dx1 f[rho(r1)] = [0.5 alpha / sqrt(rho(r1)) ] (d rho(r1) / dx1)
!
! d/dx1 (rho(r1) f[rho(r1)] = rho(r1) * d/dx1 f[rho(r1)] + f[rho(r1)] * d/dx1 rho(r1)
!!!!!!!!! rho1,rho2,rho1+rho2
call get_all_rho_grad_rho(r1,r2,rho1,rho2,grad_rho1)
rho_tot = rho1 + rho2
if(rho_tot.lt.1.d-10)rho_tot = 1.d-10
inv_rho_tot = 1.d0/rho_tot
! f(rho) = mu_r_ct * rho**beta_rho_power + mu_erf
call get_all_f_rho_simple(rho1,rho2,mu_r_ct,mu_erf,beta_rho_power,f_rho1,d_drho_f_rho1,f_rho2)
d_dx1_f_rho1(1:3) = d_drho_f_rho1 * grad_rho1(1:3)
d_dx_rho_f_rho(1:3) = rho1 * d_dx1_f_rho1(1:3) + f_rho1 * grad_rho1(1:3)
mu_val = 0.5d0 * ( f_rho1 + f_rho2)
mu_der(1:3) = d_dx_rho_f_rho(1:3)
else
print *, ' j1b_type = ', j1b_type, 'not implemented yet'
stop
@ -684,3 +776,76 @@ end function j12_mu_square
! ---
subroutine f_mu_and_deriv_mu(rho,alpha,mu0,beta,f_mu,d_drho_f_mu)
implicit none
BEGIN_DOC
! function giving mu as a function of rho
!
! f_mu = alpha * rho**beta + mu0 * exp(-rho)
!
! and its derivative with respect to rho d_drho_f_mu
END_DOC
double precision, intent(in) :: rho,alpha,mu0,beta
double precision, intent(out) :: f_mu,d_drho_f_mu
f_mu = alpha * (rho)**beta + mu0 * dexp(-rho)
d_drho_f_mu = alpha * beta * rho**(beta-1.d0) - mu0 * dexp(-rho)
end
subroutine get_all_rho_grad_rho(r1,r2,rho1,rho2,grad_rho1)
implicit none
BEGIN_DOC
! returns the density in r1,r2 and grad_rho at r1
END_DOC
double precision, intent(in) :: r1(3),r2(3)
double precision, intent(out):: grad_rho1(3),rho1,rho2
double precision :: dm_a(1), dm_b(1), grad_dm_a(3,1), grad_dm_b(3,1)
call density_and_grad_alpha_beta(r1, dm_a, dm_b, grad_dm_a, grad_dm_b)
rho1 = dm_a(1) + dm_b(1)
grad_rho1(1:3) = grad_dm_a(1:3,1) + grad_dm_b(1:3,1)
call density_and_grad_alpha_beta(r2, dm_a, dm_b, grad_dm_a, grad_dm_b)
rho2 = dm_a(1) + dm_b(1)
end
subroutine get_all_f_rho(rho1,rho2,alpha,mu0,beta,f_rho1,d_drho_f_rho1,f_rho2)
implicit none
BEGIN_DOC
! returns the values f(mu(r1)), f(mu(r2)) and d/drho(1) f(mu(r1))
END_DOC
double precision, intent(in) :: rho1,rho2,alpha,mu0,beta
double precision, intent(out):: f_rho1,d_drho_f_rho1,f_rho2
double precision :: tmp
call f_mu_and_deriv_mu(rho1,alpha,mu0,beta,f_rho1,d_drho_f_rho1)
call f_mu_and_deriv_mu(rho2,alpha,mu0,beta,f_rho2,tmp)
end
subroutine get_all_f_rho_simple(rho1,rho2,alpha,mu0,beta,f_rho1,d_drho_f_rho1,f_rho2)
implicit none
BEGIN_DOC
! returns the values f(mu(r1)), f(mu(r2)) and d/drho(1) f(mu(r1))
END_DOC
double precision, intent(in) :: rho1,rho2,alpha,mu0,beta
double precision, intent(out):: f_rho1,d_drho_f_rho1,f_rho2
double precision :: tmp
call f_mu_and_deriv_mu_simple(rho1,alpha,mu0,beta,f_rho1,d_drho_f_rho1)
call f_mu_and_deriv_mu_simple(rho2,alpha,mu0,beta,f_rho2,tmp)
end
subroutine f_mu_and_deriv_mu_simple(rho,alpha,mu0,beta,f_mu,d_drho_f_mu)
implicit none
BEGIN_DOC
! function giving mu as a function of rho
!
! f_mu = alpha * rho**beta + mu0
!
! and its derivative with respect to rho d_drho_f_mu
END_DOC
double precision, intent(in) :: rho,alpha,mu0,beta
double precision, intent(out) :: f_mu,d_drho_f_mu
f_mu = alpha * (rho)**beta + mu0
d_drho_f_mu = alpha * beta * rho**(beta-1.d0)
end

33
src/non_h_ints_mu/plot_mu_of_r.irp.f Normal file
View File

@ -0,0 +1,33 @@
program plot_mu_of_r
implicit none
read_wf = .False.
touch read_wf
call routine_print
end
subroutine routine_print
implicit none
character*(128) :: output
integer :: i_unit_output,getUnitAndOpen
output=trim(ezfio_filename)//'.mu_of_r'
i_unit_output = getUnitAndOpen(output,'w')
integer :: ipoint,nx
double precision :: xmax,xmin,r(3),dx
double precision :: mu_val, mu_der(3),dm_a,dm_b,grad
xmax = 5.D0
xmin = -5.D0
nx = 10000
dx = (xmax - xmin)/dble(nx)
r = 0.d0
r(1) = xmin
do ipoint = 1, nx
call mu_r_val_and_grad(r, r, mu_val, mu_der)
call dm_dft_alpha_beta_at_r(r,dm_a,dm_b)
grad = mu_der(1)**2 + mu_der(2)**2 + mu_der(3)**2
grad = dsqrt(grad)
write(i_unit_output,'(100(F16.7,X))')r(1),mu_val,dm_a+dm_b,grad
r(1) += dx
enddo
end

8
src/tc_bi_ortho/dressing_vectors_lr.irp.f
View File

@ -27,7 +27,7 @@ subroutine get_delta_bitc_right(psidet, psicoef, ndet, Nint, delta)
i = 1
j = 1
call htilde_mu_mat_bi_ortho(psidet(1,1,i), psidet(1,1,j), Nint, htc_mono, htc_twoe, htc_three, htc_tot)
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 hmat_bi_ortho (psidet(1,1,i), psidet(1,1,j), Nint, h_mono, h_twoe, h_tot)
delta = 0.d0
@ -39,7 +39,7 @@ subroutine get_delta_bitc_right(psidet, psicoef, ndet, Nint, delta)
do j = 1, ndet
! < I | Htilde | J >
call htilde_mu_mat_bi_ortho(psidet(1,1,i), psidet(1,1,j), Nint, htc_mono, htc_twoe, htc_three, htc_tot)
call htilde_mu_mat_bi_ortho_slow(psidet(1,1,i), psidet(1,1,j), Nint, htc_mono, htc_twoe, htc_three, htc_tot)
! < I | H | J >
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
j = 1
call htilde_mu_mat_bi_ortho(psidet(1,1,i), psidet(1,1,j), Nint, htc_mono, htc_twoe, htc_three, htc_tot)
call htilde_mu_mat_bi_ortho_slow(psidet(1,1,i), psidet(1,1,j), Nint, htc_mono, htc_twoe, htc_three, htc_tot)
delta = 0.d0
!$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
! < I | Htilde | J >
call htilde_mu_mat_bi_ortho(psidet(1,1,i), psidet(1,1,j), Nint, htc_mono, htc_twoe, htc_three, htc_tot)
call htilde_mu_mat_bi_ortho_slow(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
enddo

18
src/tc_bi_ortho/e_corr_bi_ortho.irp.f
View File

@ -2,7 +2,7 @@
BEGIN_PROVIDER [ double precision, e_tilde_00]
implicit none
double precision :: hmono,htwoe,hthree,htot
call htilde_mu_mat_bi_ortho(HF_bitmask,HF_bitmask,N_int,hmono,htwoe,hthree,htot)
call htilde_mu_mat_bi_ortho_slow(HF_bitmask,HF_bitmask,N_int,hmono,htwoe,hthree,htot)
e_tilde_00 = htot
END_PROVIDER
@ -18,11 +18,11 @@
do i = 1, N_det
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(psi_det(1,1,i),HF_bitmask,N_int,hmono,htwoe,hthree,htilde_ij)
call htilde_mu_mat_bi_ortho(psi_det(1,1,i),psi_det(1,1,i),N_int,hmono,htwoe,hthree,e_i0)
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_bi_ortho_slow(psi_det(1,1,i),psi_det(1,1,i),N_int,hmono,htwoe,hthree,e_i0)
delta_e = e_tilde_00 - e_i0
coef_pt1 = htilde_ij / delta_e
call htilde_mu_mat_bi_ortho(HF_bitmask,psi_det(1,1,i),N_int,hmono,htwoe,hthree,htilde_ij)
call htilde_mu_mat_bi_ortho_slow(HF_bitmask,psi_det(1,1,i),N_int,hmono,htwoe,hthree,htilde_ij)
e_pt2_tc_bi_orth += coef_pt1 * htilde_ij
if(degree == 1)then
e_pt2_tc_bi_orth_single += coef_pt1 * htilde_ij
@ -37,7 +37,7 @@
BEGIN_PROVIDER [ double precision, e_tilde_bi_orth_00]
implicit none
double precision :: hmono,htwoe,hthree,htilde_ij
call htilde_mu_mat_bi_ortho(HF_bitmask,HF_bitmask,N_int,hmono,htwoe,hthree,e_tilde_bi_orth_00)
call htilde_mu_mat_bi_ortho_slow(HF_bitmask,HF_bitmask,N_int,hmono,htwoe,hthree,e_tilde_bi_orth_00)
e_tilde_bi_orth_00 += nuclear_repulsion
END_PROVIDER
@ -57,7 +57,7 @@
e_corr_double_bi_orth = 0.d0
do i = 1, N_det
call get_excitation_degree(HF_bitmask,psi_det(1,1,i),degree,N_int)
call htilde_mu_mat_bi_ortho(HF_bitmask,psi_det(1,1,i),N_int,hmono,htwoe,hthree,htilde_ij)
call htilde_mu_mat_bi_ortho_slow(HF_bitmask,psi_det(1,1,i),N_int,hmono,htwoe,hthree,htilde_ij)
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_abs += dabs(reigvec_tc_bi_orth(i,1) * htilde_ij/reigvec_tc_bi_orth(1,1))
@ -80,7 +80,7 @@
do i = 1, N_det
accu += reigvec_tc_bi_orth(i,1) * leigvec_tc_bi_orth(i,1)
do j = 1, N_det
call htilde_mu_mat_bi_ortho(psi_det(1,1,j),psi_det(1,1,i),N_int,hmono,htwoe,hthree,htilde_ij)
call htilde_mu_mat_bi_ortho_slow(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)
enddo
enddo
@ -99,8 +99,8 @@ BEGIN_PROVIDER [ double precision, coef_pt1_bi_ortho, (N_det)]
if(degree==0)then
coef_pt1_bi_ortho(i) = 1.d0
else
call htilde_mu_mat_bi_ortho(psi_det(1,1,i),HF_bitmask,N_int,hmono,htwoe,hthree,htilde_ij)
call htilde_mu_mat_bi_ortho(psi_det(1,1,i),psi_det(1,1,i),N_int,hmono,htwoe,hthree,e_i0)
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_bi_ortho_slow(psi_det(1,1,i),psi_det(1,1,i),N_int,hmono,htwoe,hthree,e_i0)
delta_e = e_tilde_00 - e_i0
coef_pt1 = htilde_ij / delta_e
coef_pt1_bi_ortho(i)= coef_pt1

12
src/tc_bi_ortho/h_tc_bi_ortho_psi.irp.f
View File

@ -1,4 +1,4 @@
subroutine htc_bi_ortho_calc_tdav(v, u, N_st, sze)
subroutine htc_bi_ortho_calc_tdav_slow(v, u, N_st, sze)
use bitmasks
@ -27,7 +27,7 @@ subroutine htc_bi_ortho_calc_tdav(v, u, N_st, sze)
i = 1
j = 1
call htilde_mu_mat_bi_ortho_tot(psi_det(1,1,i), psi_det(1,1,j), N_int, htot)
call htilde_mu_mat_bi_ortho_tot_slow(psi_det(1,1,i), psi_det(1,1,j), N_int, htot)
v = 0.d0
!$OMP PARALLEL DO DEFAULT(NONE) SCHEDULE(dynamic,8) &
@ -36,7 +36,7 @@ subroutine htc_bi_ortho_calc_tdav(v, u, N_st, sze)
do istate = 1, N_st
do i = 1, sze
do j = 1, sze
call htilde_mu_mat_bi_ortho_tot(psi_det(1,1,i), psi_det(1,1,j), N_int, htot)
call htilde_mu_mat_bi_ortho_tot_slow(psi_det(1,1,i), psi_det(1,1,j), N_int, htot)
v(i,istate) = v(i,istate) + htot * u(j,istate)
enddo
enddo
@ -45,7 +45,7 @@ subroutine htc_bi_ortho_calc_tdav(v, u, N_st, sze)
end
subroutine htcdag_bi_ortho_calc_tdav(v, u, N_st, sze)
subroutine htcdag_bi_ortho_calc_tdav_slow(v, u, N_st, sze)
use bitmasks
@ -71,7 +71,7 @@ subroutine htcdag_bi_ortho_calc_tdav(v, u, N_st, sze)
i = 1
j = 1
call htilde_mu_mat_bi_ortho_tot(psi_det(1,1,i), psi_det(1,1,j), N_int, htot)
call htilde_mu_mat_bi_ortho_tot_slow(psi_det(1,1,i), psi_det(1,1,j), N_int, htot)
v = 0.d0
@ -81,7 +81,7 @@ subroutine htcdag_bi_ortho_calc_tdav(v, u, N_st, sze)
do istate = 1, N_st
do i = 1, sze
do j = 1, sze
call htilde_mu_mat_bi_ortho_tot(psi_det(1,1,j), psi_det(1,1,i), N_int, htot)
call htilde_mu_mat_bi_ortho_tot_slow(psi_det(1,1,j), psi_det(1,1,i), N_int, htot)
v(i,istate) = v(i,istate) + htot * u(j,istate)
enddo
enddo

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

@ -49,12 +49,12 @@ subroutine routine
do i = 1, N_det
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(psi_det(1,1,i),HF_bitmask,N_int,hmono,htwoe,hthree,htilde_ij)
call htilde_mu_mat_bi_ortho(psi_det(1,1,i),psi_det(1,1,i),N_int,hmono,htwoe,hthree,e_i0)
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_bi_ortho_slow(psi_det(1,1,i),psi_det(1,1,i),N_int,hmono,htwoe,hthree,e_i0)
delta_e = e_tilde_00 - e_i0
coef_pt1 = htilde_ij / delta_e
call htilde_mu_mat_bi_ortho(HF_bitmask,psi_det(1,1,i),N_int,hmono,htwoe,hthree,htilde_ij)
call htilde_mu_mat_bi_ortho_slow(HF_bitmask,psi_det(1,1,i),N_int,hmono,htwoe,hthree,htilde_ij)
contrib_pt = coef_pt1 * htilde_ij
e_pt2 += contrib_pt

8
src/tc_bi_ortho/pt2_tc_cisd.irp.f
View File

@ -36,11 +36,11 @@ subroutine routine
e_corr_abs = 0.d0
e_corr_pos = 0.d0
e_corr_neg = 0.d0
call htilde_mu_mat_bi_ortho_tot(psi_det(1,1,1), psi_det(1,1,1), N_int, e00)
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(psi_det(1,1,i), psi_det(1,1,1), N_int, hi0)
call htilde_mu_mat_bi_ortho_tot(psi_det(1,1,1), psi_det(1,1,i), N_int, h0i)
call htilde_mu_mat_bi_ortho_tot(psi_det(1,1,i), psi_det(1,1,i), N_int, ei)
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)

24
src/tc_bi_ortho/slater_tc_3e.irp.f → src/tc_bi_ortho/slater_tc_3e_slow.irp.f
View File

@ -1,23 +1,5 @@
subroutine provide_all_three_ints_bi_ortho
implicit none
BEGIN_DOC
! routine that provides all necessary three-electron integrals
END_DOC
if(three_body_h_tc)then
PROVIDE three_e_3_idx_direct_bi_ort three_e_3_idx_cycle_1_bi_ort three_e_3_idx_cycle_2_bi_ort
PROVIDE three_e_3_idx_exch23_bi_ort three_e_3_idx_exch13_bi_ort three_e_3_idx_exch12_bi_ort
PROVIDE three_e_4_idx_direct_bi_ort three_e_4_idx_cycle_1_bi_ort three_e_4_idx_cycle_2_bi_ort
PROVIDE three_e_4_idx_exch23_bi_ort three_e_4_idx_exch13_bi_ort three_e_4_idx_exch12_bi_ort
endif
if(.not.double_normal_ord)then
PROVIDE three_e_5_idx_direct_bi_ort three_e_5_idx_cycle_1_bi_ort three_e_5_idx_cycle_2_bi_ort
PROVIDE three_e_5_idx_exch23_bi_ort three_e_5_idx_exch13_bi_ort three_e_5_idx_exch12_bi_ort
else
PROVIDE normal_two_body_bi_orth
endif
end
subroutine diag_htilde_three_body_ints_bi_ort(Nint, key_i, hthree)
subroutine diag_htilde_three_body_ints_bi_ort_slow(Nint, key_i, hthree)
BEGIN_DOC
! diagonal element of htilde ONLY FOR THREE-BODY TERMS WITH BI ORTHONORMAL ORBITALS
@ -108,7 +90,7 @@ subroutine diag_htilde_three_body_ints_bi_ort(Nint, key_i, hthree)
end
subroutine single_htilde_three_body_ints_bi_ort(Nint, key_j, key_i, hthree)
subroutine single_htilde_three_body_ints_bi_ort_slow(Nint, key_j, key_i, hthree)
BEGIN_DOC
! <key_j | H_tilde | key_i> for single excitation ONLY FOR THREE-BODY TERMS WITH BI ORTHONORMAL ORBITALS
@ -203,7 +185,7 @@ end
! ---
subroutine double_htilde_three_body_ints_bi_ort(Nint, key_j, key_i, hthree)
subroutine double_htilde_three_body_ints_bi_ort_slow(Nint, key_j, key_i, hthree)
BEGIN_DOC
! <key_j | H_tilde | key_i> for double excitation ONLY FOR THREE-BODY TERMS WITH BI ORTHONORMAL ORBITALS

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

@ -1,3 +1,26 @@
subroutine provide_all_three_ints_bi_ortho
implicit none
BEGIN_DOC
! routine that provides all necessary three-electron integrals
END_DOC
if(three_body_h_tc)then
if(three_e_3_idx_term)then
PROVIDE three_e_3_idx_direct_bi_ort three_e_3_idx_cycle_1_bi_ort three_e_3_idx_cycle_2_bi_ort
PROVIDE three_e_3_idx_exch23_bi_ort three_e_3_idx_exch13_bi_ort three_e_3_idx_exch12_bi_ort
endif
if(three_e_4_idx_term)then
PROVIDE three_e_4_idx_direct_bi_ort three_e_4_idx_cycle_1_bi_ort three_e_4_idx_cycle_2_bi_ort
PROVIDE three_e_4_idx_exch23_bi_ort three_e_4_idx_exch13_bi_ort three_e_4_idx_exch12_bi_ort
endif
if(.not.double_normal_ord.and.three_e_5_idx_term)then
PROVIDE three_e_5_idx_direct_bi_ort three_e_5_idx_cycle_1_bi_ort three_e_5_idx_cycle_2_bi_ort
PROVIDE three_e_5_idx_exch23_bi_ort three_e_5_idx_exch13_bi_ort three_e_5_idx_exch12_bi_ort
elseif (double_normal_ord .and. (.not. three_e_5_idx_term))then
PROVIDE normal_two_body_bi_orth
endif
endif
end
subroutine htilde_mu_mat_opt_bi_ortho_tot(key_j, key_i, Nint, htot)
implicit none
BEGIN_DOC

8
src/tc_bi_ortho/slater_tc_opt_diag.irp.f
View File

@ -7,11 +7,11 @@
! Various component of the TC energy for the reference "HF" Slater determinant
END_DOC
double precision :: hmono, htwoe, htot, hthree
call diag_htilde_mu_mat_bi_ortho(N_int,HF_bitmask , hmono, htwoe, htot)
call diag_htilde_mu_mat_bi_ortho_slow(N_int,HF_bitmask , hmono, htwoe, htot)
ref_tc_energy_1e = hmono
ref_tc_energy_2e = htwoe
if(three_body_h_tc)then
call diag_htilde_three_body_ints_bi_ort(N_int, HF_bitmask, hthree)
call diag_htilde_three_body_ints_bi_ort_slow(N_int, HF_bitmask, hthree)
ref_tc_energy_3e = hthree
else
ref_tc_energy_3e = 0.d0
@ -156,7 +156,7 @@ subroutine ac_tc_operator(iorb,ispin,key,hmono,htwoe,hthree,Nint,na,nb)
htwoe = htwoe + mo_bi_ortho_tc_two_e_jj(occ(i,other_spin),iorb)
enddo
if(three_body_h_tc.and.elec_num.gt.2)then
if(three_body_h_tc.and.elec_num.gt.2.and.three_e_3_idx_term)then
!!!!! 3-e part
!! same-spin/same-spin
do j = 1, na
@ -243,7 +243,7 @@ subroutine a_tc_operator(iorb,ispin,key,hmono,htwoe,hthree,Nint,na,nb)
htwoe= htwoe- mo_bi_ortho_tc_two_e_jj(occ(i,other_spin),iorb)
enddo
if(three_body_h_tc.and.elec_num.gt.2)then
if(three_body_h_tc.and.elec_num.gt.2.and.three_e_3_idx_term)then
!!!!! 3-e part
!! same-spin/same-spin
do j = 1, na

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

@ -42,13 +42,13 @@ subroutine double_htilde_mu_mat_fock_bi_ortho(Nint, key_j, key_i, hmono, htwoe,
! opposite spin two-body
htwoe = mo_bi_ortho_tc_two_e(p2,p1,h2,h1)
if(three_body_h_tc.and.elec_num.gt.2)then
if(.not.double_normal_ord)then
if(.not.double_normal_ord.and.three_e_5_idx_term)then
if(degree_i>degree_j)then
call three_comp_two_e_elem(key_j,h1,h2,p1,p2,s1,s2,hthree)
else
call three_comp_two_e_elem(key_i,h1,h2,p1,p2,s1,s2,hthree)
endif
elseif(double_normal_ord.and.elec_num.gt.2)then
elseif(double_normal_ord)then
htwoe += normal_two_body_bi_orth(p2,h2,p1,h1)
endif
endif
@ -59,13 +59,13 @@ subroutine double_htilde_mu_mat_fock_bi_ortho(Nint, key_j, key_i, hmono, htwoe,
! exchange terms
htwoe -= mo_bi_ortho_tc_two_e(p1,p2,h2,h1)
if(three_body_h_tc.and.elec_num.gt.2)then
if(.not.double_normal_ord)then
if(.not.double_normal_ord.and.three_e_5_idx_term)then
if(degree_i>degree_j)then
call three_comp_two_e_elem(key_j,h1,h2,p1,p2,s1,s2,hthree)
else
call three_comp_two_e_elem(key_i,h1,h2,p1,p2,s1,s2,hthree)
endif
elseif(double_normal_ord.and.elec_num.gt.2)then
elseif(double_normal_ord)then
htwoe -= normal_two_body_bi_orth(h2,p1,h1,p2)
htwoe += normal_two_body_bi_orth(h1,p1,h2,p2)
endif

2
src/tc_bi_ortho/slater_tc_opt_single.irp.f
View File

@ -106,7 +106,7 @@ subroutine get_single_excitation_from_fock_tc(key_i,key_j,h,p,spin,phase,hmono,h
htwoe -= buffer_x(i)
enddo
hthree = 0.d0
if (three_body_h_tc.and.elec_num.gt.2)then
if (three_body_h_tc.and.elec_num.gt.2.and.three_e_4_idx_term)then
call three_comp_fock_elem(key_i,h,p,spin,hthree)
endif

47
src/tc_bi_ortho/slater_tc.irp.f → src/tc_bi_ortho/slater_tc_slow.irp.f
View File

@ -1,7 +1,7 @@
! ---
subroutine htilde_mu_mat_bi_ortho_tot(key_j, key_i, Nint, htot)
subroutine htilde_mu_mat_bi_ortho_tot_slow(key_j, key_i, Nint, htot)
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
@ -24,14 +24,14 @@ subroutine htilde_mu_mat_bi_ortho_tot(key_j, key_i, Nint, htot)
if(degree.gt.2)then
htot = 0.d0
else
call htilde_mu_mat_bi_ortho(key_j, key_i, Nint, hmono, htwoe, hthree, htot)
call htilde_mu_mat_bi_ortho_slow(key_j, key_i, Nint, hmono, htwoe, hthree, htot)
endif
end subroutine htilde_mu_mat_bi_ortho_tot
end subroutine htilde_mu_mat_bi_ortho_tot_slow
! --
subroutine htilde_mu_mat_bi_ortho(key_j, key_i, Nint, hmono, htwoe, hthree, htot)
subroutine htilde_mu_mat_bi_ortho_slow(key_j, key_i, Nint, hmono, htwoe, hthree, htot)
BEGIN_DOC
!
@ -61,22 +61,22 @@ subroutine htilde_mu_mat_bi_ortho(key_j, key_i, Nint, hmono, htwoe, hthree, htot
if(degree.gt.2) return
if(degree == 0)then
call diag_htilde_mu_mat_bi_ortho(Nint, key_i, hmono, htwoe, htot)
call diag_htilde_mu_mat_bi_ortho_slow(Nint, key_i, hmono, htwoe, htot)
else if (degree == 1)then
call single_htilde_mu_mat_bi_ortho(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
call double_htilde_mu_mat_bi_ortho(Nint, key_j, key_i, hmono, htwoe, htot)
call double_htilde_mu_mat_bi_ortho_slow(Nint, key_j, key_i, hmono, htwoe, htot)
endif
if(three_body_h_tc) then
if(degree == 2) then
if(.not.double_normal_ord) then
call double_htilde_three_body_ints_bi_ort(Nint, key_j, key_i, hthree)
if(.not.double_normal_ord.and.elec_num.gt.2.and.three_e_5_idx_term) then
call double_htilde_three_body_ints_bi_ort_slow(Nint, key_j, key_i, hthree)
endif
else if(degree == 1) then
call single_htilde_three_body_ints_bi_ort(Nint, key_j, key_i, hthree)
else if(degree == 0) then
call diag_htilde_three_body_ints_bi_ort(Nint, key_i, hthree)
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)
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)
endif
endif
@ -89,7 +89,7 @@ end
! ---
subroutine diag_htilde_mu_mat_bi_ortho(Nint, key_i, hmono, htwoe, htot)
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
@ -188,7 +188,7 @@ end
subroutine double_htilde_mu_mat_bi_ortho(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
! <key_j | H_tilde | key_i> for double excitation ONLY FOR ONE- AND TWO-BODY TERMS
@ -227,18 +227,7 @@ subroutine double_htilde_mu_mat_bi_ortho(Nint, key_j, key_i, hmono, htwoe, htot)
return
endif
! if(core_tc_op)then
! print*,'core_tc_op not already taken into account for bi ortho'
! print*,'stopping ...'
! stop
! 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
call get_double_excitation(key_i, key_j, exc, phase, Nint)
call decode_exc(exc, 2, h1, p1, h2, p2, s1, s2)
@ -246,7 +235,7 @@ subroutine double_htilde_mu_mat_bi_ortho(Nint, key_j, key_i, hmono, htwoe, htot)
! opposite spin two-body
! key_j, key_i
htwoe = mo_bi_ortho_tc_two_e(p2,p1,h2,h1)
if(double_normal_ord.and.+Ne(1).gt.2)then
if(three_body_h_tc.and.double_normal_ord.and.+Ne(1).gt.2)then
htwoe += normal_two_body_bi_orth(p2,h2,p1,h1)!!! WTF ???
endif
else
@ -255,7 +244,7 @@ subroutine double_htilde_mu_mat_bi_ortho(Nint, key_j, key_i, hmono, htwoe, htot)
htwoe = mo_bi_ortho_tc_two_e(p2,p1,h2,h1)
! exchange terms
htwoe -= mo_bi_ortho_tc_two_e(p1,p2,h2,h1)
if(double_normal_ord.and.+Ne(1).gt.2)then
if(three_body_h_tc.and.double_normal_ord.and.+Ne(1).gt.2)then
htwoe -= normal_two_body_bi_orth(h2,p1,h1,p2)!!! WTF ???
htwoe += normal_two_body_bi_orth(h1,p1,h2,p2)!!! WTF ???
endif
@ -266,7 +255,7 @@ subroutine double_htilde_mu_mat_bi_ortho(Nint, key_j, key_i, hmono, htwoe, htot)
end
subroutine single_htilde_mu_mat_bi_ortho(Nint, key_j, key_i, hmono, htwoe, htot)
subroutine single_htilde_mu_mat_bi_ortho_slow(Nint, key_j, key_i, hmono, htwoe, htot)
BEGIN_DOC
! <key_j | H_tilde | key_i> for single excitation ONLY FOR ONE- AND TWO-BODY TERMS

12
src/tc_bi_ortho/tc_cisd_sc2_utils.irp.f
View File

@ -11,10 +11,10 @@
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(psi_det(1,1,i), psi_det(1,1,i), N_int, H_jj(i))
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(HF_bitmask,psi_det(1,1,i),N_int,hmono,htwoe,hthree,h0j(i))
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
@ -29,7 +29,7 @@
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)
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
@ -41,8 +41,8 @@
it = 0
dressing_dets = 0.d0
double precision, allocatable :: H_jj(:),vec_tmp(:,:),eigval_tmp(:)
external htc_bi_ortho_calc_tdav
external htcdag_bi_ortho_calc_tdav
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
@ -66,7 +66,7 @@
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)
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

5
src/tc_bi_ortho/tc_h_eigvectors.irp.f
View File

@ -207,8 +207,6 @@ end
else ! n_det > N_det_max_full
double precision, allocatable :: H_jj(:),vec_tmp(:,:)
external htc_bi_ortho_calc_tdav
external htcdag_bi_ortho_calc_tdav
external H_tc_u_0_opt
external H_tc_dagger_u_0_opt
external H_tc_s2_dagger_u_0_opt
@ -217,7 +215,7 @@ end
allocate(H_jj(N_det),vec_tmp(N_det,n_states_diag))
do i = 1, N_det
call htilde_mu_mat_bi_ortho_tot(psi_det(1,1,i), psi_det(1,1,i), N_int, H_jj(i))
call htilde_mu_mat_opt_bi_ortho_tot(psi_det(1,1,i), psi_det(1,1,i), N_int, H_jj(i))
enddo
print*,'---------------------------------'
@ -259,7 +257,6 @@ end
do istate = N_states+1, n_states_diag
vec_tmp(istate,istate) = 1.d0
enddo
!call davidson_general_ext_rout_nonsym_b1space(vec_tmp, H_jj, eigval_right_tc_bi_orth, N_det, n_states, n_states_diag, converged, htc_bi_ortho_calc_tdav)
!call davidson_general_ext_rout_nonsym_b1space(vec_tmp, H_jj, eigval_right_tc_bi_orth, N_det, n_states, n_states_diag, converged, H_tc_u_0_opt)
converged = .False.
i_it = 0

16
src/tc_bi_ortho/tc_hmat.irp.f
View File

@ -9,33 +9,25 @@
implicit none
integer :: i, j
double precision :: hmono,htwoe,hthree,htot
double precision :: htot
PROVIDE N_int
i = 1
j = 1
call htilde_mu_mat_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_tot(psi_det(1,1,j), psi_det(1,1,i), N_int, htot)
!$OMP PARALLEL DO SCHEDULE(GUIDED) DEFAULT(NONE) PRIVATE(i,j,hmono, htwoe, hthree, htot) &
!$OMP PARALLEL DO SCHEDULE(GUIDED) DEFAULT(NONE) PRIVATE(i,j, htot) &
!$OMP SHARED (N_det, psi_det, N_int,htilde_matrix_elmt_bi_ortho)
do i = 1, N_det
do j = 1, N_det
! < J | Htilde | I >
call htilde_mu_mat_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_tot(psi_det(1,1,j), psi_det(1,1,i), N_int, htot)
!print *, ' hmono = ', hmono
!print *, ' htwoe = ', htwoe
!print *, ' hthree = ', hthree
htilde_matrix_elmt_bi_ortho(j,i) = htot
enddo
enddo
!$OMP END PARALLEL DO
! print*,'htilde_matrix_elmt_bi_ortho = '
! do i = 1, min(100,N_det)
! write(*,'(100(F16.10,X))')htilde_matrix_elmt_bi_ortho(1:min(100,N_det),i)
! enddo
END_PROVIDER

4
src/tc_bi_ortho/tc_som.irp.f
View File

@ -56,8 +56,8 @@ subroutine main()
U_SOM = 0.d0
do i = 1, N_det
if(i == i_HF) cycle
call htilde_mu_mat_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(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_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_bi_ortho_slow(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
enddo
U_SOM = 0.5d0 * U_SOM

2
src/tc_bi_ortho/tc_utils.irp.f
View File

@ -12,7 +12,7 @@ subroutine write_tc_energy()
do i = 1, N_det
do j = 1, N_det
!htot = htilde_matrix_elmt_bi_ortho(i,j)
call htilde_mu_mat_bi_ortho(psi_det(1,1,i), psi_det(1,1,j), N_int, hmono, htwoe, hthree, htot)
call htilde_mu_mat_bi_ortho_slow(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 + leigvec_tc_bi_orth(i,k) * reigvec_tc_bi_orth(j,k) * htot
enddo

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

@ -35,7 +35,7 @@ subroutine test
det_i = ref_bitmask
call do_single_excitation(det_i,h1,p1,s1,i_ok)
call do_single_excitation(det_i,h2,p2,s2,i_ok)
call htilde_mu_mat_bi_ortho(det_i,HF_bitmask,N_int,hmono,htwoe,hthree,htilde_ij)
call htilde_mu_mat_bi_ortho_slow(det_i,HF_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)
hthree *= phase
@ -67,7 +67,7 @@ do h1 = 1, elec_alpha_num
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(det_i,ref_bitmask,N_int,hmono,htwoe,hthree,htilde_ij)
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
@ -103,7 +103,7 @@ do h1 = 1, elec_beta_num
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(det_i,ref_bitmask,N_int,hmono,htwoe,hthree,htilde_ij)
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)

2
src/tc_bi_ortho/test_s2_tc.irp.f
View File

@ -91,7 +91,7 @@ subroutine routine_test_s2_davidson
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(psi_det(1,1,i), psi_det(1,1,i), N_int, H_jj(i))
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 '

14
src/tc_bi_ortho/test_tc_bi_ortho.irp.f
View File

@ -31,7 +31,7 @@ subroutine test_h_u0
u_0(i) = psi_r_coef_bi_ortho(i,1)
enddo
call H_tc_u_0_nstates_openmp(v_0_new,u_0,N_states,N_det, do_right)
call htc_bi_ortho_calc_tdav (v_0_ref,u_0,N_states,N_det)
call htc_bi_ortho_calc_tdav_slow (v_0_ref,u_0,N_states,N_det)
print*,'difference right '
accu = 0.d0
do i = 1, N_det
@ -42,7 +42,7 @@ subroutine test_h_u0
do_right = .False.
v_0_new = 0.d0
call H_tc_u_0_nstates_openmp(v_0_new,u_0,N_states,N_det, do_right)
call htcdag_bi_ortho_calc_tdav(v_0_ref_dagger,u_0,N_states,N_det, do_right)
call htcdag_bi_ortho_calc_tdav_slow(v_0_ref_dagger,u_0,N_states,N_det, do_right)
print*,'difference left'
accu = 0.d0
do i = 1, N_det
@ -63,7 +63,7 @@ subroutine test_slater_tc_opt
i_count = 0.d0
do i = 1, N_det
do j = 1,N_det
call htilde_mu_mat_bi_ortho(psi_det(1,1,j), psi_det(1,1,i), N_int, hmono, htwoe, hthree, htot)
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, hnewmono, hnewtwoe, hnewthree, hnewtot)
if(dabs(htot).gt.1.d-15)then
i_count += 1.D0
@ -99,7 +99,7 @@ subroutine timing_tot
do j = 1, N_det
! call get_excitation_degree(psi_det(1,1,j), psi_det(1,1,i),degree,N_int)
i_count += 1.d0
call htilde_mu_mat_bi_ortho(psi_det(1,1,j), psi_det(1,1,i), N_int, hmono, htwoe, hthree, htot)
call htilde_mu_mat_bi_ortho_slow(psi_det(1,1,j), psi_det(1,1,i), N_int, hmono, htwoe, hthree, htot)
enddo
enddo
call wall_time(wall1)
@ -146,7 +146,7 @@ subroutine timing_diag
do i = 1, N_det
do j = i,i
i_count += 1.d0
call htilde_mu_mat_bi_ortho(psi_det(1,1,j), psi_det(1,1,i), N_int, hmono, htwoe, hthree, htot)
call htilde_mu_mat_bi_ortho_slow(psi_det(1,1,j), psi_det(1,1,i), N_int, hmono, htwoe, hthree, htot)
enddo
enddo
call wall_time(wall1)
@ -183,7 +183,7 @@ subroutine timing_single
if(degree.ne.1)cycle
i_count += 1.d0
call wall_time(wall0)
call htilde_mu_mat_bi_ortho(psi_det(1,1,j), psi_det(1,1,i), N_int, hmono, htwoe, hthree, htot)
call htilde_mu_mat_bi_ortho_slow(psi_det(1,1,j), psi_det(1,1,i), N_int, hmono, htwoe, hthree, htot)
call wall_time(wall1)
accu += wall1 - wall0
enddo
@ -225,7 +225,7 @@ subroutine timing_double
if(degree.ne.2)cycle
i_count += 1.d0
call wall_time(wall0)
call htilde_mu_mat_bi_ortho(psi_det(1,1,j), psi_det(1,1,i), N_int, hmono, htwoe, hthree, htot)
call htilde_mu_mat_bi_ortho_slow(psi_det(1,1,j), psi_det(1,1,i), N_int, hmono, htwoe, hthree, htot)
call wall_time(wall1)
accu += wall1 - wall0
enddo

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

@ -25,8 +25,7 @@ subroutine test_3e
implicit none
double precision :: integral_aaa,integral_aab,integral_abb,integral_bbb,accu
double precision :: hmono, htwoe, hthree, htot
call htilde_mu_mat_bi_ortho(ref_bitmask, ref_bitmask, N_int, hmono, htwoe, hthree, htot)
! call diag_htilde_three_body_ints_bi_ort(N_int, ref_bitmask, hthree)
call htilde_mu_mat_bi_ortho_slow(ref_bitmask, ref_bitmask, N_int, hmono, htwoe, hthree, htot)
print*,'hmono = ',hmono
print*,'htwoe = ',htwoe
print*,'hthree= ',hthree
@ -88,7 +87,7 @@ subroutine routine_3()
print*, ' excited det'
call debug_det(det_i, N_int)
call htilde_mu_mat_bi_ortho(det_i, ref_bitmask, N_int, hmono, htwoe, hthree, htilde_ij)
call htilde_mu_mat_bi_ortho_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
@ -156,7 +155,7 @@ subroutine routine_tot()
stop
endif
call htilde_mu_mat_bi_ortho(det_i, ref_bitmask, N_int, hmono, htwoe, hthree, htilde_ij)
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'

24
src/tc_keywords/EZFIO.cfg
View File

@ -16,6 +16,24 @@ doc: If |true|, three-body terms are included
interface: ezfio,provider,ocaml
default: True
[three_e_3_idx_term]
type: logical
doc: If |true|, the diagonal 3-idx terms of the 3-e interaction are taken
interface: ezfio,provider,ocaml
default: True
[three_e_4_idx_term]
type: logical
doc: If |true|, the off-diagonal 4-idx terms of the 3-e interaction are taken
interface: ezfio,provider,ocaml
default: True
[three_e_5_idx_term]
type: logical
doc: If |true|, the off-diagonal 5-idx terms of the 3-e interaction are taken
interface: ezfio,provider,ocaml
default: True
[pure_three_body_h_tc]
type: logical
doc: If |true|, pure triple excitation three-body terms are included
@ -130,6 +148,12 @@ doc: a parameter used to define mu(r)
interface: ezfio, provider, ocaml
default: 6.203504908994001e-1
[beta_rho_power]
type: double precision
doc: a parameter used to define mu(r)
interface: ezfio, provider, ocaml
default: 0.5
[thr_degen_tc]
type: Threshold
doc: Threshold to determine if two orbitals are degenerate in TCSCF in order to avoid random quasi orthogonality between the right- and left-eigenvector for the same eigenvalue

30
src/trexio/EZFIO.cfg
View File

@ -10,11 +10,17 @@ doc: Name of the exported TREXIO file
interface: ezfio, ocaml, provider
default: None
[export_rdm]
[export_basis]
type: logical
doc: If True, export two-body reduced density matrix
doc: If True, export basis set and AOs
interface: ezfio, ocaml, provider
default: False
default: True
[export_mos]
type: logical
doc: If True, export basis set and AOs
interface: ezfio, ocaml, provider
default: True
[export_ao_one_e_ints]
type: logical
@ -22,12 +28,6 @@ doc: If True, export one-electron integrals in AO basis
interface: ezfio, ocaml, provider
default: False
[export_mo_one_e_ints]
type: logical
doc: If True, export one-electron integrals in MO basis
interface: ezfio, ocaml, provider
default: False
[export_ao_two_e_ints]
type: logical
doc: If True, export two-electron integrals in AO basis
@ -40,6 +40,12 @@ doc: If True, export Cholesky-decomposed two-electron integrals in AO basis
interface: ezfio, ocaml, provider
default: False
[export_mo_one_e_ints]
type: logical
doc: If True, export one-electron integrals in MO basis
interface: ezfio, ocaml, provider
default: False
[export_mo_two_e_ints]
type: logical
doc: If True, export two-electron integrals in MO basis
@ -52,3 +58,9 @@ doc: If True, export Cholesky-decomposed two-electron integrals in MO basis
interface: ezfio, ocaml, provider
default: False
[export_rdm]
type: logical
doc: If True, export two-body reduced density matrix
interface: ezfio, ocaml, provider
default: False

2
src/trexio/export_trexio.irp.f
View File

@ -2,6 +2,6 @@ program export_trexio_prog
implicit none
read_wf = .True.
SOFT_TOUCH read_wf
call export_trexio
call export_trexio(.False.)
end

34
src/trexio/export_trexio_routines.irp.f
View File

@ -1,15 +1,17 @@
subroutine export_trexio
subroutine export_trexio(update)
use trexio
implicit none
BEGIN_DOC
! Exports the wave function in TREXIO format
END_DOC
logical, intent(in) :: update
integer(trexio_t) :: f(N_states) ! TREXIO file handle
integer(trexio_exit_code) :: rc
integer :: k
double precision, allocatable :: factor(:)
character*(256) :: filenames(N_states)
character :: rw
filenames(1) = trexio_filename
do k=2,N_states
@ -18,15 +20,26 @@ subroutine export_trexio
do k=1,N_states
print *, 'TREXIO file : ', trim(filenames(k))
if (update) then
call system('test -f '//trim(filenames(k))//' && cp -r '//trim(filenames(k))//' '//trim(filenames(k))//'.bak')
else
call system('test -f '//trim(filenames(k))//' && mv '//trim(filenames(k))//' '//trim(filenames(k))//'.bak')
endif
enddo
print *, ''
if (update) then
rw = 'u'
else
rw = 'w'
endif
do k=1,N_states
if (backend == 0) then
f(k) = trexio_open(filenames(k), 'u', TREXIO_HDF5, rc)
f(k) = trexio_open(filenames(k), rw, TREXIO_HDF5, rc)
else if (backend == 1) then
f(k) = trexio_open(filenames(k), 'u', TREXIO_TEXT, rc)
f(k) = trexio_open(filenames(k), rw, TREXIO_TEXT, rc)
endif
if (f(k) == 0_8) then
print *, 'Unable to open TREXIO file for writing'
@ -171,12 +184,13 @@ subroutine export_trexio
endif
if (export_basis) then
! Basis
! -----
print *, 'Basis'
rc = trexio_write_basis_type(f(1), 'Gaussian', len('Gaussian'))
call trexio_assert(rc, TREXIO_SUCCESS)
@ -193,11 +207,11 @@ subroutine export_trexio
call trexio_assert(rc, TREXIO_SUCCESS)
allocate(factor(shell_num))
if (ao_normalized) then
factor(1:shell_num) = shell_normalization_factor(1:shell_num)
else
! if (ao_normalized) then
! factor(1:shell_num) = shell_normalization_factor(1:shell_num)
! else
factor(1:shell_num) = 1.d0
endif
! endif
rc = trexio_write_basis_shell_factor(f(1), factor)
call trexio_assert(rc, TREXIO_SUCCESS)
@ -258,6 +272,8 @@ subroutine export_trexio
call trexio_assert(rc, TREXIO_SUCCESS)
deallocate(factor)
endif
! One-e AO integrals
! ------------------
@ -375,6 +391,7 @@ subroutine export_trexio
! Molecular orbitals
! ------------------
if (export_mos) then
print *, 'MOs'
rc = trexio_write_mo_type(f(1), mo_label, len(trim(mo_label)))
@ -396,6 +413,7 @@ subroutine export_trexio
rc = trexio_write_mo_class(f(1), mo_class, len(mo_class(1)))
call trexio_assert(rc, TREXIO_SUCCESS)
endif
! One-e MO integrals
! ------------------

78
src/trexio/import_trexio_integrals.irp.f
View File

@ -3,6 +3,7 @@ program import_integrals_ao
implicit none
integer(trexio_t) :: f ! TREXIO file handle
integer(trexio_exit_code) :: rc
PROVIDE mo_num
f = trexio_open(trexio_filename, 'r', TREXIO_AUTO, rc)
if (f == 0_8) then
@ -42,10 +43,10 @@ subroutine run(f)
if (trexio_has_nucleus_repulsion(f) == TREXIO_SUCCESS) then
rc = trexio_read_nucleus_repulsion(f, s)
call trexio_assert(rc, TREXIO_SUCCESS)
if (rc /= TREXIO_SUCCESS) then
print *, irp_here, rc
print *, 'Error reading nuclear repulsion'
call trexio_assert(rc, TREXIO_SUCCESS)
stop -1
endif
call ezfio_set_nuclei_nuclear_repulsion(s)
@ -63,6 +64,7 @@ subroutine run(f)
if (rc /= TREXIO_SUCCESS) then
print *, irp_here
print *, 'Error reading AO overlap'
call trexio_assert(rc, TREXIO_SUCCESS)
stop -1
endif
call ezfio_set_ao_one_e_ints_ao_integrals_overlap(A)
@ -74,6 +76,7 @@ subroutine run(f)
if (rc /= TREXIO_SUCCESS) then
print *, irp_here
print *, 'Error reading AO kinetic integrals'
call trexio_assert(rc, TREXIO_SUCCESS)
stop -1
endif
call ezfio_set_ao_one_e_ints_ao_integrals_kinetic(A)
@ -85,6 +88,7 @@ subroutine run(f)
! if (rc /= TREXIO_SUCCESS) then
! print *, irp_here
! print *, 'Error reading AO ECP local integrals'
! call trexio_assert(rc, TREXIO_SUCCESS)
! stop -1
! endif
! call ezfio_set_ao_one_e_ints_ao_integrals_pseudo(A)
@ -96,6 +100,7 @@ subroutine run(f)
if (rc /= TREXIO_SUCCESS) then
print *, irp_here
print *, 'Error reading AO potential N-e integrals'
call trexio_assert(rc, TREXIO_SUCCESS)
stop -1
endif
call ezfio_set_ao_one_e_ints_ao_integrals_n_e(A)
@ -106,6 +111,10 @@ subroutine run(f)
! AO 2e integrals
! ---------------
rc = trexio_has_ao_2e_int(f)
PROVIDE ao_num
if (rc /= TREXIO_HAS_NOT) then
PROVIDE ao_integrals_map
integer*4 :: BUFSIZE
@ -143,4 +152,71 @@ subroutine run(f)
call map_save_to_disk(trim(ezfio_filename)//'/work/ao_ints',ao_integrals_map)
call ezfio_set_ao_two_e_ints_io_ao_two_e_integrals('Read')
deallocate(buffer_i, buffer_values, Vi, V)
print *, 'AO integrals read from TREXIO file'
else
print *, 'AO integrals not found in TREXIO file'
endif
! MO integrals
! ------------
allocate(A(mo_num, mo_num))
if (trexio_has_mo_1e_int_core_hamiltonian(f) == TREXIO_SUCCESS) then
rc = trexio_read_mo_1e_int_core_hamiltonian(f, A)
if (rc /= TREXIO_SUCCESS) then
print *, irp_here
print *, 'Error reading MO 1e integrals'
call trexio_assert(rc, TREXIO_SUCCESS)
stop -1
endif
call ezfio_set_mo_one_e_ints_mo_one_e_integrals(A)
call ezfio_set_mo_one_e_ints_io_mo_one_e_integrals('Read')
endif
deallocate(A)
! MO 2e integrals
! ---------------
rc = trexio_has_mo_2e_int(f)
if (rc /= TREXIO_HAS_NOT) then
BUFSIZE=mo_num**2
allocate(buffer_i(BUFSIZE), buffer_values(BUFSIZE))
allocate(Vi(4,BUFSIZE), V(BUFSIZE))
offset = 0_8
icount = BUFSIZE
rc = TREXIO_SUCCESS
do while (icount == size(V))
rc = trexio_read_mo_2e_int_eri(f, offset, icount, Vi, V)
do m=1,icount
i = Vi(1,m)
j = Vi(2,m)
k = Vi(3,m)
l = Vi(4,m)
integral = V(m)
call two_e_integrals_index(i, j, k, l, buffer_i(m) )
buffer_values(m) = integral
enddo
call map_append(mo_integrals_map, buffer_i, buffer_values, int(icount,4))
offset = offset + icount
if (rc /= TREXIO_SUCCESS) then
exit
endif
end do
n_integrals = offset
call map_sort(mo_integrals_map)
call map_unique(mo_integrals_map)
call map_save_to_disk(trim(ezfio_filename)//'/work/mo_ints',mo_integrals_map)
call ezfio_set_mo_two_e_ints_io_mo_two_e_integrals('Read')
deallocate(buffer_i, buffer_values, Vi, V)
print *, 'MO integrals read from TREXIO file'
else
print *, 'MO integrals not found in TREXIO file'
endif
end

295
src/utils/integration.irp.f
View File

@ -468,8 +468,6 @@ end subroutine
subroutine multiply_poly(b,nb,c,nc,d,nd)
implicit none
BEGIN_DOC
@ -484,33 +482,292 @@ subroutine multiply_poly(b,nb,c,nc,d,nd)
integer :: ndtmp
integer :: ib, ic, id, k
if(ior(nc,nb) >= 0) then ! True if nc>=0 and nb>=0
continue
else
if(ior(nc,nb) < 0) return !False if nc>=0 and nb>=0
select case (nb)
case (0)
call multiply_poly_b0(b,c,nc,d,nd)
return
endif
ndtmp = nb+nc
case (1)
call multiply_poly_b1(b,c,nc,d,nd)
return
case (2)
call multiply_poly_b2(b,c,nc,d,nd)
return
end select
select case (nc)
case (0)
call multiply_poly_c0(b,nb,c,d,nd)
return
case (1)
call multiply_poly_c1(b,nb,c,d,nd)
return
case (2)
call multiply_poly_c2(b,nb,c,d,nd)
return
end select
do ib=0,nb
do ic = 0,nc
d(ib+ic) = d(ib+ic) + c(ic) * b(ib)
enddo
enddo
do nd = nb+nc,0,-1
if (d(nd) /= 0.d0) exit
enddo
end
subroutine multiply_poly_b0(b,c,nc,d,nd)
implicit none
BEGIN_DOC
! Multiply two polynomials
! D(t) += B(t)*C(t)
END_DOC
integer, intent(in) :: nc
integer, intent(out) :: nd
double precision, intent(in) :: b(0:0), c(0:nc)
double precision, intent(inout) :: d(0:nc)
integer :: ndtmp
integer :: ic, id, k
if(nc < 0) return !False if nc>=0
do ic = 0,nc
d(ic) = d(ic) + c(ic) * b(0)
enddo
do ib=1,nb
d(ib) = d(ib) + c(0) * b(ib)
do ic = 1,nc
d(ib+ic) = d(ib+ic) + c(ic) * b(ib)
enddo
enddo
do nd = ndtmp,0,-1
if (d(nd) == 0.d0) then
cycle
endif
exit
do nd = nc,0,-1
if (d(nd) /= 0.d0) exit
enddo
end
subroutine multiply_poly_b1(b,c,nc,d,nd)
implicit none
BEGIN_DOC
! Multiply two polynomials
! D(t) += B(t)*C(t)
END_DOC
integer, intent(in) :: nc
integer, intent(out) :: nd
double precision, intent(in) :: b(0:1), c(0:nc)
double precision, intent(inout) :: d(0:1+nc)
integer :: ndtmp
integer :: ib, ic, id, k
if(nc < 0) return !False if nc>=0
select case (nc)
case (0)
d(0) = d(0) + c(0) * b(0)
d(1) = d(1) + c(0) * b(1)
case (1)
d(0) = d(0) + c(0) * b(0)
d(1) = d(1) + c(0) * b(1) + c(1) * b(0)
d(2) = d(2) + c(1) * b(1)
case default
d(0) = d(0) + c(0) * b(0)
do ic = 1,nc
d(ic) = d(ic) + c(ic) * b(0) + c(ic-1) * b(1)
enddo
d(nc+1) = d(nc+1) + c(nc) * b(1)
end select
do nd = 1+nc,0,-1
if (d(nd) /= 0.d0) exit
enddo
end
subroutine multiply_poly_b2(b,c,nc,d,nd)
implicit none
BEGIN_DOC
! Multiply two polynomials
! D(t) += B(t)*C(t)
END_DOC
integer, intent(in) :: nc
integer, intent(out) :: nd
double precision, intent(in) :: b(0:2), c(0:nc)
double precision, intent(inout) :: d(0:2+nc)
integer :: ndtmp
integer :: ib, ic, id, k
if(nc < 0) return !False if nc>=0
select case (nc)
case (0)
d(0) = d(0) + c(0) * b(0)
d(1) = d(1) + c(0) * b(1)
d(2) = d(2) + c(0) * b(2)
case (1)
d(0) = d(0) + c(0) * b(0)
d(1) = d(1) + c(0) * b(1) + c(1) * b(0)
d(2) = d(2) + c(0) * b(2) + c(1) * b(1)
d(3) = d(3) + c(1) * b(2)
case (2)
d(0) = d(0) + c(0) * b(0)
d(1) = d(1) + c(0) * b(1) + c(1) * b(0)
d(2) = d(2) + c(0) * b(2) + c(1) * b(1) + c(2) * b(0)
d(3) = d(3) + c(2) * b(1) + c(1) * b(2)
d(4) = d(4) + c(2) * b(2)
case default
d(0) = d(0) + c(0) * b(0)
d(1) = d(1) + c(0) * b(1) + c(1) * b(0)
do ic = 2,nc
d(ic) = d(ic) + c(ic) * b(0) + c(ic-1) * b(1) + c(ic-2) * b(2)
enddo
d(nc+1) = d(nc+1) + c(nc) * b(1) + c(nc-1) * b(2)
d(nc+2) = d(nc+2) + c(nc) * b(2)
end select
do nd = 2+nc,0,-1
if (d(nd) /= 0.d0) exit
enddo
end
subroutine multiply_poly_c0(b,nb,c,d,nd)
implicit none
BEGIN_DOC
! Multiply two polynomials
! D(t) += B(t)*C(t)
END_DOC
integer, intent(in) :: nb
integer, intent(out) :: nd
double precision, intent(in) :: b(0:nb), c(0:0)
double precision, intent(inout) :: d(0:nb)
integer :: ndtmp
integer :: ib, ic, id, k
if(nb < 0) return !False if nb>=0
do ib=0,nb
d(ib) = d(ib) + c(0) * b(ib)
enddo
do nd = nb,0,-1
if (d(nd) /= 0.d0) exit
enddo
end
subroutine multiply_poly_c1(b,nb,c,d,nd)
implicit none
BEGIN_DOC
! Multiply two polynomials
! D(t) += B(t)*C(t)
END_DOC
integer, intent(in) :: nb
integer, intent(out) :: nd
double precision, intent(in) :: b(0:nb), c(0:1)
double precision, intent(inout) :: d(0:nb+1)
integer :: ndtmp
integer :: ib, ic, id, k
if(nb < 0) return !False if nb>=0
select case (nb)
case (0)
d(0) = d(0) + c(0) * b(0)
d(1) = d(1) + c(1) * b(0)
case (1)
d(0) = d(0) + c(0) * b(0)
d(1) = d(1) + c(0) * b(1) + c(1) * b(0)
d(2) = d(2) + c(1) * b(1)
case default
d(0) = d(0) + c(0) * b(0)
do ib=1,nb
d(ib) = d(ib) + c(0) * b(ib) + c(1) * b(ib-1)
enddo
d(nb+1) = d(nb+1) + c(1) * b(nb)
end select
do nd = nb+1,0,-1
if (d(nd) /= 0.d0) exit
enddo
end
subroutine multiply_poly_c2(b,nb,c,d,nd)
implicit none
BEGIN_DOC
! Multiply two polynomials
! D(t) += B(t)*C(t)
END_DOC
integer, intent(in) :: nb
integer, intent(out) :: nd
double precision, intent(in) :: b(0:nb), c(0:2)
double precision, intent(inout) :: d(0:nb+2)
integer :: ndtmp
integer :: ib, ic, id, k
if(nb < 0) return !False if nb>=0
select case (nb)
case (0)
d(0) = d(0) + c(0) * b(0)
d(1) = d(1) + c(1) * b(0)
d(2) = d(2) + c(2) * b(0)
case (1)
d(0) = d(0) + c(0) * b(0)
d(1) = d(1) + c(0) * b(1) + c(1) * b(0)
d(2) = d(2) + c(1) * b(1) + c(2) * b(0)
d(3) = d(3) + c(2) * b(1)
case (2)
d(0) = d(0) + c(0) * b(0)
d(1) = d(1) + c(0) * b(1) + c(1) * b(0)
d(2) = d(2) + c(0) * b(2) + c(1) * b(1) + c(2) * b(0)
d(3) = d(3) + c(1) * b(2) + c(2) * b(1)
d(4) = d(4) + c(2) * b(2)
case default
d(0) = d(0) + c(0) * b(0)
d(1) = d(1) + c(0) * b(1) + c(1) * b(0)
do ib=2,nb
d(ib) = d(ib) + c(0) * b(ib) + c(1) * b(ib-1) + c(2) * b(ib-2)
enddo
d(nb+1) = d(nb+1) + c(1) * b(nb) + c(2) * b(nb-1)
d(nb+2) = d(nb+2) + c(2) * b(nb)
end select
do nd = nb+2,0,-1
if (d(nd) /= 0.d0) exit
enddo
end
subroutine multiply_poly_v(b,nb,c,nc,d,nd,n_points)
implicit none
BEGIN_DOC

24
src/utils/linear_algebra.irp.f
View File

@ -1825,39 +1825,37 @@ subroutine pivoted_cholesky( A, rank, tol, ndim, U)
!
integer :: ndim
integer, intent(inout) :: rank
double precision, dimension(ndim, ndim), intent(inout) :: A
double precision, dimension(ndim, rank), intent(out) :: U
double precision, intent(inout) :: A(ndim, ndim)
double precision, intent(out) :: U(ndim, rank)
double precision, intent(in) :: tol
integer, dimension(:), allocatable :: piv
double precision, dimension(:), allocatable :: work
character, parameter :: uplo = "U"
integer :: N, LDA
integer :: LDA
integer :: info
integer :: k, l, rank0
external :: dpstrf
rank0 = rank
N = size(A, dim=1)
LDA = N
allocate(piv(N))
allocate(work(2*N))
call dpstrf(uplo, N, A, LDA, piv, rank, tol, work, info)
LDA = ndim
allocate(piv(ndim))
allocate(work(2*ndim))
call dpstrf(uplo, ndim, A, LDA, piv, rank, tol, work, info)
if (rank > rank0) then
print *, 'Bug: rank > rank0 in pivoted cholesky. Increase rank before calling'
stop
end if
do k = 1, N
A(k+1:, k) = 0.00D+0
do k = 1, ndim
A(k+1:ndim, k) = 0.00D+0
end do
! TODO: It should be possible to use only one vector of size (1:rank) as a buffer
! to do the swapping in-place
U(:,:) = 0.00D+0
do k = 1, N
do k = 1, ndim
l = piv(k)
U(l, :) = A(1:rank, k)
U(l, 1:rank) = A(1:rank, k)
end do
end subroutine pivoted_cholesky

5
src/utils_cc/energy.irp.f
View File

@ -5,9 +5,8 @@ subroutine det_energy(det,energy)
integer(bit_kind), intent(in) :: det
double precision, intent(out) :: energy
double precision, external :: diag_H_mat_elem
call i_H_j(det,det,N_int,energy)
energy = energy + nuclear_repulsion
energy = diag_H_mat_elem(det,N_int) + nuclear_repulsion
end

69
src/utils_cc/mo_integrals_cc.irp.f
View File

@ -45,61 +45,64 @@ subroutine gen_v_space(n1,n2,n3,n4,list1,list2,list3,list4,v)
integer, intent(in) :: list1(n1),list2(n2),list3(n3),list4(n4)
double precision, intent(out) :: v(n1,n2,n3,n4)
integer :: i1,i2,i3,i4,idx1,idx2,idx3,idx4
double precision :: get_two_e_integral
PROVIDE mo_two_e_integrals_in_map
integer :: i1,i2,i3,i4,idx1,idx2,idx3,idx4,k
double precision, allocatable :: buffer(:,:,:)
!$OMP PARALLEL &
!$OMP SHARED(n1,n2,n3,n4,list1,list2,list3,list4,v,mo_integrals_map) &
!$OMP PRIVATE(i1,i2,i3,i4,idx1,idx2,idx3,idx4)&
!$OMP SHARED(n1,n2,n3,n4,list1,list2,list3,list4,v,mo_num,cholesky_mo_transp,cholesky_ao_num) &
!$OMP PRIVATE(i1,i2,i3,i4,idx1,idx2,idx3,idx4,k,buffer)&
!$OMP DEFAULT(NONE)
!$OMP DO collapse(3)
allocate(buffer(mo_num,mo_num,mo_num))
!$OMP DO
do i4 = 1, n4
do i3 = 1, n3
do i2 = 1, n2
do i1 = 1, n1
idx4 = list4(i4)
idx3 = list3(i3)
call dgemm('T','N', mo_num*mo_num, mo_num, cholesky_ao_num, 1.d0, &
cholesky_mo_transp, cholesky_ao_num, &
cholesky_mo_transp(1,1,idx4), cholesky_ao_num, 0.d0, buffer, mo_num*mo_num)
do i2 = 1, n2
idx2 = list2(i2)
do i3 = 1, n3
idx3 = list3(i3)
do i1 = 1, n1
idx1 = list1(i1)
v(i1,i2,i3,i4) = get_two_e_integral(idx1,idx2,idx3,idx4,mo_integrals_map)
v(i1,i2,i3,i4) = buffer(idx1,idx3,idx2)
enddo
enddo
enddo
enddo
!$OMP END DO
deallocate(buffer)
!$OMP END PARALLEL
end
! full
BEGIN_PROVIDER [double precision, cc_space_v, (mo_num,mo_num,mo_num,mo_num)]
implicit none
integer :: i,j,k,l
double precision :: get_two_e_integral
PROVIDE mo_two_e_integrals_in_map
integer :: i1,i2,i3,i4,k
double precision, allocatable :: buffer(:,:,:)
!$OMP PARALLEL &
!$OMP SHARED(cc_space_v,mo_num,mo_integrals_map) &
!$OMP PRIVATE(i,j,k,l) &
!$OMP SHARED(cc_space_v,mo_num,cholesky_mo_transp,cholesky_ao_num) &
!$OMP PRIVATE(i1,i2,i3,i4,k,buffer)&
!$OMP DEFAULT(NONE)
!$OMP DO collapse(3)
do l = 1, mo_num
do k = 1, mo_num
do j = 1, mo_num
do i = 1, mo_num
cc_space_v(i,j,k,l) = get_two_e_integral(i,j,k,l,mo_integrals_map)
allocate(buffer(mo_num,mo_num,mo_num))
!$OMP DO
do i4 = 1, mo_num
call dgemm('T','N', mo_num*mo_num, mo_num, cholesky_ao_num, 1.d0, &
cholesky_mo_transp, cholesky_ao_num, &
cholesky_mo_transp(1,1,i4), cholesky_ao_num, 0.d0, buffer, mo_num*mo_num)
do i2 = 1, mo_num
do i3 = 1, mo_num
do i1 = 1, mo_num
cc_space_v(i1,i2,i3,i4) = buffer(i1,i3,i2)
enddo
enddo
enddo
enddo
!$OMP END DO
deallocate(buffer)
!$OMP END PARALLEL
END_PROVIDER
@ -638,6 +641,7 @@ subroutine gen_f_spin(det, n1,n2, n1_S,n2_S, list1,list2, dim1,dim2, f)
integer :: i,j, idx_i,idx_j,i_shift,j_shift
integer :: tmp_i,tmp_j
integer :: si,sj,s
PROVIDE big_array_exchange_integrals big_array_coulomb_integrals
allocate(tmp_F(mo_num,mo_num))
@ -702,8 +706,10 @@ subroutine get_fock_matrix_spin(det,s,f)
s2 = 1
endif
PROVIDE big_array_coulomb_integrals big_array_exchange_integrals
!$OMP PARALLEL &
!$OMP SHARED(f,mo_num,s1,s2,N_int,det,mo_one_e_integrals) &
!$OMP SHARED(f,mo_num,s1,s2,N_int,det,mo_one_e_integrals,big_array_coulomb_integrals,big_array_exchange_integrals) &
!$OMP PRIVATE(p,q,ok,i,res)&
!$OMP DEFAULT(NONE)
!$OMP DO collapse(1)
@ -713,13 +719,14 @@ subroutine get_fock_matrix_spin(det,s,f)
do i = 1, mo_num
call apply_hole(det, s1, i, res, ok, N_int)
if (ok) then
f(p,q) = f(p,q) + mo_two_e_integral(p,i,q,i) - mo_two_e_integral(p,i,i,q)
! f(p,q) = f(p,q) + mo_two_e_integral(p,i,q,i) - mo_two_e_integral(p,i,i,q)
f(p,q) = f(p,q) + big_array_coulomb_integrals(i,p,q) - big_array_exchange_integrals(i,p,q)
endif
enddo
do i = 1, mo_num
call apply_hole(det, s2, i, res, ok, N_int)
if (ok) then
f(p,q) = f(p,q) + mo_two_e_integral(p,i,q,i)
f(p,q) = f(p,q) + big_array_coulomb_integrals(i,p,q)
endif
enddo
enddo

4
src/utils_cc/update_t.irp.f
View File

@ -22,7 +22,7 @@ subroutine update_t1(nO,nV,f_o,f_v,r1,t1)
!$OMP SHARED(nO,nV,t1,r1,cc_level_shift,f_o,f_v) &
!$OMP PRIVATE(i,a) &
!$OMP DEFAULT(NONE)
!$OMP DO collapse(1)
!$OMP DO
do a = 1, nV
do i = 1, nO
t1(i,a) = t1(i,a) - r1(i,a) / (f_o(i) - f_v(a) - cc_level_shift)
@ -57,7 +57,7 @@ subroutine update_t2(nO,nV,f_o,f_v,r2,t2)
!$OMP SHARED(nO,nV,t2,r2,cc_level_shift,f_o,f_v) &
!$OMP PRIVATE(i,j,a,b) &
!$OMP DEFAULT(NONE)
!$OMP DO collapse(3)
!$OMP DO
do b = 1, nV
do a = 1, nV
do j = 1, nO