subroutine gethist(x,y,z,xmin,xmax,ncoord,tcoord,icoord, 1hist,ncoord0,NTN) implicit none integer*4 ncoord0,KK,NTN integer*4 i,ncoord,hist(*),icoord(ncoord0,4),j,L0,L1,L2,L3, 1tcoord(*) real*8 x(*),y(*),z(*),xmin,xmax, 1A(3),P(3),E0(3),PI,SSPP,AD,SSAE0,SSPE0,XX(3),YY(3),SSE0E0, 1SSXXXX,SSYYYY,SSXXYY,APX,APY,APZ,SIGN,SSE0AP,CB,CA,BANGLE,ANGLE, 1xc,x12,y12,z12 include 'sizes.i' include 'bound.i' include 'cell.i' PI=3.141592653589D0 do 1 i=1,ncoord if(tcoord(i).eq.2)then x12=x(icoord(i,1))-x(icoord(i,2)) y12=y(icoord(i,1))-y(icoord(i,2)) z12=z(icoord(i,1))-z(icoord(i,2)) xc=x12**2+y12**2+z12**2 c c find the smallest distance from all cells c if(use_replica)then c do 2 j=1,ncell c call image(x12,y12,z12,j) c if(x12**2+y12**2+z12**2.lt.xc)xc=x12**2+y12**2+z12**2 c endif xc=dsqrt(xc) endif if(tcoord(i).eq.4)then L0=icoord(i,1) L1=icoord(i,2) L2=icoord(i,3) L3=icoord(i,4) A(1)=X(L2)-X(L3) A(2)=Y(L2)-Y(L3) A(3)=Z(L2)-Z(L3) P(1)=X(L0)-X(L1) P(2)=Y(L0)-Y(L1) P(3)=Z(L0)-Z(L1) E0(1)=X(L1)-X(L2) E0(2)=Y(L1)-Y(L2) E0(3)=Z(L1)-Z(L2) SSE0E0=E0(1)*E0(1)+E0(2)*E0(2)+E0(3)*E0(3) DO 8 KK=1,3 8 E0(KK)=E0(KK)/DSQRT(SSE0E0) SSPP=P(1)*P(1)+P(2)*P(2)+P(3)*P(3) AD=DSQRT(SSPP) SSAE0=A(1)*E0(1)+A(2)*E0(2)+A(3)*E0(3) SSPE0=P(1)*E0(1)+P(2)*E0(2)+P(3)*E0(3) DO 9 KK=1,3 XX(KK)=A(KK)-E0(KK)*SSAE0 9 YY(KK)=P(KK)-E0(KK)*SSPE0 SSXXXX=XX(1)*XX(1)+XX(2)*XX(2)+XX(3)*XX(3) SSYYYY=YY(1)*YY(1)+YY(2)*YY(2)+YY(3)*YY(3) SSXXYY=XX(1)*YY(1)+XX(2)*YY(2)+XX(3)*YY(3) APX=XX(2)*YY(3)-XX(3)*YY(2) APY=XX(3)*YY(1)-XX(1)*YY(3) APZ=XX(1)*YY(2)-XX(2)*YY(1) SSE0AP=APX*E0(1)+APY*E0(2)+APZ*E0(3) SIGN=1.0D0 IF (SSE0AP.NE.0.0D0)SIGN=-SSE0AP/ABS(SSE0AP) CB=0.0D0 IF (AD.NE.0.0D0)CB=-SSPE0/AD CA=0.0D0 IF (SSXXXX.NE.0.0D0.AND.SSYYYY.NE.0.0D0) 1 CA=-SSXXYY/DSQRT(SSXXXX)/DSQRT(SSYYYY) BANGLE=0.0D0 ANGLE=0.0D0 IF (CB.EQ.-1.0D0)BANGLE=PI IF (ABS(CB).LT.1.0D0)BANGLE=ACOS(CB) IF (CA.LE.-1.0D0) ANGLE=PI IF (ABS(CA).LT.1.0D0) ANGLE=ACOS(CA) ANGLE=180.0D0* ANGLE/PI BANGLE=180.0D0*BANGLE/PI ANGLE=SIGN*ANGLE if(ANGLE.lt.-180.0d0)ANGLE=360.0d0+ANGLE xc=ANGLE endif j=int((xc-xmin)/(xmax-xmin)*dble(NTN))+1 1 if(j.gt.0.and.j.lt.NTN)hist(j)=hist(j)+1 return end c wwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwww