src/gromacs/mdlib/clincs.c

   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, by the GROMACS development team, led by
   7  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   8  * and including many others, as listed in the AUTHORS file in the
   9  * top-level source directory and at http://www.gromacs.org.
  10  *
  11  * GROMACS is free software; you can redistribute it and/or
  12  * modify it under the terms of the GNU Lesser General Public License
  13  * as published by the Free Software Foundation; either version 2.1
  14  * of the License, or (at your option) any later version.
  15  *
  16  * GROMACS is distributed in the hope that it will be useful,
  17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19  * Lesser General Public License for more details.
  20  *
  21  * You should have received a copy of the GNU Lesser General Public
  22  * License along with GROMACS; if not, see
  23  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  24  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  25  *
  26  * If you want to redistribute modifications to GROMACS, please
  27  * consider that scientific software is very special. Version
  28  * control is crucial - bugs must be traceable. We will be happy to
  29  * consider code for inclusion in the official distribution, but
  30  * derived work must not be called official GROMACS. Details are found
  31  * in the README & COPYING files - if they are missing, get the
  32  * official version at http://www.gromacs.org.
  33  *
  34  * To help us fund GROMACS development, we humbly ask that you cite
  35  * the research papers on the package. Check out http://www.gromacs.org.
  36  */
  37 /* This file is completely threadsafe - keep it that way! */
  38 #ifdef HAVE_CONFIG_H
  39 #include <config.h>
  40 #endif
  41
  42 #include <math.h>
  43 #include "main.h"
  44 #include "constr.h"
  45 #include "copyrite.h"
  46 #include "physics.h"
  47 #include "vec.h"
  48 #include "pbc.h"
  49 #include "smalloc.h"
  50 #include "mdrun.h"
  51 #include "nrnb.h"
  52 #include "domdec.h"
  53 #include "mtop_util.h"
  54 #include "gmx_omp_nthreads.h"
  55
  56 #include "gromacs/fileio/gmxfio.h"
  57 #include "gromacs/utility/gmxomp.h"
  58
  59 typedef struct {
  60     int    b0;         /* first constraint for this thread */
  61     int    b1;         /* b1-1 is the last constraint for this thread */
  62     int    nind;       /* number of indices */
  63     int   *ind;        /* constraint index for updating atom data */
  64     int    nind_r;     /* number of indices */
  65     int   *ind_r;      /* constraint index for updating atom data */
  66     int    ind_nalloc; /* allocation size of ind and ind_r */
  67     tensor vir_r_m_dr; /* temporary variable for virial calculation */
  68 } lincs_thread_t;
  69
  70 typedef struct gmx_lincsdata {
  71     int             ncg;          /* the global number of constraints */
  72     int             ncg_flex;     /* the global number of flexible constraints */
  73     int             ncg_triangle; /* the global number of constraints in triangles */
  74     int             nIter;        /* the number of iterations */
  75     int             nOrder;       /* the order of the matrix expansion */
  76     int             nc;           /* the number of constraints */
  77     int             nc_alloc;     /* the number we allocated memory for */
  78     int             ncc;          /* the number of constraint connections */
  79     int             ncc_alloc;    /* the number we allocated memory for */
  80     real            matlam;       /* the FE lambda value used for filling blc and blmf */
  81     real           *bllen0;       /* the reference distance in topology A */
  82     real           *ddist;        /* the reference distance in top B - the r.d. in top A */
  83     int            *bla;          /* the atom pairs involved in the constraints */
  84     real           *blc;          /* 1/sqrt(invmass1 + invmass2) */
  85     real           *blc1;         /* as blc, but with all masses 1 */
  86     int            *blnr;         /* index into blbnb and blmf */
  87     int            *blbnb;        /* list of constraint connections */
  88     int             ntriangle;    /* the local number of constraints in triangles */
  89     int            *triangle;     /* the list of triangle constraints */
  90     int            *tri_bits;     /* the bits tell if the matrix element should be used */
  91     int             ncc_triangle; /* the number of constraint connections in triangles */
  92     gmx_bool        bCommIter;    /* communicate before each LINCS interation */
  93     real           *blmf;         /* matrix of mass factors for constraint connections */
  94     real           *blmf1;        /* as blmf, but with all masses 1 */
  95     real           *bllen;        /* the reference bond length */
  96     int             nth;          /* The number of threads doing LINCS */
  97     lincs_thread_t *th;           /* LINCS thread division */
  98     unsigned       *atf;          /* atom flags for thread parallelization */
  99     int             atf_nalloc;   /* allocation size of atf */
 100     /* arrays for temporary storage in the LINCS algorithm */
 101     rvec           *tmpv;
 102     real           *tmpncc;
 103     real           *tmp1;
 104     real           *tmp2;
 105     real           *tmp3;
 106     real           *tmp4;
 107     real           *mlambda; /* the Lagrange multipliers * -1 */
 108     /* storage for the constraint RMS relative deviation output */
 109     real            rmsd_data[3];
 110 } t_gmx_lincsdata;
 111
 112 real *lincs_rmsd_data(struct gmx_lincsdata *lincsd)
 113 {
 114     return lincsd->rmsd_data;
 115 }
 116
 117 real lincs_rmsd(struct gmx_lincsdata *lincsd, gmx_bool bSD2)
 118 {
 119     if (lincsd->rmsd_data[0] > 0)
 120     {
 121         return sqrt(lincsd->rmsd_data[bSD2 ? 2 : 1]/lincsd->rmsd_data[0]);
 122     }
 123     else
 124     {
 125         return 0;
 126     }
 127 }
 128
 129 /* Do a set of nrec LINCS matrix multiplications.
 130  * This function will return with up to date thread-local
 131  * constraint data, without an OpenMP barrier.
 132  */
 133 static void lincs_matrix_expand(const struct gmx_lincsdata *lincsd,
 134                                 int b0, int b1,
 135                                 const real *blcc,
 136                                 real *rhs1, real *rhs2, real *sol)
 137 {
 138     int        nrec, rec, b, j, n, nr0, nr1;
 139     real       mvb, *swap;
 140     int        ntriangle, tb, bits;
 141     const int *blnr     = lincsd->blnr, *blbnb = lincsd->blbnb;
 142     const int *triangle = lincsd->triangle, *tri_bits = lincsd->tri_bits;
 143
 144     ntriangle = lincsd->ntriangle;
 145     nrec      = lincsd->nOrder;
 146
 147     for (rec = 0; rec < nrec; rec++)
 148     {
 149 #pragma omp barrier
 150         for (b = b0; b < b1; b++)
 151         {
 152             mvb = 0;
 153             for (n = blnr[b]; n < blnr[b+1]; n++)
 154             {
 155                 j   = blbnb[n];
 156                 mvb = mvb + blcc[n]*rhs1[j];
 157             }
 158             rhs2[b] = mvb;
 159             sol[b]  = sol[b] + mvb;
 160         }
 161         swap = rhs1;
 162         rhs1 = rhs2;
 163         rhs2 = swap;
 164     } /* nrec*(ncons+2*nrtot) flops */
 165
 166     if (ntriangle > 0)
 167     {
 168         /* Perform an extra nrec recursions for only the constraints
 169          * involved in rigid triangles.
 170          * In this way their accuracy should come close to those of the other
 171          * constraints, since traingles of constraints can produce eigenvalues
 172          * around 0.7, while the effective eigenvalue for bond constraints
 173          * is around 0.4 (and 0.7*0.7=0.5).
 174          */
 175         /* We need to copy the temporary array, since only the elements
 176          * for constraints involved in triangles are updated and then
 177          * the pointers are swapped. This saving copying the whole arrary.
 178          * We need barrier as other threads might still be reading from rhs2.
 179          */
 180 #pragma omp barrier
 181         for (b = b0; b < b1; b++)
 182         {
 183             rhs2[b] = rhs1[b];
 184         }
 185 #pragma omp barrier
 186 #pragma omp master
 187         {
 188             for (rec = 0; rec < nrec; rec++)
 189             {
 190                 for (tb = 0; tb < ntriangle; tb++)
 191                 {
 192                     b    = triangle[tb];
 193                     bits = tri_bits[tb];
 194                     mvb  = 0;
 195                     nr0  = blnr[b];
 196                     nr1  = blnr[b+1];
 197                     for (n = nr0; n < nr1; n++)
 198                     {
 199                         if (bits & (1<<(n-nr0)))
 200                         {
 201                             j   = blbnb[n];
 202                             mvb = mvb + blcc[n]*rhs1[j];
 203                         }
 204                     }
 205                     rhs2[b] = mvb;
 206                     sol[b]  = sol[b] + mvb;
 207                 }
 208                 swap = rhs1;
 209                 rhs1 = rhs2;
 210                 rhs2 = swap;
 211             }
 212         } /* flops count is missing here */
 213
 214         /* We need a barrier here as the calling routine will continue
 215          * to operate on the thread-local constraints without barrier.
 216          */
 217 #pragma omp barrier
 218     }
 219 }
 220
 221 static void lincs_update_atoms_noind(int ncons, const int *bla,
 222                                      real prefac,
 223                                      const real *fac, rvec *r,
 224                                      const real *invmass,
 225                                      rvec *x)
 226 {
 227     int  b, i, j;
 228     real mvb, im1, im2, tmp0, tmp1, tmp2;
 229
 230     if (invmass != NULL)
 231     {
 232         for (b = 0; b < ncons; b++)
 233         {
 234             i        = bla[2*b];
 235             j        = bla[2*b+1];
 236             mvb      = prefac*fac[b];
 237             im1      = invmass[i];
 238             im2      = invmass[j];
 239             tmp0     = r[b][0]*mvb;
 240             tmp1     = r[b][1]*mvb;
 241             tmp2     = r[b][2]*mvb;
 242             x[i][0] -= tmp0*im1;
 243             x[i][1] -= tmp1*im1;
 244             x[i][2] -= tmp2*im1;
 245             x[j][0] += tmp0*im2;
 246             x[j][1] += tmp1*im2;
 247             x[j][2] += tmp2*im2;
 248         } /* 16 ncons flops */
 249     }
 250     else
 251     {
 252         for (b = 0; b < ncons; b++)
 253         {
 254             i        = bla[2*b];
 255             j        = bla[2*b+1];
 256             mvb      = prefac*fac[b];
 257             tmp0     = r[b][0]*mvb;
 258             tmp1     = r[b][1]*mvb;
 259             tmp2     = r[b][2]*mvb;
 260             x[i][0] -= tmp0;
 261             x[i][1] -= tmp1;
 262             x[i][2] -= tmp2;
 263             x[j][0] += tmp0;
 264             x[j][1] += tmp1;
 265             x[j][2] += tmp2;
 266         }
 267     }
 268 }
 269
 270 static void lincs_update_atoms_ind(int ncons, const int *ind, const int *bla,
 271                                    real prefac,
 272                                    const real *fac, rvec *r,
 273                                    const real *invmass,
 274                                    rvec *x)
 275 {
 276     int  bi, b, i, j;
 277     real mvb, im1, im2, tmp0, tmp1, tmp2;
 278
 279     if (invmass != NULL)
 280     {
 281         for (bi = 0; bi < ncons; bi++)
 282         {
 283             b        = ind[bi];
 284             i        = bla[2*b];
 285             j        = bla[2*b+1];
 286             mvb      = prefac*fac[b];
 287             im1      = invmass[i];
 288             im2      = invmass[j];
 289             tmp0     = r[b][0]*mvb;
 290             tmp1     = r[b][1]*mvb;
 291             tmp2     = r[b][2]*mvb;
 292             x[i][0] -= tmp0*im1;
 293             x[i][1] -= tmp1*im1;
 294             x[i][2] -= tmp2*im1;
 295             x[j][0] += tmp0*im2;
 296             x[j][1] += tmp1*im2;
 297             x[j][2] += tmp2*im2;
 298         } /* 16 ncons flops */
 299     }
 300     else
 301     {
 302         for (bi = 0; bi < ncons; bi++)
 303         {
 304             b        = ind[bi];
 305             i        = bla[2*b];
 306             j        = bla[2*b+1];
 307             mvb      = prefac*fac[b];
 308             tmp0     = r[b][0]*mvb;
 309             tmp1     = r[b][1]*mvb;
 310             tmp2     = r[b][2]*mvb;
 311             x[i][0] -= tmp0;
 312             x[i][1] -= tmp1;
 313             x[i][2] -= tmp2;
 314             x[j][0] += tmp0;
 315             x[j][1] += tmp1;
 316             x[j][2] += tmp2;
 317         } /* 16 ncons flops */
 318     }
 319 }
 320
 321 static void lincs_update_atoms(struct gmx_lincsdata *li, int th,
 322                                real prefac,
 323                                const real *fac, rvec *r,
 324                                const real *invmass,
 325                                rvec *x)
 326 {
 327     if (li->nth == 1)
 328     {
 329         /* Single thread, we simply update for all constraints */
 330         lincs_update_atoms_noind(li->nc, li->bla, prefac, fac, r, invmass, x);
 331     }
 332     else
 333     {
 334         /* Update the atom vector components for our thread local
 335          * constraints that only access our local atom range.
 336          * This can be done without a barrier.
 337          */
 338         lincs_update_atoms_ind(li->th[th].nind, li->th[th].ind,
 339                                li->bla, prefac, fac, r, invmass, x);
 340
 341         if (li->th[li->nth].nind > 0)
 342         {
 343             /* Update the constraints that operate on atoms
 344              * in multiple thread atom blocks on the master thread.
 345              */
 346 #pragma omp barrier
 347 #pragma omp master
 348             {
 349                 lincs_update_atoms_ind(li->th[li->nth].nind,
 350                                        li->th[li->nth].ind,
 351                                        li->bla, prefac, fac, r, invmass, x);
 352             }
 353         }
 354     }
 355 }
 356
 357 /* LINCS projection, works on derivatives of the coordinates */
 358 static void do_lincsp(rvec *x, rvec *f, rvec *fp, t_pbc *pbc,
 359                       struct gmx_lincsdata *lincsd, int th,
 360                       real *invmass,
 361                       int econq, real *dvdlambda,
 362                       gmx_bool bCalcVir, tensor rmdf)
 363 {
 364     int      b0, b1, b, i, j, k, n;
 365     real     tmp0, tmp1, tmp2, im1, im2, mvb, rlen, len, wfac, lam;
 366     rvec     dx;
 367     int     *bla, *blnr, *blbnb;
 368     rvec    *r;
 369     real    *blc, *blmf, *blcc, *rhs1, *rhs2, *sol;
 370
 371     b0 = lincsd->th[th].b0;
 372     b1 = lincsd->th[th].b1;
 373
 374     bla    = lincsd->bla;
 375     r      = lincsd->tmpv;
 376     blnr   = lincsd->blnr;
 377     blbnb  = lincsd->blbnb;
 378     if (econq != econqForce)
 379     {
 380         /* Use mass-weighted parameters */
 381         blc  = lincsd->blc;
 382         blmf = lincsd->blmf;
 383     }
 384     else
 385     {
 386         /* Use non mass-weighted parameters */
 387         blc  = lincsd->blc1;
 388         blmf = lincsd->blmf1;
 389     }
 390     blcc   = lincsd->tmpncc;
 391     rhs1   = lincsd->tmp1;
 392     rhs2   = lincsd->tmp2;
 393     sol    = lincsd->tmp3;
 394
 395     /* Compute normalized i-j vectors */
 396     if (pbc)
 397     {
 398         for (b = b0; b < b1; b++)
 399         {
 400             pbc_dx_aiuc(pbc, x[bla[2*b]], x[bla[2*b+1]], dx);
 401             unitv(dx, r[b]);
 402         }
 403     }
 404     else
 405     {
 406         for (b = b0; b < b1; b++)
 407         {
 408             rvec_sub(x[bla[2*b]], x[bla[2*b+1]], dx);
 409             unitv(dx, r[b]);
 410         } /* 16 ncons flops */
 411     }
 412
 413 #pragma omp barrier
 414     for (b = b0; b < b1; b++)
 415     {
 416         tmp0 = r[b][0];
 417         tmp1 = r[b][1];
 418         tmp2 = r[b][2];
 419         i    = bla[2*b];
 420         j    = bla[2*b+1];
 421         for (n = blnr[b]; n < blnr[b+1]; n++)
 422         {
 423             k       = blbnb[n];
 424             blcc[n] = blmf[n]*(tmp0*r[k][0] + tmp1*r[k][1] + tmp2*r[k][2]);
 425         } /* 6 nr flops */
 426         mvb = blc[b]*(tmp0*(f[i][0] - f[j][0]) +
 427                       tmp1*(f[i][1] - f[j][1]) +
 428                       tmp2*(f[i][2] - f[j][2]));
 429         rhs1[b] = mvb;
 430         sol[b]  = mvb;
 431         /* 7 flops */
 432     }
 433     /* Together: 23*ncons + 6*nrtot flops */
 434
 435     lincs_matrix_expand(lincsd, b0, b1, blcc, rhs1, rhs2, sol);
 436     /* nrec*(ncons+2*nrtot) flops */
 437
 438     if (econq == econqDeriv_FlexCon)
 439     {
 440         /* We only want to constraint the flexible constraints,
 441          * so we mask out the normal ones by setting sol to 0.
 442          */
 443         for (b = b0; b < b1; b++)
 444         {
 445             if (!(lincsd->bllen0[b] == 0 && lincsd->ddist[b] == 0))
 446             {
 447                 sol[b] = 0;
 448             }
 449         }
 450     }
 451
 452     /* We multiply sol by blc, so we can use lincs_update_atoms for OpenMP */
 453     for (b = b0; b < b1; b++)
 454     {
 455         sol[b] *= blc[b];
 456     }
 457
 458     /* When constraining forces, we should not use mass weighting,
 459      * so we pass invmass=NULL, which results in the use of 1 for all atoms.
 460      */
 461     lincs_update_atoms(lincsd, th, 1.0, sol, r,
 462                        (econq != econqForce) ? invmass : NULL, fp);
 463
 464     if (dvdlambda != NULL)
 465     {
 466 #pragma omp barrier
 467         for (b = b0; b < b1; b++)
 468         {
 469             *dvdlambda -= sol[b]*lincsd->ddist[b];
 470         }
 471         /* 10 ncons flops */
 472     }
 473
 474     if (bCalcVir)
 475     {
 476         /* Constraint virial,
 477          * determines sum r_bond x delta f,
 478          * where delta f is the constraint correction
 479          * of the quantity that is being constrained.
 480          */
 481         for (b = b0; b < b1; b++)
 482         {
 483             mvb = lincsd->bllen[b]*sol[b];
 484             for (i = 0; i < DIM; i++)
 485             {
 486                 tmp1 = mvb*r[b][i];
 487                 for (j = 0; j < DIM; j++)
 488                 {
 489                     rmdf[i][j] += tmp1*r[b][j];
 490                 }
 491             }
 492         } /* 23 ncons flops */
 493     }
 494 }
 495
 496 static void do_lincs(rvec *x, rvec *xp, matrix box, t_pbc *pbc,
 497                      struct gmx_lincsdata *lincsd, int th,
 498                      real *invmass,
 499                      t_commrec *cr,
 500                      gmx_bool bCalcLambda,
 501                      real wangle, int *warn,
 502                      real invdt, rvec *v,
 503                      gmx_bool bCalcVir, tensor vir_r_m_dr)
 504 {
 505     int      b0, b1, b, i, j, k, n, iter;
 506     real     tmp0, tmp1, tmp2, im1, im2, mvb, rlen, len, len2, dlen2, wfac;
 507     rvec     dx;
 508     int     *bla, *blnr, *blbnb;
 509     rvec    *r;
 510     real    *blc, *blmf, *bllen, *blcc, *rhs1, *rhs2, *sol, *blc_sol, *mlambda;
 511     int     *nlocat;
 512
 513     b0 = lincsd->th[th].b0;
 514     b1 = lincsd->th[th].b1;
 515
 516     bla     = lincsd->bla;
 517     r       = lincsd->tmpv;
 518     blnr    = lincsd->blnr;
 519     blbnb   = lincsd->blbnb;
 520     blc     = lincsd->blc;
 521     blmf    = lincsd->blmf;
 522     bllen   = lincsd->bllen;
 523     blcc    = lincsd->tmpncc;
 524     rhs1    = lincsd->tmp1;
 525     rhs2    = lincsd->tmp2;
 526     sol     = lincsd->tmp3;
 527     blc_sol = lincsd->tmp4;
 528     mlambda = lincsd->mlambda;
 529
 530     if (DOMAINDECOMP(cr) && cr->dd->constraints)
 531     {
 532         nlocat = dd_constraints_nlocalatoms(cr->dd);
 533     }
 534     else
 535     {
 536         nlocat = NULL;
 537     }
 538
 539     if (pbc)
 540     {
 541         /* Compute normalized i-j vectors */
 542         for (b = b0; b < b1; b++)
 543         {
 544             pbc_dx_aiuc(pbc, x[bla[2*b]], x[bla[2*b+1]], dx);
 545             unitv(dx, r[b]);
 546         }
 547 #pragma omp barrier
 548         for (b = b0; b < b1; b++)
 549         {
 550             for (n = blnr[b]; n < blnr[b+1]; n++)
 551             {
 552                 blcc[n] = blmf[n]*iprod(r[b], r[blbnb[n]]);
 553             }
 554             pbc_dx_aiuc(pbc, xp[bla[2*b]], xp[bla[2*b+1]], dx);
 555             mvb     = blc[b]*(iprod(r[b], dx) - bllen[b]);
 556             rhs1[b] = mvb;
 557             sol[b]  = mvb;
 558         }
 559     }
 560     else
 561     {
 562         /* Compute normalized i-j vectors */
 563         for (b = b0; b < b1; b++)
 564         {
 565             i       = bla[2*b];
 566             j       = bla[2*b+1];
 567             tmp0    = x[i][0] - x[j][0];
 568             tmp1    = x[i][1] - x[j][1];
 569             tmp2    = x[i][2] - x[j][2];
 570             rlen    = gmx_invsqrt(tmp0*tmp0+tmp1*tmp1+tmp2*tmp2);
 571             r[b][0] = rlen*tmp0;
 572             r[b][1] = rlen*tmp1;
 573             r[b][2] = rlen*tmp2;
 574         } /* 16 ncons flops */
 575
 576 #pragma omp barrier
 577         for (b = b0; b < b1; b++)
 578         {
 579             tmp0 = r[b][0];
 580             tmp1 = r[b][1];
 581             tmp2 = r[b][2];
 582             len  = bllen[b];
 583             i    = bla[2*b];
 584             j    = bla[2*b+1];
 585             for (n = blnr[b]; n < blnr[b+1]; n++)
 586             {
 587                 k       = blbnb[n];
 588                 blcc[n] = blmf[n]*(tmp0*r[k][0] + tmp1*r[k][1] + tmp2*r[k][2]);
 589             } /* 6 nr flops */
 590             mvb = blc[b]*(tmp0*(xp[i][0] - xp[j][0]) +
 591                           tmp1*(xp[i][1] - xp[j][1]) +
 592                           tmp2*(xp[i][2] - xp[j][2]) - len);
 593             rhs1[b] = mvb;
 594             sol[b]  = mvb;
 595             /* 10 flops */
 596         }
 597         /* Together: 26*ncons + 6*nrtot flops */
 598     }
 599
 600     lincs_matrix_expand(lincsd, b0, b1, blcc, rhs1, rhs2, sol);
 601     /* nrec*(ncons+2*nrtot) flops */
 602
 603     for (b = b0; b < b1; b++)
 604     {
 605         mlambda[b] = blc[b]*sol[b];
 606     }
 607
 608     /* Update the coordinates */
 609     lincs_update_atoms(lincsd, th, 1.0, mlambda, r, invmass, xp);
 610
 611     /*
 612      ********  Correction for centripetal effects  ********
 613      */
 614
 615     wfac = cos(DEG2RAD*wangle);
 616     wfac = wfac*wfac;
 617
 618     for (iter = 0; iter < lincsd->nIter; iter++)
 619     {
 620         if ((lincsd->bCommIter && DOMAINDECOMP(cr) && cr->dd->constraints))
 621         {
 622 #pragma omp barrier
 623 #pragma omp master
 624             {
 625                 /* Communicate the corrected non-local coordinates */
 626                 if (DOMAINDECOMP(cr))
 627                 {
 628                     dd_move_x_constraints(cr->dd, box, xp, NULL);
 629                 }
 630             }
 631         }
 632
 633 #pragma omp barrier
 634         for (b = b0; b < b1; b++)
 635         {
 636             len = bllen[b];
 637             if (pbc)
 638             {
 639                 pbc_dx_aiuc(pbc, xp[bla[2*b]], xp[bla[2*b+1]], dx);
 640             }
 641             else
 642             {
 643                 rvec_sub(xp[bla[2*b]], xp[bla[2*b+1]], dx);
 644             }
 645             len2  = len*len;
 646             dlen2 = 2*len2 - norm2(dx);
 647             if (dlen2 < wfac*len2 && (nlocat == NULL || nlocat[b]))
 648             {
 649                 *warn = b;
 650             }
 651             if (dlen2 > 0)
 652             {
 653                 mvb = blc[b]*(len - dlen2*gmx_invsqrt(dlen2));
 654             }
 655             else
 656             {
 657                 mvb = blc[b]*len;
 658             }
 659             rhs1[b] = mvb;
 660             sol[b]  = mvb;
 661         } /* 20*ncons flops */
 662
 663         lincs_matrix_expand(lincsd, b0, b1, blcc, rhs1, rhs2, sol);
 664         /* nrec*(ncons+2*nrtot) flops */
 665
 666         for (b = b0; b < b1; b++)
 667         {
 668             mvb         = blc[b]*sol[b];
 669             blc_sol[b]  = mvb;
 670             mlambda[b] += mvb;
 671         }
 672
 673         /* Update the coordinates */
 674         lincs_update_atoms(lincsd, th, 1.0, blc_sol, r, invmass, xp);
 675     }
 676     /* nit*ncons*(37+9*nrec) flops */
 677
 678     if (v != NULL)
 679     {
 680         /* Update the velocities */
 681         lincs_update_atoms(lincsd, th, invdt, mlambda, r, invmass, v);
 682         /* 16 ncons flops */
 683     }
 684
 685     if (nlocat != NULL && bCalcLambda)
 686     {
 687         /* In lincs_update_atoms thread might cross-read mlambda */
 688 #pragma omp barrier
 689
 690         /* Only account for local atoms */
 691         for (b = b0; b < b1; b++)
 692         {
 693             mlambda[b] *= 0.5*nlocat[b];
 694         }
 695     }
 696
 697     if (bCalcVir)
 698     {
 699         /* Constraint virial */
 700         for (b = b0; b < b1; b++)
 701         {
 702             tmp0 = -bllen[b]*mlambda[b];
 703             for (i = 0; i < DIM; i++)
 704             {
 705                 tmp1 = tmp0*r[b][i];
 706                 for (j = 0; j < DIM; j++)
 707                 {
 708                     vir_r_m_dr[i][j] -= tmp1*r[b][j];
 709                 }
 710             }
 711         } /* 22 ncons flops */
 712     }
 713
 714     /* Total:
 715      * 26*ncons + 6*nrtot + nrec*(ncons+2*nrtot)
 716      * + nit * (20*ncons + nrec*(ncons+2*nrtot) + 17 ncons)
 717      *
 718      * (26+nrec)*ncons + (6+2*nrec)*nrtot
 719      * + nit * ((37+nrec)*ncons + 2*nrec*nrtot)
 720      * if nit=1
 721      * (63+nrec)*ncons + (6+4*nrec)*nrtot
 722      */
 723 }
 724
 725 void set_lincs_matrix(struct gmx_lincsdata *li, real *invmass, real lambda)
 726 {
 727     int        i, a1, a2, n, k, sign, center;
 728     int        end, nk, kk;
 729     const real invsqrt2 = 0.7071067811865475244;
 730
 731     for (i = 0; (i < li->nc); i++)
 732     {
 733         a1          = li->bla[2*i];
 734         a2          = li->bla[2*i+1];
 735         li->blc[i]  = gmx_invsqrt(invmass[a1] + invmass[a2]);
 736         li->blc1[i] = invsqrt2;
 737     }
 738
 739     /* Construct the coupling coefficient matrix blmf */
 740     li->ntriangle    = 0;
 741     li->ncc_triangle = 0;
 742     for (i = 0; (i < li->nc); i++)
 743     {
 744         a1 = li->bla[2*i];
 745         a2 = li->bla[2*i+1];
 746         for (n = li->blnr[i]; (n < li->blnr[i+1]); n++)
 747         {
 748             k = li->blbnb[n];
 749             if (a1 == li->bla[2*k] || a2 == li->bla[2*k+1])
 750             {
 751                 sign = -1;
 752             }
 753             else
 754             {
 755                 sign = 1;
 756             }
 757             if (a1 == li->bla[2*k] || a1 == li->bla[2*k+1])
 758             {
 759                 center = a1;
 760                 end    = a2;
 761             }
 762             else
 763             {
 764                 center = a2;
 765                 end    = a1;
 766             }
 767             li->blmf[n]  = sign*invmass[center]*li->blc[i]*li->blc[k];
 768             li->blmf1[n] = sign*0.5;
 769             if (li->ncg_triangle > 0)
 770             {
 771                 /* Look for constraint triangles */
 772                 for (nk = li->blnr[k]; (nk < li->blnr[k+1]); nk++)
 773                 {
 774                     kk = li->blbnb[nk];
 775                     if (kk != i && kk != k &&
 776                         (li->bla[2*kk] == end || li->bla[2*kk+1] == end))
 777                     {
 778                         if (li->ntriangle == 0 ||
 779                             li->triangle[li->ntriangle-1] < i)
 780                         {
 781                             /* Add this constraint to the triangle list */
 782                             li->triangle[li->ntriangle] = i;
 783                             li->tri_bits[li->ntriangle] = 0;
 784                             li->ntriangle++;
 785                             if (li->blnr[i+1] - li->blnr[i] > sizeof(li->tri_bits[0])*8 - 1)
 786                             {
 787                                 gmx_fatal(FARGS, "A constraint is connected to %d constraints, this is more than the %d allowed for constraints participating in triangles",
 788                                           li->blnr[i+1] - li->blnr[i],
 789                                           sizeof(li->tri_bits[0])*8-1);
 790                             }
 791                         }
 792                         li->tri_bits[li->ntriangle-1] |= (1<<(n-li->blnr[i]));
 793                         li->ncc_triangle++;
 794                     }
 795                 }
 796             }
 797         }
 798     }
 799
 800     if (debug)
 801     {
 802         fprintf(debug, "Of the %d constraints %d participate in triangles\n",
 803                 li->nc, li->ntriangle);
 804         fprintf(debug, "There are %d couplings of which %d in triangles\n",
 805                 li->ncc, li->ncc_triangle);
 806     }
 807
 808     /* Set matlam,
 809      * so we know with which lambda value the masses have been set.
 810      */
 811     li->matlam = lambda;
 812 }
 813
 814 static int count_triangle_constraints(t_ilist *ilist, t_blocka *at2con)
 815 {
 816     int      ncon1, ncon_tot;
 817     int      c0, a00, a01, n1, c1, a10, a11, ac1, n2, c2, a20, a21;
 818     int      ncon_triangle;
 819     gmx_bool bTriangle;
 820     t_iatom *ia1, *ia2, *iap;
 821
 822     ncon1    = ilist[F_CONSTR].nr/3;
 823     ncon_tot = ncon1 + ilist[F_CONSTRNC].nr/3;
 824
 825     ia1 = ilist[F_CONSTR].iatoms;
 826     ia2 = ilist[F_CONSTRNC].iatoms;
 827
 828     ncon_triangle = 0;
 829     for (c0 = 0; c0 < ncon_tot; c0++)
 830     {
 831         bTriangle = FALSE;
 832         iap       = constr_iatomptr(ncon1, ia1, ia2, c0);
 833         a00       = iap[1];
 834         a01       = iap[2];
 835         for (n1 = at2con->index[a01]; n1 < at2con->index[a01+1]; n1++)
 836         {
 837             c1 = at2con->a[n1];
 838             if (c1 != c0)
 839             {
 840                 iap = constr_iatomptr(ncon1, ia1, ia2, c1);
 841                 a10 = iap[1];
 842                 a11 = iap[2];
 843                 if (a10 == a01)
 844                 {
 845                     ac1 = a11;
 846                 }
 847                 else
 848                 {
 849                     ac1 = a10;
 850                 }
 851                 for (n2 = at2con->index[ac1]; n2 < at2con->index[ac1+1]; n2++)
 852                 {
 853                     c2 = at2con->a[n2];
 854                     if (c2 != c0 && c2 != c1)
 855                     {
 856                         iap = constr_iatomptr(ncon1, ia1, ia2, c2);
 857                         a20 = iap[1];
 858                         a21 = iap[2];
 859                         if (a20 == a00 || a21 == a00)
 860                         {
 861                             bTriangle = TRUE;
 862                         }
 863                     }
 864                 }
 865             }
 866         }
 867         if (bTriangle)
 868         {
 869             ncon_triangle++;
 870         }
 871     }
 872
 873     return ncon_triangle;
 874 }
 875
 876 static gmx_bool more_than_two_sequential_constraints(const t_ilist  *ilist,
 877                                                      const t_blocka *at2con)
 878 {
 879     t_iatom  *ia1, *ia2, *iap;
 880     int       ncon1, ncon_tot, c;
 881     int       a1, a2;
 882     gmx_bool  bMoreThanTwoSequentialConstraints;
 883
 884     ncon1    = ilist[F_CONSTR].nr/3;
 885     ncon_tot = ncon1 + ilist[F_CONSTRNC].nr/3;
 886
 887     ia1 = ilist[F_CONSTR].iatoms;
 888     ia2 = ilist[F_CONSTRNC].iatoms;
 889
 890     bMoreThanTwoSequentialConstraints = FALSE;
 891     for (c = 0; c < ncon_tot && !bMoreThanTwoSequentialConstraints; c++)
 892     {
 893         iap = constr_iatomptr(ncon1, ia1, ia2, c);
 894         a1  = iap[1];
 895         a2  = iap[2];
 896         /* Check if this constraint has constraints connected at both atoms */
 897         if (at2con->index[a1+1] - at2con->index[a1] > 1 &&
 898             at2con->index[a2+1] - at2con->index[a2] > 1)
 899         {
 900             bMoreThanTwoSequentialConstraints = TRUE;
 901         }
 902     }
 903
 904     return bMoreThanTwoSequentialConstraints;
 905 }
 906
 907 static int int_comp(const void *a, const void *b)
 908 {
 909     return (*(int *)a) - (*(int *)b);
 910 }
 911
 912 gmx_lincsdata_t init_lincs(FILE *fplog, gmx_mtop_t *mtop,
 913                            int nflexcon_global, t_blocka *at2con,
 914                            gmx_bool bPLINCS, int nIter, int nProjOrder)
 915 {
 916     struct gmx_lincsdata *li;
 917     int                   mb;
 918     gmx_moltype_t        *molt;
 919
 920     if (fplog)
 921     {
 922         fprintf(fplog, "\nInitializing%s LINear Constraint Solver\n",
 923                 bPLINCS ? " Parallel" : "");
 924     }
 925
 926     snew(li, 1);
 927
 928     li->ncg      =
 929         gmx_mtop_ftype_count(mtop, F_CONSTR) +
 930         gmx_mtop_ftype_count(mtop, F_CONSTRNC);
 931     li->ncg_flex = nflexcon_global;
 932
 933     li->nIter  = nIter;
 934     li->nOrder = nProjOrder;
 935
 936     li->ncg_triangle = 0;
 937     li->bCommIter    = FALSE;
 938     for (mb = 0; mb < mtop->nmolblock; mb++)
 939     {
 940         molt              = &mtop->moltype[mtop->molblock[mb].type];
 941         li->ncg_triangle +=
 942             mtop->molblock[mb].nmol*
 943             count_triangle_constraints(molt->ilist,
 944                                        &at2con[mtop->molblock[mb].type]);
 945         if (bPLINCS && li->bCommIter == FALSE)
 946         {
 947             /* Check if we need to communicate not only before LINCS,
 948              * but also before each iteration.
 949              * The check for only two sequential constraints is only
 950              * useful for the common case of H-bond only constraints.
 951              * With more effort we could also make it useful for small
 952              * molecules with nr. sequential constraints <= nOrder-1.
 953              */
 954             li->bCommIter = (li->nOrder < 1 || more_than_two_sequential_constraints(molt->ilist, &at2con[mtop->molblock[mb].type]));
 955         }
 956     }
 957     if (debug && bPLINCS)
 958     {
 959         fprintf(debug, "PLINCS communication before each iteration: %d\n",
 960                 li->bCommIter);
 961     }
 962
 963     /* LINCS can run on any number of threads.
 964      * Currently the number is fixed for the whole simulation,
 965      * but it could be set in set_lincs().
 966      */
 967     li->nth = gmx_omp_nthreads_get(emntLINCS);
 968     if (li->nth == 1)
 969     {
 970         snew(li->th, 1);
 971     }
 972     else
 973     {
 974         /* Allocate an extra elements for "thread-overlap" constraints */
 975         snew(li->th, li->nth+1);
 976     }
 977     if (debug)
 978     {
 979         fprintf(debug, "LINCS: using %d threads\n", li->nth);
 980     }
 981
 982     if (bPLINCS || li->ncg_triangle > 0)
 983     {
 984         please_cite(fplog, "Hess2008a");
 985     }
 986     else
 987     {
 988         please_cite(fplog, "Hess97a");
 989     }
 990
 991     if (fplog)
 992     {
 993         fprintf(fplog, "The number of constraints is %d\n", li->ncg);
 994         if (bPLINCS)
 995         {
 996             fprintf(fplog, "There are inter charge-group constraints,\n"
 997                     "will communicate selected coordinates each lincs iteration\n");
 998         }
 999         if (li->ncg_triangle > 0)
1000         {
1001             fprintf(fplog,
1002                     "%d constraints are involved in constraint triangles,\n"
1003                     "will apply an additional matrix expansion of order %d for couplings\n"
1004                     "between constraints inside triangles\n",
1005                     li->ncg_triangle, li->nOrder);
1006         }
1007     }
1008
1009     return li;
1010 }
1011
1012 /* Sets up the work division over the threads */
1013 static void lincs_thread_setup(struct gmx_lincsdata *li, int natoms)
1014 {
1015     lincs_thread_t *li_m;
1016     int             th;
1017     unsigned       *atf;
1018     int             a;
1019
1020     if (natoms > li->atf_nalloc)
1021     {
1022         li->atf_nalloc = over_alloc_large(natoms);
1023         srenew(li->atf, li->atf_nalloc);
1024     }
1025
1026     atf = li->atf;
1027     /* Clear the atom flags */
1028     for (a = 0; a < natoms; a++)
1029     {
1030         atf[a] = 0;
1031     }
1032
1033     for (th = 0; th < li->nth; th++)
1034     {
1035         lincs_thread_t *li_th;
1036         int             b;
1037
1038         li_th = &li->th[th];
1039
1040         /* The constraints are divided equally over the threads */
1041         li_th->b0 = (li->nc* th   )/li->nth;
1042         li_th->b1 = (li->nc*(th+1))/li->nth;
1043
1044         if (th < sizeof(*atf)*8)
1045         {
1046             /* For each atom set a flag for constraints from each */
1047             for (b = li_th->b0; b < li_th->b1; b++)
1048             {
1049                 atf[li->bla[b*2]  ] |= (1U<<th);
1050                 atf[li->bla[b*2+1]] |= (1U<<th);
1051             }
1052         }
1053     }
1054
1055 #pragma omp parallel for num_threads(li->nth) schedule(static)
1056     for (th = 0; th < li->nth; th++)
1057     {
1058         lincs_thread_t *li_th;
1059         unsigned        mask;
1060         int             b;
1061
1062         li_th = &li->th[th];
1063
1064         if (li_th->b1 - li_th->b0 > li_th->ind_nalloc)
1065         {
1066             li_th->ind_nalloc = over_alloc_large(li_th->b1-li_th->b0);
1067             srenew(li_th->ind, li_th->ind_nalloc);
1068             srenew(li_th->ind_r, li_th->ind_nalloc);
1069         }
1070
1071         if (th < sizeof(*atf)*8)
1072         {
1073             mask = (1U<<th) - 1U;
1074
1075             li_th->nind   = 0;
1076             li_th->nind_r = 0;
1077             for (b = li_th->b0; b < li_th->b1; b++)
1078             {
1079                 /* We let the constraint with the lowest thread index
1080                  * operate on atoms with constraints from multiple threads.
1081                  */
1082                 if (((atf[li->bla[b*2]]   & mask) == 0) &&
1083                     ((atf[li->bla[b*2+1]] & mask) == 0))
1084                 {
1085                     /* Add the constraint to the local atom update index */
1086                     li_th->ind[li_th->nind++] = b;
1087                 }
1088                 else
1089                 {
1090                     /* Add the constraint to the rest block */
1091                     li_th->ind_r[li_th->nind_r++] = b;
1092                 }
1093             }
1094         }
1095         else
1096         {
1097             /* We are out of bits, assign all constraints to rest */
1098             for (b = li_th->b0; b < li_th->b1; b++)
1099             {
1100                 li_th->ind_r[li_th->nind_r++] = b;
1101             }
1102         }
1103     }
1104
1105     /* We need to copy all constraints which have not be assigned
1106      * to a thread to a separate list which will be handled by one thread.
1107      */
1108     li_m = &li->th[li->nth];
1109
1110     li_m->nind = 0;
1111     for (th = 0; th < li->nth; th++)
1112     {
1113         lincs_thread_t *li_th;
1114         int             b;
1115
1116         li_th   = &li->th[th];
1117
1118         if (li_m->nind + li_th->nind_r > li_m->ind_nalloc)
1119         {
1120             li_m->ind_nalloc = over_alloc_large(li_m->nind+li_th->nind_r);
1121             srenew(li_m->ind, li_m->ind_nalloc);
1122         }
1123
1124         for (b = 0; b < li_th->nind_r; b++)
1125         {
1126             li_m->ind[li_m->nind++] = li_th->ind_r[b];
1127         }
1128
1129         if (debug)
1130         {
1131             fprintf(debug, "LINCS thread %d: %d constraints\n",
1132                     th, li_th->nind);
1133         }
1134     }
1135
1136     if (debug)
1137     {
1138         fprintf(debug, "LINCS thread r: %d constraints\n",
1139                 li_m->nind);
1140     }
1141 }
1142
1143
1144 void set_lincs(t_idef *idef, t_mdatoms *md,
1145                gmx_bool bDynamics, t_commrec *cr,
1146                struct gmx_lincsdata *li)
1147 {
1148     int          start, natoms, nflexcon;
1149     t_blocka     at2con;
1150     t_iatom     *iatom;
1151     int          i, k, ncc_alloc, ni, con, nconnect, concon;
1152     int          type, a1, a2;
1153     real         lenA = 0, lenB;
1154     gmx_bool     bLocal;
1155
1156     li->nc  = 0;
1157     li->ncc = 0;
1158     /* Zero the thread index ranges.
1159      * Otherwise without local constraints we could return with old ranges.
1160      */
1161     for (i = 0; i < li->nth; i++)
1162     {
1163         li->th[i].b0   = 0;
1164         li->th[i].b1   = 0;
1165         li->th[i].nind = 0;
1166     }
1167     if (li->nth > 1)
1168     {
1169         li->th[li->nth].nind = 0;
1170     }
1171
1172     /* This is the local topology, so there are only F_CONSTR constraints */
1173     if (idef->il[F_CONSTR].nr == 0)
1174     {
1175         /* There are no constraints,
1176          * we do not need to fill any data structures.
1177          */
1178         return;
1179     }
1180
1181     if (debug)
1182     {
1183         fprintf(debug, "Building the LINCS connectivity\n");
1184     }
1185
1186     if (DOMAINDECOMP(cr))
1187     {
1188         if (cr->dd->constraints)
1189         {
1190             dd_get_constraint_range(cr->dd, &start, &natoms);
1191         }
1192         else
1193         {
1194             natoms = cr->dd->nat_home;
1195         }
1196         start = 0;
1197     }
1198     else
1199     {
1200         start  = 0;
1201         natoms = md->homenr;
1202     }
1203     at2con = make_at2con(start, natoms, idef->il, idef->iparams, bDynamics,
1204                          &nflexcon);
1205
1206
1207     if (idef->il[F_CONSTR].nr/3 > li->nc_alloc || li->nc_alloc == 0)
1208     {
1209         li->nc_alloc = over_alloc_dd(idef->il[F_CONSTR].nr/3);
1210         srenew(li->bllen0, li->nc_alloc);
1211         srenew(li->ddist, li->nc_alloc);
1212         srenew(li->bla, 2*li->nc_alloc);
1213         srenew(li->blc, li->nc_alloc);
1214         srenew(li->blc1, li->nc_alloc);
1215         srenew(li->blnr, li->nc_alloc+1);
1216         srenew(li->bllen, li->nc_alloc);
1217         srenew(li->tmpv, li->nc_alloc);
1218         srenew(li->tmp1, li->nc_alloc);
1219         srenew(li->tmp2, li->nc_alloc);
1220         srenew(li->tmp3, li->nc_alloc);
1221         srenew(li->tmp4, li->nc_alloc);
1222         srenew(li->mlambda, li->nc_alloc);
1223         if (li->ncg_triangle > 0)
1224         {
1225             /* This is allocating too much, but it is difficult to improve */
1226             srenew(li->triangle, li->nc_alloc);
1227             srenew(li->tri_bits, li->nc_alloc);
1228         }
1229     }
1230
1231     iatom = idef->il[F_CONSTR].iatoms;
1232
1233     ncc_alloc   = li->ncc_alloc;
1234     li->blnr[0] = 0;
1235
1236     ni = idef->il[F_CONSTR].nr/3;
1237
1238     con           = 0;
1239     nconnect      = 0;
1240     li->blnr[con] = nconnect;
1241     for (i = 0; i < ni; i++)
1242     {
1243         bLocal = TRUE;
1244         type   = iatom[3*i];
1245         a1     = iatom[3*i+1];
1246         a2     = iatom[3*i+2];
1247         lenA   = idef->iparams[type].constr.dA;
1248         lenB   = idef->iparams[type].constr.dB;
1249         /* Skip the flexible constraints when not doing dynamics */
1250         if (bDynamics || lenA != 0 || lenB != 0)
1251         {
1252             li->bllen0[con]  = lenA;
1253             li->ddist[con]   = lenB - lenA;
1254             /* Set the length to the topology A length */
1255             li->bllen[con]   = li->bllen0[con];
1256             li->bla[2*con]   = a1;
1257             li->bla[2*con+1] = a2;
1258             /* Construct the constraint connection matrix blbnb */
1259             for (k = at2con.index[a1-start]; k < at2con.index[a1-start+1]; k++)
1260             {
1261                 concon = at2con.a[k];
1262                 if (concon != i)
1263                 {
1264                     if (nconnect >= ncc_alloc)
1265                     {
1266                         ncc_alloc = over_alloc_small(nconnect+1);
1267                         srenew(li->blbnb, ncc_alloc);
1268                     }
1269                     li->blbnb[nconnect++] = concon;
1270                 }
1271             }
1272             for (k = at2con.index[a2-start]; k < at2con.index[a2-start+1]; k++)
1273             {
1274                 concon = at2con.a[k];
1275                 if (concon != i)
1276                 {
1277                     if (nconnect+1 > ncc_alloc)
1278                     {
1279                         ncc_alloc = over_alloc_small(nconnect+1);
1280                         srenew(li->blbnb, ncc_alloc);
1281                     }
1282                     li->blbnb[nconnect++] = concon;
1283                 }
1284             }
1285             li->blnr[con+1] = nconnect;
1286
1287             if (cr->dd == NULL)
1288             {
1289                 /* Order the blbnb matrix to optimize memory access */
1290                 qsort(&(li->blbnb[li->blnr[con]]), li->blnr[con+1]-li->blnr[con],
1291                       sizeof(li->blbnb[0]), int_comp);
1292             }
1293             /* Increase the constraint count */
1294             con++;
1295         }
1296     }
1297
1298     done_blocka(&at2con);
1299
1300     /* This is the real number of constraints,
1301      * without dynamics the flexible constraints are not present.
1302      */
1303     li->nc = con;
1304
1305     li->ncc = li->blnr[con];
1306     if (cr->dd == NULL)
1307     {
1308         /* Since the matrix is static, we can free some memory */
1309         ncc_alloc = li->ncc;
1310         srenew(li->blbnb, ncc_alloc);
1311     }
1312
1313     if (ncc_alloc > li->ncc_alloc)
1314     {
1315         li->ncc_alloc = ncc_alloc;
1316         srenew(li->blmf, li->ncc_alloc);
1317         srenew(li->blmf1, li->ncc_alloc);
1318         srenew(li->tmpncc, li->ncc_alloc);
1319     }
1320
1321     if (debug)
1322     {
1323         fprintf(debug, "Number of constraints is %d, couplings %d\n",
1324                 li->nc, li->ncc);
1325     }
1326
1327     if (li->nth == 1)
1328     {
1329         li->th[0].b0 = 0;
1330         li->th[0].b1 = li->nc;
1331     }
1332     else
1333     {
1334         lincs_thread_setup(li, md->nr);
1335     }
1336
1337     set_lincs_matrix(li, md->invmass, md->lambda);
1338 }
1339
1340 static void lincs_warning(FILE *fplog,
1341                           gmx_domdec_t *dd, rvec *x, rvec *xprime, t_pbc *pbc,
1342                           int ncons, int *bla, real *bllen, real wangle,
1343                           int maxwarn, int *warncount)
1344 {
1345     int  b, i, j;
1346     rvec v0, v1;
1347     real wfac, d0, d1, cosine;
1348     char buf[STRLEN];
1349
1350     wfac = cos(DEG2RAD*wangle);
1351
1352     sprintf(buf, "bonds that rotated more than %g degrees:\n"
1353             " atom 1 atom 2  angle  previous, current, constraint length\n",
1354             wangle);
1355     fprintf(stderr, "%s", buf);
1356     if (fplog)
1357     {
1358         fprintf(fplog, "%s", buf);
1359     }
1360
1361     for (b = 0; b < ncons; b++)
1362     {
1363         i = bla[2*b];
1364         j = bla[2*b+1];
1365         if (pbc)
1366         {
1367             pbc_dx_aiuc(pbc, x[i], x[j], v0);
1368             pbc_dx_aiuc(pbc, xprime[i], xprime[j], v1);
1369         }
1370         else
1371         {
1372             rvec_sub(x[i], x[j], v0);
1373             rvec_sub(xprime[i], xprime[j], v1);
1374         }
1375         d0     = norm(v0);
1376         d1     = norm(v1);
1377         cosine = iprod(v0, v1)/(d0*d1);
1378         if (cosine < wfac)
1379         {
1380             sprintf(buf, " %6d %6d  %5.1f  %8.4f %8.4f    %8.4f\n",
1381                     ddglatnr(dd, i), ddglatnr(dd, j),
1382                     RAD2DEG*acos(cosine), d0, d1, bllen[b]);
1383             fprintf(stderr, "%s", buf);
1384             if (fplog)
1385             {
1386                 fprintf(fplog, "%s", buf);
1387             }
1388             if (!gmx_isfinite(d1))
1389             {
1390                 gmx_fatal(FARGS, "Bond length not finite.");
1391             }
1392
1393             (*warncount)++;
1394         }
1395     }
1396     if (*warncount > maxwarn)
1397     {
1398         too_many_constraint_warnings(econtLINCS, *warncount);
1399     }
1400 }
1401
1402 static void cconerr(gmx_domdec_t *dd,
1403                     int ncons, int *bla, real *bllen, rvec *x, t_pbc *pbc,
1404                     real *ncons_loc, real *ssd, real *max, int *imax)
1405 {
1406     real       len, d, ma, ssd2, r2;
1407     int       *nlocat, count, b, im;
1408     rvec       dx;
1409
1410     if (dd && dd->constraints)
1411     {
1412         nlocat = dd_constraints_nlocalatoms(dd);
1413     }
1414     else
1415     {
1416         nlocat = 0;
1417     }
1418
1419     ma    = 0;
1420     ssd2  = 0;
1421     im    = 0;
1422     count = 0;
1423     for (b = 0; b < ncons; b++)
1424     {
1425         if (pbc)
1426         {
1427             pbc_dx_aiuc(pbc, x[bla[2*b]], x[bla[2*b+1]], dx);
1428         }
1429         else
1430         {
1431             rvec_sub(x[bla[2*b]], x[bla[2*b+1]], dx);
1432         }
1433         r2  = norm2(dx);
1434         len = r2*gmx_invsqrt(r2);
1435         d   = fabs(len/bllen[b]-1);
1436         if (d > ma && (nlocat == NULL || nlocat[b]))
1437         {
1438             ma = d;
1439             im = b;
1440         }
1441         if (nlocat == NULL)
1442         {
1443             ssd2 += d*d;
1444             count++;
1445         }
1446         else
1447         {
1448             ssd2  += nlocat[b]*d*d;
1449             count += nlocat[b];
1450         }
1451     }
1452
1453     *ncons_loc = (nlocat ? 0.5 : 1)*count;
1454     *ssd       = (nlocat ? 0.5 : 1)*ssd2;
1455     *max       = ma;
1456     *imax      = im;
1457 }
1458
1459 gmx_bool constrain_lincs(FILE *fplog, gmx_bool bLog, gmx_bool bEner,
1460                          t_inputrec *ir,
1461                          gmx_int64_t step,
1462                          struct gmx_lincsdata *lincsd, t_mdatoms *md,
1463                          t_commrec *cr,
1464                          rvec *x, rvec *xprime, rvec *min_proj,
1465                          matrix box, t_pbc *pbc,
1466                          real lambda, real *dvdlambda,
1467                          real invdt, rvec *v,
1468                          gmx_bool bCalcVir, tensor vir_r_m_dr,
1469                          int econq,
1470                          t_nrnb *nrnb,
1471                          int maxwarn, int *warncount)
1472 {
1473     char      buf[STRLEN], buf2[22], buf3[STRLEN];
1474     int       i, warn, p_imax, error;
1475     real      ncons_loc, p_ssd, p_max = 0;
1476     rvec      dx;
1477     gmx_bool  bOK;
1478
1479     bOK = TRUE;
1480
1481     if (lincsd->nc == 0 && cr->dd == NULL)
1482     {
1483         if (bLog || bEner)
1484         {
1485             lincsd->rmsd_data[0] = 0;
1486             if (ir->eI == eiSD2 && v == NULL)
1487             {
1488                 i = 2;
1489             }
1490             else
1491             {
1492                 i = 1;
1493             }
1494             lincsd->rmsd_data[i] = 0;
1495         }
1496
1497         return bOK;
1498     }
1499
1500     if (econq == econqCoord)
1501     {
1502         if (ir->efep != efepNO)
1503         {
1504             if (md->nMassPerturbed && lincsd->matlam != md->lambda)
1505             {
1506                 set_lincs_matrix(lincsd, md->invmass, md->lambda);
1507             }
1508
1509             for (i = 0; i < lincsd->nc; i++)
1510             {
1511                 lincsd->bllen[i] = lincsd->bllen0[i] + lambda*lincsd->ddist[i];
1512             }
1513         }
1514
1515         if (lincsd->ncg_flex)
1516         {
1517             /* Set the flexible constraint lengths to the old lengths */
1518             if (pbc != NULL)
1519             {
1520                 for (i = 0; i < lincsd->nc; i++)
1521                 {
1522                     if (lincsd->bllen[i] == 0)
1523                     {
1524                         pbc_dx_aiuc(pbc, x[lincsd->bla[2*i]], x[lincsd->bla[2*i+1]], dx);
1525                         lincsd->bllen[i] = norm(dx);
1526                     }
1527                 }
1528             }
1529             else
1530             {
1531                 for (i = 0; i < lincsd->nc; i++)
1532                 {
1533                     if (lincsd->bllen[i] == 0)
1534                     {
1535                         lincsd->bllen[i] =
1536                             sqrt(distance2(x[lincsd->bla[2*i]],
1537                                            x[lincsd->bla[2*i+1]]));
1538                     }
1539                 }
1540             }
1541         }
1542
1543         if (bLog && fplog)
1544         {
1545             cconerr(cr->dd, lincsd->nc, lincsd->bla, lincsd->bllen, xprime, pbc,
1546                     &ncons_loc, &p_ssd, &p_max, &p_imax);
1547         }
1548
1549         /* This warn var can be updated by multiple threads
1550          * at the same time. But as we only need to detect
1551          * if a warning occured or not, this is not an issue.
1552          */
1553         warn = -1;
1554
1555         /* The OpenMP parallel region of constrain_lincs for coords */
1556 #pragma omp parallel num_threads(lincsd->nth)
1557         {
1558             int th = gmx_omp_get_thread_num();
1559
1560             clear_mat(lincsd->th[th].vir_r_m_dr);
1561
1562             do_lincs(x, xprime, box, pbc, lincsd, th,
1563                      md->invmass, cr,
1564                      bCalcVir || (ir->efep != efepNO),
1565                      ir->LincsWarnAngle, &warn,
1566                      invdt, v, bCalcVir,
1567                      th == 0 ? vir_r_m_dr : lincsd->th[th].vir_r_m_dr);
1568         }
1569
1570         if (ir->efep != efepNO)
1571         {
1572             real dt_2, dvdl = 0;
1573
1574             dt_2 = 1.0/(ir->delta_t*ir->delta_t);
1575             for (i = 0; (i < lincsd->nc); i++)
1576             {
1577                 dvdl -= lincsd->mlambda[i]*dt_2*lincsd->ddist[i];
1578             }
1579             *dvdlambda += dvdl;
1580         }
1581
1582         if (bLog && fplog && lincsd->nc > 0)
1583         {
1584             fprintf(fplog, "   Rel. Constraint Deviation:  RMS         MAX     between atoms\n");
1585             fprintf(fplog, "       Before LINCS          %.6f    %.6f %6d %6d\n",
1586                     sqrt(p_ssd/ncons_loc), p_max,
1587                     ddglatnr(cr->dd, lincsd->bla[2*p_imax]),
1588                     ddglatnr(cr->dd, lincsd->bla[2*p_imax+1]));
1589         }
1590         if (bLog || bEner)
1591         {
1592             cconerr(cr->dd, lincsd->nc, lincsd->bla, lincsd->bllen, xprime, pbc,
1593                     &ncons_loc, &p_ssd, &p_max, &p_imax);
1594             /* Check if we are doing the second part of SD */
1595             if (ir->eI == eiSD2 && v == NULL)
1596             {
1597                 i = 2;
1598             }
1599             else
1600             {
1601                 i = 1;
1602             }
1603             lincsd->rmsd_data[0] = ncons_loc;
1604             lincsd->rmsd_data[i] = p_ssd;
1605         }
1606         else
1607         {
1608             lincsd->rmsd_data[0] = 0;
1609             lincsd->rmsd_data[1] = 0;
1610             lincsd->rmsd_data[2] = 0;
1611         }
1612         if (bLog && fplog && lincsd->nc > 0)
1613         {
1614             fprintf(fplog,
1615                     "        After LINCS          %.6f    %.6f %6d %6d\n\n",
1616                     sqrt(p_ssd/ncons_loc), p_max,
1617                     ddglatnr(cr->dd, lincsd->bla[2*p_imax]),
1618                     ddglatnr(cr->dd, lincsd->bla[2*p_imax+1]));
1619         }
1620
1621         if (warn >= 0)
1622         {
1623             if (maxwarn >= 0)
1624             {
1625                 cconerr(cr->dd, lincsd->nc, lincsd->bla, lincsd->bllen, xprime, pbc,
1626                         &ncons_loc, &p_ssd, &p_max, &p_imax);
1627                 if (MULTISIM(cr))
1628                 {
1629                     sprintf(buf3, " in simulation %d", cr->ms->sim);
1630                 }
1631                 else
1632                 {
1633                     buf3[0] = 0;
1634                 }
1635                 sprintf(buf, "\nStep %s, time %g (ps)  LINCS WARNING%s\n"
1636                         "relative constraint deviation after LINCS:\n"
1637                         "rms %.6f, max %.6f (between atoms %d and %d)\n",
1638                         gmx_step_str(step, buf2), ir->init_t+step*ir->delta_t,
1639                         buf3,
1640                         sqrt(p_ssd/ncons_loc), p_max,
1641                         ddglatnr(cr->dd, lincsd->bla[2*p_imax]),
1642                         ddglatnr(cr->dd, lincsd->bla[2*p_imax+1]));
1643                 if (fplog)
1644                 {
1645                     fprintf(fplog, "%s", buf);
1646                 }
1647                 fprintf(stderr, "%s", buf);
1648                 lincs_warning(fplog, cr->dd, x, xprime, pbc,
1649                               lincsd->nc, lincsd->bla, lincsd->bllen,
1650                               ir->LincsWarnAngle, maxwarn, warncount);
1651             }
1652             bOK = (p_max < 0.5);
1653         }
1654
1655         if (lincsd->ncg_flex)
1656         {
1657             for (i = 0; (i < lincsd->nc); i++)
1658             {
1659                 if (lincsd->bllen0[i] == 0 && lincsd->ddist[i] == 0)
1660                 {
1661                     lincsd->bllen[i] = 0;
1662                 }
1663             }
1664         }
1665     }
1666     else
1667     {
1668         /* The OpenMP parallel region of constrain_lincs for derivatives */
1669 #pragma omp parallel num_threads(lincsd->nth)
1670         {
1671             int th = gmx_omp_get_thread_num();
1672
1673             do_lincsp(x, xprime, min_proj, pbc, lincsd, th,
1674                       md->invmass, econq, ir->efep != efepNO ? dvdlambda : NULL,
1675                       bCalcVir, th == 0 ? vir_r_m_dr : lincsd->th[th].vir_r_m_dr);
1676         }
1677     }
1678
1679     if (bCalcVir && lincsd->nth > 1)
1680     {
1681         for (i = 1; i < lincsd->nth; i++)
1682         {
1683             m_add(vir_r_m_dr, lincsd->th[i].vir_r_m_dr, vir_r_m_dr);
1684         }
1685     }
1686
1687     /* count assuming nit=1 */
1688     inc_nrnb(nrnb, eNR_LINCS, lincsd->nc);
1689     inc_nrnb(nrnb, eNR_LINCSMAT, (2+lincsd->nOrder)*lincsd->ncc);
1690     if (lincsd->ntriangle > 0)
1691     {
1692         inc_nrnb(nrnb, eNR_LINCSMAT, lincsd->nOrder*lincsd->ncc_triangle);
1693     }
1694     if (v)
1695     {
1696         inc_nrnb(nrnb, eNR_CONSTR_V, lincsd->nc*2);
1697     }
1698     if (bCalcVir)
1699     {
1700         inc_nrnb(nrnb, eNR_CONSTR_VIR, lincsd->nc);
1701     }
1702
1703     return bOK;
1704 }