1 |
!! |
2 |
!! Copyright (c) 2005 The University of Notre Dame. All Rights Reserved. |
3 |
!! |
4 |
!! The University of Notre Dame grants you ("Licensee") a |
5 |
!! non-exclusive, royalty free, license to use, modify and |
6 |
!! redistribute this software in source and binary code form, provided |
7 |
!! that the following conditions are met: |
8 |
!! |
9 |
!! 1. Acknowledgement of the program authors must be made in any |
10 |
!! publication of scientific results based in part on use of the |
11 |
!! program. An acceptable form of acknowledgement is citation of |
12 |
!! the article in which the program was described (Matthew |
13 |
!! A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher |
14 |
!! J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented |
15 |
!! Parallel Simulation Engine for Molecular Dynamics," |
16 |
!! J. Comput. Chem. 26, pp. 252-271 (2005)) |
17 |
!! |
18 |
!! 2. Redistributions of source code must retain the above copyright |
19 |
!! notice, this list of conditions and the following disclaimer. |
20 |
!! |
21 |
!! 3. Redistributions in binary form must reproduce the above copyright |
22 |
!! notice, this list of conditions and the following disclaimer in the |
23 |
!! documentation and/or other materials provided with the |
24 |
!! distribution. |
25 |
!! |
26 |
!! This software is provided "AS IS," without a warranty of any |
27 |
!! kind. All express or implied conditions, representations and |
28 |
!! warranties, including any implied warranty of merchantability, |
29 |
!! fitness for a particular purpose or non-infringement, are hereby |
30 |
!! excluded. The University of Notre Dame and its licensors shall not |
31 |
!! be liable for any damages suffered by licensee as a result of |
32 |
!! using, modifying or distributing the software or its |
33 |
!! derivatives. In no event will the University of Notre Dame or its |
34 |
!! licensors be liable for any lost revenue, profit or data, or for |
35 |
!! direct, indirect, special, consequential, incidental or punitive |
36 |
!! damages, however caused and regardless of the theory of liability, |
37 |
!! arising out of the use of or inability to use software, even if the |
38 |
!! University of Notre Dame has been advised of the possibility of |
39 |
!! such damages. |
40 |
!! |
41 |
|
42 |
module electrostatic_module |
43 |
|
44 |
use force_globals |
45 |
use definitions |
46 |
use atype_module |
47 |
use vector_class |
48 |
use simulation |
49 |
use status |
50 |
#ifdef IS_MPI |
51 |
use mpiSimulation |
52 |
#endif |
53 |
implicit none |
54 |
|
55 |
PRIVATE |
56 |
|
57 |
real(kind=dp), parameter :: pre11 = 332.0637778_dp |
58 |
real(kind=dp), parameter :: pre12 = 69.13291783_dp |
59 |
real(kind=dp), parameter :: pre22 = 14.39289874_dp |
60 |
real(kind=dp), parameter :: pre14 = 0.0_dp |
61 |
|
62 |
public :: newElectrostaticType |
63 |
public :: setCharge |
64 |
public :: setDipoleMoment |
65 |
public :: setSplitDipoleDistance |
66 |
public :: setQuadrupoleMoments |
67 |
public :: doElectrostaticPair |
68 |
public :: getCharge |
69 |
public :: getDipoleMoment |
70 |
|
71 |
type :: Electrostatic |
72 |
integer :: c_ident |
73 |
logical :: is_Charge = .false. |
74 |
logical :: is_Dipole = .false. |
75 |
logical :: is_SplitDipole = .false. |
76 |
logical :: is_Quadrupole = .false. |
77 |
real(kind=DP) :: charge = 0.0_DP |
78 |
real(kind=DP) :: dipole_moment = 0.0_DP |
79 |
real(kind=DP) :: split_dipole_distance = 0.0_DP |
80 |
real(kind=DP), dimension(3) :: quadrupole_moments = 0.0_DP |
81 |
end type Electrostatic |
82 |
|
83 |
type(Electrostatic), dimension(:), allocatable :: ElectrostaticMap |
84 |
|
85 |
contains |
86 |
|
87 |
subroutine newElectrostaticType(c_ident, is_Charge, is_Dipole, & |
88 |
is_SplitDipole, is_Quadrupole, status) |
89 |
|
90 |
integer, intent(in) :: c_ident |
91 |
logical, intent(in) :: is_Charge |
92 |
logical, intent(in) :: is_Dipole |
93 |
logical, intent(in) :: is_SplitDipole |
94 |
logical, intent(in) :: is_Quadrupole |
95 |
integer, intent(out) :: status |
96 |
integer :: nAtypes, myATID, i, j |
97 |
|
98 |
status = 0 |
99 |
myATID = getFirstMatchingElement(atypes, "c_ident", c_ident) |
100 |
|
101 |
!! Be simple-minded and assume that we need an ElectrostaticMap that |
102 |
!! is the same size as the total number of atom types |
103 |
|
104 |
if (.not.allocated(ElectrostaticMap)) then |
105 |
|
106 |
nAtypes = getSize(atypes) |
107 |
|
108 |
if (nAtypes == 0) then |
109 |
status = -1 |
110 |
return |
111 |
end if |
112 |
|
113 |
if (.not. allocated(ElectrostaticMap)) then |
114 |
allocate(ElectrostaticMap(nAtypes)) |
115 |
endif |
116 |
|
117 |
end if |
118 |
|
119 |
if (myATID .gt. size(ElectrostaticMap)) then |
120 |
status = -1 |
121 |
return |
122 |
endif |
123 |
|
124 |
! set the values for ElectrostaticMap for this atom type: |
125 |
|
126 |
ElectrostaticMap(myATID)%c_ident = c_ident |
127 |
ElectrostaticMap(myATID)%is_Charge = is_Charge |
128 |
ElectrostaticMap(myATID)%is_Dipole = is_Dipole |
129 |
ElectrostaticMap(myATID)%is_SplitDipole = is_SplitDipole |
130 |
ElectrostaticMap(myATID)%is_Quadrupole = is_Quadrupole |
131 |
|
132 |
end subroutine newElectrostaticType |
133 |
|
134 |
subroutine setCharge(c_ident, charge, status) |
135 |
integer, intent(in) :: c_ident |
136 |
real(kind=dp), intent(in) :: charge |
137 |
integer, intent(out) :: status |
138 |
integer :: myATID |
139 |
|
140 |
status = 0 |
141 |
myATID = getFirstMatchingElement(atypes, "c_ident", c_ident) |
142 |
|
143 |
if (.not.allocated(ElectrostaticMap)) then |
144 |
call handleError("electrostatic", "no ElectrostaticMap was present before first call of setCharge!") |
145 |
status = -1 |
146 |
return |
147 |
end if |
148 |
|
149 |
if (myATID .gt. size(ElectrostaticMap)) then |
150 |
call handleError("electrostatic", "ElectrostaticMap was found to be too small during setCharge!") |
151 |
status = -1 |
152 |
return |
153 |
endif |
154 |
|
155 |
if (.not.ElectrostaticMap(myATID)%is_Charge) then |
156 |
call handleError("electrostatic", "Attempt to setCharge of an atom type that is not a charge!") |
157 |
status = -1 |
158 |
return |
159 |
endif |
160 |
|
161 |
ElectrostaticMap(myATID)%charge = charge |
162 |
end subroutine setCharge |
163 |
|
164 |
subroutine setDipoleMoment(c_ident, dipole_moment, status) |
165 |
integer, intent(in) :: c_ident |
166 |
real(kind=dp), intent(in) :: dipole_moment |
167 |
integer, intent(out) :: status |
168 |
integer :: myATID |
169 |
|
170 |
status = 0 |
171 |
myATID = getFirstMatchingElement(atypes, "c_ident", c_ident) |
172 |
|
173 |
if (.not.allocated(ElectrostaticMap)) then |
174 |
call handleError("electrostatic", "no ElectrostaticMap was present before first call of setDipoleMoment!") |
175 |
status = -1 |
176 |
return |
177 |
end if |
178 |
|
179 |
if (myATID .gt. size(ElectrostaticMap)) then |
180 |
call handleError("electrostatic", "ElectrostaticMap was found to be too small during setDipoleMoment!") |
181 |
status = -1 |
182 |
return |
183 |
endif |
184 |
|
185 |
if (.not.ElectrostaticMap(myATID)%is_Dipole) then |
186 |
call handleError("electrostatic", "Attempt to setDipoleMoment of an atom type that is not a dipole!") |
187 |
status = -1 |
188 |
return |
189 |
endif |
190 |
|
191 |
ElectrostaticMap(myATID)%dipole_moment = dipole_moment |
192 |
end subroutine setDipoleMoment |
193 |
|
194 |
subroutine setSplitDipoleDistance(c_ident, split_dipole_distance, status) |
195 |
integer, intent(in) :: c_ident |
196 |
real(kind=dp), intent(in) :: split_dipole_distance |
197 |
integer, intent(out) :: status |
198 |
integer :: myATID |
199 |
|
200 |
status = 0 |
201 |
myATID = getFirstMatchingElement(atypes, "c_ident", c_ident) |
202 |
|
203 |
if (.not.allocated(ElectrostaticMap)) then |
204 |
call handleError("electrostatic", "no ElectrostaticMap was present before first call of setSplitDipoleDistance!") |
205 |
status = -1 |
206 |
return |
207 |
end if |
208 |
|
209 |
if (myATID .gt. size(ElectrostaticMap)) then |
210 |
call handleError("electrostatic", "ElectrostaticMap was found to be too small during setSplitDipoleDistance!") |
211 |
status = -1 |
212 |
return |
213 |
endif |
214 |
|
215 |
if (.not.ElectrostaticMap(myATID)%is_SplitDipole) then |
216 |
call handleError("electrostatic", "Attempt to setSplitDipoleDistance of an atom type that is not a splitDipole!") |
217 |
status = -1 |
218 |
return |
219 |
endif |
220 |
|
221 |
ElectrostaticMap(myATID)%split_dipole_distance = split_dipole_distance |
222 |
end subroutine setSplitDipoleDistance |
223 |
|
224 |
subroutine setQuadrupoleMoments(c_ident, quadrupole_moments, status) |
225 |
integer, intent(in) :: c_ident |
226 |
real(kind=dp), intent(in), dimension(3) :: quadrupole_moments |
227 |
integer, intent(out) :: status |
228 |
integer :: myATID, i, j |
229 |
|
230 |
status = 0 |
231 |
myATID = getFirstMatchingElement(atypes, "c_ident", c_ident) |
232 |
|
233 |
if (.not.allocated(ElectrostaticMap)) then |
234 |
call handleError("electrostatic", "no ElectrostaticMap was present before first call of setQuadrupoleMoments!") |
235 |
status = -1 |
236 |
return |
237 |
end if |
238 |
|
239 |
if (myATID .gt. size(ElectrostaticMap)) then |
240 |
call handleError("electrostatic", "ElectrostaticMap was found to be too small during setQuadrupoleMoments!") |
241 |
status = -1 |
242 |
return |
243 |
endif |
244 |
|
245 |
if (.not.ElectrostaticMap(myATID)%is_Quadrupole) then |
246 |
call handleError("electrostatic", "Attempt to setQuadrupoleMoments of an atom type that is not a quadrupole!") |
247 |
status = -1 |
248 |
return |
249 |
endif |
250 |
|
251 |
do i = 1, 3 |
252 |
ElectrostaticMap(myATID)%quadrupole_moments(i) = & |
253 |
quadrupole_moments(i) |
254 |
enddo |
255 |
|
256 |
end subroutine setQuadrupoleMoments |
257 |
|
258 |
|
259 |
function getCharge(atid) result (c) |
260 |
integer, intent(in) :: atid |
261 |
integer :: localError |
262 |
real(kind=dp) :: c |
263 |
|
264 |
if (.not.allocated(ElectrostaticMap)) then |
265 |
call handleError("electrostatic", "no ElectrostaticMap was present before first call of getCharge!") |
266 |
return |
267 |
end if |
268 |
|
269 |
if (.not.ElectrostaticMap(atid)%is_Charge) then |
270 |
call handleError("electrostatic", "getCharge was called for an atom type that isn't a charge!") |
271 |
return |
272 |
endif |
273 |
|
274 |
c = ElectrostaticMap(atid)%charge |
275 |
end function getCharge |
276 |
|
277 |
function getDipoleMoment(atid) result (dm) |
278 |
integer, intent(in) :: atid |
279 |
integer :: localError |
280 |
real(kind=dp) :: dm |
281 |
|
282 |
if (.not.allocated(ElectrostaticMap)) then |
283 |
call handleError("electrostatic", "no ElectrostaticMap was present before first call of getDipoleMoment!") |
284 |
return |
285 |
end if |
286 |
|
287 |
if (.not.ElectrostaticMap(atid)%is_Dipole) then |
288 |
call handleError("electrostatic", "getDipoleMoment was called for an atom type that isn't a dipole!") |
289 |
return |
290 |
endif |
291 |
|
292 |
dm = ElectrostaticMap(atid)%dipole_moment |
293 |
end function getDipoleMoment |
294 |
|
295 |
subroutine doElectrostaticPair(atom1, atom2, d, rij, r2, sw, & |
296 |
vpair, fpair, pot, eFrame, f, t, do_pot) |
297 |
|
298 |
logical, intent(in) :: do_pot |
299 |
|
300 |
integer, intent(in) :: atom1, atom2 |
301 |
integer :: localError |
302 |
|
303 |
real(kind=dp), intent(in) :: rij, r2, sw |
304 |
real(kind=dp), intent(in), dimension(3) :: d |
305 |
real(kind=dp), intent(inout) :: vpair |
306 |
real(kind=dp), intent(inout), dimension(3) :: fpair |
307 |
|
308 |
real( kind = dp ) :: pot |
309 |
real( kind = dp ), dimension(9,nLocal) :: eFrame |
310 |
real( kind = dp ), dimension(3,nLocal) :: f |
311 |
real( kind = dp ), dimension(3,nLocal) :: t |
312 |
|
313 |
real (kind = dp), dimension(3) :: ul_i |
314 |
real (kind = dp), dimension(3) :: ul_j |
315 |
|
316 |
logical :: i_is_Charge, i_is_Dipole, i_is_SplitDipole, i_is_Quadrupole |
317 |
logical :: j_is_Charge, j_is_Dipole, j_is_SplitDipole, j_is_Quadrupole |
318 |
integer :: me1, me2, id1, id2 |
319 |
real (kind=dp) :: q_i, q_j, mu_i, mu_j, d_i, d_j |
320 |
real (kind=dp) :: ct_i, ct_j, ct_ij, a1 |
321 |
real (kind=dp) :: riji, ri, ri2, ri3, ri4 |
322 |
real (kind=dp) :: pref, vterm, epot, dudr |
323 |
real (kind=dp) :: xhat, yhat, zhat |
324 |
real (kind=dp) :: dudx, dudy, dudz |
325 |
real (kind=dp) :: drdxj, drdyj, drdzj |
326 |
real (kind=dp) :: duduix, duduiy, duduiz, dudujx, dudujy, dudujz |
327 |
real (kind=dp) :: scale, sc2, bigR |
328 |
|
329 |
if (.not.allocated(ElectrostaticMap)) then |
330 |
call handleError("electrostatic", "no ElectrostaticMap was present before first call of do_electrostatic_pair!") |
331 |
return |
332 |
end if |
333 |
|
334 |
#ifdef IS_MPI |
335 |
me1 = atid_Row(atom1) |
336 |
me2 = atid_Col(atom2) |
337 |
#else |
338 |
me1 = atid(atom1) |
339 |
me2 = atid(atom2) |
340 |
#endif |
341 |
|
342 |
!! some variables we'll need independent of electrostatic type: |
343 |
|
344 |
riji = 1.0d0 / rij |
345 |
|
346 |
xhat = d(1) * riji |
347 |
yhat = d(2) * riji |
348 |
zhat = d(3) * riji |
349 |
|
350 |
drdxj = xhat |
351 |
drdyj = yhat |
352 |
drdzj = zhat |
353 |
|
354 |
!! logicals |
355 |
|
356 |
i_is_Charge = ElectrostaticMap(me1)%is_Charge |
357 |
i_is_Dipole = ElectrostaticMap(me1)%is_Dipole |
358 |
i_is_SplitDipole = ElectrostaticMap(me1)%is_SplitDipole |
359 |
i_is_Quadrupole = ElectrostaticMap(me1)%is_Quadrupole |
360 |
|
361 |
j_is_Charge = ElectrostaticMap(me2)%is_Charge |
362 |
j_is_Dipole = ElectrostaticMap(me2)%is_Dipole |
363 |
j_is_SplitDipole = ElectrostaticMap(me2)%is_SplitDipole |
364 |
j_is_Quadrupole = ElectrostaticMap(me2)%is_Quadrupole |
365 |
|
366 |
if (i_is_Charge) then |
367 |
q_i = ElectrostaticMap(me1)%charge |
368 |
endif |
369 |
|
370 |
if (i_is_Dipole) then |
371 |
mu_i = ElectrostaticMap(me1)%dipole_moment |
372 |
#ifdef IS_MPI |
373 |
ul_i(1) = eFrame_Row(3,atom1) |
374 |
ul_i(2) = eFrame_Row(6,atom1) |
375 |
ul_i(3) = eFrame_Row(9,atom1) |
376 |
#else |
377 |
ul_i(1) = eFrame(3,atom1) |
378 |
ul_i(2) = eFrame(6,atom1) |
379 |
ul_i(3) = eFrame(9,atom1) |
380 |
#endif |
381 |
ct_i = ul_i(1)*drdxj + ul_i(2)*drdyj + ul_i(3)*drdzj |
382 |
|
383 |
if (i_is_SplitDipole) then |
384 |
d_i = ElectrostaticMap(me1)%split_dipole_distance |
385 |
endif |
386 |
|
387 |
endif |
388 |
|
389 |
if (j_is_Charge) then |
390 |
q_j = ElectrostaticMap(me2)%charge |
391 |
endif |
392 |
|
393 |
if (j_is_Dipole) then |
394 |
mu_j = ElectrostaticMap(me2)%dipole_moment |
395 |
#ifdef IS_MPI |
396 |
ul_j(1) = eFrame_Col(3,atom2) |
397 |
ul_j(2) = eFrame_Col(6,atom2) |
398 |
ul_j(3) = eFrame_Col(9,atom2) |
399 |
#else |
400 |
ul_j(1) = eFrame(3,atom2) |
401 |
ul_j(2) = eFrame(6,atom2) |
402 |
ul_j(3) = eFrame(9,atom2) |
403 |
#endif |
404 |
ct_j = ul_j(1)*drdxj + ul_j(2)*drdyj + ul_j(3)*drdzj |
405 |
|
406 |
if (j_is_SplitDipole) then |
407 |
d_j = ElectrostaticMap(me2)%split_dipole_distance |
408 |
endif |
409 |
endif |
410 |
|
411 |
epot = 0.0_dp |
412 |
dudx = 0.0_dp |
413 |
dudy = 0.0_dp |
414 |
dudz = 0.0_dp |
415 |
|
416 |
duduix = 0.0_dp |
417 |
duduiy = 0.0_dp |
418 |
duduiz = 0.0_dp |
419 |
|
420 |
dudujx = 0.0_dp |
421 |
dudujy = 0.0_dp |
422 |
dudujz = 0.0_dp |
423 |
|
424 |
if (i_is_Charge) then |
425 |
|
426 |
if (j_is_Charge) then |
427 |
|
428 |
vterm = pre11 * q_i * q_j * riji |
429 |
vpair = vpair + vterm |
430 |
epot = epot + sw*vterm |
431 |
|
432 |
dudr = - sw * vterm * riji |
433 |
|
434 |
dudx = dudx + dudr * drdxj |
435 |
dudy = dudy + dudr * drdyj |
436 |
dudz = dudz + dudr * drdzj |
437 |
|
438 |
endif |
439 |
|
440 |
if (j_is_Dipole) then |
441 |
|
442 |
if (j_is_SplitDipole) then |
443 |
BigR = sqrt(r2 + 0.25_dp * d_j * d_j) |
444 |
ri = 1.0_dp / BigR |
445 |
scale = rij * ri |
446 |
else |
447 |
ri = riji |
448 |
scale = 1.0_dp |
449 |
endif |
450 |
|
451 |
ri2 = ri * ri |
452 |
ri3 = ri2 * ri |
453 |
sc2 = scale * scale |
454 |
|
455 |
pref = pre12 * q_i * mu_j |
456 |
vterm = pref * ct_j * ri2 * scale |
457 |
vpair = vpair + vterm |
458 |
epot = epot + sw * vterm |
459 |
|
460 |
!! this has a + sign in the () because the rij vector is |
461 |
!! r_j - r_i and the charge-dipole potential takes the origin |
462 |
!! as the point dipole, which is atom j in this case. |
463 |
|
464 |
dudx = dudx + pref * sw * ri3 * ( ul_j(1) + 3.0d0*ct_j*xhat*sc2) |
465 |
dudy = dudy + pref * sw * ri3 * ( ul_j(2) + 3.0d0*ct_j*yhat*sc2) |
466 |
dudz = dudz + pref * sw * ri3 * ( ul_j(3) + 3.0d0*ct_j*zhat*sc2) |
467 |
|
468 |
dudujx = dudujx - pref * sw * ri2 * xhat * scale |
469 |
dudujy = dudujy - pref * sw * ri2 * yhat * scale |
470 |
dudujz = dudujz - pref * sw * ri2 * zhat * scale |
471 |
|
472 |
endif |
473 |
|
474 |
endif |
475 |
|
476 |
if (i_is_Dipole) then |
477 |
|
478 |
if (j_is_Charge) then |
479 |
|
480 |
if (i_is_SplitDipole) then |
481 |
BigR = sqrt(r2 + 0.25_dp * d_i * d_i) |
482 |
ri = 1.0_dp / BigR |
483 |
scale = rij * ri |
484 |
else |
485 |
ri = riji |
486 |
scale = 1.0_dp |
487 |
endif |
488 |
|
489 |
ri2 = ri * ri |
490 |
ri3 = ri2 * ri |
491 |
sc2 = scale * scale |
492 |
|
493 |
pref = pre12 * q_j * mu_i |
494 |
vterm = pref * ct_i * ri2 * scale |
495 |
vpair = vpair + vterm |
496 |
epot = epot + sw * vterm |
497 |
|
498 |
dudx = dudx + pref * sw * ri3 * ( ul_i(1) - 3.0d0 * ct_i * xhat*sc2) |
499 |
dudy = dudy + pref * sw * ri3 * ( ul_i(2) - 3.0d0 * ct_i * yhat*sc2) |
500 |
dudz = dudz + pref * sw * ri3 * ( ul_i(3) - 3.0d0 * ct_i * zhat*sc2) |
501 |
|
502 |
duduix = duduix + pref * sw * ri2 * xhat * scale |
503 |
duduiy = duduiy + pref * sw * ri2 * yhat * scale |
504 |
duduiz = duduiz + pref * sw * ri2 * zhat * scale |
505 |
endif |
506 |
|
507 |
if (j_is_Dipole) then |
508 |
|
509 |
if (i_is_SplitDipole) then |
510 |
if (j_is_SplitDipole) then |
511 |
BigR = sqrt(r2 + 0.25_dp * d_i * d_i + 0.25_dp * d_j * d_j) |
512 |
else |
513 |
BigR = sqrt(r2 + 0.25_dp * d_i * d_i) |
514 |
endif |
515 |
ri = 1.0_dp / BigR |
516 |
scale = rij * ri |
517 |
else |
518 |
if (j_is_SplitDipole) then |
519 |
BigR = sqrt(r2 + 0.25_dp * d_j * d_j) |
520 |
ri = 1.0_dp / BigR |
521 |
scale = rij * ri |
522 |
else |
523 |
ri = riji |
524 |
scale = 1.0_dp |
525 |
endif |
526 |
endif |
527 |
|
528 |
ct_ij = ul_i(1)*ul_j(1) + ul_i(2)*ul_j(2) + ul_i(3)*ul_j(3) |
529 |
|
530 |
ri2 = ri * ri |
531 |
ri3 = ri2 * ri |
532 |
ri4 = ri2 * ri2 |
533 |
sc2 = scale * scale |
534 |
|
535 |
pref = pre22 * mu_i * mu_j |
536 |
vterm = pref * ri3 * (ct_ij - 3.0d0 * ct_i * ct_j * sc2) |
537 |
vpair = vpair + vterm |
538 |
epot = epot + sw * vterm |
539 |
|
540 |
a1 = 5.0d0 * ct_i * ct_j * sc2 - ct_ij |
541 |
|
542 |
dudx=dudx+pref*sw*3.0d0*ri4*scale*(a1*xhat-ct_i*ul_j(1)-ct_j*ul_i(1)) |
543 |
dudy=dudy+pref*sw*3.0d0*ri4*scale*(a1*yhat-ct_i*ul_j(2)-ct_j*ul_i(2)) |
544 |
dudz=dudz+pref*sw*3.0d0*ri4*scale*(a1*zhat-ct_i*ul_j(3)-ct_j*ul_i(3)) |
545 |
|
546 |
duduix = duduix + pref*sw*ri3*(ul_j(1) - 3.0d0*ct_j*xhat*sc2) |
547 |
duduiy = duduiy + pref*sw*ri3*(ul_j(2) - 3.0d0*ct_j*yhat*sc2) |
548 |
duduiz = duduiz + pref*sw*ri3*(ul_j(3) - 3.0d0*ct_j*zhat*sc2) |
549 |
|
550 |
dudujx = dudujx + pref*sw*ri3*(ul_i(1) - 3.0d0*ct_i*xhat*sc2) |
551 |
dudujy = dudujy + pref*sw*ri3*(ul_i(2) - 3.0d0*ct_i*yhat*sc2) |
552 |
dudujz = dudujz + pref*sw*ri3*(ul_i(3) - 3.0d0*ct_i*zhat*sc2) |
553 |
endif |
554 |
|
555 |
endif |
556 |
|
557 |
if (do_pot) then |
558 |
#ifdef IS_MPI |
559 |
pot_row(atom1) = pot_row(atom1) + 0.5d0*epot |
560 |
pot_col(atom2) = pot_col(atom2) + 0.5d0*epot |
561 |
#else |
562 |
pot = pot + epot |
563 |
#endif |
564 |
endif |
565 |
|
566 |
#ifdef IS_MPI |
567 |
f_Row(1,atom1) = f_Row(1,atom1) + dudx |
568 |
f_Row(2,atom1) = f_Row(2,atom1) + dudy |
569 |
f_Row(3,atom1) = f_Row(3,atom1) + dudz |
570 |
|
571 |
f_Col(1,atom2) = f_Col(1,atom2) - dudx |
572 |
f_Col(2,atom2) = f_Col(2,atom2) - dudy |
573 |
f_Col(3,atom2) = f_Col(3,atom2) - dudz |
574 |
|
575 |
if (i_is_Dipole .or. i_is_Quadrupole) then |
576 |
t_Row(1,atom1) = t_Row(1,atom1) - ul_i(2)*duduiz + ul_i(3)*duduiy |
577 |
t_Row(2,atom1) = t_Row(2,atom1) - ul_i(3)*duduix + ul_i(1)*duduiz |
578 |
t_Row(3,atom1) = t_Row(3,atom1) - ul_i(1)*duduiy + ul_i(2)*duduix |
579 |
endif |
580 |
|
581 |
if (j_is_Dipole .or. j_is_Quadrupole) then |
582 |
t_Col(1,atom2) = t_Col(1,atom2) - ul_j(2)*dudujz + ul_j(3)*dudujy |
583 |
t_Col(2,atom2) = t_Col(2,atom2) - ul_j(3)*dudujx + ul_j(1)*dudujz |
584 |
t_Col(3,atom2) = t_Col(3,atom2) - ul_j(1)*dudujy + ul_j(2)*dudujx |
585 |
endif |
586 |
|
587 |
#else |
588 |
f(1,atom1) = f(1,atom1) + dudx |
589 |
f(2,atom1) = f(2,atom1) + dudy |
590 |
f(3,atom1) = f(3,atom1) + dudz |
591 |
|
592 |
f(1,atom2) = f(1,atom2) - dudx |
593 |
f(2,atom2) = f(2,atom2) - dudy |
594 |
f(3,atom2) = f(3,atom2) - dudz |
595 |
|
596 |
if (i_is_Dipole .or. i_is_Quadrupole) then |
597 |
t(1,atom1) = t(1,atom1) - ul_i(2)*duduiz + ul_i(3)*duduiy |
598 |
t(2,atom1) = t(2,atom1) - ul_i(3)*duduix + ul_i(1)*duduiz |
599 |
t(3,atom1) = t(3,atom1) - ul_i(1)*duduiy + ul_i(2)*duduix |
600 |
endif |
601 |
|
602 |
if (j_is_Dipole .or. j_is_Quadrupole) then |
603 |
t(1,atom2) = t(1,atom2) - ul_j(2)*dudujz + ul_j(3)*dudujy |
604 |
t(2,atom2) = t(2,atom2) - ul_j(3)*dudujx + ul_j(1)*dudujz |
605 |
t(3,atom2) = t(3,atom2) - ul_j(1)*dudujy + ul_j(2)*dudujx |
606 |
endif |
607 |
#endif |
608 |
|
609 |
#ifdef IS_MPI |
610 |
id1 = AtomRowToGlobal(atom1) |
611 |
id2 = AtomColToGlobal(atom2) |
612 |
#else |
613 |
id1 = atom1 |
614 |
id2 = atom2 |
615 |
#endif |
616 |
|
617 |
if (molMembershipList(id1) .ne. molMembershipList(id2)) then |
618 |
|
619 |
fpair(1) = fpair(1) + dudx |
620 |
fpair(2) = fpair(2) + dudy |
621 |
fpair(3) = fpair(3) + dudz |
622 |
|
623 |
endif |
624 |
|
625 |
return |
626 |
end subroutine doElectrostaticPair |
627 |
|
628 |
end module electrostatic_module |