*/
#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"
/****************************************************************************/
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];
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;
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;
sgmol[i] += cost2;
/* determine distribution */
- ibin = nsgbin * cost2;
+ ibin = static_cast<int>(nsgbin * cost2);
if (ibin < nsgbin)
{
sgbin[ibin]++;
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;
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;
/* 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)
{
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
}
else
{
- sdbangle += acos(Sz[axis]);
+ sdbangle += std::acos(Sz[axis]);
}
}
else
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 */
}
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
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)
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);
/* 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)
{
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 */
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++)
{
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]));
}
}
/*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;
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;
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*/
}
}
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;
}