1*     **********************************
2*     *	                               *
3*     *          BLCJ_ion_Coulomb      *
4*     *                                *
5*     **********************************
6
7*   This routine calculate the Coulomb contribution
8*  to the BLCJ water psp.
9*
10
11      subroutine BLCJ_ion_Coulomb(RO,R1,R2,
12     >                        rion,zv,
13     >                        ec)
14      implicit none
15      real*8 RO(3),R1(3),R2(3)
16      real*8 rion(3),zv
17      real*8 ec
18
19*     **** Parameters ****
20#include "waterpsp_param.fh"
21
22*     **** local variables ****
23      real*8 R3(3)
24      real*8 x1,y1,z1,d1
25      real*8 x2,y2,z2,d2
26      real*8 x3,y3,z3,d3
27
28*     **** define R3 = center of negative charge ***
29      R3(1) = RO(1) + DELTA*(R1(1)+R2(1)-2.0d0*RO(1))
30      R3(2) = RO(2) + DELTA*(R1(2)+R2(2)-2.0d0*RO(2))
31      R3(3) = RO(3) + DELTA*(R1(3)+R2(3)-2.0d0*RO(3))
32
33
34      x1 = rion(1)-R1(1)
35      x2 = rion(1)-R2(1)
36      x3 = rion(1)-R3(1)
37
38      y1 = rion(2)-R1(2)
39      y2 = rion(2)-R2(2)
40      y3 = rion(2)-R3(2)
41
42      z1 = rion(3)-R1(3)
43      z2 = rion(3)-R2(3)
44      z3 = rion(3)-R3(3)
45
46      d1 = dsqrt(x1*x1 + y1*y1 + z1*z1)
47      d2 = dsqrt(x2*x2 + y2*y2 + z2*z2)
48      d3 = dsqrt(x3*x3 + y3*y3 + z3*z3)
49
50      ec = zv*( q1/d1
51     >        + q2/d2
52     >        + q3/d3)
53
54      return
55      end
56
57
58
59*     **********************************
60*     *	                               *
61*     *     BLCJ_ion_Coulomb_Fion      *
62*     *                                *
63*     **********************************
64
65*   This routine calculate the Coulomb contribution
66*  to the BLCJ water psp.
67*
68
69      subroutine BLCJ_ion_Coulomb_Fion(RO,R1,R2,
70     >                        rion,zv,fion)
71      implicit none
72      real*8 RO(3),R1(3),R2(3)
73      real*8 rion(3),zv
74      real*8 fion(3)
75
76*     **** Parameters ****
77#include "waterpsp_param.fh"
78
79
80*     **** local variables ****
81      real*8 R3(3)
82      real*8 x1,y1,z1,d1
83      real*8 x2,y2,z2,d2
84      real*8 x3,y3,z3,d3
85      real*8 dec1,dec2,dec3
86
87
88*     **** define R3 = center of negative charge ***
89      R3(1) = RO(1) + DELTA*(R1(1)+R2(1)-2*RO(1))
90      R3(2) = RO(2) + DELTA*(R1(2)+R2(2)-2*RO(2))
91      R3(3) = RO(3) + DELTA*(R1(3)+R2(3)-2*RO(3))
92
93
94      x1 = rion(1)-R1(1)
95      x2 = rion(1)-R2(1)
96      x3 = rion(1)-R3(1)
97
98      y1 = rion(2)-R1(2)
99      y2 = rion(2)-R2(2)
100      y3 = rion(2)-R3(2)
101
102      z1 = rion(3)-R1(3)
103      z2 = rion(3)-R2(3)
104      z3 = rion(3)-R3(3)
105
106      d1 = dsqrt(x1*x1 + y1*y1 + z1*z1)
107      d2 = dsqrt(x2*x2 + y2*y2 + z2*z2)
108      d3 = dsqrt(x3*x3 + y3*y3 + z3*z3)
109
110
111      dec1 = -zv*q1/(d1*d1)
112      dec2 = -zv*q2/(d2*d2)
113      dec3 = -zv*q3/(d3*d3)
114
115      fion(1) = -(x1/d1)*dec1 - (x2/d2)*dec2 - (x3/d3)*dec3
116      fion(2) = -(y1/d1)*dec1 - (y2/d2)*dec2 - (y3/d3)*dec3
117      fion(3) = -(z1/d1)*dec1 - (z2/d2)*dec2 - (z3/d3)*dec3
118
119      return
120      end
121
122
123*     **********************************
124*     *	                               *
125*     *    BLCJ_ion_Coulomb_Fwater     *
126*     *                                *
127*     **********************************
128
129*   This routine calculate the Coulomb contribution
130*  to the BLCJ water psp.
131*
132
133      subroutine BLCJ_ion_Coulomb_Fwater(RO,R1,R2,
134     >                        rion,zv,fo,f1,f2)
135      implicit none
136      real*8 RO(3),R1(3),R2(3)
137      real*8 rion(3),zv
138      real*8 fo(3),f1(3),f2(3)
139
140*     **** Parameters ****
141#include "waterpsp_param.fh"
142
143*     **** local variables ****
144      real*8 R3(3)
145      real*8 x1,y1,z1,d1
146      real*8 x2,y2,z2,d2
147      real*8 x3,y3,z3,d3
148      real*8 dec1,dec2,dec3
149      real*8 d3do
150
151
152*     **** define R3 = center of negative charge ***
153      R3(1) = RO(1) + DELTA*(R1(1)+R2(1)-2.0d0*RO(1))
154      R3(2) = RO(2) + DELTA*(R1(2)+R2(2)-2.0d0*RO(2))
155      R3(3) = RO(3) + DELTA*(R1(3)+R2(3)-2.0d0*RO(3))
156      d3do  = (1.0d0-2.0d0*DELTA)
157
158      x1 = rion(1)-R1(1)
159      x2 = rion(1)-R2(1)
160      x3 = rion(1)-R3(1)
161
162      y1 = rion(2)-R1(2)
163      y2 = rion(2)-R2(2)
164      y3 = rion(2)-R3(2)
165
166      z1 = rion(3)-R1(3)
167      z2 = rion(3)-R2(3)
168      z3 = rion(3)-R3(3)
169
170      d1 = dsqrt(x1*x1 + y1*y1 + z1*z1)
171      d2 = dsqrt(x2*x2 + y2*y2 + z2*z2)
172      d3 = dsqrt(x3*x3 + y3*y3 + z3*z3)
173
174      dec1 = -zv*q1/(d1*d1)
175      dec2 = -zv*q2/(d2*d2)
176      dec3 = -zv*q3/(d3*d3)
177
178      fo(1) =   d3do*(x3/d3)*dec3
179      fo(2) =   d3do*(y3/d3)*dec3
180      fo(3) =   d3do*(z3/d3)*dec3
181
182      f1(1) =   (x1/d1)*dec1 + DELTA*(x3/d3)*dec3
183      f1(2) =   (y1/d1)*dec1 + DELTA*(y3/d3)*dec3
184      f1(3) =   (z1/d1)*dec1 + DELTA*(z3/d3)*dec3
185
186      f2(1) =   (x2/d2)*dec2 + DELTA*(x3/d3)*dec3
187      f2(2) =   (y2/d2)*dec2 + DELTA*(y3/d3)*dec3
188      f2(3) =   (z2/d2)*dec2 + DELTA*(z3/d3)*dec3
189
190      return
191      end
192
193
194c $Id$
195