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     if (!parse_common_args(&argc, argv, PCA_BE_NICE, NFILE, fnm, asize(pa), pa,
 405                            asize(desc), desc, asize(bugs), bugs, &oenv))
 406     {
 407         return 0;
 408     }
 409
 410     /* Check input for something sensible */
 411     if ((p_num < 0) || (n_num < 0))
 412     {
 413         gmx_fatal(FARGS, "Negative number of ions to add?");
 414     }
 415
 416     if (conc > 0 && (p_num > 0 || n_num > 0))
 417     {
 418         fprintf(stderr, "WARNING: -conc specified, overriding -nn and -np.\n");
 419     }
 420
 421     /* Read atom positions and charges */
 422     read_tps_conf(ftp2fn(efTPX, NFILE, fnm), title, &top, &ePBC, &x, &v, box, FALSE);
 423     atoms = top.atoms;
 424
 425     /* Compute total charge */
 426     qtot = 0;
 427     for (i = 0; (i < atoms.nr); i++)
 428     {
 429         qtot += atoms.atom[i].q;
 430     }
 431     iqtot = gmx_nint(qtot);
 432
 433
 434     if (conc > 0)
 435     {
 436         /* Compute number of ions to be added */
 437         vol   = det(box);
 438         nsalt = gmx_nint(conc*vol*AVOGADRO/1e24);
 439         p_num = abs(nsalt*n_q);
 440         n_num = abs(nsalt*p_q);
 441     }
 442     if (bNeutral)
 443     {
 444         int qdelta = p_num*p_q + n_num*n_q + iqtot;
 445
 446         /* Check if the system is neutralizable
 447          * is (qdelta == p_q*p_num + n_q*n_num) solvable for p_num and n_num? */
 448         int gcd = gmx_greatest_common_divisor(n_q, p_q);
 449         if ((qdelta % gcd) != 0)
 450         {
 451             gmx_fatal(FARGS, "Can't neutralize this system using -nq %d and"
 452                       " -pq %d.\n", n_q, p_q);
 453         }
 454
 455         while (qdelta != 0)
 456         {
 457             while (qdelta < 0)
 458             {
 459                 p_num++;
 460                 qdelta += p_q;
 461             }
 462             while (qdelta > 0)
 463             {
 464                 n_num++;
 465                 qdelta += n_q;
 466             }
 467         }
 468     }
 469
 470     if ((p_num == 0) && (n_num == 0))
 471     {
 472         fprintf(stderr, "No ions to add.\n");
 473         exit(0);
 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
 510     if (opt2bSet("-p", NFILE, fnm))
 511     {
 512         update_topol(opt2fn("-p", NFILE, fnm), p_num, n_num, p_name, n_name, grpname);
 513     }
 514
 515     snew(bSet, nw);
 516     snew(repl, nw);
 517
 518     snew(v, atoms.nr);
 519     snew(atoms.pdbinfo, atoms.nr);
 520
 521     set_pbc(&pbc, ePBC, box);
 522
 523     /* Now loop over the ions that have to be placed */
 524     while (p_num-- > 0)
 525     {
 526         insert_ion(nsa, &nw, bSet, repl, index, x, &pbc,
 527                    1, p_q, p_name, &atoms, rmin, &seed);
 528     }
 529     while (n_num-- > 0)
 530     {
 531         insert_ion(nsa, &nw, bSet, repl, index, x, &pbc,
 532                    -1, n_q, n_name, &atoms, rmin, &seed);
 533     }
 534     fprintf(stderr, "\n");
 535
 536     if (nw)
 537     {
 538         sort_ions(nsa, nw, repl, index, &atoms, x, p_name, n_name);
 539     }
 540
 541     sfree(atoms.pdbinfo);
 542     atoms.pdbinfo = NULL;
 543     write_sto_conf(ftp2fn(efSTO, NFILE, fnm), *top.name, &atoms, x, NULL, ePBC,
 544                    box);
 545
 546     return 0;
 547 }