/* Open output file for slab center data. Call on master only */
-static FILE *open_slab_out(const char *fn, t_rot *rot, const output_env_t oenv)
+static FILE *open_slab_out(const char *fn, t_rot *rot)
{
FILE *fp;
int g, i;
/* Call on master only */
-static FILE *open_angles_out(const char *fn, t_rot *rot, const output_env_t oenv)
+static FILE *open_angles_out(const char *fn, t_rot *rot)
{
int g, i;
FILE *fp;
/* Open torque output file and write some information about it's structure.
* Call on master only */
-static FILE *open_torque_out(const char *fn, t_rot *rot, const output_env_t oenv)
+static FILE *open_torque_out(const char *fn, t_rot *rot)
{
FILE *fp;
int g;
}
count--;
- /* Determine the max slab */
+ /* Determine the min slab */
slab = homeslab;
do
{
/* Subtract the slab center from xj */
rvec_sub(xj, xcn, tmpvec2); /* tmpvec2 = xj - xcn */
+
+ /* In rare cases, when an atom position coincides with a slab center
+ * (tmpvec2 == 0) we cannot compute the vector product for sjn.
+ * However, since the atom is located directly on the pivot, this
+ * slab's contribution to the force on that atom will be zero
+ * anyway. Therefore, we directly move on to the next slab. */
+ if ( 0 == norm(tmpvec2) )
+ {
+ continue;
+ }
/* Calculate sjn */
cprod(rotg->vec, tmpvec2, tmpvec); /* tmpvec = v x (xj - xcn) */
/* Subtract the slab center from xj */
rvec_sub(xj, xcn, xj_xcn); /* xj_xcn = xj - xcn */
+ /* In rare cases, when an atom position coincides with a slab center
+ * (xj_xcn == 0) we cannot compute the vector product for qjn.
+ * However, since the atom is located directly on the pivot, this
+ * slab's contribution to the force on that atom will be zero
+ * anyway. Therefore, we directly move on to the next slab. */
+ if ( 0 == norm(xj_xcn) )
+ {
+ continue;
+ }
+
/* Calculate qjn */
cprod(rotg->vec, tmpvec2, tmpvec); /* tmpvec= v x Omega.(yj0-ycn) */
t_rotgrp *rotg, /* The rotation group (inputrec data) */
t_gmx_enfrotgrp *erg, /* The rotation group (data only accessible in this file) */
rvec firstatom, /* First atom after sorting along the rotation vector v */
- rvec lastatom, /* Last atom along v */
- int g) /* The rotation group number */
+ rvec lastatom) /* Last atom along v */
{
erg->slab_first = get_first_slab(rotg, erg->max_beta, firstatom);
erg->slab_last = get_last_slab(rotg, erg->max_beta, lastatom);
+ /* Calculate the slab buffer size, which changes when slab_first changes */
+ erg->slab_buffer = erg->slab_first - erg->slab_first_ref;
+
/* Check whether we have reference data to compare against */
if (erg->slab_first < erg->slab_first_ref)
{
rvec x[], /* The local positions */
matrix box,
double t, /* Time in picoseconds */
- gmx_large_int_t step, /* The time step */
gmx_bool bOutstepRot, /* Output to main rotation output file */
gmx_bool bOutstepSlab) /* Output per-slab data */
{
/* Determine the first relevant slab for the first atom and the last
* relevant slab for the last atom */
- get_firstlast_slab_check(rotg, erg, erg->xc[0], erg->xc[rotg->nat-1], g);
+ get_firstlast_slab_check(rotg, erg, erg->xc[0], erg->xc[rotg->nat-1]);
/* Determine for each slab depending on the min_gaussian cutoff criterium,
* a first and a last atom index inbetween stuff needs to be calculated */
/* springs to the reference atoms. */
static void do_fixed(
t_rotgrp *rotg, /* The rotation group */
- rvec x[], /* The positions */
- matrix box, /* The simulation box */
- double t, /* Time in picoseconds */
- gmx_large_int_t step, /* The time step */
gmx_bool bOutstepRot, /* Output to main rotation output file */
gmx_bool bOutstepSlab) /* Output per-slab data */
{
/* Calculate the radial motion potential and forces */
static void do_radial_motion(
t_rotgrp *rotg, /* The rotation group */
- rvec x[], /* The positions */
- matrix box, /* The simulation box */
- double t, /* Time in picoseconds */
- gmx_large_int_t step, /* The time step */
gmx_bool bOutstepRot, /* Output to main rotation output file */
gmx_bool bOutstepSlab) /* Output per-slab data */
{
t_rotgrp *rotg, /* The rotation group */
rvec x[], /* The positions */
matrix box, /* The simulation box */
- double t, /* Time in picoseconds */
- gmx_large_int_t step, /* The time step */
gmx_bool bOutstepRot, /* Output to main rotation output file */
gmx_bool bOutstepSlab) /* Output per-slab data */
{
t_rotgrp *rotg, /* The rotation group */
rvec x[], /* The positions */
matrix box, /* The simulation box */
- double t, /* Time in picoseconds */
- gmx_large_int_t step, /* The time step */
gmx_bool bOutstepRot, /* Output to main rotation output file */
gmx_bool bOutstepSlab) /* Output per-slab data */
{
}
-/* From the extreme coordinates of the reference group, determine the first
+/* From the extreme positions of the reference group, determine the first
* and last slab of the reference. We can never have more slabs in the real
* simulation than calculated here for the reference.
*/
static void get_firstlast_slab_ref(t_rotgrp *rotg, real mc[], int ref_firstindex, int ref_lastindex)
{
gmx_enfrotgrp_t erg; /* Pointer to enforced rotation group data */
- int first, last, firststart;
+ int first, last;
rvec dummy;
erg = rotg->enfrotgrp;
first = get_first_slab(rotg, erg->max_beta, rotg->x_ref[ref_firstindex]);
last = get_last_slab( rotg, erg->max_beta, rotg->x_ref[ref_lastindex ]);
- firststart = first;
while (get_slab_weight(first, rotg, rotg->x_ref, mc, &dummy) > WEIGHT_MIN)
{
last++;
}
erg->slab_last_ref = last-1;
-
- erg->slab_buffer = firststart - erg->slab_first_ref;
}
/* Determine the smallest and largest coordinate with respect to the rotation vector */
get_firstlast_atom_ref(rotg, &ref_firstindex, &ref_lastindex);
- /* From the extreme coordinates of the reference group, determine the first
+ /* From the extreme positions of the reference group, determine the first
* and last slab of the reference. */
get_firstlast_slab_ref(rotg, erg->mc, ref_firstindex, ref_lastindex);
er->out_slabs = NULL;
if (MASTER(cr) && HaveFlexibleGroups(rot) )
{
- er->out_slabs = open_slab_out(opt2fn("-rs", nfile, fnm), rot, oenv);
+ er->out_slabs = open_slab_out(opt2fn("-rs", nfile, fnm), rot);
}
if (MASTER(cr))
{
if (HaveFlexibleGroups(rot) || HavePotFitGroups(rot) )
{
- er->out_angles = open_angles_out(opt2fn("-ra", nfile, fnm), rot, oenv);
+ er->out_angles = open_angles_out(opt2fn("-ra", nfile, fnm), rot);
}
if (HaveFlexibleGroups(rot) )
{
- er->out_torque = open_torque_out(opt2fn("-rt", nfile, fnm), rot, oenv);
+ er->out_torque = open_torque_out(opt2fn("-rt", nfile, fnm), rot);
}
}
}
-extern void finish_rot(FILE *fplog, t_rot *rot)
+extern void finish_rot(t_rot *rot)
{
gmx_enfrot_t er; /* Pointer to the enforced rotation buffer variables */
case erotgISOPF:
case erotgPM:
case erotgPMPF:
- do_fixed(rotg, x, box, t, step, outstep_rot, outstep_slab);
+ do_fixed(rotg, outstep_rot, outstep_slab);
break;
case erotgRM:
- do_radial_motion(rotg, x, box, t, step, outstep_rot, outstep_slab);
+ do_radial_motion(rotg, outstep_rot, outstep_slab);
break;
case erotgRMPF:
- do_radial_motion_pf(rotg, x, box, t, step, outstep_rot, outstep_slab);
+ do_radial_motion_pf(rotg, x, box, outstep_rot, outstep_slab);
break;
case erotgRM2:
case erotgRM2PF:
- do_radial_motion2(rotg, x, box, t, step, outstep_rot, outstep_slab);
+ do_radial_motion2(rotg, x, box, outstep_rot, outstep_slab);
break;
case erotgFLEXT:
case erotgFLEX2T:
get_center(erg->xc, erg->mc, rotg->nat, erg->xc_center);
svmul(-1.0, erg->xc_center, transvec);
translate_x(erg->xc, rotg->nat, transvec);
- do_flexible(MASTER(cr), er, rotg, g, x, box, t, step, outstep_rot, outstep_slab);
+ do_flexible(MASTER(cr), er, rotg, g, x, box, t, outstep_rot, outstep_slab);
break;
case erotgFLEX:
case erotgFLEX2:
/* Do NOT subtract the center of mass in the low level routines! */
clear_rvec(erg->xc_center);
- do_flexible(MASTER(cr), er, rotg, g, x, box, t, step, outstep_rot, outstep_slab);
+ do_flexible(MASTER(cr), er, rotg, g, x, box, t, outstep_rot, outstep_slab);
break;
default:
gmx_fatal(FARGS, "No such rotation potential.");
biod.pnpi.spb.ru / alexxy / gromacs.git / blobdiff