source/unres/src_MD/src/local_move.f

   1 c-------------------------------------------------------------
   2
   3       subroutine local_move_init(debug)
   4 crc      implicit none
   5
   6 c     Includes
   7       implicit real*8 (a-h,o-z)
   8       include 'DIMENSIONS'  ! Needed by COMMON.LOCAL
   9       include 'COMMON.GEO'  ! For pi, deg2rad
  10       include 'COMMON.LOCAL'  ! For vbl
  11       include 'COMMON.LOCMOVE'
  12
  13 c     INPUT arguments
  14       logical debug
  15
  16
  17 c     Determine wheter to do some debugging output
  18       locmove_output=debug
  19
  20 c     Set the init_called flag to 1
  21       init_called=1
  22
  23 c     The following are never changed
  24       min_theta=60.D0*deg2rad  ! (0,PI)
  25       max_theta=175.D0*deg2rad  ! (0,PI)
  26       dmin2=vbl*vbl*2.*(1.-cos(min_theta))
  27       dmax2=vbl*vbl*2.*(1.-cos(max_theta))
  28       flag=1.0D300
  29       small=1.0D-5
  30       small2=0.5*small*small
  31
  32 c     Not really necessary...
  33       a_n=0
  34       b_n=0
  35       res_n=0
  36
  37       return
  38       end
  39
  40 c-------------------------------------------------------------
  41
  42       subroutine local_move(n_start, n_end, PHImin, PHImax)
  43 c     Perform a local move between residues m and n (inclusive)
  44 c     PHImin and PHImax [0,PI] determine the size of the move
  45 c     Works on whatever structure is in the variables theta and phi,
  46 c     sidechain variables are left untouched
  47 c     The final structure is NOT minimized, but both the cartesian
  48 c     variables c and the angles are up-to-date at the end (no further
  49 c     chainbuild is required)
  50 crc      implicit none
  51
  52 c     Includes
  53       implicit real*8 (a-h,o-z)
  54       include 'DIMENSIONS'
  55       include 'COMMON.GEO'
  56       include 'COMMON.CHAIN'
  57       include 'COMMON.VAR'
  58       include 'COMMON.MINIM'
  59       include 'COMMON.SBRIDGE'
  60       include 'COMMON.LOCMOVE'
  61
  62 c     External functions
  63       integer move_res
  64       external move_res
  65       double precision ran_number
  66       external ran_number
  67
  68 c     INPUT arguments
  69       integer n_start, n_end  ! First and last residues to move
  70       double precision PHImin, PHImax  ! min/max angles [0,PI]
  71
  72 c     Local variables
  73       integer i,j
  74       double precision min,max
  75       integer iretcode
  76
  77
  78 c     Check if local_move_init was called.  This assumes that it
  79 c     would not be 1 if not explicitely initialized
  80       if (init_called.ne.1) then
  81         write(6,*)'   ***   local_move_init not called!!!'
  82         stop
  83       endif
  84
  85 c     Quick check for crazy range
  86       if (n_start.gt.n_end .or. n_start.lt.1 .or. n_end.gt.nres) then
  87         write(6,'(a,i3,a,i3)')
  88      +       '   ***   Cannot make local move between n_start = ',
  89      +       n_start,' and n_end = ',n_end
  90         return
  91       endif
  92
  93 c     Take care of end residues first...
  94       if (n_start.eq.1) then
  95 c     Move residue 1 (completely random)
  96         theta(3)=ran_number(min_theta,max_theta)
  97         phi(4)=ran_number(-PI,PI)
  98         i=2
  99       else
 100         i=n_start
 101       endif
 102       if (n_end.eq.nres) then
 103 c     Move residue nres (completely random)
 104         theta(nres)=ran_number(min_theta,max_theta)
 105         phi(nres)=ran_number(-PI,PI)
 106         j=nres-1
 107       else
 108         j=n_end
 109       endif
 110
 111 c     ...then go through all other residues one by one
 112 c     Start from the two extremes and converge
 113       call chainbuild
 114       do while (i.le.j)
 115         min=PHImin
 116         max=PHImax
 117 c$$$c     Move the first two residues by less than the others
 118 c$$$        if (i-n_start.lt.3) then
 119 c$$$          if (i-n_start.eq.0) then
 120 c$$$            min=0.4*PHImin
 121 c$$$            max=0.4*PHImax
 122 c$$$          else if (i-n_start.eq.1) then
 123 c$$$            min=0.8*PHImin
 124 c$$$            max=0.8*PHImax
 125 c$$$          else if (i-n_start.eq.2) then
 126 c$$$            min=PHImin
 127 c$$$            max=PHImax
 128 c$$$          endif
 129 c$$$        endif
 130
 131 c     The actual move, on residue i
 132         iretcode=move_res(min,max,i,c)  ! Discard iretcode
 133         i=i+1
 134
 135         if (i.le.j) then
 136           min=PHImin
 137           max=PHImax
 138 c$$$c     Move the last two residues by less than the others
 139 c$$$          if (n_end-j.lt.3) then
 140 c$$$            if (n_end-j.eq.0) then
 141 c$$$              min=0.4*PHImin
 142 c$$$              max=0.4*PHImax
 143 c$$$            else if (n_end-j.eq.1) then
 144 c$$$              min=0.8*PHImin
 145 c$$$              max=0.8*PHImax
 146 c$$$            else if (n_end-j.eq.2) then
 147 c$$$              min=PHImin
 148 c$$$              max=PHImax
 149 c$$$            endif
 150 c$$$          endif
 151
 152 c     The actual move, on residue j
 153           iretcode=move_res(min,max,j,c)  ! Discard iretcode
 154           j=j-1
 155         endif
 156       enddo
 157
 158       call int_from_cart(.false.,.false.)
 159
 160       return
 161       end
 162
 163 c-------------------------------------------------------------
 164
 165       subroutine output_tabs
 166 c     Prints out the contents of a_..., b_..., res_...
 167       implicit none
 168
 169 c     Includes
 170       include 'COMMON.GEO'
 171       include 'COMMON.LOCMOVE'
 172
 173 c     Local variables
 174       integer i,j
 175
 176
 177       write(6,*)'a_...'
 178       write(6,'(8f7.1)')(a_ang(i)*rad2deg,i=0,a_n-1)
 179       write(6,'(8(2x,3l1,2x))')((a_tab(i,j),i=0,2),j=0,a_n-1)
 180
 181       write(6,*)'b_...'
 182       write(6,'(4f7.1)')(b_ang(i)*rad2deg,i=0,b_n-1)
 183       write(6,'(4(2x,3l1,2x))')((b_tab(i,j),i=0,2),j=0,b_n-1)
 184
 185       write(6,*)'res_...'
 186       write(6,'(12f7.1)')(res_ang(i)*rad2deg,i=0,res_n-1)
 187       write(6,'(12(2x,3l1,2x))')((res_tab(0,i,j),i=0,2),j=0,res_n-1)
 188       write(6,'(12(2x,3l1,2x))')((res_tab(1,i,j),i=0,2),j=0,res_n-1)
 189       write(6,'(12(2x,3l1,2x))')((res_tab(2,i,j),i=0,2),j=0,res_n-1)
 190
 191       return
 192       end
 193
 194 c-------------------------------------------------------------
 195
 196       subroutine angles2tab(PHImin,PHImax,n,ang,tab)
 197 c     Only uses angles if [0,PI] (but PHImin cannot be 0.,
 198 c     and PHImax cannot be PI)
 199       implicit none
 200
 201 c     Includes
 202       include 'COMMON.GEO'
 203
 204 c     INPUT arguments
 205       double precision PHImin,PHImax
 206
 207 c     OUTPUT arguments
 208       integer n
 209       double precision ang(0:3)
 210       logical tab(0:2,0:3)
 211
 212
 213       if (PHImin .eq. PHImax) then
 214 c     Special case with two 010's
 215         n = 2;
 216         ang(0) = -PHImin;
 217         ang(1) = PHImin;
 218         tab(0,0) = .false.
 219         tab(2,0) = .false.
 220         tab(0,1) = .false.
 221         tab(2,1) = .false.
 222         tab(1,0) = .true.
 223         tab(1,1) = .true.
 224       else if (PHImin .eq. PI) then
 225 c     Special case with one 010
 226         n = 1
 227         ang(0) = PI
 228         tab(0,0) = .false.
 229         tab(2,0) = .false.
 230         tab(1,0) = .true.
 231       else if (PHImax .eq. 0.) then
 232 c     Special case with one 010
 233         n = 1
 234         ang(0) = 0.
 235         tab(0,0) = .false.
 236         tab(2,0) = .false.
 237         tab(1,0) = .true.
 238       else
 239 c     Standard cases
 240         n = 0
 241         if (PHImin .gt. 0.) then
 242 c     Start of range (011)
 243           ang(n) = PHImin
 244           tab(0,n) = .false.
 245           tab(1,n) = .true.
 246           tab(2,n) = .true.
 247 c     End of range (110)
 248           ang(n+1) = -PHImin
 249           tab(0,n+1) = .true.
 250           tab(1,n+1) = .true.
 251           tab(2,n+1) = .false.
 252           n = n+2
 253         endif
 254         if (PHImax .lt. PI) then
 255 c     Start of range (011)
 256           ang(n) = -PHImax
 257           tab(0,n) = .false.
 258           tab(1,n) = .true.
 259           tab(2,n) = .true.
 260 c     End of range (110)
 261           ang(n+1) = PHImax
 262           tab(0,n+1) = .true.
 263           tab(1,n+1) = .true.
 264           tab(2,n+1) = .false.
 265           n = n+2
 266         endif
 267       endif
 268
 269       return
 270       end
 271
 272 c-------------------------------------------------------------
 273
 274       subroutine minmax_angles(x,y,z,r,n,ang,tab)
 275 c     When solutions do not exist, assume all angles
 276 c     are acceptable - i.e., initial geometry must be correct
 277       implicit none
 278
 279 c     Includes
 280       include 'COMMON.GEO'
 281       include 'COMMON.LOCMOVE'
 282
 283 c     Input arguments
 284       double precision x,y,z,r
 285
 286 c     Output arguments
 287       integer n
 288       double precision ang(0:3)
 289       logical tab(0:2,0:3)
 290
 291 c     Local variables
 292       double precision num, denom, phi
 293       double precision Kmin, Kmax
 294       integer i
 295
 296
 297       num = x*x + y*y + z*z
 298       denom = x*x + y*y
 299       n = 0
 300       if (denom .gt. 0.) then
 301         phi = atan2(y,x)
 302         denom = 2.*r*sqrt(denom)
 303         num = num+r*r
 304         Kmin = (num - dmin2)/denom
 305         Kmax = (num - dmax2)/denom
 306
 307 c     Allowed values of K (else all angles are acceptable)
 308 c     -1 <= Kmin <  1
 309 c     -1 <  Kmax <= 1
 310         if (Kmin .gt. 1. .or. abs(Kmin-1.) .lt. small2) then
 311           Kmin = -flag
 312         else if (Kmin .lt. -1. .or. abs(Kmin+1.) .lt. small2) then
 313           Kmin = PI
 314         else
 315           Kmin = acos(Kmin)
 316         endif
 317
 318         if (Kmax .lt. -1. .or. abs(Kmax+1.) .lt. small2) then
 319           Kmax = flag
 320         else if (Kmax .gt. 1. .or. abs(Kmax-1.) .lt. small2) then
 321           Kmax = 0.
 322         else
 323           Kmax = acos(Kmax)
 324         endif
 325
 326         if (Kmax .lt. Kmin) Kmax = Kmin
 327
 328         call angles2tab(Kmin, Kmax, n, ang, tab)
 329
 330 c     Add phi and check that angles are within range (-PI,PI]
 331         do i=0,n-1
 332           ang(i) = ang(i)+phi
 333           if (ang(i) .le. -PI) then
 334             ang(i) = ang(i)+2.*PI
 335           else if (ang(i) .gt. PI) then
 336             ang(i) = ang(i)-2.*PI
 337           endif
 338         enddo
 339       endif
 340
 341       return
 342       end
 343
 344 c-------------------------------------------------------------
 345
 346       subroutine construct_tab
 347 c     Take a_... and b_... values and produces the results res_...
 348 c     x_ang are assumed to be all different (diff > small)
 349 c     x_tab(1,i) must be 1 for all i (i.e., all x_ang are acceptable)
 350       implicit none
 351
 352 c     Includes
 353       include 'COMMON.LOCMOVE'
 354
 355 c     Local variables
 356       integer n_max,i,j,index
 357       logical done
 358       double precision phi
 359
 360
 361       n_max = a_n + b_n
 362       if (n_max .eq. 0) then
 363         res_n = 0
 364         return
 365       endif
 366
 367       do i=0,n_max-1
 368         do j=0,1
 369           res_tab(j,0,i) = .true.
 370           res_tab(j,2,i) = .true.
 371           res_tab(j,1,i) = .false.
 372         enddo
 373       enddo
 374
 375       index = 0
 376       phi = -flag
 377       done = .false.
 378       do while (.not.done)
 379         res_ang(index) = flag
 380
 381 c     Check a first...
 382         do i=0,a_n-1
 383           if ((a_ang(i)-phi).gt.small .and.
 384      +         a_ang(i) .lt. res_ang(index)) then
 385 c     Found a lower angle
 386             res_ang(index) = a_ang(i)
 387 c     Copy the values from a_tab into res_tab(0,,)
 388             res_tab(0,0,index) = a_tab(0,i)
 389             res_tab(0,1,index) = a_tab(1,i)
 390             res_tab(0,2,index) = a_tab(2,i)
 391 c     Set default values for res_tab(1,,)
 392             res_tab(1,0,index) = .true.
 393             res_tab(1,1,index) = .false.
 394             res_tab(1,2,index) = .true.
 395           else if (abs(a_ang(i)-res_ang(index)).lt.small) then
 396 c     Found an equal angle (can only be equal to a b_ang)
 397             res_tab(0,0,index) = a_tab(0,i)
 398             res_tab(0,1,index) = a_tab(1,i)
 399             res_tab(0,2,index) = a_tab(2,i)
 400           endif
 401         enddo
 402 c     ...then check b
 403         do i=0,b_n-1
 404           if ((b_ang(i)-phi).gt.small .and.
 405      +         b_ang(i) .lt. res_ang(index)) then
 406 c     Found a lower angle
 407             res_ang(index) = b_ang(i)
 408 c     Copy the values from b_tab into res_tab(1,,)
 409             res_tab(1,0,index) = b_tab(0,i)
 410             res_tab(1,1,index) = b_tab(1,i)
 411             res_tab(1,2,index) = b_tab(2,i)
 412 c     Set default values for res_tab(0,,)
 413             res_tab(0,0,index) = .true.
 414             res_tab(0,1,index) = .false.
 415             res_tab(0,2,index) = .true.
 416           else if (abs(b_ang(i)-res_ang(index)).lt.small) then
 417 c     Found an equal angle (can only be equal to an a_ang)
 418             res_tab(1,0,index) = b_tab(0,i)
 419             res_tab(1,1,index) = b_tab(1,i)
 420             res_tab(1,2,index) = b_tab(2,i)
 421           endif
 422         enddo
 423
 424         if (res_ang(index) .eq. flag) then
 425           res_n = index
 426           done = .true.
 427         else if (index .eq. n_max-1) then
 428           res_n = n_max
 429           done = .true.
 430         else
 431           phi = res_ang(index)  ! Store previous angle
 432           index = index+1
 433         endif
 434       enddo
 435
 436 c     Fill the gaps
 437 c     First a...
 438       index = 0
 439       if (a_n .gt. 0) then
 440         do while (.not.res_tab(0,1,index))
 441           index=index+1
 442         enddo
 443         done = res_tab(0,2,index)
 444         do i=index+1,res_n-1
 445           if (res_tab(0,1,i)) then
 446             done = res_tab(0,2,i)
 447           else
 448             res_tab(0,0,i) = done
 449             res_tab(0,1,i) = done
 450             res_tab(0,2,i) = done
 451           endif
 452         enddo
 453         done = res_tab(0,0,index)
 454         do i=index-1,0,-1
 455           if (res_tab(0,1,i)) then
 456             done = res_tab(0,0,i)
 457           else
 458             res_tab(0,0,i) = done
 459             res_tab(0,1,i) = done
 460             res_tab(0,2,i) = done
 461           endif
 462         enddo
 463       else
 464         do i=0,res_n-1
 465           res_tab(0,0,i) = .true.
 466           res_tab(0,1,i) = .true.
 467           res_tab(0,2,i) = .true.
 468         enddo
 469       endif
 470 c     ...then b
 471       index = 0
 472       if (b_n .gt. 0) then
 473         do while (.not.res_tab(1,1,index))
 474           index=index+1
 475         enddo
 476         done = res_tab(1,2,index)
 477         do i=index+1,res_n-1
 478           if (res_tab(1,1,i)) then
 479             done = res_tab(1,2,i)
 480           else
 481             res_tab(1,0,i) = done
 482             res_tab(1,1,i) = done
 483             res_tab(1,2,i) = done
 484           endif
 485         enddo
 486         done = res_tab(1,0,index)
 487         do i=index-1,0,-1
 488           if (res_tab(1,1,i)) then
 489             done = res_tab(1,0,i)
 490           else
 491             res_tab(1,0,i) = done
 492             res_tab(1,1,i) = done
 493             res_tab(1,2,i) = done
 494           endif
 495         enddo
 496       else
 497         do i=0,res_n-1
 498           res_tab(1,0,i) = .true.
 499           res_tab(1,1,i) = .true.
 500           res_tab(1,2,i) = .true.
 501         enddo
 502       endif
 503
 504 c     Finally fill the last row with AND operation
 505       do i=0,res_n-1
 506         do j=0,2
 507           res_tab(2,j,i) = (res_tab(0,j,i) .and. res_tab(1,j,i))
 508         enddo
 509       enddo
 510
 511       return
 512       end
 513
 514 c-------------------------------------------------------------
 515
 516       subroutine construct_ranges(phi_n,phi_start,phi_end)
 517 c     Given the data in res_..., construct a table of
 518 c     min/max allowed angles
 519       implicit none
 520
 521 c     Includes
 522       include 'COMMON.GEO'
 523       include 'COMMON.LOCMOVE'
 524
 525 c     Output arguments
 526       integer phi_n
 527       double precision phi_start(0:11),phi_end(0:11)
 528
 529 c     Local variables
 530       logical done
 531       integer index
 532
 533
 534       if (res_n .eq. 0) then
 535 c     Any move is allowed
 536         phi_n = 1
 537         phi_start(0) = -PI
 538         phi_end(0) = PI
 539       else
 540         phi_n = 0
 541         index = 0
 542         done = .false.
 543         do while (.not.done)
 544 c     Find start of range (01x)
 545           done = .false.
 546           do while (.not.done)
 547             if (res_tab(2,0,index).or.(.not.res_tab(2,1,index))) then
 548               index=index+1
 549             else
 550               done = .true.
 551               phi_start(phi_n) = res_ang(index)
 552             endif
 553             if (index .eq. res_n) done = .true.
 554           enddo
 555 c     If a start was found (index < res_n), find the end of range (x10)
 556 c     It may not be found without wrapping around
 557           if (index .lt. res_n) then
 558             done = .false.
 559             do while (.not.done)
 560               if ((.not.res_tab(2,1,index)).or.res_tab(2,2,index)) then
 561                 index=index+1
 562               else
 563                 done = .true.
 564               endif
 565               if (index .eq. res_n) done = .true.
 566             enddo
 567             if (index .lt. res_n) then
 568 c     Found the end of the range
 569               phi_end(phi_n) = res_ang(index)
 570               phi_n=phi_n+1
 571               index=index+1
 572               if (index .eq. res_n) then
 573                 done = .true.
 574               else
 575                 done = .false.
 576               endif
 577             else
 578 c     Need to wrap around
 579               done = .true.
 580               phi_end(phi_n) = flag
 581             endif
 582           endif
 583         enddo
 584 c     Take care of the last one if need to wrap around
 585         if (phi_end(phi_n) .eq. flag) then
 586           index = 0
 587           do while ((.not.res_tab(2,1,index)).or.res_tab(2,2,index))
 588             index=index+1
 589           enddo
 590           phi_end(phi_n) = res_ang(index) + 2.*PI
 591           phi_n=phi_n+1
 592         endif
 593       endif
 594
 595       return
 596       end
 597
 598 c-------------------------------------------------------------
 599
 600       subroutine fix_no_moves(phi)
 601       implicit none
 602
 603 c     Includes
 604       include 'COMMON.GEO'
 605       include 'COMMON.LOCMOVE'
 606
 607 c     Output arguments
 608       double precision phi
 609
 610 c     Local variables
 611       integer index
 612       double precision diff,temp
 613
 614
 615 c     Look for first 01x in gammas (there MUST be at least one)
 616       diff = flag
 617       index = 0
 618       do while (res_tab(1,0,index) .or. (.not.res_tab(1,1,index)))
 619         index=index+1
 620       enddo
 621       if (res_ang(index) .le. 0.D0) then ! Make sure it's from PHImax
 622 c     Try to increase PHImax
 623         if (index .gt. 0) then
 624           phi = res_ang(index-1)
 625           diff = abs(res_ang(index) - res_ang(index-1))
 626         endif
 627 c     Look for last (corresponding) x10
 628         index = res_n - 1
 629         do while ((.not.res_tab(1,1,index)) .or. res_tab(1,2,index))
 630           index=index-1
 631         enddo
 632         if (index .lt. res_n-1) then
 633           temp = abs(res_ang(index) - res_ang(index+1))
 634           if (temp .lt. diff) then
 635             phi = res_ang(index+1)
 636             diff = temp
 637           endif
 638         endif
 639       endif
 640
 641 c     If increasing PHImax didn't work, decreasing PHImin
 642 c     will (with one exception)
 643 c     Look for first x10 (there MUST be at least one)
 644       index = 0
 645       do while ((.not.res_tab(1,1,index)) .or. res_tab(1,2,index))
 646         index=index+1
 647       enddo
 648       if (res_ang(index) .lt. 0.D0) then ! Make sure it's from PHImin
 649 c     Try to decrease PHImin
 650         if (index .lt. res_n-1) then
 651           temp = abs(res_ang(index) - res_ang(index+1))
 652           if (res_ang(index+1) .le. 0.D0 .and. temp .lt. diff) then
 653             phi = res_ang(index+1)
 654             diff = temp
 655           endif
 656         endif
 657 c     Look for last (corresponding) 01x
 658         index = res_n - 1
 659         do while (res_tab(1,0,index) .or. (.not.res_tab(1,1,index)))
 660           index=index-1
 661         enddo
 662         if (index .gt. 0) then
 663           temp = abs(res_ang(index) - res_ang(index-1))
 664           if (res_ang(index-1) .ge. 0.D0 .and. temp .lt. diff) then
 665             phi = res_ang(index-1)
 666             diff = temp
 667           endif
 668         endif
 669       endif
 670
 671 c     If it still didn't work, it must be PHImax == 0. or PHImin == PI
 672       if (diff .eq. flag) then
 673         index = 0
 674         if (res_tab(index,1,0) .or. (.not.res_tab(index,1,1)) .or.
 675      +       res_tab(index,1,2)) index = res_n - 1
 676 c     This MUST work at this point
 677         if (index .eq. 0) then
 678           phi = res_ang(1)
 679         else
 680           phi = res_ang(index - 1)
 681         endif
 682       endif
 683
 684       return
 685       end
 686
 687 c-------------------------------------------------------------
 688
 689       integer function move_res(PHImin,PHImax,i_move)
 690 c     Moves residue i_move (in array c), leaving everything else fixed
 691 c     Starting geometry is not checked, it should be correct!
 692 c     R(,i_move) is the only residue that will move, but must have
 693 c     1 < i_move < nres (i.e., cannot move ends)
 694 c     Whether any output is done is controlled by locmove_output
 695 crc      implicit none
 696
 697 c     Includes
 698       implicit real*8 (a-h,o-z)
 699       include 'DIMENSIONS'
 700       include 'COMMON.CHAIN'
 701       include 'COMMON.GEO'
 702       include 'COMMON.LOCMOVE'
 703
 704 c     External functions
 705       double precision ran_number
 706       external ran_number
 707
 708 c     Input arguments
 709       double precision PHImin,PHImax
 710       integer i_move
 711
 712 c     RETURN VALUES:
 713 c     0: move successfull
 714 c     1: Dmin or Dmax had to be modified
 715 c     2: move failed - check your input geometry
 716
 717
 718 c     Local variables
 719       double precision X(0:2),Y(0:2),Z(0:2),Orig(0:2)
 720       double precision P(0:2)
 721       logical no_moves,done
 722       integer index,i,j
 723       double precision phi,temp,radius
 724       double precision phi_start(0:11), phi_end(0:11)
 725       integer phi_n
 726
 727 c     Set up the coordinate system
 728       do i=0,2
 729         Orig(i)=0.5*(c(i+1,i_move-1)+c(i+1,i_move+1)) ! Position of origin
 730       enddo
 731
 732       do i=0,2
 733         Z(i)=c(i+1,i_move+1)-c(i+1,i_move-1)
 734       enddo
 735       temp=sqrt(Z(0)*Z(0)+Z(1)*Z(1)+Z(2)*Z(2))
 736       do i=0,2
 737         Z(i)=Z(i)/temp
 738       enddo
 739
 740       do i=0,2
 741         X(i)=c(i+1,i_move)-Orig(i)
 742       enddo
 743 c     radius is the radius of the circle on which c(,i_move) can move
 744       radius=sqrt(X(0)*X(0)+X(1)*X(1)+X(2)*X(2))
 745       do i=0,2
 746         X(i)=X(i)/radius
 747       enddo
 748
 749       Y(0)=Z(1)*X(2)-X(1)*Z(2)
 750       Y(1)=X(0)*Z(2)-Z(0)*X(2)
 751       Y(2)=Z(0)*X(1)-X(0)*Z(1)
 752
 753 c     Calculate min, max angles coming from dmin, dmax to c(,i_move-2)
 754       if (i_move.gt.2) then
 755         do i=0,2
 756           P(i)=c(i+1,i_move-2)-Orig(i)
 757         enddo
 758         call minmax_angles(P(0)*X(0)+P(1)*X(1)+P(2)*X(2),
 759      +       P(0)*Y(0)+P(1)*Y(1)+P(2)*Y(2),
 760      +       P(0)*Z(0)+P(1)*Z(1)+P(2)*Z(2),
 761      +       radius,a_n,a_ang,a_tab)
 762       else
 763         a_n=0
 764       endif
 765
 766 c     Calculate min, max angles coming from dmin, dmax to c(,i_move+2)
 767       if (i_move.lt.nres-2) then
 768         do i=0,2
 769           P(i)=c(i+1,i_move+2)-Orig(i)
 770         enddo
 771         call minmax_angles(P(0)*X(0)+P(1)*X(1)+P(2)*X(2),
 772      +       P(0)*Y(0)+P(1)*Y(1)+P(2)*Y(2),
 773      +       P(0)*Z(0)+P(1)*Z(1)+P(2)*Z(2),
 774      +       radius,b_n,b_ang,b_tab)
 775       else
 776         b_n=0
 777       endif
 778
 779 c     Construct the resulting table for alpha and beta
 780       call construct_tab()
 781
 782       if (locmove_output) then
 783         print *,'ALPHAS & BETAS TABLE'
 784         call output_tabs()
 785       endif
 786
 787 c     Check that there is at least one possible move
 788       no_moves = .true.
 789       if (res_n .eq. 0) then
 790         no_moves = .false.
 791       else
 792         index = 0
 793         do while ((index .lt. res_n) .and. no_moves)
 794           if (res_tab(2,1,index)) no_moves = .false.
 795           index=index+1
 796         enddo
 797       endif
 798       if (no_moves) then
 799         if (locmove_output) print *,'   ***   Cannot move anywhere'
 800         move_res=2
 801         return
 802       endif
 803
 804 c     Transfer res_... into a_...
 805       a_n = 0
 806       do i=0,res_n-1
 807         if ( (res_tab(2,0,i).neqv.res_tab(2,1,i)) .or.
 808      +       (res_tab(2,0,i).neqv.res_tab(2,2,i)) ) then
 809           a_ang(a_n) = res_ang(i)
 810           do j=0,2
 811             a_tab(j,a_n) = res_tab(2,j,i)
 812           enddo
 813           a_n=a_n+1
 814         endif
 815       enddo
 816
 817 c     Check that the PHI's are within [0,PI]
 818       if (PHImin .lt. 0. .or. abs(PHImin) .lt. small) PHImin = -flag
 819       if (PHImin .gt. PI .or. abs(PHImin-PI) .lt. small) PHImin = PI
 820       if (PHImax .gt. PI .or. abs(PHImax-PI) .lt. small) PHImax = flag
 821       if (PHImax .lt. 0. .or. abs(PHImax) .lt. small) PHImax = 0.
 822       if (PHImax .lt. PHImin) PHImax = PHImin
 823 c     Calculate min and max angles coming from PHImin and PHImax,
 824 c     and put them in b_...
 825       call angles2tab(PHImin, PHImax, b_n, b_ang, b_tab)
 826 c     Construct the final table
 827       call construct_tab()
 828
 829       if (locmove_output) then
 830         print *,'FINAL TABLE'
 831         call output_tabs()
 832       endif
 833
 834 c     Check that there is at least one possible move
 835       no_moves = .true.
 836       if (res_n .eq. 0) then
 837         no_moves = .false.
 838       else
 839         index = 0
 840         do while ((index .lt. res_n) .and. no_moves)
 841           if (res_tab(2,1,index)) no_moves = .false.
 842           index=index+1
 843         enddo
 844       endif
 845
 846       if (no_moves) then
 847 c     Take care of the case where no solution exists...
 848         call fix_no_moves(phi)
 849         if (locmove_output) then
 850           print *,'   ***   Had to modify PHImin or PHImax'
 851           print *,'phi: ',phi*rad2deg
 852         endif
 853         move_res=1
 854       else
 855 c     ...or calculate the solution
 856 c     Construct phi_start/phi_end arrays
 857         call construct_ranges(phi_n, phi_start, phi_end)
 858 c     Choose random angle phi in allowed range(s)
 859         temp = 0.
 860         do i=0,phi_n-1
 861           temp = temp + phi_end(i) - phi_start(i)
 862         enddo
 863         phi = ran_number(phi_start(0),phi_start(0)+temp)
 864         index = 0
 865         done = .false.
 866         do while (.not.done)
 867           if (phi .lt. phi_end(index)) then
 868             done = .true.
 869           else
 870             index=index+1
 871           endif
 872           if (index .eq. phi_n) then
 873             done = .true.
 874           else if (.not.done) then
 875             phi = phi + phi_start(index) - phi_end(index-1)
 876           endif
 877         enddo
 878         if (index.eq.phi_n) phi=phi_end(phi_n-1) ! Fix numerical errors
 879         if (phi .gt. PI) phi = phi-2.*PI
 880
 881         if (locmove_output) then
 882           print *,'ALLOWED RANGE(S)'
 883           do i=0,phi_n-1
 884             print *,phi_start(i)*rad2deg,phi_end(i)*rad2deg
 885           enddo
 886           print *,'phi: ',phi*rad2deg
 887         endif
 888         move_res=0
 889       endif
 890
 891 c     Re-use radius as temp variable
 892       temp=radius*cos(phi)
 893       radius=radius*sin(phi)
 894       do i=0,2
 895         c(i+1,i_move)=Orig(i)+temp*X(i)+radius*Y(i)
 896       enddo
 897
 898       return
 899       end
 900
 901 c-------------------------------------------------------------
 902
 903       subroutine loc_test
 904 crc      implicit none
 905
 906 c     Includes
 907       implicit real*8 (a-h,o-z)
 908       include 'DIMENSIONS'
 909       include 'COMMON.GEO'
 910       include 'COMMON.LOCAL'
 911       include 'COMMON.LOCMOVE'
 912
 913 c     External functions
 914       integer move_res
 915       external move_res
 916
 917 c     Local variables
 918       integer i,j
 919       integer phi_n
 920       double precision phi_start(0:11),phi_end(0:11)
 921       double precision phi
 922       double precision R(0:2,0:5)
 923
 924       locmove_output=.true.
 925
 926 c      call angles2tab(30.*deg2rad,70.*deg2rad,a_n,a_ang,a_tab)
 927 c      call angles2tab(80.*deg2rad,130.*deg2rad,b_n,b_ang,b_tab)
 928 c      call minmax_angles(0.D0,3.8D0,0.D0,3.8D0,b_n,b_ang,b_tab)
 929 c      call construct_tab
 930 c      call output_tabs
 931
 932 c      call construct_ranges(phi_n,phi_start,phi_end)
 933 c      do i=0,phi_n-1
 934 c        print *,phi_start(i)*rad2deg,phi_end(i)*rad2deg
 935 c      enddo
 936
 937 c      call fix_no_moves(phi)
 938 c      print *,'NO MOVES FOUND, BEST PHI IS',phi*rad2deg
 939
 940       R(0,0)=0.D0
 941       R(1,0)=0.D0
 942       R(2,0)=0.D0
 943       R(0,1)=0.D0
 944       R(1,1)=-cos(28.D0*deg2rad)
 945       R(2,1)=-0.5D0-sin(28.D0*deg2rad)
 946       R(0,2)=0.D0
 947       R(1,2)=0.D0
 948       R(2,2)=-0.5D0
 949       R(0,3)=cos(30.D0*deg2rad)
 950       R(1,3)=0.D0
 951       R(2,3)=0.D0
 952       R(0,4)=0.D0
 953       R(1,4)=0.D0
 954       R(2,4)=0.5D0
 955       R(0,5)=0.D0
 956       R(1,5)=cos(26.D0*deg2rad)
 957       R(2,5)=0.5D0+sin(26.D0*deg2rad)
 958       do i=1,5
 959         do j=0,2
 960           R(j,i)=vbl*R(j,i)
 961         enddo
 962       enddo
 963       i=move_res(R(0,1),0.D0*deg2rad,180.D0*deg2rad)
 964       print *,'RETURNED ',i
 965       print *,(R(i,3)/vbl,i=0,2)
 966
 967       return
 968       end
 969
 970 c-------------------------------------------------------------