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