API

CentralPotential(f::Function)

A concrete implementation of PotentialFunction that represents a central potential. It takes a function f that defines the potential as a function of the radial distance r.

ContactPotential1D(v0::Float64, z0::Float64)

A concrete implementation of PotentialFunction that represents a 1D Contact (Dirac) potential:

\[V(z) = v_0 \delta(z - z_0)\]

where z is the 1D coordinate.

Arguments:

• v0::Float64: The strength of the potential.
• z0::Float64: The position of the delta function.

function (gp::ContactPotential1D)(z::Float64)

Evaluates the contact potential at a given 1D coordinate z.

GaussianPotential(v0::Float64, mu_g::Float64)

A concrete implementation of PotentialFunction that represents a Gaussian potential:

\[V(r) = v_0 e^{-\mu_g r^2}\]

where r is the radial distance.

Arguments:

• v0::Float64: The strength of the potential.
• mu_g::Float64: The width parameter of the Gaussian potential.

function (gp::GaussianPotential)(r::Float64)

Evaluates the Gaussian potential at a given radial distance r.

PotentialFunction

Abstract type for potential functions used in few-body calculations. This type serves as a parent type for more specific potential implementations

PreallocStruct2B{TTV, TS, TE}

A structure that holds preallocated arrays for the GEM-2B solver, used to store basis parameters, matrices, and results for two-body calculations.

Fields

• nu_arr::Vector{TTV}: Array of nonlinear variational parameters (basis exponents).
• S::Matrix{TS}: Overlap matrix between basis functions.
• T::Matrix{TTV}: Kinetic energy matrix.
• V::Matrix{TTV}: Potential energy matrix.
• energies::Vector{TE}: Array to store computed eigenvalues (energies).
• wavefunctions::Matrix{TTV}: Matrix to store eigenvectors (wavefunctions).

Type Parameters

• TTV: Element type for kinetic, potential, and basis parameter arrays (e.g., Float64 or ComplexF64).
• TS: Element type for the overlap matrix (e.g., Float64 or ComplexF64).
• TE: Element type for the energies array (e.g., Float64 or ComplexF64).

Description

This struct is designed to minimize memory allocations and improve performance by reusing arrays during repeated GEM-2B calculations. The types and sizes of the arrays are determined by the numerical parameters and whether complex rotation or complex scaling is used.

Keyword arguments

• num_params: A named tuple containing numerical parameters, including the maximum number of basis functions (nmax).
• cr_bool: Boolean indicating whether complex range basis functions are used (1 for true, 0 for false).
• csm_bool: Boolean indicating whether complex scaling is used (1 for true, 0 for false).

Example

PreallocStruct2B(num_params, cr_bool=0, csm_bool=0) # for real basis functions and no complex scaling
PreallocStruct2B(num_params, cr_bool=1, csm_bool=0) # for complex range basis functions and no complex scaling
PreallocStruct2B(num_params, cr_bool=0, csm_bool=1) # for real basis functions with complex scaling

GEM2B_solve(phys_params, num_params; wf_bool=0, cr_bool=0, csm_bool=0)

Solves two-body quantum mechanical problems using the Gaussian Expansion Method (GEM).

Arguments

• phys_params: Physical parameters describing the two-body system:
  • hbar::Float64: reduced Planck constant
  • mur::Float64: reduced mass
  • vint_arr=Vector{Any}: a vector of interactions
  • lmax::Int: power of r^lmax in the basis functions; indicator for the angular momentum in 3D
  • dim::Int: the spatial dimension
• num_params: Numerical parameters struct containing information on the set of basis functions:
  • gem_params::NamedTuple: (number of basis functions, smallest and largest range parameters).
  • theta_csm::Float64: Complex scaling angle (in radians) for the Complex Scaling Method.
  • omega_cr::Float64: Parameter controlling the frequency for complex-ranged basis functions.
  • threshold::Float64: Numerical threshold generalized eigenvalue solver.

Keywords

• wf_bool=0: Whether to return wavefunctions (0: energies only, 1: energies and wavefunctions)
• cr_bool=0: Whether to use complex rotation method (0: no, 1: yes)
• csm_bool=0: Whether to use complex scaling method (0: no, 1: yes)

Returns

• energies: Array of energy eigenvalues
• wavefunctions: (Optional) Array of eigenvectors if wf_bool=1

Example

