
c=====================================================================
c ==================================================
c         Copyright And License Agreement
c ==================================================
c
c MMSQ Version 1.0, 4th of april 1991.
c
c You may copy and distribute copies of the complete MMSQ source
c code as you receive it, in any medium, provided that you
c conspicuously and appropriately publish on each file a valid
c copyright notice "Copyright (C) 1991 Rob Hooft, Utrecht
c University"; keep intact the notices on all files that refer to
c this License Agreement and give other recipients of MMSQ a copy
c of this license agreement along with the program. You may charge
c a distribution fee for the physical act of transferring a copy. 
c
c You may copy and distribute MMSQ in object code or executable
c form under the terms mentioned above provided that you also
c accompany it with the complete corresponding machine-readable
c source code, which must also be distributed under the terms
c mentioned above.
c
c You may not copy, sublicense, distribute or transfer MMSQ except 
c as expressly provided under this license agreement.
c
c (This license agreement has been modelled after the 
c GNU GENERAL PUBLIC LICENSE)
c
c=====================================================================
      program mmstart
      character*3 k1s,k2s
      logical not_busy1,not_busy2,not_busy3,not_busy4

c======== Variables and Readin of SQUARE.DAT parameters.==========

      include 'square.inc'

c=================Now create next job input.======================

      call waitbinfile(2)

      energy_low=max_energy
      k1=-1
      k2=-1
      kk1=-1
      kk2=-1
      do j1=0,i1
        do j2=0,i2
          nr=(i2+1)*j1+j2+1
          read(2,rec=nr) rec
          arr(j1,j2)=rec
        enddo
      enddo
      do j1=0,i1
        do j2=0,i2
c------------------------------------------------------------------------
c- what's the best l1 to go to if I want to calculate from j1
c- and what is the energy of the j1 square? Returns 1E9 if none is
c- interesting enough to try,
c- or if all the squares are already being calculated now.
c------------------------------------------------------------------------
          energy=1E9
          l1=-1
          l2=-1
c
c Periodicity of the map is to be taken into account?
c
          j1m1=-1
          if (dr_end(1)-dr_start(1)+dr_step(1).eq.360) j1m1=i1
          if (j1.gt.0) j1m1=j1-1
          j1p1=-1
          if (dr_end(1)-dr_start(1)+dr_step(1).eq.360) j1p1=0
          if (j1.lt.i1) j1p1=j1+1
          j2m1=-1
          if (dr_end(2)-dr_start(2)+dr_step(2).eq.360) j2m1=i2
          if (j2.gt.0) j2m1=j2-1
          j2p1=-1
          if (dr_end(2)-dr_start(2)+dr_step(2).eq.360) j2p1=0
          if (j2.lt.i2) j2p1=j2+1
c
c don't calculate starting from this point if the energy is 
c above requested level
c
          if (arr(j1,j2).total_energy.ge.max_energy) goto 7999
c
c don't calculate starting from this point if it is busy.
c
          if(arr(j1,j2).first_busy.or.
     1       arr(j1,j2).busy1h.or.
     1       arr(j1,j2).busy1l.or.
     1       arr(j1,j2).busy2h.or.
     1       arr(j1,j2).busy2l) goto 7999
c
c collect statistics about neighbours that are busy
c
          not_busy1=.true.
          if (j1m1.ne.-1) then
          if(arr(j1m1,j2).first_busy.or.
     1       arr(j1m1,j2).busy1h.or.
     1       arr(j1m1,j2).busy1l.or.
     1       arr(j1m1,j2).busy2h.or.
     1       arr(j1m1,j2).busy2l) not_busy1=.false.
          endif
          not_busy2=.true.
          if (j1p1.ne.-1) then
          if(arr(j1p1,j2).first_busy.or.
     1       arr(j1p1,j2).busy1h.or.
     1       arr(j1p1,j2).busy1l.or.
     1       arr(j1p1,j2).busy2h.or.
     1       arr(j1p1,j2).busy2l) not_busy2=.false.
          endif
          not_busy3=.true.
          if (j2m1.ne.-1) then
          if(arr(j1,j2m1).first_busy.or.
     1       arr(j1,j2m1).busy1h.or.
     1       arr(j1,j2m1).busy1l.or.
     1       arr(j1,j2m1).busy2h.or.
     1       arr(j1,j2m1).busy2l) not_busy3=.false.
          endif
          not_busy4=.true.
          if (j2p1.ne.-1) then
          if(arr(j1,j2p1).first_busy.or.
     1       arr(j1,j2p1).busy1h.or.
     1       arr(j1,j2p1).busy1l.or.
     1       arr(j1,j2p1).busy2h.or.
     1       arr(j1,j2p1).busy2l) not_busy4=.false.
          endif
