Subject: SUMMARY: Cartesian to internal conversion

From: Jeffrey J Gosper <[email-masked]>
Date: Wed, 30 Aug 1995 09:21:07 BST

Thanks to those who responsed to my question regarding algorithms 
and/or code for
cartesian to internal coordinate conversion. They were all very 
useful.

There were a number of references to the use of BABEL and MOPAC so 
I've only included
one of each of these.

Here are the responses I received.

********************************************************************
*****************************
the relevant references are :

(1) E. Bright Wilson, Jr., J. C. Decius, Paul C. Cross, Molecular
Vibrations, McGraw-Hill
     Book Company, New York 1955
(2) R. L. Hilderbrandt, J. Chem. Phys. 51(4) (1969) 1654-1659
(3) P. Pulay, Mol. Phys. 17(2) (1969) 197-204
(4) Peter Pulay, Geza Fogarasi, Frank Pang, James E. Boggs, J. Am. 
Chem.
Soc.
    101(10) (1979) 2550-2560
(5) Geza Fogarasi, Xuefeng Zhou, Patterson W. Taylor, Peter Pulay, 
J. Am.
Chem. Soc
     114(21) (1992) 8191-8201
(6) P. Pulay, G. Fogarasi, J. Chem. Phys. 96(4) (1992) 2856-2860
(7) Thomas H. Fischer, Jan Almloef, J. Phys. Chem. 96(24) (1992) 
9768-9774

Esp. refs (3), (4), and (6) give a description of the corresponding
algorithm. It is iterative because of the nonlinearity of this
transformation and it goes as follows :

     x(i+1) = x(i) + A ( q - q(i) )           (G1)

x(i) and x(i+1) are the old and new cartesian coordinates ( so, you 
need
reasonable start cartesians x(0) ! ), q are the given internal 
coordinates
and q(i) are the internal coordinates which correspond to x(i). A is 
the
generalized right inverse of Wilsons B-matrix ( ref. 1 ) :

    A = m BT ( B m BT )**(-1)           (G2)

B is Wilsons B-matrix, BT the transpose of B, and m is any ( 
non-uniqueness
! )
3N x 3N matrix such, that B m BT is not singular. Usually it is a 
unit
matrix. The dimensions of the matrices involved are :

   m   :    3N x 3N
   B   :     M x 3N
   BT :     3N x M
   A   :     3N x M

with M less or equal 3N-6 for nonlinear molecules. N is the number 
of atoms
in the molecule. So, the individual steps are :

     (1)   given some x(i)
     (2)   compute q(i) from the given x(i) ( a trivial task ), 
compute B (
actually this is just
            a by-product of the computing of q(i) ! )
     (3)   make B m BT and invert it
     (4)   compute A and the new cartesians x(i+1) from equations 
(G1) and
(G2) above.

Repeat steps (1) to (4) until q-q(i) is acceptable.

The ab initio packages  TEXAS ( Peter Pulay ), TURBOMOLE ( Reinhart 
Ahlrichs
), and
DISCO ( Jan Almloef ) certainly have a module to do precisely this
conversion from internal to cartesian coordinates. See esp. section 
IV of
ref. (5). Because I am working with TURBOMOLE and this is not free 
for
distribution you should probably contact one of the above mentioned 
authors
to get the corresponding modul or even better, write your own code.

Good luck
Heinz

Heinz Schiffer
Hoechst AG
Scientific Computing
65926 Frankfurt am Main
Phone   ++49-69-305-2330
Fax       ++49-69-305-81162
Email     [email-masked]
********************************************************************
*****************************
You can get the MOPAC source code at www.ccl.net and look
at the GMETRY.FOR code which has this convertion routines.

Best regards,

Edgardo Garcia
Cristol Chem & Biochem
University of Colorado
BOULDER CO 80309
USA
********************************************************************
*****************************
        Reply to:   RE>CCL:Cartesian to Internal Conversion

such algorithms are often part of standard MD programs (etc)
which have the potential in internal coords but the hard
work (integration, etc) in cartesian.