phys_params = make_phys_params2B(hbar=1.0, mur=1.0, vint_arr=[GaussianPotential(-1.0, 0.5)], lmax=0, dim=3)
num_params = make_num_params2B(gem_params=(nmax=5, r1=1.0, rnmax=10.0), omega_cr=0.5, theta_csm=0.0, threshold=1e-10)
energies = GEM2B_solve(phys_params, num_params; wf_bool=0, cr_bool=0, csm_bool=0)
# or with wavefunctions:
energies, wavefunctions = GEM2B_solve(phys_params, num_params; wf_bool=1)
# Note: The function can handle 1D, 2D, or 3D problems based on the `dim` parameter in `phys_params`.

GEM_Optim_2B(phys_params, num_params, stateindex; cr_bool=0, g_tol=1e-9)

Optimize the ranges used in the Gaussian Expansion Method (GEM) for a specific 2-body state.

Arguments

• phys_params::NamedTuple: Physical parameters such as hbar, mur, vint_arr, lmax, lmin, and dim.
• num_params::NamedTuple: Numerical parameters such as gem_params, omega_cr, theta_csm, and threshold.
• stateindex::Int: An integer specifying the index of the state to optimize (e.g., 1 for the ground state).

Keyword Arguments

• cr_bool::Int=0: Indicates whether to use complex rotation.
• g_tol::Float64=1e-9: A value specifying the tolerance for the optimizer.

Returns

• A vector containing:
  • Optimized GEM parameters: [r1, rnmax].
  • The energy value of the optimized state.

Example

r1opt,rnmaxopt,energy = GEM_Optim_2B(phys_params, num_params, 1) # optimize for the ground state
r1opt,rnmaxopt,energy = GEM_Optim_2B(phys_params, num_params, 3; cr_bool = 1) # optimize for the third state (2nd excited) using complex-ranged basis functions

make_num_params2B(; gem_params=(nmax=5, r1=1.0, rnmax=10.0), theta_csm=0.0, omega_cr=0.9, threshold=1e-8)

Create and return a named tuple containing the numerical parameters for a two-body GEM calculation.

Keyword arguments

• gem_params::NamedTuple = (nmax=5, r1=1.0, rnmax=10.0): Parameters for the Gaussian Expansion Method (number of basis functions, smallest and largest range parameters).
• theta_csm::Float64 = 0.0: Complex scaling angle (in radians) for the Complex Scaling Method.
• omega_cr::Float64 = 0.9: Parameter controlling the frequency for complex-ranged basis functions.
• threshold::Float64 = 1e-8: Numerical threshold generalized eigenvalue solver.

Returns

• NamedTuple: Named tuple with the specified numerical parameters.

Example

make_num_params2B(gem_params=(nmax=20, r1=0.1, rnmax=50.0)) # for a larger basis set
make_num_params2B(gem_params=(nmax=20, r1=0.1, rnmax=50.0), theta_csm = 10.0) # non-zero rotation angle for complex scaling method (CSM)

make_phys_params2B(; hbar=1.0, mur=1.0, vint_arr=[[GaussianPotential(-1.0, 1.0)]], lmin=0, lmax=0, dim=3)

Create and return a named tuple containing the physical parameters for a two-body system.

Keyword arguments

• hbar::Float64 = 1.0: Reduced Planck constant used in calculations.
• mur::Float64 = 1.0: Reduced mass of the two-body system.
• vint_arr::Vector{Any} = [[GaussianPotential(-1.0, 1.0)]]: Array of interaction potentials or related parameters.
• lmin::Int = 0: Minimum orbital angular momentum quantum number.
• lmax::Int = 0: Maximum orbital angular momentum quantum number.
• dim::Int = 3: Spatial dimension of the system.

Returns

• NamedTuple: Named tuple with the specified physical parameters.

Example

make_phys_params2B(vint_arr=[GaussianPotential(-1.0, 0.5)], dim=1) # 1D system with Gaussian potential
make_phys_params2B(mur=0.5, vint_arr=[r -> -1/r], lmax=2)       # 3D Coulomb potential in d-wave (l=2) with reduced mass 0.5

v0GEMOptim(phys_params, num_params, stateindex, target_e2; cr_bool=0, rtol=1e-4, atol=10*eps(), g_tol=1e-9, output=false)

Finds a value v0crit to globally scale the potential in order to achieve the target energy target_e2 for the state specified by stateindex. Additionally, performs intermediate optimization of Gaussian Expansion Method (GEM) parameters.

Arguments

• phys_params::NamedTuple: Physical parameters including hbar, mur, vint_arr, lmax, lmin, and dim.
• num_params::NamedTuple: Numerical parameters including gem_params, omega_cr, theta_csm, and threshold.
• stateindex::Int: Index of the state to optimize, where 1 indicates the ground state.
• target_e2::Float64: Target energy value.

Keyword Arguments

