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
|
subroutine laplac(r,gradup,graddn,i,imin,imax,done,
j lapl,laplup,lapldn)
c calculates laplace(density), i.e. (del^2)(density)
c spherical coordinates => (del^2)phi = 1/r^2*d/dr[r^2*d(phi)/dr]
c input r : radial coordinate (array)
c input gradup : grad(density), spin up (array)
c input graddn : grad(density), spin down (array)
c input i : counter, i.e. current position = r(i)
c input imin : min value of i
c input imax : max value of i
c in/output done : true if r^2*grad(density) is splined
c output lapl : laplace(density)
c output laplup,lapldn : laplace(density), spin up and down
implicit logical (a-z)
logical done
double precision r,gradup,graddn,lapl,laplup,lapldn,lapl2,temp
double precision rj2,r2,deriv
integer i,imin,imax,j
common/ggalap/temp,lapl2
dimension r(561),gradup(561),graddn(561)
dimension lapl2(2,561),temp(2,561)
if (.not.done) then
do 10 j = imin,imax
rj2 = r(j)**2
temp(1,j) = rj2*gradup(j)
temp(2,j) = rj2*graddn(j)
10 continue
call spline(r,imin+1,imax,temp,lapl2)
done = .true.
endif
if (i.eq.imin) then
call laplac2(r,gradup,graddn,i+1,imin,imax,done,lapl,laplup,
j lapldn)
else
r2 = r(i)**2
laplup = (deriv(r,temp,lapl2,1,i,imin,imax))/r2
lapldn = (deriv(r,temp,lapl2,-1,i,imin,imax))/r2
lapl = laplup + lapldn
endif
return
end
|
