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