src/gromacs/mdlib/shakef.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 #ifdef HAVE_CONFIG_H
  37 #include <config.h>
  38 #endif
  39
  40 #include <math.h>
  41 #include "sysstuff.h"
  42 #include "typedefs.h"
  43 #include "smalloc.h"
  44 #include "pbc.h"
  45 #include "txtdump.h"
  46 #include "vec.h"
  47 #include "nrnb.h"
  48 #include "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                int natoms, 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                  int natoms, 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, natoms, 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             dt_2 = 1/sqr(ir->delta_t);
 440             dvdl = 0;
 441             for (i = 0; i < ncons; i++)
 442             {
 443                 type  = idef->il[F_CONSTR].iatoms[3*i];
 444                 dvdl += lagr[i]*dt_2*
 445                     (idef->iparams[type].constr.dB-idef->iparams[type].constr.dA);
 446             }
 447             *dvdlambda += dvdl;
 448         }
 449     }
 450 #ifdef DEBUG
 451     fprintf(log, "tnit: %5d  omega: %10.5f\n", tnit, omega);
 452 #endif
 453     if (ir->bShakeSOR)
 454     {
 455         if (tnit > shaked->gamma)
 456         {
 457             shaked->delta *= -0.5;
 458         }
 459         shaked->omega += shaked->delta;
 460         shaked->gamma  = tnit;
 461     }
 462     inc_nrnb(nrnb, eNR_SHAKE, tnit);
 463     inc_nrnb(nrnb, eNR_SHAKE_RIJ, trij);
 464     if (v)
 465     {
 466         inc_nrnb(nrnb, eNR_CONSTR_V, trij*2);
 467     }
 468     if (bCalcVir)
 469     {
 470         inc_nrnb(nrnb, eNR_CONSTR_VIR, trij);
 471     }
 472
 473     return TRUE;
 474 }
 475
 476 void crattle(atom_id iatom[], int ncon, int *nnit, int maxnit,
 477              real dist2[], real vp[], real rij[], real m2[], real omega,
 478              real invmass[], real tt[], real lagr[], int *nerror, real invdt, t_vetavars *vetavar)
 479 {
 480     /*
 481      *     r.c. van schaik and w.f. van gunsteren
 482      *     eth zuerich
 483      *     june 1992
 484      *     Adapted for use with Gromacs by David van der Spoel november 92 and later.
 485      *     rattle added by M.R. Shirts, April 2004, from code written by Jay Ponder in TINKER
 486      *     second part of rattle algorithm
 487      */
 488
 489     const   real mytol = 1e-10;
 490
 491     int          ll, i, j, i3, j3, l3, ii;
 492     int          ix, iy, iz, jx, jy, jz;
 493     real         toler, rijd, vpijd, vx, vy, vz, diff, acor, xdotd, fac, im, jm, imdt, jmdt;
 494     real         xh, yh, zh, rijx, rijy, rijz;
 495     real         tix, tiy, tiz;
 496     real         tjx, tjy, tjz;
 497     int          nit, error, nconv;
 498     real         veta, vscale_nhc, iconvf;
 499
 500     veta       = vetavar->veta;
 501     vscale_nhc = vetavar->vscale_nhc[0];  /* for now, just use the first state */
 502
 503     error = 0;
 504     nconv = 1;
 505     for (nit = 0; (nit < maxnit) && (nconv != 0) && (error == 0); nit++)
 506     {
 507         nconv = 0;
 508         for (ll = 0; (ll < ncon) && (error == 0); ll++)
 509         {
 510             l3      = 3*ll;
 511             rijx    = rij[l3+XX];
 512             rijy    = rij[l3+YY];
 513             rijz    = rij[l3+ZZ];
 514             i       = iatom[l3+1];
 515             j       = iatom[l3+2];
 516             i3      = 3*i;
 517             j3      = 3*j;
 518             ix      = i3+XX;
 519             iy      = i3+YY;
 520             iz      = i3+ZZ;
 521             jx      = j3+XX;
 522             jy      = j3+YY;
 523             jz      = j3+ZZ;
 524             vx      = vp[ix]-vp[jx];
 525             vy      = vp[iy]-vp[jy];
 526             vz      = vp[iz]-vp[jz];
 527
 528             vpijd   = vx*rijx+vy*rijy+vz*rijz;
 529             toler   = dist2[ll];
 530             /* this is r(t+dt) \dotproduct \dot{r}(t+dt) */
 531             xdotd   = vpijd*vscale_nhc + veta*toler;
 532
 533             /* iconv is zero when the error is smaller than a bound */
 534             iconvf   = fabs(xdotd)*(tt[ll]/invdt);
 535
 536             if (iconvf > 1)
 537             {
 538                 nconv     = iconvf;
 539                 fac       = omega*2.0*m2[ll]/toler;
 540                 acor      = -fac*xdotd;
 541                 lagr[ll] += acor;
 542
 543                 xh        = rijx*acor;
 544                 yh        = rijy*acor;
 545                 zh        = rijz*acor;
 546
 547                 im        = invmass[i]/vscale_nhc;
 548                 jm        = invmass[j]/vscale_nhc;
 549
 550                 vp[ix] += xh*im;
 551                 vp[iy] += yh*im;
 552                 vp[iz] += zh*im;
 553                 vp[jx] -= xh*jm;
 554                 vp[jy] -= yh*jm;
 555                 vp[jz] -= zh*jm;
 556             }
 557         }
 558     }
 559     *nnit   = nit;
 560     *nerror = error;
 561 }