*
* Copyright (c) 1991-2000, University of Groningen, The Netherlands.
* Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2013,2014,2015,2016,2017,2018, by the GROMACS development team, led by
+ * Copyright (c) 2013,2014,2015,2016,2017,2018,2019, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
/* directions. */
/****************************************************************************/
-static void calc_h2order(const char *fn, const int index[], int ngx, rvec **slDipole,
- real **slOrder, real *slWidth, int *nslices,
- const t_topology *top, int ePBC,
- int axis, gmx_bool bMicel, int micel[], int nmic,
- const gmx_output_env_t *oenv)
+static void calc_h2order(const char* fn,
+ const int index[],
+ int ngx,
+ rvec** slDipole,
+ real** slOrder,
+ real* slWidth,
+ int* nslices,
+ const t_topology* top,
+ int ePBC,
+ int axis,
+ gmx_bool bMicel,
+ int micel[],
+ int nmic,
+ const gmx_output_env_t* oenv)
{
rvec *x0, /* coordinates with pbc */
- dipole, /* dipole moment due to one molecules */
- normal,
- com; /* center of mass of micel, with bMicel */
- rvec *dip; /* sum of dipoles, unnormalized */
+ dipole, /* dipole moment due to one molecules */
+ normal, com; /* center of mass of micel, with bMicel */
+ rvec* dip; /* sum of dipoles, unnormalized */
matrix box; /* box (3x3) */
- t_trxstatus *status;
- real t, /* time from trajectory */
- *sum, /* sum of all cosines of dipoles, per slice */
- *frame; /* order over one frame */
- int natoms, /* nr. atoms in trj */
- i, j, teller = 0,
- slice = 0, /* current slice number */
- *count; /* nr. of atoms in one slice */
- gmx_rmpbc_t gpbc = nullptr;
+ t_trxstatus* status;
+ real t, /* time from trajectory */
+ *sum, /* sum of all cosines of dipoles, per slice */
+ *frame; /* order over one frame */
+ int natoms, /* nr. atoms in trj */
+ i, j, teller = 0, slice = 0, /* current slice number */
+ *count; /* nr. of atoms in one slice */
+ gmx_rmpbc_t gpbc = nullptr;
if ((natoms = read_first_x(oenv, &status, fn, &t, &x0, box)) == 0)
{
if (!*nslices)
{
*nslices = static_cast<int>(box[axis][axis] * 10); /* default value */
-
-
}
switch (axis)
{
case 0:
- normal[0] = 1; normal[1] = 0; normal[2] = 0;
+ normal[0] = 1;
+ normal[1] = 0;
+ normal[2] = 0;
break;
case 1:
- normal[0] = 0; normal[1] = 1; normal[2] = 0;
+ normal[0] = 0;
+ normal[1] = 1;
+ normal[2] = 0;
break;
case 2:
- normal[0] = 0; normal[1] = 0; normal[2] = 1;
+ normal[0] = 0;
+ normal[1] = 0;
+ normal[2] = 1;
break;
- default:
- gmx_fatal(FARGS, "No valid value for -axis-. Exiting.\n");
+ default: gmx_fatal(FARGS, "No valid value for -axis-. Exiting.\n");
}
clear_rvec(dipole);
snew(dip, *nslices);
snew(frame, *nslices);
- *slWidth = box[axis][axis]/(*nslices);
- fprintf(stderr, "Box divided in %d slices. Initial width of slice: %f\n",
- *nslices, *slWidth);
+ *slWidth = box[axis][axis] / (*nslices);
+ fprintf(stderr, "Box divided in %d slices. Initial width of slice: %f\n", *nslices, *slWidth);
teller = 0;
/*********** Start processing trajectory ***********/
do
{
- *slWidth = box[axis][axis]/(*nslices);
+ *slWidth = box[axis][axis] / (*nslices);
teller++;
gmx_rmpbc(gpbc, natoms, box, x0);
calc_xcm(x0, nmic, micel, top->atoms.atom, com, FALSE);
}
- for (i = 0; i < ngx/3; i++)
+ for (i = 0; i < ngx / 3; i++)
{
/* put all waters in box */
for (j = 0; j < DIM; j++)
{
- if (x0[index[3*i]][j] < 0)
+ if (x0[index[3 * i]][j] < 0)
{
- x0[index[3*i]][j] += box[j][j];
- x0[index[3*i+1]][j] += box[j][j];
- x0[index[3*i+2]][j] += box[j][j];
+ x0[index[3 * i]][j] += box[j][j];
+ x0[index[3 * i + 1]][j] += box[j][j];
+ x0[index[3 * i + 2]][j] += box[j][j];
}
- if (x0[index[3*i]][j] > box[j][j])
+ if (x0[index[3 * i]][j] > box[j][j])
{
- x0[index[3*i]][j] -= box[j][j];
- x0[index[3*i+1]][j] -= box[j][j];
- x0[index[3*i+2]][j] -= box[j][j];
+ x0[index[3 * i]][j] -= box[j][j];
+ x0[index[3 * i + 1]][j] -= box[j][j];
+ x0[index[3 * i + 2]][j] -= box[j][j];
}
}
for (j = 0; j < DIM; j++)
{
- dipole[j] =
- x0[index[3*i]][j] * top->atoms.atom[index[3*i]].q +
- x0[index[3*i+1]][j] * top->atoms.atom[index[3*i+1]].q +
- x0[index[3*i+2]][j] * top->atoms.atom[index[3*i+2]].q;
+ dipole[j] = x0[index[3 * i]][j] * top->atoms.atom[index[3 * i]].q
+ + x0[index[3 * i + 1]][j] * top->atoms.atom[index[3 * i + 1]].q
+ + x0[index[3 * i + 2]][j] * top->atoms.atom[index[3 * i + 2]].q;
}
/* now we have a dipole vector. Might as well safe it. Then the
*/
if (bMicel)
- { /* this is for spherical interfaces */
- rvec_sub(com, x0[index[3*i]], normal); /* vector from Oxygen to COM */
- slice = static_cast<int>(norm(normal)/(*slWidth)); /* spherical slice */
+ { /* this is for spherical interfaces */
+ rvec_sub(com, x0[index[3 * i]], normal); /* vector from Oxygen to COM */
+ slice = static_cast<int>(norm(normal) / (*slWidth)); /* spherical slice */
- sum[slice] += iprod(dipole, normal) / (norm(dipole) * norm(normal));
+ sum[slice] += iprod(dipole, normal) / (norm(dipole) * norm(normal));
frame[slice] += iprod(dipole, normal) / (norm(dipole) * norm(normal));
count[slice]++;
-
}
else
{
/* this is for flat interfaces */
/* determine which slice atom is in */
- slice = static_cast<int>(x0[index[3*i]][axis] / (*slWidth));
+ slice = static_cast<int>(x0[index[3 * i]][axis] / (*slWidth));
if (slice < 0 || slice >= *nslices)
{
- fprintf(stderr, "Coordinate: %f ", x0[index[3*i]][axis]);
+ fprintf(stderr, "Coordinate: %f ", x0[index[3 * i]][axis]);
fprintf(stderr, "HELP PANIC! slice = %d, OUT OF RANGE!\n", slice);
}
else
{
rvec_add(dipole, dip[slice], dip[slice]);
/* Add dipole to total. mag[slice] is total dipole in axis direction */
- sum[slice] += iprod(dipole, normal)/norm(dipole);
- frame[slice] += iprod(dipole, normal)/norm(dipole);
+ sum[slice] += iprod(dipole, normal) / norm(dipole);
+ frame[slice] += iprod(dipole, normal) / norm(dipole);
/* increase count for that slice */
count[slice]++;
}
}
}
- }
- while (read_next_x(oenv, status, &t, x0, box));
+ } while (read_next_x(oenv, status, &t, x0, box));
/*********** done with status file **********/
fprintf(stderr, "\nRead trajectory. Printing parameters to file\n");
*slOrder = sum; /* copy a pointer, I hope */
*slDipole = dip;
- sfree(x0); /* free memory used by coordinate arrays */
+ sfree(x0); /* free memory used by coordinate arrays */
}
-static void h2order_plot(rvec dipole[], real order[], const char *afile,
- int nslices, real slWidth, const gmx_output_env_t *oenv)
+static void h2order_plot(rvec dipole[], real order[], const char* afile, int nslices, real slWidth, const gmx_output_env_t* oenv)
{
- FILE *ord; /* xvgr files with order parameters */
- int slice; /* loop index */
- char buf[256]; /* for xvgr title */
- real factor; /* conversion to Debye from electron*nm */
+ FILE* ord; /* xvgr files with order parameters */
+ int slice; /* loop index */
+ char buf[256]; /* for xvgr title */
+ real factor; /* conversion to Debye from electron*nm */
/* factor = 1e-9*1.60217733e-19/3.336e-30 */
- factor = 1.60217733/3.336e-2;
+ factor = 1.60217733 / 3.336e-2;
fprintf(stderr, "%d slices\n", nslices);
sprintf(buf, "Water orientation with respect to normal");
- ord = xvgropen(afile, buf,
- "box (nm)", "mu_x, mu_y, mu_z (D), cosine with normal", oenv);
+ ord = xvgropen(afile, buf, "box (nm)", "mu_x, mu_y, mu_z (D), cosine with normal", oenv);
for (slice = 0; slice < nslices; slice++)
{
- fprintf(ord, "%8.3f %8.3f %8.3f %8.3f %e\n", slWidth*slice,
- factor*dipole[slice][XX], factor*dipole[slice][YY],
- factor*dipole[slice][ZZ], order[slice]);
+ fprintf(ord, "%8.3f %8.3f %8.3f %8.3f %e\n", slWidth * slice, factor * dipole[slice][XX],
+ factor * dipole[slice][YY], factor * dipole[slice][ZZ], order[slice]);
}
xvgrclose(ord);
}
-int gmx_h2order(int argc, char *argv[])
+int gmx_h2order(int argc, char* argv[])
{
- const char *desc[] = {
+ const char* desc[] = {
"[THISMODULE] computes the orientation of water molecules with respect to the normal",
"of the box. The program determines the average cosine of the angle",
"between the dipole moment of water and an axis of the box. The box is",
"dipole and the axis from the center of mass to the oxygen is calculated",
"instead of the angle between the dipole and a box axis."
};
- static int axis = 2; /* normal to memb. default z */
- static const char *axtitle = "Z";
- static int nslices = 0; /* nr of slices defined */
- t_pargs pa[] = {
- { "-d", FALSE, etSTR, {&axtitle},
- "Take the normal on the membrane in direction X, Y or Z." },
- { "-sl", FALSE, etINT, {&nslices},
- "Calculate order parameter as function of boxlength, dividing the box"
- " in this number of slices."}
- };
- const char *bugs[] = {
+ static int axis = 2; /* normal to memb. default z */
+ static const char* axtitle = "Z";
+ static int nslices = 0; /* nr of slices defined */
+ t_pargs pa[] = { { "-d",
+ FALSE,
+ etSTR,
+ { &axtitle },
+ "Take the normal on the membrane in direction X, Y or Z." },
+ { "-sl",
+ FALSE,
+ etINT,
+ { &nslices },
+ "Calculate order parameter as function of boxlength, dividing the box"
+ " in this number of slices." } };
+ const char* bugs[] = {
"The program assigns whole water molecules to a slice, based on the first "
"atom of three in the index file group. It assumes an order O,H,H. "
"Name is not important, but the order is. If this demand is not met, "
"assigning molecules to slices is different."
};
- gmx_output_env_t *oenv;
- real *slOrder, /* av. cosine, per slice */
- slWidth = 0.0; /* width of a slice */
- rvec *slDipole;
- char *grpname, /* groupnames */
- *micname;
- int ngx, /* nr. of atomsin sol group */
- nmic = 0; /* nr. of atoms in micelle */
- t_topology *top; /* topology */
- int ePBC;
- int *index, /* indices for solvent group */
- *micelle = nullptr;
- gmx_bool bMicel = FALSE; /* think we're a micel */
- t_filenm fnm[] = { /* files for g_order */
- { efTRX, "-f", nullptr, ffREAD }, /* trajectory file */
- { efNDX, nullptr, nullptr, ffREAD }, /* index file */
- { efNDX, "-nm", nullptr, ffOPTRD }, /* index with micelle atoms */
- { efTPR, nullptr, nullptr, ffREAD }, /* topology file */
- { efXVG, "-o", "order", ffWRITE }, /* xvgr output file */
+ gmx_output_env_t* oenv;
+ real * slOrder, /* av. cosine, per slice */
+ slWidth = 0.0; /* width of a slice */
+ rvec* slDipole;
+ char *grpname, /* groupnames */
+ *micname;
+ int ngx, /* nr. of atomsin sol group */
+ nmic = 0; /* nr. of atoms in micelle */
+ t_topology* top; /* topology */
+ int ePBC;
+ int * index, /* indices for solvent group */
+ *micelle = nullptr;
+ gmx_bool bMicel = FALSE; /* think we're a micel */
+ t_filenm fnm[] = {
+ /* files for g_order */
+ { efTRX, "-f", nullptr, ffREAD }, /* trajectory file */
+ { efNDX, nullptr, nullptr, ffREAD }, /* index file */
+ { efNDX, "-nm", nullptr, ffOPTRD }, /* index with micelle atoms */
+ { efTPR, nullptr, nullptr, ffREAD }, /* topology file */
+ { efXVG, "-o", "order", ffWRITE }, /* xvgr output file */
};
#define NFILE asize(fnm)
- if (!parse_common_args(&argc, argv,
- PCA_CAN_VIEW | PCA_CAN_TIME, NFILE,
- fnm, asize(pa), pa, asize(desc), desc, asize(bugs), bugs, &oenv))
+ if (!parse_common_args(&argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME, NFILE, fnm, asize(pa), pa,
+ asize(desc), desc, asize(bugs), bugs, &oenv))
{
return 0;
}
rd_index(opt2fn("-nm", NFILE, fnm), 1, &nmic, &micelle, &micname);
}
- calc_h2order(ftp2fn(efTRX, NFILE, fnm), index, ngx, &slDipole, &slOrder,
- &slWidth, &nslices, top, ePBC, axis, bMicel, micelle, nmic,
- oenv);
+ calc_h2order(ftp2fn(efTRX, NFILE, fnm), index, ngx, &slDipole, &slOrder, &slWidth, &nslices,
+ top, ePBC, axis, bMicel, micelle, nmic, oenv);
- h2order_plot(slDipole, slOrder, opt2fn("-o", NFILE, fnm), nslices,
- slWidth, oenv);
+ h2order_plot(slDipole, slOrder, opt2fn("-o", NFILE, fnm), nslices, slWidth, oenv);
do_view(oenv, opt2fn("-o", NFILE, fnm), "-nxy"); /* view xvgr file */
biod.pnpi.spb.ru / alexxy / gromacs.git / blobdiff