src/gromacs/gmxpreprocess/genion.cpp

   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,2015,2016,2017,2018,2019, 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 "genion.h"
  40
  41 #include <cctype>
  42 #include <cmath>
  43 #include <cstdlib>
  44 #include <cstring>
  45
  46 #include "gromacs/commandline/pargs.h"
  47 #include "gromacs/fileio/confio.h"
  48 #include "gromacs/math/units.h"
  49 #include "gromacs/math/utilities.h"
  50 #include "gromacs/math/vec.h"
  51 #include "gromacs/mdlib/force.h"
  52 #include "gromacs/pbcutil/pbc.h"
  53 #include "gromacs/random/threefry.h"
  54 #include "gromacs/random/uniformintdistribution.h"
  55 #include "gromacs/topology/index.h"
  56 #include "gromacs/topology/topology.h"
  57 #include "gromacs/utility/arraysize.h"
  58 #include "gromacs/utility/cstringutil.h"
  59 #include "gromacs/utility/fatalerror.h"
  60 #include "gromacs/utility/futil.h"
  61 #include "gromacs/utility/smalloc.h"
  62
  63 static void insert_ion(int nsa, const int *nwater,
  64                        gmx_bool bSet[], int repl[], int index[],
  65                        rvec x[], t_pbc *pbc,
  66                        int sign, int q, const char *ionname,
  67                        t_atoms *atoms,
  68                        real rmin,
  69                        gmx::DefaultRandomEngine * rng)
  70 {
  71     int                                i, ei, nw;
  72     real                               rmin2;
  73     rvec                               dx;
  74     int64_t                            maxrand;
  75     gmx::UniformIntDistribution<int>   dist(0, *nwater-1);
  76
  77     nw       = *nwater;
  78     maxrand  = nw;
  79     maxrand *= 1000;
  80
  81     do
  82     {
  83         ei = dist(*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 = gmx_strdup(mname);
 130     i   = std::strlen(str)-1;
 131     while (i > 1 && (std::isdigit(str[i]) || (str[i] == '+') || (str[i] == '-')))
 132     {
 133         str[i] = '\0';
 134         i--;
 135     }
 136
 137     return str;
 138 }
 139
 140 static void sort_ions(int nsa, int nw, const int repl[], const int 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 = nullptr, **nptr = nullptr, **paptr = nullptr, **naptr = nullptr;
 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] = gmx_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] = gmx_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     FILE    *fpin, *fpout;
 242     char     buf[STRLEN], buf2[STRLEN], *temp, **mol_line = nullptr;
 243     int      line, i, nmol_line, sol_line, nsol_last;
 244     gmx_bool bMolecules;
 245     char     temporary_filename[STRLEN];
 246
 247     printf("\nProcessing topology\n");
 248     fpin  = gmx_ffopen(topinout, "r");
 249     std::strncpy(temporary_filename, "temp.topXXXXXX", STRLEN);
 250     fpout = gmx_fopen_temporary(temporary_filename);
 251
 252     line       = 0;
 253     bMolecules = FALSE;
 254     nmol_line  = 0;
 255     sol_line   = -1;
 256     nsol_last  = -1;
 257     while (fgets(buf, STRLEN, fpin))
 258     {
 259         line++;
 260         std::strcpy(buf2, buf);
 261         if ((temp = std::strchr(buf2, '\n')) != nullptr)
 262         {
 263             temp[0] = '\0';
 264         }
 265         ltrim(buf2);
 266         if (buf2[0] == '[')
 267         {
 268             buf2[0] = ' ';
 269             if ((temp = std::strchr(buf2, '\n')) != nullptr)
 270             {
 271                 temp[0] = '\0';
 272             }
 273             rtrim(buf2);
 274             if (buf2[std::strlen(buf2)-1] == ']')
 275             {
 276                 buf2[std::strlen(buf2)-1] = '\0';
 277                 ltrim(buf2);
 278                 rtrim(buf2);
 279                 bMolecules = (gmx_strcasecmp(buf2, "molecules") == 0);
 280             }
 281             fprintf(fpout, "%s", buf);
 282         }
 283         else if (!bMolecules)
 284         {
 285             fprintf(fpout, "%s", buf);
 286         }
 287         else
 288         {
 289             /* Check if this is a line with solvent molecules */
 290             sscanf(buf, "%s", buf2);
 291             if (gmx_strcasecmp(buf2, grpname) == 0)
 292             {
 293                 sol_line = nmol_line;
 294                 sscanf(buf, "%*s %d", &nsol_last);
 295             }
 296             /* Store this molecules section line */
 297             srenew(mol_line, nmol_line+1);
 298             mol_line[nmol_line] = gmx_strdup(buf);
 299             nmol_line++;
 300         }
 301     }
 302     gmx_ffclose(fpin);
 303
 304     if (sol_line == -1)
 305     {
 306         gmx_ffclose(fpout);
 307         gmx_fatal(FARGS, "No line with moleculetype '%s' found the [ molecules ] section of file '%s'", grpname, topinout);
 308     }
 309     if (nsol_last < p_num+n_num)
 310     {
 311         gmx_ffclose(fpout);
 312         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);
 313     }
 314
 315     /* Print all the molecule entries */
 316     for (i = 0; i < nmol_line; i++)
 317     {
 318         if (i != sol_line)
 319         {
 320             fprintf(fpout, "%s", mol_line[i]);
 321         }
 322         else
 323         {
 324             printf("Replacing %d solute molecules in topology file (%s) "
 325                    " by %d %s and %d %s ions.\n",
 326                    p_num+n_num, topinout, p_num, p_name, n_num, n_name);
 327             nsol_last -= p_num + n_num;
 328             if (nsol_last > 0)
 329             {
 330                 fprintf(fpout, "%-10s  %d\n", grpname, nsol_last);
 331             }
 332             if (p_num > 0)
 333             {
 334                 fprintf(fpout, "%-15s  %d\n", p_name, p_num);
 335             }
 336             if (n_num > 0)
 337             {
 338                 fprintf(fpout, "%-15s  %d\n", n_name, n_num);
 339             }
 340         }
 341     }
 342     gmx_ffclose(fpout);
 343     make_backup(topinout);
 344     gmx_file_rename(temporary_filename, topinout);
 345 }
 346
 347 int gmx_genion(int argc, char *argv[])
 348 {
 349     const char        *desc[] = {
 350         "[THISMODULE] randomly replaces solvent molecules with monoatomic ions.",
 351         "The group of solvent molecules should be continuous and all molecules",
 352         "should have the same number of atoms.",
 353         "The user should add the ion molecules to the topology file or use",
 354         "the [TT]-p[tt] option to automatically modify the topology.[PAR]",
 355         "The ion molecule type, residue and atom names in all force fields",
 356         "are the capitalized element names without sign. This molecule name",
 357         "should be given with [TT]-pname[tt] or [TT]-nname[tt], and the",
 358         "[TT][molecules][tt] section of your topology updated accordingly,",
 359         "either by hand or with [TT]-p[tt]. Do not use an atom name instead!",
 360         "[PAR]Ions which can have multiple charge states get the multiplicity",
 361         "added, without sign, for the uncommon states only.[PAR]",
 362         "For larger ions, e.g. sulfate we recommended using [gmx-insert-molecules]."
 363     };
 364     const char        *bugs[] = {
 365         "If you specify a salt concentration existing ions are not taken into "
 366         "account. In effect you therefore specify the amount of salt to be added.",
 367     };
 368     static int         p_num    = 0, n_num = 0, p_q = 1, n_q = -1;
 369     static const char *p_name   = "NA", *n_name = "CL";
 370     static real        rmin     = 0.6, conc = 0;
 371     static int         seed     = 0;
 372     static gmx_bool    bNeutral = FALSE;
 373     static t_pargs     pa[]     = {
 374         { "-np",    FALSE, etINT,  {&p_num}, "Number of positive ions"       },
 375         { "-pname", FALSE, etSTR,  {&p_name}, "Name of the positive ion"      },
 376         { "-pq",    FALSE, etINT,  {&p_q},   "Charge of the positive ion"    },
 377         { "-nn",    FALSE, etINT,  {&n_num}, "Number of negative ions"       },
 378         { "-nname", FALSE, etSTR,  {&n_name}, "Name of the negative ion"      },
 379         { "-nq",    FALSE, etINT,  {&n_q},   "Charge of the negative ion"    },
 380         { "-rmin",  FALSE, etREAL, {&rmin},  "Minimum distance between ions" },
 381         { "-seed",  FALSE, etINT,  {&seed},  "Seed for random number generator (0 means generate)" },
 382         { "-conc",  FALSE, etREAL, {&conc},
 383           "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 [REF].tpr[ref] file. Overrides the [TT]-np[tt] and [TT]-nn[tt] options." },
 384         { "-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]. "}
 385     };
 386     t_topology         top;
 387     rvec              *x, *v;
 388     real               vol;
 389     matrix             box;
 390     t_atoms            atoms;
 391     t_pbc              pbc;
 392     int               *repl, ePBC;
 393     int               *index;
 394     char              *grpname;
 395     gmx_bool          *bSet;
 396     int                i, nw, nwa, nsa, nsalt, iqtot;
 397     gmx_output_env_t  *oenv;
 398     t_filenm           fnm[] = {
 399         { efTPR, nullptr,  nullptr,      ffREAD  },
 400         { efNDX, nullptr,  nullptr,      ffOPTRD },
 401         { efSTO, "-o",  nullptr,      ffWRITE },
 402         { efTOP, "-p",  "topol",   ffOPTRW }
 403     };
 404 #define NFILE asize(fnm)
 405
 406     if (!parse_common_args(&argc, argv, 0, NFILE, fnm, asize(pa), pa,
 407                            asize(desc), desc, asize(bugs), bugs, &oenv))
 408     {
 409         return 0;
 410     }
 411
 412     /* Check input for something sensible */
 413     if ((p_num < 0) || (n_num < 0))
 414     {
 415         gmx_fatal(FARGS, "Negative number of ions to add?");
 416     }
 417
 418     if (conc > 0 && (p_num > 0 || n_num > 0))
 419     {
 420         fprintf(stderr, "WARNING: -conc specified, overriding -nn and -np.\n");
 421     }
 422
 423     /* Read atom positions and charges */
 424     read_tps_conf(ftp2fn(efTPR, NFILE, fnm), &top, &ePBC, &x, &v, box, FALSE);
 425     atoms = top.atoms;
 426
 427     /* Compute total charge */
 428     double qtot = 0;
 429     for (i = 0; (i < atoms.nr); i++)
 430     {
 431         qtot += atoms.atom[i].q;
 432     }
 433     iqtot = gmx::roundToInt(qtot);
 434
 435
 436     if (conc > 0)
 437     {
 438         /* Compute number of ions to be added */
 439         vol   = det(box);
 440         nsalt = gmx::roundToInt(conc*vol*AVOGADRO/1e24);
 441         p_num = abs(nsalt*n_q);
 442         n_num = abs(nsalt*p_q);
 443     }
 444     if (bNeutral)
 445     {
 446         int qdelta = p_num*p_q + n_num*n_q + iqtot;
 447
 448         /* Check if the system is neutralizable
 449          * is (qdelta == p_q*p_num + n_q*n_num) solvable for p_num and n_num? */
 450         int gcd = gmx_greatest_common_divisor(n_q, p_q);
 451         if ((qdelta % gcd) != 0)
 452         {
 453             gmx_fatal(FARGS, "Can't neutralize this system using -nq %d and"
 454                       " -pq %d.\n", n_q, p_q);
 455         }
 456
 457         while (qdelta != 0)
 458         {
 459             while (qdelta < 0)
 460             {
 461                 p_num++;
 462                 qdelta += p_q;
 463             }
 464             while (qdelta > 0)
 465             {
 466                 n_num++;
 467                 qdelta += n_q;
 468             }
 469         }
 470     }
 471
 472     if ((p_num == 0) && (n_num == 0))
 473     {
 474         fprintf(stderr, "No ions to add, will just copy input configuration.\n");
 475     }
 476     else
 477     {
 478         printf("Will try to add %d %s ions and %d %s ions.\n",
 479                p_num, p_name, n_num, n_name);
 480         printf("Select a continuous group of solvent molecules\n");
 481         get_index(&atoms, ftp2fn_null(efNDX, NFILE, fnm), 1, &nwa, &index, &grpname);
 482         for (i = 1; i < nwa; i++)
 483         {
 484             if (index[i] != index[i-1]+1)
 485             {
 486                 gmx_fatal(FARGS, "The solvent group %s is not continuous: "
 487                           "index[%d]=%d, index[%d]=%d",
 488                           grpname, i, index[i-1]+1, i+1, index[i]+1);
 489             }
 490         }
 491         nsa = 1;
 492         while ((nsa < nwa) &&
 493                (atoms.atom[index[nsa]].resind ==
 494                 atoms.atom[index[nsa-1]].resind))
 495         {
 496             nsa++;
 497         }
 498         if (nwa % nsa)
 499         {
 500             gmx_fatal(FARGS, "Your solvent group size (%d) is not a multiple of %d",
 501                       nwa, nsa);
 502         }
 503         nw = nwa/nsa;
 504         fprintf(stderr, "Number of (%d-atomic) solvent molecules: %d\n", nsa, nw);
 505         if (p_num+n_num > nw)
 506         {
 507             gmx_fatal(FARGS, "Not enough solvent for adding ions");
 508         }
 509
 510         if (opt2bSet("-p", NFILE, fnm))
 511         {
 512             update_topol(opt2fn("-p", NFILE, fnm), p_num, n_num, p_name, n_name, grpname);
 513         }
 514
 515         snew(bSet, nw);
 516         snew(repl, nw);
 517
 518         snew(v, atoms.nr);
 519         snew(atoms.pdbinfo, atoms.nr);
 520
 521         set_pbc(&pbc, ePBC, box);
 522
 523
 524         if (seed == 0)
 525         {
 526             // For now we make do with 32 bits to avoid changing the user input to 64 bit hex
 527             seed = static_cast<int>(gmx::makeRandomSeed());
 528         }
 529         fprintf(stderr, "Using random seed %d.\n", seed);
 530
 531         gmx::DefaultRandomEngine rng(seed);
 532
 533         /* Now loop over the ions that have to be placed */
 534         while (p_num-- > 0)
 535         {
 536             insert_ion(nsa, &nw, bSet, repl, index, x, &pbc,
 537                        1, p_q, p_name, &atoms, rmin, &rng);
 538         }
 539         while (n_num-- > 0)
 540         {
 541             insert_ion(nsa, &nw, bSet, repl, index, x, &pbc,
 542                        -1, n_q, n_name, &atoms, rmin, &rng);
 543         }
 544         fprintf(stderr, "\n");
 545
 546         if (nw)
 547         {
 548             sort_ions(nsa, nw, repl, index, &atoms, x, p_name, n_name);
 549         }
 550
 551         sfree(atoms.pdbinfo);
 552         atoms.pdbinfo = nullptr;
 553     }
 554     write_sto_conf(ftp2fn(efSTO, NFILE, fnm), *top.name, &atoms, x, nullptr, ePBC, box);
 555
 556     return 0;
 557 }