1c
2c
3c     ###################################################
4c     ##  COPYRIGHT (C)  1990  by  Jay William Ponder  ##
5c     ##              All Rights Reserved              ##
6c     ###################################################
7c
8c     ##############################################################
9c     ##                                                          ##
10c     ##  subroutine edipole1  --  dipole-dipole energy & derivs  ##
11c     ##                                                          ##
12c     ##############################################################
13c
14c
15c     "edipole1" calculates the dipole-dipole interaction energy
16c     and first derivatives with respect to Cartesian coordinates
17c
18c
19      subroutine edipole1
20      use atoms
21      use bound
22      use cell
23      use chgpot
24      use deriv
25      use dipole
26      use energi
27      use group
28      use shunt
29      use units
30      use usage
31      use virial
32      implicit none
33      integer i,j,k
34      integer i1,i2,k1,k2
35      real*8 xi,yi,zi
36      real*8 xk,yk,zk
37      real*8 xq,yq,zq
38      real*8 xr,yr,zr
39      real*8 xq1,yq1,zq1
40      real*8 xq2,yq2,zq2
41      real*8 f,fi,fik,fgrp
42      real*8 e,r2,ri2,rk2,rirkr3
43      real*8 doti,dotk,dotp
44      real*8 si1,si2,sk1,sk2
45      real*8 de,dedr,dedrirk
46      real*8 deddoti,deddotk,deddotp
47      real*8 termx,termy,termz
48      real*8 dedrirkri2,dedrirkrk2
49      real*8 termxi,termyi,termzi
50      real*8 termxk,termyk,termzk
51      real*8 dedxi1,dedyi1,dedzi1
52      real*8 dedxi2,dedyi2,dedzi2
53      real*8 dedxk1,dedyk1,dedzk1
54      real*8 dedxk2,dedyk2,dedzk2
55      real*8 r,r3,r4,r5,taper,dtaper
56      real*8 dtaperx,dtapery,dtaperz
57      real*8 vxx,vyy,vzz
58      real*8 vyx,vzx,vzy
59      logical proceed
60      character*6 mode
61c
62c
63c     zero out the overall dipole interaction energy and derivs,
64c     then set up the constants for the calculation
65c
66      ed = 0.0d0
67      do i = 1, n
68         ded(1,i) = 0.0d0
69         ded(2,i) = 0.0d0
70         ded(3,i) = 0.0d0
71      end do
72      if (ndipole .eq. 0)  return
73c
74c     set conversion factor and switching function coefficients
75c
76      f = electric / (debye**2 * dielec)
77      mode = 'DIPOLE'
78      call switch (mode)
79c
80c     compute the dipole interaction energy and first derivatives
81c
82      do i = 1, ndipole-1
83         i1 = idpl(1,i)
84         i2 = idpl(2,i)
85         si1 = 1.0d0 - sdpl(i)
86         si2 = sdpl(i)
87         xi = x(i2) - x(i1)
88         yi = y(i2) - y(i1)
89         zi = z(i2) - z(i1)
90         if (use_polymer)  call imager (xi,yi,zi,-1)
91         ri2 = xi*xi + yi*yi + zi*zi
92         xq = x(i1) + xi*si2
93         yq = y(i1) + yi*si2
94         zq = z(i1) + zi*si2
95         fi = f * bdpl(i)
96c
97c     decide whether to compute the current interaction
98c
99         do k = i+1, ndipole
100            k1 = idpl(1,k)
101            k2 = idpl(2,k)
102            proceed = .true.
103            if (use_group)  call groups (proceed,fgrp,i1,i2,k1,k2,0,0)
104            if (proceed)  proceed = (use(i1) .or. use(i2) .or.
105     &                                 use(k1) .or. use(k2))
106            if (proceed)  proceed = (k1.ne.i1 .and. k1.ne.i2 .and.
107     &                                 k2.ne.i1 .and. k2.ne.i2)
108c
109c     compute the energy contribution for this interaction
110c
111            if (proceed) then
112               sk1 = 1.0d0 - sdpl(k)
113               sk2 = sdpl(k)
114               xk = x(k2) - x(k1)
115               yk = y(k2) - y(k1)
116               zk = z(k2) - z(k1)
117               if (use_polymer)  call imager (xk,yk,zk,-1)
118               xr = xq - x(k1) - xk*sk2
119               yr = yq - y(k1) - yk*sk2
120               zr = zq - z(k1) - zk*sk2
121               call image (xr,yr,zr)
122               r2 = xr*xr + yr*yr + zr*zr
123               if (r2 .le. off2) then
124                  rk2 = xk*xk + yk*yk + zk*zk
125                  rirkr3 = sqrt(ri2*rk2*r2) * r2
126                  dotp = xi*xk + yi*yk + zi*zk
127                  doti = xi*xr + yi*yr + zi*zr
128                  dotk = xk*xr + yk*yr + zk*zr
129                  fik = fi * bdpl(k)
130c
131c     form the energy and master chain rule term for derivatives
132c
133                  e = fik * (dotp-3.0d0*doti*dotk/r2) / rirkr3
134                  de = -fik / (rirkr3*r2)
135c
136c     scale the interaction based on its group membership
137c
138                  if (use_group) then
139                     e = e * fgrp
140                     de = de * fgrp
141                  end if
142c
143c     secondary chain rule terms for derivative expressions
144c
145                  deddotp = -de * r2
146                  deddoti = de * 3.0d0*dotk
147                  deddotk = de * 3.0d0*doti
148                  dedr = de * (3.0d0*dotp-15.0d0*doti*dotk/r2)
149                  dedrirk = -e
150                  dedrirkri2 = dedrirk / ri2
151                  dedrirkrk2 = dedrirk / rk2
152c
153c     more chain rule terms for derivative expressions
154c
155                  termx = dedr*xr + deddoti*xi + deddotk*xk
156                  termy = dedr*yr + deddoti*yi + deddotk*yk
157                  termz = dedr*zr + deddoti*zi + deddotk*zk
158                  termxi = dedrirkri2*xi + deddotp*xk + deddoti*xr
159                  termyi = dedrirkri2*yi + deddotp*yk + deddoti*yr
160                  termzi = dedrirkri2*zi + deddotp*zk + deddoti*zr
161                  termxk = dedrirkrk2*xk + deddotp*xi + deddotk*xr
162                  termyk = dedrirkrk2*yk + deddotp*yi + deddotk*yr
163                  termzk = dedrirkrk2*zk + deddotp*zi + deddotk*zr
164c
165c     finally, the individual first derivative components
166c
167                  dedxi1 = si1*termx - termxi
168                  dedyi1 = si1*termy - termyi
169                  dedzi1 = si1*termz - termzi
170                  dedxi2 = si2*termx + termxi
171                  dedyi2 = si2*termy + termyi
172                  dedzi2 = si2*termz + termzi
173                  dedxk1 = -sk1*termx - termxk
174                  dedyk1 = -sk1*termy - termyk
175                  dedzk1 = -sk1*termz - termzk
176                  dedxk2 = -sk2*termx + termxk
177                  dedyk2 = -sk2*termy + termyk
178                  dedzk2 = -sk2*termz + termzk
179c
180c     use energy switching if near the cutoff distance
181c
182                  if (r2 .gt. cut2) then
183                     r = sqrt(r2)
184                     r3 = r2 * r
185                     r4 = r2 * r2
186                     r5 = r2 * r3
187                     taper = c5*r5 + c4*r4 + c3*r3
188     &                          + c2*r2 + c1*r + c0
189                     dtaper = 5.0d0*c5*r4 + 4.0d0*c4*r3
190     &                           + 3.0d0*c3*r2 + 2.0d0*c2*r + c1
191                     dtaper = dtaper * e/r
192                     dtaperx = xr * dtaper
193                     dtapery = yr * dtaper
194                     dtaperz = zr * dtaper
195                     e = e * taper
196                     dedxi1 = dedxi1*taper + si1*dtaperx
197                     dedyi1 = dedyi1*taper + si1*dtapery
198                     dedzi1 = dedzi1*taper + si1*dtaperz
199                     dedxi2 = dedxi2*taper + si2*dtaperx
200                     dedyi2 = dedyi2*taper + si2*dtapery
201                     dedzi2 = dedzi2*taper + si2*dtaperz
202                     dedxk1 = dedxk1*taper - sk1*dtaperx
203                     dedyk1 = dedyk1*taper - sk1*dtapery
204                     dedzk1 = dedzk1*taper - sk1*dtaperz
205                     dedxk2 = dedxk2*taper - sk2*dtaperx
206                     dedyk2 = dedyk2*taper - sk2*dtapery
207                     dedzk2 = dedzk2*taper - sk2*dtaperz
208                  end if
209c
210c     increment the overall energy and derivative expressions
211c
212                  ed = ed + e
213                  ded(1,i1) = ded(1,i1) + dedxi1
214                  ded(2,i1) = ded(2,i1) + dedyi1
215                  ded(3,i1) = ded(3,i1) + dedzi1
216                  ded(1,i2) = ded(1,i2) + dedxi2
217                  ded(2,i2) = ded(2,i2) + dedyi2
218                  ded(3,i2) = ded(3,i2) + dedzi2
219                  ded(1,k1) = ded(1,k1) + dedxk1
220                  ded(2,k1) = ded(2,k1) + dedyk1
221                  ded(3,k1) = ded(3,k1) + dedzk1
222                  ded(1,k2) = ded(1,k2) + dedxk2
223                  ded(2,k2) = ded(2,k2) + dedyk2
224                  ded(3,k2) = ded(3,k2) + dedzk2
225c
226c     increment the internal virial tensor components
227c
228                  xq1 = x(k1) - xq
229                  yq1 = y(k1) - yq
230                  zq1 = z(k1) - zq
231                  xq2 = x(k2) - xq
232                  yq2 = y(k2) - yq
233                  zq2 = z(k2) - zq
234                  vxx = xq1*dedxk1 + xq2*dedxk2
235                  vyx = 0.5d0 * (yq1*dedxk1 + yq2*dedxk2
236     &                              + xq1*dedyk1 + xq2*dedyk2)
237                  vzx = 0.5d0 * (zq1*dedxk1 + zq2*dedxk2
238     &                              + xq1*dedzk1 + xq2*dedzk2)
239                  vyy = yq1*dedyk1 + yq2*dedyk2
240                  vzy = 0.5d0 * (zq1*dedyk1 + zq2*dedyk2
241     &                              + yq1*dedzk1 + yq2*dedzk2)
242                  vzz = zq1*dedzk1 + zq2*dedzk2
243                  vir(1,1) = vir(1,1) + vxx
244                  vir(2,1) = vir(2,1) + vyx
245                  vir(3,1) = vir(3,1) + vzx
246                  vir(1,2) = vir(1,2) + vyx
247                  vir(2,2) = vir(2,2) + vyy
248                  vir(3,2) = vir(3,2) + vzy
249                  vir(1,3) = vir(1,3) + vzx
250                  vir(2,3) = vir(2,3) + vzy
251                  vir(3,3) = vir(3,3) + vzz
252               end if
253            end if
254         end do
255      end do
256c
257c     for periodic boundary conditions with large cutoffs
258c     neighbors must be found by the replicates method
259c
260      if (.not. use_replica)  return
261c
262c     calculate interaction energy with other unit cells
263c
264      do i = 1, ndipole
265         i1 = idpl(1,i)
266         i2 = idpl(2,i)
267         si1 = 1.0d0 - sdpl(i)
268         si2 = sdpl(i)
269         xi = x(i2) - x(i1)
270         yi = y(i2) - y(i1)
271         zi = z(i2) - z(i1)
272         if (use_polymer)  call imager (xi,yi,zi,-1)
273         ri2 = xi*xi + yi*yi + zi*zi
274         xq = x(i1) + xi*si2
275         yq = y(i1) + yi*si2
276         zq = z(i1) + zi*si2
277         fi = f * bdpl(i)
278c
279c     decide whether to compute the current interaction
280c
281         do k = i, ndipole
282            k1 = idpl(1,k)
283            k2 = idpl(2,k)
284            proceed = .true.
285            if (use_group)  call groups (proceed,fgrp,i1,i2,k1,k2,0,0)
286            if (proceed)  proceed = (use(i1) .or. use(i2) .or.
287     &                                 use(k1) .or. use(k2))
288c
289c     compute the energy contribution for this interaction
290c
291            if (proceed) then
292               sk1 = 1.0d0 - sdpl(k)
293               sk2 = sdpl(k)
294               do j = 2, ncell
295                  xk = x(k2) - x(k1)
296                  yk = y(k2) - y(k1)
297                  zk = z(k2) - z(k1)
298                  if (use_polymer)  call imager (xk,yk,zk,-1)
299                  xr = xq - x(k1) - xk*sk2
300                  yr = yq - y(k1) - yk*sk2
301                  zr = zq - z(k1) - zk*sk2
302                  call imager (xr,yr,zr,j)
303                  r2 = xr*xr + yr*yr + zr*zr
304                  if (r2 .le. off2) then
305                     rk2 = xk*xk + yk*yk + zk*zk
306                     rirkr3 = sqrt(ri2*rk2*r2) * r2
307                     dotp = xi*xk + yi*yk + zi*zk
308                     doti = xi*xr + yi*yr + zi*zr
309                     dotk = xk*xr + yk*yr + zk*zr
310                     fik = fi * bdpl(k)
311                     if (use_polymer) then
312                        if (r2 .lt. polycut2) then
313                           if (k1.eq.i1 .or. k1.eq.i2 .or.
314     &                         k2.eq.i1 .or. k2.eq.i2)  fik = 0.0d0
315                        end if
316                     end if
317c
318c     form the energy and master chain rule term for derivatives
319c
320                     e = fik * (dotp-3.0d0*doti*dotk/r2) / rirkr3
321                     de = -fik / (rirkr3*r2)
322c
323c     scale the interaction based on its group membership
324c
325                     if (use_group) then
326                        e = e * fgrp
327                        de = de * fgrp
328                     end if
329c
330c     secondary chain rule terms for derivative expressions
331c
332                     deddotp = -de * r2
333                     deddoti = de * 3.0d0*dotk
334                     deddotk = de * 3.0d0*doti
335                     dedr = de * (3.0d0*dotp-15.0d0*doti*dotk/r2)
336                     dedrirk = -e
337                     dedrirkri2 = dedrirk / ri2
338                     dedrirkrk2 = dedrirk / rk2
339c
340c     more chain rule terms for derivative expressions
341c
342                     termx = dedr*xr + deddoti*xi + deddotk*xk
343                     termy = dedr*yr + deddoti*yi + deddotk*yk
344                     termz = dedr*zr + deddoti*zi + deddotk*zk
345                     termxi = dedrirkri2*xi + deddotp*xk + deddoti*xr
346                     termyi = dedrirkri2*yi + deddotp*yk + deddoti*yr
347                     termzi = dedrirkri2*zi + deddotp*zk + deddoti*zr
348                     termxk = dedrirkrk2*xk + deddotp*xi + deddotk*xr
349                     termyk = dedrirkrk2*yk + deddotp*yi + deddotk*yr
350                     termzk = dedrirkrk2*zk + deddotp*zi + deddotk*zr
351c
352c     finally, the individual first derivative components
353c
354                     dedxi1 = si1*termx - termxi
355                     dedyi1 = si1*termy - termyi
356                     dedzi1 = si1*termz - termzi
357                     dedxi2 = si2*termx + termxi
358                     dedyi2 = si2*termy + termyi
359                     dedzi2 = si2*termz + termzi
360                     dedxk1 = -sk1*termx - termxk
361                     dedyk1 = -sk1*termy - termyk
362                     dedzk1 = -sk1*termz - termzk
363                     dedxk2 = -sk2*termx + termxk
364                     dedyk2 = -sk2*termy + termyk
365                     dedzk2 = -sk2*termz + termzk
366c
367c     use energy switching if near the cutoff distance
368c
369                     if (r2 .gt. cut2) then
370                        r = sqrt(r2)
371                        r3 = r2 * r
372                        r4 = r2 * r2
373                        r5 = r2 * r3
374                        taper = c5*r5 + c4*r4 + c3*r3
375     &                             + c2*r2 + c1*r + c0
376                        dtaper = 5.0d0*c5*r4 + 4.0d0*c4*r3
377     &                              + 3.0d0*c3*r2 + 2.0d0*c2*r + c1
378                        dtaper = dtaper * e/r
379                        dtaperx = xr * dtaper
380                        dtapery = yr * dtaper
381                        dtaperz = zr * dtaper
382                        e = e * taper
383                        dedxi1 = dedxi1*taper + si1*dtaperx
384                        dedyi1 = dedyi1*taper + si1*dtapery
385                        dedzi1 = dedzi1*taper + si1*dtaperz
386                        dedxi2 = dedxi2*taper + si2*dtaperx
387                        dedyi2 = dedyi2*taper + si2*dtapery
388                        dedzi2 = dedzi2*taper + si2*dtaperz
389                        dedxk1 = dedxk1*taper - sk1*dtaperx
390                        dedyk1 = dedyk1*taper - sk1*dtapery
391                        dedzk1 = dedzk1*taper - sk1*dtaperz
392                        dedxk2 = dedxk2*taper - sk2*dtaperx
393                        dedyk2 = dedyk2*taper - sk2*dtapery
394                        dedzk2 = dedzk2*taper - sk2*dtaperz
395                     end if
396c
397c     increment the overall energy and derivative expressions
398c
399                     if (i .eq. k)  e = 0.5d0 * e
400                     ed = ed + e
401                     ded(1,i1) = ded(1,i1) + dedxi1
402                     ded(2,i1) = ded(2,i1) + dedyi1
403                     ded(3,i1) = ded(3,i1) + dedzi1
404                     ded(1,i2) = ded(1,i2) + dedxi2
405                     ded(2,i2) = ded(2,i2) + dedyi2
406                     ded(3,i2) = ded(3,i2) + dedzi2
407                     if (i .ne. k) then
408                        ded(1,k1) = ded(1,k1) + dedxk1
409                        ded(2,k1) = ded(2,k1) + dedyk1
410                        ded(3,k1) = ded(3,k1) + dedzk1
411                        ded(1,k2) = ded(1,k2) + dedxk2
412                        ded(2,k2) = ded(2,k2) + dedyk2
413                        ded(3,k2) = ded(3,k2) + dedzk2
414                     end if
415c
416c     increment the internal virial tensor components
417c
418                     xq1 = x(k1) - xq
419                     yq1 = y(k1) - yq
420                     zq1 = z(k1) - zq
421                     xq2 = x(k2) - xq
422                     yq2 = y(k2) - yq
423                     zq2 = z(k2) - zq
424                     vxx = xq1*dedxk1 + xq2*dedxk2
425                     vyx = 0.5d0 * (yq1*dedxk1 + yq2*dedxk2
426     &                                 + xq1*dedyk1 + xq2*dedyk2)
427                     vzx = 0.5d0 * (zq1*dedxk1 + zq2*dedxk2
428     &                                 + xq1*dedzk1 + xq2*dedzk2)
429                     vyy = yq1*dedyk1 + yq2*dedyk2
430                     vzy = 0.5d0 * (zq1*dedyk1 + zq2*dedyk2
431     &                                 + yq1*dedzk1 + yq2*dedzk2)
432                     vzz = zq1*dedzk1 + zq2*dedzk2
433                     vir(1,1) = vir(1,1) + vxx
434                     vir(2,1) = vir(2,1) + vyx
435                     vir(3,1) = vir(3,1) + vzx
436                     vir(1,2) = vir(1,2) + vyx
437                     vir(2,2) = vir(2,2) + vyy
438                     vir(3,2) = vir(3,2) + vzy
439                     vir(1,3) = vir(1,3) + vzx
440                     vir(2,3) = vir(2,3) + vzy
441                     vir(3,3) = vir(3,3) + vzz
442                  end if
443               end do
444            end if
445         end do
446      end do
447      return
448      end
449