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     fprintf(log,"  %10.3f",x[m]);
  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     nconv=0;
 120     for(ll=0; (ll<ncon) && (error == 0); ll++) {
 121       l3    = 3*ll;
 122       rijx  = rij[l3+XX];
 123       rijy  = rij[l3+YY];
 124       rijz  = rij[l3+ZZ];
 125       i     = iatom[l3+1];
 126       j     = iatom[l3+2];
 127       i3    = 3*i;
 128       j3    = 3*j;
 129       ix    = i3+XX;
 130       iy    = i3+YY;
 131       iz    = i3+ZZ;
 132       jx    = j3+XX;
 133       jy    = j3+YY;
 134       jz    = j3+ZZ;
 135
 136       tx      = xp[ix]-xp[jx];
 137       ty      = xp[iy]-xp[jy];
 138       tz      = xp[iz]-xp[jz];
 139       rpij2   = tx*tx+ty*ty+tz*tz;
 140       toler   = dist2[ll];
 141       diff    = toler-rpij2;
 142
 143       /* iconvf is less than 1 when the error is smaller than a bound */
 144       /* But if tt is too big, then it will result in looping in iconv */
 145
 146       iconvf = fabs(diff)*tt[ll];
 147
 148       if (iconvf > 1) {
 149           nconv   = iconvf;
 150           rrpr    = rijx*tx+rijy*ty+rijz*tz;
 151
 152           if (rrpr < toler*mytol)
 153               error=ll;
 154           else {
 155               acor      = omega*diff*m2[ll]/rrpr;
 156               lagr[ll] += acor;
 157               xh        = rijx*acor;
 158               yh        = rijy*acor;
 159               zh        = rijz*acor;
 160               im        = invmass[i];
 161               jm        = invmass[j];
 162               xp[ix] += xh*im;
 163               xp[iy] += yh*im;
 164               xp[iz] += zh*im;
 165               xp[jx] -= xh*jm;
 166               xp[jy] -= yh*jm;
 167               xp[jz] -= zh*jm;
 168           }
 169       }
 170     }
 171   }
 172   *nnit=nit;
 173   *nerror=error;
 174 }
 175
 176 int vec_shakef(FILE *fplog,gmx_shakedata_t shaked,
 177                int natoms,real invmass[],int ncon,
 178                t_iparams ip[],t_iatom *iatom,
 179                real tol,rvec x[],rvec prime[],real omega,
 180                gmx_bool bFEP,real lambda,real lagr[],
 181                real invdt,rvec *v,
 182                gmx_bool bCalcVir,tensor rmdr,int econq,
 183                t_vetavars *vetavar)
 184 {
 185     rvec *rij;
 186     real *M2,*tt,*dist2;
 187     int     maxnit=1000;
 188     int     nit=0,ll,i,j,type;
 189     t_iatom *ia;
 190     real    L1,tol2,toler;
 191     real    mm=0.,tmp;
 192     int     error=0;
 193     real    g,vscale,rscale,rvscale;
 194
 195     if (ncon > shaked->nalloc)
 196     {
 197         shaked->nalloc = over_alloc_dd(ncon);
 198         srenew(shaked->rij,shaked->nalloc);
 199         srenew(shaked->M2,shaked->nalloc);
 200         srenew(shaked->tt,shaked->nalloc);
 201         srenew(shaked->dist2,shaked->nalloc);
 202     }
 203     rij   = shaked->rij;
 204     M2    = shaked->M2;
 205     tt    = shaked->tt;
 206     dist2 = shaked->dist2;
 207
 208     L1=1.0-lambda;
 209     tol2=2.0*tol;
 210     ia=iatom;
 211     for(ll=0; (ll<ncon); ll++,ia+=3) {
 212         type  = ia[0];
 213         i=ia[1];
 214         j=ia[2];
 215
 216         mm=2*(invmass[i]+invmass[j]);
 217         rij[ll][XX]=x[i][XX]-x[j][XX];
 218         rij[ll][YY]=x[i][YY]-x[j][YY];
 219         rij[ll][ZZ]=x[i][ZZ]-x[j][ZZ];
 220         M2[ll]=1.0/mm;
 221         if (bFEP)
 222             toler = sqr(L1*ip[type].constr.dA + lambda*ip[type].constr.dB);
 223         else
 224             toler = sqr(ip[type].constr.dA);
 225         dist2[ll] = toler;
 226         tt[ll] = 1.0/(toler*tol2);
 227     }
 228
 229     switch (econq) {
 230     case econqCoord:
 231         cshake(iatom,ncon,&nit,maxnit,dist2,prime[0],rij[0],M2,omega,invmass,tt,lagr,&error);
 232         break;
 233     case econqVeloc:
 234         crattle(iatom,ncon,&nit,maxnit,dist2,prime[0],rij[0],M2,omega,invmass,tt,lagr,&error,invdt,vetavar);
 235         break;
 236     }
 237
 238     if (nit >= maxnit) {
 239         if (fplog) {
 240             fprintf(fplog,"Shake did not converge in %d steps\n",maxnit);
 241         }
 242         fprintf(stderr,"Shake did not converge in %d steps\n",maxnit);
 243         nit=0;
 244     }
 245     else if (error != 0)
 246     {
 247         if (fplog)
 248             fprintf(fplog,"Inner product between old and new vector <= 0.0!\n"
 249                     "constraint #%d atoms %u and %u\n",
 250                     error-1,iatom[3*(error-1)+1]+1,iatom[3*(error-1)+2]+1);
 251         fprintf(stderr,"Inner product between old and new vector <= 0.0!\n"
 252                 "constraint #%d atoms %u and %u\n",
 253                 error-1,iatom[3*(error-1)+1]+1,iatom[3*(error-1)+2]+1);
 254         nit=0;
 255     }
 256
 257     /* Constraint virial and correct the lagrange multipliers for the length */
 258
 259     ia=iatom;
 260
 261     for(ll=0; (ll<ncon); ll++,ia+=3)
 262     {
 263
 264         if ((econq == econqCoord) && v!=NULL)
 265         {
 266             /* Correct the velocities */
 267             mm = lagr[ll]*invmass[ia[1]]*invdt/vetavar->rscale;
 268             for(i=0; i<DIM; i++)
 269             {
 270                 v[ia[1]][i] += mm*rij[ll][i];
 271             }
 272             mm = lagr[ll]*invmass[ia[2]]*invdt/vetavar->rscale;
 273             for(i=0; i<DIM; i++)
 274             {
 275                 v[ia[2]][i] -= mm*rij[ll][i];
 276             }
 277             /* 16 flops */
 278         }
 279
 280         /* constraint virial */
 281         if (bCalcVir)
 282         {
 283             if (econq == econqCoord)
 284             {
 285                 mm = lagr[ll]/vetavar->rvscale;
 286             }
 287             if (econq == econqVeloc)
 288             {
 289                 mm = lagr[ll]/(vetavar->vscale*vetavar->vscale_nhc[0]);
 290             }
 291             for(i=0; i<DIM; i++)
 292             {
 293                 tmp = mm*rij[ll][i];
 294                 for(j=0; j<DIM; j++)
 295                 {
 296                     rmdr[i][j] -= tmp*rij[ll][j];
 297                 }
 298             }
 299             /* 21 flops */
 300         }
 301
 302         /* Correct the lagrange multipliers for the length  */
 303         /* (more details would be useful here . . . )*/
 304
 305         type  = ia[0];
 306         if (bFEP)
 307         {
 308             toler = L1*ip[type].constr.dA + lambda*ip[type].constr.dB;
 309         }
 310         else
 311         {
 312             toler = ip[type].constr.dA;
 313             lagr[ll] *= toler;
 314         }
 315     }
 316
 317     return nit;
 318 }
 319
 320 static void check_cons(FILE *log,int nc,rvec x[],rvec prime[], rvec v[],
 321                        t_iparams ip[],t_iatom *iatom,
 322                        real invmass[], int econq)
 323 {
 324   t_iatom *ia;
 325   int     ai,aj;
 326   int     i;
 327   real    d,dp;
 328   rvec    dx,dv;
 329
 330   fprintf(log,
 331           "    i     mi      j     mj      before       after   should be\n");
 332   ia=iatom;
 333   for(i=0; (i<nc); i++,ia+=3) {
 334     ai=ia[1];
 335     aj=ia[2];
 336     rvec_sub(x[ai],x[aj],dx);
 337     d=norm(dx);
 338
 339     switch (econq) {
 340     case econqCoord:
 341         rvec_sub(prime[ai],prime[aj],dx);
 342         dp=norm(dx);
 343         fprintf(log,"%5d  %5.2f  %5d  %5.2f  %10.5f  %10.5f  %10.5f\n",
 344                 ai+1,1.0/invmass[ai],
 345                 aj+1,1.0/invmass[aj],d,dp,ip[ia[0]].constr.dA);
 346         break;
 347     case econqVeloc:
 348         rvec_sub(v[ai],v[aj],dv);
 349         d=iprod(dx,dv);
 350         rvec_sub(prime[ai],prime[aj],dv);
 351         dp=iprod(dx,dv);
 352         fprintf(log,"%5d  %5.2f  %5d  %5.2f  %10.5f  %10.5f  %10.5f\n",
 353                 ai+1,1.0/invmass[ai],
 354                 aj+1,1.0/invmass[aj],d,dp,0.);
 355         break;
 356     }
 357   }
 358 }
 359
 360 gmx_bool bshakef(FILE *log,gmx_shakedata_t shaked,
 361              int natoms,real invmass[],int nblocks,int sblock[],
 362              t_idef *idef,t_inputrec *ir,matrix box,rvec x_s[],rvec prime[],
 363              t_nrnb *nrnb,real *lagr,real lambda,real *dvdlambda,
 364              real invdt,rvec *v,gmx_bool bCalcVir,tensor rmdr,gmx_bool bDumpOnError,int econq,t_vetavars *vetavar)
 365 {
 366   t_iatom *iatoms;
 367   real    *lam,dt_2,dvdl;
 368   int     i,n0,ncons,blen,type;
 369   int     tnit=0,trij=0;
 370
 371 #ifdef DEBUG
 372   fprintf(log,"nblocks=%d, sblock[0]=%d\n",nblocks,sblock[0]);
 373 #endif
 374
 375   ncons=idef->il[F_CONSTR].nr/3;
 376
 377   for(i=0; i<ncons; i++)
 378     lagr[i] =0;
 379
 380   iatoms = &(idef->il[F_CONSTR].iatoms[sblock[0]]);
 381   lam    = lagr;
 382   for(i=0; (i<nblocks); ) {
 383     blen  = (sblock[i+1]-sblock[i]);
 384     blen /= 3;
 385     n0 = vec_shakef(log,shaked,natoms,invmass,blen,idef->iparams,
 386                     iatoms,ir->shake_tol,x_s,prime,shaked->omega,
 387                     ir->efep!=efepNO,lambda,lam,invdt,v,bCalcVir,rmdr,econq,vetavar);
 388
 389 #ifdef DEBUGSHAKE
 390     check_cons(log,blen,x_s,prime,v,idef->iparams,iatoms,invmass,econq);
 391 #endif
 392
 393     if (n0 == 0) {
 394         if (bDumpOnError && log)
 395         {
 396             {
 397                 check_cons(log,blen,x_s,prime,v,idef->iparams,iatoms,invmass,econq);
 398             }
 399         }
 400         return FALSE;
 401     }
 402     tnit   += n0*blen;
 403     trij   += blen;
 404     iatoms += 3*blen;   /* Increment pointer! */
 405     lam    += blen;
 406     i++;
 407   }
 408   /* only for position part? */
 409   if (econq == econqCoord) {
 410       if (ir->efep != efepNO) {
 411           dt_2 = 1/sqr(ir->delta_t);
 412           dvdl = 0;
 413           for(i=0; i<ncons; i++) {
 414               type = idef->il[F_CONSTR].iatoms[3*i];
 415               dvdl += lagr[i]*dt_2*
 416                   (idef->iparams[type].constr.dB-idef->iparams[type].constr.dA);
 417           }
 418           *dvdlambda += dvdl;
 419       }
 420   }
 421 #ifdef DEBUG
 422   fprintf(log,"tnit: %5d  omega: %10.5f\n",tnit,omega);
 423 #endif
 424   if (ir->bShakeSOR) {
 425     if (tnit > shaked->gamma) {
 426       shaked->delta *= -0.5;
 427     }
 428     shaked->omega += shaked->delta;
 429     shaked->gamma  = tnit;
 430   }
 431   inc_nrnb(nrnb,eNR_SHAKE,tnit);
 432   inc_nrnb(nrnb,eNR_SHAKE_RIJ,trij);
 433   if (v)
 434     inc_nrnb(nrnb,eNR_CONSTR_V,trij*2);
 435   if (bCalcVir)
 436     inc_nrnb(nrnb,eNR_CONSTR_VIR,trij);
 437
 438   return TRUE;
 439 }
 440
 441 void crattle(atom_id iatom[],int ncon,int *nnit,int maxnit,
 442              real dist2[],real vp[],real rij[],real m2[],real omega,
 443              real invmass[],real tt[],real lagr[],int *nerror,real invdt,t_vetavars *vetavar)
 444 {
 445     /*
 446      *     r.c. van schaik and w.f. van gunsteren
 447      *     eth zuerich
 448      *     june 1992
 449      *     Adapted for use with Gromacs by David van der Spoel november 92 and later.
 450      *     rattle added by M.R. Shirts, April 2004, from code written by Jay Ponder in TINKER
 451      *     second part of rattle algorithm
 452      */
 453
 454     const   real mytol=1e-10;
 455
 456     int     ll,i,j,i3,j3,l3,ii;
 457     int     ix,iy,iz,jx,jy,jz;
 458     real    toler,rijd,vpijd, vx,vy,vz,diff,acor,xdotd,fac,im,jm,imdt,jmdt;
 459     real    xh,yh,zh,rijx,rijy,rijz;
 460     real    tix,tiy,tiz;
 461     real    tjx,tjy,tjz;
 462     int     nit,error,nconv;
 463     real    veta,vscale_nhc,iconvf;
 464
 465     veta = vetavar->veta;
 466     vscale_nhc = vetavar->vscale_nhc[0];  /* for now, just use the first state */
 467
 468     error=0;
 469     nconv=1;
 470     for (nit=0; (nit<maxnit) && (nconv != 0) && (error == 0); nit++) {
 471         nconv=0;
 472         for(ll=0; (ll<ncon) && (error == 0); ll++) {
 473             l3    = 3*ll;
 474             rijx  = rij[l3+XX];
 475             rijy  = rij[l3+YY];
 476             rijz  = rij[l3+ZZ];
 477             i     = iatom[l3+1];
 478             j     = iatom[l3+2];
 479             i3    = 3*i;
 480             j3    = 3*j;
 481             ix    = i3+XX;
 482             iy    = i3+YY;
 483             iz    = i3+ZZ;
 484             jx    = j3+XX;
 485             jy    = j3+YY;
 486             jz    = j3+ZZ;
 487             vx      = vp[ix]-vp[jx];
 488             vy      = vp[iy]-vp[jy];
 489             vz      = vp[iz]-vp[jz];
 490
 491             vpijd   = vx*rijx+vy*rijy+vz*rijz;
 492             toler   = dist2[ll];
 493             /* this is r(t+dt) \dotproduct \dot{r}(t+dt) */
 494             xdotd   = vpijd*vscale_nhc + veta*toler;
 495
 496             /* iconv is zero when the error is smaller than a bound */
 497             iconvf   = fabs(xdotd)*(tt[ll]/invdt);
 498
 499             if (iconvf > 1) {
 500                 nconv     = iconvf;
 501                 fac       = omega*2.0*m2[ll]/toler;
 502                 acor      = -fac*xdotd;
 503                 lagr[ll] += acor;
 504
 505                 xh        = rijx*acor;
 506                 yh        = rijy*acor;
 507                 zh        = rijz*acor;
 508
 509                 im        = invmass[i]/vscale_nhc;
 510                 jm        = invmass[j]/vscale_nhc;
 511
 512                 vp[ix] += xh*im;
 513                 vp[iy] += yh*im;
 514                 vp[iz] += zh*im;
 515                 vp[jx] -= xh*jm;
 516                 vp[jy] -= yh*jm;
 517                 vp[jz] -= zh*jm;
 518             }
 519         }
 520     }
 521     *nnit=nit;
 522     *nerror=error;
 523 }