From chm6@midway.uchicago.edu Sun May 23 08:52:44 1993 Date: Sun, 23 May 93 13:52:44 CDT From: "charles h martin" Message-Id: <9305231852.AA06784@midway.uchicago.edu> To: CHEMISTRY@ccl.net Subject: zmat->cartesian summary Dear Folks, As there were so many request for the zmat -> cartesian transformation program, here is a summary of the responses. ---------------------- (here is a summary of a summary sent to me) Thanks to all who responded to my request concerning programs to convert cartesian coordinates to z-matrix form. I was asked to summarize responses --- from Mark A. Zottola: Gaussian 92 has a utility Newzmat that converts PDB files to z-matrix form. from Leslie Glasser (009LGZS@WITSVMA.WITS.AC.ZA): Prof. Glasser has written a MathCAD file that converts crystal to cartesian coordinates or internal coordinates. You need MathCAD (either DOS or Windows) to use the script. from Gregory L. Durst: QCPE has a program EXTOIN (QCMP070) for PC's that converts MM2 external cartesian coordinates into a z-matrix. It is in Fortran. from James Stewart: MOPAC has a subroutine XYZINT that can form the core of a program to perform this conversion. Feel free to use it. in addition, Zoran Zdravkovski reponds: MOPAC does have an option to get cartesian coordinates from Z- matrix, by using the keywords XYZ 0SCF. AMPAC probably has the same option, and you could probaly easily find the source of earlier versions of either somewhere on the network. Try using, for example, archie -c MOPAC (or AMPAC) to locate an FTP site with the programs. ------------------------- In addition, I recieved a few source code listing which I now provide: From: FREDVC@esds01.es.dupont.com Try the program appended below. It's a reworking of a QCPE code; I have used it for years. ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ FREDERIC A. VAN-CATLEDGE Scientific Computing Division || Office: (302) 695-1187 Central Research & Development Dept. || FAX: (302) 695-9658 The Du Pont company || P. O. Box 80320 || Internet: fredvc@esvax.dnet.dupont.com Wiilmington DE 19880-0320 || ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ---------------------------------- CUT HERE ---------------------------------- CFPP$ SWITCH RENUMB=100:10,FORMAT=9000:10,INDAL=2 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 ******* 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 --------------- From: GERARD@XRAY.BMC.UU.SE (Gerard Kleijwegt a.k.a. gerard@xray.bmc.uu.se) c c =========================================================================== c c Subroutines to convert internal to Cartesian coordinates c c Gerard J. Kleywegt c Department of Molecular Biology c Biomedical Centre c University of Uppsala c Uppsala, SWEDEN c c "gerard@xray.bmc.uu.se" c c =========================================================================== c c CALLS: c ====== c C2CART (NATOMS,XYZ(3,*),CORINT(3,*),IREFAT(3,*)) c C2CAR4 (NATOMS,XYZ(3,*),CORINT(3,*),IREFAT(3,*)) c c VARIABLES: c ========== c NATOMS = integer = nr of atoms c XYZ(3,*) = real = returned Cartesian coordinates c CORINT(3,*) = real = internal coordinates (distance, angle, dihedral) c IREFAT(3,*) = integer = indices of the reference atoms c c USE: c ==== c Normally, one would merely call C2CART. However, this puts the c first atom at (0,0,0). If you want to orient the atoms in a c particular way, you may set the XYZ coordinates of the first c three atoms yourself and then call C2CAR4 instead ! c c =========================================================================== c subroutine c2car4 (nat,xyz,int,ref) c c ... C2CAR4 - convert internal to Cartesian coordinates c c ... RETAINS THE COORDINATES OF THE FIRST THREE ATOMS !!! c implicit none c real twopi,pi,degtor,rtodeg c parameter (twopi=6.2831853071796) parameter (pi=0.50*twopi) parameter (degtor=twopi/360.0) parameter (rtodeg=360.0/twopi) c real xyz(3,*),int(3,*) real xa,xb,ya,yb,za,zb,xyb,xpa,zpa,xpb,zpb,costh,sinth real rbc,cosph,sinph,xqa,zqa,a,d,h,coskh,sinkh,sina,cosa real sind,cosd,xd,yd,zd,xpd,ypd,zpd,xqd,yqd,zqd,xrd,zrd real ypa,yza c integer ref(3,*),nat,i,n1,n2,n3,k c code ... c if (nat .lt. 4) return c do i=4,nat c n3 = ref(1,i) n2 = ref(2,i) n1 = ref(3,i) c d = int(1,i) a = int(2,i) c c ... the following line may need a sign-change, depending c on how you defined your dihedral angles c h = -int(3,i) c xa = xyz(1,n1) - xyz(1,n3) ya = xyz(2,n1) - xyz(2,n3) za = xyz(3,n1) - xyz(3,n3) c xb = xyz(1,n2) - xyz(1,n3) yb = xyz(2,n2) - xyz(2,n3) zb = xyz(3,n2) - xyz(3,n3) c xyb = sqrt(xb**2 + yb**2) k = 1 c if (xyb .lt. 0.1) then k = 0 c xpa = za zpa = -xa xa = xpa za = zpa c xpb = zb zpb = -xb xb = xpb zb = zpb c xyb = sqrt(xb**2 + yb**2) end if c costh = xb/xyb sinth = yb/xyb xpa = xa*costh + ya*sinth ypa = ya*costh - xa*sinth c rbc = sqrt (xb**2 + yb**2 + zb**2) sinph = zb/rbc cosph = sqrt (1.0 - sinph**2) xqa = xpa*cosph + za*sinph zqa = za*cosph - xpa*sinph c yza = sqrt (ypa**2 + zqa**2) coskh = ypa/yza sinkh = zqa/yza c sina = sin(degtor*a) cosa = cos(degtor*a) sind = sin(degtor*h) cosd = cos(degtor*h) c xd = d*cosa yd = d*sina*cosd zd = d*sina*sind c ypd = yd*coskh - zd*sinkh zpd = zd*coskh + yd*sinkh xpd = xd*cosph - zpd*sinph c zqd = zpd*cosph + xd*sinph xqd = xpd*costh - ypd*sinth yqd = ypd*costh + xpd*sinth c if (k .ne. 1) then xrd = -zqd zrd = xqd xqd = xrd zqd = zrd end if c xyz (1,i) = xyz(1,n3) + xqd xyz (2,i) = xyz(2,n3) + yqd xyz (3,i) = xyz(3,n3) + zqd c cc print *,xyz(1,i),xyz(2,i),xyz(3,i) c end do c return end c c c subroutine c2cart (nat,xyz,int,ref) c c ... C2CART - convert internal to Cartesian coordinates c implicit none c real twopi,pi,degtor,rtodeg c parameter (twopi=6.2831853071796) parameter (pi=0.50*twopi) parameter (degtor=twopi/360.0) parameter (rtodeg=360.0/twopi) c real xyz(3,*),int(3,*) c integer ref(3,*),nat c code ... c if (nat .lt. 1) return c xyz (1,1) = 0.0 xyz (2,1) = 0.0 xyz (3,1) = 0.0 if (nat .lt. 2) return c xyz (1,2) = int(1,2) xyz (2,2) = 0.0 xyz (3,2) = 0.0 if (nat .lt. 3) return c xyz (1,3) = xyz(1,2) - int(1,3)*cos(degtor*int(2,3)) xyz (2,3) = int(1,3)*sin(degtor*int(2,3)) xyz (3,3) = 0.0 if (nat .lt. 4) return c call c2car4 (nat,xyz,int,ref) c return end ----------------- I hope this helps ================================================================== Charles H. Martin email: chm6@quads.uchicago.edu U.S. nail: c/o Freed Group The James Franck Institute and The Department of Chemistry The University of Chicago 5640 South Ellis Avenue Chicago, Illinois 60637 Work: (312) 702-3457 Fax: (312) 702-6853 ==================================================================