1 /*****
2 This file is part of the Babel Program
3 Copyright (C) 1992-96 W. Patrick Walters and Matthew T. Stahl
4 All Rights Reserved
5 All Rights Reserved
6 All Rights Reserved
7 All Rights Reserved
8
9 For more information please contact :
10
11 babel@mercury.aichem.arizona.edu
12 --------------------------------------------------------------------------------
13
14 FILE : int2cart.c
15 AUTHOR(S) : Pat Walters
16 DATE : 11-92
17 PURPOSE : Routines to convert internal to cartesian coordinates
18 Most of this was adapted from Ajay's routines and those found in MOPAC 5.
19 ******/
20
21
22
23 #include "bbltyp.h"
24
25 #define CONV PI/180.0
26 static warning wstr;
27
old_int_to_cart(ums_type * mol)28 int old_int_to_cart(ums_type *mol)
29 {
30
31 double *coord[4], *geo[4], *geol[4];
32 int *na, *nb, *nc;
33
34 int i;
35 double ccos,cosa;
36 int natoms;
37 int ma,mb,mc;
38 double xa, ya, za, xb, yb, zb;
39 double rbc;
40 double xyb;
41 int k;
42 double yza;
43 double xpa, xpb, ypa, xqa, zqa;
44 double cosph, sinph, costh, sinth, coskh, sinkh;
45 double cosd, sina, sind;
46 double xd, yd, zd;
47 double xpd, ypd, zpd;
48 double xqd, yqd, zqd;
49 double xrd;
50
51
52 natoms = Atoms;
53
54
55 for (i = 0; i <= 3; i++)
56 {
57 coord[i] = (double *)malloc((natoms + 1) * sizeof(double));
58 geo[i] = (double *)malloc((natoms + 1) * sizeof(double));
59 geol[i] = (double *)malloc((natoms + 1) * sizeof(double));
60 }
61
62 na = (int *)malloc((natoms + 1) * sizeof(int));
63 nb = (int *)malloc((natoms + 1) * sizeof(int));
64 nc = (int *)malloc((natoms + 1) * sizeof(int));
65
66 for (i = 1; i <= natoms; i++)
67 {
68 geo[1][i] = mol->internal[i].r;
69 geo[2][i] = mol->internal[i].w;
70 geo[3][i] = mol->internal[i].t;
71 na[i] = mol->internal[i].na;
72 nb[i] = mol->internal[i].nb;
73 nc[i] = mol->internal[i].nc;
74 }
75
76 for (i = 1; i <= natoms; i++)
77 {
78 geol[1][i] = geo[1][i];
79 geol[2][i] = geo[2][i] * CONV;
80 geol[3][i] = geo[3][i] * CONV;
81 }
82
83 if( natoms > 1 )
84 {
85 coord[1][2] = geol[1][2];
86 coord[2][2] = 0.0;
87 coord[3][2] = 0.0;
88 }
89
90 if( natoms > 2 )
91 {
92 ccos = cos(geol[2][3]);
93 if (na[3] == 1)
94 coord[1][3] = coord[1][1] + geol[1][3] * ccos;
95 else
96 coord[1][3] = coord[1][2] - geol[1][3] * ccos;
97 coord[2][3] = geol[1][3] * sin(geol[2][3]);
98 coord[3][3] = 0.0;
99 }
100
101 for (i = 4; i <= natoms; i++)
102 {
103 cosa = cos(geol[2][i]);
104 mb = nb[i];
105 mc = na[i];
106 xb = coord[1][mb] - coord[1][mc];
107 yb = coord[2][mb] - coord[2][mc];
108 zb = coord[3][mb] - coord[3][mc];
109 rbc = 1.0/sqrt(xb*xb + yb*yb + zb*zb);
110
111 if ( fabs(cosa) >= 0.9999999991)
112 {
113 rbc = geol[1][i] * rbc * cosa;
114 coord[1][i] = coord[1][mc] + xb * rbc;
115 coord[2][i] = coord[2][mc] + yb * rbc;
116 coord[3][i] = coord[3][mc] + zb * rbc;
117 }
118 else
119 {
120 ma = nc[i];
121 xa = coord[1][ma] - coord[1][mc];
122 ya = coord[2][ma] - coord[2][mc];
123 za = coord[3][ma] - coord[3][mc];
124
125 xyb = sqrt(xb*xb + yb*yb);
126 k = -1;
127 if (xyb <= 0.10)
128 {
129 xpa = za;
130 za = -xa;
131 xa = xpa;
132 xpb = zb;
133 zb = -xb;
134 xb = xpb;
135 xyb = sqrt(xb*xb + yb*yb);
136 k = 1;
137 }
138 costh = xb/xyb;
139 sinth = yb/xyb;
140 xpa = xa * costh + ya * sinth;
141 ypa = ya * costh - xa * sinth;
142 sinph = zb * rbc;
143 cosph = sqrt(fabs(1.00 - sinph * sinph));
144 xqa = xpa * cosph + za * sinph;
145 zqa = za * cosph - xpa * sinph;
146
147 yza = sqrt(ypa*ypa + zqa*zqa);
148 if ((yza < 1.0E-1) && (yza > 1.0E-10))
149 {
150 sprintf(wstr,"THE FAULTY ATOM IS ATOM NUMBER %d",i);
151 show_warning(wstr);
152 sprintf(wstr,"ATOMS %d, %d AND %d ARE WITHIN %f ANGSTROMS OF A STRAIGHT LINE",mc,mb,ma,yza);
153 show_warning(wstr);
154 }
155 coskh = ypa/yza;
156 sinkh = zqa/yza;
157 if (yza < 1.0E-10)
158 {
159 coskh = 1.0;
160 sinkh = 0.0;
161 }
162
163 sina = sin(geol[2][i]);
164 sind = -sin(geol[3][i]);
165 cosd = cos(geol[3][i]);
166 xd = geol[1][i] * cosa;
167 yd = geol[1][i] * sina * cosd;
168 zd = geol[1][i] * sina * sind;
169
170 ypd = yd * coskh - zd * sinkh;
171 zpd = zd * coskh + yd * sinkh;
172 xpd = xd * cosph - zpd * sinph;
173 zqd = zpd * cosph + xd * sinph;
174 xqd = xpd * costh - ypd * sinth;
175 yqd = ypd * costh + xpd * sinth;
176 if (k >= 1)
177 {
178 xrd = -zqd;
179 zqd = xqd;
180 xqd = xrd;
181 }
182 coord[1][i] = xqd + coord[1][mc];
183 coord[2][i] = yqd + coord[2][mc];
184 coord[3][i] = zqd + coord[3][mc];
185 }
186 }
187 for (i = 1; i <= natoms; i++)
188 {
189 X(i) = coord[1][i];
190 Y(i) = coord[2][i];
191 Z(i) = coord[3][i];
192 }
193
194 for (i = 0; i <= 3; i++)
195 {
196 free(coord[i]);
197 free(geo[i]);
198 free(geol[i]);
199 }
200 free(na);
201 free(nb);
202 free(nc);
203 return(TRUE);
204 }
205
206