real*8 M1(3,3),M2(3,3),M3(3,3),M(3,3),ACC(3,3) REAL*8 P1(3),P2(3),T(3),x,y,z,big,theta,pi WRITE(6,*)'ENTER COORDINATES OF 2 POINTS ON THE NCS AXIS:' READ(5,*)(P1(I),I=1,3) WRITE (6,*)'P1=',(P1(I),I=1,3) READ(5,*)(P2(I),I=1,3) WRITE (6,*)'P2=',(P2(I),I=1,3) C INITIALIZE M1-M3 AS IDENTITY MATRICES: DO 80 I=1,3 DO 70 J=1,I-1 M1(I,J)=0. M1(J,I)=0. M2(I,J)=0. M2(J,I)=0. M3(I,J)=0. 70 M3(J,I)=0. M1(I,I)=1.0 M2(I,I)=1.0 80 M3(I,I)=1.0 C TRANSLATE BY -P1: DO 100 I=1,3 100 P2(I)=P2(I)-P1(I) C BETTER WAY: ROTATE ONTO WHICH AXIS HAS BIGGER COMPONENT ORIGINALLY- C L,M,N REPRESENT THE NUMBERS 1,2,3 IN ORDER OF MAG OF COMPONENT OF P1-P2 BIG=0 DO 333 I=1,3 IF (ABS(P2(I)).GE.BIG) THEN BIG=ABS(P2(I)) N1=I ENDIF 333 CONTINUE BIG=0 DO 334 I=1,3 IF (I.eq.N1) GOTO 334 IF (ABS(P2(I)).GE.BIG) THEN BIG=ABS(P2(I)) N2=I ENDIF 334 CONTINUE BIG=0 DO 335 I=1,3 IF (I.eq.N1) GOTO 335 IF (I.eq.N2) GOTO 335 C IF (ABS(P2(I)).GE.BIG) THEN C BIG=ABS(P2(I)) N3=I C ENDIF 335 CONTINUE WRITE (6,*)'P2=',(P2(I),I=1,3) WRITE(6,*)'ORDER OF AXES:',N1,N2,N3 C ROTATE P2 INTO N1-N2 PLANE- SAVE MATRIX IN M1: IF (P2(N3).NE.0.) THEN THETA=DATAN(P2(N3)/P2(N1)) M1(N1,N1)=DCOS(THETA) M1(N3,N3)=M1(N1,N1) M1(N1,N3)=DSIN(THETA) M1(N3,N1)=-M1(N1,N3) X=M1(N1,N1)*P2(N1)+M1(N1,N3)*P2(N3) Z=M1(N3,N1)*P2(N1)+M1(N3,N3)*P2(N3) P2(N1)=X P2(N3)=Z WRITE (6,*)'P2=',(P2(I),I=1,3) write(6,*)'M1:' DO 340 I=1,3 340 WRITE(6,9) (M1(I,J),J=1,3) IF (DABS(Z).GT.1E-12) STOP 'M1 MATRIX FAILED TO ZERO P2(N3)!' ENDIF C ROTATE P2 ONTO N1- SAVE MATRIX IN M2: IF (P2(N2).NE.0.) THEN THETA=ATAN(P2(N2)/P2(N1)) write(6,*)'p2(n2), p2(n1),theta:',p2(n2), p2(n1),theta M2(N1,N1)=COS(THETA) M2(N2,N2)=M2(N1,N1) M2(N1,N2)=SIN(THETA) M2(N2,N1)=-M2(N1,N2) X=M2(N1,N1)*P2(N1)+M2(N1,N2)*P2(N2) Y=M2(N2,N1)*P2(N1)+M2(N2,N2)*P2(N2) P2(N1)=X P2(N2)=Y WRITE (6,*)'P2=',(P2(I),I=1,3) write(6,*)'M2:' DO 350 I=1,3 350 WRITE(6,9) (M2(I,J),J=1,3) IF (DABS(Y).GT.1E-12) STOP 'M2 MATRIX FAILED TO ZERO P2(N2)!' ENDIF C SET UP NCS ROTATION MATRIX- IN THIS CASE 180 DEGREES - ABOUT N1 AXIS write(6,*)'Enter rotation angle:' read(5,*)theta pi=4*datan(1.0D0) write(6,*)'pi=',pi theta=pi*theta/180. M3(N2,N2)=cos(theta) M3(N3,N3)=cos(theta) M3(N2,N3)=-sin(theta) M3(N3,N2)=sin(theta) C FORM PRODUCT OF M2*M1 IN ACC DO 300 I=1,3 DO 300 K=1,3 ACC(I,K)=0 DO 300 J=1,3 300 ACC(I,K)=ACC(I,K)+M2(I,J)*M1(J,K) C COPY PRODUCT INTO M: DO 400 I=1,3 DO 400 K=1,3 400 M(I,K)=ACC(I,K) write(6,*)'M = M2*M1:' DO 360 I=1,3 360 WRITE(6,9) (M(I,J),J=1,3) C FORM PRODUCT OF M3*M IN ACC DO 500 I=1,3 DO 500 K=1,3 ACC(I,K)=0 DO 500 J=1,3 500 ACC(I,K)=ACC(I,K)+M3(I,J)*M(J,K) C COPY PRODUCT INTO M: DO 600 I=1,3 DO 600 K=1,3 600 M(I,K)=ACC(I,K) write(6,*)'M = M3*M2*M1:' DO 361 I=1,3 361 WRITE(6,9) (M(I,J),J=1,3) C FORM PRODUCT OF M2'*M IN ACC DO 700 I=1,3 DO 700 K=1,3 ACC(I,K)=0 DO 700 J=1,3 700 ACC(I,K)=ACC(I,K)+M2(J,I)*M(J,K) C COPY PRODUCT INTO M: DO 800 I=1,3 DO 800 K=1,3 800 M(I,K)=ACC(I,K) write(6,*)'M = M2t*M3*M2*M1:' DO 362 I=1,3 362 WRITE(6,9) (M(I,J),J=1,3) C FORM PRODUCT OF M1'*M IN ACC DO 900 I=1,3 DO 900 K=1,3 ACC(I,K)=0 DO 900 J=1,3 900 ACC(I,K)=ACC(I,K)+M1(J,I)*M(J,K) C COPY PRODUCT INTO M: DO 950 I=1,3 DO 950 K=1,3 950 M(I,K)=ACC(I,K) write(6,*)'M = M1t*M2t*M3*M2*M1:' DO 363 I=1,3 363 WRITE(6,9) (M(I,J),J=1,3) C OPERATE ON -P1 WITH M, PUT IN P2: P2=[I-M]*P1 DO 1000 I=1,3 P2(I)=P1(I) DO 1000 J=1,3 1000 P2(I)=P2(I)-M(I,J)*P1(J) WRITE(6,*)' NCS OPERATOR:' WRITE(6,*)' ROTATE MATRIX' DO 1100 I=1,3 1100 WRITE(6,9) (M(I,J),J=1,3) WRITE(6,*)' TRANSLATE ' WRITE(6,9) (P2(I),I=1,3) 9 FORMAT(' ',3F14.7) write(6,*)' this operator will be saved as an O datablock rt.o' WRITE(6,*)' Note this operator is for orthogonal coord, A translation' WRITE(6,*)' You may have to convert to fractional for your cell. ' open (unit=3,name='rt.o',status='unknown') WRITE(3,*) '.SPACE_GROUP_operators r 12 (3f14.7)' do 1150 j=1,3 1150 write(3,1151)(M(I,J),I=1,3) write(3,1151)(P2(I),I=1,3) 1151 FORMAT(3f14.7) END