• cr_bool::Int: Use complex rotation (default: 0).
• rtol::Float64: Relative tolerance for energy convergence (default: 1e-4).
• atol::Float64: Absolute tolerance for energy convergence (default: 10*eps()).
• g_tol::Float64: Tolerance for the optimizer (default: 1e-9).
• output::Bool: If true, prints intermediate optimization results (default: false).

Returns

• phys_params::NamedTuple: Updated physical parameters with optimized potential.
• num_params::NamedTuple: Updated numerical parameters with optimized GEM ranges.
• v0crit::Float64: The value to scale the potential with, to achieve the target energy. This is not the overall value of v0 but rather a scaling factor for the potential.

Example

phys_params_scaled,num_params_optimized,scalingfactor = v0GEMOptim(phys_params, num_params, 1, -2.5)

wavefun_arr(r_arr, phys_params, num_params, wf_arr; cr_bool=0)

Calculates the 2-body wavefunction at specified positions based on input parameters and coefficients.

Arguments

• r_arr::Vector{Float64}: Array of positions where the wavefunction is evaluated.
• phys_params::NamedTuple: Physical parameters, here we only need lmaxand dim
• num_params::NamedTuple: Numerical parameters containing the information about the Gaussian ranges
• wf_arr::Vector: Eigenvector from the diagonalization routine. Can be complex-valued.

Keyword Arguments

• cr_bool::Int=0: Determines whether to use complex-ranged Gaussians (1) or not (0). Defaults to 0.

Returns

• Vector{Float64}: The wavefunction values at the specified positions r_arr.

Example

wavefun_arr(0.0:0.1:10.0, phys_params, num_params, wf_arr; cr_bool=0)

wavefun_point(r, nu_arr, wf_arr, ll, dim)

Compute the value of the 2-body wavefunction at a given position r.

Arguments

• r::Float64: The radial position where the wavefunction is evaluated.
• nu_arr::Vector: Array of Gaussian widths.
• wf_arr::Vector: Coefficients of the wavefunction for the given state.
• ll::Int: Orbital angular momentum quantum number.
• dim::Int: Dimensionality of the system.

Returns

• Float64: The computed value of the wavefunction at the specified position r.

Example

nu_arr = [1.0, 0.0025]
wf_arr = [0.8,0.6]
psi_value = wavefun_point(1.0, nu_arr, wf_arr, 0, 3)

GEM3B1D_solve(phys_params, num_params; wf_bool=0, csm_bool=0, observ_params=(;stateindices=[],centobs_arr=[[],[],[]],R2_arr=[0,0,0]))

Solves the 1D three-body problem using the Gaussian Expansion Method (GEM).

Arguments

• phys_params: Physical parameters for the three-body system (e.g., masses, interaction potentials, etc.).
• num_params: Numerical parameters for the GEM calculation (e.g., basis size, grid parameters, etc.).
• wf_bool: (optional) If 1, also returns wavefunction-related observables. Default is 0.
• csm_bool: (optional) If 1, uses complex scaling method. Default is 0.
• observ_params: (optional) Parameters for observable calculations. Currently not supported for 1D.

Returns

• If wf_bool == 0: Returns an array of computed energies.
• If wf_bool == 1: Returns a tuple (energies, wavefunctions).

Example

phys_params = make_phys_params3B1D()
num_params = make_num_params3B1D()
energies = GEM3B3D_solve(phys_params, num_params) #solving with default parameters: three particles with the same mass and gaussian interaction

make_num_params3B1D(;lmin=0, Lmin=0, lmax=0, Lmax=0, gem_params=(nmax=5, r1=1.0, rnmax=10.0, Nmax=5, R1=1.0, RNmax=10.0), theta_csm=0.0, omega_cr=0.9, kmax_interpol=1000, threshold=10^-8)

Create and return a named tuple containing the numerical parameters for a three-body GEM calculation in 1D.

Keyword arguments

• lmin::Int = 0: Minimum power r^l used in the basis functions of the $r$ Jacobi coordinate.
• Lmin::Int = 0: Minimum power r^L used in the basis functions of the $R$ Jacobi coordinate.
• lmax::Int = 0: Maximum power r^l used in the basis functions of the $r$ Jacobi coordinate.
• Lmax::Int = 0: Maximum power r^L used in the basis functions of the $R$ Jacobi coordinate.
• gem_params::NamedTuple = (nmax=5, r1=1.0, rnmax=10.0, Nmax=5, R1=1.0, RNmax=10.0): Parameters for the Gaussian Expansion Method (number of basis functions, smallest and largest range parameters for both Jacobi coordinates).
• theta_csm::Float64 = 0.0: Complex scaling angle (in degrees) for the Complex Scaling Method.
• omega_cr::Float64 = 0.9: Parameter controlling the frequency for complex-ranged basis functions. Currently unsupported.
• kmax_interpol::Int = 1000: Number of numerical integration with effective Gaussian ranges used for interpolation.
• threshold::Float64 = 1e-8: Numerical threshold for the generalized eigenvalue solver.

