Conversion of Cartesian Coordinates to Sybyl Mol1 file type
- From: Dongchul Lim <lim |-at-| rani.chem.yale.edu>
- Subject: Conversion of Cartesian Coordinates to Sybyl Mol1 file
type
- Date: Tue, 22 Dec 92 17:18:03 EST
Dear netters,
Here's my 2 hour hack of nawk script which converts cartesian coordinates
to sybyl MOL1 file type. If you don't want to bother ftp'ing previously
announced q2m or need a source code, you may want to try this script,
although q2m will probably generate better sybyl files.
This is not tested with old awk, so it'll be safer to use nawk, mawk, gawk,
etc. If your sysmtem is a BSD or its clone, you'll need to change the first
line of the script to "#!/bin/nawk -f". SYS V users don't need to
change it,
as long as "nawk" is located at /usr/bin directory.
I'll explain briefly how it works.
1. read in a line (first line is treated as title and ignored)
2. if the line consists of the data types "integer number number
number",
the first integer is stored in "an" array, and the following three
numbers
are stored in "x", "y", and "z" arrays.
Otherwise, the program exits with error message.
3. continue until EOF is encountered, or filename changes.
4. print out sybyl MOL1 file.
1) construct bond list by comparing the distance of each atom pair and the
sum of Bragg-Slater radii of them (multiplied by "factor").
2) count the number of neighbors for each atom.
3) determine sybyl atom type based on atomic number and # of neighbors.
4) print out atomtype, coordinates, and bond lists
5. if there're more than one file, repeat the above procedure.
As you may notice, determination of sybyl atom type is not perfect and you
will perhaps have to modify the script to obtain appropriate atom types.
To use this program, cut out the script from this mail and name it
"c2m"
(or whatever you like). Then make it executable by "chmod +x c2m".
On command line, type "c2m file1 file2...", where file1, file2... are
cartesian type files (the first line is assumed to be a comment).
This program is provided AS IS, without any warranty, implied or expressed.
-Dongchul Lim
Department of Chemistry, Yale Univ.
------------------- CUT HERE ----------------------------------------------
#!/usr/bin/nawk -f
# |-at-| (#) c2m - Convert cartesian coordinates to Sybyl MOL1 file
# Usage: c2m file1 file2...
# Author: Dongchul Lim
# Department of Chemistry, Yale University, lim |-at-|
rani.chem.yale.edu
# Date: December 22, 1992
# Disclaimer: This program is provided "AS IS", without any warranty.
# The author will not be responsible for any disters which may result from
# using this program.
BEGIN {
OLDFILE = FILENAME;
factor = 1.5; # bond tolerance factor
number = "^[-+]?([0-9]+[.]?[0-9]*|[.][0-9]+)([eE][-+]?[0-9]+)?$";
integer = "^[-+]?[0-9]+$";
natoms = 0; nbonds = 0;
# Bragg-Slater Radii (change values if needed)
r[1] =0.35; r[2] =0.50; r[3] =1.45; r[4] =1.05; r[5] =0.85; r[6] =0.70; r[7]
=0.65; r[8] =0.60;
r[9] =0.55; r[10]=0.65; r[11]=1.80; r[12]=1.50; r[13]=1.25; r[14]=1.10;
r[15]=1.00; r[16]=1.00;
r[17]=1.00; r[18]=0.95; r[19]=2.20; r[20]=1.80; r[21]=1.60; r[22]=1.40;
r[23]=1.35; r[24]=1.40;
r[25]=1.40; r[26]=1.40; r[27]=1.35; r[28]=1.35; r[29]=1.16; r[30]=1.35;
r[31]=1.35; r[32]=1.30;
r[33]=1.25; r[34]=1.15; r[35]=1.15; r[36]=1.15; r[37]=1.10; r[38]=2.35;
r[39]=2.00; r[40]=1.80;
r[41]=1.55; r[42]=1.45; r[43]=1.45; r[44]=1.35; r[45]=1.30; r[46]=1.35;
r[47]=1.40; r[48]=1.60;
r[49]=1.55; r[50]=1.55; r[51]=1.45; r[52]=1.45; r[53]=1.40; r[54]=1.40;
r[55]=1.50; r[56]=1.50;
r[57]=1.50; r[58]=1.50; r[59]=1.50; r[60]=1.50; r[61]=1.50; r[62]=1.50;
r[63]=1.50; r[64]=1.50;
r[65]=1.50; r[66]=1.50; r[67]=1.50; r[68]=1.50; r[69]=1.50; r[70]=1.50;
r[71]=1.50; r[72]=1.50;
r[73]=1.50; r[74]=1.50; r[75]=1.50; r[76]=1.50; r[77]=1.50; r[78]=1.50;
r[79]=1.44; r[80]=1.50;
r[81]=1.50; r[82]=1.50; r[83]=1.50; r[84]=1.50; r[85]=1.50; r[86]=1.50;
r[87]=1.50; r[88]=1.50;
r[89]=1.50; r[90]=1.50; r[91]=1.50; r[92]=1.50; r[93]=1.50; r[94]=1.50;
r[95]=1.50; r[96]=1.50;
}
FNR == 1 {next;}
$1 !~ integer || $2 !~ number || $3 !~ number || $4 !~ number {
printf("\007Cartesian format error in line %d: %s\n", FNR, $0);
exit;
}
# read in atomic number and cartesian coordinates
{
if (FILENAME != OLDFILE) {
printSybyl(); printf("\n");
natoms = 0; nbonds = 0;
}
if ($1 <= 0+0) next; # skip dummy atoms
natoms++;
an[natoms] = $1; x[natoms] = $2; y[natoms] = $3; z[natoms] = $4;
OLDFILE = FILENAME;
}
END {
printSybyl();
}
# print Sybyl MOL1 file type
function printSybyl() {
# construct list of bonds
for(i=1;i<=natoms;i++) {
for(j=i+1;j<=natoms;j++) {
if (dist(i, j) <= r[an[i]]*factor + r[an[j]]*factor) {
nbonds++;
bond_a[nbonds] = i;
bond_b[nbonds] = j;
}
}
}
# construct list of neighbors
for(i=1;i<=natoms;i++) {
neighbors[i] = 0;
for(j=1;j<=nbonds;j++) if (bond_a[j] == i || bond_b[j] == i)
neighbors[i]++;
}
# print out sybyl MOL1 type file
printf("%4d MOL %-40s %4d\n", natoms, FILENAME,
0);
for(i=1;i<=natoms;i++) printf("%4d%4d% 9.4f% 9.4f% 9.4f\n", i,
atomtype(i), x[i], y[i], z[i]);
printf("%4d MOL\n", nbonds);
for(i=1;i<=nbonds;i++) printf("%4d%4d%4d %4d\n", i,
bond_a[i], bond_b[i], 1);
printf("%4d MOL\n", 0);
}
function dist(a1, a2) {
# get distance between a1-th and a2-th atom
d = (x[a1]-x[a2])*(x[a1]-x[a2]) + (y[a1]-y[a2])*(y[a1]-y[a2]) +
(z[a1]-z[a2])*(z[a1]-z[a2])
return sqrt(d)
}
# get sybyl atomtype of a1-th atom
function atomtype(a1) {
a = an[a1]; # atomic number
if (a == 0) return 26; # Dummy
else if (a == 1) return 13; # H
else if (a == 3) return 24; # Li
else if (a == 6) { # C
if (neighbors[a1] == 1) return 4;
else if (neighbors[a1] == 2) return 3;
else if (neighbors[a1] == 3) return 2;
else if (neighbors[a1] == 4) return 1;
else return 1;
} else if (a == 7) {
if (neighbors[a1] == 1) return 7; # N1
else if (neighbors[a1] == 2) return 6; # N2
else if (neighbors[a1] == 3) return 5; # N3
else if (neighbors[a1] == 4) return 31; # N3+
else return 5;
} else if (a == 8) {
if (neighbors[a1] == 1) return 9; # O2
else if (neighbors[a1] == 2) return 8; # O3
else return 8;
} else if (a == 9) return 16; # F
else if (a == 11) return 21; # Na
else if (a == 13) return 25; # Al
else if (a == 14) return 27; # Si
else if (a == 15) return 12; # P3
else if (a == 16) {
if (neighbors[a1] == 1) return 18; # S2
else if (neighbors[a1] == 2) return 10; # S3
else return 10;
} else if (a == 17) return 15; # Cl
else if (a == 19) return 22; # K
else if (a == 20) return 23; # Ca
else if (a == 35) return 14; # Br
else if (a == 53) return 17; # I
else return 26;
}
------------------- CUT HERE ----------------------------------------------