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