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 build_square

c==============Read Square.dat file=======================================

      include 'square.inc'

c==============Build shell script/comfile for directory creation==========

      if (drive_option(2).ne.'999') then
        call dir_mk(dr_start(1),dr_step(1),i1)
      endif

c==============Build the square-file======================================

      open(2,file='square.bin',access='direct',recl=12,
     1      form='unformatted',status='new',iostat=iostat)
      if (iostat.ne.0) stop 'Could not create SQUARE.BIN'

      if (pos1.ge.900) then
        if (drive_option(1).ne.'999') then
          call line_to_words(dr_atom(1),word,nword)
          read(word(1),*) jj1
          read(word(2),*) jj2
          read(word(3),*) jj3
          read(word(4),*) jj4
        else
          jj1=-1
        endif
        if (drive_option(2).ne.'999') then
          call line_to_words(dr_atom(2),word,nword)
          read(word(1),*) kk1
          read(word(2),*) kk2
          read(word(3),*) kk3
          read(word(4),*) kk4
        else
          kk1=-1
        endif
        call calc_now(resfile,
     1                jj1,jj2,jj3,jj4,kk1,kk2,kk3,kk4,
     2                pos1,pos2)
      endif
      write(*,*) 'Calculation starts at torsions: ',pos1,pos2
      k1=nint((pos1-dr_start(1))/dr_step(1))
      k2=nint((pos2-dr_start(2))/dr_step(2))
      if (k1.lt.0) k1=0
      if (k2.lt.0) k2=0
      if (k1.gt.i1) k1=i1
      if (k2.gt.i2) k2=i2
      rec.first_one=.false.
      rec.first_busy=.false.
      rec.first_done=.false.
      rec.busy1h=.false.
      rec.busy1l=.false.
      rec.busy2h=.false.
      rec.busy2l=.false.
      rec.done1h=.false.
      rec.done1l=.false.
      rec.done2h=.false.
      rec.done2l=.false.
      rec.total_energy=max_energy+100.0
      do j1=0,i1
        do j2=0,i2      
          nr=(i2+1)*j1+j2+1
          write(2,rec=nr) rec
        enddo
      enddo          
      write(*,*) 'Nr of records : ',(i1+1)*(i2+1)
      nr=(i2+1)*k1+k2+1
      write(*,*) 'Enabled record : ',nr
      rec.first_one=.true.
      write(2,rec=nr) rec
      close(2)
      write(*,*) 'Square-file created.'

c========Create the history file==========================================

      call newfile(2,'square.his',iostat)
      if (iostat.ne.0) stop 'Create of history file failed.'
      close(2)
      write(*,*) 'History file created.'
c      call finish(6)
      end

      subroutine calc_now(resfile,j1,j2,j3,j4,k1,k2,k3,k4,r1,r2)
      character*(*) resfile
      character*80 line,word(20)
      common /datax/ coor(3,1000)
      call readfile(45,resfile,iostat)
      if (iostat.ne.0) stop 'Could not open RESFILE'
10    continue
      read(45,'(A)',end=90) line
      if (line(5:5).ne.'.') goto 10            
      read(line,800,end=90) x1,y1,z1,i1,x2,y2,z2,i2
800   format(3f10.5,2x,i3,5x,3f10.5,2x,i3)
      if (i1.gt.0) then
        i=i+1
        coor(1,i)=x1
        coor(2,i)=y1
        coor(3,i)=z1
      endif
      if (i2.gt.0) then
        i=i+1
        coor(1,i)=x2
        coor(2,i)=y2
        coor(3,i)=z2
      endif
      goto 10
90    continue
      close(45)
      if (j1.gt.0) then
        r1=tors(j1,j2,j3,j4)
      endif
      if (k1.gt.0) then
        r2=tors(k1,k2,k3,k4)
      endif      
      return
      end

      function dist(n1,n2)
      common /datax/ coor(3,1000)
      x1=coor(1,n1)
      y1=coor(2,n1)
      z1=coor(3,n1)
      x2=coor(1,n2)
      y2=coor(2,n2)
      z2=coor(3,n2)
      dist=sqrt((x1-x2)**2+(y1-y2)**2+(z1-z2)**2)
      return
      end

      function angl(n1,n2,n3)
      common /datax/ coor(3,1000)
      x1=coor(1,n1)
      y1=coor(2,n1)
      z1=coor(3,n1)
      x2=coor(1,n2)
      y2=coor(2,n2)
      z2=coor(3,n2)
      x3=coor(1,n3)
      y3=coor(2,n3)
      z3=coor(3,n3)
c      write(*,*) x1,y1,z1,x2,y2,z2,x3,y3,z3
c- direction 1
      dx1=x2-x1
      dy1=y2-y1
      dz1=z2-z1
      d=sqrt(dx1**2+dy1**2+dz1**2)
c- direction 2
      dx2=x2-x3
      dy2=y2-y3
      dz2=z2-z3
      d=d*sqrt(dx2**2+dy2**2+dz2**2)
c- dot-product
      dot=dx1*dx2+dy1*dy2+dz1*dz2
c- angle in degrees.
      angl=acosd(dot/d)      
      return
      end

      function tors(n1,n2,n3,n4)
      common /datax/ coor(3,1000)
      x1=coor(1,n1)
      y1=coor(2,n1)
      z1=coor(3,n1)
      x2=coor(1,n2)
      y2=coor(2,n2)
      z2=coor(3,n2)
      x3=coor(1,n3)
      y3=coor(2,n3)
      z3=coor(3,n3)
      x4=coor(1,n4)
      y4=coor(2,n4)
      z4=coor(3,n4)
c- direction 1
      dx1=x2-x1
      dy1=y2-y1
      dz1=z2-z1
c- direction 2
      dx2=x2-x3
      dy2=y2-y3
      dz2=z2-z3
c- direction 3
      dx3=x3-x4
      dy3=y3-y4
      dz3=z3-z4
c- vector-products
      vx1=dy1*dz2-dz1*dy2
      vy1=dz1*dx2-dx1*dz2
      vz1=dx1*dy2-dy1*dx2
      d=sqrt(vx1**2+vy1**2+vz1**2)
      vx2=dy2*dz3-dz2*dy3
      vy2=dz2*dx3-dx2*dz3
      vz2=dx2*dy3-dy2*dx3
      d=d*sqrt(vx2**2+vy2**2+vz2**2)
c- costh
      costh=(vx1*vx2+vy1*vy2+vz1*vz2)/d
c- angle in degrees.
      if (costh.ge.1) then
        t=0
      else if (costh.le.-1) then
        t=180
      else
        t=acosd(-costh)      
      endif
c- sign is ambiguous
      vol=vx2*dx1+vy2*dy1+vz2*dz1
      if (vol.lt.0) t=-t
      tors=t
      return
      end
