Skip to content

Fix related to performing AD on molecular integrals #22

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
tansongchen wants to merge 1 commit into
base: main
Choose a base branch
from
Draft

Fix related to performing AD on molecular integrals #22

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
15 changes: 7 additions & 8 deletions src/Integrals/OneElectron.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -11,7 +11,7 @@ function nuclear!(out, BS::BasisSet{LCint}, i, j)
cint1e_nuc_sph!(out, [i,j], BS.lib)
end

# Backend: ACSint -- fall back
# Backend: ACSint -- fall back
function overlap!(out, BS::BasisSet, i, j)
generate_S_pair!(out, BS, i ,j)
end
Expand Down Expand Up @@ -62,11 +62,10 @@ function get_1e_matrix!(callback, out, BS::BasisSet, rank::Integer = 0)
# Offset list for each shell, used to map shell index to AO index
ao_offset = [sum(Nvals[1:(i-1)]) for i = 1:BS.nshells]

buf_arrays = [zeros(eltype(BS.atoms[1].xyz), 3^rank * Nmax^2) for _ = 1:Threads.nthreads()]
buf = zeros(eltype(BS.atoms[1].xyz), 3^rank * Nmax^2)

Threads.@threads :static for i in 1:BS.nshells
for i in 1:BS.nshells
Ni = Nvals[i]
buf = buf_arrays[Threads.threadid()]
ioff = ao_offset[i]
for j in i:BS.nshells
Nj = Nvals[j]
Expand All @@ -91,11 +90,11 @@ function get_1e_matrix!(callback, out, BS::BasisSet, rank::Integer = 0)
return out
end

###########################################################
###########################################################
###########################################################
###########################################################
# Mixed basis #
###########################################################
###########################################################
###########################################################
###########################################################

# Backend: ACSint -- fall back
function overlap!(out, BS1::BasisSet, BS2::BasisSet, i, j)
Expand Down
19 changes: 9 additions & 10 deletions src/Integrals/TwoElectronFourCenter.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -22,7 +22,7 @@ function sparseERI_2e4c(BS::BasisSet, cutoff = 1e-12)
# Pre allocate array to save ij pairs
ij_vals = Array{NTuple{2,Int32}}(undef, num_ij)

# Pre allocate array to save σij, that is the screening parameter for Schwarz
# Pre allocate array to save σij, that is the screening parameter for Schwarz
σvals = zeros(T, num_ij)

### Loop thorugh i and j such that i ≤ j. Save each pair into ij_vals and compute √σij for integral screening
Expand All @@ -42,14 +42,14 @@ function sparseERI_2e4c(BS::BasisSet, cutoff = 1e-12)
end

buf_arrays = [zeros(T, Nmax^4) for _ = 1:Threads.nthreads()]

# i,j,k,l => Shell indexes starting at zero
# I, J, K, L => AO indexes starting at one
@sync for ij in eachindex(ij_vals)
Threads.@spawn begin
@inbounds begin
buf = buf_arrays[Threads.threadid()]
i,j = ij_vals[ij]
i,j = ij_vals[ij]
Li, Lj = Nvals[i], Nvals[j]
Lij = Li*Lj
ioff = ao_offset[i]
Expand All @@ -59,7 +59,7 @@ function sparseERI_2e4c(BS::BasisSet, cutoff = 1e-12)
if σ < cutoff
continue
end
k,l = ij_vals[kl]
k,l = ij_vals[kl]
Lk, Ll = Nvals[k], Nvals[l]
Lijk = Lij*Lk
koff = ao_offset[k]
Expand All @@ -80,7 +80,7 @@ function sparseERI_2e4c(BS::BasisSet, cutoff = 1e-12)
L < K ? break : nothing

# L ≥ K
KL = (L * (L + 1)) >> 1 + K
KL = (L * (L + 1)) >> 1 + K
bkl = Lij*(ks-1) + bl
for js = 1:Lj
J = joff + js
Expand All @@ -91,7 +91,7 @@ function sparseERI_2e4c(BS::BasisSet, cutoff = 1e-12)

IJ = (J * (J + 1)) >> 1 + I

#KL < IJ ? continue : nothing # This restriction does not work... idk why
#KL < IJ ? continue : nothing # This restriction does not work... idk why

idx = index2(IJ,KL) + 1
out[idx] = buf[is + bjkl]
Expand Down Expand Up @@ -125,7 +125,7 @@ end

function ERI_2e4c(BS::BasisSet)
N = BS.nbas
out = zeros(N, N, N, N)
out = zeros(eltype(BS.atoms[1].xyz), N, N, N, N)
ERI_2e4c!(out, BS)
end

Expand Down Expand Up @@ -160,13 +160,12 @@ function ERI_2e4c!(out, BS::BasisSet)
end

# Initialize array for results
bufs = [zeros(Cdouble, Nmax^4) for _ = 1:Threads.nthreads()]
Threads.@threads :static for (i,j,k,l) in unique_idx
buf = zeros(eltype(out), Nmax^4)
for (i,j,k,l) in unique_idx
# Shift indexes (C starts with 0, Julia 1)
id, jd, kd, ld = i+1, j+1, k+1, l+1
Ni, Nj, Nk, Nl = Nvals[id], Nvals[jd], Nvals[kd], Nvals[ld]

buf = bufs[Threads.threadid()]
# Compute ERI
ERI_2e4c!(buf, BS, id, jd, kd, ld)

Expand Down