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