1c 2c 3c ################################################################ 4c ## COPYRIGHT (C) 2001 by ## 5c ## Michael Schnieders, Alan Grossfield & Jay William Ponder ## 6c ## All Rights Reserved ## 7c ################################################################ 8c 9c ################################################################# 10c ## ## 11c ## program monte -- Monte Carlo-Minimization search method ## 12c ## ## 13c ################################################################# 14c 15c 16c "monte" performs a Monte Carlo-Minimization conformational 17c search using Cartesian single atom or torsional move sets 18c 19c literature references: 20c 21c Z. Li and H. A. Scheraga, "Monte Carlo-Minimization Approach 22c to the Multiple-Minima Problem in Protein Folding", Proc. Natl. 23c Acad. Sci. USA, 84, 6611-6615 (1987) 24c 25c D. J. Wales, "Energy Landscapes with Applications to Clusters, 26c Biomolecules and Glasses", Cambridge University Press, 2003, 27c Section 6.7.4 28c 29c 30 program monte 31 use atoms 32 use files 33 use inform 34 use iounit 35 use omega 36 use output 37 use units 38 use usage 39 use zcoord 40 implicit none 41 integer i,k,m,next 42 integer keep,nbig 43 integer nmap,lext 44 integer istep,nstep 45 integer ixyz,freeunit 46 real*8 global,ratio 47 real*8 big,eps,size 48 real*8 grdmin,temper 49 real*8 minimum,pminimum 50 real*8 tsize,factor 51 real*8 beta,boltz 52 real*8 random,trial 53 real*8 converge,delta 54 real*8 efficient 55 real*8 vector(3) 56 real*8, allocatable :: xg(:) 57 real*8, allocatable :: yg(:) 58 real*8, allocatable :: zg(:) 59 real*8, allocatable :: xi(:) 60 real*8, allocatable :: yi(:) 61 real*8, allocatable :: zi(:) 62 real*8, allocatable :: xp(:) 63 real*8, allocatable :: yp(:) 64 real*8, allocatable :: zp(:) 65 logical exist,reset,done 66 logical torsmove 67 character*1 answer 68 character*6 status 69 character*7 ext 70 character*240 xyzfile 71 character*240 record 72 character*240 string 73 external random 74c 75c 76c set up the structure and mechanics calculation 77c 78 call initial 79 call getxyz 80 call mechanic 81c 82c initialize values of some counters and parameters 83c 84 istep = 0 85 keep = 0 86 nbig = 0 87 nmap = 0 88 delta = 0.00001d0 89 eps = 0.0001d0 90 big = 100000.0d0 91 reset = .false. 92c 93c get the desired number of Monte Carlo steps 94c 95 nstep = -1 96 call nextarg (string,exist) 97 if (exist) read (string,*,err=10,end=10) nstep 98 10 continue 99 if (nstep .le. 0) then 100 write (iout,20) 101 20 format (/,' Maximum Number of Monte Carlo Steps [1000] : ', $) 102 read (input,30) nstep 103 30 format (i10) 104 if (nstep .le. 0) nstep = 1000 105 end if 106c 107c get the search efficiency criterion for convergence 108c 109 converge = -1.0d0 110 call nextarg (string,exist) 111 if (exist) read (string,*,err=40,end=40) converge 112 40 continue 113 if (converge .lt. 0.0d0) then 114 write (iout,50) 115 50 format (/,' Enter Search Efficiency Termination Criterion', 116 & ' [0.01] : ', $) 117 read (input,60) string 118 60 format (a240) 119 read (string,*,err=70,end=70) converge 120 70 continue 121 if (converge .lt. 0.0d0) converge = 0.01 122 end if 123 converge = converge + delta 124c 125c choose either the torsional or single atom move set 126c 127 torsmove = .false. 128 call nextarg (answer, exist) 129 if (.not. exist) then 130 write (iout,80) 131 80 format (/,' Use [C]artesian or [T]orsional Moves [C] : ',$) 132 read (input,90) record 133 90 format (a240) 134 next = 1 135 call gettext (record,answer,next) 136 end if 137 call upcase (answer) 138 if (answer .eq. 'T') torsmove = .true. 139c 140c for torsional moves, generate the internal coordinates 141c 142 if (torsmove) then 143 call makeint (0) 144 call initrot 145c 146c set all atoms active to simplify torsional calculation 147c 148 nuse = n 149 do i = 1, n 150 use(i) = .true. 151 end do 152 end if 153c 154c get the desired Cartesian or torsional step size 155c 156 size = -1.0d0 157 call nextarg (string,exist) 158 if (exist) read (string,*,err=100,end=100) size 159 100 continue 160 if (size .lt. 0.0d0) then 161 if (torsmove) then 162 write (iout,110) 163 110 format (/,' Enter Maximum Step in Degrees [180.0] : ', $) 164 else 165 write (iout,120) 166 120 format (/,' Enter Maximum Step in Angstroms [3.0] : ', $) 167 end if 168 read (input,130) string 169 130 format (a240) 170 read (string,*,err=140,end=140) size 171 140 continue 172 if (size .lt. 0.0d0) then 173 if (torsmove) then 174 size = 180.0d0 175 else 176 size = 3.0d0 177 end if 178 end if 179 if (torsmove) size = min(size,180.0d0) 180 end if 181c 182c get the gradient convergence for local minimizations 183c 184 grdmin = -1.0d0 185 call nextarg (string,exist) 186 if (exist) read (string,*,err=150,end=150) grdmin 187 150 continue 188 if (grdmin .lt. 0.0d0) then 189 write (iout,160) 190 160 format (/,' Enter RMS Gradient Criterion for Minima', 191 & ' [0.01] : ', $) 192 read (input,170) string 193 170 format (a240) 194 read (string,*,err=180,end=180) grdmin 195 180 continue 196 if (grdmin .lt. 0.0d0) grdmin = 0.01 197 end if 198c 199c get the desired temperature for Metropolis criterion 200c 201 temper = -1.0d0 202 call nextarg (string,exist) 203 if (exist) read (string,*,err=190,end=190) temper 204 190 continue 205 if (temper .lt. 0.0d0) then 206 write (iout,200) 207 200 format (/,' Enter the Desired Temperature in Degrees', 208 & ' K [500] : ', $) 209 read (input,210) string 210 210 format (a240) 211 read (string,*,err=220,end=220) temper 212 220 continue 213 if (temper .lt. 0.0d0) temper = 500.0d0 214 end if 215 beta = 1.0d0 / (gasconst*temper) 216c 217c perform dynamic allocation of some local arrays 218c 219 allocate (xg(n)) 220 allocate (yg(n)) 221 allocate (zg(n)) 222 allocate (xi(n)) 223 allocate (yi(n)) 224 allocate (zi(n)) 225 allocate (xp(n)) 226 allocate (yp(n)) 227 allocate (zp(n)) 228c 229c print some information prior to initial iteration 230c 231 write (iout,230) 232 230 format (/,' Monte Carlo Minimization Global Search :') 233 write (iout,240) 234 240 format (/,' MCM Iter Current Global', 235 & ' Efficiency Accept Status',/) 236 flush (iout) 237c 238c create and open an output file if using archive mode 239c 240 if (archive) then 241 ixyz = freeunit () 242 xyzfile = filename(1:leng) 243 call suffix (xyzfile,'arc','new') 244 open (unit=ixyz,file=xyzfile,status='new') 245 close (unit=ixyz) 246 end if 247c 248c store the coordinates, then perform a minimization 249c 250 do i = 1, n 251 xi(i) = x(i) 252 yi(i) = y(i) 253 zi(i) = z(i) 254 end do 255 call mcmstep (minimum,grdmin) 256 pminimum = minimum 257 write (iout,250) 0,minimum 258 250 format (i8,3x,f12.4) 259c 260c save coordinates as the initial global minimum 261c 262 do i = 1, n 263 xg(i) = x(i) 264 yg(i) = y(i) 265 zg(i) = z(i) 266 end do 267 global = minimum 268 nmap = nmap + 1 269 lext = 3 270 call numeral (nmap,ext,lext) 271 ixyz = freeunit () 272 if (archive) then 273 xyzfile = filename(1:leng) 274 call suffix (xyzfile,'arc','old') 275 inquire (file=xyzfile,exist=exist) 276 if (exist) then 277 call openend (ixyz,xyzfile) 278 else 279 open (unit=ixyz,file=xyzfile,status='new') 280 end if 281 else 282 xyzfile = filename(1:leng)//'.'//ext(1:lext) 283 call version (xyzfile,'new') 284 open (unit=ixyz,file=xyzfile,status='new') 285 end if 286 call prtxyz (ixyz) 287 close (unit=ixyz) 288 write (iout,260) nmap,global 289 260 format (/,4x,'Minimum Energy Structure',i7,6x,f16.4,/) 290 call flush (iout) 291c 292c optionally reset coordinates to before the minimization 293c 294 if (reset) then 295 do i = 1, n 296 x(i) = xi(i) 297 y(i) = yi(i) 298 z(i) = zi(i) 299 end do 300 end if 301 if (torsmove) call makeint (2) 302c 303c store the prior coordinates to start each MCM iteration 304c 305 done = .false. 306 dowhile (.not. done) 307 istep = istep + 1 308 do i = 1, n 309 xp(i) = x(i) 310 yp(i) = y(i) 311 zp(i) = z(i) 312 end do 313c 314c generate random angle moves for a few torsions 315c 316 if (torsmove) then 317 m = int(-log(max(random(),0.0001d0))) + 1 318 do i = 1, m 319 k = int(nomega * random()) + 1 320 k = zline(k) 321 tsize = 2.0d0 * size * (random()-0.5d0) 322 ztors(k) = ztors(k) + tsize 323 if (ztors(k) .gt. 180.0d0) then 324 ztors(k) = ztors(k) - 360.0d0 325 else if (ztors(k) .lt. -180.0d0) then 326 ztors(k) = ztors(k) + 360.0d0 327 end if 328 end do 329 call makexyz 330c 331c generate a random Cartesian move for each atom 332c 333 else 334 do i = 1, nuse 335 k = iuse(i) 336 call ranvec (vector) 337 factor = size * random () 338 x(k) = x(k) + factor*vector(1) 339 y(k) = y(k) + factor*vector(2) 340 z(k) = z(k) + factor*vector(3) 341 end do 342 end if 343c 344c store the coordinates, then perform a minimization 345c 346 do i = 1, n 347 xi(i) = x(i) 348 yi(i) = y(i) 349 zi(i) = z(i) 350 end do 351 call mcmstep (minimum,grdmin) 352c 353c test for an unreasonably low energy at the minimum 354c 355 if (minimum .lt. -big) minimum = big 356c 357c step is probably degenerate if energy is identical 358c 359 if (abs(minimum-pminimum) .le. eps) then 360 status = 'Same' 361 pminimum = minimum 362c 363c accept the step if the new minimum has lower energy 364c 365 else if (minimum .le. pminimum) then 366 status = 'Accept' 367 pminimum = minimum 368c 369c if the energy increased, apply the Metropolis criterion 370c 371 else 372 boltz = exp(-beta*(minimum-pminimum)) 373 trial = random () 374c 375c reject the step if the energy increase is too large 376c 377 if (boltz .lt. trial) then 378 status = 'Reject' 379c 380c accept the step if the energy increase is small enough 381c 382 else 383 status = 'Accept' 384 pminimum = minimum 385 end if 386 end if 387c 388c save coordinates with the best energy as global minimum 389c 390 if (minimum .lt. global-eps) then 391 do i = 1, n 392 xg(i) = x(i) 393 yg(i) = y(i) 394 zg(i) = z(i) 395 end do 396 global = minimum 397 nmap = nmap + 1 398 lext = 3 399 call numeral (nmap,ext,lext) 400 ixyz = freeunit () 401 if (archive) then 402 xyzfile = filename(1:leng) 403 call suffix (xyzfile,'arc','old') 404 inquire (file=xyzfile,exist=exist) 405 if (exist) then 406 call openend (ixyz,xyzfile) 407 else 408 open (unit=ixyz,file=xyzfile,status='new') 409 end if 410 else 411 xyzfile = filename(1:leng)//'.'//ext(1:lext) 412 call version (xyzfile,'new') 413 open (unit=ixyz,file=xyzfile,status='new') 414 end if 415 call prtxyz (ixyz) 416 close (unit=ixyz) 417 write (iout,270) nmap,global 418 270 format (/,4x,'Minimum Energy Structure',i7,6x,f16.4,/) 419 flush (iout) 420 end if 421c 422c update the efficiency and Monte Carlo acceptance ratio 423c 424 efficient = dble(nmap) / dble(istep) 425 if (status .eq. 'Accept') keep = keep + 1 426 ratio = dble(keep) / dble(istep) 427c 428c print intermediate results for the current iteration 429c 430 if (istep.ne.1 .and. mod(istep,100).eq.1) then 431 write (iout,280) 432 280 format (/,' MCM Iter Current Global', 433 & ' Efficiency Accept Status',/) 434 end if 435 if (minimum .lt. big) then 436 nbig = 0 437 write (iout,290) istep,minimum,global,efficient, 438 & ratio,status 439 290 format (i8,3x,f12.4,3x,f12.4,3x,f9.4,3x,f9.4,6x,a6) 440 else 441 nbig = nbig + 1 442 write (iout,300) istep,global,efficient,ratio,status 443 300 format (i8,9x,'------',3x,f12.4,3x,f9.4,3x,f9.4,6x,a6) 444 end if 445 flush (iout) 446c 447c restore global minimum after repeated bad iterations 448c 449 if (nbig .ge. 3) then 450 nbig = 0 451 do i = 1, n 452 x(i) = xg(i) 453 y(i) = yg(i) 454 z(i) = zg(i) 455 end do 456c 457c optionally reset coordinates to before the minimization 458c 459 else if (status.eq.'Same' .or. status.eq.'Accept') then 460 if (reset) then 461 do i = 1, n 462 x(i) = xi(i) 463 y(i) = yi(i) 464 z(i) = zi(i) 465 end do 466 end if 467c 468c restore coordinates to those from the previous iteration 469c 470 else if (status .eq. 'Reject') then 471 do i = 1, n 472 x(i) = xp(i) 473 y(i) = yp(i) 474 z(i) = zp(i) 475 end do 476 end if 477c 478c update internal coordinates if using torsional moves 479c 480 if (torsmove) call makeint (2) 481c 482c check criteria based on search efficiency and step number 483c 484 if (efficient .le. converge) then 485 done = .true. 486 write (iout,310) 487 310 format (/,' MONTE -- Termination based on Overall', 488 & ' Search Efficiency') 489 end if 490 if (istep .ge. nstep) then 491 done = .true. 492 write (iout,320) 493 320 format (/,' MONTE -- Termination based on Maximum', 494 & ' MCM Step Limit') 495 end if 496 end do 497c 498c perform deallocation of some local arrays 499c 500 deallocate (xg) 501 deallocate (yg) 502 deallocate (zg) 503 deallocate (xi) 504 deallocate (yi) 505 deallocate (zi) 506 deallocate (xp) 507 deallocate (yp) 508 deallocate (zp) 509c 510c write out the final global minimum energy value 511c 512 if (digits .ge. 8) then 513 write (iout,330) global 514 330 format (/,' Global Minimum Energy Value :',2x,f18.8) 515 else if (digits .ge. 6) then 516 write (iout,340) global 517 340 format (/,' Global Minimum Energy Value :',4x,f16.6) 518 else 519 write (iout,350) global 520 350 format (/,' Global Minimum Energy Value :',6x,f14.4) 521 end if 522c 523c perform any final tasks before program exit 524c 525 call final 526 end 527c 528c 529c ############################################################### 530c ## ## 531c ## function mcmstep -- minimization phase of an MCM step ## 532c ## ## 533c ############################################################### 534c 535c 536c "mcmstep" implements the minimization phase of an MCM step 537c via Cartesian minimization following a Monte Carlo step 538c 539c 540 subroutine mcmstep (minimum,grdmin) 541 use atoms 542 use bound 543 use files 544 use inform 545 use output 546 use potent 547 use usage 548 implicit none 549 integer i,k,nvar 550 real*8 mcm1,minimum,grdmin 551 real*8, allocatable :: xx(:) 552 character*6 mode,method 553 external mcm1,mcm2,optsave 554c 555c 556c prepare for the truncated Newton minimization 557c 558 mode = 'AUTO' 559 method = 'AUTO' 560 verbose = .false. 561 iprint = 0 562 iwrite = 0 563 coordtype = 'CARTESIAN' 564c 565c perform dynamic allocation of some local arrays 566c 567 allocate (xx(3*n)) 568c 569c convert atomic coordinates to optimization parameters 570c 571 nvar = 0 572 do i = 1, nuse 573 k = iuse(i) 574 nvar = nvar + 1 575 xx(nvar) = x(k) 576 nvar = nvar + 1 577 xx(nvar) = y(k) 578 nvar = nvar + 1 579 xx(nvar) = z(k) 580 end do 581c 582c make the call to the optimization routine 583c 584 call tncg (mode,method,nvar,xx,minimum,grdmin, 585 & mcm1,mcm2,optsave) 586c 587c convert optimization parameters to atomic coordinates 588c 589 nvar = 0 590 do i = 1, nuse 591 k = iuse(i) 592 nvar = nvar + 1 593 x(k) = xx(nvar) 594 nvar = nvar + 1 595 y(k) = xx(nvar) 596 nvar = nvar + 1 597 z(k) = xx(nvar) 598 end do 599c 600c maintain any periodic boundary conditions 601c 602 if (use_bounds) call bounds 603c 604c perform deallocation of some local arrays 605c 606 deallocate (xx) 607 return 608 end 609c 610c 611c ############################################################# 612c ## ## 613c ## function mcm1 -- energy and gradient for MCM search ## 614c ## ## 615c ############################################################# 616c 617c 618c "mcm1" is a service routine that computes the energy and 619c gradient for truncated Newton optimization in Cartesian 620c coordinate space 621c 622c 623 function mcm1 (xx,g) 624 use atoms 625 use usage 626 implicit none 627 integer i,k,nvar 628 real*8 mcm1,e 629 real*8 xx(*) 630 real*8 g(*) 631 real*8, allocatable :: derivs(:,:) 632c 633c 634c convert optimization parameters to atomic coordinates 635c 636 nvar = 0 637 do i = 1, nuse 638 k = iuse(i) 639 nvar = nvar + 1 640 x(k) = xx(nvar) 641 nvar = nvar + 1 642 y(k) = xx(nvar) 643 nvar = nvar + 1 644 z(k) = xx(nvar) 645 end do 646c 647c perform dynamic allocation of some local arrays 648c 649 allocate (derivs(3,n)) 650c 651c compute and store the energy and gradient 652c 653 call gradient (e,derivs) 654 mcm1 = e 655c 656c store gradient components to optimization parameters 657c 658 nvar = 0 659 do i = 1, nuse 660 k = iuse(i) 661 nvar = nvar + 1 662 g(nvar) = derivs(1,k) 663 nvar = nvar + 1 664 g(nvar) = derivs(2,k) 665 nvar = nvar + 1 666 g(nvar) = derivs(3,k) 667 end do 668c 669c perform deallocation of some local arrays 670c 671 deallocate (derivs) 672 return 673 end 674c 675c 676c ########################################################## 677c ## ## 678c ## subroutine mcm2 -- Hessian values for MCM search ## 679c ## ## 680c ########################################################## 681c 682c 683c "mcm2" is a service routine that computes the sparse matrix 684c Hessian elements for truncated Newton optimization in Cartesian 685c coordinate space 686c 687c 688 subroutine mcm2 (mode,xx,h,hinit,hstop,hindex,hdiag) 689 use atoms 690 use usage 691 implicit none 692 integer i,j,k,nvar 693 integer hinit(*) 694 integer hstop(*) 695 integer hindex(*) 696 integer, allocatable :: hvar(:) 697 integer, allocatable :: huse(:) 698 real*8 xx(*) 699 real*8 hdiag(*) 700 real*8 h(*) 701 character*4 mode 702c 703c 704c convert optimization parameters to atomic coordinates 705c 706 if (mode .eq. 'NONE') return 707 nvar = 0 708 do i = 1, nuse 709 k = iuse(i) 710 nvar = nvar + 1 711 x(k) = xx(nvar) 712 nvar = nvar + 1 713 y(k) = xx(nvar) 714 nvar = nvar + 1 715 z(k) = xx(nvar) 716 end do 717c 718c compute and store the Hessian elements 719c 720 call hessian (h,hinit,hstop,hindex,hdiag) 721c 722c perform dynamic allocation of some local arrays 723c 724 allocate (hvar(nvar)) 725 allocate (huse(3*n)) 726c 727c transform the sparse Hessian to use only active atoms 728c 729 nvar = 0 730 if (nuse .ne. n) then 731 do i = 1, n 732 k = 3 * (i-1) 733 if (use(i)) then 734 do j = 1, 3 735 nvar = nvar + 1 736 hvar(nvar) = j + k 737 huse(j+k) = nvar 738 end do 739 else 740 do j = 1, 3 741 huse(j+k) = 0 742 end do 743 end if 744 end do 745 do i = 1, nvar 746 k = hvar(i) 747 hinit(i) = hinit(k) 748 hstop(i) = hstop(k) 749 hdiag(i) = hdiag(k) 750 do j = hinit(i), hstop(i) 751 hindex(j) = huse(hindex(j)) 752 end do 753 end do 754 end if 755c 756c convert atomic coordinates to optimization parameters 757c 758 nvar = 0 759 do i = 1, nuse 760 k = iuse(i) 761 nvar = nvar + 1 762 xx(nvar) = x(k) 763 nvar = nvar + 1 764 xx(nvar) = y(k) 765 nvar = nvar + 1 766 xx(nvar) = z(k) 767 end do 768c 769c perform deallocation of some local arrays 770c 771 deallocate (hvar) 772 deallocate (huse) 773 return 774 end 775