|
|
*******************************************************************************
* 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
|
