Merge branch 'master' into pygromacs

[alexxy/gromacs.git] / src / gromacs / gmxana / gmx_cluster.cpp

diff --git a/src/gromacs/gmxana/gmx_cluster.c b/src/gromacs/gmxana/gmx_cluster.cpp
similarity index 96%
rename from src/gromacs/gmxana/gmx_cluster.c
rename to src/gromacs/gmxana/gmx_cluster.cpp
index 011fe580e29347bb9092f6b99c17b051ca5e0bcb..5f88080bb908aae2bf34b78cb298073cad98fad3 100644 (file)
--- a/src/gromacs/gmxana/gmx_cluster.c
+++ b/src/gromacs/gmxana/gmx_cluster.cpp
@@ -36,9 +36,11 @@
  */
 #include "gmxpre.h"
 
-#include <math.h>
-#include <stdlib.h>
-#include <string.h>
+#include <cmath>
+#include <cstdlib>
+#include <cstring>
+
+#include <algorithm>
 
 #include "gromacs/commandline/pargs.h"
 #include "gromacs/fileio/confio.h"
@@ -208,8 +210,8 @@ void mc_optimize(FILE *log, t_mat *m, real *time,
         /* Generate new swapping candidates */
         do
         {
-            iswap = 1+(nn-2)*gmx_rng_uniform_real(rng);
-            jswap = 1+(nn-2)*gmx_rng_uniform_real(rng);
+            iswap = static_cast<int>(1+(nn-2)*gmx_rng_uniform_real(rng));
+            jswap = static_cast<int>(1+(nn-2)*gmx_rng_uniform_real(rng));
         }
         while ((iswap == jswap) || (iswap >= nn-1) || (jswap >= nn-1));
 
@@ -232,7 +234,7 @@ void mc_optimize(FILE *log, t_mat *m, real *time,
         else if (kT > 0)
         {
             /* Try Monte Carlo step */
-            prob = exp(-(enext-ecur)/(enorm*kT));
+            prob = std::exp(-(enext-ecur)/(enorm*kT));
         }
 
         if ((prob == 1) || (gmx_rng_uniform_real(rng) < prob))
@@ -314,7 +316,7 @@ static real rms_dist(int isize, real **d, real **d_r)
     }
     r2 /= (isize*(isize-1))/2;
 
-    return sqrt(r2);
+    return std::sqrt(r2);
 }
 
 static int rms_dist_comp(const void *a, const void *b)
