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