1c 2c 3c ################################################### 4c ## COPYRIGHT (C) 1990 by Jay William Ponder ## 5c ## All Rights Reserved ## 6c ################################################### 7c 8c ############################################################# 9c ## ## 10c ## subroutine kdipole -- assign bond dipole parameters ## 11c ## ## 12c ############################################################# 13c 14c 15c "kdipole" assigns bond dipoles to the bonds within 16c the structure and processes any new or changed values 17c 18c 19 subroutine kdipole 20 use atmlst 21 use atoms 22 use bndstr 23 use couple 24 use dipole 25 use inform 26 use iounit 27 use kdipol 28 use keys 29 use potent 30 implicit none 31 integer i,j,k 32 integer ia,ib,ita,itb 33 integer nd,nd5,nd4,nd3 34 integer iring,size,next 35 real*8 dp,ps 36 logical header 37 logical use_ring 38 character*4 pa,pb 39 character*6 label 40 character*8 blank,pt 41 character*20 keyword 42 character*240 record 43 character*240 string 44c 45c 46c process keywords containing bond dipole parameters 47c 48 blank = ' ' 49 header = .true. 50 do i = 1, nkey 51 next = 1 52 record = keyline(i) 53 call gettext (record,keyword,next) 54 call upcase (keyword) 55 iring = -1 56 if (keyword(1:7) .eq. 'DIPOLE ') iring = 0 57 if (keyword(1:8) .eq. 'DIPOLE5 ') iring = 5 58 if (keyword(1:8) .eq. 'DIPOLE4 ') iring = 4 59 if (keyword(1:8) .eq. 'DIPOLE3 ') iring = 3 60 if (iring .ge. 0) then 61 ia = 0 62 ib = 0 63 dp = 0.0d0 64 ps = 0.5d0 65 string = record(next:240) 66 read (string,*,err=10,end=10) ia,ib,dp,ps 67 10 continue 68 if (ia.gt.0 .and. ib.gt.0) then 69 if (.not. silent) then 70 if (header) then 71 header = .false. 72 write (iout,20) 73 20 format (/,' Additional Bond Dipole Moment ', 74 & 'Parameters :', 75 & //,5x,'Atom Types',13x,'Moment', 76 & 8x,'Position',/) 77 end if 78 if (iring .eq. 0) then 79 write (iout,30) ia,ib,dp,ps 80 30 format (6x,2i4,5x,2f15.3) 81 else 82 if (iring .eq. 5) label = '5-Ring' 83 if (iring .eq. 4) label = '4-Ring' 84 if (iring .eq. 3) label = '3-Ring' 85 write (iout,40) ia,ib,dp,ps,label 86 40 format (6x,2i4,5x,2f15.3,3x,a6) 87 end if 88 end if 89 size = 4 90 call numeral (ia,pa,size) 91 call numeral (ib,pb,size) 92 if (ia .le. ib) then 93 pt = pa//pb 94 else 95 pt = pb//pa 96 end if 97 if (iring .eq. 0) then 98 do j = 1, maxnd 99 if (kd(j).eq.blank .or. kd(j).eq.pt) then 100 kd(j) = pt 101 if (ia .le. ib) then 102 dpl(j) = dp 103 pos(j) = ps 104 else 105 dpl(j) = -dp 106 pos(j) = 1.0d0 - ps 107 end if 108 goto 90 109 end if 110 end do 111 write (iout,50) 112 50 format (/,' KDIPOLE -- Too many Bond Dipole', 113 & ' Moment Parameters') 114 abort = .true. 115 else if (iring .eq. 5) then 116 do j = 1, maxnd5 117 if (kd5(j).eq.blank .or. kd5(j).eq.pt) then 118 kd5(j) = pt 119 if (ia .le. ib) then 120 dpl5(j) = dp 121 pos5(j) = ps 122 else 123 dpl5(j) = -dp 124 pos5(j) = 1.0d0 - ps 125 end if 126 goto 90 127 end if 128 end do 129 write (iout,60) 130 60 format (/,' KDIPOLE -- Too many 5-Ring Bond', 131 & ' Dipole Parameters') 132 abort = .true. 133 else if (iring .eq. 4) then 134 do j = 1, maxnd4 135 if (kd4(j).eq.blank .or. kd4(j).eq.pt) then 136 kd4(j) = pt 137 if (ia .le. ib) then 138 dpl4(j) = dp 139 pos4(j) = ps 140 else 141 dpl4(j) = -dp 142 pos4(j) = 1.0d0 - ps 143 end if 144 goto 90 145 end if 146 end do 147 write (iout,70) 148 70 format (/,' KDIPOLE -- Too many 4-Ring Bond', 149 & ' Dipole Parameters') 150 abort = .true. 151 else if (iring .eq. 3) then 152 do j = 1, maxnd3 153 if (kd3(j).eq.blank .or. kd3(j).eq.pt) then 154 kd3(j) = pt 155 if (ia .le. ib) then 156 dpl3(j) = dp 157 pos3(j) = ps 158 else 159 dpl3(j) = -dp 160 pos3(j) = 1.0d0 - ps 161 end if 162 goto 90 163 end if 164 end do 165 write (iout,80) 166 80 format (/,' KDIPOLE -- Too many 3-Ring Bond', 167 & ' Dipole Parameters') 168 abort = .true. 169 end if 170 end if 171 90 continue 172 end if 173 end do 174c 175c determine the total number of forcefield parameters 176c 177 nd = maxnd 178 nd5 = maxnd5 179 nd4 = maxnd4 180 nd3 = maxnd3 181 do i = maxnd, 1, -1 182 if (kd(i) .eq. blank) nd = i - 1 183 end do 184 do i = maxnd5, 1, -1 185 if (kd5(i) .eq. blank) nd5 = i - 1 186 end do 187 do i = maxnd4, 1, -1 188 if (kd4(i) .eq. blank) nd4 = i - 1 189 end do 190 do i = maxnd3, 1, -1 191 if (kd3(i) .eq. blank) nd3 = i - 1 192 end do 193 use_ring = .false. 194 if (min(nd5,nd4,nd3) .ne. 0) use_ring = .true. 195c 196c perform dynamic allocation of some global arrays 197c 198 if (allocated(idpl)) deallocate (idpl) 199 if (allocated(bdpl)) deallocate (bdpl) 200 if (allocated(sdpl)) deallocate (sdpl) 201 allocate (idpl(2,nbond)) 202 allocate (bdpl(nbond)) 203 allocate (sdpl(nbond)) 204c 205c find and store all the bond dipole moments 206c 207 do i = 1, nbond 208 ia = ibnd(1,i) 209 ib = ibnd(2,i) 210 ita = type(ia) 211 itb = type(ib) 212 size = 4 213 call numeral (ita,pa,size) 214 call numeral (itb,pb,size) 215 if (ita .le. itb) then 216 pt = pa//pb 217 else 218 pt = pb//pa 219 end if 220 bdpl(i) = 0.0d0 221c 222c make a check for bonds contained inside small rings 223c 224 iring = 0 225 if (use_ring) then 226 call chkring (iring,ia,ib,0,0) 227 if (iring .eq. 6) iring = 0 228 if (iring.eq.5 .and. nd5.eq.0) iring = 0 229 if (iring.eq.4 .and. nd4.eq.0) iring = 0 230 if (iring.eq.3 .and. nd3.eq.0) iring = 0 231 end if 232c 233c try to assign bond dipole parameters for the bond 234c 235 if (iring .eq. 0) then 236 do j = 1, nd 237 if (kd(j) .eq. pt) then 238 if (ita .le. itb) then 239 idpl(1,i) = ia 240 idpl(2,i) = ib 241 else 242 idpl(1,i) = ib 243 idpl(2,i) = ia 244 end if 245 bdpl(i) = dpl(j) 246 sdpl(i) = pos(j) 247 goto 100 248 end if 249 end do 250 else if (iring .eq. 5) then 251 do j = 1, nd5 252 if (kd5(j) .eq. pt) then 253 if (ita .le. itb) then 254 idpl(1,i) = ia 255 idpl(2,i) = ib 256 else 257 idpl(1,i) = ib 258 idpl(2,i) = ia 259 end if 260 bdpl(i) = dpl5(j) 261 sdpl(i) = pos5(j) 262 goto 100 263 end if 264 end do 265 else if (iring .eq. 4) then 266 do j = 1, nd4 267 if (kd4(j) .eq. pt) then 268 if (ita .le. itb) then 269 idpl(1,i) = ia 270 idpl(2,i) = ib 271 else 272 idpl(1,i) = ib 273 idpl(2,i) = ia 274 end if 275 bdpl(i) = dpl4(j) 276 sdpl(i) = pos4(j) 277 goto 100 278 end if 279 end do 280 else if (iring .eq. 3) then 281 do j = 1, nd3 282 if (kd3(j) .eq. pt) then 283 if (ita .le. itb) then 284 idpl(1,i) = ia 285 idpl(2,i) = ib 286 else 287 idpl(1,i) = ib 288 idpl(2,i) = ia 289 end if 290 bdpl(i) = dpl3(j) 291 sdpl(i) = pos3(j) 292 goto 100 293 end if 294 end do 295 end if 296 100 continue 297 end do 298c 299c process keywords containing bond specific bond dipoles 300c 301 header = .true. 302 do i = 1, nkey 303 next = 1 304 record = keyline(i) 305 call gettext (record,keyword,next) 306 call upcase (keyword) 307 if (keyword(1:7) .eq. 'DIPOLE ') then 308 ia = 0 309 ib = 0 310 dp = 0.0d0 311 ps = 0.0d0 312 string = record(next:240) 313 read (string,*,err=110,end=110) ia,ib,dp,ps 314 110 continue 315 if (ia.lt.0 .and. ib.lt.0) then 316 ia = -ia 317 ib = -ib 318 if (header .and. .not.silent) then 319 header = .false. 320 write (iout,120) 321 120 format (/,' Additional Bond Dipoles for', 322 & ' Specific Bonds :', 323 & //,5x,'Bonded Atoms',11x,'Moment', 324 & 8x,'Position',/) 325 end if 326 do j = 1, n12(ia) 327 if (i12(j,ia) .eq. ib) then 328 k = bndlist(j,ia) 329 if (ps .eq. 0.0d0) ps = 0.5d0 330 if (idpl(1,k) .eq. ib) then 331 bdpl(k) = dp 332 sdpl(k) = ps 333 else 334 bdpl(k) = -dp 335 sdpl(k) = 1.0d0 - ps 336 end if 337 if (.not. silent) then 338 write (iout,130) ia,ib,dp,ps 339 130 format (4x,i5,' -',i5,3x,2f15.3) 340 end if 341 goto 140 342 end if 343 end do 344 end if 345 140 continue 346 end if 347 end do 348c 349c remove zero bond dipoles from the list of dipoles 350c 351 ndipole = 0 352 do i = 1, nbond 353 if (bdpl(i) .ne. 0.0d0) then 354 ndipole = ndipole + 1 355 idpl(1,ndipole) = idpl(1,i) 356 idpl(2,ndipole) = idpl(2,i) 357 bdpl(ndipole) = bdpl(i) 358 sdpl(ndipole) = sdpl(i) 359 end if 360 end do 361c 362c turn off dipole-dipole and charge-dipole terms if not used 363c 364 if (ndipole .eq. 0) then 365 use_dipole = .false. 366 use_chgdpl = .false. 367 end if 368 return 369 end 370