X-Git-Url: http://biod.pnpi.spb.ru/gitweb/?a=blobdiff_plain;f=src%2Fgromacs%2Fgmxana%2Fgmx_order.cpp;fp=src%2Fgromacs%2Fgmxana%2Fgmx_order.c;h=0b28c3bbd3d644c907a7b12cc490767b0a06776b;hb=c3f2d46e4047f0c465f7234b3784a2fa6f02a065;hp=c0dded8283de35e32d686594ece504779146a3eb;hpb=0595b4a4c763a0bc574658992081abf8b0abc3fe;p=alexxy%2Fgromacs.git diff --git a/src/gromacs/gmxana/gmx_order.c b/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 c0dded8283..0b28c3bbd3 100644 --- a/src/gromacs/gmxana/gmx_order.c +++ b/src/gromacs/gmxana/gmx_order.cpp @@ -36,8 +36,10 @@ */ #include "gmxpre.h" -#include -#include +#include +#include + +#include #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(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((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; }