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