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