Merge branch 'master' into pygromacs
[alexxy/gromacs.git] / src / gromacs / gmxana / gmx_editconf.cpp
similarity index 93%
rename from src/gromacs/gmxana/gmx_editconf.c
rename to src/gromacs/gmxana/gmx_editconf.cpp
index 6013b75f8e50020b3f860bea8449af440f4cb31e..c32cfb6b928c678d8f190e4848cb312d6ff4c4f5 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"
 #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"
 
-typedef struct
-{
-    char sanm[12];
-    int  natm;
-    int  nw;
-    char anm[6][12];
-} t_simat;
-
-typedef struct
-{
-    char    reso[12];
-    char    resn[12];
-    int     nsatm;
-    t_simat sat[3];
-} t_simlist;
 
 real calc_mass(t_atoms *atoms, gmx_bool bGetMass, gmx_atomprop_t aps)
 {
@@ -99,19 +87,18 @@ real calc_mass(t_atoms *atoms, gmx_bool bGetMass, gmx_atomprop_t aps)
     return tmass;
 }
 
-real calc_geom(int isize, atom_id *index, rvec *x, rvec geom_center, rvec min,
-               rvec max, gmx_bool bDiam)
+real calc_geom(int isize, atom_id *index, rvec *x, rvec geom_center, rvec minval,
+               rvec maxval, gmx_bool bDiam)
 {
     real  diam2, d;
-    char *grpnames;
     int   ii, i, j;
 
     clear_rvec(geom_center);
     diam2 = 0;
     if (isize == 0)
     {
-        clear_rvec(min);
-        clear_rvec(max);
+        clear_rvec(minval);
+        clear_rvec(maxval);
     }
     else
     {
@@ -125,7 +112,7 @@ real calc_geom(int isize, atom_id *index, rvec *x, rvec geom_center, rvec min,
         }
         for (j = 0; j < DIM; j++)
         {
-            min[j] = max[j] = x[ii][j];
+            minval[j] = maxval[j] = x[ii][j];
         }
         for (i = 0; i < isize; i++)
         {
@@ -140,13 +127,13 @@ real calc_geom(int isize, atom_id *index, rvec *x, rvec geom_center, rvec min,
             rvec_inc(geom_center, x[ii]);
             for (j = 0; j < DIM; j++)
             {
-                if (x[ii][j] < min[j])
+                if (x[ii][j] < minval[j])
                 {
-                    min[j] = x[ii][j];
+                    minval[j] = x[ii][j];
                 }
-                if (x[ii][j] > max[j])
+                if (x[ii][j] > maxval[j])
                 {
-                    max[j] = x[ii][j];
+                    maxval[j] = x[ii][j];
                 }
             }
             if (bDiam)
@@ -156,7 +143,7 @@ real calc_geom(int isize, atom_id *index, rvec *x, rvec geom_center, rvec min,
                     for (j = i + 1; j < isize; j++)
                     {
                         d     = distance2(x[ii], x[index[j]]);
-                        diam2 = max(d, diam2);
+                        diam2 = std::max(d, diam2);
                     }
                 }
                 else
@@ -164,7 +151,7 @@ real calc_geom(int isize, atom_id *index, rvec *x, rvec geom_center, rvec min,
                     for (j = i + 1; j < isize; j++)
                     {
                         d     = distance2(x[i], x[j]);
-                        diam2 = max(d, diam2);
+                        diam2 = std::max(d, diam2);
                     }
                 }
             }
@@ -172,7 +159,7 @@ real calc_geom(int isize, atom_id *index, rvec *x, rvec geom_center, rvec min,
         svmul(1.0 / isize, geom_center, geom_center);
     }
 
-    return sqrt(diam2);
+    return std::sqrt(diam2);
 }
 
 void center_conf(int natom, rvec *x, rvec center, rvec geom_cent)
@@ -232,7 +219,6 @@ void read_bfac(const char *fn, int *n_bfac, double **bfac_val, int **bfac_nr)
 void set_pdb_conf_bfac(int natoms, int nres, t_atoms *atoms, int n_bfac,
                        double *bfac, int *bfac_nr, gmx_bool peratom)
 {
-    FILE    *out;
     real     bfac_min, bfac_max;
     int      i, n;
     gmx_bool found;
@@ -269,7 +255,7 @@ void set_pdb_conf_bfac(int natoms, int nres, t_atoms *atoms, int n_bfac,
         bfac_max /= 10;
         bfac_min /= 10;
     }
-    while ((fabs(bfac_max) < 0.5) && (fabs(bfac_min) < 0.5))
+    while ((std::abs(bfac_max) < 0.5) && (std::abs(bfac_min) < 0.5))
     {
         fprintf(stderr,
                 "Range of values for B-factors too small (min %g, max %g) "
@@ -332,11 +318,11 @@ void pdb_legend(FILE *out, int natoms, int nres, t_atoms *atoms, rvec x[])
     zmin     = 1e10;
     for (i = 0; (i < natoms); i++)
     {
-        xmin     = min(xmin, x[i][XX]);
-        ymin     = min(ymin, x[i][YY]);
-        zmin     = min(zmin, x[i][ZZ]);
-        bfac_min = min(bfac_min, atoms->pdbinfo[i].bfac);
-        bfac_max = max(bfac_max, atoms->pdbinfo[i].bfac);
+        xmin     = std::min(xmin, x[i][XX]);
+        ymin     = std::min(ymin, x[i][YY]);
+        zmin     = std::min(zmin, x[i][ZZ]);
+        bfac_min = std::min(bfac_min, atoms->pdbinfo[i].bfac);
+        bfac_max = std::max(bfac_max, atoms->pdbinfo[i].bfac);
     }
     fprintf(stderr, "B-factors range from %g to %g\n", bfac_min, bfac_max);
     for (i = 1; (i < 12); i++)
@@ -394,9 +380,9 @@ void visualize_box(FILE *out, int a0, int r0, matrix box, rvec gridsize)
     a0++;
     r0++;
 
-    nx   = (int) (gridsize[XX] + 0.5);
-    ny   = (int) (gridsize[YY] + 0.5);
-    nz   = (int) (gridsize[ZZ] + 0.5);
+    nx   = static_cast<int>(gridsize[XX] + 0.5);
+    ny   = static_cast<int>(gridsize[YY] + 0.5);
+    nz   = static_cast<int>(gridsize[ZZ] + 0.5);
     nbox = nx * ny * nz;
     if (TRICLINIC(box))
     {
@@ -459,8 +445,7 @@ void visualize_box(FILE *out, int a0, int r0, matrix box, rvec gridsize)
         }
         for (i = 0; i < 24; i += 2)
         {
-            fprintf(out, "CONECT%5d%5d\n", a0 + rectedge[i], a0 + rectedge[i
-                                                                           + 1]);
+            fprintf(out, "CONECT%5d%5d\n", a0 + rectedge[i], a0 + rectedge[i + 1]);
         }
     }
 }
@@ -471,7 +456,7 @@ void calc_rotmatrix(rvec principal_axis, rvec targetvec, matrix rotmatrix)
     real ux, uy, uz, costheta, sintheta;
 
     costheta = cos_angle(principal_axis, targetvec);
-    sintheta = sqrt(1.0-costheta*costheta); /* sign is always positive since 0<theta<pi */
+    sintheta = std::sqrt(1.0-costheta*costheta); /* sign is always positive since 0<theta<pi */
 
     /* Determine rotation from cross product with target vector */
     cprod(principal_axis, targetvec, rotvec);
@@ -614,7 +599,7 @@ int gmx_editconf(int argc, char *argv[])
         "For complex molecules, the periodicity removal routine may break down, "
         "in that case you can use [gmx-trjconv]."
     };
-    static real     dist    = 0.0, rbox = 0.0, to_diam = 0.0;
+    static real     dist    = 0.0;
     static gmx_bool bNDEF   = FALSE, bRMPBC = FALSE, bCenter = FALSE, bReadVDW =
         FALSE, bCONECT      = FALSE;
     static gmx_bool peratom = FALSE, bLegend = FALSE, bOrient = FALSE, bMead =
@@ -720,15 +705,15 @@ int gmx_editconf(int argc, char *argv[])
     t_topology    *top     = NULL;
     t_atoms        atoms;
     char          *grpname, *sgrpname, *agrpname;
-    int            isize, ssize, tsize, asize;
-    atom_id       *index, *sindex, *tindex, *aindex;
-    rvec          *x, *v, gc, min, max, size;
+    int            isize, ssize, asize;
+    atom_id       *index, *sindex, *aindex;
+    rvec          *x, *v, gc, rmin, rmax, size;
     int            ePBC;
     matrix         box, rotmatrix, trans;
     rvec           princd, tmpvec;
-    gmx_bool       bIndex, bSetSize, bSetAng, bCubic, bDist, bSetCenter, bAlign;
+    gmx_bool       bIndex, bSetSize, bSetAng, bDist, bSetCenter, bAlign;
     gmx_bool       bHaveV, bScale, bRho, bTranslate, bRotate, bCalcGeom, bCalcDiam;
-    real           xs, ys, zs, xcent, ycent, zcent, diam = 0, mass = 0, d, vdw;
+    real           diam = 0, mass = 0, d, vdw;
     gmx_atomprop_t aps;
     gmx_conect     conect;
     output_env_t   oenv;
@@ -770,7 +755,10 @@ int gmx_editconf(int argc, char *argv[])
     }
     bScale    = bScale || bRho;
     bCalcGeom = bCenter || bRotate || bOrient || bScale;
-    bCalcDiam = btype[0][0] == 'c' || btype[0][0] == 'd' || btype[0][0] == 'o';
+
+    GMX_RELEASE_ASSERT(btype[0] != NULL, "Option setting inconsistency; btype[0] is NULL");
+
+    bCalcDiam = (btype[0][0] == 'c' || btype[0][0] == 'd' || btype[0][0] == 'o');
 
     infile = ftp2fn(efSTX, NFILE, fnm);
     if (bMead)
@@ -823,7 +811,7 @@ int gmx_editconf(int argc, char *argv[])
     {
         real vol = det(box);
         printf("Volume: %g nm^3, corresponds to roughly %d electrons\n",
-               vol, 100*((int)(vol*4.5)));
+               vol, 100*(static_cast<int>(vol*4.5)));
     }
 
     if (bMead || bGrasp || bCONECT)
