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                 }
1445             }
1446             else
1447             {
1448                 ref      = L1*posA + lambda*posB;
1449                 rdist[m] = 0;
1450                 dpdl[m]  = posB - posA;
1451             }
1452
1453             /* We do pbc_dx with ref+rdist,
1454              * since with only ref we can be up to half a box vector wrong.
1455              */
1456             pos[m] = ref + rdist[m];
1457         }
1458
1459         if (pbc)
1460         {
1461             pbc_dx(pbc,x[ai],pos,dx);
1462         }
1463         else
1464         {
1465             rvec_sub(x[ai],pos,dx);
1466         }
1467
1468         for (m=0; (m<DIM); m++)
1469         {
1470             kk          = L1*pr->posres.fcA[m] + lambda*pr->posres.fcB[m];
1471             fm          = -kk*dx[m];
1472             f[ai][m]   += fm;
1473             vtot       += 0.5*kk*dx[m]*dx[m];
1474             *dvdlambda +=
1475                 0.5*(pr->posres.fcB[m] - pr->posres.fcA[m])*dx[m]*dx[m]
1476                 -fm*dpdl[m];
1477
1478             /* Here we correct for the pbc_dx which included rdist */
1479             vir_diag[m] -= 0.5*(dx[m] + rdist[m])*fm;
1480         }
1481     }
1482
1483     return vtot;
1484 }
1485
1486 static real low_angres(int nbonds,
1487                        const t_iatom forceatoms[],const t_iparams forceparams[],
1488                        const rvec x[],rvec f[],rvec fshift[],
1489                        const t_pbc *pbc,const t_graph *g,
1490                        real lambda,real *dvdlambda,
1491                        gmx_bool bZAxis)
1492 {
1493   int  i,m,type,ai,aj,ak,al;
1494   int  t1,t2;
1495   real phi,cos_phi,cos_phi2,vid,vtot,dVdphi;
1496   rvec r_ij,r_kl,f_i,f_k={0,0,0};
1497   real st,sth,nrij2,nrkl2,c,cij,ckl;
1498
1499   ivec dt;
1500   t2 = 0; /* avoid warning with gcc-3.3. It is never used uninitialized */
1501
1502   vtot = 0.0;
1503   ak=al=0; /* to avoid warnings */
1504   for(i=0; i<nbonds; ) {
1505     type = forceatoms[i++];
1506     ai   = forceatoms[i++];
1507     aj   = forceatoms[i++];
1508     t1   = pbc_rvec_sub(pbc,x[aj],x[ai],r_ij);                  /*  3           */
1509     if (!bZAxis) {
1510       ak   = forceatoms[i++];
1511       al   = forceatoms[i++];
1512       t2   = pbc_rvec_sub(pbc,x[al],x[ak],r_kl);           /*  3                */
1513     } else {
1514       r_kl[XX] = 0;
1515       r_kl[YY] = 0;
1516       r_kl[ZZ] = 1;
1517     }
1518
1519     cos_phi = cos_angle(r_ij,r_kl);             /* 25           */
1520     phi     = acos(cos_phi);                    /* 10           */
1521
1522     *dvdlambda += dopdihs_min(forceparams[type].pdihs.cpA,
1523                               forceparams[type].pdihs.cpB,
1524                               forceparams[type].pdihs.phiA,
1525                               forceparams[type].pdihs.phiB,
1526                               forceparams[type].pdihs.mult,
1527                               phi,lambda,&vid,&dVdphi); /*  40  */
1528
1529     vtot += vid;
1530
1531     cos_phi2 = sqr(cos_phi);                    /*   1          */
1532     if (cos_phi2 < 1) {
1533       st  = -dVdphi*gmx_invsqrt(1 - cos_phi2);      /*  12              */
1534       sth = st*cos_phi;                         /*   1          */
1535       nrij2 = iprod(r_ij,r_ij);                 /*   5          */
1536       nrkl2 = iprod(r_kl,r_kl);                 /*   5          */
1537
1538       c   = st*gmx_invsqrt(nrij2*nrkl2);                /*  11          */
1539       cij = sth/nrij2;                          /*  10          */
1540       ckl = sth/nrkl2;                          /*  10          */
1541
1542       for (m=0; m<DIM; m++) {                   /*  18+18       */
1543         f_i[m] = (c*r_kl[m]-cij*r_ij[m]);
1544         f[ai][m] += f_i[m];
1545         f[aj][m] -= f_i[m];
1546         if (!bZAxis) {
1547           f_k[m] = (c*r_ij[m]-ckl*r_kl[m]);
1548           f[ak][m] += f_k[m];
1549           f[al][m] -= f_k[m];
1550         }
1551       }
1552
1553       if (g) {
1554         ivec_sub(SHIFT_IVEC(g,ai),SHIFT_IVEC(g,aj),dt);
1555         t1=IVEC2IS(dt);
1556       }
1557       rvec_inc(fshift[t1],f_i);
1558       rvec_dec(fshift[CENTRAL],f_i);
1559       if (!bZAxis) {
1560         if (g) {
1561           ivec_sub(SHIFT_IVEC(g,ak),SHIFT_IVEC(g,al),dt);
1562           t2=IVEC2IS(dt);
1563         }
1564         rvec_inc(fshift[t2],f_k);
1565         rvec_dec(fshift[CENTRAL],f_k);
1566       }
1567     }
1568   }
1569
1570   return vtot;  /*  184 / 157 (bZAxis)  total  */
1571 }
1572
1573 real angres(int nbonds,
1574             const t_iatom forceatoms[],const t_iparams forceparams[],
1575             const rvec x[],rvec f[],rvec fshift[],
1576             const t_pbc *pbc,const t_graph *g,
1577             real lambda,real *dvdlambda,
1578             const t_mdatoms *md,t_fcdata *fcd,
1579             int *global_atom_index)
1580 {
1581   return low_angres(nbonds,forceatoms,forceparams,x,f,fshift,pbc,g,
1582                     lambda,dvdlambda,FALSE);
1583 }
1584
1585 real angresz(int nbonds,
1586              const t_iatom forceatoms[],const t_iparams forceparams[],
1587              const rvec x[],rvec f[],rvec fshift[],
1588              const t_pbc *pbc,const t_graph *g,
1589              real lambda,real *dvdlambda,
1590              const t_mdatoms *md,t_fcdata *fcd,
1591              int *global_atom_index)
1592 {
1593   return low_angres(nbonds,forceatoms,forceparams,x,f,fshift,pbc,g,
1594                     lambda,dvdlambda,TRUE);
1595 }
1596
1597
1598 real unimplemented(int nbonds,
1599                    const t_iatom forceatoms[],const t_iparams forceparams[],
1600                    const rvec x[],rvec f[],rvec fshift[],
1601                    const t_pbc *pbc,const t_graph *g,
1602                    real lambda,real *dvdlambda,
1603                    const t_mdatoms *md,t_fcdata *fcd,
1604                    int *global_atom_index)
1605 {
1606   gmx_impl("*** you are using a not implemented function");
1607
1608   return 0.0; /* To make the compiler happy */
1609 }
1610
1611 real rbdihs(int nbonds,
1612             const t_iatom forceatoms[],const t_iparams forceparams[],
1613             const rvec x[],rvec f[],rvec fshift[],
1614             const t_pbc *pbc,const t_graph *g,
1615             real lambda,real *dvdlambda,
1616             const t_mdatoms *md,t_fcdata *fcd,
1617             int *global_atom_index)
1618 {
1619   const real c0=0.0,c1=1.0,c2=2.0,c3=3.0,c4=4.0,c5=5.0;
1620   int  type,ai,aj,ak,al,i,j;
1621   int  t1,t2,t3;
1622   rvec r_ij,r_kj,r_kl,m,n;
1623   real parmA[NR_RBDIHS];
1624   real parmB[NR_RBDIHS];
1625   real parm[NR_RBDIHS];
1626   real cos_phi,phi,rbp,rbpBA;
1627   real v,sign,ddphi,sin_phi;
1628   real cosfac,vtot;
1629   real L1   = 1.0-lambda;
1630   real dvdl=0;
1631
1632   vtot = 0.0;
1633   for(i=0; (i<nbonds); ) {
1634     type = forceatoms[i++];
1635     ai   = forceatoms[i++];
1636     aj   = forceatoms[i++];
1637     ak   = forceatoms[i++];
1638     al   = forceatoms[i++];
1639
1640       phi=dih_angle(x[ai],x[aj],x[ak],x[al],pbc,r_ij,r_kj,r_kl,m,n,
1641                     &sign,&t1,&t2,&t3);                 /*  84          */
1642
1643     /* Change to polymer convention */
1644     if (phi < c0)
1645       phi += M_PI;
1646     else
1647       phi -= M_PI;                      /*   1          */
1648
1649     cos_phi = cos(phi);
1650     /* Beware of accuracy loss, cannot use 1-sqrt(cos^2) ! */
1651     sin_phi = sin(phi);
1652
1653     for(j=0; (j<NR_RBDIHS); j++) {
1654       parmA[j] = forceparams[type].rbdihs.rbcA[j];
1655       parmB[j] = forceparams[type].rbdihs.rbcB[j];
1656       parm[j]  = L1*parmA[j]+lambda*parmB[j];
1657     }
1658     /* Calculate cosine powers */
1659     /* Calculate the energy */
1660     /* Calculate the derivative */
1661
1662     v       = parm[0];
1663     dvdl   += (parmB[0]-parmA[0]);
1664     ddphi   = c0;
1665     cosfac  = c1;
1666
1667     rbp     = parm[1];
1668     rbpBA   = parmB[1]-parmA[1];
1669     ddphi  += rbp*cosfac;
1670     cosfac *= cos_phi;
1671     v      += cosfac*rbp;
1672     dvdl   += cosfac*rbpBA;
1673     rbp     = parm[2];
1674     rbpBA   = parmB[2]-parmA[2];
1675     ddphi  += c2*rbp*cosfac;
1676     cosfac *= cos_phi;
1677     v      += cosfac*rbp;
1678     dvdl   += cosfac*rbpBA;
1679     rbp     = parm[3];
1680     rbpBA   = parmB[3]-parmA[3];
1681     ddphi  += c3*rbp*cosfac;
1682     cosfac *= cos_phi;
1683     v      += cosfac*rbp;
1684     dvdl   += cosfac*rbpBA;
1685     rbp     = parm[4];
1686     rbpBA   = parmB[4]-parmA[4];
1687     ddphi  += c4*rbp*cosfac;
1688     cosfac *= cos_phi;
1689     v      += cosfac*rbp;
1690     dvdl   += cosfac*rbpBA;
1691     rbp     = parm[5];
1692     rbpBA   = parmB[5]-parmA[5];
1693     ddphi  += c5*rbp*cosfac;
1694     cosfac *= cos_phi;
1695     v      += cosfac*rbp;
1696     dvdl   += cosfac*rbpBA;
1697
1698     ddphi = -ddphi*sin_phi;                             /*  11          */
1699
1700     do_dih_fup(ai,aj,ak,al,ddphi,r_ij,r_kj,r_kl,m,n,
1701                f,fshift,pbc,g,x,t1,t2,t3);              /* 112          */
1702     vtot += v;
1703   }
1704   *dvdlambda += dvdl;
1705
1706   return vtot;
1707 }
1708
1709 int cmap_setup_grid_index(int ip, int grid_spacing, int *ipm1, int *ipp1, int *ipp2)
1710 {
1711         int im1, ip1, ip2;
1712
1713         if(ip<0)
1714         {
1715                 ip = ip + grid_spacing - 1;
1716         }
1717         else if(ip > grid_spacing)
1718         {
1719                 ip = ip - grid_spacing - 1;
1720         }
1721
1722         im1 = ip - 1;
1723         ip1 = ip + 1;
1724         ip2 = ip + 2;
1725
1726         if(ip == 0)
1727         {
1728                 im1 = grid_spacing - 1;
1729         }
1730         else if(ip == grid_spacing-2)
1731         {
1732                 ip2 = 0;
1733         }
1734         else if(ip == grid_spacing-1)
1735         {
1736                 ip1 = 0;
1737                 ip2 = 1;
1738         }
1739
1740         *ipm1 = im1;
1741         *ipp1 = ip1;
1742         *ipp2 = ip2;
1743
1744         return ip;
1745
1746 }
1747
1748 real cmap_dihs(int nbonds,
1749                            const t_iatom forceatoms[],const t_iparams forceparams[],
1750                const gmx_cmap_t *cmap_grid,
1751                            const rvec x[],rvec f[],rvec fshift[],
1752                            const t_pbc *pbc,const t_graph *g,
1753                            real lambda,real *dvdlambda,
1754                            const t_mdatoms *md,t_fcdata *fcd,
1755                            int *global_atom_index)
1756 {
1757         int i,j,k,n,idx;
1758         int ai,aj,ak,al,am;
1759         int a1i,a1j,a1k,a1l,a2i,a2j,a2k,a2l;
1760         int type,cmapA;
1761         int t11,t21,t31,t12,t22,t32;
1762         int iphi1,ip1m1,ip1p1,ip1p2;
1763         int iphi2,ip2m1,ip2p1,ip2p2;
1764         int l1,l2,l3,l4;
1765         int pos1,pos2,pos3,pos4,tmp;
1766
1767         real ty[4],ty1[4],ty2[4],ty12[4],tc[16],tx[16];
1768         real phi1,psi1,cos_phi1,sin_phi1,sign1,xphi1;
1769         real phi2,psi2,cos_phi2,sin_phi2,sign2,xphi2;
1770         real dx,xx,tt,tu,e,df1,df2,ddf1,ddf2,ddf12,vtot;
1771         real ra21,rb21,rg21,rg1,rgr1,ra2r1,rb2r1,rabr1;
1772         real ra22,rb22,rg22,rg2,rgr2,ra2r2,rb2r2,rabr2;
1773         real fg1,hg1,fga1,hgb1,gaa1,gbb1;
1774         real fg2,hg2,fga2,hgb2,gaa2,gbb2;
1775         real fac;
1776
1777         rvec r1_ij, r1_kj, r1_kl,m1,n1;
1778         rvec r2_ij, r2_kj, r2_kl,m2,n2;
1779         rvec f1_i,f1_j,f1_k,f1_l;
1780         rvec f2_i,f2_j,f2_k,f2_l;
1781         rvec a1,b1,a2,b2;
1782         rvec f1,g1,h1,f2,g2,h2;
1783         rvec dtf1,dtg1,dth1,dtf2,dtg2,dth2;
1784         ivec jt1,dt1_ij,dt1_kj,dt1_lj;
1785         ivec jt2,dt2_ij,dt2_kj,dt2_lj;
1786
1787     const real *cmapd;
1788
1789         int loop_index[4][4] = {
1790                 {0,4,8,12},
1791                 {1,5,9,13},
1792                 {2,6,10,14},
1793                 {3,7,11,15}
1794         };
1795
1796         /* Total CMAP energy */
1797         vtot = 0;
1798
1799         for(n=0;n<nbonds; )
1800         {
1801                 /* Five atoms are involved in the two torsions */
1802                 type   = forceatoms[n++];
1803                 ai     = forceatoms[n++];
1804                 aj     = forceatoms[n++];
1805                 ak     = forceatoms[n++];
1806                 al     = forceatoms[n++];
1807                 am     = forceatoms[n++];
1808
1809                 /* Which CMAP type is this */
1810                 cmapA = forceparams[type].cmap.cmapA;
1811         cmapd = cmap_grid->cmapdata[cmapA].cmap;
1812
1813                 /* First torsion */
1814                 a1i   = ai;
1815                 a1j   = aj;
1816                 a1k   = ak;
1817                 a1l   = al;
1818
1819                 phi1  = dih_angle(x[a1i], x[a1j], x[a1k], x[a1l], pbc, r1_ij, r1_kj, r1_kl, m1, n1,
1820                                                    &sign1, &t11, &t21, &t31); /* 84 */
1821
1822         cos_phi1 = cos(phi1);
1823
1824                 a1[0] = r1_ij[1]*r1_kj[2]-r1_ij[2]*r1_kj[1];
1825                 a1[1] = r1_ij[2]*r1_kj[0]-r1_ij[0]*r1_kj[2];
1826                 a1[2] = r1_ij[0]*r1_kj[1]-r1_ij[1]*r1_kj[0]; /* 9 */
1827
1828                 b1[0] = r1_kl[1]*r1_kj[2]-r1_kl[2]*r1_kj[1];
1829                 b1[1] = r1_kl[2]*r1_kj[0]-r1_kl[0]*r1_kj[2];
1830                 b1[2] = r1_kl[0]*r1_kj[1]-r1_kl[1]*r1_kj[0]; /* 9 */
1831
1832                 tmp = pbc_rvec_sub(pbc,x[a1l],x[a1k],h1);
1833
1834                 ra21  = iprod(a1,a1);       /* 5 */
1835                 rb21  = iprod(b1,b1);       /* 5 */
1836                 rg21  = iprod(r1_kj,r1_kj); /* 5 */
1837                 rg1   = sqrt(rg21);
1838
1839                 rgr1  = 1.0/rg1;
1840                 ra2r1 = 1.0/ra21;
1841                 rb2r1 = 1.0/rb21;
1842                 rabr1 = sqrt(ra2r1*rb2r1);
1843
1844                 sin_phi1 = rg1 * rabr1 * iprod(a1,h1) * (-1);
1845
1846                 if(cos_phi1 < -0.5 || cos_phi1 > 0.5)
1847                 {
1848                         phi1 = asin(sin_phi1);
1849
1850                         if(cos_phi1 < 0)
1851                         {
1852                                 if(phi1 > 0)
1853                                 {
1854                                         phi1 = M_PI - phi1;
1855                                 }
1856                                 else
1857                                 {
1858                                         phi1 = -M_PI - phi1;
1859                                 }
1860                         }
1861                 }
1862                 else
1863                 {
1864                         phi1 = acos(cos_phi1);
1865
1866                         if(sin_phi1 < 0)
1867                         {
1868                                 phi1 = -phi1;
1869                         }
1870                 }
1871
1872                 xphi1 = phi1 + M_PI; /* 1 */
1873
1874                 /* Second torsion */
1875                 a2i   = aj;
1876                 a2j   = ak;
1877                 a2k   = al;
1878                 a2l   = am;
1879
1880                 phi2  = dih_angle(x[a2i], x[a2j], x[a2k], x[a2l], pbc, r2_ij, r2_kj, r2_kl, m2, n2,
1881                                                   &sign2, &t12, &t22, &t32); /* 84 */
1882
1883         cos_phi2 = cos(phi2);
1884
1885                 a2[0] = r2_ij[1]*r2_kj[2]-r2_ij[2]*r2_kj[1];
1886                 a2[1] = r2_ij[2]*r2_kj[0]-r2_ij[0]*r2_kj[2];
1887                 a2[2] = r2_ij[0]*r2_kj[1]-r2_ij[1]*r2_kj[0]; /* 9 */
1888
1889                 b2[0] = r2_kl[1]*r2_kj[2]-r2_kl[2]*r2_kj[1];
1890                 b2[1] = r2_kl[2]*r2_kj[0]-r2_kl[0]*r2_kj[2];
1891                 b2[2] = r2_kl[0]*r2_kj[1]-r2_kl[1]*r2_kj[0]; /* 9 */
1892
1893                 tmp = pbc_rvec_sub(pbc,x[a2l],x[a2k],h2);
1894
1895                 ra22  = iprod(a2,a2);         /* 5 */
1896                 rb22  = iprod(b2,b2);         /* 5 */
1897                 rg22  = iprod(r2_kj,r2_kj);   /* 5 */
1898                 rg2   = sqrt(rg22);
1899
1900                 rgr2  = 1.0/rg2;
1901                 ra2r2 = 1.0/ra22;
1902                 rb2r2 = 1.0/rb22;
1903                 rabr2 = sqrt(ra2r2*rb2r2);
1904
1905                 sin_phi2 = rg2 * rabr2 * iprod(a2,h2) * (-1);
1906
1907                 if(cos_phi2 < -0.5 || cos_phi2 > 0.5)
1908                 {
1909                         phi2 = asin(sin_phi2);
1910
1911                         if(cos_phi2 < 0)
1912                         {
1913                                 if(phi2 > 0)
1914                                 {
1915                                         phi2 = M_PI - phi2;
1916                                 }
1917                                 else
1918                                 {
1919                                         phi2 = -M_PI - phi2;
1920                                 }
1921                         }
1922                 }
1923                 else
1924                 {
1925                         phi2 = acos(cos_phi2);
1926
1927                         if(sin_phi2 < 0)
1928                         {
1929                                 phi2 = -phi2;
1930                         }
1931                 }
1932
1933                 xphi2 = phi2 + M_PI; /* 1 */
1934
1935                 /* Range mangling */
1936                 if(xphi1<0)
1937                 {
1938                         xphi1 = xphi1 + 2*M_PI;
1939                 }
1940                 else if(xphi1>=2*M_PI)
1941                 {
1942                         xphi1 = xphi1 - 2*M_PI;
1943                 }
1944
1945                 if(xphi2<0)
1946                 {
1947                         xphi2 = xphi2 + 2*M_PI;
1948                 }
1949                 else if(xphi2>=2*M_PI)
1950                 {
1951                         xphi2 = xphi2 - 2*M_PI;
1952                 }
1953
1954                 /* Number of grid points */
1955                 dx = 2*M_PI / cmap_grid->grid_spacing;
1956
1957                 /* Where on the grid are we */
1958                 iphi1 = (int)(xphi1/dx);
1959                 iphi2 = (int)(xphi2/dx);
1960
1961                 iphi1 = cmap_setup_grid_index(iphi1, cmap_grid->grid_spacing, &ip1m1,&ip1p1,&ip1p2);
1962                 iphi2 = cmap_setup_grid_index(iphi2, cmap_grid->grid_spacing, &ip2m1,&ip2p1,&ip2p2);
1963
1964                 pos1    = iphi1*cmap_grid->grid_spacing+iphi2;
1965                 pos2    = ip1p1*cmap_grid->grid_spacing+iphi2;
1966                 pos3    = ip1p1*cmap_grid->grid_spacing+ip2p1;
1967                 pos4    = iphi1*cmap_grid->grid_spacing+ip2p1;
1968
1969                 ty[0]   = cmapd[pos1*4];
1970                 ty[1]   = cmapd[pos2*4];
1971                 ty[2]   = cmapd[pos3*4];
1972                 ty[3]   = cmapd[pos4*4];
1973
1974                 ty1[0]   = cmapd[pos1*4+1];
1975                 ty1[1]   = cmapd[pos2*4+1];
1976                 ty1[2]   = cmapd[pos3*4+1];
1977                 ty1[3]   = cmapd[pos4*4+1];
1978
1979                 ty2[0]   = cmapd[pos1*4+2];
1980                 ty2[1]   = cmapd[pos2*4+2];
1981                 ty2[2]   = cmapd[pos3*4+2];
1982                 ty2[3]   = cmapd[pos4*4+2];
1983
1984                 ty12[0]   = cmapd[pos1*4+3];
1985                 ty12[1]   = cmapd[pos2*4+3];
1986                 ty12[2]   = cmapd[pos3*4+3];
1987                 ty12[3]   = cmapd[pos4*4+3];
1988
1989                 /* Switch to degrees */
1990                 dx = 360.0 / cmap_grid->grid_spacing;
1991                 xphi1 = xphi1 * RAD2DEG;
1992                 xphi2 = xphi2 * RAD2DEG;
1993
1994                 for(i=0;i<4;i++) /* 16 */
1995                 {
1996                         tx[i] = ty[i];
1997                         tx[i+4] = ty1[i]*dx;
1998                         tx[i+8] = ty2[i]*dx;
1999                         tx[i+12] = ty12[i]*dx*dx;
2000                 }
2001
2002                 idx=0;
2003                 for(i=0;i<4;i++) /* 1056 */
2004                 {
2005                         for(j=0;j<4;j++)
2006                         {
2007                                 xx = 0;
2008                                 for(k=0;k<16;k++)
2009                                 {
2010                                         xx = xx + cmap_coeff_matrix[k*16+idx]*tx[k];
2011                                 }
2012
2013                                 idx++;
2014                                 tc[i*4+j]=xx;
2015                         }
2016                 }
2017
2018                 tt    = (xphi1-iphi1*dx)/dx;
2019                 tu    = (xphi2-iphi2*dx)/dx;
2020
2021                 e     = 0;
2022                 df1   = 0;
2023                 df2   = 0;
2024                 ddf1  = 0;
2025                 ddf2  = 0;
2026                 ddf12 = 0;
2027
2028                 for(i=3;i>=0;i--)
2029                 {
2030                         l1 = loop_index[i][3];
2031                         l2 = loop_index[i][2];
2032                         l3 = loop_index[i][1];
2033
2034                         e     = tt * e    + ((tc[i*4+3]*tu+tc[i*4+2])*tu + tc[i*4+1])*tu+tc[i*4];
2035                         df1   = tu * df1  + (3.0*tc[l1]*tt+2.0*tc[l2])*tt+tc[l3];
2036                         df2   = tt * df2  + (3.0*tc[i*4+3]*tu+2.0*tc[i*4+2])*tu+tc[i*4+1];
2037                         ddf1  = tu * ddf1 + 2.0*3.0*tc[l1]*tt+2.0*tc[l2];
2038                         ddf2  = tt * ddf2 + 2.0*3.0*tc[4*i+3]*tu+2.0*tc[4*i+2];
2039                 }
2040
2041                 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) +
2042                 3.0*tu*tu*(tc[7]+2.0*tc[11]*tt+3.0*tc[15]*tt*tt);
2043
2044                 fac     = RAD2DEG/dx;
2045                 df1     = df1   * fac;
2046                 df2     = df2   * fac;
2047                 ddf1    = ddf1  * fac * fac;
2048                 ddf2    = ddf2  * fac * fac;
2049                 ddf12   = ddf12 * fac * fac;
2050
2051                 /* CMAP energy */
2052                 vtot += e;
2053
2054                 /* Do forces - first torsion */
2055                 fg1       = iprod(r1_ij,r1_kj);
2056                 hg1       = iprod(r1_kl,r1_kj);
2057                 fga1      = fg1*ra2r1*rgr1;
2058                 hgb1      = hg1*rb2r1*rgr1;
2059                 gaa1      = -ra2r1*rg1;
2060                 gbb1      = rb2r1*rg1;
2061
2062                 for(i=0;i<DIM;i++)
2063                 {
2064                         dtf1[i]   = gaa1 * a1[i];
2065                         dtg1[i]   = fga1 * a1[i] - hgb1 * b1[i];
2066                         dth1[i]   = gbb1 * b1[i];
2067
2068                         f1[i]     = df1  * dtf1[i];
2069                         g1[i]     = df1  * dtg1[i];
2070                         h1[i]     = df1  * dth1[i];
2071
2072                         f1_i[i]   =  f1[i];
2073                         f1_j[i]   = -f1[i] - g1[i];
2074                         f1_k[i]   =  h1[i] + g1[i];
2075                         f1_l[i]   = -h1[i];
2076
2077                         f[a1i][i] = f[a1i][i] + f1_i[i];
2078                         f[a1j][i] = f[a1j][i] + f1_j[i]; /* - f1[i] - g1[i] */
2079                         f[a1k][i] = f[a1k][i] + f1_k[i]; /* h1[i] + g1[i] */
2080                         f[a1l][i] = f[a1l][i] + f1_l[i]; /* h1[i] */
2081                 }
2082
2083                 /* Do forces - second torsion */
2084                 fg2       = iprod(r2_ij,r2_kj);
2085                 hg2       = iprod(r2_kl,r2_kj);
2086                 fga2      = fg2*ra2r2*rgr2;
2087                 hgb2      = hg2*rb2r2*rgr2;
2088                 gaa2      = -ra2r2*rg2;
2089                 gbb2      = rb2r2*rg2;
2090
2091                 for(i=0;i<DIM;i++)
2092                 {
2093                         dtf2[i]   = gaa2 * a2[i];
2094                         dtg2[i]   = fga2 * a2[i] - hgb2 * b2[i];
2095                         dth2[i]   = gbb2 * b2[i];
2096
2097                         f2[i]     = df2  * dtf2[i];
2098                         g2[i]     = df2  * dtg2[i];
2099                         h2[i]     = df2  * dth2[i];
2100
2101                         f2_i[i]   =  f2[i];
2102                         f2_j[i]   = -f2[i] - g2[i];
2103                         f2_k[i]   =  h2[i] + g2[i];
2104                         f2_l[i]   = -h2[i];
2105
2106                         f[a2i][i] = f[a2i][i] + f2_i[i]; /* f2[i] */
2107                         f[a2j][i] = f[a2j][i] + f2_j[i]; /* - f2[i] - g2[i] */
2108                         f[a2k][i] = f[a2k][i] + f2_k[i]; /* h2[i] + g2[i] */
2109                         f[a2l][i] = f[a2l][i] + f2_l[i]; /* - h2[i] */
2110                 }
2111
2112                 /* Shift forces */
2113                 if(g)
2114                 {
2115                         copy_ivec(SHIFT_IVEC(g,a1j), jt1);
2116                         ivec_sub(SHIFT_IVEC(g,a1i),  jt1,dt1_ij);
2117                         ivec_sub(SHIFT_IVEC(g,a1k),  jt1,dt1_kj);
2118                         ivec_sub(SHIFT_IVEC(g,a1l),  jt1,dt1_lj);
2119                         t11 = IVEC2IS(dt1_ij);
2120                         t21 = IVEC2IS(dt1_kj);
2121                         t31 = IVEC2IS(dt1_lj);
2122
2123                         copy_ivec(SHIFT_IVEC(g,a2j), jt2);
2124                         ivec_sub(SHIFT_IVEC(g,a2i),  jt2,dt2_ij);
2125                         ivec_sub(SHIFT_IVEC(g,a2k),  jt2,dt2_kj);
2126                         ivec_sub(SHIFT_IVEC(g,a2l),  jt2,dt2_lj);
2127                         t12 = IVEC2IS(dt2_ij);
2128                         t22 = IVEC2IS(dt2_kj);
2129                         t32 = IVEC2IS(dt2_lj);
2130                 }
2131                 else if(pbc)
2132                 {
2133                         t31 = pbc_rvec_sub(pbc,x[a1l],x[a1j],h1);
2134                         t32 = pbc_rvec_sub(pbc,x[a2l],x[a2j],h2);
2135                 }
2136                 else
2137                 {
2138                         t31 = CENTRAL;
2139                         t32 = CENTRAL;
2140                 }
2141
2142                 rvec_inc(fshift[t11],f1_i);
2143                 rvec_inc(fshift[CENTRAL],f1_j);
2144                 rvec_inc(fshift[t21],f1_k);
2145                 rvec_inc(fshift[t31],f1_l);
2146
2147                 rvec_inc(fshift[t21],f2_i);
2148                 rvec_inc(fshift[CENTRAL],f2_j);
2149                 rvec_inc(fshift[t22],f2_k);
2150                 rvec_inc(fshift[t32],f2_l);
2151         }
2152         return vtot;
2153 }
2154
2155
2156
2157 /***********************************************************
2158  *
2159  *   G R O M O S  9 6   F U N C T I O N S
2160  *
2161  ***********************************************************/
2162 real g96harmonic(real kA,real kB,real xA,real xB,real x,real lambda,
2163                  real *V,real *F)
2164 {
2165   const real half=0.5;
2166   real  L1,kk,x0,dx,dx2;
2167   real  v,f,dvdl;
2168
2169   L1    = 1.0-lambda;
2170   kk    = L1*kA+lambda*kB;
2171   x0    = L1*xA+lambda*xB;
2172
2173   dx    = x-x0;
2174   dx2   = dx*dx;
2175
2176   f     = -kk*dx;
2177   v     = half*kk*dx2;
2178   dvdl  = half*(kB-kA)*dx2 + (xA-xB)*kk*dx;
2179
2180   *F    = f;
2181   *V    = v;
2182
2183   return dvdl;
2184
2185   /* That was 21 flops */
2186 }
2187
2188 real g96bonds(int nbonds,
2189               const t_iatom forceatoms[],const t_iparams forceparams[],
2190               const rvec x[],rvec f[],rvec fshift[],
2191               const t_pbc *pbc,const t_graph *g,
2192               real lambda,real *dvdlambda,
2193               const t_mdatoms *md,t_fcdata *fcd,
2194               int *global_atom_index)
2195 {
2196   int  i,m,ki,ai,aj,type;
2197   real dr2,fbond,vbond,fij,vtot;
2198   rvec dx;
2199   ivec dt;
2200
2201   vtot = 0.0;
2202   for(i=0; (i<nbonds); ) {
2203     type = forceatoms[i++];
2204     ai   = forceatoms[i++];
2205     aj   = forceatoms[i++];
2206
2207     ki   = pbc_rvec_sub(pbc,x[ai],x[aj],dx);            /*   3          */
2208     dr2  = iprod(dx,dx);                                /*   5          */
2209
2210     *dvdlambda += g96harmonic(forceparams[type].harmonic.krA,
2211                               forceparams[type].harmonic.krB,
2212                               forceparams[type].harmonic.rA,
2213                               forceparams[type].harmonic.rB,
2214                               dr2,lambda,&vbond,&fbond);
2215
2216     vtot  += 0.5*vbond;                             /* 1*/
2217 #ifdef DEBUG
2218     if (debug)
2219       fprintf(debug,"G96-BONDS: dr = %10g  vbond = %10g  fbond = %10g\n",
2220               sqrt(dr2),vbond,fbond);
2221 #endif
2222
2223     if (g) {
2224       ivec_sub(SHIFT_IVEC(g,ai),SHIFT_IVEC(g,aj),dt);
2225       ki=IVEC2IS(dt);
2226     }
2227     for (m=0; (m<DIM); m++) {                   /*  15          */
2228       fij=fbond*dx[m];
2229       f[ai][m]+=fij;
2230       f[aj][m]-=fij;
2231       fshift[ki][m]+=fij;
2232       fshift[CENTRAL][m]-=fij;
2233     }
2234   }                                     /* 44 TOTAL     */
2235   return vtot;
2236 }
2237
2238 real g96bond_angle(const rvec xi,const rvec xj,const rvec xk,const t_pbc *pbc,
2239                    rvec r_ij,rvec r_kj,
2240                    int *t1,int *t2)
2241 /* Return value is the angle between the bonds i-j and j-k */
2242 {
2243   real costh;
2244
2245   *t1 = pbc_rvec_sub(pbc,xi,xj,r_ij);                   /*  3           */
2246   *t2 = pbc_rvec_sub(pbc,xk,xj,r_kj);                   /*  3           */
2247
2248   costh=cos_angle(r_ij,r_kj);           /* 25           */
2249                                         /* 41 TOTAL     */
2250   return costh;
2251 }
2252
2253 real g96angles(int nbonds,
2254                const t_iatom forceatoms[],const t_iparams forceparams[],
2255                const rvec x[],rvec f[],rvec fshift[],
2256                const t_pbc *pbc,const t_graph *g,
2257                real lambda,real *dvdlambda,
2258                const t_mdatoms *md,t_fcdata *fcd,
2259                int *global_atom_index)
2260 {
2261   int  i,ai,aj,ak,type,m,t1,t2;
2262   rvec r_ij,r_kj;
2263   real cos_theta,dVdt,va,vtot;
2264   real rij_1,rij_2,rkj_1,rkj_2,rijrkj_1;
2265   rvec f_i,f_j,f_k;
2266   ivec jt,dt_ij,dt_kj;
2267
2268   vtot = 0.0;
2269   for(i=0; (i<nbonds); ) {
2270     type = forceatoms[i++];
2271     ai   = forceatoms[i++];
2272     aj   = forceatoms[i++];
2273     ak   = forceatoms[i++];
2274
2275     cos_theta  = g96bond_angle(x[ai],x[aj],x[ak],pbc,r_ij,r_kj,&t1,&t2);
2276
2277     *dvdlambda += g96harmonic(forceparams[type].harmonic.krA,
2278                               forceparams[type].harmonic.krB,
2279                               forceparams[type].harmonic.rA,
2280                               forceparams[type].harmonic.rB,
2281                               cos_theta,lambda,&va,&dVdt);
2282     vtot    += va;
2283
2284     rij_1    = gmx_invsqrt(iprod(r_ij,r_ij));
2285     rkj_1    = gmx_invsqrt(iprod(r_kj,r_kj));
2286     rij_2    = rij_1*rij_1;
2287     rkj_2    = rkj_1*rkj_1;
2288     rijrkj_1 = rij_1*rkj_1;                     /* 23 */
2289
2290 #ifdef DEBUG
2291     if (debug)
2292       fprintf(debug,"G96ANGLES: costheta = %10g  vth = %10g  dV/dct = %10g\n",
2293               cos_theta,va,dVdt);
2294 #endif
2295     for (m=0; (m<DIM); m++) {                   /*  42  */
2296       f_i[m]=dVdt*(r_kj[m]*rijrkj_1 - r_ij[m]*rij_2*cos_theta);
2297       f_k[m]=dVdt*(r_ij[m]*rijrkj_1 - r_kj[m]*rkj_2*cos_theta);
2298       f_j[m]=-f_i[m]-f_k[m];
2299       f[ai][m]+=f_i[m];
2300       f[aj][m]+=f_j[m];
2301       f[ak][m]+=f_k[m];
2302     }
2303
2304     if (g) {
2305       copy_ivec(SHIFT_IVEC(g,aj),jt);
2306
2307       ivec_sub(SHIFT_IVEC(g,ai),jt,dt_ij);
2308       ivec_sub(SHIFT_IVEC(g,ak),jt,dt_kj);
2309       t1=IVEC2IS(dt_ij);
2310       t2=IVEC2IS(dt_kj);
2311     }
2312     rvec_inc(fshift[t1],f_i);
2313     rvec_inc(fshift[CENTRAL],f_j);
2314     rvec_inc(fshift[t2],f_k);               /* 9 */
2315     /* 163 TOTAL        */
2316   }
2317   return vtot;
2318 }
2319
2320 real cross_bond_bond(int nbonds,
2321                      const t_iatom forceatoms[],const t_iparams forceparams[],
2322                      const rvec x[],rvec f[],rvec fshift[],
2323                      const t_pbc *pbc,const t_graph *g,
2324                      real lambda,real *dvdlambda,
2325                      const t_mdatoms *md,t_fcdata *fcd,
2326                      int *global_atom_index)
2327 {
2328   /* Potential from Lawrence and Skimmer, Chem. Phys. Lett. 372 (2003)
2329    * pp. 842-847
2330    */
2331   int  i,ai,aj,ak,type,m,t1,t2;
2332   rvec r_ij,r_kj;
2333   real vtot,vrr,s1,s2,r1,r2,r1e,r2e,krr;
2334   rvec f_i,f_j,f_k;
2335   ivec jt,dt_ij,dt_kj;
2336
2337   vtot = 0.0;
2338   for(i=0; (i<nbonds); ) {
2339     type = forceatoms[i++];
2340     ai   = forceatoms[i++];
2341     aj   = forceatoms[i++];
2342     ak   = forceatoms[i++];
2343     r1e  = forceparams[type].cross_bb.r1e;
2344     r2e  = forceparams[type].cross_bb.r2e;
2345     krr  = forceparams[type].cross_bb.krr;
2346
2347     /* Compute distance vectors ... */
2348     t1 = pbc_rvec_sub(pbc,x[ai],x[aj],r_ij);
2349     t2 = pbc_rvec_sub(pbc,x[ak],x[aj],r_kj);
2350
2351     /* ... and their lengths */
2352     r1 = norm(r_ij);
2353     r2 = norm(r_kj);
2354
2355     /* Deviations from ideality */
2356     s1 = r1-r1e;
2357     s2 = r2-r2e;
2358
2359     /* Energy (can be negative!) */
2360     vrr   = krr*s1*s2;
2361     vtot += vrr;
2362
2363     /* Forces */
2364     svmul(-krr*s2/r1,r_ij,f_i);
2365     svmul(-krr*s1/r2,r_kj,f_k);
2366
2367     for (m=0; (m<DIM); m++) {                   /*  12  */
2368       f_j[m]    = -f_i[m] - f_k[m];
2369       f[ai][m] += f_i[m];
2370       f[aj][m] += f_j[m];
2371       f[ak][m] += f_k[m];
2372     }
2373
2374     /* Virial stuff */
2375     if (g) {
2376       copy_ivec(SHIFT_IVEC(g,aj),jt);
2377
2378       ivec_sub(SHIFT_IVEC(g,ai),jt,dt_ij);
2379       ivec_sub(SHIFT_IVEC(g,ak),jt,dt_kj);
2380       t1=IVEC2IS(dt_ij);
2381       t2=IVEC2IS(dt_kj);
2382     }
2383     rvec_inc(fshift[t1],f_i);
2384     rvec_inc(fshift[CENTRAL],f_j);
2385     rvec_inc(fshift[t2],f_k);               /* 9 */
2386     /* 163 TOTAL        */
2387   }
2388   return vtot;
2389 }
2390
2391 real cross_bond_angle(int nbonds,
2392                       const t_iatom forceatoms[],const t_iparams forceparams[],
2393                       const rvec x[],rvec f[],rvec fshift[],
2394                       const t_pbc *pbc,const t_graph *g,
2395                       real lambda,real *dvdlambda,
2396                       const t_mdatoms *md,t_fcdata *fcd,
2397                       int *global_atom_index)
2398 {
2399   /* Potential from Lawrence and Skimmer, Chem. Phys. Lett. 372 (2003)
2400    * pp. 842-847
2401    */
2402   int  i,ai,aj,ak,type,m,t1,t2,t3;
2403   rvec r_ij,r_kj,r_ik;
2404   real vtot,vrt,s1,s2,s3,r1,r2,r3,r1e,r2e,r3e,krt,k1,k2,k3;
2405   rvec f_i,f_j,f_k;
2406   ivec jt,dt_ij,dt_kj;
2407
2408   vtot = 0.0;
2409   for(i=0; (i<nbonds); ) {
2410     type = forceatoms[i++];
2411     ai   = forceatoms[i++];
2412     aj   = forceatoms[i++];
2413     ak   = forceatoms[i++];
2414     r1e  = forceparams[type].cross_ba.r1e;
2415     r2e  = forceparams[type].cross_ba.r2e;
2416     r3e  = forceparams[type].cross_ba.r3e;
2417     krt  = forceparams[type].cross_ba.krt;
2418
2419     /* Compute distance vectors ... */
2420     t1 = pbc_rvec_sub(pbc,x[ai],x[aj],r_ij);
2421     t2 = pbc_rvec_sub(pbc,x[ak],x[aj],r_kj);
2422     t3 = pbc_rvec_sub(pbc,x[ai],x[ak],r_ik);
2423
2424     /* ... and their lengths */
2425     r1 = norm(r_ij);
2426     r2 = norm(r_kj);
2427     r3 = norm(r_ik);
2428
2429     /* Deviations from ideality */
2430     s1 = r1-r1e;
2431     s2 = r2-r2e;
2432     s3 = r3-r3e;
2433
2434     /* Energy (can be negative!) */
2435     vrt   = krt*s3*(s1+s2);
2436     vtot += vrt;
2437
2438     /* Forces */
2439     k1 = -krt*(s3/r1);
2440     k2 = -krt*(s3/r2);
2441     k3 = -krt*(s1+s2)/r3;
2442     for(m=0; (m<DIM); m++) {
2443       f_i[m] = k1*r_ij[m] + k3*r_ik[m];
2444       f_k[m] = k2*r_kj[m] - k3*r_ik[m];
2445       f_j[m] = -f_i[m] - f_k[m];
2446     }
2447
2448     for (m=0; (m<DIM); m++) {                   /*  12  */
2449       f[ai][m] += f_i[m];
2450       f[aj][m] += f_j[m];
2451       f[ak][m] += f_k[m];
2452     }
2453
2454     /* Virial stuff */
2455     if (g) {
2456       copy_ivec(SHIFT_IVEC(g,aj),jt);
2457
2458       ivec_sub(SHIFT_IVEC(g,ai),jt,dt_ij);
2459       ivec_sub(SHIFT_IVEC(g,ak),jt,dt_kj);
2460       t1=IVEC2IS(dt_ij);
2461       t2=IVEC2IS(dt_kj);
2462     }
2463     rvec_inc(fshift[t1],f_i);
2464     rvec_inc(fshift[CENTRAL],f_j);
2465     rvec_inc(fshift[t2],f_k);               /* 9 */
2466     /* 163 TOTAL        */
2467   }
2468   return vtot;
2469 }
2470
2471 static real bonded_tab(const char *type,int table_nr,
2472                        const bondedtable_t *table,real kA,real kB,real r,
2473                        real lambda,real *V,real *F)
2474 {
2475   real k,tabscale,*VFtab,rt,eps,eps2,Yt,Ft,Geps,Heps2,Fp,VV,FF;
2476   int  n0,nnn;
2477   real v,f,dvdl;
2478
2479   k = (1.0 - lambda)*kA + lambda*kB;
2480
2481   tabscale = table->scale;
2482   VFtab    = table->tab;
2483
2484   rt    = r*tabscale;
2485   n0    = rt;
2486   if (n0 >= table->n) {
2487     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",
2488               type,table_nr,r,n0,n0+1,table->n);
2489   }
2490   eps   = rt - n0;
2491   eps2  = eps*eps;
2492   nnn   = 4*n0;
2493   Yt    = VFtab[nnn];
2494   Ft    = VFtab[nnn+1];
2495   Geps  = VFtab[nnn+2]*eps;
2496   Heps2 = VFtab[nnn+3]*eps2;
2497   Fp    = Ft + Geps + Heps2;
2498   VV    = Yt + Fp*eps;
2499   FF    = Fp + Geps + 2.0*Heps2;
2500
2501   *F    = -k*FF*tabscale;
2502   *V    = k*VV;
2503   dvdl  = (kB - kA)*VV;
2504
2505   return dvdl;
2506
2507   /* That was 22 flops */
2508 }
2509
2510 real tab_bonds(int nbonds,
2511                const t_iatom forceatoms[],const t_iparams forceparams[],
2512                const rvec x[],rvec f[],rvec fshift[],
2513                const t_pbc *pbc,const t_graph *g,
2514                real lambda,real *dvdlambda,
2515                const t_mdatoms *md,t_fcdata *fcd,
2516                int *global_atom_index)
2517 {
2518   int  i,m,ki,ai,aj,type,table;
2519   real dr,dr2,fbond,vbond,fij,vtot;
2520   rvec dx;
2521   ivec dt;
2522
2523   vtot = 0.0;
2524   for(i=0; (i<nbonds); ) {
2525     type = forceatoms[i++];
2526     ai   = forceatoms[i++];
2527     aj   = forceatoms[i++];
2528
2529     ki   = pbc_rvec_sub(pbc,x[ai],x[aj],dx);    /*   3          */
2530     dr2  = iprod(dx,dx);                        /*   5          */
2531     dr   = dr2*gmx_invsqrt(dr2);                        /*  10          */
2532
2533     table = forceparams[type].tab.table;
2534
2535     *dvdlambda += bonded_tab("bond",table,
2536                              &fcd->bondtab[table],
2537                              forceparams[type].tab.kA,
2538                              forceparams[type].tab.kB,
2539                              dr,lambda,&vbond,&fbond);  /*  22 */
2540
2541     if (dr2 == 0.0)
2542       continue;
2543
2544
2545     vtot  += vbond;/* 1*/
2546     fbond *= gmx_invsqrt(dr2);                  /*   6          */
2547 #ifdef DEBUG
2548     if (debug)
2549       fprintf(debug,"TABBONDS: dr = %10g  vbond = %10g  fbond = %10g\n",
2550               dr,vbond,fbond);
2551 #endif
2552     if (g) {
2553       ivec_sub(SHIFT_IVEC(g,ai),SHIFT_IVEC(g,aj),dt);
2554       ki=IVEC2IS(dt);
2555     }
2556     for (m=0; (m<DIM); m++) {                   /*  15          */
2557       fij=fbond*dx[m];
2558       f[ai][m]+=fij;
2559       f[aj][m]-=fij;
2560       fshift[ki][m]+=fij;
2561       fshift[CENTRAL][m]-=fij;
2562     }
2563   }                                     /* 62 TOTAL     */
2564   return vtot;
2565 }
2566
2567 real tab_angles(int nbonds,
2568                 const t_iatom forceatoms[],const t_iparams forceparams[],
2569                 const rvec x[],rvec f[],rvec fshift[],
2570                 const t_pbc *pbc,const t_graph *g,
2571                 real lambda,real *dvdlambda,
2572                 const t_mdatoms *md,t_fcdata *fcd,
2573                 int *global_atom_index)
2574 {
2575   int  i,ai,aj,ak,t1,t2,type,table;
2576   rvec r_ij,r_kj;
2577   real cos_theta,cos_theta2,theta,dVdt,va,vtot;
2578   ivec jt,dt_ij,dt_kj;
2579
2580   vtot = 0.0;
2581   for(i=0; (i<nbonds); ) {
2582     type = forceatoms[i++];
2583     ai   = forceatoms[i++];
2584     aj   = forceatoms[i++];
2585     ak   = forceatoms[i++];
2586
2587     theta  = bond_angle(x[ai],x[aj],x[ak],pbc,
2588                         r_ij,r_kj,&cos_theta,&t1,&t2);  /*  41          */
2589
2590     table = forceparams[type].tab.table;
2591
2592     *dvdlambda += bonded_tab("angle",table,
2593                              &fcd->angletab[table],
2594                              forceparams[type].tab.kA,
2595                              forceparams[type].tab.kB,
2596                              theta,lambda,&va,&dVdt);  /*  22  */
2597     vtot += va;
2598
2599     cos_theta2 = sqr(cos_theta);                /*   1          */
2600     if (cos_theta2 < 1) {
2601       int  m;
2602       real snt,st,sth;
2603       real cik,cii,ckk;
2604       real nrkj2,nrij2;
2605       rvec f_i,f_j,f_k;
2606
2607       st  = dVdt*gmx_invsqrt(1 - cos_theta2);   /*  12          */
2608       sth = st*cos_theta;                       /*   1          */
2609 #ifdef DEBUG
2610       if (debug)
2611         fprintf(debug,"ANGLES: theta = %10g  vth = %10g  dV/dtheta = %10g\n",
2612                 theta*RAD2DEG,va,dVdt);
2613 #endif
2614       nrkj2=iprod(r_kj,r_kj);                   /*   5          */
2615       nrij2=iprod(r_ij,r_ij);
2616
2617       cik=st*gmx_invsqrt(nrkj2*nrij2);          /*  12          */
2618       cii=sth/nrij2;                            /*  10          */
2619       ckk=sth/nrkj2;                            /*  10          */
2620
2621       for (m=0; (m<DIM); m++) {                 /*  39          */
2622         f_i[m]=-(cik*r_kj[m]-cii*r_ij[m]);
2623         f_k[m]=-(cik*r_ij[m]-ckk*r_kj[m]);
2624         f_j[m]=-f_i[m]-f_k[m];
2625         f[ai][m]+=f_i[m];
2626         f[aj][m]+=f_j[m];
2627         f[ak][m]+=f_k[m];
2628       }
2629       if (g) {
2630         copy_ivec(SHIFT_IVEC(g,aj),jt);
2631
2632         ivec_sub(SHIFT_IVEC(g,ai),jt,dt_ij);
2633         ivec_sub(SHIFT_IVEC(g,ak),jt,dt_kj);
2634         t1=IVEC2IS(dt_ij);
2635         t2=IVEC2IS(dt_kj);
2636       }
2637       rvec_inc(fshift[t1],f_i);
2638       rvec_inc(fshift[CENTRAL],f_j);
2639       rvec_inc(fshift[t2],f_k);
2640     }                                           /* 169 TOTAL    */
2641   }
2642   return vtot;
2643 }
2644
2645 real tab_dihs(int nbonds,
2646               const t_iatom forceatoms[],const t_iparams forceparams[],
2647               const rvec x[],rvec f[],rvec fshift[],
2648               const t_pbc *pbc,const t_graph *g,
2649               real lambda,real *dvdlambda,
2650               const t_mdatoms *md,t_fcdata *fcd,
2651               int *global_atom_index)
2652 {
2653   int  i,type,ai,aj,ak,al,table;
2654   int  t1,t2,t3;
2655   rvec r_ij,r_kj,r_kl,m,n;
2656   real phi,sign,ddphi,vpd,vtot;
2657
2658   vtot = 0.0;
2659   for(i=0; (i<nbonds); ) {
2660     type = forceatoms[i++];
2661     ai   = forceatoms[i++];
2662     aj   = forceatoms[i++];
2663     ak   = forceatoms[i++];
2664     al   = forceatoms[i++];
2665
2666     phi=dih_angle(x[ai],x[aj],x[ak],x[al],pbc,r_ij,r_kj,r_kl,m,n,
2667                   &sign,&t1,&t2,&t3);                   /*  84  */
2668
2669     table = forceparams[type].tab.table;
2670
2671     /* Hopefully phi+M_PI never results in values < 0 */
2672     *dvdlambda += bonded_tab("dihedral",table,
2673                              &fcd->dihtab[table],
2674                              forceparams[type].tab.kA,
2675                              forceparams[type].tab.kB,
2676                              phi+M_PI,lambda,&vpd,&ddphi);
2677
2678     vtot += vpd;
2679     do_dih_fup(ai,aj,ak,al,-ddphi,r_ij,r_kj,r_kl,m,n,
2680                f,fshift,pbc,g,x,t1,t2,t3);                      /* 112  */
2681
2682 #ifdef DEBUG
2683     fprintf(debug,"pdih: (%d,%d,%d,%d) phi=%g\n",
2684             ai,aj,ak,al,phi);
2685 #endif
2686   } /* 227 TOTAL        */
2687
2688   return vtot;
2689 }
2690
2691 void calc_bonds(FILE *fplog,const gmx_multisim_t *ms,
2692                 const t_idef *idef,
2693                 rvec x[],history_t *hist,
2694                 rvec f[],t_forcerec *fr,
2695                 const t_pbc *pbc,const t_graph *g,
2696                 gmx_enerdata_t *enerd,t_nrnb *nrnb,
2697                 real lambda,
2698                 const t_mdatoms *md,
2699                 t_fcdata *fcd,int *global_atom_index,
2700                 t_atomtypes *atype, gmx_genborn_t *born,
2701                 gmx_bool bPrintSepPot,gmx_large_int_t step)
2702 {
2703   int    ftype,nbonds,ind,nat1;
2704   real   *epot,v,dvdl;
2705   const  t_pbc *pbc_null;
2706   char   buf[22];
2707
2708   if (fr->bMolPBC)
2709     pbc_null = pbc;
2710   else
2711     pbc_null = NULL;
2712
2713   if (bPrintSepPot)
2714     fprintf(fplog,"Step %s: bonded V and dVdl for this node\n",
2715             gmx_step_str(step,buf));
2716
2717 #ifdef DEBUG
2718   if (g && debug)
2719     p_graph(debug,"Bondage is fun",g);
2720 #endif
2721
2722   epot = enerd->term;
2723
2724   /* Do pre force calculation stuff which might require communication */
2725   if (idef->il[F_ORIRES].nr) {
2726     epot[F_ORIRESDEV] = calc_orires_dev(ms,idef->il[F_ORIRES].nr,
2727                                         idef->il[F_ORIRES].iatoms,
2728                                         idef->iparams,md,(const rvec*)x,
2729                                         pbc_null,fcd,hist);
2730   }
2731   if (idef->il[F_DISRES].nr) {
2732     calc_disres_R_6(ms,idef->il[F_DISRES].nr,
2733                     idef->il[F_DISRES].iatoms,
2734                     idef->iparams,(const rvec*)x,pbc_null,
2735                     fcd,hist);
2736   }
2737
2738   /* Loop over all bonded force types to calculate the bonded forces */
2739   for(ftype=0; (ftype<F_NRE); ftype++) {
2740           if(ftype<F_GB12 || ftype>F_GB14) {
2741     if ((interaction_function[ftype].flags & IF_BOND) &&
2742         !(ftype == F_CONNBONDS || ftype == F_POSRES)) {
2743       nbonds=idef->il[ftype].nr;
2744       if (nbonds > 0) {
2745         ind  = interaction_function[ftype].nrnb_ind;
2746         nat1 = interaction_function[ftype].nratoms + 1;
2747         dvdl = 0;
2748         if (ftype < F_LJ14 || ftype > F_LJC_PAIRS_NB) {
2749                 if(ftype==F_CMAP)
2750                 {
2751                         v = cmap_dihs(nbonds,idef->il[ftype].iatoms,
2752                                                   idef->iparams,&idef->cmap_grid,
2753                                                   (const rvec*)x,f,fr->fshift,
2754                                                   pbc_null,g,lambda,&dvdl,md,fcd,
2755                                                   global_atom_index);
2756                 }
2757                 else
2758                 {
2759                         v =
2760             interaction_function[ftype].ifunc(nbonds,idef->il[ftype].iatoms,
2761                                               idef->iparams,
2762                                               (const rvec*)x,f,fr->fshift,
2763                                               pbc_null,g,lambda,&dvdl,md,fcd,
2764                                               global_atom_index);
2765                 }
2766
2767           if (bPrintSepPot) {
2768             fprintf(fplog,"  %-23s #%4d  V %12.5e  dVdl %12.5e\n",
2769                     interaction_function[ftype].longname,nbonds/nat1,v,dvdl);
2770           }
2771         } else {
2772           v = do_listed_vdw_q(ftype,nbonds,idef->il[ftype].iatoms,
2773                               idef->iparams,
2774                               (const rvec*)x,f,fr->fshift,
2775                               pbc_null,g,
2776                               lambda,&dvdl,
2777                               md,fr,&enerd->grpp,global_atom_index);
2778           if (bPrintSepPot) {
2779             fprintf(fplog,"  %-5s + %-15s #%4d                  dVdl %12.5e\n",
2780                     interaction_function[ftype].longname,
2781                     interaction_function[F_COUL14].longname,nbonds/nat1,dvdl);
2782           }
2783         }
2784         if (ind != -1)
2785           inc_nrnb(nrnb,ind,nbonds/nat1);
2786         epot[ftype]        += v;
2787         enerd->dvdl_nonlin += dvdl;
2788       }
2789     }
2790   }
2791   }
2792   /* Copy the sum of violations for the distance restraints from fcd */
2793   if (fcd)
2794     epot[F_DISRESVIOL] = fcd->disres.sumviol;
2795 }
2796
2797 void calc_bonds_lambda(FILE *fplog,
2798                        const t_idef *idef,
2799                        rvec x[],
2800                        t_forcerec *fr,
2801                        const t_pbc *pbc,const t_graph *g,
2802                        gmx_enerdata_t *enerd,t_nrnb *nrnb,
2803                        real lambda,
2804                        const t_mdatoms *md,
2805                        t_fcdata *fcd,int *global_atom_index)
2806 {
2807     int    ftype,nbonds_np,nbonds,ind, nat1;
2808   real   *epot,v,dvdl;
2809   rvec   *f,*fshift_orig;
2810   const  t_pbc *pbc_null;
2811   t_iatom *iatom_fe;
2812
2813   if (fr->bMolPBC)
2814     pbc_null = pbc;
2815   else
2816     pbc_null = NULL;
2817
2818   epot = enerd->term;
2819
2820   snew(f,fr->natoms_force);
2821   /* We want to preserve the fshift array in forcerec */
2822   fshift_orig = fr->fshift;
2823   snew(fr->fshift,SHIFTS);
2824
2825   /* Loop over all bonded force types to calculate the bonded forces */
2826   for(ftype=0; (ftype<F_NRE); ftype++) {
2827       if(ftype<F_GB12 || ftype>F_GB14) {
2828
2829           if ((interaction_function[ftype].flags & IF_BOND) &&
2830               !(ftype == F_CONNBONDS || ftype == F_POSRES))
2831           {
2832               nbonds_np = idef->il[ftype].nr_nonperturbed;
2833               nbonds    = idef->il[ftype].nr - nbonds_np;
2834               nat1 = interaction_function[ftype].nratoms + 1;
2835               if (nbonds > 0) {
2836                   ind  = interaction_function[ftype].nrnb_ind;
2837                   iatom_fe = idef->il[ftype].iatoms + nbonds_np;
2838                   dvdl = 0;
2839                   if (ftype < F_LJ14 || ftype > F_LJC_PAIRS_NB) {
2840                       v =
2841                           interaction_function[ftype].ifunc(nbonds,iatom_fe,
2842                                                             idef->iparams,
2843                                                             (const rvec*)x,f,fr->fshift,
2844                                                             pbc_null,g,lambda,&dvdl,md,fcd,
2845                                                             global_atom_index);
2846                   } else {
2847                       v = do_listed_vdw_q(ftype,nbonds,iatom_fe,
2848                                           idef->iparams,
2849                                           (const rvec*)x,f,fr->fshift,
2850                                           pbc_null,g,
2851                                           lambda,&dvdl,
2852                                           md,fr,&enerd->grpp,global_atom_index);
2853                   }
2854                   if (ind != -1)
2855                       inc_nrnb(nrnb,ind,nbonds/nat1);
2856                   epot[ftype] += v;
2857               }
2858           }
2859       }
2860   }
2861
2862   sfree(fr->fshift);
2863   fr->fshift = fshift_orig;
2864   sfree(f);
2865 }