Sort all includes in src/gromacs
[alexxy/gromacs.git] / src / gromacs / gmxpreprocess / solvate.cpp
index 3f2dd781350bbf76f64bec89e8cc383860978530..5cfbe5d03ec22cebb258509e0f01b864ba7d099d 100644 (file)
  * To help us fund GROMACS development, we humbly ask that you cite
  * the research papers on the package. Check out http://www.gromacs.org.
  */
-#include "solvate.h"
+#include "gmxpre.h"
 
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
+#include "solvate.h"
 
+#include <string.h>
 
-#include "sysstuff.h"
-#include "typedefs.h"
-#include "smalloc.h"
-#include "string2.h"
-#include "gromacs/fileio/confio.h"
-#include "macros.h"
-#include "gromacs/fileio/futil.h"
-#include "atomprop.h"
-#include "names.h"
-#include "vec.h"
-#include "gmx_fatal.h"
 #include "gromacs/commandline/pargs.h"
-#include "gromacs/gmxlib/conformation-utilities.h"
-#include "addconf.h"
-#include "read-conformation.h"
+#include "gromacs/fileio/confio.h"
 #include "gromacs/fileio/pdbio.h"
-#include "pbc.h"
+#include "gromacs/gmxlib/conformation-utilities.h"
+#include "gromacs/gmxpreprocess/addconf.h"
+#include "gromacs/gmxpreprocess/read-conformation.h"
+#include "gromacs/legacyheaders/macros.h"
+#include "gromacs/legacyheaders/names.h"
+#include "gromacs/legacyheaders/typedefs.h"
+#include "gromacs/math/vec.h"
+#include "gromacs/topology/atomprop.h"
+#include "gromacs/utility/cstringutil.h"
+#include "gromacs/utility/fatalerror.h"
+#include "gromacs/utility/futil.h"
+#include "gromacs/utility/smalloc.h"
 
 #ifdef DEBUG
 static void print_stat(rvec *x, int natoms, matrix box)
@@ -116,8 +113,10 @@ static void sort_molecule(t_atoms **atoms_solvt, rvec *x, rvec *v, real *r)
             moltp = NOTSET;
             for (j = 0; (j < nrmoltypes) && (moltp == NOTSET); j++)
             {
-                if (strcmp(*(atoms->resinfo[atoms->atom[i].resind].name),
-                           moltypes[j].name) == 0)
+                /* cppcheck-suppress nullPointer
+                 * moltypes is guaranteed to be allocated because otherwise
+                 * nrmoltypes is 0. */
+                if (strcmp(*(atoms->resinfo[atoms->atom[i].resind].name), moltypes[j].name) == 0)
                 {
                     moltp = j;
                 }
@@ -369,10 +368,10 @@ static void make_new_conformation(t_atoms *atoms, rvec *x, rvec *v, real *r, mat
     }
 }
 