for example try program VENUS, coauthored by Bill Hase and
others (including myself).  Bill Hase is at Wayne State Univ
                   [email-masked]
VENUS has scheduled to be submitted to Quantum Chem
Program Exchange this month.

Kieran Lim

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
x                                                                x
x Kieran F Lim             ph:  + [61] (3) 9344 7090  (office)   x
x School of Chemistry           + [61] (3) 9344 6481  (School)   x
x University of Melbourne       + [61] (3) 9435 0259  (home)     x
x Parkville      _.-_|\                                          x
x VIC  3052     /      \   fax: + [61] (3) 9347 5180 (preferred) x
x AUSTRALIA     \_.--.o/        + [61] (3) 9344 6233             x
x                     v                                          x
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
x                                                                x
x           [email-masked]           x
x                                                                x
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
x                                                                x
x Please note:  All Melbourne phone and fax numbers have changed x
x               Numbers  xxx-xxxx have become  9xxx-xxxx.        x
x          eg:  + [61] (3) 347 5180 becomes + [61] (3) 9347 5180 x
x                     (03) 435 0259 becomes       (03) 9435 0259 x
x                                                                x
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
********************************************************************
*****************************
Hello; Here is how to get Cart's --> Internal with Mopac.
  aigout=ab in. geom. output; noxyz =don't output any Cartesians
  Mopac will also go from internal--> Cartesians

---------------
INPUT:

am1 0scf aigout noxyz
(CN)2 CARBENE, MP2 opt. Input: SPARTAN 6-31G* geom

C      -0.25275 1  -0.33068 1  -0.68783 1
C      -0.60299 1  -0.53886 1   0.71870 1
C      -0.70363 1   0.55964 1   1.55581 1
C       1.19821 1  -0.46268 1  -0.94150 1
O      -1.05071 1  -0.16199 1  -1.54391 1
N      -0.80356 1   1.41686 1   2.31356 1
N       2.31853 1  -0.54893 1  -1.09341 1
------------------

OUTPUT:

********************************************************************
***********
 **                            MOPAC 93 (c) Fujitsu                 
          **
 
********************************************************************
***********

                                 AM1 CALCULATION RESULTS

********************************************************************
***********
 *                   MOPAC  93.01               CALC'D. Wed Aug 23 
10:58:47 1995
 *  AIGOUT   - IN ARC FILE, INCLUDE AB-INITIO GEOMETRY
 *   T=      - A TIME OF  3600.000 SECONDS REQUESTED
 *  DUMP=N   - RESTART FILE WRITTEN EVERY  3600.000 SECONDS
 *  AM1      - THE AM1 HAMILTONIAN TO BE USED
 *  NOXYZ    - CARTESIAN COORDINATES NOT TO BE PRINTED
 *  0SCF     - AFTER READING AND PRINTING DATA, STOP
 
********************************************************************
***060BY060
 AM1 0SCF AIGOUT NOXYZ
 (CN)2 CARBENE, MP2 opt. Input: SPARTAN 6-31G* geom

    ATOM   CHEMICAL  BOND LENGTH    BOND ANGLE     TWIST ANGLE
   NUMBER  SYMBOL    (ANGSTROMS)     (DEGREES)      (DEGREES)
    (I)                  NA:I          NB:NA:I      NC:NB:NA:I     
NA   NB   NC

      1      C
      2      C         1.46435  *                                  1
      3      C         1.38477  *     119.02875  *                 2 
   1
      4      C         1.47887  *     112.75119  *   94.13493  *   1 
   2    3
      5      O         1.18240  *     123.66246  *  -90.51055  *   1 
   2    3
      6      N         1.14848  *     175.76608  * -177.86070  *   3 
   2    1
      7      N         1.13386  *     177.68500  *  -30.11819  *   4 
   1    2
  C: (AM1): M.J.S. DEWAR ET AL, J. AM. CHEM. SOC. 107 3902-3909 
(1985)
  N: (AM1): M.J.S. DEWAR ET AL, J. AM. CHEM. SOC. 107 3902-3909 
