******************************************************************************* * PROGRAM DRAWCRYS.FOR * * THIS PROGRAM WILL TAKE CRYSTALLOGRAPHIC DATA * * (XYZ COORDINATES AND UNIT CELL INFORMATION) * * FROM AN ORFEE INPUT FILE AND PREPARE A FILE * * OF TRANSFORMED UNIT CELL COORDINATES * * THIS IS ONLY A PRELIMINARY VERSION OF THE PROGRAM * * PLEASE SEND SUGGESTIONS FOR IMPROVEMENT TO ME * * DAVID CLOSE * * PHYSICS DEPT. BOX 70652 * * EAST TENN. ST. UNIV. * * JOHNSON CITY, TN 37614 * * e-mail R29CLOSE@ETSU.EAST-TENN-ST.EDU * * phone 423-929-5646 * * * * THIS WORK WAS SUPPORTED BY PHS GRANT CA 36810-011, * * AWARDED BY NATIONAL CANCER INSTITUTE, DHHS * * * ******************************************************************************* * FORTRAN CODE FOR PREPARING AN INPUT FILE FOR ORFEE2 * * TWO SETS OF DATA ARE TO BE READ IN: * * 1) FIRST THE NORMAL ORFEE INPUT * * 2) SECOND, THE OUTPUT FROM THE SUBROUTINE 23INPUT.FOR * * THE OUTPUT FROM BOTH SECTIONS IS READ BY ORFEE2 * * DAVID CLOSE MARCH 1996 * ******************************************************************************* CHARACTER*64 INFIL1 CHARACTER*72 LINE(20) CHARACTER*4 BLANK CHARACTER*1 ANS DIMENSION TITLE(20) INTEGER OPTION, ORIGIN OPTION = 23 WRITE (*,100) 100 FORMAT (//,2X,'PROGRAM TO PRODUCE AN INPUT FILE FOR BABEL', 1/,2X,'THIS FILE IS THEN RUN IN ALCHEMY') 200 WRITE (*,210) 210 FORMAT (/,5X, 'INPUT FILE NAME- ',/ 1,7X,'THIS IS JUST A NORMAL ORFEE INPUT FILE ' ) READ (*,220) INFIL1 220 FORMAT (A) OPEN (5,FILE=INFIL1,ERR=230) GO TO 240 230 WRITE (*, 235) 235 FORMAT (/,2X,'CAN'T FIND ORFEE FILE', 1/,2X,'DID YOU FORGET DIRECTORY PATH? ', 2/,2X,'TRY AGAIN ? (Y/N) ') READ (*,238) ANS 238 FORMAT (A) IF (ANS .EQ. 'Y' .OR. ANS .EQ. 'y') GO TO 200 IF (ANS .EQ. 'N' .OR. ANS .EQ. 'n') GO TO 2000 240 OPEN (6,FILE='ORFEE2.TMP',STATUS= 'UNKNOWN') READ (5,250) (TITLE(I),I=1,20) 250 FORMAT (20A4) WRITE (6,250) (TITLE(I),I=1,20) READ (5,270) NUM,NSITES,NATOMS 270 FORMAT (I3,6X,2I3) WRITE (6,270) NUM,NSITES,NATOMS C**** NOW READ IN (NATOMS + NSITES + 1) LINES IN ORFEE INPUT*** N = NSITES+NATOMS+1 J = 1 280 READ (5,290) (LINE(I),I=1,20) WRITE (6,290) (LINE(I),I=1,20) 290 FORMAT (20A4) IF (J .EQ. N) GO TO 300 J = J + 1 GO TO 280 300 CONTINUE ******************************************************************************* * FORTRAN CODE FOR INPUTING DATA TO OPTION 23 OF ORFEE * * WHICH ALLOWS ONE TO GENERATE SYMMETRY RELATED * * SITES FROM THE ORIGINAL INPUT DATA * * THE NOTATION IS AS IN ORFEE. 55501 SPECIFIES * * THE ORIGIN OF THE FIRST MOLECULE. OTHER MOLECULES * * IN THE UNIT CELL ARE 5550X. MOLECULES OUT OF THE * * UNIT CELL ARE (FOR EXAMPLE ONE UNIT UP THE * * AXIS) ARE DESIGNATED 55601 * * DAVID CLOSE JULY 1995 (revised Oct. 1995) * ******************************************************************************* NUMMOLE = 0 330 CONTINUE WRITE (*,340) 340 FORMAT (///) WRITE (*,350) 350 FORMAT (5X,'MENU: ',/) WRITE (*,360) 360 FORMAT (5X,'S) INPUT SYMMETRY OPERATION (LIKE 01--0X) ',/) WRITE (*,370) 370 FORMAT (5X,'T) INPUT SPECIFIC TRANSFORMATION (LIKE 555--XXX') WRITE (*,380) 380 FORMAT (/,5X, 'Q) QUIT',/) READ (*,390) CHOICE 390 FORMAT (1A1) IF (CHOICE .EQ. 'S' .OR. CHOICE .EQ. 's') THEN GO TO 500 END IF IF (CHOICE .EQ. 'T' .OR. CHOICE .EQ. 't') THEN GO TO 800 END IF IF (CHOICE .EQ. 'Q' .OR. CHOICE .EQ. 'q') THEN GO TO 1000 END IF IF (CHOICE .NE. 'S' .OR. CHOICE .NE. 's' 1 .OR. CHOICE .NE. 'T' .OR. CHOICE .NE. 't' 2 .OR. CHOICE .NE. 'Q' .OR. CHOICE .NE. 'q') THEN GO TO 330 END IF * THIS SECTION TRANSFORMS ATOMS WITHIN THE UNIT CELL VIA THE SYMMETRY * OPERATIONS IN ORFEE. 500 WRITE (*, 550) 550 FORMAT (5X,'THIS IS OPT. (S), SO MUST SELECT SYMMETRY OPERATION' 1,/,8X,'SPECIFY SITE 1 AS 01, ETC.') READ (*,510) ISITEX 510 FORMAT (I2) IATOMS = NATOMS-4 ISITE1 = 1 WRITE (*,455) IATOMS,ISITE1,ISITEX 455 FORMAT (//,5X, 'WILL TRANSFORM',I3,1X,'ATOMS',/, 15X,'FROM SITE',I2,2X,'TO SITE',1X,I2) ORIGIN = 555 K = 0 DO 460 I = 5,NATOMS 460 WRITE (6,470) OPTION,I,ORIGIN,K,ISITE1,I,ORIGIN,K,ISITEX 470 FORMAT (I5,I5,I3,I1,I1,I5,I3,I1,I1) NUMMOLE = NUMMOLE + 1 GO TO 330 * THIS SECTION TRANSFORMS ATOMS OUTSIDE OF THE UNIT CELL VIA OPERATIONS * FROM THE MOLECULE AT THE ORIGIN (55501). 800 WRITE (*,810) 810 FORMAT (//,5X,'THIS IS OPT. (T), SO MUST SELECT A TRANSFORMATION 1',/,8X, 'OUTSIDE THE UNIT CELL',//) WRITE (*,820) 820 FORMAT (5X, 'SITES ARE DESSIGNATED HERE AS 55501,55601,55602', 1/,8X,'SO MUST SELECT BOTH SITE AND TRANSFORMATION') WRITE (*,830) 830 FORMAT (//,5X,'WHICH SITE SHOULD BE SPECIFIED ') READ (*,840) ISITEX 840 FORMAT (I2) WRITE (*,850) 850 FORMAT (5X, 'WHICH TRANSFORMATION (566?) SHOULD BE SPECIFIED') READ (*,860) ITRANS 860 FORMAT (I3) IATOMS = NATOMS-4 ORIGIN = 555 ISITE1 = 1 K = 0 WRITE (*,870) IATOMS,ISITE1,ISITEX 870 FORMAT (//,5X, 'WILL TRANSFORM',I3,1X'ATOMS',/, 15X,'FROM SITE',I4,2X,'TO SITE',1X,I4) DO 880 I = 5, NATOMS 880 WRITE (6,890) OPTION,I,ORIGIN,K,ISITE1,I,ITRANS,K,ISITEX 890 FORMAT (I5,I5,I3,I1,I1,I5,I3,I1,I1) NUMMOLE = NUMMOLE + 1 GO TO 330 1000 BLANK = '0000' WRITE (6,1010) BLANK 1010 FORMAT (A4) CLOSE (UNIT=6) WRITE (*,1020) NUMMOLE 1020 FORMAT (//,2X,'HAVE CALLED FOR',I2,' SYMMETRY OPERATIONS',//) CALL ORFEE2 CALL MAKEBAB (INFIL1,NUMMOLE) 2000 CONTINUE END ****************************************************************************** * THIS IS AN OLD COPY OF OFREE. I HAVE MODIFIED IT SEVERAL * * TIMES FOR MY WORK IN EPR/ENDOR. THERE IS A LOT CONTAINED * * HERE THAT IS ENTIRELY UNNECESSARY FOR THE PRESENT USE. * * HOWEVER, IT IS NOT TRIVIAL TO GO INTO THIS PROGRAM TO * * DELETE VARIOUS SUBROUTINES, BECAUSE MANY OF THEM ARE * * USED BY OTHER PARTS OF THE PROGRAM. SO BEWARE. * ****************************************************************************** SUBROUTINE ORFEE2 C ORFFE3 CRYSTALLOGRAPHIC FUNCTION AND ERROR PROGRAM C OAK RIDGE NATIONAL LABORATORY, OAK RIDGE, TENNESSEE 37830 C BASED ON ORFFE BY W R BUSING, K O MARTIN, AND H A LEVY C WITH MODIFICATIONS BY G M BROWN, C K JOHNSON, AND W E THIESSEN C JANUARY, 1971 VERSION C DIMENSIONED FOR 200 ATOMS, 110 VARIABLES , 160 PARAMETERS , AND C 9 SCALE FACTORS C COMMON A(6),AA(3,3),AM(21),BB(3,3),DA(6),DFDA(6),DFDP( 160) 1,DP( 160),E,E1,F,FX,P( 160),SAVEA,SAVEP ,VARA,VARP REAL*8 A,AA,AM,BB,DA,DFDA,DFDP,DP,CHEM REAL*8 E,E1,F,FX,P,SAVEA,SAVEP,VARA,VARP COMMON CD(20,2),CHEM(200),DMAX,FS(3,3,48),IAM,IH,IN(261) A,INSAVE,IPM,ITF,JB,JBP,JX,JXP,KD(20,4),KI1( 160),KI2( 160) B,KK,KKD,LNZ,NA,NG,NM,NP,NQ,NS,NSP,NV,PM( 5),ROW(9) C,SCALE,TITLE(18),TS(3,48) D,LOCX(201),ITA(200),IGM(200),INCD,IEXT,ITO,IKE(9),NPX,NI,NJ,NO DIMENSION IDL(200) NI=5 NJ=3 NO=6 OPEN(5,FILE='ORFEE2.TMP') OPEN(6,FILE='ORFEE2.OUT',STATUS='UNKNOWN') C READ AND PUT OUT TITLE AND CONTROL CARD 200 READ (NI,2)(TITLE(I),I=1,18) 2 FORMAT (18A4) WRITE(NO,4)(TITLE(I),I=1,18) 4 FORMAT (1 18A4) READ(NI,6) INCD,IPM,IAM,NS,NA,ITF 6 FORMAT (24I3) DO 205 I=1,261 205 IN(I)=0 C*****NOTE THE LIMIT HERE MUST BE SET EQUAL TO THE MAXIMUM NUMBER OF C PARAMETERS ALLOWED IN THE DIMENSION STATEMENTS ABOVE. DO 210 I=1,160 DP(I)=0.D0 KI1(I)=0 210 P(I)=0.D0 IF(INCD) 220,215,220 215 WRITE(NO,8) 8 FORMAT ('INPUT DATA TO BE READ FROM OR FLS TAPE') GO TO 230 220 WRITE(NO,10) 10 FORMAT('INPUT DATA TO BE READ FROM CARDS') IF(ITF.EQ.0) GO TO 225 WRITE(NO,12) 12 FORMAT('POSITION AND TEMPERATURE FACTOR PARAMETERS WILL BE READ') GO TO 230 225 WRITE(NO,14) 14 FORMAT('POSITION PARAMETERS ONLY WILL BE READ') 230 WRITE(NO,16) 16 FORMAT ('VARIANCE-COVARIANCE MATRIX AND PARAMETER SELECTION INF' 1'ORMATION WILL') IF(IPM)235,240,245 235 WRITE(NO,18) 18 FORMAT(' BE USED WITHOUT COVARIANCE') GO TO 250 240 WRITE(NO,20) 20 FORMAT (' NOT BE USED') GO TO 250 245 WRITE(NO,22) 22 FORMAT (' BE USED') 250 WRITE(NO,24)NS 24 FORMAT ('NUMBER OF SYMMETRY CARDS IS',I3) WRITE(NO,26) 26 FORMAT ('CELL PARAMETER ERRORS ARE') IF(IAM-1)255,260,265 255 WRITE(NO,28) 28 FORMAT (' NOT TO BE USED') GO TO 270 260 WRITE(NO,30) 30 FORMAT (' TO BE READ IN THE FORM OF STANDARD ERRORS') GO TO 270 265 WRITE(NO,32) 32 FORMAT (' TO BE READ IN THE FORM OF A VARIANCE-COVARIANCE 17H MATRIX') 270 IF(INCD)310,310,275 C READ PARAMETERS AND VARIANCE-COVARIANCE MATRIX FROM CARDS C READ ATOM PARAMETERS FROM CARDS 275 LOCX(1)=1 DO 285 I=1,NA LX=LOCX(I) LE=LX+2 READ(NI,34) CHEM(I),(P(L),L=LX,LE) 34 FORMAT(A6,21X,3F9.6) IF(ITF.NE.0) GO TO 280 LOCX(I+1)=LE+1 WRITE(NO,36)CHEM(I),I,LX,LE,(P(L),L=LX,LE) 36 FORMAT(1X,A6,1X,'ATOM',I4,2X,'PARAM',I4,1X,'-',I4,9F9.6) GO TO 285 280 LB=LX+3 LE=LB+5 READ(NI,38)(P(L),L=LB,LE),ITA(I),IGM(I) 38 FORMAT(6F9.6,9X,2I3) IF(ITA(I).EQ.0) ITA(I)=ITF IF(ITA(I).EQ.1) LE=LB LOCX(I+1)=LE+1 WRITE(NO,36)CHEM(I),I,LX,LE,(P(L),L=LX,LE) IF(IGM(I).EQ.0) GO TO 285 LB=LOCX(I+1) LE=LB+9 READ(NI,40)(P(L),L=LB,LE) 40 FORMAT(5F14.10,10X) LOCX(I+1)=LB+10 WRITE(NO,42) LB,LE,(P(L),L=LB,LE) 42 FORMAT(15X,7H PARAMI4,2H -I4,5F14.10/33X,5F14.10) 285 CONTINUE NP=LE IF(IPM)295,335,290 290 WRITE(NO,44) 44 FORMAT('THIS PROGRAM WILL NOT READ VARIANCE-COVARIANCE MATRIX FRO 1M CARDS') STOP C READ STANDARD ERRORS OF ATOM PARAMETERS FROM CARDS 295 DO 300 I=1,NA LB=LOCX(I) LE=LB+2 READ(NI,46)(DP(L),L=LB,LE) 46 FORMAT(27X,3F9.6) IF(ITF.EQ.0) GO TO 300 LB=LB+3 LE=LB+5*(ITA(I)-1) READ(NI,38)(DP(L),L=LB,LE) IF(IGM(I).EQ.0) GO TO 300 LB=LE+1 LE=LB+9 READ(NI,40)(DP(L),L=LB,LE) 300 CONTINUE NV=0 DO 305 I=1,NP IF(DP(I).EQ.0.) GO TO 305 NV=NV+1 KI1(I)=1 305 CONTINUE GO TO 335 C READ PARAMETERS AND VARIANCE-COVARIANCE MATRIX C FROM OR XFLS3 TAPE 310 READ(NJ,6) NQ,NA,NP,NV,IEXT,ITO,NPX IF(IEXT.NE.0)READ(NJ,6)(IKE(I),I=1,NQ) READ(NJ,48) (CHEM(I),I=1,NA) 48 FORMAT(12A6) READ(NJ,6) (ITA(I),IGM(I),IDL(I),I=1,NA) READ(NJ,50) (P(I),I=1,NP) 50 FORMAT (10D12.8) J=NQ+1 IF(IEXT.EQ.0) GO TO 320 DO 315 I=1,NQ J=J+1 IF(IKE(I).NE.0) J=J+5 315 CONTINUE 320 IF(ITO.NE.0) J=J+1 DO 325 I=1,NA LOCX(I)=J+3 J=J+12 IF(ITA(I).EQ.1) J=J-5 IF(IGM(I).NE.0) J=J+10 325 CONTINUE LOCX(NA+1)=J IF(IPM)330,335,330 330 READ(NJ,52)(KI1(I),I=1,NP) 52 FORMAT (72I1) NM=(NV*(NV+1))/2 READ(NJ,54)(PM(K),K=1,NM) 54 FORMAT (8E15.7) C READ AND PUT OUT CELL PARAMETERS 335 READ (NI,56)(A(I),I=1,6) 56 FORMAT (6D9.4) WRITE(NO,4)(TITLE(I),I=1,18) WRITE(NO,58)(A(I),I=1,6) 58 FORMAT(' CELL PARAMETERS'/3F11.4,3F16.8) C STORE METRIC TENSOR CALL SETP(P) CALL SETA(A) CALL STOAA C INVERT METRIC TENSOR CALL STOBB IF(IAM-1)385,340,355 C READ STANDARD ERRORS OF CELL PARAMETERS 340 DO 345 I=1,21 345 AM(I)=0.D0 READ (NI,56)AM(1),AM(7),AM(12),AM(16),AM(19),AM(21) WRITE(NO,60)AM(1),AM(7),AM(12),AM(16),AM(19),AM(21) 60 FORMAT (' STANDARD ERRORS, RESPECTIVELY, OF THE ABOVE CELL 111H PARAMETERS'/1H03F11.4,3F16.8) DO 350 I=1,21 350 AM(I)=AM(I)*AM(I) GO TO 375 C READ VARIANCE-COVARIANCE MATRIX FOR CELL PARAMETERS 355 READ (NI,62)(AM(I),I=1,21) 62 FORMAT (8D9.4) WRITE(NO,64) 64 FORMAT (' VARIANCE-COVARIANCE MATRIX FOR CELL PARAMETERS') IJ=1 DO 370 I=1,6 DO 360 J=1,6 360 ROW(J)=0.0 DO 365 J=1,6 ROW(J)=AM(IJ) 365 IJ=IJ+1 370 WRITE(NO,66)(ROW(J),J=1,6) 66 FORMAT (1X,6F11.4) C COMPUTE CELL PARAMETER INCREMENTS USED TO OBTAIN DERIVATIVES 375 K=1 L=6 DO 380 I=1,6 DA(I)=0.01D0*DSQRT(AM(K)) K=K+L 380 L=L-1 385 IF(NS)395,395,390 C READ AND PUT OUT SYMMETRY TRANSFORMATIONS 390 READ (NI,68)((TS(I,J),(FS(K,I,J),K=1,3),I=1,3),J=1,NS) 68 FORMAT (3(F15.10,3F3.0)) WRITE(NO,70) 70 FORMAT (' SYMMETRY INFORMATION'/' TRANSFORMED X' 1,20X,'TRANSFORMED Y',20X,'TRANSFORMED Z'/1X) WRITE(NO,72)((TS(I,J),(FS(K,I,J),K=1,3),I=1,3),J=1, 1NS) 72 FORMAT((1X,3(5X,F7.4,F5.1,2H*X,F5.1,2H*Y,F5.1,2H*Z))) 395 IF(IPM)415,425,400 C COMPUTE PARAMETER INCREMENTS USED TO OBTAIN DERIVATIVES 400 K=1 L=NV DO 410 I=1,NP IF(KI1(I))405,410,405 405 DP(I)= SQRT(PM(K)) K=K+L L=L-1 410 CONTINUE C PUT OUT INPUT PARAMETERS 415 WRITE(NO,74) 74 FORMAT (' INPUT DATA'/' I P(I) KI1(I) SIGMA(I)'/1H ) DO 420 I=1,NP WRITE(NO,76)I,P(I),KI1(I),DP(I) 76 FORMAT (1X,I3,1X,F11.6,I4,3X,F11.6) 420 DP(I)=(0.01D0)*DP(I) 425 WRITE(NO,78)(TITLE(I),I=1,18) 78 FORMAT (1 18A4/1H ) C READ ONE SET OF INSTRUCTIONS 430 K=1 DO 435 I=1,20 L=K+13 READ (NI,80)(IN(J),J=K,L) 80 FORMAT (14I5) IF(IN(L))435,440,435 435 K=L 440 IF(IN(1))445,455,460 445 IF(IN(1)+1)450,200,450 450 RETURN C END OF JOB 455 RETURN 460 IF(IN(1)-INSAVE)465,470,465 C PUT OUT HEADING FOR NEW TYPE OF FUNCTION 465 CALL HEDI(IN(1)) 470 INSAVE=IN(1) IF(IN(1)-100)475,475,480 C COMPUTE SINGLE-VALUED FUNCTION 475 CALL SUB19 GO TO 430 C COMPUTE MULTIPLE-VALUED FUNCTION 480 CALL SUB21 GO TO 430 END FUNCTION ARCCOS(X) C ***** ARCCOS(X) IN DEGREES ***** REAL*8 ARCCOS,X,SINE IF(1.D0-DABS(X))200,205,205 200 X=DSIGN(1.D0,X) 205 SINE=DSQRT(1.D0-X**2) IF(X)210,215,220 210 ARCCOS=57.29577951D0*DATAN(SINE/X)+180.D0 RETURN 215 ARCCOS=90.D0 RETURN 220 ARCCOS=57.29577951D0*DATAN(SINE/X) RETURN END SUBROUTINE ATOM(I,Z) C ATOM COORDINATE SUBROUTINE DIMENSION I(2),X(3),Y(3),Z(3) COMMON A(6),AA(3,3),AM(21),BB(3,3),DA(6),DFDA(6),DFDP( 160) 1,DP( 160),E,E1,F,FX,P( 160),SAVEA,SAVEP ,VARA,VARP REAL*8 A,AA,AM,BB,DA,DFDA,DFDP,DP,CHEM REAL*8 E,E1,F,FX,P,SAVEA,SAVEP,VARA,VARP COMMON CD(20,2),CHEM(200),DMAX,FS(3,3,48),IAM,IH,IN(261) A,INSAVE,IPM,ITF,JB,JBP,JX,JXP,KD(20,4),KI1( 160),KI2( 160) B,KK,KKD,LNZ,NA,NG,NM,NP,NQ,NS,NSP,NV,PM(5),ROW(9) C,SCALE,TITLE(18),TS(3,48) D,LOCX(201),ITA(200),IGM(200),INCD,IEXT,ITO,IKE(9),NPX,NI,NJ,NO REAL*8 X,Z IA=I(1) IF(IA) 210,200,205 200 X(1)=0.D0 X(2)=0.D0 X(3)=0.D0 GO TO 225 205 IF(IA-NA) 215,215,210 210 NG=5 GO TO 275 215 K=LOCX(IA) DO 220 J=1,3 X(J)=P(K) 220 K=K+1 225 KT=I(2)/100 KS=I(2)-100*KT IF(KS-NS) 235,235,230 230 NG=1 GO TO 275 235 IF(KS) 230,240,245 240 Z(1)=X(1) Z(2)=X(2) Z(3)=X(3) GO TO 255 245 DO 250 K=1,3 Z(K)=TS(K,KS) DO 250 J=1,3 250 Z(K)=Z(K)+FS(J,K,KS)*X(J) 255 IF(KT) 230,275,260 260 IF(KT-555) 270,265,270 265 I(2)=KS GO TO 275 270 K1=KT/100 K=KT-100*K1 K2=K/10 K3=K-10*K2 Z(1)=Z(1)+FLOAT(K1-5) Z(2)=Z(2)+FLOAT(K2-5) Z(3)=Z(3)+FLOAT(K3-5) 275 RETURN END SUBROUTINE AXES(U,V,X) C STORE THREE MUTUALLY PERPENDICULAR C VECTORS X(I,1), X(I,2), AND X(I,3) GIVEN C VECTORS U AND V. DIMENSIONU(3),V(3),X(3,3) COMMON A(6),AA(3,3),AM(21),BB(3,3),DA(6),DFDA(6),DFDP( 160) 1,DP( 160),E,E1,F,FX,P( 160),SAVEA,SAVEP ,VARA,VARP REAL*8 A,AA,AM,BB,DA,DFDA,DFDP,DP,CHEM REAL*8 E,E1,F,FX,P,SAVEA,SAVEP,VARA,VARP COMMON CD(20,2),CHEM(200),DMAX,FS(3,3,48),IAM,IH,IN(261) A,INSAVE,IPM,ITF,JB,JBP,JX,JXP,KD(20,4),KI1( 160),KI2( 160) B,KK,KKD,LNZ,NA,NG,NM,NP,NQ,NS,NSP,NV,PM( 5),ROW(9) C,SCALE,TITLE(18),TS(3,48) D,LOCX(201),ITA(200),IGM(200),INCD,IEXT,ITO,IKE(9),NPX,NI,NJ,NO REAL*8 T1,U,V,X,VMV DO200I=1,3 200 X(I,1)=U(I) CALLNORM(U,V,X(1,2)) CALLNORM(X(1,1),X(1,2),X(1,3)) DO 205 I=1,3 T1=DSQRT(VMV(X(1,I),AA,X(1,I))) DO 205 J=1,3 205 X(J,I)=X(J,I)/T1 RETURN END SUBROUTINE BETA(INS,Z) C STORE TRANSFORMED ANISOTROPIC TEMP FACTOR MATRIX C INS IS ATOM DESCRIPTION, Z IS TRANSFORMED MATRIX DIMENSION INS(2),Z(3 ,3),B1(3,3) COMMON A(6),AA(3,3),AM(21),BB(3,3),DA(6),DFDA(6),DFDP( 160) 1,DP( 160),E,E1,F,FX,P( 160),SAVEA,SAVEP ,VARA,VARP REAL*8 A,AA,AM,BB,DA,DFDA,DFDP,DP,CHEM REAL*8 E,E1,F,FX,P,SAVEA,SAVEP,VARA,VARP COMMON CD(20,2),CHEM(200),DMAX,FS(3,3,48),IAM,IH,IN(261) A,INSAVE,IPM,ITF,JB,JBP,JX,JXP,KD(20,4),KI1( 160),KI2( 160) B,KK,KKD,LNZ,NA,NG,NM,NP,NQ,NS,NSP,NV,PM( 5),ROW(9) C,SCALE,TITLE(18),TS(3,48) D,LOCX(201),ITA(200),IGM(200),INCD,IEXT,ITO,IKE(9),NPX,NI,NJ,NO REAL*8 B1,B2,Z IA=INS(1) IF(IA)215,200,210 200 DO 205 I=1,9 205 Z(I,1) =0.D0 GO TO 285 210 IF(IA-NA)220,220,215 215 NG=5 GO TO 285 220 IF(ITA(IA)-1)225,225,230 225 NG=4 GO TO 285 230 KT=INS(2)/100 KS=INS(2)-100*KT IF(KT-555)240,235,240 235 INS(2)=KS 240 IF(KS-NS)250,250,245 245 NG=1 GO TO 285 250 J=LOCX(IA)+3 255 B1(1,1)=P(J) B1(2,1)=P(J+3) B1(3,1)=P(J+4) B1(1,2)=P(J+3) B1(2,2)=P(J+1) B1(3,2)=P(J+5) B1(1,3)=P(J+4) B1(2,3)=P(J+5) B1(3,3)=P(J+2) IF(KS)245,260,270 260 DO 265 I=1,9 265 Z(I,1)=B1(I,1) GO TO 285 270 DO 280 I=1,3 DO 280 J=I,3 B2=0.D0 DO 275 K=1,3 DO 275 L=1,3 275 B2=B2+FS(K,I,KS)*FS(L,J,KS)*B1(K,L) Z(J,I)=B2 280 Z(I,J)=B2 285 RETURN END FUNCTION COSVV(X,Y) C COSINE OF ANGLE BETWEEN VECTORS X AND Y DIMENSIONX(3),Y(3) COMMON A(6),AA(3,3),AM(21),BB(3,3),DA(6),DFDA(6),DFDP( 160) 1,DP( 160),E,E1,F,FX,P( 160),SAVEA,SAVEP ,VARA,VARP REAL*8 A,AA,AM,BB,DA,DFDA,DFDP,DP,CHEM REAL*8 E,E1,F,FX,P,SAVEA,SAVEP,VARA,VARP COMMON CD(20,2),CHEM(200),DMAX,FS(3,3,48),IAM,IH,IN(261) A,INSAVE,IPM,ITF,JB,JBP,JX,JXP,KD(20,4),KI1( 160),KI2( 160) B,KK,KKD,LNZ,NA,NG,NM,NP,NQ,NS,NSP,NV,PM( 5),ROW(9) C,SCALE,TITLE(18),TS(3,48) D,LOCX(201),ITA(200),IGM(200),INCD,IEXT,ITO,IKE(9),NPX,NI,NJ,NO REAL*8 COSVV,X,Y,VMV,D D=DSQRT(VMV(X,AA,X)*VMV(Y,AA,Y)) IF(D)200,200,205 200 NG=9 GOTO210 205 COSVV=VMV(X,AA,Y)/D 210 RETURN END SUBROUTINE DIFV(X,Y,Z) C VECTOR - VECTOR C Z(3)=X(3)-Y(3) DIMENSIONX(3),Y(3),Z(3) REAL*8 X,Y,Z DO200I=1,3 200 Z(I)=X(I)-Y(I) RETURN END SUBROUTINE EIGVAL(W,Y) C ***** FIND EIGENVALUES Y OF MATRIX W ***** COMMON A(6),AA(3,3),AM(21),BB(3,3),DA(6),DFDA(6),DFDP( 160) 1,DP( 160),E,E1,F,FX,P( 160),SAVEA,SAVEP ,VARA,VARP REAL*8 A,AA,AM,BB,DA,DFDA,DFDP,DP,CHEM REAL*8 E,E1,F,FX,P,SAVEA,SAVEP,VARA,VARP COMMON CD(20,2),CHEM(200),DMAX,FS(3,3,48),IAM,IH,IN(261) A,INSAVE,IPM,ITF,JB,JBP,JX,JXP,KD(20,4),KI1( 160),KI2( 160) B,KK,KKD,LNZ,NA,NG,NM,NP,NQ,NS,NSP,NV,PM( 5),ROW(9) C,SCALE,TITLE(18),TS(3,48) D,LOCX(201),ITA(200),IGM(200),INCD,IEXT,ITO,IKE(9),NPX,NI,NJ,NO REAL*8 A1,A3,A27,B,B2,B4,C,HOLD,P1,P3,PHI3,Q,R,W,X,Y,Z,Z1 DIMENSION W(3,3),Y(3),X(3),Z(6,6) DO 200 J=1,3 DO 200 I=1,3 Z1=W(I,J) Z(I,J)=Z1 Z(I+3,J)=Z1 Z(I,J+3)=Z1 200 Z(I+3,J+3)=Z1 P1=0.D0 Q=0.D0 R=0.D0 DO 205 I=1,3 P1=P1-Z(I,I) Q=Q+Z(I,I)*Z(I+1,I+1)-Z(I,I+1)*Z(I+1,I) 205 R=R+Z(3,I)*Z(2,I+1)*Z(1,I+2)-Z(1,I)*Z(2,I+1)*Z(3,I+2) P3=P1/3.D0 A1=Q-P1*P3 B=2.D0*P3*P3*P3-Q*P3+R B2=B/2.D0 A3=A1/3.D0 B4=B2*B2 A27=A3*A3*A3 IF(B4+A27)225,225,210 210 IF(B4+1.0001D0*A27)220,220,215 215 NG=7 RETURN 220 A27=-B4 225 PHI3=(DATAN(DSQRT(-1.D0-A27/B4)))/3.D0 IF(B)235,230,235 230 B=0.D0 235 C=-DSIGN((2.D0*DSQRT(-A3)),B) X(1)=C*DCOS(PHI3) X(2)=C*DCOS(PHI3+4.188790205D0) X(3)=C*DCOS(PHI3+2.094395103D0) IF(B)240,245,245 240 HOLD=X(1) X(1)=X(3) X(3)=HOLD 245 DO 250 I=1,3 250 Y(I)=X(I)-P3 RETURN END SUBROUTINE EIGVEC(W,Y,Z) C COMPUTE EIGENVECTOR Z OF MATRIX C W GIVEN EIGENVALUE Y DIMENSIONW(3,3),X(6,6),Z(3),P1(3) COMMON A(6),AA(3,3),AM(21),BB(3,3),DA(6),DFDA(6),DFDP( 160) 1,DP( 160),E,E1,F,FX,P( 160),SAVEA,SAVEP ,VARA,VARP REAL*8 A,AA,AM,BB,DA,DFDA,DFDP,DP,CHEM REAL*8 E,E1,F,FX,P,SAVEA,SAVEP,VARA,VARP COMMON CD(20,2),CHEM(200),DMAX,FS(3,3,48),IAM,IH,IN(261) A,INSAVE,IPM,ITF,JB,JBP,JX,JXP,KD(20,4),KI1( 160),KI2( 160) B,KK,KKD,LNZ,NA,NG,NM,NP,NQ,NS,NSP,NV,PM( 5),ROW(9) C,SCALE,TITLE(18),TS(3,48) D,LOCX(201),ITA(200),IGM(200),INCD,IEXT,ITO,IKE(9),NPX,NI,NJ,NO REAL*8 P1,PJ,S,S1,W,X,X1,Y,Y1,Z DO200J=1,3 DO200I=1,3 X1=W(I,J) X(I,J)=X1 X(I+3,J)=X1 X(I,J+3)=X1 200 X(I+3,J+3)=X1 Y1=Y DO205I=1,3 X(I,I)=X(I,I)-Y1 X(I+3,I)=X(I+3,I)-Y1 X(I,I+3)=X(I,I+3)-Y1 205 X(I+3,I+3)=X(I+3,I+3)-Y1 S1=0.D0 DO225I=1,3 S=0.D0 DO210J=1,3 PJ=X(I,J+1)*X(I+1,J+2)-X(I,J+2)*X(I+1,J+1) P1(J)=PJ 210 S=S+PJ*PJ IF(S-S1)225,225,215 215 S1=S DO220J=1,3 220 Z(J)=P1(J) 225 CONTINUE IF(S1)230,230,235 230 NG=8 235 RETURN END FUNCTION FUNA(I) C ANGLE SUBROUTINE USED BY FUN2, FUN5, FUN6 DIMENSIONI(6),X1(3),X2(3),X3(3),V1(3),V2(3) COMMON A(6),AA(3,3),AM(21),BB(3,3),DA(6),DFDA(6),DFDP( 160) 1,DP( 160),E,E1,F,FX,P( 160),SAVEA,SAVEP ,VARA,VARP REAL*8 A,AA,AM,BB,DA,DFDA,DFDP,DP,CHEM REAL*8 E,E1,F,FX,P,SAVEA,SAVEP,VARA,VARP COMMON CD(20,2),CHEM(200),DMAX,FS(3,3,48),IAM,IH,IN(261) A,INSAVE,IPM,ITF,JB,JBP,JX,JXP,KD(20,4),KI1( 160),KI2( 160) B,KK,KKD,LNZ,NA,NG,NM,NP,NQ,NS,NSP,NV,PM( 5),ROW(9) C,SCALE,TITLE(18),TS(3,48) D,LOCX(201),ITA(200),IGM(200),INCD,IEXT,ITO,IKE(9),NPX,NI,NJ,NO REAL*8 FUNA,V1,V2,X1,X2,X3,ARCCOS,COSVV CALLATOM(I(1),X1) CALLATOM(I(3),X2) CALLATOM(I(5),X3) IF(NG)205,200,205 200 CALLDIFV(X1,X2,V1) CALLDIFV(X3,X2,V2) FUNA=ARCCOS(COSVV(V1,V2)) 205 RETURN END SUBROUTINE FUNB(W,Z,Z1) C SET UP MATRIX AND GET EIGENVALUE DIMENSIONB(3,3),W(3,3),Y(3) COMMON A(6),AA(3,3),AM(21),BB(3,3),DA(6),DFDA(6),DFDP( 160) 1,DP( 160),E,E1,F,FX,P( 160),SAVEA,SAVEP ,VARA,VARP REAL*8 A,AA,AM,BB,DA,DFDA,DFDP,DP,CHEM REAL*8 E,E1,F,FX,P,SAVEA,SAVEP,VARA,VARP COMMON CD(20,2),CHEM(200),DMAX,FS(3,3,48),IAM,IH,IN(261) A,INSAVE,IPM,ITF,JB,JBP,JX,JXP,KD(20,4),KI1( 160),KI2( 160) B,KK,KKD,LNZ,NA,NG,NM,NP,NQ,NS,NSP,NV,PM( 5),ROW(9) C,SCALE,TITLE(18),TS(3,48) D,LOCX(201),ITA(200),IGM(200),INCD,IEXT,ITO,IKE(9),NPX,NI,NJ,NO REAL*8 B,W,Y,Z,Z1 CALLBETA(IN(2),B) IF(NG)205,200,205 200 CALLMM(B,AA,W) CALLEIGVAL(W,Y) I=IN(4) Z=Y(I) Z1=DSQRT(Z*0.0506605918D0) 205 RETURN END SUBROUTINE FUNC(C,Z) C COS ANGLE OF PRINCIPAL AXIS AND VECTOR DIMENSIONW(3,3),X1(3),X2(3),V1(3),V2(3) COMMON A(6),AA(3,3),AM(21),BB(3,3),DA(6),DFDA(6),DFDP( 160) 1,DP( 160),E,E1,F,FX,P( 160),SAVEA,SAVEP ,VARA,VARP REAL*8 A,AA,AM,BB,DA,DFDA,DFDP,DP,CHEM REAL*8 E,E1,F,FX,P,SAVEA,SAVEP,VARA,VARP COMMON CD(20,2),CHEM(200),DMAX,FS(3,3,48),IAM,IH,IN(261) A,INSAVE,IPM,ITF,JB,JBP,JX,JXP,KD(20,4),KI1( 160),KI2( 160) B,KK,KKD,LNZ,NA,NG,NM,NP,NQ,NS,NSP,NV,PM( 5),ROW(9) C,SCALE,TITLE(18),TS(3,48) D,LOCX(201),ITA(200),IGM(200),INCD,IEXT,ITO,IKE(9),NPX,NI,NJ,NO REAL*8 C,V1,V2,W,X1,X2,Y,Z,COSVV CALLFUNB(W,Y,Z) IF(NG)215,200,215 200 CALLEIGVEC(W,Y,V1) IF(NG)215,205,215 205 CALLATOM(IN(5),X1) CALLATOM(IN(7),X2) IF(NG)215,210,215 210 CALLDIFV(X2,X1,V2) C=COSVV(V1,V2) 215 RETURN END SUBROUTINE FUNCR(C,R) C COMPUTE QUANTITIES FOR MEAN BOND DISTANCE DIMENSIONC(2),X(3,2),V(3) COMMON A(6),AA(3,3),AM(21),BB(3,3),DA(6),DFDA(6),DFDP( 160) 1,DP( 160),E,E1,F,FX,P( 160),SAVEA,SAVEP ,VARA,VARP REAL*8 A,AA,AM,BB,DA,DFDA,DFDP,DP,CHEM REAL*8 E,E1,F,FX,P,SAVEA,SAVEP,VARA,VARP COMMON CD(20,2),CHEM(200),DMAX,FS(3,3,48),IAM,IH,IN(261) A,INSAVE,IPM,ITF,JB,JBP,JX,JXP,KD(20,4),KI1( 160),KI2( 160) B,KK,KKD,LNZ,NA,NG,NM,NP,NQ,NS,NSP,NV,PM( 5),ROW(9) C,SCALE,TITLE(18),TS(3,48) D,LOCX(201),ITA(200),IGM(200),INCD,IEXT,ITO,IKE(9),NPX,NI,NJ,NO REAL*8 C,R,RSQ,V,X,XISQ,VMV DO200I=1,2 CALLFUNR(IN(2*I),RSQ) CALLFUNXI(IN(2*I),IN(2),XISQ) C(I)=RSQ-XISQ 200 CALLATOM(IN(2*I),X(1,I)) CALLDIFV(X(1,2),X(1,1),V) R=DSQRT(VMV(V,AA,V)) RETURN END FUNCTION FUND(I) C DISTANCE SUBROUTINE USED BY FUN1 AND FUN4 DIMENSIONI(4),X1(3),X2(3),V(3) COMMON A(6),AA(3,3),AM(21),BB(3,3),DA(6),DFDA(6),DFDP( 160) 1,DP( 160),E,E1,F,FX,P( 160),SAVEA,SAVEP ,VARA,VARP REAL*8 A,AA,AM,BB,DA,DFDA,DFDP,DP,CHEM REAL*8 E,E1,F,FX,P,SAVEA,SAVEP,VARA,VARP COMMON CD(20,2),CHEM(200),DMAX,FS(3,3,48),IAM,IH,IN(261) A,INSAVE,IPM,ITF,JB,JBP,JX,JXP,KD(20,4),KI1( 160),KI2( 160) B,KK,KKD,LNZ,NA,NG,NM,NP,NQ,NS,NSP,NV,PM( 5),ROW(9) C,SCALE,TITLE(18),TS(3,48) D,LOCX(201),ITA(200),IGM(200),INCD,IEXT,ITO,IKE(9),NPX,NI,NJ,NO REAL*8 FUND,V,X1,X2,VMV CALLATOM(I(1),X1) CALLATOM(I(3),X2) CALLDIFV(X2,X1,V) FUND=DSQRT(VMV(V,AA,V)) RETURN END SUBROUTINE FUNI(I) C SELECT THE FUN SUBROUTINE TO BE ENTERED DIMENSION X(3,6),V1(3),V2(3),V3(3),V4(3),V5(3),V6(3),CC(2),W(3,3) COMMON A(6),AA(3,3),AM(21),BB(3,3),DA(6),DFDA(6),DFDP( 160) 1,DP( 160),E,E1,F,FX,P( 160),SAVEA,SAVEP ,VARA,VARP REAL*8 A,AA,AM,BB,DA,DFDA,DFDP,DP,CHEM REAL*8 E,E1,F,FX,P,SAVEA,SAVEP,VARA,VARP COMMON CD(20,2),CHEM(200),DMAX,FS(3,3,48),IAM,IH,IN(261) A,INSAVE,IPM,ITF,JB,JBP,JX,JXP,KD(20,4),KI1( 160),KI2( 160) B,KK,KKD,LNZ,NA,NG,NM,NP,NQ,NS,NSP,NV,PM( 5),ROW(9) C,SCALE,TITLE(18),TS(3,48) D,LOCX(201),ITA(200),IGM(200),INCD,IEXT,ITO,IKE(9),NPX,NI,NJ,NO REAL*8 C,CC,R,RSQ,V1,V2,V3,V4,V5,V6,W,X,XISQ,Z REAL*8 FUND,FUNA,ARCCOS,COSVV,VMV IF(I)205,205,200 200 IF(I-15)210,210,310 205 NG=11 GO TO 315 210 GO TO (215,220,225,245,250,255,265,270,275,280,285,290,295,300,305 1),I C INTERATOMIC DISTANCE 215 FX=FUND(IN(2)) GO TO 315 C INTERATOMIC ANGLE 220 FX=FUNA(IN(2)) GO TO 315 C DIHEDRAL ANGLE 225 IF(NG)315,230,315 230 DO 235 J=1,6 235 CALL ATOM(IN(2*J),X(1,J)) IF(NG)315,240,315 240 CALL DIFV(X(1,2),X(1,1),V1) CALL DIFV(X(1,3),X(1,1),V2) CALL DIFV(X(1,5),X(1,4),V3) CALL DIFV(X(1,6),X(1,4),V4) CALL NORM(V1,V2,V5) CALL NORM(V3,V4,V6) FX=DSIGN(ARCCOS(COSVV(V5,V6)),VMV(V4,AA,V5)) GO TO 315 C DIFFERENCE BETWEEN INTERATOMIC DISTANCES 245 FX=FUND(IN(2))-FUND(IN(6)) GO TO 315 C DIFFERENCE BETWEEN INTERATOMIC ANGLES 250 FX=FUNA(IN(8))-FUNA(IN(2)) GO TO 315 C SUM OF ANGLES 255 K=IN(2) FX=0.D0 DO 260 J=1,K 260 FX=FX+FUNA(IN(6*J-3)) GO TO 315 C RMS PRINCIPAL DISPLACEMENT 265 CALL FUNB(W,Z,FX) GO TO 315 C ANGLE BETWEEN PRINCIPAL AXIS AND VECTOR 270 CALL FUNC(C,Z) FX=ARCCOS(C) GO TO 315 C PRINCIPAL AXIS PROJECTED ON VECTOR 275 CALL FUNC(C,Z) FX=C*Z GO TO 315 C ANGLE BETWEEN PRINCIPAL AND CARTESIAN AXIS 280 CALL FUNX(C,Z) FX=ARCCOS(C) GO TO 315 C PRINCIPAL AXIS PROJECTED ON CARTESIAN AXIS 285 CALL FUNX(C,Z) FX=C*Z GO TO 315 C RMS DISPLACEMENT IN GIVEN DIRECTION 290 CALL FUNXI(IN(2),IN(4),XISQ) FX=DSQRT(XISQ) GO TO 315 C RMS RADIAL DISPLACEMENT 295 CALL FUNR(IN(2),RSQ) FX=DSQRT(RSQ) GO TO 315 C MEAN BOND DISTANCE ASSUMING RIDING 300 CALL FUNCR(CC,R) FX=R+(CC(2)-CC(1))/(2.D0*R) GO TO 315 C MEAN INTERATOMIC DISTANCE ASSUMING INDEPENDENT MOTION 305 CALL FUNCR(CC,R) FX=R+(CC(2)+CC(1))/(2.D0*R) GO TO 315 310 CALL SPARE(I,3) 315 RETURN END SUBROUTINE FUNR(I,RSQ) C MEAN SQUARE RADIAL DISPLACEMENT DIMENSIONB(3,3),BAA(3,3),I(2) COMMON A(6),AA(3,3),AM(21),BB(3,3),DA(6),DFDA(6),DFDP( 160) 1,DP( 160),E,E1,F,FX,P( 160),SAVEA,SAVEP ,VARA,VARP REAL*8 A,AA,AM,BB,DA,DFDA,DFDP,DP,CHEM REAL*8 E,E1,F,FX,P,SAVEA,SAVEP,VARA,VARP COMMON CD(20,2),CHEM(200),DMAX,FS(3,3,48),IAM,IH,IN(261) A,INSAVE,IPM,ITF,JB,JBP,JX,JXP,KD(20,4),KI1( 160),KI2( 160) B,KK,KKD,LNZ,NA,NG,NM,NP,NQ,NS,NSP,NV,PM( 5),ROW(9) C,SCALE,TITLE(18),TS(3,48) D,LOCX(201),ITA(200),IGM(200),INCD,IEXT,ITO,IKE(9),NPX,NI,NJ,NO REAL*8 B,BAA,RSQ,TRACE CALLBETA(I,B) CALLMM(B,AA,BAA) RSQ=TRACE(BAA)*0.0506605918D0 RETURN END SUBROUTINE FUNX(C,Z) C COS ANGLE OF PRINCIPAL AND CARTESIAN AXES DIMENSIONW(3,3),V(3),X(3,4),V1(3),V2(3),AX(3,3) COMMON A(6),AA(3,3),AM(21),BB(3,3),DA(6),DFDA(6),DFDP( 160) 1,DP( 160),E,E1,F,FX,P( 160),SAVEA,SAVEP ,VARA,VARP REAL*8 A,AA,AM,BB,DA,DFDA,DFDP,DP,CHEM REAL*8 E,E1,F,FX,P,SAVEA,SAVEP,VARA,VARP COMMON CD(20,2),CHEM(200),DMAX,FS(3,3,48),IAM,IH,IN(261) A,INSAVE,IPM,ITF,JB,JBP,JX,JXP,KD(20,4),KI1( 160),KI2( 160) B,KK,KKD,LNZ,NA,NG,NM,NP,NQ,NS,NSP,NV,PM( 5),ROW(9) C,SCALE,TITLE(18),TS(3,48) D,LOCX(201),ITA(200),IGM(200),INCD,IEXT,ITO,IKE(9),NPX,NI,NJ,NO REAL*8 AX,C,V,V1,V2,W,X,Y,Z,COSVV CALLFUNB(W,Y,Z) IF(NG)220,200,220 200 CALLEIGVEC(W,Y,V) IF(NG)220,205,220 205 DO210I=1,4 210 CALLATOM(IN(2*I+4),X(1,I)) IF(NG)220,215,220 215 CALLDIFV(X(1,2),X(1,1),V1) CALLDIFV(X(1,4),X(1,3),V2) CALLAXES(V1,V2,AX) I=IN(5) C=COSVV(V,AX(1,I)) 220 RETURN END SUBROUTINE FUNXI(I,J,XISQ) C MEAN SQUARE DISPLACEMENT IN GIVEN DIRECTION DIMENSIONI(2),J(4),B(3,3),X1(3),X2(3),V(3),BAA(3,3),AABAA(3,3) COMMON A(6),AA(3,3),AM(21),BB(3,3),DA(6),DFDA(6),DFDP( 160) 1,DP( 160),E,E1,F,FX,P( 160),SAVEA,SAVEP ,VARA,VARP REAL*8 A,AA,AM,BB,DA,DFDA,DFDP,DP,CHEM REAL*8 E,E1,F,FX,P,SAVEA,SAVEP,VARA,VARP COMMON CD(20,2),CHEM(200),DMAX,FS(3,3,48),IAM,IH,IN(261) A,INSAVE,IPM,ITF,JB,JBP,JX,JXP,KD(20,4),KI1( 160),KI2( 160) B,KK,KKD,LNZ,NA,NG,NM,NP,NQ,NS,NSP,NV,PM( 5),ROW(9) C,SCALE,TITLE(18),TS(3,48) D,LOCX(201),ITA(200),IGM(200),INCD,IEXT,ITO,IKE(9),NPX,NI,NJ,NO REAL*8 AABAA,B,BAA,D,V,X1,X2,XISQ,VMV CALLBETA(I,B) CALLATOM(J,X1) CALLATOM(J(3),X2) IF(NG)215,200,215 200 CALLDIFV(X2,X1,V) D=VMV(V,AA,V) IF(D)205,205,210 205 NG=10 GOTO215 210 CALLMM(B,AA,BAA) CALLMM(AA,BAA,AABAA) XISQ=VMV(V,AABAA,V)*0.0506605918D0/D 215 RETURN END SUBROUTINE HEDI(I) C HEADER SUBROUTINE COMMON A(6),AA(3,3),AM(21),BB(3,3),DA(6),DFDA(6),DFDP( 160) 1,DP( 160),E,E1,F,FX,P( 160),SAVEA,SAVEP ,VARA,VARP REAL*8 A,AA,AM,BB,DA,DFDA,DFDP,DP,CHEM REAL*8 E,E1,F,FX,P,SAVEA,SAVEP,VARA,VARP COMMON CD(20,2),CHEM(200),DMAX,FS(3,3,48),IAM,IH,IN(261) A,INSAVE,IPM,ITF,JB,JBP,JX,JXP,KD(20,4),KI1( 160),KI2( 160) B,KK,KKD,LNZ,NA,NG,NM,NP,NQ,NS,NSP,NV,PM( 5),ROW(9) C,SCALE,TITLE(18),TS(3,48) D,LOCX(201),ITA(200),IGM(200),INCD,IEXT,ITO,IKE(9),NPX,NI,NJ,NO K=MOD(I,100) IF(K)290,290,200 200 IF(K-15)205,205,285 205 GO TO (210,215,220,225,230,235,240,245,250,255,260,265,270,275,280 1),K 2 FORMAT('INTERATOMIC DISTANCE IN ANGSTROMS') 210 WRITE(NO,2) GO TO 290 4 FORMAT('BOND ANGLE IN DEGREES. CENTRAL ATOM IS VERTEX') 215 WRITE(NO,4) GO TO 290 6 FORMAT('DIHEDRAL ANGLE BETWEEN PLANES EACH DEFINED BY THREE ATO 1MS') 220 WRITE(NO,6) GO TO 290 8 FORMAT('DIFFERENCE BETWEEN TWO INTERATOMIC DISTANCES') 225 WRITE(NO,8) GO TO 290 10 FORMAT('DIFFERENCE BETWEEN TWO BOND ANGLES') 230 WRITE(NO,10) GO TO 290 12 FORMAT('SUM OF SEVERAL BOND ANGLES') 235 WRITE(NO,12) GO TO 290 14 FORMAT('RMS COMPONENT OF THERMAL DISPLACEMENT ALONG PRINCIPAL A 1XIS R. ANGSTROMS'/11X,'ATOM',11X,'R') 240 WRITE(NO,14) GO TO 290 16 FORMAT('ANGLE BETWEEN PRINCIPAL AXIS R AND VECTOR DEFINED BY TW 1O ATOMS'/11X,'ATOM',11X,'R',16X,'VECTOR') 245 WRITE(NO,16) GO TO 290 18 FORMAT('RMS COMPONENT OF THERMAL DISPLACEMENT ALONG PRINCIPAL 1AXIS R PROJECTED ON VECTOR DEFINED BY TWO ATOMS. ANGSTROMS'/ 2 11X,'ATOM',11X,'R',16X,'VECTOR') 250 WRITE(NO,18) GO TO 290 20 FORMAT('ANGLE BETWEEN PRINCIPAL AXIS R AND AXIS I OF CARTESIAN 1SYSTEM DEFINED BY TWO VECTORS'/ 2 11X,'ATOM',11X,'R I',8X,'DEFINING VECTORS') 255 WRITE(NO,20) GO TO 290 22 FORMAT('RMS COMPONENT OF THERMAL DISPLACEMENT ALONG PRINCIPAL 1AXIS R PROJECTED ON AXIS I OF CARTESIAN SYSTEM'/' DEFINED BY TWO 2VECTORS. ANGSTROMS'/11X,'ATOM',11X,'R I',8X,'DEFINING VECTORS') 260 WRITE(NO,22) GO TO 290 24 FORMAT('RMS COMPONENT OF THERMAL DISPLACEMENT IN DIRECTION DEFI 1NED BY TWO ATOMS. ANGSTROMS'/11X,'ATOM',28X,'VECTOR') 265 WRITE(NO,24) GO TO 290 26 FORMAT('RMS RADIAL THERMAL DISPLACEMENT OF ATOM. ANGSTROMS') 270 WRITE(NO,26) GO TO 290 28 FORMAT('INTERATOMIC DISTANCE AVERAGED OVER THERMAL MOTION. SECO 1ND ATOM ASSUMED TO RIDE ON FIRST') 275 WRITE(NO,28) GO TO 290 30 FORMAT('INTERATOMIC DISTANCE AVERAGED OVER THERMAL MOTION. ATOM 1S ASSUMED TO MOVE INDEPENDENTLY') 280 WRITE(NO,30) GO TO 290 285 CALL SPARE(I,1) 290 RETURN END SUBROUTINE MM(X,Y,Z) C MULTIPLY TWO MATRICES C Z(3,3)=X(3,3)*Y(3,3) DIMENSIONX(3,3),Y(3,3),Z(3,3) REAL*8 X,Y,Z DO200I=1,3 DO200K=1,3 Z(I,K)=0.D0 DO200J=1,3 200 Z(I,K)=Z(I,K)+X(I,J)*Y(J,K) RETURN END SUBROUTINE MV(X,Y,Z) C MATRIX * VECTOR C Z(3)=X(3,3)*Y(3) DIMENSIONX(3,3),Y(3),Z(3) REAL*8 X,Y,Z DO200I=1,3 Z(I)=0.D0 DO200J=1,3 200 Z(I)=Z(I)+X(I,J)*Y(J) RETURN END SUBROUTINE NORM(X,Y,Z) C STORE A VECTOR Z NORMAL TO VECTORS X AND Y DIMENSIONX(3),Y(3),Z(3),X1(6),Y1(6),Z1(3) COMMON A(6),AA(3,3),AM(21),BB(3,3),DA(6),DFDA(6),DFDP( 160) 1,DP( 160),E,E1,F,FX,P( 160),SAVEA,SAVEP ,VARA,VARP REAL*8 A,AA,AM,BB,DA,DFDA,DFDP,DP,CHEM REAL*8 E,E1,F,FX,P,SAVEA,SAVEP,VARA,VARP COMMON CD(20,2),CHEM(200),DMAX,FS(3,3,48),IAM,IH,IN(261) A,INSAVE,IPM,ITF,JB,JBP,JX,JXP,KD(20,4),KI1( 160),KI2( 160) B,KK,KKD,LNZ,NA,NG,NM,NP,NQ,NS,NSP,NV,PM( 5),ROW(9) C,SCALE,TITLE(18),TS(3,48) D,LOCX(201),ITA(200),IGM(200),INCD,IEXT,ITO,IKE(9),NPX,NI,NJ,NO EQUIVALENCE (IN(1),IN1),(IN(2),IN2),(IN(3),IN3),(IN(4),IN4), 1 (IN(5),IN5),(IN(6),IN6),(IN(7),IN7),(IN(8),IN8), 1 (IN(9),IN9),(IN(10),IN10),(IN(11),IN11),(IN(12),IN12) REAL*8 X,X1,Y,Y1,Z,Z1 DO200I=1,3 X1(I)=X(I) X1(I+3)=X(I) Y1(I)=Y(I) 200 Y1(I+3)=Y(I) DO205I=1,3 205 Z1(I)=X1(I+1)*Y1(I+2)-X1(I+2)*Y1(I+1) CALLMV(BB,Z1,Z) RETURN END SUBROUTINE OUTI(I) C OUTPUT SUBROUTINE COMMON A(6),AA(3,3),AM(21),BB(3,3),DA(6),DFDA(6),DFDP( 160) 1,DP( 160),E,E1,F,FX,P( 160),SAVEA,SAVEP ,VARA,VARP REAL*8 A,AA,AM,BB,DA,DFDA,DFDP,DP,CHEM REAL*8 E,E1,F,FX,P,SAVEA,SAVEP,VARA,VARP COMMON CD(20,2),CHEM(200),DMAX,FS(3,3,48),IAM,IH,IN(261) A,INSAVE,IPM,ITF,JB,JBP,JX,JXP,KD(20,4),KI1( 160),KI2( 160) B,KK,KKD,LNZ,NA,NG,NM,NP,NQ,NS,NSP,NV,PM( 5),ROW(9) C,SCALE,TITLE(18),TS(3,48) D,LOCX(201),ITA(200),IGM(200),INCD,IEXT,ITO,IKE(9),NPX,NI,NJ,NO EQUIVALENCE (IN(1),IN1),(IN(2),IN2),(IN(3),IN3),(IN(4),IN4), 1 (IN(5),IN5),(IN(6),IN6),(IN(7),IN7),(IN(8),IN8), 1 (IN(9),IN9),(IN(10),IN10),(IN(11),IN11),(IN(12),IN12) 200 IF(I)275,275,205 205 IF(I-15)210,210,270 210 GO TO (215,220,225,230,225,235,245,250,250,255,255,260,265,215,215 1),I 215 WRITE(NO,2) CHEM(IN2),CHEM(IN4),(IN(J),J=2,5) 2 FORMAT( 5X,2(A6,1X),15X,2('(',I3,',',I5,') ')) GO TO 275 220 WRITE(NO,4) CHEM(IN2),CHEM(IN4),CHEM(IN6),(IN(J),J=2,7) 4 FORMAT( 5X,3(A6,1X), 8X,3('(',I3,',',I5,') ')) GO TO 275 225 WRITE(NO,6) CHEM(IN2),CHEM(IN4),CHEM(IN6),(IN(J),J=2,7) 1,CHEM(IN8),CHEM(IN10),CHEM(IN12),(IN(J),J=8,13) 6 FORMAT( 5X,3(A6,1X), 8X,3('(',I3,',',I5,') ')/ 1 6X,3(A6,1X), 8X,3('(',I3,',',I5,') ')) GO TO 275 230 WRITE(NO,8) CHEM(IN2),CHEM(IN4),(IN(J),J=2,5) 1, CHEM(IN6),CHEM(IN8),(IN(J),J=6,9) 8 FORMAT( 5X,2(A6,1X),15X,2('(',I3,',',I5,') ')/ 1 6X,2(A6,1X),15X,2('(',I3,',',I5,') ')) GO TO 275 235 WRITE(NO,10) 10 FORMAT(1X) NANG=IN(2) L=2 DO 240 J=1,NANG K=L+1 L=L+6 INK=IN(K) INK2=IN(K+2) INK4=IN(K+4) WRITE(NO,12) CHEM(INK),CHEM(INK2),CHEM(INK4),(IN(M),M=K,L) 12 FORMAT( 6X,3(A6,1X), 8X,3('(',I3,',',I5,') ')) 240 CONTINUE GO TO 275 245 WRITE(NO,14) CHEM(IN2),(IN(J),J=2,4) 14 FORMAT( 5X,A6, 1X,'(',I3,',',I5,')',I3) GO TO 275 250 WRITE(NO,16) CHEM(IN2),(IN(J),J=2,4),CHEM(IN 5 ),CHEM(IN 7 ) 1,(IN(J),J=5,8) 16 FORMAT( 5X,A6, 1X,'(',I3,',',I5,')', I3,6X,2(A6,1X) 1,2('(',I3,',',I5,') ')) GO TO 275 255 WRITE(NO,18) CHEM(IN2),(IN(J),J=2,5),CHEM(IN6),CHEM(IN8) 1,(IN(J),J=6,9),CHEM(IN10),CHEM(IN12),(IN(J),J=10,13) 18 FORMAT( 5X,A6, 1X,'(',I3,',',I5,')',2I3,3X,2(A6,1X) 1,2('(',I3,',',I5,') ')/33X,2(A6,1X),2('(',I3,',',I5,') ')) GO TO 275 260 WRITE(NO,20) CHEM(IN2),IN(2),IN(3),CHEM(IN4),CHEM(IN6) 1,(IN(J),J=4,7) 20 FORMAT( 5X,A6, 1X,'(',I3,',',I5,')', 9X,2(A6,1X) 1,2('(',I3,',',I5,') ')) GO TO 275 265 WRITE(NO,22)CHEM(IN2),IN(2),IN(3) 22 FORMAT( 5X,A6, 1X,'(',I3,',',I5,')') GO TO 275 270 CALL SPARE(I,4) 275 RETURN END SUBROUTINE OUTPUT(I) C SUBROUTINE FOR ADDITIONAL DATA OUTPUT RETURN END SUBROUTINE PREI(I) C SELECT THE PARAMETERS TO BE VARIED COMMON A(6),AA(3,3),AM(21),BB(3,3),DA(6),DFDA(6),DFDP( 160) 1,DP( 160),E,E1,F,FX,P( 160),SAVEA,SAVEP ,VARA,VARP REAL*8 A,AA,AM,BB,DA,DFDA,DFDP,DP,CHEM REAL*8 E,E1,F,FX,P,SAVEA,SAVEP,VARA,VARP COMMON CD(20,2),CHEM(200),DMAX,FS(3,3,48),IAM,IH,IN(261) A,INSAVE,IPM,ITF,JB,JBP,JX,JXP,KD(20,4),KI1( 160),KI2( 160) B,KK,KKD,LNZ,NA,NG,NM,NP,NQ,NS,NSP,NV,PM( 5),ROW(9) C,SCALE,TITLE(18),TS(3,48) D,LOCX(201),ITA(200),IGM(200),INCD,IEXT,ITO,IKE(9),NPX,NI,NJ,NO 200 IF(I)310,310,205 205 IF(I-15)210,210,305 210 GO TO (215,220,230,240,230,250,260,265,265,275,275,285,260,295,295 1),I 215 CALL SETKX(IN(2)) CALL SETKX(IN(4)) GO TO 310 220 DO 225 J=2,6,2 225 CALL SETKX(IN(J)) GO TO 310 230 DO 235 J=2,12,2 235 CALL SETKX(IN(J)) GO TO 310 240 DO 245 J=2,8,2 245 CALL SETKX(IN(J)) GO TO 310 250 K=IN(2)*6+1 DO 255 J=3,K,2 255 CALL SETKX(IN(J)) GO TO 310 260 CALL SETKB(IN(2)) GO TO 310 265 CALL SETKB(IN(2)) DO 270 J=5,7,2 270 CALL SETKX(IN(J)) GO TO 310 275 CALL SETKB(IN(2)) DO 280 J=6,12,2 280 CALL SETKX(IN(J)) GO TO 310 285 CALL SETKB(IN(2)) DO 290 J=4,6,2 290 CALL SETKX(IN(J)) GO TO 310 295 DO 300 J=2,4,2 CALL SETKX(IN(J)) 300 CALL SETKB(IN(J)) GO TO 310 305 CALL SPARE(I,2) 310 RETURN END SUBROUTINE SETKB(I) C SET KEY WORDS FOR ATOM BETAS C I=IN(K), THE INSTRUCTION INTEGER SPECIFYING THE ATOM NUMBER COMMON A(6),AA(3,3),AM(21),BB(3,3),DA(6),DFDA(6),DFDP( 160) 1,DP( 160),E,E1,F,FX,P( 160),SAVEA,SAVEP ,VARA,VARP REAL*8 A,AA,AM,BB,DA,DFDA,DFDP,DP,CHEM REAL*8 E,E1,F,FX,P,SAVEA,SAVEP,VARA,VARP COMMON CD(20,2),CHEM(200),DMAX,FS(3,3,48),IAM,IH,IN(261) A,INSAVE,IPM,ITF,JB,JBP,JX,JXP,KD(20,4),KI1( 160),KI2( 160) B,KK,KKD,LNZ,NA,NG,NM,NP,NQ,NS,NSP,NV,PM( 5),ROW(9) C,SCALE,TITLE(18),TS(3,48) D,LOCX(201),ITA(200),IGM(200),INCD,IEXT,ITO,IKE(9),NPX,NI,NJ,NO IF(I)210,210,200 200 J=LOCX(I)+3 DO 205 K=1,6 KI2(J)=1 205 J=J+1 210 RETURN END SUBROUTINE SETKX(I) C SELECT THE PRE SUBROUTINE TO BE ENTERED C SET KEY WORDS FOR ATOM COORDINATES C I=IN(K), THE INSTRUCTION INTEGER SPECIFYING THE ATOM NUMBER COMMON A(6),AA(3,3),AM(21),BB(3,3),DA(6),DFDA(6),DFDP( 160) 1,DP( 160),E,E1,F,FX,P( 160),SAVEA,SAVEP ,VARA,VARP REAL*8 A,AA,AM,BB,DA,DFDA,DFDP,DP,CHEM REAL*8 E,E1,F,FX,P,SAVEA,SAVEP,VARA,VARP COMMON CD(20,2),CHEM(200),DMAX,FS(3,3,48),IAM,IH,IN(261) A,INSAVE,IPM,ITF,JB,JBP,JX,JXP,KD(20,4),KI1( 160),KI2( 160) B,KK,KKD,LNZ,NA,NG,NM,NP,NQ,NS,NSP,NV,PM( 5),ROW(9) C,SCALE,TITLE(18),TS(3,48) D,LOCX(201),ITA(200),IGM(200),INCD,IEXT,ITO,IKE(9),NPX,NI,NJ,NO IF(I)205,205,200 200 J=LOCX(I) KI2(J)=1 KI2(J+1)=1 KI2(J+2)=1 205 RETURN END SUBROUTINE SPARE(I,J) C SPARE SUBROUTINE FOR NEW FUNCTIONS. C I IS IN(1), AN INTEGER WHICH DEFINES THE FUNCTION. C J INTERPRETED AS FOLLOWS- C J= 1 OUTPUT HEADING AS IN HEDI C J= 2 SPECIFY VARIABLES AS IN PREI C J= 3 PERFORM CALCULATIONS AS IN FUNI C J= 4 OUTPUT FUNCTION DESCRIPTION AS IN OUTI C FUNCTION 16 AS PROGRAMMED HERE WILL SERVE AS AN EXAMPLE. COMMON A(6),AA(3,3),AM(21),BB(3,3),DA(6),DFDA(6),DFDP( 160) 1,DP( 160),E,E1,F,FX,P( 160),SAVEA,SAVEP ,VARA,VARP REAL*8 A,AA,AM,BB,DA,DFDA,DFDP,DP,CHEM REAL*8 E,E1,F,FX,P,SAVEA,SAVEP,VARA,VARP COMMON CD(20,2),CHEM(200),DMAX,FS(3,3,48),IAM,IH,IN(261) A,INSAVE,IPM,ITF,JB,JBP,JX,JXP,KD(20,4),KI1( 160),KI2( 160) B,KK,KKD,LNZ,NA,NG,NM,NP,NQ,NS,NSP,NV,PM( 5),ROW(9) C,SCALE,TITLE(18),TS(3,48) D,LOCX(201),ITA(200),IGM(200),INCD,IEXT,ITO,IKE(9),NPX,NI,NJ,NO EQUIVALENCE (IN(1),IN1),(IN(2),IN2),(IN(3),IN3),(IN(4),IZ4), 1 (IN(5),IN5),(IN(6),IN6),(IN(7),IN7),(IN(8),IN8), 1 (IN(9),IN9),(IN(10),IN10),(IN(11),IN11),(IN(12),IN12) REAL*8 X,V1,V2,V3,V4,V5,D,VMV,COSVV,ARCCOS DIMENSION X(3,4),V1(3),V2(3),V3(3),V4(3),V5(3) DIMENSION AT(3,4),AX(3,3),FFP(3),V6(3),FFH(3),FFO(3),V7(3) DIMENSION IN4(8) DIMENSION FFPCOS(3),FFHCOS(3),FFOCOS(3) C REAL*8 QO,QH,QP REAL*8 AT,AX,FFP,V6,FFH,FFO,VNORM, V7 IF (I.EQ.22) GO TO 501 IF (I.EQ.23) GO TO 601 IF (I-21) 190,190,340 190 IF (I-20) 195,400,450 195 IF (I-18) 200,200,340 200 IF(I-17)205,260,300 C C DISTANCE OF ATOM 1 FROM PLANE DEFINED BY ATOMS 2, 3, AND 4. C ATOM 1 SPECIFIED BY IN(2) AND IN(3), ETC. C RIGHT-HAND THUMB POINTS IN POSITIVE DIRECTION IF FINGERS GO C THROUGH ATOMS 2, 3, AND 4 IN SUCCESSION. 205 GO TO (210,215,225,250),J 210 WRITE(NO,2) 2 FORMAT('DISTANCE OF ATOM FROM PLANE DEFINED BY THREE OTHERS'/ 1 11X,'ATOM',42X,'PLANE') GO TO 255 215 DO 220 K=2,8,2 CALL SETKX(IN(K)) 220 CONTINUE GO TO 255 225 DO 230 K=1,4 CALL ATOM(IN(2*K),X(1,K)) 230 CONTINUE IF(NG)255,235,255 235 CALL DIFV(X(1,1),X(1,2),V1) CALL DIFV(X(1,3),X(1,2),V2) CALL DIFV(X(1,4),X(1,2),V3) CALL NORM(V2,V3,V4) D=DSQRT(VMV(V4,AA,V4)) IF(D)240,240,245 240 NG=16 GO TO 255 245 FX=VMV(V1,AA,V4)/D GO TO 255 250 WRITE(NO,4)CHEM(IN2),IN(2),IN(3),CHEM(IZ4),IN(4),IN(5), 1 CHEM(IN6),IN(6),IN(7),CHEM(IN8),IN(8),IN(9) 4 FORMAT( 5X,A6,' (',I3,',',I5,')',28X,A6,' (',I3,',',I5,')'/ 1(52X,A6,' (',I3,',',I5,')')) 255 RETURN C C SIGNED CONFORMATION OR TORSION ANGLE 260 GO TO (265,270,280,295),J 265 WRITE(NO,6) 6 FORMAT('SIGNED CONFORMATION OR TORSION ANGLE FOR A CHAIN OF FOUR 1ATOMS') GO TO 340 270 DO 275 K=2,8,2 CALL SETKX(IN(K)) 275 CONTINUE GO TO 340 280 DO 285 K=1,4 285 CALL ATOM(IN(2*K),X(1,K)) IF(NG)340,290,340 290 CALL DIFV(X(1,2),X(1,1),V1) CALL DIFV(X(1,3),X(1,2),V2) CALL DIFV(X(1,4),X(1,3),V3) CALL NORM(V1,V2,V4) CALL NORM(V2,V3,V5) FX=DSIGN(ARCCOS(COSVV(V4,V5)),VMV(V3,AA,V4)) GO TO 340 295 WRITE(NO,8) CHEM(IN2),CHEM(IZ4),CHEM(IN6),CHEM(IN8),(IN(K),K=2,9) 8 FORMAT( 4(A6,1X),4('(',I3,',',I5,') ')) GO TO 340 C C ANGLE BETWEEN TWO VECTORS 300 GO TO (305,310,320,335),J 305 WRITE(NO,10) 10 FORMAT('ANGLE BETWEEN VECTORS EACH DEFINED BY TWO ATOMS') GO TO 340 310 DO 315 K=2,8,2 CALL SETKX(IN(K)) 315 CONTINUE GO TO 340 320 DO 325 K=1,4 325 CALL ATOM(IN(2*K),X(1,K)) IF(NG)340,330,340 330 CALL DIFV(X(1,2),X(1,1),V1) CALL DIFV(X(1,4),X(1,3),V2) FX=ARCCOS(COSVV(V1,V2)) GO TO 340 335 WRITE(NO,12) CHEM(IN2),CHEM(IZ4) ,(IN(K),K=2,5) 1, CHEM(IN6),CHEM(IN8),(IN(K),K=2,5) 12 FORMAT( 5X,2(A6,1X),15X,2('(',I3,',',I5,') ')/ 1 6X,2(A6,1X),15X,2('(',I3,',',I5,') ')) C C C TORSION ANGLE FOR NORMAL TO A PLANE OF ATOMS 1,2, AND 3 WITH THE C BOND BETWEEN ATOMS 3AND 4 400 GO TO (405,410,420,435),J 405 WRITE(NO,40) 40 FORMAT(' TORSION ANGLE FOR THE NORMAL TO THE PLANE OF THE FIRST T 1HREE ATOMS IN A CHAIN OF 4 ATOMS WITH THE BOND '/' BETWEEN ATOMS 3 2 AND 4') GO TO 340 410 DO 415 K=2,8,2 CALL SETKX(IN(K)) 415 CONTINUE GO TO 340 420 DO 425 K=1,4 425 CALL ATOM (IN(2*K),X(1,K)) IF(NG) 340,430,340 430 CALL DIFV(X(1,2),X(1,1),V1) CALL DIFV(X(1,3),X(1,2),V2) CALL DIFV(X(1,4),X(1,3),V3) CALL NORM(V1,V2,V4) CALL NORM(V2,V3,V5) CALL NORM(V5,V2,V6) FX=DSIGN(ARCCOS(COSVV(V4,V6)),VMV(V3,AA,V5)) GO TO 340 435 WRITE(NO,41)CHEM(IN2),CHEM(IZ4),CHEM(IN6),CHEM(IN8),(IN(K),K=2,9) 41 FORMAT( 4(A6,1X),4('(',I3,',',I5,') ')) GO TO 340 C C ANGLES BETWEEN NORMAL TO THE PLANE OF THREE ATOMS AND THE AXES C DEFINED BY THE FIRST FOUR ATOMS 1,2,3,& 4. WHERE AXES ARE 1-2, C 1-3, 1-4. 450 GO TO (455,460,470,490),J 455 WRITE (NO,45) 45 FORMAT('NORMAL TO THE PLANE OF ATOMS A1,A2,A3'/' IN PLANE BISECTO 1R OF THE ANGLE A1-A2-A3'/ ' NORMAL TO THE ABOVE TWO DIRECTIONS'/ 210X, ' ANGLES ARE IN REFERENCE TO AXES A,B,C WHICH ARE DEFINED 3BY THE FIRST FOUR NA ATOMS 1,2,3,&4'/ 4' A=1---2, B=1---3, C=1---4') GO TO 340 460 DO 465 K=2,6,2 CALL SETKX(IN(K)) 465 CONTINUE GO TO 340 470 DO 475 K=1,3 475 CALL ATOM(IN(2*K),X(1,K)) IF(NG) 340,480,340 480 CALL DIFV(X(1,1),X(1,2),V1) CALL DIFV(X(1,3),X(1,2),V2) CALL NORM (V1,V2,V3) VNORM= DSQRT(VMV(V1,AA,V1)) V1(1)=V1(1)/VNORM V1(2)=V1(2)/VNORM V1(3)=V1(3)/VNORM VNORM= DSQRT(VMV(V2,AA,V2)) V2(1)=V2(1)/VNORM V2(2)=V2(2)/VNORM V2(3)=V2(3)/VNORM CALL DIFV(V1(1),V2(1),V6) CALL VM(V6,AA,V7) CALL VM(V3,AA,V6) IF(V6(2).EQ.0.)GO TO 488 IF(V6(3).EQ.0.)GO TO 488 IF(V7(3).EQ.0.)GO TO 488 IF(V7(2).EQ.0.)GO TO 488 V4(1)=1.0 V4(2)=-(V6(1)/V6(3)-V7(1)/V7(3))/(V6(2)/V6(3)-V7(2)/V7(3)) V4(3)=-(V6(1)/V6(2)-V7(1)/V7(2))/(V6(3)/V6(2)-V7(3)/V7(2)) CALL NORM(V4,V3,V5) DO 483 K=1,4 IN4(2*K-1)=K 483 IN4(2*K)=0 DO 485 K=1,4 485 CALL ATOM (IN4(2*K-1), AT(1,K)) DO 486 K=2,4 486 CALL DIFV(AT(1,K),AT(1,1),AX(1,K-1)) DO 487 K=1,3 FFO(K)=ARCCOS(COSVV(AX(1,K),V5)) FFH(K)= ARCCOS(COSVV(AX(1,K),V4)) FFP(K)= ARCCOS(COSVV(AX(1,K),V3)) QO=FFO(K)*3.1415927D0/180. QH=FFH(K)*3.1415927D0/180. QP=FFP(K)*3.1415927D0/180. FFOCOS(K)=DCOS(QO) FFHCOS(K)=DCOS(QH) 487 FFPCOS(K)=DCOS(QP) NDIV0=0 GO TO 340 488 NDIV0=50 GO TO 340 490 IF(NDIV0.EQ.50) GO TO 491 GO TO 492 491 WRITE(NO,55) 55 FORMAT( T60,'ATTEMPTED DIVISION BY ZERO -DISREGARD ANSWERS') 492 WRITE(NO,44) 44 FORMAT(/ 'NORMAL TO THE PLANE',T79,'********************') LLL=1 WRITE(NO,56)CHEM(IN2),CHEM(IZ4),CHEM(IN6),(IN(K),K=2,7) DO 495 L=2,4 LL=L-1 495 WRITE(NO,46) CHEM(1),CHEM(L),LLL,L,FFP(LL) ,FFPCOS(LL) WRITE(NO,47) 47 FORMAT( 'IN THE PLANE BISECTOR OF THE ANGLE' ) WRITE(NO,56)CHEM(IN2),CHEM(IZ4),CHEM(IN6),(IN(K),K=2,7) DO 500 L=2,4 LL=L-1 500 WRITE(NO,46) CHEM(1),CHEM(L),LLL,L,FFH(LL),FFHCOS(LL) WRITE(NO,49) 49 FORMAT( 'NORMAL TO THE ABOVE TWO AXES' ) WRITE(NO,56)CHEM(IN2),CHEM(IZ4),CHEM(IN6),(IN(K),K=2,7) DO 502 L=2,4 LL=L-1 502 WRITE(NO,46) CHEM(1),CHEM(L),LLL,L,FFO(LL),FFOCOS(LL) 46 FORMAT( 5X,2(A6,1X),15X,2('(',I3,',',5X,') '),T79,F9.4,10X, 1F7.5) 56 FORMAT( 5X,3(A6,1X), 8X,3('(',I3,',',I5,')')) F=0.0 FX=0.0 C WRITE(NO,48) 48 FORMAT(/) GO TO 340 C C C ANGLES BETWEEN A VECTOR FROM A1 TO A2 AND THREE AXES DEFINED BY C THE FIRST FOUR NA ATOMS NA1, NA2, NA3, NA4. THE AXES ARE: C NA1-NA2, NA1-NA3, NA1-NA4. C 501 GO TO (505,510,520,540),J 505 WRITE(NO,50) 50 FORMAT('ANGLES BETWEEN A VECTOR BETWEEN TWO ATOMS AND THE AXES A, 1B,C'/' WHICH ARE DEFINED BY THE FIRST FOUR NA ATOMS 1, 2, 3, 4' 2' A=1---2, B=1---3, C=1---4') GO TO 340 510 DO 515 K=2,4,2 CALL SETKX(IN(K)) 515 CONTINUE GO TO 340 520 DO 525 K=1,2 525 CALL ATOM (IN(2*K),X(1,K)) IF(NG) 340,530,340 530 CALL DIFV(X(1,2),X(1,1),V1) DO 533 K=1,4 IN4(2*K-1)=K 533 IN4(2*K)=0 DO 535 K=1,4 535 CALL ATOM (IN4(2*K-1), AT(1,K)) DO 536 K=2,4 536 CALL DIFV(AT(1,K),AT(1,1),AX(1,K-1)) DO 537 K=1,3 KR=K+1- (K/3)*3 FFP(K)= ARCCOS(COSVV(AX(1,K),V1)) QP=FFP(K)*3.1415927D0/180. 537 FFPCOS(K)=DCOS(QP) F=0.0 FX=0.0 GO TO 340 540 WRITE(NO,57)CHEM(IN2),CHEM(IZ4),(IN(K),K=2,5) LLL=1 DO 545 L=2,4 LL=L-1 545 WRITE(NO,51) CHEM(1),CHEM(L),LLL,L,FFP(LL),FFPCOS(LL) 51 FORMAT( 5X,2(A6,1X),15X,2('(',I3,',',5X,') '),T79,F9.4,10X, 1F7.5) 57 FORMAT( T75,'******************'/6X,2(A6,1X),15X,2('(',I3,',',I5, 1')')) FX=0.0 C WRITE(NO,48) GO TO 340 601 CONTINUE C SUBROUTINE TO COMPUTE SYMMETRY RELATED ATOMS OF ORIGINAL XYZ COORDINATED C L = IN(5) GO TO (605,610,620,635),J 605 WRITE (NO,62) 62 FORMAT (//,2X,'SYMMETRY RELATED COORDINATE',/) WRITE (NO,608) 608 FORMAT ('********************') GO TO 340 610 DO 615 K=2,4,2 CALL SETKX(IN(K)) 615 CONTINUE GO TO 340 620 CALL ATOM (IN(4),X) WRITE (NO, 630) IN(2),IN(3),IN(4),L,X(1,1),X(2,1),X(3,1) 630 FORMAT (2X,4(I5,2X),3F10.5) F = 0.0 FX = 0.0 GO TO 340 635 CONTINUE F = 0.0 340 RETURN END SUBROUTINE STOAA C STORE METRIC TENSOR COMMON A(6),AA(3,3),AM(21),BB(3,3),DA(6),DFDA(6),DFDP( 160) 1,DP( 160),E,E1,F,FX,P( 160),SAVEA,SAVEP ,VARA,VARP REAL*8 A,AA,AM,BB,DA,DFDA,DFDP,DP,CHEM REAL*8 E,E1,F,FX,P,SAVEA,SAVEP,VARA,VARP COMMON CD(20,2),CHEM(200),DMAX,FS(3,3,48),IAM,IH,IN(261) A,INSAVE,IPM,ITF,JB,JBP,JX,JXP,KD(20,4),KI1( 160),KI2( 160) B,KK,KKD,LNZ,NA,NG,NM,NP,NQ,NS,NSP,NV,PM( 5),ROW(9) C,SCALE,TITLE(18),TS(3,48) D,LOCX(201),ITA(200),IGM(200),INCD,IEXT,ITO,IKE(9),NPX,NI,NJ,NO AA(1,1)=A(1)*A(1) AA(2,2)=A(2)*A(2) AA(3,3)=A(3)*A(3) AA(1,2)=A(1)*A(2)*A(6) AA(2,1)=AA(1,2) AA(1,3)=A(1)*A(3)*A(5) AA(3,1)=AA(1,3) AA(2,3)=A(2)*A(3)*A(4) AA(3,2)=AA(2,3) RETURN END SUBROUTINE STOBB C STORE RECIPROCAL METRIC TENSOR DIMENSIONAI(6),CI(6),BBII(3),BBJK(3) COMMON A(6),AA(3,3),AM(21),BB(3,3),DA(6),DFDA(6),DFDP( 160) 1,DP( 160),E,E1,F,FX,P( 160),SAVEA,SAVEP ,VARA,VARP REAL*8 A,AA,AM,BB,DA,DFDA,DFDP,DP,CHEM REAL*8 E,E1,F,FX,P,SAVEA,SAVEP,VARA,VARP COMMON CD(20,2),CHEM(200),DMAX,FS(3,3,48),IAM,IH,IN(261) A,INSAVE,IPM,ITF,JB,JBP,JX,JXP,KD(20,4),KI1( 160),KI2( 160) B,KK,KKD,LNZ,NA,NG,NM,NP,NQ,NS,NSP,NV,PM( 5),ROW(9) C,SCALE,TITLE(18),TS(3,48) D,LOCX(201),ITA(200),IGM(200),INCD,IEXT,ITO,IKE(9),NPX,NI,NJ,NO REAL*8 AI,CI,X,BBII,BBJK DO210I=1,3 IF(A(I))200,200,205 200 NG=6 GOTO220 205 AI(I)=A(I) AI(I+3)=A(I) CI(I)=A(I+3) 210 CI(I+3)=A(I+3) X=1.D0/(AI(1)*AI(2)*AI(3)*(1.D0-CI(1)*CI(1)-CI(2)*CI(2) 1-CI(3)*CI(3)+2.D0*CI(1)*CI(2)*CI(3))) DO215I=1,3 BBII(I)=X*(1.D0-CI(I)*CI(I))*AI(I+1)*AI(I+2)/AI(I) 215 BBJK(I)=X*AI(I)*(CI(I+1)*CI(I+2)-CI(I)) BB(1,1)=BBII(1) BB(1,2)=BBJK(3) BB(1,3)=BBJK(2) BB(2,1)=BBJK(3) BB(2,2)=BBII(2) BB(2,3)=BBJK(1) BB(3,1)=BBJK(2) BB(3,2)=BBJK(1) BB(3,3)=BBII(3) 220 RETURN END SUBROUTINE SUB10 C MULTIVALUED FUNCTIONS 7, 8, AND 9 COMMON A(6),AA(3,3),AM(21),BB(3,3),DA(6),DFDA(6),DFDP( 160) 1,DP( 160),E,E1,F,FX,P( 160),SAVEA,SAVEP ,VARA,VARP REAL*8 A,AA,AM,BB,DA,DFDA,DFDP,DP,CHEM REAL*8 E,E1,F,FX,P,SAVEA,SAVEP,VARA,VARP COMMON CD(20,2),CHEM(200),DMAX,FS(3,3,48),IAM,IH,IN(261) A,INSAVE,IPM,ITF,JB,JBP,JX,JXP,KD(20,4),KI1( 160),KI2( 160) B,KK,KKD,LNZ,NA,NG,NM,NP,NQ,NS,NSP,NV,PM( 5),ROW(9) C,SCALE,TITLE(18),TS(3,48) D,LOCX(201),ITA(200),IGM(200),INCD,IEXT,ITO,IKE(9),NPX,NI,NJ,NO C START LOOP THRU PRINCIPAL AXES DO200I=1,3 IN(4)=I C COMPUTE AND PUT OUT FUNCTION AND ERROR 200 CALLSUB19 RETURN END SUBROUTINE SUB11 C MULTIVALUED FUNCTIONS 10 AND 11 COMMON A(6),AA(3,3),AM(21),BB(3,3),DA(6),DFDA(6),DFDP( 160) 1,DP( 160),E,E1,F,FX,P( 160),SAVEA,SAVEP ,VARA,VARP REAL*8 A,AA,AM,BB,DA,DFDA,DFDP,DP,CHEM REAL*8 E,E1,F,FX,P,SAVEA,SAVEP,VARA,VARP COMMON CD(20,2),CHEM(200),DMAX,FS(3,3,48),IAM,IH,IN(261) A,INSAVE,IPM,ITF,JB,JBP,JX,JXP,KD(20,4),KI1( 160),KI2( 160) B,KK,KKD,LNZ,NA,NG,NM,NP,NQ,NS,NSP,NV,PM( 5),ROW(9) C,SCALE,TITLE(18),TS(3,48) D,LOCX(201),ITA(200),IGM(200),INCD,IEXT,ITO,IKE(9),NPX,NI,NJ,NO C START LOOP THRU PRINCIPAL AXES DO200I=1,3 IN(4)=I C START LOOP THRU REFERENCE AXES DO200J=1,3 IN(5)=J C COMPUTE AND PUT OUT FUNCTION AND ERROR 200 CALLSUB19 RETURN END SUBROUTINE SUB13 C ERROR CALCULATION AND OUTPUT COMMON A(6),AA(3,3),AM(21),BB(3,3),DA(6),DFDA(6),DFDP( 160) 1,DP( 160),E,E1,F,FX,P( 160),SAVEA,SAVEP ,VARA,VARP REAL*8 A,AA,AM,BB,DA,DFDA,DFDP,DP,CHEM REAL*8 E,E1,F,FX,P,SAVEA,SAVEP,VARA,VARP COMMON CD(20,2),CHEM(200),DMAX,FS(3,3,48),IAM,IH,IN(261) A,INSAVE,IPM,ITF,JB,JBP,JX,JXP,KD(20,4),KI1( 160),KI2( 160) B,KK,KKD,LNZ,NA,NG,NM,NP,NQ,NS,NSP,NV,PM( 5),ROW(9) C,SCALE,TITLE(18),TS(3,48) D,LOCX(201),ITA(200),IGM(200),INCD,IEXT,ITO,IKE(9),NPX,NI,NJ,NO REAL*8 C,TEMP C PUT OUT FUNCTION DESCRIPTION CALL OUTI(IN(1)) IF(NG)200,205,200 C PUT OUT ERROR INDICATOR IF NOT ZERO 200 WRITE(NO,2)NG 2 FORMAT ( 78X,'***',I3) GO TO 415 205 VARA=0.D0 VARP=0.D0 IF(IAM)210,260,210 C COMPUTE DERIVATIVES WITH RESPECT TO CELL PARAMETERS 210 DO 225 I=1,6 IF(DA(I))220,215,220 215 DFDA(I)=0.D0 GO TO 225 220 SAVEA=A(I) A(I)=A(I)+DA(I) CALL SETA(A) CALL STOAA CALL STOBB CALL FUNI(IN(1)) A(I)=SAVEA DFDA(I)=(FX-F)/DA(I) 225 CONTINUE CALL SETA(A) CALL STOAA CALL STOBB C COMPUTE VARIANCE BASED ON CELL PARAMETERS K=1 L=6 DO 255 I=1,6 IF(DFDA(I))235,230,235 230 K=K+L GO TO 255 235 C=1.D0 DO 250 J=I,6 IF(DFDA(J))240,245,240 240 VARA=VARA+C*DFDA(I)*DFDA(J)*AM(K) 245 K=K+1 250 C=2.D0 255 L=L-1 260 IF(IPM)265,370,265 C SELECT DERIVATIVES TO BE COMPUTED 265 DO 270 I=1,NP 270 KI2(I)=0 CALL PREI(IN(1)) C COMPUTE DERIVATIVES WITH RESPECT TO STRUCTURE PARAMETERS I=0 NUL=0 DO 305 J=1,NP IF(KI1(J))285,275,285 275 IF(KI2(J))280,305,280 280 NUL=NUL+1 GO TO 305 285 I=I+1 IF(KI2(J))295,290,295 290 DFDP(I)=0.D0 GO TO 305 295 IF(DP(J))290,290,300 300 SAVEP=P(J) P(J)=P(J)+DP(J) CALL SETP(P) CALL FUNI(IN(1)) P(J)=SAVEP DFDP(I)=(FX-F)/DP(J) LNZ=I 305 CONTINUE CALL SETP(P) C COMPUTE VARIANCE BASED ON STRUCTURE PARAMETERS IF(IPM)345,310,310 310 KK=1 KKD=NV DO 340 I=1,LNZ TEMP=DFDP(I) IF(TEMP)315,335,315 315 K=KK C=1.D0 DO 330 J=I,LNZ IF(DFDP(J))320,325,320 320 VARP=VARP+C*TEMP*DFDP(J)*PM(K) 325 K=K+1 330 C=2.D0 335 KK=KK+KKD 340 KKD=KKD-1 GO TO 370 C VARIANCE BASED ON DIAGONAL VARIANCE MATRIX 345 J=0 DO 365 I=1,LNZ 350 J=J+1 IF(KI1(J)) 355,350,355 355 TEMP=DFDP(I) IF(TEMP) 360,365,360 360 VARP=VARP+(TEMP*DP(J)*100.D0)**2 365 CONTINUE 370 IF(NG) 375,380,375 C PUT OUT ERROR INDICATOR IF NOT ZERO 375 WRITE(NO,4) F,NG 4 FORMAT ( 78X,F9.4,6X,'***',I3) GO TO 415 C COMPUTE STANDARD ERRORS AND PUT OUT RESULTS 380 E1=DSQRT(VARP) E=DSQRT(VARP+VARA) IF(IPM)385,390,385 385 IF(IAM)395,400,395 390 IF(IAM)400,405,400 395 WRITE(NO,6) F,E,E1,NUL 6 FORMAT( 78X,F9.4,' +OR-',F8.4,2H (F7.4,4H) *I2,2H *) GO TO 410 400 WRITE(NO,8) F,E,NUL 8 FORMAT( 78X,F9.4,6H +OR-F8.4,3H *I2,2H *) GO TO 410 405 CONTINUE IF (IN(1) .GT. 20) GO TO 410 WRITE(NO,10) F 10 FORMAT(78X,F9.4) 410 CALL OUTPUT(IN) 415 RETURN END SUBROUTINE SUB19 C FUNCTION AND ERROR CALCULATION COMMON A(6),AA(3,3),AM(21),BB(3,3),DA(6),DFDA(6),DFDP( 160) 1,DP( 160),E,E1,F,FX,P( 160),SAVEA,SAVEP ,VARA,VARP REAL*8 A,AA,AM,BB,DA,DFDA,DFDP,DP,CHEM REAL*8 E,E1,F,FX,P,SAVEA,SAVEP,VARA,VARP COMMON CD(20,2),CHEM(200),DMAX,FS(3,3,48),IAM,IH,IN(261) A,INSAVE,IPM,ITF,JB,JBP,JX,JXP,KD(20,4),KI1( 160),KI2( 160) B,KK,KKD,LNZ,NA,NG,NM,NP,NQ,NS,NSP,NV,PM( 5),ROW(9) C,SCALE,TITLE(18),TS(3,48) D,LOCX(201),ITA(200),IGM(200),INCD,IEXT,ITO,IKE(9),NPX,NI,NJ,NO NG=0 CALL FUNI(IN(1)) F=FX 200 CALL SUB13 205 RETURN END SUBROUTINE SUB21 C COMPUTE MULTIVALUED FUNCTIONS COMMON A(6),AA(3,3),AM(21),BB(3,3),DA(6),DFDA(6),DFDP( 160) 1,DP( 160),E,E1,F,FX,P( 160),SAVEA,SAVEP ,VARA,VARP REAL*8 A,AA,AM,BB,DA,DFDA,DFDP,DP,CHEM REAL*8 E,E1,F,FX,P,SAVEA,SAVEP,VARA,VARP COMMON CD(20,2),CHEM(200),DMAX,FS(3,3,48),IAM,IH,IN(261) A,INSAVE,IPM,ITF,JB,JBP,JX,JXP,KD(20,4),KI1( 160),KI2( 160) B,KK,KKD,LNZ,NA,NG,NM,NP,NQ,NS,NSP,NV,PM( 5),ROW(9) C,SCALE,TITLE(18),TS(3,48) D,LOCX(201),ITA(200),IGM(200),INCD,IEXT,ITO,IKE(9),NPX,NI,NJ,NO C UNPACK INSTRUCTION NUMBER IH=IN(1)/100 IN(1)=IN(1)-100*IH C TRANSFER TO APPROPRIATE SECTION IF(IN(1)-11)200,200,260 200 IF(IN(1)-9)205,205,240 205 IF(IN(1)-6)210,210,220 210 IF(IN(1)-1)215,215,260 215 CALL SUB 23 GOTO260 220 IF(IH-1)225,225,230 C COMPUTE FUNCTIONS INVOLVING THREE PRINCIPAL AXES 225 CALLSUB10 GOTO260 C COMPUTE FUNCTIONS INVOLVING THREE PRINCIPAL AXES C FOR ALL ATOMS 230 NA=IN(2) DO235I=1,NA IN(2)=I 235 CALLSUB10 GOTO260 240 IF(IH-1)245,245,250 C COMPUTE FOR ALL PRINCIPAL AXES AND ALL REFERENCE AXES 245 CALLSUB11 GOTO260 C COMPUTE FOR ALL PRINCIPAL AXES, ALL REFERENCE AXES, 250 NA=IN(2) DO255I=1,NA IN(2)=I 255 CALLSUB11 260 RETURN END SUBROUTINE SUB23 DIMENSION DX(3),IZ(2),KC(2),NW(6),D(2,300) DIMENSION U(3),V(3),W(2,4),WW(2,3),X(4),Y(3),Z(3),COMNT(6) DIMENSION AMINMX(4) COMMON A(6),AA(3,3),AM(21),BB(3,3),DA(6),DFDA(6),DFDP( 160) 1,DP( 160),E,E1,F,FX,P( 160),SAVEA,SAVEP ,VARA,VARP REAL*8 A,AA,AM,BB,DA,DFDA,DFDP,DP,CHEM REAL*8 E,E1,F,FX,P,SAVEA,SAVEP,VARA,VARP COMMON CD(20,2),CHEM(200),DMAX,FS(3,3,48),IAM,IH,IN(261) A,INSAVE,IPM,ITF,JB,JBP,JX,JXP,KD(20,4),KI1( 160),KI2( 160) B,KK,KKD,LNZ,NA,NG,NM,NP,NQ,NS,NSP,NV,PM( 5),ROW(9) C,SCALE,TITLE(18),TS(3,48) D,LOCX(201),ITA(200),IGM(200),INCD,IEXT,ITO,IKE(9),NPX,NI,NJ,NO REAL*8 U,V,X,Y,Z,ARCCOS,VMV,FUND,FUNA EQUIVALENCE(NW(1),LL),(NW(2),LU),(NW(3),ML) EQUIVALENCE(NW(4),MU),(NW(5),NL),(NW(6),NU) 2 FORMAT( 'FIND INTERATOMIC VECTORS') 4 FORMAT( ' FROM ATOM NO. TO ATOMS NOS. WITH LENGTHS' /' CO 1DE (MIN--MAX) (MIN--MAX) (MIN--MAX) *** COMMENTS ***') 6 FORMAT( 'SAVE VECTORS FOUND WHICH SATISIFY ONE OF FOLLOWING SELE 1CTION CODES') 8 FORMAT( 'NO FURTHER SELECTION') 10 FORMAT(5X,4I5,2F5.2,9A4) 12 FORMAT( I4,I8,I5,I8,I5,2F9.3) 14 FORMAT( I4,I8,I5,I8,I5,2F9.3,9A4) C ***** OBTAIN PROBLEM PARAMETERS ***** NG=0 ITOM1=IN(2) ITOM2=IN(3) ITAR1=IN(4) IF(ITAR1)200,200,205 200 ITAR1=1 205 ITAR2=IN(5) IF(IN(6))210,210,215 210 DMAX=4.0 GO TO 220 215 DMAX=FLOAT(IN(6))/100.0 220 KODE=IN(7) IOUT=IN(8) DO 225 I=1,4 225 AMINMX(I)=FLOAT(IN(I+8))/100. DMX=DMAX*DMAX WRITE(NO,2) WRITE(NO,4) TEM=0.01 II=IN(1)+100*IH WRITE(NO,12)II,ITOM1,ITOM2,ITAR1,ITAR2,TEM,DMAX IF(KODE)240,240,230 230 WRITE(NO,6) WRITE(NO,4) DO 235 I=1,KODE READ (NI,10)KD(I,1),KD(I,2),KD(I,3),KD(I,4),(CD(I,J),J=1 1,2),(COMNT(K),K=1,6) 235 WRITE(NO,14)I,KD(I,1),KD(I,2),KD(I,3),KD(I,4),(CD(I,J), 1J=1,2),(COMNT(K),K=1,6) GO TO 245 240 WRITE(NO,8) C ***** FIND RANGE OF X,Y,Z,X-Y ***** 245 DO 250 J=1,4 W(1,J)=99.0 250 W(2,J)=-99.0 IZ(2)=0 DO 275 I=ITAR1,ITAR2 IZ(1)=I CALL ATOM(IZ,X) IF(NG)255,255,610 255 X(4)=X(1)-X(2) DO 275 J=1,4 TEM=X(J) IF(W(2,J)-TEM)260,265,265 260 W(2,J)=TEM 265 IF(TEM-W(1,J))270,275,275 270 W(1,J)=TEM 275 CONTINUE C ***** FIND PARALLELEPIPED WHICH ENCLOSES DMAX SPHERE ***** TEM=1.0-A(4)*A(4)-A(5)*A(5)-A(6)*A(6)+2.0*A(4)*A(5)*A(6) DO 280 J=1,3 AJ3=A(J+3) AJ=A(J) 280 DX(J)=SQRT((1.0-AJ3**2)/TEM)*DMAX/AJ C ***** START SEARCH AROUND REFERENCE ATOMS ***** DO 605 ITOM=ITOM1,ITOM2 IZ(1)=ITOM IZ(2)=1 CALL ATOM(IZ,Y) C ***** K=SYMMETRY EQUIVALENT POSITION ***** NUM=0 285 DO 475 K=1,NS C ***** SUBTRACT SYMMETRY TRANSLATION FROM REF ATOM ***** DO 290 J=1,3 290 U(J)=Y(J)-TS(J,K) C ***** DETERMINE LIMITING CELLS TO BE SEARCHED ***** C ***** FIRST,MOVE THE BOX THROUGH THE SYMMETRY OPERATION ***** DO 310 J=1,3 DO 310 L=1,2 WW(L,J)=0.0 DO 310 I=1,3 TEM=FS(I,J,K) IF(TEM)295,310,300 295 N=MOD(L,2)+1 GO TO 305 300 N=L 305 WW(L,J)=WW(L,J)+W(N,I)*TEM 310 CONTINUE C ***** CHECK FOR MIXED INDEX TRANSFORMATION ***** DO 330 J=1,2 TEM=FS(1,J,K) IF(TEM+FS(2,J,K))330,315,330 315 IF(TEM)320,330,325 320 WW(1,J)=W(2,4)*TEM WW(2,J)=W(1,4)*TEM GO TO 330 325 WW(1,J)=W(1,4)*TEM WW(2,J)=W(2,4)*TEM 330 CONTINUE C ***** MOVE 4 CELLS AWAY THEN MOVE BACK UNTIL PARALLELEPIPED AROUND C REF ATOM AND BOX AROUND TRANSFORMED ASYM UNIT INTERSECT ***** N=0 DO 345 J=1,3 DO 340 I=1,2 N=N+1 TT=(U(J)-WW(I,J))*((-1.)**I)-DX(J) TEM=5.0 335 TEM=TEM-1.0 IF(TEM+TT)340,340,335 340 NW(N)=TEM*((-1.)**I)+5. C ***** IF NO POSSIBILITY OF A HIT, GO TO NEXT SYMMETRY OPER ***** IF(NW(N)-NW(N-1))475,345,345 345 CONTINUE C ***** L CELL TRANSLATIONS IN X ***** DO 470 L=LL,LU V(1)=U(1)+FLOAT(L-5) C ***** M CELL TRANSLATIONS IN Y ***** DO 470 M=ML,MU V(2)=U(2)+FLOAT(M-5) C ***** N CELL TRANSLATIONS IN Z ***** DO 470 NN=NL,NU V(3)=U(3)+FLOAT(NN-5) C ***** I = TARGET ATOM ***** DO 470 I=ITAR1,ITAR2 JJ=LOCX(I) DO 355 J=1,3 TEM=0.0 IJ=JJ DO 350 II=1,3 TEM=TEM+FS(II,J,K)*P(IJ) 350 IJ=IJ+1 C ***** SEE IF WITHIN PARALLELEPIPED***** TEM=TEM-V(J) IF(DX(J)-TEM)470,470,355 355 X(J)=TEM C ***** CALCULATE D SQ,SEE IF WITHIN SPHERE ***** DSQ=VMV(X,AA,X) IF(DMX-DSQ)470,360,360 360 IF(DSQ-.0001)470,365,365 C ***** SELECT VECTORS ACCORDING TO CODES IF ANY ***** 365 TEM=SQRT(DSQ) IF(KODE)410,410,370 370 DO 405 J=1,KODE 375 IF(ITOM-KD(J,1))405,380,380 380 IF(KD(J,2)-ITOM)405,385,385 385 IF(I-KD(J,3)) 405,390,390 390 IF(KD(J,4)-I) 405,395,395 395 IF(TEM-CD(J,1)) 405,400,400 400 IF(CD(J,2)-TEM) 405,410,410 405 CONTINUE GO TO 470 C ***** DETERMINE CORRECT POSITION IN SORTED VECTOR TABLE ***** 410 IF(NUM)460,460,415 415 DO 455 II=1,NUM TT=D(2,II)-TEM IF(ABS(TT)-0.0001)445,445,420 420 IF(TT)455,445,425 C ***** MOVE LONGER VECTORS TOWARD END OF TABLE ***** 425 IF(200-NUM)430,430,435 430 NUM=199 435 IJ=NUM DO 440 J=II,NUM D(1,IJ+1)=D(1,IJ) D(2,IJ+1)=D(2,IJ) 440 IJ=IJ-1 GO TO 465 C ***** CHECK FOR DUPLICATE VECTORS IF DISTANCES ARE EQUAL ***** 445 CALL SUB24(D(1,II),KC) CALL ATOM(KC,Z) DO 450 J=1,3 IF(DABS(X(J)+Y(J)-Z(J))-.0001D0)450,450,455 450 CONTINUE GO TO 470 455 CONTINUE IF(200-NUM)470,470,460 C ***** STORE THE RESULT IN VECTOR TABLE ***** 460 II=NUM+1 465 NUM=NUM+1 D(1,II)=100000.*FLOAT(I)+FLOAT((1110-L*100-M*10-NN)*100+K) D(2,II)=TEM 470 CONTINUE 475 CONTINUE C ***** PRINT OUT DISTANCES ***** 16 FORMAT('DISTANCES FROM ATOM ',I5,7X,' TO ATOMS ',I8,' THROUGH ',I3 1) WRITE(NO,16)ITOM,ITAR1,ITAR2 IF(NUM)605,605,480 480 IN(2)=ITOM IN(3)=1 IN(1)=1 485 DO 500 I=1,NUM CALL SUB24(D(1,I),IN(4)) F=FUND(IN(2)) NG=0 IF(IPM)495,490,495 18 FORMAT(1H 5X,2(A6,1X),2(3H (I3,1H,I5,1H)3F7.4,3X),8H D =F6.3) 490 CALL ATOM(IN(4),Z) L=IN(4) WRITE(NO,18)CHEM(ITOM),CHEM(L),IN(2),IN(3),(Y(J),J=1,3 1),IN(4),IN(5),(Z(J),J=1,3),F GO TO 500 495 CALL SUB13 500 CONTINUE C ***** CALCULATE ANGLES ABOUT REF ATOM IF CODE IS 201 ***** 505 IF(IH-2)605,510,510 20 FORMAT('ANGLES AROUND ATOM ',I6) 510 WRITE(NO,20)ITOM L=NUM-1 IF(L)605,605,515 515 IN(4)=ITOM IN(5)=1 IN(1)=2 K=1 DO 540 I=1,L CALL SUB24(D(1,I),IN(2)) LL=IN(2) M=I+1 DO 540 J=M,NUM CALL SUB24(D(1,J),IN(6)) NN=IN(6) F=FUNA(IN(2)) IF(IOUT-1)520,535,535 520 IF(IPM)530,525,530 22 FORMAT(1H 5X,3(A6,1X), 7X,3(2H (I3,1H,I5,1H)),32X,F6.2) 525 WRITE(NO,22)CHEM(LL),CHEM(ITOM),CHEM(NN),(IN(J1),J1=2, 17),F GO TO 540 530 CALL SUB13 GO TO 540 535 PM(K)=F 540 K=K+1 IF(IOUT-1)605,545,545 545 M=NUM-2 LL=2*L IF(M)605,605,550 550 DO 600 I=1,M J1=I+1 NN=((I-1)*(LL-I))/2-1 DO 600 J=J1,L IJK=NN+J X(1)=PM(IJK) IF(X(1)-AMINMX(1))600,555,555 555 IF(AMINMX(2)-X(1))600,560,560 560 K1=J+1 DO 595 K=K1,NUM IJK=NN+K X(2)=PM(IJK) IF(X(2)-AMINMX(1))595,565,565 565 IF(AMINMX(2)-X(2))595,570,570 570 IJK=((J-1)*(LL-J))/2-1+K X(3)=PM(IJK) IF(X(3)-AMINMX(1))595,575,575 575 IF(AMINMX(2)-X(3))595,580,580 580 X(4)=X(1)+X(2)+X(3) IF(X(4)-AMINMX(3))595,585,585 585 IF(AMINMX(4)-X(4))595,590,590 590 WRITE(NO,24) D(1,I),X(1),D(2,I),D(1,J),D(2,J),ITOM,X(2 1),D(2,K),X(3),X(4),D(1,K) 24 FORMAT( 50X,1H(F9.0,1H)/42X,D5.1,7X,1H*/50X,F6.3,2H A/ 152X,1H*/24X,1H(F9.0,2H)*F6.3,6H A *((I3,3H ))5X,D5.1,4H DEG/ 252X,1H*/50X,F6.3,2H A/42X,D5.1,7X,1H*/ 316H SUM OF ANGLES =D6.1,28X,1H(F9.0,1H)// 424X,43H*******************************************) 595 CONTINUE 600 CONTINUE 605 CONTINUE 610 RETURN END SUBROUTINE SUB24(D,I) C IDENTIFICATION CODE UNPACK ROUTINE DIMENSION I(2) I(1)=D/100000. I(2)=D-FLOAT(I(1))*100000. IF((I(2)/100)-555)205,200,205 200 I(2)=I(2)-30000 - 25500 205 RETURN END SUBROUTINE SUMV(X,Y,Z) C COMPUTE THE SUM OF TWO VECTORS C Z(3)=X(3)+Y(3) DIMENSIONX(3),Y(3),Z(3) REAL*8 X,Y,Z DO200I=1,3 200 Z(I)=X(I)+Y(I) RETURN END FUNCTION TRACE(X) C COMPUTE TRACE OF MATRIX X DIMENSIONX(3,3) REAL*8 TRACE,X TRACE=0.D0 DO200I=1,3 200 TRACE=TRACE+X(I,I) RETURN END SUBROUTINE VM(X,Y,Z) C TRANSPOSED VECTOR TIMES MATRIX C Z(3)=X(3)*Y(3,3) DIMENSIONX(3),Y(3,3),Z(3) REAL*8 X,Y,Z DO200J=1,3 Z(J)=0.D0 DO200I=1,3 200 Z(J)=Z(J)+X(I)*Y(I,J) RETURN END FUNCTION VMV(X1,Q,X2) C TRANSPOSED VECTOR * MATRIX * VECTOR C VMV=X1(3)*Q(3,3)*X2(3) TO EVALUATE QUADRATIC OR BILINEAR FORM DIMENSION X1(3),Q(3,3),X2(3) REAL*8 Q,T1,VMV,X1,X2 T1=0.D0 DO 200 J=1,3 200 T1=T1+X1(J)*(X2(1)*Q(J,1)+X2(2)*Q(J,2)+X2(3)*Q(J,3)) VMV=T1 RETURN END FUNCTION VV(X,Y) C TRANSPOSED VECTOR * VECTOR C VV=X(3)*Y(3) DIMENSIONX(3),Y(3) REAL*8 VV,X,Y VV=0.D0 DO200I=1,3 200 VV=VV+X(I)*Y(I) RETURN END SUBROUTINE SETA(A) C DUMMY SUBROUTINE. USED TO SET CONSTRAINTS BETWEEN LATTICE PARAM. RETURN END SUBROUTINE SETP(P) C DUMMY SUBROUTINE. USED TO SET CONSTRAINTS BETWEEN PARAMETERS. RETURN END SUBROUTINE MAKEBAB (INFIL1,NUMMOLE) ******************************************************************************* * FORTRAN CODE FOR COLLATING DATA FOR OPTION 23 OF ORFEE: * * TWO SETS OF DATA ARE TO BE READ IN: * * 1) FIRST THE ORIGINAL ORFEE INPUT * * FROM WHICH IS EXTRACTED THE ORIGINAL SITE 01 * * FRACTIONAL COORDINATES AND UNIT CELL INFORMATION * * 2) SECOND, THE OUTPUT OF OPTION 23 IS READ * * AND MERGED WITH THE SITE 01 INFORMATION * * THE OUTPUT FROM BOTH SECTIONS IS A DATA SET FOR BABEL * * DAVID CLOSE JULY 1995 * ******************************************************************************* CHARACTER*64 INFIL1,INFIL2,OUTFIL CHARACTER*72 LINE(20) CHARACTER*6 CHAR,ATOMNAME CHARACTER*1 ANS,ASYMBOL(250) DIMENSION TITLE(20) OPEN (5,FILE=INFIL1) OPEN (6,FILE='ORFEE20.TMP',STATUS= 'UNKNOWN') READ (5,250) (TITLE(I),I=1,20) 250 FORMAT (20A4) READ (5,260) NATOMS 260 FORMAT (12X,I3) WRITE (*,270) NATOMS 270 FORMAT (//,2X,'NUMBER OF ATOMS READ IN IS ',I3) NTOTAL = 0 I = 0 400 WRITE (*,410) 410 FORMAT (//,2X,'MUST NOW DETERMINE HOW MANY ATOMS THERE ARE',/, 15X,'ARE THE FIRST 4 ATOMS USED TO SPECIFY AXES? (Y/N)' ) READ (*,420) ANS 420 FORMAT (A) IF (ANS .EQ. 'Y' .OR. ANS .EQ. 'y') THEN GO TO 450 ENDIF IF (ANS .EQ. 'N' .OR. AND .EQ. 'n') THEN GO TO 480 ENDIF GO TO 400 450 NATOMS = NATOMS-4 WRITE (*,460) NATOMS 460 FORMAT (2X,'NATOMS = ',I5) READ (5,470) ATOMNAME,X,Y,Z 470 FORMAT (////,A6,20X,3F9.5) GO TO 500 480 READ (5,490) ATOMNAME,X,Y,Z 490 FORMAT (A6,20X,3F9.5) 500 CONTINUE I = I+1 C**** NOW MUST DECIDE ON A SYMBOL TO ASSIGN TO ATOMNAME**** WRITE (*,515) ATOMNAME 515 FORMAT (5X,A6) WRITE (*,520) 520 FORMAT (/,5X, 'NEED TO WRITE ONLY THE ATOMIC SYMBOL NEXT' ) READ (*,420) ASYMBOL(I) WRITE (*,530) 530 FORMAT (/) WRITE (*,550) ASYMBOL(I), X,Y,Z WRITE (6,550) ASYMBOL(I), X,Y,Z 550 FORMAT (1X,A1,2X,3(F9.5,2X)) WRITE (*,560) 560 FORMAT (/) IF (I .EQ. NATOMS) GO TO 570 GO TO 480 C**** HAVE NOW FINISHED PART I. NOW MUST READ DATA FROM FUNCTION 23*** 570 READ (5,580) A,B,C, COSA, COSB, COSC 580 FORMAT (6F9.5) CLOSE (UNIT = 5) WRITE (*,600) 600 FORMAT (2X,'FINISHED PART I. NOW GO GET DATA FROM FUNCTION 23') OPEN (5,FILE='ORFEE2.OUT') IF (NUMMOLE .GT. 1) GO TO 650 620 WRITE (*, 630) 630 FORMAT (//,2X, 'JUST ONE SYMMETRY RELATED SITE REQUESTED') GO TO 700 650 WRITE (*,660) 660 FORMAT (//,2X,'MORE THAN ONE SYMMETRY RELATED SITE REQUESTED') GO TO 800 700 I = 1 WRITE (*,710) NATOMS 710 FORMAT (/,2X, 'ONE SITE WITH',I3,1X,'ATOMS WILL BE READ IN') NTOTAL = 2*NATOMS 750 READ (5,760) CHAR 760 FORMAT (A5) IF (CHAR .EQ. '*****') GO TO 770 GO TO 750 770 READ (5,780) NUMATOM,X,Y,Z 780 FORMAT (2X,I5,25X,3F10.5) WRITE (*,790) ASYMBOL(I),X,Y,Z WRITE (6,790) ASYMBOL(I),X,Y,Z 790 FORMAT (1X,A1,2X,3(F9.5,2X)) IF (I .EQ. NATOMS) GO TO 1000 I = I + 1 GO TO 770 800 I = 1 J = 1 820 FORMAT(I3) NATOMB = NATOMS*NUMMOLE NTOTAL = NATOMS + NATOMB WRITE (*,830) NATOMB 830 FORMAT (//,2X,'TOTAL NUMBER OF ATOMS TO BE READ IN IS',I3,/) 850 READ (5,860) CHAR 860 FORMAT (A5) IF (CHAR .EQ. '*****') GO TO 870 GO TO 850 870 READ (5,780) NUMATOM,X,Y,Z WRITE (*,790) ASYMBOL(I),X,Y,Z WRITE (6,790) ASYMBOL(I),X,Y,Z IF (I .EQ. NATOMS) GO TO 880 I = I + 1 J = J + 1 GO TO 870 880 I = 1 IF (J .EQ. NATOMB) GO TO 1000 J = J + 1 GO TO 870 C***** ORFEE USES COS(ALPHA) ETC. AND BABEL NEEDS ACTUAL ANGLES**** 1000 CONTINUE PI = 3.141592653590 AA=DACOS(COSA) AA = (180.D0/PI)*AA BB=DACOS(COSB) BB = (180.D0/PI)*BB CC=DACOS(COSC) CC = (180.D0/PI)*CC CLOSE (UNIT = 5) CLOSE (UNIT = 6) OPEN (5,FILE='ORFEE20.TMP') WRITE (*,1040) 1040 FORMAT (//,2X,'INPUT FINAL OUTPUT FILENAME FOR BABEL (......BAB)') READ (*,220) OUTFIL 220 FORMAT (A) OPEN (6,FILE=OUTFIL, STATUS = 'UNKNOWN') WRITE (*,1050) NTOTAL 1050 FORMAT (//,2X,'THE CURRENT VALUE OF NATOMS IS ',I5) WRITE (*,1060) 1060 FORMAT (/,5X,'IS THIS CORRECT? (Y/N)') READ (*,420) ANS IF (ANS .EQ. 'Y' .OR. ANS .EQ. 'y') GO TO 1085 IF (ANS .EQ. 'N' .OR. ANS .EQ. 'n') GO TO 1070 1070 WRITE (*,1080) 1080 FORMAT (2X,'NTOTAL = ') READ (*,1090) NTOTAL 1090 FORMAT (I5) 1085 WRITE (*,1095) 1095 FORMAT (//,2X,'NOW JUST GO TO BABEL AND TYPE: ', 1//,5X,'BABEL -if FILENAME.BAB -ot FILENAME.ALC ') 1100 WRITE (6,1110) NTOTAL,AA,BB,CC,A,B,C 1110 FORMAT (I5,6(F9.5,1X)) WRITE (6,1120) (TITLE(I), I=1,20) 1120 FORMAT (20A4) 1150 READ (5,1120,END=2000) (LINE(I), I=1,20) WRITE (6,1120) (LINE(I), I=1,20) GO TO 1150 2000 END