********> bugfix.43
Author:   Jim Caldwell
Date:     9/21/95

Programs: Edit

Severity: Moderate

Problem:  Incorrect solvent box with GEN option.


Fix:  Make the following change to gen_sol.f:
---------------------------------------------------------------------------
Replace subroutine read_equil with the following:

      subroutine read_equil(solvs,nsolvres,nsolvat,boxx,boxy,boxz)
      parameter ( isize=500 000 )
      logical skip
      character*4 atnam
#include "sizes.h"
      dimension solvs(3,*)
      dimension scratch(3,isize)
      iat = 0
      write(6,'(t5,a,i8)')'number of atoms per solvent residue: ',
     $   nsolvat
      write(6,'(t5,a,3f10.2)')'requested sizes of solvent: ',
     $   boxx,boxy,boxz
      rewind(54)
      do jn = 1,888888

        read(54,9999,end=100)atnam,tempx,tempy,tempz
9999    format(a4,26x,3f8.3)

        if(atnam.eq.'ATOM')then
          iat = iat + 1 
          solvs(1,iat) = tempx
          solvs(2,iat) = tempy
          solvs(3,iat) = tempz
        endif
      enddo
100   continue

      if(iat.gt.MAXATM)then
          print *,'iat = ',iat,' exceeds MAXATM ',MAXATM
          call mexit(6,1)
      endif
      isolvres = iat/nsolvat
      write(6,'(t5,a,i8)') 'number of read in solvent residues',
     $                     isolvres
      write(6,'(t5,a,i8)') 'number of read in solvent atoms   ',iat

c  move the solvent box to the origin.  Note that the orientation
c  of the solvent box is assumed to be along the x,y,z axes (valid
c  if it came from an Amber run.

      call movecx(iat,solvs,0) 

c  find out how big it is

      bx = 0.0d0
      by = 0.0d0
      bz = 0.0d0
      do i = 1,iat
        tx=max(bx,solvs(1,i))
        ty=max(by,solvs(2,i))
        tz=max(bz,solvs(3,i))

        if(tx.ge.bx)bx=tx
        if(ty.ge.by)by=ty
        if(tz.ge.bz)bz=tz
      enddo
      bx = bx + bx
      by = by + by
      bz = bz + bz
 
      write(6,'(t5,a,3f10.2)') 'solvent box sizes: ',bx,by,bz
      write(6,'(t5,a,3f10.2)') 'needed box sizes:  ',boxx,boxy,boxz

c  if the box read in is too small expand it.

      if(bx.lt.boxx .and. by.lt.boxy .and. bz.lt.boxz)  then

        mat = 0
        do i = 1,3
          do j = 1,3
            do k = 1,3
               shiftx = (i-2)*bx 
               shifty = (j-2)*by
               shiftz = (k-2)*bz
               do ll = 1,iat
                  mat = mat + 1
                  scratch(1,mat) = solvs(1,ll) + shiftx
                  scratch(2,mat) = solvs(2,ll) + shifty 
                  scratch(3,mat) = solvs(3,ll) + shiftz
                  if(mat.gt.isize) then
                    write(6,*)'>>>>>>>>>overflow in scratch',mat,isize
                    call mexit(6,1)
                  endif
               enddo
            enddo
          enddo
        enddo
        newres = mat/nsolvat

c       -- now trim the cube to boxx,boxy,boxz

      else
      
         do jn = 1,iat
             scratch(1,jn) = solvs(1,jn)
             scratch(2,jn) = solvs(2,jn)
             scratch(3,jn) = solvs(3,jn)
         enddo
         newres = iat/nsolvat

      endif

      jcnt = 1
      llim = 0
      kres = 0
      idump = 0
      do i = 1,newres
         skip = .false.
         do j = 1,nsolvat
            if(abs(scratch(1,j+llim)) .gt.boxx/2.0) skip = .true. 
            if(abs(scratch(2,j+llim)) .gt.boxy/2.0) skip = .true. 
            if(abs(scratch(3,j+llim)) .gt.boxz/2.0) skip = .true. 
         enddo
         if(.not.skip) then
            kres = kres + 1
            do j = 1,nsolvat
              solvs(1,jcnt) = scratch(1,j+llim)
              solvs(2,jcnt) = scratch(2,j+llim)
              solvs(3,jcnt) = scratch(3,j+llim)
              jcnt = jcnt + 1
            enddo
            if(jcnt .gt.MAXATM) then
                print *,'jcnt =',jcnt,' exceeds  MAXATM ',MAXATM
                call mexit(6,1)
            endif
         else 
            idump = idump + 1 
         endif
         llim = llim + nsolvat
      enddo
      write(6,1002)kres
1002  format(t5,'a total of ',i5,' residues now in the solvent')
      write(6,1003)jcnt
1003  format(t5,'a total of ',i5,' atoms    now in the solvent')
      nsolvres = kres    
      return 
      end
---------------------------------------------------------------------------

Temporary workarounds: none

--