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