@@ -867,7 +855,7 @@ int gmx_editconf(int argc, char *argv[])
                     {
                         sig6 = c12/c6;
                     }
-                    vdw   = 0.5*pow(sig6, 1.0/6.0);
+                    vdw   = 0.5*std::pow(sig6, static_cast<real>(1.0/6.0));
                 }
                 else
                 {
@@ -934,8 +922,8 @@ int gmx_editconf(int argc, char *argv[])
             ssize  = atoms.nr;
             sindex = NULL;
         }
-        diam = calc_geom(ssize, sindex, x, gc, min, max, bCalcDiam);
-        rvec_sub(max, min, size);
+        diam = calc_geom(ssize, sindex, x, gc, rmin, rmax, bCalcDiam);
+        rvec_sub(rmax, rmin, size);
         printf("    system size :%7.3f%7.3f%7.3f (nm)\n",
                size[XX], size[YY], size[ZZ]);
         if (bCalcDiam)
@@ -947,11 +935,11 @@ int gmx_editconf(int argc, char *argv[])
                norm(box[XX]), norm(box[YY]), norm(box[ZZ]));
         printf("    box angles  :%7.2f%7.2f%7.2f (degrees)\n",
                norm2(box[ZZ]) == 0 ? 0 :
-               RAD2DEG*acos(cos_angle_no_table(box[YY], box[ZZ])),
+               RAD2DEG*std::acos(cos_angle_no_table(box[YY], box[ZZ])),
                norm2(box[ZZ]) == 0 ? 0 :
-               RAD2DEG*acos(cos_angle_no_table(box[XX], box[ZZ])),
+               RAD2DEG*std::acos(cos_angle_no_table(box[XX], box[ZZ])),
                norm2(box[YY]) == 0 ? 0 :
-               RAD2DEG*acos(cos_angle_no_table(box[XX], box[YY])));
+               RAD2DEG*std::acos(cos_angle_no_table(box[XX], box[YY])));
         printf("    box volume  :%7.2f               (nm^3)\n", det(box));
     }
 
