namespace gmx
{
+//! Indices of the two atoms involved in a single constraint
+struct AtomPair
+{
+ //! \brief Constructor, does not initialize to catch bugs and faster construction
+ AtomPair()
+ {
+ };
+
+ //! Index of atom 1
+ int index1;
+ //! Index of atom 2
+ int index2;
+};
+
//! Unit of work within LINCS.
struct Task
{
//! The reference distance in top B - the r.d. in top A.
std::vector < real, AlignedAllocator < real>> ddist;
//! The atom pairs involved in the constraints.
- std::vector<int> bla;
+ std::vector<AtomPair> atoms;
//! 1/sqrt(invmass1 invmass2).
std::vector < real, AlignedAllocator < real>> blc;
//! As blc, but with all masses 1.
//! Update atomic coordinates when an index is not required.
static void
lincs_update_atoms_noind(int ncons,
- gmx::ArrayRef<const int> bla,
+ gmx::ArrayRef<const AtomPair> atoms,
real preFactor,
gmx::ArrayRef<const real> fac,
gmx::ArrayRef<const gmx::RVec> r,
{
for (int b = 0; b < ncons; b++)
{
- int i = bla[2*b];
- int j = bla[2*b+1];
+ int i = atoms[b].index1;
+ int j = atoms[b].index2;
real mvb = preFactor*fac[b];
real im1 = invmass[i];
real im2 = invmass[j];
{
for (int b = 0; b < ncons; b++)
{
- int i = bla[2*b];
- int j = bla[2*b+1];
+ int i = atoms[b].index1;
+ int j = atoms[b].index2;
real mvb = preFactor*fac[b];
real tmp0 = r[b][0]*mvb;
real tmp1 = r[b][1]*mvb;
//! Update atomic coordinates when an index is required.
static void
lincs_update_atoms_ind(gmx::ArrayRef<const int> ind,
- gmx::ArrayRef<const int> bla,
+ gmx::ArrayRef<const AtomPair> atoms,
real preFactor,
gmx::ArrayRef<const real> fac,
gmx::ArrayRef<const gmx::RVec> r,
{
for (int b : ind)
{
- int i = bla[2*b];
- int j = bla[2*b+1];
+ int i = atoms[b].index1;
+ int j = atoms[b].index2;
real mvb = preFactor*fac[b];
real im1 = invmass[i];
real im2 = invmass[j];
{
for (int b : ind)
{
- int i = bla[2*b];
- int j = bla[2*b+1];
+ int i = atoms[b].index1;
+ int j = atoms[b].index2;
real mvb = preFactor*fac[b];
real tmp0 = r[b][0]*mvb;
real tmp1 = r[b][1]*mvb;
{
/* Single thread, we simply update for all constraints */
lincs_update_atoms_noind(li->nc_real,
- li->bla, preFactor, fac, r, invmass, x);
+ li->atoms, preFactor, fac, r, invmass, x);
}
else
{
* This can be done without a barrier.
*/
lincs_update_atoms_ind(li->task[th].ind,
- li->bla, preFactor, fac, r, invmass, x);
+ li->atoms, preFactor, fac, r, invmass, x);
if (!li->task[li->ntask].ind.empty())
{
#pragma omp master
{
lincs_update_atoms_ind(li->task[li->ntask].ind,
- li->bla, preFactor, fac, r, invmass, x);
+ li->atoms, preFactor, fac, r, invmass, x);
}
}
}
* This function only differs from calc_dr_x_xp_simd below in that
* no constraint length is subtracted and no PBC is used for f. */
static void gmx_simdcall
-calc_dr_x_f_simd(int b0,
- int b1,
- gmx::ArrayRef<const int> bla,
- const rvec * gmx_restrict x,
- const rvec * gmx_restrict f,
- const real * gmx_restrict blc,
- const real * pbc_simd,
- rvec * gmx_restrict r,
- real * gmx_restrict rhs,
- real * gmx_restrict sol)
+calc_dr_x_f_simd(int b0,
+ int b1,
+ gmx::ArrayRef<const AtomPair> atoms,
+ const rvec * gmx_restrict x,
+ const rvec * gmx_restrict f,
+ const real * gmx_restrict blc,
+ const real * pbc_simd,
+ rvec * gmx_restrict r,
+ real * gmx_restrict rhs,
+ real * gmx_restrict sol)
{
assert(b0 % GMX_SIMD_REAL_WIDTH == 0);
for (int i = 0; i < GMX_SIMD_REAL_WIDTH; i++)
{
- offset0[i] = bla[bs*2 + i*2];
- offset1[i] = bla[bs*2 + i*2 + 1];
+ offset0[i] = atoms[bs + i].index1;
+ offset1[i] = atoms[bs + i].index2;
}
gatherLoadUTransposeTSANSafe<3>(reinterpret_cast<const real *>(x), offset0, &x0_S, &y0_S, &z0_S);
ConstraintVariable econq, bool bCalcDHDL,
bool bCalcVir, tensor rmdf)
{
- const int b0 = lincsd->task[th].b0;
- const int b1 = lincsd->task[th].b1;
+ const int b0 = lincsd->task[th].b0;
+ const int b1 = lincsd->task[th].b1;
- gmx::ArrayRef<const int> bla = lincsd->bla;
- gmx::ArrayRef<gmx::RVec> r = lincsd->tmpv;
- gmx::ArrayRef<const int> blnr = lincsd->blnr;
- gmx::ArrayRef<const int> blbnb = lincsd->blbnb;
+ gmx::ArrayRef<const AtomPair> atoms = lincsd->atoms;
+ gmx::ArrayRef<gmx::RVec> r = lincsd->tmpv;
+ gmx::ArrayRef<const int> blnr = lincsd->blnr;
+ gmx::ArrayRef<const int> blbnb = lincsd->blbnb;
- gmx::ArrayRef<const real> blc;
- gmx::ArrayRef<const real> blmf;
+ gmx::ArrayRef<const real> blc;
+ gmx::ArrayRef<const real> blmf;
if (econq != ConstraintVariable::Force)
{
/* Use mass-weighted parameters */
/* Compute normalized x i-j vectors, store in r.
* Compute the inner product of r and xp i-j and store in rhs1.
*/
- calc_dr_x_f_simd(b0, b1, bla, x, f, blc.data(),
+ calc_dr_x_f_simd(b0, b1, atoms, x, f, blc.data(),
pbc_simd,
as_rvec_array(r.data()), rhs1.data(), sol.data());
{
rvec dx;
- pbc_dx_aiuc(pbc, x[bla[2*b]], x[bla[2*b+1]], dx);
+ pbc_dx_aiuc(pbc, x[atoms[b].index1], x[atoms[b].index2], dx);
unitv(dx, r[b]);
}
}
{
rvec dx;
- rvec_sub(x[bla[2*b]], x[bla[2*b+1]], dx);
+ rvec_sub(x[atoms[b].index1], x[atoms[b].index2], dx);
unitv(dx, r[b]);
} /* 16 ncons flops */
}
for (int b = b0; b < b1; b++)
{
- int i = bla[2*b];
- int j = bla[2*b+1];
+ int i = atoms[b].index1;
+ int j = atoms[b].index2;
real mvb = blc[b]*(r[b][0]*(f[i][0] - f[j][0]) +
r[b][1]*(f[i][1] - f[j][1]) +
r[b][2]*(f[i][2] - f[j][2]));
*
* Determine the right-hand side of the matrix equation using coordinates xp. */
static void gmx_simdcall
-calc_dr_x_xp_simd(int b0,
- int b1,
- gmx::ArrayRef<const int> bla,
- const rvec * gmx_restrict x,
- const rvec * gmx_restrict xp,
- const real * gmx_restrict bllen,
- const real * gmx_restrict blc,
- const real * pbc_simd,
- rvec * gmx_restrict r,
- real * gmx_restrict rhs,
- real * gmx_restrict sol)
+calc_dr_x_xp_simd(int b0,
+ int b1,
+ gmx::ArrayRef<const AtomPair> atoms,
+ const rvec * gmx_restrict x,
+ const rvec * gmx_restrict xp,
+ const real * gmx_restrict bllen,
+ const real * gmx_restrict blc,
+ const real * pbc_simd,
+ rvec * gmx_restrict r,
+ real * gmx_restrict rhs,
+ real * gmx_restrict sol)
{
assert(b0 % GMX_SIMD_REAL_WIDTH == 0);
alignas(GMX_SIMD_ALIGNMENT) std::int32_t offset2[GMX_SIMD_REAL_WIDTH];
for (int i = 0; i < GMX_SIMD_REAL_WIDTH; i++)
{
- offset0[i] = bla[bs*2 + i*2];
- offset1[i] = bla[bs*2 + i*2 + 1];
+ offset0[i] = atoms[bs + i].index1;
+ offset1[i] = atoms[bs + i].index2;
}
gatherLoadUTransposeTSANSafe<3>(reinterpret_cast<const real *>(x), offset0, &x0_S, &y0_S, &z0_S);
/*! \brief Determine the distances and right-hand side for the next iteration. */
gmx_unused static void calc_dist_iter(
- int b0,
- int b1,
- gmx::ArrayRef<const int> bla,
- const rvec * gmx_restrict xp,
- const real * gmx_restrict bllen,
- const real * gmx_restrict blc,
- const t_pbc * pbc,
- real wfac,
- real * gmx_restrict rhs,
- real * gmx_restrict sol,
- bool * bWarn)
+ int b0,
+ int b1,
+ gmx::ArrayRef<const AtomPair> atoms,
+ const rvec * gmx_restrict xp,
+ const real * gmx_restrict bllen,
+ const real * gmx_restrict blc,
+ const t_pbc * pbc,
+ real wfac,
+ real * gmx_restrict rhs,
+ real * gmx_restrict sol,
+ bool * bWarn)
{
- int b;
-
- for (b = b0; b < b1; b++)
+ for (int b = b0; b < b1; b++)
{
- real len, len2, dlen2, mvb;
+ real len = bllen[b];
rvec dx;
-
- len = bllen[b];
if (pbc)
{
- pbc_dx_aiuc(pbc, xp[bla[2*b]], xp[bla[2*b+1]], dx);
+ pbc_dx_aiuc(pbc, xp[atoms[b].index1], xp[atoms[b].index2], dx);
}
else
{
- rvec_sub(xp[bla[2*b]], xp[bla[2*b+1]], dx);
+ rvec_sub(xp[atoms[b].index1], xp[atoms[b].index2], dx);
}
- len2 = len*len;
- dlen2 = 2*len2 - ::norm2(dx);
+ real len2 = len*len;
+ real dlen2 = 2*len2 - ::norm2(dx);
if (dlen2 < wfac*len2)
{
/* not race free - see detailed comment in caller */
*bWarn = TRUE;
}
+ real mvb;
if (dlen2 > 0)
{
mvb = blc[b]*(len - dlen2*gmx::invsqrt(dlen2));
#if GMX_SIMD_HAVE_REAL
/*! \brief As calc_dist_iter(), but using SIMD intrinsics. */
static void gmx_simdcall
-calc_dist_iter_simd(int b0,
- int b1,
- gmx::ArrayRef<const int> bla,
- const rvec * gmx_restrict x,
- const real * gmx_restrict bllen,
- const real * gmx_restrict blc,
- const real * pbc_simd,
- real wfac,
- real * gmx_restrict rhs,
- real * gmx_restrict sol,
- bool * bWarn)
+calc_dist_iter_simd(int b0,
+ int b1,
+ gmx::ArrayRef<const AtomPair> atoms,
+ const rvec * gmx_restrict x,
+ const real * gmx_restrict bllen,
+ const real * gmx_restrict blc,
+ const real * pbc_simd,
+ real wfac,
+ real * gmx_restrict rhs,
+ real * gmx_restrict sol,
+ bool * bWarn)
{
SimdReal min_S(GMX_REAL_MIN);
SimdReal two_S(2.0);
SimdReal wfac_S(wfac);
SimdBool warn_S;
- int bs;
-
assert(b0 % GMX_SIMD_REAL_WIDTH == 0);
/* Initialize all to FALSE */
warn_S = (two_S < setZero());
- for (bs = b0; bs < b1; bs += GMX_SIMD_REAL_WIDTH)
+ for (int bs = b0; bs < b1; bs += GMX_SIMD_REAL_WIDTH)
{
SimdReal x0_S, y0_S, z0_S;
SimdReal x1_S, y1_S, z1_S;
for (int i = 0; i < GMX_SIMD_REAL_WIDTH; i++)
{
- offset0[i] = bla[bs*2 + i*2];
- offset1[i] = bla[bs*2 + i*2 + 1];
+ offset0[i] = atoms[bs + i].index1;
+ offset1[i] = atoms[bs + i].index2;
}
gatherLoadUTransposeTSANSafe<3>(reinterpret_cast<const real *>(x), offset0, &x0_S, &y0_S, &z0_S);
real invdt, rvec * gmx_restrict v,
bool bCalcVir, tensor vir_r_m_dr)
{
- const int b0 = lincsd->task[th].b0;
- const int b1 = lincsd->task[th].b1;
-
- gmx::ArrayRef<const int> bla = lincsd->bla;
- gmx::ArrayRef<gmx::RVec> r = lincsd->tmpv;
- gmx::ArrayRef<const int> blnr = lincsd->blnr;
- gmx::ArrayRef<const int> blbnb = lincsd->blbnb;
- gmx::ArrayRef<const real> blc = lincsd->blc;
- gmx::ArrayRef<const real> blmf = lincsd->blmf;
- gmx::ArrayRef<const real> bllen = lincsd->bllen;
- gmx::ArrayRef<real> blcc = lincsd->tmpncc;
- gmx::ArrayRef<real> rhs1 = lincsd->tmp1;
- gmx::ArrayRef<real> rhs2 = lincsd->tmp2;
- gmx::ArrayRef<real> sol = lincsd->tmp3;
- gmx::ArrayRef<real> blc_sol = lincsd->tmp4;
- gmx::ArrayRef<real> mlambda = lincsd->mlambda;
- gmx::ArrayRef<const int> nlocat = lincsd->nlocat;
+ const int b0 = lincsd->task[th].b0;
+ const int b1 = lincsd->task[th].b1;
+
+ gmx::ArrayRef<const AtomPair> atoms = lincsd->atoms;
+ gmx::ArrayRef<gmx::RVec> r = lincsd->tmpv;
+ gmx::ArrayRef<const int> blnr = lincsd->blnr;
+ gmx::ArrayRef<const int> blbnb = lincsd->blbnb;
+ gmx::ArrayRef<const real> blc = lincsd->blc;
+ gmx::ArrayRef<const real> blmf = lincsd->blmf;
+ gmx::ArrayRef<const real> bllen = lincsd->bllen;
+ gmx::ArrayRef<real> blcc = lincsd->tmpncc;
+ gmx::ArrayRef<real> rhs1 = lincsd->tmp1;
+ gmx::ArrayRef<real> rhs2 = lincsd->tmp2;
+ gmx::ArrayRef<real> sol = lincsd->tmp3;
+ gmx::ArrayRef<real> blc_sol = lincsd->tmp4;
+ gmx::ArrayRef<real> mlambda = lincsd->mlambda;
+ gmx::ArrayRef<const int> nlocat = lincsd->nlocat;
#if GMX_SIMD_HAVE_REAL
/* Compute normalized x i-j vectors, store in r.
* Compute the inner product of r and xp i-j and store in rhs1.
*/
- calc_dr_x_xp_simd(b0, b1, bla, x, xp, bllen.data(), blc.data(),
+ calc_dr_x_xp_simd(b0, b1, atoms, x, xp, bllen.data(), blc.data(),
pbc_simd,
as_rvec_array(r.data()), rhs1.data(), sol.data());
for (int b = b0; b < b1; b++)
{
rvec dx;
- pbc_dx_aiuc(pbc, x[bla[2*b]], x[bla[2*b+1]], dx);
+ pbc_dx_aiuc(pbc, x[atoms[b].index1], x[atoms[b].index2], dx);
unitv(dx, r[b]);
- pbc_dx_aiuc(pbc, xp[bla[2*b]], xp[bla[2*b+1]], dx);
+ pbc_dx_aiuc(pbc, xp[atoms[b].index1], xp[atoms[b].index2], dx);
real mvb = blc[b]*(::iprod(r[b], dx) - bllen[b]);
rhs1[b] = mvb;
sol[b] = mvb;
/* Compute normalized i-j vectors */
for (int b = b0; b < b1; b++)
{
- int i = bla[2*b];
- int j = bla[2*b+1];
+ int i = atoms[b].index1;
+ int j = atoms[b].index2;
real tmp0 = x[i][0] - x[j][0];
real tmp1 = x[i][1] - x[j][1];
real tmp2 = x[i][2] - x[j][2];
}
#if GMX_SIMD_HAVE_REAL
- calc_dist_iter_simd(b0, b1, bla,
+ calc_dist_iter_simd(b0, b1, atoms,
xp, bllen.data(), blc.data(), pbc_simd, wfac,
rhs1.data(), sol.data(), bWarn);
#else
- calc_dist_iter(b0, b1, bla, xp,
+ calc_dist_iter(b0, b1, atoms, xp,
bllen.data(), blc.data(), pbc, wfac,
rhs1.data(), sol.data(), bWarn);
/* 20*ncons flops */
int *ncc_triangle,
int *nCrossTaskTriangles)
{
- int i;
-
/* Construct the coupling coefficient matrix blmf */
li_task->ntriangle = 0;
*ncc_triangle = 0;
*nCrossTaskTriangles = 0;
- for (i = li_task->b0; i < li_task->b1; i++)
+ for (int i = li_task->b0; i < li_task->b1; i++)
{
- int a1, a2, n;
-
- a1 = li->bla[2*i];
- a2 = li->bla[2*i+1];
- for (n = li->blnr[i]; (n < li->blnr[i+1]); n++)
+ const int a1 = li->atoms[i].index1;
+ const int a2 = li->atoms[i].index2;
+ for (int n = li->blnr[i]; n < li->blnr[i + 1]; n++)
{
- int k, sign, center, end;
-
- k = li->blbnb[n];
+ const int k = li->blbnb[n];
/* If we are using multiple, independent tasks for LINCS,
* the calls to check_assign_connected should have
*/
assert(li->bTaskDep || (k >= li_task->b0 && k < li_task->b1));
- if (a1 == li->bla[2*k] || a2 == li->bla[2*k+1])
+ int sign;
+ if (a1 == li->atoms[k].index1 || a2 == li->atoms[k].index2)
{
sign = -1;
}
{
sign = 1;
}
- if (a1 == li->bla[2*k] || a1 == li->bla[2*k+1])
+ int center;
+ int end;
+ if (a1 == li->atoms[k].index1 || a1 == li->atoms[k].index2)
{
center = a1;
end = a2;
li->blmf1[n] = sign*0.5;
if (li->ncg_triangle > 0)
{
- int nk, kk;
-
/* Look for constraint triangles */
- for (nk = li->blnr[k]; (nk < li->blnr[k+1]); nk++)
+ for (int nk = li->blnr[k]; nk < li->blnr[k + 1]; nk++)
{
- kk = li->blbnb[nk];
+ const int kk = li->blbnb[nk];
if (kk != i && kk != k &&
- (li->bla[2*kk] == end || li->bla[2*kk+1] == end))
+ (li->atoms[kk].index1 == end ||
+ li->atoms[kk].index2 == end))
{
/* Check if the constraints in this triangle actually
* belong to a different task. We still assign them
/*! \brief Sets the elements in the LINCS matrix. */
static void set_lincs_matrix(Lincs *li, real *invmass, real lambda)
{
- int i;
const real invsqrt2 = 0.7071067811865475244;
- for (i = 0; (i < li->nc); i++)
+ for (int i = 0; (i < li->nc); i++)
{
- int a1, a2;
-
- a1 = li->bla[2*i];
- a2 = li->bla[2*i+1];
- li->blc[i] = gmx::invsqrt(invmass[a1] + invmass[a2]);
- li->blc1[i] = invsqrt2;
+ const int a1 = li->atoms[i].index1;
+ const int a2 = li->atoms[i].index2;
+ li->blc[i] = gmx::invsqrt(invmass[a1] + invmass[a2]);
+ li->blc1[i] = invsqrt2;
}
/* Construct the coupling coefficient matrix blmf */
- int th, ntriangle = 0, ncc_triangle = 0, nCrossTaskTriangles = 0;
+ int ntriangle = 0, ncc_triangle = 0, nCrossTaskTriangles = 0;
#pragma omp parallel for reduction(+: ntriangle, ncc_triangle, nCrossTaskTriangles) num_threads(li->ntask) schedule(static)
- for (th = 0; th < li->ntask; th++)
+ for (int th = 0; th < li->ntask; th++)
{
try
{
/* For each atom set a flag for constraints from each */
for (int b = li_task->b0; b < li_task->b1; b++)
{
- bitmask_set_bit(&atf[li->bla[b*2 ]], th);
- bitmask_set_bit(&atf[li->bla[b*2 + 1]], th);
+ bitmask_set_bit(&atf[li->atoms[b].index1], th);
+ bitmask_set_bit(&atf[li->atoms[b].index2], th);
}
}
/* We let the constraint with the lowest thread index
* operate on atoms with constraints from multiple threads.
*/
- if (bitmask_is_disjoint(atf[li->bla[b*2]], mask) &&
- bitmask_is_disjoint(atf[li->bla[b*2+1]], mask))
+ if (bitmask_is_disjoint(atf[li->atoms[b].index1], mask) &&
+ bitmask_is_disjoint(atf[li->atoms[b].index2], mask))
{
/* Add the constraint to the local atom update index */
li_task->ind.push_back(b);
/* Make an mapping of local topology constraint index to LINCS index */
li->con_index[constraint_index] = con;
- li->bllen0[con] = lenA;
- li->ddist[con] = lenB - lenA;
+ li->bllen0[con] = lenA;
+ li->ddist[con] = lenB - lenA;
/* Set the length to the topology A length */
- li->bllen[con] = lenA;
- li->bla[2*con] = a1;
- li->bla[2*con+1] = a2;
+ li->bllen[con] = lenA;
+ li->atoms[con].index1 = a1;
+ li->atoms[con].index2 = a2;
/* Make space in the constraint connection matrix for constraints
* connected to both end of the current constraint.
//! Sets matrix indices.
static void set_matrix_indices(Lincs *li,
- const Task *li_task,
+ const Task &li_task,
const t_blocka *at2con,
bool bSortMatrix)
{
- int b;
-
- for (b = li_task->b0; b < li_task->b1; b++)
+ for (int b = li_task.b0; b < li_task.b1; b++)
{
- int a1, a2, i, k;
+ const int a1 = li->atoms[b].index1;
+ const int a2 = li->atoms[b].index2;
- a1 = li->bla[b*2];
- a2 = li->bla[b*2 + 1];
-
- i = li->blnr[b];
- for (k = at2con->index[a1]; k < at2con->index[a1 + 1]; k++)
+ int i = li->blnr[b];
+ for (int k = at2con->index[a1]; k < at2con->index[a1 + 1]; k++)
{
- int concon;
-
- concon = li->con_index[at2con->a[k]];
+ int concon = li->con_index[at2con->a[k]];
if (concon != b)
{
li->blbnb[i++] = concon;
}
}
- for (k = at2con->index[a2]; k < at2con->index[a2 + 1]; k++)
+ for (int k = at2con->index[a2]; k < at2con->index[a2 + 1]; k++)
{
- int concon;
-
- concon = li->con_index[at2con->a[k]];
+ int concon = li->con_index[at2con->a[k]];
if (concon != b)
{
li->blbnb[i++] = concon;
li->con_index.resize(numEntries);
li->bllen0.resize(numEntries);
li->ddist.resize(numEntries);
- li->bla.resize(numEntries*2);
+ li->atoms.resize(numEntries);
li->blc.resize(numEntries);
li->blc1.resize(numEntries);
li->blnr.resize(numEntries);
* (e.g. because we are doing EM) we get imbalance, but since that doesn't
* happen during normal MD, that's ok.
*/
- int ncon_assign, ncon_target, con, th;
/* Determine the number of constraints we need to assign here */
- ncon_assign = ncon_tot;
+ int ncon_assign = ncon_tot;
if (!bDynamics)
{
/* With energy minimization, flexible constraints are ignored
/* Set the target constraint count per task to exactly uniform,
* this might be overridden below.
*/
- ncon_target = (ncon_assign + li->ntask - 1)/li->ntask;
+ int ncon_target = (ncon_assign + li->ntask - 1)/li->ntask;
/* Mark all constraints as unassigned by setting their index to -1 */
- for (con = 0; con < ncon_tot; con++)
+ for (int con = 0; con < ncon_tot; con++)
{
li->con_index[con] = -1;
}
- con = 0;
- for (th = 0; th < li->ntask; th++)
+ int con = 0;
+ for (int th = 0; th < li->ntask; th++)
{
Task *li_task;
last = li_task->b1 - 1;
for (i = li_task->b1; i < li->nc; i++)
{
- li->bla[i*2 ] = li->bla[last*2 ];
- li->bla[i*2 + 1] = li->bla[last*2 + 1];
+ li->atoms[i] = li->atoms[last];
li->bllen0[i] = li->bllen0[last];
li->ddist[i] = li->ddist[last];
li->bllen[i] = li->bllen[last];
li->blbnb.resize(li->ncc);
#pragma omp parallel for num_threads(li->ntask) schedule(static)
- for (th = 0; th < li->ntask; th++)
+ for (int th = 0; th < li->ntask; th++)
{
try
{
- Task *li_task;
-
- li_task = &li->task[th];
+ Task &li_task = li->task[th];
if (li->ncg_triangle > 0)
{
/* This is allocating too much, but it is difficult to improve */
- li_task->triangle.resize(li_task->b1 - li_task->b0);
- li_task->tri_bits.resize(li_task->b1 - li_task->b0);
+ li_task.triangle.resize(li_task.b1 - li_task.b0);
+ li_task.tri_bits.resize(li_task.b1 - li_task.b0);
}
set_matrix_indices(li, li_task, &at2con, bSortMatrix);
//! Issues a warning when LINCS constraints cannot be satisfied.
static void lincs_warning(gmx_domdec_t *dd, const rvec *x, rvec *xprime, t_pbc *pbc,
- int ncons, gmx::ArrayRef<const int> bla,
+ int ncons, gmx::ArrayRef<const AtomPair> atoms,
gmx::ArrayRef<const real> bllen, real wangle,
int maxwarn, int *warncount)
{
- int b, i, j;
- rvec v0, v1;
- real wfac, d0, d1, cosine;
-
- wfac = std::cos(DEG2RAD*wangle);
+ real wfac = std::cos(DEG2RAD*wangle);
fprintf(stderr,
"bonds that rotated more than %g degrees:\n"
" atom 1 atom 2 angle previous, current, constraint length\n",
wangle);
- for (b = 0; b < ncons; b++)
+ for (int b = 0; b < ncons; b++)
{
- i = bla[2*b];
- j = bla[2*b+1];
+ const int i = atoms[b].index1;
+ const int j = atoms[b].index2;
+ rvec v0;
+ rvec v1;
if (pbc)
{
pbc_dx_aiuc(pbc, x[i], x[j], v0);
rvec_sub(x[i], x[j], v0);
rvec_sub(xprime[i], xprime[j], v1);
}
- d0 = norm(v0);
- d1 = norm(v1);
- cosine = ::iprod(v0, v1)/(d0*d1);
+ real d0 = norm(v0);
+ real d1 = norm(v1);
+ real cosine = ::iprod(v0, v1)/(d0*d1);
if (cosine < wfac)
{
fprintf(stderr,
const rvec *x, const t_pbc *pbc,
real *ncons_loc, real *ssd, real *max, int *imax)
{
- gmx::ArrayRef<const int> bla = lincsd.bla;
- gmx::ArrayRef<const real> bllen = lincsd.bllen;
- gmx::ArrayRef<const int> nlocat = lincsd.nlocat;
+ gmx::ArrayRef<const AtomPair> atoms = lincsd.atoms;
+ gmx::ArrayRef<const real> bllen = lincsd.bllen;
+ gmx::ArrayRef<const int> nlocat = lincsd.nlocat;
real ma = 0;
real ssd2 = 0;
rvec dx;
if (pbc)
{
- pbc_dx_aiuc(pbc, x[bla[2*b]], x[bla[2*b+1]], dx);
+ pbc_dx_aiuc(pbc, x[atoms[b].index1], x[atoms[b].index2], dx);
}
else
{
- rvec_sub(x[bla[2*b]], x[bla[2*b+1]], dx);
+ rvec_sub(x[atoms[b].index1], x[atoms[b].index2], dx);
}
real r2 = ::norm2(dx);
real len = r2*gmx::invsqrt(r2);
t_nrnb *nrnb,
int maxwarn, int *warncount)
{
- gmx_bool bCalcDHDL;
- char buf2[22], buf3[STRLEN];
- int i, p_imax;
- real ncons_loc, p_ssd, p_max = 0;
- rvec dx;
- bool bOK, bWarn;
-
- bOK = TRUE;
+ bool bOK = TRUE;
/* This boolean should be set by a flag passed to this routine.
* We can also easily check if any constraint length is changed,
* if not dH/dlambda=0 and we can also set the boolean to FALSE.
*/
- bCalcDHDL = (ir.efep != efepNO && dvdlambda != nullptr);
+ bool bCalcDHDL = (ir.efep != efepNO && dvdlambda != nullptr);
if (lincsd->nc == 0 && cr->dd == nullptr)
{
set_lincs_matrix(lincsd, md.invmass, md.lambda);
}
- for (i = 0; i < lincsd->nc; i++)
+ for (int i = 0; i < lincsd->nc; i++)
{
lincsd->bllen[i] = lincsd->bllen0[i] + lambda*lincsd->ddist[i];
}
/* Set the flexible constraint lengths to the old lengths */
if (pbc != nullptr)
{
- for (i = 0; i < lincsd->nc; i++)
+ for (int i = 0; i < lincsd->nc; i++)
{
if (lincsd->bllen[i] == 0)
{
- pbc_dx_aiuc(pbc, x[lincsd->bla[2*i]], x[lincsd->bla[2*i+1]], dx);
+ rvec dx;
+ pbc_dx_aiuc(pbc, x[lincsd->atoms[i].index1], x[lincsd->atoms[i].index2], dx);
lincsd->bllen[i] = norm(dx);
}
}
}
else
{
- for (i = 0; i < lincsd->nc; i++)
+ for (int i = 0; i < lincsd->nc; i++)
{
if (lincsd->bllen[i] == 0)
{
lincsd->bllen[i] =
- std::sqrt(distance2(x[lincsd->bla[2*i]],
- x[lincsd->bla[2*i+1]]));
+ std::sqrt(distance2(x[lincsd->atoms[i].index1],
+ x[lincsd->atoms[i].index2]));
}
}
}
}
+ int p_imax;
+ real ncons_loc;
+ real p_ssd;
+ real p_max = 0;
if (debug)
{
cconerr(*lincsd, xprime, pbc,
* at the same time. But as we only need to detect
* if a warning occurred or not, this is not an issue.
*/
- bWarn = FALSE;
+ bool bWarn = FALSE;
/* The OpenMP parallel region of constrain_lincs for coords */
#pragma omp parallel num_threads(lincsd->ntask)
fprintf(debug, " Rel. Constraint Deviation: RMS MAX between atoms\n");
fprintf(debug, " Before LINCS %.6f %.6f %6d %6d\n",
std::sqrt(p_ssd/ncons_loc), p_max,
- ddglatnr(cr->dd, lincsd->bla[2*p_imax]),
- ddglatnr(cr->dd, lincsd->bla[2*p_imax+1]));
+ ddglatnr(cr->dd, lincsd->atoms[p_imax].index1),
+ ddglatnr(cr->dd, lincsd->atoms[p_imax].index2));
}
if (computeRmsd || debug)
{
fprintf(debug,
" After LINCS %.6f %.6f %6d %6d\n\n",
std::sqrt(p_ssd/ncons_loc), p_max,
- ddglatnr(cr->dd, lincsd->bla[2*p_imax]),
- ddglatnr(cr->dd, lincsd->bla[2*p_imax+1]));
+ ddglatnr(cr->dd, lincsd->atoms[p_imax].index1),
+ ddglatnr(cr->dd, lincsd->atoms[p_imax].index2));
}
if (bWarn)
{
cconerr(*lincsd, xprime, pbc,
&ncons_loc, &p_ssd, &p_max, &p_imax);
+ std::string simMesg;
if (isMultiSim(ms))
{
- sprintf(buf3, " in simulation %d", ms->sim);
- }
- else
- {
- buf3[0] = 0;
+ simMesg += gmx::formatString(" in simulation %d", ms->sim);
}
fprintf(stderr,
- "\nStep %s, time %g (ps) LINCS WARNING%s\n"
+ "\nStep %" PRId64 ", time %g (ps) LINCS WARNING%s\n"
"relative constraint deviation after LINCS:\n"
"rms %.6f, max %.6f (between atoms %d and %d)\n",
- gmx_step_str(step, buf2), ir.init_t+step*ir.delta_t,
- buf3,
+ step, ir.init_t+step*ir.delta_t,
+ simMesg.c_str(),
std::sqrt(p_ssd/ncons_loc), p_max,
- ddglatnr(cr->dd, lincsd->bla[2*p_imax]),
- ddglatnr(cr->dd, lincsd->bla[2*p_imax+1]));
+ ddglatnr(cr->dd, lincsd->atoms[p_imax].index1),
+ ddglatnr(cr->dd, lincsd->atoms[p_imax].index2));
lincs_warning(cr->dd, x, xprime, pbc,
- lincsd->nc, lincsd->bla, lincsd->bllen,
+ lincsd->nc, lincsd->atoms, lincsd->bllen,
ir.LincsWarnAngle, maxwarn, warncount);
}
bOK = (p_max < 0.5);
if (lincsd->ncg_flex)
{
- for (i = 0; (i < lincsd->nc); i++)
+ for (int i = 0; (i < lincsd->nc); i++)
{
if (lincsd->bllen0[i] == 0 && lincsd->ddist[i] == 0)
{
if (bCalcVir && lincsd->ntask > 1)
{
- for (i = 1; i < lincsd->ntask; i++)
+ for (int i = 1; i < lincsd->ntask; i++)
{
m_add(vir_r_m_dr, lincsd->task[i].vir_r_m_dr, vir_r_m_dr);
}
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff