I am trying to develop an implementation of a certain electronic structure method within Julia aimed at modelling strong light-matter interactions. The method is closely related to the widely popular Hartree-Fock method. For this, I would be needing a lot of basis-set related integrals to be computed, namely the one and two electron electron repulsion integrals (which I found can be easily gotten by specifying molecular geometry and basis set using Fermi.jl or GaussianBasis.jl) and dipole moment integrals. What I cannot understand is that how to retrieve are the dipole moment integrals in the atomic orbital (AO) basis, this is essential for the method in question to work.
Given the time constraints, obviously it is not a good idea to write the handling of basis sets and integral libraries from scratch. Is there any way I can get these integrals using some package or from somewhere else?
1 Like
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.