source/unres/src-HCD-5D/optsave.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     ##  subroutine optsave  --  save optimization info and results  ##
  11 c     ##                                                              ##
  12 c     ##################################################################
  13 c
  14 c
  15 c     "optsave" is used by the optimizers to write imtermediate
  16 c     coordinates and other relevant information; also checks for
  17 c     user requested termination of an optimization
  18 c
  19 c
  20       subroutine optsave (ncycle,f,xx)
  21       use atomid
  22       use atoms
  23       use bound
  24       use deriv
  25       use files
  26       use iounit
  27       use math
  28       use mpole
  29       use omega
  30       use output
  31       use polar
  32       use potent
  33       use scales
  34       use socket
  35       use titles
  36       use units
  37       use usage
  38       use zcoord
  39       implicit none
  40       integer i,j,k,lext
  41       integer iopt,ifrc
  42       integer iind,iend
  43       integer ncycle,nvar
  44       integer freeunit
  45       integer trimtext
  46       real*8 f,xx(*)
  47       logical exist
  48       character*7 ext
  49       character*240 optfile
  50       character*240 frcfile
  51       character*240 indfile
  52       character*240 endfile
  53 c
  54 c
  55 c     nothing to do if coordinate type is undefined
  56 c
  57       if (coordtype .eq. 'NONE')  return
  58 c
  59 c     check scaling factors for optimization parameters
  60 c
  61       if (.not. set_scale) then
  62          set_scale = .true.
  63          if (coordtype .eq. 'CARTESIAN') then
  64             if (.not. allocated(scale))  allocate (scale(3*n))
  65             do i = 1, 3*n
  66                scale(i) = 1.0d0
  67             end do
  68          else if (coordtype .eq. 'INTERNAL') then
  69             if (.not. allocated(scale))  allocate (scale(nomega))
  70             do i = 1, nomega
  71                scale(i) = 1.0d0
  72             end do
  73          end if
  74       end if
  75 c
  76 c     convert optimization parameters to atomic coordinates
  77 c
  78       if (coordtype .eq. 'CARTESIAN') then
  79          nvar = 0
  80          do i = 1, n
  81             if (use(i)) then
  82                nvar = nvar + 1
  83                x(i) = xx(nvar) / scale(nvar)
  84                nvar = nvar + 1
  85                y(i) = xx(nvar) / scale(nvar)
  86                nvar = nvar + 1
  87                z(i) = xx(nvar) / scale(nvar)
  88             end if
  89          end do
  90          if (use_bounds)  call bounds
  91       else if (coordtype .eq. 'INTERNAL') then
  92          do i = 1, nomega
  93             dihed(i) = xx(i) / scale(i)
  94             ztors(zline(i)) = dihed(i) * radian
  95          end do
  96       end if
  97 c
  98 c     get name of archive or intermediate coordinates file
  99 c
 100       iopt = freeunit ()
 101       if (cyclesave) then
 102          if (archive) then
 103             optfile = filename(1:leng)
 104             call suffix (optfile,'arc','old')
 105             inquire (file=optfile,exist=exist)
 106             if (exist) then
 107                call openend (iopt,optfile)
 108             else
 109                open (unit=iopt,file=optfile,status='new')
 110             end if
 111          else
 112             lext = 3
 113             call numeral (ncycle,ext,lext)
 114             optfile = filename(1:leng)//'.'//ext(1:lext)
 115             call version (optfile,'new')
 116             open (unit=iopt,file=optfile,status='new')
 117          end if
 118       else
 119          optfile = outfile
 120          call version (optfile,'old')
 121          open (unit=iopt,file=optfile,status='old')
 122          rewind (unit=iopt)
 123       end if
 124 c
 125 c     update intermediate file with desired coordinate type
 126 c
 127       if (coordtype .eq. 'CARTESIAN') then
 128          call prtxyz (iopt)
 129       else if (coordtype .eq. 'INTERNAL') then
 130          call prtint (iopt)
 131       else if (coordtype .eq. 'RIGIDBODY') then
 132          call prtxyz (iopt)
 133       end if
 134       close (unit=iopt)
 135 c
 136 c     save the force vector components for the current step
 137 c
 138       if (frcsave .and. coordtype.eq.'CARTESIAN') then
 139          ifrc = freeunit ()
 140          if (archive) then
 141             frcfile = filename(1:leng)
 142             call suffix (frcfile,'frc','old')
 143             inquire (file=frcfile,exist=exist)
 144             if (exist) then
 145                call openend (ifrc,frcfile)
 146             else
 147                open (unit=ifrc,file=frcfile,status='new')
 148             end if
 149          else
 150             frcfile = filename(1:leng)//'.'//ext(1:lext)//'f'
 151             call version (frcfile,'new')
 152             open (unit=ifrc,file=frcfile,status='new')
 153          end if
 154          write (ifrc,250)  n,title(1:ltitle)
 155   250    format (i6,2x,a)
 156          do i = 1, n
 157             write (ifrc,260)  i,name(i),(-desum(j,i),j=1,3)
 158   260       format (i6,2x,a3,3x,d13.6,3x,d13.6,3x,d13.6)
 159          end do
 160          close (unit=ifrc)
 161          write (iout,270)  frcfile(1:trimtext(frcfile))
 162   270    format (' Force Vector File',11x,a)
 163       end if
 164 c
 165 c     save the current induced dipole moment at each site
 166 c
 167       if (uindsave .and. use_polar .and. coordtype.eq.'CARTESIAN') then
 168          iind = freeunit ()
 169          if (archive) then
 170             indfile = filename(1:leng)
 171             call suffix (indfile,'uind','old')
 172             inquire (file=indfile,exist=exist)
 173             if (exist) then
 174                call openend (iind,indfile)
 175             else
 176                open (unit=iind,file=indfile,status='new')
 177             end if
 178          else
 179             indfile = filename(1:leng)//'.'//ext(1:lext)//'u'
 180             call version (indfile,'new')
 181             open (unit=iind,file=indfile,status='new')
 182          end if
 183          write (iind,280)  n,title(1:ltitle)
 184   280    format (i6,2x,a)
 185          do i = 1, npole
 186             if (polarity(i) .ne. 0.0d0) then
 187                k = ipole(i)
 188                write (iind,290)  k,name(k),(debye*uind(j,i),j=1,3)
 189   290          format (i6,2x,a3,3f12.6)
 190             end if
 191          end do
 192          close (unit=iind)
 193          write (iout,300)  indfile(1:trimtext(indfile))
 194   300    format (' Induced Dipole File',10x,a)
 195       end if
 196 c
 197 c     send data via external socket communication if desired
 198 c
 199       if (.not.sktstart .or. use_socket) then
 200          if (coordtype .eq. 'INTERNAL')  call makexyz
 201          call sktopt (ncycle,f)
 202       end if
 203 c
 204 c     test for requested termination of the optimization
 205 c
 206       endfile = 'tinker.end'
 207       inquire (file=endfile,exist=exist)
 208       if (.not. exist) then
 209          endfile = filename(1:leng)//'.end'
 210          inquire (file=endfile,exist=exist)
 211          if (exist) then
 212             iend = freeunit ()
 213             open (unit=iend,file=endfile,status='old')
 214             close (unit=iend,status='delete')
 215          end if
 216       end if
 217       if (exist) then
 218          write (iout,10)
 219    10    format (/,' OPTSAVE  --  Optimization Calculation Ending',
 220      &              ' due to User Request')
 221          call fatal
 222       end if
 223       return
 224       end