-
Notifications
You must be signed in to change notification settings - Fork 4
/
screen_pot_large_deltaE.py
164 lines (130 loc) · 5.14 KB
/
screen_pot_large_deltaE.py
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
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
#!/usr/bin/python
import numpy as np
import sciconv as sc
import wellenfkt as wf
import sys
import warnings
import math
# don't print warnings unless python -W ... is used
if not sys.warnoptions:
warnings.simplefilter("ignore")
#--------------------------------------------------------------------------
# input parameters
mu = wf.red_mass_au(20.1797,20.1797)
gs_de = 0.0001102
gs_a = 1.5
#gs_a = 3.8
gs_Req = 6.0
gs_const = 0.0
nmax_res = 1
nmax_fin = 1
n_vis_resstates = 2
n_vis_finstates = 4
De_min_eV = 0.01
De_max_eV = 6.0
#De_max_eV = 0.2
alpha_max = 30.0
De_step_eV = 0.1
alpha_step = 0.5
#alpha_in_step = 0.05
FCmin = 0.4
FCmin_fin = 0.2
FCmin_res = 0.2
FCmax = 1.0E-3
minEdiff_eV = 0.19
res_Req = 6.0
fin_Req = 6.0
R_min = sc.angstrom_to_bohr(1.5)
R_max = sc.angstrom_to_bohr(30.0)
#--------------------------------------------------------------------------
#convert units
De_min = sc.ev_to_hartree(De_min_eV)
De_max = sc.ev_to_hartree(De_max_eV)
minEdiff = sc.ev_to_hartree(minEdiff_eV)
De_step = sc.ev_to_hartree(De_step_eV)
#--------------------------------------------------------------------------
#definition of functions
def alpha_min(mu,De,nmax):
return 2* np.sqrt(2*mu*De) / (1 + 2*(nmax+1))
def lambda_param(mu,De,alpha):
return np.sqrt(2*mu*De) / alpha
def is_correct_nmax(l_param,nmax):
calc_n = int(l_param - 0.5)
if (calc_n == nmax):
return True
else:
return False
#--------------------------------------------------------------------------
#--------------------------------------------------------------------
# 2 lambda, viele (10) mu
#--------------------------------------------------------------------
#start with resonant state
De_res = De_min
while (De_res < De_max):
res_alpha = alpha_min(mu,De_res,nmax_res)
while (is_correct_nmax(lambda_param(mu,De_res,res_alpha+alpha_step),nmax_res)
and (res_alpha <= alpha_max)):
res_alpha = res_alpha + alpha_step
# consider only cases with minimum energy gap between vibrational states
res_evs = []
for i in range(0,nmax_res+1):
tmp = wf.eigenvalue(i,De_res,res_alpha,mu)
res_evs.append(tmp)
#if (nmax_res == 1):
# if ((res_evs[1]-res_evs[0]) < minEdiff):
# continue
# consider only cases, where the overlap of both res vib states with the
# ground state is sufficient
FCres = []
for i in range(0,nmax_res+1):
FC_tmp = wf.FC(i,res_alpha,res_Req,De_res,mu,
0,gs_a,gs_Req,gs_de,R_min,R_max)
FCres.append(FC_tmp)
if any(math.isnan(x) for x in FCres):
break
n_realistic_res_states = sum(abs(x) >= FCmin_res for x in FCres)
if (n_realistic_res_states < n_vis_resstates):
#if not all(abs(x) >= FCmin_fin for x in FCfins):
continue
De_fin = De_min
while (De_fin < De_max):
# consider only cases with minimum energy gap between vibrational states
fin_evs = []
fin_alpha = alpha_min(mu,De_fin,nmax_fin)
while (is_correct_nmax(lambda_param(mu,De_fin,fin_alpha+alpha_step),nmax_fin) and (fin_alpha <= alpha_max)):
fin_alpha = fin_alpha + alpha_step
for i in range(0,nmax_fin+1):
tmp = wf.eigenvalue(i,De_fin,fin_alpha,mu)
fin_evs.append(tmp)
if (nmax_fin == 1):
if ((fin_evs[1]-fin_evs[0]) < minEdiff):
continue
# consider only cases, where the overlap of both res vib states with the
# ground state is sufficient
FCfins = []
for i in range(0,nmax_res+1):
for j in range(0,nmax_fin+1):
FC_tmp = wf.FC(i,res_alpha,res_Req,De_res,mu,
j,fin_alpha,fin_Req,De_fin,R_min,R_max)
FCfins.append(FC_tmp)
if any(math.isnan(x) for x in FCfins):
break
#print FCfins
n_realistic_fin_states = sum(abs(x) >= FCmin_fin for x in FCfins)
if (n_realistic_fin_states < n_vis_finstates):
#if not all(abs(x) >= FCmin_fin for x in FCfins):
continue
print "n_realistic_res_states = ", n_realistic_res_states
print "n_realistic_fin_states = ", n_realistic_fin_states
print "YES!!!!!!!!!!!!!!!!!!!!!!!!"
print 'resonant', sc.hartree_to_ev(De_res), De_res, res_alpha
print 'final ', sc.hartree_to_ev(De_fin), De_fin, fin_alpha
print 'FCres', FCres
print 'FCfins', FCfins
if (nmax_res+1 >= 2):
print 'Ediff_res', sc.hartree_to_ev(res_evs[1]-res_evs[0])
if (nmax_fin+1 >= 2):
print 'Ediff_fin', sc.hartree_to_ev(fin_evs[1]-fin_evs[0])
print '---------------------------------------------------'
De_fin = De_fin + De_step
De_res = De_res + De_step