@@ -994,7 +982,7 @@ int gmx_editconf(int argc, char *argv[])
                           "zero mass (%g) or volume (%g)\n", mass, vol);
             }
 
-            scale[XX] = scale[YY] = scale[ZZ] = pow(dens/rho, 1.0/3.0);
+            scale[XX] = scale[YY] = scale[ZZ] = std::pow(dens/rho, static_cast<real>(1.0/3.0));
             fprintf(stderr, "Scaling all box vectors by %g\n", scale[XX]);
         }
         scale_conf(atoms.nr, x, box, scale);
@@ -1096,8 +1084,8 @@ int gmx_editconf(int argc, char *argv[])
     if (bCalcGeom)
     {
         /* recalc geometrical center and max and min coordinates and size */
-        calc_geom(ssize, sindex, x, gc, min, max, FALSE);
-        rvec_sub(max, min, size);
+        calc_geom(ssize, sindex, x, gc, rmin, rmax, FALSE);
+        rvec_sub(rmax, rmin, size);
         if (bScale || bOrient || bRotate)
         {
             printf("new system size : %6.3f %6.3f %6.3f\n",
@@ -1105,7 +1093,7 @@ int gmx_editconf(int argc, char *argv[])
         }
     }
 
-    if (bSetSize || bDist || (btype[0][0] == 't' && bSetAng))
+    if ((btype[0] != NULL) && (bSetSize || bDist || (btype[0][0] == 't' && bSetAng)))
     {
         ePBC = epbcXYZ;
         if (!(bSetSize || bDist))
@@ -1162,16 +1150,16 @@ int gmx_editconf(int argc, char *argv[])
                     box[YY][YY] = d;
                     box[ZZ][XX] = d/2;
                     box[ZZ][YY] = d/2;
-                    box[ZZ][ZZ] = d*sqrt(2)/2;
+                    box[ZZ][ZZ] = d*std::sqrt(2.0)/2.0;
                 }
                 else
                 {
                     box[XX][XX] = d;
                     box[YY][XX] = d/3;
-                    box[YY][YY] = d*sqrt(2)*2/3;
+                    box[YY][YY] = d*std::sqrt(2.0)*2.0/3.0;
                     box[ZZ][XX] = -d/3;
-                    box[ZZ][YY] = d*sqrt(2)/3;
-                    box[ZZ][ZZ] = d*sqrt(6)/3;
+                    box[ZZ][YY] = d*std::sqrt(2.0)/3.0;
+                    box[ZZ][ZZ] = d*std::sqrt(6.0)/3.0;
                 }
                 break;
         }
