PROGRAM VOLCAL C VOLCAL.FOR C C PROGRAM TO CALCULATE POLYHEDRAL VOLUMES AND DISTORTION PARAMETERS C WITH ERRORS (IF DESIRED). C C ORIGINAL PROGRAM BY L. W. FINGER, 9/21/71 C MODIFIED 10/24/73 BY Y. OHASHI C MODIFIED 6/20/79 BY L. W. FINGER - MADE INTERACTIVE AND ERRORS ADDED C C Modified 10/18/96 by L.W.F. to fix problem with more than 3 points on face. C C ***************************************************************** C INPUT DATA FORMAT IF READ FROM FILE C ***************************************************************** C C 1. TITLE CARD (A80) C 2. CELL CARD (6F8.0,2I4) C COLS.1-8,9-16,17-24 = A,B AND C C COLS.25-32,33-40,41-48 = ALPHA,BETA AND GAMMA IN DEG. C COLS.49-52 MP = NUMBER OF POLYHEDRA TO BE PROCESSED C COLS.53-56 IER = NON-ZERO IF ERRORS TO BE CALCULATED C 3. CELL VARIANCE-COVARIANCE MATRIX (6F10.0) (NEEDED ONLY IF IER NON-ZERO) C THIS INPUT IS SIX LINES WITH ONE ROW OF MATRIX/LINE. C 4. CENTRAL ATOM CARD (A4,A2,3F10.5,I4) C COLS.1-6 SITE NAME C COLS.7-16,17-26,27-36 = X,Y AND Z C COLS.37-40 NP = NUMBER OF COORDINATING ATOMS C 5. CENTRAL ATOM SIGMA CARDS (3F9.5) (NEEDED ONLY IF IER NON-ZERO) C COLS.1-9,10-18,19-27 = SIGMA X, Y, AND Z C 6. COORDINATING ATOM CARDS (A4,A2,3F10.5) C COLS. 1-6 SITE NAME C COLS.7-16,17-26,27-36 = X,Y AND Z C 7. COORDINATING ATOM SIGMA CARDS (3F9.5) NEEDED ONLY IF IER NON-ZERO) C COLS.1-9,10-18,19-27 = SIGMA X, Y, AND Z C C CARDS 6 AND 7 SHOULD BE REPEATED FOR EACH COORDINATING ATOM C C IF DESIRED, ANOTHER SET OF DATA CAN FOLLOW C C ****************************************************************** C DIMENSION A(6),VARA(6,6),SITE(2),XC(3),SXC(3) DIMENSION XP(3,20),SXP(3,20),DERDA(6,3),AX(3),SIGMA(3) 1,X(3),PRM(3) CHARACTER TITLE*80,IYES*1,INO*1,ILYES*1,ILNO*1,IPAGE*1,IRESP*1 DATA IYES,INO,ILYES,ILNO,IPAGE,IN/'Y','N','y','n',' ',4/ C 90 WRITE(*,1) 1 FORMAT('0POLYHEDRAL VOLUME AND DISTORTION PARAMETER PROGRAM' 1/' ARE DATA TO BE READ FROM FILE(Y,N)?',$) READ(*,'(a)')IRESP IFILE=1 IF(IRESP.EQ.IYES.OR.IRESP.EQ.ILYES)GO TO 1000 IF(IRESP.NE.INO.AND.IRESP.NE.ILNO)GO TO 90 IFILE=0 WRITE(*,2) 2 FORMAT(' TYPE TITLE?',$) READ(*,'(A)')TITLE WRITE(*,6) 6 FORMAT(' ARE ERRORS TO BE CALCULATED(Y,N)?',$) 100 READ(*,'(A)')IRESP 8 FORMAT(A1) IER=1 IF(IRESP.EQ.IYES.OR.IRESP.EQ.ILYES)GO TO 110 IF(IRESP.NE.INO.AND.IRESP.NE.ILNO)GO TO 100 IER=0 110 WRITE(*, 10) 10 FORMAT(' TYPE CELL CONSTANTS?',$) READ(*, 12)A 12 FORMAT(6F10.0) IF(IER.EQ.0)GO TO 120 WRITE(*, 14) 14 FORMAT(' TYPE VARIANCE MATRIX FOR CELL DATA (1 ROW PER LINE)') READ(*, 12)VARA 120 WRITE(*, 16) 16 FORMAT(' TYPE NO. OF COORDINATING IONS?',$) READ(*, 18)NP 18 FORMAT(I6) IF(NP.EQ.0)STOP WRITE(6,5)IPAGE,TITLE 5 FORMAT(A1,a) IPAGE='1' WRITE(6,13)A 13 FORMAT('0CELL CONSTANTS',6F10.5) IF(IER.NE.0)WRITE(6,15)VARA 15 FORMAT('0CELL COVARIANCE MATRIX'/(6E13.6)) WRITE(*, 20) 20 FORMAT(' TYPE CATION COORDS. AND SITE NAME?',$) READ(*, 22)XC,SITE 22 FORMAT(3F9.0,2A4) WRITE(6,23)SITE,XC 23 FORMAT('0CATION: ',2A4,3F10.5) IF(IER.EQ.0)GO TO 130 WRITE(*, 24) 24 FORMAT(' TYPE SIGMAS FOR SITE?',$) READ(*, 22)SXC WRITE(6,25)SXC 25 FORMAT(5X,'SIGMAS',6X,3F10.5) 130 WRITE(6,27) 27 FORMAT('0ANIONS') DO 140 I=1,NP WRITE(*, 26)I 26 FORMAT(' TYPE ANION COORDS. AND LABEL FOR NO.',I3,'?',$) READ(*, 22)(XP(J,I),J=1,3),SITE WRITE(6,29)I,SITE,(XP(J,I),J=1,3) 29 FORMAT(3X,I5,1X,2A4,3F10.5) IF(IER.EQ.0)GO TO 140 WRITE(*, 28) 28 FORMAT(' SIGMAS?',$) READ(*, 22)(SXP(J,I),J=1,3) WRITE(6,25)(SXP(J,I),J=1,3) 140 CONTINUE GO TO 1410 C C READ DATA FROM FILE C 1000 WRITE(*, 40) 40 FORMAT(' TYPE INPUT FILE NAME?',$) READ(*, 42)TITLE 42 FORMAT(A) OPEN(UNIT=IN,FILE=TITLE,STATUS='OLD') 1010 READ(IN,'(A)',END=260)TITLE READ(IN,44)A,MP,IER NCOUNT=0 44 FORMAT(6F8.0,2I4) WRITE(6,5)IPAGE,TITLE IPAGE='1' WRITE(6,13)A IF(IER.EQ.0)GO TO 1020 READ(IN,12)VARA WRITE(6,15)VARA 1020 READ(IN,46)SITE,XC,NP 46 FORMAT(A4,A2,3F10.0,I4) WRITE(6,23)SITE,XC IF(IER.EQ.0)GO TO 1030 READ(IN,22)SXC WRITE(6,25)SXC 1030 DO 1040 I=1,NP READ(IN,46)SITE,(XP(J,I),J=1,3) WRITE(6,29)I,SITE,(XP(J,I),J=1,3) IF(IER.EQ.0)GO TO 1040 READ(IN,22)(SXP(J,I),J=1,3) WRITE(6,25)(SXP(J,I),J=1,3) 1040 CONTINUE C 1410 CALL PRMCAL(0,A,A(4),XC,NP,XP,PRM(1),PRM(2)) IF(IER.EQ.0)GO TO 250 C FUNCTIONS CALULATED - NOW GET ERRORS DO 142 I=1,3 142 SIGMA(I)=0.0 DO 155 I=1,6 DX=SQRT(VARA(I,I))*.1 IF(DX.EQ.0.0)GO TO 155 A(I)=A(I)+DX CALL PRMCAL(1,A,A(4),XC,NP,XP,X(1),X(2)) A(I)=A(I)-DX DO 150 J=1,3 150 DERDA(I,J)=(X(J)-PRM(J))/DX 155 CONTINUE DO 170 I=1,6 DO 170 J=1,3 AX(I)=0.0 DO 160 K=1,6 160 AX(I)=AX(I)+DERDA(I,J)*VARA(K,I) 170 SIGMA(J)=SIGMA(J)+AX(I)*DERDA(I,J) DO 225 I=1,3 DX=XC(I)+.1*SXC(I) IF(DX.EQ.0.0)GO TO 225 XC(I)=XC(I)+DX CALL PRMCAL(1,A,A(4),XC,NP,XP,X(1),X(2)) XC(I)=XC(I)-DX DO 220 J=1,3 220 SIGMA(J)=SIGMA(J)+((X(J)-PRM(J))/DX*SXC(I))**2 225 CONTINUE DO 231 I=1,NP DO 231 J=1,3 DX=XP(J,I)+0.1*SXP(J,I) IF(DX.EQ.0.0)GO TO 231 XP(J,I)=XP(J,I)+DX CALL PRMCAL(1,A,A(4),XC,NP,XP,X(1),X(2)) XP(J,I)=XP(J,I)-DX DO 230 K=1,3 230 SIGMA(K)=SIGMA(K)+((X(K)-PRM(K))/DX*SXP(J,I))**2 231 CONTINUE DO 232 I=1,3 232 SIGMA(I)=SQRT(SIGMA(I)) IF(NP.NE.4.AND.NP.NE.6)GO TO 240 WRITE(6,30)SIGMA 30 FORMAT(' SIGMA(S)',3F12.5) GO TO 250 240 WRITE(6,30)SIGMA(1) 250 IF(IFILE.EQ.0)GO TO 120 NCOUNT=NCOUNT+1 IF(NCOUNT-MP)1020,1010,1010 260 STOP END SUBROUTINE PRMCAL(IORG,A,ANG1,XC,NP,YP,VOL,QE) DIMENSION XC(3),A(3),ANG1(3),ANG(3),XP(3,20),VXYZ(6),D(3) 1,DM(20),IH(3) DIMENSION YP(3,20) C DATA IOUT/6/ DATA PI/3.141592/ C 6 FORMAT(1H0,'POLYHEDRAL VOLUME',F13.5/5X,'QUAD. ELONG.',F14.5/5X, 1'ANGLE VARIANCE',F12.5) 8 FORMAT (1H0,'POLYHEDRAL VOLUME',F12.5) 15 FORMAT(1H0,10X,'ATOMS',13X,'DISTANCES',20X,'ANGLE'/12X,'I J', 1 7X,'0-I',9X,'0-J',9X,'I-J',7X,'I-0-J'/) 20 FORMAT(1H ,8X,2I4,3F12.5,F12.3) 25 FORMAT(1H0,12X,'FACE DEFINED',10X,'PLANE NORMAL',5X,'AREA'/13X, 1 'BY POINTS',7X,'OX',6X,'OY',6X,'OZ'/) 26 FORMAT (1H ,10X,3I3,4X,3F8.5,4X,F10.5) C DO 100 I=1,3 100 ANG(I)=COS(ANG1(I)*PI/180.0) SINA=SQRT(1.0-ANG(1)**2) ASTAR=SINA/(A(1)*SQRT(1.0-ANG(1)**2-ANG(2)**2-ANG(3)**2+2.0*ANG(1) 1 *ANG(2)*ANG(3))) C CONVERT POINTS TO ORTHOGONAL COORDS. DO 130 I=1,NP C REFER TO ORIGIN AT CENTRAL ATOM DO 120 J=1,3 120 XP(J,I)=YP(J,I)-XC(J) XP(3,I)=A(3)*XP(3,I)+XP(2,I)*A(2)*ANG(1)+XP(1,I)*A(1)*ANG(2) XP(2,I)=XP(2,I)*A(2)*SINA+XP(1,I)*A(1)*(ANG(3)-ANG(2)*ANG(1))/SINA XP(1,I)=XP(1,I)/ASTAR DM(I)=SQRT(XP(1,I)**2+XP(2,I)**2+XP(3,I)**2) 130 CONTINUE C DISTANCES AND ANGLES IF(IORG.EQ.0)WRITE (IOUT,15) DO 145 I=1,NP-1 DO 145 J=I+1,NP TEMPA=0. TEMPB=0. DO 135 K=1,3 TEMPA=TEMPA+(XP(K,I)-XP(K,J))**2 135 TEMPB=TEMPB+XP(K,I)*XP(K,J) DIS=SQRT(TEMPA) ANGL=TEMPB/(DM(I)*DM(J)) SUMTH=1.0-ANGL*ANGL IF(SUMTH.LT.0.0)SUMTH=0.0 ANGL=57.2958*ATAN2(SQRT(SUMTH),ANGL) IF(IORG.EQ.0)WRITE (IOUT,20) I,J,DM(I),DM(J),DIS,ANGL 145 CONTINUE IF(IORG.EQ.0)WRITE(IOUT,25) C FIND FACES - CHOOSE ANY THREE POINTS VOL=0.0 NP1=NP-1 NP2=NP-2 SUMTH =0.0 SUMTH2=0.0 NA=0 DO 180 I=1,NP2 IH(1)=I I1=I+1 DO 180 J=I1,NP1 J1=J+1 IH(2)=J VXYZ(1)=XP(1,J)-XP(1,I) VXYZ(2)=XP(2,J)-XP(2,I) VXYZ(3)=XP(3,J)-XP(3,I) DO 180 K=J1,NP IH(3)=K NS=0 VXYZ(4)=XP(1,K)-XP(1,I) VXYZ(5)=XP(2,K)-XP(2,I) VXYZ(6)=XP(3,K)-XP(3,I) D(1)=VXYZ(2)*VXYZ(6)-VXYZ(5)*VXYZ(3) D(2)=VXYZ(4)*VXYZ(3)-VXYZ(1)*VXYZ(6) D(3)=VXYZ(1)*VXYZ(5)-VXYZ(4)*VXYZ(2) AREA=0.5*SQRT(D(1)**2+D(2)**2+D(3)**2) Z0=0.5*(XP(1,I)*D(1)+XP(2,I)*D(2)+XP(3,I)*D(3))/AREA C CHECK FOR AND AVOID PLANE THROUGH ORIGIN IF(ABS(Z0).LT.1.0E-5) GO TO 180 Factor = 3.0 DO 170 L=1,NP IF(L.EQ.I.OR.L.EQ.J.OR.L.EQ.K) GO TO 170 C CALCULATE DISTANCE OF POINT L FROM PLANE OF IJK Z=0.5*((XP(1,I)-XP(1,L))*D(1)+(XP(2,I)-XP(2,L))*D(2)+(XP(3,I)-XP(3 1,L))*D(3))/AREA c c Code added 18-Oct-1996 C C Z and Z0 must have the same sign c if (z * z0 .lt. -0.001)go to 180 if (abs(z * z0) .lt. 0.001)then c If more than 3 corners on this face, the area will be counted twice. c Change factor to handle this case. factor = 6.0 endif c 170 CONTINUE c C ALL POINTS ON SAME SIDE, THUS IJK ARE FACE C C DIRECTION COSINES OF PLANE NORMAL D(1)=D(1)/(AREA*2.) D(2)=D(2)/(AREA*2.) D(3)=D(3)/(AREA*2.) IF(IORG.EQ.0)WRITE (IOUT,26) I,J,K,D,AREA VOL=VOL+AREA*ABS(Z0)/FACTOR DO 172 L=1,2 L1=L+1 NM=IH(L) DO 172 M=L1,3 MN=IH(M) TEMP=0.0 DO 171 N=1,3 171 TEMP=TEMP+XP(N,NM)*XP(N,MN) TEMP=TEMP/(DM(NM)*DM(MN)) TEMP=57.2958*ATAN2(SQRT(1.0-TEMP*TEMP),TEMP) SUMTH=SUMTH+TEMP 172 SUMTH2=SUMTH2+TEMP**2 NA=NA+3 180 CONTINUE IF((NP-4)*(NP-6))230,190,230 C ALL FACES FOUND, CALCULATE QUADRATIC ELONGATION FOR TETR AND OC 190 SUM=0.0 DO 200 I=1,NP 200 SUM=SUM+DM(I)**2 IF(NP.EQ.4) GO TO 210 C SET UP FOR OCTAHEDRON CONS=0.75 CONS2=90.0 GO TO 220 C SET UP FOR TETRAHEDRON 210 CONS=9.0*SQRT(3.0)/8.0 CONS2=109.45 220 VLO=EXP(2.0*ALOG(CONS*VOL)/3.0) QE=SUM/(NP*VLO) NA=(NA+1)/2 SUMTH2=0.5*SUMTH2 SUMTH =0.5*SUMTH SIGA=(SUMTH2-2.0*CONS2*SUMTH+NA*CONS2**2)/(NA-1) IF(IORG.EQ.0)WRITE(IOUT,6) VOL,QE,SIGA GO TO 240 230 IF(IORG.EQ.0)WRITE(IOUT,8) VOL 240 CONTINUE RETURN END