| 1 |
#define EXTERN extern |
| 2 |
|
| 3 |
#include "pcwin.h" |
| 4 |
#include "pcmod.h" |
| 5 |
#include "utility.h" |
| 6 |
#include "asnsym.h" |
| 7 |
|
| 8 |
double TOLER = 0.002; |
| 9 |
double TOL2 = 0.00002; |
| 10 |
|
| 11 |
double SIGN(double,double); |
| 12 |
static void vibsym(int); |
| 13 |
static int tstc3(double **,double **,int,int,int,double *); |
| 14 |
static int tstc4(double **,double **,int,int,int,double *); |
| 15 |
static int tstc5(double **,double **,int,int,int,double *); |
| 16 |
static void secmom(double **, double *, double *, double *, int *, double (*)[3]); |
| 17 |
static void orients(double **, double **,double *, double *,int,int); |
| 18 |
static void tcopy(double **, double **,int); |
| 19 |
static int findcn(double **,double **,double *, double *, int *); |
| 20 |
static int findc2(double **,double **, double **, int *); |
| 21 |
static int findv(double **,double **, double **, int *); |
| 22 |
static int findsn(double **,double **, double **, int *); |
| 23 |
static void tform(double (*)[3], double **, double **,int); |
| 24 |
static void reflect(double **, double **,int); |
| 25 |
static int tcomp(double **,double **); |
| 26 |
static void circset(double **, double **,int,int *,int *); |
| 27 |
static void rotates(double **a,double **b,int,double); |
| 28 |
static void sphere(double **,double **, double **,int *, int *); |
| 29 |
static void sym_invert(double **, double **); |
| 30 |
static void triang(double **,int,int,int,double *, double *,double *,double *,double *,double *); |
| 31 |
static void zero_array(double **,int); |
| 32 |
static void hqrii(double *,int,int,double *,double (*)[3]); |
| 33 |
static void sphset(double **, double **,int *, int *,int *); |
| 34 |
|
| 35 |
// ======================== |
| 36 |
void zero_array(double **a,int num) |
| 37 |
{ |
| 38 |
int i; |
| 39 |
for (i=1; i <= num; i++) |
| 40 |
{ |
| 41 |
a[i][0] = 0.0; |
| 42 |
a[i][1] = 0.0; |
| 43 |
a[i][2] = 0.0; |
| 44 |
} |
| 45 |
} |
| 46 |
void triang(double **a,int iat,int jat,int kat,double *alpha, double *beta, |
| 47 |
double *gamma,double *dij, double *dik, double *djk) |
| 48 |
{ |
| 49 |
double xi,yi,zi,xj,yj,zj,xk,yk,zk,dotjk,dotik,dotij; |
| 50 |
|
| 51 |
xi = a[iat][0]; |
| 52 |
yi = a[iat][1]; |
| 53 |
zi = a[iat][2]; |
| 54 |
xj = a[jat][0]; |
| 55 |
yj = a[jat][1]; |
| 56 |
zj = a[jat][2]; |
| 57 |
xk = a[kat][0]; |
| 58 |
yk = a[kat][1]; |
| 59 |
zk = a[kat][2]; |
| 60 |
*dij = sqrt( (xi-xj)*(xi-xj) + (yi-yj)*(yi-yj) + (zi-zj)*(zi-zj) ); |
| 61 |
*dik = sqrt( (xi-xk)*(xi-xk) + (yi-yk)*(yi-yk) + (zi-zk)*(zi-zk) ); |
| 62 |
*djk = sqrt( (xj-xk)*(xj-xk) + (yj-yk)*(yj-yk) + (zj-zk)*(zj-zk) ); |
| 63 |
dotjk = (xj-xi)*(xk-xi) + (yj-yi)*(yk-yi) + (zj-zi)*(zk-zi); |
| 64 |
dotik = (xi-xj)*(xk-xj) + (yi-yj)*(yk-yj) + (zi-zj)*(zk-zj); |
| 65 |
dotij = (xi-xk)*(xj-xk) + (yi-yk)*(yj-yk) + (zi-zk)*(zj-zk); |
| 66 |
*alpha = acos(dotjk/(*dij**dik)); |
| 67 |
*beta = acos(dotik/(*dij**djk)); |
| 68 |
*gamma = acos(dotij/(*dik**djk)); |
| 69 |
|
| 70 |
} |
| 71 |
/* ================================================= */ |
| 72 |
int tstc3(double **a,double **b,int iat,int jat,int kat,double *center) |
| 73 |
{ |
| 74 |
int itst, iz; |
| 75 |
double p1[3]; |
| 76 |
double phi3,theta3,fact3; |
| 77 |
double alpha,beta,gamma,dij,dik,djk; |
| 78 |
double px,py,pz; |
| 79 |
|
| 80 |
p1[0] = p1[1] = p1[2] = 0.0; |
| 81 |
itst = 0; |
| 82 |
phi3 = (8.0/3.0)*atan(1.0); |
| 83 |
theta3 = 0.5*phi3; |
| 84 |
fact3 = 2.0/3.0; |
| 85 |
|
| 86 |
triang(a,iat,jat,kat,&alpha,&beta,&gamma,&dij,&dik,&djk); |
| 87 |
|
| 88 |
if (fabs(alpha-theta3) < TOLER && fabs(dij-dik) < TOLER) |
| 89 |
{ |
| 90 |
px = 0.5*(a[jat][0] + a[kat][0]); |
| 91 |
py = 0.5*(a[jat][1] + a[kat][1]); |
| 92 |
pz = 0.5*(a[jat][2] + a[kat][2]); |
| 93 |
center[0] = a[iat][0] + fact3*(px - a[iat][0]); |
| 94 |
center[1] = a[iat][1] + fact3*(py - a[iat][1]); |
| 95 |
center[2] = a[iat][2] + fact3*(pz - a[iat][2]); |
| 96 |
|
| 97 |
orients(a,b,p1,center,3,natom); |
| 98 |
tcopy(b,a,natom); |
| 99 |
rotates(a,b,3,phi3); |
| 100 |
iz = tcomp(a,b); |
| 101 |
return (iz); |
| 102 |
} |
| 103 |
return FALSE; |
| 104 |
} |
| 105 |
int tstc4(double **a,double **b,int iat,int jat,int kat,double *center) |
| 106 |
{ |
| 107 |
int itst, iz; |
| 108 |
double p1[3]; |
| 109 |
double alpha,beta,gamma,dij,dik,djk; |
| 110 |
double halfpi; |
| 111 |
|
| 112 |
p1[0] = p1[1] = p1[2] = 0.0; |
| 113 |
itst = 0; |
| 114 |
halfpi = 2.0*atan(1.0); |
| 115 |
|
| 116 |
triang(a,iat,jat,kat,&alpha,&beta,&gamma,&dij,&dik,&djk); |
| 117 |
|
| 118 |
if (fabs(alpha-halfpi) < TOLER && fabs(dij-dik) < TOLER) |
| 119 |
{ |
| 120 |
center[0] = 0.5*(a[jat][0] + a[kat][0]); |
| 121 |
center[1] = 0.5*(a[jat][1] + a[kat][1]); |
| 122 |
center[2] = 0.5*(a[jat][2] + a[kat][2]); |
| 123 |
} else if (fabs(beta-halfpi) < TOLER && fabs(dij-djk) < TOLER) |
| 124 |
{ |
| 125 |
center[0] = 0.5*(a[iat][0] + a[kat][0]); |
| 126 |
center[1] = 0.5*(a[iat][1] + a[kat][1]); |
| 127 |
center[2] = 0.5*(a[iat][2] + a[kat][2]); |
| 128 |
} else if (fabs(gamma-halfpi) < TOLER && fabs(dik-djk) < TOLER) |
| 129 |
{ |
| 130 |
center[0] = 0.5*(a[jat][0] + a[iat][0]); |
| 131 |
center[1] = 0.5*(a[jat][1] + a[iat][1]); |
| 132 |
center[2] = 0.5*(a[jat][2] + a[iat][2]); |
| 133 |
} else if (fabs(dij-dik) < TOLER && fabs(dij-djk) < TOLER && TOLER && fabs(dik-djk)) |
| 134 |
{ |
| 135 |
center[0] = a[iat][0]; |
| 136 |
center[1] = a[iat][1]; |
| 137 |
center[2] = a[iat][2]; |
| 138 |
} else |
| 139 |
return FALSE; |
| 140 |
|
| 141 |
orients(a,b,p1,center,3,natom); |
| 142 |
tcopy(b,a,natom); |
| 143 |
rotates(a,b,3,halfpi); |
| 144 |
iz = tcomp(a,b); |
| 145 |
return (iz); |
| 146 |
|
| 147 |
} |
| 148 |
int tstc5(double **a,double **b,int iat,int jat,int kat,double *center) |
| 149 |
{ |
| 150 |
int itst, iz; |
| 151 |
double p1[3]; |
| 152 |
double piovr4, phi5,theta5,fact5; |
| 153 |
double alpha,beta,gamma,dij,dik,djk; |
| 154 |
double px,py,pz; |
| 155 |
|
| 156 |
p1[0] = p1[1] = p1[2] = 0.0; |
| 157 |
itst = 0; |
| 158 |
piovr4 = atan(1.0); |
| 159 |
phi5 = 1.6*piovr4; |
| 160 |
theta5 = 2.4*piovr4; |
| 161 |
fact5 = 1.0/(2.0*sin(0.8*piovr4)*sin(0.8*piovr4)); |
| 162 |
|
| 163 |
triang(a,iat,jat,kat,&alpha,&beta,&gamma,&dij,&dik,&djk); |
| 164 |
|
| 165 |
if (fabs(alpha-theta5) < TOLER && fabs(dij-dik) < TOLER) |
| 166 |
{ |
| 167 |
px = 0.5*(a[jat][0] + a[kat][0]); |
| 168 |
py = 0.5*(a[jat][1] + a[kat][1]); |
| 169 |
pz = 0.5*(a[jat][2] + a[kat][2]); |
| 170 |
center[0] = a[iat][0] + fact5*(px - a[iat][0]); |
| 171 |
center[1] = a[iat][1] + fact5*(py - a[iat][1]); |
| 172 |
center[2] = a[iat][2] + fact5*(pz - a[iat][2]); |
| 173 |
} else if (fabs(beta-theta5) < TOLER && fabs(dij-djk) < TOLER) |
| 174 |
{ |
| 175 |
px = 0.5*(a[iat][0] + a[kat][0]); |
| 176 |
py = 0.5*(a[iat][1] + a[kat][1]); |
| 177 |
pz = 0.5*(a[iat][2] + a[kat][2]); |
| 178 |
center[0] = a[jat][0] + fact5*(px - a[jat][0]); |
| 179 |
center[1] = a[jat][1] + fact5*(py - a[jat][1]); |
| 180 |
center[2] = a[jat][2] + fact5*(pz - a[jat][2]); |
| 181 |
} else if (fabs(gamma-theta5) < TOLER && fabs(dik-djk) < TOLER) |
| 182 |
{ |
| 183 |
px = 0.5*(a[jat][0] + a[iat][0]); |
| 184 |
py = 0.5*(a[jat][1] + a[iat][1]); |
| 185 |
pz = 0.5*(a[jat][2] + a[iat][2]); |
| 186 |
center[0] = a[kat][0] + fact5*(px - a[kat][0]); |
| 187 |
center[1] = a[kat][1] + fact5*(py - a[kat][1]); |
| 188 |
center[2] = a[kat][2] + fact5*(pz - a[kat][2]); |
| 189 |
} else |
| 190 |
return FALSE; |
| 191 |
orients(a,b,p1,center,3,natom); |
| 192 |
tcopy(b,a,natom); |
| 193 |
rotates(a,b,3,phi5); |
| 194 |
iz = tcomp(a,b); |
| 195 |
return (iz); |
| 196 |
} |
| 197 |
/* ====================================== */ |
| 198 |
void sphset(double **a, double **aset,int *npop, int *nset,int *numset) |
| 199 |
{ |
| 200 |
int i,j1,j,iset, num; |
| 201 |
double xtmp,ytmp,ztmp; |
| 202 |
|
| 203 |
for (i=1; i <= natom; i++) |
| 204 |
{ |
| 205 |
aset[i][0] = i; |
| 206 |
aset[i][1] = atom.atomwt[i]; |
| 207 |
xtmp = a[i][0]; |
| 208 |
ytmp = a[i][1]; |
| 209 |
ztmp = a[i][2]; |
| 210 |
aset[i][2] = sqrt(xtmp*xtmp+ytmp*ytmp+ztmp*ztmp); |
| 211 |
} |
| 212 |
iset = 0; |
| 213 |
for (i=1; i < natom; i++) |
| 214 |
{ |
| 215 |
num = 0; |
| 216 |
if (aset[i][0] != 0.0) |
| 217 |
{ |
| 218 |
iset++; |
| 219 |
num++; |
| 220 |
nset[i] = iset; |
| 221 |
aset[i][0] = 0.0; |
| 222 |
j1 = i + 1; |
| 223 |
for (j=j1; j <= natom; j++) |
| 224 |
{ |
| 225 |
if (aset[j][0] != 0.0) |
| 226 |
{ |
| 227 |
if (fabs(aset[j][1]-aset[i][1]) <= TOLER && |
| 228 |
fabs(aset[j][2]-aset[i][2]) <= TOLER) |
| 229 |
{ |
| 230 |
num++; |
| 231 |
nset[j] = iset; |
| 232 |
aset[j][0] = 0.0; |
| 233 |
} |
| 234 |
} |
| 235 |
} |
| 236 |
npop[iset] = num; |
| 237 |
} |
| 238 |
} |
| 239 |
*numset = iset; |
| 240 |
} |
| 241 |
void sphere(double **a, double **b, double **c, int *nset, int *norder) |
| 242 |
{ |
| 243 |
int i,j,k,keyset, numset, minnum; |
| 244 |
int mpop, num3, i2, j1, j2, k1, iat=0,jat,kat; |
| 245 |
int iz,iz1,iz2,key,ihop,ixyz; |
| 246 |
double savez,curz,x,y,z,cmax,theta,tst; |
| 247 |
int *npop; |
| 248 |
double halfpi,pi; |
| 249 |
double p1[3],center[3],save1[3],save2[3]; |
| 250 |
|
| 251 |
halfpi = 2.0*atan(1.0); |
| 252 |
pi = 2.0*halfpi; |
| 253 |
p1[0] = p1[1] = p1[2] = 0.0; |
| 254 |
save2[0] = save2[1] = save2[2] = 0.0; |
| 255 |
savez = 0.0; |
| 256 |
jat = 0; |
| 257 |
keyset = 0; |
| 258 |
npop = ivector(0,natom+1); |
| 259 |
sphset(a,c,npop,nset,&numset); |
| 260 |
|
| 261 |
minnum = natom; |
| 262 |
*norder = 0; |
| 263 |
for (i=1; i <= numset; i++) |
| 264 |
{ |
| 265 |
j = npop[i]; |
| 266 |
if (npop[i] <= minnum && npop[i] > 3) |
| 267 |
{ |
| 268 |
minnum = npop[i]; |
| 269 |
keyset = i; |
| 270 |
} |
| 271 |
} |
| 272 |
|
| 273 |
key = natom; |
| 274 |
mpop = 0; |
| 275 |
for (i=1; i <= natom; i++) |
| 276 |
{ |
| 277 |
if (nset[i] == keyset) |
| 278 |
{ |
| 279 |
mpop++; |
| 280 |
if (key < i) |
| 281 |
key = i; |
| 282 |
nset[mpop] = i; |
| 283 |
} |
| 284 |
} |
| 285 |
|
| 286 |
if (mpop < 3) |
| 287 |
{ |
| 288 |
*norder = 0; |
| 289 |
free_ivector(npop,0,natom+1); |
| 290 |
return; |
| 291 |
} |
| 292 |
num3 = 0; |
| 293 |
i2 = mpop-2; |
| 294 |
j2 = mpop-1; |
| 295 |
for (i=1; i <= i2; i++) |
| 296 |
{ |
| 297 |
iat = nset[i]; |
| 298 |
j1 = i + 1; |
| 299 |
for (j=j1; j <= j2; j++) |
| 300 |
{ |
| 301 |
jat = nset[j]; |
| 302 |
k1 = j+1; |
| 303 |
for (k=k1; k <= mpop; k++) |
| 304 |
{ |
| 305 |
kat = nset[k]; |
| 306 |
if (*norder <= 3) |
| 307 |
{ |
| 308 |
tcopy(a,c,natom); |
| 309 |
iz = tstc5(c,b,iat,jat,kat,center); |
| 310 |
if (iz == TRUE) |
| 311 |
{ |
| 312 |
*norder = 5; |
| 313 |
savez = fabs(c[key][2]); |
| 314 |
save1[0] = center[0]; |
| 315 |
save1[1] = center[1]; |
| 316 |
save1[2] = center[2]; |
| 317 |
}else |
| 318 |
{ |
| 319 |
tcopy(a,c,natom); |
| 320 |
zero_array(b,natom); |
| 321 |
iz1 = tstc4(c,b,iat,jat,kat,center); |
| 322 |
if (iz1 == TRUE) |
| 323 |
{ |
| 324 |
*norder = 4; |
| 325 |
tcopy(c,a,natom); |
| 326 |
} else |
| 327 |
{ |
| 328 |
if (num3 != 2) |
| 329 |
{ |
| 330 |
tcopy(a,c,natom); |
| 331 |
iz2 = tstc3(c,b,iat,jat,kat,center); |
| 332 |
if (iz2 == TRUE) |
| 333 |
{ |
| 334 |
ihop = num3+1; |
| 335 |
if (ihop == 1) |
| 336 |
{ |
| 337 |
*norder = 3; |
| 338 |
num3 = 1; |
| 339 |
save1[0] = center[0]; |
| 340 |
save1[1] = center[1]; |
| 341 |
save1[2] = center[2]; |
| 342 |
} else if (ihop == 2) |
| 343 |
{ |
| 344 |
num3 = 2; |
| 345 |
save2[0] = center[0]; |
| 346 |
save2[1] = center[1]; |
| 347 |
save2[2] = center[2]; |
| 348 |
} |
| 349 |
} |
| 350 |
} |
| 351 |
} |
| 352 |
} |
| 353 |
} else if (*norder == 5) |
| 354 |
{ |
| 355 |
tcopy(a,c,natom); |
| 356 |
iz = tstc5(c,b,iat,jat,kat,center); |
| 357 |
if (iz == TRUE) |
| 358 |
{ |
| 359 |
curz = fabs(c[key][2]); |
| 360 |
if (fabs(curz-savez) > TOLER && (savez < curz )) |
| 361 |
{ |
| 362 |
savez = curz; |
| 363 |
save1[0] = center[0]; |
| 364 |
save1[1] = center[1]; |
| 365 |
save1[2] = center[2]; |
| 366 |
} |
| 367 |
} |
| 368 |
} else if (*norder == 4) |
| 369 |
{ |
| 370 |
tcopy(a,c,natom); |
| 371 |
iz = tstc4(c,b,iat,jat,kat,center); |
| 372 |
if (iz == TRUE) |
| 373 |
{ |
| 374 |
if (fabs(center[2]) < TOLER) |
| 375 |
{ |
| 376 |
x = center[0]; |
| 377 |
y = center[1]; |
| 378 |
theta = halfpi; |
| 379 |
if (fabs(y) > TOLER) theta = atan(x/y); |
| 380 |
rotates(c,a,3,theta); |
| 381 |
goto L_5; |
| 382 |
} |
| 383 |
} |
| 384 |
} |
| 385 |
} |
| 386 |
} |
| 387 |
} |
| 388 |
L_5: |
| 389 |
if (*norder == 0) |
| 390 |
{ |
| 391 |
center[0] = 0.5*(a[iat][0] + a[jat][0]); |
| 392 |
center[1] = 0.5*(a[iat][1] + a[jat][1]); |
| 393 |
center[2] = 0.5*(a[iat][2] + a[jat][2]); |
| 394 |
orients(a,b,p1,center,3,natom); |
| 395 |
x = b[iat][0]; |
| 396 |
y = b[iat][1]; |
| 397 |
theta = halfpi; |
| 398 |
if (fabs(y) > TOLER) theta = atan(x/y); |
| 399 |
rotates(b,a,3,theta); |
| 400 |
} else if (*norder == 3) |
| 401 |
{ |
| 402 |
center[0] = 0.5*(save1[0] + save2[0]); |
| 403 |
center[1] = 0.5*(save1[1] + save2[1]); |
| 404 |
center[2] = 0.5*(save1[2] + save2[2]); |
| 405 |
orients(a,b,p1,center,3,natom); |
| 406 |
tcopy(b,a,natom); |
| 407 |
b[0][0] = save1[0]; |
| 408 |
b[0][1] = save1[1]; |
| 409 |
b[0][2] = save1[2]; |
| 410 |
b[1][0] = save2[0]; |
| 411 |
b[1][1] = save2[1]; |
| 412 |
b[1][2] = save2[2]; |
| 413 |
orients(b,c,p1,center,3,2); |
| 414 |
tcopy(c,b,2); |
| 415 |
save1[0] = b[0][0]; |
| 416 |
save1[1] = b[0][1]; |
| 417 |
save1[2] = b[0][2]; |
| 418 |
save2[0] = b[1][0]; |
| 419 |
save2[1] = b[1][1]; |
| 420 |
save2[2] = b[1][2]; |
| 421 |
x = 0.5*(save1[0]-save2[1]); |
| 422 |
y = 0.5*(save1[1]+save2[0]); |
| 423 |
theta = halfpi; |
| 424 |
if (fabs(y) > TOLER) theta = atan(x/y); |
| 425 |
rotates(a,b,3,theta); |
| 426 |
tcopy(b,a,natom); |
| 427 |
|
| 428 |
} else if (*norder == 4) |
| 429 |
{ |
| 430 |
x = a[key][0]; |
| 431 |
y = a[key][1]; |
| 432 |
z = a[key][2]; |
| 433 |
cmax = fabs(z); |
| 434 |
ixyz = 3; |
| 435 |
for (i=1; i <= 2; i++) |
| 436 |
{ |
| 437 |
tst = fabs(a[key][i]); |
| 438 |
if (fabs(tst-cmax) > TOLER && (cmax <= tst)) |
| 439 |
{ |
| 440 |
ixyz = i; |
| 441 |
cmax = tst; |
| 442 |
} |
| 443 |
} |
| 444 |
if (ixyz < 3) |
| 445 |
{ |
| 446 |
if (ixyz == 2) |
| 447 |
ixyz = 1; |
| 448 |
else |
| 449 |
ixyz = 2; |
| 450 |
rotates(a,b,ixyz,halfpi); |
| 451 |
tcopy(b,a,natom); |
| 452 |
} |
| 453 |
if (a[key][3] < 0.0) |
| 454 |
{ |
| 455 |
rotates(a,b,1,pi); |
| 456 |
tcopy(b,a,natom); |
| 457 |
} |
| 458 |
} else if (*norder == 5) |
| 459 |
{ |
| 460 |
orients(a,b,p1,save1,3,natom); |
| 461 |
tcopy(b,a,natom); |
| 462 |
if (a[key][3] < 0.0) |
| 463 |
{ |
| 464 |
rotates(a,b,1,pi); |
| 465 |
tcopy(b,a,natom); |
| 466 |
} |
| 467 |
} |
| 468 |
|
| 469 |
|
| 470 |
free_ivector(npop,0,natom+1); |
| 471 |
} |
| 472 |
void sym_invert(double **a, double **b) |
| 473 |
{ |
| 474 |
int i,j; |
| 475 |
double t[3][3]; |
| 476 |
for (i=0; i < 3; i++) |
| 477 |
{ |
| 478 |
for (j=0; j <3; j++) |
| 479 |
t[i][j] = 0.0; |
| 480 |
} |
| 481 |
t[0][0] = -1.0; |
| 482 |
t[1][1] = -1.0; |
| 483 |
t[2][2] = -1.0; |
| 484 |
tform(t,a,b,natom); |
| 485 |
} |
| 486 |
void reflect(double **a, double **b, int iaxis) |
| 487 |
{ |
| 488 |
int i,j; |
| 489 |
double t[3][3]; |
| 490 |
for (i=0; i < 3; i++) |
| 491 |
{ |
| 492 |
for (j=0; j <3; j++) |
| 493 |
t[i][j] = 0.0; |
| 494 |
t[i][i] = 1.0; |
| 495 |
} |
| 496 |
t[iaxis-1][iaxis-1] = -1.0; |
| 497 |
tform(t,a,b,natom); |
| 498 |
} |
| 499 |
/* ==================================== */ |
| 500 |
int tcomp(double **a, double **b) |
| 501 |
{ |
| 502 |
int i,j, k; |
| 503 |
int *nrs; |
| 504 |
double dsum, x1, x2; |
| 505 |
|
| 506 |
nrs = ivector(0,natom+1); |
| 507 |
for (i=0; i <= natom; i++) |
| 508 |
nrs[i] = 0; |
| 509 |
|
| 510 |
for (i=1; i <= natom; i++) |
| 511 |
{ |
| 512 |
for (j=1; j <= natom; j++) |
| 513 |
{ |
| 514 |
if (fabs(atom.atomwt[i]-atom.atomwt[j]) <= TOLER) |
| 515 |
{ |
| 516 |
if (nrs[i] == 0) |
| 517 |
{ |
| 518 |
dsum = 0; |
| 519 |
for (k=0; k < 3; k++) |
| 520 |
{ |
| 521 |
x1 = a[i][k]; |
| 522 |
x2 = b[j][k]; |
| 523 |
dsum += (a[i][k]-b[j][k])*(a[i][k]-b[j][k]); |
| 524 |
} |
| 525 |
if (dsum <= TOLER) |
| 526 |
{ |
| 527 |
nrs[i] = j; |
| 528 |
goto L_5; |
| 529 |
} |
| 530 |
} |
| 531 |
} |
| 532 |
} |
| 533 |
free_ivector(nrs,0,natom+1); |
| 534 |
return FALSE; |
| 535 |
L_5: |
| 536 |
continue; |
| 537 |
} |
| 538 |
free_ivector(nrs,0,natom+1); |
| 539 |
return TRUE; |
| 540 |
} |
| 541 |
void tcopy(double **a, double **b, int num) |
| 542 |
{ |
| 543 |
int i,j; |
| 544 |
for (i=1; i <= num; i++) |
| 545 |
{ |
| 546 |
for (j=0; j < 3; j++) |
| 547 |
b[i][j] = a[i][j]; |
| 548 |
} |
| 549 |
} |
| 550 |
void tform(double t[][3], double **a, double **b, int num) |
| 551 |
{ |
| 552 |
int i; |
| 553 |
for (i=1; i <= num; i++) |
| 554 |
{ |
| 555 |
b[i][0] = t[0][0]*a[i][0] + t[0][1]*a[i][1] + t[0][2]*a[i][2]; |
| 556 |
b[i][1] = t[1][0]*a[i][0] + t[1][1]*a[i][1] + t[1][2]*a[i][2]; |
| 557 |
b[i][2] = t[2][0]*a[i][0] + t[2][1]*a[i][1] + t[2][2]*a[i][2]; |
| 558 |
} |
| 559 |
} |
| 560 |
/* ==================================== */ |
| 561 |
void circset(double **a, double **aset,int ixyz,int *nset,int *numset) |
| 562 |
{ |
| 563 |
int i,j, i1, i2, i3, iset,j1; |
| 564 |
double q2, q3, an,ap,ad; |
| 565 |
|
| 566 |
i1 = ixyz; |
| 567 |
i2 = ixyz%3+1; |
| 568 |
i3 = i2%3 +1; |
| 569 |
i1--; |
| 570 |
i2--; |
| 571 |
i3--; |
| 572 |
for (i=1; i <= natom; i++) |
| 573 |
{ |
| 574 |
aset[i][0] = atom.atomwt[i]; |
| 575 |
aset[i][1] = a[i][i1]; |
| 576 |
q2 = a[i][i2]; |
| 577 |
q3 = a[i][i3]; |
| 578 |
aset[i][2] = sqrt(q2*q2+q3*q3); |
| 579 |
} |
| 580 |
// set 0 on axis atoms |
| 581 |
for (i=1; i <= natom; i++) |
| 582 |
{ |
| 583 |
if (fabs(aset[i][2]) <= TOLER) |
| 584 |
{ |
| 585 |
nset[i] = 0; |
| 586 |
aset[i][0] = 0.0; |
| 587 |
} |
| 588 |
} |
| 589 |
// remaining sets |
| 590 |
iset = 0; |
| 591 |
for (i=1; i < natom; i++) // loop is natom-1 |
| 592 |
{ |
| 593 |
if (aset[i][0] != 0.0) |
| 594 |
{ |
| 595 |
iset++; |
| 596 |
nset[i] = iset; |
| 597 |
an = aset[i][0]; |
| 598 |
ap = aset[i][1]; |
| 599 |
ad = aset[i][2]; |
| 600 |
aset[i][0] = 0.0; |
| 601 |
j1 = i+1; |
| 602 |
for (j=j1; j <= natom; j++) |
| 603 |
{ |
| 604 |
if (fabs(aset[j][0]-an) <= TOL2 && |
| 605 |
fabs(aset[j][1]-ap) <= TOLER && |
| 606 |
fabs(aset[j][2]-ad) <= TOLER) |
| 607 |
{ |
| 608 |
nset[j] = iset; |
| 609 |
aset[j][0] = 0.0; |
| 610 |
} |
| 611 |
} |
| 612 |
} |
| 613 |
} |
| 614 |
*numset = iset; |
| 615 |
for (i=1; i <= natom; i++) |
| 616 |
aset[i][0] = atom.atomwt[i]; |
| 617 |
} |
| 618 |
/* =============================================== */ |
| 619 |
void rotates(double **a,double **b,int iaxis,double theta) |
| 620 |
{ |
| 621 |
int i,j,i1,i2,i3; |
| 622 |
double s,c; |
| 623 |
double t[3][3]; |
| 624 |
|
| 625 |
i1 = iaxis; |
| 626 |
i2 = iaxis%3 +1; |
| 627 |
i3 = i2%3 + 1; |
| 628 |
i1--; |
| 629 |
i2--; |
| 630 |
i3--; |
| 631 |
s = sin(theta); |
| 632 |
c = cos(theta); |
| 633 |
for (i=0; i < 3; i++) |
| 634 |
{ |
| 635 |
for (j=0; j < 3; j++) |
| 636 |
t[i][j] = 0.0; |
| 637 |
} |
| 638 |
t[i1][i1] = 1.0; |
| 639 |
t[i2][i2] = c; |
| 640 |
t[i2][i3] = s; |
| 641 |
t[i3][i2] = -s; |
| 642 |
t[i3][i3] = c; |
| 643 |
tform(t,a,b,natom); |
| 644 |
} |
| 645 |
/* =============================================== */ |
| 646 |
void orients(double **a, double **b, double *p1, double *p2, int nset, int num) |
| 647 |
{ |
| 648 |
int i,j, i1, i2, i3; |
| 649 |
double trvec[3], t[3][3]; |
| 650 |
double v1, v2, v3, vnorm; |
| 651 |
double alph, beta, gamm, v2v2, v3v3, v2233; |
| 652 |
double tol = 1.0e-17; |
| 653 |
|
| 654 |
trvec[0] = p1[0]; |
| 655 |
trvec[1] = p1[1]; |
| 656 |
trvec[2] = p1[2]; |
| 657 |
|
| 658 |
for (i=1; i <= num; i++) |
| 659 |
{ |
| 660 |
for (j=0; j < 3; j++) |
| 661 |
a[i][j] -= trvec[j]; |
| 662 |
} |
| 663 |
|
| 664 |
i1 = nset; |
| 665 |
i2 = i1%3 + 1; |
| 666 |
i3 = i2%3 + 1; |
| 667 |
i1--; |
| 668 |
i2--; |
| 669 |
i3--; |
| 670 |
v1 = p2[i1] - p1[i1]; |
| 671 |
v2 = p2[i2] - p1[i2]; |
| 672 |
v3 = p2[i3] - p1[i3]; |
| 673 |
vnorm = sqrt(v1*v1 + v2*v2 + v3*v3); |
| 674 |
if (vnorm == 0.0) |
| 675 |
return; |
| 676 |
if ( fabs(v1) < TOLER && fabs(v2) < TOLER && fabs(v3) < TOLER) |
| 677 |
return; |
| 678 |
if (fabs(v2) < tol && fabs(v3) < tol) |
| 679 |
{ |
| 680 |
if (v1 > 0.0) |
| 681 |
{ |
| 682 |
t[i1][i1] = 1.0; |
| 683 |
t[i2][i2] = 1.0; |
| 684 |
t[i3][i3] = 1.0; |
| 685 |
t[i1][i2] = t[i2][i1] = 0.0; |
| 686 |
t[i1][i3] = t[i3][i1] = 0.0; |
| 687 |
t[i3][i2] = t[i2][i3] = 0.0; |
| 688 |
} else |
| 689 |
{ |
| 690 |
t[i1][i1] = -1.0; |
| 691 |
t[i2][i2] = 1.0; |
| 692 |
t[i3][i3] = 1.0; |
| 693 |
t[i1][i2] = t[i2][i1] = 0.0; |
| 694 |
t[i1][i3] = t[i3][i1] = 0.0; |
| 695 |
t[i3][i2] = t[i2][i3] = 0.0; |
| 696 |
} |
| 697 |
} else |
| 698 |
{ |
| 699 |
alph = v1/vnorm; |
| 700 |
beta = v2/vnorm; |
| 701 |
gamm = v3/vnorm; |
| 702 |
v2v2 = v2*v2; |
| 703 |
v3v3 = v3*v3; |
| 704 |
v2233 = 1.0/(v2v2+v3v3); |
| 705 |
t[i1][i1] = alph; |
| 706 |
t[i1][i2] = beta; |
| 707 |
t[i1][i3] = gamm; |
| 708 |
t[i2][i1] = -t[i1][i2]; |
| 709 |
t[i3][i1] = -t[i1][i3]; |
| 710 |
t[i2][i3] = v2*v3*(alph-1.0)*v2233; |
| 711 |
t[i3][i2] = t[i2][i3]; |
| 712 |
t[i2][i2] = (v2v2*alph +v3v3)*v2233; |
| 713 |
t[i3][i3] = (v3v3*alph +v2v2)*v2233; |
| 714 |
} |
| 715 |
tform(t,a,b,num); |
| 716 |
for (i=1; i <= num; i++) |
| 717 |
{ |
| 718 |
for (j=0; j < 3; j++) |
| 719 |
b[i][j] += trvec[j]; |
| 720 |
} |
| 721 |
} |
| 722 |
/* ================================== */ |
| 723 |
int findi(double **a, double **b) |
| 724 |
{ |
| 725 |
int iz; |
| 726 |
sym_invert(a,b); |
| 727 |
iz = tcomp(a,b); |
| 728 |
return (iz); |
| 729 |
} |
| 730 |
int findh(double **a, double **b, int nset) |
| 731 |
{ |
| 732 |
int iz; |
| 733 |
reflect(a,b,nset); |
| 734 |
iz = tcomp(a,b); |
| 735 |
return (iz); |
| 736 |
} |
| 737 |
int findcn(double **a, double **b,double *p1, double *p2, int *maxcn) |
| 738 |
{ |
| 739 |
int i,j,k,l, iz, iaxis; |
| 740 |
int i1, i2; |
| 741 |
double halfpi, pi,theta; |
| 742 |
double t[3][3]; |
| 743 |
|
| 744 |
halfpi = 2.0*atan(1.0); |
| 745 |
pi = 2.0*halfpi; |
| 746 |
iaxis = 0; |
| 747 |
for (i=6; i >= 2; i--) |
| 748 |
{ |
| 749 |
theta = 2.0*pi/(float)i; |
| 750 |
for (j=0; j < 3; j++) |
| 751 |
{ |
| 752 |
for (k=0; k < 3; k++) |
| 753 |
{ |
| 754 |
for (l=0; l < 3; l++) |
| 755 |
{ |
| 756 |
t[k][l] = 0.0; |
| 757 |
} |
| 758 |
} |
| 759 |
i1 = j+1; |
| 760 |
i2 = j+2; |
| 761 |
if (i1 > 2) i1 -= 3; |
| 762 |
if (i2 > 2) i2 -= 3; |
| 763 |
t[i1][i1] = cos(theta); |
| 764 |
t[i2][i2] = cos(theta); |
| 765 |
t[i1][i2] = sin(theta); |
| 766 |
t[i2][i1] = -sin(theta); |
| 767 |
t[j][j] = 1.0; |
| 768 |
tform(t,a,b,natom); |
| 769 |
iz = tcomp(a,b); |
| 770 |
if (iz == TRUE) |
| 771 |
{ |
| 772 |
*maxcn = i; |
| 773 |
for (k=0; k < 3; k++) |
| 774 |
{ |
| 775 |
p1[k] = 0.0; |
| 776 |
p2[k] = 0.0; |
| 777 |
if (k == j) p2[k] = 1.0; |
| 778 |
} |
| 779 |
orients(b,a,p1,p2,3,natom); |
| 780 |
return TRUE; |
| 781 |
} |
| 782 |
} |
| 783 |
} |
| 784 |
*maxcn = 1; |
| 785 |
return FALSE; |
| 786 |
} |
| 787 |
int findsn(double **a, double **b, double **c, int *maxcn) |
| 788 |
{ |
| 789 |
int iz; |
| 790 |
double pi,theta; |
| 791 |
pi = 4.0*atan(1.0); |
| 792 |
theta = pi/(*maxcn); |
| 793 |
rotates(a,b,3,theta); |
| 794 |
reflect(b,c,3); |
| 795 |
iz = tcomp(a,c); |
| 796 |
return iz; |
| 797 |
} |
| 798 |
int findc2(double **a, double **b, double **aset, int *nset) |
| 799 |
{ |
| 800 |
int i,j,k, j1, numset, iset, jset, iz; |
| 801 |
double halfpi, pi, x,y,theta, xi,yi; |
| 802 |
double proi, ani,disi; |
| 803 |
|
| 804 |
halfpi = 2.0*atan(1.0); |
| 805 |
pi = 2.0*halfpi; |
| 806 |
|
| 807 |
circset(a,aset,3,nset,&numset); |
| 808 |
for (i=1; i <= numset; i++) |
| 809 |
{ |
| 810 |
for (j=1; j <= natom; j++) |
| 811 |
{ |
| 812 |
if(nset[j] == i) // atom belongs to set |
| 813 |
{ |
| 814 |
if ( fabs(aset[j][1]) < TOLER ) |
| 815 |
{ |
| 816 |
j1 = j+1; |
| 817 |
for (k=j1; k <= natom; k++) |
| 818 |
{ |
| 819 |
if (nset[k] == i) |
| 820 |
{ |
| 821 |
x = (a[j][0]+a[k][0])*0.5; |
| 822 |
y = (a[j][1]+a[k][1])*0.5; |
| 823 |
theta = halfpi; |
| 824 |
if (fabs(y) > TOLER) theta = -atan(x/y); |
| 825 |
rotates(a,b,3,theta); |
| 826 |
rotates(b,aset,2,pi); |
| 827 |
iz = tcomp(b,aset); |
| 828 |
if (iz == TRUE) |
| 829 |
{ |
| 830 |
tcopy(b,a,natom); |
| 831 |
return TRUE; |
| 832 |
} |
| 833 |
} |
| 834 |
} |
| 835 |
} |
| 836 |
} |
| 837 |
} |
| 838 |
} |
| 839 |
// |
| 840 |
circset(a,aset,3,nset,&numset); |
| 841 |
iset = 1; |
| 842 |
for (i=1; i <= natom; i++) |
| 843 |
if (nset[i] == iset) goto L_10; |
| 844 |
return FALSE; |
| 845 |
// |
| 846 |
L_10: |
| 847 |
proi = aset[i][1]; |
| 848 |
ani = aset[i][0]; |
| 849 |
disi = aset[i][2]; |
| 850 |
xi = a[i][0]; |
| 851 |
yi = a[i][1]; |
| 852 |
j1 = iset+1; |
| 853 |
for (j=j1; j <= numset; j++) |
| 854 |
{ |
| 855 |
for (k=1; k <= natom; k++) |
| 856 |
{ |
| 857 |
if (nset[k] == j) |
| 858 |
{ |
| 859 |
jset = j; |
| 860 |
goto L_20; |
| 861 |
} |
| 862 |
} |
| 863 |
return FALSE; |
| 864 |
L_20: |
| 865 |
if ( fabs(proi+aset[k][1]) > TOLER || fabs(ani-aset[k][0]) > TOL2 || |
| 866 |
fabs(disi-aset[k][2]) > TOLER ) |
| 867 |
goto L_30; |
| 868 |
|
| 869 |
for (k=1; k <= natom; k++) |
| 870 |
{ |
| 871 |
if (nset[k] == jset) |
| 872 |
{ |
| 873 |
x = (xi+a[k][0])*0.5; |
| 874 |
y = (yi+a[k][1])*0.5; |
| 875 |
theta = halfpi; |
| 876 |
if (fabs(y) > TOL2) theta = -atan(x/y); |
| 877 |
rotates(a,b,3,theta); |
| 878 |
rotates(b,aset,2,pi); |
| 879 |
iz = tcomp(b,aset); |
| 880 |
if (iz == TRUE) |
| 881 |
{ |
| 882 |
tcopy(b,a,natom); |
| 883 |
return TRUE; |
| 884 |
} |
| 885 |
} |
| 886 |
} |
| 887 |
L_30: |
| 888 |
continue; |
| 889 |
} |
| 890 |
return FALSE; |
| 891 |
} |
| 892 |
int findv(double **a, double **b, double **aset, int *nset) |
| 893 |
{ |
| 894 |
int i,j, j1, numset, iset, jset, iz; |
| 895 |
double halfpi, pi, x,y,theta; |
| 896 |
|
| 897 |
halfpi = 2.0*atan(1.0); |
| 898 |
pi = 2.0*halfpi; |
| 899 |
|
| 900 |
circset(a,aset,3,nset,&numset); |
| 901 |
iset = 1; |
| 902 |
for (i=1; i < natom; i++) |
| 903 |
{ |
| 904 |
if (nset[i] == iset) |
| 905 |
{ |
| 906 |
jset = i; |
| 907 |
goto L_5; |
| 908 |
} |
| 909 |
} |
| 910 |
return FALSE; |
| 911 |
L_5: |
| 912 |
j1 = jset + 1; |
| 913 |
for (j=j1; j <= natom; j++) |
| 914 |
{ |
| 915 |
if (nset[j] == iset) |
| 916 |
{ |
| 917 |
x = (a[i][0]+a[j][0])*0.5; |
| 918 |
y = (a[i][1]+a[j][1])*0.5; |
| 919 |
if (fabs(x) <= TOLER && fabs(y) <= TOLER) |
| 920 |
{ |
| 921 |
x = 0.5*a[j][0]; |
| 922 |
y = 0.5*a[j][1]; |
| 923 |
} |
| 924 |
theta = halfpi; |
| 925 |
if (fabs(y) > TOLER) theta = -atan(x/y); |
| 926 |
rotates(a,b,3,theta); |
| 927 |
reflect(b,aset,1); |
| 928 |
iz = tcomp(b,aset); |
| 929 |
if (iz == TRUE) |
| 930 |
{ |
| 931 |
tcopy(b,a,natom); |
| 932 |
return TRUE; |
| 933 |
} |
| 934 |
} |
| 935 |
} |
| 936 |
return FALSE; |
| 937 |
} |
| 938 |
/* ======================================================= */ |
| 939 |
// compute moments of inertia |
| 940 |
void secmom(double **xyz, double *xi,double *yi, double *zi, int *ishape,double t[3][3]) |
| 941 |
{ |
| 942 |
int i,j,k,ii; |
| 943 |
float b[3]; |
| 944 |
double a[6], eigval[3], eigvec[3][3], atemp; |
| 945 |
double summ, wti, xx,yy,zz; |
| 946 |
|
| 947 |
for (i=0; i < 3; i++) |
| 948 |
{ |
| 949 |
b[i] = 0.0; |
| 950 |
eigval[i] = 0.0; |
| 951 |
for (j=0; j < 3; j++) |
| 952 |
eigvec[i][j] = 0.0; |
| 953 |
} |
| 954 |
for (i=0; i < 6; i++) |
| 955 |
a[i] = 0.0; |
| 956 |
|
| 957 |
summ = 0.0; |
| 958 |
for (i=1; i <= natom; i++) |
| 959 |
{ |
| 960 |
wti = atom.atomwt[i]; |
| 961 |
for (j=0; j < 3; j++) |
| 962 |
b[j] += xyz[i][j]*wti; |
| 963 |
summ += wti; |
| 964 |
} |
| 965 |
for (i=0; i < 3; i++) |
| 966 |
b[i] /= summ; |
| 967 |
for (i=1; i <= natom; i++) |
| 968 |
{ |
| 969 |
for (j=0; j < 3; j++) |
| 970 |
xyz[i][j] -= b[j]; |
| 971 |
} |
| 972 |
// compute moments of inertia , sort and place smallest moment on x axis |
| 973 |
for (i=1; i <= natom; i++) |
| 974 |
{ |
| 975 |
wti = atom.atomwt[i]; |
| 976 |
*xi = xyz[i][0]; |
| 977 |
*yi = xyz[i][1]; |
| 978 |
*zi = xyz[i][2]; |
| 979 |
xx = *xi* *xi; |
| 980 |
yy = *yi* *yi; |
| 981 |
zz = *zi* *zi; |
| 982 |
a[0] += wti*(yy+zz); |
| 983 |
a[2] += wti*(zz+xx); |
| 984 |
a[5] += wti*(xx+yy); |
| 985 |
a[1] -= wti**xi**yi; |
| 986 |
a[3] -= wti**xi**zi; |
| 987 |
a[4] -= wti**yi**zi; |
| 988 |
} |
| 989 |
summ = fabs(a[1]) + fabs(a[3]) + fabs(a[4]); |
| 990 |
if (summ <= 1.0e-15) |
| 991 |
{ |
| 992 |
eigval[0] = a[0]; |
| 993 |
eigval[1] = a[2]; |
| 994 |
eigval[2] = a[5]; |
| 995 |
eigvec[0][0] = 1.0; |
| 996 |
eigvec[1][1] = 1.0; |
| 997 |
eigvec[2][2] = 1.0; |
| 998 |
} else |
| 999 |
{ |
| 1000 |
hqrii(a,3,3,eigval,eigvec); |
| 1001 |
} |
| 1002 |
|
| 1003 |
for (i=0; i < 2; i++) |
| 1004 |
{ |
| 1005 |
ii = i + 1; |
| 1006 |
for (j=ii; j < 3; j++) |
| 1007 |
{ |
| 1008 |
if (eigval[i] > eigval[j]) |
| 1009 |
{ |
| 1010 |
atemp = eigval[i]; |
| 1011 |
eigval[i] = eigval[j]; |
| 1012 |
eigval[j] = atemp; |
| 1013 |
for (k=0; k < 3; k++) |
| 1014 |
{ |
| 1015 |
atemp = eigvec[k][i]; |
| 1016 |
eigvec[k][i] = eigvec[k][j]; |
| 1017 |
eigvec[k][j] = atemp; |
| 1018 |
} |
| 1019 |
} |
| 1020 |
} |
| 1021 |
} |
| 1022 |
*xi = eigval[0] / 6.0228; |
| 1023 |
*yi = eigval[1] / 6.0228; |
| 1024 |
*zi = eigval[2] / 6.0228; |
| 1025 |
// now rotate the molecule |
| 1026 |
for (i=1; i <= natom; i++) |
| 1027 |
{ |
| 1028 |
for (j=0; j < 3; j++) |
| 1029 |
{ |
| 1030 |
b[j] = 0.0; |
| 1031 |
for (k=0; k < 3; k++) |
| 1032 |
b[j] += eigvec[j][k]*xyz[i][k]; |
| 1033 |
} |
| 1034 |
xyz[i][0] = b[0]; |
| 1035 |
xyz[i][1] = b[1]; |
| 1036 |
xyz[i][2] = b[2]; |
| 1037 |
} |
| 1038 |
// assign shape |
| 1039 |
*ishape = 0; |
| 1040 |
if (*xi <= 0.01) |
| 1041 |
*ishape = 1; |
| 1042 |
else |
| 1043 |
{ |
| 1044 |
if (fabs(*xi+*yi-*zi) <= 0.01) |
| 1045 |
*ishape = 6; |
| 1046 |
else |
| 1047 |
{ |
| 1048 |
if (fabs( (*zi/(*yi))-(*yi/(*xi))) <= 0.01) |
| 1049 |
*ishape = 5; |
| 1050 |
else |
| 1051 |
{ |
| 1052 |
atemp = (2.0/(*yi) - 1.0/(*zi) - 1.0/(*xi))/ (1.0/(*xi) - 1.0/(*zi)); |
| 1053 |
if (atemp <= -0.999) |
| 1054 |
*ishape = 3; |
| 1055 |
else |
| 1056 |
{ |
| 1057 |
if (atemp >= 0.999) |
| 1058 |
*ishape = 4; |
| 1059 |
else |
| 1060 |
*ishape = 2; |
| 1061 |
} |
| 1062 |
} |
| 1063 |
} |
| 1064 |
} |
| 1065 |
} |
| 1066 |
// assign point group, symmetry number and chirality |
| 1067 |
void vibsym(int mode) |
| 1068 |
{ |
| 1069 |
int i; |
| 1070 |
double **a; |
| 1071 |
char ngrp[4]; |
| 1072 |
|
| 1073 |
a = dmatrix (0, natom+1, 0,3); |
| 1074 |
for (i=1; i <= natom; i++) |
| 1075 |
{ |
| 1076 |
a[i][0] = atom.x[i]; |
| 1077 |
a[i][1] = atom.y[i]; |
| 1078 |
a[i][2] = atom.z[i]; |
| 1079 |
} |
| 1080 |
|
| 1081 |
if (mode == 0) |
| 1082 |
{ |
| 1083 |
strcpy(ngrp,""); |
| 1084 |
ptgrp(mode, ngrp,a); |
| 1085 |
} |
| 1086 |
free_dmatrix(a ,0, natom+1, 0,3); |
| 1087 |
} |
| 1088 |
|
| 1089 |
void ptgrp(int mode,char *ngrp,double **a) |
| 1090 |
{ |
| 1091 |
int i, j, k, ishape, numsym, iaxis, maxcn, xxx; |
| 1092 |
int iz, iz1, iz2, iopt; |
| 1093 |
int chiral, *nset; |
| 1094 |
double **b, **c; |
| 1095 |
double t[3][3], p1[3], p2[3]; |
| 1096 |
double xi, yi, zi; |
| 1097 |
float halfpi, pi; |
| 1098 |
char string[256]; |
| 1099 |
|
| 1100 |
numsym = 0; |
| 1101 |
halfpi = 2.0*atan(1.0); |
| 1102 |
pi = 2.0*halfpi; |
| 1103 |
chiral = FALSE; |
| 1104 |
strcpy(ngrp,"Cs"); |
| 1105 |
|
| 1106 |
|
| 1107 |
b = dmatrix (0, natom+1, 0,3); |
| 1108 |
c = dmatrix (0, natom+1, 0,3); |
| 1109 |
nset = ivector(0,natom+1); |
| 1110 |
|
| 1111 |
secmom(a, &xi, &yi, &zi, &ishape, t); |
| 1112 |
|
| 1113 |
if (ishape == 1) // linear |
| 1114 |
{ |
| 1115 |
p1[0] = 0.0; |
| 1116 |
p1[1] = 0.0; |
| 1117 |
p1[2] = 0.0; |
| 1118 |
p2[0] = 1.0; |
| 1119 |
p2[1] = 0.0; |
| 1120 |
p2[1] = 0.0; |
| 1121 |
orients(a,b,p1,p2,3,natom); |
| 1122 |
tcopy(b,a,natom); |
| 1123 |
iaxis = findh(a,b,3); |
| 1124 |
if (iaxis == TRUE) |
| 1125 |
{ |
| 1126 |
strcpy(ngrp,"D*h"); |
| 1127 |
numsym = 2; |
| 1128 |
} else |
| 1129 |
{ |
| 1130 |
strcpy(ngrp,"C*v"); |
| 1131 |
numsym = 1; |
| 1132 |
} |
| 1133 |
} else if (ishape == 2) // asymmetric |
| 1134 |
/*------------------------* |
| 1135 |
ASyMMETRiC TOP MOLECULES |
| 1136 |
*------------------------* |
| 1137 |
|
| 1138 |
These molecules can have no axes of order greater than 2. Thus |
| 1139 |
the possible point groups are: D2H, D2, C2V, C2H, C2, Ci, CS, |
| 1140 |
and C1. FiNDCN routine finds the principal axes and align the |
| 1141 |
principal axes with the z-axis. */ |
| 1142 |
{ |
| 1143 |
iaxis = findcn(a,b,p1,p2,&maxcn); |
| 1144 |
if (maxcn == 1) |
| 1145 |
{ |
| 1146 |
j = findh(a,b,3); |
| 1147 |
if (j == TRUE) |
| 1148 |
{ |
| 1149 |
strcpy(ngrp,"Cs"); |
| 1150 |
numsym = 1; |
| 1151 |
}else |
| 1152 |
{ |
| 1153 |
i = findh(a,b,2); |
| 1154 |
if (i == TRUE) |
| 1155 |
{ |
| 1156 |
strcpy(ngrp,"Cs"); |
| 1157 |
numsym = 1; |
| 1158 |
} else |
| 1159 |
{ |
| 1160 |
k = findh(a,b,1); |
| 1161 |
if (k == TRUE) |
| 1162 |
{ |
| 1163 |
strcpy(ngrp,"Cs"); |
| 1164 |
numsym = 1; |
| 1165 |
} else |
| 1166 |
{ |
| 1167 |
xxx = findi(a,b); |
| 1168 |
if (xxx == TRUE) |
| 1169 |
{ |
| 1170 |
strcpy(ngrp,"Ci"); |
| 1171 |
numsym = 1; |
| 1172 |
} else |
| 1173 |
{ |
| 1174 |
strcpy(ngrp,"C1"); |
| 1175 |
numsym = 1; |
| 1176 |
chiral = TRUE; |
| 1177 |
} |
| 1178 |
} |
| 1179 |
} |
| 1180 |
} |
| 1181 |
} else // maxcn == 2 |
| 1182 |
{ |
| 1183 |
iaxis = findc2(a,b,c,nset); |
| 1184 |
if (iaxis == FALSE) |
| 1185 |
{ |
| 1186 |
i = findv(a,b,c,nset); |
| 1187 |
if (i == TRUE) |
| 1188 |
{ |
| 1189 |
strcpy(ngrp,"C2v"); |
| 1190 |
numsym = 2; |
| 1191 |
} else |
| 1192 |
{ |
| 1193 |
j = findh(a,b,3); |
| 1194 |
if (j == TRUE) |
| 1195 |
{ |
| 1196 |
strcpy(ngrp,"C2h"); |
| 1197 |
numsym = 2; |
| 1198 |
} else |
| 1199 |
{ |
| 1200 |
strcpy(ngrp,"C2"); |
| 1201 |
numsym = 2; |
| 1202 |
chiral = TRUE; |
| 1203 |
} |
| 1204 |
} |
| 1205 |
} else |
| 1206 |
{ |
| 1207 |
i = findh(a,b,3); |
| 1208 |
if (i == TRUE) |
| 1209 |
{ |
| 1210 |
strcpy(ngrp,"D2h"); |
| 1211 |
numsym = 4; |
| 1212 |
} else |
| 1213 |
{ |
| 1214 |
j = findv(a,b,c,nset); |
| 1215 |
if (j == TRUE) |
| 1216 |
{ |
| 1217 |
strcpy(ngrp,"D2d"); |
| 1218 |
numsym = 4; |
| 1219 |
} else |
| 1220 |
{ |
| 1221 |
strcpy(ngrp,"D2"); |
| 1222 |
numsym = 4; |
| 1223 |
chiral = TRUE; |
| 1224 |
} |
| 1225 |
} |
| 1226 |
} |
| 1227 |
} |
| 1228 |
} else if (ishape == 3 || ishape == 4) //prolate |
| 1229 |
{ |
| 1230 |
/*-----------------------* |
| 1231 |
SyMMETRiC TOP MOLECULES |
| 1232 |
*-----------------------* |
| 1233 |
|
| 1234 |
These molecules can belong to any axial point group, thus only |
| 1235 |
the cubic point groups (T, Td, Th, O, Oh, i, ih) are impossible. |
| 1236 |
However, except in rare cases the unique axis is a rotation axis |
| 1237 |
of order 3 or greater. */ |
| 1238 |
|
| 1239 |
iz = findcn(a,b,p1,p2,&maxcn); |
| 1240 |
iz = findsn(a,b,c,&maxcn); |
| 1241 |
if (iz == TRUE) |
| 1242 |
{ |
| 1243 |
iz1 = findv(a,b,c,nset); |
| 1244 |
if (iz1 == FALSE) |
| 1245 |
{ |
| 1246 |
sprintf(ngrp,"S%1d",2*maxcn); |
| 1247 |
numsym = maxcn; |
| 1248 |
goto L_2000; |
| 1249 |
} |
| 1250 |
} |
| 1251 |
iz = findc2(a,b,c,nset); |
| 1252 |
if (iz == FALSE) |
| 1253 |
{ |
| 1254 |
iz1 = findv(a,b,c,nset); |
| 1255 |
if (iz1 == TRUE) |
| 1256 |
{ |
| 1257 |
sprintf(ngrp,"C%1dv",maxcn); |
| 1258 |
numsym = maxcn; |
| 1259 |
} else |
| 1260 |
{ |
| 1261 |
iz2 = findh(a,b,3); |
| 1262 |
if (iz2 == TRUE) |
| 1263 |
{ |
| 1264 |
sprintf(ngrp,"C%1dh",maxcn); |
| 1265 |
numsym = maxcn; |
| 1266 |
} else |
| 1267 |
{ |
| 1268 |
sprintf(ngrp,"C%1d",maxcn); |
| 1269 |
numsym = maxcn; |
| 1270 |
chiral = TRUE; |
| 1271 |
} |
| 1272 |
} |
| 1273 |
}else |
| 1274 |
{ |
| 1275 |
iz1 = findh(a,b,3); |
| 1276 |
if (iz1 == TRUE) |
| 1277 |
{ |
| 1278 |
sprintf(ngrp,"D%1dh",maxcn); |
| 1279 |
numsym = 2*maxcn; |
| 1280 |
} else |
| 1281 |
{ |
| 1282 |
iz2 = findv(a,b,c,nset); |
| 1283 |
if (iz2 == TRUE) |
| 1284 |
{ |
| 1285 |
sprintf(ngrp,"D%1dd",maxcn); |
| 1286 |
numsym = 2*maxcn; |
| 1287 |
} else |
| 1288 |
{ |
| 1289 |
sprintf(ngrp,"D%1d",maxcn); |
| 1290 |
numsym = 2*maxcn; |
| 1291 |
chiral = TRUE; |
| 1292 |
} |
| 1293 |
} |
| 1294 |
} |
| 1295 |
// SPHERiCAL TOP MOLECULES |
| 1296 |
} else if (ishape == 5) // spherical |
| 1297 |
{ |
| 1298 |
sphere(a,b,c,nset,&maxcn); |
| 1299 |
iopt = maxcn - 2; |
| 1300 |
if (iopt == 1) |
| 1301 |
{ |
| 1302 |
iz = findi(a,b); |
| 1303 |
if (iz == TRUE) |
| 1304 |
{ |
| 1305 |
sprintf(ngrp,"Th"); |
| 1306 |
numsym = 12; |
| 1307 |
} else |
| 1308 |
{ |
| 1309 |
iz1 = findv(a,b,c,nset); |
| 1310 |
if (iz1 == TRUE) |
| 1311 |
{ |
| 1312 |
sprintf(ngrp,"Td"); |
| 1313 |
numsym = 12; |
| 1314 |
} else |
| 1315 |
{ |
| 1316 |
sprintf(ngrp,"T"); |
| 1317 |
numsym = 12; |
| 1318 |
chiral = TRUE; |
| 1319 |
} |
| 1320 |
} |
| 1321 |
} else if (iopt == 2) |
| 1322 |
{ |
| 1323 |
iz = findi(a,b); |
| 1324 |
if (iz == TRUE) |
| 1325 |
{ |
| 1326 |
sprintf(ngrp,"Oh"); |
| 1327 |
numsym = 24; |
| 1328 |
} else |
| 1329 |
{ |
| 1330 |
sprintf(ngrp,"O"); |
| 1331 |
numsym = 24; |
| 1332 |
chiral = TRUE; |
| 1333 |
} |
| 1334 |
} else if (iopt == 3) |
| 1335 |
{ |
| 1336 |
iz = findi(a,b); |
| 1337 |
if (iz == TRUE) |
| 1338 |
{ |
| 1339 |
sprintf(ngrp,"Ih"); |
| 1340 |
numsym = 60; |
| 1341 |
} else |
| 1342 |
{ |
| 1343 |
sprintf(ngrp,"I"); |
| 1344 |
numsym = 60; |
| 1345 |
chiral = TRUE; |
| 1346 |
} |
| 1347 |
} |
| 1348 |
} else if (ishape == 6) // planar |
| 1349 |
{ |
| 1350 |
iz = findh(a,b,1); |
| 1351 |
iz1 = findh(a,b,2); |
| 1352 |
iz2 = findcn(a,b,p1,p2,&maxcn); |
| 1353 |
if (maxcn == 2) |
| 1354 |
{ |
| 1355 |
if (iz == TRUE && iz1 == TRUE) |
| 1356 |
{ |
| 1357 |
sprintf(ngrp,"D%1dh",maxcn); |
| 1358 |
numsym = maxcn*2; |
| 1359 |
} else |
| 1360 |
{ |
| 1361 |
if (iz == TRUE || iz1 == TRUE) |
| 1362 |
{ |
| 1363 |
sprintf(ngrp,"C%1dv",maxcn); |
| 1364 |
numsym = maxcn; |
| 1365 |
}else |
| 1366 |
{ |
| 1367 |
i = findh(a,b,3); |
| 1368 |
if (i == TRUE) |
| 1369 |
{ |
| 1370 |
sprintf(ngrp,"C%1dh",maxcn); |
| 1371 |
numsym = maxcn; |
| 1372 |
} else |
| 1373 |
{ |
| 1374 |
sprintf(ngrp,"C%1d",maxcn); |
| 1375 |
numsym = maxcn; |
| 1376 |
chiral = TRUE; |
| 1377 |
} |
| 1378 |
} |
| 1379 |
} |
| 1380 |
} else |
| 1381 |
{ |
| 1382 |
iz = findsn(a,b,c,&maxcn); |
| 1383 |
if (iz == TRUE) |
| 1384 |
{ |
| 1385 |
iz1 = findv(a,b,c,nset); |
| 1386 |
if (iz1 == FALSE) |
| 1387 |
{ |
| 1388 |
sprintf(ngrp,"S%1d",2*maxcn); |
| 1389 |
numsym = maxcn; |
| 1390 |
goto L_2000; |
| 1391 |
} |
| 1392 |
} |
| 1393 |
iz = findc2(a,b,c,nset); |
| 1394 |
if (iz == FALSE) |
| 1395 |
{ |
| 1396 |
iz1 = findv(a,b,c,nset); |
| 1397 |
if (iz1 == TRUE) |
| 1398 |
{ |
| 1399 |
sprintf(ngrp,"C%1dv",maxcn); |
| 1400 |
numsym = maxcn; |
| 1401 |
} else |
| 1402 |
{ |
| 1403 |
iz2 = findh(a,b,3); |
| 1404 |
if (iz2 == TRUE) |
| 1405 |
{ |
| 1406 |
sprintf(ngrp,"C%1dh",maxcn); |
| 1407 |
numsym = maxcn; |
| 1408 |
} else |
| 1409 |
{ |
| 1410 |
sprintf(ngrp,"C%1d",maxcn); |
| 1411 |
numsym = maxcn; |
| 1412 |
chiral = TRUE; |
| 1413 |
} |
| 1414 |
} |
| 1415 |
} else |
| 1416 |
{ |
| 1417 |
iz1 = findh(a,b,3); |
| 1418 |
if (iz1 == TRUE) |
| 1419 |
{ |
| 1420 |
sprintf(ngrp,"D%1dh",maxcn); |
| 1421 |
numsym = 2*maxcn; |
| 1422 |
} else |
| 1423 |
{ |
| 1424 |
iz2 = findv(a,b,c,nset); |
| 1425 |
if (iz2 == TRUE) |
| 1426 |
{ |
| 1427 |
sprintf(ngrp,"D%1dd",maxcn); |
| 1428 |
numsym = 2*maxcn; |
| 1429 |
} else |
| 1430 |
{ |
| 1431 |
sprintf(ngrp,"D%1d",maxcn); |
| 1432 |
numsym = 2*maxcn; |
| 1433 |
chiral = TRUE; |
| 1434 |
} |
| 1435 |
} |
| 1436 |
} |
| 1437 |
} |
| 1438 |
} |
| 1439 |
L_2000: |
| 1440 |
free_dmatrix(b ,0, natom+1, 0,3); |
| 1441 |
free_dmatrix(c ,0, natom+1, 0,3); |
| 1442 |
free_ivector(nset, 0, natom+1); |
| 1443 |
|
| 1444 |
if (mode == 2) |
| 1445 |
{ |
| 1446 |
symmetry.chiral = chiral; |
| 1447 |
symmetry.numsym = numsym; |
| 1448 |
symmetry.ishape = ishape; |
| 1449 |
symmetry.xi = xi; |
| 1450 |
symmetry.yi = yi; |
| 1451 |
symmetry.zi = zi; |
| 1452 |
} |
| 1453 |
|
| 1454 |
if (mode == 0) |
| 1455 |
{ |
| 1456 |
sprintf(string,"Symmetry group: %s\nMoments of Inertia\nIx: %7.2f\nIy: %7.2f\nIz: %7.2f",ngrp,xi,yi,zi); |
| 1457 |
message_alert(string,"Symmetry Information"); |
| 1458 |
} |
| 1459 |
} |
| 1460 |
/* ========================================== */ |
| 1461 |
void hqrii(double *a,int n,int m,double *e,double v[3][3]) |
| 1462 |
{ |
| 1463 |
#define V(I_,J_) (*(v+(I_)*(n)+(J_))) |
| 1464 |
int _do0, _do1, _do10, _do11, _do12, _do13, _do14, _do15, |
| 1465 |
_do16, _do17, _do18, _do19, _do2, _do20, _do21, _do22, _do3, |
| 1466 |
_do4, _do5, _do6, _do7, _do8, _do9, i, i_, ig, ii, im1, iord, |
| 1467 |
ip1, irank, itere, itere_, j, j_, jrank, k, k_, kk, kk_, kp1, |
| 1468 |
kpiv, krank, l, ll, nm1, nm2=0; |
| 1469 |
double c, del, ee, eps, eps1, eps2, eps3, ff, fn, gersch, h, r, |
| 1470 |
ra, rn, s, seps, sinv, sorter, sum, summ, t, u, vn, w[5][5], |
| 1471 |
ww, z, zero; |
| 1472 |
|
| 1473 |
double *const A = &a[0] - 1; |
| 1474 |
double *const E = &e[0] - 1; |
| 1475 |
|
| 1476 |
/************************************************************* |
| 1477 |
* |
| 1478 |
* HQRII IS A DIAGONALISATION ROUTINE, WRITTEN BY YOSHITAKA BEPPU OF |
| 1479 |
* NAGOYA UNIVERSITY, JAPAN. |
| 1480 |
* FOR DETAILS SEE 'COMPUTERS & CHEMISTRY' VOL.6 1982. PAGE 000. |
| 1481 |
* |
| 1482 |
* ON INPUT A = MATRIX TO BE DIAGONALISED (PACKED CANONICAL) |
| 1483 |
* N = SIZE OF MATRIX TO BE DIAGONALISED. |
| 1484 |
* M = NUMBER OF EIGENVECTORS NEEDED. |
| 1485 |
* E = ARRAY OF SIZE AT LEAST N |
| 1486 |
* V = ARRAY OF SIZE AT LEAST NMAX*M |
| 1487 |
* |
| 1488 |
* ON OUTPUT E = EIGENVALUES |
| 1489 |
* V = EIGENVECTORS IN ARRAY OF SIZE NMAX*M |
| 1490 |
* |
| 1491 |
************************************************************************ */ |
| 1492 |
|
| 1493 |
/* EPS3 AND EPS ARE MACHINE-PRECISION DEPENDENT |
| 1494 |
* */ |
| 1495 |
eps3 = 1.e-30; |
| 1496 |
zero = 0.e0; |
| 1497 |
ll = (n*(n + 1))/2 + 1; |
| 1498 |
eps = 1.e-8; |
| 1499 |
iord = -1; |
| 1500 |
nm1 = n - 1; |
| 1501 |
if( n == 2 ) |
| 1502 |
goto L_90; |
| 1503 |
nm2 = n - 2; |
| 1504 |
krank = 0; |
| 1505 |
/* HOUSEHOLDER TRANSFORMATION */ |
| 1506 |
for( k = 1, _do0 = nm2; k <= _do0; k++ ) |
| 1507 |
{ |
| 1508 |
k_ = k - 1; |
| 1509 |
kp1 = k + 1; |
| 1510 |
krank += k; |
| 1511 |
w[k_][1] = A[krank]; |
| 1512 |
sum = 0.; |
| 1513 |
jrank = krank; |
| 1514 |
for( j = kp1, _do1 = n; j <= _do1; j++ ) |
| 1515 |
{ |
| 1516 |
j_ = j - 1; |
| 1517 |
w[j_][1] = A[jrank + k]; |
| 1518 |
jrank += j; |
| 1519 |
sum += (w[j_][1])*(w[j_][1]); |
| 1520 |
} |
| 1521 |
s = SIGN( sqrt( sum ), w[kp1 - 1][1] ); |
| 1522 |
w[k_][0] = -s; |
| 1523 |
w[kp1 - 1][1] += s; |
| 1524 |
A[k + krank] = w[kp1 - 1][1]; |
| 1525 |
h = w[kp1 - 1][1]*s; |
| 1526 |
if( fabs( h ) < eps3 ) |
| 1527 |
goto L_80; |
| 1528 |
summ = 0.e0; |
| 1529 |
irank = krank; |
| 1530 |
for( i = kp1, _do2 = n; i <= _do2; i++ ) |
| 1531 |
{ |
| 1532 |
i_ = i - 1; |
| 1533 |
sum = 0.e0; |
| 1534 |
for( j = kp1, _do3 = i; j <= _do3; j++ ) |
| 1535 |
{ |
| 1536 |
j_ = j - 1; |
| 1537 |
sum += A[j + irank]*w[j_][1]; |
| 1538 |
} |
| 1539 |
if( i >= n ) |
| 1540 |
goto L_40; |
| 1541 |
ip1 = i + 1; |
| 1542 |
jrank = i*(i + 3)/2; |
| 1543 |
for( j = ip1, _do4 = n; j <= _do4; j++ ) |
| 1544 |
{ |
| 1545 |
j_ = j - 1; |
| 1546 |
sum += A[jrank]*w[j_][1]; |
| 1547 |
jrank += j; |
| 1548 |
} |
| 1549 |
L_40: |
| 1550 |
w[i_][0] = sum/h; |
| 1551 |
irank += i; |
| 1552 |
summ += w[i_][0]*w[i_][1]; |
| 1553 |
} |
| 1554 |
u = summ*0.5e0/h; |
| 1555 |
jrank = krank; |
| 1556 |
for( j = kp1, _do5 = n; j <= _do5; j++ ) |
| 1557 |
{ |
| 1558 |
j_ = j - 1; |
| 1559 |
w[j_][0] = w[j_][1]*u - w[j_][0]; |
| 1560 |
for( i = kp1, _do6 = j; i <= _do6; i++ ) |
| 1561 |
{ |
| 1562 |
i_ = i - 1; |
| 1563 |
A[i + jrank] += w[i_][0]*w[j_][1] + w[j_][0]*w[i_][1]; |
| 1564 |
} |
| 1565 |
jrank += j; |
| 1566 |
} |
| 1567 |
L_80: |
| 1568 |
A[krank] = h; |
| 1569 |
} |
| 1570 |
L_90: |
| 1571 |
w[nm1 - 1][1] = A[(nm1*(nm1 + 1))/2]; |
| 1572 |
w[n - 1][1] = A[(n*(n + 1))/2]; |
| 1573 |
w[nm1 - 1][0] = A[nm1 + (n*(n - 1))/2]; |
| 1574 |
w[n - 1][0] = 0.e0; |
| 1575 |
gersch = fabs( w[0][1] ) + fabs( w[0][0] ); |
| 1576 |
for( i = 1, _do7 = nm1; i <= _do7; i++ ) |
| 1577 |
{ |
| 1578 |
i_ = i - 1; |
| 1579 |
if ( fabs( w[i_ + 1][1] ) + fabs( w[i_][0] ) + fabs( w[i_ + 1][0] ) > gersch) |
| 1580 |
gersch = fabs( w[i_ + 1][1] ) + fabs( w[i_][0] ) + fabs( w[i_ + 1][0] ); |
| 1581 |
} |
| 1582 |
del = eps*gersch; |
| 1583 |
for( i = 1, _do8 = n; i <= _do8; i++ ) |
| 1584 |
{ |
| 1585 |
i_ = i - 1; |
| 1586 |
w[i_][2] = w[i_][0]; |
| 1587 |
E[i] = w[i_][1]; |
| 1588 |
v[m-1][i_] = E[i]; |
| 1589 |
} |
| 1590 |
if( fabs( del ) < eps3 ) |
| 1591 |
goto L_220; |
| 1592 |
/* QR-METHOD WITH ORIGIN SHIFT */ |
| 1593 |
k = n; |
| 1594 |
L_120: |
| 1595 |
l = k; |
| 1596 |
L_130: |
| 1597 |
if( fabs( w[l - 2][2] ) < del ) |
| 1598 |
goto L_140; |
| 1599 |
l -= 1; |
| 1600 |
if( l > 1 ) |
| 1601 |
goto L_130; |
| 1602 |
L_140: |
| 1603 |
if( l == k ) |
| 1604 |
goto L_170; |
| 1605 |
ww = (E[k - 1] + E[k])*0.5e0; |
| 1606 |
r = E[k] - ww; |
| 1607 |
z = SIGN( sqrt( (w[k - 2][2])*(w[k - 2][2]) + r*r ), r ) + ww; |
| 1608 |
ee = E[l] - z; |
| 1609 |
E[l] = ee; |
| 1610 |
ff = w[l - 1][2]; |
| 1611 |
r = sqrt( ee*ee + ff*ff ); |
| 1612 |
j = l; |
| 1613 |
goto L_160; |
| 1614 |
L_150: |
| 1615 |
r = sqrt( (E[j])*(E[j]) + (w[j - 1][2])*(w[j - 1][2]) ); |
| 1616 |
w[j - 2][2] = s*r; |
| 1617 |
ee = E[j]*c; |
| 1618 |
ff = w[j - 1][2]*c; |
| 1619 |
L_160: |
| 1620 |
c = E[j]/r; |
| 1621 |
s = w[j - 1][2]/r; |
| 1622 |
ww = E[j + 1] - z; |
| 1623 |
E[j] = (ff*c + ww*s)*s + ee + z; |
| 1624 |
E[j + 1] = c*ww - s*ff; |
| 1625 |
j += 1; |
| 1626 |
if( j < k ) |
| 1627 |
goto L_150; |
| 1628 |
w[k - 2][2] = E[k]*s; |
| 1629 |
E[k] = E[k]*c + z; |
| 1630 |
goto L_120; |
| 1631 |
L_170: |
| 1632 |
k -= 1; |
| 1633 |
if( k > 1 ) |
| 1634 |
goto L_120; |
| 1635 |
/* * * * * * * * * * * * * |
| 1636 |
* |
| 1637 |
* AT THIS POINT THE ARRAY 'E' CONTAINS THE UN-ORDERED EIGENVALUES |
| 1638 |
* |
| 1639 |
* * * * * * * * * * * * * |
| 1640 |
* STRAIGHT SELECTION SORT OF EIGENVALUES */ |
| 1641 |
sorter = 1.e0; |
| 1642 |
if( iord < 0 ) |
| 1643 |
sorter = -1.e0; |
| 1644 |
j = n; |
| 1645 |
L_180: |
| 1646 |
l = 1; |
| 1647 |
ii = 1; |
| 1648 |
ll = 1; |
| 1649 |
for( i = 2, _do9 = j; i <= _do9; i++ ) |
| 1650 |
{ |
| 1651 |
i_ = i - 1; |
| 1652 |
if( (E[i] - E[l])*sorter > 0.e0 ) |
| 1653 |
goto L_190; |
| 1654 |
l = i; |
| 1655 |
goto L_200; |
| 1656 |
L_190: |
| 1657 |
ii = i; |
| 1658 |
ll = l; |
| 1659 |
L_200: |
| 1660 |
; |
| 1661 |
} |
| 1662 |
if( ii == ll ) |
| 1663 |
goto L_210; |
| 1664 |
ww = E[ll]; |
| 1665 |
E[ll] = E[ii]; |
| 1666 |
E[ii] = ww; |
| 1667 |
L_210: |
| 1668 |
j = ii - 1; |
| 1669 |
if( j >= 2 ) |
| 1670 |
goto L_180; |
| 1671 |
L_220: |
| 1672 |
if( !m ) |
| 1673 |
return; |
| 1674 |
/*************** |
| 1675 |
* ORDERING OF EIGENVALUES COMPLETE. |
| 1676 |
*************** |
| 1677 |
* INVERSE-ITERATION FOR EIGENVECTORS */ |
| 1678 |
fn = (float)( n ); |
| 1679 |
eps1 = 1.e-5; |
| 1680 |
seps = sqrt( eps ); |
| 1681 |
eps2 = 0.05e0; |
| 1682 |
rn = 0.e0; |
| 1683 |
ra = eps*0.6180339887485e0; |
| 1684 |
/* 0.618... IS THE FIBONACCI NUMBER (-1+SQRT(5))/2. */ |
| 1685 |
ig = 1; |
| 1686 |
for( i = 1, _do10 = m; i <= _do10; i++ ) |
| 1687 |
{ |
| 1688 |
i_ = i - 1; |
| 1689 |
im1 = i - 1; |
| 1690 |
for( j = 1, _do11 = n; j <= _do11; j++ ) |
| 1691 |
{ |
| 1692 |
j_ = j - 1; |
| 1693 |
w[j_][2] = 0.e0; |
| 1694 |
w[j_][3] = w[j_][0]; |
| 1695 |
w[j_][4] = v[m-1][j_] - E[i]; |
| 1696 |
rn += ra; |
| 1697 |
if( rn >= eps ) |
| 1698 |
rn -= eps; |
| 1699 |
v[i_][j_] = rn; |
| 1700 |
} |
| 1701 |
for( j = 1, _do12 = nm1; j <= _do12; j++ ) |
| 1702 |
{ |
| 1703 |
j_ = j - 1; |
| 1704 |
if( fabs( w[j_][4] ) >= fabs( w[j_][0] ) ) |
| 1705 |
goto L_240; |
| 1706 |
w[j_][1] = -w[j_][4]/w[j_][0]; |
| 1707 |
w[j_][4] = w[j_][0]; |
| 1708 |
t = w[j_ + 1][4]; |
| 1709 |
w[j_ + 1][4] = w[j_][3]; |
| 1710 |
w[j_][3] = t; |
| 1711 |
w[j_][2] = w[j_ + 1][3]; |
| 1712 |
if( fabs( w[j_][2] ) < eps3 ) |
| 1713 |
w[j_][2] = del; |
| 1714 |
w[j_ + 1][3] = 0.e0; |
| 1715 |
goto L_250; |
| 1716 |
L_240: |
| 1717 |
if( fabs( w[j_][4] ) < eps3 ) |
| 1718 |
w[j_][4] = del; |
| 1719 |
w[j_][1] = -w[j_][0]/w[j_][4]; |
| 1720 |
L_250: |
| 1721 |
w[j_ + 1][3] += w[j_][2]*w[j_][1]; |
| 1722 |
w[j_ + 1][4] += w[j_][3]*w[j_][1]; |
| 1723 |
} |
| 1724 |
if( fabs( w[n - 1][4] ) < eps3 ) |
| 1725 |
w[n - 1][4] = del; |
| 1726 |
for( itere = 1; itere <= 5; itere++ ) |
| 1727 |
{ |
| 1728 |
itere_ = itere - 1; |
| 1729 |
if( itere == 1 ) |
| 1730 |
goto L_280; |
| 1731 |
for( j = 1, _do13 = nm1; j <= _do13; j++ ) |
| 1732 |
{ |
| 1733 |
j_ = j - 1; |
| 1734 |
if( fabs( w[j_][2] ) < eps3 ) |
| 1735 |
goto L_270; |
| 1736 |
t = v[i_][j_]; |
| 1737 |
v[i_][j_] = v[i_][j_+1]; |
| 1738 |
v[i_][j_+1] = t; |
| 1739 |
L_270: |
| 1740 |
v[i_][j_+1] += v[i_][j_]*w[j_][1]; |
| 1741 |
} |
| 1742 |
L_280: |
| 1743 |
v[i_][n-1] /= w[n - 1][4]; |
| 1744 |
v[i_][nm1 - 1] = (v[i_][nm1 - 1] - v[i_][n-1]*w[nm1 - 1][3])/w[nm1 - 1][4]; |
| 1745 |
vn = 1.0e-20; |
| 1746 |
if (fabs( v[i_][n-1] ) > vn) |
| 1747 |
vn = fabs( v[i_][n-1] ); |
| 1748 |
if ( fabs( v[i_][nm1 - 1] ) > vn) |
| 1749 |
vn = fabs( v[i_][nm1 - 1] ); |
| 1750 |
if( n == 2 ) |
| 1751 |
goto L_300; |
| 1752 |
k = nm2; |
| 1753 |
L_290: |
| 1754 |
v[i_][k - 1] = (v[i_][k - 1] - v[i_][k]*w[k - 1][3] - v[i_][k + 1]* |
| 1755 |
w[k - 1][2])/w[k - 1][4]; |
| 1756 |
if (vn < 1.0e-20) |
| 1757 |
vn = 1.0e-20; |
| 1758 |
if (fabs( v[i_][k - 1] ) > vn) |
| 1759 |
vn = fabs( v[i_][n-1] ); |
| 1760 |
k -= 1; |
| 1761 |
if( k >= 1 ) |
| 1762 |
goto L_290; |
| 1763 |
L_300: |
| 1764 |
s = eps1/vn; |
| 1765 |
for( j = 1, _do14 = n; j <= _do14; j++ ) |
| 1766 |
{ |
| 1767 |
j_ = j - 1; |
| 1768 |
v[i_][j_] *= s; |
| 1769 |
} |
| 1770 |
if( itere > 1 && vn > 1 ) |
| 1771 |
goto L_330; |
| 1772 |
} |
| 1773 |
/* TRANSFORMATION OF EIGENVECTORS */ |
| 1774 |
L_330: |
| 1775 |
if( n == 2 ) |
| 1776 |
goto L_380; |
| 1777 |
krank = nm2*(n + 1)/2; |
| 1778 |
kpiv = nm2*nm1/2; |
| 1779 |
for( k = nm2; k >= 1; k-- ) |
| 1780 |
{ |
| 1781 |
k_ = k - 1; |
| 1782 |
kp1 = k + 1; |
| 1783 |
if( fabs( A[kpiv] ) <= eps3 ) |
| 1784 |
goto L_360; |
| 1785 |
sum = 0.e0; |
| 1786 |
for( kk = kp1, _do15 = n; kk <= _do15; kk++ ) |
| 1787 |
{ |
| 1788 |
kk_ = kk - 1; |
| 1789 |
sum += A[krank]*v[i_][kk_]; |
| 1790 |
krank += kk; |
| 1791 |
} |
| 1792 |
s = -sum/A[kpiv]; |
| 1793 |
for( kk = n, _do16 = kp1; kk >= _do16; kk-- ) |
| 1794 |
{ |
| 1795 |
kk_ = kk - 1; |
| 1796 |
krank -= kk; |
| 1797 |
v[i_][kk_] += A[krank]*s; |
| 1798 |
} |
| 1799 |
L_360: |
| 1800 |
kpiv -= k; |
| 1801 |
krank -= kp1; |
| 1802 |
} |
| 1803 |
L_380: |
| 1804 |
for( j = ig, _do17 = i; j <= _do17; j++ ) |
| 1805 |
{ |
| 1806 |
j_ = j - 1; |
| 1807 |
if( fabs( E[j] - E[i] ) < eps2 ) |
| 1808 |
goto L_400; |
| 1809 |
} |
| 1810 |
j = i; |
| 1811 |
L_400: |
| 1812 |
ig = j; |
| 1813 |
if( ig == i ) |
| 1814 |
goto L_430; |
| 1815 |
/* RE-ORTHOGONALISATION */ |
| 1816 |
for( k = ig, _do18 = im1; k <= _do18; k++ ) |
| 1817 |
{ |
| 1818 |
k_ = k - 1; |
| 1819 |
sum = 0.e0; |
| 1820 |
for( j = 1, _do19 = n; j <= _do19; j++ ) |
| 1821 |
{ |
| 1822 |
j_ = j - 1; |
| 1823 |
sum += v[k_][j_]*v[i_][j_]; |
| 1824 |
} |
| 1825 |
s = -sum; |
| 1826 |
for( j = 1, _do20 = n; j <= _do20; j++ ) |
| 1827 |
{ |
| 1828 |
j_ = j - 1; |
| 1829 |
v[i_][j_] += v[k_][j_]*s; |
| 1830 |
} |
| 1831 |
} |
| 1832 |
/* NORMALISATION */ |
| 1833 |
L_430: |
| 1834 |
sum = 1.e-24; |
| 1835 |
for( j = 1, _do21 = n; j <= _do21; j++ ) |
| 1836 |
{ |
| 1837 |
j_ = j - 1; |
| 1838 |
sum += (v[i_][j_])*(v[i_][j_]); |
| 1839 |
} |
| 1840 |
sinv = 1.e0/sqrt( sum ); |
| 1841 |
for( j = 1, _do22 = n; j <= _do22; j++ ) |
| 1842 |
{ |
| 1843 |
j_ = j - 1; |
| 1844 |
v[i_][j_] *= sinv; |
| 1845 |
} |
| 1846 |
} |
| 1847 |
return; |
| 1848 |
#undef V |
| 1849 |
} |
