<jupyter_start><jupyter_text>CS/GTO: step 1 and 2<jupyter_code>show_plots=False
import numpy as np
from scipy.linalg import eig, eigh, eigvals, eigvalsh
from scipy.optimize import minimize
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')
from captools import simple_traj_der, five_point_derivative 
from GTO_basis import GBR
from jolanta import Jolanta_3D
amu_to_au=1822.888486192
au2cm=219474.63068
au2eV=27.211386027
Angs2Bohr=1.8897259886
#
# Jolanata parameters a, b, c:
#
#   CS-DVR:   
#      bound state:  -7.17051 eV
#      resonance (3.1729556 - 0.16085j) eV
#
jparam=(0.028, 1.0, 0.028)<jupyter_output><empty_output><jupyter_text>* Create a valence set $[\alpha_0, \alpha_0/s, \alpha_0/s^2, ..., \alpha_N]$* Diagonalize **H** to compare with $E_0^{DVR}$* Add diffuse functions $[\alpha_N/s_{df}, ...]$* Diagonalize **H** again<jupyter_code>sets=['GTO_unc', 'GTO_DZ', 'GTO_TZ']
bas = sets[0]
nval=10
a0=17.0
s=2
ndf=4
sdf=1.5
if bas == 'GTO_unc':
    contract = (0,0)
elif bas == 'GTO_DZ':
    contract = (1,1)  # one contracted, one uncontracted function
elif bas == 'GTO_TZ':
    contract = (1,2)  # one contracted, two uncontracted function
else:
    print('No such basis.')

fname='Traj_' + bas + '.csv'
print(fname)<jupyter_output>Traj_GTO_unc.csv<jupyter_text>Valence setCompare the bound state with DVR: $E_0 = -7.17051$ eV<jupyter_code>alpha_val=[a0]
for i in range(nval-1):
    alpha_val.append(alpha_val[-1]/s)
Val = GBR(alpha_val, jparam, contract=contract, diffuse=(0,0))
S, T, V = Val.STV()
Es, cs = eigh(T+V, b=S)
print(f'E0 = {Es[0]*au2eV:.6f}   Emax = {Es[-1]*au2eV:.6f}')
Val.print_exp()

if show_plots:
    scale=10
    xmax=15
    xs=np.linspace(0.1,xmax,200)
    Vs=Jolanta_3D(xs, jparam)
    plt.cla()
    plt.plot(xs,Vs*au2eV, '-', color="blue")
    for i in range(len(Es)):
        ys=Val.eval_vector(cs[:,i], xs)
        plt.plot(xs,scale*ys**2+Es[i]*au2eV, '-')
    plt.ylim(-8,10)
    plt.show()<jupyter_output>E0 = -7.170439   Emax = 2419.320761
      alpha         r0=1/sqrt(alpha)     Norm
   1.70000000e+01    2.42535625e-01    1.0070e+02
   8.50000000e+00    3.42997170e-01    4.2341e+01
   4.25000000e+00    4.85071250e-01    1.7802e+01
   2.12500000e+00    6.85994341e-01    7.4848e+00
   1.06250000e+00    9.70142500e-01    3.1470e+00
   5.31250000e-01    1.37198868e+00    1.3231e+00
   2.65625000e-01    1.94028500e+00    5.5631e-01
   1.32812500e-01    2.74397736e+00    2.3390e-01
   6.64062500e-02    3.88057000e+00    9.8343e-02
   3.32031250e-02    5.48795472e+00    4.1348e-02<jupyter_text>Extend the basis by a diffuse set<jupyter_code>Bas = GBR(alpha_val, jparam, contract=contract, diffuse=(ndf,sdf))
S, T, V = Bas.STV()
Es, cs = eigh(T+V, b=S)
nEs = len(Es)
print(f'E0 = {Es[0]*au2eV:.6f}   Emax = {Es[-1]*au2eV:.6f}')

