1c
2c
3c     ##########################################################
4c     ##  COPYRIGHT (C) 2020 by Chengwen Liu & Jay W. Ponder  ##
5c     ##                 All Rights Reserved                  ##
6c     ##########################################################
7c
8c     ################################################################
9c     ##                                                            ##
10c     ##  subroutine kchgflx  --  charge flux parameter assignment  ##
11c     ##                                                            ##
12c     ################################################################
13c
14c
15c     "kchgflx" assigns bond stretch and angle bend charge flux
16c     correction values and processes any new or changed values
17c     for these parameters
18c
19c
20      subroutine kchgflx
21      use sizes
22      use angbnd
23      use atmlst
24      use atomid
25      use atoms
26      use bndstr
27      use cflux
28      use couple
29      use inform
30      use iounit
31      use kangs
32      use kbonds
33      use kcflux
34      use keys
35      use potent
36      use usage
37      implicit none
38      integer i,j,k
39      integer ia,ib,ic
40      integer ita,itb,itc
41      integer na,nb
42      integer size,next
43      real*8 fc,bd,cfb
44      real*8 cfa1,cfa2
45      real*8 cfb1,cfb2
46      logical headerb
47      logical headera
48      character*4 pa,pb,pc
49      character*8 blank8
50      character*8 pt,pt2
51      character*8 ptab,ptbc
52      character*12 blank12,pt3
53      character*20 keyword
54      character*240 record
55      character*240 string
56c
57c
58c     process keywords containing charge flux parameters
59c
60      blank8 = '        '
61      blank12 = '            '
62      size = 4
63      headerb = .true.
64      headera = .true.
65      do i = 1, nkey
66         next = 1
67         record = keyline(i)
68         call gettext (record,keyword,next)
69         call upcase (keyword)
70         if (keyword(1:9) .eq. 'BNDCFLUX ') then
71            ia = 0
72            ib = 0
73            cfb = 0.0d0
74            string = record(next:240)
75            read (string,*,err=10,end=10) ia,ib,cfb
76   10       continue
77            if (headerb .and. .not.silent) then
78               headerb = .false.
79               write (iout,20)
80   20          format (/,' Additional Bond Charge Flux Parameters :',
81     &                 //,5x,'Atom Classes',19x,'K(CFB)',/)
82            end if
83            if (.not. silent) then
84               write (iout,30)  ia,ib,cfb
85   30          format (6x,2i4,13x,f15.6)
86            end if
87            call numeral (ia,pa,size)
88            call numeral (ib,pb,size)
89            if (ia .le. ib) then
90               pt2 = pa//pb
91            else
92               pt2 = pb//pa
93            end if
94            do j = 1, maxncfb
95               if (kcfb(j).eq.blank8 .or. kcfb(j).eq.pt2) then
96                  kcfb(j) = pt2
97                  if (ia .lt. ib) then
98                     cflb(j) = cfb
99                  else if (ib .lt. ia) then
100                     cflb(j) = -cfb
101                  else
102                     cflb(j) = 0.0d0
103                     write (iout,40)
104   40                format (/,' KCHGFLX  --  Bond Charge Flux for',
105     &                          ' Identical Classes Set to Zero')
106                  end if
107                  goto 50
108               end if
109            end do
110   50       continue
111         else if (keyword(1:9) .eq. 'ANGCFLUX ') then
112            ia = 0
113            ib = 0
114            ic = 0
115            cfa1 = 0.0d0
116            cfa2 = 0.0d0
117            cfb1 = 0.0d0
118            cfb2 = 0.0d0
119            string = record(next:240)
120            read (string,*,err=60,end=60) ia,ib,ic,cfa1,cfa2,cfb1,cfb2
121   60       continue
122            if (headera .and. .not.silent) then
123               headera = .false.
124               write (iout,70)
125   70          format (/,' Additional Angle Charge Flux Parameters :',
126     &                 //,5x,'Atom Classes',10x,'K(CFA1)',
127     &                    7x,'K(CFA2)',7x,'K(CFB1)',7x,'K(CFB2)',/)
128            end if
129            if (.not. silent) then
130               write (iout,80)  ia,ib,ic,cfa1,cfa2,cfb1,cfb2
131   80          format (4x,3i4,4x,4f14.6)
132            end if
133            call numeral (ia,pa,size)
134            call numeral (ib,pb,size)
135            call numeral (ic,pc,size)
136            if (ia .le. ic) then
137               pt3 = pa//pb//pc
138            else
139               pt3 = pc//pb//pa
140            end if
141            do j = 1, maxncfa
142               if (kcfa(j).eq.blank12 .or. kcfa(j).eq.pt3) then
143                  kcfa(j) = pt3
144                  cfla(1,j) = cfa1
145                  cfla(2,j) = cfa2
146                  cflab(1,j) = cfb1
147                  cflab(2,j) = cfb2
148                  goto 90
149               end if
150            end do
151   90       continue
152         end if
153      end do
154c
155c     determine the total number of forcefield parameters
156c
157      nb = maxncfb
158      do i = maxncfb, 1, -1
159         if (kcfb(i) .eq. blank8)  nb = i - 1
160      end do
161      na = maxncfa
162      do i = maxncfa, 1, -1
163         if (kcfa(i) .eq. blank12)  na = i - 1
164      end do
165c
166c     perform dynamic allocation of some global arrays
167c
168      if (allocated(bflx))  deallocate (bflx)
169      if (allocated(aflx))  deallocate (aflx)
170      if (allocated(abflx))  deallocate (abflx)
171      allocate (bflx(nbond))
172      allocate (aflx(2,nangle))
173      allocate (abflx(2,nangle))
174c
175c     assign bond charge flux parameters for each bond
176c
177      nbflx = 0
178      do i = 1, nbond
179         ia = ibnd(1,i)
180         ib = ibnd(2,i)
181         ita = class(ia)
182         itb = class(ib)
183         size = 4
184         call numeral (ita,pa,size)
185         call numeral (itb,pb,size)
186         if (ita .le. itb) then
187            pt2 = pa//pb
188         else
189            pt2 = pb//pa
190         end if
191         bflx(i) = 0.0d0
192         do j = 1, nb
193            if (kcfb(j) .eq. pt2) then
194               nbflx = nbflx + 1
195               if (ita .le. itb) then
196                  bflx(i) = cflb(j)
197               else
198                  bflx(i) = -cflb(j)
199               end if
200            end if
201         end do
202      end do
203c
204c    assign angle charge flux parameters for each angle
205c
206      naflx = 0
207      do i = 1, nangle
208         ia = iang(1,i)
209         ib = iang(2,i)
210         ic = iang(3,i)
211         ita = class(ia)
212         itb = class(ib)
213         itc = class(ic)
214         call numeral (ita,pa,size)
215         call numeral (itb,pb,size)
216         call numeral (itc,pc,size)
217         if (ita .le. itc) then
218            pt3 = pa//pb//pc
219         else
220            pt3 = pc//pb//pa
221         end if
222         if (ita .le. itb) then
223            ptab = pa//pb
224         else
225            ptab = pb//pa
226         end if
227         if (itb .le. itc) then
228            ptbc = pb//pc
229         else
230            ptbc = pc//pb
231         end if
232         aflx(1,i) = 0.0d0
233         aflx(2,i) = 0.0d0
234         abflx(1,i) = 0.0d0
235         abflx(2,i) = 0.0d0
236         do j = 1, na
237            if (kcfa(j) .eq. pt3) then
238               naflx = naflx + 1
239               aflx(1,i) = cfla(1,j)
240               aflx(2,i) = cfla(2,j)
241               abflx(1,i) = cflab(1,j)
242               abflx(2,i) = cflab(2,j)
243            end if
244         end do
245      end do
246c
247c     turn off bond and angle charge flux term if not used
248c
249      if (nbflx.eq.0 .and. naflx.eq.0)  use_chgflx = .false.
250      if (.not.use_charge .and. .not.use_mpole
251     &        .and. .not.use_polar)  use_chgflx = .false.
252      return
253      end
254