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,2015, 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     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     char     temporary_filename[STRLEN];
 246
 247     printf("\nProcessing topology\n");
 248     fpin  = gmx_ffopen(topinout, "r");
 249     strncpy(temporary_filename, "temp.topXXXXXX", STRLEN);
 250     gmx_tmpnam(temporary_filename);
 251     fpout = gmx_ffopen(temporary_filename, "w");
 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         strcpy(buf2, buf);
 262         if ((temp = strchr(buf2, '\n')) != NULL)
 263         {
 264             temp[0] = '\0';
 265         }
 266         ltrim(buf2);
 267         if (buf2[0] == '[')
 268         {
 269             buf2[0] = ' ';
 270             if ((temp = strchr(buf2, '\n')) != NULL)
 271             {
 272                 temp[0] = '\0';
 273             }
 274             rtrim(buf2);
 275             if (buf2[strlen(buf2)-1] == ']')
 276             {
 277                 buf2[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] = 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     /* use gmx_ffopen to generate backup of topinout */
 345     fpout = gmx_ffopen(topinout, "w");
 346     gmx_ffclose(fpout);
 347     rename(temporary_filename, topinout);
 348 }
 349
 350 int gmx_genion(int argc, char *argv[])
 351 {
 352     const char        *desc[] = {
 353         "[THISMODULE] randomly replaces solvent molecules with monoatomic ions.",
 354         "The group of solvent molecules should be continuous and all molecules",
 355         "should have the same number of atoms.",
 356         "The user should add the ion molecules to the topology file or use",
 357         "the [TT]-p[tt] option to automatically modify the topology.[PAR]",
 358         "The ion molecule type, residue and atom names in all force fields",
 359         "are the capitalized element names without sign. This molecule name",
 360         "should be given with [TT]-pname[tt] or [TT]-nname[tt], and the",
 361         "[TT][molecules][tt] section of your topology updated accordingly,",
 362         "either by hand or with [TT]-p[tt]. Do not use an atom name instead!",
 363         "[PAR]Ions which can have multiple charge states get the multiplicity",
 364         "added, without sign, for the uncommon states only.[PAR]",
 365         "For larger ions, e.g. sulfate we recommended using [gmx-insert-molecules]."
 366     };
 367     const char        *bugs[] = {
 368         "If you specify a salt concentration existing ions are not taken into "
 369         "account. In effect you therefore specify the amount of salt to be added.",
 370     };
 371     static int         p_num    = 0, n_num = 0, p_q = 1, n_q = -1;
 372     static const char *p_name   = "NA", *n_name = "CL";
 373     static real        rmin     = 0.6, conc = 0;
 374     static int         seed     = 1993;
 375     static gmx_bool    bNeutral = FALSE;
 376     static t_pargs     pa[]     = {
 377         { "-np",    FALSE, etINT,  {&p_num}, "Number of positive ions"       },
 378         { "-pname", FALSE, etSTR,  {&p_name}, "Name of the positive ion"      },
 379         { "-pq",    FALSE, etINT,  {&p_q},   "Charge of the positive ion"    },
 380         { "-nn",    FALSE, etINT,  {&n_num}, "Number of negative ions"       },
 381         { "-nname", FALSE, etSTR,  {&n_name}, "Name of the negative ion"      },
 382         { "-nq",    FALSE, etINT,  {&n_q},   "Charge of the negative ion"    },
 383         { "-rmin",  FALSE, etREAL, {&rmin},  "Minimum distance between ions" },
 384         { "-seed",  FALSE, etINT,  {&seed},  "Seed for random number generator" },
 385         { "-conc",  FALSE, etREAL, {&conc},
 386           "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." },
 387         { "-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]. "}
 388     };
 389     t_topology         top;
 390     rvec              *x, *v;
 391     real               vol, qtot;
 392     matrix             box;
 393     t_atoms            atoms;
 394     t_pbc              pbc;
 395     int               *repl, ePBC;
 396     atom_id           *index;
 397     char              *grpname, title[STRLEN];
 398     gmx_bool          *bSet;
 399     int                i, nw, nwa, nsa, nsalt, iqtot;
 400     output_env_t       oenv;
 401     gmx_rng_t          rng;
 402     t_filenm           fnm[] = {
 403         { efTPX, NULL,  NULL,      ffREAD  },
 404         { efNDX, NULL,  NULL,      ffOPTRD },
 405         { efSTO, "-o",  NULL,      ffWRITE },
 406         { efTOP, "-p",  "topol",   ffOPTRW }
 407     };
 408 #define NFILE asize(fnm)
 409
 410     if (!parse_common_args(&argc, argv, PCA_BE_NICE, NFILE, fnm, asize(pa), pa,
 411                            asize(desc), desc, asize(bugs), bugs, &oenv))
 412     {
 413         return 0;
 414     }
 415
 416     /* Check input for something sensible */
 417     if ((p_num < 0) || (n_num < 0))
 418     {
 419         gmx_fatal(FARGS, "Negative number of ions to add?");
 420     }
 421
 422     if (conc > 0 && (p_num > 0 || n_num > 0))
 423     {
 424         fprintf(stderr, "WARNING: -conc specified, overriding -nn and -np.\n");
 425     }
 426
 427     /* Read atom positions and charges */
 428     read_tps_conf(ftp2fn(efTPX, NFILE, fnm), title, &top, &ePBC, &x, &v, box, FALSE);
 429     atoms = top.atoms;
 430
 431     /* Compute total charge */
 432     qtot = 0;
 433     for (i = 0; (i < atoms.nr); i++)
 434     {
 435         qtot += atoms.atom[i].q;
 436     }
 437     iqtot = gmx_nint(qtot);
 438
 439
 440     if (conc > 0)
 441     {
 442         /* Compute number of ions to be added */
 443         vol   = det(box);
 444         nsalt = gmx_nint(conc*vol*AVOGADRO/1e24);
 445         p_num = abs(nsalt*n_q);
 446         n_num = abs(nsalt*p_q);
 447     }
 448     if (bNeutral)
 449     {
 450         int qdelta = p_num*p_q + n_num*n_q + iqtot;
 451
 452         /* Check if the system is neutralizable
 453          * is (qdelta == p_q*p_num + n_q*n_num) solvable for p_num and n_num? */
 454         int gcd = gmx_greatest_common_divisor(n_q, p_q);
 455         if ((qdelta % gcd) != 0)
 456         {
 457             gmx_fatal(FARGS, "Can't neutralize this system using -nq %d and"
 458                       " -pq %d.\n", n_q, p_q);
 459         }
 460
 461         while (qdelta != 0)
 462         {
 463             while (qdelta < 0)
 464             {
 465                 p_num++;
 466                 qdelta += p_q;
 467             }
 468             while (qdelta > 0)
 469             {
 470                 n_num++;
 471                 qdelta += n_q;
 472             }
 473         }
 474     }
 475
 476     if ((p_num == 0) && (n_num == 0))
 477     {
 478         fprintf(stderr, "No ions to add, will just copy input configuration.\n");
 479     }
 480     else
 481     {
 482         printf("Will try to add %d %s ions and %d %s ions.\n",
 483                p_num, p_name, n_num, n_name);
 484         printf("Select a continuous group of solvent molecules\n");
 485         get_index(&atoms, ftp2fn_null(efNDX, NFILE, fnm), 1, &nwa, &index, &grpname);
 486         for (i = 1; i < nwa; i++)
 487         {
 488             if (index[i] != index[i-1]+1)
 489             {
 490                 gmx_fatal(FARGS, "The solvent group %s is not continuous: "
 491                           "index[%d]=%d, index[%d]=%d",
 492                           grpname, i, index[i-1]+1, i+1, index[i]+1);
 493             }
 494         }
 495         nsa = 1;
 496         while ((nsa < nwa) &&
 497                (atoms.atom[index[nsa]].resind ==
 498                 atoms.atom[index[nsa-1]].resind))
 499         {
 500             nsa++;
 501         }
 502         if (nwa % nsa)
 503         {
 504             gmx_fatal(FARGS, "Your solvent group size (%d) is not a multiple of %d",
 505                       nwa, nsa);
 506         }
 507         nw = nwa/nsa;
 508         fprintf(stderr, "Number of (%d-atomic) solvent molecules: %d\n", nsa, nw);
 509         if (p_num+n_num > nw)
 510         {
 511             gmx_fatal(FARGS, "Not enough solvent for adding ions");
 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     }
 557     write_sto_conf(ftp2fn(efSTO, NFILE, fnm), *top.name, &atoms, x, NULL, ePBC, box);
 558
 559     return 0;
 560 }