|
|
From Jeffrey.Gosper@brunel.ac.uk Wed Aug 30 04:21:52 1995
Received: from monge.brunel.ac.uk for Jeffrey.Gosper@brunel.ac.uk
by www.ccl.net (8.6.10/950822.1) id EAA13065; Wed, 30 Aug 1995 04:21:48 -0400
Received: from chem-pc-18.brunel.ac.uk by monge.brunel.ac.uk with SMTP (PP)
id <07531-0@monge.brunel.ac.uk>; Wed, 30 Aug 1995 09:21:15 +0100
Date: Wed, 30 Aug 1995 09:21:07 BST
From: Jeffrey J Gosper
Reply-To: Jeffrey.Gosper@brunel.ac.uk
Subject: SUMMARY: Cartesian to internal conversion
To: chemistry@www.ccl.net
Message-ID:
Priority: Normal
MIME-Version: 1.0
Content-Type: TEXT/PLAIN; CHARSET=US-ASCII
Status: R
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 schiffer@msmwia.hoechst.hoechst-ag.d400.de
********************************************************************
*****************************
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
hase@sun.chem.wayne.edu
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 kieran_lim.chemistry@muwayf.unimelb.edu.au 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:
Karl.F.Moschner@urlus.sprint.com |
| 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 babel@mercury.aichem.arizona.edu
I believe it ftp's from there.
John Reissner Pembroke State University Pembroke NC
28372 USA
reissner@pembvax1.pembroke.edu 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 psh@dl.ac.uk
I hope this helps
All the best
Alex
+--------------------------------------------------+
|Alternate E-mail A.J.Turner@Bath.ac.uk |
|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: eisenmen@orion.rz.mdc-berlin.de
|
+-------------------------------------------------------------------
-----+
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 <20604-0@monge.brunel.ac.uk>; 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: lcaogto@n32.rhea.cnusc.fr (\(SP_6\))
Message-Id: <9508231151.AA30959@n32.rhea.cnusc.fr>
********************************************************************
***********
/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\
Dr. Jeff Gosper
Dept. of Chemistry
BRUNEL University
Uxbridge Middx UB8 3PH, UK
voice: 01895 274000 x2187
facsim: 01895 256844
internet/email/work: Jeffrey.Gosper@brunel.ac.uk
internet/WWW: http://http2.brunel.ac.uk:8080/~castjjg
\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/
|
