src/gromacs/math/do_fit.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team.
   6  * Copyright (c) 2012,2014, by the GROMACS development team, led by
   7  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   8  * and including many others, as listed in the AUTHORS file in the
   9  * top-level source directory and at http://www.gromacs.org.
  10  *
  11  * GROMACS is free software; you can redistribute it and/or
  12  * modify it under the terms of the GNU Lesser General Public License
  13  * as published by the Free Software Foundation; either version 2.1
  14  * of the License, or (at your option) any later version.
  15  *
  16  * GROMACS is distributed in the hope that it will be useful,
  17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19  * Lesser General Public License for more details.
  20  *
  21  * You should have received a copy of the GNU Lesser General Public
  22  * License along with GROMACS; if not, see
  23  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  24  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  25  *
  26  * If you want to redistribute modifications to GROMACS, please
  27  * consider that scientific software is very special. Version
  28  * control is crucial - bugs must be traceable. We will be happy to
  29  * consider code for inclusion in the official distribution, but
  30  * derived work must not be called official GROMACS. Details are found
  31  * in the README & COPYING files - if they are missing, get the
  32  * official version at http://www.gromacs.org.
  33  *
  34  * To help us fund GROMACS development, we humbly ask that you cite
  35  * the research papers on the package. Check out http://www.gromacs.org.
  36  */
  37 #include "do_fit.h"
  38
  39 #include "gromacs/linearalgebra/nrjac.h"
  40 #include "gromacs/math/utilities.h"
  41 #include "gromacs/math/vec.h"
  42 #include "gromacs/utility/fatalerror.h"
  43 #include "gromacs/utility/smalloc.h"
  44
  45 real calc_similar_ind(gmx_bool bRho, int nind, atom_id *index, real mass[],
  46                       rvec x[], rvec xp[])
  47 {
  48     int  i, j, d;
  49     real m, tm, xs, xd, rs, rd;
  50
  51     tm = 0;
  52     rs = 0;
  53     rd = 0;
  54     for (j = 0; j < nind; j++)
  55     {
  56         if (index)
  57         {
  58             i = index[j];
  59         }
  60         else
  61         {
  62             i = j;
  63         }
  64         m   = mass[i];
  65         tm += m;
  66         for (d = 0; d < DIM; d++)
  67         {
  68             xd  = x[i][d] - xp[i][d];
  69             rd += m * sqr(xd);
  70             if (bRho)
  71             {
  72                 xs  = x[i][d] + xp[i][d];
  73                 rs += m * sqr(xs);
  74             }
  75         }
  76     }
  77     if (bRho)
  78     {
  79         return 2*sqrt(rd/rs);
  80     }
  81     else
  82     {
  83         return sqrt(rd/tm);
  84     }
  85 }
  86
  87 real rmsdev_ind(int nind, atom_id index[], real mass[], rvec x[], rvec xp[])
  88 {
  89     return calc_similar_ind(FALSE, nind, index, mass, x, xp);
  90 }
  91
  92 real rmsdev(int natoms, real mass[], rvec x[], rvec xp[])
  93 {
  94     return calc_similar_ind(FALSE, natoms, NULL, mass, x, xp);
  95 }
  96
  97 real rhodev_ind(int nind, atom_id index[], real mass[], rvec x[], rvec xp[])
  98 {
  99     return calc_similar_ind(TRUE, nind, index, mass, x, xp);
 100 }
 101
 102 real rhodev(int natoms, real mass[], rvec x[], rvec xp[])
 103 {
 104     return calc_similar_ind(TRUE, natoms, NULL, mass, x, xp);
 105 }
 106
 107 void calc_fit_R(int ndim, int natoms, real *w_rls, rvec *xp, rvec *x, matrix R)
 108 {
 109     int      c, r, n, j, m, i, irot, s;
 110     double **omega, **om;
 111     double   d[2*DIM], xnr, xpc;
 112     matrix   vh, vk, u;
 113     real     mn;
 114     int      index;
 115     real     max_d;
 116
 117     if (ndim != 3 && ndim != 2)
 118     {
 119         gmx_fatal(FARGS, "calc_fit_R called with ndim=%d instead of 3 or 2", ndim);
 120     }
 121
 122     snew(omega, 2*ndim);
 123     snew(om, 2*ndim);
 124     for (i = 0; i < 2*ndim; i++)
 125     {
 126         snew(omega[i], 2*ndim);
 127         snew(om[i], 2*ndim);
 128     }
 129
 130     for (i = 0; i < 2*ndim; i++)
 131     {
 132         d[i] = 0;
 133         for (j = 0; j < 2*ndim; j++)
 134         {
 135             omega[i][j] = 0;
 136             om[i][j]    = 0;
 137         }
 138     }
 139
 140     /*calculate the matrix U*/
 141     clear_mat(u);
 142     for (n = 0; (n < natoms); n++)
 143     {
 144         if ((mn = w_rls[n]) != 0.0)
 145         {
 146             for (c = 0; (c < ndim); c++)
 147             {
 148                 xpc = xp[n][c];
 149                 for (r = 0; (r < ndim); r++)
 150                 {
 151                     xnr      = x[n][r];
 152                     u[c][r] += mn*xnr*xpc;
 153                 }
 154             }
 155         }
 156     }
 157
 158     /*construct omega*/
 159     /*omega is symmetric -> omega==omega' */
 160     for (r = 0; r < 2*ndim; r++)
 161     {
 162         for (c = 0; c <= r; c++)
 163         {
 164             if (r >= ndim && c < ndim)
 165             {
 166                 omega[r][c] = u[r-ndim][c];
 167                 omega[c][r] = u[r-ndim][c];
 168             }
 169             else
 170             {
 171                 omega[r][c] = 0;
 172                 omega[c][r] = 0;
 173             }
 174         }
 175     }
 176
 177     /*determine h and k*/
 178     jacobi(omega, 2*ndim, d, om, &irot);
 179     /*real   **omega = input matrix a[0..n-1][0..n-1] must be symmetric
 180      * int     natoms = number of rows and columns
 181      * real      NULL = d[0]..d[n-1] are the eigenvalues of a[][]
 182      * real       **v = v[0..n-1][0..n-1] contains the vectors in columns
 183      * int      *irot = number of jacobi rotations
 184      */
 185
 186     if (debug && irot == 0)
 187     {
 188         fprintf(debug, "IROT=0\n");
 189     }
 190
 191     index = 0; /* For the compiler only */
 192
 193     /* Copy only the first ndim-1 eigenvectors */
 194     for (j = 0; j < ndim-1; j++)
 195     {
 196         max_d = -1000;
 197         for (i = 0; i < 2*ndim; i++)
 198         {
 199             if (d[i] > max_d)
 200             {
 201                 max_d = d[i];
 202                 index = i;
 203             }
 204         }
 205         d[index] = -10000;
 206         for (i = 0; i < ndim; i++)
 207         {
 208             vh[j][i] = M_SQRT2*om[i][index];
 209             vk[j][i] = M_SQRT2*om[i+ndim][index];
 210         }
 211     }
 212     if (ndim == 3)
 213     {
 214         /* Calculate the last eigenvector as the outer-product of the first two.
 215          * This insures that the conformation is not mirrored and
 216          * prevents problems with completely flat reference structures.
 217          */
 218         cprod(vh[0], vh[1], vh[2]);
 219         cprod(vk[0], vk[1], vk[2]);
 220     }
 221     else if (ndim == 2)
 222     {
 223         /* Calculate the last eigenvector from the first one */
 224         vh[1][XX] = -vh[0][YY];
 225         vh[1][YY] =  vh[0][XX];
 226         vk[1][XX] = -vk[0][YY];
 227         vk[1][YY] =  vk[0][XX];
 228     }
 229
 230     /* determine R */
 231     clear_mat(R);
 232     for (r = 0; r < ndim; r++)
 233     {
 234         for (c = 0; c < ndim; c++)
 235         {
 236             for (s = 0; s < ndim; s++)
 237             {
 238                 R[r][c] += vk[s][r]*vh[s][c];
 239             }
 240         }
 241     }
 242     for (r = ndim; r < DIM; r++)
 243     {
 244         R[r][r] = 1;
 245     }
 246
 247     for (i = 0; i < 2*ndim; i++)
 248     {
 249         sfree(omega[i]);
 250         sfree(om[i]);
 251     }
 252     sfree(omega);
 253     sfree(om);
 254 }
 255
 256 void do_fit_ndim(int ndim, int natoms, real *w_rls, rvec *xp, rvec *x)
 257 {
 258     int    i, j, m, r, c;
 259     matrix R;
 260     rvec   x_old;
 261
 262     /* Calculate the rotation matrix R */
 263     calc_fit_R(ndim, natoms, w_rls, xp, x, R);
 264
 265     /*rotate X*/
 266     for (j = 0; j < natoms; j++)
 267     {
 268         for (m = 0; m < DIM; m++)
 269         {
 270             x_old[m] = x[j][m];
 271         }
 272         for (r = 0; r < DIM; r++)
 273         {
 274             x[j][r] = 0;
 275             for (c = 0; c < DIM; c++)
 276             {
 277                 x[j][r] += R[r][c]*x_old[c];
 278             }
 279         }
 280     }
 281 }
 282
 283 void do_fit(int natoms, real *w_rls, rvec *xp, rvec *x)
 284 {
 285     do_fit_ndim(3, natoms, w_rls, xp, x);
 286 }
 287
 288 void reset_x_ndim(int ndim, int ncm, const atom_id *ind_cm,
 289                   int nreset, const atom_id *ind_reset,
 290                   rvec x[], const real mass[])
 291 {
 292     int  i, m, ai;
 293     rvec xcm;
 294     real tm, mm;
 295
 296     if (ndim > DIM)
 297     {
 298         gmx_incons("More than 3 dimensions not supported.");
 299     }
 300     tm = 0.0;
 301     clear_rvec(xcm);
 302     if (ind_cm != NULL)
 303     {
 304         for (i = 0; i < ncm; i++)
 305         {
 306             ai = ind_cm[i];
 307             mm = mass[ai];
 308             for (m = 0; m < ndim; m++)
 309             {
 310                 xcm[m] += mm*x[ai][m];
 311             }
 312             tm += mm;
 313         }
 314     }
 315     else
 316     {
 317         for (i = 0; i < ncm; i++)
 318         {
 319             mm = mass[i];
 320             for (m = 0; m < ndim; m++)
 321             {
 322                 xcm[m] += mm*x[i][m];
 323             }
 324             tm += mm;
 325         }
 326     }
 327     for (m = 0; m < ndim; m++)
 328     {
 329         xcm[m] /= tm;
 330     }
 331
 332     if (ind_reset != NULL)
 333     {
 334         for (i = 0; i < nreset; i++)
 335         {
 336             rvec_dec(x[ind_reset[i]], xcm);
 337         }
 338     }
 339     else
 340     {
 341         for (i = 0; i < nreset; i++)
 342         {
 343             rvec_dec(x[i], xcm);
 344         }
 345     }
 346 }
 347
 348 void reset_x(int ncm, const atom_id *ind_cm,
 349              int nreset, const atom_id *ind_reset,
 350              rvec x[], const real mass[])
 351 {
 352     reset_x_ndim(3, ncm, ind_cm, nreset, ind_reset, x, mass);
 353 }