1c
2c
3c     ############################################################
4c     ##  COPYRIGHT (C) 2018 by Joshua Rackers & Jay W. Ponder  ##
5c     ##                   All Rights Reserved                  ##
6c     ############################################################
7c
8c     ################################################################
9c     ##                                                            ##
10c     ##  subroutine echgtrn1  --  charge transfer energy & derivs  ##
11c     ##                                                            ##
12c     ################################################################
13c
14c
15c     "echgtrn1" calculates the charge transfer energy and first
16c     derivatives with respect to Cartesian coordinates
17c
18c
19      subroutine echgtrn1
20      use limits
21      implicit none
22c
23c
24c     choose method for summing over charge transfer interactions
25c
26      if (use_mlist) then
27         call echgtrn1b
28      else
29         call echgtrn1a
30      end if
31      return
32      end
33c
34c
35c     #################################################################
36c     ##                                                             ##
37c     ##  subroutine echgtrn1a  --  charge transfer derivs via loop  ##
38c     ##                                                             ##
39c     #################################################################
40c
41c
42c     "echgtrn1a" calculates the charge transfer interaction energy
43c     and first derivatives using a double loop
44c
45c
46      subroutine echgtrn1a
47      use atoms
48      use bound
49      use chgpot
50      use chgtrn
51      use cell
52      use couple
53      use ctrpot
54      use deriv
55      use energi
56      use group
57      use mplpot
58      use mpole
59      use shunt
60      use usage
61      use virial
62      implicit none
63      integer i,j,k
64      integer ii,kk
65      integer jcell
66      real*8 e,de,f,fgrp
67      real*8 rr1,r,r2
68      real*8 r3,r4,r5
69      real*8 xi,yi,zi
70      real*8 xr,yr,zr
71      real*8 chgi,chgk
72      real*8 chgik
73      real*8 alphai,alphak
74      real*8 alphaik
75      real*8 expi,expk
76      real*8 expik
77      real*8 frcx,frcy,frcz
78      real*8 vxx,vyy,vzz
79      real*8 vxy,vxz,vyz
80      real*8 taper,dtaper
81      real*8, allocatable :: mscale(:)
82      logical proceed,usei
83      character*6 mode
84c
85c
86c     zero out the charge transfer energy and first derivatives
87c
88      ect = 0.0d0
89      do i = 1, n
90         dect(1,i) = 0.0d0
91         dect(2,i) = 0.0d0
92         dect(3,i) = 0.0d0
93      end do
94      if (nct .eq. 0)  return
95c
96c     perform dynamic allocation of some local arrays
97c
98      allocate (mscale(n))
99c
100c     initialize connected atom exclusion coefficients
101c
102      do i = 1, n
103         mscale(i) = 1.0d0
104      end do
105c
106c     set conversion factor, cutoff and switching coefficients
107c
108      f = electric / dielec
109      mode = 'CHGTRN'
110      call switch (mode)
111c
112c     calculate the charge transfer energy and derivatives
113c
114      do ii = 1, npole-1
115         i = ipole(ii)
116         xi = x(i)
117         yi = y(i)
118         zi = z(i)
119         chgi = chgct(ii)
120         alphai = dmpct(ii)
121         if (alphai .eq. 0.0d0)  alphai = 1000.0d0
122         usei = use(i)
123c
124c     set exclusion coefficients for connected atoms
125c
126         do j = 1, n12(i)
127            mscale(i12(j,i)) = m2scale
128         end do
129         do j = 1, n13(i)
130            mscale(i13(j,i)) = m3scale
131         end do
132         do j = 1, n14(i)
133            mscale(i14(j,i)) = m4scale
134         end do
135         do j = 1, n15(i)
136            mscale(i15(j,i)) = m5scale
137         end do
138c
139c     evaluate all sites within the cutoff distance
140c
141         do kk = ii+1, npole
142            k = ipole(kk)
143            proceed = .true.
144            if (use_group)  call groups (proceed,fgrp,i,k,0,0,0,0)
145            if (.not. use_intra)  proceed = .true.
146            if (proceed)  proceed = (usei .or. use(k))
147            if (proceed) then
148               xr = x(k) - xi
149               yr = y(k) - yi
150               zr = z(k) - zi
151               if (use_bounds)  call image (xr,yr,zr)
152               r2 = xr*xr + yr* yr + zr*zr
153               if (r2 .le. off2) then
154                  r = sqrt(r2)
155                  rr1 = 1.0d0 / r
156                  chgk = chgct(kk)
157                  alphak = dmpct(kk)
158                  if (alphak .eq. 0.0d0)  alphak = 1000.0d0
159                  if (ctrntyp .eq. 'SEPARATE') then
160                     expi = exp(-alphai*r)
161                     expk = exp(-alphak*r)
162                     e = -chgi*expk - chgk*expi
163                     de = chgi*expk*alphak + chgk*expi*alphai
164                  else
165                     chgik = sqrt(abs(chgi*chgk))
166                     alphaik = 0.5d0 * (alphai+alphak)
167                     expik = exp(-alphaik*r)
168                     e = -chgik * expik
169                     de = -e * alphaik
170                  end if
171                  e = f * e * mscale(k)
172                  de = f * de * mscale(k)
173c
174c     use energy switching if near the cutoff distance
175c
176                  if (r2 .gt. cut2) then
177                     r3 = r2 * r
178                     r4 = r2 * r2
179                     r5 = r2 * r3
180                     taper = c5*r5 + c4*r4 + c3*r3
181     &                          + c2*r2 + c1*r + c0
182                     dtaper = 5.0d0*c5*r4 + 4.0d0*c4*r3
183     &                           + 3.0d0*c3*r2 + 2.0d0*c2*r + c1
184                     de = e*dtaper + de*taper
185                     e = e * taper
186                  end if
187c
188c     scale the interaction based on its group membership
189c
190                  if (use_group) then
191                     e = e * fgrp
192                     de = de * fgrp
193                  end if
194c
195c     compute the force components for this interaction
196c
197                  frcx = de * xr * rr1
198                  frcy = de * yr * rr1
199                  frcz = de * zr * rr1
200c
201c     increment the total charge transfer energy and derivatives
202c
203                  ect = ect + e
204                  dect(1,i) = dect(1,i) - frcx
205                  dect(2,i) = dect(2,i) - frcy
206                  dect(3,i) = dect(3,i) - frcz
207                  dect(1,k) = dect(1,k) + frcx
208                  dect(2,k) = dect(2,k) + frcy
209                  dect(3,k) = dect(3,k) + frcz
210c
211c     increment the internal virial tensor components
212c
213                  vxx = xr * frcx
214                  vxy = yr * frcx
215                  vxz = zr * frcx
216                  vyy = yr * frcy
217                  vyz = zr * frcy
218                  vzz = zr * frcz
219                  vir(1,1) = vir(1,1) + vxx
220                  vir(2,1) = vir(2,1) + vxy
221                  vir(3,1) = vir(3,1) + vxz
222                  vir(1,2) = vir(1,2) + vxy
223                  vir(2,2) = vir(2,2) + vyy
224                  vir(3,2) = vir(3,2) + vyz
225                  vir(1,3) = vir(1,3) + vxz
226                  vir(2,3) = vir(2,3) + vyz
227                  vir(3,3) = vir(3,3) + vzz
228               end if
229            end if
230         end do
231c
232c     reset exclusion coefficients for connected atoms
233c
234         do j = 1, n12(i)
235            mscale(i12(j,i)) = 1.0d0
236         end do
237         do j = 1, n13(i)
238            mscale(i13(j,i)) = 1.0d0
239         end do
240         do j = 1, n14(i)
241            mscale(i14(j,i)) = 1.0d0
242         end do
243         do j = 1, n15(i)
244            mscale(i15(j,i)) = 1.0d0
245         end do
246      end do
247c
248c     for periodic boundary conditions with large cutoffs
249c     neighbors must be found by the replicates method
250c
251      if (use_replica) then
252c
253c     calculate interaction components with other unit cells
254c
255         do ii = 1, npole
256            i = ipole(ii)
257            xi = x(i)
258            yi = y(i)
259            zi = z(i)
260            chgi = chgct(ii)
261            alphai = dmpct(ii)
262            if (alphai .eq. 0.0d0)  alphai = 1000.0d0
263            usei = use(i)
264c
265c     set exclusion coefficients for connected atoms
266c
267            do j = 1, n12(i)
268               mscale(i12(j,i)) = m2scale
269            end do
270            do j = 1, n13(i)
271               mscale(i13(j,i)) = m3scale
272            end do
273            do j = 1, n14(i)
274               mscale(i14(j,i)) = m4scale
275            end do
276            do j = 1, n15(i)
277               mscale(i15(j,i)) = m5scale
278            end do
279c
280c     evaluate all sites within the cutoff distance
281c
282            do kk = ii, npole
283               k = ipole(kk)
284               proceed = .true.
285               if (use_group)  call groups (proceed,fgrp,i,k,0,0,0,0)
286               if (.not. use_intra)  proceed = .true.
287               if (proceed)  proceed = (usei .or. use(k))
288               if (proceed) then
289                  do jcell = 2, ncell
290                     xr = x(k) - xi
291                     yr = y(k) - yi
292                     zr = z(k) - zi
293                     call imager (xr,yr,zr,jcell)
294                     r2 = xr*xr + yr* yr + zr*zr
295                     if (r2 .le. off2) then
296                        r = sqrt(r2)
297                        rr1 = 1.0d0 / r
298                        chgk = chgct(kk)
299                        alphak = dmpct(kk)
300                        if (alphak .eq. 0.0d0)  alphak = 1000.0d0
301                        if (ctrntyp .eq. 'SEPARATE') then
302                           expi = exp(-alphai*r)
303                           expk = exp(-alphak*r)
304                           e = -chgi*expk - chgk*expi
305                           de = chgi*expk*alphak + chgk*expi*alphai
306                        else
307                           chgik = sqrt(abs(chgi*chgk))
308                           alphaik = 0.5d0 * (alphai+alphak)
309                           expik = exp(-alphaik*r)
310                           e = -chgik * expik
311                           de = -e * alphaik
312                        end if
313                        if (use_group) then
314                           e = e * fgrp
315                           de = de * fgrp
316                        end if
317                        e = f * e * mscale(k)
318                        de = f * de * mscale(k)
319                        if (i .eq. k) then
320                           e = 0.5d0 * e
321                           de = 0.5d0 * de
322                        end if
323c
324c     use energy switching if near the cutoff distance
325c
326                        if (r2 .gt. cut2) then
327                           r3 = r2 * r
328                           r4 = r2 * r2
329                           r5 = r2 * r3
330                           taper = c5*r5 + c4*r4 + c3*r3
331     &                                + c2*r2 + c1*r + c0
332                           dtaper = 5.0d0*c5*r4 + 4.0d0*c4*r3
333     &                                 + 3.0d0*c3*r2 + 2.0d0*c2*r + c1
334                           de = e*dtaper + de*taper
335                           e = e * taper
336                        end if
337c
338c     scale the interaction based on its group membership
339c
340                        if (use_group) then
341                           e = e * fgrp
342                           de = de * fgrp
343                        end if
344c
345c     compute the force components for this interaction
346c
347                        frcx = de * xr * rr1
348                        frcy = de * yr * rr1
349                        frcz = de * zr * rr1
350c
351c     increment the total charge transfer energy and derivatives
352c
353                        ect = ect + e
354                        dect(1,i) = dect(1,i) - frcx
355                        dect(2,i) = dect(2,i) - frcy
356                        dect(3,i) = dect(3,i) - frcz
357                        dect(1,k) = dect(1,k) + frcx
358                        dect(2,k) = dect(2,k) + frcy
359                        dect(3,k) = dect(3,k) + frcz
360c
361c     increment the internal virial tensor components
362c
363                        vxx = xr * frcx
364                        vxy = yr * frcx
365                        vxz = zr * frcx
366                        vyy = yr * frcy
367                        vyz = zr * frcy
368                        vzz = zr * frcz
369                        vir(1,1) = vir(1,1) + vxx
370                        vir(2,1) = vir(2,1) + vxy
371                        vir(3,1) = vir(3,1) + vxz
372                        vir(1,2) = vir(1,2) + vxy
373                        vir(2,2) = vir(2,2) + vyy
374                        vir(3,2) = vir(3,2) + vyz
375                        vir(1,3) = vir(1,3) + vxz
376                        vir(2,3) = vir(2,3) + vyz
377                        vir(3,3) = vir(3,3) + vzz
378                     end if
379                  end do
380               end if
381            end do
382c
383c     reset exclusion coefficients for connected atoms
384c
385            do j = 1, n12(i)
386               mscale(i12(j,i)) = 1.0d0
387            end do
388            do j = 1, n13(i)
389               mscale(i13(j,i)) = 1.0d0
390            end do
391            do j = 1, n14(i)
392               mscale(i14(j,i)) = 1.0d0
393            end do
394            do j = 1, n15(i)
395               mscale(i15(j,i)) = 1.0d0
396            end do
397         end do
398      end if
399c
400c     perform deallocation of some local arrays
401c
402      deallocate (mscale)
403      return
404      end
405c
406c
407c     #################################################################
408c     ##                                                             ##
409c     ##  subroutine echgtrn1b  --  charge transfer derivs via list  ##
410c     ##                                                             ##
411c     #################################################################
412c
413c
414c     "echgtrn1b" calculates the charge transfer energy and first
415c     derivatives using a pairwise neighbor list
416c
417c
418      subroutine echgtrn1b
419      use atoms
420      use bound
421      use chgpot
422      use chgtrn
423      use cell
424      use couple
425      use ctrpot
426      use deriv
427      use energi
428      use group
429      use mplpot
430      use mpole
431      use neigh
432      use shunt
433      use usage
434      use virial
435      implicit none
436      integer i,j,k
437      integer ii,kk,kkk
438      real*8 e,de,f,fgrp
439      real*8 rr1,r,r2
440      real*8 r3,r4,r5
441      real*8 xi,yi,zi
442      real*8 xr,yr,zr
443      real*8 chgi,chgk
444      real*8 chgik
445      real*8 alphai,alphak
446      real*8 alphaik
447      real*8 expi,expk
448      real*8 expik
449      real*8 frcx,frcy,frcz
450      real*8 vxx,vyy,vzz
451      real*8 vxy,vxz,vyz
452      real*8 taper,dtaper
453      real*8, allocatable :: mscale(:)
454      logical proceed,usei
455      character*6 mode
456c
457c
458c     zero out the charge transfer energy and first derivatives
459c
460      ect = 0.0d0
461      do i = 1, n
462         dect(1,i) = 0.0d0
463         dect(2,i) = 0.0d0
464         dect(3,i) = 0.0d0
465      end do
466      if (nct .eq. 0)  return
467c
468c     perform dynamic allocation of some local arrays
469c
470      allocate (mscale(n))
471c
472c     initialize connected atom exclusion coefficients
473c
474      do i = 1, n
475         mscale(i) = 1.0d0
476      end do
477c
478c     set conversion factor, cutoff and switching coefficients
479c
480      f = electric / dielec
481      mode = 'CHGTRN'
482      call switch (mode)
483c
484c     OpenMP directives for the major loop structure
485c
486!$OMP PARALLEL default(private)
487!$OMP& shared(npole,ipole,x,y,z,chgct,dmpct,n12,i12,n13,i13,
488!$OMP& n14,i14,n15,i15,m2scale,m3scale,m4scale,m5scale,nelst,
489!$OMP& elst,use,use_group,use_intra,use_bounds,ctrntyp,f,off2,
490!$OMP& cut2,c0,c1,c2,c3,c4,c5)
491!$OMP& firstprivate(mscale) shared(ect,dect,vir)
492!$OMP DO reduction(+:ect,dect,vir) schedule(guided)
493c
494c     compute the charge transfer energy and derivatives
495c
496      do ii = 1, npole
497         i = ipole(ii)
498         xi = x(i)
499         yi = y(i)
500         zi = z(i)
501         chgi = chgct(ii)
502         alphai = dmpct(ii)
503         if (alphai .eq. 0.0d0)  alphai = 1000.0d0
504         usei = use(i)
505c
506c     set exclusion coefficients for connected atoms
507c
508         do j = 1, n12(i)
509            mscale(i12(j,i)) = m2scale
510         end do
511         do j = 1, n13(i)
512            mscale(i13(j,i)) = m3scale
513         end do
514         do j = 1, n14(i)
515            mscale(i14(j,i)) = m4scale
516         end do
517         do j = 1, n15(i)
518            mscale(i15(j,i)) = m5scale
519         end do
520c
521c     evaluate all sites within the cutoff distance
522c
523         do kkk = 1, nelst(ii)
524            kk = elst(kkk,ii)
525            k = ipole(kk)
526            proceed = .true.
527            if (use_group)  call groups (proceed,fgrp,i,k,0,0,0,0)
528            if (.not. use_intra)  proceed = .true.
529            if (proceed)  proceed = (usei .or. use(k))
530            if (proceed) then
531               xr = x(k) - xi
532               yr = y(k) - yi
533               zr = z(k) - zi
534               if (use_bounds)  call image (xr,yr,zr)
535               r2 = xr*xr + yr* yr + zr*zr
536               if (r2 .le. off2) then
537                  r = sqrt(r2)
538                  rr1 = 1.0d0 / r
539                  chgk = chgct(kk)
540                  alphak = dmpct(kk)
541                  if (alphak .eq. 0.0d0)  alphak = 1000.0d0
542                  if (ctrntyp .eq. 'SEPARATE') then
543                     expi = exp(-alphai*r)
544                     expk = exp(-alphak*r)
545                     e = -chgi*expk - chgk*expi
546                     de = chgi*expk*alphak + chgk*expi*alphai
547                  else
548                     chgik = sqrt(abs(chgi*chgk))
549                     alphaik = 0.5d0 * (alphai+alphak)
550                     expik = exp(-alphaik*r)
551                     e = -chgik * expik
552                     de = -e * alphaik
553                  end if
554                  e = f * e * mscale(k)
555                  de = f * de * mscale(k)
556c
557c     use energy switching if near the cutoff distance
558c
559                  if (r2 .gt. cut2) then
560                     r3 = r2 * r
561                     r4 = r2 * r2
562                     r5 = r2 * r3
563                     taper = c5*r5 + c4*r4 + c3*r3
564     &                          + c2*r2 + c1*r + c0
565                     dtaper = 5.0d0*c5*r4 + 4.0d0*c4*r3
566     &                           + 3.0d0*c3*r2 + 2.0d0*c2*r + c1
567                     de = e*dtaper + de*taper
568                     e = e * taper
569                  end if
570c
571c     scale the interaction based on its group membership
572c
573                  if (use_group) then
574                     e = e * fgrp
575                     de = de * fgrp
576                  end if
577c
578c     compute the force components for this interaction
579c
580                  frcx = de * xr * rr1
581                  frcy = de * yr * rr1
582                  frcz = de * zr * rr1
583c
584c     increment the total charge transfer energy and derivatives
585c
586                  ect = ect + e
587                  dect(1,i) = dect(1,i) - frcx
588                  dect(2,i) = dect(2,i) - frcy
589                  dect(3,i) = dect(3,i) - frcz
590                  dect(1,k) = dect(1,k) + frcx
591                  dect(2,k) = dect(2,k) + frcy
592                  dect(3,k) = dect(3,k) + frcz
593c
594c     increment the internal virial tensor components
595c
596                  vxx = xr * frcx
597                  vxy = yr * frcx
598                  vxz = zr * frcx
599                  vyy = yr * frcy
600                  vyz = zr * frcy
601                  vzz = zr * frcz
602                  vir(1,1) = vir(1,1) + vxx
603                  vir(2,1) = vir(2,1) + vxy
604                  vir(3,1) = vir(3,1) + vxz
605                  vir(1,2) = vir(1,2) + vxy
606                  vir(2,2) = vir(2,2) + vyy
607                  vir(3,2) = vir(3,2) + vyz
608                  vir(1,3) = vir(1,3) + vxz
609                  vir(2,3) = vir(2,3) + vyz
610                  vir(3,3) = vir(3,3) + vzz
611               end if
612            end if
613         end do
614c
615c     reset exclusion coefficients for connected atoms
616c
617         do j = 1, n12(i)
618            mscale(i12(j,i)) = 1.0d0
619         end do
620         do j = 1, n13(i)
621            mscale(i13(j,i)) = 1.0d0
622         end do
623         do j = 1, n14(i)
624            mscale(i14(j,i)) = 1.0d0
625         end do
626         do j = 1, n15(i)
627            mscale(i15(j,i)) = 1.0d0
628         end do
629      end do
630c
631c     OpenMP directives for the major loop structure
632c
633!$OMP END DO
634!$OMP END PARALLEL
635c
636c     perform deallocation of some local arrays
637c
638      deallocate (mscale)
639      return
640      end
641