-static void add_solv(const char *fn, t_atoms *atoms, rvec **x, rvec **v, real **r,
+static void add_solv(const char *fn, t_atoms *atoms, rvec **x, rvec **v, real **exclusionDistances,
                      int ePBC, matrix box,
                      gmx_atomprop_t aps,
-                     real r_distance, real r_scale, int *atoms_added,
+                     real defaultDistance, real scaleFactor, int *atoms_added,
                      int *residues_added, real rshell, int max_sol,
                      const output_env_t oenv)
 {
@@ -382,7 +381,7 @@ static void add_solv(const char *fn, t_atoms *atoms, rvec **x, rvec **v, real **
     char     title_solvt[STRLEN];
     t_atoms *atoms_solvt;
     rvec    *x_solvt, *v_solvt = NULL;
-    real    *r_solvt;
+    real    *exclusionDistances_solvt;
     int      ePBC_solvt;
     matrix   box_solvt;
     int      onr, onres;
@@ -391,18 +390,22 @@ static void add_solv(const char *fn, t_atoms *atoms, rvec **x, rvec **v, real **
     lfn = gmxlibfn(fn);
     strncpy(filename, lfn, STRLEN);
     sfree(lfn);
-    snew(atoms_solvt, 1);
-    get_stx_coordnum(filename, &(atoms_solvt->nr));
-    if (atoms_solvt->nr == 0)
     {
-        gmx_fatal(FARGS, "No solvent in %s, please check your input\n", filename);
+        int natoms;
+        get_stx_coordnum(filename, &natoms);
+        if (0 == natoms)
+        {
+            gmx_fatal(FARGS, "No solvent in %s, please check your input\n", filename);
+        }
+        snew(atoms_solvt, 1);
+        init_t_atoms(atoms_solvt, natoms, FALSE);
     }
     snew(x_solvt, atoms_solvt->nr);
     if (v)
     {
         snew(v_solvt, atoms_solvt->nr);
     }
-    snew(r_solvt, atoms_solvt->nr);
+    snew(exclusionDistances_solvt, atoms_solvt->nr);
     snew(atoms_solvt->resinfo, atoms_solvt->nr);
     snew(atoms_solvt->atomname, atoms_solvt->nr);
     snew(atoms_solvt->atom, atoms_solvt->nr);
@@ -419,8 +422,8 @@ static void add_solv(const char *fn, t_atoms *atoms, rvec **x, rvec **v, real **
     /* apply pbc for solvent configuration for whole molecules */
     rm_res_pbc(atoms_solvt, x_solvt, box_solvt);
 
-    /* initialise van der waals arrays for solvent configuration */
-    mk_vdw(atoms_solvt, r_solvt, aps, r_distance, r_scale);
+    /* initialise distance arrays for solvent configuration */
+    exclusionDistances_solvt = makeExclusionDistances(atoms_solvt, aps, defaultDistance, scaleFactor);
 
     /* calculate the box multiplication factors n_box[0...DIM] */
     nmol = 1;
@@ -445,10 +448,10 @@ static void add_solv(const char *fn, t_atoms *atoms, rvec **x, rvec **v, real **
     {
         srenew(v_solvt, atoms_solvt->nr*nmol);
     }
-    srenew(r_solvt, atoms_solvt->nr*nmol);
+    srenew(exclusionDistances_solvt, atoms_solvt->nr*nmol);
 
     /* generate a new solvent configuration */
-    make_new_conformation(atoms_solvt, x_solvt, v_solvt, r_solvt, box_solvt, n_box);
+    make_new_conformation(atoms_solvt, x_solvt, v_solvt, exclusionDistances_solvt, box_solvt, n_box);
 
 #ifdef DEBUG
     print_stat(x_solvt, atoms_solvt->nr, box_solvt);
@@ -458,18 +461,20 @@ static void add_solv(const char *fn, t_atoms *atoms, rvec **x, rvec **v, real **
     print_stat(x_solvt, atoms_solvt->nr, box_solvt);
 #endif
     /* Sort the solvent mixture, not the protein... */
-    sort_molecule(&atoms_solvt, x_solvt, v_solvt, r_solvt);
+    sort_molecule(&atoms_solvt, x_solvt, v_solvt, exclusionDistances_solvt);
 
     /* add the two configurations */
     onr   = atoms->nr;
     onres = atoms->nres;
-    add_conf(atoms, x, v, r, TRUE, ePBC, box, FALSE,
-             atoms_solvt, x_solvt, v_solvt, r_solvt, TRUE, rshell, max_sol, oenv);
+    add_conf(atoms, x, v, exclusionDistances, TRUE, ePBC, box, FALSE,
+             atoms_solvt, x_solvt, v_solvt, exclusionDistances_solvt, TRUE, rshell, max_sol, oenv);
     *atoms_added    = atoms->nr-onr;
     *residues_added = atoms->nres-onres;
 
     sfree(x_solvt);
-    sfree(r_solvt);
+    sfree(exclusionDistances_solvt);
+    done_atom(atoms_solvt);
+    sfree(atoms_solvt);
 
     fprintf(stderr, "Generated solvent containing %d atoms in %d residues\n",
             *atoms_added, *residues_added);
@@ -519,8 +524,8 @@ static void update_top(t_atoms *atoms, matrix box, int NFILE, t_filenm fnm[],
     if (ftp2bSet(efTOP, NFILE, fnm) )
     {
         fprintf(stderr, "Processing topology\n");
-        fpin    = ffopen(topinout, "r");
-        fpout   = ffopen(TEMP_FILENM, "w");
+        fpin    = gmx_ffopen(topinout, "r");
+        fpout   = gmx_ffopen(TEMP_FILENM, "w");
         line    = 0;
         bSystem = bMolecules = FALSE;
         while (fgets(buf, STRLEN, fpin))
@@ -589,17 +594,17 @@ static void update_top(t_atoms *atoms, matrix box, int NFILE, t_filenm fnm[],
                 fprintf(fpout, "%s", buf);
             }
         }
-        ffclose(fpin);
+        gmx_ffclose(fpin);
         if (nsol)
         {
             fprintf(stdout, "Adding line for %d solvent molecules to "
                     "topology file (%s)\n", nsol, topinout);
             fprintf(fpout, "%-15s %5d\n", "SOL", nsol);
         }
-        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
@@ -623,13 +628,11 @@ int gmx_solvate(int argc, char *argv[])
         "to change the box dimensions and center the solute.",
         "Solvent molecules are removed from the box where the ",
         "distance between any atom of the solute molecule(s) and any atom of ",
-        "the solvent molecule is less than the sum of the van der Waals radii of ",
-        "both atoms. A database ([TT]vdwradii.dat[tt]) of van der Waals radii is ",
-        "read by the program, and atoms not in the database are ",
-        "assigned a default distance [TT]-vdwd[tt].",
-        "Note that this option will also influence the distances between",
-        "solvent molecules if they contain atoms that are not in the database.",
-        "[PAR]",
+        "the solvent molecule is less than the sum of the scaled van der Waals",
+        "radii of both atoms. A database ([TT]vdwradii.dat[tt]) of van der",
+        "Waals radii is read by the program, and the resulting radii scaled",
+        "by [TT]-scale[tt]. If radii are not found in the database, those"
+        "atoms are assigned the (pre-scaled) distance [TT]-radius[tt].[PAR]",
 
         "The default solvent is Simple Point Charge water (SPC), with coordinates ",
         "from [TT]$GMXLIB/spc216.gro[tt]. These coordinates can also be used",
@@ -675,13 +678,13 @@ int gmx_solvate(int argc, char *argv[])
     /* parameter data */
     gmx_bool       bProt, bBox;
     const char    *conf_prot, *confout;
-    real          *r;
+    real          *exclusionDistances = NULL;
     gmx_atomprop_t aps;
 
     /* protein configuration data */
     char    *title = NULL;
-    t_atoms  atoms;
-    rvec    *x, *v = NULL;
+    t_atoms *atoms;
+    rvec    *x    = NULL, *v = NULL;
     int      ePBC = -1;
     matrix   box;
 
@@ -696,18 +699,18 @@ int gmx_solvate(int argc, char *argv[])
     };
 #define NFILE asize(fnm)
 
-    static real     r_distance = 0.105, r_shell = 0, r_scale = 0.57;
-    static rvec     new_box    = {0.0, 0.0, 0.0};
-    static gmx_bool bReadV     = FALSE;
-    static int      max_sol    = 0;
+    static real     defaultDistance = 0.105, r_shell = 0, scaleFactor = 0.57;
+    static rvec     new_box         = {0.0, 0.0, 0.0};
+    static gmx_bool bReadV          = FALSE;
+    static int      max_sol         = 0;
     output_env_t    oenv;
     t_pargs         pa[]              = {
         { "-box",    FALSE, etRVEC, {new_box},
           "Box size (in nm)" },
-        { "-vdwd",   FALSE, etREAL, {&r_distance},
+        { "-radius",   FALSE, etREAL, {&defaultDistance},
           "Default van der Waals distance"},
-        { "-vdwscale", FALSE, etREAL, {&r_scale},
-          "HIDDENScale factor to multiply Van der Waals radii from the database in share/gromacs/top/vdwradii.dat. The default value of 0.57 yields density close to 1000 g/l for proteins in water." },
+        { "-scale", FALSE, etREAL, {&scaleFactor},
+          "Scale factor to multiply Van der Waals radii from the database in share/gromacs/top/vdwradii.dat. The default value of 0.57 yields density close to 1000 g/l for proteins in water." },
         { "-shell",  FALSE, etREAL, {&r_shell},
           "Thickness of optional water layer around solute" },
         { "-maxsol", FALSE, etINT,  {&max_sol},
@@ -716,7 +719,7 @@ int gmx_solvate(int argc, char *argv[])
           "Keep velocities from input solute and solvent" },
     };
 
-    if (!parse_common_args(&argc, argv, PCA_BE_NICE, NFILE, fnm, asize(pa), pa,
+    if (!parse_common_args(&argc, argv, 0, NFILE, fnm, asize(pa), pa,
                            asize(desc), desc, asize(bugs), bugs, &oenv))
     {
         return 0;
@@ -735,33 +738,26 @@ int gmx_solvate(int argc, char *argv[])
 
     aps = gmx_atomprop_init();
 
+    snew(atoms, 1);
+    init_t_atoms(atoms, 0, FALSE);
     if (bProt)
     {
         /* Generate a solute configuration */
         conf_prot = opt2fn("-cp", NFILE, fnm);
-        title     = read_conformation(conf_prot, &atoms, &x, bReadV ? &v : NULL, &r, &ePBC, box,
-                                      aps, r_distance, r_scale);
+        title     = readConformation(conf_prot, atoms, &x,
+                                     bReadV ? &v : NULL, &ePBC, box);
+        exclusionDistances = makeExclusionDistances(atoms, aps, defaultDistance, scaleFactor);
+
         if (bReadV && !v)
         {
             fprintf(stderr, "Note: no velocities found\n");
         }
-        if (atoms.nr == 0)
+        if (atoms->nr == 0)
         {
             fprintf(stderr, "Note: no atoms in %s\n", conf_prot);
             bProt = FALSE;
         }
     }
-    else
-    {
-        atoms.nr       = 0;
-        atoms.nres     = 0;
-        atoms.resinfo  = NULL;
-        atoms.atomname = NULL;
-        atoms.atom     = NULL;
-        atoms.pdbinfo  = NULL;
-        x              = NULL;
-        r              = NULL;
-    }
     if (bBox)
     {
         ePBC = epbcXYZ;
@@ -776,8 +772,8 @@ int gmx_solvate(int argc, char *argv[])
                   "or give explicit -box command line option");
     }
 
-    add_solv(solventFileName, &atoms, &x, v ? &v : NULL, &r, ePBC, box,
-             aps, r_distance, r_scale, &atoms_added, &residues_added, r_shell, max_sol,
+    add_solv(solventFileName, atoms, &x, v ? &v : NULL, &exclusionDistances, ePBC, box,
+             aps, defaultDistance, scaleFactor, &atoms_added, &residues_added, r_shell, max_sol,
              oenv);
 
     /* write new configuration 1 to file confout */
@@ -785,21 +781,26 @@ int gmx_solvate(int argc, char *argv[])
     fprintf(stderr, "Writing generated configuration to %s\n", confout);
     if (bProt)
     {
-        write_sto_conf(confout, title, &atoms, x, v, ePBC, box);
+        write_sto_conf(confout, title, atoms, x, v, ePBC, box);
         /* print box sizes and box type to stderr */
         fprintf(stderr, "%s\n", title);
     }
     else
     {
-        write_sto_conf(confout, "Generated by gmx solvate", &atoms, x, v, ePBC, box);
+        write_sto_conf(confout, "Generated by gmx solvate", atoms, x, v, ePBC, box);
     }
 
     /* print size of generated configuration */
     fprintf(stderr, "\nOutput configuration contains %d atoms in %d residues\n",
-            atoms.nr, atoms.nres);
-    update_top(&atoms, box, NFILE, fnm, aps);
+            atoms->nr, atoms->nres);
+    update_top(atoms, box, NFILE, fnm, aps);
 
     gmx_atomprop_destroy(aps);
+    sfree(exclusionDistances);
+    sfree(x);
+    sfree(v);
+    done_atom(atoms);
+    sfree(atoms);
 
     return 0;
 }