(1985)
  O: (AM1): M.J.S. DEWAR ET AL, J. AM. CHEM. SOC. 107 3902-3909 
(1985)

      MOLECULAR POINT GROUP   :   C1

      RHF CALCULATION, NO. OF DOUBLY OCCUPIED LEVELS = 16

            INTERATOMIC DISTANCES
0
                  C  1       C  2       C  3       C  4       O  5  
     N  6
 
--------------------------------------------------------------------
----------
     C    1    .000000
     C    2   1.464354    .000000
     C    3   2.455582   1.384769    .000000
     C    4   1.478870   2.450794   3.301317    .000000
     O    5   1.182400   2.337068   3.201481   2.347542    .000000
     N    6   3.516480   2.531531   1.148476   4.258541   4.175395  
  .000000
     N    7   2.612204   3.437895   4.169022   1.133857   3.421177  
 5.021874
0
                  N  7
 ------------------
     N    7    .000000
 GEOMETRY IN MOPAC Z-MATRIX FORMAT
 AM1 0SCF AIGOUT NOXYZ
 (CN)2 CARBENE, MP2 opt. Input: SPARTAN 6-31G* geom

  C     .00000000  0     .0000000  0     .0000000  0    0    0    0 
     4.0000
  C    1.46435433  1     .0000000  0     .0000000  0    1    0    0 
     4.0000
  C    1.38476850  1  119.0287496  1     .0000000  0    2    1    0 
     4.0000
  C    1.47887031  1  112.7511903  1   94.1349333  1    1    2    3 
     4.0000
  O    1.18239987  1  123.6624556  1  -90.5105473  1    1    2    3 
     6.0000
  N    1.14847603  1  175.7660839  1 -177.8607024  1    3    2    1 
     5.0000
  N    1.13385740  1  177.6850027  1  -30.1181886  1    4    1    2 
     5.0000

  GEOMETRY IN GAUSSIAN Z-MATRIX FORMAT
 AM1 0SCF AIGOUT NOXYZ
 (CN)2 CARBENE, MP2 opt. Input: SPARTAN 6-31G* geom

  C
  C     1   r21
  C     2   r32            1   a321
  C     1   r41            2   a412           3   d4123          0
  O     1   r51            2   a512           3   d5123          0
  N     3   r63            2   a632           1   d6321          0
  N     4   r74            1   a741           2   d7412          0

   r21            1.464354
   r32            1.384769
   r41            1.478870
   r51            1.182400
   r63            1.148476
   r74            1.133857
   a321         119.028750
   a412         112.751190
   a512         123.662456
   a632         175.766084
   a741         177.685003
   d4123         94.134933
   d5123        269.489453
   d6321        182.139298
   d7412        329.881811

 TOTAL CPU TIME:              .03 SECONDS

 == MOPAC DONE ==
 JOB FINISHED
********************************************************************
*************
The Gaussian utility newzmat can do the job.

In general, I am rarely satisfied with automatically generated 
Z-matrices,
whether from Gaussian newzmat, MOPAC or other software.  At best 
they give
a fair starting point.  At worst you end up doing the Z-matrix by 
hand
anyway, especially if you want to impose specific symmetric 
constraints or
a define particular reaction coordinate.  Alternative methods such
Spartan's Cartessian coordinate optimization which supports 
"internal"
constraints and Gaussian 94 use of redundant internal coordinates 
are
overall more convenient and, in my opinion, superior.

Regards, Karl
 
____________________________________________________________________
___
/                                                                   
    \
| Comments are those of the author and not Unilever Research U. S.  
    |
|                                                                   
    |
| Karl F. Moschner, Ph. D.                                          
    |
|                                                                   
    |
| Unilever Research U. S.      e-mail: 
[email-masked] |
| 45 River Road                Phone:  (201) 943-7100 x2629         
    |
| Edgewater, NJ 07020          FAX:    (201) 943-5653               
    |
