1c 2c 3c ################################################### 4c ## COPYRIGHT (C) 1990 by Jay William Ponder ## 5c ## All Rights Reserved ## 6c ################################################### 7c 8c ############################################################## 9c ## ## 10c ## subroutine kangle -- angle bend parameter assignment ## 11c ## ## 12c ############################################################## 13c 14c 15c "kangle" assigns the force constants and ideal angles for 16c the bond angles; also processes new or changed parameters 17c 18c 19 subroutine kangle 20 use angbnd 21 use angpot 22 use atomid 23 use atoms 24 use couple 25 use fields 26 use inform 27 use iounit 28 use kangs 29 use keys 30 use potent 31 use usage 32 implicit none 33 integer i,j 34 integer ia,ib,ic 35 integer ita,itb,itc 36 integer na,nap,naf 37 integer na3,na4,na5 38 integer jen,ih,nh 39 integer next,size 40 integer minat,iring 41 real*8 fc,an,pr 42 real*8 an1,an2,an3 43 logical header,done 44 logical use_ring 45 character*4 pa,pb,pc 46 character*6 label 47 character*12 blank,pt 48 character*20 keyword 49 character*240 record 50 character*240 string 51c 52c 53c process keywords containing angle bending parameters 54c 55 blank = ' ' 56 header = .true. 57 do i = 1, nkey 58 next = 1 59 record = keyline(i) 60 call gettext (record,keyword,next) 61 call upcase (keyword) 62 iring = -1 63 if (keyword(1:6) .eq. 'ANGLE ') iring = 0 64 if (keyword(1:7) .eq. 'ANGLE5 ') iring = 5 65 if (keyword(1:7) .eq. 'ANGLE4 ') iring = 4 66 if (keyword(1:7) .eq. 'ANGLE3 ') iring = 3 67 if (iring .ge. 0) then 68 ia = 0 69 ib = 0 70 ic = 0 71 fc = 0.0d0 72 an1 = 0.0d0 73 an2 = 0.0d0 74 an3 = 0.0d0 75 jen = 0 76 string = record(next:240) 77 read (string,*,err=10,end=10) ia,ib,ic,fc,an1,an2,an3 78 10 continue 79 if (an2.ne.0.0d0 .or. an3.ne.0.0d0) jen = 1 80 if (.not. silent) then 81 if (header) then 82 header = .false. 83 write (iout,20) 84 20 format (/,' Additional Angle Bending Parameters :', 85 & //,5x,'Atom Classes',13x,'K(B)',10x,'Angle',/) 86 end if 87 if (iring .eq. 0) then 88 if (jen .eq. 0) then 89 write (iout,30) ia,ib,ic,fc,an1 90 30 format (4x,3i4,3x,2f15.3) 91 else if (an1 .ne. 0.0d0) then 92 write (iout,40) ia,ib,ic,fc,an1 93 40 format (4x,3i4,3x,2f15.3,3x,'0-H''s') 94 end if 95 if (an2 .ne. 0.0d0) then 96 write (iout,50) ia,ib,ic,fc,an2 97 50 format (4x,3i4,3x,2f15.3,3x,'1-H''s') 98 end if 99 if (an3 .ne. 0.0d0) then 100 write (iout,60) ia,ib,ic,fc,an3 101 60 format (4x,3i4,3x,2f15.3,3x,'2-H''s') 102 end if 103 else 104 if (iring .eq. 5) label = '5-Ring' 105 if (iring .eq. 4) label = '4-Ring' 106 if (iring .eq. 3) label = '3-Ring' 107 if (jen .eq. 0) then 108 write (iout,70) ia,ib,ic,fc,an1,label 109 70 format (4x,3i4,3x,2f15.3,3x,a6) 110 else if (an1 .ne. 0.0d0) then 111 write (iout,80) ia,ib,ic,fc,an1,label 112 80 format (4x,3i4,3x,2f15.3,3x,a6,3x,'0-H''s') 113 end if 114 if (an2 .ne. 0.0d0) then 115 write (iout,90) ia,ib,ic,fc,an2,label 116 90 format (4x,3i4,3x,2f15.3,3x,a6,3x,'1-H''s') 117 end if 118 if (an3 .ne. 0.0d0) then 119 write (iout,100) ia,ib,ic,fc,an3,label 120 100 format (4x,3i4,3x,2f15.3,3x,a6,3x,'2-H''s') 121 end if 122 end if 123 end if 124 size = 4 125 call numeral (ia,pa,size) 126 call numeral (ib,pb,size) 127 call numeral (ic,pc,size) 128 if (ia .le. ic) then 129 pt = pa//pb//pc 130 else 131 pt = pc//pb//pa 132 end if 133 if (iring .eq. 0) then 134 do j = 1, maxna 135 if (ka(j).eq.blank .or. ka(j).eq.pt) then 136 ka(j) = pt 137 acon(j) = fc 138 ang(1,j) = an1 139 ang(2,j) = an2 140 ang(3,j) = an3 141 goto 150 142 end if 143 end do 144 write (iout,110) 145 110 format (/,' KANGLE -- Too many Bond Angle', 146 & ' Bending Parameters') 147 abort = .true. 148 else if (iring .eq. 5) then 149 do j = 1, maxna5 150 if (ka5(j).eq.blank .or. ka5(j).eq.pt) then 151 ka5(j) = pt 152 acon5(j) = fc 153 ang5(1,j) = an1 154 ang5(2,j) = an2 155 ang5(3,j) = an3 156 goto 150 157 end if 158 end do 159 write (iout,120) 160 120 format (/,' KANGLE -- Too many 5-Ring Angle', 161 & ' Bending Parameters') 162 abort = .true. 163 else if (iring .eq. 4) then 164 do j = 1, maxna4 165 if (ka4(j).eq.blank .or. ka4(j).eq.pt) then 166 ka4(j) = pt 167 acon4(j) = fc 168 ang4(1,j) = an1 169 ang4(2,j) = an2 170 ang4(3,j) = an3 171 goto 150 172 end if 173 end do 174 write (iout,130) 175 130 format (/,' KANGLE -- Too many 4-Ring Angle', 176 & ' Bending Parameters') 177 abort = .true. 178 else if (iring .eq. 3) then 179 do j = 1, maxna3 180 if (ka3(j).eq.blank .or. ka3(j).eq.pt) then 181 ka3(j) = pt 182 acon3(j) = fc 183 ang3(1,j) = an1 184 ang3(2,j) = an2 185 ang3(3,j) = an3 186 goto 150 187 end if 188 end do 189 write (iout,140) 190 140 format (/,' KANGLE -- Too many 3-Ring Angle', 191 & ' Bending Parameters') 192 abort = .true. 193 end if 194 150 continue 195 end if 196 end do 197c 198c process keywords containing in-plane angle bending parameters 199c 200 header = .true. 201 do i = 1, nkey 202 next = 1 203 record = keyline(i) 204 call gettext (record,keyword,next) 205 call upcase (keyword) 206 if (keyword(1:7) .eq. 'ANGLEP ') then 207 ia = 0 208 ib = 0 209 ic = 0 210 fc = 0.0d0 211 an1 = 0.0d0 212 an2 = 0.0d0 213 jen = 0 214 string = record(next:240) 215 read (string,*,err=160,end=160) ia,ib,ic,fc,an1,an2 216 160 continue 217 if (an2 .ne. 0.0d0) jen = 1 218 if (.not. silent) then 219 if (header) then 220 header = .false. 221 write (iout,170) 222 170 format (/,' Additional In-Plane Angle Bending', 223 & ' Parameters :', 224 & //,5x,'Atom Classes',13x,'K(B)',10x,'Angle',/) 225 end if 226 if (jen .eq. 0) then 227 write (iout,180) ia,ib,ic,fc,an1 228 180 format (4x,3i4,3x,2f15.3) 229 else if (an1 .ne. 0.0d0) then 230 write (iout,190) ia,ib,ic,fc,an1 231 190 format (4x,3i4,3x,2f15.3,3x,'0-H''s') 232 end if 233 if (an2 .ne. 0.0d0) then 234 write (iout,200) ia,ib,ic,fc,an2 235 200 format (4x,3i4,3x,2f15.3,3x,'1-H''s') 236 end if 237 end if 238 size = 4 239 call numeral (ia,pa,size) 240 call numeral (ib,pb,size) 241 call numeral (ic,pc,size) 242 if (ia .le. ic) then 243 pt = pa//pb//pc 244 else 245 pt = pc//pb//pa 246 end if 247 do j = 1, maxnap 248 if (kap(j).eq.blank .or. kap(j).eq.pt) then 249 kap(j) = pt 250 aconp(j) = fc 251 angp(1,j) = an1 252 angp(2,j) = an2 253 goto 230 254 end if 255 end do 256 write (iout,220) 257 220 format (/,' KANGLE -- Too many In-Plane Angle', 258 & ' Bending Parameters') 259 abort = .true. 260 230 continue 261 end if 262 end do 263c 264c process keywords containing Fourier angle bending parameters 265c 266 header = .true. 267 do i = 1, nkey 268 next = 1 269 record = keyline(i) 270 call gettext (record,keyword,next) 271 call upcase (keyword) 272 if (keyword(1:7) .eq. 'ANGLEF ') then 273 ia = 0 274 ib = 0 275 ic = 0 276 fc = 0.0d0 277 an = 0.0d0 278 pr = 0.0d0 279 string = record(next:240) 280 read (string,*,err=240,end=240) ia,ib,ic,fc,an,pr 281 240 continue 282 if (.not. silent) then 283 if (header) then 284 header = .false. 285 write (iout,250) 286 250 format (/,' Additional Fourier Angle Bending', 287 & ' Parameters :', 288 & //,5x,'Atom Classes',13x,'K(B)',10x,'Shift', 289 & 9x,'Period',/) 290 end if 291 write (iout,260) ia,ib,ic,fc,an,pr 292 260 format (4x,3i4,3x,3f15.3) 293 end if 294 size = 4 295 call numeral (ia,pa,size) 296 call numeral (ib,pb,size) 297 call numeral (ic,pc,size) 298 if (ia .le. ic) then 299 pt = pa//pb//pc 300 else 301 pt = pc//pb//pa 302 end if 303 do j = 1, maxnaf 304 if (kaf(j).eq.blank .or. kaf(j).eq.pt) then 305 kaf(j) = pt 306 aconf(j) = fc 307 angf(1,j) = an 308 angf(2,j) = pr 309 goto 280 310 end if 311 end do 312 write (iout,270) 313 270 format (/,' KANGLE -- Too many Fourier Angle', 314 & ' Bending Parameters') 315 abort = .true. 316 280 continue 317 end if 318 end do 319c 320c determine the total number of forcefield parameters 321c 322 na = maxna 323 na5 = maxna5 324 na4 = maxna4 325 na3 = maxna3 326 nap = maxnap 327 naf = maxnaf 328 do i = maxna, 1, -1 329 if (ka(i) .eq. blank) na = i - 1 330 end do 331 do i = maxna5, 1, -1 332 if (ka5(i) .eq. blank) na5 = i - 1 333 end do 334 do i = maxna4, 1, -1 335 if (ka4(i) .eq. blank) na4 = i - 1 336 end do 337 do i = maxna3, 1, -1 338 if (ka3(i) .eq. blank) na3 = i - 1 339 end do 340 do i = maxnap, 1, -1 341 if (kap(i) .eq. blank) nap = i - 1 342 end do 343 do i = maxnaf, 1, -1 344 if (kaf(i) .eq. blank) naf = i - 1 345 end do 346 use_ring = .false. 347 if (min(na5,na4,na3) .ne. 0) use_ring = .true. 348c 349c set generic parameters for use with any number of hydrogens 350c 351 do i = 1, na 352 if (ang(2,i).eq.0.0d0 .and. ang(3,i).eq.0.0d0) then 353 ang(2,i) = ang(1,i) 354 ang(3,i) = ang(1,i) 355 end if 356 end do 357 do i = 1, na5 358 if (ang5(2,i).eq.0.0d0 .and. ang5(3,i).eq.0.0d0) then 359 ang5(2,i) = ang5(1,i) 360 ang5(3,i) = ang5(1,i) 361 end if 362 end do 363 do i = 1, na4 364 if (ang4(2,i).eq.0.0d0 .and. ang4(3,i).eq.0.0d0) then 365 ang4(2,i) = ang4(1,i) 366 ang4(3,i) = ang4(1,i) 367 end if 368 end do 369 do i = 1, na3 370 if (ang3(2,i).eq.0.0d0 .and. ang3(3,i).eq.0.0d0) then 371 ang3(2,i) = ang3(1,i) 372 ang3(3,i) = ang3(1,i) 373 end if 374 end do 375 do i = 1, nap 376 if (angp(2,i) .eq. 0.0d0) then 377 angp(2,i) = angp(1,i) 378 end if 379 end do 380c 381c perform dynamic allocation of some global arrays 382c 383 if (allocated(ak)) deallocate (ak) 384 if (allocated(anat)) deallocate (anat) 385 if (allocated(afld)) deallocate (afld) 386 if (allocated(angtyp)) deallocate (angtyp) 387 allocate (ak(nangle)) 388 allocate (anat(nangle)) 389 allocate (afld(nangle)) 390 allocate (angtyp(nangle)) 391c 392c use special angle parameter assignment method for MMFF 393c 394 if (forcefield .eq. 'MMFF94') then 395 call kanglem 396 return 397 end if 398c 399c assign ideal bond angle and force constant for each angle 400c 401 header = .true. 402 do i = 1, nangle 403 ia = iang(1,i) 404 ib = iang(2,i) 405 ic = iang(3,i) 406 ita = class(ia) 407 itb = class(ib) 408 itc = class(ic) 409 size = 4 410 call numeral (ita,pa,size) 411 call numeral (itb,pb,size) 412 call numeral (itc,pc,size) 413 if (ita .le. itc) then 414 pt = pa//pb//pc 415 else 416 pt = pc//pb//pa 417 end if 418 ak(i) = 0.0d0 419 anat(i) = 0.0d0 420 afld(i) = 0.0d0 421 angtyp(i) = 'HARMONIC' 422 done = .false. 423c 424c count number of non-angle hydrogens on the central atom 425c 426 nh = 1 427 do j = 1, n12(ib) 428 ih = i12(j,ib) 429 if (ih.ne.ia .and. ih.ne.ic .and. atomic(ih).eq.1) 430 & nh = nh + 1 431 end do 432c 433c make a check for bond angles contained inside small rings 434c 435 iring = 0 436 if (use_ring) then 437 call chkring (iring,ia,ib,ic,0) 438 if (iring .eq. 6) iring = 0 439 if (iring.eq.5 .and. na5.eq.0) iring = 0 440 if (iring.eq.4 .and. na4.eq.0) iring = 0 441 if (iring.eq.3 .and. na3.eq.0) iring = 0 442 end if 443c 444c assign angle bending parameters for bond angles 445c 446 if (iring .eq. 0) then 447 do j = 1, na 448 if (ka(j).eq.pt .and. ang(nh,j).ne.0.0d0) then 449 ak(i) = acon(j) 450 anat(i) = ang(nh,j) 451 done = .true. 452 goto 290 453 end if 454 end do 455c 456c assign bending parameters for 5-membered ring angles 457c 458 else if (iring .eq. 5) then 459 do j = 1, na5 460 if (ka5(j).eq.pt .and. ang5(nh,j).ne.0.0d0) then 461 ak(i) = acon5(j) 462 anat(i) = ang5(nh,j) 463 done = .true. 464 goto 290 465 end if 466 end do 467c 468c assign bending parameters for 4-membered ring angles 469c 470 else if (iring .eq. 4) then 471 do j = 1, na4 472 if (ka4(j).eq.pt .and. ang4(nh,j).ne.0.0d0) then 473 ak(i) = acon4(j) 474 anat(i) = ang4(nh,j) 475 done = .true. 476 goto 290 477 end if 478 end do 479c 480c assign bending parameters for 3-membered ring angles 481c 482 else if (iring .eq. 3) then 483 do j = 1, na3 484 if (ka3(j).eq.pt .and. ang3(nh,j).ne.0.0d0) then 485 ak(i) = acon3(j) 486 anat(i) = ang3(nh,j) 487 done = .true. 488 goto 290 489 end if 490 end do 491 end if 492c 493c assign in-plane angle bending parameters for bond angles 494c 495 if (.not.done .and. n12(ib).eq.3) then 496 do j = 1, nap 497 if (kap(j).eq.pt .and. angp(nh,j).ne.0.0d0) then 498 ak(i) = aconp(j) 499 anat(i) = angp(nh,j) 500 angtyp(i) = 'IN-PLANE' 501 done = .true. 502 goto 290 503 end if 504 end do 505 end if 506c 507c assign Fourier angle bending parameters for bond angles 508c 509 if (.not. done) then 510 do j = 1, naf 511 if (kaf(j) .eq. pt) then 512 ak(i) = aconf(j) 513 anat(i) = angf(1,j) 514 afld(i) = angf(2,j) 515 angtyp(i) = 'FOURIER' 516 done = .true. 517 goto 290 518 end if 519 end do 520 end if 521c 522c warning if suitable angle bending parameter not found 523c 524 290 continue 525 minat = min(atomic(ia),atomic(ib),atomic(ic)) 526 if (minat .eq. 0) done = .true. 527 if (use_angle .and. .not.done) then 528 if (use(ia) .or. use(ib) .or. use(ic)) abort = .true. 529 if (header) then 530 header = .false. 531 write (iout,300) 532 300 format (/,' Undefined Angle Bending Parameters :', 533 & //,' Type',18x,'Atom Names',19x, 534 & 'Atom Classes',/) 535 end if 536 label = 'Angle ' 537 if (iring .eq. 5) label = '5-Ring' 538 if (iring .eq. 4) label = '4-Ring' 539 if (iring .eq. 3) label = '3-Ring' 540 write (iout,310) label,ia,name(ia),ib,name(ib), 541 & ic,name(ic),ita,itb,itc 542 310 format (1x,a6,5x,3(i6,'-',a3),7x,3i5) 543 end if 544 end do 545c 546c turn off the angle bending potential if it is not used 547c 548 if (nangle .eq. 0) use_angle = .false. 549 return 550 end 551c 552c 553c ############################################################### 554c ## ## 555c ## subroutine kanglem -- MMFF angle parameter assignment ## 556c ## ## 557c ############################################################### 558c 559c 560c "kanglem" assigns the force constants and ideal angles for 561c bond angles according to the Merck Molecular Force Field (MMFF) 562c 563c literature reference: 564c 565c T. A. Halgren, "Merck Molecular Force Field. I. Basis, Form, 566c Scope, Parametrization, and Performance of MMFF94", Journal of 567c Computational Chemistry, 17, 490-519 (1995) 568c 569c T. A. Halgren, "Merck Molecular Force Field. V. Extension of 570c MMFF94 Using Experimental Data, Additional Computational Data, 571c and Empirical Rules", Journal of Computational Chemistry, 17, 572c 616-641 (1995) 573c 574c 575 subroutine kanglem 576 use angbnd 577 use angpot 578 use atomid 579 use atoms 580 use bndstr 581 use merck 582 use potent 583 use ring 584 implicit none 585 integer i,j,k,l,m 586 integer ia,ib,ic 587 integer ita,itb,itc 588 integer ina,inb,inc 589 integer itta,ittb,ittc 590 integer bnd_ab,bnd_bc 591 integer at,minat 592 integer mclass 593 real*8 d,beta 594 real*8 z2(100),c(100) 595 logical done 596 logical ring3,ring4 597c 598c 599c set empirical rule parameters for some common elements 600c 601 do i = 1, 100 602 z2(i) = 1000.0d0 603 c(i) = 1000.0d0 604 end do 605 z2(1) = 1.395d0 606 z2(5) = 0.0d0 607 z2(6) = 2.494d0 608 z2(7) = 2.711d0 609 z2(8) = 3.045d0 610 z2(9) = 2.847d0 611 z2(14) = 2.350d0 612 z2(15) = 2.350d0 613 z2(16) = 2.980d0 614 z2(17) = 2.909d0 615 z2(35) = 3.017d0 616 z2(33) = 0.0d0 617 z2(53) = 3.086d0 618 c(1) = 0.0d0 619 c(5) = 0.704d0 620 c(6) = 1.016d0 621 c(7) = 1.113d0 622 c(8) = 1.337d0 623 c(9) = 0.0d0 624 c(14) = 0.811d0 625 c(15) = 1.068d0 626 c(16) = 1.249d0 627 c(17) = 1.078d0 628 c(35) = 0.0d0 629 c(33) = 0.825d0 630 c(53) = 0.0d0 631c 632c assign MMFF bond angle and force constant for each angle 633c 634 do i = 1, nangle 635 ia = iang(1,i) 636 ib = iang(2,i) 637 ic = iang(3,i) 638 ita = class(ia) 639 itb = class(ib) 640 itc = class(ic) 641 itta = type(ia) 642 ittb = type(ib) 643 ittc = type(ic) 644 ina = atomic(ia) 645 inb = atomic(ib) 646 inc = atomic(ic) 647c 648c set angle index value, accounting for MMFF bond type = 1 649c 650 at = 0 651 do j = 1, nligne 652 if ((ia.eq.bt_1(j,1) .and. ib.eq.bt_1(j,2)) .or. 653 & (ib.eq.bt_1(j,1) .and. ia.eq.bt_1(j,2))) then 654 at = at + 1 655 end if 656 if ((ic.eq.bt_1(j,1) .and. ib.eq.bt_1(j,2)) .or. 657 & (ib.eq.bt_1(j,1) .and. ic.eq.bt_1(j,2))) then 658 at = at + 1 659 end if 660 end do 661c 662c determine if the atoms belong to a 3- or 4-membered ring 663c 664 ring3 = .false. 665 ring4 = .false. 666 do j = 1, nring3 667 do k = 1, 3 668 if (ia .eq. iring3(k,j)) then 669 do l = 1, 3 670 if (ib .eq. iring3(l,j)) then 671 do m = 1, 3 672 if (ic .eq. iring3(m,j)) ring3 = .true. 673 end do 674 end if 675 end do 676 end if 677 end do 678 end do 679 if (.not. ring3) then 680 do j = 1, nring4 681 do k = 1, 4 682 if (ia .eq. iring4(k,j)) then 683 do l = 1, 4 684 if (ib .eq. iring4(l,j)) then 685 do m = 1, 4 686 if (ic .eq. iring4(m,j)) ring4 = .true. 687 end do 688 end if 689 end do 690 end if 691 end do 692 end do 693 end if 694c 695c set special index value when 3- or 4-rings are present 696c 697 if (at.eq.0 .and. ring4) then 698 at = 4 699 else if (at.eq.1 .and. ring4) then 700 at = 7 701 else if (at.eq.2 .and. ring4) then 702 at = 8 703 else if (at.eq.0 .and. ring3) then 704 at = 3 705 else if (at.eq.1 .and. ring3) then 706 at = 5 707 else if (at.eq.2 .and. ring3) then 708 at = 6 709 end if 710c 711c setup the atom class equivalencies assignment 712c 713 mclass = 0 714 10 continue 715 mclass = mclass + 1 716 if (mclass .eq. 1) then 717 ita = eqclass(itta,1) 718 itb = eqclass(ittb,1) 719 itc = eqclass(ittc,1) 720 else if (mclass .eq. 2) then 721 ita = eqclass(itta,2) 722 itb = eqclass(ittb,2) 723 itc = eqclass(ittc,2) 724 else if (mclass .eq. 3) then 725 ita = eqclass(itta,3) 726 itb = eqclass(ittb,2) 727 itc = eqclass(ittc,3) 728 else if (mclass .eq. 4) then 729 ita = eqclass(itta,4) 730 itb = eqclass(ittb,2) 731 itc = eqclass(ittc,4) 732 else if (mclass .eq. 5) then 733 ita = eqclass(itta,5) 734 itb = eqclass(ittb,2) 735 itc = eqclass(ittc,5) 736 end if 737 if (mclass .gt. 5) then 738 goto 20 739 else 740 if (at .eq. 0) then 741 ak(i) = mmff_ka(ita,itb,itc) 742 anat(i) = mmff_ang0(ita,itb,itc) 743 else if (at .eq. 1) then 744 ak(i) = mmff_ka1(ita,itb,itc) 745 anat(i) = mmff_ang1(ita,itb,itc) 746 else if (at .eq. 2) then 747 ak(i) = mmff_ka2(ita,itb,itc) 748 anat(i) = mmff_ang2(ita,itb,itc) 749 else if (at .eq. 3) then 750 ak(i) = mmff_ka3(ita,itb,itc) 751 anat(i) = mmff_ang3(ita,itb,itc) 752 else if (at .eq. 4) then 753 ak(i) = mmff_ka4(ita,itb,itc) 754 anat(i) = mmff_ang4(ita,itb,itc) 755 else if (at .eq. 5) then 756 ak(i) = mmff_ka5(ita,itb,itc) 757 anat(i) = mmff_ang5(ita,itb,itc) 758 else if (at .eq. 6) then 759 ak(i) = mmff_ka6(ita,itb,itc) 760 anat(i) = mmff_ang6(ita,itb,itc) 761 else if (at .eq. 7) then 762 ak(i) = mmff_ka7(ita,itb,itc) 763 anat(i) = mmff_ang7(ita,itb,itc) 764 else if (at .eq. 8) then 765 ak(i) = mmff_ka8(ita,itb,itc) 766 anat(i) = mmff_ang8(ita,itb,itc) 767 end if 768c 769c use empirical rule to calculate the force constant 770c 771 if (mclass .eq. 5) then 772 if (z2(ina) .eq. 1000.0d0) goto 20 773 if (z2(inb) .eq. 1000.0d0) goto 20 774 if (z2(inc) .eq. 1000.0d0) goto 20 775 if (c(ina) .eq. 1000.0d0) goto 20 776 if (c(inb) .eq. 1000.0d0) goto 20 777 if (c(inc) .eq. 1000.0d0) goto 20 778 do k = 1, nbond 779 if ((min(ia,ib).eq.ibnd(1,k)) .and. 780 & (max(ia,ib).eq.ibnd(2,k))) then 781 bnd_ab = k 782 end if 783 if ((min(ic,ib).eq.ibnd(1,k)) .and. 784 & (max(ic,ib).eq.ibnd(2,k))) then 785 bnd_bc = k 786 end if 787 end do 788 d = (bl(bnd_ab)-bl(bnd_bc))**2 789 & / (bl(bnd_ab)+bl(bnd_bc))**2 790 beta = 1.0d0 791 if (ring4) beta = 0.85d0 792 if (ring3) beta = 0.05d0 793 ak(i) = beta*1.75d0*z2(ina)*z2(inc)*c(inb) 794 & / ((0.01745329252d0*anat(i))**2 795 & *(bl(bnd_ab)+bl(bnd_bc))*exp(2.0d0*d)) 796 end if 797 done = .true. 798 if (ak(i) .eq. 1000.0d0) done = .false. 799 if (anat(i) .eq. 1000.0d0) done = .false. 800 if (.not. done) goto 10 801 goto 20 802 end if 803c 804c use empirical rule for ideal angle and force constant 805c 806 20 continue 807 minat = min(ina,inb,inc) 808 if (minat .eq. 0) done = .true. 809 if (.not. done) then 810 if (use_angle) then 811 anat(i) = 120.0d0 812 if (crd(itb) .eq. 4) anat(i) = 109.45d0 813 if (crd(itb) .eq. 2) then 814 if (inb .eq. 8) then 815 anat(i) = 105.0d0 816 else if (inb .gt. 10) then 817 anat(i) = 95.0d0 818 else if (lin(itb) .eq. 1) then 819 anat(i) = 180.0d0 820 end if 821 end if 822 if (crd(itb).eq.3 .and. val(itb).eq.3 823 & .and. mltb(itb).eq.0) then 824 if (inb .eq. 7) then 825 anat(i) = 107.0d0 826 else 827 anat(i) = 92.0d0 828 end if 829 end if 830 if (ring3) anat(i) = 60.0d0 831 if (ring4) anat(i) = 90.0d0 832 do k = 1, nbond 833 if ((min(ia,ib).eq.ibnd(1,k)) .and. 834 & (max(ia,ib).eq.ibnd(2,k))) then 835 bnd_ab = k 836 end if 837 if ((min(ic,ib).eq.ibnd(1,k)) .and. 838 & (max(ic,ib).eq.ibnd(2,k))) then 839 bnd_bc = k 840 end if 841 end do 842 d = (bl(bnd_ab)-bl(bnd_bc))**2 843 & / (bl(bnd_ab)+bl(bnd_bc))**2 844 beta = 1.0d0 845 if (ring4) beta = 0.85d0 846 if (ring3) beta = 0.05d0 847 ak(i) = beta*1.75d0*z2(ina)*z2(inc)*c(inb) 848 & / ((0.01745329252d0*anat(i))**2 849 & *(bl(bnd_ab)+bl(bnd_bc))*exp(2.0d0*d)) 850 end if 851 end if 852 angtyp(i) = 'HARMONIC' 853 if (anat(i) .eq. 180.0d0) angtyp(i) = 'LINEAR' 854 end do 855c 856c turn off the angle bending potential if it is not used 857c 858 if (nangle .eq. 0) use_angle = .false. 859 return 860 end 861