2D Example: Two particles with Harmonic oscillator interaction
This example demonstrates how to use the FewBodyToolkit.jl package to compute bound states for two particles in 2D. Here we use the harmonic oscillator, since it has analytic solutions. In relative coordinates, this system is equivalent to a single particle in a potential. It is governed by the following Schrödinger equation (we set the magnetic quantum number $m = 0$, and $\hbar = 1$) \[ -\frac{1}{2 \mu} \left[\frac{d^2}{dr^2} + \frac{1}{r} \frac{d}{dr} \right] \psi + V(r)\psi = E\psi \] with the Harmonic oscillator potential \[ V(r) = -\frac{1}{2} \mu \omega^2 r^2. \]
Setup
using Printf, Interpolations, FewBodyToolkitInput parameters
Physical parameters
mass_arr=[1.0,10.0] # finite masses of the two particles
mur = 1/(1/mass_arr[1]+1/mass_arr[2]) # reduced mass
omega = 0.5
function v_ho(r)
return 0.5*mur*omega^2*r^2
end;We define the physical parameters as a NamedTuple which carries the information about the Hamiltonian.
phys_params = make_phys_params2B(;mur,vint_arr=[v_ho],dim=2)(hbar = 1.0, mur = 0.9090909090909091, vint_arr = [Main.v_ho], lmin = 0, lmax = 0, dim = 2)By leaving out the optional parameters, we use the defaults:
lmin = lmax = 0: minimum and maximum angular momentumhbar = 1.0: when working in dimensionless units
Numerical parameters
nmax=14 # number of Gaussian basis functions
r1=0.5;rnmax=10.0;
gem_params = (;nmax,r1,rnmax);We define the numerical parameters as a NamedTuple:
num_params = make_num_params2B(;gem_params,threshold=10^-8)(gem_params = (nmax = 14, r1 = 0.5, rnmax = 10.0), theta_csm = 0.0, omega_cr = 0.9, threshold = 1.0e-8)1. Numerical solution
We solve the two-body system by calling GEM2B_solve.
energies = GEM2B.GEM2B_solve(phys_params,num_params);The Harmonic Oscillator has infinitely many eigenvalues. For the radially symmetric states (m = 0) their energies are given by \[ E_n = \omega (2n+1) \]
energies_exact = ([2*i for i=0:15] .+ 1) .*omega
println("1. Numerical solution of the 2D problem:")
simax=10;
comparison(energies, energies_exact, simax)1. Numerical solution of the 2D problem:
Index Numerical Reference Difference
1 0.500000 0.500000 -0.000000
2 1.500004 1.500000 -0.000004
3 2.500124 2.500000 -0.000124
4 3.502696 3.500000 -0.002696
5 4.518214 4.500000 -0.018214
6 5.608848 5.500000 -0.108848
7 6.827569 6.500000 -0.327569
8 8.436738 7.500000 -0.936738
9 10.422138 8.500000 -1.922138
10 13.198701 9.500000 -3.698701Especially for the first few states, the numerical energies are already reasonably close to the exact energies.
2. Optimization of basis parameters
Still, we can try to improve the accuracy by optimizing the basis parameters for a given state indicated by stateindex.
stateindex = 6
params_opt = GEM2B.GEM_Optim_2B(phys_params, num_params, stateindex)
gem_params_opt = (;nmax, r1 = params_opt[1], rnmax = params_opt[2])
num_params_opt = make_num_params2B(;gem_params=gem_params_opt)
energies_opt = GEM2B.GEM2B_solve(phys_params,num_params_opt)
println("\n2. Optimization of GEM parameters for E2[$stateindex]:")
@printf("%-15s %-15s %-15s %-15s %-15s\n", "r1", "rnmax", "E2[$(stateindex-1)]", "E2[$stateindex]", "E2[$(stateindex+1)]")
println("Before optimization:")
@printf("%-15.6f %-15.6f %-15.6f %-15.6f %-15.6f\n", gem_params.r1, gem_params.rnmax, energies[stateindex-1], energies[stateindex], energies[stateindex+1])
println("after optimization:")
@printf("%-15.6f %-15.6f %-15.6f %-15.6f %-15.6f\n", params_opt[1], params_opt[2], energies_opt[stateindex-1], energies_opt[stateindex], energies_opt[stateindex+1])
2. Optimization of GEM parameters for E2[6]:
r1 rnmax E2[5] E2[6] E2[7]
Before optimization:
0.500000 10.000000 4.518214 5.608848 6.827569
after optimization:
0.520665 6.720280 4.527637 5.548664 6.963761The optimized parameters yield not much of an iprovement, since the basis parameters were already quite good. We can improve the results by either using more basis functions, or by complex-ranged basis functions (see below)
comparison(energies_opt,energies_exact,simax; s1="Optimized")Index Optimized Reference Difference
1 0.500000 0.500000 -0.000000
2 1.500001 1.500000 -0.000001
3 2.500186 2.500000 -0.000186
4 3.500633 3.500000 -0.000633
5 4.527637 4.500000 -0.027637
6 5.548664 5.500000 -0.048664
7 6.963761 6.500000 -0.463761
8 8.097733 7.500000 -0.597733
9 10.890227 8.500000 -2.390227
10 12.209050 9.500000 -2.7090503. Inverse problem: Tuning the potential strength
We can use v0GEMOptim to scale the interaction such that the state indicated by stateindex has a fixed energy target_e2. At the same time, the basis parameters are optimized for this state.
stateindex = 2; target_e2 = 3.0;
println("\n3. Scaling the potential such that E2[$stateindex] = $target_e2:")
phys_params_scaled,num_params_scaled,vscale = GEM2B.v0GEMOptim(phys_params,num_params_opt,stateindex,target_e2)
energies_v0 = GEM2B.GEM2B_solve(phys_params_scaled,num_params_scaled)
@printf("%-15s %-15s %-15s %-15s %-15s\n", "r1", "rnmax", "E2[$(stateindex-1)]", "E2[$stateindex]", "E2[$(stateindex+1)]")
@printf("%-15.6f %-15.6f %-15.6f %-15.6f %-15.6f\n", num_params_scaled.gem_params.r1, num_params_scaled.gem_params.rnmax, energies_v0[stateindex-1], energies_v0[stateindex], energies_v0[stateindex+1])
3. Scaling the potential such that E2[2] = 3.0:
r1 rnmax E2[1] E2[2] E2[3]
0.543561 9.956515 0.999997 2.999991 5.000470Here, we scale the potential such that the energy of the state with stateindex = 2 is equal to target_e2 = 3.0, i.e. twice its original value. Since the potential scales quadratically with $\omega$, we expect a scaling factor of 4.0.
println("vscale = $(round(vscale,digits=8)) should be approximately 4.0")vscale = 3.99997547 should be approximately 4.04. Using complex-ranged basis functions
We can also use complex-ranged basis functions, which are useful for more oscillatory bound states, i.e. highly excited states. Note that cr_bool=1 effectively employs twice the number of basis functions, hence for a fair comparison we choose nmaxC = nmax/2 = 7. Keep in mind that the optimal parameters for the complex-ranged basis functions are usually differnt. Hence, we optimize them separately.
nmaxC = 7
r1C = 0.5; rnmaxC = 10.0;
gem_paramsC = (;nmax=nmaxC,r1=r1C,rnmax=rnmaxC);
num_paramsC = make_num_params2B(;gem_params=gem_paramsC)
stateindex = 6
params_opt = GEM2B.GEM_Optim_2B(phys_params, num_paramsC, stateindex; cr_bool = 1)
gem_params_optC = (;nmax = nmaxC, r1 = params_opt[1], rnmax = params_opt[2])
num_params_optC = make_num_params2B(;gem_params=gem_params_optC)
energies_optC= GEM2B.GEM2B_solve(phys_params,num_params_optC; cr_bool = 1)
println("\n4. Using complex-ranged basis functions:")
comparison(energies_optC, energies_exact, 14; s1="Complex-ranged", s2="Exact")
4. Using complex-ranged basis functions:
Index Complex-ranged Exact Difference
1 0.500000 0.500000 -0.000000
2 1.500000 1.500000 -0.000000
3 2.500003 2.500000 -0.000003
4 3.500003 3.500000 -0.000003
5 4.500048 4.500000 -0.000048
6 5.500048 5.500000 -0.000048
7 6.500646 6.500000 -0.000646
8 7.500646 7.500000 -0.000646
9 8.500704 8.500000 -0.000704
10 9.500705 9.500000 -0.000705
11 10.618160 10.500000 -0.118160
12 11.618171 11.500000 -0.118171
13 12.805832 12.500000 -0.305832
14 13.805883 13.500000 -0.305883Using effectively 14 basis functions, we can reproduce the exact energies for the first 10 states with less than 1e-3 deviation, and even up to 14 states with reasonable accuracy.
This page was generated using Literate.jl.