\___________________________________________________________________
____/
********************************************************************
Pat Walters' group at Tucson offers Babel, which goes "both ways,"
and runs on dos and unix.
The address is [email-masked]
I believe it ftp's from there.

John Reissner         Pembroke State University     Pembroke NC  
28372  USA
[email-masked]     vox: (910)521-6425    fax: 
(910)521-6649
********************************************************************
**********************
Dear Jeff,

on April 27 I posted the same question the other way round, i.e.
internal to cartesian. I can send you my summary again if you wish.

For your problem, you should first take a look at BABEL which is
a nice tool for conversion between several chemical file formats.
I think you can get the BABEL package from ftp, maybe from the Ohio
chemistry software server

In general I think the this conversion isn't quite hard. Bonds are
computed as Euclidean distances. Angles by the equation

                         (x0-x1) o (x2-x1)
      cos phi  =        -------------------
                         |x0-x1| |x2-x1|

if x1 is the top of the angle and "o" the standard scalar product.
Torsion angles are defined as the angle between the planes made up
by (x0,x1,x2) on one hand and (x1,x2,x3) on the other. So compute
these planes and their normal vectors and then calculate the
angle between them by the equation above. Ask me if you need more 
help.

By the way, I think I also read a solution to this problem in the 
following
book : Tim Clark: A handbook of computational chemistry, Wiley, New 
York
1985, 4. ed.

Best wishes,

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

P.S.: Take a look at MOLGEN:
      http://btm2xd.mat.uni-bayreuth.de/molgen
********************************************************************
***
In reponce to your request I have written a TCL code for this 
purpose.
I am glad to sent it to you - but it is still a little developmental 
- in
that I have not checked the patch I put on the four quadrant 
correector yet.

A more elegant solution might be a developmental program called 
chemsh.
This is being authored by Paul Sherwood (et all) at Daresbury.  I 
cannot say
if he would be willing to release a version to you, but if you were 
interested
in testing it - it does have just about any z-matrix conversion you 
could
want.

Try [email-masked]

I hope this helps

All the best

Alex

+--------------------------------------------------+
|Alternate E-mail [email-masked]            |
|www home @ http://www.bath.ac.uk/~chpajt/home.html|
+--------------------------------------------------+
********************************************************************
**********************
provided, the cartesian co-ordinates of atoms are given in vectors
XAT,YAT,ZAT then you can use the following routines (written by me 
using
algorithms described by many people, so I do not remember).

Good luck,

Frank Eisenmenger

+-------------------------------------------------------------------
-----+
| Institute of Biochemistry                                         
     |
| Medical Faculty of the Humboldt-University (Charite)              
     |
|                                                                   
     |
| Hessische Str. 3-4                                                
     |
| 10115 Berlin, Germany                                             
     |
|                                                                   
     |
| Phone: +49-30-28468-612 or -239       FAX: +49-30-28468-600       
     |
| E-mail: [email-masked]                           
     |
+-------------------------------------------------------------------
-----+

c *******************************
      real function bndlen(i1,i2)

