src/gromacs/mdlib/shakef.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 #include "gmxpre.h"
  38
  39 #include "config.h"
  40
  41 #include <math.h>
  42
  43 #include "gromacs/legacyheaders/typedefs.h"
  44 #include "gromacs/utility/smalloc.h"
  45 #include "gromacs/legacyheaders/txtdump.h"
  46 #include "gromacs/math/vec.h"
  47 #include "gromacs/legacyheaders/nrnb.h"
  48 #include "gromacs/legacyheaders/constr.h"
  49
  50 typedef struct gmx_shakedata
  51 {
  52     rvec *rij;
  53     real *M2;
  54     real *tt;
  55     real *dist2;
  56     int   nalloc;
  57     /* SOR stuff */
  58     real  delta;
  59     real  omega;
  60     real  gamma;
  61 } t_gmx_shakedata;
  62
  63 gmx_shakedata_t shake_init()
  64 {
  65     gmx_shakedata_t d;
  66
  67     snew(d, 1);
  68
  69     d->nalloc = 0;
  70     d->rij    = NULL;
  71     d->M2     = NULL;
  72     d->tt     = NULL;
  73     d->dist2  = NULL;
  74
  75     /* SOR initialization */
  76     d->delta = 0.1;
  77     d->omega = 1.0;
  78     d->gamma = 1000000;
  79
  80     return d;
  81 }
  82
  83 static void pv(FILE *log, char *s, rvec x)
  84 {
  85     int m;
  86
  87     fprintf(log, "%5s:", s);
  88     for (m = 0; (m < DIM); m++)
  89     {
  90         fprintf(log, "  %10.3f", x[m]);
  91     }
  92     fprintf(log, "\n");
  93     fflush(log);
  94 }
  95
  96 void cshake(atom_id iatom[], int ncon, int *nnit, int maxnit,
  97             real dist2[], real xp[], real rij[], real m2[], real omega,
  98             real invmass[], real tt[], real lagr[], int *nerror)
  99 {
 100     /*
 101      *     r.c. van schaik and w.f. van gunsteren
 102      *     eth zuerich
 103      *     june 1992
 104      *     Adapted for use with Gromacs by David van der Spoel november 92 and later.
 105      */
 106     /* default should be increased! MRS 8/4/2009 */
 107     const   real mytol = 1e-10;
 108
 109     int          ll, i, j, i3, j3, l3;
 110     int          ix, iy, iz, jx, jy, jz;
 111     real         toler, rpij2, rrpr, tx, ty, tz, diff, acor, im, jm;
 112     real         xh, yh, zh, rijx, rijy, rijz;
 113     real         tix, tiy, tiz;
 114     real         tjx, tjy, tjz;
 115     int          nit, error, nconv;
 116     real         iconvf;
 117
 118     error = 0;
 119     nconv = 1;
 120     for (nit = 0; (nit < maxnit) && (nconv != 0) && (error == 0); nit++)
 121     {
 122         nconv = 0;
 123         for (ll = 0; (ll < ncon) && (error == 0); ll++)
 124         {
 125             l3    = 3*ll;
 126             rijx  = rij[l3+XX];
 127             rijy  = rij[l3+YY];
 128             rijz  = rij[l3+ZZ];
 129             i     = iatom[l3+1];
 130             j     = iatom[l3+2];
 131             i3    = 3*i;
 132             j3    = 3*j;
 133             ix    = i3+XX;
 134             iy    = i3+YY;
 135             iz    = i3+ZZ;
 136             jx    = j3+XX;
 137             jy    = j3+YY;
 138             jz    = j3+ZZ;
 139
 140             tx      = xp[ix]-xp[jx];
 141             ty      = xp[iy]-xp[jy];
 142             tz      = xp[iz]-xp[jz];
 143             rpij2   = tx*tx+ty*ty+tz*tz;
 144             toler   = dist2[ll];
 145             diff    = toler-rpij2;
 146
 147             /* iconvf is less than 1 when the error is smaller than a bound */
 148             /* But if tt is too big, then it will result in looping in iconv */
 149
 150             iconvf = fabs(diff)*tt[ll];
 151
 152             if (iconvf > 1)
 153             {
 154                 nconv   = iconvf;
 155                 rrpr    = rijx*tx+rijy*ty+rijz*tz;
 156
 157                 if (rrpr < toler*mytol)
 158                 {
 159                     error = ll+1;
 160                 }
 161                 else
 162                 {
 163                     acor      = omega*diff*m2[ll]/rrpr;
 164                     lagr[ll] += acor;
 165                     xh        = rijx*acor;
 166                     yh        = rijy*acor;
 167                     zh        = rijz*acor;
 168                     im        = invmass[i];
 169                     jm        = invmass[j];
 170                     xp[ix]   += xh*im;
 171                     xp[iy]   += yh*im;
 172                     xp[iz]   += zh*im;
 173                     xp[jx]   -= xh*jm;
 174                     xp[jy]   -= yh*jm;
 175                     xp[jz]   -= zh*jm;
 176                 }
 177             }
 178         }
 179     }
 180     *nnit   = nit;
 181     *nerror = error;
 182 }
 183
 184 int vec_shakef(FILE *fplog, gmx_shakedata_t shaked,
 185                real invmass[], int ncon,
 186                t_iparams ip[], t_iatom *iatom,
 187                real tol, rvec x[], rvec prime[], real omega,
 188                gmx_bool bFEP, real lambda, real lagr[],
 189                real invdt, rvec *v,
 190                gmx_bool bCalcVir, tensor vir_r_m_dr, int econq,
 191                t_vetavars *vetavar)
 192 {
 193     rvec    *rij;
 194     real    *M2, *tt, *dist2;
 195     int      maxnit = 1000;
 196     int      nit    = 0, ll, i, j, type;
 197     t_iatom *ia;
 198     real     L1, tol2, toler;
 199     real     mm    = 0., tmp;
 200     int      error = 0;
 201     real     g, vscale, rscale, rvscale;
 202
 203     if (ncon > shaked->nalloc)
 204     {
 205         shaked->nalloc = over_alloc_dd(ncon);
 206         srenew(shaked->rij, shaked->nalloc);
 207         srenew(shaked->M2, shaked->nalloc);
 208         srenew(shaked->tt, shaked->nalloc);
 209         srenew(shaked->dist2, shaked->nalloc);
 210     }
 211     rij   = shaked->rij;
 212     M2    = shaked->M2;
 213     tt    = shaked->tt;
 214     dist2 = shaked->dist2;
 215
 216     L1   = 1.0-lambda;
 217     tol2 = 2.0*tol;
 218     ia   = iatom;
 219     for (ll = 0; (ll < ncon); ll++, ia += 3)
 220     {
 221         type  = ia[0];
 222         i     = ia[1];
 223         j     = ia[2];
 224
 225         mm          = 2*(invmass[i]+invmass[j]);
 226         rij[ll][XX] = x[i][XX]-x[j][XX];
 227         rij[ll][YY] = x[i][YY]-x[j][YY];
 228         rij[ll][ZZ] = x[i][ZZ]-x[j][ZZ];
 229         M2[ll]      = 1.0/mm;
 230         if (bFEP)
 231         {
 232             toler = sqr(L1*ip[type].constr.dA + lambda*ip[type].constr.dB);
 233         }
 234         else
 235         {
 236             toler = sqr(ip[type].constr.dA);
 237         }
 238         dist2[ll] = toler;
 239         tt[ll]    = 1.0/(toler*tol2);
 240     }
 241
 242     switch (econq)
 243     {
 244         case econqCoord:
 245             cshake(iatom, ncon, &nit, maxnit, dist2, prime[0], rij[0], M2, omega, invmass, tt, lagr, &error);
 246             break;
 247         case econqVeloc:
 248             crattle(iatom, ncon, &nit, maxnit, dist2, prime[0], rij[0], M2, omega, invmass, tt, lagr, &error, invdt, vetavar);
 249             break;
 250     }
 251
 252     if (nit >= maxnit)
 253     {
 254         if (fplog)
 255         {
 256             fprintf(fplog, "Shake did not converge in %d steps\n", maxnit);
 257         }
 258         fprintf(stderr, "Shake did not converge in %d steps\n", maxnit);
 259         nit = 0;
 260     }
 261     else if (error != 0)
 262     {
 263         if (fplog)
 264         {
 265             fprintf(fplog, "Inner product between old and new vector <= 0.0!\n"
 266                     "constraint #%d atoms %u and %u\n",
 267                     error-1, iatom[3*(error-1)+1]+1, iatom[3*(error-1)+2]+1);
 268         }
 269         fprintf(stderr, "Inner product between old and new vector <= 0.0!\n"
 270                 "constraint #%d atoms %u and %u\n",
 271                 error-1, iatom[3*(error-1)+1]+1, iatom[3*(error-1)+2]+1);
 272         nit = 0;
 273     }
 274
 275     /* Constraint virial and correct the lagrange multipliers for the length */
 276
 277     ia = iatom;
 278
 279     for (ll = 0; (ll < ncon); ll++, ia += 3)
 280     {
 281
 282         if ((econq == econqCoord) && v != NULL)
 283         {
 284             /* Correct the velocities */
 285             mm = lagr[ll]*invmass[ia[1]]*invdt/vetavar->rscale;
 286             for (i = 0; i < DIM; i++)
 287             {
 288                 v[ia[1]][i] += mm*rij[ll][i];
 289             }
 290             mm = lagr[ll]*invmass[ia[2]]*invdt/vetavar->rscale;
 291             for (i = 0; i < DIM; i++)
 292             {
 293                 v[ia[2]][i] -= mm*rij[ll][i];
 294             }
 295             /* 16 flops */
 296         }
 297
 298         /* constraint virial */
 299         if (bCalcVir)
 300         {
 301             if (econq == econqCoord)
 302             {
 303                 mm = lagr[ll]/vetavar->rvscale;
 304             }
 305             if (econq == econqVeloc)
 306             {
 307                 mm = lagr[ll]/(vetavar->vscale*vetavar->vscale_nhc[0]);
 308             }
 309             for (i = 0; i < DIM; i++)
 310             {
 311                 tmp = mm*rij[ll][i];
 312                 for (j = 0; j < DIM; j++)
 313                 {
 314                     vir_r_m_dr[i][j] -= tmp*rij[ll][j];
 315                 }
 316             }
 317             /* 21 flops */
 318         }
 319
 320         /* Correct the lagrange multipliers for the length  */
 321         /* (more details would be useful here . . . )*/
 322
 323         type  = ia[0];
 324         if (bFEP)
 325         {
 326             toler = L1*ip[type].constr.dA + lambda*ip[type].constr.dB;
 327         }
 328         else
 329         {
 330             toler     = ip[type].constr.dA;
 331             lagr[ll] *= toler;
 332         }
 333     }
 334
 335     return nit;
 336 }
 337
 338 static void check_cons(FILE *log, int nc, rvec x[], rvec prime[], rvec v[],
 339                        t_iparams ip[], t_iatom *iatom,
 340                        real invmass[], int econq)
 341 {
 342     t_iatom *ia;
 343     int      ai, aj;
 344     int      i;
 345     real     d, dp;
 346     rvec     dx, dv;
 347
 348     fprintf(log,
 349             "    i     mi      j     mj      before       after   should be\n");
 350     ia = iatom;
 351     for (i = 0; (i < nc); i++, ia += 3)
 352     {
 353         ai = ia[1];
 354         aj = ia[2];
 355         rvec_sub(x[ai], x[aj], dx);
 356         d = norm(dx);
 357
 358         switch (econq)
 359         {
 360             case econqCoord:
 361                 rvec_sub(prime[ai], prime[aj], dx);
 362                 dp = norm(dx);
 363                 fprintf(log, "%5d  %5.2f  %5d  %5.2f  %10.5f  %10.5f  %10.5f\n",
 364                         ai+1, 1.0/invmass[ai],
 365                         aj+1, 1.0/invmass[aj], d, dp, ip[ia[0]].constr.dA);
 366                 break;
 367             case econqVeloc:
 368                 rvec_sub(v[ai], v[aj], dv);
 369                 d = iprod(dx, dv);
 370                 rvec_sub(prime[ai], prime[aj], dv);
 371                 dp = iprod(dx, dv);
 372                 fprintf(log, "%5d  %5.2f  %5d  %5.2f  %10.5f  %10.5f  %10.5f\n",
 373                         ai+1, 1.0/invmass[ai],
 374                         aj+1, 1.0/invmass[aj], d, dp, 0.);
 375                 break;
 376         }
 377     }
 378 }
 379
 380 gmx_bool bshakef(FILE *log, gmx_shakedata_t shaked,
 381                  real invmass[], int nblocks, int sblock[],
 382                  t_idef *idef, t_inputrec *ir, rvec x_s[], rvec prime[],
 383                  t_nrnb *nrnb, real *lagr, real lambda, real *dvdlambda,
 384                  real invdt, rvec *v, gmx_bool bCalcVir, tensor vir_r_m_dr,
 385                  gmx_bool bDumpOnError, int econq, t_vetavars *vetavar)
 386 {
 387     t_iatom *iatoms;
 388     real    *lam, dt_2, dvdl;
 389     int      i, n0, ncons, blen, type;
 390     int      tnit = 0, trij = 0;
 391
 392 #ifdef DEBUG
 393     fprintf(log, "nblocks=%d, sblock[0]=%d\n", nblocks, sblock[0]);
 394 #endif
 395
 396     ncons = idef->il[F_CONSTR].nr/3;
 397
 398     for (i = 0; i < ncons; i++)
 399     {
 400         lagr[i] = 0;
 401     }
 402
 403     iatoms = &(idef->il[F_CONSTR].iatoms[sblock[0]]);
 404     lam    = lagr;
 405     for (i = 0; (i < nblocks); )
 406     {
 407         blen  = (sblock[i+1]-sblock[i]);
 408         blen /= 3;
 409         n0    = vec_shakef(log, shaked, invmass, blen, idef->iparams,
 410                            iatoms, ir->shake_tol, x_s, prime, shaked->omega,
 411                            ir->efep != efepNO, lambda, lam, invdt, v, bCalcVir, vir_r_m_dr,
 412                            econq, vetavar);
 413
 414 #ifdef DEBUGSHAKE
 415         check_cons(log, blen, x_s, prime, v, idef->iparams, iatoms, invmass, econq);
 416 #endif
 417
 418         if (n0 == 0)
 419         {
 420             if (bDumpOnError && log)
 421             {
 422                 {
 423                     check_cons(log, blen, x_s, prime, v, idef->iparams, iatoms, invmass, econq);
 424                 }
 425             }
 426             return FALSE;
 427         }
 428         tnit   += n0*blen;
 429         trij   += blen;
 430         iatoms += 3*blen; /* Increment pointer! */
 431         lam    += blen;
 432         i++;
 433     }
 434     /* only for position part? */
 435     if (econq == econqCoord)
 436     {
 437         if (ir->efep != efepNO)
 438         {
 439             real bondA, bondB;
 440             dt_2 = 1/sqr(ir->delta_t);
 441             dvdl = 0;
 442             for (i = 0; i < ncons; i++)
 443             {
 444                 type  = idef->il[F_CONSTR].iatoms[3*i];
 445
 446                 /* dh/dl contribution from constraint force is  dh/dr (constraint force) dot dr/dl */
 447                 /* constraint force is -\sum_i lagr_i* d(constraint)/dr, with constrant = r^2-d^2  */
 448                 /* constraint force is -\sum_i lagr_i* 2 r  */
 449                 /* so dh/dl = -\sum_i lagr_i* 2 r * dr/dl */
 450                 /* However, by comparison with lincs and with
 451                    comparison with a full thermodynamics cycle (see
 452                    redmine issue #1255), this is off by a factor of
 453                    two -- the 2r should apparently just be r.  Further
 454                    investigation should be done at some point to
 455                    understand why and see if there is something deeper
 456                    we are missing */
 457
 458                 bondA = idef->iparams[type].constr.dA;
 459                 bondB = idef->iparams[type].constr.dB;
 460                 dvdl += lagr[i] * dt_2 * ((1.0-lambda)*bondA + lambda*bondB) * (bondB-bondA);
 461             }
 462             *dvdlambda += dvdl;
 463         }
 464     }
 465 #ifdef DEBUG
 466     fprintf(log, "tnit: %5d  omega: %10.5f\n", tnit, omega);
 467 #endif
 468     if (ir->bShakeSOR)
 469     {
 470         if (tnit > shaked->gamma)
 471         {
 472             shaked->delta *= -0.5;
 473         }
 474         shaked->omega += shaked->delta;
 475         shaked->gamma  = tnit;
 476     }
 477     inc_nrnb(nrnb, eNR_SHAKE, tnit);
 478     inc_nrnb(nrnb, eNR_SHAKE_RIJ, trij);
 479     if (v)
 480     {
 481         inc_nrnb(nrnb, eNR_CONSTR_V, trij*2);
 482     }
 483     if (bCalcVir)
 484     {
 485         inc_nrnb(nrnb, eNR_CONSTR_VIR, trij);
 486     }
 487
 488     return TRUE;
 489 }
 490
 491 void crattle(atom_id iatom[], int ncon, int *nnit, int maxnit,
 492              real dist2[], real vp[], real rij[], real m2[], real omega,
 493              real invmass[], real tt[], real lagr[], int *nerror, real invdt, t_vetavars *vetavar)
 494 {
 495     /*
 496      *     r.c. van schaik and w.f. van gunsteren
 497      *     eth zuerich
 498      *     june 1992
 499      *     Adapted for use with Gromacs by David van der Spoel november 92 and later.
 500      *     rattle added by M.R. Shirts, April 2004, from code written by Jay Ponder in TINKER
 501      *     second part of rattle algorithm
 502      */
 503
 504     const   real mytol = 1e-10;
 505
 506     int          ll, i, j, i3, j3, l3, ii;
 507     int          ix, iy, iz, jx, jy, jz;
 508     real         toler, rijd, vpijd, vx, vy, vz, diff, acor, xdotd, fac, im, jm, imdt, jmdt;
 509     real         xh, yh, zh, rijx, rijy, rijz;
 510     real         tix, tiy, tiz;
 511     real         tjx, tjy, tjz;
 512     int          nit, error, nconv;
 513     real         veta, vscale_nhc, iconvf;
 514
 515     veta       = vetavar->veta;
 516     vscale_nhc = vetavar->vscale_nhc[0];  /* for now, just use the first state */
 517
 518     error = 0;
 519     nconv = 1;
 520     for (nit = 0; (nit < maxnit) && (nconv != 0) && (error == 0); nit++)
 521     {
 522         nconv = 0;
 523         for (ll = 0; (ll < ncon) && (error == 0); ll++)
 524         {
 525             l3      = 3*ll;
 526             rijx    = rij[l3+XX];
 527             rijy    = rij[l3+YY];
 528             rijz    = rij[l3+ZZ];
 529             i       = iatom[l3+1];
 530             j       = iatom[l3+2];
 531             i3      = 3*i;
 532             j3      = 3*j;
 533             ix      = i3+XX;
 534             iy      = i3+YY;
 535             iz      = i3+ZZ;
 536             jx      = j3+XX;
 537             jy      = j3+YY;
 538             jz      = j3+ZZ;
 539             vx      = vp[ix]-vp[jx];
 540             vy      = vp[iy]-vp[jy];
 541             vz      = vp[iz]-vp[jz];
 542
 543             vpijd   = vx*rijx+vy*rijy+vz*rijz;
 544             toler   = dist2[ll];
 545             /* this is r(t+dt) \dotproduct \dot{r}(t+dt) */
 546             xdotd   = vpijd*vscale_nhc + veta*toler;
 547
 548             /* iconv is zero when the error is smaller than a bound */
 549             iconvf   = fabs(xdotd)*(tt[ll]/invdt);
 550
 551             if (iconvf > 1)
 552             {
 553                 nconv     = iconvf;
 554                 fac       = omega*2.0*m2[ll]/toler;
 555                 acor      = -fac*xdotd;
 556                 lagr[ll] += acor;
 557
 558                 xh        = rijx*acor;
 559                 yh        = rijy*acor;
 560                 zh        = rijz*acor;
 561
 562                 im        = invmass[i]/vscale_nhc;
 563                 jm        = invmass[j]/vscale_nhc;
 564
 565                 vp[ix] += xh*im;
 566                 vp[iy] += yh*im;
 567                 vp[iz] += zh*im;
 568                 vp[jx] -= xh*jm;
 569                 vp[jy] -= yh*jm;
 570                 vp[jz] -= zh*jm;
 571             }
 572         }
 573     }
 574     *nnit   = nit;
 575     *nerror = error;
 576 }