if(GMX_BUILD_HOSTNAME AND NOT "${GMX_BUILD_HOSTNAME}" STREQUAL "${TMP_HOSTNAME}")
message(WARNING "
The CMake cache, probably generated on a different host (${GMX_BUILD_HOSTNAME}),
- is being reused! This could lead to inconsitencies; therefore, it is
+ is being reused! This could lead to inconsistencies; therefore, it is
recommended to regenerate the cache!")
endif()
set(GMX_BUILD_HOSTNAME "${TMP_HOSTNAME}" CACHE INTERNAL
from the whole bilayer. The whole bilayer can be taken as reference
group in that case, but it might also be of interest to define the
reaction coordinate for the PMF more locally. The {\tt .mdp} option
-{\tt pull_geometry = cylinder} does not
+{\tt pull-geometry = cylinder} does not
use all the atoms of the reference group, but instead dynamically only those
-within a cylinder with radius {\tt r_1} around the pull vector going
+within a cylinder with radius {\tt pull-r1} around the pull vector going
through the pull group. This only
works for distances defined in one dimension, and the cylinder is
oriented with its long axis along this one dimension. A second cylinder
-can be defined with {\tt r_0}, with a linear switch function that weighs
-the contribution of atoms between {\tt r_0} and {\tt r_1} with
+can be defined with {\tt pull-r0}, with a linear switch function that weighs
+the contribution of atoms between {\tt pull-r0} and {\tt pull-r1} with
distance. This smooths the effects of atoms moving in and out of the
cylinder (which causes jumps in the pull forces).
\subsection{Usage}
To apply enforced rotation, the particles $i$ that are to
be subjected to one of the rotation potentials are defined via index groups
-{\tt rot\_group0}, {\tt rot\_group1}, etc., in the {\tt .mdp} input file.
+{\tt rot-group0}, {\tt rot-group1}, etc., in the {\tt .mdp} input file.
The reference positions $\ve{y}_i^0$ are
read from a special {\tt .trr} file provided to {\tt grompp}. If no such file is found,
$\ve{x}_i(t=0)$ are used as reference positions and written to {\tt .trr} such
that they can be used for subsequent setups. All parameters of the potentials
such as $k$, $\epsilon'$, etc. (\tabref{vars}) are provided as {\tt .mdp}
-parameters; {\tt rot\_type} selects the type of the potential.
-The option {\tt rot\_massw} allows to choose whether or not to use
+parameters; {\tt rot-type} selects the type of the potential.
+The option {\tt rot-massw} allows to choose whether or not to use
mass-weighted averaging.
For the flexible potentials, a cutoff value $g_n^\mathrm{min}$
(typically $g_n^\mathrm{min}=0.001$) makes shure that only
\begin{tabular}{l>{$}l<{$}rccccccc}
\hline
parameter & & & $k$ & $\hat{\ve{v}}$ & $\ve{u}$ & $\omega$ & $\epsilon'$ & $\Delta x$ & $g_n^\mathrm{min}$ \\
-\multicolumn{3}{l}{{\tt .mdp} input variable name} & \smtt{k} & \smtt{vec} & \smtt{pivot} & \smtt{rate} & \smtt{eps} & \smtt{slab\_dist} & \smtt{min\_gauss} \\
+\multicolumn{3}{l}{{\tt .mdp} input variable name} & \smtt{k} & \smtt{vec} & \smtt{pivot} & \smtt{rate} & \smtt{eps} & \smtt{slab-dist} & \smtt{min-gauss} \\
unit & & & \kunit & - & nm & $^\circ$/ps & nm$^2$ & nm & \,-\, \\[1mm]
\hline \multicolumn{2}{l}{fixed axis potentials:} & eqn.\\
isotropic & V\rotiso & (\ref{eqn:potiso}) & {\sf x} & {\sf x} & {\sf x} & {\sf x} & - & - & - \\
$\theta_\mathrm{fit}(t,n)$ is calculated as in \eqnref{rmsdfit}) from the
Gaussian-weighted positions.
-For all angles, the {\tt .mdp} input option {\tt rot\_fit\_method} controls
+For all angles, the {\tt .mdp} input option {\tt rot-fit-method} controls
whether a normal RMSD fit is performed or whether for the fit each
position $\ve{x}_i$ is put at the same distance to the rotation axis as its
reference counterpart $\ve{y}_i^0$. In the latter case, the RMSD
\subsubsection*{Angle Determination by Searching the Energy Minimum}
-Alternatively, for {\tt rot\_fit\_method = potential}, the angle of the rotation
+Alternatively, for {\tt rot-fit-method = potential}, the angle of the rotation
group is determined as the angle for which the rotation potential energy is minimal.
Therefore, the used rotation potential is additionally evaluated for a set of angles
around the current reference angle. In this case, the {\tt rotangles.log} output file
{\small
\begin{verbatim}
swapcoords = Z ; Swap positions: no, X, Y, Z
-swap_frequency = 100 ; Swap attempt frequency
+swap-frequency = 100 ; Swap attempt frequency
\end{verbatim}}
Choose {\tt Z} if your membrane is in the $xy$-plane (\figref{compelsetup} A, B).
Ions will be exchanged between compartments depending on their $z$-positions alone.
-{\tt swap_frequency} determines how often a swap attempt will be made.
+{\tt swap-frequency} determines how often a swap attempt will be made.
This step requires that the positions of the ions, solvent, and swap groups are
-communicated between the parallel ranks, so if chosen too small it can decrease the simulation
+communicated between the parallel processes, so if chosen too small it can decrease the simulation
performance.
{\small
\begin{verbatim}
-split_group0 = channel0 ; Defines compartment boundary
-split_group1 = channel1 ; Defines other compartment boundary
-massw_split0 = no ; use mass-weighted center?
-massw_split1 = no
+split-group0 = channel0 ; Defines compartment boundary
+split-group1 = channel1 ; Defines other compartment boundary
+massw-split0 = no ; use mass-weighted center?
+massw-split1 = no
\end{verbatim}}
-{\tt split_group0} and {\tt split_group1} are two index groups that define the boundaries
+{\tt split-group0} and {\tt split-group1} are two index groups that define the boundaries
between the two compartments, which are usually the centers of the channels.
-If {\tt massw_split0} or {\tt massw_split1} are set to {\tt yes}, the center of mass
+If {\tt massw-split0} or {\tt massw-split1} are set to {\tt yes}, the center of mass
of each index group is used as boundary, here in $z$-direction. Otherwise, the
geometrical centers will be used ($\times$ in \figref{compelsetup} A). If, such as here, a membrane
channel is selected as split group, the center of the channel will define the dividing
must be defined in the index file.
{\small
\begin{verbatim}
-swap_group = NA+_CL- ; Ions to be included in exchange
-solvent_group = SOL ; Group name of solvent molecules
-cyl0_r = 5.0 ; Split cylinder 0: pore radius (nm)
-cyl0_up = 0.75 ; Split cylinder 0 upper extension (nm)
-cyl0_down = 0.75 ; Split cylinder 0 lower extension (nm)
-cyl1_r = 5.0 ; same for other channel
-cyl1_up = 0.75
-cyl1_down = 0.75
-coupl_steps = 10 ; Average over these many swap steps
+swap-group = NA+_CL- ; Ions to be included in exchange
+solvent-group = SOL ; Group name of solvent molecules
+cyl0-r = 5.0 ; Split cylinder 0: pore radius (nm)
+cyl0-up = 0.75 ; Split cylinder 0 upper extension (nm)
+cyl0-down = 0.75 ; Split cylinder 0 lower extension (nm)
+cyl1-r = 5.0 ; same for other channel
+cyl1-up = 0.75
+cyl1-down = 0.75
+coupl-steps = 10 ; Average over these many swap steps
threshold = 1 ; Do not swap if < threshold
\end{verbatim}}
-{\tt swap_group} identifies the index group of ions that
-should be involved in the flux and exchange cycles, {\tt solvent_group} defines the solvent
+{\tt swap-group} identifies the index group of ions that
+should be involved in the flux and exchange cycles, {\tt solvent-group} defines the solvent
group with which they are swapped. The cylinder options only influence the counting of
ions, i.e., ions will be counted as having traveled through either channel 0 or channel 1
according to the definition of (channel) cylinder radius, upper and lower extension,
flux or exchange, but will provide you with the ion permeation numbers across
each of the channels. Note that an ion can only be counted as passing through a particular
channel if it is detected \emph{within} the defined split cylinder in a swap step.
-If {\tt swap_frequency} is chosen too high, a particular ion may be detected in compartment {\bf A}
+If {\tt swap-frequency} is chosen too high, a particular ion may be detected in compartment {\bf A}
in one swap step, and in compartment {\bf B} in the following swap step, so it will be unclear
through which of the channels it has passed.
-{\tt coupl_steps} sets the number of swap attempt steps. A discrepancy between
+{\tt coupl-steps} sets the number of swap attempt steps. A discrepancy between
actual and reference ion numbers in each compartment must persist over this many attempts
-before an actual exchange takes place. If {\tt coupl_steps} is set to 1, then the momentary ion distribution determines
-whether ions are exchanged. {\tt coupl_steps} \textgreater\ 1 will use the time-average
+before an actual exchange takes place. If {\tt coupl-steps} is set to 1, then the momentary ion distribution determines
+whether ions are exchanged. {\tt coupl-steps} \textgreater\ 1 will use the time-average
of ion distributions over the selected number of attempt steps instead. This can be useful, for example,
when ions diffuse near compartment boundaries, which would lead to numerous unproductive
ion exchanges. A {\tt threshold} of 1 means that a swap is performed if the average ion
[ dihedrals ]
CA C +N +H
O C +N +H
- -C N CA CB
-C N CA C
CA C +N +CA
O C +N +CA
comm->cutoff = std::max(comm->cutoff, comm->cutoff_mbody);
}
}
- comm->cellsize_limit = std::max(comm->cellsize_limit, r_bonded_limit);
if (fplog)
{
fprintf(fplog,
"Minimum cell size due to bonded interactions: %.3f nm\n",
- comm->cellsize_limit);
+ r_bonded_limit);
}
+ comm->cellsize_limit = std::max(comm->cellsize_limit, r_bonded_limit);
}
if (dd->bInterCGcons && rconstr <= 0)
{
if (ftpout == efTNG)
{
+ if (ftpout != ftpin)
+ {
+ gmx_fatal(FARGS, "When writing TNG the input file format must also be TNG");
+ }
if (bIndex)
{
trjtools_gmx_prepare_tng_writing(out_file, 'w', NULL, &trxout,
PR("rot-min-gauss", rotg->min_gaussian);
PR("rot-eps", rotg->eps);
PS("rot-fit-method", EROTFIT(rotg->eFittype));
- PI("rot_potfit_nstep", rotg->PotAngle_nstep);
- PR("rot_potfit_step", rotg->PotAngle_step);
+ PI("rot-potfit-nstep", rotg->PotAngle_nstep);
+ PR("rot-potfit-step", rotg->PotAngle_step);
}
static void pr_rot(FILE *fp, int indent, t_rot *rot)
real r_buffer,
real rlist, real boxvol)
{
- double drift_tot, pot1, pot2, pot3, pot;
- int i, j;
- real s2i_2d, s2i_3d, s2j_2d, s2j_3d, s2, s;
- int ti, tj;
- real md1, d2, md3;
- real sc_fac, rsh, rsh2;
- double c_exp, c_erfc;
+ /* Erfc(8)=1e-29, use this limit so we have some space for arithmetic
+ * on the result when using float precision.
+ */
+ const real erfc_arg_max = 8.0;
+
+ double drift_tot, pot1, pot2, pot3, pot;
+ int i, j;
+ real s2i_2d, s2i_3d, s2j_2d, s2j_3d, s2, s;
+ int ti, tj;
+ real md1, d2, md3;
+ real sc_fac, rsh, rsh2;
+ double c_exp, c_erfc;
drift_tot = 0;
rsh = r_buffer;
sc_fac = 1.0;
- /* For constraints: adapt r and scaling for the Gaussian */
- if (att[i].prop.bConstr)
- {
- real sh, sc;
- approx_2dof(s2i_2d, r_buffer*s2i_2d/s2, &sh, &sc);
- rsh += sh;
- sc_fac *= sc;
+ if (rsh*rsh > 2*s2*erfc_arg_max*erfc_arg_max)
+ {
+ /* Erfc might run out of float and become 0, somewhat before
+ * c_exp becomes 0. To avoid this and to avoid NaN in
+ * approx_2dof, we set both c_expc and c_erfc to zero.
+ * In any relevant case this has no effect on the results,
+ * since c_exp < 6e-29, so the displacement is completely
+ * negligible for such atom pairs (and an overestimate).
+ * In nearly all use cases, there will be other atom
+ * pairs that contribute much more to the total, so zeroing
+ * this particular contribution has no effect at all.
+ */
+ c_exp = 0;
+ c_erfc = 0;
}
- if (att[j].prop.bConstr)
+ else
{
- real sh, sc;
+ /* For constraints: adapt r and scaling for the Gaussian */
+ if (att[i].prop.bConstr)
+ {
+ real sh, sc;
- approx_2dof(s2j_2d, r_buffer*s2j_2d/s2, &sh, &sc);
- rsh += sh;
- sc_fac *= sc;
- }
+ approx_2dof(s2i_2d, r_buffer*s2i_2d/s2, &sh, &sc);
+ rsh += sh;
+ sc_fac *= sc;
+ }
+ if (att[j].prop.bConstr)
+ {
+ real sh, sc;
- /* Exact contribution of an atom pair with Gaussian displacement
- * with sigma s to the energy drift for a potential with
- * derivative -md and second derivative dd at the cut-off.
- * The only catch is that for potentials that change sign
- * near the cut-off there could be an unlucky compensation
- * of positive and negative energy drift.
- * Such potentials are extremely rare though.
- *
- * Note that pot has unit energy*length, as the linear
- * atom density still needs to be put in.
- */
- c_exp = exp(-rsh*rsh/(2*s2))/sqrt(2*M_PI);
- c_erfc = 0.5*gmx_erfc(rsh/(sqrt(2*s2)));
+ approx_2dof(s2j_2d, r_buffer*s2j_2d/s2, &sh, &sc);
+ rsh += sh;
+ sc_fac *= sc;
+ }
+
+ /* Exact contribution of an atom pair with Gaussian displacement
+ * with sigma s to the energy drift for a potential with
+ * derivative -md and second derivative dd at the cut-off.
+ * The only catch is that for potentials that change sign
+ * near the cut-off there could be an unlucky compensation
+ * of positive and negative energy drift.
+ * Such potentials are extremely rare though.
+ *
+ * Note that pot has unit energy*length, as the linear
+ * atom density still needs to be put in.
+ */
+ c_exp = exp(-rsh*rsh/(2*s2))/sqrt(2*M_PI);
+ c_erfc = 0.5*gmx_erfc(rsh/(sqrt(2*s2)));
+ }
s = sqrt(s2);
rsh2 = rsh*rsh;
md1/2*((rsh2 + s2)*c_erfc - rsh*s*c_exp);
pot2 = sc_fac*
d2/6*(s*(rsh2 + 2*s2)*c_exp - rsh*(rsh2 + 3*s2)*c_erfc);
- pot3 =
+ pot3 = sc_fac*
md3/24*((rsh2*rsh2 + 6*rsh2*s2 + 3*s2*s2)*c_erfc - rsh*s*(rsh2 + 5*s2)*c_exp);
pot = pot1 + pot2 + pot3;
if (debug)
{
- fprintf(debug, "ib %3d %3d %3d rb %.3f %dx%d fac %.3f drift %f\n",
+ fprintf(debug, "ib %3d %3d %3d rb %.3f %dx%d fac %.3f drift %.1e\n",
ib0, ib, ib1, rb,
list_setup->cluster_size_i, list_setup->cluster_size_j,
nb_clust_frac_pairs_not_in_list_at_cutoff,
/* read rotation parameters */
CTYPE("Output frequency for angle, torque and rotation potential energy for the whole group");
- ITYPE("rot_nstrout", rot->nstrout, 100);
+ ITYPE("rot-nstrout", rot->nstrout, 100);
CTYPE("Output frequency for per-slab data (angles, torques and slab centers)");
- ITYPE("rot_nstsout", rot->nstsout, 1000);
+ ITYPE("rot-nstsout", rot->nstsout, 1000);
CTYPE("Number of rotation groups");
- ITYPE("rot_ngroups", rot->ngrp, 1);
+ ITYPE("rot-ngroups", rot->ngrp, 1);
if (rot->ngrp < 1)
{
- gmx_fatal(FARGS, "rot_ngroups should be >= 1");
+ gmx_fatal(FARGS, "rot-ngroups should be >= 1");
}
snew(rot->grp, rot->ngrp);
rotg = &rot->grp[g];
snew(grpbuf[g], STRLEN);
CTYPE("Rotation group name");
- sprintf(buf, "rot_group%d", g);
+ sprintf(buf, "rot-group%d", g);
STYPE(buf, grpbuf[g], "");
CTYPE("Rotation potential. Can be iso, iso-pf, pm, pm-pf, rm, rm-pf, rm2, rm2-pf, flex, flex-t, flex2, flex2-t");
- sprintf(buf, "rot_type%d", g);
+ sprintf(buf, "rot-type%d", g);
ETYPE(buf, rotg->eType, erotg_names);
CTYPE("Use mass-weighting of the rotation group positions");
- sprintf(buf, "rot_massw%d", g);
+ sprintf(buf, "rot-massw%d", g);
ETYPE(buf, rotg->bMassW, yesno_names);
CTYPE("Rotation vector, will get normalized");
- sprintf(buf, "rot_vec%d", g);
+ sprintf(buf, "rot-vec%d", g);
STYPE(buf, s_vec, "1.0 0.0 0.0");
string2dvec(s_vec, vec);
/* Normalize the rotation vector */
}
else
{
- sprintf(warn_buf, "rot_vec%d = 0", g);
+ sprintf(warn_buf, "rot-vec%d = 0", g);
warning_error(wi, warn_buf);
}
fprintf(stderr, "%s Group %d (%s) normalized rot. vector: %f %f %f\n",
}
CTYPE("Pivot point for the potentials iso, pm, rm, and rm2 (nm)");
- sprintf(buf, "rot_pivot%d", g);
+ sprintf(buf, "rot-pivot%d", g);
STYPE(buf, s_vec, "0.0 0.0 0.0");
clear_dvec(vec);
if ( (rotg->eType == erotgISO) || (rotg->eType == erotgPM) || (rotg->eType == erotgRM) || (rotg->eType == erotgRM2) )
}
CTYPE("Rotation rate (degree/ps) and force constant (kJ/(mol*nm^2))");
- sprintf(buf, "rot_rate%d", g);
+ sprintf(buf, "rot-rate%d", g);
RTYPE(buf, rotg->rate, 0.0);
- sprintf(buf, "rot_k%d", g);
+ sprintf(buf, "rot-k%d", g);
RTYPE(buf, rotg->k, 0.0);
if (rotg->k <= 0.0)
{
- sprintf(warn_buf, "rot_k%d <= 0", g);
+ sprintf(warn_buf, "rot-k%d <= 0", g);
warning_note(wi, warn_buf);
}
CTYPE("Slab distance for flexible axis rotation (nm)");
- sprintf(buf, "rot_slab_dist%d", g);
+ sprintf(buf, "rot-slab-dist%d", g);
RTYPE(buf, rotg->slab_dist, 1.5);
if (rotg->slab_dist <= 0.0)
{
- sprintf(warn_buf, "rot_slab_dist%d <= 0", g);
+ sprintf(warn_buf, "rot-slab-dist%d <= 0", g);
warning_error(wi, warn_buf);
}
CTYPE("Minimum value of Gaussian function for the force to be evaluated (for flex* potentials)");
- sprintf(buf, "rot_min_gauss%d", g);
+ sprintf(buf, "rot-min-gauss%d", g);
RTYPE(buf, rotg->min_gaussian, 1e-3);
if (rotg->min_gaussian <= 0.0)
{
- sprintf(warn_buf, "rot_min_gauss%d <= 0", g);
+ sprintf(warn_buf, "rot-min-gauss%d <= 0", g);
warning_error(wi, warn_buf);
}
CTYPE("Value of additive constant epsilon' (nm^2) for rm2* and flex2* potentials");
- sprintf(buf, "rot_eps%d", g);
+ sprintf(buf, "rot-eps%d", g);
RTYPE(buf, rotg->eps, 1e-4);
if ( (rotg->eps <= 0.0) && (rotg->eType == erotgRM2 || rotg->eType == erotgFLEX2) )
{
- sprintf(warn_buf, "rot_eps%d <= 0", g);
+ sprintf(warn_buf, "rot-eps%d <= 0", g);
warning_error(wi, warn_buf);
}
CTYPE("Fitting method to determine angle of rotation group (rmsd, norm, or potential)");
- sprintf(buf, "rot_fit_method%d", g);
+ sprintf(buf, "rot-fit-method%d", g);
ETYPE(buf, rotg->eFittype, erotg_fitnames);
CTYPE("For fit type 'potential', nr. of angles around the reference for which the pot. is evaluated");
- sprintf(buf, "rot_potfit_nsteps%d", g);
+ sprintf(buf, "rot-potfit-nsteps%d", g);
ITYPE(buf, rotg->PotAngle_nstep, 21);
if ( (rotg->eFittype == erotgFitPOT) && (rotg->PotAngle_nstep < 1) )
{
- sprintf(warn_buf, "rot_potfit_nsteps%d < 1", g);
+ sprintf(warn_buf, "rot-potfit-nsteps%d < 1", g);
warning_error(wi, warn_buf);
}
CTYPE("For fit type 'potential', distance in degrees between two consecutive angles");
- sprintf(buf, "rot_potfit_step%d", g);
+ sprintf(buf, "rot-potfit-step%d", g);
RTYPE(buf, rotg->PotAngle_step, 0.25);
}
typedef struct {
rvec c0; /* The lower corner of the (local) grid */
rvec c1; /* The upper corner of the (local) grid */
+ rvec size; /* c1 - c0 */
real atom_density; /* The atom number density for the local grid */
gmx_bool bSimple; /* Is this grid simple or super/sub */
copy_rvec(corner0, grid->c0);
copy_rvec(corner1, grid->c1);
+ copy_rvec(size, grid->size);
return nc_max;
}
/* Prints the average bb size, used for debug output */
static void print_bbsizes_simple(FILE *fp,
- const nbnxn_search_t nbs,
const nbnxn_grid_t *grid)
{
int c, d;
}
dsvmul(1.0/grid->nc, ba, ba);
- fprintf(fp, "ns bb: %4.2f %4.2f %4.2f %4.2f %4.2f %4.2f rel %4.2f %4.2f %4.2f\n",
- nbs->box[XX][XX]/grid->ncx,
- nbs->box[YY][YY]/grid->ncy,
- nbs->box[ZZ][ZZ]*grid->ncx*grid->ncy/grid->nc,
+ fprintf(fp, "ns bb: grid %4.2f %4.2f %4.2f abs %4.2f %4.2f %4.2f rel %4.2f %4.2f %4.2f\n",
+ grid->sx,
+ grid->sy,
+ grid->na_c/(grid->atom_density*grid->sx*grid->sy),
ba[XX], ba[YY], ba[ZZ],
- ba[XX]*grid->ncx/nbs->box[XX][XX],
- ba[YY]*grid->ncy/nbs->box[YY][YY],
- ba[ZZ]*grid->nc/(grid->ncx*grid->ncy*nbs->box[ZZ][ZZ]));
+ ba[XX]/grid->sx,
+ ba[YY]/grid->sy,
+ ba[ZZ]/(grid->na_c/(grid->atom_density*grid->sx*grid->sy)));
}
/* Prints the average bb size, used for debug output */
static void print_bbsizes_supersub(FILE *fp,
- const nbnxn_search_t nbs,
const nbnxn_grid_t *grid)
{
int ns, c, s;
}
dsvmul(1.0/ns, ba, ba);
- fprintf(fp, "ns bb: %4.2f %4.2f %4.2f %4.2f %4.2f %4.2f rel %4.2f %4.2f %4.2f\n",
- nbs->box[XX][XX]/(grid->ncx*GPU_NSUBCELL_X),
- nbs->box[YY][YY]/(grid->ncy*GPU_NSUBCELL_Y),
- nbs->box[ZZ][ZZ]*grid->ncx*grid->ncy/(grid->nc*GPU_NSUBCELL_Z),
+ fprintf(fp, "ns bb: grid %4.2f %4.2f %4.2f abs %4.2f %4.2f %4.2f rel %4.2f %4.2f %4.2f\n",
+ grid->sx/GPU_NSUBCELL_X,
+ grid->sy/GPU_NSUBCELL_Y,
+ grid->na_sc/(grid->atom_density*grid->sx*grid->sy*GPU_NSUBCELL_Z),
ba[XX], ba[YY], ba[ZZ],
- ba[XX]*grid->ncx*GPU_NSUBCELL_X/nbs->box[XX][XX],
- ba[YY]*grid->ncy*GPU_NSUBCELL_Y/nbs->box[YY][YY],
- ba[ZZ]*grid->nc*GPU_NSUBCELL_Z/(grid->ncx*grid->ncy*nbs->box[ZZ][ZZ]));
+ ba[XX]*GPU_NSUBCELL_X/grid->sx,
+ ba[YY]*GPU_NSUBCELL_Y/grid->sy,
+ ba[ZZ]/(grid->na_sc/(grid->atom_density*grid->sx*grid->sy*GPU_NSUBCELL_Z)));
}
/* Potentially sorts atoms on LJ coefficients !=0 and ==0.
sort_atoms(ZZ, FALSE, dd_zone,
nbs->a+ash, na, x,
grid->c0[ZZ],
- 1.0/nbs->box[ZZ][ZZ], ncz*grid->na_sc,
+ 1.0/grid->size[ZZ], ncz*grid->na_sc,
sort_work);
/* Fill the ncz cells in this column */
sort_atoms(ZZ, FALSE, dd_zone,
nbs->a+ash, na, x,
grid->c0[ZZ],
- 1.0/nbs->box[ZZ][ZZ], ncz*grid->na_sc,
+ 1.0/grid->size[ZZ], ncz*grid->na_sc,
sort_work);
/* This loop goes over the supercells and subcells along z at once */
{
if (grid->bSimple)
{
- print_bbsizes_simple(debug, nbs, grid);
+ print_bbsizes_simple(debug, grid);
}
else
{
fprintf(debug, "ns non-zero sub-cells: %d average atoms %.2f\n",
grid->nsubc_tot, (a1-a0)/(double)grid->nsubc_tot);
- print_bbsizes_supersub(debug, nbs, grid);
+ print_bbsizes_supersub(debug, grid);
}
}
}
fprintf(fp, "nbl na_sc %d rl %g ncp %d per cell %.1f atoms %.1f ratio %.2f\n",
nbl->na_sc, rl, nbl->ncj, nbl->ncj/(double)grid->nc,
nbl->ncj/(double)grid->nc*grid->na_sc,
- nbl->ncj/(double)grid->nc*grid->na_sc/(0.5*4.0/3.0*M_PI*rl*rl*rl*grid->nc*grid->na_sc/det(nbs->box)));
+ nbl->ncj/(double)grid->nc*grid->na_sc/(0.5*4.0/3.0*M_PI*rl*rl*rl*grid->nc*grid->na_sc/(grid->size[XX]*grid->size[YY]*grid->size[ZZ])));
fprintf(fp, "nbl average j cell list length %.1f\n",
0.25*nbl->ncj/(double)nbl->nci);
fprintf(fp, "nbl na_c %d rl %g ncp %d per cell %.1f atoms %.1f ratio %.2f\n",
nbl->na_ci, rl, nbl->nci_tot, nbl->nci_tot/(double)grid->nsubc_tot,
nbl->nci_tot/(double)grid->nsubc_tot*grid->na_c,
- nbl->nci_tot/(double)grid->nsubc_tot*grid->na_c/(0.5*4.0/3.0*M_PI*rl*rl*rl*grid->nsubc_tot*grid->na_c/det(nbs->box)));
+ nbl->nci_tot/(double)grid->nsubc_tot*grid->na_c/(0.5*4.0/3.0*M_PI*rl*rl*rl*grid->nsubc_tot*grid->na_c/(grid->size[XX]*grid->size[YY]*grid->size[ZZ])));
fprintf(fp, "nbl average j super cell list length %.1f\n",
0.25*nbl->ncj4/(double)nbl->nsci);
{
xp[i][d] = state->x[i][d] + fac*f_lr[i][d]*md->invmass[i];
}
+ else
+ {
+ xp[i][d] = state->x[i][d];
+ }
}
}
constrain(NULL, FALSE, FALSE, constr, idef, ir, NULL, cr, step, 0, 1.0, md,
fp = xvgropen(fn, "Rotation angles and energy", "Time (ps)", "angles (degrees) and energies (kJ/mol)", oenv);
fprintf(fp, "# Output of enforced rotation data is written in intervals of %d time step%s.\n#\n", rot->nstrout, rot->nstrout > 1 ? "s" : "");
fprintf(fp, "# The scalar tau is the torque (kJ/mol) in the direction of the rotation vector v.\n");
- fprintf(fp, "# To obtain the vectorial torque, multiply tau with the group's rot_vec.\n");
+ fprintf(fp, "# To obtain the vectorial torque, multiply tau with the group's rot-vec.\n");
fprintf(fp, "# For flexible groups, tau(t,n) from all slabs n have been summed in a single value tau(t) here.\n");
fprintf(fp, "# The torques tau(t,n) are found in the rottorque.log (-rt) output file\n");
fprintf(fp, "#\n");
fprintf(fp, "# ROTATION GROUP %d, potential type '%s':\n", g, erotg_names[rotg->eType]);
- fprintf(fp, "# rot_massw%d %s\n", g, yesno_names[rotg->bMassW]);
- fprintf(fp, "# rot_vec%d %12.5e %12.5e %12.5e\n", g, rotg->vec[XX], rotg->vec[YY], rotg->vec[ZZ]);
- fprintf(fp, "# rot_rate%d %12.5e degrees/ps\n", g, rotg->rate);
- fprintf(fp, "# rot_k%d %12.5e kJ/(mol*nm^2)\n", g, rotg->k);
+ fprintf(fp, "# rot-massw%d %s\n", g, yesno_names[rotg->bMassW]);
+ fprintf(fp, "# rot-vec%d %12.5e %12.5e %12.5e\n", g, rotg->vec[XX], rotg->vec[YY], rotg->vec[ZZ]);
+ fprintf(fp, "# rot-rate%d %12.5e degrees/ps\n", g, rotg->rate);
+ fprintf(fp, "# rot-k%d %12.5e kJ/(mol*nm^2)\n", g, rotg->k);
if (rotg->eType == erotgISO || rotg->eType == erotgPM || rotg->eType == erotgRM || rotg->eType == erotgRM2)
{
- fprintf(fp, "# rot_pivot%d %12.5e %12.5e %12.5e nm\n", g, rotg->pivot[XX], rotg->pivot[YY], rotg->pivot[ZZ]);
+ fprintf(fp, "# rot-pivot%d %12.5e %12.5e %12.5e nm\n", g, rotg->pivot[XX], rotg->pivot[YY], rotg->pivot[ZZ]);
}
if (bFlex)
{
- fprintf(fp, "# rot_slab_distance%d %f nm\n", g, rotg->slab_dist);
- fprintf(fp, "# rot_min_gaussian%d %12.5e\n", g, rotg->min_gaussian);
+ fprintf(fp, "# rot-slab-distance%d %f nm\n", g, rotg->slab_dist);
+ fprintf(fp, "# rot-min-gaussian%d %12.5e\n", g, rotg->min_gaussian);
}
/* Output the centers of the rotation groups for the pivot-free potentials */
if ( (rotg->eType == erotgRM2) || (rotg->eType == erotgFLEX2) || (rotg->eType == erotgFLEX2T) )
{
- fprintf(fp, "# rot_eps%d %12.5e nm^2\n", g, rotg->eps);
+ fprintf(fp, "# rot-eps%d %12.5e nm^2\n", g, rotg->eps);
}
if (erotgFitPOT == rotg->eFittype)
{
fprintf(fp, "# Rotation group %d (%s), slab distance %f nm.\n", g, erotg_names[rotg->eType], rotg->slab_dist);
fprintf(fp, "# The scalar tau is the torque (kJ/mol) in the direction of the rotation vector.\n");
fprintf(fp, "# To obtain the vectorial torque, multiply tau with\n");
- fprintf(fp, "# rot_vec%d %10.3e %10.3e %10.3e\n", g, rotg->vec[XX], rotg->vec[YY], rotg->vec[ZZ]);
+ fprintf(fp, "# rot-vec%d %10.3e %10.3e %10.3e\n", g, rotg->vec[XX], rotg->vec[YY], rotg->vec[ZZ]);
fprintf(fp, "#\n");
}
}
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff