src/mdlib/dummies.c

   1 /*
   2  * $Id$
   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 2.0
  11  *
  12  * Copyright (c) 1991-1999
  13  * BIOSON Research Institute, Dept. of Biophysical Chemistry
  14  * University of Groningen, The Netherlands
  15  *
  16  * Please refer to:
  17  * GROMACS: A message-passing parallel molecular dynamics implementation
  18  * H.J.C. Berendsen, D. van der Spoel and R. van Drunen
  19  * Comp. Phys. Comm. 91, 43-56 (1995)
  20  *
  21  * Also check out our WWW page:
  22  * http://md.chem.rug.nl/~gmx
  23  * or e-mail to:
  24  * gromacs@chem.rug.nl
  25  *
  26  * And Hey:
  27  * GRowing Old MAkes el Chrono Sweat
  28  */
  29 static char *SRCID_dummies_c = "$Id$";
  30
  31 #include <stdio.h>
  32 #include "typedefs.h"
  33 #include "assert.h"
  34 #include "dummies.h"
  35 #include "macros.h"
  36 #include "smalloc.h"
  37 #include "nrnb.h"
  38 #include "vec.h"
  39
  40 static void constr_dum2(rvec xi,rvec xj,rvec x,real a)
  41 {
  42   real b;
  43
  44   b=1.0-a;
  45   /* 1 flop */
  46
  47   x[XX]=b*xi[XX]+a*xj[XX];
  48   x[YY]=b*xi[YY]+a*xj[YY];
  49   x[ZZ]=b*xi[ZZ]+a*xj[ZZ];
  50   /* 9 Flops */
  51
  52   /* TOTAL: 10 flops */
  53 }
  54
  55 static void constr_dum3(rvec xi,rvec xj,rvec xk,rvec x,real a,real b)
  56 {
  57   real c;
  58
  59   c=1.0-a-b;
  60   /* 2 flops */
  61
  62   x[XX] = c*xi[XX] + a*xj[XX] + b*xk[XX];
  63   x[YY] = c*xi[YY] + a*xj[YY] + b*xk[YY];
  64   x[ZZ] = c*xi[ZZ] + a*xj[ZZ] + b*xk[ZZ];
  65   /* 15 Flops */
  66
  67   /* TOTAL: 17 flops */
  68 }
  69
  70 static void constr_dum3FD(rvec xi,rvec xj,rvec xk,rvec x,real a,real b)
  71 {
  72   rvec xij,xjk,temp;
  73   real c;
  74
  75   rvec_sub(xj,xi,xij);
  76   rvec_sub(xk,xj,xjk);
  77   /* 6 flops */
  78
  79   /* temp goes from i to a point on the line jk */
  80   temp[XX] = xij[XX] + a*xjk[XX];
  81   temp[YY] = xij[YY] + a*xjk[YY];
  82   temp[ZZ] = xij[ZZ] + a*xjk[ZZ];
  83   /* 6 flops */
  84
  85   c=b*invsqrt(iprod(temp,temp));
  86   /* 6 + 10 flops */
  87
  88   x[XX] = xi[XX] + c*temp[XX];
  89   x[YY] = xi[YY] + c*temp[YY];
  90   x[ZZ] = xi[ZZ] + c*temp[ZZ];
  91   /* 6 Flops */
  92
  93   /* TOTAL: 34 flops */
  94 }
  95
  96 static void constr_dum3FAD(rvec xi,rvec xj,rvec xk,rvec x,real a,real b)
  97 {
  98   rvec xij,xjk,xp;
  99   real a1,b1,c1,invdij;
 100
 101   rvec_sub(xj,xi,xij);
 102   rvec_sub(xk,xj,xjk);
 103   /* 6 flops */
 104
 105   invdij = invsqrt(iprod(xij,xij));
 106   c1 = invdij * invdij * iprod(xij,xjk);
 107   xp[XX] = xjk[XX] - c1*xij[XX];
 108   xp[YY] = xjk[YY] - c1*xij[YY];
 109   xp[ZZ] = xjk[ZZ] - c1*xij[ZZ];
 110   a1 = a*invdij;
 111   b1 = b*invsqrt(iprod(xp,xp));
 112   /* 45 */
 113
 114   x[XX] = xi[XX] + a1*xij[XX] + b1*xp[XX];
 115   x[YY] = xi[YY] + a1*xij[YY] + b1*xp[YY];
 116   x[ZZ] = xi[ZZ] + a1*xij[ZZ] + b1*xp[ZZ];
 117   /* 12 Flops */
 118
 119   /* TOTAL: 63 flops */
 120 }
 121
 122 static void constr_dum3OUT(rvec xi,rvec xj,rvec xk,rvec x,real a,real b,real c)
 123 {
 124   rvec xij,xik,temp;
 125
 126   rvec_sub(xj,xi,xij);
 127   rvec_sub(xk,xi,xik);
 128   oprod(xij,xik,temp);
 129   /* 15 Flops */
 130
 131   x[XX] = xi[XX] + a*xij[XX] + b*xik[XX] + c*temp[XX];
 132   x[YY] = xi[YY] + a*xij[YY] + b*xik[YY] + c*temp[YY];
 133   x[ZZ] = xi[ZZ] + a*xij[ZZ] + b*xik[ZZ] + c*temp[ZZ];
 134   /* 18 Flops */
 135
 136   /* TOTAL: 33 flops */
 137 }
 138
 139 static void constr_dum4FD(rvec xi,rvec xj,rvec xk,rvec xl,rvec x,
 140                           real a,real b,real c)
 141 {
 142   rvec xij,xjk,xjl,temp;
 143   real d;
 144
 145   rvec_sub(xj,xi,xij);
 146   rvec_sub(xk,xj,xjk);
 147   rvec_sub(xl,xj,xjl);
 148   /* 9 flops */
 149
 150   /* temp goes from i to a point on the plane jkl */
 151   temp[XX] = xij[XX] + a*xjk[XX] + b*xjl[XX];
 152   temp[YY] = xij[YY] + a*xjk[YY] + b*xjl[YY];
 153   temp[ZZ] = xij[ZZ] + a*xjk[ZZ] + b*xjl[ZZ];
 154   /* 12 flops */
 155
 156   d=c*invsqrt(iprod(temp,temp));
 157   /* 6 + 10 flops */
 158
 159   x[XX] = xi[XX] + d*temp[XX];
 160   x[YY] = xi[YY] + d*temp[YY];
 161   x[ZZ] = xi[ZZ] + d*temp[ZZ];
 162   /* 6 Flops */
 163
 164   /* TOTAL: 43 flops */
 165 }
 166
 167 void construct_dummies(FILE *log,rvec x[],t_nrnb *nrnb,real dt,
 168                        rvec *v,t_idef *idef)
 169 {
 170   rvec      xd,vv;
 171   real      a1,b1,c1,inv_dt;
 172   int       i,nra,nrd,tp,ftype;
 173   t_iatom   adum,ai,aj,ak,al;
 174   t_iatom   *ia;
 175   t_iparams *ip;
 176
 177   ip     = idef->iparams;
 178   if (v)
 179     inv_dt = 1.0/dt;
 180   else
 181     inv_dt = 1.0;
 182
 183   for(ftype=0; (ftype<F_NRE); ftype++) {
 184     if (interaction_function[ftype].flags & IF_DUMMY) {
 185       nra    = interaction_function[ftype].nratoms;
 186       nrd    = idef->il[ftype].nr;
 187       ia     = idef->il[ftype].iatoms;
 188
 189       for(i=0; (i<nrd); ) {
 190         tp   = ia[0];
 191         assert(ftype == idef->functype[tp]);
 192
 193         /* The dummy and constructing atoms */
 194         adum = ia[1];
 195         ai   = ia[2];
 196         aj   = ia[3];
 197
 198         /* Constants for constructing dummies */
 199         a1   = ip[tp].dummy.a;
 200
 201         /* Copy the old position */
 202         copy_rvec(x[adum],xd);
 203
 204         /* Construct the dummy depending on type */
 205         switch (ftype) {
 206         case F_DUMMY2:
 207           constr_dum2(x[ai],x[aj],x[adum],a1);
 208           break;
 209         case F_DUMMY3:
 210           ak = ia[4];
 211           b1 = ip[tp].dummy.b;
 212           constr_dum3(x[ai],x[aj],x[ak],x[adum],a1,b1);
 213           break;
 214         case F_DUMMY3FD:
 215           ak = ia[4];
 216           b1 = ip[tp].dummy.b;
 217           constr_dum3FD(x[ai],x[aj],x[ak],x[adum],a1,b1);
 218           break;
 219         case F_DUMMY3FAD:
 220           ak = ia[4];
 221           b1 = ip[tp].dummy.b;
 222           constr_dum3FAD(x[ai],x[aj],x[ak],x[adum],a1,b1);
 223           break;
 224         case F_DUMMY3OUT:
 225           ak = ia[4];
 226           b1 = ip[tp].dummy.b;
 227           c1 = ip[tp].dummy.c;
 228           constr_dum3OUT(x[ai],x[aj],x[ak],x[adum],a1,b1,c1);
 229           break;
 230         case F_DUMMY4FD:
 231           ak = ia[4];
 232           al = ia[5];
 233           b1 = ip[tp].dummy.b;
 234           c1 = ip[tp].dummy.c;
 235           constr_dum4FD(x[ai],x[aj],x[ak],x[al],x[adum],a1,b1,c1);
 236           break;
 237         default:
 238           fatal_error(0,"No such dummy type %d in %s, line %d",
 239                       ftype,__FILE__,__LINE__);
 240         }
 241         if (v) {
 242           /* Calculate velocity of dummy... */
 243           rvec_sub(x[adum],xd,vv);
 244           svmul(inv_dt,vv,v[adum]);
 245         }
 246         /* Increment loop variables */
 247         i  += nra+1;
 248         ia += nra+1;
 249       }
 250     }
 251   }
 252 }
 253
 254 static void spread_dum2(rvec fi,rvec fj,rvec f,real a)
 255 {
 256   real fx,fy,fz,b;
 257
 258   b=1.0-a;
 259   /* 1 flop */
 260
 261   fx=f[XX];
 262   fy=f[YY];
 263   fz=f[ZZ];
 264   fi[XX]+=b*fx;
 265   fi[YY]+=b*fy;
 266   fi[ZZ]+=b*fz;
 267   fj[XX]+=a*fx;
 268   fj[YY]+=a*fy;
 269   fj[ZZ]+=a*fz;
 270   /* 6 Flops */
 271
 272   /* TOTAL: 7 flops */
 273 }
 274
 275 static void spread_dum3(rvec fi,rvec fj,rvec fk,rvec f,real a,real b)
 276 {
 277   real fx,fy,fz,c;
 278
 279   c=1.0-a-b;
 280   /* 2 flops */
 281
 282   fx=f[XX];
 283   fy=f[YY];
 284   fz=f[ZZ];
 285   fi[XX]+=c*fx;
 286   fi[YY]+=c*fy;
 287   fi[ZZ]+=c*fz;
 288   fj[XX]+=a*fx;
 289   fj[YY]+=a*fy;
 290   fj[ZZ]+=a*fz;
 291   fk[XX]+=b*fx;
 292   fk[YY]+=b*fy;
 293   fk[ZZ]+=b*fz;
 294   /* 9 Flops */
 295
 296   /* TOTAL: 11 flops */
 297 }
 298
 299 static void spread_dum3FD(rvec xi,rvec xj,rvec xk,
 300                           rvec fi,rvec fj,rvec fk,rvec f,real a,real b)
 301 {
 302   real fx,fy,fz,c,invl,fproj,a1;
 303   rvec xij,xjk,xix,temp;
 304
 305   rvec_sub(xj,xi,xij);
 306   rvec_sub(xk,xj,xjk);
 307   /* 6 flops */
 308
 309   /* xix goes from i to point x on the line jk */
 310   xix[XX]=xij[XX]+a*xjk[XX];
 311   xix[YY]=xij[YY]+a*xjk[YY];
 312   xix[ZZ]=xij[ZZ]+a*xjk[ZZ];
 313   /* 6 flops */
 314
 315   invl=invsqrt(iprod(xix,xix));
 316   c=b*invl;
 317   /* 4 + ?10? flops */
 318
 319   fproj=iprod(xix,f)*invl*invl; /* = (xix . f)/(xix . xix) */
 320
 321   fx=f[XX];
 322   fy=f[YY];
 323   fz=f[ZZ];
 324
 325   temp[XX]=c*(fx-fproj*xix[XX]);
 326   temp[YY]=c*(fy-fproj*xix[YY]);
 327   temp[ZZ]=c*(fz-fproj*xix[ZZ]);
 328   /* 16 */
 329
 330   /* c is already calculated in constr_dum3FD
 331      storing c somewhere will save 26 flops!     */
 332
 333   a1=1-a;
 334   fi[XX]+=fx-temp[XX];
 335   fi[YY]+=fy-temp[YY];
 336   fi[ZZ]+=fz-temp[ZZ];
 337   fj[XX]+=a1*temp[XX];
 338   fj[YY]+=a1*temp[YY];
 339   fj[ZZ]+=a1*temp[ZZ];
 340   fk[XX]+= a*temp[XX];
 341   fk[YY]+= a*temp[YY];
 342   fk[ZZ]+= a*temp[ZZ];
 343   /* 19 Flops */
 344
 345   /* TOTAL: 61 flops */
 346 }
 347
 348 static void spread_dum3FAD(rvec xi,rvec xj,rvec xk,
 349                            rvec fi,rvec fj,rvec fk,rvec f,real a,real b)
 350 {
 351   rvec xij,xjk,xperp,Fpij,Fppp,f1,f2,f3;
 352   real a1,b1,c1,c2,invdij,invdij2,invdp,fproj;
 353   int d;
 354
 355   rvec_sub(xj,xi,xij);
 356   rvec_sub(xk,xj,xjk);
 357   /* 6 flops */
 358
 359   invdij = invsqrt(iprod(xij,xij));
 360   invdij2 = invdij * invdij;
 361   c1 = iprod(xij,xjk) * invdij2;
 362   xperp[XX] = xjk[XX] - c1*xij[XX];
 363   xperp[YY] = xjk[YY] - c1*xij[YY];
 364   xperp[ZZ] = xjk[ZZ] - c1*xij[ZZ];
 365   /* xperp in plane ijk, perp. to ij */
 366   invdp = invsqrt(iprod(xperp,xperp));
 367   a1 = a*invdij;
 368   b1 = b*invdp;
 369   /* 45 flops */
 370
 371   /* a1, b1 and c1 are already calculated in constr_dum3FAD
 372      storing them somewhere will save 45 flops!     */
 373
 374   fproj=iprod(xij  ,f)*invdij2;
 375   svmul(fproj,                     xij,  Fpij); /* proj. f on xij */
 376   svmul(iprod(xperp,f)*invdp*invdp,xperp,Fppp); /* proj. f on xperp */
 377   svmul(b1*fproj,                  xperp,f3);
 378   /* 23 flops */
 379
 380   rvec_sub(f,Fpij,f1);  /* f1 = f - Fpij */
 381   rvec_sub(f1,Fppp,f2); /* f2 = f - Fpij - Fppp */
 382   for (d=0; (d<DIM); d++) {
 383     f1[d]*=a1;
 384     f2[d]*=b1;
 385   }
 386   /* 12 flops */
 387
 388   c2=1+c1;
 389   fi[XX] += f[XX] - f1[XX] + c1*f2[XX] + f3[XX];
 390   fi[YY] += f[YY] - f1[YY] + c1*f2[YY] + f3[YY];
 391   fi[ZZ] += f[ZZ] - f1[ZZ] + c1*f2[ZZ] + f3[ZZ];
 392   fj[XX] +=         f1[XX] - c2*f2[XX] - f3[XX];
 393   fj[YY] +=         f1[YY] - c2*f2[YY] - f3[YY];
 394   fj[ZZ] +=         f1[ZZ] - c2*f2[ZZ] - f3[ZZ];
 395   fk[XX] +=                     f2[XX];
 396   fk[YY] +=                     f2[YY];
 397   fk[ZZ] +=                     f2[ZZ];
 398   /* 30 Flops */
 399
 400   /* TOTAL: 113 flops */
 401 }
 402
 403 static void spread_dum3OUT(rvec xi,rvec xj,rvec xk,
 404                         rvec fi,rvec fj,rvec fk,rvec f,real a,real b,real c)
 405 {
 406   rvec xij,xik,ffj,ffk;
 407   real cfx,cfy,cfz;
 408   int  m;
 409
 410   rvec_sub(xj,xi,xij);
 411   rvec_sub(xk,xi,xik);
 412   /* 6 Flops */
 413
 414   cfx=c*f[XX];
 415   cfy=c*f[YY];
 416   cfz=c*f[ZZ];
 417   /* 3 Flops */
 418
 419   ffj[XX] =  a*f[XX]     - xik[ZZ]*cfy + xik[YY]*cfz;
 420   ffj[YY] =  xik[ZZ]*cfx + a*f[YY]     - xik[XX]*cfz;
 421   ffj[ZZ] = -xik[YY]*cfx + xik[XX]*cfy + a*f[ZZ];
 422
 423   ffk[XX] =  b*f[XX]     + xij[ZZ]*cfy - xij[YY]*cfz;
 424   ffk[YY] = -xij[ZZ]*cfx + b*f[YY]     + xij[XX]*cfz;
 425   ffk[ZZ] =  xij[YY]*cfx - xij[XX]*cfy + b*f[ZZ];
 426   /* 30 Flops */
 427
 428   for(m=0; (m<DIM); m++) {
 429     fi[m]+=f[m]-ffj[m]-ffk[m];
 430     fj[m]+=ffj[m];
 431     fk[m]+=ffk[m];
 432   }
 433   /* 15 Flops */
 434
 435   /* TOTAL: 54 flops */
 436 }
 437
 438 static void spread_dum4FD(rvec xi,rvec xj,rvec xk,rvec xl,
 439                           rvec fi,rvec fj,rvec fk,rvec fl,rvec f,
 440                           real a,real b,real c)
 441 {
 442   real fx,fy,fz,d,invl,fproj,a1;
 443   rvec xij,xjk,xjl,xix,temp;
 444
 445   rvec_sub(xj,xi,xij);
 446   rvec_sub(xk,xj,xjk);
 447   rvec_sub(xl,xj,xjl);
 448   /* 9 flops */
 449
 450   /* xix goes from i to point x on the plane jkl */
 451   xix[XX]=xij[XX]+a*xjk[XX]+b*xjl[XX];
 452   xix[YY]=xij[YY]+a*xjk[YY]+b*xjl[YY];
 453   xix[ZZ]=xij[ZZ]+a*xjk[ZZ]+b*xjl[ZZ];
 454   /* 12 flops */
 455
 456   invl=invsqrt(iprod(xix,xix));
 457   d=c*invl;
 458   /* 4 + ?10? flops */
 459
 460   fproj=iprod(xix,f)*invl*invl; /* = (xix . f)/(xix . xix) */
 461
 462   fx=f[XX];
 463   fy=f[YY];
 464   fz=f[ZZ];
 465
 466   temp[XX]=d*(fx-fproj*xix[XX]);
 467   temp[YY]=d*(fy-fproj*xix[YY]);
 468   temp[ZZ]=d*(fz-fproj*xix[ZZ]);
 469   /* 16 */
 470
 471   /* c is already calculated in constr_dum3FD
 472      storing c somewhere will save 35 flops!     */
 473
 474   a1=1-a-b;
 475   fi[XX]+=fx-temp[XX];
 476   fi[YY]+=fy-temp[YY];
 477   fi[ZZ]+=fz-temp[ZZ];
 478   fj[XX]+=a1*temp[XX];
 479   fj[YY]+=a1*temp[YY];
 480   fj[ZZ]+=a1*temp[ZZ];
 481   fk[XX]+= a*temp[XX];
 482   fk[YY]+= a*temp[YY];
 483   fk[ZZ]+= a*temp[ZZ];
 484   fl[XX]+= b*temp[XX];
 485   fl[YY]+= b*temp[YY];
 486   fl[ZZ]+= b*temp[ZZ];
 487   /* 26 Flops */
 488
 489   /* TOTAL: 77 flops */
 490 }
 491
 492 void spread_dummy_f(FILE *log,rvec x[],rvec f[],t_nrnb *nrnb,t_idef *idef)
 493 {
 494   real      a1,b1,c1;
 495   int       i,nra,nrd,tp,ftype;
 496   int       nd2,nd3,nd3FD,nd3FAD,nd3OUT,nd4FD;
 497   t_iatom   adum,ai,aj,ak,al;
 498   t_iatom   *ia;
 499   t_iparams *ip;
 500
 501   ip     = idef->iparams;
 502
 503   nd2    = 0;
 504   nd3    = 0;
 505   nd3FD  = 0;
 506   nd3FAD = 0;
 507   nd3OUT = 0;
 508   nd4FD  = 0;
 509   /* this loop goes backwards to be able to build *
 510    * higher type dummies from lower types         */
 511   for(ftype=F_NRE-1; (ftype>=0); ftype--) {
 512     if (interaction_function[ftype].flags & IF_DUMMY) {
 513       nra    = interaction_function[ftype].nratoms;
 514       nrd    = idef->il[ftype].nr;
 515       ia     = idef->il[ftype].iatoms;
 516
 517       for(i=0; (i<nrd); ) {
 518         tp   = ia[0];
 519         assert(ftype == idef->functype[tp]);
 520
 521         /* The dummy and constructing atoms */
 522         adum = ia[1];
 523         ai   = ia[2];
 524         aj   = ia[3];
 525
 526         /* Constants for constructing */
 527         a1   = ip[tp].dummy.a;
 528
 529         /* Construct the dummy depending on type */
 530         switch (ftype) {
 531         case F_DUMMY2:
 532           spread_dum2(f[ai],f[aj],f[adum],a1);
 533           nd2++;
 534           break;
 535         case F_DUMMY3:
 536           ak = ia[4];
 537           b1 = ip[tp].dummy.b;
 538           spread_dum3(f[ai],f[aj],f[ak],f[adum],a1,b1);
 539           nd3++;
 540           break;
 541         case F_DUMMY3FD:
 542           ak = ia[4];
 543           b1 = ip[tp].dummy.b;
 544           spread_dum3FD(x[ai],x[aj],x[ak],f[ai],f[aj],f[ak],f[adum],a1,b1);
 545           nd3FD++;
 546           break;
 547         case F_DUMMY3FAD:
 548           ak = ia[4];
 549           b1 = ip[tp].dummy.b;
 550           spread_dum3FAD(x[ai],x[aj],x[ak],f[ai],f[aj],f[ak],f[adum],a1,b1);
 551           nd3FAD++;
 552           break;
 553         case F_DUMMY3OUT:
 554           ak = ia[4];
 555           b1 = ip[tp].dummy.b;
 556           c1 = ip[tp].dummy.c;
 557           spread_dum3OUT(x[ai],x[aj],x[ak],f[ai],f[aj],f[ak],f[adum],a1,b1,c1);
 558           nd3OUT++;
 559           break;
 560         case F_DUMMY4FD:
 561           ak = ia[4];
 562           al = ia[5];
 563           b1 = ip[tp].dummy.b;
 564           c1 = ip[tp].dummy.c;
 565           spread_dum4FD(x[ai],x[aj],x[ak],x[al],
 566                         f[ai],f[aj],f[ak],f[al],f[adum],a1,b1,c1);
 567           nd4FD++;
 568           break;
 569         default:
 570           fatal_error(0,"No such dummy type %d in %s, line %d",
 571                       ftype,__FILE__,__LINE__);
 572         }
 573         clear_rvec(f[adum]);
 574
 575         /* Increment loop variables */
 576         i  += nra+1;
 577         ia += nra+1;
 578       }
 579     }
 580   }
 581   inc_nrnb(nrnb,eNR_DUM2,   nd2   );
 582   inc_nrnb(nrnb,eNR_DUM3,   nd3   );
 583   inc_nrnb(nrnb,eNR_DUM3FD, nd3FD );
 584   inc_nrnb(nrnb,eNR_DUM3FAD,nd3FAD);
 585   inc_nrnb(nrnb,eNR_DUM3OUT,nd3OUT);
 586   inc_nrnb(nrnb,eNR_DUM4FD, nd4FD );
 587 }
 588