src/gmxlib/do_fit.c

   1 /*
   2  * $Id$
   3  *
   4  *       This source code is part of
   5  *
   6  *        G   R   O   M   A   C   S
   7  *
   8  * GROningen MAchine for Chemical Simulations
   9  *
  10  *               VERSION 2.0
  11  *
  12  * Copyright (c) 1991-1999
  13  * BIOSON Research Institute, Dept. of Biophysical Chemistry
  14  * University of Groningen, The Netherlands
  15  *
  16  * Please refer to:
  17  * GROMACS: A message-passing parallel molecular dynamics implementation
  18  * H.J.C. Berendsen, D. van der Spoel and R. van Drunen
  19  * Comp. Phys. Comm. 91, 43-56 (1995)
  20  *
  21  * Also check out our WWW page:
  22  * http://md.chem.rug.nl/~gmx
  23  * or e-mail to:
  24  * gromacs@chem.rug.nl
  25  *
  26  * And Hey:
  27  * Green Red Orange Magenta Azure Cyan Skyblue
  28  */
  29 static char *SRCID_do_fit_c = "$Id$";
  30
  31 #include "maths.h"
  32 #include "sysstuff.h"
  33 #include "typedefs.h"
  34 #include "nrjac.h"
  35 #include "vec.h"
  36 #include "txtdump.h"
  37
  38 #define EPS  1.0e-09
  39
  40 real rmsdev_ind(int nind,atom_id index[],real mass[],rvec x[],rvec xp[])
  41 {
  42   int  i,j,m;
  43   real sqd,totmass;
  44
  45   sqd = 0;
  46   totmass = 0;
  47   for(j=0; j<nind; j++) {
  48     i=index[j];
  49     totmass += mass[i];
  50     for(m=0; m<DIM; m++)
  51       sqd += mass[i]*(x[i][m]-xp[i][m])*(x[i][m]-xp[i][m]);
  52   }
  53
  54   return sqrt(sqd/totmass);
  55 }
  56
  57 real rmsdev(int natoms,real mass[],rvec x[],rvec xp[])
  58 {
  59   int  i,m;
  60   real sqd,totmass;
  61
  62   sqd = 0;
  63   totmass = 0;
  64   for(i=0; i<natoms; i++) {
  65     totmass += mass[i];
  66     for(m=0; m<DIM; m++)
  67       sqd += mass[i]*(x[i][m]-xp[i][m])*(x[i][m]-xp[i][m]);
  68   }
  69
  70   return sqrt(sqd/totmass);
  71 }
  72
  73 void do_fit(int natoms,real *w_rls,rvec *xp,rvec *x)
  74 {
  75   int    c,r,n,j,m,i,irot;
  76   double omega[7][7],om[7][7],d[7],xnr,xpc;
  77   matrix vh,vk,R,u;
  78   /*
  79   matrix vh,vk,R,vh_d,vk_d,u;
  80   */
  81   real   mn;
  82   int    index;
  83   real   max_d;
  84   rvec   x_old;
  85
  86   for(i=0;(i<7);i++) {
  87     d[i]=0;
  88     for(j=0;(j<7);j++) {
  89       omega[i][j]=0;
  90       om[i][j]=0;
  91     }
  92   }
  93
  94   /*calculate the matrix U*/
  95   clear_mat(u);
  96   for(n=0;(n<natoms);n++) {
  97     if ((mn = w_rls[n]) != 0.0) {
  98       for(c=0; (c<DIM); c++) {
  99         xpc=xp[n][c];
 100         for(r=0; (r<DIM); r++) {
 101           xnr=x[n][r];
 102           u[c][r]+=mn*xnr*xpc;
 103         }
 104       }
 105     }
 106   }
 107
 108   /*construct omega*/
 109   /*omega is symmetric -> omega==omega' */
 110   for(r=0;(r<6);r++)
 111     for(c=0;(c<=r);c++)
 112       if ((r>=3) && (c<3)) {
 113         omega[r+1][c+1]=u[r-3][c];
 114         omega[c+1][r+1]=u[r-3][c];
 115       }
 116       else {
 117         omega[r+1][c+1]=0;
 118         omega[c+1][r+1]=0;
 119       }
 120
 121   /*determine h and k*/
 122   jacobi(omega,6,d,om,&irot);
 123   /*real   **omega = input matrix a[1..n][1..n] must be symmetric
 124    *int     natoms = number of rows and columns
 125    *real      NULL = d[1]..d[n] are the eigenvalues of a[][]
 126    *real       **v = v[1..n][1..n] contains the vectors in columns
 127    *int      *irot = number of jacobi rotations
 128    */
 129
 130   if (irot==0) {
 131     fprintf(stderr,"IROT=0\n");
 132   }
 133
 134   index=0; /* For the compiler only */
 135
 136   /* Copy only the first two eigenvectors */
 137   for(j=0;(j<2);j++) {
 138     max_d=-1000;
 139     for(i=0;(i<6);i++)
 140       if (d[i+1]>max_d) {
 141         max_d=d[i+1];
 142         index=i;
 143       }
 144     d[index+1]=-10000;
 145     for(i=0;(i<3);i++) {
 146       vh[j][i]=M_SQRT2*om[i+1][index+1];
 147       vk[j][i]=M_SQRT2*om[i+4][index+1];
 148     }
 149   }
 150   /* Calculate the last eigenvector as the outer-product of the first two.
 151    * This insures that the conformation is not mirrored and
 152    * prevents problems with completely flat reference structures.
 153    */
 154   oprod(vh[0],vh[1],vh[2]);
 155   oprod(vk[0],vk[1],vk[2]);
 156
 157   /*determine R*/
 158   for(c=0;(c<3);c++)
 159     for(r=0;(r<3);r++)
 160       R[c][r]=vk[0][r]*vh[0][c]+
 161               vk[1][r]*vh[1][c]+
 162               vk[2][r]*vh[2][c];
 163
 164   /*rotate X*/
 165   for(j=0;(j<natoms);j++) {
 166     for(m=0;(m<3);m++)
 167       x_old[m]=x[j][m];
 168     for(r=0;(r<3);r++) {
 169       x[j][r]=0;
 170       for(c=0;(c<3);c++)
 171         x[j][r]+=R[c][r]*x_old[c];
 172     }
 173   }
 174 }
 175
 176 void reset_x(int ncm,atom_id ind_cm[],
 177              int nrms,atom_id ind_rms[],
 178              rvec x[],real mass[])
 179 {
 180   int  i,m,ai;
 181   rvec xcm;
 182   real tm,mm;
 183
 184   tm=0.0;
 185   clear_rvec(xcm);
 186   for(i=0; (i<ncm); i++) {
 187     ai=ind_cm[i];
 188     mm=mass[ai];
 189     for(m=0; (m<DIM); m++)
 190       xcm[m]+=mm*x[ai][m];
 191     tm+=mm;
 192   }
 193   for(m=0; (m<DIM); m++)
 194     xcm[m]/=tm;
 195
 196   for(i=0; (i<nrms); i++)
 197     rvec_dec(x[ind_rms[i]],xcm);
 198 }
 199