src/tools/gmx_sgangle.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,2013, 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 #ifdef HAVE_CONFIG_H
  39 #include <config.h>
  40 #endif
  41 #include <math.h>
  42
  43 #include "sysstuff.h"
  44 #include <string.h>
  45 #include "typedefs.h"
  46 #include "smalloc.h"
  47 #include "macros.h"
  48 #include "rmpbc.h"
  49 #include "vec.h"
  50 #include "xvgr.h"
  51 #include "pbc.h"
  52 #include "copyrite.h"
  53 #include "futil.h"
  54 #include "statutil.h"
  55 #include "index.h"
  56 #include "gmx_ana.h"
  57
  58
  59 /* this version only works correctly if one of the entries in the index file
  60    is a plane (three atoms specified) and the other a vector. Distance
  61    is calculated from the center of the plane to both atoms of the vector */
  62
  63 static void print_types(atom_id index1[], int gnx1, char *group1,
  64                         atom_id index2[], int gnx2, char *group2,
  65                         t_topology *top)
  66 {
  67     int i, j;
  68
  69     fprintf(stderr, "\n");
  70     fprintf(stderr, "Group %s contains the following atoms: \n", group1);
  71     for (i = 0; i < gnx1; i++)
  72     {
  73         fprintf(stderr, "Atomname %d: %s\n", i, *(top->atoms.atomname[index1[i]]));
  74     }
  75     fprintf(stderr, "\n");
  76
  77     fprintf(stderr, "Group %s contains the following atoms: \n", group2);
  78     for (j = 0; j < gnx2; j++)
  79     {
  80         fprintf(stderr, "Atomname %d: %s\n", j, *(top->atoms.atomname[index2[j]]));
  81     }
  82     fprintf(stderr, "\n");
  83
  84     fprintf(stderr, "Careful: distance only makes sense in some situations.\n\n");
  85 }
  86
  87 static void calculate_normal(atom_id index[], rvec x[], rvec result, rvec center)
  88 {
  89     rvec c1, c2;
  90     int  i;
  91
  92     /* calculate centroid of triangle spanned by the three points */
  93     for (i = 0; i < 3; i++)
  94     {
  95         center[i] = (x[index[0]][i] + x[index[1]][i] + x[index[2]][i])/3;
  96     }
  97
  98     /* use P1P2 x P1P3 to calculate normal, given three points P1-P3 */
  99     rvec_sub(x[index[1]], x[index[0]], c1); /* find two vectors */
 100     rvec_sub(x[index[2]], x[index[0]], c2);
 101
 102     cprod(c1, c2, result);                /* take crossproduct between these */
 103 }
 104
 105 /* calculate the angle and distance between the two groups */
 106 static void calc_angle(int ePBC, matrix box, rvec x[], atom_id index1[],
 107                        atom_id index2[], int gnx1, int gnx2,
 108                        real *angle,      real *distance,
 109                        real *distance1,  real *distance2)
 110
 111 /* distance is distance between centers, distance 1 between center of plane
 112    and atom one of vector, distance 2 same for atom two
 113  */
 114
 115 {
 116     rvec
 117         normal1, normal2,   /* normals on planes of interest */
 118         center1, center2,   /* center of triangle of points given to define plane,*/
 119                             /* or center of vector if a vector is given */
 120         h1, h2, h3, h4, h5; /* temp. vectors */
 121     t_pbc pbc;
 122
 123     set_pbc(&pbc, ePBC, box);
 124
 125     switch (gnx1)
 126     {
 127         case 3:       /* group 1 defines plane */
 128             calculate_normal(index1, x, normal1, center1);
 129             break;
 130         case 2:       /* group 1 defines vector */
 131             rvec_sub(x[index1[0]], x[index1[1]], normal1);
 132             rvec_add(x[index1[0]], x[index1[1]], h1);
 133             svmul(0.5, h1, center1); /* center is geometric mean */
 134             break;
 135         default:                     /* group 1 does none of the above */
 136             gmx_fatal(FARGS, "Something wrong with contents of index file. Groups should contain 2 or 3 atoms.\n");
 137     }
 138
 139     switch (gnx2)
 140     {
 141         case 3:      /* group 2 defines plane */
 142             calculate_normal(index2, x, normal2, center2);
 143             break;
 144         case 2:      /* group 2 defines vector */
 145             rvec_sub(x[index2[0]], x[index2[1]], normal2);
 146             rvec_add(x[index2[0]], x[index2[1]], h2);
 147             svmul(0.5, h2, center2); /* center is geometric mean */
 148             break;
 149         case 0:
 150             normal2[XX] = 0;
 151             normal2[YY] = 0;
 152             normal2[ZZ] = 1;
 153             center2[XX] = 0;
 154             center2[YY] = 0;
 155             center2[ZZ] = 0;
 156             break;
 157         default:     /* group 2 does none of the above */
 158             gmx_fatal(FARGS, "Something wrong with contents of index file.\n");
 159     }
 160
 161     *angle = cos_angle(normal1, normal2);
 162
 163     if (box)
 164     {
 165         pbc_dx(&pbc, center1, center2, h3);
 166     }
 167     else
 168     {
 169         rvec_sub(center1, center2, h3);
 170     }
 171     *distance = norm(h3);
 172
 173     if (gnx1 == 3 && gnx2 == 2)
 174     {
 175         if (box)
 176         {
 177             pbc_dx(&pbc, center1, x[index2[0]], h4);
 178             pbc_dx(&pbc, center1, x[index2[1]], h5);
 179         }
 180         else
 181         {
 182             rvec_sub(center1, x[index2[0]], h4);
 183             rvec_sub(center1, x[index2[1]], h5);
 184         }
 185         *distance1 = norm(h4);
 186         *distance2 = norm(h5);
 187     }
 188     else if (gnx1 == 2 && gnx2 == 3)
 189     {
 190         rvec_sub(center1, x[index1[0]], h4);
 191         rvec_sub(center1, x[index1[1]], h5);
 192         *distance1 = norm(h4);
 193         *distance2 = norm(h5);
 194     }
 195     else
 196     {
 197         *distance1 = 0; *distance2 = 0;
 198     }
 199 }
 200
 201 void sgangle_plot(const char *fn, const char *afile, const char *dfile,
 202                   const char *d1file, const char *d2file,
 203                   atom_id index1[], int gnx1, char *grpn1,
 204                   atom_id index2[], int gnx2, char *grpn2,
 205                   t_topology *top, int ePBC, const output_env_t oenv)
 206 {
 207     FILE
 208     *sg_angle,            /* xvgr file with angles */
 209     *sg_distance  = NULL, /* xvgr file with distances */
 210     *sg_distance1 = NULL, /* xvgr file with distance between plane and atom */
 211     *sg_distance2 = NULL; /* xvgr file with distance between plane and atom2 */
 212     real
 213         t,                /* time */
 214         angle,            /* cosine of angle between two groups */
 215         distance,         /* distance between two groups. */
 216         distance1,        /* distance between plane and one of two atoms */
 217         distance2;        /* same for second of two atoms */
 218     t_trxstatus *status;
 219     int          natoms, teller = 0;
 220     rvec        *x0;       /* coordinates, and coordinates corrected for pb */
 221     matrix       box;
 222     char         buf[256]; /* for xvgr title */
 223     gmx_rmpbc_t  gpbc = NULL;
 224     const char*  aleg[2] = { "cos(Angle)", "Angle (degrees)" };     /* legends for sg_angle output file */
 225
 226     if ((natoms = read_first_x(oenv, &status, fn, &t, &x0, box)) == 0)
 227     {
 228         gmx_fatal(FARGS, "Could not read coordinates from statusfile\n");
 229     }
 230
 231     sprintf(buf, "Angle between %s and %s", grpn1, grpn2);
 232     sg_angle = xvgropen(afile, buf, "Time (ps)", "Angle (degrees)", oenv);
 233     xvgr_legend(sg_angle, 2, aleg, oenv);
 234
 235     if (dfile)
 236     {
 237         sprintf(buf, "Distance between %s and %s", grpn1, grpn2);
 238         sg_distance = xvgropen(dfile, buf, "Time (ps)", "Distance (nm)", oenv);
 239     }
 240
 241     if (d1file)
 242     {
 243         sprintf(buf, "Distance between plane and first atom of vector");
 244         sg_distance1 = xvgropen(d1file, buf, "Time (ps)", "Distance (nm)", oenv);
 245     }
 246
 247     if (d2file)
 248     {
 249         sprintf(buf, "Distance between plane and second atom of vector");
 250         sg_distance2 = xvgropen(d2file, buf, "Time (ps", "Distance (nm)", oenv);
 251     }
 252
 253     gpbc = gmx_rmpbc_init(&(top->idef), ePBC, natoms, box);
 254
 255     do
 256     {
 257         teller++;
 258
 259         gmx_rmpbc(gpbc, natoms, box, x0);
 260
 261         calc_angle(ePBC, box, x0, index1, index2, gnx1, gnx2, &angle,
 262                    &distance, &distance1, &distance2);
 263
 264         fprintf(sg_angle, "%12g  %12g  %12g\n", t, angle, acos(angle)*180.0/M_PI);
 265         if (dfile)
 266         {
 267             fprintf(sg_distance, "%12g  %12g\n", t, distance);
 268         }
 269         if (d1file)
 270         {
 271             fprintf(sg_distance1, "%12g  %12g\n", t, distance1);
 272         }
 273         if (d2file)
 274         {
 275             fprintf(sg_distance2, "%12g  %12g\n", t, distance1);
 276         }
 277
 278     }
 279     while (read_next_x(oenv, status, &t, natoms, x0, box));
 280
 281     gmx_rmpbc_done(gpbc);
 282
 283     fprintf(stderr, "\n");
 284     close_trj(status);
 285     ffclose(sg_angle);
 286     if (dfile)
 287     {
 288         ffclose(sg_distance);
 289     }
 290     if (d1file)
 291     {
 292         ffclose(sg_distance1);
 293     }
 294     if (d2file)
 295     {
 296         ffclose(sg_distance2);
 297     }
 298
 299     sfree(x0);
 300 }
 301
 302 static void calc_angle_single(int     ePBC,
 303                               matrix  box,
 304                               rvec    xzero[],
 305                               rvec    x[],
 306                               atom_id index1[],
 307                               atom_id index2[],
 308                               int     gnx1,
 309                               int     gnx2,
 310                               real   *angle,
 311                               real   *distance,
 312                               real   *distance1,
 313                               real   *distance2)
 314 {
 315     t_pbc pbc;
 316
 317     /* distance is distance between centers, distance 1 between center of plane
 318        and atom one of vector, distance 2 same for atom two
 319      */
 320
 321     rvec  normal1, normal2; /* normals on planes of interest */
 322     rvec  center1, center2;
 323     /* center of triangle of pts to define plane,
 324      * or center of vector if a vector is given
 325      */
 326     rvec  h1, h2, h3, h4, h5; /* temp. vectors */
 327
 328     if (box)
 329     {
 330         set_pbc(&pbc, ePBC, box);
 331     }
 332
 333     switch (gnx1)
 334     {
 335         case 3:     /* group 1 defines plane */
 336             calculate_normal(index1, xzero, normal1, center1);
 337             break;
 338         case 2:     /* group 1 defines vector */
 339             rvec_sub(xzero[index1[0]], xzero[index1[1]], normal1);
 340             rvec_add(xzero[index1[0]], xzero[index1[1]], h1);
 341             svmul(0.5, h1, center1); /* center is geometric mean */
 342             break;
 343         default:                     /* group 1 does none of the above */
 344             gmx_fatal(FARGS, "Something wrong with contents of index file.\n");
 345     }
 346
 347     switch (gnx2)
 348     {
 349         case 3:    /* group 2 defines plane */
 350             calculate_normal(index2, x, normal2, center2);
 351             break;
 352         case 2:    /* group 2 defines vector */
 353             rvec_sub(x[index2[0]], x[index2[1]], normal2);
 354             rvec_add(x[index2[0]], x[index2[1]], h2);
 355             svmul(0.5, h2, center2); /* center is geometric mean */
 356             break;
 357         default:                     /* group 2 does none of the above */
 358             gmx_fatal(FARGS, "Something wrong with contents of index file.\n");
 359     }
 360
 361     *angle = cos_angle(normal1, normal2);
 362
 363     if (box)
 364     {
 365         pbc_dx(&pbc, center1, center2, h3);
 366     }
 367     else
 368     {
 369         rvec_sub(center1, center2, h3);
 370     }
 371     *distance = norm(h3);
 372
 373     if (gnx1 == 3 && gnx2 == 2)
 374     {
 375         if (box)
 376         {
 377             pbc_dx(&pbc, center1, x[index2[0]], h4);
 378             pbc_dx(&pbc, center1, x[index2[1]], h5);
 379         }
 380         else
 381         {
 382             rvec_sub(center1, x[index2[0]], h4);
 383             rvec_sub(center1, x[index2[1]], h5);
 384         }
 385         *distance1 = norm(h4);
 386         *distance2 = norm(h5);
 387     }
 388     else if (gnx1 == 2 && gnx2 == 3)
 389     {
 390         rvec_sub(center1, xzero[index1[0]], h4);
 391         rvec_sub(center1, xzero[index1[1]], h5);
 392         *distance1 = norm(h4);
 393         *distance2 = norm(h5);
 394     }
 395     else
 396     {
 397         *distance1 = 0; *distance2 = 0;
 398     }
 399 }
 400
 401
 402 void sgangle_plot_single(const char *fn, const char *afile, const char *dfile,
 403                          const char *d1file, const char *d2file,
 404                          atom_id index1[], int gnx1, char *grpn1,
 405                          atom_id index2[], int gnx2, char *grpn2,
 406                          t_topology *top, int ePBC, const output_env_t oenv)
 407 {
 408     FILE
 409     *sg_angle,            /* xvgr file with angles */
 410     *sg_distance  = NULL, /* xvgr file with distances */
 411     *sg_distance1 = NULL, /* xvgr file with distance between plane and atom */
 412     *sg_distance2 = NULL; /* xvgr file with distance between plane and atom2 */
 413     real
 414         t,                /* time */
 415         angle,            /* cosine of angle between two groups */
 416         distance,         /* distance between two groups. */
 417         distance1,        /* distance between plane and one of two atoms */
 418         distance2;        /* same for second of two atoms */
 419     t_trxstatus *status;
 420     int          natoms, teller = 0;
 421     int          i;
 422     rvec        *x0; /* coordinates, and coordinates corrected for pb */
 423     rvec        *xzero;
 424     matrix       box;
 425     char         buf[256]; /* for xvgr title */
 426     gmx_rmpbc_t  gpbc = NULL;
 427
 428
 429     if ((natoms = read_first_x(oenv, &status, fn, &t, &x0, box)) == 0)
 430     {
 431         gmx_fatal(FARGS, "Could not read coordinates from statusfile\n");
 432     }
 433
 434     sprintf(buf, "Angle between %s and %s", grpn1, grpn2);
 435     sg_angle = xvgropen(afile, buf, "Time (ps)", "Cos(angle) ", oenv);
 436
 437     if (dfile)
 438     {
 439         sprintf(buf, "Distance between %s and %s", grpn1, grpn2);
 440         sg_distance = xvgropen(dfile, buf, "Time (ps)", "Distance (nm)", oenv);
 441     }
 442
 443     if (d1file)
 444     {
 445         sprintf(buf, "Distance between plane and first atom of vector");
 446         sg_distance1 = xvgropen(d1file, buf, "Time (ps)", "Distance (nm)", oenv);
 447     }
 448
 449     if (d2file)
 450     {
 451         sprintf(buf, "Distance between plane and second atom of vector");
 452         sg_distance2 = xvgropen(d2file, buf, "Time (ps", "Distance (nm)", oenv);
 453     }
 454
 455     snew(xzero, natoms);
 456     gpbc = gmx_rmpbc_init(&top->idef, ePBC, natoms, box);
 457
 458     do
 459     {
 460         teller++;
 461
 462         gmx_rmpbc(gpbc, natoms, box, x0);
 463         if (teller == 1)
 464         {
 465             for (i = 0; i < natoms; i++)
 466             {
 467                 copy_rvec(x0[i], xzero[i]);
 468             }
 469         }
 470
 471
 472         calc_angle_single(ePBC, box,
 473                           xzero, x0, index1, index2, gnx1, gnx2, &angle,
 474                           &distance, &distance1, &distance2);
 475
 476         fprintf(sg_angle, "%12g  %12g  %12g\n", t, angle, acos(angle)*180.0/M_PI);
 477         if (dfile)
 478         {
 479             fprintf(sg_distance, "%12g  %12g\n", t, distance);
 480         }
 481         if (d1file)
 482         {
 483             fprintf(sg_distance1, "%12g  %12g\n", t, distance1);
 484         }
 485         if (d2file)
 486         {
 487             fprintf(sg_distance2, "%12g  %12g\n", t, distance1);
 488         }
 489
 490     }
 491     while (read_next_x(oenv, status, &t, natoms, x0, box));
 492     gmx_rmpbc_done(gpbc);
 493
 494     fprintf(stderr, "\n");
 495     close_trj(status);
 496     ffclose(sg_angle);
 497     if (dfile)
 498     {
 499         ffclose(sg_distance);
 500     }
 501     if (d1file)
 502     {
 503         ffclose(sg_distance1);
 504     }
 505     if (d2file)
 506     {
 507         ffclose(sg_distance2);
 508     }
 509
 510     sfree(x0);
 511 }
 512
 513
 514
 515 int gmx_sgangle(int argc, char *argv[])
 516 {
 517     const char     *desc[] = {
 518         "Compute the angle and distance between two groups. ",
 519         "The groups are defined by a number of atoms given in an index file and",
 520         "may be two or three atoms in size.",
 521         "If [TT]-one[tt] is set, only one group should be specified in the index",
 522         "file and the angle between this group at time 0 and t will be computed.",
 523         "The angles calculated depend on the order in which the atoms are ",
 524         "given. Giving, for instance, 5 6 will rotate the vector 5-6 with ",
 525         "180 degrees compared to giving 6 5. [PAR]If three atoms are given, ",
 526         "the normal on the plane spanned by those three atoms will be",
 527         "calculated, using the formula  P1P2 x P1P3.",
 528         "The cos of the angle is calculated, using the inproduct of the two",
 529         "normalized vectors.[PAR]",
 530         "Here is what some of the file options do:[BR]",
 531         "[TT]-oa[tt]: Angle between the two groups specified in the index file. If a group contains three atoms the normal to the plane defined by those three atoms will be used. If a group contains two atoms, the vector defined by those two atoms will be used.[BR]",
 532         "[TT]-od[tt]: Distance between two groups. Distance is taken from the center of one group to the center of the other group.[BR]",
 533         "[TT]-od1[tt]: If one plane and one vector is given, the distances for each of the atoms from the center of the plane is given separately.[BR]",
 534         "[TT]-od2[tt]: For two planes this option has no meaning."
 535     };
 536
 537     output_env_t    oenv;
 538     const char     *fna, *fnd, *fnd1, *fnd2;
 539     char     *      grpname[2];         /* name of the two groups */
 540     int             gnx[2];             /* size of the two groups */
 541     t_topology     *top;                /* topology         */
 542     int             ePBC;
 543     atom_id        *index[2];
 544     static gmx_bool bOne = FALSE, bZ = FALSE;
 545     t_pargs         pa[] = {
 546         { "-one", FALSE, etBOOL, {&bOne},
 547           "Only one group compute angle between vector at time zero and time t"},
 548         { "-z", FALSE, etBOOL, {&bZ},
 549           "Use the [IT]z[it]-axis as reference" }
 550     };
 551 #define NPA asize(pa)
 552
 553     t_filenm  fnm[] = {                         /* files for g_sgangle  */
 554         { efTRX, "-f", NULL,  ffREAD },         /* trajectory file  */
 555         { efNDX, NULL, NULL,  ffREAD },         /* index file       */
 556         { efTPX, NULL, NULL,  ffREAD },         /* topology file    */
 557         { efXVG, "-oa", "sg_angle", ffWRITE },  /* xvgr output file     */
 558         { efXVG, "-od", "sg_dist", ffOPTWR },   /* xvgr output file     */
 559         { efXVG, "-od1", "sg_dist1", ffOPTWR }, /* xvgr output file     */
 560         { efXVG, "-od2", "sg_dist2", ffOPTWR }  /* xvgr output file     */
 561     };
 562
 563 #define NFILE asize(fnm)
 564
 565     CopyRight(stderr, argv[0]);
 566     parse_common_args(&argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME | PCA_BE_NICE,
 567                       NFILE, fnm, NPA, pa, asize(desc), desc, 0, NULL, &oenv);
 568
 569
 570     top = read_top(ftp2fn(efTPX, NFILE, fnm), &ePBC); /* read topology file */
 571
 572     fna  = opt2fn("-oa", NFILE, fnm);
 573     fnd  = opt2fn_null("-od", NFILE, fnm);
 574     fnd1 = opt2fn_null("-od1", NFILE, fnm);
 575     fnd2 = opt2fn_null("-od2", NFILE, fnm);
 576
 577     /* read index file. */
 578     if (bOne)
 579     {
 580         rd_index(ftp2fn(efNDX, NFILE, fnm), 1, gnx, index, grpname);
 581         print_types(index[0], gnx[0], grpname[0],
 582                     index[0], gnx[0], grpname[0], top);
 583
 584         sgangle_plot_single(ftp2fn(efTRX, NFILE, fnm), fna, fnd, fnd1, fnd2,
 585                             index[0], gnx[0], grpname[0],
 586                             index[0], gnx[0], grpname[0],
 587                             top, ePBC, oenv);
 588     }
 589     else
 590     {
 591         rd_index(ftp2fn(efNDX, NFILE, fnm), bZ ? 1 : 2, gnx, index, grpname);
 592         if (!bZ)
 593         {
 594             print_types(index[0], gnx[0], grpname[0],
 595                         index[1], gnx[1], grpname[1], top);
 596         }
 597         else
 598         {
 599             gnx[1]     = 0;
 600             grpname[1] = strdup("Z-axis");
 601         }
 602         sgangle_plot(ftp2fn(efTRX, NFILE, fnm), fna, fnd, fnd1, fnd2,
 603                      index[0], gnx[0], grpname[0],
 604                      index[1], gnx[1], grpname[1],
 605                      top, ePBC, oenv);
 606     }
 607
 608     do_view(oenv, fna, "-nxy");  /* view xvgr file */
 609     do_view(oenv, fnd, "-nxy");  /* view xvgr file */
 610     do_view(oenv, fnd1, "-nxy"); /* view xvgr file */
 611     do_view(oenv, fnd2, "-nxy"); /* view xvgr file */
 612
 613     thanx(stderr);
 614     return 0;
 615 }