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 #include "gmxpre.h"
  38
  39 #include "config.h"
  40
  41 #include <ctype.h>
  42 #include <stdlib.h>
  43 #include <string.h>
  44
  45 #include "gromacs/utility/cstringutil.h"
  46 #include "gromacs/utility/smalloc.h"
  47 #include "gromacs/fileio/confio.h"
  48 #include "gromacs/commandline/pargs.h"
  49 #include "gromacs/pbcutil/pbc.h"
  50 #include "gromacs/legacyheaders/force.h"
  51 #include "gromacs/utility/fatalerror.h"
  52 #include "gromacs/utility/futil.h"
  53 #include "gromacs/math/utilities.h"
  54 #include "gromacs/legacyheaders/macros.h"
  55 #include "gromacs/math/vec.h"
  56 #include "gromacs/fileio/tpxio.h"
  57 #include "gromacs/legacyheaders/mdrun.h"
  58 #include "gromacs/random/random.h"
  59 #include "gromacs/topology/index.h"
  60 #include "gmx_ana.h"
  61
  62 static void insert_ion(int nsa, int *nwater,
  63                        gmx_bool bSet[], int repl[], atom_id index[],
  64                        rvec x[], t_pbc *pbc,
  65                        int sign, int q, const char *ionname,
  66                        t_atoms *atoms,
  67                        real rmin, gmx_rng_t rng)
  68 {
  69     int             i, ei, nw;
  70     real            rmin2;
  71     rvec            dx;
  72     gmx_int64_t     maxrand;
  73
  74     ei       = -1;
  75     nw       = *nwater;
  76     maxrand  = nw;
  77     maxrand *= 1000;
  78
  79     do
  80     {
  81         ei = nw*gmx_rng_uniform_real(rng);
  82         maxrand--;
  83     }
  84     while (bSet[ei] && (maxrand > 0));
  85     if (bSet[ei])
  86     {
  87         gmx_fatal(FARGS, "No more replaceable solvent!");
  88     }
  89
  90     fprintf(stderr, "Replacing solvent molecule %d (atom %d) with %s\n",
  91             ei, index[nsa*ei], ionname);
  92
  93     /* Replace solvent molecule charges with ion charge */
  94     bSet[ei] = TRUE;
  95     repl[ei] = sign;
  96
  97     atoms->atom[index[nsa*ei]].q = q;
  98     for (i = 1; i < nsa; i++)
  99     {
 100         atoms->atom[index[nsa*ei+i]].q = 0;
 101     }
 102
 103     /* Mark all solvent molecules within rmin as unavailable for substitution */
 104     if (rmin > 0)
 105     {
 106         rmin2 = rmin*rmin;
 107         for (i = 0; (i < nw); i++)
 108         {
 109             if (!bSet[i])
 110             {
 111                 pbc_dx(pbc, x[index[nsa*ei]], x[index[nsa*i]], dx);
 112                 if (iprod(dx, dx) < rmin2)
 113                 {
 114                     bSet[i] = TRUE;
 115                 }
 116             }
 117         }
 118     }
 119 }
 120
 121
 122 static char *aname(const char *mname)
 123 {
 124     char *str;
 125     int   i;
 126
 127     str = gmx_strdup(mname);
 128     i   = strlen(str)-1;
 129     while (i > 1 && (isdigit(str[i]) || (str[i] == '+') || (str[i] == '-')))
 130     {
 131         str[i] = '\0';
 132         i--;
 133     }
 134
 135     return str;
 136 }
 137
 138 void sort_ions(int nsa, int nw, int repl[], atom_id index[],
 139                t_atoms *atoms, rvec x[],
 140                const char *p_name, const char *n_name)
 141 {
 142     int    i, j, k, r, np, nn, starta, startr, npi, nni;
 143     rvec  *xt;
 144     char **pptr = NULL, **nptr = NULL, **paptr = NULL, **naptr = NULL;
 145
 146     snew(xt, atoms->nr);
 147
 148     /* Put all the solvent in front and count the added ions */
 149     np = 0;
 150     nn = 0;
 151     j  = index[0];
 152     for (i = 0; i < nw; i++)
 153     {
 154         r = repl[i];
 155         if (r == 0)
 156         {
 157             for (k = 0; k < nsa; k++)
 158             {
 159                 copy_rvec(x[index[nsa*i+k]], xt[j++]);
 160             }
 161         }
 162         else if (r > 0)
 163         {
 164             np++;
 165         }
 166         else if (r < 0)
 167         {
 168             nn++;
 169         }
 170     }
 171
 172     if (np+nn > 0)
 173     {
 174         /* Put the positive and negative ions at the end */
 175         starta = index[nsa*(nw - np - nn)];
 176         startr = atoms->atom[starta].resind;
 177
 178         if (np)
 179         {
 180             snew(pptr, 1);
 181             pptr[0] = gmx_strdup(p_name);
 182             snew(paptr, 1);
 183             paptr[0] = aname(p_name);
 184         }
 185         if (nn)
 186         {
 187             snew(nptr, 1);
 188             nptr[0] = gmx_strdup(n_name);
 189             snew(naptr, 1);
 190             naptr[0] = aname(n_name);
 191         }
 192         npi = 0;
 193         nni = 0;
 194         for (i = 0; i < nw; i++)
 195         {
 196             r = repl[i];
 197             if (r > 0)
 198             {
 199                 j = starta+npi;
 200                 k = startr+npi;
 201                 copy_rvec(x[index[nsa*i]], xt[j]);
 202                 atoms->atomname[j]     = paptr;
 203                 atoms->atom[j].resind  = k;
 204                 atoms->resinfo[k].name = pptr;
 205                 npi++;
 206             }
 207             else if (r < 0)
 208             {
 209                 j = starta+np+nni;
 210                 k = startr+np+nni;
 211                 copy_rvec(x[index[nsa*i]], xt[j]);
 212                 atoms->atomname[j]     = naptr;
 213                 atoms->atom[j].resind  = k;
 214                 atoms->resinfo[k].name = nptr;
 215                 nni++;
 216             }
 217         }
 218         for (i = index[nsa*nw-1]+1; i < atoms->nr; i++)
 219         {
 220             j                  = i-(nsa-1)*(np+nn);
 221             atoms->atomname[j] = atoms->atomname[i];
 222             atoms->atom[j]     = atoms->atom[i];
 223             copy_rvec(x[i], xt[j]);
 224         }
 225         atoms->nr -= (nsa-1)*(np+nn);
 226
 227         /* Copy the new positions back */
 228         for (i = index[0]; i < atoms->nr; i++)
 229         {
 230             copy_rvec(xt[i], x[i]);
 231         }
 232         sfree(xt);
 233     }
 234 }
 235
 236 static void update_topol(const char *topinout, int p_num, int n_num,
 237                          const char *p_name, const char *n_name, char *grpname)
 238 {
 239 #define TEMP_FILENM "temp.top"
 240     FILE    *fpin, *fpout;
 241     char     buf[STRLEN], buf2[STRLEN], *temp, **mol_line = NULL;
 242     int      line, i, nsol, nmol_line, sol_line, nsol_last;
 243     gmx_bool bMolecules;
 244
 245     printf("\nProcessing topology\n");
 246     fpin  = gmx_ffopen(topinout, "r");
 247     fpout = gmx_ffopen(TEMP_FILENM, "w");
 248
 249     line       = 0;
 250     bMolecules = FALSE;
 251     nmol_line  = 0;
 252     sol_line   = -1;
 253     nsol_last  = -1;
 254     while (fgets(buf, STRLEN, fpin))
 255     {
 256         line++;
 257         strcpy(buf2, buf);
 258         if ((temp = strchr(buf2, '\n')) != NULL)
 259         {
 260             temp[0] = '\0';
 261         }
 262         ltrim(buf2);
 263         if (buf2[0] == '[')
 264         {
 265             buf2[0] = ' ';
 266             if ((temp = strchr(buf2, '\n')) != NULL)
 267             {
 268                 temp[0] = '\0';
 269             }
 270             rtrim(buf2);
 271             if (buf2[strlen(buf2)-1] == ']')
 272             {
 273                 buf2[strlen(buf2)-1] = '\0';
 274                 ltrim(buf2);
 275                 rtrim(buf2);
 276                 bMolecules = (gmx_strcasecmp(buf2, "molecules") == 0);
 277             }
 278             fprintf(fpout, "%s", buf);
 279         }
 280         else if (!bMolecules)
 281         {
 282             fprintf(fpout, "%s", buf);
 283         }
 284         else
 285         {
 286             /* Check if this is a line with solvent molecules */
 287             sscanf(buf, "%s", buf2);
 288             if (gmx_strcasecmp(buf2, grpname) == 0)
 289             {
 290                 sol_line = nmol_line;
 291                 sscanf(buf, "%*s %d", &nsol_last);
 292             }
 293             /* Store this molecules section line */
 294             srenew(mol_line, nmol_line+1);
 295             mol_line[nmol_line] = gmx_strdup(buf);
 296             nmol_line++;
 297         }
 298     }
 299     gmx_ffclose(fpin);
 300
 301     if (sol_line == -1)
 302     {
 303         gmx_ffclose(fpout);
 304         gmx_fatal(FARGS, "No line with moleculetype '%s' found the [ molecules ] section of file '%s'", grpname, topinout);
 305     }
 306     if (nsol_last < p_num+n_num)
 307     {
 308         gmx_ffclose(fpout);
 309         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);
 310     }
 311
 312     /* Print all the molecule entries */
 313     for (i = 0; i < nmol_line; i++)
 314     {
 315         if (i != sol_line)
 316         {
 317             fprintf(fpout, "%s", mol_line[i]);
 318         }
 319         else
 320         {
 321             printf("Replacing %d solute molecules in topology file (%s) "
 322                    " by %d %s and %d %s ions.\n",
 323                    p_num+n_num, topinout, p_num, p_name, n_num, n_name);
 324             nsol_last -= p_num + n_num;
 325             if (nsol_last > 0)
 326             {
 327                 fprintf(fpout, "%-10s  %d\n", grpname, nsol_last);
 328             }
 329             if (p_num > 0)
 330             {
 331                 fprintf(fpout, "%-15s  %d\n", p_name, p_num);
 332             }
 333             if (n_num > 0)
 334             {
 335                 fprintf(fpout, "%-15s  %d\n", n_name, n_num);
 336             }
 337         }
 338     }
 339     gmx_ffclose(fpout);
 340     /* use gmx_ffopen to generate backup of topinout */
 341     fpout = gmx_ffopen(topinout, "w");
 342     gmx_ffclose(fpout);
 343     rename(TEMP_FILENM, topinout);
 344 #undef TEMP_FILENM
 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     = 1993;
 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" },
 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 [TT].tpr[tt] 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, qtot;
 389     matrix             box;
 390     t_atoms            atoms;
 391     t_pbc              pbc;
 392     int               *repl, ePBC;
 393     atom_id           *index;
 394     char              *grpname, title[STRLEN];
 395     gmx_bool          *bSet;
 396     int                i, nw, nwa, nsa, nsalt, iqtot;
 397     output_env_t       oenv;
 398     gmx_rng_t          rng;
 399     t_filenm           fnm[] = {
 400         { efTPX, NULL,  NULL,      ffREAD  },
 401         { efNDX, NULL,  NULL,      ffOPTRD },
 402         { efSTO, "-o",  NULL,      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(efTPX, NFILE, fnm), title, &top, &ePBC, &x, &v, box, FALSE);
 426     atoms = top.atoms;
 427
 428     /* Compute total charge */
 429     qtot = 0;
 430     for (i = 0; (i < atoms.nr); i++)
 431     {
 432         qtot += atoms.atom[i].q;
 433     }
 434     iqtot = gmx_nint(qtot);
 435
 436
 437     if (conc > 0)
 438     {
 439         /* Compute number of ions to be added */
 440         vol   = det(box);
 441         nsalt = gmx_nint(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         if (seed == 0)
 525         {
 526             rng = gmx_rng_init(gmx_rng_make_seed());
 527         }
 528         else
 529         {
 530             rng = gmx_rng_init(seed);
 531         }
 532         /* Now loop over the ions that have to be placed */
 533         while (p_num-- > 0)
 534         {
 535             insert_ion(nsa, &nw, bSet, repl, index, x, &pbc,
 536                        1, p_q, p_name, &atoms, rmin, rng);
 537         }
 538         while (n_num-- > 0)
 539         {
 540             insert_ion(nsa, &nw, bSet, repl, index, x, &pbc,
 541                        -1, n_q, n_name, &atoms, rmin, rng);
 542         }
 543         gmx_rng_destroy(rng);
 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 = NULL;
 553     }
 554     write_sto_conf(ftp2fn(efSTO, NFILE, fnm), *top.name, &atoms, x, NULL, ePBC, box);
 555
 556     return 0;
 557 }