Subject: [Summary] Methods for converting internal to external coordinates

From: [email-masked] (Thomas Wieland)
Date: Fri, 5 May 95 13:40:50 +0200

Dear Netters,

last week I posted a question about methods for converting internal to external 
coordinates. I want to thank all peole who replied and the OSC for providing 
thisg reat and really useful list. Below is a summary of the messages I got. 
Some replies are in German; if translations are needed, please refer to me 
or the original senders.

So long,

Thomas Wieland               +---------------+
Dipl. Math.                  |+----    +----+|
Lehrstuhl II f. Mathematik   |\    \   |    ||
Universitaet Bayreuth        | \    \  |    ||
                             |  \    \ |    ||
95440 Bayreuth               |   \    \\    ||
Germany                      |    \    \\   ||
Tel. +49 (921) 553386        |     \     \\ ||
Fax  +49 (921) 553385        |      \-------||
                             +---------------+

Original question:
--------------------
A colleague of mine is looking for conversion methods from internal to
external 3D-coordinates. Given a set of points, their distances,
angles and dihedral angles, he wants to compute Cartesian
coordinates in space. We know that this routine is implemented in many
packages and to avoid reinventing the wheel we wondered if somebody
out there could point us to a description of the method (no code or
complete program since it shall be adapted to our own data structure).

What we already know is MOPAC on the one hand, which, however, has a 
slightly different approach, and BABEL on the other hand, which 
also tells nothing about the underlying methods.

From [email-masked]  
-----------------------------------------------
ich habe die Konvertierung von Z-Matrix (innere Koordinaten) zu kartesischen 
Koordinaten so geloest, das ich das Molekuel nach und
nach aufbaue. 
1.Atom (0,0,0);
2.Atom (d,0,0);
3.Atom (x2+d*cos(180-w), d*sin(180-w), 0);
 alle naechsten Atome lege ich mit dem Abstand zum vorhergenden Atom
 parallel zur X-Achse (x-Koordinate die selbe wie Nachbar). Da jetzt der 
 Bindungswinkel und der Torisonswinkel wilkuerlich sind, aendere ich
 dann den Winkel und den Torionswinkel so, das der in der Z-Matrix definierte
 Wert erhalten wird.

 der Code ist in C in meinem Programm implementiert und benoetigt Routinen 
 zum Berechnen und Aendern von Bindungslaengen, Winkeln und Torsionswinkeln. 

 Wenn Du willst, stelle ich Dir die Routinen zusammen. Es ist allerdings
 keine mathematisch analytische Loesung des Problemes im Sinne von nur
 einer Formel  P(x,y,z) = f(d,w,t)
 der Algorithmus laeuft aber seit laengerem stabil in meinem Programm.

folgendes Beispiel:  Methan  

           H3
         /
   H1--C2
        \ \
        H4  H5

alle distancen 1  alle winkel 109x  Torsion 4-2-1-3  120x
                                            5-2-1-3  240x

H1    0  0  0
C2    1  0  0
H3    1+cos(71)  sin(71)  0
jetzt:

H4 zunaechst y wie Nachbar(C2) = 0
             x = x(2)+d        = 2
             z wie Nachbar(C2) = 0

daraus ergeben sich ein Winkel(4-2-1) von 180x
                                           der wird nun korrigiert zu 109x
der Torsionwinkel ist nun zufllig und mua zu 120x korrigiert werden
usw.

From [email-masked]  
---------------------------------------
Try QCPE.  It's a nice place to "shop".  There are a series of programs for 
calculation of atomic coordinates.  The most likely candidate is COORD 
QCPE#419, but there may also be a newer version.

     QCPE
     Creative Arts Building 181
     Indiana University
     Bloomington, IN 47405
     USA

     Phone: (812) 855-4784
     FAX:   (812) 855-5539

     E-mail: [email-masked]

     For information and catalogues:

     login:     [email-masked]
     password:  guest
                                                                       
 Karl F. Moschner, Ph. D.                                               
 Unilever Research U. S.      

