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