source/unres/optsave.F90

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