From [email-masked]  
-------------------------------------

        Despite your reluctance to deal with existing code, the easiest way to 
get something working is to adapt an existing program to read/write the data 
files/structures that you prefer... which is what I did some time ago.  Also,
the methodology is quite transparent in a well-documented program.

        Appended is my current version of COORD, adapted from something that
QCPE distributes.  I have, infact, modified it on two separate occasions
because of changes in the form of the output file that were desired.  Feel free
to do whatever you need to do to make it work in your environment.

                           FREDERIC A. VAN-CATLEDGE
The DuPont Company
 ------------------------------ 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

From: David Heisterberg <[email-masked]>
------------------

The following is some Fortran code to do the conversion, but it explains
the method pretty well.
--------
C
C Z2C
C Convert Z matrix to cartesian coordinates.
C
C The Z matrix coordinates are simply local polar coordinates in which
C   c0 - defines the local origin
C   c1 - defines the local z axis
C   c2 - defines the local +zx quadrant
C
      SUBROUTINE z2c (n, nc, zc, cc, q, lerr)
      IMPLICIT CHARACTER (A-Z)
C
C n      - number of centers
C nc     - neighbor (or connectivity) indices
C zc     - zmatrix (r, a, d) coordinates
C cc     - cartesian coordinates
C q      - quality of local coordinates (best if < 1)
C lerr   - lun on which to report fatal errors
C          fatal errors are
C            forward referenced neighbors
C            numerically degenerate local axes
C
      INTEGER n
      INTEGER nc (3, *)
      DOUBLEPRECISION zc (3, *)
      DOUBLEPRECISION cc (3, *)
      DOUBLEPRECISION q (*)
      INTEGER lerr
C
      DOUBLEPRECISION qeps
      PARAMETER (qeps = 0.001)
C
      DOUBLEPRECISION c0 (3), c2 (3)
      DOUBLEPRECISION x (3), y (3), z (3)
      DOUBLEPRECISION xnorm, ynorm, znorm
      DOUBLEPRECISION sint, cost
      DOUBLEPRECISION xfac, yfac, zfac
      DOUBLEPRECISION r, a, d
      INTEGER i, i0, i1, i2
C
C center 1 at the origin
C
      IF (n .LT. 1) RETURN
      cc (1, 1) = 0.0D0
      cc (2, 1) = 0.0D0
      cc (3, 1) = 0.0D0
      q (1) = 0.0D0
C
C center 2 on z axis
C
      IF (n .LT. 2) RETURN
      r = zc (1, 2)
      cc (1, 2) = 0.0D0
      cc (2, 2) = 0.0D0
      cc (3, 2) = r
      q (2) = 0.0D0
C
C center 3 in zx plane
C
      IF (n .LT. 3) RETURN
      r = zc (1, 3)
      a = zc (2, 3)
      i0 = nc (1, 3)
      i1 = nc (2, 3)
C
C check i0 and i1
C
      IF ((i0 .LE. 0) .OR. (i0 .GE. 3)) THEN
        WRITE (lerr, 1000) i0, i
        STOP
      ELSEIF ((i1 .LE. 0) .OR. (i1 .GE. 3)) THEN
        WRITE (lerr, 1001) i1, i
        STOP
      ENDIF
C
C find the local origin and +z axis
C
      c0 (3) = cc (3, i0)
      z (3) = SIGN (1.0D0, cc (3, i1) - c0 (3))
C
C target coordinates are in local +xz quadrant
C
      cc (1, 3) = r * sin (a) * z (3)
      cc (2, 3) = 0.0D0
      cc (3, 3) = r * cos (a) * z (3) + c0 (3)
      q (3) = 0.0D0
C
C centers 4 to n
C
      DO 10000 i = 4, n
        r = zc (1, i)
        a = zc (2, i)
        d = zc (3, i)
        i0 = nc (1, i)
        i1 = nc (2, i)
        i2 = nc (3, i)
