Remove .tpa, .tpb, .tpx, .trj files. Part of #1500.
[alexxy/gromacs.git] / src / gromacs / gmxana / gmx_genion.c
index 8c521daae19228d63d81cb4cd64709c49922e928..1301d8e01a92c7de5b5990531472fcb3c2d5c25d 100644 (file)
@@ -1,85 +1,73 @@
 /*
+ * This file is part of the GROMACS molecular simulation package.
  *
- *                This source code is part of
- *
- *                 G   R   O   M   A   C   S
- *
- *          GROningen MAchine for Chemical Simulations
- *
- *                        VERSION 3.2.0
- * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
- * Copyright (c) 2001-2004, The GROMACS development team,
- * check out http://www.gromacs.org for more information.
-
- * This program is free software; you can redistribute it and/or
- * modify it under the terms of the GNU General Public License
- * as published by the Free Software Foundation; either version 2
+ * Copyright (c) 2001-2004, The GROMACS development team.
+ * Copyright (c) 2013,2014, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
  * of the License, or (at your option) any later version.
  *
- * If you want to redistribute modifications, please consider that
- * scientific software is very special. Version control is crucial -
- * bugs must be traceable. We will be happy to consider code for
- * inclusion in the official distribution, but derived work must not
- * be called official GROMACS. Details are found in the README & COPYING
- * files - if they are missing, get the official version at www.gromacs.org.
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * Lesser General Public License for more details.
  *
- * To help us fund GROMACS development, we humbly ask that you cite
- * the papers on the package - you can find them in the top README file.
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  *
- * For more info, check our website at http://www.gromacs.org
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
  *
- * And Hey:
- * Green Red Orange Magenta Azure Cyan Skyblue
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
  */
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
+#include "gmxpre.h"
 
 #include <ctype.h>
-#include "string2.h"
-#include "smalloc.h"
-#include "sysstuff.h"
-#include "confio.h"
-#include "statutil.h"
-#include "pbc.h"
-#include "force.h"
-#include "gmx_fatal.h"
-#include "futil.h"
-#include "maths.h"
-#include "macros.h"
-#include "vec.h"
-#include "tpxio.h"
-#include "mdrun.h"
-#include "main.h"
-#include "random.h"
-#include "index.h"
-#include "mtop_util.h"
+#include <stdlib.h>
+#include <string.h>
+
+#include "gromacs/utility/cstringutil.h"
+#include "gromacs/utility/smalloc.h"
+#include "gromacs/fileio/confio.h"
+#include "gromacs/commandline/pargs.h"
+#include "gromacs/pbcutil/pbc.h"
+#include "gromacs/legacyheaders/force.h"
+#include "gromacs/utility/fatalerror.h"
+#include "gromacs/utility/futil.h"
+#include "gromacs/math/utilities.h"
+#include "gromacs/legacyheaders/macros.h"
+#include "gromacs/math/vec.h"
+#include "gromacs/fileio/tpxio.h"
+#include "gromacs/legacyheaders/mdrun.h"
+#include "gromacs/random/random.h"
+#include "gromacs/topology/index.h"
 #include "gmx_ana.h"
 
-static int greatest_common_divisor(int p, int q)
-{
-    int tmp;
-    while (q != 0)
-    {
-        tmp = q;
-        q = p % q;
-        p = tmp;
-    }
-    return p;
-}
-
 static void insert_ion(int nsa, int *nwater,
                        gmx_bool bSet[], int repl[], atom_id index[],
                        rvec x[], t_pbc *pbc,
                        int sign, int q, const char *ionname,
                        t_atoms *atoms,
-                       real rmin, int *seed)
+                       real rmin, gmx_rng_t rng)
 {
-    int             i, ei,nw;
+    int             i, ei, nw;
     real            rmin2;
     rvec            dx;
-    gmx_large_int_t maxrand;
+    gmx_int64_t     maxrand;
 
     ei       = -1;
     nw       = *nwater;
@@ -88,7 +76,7 @@ static void insert_ion(int nsa, int *nwater,
 
     do
     {
-        ei = nw*rando(seed);
+        ei = nw*gmx_rng_uniform_real(rng);
         maxrand--;
     }
     while (bSet[ei] && (maxrand > 0));
@@ -134,7 +122,7 @@ static char *aname(const char *mname)
     char *str;
     int   i;
 
-    str = strdup(mname);
+    str = gmx_strdup(mname);
     i   = strlen(str)-1;
     while (i > 1 && (isdigit(str[i]) || (str[i] == '+') || (str[i] == '-')))
     {
@@ -188,14 +176,14 @@ void sort_ions(int nsa, int nw, int repl[], atom_id index[],
         if (np)
         {
             snew(pptr, 1);
-            pptr[0] = strdup(p_name);
+            pptr[0] = gmx_strdup(p_name);
             snew(paptr, 1);
             paptr[0] = aname(p_name);
         }
         if (nn)
         {
             snew(nptr, 1);
-            nptr[0] = strdup(n_name);
+            nptr[0] = gmx_strdup(n_name);
             snew(naptr, 1);
             naptr[0] = aname(n_name);
         }
@@ -253,8 +241,8 @@ static void update_topol(const char *topinout, int p_num, int n_num,
     gmx_bool bMolecules;
 
     printf("\nProcessing topology\n");
-    fpin  = ffopen(topinout, "r");
-    fpout = ffopen(TEMP_FILENM, "w");
+    fpin  = gmx_ffopen(topinout, "r");
+    fpout = gmx_ffopen(TEMP_FILENM, "w");
 
     line       = 0;
     bMolecules = FALSE;
@@ -302,20 +290,20 @@ static void update_topol(const char *topinout, int p_num, int n_num,
             }
             /* Store this molecules section line */
             srenew(mol_line, nmol_line+1);
-            mol_line[nmol_line] = strdup(buf);
+            mol_line[nmol_line] = gmx_strdup(buf);
             nmol_line++;
         }
     }
-    ffclose(fpin);
+    gmx_ffclose(fpin);
 
     if (sol_line == -1)
     {
-        ffclose(fpout);
+        gmx_ffclose(fpout);
         gmx_fatal(FARGS, "No line with moleculetype '%s' found the [ molecules ] section of file '%s'", grpname, topinout);
     }
     if (nsol_last < p_num+n_num)
     {
-        ffclose(fpout);
+        gmx_ffclose(fpout);
         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);
     }
 
@@ -346,10 +334,10 @@ static void update_topol(const char *topinout, int p_num, int n_num,
             }
         }
     }
