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