C
C check i0, i1, i2
C
        IF ((i0 .LE. 0) .OR. (i0 .GE. i)) THEN
          WRITE (lerr, 1000) i0, i
          STOP
        ELSEIF ((i1 .LE. 0) .OR. (i1 .GE. i)) THEN
          WRITE (lerr, 1001) i1, i
          STOP
        ELSEIF ((i2 .LE. 0) .OR. (i2 .GE. i)) THEN
          WRITE (lerr, 1002) i2, i
          STOP
        ENDIF
C
C c0 is the local origin
C
        c0 (1) = cc (1, i0)
        c0 (2) = cc (2, i0)
        c0 (3) = cc (3, i0)
C
C the local z axis
C check for degeneracy
C
        z (1) = cc (1, i1) - c0 (1)
        z (2) = cc (2, i1) - c0 (2)
        z (3) = cc (3, i1) - c0 (3)
        znorm = SQRT (z (1) * z (1) + z (2) * z (2) + z (3) * z (3))
        IF (znorm .LE. 0.0D0) THEN
          WRITE (lerr, 1011) 'Z', i
          STOP
          ENDIF
C
C c2 defines the local +zx quadrant
C
        c2 (1) = cc (1, i2) - c0 (1)
        c2 (2) = cc (2, i2) - c0 (2)
        c2 (3) = cc (3, i2) - c0 (3)
C
C the local y axis is z cross c2
C check for degeneracy
C
        y (1) = z (2) * c2 (3) - z (3) * c2 (2)
        y (2) = z (3) * c2 (1) - z (1) * c2 (3)
        y (3) = z (1) * c2 (2) - z (2) * c2 (1)
        ynorm = SQRT (y (1) * y (1) + y (2) * y (2) + y (3) * y (3))
        IF (ynorm .LE. 0.0D0) THEN
          WRITE (lerr, 1011) 'Y', i
          STOP
          ENDIF
C
C the local x axis is y cross z
C check for degeneracy
C
        x (1) = y (2) * z (3) - y (3) * z (2)
        x (2) = y (3) * z (1) - y (1) * z (3)
        x (3) = y (1) * z (2) - y (2) * z (1)
        xnorm = SQRT (x (1) * x (1) + x (2) * x (2) + x (3) * x (3))
        IF (xnorm .LE. 0.0D0) THEN
          WRITE (lerr, 1011) 'X', i
          STOP
          ENDIF
C
C cost = ABS (z dot c2)
C sint = ynorm
C if sint is very small compared to cost, "floor" it at qeps * cost
C
        cost = ABS (z (1) * c2 (1) + z (2) * c2 (2) + z (3) * c2 (3))
        sint = ynorm
        IF (sint .LT. qeps * cost) sint = qeps * cost
        q (i) = cost / sint
C
C normalize the axes
C (do this as late as possible to preserve accuracy, and do
C  it as a separate step because the magnitude of the norms
C  are only known to be commensurate with the axis components)
C
        x (1) = x (1) / xnorm
        x (2) = x (2) / xnorm
        x (3) = x (3) / xnorm
        y (1) = y (1) / ynorm
        y (2) = y (2) / ynorm
        y (3) = y (3) / ynorm
        z (1) = z (1) / znorm
        z (2) = z (2) / znorm
        z (3) = z (3) / znorm
C
C fold r, a, d into the 'fac'tors
C
        xfac = r * SIN (a) * COS (d)
        yfac = r * SIN (a) * SIN (d)
        zfac = r * COS (a)
C
C now the target coordinates are a simple affine transformation of
C the local axes
C
        cc (1, i) = xfac * x (1) + yfac * y (1) + zfac * z (1) + c0 (1)
        cc (2, i) = xfac * x (2) + yfac * y (2) + zfac * z (2) + c0 (2)
        cc (3, i) = xfac * x (3) + yfac * y (3) + zfac * z (3) + c0 (3)
10000   CONTINUE
C
      RETURN
