src/gromacs/gmxana/gmx_saltbr.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  * Copyright (c) 2013,2014, by the GROMACS development team, led by
   7  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   8  * and including many others, as listed in the AUTHORS file in the
   9  * top-level source directory and at http://www.gromacs.org.
  10  *
  11  * GROMACS is free software; you can redistribute it and/or
  12  * modify it under the terms of the GNU Lesser General Public License
  13  * as published by the Free Software Foundation; either version 2.1
  14  * of the License, or (at your option) any later version.
  15  *
  16  * GROMACS is distributed in the hope that it will be useful,
  17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19  * Lesser General Public License for more details.
  20  *
  21  * You should have received a copy of the GNU Lesser General Public
  22  * License along with GROMACS; if not, see
  23  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  24  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  25  *
  26  * If you want to redistribute modifications to GROMACS, please
  27  * consider that scientific software is very special. Version
  28  * control is crucial - bugs must be traceable. We will be happy to
  29  * consider code for inclusion in the official distribution, but
  30  * derived work must not be called official GROMACS. Details are found
  31  * in the README & COPYING files - if they are missing, get the
  32  * official version at http://www.gromacs.org.
  33  *
  34  * To help us fund GROMACS development, we humbly ask that you cite
  35  * the research papers on the package. Check out http://www.gromacs.org.
  36  */
  37 #include "gmxpre.h"
  38
  39 #include <math.h>
  40 #include <string.h>
  41
  42 #include "gromacs/commandline/pargs.h"
  43 #include "gromacs/fileio/filenm.h"
  44 #include "gromacs/fileio/trxio.h"
  45 #include "gromacs/fileio/xvgr.h"
  46 #include "gromacs/gmxana/gmx_ana.h"
  47 #include "gromacs/legacyheaders/macros.h"
  48 #include "gromacs/legacyheaders/typedefs.h"
  49 #include "gromacs/math/vec.h"
  50 #include "gromacs/pbcutil/pbc.h"
  51 #include "gromacs/utility/fatalerror.h"
  52 #include "gromacs/utility/futil.h"
  53 #include "gromacs/utility/smalloc.h"
  54
  55 typedef struct {
  56     char *label;
  57     int   cg;
  58     real  q;
  59 } t_charge;
  60
  61 static t_charge *mk_charge(t_atoms *atoms, t_block *cgs, int *nncg)
  62 {
  63     t_charge *cg = NULL;
  64     char      buf[32];
  65     int       i, j, ncg, resnr, anr;
  66     real      qq;
  67
  68     /* Find the charged groups */
  69     ncg = 0;
  70     for (i = 0; (i < cgs->nr); i++)
  71     {
  72         qq = 0.0;
  73         for (j = cgs->index[i]; (j < cgs->index[i+1]); j++)
  74         {
  75             qq += atoms->atom[j].q;
  76         }
  77         if (fabs(qq) > 1.0e-5)
  78         {
  79             srenew(cg, ncg+1);
  80             cg[ncg].q  = qq;
  81             cg[ncg].cg = i;
  82             anr        = cgs->index[i];
  83             resnr      = atoms->atom[anr].resind;
  84             sprintf(buf, "%s%d-%d",
  85                     *(atoms->resinfo[resnr].name),
  86                     atoms->resinfo[resnr].nr,
  87                     anr+1);
  88             cg[ncg].label = gmx_strdup(buf);
  89             ncg++;
  90         }
  91     }
  92     *nncg = ncg;
  93
  94     for (i = 0; (i < ncg); i++)
  95     {
  96         printf("CG: %10s Q: %6g  Atoms:",
  97                cg[i].label, cg[i].q);
  98         for (j = cgs->index[cg[i].cg]; (j < cgs->index[cg[i].cg+1]); j++)
  99         {
 100             printf(" %4d", j);
 101         }
 102         printf("\n");
 103     }
 104
 105     return cg;
 106 }
 107
 108 static real calc_dist(t_pbc *pbc, rvec x[], t_block *cgs, int icg, int jcg)
 109 {
 110     int  i, j;
 111     rvec dx;
 112     real d2, mindist2 = 1000;
 113
 114     for (i = cgs->index[icg]; (i < cgs->index[icg+1]); i++)
 115     {
 116         for (j = cgs->index[jcg]; (j < cgs->index[jcg+1]); j++)
 117         {
 118             pbc_dx(pbc, x[i], x[j], dx);
 119             d2 = norm2(dx);
 120             if (d2 < mindist2)
 121             {
 122                 mindist2 = d2;
 123             }
 124         }
 125     }
 126     return sqrt(mindist2);
 127 }
 128
 129 int gmx_saltbr(int argc, char *argv[])
 130 {
 131     const char     *desc[] = {
 132         "[THISMODULE] plots the distance between all combination of charged groups",
 133         "as a function of time. The groups are combined in different ways.",
 134         "A minimum distance can be given (i.e. a cut-off), such that groups",
 135         "that are never closer than that distance will not be plotted.[PAR]",
 136         "Output will be in a number of fixed filenames, [TT]min-min.xvg[tt], [TT]plus-min.xvg[tt]",
 137         "and [TT]plus-plus.xvg[tt], or files for every individual ion pair if the [TT]-sep[tt]",
 138         "option is selected. In this case, files are named as [TT]sb-(Resname)(Resnr)-(Atomnr)[tt].",
 139         "There may be [BB]many[bb] such files."
 140     };
 141     static gmx_bool bSep     = FALSE;
 142     static real     truncate = 1000.0;
 143     t_pargs         pa[]     = {
 144         { "-t",   FALSE, etREAL, {&truncate},
 145           "Groups that are never closer than this distance are not plotted" },
 146         { "-sep", FALSE, etBOOL, {&bSep},
 147           "Use separate files for each interaction (may be MANY)" }
 148     };
 149     t_filenm        fnm[] = {
 150         { efTRX, "-f",  NULL, ffREAD },
 151         { efTPR, NULL,  NULL, ffREAD },
 152     };
 153 #define NFILE asize(fnm)
 154
 155     FILE              *out[3], *fp;
 156     static const char *title[3] = {
 157         "Distance between positively charged groups",
 158         "Distance between negatively charged groups",
 159         "Distance between oppositely charged groups"
 160     };
 161     static const char *fn[3] = {
 162         "plus-plus.xvg",
 163         "min-min.xvg",
 164         "plus-min.xvg"
 165     };
 166     int                nset[3] = {0, 0, 0};
 167
 168     t_topology        *top;
 169     int                ePBC;
 170     char              *buf;
 171     t_trxstatus       *status;
 172     int                i, j, k, m, nnn, teller, ncg, n1, n2, n3, natoms;
 173     real               t, *time, qi, qj;
 174     t_charge          *cg;
 175     real            ***cgdist;
 176     int              **nWithin;
 177
 178     double             t0, dt;
 179     char               label[234];
 180     t_pbc              pbc;
 181     rvec              *x;
 182     matrix             box;
 183     output_env_t       oenv;
 184
 185     if (!parse_common_args(&argc, argv, PCA_CAN_TIME,
 186                            NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, NULL, &oenv))
 187     {
 188         return 0;
 189     }
 190
 191     top = read_top(ftp2fn(efTPR, NFILE, fnm), &ePBC);
 192     cg  = mk_charge(&top->atoms, &(top->cgs), &ncg);
 193     snew(cgdist, ncg);
 194     snew(nWithin, ncg);
 195     for (i = 0; (i < ncg); i++)
 196     {
 197         snew(cgdist[i], ncg);
 198         snew(nWithin[i], ncg);
 199     }
 200
 201     natoms = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);
 202
 203     teller = 0;
 204     time   = NULL;
 205     do
 206     {
 207         srenew(time, teller+1);
 208         time[teller] = t;
 209
 210         set_pbc(&pbc, ePBC, box);
 211
 212         for (i = 0; (i < ncg); i++)
 213         {
 214             for (j = i+1; (j < ncg); j++)
 215             {
 216                 srenew(cgdist[i][j], teller+1);
 217                 cgdist[i][j][teller] =
 218                     calc_dist(&pbc, x, &(top->cgs), cg[i].cg, cg[j].cg);
 219                 if (cgdist[i][j][teller] < truncate)
 220                 {
 221                     nWithin[i][j] = 1;
 222                 }
 223             }
 224         }
 225
 226         teller++;
 227     }
 228     while (read_next_x(oenv, status, &t, x, box));
 229     fprintf(stderr, "\n");
 230     close_trj(status);
 231
 232     if (bSep)
 233     {
 234         snew(buf, 256);
 235         for (i = 0; (i < ncg); i++)
 236         {
 237             for (j = i+1; (j < ncg); j++)
 238             {
 239                 if (nWithin[i][j])
 240                 {
 241                     sprintf(buf, "sb-%s:%s.xvg", cg[i].label, cg[j].label);
 242                     fp = xvgropen(buf, buf, "Time (ps)", "Distance (nm)", oenv);
 243                     for (k = 0; (k < teller); k++)
 244                     {
 245                         fprintf(fp, "%10g  %10g\n", time[k], cgdist[i][j][k]);
 246                     }
 247                     gmx_ffclose(fp);
 248                 }
 249             }
 250         }
 251         sfree(buf);
 252     }
 253     else
 254     {
 255
 256         for (m = 0; (m < 3); m++)
 257         {
 258             out[m] = xvgropen(fn[m], title[m], "Time (ps)", "Distance (nm)", oenv);
 259         }
 260
 261         snew(buf, 256);
 262         for (i = 0; (i < ncg); i++)
 263         {
 264             qi = cg[i].q;
 265             for (j = i+1; (j < ncg); j++)
 266             {
 267                 qj = cg[j].q;
 268                 if (nWithin[i][j])
 269                 {
 270                     sprintf(buf, "%s:%s", cg[i].label, cg[j].label);
 271                     if (qi*qj < 0)
 272                     {
 273                         nnn = 2;
 274                     }
 275                     else if (qi+qj > 0)
 276                     {
 277                         nnn = 0;
 278                     }
 279                     else
 280                     {
 281                         nnn = 1;
 282                     }
 283
 284                     if (nset[nnn] == 0)
 285                     {
 286                         xvgr_legend(out[nnn], 1, (const char**)&buf, oenv);
 287                     }
 288                     else
 289                     {
 290                         if (output_env_get_xvg_format(oenv) == exvgXMGR)
 291                         {
 292                             fprintf(out[nnn], "@ legend string %d \"%s\"\n", nset[nnn], buf);
 293                         }
 294                         else if (output_env_get_xvg_format(oenv) == exvgXMGRACE)
 295                         {
 296                             fprintf(out[nnn], "@ s%d legend \"%s\"\n", nset[nnn], buf);
 297                         }
 298                     }
 299                     nset[nnn]++;
 300                     nWithin[i][j] = nnn+1;
 301                 }
 302             }
 303         }
 304         for (k = 0; (k < teller); k++)
 305         {
 306             for (m = 0; (m < 3); m++)
 307             {
 308                 fprintf(out[m], "%10g", time[k]);
 309             }
 310
 311             for (i = 0; (i < ncg); i++)
 312             {
 313                 for (j = i+1; (j < ncg); j++)
 314                 {
 315                     nnn = nWithin[i][j];
 316                     if (nnn > 0)
 317                     {
 318                         fprintf(out[nnn-1], "  %10g", cgdist[i][j][k]);
 319                     }
 320                 }
 321             }
 322             for (m = 0; (m < 3); m++)
 323             {
 324                 fprintf(out[m], "\n");
 325             }
 326         }
 327         for (m = 0; (m < 3); m++)
 328         {
 329             gmx_ffclose(out[m]);
 330             if (nset[m] == 0)
 331             {
 332                 remove(fn[m]);
 333             }
 334         }
 335     }
 336
 337     return 0;
 338 }