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
 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
 234     if (dfile)
 235     {
 236         sprintf(buf, "Distance between %s and %s", grpn1, grpn2);
 237         sg_distance = xvgropen(dfile, buf, "Time (ps)", "Distance (nm)", oenv);
 238     }
 239
 240     if (d1file)
 241     {
 242         sprintf(buf, "Distance between plane and first atom of vector");
 243         sg_distance1 = xvgropen(d1file, buf, "Time (ps)", "Distance (nm)", oenv);
 244     }
 245
 246     if (d2file)
 247     {
 248         sprintf(buf, "Distance between plane and second atom of vector");
 249         sg_distance2 = xvgropen(d2file, buf, "Time (ps", "Distance (nm)", oenv);
 250     }
 251
 252     gpbc = gmx_rmpbc_init(&(top->idef), ePBC, natoms, box);
 253
 254     do
 255     {
 256         teller++;
 257
 258         gmx_rmpbc(gpbc, natoms, box, x0);
 259
 260         calc_angle(ePBC, box, x0, index1, index2, gnx1, gnx2, &angle,
 261                    &distance, &distance1, &distance2);
 262
 263         fprintf(sg_angle, "%12g  %12g  %12g\n", t, angle, acos(angle)*180.0/M_PI);
 264         if (dfile)
 265         {
 266             fprintf(sg_distance, "%12g  %12g\n", t, distance);
 267         }
 268         if (d1file)
 269         {
 270             fprintf(sg_distance1, "%12g  %12g\n", t, distance1);
 271         }
 272         if (d2file)
 273         {
 274             fprintf(sg_distance2, "%12g  %12g\n", t, distance1);
 275         }
 276
 277     }
 278     while (read_next_x(oenv, status, &t, natoms, x0, box));
 279
 280     gmx_rmpbc_done(gpbc);
 281
 282     fprintf(stderr, "\n");
 283     close_trj(status);
 284     ffclose(sg_angle);
 285     if (dfile)
 286     {
 287         ffclose(sg_distance);
 288     }
 289     if (d1file)
 290     {
 291         ffclose(sg_distance1);
 292     }
 293     if (d2file)
 294     {
 295         ffclose(sg_distance2);
 296     }
 297
 298     sfree(x0);
 299 }
 300
 301 static void calc_angle_single(int     ePBC,
 302                               matrix  box,
 303                               rvec    xzero[],
 304                               rvec    x[],
 305                               atom_id index1[],
 306                               atom_id index2[],
 307                               int     gnx1,
 308                               int     gnx2,
 309                               real   *angle,
 310                               real   *distance,
 311                               real   *distance1,
 312                               real   *distance2)
 313 {
 314     t_pbc pbc;
 315
 316     /* distance is distance between centers, distance 1 between center of plane
 317        and atom one of vector, distance 2 same for atom two
 318      */
 319
 320     rvec  normal1, normal2; /* normals on planes of interest */
 321     rvec  center1, center2;
 322     /* center of triangle of pts to define plane,
 323      * or center of vector if a vector is given
 324      */
 325     rvec  h1, h2, h3, h4, h5; /* temp. vectors */
 326
 327     if (box)
 328     {
 329         set_pbc(&pbc, ePBC, box);
 330     }
 331
 332     switch (gnx1)
 333     {
 334         case 3:     /* group 1 defines plane */
 335             calculate_normal(index1, xzero, normal1, center1);
 336             break;
 337         case 2:     /* group 1 defines vector */
 338             rvec_sub(xzero[index1[0]], xzero[index1[1]], normal1);
 339             rvec_add(xzero[index1[0]], xzero[index1[1]], h1);
 340             svmul(0.5, h1, center1); /* center is geometric mean */
 341             break;
 342         default:                     /* group 1 does none of the above */
 343             gmx_fatal(FARGS, "Something wrong with contents of index file.\n");
 344     }
 345
 346     switch (gnx2)
 347     {
 348         case 3:    /* group 2 defines plane */
 349             calculate_normal(index2, x, normal2, center2);
 350             break;
 351         case 2:    /* group 2 defines vector */
 352             rvec_sub(x[index2[0]], x[index2[1]], normal2);
 353             rvec_add(x[index2[0]], x[index2[1]], h2);
 354             svmul(0.5, h2, center2); /* center is geometric mean */
 355             break;
 356         default:                     /* group 2 does none of the above */
 357             gmx_fatal(FARGS, "Something wrong with contents of index file.\n");
 358     }
 359
 360     *angle = cos_angle(normal1, normal2);
 361
 362     if (box)
 363     {
 364         pbc_dx(&pbc, center1, center2, h3);
 365     }
 366     else
 367     {
 368         rvec_sub(center1, center2, h3);
 369     }
 370     *distance = norm(h3);
 371
 372     if (gnx1 == 3 && gnx2 == 2)
 373     {
 374         if (box)
 375         {
 376             pbc_dx(&pbc, center1, x[index2[0]], h4);
 377             pbc_dx(&pbc, center1, x[index2[1]], h5);
 378         }
 379         else
 380         {
 381             rvec_sub(center1, x[index2[0]], h4);
 382             rvec_sub(center1, x[index2[1]], h5);
 383         }
 384         *distance1 = norm(h4);
 385         *distance2 = norm(h5);
 386     }
 387     else if (gnx1 == 2 && gnx2 == 3)
 388     {
 389         rvec_sub(center1, xzero[index1[0]], h4);
 390         rvec_sub(center1, xzero[index1[1]], h5);
 391         *distance1 = norm(h4);
 392         *distance2 = norm(h5);
 393     }
 394     else
 395     {
 396         *distance1 = 0; *distance2 = 0;
 397     }
 398 }
 399
 400
 401 void sgangle_plot_single(const char *fn, const char *afile, const char *dfile,
 402                          const char *d1file, const char *d2file,
 403                          atom_id index1[], int gnx1, char *grpn1,
 404                          atom_id index2[], int gnx2, char *grpn2,
 405                          t_topology *top, int ePBC, const output_env_t oenv)
 406 {
 407     FILE
 408     *sg_angle,            /* xvgr file with angles */
 409     *sg_distance  = NULL, /* xvgr file with distances */
 410     *sg_distance1 = NULL, /* xvgr file with distance between plane and atom */
 411     *sg_distance2 = NULL; /* xvgr file with distance between plane and atom2 */
 412     real
 413         t,                /* time */
 414         angle,            /* cosine of angle between two groups */
 415         distance,         /* distance between two groups. */
 416         distance1,        /* distance between plane and one of two atoms */
 417         distance2;        /* same for second of two atoms */
 418     t_trxstatus *status;
 419     int          natoms, teller = 0;
 420     int          i;
 421     rvec        *x0; /* coordinates, and coordinates corrected for pb */
 422     rvec        *xzero;
 423     matrix       box;
 424     char         buf[256]; /* for xvgr title */
 425     gmx_rmpbc_t  gpbc = NULL;
 426
 427
 428     if ((natoms = read_first_x(oenv, &status, fn, &t, &x0, box)) == 0)
 429     {
 430         gmx_fatal(FARGS, "Could not read coordinates from statusfile\n");
 431     }
 432
 433     sprintf(buf, "Angle between %s and %s", grpn1, grpn2);
 434     sg_angle = xvgropen(afile, buf, "Time (ps)", "Cos(angle) ", oenv);
 435
 436     if (dfile)
 437     {
 438         sprintf(buf, "Distance between %s and %s", grpn1, grpn2);
 439         sg_distance = xvgropen(dfile, buf, "Time (ps)", "Distance (nm)", oenv);
 440     }
 441
 442     if (d1file)
 443     {
 444         sprintf(buf, "Distance between plane and first atom of vector");
 445         sg_distance1 = xvgropen(d1file, buf, "Time (ps)", "Distance (nm)", oenv);
 446     }
 447
 448     if (d2file)
 449     {
 450         sprintf(buf, "Distance between plane and second atom of vector");
 451         sg_distance2 = xvgropen(d2file, buf, "Time (ps", "Distance (nm)", oenv);
 452     }
 453
 454     snew(xzero, natoms);
 455     gpbc = gmx_rmpbc_init(&top->idef, ePBC, natoms, box);
 456
 457     do
 458     {
 459         teller++;
 460
 461         gmx_rmpbc(gpbc, natoms, box, x0);
 462         if (teller == 1)
 463         {
 464             for (i = 0; i < natoms; i++)
 465             {
 466                 copy_rvec(x0[i], xzero[i]);
 467             }
 468         }
 469
 470
 471         calc_angle_single(ePBC, box,
 472                           xzero, x0, index1, index2, gnx1, gnx2, &angle,
 473                           &distance, &distance1, &distance2);
 474
 475         fprintf(sg_angle, "%12g  %12g  %12g\n", t, angle, acos(angle)*180.0/M_PI);
 476         if (dfile)
 477         {
 478             fprintf(sg_distance, "%12g  %12g\n", t, distance);
 479         }
 480         if (d1file)
 481         {
 482             fprintf(sg_distance1, "%12g  %12g\n", t, distance1);
 483         }
 484         if (d2file)
 485         {
 486             fprintf(sg_distance2, "%12g  %12g\n", t, distance1);
 487         }
 488
 489     }
 490     while (read_next_x(oenv, status, &t, natoms, x0, box));
 491     gmx_rmpbc_done(gpbc);
 492
 493     fprintf(stderr, "\n");
 494     close_trj(status);
 495     ffclose(sg_angle);
 496     if (dfile)
 497     {
 498         ffclose(sg_distance);
 499     }
 500     if (d1file)
 501     {
 502         ffclose(sg_distance1);
 503     }
 504     if (d2file)
 505     {
 506         ffclose(sg_distance2);
 507     }
 508
 509     sfree(x0);
 510 }
 511
 512
 513
 514 int gmx_sgangle(int argc, char *argv[])
 515 {
 516     const char     *desc[] = {
 517         "Compute the angle and distance between two groups. ",
 518         "The groups are defined by a number of atoms given in an index file and",
 519         "may be two or three atoms in size.",
 520         "If [TT]-one[tt] is set, only one group should be specified in the index",
 521         "file and the angle between this group at time 0 and t will be computed.",
 522         "The angles calculated depend on the order in which the atoms are ",
 523         "given. Giving, for instance, 5 6 will rotate the vector 5-6 with ",
 524         "180 degrees compared to giving 6 5. [PAR]If three atoms are given, ",
 525         "the normal on the plane spanned by those three atoms will be",
 526         "calculated, using the formula  P1P2 x P1P3.",
 527         "The cos of the angle is calculated, using the inproduct of the two",
 528         "normalized vectors.[PAR]",
 529         "Here is what some of the file options do:[BR]",
 530         "[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]",
 531         "[TT]-od[tt]: Distance between two groups. Distance is taken from the center of one group to the center of the other group.[BR]",
 532         "[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]",
 533         "[TT]-od2[tt]: For two planes this option has no meaning."
 534     };
 535
 536     output_env_t    oenv;
 537     const char     *fna, *fnd, *fnd1, *fnd2;
 538     char     *      grpname[2];         /* name of the two groups */
 539     int             gnx[2];             /* size of the two groups */
 540     t_topology     *top;                /* topology         */
 541     int             ePBC;
 542     atom_id        *index[2];
 543     static gmx_bool bOne = FALSE, bZ = FALSE;
 544     t_pargs         pa[] = {
 545         { "-one", FALSE, etBOOL, {&bOne},
 546           "Only one group compute angle between vector at time zero and time t"},
 547         { "-z", FALSE, etBOOL, {&bZ},
 548           "Use the [IT]z[it]-axis as reference" }
 549     };
 550 #define NPA asize(pa)
 551
 552     t_filenm  fnm[] = {                         /* files for g_sgangle  */
 553         { efTRX, "-f", NULL,  ffREAD },         /* trajectory file  */
 554         { efNDX, NULL, NULL,  ffREAD },         /* index file       */
 555         { efTPX, NULL, NULL,  ffREAD },         /* topology file    */
 556         { efXVG, "-oa", "sg_angle", ffWRITE },  /* xvgr output file     */
 557         { efXVG, "-od", "sg_dist", ffOPTWR },   /* xvgr output file     */
 558         { efXVG, "-od1", "sg_dist1", ffOPTWR }, /* xvgr output file     */
 559         { efXVG, "-od2", "sg_dist2", ffOPTWR }  /* xvgr output file     */
 560     };
 561
 562 #define NFILE asize(fnm)
 563
 564     CopyRight(stderr, argv[0]);
 565     parse_common_args(&argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME | PCA_BE_NICE,
 566                       NFILE, fnm, NPA, pa, asize(desc), desc, 0, NULL, &oenv);
 567
 568
 569     top = read_top(ftp2fn(efTPX, NFILE, fnm), &ePBC); /* read topology file */
 570
 571     fna  = opt2fn("-oa", NFILE, fnm);
 572     fnd  = opt2fn_null("-od", NFILE, fnm);
 573     fnd1 = opt2fn_null("-od1", NFILE, fnm);
 574     fnd2 = opt2fn_null("-od2", NFILE, fnm);
 575
 576     /* read index file. */
 577     if (bOne)
 578     {
 579         rd_index(ftp2fn(efNDX, NFILE, fnm), 1, gnx, index, grpname);
 580         print_types(index[0], gnx[0], grpname[0],
 581                     index[0], gnx[0], grpname[0], top);
 582
 583         sgangle_plot_single(ftp2fn(efTRX, NFILE, fnm), fna, fnd, fnd1, fnd2,
 584                             index[0], gnx[0], grpname[0],
 585                             index[0], gnx[0], grpname[0],
 586                             top, ePBC, oenv);
 587     }
 588     else
 589     {
 590         rd_index(ftp2fn(efNDX, NFILE, fnm), bZ ? 1 : 2, gnx, index, grpname);
 591         if (!bZ)
 592         {
 593             print_types(index[0], gnx[0], grpname[0],
 594                         index[1], gnx[1], grpname[1], top);
 595         }
 596         else
 597         {
 598             gnx[1]     = 0;
 599             grpname[1] = strdup("Z-axis");
 600         }
 601         sgangle_plot(ftp2fn(efTRX, NFILE, fnm), fna, fnd, fnd1, fnd2,
 602                      index[0], gnx[0], grpname[0],
 603                      index[1], gnx[1], grpname[1],
 604                      top, ePBC, oenv);
 605     }
 606
 607     do_view(oenv, fna, "-nxy");  /* view xvgr file */
 608     do_view(oenv, fnd, "-nxy");  /* view xvgr file */
 609     do_view(oenv, fnd1, "-nxy"); /* view xvgr file */
 610     do_view(oenv, fnd2, "-nxy"); /* view xvgr file */
 611
 612     thanx(stderr);
 613     return 0;
 614 }