-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathsbasis.f90
151 lines (146 loc) · 4.09 KB
/
sbasis.f90
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
! File created at Fri Jun 5 21:58:57 PDT 2020
! Original source code: sbasis.f
subroutine sbasis (rho,cvi,eint,ilev,jlev,klev,nlev,mode)
implicit double precision (a-h,o-z)
!
! -----------------------------------------------------------------
! This subroutine constructs a multiple-arrangement
! hyperspherical basis set at hyperradius rho.
! -----------------------------------------------------------------
!
! common blocks
!
common /arrays/ mro,mvi,nvi,n
common /ranges/ rmin,rmax,smax
common /rotors/ jtot,ipar,jpar,jmax,kmin,kmax
!
! input arrays
!
dimension cvi(nvi,n)
dimension eint(n),ilev(n),jlev(n),klev(n),nlev(n)
!
! local arrays
!
! dimension h(nvi,nvi),t(nvi),e(nvi)
allocatable h(:,:),t(:),e(:)
! dimension c(0:2*nvi),s(0:2*nvi),v(0:2*nvi)
allocatable c(:),s(:),v(:)
! dimension wvi(mvi),xvi(mvi)
allocatable wvi(:),xvi(:)
allocate (h(nvi,nvi),t(nvi),e(nvi))
allocate (c(0:2*nvi),s(0:2*nvi),v(0:2*nvi))
allocate (wvi(mvi),xvi(mvi))
!
! vibrational quadrature rule
!
tmin = 0.d0
tmax = asin(min(1.d0,smax/rho))
call qvib (tmin,tmax,mvi,wvi,xvi)
!
! kinetic energy
!
pi = acos(-1.d0)
tscale = pi/tmax
do jvi = 1,nvi
t(jvi) = ((jvi*tscale)**2-0.25d0)/rho**2
enddo
do jvi = 0,2*nvi
c(jvi) = 0.d0
s(jvi) = 0.d0
enddo
do kvi = 1,mvi
theta = xvi(kvi)
weight = wvi(kvi)/tmax
oncos2 = weight/cos(theta)**2
onsin2 = weight/sin(theta)**2
arg = tscale*theta
do jvi = 0,2*nvi
cosine = cos(jvi*arg)
c(jvi) = c(jvi)+cosine*oncos2
s(jvi) = s(jvi)+cosine*onsin2
enddo
enddo
!
! loop over arrangements
!
na = 3-iabs(jpar)
do ia = 1,na
!
! reference potential
!
do jvi = 0,2*nvi
v(jvi) = 0.d0
enddo
do kvi = 1,mvi
theta = xvi(kvi)
weight = wvi(kvi)/tmax
sa = rho*sin(theta)
call potenl (100.d0,sa,0.d0,va,ia)
vtheta = weight*va
arg = tscale*theta
do jvi = 0,2*nvi
cosine = cos(jvi*arg)
v(jvi) = v(jvi)+cosine*vtheta
enddo
enddo
!
! loop over rotational quantum numbers
!
if (ia .eq. 1) then
jmin = (1-jpar)/2
jinc = 1+iabs(jpar)
else
jmin = 0
jinc = 1
endif
do j = jmin,jmax,jinc
kmaxj = min(j,kmax)
do k = kmin,kmaxj
!
! vibrational eigenvalue problem
!
if (mode .eq. 0) then
!
! k-dependent vibrational functions
! in the interaction region
!
ctot = (jtot*(jtot+1)+j*(j+1)-2*k*k)/rho**2
else
!
! k-independent vibrational functions
! in the asymptotic region
!
ctot = 0.d0
endif
cent = j*(j+1)/rho**2
do jvi = 1,nvi
do ivi = 1,jvi
h(ivi,jvi) = ctot*(c(jvi-ivi)-c(jvi+ivi)) &
+ cent*(s(jvi-ivi)-s(jvi+ivi)) &
+ (v(jvi-ivi)-v(jvi+ivi))
h(jvi,ivi) = h(ivi,jvi)
enddo
h(jvi,jvi) = h(jvi,jvi)+t(jvi)
enddo
call symevp (h,nvi,nvi,e,ierr)
if (ierr .ne. 0) stop 'sbasis 1'
!
! specified surface basis functions
!
do kvi = 1,nvi
do 1 i = 1,n
if (ilev(i) .ne. ia) go to 1
if (jlev(i) .ne. j) go to 1
if (klev(i) .ne. k) go to 1
if (nlev(i) .ne. kvi-1) go to 1
do jvi = 1,nvi
cvi(jvi,i) = h(jvi,kvi)
enddo
eint(i) = e(kvi)
1 continue
enddo
enddo
enddo
enddo
return
end