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 |
|
61 |
public :: newElectrostaticType |
62 |
public :: setCharge |
63 |
public :: setDipoleMoment |
64 |
public :: setSplitDipoleDistance |
65 |
public :: setQuadrupoleMoments |
66 |
public :: doElectrostaticPair |
67 |
public :: getCharge |
68 |
public :: getDipoleMoment |
69 |
|
70 |
type :: Electrostatic |
71 |
integer :: c_ident |
72 |
logical :: is_Charge = .false. |
73 |
logical :: is_Dipole = .false. |
74 |
logical :: is_SplitDipole = .false. |
75 |
logical :: is_Quadrupole = .false. |
76 |
real(kind=DP) :: charge = 0.0_DP |
77 |
real(kind=DP) :: dipole_moment = 0.0_DP |
78 |
real(kind=DP) :: split_dipole_distance = 0.0_DP |
79 |
real(kind=DP), dimension(3) :: quadrupole_moments = 0.0_DP |
80 |
end type Electrostatic |
81 |
|
82 |
type(Electrostatic), dimension(:), allocatable :: ElectrostaticMap |
83 |
|
84 |
contains |
85 |
|
86 |
subroutine newElectrostaticType(c_ident, is_Charge, is_Dipole, & |
87 |
is_SplitDipole, is_Quadrupole, status) |
88 |
|
89 |
integer, intent(in) :: c_ident |
90 |
logical, intent(in) :: is_Charge |
91 |
logical, intent(in) :: is_Dipole |
92 |
logical, intent(in) :: is_SplitDipole |
93 |
logical, intent(in) :: is_Quadrupole |
94 |
integer, intent(out) :: status |
95 |
integer :: nAtypes, myATID, i, j |
96 |
|
97 |
status = 0 |
98 |
myATID = getFirstMatchingElement(atypes, "c_ident", c_ident) |
99 |
|
100 |
!! Be simple-minded and assume that we need an ElectrostaticMap that |
101 |
!! is the same size as the total number of atom types |
102 |
|
103 |
if (.not.allocated(ElectrostaticMap)) then |
104 |
|
105 |
nAtypes = getSize(atypes) |
106 |
|
107 |
if (nAtypes == 0) then |
108 |
status = -1 |
109 |
return |
110 |
end if |
111 |
|
112 |
if (.not. allocated(ElectrostaticMap)) then |
113 |
allocate(ElectrostaticMap(nAtypes)) |
114 |
endif |
115 |
|
116 |
end if |
117 |
|
118 |
if (myATID .gt. size(ElectrostaticMap)) then |
119 |
status = -1 |
120 |
return |
121 |
endif |
122 |
|
123 |
! set the values for ElectrostaticMap for this atom type: |
124 |
|
125 |
ElectrostaticMap(myATID)%c_ident = c_ident |
126 |
ElectrostaticMap(myATID)%is_Charge = is_Charge |
127 |
ElectrostaticMap(myATID)%is_Dipole = is_Dipole |
128 |
ElectrostaticMap(myATID)%is_SplitDipole = is_SplitDipole |
129 |
ElectrostaticMap(myATID)%is_Quadrupole = is_Quadrupole |
130 |
|
131 |
end subroutine newElectrostaticType |
132 |
|
133 |
subroutine setCharge(c_ident, charge, status) |
134 |
integer, intent(in) :: c_ident |
135 |
real(kind=dp), intent(in) :: charge |
136 |
integer, intent(out) :: status |
137 |
integer :: myATID |
138 |
|
139 |
status = 0 |
140 |
myATID = getFirstMatchingElement(atypes, "c_ident", c_ident) |
141 |
|
142 |
if (.not.allocated(ElectrostaticMap)) then |
143 |
call handleError("electrostatic", "no ElectrostaticMap was present before first call of setCharge!") |
144 |
status = -1 |
145 |
return |
146 |
end if |
147 |
|
148 |
if (myATID .gt. size(ElectrostaticMap)) then |
149 |
call handleError("electrostatic", "ElectrostaticMap was found to be too small during setCharge!") |
150 |
status = -1 |
151 |
return |
152 |
endif |
153 |
|
154 |
if (.not.ElectrostaticMap(myATID)%is_Charge) then |
155 |
call handleError("electrostatic", "Attempt to setCharge of an atom type that is not a charge!") |
156 |
status = -1 |
157 |
return |
158 |
endif |
159 |
|
160 |
ElectrostaticMap(myATID)%charge = charge |
161 |
end subroutine setCharge |
162 |
|
163 |
subroutine setDipoleMoment(c_ident, dipole_moment, status) |
164 |
integer, intent(in) :: c_ident |
165 |
real(kind=dp), intent(in) :: dipole_moment |
166 |
integer, intent(out) :: status |
167 |
integer :: myATID |
168 |
|
169 |
status = 0 |
170 |
myATID = getFirstMatchingElement(atypes, "c_ident", c_ident) |
171 |
|
172 |
if (.not.allocated(ElectrostaticMap)) then |
173 |
call handleError("electrostatic", "no ElectrostaticMap was present before first call of setDipoleMoment!") |
174 |
status = -1 |
175 |
return |
176 |
end if |
177 |
|
178 |
if (myATID .gt. size(ElectrostaticMap)) then |
179 |
call handleError("electrostatic", "ElectrostaticMap was found to be too small during setDipoleMoment!") |
180 |
status = -1 |
181 |
return |
182 |
endif |
183 |
|
184 |
if (.not.ElectrostaticMap(myATID)%is_Dipole) then |
185 |
call handleError("electrostatic", "Attempt to setDipoleMoment of an atom type that is not a dipole!") |
186 |
status = -1 |
187 |
return |
188 |
endif |
189 |
|
190 |
ElectrostaticMap(myATID)%dipole_moment = dipole_moment |
191 |
end subroutine setDipoleMoment |
192 |
|
193 |
subroutine setSplitDipoleDistance(c_ident, split_dipole_distance, status) |
194 |
integer, intent(in) :: c_ident |
195 |
real(kind=dp), intent(in) :: split_dipole_distance |
196 |
integer, intent(out) :: status |
197 |
integer :: myATID |
198 |
|
199 |
status = 0 |
200 |
myATID = getFirstMatchingElement(atypes, "c_ident", c_ident) |
201 |
|
202 |
if (.not.allocated(ElectrostaticMap)) then |
203 |
call handleError("electrostatic", "no ElectrostaticMap was present before first call of setSplitDipoleDistance!") |
204 |
status = -1 |
205 |
return |
206 |
end if |
207 |
|
208 |
if (myATID .gt. size(ElectrostaticMap)) then |
209 |
call handleError("electrostatic", "ElectrostaticMap was found to be too small during setSplitDipoleDistance!") |
210 |
status = -1 |
211 |
return |
212 |
endif |
213 |
|
214 |
if (.not.ElectrostaticMap(myATID)%is_SplitDipole) then |
215 |
call handleError("electrostatic", "Attempt to setSplitDipoleDistance of an atom type that is not a splitDipole!") |
216 |
status = -1 |
217 |
return |
218 |
endif |
219 |
|
220 |
ElectrostaticMap(myATID)%split_dipole_distance = split_dipole_distance |
221 |
end subroutine setSplitDipoleDistance |
222 |
|
223 |
subroutine setQuadrupoleMoments(c_ident, quadrupole_moments, status) |
224 |
integer, intent(in) :: c_ident |
225 |
real(kind=dp), intent(in), dimension(3) :: quadrupole_moments |
226 |
integer, intent(out) :: status |
227 |
integer :: myATID, i, j |
228 |
|
229 |
status = 0 |
230 |
myATID = getFirstMatchingElement(atypes, "c_ident", c_ident) |
231 |
|
232 |
if (.not.allocated(ElectrostaticMap)) then |
233 |
call handleError("electrostatic", "no ElectrostaticMap was present before first call of setQuadrupoleMoments!") |
234 |
status = -1 |
235 |
return |
236 |
end if |
237 |
|
238 |
if (myATID .gt. size(ElectrostaticMap)) then |
239 |
call handleError("electrostatic", "ElectrostaticMap was found to be too small during setQuadrupoleMoments!") |
240 |
status = -1 |
241 |
return |
242 |
endif |
243 |
|
244 |
if (.not.ElectrostaticMap(myATID)%is_Quadrupole) then |
245 |
call handleError("electrostatic", "Attempt to setQuadrupoleMoments of an atom type that is not a quadrupole!") |
246 |
status = -1 |
247 |
return |
248 |
endif |
249 |
|
250 |
do i = 1, 3 |
251 |
ElectrostaticMap(myATID)%quadrupole_moments(i) = & |
252 |
quadrupole_moments(i) |
253 |
enddo |
254 |
|
255 |
end subroutine setQuadrupoleMoments |
256 |
|
257 |
|
258 |
function getCharge(atid) result (c) |
259 |
integer, intent(in) :: atid |
260 |
integer :: localError |
261 |
real(kind=dp) :: c |
262 |
|
263 |
if (.not.allocated(ElectrostaticMap)) then |
264 |
call handleError("electrostatic", "no ElectrostaticMap was present before first call of getCharge!") |
265 |
return |
266 |
end if |
267 |
|
268 |
if (.not.ElectrostaticMap(atid)%is_Charge) then |
269 |
call handleError("electrostatic", "getCharge was called for an atom type that isn't a charge!") |
270 |
return |
271 |
endif |
272 |
|
273 |
c = ElectrostaticMap(atid)%charge |
274 |
end function getCharge |
275 |
|
276 |
function getDipoleMoment(atid) result (dm) |
277 |
integer, intent(in) :: atid |
278 |
integer :: localError |
279 |
real(kind=dp) :: dm |
280 |
|
281 |
if (.not.allocated(ElectrostaticMap)) then |
282 |
call handleError("electrostatic", "no ElectrostaticMap was present before first call of getDipoleMoment!") |
283 |
return |
284 |
end if |
285 |
|
286 |
if (.not.ElectrostaticMap(atid)%is_Dipole) then |
287 |
call handleError("electrostatic", "getDipoleMoment was called for an atom type that isn't a dipole!") |
288 |
return |
289 |
endif |
290 |
|
291 |
dm = ElectrostaticMap(atid)%dipole_moment |
292 |
end function getDipoleMoment |
293 |
|
294 |
subroutine doElectrostaticPair(atom1, atom2, d, rij, r2, sw, & |
295 |
vpair, fpair, pot, eFrame, f, t, do_pot) |
296 |
|
297 |
logical, intent(in) :: do_pot |
298 |
|
299 |
integer, intent(in) :: atom1, atom2 |
300 |
integer :: localError |
301 |
|
302 |
real(kind=dp), intent(in) :: rij, r2, sw |
303 |
real(kind=dp), intent(in), dimension(3) :: d |
304 |
real(kind=dp), intent(inout) :: vpair |
305 |
real(kind=dp), intent(inout), dimension(3) :: fpair |
306 |
|
307 |
real( kind = dp ) :: pot |
308 |
real( kind = dp ), dimension(9,nLocal) :: eFrame |
309 |
real( kind = dp ), dimension(3,nLocal) :: f |
310 |
real( kind = dp ), dimension(3,nLocal) :: t |
311 |
|
312 |
real (kind = dp), dimension(3) :: ul_i |
313 |
real (kind = dp), dimension(3) :: ul_j |
314 |
|
315 |
logical :: i_is_Charge, i_is_Dipole, i_is_SplitDipole, i_is_Quadrupole |
316 |
logical :: j_is_Charge, j_is_Dipole, j_is_SplitDipole, j_is_Quadrupole |
317 |
integer :: me1, me2, id1, id2 |
318 |
real (kind=dp) :: q_i, q_j, mu_i, mu_j, d_i, d_j |
319 |
real (kind=dp) :: ct_i, ct_j, ct_ij, a1 |
320 |
real (kind=dp) :: riji, ri2, ri3, ri4 |
321 |
real (kind=dp) :: pref, vterm, epot, dudr |
322 |
real (kind=dp) :: dudx, dudy, dudz |
323 |
real (kind=dp) :: drdxj, drdyj, drdzj |
324 |
real (kind=dp) :: duduix, duduiy, duduiz, dudujx, dudujy, dudujz |
325 |
|
326 |
|
327 |
if (.not.allocated(ElectrostaticMap)) then |
328 |
call handleError("electrostatic", "no ElectrostaticMap was present before first call of do_electrostatic_pair!") |
329 |
return |
330 |
end if |
331 |
|
332 |
#ifdef IS_MPI |
333 |
me1 = atid_Row(atom1) |
334 |
me2 = atid_Col(atom2) |
335 |
#else |
336 |
me1 = atid(atom1) |
337 |
me2 = atid(atom2) |
338 |
#endif |
339 |
|
340 |
!! some variables we'll need independent of electrostatic type: |
341 |
|
342 |
riji = 1.0d0 / rij |
343 |
|
344 |
!! these are also useful as the unit vector of \vec{r} |
345 |
!! \hat{r} = \vec{r} / r = {(x_j-x_i) / r, (y_j-y_i)/r, (z_j-z_i)/r} |
346 |
|
347 |
drdxj = d(1) * riji |
348 |
drdyj = d(2) * riji |
349 |
drdzj = d(3) * riji |
350 |
|
351 |
!! logicals |
352 |
|
353 |
i_is_Charge = ElectrostaticMap(me1)%is_Charge |
354 |
i_is_Dipole = ElectrostaticMap(me1)%is_Dipole |
355 |
i_is_SplitDipole = ElectrostaticMap(me1)%is_SplitDipole |
356 |
i_is_Quadrupole = ElectrostaticMap(me1)%is_Quadrupole |
357 |
|
358 |
j_is_Charge = ElectrostaticMap(me2)%is_Charge |
359 |
j_is_Dipole = ElectrostaticMap(me2)%is_Dipole |
360 |
j_is_SplitDipole = ElectrostaticMap(me2)%is_SplitDipole |
361 |
j_is_Quadrupole = ElectrostaticMap(me2)%is_Quadrupole |
362 |
|
363 |
if (i_is_Charge) then |
364 |
q_i = ElectrostaticMap(me1)%charge |
365 |
endif |
366 |
|
367 |
if (i_is_Dipole) then |
368 |
mu_i = ElectrostaticMap(me1)%dipole_moment |
369 |
#ifdef IS_MPI |
370 |
ul_i(1) = eFrame_Row(3,atom1) |
371 |
ul_i(2) = eFrame_Row(6,atom1) |
372 |
ul_i(3) = eFrame_Row(9,atom1) |
373 |
#else |
374 |
ul_i(1) = eFrame(3,atom1) |
375 |
ul_i(2) = eFrame(6,atom1) |
376 |
ul_i(3) = eFrame(9,atom1) |
377 |
#endif |
378 |
ct_i = ul_i(1)*drdxj + ul_i(2)*drdyj + ul_i(3)*drdzj |
379 |
|
380 |
if (i_is_SplitDipole) then |
381 |
d_i = ElectrostaticMap(me1)%split_dipole_distance |
382 |
endif |
383 |
|
384 |
endif |
385 |
|
386 |
if (j_is_Charge) then |
387 |
q_j = ElectrostaticMap(me2)%charge |
388 |
endif |
389 |
|
390 |
if (j_is_Dipole) then |
391 |
mu_j = ElectrostaticMap(me2)%dipole_moment |
392 |
#ifdef IS_MPI |
393 |
ul_j(1) = eFrame_Col(3,atom2) |
394 |
ul_j(2) = eFrame_Col(6,atom2) |
395 |
ul_j(3) = eFrame_Col(9,atom2) |
396 |
#else |
397 |
ul_j(1) = eFrame(3,atom2) |
398 |
ul_j(2) = eFrame(6,atom2) |
399 |
ul_j(3) = eFrame(9,atom2) |
400 |
#endif |
401 |
ct_j = ul_j(1)*drdxj + ul_j(2)*drdyj + ul_j(3)*drdzj |
402 |
|
403 |
if (j_is_SplitDipole) then |
404 |
d_j = ElectrostaticMap(me2)%split_dipole_distance |
405 |
endif |
406 |
endif |
407 |
|
408 |
epot = 0.0_dp |
409 |
dudx = 0.0_dp |
410 |
dudy = 0.0_dp |
411 |
dudz = 0.0_dp |
412 |
|
413 |
duduix = 0.0_dp |
414 |
duduiy = 0.0_dp |
415 |
duduiz = 0.0_dp |
416 |
|
417 |
dudujx = 0.0_dp |
418 |
dudujy = 0.0_dp |
419 |
dudujz = 0.0_dp |
420 |
|
421 |
if (i_is_Charge) then |
422 |
|
423 |
if (j_is_Charge) then |
424 |
|
425 |
vterm = pre11 * q_i * q_j * riji |
426 |
vpair = vpair + vterm |
427 |
epot = epot + sw*vterm |
428 |
|
429 |
dudr = - sw * vterm * riji |
430 |
|
431 |
dudx = dudx + dudr * drdxj |
432 |
dudy = dudy + dudr * drdyj |
433 |
dudz = dudz + dudr * drdzj |
434 |
|
435 |
endif |
436 |
|
437 |
if (j_is_Dipole) then |
438 |
|
439 |
ri2 = riji * riji |
440 |
ri3 = ri2 * riji |
441 |
|
442 |
pref = pre12 * q_i * mu_j |
443 |
vterm = pref * ct_j * riji * riji |
444 |
vpair = vpair + vterm |
445 |
epot = epot + sw * vterm |
446 |
|
447 |
dudx = dudx + pref * sw * ri3 * ( ul_j(1) + 3.0d0 * ct_j * drdxj) |
448 |
dudy = dudy + pref * sw * ri3 * ( ul_j(2) + 3.0d0 * ct_j * drdyj) |
449 |
dudz = dudz + pref * sw * ri3 * ( ul_j(3) + 3.0d0 * ct_j * drdzj) |
450 |
|
451 |
dudujx = dudujx - pref * sw * ri2 * drdxj |
452 |
dudujy = dudujy - pref * sw * ri2 * drdyj |
453 |
dudujz = dudujz - pref * sw * ri2 * drdzj |
454 |
|
455 |
endif |
456 |
endif |
457 |
|
458 |
if (i_is_Dipole) then |
459 |
|
460 |
if (j_is_Charge) then |
461 |
|
462 |
ri2 = riji * riji |
463 |
ri3 = ri2 * riji |
464 |
|
465 |
pref = pre12 * q_j * mu_i |
466 |
vterm = pref * ct_i * riji * riji |
467 |
vpair = vpair + vterm |
468 |
epot = epot + sw * vterm |
469 |
|
470 |
dudx = dudx + pref * sw * ri3 * ( ul_i(1) - 3.0d0 * ct_i * drdxj) |
471 |
dudy = dudy + pref * sw * ri3 * ( ul_i(2) - 3.0d0 * ct_i * drdyj) |
472 |
dudz = dudz + pref * sw * ri3 * ( ul_i(3) - 3.0d0 * ct_i * drdzj) |
473 |
|
474 |
duduix = duduix + pref * sw * ri2 * drdxj |
475 |
duduiy = duduiy + pref * sw * ri2 * drdyj |
476 |
duduiz = duduiz + pref * sw * ri2 * drdzj |
477 |
endif |
478 |
|
479 |
if (j_is_Dipole) then |
480 |
|
481 |
ct_ij = ul_i(1)*ul_j(1) + ul_i(2)*ul_j(2) + ul_i(3)*ul_j(3) |
482 |
ri2 = riji * riji |
483 |
ri3 = ri2 * riji |
484 |
ri4 = ri2 * ri2 |
485 |
|
486 |
pref = pre22 * mu_i * mu_j |
487 |
vterm = pref * ri3 * (ct_ij - 3.0d0 * ct_i * ct_j) |
488 |
vpair = vpair + vterm |
489 |
epot = epot + sw * vterm |
490 |
|
491 |
a1 = 5.0d0 * ct_i * ct_j - ct_ij |
492 |
|
493 |
dudx = dudx + pref*sw*3.0d0*ri4*(a1*drdxj-ct_i*ul_j(1)-ct_j*ul_i(1)) |
494 |
dudy = dudy + pref*sw*3.0d0*ri4*(a1*drdyj-ct_i*ul_j(2)-ct_j*ul_i(2)) |
495 |
dudz = dudz + pref*sw*3.0d0*ri4*(a1*drdzj-ct_i*ul_j(3)-ct_j*ul_i(3)) |
496 |
|
497 |
duduix = duduix + pref*sw*ri3*(ul_j(1) - 3.0d0*ct_j*drdxj) |
498 |
duduiy = duduiy + pref*sw*ri3*(ul_j(2) - 3.0d0*ct_j*drdyj) |
499 |
duduiz = duduiz + pref*sw*ri3*(ul_j(3) - 3.0d0*ct_j*drdzj) |
500 |
|
501 |
dudujx = dudujx + pref*sw*ri3*(ul_i(1) - 3.0d0*ct_i*drdxj) |
502 |
dudujy = dudujy + pref*sw*ri3*(ul_i(2) - 3.0d0*ct_i*drdyj) |
503 |
dudujz = dudujz + pref*sw*ri3*(ul_i(3) - 3.0d0*ct_i*drdzj) |
504 |
endif |
505 |
|
506 |
endif |
507 |
|
508 |
if (do_pot) then |
509 |
#ifdef IS_MPI |
510 |
pot_row(atom1) = pot_row(atom1) + 0.5d0*epot |
511 |
pot_col(atom2) = pot_col(atom2) + 0.5d0*epot |
512 |
#else |
513 |
pot = pot + epot |
514 |
#endif |
515 |
endif |
516 |
|
517 |
#ifdef IS_MPI |
518 |
f_Row(1,atom1) = f_Row(1,atom1) + dudx |
519 |
f_Row(2,atom1) = f_Row(2,atom1) + dudy |
520 |
f_Row(3,atom1) = f_Row(3,atom1) + dudz |
521 |
|
522 |
f_Col(1,atom2) = f_Col(1,atom2) - dudx |
523 |
f_Col(2,atom2) = f_Col(2,atom2) - dudy |
524 |
f_Col(3,atom2) = f_Col(3,atom2) - dudz |
525 |
|
526 |
if (i_is_Dipole .or. i_is_Quadrupole) then |
527 |
t_Row(1,atom1) = t_Row(1,atom1) - ul_i(2)*duduiz + ul_i(3)*duduiy |
528 |
t_Row(2,atom1) = t_Row(2,atom1) - ul_i(3)*duduix + ul_i(1)*duduiz |
529 |
t_Row(3,atom1) = t_Row(3,atom1) - ul_i(1)*duduiy + ul_i(2)*duduix |
530 |
endif |
531 |
|
532 |
if (j_is_Dipole .or. j_is_Quadrupole) then |
533 |
t_Col(1,atom2) = t_Col(1,atom2) - ul_j(2)*dudujz + ul_j(3)*dudujy |
534 |
t_Col(2,atom2) = t_Col(2,atom2) - ul_j(3)*dudujx + ul_j(1)*dudujz |
535 |
t_Col(3,atom2) = t_Col(3,atom2) - ul_j(1)*dudujy + ul_j(2)*dudujx |
536 |
endif |
537 |
|
538 |
#else |
539 |
f(1,atom1) = f(1,atom1) + dudx |
540 |
f(2,atom1) = f(2,atom1) + dudy |
541 |
f(3,atom1) = f(3,atom1) + dudz |
542 |
|
543 |
f(1,atom2) = f(1,atom2) - dudx |
544 |
f(2,atom2) = f(2,atom2) - dudy |
545 |
f(3,atom2) = f(3,atom2) - dudz |
546 |
|
547 |
if (i_is_Dipole .or. i_is_Quadrupole) then |
548 |
t(1,atom1) = t(1,atom1) - ul_i(2)*duduiz + ul_i(3)*duduiy |
549 |
t(2,atom1) = t(2,atom1) - ul_i(3)*duduix + ul_i(1)*duduiz |
550 |
t(3,atom1) = t(3,atom1) - ul_i(1)*duduiy + ul_i(2)*duduix |
551 |
endif |
552 |
|
553 |
if (j_is_Dipole .or. j_is_Quadrupole) then |
554 |
t(1,atom2) = t(1,atom2) - ul_j(2)*dudujz + ul_j(3)*dudujy |
555 |
t(2,atom2) = t(2,atom2) - ul_j(3)*dudujx + ul_j(1)*dudujz |
556 |
t(3,atom2) = t(3,atom2) - ul_j(1)*dudujy + ul_j(2)*dudujx |
557 |
endif |
558 |
#endif |
559 |
|
560 |
#ifdef IS_MPI |
561 |
id1 = AtomRowToGlobal(atom1) |
562 |
id2 = AtomColToGlobal(atom2) |
563 |
#else |
564 |
id1 = atom1 |
565 |
id2 = atom2 |
566 |
#endif |
567 |
|
568 |
if (molMembershipList(id1) .ne. molMembershipList(id2)) then |
569 |
|
570 |
fpair(1) = fpair(1) + dudx |
571 |
fpair(2) = fpair(2) + dudy |
572 |
fpair(3) = fpair(3) + dudz |
573 |
|
574 |
endif |
575 |
|
576 |
return |
577 |
end subroutine doElectrostaticPair |
578 |
|
579 |
end module electrostatic_module |
580 |
|