Merge branch 'master' into pygromacs
[alexxy/gromacs.git] / src / gromacs / gmxana / gmx_order.cpp
similarity index 95%
rename from src/gromacs/gmxana/gmx_order.c
rename to src/gromacs/gmxana/gmx_order.cpp
index c0dded8283de35e32d686594ece504779146a3eb..0b28c3bbd3d644c907a7b12cc490767b0a06776b 100644 (file)
  */
 #include "gmxpre.h"
 
-#include <math.h>
-#include <string.h>
+#include <cmath>
+#include <cstring>
+
+#include <algorithm>
 
 #include "gromacs/commandline/pargs.h"
 #include "gromacs/fileio/confio.h"
@@ -56,6 +58,7 @@
 #include "gromacs/utility/cstringutil.h"
 #include "gromacs/utility/fatalerror.h"
 #include "gromacs/utility/futil.h"
+#include "gromacs/utility/gmxassert.h"
 #include "gromacs/utility/smalloc.h"
 
 /****************************************************************************/
@@ -77,17 +80,13 @@ static void find_nearest_neighbours(int ePBC,
                                     real sgslice[], real skslice[],
                                     gmx_rmpbc_t gpbc)
 {
-    FILE    *fpoutdist;
-    char     fnsgdist[32];
     int      ix, jx, nsgbin, *sgbin;
-    int      i1, i2, i, ibin, j, k, l, n, *nn[4];
-    rvec     dx, dx1, dx2, rj, rk, urk, urj;
+    int      i, ibin, j, k, l, n, *nn[4];
+    rvec     dx, rj, rk, urk, urj;
     real     cost, cost2, *sgmol, *skmol, rmean, rmean2, r2, box2, *r_nn[4];
     t_pbc    pbc;
-    t_mat   *dmat;
-    t_dist  *d;
-    int      m1, mm, sl_index;
-    int    **nnb, *sl_count;
+    int      sl_index;
+    int     *sl_count;
     real     onethird = 1.0/3.0;
     /*  dmat = init_mat(maxidx, FALSE); */
     box2 = box[XX][XX] * box[XX][XX];
@@ -111,7 +110,7 @@ static void find_nearest_neighbours(int ePBC,
 
     gmx_rmpbc(gpbc, natoms, box, x);
 
-    nsgbin = 1 + 1/0.0005;
+    nsgbin = 2001; // Calculated as (1 + 1/0.0005)
     snew(sgbin, nsgbin);
 
     *sgmean = 0.0;
@@ -164,7 +163,7 @@ static void find_nearest_neighbours(int ePBC,
         rmean = 0;
         for (j = 0; (j < 4); j++)
         {
-            r_nn[j][i] = sqrt(r_nn[j][i]);
+            r_nn[j][i] = std::sqrt(r_nn[j][i]);
             rmean     += r_nn[j][i];
         }
         rmean /= 4;
@@ -192,7 +191,7 @@ static void find_nearest_neighbours(int ePBC,
                 sgmol[i] += cost2;
 
                 /* determine distribution */
-                ibin = nsgbin * cost2;
+                ibin = static_cast<int>(nsgbin * cost2);
                 if (ibin < nsgbin)
                 {
                     sgbin[ibin]++;
@@ -259,7 +258,7 @@ static void calc_tetra_order_parm(const char *fnNDX, const char *fnTPS,
     FILE        *fpsg = NULL, *fpsk = NULL;
     t_topology   top;
     int          ePBC;
-    char         title[STRLEN], fn[STRLEN], subtitle[STRLEN];
+    char         title[STRLEN];
     t_trxstatus *status;
     int          natoms;
     real         t;
@@ -396,12 +395,11 @@ void calc_order(const char *fn, atom_id *index, atom_id *a, rvec **order,
     int natoms,                                  /* nr. atoms in trj                               */
         nr_tails,                                /* nr tails, to check if index file is correct    */
         size = 0,                                /* nr. of atoms in group. same as nr_tails        */
-        i, j, m, k, l, teller = 0,
+        i, j, m, k, teller = 0,
         slice,                                   /* current slice number                           */
         nr_frames = 0;
     int         *slCount;                        /* nr. of atoms in one slice                      */
-    real         dbangle                = 0,     /* angle between double bond and  axis            */
-                 sdbangle               = 0;     /* sum of these angles                            */
+    real         sdbangle               = 0;     /* sum of these angles                            */
     gmx_bool     use_unitvector         = FALSE; /* use a specified unit vector instead of axis to specify unit normal*/
     rvec         direction, com, dref, dvec;
     int          comsize, distsize;
@@ -413,7 +411,7 @@ void calc_order(const char *fn, atom_id *index, atom_id *a, rvec **order,
 
     /* PBC added for center-of-mass vector*/
     /* Initiate the pbc structure */
-    memset(&pbc, 0, sizeof(pbc));
+    std::memset(&pbc, 0, sizeof(pbc));
 
     if ((natoms = read_first_x(oenv, &status, fn, &t, &x0, box)) == 0)
     {
@@ -567,7 +565,7 @@ void calc_order(const char *fn, atom_id *index, atom_id *a, rvec **order,
                     rvec_sub(x1[a[index[i+1]+j]], x1[a[index[i]+j]], dist);
                     length = norm(dist);
                     check_length(length, a[index[i]+j], a[index[i+1]+j]);
-                    svmul(1/length, dist, Sz);
+                    svmul(1.0/length, dist, Sz);
 
                     /* this is actually the cosine of the angle between the double bond
                        and axis, because Sz is normalized and the two other components of
@@ -578,7 +576,7 @@ void calc_order(const char *fn, atom_id *index, atom_id *a, rvec **order,
                     }
                     else
                     {
-                        sdbangle += acos(Sz[axis]);
+                        sdbangle += std::acos(Sz[axis]);
                     }
                 }
                 else
@@ -588,7 +586,7 @@ void calc_order(const char *fn, atom_id *index, atom_id *a, rvec **order,
                     length = norm(dist); /* determine distance between two atoms */
                     check_length(length, a[index[i-1]+j], a[index[i+1]+j]);
 
-                    svmul(1/length, dist, Sz);
+                    svmul(1.0/length, dist, Sz);
                     /* Sz is now the molecular axis Sz, normalized and all that */
                 }
 
@@ -597,11 +595,11 @@ void calc_order(const char *fn, atom_id *index, atom_id *a, rvec **order,
                 rvec_sub(x1[a[index[i+1]+j]], x1[a[index[i]+j]], tmp1);
                 rvec_sub(x1[a[index[i-1]+j]], x1[a[index[i]+j]], tmp2);
                 cprod(tmp1, tmp2, Sx);
-                svmul(1/norm(Sx), Sx, Sx);
+                svmul(1.0/norm(Sx), Sx, Sx);
 
                 /* now we can get Sy from the outer product of Sx and Sz   */
                 cprod(Sz, Sx, Sy);
-                svmul(1/norm(Sy), Sy, Sy);
+                svmul(1.0/norm(Sy), Sy, Sy);
 
                 /* the square of cosine of the angle between dist and the axis.
                    Using the innerproduct, but two of the three elements are zero
@@ -622,7 +620,7 @@ void calc_order(const char *fn, atom_id *index, atom_id *a, rvec **order,
 
                 for (m = 0; m < DIM; m++)
                 {
-                    frameorder[m] += 0.5 * (3 * cossum[m] - 1);
+                    frameorder[m] += 0.5 * (3.0 * cossum[m] - 1.0);
                 }
 
                 if (bSliced)
@@ -646,7 +644,7 @@ void calc_order(const char *fn, atom_id *index, atom_id *a, rvec **order,
                         z_ave -= box[axis][axis];
                     }
 
-                    slice  = (int)(0.5 + (z_ave / (*slWidth))) - 1;
+                    slice  = static_cast<int>((0.5 + (z_ave / (*slWidth))) - 1);
                     slCount[slice]++;     /* determine slice, increase count */
 
                     slFrameorder[slice] += 0.5 * (3 * cossum[axis] - 1);
@@ -656,7 +654,7 @@ void calc_order(const char *fn, atom_id *index, atom_id *a, rvec **order,
                     /*  store per-molecule order parameter
                      *  To just track single-axis order: (*slOrder)[j][i] += 0.5 * (3 * iprod(cossum,direction) - 1);
                      *  following is for Scd order: */
-                    (*slOrder)[j][i] += -1* (0.3333 * (3 * cossum[XX] - 1) + 0.3333 * 0.5 * (3 * cossum[YY] - 1));
+                    (*slOrder)[j][i] += -1* (1.0/3.0 * (3 * cossum[XX] - 1) + 1.0/3.0 * 0.5 * (3.0 * cossum[YY] - 1));
                 }
                 if (distcalc)
                 {
@@ -675,10 +673,10 @@ void calc_order(const char *fn, atom_id *index, atom_id *a, rvec **order,
                             pbc_dx(&pbc, x1[distidx[k]], x1[a[index[i]+j]], dvec);
                             /* at the moment, just remove dvec[axis] */
                             dvec[axis] = 0;
-                            tmpdist    = min(tmpdist, norm2(dvec));
+                            tmpdist    = std::min(tmpdist, norm2(dvec));
                         }
                         //fprintf(stderr, "Min dist %f; trace %f\n", tmpdist, trace(box));
-                        (*distvals)[j][i] += sqrt(tmpdist);
+                        (*distvals)[j][i] += std::sqrt(tmpdist);
                     }
                 }
             } /* end loop j, over all atoms in group */
@@ -771,8 +769,8 @@ void order_plot(rvec order[], real *slOrder[], const char *afile, const char *bf
         slOrd = xvgropen(bfile, buf, "Molecule", "S", oenv);
         for (atom = 1; atom < ngrps - 1; atom++)
         {
-            fprintf(ord, "%12d   %12g\n", atom, -1 * (0.6667 * order[atom][XX] +
-                                                      0.333 * order[atom][YY]));
+            fprintf(ord, "%12d   %12g\n", atom, -1.0 * (2.0/3.0 * order[atom][XX] +
+                                                        1.0/3.0 * order[atom][YY]));
         }
 
         for (slice = 0; slice < nslices; slice++)
@@ -826,8 +824,8 @@ void order_plot(rvec order[], real *slOrder[], const char *afile, const char *bf
         {
             fprintf(ord, "%12d   %12g   %12g   %12g\n", atom, order[atom][XX],
                     order[atom][YY], order[atom][ZZ]);
-            fprintf(slOrd, "%12d   %12g\n", atom, -1 * (0.6667 * order[atom][XX] +
-                                                        0.333 * order[atom][YY]));
+            fprintf(slOrd, "%12d   %12g\n", atom, -1.0 * (2.0/3.0 * order[atom][XX] +
+                                                          1.0/3.0 * order[atom][YY]));
         }
 
     }
@@ -840,7 +838,6 @@ void write_bfactors(t_filenm  *fnm, int nfile, atom_id *index, atom_id *a, int n
     /*function to write order parameters as B factors in PDB file using
           first frame of trajectory*/
     t_trxstatus *status;
-    int          natoms;
     t_trxframe   fr, frout;
     t_atoms      useatoms;
     int          i, j, ctr, nout;
@@ -848,8 +845,8 @@ void write_bfactors(t_filenm  *fnm, int nfile, atom_id *index, atom_id *a, int n
     ngrps -= 2;  /*we don't have an order parameter for the first or
                        last atom in each chain*/
     nout   = nslices*ngrps;
-    natoms = read_first_frame(oenv, &status, ftp2fn(efTRX, nfile, fnm), &fr,
-                              TRX_NEED_X);
+    read_first_frame(oenv, &status, ftp2fn(efTRX, nfile, fnm), &fr, TRX_NEED_X);
+
     close_trj(status);
     frout        = fr;
     frout.natoms = nout;
@@ -883,7 +880,7 @@ void write_bfactors(t_filenm  *fnm, int nfile, atom_id *index, atom_id *a, int n
             copy_rvec(fr.x[a[index[i+1]+j]], frout.x[ctr]);
             useatoms.atomname[ctr] = top->atoms.atomname[a[index[i+1]+j]];
             useatoms.atom[ctr]     = top->atoms.atom[a[index[i+1]+j]];
-            useatoms.nres          = max(useatoms.nres, useatoms.atom[ctr].resind+1);
+            useatoms.nres          = std::max(useatoms.nres, useatoms.atom[ctr].resind+1);
             useatoms.resinfo[useatoms.atom[ctr].resind] = top->atoms.resinfo[useatoms.atom[ctr].resind]; /*copy resinfo*/
         }
     }
@@ -991,15 +988,16 @@ int gmx_order(int argc, char *argv[])
     trxfnm = ftp2fn(efTRX, NFILE, fnm);
 
     /* Calculate axis */
-    if (strcmp(normal_axis[0], "x") == 0)
+    GMX_RELEASE_ASSERT(normal_axis[0] != NULL, "Options inconsistency; normal_axis[0] is NULL");
+    if (std::strcmp(normal_axis[0], "x") == 0)
     {
         axis = XX;
     }
-    else if (strcmp(normal_axis[0], "y") == 0)
+    else if (std::strcmp(normal_axis[0], "y") == 0)
     {
         axis = YY;
     }
-    else if (strcmp(normal_axis[0], "z") == 0)
+    else if (std::strcmp(normal_axis[0], "z") == 0)
     {
         axis = ZZ;
     }