How to get dipole moments for a molecule in Atomic Orbital basis within the Julia ecosystem?

Till now, I best solution I found was to use PySCF through PyCall.jl package. After defining the PySCF molecule object through specifying the geometry and the basis set, the integrals in question can be extracted as the follows

julia> using PyCall

julia> gto = pyimport("pyscf.gto")
PyObject <module 'pyscf.gto' from '/home/anik/.local/lib/python3.10/site-packages/pyscf/gto/__init__.py'>

julia> scf = pyimport("pyscf.scf")
PyObject <module 'pyscf.scf' from '/home/anik/.local/lib/python3.10/site-packages/pyscf/scf/__init__.py'>

julia> mol = gto.Mole()
PyObject <pyscf.gto.mole.Mole object at 0xfeb583254b0>

julia> mol.atom = "O 0 0 0; H 0 1 0; H 0 0 1"
"O 0 0 0; H 0 1 0; H 0 0 1"

julia> mol.basis="sto-3g"
"sto-3g"

julia> mol.build()
PyObject <pyscf.gto.mole.Mole object at 0xfeb583254b0>

julia> mol.atom
"O 0 0 0; H 0 1 0; H 0 0 1"

julia> mol.basis
"sto-3g"

julia> hf = scf.RHF(mol)
PyObject <pyscf.scf.hf.RHF object at 0xfeb58907c40>

julia> hf.kernel()
converged SCF energy = -74.9611711378677
-74.96117113786772

julia> mol.set_common_origin([0,0,0])
PyObject <pyscf.gto.mole.Mole object at 0xfeb583254b0>

julia> dipole_ao = mol.intor('int1e_r')
ERROR: ParseError:
# Error @ REPL[13]:1:24
dipole_ao = mol.intor('int1e_r')
#                      └─────┘ ── character literal contains multiple characters
Stacktrace:
 [1] top-level scope
   @ none:1

julia> dipole_ao = mol.intor("int1e_r")
3×7×7 Array{Float64, 3}:
[:, :, 1] =
 0.0  0.0  0.0507919  0.0        0.0        0.0        0.0
 0.0  0.0  0.0        0.0507919  0.0        0.0035842  0.0
 0.0  0.0  0.0        0.0        0.0507919  0.0        0.0035842

[:, :, 2] =
 0.0  0.0  0.641173  0.0       0.0       0.0       0.0
 0.0  0.0  0.0       0.641173  0.0       0.357697  0.0
 0.0  0.0  0.0       0.0       0.641173  0.0       0.357697

[:, :, 3] =
 0.0507919  0.641173  0.0  0.0  0.0  0.287551  0.287551
 0.0        0.0       0.0  0.0  0.0  0.0       0.0
 0.0        0.0       0.0  0.0  0.0  0.0       0.0

This seems to give the output automatically as Julia’s tensor data type, so there seems to be no need to parse anything. If anyone has a better solution, please let me know.