src/gromacs/gmxana/gmx_sorient.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
  39 #include "macros.h"
  40 #include "statutil.h"
  41 #include "smalloc.h"
  42 #include "copyrite.h"
  43 #include "gstat.h"
  44 #include "vec.h"
  45 #include "xvgr.h"
  46 #include "pbc.h"
  47 #include "index.h"
  48 #include "tpxio.h"
  49 #include "gmx_ana.h"
  50
  51
  52 static void calc_com_pbc(int nrefat, t_topology *top, rvec x[], t_pbc *pbc,
  53                          atom_id index[], rvec xref, gmx_bool bPBC, matrix box)
  54 {
  55     const real tol = 1e-4;
  56     gmx_bool   bChanged;
  57     int        m, j, ai, iter;
  58     real       mass, mtot;
  59     rvec       dx, xtest;
  60
  61     /* First simple calculation */
  62     clear_rvec(xref);
  63     mtot = 0;
  64     for (m = 0; (m < nrefat); m++)
  65     {
  66         ai   = index[m];
  67         mass = top->atoms.atom[ai].m;
  68         for (j = 0; (j < DIM); j++)
  69         {
  70             xref[j] += mass*x[ai][j];
  71         }
  72         mtot += mass;
  73     }
  74     svmul(1/mtot, xref, xref);
  75     /* Now check if any atom is more than half the box from the COM */
  76     if (bPBC)
  77     {
  78         iter = 0;
  79         do
  80         {
  81             bChanged = FALSE;
  82             for (m = 0; (m < nrefat); m++)
  83             {
  84                 ai   = index[m];
  85                 mass = top->atoms.atom[ai].m/mtot;
  86                 pbc_dx(pbc, x[ai], xref, dx);
  87                 rvec_add(xref, dx, xtest);
  88                 for (j = 0; (j < DIM); j++)
  89                 {
  90                     if (fabs(xtest[j]-x[ai][j]) > tol)
  91                     {
  92                         /* Here we have used the wrong image for contributing to the COM */
  93                         xref[j] += mass*(xtest[j]-x[ai][j]);
  94                         x[ai][j] = xtest[j];
  95                         bChanged = TRUE;
  96                     }
  97                 }
  98             }
  99             if (bChanged)
 100             {
 101                 printf("COM: %8.3f  %8.3f  %8.3f  iter = %d\n", xref[XX], xref[YY], xref[ZZ], iter);
 102             }
 103             iter++;
 104         }
 105         while (bChanged);
 106     }
 107 }
 108
 109 int gmx_sorient(int argc, char *argv[])
 110 {
 111     t_topology      top;
 112     int             ePBC = -1;
 113     char            title[STRLEN];
 114     t_trxstatus    *status;
 115     int             natoms;
 116     real            t;
 117     rvec           *xtop, *x;
 118     matrix          box;
 119
 120     FILE           *fp;
 121     int             i, j, p, sa0, sa1, sa2, n, ntot, nf, m, *hist1, *hist2, *histn, nbin1, nbin2, nrbin;
 122     real           *histi1, *histi2, invbw, invrbw;
 123     double          sum1, sum2;
 124     int            *isize, nrefgrp, nrefat;
 125     atom_id       **index;
 126     char          **grpname;
 127     real            inp, outp, two_pi, nav, normfac, rmin2, rmax2, rcut, rcut2, r2, r, mass, mtot;
 128     real            c1, c2;
 129     char            str[STRLEN];
 130     gmx_bool        bTPS;
 131     rvec            xref, dx, dxh1, dxh2, outer;
 132     gmx_rmpbc_t     gpbc = NULL;
 133     t_pbc           pbc;
 134     const char     *legr[] = {
 135         "<cos(\\8q\\4\\s1\\N)>",
 136         "<3cos\\S2\\N(\\8q\\4\\s2\\N)-1>"
 137     };
 138     const char     *legc[] = {
 139         "cos(\\8q\\4\\s1\\N)",
 140         "3cos\\S2\\N(\\8q\\4\\s2\\N)-1"
 141     };
 142
 143     const char     *desc[] = {
 144         "[TT]g_sorient[tt] analyzes solvent orientation around solutes.",
 145         "It calculates two angles between the vector from one or more",
 146         "reference positions to the first atom of each solvent molecule:[PAR]",
 147         "[GRK]theta[grk][SUB]1[sub]: the angle with the vector from the first atom of the solvent",
 148         "molecule to the midpoint between atoms 2 and 3.[BR]",
 149         "[GRK]theta[grk][SUB]2[sub]: the angle with the normal of the solvent plane, defined by the",
 150         "same three atoms, or, when the option [TT]-v23[tt] is set, ",
 151         "the angle with the vector between atoms 2 and 3.[PAR]",
 152         "The reference can be a set of atoms or",
 153         "the center of mass of a set of atoms. The group of solvent atoms should",
 154         "consist of 3 atoms per solvent molecule.",
 155         "Only solvent molecules between [TT]-rmin[tt] and [TT]-rmax[tt] are",
 156         "considered for [TT]-o[tt] and [TT]-no[tt] each frame.[PAR]",
 157         "[TT]-o[tt]: distribtion of [MATH][COS][GRK]theta[grk][SUB]1[sub][cos][math] for rmin<=r<=rmax.[PAR]",
 158         "[TT]-no[tt]: distribution of [MATH][COS][GRK]theta[grk][SUB]2[sub][cos][math] for rmin<=r<=rmax.[PAR]",
 159         "[TT]-ro[tt]: [MATH][CHEVRON][COS][GRK]theta[grk][SUB]1[sub][cos][chevron][math] and [MATH][CHEVRON]3[COS]^2[GRK]theta[grk][SUB]2[sub][cos]-1[chevron][math] as a function of the",
 160         "distance.[PAR]",
 161         "[TT]-co[tt]: the sum over all solvent molecules within distance r",
 162         "of [MATH][COS][GRK]theta[grk][SUB]1[sub][cos][math] and [MATH]3[COS]^2([GRK]theta[grk][SUB]2[sub])-1[cos][math] as a function of r.[PAR]",
 163         "[TT]-rc[tt]: the distribution of the solvent molecules as a function of r"
 164     };
 165
 166     output_env_t    oenv;
 167     static gmx_bool bCom = FALSE, bVec23 = FALSE, bPBC = FALSE;
 168     static real     rmin = 0.0, rmax = 0.5, binwidth = 0.02, rbinw = 0.02;
 169     t_pargs         pa[] = {
 170         { "-com",  FALSE, etBOOL,  {&bCom},
 171           "Use the center of mass as the reference postion" },
 172         { "-v23",  FALSE, etBOOL,  {&bVec23},
 173           "Use the vector between atoms 2 and 3" },
 174         { "-rmin",  FALSE, etREAL, {&rmin}, "Minimum distance (nm)" },
 175         { "-rmax",  FALSE, etREAL, {&rmax}, "Maximum distance (nm)" },
 176         { "-cbin",  FALSE, etREAL, {&binwidth}, "Binwidth for the cosine" },
 177         { "-rbin",  FALSE, etREAL, {&rbinw}, "Binwidth for r (nm)" },
 178         { "-pbc",   FALSE, etBOOL, {&bPBC}, "Check PBC for the center of mass calculation. Only necessary when your reference group consists of several molecules." }
 179     };
 180
 181     t_filenm        fnm[] = {
 182         { efTRX, NULL,  NULL,  ffREAD },
 183         { efTPS, NULL,  NULL,  ffREAD },
 184         { efNDX, NULL,  NULL,  ffOPTRD },
 185         { efXVG, NULL,  "sori.xvg",  ffWRITE },
 186         { efXVG, "-no", "snor.xvg",  ffWRITE },
 187         { efXVG, "-ro", "sord.xvg",  ffWRITE },
 188         { efXVG, "-co", "scum.xvg",  ffWRITE },
 189         { efXVG, "-rc", "scount.xvg",  ffWRITE }
 190     };
 191 #define NFILE asize(fnm)
 192
 193     parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_CAN_VIEW | PCA_BE_NICE,
 194                       NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, NULL, &oenv);
 195
 196     two_pi = 2/M_PI;
 197
 198     bTPS = (opt2bSet("-s", NFILE, fnm) || !opt2bSet("-n", NFILE, fnm) || bCom);
 199     if (bTPS)
 200     {
 201         read_tps_conf(ftp2fn(efTPS, NFILE, fnm), title, &top, &ePBC, &xtop, NULL, box,
 202                       bCom);
 203     }
 204
 205     /* get index groups */
 206     printf("Select a group of reference particles and a solvent group:\n");
 207     snew(grpname, 2);
 208     snew(index, 2);
 209     snew(isize, 2);
 210     if (bTPS)
 211     {
 212         get_index(&top.atoms, ftp2fn_null(efNDX, NFILE, fnm), 2, isize, index, grpname);
 213     }
 214     else
 215     {
 216         get_index(NULL, ftp2fn(efNDX, NFILE, fnm), 2, isize, index, grpname);
 217     }
 218
 219     if (bCom)
 220     {
 221         nrefgrp = 1;
 222         nrefat  = isize[0];
 223     }
 224     else
 225     {
 226         nrefgrp = isize[0];
 227         nrefat  = 1;
 228     }
 229
 230     if (isize[1] % 3)
 231     {
 232         gmx_fatal(FARGS, "The number of solvent atoms (%d) is not a multiple of 3",
 233                   isize[1]);
 234     }
 235
 236     /* initialize reading trajectory:                         */
 237     natoms = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);
 238
 239     rmin2 = sqr(rmin);
 240     rmax2 = sqr(rmax);
 241     rcut  = 0.99*sqrt(max_cutoff2(guess_ePBC(box), box));
 242     if (rcut == 0)
 243     {
 244         rcut = 10*rmax;
 245     }
 246     rcut2 = sqr(rcut);
 247
 248     invbw = 1/binwidth;
 249     nbin1 = 1+(int)(2*invbw + 0.5);
 250     nbin2 = 1+(int)(invbw + 0.5);
 251
 252     invrbw = 1/rbinw;
 253
 254     snew(hist1, nbin1);
 255     snew(hist2, nbin2);
 256     nrbin = 1+(int)(rcut/rbinw);
 257     if (nrbin == 0)
 258     {
 259         nrbin = 1;
 260     }
 261     snew(histi1, nrbin);
 262     snew(histi2, nrbin);
 263     snew(histn, nrbin);
 264
 265     ntot = 0;
 266     nf   = 0;
 267     sum1 = 0;
 268     sum2 = 0;
 269
 270     if (bTPS)
 271     {
 272         /* make molecules whole again */
 273         gpbc = gmx_rmpbc_init(&top.idef, ePBC, natoms, box);
 274     }
 275     /* start analysis of trajectory */
 276     do
 277     {
 278         if (bTPS)
 279         {
 280             /* make molecules whole again */
 281             gmx_rmpbc(gpbc, natoms, box, x);
 282         }
 283
 284         set_pbc(&pbc, ePBC, box);
 285         n    = 0;
 286         inp  = 0;
 287         outp = 0;
 288         for (p = 0; (p < nrefgrp); p++)
 289         {
 290             if (bCom)
 291             {
 292                 calc_com_pbc(nrefat, &top, x, &pbc, index[0], xref, bPBC, box);
 293             }
 294             else
 295             {
 296                 copy_rvec(x[index[0][p]], xref);
 297             }
 298
 299             for (m = 0; m < isize[1]; m += 3)
 300             {
 301                 sa0 = index[1][m];
 302                 sa1 = index[1][m+1];
 303                 sa2 = index[1][m+2];
 304                 range_check(sa0, 0, natoms);
 305                 range_check(sa1, 0, natoms);
 306                 range_check(sa2, 0, natoms);
 307                 pbc_dx(&pbc, x[sa0], xref, dx);
 308                 r2  = norm2(dx);
 309                 if (r2 < rcut2)
 310                 {
 311                     r = sqrt(r2);
 312                     if (!bVec23)
 313                     {
 314                         /* Determine the normal to the plain */
 315                         rvec_sub(x[sa1], x[sa0], dxh1);
 316                         rvec_sub(x[sa2], x[sa0], dxh2);
 317                         rvec_inc(dxh1, dxh2);
 318                         svmul(1/r, dx, dx);
 319                         unitv(dxh1, dxh1);
 320                         inp = iprod(dx, dxh1);
 321                         cprod(dxh1, dxh2, outer);
 322                         unitv(outer, outer);
 323                         outp = iprod(dx, outer);
 324                     }
 325                     else
 326                     {
 327                         /* Use the vector between the 2nd and 3rd atom */
 328                         rvec_sub(x[sa2], x[sa1], dxh2);
 329                         unitv(dxh2, dxh2);
 330                         outp = iprod(dx, dxh2)/r;
 331                     }
 332                     {
 333                         int ii = (int)(invrbw*r);
 334                         range_check(ii, 0, nrbin);
 335                         histi1[ii] += inp;
 336                         histi2[ii] += 3*sqr(outp) - 1;
 337                         histn[ii]++;
 338                     }
 339                     if ((r2 >= rmin2) && (r2 < rmax2))
 340                     {
 341                         int ii1 = (int)(invbw*(inp + 1));
 342                         int ii2 = (int)(invbw*fabs(outp));
 343
 344                         range_check(ii1, 0, nbin1);
 345                         range_check(ii2, 0, nbin2);
 346                         hist1[ii1]++;
 347                         hist2[ii2]++;
 348                         sum1 += inp;
 349                         sum2 += outp;
 350                         n++;
 351                     }
 352                 }
 353             }
 354         }
 355         ntot += n;
 356         nf++;
 357
 358     }
 359     while (read_next_x(oenv, status, &t, natoms, x, box));
 360
 361     /* clean up */
 362     sfree(x);
 363     close_trj(status);
 364     gmx_rmpbc_done(gpbc);
 365
 366     /* Add the bin for the exact maximum to the previous bin */
 367     hist1[nbin1-1] += hist1[nbin1];
 368     hist2[nbin2-1] += hist2[nbin2];
 369
 370     nav     = (real)ntot/(nrefgrp*nf);
 371     normfac = invbw/ntot;
 372
 373     fprintf(stderr,  "Average nr of molecules between %g and %g nm: %.1f\n",
 374             rmin, rmax, nav);
 375     if (ntot > 0)
 376     {
 377         sum1 /= ntot;
 378         sum2 /= ntot;
 379         fprintf(stderr, "Average cos(theta1)     between %g and %g nm: %6.3f\n",
 380                 rmin, rmax, sum1);
 381         fprintf(stderr, "Average 3cos2(theta2)-1 between %g and %g nm: %6.3f\n",
 382                 rmin, rmax, sum2);
 383     }
 384
 385     sprintf(str, "Solvent orientation between %g and %g nm", rmin, rmax);
 386     fp = xvgropen(opt2fn("-o", NFILE, fnm), str, "cos(\\8q\\4\\s1\\N)", "", oenv);
 387     if (output_env_get_print_xvgr_codes(oenv))
 388     {
 389         fprintf(fp, "@ subtitle \"average shell size %.1f molecules\"\n", nav);
 390     }
 391     for (i = 0; i < nbin1; i++)
 392     {
 393         fprintf(fp, "%g %g\n", (i+0.5)*binwidth-1, 2*normfac*hist1[i]);
 394     }
 395     ffclose(fp);
 396
 397     sprintf(str, "Solvent normal orientation between %g and %g nm", rmin, rmax);
 398     fp = xvgropen(opt2fn("-no", NFILE, fnm), str, "cos(\\8q\\4\\s2\\N)", "", oenv);
 399     if (output_env_get_print_xvgr_codes(oenv))
 400     {
 401         fprintf(fp, "@ subtitle \"average shell size %.1f molecules\"\n", nav);
 402     }
 403     for (i = 0; i < nbin2; i++)
 404     {
 405         fprintf(fp, "%g %g\n", (i+0.5)*binwidth, normfac*hist2[i]);
 406     }
 407     ffclose(fp);
 408
 409
 410     sprintf(str, "Solvent orientation");
 411     fp = xvgropen(opt2fn("-ro", NFILE, fnm), str, "r (nm)", "", oenv);
 412     if (output_env_get_print_xvgr_codes(oenv))
 413     {
 414         fprintf(fp, "@ subtitle \"as a function of distance\"\n");
 415     }
 416     xvgr_legend(fp, 2, legr, oenv);
 417     for (i = 0; i < nrbin; i++)
 418     {
 419         fprintf(fp, "%g %g %g\n", (i+0.5)*rbinw,
 420                 histn[i] ? histi1[i]/histn[i] : 0,
 421                 histn[i] ? histi2[i]/histn[i] : 0);
 422     }
 423     ffclose(fp);
 424
 425     sprintf(str, "Cumulative solvent orientation");
 426     fp = xvgropen(opt2fn("-co", NFILE, fnm), str, "r (nm)", "", oenv);
 427     if (output_env_get_print_xvgr_codes(oenv))
 428     {
 429         fprintf(fp, "@ subtitle \"as a function of distance\"\n");
 430     }
 431     xvgr_legend(fp, 2, legc, oenv);
 432     normfac = 1.0/(nrefgrp*nf);
 433     c1      = 0;
 434     c2      = 0;
 435     fprintf(fp, "%g %g %g\n", 0.0, c1, c2);
 436     for (i = 0; i < nrbin; i++)
 437     {
 438         c1 += histi1[i]*normfac;
 439         c2 += histi2[i]*normfac;
 440         fprintf(fp, "%g %g %g\n", (i+1)*rbinw, c1, c2);
 441     }
 442     ffclose(fp);
 443
 444     sprintf(str, "Solvent distribution");
 445     fp = xvgropen(opt2fn("-rc", NFILE, fnm), str, "r (nm)", "molecules/nm", oenv);
 446     if (output_env_get_print_xvgr_codes(oenv))
 447     {
 448         fprintf(fp, "@ subtitle \"as a function of distance\"\n");
 449     }
 450     normfac = 1.0/(rbinw*nf);
 451     for (i = 0; i < nrbin; i++)
 452     {
 453         fprintf(fp, "%g %g\n", (i+0.5)*rbinw, histn[i]*normfac);
 454     }
 455     ffclose(fp);
 456
 457     do_view(oenv, opt2fn("-o", NFILE, fnm), NULL);
 458     do_view(oenv, opt2fn("-no", NFILE, fnm), NULL);
 459     do_view(oenv, opt2fn("-ro", NFILE, fnm), "-nxy");
 460     do_view(oenv, opt2fn("-co", NFILE, fnm), "-nxy");
 461
 462     thanx(stderr);
 463
 464     return 0;
 465 }