src/gromacs/gmxana/gmx_genion.c

   1 /*
   2  *
   3  *                This source code is part of
   4  *
   5  *                 G   R   O   M   A   C   S
   6  *
   7  *          GROningen MAchine for Chemical Simulations
   8  *
   9  *                        VERSION 3.2.0
  10  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  11  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  12  * Copyright (c) 2001-2004, The GROMACS development team,
  13  * check out http://www.gromacs.org for more information.
  14
  15  * This program is free software; you can redistribute it and/or
  16  * modify it under the terms of the GNU General Public License
  17  * as published by the Free Software Foundation; either version 2
  18  * of the License, or (at your option) any later version.
  19  *
  20  * If you want to redistribute modifications, please consider that
  21  * scientific software is very special. Version control is crucial -
  22  * bugs must be traceable. We will be happy to consider code for
  23  * inclusion in the official distribution, but derived work must not
  24  * be called official GROMACS. Details are found in the README & COPYING
  25  * files - if they are missing, get the official version at www.gromacs.org.
  26  *
  27  * To help us fund GROMACS development, we humbly ask that you cite
  28  * the papers on the package - you can find them in the top README file.
  29  *
  30  * For more info, check our website at http://www.gromacs.org
  31  *
  32  * And Hey:
  33  * Green Red Orange Magenta Azure Cyan Skyblue
  34  */
  35 #ifdef HAVE_CONFIG_H
  36 #include <config.h>
  37 #endif
  38
  39 #include <ctype.h>
  40 #include "string2.h"
  41 #include "smalloc.h"
  42 #include "sysstuff.h"
  43 #include "confio.h"
  44 #include "statutil.h"
  45 #include "pbc.h"
  46 #include "force.h"
  47 #include "gmx_fatal.h"
  48 #include "futil.h"
  49 #include "maths.h"
  50 #include "macros.h"
  51 #include "vec.h"
  52 #include "tpxio.h"
  53 #include "mdrun.h"
  54 #include "main.h"
  55 #include "random.h"
  56 #include "index.h"
  57 #include "mtop_util.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, int *seed)
  66 {
  67     int             i, ei,nw;
  68     real            rmin2;
  69     rvec            dx;
  70     gmx_large_int_t maxrand;
  71
  72     ei       = -1;
  73     nw       = *nwater;
  74     maxrand  = nw;
  75     maxrand *= 1000;
  76
  77     do
  78     {
  79         ei = nw*rando(seed);
  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 = 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] = 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] = 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  = ffopen(topinout, "r");
 245     fpout = 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] = strdup(buf);
 294             nmol_line++;
 295         }
 296     }
 297     ffclose(fpin);
 298
 299     if (sol_line == -1)
 300     {
 301         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         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     ffclose(fpout);
 338     /* use ffopen to generate backup of topinout */
 339     fpout = ffopen(topinout, "w");
 340     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         "[TT]genion[tt] 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 [TT]genbox[tt]."
 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     t_filenm           fnm[] = {
 397         { efTPX, NULL,  NULL,      ffREAD  },
 398         { efNDX, NULL,  NULL,      ffOPTRD },
 399         { efSTO, "-o",  NULL,      ffWRITE },
 400         { efTOP, "-p",  "topol",   ffOPTRW }
 401     };
 402 #define NFILE asize(fnm)
 403
 404     parse_common_args(&argc, argv, PCA_BE_NICE, NFILE, fnm, asize(pa), pa,
 405                       asize(desc), desc, asize(bugs), bugs, &oenv);
 406
 407     /* Check input for something sensible */
 408     if ((p_num < 0) || (n_num < 0))
 409     {
 410         gmx_fatal(FARGS, "Negative number of ions to add?");
 411     }
 412
 413     if (conc > 0 && (p_num > 0 || n_num > 0))
 414     {
 415         fprintf(stderr, "WARNING: -conc specified, overriding -nn and -np.\n");
 416     }
 417
 418     /* Read atom positions and charges */
 419     read_tps_conf(ftp2fn(efTPX, NFILE, fnm), title, &top, &ePBC, &x, &v, box, FALSE);
 420     atoms = top.atoms;
 421
 422     /* Compute total charge */
 423     qtot = 0;
 424     for (i = 0; (i < atoms.nr); i++)
 425     {
 426         qtot += atoms.atom[i].q;
 427     }
 428     iqtot = gmx_nint(qtot);
 429
 430
 431     if (conc > 0)
 432     {
 433         /* Compute number of ions to be added */
 434         vol = det(box);
 435         nsalt = gmx_nint(conc*vol*AVOGADRO/1e24);
 436         p_num = abs(nsalt*n_q);
 437         n_num = abs(nsalt*p_q);
 438     }
 439     if (bNeutral)
 440     {
 441         int qdelta = p_num*p_q + n_num*n_q + iqtot;
 442
 443         /* Check if the system is neutralizable
 444          * is (qdelta == p_q*p_num + n_q*n_num) solvable for p_num and n_num? */
 445         int gcd = gmx_greatest_common_divisor(n_q, p_q);
 446         if ((qdelta % gcd) != 0)
 447         {
 448             gmx_fatal(FARGS, "Can't neutralize this system using -nq %d and"
 449                     " -pq %d.\n", n_q, p_q);
 450         }
 451
 452         while (qdelta != 0)
 453         {
 454             while (qdelta < 0)
 455             {
 456                 p_num++;
 457                 qdelta += p_q;
 458             }
 459             while (qdelta > 0)
 460             {
 461                 n_num++;
 462                 qdelta += n_q;
 463             }
 464         }
 465     }
 466
 467     if ((p_num == 0) && (n_num == 0))
 468     {
 469         fprintf(stderr, "No ions to add.\n");
 470         exit(0);
 471     }
 472     else
 473     {
 474         printf("Will try to add %d %s ions and %d %s ions.\n",
 475                p_num, p_name, n_num, n_name);
 476         printf("Select a continuous group of solvent molecules\n");
 477         get_index(&atoms, ftp2fn_null(efNDX, NFILE, fnm), 1, &nwa, &index, &grpname);
 478         for (i = 1; i < nwa; i++)
 479         {
 480             if (index[i] != index[i-1]+1)
 481             {
 482                 gmx_fatal(FARGS, "The solvent group %s is not continuous: "
 483                           "index[%d]=%d, index[%d]=%d",
 484                           grpname, i, index[i-1]+1, i+1, index[i]+1);
 485             }
 486         }
 487         nsa = 1;
 488         while ((nsa < nwa) &&
 489                (atoms.atom[index[nsa]].resind ==
 490                 atoms.atom[index[nsa-1]].resind))
 491         {
 492             nsa++;
 493         }
 494         if (nwa % nsa)
 495         {
 496             gmx_fatal(FARGS, "Your solvent group size (%d) is not a multiple of %d",
 497                       nwa, nsa);
 498         }
 499         nw = nwa/nsa;
 500         fprintf(stderr, "Number of (%d-atomic) solvent molecules: %d\n", nsa, nw);
 501         if (p_num+n_num > nw)
 502         {
 503             gmx_fatal(FARGS, "Not enough solvent for adding ions");
 504         }
 505     }
 506
 507     if (opt2bSet("-p", NFILE, fnm))
 508     {
 509         update_topol(opt2fn("-p", NFILE, fnm), p_num, n_num, p_name, n_name, grpname);
 510     }
 511
 512     snew(bSet, nw);
 513     snew(repl, nw);
 514
 515     snew(v, atoms.nr);
 516     snew(atoms.pdbinfo, atoms.nr);
 517
 518     set_pbc(&pbc, ePBC, box);
 519
 520     /* Now loop over the ions that have to be placed */
 521     while (p_num-- > 0)
 522     {
 523         insert_ion(nsa, &nw, bSet, repl, index, x, &pbc,
 524                    1, p_q, p_name, &atoms, rmin, &seed);
 525     }
 526     while (n_num-- > 0)
 527     {
 528         insert_ion(nsa, &nw, bSet, repl, index, x, &pbc,
 529                    -1, n_q, n_name, &atoms, rmin, &seed);
 530     }
 531     fprintf(stderr, "\n");
 532
 533     if (nw)
 534     {
 535         sort_ions(nsa, nw, repl, index, &atoms, x, p_name, n_name);
 536     }
 537
 538     sfree(atoms.pdbinfo);
 539     atoms.pdbinfo = NULL;
 540     write_sto_conf(ftp2fn(efSTO, NFILE, fnm), *top.name, &atoms, x, NULL, ePBC,
 541                    box);
 542
 543     return 0;
 544 }