-    ffclose(fpout);
-    /* use ffopen to generate backup of topinout */
-    fpout = ffopen(topinout, "w");
-    ffclose(fpout);
+    gmx_ffclose(fpout);
+    /* use gmx_ffopen to generate backup of topinout */
+    fpout = gmx_ffopen(topinout, "w");
+    gmx_ffclose(fpout);
     rename(TEMP_FILENM, topinout);
 #undef TEMP_FILENM
 }
@@ -357,7 +345,7 @@ static void update_topol(const char *topinout, int p_num, int n_num,
 int gmx_genion(int argc, char *argv[])
 {
     const char        *desc[] = {
-        "[TT]genion[tt] randomly replaces solvent molecules with monoatomic ions.",
+        "[THISMODULE] randomly replaces solvent molecules with monoatomic ions.",
         "The group of solvent molecules should be continuous and all molecules",
         "should have the same number of atoms.",
         "The user should add the ion molecules to the topology file or use",
@@ -369,18 +357,18 @@ int gmx_genion(int argc, char *argv[])
         "either by hand or with [TT]-p[tt]. Do not use an atom name instead!",
         "[PAR]Ions which can have multiple charge states get the multiplicity",
         "added, without sign, for the uncommon states only.[PAR]",
-        "For larger ions, e.g. sulfate we recommended using [TT]genbox[tt]."
+        "For larger ions, e.g. sulfate we recommended using [gmx-insert-molecules]."
     };
     const char        *bugs[] = {
         "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.",
     };
-    static int         p_num   = 0, n_num = 0, p_q = 1, n_q = -1;
-    static const char *p_name  = "NA", *n_name = "CL";
-    static real        rmin    = 0.6, conc = 0;
-    static int         seed    = 1993;
+    static int         p_num    = 0, n_num = 0, p_q = 1, n_q = -1;
+    static const char *p_name   = "NA", *n_name = "CL";
+    static real        rmin     = 0.6, conc = 0;
+    static int         seed     = 1993;
     static gmx_bool    bNeutral = FALSE;
-    static t_pargs     pa[]    = {
+    static t_pargs     pa[]     = {
         { "-np",    FALSE, etINT,  {&p_num}, "Number of positive ions"       },
         { "-pname", FALSE, etSTR,  {&p_name}, "Name of the positive ion"      },
         { "-pq",    FALSE, etINT,  {&p_q},   "Charge of the positive ion"    },
@@ -393,7 +381,7 @@ int gmx_genion(int argc, char *argv[])
           "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." },
         { "-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]. "}
     };
-    t_topology        top;
+    t_topology         top;
     rvec              *x, *v;
     real               vol, qtot;
     matrix             box;
@@ -405,16 +393,20 @@ int gmx_genion(int argc, char *argv[])
     gmx_bool          *bSet;
     int                i, nw, nwa, nsa, nsalt, iqtot;
     output_env_t       oenv;
+    gmx_rng_t          rng;
     t_filenm           fnm[] = {
-        { efTPX, NULL,  NULL,      ffREAD  },
+        { efTPR, NULL,  NULL,      ffREAD  },
         { efNDX, NULL,  NULL,      ffOPTRD },
         { efSTO, "-o",  NULL,      ffWRITE },
         { efTOP, "-p",  "topol",   ffOPTRW }
     };
 #define NFILE asize(fnm)
 
-    parse_common_args(&argc, argv, PCA_BE_NICE, NFILE, fnm, asize(pa), pa,
-                      asize(desc), desc, asize(bugs), bugs, &oenv);
+    if (!parse_common_args(&argc, argv, 0, NFILE, fnm, asize(pa), pa,
+                           asize(desc), desc, asize(bugs), bugs, &oenv))
+    {
+        return 0;
+    }
 
     /* Check input for something sensible */
     if ((p_num < 0) || (n_num < 0))
@@ -428,7 +420,7 @@ int gmx_genion(int argc, char *argv[])
     }
 
     /* Read atom positions and charges */
-    read_tps_conf(ftp2fn(efTPX, NFILE, fnm), title, &top, &ePBC, &x, &v, box, FALSE);
+    read_tps_conf(ftp2fn(efTPR, NFILE, fnm), title, &top, &ePBC, &x, &v, box, FALSE);
     atoms = top.atoms;
 
     /* Compute total charge */
@@ -439,11 +431,11 @@ int gmx_genion(int argc, char *argv[])
     }
     iqtot = gmx_nint(qtot);
 
