/* The program to superimpose atoms of two molecules by quaternion method David J. Heisterberg The Ohio Supercomputer Center 1224 Kinnear Rd. Columbus, OH 43212-1163 (614)292-6036 djh@ccl.net djh@ohstpy.bitnet ohstpy::djh Translated to C from fitest.f program and interfaced with Xmol program by Jan Labanowski, jkl@ccl.net jkl@ohstpy.bitnet ohstpy::jkl To complile: cc -o quatfit quatfit.c -lm Copyright: Ohio Supercomputer Center, David J. Heisterberg, 1990. The program can be copied and distributed freely, provided that this copyright in not removed. You may acknowledge the use of the program in published material as: David J. Heisterberg, 1990, unpublished results. */ /* there was a typo in the formula for RMS. In the part calculating the RMS, there was a line: wnorm += s; while it should be: wnorm += s*s; Martin Lema mlema[-at-]unq.edu.ar was kind to find it. THANKS!!! jkl. It did not affect fitting too much. But the RMS value was not right when weights were different than 1. */ #include #include #include #include #include #define MAXPOINTS 400 #define MAXLINELEN 250 #define OPTIONS "r:f:p:o:s:" /* options -r refmol reference molecule xmol file. If this option is not given, the information is read from standard input. -f fitmol fitted molecule input xmol file If this option is not given, the information is read from standard input. -p pairs file with the list of fitted atom pairs and weights. If this option is not specified, pairs and weights are taken from stdin. If file name "none" is used (i.e. -p none), atoms of the fitted molecule are fitted to atoms of the reference molecule with the same number, and all weights are assumed 1.0. If molecules do not have the same number of atoms, the smaller number of atoms is fitted. -o outmol output file for transformed fitted molecule. If this option is not given, the information is written to standard output. -s statfile file with fit statistics (distances between fitted atoms, etc). If this option is not given, the information is written to standard output. If any files are read from stdin, the order is: refmol, fitmol, pairs. If any files are written to stdout, the order is: outmol, statfile. The file formats are: The refmol, fitmol and outmol files are in the XYZ format used with xmol program from Minnesota Supercomputer Institute. The format is: 1st line: number of atoms 2nd line: title 3rd and next lines have the format depending on the kind of information: T X Y Z (total of 4 columns) T X Y Z C (total of 5 columns) T X Y Z Mx My Mz (total of 7 columns) T X Y Z C Mx My Mz (total of 8 columns) where T is atom type (usually elemnt symbol), X, Y, Z are cartesian coordinates of an atom in Angstroms, C is atomic charge, and Mx, My, Mz are normal modes. The pairs file format is: 1st line: number of pairs 2nd and next lines: Ar Af W where Ar is the atom number of the reference molecule, Af is the atom number of fitted molecule, w is the statistical weight. Weights W are related to the square of expected deviation "sigma" between the reference and fitted molecule atoms and allow to make fit of some atom pairs more tight. W is proportional to 1/sigma^2. The larger the weight, the more tight will be the resulting fit for the given pair. The statfile lists results of the fit with explanation. */ /*==================================================================== CENTER center or translate a molecule. n - number of atoms x - on input - original xyz coordinates of a molecule on output - moved xyz coordinates (see io for modes). w - if io=1, weights of atoms if io=2 or 3, unused io - 1 weighted geometric center of the molecule will be at (0,0,0) 2 molecule will be moved by a vector -o (i.e., components of a vector o will be subtracted from atom coordinates). 3 molecule will be moved by a vector +o (i.e., components of a vector o will be added atom coordinates). o - if io=1, output, center of original coordinates if io=2, input, vector o will be subtracted from atomic coordinates if io=3, input, vector o will be added to atomic coordinates =====================================================================*/ void center(n, x, w, io, o) int n; double x[4][MAXPOINTS]; double w[MAXPOINTS]; int io; double o[4]; { double wnorm, modif; int i; if (io == 2) { modif = -1.0; } else if (io == 3) { modif = 1.0; } else { modif = -1.0; o[1] = 0.0; o[2] = 0.0; o[3] = 0.0; wnorm = 0.0; for (i = 1; i <= n; i++) { o[1] = o[1] + x[1][i] * sqrt(w[i]); o[2] = o[2] + x[2][i] * sqrt(w[i]); o[3] = o[3] + x[3][i] * sqrt(w[i]); wnorm = wnorm + sqrt(w[i]); } o[1] = o[1] / wnorm; o[2] = o[2] / wnorm; o[3] = o[3] / wnorm; } for (i = 1; i <= n; i++) { x[1][i] = x[1][i] + modif*o[1]; x[2][i] = x[2][i] + modif*o[2]; x[3][i] = x[3][i] + modif*o[3]; } } /*================================ ROTMOL rotate a molecule n - number of atoms x - input coordinates y - rotated coordinates y = u * x u - left rotation matrix ==================================*/ void rotmol (n, x, y, u) int n; double x[4][MAXPOINTS]; double y[4][MAXPOINTS]; double u[4][4]; { double yx, yy, yz; int i; for (i = 1; i <= n; i++) { yx = u[1][1] * x[1][i] + u[1][2] * x[2][i] + u[1][3] * x[3][i]; yy = u[2][1] * x[1][i] + u[2][2] * x[2][i] + u[2][3] * x[3][i]; yz = u[3][1] * x[1][i] + u[3][2] * x[2][i] + u[3][3] * x[3][i]; y[1][i] = yx; y[2][i] = yy; y[3][i] = yz; } } /*======================================================= JACOBI Jacobi diagonalizer with sorted output. It is only good for 4x4 matrices. (was too lazy to do pointers...) a - input: matrix to diagonalize v - output: eigenvectors d - output: eigenvalues nrot - input: maximum number of sweeps =========================================================*/ void jacobi (a, d, v, nrot) int nrot; double a[4][4]; double d[4]; double v[4][4]; { double onorm, dnorm; double b, dma, q, t, c, s; double atemp, vtemp, dtemp; int i, j, k, l; for (j = 0; j <= 3; j++) { for (i = 0; i <= 3; i++) { v[i][j] = 0.0; } v[j][j] = 1.0; d[j] = a[j][j]; } for (l = 1; l <= nrot; l++) { dnorm = 0.0; onorm = 0.0; for (j = 0; j <= 3; j++) { dnorm = dnorm + fabs(d[j]); for (i = 0; i <= j - 1; i++) { onorm = onorm + fabs(a[i][j]); } } if((onorm/dnorm) <= 1.0e-12) goto Exit_now; for (j = 1; j <= 3; j++) { for (i = 0; i <= j - 1; i++) { b = a[i][j]; if(fabs(b) > 0.0) { dma = d[j] - d[i]; if((fabs(dma) + fabs(b)) <= fabs(dma)) { t = b / dma; } else { q = 0.5 * dma / b; t = 1.0/(fabs(q) + sqrt(1.0+q*q)); if(q < 0.0) { t = -t; } } c = 1.0/sqrt(t * t + 1.0); s = t * c; a[i][j] = 0.0; for (k = 0; k <= i-1; k++) { atemp = c * a[k][i] - s * a[k][j]; a[k][j] = s * a[k][i] + c * a[k][j]; a[k][i] = atemp; } for (k = i+1; k <= j-1; k++) { atemp = c * a[i][k] - s * a[k][j]; a[k][j] = s * a[i][k] + c * a[k][j]; a[i][k] = atemp; } for (k = j+1; k <= 3; k++) { atemp = c * a[i][k] - s * a[j][k]; a[j][k] = s * a[i][k] + c * a[j][k]; a[i][k] = atemp; } for (k = 0; k <= 3; k++) { vtemp = c * v[k][i] - s * v[k][j]; v[k][j] = s * v[k][i] + c * v[k][j]; v[k][i] = vtemp; } dtemp = c*c*d[i] + s*s*d[j] - 2.0*c*s*b; d[j] = s*s*d[i] + c*c*d[j] + 2.0*c*s*b; d[i] = dtemp; } /* end if */ } /* end for i */ } /* end for j */ } /* end for l */ Exit_now: nrot = l; for (j = 0; j <= 2; j++) { k = j; dtemp = d[k]; for (i = j+1; i <= 3; i++) { if(d[i] < dtemp) { k = i; dtemp = d[k]; } } if(k > j) { d[k] = d[j]; d[j] = dtemp; for (i = 0; i <= 3; i++) { dtemp = v[i][k]; v[i][k] = v[i][j]; v[i][j] = dtemp; } } } } /*========================================== Q2MAT Generate a left rotation matrix from a normalized quaternion INPUT q - normalized quaternion OUTPUT u - the rotation matrix ===========================================*/ void q2mat (q, u) double q[4]; double u[4][4]; { u[1][1] = q[0]*q[0] + q[1]*q[1] - q[2]*q[2] - q[3]*q[3]; u[2][1] = 2.0 * (q[1] * q[2] - q[0] * q[3]); u[3][1] = 2.0 * (q[1] * q[3] + q[0] * q[2]); u[1][2] = 2.0 * (q[2] * q[1] + q[0] * q[3]); u[2][2] = q[0]*q[0] - q[1]*q[1] + q[2]*q[2] - q[3]*q[3]; u[3][2] = 2.0 * (q[2] * q[3] - q[0] * q[1]); u[1][3] = 2.0 *(q[3] * q[1] - q[0] * q[2]); u[2][3] = 2.0 * (q[3] * q[2] + q[0] * q[1]); u[3][3] = q[0]*q[0] - q[1]*q[1] - q[2]*q[2] + q[3]*q[3]; } /*========================================== QTRFIT Find the quaternion, q,[and left rotation matrix, u] that minimizes |qTXq - Y| ^ 2 [|uX - Y| ^ 2] This is equivalent to maximizing Re (qTXTqY). This is equivalent to finding the largest eigenvalue and corresponding eigenvector of the matrix [A2 AUx AUy AUz ] [AUx Ux2 UxUy UzUx] [AUy UxUy Uy2 UyUz] [AUz UzUx UyUz Uz2 ] where A2 = Xx Yx + Xy Yy + Xz Yz Ux2 = Xx Yx - Xy Yy - Xz Yz Uy2 = Xy Yy - Xz Yz - Xx Yx Uz2 = Xz Yz - Xx Yx - Xy Yy AUx = Xz Yy - Xy Yz AUy = Xx Yz - Xz Yx AUz = Xy Yx - Xx Yy UxUy = Xx Yy + Xy Yx UyUz = Xy Yz + Xz Yy UzUx = Xz Yx + Xx Yz The left rotation matrix, u, is obtained from q by u = qT1q INPUT n - number of points x - fitted molecule coordinates y - reference molecule coordinates w - weights OUTPUT q - the best-fit quaternion u - the best-fit left rotation matrix nr - max number of jacobi sweeps =====================================*/ void qtrfit (n, x, y, w, q, u, nr) int n; double x[4][MAXPOINTS]; double y[4][MAXPOINTS]; double w[MAXPOINTS]; double q[4]; double u[4][4]; int nr; { double xxyx, xxyy, xxyz; double xyyx, xyyy, xyyz; double xzyx, xzyy, xzyz; double c[4][4], v[4][4]; double d[4]; int i, j; /* generate the upper triangle of the quadratic form matrix */ xxyx = 0.0; xxyy = 0.0; xxyz = 0.0; xyyx = 0.0; xyyy = 0.0; xyyz = 0.0; xzyx = 0.0; xzyy = 0.0; xzyz = 0.0; for (i = 1; i <= n; i++) { xxyx = xxyx + x[1][i] * y[1][i] * w[i]; xxyy = xxyy + x[1][i] * y[2][i] * w[i]; xxyz = xxyz + x[1][i] * y[3][i] * w[i]; xyyx = xyyx + x[2][i] * y[1][i] * w[i]; xyyy = xyyy + x[2][i] * y[2][i] * w[i]; xyyz = xyyz + x[2][i] * y[3][i] * w[i]; xzyx = xzyx + x[3][i] * y[1][i] * w[i]; xzyy = xzyy + x[3][i] * y[2][i] * w[i]; xzyz = xzyz + x[3][i] * y[3][i] * w[i]; } for(i = 0; i <= 3; i++) { for(j = 0; j <= 3; j++) { c[i][j] = 0.0; } } c[0][0] = xxyx + xyyy + xzyz; c[0][1] = xzyy - xyyz; c[1][1] = xxyx - xyyy - xzyz; c[0][2] = xxyz - xzyx; c[1][2] = xxyy + xyyx; c[2][2] = xyyy - xzyz - xxyx; c[0][3] = xyyx - xxyy; c[1][3] = xzyx + xxyz; c[2][3] = xyyz + xzyy; c[3][3] = xzyz - xxyx - xyyy; /* diagonalize c */ jacobi (c, d, v, nr); /* extract the desired quaternion */ q[0] = v[0][3]; q[1] = v[1][3]; q[2] = v[2][3]; q[3] = v[3][3]; /* generate the rotation matrix */ q2mat (q, u); } /*============================================ scan_line - scans line of XYZ file of Xmol and parses it into individual components. ==============================================*/ int scan_line(file, n_fields, symb, x, y, z, charge, mx, my, mz) FILE *file; int *n_fields; char *symb; double *x, *y, *z, *charge, *mx, *my, *mz; { int n, i; char line[MAXLINELEN]; *x = 0.0; *y = 0.0; *z = 0.0; *charge = 0.0; *mx = 0.0; *my = 0.0; *mz = 0.0; /* get a line of text from the file, exit EOF reached */ if(fgets(line, MAXLINELEN, file) == NULL) { return((int)-1); } /* find number of fields in the line */ i = 1; n = 0; while (line[i] != '\0') { if((!isspace(line[i-1])) && isspace(line[i])) { n++; } i++; } /* check if n_fields is set already */ if(*n_fields == 0) { *n_fields = n; } /* check if n_fields is correct */ if(*n_fields != n) { return((int)-2); } switch (n) { case 4: if(sscanf(line, "%s%le%le%le", symb, x, y, z) != 4) { return((int)-3); } else { return((int)0); } case 5: if(sscanf(line, "%s%le%le%le%le", symb, x, y, z, charge) != 5) { return((int)-3); } else { return((int)0); } case 7: if(sscanf(line, "%s%le%le%le%le%le%le", symb, x, y, z, mx, my, mz) != 7) { return((int)-3); } else { return((int)0); } case 8: if(sscanf(line, "%s%le%le%le%le%le%le%le", symb, x, y, z, charge, mx, my, mz) != 8) { return((int)-3); } else { return((int)0); } default: return((int)-3); } /* end switch */ } /*============================================= write one line of XYZ xmol coordinate file. n_fields - number of columns in the XYZ file of Xmol. 4 At.Symb X Y Z 5 At.Symb X Y Z Charge 7 At.Symb X Y Z X-amplitude Y-amplitude Z-amplitude 8 At.Symb X Y Z Charge X-ampl. Y-ampl. Z-ampl. ===============================================*/ int write_line(file, n_fields, symb, x, y, z, charge, mx, my, mz) FILE *file; int n_fields; char *symb; double x, y, z, charge, mx, my, mz; { switch (n_fields) { case 4: fprintf(file, "%10s%11.6f%11.6f%11.6f\n", symb, x, y, z); break; case 5: fprintf(file, "%10s%11.6f%11.6f%11.6f%11.6f\n", symb, x, y, z, charge); break; case 7: fprintf(file, "%10s%11.6f%11.6f%11.6f%11.6f%11.6f%11.6f\n", symb, x, y, z, mx, my, mz); break; case 8: fprintf(file, "%10s%11.6f%11.6f%11.6f%11.6f%11.6f%11.6f%11.6f\n", symb, x, y, z, charge, mx, my, mz); break; default: return((int)-3); } /* end switch */ return((int)0); } /*=======================================================*/ /* FITEST rigid fit test driver reads in data, fits, and writes out */ int main(argc, argv) int argc; char *argv[]; { int n_fields_r; /* no of fields in xmol file for ref. molec */ int n_fields_f; /* no of fields in xmol file for fit. molec */ int nat_r; /* number of all atoms in reference molecule */ int nat_f; /* number of all atoms in fitted molecule */ char title_r[MAXLINELEN]; /* title line in reference mol. xmol file */ char title_f[MAXLINELEN]; /* title line in fitted mol. xmol file */ char symb_r[MAXPOINTS][10]; /* atom type symbols for reference molecule */ char symb_f[MAXPOINTS][10]; /* atom type symbols for fitted molecule */ char line[MAXLINELEN]; /* scratch space for line */ double xyz_r[4][MAXPOINTS]; /* coordinates for reference molecule */ double xyz_f[4][MAXPOINTS]; /* coordinates for fitted molecule */ double charge_r[MAXPOINTS]; /* charges for reference molecule */ double charge_f[MAXPOINTS]; /* charges for fitted molecule */ double modes_r[4][MAXPOINTS]; /* normal modes for reference */ double modes_f[4][MAXPOINTS]; /* normal modes for fitted */ int npairs; /* no of fitted atom pairs */ int atoms_r[MAXPOINTS]; /* atoms of ref. molecule to be superimposed */ int atoms_f[MAXPOINTS]; /* atoms of fit. molecule to be superimposed */ double ref_xyz[4][MAXPOINTS]; /* ref. molecule atom coordinates to fit */ double fit_xyz[4][MAXPOINTS]; /* fit. molecule atom coordinates to fit */ double weight[MAXPOINTS]; /* fitted atom pair weights */ double ref_center[4]; /* center of ref. molecule fitted atoms */ double fit_center[4]; /* center of ref. molecule fitted atoms */ double q[4]; /* quaternion */ double u[4][4]; /* left rotation matrix for coordinates */ int i, j, ch; /* aux variables */ double s, d, wd, rms, wnorm, dotm; /* aux variables */ FILE *refmol; /* file variable for reference mol xmol file */ FILE *fitmol; /* file variable for reference mol xmol file */ FILE *pairinp; /* file with fitted atom pairs and weights */ FILE *outfile; /* rotated and translated fit. mol. xmol file*/ FILE *statfile; /* file with fit goodness values */ int max_sweeps; /* max number of iterations in jacobi */ int opt; /* option letter */ int read_pairs; /* 1 - read pairs, 0 - make pairs, weights */ static char usage[]= "Usage : quatfit [-r ref] [-f fit] [-p pairs] [-o out] [-s stat]\n"; extern char *optarg; /* option argument from getopt */ /* set defaults */ max_sweeps = 30; read_pairs = 1; refmol = stdin; fitmol = stdin; pairinp = stdin; outfile = stdout; statfile = stdout; while ((opt = getopt(argc, argv, OPTIONS)) != EOF) { switch (opt) { case 'r': if((refmol = fopen(optarg, "r")) == NULL) { fprintf(stderr,"Error: Could not find ref-mol-file: %s\n", optarg); return(1); } break; case 'f': if((fitmol = fopen(optarg, "r")) == NULL) { fprintf(stderr,"Error: Could not find fit-mol-file: %s\n", optarg); return(2); } break; case 'p': if(strcmp(optarg,"none") == 0) { /* if argument is "none" */ read_pairs = 0; break; } if((pairinp = fopen(optarg, "r")) == NULL) { fprintf(stderr, "Error: Could not find pairs-wieight-file: %s\n", optarg); return(3); } break; case 's': if((statfile = fopen(optarg, "w")) == NULL) { fprintf(stderr, "Error: Could not open out-file %s\n", optarg); return(4); } break; case 'o': if((outfile = fopen(optarg, "w")) == NULL) { fprintf(stderr, "Error: Could not open fitted-mol-out-file: %s\n", optarg); return(5); } break; case '?': fprintf(stderr,"Error: %s\n", usage); return(6); } /* end switch */ } /* end while */ /* Now read in the ref molecule */ while (fgets(line, MAXLINELEN, refmol) != NULL) { /* scan until n_at found */ i = 0; while ((ch = line[i++]) != '\0') { if(!isspace(ch)) { i = -1; break; } } if(i == -1) { break; } } if(i != -1) { /* if white space only to the end */ fprintf(stderr, "Error: Error in line 1 of refmol file\n"); return(7); } if(sscanf(line, "%d\n", &nat_r) != 1) { /* n atoms */ fprintf(stderr, "Error: Error in line 1 of refmol file.\n"); return(8); } if(nat_r >= MAXPOINTS) { fprintf(stderr, "Error: Molecule too big. Recompile program with larger MAXPOINTS\n"); return(8); } if(fgets(title_r, MAXLINELEN, refmol) == NULL) { /* title line */ fprintf(stderr, "Error: Error in line 2 of refmol file.\n"); return(9); } /* read coordinates of ref molecule */ n_fields_r = 0; for (i = 1; i <= nat_r; i++) { if(scan_line(refmol, &n_fields_r, &symb_r[i][0], &xyz_r[1][i], &xyz_r[2][i], &xyz_r[3][i], &charge_r[i], &modes_r[1][i], &modes_r[2][i], &modes_r[3][i]) != 0) { fprintf(stderr,"Error: Error in line %d of refmol file.\n", i+2); return(10); } } /* Now read in the fitted molecule */ while (fgets(line, MAXLINELEN, fitmol) != NULL) { /* scan until n_at found */ i = 0; while ((ch = line[i++]) != '\0') { if(!isspace(ch)) { i = -1; break; } } if(i == -1) { break; } } if(i != -1) { /* if white space only to the end */ fprintf(stderr, "Error: Error in line 1 of fitmol file\n"); return(11); } if(sscanf(line, "%d\n", &nat_f) != 1) { fprintf(stderr, "Error: Error in line 1 of fitmol file.\n"); return(12); } if(nat_f >= MAXPOINTS) { fprintf(stderr, "Error: Molecule too big. Recompile program with larger MAXPOINTS\n"); return(13); } if(fgets(title_f, MAXLINELEN, fitmol) == NULL) { fprintf(stderr, "Error: Error in line 2 of fitmol file\n"); return(13); } /* read coordinates of fitted molecule */ n_fields_f = 0; for (i = 1; i <= nat_f; i++) { if(scan_line(fitmol, &n_fields_f, &symb_f[i][0], &xyz_f[1][i], &xyz_f[2][i], &xyz_f[3][i], &charge_f[i], &modes_f[1][i], &modes_f[2][i], &modes_f[3][i]) != 0) { fprintf(stderr,"Error: Error in line %d of fitmol file.\n", i+2); return(14); } } /* Read or set pairs and weights */ if(read_pairs == 1) { while (fgets(line, MAXLINELEN, pairinp) != NULL) { /* skip white space */ i = 0; while ((ch = line[i++]) != '\0') { if(!isspace(ch)) { i = -1; break; } } if(i == -1) { break; } } if(i != -1) { /* if white space only to the end */ fprintf(stderr, "Error: Error in line 1 of pairs and weights file\n"); return(15); } if(sscanf(line, "%d", &npairs) != 1) { fprintf(stderr, "Error: Error in line 1 of file with pairs and weights\n"); return(16); } if(npairs < 2) { fprintf(stderr, "Error: Cannot fit a single atom. Need at least 2\n"); return(16); } for (i = 1; i <= npairs; i++) { if(fscanf(pairinp,"%d%d%le",&atoms_r[i], &atoms_f[i], &weight[i]) != 3) { fprintf(stderr, "Error: Error in line %d of pairs and weights file.\n", i+2); return(17); } if((atoms_r[i] < 1) || (atoms_r[i] > nat_r) || (atoms_f[i] < 1) || (atoms_f[i] > nat_f) ) { fprintf(stderr, "Error: Error in line %d of pairs and weights (number out of range)\n", i+2); return(18); } } } else { /* initialize pairs to consectutive numbers */ npairs = nat_r; if(nat_r > nat_f) { npairs = nat_f; } for(i = 1; i <= npairs; i++) { atoms_r[i] = i; atoms_f[i] = i; weight[i] = 1.0; } } /* extract fitted atoms to tables */ for (i = 1; i <= npairs; i++) { for (j = 1; j <= 3; j++) { ref_xyz[j][i] = xyz_r[j][atoms_r[i]]; fit_xyz[j][i] = xyz_f[j][atoms_f[i]]; } } /* === Atom coordinates are fit in both modes === */ /* center ref molecule fitted atoms around (0,0,0) */ center (npairs, ref_xyz, weight, 1, ref_center); /* center fitted molecule fitted atoms around (0,0,0) */ center (npairs, fit_xyz, weight, 1, fit_center); /* fit specified atoms of fit_molecule to those of ref_molecule */ qtrfit(npairs, fit_xyz, ref_xyz, weight, q, u, max_sweeps); /* subtract coordinates of the center of fitted atoms of the fitted molecule from all atom coordinates of the fitted molecule (note that weight is a dummy parameter) */ center(nat_f, xyz_f, weight, 2, fit_center); /* rotate the fitted molecule by the rotation matrix u */ rotmol(nat_f, xyz_f, xyz_f, u); /* same with set of fitted atoms of the fitted molecule */ rotmol(npairs, fit_xyz, fit_xyz, u); /* if modes given in fitted molecule, rotate the modes too */ if(n_fields_f > 4) { rotmol(nat_f, modes_f, modes_f, u); /* calculate dot product of reference and fitted molecule modes */ if(n_fields_r > 4) { dotm = 0.0; for (i = 1; i <= npairs; i++) { for (j = 1; j <= 3; j++) { dotm += modes_r[j][atoms_r[i]]*modes_f[j][atoms_f[i]]; } } } } /* translate atoms of the fitted molecule to the center of fitted atoms of the reference molecule */ center(nat_f, xyz_f, weight, 3, ref_center); /* same with set of fitted atoms of the fitted molecule */ center(npairs, fit_xyz, weight, 3, ref_center); /* translate fitted atoms of reference molecule to their orig. location */ center(npairs, ref_xyz, weight, 3, ref_center); /* write modified XYZ file for fitted molecule */ fprintf(outfile,"%d\n%s", nat_f, title_f); for (i = 1; i <= nat_f; i++) { write_line(outfile, n_fields_f, symb_f[i], xyz_f[1][i], xyz_f[2][i], xyz_f[3][i], charge_f[i], modes_f[1][i], modes_f[2][i], modes_f[3][i]); } fflush(outfile); /* find distances between fitted and reference atoms and print them in out file */ fprintf(statfile, "\nDistances and weighted distances between fitted atoms\n"); fprintf(statfile,"Ref.At. Fit.At. Distance Dist*sqrt(weight) weight\n"); rms = 0.0; wnorm = 0.0; for (i = 1; i <= npairs; i++) { d = 0.0; for (j = 1; j <= 3; j++) { s = ref_xyz[j][i] - fit_xyz[j][i]; d += s*s; } d = sqrt(d); s = sqrt(weight[i]); wd = s*d; rms += wd*wd; /* wnorm += s; Bug!!! was fixed by Martin Lema */ wnorm += s*s; fprintf(statfile, " %3d %3d %11.6f %11.6f %11.6f\n", atoms_r[i], atoms_f[i], d, wd, weight[i]); } rms = sqrt(rms/wnorm); fprintf(statfile, "\n\nWeighted root mean square=%10.6f\n\n", rms); fprintf(statfile, "\n\nCenter of reference molecule fitted atoms\n"); fprintf(statfile, "Xc = %11.6f Yc = %11.6f Zc = %11.6f\n", ref_center[1], ref_center[2], ref_center[3]); fprintf(statfile, "\n\nCenter of fitted molecule fitted atoms\n"); fprintf(statfile, "Xc = %11.6f Yc = %11.6f Zc = %11.6f\n", fit_center[1], fit_center[2], fit_center[3]); fprintf(statfile,"\n\nLeft rotation matrix\n"); for (i = 1; i <= 3; i++) { fprintf(statfile, " %11.6f %11.6f %11.6f\n", u[1][i], u[2][i], u[3][i]); } if((n_fields_f > 4) && (n_fields_r > 4)) { fprintf(statfile, "\nDot product of normal modes on fitted atom pairs =%11.6f\n", dotm); } fclose(outfile); fclose(statfile); fclose(pairinp); fclose(fitmol); fclose(refmol); return(0); }