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