if show_plots:
    Emax=10 # eV
    plt.cla()
    plt.plot(xs,Vs*au2eV, '-', color="blue")
    for i, E in enumerate(Es):
        ys=Bas.eval_vector(cs[:,i], xs)
        plt.plot(xs,scale*ys**2+E*au2eV, '-')
        if E*au2eV > Emax:
            break

    plt.ylim(-10,Emax+1)
    plt.show()<jupyter_output>E0 = -7.170447   Emax = 2421.905980<jupyter_text>*** CSExample for testing library. With:`nval=10, a0=17.0, s=2, ndf=4, sdf=1.5, theta = 10 deg`1. `(-7.171756, 0.000583)`2. `(0.327329, -0.171381)`3. `(1.228884, -0.692745)`4. `(3.165013, -0.147167)`5. `(3.431676, -1.806220)`<jupyter_code>if True:
    theta=8.0/180.0*np.pi
    print("theta = %f" % (theta))
    H_theta = Bas.H_theta(theta+0.00005, 1.00000)
    energies = eigvals(H_theta, b=S)
    energies.sort()
    energies*=au2eV
    for e in energies:
        print("(%f, %f)" % (e.real, e.imag))<jupyter_output>theta = 0.139626
(-7.171084, 0.000342)
(0.344332, -0.140453)
(1.304015, -0.565986)
(3.152077, -0.151021)
(3.591431, -1.412265)
(7.487083, -2.093086)
(12.154227, -3.213948)
(21.924637, -7.105574)
(47.856322, -17.015877)
(108.475496, -36.285772)
(236.367357, -74.120776)
(500.570665, -150.587619)
(1056.707013, -310.511939)
(2326.966794, -675.120632)<jupyter_text>$\theta$-run<jupyter_code>n_keep=nEs

theta_min=0
theta_max=16
n_theta=101
thetas=np.linspace(theta_min, theta_max, num=n_theta)
run_data = np.zeros((n_theta,n_keep), complex)  # array used to collect all theta-run data

for i, tdeg in enumerate(thetas):
    theta = tdeg/180.0*np.pi
    H_theta = Bas.H_theta(theta, 1.0)
    energies = au2eV * eigvals(H_theta, b=S)
    energies.sort()
    run_data[i,:] = energies[0:n_keep]
    print(i+1, end=" ")
    if (i+1)%10==0:
        print()<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<jupyter_text>Raw $\eta$ trajectories<jupyter_code>plt.cla()
for i in range(0, n_keep):
    plt.plot(run_data[:,i].real,  run_data[:,i].imag, 'o')
plt.xlim(0,8)
plt.ylim(-2,0.5)
plt.show()<jupyter_output><empty_output><jupyter_text>Get the resonance trajectory by naive nearest followWe start with the energy nearest `follow`<jupyter_code>#follow=3.7
follow=3.0
es=np.zeros(n_theta,complex)

for j in range(0,n_theta):
    i = np.argmin(abs(run_data[j,:]-follow))
    es[j] = run_data[j,i]
    follow = es[j]
plt.cla()
plt.ylim(-1,0.0)
plt.plot(es.real, es.imag, 'o-')
plt.show()
abs_der = np.abs(simple_traj_der(thetas, es))
plt.cla()
plt.plot(thetas, np.log10(abs_der), 'o-')
plt.xlabel(r'$\theta$ [deg.]')
plt.ylabel(r'$\log \vert dE/d \theta{} \vert $')
plt.show()
df = pd.DataFrame({"theta": thetas, 
                   "ReE":es.real, 
                   "ImE":es.imag,
                   "der": abs_der})
fname = 'Traj_' + bas + '.csv'
df.to_csv(fname, index=False)<jupyter_output><empty_output><jupyter_text>Energy, $\theta_{opt} $, and stability for $\alpha=1$<jupyter_code>j_opt = np.argmin(abs_der)
if j_opt == len(thetas)-1:
    j_opt -= 1