@@ -1192,7 +1180,7 @@ int gmx_editconf(int argc, char *argv[])
     /* print some */
     if (bCalcGeom)
     {
-        calc_geom(ssize, sindex, x, gc, min, max, FALSE);
+        calc_geom(ssize, sindex, x, gc, rmin, rmax, FALSE);
         printf("new center      :%7.3f%7.3f%7.3f (nm)\n", gc[XX], gc[YY], gc[ZZ]);
     }
     if (bOrient || bScale || bDist || bSetSize)
@@ -1201,11 +1189,11 @@ int gmx_editconf(int argc, char *argv[])
                norm(box[XX]), norm(box[YY]), norm(box[ZZ]));
         printf("new box angles  :%7.2f%7.2f%7.2f (degrees)\n",
                norm2(box[ZZ]) == 0 ? 0 :
-               RAD2DEG*acos(cos_angle_no_table(box[YY], box[ZZ])),
+               RAD2DEG*std::acos(cos_angle_no_table(box[YY], box[ZZ])),
                norm2(box[ZZ]) == 0 ? 0 :
-               RAD2DEG*acos(cos_angle_no_table(box[XX], box[ZZ])),
+               RAD2DEG*std::acos(cos_angle_no_table(box[XX], box[ZZ])),
                norm2(box[YY]) == 0 ? 0 :
-               RAD2DEG*acos(cos_angle_no_table(box[XX], box[YY])));
+               RAD2DEG*std::acos(cos_angle_no_table(box[XX], box[YY])));
         printf("new box volume  :%7.2f               (nm^3)\n", det(box));
     }