From 9cf881c45297c78b282b14e8ff3da563b7644a3f Mon Sep 17 00:00:00 2001
From: Daniel Menendez Crespo <danielmail7@gmail.com>
Date: Sun, 27 Sep 2020 12:31:20 +0200
Subject: [PATCH] An example of Julia calling libcint.

---
 examples/julia_call.jl | 117 +++++++++++++++++++++++++++++++++++++++++
 1 file changed, 117 insertions(+)
 create mode 100644 examples/julia_call.jl

diff --git a/examples/julia_call.jl b/examples/julia_call.jl
new file mode 100644
index 00000000..2c2166f8
--- /dev/null
+++ b/examples/julia_call.jl
@@ -0,0 +1,117 @@
+
+# general contracted DZ basis [3s1p/2s1p] for H₂
+#     exponents    contract-coeff
+# S   6.0          0.7               0.4
+#     2.0          0.6               0.3
+#     0.8          0.5               0.2
+# P   0.9          1.
+
+function gto_norm(n, a)
+   # normalization factor of function rⁿ exp(-ar²)
+    s = 2^(2n+3) * factorial(n+1) * (2a)^(n+1.5) / (factorial(2n+2) * √π)
+    sqrt(s)
+end
+
+CHARGE_OF  = 1
+PTR_COORD  = 2
+NUC_MOD_OF = 3
+PTR_ZETA   = 4
+ATM_SLOTS  = 6
+
+ATOM_OF    = 1
+ANG_OF     = 2
+NPRIM_OF   = 3
+NCTR_OF    = 4
+KAPPA_OF   = 5
+PTR_EXP    = 6
+PTR_COEFF  = 7
+BAS_SLOTS  = 8
+
+ptr_env = 20
+env = fill(0.0, ptr_env)
+
+i = 0
+#                     CHARGE_OF,PTR_COORD
+atm =           Int32[1         ptr_env  0 0 0 0]
+ptr_env += 3
+atm = vcat(atm, Int32[1         ptr_env  0 0 0 0])
+#             x  y   z (Bohr)
+append!(env, [0, 0, -0.8])
+append!(env, [0, 0,  0.8])
+
+# basis for atom #0
+# 3s -> 2s
+append!(env, [6., 2., .8])
+append!(env, [.7, .6, .5, .4, .3, .2])
+#                     ATOM_OF, ANG_OF, NPRIM_OF, NCTR_OF, KAPPA_OF, PTR_EXP, PTR_COEFF
+bas =           Int32[0        0       3         2         0        ptr_env  ptr_env+3  0]
+ptr_env += 9
+bas = vcat(bas, Int32[0        1       1         1        0         ptr_env  ptr_env+1  0])
+append!(env, [.9])
+append!(env, [1.])
+ptr_env += 1
+
+# basis functions for atom #1, they are the same to thoes of atom #0
+bas = vcat(bas, bas)
+
+shls = Int32[0, 1]
+
+natm = length(atm)
+nbas = length(bas)
+
+@info "Input" atm bas #env
+
+bas = permutedims(bas, (2,1))   # Column (Julia) -> Row (C) order
+atm = permutedims(atm, (2,1))   # Column (Julia) -> Row (C) order
+
+# Libcint API
+const LIBCINT = joinpath(@__DIR__, "path/to/libcint.so/dylib/dll")
+
+function CINTcgtos_spheric(bas_id, bas)
+   @static if VERSION >= v"1.5.0" # clearer syntax
+      @ccall LIBCINT.CINTcgtos_spheric(bas_id::Cint, bas::Ptr{Cint})::Cint
+   else # works anyway
+      ccall( (:CINTcgtos_spheric, LIBCINT), Cint, (Cint, Ptr{Cint},), bas_id, bas)
+   end
+end
+
+function cint1e_ipnuc_sph(buf, shls, atm, natm, bas, nbas, env)
+   @static if VERSION >= v"1.5.0" # clearer syntax
+      @ccall LIBCINT.cint1e_ipnuc_sph(
+                                      buf  :: Ptr{Cdouble},
+                                      shls :: Ptr{Cint},
+                                      atm  :: Ptr{Cint},
+                                      natm :: Cint,
+                                      bas  :: Ptr{Cint},
+                                      nbas :: Cint,
+                                      env  :: Ptr{Cdouble}
+                                     )::Cvoid
+   else # works anyway
+      ccall( (:cint1e_ipnuc_sph, LIBCINT), Cvoid, (
+                                      Ptr{Cdouble},
+                                      Ptr{Cint},
+                                      Ptr{Cint},
+                                      Cint,
+                                      Ptr{Cint},
+                                      Cint,
+                                      Ptr{Cdouble}, ),
+                                      buf  ,
+                                      shls ,
+                                      atm  ,
+                                      natm ,
+                                      bas  ,
+                                      nbas ,
+                                      env  ,
+                                     )
+   end
+end
+
+di = CINTcgtos_spheric(0, bas) # C uses 0-indexing!
+dj = CINTcgtos_spheric(1, bas) # C uses 0-indexing!
+@assert di == (2bas[ANG_OF,1]+1) * bas[NCTR_OF,1]
+@assert dj == (2bas[ANG_OF,2]+1) * bas[NCTR_OF,2]
+
+buf = zeros(3, dj, di) # Column (Julia) -> Row (C) order
+
+cint1e_ipnuc_sph(buf, shls, atm, natm, bas, nbas, env)
+@info "∇⟨i|nuc|j⟩ = " permutedims(buf, (3,2,1))