src/gromacs/gmxlib/bondfree.c

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