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