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