@@ -492,7 +494,7 @@ static void jarvis_patrick(int n1, real **mat, int M, int P,
     t_dist     *row;
     t_clustid  *c;
     int       **nnb;
-    int         i, j, k, cid, diff, max;
+    int         i, j, k, cid, diff, maxval;
     gmx_bool    bChange;
     real      **mcpy = NULL;
 
@@ -531,24 +533,24 @@ static void jarvis_patrick(int n1, real **mat, int M, int P,
         else
         {
             /* Put all neighbors nearer than rmsdcut in the list */
-            max = 0;
-            k   = 0;
+            maxval = 0;
+            k      = 0;
             for (j = 0; (j < n1) && (mat[i][row[j].j] < rmsdcut); j++)
             {
                 if (row[j].j != i)
                 {
-                    if (k >= max)
+                    if (k >= maxval)
                     {
-                        max += 10;
-                        srenew(nnb[i], max);
+                        maxval += 10;
+                        srenew(nnb[i], maxval);
                     }
                     nnb[i][k] = row[j].j;
                     k++;
                 }
             }
-            if (k == max)
+            if (k == maxval)
             {
-                srenew(nnb[i], max+1);
+                srenew(nnb[i], maxval+1);
             }
             nnb[i][k] = -1;
         }
@@ -641,15 +643,6 @@ static void jarvis_patrick(int n1, real **mat, int M, int P,
         }
     }
 
-/* Again, I don't see the point in this... (AF) */
-/*   for(i=0; (i<n1); i++) { */
-/*     for(j=0; (j<n1); j++) */
-/*       mcpy[c[i].conf][c[j].conf] = mat[i][j]; */
-/*   } */
-/*   for(i=0; (i<n1); i++) { */
-/*     for(j=0; (j<n1); j++) */
-/*       mat[i][j] = mcpy[i][j]; */
-/*   } */
     done_matrix(n1, &mcpy);
 
     sfree(c);
@@ -681,7 +674,7 @@ static void gromos(int n1, real **mat, real rmsdcut, t_clusters *clust)
 {
     t_dist *row;
     t_nnb  *nnb;
-    int     i, j, k, j1, max;
+    int     i, j, k, j1, maxval;
 
     /* Put all neighbors nearer than rmsdcut in the list */
     fprintf(stderr, "Making list of neighbors within cutoff ");
@@ -689,17 +682,17 @@ static void gromos(int n1, real **mat, real rmsdcut, t_clusters *clust)
     snew(row, n1);
     for (i = 0; (i < n1); i++)
     {
-        max = 0;
-        k   = 0;
+        maxval = 0;
+        k      = 0;
         /* put all neighbors within cut-off in list */
         for (j = 0; j < n1; j++)
         {
             if (mat[i][j] < rmsdcut)
             {
-                if (k >= max)
+                if (k >= maxval)
                 {
-                    max += 10;
-                    srenew(nnb[i].nb, max);
+                    maxval += 10;
+                    srenew(nnb[i].nb, maxval);
                 }
                 nnb[i].nb[k] = j;
                 k++;
@@ -892,7 +885,7 @@ static int plot_clusters(int nf, real **mat, t_clusters *clust,
 
 static void mark_clusters(int nf, real **mat, real val, t_clusters *clust)
 {
-    int i, j, v;
+    int i, j;
 
     for (i = 0; i < nf; i++)
     {
@@ -917,14 +910,14 @@ static char *parse_filename(const char *fn, int maxnr)
     const char *ext;
     char        buf[STRLEN];
 
-    if (strchr(fn, '%'))
+    if (std::strchr(fn, '%'))
     {
         gmx_fatal(FARGS, "will not number filename %s containing '%c'", fn, '%');
     }
     /* number of digits needed in numbering */
-    i = (int)(log(maxnr)/log(10)) + 1;
+    i = static_cast<int>((std::log(static_cast<real>(maxnr))/std::log(10.0)) + 1);
     /* split fn and ext */
-    ext = strrchr(fn, '.');
+    ext = std::strrchr(fn, '.');
     if (!ext)
     {
         gmx_fatal(FARGS, "cannot separate extension in filename %s", fn);
@@ -932,8 +925,8 @@ static char *parse_filename(const char *fn, int maxnr)
     ext++;
     /* insert e.g. '%03d' between fn and ext */
     sprintf(buf, "%s%%0%dd.%s", fn, i, ext);
-    snew(fnout, strlen(buf)+1);
-    strcpy(fnout, buf);
+    snew(fnout, std::strlen(buf)+1);
+    std::strcpy(fnout, buf);
 
     return fnout;
 }
@@ -966,7 +959,7 @@ static void ana_trans(t_clusters *clust, int nf,
             ntrans[clust->cl[i-1]-1]++;
             ntrans[clust->cl[i]-1]++;
             trans[clust->cl[i-1]-1][clust->cl[i]-1]++;
-            maxtrans = max(maxtrans, trans[clust->cl[i]-1][clust->cl[i-1]-1]);
+            maxtrans = static_cast<int>(std::max(static_cast<real>(maxtrans), trans[clust->cl[i]-1][clust->cl[i-1]-1]));
         }
     }
     ffprintf_dd(stderr, log, buf, "Counted %d transitions in total, "
@@ -974,7 +967,7 @@ static void ana_trans(t_clusters *clust, int nf,
     if (transfn)
     {
         fp = gmx_ffopen(transfn, "w");
-        i  = min(maxtrans+1, 80);
+        i  = std::min(maxtrans+1, 80);
         write_xpm(fp, 0, "Cluster Transitions", "# transitions",
                   "from cluster", "to cluster",
                   clust->ncl, clust->ncl, axis, axis, trans,
@@ -1031,7 +1024,7 @@ static void analyze_clusters(int nf, t_clusters *clust, real **rmsd,
         /* do we write all structures? */
         if (write_ncl)
         {
-            trxsfn = parse_filename(trxfn, max(write_ncl, clust->ncl));
+            trxsfn = parse_filename(trxfn, std::max(write_ncl, clust->ncl));
             snew(bWrite, nf);
         }
         ffprintf_ss(stderr, log, buf, "Writing %s structure for each cluster to %s\n",
@@ -1281,7 +1274,6 @@ static void analyze_clusters(int nf, t_clusters *clust, real **rmsd,
             {
                 do_fit(natom, mass, xtps, xav);
             }
-            r = cl;
             write_trx(trxout, iosize, outidx, atoms, cl, time[midstr], zerobox, xav, NULL, NULL);
         }
     }
@@ -1325,10 +1317,10 @@ static void convert_mat(t_matrix *mat, t_mat *rms)
         for (j = i; j < mat->nx; j++)
         {
             rms->sumrms += rms->mat[i][j];
-            rms->maxrms  = max(rms->maxrms, rms->mat[i][j]);
+            rms->maxrms  = std::max(rms->maxrms, rms->mat[i][j]);
             if (j != i)
             {
-                rms->minrms = min(rms->minrms, rms->mat[i][j]);
+                rms->minrms = std::min(rms->minrms, rms->mat[i][j]);
             }
         }
     }
@@ -1429,7 +1421,7 @@ int gmx_cluster(int argc, char *argv[])
     real              *eigenvectors;
 
     int                isize = 0, ifsize = 0, iosize = 0;
-    atom_id           *index = NULL, *fitidx, *outidx;
+    atom_id           *index = NULL, *fitidx = NULL, *outidx = NULL;
     char              *grpname;
     real               rmsd, **d1, **d2, *time = NULL, time_invfac, *mass = NULL;
     char               buf[STRLEN], buf1[80], title[STRLEN];
@@ -1749,7 +1741,7 @@ int gmx_cluster(int argc, char *argv[])
     else   /* !bReadMat */
     {
         rms  = init_mat(nf, method == m_diagonalize);
-        nrms = ((gmx_int64_t)nf*((gmx_int64_t)nf-1))/2;
+        nrms = (static_cast<gmx_int64_t>(nf)*static_cast<gmx_int64_t>(nf-1))/2;
         if (!bRMSdist)
         {
             fprintf(stderr, "Computing %dx%d RMS deviation matrix\n", nf, nf);
@@ -1770,8 +1762,8 @@ int gmx_cluster(int argc, char *argv[])
                     rmsd = rmsdev(isize, mass, xx[i2], x1);
                     set_mat_entry(rms, i1, i2, rmsd);
                 }
-                nrms -= (gmx_int64_t) (nf-i1-1);
-                fprintf(stderr, "\r# RMSD calculations left: " "%"GMX_PRId64 "   ", nrms);
+                nrms -= nf-i1-1;
+                fprintf(stderr, "\r# RMSD calculations left: " "%" GMX_PRId64 "   ", nrms);
             }
             sfree(x1);
         }
@@ -1795,8 +1787,8 @@ int gmx_cluster(int argc, char *argv[])
                     calc_dist(isize, xx[i2], d2);
                     set_mat_entry(rms, i1, i2, rms_dist(isize, d1, d2));
                 }
-                nrms -= (nf-i1-1);
-                fprintf(stderr, "\r# RMSD calculations left: " "%"GMX_PRId64 "   ", nrms);
+                nrms -= nf-i1-1;
+                fprintf(stderr, "\r# RMSD calculations left: " "%" GMX_PRId64 "   ", nrms);
             }
             /* Clean up work arrays */
             for (i = 0; (i < isize); i++)
@@ -1862,7 +1854,7 @@ int gmx_cluster(int argc, char *argv[])
             /* Do a diagonalization */
             snew(eigenvalues, nf);
             snew(eigenvectors, nf*nf);
-            memcpy(eigenvectors, rms->mat[0], nf*nf*sizeof(real));
+            std::memcpy(eigenvectors, rms->mat[0], nf*nf*sizeof(real));
             eigensolver(eigenvectors, nf, 0, nf, eigenvalues, rms->mat[0]);
             sfree(eigenvectors);
 
@@ -1914,7 +1906,7 @@ int gmx_cluster(int argc, char *argv[])
         {
             useatoms.atomname[i]    = top.atoms.atomname[index[i]];
             useatoms.atom[i].resind = top.atoms.atom[index[i]].resind;
-            useatoms.nres           = max(useatoms.nres, useatoms.atom[i].resind+1);
+            useatoms.nres           = std::max(useatoms.nres, useatoms.atom[i].resind+1);
             copy_rvec(xtps[index[i]], usextps[i]);
         }
         useatoms.nr = isize;
@@ -1962,7 +1954,7 @@ int gmx_cluster(int argc, char *argv[])
         {
             write_xpm_split(fp, 0, title, "RMSD (nm)", buf, buf,
                             nf, nf, time, time, rms->mat, 0.0, rms->maxrms, &nlevels,
-                            rlo_top, rhi_top, 0.0, (real) ncluster,
+                            rlo_top, rhi_top, 0.0, ncluster,
                             &ncluster, TRUE, rlo_bot, rhi_bot);
         }
         else