src/gromacs/gmxlib/bondfree.c

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