src/gmxlib/bondfree.c

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