2121! PCMSolver API, see: <http://pcmsolver.github.io/pcmsolver-doc>
2222! pcmsolver_copyright_end
2323
24- module metal_sphere
24+ module metal_sphere
2525
26- use , intrinsic :: iso_c_binding
26+ use , intrinsic :: iso_c_binding
2727
28- implicit none
28+ implicit none
2929
30- public greens_function
31-
32- private
30+ public greens_function
31+
32+ private
33+ ! integer(kind=regint_k) types
34+ ! 32-bit integer(kind=regint_k)s
35+ integer , parameter :: regint_k = selected_int_kind (8 )
36+ ! 64-bit integer(kind=regint_k)s
37+ integer , parameter :: largeint_k = selected_int_kind (18 )
38+
39+ ! Real types
40+ ! Single-precision real
41+ integer , parameter :: sp = kind (1.0 )
42+ ! Double-precision real
43+ integer , parameter :: dp = selected_real_kind (2 * precision (1.0_sp ))
3344
34- contains
45+ contains
3546
3647 subroutine greens_function (epssol , epsre , epsim , radius , &
3748 ps , p1 , p2 , greenre , greenim ) bind(c, name= ' greens_function' )
3849
39- ! Passed variables
50+ ! Passed variables
4051 real (c_double), intent (in ) :: epssol, epsre, epsim, radius
4152 real (c_double), intent (in ) :: p1(* ), p2(* ), ps(* )
4253 real (c_double), intent (out ) :: greenre, greenim
43- ! Local variables
44- complex (kind= 8 ) :: green, eps2, ui
54+ ! Local variables
55+ complex (kind= dp ) :: green, eps2, ui
4556
46- ui = (0.0 , 1.0 )
57+ ui = (0.0_dp , 1.0_dp )
4758 eps2 = epsre + epsim * ui
4859 green = gsfera(epssol, eps2, ps(1 ), ps(2 ), ps(3 ), radius, &
4960 p1(1 ), p1(2 ), p1(3 ), &
@@ -53,25 +64,25 @@ subroutine greens_function(epssol, epsre, epsim, radius, &
5364
5465 end subroutine greens_function
5566
56- complex (kind= 8 ) function gsfera(eps, eps2, xs, ys, zs, rs, &
67+ complex (kind= dp ) function gsfera(eps, eps2, xs, ys, zs, rs, &
5768 xi, yi, zi, xj, yj, zj)
58- ! Passed variables
59- real (8 ), intent (in ) :: eps
60- complex (kind= 8 ), intent (in ) :: eps2
61- real (8 ), intent (in ) :: xs, ys, zs, rs ! Sphere center and radius
62- real (8 ), intent (in ) :: xi, yi, zi ! Source point
63- real (8 ), intent (in ) :: xj, yj, zj ! Probe point
64- ! Local variables
65- complex (kind= 8 ) :: coefl
66- real (8 ) :: di, dj
67- real (8 ) :: dim, xim, yim, zim, qim ! Image quantities
68- real (8 ) :: gc, cost, aa, arg, argl, gsfera2
69- integer :: maxl, l, m
69+ ! Passed variables
70+ real (kind = dp ), intent (in ) :: eps
71+ complex (kind= dp ), intent (in ) :: eps2
72+ real (kind = dp ), intent (in ) :: xs, ys, zs, rs ! Sphere center and radius
73+ real (kind = dp ), intent (in ) :: xi, yi, zi ! Source point
74+ real (kind = dp ), intent (in ) :: xj, yj, zj ! Probe point
75+ ! Local variables
76+ complex (kind= dp) :: coefl
77+ real (kind = dp ) :: di, dj
78+ real (kind = dp ) :: dim, xim, yim, zim, qim ! Image quantities
79+ real (kind = dp ) :: gc, cost, aa, arg, argl, gsfera2
80+ integer (kind = regint_k) :: maxl, l, m
7081
71- ! lf print *, eps,eps2
72- ! lf print *, 'Sphere', xs,ys,zs,rs
73- ! lf print *, 'Source', xi,yi,zi
74- ! lf print *, 'Probe ', xj,yj,zj
82+ ! lf print *, eps,eps2
83+ ! lf print *, 'Sphere', xs,ys,zs,rs
84+ ! lf print *, 'Source', xi,yi,zi
85+ ! lf print *, 'Probe ', xj,yj,zj
7586 di = sqrt ((xi - xs)** 2 + (yi - ys)** 2 + (zi - zs)** 2 )
7687 dj = sqrt ((xj - xs)** 2 + (yj - ys)** 2 + (zj - zs)** 2 )
7788 dim = rs** 2 / di
@@ -86,41 +97,41 @@ complex(kind=8) function gsfera(eps, eps2, xs, ys, zs, rs, &
8697 aa = (rs * rs / (di * dj))** 4 .
8798 maxl = int (800 . / (abs (eps2 + eps))** 0.4 * aa)
8899 maxl = max (maxl, 10 )
89- ! lf write (6,*) "maxl",maxl
100+ ! lf write (6,*) "maxl",maxl
90101 arg = rs * rs / (di * dj)
91102 argl = rs / (di * dj)
92- gsfera2 = 0.0d0
103+ gsfera2 = 0.0_dp
93104 m = 0
94105 do l = 1 , maxl
95106 argl = argl * arg
96- coefl = (eps2 - eps) * float(l) / ((eps2 + eps) * float(l) + 1.0d0 )
107+ coefl = (eps2 - eps) * float(l) / ((eps2 + eps) * float(l) + 1.0_dp )
97108 coefl = coefl - (eps2 - eps) / (eps2 + eps)
98109 gsfera = gsfera + coefl * argl * legendre_polynomial(l, m, cost)
99110 enddo
100111 gsfera = gsfera + gc * (eps2- eps) / (eps2 + eps)
101- ! lf: WARNING: sign change and divided by epsilon of solvent!!!!!!
112+ ! lf: WARNING: sign change and divided by epsilon of solvent!!!!!!
102113 gsfera = - gsfera / eps
103- ! DEBUGGING STUFF
104- ! arg = rs * rs / (di * dj)
105- ! argl = rs / (di * dj)
106- ! gsfera = 0.0d0
107- ! m = 0
108- ! do l = 1, 5000
109- ! argl = argl * arg
110- ! coefl = (eps2 - eps) * float(l) / ((eps2 + eps) * float(l) + 1.0d0)
111- ! gsfera = gsfera + coefl * argl * legendre_polynomial(l, m, cost)
112- ! enddo
113- ! lf write (6, *) "gsfera2", gsfera
114+ ! DEBUGGING STUFF
115+ ! arg = rs * rs / (di * dj)
116+ ! argl = rs / (di * dj)
117+ ! gsfera = 0.0d0
118+ ! m = 0
119+ ! do l = 1, 5000
120+ ! argl = argl * arg
121+ ! coefl = (eps2 - eps) * float(l) / ((eps2 + eps) * float(l) + 1.0d0)
122+ ! gsfera = gsfera + coefl * argl * legendre_polynomial(l, m, cost)
123+ ! enddo
124+ ! lf write (6, *) "gsfera2", gsfera
114125 end function gsfera
115126
116- real (8 ) function legendre_polynomial(l, m, x)
117- ! Computes the associated Legendre polynomial P_l^m(x)
118- ! Passed variables
119- integer , intent (in ) :: l, m
120- real (8 ), intent (in ) :: x
121- ! Local variables
122- integer :: i, ll
123- real (8 ) :: fact, pll, pmm, pmmp1, somx2
127+ real (kind = dp ) function legendre_polynomial(l, m, x)
128+ ! Computes the associated Legendre polynomial P_l^m(x)
129+ ! Passed variables
130+ integer (kind = regint_k) , intent (in ) :: l, m
131+ real (kind = dp ), intent (in ) :: x
132+ ! Local variables
133+ integer (kind = regint_k) :: i, ll
134+ real (kind = dp ) :: fact, pll, pmm, pmmp1, somx2
124135
125136 pmm = 1.0d0
126137 if (m .gt. 0 ) then
@@ -149,4 +160,4 @@ real(8) function legendre_polynomial(l, m, x)
149160
150161 end function legendre_polynomial
151162
152- end module metal_sphere
163+ end module metal_sphere
0 commit comments