<jupyter_start><jupyter_text>Stabilization plot via DVR box-size<jupyter_code>import numpy as np
from numpy.linalg import eigvals
import matplotlib
import matplotlib.pyplot as plt
matplotlib.use('Qt5Agg')
%matplotlib qt5
import pandas as pd
#
#  extend path by location of the dvr package
#
import sys
sys.path.append('../../Python_libs')
#sys.path.append('/home/thomas/Current_Work/Jolanta-by-dvr/Python_libs')
import dvr
from jolanta import Jolanta_3D
amu_to_au=1822.888486192
au2cm=219474.63068
au2eV=27.211386027
Angs2Bohr=1.8897259886
#
# Jolanata parameters a, b, c:
#
#    bound state:  -12.26336 eV
#    resonance: (3.279526396 - 0.2079713j)  eV
jparam=(0.028, 1.0, 0.028)
#
#  compute DVR of T and V
#  then show the density of states
#  in a potential + energy-levels plot
#
xmin=0
xmax=20      # grid from -xmax to xmax
thresh = 5   # maximum energy for plot
ppB = 15     # grid points per Bohr

nGrid=int((xmax-xmin)*ppB)
xs = dvr.DVRGrid(xmin, xmax, nGrid)
Ts = dvr.KineticEnergy(1, xmin, xmax, nGrid)
Vs=Jolanta_3D(xs, jparam)
[energy, wf] = dvr.DVRDiag2(nGrid, Ts, Vs, wf=True)

n_ene=0
for i in range(nGrid):
    print("%3d  %12.8f au = %12.5f eV" % (i+1, energy[i], energy[i]*au2eV))
    n_ene += 1
    if energy[i]*au2eV > thresh:
        break

# "DVR normalization", sum(wf[:,0]**2)
# this is correct for plotting

c=["orange", "blue"]
h=float(xmax) / (nGrid+1.0)
scale=400

plt.plot(xs,Vs*au2eV, '-', color="black")
for i in range(n_ene):
    plt.plot(xs, scale*wf[:,i]**2+energy[i]*au2eV, '-', color=c[i%len(c)])
#plt.ylim(0, 20)
plt.xlabel('$x$ [Bohr]')
plt.ylabel('$E$ [eV]')
plt.ylim(-8,8)
plt.show()<jupyter_output>1   -0.26351131 au =     -7.17051 eV
  2    0.06066341 au =      1.65074 eV
  3    0.11675239 au =      3.17699 eV
  4    0.17971292 au =      4.89024 eV
  5    0.29854286 au =      8.12377 eV<jupyter_text>Stabilization loop<jupyter_code>L_min=20        # smallest grid extent (box size)
L_max=45       # largest grid extent (box size)
nEs_keep=30    # how many energies are kept

s_min=1/L_max**2
s_max=1/L_min**2


n_s=121  
slist=np.linspace(s_min, s_max, num=n_s, endpoint=True)
Llist=1/np.sqrt(slist)

run_data = np.zeros((n_s, nEs_keep))  # array used to collect the energies
dc_data  = np.zeros((n_s, nEs_keep))  # array used to collect the DC energies


for i in range(n_s):
    xmax=1/np.sqrt(slist[i])
    xmin=0
    nGrid=int((xmax-xmin)*ppB)
    xs = dvr.DVRGrid(xmin, xmax, nGrid)
    Vs = Jolanta_3D(xs, jparam)
    Ts = dvr.KineticEnergy(1, xmin, xmax, nGrid)
    energy, wf= dvr.DVRDiag2(nGrid, Ts, Vs, wf=True)
    run_data[i,:] = energy[0:nEs_keep]        
    print(i+1, end=" ")
    if (i+1)%10==0:
        print()

run_data *= au2eV
dc_data *= au2eV
slist*=1000<jupyter_output>1 2 3 4 5 6 7 8 9 10 
11 12 13 14 15 16 17 18 19 20 
21 22 23 24 25 26 27 28 29 30 
31 32 33 34 35 36 37 38 39 40 
41 42 43 44 45 46 47 48 49 50 
51 52 53 54 55 56 57 58 59 60 
61 62 63 64 65 66 67 68 69 70 
71 72 73 74 75 76 77 78 79 80 
81 82 83 84 85 86 87 88 89 90 
91 92 93 94 95 96 97 98 99 100 
101 102 103 104 105 106 107 108 109 110 
111 112 113 114 115 116 117 118 119 120 
121<jupyter_text>$E(L)$ and $E(\alpha)$ plots<jupyter_code>plt.cla()
for i in range(0, nEs_keep):
    plt.plot(slist,run_data[:,i], 'o-', color='blue')
    if do_DC:
        plt.plot(slist,dc_data[:,i], '-', color='grey')
    plt.ylim(0,7)
plt.xscale('linear')
plt.show()
#
# copy the nEs lowest states from run_data to stab_data
#
nEs=6
stab_data = np.zeros((n_s,nEs+2))
i = 0
stab_data[:,0] = 1/np.sqrt(slist)
stab_data[:,1] = slist
for i in range(nEs):
    stab_data[:,i+2] = run_data[:,i]

# for checking what we got
plt.cla()
for i in range(nEs):
    plt.plot(stab_data[:,0],stab_data[:,i+2], '-', color='blue')
plt.ylim(0,10)
plt.show()
cols = ['L', 'z']
for i in range(nEs):
    cols.append('E'+str(i+1))
df = pd.DataFrame(stab_data, columns=cols)
df.head(5)
df.to_csv('DVR_stab_plot.csv', index=False)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1491670873296315, 'fraction_numerical': 0.12771327612317013, 'mean_word_length': 3.6405152224824358, 'pattern_counts': {'":': 0, '<': 11, '<?xml version=': 0, '>': 12, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '3', 'dir': 'jupyter-structured-clean-dedup', 'id': '247299', 'license': 'MIT', 'n_tokens_mistral': 2149, 'n_tokens_neox': 1746, 'n_words': 462, 'path': 'notebooks/Stab_3D/.ipynb_checkpoints/Stab-DVR-checkpoint.ipynb', 'repo_name': 'tsommerfeld/L2-methods_for_resonances', 'size': '3962'}