c ..................................................
c  PURPOSE: return bond length between atoms i1 & i2
c
c  INPUT:   i1,i2 - indices of 2 atoms
c
c ..................................................

      COMMON /COOR/ XAT(1000),YAT(1000),ZAT(1000)

      bndlen=sqrt( (xat(i1)-xat(i2))**2
     #            +(yat(i1)-yat(i2))**2+(zat(i1)-zat(i2))**2 )

      return
      end
c **********************************
      real function valang(i1,i2,i3)

c .........................................
c  PURPOSE: return valence angle (i1,i2,i3)
c           [in rad.] with 'i2' as vertex
c
c  INPUT:   i1,i2,i3 - indices of 3 atoms
c .........................................

      COMMON /COOR/ XAT(1000),YAT(1000),ZAT(1000)

      h1=xat(i2)
      h2=yat(i2)
      h3=zat(i2)
      x1=xat(i1)-h1
      x2=yat(i1)-h2
      x3=zat(i1)-h3
      y1=xat(i3)-h1
      y2=yat(i3)-h2
      y3=zat(i3)-h3

      x=x1*x1+x2*x2+x3*x3
      y=y1*y1+y2*y2+y3*y3
      u=x*y

      if (u.ne.0.0) then

        a=(x1*y1+x2*y2+x3*y3)/sqrt(u)
        a=max(a,-1.0)
        a=min(a,1.0)
        valang=acos(a)
        return

      else
        write (*,'(a,3i5)')' valang> Error in coordinates of atoms 
#: '
     #                     ,i1,i2,i3
        stop
      endif

      end
c *************************************
      real function dihedr(i1,i2,i3,i4)

c .............................................
c  PURPOSE: return dihedral angle (i1,i2,i3,i4)
c           [in rad.]
c
c  INPUT:   i1,i2,i3,i4 - indices of four atoms
c
c  CALLS:   none
c .............................................

      COMMON /COOR/ XAT(1000),YAT(1000),ZAT(1000)

      x1=xat(i2)-xat(i1)
      y1=yat(i2)-yat(i1)
      z1=zat(i2)-zat(i1)
      x2=xat(i3)-xat(i2)
      y2=yat(i3)-yat(i2)
      z2=zat(i3)-zat(i2)
      ux1=y1*z2-z1*y2
      uy1=z1*x2-x1*z2
      uz1=x1*y2-y1*x2
      x1=xat(i4)-xat(i3)
      y1=yat(i4)-yat(i3)
      z1=zat(i4)-zat(i3)
      ux2=z1*y2-y1*z2
      uy2=x1*z2-z1*x2
      uz2=y1*x2-x1*y2

      u1=ux1*ux1+uy1*uy1+uz1*uz1
      u2=ux2*ux2+uy2*uy2+uz2*uz2
      u=u1*u2

      if (u.ne.0.0) then
        a=(ux1*ux2+uy1*uy2+uz1*uz2)/sqrt(u)
        a=max(a,-1.0)
        a=min(a,1.0)
        dihedr=acos(a)
        if (ux1*(uy2*z2-uz2*y2)+uy1*(uz2*x2-ux2*z2)+
     #      uz1*(ux2*y2-uy2*x2).lt.0.0) dihedr =-dihedr
        return
      else
        write (*,'(a,4i5)')' dihedr> Error in coordinates of atoms 
#: '
     #                     ,i1,i2,i3,i4
        stop
      endif

      end

********************************************************************
Dear Jeff,

Many years ago when I was young and writing my programs myself I 
wrote such
a programm which converts cartesian coordinates into internal ones. 
It worked
in that years so I hope it will work now...

If you will have problems with it - ask me.

Cheers, Vlad

c
c
      subroutine fmintc(cartes,cint,int,natoms)
c
c Generate internal coordinates from cartesian
c INPUT:
c cartes(3,*) - initial cartesian coordinates
c int(3,*) - connectivity of a given atom with others (bond length 
-int(1,*),
c            valence angle -int(2,i), torsion- int(3,*)
c natoms - number of atoms
c
c OUTPUT:
c cint(3,*) - Z-matrix
c
      implicit real*8 (a-h,o-z)
      dimension cint(3,*),int(3,*),cartes(3,*),ext(3,4),v1(3),v2(3)
      data rad_to_degree/57.2958/
c
      do i=1,3
        do j=1,3
          cint(j,i) = 0.0
        enddo
      enddo
c
      cint(1,2) = sqrt( (cartes(1,1)-cartes(1,2))**2 +
     &                  (cartes(2,1)-cartes(2,2))**2 +
     &                  (cartes(3,1)-cartes(3,2))**2 )
      cint(1,3) = sqrt( (cartes(1,3)-cartes(1,2))**2 +
     &                  (cartes(2,3)-cartes(2,2))**2 +
     &                  (cartes(3,3)-cartes(3,2))**2 )
      do 2 i=1,3
        v1(i) = cartes(i,1) - cartes(i,2)
 2      v2(i) = cartes(i,3) - cartes(i,2)
      cint(2,3) = acos( (v1(1)*v2(1) + v1(2)*v2(2) + v1(3)*v2(3))/
     &            (cint(1,2)*cint(1,3)) ) * rad_to_degree
c
      do 1 i=4,natoms
        do 3 j=1,3
          ext(j,4) = cartes(j,i)
          kk = 4 - j
          do 3 k=1,3
 3          ext(k,j) = cartes(k,int(kk,i))
 1      call intrnl(ext,cint(1,i) )
c
      return
      end
c
c
      subroutine intrnl(ext,cint )
c
      implicit real*8 (a-h,o-z)
      dimension ext(3,*),cint(*),h1(3),h2(3),vnorm(3),unitar(3,3),
     &          vmidl(3,2)
      data pi/3.14159265453/,pp/57.2958/
c
      do i=1,3
        h1(i)=ext(i,1)-ext(i,2)
        h2(i)=ext(i,3)-ext(i,2)
      enddo
c
      call vecmul(h2,h1,vnorm)
      call rulcos(vnorm,unitar(1,3) )
      call rulcos(h2,unitar(1,1) )
      call vecmul(unitar(1,3),unitar(1,1),unitar(1,2) )
c
      do 2 i=1,2
        do 2 j=1,3
 2        vmidl(j,i) = ext(1,i+2) * unitar(1,j) +
     &                 ext(2,i+2) * unitar(2,j) +
     &                 ext(3,i+2) * unitar(3,j)
c
      cint(1) = sqrt( (ext(1,4)-ext(1,3))**2 +
     &                (ext(2,4)-ext(2,3))**2 +
     &                (ext(3,4)-ext(3,3))**2 )
      cint(2) = acos( (vmidl(1,1)-vmidl(1,2))/cint(1) ) * pp
      subd=sqrt( (vmidl(2,2)-vmidl(2,1))**2 +
     &           (vmidl(3,2)-vmidl(3,1))**2 )
      arg = (vmidl(2,2)-vmidl(2,1))/subd
      if ( abs(arg).gt. 1.D00 ) arg = sign( 0.9999999D00, arg )
CC    cint(3) = -acos( arg ) * pp
      cint(3) =  acos( arg ) * pp
      if ( vmidl(3,2)-vmidl(3,1).lt.0.0 ) cint(3)= -cint(3)
c
      return
      end
c
c
       subroutine rulcos( ventry,rcos )
c
      implicit real*8 (a-h,o-z)
      dimension ventry(*),rcos(*)
c
       s=0.0
       do 1 i=1,3
 1       s=s+ventry(i)**2
       s=sqrt(s)
       do 2 i=1,3
 2       rcos(i)=ventry(i)/s
c
       return
       end
c
c
      subroutine vecmul(v1,v2,vnorm)
c
      implicit real*8 (a-h,o-z)
      dimension v1(*),v2(*),vnorm(*)
c
      vnorm(1) = v1(2)*v2(3) - v1(3)*v2(2)
      vnorm(2) = v1(3)*v2(1) - v1(1)*v2(3)
      vnorm(3) = v1(1)*v2(2) - v1(2)*v2(1)
c
      return
      end
Received: from n32.rhea.cnusc.fr by monge.brunel.ac.uk with SMTP 
(PP)
          id <[email-masked]>; Wed, 23 Aug 1995 12:54:21 
+0100
Received: by n32.rhea.cnusc.fr (AIX 3.2/UCB 5.64/4.03) id AA30959;
          Wed, 23 Aug 1995 13:51:50 +0200
Date: Wed, 23 Aug 1995 13:51:50 +0200
From: [email-masked] (\(SP_6\))
Message-Id: <[email-masked]>
********************************************************************
***********

/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\
 Dr. Jeff Gosper                                         
 Dept. of Chemistry		                        
 BRUNEL University                                     
 Uxbridge Middx UB8 3PH, UK                            
 voice:  01895 274000 x2187                            
 facsim: 01895 256844                                  
 internet/email/work:   [email-masked]     
 internet/WWW: http://http2.brunel.ac.uk:8080/~castjjg 
\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/