|
|
C >From fredvc@esvax.dnet.dupont.com
C -------------------------------------
C
C Despite your reluctance to deal with existing code, the easiest way to
C get something working is to adapt an existing program to read/write the data
C files/structures that you prefer... which is what I did some time ago. Also,
C the methodology is quite transparent in a well-documented program.
C
C Appended is my current version of COORD, adapted from something that
C QCPE distributes. I have, infact, modified it on two separate occasions
C because of changes in the form of the output file that were desired. Feel free
C to do whatever you need to do to make it work in your environment.
C
C FREDERIC A. VAN-CATLEDGE
C The DuPont Company
C ------------------------------ CUT HERE -------------------------------------
PROGRAM COORD_XYZ
C**************************************************************************
C
C
C COORD CALCULATES THE COORDINATES IN 3-DIMENSIONAL MOLECULES,GIVEN
C THE BOND LENGTHS,BOND ANGLES AND DIHEDRAL ANGLES
C
C THE CARD OUTPUT FROM COORD CAN BE USED DIRECTLY AS INPUT TO MINDO3
C OR REINER SUSTMANN'S NDDO PROGRAMS
C
C
C ******* INPUT DATA FOR EACH MOLECULE *****
C FIRST CARD HAS MOLECULE CHARGE IN COLS.1-2,AND HAS THE MOLECULE
C NAME IN COLS.3-63
C SECOND CARD. NOAT I2
C IZAT(1) I2
C IZAT(2) I2
C IZAT(3) I2
C KWIK I1
C R12 F7.4
C R23 F7.4
C THETA F14.7
C EACH SUCCESSIVE CARD. NA I2
C NB I2
C NC I2
C ND I2
C IZAT(ND) I2
C ILAZY I1
C RCD F7.4
C THBCD F14.7
C PHABCD F14.7
C
C A CARD WITH 99 IN COLS.1-2 MUST BE ADDED AT THE END OF
C THE ENTIRE DECK OF MOLECULES TO TERMINATE THE PROGRAM
C
C CARDS FOR ATOMS WITH IZAT>99 ARE NOT PUNCHED
C
C**************************************************************************
C
IMPLICIT REAL*8 (A-H,O-Z)
DIMENSION IAT(100), NAT(100,10),COVRAD(99),ATWT(99)
DIMENSION X(100),Y(100),Z(100),RIJ(100,100),IZAT(100),NAME(15)
LOGICAL QDONE(100),QFILE
CHARACTER*2 EL(99)
CHARACTER*4 SCREXT,CEXT,OUTEXT,INEXT
CHARACTER*20 OFIL
CHARACTER*24 NEWFIL
C
COMMON/COVRAD/COVRAD,ATWT,EL
C
COMMON /WORK/ X, Y, Z, DEGRAD
C
DATA DEGRAD/1.74532925199432D-02/
DATA SCREXT/'.SCR'/,CEXT/'.XYZ'/,OUTEXT/'.OUT'/,INEXT/'.Q'/
C
SQRT2 = DSQRT ( 2.0D0 )
SQRT3 = DSQRT ( 3.0D0 )
THTET = DACOS(-1.D0/3.D0)/DEGRAD
THTRI = 120.
C
D OPEN(UNIT=46,FILE='COORD.TRK',TYPE='NEW')
C
100 CONTINUE
WRITE ( 6, 9000 )
READ ( 5, 9010, END=250 ) OFIL
NCH = INDEX(OFIL, ' ') -1
C
NEWFIL = OFIL(1:NCH) // SCREXT
WRITE ( 6, 9020 ) NEWFIL
OPEN(UNIT=23,NAME=NEWFIL,TYPE='NEW',CARRIAGECONTROL='LIST')
C
NEWFIL = OFIL(1:NCH) // CEXT
WRITE ( 6, 9030 ) NEWFIL
OPEN(UNIT=21,NAME=NEWFIL,TYPE='NEW',CARRIAGECONTROL='LIST')
C
NEWFIL = OFIL(1:NCH) // OUTEXT
OPEN(UNIT=16,NAME=NEWFIL,TYPE='NEW',CARRIAGECONTROL='FORTRAN')
WRITE ( 6, 9040 ) NEWFIL
C
NEWFIL = OFIL(1:NCH) // INEXT
INQUIRE(FILE=NEWFIL,EXIST=QFILE)
IF (QFILE) THEN
WRITE ( 6, 9045 ) NEWFIL
IGET = 15
IPRMPT = 36
OPEN(UNIT=IGET,TYPE='OLD',NAME=NEWFIL,READONLY)
OPEN(UNIT=IPRMPT,NAME='NLA0:',TYPE='NEW')
ELSE
IGET = 5
IPRMPT = 6
END IF
C
C************************************************************************
C
WRITE ( IPRMPT, 9050 )
READ ( IGET, 9060 ) ICHG, ( NAME(I), I = 1, 15 )
IF (ICHG .EQ. 99) GO TO 240
C
OPEN(UNIT=22,TYPE='SCRATCH',FORM='UNFORMATTED')
C
WRITE ( 23, 9070 ) ICHG, ( NAME(I), I = 1, 15 )
C
C>>>>>NOAT IS THE NUMBER OF ATOMS. IZAT(I) IS THE ATOMIC NUMBER OF ATOM
C NUMBER I. KWIK ALLOWS AUTOMATIC CALCULATION OF COORDINATES OF
C ATOMS 1, 2, 3 IN SIMPLE CASES, KWIK = 0, TETRAHEDRAL, KWIK = 1,
C ANGLE 120 DEGREES, KWIK = 2, ANGLE SUPPLIED AS DATAUM.
C R12, R23 ARE BOND LENGTHS. THETA IS THE 123 BOND ANGLE.
C
WRITE ( IPRMPT, 9080 )
READ(IGET,9090)NOAT,(IZAT(I),I=1,3),KWIK,R12,R23,THETA
IF (KWIK .EQ. 0) THEN
THETA = THTET
ELSE IF (KWIK .EQ. 1) THEN
THETA = THTRI
ENDIF
C
C>>>>>WRITE INPUT TO VARIOUS FILES
WRITE(23,9100)NOAT,(IZAT(K),K=1,3),KWIK,R12,R23,THETA
WRITE ( 16, 9240 ) ( NAME(I), I = 1, 15 ) , ICHG
WRITE ( 16, 9250 ) NOAT, ( IZAT(I), I = 1, 3 ) , KWIK
IF (KWIK .LT. 2) THEN
WRITE ( 16, 9260 ) R12, R23, THETA
ELSE
WRITE ( 16, 9265 ) R12, R23, THETA
END IF
WRITE ( 16, 9270 )
C
C>>>>>CALC TRIG FCNS OF BOND ANGLE
THETAK = THETA * DEGRAD
CCOS = DCOS ( THETAK )
SSIN = DSIN ( THETAK )
C
C>>>>>CALC XYZ FOR ATOMS 1, 2, & 3
DO 110 I = 1, 3
X(I) = 0.0
Y(I) = 0.0
Z(I) = 0.0
QDONE(I) = .TRUE.
110 CONTINUE
C
X(2) = R12
X(3) = R12 - R23 * CCOS
Y(3) = R23 * SSIN
C
C>>>>>SET COORDINATE FLAG FOR REMAINING ATOMS
DO 120 I = 4, NOAT
QDONE(I) = .FALSE.
120 CONTINUE
C
C>>>>>LOOP TO CALC REMAINING XYZ COORDINATES
DO 160 I = 4, NOAT
C
C>>>>>>>ATOMS NA, NB, NC, HAVE KNOWN COORDINATES AND ARE NOT COLLINEAR.
C IZAT(ND) IS THE ATOMIC NUMBER OF ATOM ND. THBCD IS THE BCD BOND
C ANGLEIN DEGREES AND PHABCD THE DIHEDRAL ANGLE OF CD RELATIVE TO
C AB, MEASURED CLOCKWISE ALONG THE DIRECTION B TO C. ILAZY ALLOWS
C AUTOMATIC CALCULATION OF ANGLES IN NORMAL TETRAHEDRAL AND PLANAR
C SYSTEMS. ILAZY = 0, 1, 2, 3, 4, 5 TETRAHEDRAL WITH DIHEDRAL
C ANGLES OF RESPECTIVELY 0, 60, 120, 180, 240, 310 DEGREES.
C ILAZY = 6, 7 PLANAR CIS, TRANS RESPECTIVELY. ILAZY = 8, ATOMS
C B, C, D COLLINEAR. ILAZY = 9, ANGLES FROM DATA
C
WRITE ( IPRMPT, 9110 )
READ(IGET,9120)NA,NB,NC,ND,IZAT(ND),ILAZY,RCD,THBCD,PHABCD
IF (ILAZY .LT. 6) THEN
THBCD = THTET
PHABCD = ILAZY * 60.
ELSE IF (ILAZY .LT. 8) THEN
THBCD = THTRI
PHABCD = (ILAZY - 6) * 180.
ELSE IF (ILAZY .EQ. 8) THEN
THBCD = 180.
PHABCD = 0.
ELSE
IF (DABS(THBCD - 180.) .LT. 0.01) THEN
THBCD = 180.
PHABCD = 0.
END IF
ENDIF
C
C>>>>>>>WRITE INPUT TO .SCR FILE
WRITE(23,9130)NA,NB,NC,ND,IZAT(ND),ILAZY,RCD,THBCD,PHABCD
C
C>>>>>>>CHECK THAT COORDINATES OF ATOMS NA, NB, NC HAVE BEEN CALCULATED
130 CONTINUE
IF(.NOT.(QDONE(NA).AND.QDONE(NB).AND.QDONE(NC)))GO TO 230
C
140 CONTINUE
WRITE ( 22 ) NA, NB, NC, ND, IZAT(ND), ILAZY, RCD, THBCD, PHABCD
C
THETA = THBCD
PHI = PHABCD
R = RCD
CALL GETXYZ(NA, NB, NC, ND, R, THETA, PHI)
QDONE(ND) = .TRUE.
160 CONTINUE
C
REWIND 22
DO 180 I = 4, NOAT
170 CONTINUE
READ ( 22 ) NA, NB, NC, ND, IZ, ILAZY, RCD, THBCD, PHABCD
IF (ILAZY .LT. 9) THEN
WRITE(16,9140)NA,NB,NC,ND,IZ,ILAZY,RCD,THBCD,PHABCD
ELSE
WRITE(16,9145)NA,NB,NC,ND,IZ,ILAZY,RCD,THBCD,PHABCD
END IF
180 CONTINUE
CLOSE ( UNIT=22 )
C
WRITE ( 16, 9150 )
C
NCON = NOAT
DO 200 I = 1, NOAT
IF (IZAT(I) .GT. 99) NCON = NCON - 1
DO 190 J = 1, NOAT
RIJ(I,J)=DSQRT((X(I)-X(J))**2+(Y(I)-Y(J))**2+(Z(I)-Z(J))**2)
190 CONTINUE
200 CONTINUE
C
C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
C CODE FOR READING A .C FILE
C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
C READ ( 21, 9000 ) ID
C READ ( 21, 9010 ) NATOMS, ICHG, IMULT, IANG
C BOHR = IANG .EQ. 0
C NEL = -ICHG
C DO 100 K = 1, NATOMS
C READ ( 21, 9020 ) IZ, X(K), Y(K), Z(K)
C 100 CONTINUE
C9000 FORMAT(A80)
C9010 FORMAT(4I)
C9020 FORMAT(I,3F)
C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
WRITE ( 21, 9160 ) ( NAME(I), I = 1, 15 )
WRITE ( 21, 9170 ) NOAT, ICHG
SX = 0.
SY = 0.
SZ = 0.
SW = 0.
DO 220 I = 1, NOAT
IF (IZAT(I) .LE. 99) THEN
IT = IZAT(I)
W = ATWT(IT)
SW = SW + W
SX = SX + W*X(I)
SY = SY + W*Y(I)
SZ = SZ + W*Z(I)
IF (COVRAD(IT) .EQ. 0.) THEN
RAD1 = 1.65
ELSE
RAD1 = 1.15 * COVRAD(IT)
ENDIF
IAT(I) = 0
DO 210 J = 1, NOAT
IF (I .NE. J) THEN
IF (IZAT(J) .LE. 99) THEN
JT = IZAT(J)
IF (COVRAD(JT) .EQ. 0.) THEN
RAD2 = 1.65
ELSE
RAD2 = 1.15 * COVRAD(JT)
ENDIF
IF (RIJ(I, J) .LE. RAD1+RAD2) THEN
IAT(I) = IAT(I) + 1
NAT(I,IAT(I)) = J
ENDIF
ENDIF
ENDIF
210 CONTINUE
ENDIF
220 CONTINUE
C
C.....MOVE TO CENTER-OF-MASS COORDINATES
XO = SX/SW
YO = SY/SW
ZO = SZ/SW
DO 225 I = 1, NOAT
X(I) = X(I) - XO
Y(I) = Y(I) - YO
Z(I) = Z(I) - ZO
IF (IZAT(I) .LE. 99) THEN
WRITE(21,9180)IZAT(I),X(I),Y(I),Z(I),I,IAT(I),
& (NAT(I,K),K=1,IAT(I))
WRITE(16,9190)I,X(I),Y(I),Z(I),(NAT(I,L),L=1,IAT(I))
END IF
225 CONTINUE
C
WRITE ( 16, 9240 ) ( NAME(K), K = 1, 15 ) , ICHG
WRITE ( 16, 9200 )
CALL NATPRT ( RIJ, NOAT, NOAT, 100, 10, 0 )
GO TO 100
C
230 CONTINUE
WRITE ( 6, 9210 ) NA, NB, NC, ND
WRITE ( 23, 9210 ) NA, NB, NB, ND
IF (.NOT.QDONE(NA)) WRITE ( 23, 9220 ) NA
IF (.NOT.QDONE(NB)) WRITE ( 23, 9220 ) NB
IF (.NOT.QDONE(NC)) WRITE ( 23, 9220 ) NC
240 CONTINUE
WRITE ( 23, 9230 ) ICHG
250 CONTINUE
STOP
9000 FORMAT(//1H$,7X,'ENTER FILE ID (CTRL/Z TO END)... ')
9010 FORMAT(A)
9020 FORMAT('0 INPUT IS SAVED IN FILE ',A24)
9030 FORMAT('0 CONNECTIVITY FILE IS ',A24)
9040 FORMAT('0 SUMMARY IS IN FILE ',A24,', WHICH CAN BE PRINTED')
9045 FORMAT('0(INPUT IS TAKEN FROM FILE ',A24,')')
9050 FORMAT(1H0,'CHARGE, NAME = (I2,15A4)'/)
9060 FORMAT (I,15A4)
9070 FORMAT(I2,',',15A4)
9080 FORMAT(1H0,'# ATOMS, Z1, Z2, Z3, KWIK, R12, R23, THETA'/,
1 ' KWIK = 0,1,2 FOR TET, 120 OR THETA SUPPLIED'/)
9090 FORMAT(5I,3F)
9100 FORMAT(I4,4(',',I4),3(',',F12.5))
9110 FORMAT(1H0,'NA, NB, NC, ND, Z(ND), LAZY, RCD,'
1 ' THBCD, PHABCD'/' LAZY=0,1,2,3,4,5 FOR TET WITH',
2 ' PHI=0,60,120,180,240,310: 6,7 FOR CIS,TRANS: '
3 '8 COLINEAR: 9 SUPPLIED'/)
9120 FORMAT(6I,3F)
9130 FORMAT(I4,5(',',I4),3(',',F12.5))
9140 FORMAT(1H ,4I5,2I8,F12.5,2F12.3,3X,'(STD. ANGLES)')
9145 FORMAT(1H ,4I5,2I8,F12.5,2F12.3)
9150 FORMAT(1H0,2X,'NO. OF ATOM',5X,'X-COORDINATE',5X,
1 'Y-COORDINATE',5X,'Z-COORDINATE'/)
9160 FORMAT(15A4,20X)
9170 FORMAT(I3,', ',I2,', 1, 1, <- (ANGSTROMS)')
9180 FORMAT(I3,', ',F15.10,', ',F15.10,', ',F15.10,', ',
1 I3,', ',I1,', ',9(I3,', '))
9190 FORMAT(1H ,8X,I2,3X,3F17.6,10I5)
9200 FORMAT(1H0,21HINTERATOMIC DISTANCES,//)
9210 FORMAT(1H0,38HCOORDS.OF 1 REFERENCE ATOM UNAVAILABLE,4I5)
9220 FORMAT('0 ATOM',I4.' XYZ NOT CALC-D')
9230 FORMAT(I2)
9240 FORMAT(1H1,15A4,7HCHARGE=,I3)
9250 FORMAT (8H0NOAT = I2, 14H IZAT(1) = I2, 14H IZAT(2) = I2,
114H IZAT(3) = I2, 11H KWIK = I1)
9260 FORMAT (' R12 = ',F7.4,4X,'R23 = ',F7.4,4X,'THETA = ',F10.3,3X,
& '(STD. ANGLE)')
9265 FORMAT (' R12 = ',F7.4,4X,'R23 = ',F7.4,4X,'THETA = ',F10.3)
9270 FORMAT(1H0,2X,'NA',3X,'NB',3X,'NC',3X,'ND',3X,
1 'IZAT(ND)',2X,'ILAZY',4X,'RCD',7X,'THBCD',
2 6X,'PHABCD'/)
END
SUBROUTINE NATPRT(A,N,M,MA,NC,II)
IMPLICIT REAL*8 (A-H,O-Z)
C...Translated by FPP 5.0 (3.03N9) 08/16/92 22:09:51
C...Switches: format=9000:10,indal=2
C**************************************************************
C *
C SUBROUTINES AND FUNCTIONS CALLED FROM THIS ROUTINE *
C *
C NO SUBROUTINES OR FUNCTIONS CALLED FROM THIS PROGRAM UNIT *
C *
C**************************************************************
DIMENSION A(MA,MA)
C MATPRT PRINTS MATRICES** FROM KLOPMANS PROGRAM SCF
KK = 0
NCM1 = NC - 1
J = 0
L = 1
IF (II - 1 .GT. 0) THEN
L = II - 10
II = 0
KK = 5
ENDIF
DO 150 IZ = L, M, NC
NIF = IZ + NCM1
NIF = MIN0 ( M,NIF )
J = J + N - II * (IZ - 1 )
IF (J - 52 .GE. 0) THEN
I = 0
J = 0
ELSE
I = 1
ENDIF
IF (I + KK - 1 .LT. 0) GO TO 100
IF (I + KK - 1 .EQ. 0) GO TO 110
GO TO 120
100 CONTINUE
WRITE ( 16, 9000 )
110 CONTINUE
WRITE ( 16, 9010 ) ( K, K = IZ, NIF )
120 CONTINUE
IJ = 2 * (NIF - IZ + 1 ) + 1
IF (II .GT. 0) THEN
DO 130 IR = IZ, N
JJ = IR
JJ = MIN0 ( NIF,JJ )
WRITE ( 16, 9020 ) IR, ( A(IR,IC), IC = IZ, JJ )
130 CONTINUE
ELSE
DO 140 IR = 1, N
WRITE ( 16, 9020 ) IR, ( A(IR,IC), IC = IZ, NIF )
140 CONTINUE
ENDIF
150 CONTINUE
RETURN
9000 FORMAT (1H1)
9010 FORMAT(1H0,/10I12)
9020 FORMAT(1H I2,10F12.5)
END
BLOCK DATA ATOMS
IMPLICIT REAL*8 (A-H,O-Z)
DIMENSION COVRAD(99),ATWT(99)
CHARACTER*2 EL(99)
C
COMMON/COVRAD/COVRAD,ATWT,EL
C
DATA COVRAD/
1 0.3200, 0.9300, 1.2300, 0.9000, 0.8200, 0.7700,
2 0.7500, 0.7300, 0.7200, 0.7100, 1.5400, 1.3600,
3 1.1800, 1.1100, 1.0600, 1.0200, 0.9900, 0.9800,
4 2.0300, 1.7400, 1.4400, 1.3200, 1.2200, 1.1800,
5 1.1700, 1.1700, 1.1600, 1.1500, 1.1700, 1.2500,
6 1.2600, 1.2200, 1.2000, 1.1600, 1.1400, 1.1200,
7 2.1600, 1.9100, 1.6200, 1.4500, 1.3400, 1.3000,
8 1.2700, 1.2500, 1.2500, 1.2800, 1.3400, 1.4800,
9 1.4400, 1.4100, 1.4000, 1.3600, 1.3300, 1.3100,
1 2.3500, 1.9800, 1.6900, 1.6500, 1.6500, 1.6400,
2 1.6300, 1.6200, 1.8500, 1.6100, 1.5900, 1.5900,
3 1.5800, 1.5700, 1.5600, 1.5600, 1.5600, 1.4400,
4 1.3400, 1.3000, 1.2800, 1.2600, 1.2700, 1.3000,
5 1.3400, 1.4900, 1.4800, 1.4700, 1.4600, 1.4600,
6 1.4500, 0.0000, 0.0000, 0.0000, 0.0000, 1.6500,
7 0.0000, 1.4200, 0.0000, 0.0000, 0.0000, 0.0000,
8 0.0000, 0.0000, 0.0000/
C
DATA ATWT/
& 1.00797, 4.00260, 6.93900, 9.01220, 10.81100,
& 12.01115, 14.00670, 15.99940, 18.99840, 20.18300,
& 22.98980, 24.31200, 26.98150, 28.08600, 30.97380,
& 32.06400, 35.45300, 39.94800, 39.10200, 40.08000,
& 44.95600, 47.90000, 50.94200, 51.99600, 54.93800,
& 55.84700, 58.93300, 58.71000, 63.54000, 65.37000,
& 69.72000, 72.59000, 74.92200, 78.96000, 79.90900,
& 83.80000, 85.47000, 87.62000, 88.90500, 91.22000,
& 92.90600, 95.94000, 98.00000, 101.07000, 102.90500,
& 106.40000, 107.87000, 112.40000, 114.82000, 118.69000,
& 121.75000, 127.60000, 126.90400, 131.30000, 132.90500,
& 137.34000, 138.91000, 140.12000, 140.90700, 144.24000,
& 147.00000, 150.35000, 151.96000, 157.25000, 158.92400,
& 162.50000, 164.93000, 167.26000, 168.93400, 173.04000,
& 174.97000, 178.49000, 180.94800, 183.85000, 186.20000,
& 190.20000, 192.20000, 195.09000, 196.96700, 200.59000,
& 204.37000, 207.19000, 208.98000, 210.00000, 210.00000,
& 222.00000, 13*250.0/
C
DATA EL/' H','HE','LI','BE',' B',' C',' N',' O',
1 ' F','NE','NA','MG','AL','SI',' P',' S','CL',
2 'AR',' K','CA','SC','TI',' V','CR','MN','FE',
3 'CO','NI','CU','ZN','GA','GE','AS','SE','BR',
4 'KR','RB','SR',' Y','ZR','NB','MO','TC','RU','RH',
5 'PD','AG','CD','IN','SN','SB','TE',' I','XE',
6 'CS','BA','LA','CE','PR','ND','PM','SM','EU','GD',
7 'TB','DY','HO','ER','TM','YB','LU','HF','TA',' W',
8 'RE','OS','IR','PT','AU','HG','TL','PB','BI','PO',
9 'AT','RN','FR','RA','AC','TH','PA',' U','NP','PU','AM','CM',
1'BK','CF','ES'/
END
SUBROUTINE GETXYZ(NA, NB, NC, ND, RIJ, THETA, PHI)
C**************************************************************************
C IMSL ROUTINES CALLED: DMRRRR - MATRIX MULTIPLY (REAL * REAL)
C DCRGRG - MATRIX COPY (REAL TO REAL)
C DLINRG - MATRIX INVERSION (REAL)
C**************************************************************************
IMPLICIT REAL*8 (A-H,O-Z)
DIMENSION X(100),Y(100),Z(100)
COMMON /WORK/ X, Y, Z, DEGRAD
C
DIMENSION Q(4,4), R(4,4), RTOT(4,4), TMP(4,4), E(4,4)
C
DATA EPS/1.0D-08/
C
C-----LOAD Q ARRAY
DO 5 J = 1, 4
DO 4 I = 1, 4
RTOT(I,J) = 0.
Q(I,J) = 0.
E(I,J) = 0.
4 CONTINUE
RTOT(J,J) = 1.
E(J,J) = 1.
Q(4,J) = 1.
5 CONTINUE
Q(1,1) = X(NA)
Q(2,1) = Y(NA)
Q(3,1) = Z(NA)
Q(1,2) = X(NB)
Q(2,2) = Y(NB)
Q(3,2) = Z(NB)
Q(1,3) = X(NC)
Q(2,3) = Y(NC)
Q(3,3) = Z(NC)
Q(4,4) = 0.
C
D WRITE(46,7999)NA, NB, NC, ND, RIJ, THETA, PHI
D7999 FORMAT('1ENTERING *GETXYZ*; ATOM SEQUENCE, R/TH/PH = ',4I5,
D & 3F12.5)
D ASSIGN 9 TO IQW
D GO TO 1000
9 CONTINUE
C
C-----TRANSLATE ORIGIN TO ATOM B
10 CONTINUE
CALL DCRGRG(4, E, 4, R, 4)
DO 15 I = 1, 3
R(I,4) = -Q(I,2)
15 CONTINUE
C
D WRITE(46,7998)
D7998 FORMAT(////'0TRANSLATE ORIGIN TO ATOM B')
D ASSIGN 16 TO IRQW
D GO TO 2000
C
16 CONTINUE
CALL DMRRRR(4, 4, R, 4, 4, 4, Q, 4, 4, 4, TMP, 4)
CALL DCRGRG(4, TMP, 4, Q, 4)
CALL DMRRRR(4, 4, R, 4, 4, 4, RTOT, 4, 4, 4, TMP, 4)
CALL DCRGRG(4, TMP, 4, RTOT, 4)
C
D ASSIGN 19 TO IQW
D GO TO 1000
C
19 CONTINUE
C
C-----Z-ROTATION: ATOM C INTO XZ-PLANE (YC->0)
20 CONTINUE
CALL DCRGRG(4, E, 4, R, 4)
DXY = DSQRT(Q(1,3)**2 + Q(2,3)**2)
IF (DXY .LT. EPS) GO TO 29
R(1,1) = Q(1,3)/DXY
R(2,2) = R(1,1)
R(1,2) = Q(2,3)/DXY
R(2,1) = -R(1,2)
C
D WRITE(46,7997)
D7997 FORMAT(////'0Z-ROTATION: ATOM C INTO XZ-PLANE (YC->0)')
D ASSIGN 26 TO IRQW
D GO TO 2000
26 CONTINUE
C
CALL DMRRRR(4, 4, R, 4, 4, 4, Q, 4, 4, 4, TMP, 4)
CALL DCRGRG(4, TMP, 4, Q, 4)
CALL DMRRRR(4, 4, R, 4, 4, 4, RTOT, 4, 4, 4, TMP, 4)
CALL DCRGRG(4, TMP, 4, RTOT, 4)
C
D ASSIGN 29 TO IQW
D GO TO 1000
29 CONTINUE
C
C-----Y-ROTATION: ATOM C ONTO Z-AXIS (XC->0)
30 CONTINUE
CALL DCRGRG(4, E, 4, R, 4)
DXZ = DSQRT(Q(3,3)**2 + Q(1,3)**2)
IF (DXZ .LT. EPS) GO TO 39
R(3,3) = Q(3,3)/DXZ
R(1,1) = R(3,3)
R(3,1) = Q(1,3)/DXZ
R(1,3) = -R(3,1)
C
D WRITE(46,7996)
D7996 FORMAT(////'0Y-ROTATION: ATOM C ONTO Z-AXIS (XC->0)')
D ASSIGN 36 TO IRQW
D GO TO 2000
C
36 CONTINUE
CALL DMRRRR(4, 4, R, 4, 4, 4, Q, 4, 4, 4, TMP, 4)
CALL DCRGRG(4, TMP, 4, Q, 4)
CALL DMRRRR(4, 4, R, 4, 4, 4, RTOT, 4, 4, 4, TMP, 4)
CALL DCRGRG(4, TMP, 4, RTOT, 4)
C
D ASSIGN 39 TO IQW
D GO TO 1000
C
39 CONTINUE
C
C-----Z-ROTATION: ATOM A INTO XZ-PLANE (YA = 0)
40 CONTINUE
CALL DCRGRG(4, E, 4, R, 4)
DXY = DSQRT(Q(1,1)**2 + Q(2,1)**2)
IF (DXY .LT. EPS) GO TO 49
R(1,1) = Q(1,1)/DXY
R(2,2) = R(1,1)
R(1,2) = Q(2,1)/DXY
R(2,1) = -R(1,2)
C
D WRITE(46,7995)
D7995 FORMAT(////'0Z-ROTATION: ATOM Z INTO XZ-PLANE (YA = 0)')
D ASSIGN 46 TO IRQW
D GO TO 2000
C
46 CONTINUE
CALL DMRRRR(4, 4, R, 4, 4, 4, Q, 4, 4, 4, TMP, 4)
CALL DCRGRG(4, TMP, 4, Q, 4)
CALL DMRRRR(4, 4, R, 4, 4, 4, RTOT, 4, 4, 4, TMP, 4)
CALL DCRGRG(4, TMP, 4, RTOT, 4)
C
D ASSIGN 49 TO IQW
D GO TO 1000
C
49 CONTINUE
C
C-----ADD ATOM D IN LOCAL COORDINATES
IF (THETA .EQ. 180.) THEN
Q(3,4) = Q(3,3) + RIJ
ELSE
THBCD = THETA * DEGRAD
PHABCD = PHI * DEGRAD
DZ = RIJ * DCOS(THBCD)
DXY = RIJ * DSIN(THBCD)
Q(3,4) = Q(3,3) - DZ
Q(1,4) = DXY * DCOS(PHABCD)
Q(2,4) = DXY * DSIN(PHABCD)
END IF
Q(4,4) = 1.
D WRITE(46,6995)
D6995 FORMAT(////'0ADD ATOM D IN LOCAL COORDINATES')
D ASSIGN 59 TO IQW
D GO TO 1000
59 CONTINUE
C
C-----PUT INVERSE OF RTOT-MATRIX INTO R-MATRIX; APPLY INVERSE TRANSFORM
CALL DLINRG(4, RTOT, 4, R, 4)
WRITE(46,7994)
7994 FORMAT(////'0APPLY INVERSE TRANSFORM')
D ASSIGN 96 TO IRQW
D GO TO 2000
96 CONTINUE
CALL DMRRRR(4, 4, R, 4, 4, 4, Q, 4, 4, 4, TMP, 4)
CALL DCRGRG(4, TMP, 4, Q, 4)
D ASSIGN 99 TO IQW
D GO TO 1000
C
99 CONTINUE
X(ND) = Q(1,4)
Y(ND) = Q(2,4)
Z(ND) = Q(3,4)
C
RETURN
C============================================================================
1000 CONTINUE
D WRITE(46,9999)
D9999 FORMAT('0 CURRENT (Q MATRIX)'/)
D DO 1010 I = 1, 4
D WRITE(46,9998)(Q(I,J),J = 1, 4)
D1010 CONTINUE
D9998 FORMAT('|',4F15.9,'|')
D GO TO IQW
C============================================================================
2000 CONTINUE
D WRITE(46,8999)
D8999 FORMAT('0TRANSFORMATION: (R-MATRIX) * (Q-MATRIX)'/)
D DO 2010 I = 1, 4
D WRITE(46,8998)(R(I,J),J=1,4), (Q(I,K),K=1,4)
D2010 CONTINUE
D8998 FORMAT(' |',4F15.9,'|',6X,'|',4F15.9,'|')
D GO TO IRQW
C============================================================================
C
END
|
