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 "gromacs/utility/fatalerror.h"
  52 #include "gromacs/utility/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 "gromacs/random/random.h"
  59 #include "index.h"
  60 #include "mtop_util.h"
  61 #include "gmx_ana.h"
  62
  63 static void insert_ion(int nsa, int *nwater,
  64                        gmx_bool bSet[], int repl[], atom_id index[],
  65                        rvec x[], t_pbc *pbc,
  66                        int sign, int q, const char *ionname,
  67                        t_atoms *atoms,
  68                        real rmin, gmx_rng_t rng)
  69 {
  70     int             i, ei, nw;
  71     real            rmin2;
  72     rvec            dx;
  73     gmx_int64_t     maxrand;
  74
  75     ei       = -1;
  76     nw       = *nwater;
  77     maxrand  = nw;
  78     maxrand *= 1000;
  79
  80     do
  81     {
  82         ei = nw*gmx_rng_uniform_real(rng);
  83         maxrand--;
  84     }
  85     while (bSet[ei] && (maxrand > 0));
  86     if (bSet[ei])
  87     {
  88         gmx_fatal(FARGS, "No more replaceable solvent!");
  89     }
  90
  91     fprintf(stderr, "Replacing solvent molecule %d (atom %d) with %s\n",
  92             ei, index[nsa*ei], ionname);
  93
  94     /* Replace solvent molecule charges with ion charge */
  95     bSet[ei] = TRUE;
  96     repl[ei] = sign;
  97
  98     atoms->atom[index[nsa*ei]].q = q;
  99     for (i = 1; i < nsa; i++)
 100     {
 101         atoms->atom[index[nsa*ei+i]].q = 0;
 102     }
 103
 104     /* Mark all solvent molecules within rmin as unavailable for substitution */
 105     if (rmin > 0)
 106     {
 107         rmin2 = rmin*rmin;
 108         for (i = 0; (i < nw); i++)
 109         {
 110             if (!bSet[i])
 111             {
 112                 pbc_dx(pbc, x[index[nsa*ei]], x[index[nsa*i]], dx);
 113                 if (iprod(dx, dx) < rmin2)
 114                 {
 115                     bSet[i] = TRUE;
 116                 }
 117             }
 118         }
 119     }
 120 }
 121
 122
 123 static char *aname(const char *mname)
 124 {
 125     char *str;
 126     int   i;
 127
 128     str = strdup(mname);
 129     i   = strlen(str)-1;
 130     while (i > 1 && (isdigit(str[i]) || (str[i] == '+') || (str[i] == '-')))
 131     {
 132         str[i] = '\0';
 133         i--;
 134     }
 135
 136     return str;
 137 }
 138
 139 void sort_ions(int nsa, int nw, int repl[], atom_id index[],
 140                t_atoms *atoms, rvec x[],
 141                const char *p_name, const char *n_name)
 142 {
 143     int    i, j, k, r, np, nn, starta, startr, npi, nni;
 144     rvec  *xt;
 145     char **pptr = NULL, **nptr = NULL, **paptr = NULL, **naptr = NULL;
 146
 147     snew(xt, atoms->nr);
 148
 149     /* Put all the solvent in front and count the added ions */
 150     np = 0;
 151     nn = 0;
 152     j  = index[0];
 153     for (i = 0; i < nw; i++)
 154     {
 155         r = repl[i];
 156         if (r == 0)
 157         {
 158             for (k = 0; k < nsa; k++)
 159             {
 160                 copy_rvec(x[index[nsa*i+k]], xt[j++]);
 161             }
 162         }
 163         else if (r > 0)
 164         {
 165             np++;
 166         }
 167         else if (r < 0)
 168         {
 169             nn++;
 170         }
 171     }
 172
 173     if (np+nn > 0)
 174     {
 175         /* Put the positive and negative ions at the end */
 176         starta = index[nsa*(nw - np - nn)];
 177         startr = atoms->atom[starta].resind;
 178
 179         if (np)
 180         {
 181             snew(pptr, 1);
 182             pptr[0] = strdup(p_name);
 183             snew(paptr, 1);
 184             paptr[0] = aname(p_name);
 185         }
 186         if (nn)
 187         {
 188             snew(nptr, 1);
 189             nptr[0] = strdup(n_name);
 190             snew(naptr, 1);
 191             naptr[0] = aname(n_name);
 192         }
 193         npi = 0;
 194         nni = 0;
 195         for (i = 0; i < nw; i++)
 196         {
 197             r = repl[i];
 198             if (r > 0)
 199             {
 200                 j = starta+npi;
 201                 k = startr+npi;
 202                 copy_rvec(x[index[nsa*i]], xt[j]);
 203                 atoms->atomname[j]     = paptr;
 204                 atoms->atom[j].resind  = k;
 205                 atoms->resinfo[k].name = pptr;
 206                 npi++;
 207             }
 208             else if (r < 0)
 209             {
 210                 j = starta+np+nni;
 211                 k = startr+np+nni;
 212                 copy_rvec(x[index[nsa*i]], xt[j]);
 213                 atoms->atomname[j]     = naptr;
 214                 atoms->atom[j].resind  = k;
 215                 atoms->resinfo[k].name = nptr;
 216                 nni++;
 217             }
 218         }
 219         for (i = index[nsa*nw-1]+1; i < atoms->nr; i++)
 220         {
 221             j                  = i-(nsa-1)*(np+nn);
 222             atoms->atomname[j] = atoms->atomname[i];
 223             atoms->atom[j]     = atoms->atom[i];
 224             copy_rvec(x[i], xt[j]);
 225         }
 226         atoms->nr -= (nsa-1)*(np+nn);
 227
 228         /* Copy the new positions back */
 229         for (i = index[0]; i < atoms->nr; i++)
 230         {
 231             copy_rvec(xt[i], x[i]);
 232         }
 233         sfree(xt);
 234     }
 235 }
 236
 237 static void update_topol(const char *topinout, int p_num, int n_num,
 238                          const char *p_name, const char *n_name, char *grpname)
 239 {
 240 #define TEMP_FILENM "temp.top"
 241     FILE    *fpin, *fpout;
 242     char     buf[STRLEN], buf2[STRLEN], *temp, **mol_line = NULL;
 243     int      line, i, nsol, nmol_line, sol_line, nsol_last;
 244     gmx_bool bMolecules;
 245
 246     printf("\nProcessing topology\n");
 247     fpin  = gmx_ffopen(topinout, "r");
 248     fpout = gmx_ffopen(TEMP_FILENM, "w");
 249
 250     line       = 0;
 251     bMolecules = FALSE;
 252     nmol_line  = 0;
 253     sol_line   = -1;
 254     nsol_last  = -1;
 255     while (fgets(buf, STRLEN, fpin))
 256     {
 257         line++;
 258         strcpy(buf2, buf);
 259         if ((temp = strchr(buf2, '\n')) != NULL)
 260         {
 261             temp[0] = '\0';
 262         }
 263         ltrim(buf2);
 264         if (buf2[0] == '[')
 265         {
 266             buf2[0] = ' ';
 267             if ((temp = strchr(buf2, '\n')) != NULL)
 268             {
 269                 temp[0] = '\0';
 270             }
 271             rtrim(buf2);
 272             if (buf2[strlen(buf2)-1] == ']')
 273             {
 274                 buf2[strlen(buf2)-1] = '\0';
 275                 ltrim(buf2);
 276                 rtrim(buf2);
 277                 bMolecules = (gmx_strcasecmp(buf2, "molecules") == 0);
 278             }
 279             fprintf(fpout, "%s", buf);
 280         }
 281         else if (!bMolecules)
 282         {
 283             fprintf(fpout, "%s", buf);
 284         }
 285         else
 286         {
 287             /* Check if this is a line with solvent molecules */
 288             sscanf(buf, "%s", buf2);
 289             if (gmx_strcasecmp(buf2, grpname) == 0)
 290             {
 291                 sol_line = nmol_line;
 292                 sscanf(buf, "%*s %d", &nsol_last);
 293             }
 294             /* Store this molecules section line */
 295             srenew(mol_line, nmol_line+1);
 296             mol_line[nmol_line] = strdup(buf);
 297             nmol_line++;
 298         }
 299     }
 300     gmx_ffclose(fpin);
 301
 302     if (sol_line == -1)
 303     {
 304         gmx_ffclose(fpout);
 305         gmx_fatal(FARGS, "No line with moleculetype '%s' found the [ molecules ] section of file '%s'", grpname, topinout);
 306     }
 307     if (nsol_last < p_num+n_num)
 308     {
 309         gmx_ffclose(fpout);
 310         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);
 311     }
 312
 313     /* Print all the molecule entries */
 314     for (i = 0; i < nmol_line; i++)
 315     {
 316         if (i != sol_line)
 317         {
 318             fprintf(fpout, "%s", mol_line[i]);
 319         }
 320         else
 321         {
 322             printf("Replacing %d solute molecules in topology file (%s) "
 323                    " by %d %s and %d %s ions.\n",
 324                    p_num+n_num, topinout, p_num, p_name, n_num, n_name);
 325             nsol_last -= p_num + n_num;
 326             if (nsol_last > 0)
 327             {
 328                 fprintf(fpout, "%-10s  %d\n", grpname, nsol_last);
 329             }
 330             if (p_num > 0)
 331             {
 332                 fprintf(fpout, "%-15s  %d\n", p_name, p_num);
 333             }
 334             if (n_num > 0)
 335             {
 336                 fprintf(fpout, "%-15s  %d\n", n_name, n_num);
 337             }
 338         }
 339     }
 340     gmx_ffclose(fpout);
 341     /* use gmx_ffopen to generate backup of topinout */
 342     fpout = gmx_ffopen(topinout, "w");
 343     gmx_ffclose(fpout);
 344     rename(TEMP_FILENM, topinout);
 345 #undef TEMP_FILENM
 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     = 1993;
 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" },
 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 [TT].tpr[tt] 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, qtot;
 390     matrix             box;
 391     t_atoms            atoms;
 392     t_pbc              pbc;
 393     int               *repl, ePBC;
 394     atom_id           *index;
 395     char              *grpname, title[STRLEN];
 396     gmx_bool          *bSet;
 397     int                i, nw, nwa, nsa, nsalt, iqtot;
 398     output_env_t       oenv;
 399     gmx_rng_t          rng;
 400     t_filenm           fnm[] = {
 401         { efTPX, NULL,  NULL,      ffREAD  },
 402         { efNDX, NULL,  NULL,      ffOPTRD },
 403         { efSTO, "-o",  NULL,      ffWRITE },
 404         { efTOP, "-p",  "topol",   ffOPTRW }
 405     };
 406 #define NFILE asize(fnm)
 407
 408     if (!parse_common_args(&argc, argv, PCA_BE_NICE, NFILE, fnm, asize(pa), pa,
 409                            asize(desc), desc, asize(bugs), bugs, &oenv))
 410     {
 411         return 0;
 412     }
 413
 414     /* Check input for something sensible */
 415     if ((p_num < 0) || (n_num < 0))
 416     {
 417         gmx_fatal(FARGS, "Negative number of ions to add?");
 418     }
 419
 420     if (conc > 0 && (p_num > 0 || n_num > 0))
 421     {
 422         fprintf(stderr, "WARNING: -conc specified, overriding -nn and -np.\n");
 423     }
 424
 425     /* Read atom positions and charges */
 426     read_tps_conf(ftp2fn(efTPX, NFILE, fnm), title, &top, &ePBC, &x, &v, box, FALSE);
 427     atoms = top.atoms;
 428
 429     /* Compute total charge */
 430     qtot = 0;
 431     for (i = 0; (i < atoms.nr); i++)
 432     {
 433         qtot += atoms.atom[i].q;
 434     }
 435     iqtot = gmx_nint(qtot);
 436
 437
 438     if (conc > 0)
 439     {
 440         /* Compute number of ions to be added */
 441         vol   = det(box);
 442         nsalt = gmx_nint(conc*vol*AVOGADRO/1e24);
 443         p_num = abs(nsalt*n_q);
 444         n_num = abs(nsalt*p_q);
 445     }
 446     if (bNeutral)
 447     {
 448         int qdelta = p_num*p_q + n_num*n_q + iqtot;
 449
 450         /* Check if the system is neutralizable
 451          * is (qdelta == p_q*p_num + n_q*n_num) solvable for p_num and n_num? */
 452         int gcd = gmx_greatest_common_divisor(n_q, p_q);
 453         if ((qdelta % gcd) != 0)
 454         {
 455             gmx_fatal(FARGS, "Can't neutralize this system using -nq %d and"
 456                       " -pq %d.\n", n_q, p_q);
 457         }
 458
 459         while (qdelta != 0)
 460         {
 461             while (qdelta < 0)
 462             {
 463                 p_num++;
 464                 qdelta += p_q;
 465             }
 466             while (qdelta > 0)
 467             {
 468                 n_num++;
 469                 qdelta += n_q;
 470             }
 471         }
 472     }
 473
 474     if ((p_num == 0) && (n_num == 0))
 475     {
 476         fprintf(stderr, "No ions to add.\n");
 477         exit(0);
 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
 514     if (opt2bSet("-p", NFILE, fnm))
 515     {
 516         update_topol(opt2fn("-p", NFILE, fnm), p_num, n_num, p_name, n_name, grpname);
 517     }
 518
 519     snew(bSet, nw);
 520     snew(repl, nw);
 521
 522     snew(v, atoms.nr);
 523     snew(atoms.pdbinfo, atoms.nr);
 524
 525     set_pbc(&pbc, ePBC, box);
 526
 527     if (seed == 0)
 528     {
 529         rng = gmx_rng_init(gmx_rng_make_seed());
 530     }
 531     else
 532     {
 533         rng = gmx_rng_init(seed);
 534     }
 535     /* Now loop over the ions that have to be placed */
 536     while (p_num-- > 0)
 537     {
 538         insert_ion(nsa, &nw, bSet, repl, index, x, &pbc,
 539                    1, p_q, p_name, &atoms, rmin, rng);
 540     }
 541     while (n_num-- > 0)
 542     {
 543         insert_ion(nsa, &nw, bSet, repl, index, x, &pbc,
 544                    -1, n_q, n_name, &atoms, rmin, rng);
 545     }
 546     gmx_rng_destroy(rng);
 547     fprintf(stderr, "\n");
 548
 549     if (nw)
 550     {
 551         sort_ions(nsa, nw, repl, index, &atoms, x, p_name, n_name);
 552     }
 553
 554     sfree(atoms.pdbinfo);
 555     atoms.pdbinfo = NULL;
 556     write_sto_conf(ftp2fn(efSTO, NFILE, fnm), *top.name, &atoms, x, NULL, ePBC,
 557                    box);
 558
 559     return 0;
 560 }