source/unres/Tinker-minimize/minimize.f

   1 c
   2 c
   3 c     ###################################################
   4 c     ##  COPYRIGHT (C)  1990  by  Jay William Ponder  ##
   5 c     ##              All Rights Reserved              ##
   6 c     ###################################################
   7 c
   8 c     ##################################################################
   9 c     ##                                                              ##
  10 c     ##  program minimize  --  low storage BFGS Cartesian optimizer  ##
  11 c     ##                                                              ##
  12 c     ##################################################################
  13 c
  14 c
  15 c     "minimize" performs energy minimization in Cartesian coordinate
  16 c     space using a low storage BFGS nonlinear optimization
  17 c
  18 c
  19       program minimize
  20       use atoms
  21       use bound
  22       use files
  23       use freeze
  24       use inform
  25       use iounit
  26       use scales
  27       use usage
  28       implicit none
  29       integer i,j
  30       integer imin,nvar
  31       integer freeunit
  32       real*8 minimum,minimiz1
  33       real*8 grdmin,gnorm,grms
  34       real*8 energy,eps
  35       real*8, allocatable :: xx(:)
  36       real*8, allocatable :: derivs(:,:)
  37       logical exist,analytic
  38       character*240 minfile
  39       character*240 string
  40       external energy
  41       external minimiz1
  42       external optsave
  43 c
  44 c
  45 c     set up the structure and mechanics calculation
  46 c
  47       call initial
  48       call getxyz
  49       call mechanic
  50 c
  51 c     perform the setup functions needed for optimization
  52 c
  53       call optinit
  54 c
  55 c     use either analytical or numerical gradients
  56 c
  57       analytic = .true.
  58       eps = 0.00001d0
  59 c
  60 c     get termination criterion as RMS gradient per atom
  61 c
  62       grdmin = -1.0d0
  63       call nextarg (string,exist)
  64       if (exist)  read (string,*,err=10,end=10)  grdmin
  65    10 continue
  66       if (grdmin .le. 0.0d0) then
  67          write (iout,20)
  68    20    format (/,' Enter RMS Gradient per Atom Criterion',
  69      &              ' [0.01] :  ',$)
  70          read (input,30)  grdmin
  71    30    format (f20.0)
  72       end if
  73       if (grdmin .le. 0.0d0)  grdmin = 0.01d0
  74 c
  75 c     write out a copy of coordinates for later update
  76 c
  77       imin = freeunit ()
  78       minfile = filename(1:leng)//'.xyz'
  79       call version (minfile,'new')
  80       open (unit=imin,file=minfile,status='new')
  81       call prtxyz (imin)
  82       close (unit=imin)
  83       outfile = minfile
  84 c
  85 c     perform dynamic allocation of some global arrays
  86 c
  87       if (.not. allocated(scale))  allocate (scale(3*n))
  88 c
  89 c     set scaling parameter for function and derivative values;
  90 c     use square root of median eigenvalue of typical Hessian
  91 c
  92       set_scale = .true.
  93       nvar = 0
  94       do i = 1, n
  95          if (use(i)) then
  96             do j = 1, 3
  97                nvar = nvar + 1
  98                scale(nvar) = 12.0d0
  99             end do
 100          end if
 101       end do
 102 c
 103 c     perform dynamic allocation of some local arrays
 104 c
 105       allocate (xx(nvar))
 106       allocate (derivs(3,n))
 107 c
 108 c     convert atomic coordinates to optimization parameters
 109 c
 110       nvar = 0
 111       do i = 1, n
 112          if (use(i)) then
 113             nvar = nvar + 1
 114             xx(nvar) = x(i) * scale(nvar)
 115             nvar = nvar + 1
 116             xx(nvar) = y(i) * scale(nvar)
 117             nvar = nvar + 1
 118             xx(nvar) = z(i) * scale(nvar)
 119          end if
 120       end do
 121 c
 122 c     make the call to the optimization routine
 123 c
 124       call lbfgs (nvar,xx,minimum,grdmin,minimiz1,optsave)
 125 c
 126 c     convert optimization parameters to atomic coordinates
 127 c
 128       nvar = 0
 129       do i = 1, n
 130          if (use(i)) then
 131             nvar = nvar + 1
 132             x(i) = xx(nvar) / scale(nvar)
 133             nvar = nvar + 1
 134             y(i) = xx(nvar) / scale(nvar)
 135             nvar = nvar + 1
 136             z(i) = xx(nvar) / scale(nvar)
 137          end if
 138       end do
 139 c
 140 c     compute the final function and RMS gradient values
 141 c
 142       if (analytic) then
 143          call gradient (minimum,derivs)
 144       else
 145          minimum = energy ()
 146          call numgrad (energy,derivs,eps)
 147       end if
 148       if (use_rattle)  call shakef (derivs)
 149       gnorm = 0.0d0
 150       do i = 1, n
 151          if (use(i)) then
 152             do j = 1, 3
 153                gnorm = gnorm + derivs(j,i)**2
 154             end do
 155          end if
 156       end do
 157       gnorm = sqrt(gnorm)
 158       grms = gnorm / sqrt(dble(nvar/3))
 159 c
 160 c     perform deallocation of some local arrays
 161 c
 162       deallocate (xx)
 163       deallocate (derivs)
 164 c
 165 c     write out the final function and gradient values
 166 c
 167       if (digits .ge. 8) then
 168          if (grms .gt. 1.0d-8) then
 169             write (iout,40)  minimum,grms,gnorm
 170    40       format (/,' Final Function Value :',2x,f20.8,
 171      &              /,' Final RMS Gradient :',4x,f20.8,
 172      &              /,' Final Gradient Norm :',3x,f20.8)
 173          else
 174             write (iout,50)  minimum,grms,gnorm
 175    50       format (/,' Final Function Value :',2x,f20.8,
 176      &              /,' Final RMS Gradient :',4x,d20.8,
 177      &              /,' Final Gradient Norm :',3x,d20.8)
 178          end if
 179       else if (digits .ge. 6) then
 180          if (grms .gt. 1.0d-6) then
 181             write (iout,60)  minimum,grms,gnorm
 182    60       format (/,' Final Function Value :',2x,f18.6,
 183      &              /,' Final RMS Gradient :',4x,f18.6,
 184      &              /,' Final Gradient Norm :',3x,f18.6)
 185          else
 186             write (iout,70)  minimum,grms,gnorm
 187    70       format (/,' Final Function Value :',2x,f18.6,
 188      &              /,' Final RMS Gradient :',4x,d18.6,
 189      &              /,' Final Gradient Norm :',3x,d18.6)
 190          end if
 191       else
 192          if (grms .gt. 1.0d-4) then
 193             write (iout,80)  minimum,grms,gnorm
 194    80       format (/,' Final Function Value :',2x,f16.4,
 195      &              /,' Final RMS Gradient :',4x,f16.4,
 196      &              /,' Final Gradient Norm :',3x,f16.4)
 197          else
 198             write (iout,90)  minimum,grms,gnorm
 199    90       format (/,' Final Function Value :',2x,f16.4,
 200      &              /,' Final RMS Gradient :',4x,d16.4,
 201      &              /,' Final Gradient Norm :',3x,d16.4)
 202          end if
 203       end if
 204 c
 205 c     write the final coordinates into a file
 206 c
 207       if (use_bounds)  call bounds
 208       imin = freeunit ()
 209       open (unit=imin,file=minfile,status='old')
 210       rewind (unit=imin)
 211       call prtxyz (imin)
 212       close (unit=imin)
 213 c
 214 c     perform any final tasks before program exit
 215 c
 216       call final
 217       end
 218 c
 219 c
 220 c     ###############################################################
 221 c     ##                                                           ##
 222 c     ##  function minimiz1  --  energy and gradient for minimize  ##
 223 c     ##                                                           ##
 224 c     ###############################################################
 225 c
 226 c
 227 c     "minimiz1" is a service routine that computes the energy and
 228 c     gradient for a low storage BFGS optimization in Cartesian
 229 c     coordinate space
 230 c
 231 c
 232       function minimiz1 (xx,g)
 233       use atoms
 234       use freeze
 235       use scales
 236       use usage
 237       implicit none
 238       integer i,nvar
 239       real*8 minimiz1,e
 240       real*8 energy,eps
 241       real*8 xx(*)
 242       real*8 g(*)
 243       real*8, allocatable :: derivs(:,:)
 244       logical analytic
 245       external energy
 246 c
 247 c
 248 c     use either analytical or numerical gradients
 249 c
 250       analytic = .true.
 251       eps = 0.00001d0
 252 c
 253 c     convert optimization parameters to atomic coordinates
 254 c
 255       nvar = 0
 256       do i = 1, n
 257          if (use(i)) then
 258             nvar = nvar + 1
 259             x(i) = xx(nvar) / scale(nvar)
 260             nvar = nvar + 1
 261             y(i) = xx(nvar) / scale(nvar)
 262             nvar = nvar + 1
 263             z(i) = xx(nvar) / scale(nvar)
 264          end if
 265       end do
 266 c
 267 c     adjust atomic coordinates to satisfy distance constraints
 268 c
 269       if (use_rattle)  call shake (x,y,z)
 270 c
 271 c     perform dynamic allocation of some local arrays
 272 c
 273       allocate (derivs(3,n))
 274 c
 275 c     compute and store the energy and gradient
 276 c
 277       if (analytic) then
 278          call gradient (e,derivs)
 279       else
 280          e = energy ()
 281          call numgrad (energy,derivs,eps)
 282       end if
 283       minimiz1 = e
 284 c
 285 c     adjust gradient to remove components along constraints
 286 c
 287       if (use_rattle)  call shakef (derivs)
 288 c
 289 c     convert coordinates and gradient to optimization parameters
 290 c
 291       nvar = 0
 292       do i = 1, n
 293          if (use(i)) then
 294             nvar = nvar + 1
 295             xx(nvar) = x(i) * scale(nvar)
 296             g(nvar) = derivs(1,i) / scale(nvar)
 297             nvar = nvar + 1
 298             xx(nvar) = y(i) * scale(nvar)
 299             g(nvar) = derivs(2,i) / scale(nvar)
 300             nvar = nvar + 1
 301             xx(nvar) = z(i) * scale(nvar)
 302             g(nvar) = derivs(3,i) / scale(nvar)
 303          end if
 304       end do
 305 c
 306 c     perform deallocation of some local arrays
 307 c
 308       deallocate (derivs)
 309       return
 310       end