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