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