biod.pnpi.spb.ru
/
alexxy
/
gromacs.git
/ blobdiff
commit
grep
author
committer
pickaxe
?
search:
re
summary
|
shortlog
|
log
|
commit
|
commitdiff
|
tree
raw
|
inline
| side by side
Merge branch 'master' into pygromacs
[alexxy/gromacs.git]
/
src
/
gromacs
/
gmxana
/
gmx_order.cpp
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 c0dded8283de35e32d686594ece504779146a3eb..0b28c3bbd3d644c907a7b12cc490767b0a06776b 100644
(file)
--- a/
src/gromacs/gmxana/gmx_order.c
+++ b/
src/gromacs/gmxana/gmx_order.cpp
@@
-36,8
+36,10
@@
*/
#include "gmxpre.h"
*/
#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/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/cstringutil.h"
#include "gromacs/utility/fatalerror.h"
#include "gromacs/utility/futil.h"
+#include "gromacs/utility/gmxassert.h"
#include "gromacs/utility/smalloc.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)
{
real sgslice[], real skslice[],
gmx_rmpbc_t gpbc)
{
- FILE *fpoutdist;
- char fnsgdist[32];
int ix, jx, nsgbin, *sgbin;
int ix, jx, nsgbin, *sgbin;
- int i
1, 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;
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];
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);
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;
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++)
{
rmean = 0;
for (j = 0; (j < 4); j++)
{
- r_nn[j][i] = sqrt(r_nn[j][i]);
+ r_nn[j][i] = s
td::s
qrt(r_nn[j][i]);
rmean += r_nn[j][i];
}
rmean /= 4;
rmean += r_nn[j][i];
}
rmean /= 4;
@@
-192,7
+191,7
@@
static void find_nearest_neighbours(int ePBC,
sgmol[i] += cost2;
/* determine distribution */
sgmol[i] += cost2;
/* determine distribution */
- ibin =
nsgbin * cost2
;
+ ibin =
static_cast<int>(nsgbin * cost2)
;
if (ibin < nsgbin)
{
sgbin[ibin]++;
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;
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;
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 */
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 */
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;
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 */
/* 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)
{
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]);
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
/* 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
{
}
else
{
- sdbangle += acos(Sz[axis]);
+ sdbangle +=
std::
acos(Sz[axis]);
}
}
else
}
}
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]);
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 */
}
/* 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);
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);
/* 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
/* 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++)
{
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)
}
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];
}
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);
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: */
/* 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)
{
}
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;
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));
}
//fprintf(stderr, "Min dist %f; trace %f\n", tmpdist, trace(box));
- (*distvals)[j][i] += sqrt(tmpdist);
+ (*distvals)[j][i] += s
td::s
qrt(tmpdist);
}
}
} /* end loop j, over all atoms in group */
}
}
} /* 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++)
{
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++)
}
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(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;
/*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;
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;
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;
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]];
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*/
}
}
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 */
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;
}
{
axis = XX;
}
- else if (strcmp(normal_axis[0], "y") == 0)
+ else if (st
d::st
rcmp(normal_axis[0], "y") == 0)
{
axis = YY;
}
{
axis = YY;
}
- else if (strcmp(normal_axis[0], "z") == 0)
+ else if (st
d::st
rcmp(normal_axis[0], "z") == 0)
{
axis = ZZ;
}
{
axis = ZZ;
}