src/gromacs/gmxana/fitahx.c

   1 /*
   2  *
   3  *                This source code is part of
   4  *
   5  *                 G   R   O   M   A   C   S
   6  *
   7  *          GROningen MAchine for Chemical Simulations
   8  *
   9  *                        VERSION 3.2.0
  10  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  11  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  12  * Copyright (c) 2001-2004, The GROMACS development team,
  13  * check out http://www.gromacs.org for more information.
  14
  15  * This program is free software; you can redistribute it and/or
  16  * modify it under the terms of the GNU General Public License
  17  * as published by the Free Software Foundation; either version 2
  18  * of the License, or (at your option) any later version.
  19  *
  20  * If you want to redistribute modifications, please consider that
  21  * scientific software is very special. Version control is crucial -
  22  * bugs must be traceable. We will be happy to consider code for
  23  * inclusion in the official distribution, but derived work must not
  24  * be called official GROMACS. Details are found in the README & COPYING
  25  * files - if they are missing, get the official version at www.gromacs.org.
  26  *
  27  * To help us fund GROMACS development, we humbly ask that you cite
  28  * the papers on the package - you can find them in the top README file.
  29  *
  30  * For more info, check our website at http://www.gromacs.org
  31  *
  32  * And Hey:
  33  * Green Red Orange Magenta Azure Cyan Skyblue
  34  */
  35 #ifdef HAVE_CONFIG_H
  36 #include <config.h>
  37 #endif
  38
  39 #include "fitahx.h"
  40 #include "vec.h"
  41 #include "do_fit.h"
  42 #include "smalloc.h"
  43
  44 static void my_calc_xcm(int nbb, atom_id bbind[], rvec x[], rvec xcm)
  45 {
  46     int    i, m, ai;
  47
  48     clear_rvec(xcm);
  49     for (i = 0; (i < nbb); i++)
  50     {
  51         ai = bbind[i];
  52         rvec_inc(xcm, x[ai]);
  53     }
  54     for (m = 0; (m < DIM); m++)
  55     {
  56         xcm[m] /= (nbb);
  57     }
  58 }
  59
  60 static void my_sub_xcm(int nbb, atom_id bbind[], rvec x[], rvec xcm)
  61 {
  62     int i, ai;
  63
  64     for (i = 0; (i < nbb); i++)
  65     {
  66         ai = bbind[i];
  67         rvec_dec(x[ai], xcm);
  68     }
  69 }
  70
  71 real fit_ahx(int nres, t_bb bb[], int natoms, int nall, atom_id allindex[],
  72              rvec x[], int nca,
  73              atom_id caindex[], gmx_bool bFit)
  74 {
  75     static rvec *xref = NULL;
  76     static real *mass = NULL;
  77     const  real  d    = 0.15;  /* Rise per residue (nm)    */
  78     const  real  tw   = 1.745; /* Twist per residue (rad)  */
  79     const  real  rad  = 0.23;  /* Radius of the helix (nm) */
  80     real         phi0, trms, rms;
  81     rvec         dx, xcm;
  82     int          ai, i, nmass;
  83
  84     if (nca < 3)
  85     {
  86         gmx_fatal(FARGS, "Need at least 3 Calphas to fit to, (now %d)...\n", nca);
  87     }
  88
  89     if (xref == NULL)
  90     {
  91         snew(xref, natoms);
  92         snew(mass, natoms);
  93     }
  94     phi0 = 0;
  95     for (i = 0; (i < nca); i++)
  96     {
  97         ai           = caindex[i];
  98         xref[ai][XX] = rad*cos(phi0);
  99         xref[ai][YY] = rad*sin(phi0);
 100         xref[ai][ZZ] = d*i;
 101
 102         /* Set the mass to select that this Calpha contributes to fitting */
 103         mass[ai] = 10.0;
 104 /*#define DEBUG*/
 105 #ifdef DEBUG
 106         fprintf(stderr, "%5d  %8.3f  %8.3f  %8.3f  %8.3f  %8.3f  %8.3f\n", ai,
 107                 x[ai][XX], x[ai][YY], x[ai][ZZ],
 108                 xref[ai][XX], xref[ai][YY], xref[ai][ZZ]);
 109 #endif
 110         phi0 += tw;
 111     }
 112
 113     /* Center the referece around the origin */
 114     my_calc_xcm(nca, caindex, xref, xcm);
 115     my_sub_xcm(nca, caindex, xref, xcm);
 116
 117     if (bFit)
 118     {
 119         /* Now center the to-be-fitted coords around the origin */
 120         my_calc_xcm(nca, caindex, x, xcm);
 121         my_sub_xcm(nall, allindex, x, xcm);
 122     }
 123 #ifdef DEBUG
 124     dump_ahx(nres, bb, xref, box, 0);
 125 #endif
 126
 127     nmass = 0;
 128     for (i = 0; (i < natoms); i++)
 129     {
 130         if (mass[i] > 0)
 131         {
 132             nmass++;
 133         }
 134     }
 135     if (nmass != nca)
 136     {
 137         gmx_fatal(FARGS, "nmass=%d, nca=%d\n", nmass, nca);
 138     }
 139
 140     /* Now call the fitting routine */
 141     if (bFit)
 142     {
 143         do_fit(natoms, mass, xref, x);
 144     }
 145
 146     /* Reset masses and calc rms */
 147     trms = 0.0;
 148     for (i = 0; (i < nres); i++)
 149     {
 150         ai = bb[i].CA;
 151
 152         if (mass[ai] > 0.0)
 153         {
 154             rvec_sub(x[ai], xref[ai], dx);
 155             rms         = iprod(dx, dx);
 156             bb[i].rmsa += sqrt(rms);
 157             bb[i].nrms++;
 158             trms    += rms;
 159             mass[ai] = 0.0;
 160         }
 161     }
 162     return sqrt(trms/nca);
 163 }