c
c ok now the reasoning starts.
c select the neighbouring point with the highest energy, save in l1
c it must have an energy 'min_jump' higher than this point itself.
c
          energy=arr(j1,j2).total_energy+min_jump
          l1=-1
          if (j1m1.ne.-1) then
            if ((.not.arr(j1m1,j2).done1h).and.not_busy1.and.
     1           arr(j1m1,j2).total_energy.gt.energy) then
              l1=j1m1
              l2=j2
              energy=arr(l1,l2).total_energy
            endif
          endif         
          if (j1p1.ne.-1) then
            if ((.not.arr(j1p1,j2).done1l).and.not_busy2.and.
     1           arr(j1p1,j2).total_energy.gt.energy) then
              l1=j1p1
              l2=j2
              energy=arr(l1,l2).total_energy
            endif
          endif         
          if (j2m1.ne.-1) then
            if ((.not.arr(j1,j2m1).done2h).and.not_busy3.and.
     1           arr(j1,j2m1).total_energy.gt.energy) then
              l1=j1
              l2=j2m1
              energy=arr(l1,l2).total_energy
            endif
          endif         
          if (j2p1.ne.-1) then
            if ((.not.arr(j1,j2p1).done2l).and.not_busy4.and.
     1           arr(j1,j2p1).total_energy.gt.energy) then
              l1=j1
              l2=j2p1
              energy=arr(l1,l2).total_energy
            endif
          endif         
          energy=arr(j1,j2).total_energy
c
c did we find a point that we can improve on? no? 
c then change energy to 1E9 again.
c
          if (l1.eq.-1) energy=1E9

7999      continue
c
c is this an improvement on the best point found?
c
          if (energy.lt.energy_low) then
            energy_low=energy
            k1=l1
            k2=l2
            kk1=j1
            kk2=j2
            val1=dr_start(1)+kk1*dr_step(1)
            val2=dr_start(2)+kk2*dr_step(2)
            n1=nint(val1)
            n2=nint(val2)
            if (n1.lt.0) n1=n1+360
            if (n2.lt.0) n2=n2+360
            call newjob(resfile,n1,n2,i2)
            resfile=resfile(1:lenstr(resfile))//'.res'
c
c did we run the first calculation yet?
c (This could have been tested earlier, but who cares...)
c
          else if (arr(j1,j2).first_one) then
            if (.not.arr(j1,j2).first_busy.and.
     1          .not.arr(j1,j2).first_done) then
              k1=j1
              k2=j2
            endif
          endif              
        enddo
      enddo
c
c Did we assign a possible calculation?
c
      if (k1.eq.-1) stop 'nothing left to do'
      write(*,678) ' Now calculating ',k1,k2,' from ',kk1,kk2
678   format(2(a,2i3))
c
c Modify information in the binary file.
c
      nr=(i2+1)*k1+k2+1
      read(2,rec=nr) rec
      if (kk1.eq.-1) then
        rec.first_busy=.true.
      else if (kk2.eq.k2+1) then
        rec.busy2h=.true.
      else if (kk2.eq.k2-1) then
        rec.busy2l=.true.
      else if (kk1.eq.k1+1) then
        rec.busy1h=.true.
      else if (kk1.eq.k1-1) then
        rec.busy1l=.true.
      else if (kk1.eq.0.and.k1.eq.i1) then
        rec.busy1h=.true.
      else if (kk2.eq.0.and.k2.eq.i2) then
        rec.busy2h=.true.
      else if (k1.eq.0.and.kk1.eq.i1) then
        rec.busy1l=.true.
      else if (k2.eq.0.and.kk2.eq.i2) then
        rec.busy2l=.true.
      else
        print *,'There is something really wrong.'
        print *,'I couldn''t assign a direction!'
c
c- create a nice VMS fatal error
c
        call exit(1572)
      endif
      write(2,rec=nr) rec
      close(2)
