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 "physics.h"
  53 #include "vec.h"
  54 #include "tpxio.h"
  55 #include "mdrun.h"
  56 #include "calcpot.h"
  57 #include "main.h"
  58 #include "random.h"
  59 #include "index.h"
  60 #include "mtop_util.h"
  61 #include "gmx_ana.h"
  62
  63 static void insert_ion(int nsa, int *nwater,
  64                        gmx_bool bSet[], int repl[], atom_id index[],
  65                        real pot[], rvec x[], t_pbc *pbc,
  66                        int sign, int q, const char *ionname,
  67                        t_mdatoms *mdatoms,
  68                        real rmin, gmx_bool bRandom, int *seed)
  69 {
  70     int             i, ii, ei, owater, wlast, m, nw;
  71     real            extr_e, poti, rmin2;
  72     rvec            xei, dx;
  73     gmx_bool        bSub = FALSE;
  74     gmx_large_int_t maxrand;
  75
  76     ei       = -1;
  77     nw       = *nwater;
  78     maxrand  = nw;
  79     maxrand *= 1000;
  80     if (bRandom)
  81     {
  82         do
  83         {
  84             ei = nw*rando(seed);
  85             maxrand--;
  86         }
  87         while (bSet[ei] && (maxrand > 0));
  88         if (bSet[ei])
  89         {
  90             gmx_fatal(FARGS, "No more replaceable solvent!");
  91         }
  92     }
  93     else
  94     {
  95         extr_e = 0;
  96         for (i = 0; (i < nw); i++)
  97         {
  98             if (!bSet[i])
  99             {
 100                 ii   = index[nsa*i];
 101                 poti = pot[ii];
 102                 if (q > 0)
 103                 {
 104                     if ((poti <= extr_e) || !bSub)
 105                     {
 106                         extr_e = poti;
 107                         ei     = i;
 108                         bSub   = TRUE;
 109                     }
 110                 }
 111                 else
 112                 {
 113                     if ((poti >= extr_e) || !bSub)
 114                     {
 115                         extr_e = poti;
 116                         ei     = i;
 117                         bSub   = TRUE;
 118                     }
 119                 }
 120             }
 121         }
 122         if (ei == -1)
 123         {
 124             gmx_fatal(FARGS, "No more replaceable solvent!");
 125         }
 126     }
 127     fprintf(stderr, "Replacing solvent molecule %d (atom %d) with %s\n",
 128             ei, index[nsa*ei], ionname);
 129
 130     /* Replace solvent molecule charges with ion charge */
 131     bSet[ei] = TRUE;
 132     repl[ei] = sign;
 133     mdatoms->chargeA[index[nsa*ei]] = q;
 134     for (i = 1; i < nsa; i++)
 135     {
 136         mdatoms->chargeA[index[nsa*ei+i]] = 0;
 137     }
 138
 139     /* Mark all solvent molecules within rmin as unavailable for substitution */
 140     if (rmin > 0)
 141     {
 142         rmin2 = rmin*rmin;
 143         for (i = 0; (i < nw); i++)
 144         {
 145             if (!bSet[i])
 146             {
 147                 pbc_dx(pbc, x[index[nsa*ei]], x[index[nsa*i]], dx);
 148                 if (iprod(dx, dx) < rmin2)
 149                 {
 150                     bSet[i] = TRUE;
 151                 }
 152             }
 153         }
 154     }
 155 }
 156
 157 static char *aname(const char *mname)
 158 {
 159     char *str;
 160     int   i;
 161
 162     str = strdup(mname);
 163     i   = strlen(str)-1;
 164     while (i > 1 && (isdigit(str[i]) || (str[i] == '+') || (str[i] == '-')))
 165     {
 166         str[i] = '\0';
 167         i--;
 168     }
 169
 170     return str;
 171 }
 172
 173 void sort_ions(int nsa, int nw, int repl[], atom_id index[],
 174                t_atoms *atoms, rvec x[],
 175                const char *p_name, const char *n_name)
 176 {
 177     int    i, j, k, r, np, nn, starta, startr, npi, nni;
 178     rvec  *xt;
 179     char **pptr = NULL, **nptr = NULL, **paptr = NULL, **naptr = NULL;
 180
 181     snew(xt, atoms->nr);
 182
 183     /* Put all the solvent in front and count the added ions */
 184     np = 0;
 185     nn = 0;
 186     j  = index[0];
 187     for (i = 0; i < nw; i++)
 188     {
 189         r = repl[i];
 190         if (r == 0)
 191         {
 192             for (k = 0; k < nsa; k++)
 193             {
 194                 copy_rvec(x[index[nsa*i+k]], xt[j++]);
 195             }
 196         }
 197         else if (r > 0)
 198         {
 199             np++;
 200         }
 201         else if (r < 0)
 202         {
 203             nn++;
 204         }
 205     }
 206
 207     if (np+nn > 0)
 208     {
 209         /* Put the positive and negative ions at the end */
 210         starta = index[nsa*(nw - np - nn)];
 211         startr = atoms->atom[starta].resind;
 212
 213         if (np)
 214         {
 215             snew(pptr, 1);
 216             pptr[0] = strdup(p_name);
 217             snew(paptr, 1);
 218             paptr[0] = aname(p_name);
 219         }
 220         if (nn)
 221         {
 222             snew(nptr, 1);
 223             nptr[0] = strdup(n_name);
 224             snew(naptr, 1);
 225             naptr[0] = aname(n_name);
 226         }
 227         npi = 0;
 228         nni = 0;
 229         for (i = 0; i < nw; i++)
 230         {
 231             r = repl[i];
 232             if (r > 0)
 233             {
 234                 j = starta+npi;
 235                 k = startr+npi;
 236                 copy_rvec(x[index[nsa*i]], xt[j]);
 237                 atoms->atomname[j]     = paptr;
 238                 atoms->atom[j].resind  = k;
 239                 atoms->resinfo[k].name = pptr;
 240                 npi++;
 241             }
 242             else if (r < 0)
 243             {
 244                 j = starta+np+nni;
 245                 k = startr+np+nni;
 246                 copy_rvec(x[index[nsa*i]], xt[j]);
 247                 atoms->atomname[j]     = naptr;
 248                 atoms->atom[j].resind  = k;
 249                 atoms->resinfo[k].name = nptr;
 250                 nni++;
 251             }
 252         }
 253         for (i = index[nsa*nw-1]+1; i < atoms->nr; i++)
 254         {
 255             j                  = i-(nsa-1)*(np+nn);
 256             atoms->atomname[j] = atoms->atomname[i];
 257             atoms->atom[j]     = atoms->atom[i];
 258             copy_rvec(x[i], xt[j]);
 259         }
 260         atoms->nr -= (nsa-1)*(np+nn);
 261
 262         /* Copy the new positions back */
 263         for (i = index[0]; i < atoms->nr; i++)
 264         {
 265             copy_rvec(xt[i], x[i]);
 266         }
 267         sfree(xt);
 268     }
 269 }
 270
 271 static void update_topol(const char *topinout, int p_num, int n_num,
 272                          const char *p_name, const char *n_name, char *grpname)
 273 {
 274 #define TEMP_FILENM "temp.top"
 275     FILE    *fpin, *fpout;
 276     char     buf[STRLEN], buf2[STRLEN], *temp, **mol_line = NULL;
 277     int      line, i, nsol, nmol_line, sol_line, nsol_last;
 278     gmx_bool bMolecules;
 279
 280     printf("\nProcessing topology\n");
 281     fpin  = ffopen(topinout, "r");
 282     fpout = ffopen(TEMP_FILENM, "w");
 283
 284     line       = 0;
 285     bMolecules = FALSE;
 286     nmol_line  = 0;
 287     sol_line   = -1;
 288     nsol_last  = -1;
 289     while (fgets(buf, STRLEN, fpin))
 290     {
 291         line++;
 292         strcpy(buf2, buf);
 293         if ((temp = strchr(buf2, '\n')) != NULL)
 294         {
 295             temp[0] = '\0';
 296         }
 297         ltrim(buf2);
 298         if (buf2[0] == '[')
 299         {
 300             buf2[0] = ' ';
 301             if ((temp = strchr(buf2, '\n')) != NULL)
 302             {
 303                 temp[0] = '\0';
 304             }
 305             rtrim(buf2);
 306             if (buf2[strlen(buf2)-1] == ']')
 307             {
 308                 buf2[strlen(buf2)-1] = '\0';
 309                 ltrim(buf2);
 310                 rtrim(buf2);
 311                 bMolecules = (gmx_strcasecmp(buf2, "molecules") == 0);
 312             }
 313             fprintf(fpout, "%s", buf);
 314         }
 315         else if (!bMolecules)
 316         {
 317             fprintf(fpout, "%s", buf);
 318         }
 319         else
 320         {
 321             /* Check if this is a line with solvent molecules */
 322             sscanf(buf, "%s", buf2);
 323             if (gmx_strcasecmp(buf2, grpname) == 0)
 324             {
 325                 sol_line = nmol_line;
 326                 sscanf(buf, "%*s %d", &nsol_last);
 327             }
 328             /* Store this molecules section line */
 329             srenew(mol_line, nmol_line+1);
 330             mol_line[nmol_line] = strdup(buf);
 331             nmol_line++;
 332         }
 333     }
 334     ffclose(fpin);
 335
 336     if (sol_line == -1)
 337     {
 338         ffclose(fpout);
 339         gmx_fatal(FARGS, "No line with moleculetype '%s' found the [ molecules ] section of file '%s'", grpname, topinout);
 340     }
 341     if (nsol_last < p_num+n_num)
 342     {
 343         ffclose(fpout);
 344         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);
 345     }
 346
 347     /* Print all the molecule entries */
 348     for (i = 0; i < nmol_line; i++)
 349     {
 350         if (i != sol_line)
 351         {
 352             fprintf(fpout, "%s", mol_line[i]);
 353         }
 354         else
 355         {
 356             printf("Replacing %d solute molecules in topology file (%s) "
 357                    " by %d %s and %d %s ions.\n",
 358                    p_num+n_num, topinout, p_num, p_name, n_num, n_name);
 359             nsol_last -= p_num + n_num;
 360             if (nsol_last > 0)
 361             {
 362                 fprintf(fpout, "%-10s  %d\n", grpname, nsol_last);
 363             }
 364             if (p_num > 0)
 365             {
 366                 fprintf(fpout, "%-15s  %d\n", p_name, p_num);
 367             }
 368             if (n_num > 0)
 369             {
 370                 fprintf(fpout, "%-15s  %d\n", n_name, n_num);
 371             }
 372         }
 373     }
 374     ffclose(fpout);
 375     /* use ffopen to generate backup of topinout */
 376     fpout = ffopen(topinout, "w");
 377     ffclose(fpout);
 378     rename(TEMP_FILENM, topinout);
 379 #undef TEMP_FILENM
 380 }
 381
 382 int gmx_genion(int argc, char *argv[])
 383 {
 384     const char        *desc[] = {
 385         "[TT]genion[tt] replaces solvent molecules by monoatomic ions at",
 386         "the position of the first atoms with the most favorable electrostatic",
 387         "potential or at random. The potential is calculated on all atoms, using",
 388         "normal GROMACS particle-based methods (in contrast to other methods",
 389         "based on solving the Poisson-Boltzmann equation).",
 390         "The potential is recalculated after every ion insertion.",
 391         "If specified in the run input file, a reaction field, shift function",
 392         "or user function can be used. For the user function a table file",
 393         "can be specified with the option [TT]-table[tt].",
 394         "The group of solvent molecules should be continuous and all molecules",
 395         "should have the same number of atoms.",
 396         "The user should add the ion molecules to the topology file or use",
 397         "the [TT]-p[tt] option to automatically modify the topology.[PAR]",
 398         "The ion molecule type, residue and atom names in all force fields",
 399         "are the capitalized element names without sign. This molecule name",
 400         "should be given with [TT]-pname[tt] or [TT]-nname[tt], and the",
 401         "[TT][molecules][tt] section of your topology updated accordingly,",
 402         "either by hand or with [TT]-p[tt]. Do not use an atom name instead!",
 403         "[PAR]Ions which can have multiple charge states get the multiplicity",
 404         "added, without sign, for the uncommon states only.[PAR]",
 405         "With the option [TT]-pot[tt] the potential can be written as B-factors",
 406         "in a [TT].pdb[tt] file (for visualisation using e.g. Rasmol).",
 407         "The unit of the potential is 1000 kJ/(mol e), the scaling be changed",
 408         "with the [TT]-scale[tt] option.[PAR]",
 409         "For larger ions, e.g. sulfate we recommended using [TT]genbox[tt]."
 410     };
 411     const char        *bugs[] = {
 412         "Calculation of the potential is not reliable, therefore the [TT]-random[tt] option is now turned on by default.",
 413         "If you specify a salt concentration existing ions are not taken into account. In effect you therefore specify the amount of salt to be added."
 414     };
 415     static int         p_num   = 0, n_num = 0, p_q = 1, n_q = -1;
 416     static const char *p_name  = "NA", *n_name = "CL";
 417     static real        rmin    = 0.6, scale = 0.001, conc = 0;
 418     static int         seed    = 1993;
 419     static gmx_bool    bRandom = TRUE, bNeutral = FALSE;
 420     static t_pargs     pa[]    = {
 421         { "-np",    FALSE, etINT,  {&p_num}, "Number of positive ions"       },
 422         { "-pname", FALSE, etSTR,  {&p_name}, "Name of the positive ion"      },
 423         { "-pq",    FALSE, etINT,  {&p_q},   "Charge of the positive ion"    },
 424         { "-nn",    FALSE, etINT,  {&n_num}, "Number of negative ions"       },
 425         { "-nname", FALSE, etSTR,  {&n_name}, "Name of the negative ion"      },
 426         { "-nq",    FALSE, etINT,  {&n_q},   "Charge of the negative ion"    },
 427         { "-rmin",  FALSE, etREAL, {&rmin},  "Minimum distance between ions" },
 428         { "-random", FALSE, etBOOL, {&bRandom}, "Use random placement of ions instead of based on potential. The rmin option should still work" },
 429         { "-seed",  FALSE, etINT,  {&seed},  "Seed for random number generator" },
 430         { "-scale", FALSE, etREAL, {&scale}, "Scaling factor for the potential for [TT]-pot[tt]" },
 431         { "-conc",  FALSE, etREAL, {&conc},
 432           "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." },
 433         { "-neutral", FALSE, etBOOL, {&bNeutral},
 434           "This option will add enough ions to neutralize the system. In combination with the concentration option a neutral system at a given salt concentration will be generated." }
 435     };
 436     gmx_mtop_t        *mtop;
 437     gmx_localtop_t    *top;
 438     t_inputrec         inputrec;
 439     t_commrec         *cr;
 440     t_mdatoms         *mdatoms;
 441     gmx_enerdata_t     enerd;
 442     t_graph           *graph;
 443     t_forcerec        *fr;
 444     rvec              *x, *v;
 445     real              *pot, vol, qtot;
 446     matrix             box;
 447     t_atoms            atoms;
 448     t_pbc              pbc;
 449     int               *repl;
 450     atom_id           *index;
 451     char              *grpname;
 452     gmx_bool          *bSet, bPDB;
 453     int                i, nw, nwa, nsa, nsalt, iqtot;
 454     FILE              *fplog;
 455     output_env_t       oenv;
 456     t_filenm           fnm[] = {
 457         { efTPX, NULL,  NULL,      ffREAD  },
 458         { efXVG, "-table", "table", ffOPTRD },
 459         { efNDX, NULL,  NULL,      ffOPTRD },
 460         { efSTO, "-o",  NULL,      ffWRITE },
 461         { efLOG, "-g",  "genion",  ffWRITE },
 462         { efPDB, "-pot", "pot",    ffOPTWR },
 463         { efTOP, "-p",  "topol",   ffOPTRW }
 464     };
 465 #define NFILE asize(fnm)
 466
 467     parse_common_args(&argc, argv, PCA_BE_NICE, NFILE, fnm, asize(pa), pa,
 468                       asize(desc), desc, asize(bugs), bugs, &oenv);
 469     bPDB = ftp2bSet(efPDB, NFILE, fnm);
 470     if (bRandom && bPDB)
 471     {
 472         fprintf(stderr, "Not computing potential with random option!\n");
 473         bPDB = FALSE;
 474     }
 475
 476     /* Check input for something sensible */
 477     if ((p_num < 0) || (n_num < 0))
 478     {
 479         gmx_fatal(FARGS, "Negative number of ions to add?");
 480     }
 481
 482     snew(mtop, 1);
 483     snew(top, 1);
 484     fplog = init_calcpot(ftp2fn(efLOG, NFILE, fnm), ftp2fn(efTPX, NFILE, fnm),
 485                          opt2fn("-table", NFILE, fnm), mtop, top, &inputrec, &cr,
 486                          &graph, &mdatoms, &fr, &enerd, &pot, box, &x, oenv);
 487
 488     atoms = gmx_mtop_global_atoms(mtop);
 489
 490     qtot = 0;
 491     for (i = 0; (i < atoms.nr); i++)
 492     {
 493         qtot += atoms.atom[i].q;
 494     }
 495     iqtot = gmx_nint(qtot);
 496
 497     if ((conc > 0) || bNeutral)
 498     {
 499         /* Compute number of ions to be added */
 500         vol = det(box);
 501         if (conc > 0)
 502         {
 503             nsalt = gmx_nint(conc*vol*AVOGADRO/1e24);
 504             p_num = abs(nsalt*n_q);
 505             n_num = abs(nsalt*p_q);
 506             if (bNeutral)
 507             {
 508                 int qdelta = 0;
 509                 do
 510                 {
 511                     qdelta = (p_num*p_q + n_num*n_q + iqtot);
 512                     if (qdelta < 0)
 513                     {
 514                         p_num  += abs(qdelta/p_q);
 515                         qdelta  = (p_num*p_q + n_num*n_q + iqtot);
 516                     }
 517                     if (qdelta > 0)
 518                     {
 519                         n_num  += abs(qdelta/n_q);
 520                         qdelta  = (p_num*p_q + n_num*n_q + iqtot);
 521                     }
 522                 }
 523                 while (qdelta != 0);
 524             }
 525         }
 526     }
 527
 528     if ((p_num == 0) && (n_num == 0))
 529     {
 530         if (!bPDB)
 531         {
 532             fprintf(stderr, "No ions to add and no potential to calculate.\n");
 533             exit(0);
 534         }
 535         nw  = 0;
 536         nsa = 0; /* to keep gcc happy */
 537     }
 538     else
 539     {
 540         printf("Will try to add %d %s ions and %d %s ions.\n",
 541                p_num, p_name, n_num, n_name);
 542         printf("Select a continuous group of solvent molecules\n");
 543         get_index(&atoms, ftp2fn_null(efNDX, NFILE, fnm), 1, &nwa, &index, &grpname);
 544         for (i = 1; i < nwa; i++)
 545         {
 546             if (index[i] != index[i-1]+1)
 547             {
 548                 gmx_fatal(FARGS, "The solvent group %s is not continuous: "
 549                           "index[%d]=%d, index[%d]=%d",
 550                           grpname, i, index[i-1]+1, i+1, index[i]+1);
 551             }
 552         }
 553         nsa = 1;
 554         while ((nsa < nwa) &&
 555                (atoms.atom[index[nsa]].resind ==
 556                 atoms.atom[index[nsa-1]].resind))
 557         {
 558             nsa++;
 559         }
 560         if (nwa % nsa)
 561         {
 562             gmx_fatal(FARGS, "Your solvent group size (%d) is not a multiple of %d",
 563                       nwa, nsa);
 564         }
 565         nw = nwa/nsa;
 566         fprintf(stderr, "Number of (%d-atomic) solvent molecules: %d\n", nsa, nw);
 567         if (p_num+n_num > nw)
 568         {
 569             gmx_fatal(FARGS, "Not enough solvent for adding ions");
 570         }
 571     }
 572
 573     if (opt2bSet("-p", NFILE, fnm))
 574     {
 575         update_topol(opt2fn("-p", NFILE, fnm), p_num, n_num, p_name, n_name, grpname);
 576     }
 577
 578     snew(bSet, nw);
 579     snew(repl, nw);
 580
 581     snew(v, atoms.nr);
 582     snew(atoms.pdbinfo, atoms.nr);
 583
 584     set_pbc(&pbc, inputrec.ePBC, box);
 585
 586     /* Now loop over the ions that have to be placed */
 587     do
 588     {
 589         if (!bRandom)
 590         {
 591             calc_pot(fplog, cr, mtop, &inputrec, top, x, fr, &enerd, mdatoms, pot, box, graph);
 592             if (bPDB || debug)
 593             {
 594                 char buf[STRLEN];
 595
 596                 if (debug)
 597                 {
 598                     sprintf(buf, "%d_%s", p_num+n_num, ftp2fn(efPDB, NFILE, fnm));
 599                 }
 600                 else
 601                 {
 602                     strcpy(buf, ftp2fn(efPDB, NFILE, fnm));
 603                 }
 604                 for (i = 0; (i < atoms.nr); i++)
 605                 {
 606                     atoms.pdbinfo[i].bfac = pot[i]*scale;
 607                 }
 608                 write_sto_conf(buf, "Potential calculated by genion",
 609                                &atoms, x, v, inputrec.ePBC, box);
 610                 bPDB = FALSE;
 611             }
 612         }
 613         if ((p_num > 0) && (p_num >= n_num))
 614         {
 615             insert_ion(nsa, &nw, bSet, repl, index, pot, x, &pbc,
 616                        1, p_q, p_name, mdatoms, rmin, bRandom, &seed);
 617             p_num--;
 618         }
 619         else if (n_num > 0)
 620         {
 621             insert_ion(nsa, &nw, bSet, repl, index, pot, x, &pbc,
 622                        -1, n_q, n_name, mdatoms, rmin, bRandom, &seed);
 623             n_num--;
 624         }
 625     }
 626     while (p_num+n_num > 0);
 627     fprintf(stderr, "\n");
 628
 629     if (nw)
 630     {
 631         sort_ions(nsa, nw, repl, index, &atoms, x, p_name, n_name);
 632     }
 633
 634     sfree(atoms.pdbinfo);
 635     atoms.pdbinfo = NULL;
 636     write_sto_conf(ftp2fn(efSTO, NFILE, fnm), *mtop->name, &atoms, x, NULL,
 637                    inputrec.ePBC, box);
 638
 639     thanx(stderr);
 640
 641     gmx_log_close(fplog);
 642
 643     return 0;
 644 }