1c 2c 3c ############################################################## 4c ## COPYRIGHT (C) 1998 by Rohit Pappu & Jay William Ponder ## 5c ## All Rights Reserved ## 6c ############################################################## 7c 8c ################################################################# 9c ## ## 10c ## program scan -- maps minima on potential energy surface ## 11c ## ## 12c ################################################################# 13c 14c 15c "scan" attempts to find all the local minima on a potential 16c energy surface via an iterative series of local searches along 17c normal mode directions 18c 19c literature reference: 20c 21c I. Kolossvary and W. C. Guida, "Low-Mode Conformational Search 22c Elucidated: Application to C39H80 and Flexible Docking of 23c 9-Deazaguanine Inhibitors into PNP, Journal of Computational 24c Chemistry, 20, 1671-1684 (1999) 25c 26c 27 program scan 28 use files 29 use inform 30 use iounit 31 use omega 32 use output 33 implicit none 34 integer maxmap 35 parameter (maxmap=100000) 36 integer i,ixyz 37 integer lext,freeunit 38 integer nmap,niter 39 integer nvec,neigen 40 real*8 minimum,grdmin,range 41 real*8 emap(maxmap) 42 logical exist 43 character*7 ext 44 character*240 xyzfile 45 character*240 string 46c 47c 48c set up the structure and mechanics calculation 49c 50 call initial 51 call getxyz 52 call mechanic 53c 54c initialize the number of minima and coordinate type 55c 56 nmap = 0 57 coordtype = 'CARTESIAN' 58c 59c get the rotatable bonds for torsional local search 60c 61 call makeint (0) 62 call initrot 63 call active 64c 65c get the number of eigenvectors to use for the local search 66c 67 neigen = -1 68 call nextarg (string,exist) 69 if (exist) read (string,*,err=10,end=10) neigen 70 10 continue 71 nvec = min(nomega,5) 72 if (neigen .le. 0) then 73 write (iout,20) nvec 74 20 format(/,' Enter the Number of Eigenvectors for Local', 75 & ' Search [',i1,'] : ',$) 76 read (input,30) neigen 77 30 format (i10) 78 if (neigen .le. 0) neigen = nvec 79 end if 80 neigen = min(neigen,nvec) 81c 82c get the energy threshold criterion for map membership 83c 84 range = -1.0d0 85 call nextarg (string,exist) 86 if (exist) read (string,*,err=40,end=40) range 87 40 continue 88 if (range .le. 0.0d0) then 89 write (iout,50) 90 50 format (/,' Enter the Energy Threshold for Local Minima', 91 & ' [100.0] : ',$) 92 read (input,60) range 93 60 format (f20.0) 94 end if 95 if (range .le. 0.0d0) range = 100.0d0 96c 97c get the termination criterion as RMS gradient per atom 98c 99 grdmin = -1.0d0 100 call nextarg (string,exist) 101 if (exist) read (string,*,err=70,end=70) grdmin 102 70 continue 103 if (grdmin .le. 0.0d0) then 104 write (iout,80) 105 80 format (/,' Enter RMS Gradient per Atom Criterion', 106 & ' [0.0001] : ',$) 107 read (input,90) grdmin 108 90 format (f20.0) 109 end if 110 if (grdmin .le. 0.0d0) grdmin = 0.0001d0 111c 112c set the energy output precision via convergence criterion 113c 114 if (grdmin .le. 0.000001d0) digits = 6 115 if (grdmin .le. 0.00000001d0) digits = 8 116c 117c create and open an output file if using archive mode 118c 119 if (archive) then 120 ixyz = freeunit () 121 xyzfile = filename(1:leng) 122 call suffix (xyzfile,'arc','new') 123 open (unit=ixyz,file=xyzfile,status='new') 124 close (unit=ixyz) 125 end if 126c 127c find the first map point from the input structure 128c 129 write (iout,100) 130 100 format (/,' Generating Seed Point for Potential Energy', 131 & ' Surface Scan',/) 132 call localmin (minimum,grdmin) 133 call mapcheck (nmap,emap,range,minimum,grdmin) 134c 135c use normal mode local search to explore adjacent minima 136c 137 niter = 0 138 do while (niter .lt. nmap) 139 niter = niter + 1 140 write (iout,110) niter 141 110 format (/,' Normal Mode Local Search',7x,'Minimum',i7,/) 142 ixyz = freeunit () 143 if (archive) then 144 xyzfile = filename(1:leng) 145 call suffix (xyzfile,'arc','old') 146 open (unit=ixyz,file=xyzfile,status='old') 147 do i = 1, niter-1 148 call readxyz (ixyz) 149 end do 150 else 151 lext = 3 152 call numeral (niter,ext,lext) 153 xyzfile = filename(1:leng)//'.'//ext(1:lext) 154 call version (xyzfile,'old') 155 open (unit=ixyz,file=xyzfile,status='old') 156 end if 157 call readxyz (ixyz) 158 close (unit=ixyz) 159 call modesrch (nmap,emap,range,neigen,grdmin) 160 end do 161c 162c perform any final tasks before program exit 163c 164 call final 165 end 166c 167c 168c ############################################################### 169c ## ## 170c ## subroutine mapcheck -- addition to local minimum list ## 171c ## ## 172c ############################################################### 173c 174c 175c "mapcheck" checks the current minimum energy structure 176c for possible addition to the master list of local minima 177c 178c 179 subroutine mapcheck (nmap,emap,range,minimum,grdmin) 180 use files 181 use inform 182 use iounit 183 use output 184 implicit none 185 integer i,ixyz,lext 186 integer nmap,freeunit 187 real*8 minimum,grdmin 188 real*8 delta,eps,range 189 real*8 emap(*) 190 logical unique,exist 191 character*7 ext 192 character*240 xyzfile 193c 194c 195c check to see if the current minimum was previously found 196c 197 eps = grdmin 198 unique = .true. 199 do i = 1, nmap 200 delta = minimum - emap(i) 201 if (abs(delta) .lt. eps) unique = .false. 202 if (delta .gt. range) unique = .false. 203 end do 204c 205c add minimum to master list if it was not previously known 206c 207 if (unique) then 208 nmap = nmap + 1 209 emap(nmap) = minimum 210 if (digits .ge. 8) then 211 write (iout,10) nmap,minimum 212 10 format (/,4x,'Potential Surface Map',7x,'Minimum', 213 & i7,6x,f20.8,/) 214 else if (digits .ge. 6) then 215 write (iout,20) nmap,minimum 216 20 format (/,4x,'Potential Surface Map',7x,'Minimum', 217 & i7,6x,f18.6,/) 218 else 219 write (iout,30) nmap,minimum 220 30 format (/,4x,'Potential Surface Map',7x,'Minimum', 221 & i7,6x,f16.4,/) 222 end if 223c 224c write the coordinates of the new minimum to a file 225c 226 ixyz = freeunit () 227 if (archive) then 228 xyzfile = filename(1:leng) 229 call suffix (xyzfile,'arc','old') 230 inquire (file=xyzfile,exist=exist) 231 if (exist) then 232 call openend (ixyz,xyzfile) 233 else 234 open (unit=ixyz,file=xyzfile,status='new') 235 end if 236 else 237 lext = 3 238 call numeral (nmap,ext,lext) 239 xyzfile = filename(1:leng)//'.'//ext(1:lext) 240 call version (xyzfile,'new') 241 open (unit=ixyz,file=xyzfile,status='new') 242 end if 243 call prtxyz (ixyz) 244 close (unit=ixyz) 245 end if 246 return 247 end 248c 249c 250c ############################################################### 251c ## ## 252c ## function scan1 -- energy and gradient values for scan ## 253c ## ## 254c ############################################################### 255c 256c 257c "scan1" is a service routine that computes the energy and 258c gradient during exploration of a potential energy surface 259c via iterative local search 260c 261c 262 function scan1 (xx,g) 263 use atoms 264 implicit none 265 integer i,nvar 266 real*8 scan1,e 267 real*8 xx(*) 268 real*8 g(*) 269 real*8, allocatable :: derivs(:,:) 270c 271c 272c convert optimization parameters to atomic coordinates 273c 274 nvar = 0 275 do i = 1, n 276 nvar = nvar + 1 277 x(i) = xx(nvar) 278 nvar = nvar + 1 279 y(i) = xx(nvar) 280 nvar = nvar + 1 281 z(i) = xx(nvar) 282 end do 283c 284c perform dynamic allocation of some local arrays 285c 286 allocate (derivs(3,n)) 287c 288c compute and store the energy and gradient 289c 290 call gradient (e,derivs) 291 scan1 = e 292c 293c convert gradient components to optimization parameters 294c 295 nvar = 0 296 do i = 1, n 297 nvar = nvar + 1 298 g(nvar) = derivs(1,i) 299 nvar = nvar + 1 300 g(nvar) = derivs(2,i) 301 nvar = nvar + 1 302 g(nvar) = derivs(3,i) 303 end do 304c 305c perform deallocation of some local arrays 306c 307 deallocate (derivs) 308 return 309 end 310c 311c 312c ############################################################ 313c ## ## 314c ## subroutine scan2 -- Hessian matrix values for scan ## 315c ## ## 316c ############################################################ 317c 318c 319c "scan2" is a service routine that computes the sparse matrix 320c Hessian elements during exploration of a potential energy 321c surface via iterative local search 322c 323c 324 subroutine scan2 (mode,xx,h,hinit,hstop,hindex,hdiag) 325 use atoms 326 implicit none 327 integer i,nvar 328 integer hinit(*) 329 integer hstop(*) 330 integer hindex(*) 331 real*8 xx(*) 332 real*8 hdiag(*) 333 real*8 h(*) 334 character*4 mode 335c 336c 337c translate optimization parameters to atomic coordinates 338c 339 if (mode .eq. 'NONE') return 340 nvar = 0 341 do i = 1, n 342 nvar = nvar + 1 343 x(i) = xx(nvar) 344 nvar = nvar + 1 345 y(i) = xx(nvar) 346 nvar = nvar + 1 347 z(i) = xx(nvar) 348 end do 349c 350c compute and store the Hessian elements 351c 352 call hessian (h,hinit,hstop,hindex,hdiag) 353 return 354 end 355c 356c 357c ######################################################### 358c ## ## 359c ## subroutine modesrch -- normal mode local search ## 360c ## ## 361c ######################################################### 362c 363c 364 subroutine modesrch (nmap,emap,range,neigen,grdmin) 365 use iounit 366 use omega 367 implicit none 368 integer i,k,nsearch 369 integer nmap,neigen 370 real*8 minimum,grdmin,range 371 real*8 emap(*) 372 real*8, allocatable :: step(:) 373 real*8, allocatable :: eigen(:) 374 real*8, allocatable :: vects(:,:) 375c 376c 377c store the current coordinates as the reference set 378c 379 call makeref (1) 380c 381c perform dynamic allocation of some local arrays 382c 383 allocate (step(nomega)) 384 allocate (eigen(nomega)) 385 allocate (vects(nomega,nomega)) 386c 387c convert to internal coordinates and find torsional modes 388c 389 call makeint (0) 390 call eigenrot (eigen,vects) 391c 392c search both directions along each torsional eigenvector 393c 394 nsearch = 0 395 do i = 1, neigen 396 do k = 1, nomega 397 step(k) = vects(k,nomega-i+1) 398 end do 399 nsearch = nsearch + 1 400 call climber (nsearch,minimum,step,grdmin) 401 call mapcheck (nmap,emap,range,minimum,grdmin) 402 do k = 1, nomega 403 step(k) = -vects(k,nomega-i+1) 404 end do 405 nsearch = nsearch + 1 406 call climber (nsearch,minimum,step,grdmin) 407 call mapcheck (nmap,emap,range,minimum,grdmin) 408 end do 409c 410c perform deallocation of some local arrays 411c 412 deallocate (step) 413 deallocate (eigen) 414 deallocate (vects) 415 return 416 end 417c 418c 419c ############################################################### 420c ## ## 421c ## subroutine eigenrot -- torsional Hessian eigenvectors ## 422c ## ## 423c ############################################################### 424c 425c 426 subroutine eigenrot (eigen,vects) 427 use atoms 428 use omega 429 implicit none 430 integer i,j,ihess 431 real*8 vnorm 432 real*8 eigen(*) 433 real*8, allocatable :: matrix(:) 434 real*8 vects(nomega,*) 435 real*8, allocatable :: hrot(:,:) 436c 437c 438c perform dynamic allocation of some local arrays 439c 440 allocate (matrix(nomega*(nomega+1)/2)) 441 allocate (hrot(nomega,nomega)) 442c 443c compute the Hessian in torsional space 444c 445 call hessrot ('FULL',hrot) 446c 447c place Hessian elements into triangular form 448c 449 ihess = 0 450 do i = 1, nomega 451 do j = i, nomega 452 ihess = ihess + 1 453 matrix(ihess) = hrot(i,j) 454 end do 455 end do 456c 457c diagonalize the Hessian to obtain eigenvalues 458c 459 call diagq (nomega,nomega,matrix,eigen,vects) 460c 461c perform deallocation of some local arrays 462c 463 deallocate (matrix) 464 deallocate (hrot) 465c 466c normalize the torsional Hessian eigenvectors 467c 468 do i = 1, nomega 469 vnorm = 0.0d0 470 do j = 1, nomega 471 vnorm = vnorm + vects(j,i)**2 472 end do 473 vnorm = sqrt(vnorm) 474 do j = 1, nomega 475 vects(j,i) = vects(j,i) / vnorm 476 end do 477 end do 478 return 479 end 480c 481c 482c ############################################################### 483c ## ## 484c ## subroutine climber -- explore single search direction ## 485c ## ## 486c ############################################################### 487c 488c 489 subroutine climber (nsearch,minimum,step,grdmin) 490 use inform 491 use iounit 492 use math 493 use omega 494 use potent 495 use zcoord 496 implicit none 497 integer maxstep 498 parameter (maxstep=500) 499 integer i,kstep 500 integer nstep,nsearch 501 real*8 minimum,grdmin 502 real*8 big,energy,size 503 real*8 estep(0:maxstep) 504 real*8 step(*) 505 logical done 506 logical oldpolar 507c 508c 509c convert current reference coordinates to a Z-matrix 510c 511 call getref (1) 512 call makeint (0) 513c 514c set the maximum number of steps and the step size 515c 516 done = .false. 517 big = 100000.0d0 518 minimum = big 519 kstep = 0 520 nstep = 65 521 size = 0.1d0 * radian 522 do i = 1, nomega 523 step(i) = size * step(i) 524 end do 525c 526c scan the search direction for a minimization candidate 527c 528 do while (.not. done) 529 if (kstep .ne. 0) then 530 do i = 1, nomega 531 ztors(zline(i)) = ztors(zline(i)) + step(i) 532 end do 533 end if 534 call makexyz 535 oldpolar = use_polar 536 use_polar = .false. 537 estep(kstep) = energy () 538 use_polar = oldpolar 539 if (kstep .ge. 2) then 540 if (estep(kstep) .lt. estep(kstep-2) .and. 541 & estep(kstep-1) .lt. estep(kstep-2)) then 542 done = .true. 543 do i = 1, nomega 544 ztors(zline(i)) = ztors(zline(i)) - step(i) 545 end do 546 call makexyz 547 call localmin (minimum,grdmin) 548 if (minimum .le. -big) then 549 minimum = big 550 write (iout,10) nsearch 551 10 format (4x,'Search Direction',i4,38x,'<<<<<<') 552 else if (minimum .ge. big) then 553 minimum = big 554 write (iout,20) nsearch 555 20 format (4x,'Search Direction',i4,38x,'>>>>>>') 556 else 557 if (digits .ge. 8) then 558 write (iout,30) nsearch,kstep-1,minimum 559 30 format (4x,'Search Direction',i4,11x,'Step', 560 & i7,6x,f20.8) 561 else if (digits .ge. 6) then 562 write (iout,40) nsearch,kstep-1,minimum 563 40 format (4x,'Search Direction',i4,11x,'Step', 564 & i7,6x,f18.6) 565 else 566 write (iout,50) nsearch,kstep-1,minimum 567 50 format (4x,'Search Direction',i4,11x,'Step', 568 & i7,6x,f16.4) 569 end if 570 end if 571 end if 572 end if 573 if (kstep.ge.nstep .and. .not.done) then 574 done = .true. 575 write (iout,60) nsearch 576 60 format (4x,'Search Direction',i4,38x,'------') 577 end if 578 kstep = kstep + 1 579 end do 580 return 581 end 582c 583c 584c ################################################################ 585c ## ## 586c ## subroutine localmin -- optimize local search candidate ## 587c ## ## 588c ################################################################ 589c 590c 591c "localmin" is used during normal mode local search to 592c perform a Cartesian coordinate energy minimization 593c 594c 595 subroutine localmin (minimum,grdmin) 596 use atoms 597 use inform 598 use minima 599 use output 600 use potent 601 use scales 602 implicit none 603 integer i,j,nvar 604 real*8 minimum,scan1 605 real*8 grdmin,oldgrd 606 real*8 gnorm,grms,big 607 real*8, allocatable :: xx(:) 608 real*8, allocatable :: derivs(:,:) 609 logical oldverb,oldpolar 610 character*6 mode,method 611 external scan1,scan2 612 external optsave 613c 614c 615c initialize optimization output and maximum energy 616c 617 iwrite = 0 618 iprint = 0 619 big = 100000.0d0 620c 621c perform dynamic allocation of some local arrays 622c 623 allocate (xx(3*n)) 624c 625c convert atomic coordinates to optimization parameters 626c 627 nvar = 0 628 do i = 1, n 629 nvar = nvar + 1 630 xx(nvar) = x(i) 631 nvar = nvar + 1 632 xx(nvar) = y(i) 633 nvar = nvar + 1 634 xx(nvar) = z(i) 635 end do 636c 637c perform dynamic allocation of some global arrays 638c 639 if (.not. set_scale) then 640 if (.not. allocated(scale)) allocate (scale(nvar)) 641c 642c set scaling parameters to unity due to mixed optimization 643c 644 set_scale = .true. 645 do i = 1, nvar 646 scale(i) = 1.0d0 647 end do 648 end if 649c 650c adjust polarization and set initial optimization values 651c 652 oldverb = verbose 653 oldpolar = use_polar 654 oldgrd = grdmin 655 verbose = .false. 656 use_polar = .false. 657 grdmin = 3.0 658c 659c initial optimizaton to get close to approximate minimum 660c 661 call lbfgs (nvar,xx,minimum,grdmin,scan1,optsave) 662c 663c secondary optimization to reach the exact local minimum 664c 665 use_polar = oldpolar 666 grdmin = oldgrd 667 mode = 'AUTO' 668 method = 'AUTO' 669 call tncg (mode,method,nvar,xx,minimum, 670 & grdmin,scan1,scan2,optsave) 671 verbose = oldverb 672c 673c convert optimization parameters to atomic coordinates 674c 675 nvar = 0 676 do i = 1, n 677 nvar = nvar + 1 678 x(i) = xx(nvar) 679 nvar = nvar + 1 680 y(i) = xx(nvar) 681 nvar = nvar + 1 682 z(i) = xx(nvar) 683 end do 684c 685c perform deallocation of some local arrays 686c 687 deallocate (xx) 688c 689c perform dynamic allocation of some local arrays 690c 691 allocate (derivs(3,n)) 692c 693c independently check the gradient convergence criterion 694c 695 call gradient (minimum,derivs) 696 gnorm = 0.0d0 697 do i = 1, n 698 do j = 1, 3 699 gnorm = gnorm + derivs(j,i)**2 700 end do 701 end do 702 gnorm = sqrt(gnorm) 703 grms = gnorm / sqrt(dble(n)) 704 if (grms .gt. grdmin) minimum = big 705c 706c perform deallocation of some local arrays 707c 708 deallocate (derivs) 709 return 710 end 711