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