3 * This source code is part of
7 * GROningen MAchine for Chemical Simulations
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.
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.
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.
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.
30 * For more info, check our website at http://www.gromacs.org
33 * Green Red Orange Magenta Azure Cyan Skyblue
55 #define TOLERANCE 1.0E-8
57 #define e2d(x) ENM2DEBYE*(x)
58 #define delta(a,b) (( a == b ) ? 1.0 : 0.0)
60 #define NDIM 3 /* We will be using a numerical recipes routine */
62 static char dim[DIM+1] = "XYZ";
64 typedef real tensor3[DIM][DIM][DIM]; /* 3 rank tensor */
65 typedef real tensor4[DIM][DIM][DIM][DIM]; /* 4 rank tensor */
68 void pr_coord(int k0,int k1,atom_id index[],rvec x[],char *msg)
72 fprintf(stdout,"Coordinates in nm (%s)\n",msg);
73 for(k=k0; (k<k1); k++) {
75 fprintf(stdout,"Atom %d, %15.10f %15.10f %15.10f\n",
76 kk,x[kk][XX],x[kk][YY],x[kk][ZZ]);
81 static void clear_tensor3(tensor3 a)
86 for(i=0; (i<DIM); i++)
87 for(j=0; (j<DIM); j++)
88 for(k=0; (k<DIM); k++)
92 static void clear_tensor4(tensor4 a)
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++)
104 void rotate_mol(int k0,int k1,atom_id index[],rvec x[],matrix trans)
109 for(k=k0; (k<k1); k++) {
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;
120 /* the following routines are heavily inspired by the Gaussian 94 source
125 Make the rotation matrix for angle Theta counterclockwise about
129 void make_rot_mat(int axis,real theta,matrix t_mat){
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;
156 gmx_bool test_linear_mol(rvec d)
158 /* d is sorted in descending order */
159 if ( (d[ZZ] < TOLERANCE) && (d[XX]-d[YY]) < TOLERANCE ) {
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){
172 for(k=k0; (k<k1); k++) {
175 xcq+=q*x[kk][axis]*x[kk][axis]*x[kk][axis];
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){
188 for(k=k0; (k<k1); k++) {
191 xcm+=m*x[kk][axis]*x[kk][axis];
197 real calc_xcm_mol(int k0,int k1,atom_id index[],t_atom atom[],rvec x[],
203 /* Compute the center of mass */
206 for(k=k0; (k<k1); k++) {
210 for(m=0; (m<DIM); m++)
213 for(m=0; (m<DIM); m++)
216 /* And make it the origin */
217 for(k=k0; (k<k1); k++) {
225 /* Retruns the center of charge */
226 real calc_xcq_mol(int k0,int k1,atom_id index[],t_atom atom[],rvec x[],
234 for(k=k0; (k<k1); k++) {
238 fprintf(stdout,"tq: %f, q0: %f\n",tq,q0);
239 for(m=0; (m<DIM); m++)
243 for(m=0; (m<DIM); m++)
246 for(k=k0; (k<k1); k++) {
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)
261 for(n=n0; (n<n1); n++) {
264 for(m=0; (m<DIM); m++)
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)
276 for(n=n0; (n<n1); n++) {
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;
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)
293 for(n=n0; (n<n1); n++) {
297 for(i=0; (i<DIM); i++)
298 for(j=0; (j<DIM); j++)
299 for(k=0; (k<DIM); 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;
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)
315 for(n=n0; (n<n1); n++) {
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++)
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;
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)
343 fprintf(fp,"Molecule: %d @ t= %f ps\n",mol,time);
345 m1_tot = sqrt(m1[XX]*m1[XX]+m1[YY]*m1[YY]+m1[ZZ]*m1[ZZ]);
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);
352 /* Print the quadrupole moment components */
353 void pr_M2(FILE *fp,char *msg,tensor m2,gmx_bool bFull)
357 fprintf(fp,"Quadrupole Moment %s(Debye-Ang):\n",msg);
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]);
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]);
373 /* Print the octopole moment components */
374 void pr_M3(FILE *fp,char *msg,tensor3 m3,gmx_bool bFull)
378 fprintf(fp,"Octopole Moment %s(Debye-Ang^2):\n",msg);
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]);
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]);
398 /* Print the hexadecapole moment components */
399 void pr_M4(FILE *fp,char *msg,tensor4 m4,gmx_bool bFull)
403 fprintf(fp,"Hexadecapole Moment %s(Debye-Ang^3):\n",msg);
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]);
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],
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[],
438 double **inten,dd[NDIM],tvec[NDIM],**ev;
443 for(i=0; (i<NDIM); i++) {
449 for(i=0; (i<NDIM); i++)
450 for(m=0; (m<NDIM); m++)
453 for(i=k0; (i<k1); i++) {
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);
466 inten[0][1]=inten[1][0];
467 inten[0][2]=inten[2][0];
468 inten[1][2]=inten[2][1];
470 /* Call numerical recipe routines */
471 jacobi(inten,3,dd,ev,&nrot);
473 /* Sort eigenvalues in descending order */
475 if (fabs(dd[i+1]) > fabs(dd[i])) { \
477 for(j=0; (j<NDIM); j++) tvec[j]=ev[j][i];\
479 for(j=0; (j<NDIM); j++) ev[j][i]=ev[j][i+1]; \
481 for(j=0; (j<NDIM); j++) ev[j][i+1]=tvec[j]; \
487 for(i=0; (i<DIM); i++) {
489 for(m=0; (m<DIM); m++)
490 trans[i][m]=ev[m][i];
493 for(i=0; (i<NDIM); i++) {
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)
516 /* Compute the center of mass of the molecule and make it the origin */
517 calc_xcm_mol(k0,k1,index,atom,x,xcm);
519 /* Compute the inertia moment tensor of a molecule */
520 principal_comp_mol(k0,k1,index,atom,x,trans,d);
522 /* Rotate molecule to align with principal axes */
523 rotate_mol(k0,k1,index,x,trans);
526 /* If one of the moments is zero and the other two are equal, the
530 if (test_linear_mol(d)) {
531 fprintf(stdout,"This molecule is linear\n");
533 fprintf(stdout,"This molecule is not linear\n");
535 make_rot_mat(ZZ,-0.5*M_PI,r_mat);
536 rotate_mol(k0,k1,index,x,r_mat);
540 /* Now check if the center of charge now lies on the Z-axis
541 * If not, rotate molecule so that it does.
543 for(i=0; (i<DIM); i++) {
544 xcq[i]=test_qmom3(k0,k1,index,atom,x,i);
547 if ((fabs(xcq[ZZ]) - TOLERANCE) < 0.0) {
551 fprintf(stdout,"Center of charge on Z-axis: %f\n",xcq[ZZ]);
554 make_rot_mat(XX,M_PI,r_mat);
555 rotate_mol(k0,k1,index,x,r_mat);
559 if ((fabs(xcq[XX]) - TOLERANCE) < 0.0) {
563 fprintf(stdout,"Center of charge on X-axis: %f\n",xcq[XX]);
566 make_rot_mat(YY,0.5*M_PI,r_mat);
567 rotate_mol(k0,k1,index,x,r_mat);
569 make_rot_mat(YY,-0.5*M_PI,r_mat);
570 rotate_mol(k0,k1,index,x,r_mat);
574 if ((fabs(xcq[YY]) - TOLERANCE) < 0.0) {
578 fprintf(stdout,"Center of charge on Y-axis: %f\n",xcq[YY]);
581 make_rot_mat(XX,-0.5*M_PI,r_mat);
582 rotate_mol(k0,k1,index,x,r_mat);
584 make_rot_mat(XX,0.5*M_PI,r_mat);
585 rotate_mol(k0,k1,index,x,r_mat);
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);
595 fprintf(stdout,"xcm: %f %f %f\n",xcm[XX],xcm[YY],xcm[ZZ]);
598 /* Check if X-component of inertia tensor is zero, if not
600 * This probably only works for water!!! PvM
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);
609 /* Does the real work */
610 void do_multipoles(char *trjfn,char *topfn,char *molndxfn,gmx_bool bFull)
631 gmx_rmpbc_t gpbc=NULL;
633 top = read_top(topfn);
634 rd_index(molndxfn,1,&gnx,&grpindex,&grpname);
635 atom = top->atoms.atom;
636 mols = &(top->blocks[ebMOLS]);
638 natoms = read_first_x(&status,trjfn,&t,&x,box);
644 gpbc = gmx_rmpbc_init(&top->idef,ePBC,top->atoms.nr,box);
646 /* Start while loop over frames */
648 /* PvM, bug in rm_pbc??? Does not work for virtual sites ...
649 gmx_rmpbc(gpbc,box,x,x_s); */
651 /* Begin loop of all molecules in index file */
652 for(i=0; (i<gnx); i++) {
653 int gi = grpindex[i];
655 rot_mol_to_std_orient(mols->index[gi],mols->index[gi+1],mols->a,atom,x,
658 /* Rotate the molecule along the principal moments axes */
659 /* rotate_mol(mols->index[gi],mols->index[gi+1],mols->a,x,trans); */
661 /* Compute the multipole moments */
662 mol_M1(mols->index[gi],mols->index[gi+1],mols->a,x,atom,m1[i]);
663 mol_M2(mols->index[gi],mols->index[gi+1],mols->a,x,atom,m2[i]);
664 mol_M3(mols->index[gi],mols->index[gi+1],mols->a,x,atom,m3[i]);
665 mol_M4(mols->index[gi],mols->index[gi+1],mols->a,x,atom,m4[i]);
668 pr_M1(stdout,"",i,m1[i],t);
669 pr_M2(stdout,"",m2[i],bFull);
670 pr_M3(stdout,"",m3[i],bFull);
671 pr_M4(stdout,"",m4[i],bFull);
674 } /* End loop of all molecules in index file */
676 bCont = read_next_x(status,&t,natoms,x,box);
678 gmx_rmpbc_done(gpbc);
684 int gmx_multipoles(int argc,char *argv[])
686 const char *desc[] = {
687 "g_multipoles computes the electric multipole moments of",
688 "molecules selected by a molecular index file.",
689 "The center of mass of the molecule is used as the origin"
692 static gmx_bool bFull = FALSE;
695 { "-boxtype",FALSE,etINT,&ntb, "HIDDENbox type 0=rectangular; 1=truncated octahedron (only rectangular boxes are fully implemented)"},
696 { "-full", FALSE, etBOOL, &bFull,
697 "Print all compononents of all multipoles instead of just the interesting ones" }
703 { efTPX, NULL, NULL, ffREAD },
704 { efTRX, "-f", NULL, ffREAD },
705 { efNDX, NULL, NULL, ffREAD }
707 #define NFILE asize(fnm)
722 real mtot,alfa,beta,gamma;
723 rvec CoM,*xCoM,angle;
726 CopyRight(stderr,argv[0]);
728 parse_common_args(&argc,argv,PCA_CAN_TIME | PCA_BE_NICE,
729 NFILE,fnm,asize(pa),pa,asize(desc),desc,0,NULL);
731 do_multipoles(ftp2fn(efTRX,NFILE,fnm),ftp2fn(efTPX,NFILE,fnm),
732 ftp2fn(efNDX,NFILE,fnm),bFull);