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