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 #include "gmxpre.h"
  38
  39 #include <ctype.h>
  40 #include <stdlib.h>
  41 #include <string.h>
  42
  43 #include "gromacs/commandline/pargs.h"
  44 #include "gromacs/fileio/confio.h"
  45 #include "gromacs/fileio/tpxio.h"
  46 #include "gromacs/gmxana/gmx_ana.h"
  47 #include "gromacs/legacyheaders/force.h"
  48 #include "gromacs/legacyheaders/macros.h"
  49 #include "gromacs/legacyheaders/mdrun.h"
  50 #include "gromacs/math/utilities.h"
  51 #include "gromacs/math/vec.h"
  52 #include "gromacs/pbcutil/pbc.h"
  53 #include "gromacs/random/random.h"
  54 #include "gromacs/topology/index.h"
  55 #include "gromacs/utility/cstringutil.h"
  56 #include "gromacs/utility/fatalerror.h"
  57 #include "gromacs/utility/futil.h"
  58 #include "gromacs/utility/smalloc.h"
  59
  60 static void insert_ion(int nsa, int *nwater,
  61                        gmx_bool bSet[], int repl[], atom_id index[],
  62                        rvec x[], t_pbc *pbc,
  63                        int sign, int q, const char *ionname,
  64                        t_atoms *atoms,
  65                        real rmin, gmx_rng_t rng)
  66 {
  67     int             i, ei, nw;
  68     real            rmin2;
  69     rvec            dx;
  70     gmx_int64_t     maxrand;
  71
  72     ei       = -1;
  73     nw       = *nwater;
  74     maxrand  = nw;
  75     maxrand *= 1000;
  76
  77     do
  78     {
  79         ei = nw*gmx_rng_uniform_real(rng);
  80         maxrand--;
  81     }
  82     while (bSet[ei] && (maxrand > 0));
  83     if (bSet[ei])
  84     {
  85         gmx_fatal(FARGS, "No more replaceable solvent!");
  86     }
  87
  88     fprintf(stderr, "Replacing solvent molecule %d (atom %d) with %s\n",
  89             ei, index[nsa*ei], ionname);
  90
  91     /* Replace solvent molecule charges with ion charge */
  92     bSet[ei] = TRUE;
  93     repl[ei] = sign;
  94
  95     atoms->atom[index[nsa*ei]].q = q;
  96     for (i = 1; i < nsa; i++)
  97     {
  98         atoms->atom[index[nsa*ei+i]].q = 0;
  99     }
 100
 101     /* Mark all solvent molecules within rmin as unavailable for substitution */
 102     if (rmin > 0)
 103     {
 104         rmin2 = rmin*rmin;
 105         for (i = 0; (i < nw); i++)
 106         {
 107             if (!bSet[i])
 108             {
 109                 pbc_dx(pbc, x[index[nsa*ei]], x[index[nsa*i]], dx);
 110                 if (iprod(dx, dx) < rmin2)
 111                 {
 112                     bSet[i] = TRUE;
 113                 }
 114             }
 115         }
 116     }
 117 }
 118
 119
 120 static char *aname(const char *mname)
 121 {
 122     char *str;
 123     int   i;
 124
 125     str = gmx_strdup(mname);
 126     i   = strlen(str)-1;
 127     while (i > 1 && (isdigit(str[i]) || (str[i] == '+') || (str[i] == '-')))
 128     {
 129         str[i] = '\0';
 130         i--;
 131     }
 132
 133     return str;
 134 }
 135
 136 void sort_ions(int nsa, int nw, int repl[], atom_id index[],
 137                t_atoms *atoms, rvec x[],
 138                const char *p_name, const char *n_name)
 139 {
 140     int    i, j, k, r, np, nn, starta, startr, npi, nni;
 141     rvec  *xt;
 142     char **pptr = NULL, **nptr = NULL, **paptr = NULL, **naptr = NULL;
 143
 144     snew(xt, atoms->nr);
 145
 146     /* Put all the solvent in front and count the added ions */
 147     np = 0;
 148     nn = 0;
 149     j  = index[0];
 150     for (i = 0; i < nw; i++)
 151     {
 152         r = repl[i];
 153         if (r == 0)
 154         {
 155             for (k = 0; k < nsa; k++)
 156             {
 157                 copy_rvec(x[index[nsa*i+k]], xt[j++]);
 158             }
 159         }
 160         else if (r > 0)
 161         {
 162             np++;
 163         }
 164         else if (r < 0)
 165         {
 166             nn++;
 167         }
 168     }
 169
 170     if (np+nn > 0)
 171     {
 172         /* Put the positive and negative ions at the end */
 173         starta = index[nsa*(nw - np - nn)];
 174         startr = atoms->atom[starta].resind;
 175
 176         if (np)
 177         {
 178             snew(pptr, 1);
 179             pptr[0] = gmx_strdup(p_name);
 180             snew(paptr, 1);
 181             paptr[0] = aname(p_name);
 182         }
 183         if (nn)
 184         {
 185             snew(nptr, 1);
 186             nptr[0] = gmx_strdup(n_name);
 187             snew(naptr, 1);
 188             naptr[0] = aname(n_name);
 189         }
 190         npi = 0;
 191         nni = 0;
 192         for (i = 0; i < nw; i++)
 193         {
 194             r = repl[i];
 195             if (r > 0)
 196             {
 197                 j = starta+npi;
 198                 k = startr+npi;
 199                 copy_rvec(x[index[nsa*i]], xt[j]);
 200                 atoms->atomname[j]     = paptr;
 201                 atoms->atom[j].resind  = k;
 202                 atoms->resinfo[k].name = pptr;
 203                 npi++;
 204             }
 205             else if (r < 0)
 206             {
 207                 j = starta+np+nni;
 208                 k = startr+np+nni;
 209                 copy_rvec(x[index[nsa*i]], xt[j]);
 210                 atoms->atomname[j]     = naptr;
 211                 atoms->atom[j].resind  = k;
 212                 atoms->resinfo[k].name = nptr;
 213                 nni++;
 214             }
 215         }
 216         for (i = index[nsa*nw-1]+1; i < atoms->nr; i++)
 217         {
 218             j                  = i-(nsa-1)*(np+nn);
 219             atoms->atomname[j] = atoms->atomname[i];
 220             atoms->atom[j]     = atoms->atom[i];
 221             copy_rvec(x[i], xt[j]);
 222         }
 223         atoms->nr -= (nsa-1)*(np+nn);
 224
 225         /* Copy the new positions back */
 226         for (i = index[0]; i < atoms->nr; i++)
 227         {
 228             copy_rvec(xt[i], x[i]);
 229         }
 230         sfree(xt);
 231     }
 232 }
 233
 234 static void update_topol(const char *topinout, int p_num, int n_num,
 235                          const char *p_name, const char *n_name, char *grpname)
 236 {
 237     FILE    *fpin, *fpout;
 238     char     buf[STRLEN], buf2[STRLEN], *temp, **mol_line = NULL;
 239     int      line, i, nsol, nmol_line, sol_line, nsol_last;
 240     gmx_bool bMolecules;
 241     char     temporary_filename[STRLEN];
 242
 243     printf("\nProcessing topology\n");
 244     fpin  = gmx_ffopen(topinout, "r");
 245     strncpy(temporary_filename, "temp.topXXXXXX", STRLEN);
 246     gmx_tmpnam(temporary_filename);
 247     fpout = gmx_ffopen(temporary_filename, "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(temporary_filename, topinout);
 344 }
 345
 346 int gmx_genion(int argc, char *argv[])
 347 {
 348     const char        *desc[] = {
 349         "[THISMODULE] randomly replaces solvent molecules with monoatomic ions.",
 350         "The group of solvent molecules should be continuous and all molecules",
 351         "should have the same number of atoms.",
 352         "The user should add the ion molecules to the topology file or use",
 353         "the [TT]-p[tt] option to automatically modify the topology.[PAR]",
 354         "The ion molecule type, residue and atom names in all force fields",
 355         "are the capitalized element names without sign. This molecule name",
 356         "should be given with [TT]-pname[tt] or [TT]-nname[tt], and the",
 357         "[TT][molecules][tt] section of your topology updated accordingly,",
 358         "either by hand or with [TT]-p[tt]. Do not use an atom name instead!",
 359         "[PAR]Ions which can have multiple charge states get the multiplicity",
 360         "added, without sign, for the uncommon states only.[PAR]",
 361         "For larger ions, e.g. sulfate we recommended using [gmx-insert-molecules]."
 362     };
 363     const char        *bugs[] = {
 364         "If you specify a salt concentration existing ions are not taken into "
 365         "account. In effect you therefore specify the amount of salt to be added.",
 366     };
 367     static int         p_num    = 0, n_num = 0, p_q = 1, n_q = -1;
 368     static const char *p_name   = "NA", *n_name = "CL";
 369     static real        rmin     = 0.6, conc = 0;
 370     static int         seed     = 1993;
 371     static gmx_bool    bNeutral = FALSE;
 372     static t_pargs     pa[]     = {
 373         { "-np",    FALSE, etINT,  {&p_num}, "Number of positive ions"       },
 374         { "-pname", FALSE, etSTR,  {&p_name}, "Name of the positive ion"      },
 375         { "-pq",    FALSE, etINT,  {&p_q},   "Charge of the positive ion"    },
 376         { "-nn",    FALSE, etINT,  {&n_num}, "Number of negative ions"       },
 377         { "-nname", FALSE, etSTR,  {&n_name}, "Name of the negative ion"      },
 378         { "-nq",    FALSE, etINT,  {&n_q},   "Charge of the negative ion"    },
 379         { "-rmin",  FALSE, etREAL, {&rmin},  "Minimum distance between ions" },
 380         { "-seed",  FALSE, etINT,  {&seed},  "Seed for random number generator" },
 381         { "-conc",  FALSE, etREAL, {&conc},
 382           "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." },
 383         { "-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]. "}
 384     };
 385     t_topology         top;
 386     rvec              *x, *v;
 387     real               vol, qtot;
 388     matrix             box;
 389     t_atoms            atoms;
 390     t_pbc              pbc;
 391     int               *repl, ePBC;
 392     atom_id           *index;
 393     char              *grpname, title[STRLEN];
 394     gmx_bool          *bSet;
 395     int                i, nw, nwa, nsa, nsalt, iqtot;
 396     output_env_t       oenv;
 397     gmx_rng_t          rng;
 398     t_filenm           fnm[] = {
 399         { efTPR, NULL,  NULL,      ffREAD  },
 400         { efNDX, NULL,  NULL,      ffOPTRD },
 401         { efSTO, "-o",  NULL,      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), title, &top, &ePBC, &x, &v, box, FALSE);
 425     atoms = top.atoms;
 426
 427     /* Compute total charge */
 428     qtot = 0;
 429     for (i = 0; (i < atoms.nr); i++)
 430     {
 431         qtot += atoms.atom[i].q;
 432     }
 433     iqtot = gmx_nint(qtot);
 434
 435
 436     if (conc > 0)
 437     {
 438         /* Compute number of ions to be added */
 439         vol   = det(box);
 440         nsalt = gmx_nint(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         if (seed == 0)
 524         {
 525             rng = gmx_rng_init(gmx_rng_make_seed());
 526         }
 527         else
 528         {
 529             rng = gmx_rng_init(seed);
 530         }
 531         /* Now loop over the ions that have to be placed */
 532         while (p_num-- > 0)
 533         {
 534             insert_ion(nsa, &nw, bSet, repl, index, x, &pbc,
 535                        1, p_q, p_name, &atoms, rmin, rng);
 536         }
 537         while (n_num-- > 0)
 538         {
 539             insert_ion(nsa, &nw, bSet, repl, index, x, &pbc,
 540                        -1, n_q, n_name, &atoms, rmin, rng);
 541         }
 542         gmx_rng_destroy(rng);
 543         fprintf(stderr, "\n");
 544
 545         if (nw)
 546         {
 547             sort_ions(nsa, nw, repl, index, &atoms, x, p_name, n_name);
 548         }
 549
 550         sfree(atoms.pdbinfo);
 551         atoms.pdbinfo = NULL;
 552     }
 553     write_sto_conf(ftp2fn(efSTO, NFILE, fnm), *top.name, &atoms, x, NULL, ePBC, box);
 554
 555     return 0;
 556 }