c
c Ok. We found something to do, now we create the job-file
c
 600  continue
      val1=dr_start(1)+k1*dr_step(1)
      val2=dr_start(2)+k2*dr_step(2)
      n1=nint(val1)
      n2=nint(val2)
      if (n1.lt.0) n1=n1+360
      if (n2.lt.0) n2=n2+360
      write(k1s,'(I3.3)') k1
      write(k2s,'(I3.3)') k2
      call newjob(jobnam,n1,n2,i2)
      call set_symbol('newjob',jobnam(1:lenstr(jobnam)))
      l=lenstr(jobnam)
      call set_symbol('jobi',k1s)
      call set_symbol('jobj',k2s)

      OPEN (1,FILE=resfile,READONLY,STATUS='OLD',ERR=999)
      OPEN (2,FILE=datfile,READONLY,STATUS='OLD',ERR=997)
      call newfile(3,jobnam(1:l)//'.dat',iostat)
      if (iostat.ne.0) goto 996
c
c first copy the DATFILE to the new job.
c
110   continue
      read(2,3456,end=120) line
3456  format(A)
      call upcase(line,line2)
      if (line.eq.'END') goto 120 ! suppress END card.
      if (line2(1:4).ne.'ATOM') then
        write(3,3456) line(1:lenstr(line))
        goto 110
      endif
120   continue
c
c second: write the DRIVE command generated.
c
      write(3,3457) dr_atom(1)(1:lenstr(dr_atom(1))),
     1        val1,val1,drive_option(1)
      if (drive_option(2).ne.'999') then
        write(3,3457) dr_atom(2)(1:lenstr(dr_atom(2))),
     1        val2,val2,drive_option(2)
      endif
 3457 format('DRIVE ',a,2f10.5,' 10. ',a)
      i=0
c
c Third: copy atomic coordinates from the selected RESfile.
c copy charges from the DATFILE. (also when they're not there)
c this is why NDC 0 OVERRULE should allways be specified when
c no charges are desired.
c
      read(1,3456,end=40) line1
30    continue
      read(1,3456,end=40,err=30) resline
c
c skip non-atom records in selected RESfile
c
      if (resline(5:5).ne.'.') goto 30
c
c atom record can contain 1 or 2 consecutive atoms.
c this is determined by testing value of i1 and i2.
c
      read(resline,800,err=30) x1,y1,z1,i1,x2,y2,z2,i2
      if (i1.gt.0) then
        i=i+1
        call line_to_words(line,word,nword)
        read(word(2),*,iostat=iostat) rr
c
c get charge for this atom, if applicable
c
        if (iostat.eq.0) then
          read(word(6),*,iostat=iostat) charge
          if (nword.lt.6) charge=0.0
          if (iostat.ne.0) charge=0.0
        else
          read(word(7),*,iostat=iostat) charge
          if (nword.lt.7) charge=0.0
          if (iostat.ne.0) charge=0.0
        endif
c
c find next atom in DATfile
c
 44     continue
        read(2,3456,iostat=iostat) line
        if (iostat.eq.0.and.line(1:4).ne.'ATOM') goto 44
c
c put info for this (second) atom to DAT file
c
        la=lab(i1,i)
        write(3,801) la,x1,y1,z1,i1,charge
      endif
      if (i2.gt.0) then
        i=i+1
        call line_to_words(line,word,nword)
c
c get charge for this (second) atom, if applicable
c
        read(word(2),*,iostat=iostat) rr
        if (iostat.eq.0) then
          read(word(6),*,iostat=iostat) charge
          if (iostat.ne.0) charge=0.0
        else
          read(word(7),*,iostat=iostat) charge
          if (iostat.ne.0) charge=0.0
        endif
c
c find next atom in DATfile
c
 46     continue
        read(2,3456,iostat=iostat) line
        if (iostat.eq.0.and.line(1:4).ne.'ATOM') goto 46
c
c put info for this (second) atom to DAT file
c
        la=lab(i2,i)
        write(3,801) la,x2,y2,z2,i2,charge
      endif
800   format(3f10.5,2x,i3,5x,3f10.5,2x,i3)
801   format('ATOM ',a5,3(x,f10.5),x,i5,x,f10.5)
      goto 30
40    continue
      close(1)
      close(2)
      close(3)
c      call finish(6)
      call exit(1)
999   stop 'Unable to open selected RESfile'
997   stop 'Unable to open DATFILE'
996   stop 'Unable to open the new jobfile'
      end

      character*5 function lab(ityp,i)
      parameter(maxtyp=28)
      character*2 XTYP(Maxtyp)
      data xtyp/'C ','C ','C ','C ','H ','O ','O ','N ','N ','N ',
     1          'F ','CL','BR','I ','S ','S ','S ','S ','SI','LP',
     2          'H ','C ','H ','H ','P ','B ','B ','H '/
      character*5 c      
      write(c,'(I5)') i
      if (ityp.le.maxtyp) then
        c(1:2)=xtyp(ityp)
      else
        c(1:2)='??'
      endif
      do k=4,2,-1
        if (c(k:k).eq.' ') c=c(1:k-1)//c(k+1:)
      enddo
      lab=c
      return
      end      