Returns

• NamedTuple: Named tuple with the specified numerical parameters.

Example

make_num_params3B1D() # for the default basis set
make_num_params3B1D(gem_params=(nmax=10, r1=0.5, rnmax=20.0, Nmax=15, R1=1.0, RNmax=100.0), theta_csm=10.0) for a larger basis set and non-zero rotation angle for complex scaling method

make_phys_params3B1D(; hbar=1.0, mass_arr=[1.0,1.0,1.0], svals=["x","y","z"], vint_arr=[[GaussianPotential(-1.0, 1.0)],[GaussianPotential(-1.0, 1.0)],[GaussianPotential(-1.0, 1.0)]], parity=+1)

Create and return a named tuple containing the physical parameters for a three-body system in 1D.

Keyword arguments

• hbar::Float64 = 1.0: Reduced Planck constant used in calculations.
• mass_arr::Vector{Float64} = [1.0,1.0,1.0]: 3-element vector containing the masses of the three particles.
• svals::Vector{String} = ["x","y","z"]: 3-element vector containing the types of particles, used for automatic (anti-)symmetrization. "b" ("f") for identical bosons (fermions). Either 0, 2, or 3 identical particles.
• vint_arr::Vector{Vector{Any}} = [[GaussianPotential(-1.0, 1.0)],[GaussianPotential(-1.0, 1.0)],[GaussianPotential(-1.0, 1.0)]]: Vector of Vector of interaction potentials for each pair of particles: [[v231,v232,...],[v311,v312,...],[v121,v122,...]].
• parity::Int = 1: Parity parity=(-1)^(l+L) of the wave function. Possible values: +1,-1 for positive/negative parity, 0 for parity-violating potentials.

Returns

• NamedTuple: Named tuple with the specified physical parameters.

Example

make_phys_params3B1D() # default: system of three different particles with the same mass and Gaussian interactions
make_phys_params3B1D(mass_arr=[1.0,10.0,20.0], svals=["i","j","k"], vint_arr=[[v23],[v31],[v12_1,v12_2]]) # system of three different particles with interactions v23 (between particles 2 and 3), v31 (between particles 3 and 1), and v12_1, v12_2 (between particles 1 and 2). The interactions need to be defined before this call.

ISGL_solve(phys_params, num_params; wf_bool=0, csm_bool=0, observ_params=(;stateindices=[],centobs_arr=[[],[],[]],R2_arr=[0,0,0]))

Solves the 3D three-body problem using the Gaussian Expansion Method (GEM).

Arguments

• phys_params: Physical parameters for the three-body system (e.g., masses, interaction potentials, etc.).
• num_params: Numerical parameters for the GEM calculation (e.g., basis size, grid parameters, etc.).
• wf_bool: (optional) If 1, also returns wavefunction-related observables. Default is 0.
• csm_bool: (optional) If 1, uses complex scaling method. Default is 0.
• observ_params: (optional) Parameters for observable calculations.
  • stateindices: Indices of states for which observables are calculated.
  • centobs_arr: Array of central (only dependent on $r$; must be defined as functions) observables, for each Jacobi set (similar to vint_arr in phys_params).
  • R2_arr: Array which indicates whether the observable $\langle R^2 \rangle$ should be calculated (1) for any of the three Jacobi sets, or not (0).

Returns

• If wf_bool == 0: Returns an array of computed energies.
• If wf_bool == 1: Returns a tuple (energies, wavefunctions, centobs_output, R2_output).
  • energies: Vector of computed energies.
  • wavefunctions: Matrix of eigenvectors (column-wise) which contain the coefficients of the basis functions.
  • centobs_output: Mean values of central observables for the specified states. The first dimension corresponds to the Jacobi sets, the second to the observables, and the third to the states.
  • R2_output: Mean squared radii for the R-coordinate. Rows corresond to the Jacobi sets, columns to the states.

Example

phys_params = make_phys_params3B3D()
num_params = make_num_params3B3D()
energies = ISGL_solve(phys_params, num_params) #solving with default parameters: three particles with the same mass and gaussian interaction

csmgaussopt(gaussopt, csm_bool, theta_csm)

If csm_bool == 1, returns a new Vector{Vector{Tuple{Float64,ComplexF64}}} where each mu has been scaled by exp(2im * theta_csm * pi/180). Otherwise returns the original gaussopt unchanged.

