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 #ifdef HAVE_CONFIG_H
  38 #include <config.h>
  39 #endif
  40 #include <math.h>
  41 #include <string.h>
  42
  43 #include "macros.h"
  44 #include "vec.h"
  45 #include "sysstuff.h"
  46 #include "typedefs.h"
  47 #include "gromacs/fileio/filenm.h"
  48 #include "gromacs/fileio/trxio.h"
  49 #include "gromacs/commandline/pargs.h"
  50 #include "gromacs/fileio/futil.h"
  51 #include "gmx_fatal.h"
  52 #include "gromacs/utility/smalloc.h"
  53 #include "pbc.h"
  54 #include "xvgr.h"
  55 #include "gmx_ana.h"
  56
  57 #include "gromacs/legacyheaders/gmx_fatal.h"
  58
  59 typedef struct {
  60     char *label;
  61     int   cg;
  62     real  q;
  63 } t_charge;
  64
  65 static t_charge *mk_charge(t_atoms *atoms, t_block *cgs, int *nncg)
  66 {
  67     t_charge *cg = NULL;
  68     char      buf[32];
  69     int       i, j, ncg, resnr, anr;
  70     real      qq;
  71
  72     /* Find the charged groups */
  73     ncg = 0;
  74     for (i = 0; (i < cgs->nr); i++)
  75     {
  76         qq = 0.0;
  77         for (j = cgs->index[i]; (j < cgs->index[i+1]); j++)
  78         {
  79             qq += atoms->atom[j].q;
  80         }
  81         if (fabs(qq) > 1.0e-5)
  82         {
  83             srenew(cg, ncg+1);
  84             cg[ncg].q  = qq;
  85             cg[ncg].cg = i;
  86             anr        = cgs->index[i];
  87             resnr      = atoms->atom[anr].resind;
  88             sprintf(buf, "%s%d-%d",
  89                     *(atoms->resinfo[resnr].name),
  90                     atoms->resinfo[resnr].nr,
  91                     anr+1);
  92             cg[ncg].label = strdup(buf);
  93             ncg++;
  94         }
  95     }
  96     *nncg = ncg;
  97
  98     for (i = 0; (i < ncg); i++)
  99     {
 100         printf("CG: %10s Q: %6g  Atoms:",
 101                cg[i].label, cg[i].q);
 102         for (j = cgs->index[cg[i].cg]; (j < cgs->index[cg[i].cg+1]); j++)
 103         {
 104             printf(" %4d", j);
 105         }
 106         printf("\n");
 107     }
 108
 109     return cg;
 110 }
 111
 112 static real calc_dist(t_pbc *pbc, rvec x[], t_block *cgs, int icg, int jcg)
 113 {
 114     int  i, j;
 115     rvec dx;
 116     real d2, mindist2 = 1000;
 117
 118     for (i = cgs->index[icg]; (i < cgs->index[icg+1]); i++)
 119     {
 120         for (j = cgs->index[jcg]; (j < cgs->index[jcg+1]); j++)
 121         {
 122             pbc_dx(pbc, x[i], x[j], dx);
 123             d2 = norm2(dx);
 124             if (d2 < mindist2)
 125             {
 126                 mindist2 = d2;
 127             }
 128         }
 129     }
 130     return sqrt(mindist2);
 131 }
 132
 133 int gmx_saltbr(int argc, char *argv[])
 134 {
 135     const char     *desc[] = {
 136         "[THISMODULE] plots the distance between all combination of charged groups",
 137         "as a function of time. The groups are combined in different ways.",
 138         "A minimum distance can be given (i.e. a cut-off), such that groups",
 139         "that are never closer than that distance will not be plotted.[PAR]",
 140         "Output will be in a number of fixed filenames, [TT]min-min.xvg[tt], [TT]plus-min.xvg[tt]",
 141         "and [TT]plus-plus.xvg[tt], or files for every individual ion pair if the [TT]-sep[tt]",
 142         "option is selected. In this case, files are named as [TT]sb-(Resname)(Resnr)-(Atomnr)[tt].",
 143         "There may be [BB]many[bb] such files."
 144     };
 145     static gmx_bool bSep     = FALSE;
 146     static real     truncate = 1000.0;
 147     t_pargs         pa[]     = {
 148         { "-t",   FALSE, etREAL, {&truncate},
 149           "Groups that are never closer than this distance are not plotted" },
 150         { "-sep", FALSE, etBOOL, {&bSep},
 151           "Use separate files for each interaction (may be MANY)" }
 152     };
 153     t_filenm        fnm[] = {
 154         { efTRX, "-f",  NULL, ffREAD },
 155         { efTPX, NULL,  NULL, ffREAD },
 156     };
 157 #define NFILE asize(fnm)
 158
 159     FILE              *out[3], *fp;
 160     static const char *title[3] = {
 161         "Distance between positively charged groups",
 162         "Distance between negatively charged groups",
 163         "Distance between oppositely charged groups"
 164     };
 165     static const char *fn[3] = {
 166         "plus-plus.xvg",
 167         "min-min.xvg",
 168         "plus-min.xvg"
 169     };
 170     int                nset[3] = {0, 0, 0};
 171
 172     t_topology        *top;
 173     int                ePBC;
 174     char              *buf;
 175     t_trxstatus       *status;
 176     int                i, j, k, m, nnn, teller, ncg, n1, n2, n3, natoms;
 177     real               t, *time, qi, qj;
 178     t_charge          *cg;
 179     real            ***cgdist;
 180     int              **nWithin;
 181
 182     double             t0, dt;
 183     char               label[234];
 184     t_pbc              pbc;
 185     rvec              *x;
 186     matrix             box;
 187     output_env_t       oenv;
 188
 189     if (!parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_BE_NICE,
 190                            NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, NULL, &oenv))
 191     {
 192         return 0;
 193     }
 194
 195     top = read_top(ftp2fn(efTPX, NFILE, fnm), &ePBC);
 196     cg  = mk_charge(&top->atoms, &(top->cgs), &ncg);
 197     snew(cgdist, ncg);
 198     snew(nWithin, ncg);
 199     for (i = 0; (i < ncg); i++)
 200     {
 201         snew(cgdist[i], ncg);
 202         snew(nWithin[i], ncg);
 203     }
 204
 205     natoms = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);
 206
 207     teller = 0;
 208     time   = NULL;
 209     do
 210     {
 211         srenew(time, teller+1);
 212         time[teller] = t;
 213
 214         set_pbc(&pbc, ePBC, box);
 215
 216         for (i = 0; (i < ncg); i++)
 217         {
 218             for (j = i+1; (j < ncg); j++)
 219             {
 220                 srenew(cgdist[i][j], teller+1);
 221                 cgdist[i][j][teller] =
 222                     calc_dist(&pbc, x, &(top->cgs), cg[i].cg, cg[j].cg);
 223                 if (cgdist[i][j][teller] < truncate)
 224                 {
 225                     nWithin[i][j] = 1;
 226                 }
 227             }
 228         }
 229
 230         teller++;
 231     }
 232     while (read_next_x(oenv, status, &t, x, box));
 233     fprintf(stderr, "\n");
 234     close_trj(status);
 235
 236     if (bSep)
 237     {
 238         snew(buf, 256);
 239         for (i = 0; (i < ncg); i++)
 240         {
 241             for (j = i+1; (j < ncg); j++)
 242             {
 243                 if (nWithin[i][j])
 244                 {
 245                     sprintf(buf, "sb-%s:%s.xvg", cg[i].label, cg[j].label);
 246                     fp = xvgropen(buf, buf, "Time (ps)", "Distance (nm)", oenv);
 247                     for (k = 0; (k < teller); k++)
 248                     {
 249                         fprintf(fp, "%10g  %10g\n", time[k], cgdist[i][j][k]);
 250                     }
 251                     gmx_ffclose(fp);
 252                 }
 253             }
 254         }
 255         sfree(buf);
 256     }
 257     else
 258     {
 259
 260         for (m = 0; (m < 3); m++)
 261         {
 262             out[m] = xvgropen(fn[m], title[m], "Time (ps)", "Distance (nm)", oenv);
 263         }
 264
 265         snew(buf, 256);
 266         for (i = 0; (i < ncg); i++)
 267         {
 268             qi = cg[i].q;
 269             for (j = i+1; (j < ncg); j++)
 270             {
 271                 qj = cg[j].q;
 272                 if (nWithin[i][j])
 273                 {
 274                     sprintf(buf, "%s:%s", cg[i].label, cg[j].label);
 275                     if (qi*qj < 0)
 276                     {
 277                         nnn = 2;
 278                     }
 279                     else if (qi+qj > 0)
 280                     {
 281                         nnn = 0;
 282                     }
 283                     else
 284                     {
 285                         nnn = 1;
 286                     }
 287
 288                     if (nset[nnn] == 0)
 289                     {
 290                         xvgr_legend(out[nnn], 1, (const char**)&buf, oenv);
 291                     }
 292                     else
 293                     {
 294                         if (output_env_get_xvg_format(oenv) == exvgXMGR)
 295                         {
 296                             fprintf(out[nnn], "@ legend string %d \"%s\"\n", nset[nnn], buf);
 297                         }
 298                         else if (output_env_get_xvg_format(oenv) == exvgXMGRACE)
 299                         {
 300                             fprintf(out[nnn], "@ s%d legend \"%s\"\n", nset[nnn], buf);
 301                         }
 302                     }
 303                     nset[nnn]++;
 304                     nWithin[i][j] = nnn+1;
 305                 }
 306             }
 307         }
 308         for (k = 0; (k < teller); k++)
 309         {
 310             for (m = 0; (m < 3); m++)
 311             {
 312                 fprintf(out[m], "%10g", time[k]);
 313             }
 314
 315             for (i = 0; (i < ncg); i++)
 316             {
 317                 for (j = i+1; (j < ncg); j++)
 318                 {
 319                     nnn = nWithin[i][j];
 320                     if (nnn > 0)
 321                     {
 322                         fprintf(out[nnn-1], "  %10g", cgdist[i][j][k]);
 323                     }
 324                 }
 325             }
 326             for (m = 0; (m < 3); m++)
 327             {
 328                 fprintf(out[m], "\n");
 329             }
 330         }
 331         for (m = 0; (m < 3); m++)
 332         {
 333             gmx_ffclose(out[m]);
 334             if (nset[m] == 0)
 335             {
 336                 remove(fn[m]);
 337             }
 338         }
 339     }
 340
 341     return 0;
 342 }