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