X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD-NEWSC%2Flocal_move.f;fp=source%2Funres%2Fsrc_MD-NEWSC%2Flocal_move.f;h=7a7e125ae13acf10b51e20f0ec995ebcb799196c;hb=7308760ff07636ef6b1ee28d8c3a67a23c14b34b;hp=0000000000000000000000000000000000000000;hpb=9a54ab407f6d0d9d564d52763b3e2136450b9ffc;p=unres.git

diff --git a/source/unres/src_MD-NEWSC/local_move.f b/source/unres/src_MD-NEWSC/local_move.f
new file mode 100644
index 0000000..7a7e125
--- /dev/null
+++ b/source/unres/src_MD-NEWSC/local_move.f
@@ -0,0 +1,972 @@
+c-------------------------------------------------------------
+
+      subroutine local_move_init(debug)
+crc      implicit none
+
+c     Includes
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'  ! Needed by COMMON.LOCAL
+      include 'COMMON.GEO'  ! For pi, deg2rad
+      include 'COMMON.LOCAL'  ! For vbl
+      include 'COMMON.LOCMOVE'
+
+c     INPUT arguments
+      logical debug
+
+
+c     Determine wheter to do some debugging output
+      locmove_output=debug
+
+c     Set the init_called flag to 1
+      init_called=1
+
+c     The following are never changed
+      min_theta=60.D0*deg2rad  ! (0,PI)
+      max_theta=175.D0*deg2rad  ! (0,PI)
+      dmin2=vbl*vbl*2.*(1.-cos(min_theta))
+      dmax2=vbl*vbl*2.*(1.-cos(max_theta))
+      flag=1.0D300
+      small=1.0D-5
+      small2=0.5*small*small
+
+c     Not really necessary...
+      a_n=0
+      b_n=0
+      res_n=0
+
+      return
+      end
+
+c-------------------------------------------------------------
+
+      subroutine local_move(n_start, n_end, PHImin, PHImax)
+c     Perform a local move between residues m and n (inclusive)
+c     PHImin and PHImax [0,PI] determine the size of the move
+c     Works on whatever structure is in the variables theta and phi,
+c     sidechain variables are left untouched
+c     The final structure is NOT minimized, but both the cartesian
+c     variables c and the angles are up-to-date at the end (no further
+c     chainbuild is required)
+crc      implicit none
+
+c     Includes
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+      include 'COMMON.GEO'
+      include 'COMMON.CHAIN'
+      include 'COMMON.VAR'
+      include 'COMMON.MINIM'
+      include 'COMMON.SBRIDGE'
+      include 'COMMON.LOCMOVE'
+
+c     External functions
+      integer move_res
+      external move_res
+      double precision ran_number
+      external ran_number
+
+c     INPUT arguments
+      integer n_start, n_end  ! First and last residues to move
+      double precision PHImin, PHImax  ! min/max angles [0,PI]
+
+c     Local variables
+      integer i,j
+      double precision min,max
+      integer iretcode
+
+
+c     Check if local_move_init was called.  This assumes that it
+c     would not be 1 if not explicitely initialized
+      if (init_called.ne.1) then
+        write(6,*)'   ***   local_move_init not called!!!'
+        stop
+      endif
+
+c     Quick check for crazy range
+      if (n_start.gt.n_end .or. n_start.lt.1 .or. n_end.gt.nres) then
+        write(6,'(a,i3,a,i3)')
+     +       '   ***   Cannot make local move between n_start = ',
+     +       n_start,' and n_end = ',n_end
+        return
+      endif
+
+c     Take care of end residues first...
+      if (n_start.eq.1) then
+c     Move residue 1 (completely random)
+        theta(3)=ran_number(min_theta,max_theta)
+        phi(4)=ran_number(-PI,PI)
+        i=2
+      else
+        i=n_start
+      endif
+      if (n_end.eq.nres) then
+c     Move residue nres (completely random)
+        theta(nres)=ran_number(min_theta,max_theta)
+        phi(nres)=ran_number(-PI,PI)
+        j=nres-1
+      else
+        j=n_end
+      endif
+
+c     ...then go through all other residues one by one
+c     Start from the two extremes and converge
+      call chainbuild
+      do while (i.le.j)
+        min=PHImin
+        max=PHImax
+c$$$c     Move the first two residues by less than the others
+c$$$        if (i-n_start.lt.3) then
+c$$$          if (i-n_start.eq.0) then
+c$$$            min=0.4*PHImin
+c$$$            max=0.4*PHImax
+c$$$          else if (i-n_start.eq.1) then
+c$$$            min=0.8*PHImin
+c$$$            max=0.8*PHImax
+c$$$          else if (i-n_start.eq.2) then
+c$$$            min=PHImin
+c$$$            max=PHImax
+c$$$          endif
+c$$$        endif
+
+c     The actual move, on residue i
+        iretcode=move_res(min,max,i)  ! Discard iretcode
+        i=i+1
+
+        if (i.le.j) then
+          min=PHImin
+          max=PHImax
+c$$$c     Move the last two residues by less than the others
+c$$$          if (n_end-j.lt.3) then
+c$$$            if (n_end-j.eq.0) then
+c$$$              min=0.4*PHImin
+c$$$              max=0.4*PHImax
+c$$$            else if (n_end-j.eq.1) then
+c$$$              min=0.8*PHImin
+c$$$              max=0.8*PHImax
+c$$$            else if (n_end-j.eq.2) then
+c$$$              min=PHImin
+c$$$              max=PHImax
+c$$$            endif
+c$$$          endif
+
+c     The actual move, on residue j
+          iretcode=move_res(min,max,j)  ! Discard iretcode
+          j=j-1
+        endif
+      enddo
+
+      call int_from_cart(.false.,.false.)
+
+      return
+      end
+
+c-------------------------------------------------------------
+
+      subroutine output_tabs
+c     Prints out the contents of a_..., b_..., res_...
+      implicit none
+
+c     Includes
+      include 'COMMON.GEO'
+      include 'COMMON.LOCMOVE'
+
+c     Local variables
+      integer i,j
+
+
+      write(6,*)'a_...'
+      write(6,'(8f7.1)')(a_ang(i)*rad2deg,i=0,a_n-1)
+      write(6,'(8(2x,3l1,2x))')((a_tab(i,j),i=0,2),j=0,a_n-1)
+
+      write(6,*)'b_...'
+      write(6,'(4f7.1)')(b_ang(i)*rad2deg,i=0,b_n-1)
+      write(6,'(4(2x,3l1,2x))')((b_tab(i,j),i=0,2),j=0,b_n-1)
+
+      write(6,*)'res_...'
+      write(6,'(12f7.1)')(res_ang(i)*rad2deg,i=0,res_n-1)
+      write(6,'(12(2x,3l1,2x))')((res_tab(0,i,j),i=0,2),j=0,res_n-1)
+      write(6,'(12(2x,3l1,2x))')((res_tab(1,i,j),i=0,2),j=0,res_n-1)
+      write(6,'(12(2x,3l1,2x))')((res_tab(2,i,j),i=0,2),j=0,res_n-1)
+
+      return
+      end
+
+c-------------------------------------------------------------
+
+      subroutine angles2tab(PHImin,PHImax,n,ang,tab)
+c     Only uses angles if [0,PI] (but PHImin cannot be 0.,
+c     and PHImax cannot be PI)
+      implicit none
+
+c     Includes
+      include 'COMMON.GEO'
+
+c     INPUT arguments
+      double precision PHImin,PHImax
+
+c     OUTPUT arguments
+      integer n
+      double precision ang(0:3)
+      logical tab(0:2,0:3)
+
+
+      if (PHImin .eq. PHImax) then
+c     Special case with two 010's
+        n = 2;
+        ang(0) = -PHImin;
+        ang(1) = PHImin;
+        tab(0,0) = .false.
+        tab(2,0) = .false.
+        tab(0,1) = .false.
+        tab(2,1) = .false.
+        tab(1,0) = .true.
+        tab(1,1) = .true.
+      else if (PHImin .eq. PI) then
+c     Special case with one 010
+        n = 1
+        ang(0) = PI
+        tab(0,0) = .false.
+        tab(2,0) = .false.
+        tab(1,0) = .true.
+      else if (PHImax .eq. 0.) then
+c     Special case with one 010
+        n = 1
+        ang(0) = 0.
+        tab(0,0) = .false.
+        tab(2,0) = .false.
+        tab(1,0) = .true.
+      else
+c     Standard cases
+        n = 0
+        if (PHImin .gt. 0.) then
+c     Start of range (011)
+          ang(n) = PHImin
+          tab(0,n) = .false.
+          tab(1,n) = .true.
+          tab(2,n) = .true.
+c     End of range (110)
+          ang(n+1) = -PHImin
+          tab(0,n+1) = .true.
+          tab(1,n+1) = .true.
+          tab(2,n+1) = .false.
+          n = n+2
+        endif
+        if (PHImax .lt. PI) then
+c     Start of range (011)
+          ang(n) = -PHImax
+          tab(0,n) = .false.
+          tab(1,n) = .true.
+          tab(2,n) = .true.
+c     End of range (110)
+          ang(n+1) = PHImax
+          tab(0,n+1) = .true.
+          tab(1,n+1) = .true.
+          tab(2,n+1) = .false.
+          n = n+2
+        endif
+      endif
+
+      return
+      end
+
+c-------------------------------------------------------------
+
+      subroutine minmax_angles(x,y,z,r,n,ang,tab)
+c     When solutions do not exist, assume all angles
+c     are acceptable - i.e., initial geometry must be correct
+      implicit none
+
+c     Includes
+      include 'COMMON.GEO'
+      include 'COMMON.LOCMOVE'
+
+c     Input arguments
+      double precision x,y,z,r
+
+c     Output arguments
+      integer n
+      double precision ang(0:3)
+      logical tab(0:2,0:3)
+
+c     Local variables
+      double precision num, denom, phi
+      double precision Kmin, Kmax
+      integer i
+
+
+      num = x*x + y*y + z*z
+      denom = x*x + y*y
+      n = 0
+      if (denom .gt. 0.) then
+        phi = atan2(y,x)
+        denom = 2.*r*sqrt(denom)
+        num = num+r*r
+        Kmin = (num - dmin2)/denom
+        Kmax = (num - dmax2)/denom
+
+c     Allowed values of K (else all angles are acceptable)
+c     -1 <= Kmin <  1
+c     -1 <  Kmax <= 1
+        if (Kmin .gt. 1. .or. abs(Kmin-1.) .lt. small2) then
+          Kmin = -flag
+        else if (Kmin .lt. -1. .or. abs(Kmin+1.) .lt. small2) then
+          Kmin = PI
+        else
+          Kmin = acos(Kmin)
+        endif
+
+        if (Kmax .lt. -1. .or. abs(Kmax+1.) .lt. small2) then
+          Kmax = flag
+        else if (Kmax .gt. 1. .or. abs(Kmax-1.) .lt. small2) then
+          Kmax = 0.
+        else
+          Kmax = acos(Kmax)
+        endif
+
+        if (Kmax .lt. Kmin) Kmax = Kmin
+
+        call angles2tab(Kmin, Kmax, n, ang, tab)
+
+c     Add phi and check that angles are within range (-PI,PI]
+        do i=0,n-1
+          ang(i) = ang(i)+phi
+          if (ang(i) .le. -PI) then
+            ang(i) = ang(i)+2.*PI
+          else if (ang(i) .gt. PI) then
+            ang(i) = ang(i)-2.*PI
+          endif
+        enddo
+      endif
+
+      return
+      end
+
+c-------------------------------------------------------------
+
+      subroutine construct_tab
+c     Take a_... and b_... values and produces the results res_...
+c     x_ang are assumed to be all different (diff > small)
+c     x_tab(1,i) must be 1 for all i (i.e., all x_ang are acceptable)
+      implicit none
+
+c     Includes
+      include 'COMMON.LOCMOVE'
+
+c     Local variables
+      integer n_max,i,j,index
+      logical done
+      double precision phi
+
+
+      n_max = a_n + b_n
+      if (n_max .eq. 0) then
+        res_n = 0
+        return
+      endif
+
+      do i=0,n_max-1
+        do j=0,1
+          res_tab(j,0,i) = .true.
+          res_tab(j,2,i) = .true.
+          res_tab(j,1,i) = .false.
+        enddo
+      enddo
+
+      index = 0
+      phi = -flag
+      done = .false.
+      do while (.not.done)
+        res_ang(index) = flag
+
+c     Check a first...
+        do i=0,a_n-1
+          if ((a_ang(i)-phi).gt.small .and.
+     +         a_ang(i) .lt. res_ang(index)) then
+c     Found a lower angle
+            res_ang(index) = a_ang(i)
+c     Copy the values from a_tab into res_tab(0,,)
+            res_tab(0,0,index) = a_tab(0,i)
+            res_tab(0,1,index) = a_tab(1,i)
+            res_tab(0,2,index) = a_tab(2,i)
+c     Set default values for res_tab(1,,)
+            res_tab(1,0,index) = .true.
+            res_tab(1,1,index) = .false.
+            res_tab(1,2,index) = .true.
+          else if (abs(a_ang(i)-res_ang(index)).lt.small) then
+c     Found an equal angle (can only be equal to a b_ang)
+            res_tab(0,0,index) = a_tab(0,i)
+            res_tab(0,1,index) = a_tab(1,i)
+            res_tab(0,2,index) = a_tab(2,i)
+          endif
+        enddo
+c     ...then check b
+        do i=0,b_n-1
+          if ((b_ang(i)-phi).gt.small .and.
+     +         b_ang(i) .lt. res_ang(index)) then
+c     Found a lower angle
+            res_ang(index) = b_ang(i)
+c     Copy the values from b_tab into res_tab(1,,)
+            res_tab(1,0,index) = b_tab(0,i)
+            res_tab(1,1,index) = b_tab(1,i)
+            res_tab(1,2,index) = b_tab(2,i)
+c     Set default values for res_tab(0,,)
+            res_tab(0,0,index) = .true.
+            res_tab(0,1,index) = .false.
+            res_tab(0,2,index) = .true.
+          else if (abs(b_ang(i)-res_ang(index)).lt.small) then
+c     Found an equal angle (can only be equal to an a_ang)
+            res_tab(1,0,index) = b_tab(0,i)
+            res_tab(1,1,index) = b_tab(1,i)
+            res_tab(1,2,index) = b_tab(2,i)
+          endif
+        enddo
+
+        if (res_ang(index) .eq. flag) then
+          res_n = index
+          done = .true.
+        else if (index .eq. n_max-1) then
+          res_n = n_max
+          done = .true.
+        else
+          phi = res_ang(index)  ! Store previous angle
+          index = index+1
+        endif
+      enddo
+
+c     Fill the gaps
+c     First a...
+      index = 0
+      if (a_n .gt. 0) then
+        do while (.not.res_tab(0,1,index))
+          index=index+1
+        enddo
+        done = res_tab(0,2,index)
+        do i=index+1,res_n-1
+          if (res_tab(0,1,i)) then
+            done = res_tab(0,2,i)
+          else
+            res_tab(0,0,i) = done
+            res_tab(0,1,i) = done
+            res_tab(0,2,i) = done
+          endif
+        enddo
+        done = res_tab(0,0,index)
+        do i=index-1,0,-1
+          if (res_tab(0,1,i)) then
+            done = res_tab(0,0,i)
+          else
+            res_tab(0,0,i) = done
+            res_tab(0,1,i) = done
+            res_tab(0,2,i) = done
+          endif
+        enddo
+      else
+        do i=0,res_n-1
+          res_tab(0,0,i) = .true.
+          res_tab(0,1,i) = .true.
+          res_tab(0,2,i) = .true.
+        enddo
+      endif
+c     ...then b
+      index = 0
+      if (b_n .gt. 0) then
+        do while (.not.res_tab(1,1,index))
+          index=index+1
+        enddo
+        done = res_tab(1,2,index)
+        do i=index+1,res_n-1
+          if (res_tab(1,1,i)) then
+            done = res_tab(1,2,i)
+          else
+            res_tab(1,0,i) = done
+            res_tab(1,1,i) = done
+            res_tab(1,2,i) = done
+          endif
+        enddo
+        done = res_tab(1,0,index)
+        do i=index-1,0,-1
+          if (res_tab(1,1,i)) then
+            done = res_tab(1,0,i)
+          else
+            res_tab(1,0,i) = done
+            res_tab(1,1,i) = done
+            res_tab(1,2,i) = done
+          endif
+        enddo
+      else
+        do i=0,res_n-1
+          res_tab(1,0,i) = .true.
+          res_tab(1,1,i) = .true.
+          res_tab(1,2,i) = .true.
+        enddo
+      endif
+
+c     Finally fill the last row with AND operation
+      do i=0,res_n-1
+        do j=0,2
+          res_tab(2,j,i) = (res_tab(0,j,i) .and. res_tab(1,j,i))
+        enddo
+      enddo
+
+      return 
+      end
+
+c-------------------------------------------------------------
+
+      subroutine construct_ranges(phi_n,phi_start,phi_end)
+c     Given the data in res_..., construct a table of 
+c     min/max allowed angles
+      implicit none
+
+c     Includes
+      include 'COMMON.GEO'
+      include 'COMMON.LOCMOVE'
+
+c     Output arguments
+      integer phi_n
+      double precision phi_start(0:11),phi_end(0:11)
+
+c     Local variables
+      logical done
+      integer index
+
+
+      if (res_n .eq. 0) then
+c     Any move is allowed
+        phi_n = 1
+        phi_start(0) = -PI
+        phi_end(0) = PI
+      else
+        phi_n = 0
+        index = 0
+        done = .false.
+        do while (.not.done)
+c     Find start of range (01x)
+          done = .false.
+          do while (.not.done)
+            if (res_tab(2,0,index).or.(.not.res_tab(2,1,index))) then
+              index=index+1
+            else
+              done = .true.
+              phi_start(phi_n) = res_ang(index)
+            endif
+            if (index .eq. res_n) done = .true.
+          enddo
+c     If a start was found (index < res_n), find the end of range (x10)
+c     It may not be found without wrapping around
+          if (index .lt. res_n) then
+            done = .false.
+            do while (.not.done)
+              if ((.not.res_tab(2,1,index)).or.res_tab(2,2,index)) then
+                index=index+1
+              else
+                done = .true.
+              endif
+              if (index .eq. res_n) done = .true.
+            enddo
+            if (index .lt. res_n) then
+c     Found the end of the range
+              phi_end(phi_n) = res_ang(index)
+              phi_n=phi_n+1
+              index=index+1
+              if (index .eq. res_n) then
+                done = .true.
+              else
+                done = .false.
+              endif
+            else
+c     Need to wrap around
+              done = .true.
+              phi_end(phi_n) = flag
+            endif
+          endif
+        enddo
+c     Take care of the last one if need to wrap around
+        if (phi_end(phi_n) .eq. flag) then
+          index = 0
+          do while ((.not.res_tab(2,1,index)).or.res_tab(2,2,index))
+            index=index+1
+          enddo
+          phi_end(phi_n) = res_ang(index) + 2.*PI
+          phi_n=phi_n+1
+        endif
+      endif
+
+      return
+      end
+
+c-------------------------------------------------------------
+
+      subroutine fix_no_moves(phi)
+      implicit none
+
+c     Includes
+      include 'COMMON.GEO'
+      include 'COMMON.LOCMOVE'
+
+c     Output arguments
+      double precision phi
+
+c     Local variables
+      integer index
+      double precision diff,temp
+
+
+c     Look for first 01x in gammas (there MUST be at least one)
+      diff = flag
+      index = 0
+      do while (res_tab(1,0,index) .or. (.not.res_tab(1,1,index)))
+        index=index+1
+      enddo
+      if (res_ang(index) .le. 0.D0) then ! Make sure it's from PHImax
+c     Try to increase PHImax
+        if (index .gt. 0) then
+          phi = res_ang(index-1)
+          diff = abs(res_ang(index) - res_ang(index-1))
+        endif
+c     Look for last (corresponding) x10
+        index = res_n - 1
+        do while ((.not.res_tab(1,1,index)) .or. res_tab(1,2,index))
+          index=index-1
+        enddo
+        if (index .lt. res_n-1) then
+          temp = abs(res_ang(index) - res_ang(index+1))
+          if (temp .lt. diff) then
+            phi = res_ang(index+1)
+            diff = temp
+          endif
+        endif
+      endif
+
+c     If increasing PHImax didn't work, decreasing PHImin
+c     will (with one exception)
+c     Look for first x10 (there MUST be at least one)
+      index = 0
+      do while ((.not.res_tab(1,1,index)) .or. res_tab(1,2,index))
+        index=index+1
+      enddo
+      if (res_ang(index) .lt. 0.D0) then ! Make sure it's from PHImin
+c     Try to decrease PHImin
+        if (index .lt. res_n-1) then
+          temp = abs(res_ang(index) - res_ang(index+1))
+          if (res_ang(index+1) .le. 0.D0 .and. temp .lt. diff) then
+            phi = res_ang(index+1)
+            diff = temp
+          endif
+        endif
+c     Look for last (corresponding) 01x
+        index = res_n - 1
+        do while (res_tab(1,0,index) .or. (.not.res_tab(1,1,index)))
+          index=index-1
+        enddo
+        if (index .gt. 0) then
+          temp = abs(res_ang(index) - res_ang(index-1))
+          if (res_ang(index-1) .ge. 0.D0 .and. temp .lt. diff) then
+            phi = res_ang(index-1)
+            diff = temp
+          endif
+        endif
+      endif
+
+c     If it still didn't work, it must be PHImax == 0. or PHImin == PI
+      if (diff .eq. flag) then
+        index = 0
+        if (res_tab(index,1,0) .or. (.not.res_tab(index,1,1)) .or.
+     +       res_tab(index,1,2)) index = res_n - 1
+c     This MUST work at this point
+        if (index .eq. 0) then
+          phi = res_ang(1)
+        else
+          phi = res_ang(index - 1)
+        endif
+      endif
+
+      return
+      end
+
+c-------------------------------------------------------------
+
+      integer function move_res(PHImin,PHImax,i_move)
+c     Moves residue i_move (in array c), leaving everything else fixed
+c     Starting geometry is not checked, it should be correct!
+c     R(,i_move) is the only residue that will move, but must have
+c     1 < i_move < nres (i.e., cannot move ends)
+c     Whether any output is done is controlled by locmove_output
+crc      implicit none
+
+c     Includes
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+      include 'COMMON.CHAIN'
+      include 'COMMON.GEO'
+      include 'COMMON.LOCMOVE'
+
+c     External functions
+      double precision ran_number
+      external ran_number
+
+c     Input arguments
+      double precision PHImin,PHImax
+      integer i_move
+
+c     RETURN VALUES:
+c     0: move successfull
+c     1: Dmin or Dmax had to be modified
+c     2: move failed - check your input geometry
+
+
+c     Local variables
+      double precision X(0:2),Y(0:2),Z(0:2),Orig(0:2)
+      double precision P(0:2)
+      logical no_moves,done
+      integer index,i,j
+      double precision phi,temp,radius
+      double precision phi_start(0:11), phi_end(0:11)
+      integer phi_n
+
+c     Set up the coordinate system
+      do i=0,2
+        Orig(i)=0.5*(c(i+1,i_move-1)+c(i+1,i_move+1)) ! Position of origin
+      enddo
+
+      do i=0,2
+        Z(i)=c(i+1,i_move+1)-c(i+1,i_move-1)
+      enddo
+      temp=sqrt(Z(0)*Z(0)+Z(1)*Z(1)+Z(2)*Z(2))
+      do i=0,2
+        Z(i)=Z(i)/temp
+      enddo
+
+      do i=0,2
+        X(i)=c(i+1,i_move)-Orig(i)
+      enddo
+c     radius is the radius of the circle on which c(,i_move) can move
+      radius=sqrt(X(0)*X(0)+X(1)*X(1)+X(2)*X(2))
+      do i=0,2
+        X(i)=X(i)/radius
+      enddo
+
+      Y(0)=Z(1)*X(2)-X(1)*Z(2)
+      Y(1)=X(0)*Z(2)-Z(0)*X(2)
+      Y(2)=Z(0)*X(1)-X(0)*Z(1)
+
+c     Calculate min, max angles coming from dmin, dmax to c(,i_move-2)
+      if (i_move.gt.2) then
+        do i=0,2
+          P(i)=c(i+1,i_move-2)-Orig(i)
+        enddo
+        call minmax_angles(P(0)*X(0)+P(1)*X(1)+P(2)*X(2),
+     +       P(0)*Y(0)+P(1)*Y(1)+P(2)*Y(2),
+     +       P(0)*Z(0)+P(1)*Z(1)+P(2)*Z(2),
+     +       radius,a_n,a_ang,a_tab)
+      else
+        a_n=0
+      endif
+
+c     Calculate min, max angles coming from dmin, dmax to c(,i_move+2)
+      if (i_move.lt.nres-2) then
+        do i=0,2
+          P(i)=c(i+1,i_move+2)-Orig(i)
+        enddo
+        call minmax_angles(P(0)*X(0)+P(1)*X(1)+P(2)*X(2),
+     +       P(0)*Y(0)+P(1)*Y(1)+P(2)*Y(2),
+     +       P(0)*Z(0)+P(1)*Z(1)+P(2)*Z(2),
+     +       radius,b_n,b_ang,b_tab)
+      else
+        b_n=0
+      endif
+
+c     Construct the resulting table for alpha and beta
+      call construct_tab()
+
+      if (locmove_output) then
+        print *,'ALPHAS & BETAS TABLE'
+        call output_tabs()
+      endif
+
+c     Check that there is at least one possible move
+      no_moves = .true.
+      if (res_n .eq. 0) then
+        no_moves = .false.
+      else
+        index = 0
+        do while ((index .lt. res_n) .and. no_moves)
+          if (res_tab(2,1,index)) no_moves = .false.
+          index=index+1
+        enddo
+      endif
+      if (no_moves) then
+        if (locmove_output) print *,'   ***   Cannot move anywhere'
+        move_res=2
+        return
+      endif
+
+c     Transfer res_... into a_...
+      a_n = 0
+      do i=0,res_n-1
+        if ( (res_tab(2,0,i).neqv.res_tab(2,1,i)) .or.
+     +       (res_tab(2,0,i).neqv.res_tab(2,2,i)) ) then
+          a_ang(a_n) = res_ang(i)
+          do j=0,2
+            a_tab(j,a_n) = res_tab(2,j,i)
+          enddo
+          a_n=a_n+1
+        endif
+      enddo
+
+c     Check that the PHI's are within [0,PI]
+      if (PHImin .lt. 0. .or. abs(PHImin) .lt. small) PHImin = -flag
+      if (PHImin .gt. PI .or. abs(PHImin-PI) .lt. small) PHImin = PI
+      if (PHImax .gt. PI .or. abs(PHImax-PI) .lt. small) PHImax = flag
+      if (PHImax .lt. 0. .or. abs(PHImax) .lt. small) PHImax = 0.
+      if (PHImax .lt. PHImin) PHImax = PHImin
+c     Calculate min and max angles coming from PHImin and PHImax,
+c     and put them in b_...
+      call angles2tab(PHImin, PHImax, b_n, b_ang, b_tab)
+c     Construct the final table
+      call construct_tab()
+
+      if (locmove_output) then
+        print *,'FINAL TABLE'
+        call output_tabs()
+      endif
+
+c     Check that there is at least one possible move
+      no_moves = .true.
+      if (res_n .eq. 0) then
+        no_moves = .false.
+      else
+        index = 0
+        do while ((index .lt. res_n) .and. no_moves)
+          if (res_tab(2,1,index)) no_moves = .false.
+          index=index+1
+        enddo
+      endif
+
+      if (no_moves) then
+c     Take care of the case where no solution exists...
+        call fix_no_moves(phi)
+        if (locmove_output) then
+          print *,'   ***   Had to modify PHImin or PHImax'
+          print *,'phi: ',phi*rad2deg
+        endif
+        move_res=1
+      else
+c     ...or calculate the solution
+c     Construct phi_start/phi_end arrays
+        call construct_ranges(phi_n, phi_start, phi_end)
+c     Choose random angle phi in allowed range(s)
+        temp = 0.
+        do i=0,phi_n-1
+          temp = temp + phi_end(i) - phi_start(i)
+        enddo
+        phi = ran_number(phi_start(0),phi_start(0)+temp)
+        index = 0
+        done = .false.
+        do while (.not.done)
+          if (phi .lt. phi_end(index)) then
+            done = .true.
+          else
+            index=index+1
+          endif
+          if (index .eq. phi_n) then
+            done = .true.
+          else if (.not.done) then
+            phi = phi + phi_start(index) - phi_end(index-1)
+          endif
+        enddo
+        if (index.eq.phi_n) phi=phi_end(phi_n-1) ! Fix numerical errors
+        if (phi .gt. PI) phi = phi-2.*PI
+
+        if (locmove_output) then
+          print *,'ALLOWED RANGE(S)'
+          do i=0,phi_n-1
+            print *,phi_start(i)*rad2deg,phi_end(i)*rad2deg
+          enddo
+          print *,'phi: ',phi*rad2deg
+        endif
+        move_res=0
+      endif
+
+c     Re-use radius as temp variable
+      temp=radius*cos(phi)
+      radius=radius*sin(phi)
+      do i=0,2
+        c(i+1,i_move)=Orig(i)+temp*X(i)+radius*Y(i)
+      enddo
+
+      return
+      end
+
+c-------------------------------------------------------------
+
+      subroutine loc_test
+crc      implicit none
+
+c     Includes
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+      include 'COMMON.GEO'
+      include 'COMMON.LOCAL'
+      include 'COMMON.LOCMOVE'
+
+c     External functions
+      integer move_res
+      external move_res
+
+c     Local variables
+      integer i,j
+      integer phi_n
+      double precision phi_start(0:11),phi_end(0:11)
+      double precision phi
+      double precision R(0:2,0:5)
+
+      locmove_output=.true.
+
+c      call angles2tab(30.*deg2rad,70.*deg2rad,a_n,a_ang,a_tab)
+c      call angles2tab(80.*deg2rad,130.*deg2rad,b_n,b_ang,b_tab)
+c      call minmax_angles(0.D0,3.8D0,0.D0,3.8D0,b_n,b_ang,b_tab)
+c      call construct_tab
+c      call output_tabs
+
+c      call construct_ranges(phi_n,phi_start,phi_end)
+c      do i=0,phi_n-1
+c        print *,phi_start(i)*rad2deg,phi_end(i)*rad2deg
+c      enddo
+
+c      call fix_no_moves(phi)
+c      print *,'NO MOVES FOUND, BEST PHI IS',phi*rad2deg
+
+      R(0,0)=0.D0
+      R(1,0)=0.D0
+      R(2,0)=0.D0
+      R(0,1)=0.D0
+      R(1,1)=-cos(28.D0*deg2rad)
+      R(2,1)=-0.5D0-sin(28.D0*deg2rad)
+      R(0,2)=0.D0
+      R(1,2)=0.D0
+      R(2,2)=-0.5D0
+      R(0,3)=cos(30.D0*deg2rad)
+      R(1,3)=0.D0
+      R(2,3)=0.D0
+      R(0,4)=0.D0
+      R(1,4)=0.D0
+      R(2,4)=0.5D0
+      R(0,5)=0.D0
+      R(1,5)=cos(26.D0*deg2rad)
+      R(2,5)=0.5D0+sin(26.D0*deg2rad)
+      do i=1,5
+        do j=0,2
+          R(j,i)=vbl*R(j,i)
+        enddo
+      enddo
+c      i=move_res(R(0,1),0.D0*deg2rad,180.D0*deg2rad)
+      imov=nnt
+      i=move_res(0.D0*deg2rad,180.D0*deg2rad,imov)
+      print *,'RETURNED ',i
+      print *,(R(i,3)/vbl,i=0,2)
+
+      return
+      end
+
+c-------------------------------------------------------------