make_num_params3B3D(;lmin=0, Lmin=0, lmax=0, Lmax=0, gem_params=(nmax=10, r1=0.2, rnmax=10.0, Nmax=10, R1=0.2, RNmax=20.0), theta_csm=0.0, omega_cr=0.9, mu0=0.08, c_shoulder=1.6, kmax_interpol=1000, threshold=10^-8)

Create and return a named tuple containing the numerical parameters for a three-body GEM calculation in 3D.

Keyword arguments

• lmin::Int = 0: Minimum power r^l used in the basis functions of the $r$ Jacobi coordinate.
• Lmin::Int = 0: Minimum power r^L used in the basis functions of the $R$ Jacobi coordinate.
• lmax::Int = 0: Maximum power r^l used in the basis functions of the $r$ Jacobi coordinate.
• Lmax::Int = 0: Maximum power r^L used in the basis functions of the $R$ Jacobi coordinate.
• gem_params::NamedTuple = (nmax=5, r1=1.0, rnmax=10.0, Nmax=5, R1=1.0, RNmax=10.0): Parameters for the Gaussian Expansion Method (number of basis functions, smallest and largest range parameters for both Jacobi coordinates).
• theta_csm::Float64 = 0.0: Complex scaling angle (in degrees) for the Complex Scaling Method.
• omega_cr::Float64 = 0.9: Parameter controlling the frequency for complex-ranged basis functions. Currently unsupported.
• mu0::Float64 = 0.08: Parameter (prefactor) for the ISGL method.
• c_shoulder::Float64 = 1.6: Parameter (base) for the ISGL method.
• kmax_interpol::Int = 1000: Number of numerical integration with effective Gaussian ranges used for interpolation.
• threshold::Float64 = 1e-8: Numerical threshold for the generalized eigenvalue solver.

Returns

• NamedTuple: Named tuple with the specified numerical parameters.

Example

make_num_params3B3D() # for the default basis set
make_num_params3B3D(gem_params=(nmax=10, r1=0.5, rnmax=20.0, Nmax=15, R1=1.0, RNmax=100.0), kmax_interpol=5000) for a larger basis set, and a finer interpolation grid for matrix-element calculations of interactions with many features.

make_phys_params3B3D(; hbar=1.0, mass_arr=[1.0,1.0,1.0], svals=["x","y","z"], vint_arr=[[GaussianPotential(-1.0, 1.0)],[GaussianPotential(-1.0, 1.0)],[GaussianPotential(-1.0, 1.0)]], J_tot=0, parity=+1, spin_arr=[0,0,0])

Create and return a named tuple containing the physical parameters for a three-body system in 3D.

Keyword arguments

• hbar::Float64 = 1.0: Reduced Planck constant used in calculations.
• mass_arr::Vector{Float64} = [1.0,1.0,1.0]: 3-element vector containing the masses of the three particles.
• svals::Vector{String} = ["x","y","z"]: 3-element vector containing the types of particles, used for automatic (anti-)symmetrization. "b" ("f") for identical bosons (fermions). Either 0, 2, or 3 identical particles.
• vint_arr::Vector{Vector{Any}} = [[GaussianPotential(-1.0, 1.0)],[GaussianPotential(-1.0, 1.0)],[GaussianPotential(-1.0, 1.0)]]: Vector of Vector of interaction potentials for each pair of particles: [[v231,v232,...],[v311,v312,...],[v121,v122,...]].
• J_tot::Int = 0: Total angular momentum of the system. Can take on half-integer values when spins are involved.
• parity::Int = 1: Parity parity=(-1)^(l+L) of the wave function. Possible values: +1,-1 for positive/negative parity, 0 for parity-violating potentials.
• spin_arr::Vector{Float64} = [0,0,0]: 3-element vector containing the spins [z1,z2,z_3] of the three particles. Used for automatic (anti-)symmetrization.

Returns

• NamedTuple: Named tuple with the specified physical parameters.

Example

make_phys_params3B3D() # default: system of three different particles with the same mass and Gaussian interactions
make_phys_params3B3D(mass_arr=[10.0,10.0,20.0], svals=["i","j","k"], vint_arr=[[v23],[v31],[v12_1,v12_2]], J_tot=3/2, spin_arr=[1/2,1/2,1/2],parity=-1, lmax=1, Lmax=1) # system of three different particles with half-integer spin, and interactions v23 (between particles 2 and 3), v31 (between particles 3 and 1), and v12_1, v12_2 (between particles 1 and 2). The interactions need to be defined before this call.