src/gromacs/mdlib/clincs.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team.
   6  * Copyright (c) 2013,2014,2015,2016, by the GROMACS development team, led by
   7  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   8  * and including many others, as listed in the AUTHORS file in the
   9  * top-level source directory and at http://www.gromacs.org.
  10  *
  11  * GROMACS is free software; you can redistribute it and/or
  12  * modify it under the terms of the GNU Lesser General Public License
  13  * as published by the Free Software Foundation; either version 2.1
  14  * of the License, or (at your option) any later version.
  15  *
  16  * GROMACS is distributed in the hope that it will be useful,
  17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19  * Lesser General Public License for more details.
  20  *
  21  * You should have received a copy of the GNU Lesser General Public
  22  * License along with GROMACS; if not, see
  23  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  24  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  25  *
  26  * If you want to redistribute modifications to GROMACS, please
  27  * consider that scientific software is very special. Version
  28  * control is crucial - bugs must be traceable. We will be happy to
  29  * consider code for inclusion in the official distribution, but
  30  * derived work must not be called official GROMACS. Details are found
  31  * in the README & COPYING files - if they are missing, get the
  32  * official version at http://www.gromacs.org.
  33  *
  34  * To help us fund GROMACS development, we humbly ask that you cite
  35  * the research papers on the package. Check out http://www.gromacs.org.
  36  */
  37 /* This file is completely threadsafe - keep it that way! */
  38 #include "gmxpre.h"
  39
  40 #include "config.h"
  41
  42 #include <assert.h>
  43 #include <math.h>
  44 #include <stdlib.h>
  45
  46 #include <algorithm>
  47
  48 #include "gromacs/domdec/domdec.h"
  49 #include "gromacs/fileio/gmxfio.h"
  50 #include "gromacs/legacyheaders/constr.h"
  51 #include "gromacs/legacyheaders/copyrite.h"
  52 #include "gromacs/legacyheaders/gmx_omp_nthreads.h"
  53 #include "gromacs/legacyheaders/mdrun.h"
  54 #include "gromacs/legacyheaders/nrnb.h"
  55 #include "gromacs/legacyheaders/types/commrec.h"
  56 #include "gromacs/math/units.h"
  57 #include "gromacs/math/vec.h"
  58 #include "gromacs/pbcutil/pbc-simd.h"
  59 #include "gromacs/pbcutil/pbc.h"
  60 #include "gromacs/simd/simd.h"
  61 #include "gromacs/simd/simd_math.h"
  62 #include "gromacs/simd/vector_operations.h"
  63 #include "gromacs/topology/block.h"
  64 #include "gromacs/topology/mtop_util.h"
  65 #include "gromacs/utility/bitmask.h"
  66 #include "gromacs/utility/fatalerror.h"
  67 #include "gromacs/utility/gmxomp.h"
  68 #include "gromacs/utility/smalloc.h"
  69
  70 /* MSVC 2010 produces buggy SIMD PBC code, disable SIMD for MSVC <= 2010 */
  71 #if defined GMX_SIMD_HAVE_REAL && !(defined _MSC_VER && _MSC_VER < 1700) && !defined(__ICL)
  72 #define LINCS_SIMD
  73 #endif
  74
  75
  76 #if defined(GMX_SIMD_X86_AVX_256) || defined(GMX_SIMD_X86_AVX2_256)
  77
  78 // This was originally work-in-progress for augmenting the SIMD module with
  79 // masked load/store operations. Instead, that turned into and extended SIMD
  80 // interface that supports gather/scatter in all platforms, which will be
  81 // part of a future Gromacs version. However, since the code for bonded
  82 // interactions and LINCS was already written it would be a pity not to get
  83 // the performance gains in Gromacs-5.1. For this reason we have added it as
  84 // a bit of a hack in the two files that use it. It will be replaced with the
  85 // new generic functionality after version 5.1
  86
  87 #    ifdef GMX_DOUBLE
  88 static gmx_inline void gmx_simdcall
  89 gmx_hack_simd_transpose4_r(gmx_simd_double_t *row0,
  90                            gmx_simd_double_t *row1,
  91                            gmx_simd_double_t *row2,
  92                            gmx_simd_double_t *row3)
  93 {
  94     __m256d tmp0, tmp1, tmp2, tmp3;
  95
  96     tmp0  = _mm256_unpacklo_pd(*row0, *row1);
  97     tmp2  = _mm256_unpacklo_pd(*row2, *row3);
  98     tmp1  = _mm256_unpackhi_pd(*row0, *row1);
  99     tmp3  = _mm256_unpackhi_pd(*row2, *row3);
 100     *row0 = _mm256_permute2f128_pd(tmp0, tmp2, 0x20);
 101     *row1 = _mm256_permute2f128_pd(tmp1, tmp3, 0x20);
 102     *row2 = _mm256_permute2f128_pd(tmp0, tmp2, 0x31);
 103     *row3 = _mm256_permute2f128_pd(tmp1, tmp3, 0x31);
 104 }
 105
 106 static gmx_inline void gmx_simdcall
 107 gmx_hack_simd4_transpose_to_simd_r(const gmx_simd4_double_t *a,
 108                                    gmx_simd_double_t        *row0,
 109                                    gmx_simd_double_t        *row1,
 110                                    gmx_simd_double_t        *row2,
 111                                    gmx_simd_double_t        *row3)
 112 {
 113     *row0 = a[0];
 114     *row1 = a[1];
 115     *row2 = a[2];
 116     *row3 = a[3];
 117
 118     gmx_hack_simd_transpose4_r(row0, row1, row2, row3);
 119 }
 120
 121 static gmx_inline void gmx_simdcall
 122 gmx_hack_simd_transpose_to_simd4_r(gmx_simd_double_t   row0,
 123                                    gmx_simd_double_t   row1,
 124                                    gmx_simd_double_t   row2,
 125                                    gmx_simd_double_t   row3,
 126                                    gmx_simd4_double_t *a)
 127 {
 128     a[0] = row0;
 129     a[1] = row1;
 130     a[2] = row2;
 131     a[3] = row3;
 132
 133     gmx_hack_simd_transpose4_r(&a[0], &a[1], &a[2], &a[3]);
 134 }
 135
 136
 137 #    ifdef GMX_SIMD_X86_AVX_GCC_MASKLOAD_BUG
 138 #        define gmx_hack_simd4_load3_r(mem)      _mm256_maskload_pd((mem), _mm_castsi128_ps(_mm256_set_epi32(0, 0, -1, -1, -1, -1, -1, -1)))
 139 #        define gmx_hack_simd4_store3_r(mem, x)   _mm256_maskstore_pd((mem), _mm_castsi128_ps(_mm256_set_epi32(0, 0, -1, -1, -1, -1, -1, -1)), (x))
 140 #    else
 141 #        define gmx_hack_simd4_load3_r(mem)      _mm256_maskload_pd((mem), _mm256_set_epi32(0, 0, -1, -1, -1, -1, -1, -1))
 142 #        define gmx_hack_simd4_store3_r(mem, x)   _mm256_maskstore_pd((mem), _mm256_set_epi32(0, 0, -1, -1, -1, -1, -1, -1), (x))
 143 #    endif
 144
 145 #    else /* single instead of double */
 146 static gmx_inline void gmx_simdcall
 147 gmx_hack_simd_transpose4_r(gmx_simd_float_t *row0,
 148                            gmx_simd_float_t *row1,
 149                            gmx_simd_float_t *row2,
 150                            gmx_simd_float_t *row3)
 151 {
 152     __m256 tmp0, tmp1, tmp2, tmp3;
 153
 154     tmp0  = _mm256_unpacklo_ps(*row0, *row1);
 155     tmp2  = _mm256_unpacklo_ps(*row2, *row3);
 156     tmp1  = _mm256_unpackhi_ps(*row0, *row1);
 157     tmp3  = _mm256_unpackhi_ps(*row2, *row3);
 158     *row0 = _mm256_shuffle_ps(tmp0, tmp2, 0x44);
 159     *row1 = _mm256_shuffle_ps(tmp0, tmp2, 0xEE);
 160     *row2 = _mm256_shuffle_ps(tmp1, tmp3, 0x44);
 161     *row3 = _mm256_shuffle_ps(tmp1, tmp3, 0xEE);
 162 }
 163
 164 static gmx_inline void gmx_simdcall
 165 gmx_hack_simd4_transpose_to_simd_r(const gmx_simd4_float_t *a,
 166                                    gmx_simd_float_t        *row0,
 167                                    gmx_simd_float_t        *row1,
 168                                    gmx_simd_float_t        *row2,
 169                                    gmx_simd_float_t        *row3)
 170 {
 171     *row0 = _mm256_insertf128_ps(_mm256_castps128_ps256(a[0]), a[4], 1);
 172     *row1 = _mm256_insertf128_ps(_mm256_castps128_ps256(a[1]), a[5], 1);
 173     *row2 = _mm256_insertf128_ps(_mm256_castps128_ps256(a[2]), a[6], 1);
 174     *row3 = _mm256_insertf128_ps(_mm256_castps128_ps256(a[3]), a[7], 1);
 175
 176     gmx_hack_simd_transpose4_r(row0, row1, row2, row3);
 177 }
 178
 179 static gmx_inline void gmx_simdcall
 180 gmx_hack_simd_transpose_to_simd4_r(gmx_simd_float_t   row0,
 181                                    gmx_simd_float_t   row1,
 182                                    gmx_simd_float_t   row2,
 183                                    gmx_simd_float_t   row3,
 184                                    gmx_simd4_float_t *a)
 185 {
 186     gmx_hack_simd_transpose4_r(&row0, &row1, &row2, &row3);
 187
 188     a[0] = _mm256_extractf128_ps(row0, 0);
 189     a[1] = _mm256_extractf128_ps(row1, 0);
 190     a[2] = _mm256_extractf128_ps(row2, 0);
 191     a[3] = _mm256_extractf128_ps(row3, 0);
 192     a[4] = _mm256_extractf128_ps(row0, 1);
 193     a[5] = _mm256_extractf128_ps(row1, 1);
 194     a[6] = _mm256_extractf128_ps(row2, 1);
 195     a[7] = _mm256_extractf128_ps(row3, 1);
 196 }
 197 #ifdef GMX_SIMD_X86_AVX_GCC_MASKLOAD_BUG
 198 #        define gmx_hack_simd4_load3_r(mem)      _mm_maskload_ps((mem), _mm_castsi256_pd(_mm_set_epi32(0, -1, -1, -1)))
 199 #        define gmx_hack_simd4_store3_r(mem, x)   _mm_maskstore_ps((mem), _mm_castsi256_pd(_mm_set_epi32(0, -1, -1, -1)), (x))
 200 #else
 201 #        define gmx_hack_simd4_load3_r(mem)      _mm_maskload_ps((mem), _mm_set_epi32(0, -1, -1, -1))
 202 #        define gmx_hack_simd4_store3_r(mem, x)   _mm_maskstore_ps((mem), _mm_set_epi32(0, -1, -1, -1), (x))
 203 #endif
 204
 205 #endif /* double */
 206
 207 #endif /* AVX */
 208
 209 #ifdef GMX_SIMD_HAVE_REAL
 210 /*! \brief Store differences between indexed rvecs in SIMD registers.
 211  *
 212  * Returns SIMD register with the difference vectors:
 213  *     v[pair_index[i*2]] - v[pair_index[i*2 + 1]]
 214  *
 215  * \param[in]     v           Array of rvecs
 216  * \param[in]     pair_index  Index pairs for GMX_SIMD_REAL_WIDTH vector pairs
 217  * \param[in,out] buf         Aligned tmp buffer of size 3*GMX_SIMD_REAL_WIDTH
 218  * \param[out]    dx          SIMD register with x difference
 219  * \param[out]    dy          SIMD register with y difference
 220  * \param[out]    dz          SIMD register with z difference
 221  */
 222 static gmx_inline void gmx_simdcall
 223 gmx_hack_simd_gather_rvec_dist_pair_index(const rvec      *v,
 224                                           const int       *pair_index,
 225                                           real gmx_unused *buf,
 226                                           gmx_simd_real_t *dx,
 227                                           gmx_simd_real_t *dy,
 228                                           gmx_simd_real_t *dz)
 229 {
 230 #if defined(GMX_SIMD_X86_AVX_256) || defined(GMX_SIMD_X86_AVX2_256)
 231     int              i;
 232     gmx_simd4_real_t d[GMX_SIMD_REAL_WIDTH];
 233     gmx_simd_real_t  tmp;
 234
 235     for (i = 0; i < GMX_SIMD_REAL_WIDTH; i++)
 236     {
 237         d[i] = gmx_simd4_sub_r(gmx_hack_simd4_load3_r(&(v[pair_index[i*2 + 0]][0])),
 238                                gmx_hack_simd4_load3_r(&(v[pair_index[i*2 + 1]][0])));
 239     }
 240
 241     gmx_hack_simd4_transpose_to_simd_r(d, dx, dy, dz, &tmp);
 242 #else
 243 #if GMX_ALIGNMENT
 244     GMX_ALIGNED(real, GMX_SIMD_REAL_WIDTH) buf_aligned[3*GMX_SIMD_REAL_WIDTH];
 245 #else
 246     real* buf_aligned = buf;
 247 #endif
 248
 249     int i, m;
 250
 251     for (i = 0; i < GMX_SIMD_REAL_WIDTH; i++)
 252     {
 253         /* Store the distances packed and aligned */
 254         for (m = 0; m < DIM; m++)
 255         {
 256             buf_aligned[m*GMX_SIMD_REAL_WIDTH + i] =
 257                 v[pair_index[i*2]][m] - v[pair_index[i*2 + 1]][m];
 258         }
 259     }
 260     *dx = gmx_simd_load_r(buf_aligned + 0*GMX_SIMD_REAL_WIDTH);
 261     *dy = gmx_simd_load_r(buf_aligned + 1*GMX_SIMD_REAL_WIDTH);
 262     *dz = gmx_simd_load_r(buf_aligned + 2*GMX_SIMD_REAL_WIDTH);
 263 #endif
 264 }
 265
 266
 267 /*! \brief Stores SIMD vector into multiple rvecs.
 268  *
 269  * \param[in]     x           SIMD register with x-components of the vectors
 270  * \param[in]     y           SIMD register with y-components of the vectors
 271  * \param[in]     z           SIMD register with z-components of the vectors
 272  * \param[in,out] buf         Aligned tmp buffer of size 3*GMX_SIMD_REAL_WIDTH
 273  * \param[out]    v           Array of GMX_SIMD_REAL_WIDTH rvecs
 274  */
 275 static gmx_inline void gmx_simdcall
 276 gmx_simd_store_vec_to_rvec(gmx_simd_real_t  x,
 277                            gmx_simd_real_t  y,
 278                            gmx_simd_real_t  z,
 279                            real gmx_unused *buf,
 280                            rvec            *v)
 281 {
 282 #if defined(GMX_SIMD_X86_AVX_256) || defined(GMX_SIMD_X86_AVX2_256)
 283     int              i;
 284     gmx_simd4_real_t s4[GMX_SIMD_REAL_WIDTH];
 285     gmx_simd_real_t  zero = gmx_simd_setzero_r();
 286
 287     gmx_hack_simd_transpose_to_simd4_r(x, y, z, zero, s4);
 288
 289     for (i = 0; i < GMX_SIMD_REAL_WIDTH; i++)
 290     {
 291         gmx_hack_simd4_store3_r(v[i], s4[i]);
 292     }
 293 #else
 294 #if GMX_ALIGNMENT
 295     GMX_ALIGNED(real, GMX_SIMD_REAL_WIDTH) buf_aligned[3*GMX_SIMD_REAL_WIDTH];
 296 #else
 297     real* buf_aligned = buf;
 298 #endif
 299
 300     int i, m;
 301
 302     gmx_simd_store_r(buf_aligned + 0*GMX_SIMD_REAL_WIDTH, x);
 303     gmx_simd_store_r(buf_aligned + 1*GMX_SIMD_REAL_WIDTH, y);
 304     gmx_simd_store_r(buf_aligned + 2*GMX_SIMD_REAL_WIDTH, z);
 305
 306     for (i = 0; i < GMX_SIMD_REAL_WIDTH; i++)
 307     {
 308         for (m = 0; m < DIM; m++)
 309         {
 310             v[i][m] = buf_aligned[m*GMX_SIMD_REAL_WIDTH + i];
 311         }
 312     }
 313 #endif
 314 }
 315 #endif /* GMX_SIMD_HAVE_REAL */
 316
 317
 318 typedef struct {
 319     int    b0;         /* first constraint for this task */
 320     int    b1;         /* b1-1 is the last constraint for this task */
 321     int    ntriangle;  /* the number of constraints in triangles */
 322     int   *triangle;   /* the list of triangle constraints */
 323     int   *tri_bits;   /* the bits tell if the matrix element should be used */
 324     int    tri_alloc;  /* allocation size of triangle and tri_bits */
 325     int    nind;       /* number of indices */
 326     int   *ind;        /* constraint index for updating atom data */
 327     int    nind_r;     /* number of indices */
 328     int   *ind_r;      /* constraint index for updating atom data */
 329     int    ind_nalloc; /* allocation size of ind and ind_r */
 330     tensor vir_r_m_dr; /* temporary variable for virial calculation */
 331     real   dhdlambda;  /* temporary variable for lambda derivative */
 332     real  *simd_buf;   /* Aligned work array for SIMD */
 333 } lincs_task_t;
 334
 335 typedef struct gmx_lincsdata {
 336     int             ncg;          /* the global number of constraints */
 337     int             ncg_flex;     /* the global number of flexible constraints */
 338     int             ncg_triangle; /* the global number of constraints in triangles */
 339     int             nIter;        /* the number of iterations */
 340     int             nOrder;       /* the order of the matrix expansion */
 341     int             max_connect;  /* the maximum number of constrains connected to a single atom */
 342
 343     int             nc_real;      /* the number of real constraints */
 344     int             nc;           /* the number of constraints including padding for SIMD */
 345     int             nc_alloc;     /* the number we allocated memory for */
 346     int             ncc;          /* the number of constraint connections */
 347     int             ncc_alloc;    /* the number we allocated memory for */
 348     real            matlam;       /* the FE lambda value used for filling blc and blmf */
 349     int            *con_index;    /* mapping from topology to LINCS constraints */
 350     real           *bllen0;       /* the reference distance in topology A */
 351     real           *ddist;        /* the reference distance in top B - the r.d. in top A */
 352     int            *bla;          /* the atom pairs involved in the constraints */
 353     real           *blc;          /* 1/sqrt(invmass1 + invmass2) */
 354     real           *blc1;         /* as blc, but with all masses 1 */
 355     int            *blnr;         /* index into blbnb and blmf */
 356     int            *blbnb;        /* list of constraint connections */
 357     int             ntriangle;    /* the local number of constraints in triangles */
 358     int             ncc_triangle; /* the number of constraint connections in triangles */
 359     gmx_bool        bCommIter;    /* communicate before each LINCS interation */
 360     real           *blmf;         /* matrix of mass factors for constraint connections */
 361     real           *blmf1;        /* as blmf, but with all masses 1 */
 362     real           *bllen;        /* the reference bond length */
 363     int            *nlocat;       /* the local atom count per constraint, can be NULL */
 364
 365     int             ntask;        /* The number of tasks = #threads for LINCS */
 366     lincs_task_t   *task;         /* LINCS thread division */
 367     gmx_bitmask_t  *atf;          /* atom flags for thread parallelization */
 368     int             atf_nalloc;   /* allocation size of atf */
 369     gmx_bool        bTaskDep;     /* are the LINCS tasks interdependent? */
 370     gmx_bool        bTaskDepTri;  /* are there triangle constraints that cross task borders? */
 371     /* arrays for temporary storage in the LINCS algorithm */
 372     rvec           *tmpv;
 373     real           *tmpncc;
 374     real           *tmp1;
 375     real           *tmp2;
 376     real           *tmp3;
 377     real           *tmp4;
 378     real           *mlambda; /* the Lagrange multipliers * -1 */
 379     /* storage for the constraint RMS relative deviation output */
 380     real            rmsd_data[3];
 381 } t_gmx_lincsdata;
 382
 383 /* Define simd_width for memory allocation used for SIMD code */
 384 #ifdef LINCS_SIMD
 385 static const int simd_width = GMX_SIMD_REAL_WIDTH;
 386 #else
 387 static const int simd_width = 1;
 388 #endif
 389 /* We can't use small memory alignments on many systems, so use min 64 bytes*/
 390 static const int align_bytes = std::max<int>(64, simd_width*sizeof(real));
 391
 392 real *lincs_rmsd_data(struct gmx_lincsdata *lincsd)
 393 {
 394     return lincsd->rmsd_data;
 395 }
 396
 397 real lincs_rmsd(struct gmx_lincsdata *lincsd, gmx_bool bSD2)
 398 {
 399     if (lincsd->rmsd_data[0] > 0)
 400     {
 401         return sqrt(lincsd->rmsd_data[bSD2 ? 2 : 1]/lincsd->rmsd_data[0]);
 402     }
 403     else
 404     {
 405         return 0;
 406     }
 407 }
 408
 409 /* Do a set of nrec LINCS matrix multiplications.
 410  * This function will return with up to date thread-local
 411  * constraint data, without an OpenMP barrier.
 412  */
 413 static void lincs_matrix_expand(const struct gmx_lincsdata *lincsd,
 414                                 const lincs_task_t *li_task,
 415                                 const real *blcc,
 416                                 real *rhs1, real *rhs2, real *sol)
 417 {
 418     int        b0, b1, nrec, rec;
 419     const int *blnr  = lincsd->blnr;
 420     const int *blbnb = lincsd->blbnb;
 421
 422     b0   = li_task->b0;
 423     b1   = li_task->b1;
 424     nrec = lincsd->nOrder;
 425
 426     for (rec = 0; rec < nrec; rec++)
 427     {
 428         int b;
 429
 430         if (lincsd->bTaskDep)
 431         {
 432 #pragma omp barrier
 433         }
 434         for (b = b0; b < b1; b++)
 435         {
 436             real mvb;
 437             int  n;
 438
 439             mvb = 0;
 440             for (n = blnr[b]; n < blnr[b+1]; n++)
 441             {
 442                 mvb = mvb + blcc[n]*rhs1[blbnb[n]];
 443             }
 444             rhs2[b] = mvb;
 445             sol[b]  = sol[b] + mvb;
 446         }
 447
 448         real *swap;
 449
 450         swap = rhs1;
 451         rhs1 = rhs2;
 452         rhs2 = swap;
 453     } /* nrec*(ncons+2*nrtot) flops */
 454
 455     if (lincsd->ntriangle > 0)
 456     {
 457         /* Perform an extra nrec recursions for only the constraints
 458          * involved in rigid triangles.
 459          * In this way their accuracy should come close to those of the other
 460          * constraints, since traingles of constraints can produce eigenvalues
 461          * around 0.7, while the effective eigenvalue for bond constraints
 462          * is around 0.4 (and 0.7*0.7=0.5).
 463          */
 464
 465         if (lincsd->bTaskDep)
 466         {
 467             /* We need a barrier here, since other threads might still be
 468              * reading the contents of rhs1 and/o rhs2.
 469              * We could avoid this barrier by introducing two extra rhs
 470              * arrays for the triangle constraints only.
 471              */
 472 #pragma omp barrier
 473         }
 474
 475         /* Constraints involved in a triangle are ensured to be in the same
 476          * LINCS task. This means no barriers are required during the extra
 477          * iterations for the triangle constraints.
 478          */
 479         const int *triangle = li_task->triangle;
 480         const int *tri_bits = li_task->tri_bits;
 481
 482         for (rec = 0; rec < nrec; rec++)
 483         {
 484             int tb;
 485
 486             for (tb = 0; tb < li_task->ntriangle; tb++)
 487             {
 488                 int  b, bits, nr0, nr1, n;
 489                 real mvb;
 490
 491                 b    = triangle[tb];
 492                 bits = tri_bits[tb];
 493                 mvb  = 0;
 494                 nr0  = blnr[b];
 495                 nr1  = blnr[b+1];
 496                 for (n = nr0; n < nr1; n++)
 497                 {
 498                     if (bits & (1 << (n - nr0)))
 499                     {
 500                         mvb = mvb + blcc[n]*rhs1[blbnb[n]];
 501                     }
 502                 }
 503                 rhs2[b] = mvb;
 504                 sol[b]  = sol[b] + mvb;
 505             }
 506
 507             real *swap;
 508
 509             swap = rhs1;
 510             rhs1 = rhs2;
 511             rhs2 = swap;
 512         } /* nrec*(ntriangle + ncc_triangle*2) flops */
 513
 514         if (lincsd->bTaskDepTri)
 515         {
 516             /* The constraints triangles are decoupled from each other,
 517              * but constraints in one triangle cross thread task borders.
 518              * We could probably avoid this with more advanced setup code.
 519              */
 520 #pragma omp barrier
 521         }
 522     }
 523 }
 524
 525 static void lincs_update_atoms_noind(int ncons, const int *bla,
 526                                      real prefac,
 527                                      const real *fac, rvec *r,
 528                                      const real *invmass,
 529                                      rvec *x)
 530 {
 531     int  b, i, j;
 532     real mvb, im1, im2, tmp0, tmp1, tmp2;
 533
 534     if (invmass != NULL)
 535     {
 536         for (b = 0; b < ncons; b++)
 537         {
 538             i        = bla[2*b];
 539             j        = bla[2*b+1];
 540             mvb      = prefac*fac[b];
 541             im1      = invmass[i];
 542             im2      = invmass[j];
 543             tmp0     = r[b][0]*mvb;
 544             tmp1     = r[b][1]*mvb;
 545             tmp2     = r[b][2]*mvb;
 546             x[i][0] -= tmp0*im1;
 547             x[i][1] -= tmp1*im1;
 548             x[i][2] -= tmp2*im1;
 549             x[j][0] += tmp0*im2;
 550             x[j][1] += tmp1*im2;
 551             x[j][2] += tmp2*im2;
 552         } /* 16 ncons flops */
 553     }
 554     else
 555     {
 556         for (b = 0; b < ncons; b++)
 557         {
 558             i        = bla[2*b];
 559             j        = bla[2*b+1];
 560             mvb      = prefac*fac[b];
 561             tmp0     = r[b][0]*mvb;
 562             tmp1     = r[b][1]*mvb;
 563             tmp2     = r[b][2]*mvb;
 564             x[i][0] -= tmp0;
 565             x[i][1] -= tmp1;
 566             x[i][2] -= tmp2;
 567             x[j][0] += tmp0;
 568             x[j][1] += tmp1;
 569             x[j][2] += tmp2;
 570         }
 571     }
 572 }
 573
 574 static void lincs_update_atoms_ind(int ncons, const int *ind, const int *bla,
 575                                    real prefac,
 576                                    const real *fac, rvec *r,
 577                                    const real *invmass,
 578                                    rvec *x)
 579 {
 580     int  bi, b, i, j;
 581     real mvb, im1, im2, tmp0, tmp1, tmp2;
 582
 583     if (invmass != NULL)
 584     {
 585         for (bi = 0; bi < ncons; bi++)
 586         {
 587             b        = ind[bi];
 588             i        = bla[2*b];
 589             j        = bla[2*b+1];
 590             mvb      = prefac*fac[b];
 591             im1      = invmass[i];
 592             im2      = invmass[j];
 593             tmp0     = r[b][0]*mvb;
 594             tmp1     = r[b][1]*mvb;
 595             tmp2     = r[b][2]*mvb;
 596             x[i][0] -= tmp0*im1;
 597             x[i][1] -= tmp1*im1;
 598             x[i][2] -= tmp2*im1;
 599             x[j][0] += tmp0*im2;
 600             x[j][1] += tmp1*im2;
 601             x[j][2] += tmp2*im2;
 602         } /* 16 ncons flops */
 603     }
 604     else
 605     {
 606         for (bi = 0; bi < ncons; bi++)
 607         {
 608             b        = ind[bi];
 609             i        = bla[2*b];
 610             j        = bla[2*b+1];
 611             mvb      = prefac*fac[b];
 612             tmp0     = r[b][0]*mvb;
 613             tmp1     = r[b][1]*mvb;
 614             tmp2     = r[b][2]*mvb;
 615             x[i][0] -= tmp0;
 616             x[i][1] -= tmp1;
 617             x[i][2] -= tmp2;
 618             x[j][0] += tmp0;
 619             x[j][1] += tmp1;
 620             x[j][2] += tmp2;
 621         } /* 16 ncons flops */
 622     }
 623 }
 624
 625 static void lincs_update_atoms(struct gmx_lincsdata *li, int th,
 626                                real prefac,
 627                                const real *fac, rvec *r,
 628                                const real *invmass,
 629                                rvec *x)
 630 {
 631     if (li->ntask == 1)
 632     {
 633         /* Single thread, we simply update for all constraints */
 634         lincs_update_atoms_noind(li->nc_real,
 635                                  li->bla, prefac, fac, r, invmass, x);
 636     }
 637     else
 638     {
 639         /* Update the atom vector components for our thread local
 640          * constraints that only access our local atom range.
 641          * This can be done without a barrier.
 642          */
 643         lincs_update_atoms_ind(li->task[th].nind, li->task[th].ind,
 644                                li->bla, prefac, fac, r, invmass, x);
 645
 646         if (li->task[li->ntask].nind > 0)
 647         {
 648             /* Update the constraints that operate on atoms
 649              * in multiple thread atom blocks on the master thread.
 650              */
 651 #pragma omp barrier
 652 #pragma omp master
 653             {
 654                 lincs_update_atoms_ind(li->task[li->ntask].nind,
 655                                        li->task[li->ntask].ind,
 656                                        li->bla, prefac, fac, r, invmass, x);
 657             }
 658         }
 659     }
 660 }
 661
 662 #ifdef LINCS_SIMD
 663 /* Calculate the constraint distance vectors r to project on from x.
 664  * Determine the right-hand side of the matrix equation using quantity f.
 665  * This function only differs from calc_dr_x_xp_simd below in that
 666  * no constraint length is subtracted and no PBC is used for f.
 667  */
 668 static void gmx_simdcall
 669 calc_dr_x_f_simd(int                       b0,
 670                  int                       b1,
 671                  const int *               bla,
 672                  const rvec * gmx_restrict x,
 673                  const rvec * gmx_restrict f,
 674                  const real * gmx_restrict blc,
 675                  const pbc_simd_t *        pbc_simd,
 676                  real * gmx_restrict       vbuf1,
 677                  real * gmx_restrict       vbuf2,
 678                  rvec * gmx_restrict       r,
 679                  real * gmx_restrict       rhs,
 680                  real * gmx_restrict       sol)
 681 {
 682     int bs;
 683
 684     assert(b0 % GMX_SIMD_REAL_WIDTH == 0);
 685
 686     for (bs = b0; bs < b1; bs += GMX_SIMD_REAL_WIDTH)
 687     {
 688         gmx_simd_real_t rx_S, ry_S, rz_S, n2_S, il_S;
 689         gmx_simd_real_t fx_S, fy_S, fz_S, ip_S, rhs_S;
 690
 691         gmx_hack_simd_gather_rvec_dist_pair_index(x, bla + bs*2, vbuf1,
 692                                                   &rx_S, &ry_S, &rz_S);
 693
 694         pbc_correct_dx_simd(&rx_S, &ry_S, &rz_S, pbc_simd);
 695
 696         n2_S  = gmx_simd_norm2_r(rx_S, ry_S, rz_S);
 697         il_S  = gmx_simd_invsqrt_r(n2_S);
 698
 699         rx_S  = gmx_simd_mul_r(rx_S, il_S);
 700         ry_S  = gmx_simd_mul_r(ry_S, il_S);
 701         rz_S  = gmx_simd_mul_r(rz_S, il_S);
 702
 703         gmx_simd_store_vec_to_rvec(rx_S, ry_S, rz_S, vbuf1, r + bs);
 704
 705         gmx_hack_simd_gather_rvec_dist_pair_index(f, bla + bs*2, vbuf2,
 706                                                   &fx_S, &fy_S, &fz_S);
 707
 708         ip_S  = gmx_simd_iprod_r(rx_S, ry_S, rz_S,
 709                                  fx_S, fy_S, fz_S);
 710
 711         rhs_S = gmx_simd_mul_r(gmx_simd_load_r(blc + bs), ip_S);
 712
 713         gmx_simd_store_r(rhs + bs, rhs_S);
 714         gmx_simd_store_r(sol + bs, rhs_S);
 715     }
 716 }
 717 #endif /* LINCS_SIMD */
 718
 719 /* LINCS projection, works on derivatives of the coordinates */
 720 static void do_lincsp(rvec *x, rvec *f, rvec *fp, t_pbc *pbc,
 721                       struct gmx_lincsdata *lincsd, int th,
 722                       real *invmass,
 723                       int econq, gmx_bool bCalcDHDL,
 724                       gmx_bool bCalcVir, tensor rmdf)
 725 {
 726     int      b0, b1, b;
 727     int     *bla, *blnr, *blbnb;
 728     rvec    *r;
 729     real    *blc, *blmf, *blcc, *rhs1, *rhs2, *sol;
 730
 731     b0 = lincsd->task[th].b0;
 732     b1 = lincsd->task[th].b1;
 733
 734     bla    = lincsd->bla;
 735     r      = lincsd->tmpv;
 736     blnr   = lincsd->blnr;
 737     blbnb  = lincsd->blbnb;
 738     if (econq != econqForce)
 739     {
 740         /* Use mass-weighted parameters */
 741         blc  = lincsd->blc;
 742         blmf = lincsd->blmf;
 743     }
 744     else
 745     {
 746         /* Use non mass-weighted parameters */
 747         blc  = lincsd->blc1;
 748         blmf = lincsd->blmf1;
 749     }
 750     blcc   = lincsd->tmpncc;
 751     rhs1   = lincsd->tmp1;
 752     rhs2   = lincsd->tmp2;
 753     sol    = lincsd->tmp3;
 754
 755 #ifdef LINCS_SIMD
 756
 757     /* This SIMD code does the same as the plain-C code after the #else.
 758      * The only difference is that we always call pbc code, as with SIMD
 759      * the overhead of pbc computation (when not needed) is small.
 760      */
 761     pbc_simd_t pbc_simd;
 762
 763     /* Convert the pbc struct for SIMD */
 764     set_pbc_simd(pbc, &pbc_simd);
 765
 766     /* Compute normalized x i-j vectors, store in r.
 767      * Compute the inner product of r and xp i-j and store in rhs1.
 768      */
 769     calc_dr_x_f_simd(b0, b1, bla, x, f, blc,
 770                      &pbc_simd,
 771                      lincsd->task[th].simd_buf,
 772                      lincsd->task[th].simd_buf + GMX_SIMD_REAL_WIDTH*DIM,
 773                      r, rhs1, sol);
 774
 775 #else /* LINCS_SIMD */
 776
 777     /* Compute normalized i-j vectors */
 778     if (pbc)
 779     {
 780         for (b = b0; b < b1; b++)
 781         {
 782             rvec dx;
 783
 784             pbc_dx_aiuc(pbc, x[bla[2*b]], x[bla[2*b+1]], dx);
 785             unitv(dx, r[b]);
 786         }
 787     }
 788     else
 789     {
 790         for (b = b0; b < b1; b++)
 791         {
 792             rvec dx;
 793
 794             rvec_sub(x[bla[2*b]], x[bla[2*b+1]], dx);
 795             unitv(dx, r[b]);
 796         } /* 16 ncons flops */
 797     }
 798
 799     for (b = b0; b < b1; b++)
 800     {
 801         int  i, j;
 802         real mvb;
 803
 804         i       = bla[2*b];
 805         j       = bla[2*b+1];
 806         mvb     = blc[b]*(r[b][0]*(f[i][0] - f[j][0]) +
 807                           r[b][1]*(f[i][1] - f[j][1]) +
 808                           r[b][2]*(f[i][2] - f[j][2]));
 809         rhs1[b] = mvb;
 810         sol[b]  = mvb;
 811         /* 7 flops */
 812     }
 813
 814 #endif /* LINCS_SIMD */
 815
 816     if (lincsd->bTaskDep)
 817     {
 818         /* We need a barrier, since the matrix construction below
 819          * can access entries in r of other threads.
 820          */
 821 #pragma omp barrier
 822     }
 823
 824     /* Construct the (sparse) LINCS matrix */
 825     for (b = b0; b < b1; b++)
 826     {
 827         int n;
 828
 829         for (n = blnr[b]; n < blnr[b+1]; n++)
 830         {
 831             blcc[n] = blmf[n]*iprod(r[b], r[blbnb[n]]);
 832         } /* 6 nr flops */
 833     }
 834     /* Together: 23*ncons + 6*nrtot flops */
 835
 836     lincs_matrix_expand(lincsd, &lincsd->task[th], blcc, rhs1, rhs2, sol);
 837     /* nrec*(ncons+2*nrtot) flops */
 838
 839     if (econq == econqDeriv_FlexCon)
 840     {
 841         /* We only want to constraint the flexible constraints,
 842          * so we mask out the normal ones by setting sol to 0.
 843          */
 844         for (b = b0; b < b1; b++)
 845         {
 846             if (!(lincsd->bllen0[b] == 0 && lincsd->ddist[b] == 0))
 847             {
 848                 sol[b] = 0;
 849             }
 850         }
 851     }
 852
 853     /* We multiply sol by blc, so we can use lincs_update_atoms for OpenMP */
 854     for (b = b0; b < b1; b++)
 855     {
 856         sol[b] *= blc[b];
 857     }
 858
 859     /* When constraining forces, we should not use mass weighting,
 860      * so we pass invmass=NULL, which results in the use of 1 for all atoms.
 861      */
 862     lincs_update_atoms(lincsd, th, 1.0, sol, r,
 863                        (econq != econqForce) ? invmass : NULL, fp);
 864
 865     if (bCalcDHDL)
 866     {
 867         real dhdlambda;
 868
 869         dhdlambda = 0;
 870         for (b = b0; b < b1; b++)
 871         {
 872             dhdlambda -= sol[b]*lincsd->ddist[b];
 873         }
 874
 875         lincsd->task[th].dhdlambda = dhdlambda;
 876     }
 877
 878     if (bCalcVir)
 879     {
 880         /* Constraint virial,
 881          * determines sum r_bond x delta f,
 882          * where delta f is the constraint correction
 883          * of the quantity that is being constrained.
 884          */
 885         for (b = b0; b < b1; b++)
 886         {
 887             real mvb, tmp1;
 888             int  i, j;
 889
 890             mvb = lincsd->bllen[b]*sol[b];
 891             for (i = 0; i < DIM; i++)
 892             {
 893                 tmp1 = mvb*r[b][i];
 894                 for (j = 0; j < DIM; j++)
 895                 {
 896                     rmdf[i][j] += tmp1*r[b][j];
 897                 }
 898             }
 899         } /* 23 ncons flops */
 900     }
 901 }
 902
 903 #ifdef LINCS_SIMD
 904 /* Calculate the constraint distance vectors r to project on from x.
 905  * Determine the right-hand side of the matrix equation using coordinates xp.
 906  */
 907 static void gmx_simdcall
 908 calc_dr_x_xp_simd(int                       b0,
 909                   int                       b1,
 910                   const int *               bla,
 911                   const rvec * gmx_restrict x,
 912                   const rvec * gmx_restrict xp,
 913                   const real * gmx_restrict bllen,
 914                   const real * gmx_restrict blc,
 915                   const pbc_simd_t *        pbc_simd,
 916                   real * gmx_restrict       vbuf1,
 917                   real * gmx_restrict       vbuf2,
 918                   rvec * gmx_restrict       r,
 919                   real * gmx_restrict       rhs,
 920                   real * gmx_restrict       sol)
 921 {
 922     int bs;
 923
 924     assert(b0 % GMX_SIMD_REAL_WIDTH == 0);
 925
 926     for (bs = b0; bs < b1; bs += GMX_SIMD_REAL_WIDTH)
 927     {
 928         gmx_simd_real_t rx_S, ry_S, rz_S, n2_S, il_S;
 929         gmx_simd_real_t rxp_S, ryp_S, rzp_S, ip_S, rhs_S;
 930
 931         gmx_hack_simd_gather_rvec_dist_pair_index(x, bla + bs*2, vbuf1,
 932                                                   &rx_S, &ry_S, &rz_S);
 933
 934         pbc_correct_dx_simd(&rx_S, &ry_S, &rz_S, pbc_simd);
 935
 936         n2_S  = gmx_simd_norm2_r(rx_S, ry_S, rz_S);
 937         il_S  = gmx_simd_invsqrt_r(n2_S);
 938
 939         rx_S  = gmx_simd_mul_r(rx_S, il_S);
 940         ry_S  = gmx_simd_mul_r(ry_S, il_S);
 941         rz_S  = gmx_simd_mul_r(rz_S, il_S);
 942
 943         gmx_simd_store_vec_to_rvec(rx_S, ry_S, rz_S, vbuf1, r + bs);
 944
 945         gmx_hack_simd_gather_rvec_dist_pair_index(xp, bla + bs*2, vbuf2,
 946                                                   &rxp_S, &ryp_S, &rzp_S);
 947
 948         pbc_correct_dx_simd(&rxp_S, &ryp_S, &rzp_S, pbc_simd);
 949
 950         ip_S  = gmx_simd_iprod_r(rx_S, ry_S, rz_S,
 951                                  rxp_S, ryp_S, rzp_S);
 952
 953         rhs_S = gmx_simd_mul_r(gmx_simd_load_r(blc + bs),
 954                                gmx_simd_sub_r(ip_S, gmx_simd_load_r(bllen + bs)));
 955
 956         gmx_simd_store_r(rhs + bs, rhs_S);
 957         gmx_simd_store_r(sol + bs, rhs_S);
 958     }
 959 }
 960 #endif /* LINCS_SIMD */
 961
 962 /* Determine the distances and right-hand side for the next iteration */
 963 static void calc_dist_iter(int                       b0,
 964                            int                       b1,
 965                            const int *               bla,
 966                            const rvec * gmx_restrict xp,
 967                            const real * gmx_restrict bllen,
 968                            const real * gmx_restrict blc,
 969                            const t_pbc *             pbc,
 970                            real                      wfac,
 971                            real * gmx_restrict       rhs,
 972                            real * gmx_restrict       sol,
 973                            gmx_bool *                bWarn)
 974 {
 975     int b;
 976
 977     for (b = b0; b < b1; b++)
 978     {
 979         real len, len2, dlen2, mvb;
 980         rvec dx;
 981
 982         len = bllen[b];
 983         if (pbc)
 984         {
 985             pbc_dx_aiuc(pbc, xp[bla[2*b]], xp[bla[2*b+1]], dx);
 986         }
 987         else
 988         {
 989             rvec_sub(xp[bla[2*b]], xp[bla[2*b+1]], dx);
 990         }
 991         len2  = len*len;
 992         dlen2 = 2*len2 - norm2(dx);
 993         if (dlen2 < wfac*len2)
 994         {
 995             /* not race free - see detailed comment in caller */
 996             *bWarn = TRUE;
 997         }
 998         if (dlen2 > 0)
 999         {
1000             mvb = blc[b]*(len - dlen2*gmx_invsqrt(dlen2));
1001         }
1002         else
1003         {
1004             mvb = blc[b]*len;
1005         }
1006         rhs[b]  = mvb;
1007         sol[b]  = mvb;
1008     } /* 20*ncons flops */
1009 }
1010
1011 #ifdef LINCS_SIMD
1012 /* As the function above, but using SIMD intrinsics */
1013 static void gmx_simdcall
1014 calc_dist_iter_simd(int                       b0,
1015                     int                       b1,
1016                     const int *               bla,
1017                     const rvec * gmx_restrict x,
1018                     const real * gmx_restrict bllen,
1019                     const real * gmx_restrict blc,
1020                     const pbc_simd_t *        pbc_simd,
1021                     real                      wfac,
1022                     real * gmx_restrict       vbuf,
1023                     real * gmx_restrict       rhs,
1024                     real * gmx_restrict       sol,
1025                     gmx_bool *                bWarn)
1026 {
1027     gmx_simd_real_t min_S  = gmx_simd_set1_r(GMX_REAL_MIN);
1028     gmx_simd_real_t two_S  = gmx_simd_set1_r(2.0);
1029     gmx_simd_real_t wfac_S = gmx_simd_set1_r(wfac);
1030     gmx_simd_bool_t warn_S;
1031
1032     int             bs;
1033
1034     assert(b0 % GMX_SIMD_REAL_WIDTH == 0);
1035
1036     /* Initialize all to FALSE */
1037     warn_S = gmx_simd_cmplt_r(two_S, gmx_simd_setzero_r());
1038
1039     for (bs = b0; bs < b1; bs += GMX_SIMD_REAL_WIDTH)
1040     {
1041         gmx_simd_real_t rx_S, ry_S, rz_S, n2_S;
1042         gmx_simd_real_t len_S, len2_S, dlen2_S, lc_S, blc_S;
1043
1044         gmx_hack_simd_gather_rvec_dist_pair_index(x, bla + bs*2, vbuf,
1045                                                   &rx_S, &ry_S, &rz_S);
1046
1047         pbc_correct_dx_simd(&rx_S, &ry_S, &rz_S, pbc_simd);
1048
1049         n2_S    = gmx_simd_norm2_r(rx_S, ry_S, rz_S);
1050
1051         len_S   = gmx_simd_load_r(bllen + bs);
1052         len2_S  = gmx_simd_mul_r(len_S, len_S);
1053
1054         dlen2_S = gmx_simd_fmsub_r(two_S, len2_S, n2_S);
1055
1056         warn_S  = gmx_simd_or_b(warn_S,
1057                                 gmx_simd_cmplt_r(dlen2_S,
1058                                                  gmx_simd_mul_r(wfac_S, len2_S)));
1059
1060         /* Avoid 1/0 by taking the max with REAL_MIN.
1061          * Note: when dlen2 is close to zero (90 degree constraint rotation),
1062          * the accuracy of the algorithm is no longer relevant.
1063          */
1064         dlen2_S = gmx_simd_max_r(dlen2_S, min_S);
1065
1066         lc_S    = gmx_simd_fnmadd_r(dlen2_S, gmx_simd_invsqrt_r(dlen2_S), len_S);
1067
1068         blc_S   = gmx_simd_load_r(blc + bs);
1069
1070         lc_S    = gmx_simd_mul_r(blc_S, lc_S);
1071
1072         gmx_simd_store_r(rhs + bs, lc_S);
1073         gmx_simd_store_r(sol + bs, lc_S);
1074     }
1075
1076     if (gmx_simd_anytrue_b(warn_S))
1077     {
1078         *bWarn = TRUE;
1079     }
1080 }
1081 #endif /* LINCS_SIMD */
1082
1083 static void do_lincs(rvec *x, rvec *xp, matrix box, t_pbc *pbc,
1084                      struct gmx_lincsdata *lincsd, int th,
1085                      const real *invmass,
1086                      t_commrec *cr,
1087                      gmx_bool bCalcDHDL,
1088                      real wangle, gmx_bool *bWarn,
1089                      real invdt, rvec * gmx_restrict v,
1090                      gmx_bool bCalcVir, tensor vir_r_m_dr)
1091 {
1092     int      b0, b1, b, i, j, n, iter;
1093     int     *bla, *blnr, *blbnb;
1094     rvec    *r;
1095     real    *blc, *blmf, *bllen, *blcc, *rhs1, *rhs2, *sol, *blc_sol, *mlambda;
1096     int     *nlocat;
1097
1098     b0 = lincsd->task[th].b0;
1099     b1 = lincsd->task[th].b1;
1100
1101     bla     = lincsd->bla;
1102     r       = lincsd->tmpv;
1103     blnr    = lincsd->blnr;
1104     blbnb   = lincsd->blbnb;
1105     blc     = lincsd->blc;
1106     blmf    = lincsd->blmf;
1107     bllen   = lincsd->bllen;
1108     blcc    = lincsd->tmpncc;
1109     rhs1    = lincsd->tmp1;
1110     rhs2    = lincsd->tmp2;
1111     sol     = lincsd->tmp3;
1112     blc_sol = lincsd->tmp4;
1113     mlambda = lincsd->mlambda;
1114     nlocat  = lincsd->nlocat;
1115
1116 #ifdef LINCS_SIMD
1117
1118     /* This SIMD code does the same as the plain-C code after the #else.
1119      * The only difference is that we always call pbc code, as with SIMD
1120      * the overhead of pbc computation (when not needed) is small.
1121      */
1122     pbc_simd_t pbc_simd;
1123
1124     /* Convert the pbc struct for SIMD */
1125     set_pbc_simd(pbc, &pbc_simd);
1126
1127     /* Compute normalized x i-j vectors, store in r.
1128      * Compute the inner product of r and xp i-j and store in rhs1.
1129      */
1130     calc_dr_x_xp_simd(b0, b1, bla, x, xp, bllen, blc,
1131                       &pbc_simd,
1132                       lincsd->task[th].simd_buf,
1133                       lincsd->task[th].simd_buf + GMX_SIMD_REAL_WIDTH*DIM,
1134                       r, rhs1, sol);
1135
1136 #else /* LINCS_SIMD */
1137
1138     if (pbc)
1139     {
1140         /* Compute normalized i-j vectors */
1141         for (b = b0; b < b1; b++)
1142         {
1143             rvec dx;
1144             real mvb;
1145
1146             pbc_dx_aiuc(pbc, x[bla[2*b]], x[bla[2*b+1]], dx);
1147             unitv(dx, r[b]);
1148
1149             pbc_dx_aiuc(pbc, xp[bla[2*b]], xp[bla[2*b+1]], dx);
1150             mvb     = blc[b]*(iprod(r[b], dx) - bllen[b]);
1151             rhs1[b] = mvb;
1152             sol[b]  = mvb;
1153         }
1154     }
1155     else
1156     {
1157         /* Compute normalized i-j vectors */
1158         for (b = b0; b < b1; b++)
1159         {
1160             real tmp0, tmp1, tmp2, rlen, mvb;
1161
1162             i       = bla[2*b];
1163             j       = bla[2*b+1];
1164             tmp0    = x[i][0] - x[j][0];
1165             tmp1    = x[i][1] - x[j][1];
1166             tmp2    = x[i][2] - x[j][2];
1167             rlen    = gmx_invsqrt(tmp0*tmp0 + tmp1*tmp1 + tmp2*tmp2);
1168             r[b][0] = rlen*tmp0;
1169             r[b][1] = rlen*tmp1;
1170             r[b][2] = rlen*tmp2;
1171             /* 16 ncons flops */
1172
1173             i       = bla[2*b];
1174             j       = bla[2*b+1];
1175             mvb     = blc[b]*(r[b][0]*(xp[i][0] - xp[j][0]) +
1176                               r[b][1]*(xp[i][1] - xp[j][1]) +
1177                               r[b][2]*(xp[i][2] - xp[j][2]) - bllen[b]);
1178             rhs1[b] = mvb;
1179             sol[b]  = mvb;
1180             /* 10 flops */
1181         }
1182         /* Together: 26*ncons + 6*nrtot flops */
1183     }
1184
1185 #endif /* LINCS_SIMD */
1186
1187     if (lincsd->bTaskDep)
1188     {
1189         /* We need a barrier, since the matrix construction below
1190          * can access entries in r of other threads.
1191          */
1192 #pragma omp barrier
1193     }
1194
1195     /* Construct the (sparse) LINCS matrix */
1196     for (b = b0; b < b1; b++)
1197     {
1198         for (n = blnr[b]; n < blnr[b+1]; n++)
1199         {
1200             blcc[n] = blmf[n]*iprod(r[b], r[blbnb[n]]);
1201         }
1202     }
1203     /* Together: 26*ncons + 6*nrtot flops */
1204
1205     lincs_matrix_expand(lincsd, &lincsd->task[th], blcc, rhs1, rhs2, sol);
1206     /* nrec*(ncons+2*nrtot) flops */
1207
1208 #ifndef LINCS_SIMD
1209     for (b = b0; b < b1; b++)
1210     {
1211         mlambda[b] = blc[b]*sol[b];
1212     }
1213 #else
1214     for (b = b0; b < b1; b += GMX_SIMD_REAL_WIDTH)
1215     {
1216         gmx_simd_store_r(mlambda + b,
1217                          gmx_simd_mul_r(gmx_simd_load_r(blc + b),
1218                                         gmx_simd_load_r(sol + b)));
1219     }
1220 #endif
1221
1222     /* Update the coordinates */
1223     lincs_update_atoms(lincsd, th, 1.0, mlambda, r, invmass, xp);
1224
1225     /*
1226      ********  Correction for centripetal effects  ********
1227      */
1228
1229     real wfac;
1230
1231     wfac = cos(DEG2RAD*wangle);
1232     wfac = wfac*wfac;
1233
1234     for (iter = 0; iter < lincsd->nIter; iter++)
1235     {
1236         if ((lincsd->bCommIter && DOMAINDECOMP(cr) && cr->dd->constraints))
1237         {
1238 #pragma omp barrier
1239 #pragma omp master
1240             {
1241                 /* Communicate the corrected non-local coordinates */
1242                 if (DOMAINDECOMP(cr))
1243                 {
1244                     dd_move_x_constraints(cr->dd, box, xp, NULL, FALSE);
1245                 }
1246             }
1247 #pragma omp barrier
1248         }
1249         else if (lincsd->bTaskDep)
1250         {
1251 #pragma omp barrier
1252         }
1253
1254 #ifdef LINCS_SIMD
1255         calc_dist_iter_simd(b0, b1, bla, xp, bllen, blc, &pbc_simd, wfac,
1256                             lincsd->task[th].simd_buf, rhs1, sol, bWarn);
1257 #else
1258         calc_dist_iter(b0, b1, bla, xp, bllen, blc, pbc, wfac,
1259                        rhs1, sol, bWarn);
1260         /* 20*ncons flops */
1261 #endif
1262
1263         lincs_matrix_expand(lincsd, &lincsd->task[th], blcc, rhs1, rhs2, sol);
1264         /* nrec*(ncons+2*nrtot) flops */
1265
1266 #ifndef LINCS_SIMD
1267         for (b = b0; b < b1; b++)
1268         {
1269             real mvb;
1270
1271             mvb         = blc[b]*sol[b];
1272             blc_sol[b]  = mvb;
1273             mlambda[b] += mvb;
1274         }
1275 #else
1276         for (b = b0; b < b1; b += GMX_SIMD_REAL_WIDTH)
1277         {
1278             gmx_simd_real_t mvb;
1279
1280             mvb = gmx_simd_mul_r(gmx_simd_load_r(blc + b),
1281                                  gmx_simd_load_r(sol + b));
1282             gmx_simd_store_r(blc_sol + b, mvb);
1283             gmx_simd_store_r(mlambda + b,
1284                              gmx_simd_add_r(gmx_simd_load_r(mlambda + b), mvb));
1285         }
1286 #endif
1287
1288         /* Update the coordinates */
1289         lincs_update_atoms(lincsd, th, 1.0, blc_sol, r, invmass, xp);
1290     }
1291     /* nit*ncons*(37+9*nrec) flops */
1292
1293     if (v != NULL)
1294     {
1295         /* Update the velocities */
1296         lincs_update_atoms(lincsd, th, invdt, mlambda, r, invmass, v);
1297         /* 16 ncons flops */
1298     }
1299
1300     if (nlocat != NULL && (bCalcDHDL || bCalcVir))
1301     {
1302         if (lincsd->bTaskDep)
1303         {
1304             /* In lincs_update_atoms threads might cross-read mlambda */
1305 #pragma omp barrier
1306         }
1307
1308         /* Only account for local atoms */
1309         for (b = b0; b < b1; b++)
1310         {
1311             mlambda[b] *= 0.5*nlocat[b];
1312         }
1313     }
1314
1315     if (bCalcDHDL)
1316     {
1317         real dhdl;
1318
1319         dhdl = 0;
1320         for (b = b0; b < b1; b++)
1321         {
1322             /* Note that this this is dhdl*dt^2, the dt^2 factor is corrected
1323              * later after the contributions are reduced over the threads.
1324              */
1325             dhdl -= lincsd->mlambda[b]*lincsd->ddist[b];
1326         }
1327         lincsd->task[th].dhdlambda = dhdl;
1328     }
1329
1330     if (bCalcVir)
1331     {
1332         /* Constraint virial */
1333         for (b = b0; b < b1; b++)
1334         {
1335             real tmp0, tmp1;
1336
1337             tmp0 = -bllen[b]*mlambda[b];
1338             for (i = 0; i < DIM; i++)
1339             {
1340                 tmp1 = tmp0*r[b][i];
1341                 for (j = 0; j < DIM; j++)
1342                 {
1343                     vir_r_m_dr[i][j] -= tmp1*r[b][j];
1344                 }
1345             }
1346         } /* 22 ncons flops */
1347     }
1348
1349     /* Total:
1350      * 26*ncons + 6*nrtot + nrec*(ncons+2*nrtot)
1351      * + nit * (20*ncons + nrec*(ncons+2*nrtot) + 17 ncons)
1352      *
1353      * (26+nrec)*ncons + (6+2*nrec)*nrtot
1354      * + nit * ((37+nrec)*ncons + 2*nrec*nrtot)
1355      * if nit=1
1356      * (63+nrec)*ncons + (6+4*nrec)*nrtot
1357      */
1358 }
1359
1360 /* Sets the elements in the LINCS matrix for task li_task */
1361 static void set_lincs_matrix_task(struct gmx_lincsdata *li,
1362                                   lincs_task_t         *li_task,
1363                                   const real           *invmass,
1364                                   int                  *ncc_triangle,
1365                                   int                  *nCrossTaskTriangles)
1366 {
1367     int        i;
1368
1369     /* Construct the coupling coefficient matrix blmf */
1370     li_task->ntriangle   = 0;
1371     *ncc_triangle        = 0;
1372     *nCrossTaskTriangles = 0;
1373     for (i = li_task->b0; i < li_task->b1; i++)
1374     {
1375         int a1, a2, n;
1376
1377         a1 = li->bla[2*i];
1378         a2 = li->bla[2*i+1];
1379         for (n = li->blnr[i]; (n < li->blnr[i+1]); n++)
1380         {
1381             int k, sign, center, end;
1382
1383             k = li->blbnb[n];
1384
1385             /* If we are using multiple, independent tasks for LINCS,
1386              * the calls to check_assign_connected should have
1387              * put all connected constraints in our task.
1388              */
1389             assert(li->bTaskDep || (k >= li_task->b0 && k < li_task->b1));
1390
1391             if (a1 == li->bla[2*k] || a2 == li->bla[2*k+1])
1392             {
1393                 sign = -1;
1394             }
1395             else
1396             {
1397                 sign = 1;
1398             }
1399             if (a1 == li->bla[2*k] || a1 == li->bla[2*k+1])
1400             {
1401                 center = a1;
1402                 end    = a2;
1403             }
1404             else
1405             {
1406                 center = a2;
1407                 end    = a1;
1408             }
1409             li->blmf[n]  = sign*invmass[center]*li->blc[i]*li->blc[k];
1410             li->blmf1[n] = sign*0.5;
1411             if (li->ncg_triangle > 0)
1412             {
1413                 int nk, kk;
1414
1415                 /* Look for constraint triangles */
1416                 for (nk = li->blnr[k]; (nk < li->blnr[k+1]); nk++)
1417                 {
1418                     kk = li->blbnb[nk];
1419                     if (kk != i && kk != k &&
1420                         (li->bla[2*kk] == end || li->bla[2*kk+1] == end))
1421                     {
1422                         /* Check if the constraints in this triangle actually
1423                          * belong to a different task. We still assign them
1424                          * here, since it's convenient for the triangle
1425                          * iterations, but we then need an extra barrier.
1426                          */
1427                         if (k  < li_task->b0 || k  >= li_task->b1 ||
1428                             kk < li_task->b0 || kk >= li_task->b1)
1429                         {
1430                             (*nCrossTaskTriangles)++;
1431                         }
1432
1433                         if (li_task->ntriangle == 0 ||
1434                             li_task->triangle[li_task->ntriangle - 1] < i)
1435                         {
1436                             /* Add this constraint to the triangle list */
1437                             li_task->triangle[li_task->ntriangle] = i;
1438                             li_task->tri_bits[li_task->ntriangle] = 0;
1439                             li_task->ntriangle++;
1440                             if (li->blnr[i+1] - li->blnr[i] > static_cast<int>(sizeof(li_task->tri_bits[0])*8 - 1))
1441                             {
1442                                 gmx_fatal(FARGS, "A constraint is connected to %d constraints, this is more than the %d allowed for constraints participating in triangles",
1443                                           li->blnr[i+1] - li->blnr[i],
1444                                           sizeof(li_task->tri_bits[0])*8-1);
1445                             }
1446                         }
1447                         li_task->tri_bits[li_task->ntriangle-1] |= (1 << (n - li->blnr[i]));
1448                         (*ncc_triangle)++;
1449                     }
1450                 }
1451             }
1452         }
1453     }
1454 }
1455
1456 /* Sets the elements in the LINCS matrix */
1457 void set_lincs_matrix(struct gmx_lincsdata *li, real *invmass, real lambda)
1458 {
1459     int        i;
1460     const real invsqrt2 = 0.7071067811865475244;
1461
1462     for (i = 0; (i < li->nc); i++)
1463     {
1464         int a1, a2;
1465
1466         a1          = li->bla[2*i];
1467         a2          = li->bla[2*i+1];
1468         li->blc[i]  = gmx_invsqrt(invmass[a1] + invmass[a2]);
1469         li->blc1[i] = invsqrt2;
1470     }
1471
1472     /* Construct the coupling coefficient matrix blmf */
1473     int th, ntriangle = 0, ncc_triangle = 0, nCrossTaskTriangles = 0;
1474 #pragma omp parallel for reduction(+: ntriangle, ncc_triangle, nCrossTaskTriangles) num_threads(li->ntask) schedule(static)
1475     for (th = 0; th < li->ntask; th++)
1476     {
1477         set_lincs_matrix_task(li, &li->task[th], invmass,
1478                               &ncc_triangle, &nCrossTaskTriangles);
1479         ntriangle = li->task[th].ntriangle;
1480     }
1481     li->ntriangle    = ntriangle;
1482     li->ncc_triangle = ncc_triangle;
1483     li->bTaskDepTri  = (nCrossTaskTriangles > 0);
1484
1485     if (debug)
1486     {
1487         fprintf(debug, "The %d constraints participate in %d triangles\n",
1488                 li->nc, li->ntriangle);
1489         fprintf(debug, "There are %d constraint couplings, of which %d in triangles\n",
1490                 li->ncc, li->ncc_triangle);
1491         if (li->ntriangle > 0 && li->ntask > 1)
1492         {
1493             fprintf(debug, "%d constraint triangles contain constraints assigned to different tasks\n",
1494                     nCrossTaskTriangles);
1495         }
1496     }
1497
1498     /* Set matlam,
1499      * so we know with which lambda value the masses have been set.
1500      */
1501     li->matlam = lambda;
1502 }
1503
1504 static int count_triangle_constraints(t_ilist *ilist, t_blocka *at2con)
1505 {
1506     int      ncon1, ncon_tot;
1507     int      c0, a00, a01, n1, c1, a10, a11, ac1, n2, c2, a20, a21;
1508     int      ncon_triangle;
1509     gmx_bool bTriangle;
1510     t_iatom *ia1, *ia2, *iap;
1511
1512     ncon1    = ilist[F_CONSTR].nr/3;
1513     ncon_tot = ncon1 + ilist[F_CONSTRNC].nr/3;
1514
1515     ia1 = ilist[F_CONSTR].iatoms;
1516     ia2 = ilist[F_CONSTRNC].iatoms;
1517
1518     ncon_triangle = 0;
1519     for (c0 = 0; c0 < ncon_tot; c0++)
1520     {
1521         bTriangle = FALSE;
1522         iap       = constr_iatomptr(ncon1, ia1, ia2, c0);
1523         a00       = iap[1];
1524         a01       = iap[2];
1525         for (n1 = at2con->index[a01]; n1 < at2con->index[a01+1]; n1++)
1526         {
1527             c1 = at2con->a[n1];
1528             if (c1 != c0)
1529             {
1530                 iap = constr_iatomptr(ncon1, ia1, ia2, c1);
1531                 a10 = iap[1];
1532                 a11 = iap[2];
1533                 if (a10 == a01)
1534                 {
1535                     ac1 = a11;
1536                 }
1537                 else
1538                 {
1539                     ac1 = a10;
1540                 }
1541                 for (n2 = at2con->index[ac1]; n2 < at2con->index[ac1+1]; n2++)
1542                 {
1543                     c2 = at2con->a[n2];
1544                     if (c2 != c0 && c2 != c1)
1545                     {
1546                         iap = constr_iatomptr(ncon1, ia1, ia2, c2);
1547                         a20 = iap[1];
1548                         a21 = iap[2];
1549                         if (a20 == a00 || a21 == a00)
1550                         {
1551                             bTriangle = TRUE;
1552                         }
1553                     }
1554                 }
1555             }
1556         }
1557         if (bTriangle)
1558         {
1559             ncon_triangle++;
1560         }
1561     }
1562
1563     return ncon_triangle;
1564 }
1565
1566 static gmx_bool more_than_two_sequential_constraints(const t_ilist  *ilist,
1567                                                      const t_blocka *at2con)
1568 {
1569     t_iatom  *ia1, *ia2, *iap;
1570     int       ncon1, ncon_tot, c;
1571     int       a1, a2;
1572     gmx_bool  bMoreThanTwoSequentialConstraints;
1573
1574     ncon1    = ilist[F_CONSTR].nr/3;
1575     ncon_tot = ncon1 + ilist[F_CONSTRNC].nr/3;
1576
1577     ia1 = ilist[F_CONSTR].iatoms;
1578     ia2 = ilist[F_CONSTRNC].iatoms;
1579
1580     bMoreThanTwoSequentialConstraints = FALSE;
1581     for (c = 0; c < ncon_tot && !bMoreThanTwoSequentialConstraints; c++)
1582     {
1583         iap = constr_iatomptr(ncon1, ia1, ia2, c);
1584         a1  = iap[1];
1585         a2  = iap[2];
1586         /* Check if this constraint has constraints connected at both atoms */
1587         if (at2con->index[a1+1] - at2con->index[a1] > 1 &&
1588             at2con->index[a2+1] - at2con->index[a2] > 1)
1589         {
1590             bMoreThanTwoSequentialConstraints = TRUE;
1591         }
1592     }
1593
1594     return bMoreThanTwoSequentialConstraints;
1595 }
1596
1597 static int int_comp(const void *a, const void *b)
1598 {
1599     return (*(int *)a) - (*(int *)b);
1600 }
1601
1602 gmx_lincsdata_t init_lincs(FILE *fplog, gmx_mtop_t *mtop,
1603                            int nflexcon_global, t_blocka *at2con,
1604                            gmx_bool bPLINCS, int nIter, int nProjOrder)
1605 {
1606     struct gmx_lincsdata *li;
1607     int                   mt, mb;
1608     gmx_moltype_t        *molt;
1609     gmx_bool              bMoreThanTwoSeq;
1610
1611     if (fplog)
1612     {
1613         fprintf(fplog, "\nInitializing%s LINear Constraint Solver\n",
1614                 bPLINCS ? " Parallel" : "");
1615     }
1616
1617     snew(li, 1);
1618
1619     li->ncg      =
1620         gmx_mtop_ftype_count(mtop, F_CONSTR) +
1621         gmx_mtop_ftype_count(mtop, F_CONSTRNC);
1622     li->ncg_flex = nflexcon_global;
1623
1624     li->nIter  = nIter;
1625     li->nOrder = nProjOrder;
1626
1627     li->max_connect = 0;
1628     for (mt = 0; mt < mtop->nmoltype; mt++)
1629     {
1630         int a;
1631
1632         molt = &mtop->moltype[mt];
1633         for (a = 0; a < molt->atoms.nr; a++)
1634         {
1635             li->max_connect = std::max(li->max_connect,
1636                                        at2con[mt].index[a + 1] - at2con[mt].index[a]);
1637         }
1638     }
1639
1640     li->ncg_triangle = 0;
1641     bMoreThanTwoSeq  = FALSE;
1642     for (mb = 0; mb < mtop->nmolblock; mb++)
1643     {
1644         molt              = &mtop->moltype[mtop->molblock[mb].type];
1645
1646         li->ncg_triangle +=
1647             mtop->molblock[mb].nmol*
1648             count_triangle_constraints(molt->ilist,
1649                                        &at2con[mtop->molblock[mb].type]);
1650
1651         if (!bMoreThanTwoSeq &&
1652             more_than_two_sequential_constraints(molt->ilist, &at2con[mtop->molblock[mb].type]))
1653         {
1654             bMoreThanTwoSeq = TRUE;
1655         }
1656     }
1657
1658     /* Check if we need to communicate not only before LINCS,
1659      * but also before each iteration.
1660      * The check for only two sequential constraints is only
1661      * useful for the common case of H-bond only constraints.
1662      * With more effort we could also make it useful for small
1663      * molecules with nr. sequential constraints <= nOrder-1.
1664      */
1665     li->bCommIter = (bPLINCS && (li->nOrder < 1 || bMoreThanTwoSeq));
1666
1667     if (debug && bPLINCS)
1668     {
1669         fprintf(debug, "PLINCS communication before each iteration: %d\n",
1670                 li->bCommIter);
1671     }
1672
1673     /* LINCS can run on any number of threads.
1674      * Currently the number is fixed for the whole simulation,
1675      * but it could be set in set_lincs().
1676      * The current constraint to task assignment code can create independent
1677      * tasks only when not more than two constraints are connected sequentially.
1678      */
1679     li->ntask    = gmx_omp_nthreads_get(emntLINCS);
1680     li->bTaskDep = (li->ntask > 1 && bMoreThanTwoSeq);
1681     if (debug)
1682     {
1683         fprintf(debug, "LINCS: using %d threads, tasks are %sdependent\n",
1684                 li->ntask, li->bTaskDep ? "" : "in");
1685     }
1686     if (li->ntask == 1)
1687     {
1688         snew(li->task, 1);
1689     }
1690     else
1691     {
1692         /* Allocate an extra elements for "task-overlap" constraints */
1693         snew(li->task, li->ntask + 1);
1694     }
1695     int th;
1696 #pragma omp parallel for num_threads(li->ntask)
1697     for (th = 0; th < li->ntask; th++)
1698     {
1699         /* Per thread SIMD load buffer for loading 2 simd_width rvecs */
1700         snew_aligned(li->task[th].simd_buf, 2*simd_width*DIM,
1701                      align_bytes);
1702     }
1703
1704     if (bPLINCS || li->ncg_triangle > 0)
1705     {
1706         please_cite(fplog, "Hess2008a");
1707     }
1708     else
1709     {
1710         please_cite(fplog, "Hess97a");
1711     }
1712
1713     if (fplog)
1714     {
1715         fprintf(fplog, "The number of constraints is %d\n", li->ncg);
1716         if (bPLINCS)
1717         {
1718             fprintf(fplog, "There are inter charge-group constraints,\n"
1719                     "will communicate selected coordinates each lincs iteration\n");
1720         }
1721         if (li->ncg_triangle > 0)
1722         {
1723             fprintf(fplog,
1724                     "%d constraints are involved in constraint triangles,\n"
1725                     "will apply an additional matrix expansion of order %d for couplings\n"
1726                     "between constraints inside triangles\n",
1727                     li->ncg_triangle, li->nOrder);
1728         }
1729     }
1730
1731     return li;
1732 }
1733
1734 /* Sets up the work division over the threads */
1735 static void lincs_thread_setup(struct gmx_lincsdata *li, int natoms)
1736 {
1737     lincs_task_t   *li_m;
1738     int             th;
1739     gmx_bitmask_t  *atf;
1740     int             a;
1741
1742     if (natoms > li->atf_nalloc)
1743     {
1744         li->atf_nalloc = over_alloc_large(natoms);
1745         srenew(li->atf, li->atf_nalloc);
1746     }
1747
1748     atf = li->atf;
1749     /* Clear the atom flags */
1750     for (a = 0; a < natoms; a++)
1751     {
1752         bitmask_clear(&atf[a]);
1753     }
1754
1755     if (li->ntask > BITMASK_SIZE)
1756     {
1757         gmx_fatal(FARGS, "More than %d threads is not supported for LINCS.", BITMASK_SIZE);
1758     }
1759
1760     for (th = 0; th < li->ntask; th++)
1761     {
1762         lincs_task_t *li_task;
1763         int           b;
1764
1765         li_task = &li->task[th];
1766
1767         /* For each atom set a flag for constraints from each */
1768         for (b = li_task->b0; b < li_task->b1; b++)
1769         {
1770             bitmask_set_bit(&atf[li->bla[b*2    ]], th);
1771             bitmask_set_bit(&atf[li->bla[b*2 + 1]], th);
1772         }
1773     }
1774
1775 #pragma omp parallel for num_threads(li->ntask) schedule(static)
1776     for (th = 0; th < li->ntask; th++)
1777     {
1778         lincs_task_t  *li_task;
1779         gmx_bitmask_t  mask;
1780         int            b;
1781
1782         li_task = &li->task[th];
1783
1784         if (li_task->b1 - li_task->b0 > li_task->ind_nalloc)
1785         {
1786             li_task->ind_nalloc = over_alloc_large(li_task->b1-li_task->b0);
1787             srenew(li_task->ind, li_task->ind_nalloc);
1788             srenew(li_task->ind_r, li_task->ind_nalloc);
1789         }
1790
1791         bitmask_init_low_bits(&mask, th);
1792
1793         li_task->nind   = 0;
1794         li_task->nind_r = 0;
1795         for (b = li_task->b0; b < li_task->b1; b++)
1796         {
1797             /* We let the constraint with the lowest thread index
1798              * operate on atoms with constraints from multiple threads.
1799              */
1800             if (bitmask_is_disjoint(atf[li->bla[b*2]], mask) &&
1801                 bitmask_is_disjoint(atf[li->bla[b*2+1]], mask))
1802             {
1803                 /* Add the constraint to the local atom update index */
1804                 li_task->ind[li_task->nind++] = b;
1805             }
1806             else
1807             {
1808                 /* Add the constraint to the rest block */
1809                 li_task->ind_r[li_task->nind_r++] = b;
1810             }
1811         }
1812     }
1813
1814     /* We need to copy all constraints which have not be assigned
1815      * to a thread to a separate list which will be handled by one thread.
1816      */
1817     li_m = &li->task[li->ntask];
1818
1819     li_m->nind = 0;
1820     for (th = 0; th < li->ntask; th++)
1821     {
1822         lincs_task_t *li_task;
1823         int           b;
1824
1825         li_task   = &li->task[th];
1826
1827         if (li_m->nind + li_task->nind_r > li_m->ind_nalloc)
1828         {
1829             li_m->ind_nalloc = over_alloc_large(li_m->nind+li_task->nind_r);
1830             srenew(li_m->ind, li_m->ind_nalloc);
1831         }
1832
1833         for (b = 0; b < li_task->nind_r; b++)
1834         {
1835             li_m->ind[li_m->nind++] = li_task->ind_r[b];
1836         }
1837
1838         if (debug)
1839         {
1840             fprintf(debug, "LINCS thread %d: %d constraints\n",
1841                     th, li_task->nind);
1842         }
1843     }
1844
1845     if (debug)
1846     {
1847         fprintf(debug, "LINCS thread r: %d constraints\n",
1848                 li_m->nind);
1849     }
1850 }
1851
1852 /* There is no realloc with alignment, so here we make one for reals.
1853  * Note that this function does not preserve the contents of the memory.
1854  */
1855 static void resize_real_aligned(real **ptr, int nelem)
1856 {
1857     sfree_aligned(*ptr);
1858     snew_aligned(*ptr, nelem, align_bytes);
1859 }
1860
1861 static void assign_constraint(struct gmx_lincsdata *li,
1862                               int constraint_index,
1863                               int a1, int a2,
1864                               real lenA, real lenB,
1865                               const t_blocka *at2con)
1866 {
1867     int con;
1868
1869     con = li->nc;
1870
1871     /* Make an mapping of local topology constraint index to LINCS index */
1872     li->con_index[constraint_index] = con;
1873
1874     li->bllen0[con]  = lenA;
1875     li->ddist[con]   = lenB - lenA;
1876     /* Set the length to the topology A length */
1877     li->bllen[con]   = lenA;
1878     li->bla[2*con]   = a1;
1879     li->bla[2*con+1] = a2;
1880
1881     /* Make space in the constraint connection matrix for constraints
1882      * connected to both end of the current constraint.
1883      */
1884     li->ncc +=
1885         at2con->index[a1 + 1] - at2con->index[a1] - 1 +
1886         at2con->index[a2 + 1] - at2con->index[a2] - 1;
1887
1888     li->blnr[con + 1] = li->ncc;
1889
1890     /* Increase the constraint count */
1891     li->nc++;
1892 }
1893
1894 /* Check if constraint with topology index constraint_index is connected
1895  * to other constraints, and if so add those connected constraints to our task.
1896  */
1897 static void check_assign_connected(struct gmx_lincsdata *li,
1898                                    const t_iatom *iatom,
1899                                    const t_idef *idef,
1900                                    int bDynamics,
1901                                    int a1, int a2,
1902                                    const t_blocka *at2con)
1903 {
1904     /* Currently this function only supports constraint groups
1905      * in which all constraints share at least one atom
1906      * (e.g. H-bond constraints).
1907      * Check both ends of the current constraint for
1908      * connected constraints. We need to assign those
1909      * to the same task.
1910      */
1911     int end;
1912
1913     for (end = 0; end < 2; end++)
1914     {
1915         int a, k;
1916
1917         a = (end == 0 ? a1 : a2);
1918
1919         for (k = at2con->index[a]; k < at2con->index[a + 1]; k++)
1920         {
1921             int cc;
1922
1923             cc = at2con->a[k];
1924             /* Check if constraint cc has not yet been assigned */
1925             if (li->con_index[cc] == -1)
1926             {
1927                 int  type;
1928                 real lenA, lenB;
1929
1930                 type = iatom[cc*3];
1931                 lenA = idef->iparams[type].constr.dA;
1932                 lenB = idef->iparams[type].constr.dB;
1933
1934                 if (bDynamics || lenA != 0 || lenB != 0)
1935                 {
1936                     assign_constraint(li, cc, iatom[3*cc + 1], iatom[3*cc + 2], lenA, lenB, at2con);
1937                 }
1938             }
1939         }
1940     }
1941 }
1942
1943 /* Check if constraint with topology index constraint_index is involved
1944  * in a constraint triangle, and if so add the other two constraints
1945  * in the triangle to our task.
1946  */
1947 static void check_assign_triangle(struct gmx_lincsdata *li,
1948                                   const t_iatom *iatom,
1949                                   const t_idef *idef,
1950                                   int bDynamics,
1951                                   int constraint_index,
1952                                   int a1, int a2,
1953                                   const t_blocka *at2con)
1954 {
1955     int nca, cc[32], ca[32], k;
1956     int c_triangle[2] = { -1, -1 };
1957
1958     nca = 0;
1959     for (k = at2con->index[a1]; k < at2con->index[a1 + 1]; k++)
1960     {
1961         int c;
1962
1963         c = at2con->a[k];
1964         if (c != constraint_index)
1965         {
1966             int aa1, aa2;
1967
1968             aa1 = iatom[c*3 + 1];
1969             aa2 = iatom[c*3 + 2];
1970             if (aa1 != a1)
1971             {
1972                 cc[nca] = c;
1973                 ca[nca] = aa1;
1974                 nca++;
1975             }
1976             if (aa2 != a1)
1977             {
1978                 cc[nca] = c;
1979                 ca[nca] = aa2;
1980                 nca++;
1981             }
1982         }
1983     }
1984
1985     for (k = at2con->index[a2]; k < at2con->index[a2 + 1]; k++)
1986     {
1987         int c;
1988
1989         c = at2con->a[k];
1990         if (c != constraint_index)
1991         {
1992             int aa1, aa2, i;
1993
1994             aa1 = iatom[c*3 + 1];
1995             aa2 = iatom[c*3 + 2];
1996             if (aa1 != a2)
1997             {
1998                 for (i = 0; i < nca; i++)
1999                 {
2000                     if (aa1 == ca[i])
2001                     {
2002                         c_triangle[0] = cc[i];
2003                         c_triangle[1] = c;
2004                     }
2005                 }
2006             }
2007             if (aa2 != a2)
2008             {
2009                 for (i = 0; i < nca; i++)
2010                 {
2011                     if (aa2 == ca[i])
2012                     {
2013                         c_triangle[0] = cc[i];
2014                         c_triangle[1] = c;
2015                     }
2016                 }
2017             }
2018         }
2019     }
2020
2021     if (c_triangle[0] >= 0)
2022     {
2023         int end;
2024
2025         for (end = 0; end < 2; end++)
2026         {
2027             /* Check if constraint c_triangle[end] has not yet been assigned */
2028             if (li->con_index[c_triangle[end]] == -1)
2029             {
2030                 int  i, type;
2031                 real lenA, lenB;
2032
2033                 i    = c_triangle[end]*3;
2034                 type = iatom[i];
2035                 lenA = idef->iparams[type].constr.dA;
2036                 lenB = idef->iparams[type].constr.dB;
2037
2038                 if (bDynamics || lenA != 0 || lenB != 0)
2039                 {
2040                     assign_constraint(li, c_triangle[end], iatom[i + 1], iatom[i + 2], lenA, lenB, at2con);
2041                 }
2042             }
2043         }
2044     }
2045 }
2046
2047 static void set_matrix_indices(struct gmx_lincsdata *li,
2048                                const lincs_task_t   *li_task,
2049                                const t_blocka       *at2con,
2050                                gmx_bool              bSortMatrix)
2051 {
2052     int b;
2053
2054     for (b = li_task->b0; b < li_task->b1; b++)
2055     {
2056         int a1, a2, i, k;
2057
2058         a1 = li->bla[b*2];
2059         a2 = li->bla[b*2 + 1];
2060
2061         i = li->blnr[b];
2062         for (k = at2con->index[a1]; k < at2con->index[a1 + 1]; k++)
2063         {
2064             int concon;
2065
2066             concon = li->con_index[at2con->a[k]];
2067             if (concon != b)
2068             {
2069                 li->blbnb[i++] = concon;
2070             }
2071         }
2072         for (k = at2con->index[a2]; k < at2con->index[a2 + 1]; k++)
2073         {
2074             int concon;
2075
2076             concon = li->con_index[at2con->a[k]];
2077             if (concon != b)
2078             {
2079                 li->blbnb[i++] = concon;
2080             }
2081         }
2082
2083         if (bSortMatrix)
2084         {
2085             /* Order the blbnb matrix to optimize memory access */
2086             qsort(&(li->blbnb[li->blnr[b]]), li->blnr[b + 1] - li->blnr[b],
2087                   sizeof(li->blbnb[0]), int_comp);
2088         }
2089     }
2090 }
2091
2092 void set_lincs(t_idef *idef, t_mdatoms *md,
2093                gmx_bool bDynamics, t_commrec *cr,
2094                struct gmx_lincsdata *li)
2095 {
2096     int          natoms, nflexcon;
2097     t_blocka     at2con;
2098     t_iatom     *iatom;
2099     int          i, ncc_alloc_old, ncon_tot;
2100
2101     li->nc_real = 0;
2102     li->nc      = 0;
2103     li->ncc     = 0;
2104     /* Zero the thread index ranges.
2105      * Otherwise without local constraints we could return with old ranges.
2106      */
2107     for (i = 0; i < li->ntask; i++)
2108     {
2109         li->task[i].b0   = 0;
2110         li->task[i].b1   = 0;
2111         li->task[i].nind = 0;
2112     }
2113     if (li->ntask > 1)
2114     {
2115         li->task[li->ntask].nind = 0;
2116     }
2117
2118     /* This is the local topology, so there are only F_CONSTR constraints */
2119     if (idef->il[F_CONSTR].nr == 0)
2120     {
2121         /* There are no constraints,
2122          * we do not need to fill any data structures.
2123          */
2124         return;
2125     }
2126
2127     if (debug)
2128     {
2129         fprintf(debug, "Building the LINCS connectivity\n");
2130     }
2131
2132     if (DOMAINDECOMP(cr))
2133     {
2134         if (cr->dd->constraints)
2135         {
2136             int start;
2137
2138             dd_get_constraint_range(cr->dd, &start, &natoms);
2139         }
2140         else
2141         {
2142             natoms = cr->dd->nat_home;
2143         }
2144     }
2145     else
2146     {
2147         natoms = md->homenr;
2148     }
2149     at2con = make_at2con(0, natoms, idef->il, idef->iparams, bDynamics,
2150                          &nflexcon);
2151
2152     ncon_tot = idef->il[F_CONSTR].nr/3;
2153
2154     /* Ensure we have enough padding for aligned loads for each thread */
2155     if (ncon_tot + li->ntask*simd_width > li->nc_alloc || li->nc_alloc == 0)
2156     {
2157         li->nc_alloc = over_alloc_dd(ncon_tot + li->ntask*simd_width);
2158         srenew(li->con_index, li->nc_alloc);
2159         resize_real_aligned(&li->bllen0, li->nc_alloc);
2160         resize_real_aligned(&li->ddist, li->nc_alloc);
2161         srenew(li->bla, 2*li->nc_alloc);
2162         resize_real_aligned(&li->blc, li->nc_alloc);
2163         resize_real_aligned(&li->blc1, li->nc_alloc);
2164         srenew(li->blnr, li->nc_alloc + 1);
2165         resize_real_aligned(&li->bllen, li->nc_alloc);
2166         srenew(li->tmpv, li->nc_alloc);
2167         if (DOMAINDECOMP(cr))
2168         {
2169             srenew(li->nlocat, li->nc_alloc);
2170         }
2171         resize_real_aligned(&li->tmp1, li->nc_alloc);
2172         resize_real_aligned(&li->tmp2, li->nc_alloc);
2173         resize_real_aligned(&li->tmp3, li->nc_alloc);
2174         resize_real_aligned(&li->tmp4, li->nc_alloc);
2175         resize_real_aligned(&li->mlambda, li->nc_alloc);
2176     }
2177
2178     iatom = idef->il[F_CONSTR].iatoms;
2179
2180     ncc_alloc_old = li->ncc_alloc;
2181     li->blnr[0]   = li->ncc;
2182
2183     /* Assign the constraints for li->ntask LINCS tasks.
2184      * We target a uniform distribution of constraints over the tasks.
2185      * Note that when flexible constraints are present, but are removed here
2186      * (e.g. because we are doing EM) we get imbalance, but since that doesn't
2187      * happen during normal MD, that's ok.
2188      */
2189     int ncon_assign, ncon_target, con, th;
2190
2191     /* Determine the number of constraints we need to assign here */
2192     ncon_assign      = ncon_tot;
2193     if (!bDynamics)
2194     {
2195         /* With energy minimization, flexible constraints are ignored
2196          * (and thus minimized, as they should be).
2197          */
2198         ncon_assign -= nflexcon;
2199     }
2200
2201     /* Set the target constraint count per task to exactly uniform,
2202      * this might be overridden below.
2203      */
2204     ncon_target = (ncon_assign + li->ntask - 1)/li->ntask;
2205
2206     /* Mark all constraints as unassigned by setting their index to -1 */
2207     for (con = 0; con < ncon_tot; con++)
2208     {
2209         li->con_index[con] = -1;
2210     }
2211
2212     con = 0;
2213     for (th = 0; th < li->ntask; th++)
2214     {
2215         lincs_task_t *li_task;
2216
2217         li_task = &li->task[th];
2218
2219 #ifdef LINCS_SIMD
2220         /* With indepedent tasks we likely have H-bond constraints or constraint
2221          * pairs. The connected constraints will be pulled into the task, so the
2222          * constraints per task will often exceed ncon_target.
2223          * Triangle constraints can also increase the count, but there are
2224          * relatively few of those, so we usually expect to get ncon_target.
2225          */
2226         if (li->bTaskDep)
2227         {
2228             /* We round ncon_target to a multiple of GMX_SIMD_WIDTH,
2229              * since otherwise a lot of operations can be wasted.
2230              * There are several ways to round here, we choose the one
2231              * that alternates block sizes, which helps with Intel HT.
2232              */
2233             ncon_target = ((ncon_assign*(th + 1))/li->ntask - li->nc_real + GMX_SIMD_REAL_WIDTH - 1) & ~(GMX_SIMD_REAL_WIDTH - 1);
2234         }
2235 #endif
2236
2237         /* Continue filling the arrays where we left off with the previous task,
2238          * including padding for SIMD.
2239          */
2240         li_task->b0 = li->nc;
2241
2242         while (con < ncon_tot && li->nc - li_task->b0 < ncon_target)
2243         {
2244             if (li->con_index[con] == -1)
2245             {
2246                 int  type, a1, a2;
2247                 real lenA, lenB;
2248
2249                 type   = iatom[3*con];
2250                 a1     = iatom[3*con + 1];
2251                 a2     = iatom[3*con + 2];
2252                 lenA   = idef->iparams[type].constr.dA;
2253                 lenB   = idef->iparams[type].constr.dB;
2254                 /* Skip the flexible constraints when not doing dynamics */
2255                 if (bDynamics || lenA != 0 || lenB != 0)
2256                 {
2257                     assign_constraint(li, con, a1, a2, lenA, lenB, &at2con);
2258
2259                     if (li->ntask > 1 && !li->bTaskDep)
2260                     {
2261                         /* We can generate independent tasks. Check if we
2262                          * need to assign connected constraints to our task.
2263                          */
2264                         check_assign_connected(li, iatom, idef, bDynamics,
2265                                                a1, a2, &at2con);
2266                     }
2267                     if (li->ntask > 1 && li->ncg_triangle > 0)
2268                     {
2269                         /* Ensure constraints in one triangle are assigned
2270                          * to the same task.
2271                          */
2272                         check_assign_triangle(li, iatom, idef, bDynamics,
2273                                               con, a1, a2, &at2con);
2274                     }
2275                 }
2276             }
2277
2278             con++;
2279         }
2280
2281         li_task->b1 = li->nc;
2282
2283         if (simd_width > 1)
2284         {
2285             /* Copy the last atom pair indices and lengths for constraints
2286              * up to a multiple of simd_width, such that we can do all
2287              * SIMD operations without having to worry about end effects.
2288              */
2289             int i, last;
2290
2291             li->nc = ((li_task->b1 + simd_width - 1)/simd_width)*simd_width;
2292             last   = li_task->b1 - 1;
2293             for (i = li_task->b1; i < li->nc; i++)
2294             {
2295                 li->bla[i*2    ] = li->bla[last*2    ];
2296                 li->bla[i*2 + 1] = li->bla[last*2 + 1];
2297                 li->bllen0[i]    = li->bllen0[last];
2298                 li->ddist[i]     = li->ddist[last];
2299                 li->bllen[i]     = li->bllen[last];
2300                 li->blnr[i + 1]  = li->blnr[last + 1];
2301             }
2302         }
2303
2304         /* Keep track of how many constraints we assigned */
2305         li->nc_real += li_task->b1 - li_task->b0;
2306
2307         if (debug)
2308         {
2309             fprintf(debug, "LINCS task %d constraints %d - %d\n",
2310                     th, li_task->b0, li_task->b1);
2311         }
2312     }
2313
2314     assert(li->nc_real == ncon_assign);
2315
2316     gmx_bool bSortMatrix;
2317
2318     /* Without DD we order the blbnb matrix to optimize memory access.
2319      * With DD the overhead of sorting is more than the gain during access.
2320      */
2321     bSortMatrix = !DOMAINDECOMP(cr);
2322
2323     if (li->ncc > li->ncc_alloc)
2324     {
2325         li->ncc_alloc = over_alloc_small(li->ncc);
2326         srenew(li->blbnb, li->ncc_alloc);
2327     }
2328
2329 #pragma omp parallel for num_threads(li->ntask) schedule(static)
2330     for (th = 0; th < li->ntask; th++)
2331     {
2332         lincs_task_t *li_task;
2333
2334         li_task = &li->task[th];
2335
2336         if (li->ncg_triangle > 0 &&
2337             li_task->b1 - li_task->b0 > li_task->tri_alloc)
2338         {
2339             /* This is allocating too much, but it is difficult to improve */
2340             li_task->tri_alloc = over_alloc_dd(li_task->b1 - li_task->b0);
2341             srenew(li_task->triangle, li_task->tri_alloc);
2342             srenew(li_task->tri_bits, li_task->tri_alloc);
2343         }
2344
2345         set_matrix_indices(li, li_task, &at2con, bSortMatrix);
2346     }
2347
2348     done_blocka(&at2con);
2349
2350     if (cr->dd == NULL)
2351     {
2352         /* Since the matrix is static, we should free some memory */
2353         li->ncc_alloc = li->ncc;
2354         srenew(li->blbnb, li->ncc_alloc);
2355     }
2356
2357     if (li->ncc_alloc > ncc_alloc_old)
2358     {
2359         srenew(li->blmf, li->ncc_alloc);
2360         srenew(li->blmf1, li->ncc_alloc);
2361         srenew(li->tmpncc, li->ncc_alloc);
2362     }
2363
2364     if (DOMAINDECOMP(cr) && dd_constraints_nlocalatoms(cr->dd) != NULL)
2365     {
2366         int *nlocat_dd;
2367
2368         nlocat_dd = dd_constraints_nlocalatoms(cr->dd);
2369
2370         /* Convert nlocat from local topology to LINCS constraint indexing */
2371         for (con = 0; con < ncon_tot; con++)
2372         {
2373             li->nlocat[li->con_index[con]] = nlocat_dd[con];
2374         }
2375     }
2376     else
2377     {
2378         li->nlocat = NULL;
2379     }
2380
2381     if (debug)
2382     {
2383         fprintf(debug, "Number of constraints is %d, padded %d, couplings %d\n",
2384                 li->nc_real, li->nc, li->ncc);
2385     }
2386
2387     if (li->ntask > 1)
2388     {
2389         lincs_thread_setup(li, md->nr);
2390     }
2391
2392     set_lincs_matrix(li, md->invmass, md->lambda);
2393 }
2394
2395 static void lincs_warning(FILE *fplog,
2396                           gmx_domdec_t *dd, rvec *x, rvec *xprime, t_pbc *pbc,
2397                           int ncons, int *bla, real *bllen, real wangle,
2398                           int maxwarn, int *warncount)
2399 {
2400     int  b, i, j;
2401     rvec v0, v1;
2402     real wfac, d0, d1, cosine;
2403     char buf[STRLEN];
2404
2405     wfac = cos(DEG2RAD*wangle);
2406
2407     sprintf(buf, "bonds that rotated more than %g degrees:\n"
2408             " atom 1 atom 2  angle  previous, current, constraint length\n",
2409             wangle);
2410     fprintf(stderr, "%s", buf);
2411     if (fplog)
2412     {
2413         fprintf(fplog, "%s", buf);
2414     }
2415
2416     for (b = 0; b < ncons; b++)
2417     {
2418         i = bla[2*b];
2419         j = bla[2*b+1];
2420         if (pbc)
2421         {
2422             pbc_dx_aiuc(pbc, x[i], x[j], v0);
2423             pbc_dx_aiuc(pbc, xprime[i], xprime[j], v1);
2424         }
2425         else
2426         {
2427             rvec_sub(x[i], x[j], v0);
2428             rvec_sub(xprime[i], xprime[j], v1);
2429         }
2430         d0     = norm(v0);
2431         d1     = norm(v1);
2432         cosine = iprod(v0, v1)/(d0*d1);
2433         if (cosine < wfac)
2434         {
2435             sprintf(buf, " %6d %6d  %5.1f  %8.4f %8.4f    %8.4f\n",
2436                     ddglatnr(dd, i), ddglatnr(dd, j),
2437                     RAD2DEG*acos(cosine), d0, d1, bllen[b]);
2438             fprintf(stderr, "%s", buf);
2439             if (fplog)
2440             {
2441                 fprintf(fplog, "%s", buf);
2442             }
2443             if (!gmx_isfinite(d1))
2444             {
2445                 gmx_fatal(FARGS, "Bond length not finite.");
2446             }
2447
2448             (*warncount)++;
2449         }
2450     }
2451     if (*warncount > maxwarn)
2452     {
2453         too_many_constraint_warnings(econtLINCS, *warncount);
2454     }
2455 }
2456
2457 static void cconerr(const struct gmx_lincsdata *lincsd,
2458                     rvec *x, t_pbc *pbc,
2459                     real *ncons_loc, real *ssd, real *max, int *imax)
2460 {
2461     const int  *bla, *nlocat;
2462     const real *bllen;
2463     real        ma, ssd2;
2464     int         count, im, task;
2465
2466     bla    = lincsd->bla;
2467     bllen  = lincsd->bllen;
2468     nlocat = lincsd->nlocat;
2469
2470     ma    = 0;
2471     ssd2  = 0;
2472     im    = 0;
2473     count = 0;
2474     for (task = 0; task < lincsd->ntask; task++)
2475     {
2476         int b;
2477
2478         for (b = lincsd->task[task].b0; b < lincsd->task[task].b1; b++)
2479         {
2480             real len, d, r2;
2481             rvec dx;
2482
2483             if (pbc)
2484             {
2485                 pbc_dx_aiuc(pbc, x[bla[2*b]], x[bla[2*b+1]], dx);
2486             }
2487             else
2488             {
2489                 rvec_sub(x[bla[2*b]], x[bla[2*b+1]], dx);
2490             }
2491             r2  = norm2(dx);
2492             len = r2*gmx_invsqrt(r2);
2493             d   = fabs(len/bllen[b]-1);
2494             if (d > ma && (nlocat == NULL || nlocat[b]))
2495             {
2496                 ma = d;
2497                 im = b;
2498             }
2499             if (nlocat == NULL)
2500             {
2501                 ssd2 += d*d;
2502                 count++;
2503             }
2504             else
2505             {
2506                 ssd2  += nlocat[b]*d*d;
2507                 count += nlocat[b];
2508             }
2509         }
2510     }
2511
2512     *ncons_loc = (nlocat ? 0.5 : 1)*count;
2513     *ssd       = (nlocat ? 0.5 : 1)*ssd2;
2514     *max       = ma;
2515     *imax      = im;
2516 }
2517
2518 gmx_bool constrain_lincs(FILE *fplog, gmx_bool bLog, gmx_bool bEner,
2519                          t_inputrec *ir,
2520                          gmx_int64_t step,
2521                          struct gmx_lincsdata *lincsd, t_mdatoms *md,
2522                          t_commrec *cr,
2523                          rvec *x, rvec *xprime, rvec *min_proj,
2524                          matrix box, t_pbc *pbc,
2525                          real lambda, real *dvdlambda,
2526                          real invdt, rvec *v,
2527                          gmx_bool bCalcVir, tensor vir_r_m_dr,
2528                          int econq,
2529                          t_nrnb *nrnb,
2530                          int maxwarn, int *warncount)
2531 {
2532     gmx_bool  bCalcDHDL;
2533     char      buf[STRLEN], buf2[22], buf3[STRLEN];
2534     int       i, p_imax;
2535     real      ncons_loc, p_ssd, p_max = 0;
2536     rvec      dx;
2537     gmx_bool  bOK, bWarn;
2538
2539     bOK = TRUE;
2540
2541     /* This boolean should be set by a flag passed to this routine.
2542      * We can also easily check if any constraint length is changed,
2543      * if not dH/dlambda=0 and we can also set the boolean to FALSE.
2544      */
2545     bCalcDHDL = (ir->efep != efepNO && dvdlambda != NULL);
2546
2547     if (lincsd->nc == 0 && cr->dd == NULL)
2548     {
2549         if (bLog || bEner)
2550         {
2551             lincsd->rmsd_data[0] = 0;
2552             if (ir->eI == eiSD2 && v == NULL)
2553             {
2554                 i = 2;
2555             }
2556             else
2557             {
2558                 i = 1;
2559             }
2560             lincsd->rmsd_data[i] = 0;
2561         }
2562
2563         return bOK;
2564     }
2565
2566     if (econq == econqCoord)
2567     {
2568         /* We can't use bCalcDHDL here, since NULL can be passed for dvdlambda
2569          * also with efep!=fepNO.
2570          */
2571         if (ir->efep != efepNO)
2572         {
2573             if (md->nMassPerturbed && lincsd->matlam != md->lambda)
2574             {
2575                 set_lincs_matrix(lincsd, md->invmass, md->lambda);
2576             }
2577
2578             for (i = 0; i < lincsd->nc; i++)
2579             {
2580                 lincsd->bllen[i] = lincsd->bllen0[i] + lambda*lincsd->ddist[i];
2581             }
2582         }
2583
2584         if (lincsd->ncg_flex)
2585         {
2586             /* Set the flexible constraint lengths to the old lengths */
2587             if (pbc != NULL)
2588             {
2589                 for (i = 0; i < lincsd->nc; i++)
2590                 {
2591                     if (lincsd->bllen[i] == 0)
2592                     {
2593                         pbc_dx_aiuc(pbc, x[lincsd->bla[2*i]], x[lincsd->bla[2*i+1]], dx);
2594                         lincsd->bllen[i] = norm(dx);
2595                     }
2596                 }
2597             }
2598             else
2599             {
2600                 for (i = 0; i < lincsd->nc; i++)
2601                 {
2602                     if (lincsd->bllen[i] == 0)
2603                     {
2604                         lincsd->bllen[i] =
2605                             sqrt(distance2(x[lincsd->bla[2*i]],
2606                                            x[lincsd->bla[2*i+1]]));
2607                     }
2608                 }
2609             }
2610         }
2611
2612         if (bLog && fplog)
2613         {
2614             cconerr(lincsd, xprime, pbc,
2615                     &ncons_loc, &p_ssd, &p_max, &p_imax);
2616         }
2617
2618         /* This bWarn var can be updated by multiple threads
2619          * at the same time. But as we only need to detect
2620          * if a warning occured or not, this is not an issue.
2621          */
2622         bWarn = FALSE;
2623
2624         /* The OpenMP parallel region of constrain_lincs for coords */
2625 #pragma omp parallel num_threads(lincsd->ntask)
2626         {
2627             int th = gmx_omp_get_thread_num();
2628
2629             clear_mat(lincsd->task[th].vir_r_m_dr);
2630
2631             do_lincs(x, xprime, box, pbc, lincsd, th,
2632                      md->invmass, cr,
2633                      bCalcDHDL,
2634                      ir->LincsWarnAngle, &bWarn,
2635                      invdt, v, bCalcVir,
2636                      th == 0 ? vir_r_m_dr : lincsd->task[th].vir_r_m_dr);
2637         }
2638
2639         if (bLog && fplog && lincsd->nc > 0)
2640         {
2641             fprintf(fplog, "   Rel. Constraint Deviation:  RMS         MAX     between atoms\n");
2642             fprintf(fplog, "       Before LINCS          %.6f    %.6f %6d %6d\n",
2643                     sqrt(p_ssd/ncons_loc), p_max,
2644                     ddglatnr(cr->dd, lincsd->bla[2*p_imax]),
2645                     ddglatnr(cr->dd, lincsd->bla[2*p_imax+1]));
2646         }
2647         if (bLog || bEner)
2648         {
2649             cconerr(lincsd, xprime, pbc,
2650                     &ncons_loc, &p_ssd, &p_max, &p_imax);
2651             /* Check if we are doing the second part of SD */
2652             if (ir->eI == eiSD2 && v == NULL)
2653             {
2654                 i = 2;
2655             }
2656             else
2657             {
2658                 i = 1;
2659             }
2660             lincsd->rmsd_data[0] = ncons_loc;
2661             lincsd->rmsd_data[i] = p_ssd;
2662         }
2663         else
2664         {
2665             lincsd->rmsd_data[0] = 0;
2666             lincsd->rmsd_data[1] = 0;
2667             lincsd->rmsd_data[2] = 0;
2668         }
2669         if (bLog && fplog && lincsd->nc > 0)
2670         {
2671             fprintf(fplog,
2672                     "        After LINCS          %.6f    %.6f %6d %6d\n\n",
2673                     sqrt(p_ssd/ncons_loc), p_max,
2674                     ddglatnr(cr->dd, lincsd->bla[2*p_imax]),
2675                     ddglatnr(cr->dd, lincsd->bla[2*p_imax+1]));
2676         }
2677
2678         if (bWarn)
2679         {
2680             if (maxwarn >= 0)
2681             {
2682                 cconerr(lincsd, xprime, pbc,
2683                         &ncons_loc, &p_ssd, &p_max, &p_imax);
2684                 if (MULTISIM(cr))
2685                 {
2686                     sprintf(buf3, " in simulation %d", cr->ms->sim);
2687                 }
2688                 else
2689                 {
2690                     buf3[0] = 0;
2691                 }
2692                 sprintf(buf, "\nStep %s, time %g (ps)  LINCS WARNING%s\n"
2693                         "relative constraint deviation after LINCS:\n"
2694                         "rms %.6f, max %.6f (between atoms %d and %d)\n",
2695                         gmx_step_str(step, buf2), ir->init_t+step*ir->delta_t,
2696                         buf3,
2697                         sqrt(p_ssd/ncons_loc), p_max,
2698                         ddglatnr(cr->dd, lincsd->bla[2*p_imax]),
2699                         ddglatnr(cr->dd, lincsd->bla[2*p_imax+1]));
2700                 if (fplog)
2701                 {
2702                     fprintf(fplog, "%s", buf);
2703                 }
2704                 fprintf(stderr, "%s", buf);
2705                 lincs_warning(fplog, cr->dd, x, xprime, pbc,
2706                               lincsd->nc, lincsd->bla, lincsd->bllen,
2707                               ir->LincsWarnAngle, maxwarn, warncount);
2708             }
2709             bOK = (p_max < 0.5);
2710         }
2711
2712         if (lincsd->ncg_flex)
2713         {
2714             for (i = 0; (i < lincsd->nc); i++)
2715             {
2716                 if (lincsd->bllen0[i] == 0 && lincsd->ddist[i] == 0)
2717                 {
2718                     lincsd->bllen[i] = 0;
2719                 }
2720             }
2721         }
2722     }
2723     else
2724     {
2725         /* The OpenMP parallel region of constrain_lincs for derivatives */
2726 #pragma omp parallel num_threads(lincsd->ntask)
2727         {
2728             int th = gmx_omp_get_thread_num();
2729
2730             do_lincsp(x, xprime, min_proj, pbc, lincsd, th,
2731                       md->invmass, econq, bCalcDHDL,
2732                       bCalcVir, th == 0 ? vir_r_m_dr : lincsd->task[th].vir_r_m_dr);
2733         }
2734     }
2735
2736     if (bCalcDHDL)
2737     {
2738         /* Reduce the dH/dlambda contributions over the threads */
2739         real dhdlambda;
2740         int  th;
2741
2742         dhdlambda = 0;
2743         for (th = 0; th < lincsd->ntask; th++)
2744         {
2745             dhdlambda += lincsd->task[th].dhdlambda;
2746         }
2747         if (econq == econqCoord)
2748         {
2749             /* dhdlambda contains dH/dlambda*dt^2, correct for this */
2750             /* TODO This should probably use invdt, so that sd integrator scaling works properly */
2751             dhdlambda /= ir->delta_t*ir->delta_t;
2752         }
2753         *dvdlambda += dhdlambda;
2754     }
2755
2756     if (bCalcVir && lincsd->ntask > 1)
2757     {
2758         for (i = 1; i < lincsd->ntask; i++)
2759         {
2760             m_add(vir_r_m_dr, lincsd->task[i].vir_r_m_dr, vir_r_m_dr);
2761         }
2762     }
2763
2764     /* count assuming nit=1 */
2765     inc_nrnb(nrnb, eNR_LINCS, lincsd->nc_real);
2766     inc_nrnb(nrnb, eNR_LINCSMAT, (2+lincsd->nOrder)*lincsd->ncc);
2767     if (lincsd->ntriangle > 0)
2768     {
2769         inc_nrnb(nrnb, eNR_LINCSMAT, lincsd->nOrder*lincsd->ncc_triangle);
2770     }
2771     if (v)
2772     {
2773         inc_nrnb(nrnb, eNR_CONSTR_V, lincsd->nc_real*2);
2774     }
2775     if (bCalcVir)
2776     {
2777         inc_nrnb(nrnb, eNR_CONSTR_VIR, lincsd->nc_real);
2778     }
2779
2780     return bOK;
2781 }