1c 2c 3c ################################################### 4c ## COPYRIGHT (C) 1993 by Jay William Ponder ## 5c ## All Rights Reserved ## 6c ################################################### 7c 8c ############################################################### 9c ## ## 10c ## subroutine kopbend -- out-of-plane bending parameters ## 11c ## ## 12c ############################################################### 13c 14c 15c "kopbend" assigns the force constants for out-of-plane bends 16c at trigonal centers via Wilson-Decius-Cross or Allinger angles; 17c also processes any new or changed parameter values 18c 19c 20 subroutine kopbend 21 use angbnd 22 use atomid 23 use atoms 24 use couple 25 use fields 26 use inform 27 use iounit 28 use keys 29 use kopbnd 30 use opbend 31 use potent 32 use usage 33 implicit none 34 integer i,j,it 35 integer ia,ib,ic,id 36 integer ita,itb,itc,itd 37 integer nopb,size 38 integer next,number 39 real*8 fopb 40 logical header,done 41 logical, allocatable :: jopb(:) 42 character*4 pa,pb,pc,pd 43 character*4 zero4 44 character*8 zero8 45 character*16 blank,pt 46 character*16 pt0,pt1 47 character*20 keyword 48 character*240 record 49 character*240 string 50c 51c 52c process keywords containing out-of-plane bend parameters 53c 54 blank = ' ' 55 zero4 = '0000' 56 zero8 = '00000000' 57 header = .true. 58 do i = 1, nkey 59 next = 1 60 record = keyline(i) 61 call gettext (record,keyword,next) 62 call upcase (keyword) 63 if (keyword(1:7) .eq. 'OPBEND ') then 64 ia = 0 65 ib = 0 66 ic = 0 67 id = 0 68 fopb = 0.0d0 69 string = record(next:240) 70 read (string,*,err=10,end=10) ia,ib,ic,id,fopb 71 10 continue 72 size = 4 73 call numeral (ia,pa,size) 74 call numeral (ib,pb,size) 75 call numeral (ic,pc,size) 76 call numeral (id,pd,size) 77 if (ic .le. id) then 78 pt = pa//pb//pc//pd 79 else 80 pt = pa//pb//pd//pc 81 end if 82 if (.not. silent) then 83 if (header) then 84 header = .false. 85 write (iout,20) 86 20 format (/,' Additional Out-of-Plane Bend', 87 & ' Parameters :', 88 & //,5x,'Atom Classes',19x,'K(OPB)',/) 89 end if 90 write (iout,30) ia,ib,ic,id,fopb 91 30 format (4x,4i4,10x,f12.3) 92 end if 93 size = 4 94 do j = 1, maxnopb 95 if (kopb(j).eq.blank .or. kopb(j).eq.pt) then 96 kopb(j) = pt 97 opbn(j) = fopb 98 goto 50 99 end if 100 end do 101 write (iout,40) 102 40 format (/,' KOPBEND -- Too many Out-of-Plane', 103 & ' Angle Bending Parameters') 104 abort = .true. 105 50 continue 106 end if 107 end do 108c 109c perform dynamic allocation of some global arrays 110c 111 if (allocated(iopb)) deallocate (iopb) 112 if (allocated(opbk)) deallocate (opbk) 113 allocate (iopb(nangle)) 114 allocate (opbk(nangle)) 115c 116c use special out-of-plane bend parameter assignment for MMFF 117c 118 if (forcefield .eq. 'MMFF94') then 119 call kopbendm 120 return 121 end if 122c 123c determine the total number of forcefield parameters 124c 125 nopb = maxnopb 126 do i = maxnopb, 1, -1 127 if (kopb(i) .eq. blank) nopb = i - 1 128 end do 129c 130c perform dynamic allocation of some local arrays 131c 132 allocate (jopb(maxclass)) 133c 134c make list of atom classes using out-of-plane bending 135c 136 do i = 1, maxclass 137 jopb(i) = .false. 138 end do 139 do i = 1, maxnopb 140 if (kopb(i) .eq. blank) goto 60 141 it = number(kopb(i)(5:8)) 142 jopb(it) = .true. 143 end do 144 60 continue 145c 146c assign out-of-plane bending parameters for each angle 147c 148 nopbend = 0 149 if (nopb .ne. 0) then 150 header = .true. 151 do i = 1, nangle 152 ib = iang(2,i) 153 itb = class(ib) 154 if (jopb(itb) .and. n12(ib).eq.3) then 155 ia = iang(1,i) 156 ita = class(ia) 157 ic = iang(3,i) 158 itc = class(ic) 159 id = iang(4,i) 160 itd = class(id) 161 size = 4 162 call numeral (ita,pa,size) 163 call numeral (itb,pb,size) 164 call numeral (itc,pc,size) 165 call numeral (itd,pd,size) 166 if (ita .le. itc) then 167 pt = pd//pb//pa//pc 168 else 169 pt = pd//pb//pc//pa 170 end if 171 pt1 = pd//pb//zero8 172 pt0 = zero4//pb//zero8 173 done = .false. 174 do j = 1, nopb 175 if (kopb(j) .eq. pt) then 176 nopbend = nopbend + 1 177 iopb(nopbend) = i 178 opbk(nopbend) = opbn(j) 179 done = .true. 180 goto 70 181 end if 182 end do 183 do j = 1, nopb 184 if (kopb(j) .eq. pt1) then 185 nopbend = nopbend + 1 186 iopb(nopbend) = i 187 opbk(nopbend) = opbn(j) 188 done = .true. 189 goto 70 190 end if 191 end do 192 do j = 1, nopb 193 if (kopb(j) .eq. pt0) then 194 nopbend = nopbend + 1 195 iopb(nopbend) = i 196 opbk(nopbend) = opbn(j) 197 done = .true. 198 goto 70 199 end if 200 end do 201 70 continue 202 if (use_opbend .and. .not.done) then 203 if (use(ia) .or. use(ib) .or. use(ic) .or. use(id)) 204 & abort = .true. 205 if (header) then 206 header = .false. 207 write (iout,80) 208 80 format (/,' Undefined Out-of-Plane Bend', 209 & ' Parameters :', 210 & //,' Type',24x,'Atom Names',24x, 211 & 'Atom Classes',/) 212 end if 213 write (iout,90) id,name(id),ib,name(ib),ia,name(ia), 214 & ic,name(ic),itd,itb,ita,itc 215 90 format (' Angle-OP',3x,4(i6,'-',a3),5x,4i5) 216 end if 217 else 218 iang(4,i) = ib 219 end if 220 end do 221 end if 222c 223c perform deallocation of some local arrays 224c 225 deallocate (jopb) 226c 227c turn off the out-of-plane bending term if it is not used 228c 229 if (nopbend .eq. 0) use_opbend = .false. 230 return 231 end 232c 233c 234c ################################################################## 235c ## ## 236c ## subroutine kopbendm -- MMFF out-of-plane bend parameters ## 237c ## ## 238c ################################################################## 239c 240c 241c "kopbendm" assigns the force constants for out-of-plane bends 242c according to the Merck Molecular Force Field (MMFF) 243c 244c 245 subroutine kopbendm 246 use angbnd 247 use atomid 248 use atoms 249 use kopbnd 250 use merck 251 use opbend 252 use potent 253 implicit none 254 integer i,j,m 255 integer nopb,size 256 integer ia,ib,ic,id 257 integer ita,itb,itc,itd 258 integer itta,ittb 259 integer ittc,ittd 260 character*4 pa,pb,pc,pd 261 character*16 blank,pt 262c 263c 264c determine the total number of forcefield parameters 265c 266 blank = ' ' 267 nopb = maxnopb 268 do i = maxnopb, 1, -1 269 if (kopb(i) .eq. blank) nopb = i - 1 270 end do 271c 272c assign MMFF out-of-plane bending parameter values 273c 274 nopbend = 0 275 if (nopb .ne. 0) then 276 do i = 1, nangle 277 ia = iang(1,i) 278 ib = iang(2,i) 279 ic = iang(3,i) 280 id = iang(4,i) 281 if (min(ia,ib,ic,id) .gt. 0) then 282 itta = type(ia) 283 ittb = type(ib) 284 ittc = type(ic) 285 ittd = type(id) 286 m = 0 287 10 continue 288 m = m + 1 289 if (m .eq. 1) then 290 ita = eqclass(itta,1) 291 itb = eqclass(ittb,1) 292 itc = eqclass(ittc,1) 293 itd = eqclass(ittd,1) 294 else if (m .eq. 2) then 295 ita = eqclass(itta,2) 296 itb = eqclass(ittb,2) 297 itc = eqclass(ittc,2) 298 itd = eqclass(ittd,2) 299 else if (m .eq. 3) then 300 ita = eqclass(itta,3) 301 itb = eqclass(ittb,2) 302 itc = eqclass(ittc,3) 303 itd = eqclass(ittd,3) 304 else if (m .eq. 4) then 305 ita = eqclass(itta,4) 306 itb = eqclass(ittb,2) 307 itc = eqclass(ittc,4) 308 itd = eqclass(ittd,4) 309 else if (m .eq. 5) then 310 ita = eqclass(itta,5) 311 itb = eqclass(ittb,2) 312 itc = eqclass(ittc,5) 313 itd = eqclass(ittd,5) 314 end if 315 if (m .gt. 5) then 316 nopbend = nopbend + 1 317 iopb(nopbend) = i 318 opbk(nopbend) = 0.0d0 319 else 320 size = 4 321 call numeral (ita,pa,size) 322 call numeral (itb,pb,size) 323 call numeral (itc,pc,size) 324 call numeral (itd,pd,size) 325 if (itd.le.ita .and. itd.le.itc) then 326 if (ita .le. itc) then 327 pt = pd//pb//pa//pc 328 else 329 pt = pd//pb//pc//pa 330 end if 331 else if (ita.le.itc .and. ita.le.itd) then 332 if (itd .le. itc) then 333 pt = pa//pb//pd//pc 334 else 335 pt = pa//pb//pc//pd 336 end if 337 else if (itc.le.ita .and. itc.le.itd) then 338 if (ita .le. itd) then 339 pt = pc//pb//pa//pd 340 else 341 pt = pc//pb//pd//pa 342 end if 343 end if 344 do j = 1, nopb 345 if (kopb(j) .eq. pt) then 346 nopbend = nopbend + 1 347 iopb(nopbend) = i 348 opbk(nopbend) = opbn(j) 349 goto 20 350 end if 351 end do 352 if (class(ib).eq.8 .or. class(ib).eq.17 .or. 353 & class(ib).eq.26 .or. class(ib).eq.43 .or. 354 & class(ib).eq.49 .or. class(ib).eq.73 .or. 355 & class(ib).eq.82) then 356 nopbend = nopbend + 1 357 iopb(nopbend) = i 358 opbk(nopbend) = 0.0d0 359 goto 20 360 end if 361 goto 10 362 20 continue 363 end if 364 end if 365 end do 366 end if 367c 368c turn off the out-of-plane bending term if it is not used 369c 370 if (nopbend .eq. 0) use_opbend = .false. 371 return 372 end 373