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