-    
+
     if (conc > 0)
     {
         /* Compute number of ions to be added */
-        vol = det(box);
+        vol   = det(box);
         nsalt = gmx_nint(conc*vol*AVOGADRO/1e24);
         p_num = abs(nsalt*n_q);
         n_num = abs(nsalt*p_q);
@@ -454,13 +446,13 @@ int gmx_genion(int argc, char *argv[])
 
         /* Check if the system is neutralizable
          * is (qdelta == p_q*p_num + n_q*n_num) solvable for p_num and n_num? */
-        int gcd = greatest_common_divisor(n_q, p_q);
+        int gcd = gmx_greatest_common_divisor(n_q, p_q);
         if ((qdelta % gcd) != 0)
         {
             gmx_fatal(FARGS, "Can't neutralize this system using -nq %d and"
-                    " -pq %d.\n", n_q, p_q);
+                      " -pq %d.\n", n_q, p_q);
         }
-        
+
         while (qdelta != 0)
         {
             while (qdelta < 0)
@@ -478,8 +470,7 @@ int gmx_genion(int argc, char *argv[])
 
     if ((p_num == 0) && (n_num == 0))
     {
-        fprintf(stderr, "No ions to add.\n");
-        exit(0);
+        fprintf(stderr, "No ions to add, will just copy input configuration.\n");
     }
     else
     {
@@ -514,43 +505,51 @@ int gmx_genion(int argc, char *argv[])
         {
             gmx_fatal(FARGS, "Not enough solvent for adding ions");
         }
-    }
 
-    if (opt2bSet("-p", NFILE, fnm))
-    {
-        update_topol(opt2fn("-p", NFILE, fnm), p_num, n_num, p_name, n_name, grpname);
-    }
+        if (opt2bSet("-p", NFILE, fnm))
+        {
+            update_topol(opt2fn("-p", NFILE, fnm), p_num, n_num, p_name, n_name, grpname);
+        }
 
-    snew(bSet, nw);
-    snew(repl, nw);
+        snew(bSet, nw);
+        snew(repl, nw);
 
-    snew(v, atoms.nr);
-    snew(atoms.pdbinfo, atoms.nr);
+        snew(v, atoms.nr);
+        snew(atoms.pdbinfo, atoms.nr);
 
-    set_pbc(&pbc, ePBC, box);
+        set_pbc(&pbc, ePBC, box);
 
-    /* Now loop over the ions that have to be placed */
-    while (p_num-- > 0)
-    {
-        insert_ion(nsa, &nw, bSet, repl, index, x, &pbc,
-                   1, p_q, p_name, &atoms, rmin, &seed);
-    }
-    while (n_num-- > 0)
-    {
-        insert_ion(nsa, &nw, bSet, repl, index, x, &pbc,
-                   -1, n_q, n_name, &atoms, rmin, &seed);
-    }
-    fprintf(stderr, "\n");
+        if (seed == 0)
+        {
+            rng = gmx_rng_init(gmx_rng_make_seed());
+        }
+        else
+        {
+            rng = gmx_rng_init(seed);
+        }
+        /* Now loop over the ions that have to be placed */
+        while (p_num-- > 0)
+        {
+            insert_ion(nsa, &nw, bSet, repl, index, x, &pbc,
+                       1, p_q, p_name, &atoms, rmin, rng);
+        }
+        while (n_num-- > 0)
+        {
+            insert_ion(nsa, &nw, bSet, repl, index, x, &pbc,
+                       -1, n_q, n_name, &atoms, rmin, rng);
+        }
+        gmx_rng_destroy(rng);
+        fprintf(stderr, "\n");
 
-    if (nw)
-    {
-        sort_ions(nsa, nw, repl, index, &atoms, x, p_name, n_name);
-    }
+        if (nw)
+        {
+            sort_ions(nsa, nw, repl, index, &atoms, x, p_name, n_name);
+        }
 
-    sfree(atoms.pdbinfo);
-    atoms.pdbinfo = NULL;
-    write_sto_conf(ftp2fn(efSTO, NFILE, fnm), *top.name, &atoms, x, NULL, ePBC,
-                   box);
+        sfree(atoms.pdbinfo);
+        atoms.pdbinfo = NULL;
+    }
+    write_sto_conf(ftp2fn(efSTO, NFILE, fnm), *top.name, &atoms, x, NULL, ePBC, box);
 
     return 0;
 }