src/gromacs/mdlib/clincs.c

   1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
   2  *
   3  *
   4  *                This source code is part of
   5  *
   6  *                 G   R   O   M   A   C   S
   7  *
   8  *          GROningen MAchine for Chemical Simulations
   9  *
  10  *                        VERSION 3.2.0
  11  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  12  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  13  * Copyright (c) 2001-2004, The GROMACS development team,
  14  * check out http://www.gromacs.org for more information.
  15  *
  16  * This program is free software; you can redistribute it and/or
  17  * modify it under the terms of the GNU General Public License
  18  * as published by the Free Software Foundation; either version 2
  19  * of the License, or (at your option) any later version.
  20  *
  21  * If you want to redistribute modifications, please consider that
  22  * scientific software is very special. Version control is crucial -
  23  * bugs must be traceable. We will be happy to consider code for
  24  * inclusion in the official distribution, but derived work must not
  25  * be called official GROMACS. Details are found in the README & COPYING
  26  * files - if they are missing, get the official version at www.gromacs.org.
  27  *
  28  * To help us fund GROMACS development, we humbly ask that you cite
  29  * the papers on the package - you can find them in the top README file.
  30  *
  31  * For more info, check our website at http://www.gromacs.org
  32  *
  33  * And Hey:
  34  * GROwing Monsters And Cloning Shrimps
  35  */
  36 /* This file is completely threadsafe - keep it that way! */
  37 #ifdef HAVE_CONFIG_H
  38 #include <config.h>
  39 #endif
  40
  41 #include <math.h>
  42 #include "main.h"
  43 #include "constr.h"
  44 #include "copyrite.h"
  45 #include "physics.h"
  46 #include "vec.h"
  47 #include "pbc.h"
  48 #include "smalloc.h"
  49 #include "mdrun.h"
  50 #include "nrnb.h"
  51 #include "domdec.h"
  52 #include "partdec.h"
  53 #include "mtop_util.h"
  54 #include "gmxfio.h"
  55
  56 typedef struct gmx_lincsdata {
  57     int  ncg;         /* the global number of constraints */
  58     int  ncg_flex;    /* the global number of flexible constraints */
  59     int  ncg_triangle;/* the global number of constraints in triangles */
  60     int  nIter;       /* the number of iterations */
  61     int  nOrder;      /* the order of the matrix expansion */
  62     int  nc;          /* the number of constraints */
  63     int  nc_alloc;    /* the number we allocated memory for */
  64     int  ncc;         /* the number of constraint connections */
  65     int  ncc_alloc;   /* the number we allocated memory for */
  66     real matlam;      /* the FE lambda value used for filling blc and blmf */
  67     real *bllen0;     /* the reference distance in topology A */
  68     real *ddist;      /* the reference distance in top B - the r.d. in top A */
  69     int  *bla;        /* the atom pairs involved in the constraints */
  70     real *blc;        /* 1/sqrt(invmass1 + invmass2) */
  71     real *blc1;       /* as blc, but with all masses 1 */
  72     int  *blnr;       /* index into blbnb and blmf */
  73     int  *blbnb;      /* list of constraint connections */
  74     int  ntriangle;   /* the local number of constraints in triangles */
  75     int  *triangle;   /* the list of triangle constraints */
  76     int  *tri_bits;   /* the bits tell if the matrix element should be used */
  77     int  ncc_triangle;/* the number of constraint connections in triangles */
  78     real *blmf;       /* matrix of mass factors for constraint connections */
  79     real *blmf1;      /* as blmf, but with all masses 1 */
  80     real *bllen;      /* the reference bond length */
  81     /* arrays for temporary storage in the LINCS algorithm */
  82     rvec *tmpv;
  83     real *tmpncc;
  84     real *tmp1;
  85     real *tmp2;
  86     real *tmp3;
  87     real *lambda;  /* the Lagrange multipliers */
  88     /* storage for the constraint RMS relative deviation output */
  89     real rmsd_data[3];
  90 } t_gmx_lincsdata;
  91
  92 real *lincs_rmsd_data(struct gmx_lincsdata *lincsd)
  93 {
  94     return lincsd->rmsd_data;
  95 }
  96
  97 real lincs_rmsd(struct gmx_lincsdata *lincsd,gmx_bool bSD2)
  98 {
  99     if (lincsd->rmsd_data[0] > 0)
 100     {
 101         return sqrt(lincsd->rmsd_data[bSD2 ? 2 : 1]/lincsd->rmsd_data[0]);
 102     }
 103     else
 104     {
 105         return 0;
 106     }
 107 }
 108
 109 static void lincs_matrix_expand(const struct gmx_lincsdata *lincsd,
 110                                 const real *blcc,
 111                                 real *rhs1,real *rhs2,real *sol)
 112 {
 113     int  nrec,rec,ncons,b,j,n,nr0,nr1;
 114     real mvb,*swap;
 115     int  ntriangle,tb,bits;
 116     const int *blnr=lincsd->blnr,*blbnb=lincsd->blbnb;
 117     const int *triangle=lincsd->triangle,*tri_bits=lincsd->tri_bits;
 118
 119     ncons     = lincsd->nc;
 120     ntriangle = lincsd->ntriangle;
 121     nrec      = lincsd->nOrder;
 122
 123     for(rec=0; rec<nrec; rec++)
 124     {
 125         for(b=0; b<ncons; b++)
 126         {
 127             mvb = 0;
 128             for(n=blnr[b]; n<blnr[b+1]; n++)
 129             {
 130                 j = blbnb[n];
 131                 mvb = mvb + blcc[n]*rhs1[j];
 132             }
 133             rhs2[b] = mvb;
 134             sol[b]  = sol[b] + mvb;
 135         }
 136         swap = rhs1;
 137         rhs1 = rhs2;
 138         rhs2 = swap;
 139     } /* nrec*(ncons+2*nrtot) flops */
 140
 141     if (ntriangle > 0)
 142     {
 143         /* Perform an extra nrec recursions for only the constraints
 144          * involved in rigid triangles.
 145          * In this way their accuracy should come close to those of the other
 146          * constraints, since traingles of constraints can produce eigenvalues
 147          * around 0.7, while the effective eigenvalue for bond constraints
 148          * is around 0.4 (and 0.7*0.7=0.5).
 149          */
 150         /* We need to copy the temporary array, since only the elements
 151          * for constraints involved in triangles are updated
 152          * and then the pointers are swapped.
 153          */
 154         for(b=0; b<ncons; b++)
 155         {
 156             rhs2[b] = rhs1[b];
 157         }
 158         for(rec=0; rec<nrec; rec++)
 159         {
 160             for(tb=0; tb<ntriangle; tb++)
 161             {
 162                 b    = triangle[tb];
 163                 bits = tri_bits[tb];
 164                 mvb = 0;
 165                 nr0 = blnr[b];
 166                 nr1 = blnr[b+1];
 167                 for(n=nr0; n<nr1; n++)
 168                 {
 169                     if (bits & (1<<(n-nr0)))
 170                     {
 171                         j = blbnb[n];
 172                         mvb = mvb + blcc[n]*rhs1[j];
 173                     }
 174                 }
 175                 rhs2[b] = mvb;
 176                 sol[b]  = sol[b] + mvb;
 177             }
 178             swap = rhs1;
 179             rhs1 = rhs2;
 180             rhs2 = swap;
 181         } /* flops count is missing here */
 182     }
 183 }
 184
 185 static void do_lincsp(rvec *x,rvec *f,rvec *fp,t_pbc *pbc,
 186                       struct gmx_lincsdata *lincsd,real *invmass,
 187                       int econq,real *dvdlambda,
 188                       gmx_bool bCalcVir,tensor rmdf)
 189 {
 190     int     b,i,j,k,n;
 191     real    tmp0,tmp1,tmp2,im1,im2,mvb,rlen,len,wfac,lam;
 192     rvec    dx;
 193     int     ncons,*bla,*blnr,*blbnb;
 194     rvec    *r;
 195     real    *blc,*blmf,*blcc,*rhs1,*rhs2,*sol;
 196
 197     ncons  = lincsd->nc;
 198     bla    = lincsd->bla;
 199     r      = lincsd->tmpv;
 200     blnr   = lincsd->blnr;
 201     blbnb  = lincsd->blbnb;
 202     if (econq != econqForce)
 203     {
 204         /* Use mass-weighted parameters */
 205         blc  = lincsd->blc;
 206         blmf = lincsd->blmf;
 207     }
 208     else
 209     {
 210         /* Use non mass-weighted parameters */
 211         blc  = lincsd->blc1;
 212         blmf = lincsd->blmf1;
 213     }
 214     blcc   = lincsd->tmpncc;
 215     rhs1   = lincsd->tmp1;
 216     rhs2   = lincsd->tmp2;
 217     sol    = lincsd->tmp3;
 218
 219     if (econq != econqForce)
 220     {
 221         dvdlambda = NULL;
 222     }
 223
 224     /* Compute normalized i-j vectors */
 225     if (pbc)
 226     {
 227         for(b=0; b<ncons; b++)
 228         {
 229             pbc_dx_aiuc(pbc,x[bla[2*b]],x[bla[2*b+1]],dx);
 230             unitv(dx,r[b]);
 231         }
 232     }
 233     else
 234     {
 235         for(b=0; b<ncons; b++)
 236         {
 237             rvec_sub(x[bla[2*b]],x[bla[2*b+1]],dx);
 238             unitv(dx,r[b]);
 239         } /* 16 ncons flops */
 240     }
 241
 242     for(b=0; b<ncons; b++)
 243     {
 244         tmp0 = r[b][0];
 245         tmp1 = r[b][1];
 246         tmp2 = r[b][2];
 247         i = bla[2*b];
 248         j = bla[2*b+1];
 249         for(n=blnr[b]; n<blnr[b+1]; n++)
 250         {
 251             k = blbnb[n];
 252             blcc[n] = blmf[n]*(tmp0*r[k][0] + tmp1*r[k][1] + tmp2*r[k][2]);
 253         } /* 6 nr flops */
 254         mvb = blc[b]*(tmp0*(f[i][0] - f[j][0]) +
 255                       tmp1*(f[i][1] - f[j][1]) +
 256                       tmp2*(f[i][2] - f[j][2]));
 257         rhs1[b] = mvb;
 258         sol[b]  = mvb;
 259         /* 7 flops */
 260     }
 261     /* Together: 23*ncons + 6*nrtot flops */
 262
 263     lincs_matrix_expand(lincsd,blcc,rhs1,rhs2,sol);
 264     /* nrec*(ncons+2*nrtot) flops */
 265
 266     if (econq != econqForce)
 267     {
 268         for(b=0; b<ncons; b++)
 269         {
 270             /* With econqDeriv_FlexCon only use the flexible constraints */
 271             if (econq != econqDeriv_FlexCon ||
 272                 (lincsd->bllen0[b] == 0 && lincsd->ddist[b] == 0))
 273             {
 274                 i = bla[2*b];
 275                 j = bla[2*b+1];
 276                 mvb = blc[b]*sol[b];
 277                 im1 = invmass[i];
 278                 im2 = invmass[j];
 279                 tmp0 = r[b][0]*mvb;
 280                 tmp1 = r[b][1]*mvb;
 281                 tmp2 = r[b][2]*mvb;
 282                 fp[i][0] -= tmp0*im1;
 283                 fp[i][1] -= tmp1*im1;
 284                 fp[i][2] -= tmp2*im1;
 285                 fp[j][0] += tmp0*im2;
 286                 fp[j][1] += tmp1*im2;
 287                 fp[j][2] += tmp2*im2;
 288                 if (dvdlambda)
 289                 {
 290                     /* This is only correct with forces and invmass=1 */
 291                     *dvdlambda -= mvb*lincsd->ddist[b];
 292                 }
 293             }
 294         } /* 16 ncons flops */
 295     }
 296     else
 297     {
 298         for(b=0; b<ncons; b++)
 299         {
 300             i = bla[2*b];
 301             j = bla[2*b+1];
 302             mvb = blc[b]*sol[b];
 303             tmp0 = r[b][0]*mvb;
 304             tmp1 = r[b][1]*mvb;
 305             tmp2 = r[b][2]*mvb;
 306             fp[i][0] -= tmp0;
 307             fp[i][1] -= tmp1;
 308             fp[i][2] -= tmp2;
 309             fp[j][0] += tmp0;
 310             fp[j][1] += tmp1;
 311             fp[j][2] += tmp2;
 312             if (dvdlambda)
 313             {
 314                 *dvdlambda -= mvb*lincsd->ddist[b];
 315             }
 316         }
 317         /* 10 ncons flops */
 318     }
 319
 320     if (bCalcVir)
 321     {
 322         /* Constraint virial,
 323          * determines sum r_bond x delta f,
 324          * where delta f is the constraint correction
 325          * of the quantity that is being constrained.
 326          */
 327         for(b=0; b<ncons; b++)
 328         {
 329             mvb = lincsd->bllen[b]*blc[b]*sol[b];
 330             for(i=0; i<DIM; i++)
 331             {
 332                 tmp1 = mvb*r[b][i];
 333                 for(j=0; j<DIM; j++)
 334                 {
 335                     rmdf[i][j] += tmp1*r[b][j];
 336                 }
 337             }
 338         } /* 23 ncons flops */
 339     }
 340 }
 341
 342 static void do_lincs(rvec *x,rvec *xp,matrix box,t_pbc *pbc,
 343                      struct gmx_lincsdata *lincsd,real *invmass,
 344                                          t_commrec *cr,
 345                      real wangle,int *warn,
 346                      real invdt,rvec *v,
 347                      gmx_bool bCalcVir,tensor rmdr)
 348 {
 349     int     b,i,j,k,n,iter;
 350     real    tmp0,tmp1,tmp2,im1,im2,mvb,rlen,len,len2,dlen2,wfac,lam;
 351     rvec    dx;
 352     int     ncons,*bla,*blnr,*blbnb;
 353     rvec    *r;
 354     real    *blc,*blmf,*bllen,*blcc,*rhs1,*rhs2,*sol,*lambda;
 355     int     *nlocat;
 356
 357     ncons  = lincsd->nc;
 358     bla    = lincsd->bla;
 359     r      = lincsd->tmpv;
 360     blnr   = lincsd->blnr;
 361     blbnb  = lincsd->blbnb;
 362     blc    = lincsd->blc;
 363     blmf   = lincsd->blmf;
 364     bllen  = lincsd->bllen;
 365     blcc   = lincsd->tmpncc;
 366     rhs1   = lincsd->tmp1;
 367     rhs2   = lincsd->tmp2;
 368     sol    = lincsd->tmp3;
 369     lambda = lincsd->lambda;
 370
 371     if (DOMAINDECOMP(cr) && cr->dd->constraints)
 372     {
 373         nlocat = dd_constraints_nlocalatoms(cr->dd);
 374     }
 375     else if (PARTDECOMP(cr))
 376     {
 377         nlocat = pd_constraints_nlocalatoms(cr->pd);
 378     }
 379     else
 380     {
 381         nlocat = NULL;
 382     }
 383
 384     *warn = 0;
 385
 386     if (pbc)
 387     {
 388         /* Compute normalized i-j vectors */
 389         for(b=0; b<ncons; b++)
 390         {
 391             pbc_dx_aiuc(pbc,x[bla[2*b]],x[bla[2*b+1]],dx);
 392             unitv(dx,r[b]);
 393         }
 394         for(b=0; b<ncons; b++)
 395         {
 396             for(n=blnr[b]; n<blnr[b+1]; n++)
 397             {
 398                 blcc[n] = blmf[n]*iprod(r[b],r[blbnb[n]]);
 399             }
 400             pbc_dx_aiuc(pbc,xp[bla[2*b]],xp[bla[2*b+1]],dx);
 401             mvb = blc[b]*(iprod(r[b],dx) - bllen[b]);
 402             rhs1[b] = mvb;
 403             sol[b]  = mvb;
 404         }
 405     }
 406     else
 407     {
 408         /* Compute normalized i-j vectors */
 409         for(b=0; b<ncons; b++)
 410         {
 411             i = bla[2*b];
 412             j = bla[2*b+1];
 413             tmp0 = x[i][0] - x[j][0];
 414             tmp1 = x[i][1] - x[j][1];
 415             tmp2 = x[i][2] - x[j][2];
 416             rlen = gmx_invsqrt(tmp0*tmp0+tmp1*tmp1+tmp2*tmp2);
 417             r[b][0] = rlen*tmp0;
 418             r[b][1] = rlen*tmp1;
 419             r[b][2] = rlen*tmp2;
 420         } /* 16 ncons flops */
 421
 422         for(b=0; b<ncons; b++)
 423         {
 424             tmp0 = r[b][0];
 425             tmp1 = r[b][1];
 426             tmp2 = r[b][2];
 427             len = bllen[b];
 428             i = bla[2*b];
 429             j = bla[2*b+1];
 430             for(n=blnr[b]; n<blnr[b+1]; n++)
 431             {
 432                 k = blbnb[n];
 433                 blcc[n] = blmf[n]*(tmp0*r[k][0] + tmp1*r[k][1] + tmp2*r[k][2]);
 434             } /* 6 nr flops */
 435             mvb = blc[b]*(tmp0*(xp[i][0] - xp[j][0]) +
 436                           tmp1*(xp[i][1] - xp[j][1]) +
 437                           tmp2*(xp[i][2] - xp[j][2]) - len);
 438             rhs1[b] = mvb;
 439             sol[b]  = mvb;
 440             /* 10 flops */
 441         }
 442         /* Together: 26*ncons + 6*nrtot flops */
 443     }
 444
 445     lincs_matrix_expand(lincsd,blcc,rhs1,rhs2,sol);
 446     /* nrec*(ncons+2*nrtot) flops */
 447
 448     for(b=0; b<ncons; b++)
 449     {
 450         i = bla[2*b];
 451         j = bla[2*b+1];
 452         mvb = blc[b]*sol[b];
 453         lambda[b] = -mvb;
 454         im1 = invmass[i];
 455         im2 = invmass[j];
 456         tmp0 = r[b][0]*mvb;
 457         tmp1 = r[b][1]*mvb;
 458         tmp2 = r[b][2]*mvb;
 459         xp[i][0] -= tmp0*im1;
 460         xp[i][1] -= tmp1*im1;
 461         xp[i][2] -= tmp2*im1;
 462         xp[j][0] += tmp0*im2;
 463         xp[j][1] += tmp1*im2;
 464         xp[j][2] += tmp2*im2;
 465     } /* 16 ncons flops */
 466
 467
 468     /*
 469      ********  Correction for centripetal effects  ********
 470      */
 471
 472     wfac = cos(DEG2RAD*wangle);
 473     wfac = wfac*wfac;
 474
 475     for(iter=0; iter<lincsd->nIter; iter++)
 476     {
 477         if (DOMAINDECOMP(cr) && cr->dd->constraints)
 478         {
 479             /* Communicate the corrected non-local coordinates */
 480             dd_move_x_constraints(cr->dd,box,xp,NULL);
 481         }
 482                 else if (PARTDECOMP(cr))
 483                 {
 484                         pd_move_x_constraints(cr,xp,NULL);
 485                 }
 486
 487         for(b=0; b<ncons; b++)
 488         {
 489             len = bllen[b];
 490             if (pbc)
 491             {
 492                 pbc_dx_aiuc(pbc,xp[bla[2*b]],xp[bla[2*b+1]],dx);
 493             }
 494             else
 495             {
 496                 rvec_sub(xp[bla[2*b]],xp[bla[2*b+1]],dx);
 497             }
 498             len2 = len*len;
 499             dlen2 = 2*len2 - norm2(dx);
 500             if (dlen2 < wfac*len2 && (nlocat==NULL || nlocat[b]))
 501             {
 502                 *warn = b;
 503             }
 504             if (dlen2 > 0)
 505             {
 506                 mvb = blc[b]*(len - dlen2*gmx_invsqrt(dlen2));
 507             }
 508             else
 509             {
 510                 mvb = blc[b]*len;
 511             }
 512             rhs1[b] = mvb;
 513             sol[b]  = mvb;
 514         } /* 20*ncons flops */
 515
 516         lincs_matrix_expand(lincsd,blcc,rhs1,rhs2,sol);
 517         /* nrec*(ncons+2*nrtot) flops */
 518
 519         for(b=0; b<ncons; b++)
 520         {
 521             i = bla[2*b];
 522             j = bla[2*b+1];
 523             lam = lambda[b];
 524             mvb = blc[b]*sol[b];
 525             lambda[b] = lam - mvb;
 526             im1 = invmass[i];
 527             im2 = invmass[j];
 528             tmp0 = r[b][0]*mvb;
 529             tmp1 = r[b][1]*mvb;
 530             tmp2 = r[b][2]*mvb;
 531             xp[i][0] -= tmp0*im1;
 532             xp[i][1] -= tmp1*im1;
 533             xp[i][2] -= tmp2*im1;
 534             xp[j][0] += tmp0*im2;
 535             xp[j][1] += tmp1*im2;
 536             xp[j][2] += tmp2*im2;
 537         } /* 17 ncons flops */
 538     } /* nit*ncons*(37+9*nrec) flops */
 539
 540     if (v)
 541     {
 542         /* Correct the velocities */
 543         for(b=0; b<ncons; b++)
 544         {
 545             i = bla[2*b];
 546             j = bla[2*b+1];
 547             im1 = invmass[i]*lambda[b]*invdt;
 548             im2 = invmass[j]*lambda[b]*invdt;
 549             v[i][0] += im1*r[b][0];
 550             v[i][1] += im1*r[b][1];
 551             v[i][2] += im1*r[b][2];
 552             v[j][0] -= im2*r[b][0];
 553             v[j][1] -= im2*r[b][1];
 554             v[j][2] -= im2*r[b][2];
 555         } /* 16 ncons flops */
 556     }
 557
 558     if (nlocat)
 559     {
 560         /* Only account for local atoms */
 561         for(b=0; b<ncons; b++)
 562         {
 563             lambda[b] *= 0.5*nlocat[b];
 564         }
 565     }
 566
 567     if (bCalcVir)
 568     {
 569         /* Constraint virial */
 570         for(b=0; b<ncons; b++)
 571         {
 572             tmp0 = bllen[b]*lambda[b];
 573             for(i=0; i<DIM; i++)
 574             {
 575                 tmp1 = tmp0*r[b][i];
 576                 for(j=0; j<DIM; j++)
 577                 {
 578                     rmdr[i][j] -= tmp1*r[b][j];
 579                 }
 580             }
 581         } /* 22 ncons flops */
 582     }
 583
 584     /* Total:
 585      * 26*ncons + 6*nrtot + nrec*(ncons+2*nrtot)
 586      * + nit * (20*ncons + nrec*(ncons+2*nrtot) + 17 ncons)
 587      *
 588      * (26+nrec)*ncons + (6+2*nrec)*nrtot
 589      * + nit * ((37+nrec)*ncons + 2*nrec*nrtot)
 590      * if nit=1
 591      * (63+nrec)*ncons + (6+4*nrec)*nrtot
 592      */
 593 }
 594
 595 void set_lincs_matrix(struct gmx_lincsdata *li,real *invmass,real lambda)
 596 {
 597     int i,a1,a2,n,k,sign,center;
 598     int end,nk,kk;
 599     const real invsqrt2=0.7071067811865475244;
 600
 601     for(i=0; (i<li->nc); i++)
 602     {
 603         a1 = li->bla[2*i];
 604         a2 = li->bla[2*i+1];
 605         li->blc[i]  = gmx_invsqrt(invmass[a1] + invmass[a2]);
 606         li->blc1[i] = invsqrt2;
 607     }
 608
 609     /* Construct the coupling coefficient matrix blmf */
 610     li->ntriangle = 0;
 611     li->ncc_triangle = 0;
 612     for(i=0; (i<li->nc); i++)
 613     {
 614         a1 = li->bla[2*i];
 615         a2 = li->bla[2*i+1];
 616         for(n=li->blnr[i]; (n<li->blnr[i+1]); n++)
 617         {
 618             k = li->blbnb[n];
 619             if (a1 == li->bla[2*k] || a2 == li->bla[2*k+1])
 620             {
 621                 sign = -1;
 622             }
 623             else
 624             {
 625                 sign = 1;
 626             }
 627             if (a1 == li->bla[2*k] || a1 == li->bla[2*k+1])
 628             {
 629                 center = a1;
 630                 end    = a2;
 631             }
 632             else
 633             {
 634                 center = a2;
 635                 end    = a1;
 636             }
 637             li->blmf[n]  = sign*invmass[center]*li->blc[i]*li->blc[k];
 638             li->blmf1[n] = sign*0.5;
 639             if (li->ncg_triangle > 0)
 640             {
 641                 /* Look for constraint triangles */
 642                 for(nk=li->blnr[k]; (nk<li->blnr[k+1]); nk++)
 643                 {
 644                     kk = li->blbnb[nk];
 645                     if (kk != i && kk != k &&
 646                         (li->bla[2*kk] == end || li->bla[2*kk+1] == end))
 647                     {
 648                         if (li->ntriangle == 0 ||
 649                             li->triangle[li->ntriangle-1] < i)
 650                         {
 651                             /* Add this constraint to the triangle list */
 652                             li->triangle[li->ntriangle] = i;
 653                             li->tri_bits[li->ntriangle] = 0;
 654                             li->ntriangle++;
 655                             if (li->blnr[i+1] - li->blnr[i] > sizeof(li->tri_bits[0])*8 - 1)
 656                             {
 657                                 gmx_fatal(FARGS,"A constraint is connected to %d constraints, this is more than the %d allowed for constraints participating in triangles",
 658                                           li->blnr[i+1] - li->blnr[i],
 659                                           sizeof(li->tri_bits[0])*8-1);
 660                             }
 661                         }
 662                         li->tri_bits[li->ntriangle-1] |= (1<<(n-li->blnr[i]));
 663                         li->ncc_triangle++;
 664                     }
 665                 }
 666             }
 667         }
 668     }
 669
 670     if (debug)
 671     {
 672         fprintf(debug,"Of the %d constraints %d participate in triangles\n",
 673                 li->nc,li->ntriangle);
 674         fprintf(debug,"There are %d couplings of which %d in triangles\n",
 675                 li->ncc,li->ncc_triangle);
 676     }
 677
 678     /* Set matlam,
 679      * so we know with which lambda value the masses have been set.
 680      */
 681     li->matlam = lambda;
 682 }
 683
 684 static int count_triangle_constraints(t_ilist *ilist,t_blocka *at2con)
 685 {
 686     int  ncon1,ncon_tot;
 687     int  c0,a00,a01,n1,c1,a10,a11,ac1,n2,c2,a20,a21;
 688     int  ncon_triangle;
 689     gmx_bool bTriangle;
 690     t_iatom *ia1,*ia2,*iap;
 691
 692     ncon1    = ilist[F_CONSTR].nr/3;
 693     ncon_tot = ncon1 + ilist[F_CONSTRNC].nr/3;
 694
 695     ia1 = ilist[F_CONSTR].iatoms;
 696     ia2 = ilist[F_CONSTRNC].iatoms;
 697
 698     ncon_triangle = 0;
 699     for(c0=0; c0<ncon_tot; c0++)
 700     {
 701         bTriangle = FALSE;
 702         iap = constr_iatomptr(ncon1,ia1,ia2,c0);
 703         a00 = iap[1];
 704         a01 = iap[2];
 705         for(n1=at2con->index[a01]; n1<at2con->index[a01+1]; n1++)
 706         {
 707             c1 = at2con->a[n1];
 708             if (c1 != c0)
 709             {
 710                 iap = constr_iatomptr(ncon1,ia1,ia2,c1);
 711                 a10 = iap[1];
 712                 a11 = iap[2];
 713                 if (a10 == a01)
 714                 {
 715                     ac1 = a11;
 716                 }
 717                 else
 718                 {
 719                     ac1 = a10;
 720                 }
 721                 for(n2=at2con->index[ac1]; n2<at2con->index[ac1+1]; n2++)
 722                 {
 723                     c2 = at2con->a[n2];
 724                     if (c2 != c0 && c2 != c1)
 725                     {
 726                         iap = constr_iatomptr(ncon1,ia1,ia2,c2);
 727                         a20 = iap[1];
 728                         a21 = iap[2];
 729                         if (a20 == a00 || a21 == a00)
 730                         {
 731                             bTriangle = TRUE;
 732                         }
 733                     }
 734                 }
 735             }
 736         }
 737         if (bTriangle)
 738         {
 739             ncon_triangle++;
 740         }
 741     }
 742
 743     return ncon_triangle;
 744 }
 745
 746 static int int_comp(const void *a,const void *b)
 747 {
 748     return (*(int *)a) - (*(int *)b);
 749 }
 750
 751 gmx_lincsdata_t init_lincs(FILE *fplog,gmx_mtop_t *mtop,
 752                            int nflexcon_global,t_blocka *at2con,
 753                            gmx_bool bPLINCS,int nIter,int nProjOrder)
 754 {
 755     struct gmx_lincsdata *li;
 756     int mb;
 757     gmx_moltype_t *molt;
 758
 759     if (fplog)
 760     {
 761         fprintf(fplog,"\nInitializing%s LINear Constraint Solver\n",
 762                 bPLINCS ? " Parallel" : "");
 763     }
 764
 765     snew(li,1);
 766
 767     li->ncg      =
 768         gmx_mtop_ftype_count(mtop,F_CONSTR) +
 769         gmx_mtop_ftype_count(mtop,F_CONSTRNC);
 770     li->ncg_flex = nflexcon_global;
 771
 772     li->ncg_triangle = 0;
 773     for(mb=0; mb<mtop->nmolblock; mb++)
 774     {
 775         molt = &mtop->moltype[mtop->molblock[mb].type];
 776         li->ncg_triangle +=
 777             mtop->molblock[mb].nmol*
 778             count_triangle_constraints(molt->ilist,
 779                                        &at2con[mtop->molblock[mb].type]);
 780     }
 781
 782     li->nIter  = nIter;
 783     li->nOrder = nProjOrder;
 784
 785     if (bPLINCS || li->ncg_triangle > 0)
 786     {
 787         please_cite(fplog,"Hess2008a");
 788     }
 789     else
 790     {
 791         please_cite(fplog,"Hess97a");
 792     }
 793
 794     if (fplog)
 795     {
 796         fprintf(fplog,"The number of constraints is %d\n",li->ncg);
 797         if (bPLINCS)
 798         {
 799             fprintf(fplog,"There are inter charge-group constraints,\n"
 800                     "will communicate selected coordinates each lincs iteration\n");
 801         }
 802         if (li->ncg_triangle > 0)
 803         {
 804             fprintf(fplog,
 805                     "%d constraints are involved in constraint triangles,\n"
 806                     "will apply an additional matrix expansion of order %d for couplings\n"
 807                     "between constraints inside triangles\n",
 808                     li->ncg_triangle,li->nOrder);
 809         }
 810     }
 811
 812     return li;
 813 }
 814
 815 void set_lincs(t_idef *idef,t_mdatoms *md,
 816                gmx_bool bDynamics,t_commrec *cr,
 817                struct gmx_lincsdata *li)
 818 {
 819     int      start,natoms,nflexcon;
 820     t_blocka at2con;
 821     t_iatom  *iatom;
 822     int      i,k,ncc_alloc,ni,con,nconnect,concon;
 823     int      type,a1,a2;
 824     real     lenA=0,lenB;
 825     gmx_bool     bLocal;
 826
 827     li->nc = 0;
 828     li->ncc = 0;
 829
 830     /* This is the local topology, so there are only F_CONSTR constraints */
 831     if (idef->il[F_CONSTR].nr == 0)
 832     {
 833         /* There are no constraints,
 834          * we do not need to fill any data structures.
 835          */
 836         return;
 837     }
 838
 839     if (debug)
 840     {
 841         fprintf(debug,"Building the LINCS connectivity\n");
 842     }
 843
 844     if (DOMAINDECOMP(cr))
 845     {
 846         if (cr->dd->constraints)
 847         {
 848             dd_get_constraint_range(cr->dd,&start,&natoms);
 849         }
 850         else
 851         {
 852             natoms = cr->dd->nat_home;
 853         }
 854         start = 0;
 855     }
 856     else if(PARTDECOMP(cr))
 857         {
 858                 pd_get_constraint_range(cr->pd,&start,&natoms);
 859         }
 860         else
 861     {
 862         start  = md->start;
 863         natoms = md->homenr;
 864     }
 865     at2con = make_at2con(start,natoms,idef->il,idef->iparams,bDynamics,
 866                          &nflexcon);
 867
 868
 869     if (idef->il[F_CONSTR].nr/3 > li->nc_alloc || li->nc_alloc == 0)
 870     {
 871         li->nc_alloc = over_alloc_dd(idef->il[F_CONSTR].nr/3);
 872         srenew(li->bllen0,li->nc_alloc);
 873         srenew(li->ddist,li->nc_alloc);
 874         srenew(li->bla,2*li->nc_alloc);
 875         srenew(li->blc,li->nc_alloc);
 876         srenew(li->blc1,li->nc_alloc);
 877         srenew(li->blnr,li->nc_alloc+1);
 878         srenew(li->bllen,li->nc_alloc);
 879         srenew(li->tmpv,li->nc_alloc);
 880         srenew(li->tmp1,li->nc_alloc);
 881         srenew(li->tmp2,li->nc_alloc);
 882         srenew(li->tmp3,li->nc_alloc);
 883         srenew(li->lambda,li->nc_alloc);
 884         if (li->ncg_triangle > 0)
 885         {
 886             /* This is allocating too much, but it is difficult to improve */
 887             srenew(li->triangle,li->nc_alloc);
 888             srenew(li->tri_bits,li->nc_alloc);
 889         }
 890     }
 891
 892     iatom = idef->il[F_CONSTR].iatoms;
 893
 894     ncc_alloc = li->ncc_alloc;
 895     li->blnr[0] = 0;
 896
 897     ni = idef->il[F_CONSTR].nr/3;
 898
 899     con = 0;
 900     nconnect = 0;
 901     li->blnr[con] = nconnect;
 902     for(i=0; i<ni; i++)
 903     {
 904         bLocal = TRUE;
 905         type = iatom[3*i];
 906         a1   = iatom[3*i+1];
 907         a2   = iatom[3*i+2];
 908         lenA = idef->iparams[type].constr.dA;
 909         lenB = idef->iparams[type].constr.dB;
 910         /* Skip the flexible constraints when not doing dynamics */
 911         if (bDynamics || lenA!=0 || lenB!=0)
 912         {
 913             li->bllen0[con]  = lenA;
 914             li->ddist[con]   = lenB - lenA;
 915             /* Set the length to the topology A length */
 916             li->bllen[con]   = li->bllen0[con];
 917             li->bla[2*con]   = a1;
 918             li->bla[2*con+1] = a2;
 919             /* Construct the constraint connection matrix blbnb */
 920             for(k=at2con.index[a1-start]; k<at2con.index[a1-start+1]; k++)
 921             {
 922                 concon = at2con.a[k];
 923                 if (concon != i)
 924                 {
 925                     if (nconnect >= ncc_alloc)
 926                     {
 927                         ncc_alloc = over_alloc_small(nconnect+1);
 928                         srenew(li->blbnb,ncc_alloc);
 929                     }
 930                     li->blbnb[nconnect++] = concon;
 931                 }
 932             }
 933             for(k=at2con.index[a2-start]; k<at2con.index[a2-start+1]; k++)
 934             {
 935                 concon = at2con.a[k];
 936                 if (concon != i)
 937                 {
 938                     if (nconnect+1 > ncc_alloc)
 939                     {
 940                         ncc_alloc = over_alloc_small(nconnect+1);
 941                         srenew(li->blbnb,ncc_alloc);
 942                     }
 943                     li->blbnb[nconnect++] = concon;
 944                 }
 945             }
 946             li->blnr[con+1] = nconnect;
 947
 948             if (cr->dd == NULL)
 949             {
 950                 /* Order the blbnb matrix to optimize memory access */
 951                 qsort(&(li->blbnb[li->blnr[con]]),li->blnr[con+1]-li->blnr[con],
 952                       sizeof(li->blbnb[0]),int_comp);
 953             }
 954             /* Increase the constraint count */
 955             con++;
 956         }
 957     }
 958
 959     done_blocka(&at2con);
 960
 961     /* This is the real number of constraints,
 962      * without dynamics the flexible constraints are not present.
 963      */
 964     li->nc = con;
 965
 966     li->ncc = li->blnr[con];
 967     if (cr->dd == NULL)
 968     {
 969         /* Since the matrix is static, we can free some memory */
 970         ncc_alloc = li->ncc;
 971         srenew(li->blbnb,ncc_alloc);
 972     }
 973
 974     if (ncc_alloc > li->ncc_alloc)
 975     {
 976         li->ncc_alloc = ncc_alloc;
 977         srenew(li->blmf,li->ncc_alloc);
 978         srenew(li->blmf1,li->ncc_alloc);
 979         srenew(li->tmpncc,li->ncc_alloc);
 980     }
 981
 982     if (debug)
 983     {
 984         fprintf(debug,"Number of constraints is %d, couplings %d\n",
 985                 li->nc,li->ncc);
 986     }
 987
 988     set_lincs_matrix(li,md->invmass,md->lambda);
 989 }
 990
 991 static void lincs_warning(FILE *fplog,
 992                           gmx_domdec_t *dd,rvec *x,rvec *xprime,t_pbc *pbc,
 993                           int ncons,int *bla,real *bllen,real wangle,
 994                           int maxwarn,int *warncount)
 995 {
 996     int b,i,j;
 997     rvec v0,v1;
 998     real wfac,d0,d1,cosine;
 999     char buf[STRLEN];
1000
1001     wfac=cos(DEG2RAD*wangle);
1002
1003     sprintf(buf,"bonds that rotated more than %g degrees:\n"
1004             " atom 1 atom 2  angle  previous, current, constraint length\n",
1005             wangle);
1006     fprintf(stderr,"%s",buf);
1007     if (fplog)
1008     {
1009         fprintf(fplog,"%s",buf);
1010     }
1011
1012     for(b=0;b<ncons;b++)
1013     {
1014         i = bla[2*b];
1015         j = bla[2*b+1];
1016         if (pbc)
1017         {
1018             pbc_dx_aiuc(pbc,x[i],x[j],v0);
1019             pbc_dx_aiuc(pbc,xprime[i],xprime[j],v1);
1020         }
1021         else
1022         {
1023             rvec_sub(x[i],x[j],v0);
1024             rvec_sub(xprime[i],xprime[j],v1);
1025         }
1026         d0 = norm(v0);
1027         d1 = norm(v1);
1028         cosine = iprod(v0,v1)/(d0*d1);
1029         if (cosine < wfac)
1030         {
1031             sprintf(buf," %6d %6d  %5.1f  %8.4f %8.4f    %8.4f\n",
1032                     ddglatnr(dd,i),ddglatnr(dd,j),
1033                     RAD2DEG*acos(cosine),d0,d1,bllen[b]);
1034             fprintf(stderr,"%s",buf);
1035             if (fplog)
1036             {
1037                 fprintf(fplog,"%s",buf);
1038             }
1039             if (!gmx_isfinite(d1))
1040             {
1041                 gmx_fatal(FARGS,"Bond length not finite.");
1042             }
1043
1044             (*warncount)++;
1045         }
1046     }
1047     if (*warncount > maxwarn)
1048     {
1049         too_many_constraint_warnings(econtLINCS,*warncount);
1050     }
1051 }
1052
1053 static void cconerr(gmx_domdec_t *dd,
1054                     int ncons,int *bla,real *bllen,rvec *x,t_pbc *pbc,
1055                     real *ncons_loc,real *ssd,real *max,int *imax)
1056 {
1057     real      len,d,ma,ssd2,r2;
1058     int       *nlocat,count,b,im;
1059     rvec      dx;
1060
1061     if (dd && dd->constraints)
1062     {
1063         nlocat = dd_constraints_nlocalatoms(dd);
1064     }
1065     else
1066     {
1067         nlocat = 0;
1068     }
1069
1070     ma = 0;
1071     ssd2 = 0;
1072     im = 0;
1073     count = 0;
1074     for(b=0;b<ncons;b++)
1075     {
1076         if (pbc)
1077         {
1078             pbc_dx_aiuc(pbc,x[bla[2*b]],x[bla[2*b+1]],dx);
1079         }
1080         else {
1081             rvec_sub(x[bla[2*b]],x[bla[2*b+1]],dx);
1082         }
1083         r2 = norm2(dx);
1084         len = r2*gmx_invsqrt(r2);
1085         d = fabs(len/bllen[b]-1);
1086         if (d > ma && (nlocat==NULL || nlocat[b]))
1087         {
1088             ma = d;
1089             im = b;
1090         }
1091         if (nlocat == NULL)
1092         {
1093             ssd2 += d*d;
1094             count++;
1095         }
1096         else
1097         {
1098             ssd2 += nlocat[b]*d*d;
1099             count += nlocat[b];
1100         }
1101     }
1102
1103     *ncons_loc = (nlocat ? 0.5 : 1)*count;
1104     *ssd       = (nlocat ? 0.5 : 1)*ssd2;
1105     *max = ma;
1106     *imax = im;
1107 }
1108
1109 static void dump_conf(gmx_domdec_t *dd,struct gmx_lincsdata *li,
1110                       t_blocka *at2con,
1111                       char *name,gmx_bool bAll,rvec *x,matrix box)
1112 {
1113     char str[STRLEN];
1114     FILE *fp;
1115     int  ac0,ac1,i;
1116
1117     dd_get_constraint_range(dd,&ac0,&ac1);
1118
1119     sprintf(str,"%s_%d_%d_%d.pdb",name,dd->ci[XX],dd->ci[YY],dd->ci[ZZ]);
1120     fp = gmx_fio_fopen(str,"w");
1121     fprintf(fp,"CRYST1%9.3f%9.3f%9.3f%7.2f%7.2f%7.2f P 1           1\n",
1122             10*norm(box[XX]),10*norm(box[YY]),10*norm(box[ZZ]),
1123             90.0,90.0,90.0);
1124     for(i=0; i<ac1; i++)
1125     {
1126         if (i < dd->nat_home || (bAll && i >= ac0 && i < ac1))
1127         {
1128             fprintf(fp,"%-6s%5u  %-4.4s%3.3s %c%4d    %8.3f%8.3f%8.3f%6.2f%6.2f\n",
1129                     "ATOM",ddglatnr(dd,i),"C","ALA",' ',i+1,
1130                     10*x[i][XX],10*x[i][YY],10*x[i][ZZ],
1131                     1.0,i<dd->nat_tot ? 0.0 : 1.0);
1132         }
1133     }
1134     if (bAll)
1135     {
1136         for(i=0; i<li->nc; i++)
1137         {
1138             fprintf(fp,"CONECT%5d%5d\n",
1139                     ddglatnr(dd,li->bla[2*i]),
1140                     ddglatnr(dd,li->bla[2*i+1]));
1141         }
1142     }
1143     gmx_fio_fclose(fp);
1144 }
1145
1146 gmx_bool constrain_lincs(FILE *fplog,gmx_bool bLog,gmx_bool bEner,
1147                      t_inputrec *ir,
1148                      gmx_large_int_t step,
1149                      struct gmx_lincsdata *lincsd,t_mdatoms *md,
1150                      t_commrec *cr,
1151                      rvec *x,rvec *xprime,rvec *min_proj,matrix box,
1152                      real lambda,real *dvdlambda,
1153                      real invdt,rvec *v,
1154                      gmx_bool bCalcVir,tensor rmdr,
1155                      int econq,
1156                      t_nrnb *nrnb,
1157                      int maxwarn,int *warncount)
1158 {
1159     char  buf[STRLEN],buf2[22],buf3[STRLEN];
1160     int   i,warn,p_imax,error;
1161     real  ncons_loc,p_ssd,p_max=0;
1162     t_pbc pbc,*pbc_null;
1163     rvec  dx;
1164     gmx_bool  bOK;
1165
1166     bOK = TRUE;
1167
1168     if (lincsd->nc == 0 && cr->dd == NULL)
1169     {
1170         if (bLog || bEner)
1171         {
1172             lincsd->rmsd_data[0] = 0;
1173             if (ir->eI == eiSD2 && v == NULL)
1174             {
1175                 i = 2;
1176             }
1177             else
1178             {
1179                 i = 1;
1180             }
1181             lincsd->rmsd_data[i] = 0;
1182         }
1183
1184         return bOK;
1185     }
1186
1187     /* We do not need full pbc when constraints do not cross charge groups,
1188      * i.e. when dd->constraint_comm==NULL
1189      */
1190     if ((cr->dd || ir->bPeriodicMols) && !(cr->dd && cr->dd->constraint_comm==NULL))
1191     {
1192         /* With pbc=screw the screw has been changed to a shift
1193          * by the constraint coordinate communication routine,
1194          * so that here we can use normal pbc.
1195          */
1196         pbc_null = set_pbc_dd(&pbc,ir->ePBC,cr->dd,FALSE,box);
1197     }
1198     else
1199     {
1200         pbc_null = NULL;
1201     }
1202     if (cr->dd)
1203     {
1204         /* Communicate the coordinates required for the non-local constraints */
1205         dd_move_x_constraints(cr->dd,box,x,xprime);
1206         /* dump_conf(dd,lincsd,NULL,"con",TRUE,xprime,box); */
1207     }
1208         else if (PARTDECOMP(cr))
1209         {
1210                 pd_move_x_constraints(cr,x,xprime);
1211         }
1212
1213     if (econq == econqCoord)
1214     {
1215         if (ir->efep != efepNO)
1216         {
1217             if (md->nMassPerturbed && lincsd->matlam != md->lambda)
1218             {
1219                 set_lincs_matrix(lincsd,md->invmass,md->lambda);
1220             }
1221
1222             for(i=0; i<lincsd->nc; i++)
1223             {
1224                 lincsd->bllen[i] = lincsd->bllen0[i] + lambda*lincsd->ddist[i];
1225             }
1226         }
1227
1228         if (lincsd->ncg_flex)
1229         {
1230             /* Set the flexible constraint lengths to the old lengths */
1231             if (pbc_null)
1232             {
1233                 for(i=0; i<lincsd->nc; i++)
1234                 {
1235                     if (lincsd->bllen[i] == 0) {
1236                         pbc_dx_aiuc(pbc_null,x[lincsd->bla[2*i]],x[lincsd->bla[2*i+1]],dx);
1237                         lincsd->bllen[i] = norm(dx);
1238                     }
1239                 }
1240             }
1241             else
1242             {
1243                 for(i=0; i<lincsd->nc; i++)
1244                 {
1245                     if (lincsd->bllen[i] == 0)
1246                     {
1247                         lincsd->bllen[i] =
1248                             sqrt(distance2(x[lincsd->bla[2*i]],
1249                                            x[lincsd->bla[2*i+1]]));
1250                     }
1251                 }
1252             }
1253         }
1254
1255         if (bLog && fplog)
1256         {
1257             cconerr(cr->dd,lincsd->nc,lincsd->bla,lincsd->bllen,xprime,pbc_null,
1258                     &ncons_loc,&p_ssd,&p_max,&p_imax);
1259         }
1260
1261         do_lincs(x,xprime,box,pbc_null,lincsd,md->invmass,cr,
1262                  ir->LincsWarnAngle,&warn,
1263                  invdt,v,bCalcVir,rmdr);
1264
1265         if (ir->efep != efepNO)
1266         {
1267             real dt_2,dvdl=0;
1268
1269             dt_2 = 1.0/(ir->delta_t*ir->delta_t);
1270             for(i=0; (i<lincsd->nc); i++)
1271             {
1272                 dvdl += lincsd->lambda[i]*dt_2*lincsd->ddist[i];
1273             }
1274             *dvdlambda += dvdl;
1275                 }
1276
1277         if (bLog && fplog && lincsd->nc > 0)
1278         {
1279             fprintf(fplog,"   Rel. Constraint Deviation:  RMS         MAX     between atoms\n");
1280             fprintf(fplog,"       Before LINCS          %.6f    %.6f %6d %6d\n",
1281                     sqrt(p_ssd/ncons_loc),p_max,
1282                     ddglatnr(cr->dd,lincsd->bla[2*p_imax]),
1283                     ddglatnr(cr->dd,lincsd->bla[2*p_imax+1]));
1284         }
1285         if (bLog || bEner)
1286         {
1287             cconerr(cr->dd,lincsd->nc,lincsd->bla,lincsd->bllen,xprime,pbc_null,
1288                     &ncons_loc,&p_ssd,&p_max,&p_imax);
1289             /* Check if we are doing the second part of SD */
1290             if (ir->eI == eiSD2 && v == NULL)
1291             {
1292                 i = 2;
1293             }
1294             else
1295             {
1296                 i = 1;
1297             }
1298             lincsd->rmsd_data[0] = ncons_loc;
1299             lincsd->rmsd_data[i] = p_ssd;
1300         }
1301         else
1302         {
1303             lincsd->rmsd_data[0] = 0;
1304             lincsd->rmsd_data[1] = 0;
1305             lincsd->rmsd_data[2] = 0;
1306         }
1307         if (bLog && fplog && lincsd->nc > 0)
1308         {
1309             fprintf(fplog,
1310                     "        After LINCS          %.6f    %.6f %6d %6d\n\n",
1311                     sqrt(p_ssd/ncons_loc),p_max,
1312                     ddglatnr(cr->dd,lincsd->bla[2*p_imax]),
1313                     ddglatnr(cr->dd,lincsd->bla[2*p_imax+1]));
1314         }
1315
1316         if (warn > 0)
1317         {
1318             if (maxwarn >= 0)
1319             {
1320                 cconerr(cr->dd,lincsd->nc,lincsd->bla,lincsd->bllen,xprime,pbc_null,
1321                         &ncons_loc,&p_ssd,&p_max,&p_imax);
1322                 if (MULTISIM(cr))
1323                 {
1324                     sprintf(buf3," in simulation %d", cr->ms->sim);
1325                 }
1326                 else
1327                 {
1328                     buf3[0] = 0;
1329                 }
1330                 sprintf(buf,"\nStep %s, time %g (ps)  LINCS WARNING%s\n"
1331                         "relative constraint deviation after LINCS:\n"
1332                         "rms %.6f, max %.6f (between atoms %d and %d)\n",
1333                         gmx_step_str(step,buf2),ir->init_t+step*ir->delta_t,
1334                         buf3,
1335                         sqrt(p_ssd/ncons_loc),p_max,
1336                         ddglatnr(cr->dd,lincsd->bla[2*p_imax]),
1337                         ddglatnr(cr->dd,lincsd->bla[2*p_imax+1]));
1338                 if (fplog)
1339                 {
1340                     fprintf(fplog,"%s",buf);
1341                 }
1342                 fprintf(stderr,"%s",buf);
1343                 lincs_warning(fplog,cr->dd,x,xprime,pbc_null,
1344                               lincsd->nc,lincsd->bla,lincsd->bllen,
1345                               ir->LincsWarnAngle,maxwarn,warncount);
1346             }
1347             bOK = (p_max < 0.5);
1348         }
1349
1350         if (lincsd->ncg_flex) {
1351             for(i=0; (i<lincsd->nc); i++)
1352                 if (lincsd->bllen0[i] == 0 && lincsd->ddist[i] == 0)
1353                     lincsd->bllen[i] = 0;
1354         }
1355     }
1356     else
1357     {
1358         do_lincsp(x,xprime,min_proj,pbc_null,lincsd,md->invmass,econq,dvdlambda,
1359                   bCalcVir,rmdr);
1360     }
1361
1362     /* count assuming nit=1 */
1363     inc_nrnb(nrnb,eNR_LINCS,lincsd->nc);
1364     inc_nrnb(nrnb,eNR_LINCSMAT,(2+lincsd->nOrder)*lincsd->ncc);
1365     if (lincsd->ntriangle > 0)
1366     {
1367         inc_nrnb(nrnb,eNR_LINCSMAT,lincsd->nOrder*lincsd->ncc_triangle);
1368     }
1369     if (v)
1370     {
1371         inc_nrnb(nrnb,eNR_CONSTR_V,lincsd->nc*2);
1372     }
1373     if (bCalcVir)
1374     {
1375         inc_nrnb(nrnb,eNR_CONSTR_VIR,lincsd->nc);
1376     }
1377
1378     return bOK;
1379 }