theta_opt = thetas[j_opt]
Eres=es[j_opt]
print('theta_opt', thetas[j_opt])
print('E_res', Eres)
print('dE/dtheta', abs_der[j_opt])
print('delta-E down', es[j_opt-1]-es[j_opt])
print('delta-E up  ', es[j_opt+1]-es[j_opt])<jupyter_output>theta_opt 11.040000000000001
E_res (3.1692295664154315-0.14411821624242616j)
dE/dtheta 0.004881233243806652
delta-E down (-0.0005794423188767439-0.0005237844928482893j)
delta-E up   (0.0005591583569932723+0.0005466002429930505j)<jupyter_text>Don't optimize $\eta=\alpha\,\exp(i\theta)$**For GTOs the primitive numerical derivative is unstable.**<jupyter_code>def objective(eta, dx, E_near, verbose=False):
    """
    Input: 
        eta: (alpha, theta); z = alpha*exp(i*theta)
        dx:  numerical derivatives using z +/- dx or idx
        E_near: select the eigenvalue nearest to this energy  
    Output: 
        dE/dz along Re(z), dE/dz along Im(z) for the state near E_near
    This requires four single points.
    """
    a, t = eta
    z0 = a*np.exp(1j*t)
    es = np.zeros(2, complex)
    zs = np.array([z0-dx, z0+dx])
    for i, z in enumerate(zs):
        a, t = np.abs(z), np.angle(z)
        es[i] = near_E(a, t, E_near)
    dEdz_re = (es[1] - es[0])/(2*dx)
    zs = [z0-1j*dx, z0+1j*dx]
    for i, z in enumerate(zs):
        a, t = np.abs(z), np.angle(z)
        es[i] = near_E(a, t, E_near)
    dEdz_im = (es[1] - es[0])/(2j*dx)
    if verbose:
        print(dEdz_re)
        print(dEdz_im)
    return np.abs(0.5*(dEdz_re + dEdz_im)) 

def ECS(alpha, theta):
    """ diag H(alpha*exp(i*theta)) and return eigenvalues """
    H_theta = Bas.H_theta(theta, alpha)
    return eigvals(H_theta, b=S)

def near_E(alpha, theta, E_near):
    """ diag H(alpha*exp(i*theta)) and return eigenvalue near E_near """
    Es = ECS(alpha, theta)
    j = np.argmin(np.abs(Es - E_near))
    return Es[j]
alpha_curr=1.0
theta_curr=theta_opt/180.0*np.pi
dx=0.001
E_find = Eres/au2eV
objective((alpha_curr, theta_curr), dx, E_find, verbose=True)
args=(dx, E_find)
p0=[alpha_curr, theta_curr]
print('p0=', p0)
print('Eres0=', near_E(alpha_curr, theta_curr, E_find)*au2eV)
print('f0=', objective(p0, dx, E_find))
res=minimize(objective, p0, args=args, method='Nelder-Mead')
print(res.message)
print(res.x)
print(res.fun)
Eopt=near_E(res.x[0], res.x[1], E_find)*au2eV
print('Eres=', Eopt)
print('Energy change=', Eopt-Eres)<jupyter_output>p0= [1.0, 0.19268434942017398]
Eres0= (3.1692295664154315-0.14411821624242616j)
f0= 0.4828156933534529
Optimization terminated successfully.
[1.15969003 0.20451713]
0.0246259033514554
Eres= (3.64741613569381+0.8421105067721848j)
Energy change= (0.47818656927837866+0.986228723014611j)
{'fraction_non_alphanumeric': 0.1669133553455539, 'fraction_numerical': 0.16156210861892292, 'mean_word_length': 3.730210016155089, 'pattern_counts': {'":': 4, '<': 34, '<?xml version=': 0, '>': 35, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 3, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '10', 'dir': 'jupyter-structured-clean-dedup', 'id': '5618', 'license': 'MIT', 'n_tokens_mistral': 4881, 'n_tokens_neox': 3984, 'n_words': 789, 'path': 'notebooks/CS/.ipynb_checkpoints/CS_J3D-GTO-checkpoint.ipynb', 'repo_name': 'tsommerfeld/L2-methods_for_resonances', 'size': '8783'}