src/tools/gmx_multipoles.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 <math.h>
  40
  41 #include "statutil.h"
  42 #include "macros.h"
  43 #include "tpxio.h"
  44 #include "smalloc.h"
  45 #include "physics.h"
  46 #include "vec.h"
  47 #include "gstat.h"
  48 #include "nrjac.h"
  49 #include "copyrite.h"
  50 #include "index.h"
  51 #include "gmx_ana.h"
  52
  53
  54 #define NM2ANG 10
  55 #define TOLERANCE 1.0E-8
  56
  57 #define e2d(x) ENM2DEBYE*(x)
  58 #define delta(a,b) (( a == b ) ? 1.0 : 0.0)
  59
  60 #define NDIM 3          /* We will be using a numerical recipes routine */
  61
  62 static char dim[DIM+1] = "XYZ";
  63
  64 typedef real            tensor3[DIM][DIM][DIM];      /* 3 rank tensor */
  65 typedef real            tensor4[DIM][DIM][DIM][DIM]; /* 4 rank tensor */
  66
  67
  68 void pr_coord(int k0,int k1,atom_id index[],rvec x[],char *msg)
  69 {
  70   int k,kk;
  71
  72   fprintf(stdout,"Coordinates in nm (%s)\n",msg);
  73   for(k=k0; (k<k1); k++) {
  74     kk=index[k];
  75     fprintf(stdout,"Atom %d, %15.10f %15.10f %15.10f\n",
  76             kk,x[kk][XX],x[kk][YY],x[kk][ZZ]);
  77   }
  78   fprintf(stdout,"\n");
  79 }
  80
  81 static void clear_tensor3(tensor3 a)
  82 {
  83   int i,j,k;
  84   const real nul=0.0;
  85
  86   for(i=0; (i<DIM); i++)
  87     for(j=0; (j<DIM); j++)
  88       for(k=0; (k<DIM); k++)
  89         a[i][j][k]=nul;
  90 }
  91
  92 static void clear_tensor4(tensor4 a)
  93 {
  94   int i,j,k,l;
  95   const real nul=0.0;
  96
  97   for(i=0; (i<DIM); i++)
  98     for(j=0; (j<DIM); j++)
  99       for(k=0; (k<DIM); k++)
 100         for(l=0; (l<DIM); l++)
 101           a[i][j][k][l]=nul;
 102 }
 103
 104 void rotate_mol(int k0,int k1,atom_id index[],rvec x[],matrix trans)
 105 {
 106   real   xt,yt,zt;
 107   int    k,kk;
 108
 109   for(k=k0; (k<k1); k++) {
 110     kk=index[k];
 111     xt=x[kk][XX];
 112     yt=x[kk][YY];
 113     zt=x[kk][ZZ];
 114     x[kk][XX]=trans[XX][XX]*xt+trans[XX][YY]*yt+trans[XX][ZZ]*zt;
 115     x[kk][YY]=trans[YY][XX]*xt+trans[YY][YY]*yt+trans[YY][ZZ]*zt;
 116     x[kk][ZZ]=trans[ZZ][XX]*xt+trans[ZZ][YY]*yt+trans[ZZ][ZZ]*zt;
 117   }
 118 }
 119
 120 /* the following routines are heavily inspired by the Gaussian 94 source
 121  * code
 122  */
 123
 124 /*
 125      Make the rotation matrix for angle Theta counterclockwise about
 126      axis IXYZ.
 127 */
 128
 129 void make_rot_mat(int axis,real theta,matrix t_mat){
 130
 131   ivec i;
 132   real s,c;
 133
 134
 135   i[XX]=axis + 1;
 136   i[YY]=1 + i[XX] % 3;
 137   i[ZZ]=1 + i[YY] % 3;
 138
 139   i[XX]-=1;
 140   i[YY]-=1;
 141   i[ZZ]-=1;
 142
 143   s=sin(theta);
 144   c=cos(theta);
 145   t_mat[i[XX]][i[XX]]=1.0;
 146   t_mat[i[XX]][i[YY]]=0.0;
 147   t_mat[i[XX]][i[ZZ]]=0.0;
 148   t_mat[i[YY]][i[XX]]=0.0;
 149   t_mat[i[YY]][i[YY]]=c;
 150   t_mat[i[YY]][i[ZZ]]=s;
 151   t_mat[i[ZZ]][i[XX]]=0.0;
 152   t_mat[i[ZZ]][i[YY]]=-s;
 153   t_mat[i[ZZ]][i[ZZ]]=c;
 154 }
 155
 156 bool test_linear_mol(rvec d)
 157 {
 158   /* d is sorted in descending order */
 159   if ( (d[ZZ] < TOLERANCE) && (d[XX]-d[YY]) < TOLERANCE ) {
 160     return TRUE;
 161   } else
 162     return FALSE;
 163 }
 164
 165 /* Returns the third moment of charge along an axis */
 166 real test_qmom3(int k0,int k1,atom_id index[],t_atom atom[],rvec x[],int axis){
 167
 168   int k,kk;
 169   real xcq,q;
 170
 171   xcq=0.0;
 172   for(k=k0; (k<k1); k++) {
 173     kk=index[k];
 174     q=fabs(atom[kk].q);
 175     xcq+=q*x[kk][axis]*x[kk][axis]*x[kk][axis];
 176   }
 177
 178   return xcq;
 179 }
 180
 181 /* Returns the second moment of mass along an axis */
 182 real test_mmom2(int k0,int k1,atom_id index[],t_atom atom[],rvec x[],int axis){
 183
 184   int k,kk;
 185   real xcm,m;
 186
 187   xcm=0.0;
 188   for(k=k0; (k<k1); k++) {
 189     kk=index[k];
 190     m=atom[kk].m;
 191     xcm+=m*x[kk][axis]*x[kk][axis];
 192   }
 193
 194   return xcm;
 195 }
 196
 197 real calc_xcm_mol(int k0,int k1,atom_id index[],t_atom atom[],rvec x[],
 198                   rvec xcm)
 199 {
 200   int  k,kk,m;
 201   real m0,tm;
 202
 203   /* Compute the center of mass */
 204   clear_rvec(xcm);
 205   tm=0.0;
 206   for(k=k0; (k<k1); k++) {
 207     kk=index[k];
 208     m0=atom[kk].m;
 209     tm+=m0;
 210     for(m=0; (m<DIM); m++)
 211       xcm[m]+=m0*x[kk][m];
 212   }
 213   for(m=0; (m<DIM); m++)
 214     xcm[m]/=tm;
 215
 216   /* And make it the origin */
 217   for(k=k0; (k<k1); k++) {
 218     kk=index[k];
 219     rvec_dec(x[kk],xcm);
 220   }
 221
 222   return tm;
 223 }
 224
 225 /* Retruns the center of charge */
 226 real calc_xcq_mol(int k0,int k1,atom_id index[],t_atom atom[],rvec x[],
 227                   rvec xcq)
 228 {
 229   int  k,kk,m;
 230   real q0,tq;
 231
 232   clear_rvec(xcq);
 233   tq=0.0;
 234   for(k=k0; (k<k1); k++) {
 235     kk=index[k];
 236     q0=fabs(atom[kk].q);
 237     tq+=q0;
 238     fprintf(stdout,"tq: %f, q0: %f\n",tq,q0);
 239     for(m=0; (m<DIM); m++)
 240       xcq[m]+=q0*x[kk][m];
 241   }
 242
 243   for(m=0; (m<DIM); m++)
 244     xcq[m]/=tq;
 245   /*
 246   for(k=k0; (k<k1); k++) {
 247     kk=index[k];
 248     rvec_dec(x[kk],xcq);
 249   }
 250   */
 251   return tq;
 252 }
 253
 254 /* Returns in m1 the dipole moment */
 255 void mol_M1(int n0,int n1,atom_id ma[],rvec x[],t_atom atom[],rvec m1)
 256 {
 257   int  m,n,nn;
 258   real q;
 259
 260   clear_rvec(m1);
 261   for(n=n0; (n<n1); n++) {
 262     nn = ma[n];
 263     q  = e2d(atom[nn].q);
 264     for(m=0; (m<DIM); m++)
 265       m1[m] += q*x[nn][m];
 266   }
 267 }
 268
 269 /* returns in m2 the quadrupole moment */
 270 void mol_M2(int n0,int n1,atom_id ma[],rvec x[],t_atom atom[],tensor m2)
 271 {
 272   int  n,nn,i,j;
 273   real q,r2;
 274
 275   clear_mat(m2);
 276   for(n=n0; (n<n1); n++) {
 277     nn = ma[n];
 278     q  = e2d(atom[nn].q);
 279     r2 = norm2(x[nn]);
 280     for(i=0; (i<DIM); i++)
 281       for(j=0; (j<DIM); j++)
 282         m2[i][j] += 0.5*q*(3.0*x[nn][i]*x[nn][j] - r2*delta(i,j))*NM2ANG;
 283   }
 284 }
 285
 286 /* Returns in m3 the octopole moment */
 287 void mol_M3(int n0,int n1,atom_id ma[],rvec x[],t_atom atom[],tensor3 m3)
 288 {
 289   int  i,j,k,n,nn;
 290   real q,r2;
 291
 292   clear_tensor3(m3);
 293   for(n=n0; (n<n1); n++) {
 294     nn = ma[n];
 295     q  = e2d(atom[nn].q);
 296     r2 = norm2(x[nn]);
 297     for(i=0; (i<DIM); i++)
 298       for(j=0; (j<DIM); j++)
 299         for(k=0; (k<DIM); k++)
 300           m3[i][j][k] +=
 301             0.5*q*(5.0*x[nn][i]*x[nn][j]*x[nn][k]
 302                    - r2*(x[nn][i]*delta(j,k) +
 303                          x[nn][j]*delta(k,i) +
 304                          x[nn][k]*delta(i,j)))*NM2ANG*NM2ANG;
 305   }
 306 }
 307
 308 /* Returns in m4 the hexadecapole moment */
 309 void mol_M4(int n0,int n1,atom_id ma[],rvec x[],t_atom atom[],tensor4 m4)
 310 {
 311   int  i,j,k,l,n,nn;
 312   real q,r2;
 313
 314   clear_tensor4(m4);
 315   for(n=n0; (n<n1); n++) {
 316     nn = ma[n];
 317     q  = e2d(atom[nn].q);
 318     r2 = norm2(x[nn]);
 319     for(i=0; (i<DIM); i++)
 320       for(j=0; (j<DIM); j++)
 321         for(k=0; (k<DIM); k++)
 322           for(l=0; (l<DIM); l++)
 323             m4[i][j][k][l] +=
 324               0.125*q*(35.0*x[nn][i]*x[nn][j]*x[nn][k]*x[nn][l]
 325                      - 5.0*r2*(x[nn][i]*x[nn][j]*delta(k,l) +
 326                                x[nn][i]*x[nn][k]*delta(j,l) +
 327                                x[nn][i]*x[nn][l]*delta(j,k) +
 328                                x[nn][j]*x[nn][k]*delta(i,l) +
 329                                x[nn][j]*x[nn][l]*delta(i,k) +
 330                                x[nn][k]*x[nn][l]*delta(i,j))
 331                       + r2*r2*(delta(i,j)*delta(k,l) +
 332                                delta(i,k)*delta(j,l) +
 333                                delta(i,l)*delta(j,k)))*NM2ANG*NM2ANG*NM2ANG;
 334   }
 335 }
 336
 337 /* Print the dipole moment components and the total dipole moment */
 338 void pr_M1(FILE *fp,char *msg,int mol,rvec m1,real time)
 339 {
 340   int  i;
 341   real m1_tot;
 342
 343   fprintf(fp,"Molecule: %d @ t= %f ps\n",mol,time);
 344
 345   m1_tot = sqrt(m1[XX]*m1[XX]+m1[YY]*m1[YY]+m1[ZZ]*m1[ZZ]);
 346
 347   fprintf(stdout,"Dipole Moment %s(Debye):\n",msg);
 348   fprintf(stdout,"X= %10.5f Y= %10.5f Z= %10.5f Tot= %10.5f\n",
 349           m1[XX],m1[YY],m1[ZZ],m1_tot);
 350 }
 351
 352 /* Print the quadrupole moment components */
 353 void pr_M2(FILE *fp,char *msg,tensor m2,bool bFull)
 354 {
 355   int i,j;
 356
 357   fprintf(fp,"Quadrupole Moment %s(Debye-Ang):\n",msg);
 358   if (!bFull) {
 359     fprintf(fp,"XX= %10.5f YY= %10.5f ZZ= %10.5f\n",
 360             m2[XX][XX],m2[YY][YY],m2[ZZ][ZZ]);
 361     fprintf(fp,"XY= %10.5f XZ= %10.5f YZ= %10.5f\n",
 362             m2[XX][YY],m2[XX][ZZ],m2[YY][ZZ]);
 363   }
 364   else {
 365     for(i=0; (i<DIM); i++) {
 366       for(j=0; (j<DIM); j++)
 367         fprintf(fp,"  %c%c= %10.4f",dim[i],dim[j],m2[i][j]);
 368       fprintf(fp,"\n");
 369     }
 370   }
 371 }
 372
 373 /* Print the octopole moment components */
 374 void pr_M3(FILE *fp,char *msg,tensor3 m3,bool bFull)
 375 {
 376   int i,j,k;
 377
 378   fprintf(fp,"Octopole Moment %s(Debye-Ang^2):\n",msg);
 379   if (!bFull) {
 380     fprintf(fp,"XXX= %10.5f YYY= %10.5f ZZZ= %10.5f XYY= %10.5f\n",
 381             m3[XX][XX][XX],m3[YY][YY][YY],m3[ZZ][ZZ][ZZ],m3[XX][YY][YY]);
 382     fprintf(fp,"XXY= %10.5f XXZ= %10.5f XZZ= %10.5f YZZ= %10.5f\n",
 383             m3[XX][XX][YY],m3[XX][XX][ZZ],m3[XX][ZZ][ZZ],m3[YY][ZZ][ZZ]);
 384     fprintf(fp,"YYZ= %10.5f XYZ= %10.5f\n",
 385             m3[YY][YY][ZZ],m3[XX][YY][ZZ]);
 386   }
 387   else {
 388     for(i=0; (i<DIM); i++) {
 389       for(j=0; (j<DIM); j++) {
 390         for(k=0; (k<DIM); k++)
 391           fprintf(fp,"  %c%c%c= %10.4f",dim[i],dim[j],dim[k],m3[i][j][k]);
 392         fprintf(fp,"\n");
 393       }
 394     }
 395   }
 396 }
 397
 398 /* Print the hexadecapole moment components */
 399 void pr_M4(FILE *fp,char *msg,tensor4 m4,bool bFull)
 400 {
 401   int i,j,k,l;
 402
 403   fprintf(fp,"Hexadecapole Moment %s(Debye-Ang^3):\n",msg);
 404   if (!bFull) {
 405     fprintf(fp,"XXXX= %10.5f YYYY= %10.5f ZZZZ= %10.5f XXXY= %10.5f\n",
 406             m4[XX][XX][XX][XX],m4[YY][YY][YY][YY],
 407             m4[ZZ][ZZ][ZZ][ZZ],m4[XX][XX][XX][YY]);
 408     fprintf(fp,"XXXZ= %10.5f YYYX= %10.5f YYYZ= %10.5f ZZZX= %10.5f\n",
 409             m4[XX][XX][XX][ZZ],m4[YY][YY][YY][XX],
 410             m4[YY][YY][YY][ZZ],m4[ZZ][ZZ][ZZ][XX]);
 411     fprintf(fp,"ZZZY= %10.5f XXYY= %10.5f XXZZ= %10.5f YYZZ= %10.5f\n",
 412             m4[ZZ][ZZ][ZZ][YY],m4[XX][XX][YY][YY],
 413             m4[XX][XX][ZZ][ZZ],m4[YY][YY][ZZ][ZZ]);
 414     fprintf(fp,"XXYZ= %10.5f YYXZ= %10.5f ZZXY= %10.5f\n\n",
 415             m4[XX][XX][YY][ZZ],m4[YY][YY][XX][ZZ],m4[ZZ][ZZ][XX][YY]);
 416   }
 417   else {
 418     for(i=0; (i<DIM); i++) {
 419       for(j=0; (j<DIM); j++) {
 420         for(k=0; (k<DIM); k++) {
 421           for(l=0; (l<DIM); l++)
 422             fprintf(fp,"  %c%c%c%c = %10.4f",dim[i],dim[j],dim[k],dim[l],
 423                     m4[i][j][k][l]);
 424           fprintf(fp,"\n");
 425         }
 426       }
 427     }
 428   }
 429 }
 430
 431 /* Compute the inertia tensor and returns in trans a matrix which rotates
 432  * the molecules along the principal axes system */
 433 void principal_comp_mol(int k0,int k1,atom_id index[],t_atom atom[],rvec x[],
 434                         matrix trans,rvec d)
 435 {
 436   int  i,j,ai,m,nrot;
 437   real mm,rx,ry,rz;
 438   double **inten,dd[NDIM],tvec[NDIM],**ev;
 439   real temp;
 440
 441   snew(inten,NDIM);
 442   snew(ev,NDIM);
 443   for(i=0; (i<NDIM); i++) {
 444     snew(inten[i],NDIM);
 445     snew(ev[i],NDIM);
 446     dd[i]=0.0;
 447   }
 448
 449   for(i=0; (i<NDIM); i++)
 450     for(m=0; (m<NDIM); m++)
 451       inten[i][m]=0;
 452
 453   for(i=k0; (i<k1); i++) {
 454     ai=index[i];
 455     mm=atom[ai].m;
 456     rx=x[ai][XX];
 457     ry=x[ai][YY];
 458     rz=x[ai][ZZ];
 459     inten[0][0]+=mm*(sqr(ry)+sqr(rz));
 460     inten[1][1]+=mm*(sqr(rx)+sqr(rz));
 461     inten[2][2]+=mm*(sqr(rx)+sqr(ry));
 462     inten[1][0]-=mm*(ry*rx);
 463     inten[2][0]-=mm*(rx*rz);
 464     inten[2][1]-=mm*(rz*ry);
 465   }
 466   inten[0][1]=inten[1][0];
 467   inten[0][2]=inten[2][0];
 468   inten[1][2]=inten[2][1];
 469
 470   /* Call numerical recipe routines */
 471   jacobi(inten,3,dd,ev,&nrot);
 472
 473   /* Sort eigenvalues in descending order */
 474 #define SWAPPER(i)                      \
 475   if (fabs(dd[i+1]) > fabs(dd[i])) {    \
 476     temp=dd[i];                 \
 477     for(j=0; (j<NDIM); j++) tvec[j]=ev[j][i];\
 478     dd[i]=dd[i+1];                      \
 479     for(j=0; (j<NDIM); j++) ev[j][i]=ev[j][i+1];                \
 480     dd[i+1]=temp;                       \
 481     for(j=0; (j<NDIM); j++) ev[j][i+1]=tvec[j];                 \
 482   }
 483   SWAPPER(0)
 484   SWAPPER(1)
 485   SWAPPER(0)
 486
 487   for(i=0; (i<DIM); i++) {
 488     d[i]=dd[i];
 489     for(m=0; (m<DIM); m++)
 490       trans[i][m]=ev[m][i];
 491   }
 492
 493   for(i=0; (i<NDIM); i++) {
 494     sfree(inten[i]);
 495     sfree(ev[i]);
 496   }
 497   sfree(inten);
 498   sfree(ev);
 499 }
 500
 501
 502 /* WARNING WARNING WARNING
 503  * This routine rotates a molecule (I have checked this for water, PvM)
 504  * in the standard orientation used for water by researchers in the field.
 505  * This is different from the orientation used by Gray and Gubbins,
 506  * so be careful, with molecules other than water */
 507 void rot_mol_to_std_orient(int k0,int k1,atom_id index[],t_atom atom[],
 508                            rvec x[],matrix trans)
 509 {
 510   int  i;
 511   rvec xcm,xcq,d;
 512   matrix r_mat;
 513
 514   clear_rvec(xcm);
 515
 516   /* Compute the center of mass of the molecule and make it the origin */
 517   calc_xcm_mol(k0,k1,index,atom,x,xcm);
 518
 519   /* Compute the inertia moment tensor of a molecule */
 520   principal_comp_mol(k0,k1,index,atom,x,trans,d);
 521
 522   /* Rotate molecule to align with principal axes */
 523   rotate_mol(k0,k1,index,x,trans);
 524
 525
 526   /* If one of the moments is zero and the other two are equal, the
 527    * molecule is linear
 528    */
 529
 530   if (test_linear_mol(d)) {
 531     fprintf(stdout,"This molecule is linear\n");
 532   } else {
 533     fprintf(stdout,"This molecule is not linear\n");
 534
 535 make_rot_mat(ZZ,-0.5*M_PI,r_mat);
 536 rotate_mol(k0,k1,index,x,r_mat);
 537
 538
 539
 540     /* Now check if the center of charge now lies on the Z-axis
 541      * If not, rotate molecule so that it does.
 542      */
 543     for(i=0; (i<DIM); i++) {
 544       xcq[i]=test_qmom3(k0,k1,index,atom,x,i);
 545     }
 546
 547     if ((fabs(xcq[ZZ]) - TOLERANCE) < 0.0) {
 548       xcq[ZZ]=0.0;
 549     } else {
 550 #ifdef DEBUG
 551       fprintf(stdout,"Center of charge on Z-axis: %f\n",xcq[ZZ]);
 552 #endif
 553       if (xcq[ZZ] > 0.0) {
 554         make_rot_mat(XX,M_PI,r_mat);
 555         rotate_mol(k0,k1,index,x,r_mat);
 556       }
 557     }
 558
 559     if ((fabs(xcq[XX]) - TOLERANCE) < 0.0) {
 560       xcq[XX]=0.0;
 561     } else {
 562 #ifdef DEBUG
 563       fprintf(stdout,"Center of charge on X-axis: %f\n",xcq[XX]);
 564 #endif
 565       if (xcq[XX] < 0.0) {
 566         make_rot_mat(YY,0.5*M_PI,r_mat);
 567         rotate_mol(k0,k1,index,x,r_mat);
 568       } else {
 569         make_rot_mat(YY,-0.5*M_PI,r_mat);
 570         rotate_mol(k0,k1,index,x,r_mat);
 571       }
 572     }
 573
 574     if ((fabs(xcq[YY]) - TOLERANCE) < 0.0) {
 575       xcq[YY]=0.0;
 576     } else {
 577 #ifdef DEBUG
 578       fprintf(stdout,"Center of charge on Y-axis: %f\n",xcq[YY]);
 579 #endif
 580       if (xcq[YY] < 0.0) {
 581         make_rot_mat(XX,-0.5*M_PI,r_mat);
 582         rotate_mol(k0,k1,index,x,r_mat);
 583       } else {
 584         make_rot_mat(XX,0.5*M_PI,r_mat);
 585         rotate_mol(k0,k1,index,x,r_mat);
 586       }
 587     }
 588
 589     /* Now check the trace of the inertia tensor.
 590      * I want the water molecule in the YZ-plane */
 591     for(i=0; (i<DIM); i++) {
 592       xcm[i]=test_mmom2(k0,k1,index,atom,x,i);
 593     }
 594 #ifdef DEBUG
 595     fprintf(stdout,"xcm: %f %f %f\n",xcm[XX],xcm[YY],xcm[ZZ]);
 596 #endif
 597
 598     /* Check if X-component of inertia tensor is zero, if not
 599      * rotate molecule
 600      * This probably only works for water!!! PvM
 601      */
 602     if ((xcm[XX] - TOLERANCE) > 0.0) {
 603       make_rot_mat(ZZ,-0.5*M_PI,r_mat);
 604       rotate_mol(k0,k1,index,x,r_mat);
 605     }
 606   }
 607 }
 608
 609 /* Does the real work */
 610 void do_multipoles(char *trjfn,char *topfn,char *molndxfn,bool bFull)
 611 {
 612   int        i;
 613   int        gnx;
 614   atom_id    *grpindex;
 615   char       *grpname;
 616   t_topology *top;
 617   t_atom     *atom;
 618   t_block    *mols;
 619   int        natoms,status;
 620   real       t;
 621   matrix     box;
 622   real       t0,t1,tq;
 623   int        teller;
 624   bool       bCont;
 625
 626   rvec       *x,*m1;
 627   tensor     *m2;
 628   tensor3    *m3;
 629   tensor4    *m4;
 630   matrix     trans;
 631
 632   top  = read_top(topfn);
 633   rd_index(molndxfn,1,&gnx,&grpindex,&grpname);
 634   atom = top->atoms.atom;
 635   mols = &(top->blocks[ebMOLS]);
 636
 637   natoms  = read_first_x(&status,trjfn,&t,&x,box);
 638   snew(m1,gnx);
 639   snew(m2,gnx);
 640   snew(m3,gnx);
 641   snew(m4,gnx);
 642
 643   /* Start while loop over frames */
 644   do {
 645     /* PvM, bug in rm_pbc??? Does not work for virtual sites ...
 646     rm_pbc(&(top->idef),top->atoms.nr,box,x,x_s); */
 647
 648
 649     /* Begin loop of all molecules in index file */
 650     for(i=0; (i<gnx); i++) {
 651       int gi = grpindex[i];
 652
 653       rot_mol_to_std_orient(mols->index[gi],mols->index[gi+1],mols->a,atom,x,
 654                             trans);
 655
 656       /* Rotate the molecule along the principal moments axes */
 657       /* rotate_mol(mols->index[gi],mols->index[gi+1],mols->a,x,trans); */
 658
 659       /* Compute the multipole moments */
 660       mol_M1(mols->index[gi],mols->index[gi+1],mols->a,x,atom,m1[i]);
 661       mol_M2(mols->index[gi],mols->index[gi+1],mols->a,x,atom,m2[i]);
 662       mol_M3(mols->index[gi],mols->index[gi+1],mols->a,x,atom,m3[i]);
 663       mol_M4(mols->index[gi],mols->index[gi+1],mols->a,x,atom,m4[i]);
 664
 665       /* Spit it out */
 666       pr_M1(stdout,"",i,m1[i],t);
 667       pr_M2(stdout,"",m2[i],bFull);
 668       pr_M3(stdout,"",m3[i],bFull);
 669       pr_M4(stdout,"",m4[i],bFull);
 670
 671
 672     } /* End loop of all molecules in index file */
 673
 674     bCont = read_next_x(status,&t,natoms,x,box);
 675   } while(bCont);
 676
 677 }
 678
 679 int gmx_multipoles(int argc,char *argv[])
 680 {
 681   const char *desc[] = {
 682     "g_multipoles computes the electric multipole moments of",
 683     "molecules selected by a molecular index file.",
 684     "The center of mass of the molecule is used as the origin"
 685   };
 686
 687   static bool bFull = FALSE;
 688   static int  ntb=0;
 689   t_pargs pa[] = {
 690     { "-boxtype",FALSE,etINT,&ntb, "HIDDENbox type 0=rectangular; 1=truncated octahedron (only rectangular boxes are fully implemented)"},
 691     { "-full",   FALSE, etBOOL, &bFull,
 692       "Print all compononents of all multipoles instead of just the interesting ones" }
 693   };
 694   int          gnx;
 695   atom_id      *index;
 696   char         *grpname;
 697   t_filenm fnm[] = {
 698     { efTPX, NULL, NULL,      ffREAD },
 699     { efTRX, "-f", NULL,      ffREAD },
 700     { efNDX, NULL, NULL,      ffREAD }
 701   };
 702 #define NFILE asize(fnm)
 703   int     npargs;
 704   t_pargs *ppa;
 705
 706   int         i,j,k;
 707   int         natoms;
 708   int         step;
 709   real        t,lambda;
 710
 711   t_tpxheader tpx;
 712   t_topology  top;
 713   rvec        *x,*xnew;
 714   matrix      box;
 715   t_atom      *atom;
 716   rvec        dipole,dipole2;
 717   real        mtot,alfa,beta,gamma;
 718   rvec        CoM,*xCoM,angle;
 719   real        *xSqrCoM;
 720
 721   CopyRight(stderr,argv[0]);
 722
 723   parse_common_args(&argc,argv,PCA_CAN_TIME | PCA_BE_NICE,
 724                     NFILE,fnm,asize(pa),pa,asize(desc),desc,0,NULL);
 725
 726   do_multipoles(ftp2fn(efTRX,NFILE,fnm),ftp2fn(efTPX,NFILE,fnm),
 727                 ftp2fn(efNDX,NFILE,fnm),bFull);
 728
 729   return 0;
 730 }