--- /dev/null
+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-------------------------------------------------------------
projects / unres.git / blobdiff