* the research papers on the package. Check out http://www.gromacs.org.
*/
/* This file is completely threadsafe - keep it that way! */
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
+#include "gmxpre.h"
#include <math.h>
-#include "main.h"
-#include "constr.h"
-#include "copyrite.h"
-#include "physics.h"
-#include "vec.h"
-#include "pbc.h"
-#include "gromacs/utility/smalloc.h"
-#include "mdrun.h"
-#include "nrnb.h"
-#include "domdec.h"
-#include "mtop_util.h"
-#include "gmx_omp_nthreads.h"
+#include <stdlib.h>
+#include "gromacs/domdec/domdec.h"
#include "gromacs/fileio/gmxfio.h"
-#include "gmx_fatal.h"
+#include "gromacs/legacyheaders/constr.h"
+#include "gromacs/legacyheaders/copyrite.h"
+#include "gromacs/legacyheaders/gmx_omp_nthreads.h"
+#include "gromacs/legacyheaders/mdrun.h"
+#include "gromacs/legacyheaders/nrnb.h"
+#include "gromacs/legacyheaders/types/commrec.h"
+#include "gromacs/math/units.h"
+#include "gromacs/math/vec.h"
+#include "gromacs/mdlib/bitmask.h"
+#include "gromacs/pbcutil/pbc.h"
+#include "gromacs/topology/block.h"
+#include "gromacs/topology/mtop_util.h"
+#include "gromacs/utility/fatalerror.h"
#include "gromacs/utility/gmxomp.h"
+#include "gromacs/utility/smalloc.h"
typedef struct {
int b0; /* first constraint for this thread */
real *bllen; /* the reference bond length */
int nth; /* The number of threads doing LINCS */
lincs_thread_t *th; /* LINCS thread division */
- unsigned *atf; /* atom flags for thread parallelization */
+ gmx_bitmask_t *atf; /* atom flags for thread parallelization */
int atf_nalloc; /* allocation size of atf */
/* arrays for temporary storage in the LINCS algorithm */
rvec *tmpv;
gmx_bool bCalcVir, tensor rmdf)
{
int b0, b1, b, i, j, k, n;
- real tmp0, tmp1, tmp2, im1, im2, mvb, rlen, len, wfac, lam;
+ real tmp0, tmp1, tmp2, mvb;
rvec dx;
int *bla, *blnr, *blbnb;
rvec *r;
gmx_bool bCalcVir, tensor vir_r_m_dr)
{
int b0, b1, b, i, j, k, n, iter;
- real tmp0, tmp1, tmp2, im1, im2, mvb, rlen, len, len2, dlen2, wfac;
+ real tmp0, tmp1, tmp2, mvb, rlen, len, len2, dlen2, wfac;
rvec dx;
int *bla, *blnr, *blbnb;
rvec *r;
dlen2 = 2*len2 - norm2(dx);
if (dlen2 < wfac*len2 && (nlocat == NULL || nlocat[b]))
{
+ /* not race free - see detailed comment in caller */
*warn = b;
}
if (dlen2 > 0)
li->triangle[li->ntriangle] = i;
li->tri_bits[li->ntriangle] = 0;
li->ntriangle++;
- if (li->blnr[i+1] - li->blnr[i] > sizeof(li->tri_bits[0])*8 - 1)
+ if (li->blnr[i+1] - li->blnr[i] > static_cast<int>(sizeof(li->tri_bits[0])*8 - 1))
{
gmx_fatal(FARGS, "A constraint is connected to %d constraints, this is more than the %d allowed for constraints participating in triangles",
li->blnr[i+1] - li->blnr[i],
{
lincs_thread_t *li_m;
int th;
- unsigned *atf;
+ gmx_bitmask_t *atf;
int a;
if (natoms > li->atf_nalloc)
/* Clear the atom flags */
for (a = 0; a < natoms; a++)
{
- atf[a] = 0;
+ bitmask_clear(&atf[a]);
+ }
+
+ if (li->nth > BITMASK_SIZE)
+ {
+ gmx_fatal(FARGS, "More than %d threads is not supported for LINCS.", BITMASK_SIZE);
}
for (th = 0; th < li->nth; th++)
li_th->b0 = (li->nc* th )/li->nth;
li_th->b1 = (li->nc*(th+1))/li->nth;
- if (th < sizeof(*atf)*8)
+ /* For each atom set a flag for constraints from each */
+ for (b = li_th->b0; b < li_th->b1; b++)
{
- /* For each atom set a flag for constraints from each */
- for (b = li_th->b0; b < li_th->b1; b++)
- {
- atf[li->bla[b*2] ] |= (1U<<th);
- atf[li->bla[b*2+1]] |= (1U<<th);
- }
+ bitmask_set_bit(&atf[li->bla[b*2]], th);
+ bitmask_set_bit(&atf[li->bla[b*2+1]], th);
}
}
for (th = 0; th < li->nth; th++)
{
lincs_thread_t *li_th;
- unsigned mask;
+ gmx_bitmask_t mask;
int b;
li_th = &li->th[th];
srenew(li_th->ind_r, li_th->ind_nalloc);
}
- if (th < sizeof(*atf)*8)
- {
- mask = (1U<<th) - 1U;
+ bitmask_init_low_bits(&mask, th);
- li_th->nind = 0;
- li_th->nind_r = 0;
- for (b = li_th->b0; b < li_th->b1; b++)
+ li_th->nind = 0;
+ li_th->nind_r = 0;
+ for (b = li_th->b0; b < li_th->b1; b++)
+ {
+ /* 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))
{
- /* We let the constraint with the lowest thread index
- * operate on atoms with constraints from multiple threads.
- */
- if (((atf[li->bla[b*2]] & mask) == 0) &&
- ((atf[li->bla[b*2+1]] & mask) == 0))
- {
- /* Add the constraint to the local atom update index */
- li_th->ind[li_th->nind++] = b;
- }
- else
- {
- /* Add the constraint to the rest block */
- li_th->ind_r[li_th->nind_r++] = b;
- }
+ /* Add the constraint to the local atom update index */
+ li_th->ind[li_th->nind++] = b;
}
- }
- else
- {
- /* We are out of bits, assign all constraints to rest */
- for (b = li_th->b0; b < li_th->b1; b++)
+ else
{
+ /* Add the constraint to the rest block */
li_th->ind_r[li_th->nind_r++] = b;
}
}
int i, k, ncc_alloc, ni, con, nconnect, concon;
int type, a1, a2;
real lenA = 0, lenB;
- gmx_bool bLocal;
li->nc = 0;
li->ncc = 0;
li->blnr[con] = nconnect;
for (i = 0; i < ni; i++)
{
- bLocal = TRUE;
type = iatom[3*i];
a1 = iatom[3*i+1];
a2 = iatom[3*i+2];
{
gmx_bool bCalcDHDL;
char buf[STRLEN], buf2[22], buf3[STRLEN];
- int i, warn, p_imax, error;
+ int i, warn, p_imax;
real ncons_loc, p_ssd, p_max = 0;
rvec dx;
gmx_bool bOK;
if (econqCoord)
{
/* dhdlambda contains dH/dlambda*dt^2, correct for this */
+ /* TODO This should probably use invdt, so that sd integrator scaling works properly */
dhdlambda /= ir->delta_t*ir->delta_t;
}
*dvdlambda += dhdlambda;
biod.pnpi.spb.ru / alexxy / gromacs.git / blobdiff