1      subroutine ana_groups(card,sgmnam,imol,isel,wt,x)
2c
3c $Id$
4c
5      implicit none
6c
7#include "ana_common.fh"
8#include "global.fh"
9#include "mafdecls.fh"
10#include "msgids.fh"
11#include "rtdb.fh"
12c
13      character*80 card
14      character*16 sgmnam(nsa)
15      real*8 wt(nsa)
16      integer isel(nsa),imol(msa)
17      real*8 x(nsa,3)
18c
19      integer i
20c
21      ngroups=ngroups+1
22      if(ngroups.gt.maxgrp) call md_abort('Increase maxgrp',maxgrp)
23c
24      if(card(1:6).eq.'groups') then
25      read(card(8:67),1001) (igroups(ngroups,i),i=1,6),
26     + (rgroups(ngroups,i),i=1,2)
27 1001 format(4i7,i5,i3,2f12.6)
28      endif
29      print*,'++++++++++',igroups(ngroups,5)
30c
31      filgrp=card(68:80)
32      if(filgrp(1:1).ne.' ') then
33      open(unit=lfngrp,file=filgrp(1:index(filgrp,' ')-1),
34     + form='formatted',status='unknown')
35      if(igroups(ngroups,5).eq.1)
36     + call ana_grpcogdis(sgmnam,imol,isel,wt,x,ngroups)
37      if(igroups(ngroups,5).eq.2)
38     + call ana_grpdis(sgmnam,imol,isel,wt,x,ngroups)
39      if(igroups(ngroups,5).eq.4)
40     + call ana_grpcogang(sgmnam,imol,isel,wt,x,ngroups)
41      if(igroups(ngroups,5).eq.5)
42     + call ana_grpvectors(sgmnam,imol,isel,wt,x,ngroups)
43      ngroups=ngroups-1
44      close(unit=lfngrp)
45      endif
46c
47      return
48      end
49      subroutine ana_group(card,sgmnam,imol,isel,wt,x,iwrk)
50c
51c $Id$
52c
53      implicit none
54c
55#include "ana_common.fh"
56#include "global.fh"
57#include "mafdecls.fh"
58#include "msgids.fh"
59#include "rtdb.fh"
60c
61      character*80 card
62      character*16 sgmnam(nsa)
63      real*8 wt(nsa)
64      integer isel(nsa),imol(msa),iwrk(mxdef,mxnum,maxgrp)
65      real*8 x(nsa,3)
66c
67      integer i,j
68c
69      ngroup=ngroup+1
70      if(ngroup.gt.maxgrp) call md_abort('Increase maxgrp',maxgrp)
71c
72      read(card(8:46),1000) (igroup(ngroup,i),i=1,3),
73     + (rgroup(ngroup,i),i=1,2)
74 1000 format(i7,i5,i3,2f12.6)
75c
76      do 1 i=1,mxdef
77      do 2 j=1,mxnum
78      iwrk(i,j,ngroup)=0
79    2 continue
80    1 continue
81c
82      return
83      end
84      subroutine ana_grpdis(sgmnam,imol,isel,wt,x,igr)
85c
86      implicit none
87c
88#include "ana_common.fh"
89#include "global.fh"
90#include "mafdecls.fh"
91#include "msgids.fh"
92c
93      character*16 sgmnam(nsa)
94      real*8 wt(nsa)
95      integer isel(nsa),igr,imol(msa)
96      real*8 x(nsa,3)
97c
98      integer igrp,jgrp
99      integer i,j,jfrom,ia,ja,number
100      real*8 dist,dx,dy,dz
101c
102      igrp=igroups(igr,1)
103      jgrp=igroups(igr,2)
104c
105      number=0
106c
107      if(ldef(igrp).lt.0) return
108c
109      write(lfngrp,1000) igr,igroups(igr,5)
110 1000 format(2i5)
111c
112      ito=ldef(igrp)
113      if(igrp.eq.jgrp) ito=ito-1
114      do 1 i=1,ito
115      ia=idef(igrp,i)
116      jfrom=1
117      if(igrp.eq.jgrp) jfrom=i+1
118      do 2 j=jfrom,ldef(jgrp)
119      ja=idef(jgrp,j)
120      dx=abs(x(ia,1)-x(ja,1))
121      dy=abs(x(ia,2)-x(ja,2))
122      dz=abs(x(ia,3)-x(ja,3))
123      if(igroups(igr,6).eq.1) then
124      if(dz.gt.box(3)) dz=dz-box(3)
125      elseif(igroups(igr,6).eq.2) then
126      if(dx.gt.box(1)) dx=dx-box(1)
127      if(dy.gt.box(2)) dy=dy-box(2)
128      elseif(igroups(igr,6).eq.3) then
129      if(dx.gt.box(1)) dx=dx-box(1)
130      if(dy.gt.box(2)) dy=dy-box(2)
131      if(dz.gt.box(3)) dz=dz-box(3)
132      endif
133      dist=sqrt(dx*dx+dy*dy+dz*dz)
134      if(dist.ge.rgroups(igr,1).and.dist.le.rgroups(igr,2)) then
135      write(lfngrp,1001)
136     + imol(ia),sgmnam(ia)(11:16),sgmnam(ia)(1:5),sgmnam(ia)(6:10),
137     + imol(ja),sgmnam(ja)(11:16),sgmnam(ja)(1:5),sgmnam(ja)(6:10),dist
138 1001 format(2(i5,a6,' ',a5,':',a5,' '),f12.6)
139      number=number+1
140      endif
141    2 continue
142    1 continue
143c
144      write(lfngrp,1002) 0
145 1002 format(i5)
146c
147      return
148      end
149      subroutine ana_grpvectors(sgmnam,imol,isel,wt,x,igr)
150c
151      implicit none
152c
153#include "ana_common.fh"
154#include "global.fh"
155#include "mafdecls.fh"
156#include "msgids.fh"
157c
158      character*16 sgmnam(nsa)
159      real*8 wt(nsa)
160      integer isel(nsa),igr,imol(msa)
161      real*8 x(nsa,3),vi(3),vj(3)
162c
163      integer igrp,jgrp
164      integer i,j,jfrom,ia,ja,number,nsum
165      real*8 dist,dx,dy,dz,dot,dotsum
166c
167      igrp=igroups(igr,1)
168      jgrp=igroups(igr,2)
169      print*,'igrp,jgrp=',igrp,jgrp
170c
171      number=0
172c
173      if(ldef(igrp).lt.0) return
174c
175      write(lfngrp,1000) igr,igroups(igr,5)
176 1000 format(2i5)
177c
178      if(igrp.eq.jgrp)
179     + call md_abort('vectors: single atom list',0)
180      if(ldef(igrp).ne.ldef(jgrp))
181     + call md_abort('vectors: unequal atom lists',0)
182      ito=ldef(igrp)
183      do 1 i=1,ito
184      ia=idef(igrp,i)
185      ja=idef(jgrp,i)
186      dx=abs(x(ia,1)-x(ja,1))
187      dy=abs(x(ia,2)-x(ja,2))
188      dz=abs(x(ia,3)-x(ja,3))
189      if(igroups(igr,6).eq.1) then
190      if(dz.gt.box(3)) dz=dz-box(3)
191      elseif(igroups(igr,6).eq.2) then
192      if(dx.gt.box(1)) dx=dx-box(1)
193      if(dy.gt.box(2)) dy=dy-box(2)
194      elseif(igroups(igr,6).eq.3) then
195      if(dx.gt.box(1)) dx=dx-box(1)
196      if(dy.gt.box(2)) dy=dy-box(2)
197      if(dz.gt.box(3)) dz=dz-box(3)
198      endif
199      dist=sqrt(dx*dx+dy*dy+dz*dz)
200c      write(lfngrp,1001)
201c     + imol(ia),sgmnam(ia)(11:16),sgmnam(ia)(1:5),sgmnam(ia)(6:10),
202c     + imol(ja),sgmnam(ja)(11:16),sgmnam(ja)(1:5),sgmnam(ja)(6:10),dist
203c 1001 format(2(i5,a6,' ',a5,':',a5,' '),f12.6)
204      write(lfngrp,1003) i,(x(ia,j),j=1,3),(x(ja,j)-x(ia,j),j=1,3),dist
205 1003 format(i5,7f12.6)
206      number=number+1
207    1 continue
208      dotsum=0.0d0
209      nsum=0
210      do 2 i=1,ito
211      ia=idef(igrp,i)
212      ja=idef(jgrp,i)
213      vi(1)=x(ia,1)-x(ja,1)
214      vi(2)=x(ia,2)-x(ja,2)
215      vi(3)=x(ia,3)-x(ja,3)
216      do 3 j=i+1,ito
217      ia=idef(igrp,j)
218      ja=idef(jgrp,j)
219      vj(1)=x(ia,1)-x(ja,1)
220      vj(2)=x(ia,2)-x(ja,2)
221      vj(3)=x(ia,3)-x(ja,3)
222      dot=(vi(1)*vj(1)+vi(2)*vj(2)+vi(3)*vj(3))/
223     + (sqrt(vi(1)**2+vi(2)**2+vi(3)**2)*
224     +  sqrt(vj(1)**2+vj(2)**2+vj(3)**2))
225c      write(lfngrp,1004) i,j,dot
226c 1004 format(2i5,f12.6)
227      dotsum=dotsum+dot
228      nsum=nsum+1
229    3 continue
230    2 continue
231      write(lfngrp,1005) time,dotsum/dble(nsum)
232 1005 format('da ',2f12.6)
233c    1 continue
234c
235      write(lfngrp,1002) 0
236 1002 format(i5)
237c
238      return
239      end
240      subroutine ana_lochdr(sgmnam,imol,isel)
241c
242      implicit none
243c
244#include "ana_common.fh"
245#include "global.fh"
246#include "mafdecls.fh"
247#include "msgids.fh"
248c
249      character*16 sgmnam(nsa)
250      integer imol(msa),isel(nsa)
251c
252      integer i,num,igrp,ia,j
253c
254      num=0
255      do 1 i=1,nsa
256      if(isel(i).ne.0) num=num+1
257    1 continue
258c
259      write(lfnloc,1000) ngroup
260 1000 format(2i7)
261      do 4 j=1,ngroup
262      igrp=igroup(j,1)
263      ito=ldef(igrp)
264      write(lfnloc,1000) ito,nsa
265      do 2 i=1,ito
266      ia=idef(igrp,i)
267      write(lfnloc,1001) ia,imol(ia),sgmnam(ia)(11:16),sgmnam(ia)(1:5),
268     + sgmnam(ia)(6:10)
269    2 continue
270    4 continue
271      write(lfnloc,1000) num,nsa
272      do 3 i=1,nsa
273      if(isel(i).ne.0) then
274      write(lfnloc,1001) i,imol(i),sgmnam(i)(11:16),sgmnam(i)(1:5),
275     + sgmnam(i)(6:10)
276 1001 format(i7,i5,a6,' ',a5,':',a5)
277      endif
278    3 continue
279c
280      return
281      end
282      subroutine ana_grploc(sgmnam,imol,isel,wt,x,igr,iwrk)
283c
284      implicit none
285c
286#include "ana_common.fh"
287#include "global.fh"
288#include "mafdecls.fh"
289#include "msgids.fh"
290c
291      character*16 sgmnam(nsa)
292      real*8 wt(nsa)
293      integer isel(nsa),igr,imol(msa),iwrk(mxdef,mxnum,maxgrp)
294      real*8 x(nsa,3)
295c
296      integer i,j,igrp,ia,ja,k,ntmp
297      real*8 dx,dy,dz,dist
298c
299      integer num,ndx(100)
300      real*8 rd(100)
301c
302      igrp=igroup(igr,1)
303      ito=ldef(igrp)
304c
305c      write(lfngrp,1000) time
306c 1000 format('Atom distances for selected atoms at time ',f12.6)
307c
308      do 1 i=1,ito
309      ia=idef(igrp,i)
310      num=0
311      do 2 ja=1,nsa
312      if(ia.ne.ja.and.isel(ja).ne.0) then
313      dx=abs(x(ia,1)-x(ja,1))
314      dy=abs(x(ia,2)-x(ja,2))
315      dz=abs(x(ia,3)-x(ja,3))
316      if(igroup(igr,3).eq.1) then
317      if(dz.gt.box(3)) dz=dz-box(3)
318      elseif(igroup(igr,3).eq.2) then
319      if(dx.gt.box(1)) dx=dx-box(1)
320      if(dy.gt.box(2)) dy=dy-box(2)
321      elseif(igroup(igr,3).eq.3) then
322      if(dx.gt.box(1)) dx=dx-box(1)
323      if(dy.gt.box(2)) dy=dy-box(2)
324      if(dz.gt.box(3)) dz=dz-box(3)
325      endif
326      dist=sqrt(dx*dx+dy*dy+dz*dz)
327      if(dist.ge.rgroup(igr,1).and.dist.le.rgroup(igr,2)) then
328      if(num.lt.100) then
329      num=num+1
330      ndx(num)=ja
331      rd(num)=dist
332      endif
333      endif
334      endif
335    2 continue
336      if(num.gt.0) then
337      do 3 j=1,num-1
338      do 4 k=j+1,num
339      if(ndx(j).gt.ndx(k)) then
340      ntmp=ndx(j)
341      ndx(j)=ndx(k)
342      ndx(k)=ntmp
343      endif
344    4 continue
345    3 continue
346      do 5 j=1,num
347      if(ndx(j).eq.0) goto 6
348      if(ndx(j).ne.iwrk(igrp,j,i)) goto 6
349    5 continue
350      if(iwrk(igrp,num+1,i).eq.0) goto 1
351    6 continue
352      iwrk(igrp,num+1,i)=0
353      do 7 j=1,num
354      iwrk(igrp,j,i)=ndx(j)
355    7 continue
356      if(num.lt.11) then
357      write(lfnloc,1001) time,ia,(ndx(j),j=1,num)
358 1001 format(f12.6,11i6)
359      else
360      write(lfnloc,1002) time,ia,(ndx(j),j=1,num)
361 1002 format(f12.6,11i6,/,(18x,10i6))
362      endif
363c      write(lfnloc,1001)
364c     + imol(ia),sgmnam(ia)(11:16),sgmnam(ia)(1:5),sgmnam(ia)(6:10),
365c     + (imol(ndx(j)),sgmnam(ndx(j))(11:16),sgmnam(ndx(j))(1:5),
366c     + sgmnam(ndx(j))(6:10),j=1,num)
367c 1001 format(2(i5,a6,' ',a5,':',a5,' '),/,(t25,i5,a6,' ',a5,':',a5,' '))
368      endif
369    1 continue
370c
371      return
372      end
373      subroutine ana_histo(x,w,ih)
374c
375      implicit none
376c
377#include "ana_common.fh"
378#include "global.fh"
379#include "mafdecls.fh"
380#include "msgids.fh"
381c
382      integer ih
383      real*8 x(nsa,3),w(mwm,mwa,3)
384c
385      integer i,ia,j,ndx
386c
387      dhis=1.2d0*box(3)/real(idhis(ih,2))
388c
389      ito=ldef(idhis(ih,1))
390      if(ito.gt.0) then
391      do 1 i=1,ito
392      ia=idef(idhis(ih,1),i)
393      ndx=(x(ia,3)-rhis)/dhis
394      ihis(ndx,ih)=ihis(ndx,ih)+1
395    1 continue
396      else
397      do 2 j=1,nwm
398      do 3 i=1,-ito
399      ia=idef(idhis(ih,1),i)
400      ndx=(w(j,ia,3)-rhis)/dhis
401      ihis(ndx,ih)=ihis(ndx,ih)+1
402    3 continue
403    2 continue
404      endif
405c
406      return
407      end
408      subroutine ana_order(x)
409c
410      implicit none
411c
412#include "ana_common.fh"
413#include "ana_params.fh"
414#include "global.fh"
415#include "mafdecls.fh"
416#include "msgids.fh"
417c
418      real*8 x(nsa,3)
419c
420      integer i,ia,n,igrp,j
421      real*8 xsijx,xsijy,xsijz,xskjx,xskjy,xskjz
422      real*8 rsij2,rskj2,rsij2i,rskj2i,rsikji,cphi
423c
424      real*8 xa(3),xb(3),xc(3)
425c
426      do 1 i=1,nord
427      igrp=idord(i,1)
428      n=ldef(igrp)
429      do 2 j=1,n
430      ia=idef(igrp,j)
431      xa(1)=xa(1)+x(ia,1)
432      xa(2)=xa(2)+x(ia,2)
433      xa(3)=xa(3)+x(ia,3)
434    2 continue
435      xa(1)=xa(1)/dble(n)-x(idord(i,3),1)
436      xa(2)=xa(2)/dble(n)-x(idord(i,3),2)
437      xa(3)=xa(3)/dble(n)-x(idord(i,3),3)
438      igrp=idord(i,2)
439      n=ldef(igrp)
440      do 3 j=1,n
441      ia=idef(igrp,j)
442      xb(1)=xb(1)+x(ia,1)
443      xb(2)=xb(2)+x(ia,2)
444      xb(3)=xb(3)+x(ia,3)
445    3 continue
446      xb(1)=xb(1)/dble(n)-x(idord(i,3),1)
447      xb(2)=xb(2)/dble(n)-x(idord(i,3),2)
448      xb(3)=xb(3)/dble(n)-x(idord(i,3),3)
449      xc(1)=x(idord(i,4),1)-x(idord(i,3),1)
450      xc(2)=x(idord(i,4),2)-x(idord(i,3),2)
451      xc(3)=x(idord(i,4),3)-x(idord(i,3),3)
452c
453      xsijx=xa(1)-xb(1)
454      xskjx=xc(1)-xb(1)
455      xsijy=xa(2)-xb(2)
456      xskjy=xc(2)-xb(2)
457      xsijz=xa(3)-xb(3)
458      xskjz=xc(3)-xb(3)
459c
460      rsij2=xsijx*xsijx+xsijy*xsijy+xsijz*xsijz
461      rskj2=xskjx*xskjx+xskjy*xskjy+xskjz*xskjz
462      cphi=xsijx*xskjx+xsijy*xskjy+xsijz*xskjz
463      rsij2i=one/rsij2
464      rskj2i=one/rskj2
465      rsikji=one/sqrt(rsij2*rskj2)
466      cphi=cphi*rsikji
467      if(cphi.lt.-one) cphi=-one
468      if(cphi.gt. one) cphi= one
469c
470      rord(i,1)=rord(i,1)+acos(cphi)
471      rord(i,2)=rord(i,2)+cphi*cphi
472c
473    1 continue
474c
475      return
476      end
477      subroutine ana_grpcogdis(sgmnam,imol,isel,wt,x,igr)
478c
479      implicit none
480c
481#include "ana_common.fh"
482#include "ana_params.fh"
483#include "global.fh"
484#include "mafdecls.fh"
485#include "msgids.fh"
486c
487      character*16 sgmnam(nsa)
488      real*8 wt(nsa)
489      integer isel(nsa),igr,imol(msa)
490      real*8 x(nsa,3)
491c
492      integer igrp,jgrp
493      integer i,ia,ja
494      real*8 dist,dx,dy,dz
495      real*8 cogi(3),cogj(3),factor
496      real*8 boxh(3)
497c
498      boxh(1)=half*box(1)
499      boxh(2)=half*box(2)
500      boxh(3)=half*box(3)
501c
502      igrp=igroups(igr,1)
503      jgrp=igroups(igr,2)
504c
505      if(ldef(igrp).lt.0) return
506      if(ldef(jgrp).lt.0) return
507c
508      do 1 i=1,3
509      cogi(i)=zero
510      cogj(i)=zero
511    1 continue
512c
513      do 2 i=1,ldef(igrp)
514      ia=idef(igrp,i)
515      cogi(1)=cogi(1)+x(ia,1)
516      cogi(2)=cogi(2)+x(ia,2)
517      cogi(3)=cogi(3)+x(ia,3)
518    2 continue
519      factor=one/dble(ldef(igrp))
520      cogi(1)=cogi(1)*factor
521      cogi(2)=cogi(2)*factor
522      cogi(3)=cogi(3)*factor
523c
524      do 3 i=1,ldef(jgrp)
525      ja=idef(jgrp,i)
526      cogj(1)=cogj(1)+x(ja,1)
527      cogj(2)=cogj(2)+x(ja,2)
528      cogj(3)=cogj(3)+x(ja,3)
529    3 continue
530      factor=one/dble(ldef(jgrp))
531      cogj(1)=cogj(1)*factor
532      cogj(2)=cogj(2)*factor
533      cogj(3)=cogj(3)*factor
534c
535      dx=abs(cogi(1)-cogj(1))
536      dy=abs(cogi(2)-cogj(2))
537      dz=abs(cogi(3)-cogj(3))
538c
539      if(igroups(igr,6).eq.1) then
540      if(dz.gt.boxh(3)) dz=dz-box(3)
541      elseif(igroups(igr,6).eq.2) then
542      if(dx.gt.boxh(1)) dx=dx-box(1)
543      if(dy.gt.boxh(2)) dy=dy-box(2)
544      elseif(igroups(igr,6).eq.3) then
545      if(dx.gt.boxh(1)) dx=dx-box(1)
546      if(dy.gt.boxh(2)) dy=dy-box(2)
547      if(dz.gt.boxh(3)) dz=dz-box(3)
548      endif
549      dist=sqrt(dx*dx+dy*dy+dz*dz)
550c
551      write(lfngrp,1001) igr,igroups(igr,5),dist,dx,dy,dz
552 1001 format(2i5,4f12.6)
553c
554      return
555      end
556      subroutine ana_grpcogang(sgmnam,imol,isel,wt,x,igr)
557c
558      implicit none
559c
560#include "ana_common.fh"
561#include "ana_params.fh"
562#include "global.fh"
563#include "mafdecls.fh"
564#include "msgids.fh"
565c
566      character*16 sgmnam(nsa)
567      real*8 wt(nsa)
568      integer isel(nsa),igr,imol(msa)
569      real*8 x(nsa,3)
570c
571      integer igrp,jgrp,kgrp
572      integer i,ia,ja,ka
573      real*8 dx
574      real*8 cogi(3),cogj(3),cogk(3),factor
575      real*8 boxh(3)
576      real*8 xsijx,xskjx,xsijy,xskjy,xsijz,xskjz,rsij2,rskj2
577      real*8 cphi,rsij2i,rskj2i,rsikji,phi
578c
579      boxh(1)=half*box(1)
580      boxh(2)=half*box(2)
581      boxh(3)=half*box(3)
582c
583      igrp=igroups(igr,1)
584      jgrp=igroups(igr,2)
585      kgrp=igroups(igr,3)
586c
587      if(ldef(igrp).lt.0) return
588      if(ldef(jgrp).lt.0) return
589      if(ldef(kgrp).lt.0) return
590c
591      do 1 i=1,3
592      cogi(i)=zero
593      cogj(i)=zero
594      cogk(i)=zero
595    1 continue
596c
597      do 2 i=1,ldef(igrp)
598      ia=idef(igrp,i)
599      cogi(1)=cogi(1)+x(ia,1)
600      cogi(2)=cogi(2)+x(ia,2)
601      cogi(3)=cogi(3)+x(ia,3)
602    2 continue
603      factor=one/dble(ldef(igrp))
604      cogi(1)=cogi(1)*factor
605      cogi(2)=cogi(2)*factor
606      cogi(3)=cogi(3)*factor
607c
608      do 3 i=1,ldef(jgrp)
609      ja=idef(jgrp,i)
610      cogj(1)=cogj(1)+x(ja,1)
611      cogj(2)=cogj(2)+x(ja,2)
612      cogj(3)=cogj(3)+x(ja,3)
613    3 continue
614      factor=one/dble(ldef(jgrp))
615      cogj(1)=cogj(1)*factor
616      cogj(2)=cogj(2)*factor
617      cogj(3)=cogj(3)*factor
618c
619      do 4 i=1,ldef(kgrp)
620      ka=idef(kgrp,i)
621      cogk(1)=cogk(1)+x(ka,1)
622      cogk(2)=cogk(2)+x(ka,2)
623      cogk(3)=cogk(3)+x(ka,3)
624    4 continue
625      factor=one/dble(ldef(kgrp))
626      cogk(1)=cogk(1)*factor
627      cogk(2)=cogk(2)*factor
628      cogk(3)=cogk(3)*factor
629c
630      if(igroups(igr,6).gt.0) then
631      do 6 i=1,igroups(igr,6)
632      dx=cogi(i)-cogj(i)
633      if(dx.lt.-boxh(i)) cogi(i)=cogi(i)+box(i)
634      if(dx.gt.boxh(i)) cogi(i)=cogi(i)-box(i)
635      dx=cogk(i)-cogk(i)
636      if(dx.lt.-boxh(i)) cogk(i)=cogk(i)+box(i)
637      if(dx.gt.boxh(i)) cogk(i)=cogk(i)-box(i)
638    6 continue
639      endif
640c
641      xsijx=cogi(1)-cogj(1)
642      xskjx=cogk(1)-cogj(1)
643      xsijy=cogi(2)-cogj(2)
644      xskjy=cogk(2)-cogj(2)
645      xsijz=cogi(3)-cogj(3)
646      xskjz=cogk(3)-cogj(3)
647c
648      rsij2=xsijx*xsijx+xsijy*xsijy+xsijz*xsijz
649      rskj2=xskjx*xskjx+xskjy*xskjy+xskjz*xskjz
650      cphi=xsijx*xskjx+xsijy*xskjy+xsijz*xskjz
651      rsij2i=one/rsij2
652      rskj2i=one/rskj2
653      rsikji=one/sqrt(rsij2*rskj2)
654      cphi=cphi*rsikji
655      if(cphi.lt.-one) cphi=-one
656      if(cphi.gt. one) cphi= one
657      phi=acos(cphi)
658c
659      write(lfngrp,1001) igr,igroups(igr,5),phi
660 1001 format(2i5,4f12.6)
661c
662      return
663      end
664
665