<jupyter_start><jupyter_text>For small $\theta$ everything is fine, but for larger angles basis set artifacts become apparent = there are eigenvalues off the rotation string with small real and hige imaginary part. Thus, above some $\theta$, sorting the eigenvalue by low real part stops working as a filter.<jupyter_code>#
#  complex scaling loop: 
#
#  start on the real axis (theta=0) and rotate to theta = theta_max 
#
#  we keep n_keep energies with the lowest real part 
#

n_keep=30
theta_min=0
theta_max=40
n_theta=theta_max - theta_min + 1
thetas=[t for t in range(theta_min, theta_max+1)]

run_data = np.zeros((n_theta,n_keep), complex)  # array used to collect all theta-run data
#run_data = np.zeros((n_keep,n_theta), complex)  # array used to collect all theta-run data

for i_theta in range(n_theta):
    theta=thetas[i_theta]/180.0*np.pi
    Vs = Jolanta_1Db(xs*np.exp(1j*complex(theta)), jparam)
    H_theta = np.exp(-2j*complex(theta)) * Ts + np.diag(Vs)
    energies = eigvals(H_theta)
    energies.sort()
    run_data[i_theta,:] = energies[0:n_keep]
    #run_data[:,i_theta] = energies[0:n_keep]
    print(i_theta+1, end=" ")
    if (i_theta+1)%10==0:
        print()

run_data *= au2eV
#
# useful piece of the complex plane 
# (if unknown, plot all and zoom with matplotlib)
#
plt.cla()
for i in range(0, n_keep):
    plt.plot(run_data[:,i].real,  run_data[:,i].imag, 'o')
plt.xlim(-1,10)
plt.ylim(-1,0)
plt.show()
#
#
#  two follow ideas:
#
#   - at the last theta compute the angles and compare with 2*theta
#     if significantly smaller, then resonance
#
#   - once a trajectory has five values, use them
#     to establish a Pade[2,2]
#     then predict for the next theta
#
def follow_nearest(follow, es):
    """
    follow the energy closet to e0 from the real axis into the complex plane
    es is a table of theta-run data es[i_theta,j_energies]
    the algorithm used is simply nearest to the old energy
    """
    (n_thetas, n_energies) = es.shape
    trajectory = np.zeros(n_thetas,complex)
    for j in range(0,n_thetas):
        i = np.argmin(abs(es[j,:]-follow))
        follow = trajectory[j] = es[j,i]
    return trajectory
n_save = n_keep//2
trajectories = np.zeros((n_theta, n_save), complex)
for j in range(n_save):
    trajectories[:,j] = follow_nearest(run_data[0,j], run_data)
for i in range(0, n_save):
    plt.plot(trajectories[:,i].real,  trajectories[:,i].imag, '-')
plt.show()
#
#  save n_save trajectories to file
#  csv as real and imag 
#  (at the moment easier than csv with complex)
#  also, include no header, because the energies need to be sorted
#  into trajectories first
#
#fname="complex_scaling_rmax."+str(int(rmax))+"_ppB."+str(ppB)+".csv"
#read_write.write_theta_run(fname,thetas,trajectories)
#
# regarding the question of automatization:
# the resonance clearly stands out in CS/DVR
# use to make a loop over the Jolanata parameters 
# and map Eres(a,b,c) in a useful range: 0.1 to 8 eV
#
for i in range(n_save):
    print(i, np.angle(trajectories[-1,i],deg=True))
res_traj = trajectories[:,3]
def naive_derivative(xs, ys):
    """ naive forward or backward derivative """
    return (ys[1]-ys[0])/(xs[1]-xs[0])

def central_derivative(xs, ys):
    """ central derivative at x[1] """
    return (ys[2]-ys[0])/(xs[2]-xs[0])

def five_point_derivative(xs, ys):
    """ five-point derivative at x[2] """
    """ (-ys[0] + 8*ys[1] - 8*ys[3] + ys[4])/(12*h)  """
    return (-ys[0] + 8*ys[1] - 8*ys[3] + ys[4])/(xs[4]-xs[0])/3
abs_der = np.zeros(n_theta)
abs_der[0] = abs(naive_derivative(thetas[0:2], res_traj[0:2]))
abs_der[1] = abs(central_derivative(thetas[0:3], res_traj[0:3]))
for k in range(2,n_theta-2):
    abs_der[k] = abs(five_point_derivative(thetas[k-2:k+3], res_traj[k-2:k+3]))
abs_der[-2] = abs(naive_derivative(thetas[-3:], res_traj[-3:]))
abs_der[-1] = abs(naive_derivative(thetas[-2:], res_traj[-2:]))
plt.cla()
plt.plot(thetas, np.log(abs_der), 'o-')
plt.show()
#
#  get a feeling for the stabilitity of the value
#
j_opt = np.argmin(abs_der)-2
print(j_opt, thetas[j_opt])
print(res_traj[j_opt])
print(res_traj[j_opt-1]-res_traj[j_opt])
print(res_traj[j_opt+1]-res_traj[j_opt])<jupyter_output>38 38
(0.8696689855466049-0.14183784022465826j)
(1.969203050256052e-05+5.046653604012796e-05j)
(-2.0663096704232053e-05-4.001912573572963e-05j)
{'fraction_non_alphanumeric': 0.15010423905489922, 'fraction_numerical': 0.047718322909427846, 'mean_word_length': 4.47275031685678, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '1', 'dir': 'jupyter-structured-clean-dedup', 'id': '642809', 'license': 'MIT', 'n_tokens_mistral': 1841, 'n_tokens_neox': 1720, 'n_words': 416, 'path': 'notebooks/CS/.ipynb_checkpoints/CS-J3D-checkpoint.ipynb', 'repo_name': 'tsommerfeld/L2-methods_for_resonances', 'size': '4317'}