#include #include #include #include #include #include "nrutil.h" #include "match.h" #define XEPS 1.0e-7 void ls_fit(double **sxyz, double **exyz, long n, double *rms_value, double **fitted_xyz, double **R, double *orgi) /* least-squares fitting between two structures: input xyz pairs the pair number, get the rest. */ { double temp; double ave_exyz[4], ave_sxyz[4], D[5]; double **N, **U, **V; long i, j; /* for(i=1; i<=n ; i++) { fprintf(fout,"%5d %8.3f%8.3f%8.3f %8.3f%8.3f%8.3f\n", i, exyz[i][1],exyz[i][2],exyz[i][3], sxyz[i][1],sxyz[i][2],sxyz[i][3]); } */ if (n < 3){ nrerror("too few atoms for least-squares fitting"); *rms_value= 999.0; return; } /* get the covariance matrix U */ U = dmatrix(1, 3, 1, 3); cov_matrix(sxyz, exyz, n, 3, U); /* get 4-by-4 symmetric matrix N */ N = dmatrix(1, 4, 1, 4); N[1][1] = U[1][1] + U[2][2] + U[3][3]; N[2][2] = U[1][1] - U[2][2] - U[3][3]; N[3][3] = -U[1][1] + U[2][2] - U[3][3]; N[4][4] = -U[1][1] - U[2][2] + U[3][3]; N[1][2] = U[2][3] - U[3][2]; N[2][1] = N[1][2]; N[1][3] = U[3][1] - U[1][3]; N[3][1] = N[1][3]; N[1][4] = U[1][2] - U[2][1]; N[4][1] = N[1][4]; N[2][3] = U[1][2] + U[2][1]; N[3][2] = N[2][3]; N[2][4] = U[3][1] + U[1][3]; N[4][2] = N[2][4]; N[3][4] = U[2][3] + U[3][2]; N[4][3] = N[3][4]; /* get N's eigenvalues and eigenvectors */ V = dmatrix(1, 4, 1, 4); jacobi(N, 4, D, V); /* get the rotation matrix */ for (i = 1; i <= 4; i++) for (j = 1; j <= 4; j++) N[i][j] = V[i][4] * V[j][4]; R[1][1] = N[1][1] + N[2][2] - N[3][3] - N[4][4]; R[1][2] = 2 * (N[2][3] - N[1][4]); R[1][3] = 2 * (N[2][4] + N[1][3]); R[2][1] = 2 * (N[3][2] + N[1][4]); R[2][2] = N[1][1] - N[2][2] + N[3][3] - N[4][4]; R[2][3] = 2 * (N[3][4] - N[1][2]); R[3][1] = 2 * (N[4][2] - N[1][3]); R[3][2] = 2 * (N[4][3] + N[1][2]); R[3][3] = N[1][1] - N[2][2] - N[3][3] + N[4][4]; ave_dmatrix(sxyz, n, 3, ave_sxyz); ave_dmatrix(exyz, n, 3, ave_exyz); /* fitted sxyz origin */ for (i = 1; i <= 3; i++) orgi[i] = ave_exyz[i] - dot(ave_sxyz, R[i]); /* fitted sxyz coordinates */ for (i = 1; i <= n; i++) for (j = 1; j <= 3; j++) fitted_xyz[i][j] = dot(sxyz[i], R[j]) + orgi[j]; /* rms deviation */ temp = 0.0; for (i = 1; i <= n; i++) { for (j = 1; j <= 3; j++) D[j] = exyz[i][j] - fitted_xyz[i][j]; temp += dot(D, D); } *rms_value = sqrt(temp / n); free_dmatrix(U, 1, 3, 1, 3); free_dmatrix(N, 1, 4, 1, 4); free_dmatrix(V, 1, 4, 1, 4); } void cov_matrix(double **a, double **b, long nr, long nc, double **cmtx) /* calculate the covariance matrix between two matrices */ { double ave_a[4], ave_b[4]; double **ta, **ta_x_b; long i, j; ave_dmatrix(a, nr, nc, ave_a); ave_dmatrix(b, nr, nc, ave_b); ta = dmatrix(1, nc, 1, nr); /* transpose of a */ ta_x_b = dmatrix(1, nc, 1, nc); /* transpose-a multiply b */ transpose_matrix(a, nr, nc, ta); multi_matrix(ta, nc, nr, b, nr, nc, ta_x_b); for (i = 1; i <= nc; i++) for (j = 1; j <= nc; j++) cmtx[i][j] = (ta_x_b[i][j] - ave_a[i] * ave_b[j] * nr) / (nr - 1); free_dmatrix(ta, 1, nc, 1, nr); free_dmatrix(ta_x_b, 1, nc, 1, nc); } void ave_dmatrix(double **d, long nr, long nc, double *avedm) { long i, j; for (i = 1; i <= nc; i++) { avedm[i] = 0.0; for (j = 1; j <= nr; j++) avedm[i] += d[j][i]; avedm[i] /= nr; } } void transpose_matrix(double **a, long nr, long nc, double **o) { long i, j; for (i = 1; i <= nc; i++) for (j = 1; j <= nr; j++) o[i][j] = a[j][i]; } void multi_matrix(double **a, long nra, long nca, double **b, long nrb, long ncb, double **o) { long i, j, k; if (nca != nrb) nrerror("matrices a and b do not conform"); for (i = 1; i <= nra; i++) { for (j = 1; j <= ncb; j++) { o[i][j] = 0.0; for (k = 1; k <= nca; k++) o[i][j] += a[i][k] * b[k][j]; } } } double dot(double *va, double *vb) /* dot product between two 1-by-3 vectors */ { double dsum = 0.0; long i; for (i = 1; i <= 3; i++) dsum += va[i] * vb[i]; return dsum; } void jacobi(double **a, long n, double *d, double **v) { long i, j, iq, ip; double tresh, theta, tau, t, sm, s, h, g, c, *b, *z; b = dvector(1, n); z = dvector(1, n); identity_matrix(v, n); for (ip = 1; ip <= n; ip++) { b[ip] = d[ip] = a[ip][ip]; z[ip] = 0.0; } for (i = 1; i <= 100; i++) { sm = 0.0; for (ip = 1; ip <= n - 1; ip++) { for (iq = ip + 1; iq <= n; iq++) sm += fabs(a[ip][iq]); } if (sm < XEPS) { free_dvector(z, 1, n); free_dvector(b, 1, n); eigsrt(d, v, n); return; } if (i < 4) tresh = 0.2 * sm / (n * n); else tresh = 0.0; for (ip = 1; ip <= n - 1; ip++) { for (iq = ip + 1; iq <= n; iq++) { g = 100.0 * fabs(a[ip][iq]); if (i > 4 && (fabs(d[ip]) + g) == fabs(d[ip]) && (fabs(d[iq]) + g) == fabs(d[iq])) a[ip][iq] = 0.0; else if (fabs(a[ip][iq]) > tresh) { h = d[iq] - d[ip]; if ((fabs(h) + g) == fabs(h)) t = a[ip][iq] / h; else { theta = 0.5 * h / a[ip][iq]; t = 1.0 / (fabs(theta) + sqrt(1.0 + theta * theta)); if (theta < 0.0) t = -t; } c = 1.0 / sqrt(1 + t * t); s = t * c; tau = s / (1.0 + c); h = t * a[ip][iq]; z[ip] -= h; z[iq] += h; d[ip] -= h; d[iq] += h; a[ip][iq] = 0.0; for (j = 1; j <= ip - 1; j++) rotate(a, j, ip, j, iq, &g, &h, s, tau); for (j = ip + 1; j <= iq - 1; j++) rotate(a, ip, j, j, iq, &g, &h, s, tau); for (j = iq + 1; j <= n; j++) rotate(a, ip, j, iq, j, &g, &h, s, tau); for (j = 1; j <= n; j++) rotate(v, j, ip, j, iq, &g, &h, s, tau); } } } for (ip = 1; ip <= n; ip++) { b[ip] += z[ip]; d[ip] = b[ip]; z[ip] = 0.0; } } nrerror("too many iterations"); } void identity_matrix(double **d, long n) { long i, j; for (i = 1; i <= n; i++) { for (j = 1; j <= n; j++) d[i][j] = 0.0; d[i][i] = 1.0; } } void rotate(double **a, long i, long j, long k, long l, double *g, double *h, double s, double tau) { *g = a[i][j]; *h = a[k][l]; a[i][j] = *g - s * (*h + *g * tau); a[k][l] = *h + s * (*g - *h * tau); } void eigsrt(double *d, double **v, long n) /* sort eigenvalues into ascending order and rearrange eigenvectors */ { double p; long i, j, k; for (i = 1; i < n; i++) { p = d[k = i]; for (j = i + 1; j <= n; j++) if (d[j] < p) p = d[k = j]; if (k != i) { d[k] = d[i]; d[i] = p; for (j = 1; j <= n; j++) dswap(&v[j][i], &v[j][k]); } } } void dswap(double *pa, double *pb) { double temp; temp = *pa; *pa = *pb; *pb = temp; }