src/gromacs/gmxana/gmx_genion.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
  41 #include <ctype.h>
  42 #include <string.h>
  43
  44 #include "gromacs/utility/cstringutil.h"
  45 #include "gromacs/utility/smalloc.h"
  46 #include "sysstuff.h"
  47 #include "gromacs/fileio/confio.h"
  48 #include "gromacs/commandline/pargs.h"
  49 #include "pbc.h"
  50 #include "force.h"
  51 #include "gmx_fatal.h"
  52 #include "gromacs/fileio/futil.h"
  53 #include "gromacs/math/utilities.h"
  54 #include "macros.h"
  55 #include "vec.h"
  56 #include "gromacs/fileio/tpxio.h"
  57 #include "mdrun.h"
  58 #include "main.h"
  59 #include "gromacs/random/random.h"
  60 #include "index.h"
  61 #include "mtop_util.h"
  62 #include "gmx_ana.h"
  63
  64 static void insert_ion(int nsa, int *nwater,
  65                        gmx_bool bSet[], int repl[], atom_id index[],
  66                        rvec x[], t_pbc *pbc,
  67                        int sign, int q, const char *ionname,
  68                        t_atoms *atoms,
  69                        real rmin, gmx_rng_t rng)
  70 {
  71     int             i, ei, nw;
  72     real            rmin2;
  73     rvec            dx;
  74     gmx_int64_t     maxrand;
  75
  76     ei       = -1;
  77     nw       = *nwater;
  78     maxrand  = nw;
  79     maxrand *= 1000;
  80
  81     do
  82     {
  83         ei = nw*gmx_rng_uniform_real(rng);
  84         maxrand--;
  85     }
  86     while (bSet[ei] && (maxrand > 0));
  87     if (bSet[ei])
  88     {
  89         gmx_fatal(FARGS, "No more replaceable solvent!");
  90     }
  91
  92     fprintf(stderr, "Replacing solvent molecule %d (atom %d) with %s\n",
  93             ei, index[nsa*ei], ionname);
  94
  95     /* Replace solvent molecule charges with ion charge */
  96     bSet[ei] = TRUE;
  97     repl[ei] = sign;
  98
  99     atoms->atom[index[nsa*ei]].q = q;
 100     for (i = 1; i < nsa; i++)
 101     {
 102         atoms->atom[index[nsa*ei+i]].q = 0;
 103     }
 104
 105     /* Mark all solvent molecules within rmin as unavailable for substitution */
 106     if (rmin > 0)
 107     {
 108         rmin2 = rmin*rmin;
 109         for (i = 0; (i < nw); i++)
 110         {
 111             if (!bSet[i])
 112             {
 113                 pbc_dx(pbc, x[index[nsa*ei]], x[index[nsa*i]], dx);
 114                 if (iprod(dx, dx) < rmin2)
 115                 {
 116                     bSet[i] = TRUE;
 117                 }
 118             }
 119         }
 120     }
 121 }
 122
 123
 124 static char *aname(const char *mname)
 125 {
 126     char *str;
 127     int   i;
 128
 129     str = strdup(mname);
 130     i   = strlen(str)-1;
 131     while (i > 1 && (isdigit(str[i]) || (str[i] == '+') || (str[i] == '-')))
 132     {
 133         str[i] = '\0';
 134         i--;
 135     }
 136
 137     return str;
 138 }
 139
 140 void sort_ions(int nsa, int nw, int repl[], atom_id index[],
 141                t_atoms *atoms, rvec x[],
 142                const char *p_name, const char *n_name)
 143 {
 144     int    i, j, k, r, np, nn, starta, startr, npi, nni;
 145     rvec  *xt;
 146     char **pptr = NULL, **nptr = NULL, **paptr = NULL, **naptr = NULL;
 147
 148     snew(xt, atoms->nr);
 149
 150     /* Put all the solvent in front and count the added ions */
 151     np = 0;
 152     nn = 0;
 153     j  = index[0];
 154     for (i = 0; i < nw; i++)
 155     {
 156         r = repl[i];
 157         if (r == 0)
 158         {
 159             for (k = 0; k < nsa; k++)
 160             {
 161                 copy_rvec(x[index[nsa*i+k]], xt[j++]);
 162             }
 163         }
 164         else if (r > 0)
 165         {
 166             np++;
 167         }
 168         else if (r < 0)
 169         {
 170             nn++;
 171         }
 172     }
 173
 174     if (np+nn > 0)
 175     {
 176         /* Put the positive and negative ions at the end */
 177         starta = index[nsa*(nw - np - nn)];
 178         startr = atoms->atom[starta].resind;
 179
 180         if (np)
 181         {
 182             snew(pptr, 1);
 183             pptr[0] = strdup(p_name);
 184             snew(paptr, 1);
 185             paptr[0] = aname(p_name);
 186         }
 187         if (nn)
 188         {
 189             snew(nptr, 1);
 190             nptr[0] = strdup(n_name);
 191             snew(naptr, 1);
 192             naptr[0] = aname(n_name);
 193         }
 194         npi = 0;
 195         nni = 0;
 196         for (i = 0; i < nw; i++)
 197         {
 198             r = repl[i];
 199             if (r > 0)
 200             {
 201                 j = starta+npi;
 202                 k = startr+npi;
 203                 copy_rvec(x[index[nsa*i]], xt[j]);
 204                 atoms->atomname[j]     = paptr;
 205                 atoms->atom[j].resind  = k;
 206                 atoms->resinfo[k].name = pptr;
 207                 npi++;
 208             }
 209             else if (r < 0)
 210             {
 211                 j = starta+np+nni;
 212                 k = startr+np+nni;
 213                 copy_rvec(x[index[nsa*i]], xt[j]);
 214                 atoms->atomname[j]     = naptr;
 215                 atoms->atom[j].resind  = k;
 216                 atoms->resinfo[k].name = nptr;
 217                 nni++;
 218             }
 219         }
 220         for (i = index[nsa*nw-1]+1; i < atoms->nr; i++)
 221         {
 222             j                  = i-(nsa-1)*(np+nn);
 223             atoms->atomname[j] = atoms->atomname[i];
 224             atoms->atom[j]     = atoms->atom[i];
 225             copy_rvec(x[i], xt[j]);
 226         }
 227         atoms->nr -= (nsa-1)*(np+nn);
 228
 229         /* Copy the new positions back */
 230         for (i = index[0]; i < atoms->nr; i++)
 231         {
 232             copy_rvec(xt[i], x[i]);
 233         }
 234         sfree(xt);
 235     }
 236 }
 237
 238 static void update_topol(const char *topinout, int p_num, int n_num,
 239                          const char *p_name, const char *n_name, char *grpname)
 240 {
 241 #define TEMP_FILENM "temp.top"
 242     FILE    *fpin, *fpout;
 243     char     buf[STRLEN], buf2[STRLEN], *temp, **mol_line = NULL;
 244     int      line, i, nsol, nmol_line, sol_line, nsol_last;
 245     gmx_bool bMolecules;
 246
 247     printf("\nProcessing topology\n");
 248     fpin  = gmx_ffopen(topinout, "r");
 249     fpout = gmx_ffopen(TEMP_FILENM, "w");
 250
 251     line       = 0;
 252     bMolecules = FALSE;
 253     nmol_line  = 0;
 254     sol_line   = -1;
 255     nsol_last  = -1;
 256     while (fgets(buf, STRLEN, fpin))
 257     {
 258         line++;
 259         strcpy(buf2, buf);
 260         if ((temp = strchr(buf2, '\n')) != NULL)
 261         {
 262             temp[0] = '\0';
 263         }
 264         ltrim(buf2);
 265         if (buf2[0] == '[')
 266         {
 267             buf2[0] = ' ';
 268             if ((temp = strchr(buf2, '\n')) != NULL)
 269             {
 270                 temp[0] = '\0';
 271             }
 272             rtrim(buf2);
 273             if (buf2[strlen(buf2)-1] == ']')
 274             {
 275                 buf2[strlen(buf2)-1] = '\0';
 276                 ltrim(buf2);
 277                 rtrim(buf2);
 278                 bMolecules = (gmx_strcasecmp(buf2, "molecules") == 0);
 279             }
 280             fprintf(fpout, "%s", buf);
 281         }
 282         else if (!bMolecules)
 283         {
 284             fprintf(fpout, "%s", buf);
 285         }
 286         else
 287         {
 288             /* Check if this is a line with solvent molecules */
 289             sscanf(buf, "%s", buf2);
 290             if (gmx_strcasecmp(buf2, grpname) == 0)
 291             {
 292                 sol_line = nmol_line;
 293                 sscanf(buf, "%*s %d", &nsol_last);
 294             }
 295             /* Store this molecules section line */
 296             srenew(mol_line, nmol_line+1);
 297             mol_line[nmol_line] = strdup(buf);
 298             nmol_line++;
 299         }
 300     }
 301     gmx_ffclose(fpin);
 302
 303     if (sol_line == -1)
 304     {
 305         gmx_ffclose(fpout);
 306         gmx_fatal(FARGS, "No line with moleculetype '%s' found the [ molecules ] section of file '%s'", grpname, topinout);
 307     }
 308     if (nsol_last < p_num+n_num)
 309     {
 310         gmx_ffclose(fpout);
 311         gmx_fatal(FARGS, "The last entry for moleculetype '%s' in the [ molecules ] section of file '%s' has less solvent molecules (%d) than were replaced (%d)", grpname, topinout, nsol_last, p_num+n_num);
 312     }
 313
 314     /* Print all the molecule entries */
 315     for (i = 0; i < nmol_line; i++)
 316     {
 317         if (i != sol_line)
 318         {
 319             fprintf(fpout, "%s", mol_line[i]);
 320         }
 321         else
 322         {
 323             printf("Replacing %d solute molecules in topology file (%s) "
 324                    " by %d %s and %d %s ions.\n",
 325                    p_num+n_num, topinout, p_num, p_name, n_num, n_name);
 326             nsol_last -= p_num + n_num;
 327             if (nsol_last > 0)
 328             {
 329                 fprintf(fpout, "%-10s  %d\n", grpname, nsol_last);
 330             }
 331             if (p_num > 0)
 332             {
 333                 fprintf(fpout, "%-15s  %d\n", p_name, p_num);
 334             }
 335             if (n_num > 0)
 336             {
 337                 fprintf(fpout, "%-15s  %d\n", n_name, n_num);
 338             }
 339         }
 340     }
 341     gmx_ffclose(fpout);
 342     /* use gmx_ffopen to generate backup of topinout */
 343     fpout = gmx_ffopen(topinout, "w");
 344     gmx_ffclose(fpout);
 345     rename(TEMP_FILENM, topinout);
 346 #undef TEMP_FILENM
 347 }
 348
 349 int gmx_genion(int argc, char *argv[])
 350 {
 351     const char        *desc[] = {
 352         "[THISMODULE] randomly replaces solvent molecules with monoatomic ions.",
 353         "The group of solvent molecules should be continuous and all molecules",
 354         "should have the same number of atoms.",
 355         "The user should add the ion molecules to the topology file or use",
 356         "the [TT]-p[tt] option to automatically modify the topology.[PAR]",
 357         "The ion molecule type, residue and atom names in all force fields",
 358         "are the capitalized element names without sign. This molecule name",
 359         "should be given with [TT]-pname[tt] or [TT]-nname[tt], and the",
 360         "[TT][molecules][tt] section of your topology updated accordingly,",
 361         "either by hand or with [TT]-p[tt]. Do not use an atom name instead!",
 362         "[PAR]Ions which can have multiple charge states get the multiplicity",
 363         "added, without sign, for the uncommon states only.[PAR]",
 364         "For larger ions, e.g. sulfate we recommended using [gmx-insert-molecules]."
 365     };
 366     const char        *bugs[] = {
 367         "If you specify a salt concentration existing ions are not taken into "
 368         "account. In effect you therefore specify the amount of salt to be added.",
 369     };
 370     static int         p_num    = 0, n_num = 0, p_q = 1, n_q = -1;
 371     static const char *p_name   = "NA", *n_name = "CL";
 372     static real        rmin     = 0.6, conc = 0;
 373     static int         seed     = 1993;
 374     static gmx_bool    bNeutral = FALSE;
 375     static t_pargs     pa[]     = {
 376         { "-np",    FALSE, etINT,  {&p_num}, "Number of positive ions"       },
 377         { "-pname", FALSE, etSTR,  {&p_name}, "Name of the positive ion"      },
 378         { "-pq",    FALSE, etINT,  {&p_q},   "Charge of the positive ion"    },
 379         { "-nn",    FALSE, etINT,  {&n_num}, "Number of negative ions"       },
 380         { "-nname", FALSE, etSTR,  {&n_name}, "Name of the negative ion"      },
 381         { "-nq",    FALSE, etINT,  {&n_q},   "Charge of the negative ion"    },
 382         { "-rmin",  FALSE, etREAL, {&rmin},  "Minimum distance between ions" },
 383         { "-seed",  FALSE, etINT,  {&seed},  "Seed for random number generator" },
 384         { "-conc",  FALSE, etREAL, {&conc},
 385           "Specify salt concentration (mol/liter). This will add sufficient ions to reach up to the specified concentration as computed from the volume of the cell in the input [TT].tpr[tt] file. Overrides the [TT]-np[tt] and [TT]-nn[tt] options." },
 386         { "-neutral", FALSE, etBOOL, {&bNeutral}, "This option will add enough ions to neutralize the system. These ions are added on top of those specified with [TT]-np[tt]/[TT]-nn[tt] or [TT]-conc[tt]. "}
 387     };
 388     t_topology         top;
 389     rvec              *x, *v;
 390     real               vol, qtot;
 391     matrix             box;
 392     t_atoms            atoms;
 393     t_pbc              pbc;
 394     int               *repl, ePBC;
 395     atom_id           *index;
 396     char              *grpname, title[STRLEN];
 397     gmx_bool          *bSet;
 398     int                i, nw, nwa, nsa, nsalt, iqtot;
 399     output_env_t       oenv;
 400     gmx_rng_t          rng;
 401     t_filenm           fnm[] = {
 402         { efTPX, NULL,  NULL,      ffREAD  },
 403         { efNDX, NULL,  NULL,      ffOPTRD },
 404         { efSTO, "-o",  NULL,      ffWRITE },
 405         { efTOP, "-p",  "topol",   ffOPTRW }
 406     };
 407 #define NFILE asize(fnm)
 408
 409     if (!parse_common_args(&argc, argv, PCA_BE_NICE, NFILE, fnm, asize(pa), pa,
 410                            asize(desc), desc, asize(bugs), bugs, &oenv))
 411     {
 412         return 0;
 413     }
 414
 415     /* Check input for something sensible */
 416     if ((p_num < 0) || (n_num < 0))
 417     {
 418         gmx_fatal(FARGS, "Negative number of ions to add?");
 419     }
 420
 421     if (conc > 0 && (p_num > 0 || n_num > 0))
 422     {
 423         fprintf(stderr, "WARNING: -conc specified, overriding -nn and -np.\n");
 424     }
 425
 426     /* Read atom positions and charges */
 427     read_tps_conf(ftp2fn(efTPX, NFILE, fnm), title, &top, &ePBC, &x, &v, box, FALSE);
 428     atoms = top.atoms;
 429
 430     /* Compute total charge */
 431     qtot = 0;
 432     for (i = 0; (i < atoms.nr); i++)
 433     {
 434         qtot += atoms.atom[i].q;
 435     }
 436     iqtot = gmx_nint(qtot);
 437
 438
 439     if (conc > 0)
 440     {
 441         /* Compute number of ions to be added */
 442         vol   = det(box);
 443         nsalt = gmx_nint(conc*vol*AVOGADRO/1e24);
 444         p_num = abs(nsalt*n_q);
 445         n_num = abs(nsalt*p_q);
 446     }
 447     if (bNeutral)
 448     {
 449         int qdelta = p_num*p_q + n_num*n_q + iqtot;
 450
 451         /* Check if the system is neutralizable
 452          * is (qdelta == p_q*p_num + n_q*n_num) solvable for p_num and n_num? */
 453         int gcd = gmx_greatest_common_divisor(n_q, p_q);
 454         if ((qdelta % gcd) != 0)
 455         {
 456             gmx_fatal(FARGS, "Can't neutralize this system using -nq %d and"
 457                       " -pq %d.\n", n_q, p_q);
 458         }
 459
 460         while (qdelta != 0)
 461         {
 462             while (qdelta < 0)
 463             {
 464                 p_num++;
 465                 qdelta += p_q;
 466             }
 467             while (qdelta > 0)
 468             {
 469                 n_num++;
 470                 qdelta += n_q;
 471             }
 472         }
 473     }
 474
 475     if ((p_num == 0) && (n_num == 0))
 476     {
 477         fprintf(stderr, "No ions to add, will just copy input configuration.\n");
 478     }
 479     else
 480     {
 481         printf("Will try to add %d %s ions and %d %s ions.\n",
 482                p_num, p_name, n_num, n_name);
 483         printf("Select a continuous group of solvent molecules\n");
 484         get_index(&atoms, ftp2fn_null(efNDX, NFILE, fnm), 1, &nwa, &index, &grpname);
 485         for (i = 1; i < nwa; i++)
 486         {
 487             if (index[i] != index[i-1]+1)
 488             {
 489                 gmx_fatal(FARGS, "The solvent group %s is not continuous: "
 490                           "index[%d]=%d, index[%d]=%d",
 491                           grpname, i, index[i-1]+1, i+1, index[i]+1);
 492             }
 493         }
 494         nsa = 1;
 495         while ((nsa < nwa) &&
 496                (atoms.atom[index[nsa]].resind ==
 497                 atoms.atom[index[nsa-1]].resind))
 498         {
 499             nsa++;
 500         }
 501         if (nwa % nsa)
 502         {
 503             gmx_fatal(FARGS, "Your solvent group size (%d) is not a multiple of %d",
 504                       nwa, nsa);
 505         }
 506         nw = nwa/nsa;
 507         fprintf(stderr, "Number of (%d-atomic) solvent molecules: %d\n", nsa, nw);
 508         if (p_num+n_num > nw)
 509         {
 510             gmx_fatal(FARGS, "Not enough solvent for adding ions");
 511         }
 512
 513         if (opt2bSet("-p", NFILE, fnm))
 514         {
 515             update_topol(opt2fn("-p", NFILE, fnm), p_num, n_num, p_name, n_name, grpname);
 516         }
 517
 518         snew(bSet, nw);
 519         snew(repl, nw);
 520
 521         snew(v, atoms.nr);
 522         snew(atoms.pdbinfo, atoms.nr);
 523
 524         set_pbc(&pbc, ePBC, box);
 525
 526         if (seed == 0)
 527         {
 528             rng = gmx_rng_init(gmx_rng_make_seed());
 529         }
 530         else
 531         {
 532             rng = gmx_rng_init(seed);
 533         }
 534         /* Now loop over the ions that have to be placed */
 535         while (p_num-- > 0)
 536         {
 537             insert_ion(nsa, &nw, bSet, repl, index, x, &pbc,
 538                        1, p_q, p_name, &atoms, rmin, rng);
 539         }
 540         while (n_num-- > 0)
 541         {
 542             insert_ion(nsa, &nw, bSet, repl, index, x, &pbc,
 543                        -1, n_q, n_name, &atoms, rmin, rng);
 544         }
 545         gmx_rng_destroy(rng);
 546         fprintf(stderr, "\n");
 547
 548         if (nw)
 549         {
 550             sort_ions(nsa, nw, repl, index, &atoms, x, p_name, n_name);
 551         }
 552
 553         sfree(atoms.pdbinfo);
 554         atoms.pdbinfo = NULL;
 555     }
 556     write_sto_conf(ftp2fn(efSTO, NFILE, fnm), *top.name, &atoms, x, NULL, ePBC, box);
 557
 558     return 0;
 559 }