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 <ctype.h>
  40 #include <stdlib.h>
  41 #include <string.h>
  42
  43 #include "gromacs/utility/cstringutil.h"
  44 #include "gromacs/utility/smalloc.h"
  45 #include "gromacs/fileio/confio.h"
  46 #include "gromacs/commandline/pargs.h"
  47 #include "gromacs/pbcutil/pbc.h"
  48 #include "gromacs/legacyheaders/force.h"
  49 #include "gromacs/utility/fatalerror.h"
  50 #include "gromacs/utility/futil.h"
  51 #include "gromacs/math/utilities.h"
  52 #include "gromacs/legacyheaders/macros.h"
  53 #include "gromacs/math/vec.h"
  54 #include "gromacs/fileio/tpxio.h"
  55 #include "gromacs/legacyheaders/mdrun.h"
  56 #include "gromacs/random/random.h"
  57 #include "gromacs/topology/index.h"
  58 #include "gmx_ana.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 #define TEMP_FILENM "temp.top"
 238     FILE    *fpin, *fpout;
 239     char     buf[STRLEN], buf2[STRLEN], *temp, **mol_line = NULL;
 240     int      line, i, nsol, nmol_line, sol_line, nsol_last;
 241     gmx_bool bMolecules;
 242
 243     printf("\nProcessing topology\n");
 244     fpin  = gmx_ffopen(topinout, "r");
 245     fpout = gmx_ffopen(TEMP_FILENM, "w");
 246
 247     line       = 0;
 248     bMolecules = FALSE;
 249     nmol_line  = 0;
 250     sol_line   = -1;
 251     nsol_last  = -1;
 252     while (fgets(buf, STRLEN, fpin))
 253     {
 254         line++;
 255         strcpy(buf2, buf);
 256         if ((temp = strchr(buf2, '\n')) != NULL)
 257         {
 258             temp[0] = '\0';
 259         }
 260         ltrim(buf2);
 261         if (buf2[0] == '[')
 262         {
 263             buf2[0] = ' ';
 264             if ((temp = strchr(buf2, '\n')) != NULL)
 265             {
 266                 temp[0] = '\0';
 267             }
 268             rtrim(buf2);
 269             if (buf2[strlen(buf2)-1] == ']')
 270             {
 271                 buf2[strlen(buf2)-1] = '\0';
 272                 ltrim(buf2);
 273                 rtrim(buf2);
 274                 bMolecules = (gmx_strcasecmp(buf2, "molecules") == 0);
 275             }
 276             fprintf(fpout, "%s", buf);
 277         }
 278         else if (!bMolecules)
 279         {
 280             fprintf(fpout, "%s", buf);
 281         }
 282         else
 283         {
 284             /* Check if this is a line with solvent molecules */
 285             sscanf(buf, "%s", buf2);
 286             if (gmx_strcasecmp(buf2, grpname) == 0)
 287             {
 288                 sol_line = nmol_line;
 289                 sscanf(buf, "%*s %d", &nsol_last);
 290             }
 291             /* Store this molecules section line */
 292             srenew(mol_line, nmol_line+1);
 293             mol_line[nmol_line] = gmx_strdup(buf);
 294             nmol_line++;
 295         }
 296     }
 297     gmx_ffclose(fpin);
 298
 299     if (sol_line == -1)
 300     {
 301         gmx_ffclose(fpout);
 302         gmx_fatal(FARGS, "No line with moleculetype '%s' found the [ molecules ] section of file '%s'", grpname, topinout);
 303     }
 304     if (nsol_last < p_num+n_num)
 305     {
 306         gmx_ffclose(fpout);
 307         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);
 308     }
 309
 310     /* Print all the molecule entries */
 311     for (i = 0; i < nmol_line; i++)
 312     {
 313         if (i != sol_line)
 314         {
 315             fprintf(fpout, "%s", mol_line[i]);
 316         }
 317         else
 318         {
 319             printf("Replacing %d solute molecules in topology file (%s) "
 320                    " by %d %s and %d %s ions.\n",
 321                    p_num+n_num, topinout, p_num, p_name, n_num, n_name);
 322             nsol_last -= p_num + n_num;
 323             if (nsol_last > 0)
 324             {
 325                 fprintf(fpout, "%-10s  %d\n", grpname, nsol_last);
 326             }
 327             if (p_num > 0)
 328             {
 329                 fprintf(fpout, "%-15s  %d\n", p_name, p_num);
 330             }
 331             if (n_num > 0)
 332             {
 333                 fprintf(fpout, "%-15s  %d\n", n_name, n_num);
 334             }
 335         }
 336     }
 337     gmx_ffclose(fpout);
 338     /* use gmx_ffopen to generate backup of topinout */
 339     fpout = gmx_ffopen(topinout, "w");
 340     gmx_ffclose(fpout);
 341     rename(TEMP_FILENM, topinout);
 342 #undef TEMP_FILENM
 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 [TT].tpr[tt] 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         { efTPX, 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(efTPX, 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 }