C
 1000 FORMAT ('z2c: invalid c0 (', I6, ') for center ', I6)
 1001 FORMAT ('z2c: invalid c1 (', I6, ') for center ', I6)
 1002 FORMAT ('z2c: invalid c2 (', I6, ') for center ', I6)
 1011 FORMAT ('z2c: degenerate local axis (', A, ') for center ', I6)
C
      END

From: tj ODonnell <[email-masked]>
----------------------------------------

This is the way I do it.  With n atoms, there are 3n Cartesian
coordinates needed to descibe the molecule in space.  However,
there are only 3n-6 internal coordinates available.  So we
make the assumption that atom 1 (which has NO internal coordinates)
is at (0,0,0).  Atom 2 has 1 internal coord (a bond distance) and
is assumed to be along the x axis.  Atom 3 has 2 internals (a bond
distance and angle) and is assumed to be in the xy plane.  From this
point on, it is simply a matter of geometry to find the location
of atoms 4...n (having distance, bond angle and torsion (or plane) angle).
Any atom d can be thought of as related to atoms c-b-a such that
d-c-b-a forms two planes dcb and cba.  The torsion angle is the
angle between these planes.  The bond angle is the dcb angle and the
bond distance is the dc distance.
IT IS IMPORTANT TO KNOW THE CONVENTION FOR THE SIGNS OF THE TORSION
ANGLES OR YOU MAY GET THE WRONG ENANTIOMER OF THE MOLECULE.  Unfortunately,
people make different assumptions and if it is not known, there is
nothing you can do but guess.  The only two assumptions are that
a positive torsion angle is either a clockwise or counter-clockwise
rotation of plane dcb relative to plane cba.

Hope this helps.  Let me know if I can help any more.

TJ O'Donnell
http://www.eecs.uic.edu/~tj/

From: [email-masked] (Hinsen Konrad)
--------------------------------------------------
Wenn die internen Koordinaten in Form einer Z-Matrix vorliegen, ist
die Sache einfach: jede Zeile enthaelt die komplette Information
fuer die Position eines Atoms relativ zu den vorherigen, deren
Positionen bereits bekannt sind. Die drei Koordinaten jeder
Zeile sind Kugelkoordinaten in einem Koordinatensystem, das
durch die Positionen der drei Referenzatome definiert ist.
Nur die drei ersten Atome erfordern eine Sonderbehandlung, da
einige Koordinaten nicht festgelegt sind und damit beliebig
gewaehlt werden koennen. Nicht vergessen: Wenn der Winkel
180 Grad betraegt, sollte der Diederwinkel ignoriert werden.
Es hat keinen Sinn, lineare Winkel zu verbieten, wie es etwa
Gaussian macht.

Wenn die Eingabe etwas flexibler sein soll (was ich persoenlich
bevorzuge), wird die Strategie etwas komplizierter. In jedem
Schritt muss man zunaechst ein Atom suchen, fuer das drei
interne Koordinaten relativ zu bereits fixierten Atomen
vorhanden sind. Danach geht es weiter wie oben.

Bei Interesse kann ich einen relativ lesbaren C++-Code anbieten,
der natuerlich eine etwas praezisere Beschreibung darstellt.

Konrad Hinsen                     | E-Mail: [email-masked]
Universite de Montreal            | Fax:  +1-514-343-7586

From: [email-masked] (lee wilson)
---------------------------------------------
Dear Thomas;
        Not sure what you mean about MOPAC6 being different. If you have a
set of internal coordinates: bond distances, bond angles and dihedral
angles- it will give readily the xyz coordinates. That is given a Z-matrix
of internal coordinates  in a MOPAC input file (.MNI or .DAT) with XYZ and
0SCF (either one or both, think) would give a MOPAC output file which
displays the cartesion coordinates. If you want to get the coordinates that
will do it. However I gather that you are interested in incorporating a
routine into software. Seems that you could adapt an appropriate MOPAC routine.

        lee wilson