src/gromacs/gmxlib/bondfree.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  * GROningen Mixture of Alchemy and Childrens' Stories
  35  */
  36 #ifdef HAVE_CONFIG_H
  37 #include <config.h>
  38 #endif
  39
  40 #include <math.h>
  41 #include "physics.h"
  42 #include "vec.h"
  43 #include "maths.h"
  44 #include "txtdump.h"
  45 #include "bondf.h"
  46 #include "smalloc.h"
  47 #include "pbc.h"
  48 #include "ns.h"
  49 #include "macros.h"
  50 #include "names.h"
  51 #include "gmx_fatal.h"
  52 #include "mshift.h"
  53 #include "main.h"
  54 #include "disre.h"
  55 #include "orires.h"
  56 #include "force.h"
  57 #include "nonbonded.h"
  58 #include "mdrun.h"
  59
  60 /* Find a better place for this? */
  61 const int cmap_coeff_matrix[] = {
  62 1, 0, -3,  2, 0, 0,  0,  0, -3,  0,  9, -6,  2,  0, -6,  4 ,
  63 0, 0,  0,  0, 0, 0,  0,  0,  3,  0, -9,  6, -2,  0,  6, -4,
  64 0, 0,  0,  0, 0, 0,  0,  0,  0,  0,  9, -6,  0,  0, -6,  4 ,
  65 0, 0,  3, -2, 0, 0,  0,  0,  0,  0, -9,  6,  0,  0,  6, -4,
  66 0, 0,  0,  0, 1, 0, -3,  2, -2,  0,  6, -4,  1,  0, -3,  2 ,
  67 0, 0,  0,  0, 0, 0,  0,  0, -1,  0,  3, -2,  1,  0, -3,  2 ,
  68 0, 0,  0,  0, 0, 0,  0,  0,  0,  0, -3,  2,  0,  0,  3, -2,
  69 0, 0,  0,  0, 0, 0,  3, -2,  0,  0, -6,  4,  0,  0,  3, -2,
  70 0, 1, -2,  1, 0, 0,  0,  0,  0, -3,  6, -3,  0,  2, -4,  2 ,
  71 0, 0,  0,  0, 0, 0,  0,  0,  0,  3, -6,  3,  0, -2,  4, -2,
  72 0, 0,  0,  0, 0, 0,  0,  0,  0,  0, -3,  3,  0,  0,  2, -2,
  73 0, 0, -1,  1, 0, 0,  0,  0,  0,  0,  3, -3,  0,  0, -2,  2 ,
  74 0, 0,  0,  0, 0, 1, -2,  1,  0, -2,  4, -2,  0,  1, -2,  1,
  75 0, 0,  0,  0, 0, 0,  0,  0,  0, -1,  2, -1,  0,  1, -2,  1,
  76 0, 0,  0,  0, 0, 0,  0,  0,  0,  0,  1, -1,  0,  0, -1,  1,
  77 0, 0,  0,  0, 0, 0, -1,  1,  0,  0,  2, -2,  0,  0, -1,  1
  78 };
  79
  80
  81
  82 int glatnr(int *global_atom_index,int i)
  83 {
  84     int atnr;
  85
  86     if (global_atom_index == NULL) {
  87         atnr = i + 1;
  88     } else {
  89         atnr = global_atom_index[i] + 1;
  90     }
  91
  92     return atnr;
  93 }
  94
  95 static int pbc_rvec_sub(const t_pbc *pbc,const rvec xi,const rvec xj,rvec dx)
  96 {
  97   if (pbc) {
  98     return pbc_dx_aiuc(pbc,xi,xj,dx);
  99   }
 100   else {
 101     rvec_sub(xi,xj,dx);
 102     return CENTRAL;
 103   }
 104 }
 105
 106 /*
 107  * Morse potential bond by Frank Everdij
 108  *
 109  * Three parameters needed:
 110  *
 111  * b0 = equilibrium distance in nm
 112  * be = beta in nm^-1 (actually, it's nu_e*Sqrt(2*pi*pi*mu/D_e))
 113  * cb = well depth in kJ/mol
 114  *
 115  * Note: the potential is referenced to be +cb at infinite separation
 116  *       and zero at the equilibrium distance!
 117  */
 118
 119 real morse_bonds(int nbonds,
 120                  const t_iatom forceatoms[],const t_iparams forceparams[],
 121                  const rvec x[],rvec f[],rvec fshift[],
 122                  const t_pbc *pbc,const t_graph *g,
 123                  real lambda,real *dvdlambda,
 124                  const t_mdatoms *md,t_fcdata *fcd,
 125                  int *global_atom_index)
 126 {
 127   const real one=1.0;
 128   const real two=2.0;
 129   real  dr,dr2,temp,omtemp,cbomtemp,fbond,vbond,fij,vtot;
 130   real  b0,be,cb,b0A,beA,cbA,b0B,beB,cbB,L1;
 131   rvec  dx;
 132   int   i,m,ki,type,ai,aj;
 133   ivec  dt;
 134
 135   vtot = 0.0;
 136   for(i=0; (i<nbonds); ) {
 137     type = forceatoms[i++];
 138     ai   = forceatoms[i++];
 139     aj   = forceatoms[i++];
 140
 141     b0A   = forceparams[type].morse.b0A;
 142     beA   = forceparams[type].morse.betaA;
 143     cbA   = forceparams[type].morse.cbA;
 144
 145     b0B   = forceparams[type].morse.b0B;
 146     beB   = forceparams[type].morse.betaB;
 147     cbB   = forceparams[type].morse.cbB;
 148
 149     L1 = one-lambda;                      /* 1 */
 150     b0 = L1*b0A + lambda*b0B;             /* 3 */
 151     be = L1*beA + lambda*beB;             /* 3 */
 152     cb = L1*cbA + lambda*cbB;             /* 3 */
 153
 154     ki   = pbc_rvec_sub(pbc,x[ai],x[aj],dx);            /*   3          */
 155     dr2  = iprod(dx,dx);                            /*   5          */
 156     dr   = dr2*gmx_invsqrt(dr2);                        /*  10          */
 157     temp = exp(-be*(dr-b0));                        /*  12          */
 158
 159     if (temp == one)
 160     {
 161         /* bonds are constrainted. This may _not_ include bond constraints if they are lambda dependent */
 162         *dvdlambda = cbB-cbA;
 163         continue;
 164     }
 165
 166     omtemp   = one-temp;                               /*   1          */
 167     cbomtemp = cb*omtemp;                              /*   1          */
 168     vbond    = cbomtemp*omtemp;                        /*   1          */
 169     fbond    = -two*be*temp*cbomtemp*gmx_invsqrt(dr2); /*   9          */
 170     vtot     += vbond;                                 /*   1          */
 171
 172     *dvdlambda += (cbB - cbA) * omtemp * omtemp - (2-2*omtemp)*omtemp * cb * ((b0B-b0A)*be - (beB-beA)*(dr-b0)); /* 15 */
 173
 174     if (g) {
 175       ivec_sub(SHIFT_IVEC(g,ai),SHIFT_IVEC(g,aj),dt);
 176       ki = IVEC2IS(dt);
 177     }
 178
 179     for (m=0; (m<DIM); m++) {                          /*  15          */
 180       fij=fbond*dx[m];
 181       f[ai][m]+=fij;
 182       f[aj][m]-=fij;
 183       fshift[ki][m]+=fij;
 184       fshift[CENTRAL][m]-=fij;
 185     }
 186   }                                           /*  83 TOTAL    */
 187   return vtot;
 188 }
 189
 190 real cubic_bonds(int nbonds,
 191                  const t_iatom forceatoms[],const t_iparams forceparams[],
 192                  const rvec x[],rvec f[],rvec fshift[],
 193                  const t_pbc *pbc,const t_graph *g,
 194                  real lambda,real *dvdlambda,
 195                  const t_mdatoms *md,t_fcdata *fcd,
 196                  int *global_atom_index)
 197 {
 198   const real three = 3.0;
 199   const real two   = 2.0;
 200   real  kb,b0,kcub;
 201   real  dr,dr2,dist,kdist,kdist2,fbond,vbond,fij,vtot;
 202   rvec  dx;
 203   int   i,m,ki,type,ai,aj;
 204   ivec  dt;
 205
 206   vtot = 0.0;
 207   for(i=0; (i<nbonds); ) {
 208     type = forceatoms[i++];
 209     ai   = forceatoms[i++];
 210     aj   = forceatoms[i++];
 211
 212     b0   = forceparams[type].cubic.b0;
 213     kb   = forceparams[type].cubic.kb;
 214     kcub = forceparams[type].cubic.kcub;
 215
 216     ki   = pbc_rvec_sub(pbc,x[ai],x[aj],dx);                /*   3          */
 217     dr2  = iprod(dx,dx);                                /*   5          */
 218
 219     if (dr2 == 0.0)
 220       continue;
 221
 222     dr         = dr2*gmx_invsqrt(dr2);                      /*  10          */
 223     dist       = dr-b0;
 224     kdist      = kb*dist;
 225     kdist2     = kdist*dist;
 226
 227     vbond      = kdist2 + kcub*kdist2*dist;
 228     fbond      = -(two*kdist + three*kdist2*kcub)/dr;
 229
 230     vtot      += vbond;       /* 21 */
 231
 232     if (g) {
 233       ivec_sub(SHIFT_IVEC(g,ai),SHIFT_IVEC(g,aj),dt);
 234       ki=IVEC2IS(dt);
 235     }
 236     for (m=0; (m<DIM); m++) {                          /*  15          */
 237       fij=fbond*dx[m];
 238       f[ai][m]+=fij;
 239       f[aj][m]-=fij;
 240       fshift[ki][m]+=fij;
 241       fshift[CENTRAL][m]-=fij;
 242     }
 243   }                                           /*  54 TOTAL    */
 244   return vtot;
 245 }
 246
 247 real FENE_bonds(int nbonds,
 248                 const t_iatom forceatoms[],const t_iparams forceparams[],
 249                 const rvec x[],rvec f[],rvec fshift[],
 250                 const t_pbc *pbc,const t_graph *g,
 251                 real lambda,real *dvdlambda,
 252                 const t_mdatoms *md,t_fcdata *fcd,
 253                 int *global_atom_index)
 254 {
 255   const real half=0.5;
 256   const real one=1.0;
 257   real  bm,kb;
 258   real  dr,dr2,bm2,omdr2obm2,fbond,vbond,fij,vtot;
 259   rvec  dx;
 260   int   i,m,ki,type,ai,aj;
 261   ivec  dt;
 262
 263   vtot = 0.0;
 264   for(i=0; (i<nbonds); ) {
 265     type = forceatoms[i++];
 266     ai   = forceatoms[i++];
 267     aj   = forceatoms[i++];
 268
 269     bm   = forceparams[type].fene.bm;
 270     kb   = forceparams[type].fene.kb;
 271
 272     ki   = pbc_rvec_sub(pbc,x[ai],x[aj],dx);            /*   3          */
 273     dr2  = iprod(dx,dx);                                /*   5          */
 274
 275     if (dr2 == 0.0)
 276       continue;
 277
 278     bm2 = bm*bm;
 279
 280     if (dr2 >= bm2)
 281       gmx_fatal(FARGS,
 282                 "r^2 (%f) >= bm^2 (%f) in FENE bond between atoms %d and %d",
 283                 dr2,bm2,
 284                 glatnr(global_atom_index,ai),
 285                 glatnr(global_atom_index,aj));
 286
 287     omdr2obm2  = one - dr2/bm2;
 288
 289     vbond      = -half*kb*bm2*log(omdr2obm2);
 290     fbond      = -kb/omdr2obm2;
 291
 292     vtot      += vbond;       /* 35 */
 293
 294     if (g) {
 295       ivec_sub(SHIFT_IVEC(g,ai),SHIFT_IVEC(g,aj),dt);
 296       ki=IVEC2IS(dt);
 297     }
 298     for (m=0; (m<DIM); m++) {                          /*  15          */
 299       fij=fbond*dx[m];
 300       f[ai][m]+=fij;
 301       f[aj][m]-=fij;
 302       fshift[ki][m]+=fij;
 303       fshift[CENTRAL][m]-=fij;
 304     }
 305   }                                           /*  58 TOTAL    */
 306   return vtot;
 307 }
 308
 309 real harmonic(real kA,real kB,real xA,real xB,real x,real lambda,
 310               real *V,real *F)
 311 {
 312   const real half=0.5;
 313   real  L1,kk,x0,dx,dx2;
 314   real  v,f,dvdlambda;
 315
 316   L1    = 1.0-lambda;
 317   kk    = L1*kA+lambda*kB;
 318   x0    = L1*xA+lambda*xB;
 319
 320   dx    = x-x0;
 321   dx2   = dx*dx;
 322
 323   f     = -kk*dx;
 324   v     = half*kk*dx2;
 325   dvdlambda  = half*(kB-kA)*dx2 + (xA-xB)*kk*dx;
 326
 327   *F    = f;
 328   *V    = v;
 329
 330   return dvdlambda;
 331
 332   /* That was 19 flops */
 333 }
 334
 335
 336 real bonds(int nbonds,
 337            const t_iatom forceatoms[],const t_iparams forceparams[],
 338            const rvec x[],rvec f[],rvec fshift[],
 339            const t_pbc *pbc,const t_graph *g,
 340            real lambda,real *dvdlambda,
 341            const t_mdatoms *md,t_fcdata *fcd,
 342            int *global_atom_index)
 343 {
 344   int  i,m,ki,ai,aj,type;
 345   real dr,dr2,fbond,vbond,fij,vtot;
 346   rvec dx;
 347   ivec dt;
 348
 349   vtot = 0.0;
 350   for(i=0; (i<nbonds); ) {
 351     type = forceatoms[i++];
 352     ai   = forceatoms[i++];
 353     aj   = forceatoms[i++];
 354
 355     ki   = pbc_rvec_sub(pbc,x[ai],x[aj],dx);    /*   3          */
 356     dr2  = iprod(dx,dx);                        /*   5          */
 357     dr   = dr2*gmx_invsqrt(dr2);                        /*  10          */
 358
 359     *dvdlambda += harmonic(forceparams[type].harmonic.krA,
 360                            forceparams[type].harmonic.krB,
 361                            forceparams[type].harmonic.rA,
 362                            forceparams[type].harmonic.rB,
 363                            dr,lambda,&vbond,&fbond);  /*  19  */
 364
 365     if (dr2 == 0.0)
 366       continue;
 367
 368
 369     vtot  += vbond;/* 1*/
 370     fbond *= gmx_invsqrt(dr2);                  /*   6          */
 371 #ifdef DEBUG
 372     if (debug)
 373       fprintf(debug,"BONDS: dr = %10g  vbond = %10g  fbond = %10g\n",
 374               dr,vbond,fbond);
 375 #endif
 376     if (g) {
 377       ivec_sub(SHIFT_IVEC(g,ai),SHIFT_IVEC(g,aj),dt);
 378       ki=IVEC2IS(dt);
 379     }
 380     for (m=0; (m<DIM); m++) {                   /*  15          */
 381       fij=fbond*dx[m];
 382       f[ai][m]+=fij;
 383       f[aj][m]-=fij;
 384       fshift[ki][m]+=fij;
 385       fshift[CENTRAL][m]-=fij;
 386     }
 387   }                                     /* 59 TOTAL     */
 388   return vtot;
 389 }
 390
 391 real restraint_bonds(int nbonds,
 392                      const t_iatom forceatoms[],const t_iparams forceparams[],
 393                      const rvec x[],rvec f[],rvec fshift[],
 394                      const t_pbc *pbc,const t_graph *g,
 395                      real lambda,real *dvdlambda,
 396                      const t_mdatoms *md,t_fcdata *fcd,
 397                      int *global_atom_index)
 398 {
 399     int  i,m,ki,ai,aj,type;
 400     real dr,dr2,fbond,vbond,fij,vtot;
 401     real L1;
 402     real low,dlow,up1,dup1,up2,dup2,k,dk;
 403     real drh,drh2;
 404     rvec dx;
 405     ivec dt;
 406
 407     L1   = 1.0 - lambda;
 408
 409     vtot = 0.0;
 410     for(i=0; (i<nbonds); )
 411     {
 412         type = forceatoms[i++];
 413         ai   = forceatoms[i++];
 414         aj   = forceatoms[i++];
 415
 416         ki   = pbc_rvec_sub(pbc,x[ai],x[aj],dx);        /*   3          */
 417         dr2  = iprod(dx,dx);                            /*   5          */
 418         dr   = dr2*gmx_invsqrt(dr2);                    /*  10          */
 419
 420         low  = L1*forceparams[type].restraint.lowA + lambda*forceparams[type].restraint.lowB;
 421         dlow =   -forceparams[type].restraint.lowA +        forceparams[type].restraint.lowB;
 422         up1  = L1*forceparams[type].restraint.up1A + lambda*forceparams[type].restraint.up1B;
 423         dup1 =   -forceparams[type].restraint.up1A +        forceparams[type].restraint.up1B;
 424         up2  = L1*forceparams[type].restraint.up2A + lambda*forceparams[type].restraint.up2B;
 425         dup2 =   -forceparams[type].restraint.up2A +        forceparams[type].restraint.up2B;
 426         k    = L1*forceparams[type].restraint.kA   + lambda*forceparams[type].restraint.kB;
 427         dk   =   -forceparams[type].restraint.kA   +        forceparams[type].restraint.kB;
 428         /* 24 */
 429
 430         if (dr < low)
 431         {
 432             drh   = dr - low;
 433             drh2  = drh*drh;
 434             vbond = 0.5*k*drh2;
 435             fbond = -k*drh;
 436             *dvdlambda += 0.5*dk*drh2 - k*dlow*drh;
 437         } /* 11 */
 438         else if (dr <= up1)
 439         {
 440             vbond = 0;
 441             fbond = 0;
 442         }
 443         else if (dr <= up2)
 444         {
 445             drh   = dr - up1;
 446             drh2  = drh*drh;
 447             vbond = 0.5*k*drh2;
 448             fbond = -k*drh;
 449             *dvdlambda += 0.5*dk*drh2 - k*dup1*drh;
 450         } /* 11 */
 451         else
 452         {
 453             drh   = dr - up2;
 454             vbond = k*(up2 - up1)*(0.5*(up2 - up1) + drh);
 455             fbond = -k*(up2 - up1);
 456             *dvdlambda += dk*(up2 - up1)*(0.5*(up2 - up1) + drh)
 457                 + k*(dup2 - dup1)*(up2 - up1 + drh)
 458                 - k*(up2 - up1)*dup2;
 459         }
 460
 461         if (dr2 == 0.0)
 462             continue;
 463
 464         vtot  += vbond;/* 1*/
 465         fbond *= gmx_invsqrt(dr2);                      /*   6          */
 466 #ifdef DEBUG
 467         if (debug)
 468             fprintf(debug,"BONDS: dr = %10g  vbond = %10g  fbond = %10g\n",
 469                     dr,vbond,fbond);
 470 #endif
 471         if (g) {
 472             ivec_sub(SHIFT_IVEC(g,ai),SHIFT_IVEC(g,aj),dt);
 473             ki=IVEC2IS(dt);
 474         }
 475         for (m=0; (m<DIM); m++) {                       /*  15          */
 476             fij=fbond*dx[m];
 477             f[ai][m]+=fij;
 478             f[aj][m]-=fij;
 479             fshift[ki][m]+=fij;
 480             fshift[CENTRAL][m]-=fij;
 481         }
 482     }                                   /* 59 TOTAL     */
 483
 484     return vtot;
 485 }
 486
 487 real polarize(int nbonds,
 488               const t_iatom forceatoms[],const t_iparams forceparams[],
 489               const rvec x[],rvec f[],rvec fshift[],
 490               const t_pbc *pbc,const t_graph *g,
 491               real lambda,real *dvdlambda,
 492               const t_mdatoms *md,t_fcdata *fcd,
 493               int *global_atom_index)
 494 {
 495   int  i,m,ki,ai,aj,type;
 496   real dr,dr2,fbond,vbond,fij,vtot,ksh;
 497   rvec dx;
 498   ivec dt;
 499
 500   vtot = 0.0;
 501   for(i=0; (i<nbonds); ) {
 502     type = forceatoms[i++];
 503     ai   = forceatoms[i++];
 504     aj   = forceatoms[i++];
 505     ksh  = sqr(md->chargeA[aj])*ONE_4PI_EPS0/forceparams[type].polarize.alpha;
 506     if (debug)
 507       fprintf(debug,"POL: local ai = %d aj = %d ksh = %.3f\n",ai,aj,ksh);
 508
 509     ki   = pbc_rvec_sub(pbc,x[ai],x[aj],dx);    /*   3          */
 510     dr2  = iprod(dx,dx);                        /*   5          */
 511     dr   = dr2*gmx_invsqrt(dr2);                        /*  10          */
 512
 513     *dvdlambda += harmonic(ksh,ksh,0,0,dr,lambda,&vbond,&fbond);  /*  19  */
 514
 515     if (dr2 == 0.0)
 516       continue;
 517
 518     vtot  += vbond;/* 1*/
 519     fbond *= gmx_invsqrt(dr2);                  /*   6          */
 520
 521     if (g) {
 522       ivec_sub(SHIFT_IVEC(g,ai),SHIFT_IVEC(g,aj),dt);
 523       ki=IVEC2IS(dt);
 524     }
 525     for (m=0; (m<DIM); m++) {                   /*  15          */
 526       fij=fbond*dx[m];
 527       f[ai][m]+=fij;
 528       f[aj][m]-=fij;
 529       fshift[ki][m]+=fij;
 530       fshift[CENTRAL][m]-=fij;
 531     }
 532   }                                     /* 59 TOTAL     */
 533   return vtot;
 534 }
 535
 536 real anharm_polarize(int nbonds,
 537                      const t_iatom forceatoms[],const t_iparams forceparams[],
 538                      const rvec x[],rvec f[],rvec fshift[],
 539                      const t_pbc *pbc,const t_graph *g,
 540                      real lambda,real *dvdlambda,
 541                      const t_mdatoms *md,t_fcdata *fcd,
 542                      int *global_atom_index)
 543 {
 544   int  i,m,ki,ai,aj,type;
 545   real dr,dr2,fbond,vbond,fij,vtot,ksh,khyp,drcut,ddr,ddr3;
 546   rvec dx;
 547   ivec dt;
 548
 549   vtot = 0.0;
 550   for(i=0; (i<nbonds); ) {
 551     type  = forceatoms[i++];
 552     ai    = forceatoms[i++];
 553     aj    = forceatoms[i++];
 554     ksh   = sqr(md->chargeA[aj])*ONE_4PI_EPS0/forceparams[type].anharm_polarize.alpha; /* 7*/
 555     khyp  = forceparams[type].anharm_polarize.khyp;
 556     drcut = forceparams[type].anharm_polarize.drcut;
 557     if (debug)
 558       fprintf(debug,"POL: local ai = %d aj = %d ksh = %.3f\n",ai,aj,ksh);
 559
 560     ki   = pbc_rvec_sub(pbc,x[ai],x[aj],dx);    /*   3          */
 561     dr2  = iprod(dx,dx);                        /*   5          */
 562     dr   = dr2*gmx_invsqrt(dr2);                        /*  10          */
 563
 564     *dvdlambda += harmonic(ksh,ksh,0,0,dr,lambda,&vbond,&fbond);  /*  19  */
 565
 566     if (dr2 == 0.0)
 567       continue;
 568
 569     if (dr > drcut) {
 570         ddr    = dr-drcut;
 571         ddr3   = ddr*ddr*ddr;
 572         vbond += khyp*ddr*ddr3;
 573         fbond -= 4*khyp*ddr3;
 574     }
 575     fbond *= gmx_invsqrt(dr2);                  /*   6          */
 576     vtot  += vbond;/* 1*/
 577
 578     if (g) {
 579       ivec_sub(SHIFT_IVEC(g,ai),SHIFT_IVEC(g,aj),dt);
 580       ki=IVEC2IS(dt);
 581     }
 582     for (m=0; (m<DIM); m++) {                   /*  15          */
 583       fij=fbond*dx[m];
 584       f[ai][m]+=fij;
 585       f[aj][m]-=fij;
 586       fshift[ki][m]+=fij;
 587       fshift[CENTRAL][m]-=fij;
 588     }
 589   }                                     /* 72 TOTAL     */
 590   return vtot;
 591 }
 592
 593 real water_pol(int nbonds,
 594                const t_iatom forceatoms[],const t_iparams forceparams[],
 595                const rvec x[],rvec f[],rvec fshift[],
 596                const t_pbc *pbc,const t_graph *g,
 597                real lambda,real *dvdlambda,
 598                const t_mdatoms *md,t_fcdata *fcd,
 599                int *global_atom_index)
 600 {
 601   /* This routine implements anisotropic polarizibility for water, through
 602    * a shell connected to a dummy with spring constant that differ in the
 603    * three spatial dimensions in the molecular frame.
 604    */
 605   int  i,m,aO,aH1,aH2,aD,aS,type,type0;
 606   rvec dOH1,dOH2,dHH,dOD,dDS,nW,kk,dx,kdx,proj;
 607 #ifdef DEBUG
 608   rvec df;
 609 #endif
 610   real vtot,fij,r_HH,r_OD,r_nW,tx,ty,tz,qS;
 611
 612   vtot = 0.0;
 613   if (nbonds > 0) {
 614     type0  = forceatoms[0];
 615     aS     = forceatoms[5];
 616     qS     = md->chargeA[aS];
 617     kk[XX] = sqr(qS)*ONE_4PI_EPS0/forceparams[type0].wpol.al_x;
 618     kk[YY] = sqr(qS)*ONE_4PI_EPS0/forceparams[type0].wpol.al_y;
 619     kk[ZZ] = sqr(qS)*ONE_4PI_EPS0/forceparams[type0].wpol.al_z;
 620     r_HH   = 1.0/forceparams[type0].wpol.rHH;
 621     r_OD   = 1.0/forceparams[type0].wpol.rOD;
 622     if (debug) {
 623       fprintf(debug,"WPOL: qS  = %10.5f aS = %5d\n",qS,aS);
 624       fprintf(debug,"WPOL: kk  = %10.3f        %10.3f        %10.3f\n",
 625               kk[XX],kk[YY],kk[ZZ]);
 626       fprintf(debug,"WPOL: rOH = %10.3f  rHH = %10.3f  rOD = %10.3f\n",
 627               forceparams[type0].wpol.rOH,
 628               forceparams[type0].wpol.rHH,
 629               forceparams[type0].wpol.rOD);
 630     }
 631     for(i=0; (i<nbonds); i+=6) {
 632       type = forceatoms[i];
 633       if (type != type0)
 634         gmx_fatal(FARGS,"Sorry, type = %d, type0 = %d, file = %s, line = %d",
 635                     type,type0,__FILE__,__LINE__);
 636       aO   = forceatoms[i+1];
 637       aH1  = forceatoms[i+2];
 638       aH2  = forceatoms[i+3];
 639       aD   = forceatoms[i+4];
 640       aS   = forceatoms[i+5];
 641
 642       /* Compute vectors describing the water frame */
 643       rvec_sub(x[aH1],x[aO], dOH1);
 644       rvec_sub(x[aH2],x[aO], dOH2);
 645       rvec_sub(x[aH2],x[aH1],dHH);
 646       rvec_sub(x[aD], x[aO], dOD);
 647       rvec_sub(x[aS], x[aD], dDS);
 648       cprod(dOH1,dOH2,nW);
 649
 650       /* Compute inverse length of normal vector
 651        * (this one could be precomputed, but I'm too lazy now)
 652        */
 653       r_nW = gmx_invsqrt(iprod(nW,nW));
 654       /* This is for precision, but does not make a big difference,
 655        * it can go later.
 656        */
 657       r_OD = gmx_invsqrt(iprod(dOD,dOD));
 658
 659       /* Normalize the vectors in the water frame */
 660       svmul(r_nW,nW,nW);
 661       svmul(r_HH,dHH,dHH);
 662       svmul(r_OD,dOD,dOD);
 663
 664       /* Compute displacement of shell along components of the vector */
 665       dx[ZZ] = iprod(dDS,dOD);
 666       /* Compute projection on the XY plane: dDS - dx[ZZ]*dOD */
 667       for(m=0; (m<DIM); m++)
 668         proj[m] = dDS[m]-dx[ZZ]*dOD[m];
 669
 670       /*dx[XX] = iprod(dDS,nW);
 671         dx[YY] = iprod(dDS,dHH);*/
 672       dx[XX] = iprod(proj,nW);
 673       for(m=0; (m<DIM); m++)
 674         proj[m] -= dx[XX]*nW[m];
 675       dx[YY] = iprod(proj,dHH);
 676       /*#define DEBUG*/
 677 #ifdef DEBUG
 678       if (debug) {
 679         fprintf(debug,"WPOL: dx2=%10g  dy2=%10g  dz2=%10g  sum=%10g  dDS^2=%10g\n",
 680                 sqr(dx[XX]),sqr(dx[YY]),sqr(dx[ZZ]),iprod(dx,dx),iprod(dDS,dDS));
 681         fprintf(debug,"WPOL: dHH=(%10g,%10g,%10g)\n",dHH[XX],dHH[YY],dHH[ZZ]);
 682         fprintf(debug,"WPOL: dOD=(%10g,%10g,%10g), 1/r_OD = %10g\n",
 683                 dOD[XX],dOD[YY],dOD[ZZ],1/r_OD);
 684         fprintf(debug,"WPOL: nW =(%10g,%10g,%10g), 1/r_nW = %10g\n",
 685                 nW[XX],nW[YY],nW[ZZ],1/r_nW);
 686         fprintf(debug,"WPOL: dx  =%10g, dy  =%10g, dz  =%10g\n",
 687                 dx[XX],dx[YY],dx[ZZ]);
 688         fprintf(debug,"WPOL: dDSx=%10g, dDSy=%10g, dDSz=%10g\n",
 689                 dDS[XX],dDS[YY],dDS[ZZ]);
 690       }
 691 #endif
 692       /* Now compute the forces and energy */
 693       kdx[XX] = kk[XX]*dx[XX];
 694       kdx[YY] = kk[YY]*dx[YY];
 695       kdx[ZZ] = kk[ZZ]*dx[ZZ];
 696       vtot   += iprod(dx,kdx);
 697       for(m=0; (m<DIM); m++) {
 698         /* This is a tensor operation but written out for speed */
 699         tx        =  nW[m]*kdx[XX];
 700         ty        = dHH[m]*kdx[YY];
 701         tz        = dOD[m]*kdx[ZZ];
 702         fij       = -tx-ty-tz;
 703 #ifdef DEBUG
 704         df[m] = fij;
 705 #endif
 706         f[aS][m] += fij;
 707         f[aD][m] -= fij;
 708       }
 709 #ifdef DEBUG
 710       if (debug) {
 711         fprintf(debug,"WPOL: vwpol=%g\n",0.5*iprod(dx,kdx));
 712         fprintf(debug,"WPOL: df = (%10g, %10g, %10g)\n",df[XX],df[YY],df[ZZ]);
 713       }
 714 #endif
 715     }
 716   }
 717   return 0.5*vtot;
 718 }
 719
 720 static real do_1_thole(const rvec xi,const rvec xj,rvec fi,rvec fj,
 721                        const t_pbc *pbc,real qq,
 722                        rvec fshift[],real afac)
 723 {
 724   rvec r12;
 725   real r12sq,r12_1,r12n,r12bar,v0,v1,fscal,ebar,fff;
 726   int  m,t;
 727
 728   t      = pbc_rvec_sub(pbc,xi,xj,r12); /*  3 */
 729
 730   r12sq  = iprod(r12,r12);              /*  5 */
 731   r12_1  = gmx_invsqrt(r12sq);              /*  5 */
 732   r12bar = afac/r12_1;                  /*  5 */
 733   v0     = qq*ONE_4PI_EPS0*r12_1;       /*  2 */
 734   ebar   = exp(-r12bar);                /*  5 */
 735   v1     = (1-(1+0.5*r12bar)*ebar);     /*  4 */
 736   fscal  = ((v0*r12_1)*v1 - v0*0.5*afac*ebar*(r12bar+1))*r12_1; /* 9 */
 737   if (debug)
 738     fprintf(debug,"THOLE: v0 = %.3f v1 = %.3f r12= % .3f r12bar = %.3f fscal = %.3f  ebar = %.3f\n",v0,v1,1/r12_1,r12bar,fscal,ebar);
 739
 740   for(m=0; (m<DIM); m++) {
 741     fff    = fscal*r12[m];
 742     fi[m] += fff;
 743     fj[m] -= fff;
 744     fshift[t][m]       += fff;
 745     fshift[CENTRAL][m] -= fff;
 746   } /* 15 */
 747
 748   return v0*v1; /* 1 */
 749   /* 54 */
 750 }
 751
 752 real thole_pol(int nbonds,
 753                const t_iatom forceatoms[],const t_iparams forceparams[],
 754                const rvec x[],rvec f[],rvec fshift[],
 755                const t_pbc *pbc,const t_graph *g,
 756                real lambda,real *dvdlambda,
 757                const t_mdatoms *md,t_fcdata *fcd,
 758                int *global_atom_index)
 759 {
 760   /* Interaction between two pairs of particles with opposite charge */
 761   int i,type,a1,da1,a2,da2;
 762   real q1,q2,qq,a,al1,al2,afac;
 763   real V=0;
 764
 765   for(i=0; (i<nbonds); ) {
 766     type  = forceatoms[i++];
 767     a1    = forceatoms[i++];
 768     da1   = forceatoms[i++];
 769     a2    = forceatoms[i++];
 770     da2   = forceatoms[i++];
 771     q1    = md->chargeA[da1];
 772     q2    = md->chargeA[da2];
 773     a     = forceparams[type].thole.a;
 774     al1   = forceparams[type].thole.alpha1;
 775     al2   = forceparams[type].thole.alpha2;
 776     qq    = q1*q2;
 777     afac  = a*pow(al1*al2,-1.0/6.0);
 778     V += do_1_thole(x[a1], x[a2], f[a1], f[a2], pbc, qq,fshift,afac);
 779     V += do_1_thole(x[da1],x[a2], f[da1],f[a2], pbc,-qq,fshift,afac);
 780     V += do_1_thole(x[a1], x[da2],f[a1], f[da2],pbc,-qq,fshift,afac);
 781     V += do_1_thole(x[da1],x[da2],f[da1],f[da2],pbc, qq,fshift,afac);
 782   }
 783   /* 290 flops */
 784   return V;
 785 }
 786
 787 real bond_angle(const rvec xi,const rvec xj,const rvec xk,const t_pbc *pbc,
 788                 rvec r_ij,rvec r_kj,real *costh,
 789                 int *t1,int *t2)
 790 /* Return value is the angle between the bonds i-j and j-k */
 791 {
 792   /* 41 FLOPS */
 793   real th;
 794
 795   *t1 = pbc_rvec_sub(pbc,xi,xj,r_ij);                   /*  3           */
 796   *t2 = pbc_rvec_sub(pbc,xk,xj,r_kj);                   /*  3           */
 797
 798   *costh=cos_angle(r_ij,r_kj);          /* 25           */
 799   th=acos(*costh);                      /* 10           */
 800                                         /* 41 TOTAL     */
 801   return th;
 802 }
 803
 804 real angles(int nbonds,
 805             const t_iatom forceatoms[],const t_iparams forceparams[],
 806             const rvec x[],rvec f[],rvec fshift[],
 807             const t_pbc *pbc,const t_graph *g,
 808             real lambda,real *dvdlambda,
 809             const t_mdatoms *md,t_fcdata *fcd,
 810             int *global_atom_index)
 811 {
 812   int  i,ai,aj,ak,t1,t2,type;
 813   rvec r_ij,r_kj;
 814   real cos_theta,cos_theta2,theta,dVdt,va,vtot;
 815   ivec jt,dt_ij,dt_kj;
 816
 817   vtot = 0.0;
 818   for(i=0; (i<nbonds); ) {
 819     type = forceatoms[i++];
 820     ai   = forceatoms[i++];
 821     aj   = forceatoms[i++];
 822     ak   = forceatoms[i++];
 823
 824     theta  = bond_angle(x[ai],x[aj],x[ak],pbc,
 825                         r_ij,r_kj,&cos_theta,&t1,&t2);  /*  41          */
 826
 827     *dvdlambda += harmonic(forceparams[type].harmonic.krA,
 828                            forceparams[type].harmonic.krB,
 829                            forceparams[type].harmonic.rA*DEG2RAD,
 830                            forceparams[type].harmonic.rB*DEG2RAD,
 831                            theta,lambda,&va,&dVdt);  /*  21  */
 832     vtot += va;
 833
 834     cos_theta2 = sqr(cos_theta);
 835     if (cos_theta2 < 1) {
 836       int  m;
 837       real st,sth;
 838       real cik,cii,ckk;
 839       real nrkj2,nrij2;
 840       rvec f_i,f_j,f_k;
 841
 842       st  = dVdt*gmx_invsqrt(1 - cos_theta2);   /*  12          */
 843       sth = st*cos_theta;                       /*   1          */
 844 #ifdef DEBUG
 845       if (debug)
 846         fprintf(debug,"ANGLES: theta = %10g  vth = %10g  dV/dtheta = %10g\n",
 847                 theta*RAD2DEG,va,dVdt);
 848 #endif
 849       nrkj2=iprod(r_kj,r_kj);                   /*   5          */
 850       nrij2=iprod(r_ij,r_ij);
 851
 852       cik=st*gmx_invsqrt(nrkj2*nrij2);          /*  12          */
 853       cii=sth/nrij2;                            /*  10          */
 854       ckk=sth/nrkj2;                            /*  10          */
 855
 856       for (m=0; (m<DIM); m++) {                 /*  39          */
 857         f_i[m]=-(cik*r_kj[m]-cii*r_ij[m]);
 858         f_k[m]=-(cik*r_ij[m]-ckk*r_kj[m]);
 859         f_j[m]=-f_i[m]-f_k[m];
 860         f[ai][m]+=f_i[m];
 861         f[aj][m]+=f_j[m];
 862         f[ak][m]+=f_k[m];
 863       }
 864       if (g) {
 865         copy_ivec(SHIFT_IVEC(g,aj),jt);
 866
 867         ivec_sub(SHIFT_IVEC(g,ai),jt,dt_ij);
 868         ivec_sub(SHIFT_IVEC(g,ak),jt,dt_kj);
 869         t1=IVEC2IS(dt_ij);
 870         t2=IVEC2IS(dt_kj);
 871       }
 872       rvec_inc(fshift[t1],f_i);
 873       rvec_inc(fshift[CENTRAL],f_j);
 874       rvec_inc(fshift[t2],f_k);
 875     }                                           /* 161 TOTAL    */
 876   }
 877   return vtot;
 878 }
 879
 880 real linear_angles(int nbonds,
 881                    const t_iatom forceatoms[],const t_iparams forceparams[],
 882                    const rvec x[],rvec f[],rvec fshift[],
 883                    const t_pbc *pbc,const t_graph *g,
 884                    real lambda,real *dvdlambda,
 885                    const t_mdatoms *md,t_fcdata *fcd,
 886                    int *global_atom_index)
 887 {
 888   int  i,m,ai,aj,ak,t1,t2,type;
 889   rvec f_i,f_j,f_k;
 890   real L1,kA,kB,aA,aB,dr,dr2,va,vtot,a,b,klin;
 891   ivec jt,dt_ij,dt_kj;
 892   rvec r_ij,r_kj,r_ik,dx;
 893
 894   L1   = 1-lambda;
 895   vtot = 0.0;
 896   for(i=0; (i<nbonds); ) {
 897     type = forceatoms[i++];
 898     ai   = forceatoms[i++];
 899     aj   = forceatoms[i++];
 900     ak   = forceatoms[i++];
 901
 902     kA = forceparams[type].linangle.klinA;
 903     kB = forceparams[type].linangle.klinB;
 904     klin = L1*kA + lambda*kB;
 905
 906     aA   = forceparams[type].linangle.aA;
 907     aB   = forceparams[type].linangle.aB;
 908     a    = L1*aA+lambda*aB;
 909     b    = 1-a;
 910
 911     t1 = pbc_rvec_sub(pbc,x[ai],x[aj],r_ij);
 912     t2 = pbc_rvec_sub(pbc,x[ak],x[aj],r_kj);
 913     rvec_sub(r_ij,r_kj,r_ik);
 914
 915     dr2 = 0;
 916     for(m=0; (m<DIM); m++)
 917     {
 918         dr     = - a * r_ij[m] - b * r_kj[m];
 919         dr2   += dr*dr;
 920         dx[m]  = dr;
 921         f_i[m] = a*klin*dr;
 922         f_k[m] = b*klin*dr;
 923         f_j[m] = -(f_i[m]+f_k[m]);
 924         f[ai][m] += f_i[m];
 925         f[aj][m] += f_j[m];
 926         f[ak][m] += f_k[m];
 927     }
 928     va    = 0.5*klin*dr2;
 929     *dvdlambda += 0.5*(kB-kA)*dr2 + klin*(aB-aA)*iprod(dx,r_ik);
 930
 931     vtot += va;
 932
 933     if (g) {
 934         copy_ivec(SHIFT_IVEC(g,aj),jt);
 935
 936         ivec_sub(SHIFT_IVEC(g,ai),jt,dt_ij);
 937         ivec_sub(SHIFT_IVEC(g,ak),jt,dt_kj);
 938         t1=IVEC2IS(dt_ij);
 939         t2=IVEC2IS(dt_kj);
 940     }
 941     rvec_inc(fshift[t1],f_i);
 942     rvec_inc(fshift[CENTRAL],f_j);
 943     rvec_inc(fshift[t2],f_k);
 944   }                                           /* 57 TOTAL       */
 945   return vtot;
 946 }
 947
 948 real urey_bradley(int nbonds,
 949                   const t_iatom forceatoms[],const t_iparams forceparams[],
 950                   const rvec x[],rvec f[],rvec fshift[],
 951                   const t_pbc *pbc,const t_graph *g,
 952                   real lambda,real *dvdlambda,
 953                   const t_mdatoms *md,t_fcdata *fcd,
 954                   int *global_atom_index)
 955 {
 956   int  i,m,ai,aj,ak,t1,t2,type,ki;
 957   rvec r_ij,r_kj,r_ik;
 958   real cos_theta,cos_theta2,theta;
 959   real dVdt,va,vtot,dr,dr2,vbond,fbond,fik;
 960   real kthA,th0A,kUBA,r13A,kthB,th0B,kUBB,r13B;
 961   ivec jt,dt_ij,dt_kj,dt_ik;
 962
 963   vtot = 0.0;
 964   for(i=0; (i<nbonds); ) {
 965     type = forceatoms[i++];
 966     ai   = forceatoms[i++];
 967     aj   = forceatoms[i++];
 968     ak   = forceatoms[i++];
 969     th0A  = forceparams[type].u_b.thetaA*DEG2RAD;
 970     kthA  = forceparams[type].u_b.kthetaA;
 971     r13A  = forceparams[type].u_b.r13A;
 972     kUBA  = forceparams[type].u_b.kUBA;
 973     th0B  = forceparams[type].u_b.thetaB*DEG2RAD;
 974     kthB  = forceparams[type].u_b.kthetaB;
 975     r13B  = forceparams[type].u_b.r13B;
 976     kUBB  = forceparams[type].u_b.kUBB;
 977
 978     theta  = bond_angle(x[ai],x[aj],x[ak],pbc,
 979                         r_ij,r_kj,&cos_theta,&t1,&t2);  /*  41          */
 980
 981     *dvdlambda += harmonic(kthA,kthB,th0A,th0B,theta,lambda,&va,&dVdt);  /*  21  */
 982     vtot += va;
 983
 984     ki   = pbc_rvec_sub(pbc,x[ai],x[ak],r_ik);  /*   3          */
 985     dr2  = iprod(r_ik,r_ik);                    /*   5          */
 986     dr   = dr2*gmx_invsqrt(dr2);                        /*  10          */
 987
 988     *dvdlambda += harmonic(kUBA,kUBB,r13A,r13B,dr,lambda,&vbond,&fbond); /*  19  */
 989
 990     cos_theta2 = sqr(cos_theta);                /*   1          */
 991     if (cos_theta2 < 1) {
 992       real st,sth;
 993       real cik,cii,ckk;
 994       real nrkj2,nrij2;
 995       rvec f_i,f_j,f_k;
 996
 997       st  = dVdt*gmx_invsqrt(1 - cos_theta2);   /*  12          */
 998       sth = st*cos_theta;                       /*   1          */
 999 #ifdef DEBUG
1000       if (debug)
1001         fprintf(debug,"ANGLES: theta = %10g  vth = %10g  dV/dtheta = %10g\n",
1002                 theta*RAD2DEG,va,dVdt);
1003 #endif
1004       nrkj2=iprod(r_kj,r_kj);                   /*   5          */
1005       nrij2=iprod(r_ij,r_ij);
1006
1007       cik=st*gmx_invsqrt(nrkj2*nrij2);          /*  12          */
1008       cii=sth/nrij2;                            /*  10          */
1009       ckk=sth/nrkj2;                            /*  10          */
1010
1011       for (m=0; (m<DIM); m++) {                 /*  39          */
1012         f_i[m]=-(cik*r_kj[m]-cii*r_ij[m]);
1013         f_k[m]=-(cik*r_ij[m]-ckk*r_kj[m]);
1014         f_j[m]=-f_i[m]-f_k[m];
1015         f[ai][m]+=f_i[m];
1016         f[aj][m]+=f_j[m];
1017         f[ak][m]+=f_k[m];
1018       }
1019       if (g) {
1020         copy_ivec(SHIFT_IVEC(g,aj),jt);
1021
1022         ivec_sub(SHIFT_IVEC(g,ai),jt,dt_ij);
1023         ivec_sub(SHIFT_IVEC(g,ak),jt,dt_kj);
1024         t1=IVEC2IS(dt_ij);
1025         t2=IVEC2IS(dt_kj);
1026       }
1027       rvec_inc(fshift[t1],f_i);
1028       rvec_inc(fshift[CENTRAL],f_j);
1029       rvec_inc(fshift[t2],f_k);
1030     }                                           /* 161 TOTAL    */
1031     /* Time for the bond calculations */
1032     if (dr2 == 0.0)
1033       continue;
1034
1035     vtot  += vbond;  /* 1*/
1036     fbond *= gmx_invsqrt(dr2);                  /*   6          */
1037
1038     if (g) {
1039       ivec_sub(SHIFT_IVEC(g,ai),SHIFT_IVEC(g,ak),dt_ik);
1040       ki=IVEC2IS(dt_ik);
1041     }
1042     for (m=0; (m<DIM); m++) {                   /*  15          */
1043       fik=fbond*r_ik[m];
1044       f[ai][m]+=fik;
1045       f[ak][m]-=fik;
1046       fshift[ki][m]+=fik;
1047       fshift[CENTRAL][m]-=fik;
1048     }
1049   }
1050   return vtot;
1051 }
1052
1053 real quartic_angles(int nbonds,
1054                     const t_iatom forceatoms[],const t_iparams forceparams[],
1055                     const rvec x[],rvec f[],rvec fshift[],
1056                     const t_pbc *pbc,const t_graph *g,
1057                     real lambda,real *dvdlambda,
1058                     const t_mdatoms *md,t_fcdata *fcd,
1059                     int *global_atom_index)
1060 {
1061   int  i,j,ai,aj,ak,t1,t2,type;
1062   rvec r_ij,r_kj;
1063   real cos_theta,cos_theta2,theta,dt,dVdt,va,dtp,c,vtot;
1064   ivec jt,dt_ij,dt_kj;
1065
1066   vtot = 0.0;
1067   for(i=0; (i<nbonds); ) {
1068     type = forceatoms[i++];
1069     ai   = forceatoms[i++];
1070     aj   = forceatoms[i++];
1071     ak   = forceatoms[i++];
1072
1073     theta  = bond_angle(x[ai],x[aj],x[ak],pbc,
1074                         r_ij,r_kj,&cos_theta,&t1,&t2);  /*  41          */
1075
1076     dt = theta - forceparams[type].qangle.theta*DEG2RAD; /* 2          */
1077
1078     dVdt = 0;
1079     va = forceparams[type].qangle.c[0];
1080     dtp = 1.0;
1081     for(j=1; j<=4; j++) {
1082       c = forceparams[type].qangle.c[j];
1083       dVdt -= j*c*dtp;
1084       dtp *= dt;
1085       va += c*dtp;
1086     }
1087     /* 20 */
1088
1089     vtot += va;
1090
1091     cos_theta2 = sqr(cos_theta);                /*   1          */
1092     if (cos_theta2 < 1) {
1093       int  m;
1094       real st,sth;
1095       real cik,cii,ckk;
1096       real nrkj2,nrij2;
1097       rvec f_i,f_j,f_k;
1098
1099       st  = dVdt*gmx_invsqrt(1 - cos_theta2);           /*  12          */
1100       sth = st*cos_theta;                       /*   1          */
1101 #ifdef DEBUG
1102       if (debug)
1103         fprintf(debug,"ANGLES: theta = %10g  vth = %10g  dV/dtheta = %10g\n",
1104                 theta*RAD2DEG,va,dVdt);
1105 #endif
1106       nrkj2=iprod(r_kj,r_kj);                   /*   5          */
1107       nrij2=iprod(r_ij,r_ij);
1108
1109       cik=st*gmx_invsqrt(nrkj2*nrij2);          /*  12          */
1110       cii=sth/nrij2;                            /*  10          */
1111       ckk=sth/nrkj2;                            /*  10          */
1112
1113       for (m=0; (m<DIM); m++) {                 /*  39          */
1114         f_i[m]=-(cik*r_kj[m]-cii*r_ij[m]);
1115         f_k[m]=-(cik*r_ij[m]-ckk*r_kj[m]);
1116         f_j[m]=-f_i[m]-f_k[m];
1117         f[ai][m]+=f_i[m];
1118         f[aj][m]+=f_j[m];
1119         f[ak][m]+=f_k[m];
1120       }
1121       if (g) {
1122         copy_ivec(SHIFT_IVEC(g,aj),jt);
1123
1124         ivec_sub(SHIFT_IVEC(g,ai),jt,dt_ij);
1125         ivec_sub(SHIFT_IVEC(g,ak),jt,dt_kj);
1126         t1=IVEC2IS(dt_ij);
1127         t2=IVEC2IS(dt_kj);
1128       }
1129       rvec_inc(fshift[t1],f_i);
1130       rvec_inc(fshift[CENTRAL],f_j);
1131       rvec_inc(fshift[t2],f_k);
1132     }                                           /* 153 TOTAL    */
1133   }
1134   return vtot;
1135 }
1136
1137 real dih_angle(const rvec xi,const rvec xj,const rvec xk,const rvec xl,
1138                const t_pbc *pbc,
1139                rvec r_ij,rvec r_kj,rvec r_kl,rvec m,rvec n,
1140                real *sign,int *t1,int *t2,int *t3)
1141 {
1142   real ipr,phi;
1143
1144   *t1 = pbc_rvec_sub(pbc,xi,xj,r_ij);                   /*  3           */
1145   *t2 = pbc_rvec_sub(pbc,xk,xj,r_kj);                   /*  3           */
1146   *t3 = pbc_rvec_sub(pbc,xk,xl,r_kl);                   /*  3           */
1147
1148   cprod(r_ij,r_kj,m);                   /*  9           */
1149   cprod(r_kj,r_kl,n);                   /*  9           */
1150   phi=gmx_angle(m,n);                   /* 49 (assuming 25 for atan2) */
1151   ipr=iprod(r_ij,n);                    /*  5           */
1152   (*sign)=(ipr<0.0)?-1.0:1.0;
1153   phi=(*sign)*phi;                      /*  1           */
1154                                         /* 82 TOTAL     */
1155   return phi;
1156 }
1157
1158
1159
1160 void do_dih_fup(int i,int j,int k,int l,real ddphi,
1161                 rvec r_ij,rvec r_kj,rvec r_kl,
1162                 rvec m,rvec n,rvec f[],rvec fshift[],
1163                 const t_pbc *pbc,const t_graph *g,
1164                 const rvec x[],int t1,int t2,int t3)
1165 {
1166   /* 143 FLOPS */
1167   rvec f_i,f_j,f_k,f_l;
1168   rvec uvec,vvec,svec,dx_jl;
1169   real iprm,iprn,nrkj,nrkj2;
1170   real a,p,q,toler;
1171   ivec jt,dt_ij,dt_kj,dt_lj;
1172
1173   iprm  = iprod(m,m);           /*  5   */
1174   iprn  = iprod(n,n);           /*  5   */
1175   nrkj2 = iprod(r_kj,r_kj);     /*  5   */
1176   toler = nrkj2*GMX_REAL_EPS;
1177   if ((iprm > toler) && (iprn > toler)) {
1178     nrkj  = nrkj2*gmx_invsqrt(nrkj2);   /* 10   */
1179     a     = -ddphi*nrkj/iprm;   /* 11   */
1180     svmul(a,m,f_i);             /*  3   */
1181     a     = ddphi*nrkj/iprn;    /* 11   */
1182     svmul(a,n,f_l);             /*  3   */
1183     p     = iprod(r_ij,r_kj);   /*  5   */
1184     p    /= nrkj2;              /* 10   */
1185     q     = iprod(r_kl,r_kj);   /*  5   */
1186     q    /= nrkj2;              /* 10   */
1187     svmul(p,f_i,uvec);          /*  3   */
1188     svmul(q,f_l,vvec);          /*  3   */
1189     rvec_sub(uvec,vvec,svec);   /*  3   */
1190     rvec_sub(f_i,svec,f_j);     /*  3   */
1191     rvec_add(f_l,svec,f_k);     /*  3   */
1192     rvec_inc(f[i],f_i);         /*  3   */
1193     rvec_dec(f[j],f_j);         /*  3   */
1194     rvec_dec(f[k],f_k);         /*  3   */
1195     rvec_inc(f[l],f_l);         /*  3   */
1196
1197     if (g) {
1198       copy_ivec(SHIFT_IVEC(g,j),jt);
1199       ivec_sub(SHIFT_IVEC(g,i),jt,dt_ij);
1200       ivec_sub(SHIFT_IVEC(g,k),jt,dt_kj);
1201       ivec_sub(SHIFT_IVEC(g,l),jt,dt_lj);
1202       t1=IVEC2IS(dt_ij);
1203       t2=IVEC2IS(dt_kj);
1204       t3=IVEC2IS(dt_lj);
1205     } else if (pbc) {
1206       t3 = pbc_rvec_sub(pbc,x[l],x[j],dx_jl);
1207     } else {
1208       t3 = CENTRAL;
1209     }
1210
1211     rvec_inc(fshift[t1],f_i);
1212     rvec_dec(fshift[CENTRAL],f_j);
1213     rvec_dec(fshift[t2],f_k);
1214     rvec_inc(fshift[t3],f_l);
1215   }
1216   /* 112 TOTAL  */
1217 }
1218
1219
1220 real dopdihs(real cpA,real cpB,real phiA,real phiB,int mult,
1221              real phi,real lambda,real *V,real *F)
1222 {
1223   real v,dvdlambda,mdphi,v1,sdphi,ddphi;
1224   real L1   = 1.0 - lambda;
1225   real ph0  = (L1*phiA + lambda*phiB)*DEG2RAD;
1226   real dph0 = (phiB - phiA)*DEG2RAD;
1227   real cp   = L1*cpA + lambda*cpB;
1228
1229   mdphi =  mult*phi - ph0;
1230   sdphi = sin(mdphi);
1231   ddphi = -cp*mult*sdphi;
1232   v1    = 1.0 + cos(mdphi);
1233   v     = cp*v1;
1234
1235   dvdlambda  = (cpB - cpA)*v1 + cp*dph0*sdphi;
1236
1237   *V = v;
1238   *F = ddphi;
1239
1240   return dvdlambda;
1241
1242   /* That was 40 flops */
1243 }
1244
1245 static real dopdihs_min(real cpA,real cpB,real phiA,real phiB,int mult,
1246                         real phi,real lambda,real *V,real *F)
1247      /* similar to dopdihs, except for a minus sign  *
1248       * and a different treatment of mult/phi0       */
1249 {
1250   real v,dvdlambda,mdphi,v1,sdphi,ddphi;
1251   real L1   = 1.0 - lambda;
1252   real ph0  = (L1*phiA + lambda*phiB)*DEG2RAD;
1253   real dph0 = (phiB - phiA)*DEG2RAD;
1254   real cp   = L1*cpA + lambda*cpB;
1255
1256   mdphi = mult*(phi-ph0);
1257   sdphi = sin(mdphi);
1258   ddphi = cp*mult*sdphi;
1259   v1    = 1.0-cos(mdphi);
1260   v     = cp*v1;
1261
1262   dvdlambda  = (cpB-cpA)*v1 + cp*dph0*sdphi;
1263
1264   *V = v;
1265   *F = ddphi;
1266
1267   return dvdlambda;
1268
1269   /* That was 40 flops */
1270 }
1271
1272 real pdihs(int nbonds,
1273            const t_iatom forceatoms[],const t_iparams forceparams[],
1274            const rvec x[],rvec f[],rvec fshift[],
1275            const t_pbc *pbc,const t_graph *g,
1276            real lambda,real *dvdlambda,
1277            const t_mdatoms *md,t_fcdata *fcd,
1278            int *global_atom_index)
1279 {
1280   int  i,type,ai,aj,ak,al;
1281   int  t1,t2,t3;
1282   rvec r_ij,r_kj,r_kl,m,n;
1283   real phi,sign,ddphi,vpd,vtot;
1284
1285   vtot = 0.0;
1286
1287   for(i=0; (i<nbonds); ) {
1288     type = forceatoms[i++];
1289     ai   = forceatoms[i++];
1290     aj   = forceatoms[i++];
1291     ak   = forceatoms[i++];
1292     al   = forceatoms[i++];
1293
1294     phi=dih_angle(x[ai],x[aj],x[ak],x[al],pbc,r_ij,r_kj,r_kl,m,n,
1295                   &sign,&t1,&t2,&t3);                   /*  84          */
1296     *dvdlambda += dopdihs(forceparams[type].pdihs.cpA,
1297                           forceparams[type].pdihs.cpB,
1298                           forceparams[type].pdihs.phiA,
1299                           forceparams[type].pdihs.phiB,
1300                           forceparams[type].pdihs.mult,
1301                           phi,lambda,&vpd,&ddphi);
1302
1303     vtot += vpd;
1304     do_dih_fup(ai,aj,ak,al,ddphi,r_ij,r_kj,r_kl,m,n,
1305                f,fshift,pbc,g,x,t1,t2,t3);                      /* 112          */
1306
1307 #ifdef DEBUG
1308     fprintf(debug,"pdih: (%d,%d,%d,%d) phi=%g\n",
1309             ai,aj,ak,al,phi);
1310 #endif
1311   } /* 223 TOTAL        */
1312
1313   return vtot;
1314 }
1315
1316 void make_dp_periodic(real *dp)  /* 1 flop? */
1317 {
1318     /* dp cannot be outside (-pi,pi) */
1319     if (*dp >= M_PI)
1320     {
1321         *dp -= 2*M_PI;
1322     }
1323     else if (*dp < -M_PI)
1324     {
1325         *dp += 2*M_PI;
1326     }
1327     return;
1328 }
1329
1330
1331 real idihs(int nbonds,
1332            const t_iatom forceatoms[],const t_iparams forceparams[],
1333            const rvec x[],rvec f[],rvec fshift[],
1334            const t_pbc *pbc,const t_graph *g,
1335            real lambda,real *dvdlambda,
1336            const t_mdatoms *md,t_fcdata *fcd,
1337            int *global_atom_index)
1338 {
1339   int  i,type,ai,aj,ak,al;
1340   int  t1,t2,t3;
1341   real phi,phi0,dphi0,ddphi,sign,vtot;
1342   rvec r_ij,r_kj,r_kl,m,n;
1343   real L1,kk,dp,dp2,kA,kB,pA,pB,dvdl_term;
1344
1345   L1 = 1.0-lambda;
1346   dvdl_term = 0;
1347   vtot = 0.0;
1348   for(i=0; (i<nbonds); ) {
1349     type = forceatoms[i++];
1350     ai   = forceatoms[i++];
1351     aj   = forceatoms[i++];
1352     ak   = forceatoms[i++];
1353     al   = forceatoms[i++];
1354
1355     phi=dih_angle(x[ai],x[aj],x[ak],x[al],pbc,r_ij,r_kj,r_kl,m,n,
1356                   &sign,&t1,&t2,&t3);                   /*  84          */
1357
1358     /* phi can jump if phi0 is close to Pi/-Pi, which will cause huge
1359      * force changes if we just apply a normal harmonic.
1360      * Instead, we first calculate phi-phi0 and take it modulo (-Pi,Pi).
1361      * This means we will never have the periodicity problem, unless
1362      * the dihedral is Pi away from phiO, which is very unlikely due to
1363      * the potential.
1364      */
1365     kA = forceparams[type].harmonic.krA;
1366     kB = forceparams[type].harmonic.krB;
1367     pA = forceparams[type].harmonic.rA;
1368     pB = forceparams[type].harmonic.rB;
1369
1370     kk    = L1*kA + lambda*kB;
1371     phi0  = (L1*pA + lambda*pB)*DEG2RAD;
1372     dphi0 = (pB - pA)*DEG2RAD;
1373
1374     dp = phi-phi0;
1375
1376     make_dp_periodic(&dp);
1377
1378     dp2 = dp*dp;
1379
1380     vtot += 0.5*kk*dp2;
1381     ddphi = -kk*dp;
1382
1383     dvdl_term += 0.5*(kB - kA)*dp2 - kk*dphi0*dp;
1384
1385     do_dih_fup(ai,aj,ak,al,(real)(-ddphi),r_ij,r_kj,r_kl,m,n,
1386                f,fshift,pbc,g,x,t1,t2,t3);                      /* 112          */
1387     /* 218 TOTAL        */
1388 #ifdef DEBUG
1389     if (debug)
1390       fprintf(debug,"idih: (%d,%d,%d,%d) phi=%g\n",
1391               ai,aj,ak,al,phi);
1392 #endif
1393   }
1394
1395   *dvdlambda += dvdl_term;
1396   return vtot;
1397 }
1398
1399
1400 real posres(int nbonds,
1401             const t_iatom forceatoms[],const t_iparams forceparams[],
1402             const rvec x[],rvec f[],rvec vir_diag,
1403             t_pbc *pbc,
1404             real lambda,real *dvdlambda,
1405             int refcoord_scaling,int ePBC,rvec comA,rvec comB)
1406 {
1407     int  i,ai,m,d,type,ki,npbcdim=0;
1408     const t_iparams *pr;
1409     real L1;
1410     real vtot,kk,fm;
1411     real posA,posB,ref=0;
1412     rvec comA_sc,comB_sc,rdist,dpdl,pos,dx;
1413     gmx_bool bForceValid = TRUE;
1414
1415     if ((f==NULL) || (vir_diag==NULL)) {  /* should both be null together! */
1416         bForceValid = FALSE;
1417     }
1418
1419     npbcdim = ePBC2npbcdim(ePBC);
1420
1421     if (refcoord_scaling == erscCOM)
1422     {
1423         clear_rvec(comA_sc);
1424         clear_rvec(comB_sc);
1425         for(m=0; m<npbcdim; m++)
1426         {
1427             for(d=m; d<npbcdim; d++)
1428             {
1429                 comA_sc[m] += comA[d]*pbc->box[d][m];
1430                 comB_sc[m] += comB[d]*pbc->box[d][m];
1431             }
1432         }
1433     }
1434
1435     L1 = 1.0 - lambda;
1436
1437     vtot = 0.0;
1438     for(i=0; (i<nbonds); )
1439     {
1440         type = forceatoms[i++];
1441         ai   = forceatoms[i++];
1442         pr   = &forceparams[type];
1443
1444         for(m=0; m<DIM; m++)
1445         {
1446             posA = forceparams[type].posres.pos0A[m];
1447             posB = forceparams[type].posres.pos0B[m];
1448             if (m < npbcdim)
1449             {
1450                 switch (refcoord_scaling)
1451                 {
1452                 case erscNO:
1453                     ref      = 0;
1454                     rdist[m] = L1*posA + lambda*posB;
1455                     dpdl[m]  = posB - posA;
1456                     break;
1457                 case erscALL:
1458                     /* Box relative coordinates are stored for dimensions with pbc */
1459                     posA *= pbc->box[m][m];
1460                     posB *= pbc->box[m][m];
1461                     for(d=m+1; d<npbcdim; d++)
1462                     {
1463                         posA += forceparams[type].posres.pos0A[d]*pbc->box[d][m];
1464                         posB += forceparams[type].posres.pos0B[d]*pbc->box[d][m];
1465                     }
1466                     ref      = L1*posA + lambda*posB;
1467                     rdist[m] = 0;
1468                     dpdl[m]  = posB - posA;
1469                     break;
1470                 case erscCOM:
1471                     ref      = L1*comA_sc[m] + lambda*comB_sc[m];
1472                     rdist[m] = L1*posA       + lambda*posB;
1473                     dpdl[m]  = comB_sc[m] - comA_sc[m] + posB - posA;
1474                     break;
1475                 default:
1476                     gmx_fatal(FARGS, "No such scaling method implemented");
1477                 }
1478             }
1479             else
1480             {
1481                 ref      = L1*posA + lambda*posB;
1482                 rdist[m] = 0;
1483                 dpdl[m]  = posB - posA;
1484             }
1485
1486             /* We do pbc_dx with ref+rdist,
1487              * since with only ref we can be up to half a box vector wrong.
1488              */
1489             pos[m] = ref + rdist[m];
1490         }
1491
1492         if (pbc)
1493         {
1494             pbc_dx(pbc,x[ai],pos,dx);
1495         }
1496         else
1497         {
1498             rvec_sub(x[ai],pos,dx);
1499         }
1500
1501         for (m=0; (m<DIM); m++)
1502         {
1503             kk          = L1*pr->posres.fcA[m] + lambda*pr->posres.fcB[m];
1504             fm          = -kk*dx[m];
1505             vtot       += 0.5*kk*dx[m]*dx[m];
1506             *dvdlambda +=
1507                 0.5*(pr->posres.fcB[m] - pr->posres.fcA[m])*dx[m]*dx[m]
1508                 -fm*dpdl[m];
1509
1510             /* Here we correct for the pbc_dx which included rdist */
1511             if (bForceValid) {
1512                 f[ai][m]   += fm;
1513                 vir_diag[m] -= 0.5*(dx[m] + rdist[m])*fm;
1514             }
1515         }
1516     }
1517
1518     return vtot;
1519 }
1520
1521 static real low_angres(int nbonds,
1522                        const t_iatom forceatoms[],const t_iparams forceparams[],
1523                        const rvec x[],rvec f[],rvec fshift[],
1524                        const t_pbc *pbc,const t_graph *g,
1525                        real lambda,real *dvdlambda,
1526                        gmx_bool bZAxis)
1527 {
1528   int  i,m,type,ai,aj,ak,al;
1529   int  t1,t2;
1530   real phi,cos_phi,cos_phi2,vid,vtot,dVdphi;
1531   rvec r_ij,r_kl,f_i,f_k={0,0,0};
1532   real st,sth,nrij2,nrkl2,c,cij,ckl;
1533
1534   ivec dt;
1535   t2 = 0; /* avoid warning with gcc-3.3. It is never used uninitialized */
1536
1537   vtot = 0.0;
1538   ak=al=0; /* to avoid warnings */
1539   for(i=0; i<nbonds; ) {
1540     type = forceatoms[i++];
1541     ai   = forceatoms[i++];
1542     aj   = forceatoms[i++];
1543     t1   = pbc_rvec_sub(pbc,x[aj],x[ai],r_ij);                  /*  3           */
1544     if (!bZAxis) {
1545       ak   = forceatoms[i++];
1546       al   = forceatoms[i++];
1547       t2   = pbc_rvec_sub(pbc,x[al],x[ak],r_kl);           /*  3                */
1548     } else {
1549       r_kl[XX] = 0;
1550       r_kl[YY] = 0;
1551       r_kl[ZZ] = 1;
1552     }
1553
1554     cos_phi = cos_angle(r_ij,r_kl);             /* 25           */
1555     phi     = acos(cos_phi);                    /* 10           */
1556
1557     *dvdlambda += dopdihs_min(forceparams[type].pdihs.cpA,
1558                               forceparams[type].pdihs.cpB,
1559                               forceparams[type].pdihs.phiA,
1560                               forceparams[type].pdihs.phiB,
1561                               forceparams[type].pdihs.mult,
1562                               phi,lambda,&vid,&dVdphi); /*  40  */
1563
1564     vtot += vid;
1565
1566     cos_phi2 = sqr(cos_phi);                    /*   1          */
1567     if (cos_phi2 < 1) {
1568       st  = -dVdphi*gmx_invsqrt(1 - cos_phi2);      /*  12              */
1569       sth = st*cos_phi;                         /*   1          */
1570       nrij2 = iprod(r_ij,r_ij);                 /*   5          */
1571       nrkl2 = iprod(r_kl,r_kl);                 /*   5          */
1572
1573       c   = st*gmx_invsqrt(nrij2*nrkl2);                /*  11          */
1574       cij = sth/nrij2;                          /*  10          */
1575       ckl = sth/nrkl2;                          /*  10          */
1576
1577       for (m=0; m<DIM; m++) {                   /*  18+18       */
1578         f_i[m] = (c*r_kl[m]-cij*r_ij[m]);
1579         f[ai][m] += f_i[m];
1580         f[aj][m] -= f_i[m];
1581         if (!bZAxis) {
1582           f_k[m] = (c*r_ij[m]-ckl*r_kl[m]);
1583           f[ak][m] += f_k[m];
1584           f[al][m] -= f_k[m];
1585         }
1586       }
1587
1588       if (g) {
1589         ivec_sub(SHIFT_IVEC(g,ai),SHIFT_IVEC(g,aj),dt);
1590         t1=IVEC2IS(dt);
1591       }
1592       rvec_inc(fshift[t1],f_i);
1593       rvec_dec(fshift[CENTRAL],f_i);
1594       if (!bZAxis) {
1595         if (g) {
1596           ivec_sub(SHIFT_IVEC(g,ak),SHIFT_IVEC(g,al),dt);
1597           t2=IVEC2IS(dt);
1598         }
1599         rvec_inc(fshift[t2],f_k);
1600         rvec_dec(fshift[CENTRAL],f_k);
1601       }
1602     }
1603   }
1604
1605   return vtot;  /*  184 / 157 (bZAxis)  total  */
1606 }
1607
1608 real angres(int nbonds,
1609             const t_iatom forceatoms[],const t_iparams forceparams[],
1610             const rvec x[],rvec f[],rvec fshift[],
1611             const t_pbc *pbc,const t_graph *g,
1612             real lambda,real *dvdlambda,
1613             const t_mdatoms *md,t_fcdata *fcd,
1614             int *global_atom_index)
1615 {
1616   return low_angres(nbonds,forceatoms,forceparams,x,f,fshift,pbc,g,
1617                     lambda,dvdlambda,FALSE);
1618 }
1619
1620 real angresz(int nbonds,
1621              const t_iatom forceatoms[],const t_iparams forceparams[],
1622              const rvec x[],rvec f[],rvec fshift[],
1623              const t_pbc *pbc,const t_graph *g,
1624              real lambda,real *dvdlambda,
1625              const t_mdatoms *md,t_fcdata *fcd,
1626              int *global_atom_index)
1627 {
1628   return low_angres(nbonds,forceatoms,forceparams,x,f,fshift,pbc,g,
1629                     lambda,dvdlambda,TRUE);
1630 }
1631
1632 real dihres(int nbonds,
1633             const t_iatom forceatoms[],const t_iparams forceparams[],
1634             const rvec x[],rvec f[],rvec fshift[],
1635             const t_pbc *pbc,const t_graph *g,
1636             real lambda,real *dvdlambda,
1637             const t_mdatoms *md,t_fcdata *fcd,
1638             int *global_atom_index)
1639 {
1640     real vtot = 0;
1641     int  ai,aj,ak,al,i,k,type,t1,t2,t3;
1642     real phi0A,phi0B,dphiA,dphiB,kfacA,kfacB,phi0,dphi,kfac;
1643     real phi,ddphi,ddp,ddp2,dp,sign,d2r,fc,L1;
1644     rvec r_ij,r_kj,r_kl,m,n;
1645
1646     L1 = 1.0-lambda;
1647
1648     d2r = DEG2RAD;
1649     k   = 0;
1650
1651     for (i=0; (i<nbonds); )
1652     {
1653         type = forceatoms[i++];
1654         ai   = forceatoms[i++];
1655         aj   = forceatoms[i++];
1656         ak   = forceatoms[i++];
1657         al   = forceatoms[i++];
1658
1659         phi0A  = forceparams[type].dihres.phiA*d2r;
1660         dphiA  = forceparams[type].dihres.dphiA*d2r;
1661         kfacA  = forceparams[type].dihres.kfacA;
1662
1663         phi0B  = forceparams[type].dihres.phiB*d2r;
1664         dphiB  = forceparams[type].dihres.dphiB*d2r;
1665         kfacB  = forceparams[type].dihres.kfacB;
1666
1667         phi0  = L1*phi0A + lambda*phi0B;
1668         dphi  = L1*dphiA + lambda*dphiB;
1669         kfac = L1*kfacA + lambda*kfacB;
1670
1671         phi = dih_angle(x[ai],x[aj],x[ak],x[al],pbc,r_ij,r_kj,r_kl,m,n,
1672                         &sign,&t1,&t2,&t3);
1673         /* 84 flops */
1674
1675         if (debug)
1676         {
1677             fprintf(debug,"dihres[%d]: %d %d %d %d : phi=%f, dphi=%f, kfac=%f\n",
1678                     k++,ai,aj,ak,al,phi0,dphi,kfac);
1679         }
1680         /* phi can jump if phi0 is close to Pi/-Pi, which will cause huge
1681          * force changes if we just apply a normal harmonic.
1682          * Instead, we first calculate phi-phi0 and take it modulo (-Pi,Pi).
1683          * This means we will never have the periodicity problem, unless
1684          * the dihedral is Pi away from phiO, which is very unlikely due to
1685          * the potential.
1686          */
1687         dp = phi-phi0;
1688         make_dp_periodic(&dp);
1689
1690         if (dp > dphi)
1691         {
1692             ddp = dp-dphi;
1693         }
1694         else if (dp < -dphi)
1695         {
1696             ddp = dp+dphi;
1697         }
1698         else
1699         {
1700             ddp = 0;
1701         }
1702
1703         if (ddp != 0.0)
1704         {
1705             ddp2 = ddp*ddp;
1706             vtot += 0.5*kfac*ddp2;
1707             ddphi = kfac*ddp;
1708
1709             *dvdlambda += 0.5*(kfacB - kfacA)*ddp2;
1710             /* lambda dependence from changing restraint distances */
1711             if (ddp > 0)
1712             {
1713                 *dvdlambda -= kfac*ddp*((dphiB - dphiA)+(phi0B - phi0A));
1714             }
1715             else if (ddp < 0)
1716             {
1717                 *dvdlambda += kfac*ddp*((dphiB - dphiA)-(phi0B - phi0A));
1718             }
1719             do_dih_fup(ai,aj,ak,al,ddphi,r_ij,r_kj,r_kl,m,n,
1720                        f,fshift,pbc,g,x,t1,t2,t3);              /* 112          */
1721         }
1722     }
1723     return vtot;
1724 }
1725
1726
1727 real unimplemented(int nbonds,
1728                    const t_iatom forceatoms[],const t_iparams forceparams[],
1729                    const rvec x[],rvec f[],rvec fshift[],
1730                    const t_pbc *pbc,const t_graph *g,
1731                    real lambda,real *dvdlambda,
1732                    const t_mdatoms *md,t_fcdata *fcd,
1733                    int *global_atom_index)
1734 {
1735   gmx_impl("*** you are using a not implemented function");
1736
1737   return 0.0; /* To make the compiler happy */
1738 }
1739
1740 real rbdihs(int nbonds,
1741             const t_iatom forceatoms[],const t_iparams forceparams[],
1742             const rvec x[],rvec f[],rvec fshift[],
1743             const t_pbc *pbc,const t_graph *g,
1744             real lambda,real *dvdlambda,
1745             const t_mdatoms *md,t_fcdata *fcd,
1746             int *global_atom_index)
1747 {
1748   const real c0=0.0,c1=1.0,c2=2.0,c3=3.0,c4=4.0,c5=5.0;
1749   int  type,ai,aj,ak,al,i,j;
1750   int  t1,t2,t3;
1751   rvec r_ij,r_kj,r_kl,m,n;
1752   real parmA[NR_RBDIHS];
1753   real parmB[NR_RBDIHS];
1754   real parm[NR_RBDIHS];
1755   real cos_phi,phi,rbp,rbpBA;
1756   real v,sign,ddphi,sin_phi;
1757   real cosfac,vtot;
1758   real L1   = 1.0-lambda;
1759   real dvdl_term=0;
1760
1761   vtot = 0.0;
1762   for(i=0; (i<nbonds); ) {
1763     type = forceatoms[i++];
1764     ai   = forceatoms[i++];
1765     aj   = forceatoms[i++];
1766     ak   = forceatoms[i++];
1767     al   = forceatoms[i++];
1768
1769       phi=dih_angle(x[ai],x[aj],x[ak],x[al],pbc,r_ij,r_kj,r_kl,m,n,
1770                     &sign,&t1,&t2,&t3);                 /*  84          */
1771
1772     /* Change to polymer convention */
1773     if (phi < c0)
1774       phi += M_PI;
1775     else
1776       phi -= M_PI;                      /*   1          */
1777
1778     cos_phi = cos(phi);
1779     /* Beware of accuracy loss, cannot use 1-sqrt(cos^2) ! */
1780     sin_phi = sin(phi);
1781
1782     for(j=0; (j<NR_RBDIHS); j++) {
1783       parmA[j] = forceparams[type].rbdihs.rbcA[j];
1784       parmB[j] = forceparams[type].rbdihs.rbcB[j];
1785       parm[j]  = L1*parmA[j]+lambda*parmB[j];
1786     }
1787     /* Calculate cosine powers */
1788     /* Calculate the energy */
1789     /* Calculate the derivative */
1790
1791     v       = parm[0];
1792     dvdl_term   += (parmB[0]-parmA[0]);
1793     ddphi   = c0;
1794     cosfac  = c1;
1795
1796     rbp     = parm[1];
1797     rbpBA   = parmB[1]-parmA[1];
1798     ddphi  += rbp*cosfac;
1799     cosfac *= cos_phi;
1800     v      += cosfac*rbp;
1801     dvdl_term   += cosfac*rbpBA;
1802     rbp     = parm[2];
1803     rbpBA   = parmB[2]-parmA[2];
1804     ddphi  += c2*rbp*cosfac;
1805     cosfac *= cos_phi;
1806     v      += cosfac*rbp;
1807     dvdl_term   += cosfac*rbpBA;
1808     rbp     = parm[3];
1809     rbpBA   = parmB[3]-parmA[3];
1810     ddphi  += c3*rbp*cosfac;
1811     cosfac *= cos_phi;
1812     v      += cosfac*rbp;
1813     dvdl_term   += cosfac*rbpBA;
1814     rbp     = parm[4];
1815     rbpBA   = parmB[4]-parmA[4];
1816     ddphi  += c4*rbp*cosfac;
1817     cosfac *= cos_phi;
1818     v      += cosfac*rbp;
1819     dvdl_term   += cosfac*rbpBA;
1820     rbp     = parm[5];
1821     rbpBA   = parmB[5]-parmA[5];
1822     ddphi  += c5*rbp*cosfac;
1823     cosfac *= cos_phi;
1824     v      += cosfac*rbp;
1825     dvdl_term   += cosfac*rbpBA;
1826
1827     ddphi = -ddphi*sin_phi;                             /*  11          */
1828
1829     do_dih_fup(ai,aj,ak,al,ddphi,r_ij,r_kj,r_kl,m,n,
1830                f,fshift,pbc,g,x,t1,t2,t3);              /* 112          */
1831     vtot += v;
1832   }
1833   *dvdlambda += dvdl_term;
1834
1835   return vtot;
1836 }
1837
1838 int cmap_setup_grid_index(int ip, int grid_spacing, int *ipm1, int *ipp1, int *ipp2)
1839 {
1840         int im1, ip1, ip2;
1841
1842         if(ip<0)
1843         {
1844                 ip = ip + grid_spacing - 1;
1845         }
1846         else if(ip > grid_spacing)
1847         {
1848                 ip = ip - grid_spacing - 1;
1849         }
1850
1851         im1 = ip - 1;
1852         ip1 = ip + 1;
1853         ip2 = ip + 2;
1854
1855         if(ip == 0)
1856         {
1857                 im1 = grid_spacing - 1;
1858         }
1859         else if(ip == grid_spacing-2)
1860         {
1861                 ip2 = 0;
1862         }
1863         else if(ip == grid_spacing-1)
1864         {
1865                 ip1 = 0;
1866                 ip2 = 1;
1867         }
1868
1869         *ipm1 = im1;
1870         *ipp1 = ip1;
1871         *ipp2 = ip2;
1872
1873         return ip;
1874
1875 }
1876
1877 real cmap_dihs(int nbonds,
1878                            const t_iatom forceatoms[],const t_iparams forceparams[],
1879                const gmx_cmap_t *cmap_grid,
1880                            const rvec x[],rvec f[],rvec fshift[],
1881                            const t_pbc *pbc,const t_graph *g,
1882                            real lambda,real *dvdlambda,
1883                            const t_mdatoms *md,t_fcdata *fcd,
1884                            int *global_atom_index)
1885 {
1886         int i,j,k,n,idx;
1887         int ai,aj,ak,al,am;
1888         int a1i,a1j,a1k,a1l,a2i,a2j,a2k,a2l;
1889         int type,cmapA;
1890         int t11,t21,t31,t12,t22,t32;
1891         int iphi1,ip1m1,ip1p1,ip1p2;
1892         int iphi2,ip2m1,ip2p1,ip2p2;
1893         int l1,l2,l3,l4;
1894         int pos1,pos2,pos3,pos4,tmp;
1895
1896         real ty[4],ty1[4],ty2[4],ty12[4],tc[16],tx[16];
1897         real phi1,psi1,cos_phi1,sin_phi1,sign1,xphi1;
1898         real phi2,psi2,cos_phi2,sin_phi2,sign2,xphi2;
1899         real dx,xx,tt,tu,e,df1,df2,ddf1,ddf2,ddf12,vtot;
1900         real ra21,rb21,rg21,rg1,rgr1,ra2r1,rb2r1,rabr1;
1901         real ra22,rb22,rg22,rg2,rgr2,ra2r2,rb2r2,rabr2;
1902         real fg1,hg1,fga1,hgb1,gaa1,gbb1;
1903         real fg2,hg2,fga2,hgb2,gaa2,gbb2;
1904         real fac;
1905
1906         rvec r1_ij, r1_kj, r1_kl,m1,n1;
1907         rvec r2_ij, r2_kj, r2_kl,m2,n2;
1908         rvec f1_i,f1_j,f1_k,f1_l;
1909         rvec f2_i,f2_j,f2_k,f2_l;
1910         rvec a1,b1,a2,b2;
1911         rvec f1,g1,h1,f2,g2,h2;
1912         rvec dtf1,dtg1,dth1,dtf2,dtg2,dth2;
1913         ivec jt1,dt1_ij,dt1_kj,dt1_lj;
1914         ivec jt2,dt2_ij,dt2_kj,dt2_lj;
1915
1916     const real *cmapd;
1917
1918         int loop_index[4][4] = {
1919                 {0,4,8,12},
1920                 {1,5,9,13},
1921                 {2,6,10,14},
1922                 {3,7,11,15}
1923         };
1924
1925         /* Total CMAP energy */
1926         vtot = 0;
1927
1928         for(n=0;n<nbonds; )
1929         {
1930                 /* Five atoms are involved in the two torsions */
1931                 type   = forceatoms[n++];
1932                 ai     = forceatoms[n++];
1933                 aj     = forceatoms[n++];
1934                 ak     = forceatoms[n++];
1935                 al     = forceatoms[n++];
1936                 am     = forceatoms[n++];
1937
1938                 /* Which CMAP type is this */
1939                 cmapA = forceparams[type].cmap.cmapA;
1940         cmapd = cmap_grid->cmapdata[cmapA].cmap;
1941
1942                 /* First torsion */
1943                 a1i   = ai;
1944                 a1j   = aj;
1945                 a1k   = ak;
1946                 a1l   = al;
1947
1948                 phi1  = dih_angle(x[a1i], x[a1j], x[a1k], x[a1l], pbc, r1_ij, r1_kj, r1_kl, m1, n1,
1949                                                    &sign1, &t11, &t21, &t31); /* 84 */
1950
1951         cos_phi1 = cos(phi1);
1952
1953                 a1[0] = r1_ij[1]*r1_kj[2]-r1_ij[2]*r1_kj[1];
1954                 a1[1] = r1_ij[2]*r1_kj[0]-r1_ij[0]*r1_kj[2];
1955                 a1[2] = r1_ij[0]*r1_kj[1]-r1_ij[1]*r1_kj[0]; /* 9 */
1956
1957                 b1[0] = r1_kl[1]*r1_kj[2]-r1_kl[2]*r1_kj[1];
1958                 b1[1] = r1_kl[2]*r1_kj[0]-r1_kl[0]*r1_kj[2];
1959                 b1[2] = r1_kl[0]*r1_kj[1]-r1_kl[1]*r1_kj[0]; /* 9 */
1960
1961                 tmp = pbc_rvec_sub(pbc,x[a1l],x[a1k],h1);
1962
1963                 ra21  = iprod(a1,a1);       /* 5 */
1964                 rb21  = iprod(b1,b1);       /* 5 */
1965                 rg21  = iprod(r1_kj,r1_kj); /* 5 */
1966                 rg1   = sqrt(rg21);
1967
1968                 rgr1  = 1.0/rg1;
1969                 ra2r1 = 1.0/ra21;
1970                 rb2r1 = 1.0/rb21;
1971                 rabr1 = sqrt(ra2r1*rb2r1);
1972
1973                 sin_phi1 = rg1 * rabr1 * iprod(a1,h1) * (-1);
1974
1975                 if(cos_phi1 < -0.5 || cos_phi1 > 0.5)
1976                 {
1977                         phi1 = asin(sin_phi1);
1978
1979                         if(cos_phi1 < 0)
1980                         {
1981                                 if(phi1 > 0)
1982                                 {
1983                                         phi1 = M_PI - phi1;
1984                                 }
1985                                 else
1986                                 {
1987                                         phi1 = -M_PI - phi1;
1988                                 }
1989                         }
1990                 }
1991                 else
1992                 {
1993                         phi1 = acos(cos_phi1);
1994
1995                         if(sin_phi1 < 0)
1996                         {
1997                                 phi1 = -phi1;
1998                         }
1999                 }
2000
2001                 xphi1 = phi1 + M_PI; /* 1 */
2002
2003                 /* Second torsion */
2004                 a2i   = aj;
2005                 a2j   = ak;
2006                 a2k   = al;
2007                 a2l   = am;
2008
2009                 phi2  = dih_angle(x[a2i], x[a2j], x[a2k], x[a2l], pbc, r2_ij, r2_kj, r2_kl, m2, n2,
2010                                                   &sign2, &t12, &t22, &t32); /* 84 */
2011
2012         cos_phi2 = cos(phi2);
2013
2014                 a2[0] = r2_ij[1]*r2_kj[2]-r2_ij[2]*r2_kj[1];
2015                 a2[1] = r2_ij[2]*r2_kj[0]-r2_ij[0]*r2_kj[2];
2016                 a2[2] = r2_ij[0]*r2_kj[1]-r2_ij[1]*r2_kj[0]; /* 9 */
2017
2018                 b2[0] = r2_kl[1]*r2_kj[2]-r2_kl[2]*r2_kj[1];
2019                 b2[1] = r2_kl[2]*r2_kj[0]-r2_kl[0]*r2_kj[2];
2020                 b2[2] = r2_kl[0]*r2_kj[1]-r2_kl[1]*r2_kj[0]; /* 9 */
2021
2022                 tmp = pbc_rvec_sub(pbc,x[a2l],x[a2k],h2);
2023
2024                 ra22  = iprod(a2,a2);         /* 5 */
2025                 rb22  = iprod(b2,b2);         /* 5 */
2026                 rg22  = iprod(r2_kj,r2_kj);   /* 5 */
2027                 rg2   = sqrt(rg22);
2028
2029                 rgr2  = 1.0/rg2;
2030                 ra2r2 = 1.0/ra22;
2031                 rb2r2 = 1.0/rb22;
2032                 rabr2 = sqrt(ra2r2*rb2r2);
2033
2034                 sin_phi2 = rg2 * rabr2 * iprod(a2,h2) * (-1);
2035
2036                 if(cos_phi2 < -0.5 || cos_phi2 > 0.5)
2037                 {
2038                         phi2 = asin(sin_phi2);
2039
2040                         if(cos_phi2 < 0)
2041                         {
2042                                 if(phi2 > 0)
2043                                 {
2044                                         phi2 = M_PI - phi2;
2045                                 }
2046                                 else
2047                                 {
2048                                         phi2 = -M_PI - phi2;
2049                                 }
2050                         }
2051                 }
2052                 else
2053                 {
2054                         phi2 = acos(cos_phi2);
2055
2056                         if(sin_phi2 < 0)
2057                         {
2058                                 phi2 = -phi2;
2059                         }
2060                 }
2061
2062                 xphi2 = phi2 + M_PI; /* 1 */
2063
2064                 /* Range mangling */
2065                 if(xphi1<0)
2066                 {
2067                         xphi1 = xphi1 + 2*M_PI;
2068                 }
2069                 else if(xphi1>=2*M_PI)
2070                 {
2071                         xphi1 = xphi1 - 2*M_PI;
2072                 }
2073
2074                 if(xphi2<0)
2075                 {
2076                         xphi2 = xphi2 + 2*M_PI;
2077                 }
2078                 else if(xphi2>=2*M_PI)
2079                 {
2080                         xphi2 = xphi2 - 2*M_PI;
2081                 }
2082
2083                 /* Number of grid points */
2084                 dx = 2*M_PI / cmap_grid->grid_spacing;
2085
2086                 /* Where on the grid are we */
2087                 iphi1 = (int)(xphi1/dx);
2088                 iphi2 = (int)(xphi2/dx);
2089
2090                 iphi1 = cmap_setup_grid_index(iphi1, cmap_grid->grid_spacing, &ip1m1,&ip1p1,&ip1p2);
2091                 iphi2 = cmap_setup_grid_index(iphi2, cmap_grid->grid_spacing, &ip2m1,&ip2p1,&ip2p2);
2092
2093                 pos1    = iphi1*cmap_grid->grid_spacing+iphi2;
2094                 pos2    = ip1p1*cmap_grid->grid_spacing+iphi2;
2095                 pos3    = ip1p1*cmap_grid->grid_spacing+ip2p1;
2096                 pos4    = iphi1*cmap_grid->grid_spacing+ip2p1;
2097
2098                 ty[0]   = cmapd[pos1*4];
2099                 ty[1]   = cmapd[pos2*4];
2100                 ty[2]   = cmapd[pos3*4];
2101                 ty[3]   = cmapd[pos4*4];
2102
2103                 ty1[0]   = cmapd[pos1*4+1];
2104                 ty1[1]   = cmapd[pos2*4+1];
2105                 ty1[2]   = cmapd[pos3*4+1];
2106                 ty1[3]   = cmapd[pos4*4+1];
2107
2108                 ty2[0]   = cmapd[pos1*4+2];
2109                 ty2[1]   = cmapd[pos2*4+2];
2110                 ty2[2]   = cmapd[pos3*4+2];
2111                 ty2[3]   = cmapd[pos4*4+2];
2112
2113                 ty12[0]   = cmapd[pos1*4+3];
2114                 ty12[1]   = cmapd[pos2*4+3];
2115                 ty12[2]   = cmapd[pos3*4+3];
2116                 ty12[3]   = cmapd[pos4*4+3];
2117
2118                 /* Switch to degrees */
2119                 dx = 360.0 / cmap_grid->grid_spacing;
2120                 xphi1 = xphi1 * RAD2DEG;
2121                 xphi2 = xphi2 * RAD2DEG;
2122
2123                 for(i=0;i<4;i++) /* 16 */
2124                 {
2125                         tx[i] = ty[i];
2126                         tx[i+4] = ty1[i]*dx;
2127                         tx[i+8] = ty2[i]*dx;
2128                         tx[i+12] = ty12[i]*dx*dx;
2129                 }
2130
2131                 idx=0;
2132                 for(i=0;i<4;i++) /* 1056 */
2133                 {
2134                         for(j=0;j<4;j++)
2135                         {
2136                                 xx = 0;
2137                                 for(k=0;k<16;k++)
2138                                 {
2139                                         xx = xx + cmap_coeff_matrix[k*16+idx]*tx[k];
2140                                 }
2141
2142                                 idx++;
2143                                 tc[i*4+j]=xx;
2144                         }
2145                 }
2146
2147                 tt    = (xphi1-iphi1*dx)/dx;
2148                 tu    = (xphi2-iphi2*dx)/dx;
2149
2150                 e     = 0;
2151                 df1   = 0;
2152                 df2   = 0;
2153                 ddf1  = 0;
2154                 ddf2  = 0;
2155                 ddf12 = 0;
2156
2157                 for(i=3;i>=0;i--)
2158                 {
2159                         l1 = loop_index[i][3];
2160                         l2 = loop_index[i][2];
2161                         l3 = loop_index[i][1];
2162
2163                         e     = tt * e    + ((tc[i*4+3]*tu+tc[i*4+2])*tu + tc[i*4+1])*tu+tc[i*4];
2164                         df1   = tu * df1  + (3.0*tc[l1]*tt+2.0*tc[l2])*tt+tc[l3];
2165                         df2   = tt * df2  + (3.0*tc[i*4+3]*tu+2.0*tc[i*4+2])*tu+tc[i*4+1];
2166                         ddf1  = tu * ddf1 + 2.0*3.0*tc[l1]*tt+2.0*tc[l2];
2167                         ddf2  = tt * ddf2 + 2.0*3.0*tc[4*i+3]*tu+2.0*tc[4*i+2];
2168                 }
2169
2170                 ddf12 = tc[5] + 2.0*tc[9]*tt + 3.0*tc[13]*tt*tt + 2.0*tu*(tc[6]+2.0*tc[10]*tt+3.0*tc[14]*tt*tt) +
2171                 3.0*tu*tu*(tc[7]+2.0*tc[11]*tt+3.0*tc[15]*tt*tt);
2172
2173                 fac     = RAD2DEG/dx;
2174                 df1     = df1   * fac;
2175                 df2     = df2   * fac;
2176                 ddf1    = ddf1  * fac * fac;
2177                 ddf2    = ddf2  * fac * fac;
2178                 ddf12   = ddf12 * fac * fac;
2179
2180                 /* CMAP energy */
2181                 vtot += e;
2182
2183                 /* Do forces - first torsion */
2184                 fg1       = iprod(r1_ij,r1_kj);
2185                 hg1       = iprod(r1_kl,r1_kj);
2186                 fga1      = fg1*ra2r1*rgr1;
2187                 hgb1      = hg1*rb2r1*rgr1;
2188                 gaa1      = -ra2r1*rg1;
2189                 gbb1      = rb2r1*rg1;
2190
2191                 for(i=0;i<DIM;i++)
2192                 {
2193                         dtf1[i]   = gaa1 * a1[i];
2194                         dtg1[i]   = fga1 * a1[i] - hgb1 * b1[i];
2195                         dth1[i]   = gbb1 * b1[i];
2196
2197                         f1[i]     = df1  * dtf1[i];
2198                         g1[i]     = df1  * dtg1[i];
2199                         h1[i]     = df1  * dth1[i];
2200
2201                         f1_i[i]   =  f1[i];
2202                         f1_j[i]   = -f1[i] - g1[i];
2203                         f1_k[i]   =  h1[i] + g1[i];
2204                         f1_l[i]   = -h1[i];
2205
2206                         f[a1i][i] = f[a1i][i] + f1_i[i];
2207                         f[a1j][i] = f[a1j][i] + f1_j[i]; /* - f1[i] - g1[i] */
2208                         f[a1k][i] = f[a1k][i] + f1_k[i]; /* h1[i] + g1[i] */
2209                         f[a1l][i] = f[a1l][i] + f1_l[i]; /* h1[i] */
2210                 }
2211
2212                 /* Do forces - second torsion */
2213                 fg2       = iprod(r2_ij,r2_kj);
2214                 hg2       = iprod(r2_kl,r2_kj);
2215                 fga2      = fg2*ra2r2*rgr2;
2216                 hgb2      = hg2*rb2r2*rgr2;
2217                 gaa2      = -ra2r2*rg2;
2218                 gbb2      = rb2r2*rg2;
2219
2220                 for(i=0;i<DIM;i++)
2221                 {
2222                         dtf2[i]   = gaa2 * a2[i];
2223                         dtg2[i]   = fga2 * a2[i] - hgb2 * b2[i];
2224                         dth2[i]   = gbb2 * b2[i];
2225
2226                         f2[i]     = df2  * dtf2[i];
2227                         g2[i]     = df2  * dtg2[i];
2228                         h2[i]     = df2  * dth2[i];
2229
2230                         f2_i[i]   =  f2[i];
2231                         f2_j[i]   = -f2[i] - g2[i];
2232                         f2_k[i]   =  h2[i] + g2[i];
2233                         f2_l[i]   = -h2[i];
2234
2235                         f[a2i][i] = f[a2i][i] + f2_i[i]; /* f2[i] */
2236                         f[a2j][i] = f[a2j][i] + f2_j[i]; /* - f2[i] - g2[i] */
2237                         f[a2k][i] = f[a2k][i] + f2_k[i]; /* h2[i] + g2[i] */
2238                         f[a2l][i] = f[a2l][i] + f2_l[i]; /* - h2[i] */
2239                 }
2240
2241                 /* Shift forces */
2242                 if(g)
2243                 {
2244                         copy_ivec(SHIFT_IVEC(g,a1j), jt1);
2245                         ivec_sub(SHIFT_IVEC(g,a1i),  jt1,dt1_ij);
2246                         ivec_sub(SHIFT_IVEC(g,a1k),  jt1,dt1_kj);
2247                         ivec_sub(SHIFT_IVEC(g,a1l),  jt1,dt1_lj);
2248                         t11 = IVEC2IS(dt1_ij);
2249                         t21 = IVEC2IS(dt1_kj);
2250                         t31 = IVEC2IS(dt1_lj);
2251
2252                         copy_ivec(SHIFT_IVEC(g,a2j), jt2);
2253                         ivec_sub(SHIFT_IVEC(g,a2i),  jt2,dt2_ij);
2254                         ivec_sub(SHIFT_IVEC(g,a2k),  jt2,dt2_kj);
2255                         ivec_sub(SHIFT_IVEC(g,a2l),  jt2,dt2_lj);
2256                         t12 = IVEC2IS(dt2_ij);
2257                         t22 = IVEC2IS(dt2_kj);
2258                         t32 = IVEC2IS(dt2_lj);
2259                 }
2260                 else if(pbc)
2261                 {
2262                         t31 = pbc_rvec_sub(pbc,x[a1l],x[a1j],h1);
2263                         t32 = pbc_rvec_sub(pbc,x[a2l],x[a2j],h2);
2264                 }
2265                 else
2266                 {
2267                         t31 = CENTRAL;
2268                         t32 = CENTRAL;
2269                 }
2270
2271                 rvec_inc(fshift[t11],f1_i);
2272                 rvec_inc(fshift[CENTRAL],f1_j);
2273                 rvec_inc(fshift[t21],f1_k);
2274                 rvec_inc(fshift[t31],f1_l);
2275
2276                 rvec_inc(fshift[t21],f2_i);
2277                 rvec_inc(fshift[CENTRAL],f2_j);
2278                 rvec_inc(fshift[t22],f2_k);
2279                 rvec_inc(fshift[t32],f2_l);
2280         }
2281         return vtot;
2282 }
2283
2284
2285
2286 /***********************************************************
2287  *
2288  *   G R O M O S  9 6   F U N C T I O N S
2289  *
2290  ***********************************************************/
2291 real g96harmonic(real kA,real kB,real xA,real xB,real x,real lambda,
2292                  real *V,real *F)
2293 {
2294   const real half=0.5;
2295   real  L1,kk,x0,dx,dx2;
2296   real  v,f,dvdlambda;
2297
2298   L1    = 1.0-lambda;
2299   kk    = L1*kA+lambda*kB;
2300   x0    = L1*xA+lambda*xB;
2301
2302   dx    = x-x0;
2303   dx2   = dx*dx;
2304
2305   f     = -kk*dx;
2306   v     = half*kk*dx2;
2307   dvdlambda  = half*(kB-kA)*dx2 + (xA-xB)*kk*dx;
2308
2309   *F    = f;
2310   *V    = v;
2311
2312   return dvdlambda;
2313
2314   /* That was 21 flops */
2315 }
2316
2317 real g96bonds(int nbonds,
2318               const t_iatom forceatoms[],const t_iparams forceparams[],
2319               const rvec x[],rvec f[],rvec fshift[],
2320               const t_pbc *pbc,const t_graph *g,
2321               real lambda,real *dvdlambda,
2322               const t_mdatoms *md,t_fcdata *fcd,
2323               int *global_atom_index)
2324 {
2325   int  i,m,ki,ai,aj,type;
2326   real dr2,fbond,vbond,fij,vtot;
2327   rvec dx;
2328   ivec dt;
2329
2330   vtot = 0.0;
2331   for(i=0; (i<nbonds); ) {
2332     type = forceatoms[i++];
2333     ai   = forceatoms[i++];
2334     aj   = forceatoms[i++];
2335
2336     ki   = pbc_rvec_sub(pbc,x[ai],x[aj],dx);            /*   3          */
2337     dr2  = iprod(dx,dx);                                /*   5          */
2338
2339     *dvdlambda += g96harmonic(forceparams[type].harmonic.krA,
2340                               forceparams[type].harmonic.krB,
2341                               forceparams[type].harmonic.rA,
2342                               forceparams[type].harmonic.rB,
2343                               dr2,lambda,&vbond,&fbond);
2344
2345     vtot  += 0.5*vbond;                             /* 1*/
2346 #ifdef DEBUG
2347     if (debug)
2348       fprintf(debug,"G96-BONDS: dr = %10g  vbond = %10g  fbond = %10g\n",
2349               sqrt(dr2),vbond,fbond);
2350 #endif
2351
2352     if (g) {
2353       ivec_sub(SHIFT_IVEC(g,ai),SHIFT_IVEC(g,aj),dt);
2354       ki=IVEC2IS(dt);
2355     }
2356     for (m=0; (m<DIM); m++) {                   /*  15          */
2357       fij=fbond*dx[m];
2358       f[ai][m]+=fij;
2359       f[aj][m]-=fij;
2360       fshift[ki][m]+=fij;
2361       fshift[CENTRAL][m]-=fij;
2362     }
2363   }                                     /* 44 TOTAL     */
2364   return vtot;
2365 }
2366
2367 real g96bond_angle(const rvec xi,const rvec xj,const rvec xk,const t_pbc *pbc,
2368                    rvec r_ij,rvec r_kj,
2369                    int *t1,int *t2)
2370 /* Return value is the angle between the bonds i-j and j-k */
2371 {
2372   real costh;
2373
2374   *t1 = pbc_rvec_sub(pbc,xi,xj,r_ij);                   /*  3           */
2375   *t2 = pbc_rvec_sub(pbc,xk,xj,r_kj);                   /*  3           */
2376
2377   costh=cos_angle(r_ij,r_kj);           /* 25           */
2378                                         /* 41 TOTAL     */
2379   return costh;
2380 }
2381
2382 real g96angles(int nbonds,
2383                const t_iatom forceatoms[],const t_iparams forceparams[],
2384                const rvec x[],rvec f[],rvec fshift[],
2385                const t_pbc *pbc,const t_graph *g,
2386                real lambda,real *dvdlambda,
2387                const t_mdatoms *md,t_fcdata *fcd,
2388                int *global_atom_index)
2389 {
2390   int  i,ai,aj,ak,type,m,t1,t2;
2391   rvec r_ij,r_kj;
2392   real cos_theta,dVdt,va,vtot;
2393   real rij_1,rij_2,rkj_1,rkj_2,rijrkj_1;
2394   rvec f_i,f_j,f_k;
2395   ivec jt,dt_ij,dt_kj;
2396
2397   vtot = 0.0;
2398   for(i=0; (i<nbonds); ) {
2399     type = forceatoms[i++];
2400     ai   = forceatoms[i++];
2401     aj   = forceatoms[i++];
2402     ak   = forceatoms[i++];
2403
2404     cos_theta  = g96bond_angle(x[ai],x[aj],x[ak],pbc,r_ij,r_kj,&t1,&t2);
2405
2406     *dvdlambda += g96harmonic(forceparams[type].harmonic.krA,
2407                               forceparams[type].harmonic.krB,
2408                               forceparams[type].harmonic.rA,
2409                               forceparams[type].harmonic.rB,
2410                               cos_theta,lambda,&va,&dVdt);
2411     vtot    += va;
2412
2413     rij_1    = gmx_invsqrt(iprod(r_ij,r_ij));
2414     rkj_1    = gmx_invsqrt(iprod(r_kj,r_kj));
2415     rij_2    = rij_1*rij_1;
2416     rkj_2    = rkj_1*rkj_1;
2417     rijrkj_1 = rij_1*rkj_1;                     /* 23 */
2418
2419 #ifdef DEBUG
2420     if (debug)
2421       fprintf(debug,"G96ANGLES: costheta = %10g  vth = %10g  dV/dct = %10g\n",
2422               cos_theta,va,dVdt);
2423 #endif
2424     for (m=0; (m<DIM); m++) {                   /*  42  */
2425       f_i[m]=dVdt*(r_kj[m]*rijrkj_1 - r_ij[m]*rij_2*cos_theta);
2426       f_k[m]=dVdt*(r_ij[m]*rijrkj_1 - r_kj[m]*rkj_2*cos_theta);
2427       f_j[m]=-f_i[m]-f_k[m];
2428       f[ai][m]+=f_i[m];
2429       f[aj][m]+=f_j[m];
2430       f[ak][m]+=f_k[m];
2431     }
2432
2433     if (g) {
2434       copy_ivec(SHIFT_IVEC(g,aj),jt);
2435
2436       ivec_sub(SHIFT_IVEC(g,ai),jt,dt_ij);
2437       ivec_sub(SHIFT_IVEC(g,ak),jt,dt_kj);
2438       t1=IVEC2IS(dt_ij);
2439       t2=IVEC2IS(dt_kj);
2440     }
2441     rvec_inc(fshift[t1],f_i);
2442     rvec_inc(fshift[CENTRAL],f_j);
2443     rvec_inc(fshift[t2],f_k);               /* 9 */
2444     /* 163 TOTAL        */
2445   }
2446   return vtot;
2447 }
2448
2449 real cross_bond_bond(int nbonds,
2450                      const t_iatom forceatoms[],const t_iparams forceparams[],
2451                      const rvec x[],rvec f[],rvec fshift[],
2452                      const t_pbc *pbc,const t_graph *g,
2453                      real lambda,real *dvdlambda,
2454                      const t_mdatoms *md,t_fcdata *fcd,
2455                      int *global_atom_index)
2456 {
2457   /* Potential from Lawrence and Skimmer, Chem. Phys. Lett. 372 (2003)
2458    * pp. 842-847
2459    */
2460   int  i,ai,aj,ak,type,m,t1,t2;
2461   rvec r_ij,r_kj;
2462   real vtot,vrr,s1,s2,r1,r2,r1e,r2e,krr;
2463   rvec f_i,f_j,f_k;
2464   ivec jt,dt_ij,dt_kj;
2465
2466   vtot = 0.0;
2467   for(i=0; (i<nbonds); ) {
2468     type = forceatoms[i++];
2469     ai   = forceatoms[i++];
2470     aj   = forceatoms[i++];
2471     ak   = forceatoms[i++];
2472     r1e  = forceparams[type].cross_bb.r1e;
2473     r2e  = forceparams[type].cross_bb.r2e;
2474     krr  = forceparams[type].cross_bb.krr;
2475
2476     /* Compute distance vectors ... */
2477     t1 = pbc_rvec_sub(pbc,x[ai],x[aj],r_ij);
2478     t2 = pbc_rvec_sub(pbc,x[ak],x[aj],r_kj);
2479
2480     /* ... and their lengths */
2481     r1 = norm(r_ij);
2482     r2 = norm(r_kj);
2483
2484     /* Deviations from ideality */
2485     s1 = r1-r1e;
2486     s2 = r2-r2e;
2487
2488     /* Energy (can be negative!) */
2489     vrr   = krr*s1*s2;
2490     vtot += vrr;
2491
2492     /* Forces */
2493     svmul(-krr*s2/r1,r_ij,f_i);
2494     svmul(-krr*s1/r2,r_kj,f_k);
2495
2496     for (m=0; (m<DIM); m++) {                   /*  12  */
2497       f_j[m]    = -f_i[m] - f_k[m];
2498       f[ai][m] += f_i[m];
2499       f[aj][m] += f_j[m];
2500       f[ak][m] += f_k[m];
2501     }
2502
2503     /* Virial stuff */
2504     if (g) {
2505       copy_ivec(SHIFT_IVEC(g,aj),jt);
2506
2507       ivec_sub(SHIFT_IVEC(g,ai),jt,dt_ij);
2508       ivec_sub(SHIFT_IVEC(g,ak),jt,dt_kj);
2509       t1=IVEC2IS(dt_ij);
2510       t2=IVEC2IS(dt_kj);
2511     }
2512     rvec_inc(fshift[t1],f_i);
2513     rvec_inc(fshift[CENTRAL],f_j);
2514     rvec_inc(fshift[t2],f_k);               /* 9 */
2515     /* 163 TOTAL        */
2516   }
2517   return vtot;
2518 }
2519
2520 real cross_bond_angle(int nbonds,
2521                       const t_iatom forceatoms[],const t_iparams forceparams[],
2522                       const rvec x[],rvec f[],rvec fshift[],
2523                       const t_pbc *pbc,const t_graph *g,
2524                       real lambda,real *dvdlambda,
2525                       const t_mdatoms *md,t_fcdata *fcd,
2526                       int *global_atom_index)
2527 {
2528   /* Potential from Lawrence and Skimmer, Chem. Phys. Lett. 372 (2003)
2529    * pp. 842-847
2530    */
2531   int  i,ai,aj,ak,type,m,t1,t2,t3;
2532   rvec r_ij,r_kj,r_ik;
2533   real vtot,vrt,s1,s2,s3,r1,r2,r3,r1e,r2e,r3e,krt,k1,k2,k3;
2534   rvec f_i,f_j,f_k;
2535   ivec jt,dt_ij,dt_kj;
2536
2537   vtot = 0.0;
2538   for(i=0; (i<nbonds); ) {
2539     type = forceatoms[i++];
2540     ai   = forceatoms[i++];
2541     aj   = forceatoms[i++];
2542     ak   = forceatoms[i++];
2543     r1e  = forceparams[type].cross_ba.r1e;
2544     r2e  = forceparams[type].cross_ba.r2e;
2545     r3e  = forceparams[type].cross_ba.r3e;
2546     krt  = forceparams[type].cross_ba.krt;
2547
2548     /* Compute distance vectors ... */
2549     t1 = pbc_rvec_sub(pbc,x[ai],x[aj],r_ij);
2550     t2 = pbc_rvec_sub(pbc,x[ak],x[aj],r_kj);
2551     t3 = pbc_rvec_sub(pbc,x[ai],x[ak],r_ik);
2552
2553     /* ... and their lengths */
2554     r1 = norm(r_ij);
2555     r2 = norm(r_kj);
2556     r3 = norm(r_ik);
2557
2558     /* Deviations from ideality */
2559     s1 = r1-r1e;
2560     s2 = r2-r2e;
2561     s3 = r3-r3e;
2562
2563     /* Energy (can be negative!) */
2564     vrt   = krt*s3*(s1+s2);
2565     vtot += vrt;
2566
2567     /* Forces */
2568     k1 = -krt*(s3/r1);
2569     k2 = -krt*(s3/r2);
2570     k3 = -krt*(s1+s2)/r3;
2571     for(m=0; (m<DIM); m++) {
2572       f_i[m] = k1*r_ij[m] + k3*r_ik[m];
2573       f_k[m] = k2*r_kj[m] - k3*r_ik[m];
2574       f_j[m] = -f_i[m] - f_k[m];
2575     }
2576
2577     for (m=0; (m<DIM); m++) {                   /*  12  */
2578       f[ai][m] += f_i[m];
2579       f[aj][m] += f_j[m];
2580       f[ak][m] += f_k[m];
2581     }
2582
2583     /* Virial stuff */
2584     if (g) {
2585       copy_ivec(SHIFT_IVEC(g,aj),jt);
2586
2587       ivec_sub(SHIFT_IVEC(g,ai),jt,dt_ij);
2588       ivec_sub(SHIFT_IVEC(g,ak),jt,dt_kj);
2589       t1=IVEC2IS(dt_ij);
2590       t2=IVEC2IS(dt_kj);
2591     }
2592     rvec_inc(fshift[t1],f_i);
2593     rvec_inc(fshift[CENTRAL],f_j);
2594     rvec_inc(fshift[t2],f_k);               /* 9 */
2595     /* 163 TOTAL        */
2596   }
2597   return vtot;
2598 }
2599
2600 static real bonded_tab(const char *type,int table_nr,
2601                        const bondedtable_t *table,real kA,real kB,real r,
2602                        real lambda,real *V,real *F)
2603 {
2604   real k,tabscale,*VFtab,rt,eps,eps2,Yt,Ft,Geps,Heps2,Fp,VV,FF;
2605   int  n0,nnn;
2606   real v,f,dvdlambda;
2607
2608   k = (1.0 - lambda)*kA + lambda*kB;
2609
2610   tabscale = table->scale;
2611   VFtab    = table->tab;
2612
2613   rt    = r*tabscale;
2614   n0    = rt;
2615   if (n0 >= table->n) {
2616     gmx_fatal(FARGS,"A tabulated %s interaction table number %d is out of the table range: r %f, between table indices %d and %d, table length %d",
2617               type,table_nr,r,n0,n0+1,table->n);
2618   }
2619   eps   = rt - n0;
2620   eps2  = eps*eps;
2621   nnn   = 4*n0;
2622   Yt    = VFtab[nnn];
2623   Ft    = VFtab[nnn+1];
2624   Geps  = VFtab[nnn+2]*eps;
2625   Heps2 = VFtab[nnn+3]*eps2;
2626   Fp    = Ft + Geps + Heps2;
2627   VV    = Yt + Fp*eps;
2628   FF    = Fp + Geps + 2.0*Heps2;
2629
2630   *F    = -k*FF*tabscale;
2631   *V    = k*VV;
2632   dvdlambda  = (kB - kA)*VV;
2633
2634   return dvdlambda;
2635
2636   /* That was 22 flops */
2637 }
2638
2639 real tab_bonds(int nbonds,
2640                const t_iatom forceatoms[],const t_iparams forceparams[],
2641                const rvec x[],rvec f[],rvec fshift[],
2642                const t_pbc *pbc,const t_graph *g,
2643                real lambda,real *dvdlambda,
2644                const t_mdatoms *md,t_fcdata *fcd,
2645                int *global_atom_index)
2646 {
2647   int  i,m,ki,ai,aj,type,table;
2648   real dr,dr2,fbond,vbond,fij,vtot;
2649   rvec dx;
2650   ivec dt;
2651
2652   vtot = 0.0;
2653   for(i=0; (i<nbonds); ) {
2654     type = forceatoms[i++];
2655     ai   = forceatoms[i++];
2656     aj   = forceatoms[i++];
2657
2658     ki   = pbc_rvec_sub(pbc,x[ai],x[aj],dx);    /*   3          */
2659     dr2  = iprod(dx,dx);                        /*   5          */
2660     dr   = dr2*gmx_invsqrt(dr2);                        /*  10          */
2661
2662     table = forceparams[type].tab.table;
2663
2664     *dvdlambda += bonded_tab("bond",table,
2665                              &fcd->bondtab[table],
2666                              forceparams[type].tab.kA,
2667                              forceparams[type].tab.kB,
2668                              dr,lambda,&vbond,&fbond);  /*  22 */
2669
2670     if (dr2 == 0.0)
2671       continue;
2672
2673
2674     vtot  += vbond;/* 1*/
2675     fbond *= gmx_invsqrt(dr2);                  /*   6          */
2676 #ifdef DEBUG
2677     if (debug)
2678       fprintf(debug,"TABBONDS: dr = %10g  vbond = %10g  fbond = %10g\n",
2679               dr,vbond,fbond);
2680 #endif
2681     if (g) {
2682       ivec_sub(SHIFT_IVEC(g,ai),SHIFT_IVEC(g,aj),dt);
2683       ki=IVEC2IS(dt);
2684     }
2685     for (m=0; (m<DIM); m++) {                   /*  15          */
2686       fij=fbond*dx[m];
2687       f[ai][m]+=fij;
2688       f[aj][m]-=fij;
2689       fshift[ki][m]+=fij;
2690       fshift[CENTRAL][m]-=fij;
2691     }
2692   }                                     /* 62 TOTAL     */
2693   return vtot;
2694 }
2695
2696 real tab_angles(int nbonds,
2697                 const t_iatom forceatoms[],const t_iparams forceparams[],
2698                 const rvec x[],rvec f[],rvec fshift[],
2699                 const t_pbc *pbc,const t_graph *g,
2700                 real lambda,real *dvdlambda,
2701                 const t_mdatoms *md,t_fcdata *fcd,
2702                 int *global_atom_index)
2703 {
2704   int  i,ai,aj,ak,t1,t2,type,table;
2705   rvec r_ij,r_kj;
2706   real cos_theta,cos_theta2,theta,dVdt,va,vtot;
2707   ivec jt,dt_ij,dt_kj;
2708
2709   vtot = 0.0;
2710   for(i=0; (i<nbonds); ) {
2711     type = forceatoms[i++];
2712     ai   = forceatoms[i++];
2713     aj   = forceatoms[i++];
2714     ak   = forceatoms[i++];
2715
2716     theta  = bond_angle(x[ai],x[aj],x[ak],pbc,
2717                         r_ij,r_kj,&cos_theta,&t1,&t2);  /*  41          */
2718
2719     table = forceparams[type].tab.table;
2720
2721     *dvdlambda += bonded_tab("angle",table,
2722                              &fcd->angletab[table],
2723                              forceparams[type].tab.kA,
2724                              forceparams[type].tab.kB,
2725                              theta,lambda,&va,&dVdt);  /*  22  */
2726     vtot += va;
2727
2728     cos_theta2 = sqr(cos_theta);                /*   1          */
2729     if (cos_theta2 < 1) {
2730       int  m;
2731       real snt,st,sth;
2732       real cik,cii,ckk;
2733       real nrkj2,nrij2;
2734       rvec f_i,f_j,f_k;
2735
2736       st  = dVdt*gmx_invsqrt(1 - cos_theta2);   /*  12          */
2737       sth = st*cos_theta;                       /*   1          */
2738 #ifdef DEBUG
2739       if (debug)
2740         fprintf(debug,"ANGLES: theta = %10g  vth = %10g  dV/dtheta = %10g\n",
2741                 theta*RAD2DEG,va,dVdt);
2742 #endif
2743       nrkj2=iprod(r_kj,r_kj);                   /*   5          */
2744       nrij2=iprod(r_ij,r_ij);
2745
2746       cik=st*gmx_invsqrt(nrkj2*nrij2);          /*  12          */
2747       cii=sth/nrij2;                            /*  10          */
2748       ckk=sth/nrkj2;                            /*  10          */
2749
2750       for (m=0; (m<DIM); m++) {                 /*  39          */
2751         f_i[m]=-(cik*r_kj[m]-cii*r_ij[m]);
2752         f_k[m]=-(cik*r_ij[m]-ckk*r_kj[m]);
2753         f_j[m]=-f_i[m]-f_k[m];
2754         f[ai][m]+=f_i[m];
2755         f[aj][m]+=f_j[m];
2756         f[ak][m]+=f_k[m];
2757       }
2758       if (g) {
2759         copy_ivec(SHIFT_IVEC(g,aj),jt);
2760
2761         ivec_sub(SHIFT_IVEC(g,ai),jt,dt_ij);
2762         ivec_sub(SHIFT_IVEC(g,ak),jt,dt_kj);
2763         t1=IVEC2IS(dt_ij);
2764         t2=IVEC2IS(dt_kj);
2765       }
2766       rvec_inc(fshift[t1],f_i);
2767       rvec_inc(fshift[CENTRAL],f_j);
2768       rvec_inc(fshift[t2],f_k);
2769     }                                           /* 169 TOTAL    */
2770   }
2771   return vtot;
2772 }
2773
2774 real tab_dihs(int nbonds,
2775               const t_iatom forceatoms[],const t_iparams forceparams[],
2776               const rvec x[],rvec f[],rvec fshift[],
2777               const t_pbc *pbc,const t_graph *g,
2778               real lambda,real *dvdlambda,
2779               const t_mdatoms *md,t_fcdata *fcd,
2780               int *global_atom_index)
2781 {
2782   int  i,type,ai,aj,ak,al,table;
2783   int  t1,t2,t3;
2784   rvec r_ij,r_kj,r_kl,m,n;
2785   real phi,sign,ddphi,vpd,vtot;
2786
2787   vtot = 0.0;
2788   for(i=0; (i<nbonds); ) {
2789     type = forceatoms[i++];
2790     ai   = forceatoms[i++];
2791     aj   = forceatoms[i++];
2792     ak   = forceatoms[i++];
2793     al   = forceatoms[i++];
2794
2795     phi=dih_angle(x[ai],x[aj],x[ak],x[al],pbc,r_ij,r_kj,r_kl,m,n,
2796                   &sign,&t1,&t2,&t3);                   /*  84  */
2797
2798     table = forceparams[type].tab.table;
2799
2800     /* Hopefully phi+M_PI never results in values < 0 */
2801     *dvdlambda += bonded_tab("dihedral",table,
2802                              &fcd->dihtab[table],
2803                              forceparams[type].tab.kA,
2804                              forceparams[type].tab.kB,
2805                              phi+M_PI,lambda,&vpd,&ddphi);
2806
2807     vtot += vpd;
2808     do_dih_fup(ai,aj,ak,al,-ddphi,r_ij,r_kj,r_kl,m,n,
2809                f,fshift,pbc,g,x,t1,t2,t3);                      /* 112  */
2810
2811 #ifdef DEBUG
2812     fprintf(debug,"pdih: (%d,%d,%d,%d) phi=%g\n",
2813             ai,aj,ak,al,phi);
2814 #endif
2815   } /* 227 TOTAL        */
2816
2817   return vtot;
2818 }
2819
2820 real calc_one_bond(FILE *fplog,int ftype, const t_idef *idef,
2821                    rvec x[], rvec f[], t_forcerec *fr,
2822                    const t_pbc *pbc,const t_graph *g,
2823                    gmx_enerdata_t *enerd, t_nrnb *nrnb,
2824                    real *lambda, real *dvdl,
2825                    const t_mdatoms *md,t_fcdata *fcd,
2826                    int *global_atom_index, gmx_bool bPrintSepPot)
2827 {
2828     int ind,nat1,nbonds,efptFTYPE;
2829     real v=0;
2830     t_iatom *iatoms;
2831
2832     if (IS_RESTRAINT_TYPE(ftype))
2833     {
2834         efptFTYPE = efptRESTRAINT;
2835     }
2836     else
2837     {
2838         efptFTYPE = efptBONDED;
2839     }
2840
2841     if (ftype<F_GB12 || ftype>F_GB14)
2842     {
2843         if (interaction_function[ftype].flags & IF_BOND &&
2844             !(ftype == F_CONNBONDS || ftype == F_POSRES))
2845         {
2846             ind  = interaction_function[ftype].nrnb_ind;
2847             nat1 = interaction_function[ftype].nratoms+1;
2848             nbonds    = idef->il[ftype].nr;
2849             iatoms    = idef->il[ftype].iatoms;
2850             if (nbonds > 0)
2851             {
2852                 if (ftype < F_LJ14 || ftype > F_LJC_PAIRS_NB)
2853                 {
2854                     if(ftype==F_CMAP)
2855                     {
2856                         v = cmap_dihs(nbonds,iatoms,
2857                                       idef->iparams,&idef->cmap_grid,
2858                                       (const rvec*)x,f,fr->fshift,
2859                                       pbc,g,lambda[efptFTYPE],&(dvdl[efptFTYPE]),
2860                                       md,fcd,global_atom_index);
2861                     }
2862                     else
2863                     {
2864                         v =         interaction_function[ftype].ifunc(nbonds,iatoms,
2865                                                                   idef->iparams,
2866                                                                   (const rvec*)x,f,fr->fshift,
2867                                                                   pbc,g,lambda[efptFTYPE],&(dvdl[efptFTYPE]),
2868                                                                   md,fcd,global_atom_index);
2869                     }
2870                     enerd->dvdl_nonlin[efptFTYPE] += dvdl[efptFTYPE];
2871                     if (bPrintSepPot)
2872                     {
2873                         fprintf(fplog,"  %-23s #%4d  V %12.5e  dVdl %12.5e\n",
2874                                 interaction_function[ftype].longname,
2875                                 nbonds/nat1,v,lambda[efptFTYPE]);
2876                     }
2877                 }
2878                 else
2879                 {
2880                     v = do_listed_vdw_q(ftype,nbonds,iatoms,
2881                                         idef->iparams,
2882                                         (const rvec*)x,f,fr->fshift,
2883                                         pbc,g,lambda,dvdl,
2884                                         md,fr,&enerd->grpp,global_atom_index);
2885                     enerd->dvdl_nonlin[efptCOUL] += dvdl[efptCOUL];
2886                     enerd->dvdl_nonlin[efptVDW] += dvdl[efptVDW];
2887
2888                     if (bPrintSepPot)
2889                     {
2890                         fprintf(fplog,"  %-5s + %-15s #%4d                  dVdl %12.5e\n",
2891                                 interaction_function[ftype].longname,
2892                                 interaction_function[F_LJ14].longname,nbonds/nat1,dvdl[efptVDW]);
2893                         fprintf(fplog,"  %-5s + %-15s #%4d                  dVdl %12.5e\n",
2894                                 interaction_function[ftype].longname,
2895                                 interaction_function[F_COUL14].longname,nbonds/nat1,dvdl[efptCOUL]);
2896                     }
2897                 }
2898                 if (ind != -1)
2899                 {
2900                     inc_nrnb(nrnb,ind,nbonds/nat1);
2901                 }
2902             }
2903         }
2904     }
2905     return v;
2906 }
2907
2908 /* WARNING!  THIS FUNCTION MUST EXACTLY TRACK THE calc_one_bond
2909    function, or horrible things will happen when doing free energy
2910    calculations!  In a good coding world, this would not be a
2911    different function, but for speed reasons, it needs to be made a
2912    separate function.  TODO for 5.0 - figure out a way to reorganize
2913    to reduce duplication.
2914 */
2915
2916 real calc_one_bond_foreign(FILE *fplog,int ftype, const t_idef *idef,
2917                            rvec x[], rvec f[], t_forcerec *fr,
2918                            const t_pbc *pbc,const t_graph *g,
2919                            gmx_enerdata_t *enerd, t_nrnb *nrnb,
2920                            real *lambda, real *dvdl,
2921                            const t_mdatoms *md,t_fcdata *fcd,
2922                            int *global_atom_index, gmx_bool bPrintSepPot)
2923 {
2924     int ind,nat1,nbonds,efptFTYPE,nbonds_np;
2925     real v=0;
2926     t_iatom *iatoms;
2927
2928     if (IS_RESTRAINT_TYPE(ftype))
2929     {
2930         efptFTYPE = efptRESTRAINT;
2931     }
2932     else
2933     {
2934         efptFTYPE = efptBONDED;
2935     }
2936
2937     if (ftype<F_GB12 || ftype>F_GB14)
2938     {
2939         if (interaction_function[ftype].flags & IF_BOND &&
2940             !(ftype == F_CONNBONDS || ftype == F_POSRES))
2941         {
2942             ind  = interaction_function[ftype].nrnb_ind;
2943             nat1 = interaction_function[ftype].nratoms+1;
2944             nbonds_np = idef->il[ftype].nr_nonperturbed;
2945             nbonds    = idef->il[ftype].nr - nbonds_np;
2946             iatoms    = idef->il[ftype].iatoms + nbonds_np;
2947             if (nbonds > 0)
2948             {
2949                 if (ftype < F_LJ14 || ftype > F_LJC_PAIRS_NB)
2950                 {
2951                     if(ftype==F_CMAP)
2952                     {
2953                         v = cmap_dihs(nbonds,iatoms,
2954                                       idef->iparams,&idef->cmap_grid,
2955                                       (const rvec*)x,f,fr->fshift,
2956                                       pbc,g,lambda[efptFTYPE],&(dvdl[efptFTYPE]),md,fcd,
2957                                       global_atom_index);
2958                     }
2959                     else
2960                     {
2961                         v =         interaction_function[ftype].ifunc(nbonds,iatoms,
2962                                                                   idef->iparams,
2963                                                                   (const rvec*)x,f,fr->fshift,
2964                                                                   pbc,g,lambda[efptFTYPE],&dvdl[efptFTYPE],
2965                                                                   md,fcd,global_atom_index);
2966                     }
2967                 }
2968                 else
2969                 {
2970                     v = do_listed_vdw_q(ftype,nbonds,iatoms,
2971                                         idef->iparams,
2972                                         (const rvec*)x,f,fr->fshift,
2973                                         pbc,g,lambda,dvdl,
2974                                         md,fr,&enerd->grpp,global_atom_index);
2975                 }
2976                 if (ind != -1)
2977                 {
2978                     inc_nrnb(nrnb,ind,nbonds/nat1);
2979                 }
2980             }
2981         }
2982     }
2983     return v;
2984 }
2985
2986 void calc_bonds(FILE *fplog,const gmx_multisim_t *ms,
2987                 const t_idef *idef,
2988                 rvec x[],history_t *hist,
2989                 rvec f[],t_forcerec *fr,
2990                 const t_pbc *pbc,const t_graph *g,
2991                 gmx_enerdata_t *enerd,t_nrnb *nrnb,
2992                 real *lambda,
2993                 const t_mdatoms *md,
2994                 t_fcdata *fcd,int *global_atom_index,
2995                 t_atomtypes *atype, gmx_genborn_t *born,
2996                 gmx_bool bPrintSepPot,gmx_large_int_t step)
2997 {
2998     int    i,ftype,nbonds,ind,nat;
2999     real   v,dvdl[efptNR],dvdl_dum[efptNR]; /* The dummy array is to have a place to store the dhdl at other values
3000                                                of lambda, which will be thrown away in the end*/
3001     real   *epot;
3002     const  t_pbc *pbc_null;
3003     char   buf[22];
3004
3005     for (i=0;i<efptNR;i++)
3006     {
3007         dvdl[i] = 0.0;
3008     }
3009     if (fr->bMolPBC)
3010     {
3011         pbc_null = pbc;
3012     }
3013     else
3014     {
3015         pbc_null = NULL;
3016     }
3017     if (bPrintSepPot)
3018     {
3019         fprintf(fplog,"Step %s: bonded V and dVdl for this node\n",
3020                 gmx_step_str(step,buf));
3021     }
3022
3023 #ifdef DEBUG
3024     if (g && debug)
3025     {
3026         p_graph(debug,"Bondage is fun",g);
3027     }
3028 #endif
3029
3030     epot = enerd->term;
3031
3032     /* Do pre force calculation stuff which might require communication */
3033     if (idef->il[F_ORIRES].nr) {
3034         epot[F_ORIRESDEV] = calc_orires_dev(ms,idef->il[F_ORIRES].nr,
3035                                             idef->il[F_ORIRES].iatoms,
3036                                             idef->iparams,md,(const rvec*)x,
3037                                             pbc_null,fcd,hist);
3038     }
3039     if (idef->il[F_DISRES].nr) {
3040         calc_disres_R_6(ms,idef->il[F_DISRES].nr,
3041                         idef->il[F_DISRES].iatoms,
3042                         idef->iparams,(const rvec*)x,pbc_null,
3043                         fcd,hist);
3044     }
3045
3046     /* Loop over all bonded force types to calculate the bonded forces */
3047     for(ftype=0; (ftype<F_NRE); ftype++)
3048     {
3049         v = calc_one_bond(fplog,ftype,idef,x,
3050                           f,fr,pbc_null,g,enerd,nrnb,lambda,dvdl,
3051                           md,fcd,global_atom_index,bPrintSepPot);
3052         epot[ftype]        += v;
3053     }
3054     /* Copy the sum of violations for the distance restraints from fcd */
3055     if (fcd)
3056     {
3057         epot[F_DISRESVIOL] = fcd->disres.sumviol;
3058     }
3059 }
3060
3061 void calc_bonds_lambda(FILE *fplog,
3062                        const t_idef *idef,
3063                        rvec x[],
3064                        t_forcerec *fr,
3065                        const t_pbc *pbc,const t_graph *g,
3066                        gmx_enerdata_t *enerd,t_nrnb *nrnb,
3067                        real *lambda,
3068                        const t_mdatoms *md,
3069                        t_fcdata *fcd,
3070                        int *global_atom_index)
3071 {
3072     int    i,ftype,nbonds_np,nbonds,ind,nat;
3073     real   v,dr,dr2,*epot;
3074     real   dvdl_dum[efptNR];
3075     rvec   *f,*fshift_orig;
3076     const  t_pbc *pbc_null;
3077     t_iatom *iatom_fe;
3078
3079     if (fr->bMolPBC)
3080     {
3081         pbc_null = pbc;
3082     }
3083     else
3084     {
3085         pbc_null = NULL;
3086     }
3087
3088     epot = enerd->term;
3089
3090     snew(f,fr->natoms_force);
3091     /* We want to preserve the fshift array in forcerec */
3092     fshift_orig = fr->fshift;
3093     snew(fr->fshift,SHIFTS);
3094
3095     /* Loop over all bonded force types to calculate the bonded forces */
3096     for(ftype=0; (ftype<F_NRE); ftype++)
3097     {
3098         v = calc_one_bond_foreign(fplog,ftype,idef,x,
3099                                   f,fr,pbc_null,g,enerd,nrnb,lambda,dvdl_dum,
3100                                   md,fcd,global_atom_index,FALSE);
3101         epot[ftype] += v;
3102     }
3103
3104     sfree(fr->fshift);
3105     fr->fshift = fshift_orig;
3106     sfree(f);
3107 }