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