Zmina 21 na ntyp1 w wham src i src-M
[unres.git] / source / unres / src_CSA / 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-------------------------------------------------------------