src/gmxlib/princ.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  * check out http://www.gromacs.org for more information.
   7  * Copyright (c) 2012, by the GROMACS development team, led by
   8  * David van der Spoel, Berk Hess, Erik Lindahl, and including many
   9  * others, as listed in the AUTHORS file in the top-level source
  10  * directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38 /* This file is completely threadsafe - keep it that way! */
  39 #ifdef HAVE_CONFIG_H
  40 #include <config.h>
  41 #endif
  42
  43 #include "typedefs.h"
  44 #include "vec.h"
  45 #include "smalloc.h"
  46 #include "gstat.h"
  47 #include "nrjac.h"
  48 #include "txtdump.h"
  49 #include "princ.h"
  50
  51 static void m_op(matrix mat,rvec x)
  52 {
  53   rvec xp;
  54   int  m;
  55
  56   for(m=0; (m<DIM); m++)
  57     xp[m]=mat[m][XX]*x[XX]+mat[m][YY]*x[YY]+mat[m][ZZ]*x[ZZ];
  58   fprintf(stderr,"x    %8.3f  %8.3f  %8.3f\n",x[XX],x[YY],x[ZZ]);
  59   fprintf(stderr,"xp   %8.3f  %8.3f  %8.3f\n",xp[XX],xp[YY],xp[ZZ]);
  60   fprintf(stderr,"fac  %8.3f  %8.3f  %8.3f\n",xp[XX]/x[XX],xp[YY]/x[YY],
  61           xp[ZZ]/x[ZZ]);
  62 }
  63
  64 #define NDIM 4
  65
  66 #ifdef DEBUG
  67 static void ptrans(char *s,real **inten,real d[],real e[])
  68 {
  69   int  m;
  70   real n,x,y,z;
  71   for(m=1; (m<NDIM); m++) {
  72     x=inten[m][1];
  73     y=inten[m][2];
  74     z=inten[m][3];
  75     n=x*x+y*y+z*z;
  76     fprintf(stderr,"%8s %8.3f %8.3f %8.3f, norm:%8.3f, d:%8.3f, e:%8.3f\n",
  77             s,x,y,z,sqrt(n),d[m],e[m]);
  78   }
  79   fprintf(stderr,"\n");
  80 }
  81
  82 void t_trans(matrix trans,real d[],real **ev)
  83 {
  84   rvec x;
  85   int  j;
  86   for(j=0; (j<DIM); j++) {
  87     x[XX]=ev[1][j+1];
  88     x[YY]=ev[2][j+1];
  89     x[ZZ]=ev[3][j+1];
  90     m_op(trans,x);
  91     fprintf(stderr,"d[%d]=%g\n",j,d[j+1]);
  92   }
  93 }
  94 #endif
  95
  96 void principal_comp(int n,atom_id index[],t_atom atom[],rvec x[],
  97                     matrix trans,rvec d)
  98 {
  99   int  i,j,ai,m,nrot;
 100   real mm,rx,ry,rz;
 101   double **inten,dd[NDIM],tvec[NDIM],**ev;
 102 #ifdef DEBUG
 103   real e[NDIM];
 104 #endif
 105   real temp;
 106
 107   snew(inten,NDIM);
 108   snew(ev,NDIM);
 109   for(i=0; (i<NDIM); i++) {
 110     snew(inten[i],NDIM);
 111     snew(ev[i],NDIM);
 112     dd[i]=0.0;
 113 #ifdef DEBUG
 114     e[i]=0.0;
 115 #endif
 116   }
 117
 118   for(i=0; (i<NDIM); i++)
 119     for(m=0; (m<NDIM); m++)
 120       inten[i][m]=0;
 121   for(i=0; (i<n); i++) {
 122     ai=index[i];
 123     mm=atom[ai].m;
 124     rx=x[ai][XX];
 125     ry=x[ai][YY];
 126     rz=x[ai][ZZ];
 127     inten[0][0]+=mm*(sqr(ry)+sqr(rz));
 128     inten[1][1]+=mm*(sqr(rx)+sqr(rz));
 129     inten[2][2]+=mm*(sqr(rx)+sqr(ry));
 130     inten[1][0]-=mm*(ry*rx);
 131     inten[2][0]-=mm*(rx*rz);
 132     inten[2][1]-=mm*(rz*ry);
 133   }
 134   inten[0][1]=inten[1][0];
 135   inten[0][2]=inten[2][0];
 136   inten[1][2]=inten[2][1];
 137 #ifdef DEBUG
 138   ptrans("initial",inten,dd,e);
 139 #endif
 140
 141   for(i=0; (i<DIM); i++) {
 142     for(m=0; (m<DIM); m++)
 143       trans[i][m]=inten[i][m];
 144   }
 145
 146   /* Call numerical recipe routines */
 147   jacobi(inten,3,dd,ev,&nrot);
 148 #ifdef DEBUG
 149   ptrans("jacobi",ev,dd,e);
 150 #endif
 151
 152   /* Sort eigenvalues in ascending order */
 153 #define SWAPPER(i)                      \
 154   if (fabs(dd[i+1]) < fabs(dd[i])) {    \
 155     temp=dd[i];                 \
 156     for(j=0; (j<NDIM); j++) tvec[j]=ev[j][i];\
 157     dd[i]=dd[i+1];                      \
 158     for(j=0; (j<NDIM); j++) ev[j][i]=ev[j][i+1];                \
 159     dd[i+1]=temp;                       \
 160     for(j=0; (j<NDIM); j++) ev[j][i+1]=tvec[j];                 \
 161   }
 162   SWAPPER(0)
 163   SWAPPER(1)
 164   SWAPPER(0)
 165 #ifdef DEBUG
 166   ptrans("swap",ev,dd,e);
 167   t_trans(trans,dd,ev);
 168 #endif
 169
 170   for(i=0; (i<DIM); i++) {
 171     d[i]=dd[i];
 172     for(m=0; (m<DIM); m++)
 173       trans[i][m]=ev[m][i];
 174   }
 175
 176   for(i=0; (i<NDIM); i++) {
 177     sfree(inten[i]);
 178     sfree(ev[i]);
 179   }
 180   sfree(inten);
 181   sfree(ev);
 182 }
 183
 184 void rotate_atoms(int gnx,atom_id *index,rvec x[],matrix trans)
 185 {
 186   real   xt,yt,zt;
 187   int    i,ii;
 188
 189   for(i=0; (i<gnx); i++) {
 190     ii=index ? index[i] : i;
 191     xt=x[ii][XX];
 192     yt=x[ii][YY];
 193     zt=x[ii][ZZ];
 194     x[ii][XX]=trans[XX][XX]*xt+trans[XX][YY]*yt+trans[XX][ZZ]*zt;
 195     x[ii][YY]=trans[YY][XX]*xt+trans[YY][YY]*yt+trans[YY][ZZ]*zt;
 196     x[ii][ZZ]=trans[ZZ][XX]*xt+trans[ZZ][YY]*yt+trans[ZZ][ZZ]*zt;
 197   }
 198 }
 199
 200 real calc_xcm(rvec x[],int gnx,atom_id *index,t_atom *atom,rvec xcm,
 201               gmx_bool bQ)
 202 {
 203   int  i,ii,m;
 204   real m0,tm;
 205
 206   clear_rvec(xcm);
 207   tm=0;
 208   for(i=0; (i<gnx); i++) {
 209     ii=index ? index[i] : i;
 210     if (atom) {
 211       if (bQ)
 212         m0=fabs(atom[ii].q);
 213       else
 214         m0=atom[ii].m;
 215     }
 216     else
 217       m0 = 1;
 218     tm+=m0;
 219     for(m=0; (m<DIM); m++)
 220       xcm[m]+=m0*x[ii][m];
 221   }
 222   for(m=0; (m<DIM); m++)
 223     xcm[m]/=tm;
 224
 225   return tm;
 226 }
 227
 228 real sub_xcm(rvec x[],int gnx,atom_id *index,t_atom atom[],rvec xcm,
 229              gmx_bool bQ)
 230 {
 231   int  i,ii;
 232   real tm;
 233
 234   tm=calc_xcm(x,gnx,index,atom,xcm,bQ);
 235   for(i=0; (i<gnx); i++) {
 236     ii=index ? index[i] : i;
 237     rvec_dec(x[ii],xcm);
 238   }
 239   return tm;
 240 }
 241
 242 void add_xcm(rvec x[],int gnx,atom_id *index,rvec xcm)
 243 {
 244   int  i,ii;
 245
 246   for(i=0; (i<gnx); i++) {
 247     ii=index ? index[i] : i;
 248     rvec_inc(x[ii],xcm);
 249   }
 250 }
 251
 252 void orient_princ(t_atoms *atoms,int isize,atom_id *index,
 253                   int natoms, rvec x[], rvec *v, rvec d)
 254 {
 255   int     i,m;
 256   rvec    xcm,prcomp;
 257   matrix  trans;
 258
 259   calc_xcm(x,isize,index,atoms->atom,xcm,FALSE);
 260   for(i=0; i<natoms; i++)
 261     rvec_dec(x[i],xcm);
 262   principal_comp(isize,index,atoms->atom,x,trans,prcomp);
 263   if (d)
 264     copy_rvec(prcomp, d);
 265
 266   /* Check whether this trans matrix mirrors the molecule */
 267   if (det(trans) < 0) {
 268     for(m=0; (m<DIM); m++)
 269       trans[ZZ][m] = -trans[ZZ][m];
 270   }
 271   rotate_atoms(natoms,NULL,x,trans);
 272   if (v) rotate_atoms(natoms,NULL,v,trans);
 273
 274   for(i=0; i<natoms; i++)
 275     rvec_inc(x[i],xcm);
 276 }
 277