/* _______________________________________________________________________ * * RDPARM/PTRAJ: APR 2000 * _______________________________________________________________________ * * $Header: /thr/gamow/cvsroot/amber7/src/ptraj/rms.c,v 1.6 2000/05/10 21:53:20 cheatham Exp $ * * Revision: $Revision: 1.6 $ * Date: $Date: 2000/05/10 21:53:20 $ * Last checked in by $Author: cheatham $ * * This source is now archived under CVS at Scripps in the amber7 tree. * * NOTE: this is a beta, pre-release version, is under constant development, * probably contains some bugs and is under constant revision; therefore * take care. Please report bugs (and fixes!) to cheatham@chpc.utah.edu or * cheatham@cgl.ucsf.edu * * Do not distribute this code without explicit permission. * Do not incorporate into other codes without explicit permission. * * Current contact information: * * Thomas Cheatham, III * 2000 East, 30 South, Skaggs Hall 201 * Department of Medicinal Chemistry * University of Utah * Salt Lake City, UT 84112-5820 * cheatham@chpc.utah.edu * (801) 587-9653 * FAX: (801) 585-5366 * * Other contributors: * * David Case (Scripps) * Michael Crowley (Scripps) * Jed Pitera (UCSF) * Vickie Tsui (Scripps) * _______________________________________________________________________ */ #include #include #include "utility.h" #include "vector.h" /* * Define the RMS routines; other ptraj functionality isn't necessary! */ void normalize( double a[3] ) { double b; b = 1.0/sqrt((double)(a[0]*a[0] + a[1]*a[1] + a[2]*a[2])); a[0] *= b; a[1] *= b; a[2] *= b; } int diagEsort(double *mat, double *Emat, double *Evec[], double *Eigenvalue) { int njrot; int i, j, k, i3; double eigenvector[9], *eA, v; if (!jacobi3(mat, Eigenvalue, eigenvector, &njrot)) { printf("convergence failed\n"); return(0); } for (i=i3=0; i<3; i++, i3+=3) for (j=0; j<3; j++) Emat[i3+j] = eigenvector[j*3+i]; for (i=0; i<3; i++) Evec[i] = (double *) &Emat[i*3]; for (i=0; i<2; i++) { v = Eigenvalue[k=i]; for (j=i+1; j<3; j++) if (Eigenvalue[j] > v) v = Eigenvalue[k=j]; if (k != i) { Eigenvalue[k] = Eigenvalue[i]; Eigenvalue[i] = v; eA = Evec[i]; Evec[i] = Evec[k]; Evec[k] = eA; } } return(1); } /* * jacobi3() - get jacobian of 3x3 matrix */ #define ROTATE(ARR,MAJ1,MIN1,MAJ2,MIN2) { \ g = ARR[MAJ1 + MIN1]; \ h = ARR[MAJ2 + MIN2]; \ ARR[MAJ1 + MIN1] = g - s*(h+g*tau); \ ARR[MAJ2 + MIN2] = h + s*(g-h*tau); } int jacobi3(double *a, double *d, double *v, int *nrot) /* n must be 3. see b[3] and z[3] below */ { int i, j, ip, iq, p3, j3; double tresh, theta, tau, t, sm, s, h, g, c, b[3], z[3]; for (ip=p3=0; ip<3; ip++,p3+=3) { /* initialize the identity matrix */ for (iq=0; iq<3; iq++) v[p3 + iq] = 0.0; v[p3 + ip] = 1.0; /* initialize b and d to diagonal of a */ b[ip] = d[ip] = a[p3 + ip]; z[ip] = 0.0; } *nrot = 0; for (i=0; i<50; i++) { /* 50 tries */ sm = 0.0; for (ip=p3=0; ip<2; ip++,p3+=3) { for (iq=ip+1; iq<3; iq++) sm += fabs(a[p3 + iq]); } if (sm == 0.0) { return(1); } if (i < 3) tresh = sm * 0.2 / 9.0; /* on 1st three sweeps... */ else tresh = 0.0; /* thereafter... */ for (ip=p3=0; ip<2; ip++,p3+=3) { for (iq=ip+1; iq<3; iq++) { g = 100.0 * fabs(a[p3 + iq]); if ( i > 3 && fabs(d[ip])+g == fabs(d[ip]) && fabs(d[iq])+g == fabs(d[iq])) { a[p3 + iq] = 0.0; } else if (fabs(a[p3 + iq]) > tresh) { h = d[iq]-d[ip]; if (fabs(h)+g==fabs(h)) t = a[p3 + iq] / h; else { theta = 0.5 * h / a[p3 + iq]; t = 1.0 / (fabs(theta)+ (double)sqrt(1.0+theta*theta)); if (theta < 0.0) t = -t; } c = 1.0 / (double)sqrt(1+t*t); s = t * c; tau = s / (1.0+c); h = t * a[p3 + iq]; z[ip] -= h; z[iq] += h; d[ip] -= h; d[iq] += h; a[p3 + iq] = 0.0; for (j=j3=0; j<=ip-1; j++,j3+=3) ROTATE(a,j3,ip,j3,iq) for (j=ip+1; j<=iq-1; j++) ROTATE(a,p3,j,j*3,iq) for (j=iq+1; j<3; j++) ROTATE(a,p3,j,iq*3,j) for (j3=0; j3<9; j3+=3) ROTATE(v,j3,ip,j3,iq) ++(*nrot); } } } for (ip=0; ip<3; ip++) { b[ip] += z[ip]; d[ip] = b[ip]; z[ip] = 0.0; } } printf("Too many iterations in routine JACOBI\n"); return(0); } double rms(int n, int mode, double *mass, int *mask, double *toFitX, double *toFitY, double *toFitZ, double *X, double *Y, double *Z, double rotation[3][3], double translation[3], int fit) { int ierr=0; int i, modifiedCount; char *err; double rms_return; double *weights; double rot[9], rtr[9]; int i3, k3, j, k; double mwss; double b[9], U[9]; double *Evector[3], Eigenvalue[3], Emat[9]; double x, y, z, xx, yy, zz; double total_mass; double sig3; double cp[3]; double cofmX, cofmY, cofmZ; double cofmX1, cofmY1, cofmZ1; #ifdef DEBUG_RMS double cofmXn, cofmYn, cofmZn; double cofmX1n, cofmY1n, cofmZ1n; #endif weights = safe_malloc( sizeof(double) * n ); total_mass = 0.0; for (i=0, modifiedCount=n; i < n; i++) { if ( mask != NULL && mask[i] == 0 ) { weights[i] = 0.0; modifiedCount--; } else if (mass != NULL) weights[i] = mass[i]; else weights[i] = 1.0; total_mass += weights[i]; } if ( mode ) if ( rotation == NULL || translation == NULL ) error("rms", "rotation matrix and translation vector are NULL?"); if ( modifiedCount > 2 ) { bzero(rot, sizeof(double) * 9); bzero(rtr, sizeof(double) * 9); bzero(U, sizeof(double) * 9); cofmX = 0.0; cofmY = 0.0; cofmZ = 0.0; cofmX1 = 0.0; cofmY1 = 0.0; cofmZ1 = 0.0; /* * first shift the center of mass of all the atoms to be fit * to the origin... */ for (k=0; k < n; k++) { cofmX += weights[k] * toFitX[k]; cofmY += weights[k] * toFitY[k]; cofmZ += weights[k] * toFitZ[k]; cofmX1 += weights[k] * X[k]; cofmY1 += weights[k] * Y[k]; cofmZ1 += weights[k] * Z[k]; } cofmX /= total_mass; cofmY /= total_mass; cofmZ /= total_mass; cofmX1 /= total_mass; cofmY1 /= total_mass; cofmZ1 /= total_mass; if (fit) { for (k=0; k < n; k++) { toFitX[k] -= cofmX; toFitY[k] -= cofmY; toFitZ[k] -= cofmZ; X[k] -= cofmX1; Y[k] -= cofmY1; Z[k] -= cofmZ1; } cofmX = 0.0; cofmY = 0.0; cofmZ = 0.0; cofmX1 = 0.0; cofmY1 = 0.0; cofmZ1 = 0.0; } mwss = 0.0; for (k=0; k < n; k++) { x = toFitX[k]-cofmX; y = toFitY[k]-cofmY; z = toFitZ[k]-cofmZ; xx = X[k]-cofmX1; yy = Y[k]-cofmY1; zz = Z[k]-cofmZ1; mwss += weights[k] * ( x*x + y*y + z*z + xx*xx + yy*yy + zz*zz ); #ifdef DEBUG_RMS_SERIOUS if ( weights[k] > 0.0 ) printf(" *** MWSS is %g\n", mwss); #endif rot[0] += weights[k] * x * xx; rot[1] += weights[k] * x * yy; rot[2] += weights[k] * x * zz; rot[3] += weights[k] * y * xx; rot[4] += weights[k] * y * yy; rot[5] += weights[k] * y * zz; rot[6] += weights[k] * z * xx; rot[7] += weights[k] * z * yy; rot[8] += weights[k] * z * zz; } mwss *= 0.5; #ifdef DEBUG_RMS_SERIOUS fprintf(stdout, "\n***RMS: mass weighted sum of squares is %g***\n", mwss); fprintf(stdout, " ROT:\n"); fprintf(stdout, " %g %g %g\n", rot[0], rot[1], rot[2]); fprintf(stdout, " %g %g %g\n", rot[3], rot[4], rot[5]); fprintf(stdout, " %g %g %g\n", rot[6], rot[7], rot[8]); #endif rtr[0] = rot[0]*rot[0] + rot[1]*rot[1] + rot[2]*rot[2]; rtr[1] = rot[0]*rot[3] + rot[1]*rot[4] + rot[2]*rot[5]; rtr[2] = rot[0]*rot[6] + rot[1]*rot[7] + rot[2]*rot[8]; rtr[3] = rot[3]*rot[0] + rot[4]*rot[1] + rot[5]*rot[2]; rtr[4] = rot[3]*rot[3] + rot[4]*rot[4] + rot[5]*rot[5]; rtr[5] = rot[3]*rot[6] + rot[4]*rot[7] + rot[5]*rot[8]; rtr[6] = rot[6]*rot[0] + rot[7]*rot[1] + rot[8]*rot[2]; rtr[7] = rot[6]*rot[3] + rot[7]*rot[4] + rot[8]*rot[5]; rtr[8] = rot[6]*rot[6] + rot[7]*rot[7] + rot[8]*rot[8]; if (!diagEsort(rtr, Emat, Evector, Eigenvalue)) return(0); #ifdef DEBUG_RMS_SERIOUS fprintf(stdout, "\n***RMS: post diagEsort\n Eigenvalues are %g %g %g\n", Eigenvalue[0], Eigenvalue[1], Eigenvalue[2]); fprintf(stdout, "Emat"); fprintf(stdout, " %10.8f %10.8f %10.8f\n", Emat[0], Emat[1], Emat[2]); fprintf(stdout, " %10.8f %10.8f %10.8f\n", Emat[3], Emat[4], Emat[5]); fprintf(stdout, " %10.8f %10.8f %10.8f\n", Emat[6], Emat[7], Emat[8]); #endif VOP_3D_COORDS_CROSS_PRODUCT(Evector[2][0], Evector[2][1], Evector[2][2], Evector[0][0], Evector[0][1], Evector[0][2], Evector[1][0], Evector[1][1], Evector[1][2]); #ifdef DEBUG_RMS_SERIOUS fprintf(stdout, "\n***RMS: Evector\n"); fprintf(stdout, " %10.8f %10.8f %10.8f\n", Evector[0][0], Evector[0][1], Evector[0][2]); fprintf(stdout, " %10.8f %10.8f %10.8f\n", Evector[1][0], Evector[1][1], Evector[1][2]); fprintf(stdout, " %10.8f %10.8f %10.8f\n", Evector[2][0], Evector[2][1], Evector[2][2]); #endif /* Evector dot transpose rot */ b[0] = Evector[0][0] * rot[0] + Evector[0][1] * rot[3] + Evector[0][2] * rot[6]; b[1] = Evector[0][0] * rot[1] + Evector[0][1] * rot[4] + Evector[0][2] * rot[7]; b[2] = Evector[0][0] * rot[2] + Evector[0][1] * rot[5] + Evector[0][2] * rot[8]; normalize(&b[0]); b[3] = Evector[1][0] * rot[0] + Evector[1][1] * rot[3] + Evector[1][2] * rot[6]; b[4] = Evector[1][0] * rot[1] + Evector[1][1] * rot[4] + Evector[1][2] * rot[7]; b[5] = Evector[1][0] * rot[2] + Evector[1][1] * rot[5] + Evector[1][2] * rot[8]; normalize(&b[3]); b[6] = Evector[2][0] * rot[0] + Evector[2][1] * rot[3] + Evector[2][2] * rot[6]; b[7] = Evector[2][0] * rot[1] + Evector[2][1] * rot[4] + Evector[2][2] * rot[7]; b[8] = Evector[2][0] * rot[2] + Evector[2][1] * rot[5] + Evector[2][2] * rot[8]; normalize(&b[6]); VOP_3D_COORDS_CROSS_PRODUCT(cp[0], cp[1], cp[2], b[0], b[1], b[2], b[3], b[4], b[5]); if ( (cp[0] * b[6] + cp[1] * b[7] + cp[2] * b[8]) < 0.0 ) sig3 = -1.0; else sig3 = 1.0; b[6] = cp[0]; b[7] = cp[1]; b[8] = cp[2]; for (k=k3=0; k<3; k++,k3+=3) for (i=i3=0;i<3; i++,i3+=3) for (j=0; j<3; j++) { U[i3 + j] += Evector[k][j] * b[k3 + i]; } #ifdef DEBUG_RMS_SERIOUS fprintf(stdout, "\n***RMS: Evector\n"); fprintf(stdout, " %10.8f %10.8f %10.8f\n", Evector[0][0], Evector[0][1], Evector[0][2]); fprintf(stdout, " %10.8f %10.8f %10.8f\n", Evector[1][0], Evector[1][1], Evector[1][2]); fprintf(stdout, " %10.8f %10.8f %10.8f\n", Evector[2][0], Evector[2][1], Evector[2][2]); printf("***RMS: B values are:\n %g %g %g\n %g %g %g\n %g %g %g\n", b[0], b[1], b[2], b[3], b[4], b[5], b[6], b[7], b[8]); fprintf(stdout, "\nDumping rotation matrix in RMS()\n"); fprintf(stdout, " %10.8f %10.8f %10.8f\n", U[0], U[1], U[2]); fprintf(stdout, " %10.8f %10.8f %10.8f\n", U[3], U[4], U[5]); fprintf(stdout, " %10.8f %10.8f %10.8f\n", U[6], U[7], U[8]); #endif rms_return = mwss - sqrt(fabs(Eigenvalue[0])) - sqrt(fabs(Eigenvalue[1])) - sig3 * sqrt(fabs(Eigenvalue[2])); if ( rms_return < 0 ) { #ifdef DEBUG_RMS printf("RMS returned is < 0 before sqrt, setting to zero (%f)\n", rms_return); #endif rms_return = 0.0; } else { rms_return = sqrt( (2.0 * rms_return) / total_mass); } #ifdef DEBUG_RMS printf("RMS returned is %lf, n is %i\n", rms_return, n); #endif if (mode) { rotation[0][0] = U[0]; rotation[0][1] = U[1]; rotation[0][2] = U[2]; rotation[1][0] = U[3]; rotation[1][1] = U[4]; rotation[1][2] = U[5]; rotation[2][0] = U[6]; rotation[2][1] = U[7]; rotation[2][2] = U[8]; translation[0] = cofmX1 - rotation[0][0]*cofmX - rotation[0][1]*cofmY - rotation[0][2]*cofmZ; translation[1] = cofmY1 - rotation[1][0]*cofmX - rotation[1][1]*cofmY - rotation[1][2]*cofmZ; translation[2] = cofmZ1 - rotation[2][0]*cofmX - rotation[2][1]*cofmY - rotation[2][2]*cofmZ; translation[0] = cofmX; translation[1] = cofmY; translation[2] = cofmZ; } } else ierr = -1; if (ierr != 0) { switch (ierr) { case -1: err = "Number of atoms less than 2"; break; case -2: err = "Illegal weights"; break; default: err = "Unknown error"; break; } error("rms", "KRMS_ reported %s\n", err); } safe_free(weights); return (double) rms_return; } double rmsf(int n, int mode, double *mass, int *mask, float *toFitX, float *toFitY, float *toFitZ, float *X, float *Y, float *Z, float rotation[3][3], float translation[3], int fit) { int ierr=0; int i, modifiedCount; char *err; double rms_return; double *weights; double rot[9], rtr[9]; int i3, k3, j, k; double mwss; double b[9], U[9]; double *Evector[3], Eigenvalue[3], Emat[9]; double x, y, z, xx, yy, zz; double total_mass; double sig3; double cp[3]; double cofmX, cofmY, cofmZ; double cofmX1, cofmY1, cofmZ1; #ifdef DEBUG_RMS double cofmXn, cofmYn, cofmZn; double cofmX1n, cofmY1n, cofmZ1n; #endif weights = safe_malloc( sizeof(double) * n ); total_mass = 0.0; for (i=0, modifiedCount=n; i < n; i++) { if ( mask != NULL && mask[i] == 0 ) { weights[i] = 0.0; modifiedCount--; } else if (mass != NULL) weights[i] = mass[i]; else weights[i] = 1.0; total_mass += weights[i]; } if ( mode ) if ( rotation == NULL || translation == NULL ) error("rms", "rotation matrix and translation vector are NULL?"); if ( modifiedCount > 2 ) { bzero(rot, sizeof(double) * 9); bzero(rtr, sizeof(double) * 9); bzero(U, sizeof(double) * 9); cofmX = 0.0; cofmY = 0.0; cofmZ = 0.0; cofmX1 = 0.0; cofmY1 = 0.0; cofmZ1 = 0.0; /* * first shift the center of mass of all the atoms to be fit * to the origin... */ for (k=0; k < n; k++) { cofmX += weights[k] * toFitX[k]; cofmY += weights[k] * toFitY[k]; cofmZ += weights[k] * toFitZ[k]; cofmX1 += weights[k] * X[k]; cofmY1 += weights[k] * Y[k]; cofmZ1 += weights[k] * Z[k]; } cofmX /= total_mass; cofmY /= total_mass; cofmZ /= total_mass; cofmX1 /= total_mass; cofmY1 /= total_mass; cofmZ1 /= total_mass; if (fit) { for (k=0; k < n; k++) { toFitX[k] -= cofmX; toFitY[k] -= cofmY; toFitZ[k] -= cofmZ; X[k] -= cofmX1; Y[k] -= cofmY1; Z[k] -= cofmZ1; } cofmX = 0.0; cofmY = 0.0; cofmZ = 0.0; cofmX1 = 0.0; cofmY1 = 0.0; cofmZ1 = 0.0; } #ifdef DEBUG_RMS /* * recalculated to COM */ cofmXn = 0.0; cofmYn = 0.0; cofmZn = 0.0; cofmX1n = 0.0; cofmY1n = 0.0; cofmZ1n = 0.0; for (k=0; k < n; k++) { cofmXn += weights[k] * (toFitX[k]-cofmX); cofmYn += weights[k] * (toFitY[k]-cofmY); cofmZn += weights[k] * (toFitZ[k]-cofmZ); cofmX1n += weights[k] * (X[k]-cofmX1); cofmY1n += weights[k] * (Y[k]-cofmY1); cofmZ1n += weights[k] * (Z[k]-cofmZ1); } cofmXn /= total_mass; cofmYn /= total_mass; cofmZn /= total_mass; cofmX1n /= total_mass; cofmY1n /= total_mass; cofmZ1n /= total_mass; printf("To fit new COM is %f %f %f\n", cofmXn, cofmYn, cofmZn); printf("X,Y,Z new COM is %f %f %f\n", cofmX1n, cofmY1n, cofmZ1n); #endif mwss = 0.0; for (k=0; k < n; k++) { x = toFitX[k]-cofmX; y = toFitY[k]-cofmY; z = toFitZ[k]-cofmZ; xx = X[k]-cofmX1; yy = Y[k]-cofmY1; zz = Z[k]-cofmZ1; mwss += weights[k] * ( x*x + y*y + z*z + xx*xx + yy*yy + zz*zz ); #ifdef DEBUG_RMS_SERIOUS if ( weights[k] > 0.0 ) printf(" *** MWSS is %g\n", mwss); #endif rot[0] += weights[k] * x * xx; rot[1] += weights[k] * x * yy; rot[2] += weights[k] * x * zz; rot[3] += weights[k] * y * xx; rot[4] += weights[k] * y * yy; rot[5] += weights[k] * y * zz; rot[6] += weights[k] * z * xx; rot[7] += weights[k] * z * yy; rot[8] += weights[k] * z * zz; } mwss *= 0.5; #ifdef DEBUG_RMS_SERIOUS fprintf(stdout, "\n***RMS: mass weighted sum of squares is %g***\n", mwss); fprintf(stdout, " ROT:\n"); fprintf(stdout, " %g %g %g\n", rot[0], rot[1], rot[2]); fprintf(stdout, " %g %g %g\n", rot[3], rot[4], rot[5]); fprintf(stdout, " %g %g %g\n", rot[6], rot[7], rot[8]); #endif rtr[0] = rot[0]*rot[0] + rot[1]*rot[1] + rot[2]*rot[2]; rtr[1] = rot[0]*rot[3] + rot[1]*rot[4] + rot[2]*rot[5]; rtr[2] = rot[0]*rot[6] + rot[1]*rot[7] + rot[2]*rot[8]; rtr[3] = rot[3]*rot[0] + rot[4]*rot[1] + rot[5]*rot[2]; rtr[4] = rot[3]*rot[3] + rot[4]*rot[4] + rot[5]*rot[5]; rtr[5] = rot[3]*rot[6] + rot[4]*rot[7] + rot[5]*rot[8]; rtr[6] = rot[6]*rot[0] + rot[7]*rot[1] + rot[8]*rot[2]; rtr[7] = rot[6]*rot[3] + rot[7]*rot[4] + rot[8]*rot[5]; rtr[8] = rot[6]*rot[6] + rot[7]*rot[7] + rot[8]*rot[8]; if (!diagEsort(rtr, Emat, Evector, Eigenvalue)) return(0); #ifdef DEBUG_RMS_SERIOUS fprintf(stdout, "\n***RMS: post diagEsort\n Eigenvalues are %g %g %g\n", Eigenvalue[0], Eigenvalue[1], Eigenvalue[2]); fprintf(stdout, "Emat"); fprintf(stdout, " %10.8f %10.8f %10.8f\n", Emat[0], Emat[1], Emat[2]); fprintf(stdout, " %10.8f %10.8f %10.8f\n", Emat[3], Emat[4], Emat[5]); fprintf(stdout, " %10.8f %10.8f %10.8f\n", Emat[6], Emat[7], Emat[8]); #endif VOP_3D_COORDS_CROSS_PRODUCT(Evector[2][0], Evector[2][1], Evector[2][2], Evector[0][0], Evector[0][1], Evector[0][2], Evector[1][0], Evector[1][1], Evector[1][2]); #ifdef DEBUG_RMS_SERIOUS fprintf(stdout, "\n***RMS: Evector\n"); fprintf(stdout, " %10.8f %10.8f %10.8f\n", Evector[0][0], Evector[0][1], Evector[0][2]); fprintf(stdout, " %10.8f %10.8f %10.8f\n", Evector[1][0], Evector[1][1], Evector[1][2]); fprintf(stdout, " %10.8f %10.8f %10.8f\n", Evector[2][0], Evector[2][1], Evector[2][2]); #endif /* Evector dot transpose rot */ b[0] = Evector[0][0] * rot[0] + Evector[0][1] * rot[3] + Evector[0][2] * rot[6]; b[1] = Evector[0][0] * rot[1] + Evector[0][1] * rot[4] + Evector[0][2] * rot[7]; b[2] = Evector[0][0] * rot[2] + Evector[0][1] * rot[5] + Evector[0][2] * rot[8]; normalize(&b[0]); b[3] = Evector[1][0] * rot[0] + Evector[1][1] * rot[3] + Evector[1][2] * rot[6]; b[4] = Evector[1][0] * rot[1] + Evector[1][1] * rot[4] + Evector[1][2] * rot[7]; b[5] = Evector[1][0] * rot[2] + Evector[1][1] * rot[5] + Evector[1][2] * rot[8]; normalize(&b[3]); b[6] = Evector[2][0] * rot[0] + Evector[2][1] * rot[3] + Evector[2][2] * rot[6]; b[7] = Evector[2][0] * rot[1] + Evector[2][1] * rot[4] + Evector[2][2] * rot[7]; b[8] = Evector[2][0] * rot[2] + Evector[2][1] * rot[5] + Evector[2][2] * rot[8]; normalize(&b[6]); VOP_3D_COORDS_CROSS_PRODUCT(cp[0], cp[1], cp[2], b[0], b[1], b[2], b[3], b[4], b[5]); if ( (cp[0] * b[6] + cp[1] * b[7] + cp[2] * b[8]) < 0.0 ) sig3 = -1.0; else sig3 = 1.0; b[6] = cp[0]; b[7] = cp[1]; b[8] = cp[2]; for (k=k3=0; k<3; k++,k3+=3) for (i=i3=0;i<3; i++,i3+=3) for (j=0; j<3; j++) { U[i3 + j] += Evector[k][j] * b[k3 + i]; } #ifdef DEBUG_RMS_SERIOUS fprintf(stdout, "\n***RMS: Evector\n"); fprintf(stdout, " %10.8f %10.8f %10.8f\n", Evector[0][0], Evector[0][1], Evector[0][2]); fprintf(stdout, " %10.8f %10.8f %10.8f\n", Evector[1][0], Evector[1][1], Evector[1][2]); fprintf(stdout, " %10.8f %10.8f %10.8f\n", Evector[2][0], Evector[2][1], Evector[2][2]); printf("***RMS: B values are:\n %g %g %g\n %g %g %g\n %g %g %g\n", b[0], b[1], b[2], b[3], b[4], b[5], b[6], b[7], b[8]); fprintf(stdout, "\nDumping rotation matrix in RMS()\n"); fprintf(stdout, " %10.8f %10.8f %10.8f\n", U[0], U[1], U[2]); fprintf(stdout, " %10.8f %10.8f %10.8f\n", U[3], U[4], U[5]); fprintf(stdout, " %10.8f %10.8f %10.8f\n", U[6], U[7], U[8]); #endif rms_return = mwss - sqrt(fabs(Eigenvalue[0])) - sqrt(fabs(Eigenvalue[1])) - sig3 * sqrt(fabs(Eigenvalue[2])); if ( rms_return < 0 ) { #ifdef DEBUG_RMS printf("RMS returned is < 0 before sqrt, setting to zero (%f)\n", rms_return); #endif rms_return = 0.0; } else { rms_return = sqrt( (2.0 * rms_return) / total_mass); } #ifdef DEBUG_RMS printf("RMS returned is %lf, n is %i\n", rms_return, n); #endif if (mode) { rotation[0][0] = U[0]; rotation[0][1] = U[1]; rotation[0][2] = U[2]; rotation[1][0] = U[3]; rotation[1][1] = U[4]; rotation[1][2] = U[5]; rotation[2][0] = U[6]; rotation[2][1] = U[7]; rotation[2][2] = U[8]; translation[0] = cofmX1 - rotation[0][0]*cofmX - rotation[0][1]*cofmY - rotation[0][2]*cofmZ; translation[1] = cofmY1 - rotation[1][0]*cofmX - rotation[1][1]*cofmY - rotation[1][2]*cofmZ; translation[2] = cofmZ1 - rotation[2][0]*cofmX - rotation[2][1]*cofmY - rotation[2][2]*cofmZ; } } else ierr = -1; if (ierr != 0) { switch (ierr) { case -1: err = "Number of atoms less than 2"; break; case -2: err = "Illegal weights"; break; default: err = "Unknown error"; break; } error("rms", "KRMS_ reported %s\n", err); } safe_free(weights); return (double) rms_return; }