| 1 |
subroutine doMultipolePair(atom1, atom2, d, rij, r2, rcut, sw, & |
| 2 |
vpair, fpair, pot, eFrame, f, t, do_pot) |
| 3 |
|
| 4 |
!*************************************************************** |
| 5 |
! doMultipolePair evaluates the potential, forces, and torques |
| 6 |
! between point multipoles. It is based on the ewald2 real |
| 7 |
! space routine in the MDMULP code written by W. Smith and |
| 8 |
! available from the CCP5 website. |
| 9 |
! |
| 10 |
! We're using the damped real space portion of the Ewald sum |
| 11 |
! as the entire interaction. Details on this portion of the |
| 12 |
! sum can be found in: |
| 13 |
! |
| 14 |
! "Point Multipoles in the Ewald Summation (Revisited)," by |
| 15 |
! W. Smith, CCP5 Newsletter, 46, pp. 18-30 (1998). |
| 16 |
! |
| 17 |
!************************************************************** |
| 18 |
|
| 19 |
logical, intent(in) :: do_pot |
| 20 |
|
| 21 |
integer, intent(in) :: atom1, atom2 |
| 22 |
real(kind=dp), intent(in) :: rij, r2, sw, rcut |
| 23 |
real(kind=dp), intent(in), dimension(3) :: d |
| 24 |
real(kind=dp), intent(inout) :: vpair |
| 25 |
real(kind=dp), intent(inout), dimension(3) :: fpair |
| 26 |
|
| 27 |
real( kind = dp ) :: pot |
| 28 |
real( kind = dp ), dimension(9,nLocal) :: eFrame |
| 29 |
real( kind = dp ), dimension(3,nLocal) :: f |
| 30 |
real( kind = dp ), dimension(3,nLocal) :: t |
| 31 |
logical :: i_is_Charge, i_is_Dipole, i_is_SplitDipole, i_is_Quadrupole |
| 32 |
logical :: j_is_Charge, j_is_Dipole, j_is_SplitDipole, j_is_Quadrupole |
| 33 |
integer :: me1, me2, id1, id2 |
| 34 |
real (kind=dp) :: q_i, q_j, mu_i, mu_j, d_i, d_j |
| 35 |
real (kind=dp) :: qxx_i, qyy_i, qzz_i |
| 36 |
real (kind=dp) :: qxx_j, qyy_j, qzz_j |
| 37 |
|
| 38 |
real (kind=dp) :: epot, qii, alsq2, alsq2n, exp2a, ralpha |
| 39 |
real (kind=dp), dimension(3) :: ftm, ttm_i, ttm_j |
| 40 |
real (kind=dp), dimension(3) :: qir, qjr, qidj, qjdi, qiqjr, qjqir |
| 41 |
real (kind=dp), dimension(3) :: dixdj, qixqj, dixr, djxr, rxqir, rxqjr |
| 42 |
real (kind=dp), dimension(3) :: dixqjr, djxqir, rxqidj, rxqjdi, rxqijr |
| 43 |
real (kind=dp), dimension(3) :: rxqjir, qjrxqir |
| 44 |
real (kind=dp), dimension(6) :: bn |
| 45 |
real (kind=dp), dimension(10) :: sc |
| 46 |
real (kind=dp), dimension(9) :: gl |
| 47 |
|
| 48 |
real(kind=dp) :: a1, a2, a3, a4, a5, p |
| 49 |
|
| 50 |
DATA A1,A2,A3/0.254829592,-0.284496736,1.421413741/ |
| 51 |
DATA A4,A5,P/-1.453152027,1.061405429,0.3275911/ |
| 52 |
|
| 53 |
if (.not.summationMethodChecked) then |
| 54 |
call checkSummationMethod() |
| 55 |
endif |
| 56 |
|
| 57 |
#ifdef IS_MPI |
| 58 |
me1 = atid_Row(atom1) |
| 59 |
me2 = atid_Col(atom2) |
| 60 |
#else |
| 61 |
me1 = atid(atom1) |
| 62 |
me2 = atid(atom2) |
| 63 |
#endif |
| 64 |
|
| 65 |
! logicals |
| 66 |
i_is_Charge = ElectrostaticMap(me1)%is_Charge |
| 67 |
i_is_Dipole = ElectrostaticMap(me1)%is_Dipole |
| 68 |
i_is_Quadrupole = ElectrostaticMap(me1)%is_Quadrupole |
| 69 |
|
| 70 |
j_is_Charge = ElectrostaticMap(me2)%is_Charge |
| 71 |
j_is_Dipole = ElectrostaticMap(me2)%is_Dipole |
| 72 |
j_is_Quadrupole = ElectrostaticMap(me2)%is_Quadrupole |
| 73 |
|
| 74 |
if (i_is_Charge) then |
| 75 |
q_i = ElectrostaticMap(me1)%charge |
| 76 |
else |
| 77 |
q_i = 0.0_dp |
| 78 |
endif |
| 79 |
|
| 80 |
if (j_is_Charge) then |
| 81 |
q_j = ElectrostaticMap(me2)%charge |
| 82 |
else |
| 83 |
q_j = 0.0_dp |
| 84 |
endif |
| 85 |
|
| 86 |
if (i_is_Dipole) then |
| 87 |
mu_i = ElectrostaticMap(me1)%dipole_moment |
| 88 |
#ifdef IS_MPI |
| 89 |
d_i(1) = eFrame_Row(3,atom1) |
| 90 |
d_i(2) = eFrame_Row(6,atom1) |
| 91 |
d_i(3) = eFrame_Row(9,atom1) |
| 92 |
#else |
| 93 |
d_i(1) = eFrame(3,atom1) |
| 94 |
d_i(2) = eFrame(6,atom1) |
| 95 |
d_i(3) = eFrame(9,atom1) |
| 96 |
#endif |
| 97 |
d_i = d_i * mu_i |
| 98 |
else |
| 99 |
d_i = 0.0_dp |
| 100 |
endif |
| 101 |
|
| 102 |
if (j_is_Dipole) then |
| 103 |
mu_j = ElectrostaticMap(me2)%dipole_moment |
| 104 |
#ifdef IS_MPI |
| 105 |
d_j(1) = eFrame_Col(3,atom2) |
| 106 |
d_j(2) = eFrame_Col(6,atom2) |
| 107 |
d_j(3) = eFrame_Col(9,atom2) |
| 108 |
#else |
| 109 |
d_j(1) = eFrame(3,atom2) |
| 110 |
d_j(2) = eFrame(6,atom2) |
| 111 |
d_j(3) = eFrame(9,atom2) |
| 112 |
#endif |
| 113 |
d_j = d_j * mu_j |
| 114 |
else |
| 115 |
d_j = 0.0_dp |
| 116 |
endif |
| 117 |
|
| 118 |
if (i_is_Quadrupole) then |
| 119 |
qxx_i = ElectrostaticMap(me1)%quadrupole_moments(1) |
| 120 |
qyy_i = ElectrostaticMap(me1)%quadrupole_moments(2) |
| 121 |
qzz_i = ElectrostaticMap(me1)%quadrupole_moments(3) |
| 122 |
#ifdef IS_MPI |
| 123 |
qpole_i(1) = qxx_i * eFrame_Row(1,atom1) |
| 124 |
qpole_i(2) = qxx_i * eFrame_Row(4,atom1) |
| 125 |
qpole_i(3) = qxx_i * eFrame_Row(7,atom1) |
| 126 |
qpole_i(4) = qyy_i * eFrame_Row(2,atom1) |
| 127 |
qpole_i(5) = qyy_i * eFrame_Row(5,atom1) |
| 128 |
qpole_i(6) = qyy_i * eFrame_Row(8,atom1) |
| 129 |
qpole_i(7) = qzz_i * eFrame_Row(3,atom1) |
| 130 |
qpole_i(8) = qzz_i * eFrame_Row(6,atom1) |
| 131 |
qpole_i(9) = qzz_i * eFrame_Row(9,atom1) |
| 132 |
#else |
| 133 |
qpole_i(1) = qxx_i * eFrame(1,atom1) |
| 134 |
qpole_i(2) = qxx_i * eFrame(4,atom1) |
| 135 |
qpole_i(3) = qxx_i * eFrame(7,atom1) |
| 136 |
qpole_i(4) = qyy_i * eFrame(2,atom1) |
| 137 |
qpole_i(5) = qyy_i * eFrame(5,atom1) |
| 138 |
qpole_i(6) = qyy_i * eFrame(8,atom1) |
| 139 |
qpole_i(7) = qzz_i * eFrame(3,atom1) |
| 140 |
qpole_i(8) = qzz_i * eFrame(6,atom1) |
| 141 |
qpole_i(9) = qzz_i * eFrame(9,atom1) |
| 142 |
#endif |
| 143 |
else |
| 144 |
qpole_i = 0.0_dp |
| 145 |
endif |
| 146 |
|
| 147 |
if (j_is_Quadrupole) then |
| 148 |
qxx_j = ElectrostaticMap(me2)%quadrupole_moments(1) |
| 149 |
qyy_j = ElectrostaticMap(me2)%quadrupole_moments(2) |
| 150 |
qzz_j = ElectrostaticMap(me2)%quadrupole_moments(3) |
| 151 |
#ifdef IS_MPI |
| 152 |
qpole_j(1) = qxx_j * eFrame_Col(1,atom2) |
| 153 |
qpole_j(2) = qxx_j * eFrame_Col(4,atom2) |
| 154 |
qpole_j(3) = qxx_j * eFrame_Col(7,atom2) |
| 155 |
qpole_j(4) = qyy_j * eFrame_Col(2,atom2) |
| 156 |
qpole_j(5) = qyy_j * eFrame_Col(5,atom2) |
| 157 |
qpole_j(6) = qyy_j * eFrame_Col(8,atom2) |
| 158 |
qpole_j(7) = qzz_j * eFrame_Col(3,atom2) |
| 159 |
qpole_j(8) = qzz_j * eFrame_Col(6,atom2) |
| 160 |
qpole_j(9) = qzz_j * eFrame_Col(9,atom2) |
| 161 |
#else |
| 162 |
qpole_j(1) = qxx_j * eFrame(1,atom2) |
| 163 |
qpole_j(2) = qxx_j * eFrame(4,atom2) |
| 164 |
qpole_j(3) = qxx_j * eFrame(7,atom2) |
| 165 |
qpole_j(4) = qyy_j * eFrame(2,atom2) |
| 166 |
qpole_j(5) = qyy_j * eFrame(5,atom2) |
| 167 |
qpole_j(6) = qyy_j * eFrame(8,atom2) |
| 168 |
qpole_j(7) = qzz_j * eFrame(3,atom2) |
| 169 |
qpole_j(8) = qzz_j * eFrame(6,atom2) |
| 170 |
qpole_j(9) = qzz_j * eFrame(9,atom2) |
| 171 |
#endif |
| 172 |
else |
| 173 |
qpole_j = 0.0_dp |
| 174 |
endif |
| 175 |
|
| 176 |
! |
| 177 |
! INITIALISE TEMPORARY FORCE AND TORQUE ACCUMULATORS |
| 178 |
! |
| 179 |
|
| 180 |
ftm = 0.0 |
| 181 |
ttm_i = 0.0 |
| 182 |
ttm_j = 0.0 |
| 183 |
|
| 184 |
ri = 1.0_dp / rij |
| 185 |
ri2 = ri * ri |
| 186 |
|
| 187 |
! |
| 188 |
! CALCULATE BN FUNCTIONS |
| 189 |
! |
| 190 |
if (screeningMethod .eq. DAMPED) then |
| 191 |
! assemble the damping variables |
| 192 |
call lookupUniformSpline1d(erfcSpline, rij, erfcVal, derfcVal) |
| 193 |
else |
| 194 |
dampingAlpha = 0.0_dp |
| 195 |
erfcVal = 1.0_dp |
| 196 |
endif |
| 197 |
|
| 198 |
rrtpi = 1.0 / sqrt(pi) |
| 199 |
alsq2 = 2.0 * dampingAlpha**2 |
| 200 |
ralpi = 0.0 |
| 201 |
if(dampingAlpha .gt. 0.0) ralpi = rrtpi / dampingAlpha |
| 202 |
|
| 203 |
alphar = dampingAlpha * rij |
| 204 |
t = 1.0 / (1.0 + p*alphar) |
| 205 |
exp2a = exp(-alphar**2) |
| 206 |
bn(1) = ((((a5*t+a4)*t+a3)*t+a2)*t+a1) * t * exp2a * ri |
| 207 |
alsq2n = ralpi |
| 208 |
|
| 209 |
do n = 1, 5 |
| 210 |
bfac = real(n+n-1, kind=dp) |
| 211 |
alsq2n = alsq2*alsq2n |
| 212 |
bn(n+1) = ri2 * (bfac*bn(n) + alsq2n*exp2a) |
| 213 |
enddo |
| 214 |
|
| 215 |
ralpha = dampingAlpha * rij |
| 216 |
bn(0) = erfcVal * ri |
| 217 |
alsq2 = 2.0d0 * dampingAlpha**2 |
| 218 |
alsq2n = 0.0d0 |
| 219 |
if (dampingAlpha .gt. 0.0_dp) alsq2n = 1.0_dp / (sqrtpi*dampingAlpha) |
| 220 |
exp2a = exp(-ralpha**2) |
| 221 |
do m = 1, 5 |
| 222 |
bfac = real(m+m-1, kind=dp) |
| 223 |
alsq2n = alsq2 * alsq2n |
| 224 |
bn(m) = (bfac*bn(m-1)+alsq2n*exp2a) / r2 |
| 225 |
end do |
| 226 |
|
| 227 |
! |
| 228 |
! CONSTRUCT AUXILIARY VECTORS |
| 229 |
! |
| 230 |
|
| 231 |
dixdj(1) = d_i(2)*d_j(3) - d_i(3)*d_j(2) |
| 232 |
dixdj(2) = d_i(3)*d_j(1) - d_i(1)*d_j(3) |
| 233 |
dixdj(3) = d_i(1)*d_j(2) - d_i(2)*d_j(1) |
| 234 |
|
| 235 |
dixr(1) = d_i(2)*d(3) - d_i(3)*d(2) |
| 236 |
dixr(2) = d_i(3)*d(1) - d_i(1)*d(3) |
| 237 |
dixr(3) = d_i(1)*d(2) - d_i(2)*d(1) |
| 238 |
|
| 239 |
djxr(1) = d_j(2)*d(3) - d_j(3)*d(2) |
| 240 |
djxr(2) = d_j(3)*d(1) - d_j(1)*d(3) |
| 241 |
djxr(3) = d_j(1)*d(2) - d_j(2)*d(1) |
| 242 |
|
| 243 |
qir(1) = qpole_i(1)*d(1) + qpole_i(4)*d(2) + qpole_i(7)*d(3) |
| 244 |
qir(2) = qpole_i(2)*d(1) + qpole_i(5)*d(2) + qpole_i(8)*d(3) |
| 245 |
qir(3) = qpole_i(3)*d(1) + qpole_i(6)*d(2) + qpole_i(9)*d(3) |
| 246 |
|
| 247 |
qjr(1) = qpole_j(1)*d(1) + qpole_j(4)*d(2) + qpole_j(7)*d(3) |
| 248 |
qjr(2) = qpole_j(2)*d(1) + qpole_j(5)*d(2) + qpole_j(8)*d(3) |
| 249 |
qjr(3) = qpole_j(3)*d(1) + qpole_j(6)*d(2) + qpole_j(9)*d(3) |
| 250 |
|
| 251 |
qiqjr(1) = qpole_i(1)*qjr(1) + qpole_i(4)*qjr(2) + qpole_i(7)*qjr(3) |
| 252 |
qiqjr(2) = qpole_i(2)*qjr(1) + qpole_i(5)*qjr(2) + qpole_i(8)*qjr(3) |
| 253 |
qiqjr(3) = qpole_i(3)*qjr(1) + qpole_i(6)*qjr(2) + qpole_i(9)*qjr(3) |
| 254 |
|
| 255 |
|
| 256 |
qjqir(1) = qpole_j(1)*QIR(1) + qpole_j(4)*QIR(2) + qpole_j(7)*QIR(3) |
| 257 |
qjqir(2) = qpole_j(2)*QIR(1) + qpole_j(5)*QIR(2) + qpole_j(8)*QIR(3) |
| 258 |
qjqir(3) = qpole_j(3)*QIR(1) + qpole_j(6)*QIR(2) + qpole_j(9)*QIR(3) |
| 259 |
|
| 260 |
qixqj(1) = qpole_i(2)*qpole_j(3) + qpole_i(5)*qpole_j(6) + & |
| 261 |
qpole_i(8)*qpole_j(9) - qpole_i(3)*qpole_j(2) - & |
| 262 |
qpole_i(6)*qpole_j(5) - qpole_i(9)*qpole_j(8) |
| 263 |
|
| 264 |
qixqj(2) = qpole_i(3)*qpole_j(1) + qpole_i(6)*qpole_j(4) + & |
| 265 |
qpole_i(9)*qpole_j(7) - qpole_i(1)*qpole_j(3) - & |
| 266 |
qpole_i(4)*qpole_j(6) - qpole_i(7)*qpole_j(9) |
| 267 |
|
| 268 |
qixqj(3) = qpole_i(1)*qpole_j(2) + qpole_i(4)*qpole_j(5) + & |
| 269 |
qpole_i(7)*qpole_j(8) - qpole_i(2)*qpole_j(1) - & |
| 270 |
qpole_i(5)*qpole_j(4) - qpole_i(8)*qpole_j(7) |
| 271 |
|
| 272 |
rxqir(1) = d(2)*qir(3) - d(3)*qir(2) |
| 273 |
rxqir(2) = d(3)*qir(1) - d(1)*qir(3) |
| 274 |
rxqir(3) = d(1)*qir(2) - d(2)*qir(1) |
| 275 |
|
| 276 |
rxqjr(1) = d(2)*qjr(3) - d(3)*qjr(2) |
| 277 |
rxqjr(2) = d(3)*qjr(1) - d(1)*qjr(3) |
| 278 |
rxqjr(3) = d(1)*qjr(2) - d(2)*qjr(1) |
| 279 |
|
| 280 |
rxqijr(1) = d(2)*qiqjr(3) - d(3)*qiqjr(2) |
| 281 |
rxqijr(2) = d(3)*qiqjr(1) - d(1)*qiqjr(3) |
| 282 |
rxqijr(3) = d(1)*qiqjr(2) - d(2)*qiqjr(1) |
| 283 |
|
| 284 |
rxqjir(1) = d(2)*qjqir(3) - d(3)*qjqir(2) |
| 285 |
rxqjir(2) = d(3)*qjqir(1) - d(1)*qjqir(3) |
| 286 |
rxqjir(3) = d(1)*qjqir(2) - d(2)*qjqir(1) |
| 287 |
|
| 288 |
qjrxqir(1) = qjr(2)*qir(3) - qjr(3)*qir(2) |
| 289 |
qjrxqir(2) = qjr(3)*qir(1) - qjr(1)*qir(3) |
| 290 |
qjrxqir(3) = qjr(1)*qir(2) - qjr(2)*qir(1) |
| 291 |
|
| 292 |
qidj(1) = qpole_i(1)*d_j(1) + qpole_i(4)*d_j(2) + qpole_i(7)*d_j(3) |
| 293 |
qidj(2) = qpole_i(2)*d_j(1) + qpole_i(5)*d_j(2) + qpole_i(8)*d_j(3) |
| 294 |
qidj(3) = qpole_i(3)*d_j(1) + qpole_i(6)*d_j(2) + qpole_i(9)*d_j(3) |
| 295 |
|
| 296 |
qjdi(1) = qpole_j(1)*d_i(1) + qpole_j(4)*d_i(2) + qpole_j(7)*d_i(3) |
| 297 |
qjdi(2) = qpole_j(2)*d_i(1) + qpole_j(5)*d_i(2) + qpole_j(8)*d_i(3) |
| 298 |
qjdi(3) = qpole_j(3)*d_i(1) + qpole_j(6)*d_i(2) + qpole_j(9)*d_i(3) |
| 299 |
|
| 300 |
dixqjr(1) = d_i(2)*qjr(3) - d_i(3)*qjr(2) |
| 301 |
dixqjr(2) = d_i(3)*qjr(1) - d_i(1)*qjr(3) |
| 302 |
dixqjr(3) = d_i(1)*qjr(2) - d_i(2)*qjr(1) |
| 303 |
|
| 304 |
djxqir(1) = d_j(2)*qir(3) - d_j(3)*qir(2) |
| 305 |
djxqir(2) = d_j(3)*qir(1) - d_j(1)*qir(3) |
| 306 |
djxqir(3) = d_j(1)*qir(2) - d_j(2)*qir(1) |
| 307 |
|
| 308 |
rxqidj(1) = d(2)*qidj(3) - d(3)*qidj(2) |
| 309 |
rxqidj(2) = d(3)*qidj(1) - d(1)*qidj(3) |
| 310 |
rxqidj(3) = d(1)*qidj(2) - d(2)*qidj(1) |
| 311 |
|
| 312 |
rxqjdi(1) = d(2)*qjdi(3) - d(3)*qjdi(2) |
| 313 |
rxqjdi(2) = d(3)*qjdi(1) - d(1)*qjdi(3) |
| 314 |
rxqjdi(3) = d(1)*qjdi(2) - d(2)*qjdi(1) |
| 315 |
|
| 316 |
! |
| 317 |
! CALCULATE SCALAR PRODUCTS |
| 318 |
! |
| 319 |
|
| 320 |
qii = qpole_i(1) + qpole_i(5) + qpole_i(9) |
| 321 |
|
| 322 |
sc(1) = qpole_j(1) + qpole_j(5) + qpole_j(9) |
| 323 |
sc(2) = d_i(1)*d_j(1) + d_i(2)*d_j(2) + d_i(3)*d_j(3) |
| 324 |
sc(3) = d_i(1)*d(1) + d_i(2)*d(2) + d_i(3)*d(3) |
| 325 |
sc(4) = d_j(1)*d(1) + d_j(2)*d(2) + d_j(3)*d(3) |
| 326 |
sc(5) = qir(1)*d(1) + qir(2)*d(2) + qir(3)*d(3) |
| 327 |
sc(6) = qjr(1)*d(1) + qjr(2)*d(2) + qjr(3)*d(3) |
| 328 |
sc(7) = qir(1)*d_j(1) + qir(2)*d_j(2) + qir(3)*d_j(3) |
| 329 |
sc(8) = qjr(1)*d_i(1) + qjr(2)*d_i(2) + qjr(3)*d_i(3) |
| 330 |
sc(9) = qir(1)*qjr(1) + qir(2)*qjr(2) + qir(3)*qjr(3) |
| 331 |
sc(10) = qpole_i(1)*qpole_j(1) + qpole_i(2)*qpole_j(2) + & |
| 332 |
qpole_i(3)*qpole_j(3) + qpole_i(4)*qpole_j(4) + & |
| 333 |
qpole_i(5)*qpole_j(5) + qpole_i(6)*qpole_j(6) + & |
| 334 |
qpole_i(7)*qpole_j(7) + qpole_i(8)*qpole_j(8) + & |
| 335 |
qpole_i(9)*qpole_j(9) |
| 336 |
|
| 337 |
! |
| 338 |
! CALCULATE GL FUNCTIONS FOR POTENTIAL |
| 339 |
! |
| 340 |
|
| 341 |
gl(1) = q_i*q_j |
| 342 |
gl(2) = q_j*sc(3) - q_i*sc(4) |
| 343 |
gl(3) = q_i*sc(6) + q_j*sc(5) - sc(3)*sc(4) |
| 344 |
gl(4) = sc(3)*sc(6) - sc(4)*sc(5) |
| 345 |
gl(5) = sc(5)*sc(6) |
| 346 |
gl(7) = sc(2) - q_j*qii - q_i*sc(1) |
| 347 |
gl(8) = sc(4)*qii - sc(1)*sc(3) + 2.0*(sc(7) - sc(8)) |
| 348 |
gl(9) = 2.0*sc(10) + sc(1)*qii |
| 349 |
gl(6) = -(sc(1)*sc(5) + sc(6)*qii + 4.0*sc(9)) |
| 350 |
|
| 351 |
! |
| 352 |
! CALCULATE POTENTIAL AND VIRIAL |
| 353 |
! |
| 354 |
|
| 355 |
epot = bn(1)*gl(1) + bn(2)*(gl(2) + gl(7)) + bn(3)*(gl(3) & |
| 356 |
+ gl(8) + gl(9)) + bn(4)*(gl(4) + gl(6)) + bn(5)*gl(5) |
| 357 |
|
| 358 |
vpair = vpair + epot |
| 359 |
|
| 360 |
if (do_pot) then |
| 361 |
#ifdef IS_MPI |
| 362 |
pot_row(ELECTROSTATIC_POT,atom1) = pot_row(ELECTROSTATIC_POT,atom1) + & |
| 363 |
0.5_dp*epot*sw |
| 364 |
pot_col(ELECTROSTATIC_POT,atom2) = pot_col(ELECTROSTATIC_POT,atom2) + & |
| 365 |
0.5_dp*epot*sw |
| 366 |
#else |
| 367 |
pot = pot + epot*sw |
| 368 |
#endif |
| 369 |
endif |
| 370 |
|
| 371 |
! |
| 372 |
! CALCULATE FORCE AND TORQUE COEFFICIENTS |
| 373 |
! |
| 374 |
|
| 375 |
gl(1) = bn(2)*gl(1) + bn(3)*(gl(2) + gl(7)) + bn(4)*(gl(3) & |
| 376 |
+ gl(8) + gl(9)) + bn(5)*(gl(4) + gl(6)) + bn(6)*gl(5) |
| 377 |
gl(2) = -q_j*bn(2) + (sc(4) + sc(1))*bn(3) - sc(6)*bn(4) |
| 378 |
gl(3) = q_i*bn(2) + (sc(3) - qii)*bn(3) + sc(5)*bn(4) |
| 379 |
gl(4) = 2.0*bn(3) |
| 380 |
gl(5) = 2.0*( - q_j*bn(3) + (sc(1) + sc(4))*bn(4) - sc(6)* bn(5)) |
| 381 |
gl(6) = 2.0*( - q_i*bn(3) + (qii - sc(3))*bn(4) - sc(5)*bn(5)) |
| 382 |
gl(7) = 4.0*bn(4) |
| 383 |
|
| 384 |
! |
| 385 |
! CALCULATE FORCES BETWEEN MULTIPOLES I AND J |
| 386 |
! |
| 387 |
|
| 388 |
ftm(1) = gl(1)*d(1) + gl(2)*d_i(1) + gl(3)*d_j(1) + & |
| 389 |
gl(4)*(qjdi(1) - qidj(1)) + gl(5)*qir(1) + & |
| 390 |
gl(6)*qjr(1) + gl(7)*(qiqjr(1) + qjqir(1)) |
| 391 |
ftm(2) = gl(1)*d(2) + gl(2)*d_i(2) + gl(3)*d_j(2) + & |
| 392 |
gl(4)*(qjdi(2) - qidj(2)) + gl(5)*qir(2) + & |
| 393 |
gl(6)*qjr(2) + gl(7)*(qiqjr(2) + qjqir(2)) |
| 394 |
ftm(3) = gl(1)*d(3) + gl(2)*d_i(3) + gl(3)*d_j(3) + & |
| 395 |
gl(4)*(qjdi(3) - qidj(3)) + gl(5)*qir(3) + & |
| 396 |
gl(6)*qjr(3) + gl(7)*(qiqjr(3) + qjqir(3)) |
| 397 |
|
| 398 |
! |
| 399 |
! CALCULATE TORQUES BETWEEN MULTIPOLES I AND J |
| 400 |
! |
| 401 |
|
| 402 |
ttm_i(1) = - bn(2)*dixdj(1) + gl(2)*dixr(1) + gl(4)* & |
| 403 |
(dixqjr(1) + djxqir(1) + rxqidj(1) - 2.0*qixqj(1)) - & |
| 404 |
gl(5)*rxqir(1) - gl(7)*(rxqijr(1) + qjrxqir(1)) |
| 405 |
ttm_i(2) = - bn(2)*dixdj(2) + gl(2)*dixr(2) + gl(4)* & |
| 406 |
(dixqjr(2) + djxqir(2) + rxqidj(2) - 2.0*qixqj(2)) - & |
| 407 |
gl(5)*rxqir(2) - gl(7)*(rxqijr(2) + qjrxqir(2)) |
| 408 |
ttm_i(3) = - bn(2)*dixdj(3) + gl(2)*dixr(3) + gl(4)* & |
| 409 |
(dixqjr(3) + djxqir(3) + rxqidj(3) - 2.0*qixqj(3)) - & |
| 410 |
gl(5)*rxqir(3) - gl(7)*(rxqijr(3) + qjrxqir(3)) |
| 411 |
ttm_j(1) = bn(2)*dixdj(1) + gl(3)*djxr(1) - gl(4)* & |
| 412 |
(dixqjr(1) + djxqir(1) + rxqjdi(1) - 2.0*qixqj(1)) - & |
| 413 |
gl(6)*rxqjr(1) - gl(7)*(rxqjir(1) - qjrxqir(1)) |
| 414 |
ttm_j(2) = bn(2)*dixdj(2) + gl(3)*djxr(2) - gl(4)* & |
| 415 |
(dixqjr(2) + djxqir(2) + rxqjdi(2) - 2.0*qixqj(2)) - & |
| 416 |
gl(6)*rxqjr(2) - gl(7)*(rxqjir(2) - qjrxqir(2)) |
| 417 |
ttm_j(3) = bn(2)*dixdj(3) + gl(3)*djxr(3) - gl(4)* & |
| 418 |
(dixqjr(3) + djxqir(3) + rxqjdi(3) - 2.0*qixqj(3)) - & |
| 419 |
gl(6)*rxqjr(3) - gl(7)*(rxqjir(3) - qjrxqir(3)) |
| 420 |
|
| 421 |
ftm = ftm*sw |
| 422 |
ttm_i = ttm_i*sw |
| 423 |
ttm_j = ttm_j*sw |
| 424 |
|
| 425 |
#ifdef IS_MPI |
| 426 |
f_Row(1,atom1) = f_Row(1,atom1) + ftm(1) |
| 427 |
f_Row(2,atom1) = f_Row(2,atom1) + ftm(2) |
| 428 |
f_Row(3,atom1) = f_Row(3,atom1) + ftm(3) |
| 429 |
|
| 430 |
f_Col(1,atom2) = f_Col(1,atom2) - ftm(1) |
| 431 |
f_Col(2,atom2) = f_Col(2,atom2) - ftm(2) |
| 432 |
f_Col(3,atom2) = f_Col(3,atom2) - ftm(3) |
| 433 |
|
| 434 |
t_Row(1,atom1) = t_Row(1,atom1) + ttm_i(1) |
| 435 |
t_Row(2,atom1) = t_Row(2,atom1) + ttm_i(2) |
| 436 |
t_Row(3,atom1) = t_Row(3,atom1) + ttm_i(3) |
| 437 |
|
| 438 |
t_Col(1,atom2) = t_Col(1,atom2) + ttm_j(1) |
| 439 |
t_Col(2,atom2) = t_Col(2,atom2) + ttm_j(2) |
| 440 |
t_Col(3,atom2) = t_Col(3,atom2) + ttm_j(3) |
| 441 |
#else |
| 442 |
f(1,atom1) = f(1,atom1) + ftm(1) |
| 443 |
f(2,atom1) = f(2,atom1) + ftm(2) |
| 444 |
f(3,atom1) = f(3,atom1) + ftm(3) |
| 445 |
|
| 446 |
f(1,atom2) = f(1,atom2) - ftm(1) |
| 447 |
f(2,atom2) = f(2,atom2) - ftm(2) |
| 448 |
f(3,atom2) = f(3,atom2) - ftm(3) |
| 449 |
|
| 450 |
t(1,atom1) = t(1,atom1) + ttm_i(1) |
| 451 |
t(2,atom1) = t(2,atom1) + ttm_i(2) |
| 452 |
t(3,atom1) = t(3,atom1) + ttm_i(3) |
| 453 |
|
| 454 |
t(1,atom2) = t(1,atom2) + ttm_j(1) |
| 455 |
t(2,atom2) = t(2,atom2) + ttm_j(2) |
| 456 |
t(3,atom2) = t(3,atom2) + ttm_j(3) |
| 457 |
#endif |
| 458 |
|
| 459 |
#ifdef IS_MPI |
| 460 |
id1 = AtomRowToGlobal(atom1) |
| 461 |
id2 = AtomColToGlobal(atom2) |
| 462 |
#else |
| 463 |
id1 = atom1 |
| 464 |
id2 = atom2 |
| 465 |
#endif |
| 466 |
|
| 467 |
if (molMembershipList(id1) .ne. molMembershipList(id2)) then |
| 468 |
|
| 469 |
fpair(1) = fpair(1) + ftm(1) |
| 470 |
fpair(2) = fpair(2) + ftm(2) |
| 471 |
fpair(3) = fpair(3) + ftm(3) |
| 472 |
|
| 473 |
endif |
| 474 |
|
| 475 |
return |
| 476 |
|
| 477 |
end subroutine domultipolepair |
