src/gromacs/mdlib/genborn.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2008, The GROMACS development team.
   6  * Copyright (c) 2013,2014,2015, by the GROMACS development team, led by
   7  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   8  * and including many others, as listed in the AUTHORS file in the
   9  * top-level source directory and at http://www.gromacs.org.
  10  *
  11  * GROMACS is free software; you can redistribute it and/or
  12  * modify it under the terms of the GNU Lesser General Public License
  13  * as published by the Free Software Foundation; either version 2.1
  14  * of the License, or (at your option) any later version.
  15  *
  16  * GROMACS is distributed in the hope that it will be useful,
  17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19  * Lesser General Public License for more details.
  20  *
  21  * You should have received a copy of the GNU Lesser General Public
  22  * License along with GROMACS; if not, see
  23  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  24  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  25  *
  26  * If you want to redistribute modifications to GROMACS, please
  27  * consider that scientific software is very special. Version
  28  * control is crucial - bugs must be traceable. We will be happy to
  29  * consider code for inclusion in the official distribution, but
  30  * derived work must not be called official GROMACS. Details are found
  31  * in the README & COPYING files - if they are missing, get the
  32  * official version at http://www.gromacs.org.
  33  *
  34  * To help us fund GROMACS development, we humbly ask that you cite
  35  * the research papers on the package. Check out http://www.gromacs.org.
  36  */
  37
  38 #include "gmxpre.h"
  39
  40 #include "gromacs/legacyheaders/genborn.h"
  41
  42 #include <math.h>
  43 #include <string.h>
  44
  45 #include "gromacs/domdec/domdec.h"
  46 #include "gromacs/fileio/pdbio.h"
  47 #include "gromacs/legacyheaders/names.h"
  48 #include "gromacs/legacyheaders/network.h"
  49 #include "gromacs/legacyheaders/nrnb.h"
  50 #include "gromacs/legacyheaders/typedefs.h"
  51 #include "gromacs/legacyheaders/types/commrec.h"
  52 #include "gromacs/math/units.h"
  53 #include "gromacs/math/vec.h"
  54 #include "gromacs/mdlib/genborn_allvsall.h"
  55 #include "gromacs/pbcutil/ishift.h"
  56 #include "gromacs/pbcutil/mshift.h"
  57 #include "gromacs/pbcutil/pbc.h"
  58 #include "gromacs/topology/mtop_util.h"
  59 #include "gromacs/utility/fatalerror.h"
  60 #include "gromacs/utility/gmxmpi.h"
  61 #include "gromacs/utility/smalloc.h"
  62
  63
  64 typedef struct {
  65     int  shift;
  66     int  naj;
  67     int *aj;
  68     int  aj_nalloc;
  69 } gbtmpnbl_t;
  70
  71 typedef struct gbtmpnbls {
  72     int         nlist;
  73     gbtmpnbl_t *list;
  74     int         list_nalloc;
  75 } t_gbtmpnbls;
  76
  77 /* This function is exactly the same as the one in listed-forces/bonded.cpp. The reason
  78  * it is copied here is that the bonded gb-interactions are evaluated
  79  * not in calc_bonds, but rather in calc_gb_forces
  80  */
  81 static int pbc_rvec_sub(const t_pbc *pbc, const rvec xi, const rvec xj, rvec dx)
  82 {
  83     if (pbc)
  84     {
  85         return pbc_dx_aiuc(pbc, xi, xj, dx);
  86     }
  87     else
  88     {
  89         rvec_sub(xi, xj, dx);
  90         return CENTRAL;
  91     }
  92 }
  93
  94 static int init_gb_nblist(int natoms, t_nblist *nl)
  95 {
  96     nl->maxnri      = natoms*4;
  97     nl->maxnrj      = 0;
  98     nl->nri         = 0;
  99     nl->nrj         = 0;
 100     nl->iinr        = NULL;
 101     nl->gid         = NULL;
 102     nl->shift       = NULL;
 103     nl->jindex      = NULL;
 104     nl->jjnr        = NULL;
 105     /*nl->nltype      = nltype;*/
 106
 107     srenew(nl->iinr,   nl->maxnri);
 108     srenew(nl->gid,    nl->maxnri);
 109     srenew(nl->shift,  nl->maxnri);
 110     srenew(nl->jindex, nl->maxnri+1);
 111
 112     nl->jindex[0] = 0;
 113
 114     return 0;
 115 }
 116
 117
 118 static int init_gb_still(const t_atomtypes *atype, t_idef *idef, t_atoms *atoms,
 119                          gmx_genborn_t *born, int natoms)
 120 {
 121
 122     int   i, j, i1, i2, k, m, nbond, nang, ia, ib, ic, id, nb, idx, idx2, at;
 123     int   iam, ibm;
 124     int   at0, at1;
 125     real  length, angle;
 126     real  r, ri, rj, ri2, ri3, rj2, r2, r3, r4, rk, ratio, term, h, doffset;
 127     real  p1, p2, p3, factor, cosine, rab, rbc;
 128
 129     real *vsol;
 130     real *gp;
 131
 132     snew(vsol, natoms);
 133     snew(gp, natoms);
 134     snew(born->gpol_still_work, natoms+3);
 135
 136     at0 = 0;
 137     at1 = natoms;
 138
 139     doffset = born->gb_doffset;
 140
 141     for (i = 0; i < natoms; i++)
 142     {
 143         born->gpol_globalindex[i]              = born->vsolv_globalindex[i] =
 144                 born->gb_radius_globalindex[i] = 0;
 145     }
 146
 147     /* Compute atomic solvation volumes for Still method */
 148     for (i = 0; i < natoms; i++)
 149     {
 150         ri = atype->gb_radius[atoms->atom[i].type];
 151         born->gb_radius_globalindex[i] = ri;
 152         r3 = ri*ri*ri;
 153         born->vsolv_globalindex[i] = (4*M_PI/3)*r3;
 154     }
 155
 156     for (j = 0; j < idef->il[F_GB12].nr; j += 3)
 157     {
 158         m  = idef->il[F_GB12].iatoms[j];
 159         ia = idef->il[F_GB12].iatoms[j+1];
 160         ib = idef->il[F_GB12].iatoms[j+2];
 161
 162         r = 1.01*idef->iparams[m].gb.st;
 163
 164         ri   = atype->gb_radius[atoms->atom[ia].type];
 165         rj   = atype->gb_radius[atoms->atom[ib].type];
 166
 167         ri2  = ri*ri;
 168         ri3  = ri2*ri;
 169         rj2  = rj*rj;
 170
 171         ratio  = (rj2-ri2-r*r)/(2*ri*r);
 172         h      = ri*(1+ratio);
 173         term   = (M_PI/3.0)*h*h*(3.0*ri-h);
 174
 175         born->vsolv_globalindex[ia] -= term;
 176
 177         ratio  = (ri2-rj2-r*r)/(2*rj*r);
 178         h      = rj*(1+ratio);
 179         term   = (M_PI/3.0)*h*h*(3.0*rj-h);
 180
 181         born->vsolv_globalindex[ib] -= term;
 182     }
 183
 184     /* Get the self-, 1-2 and 1-3 polarization energies for analytical Still
 185        method */
 186     /* Self */
 187     for (j = 0; j < natoms; j++)
 188     {
 189         if (born->use_globalindex[j] == 1)
 190         {
 191             born->gpol_globalindex[j] = -0.5*ONE_4PI_EPS0/
 192                 (atype->gb_radius[atoms->atom[j].type]-doffset+STILL_P1);
 193         }
 194     }
 195
 196     /* 1-2 */
 197     for (j = 0; j < idef->il[F_GB12].nr; j += 3)
 198     {
 199         m  = idef->il[F_GB12].iatoms[j];
 200         ia = idef->il[F_GB12].iatoms[j+1];
 201         ib = idef->il[F_GB12].iatoms[j+2];
 202
 203         r = idef->iparams[m].gb.st;
 204
 205         r4 = r*r*r*r;
 206
 207         born->gpol_globalindex[ia] = born->gpol_globalindex[ia]+
 208             STILL_P2*born->vsolv_globalindex[ib]/r4;
 209         born->gpol_globalindex[ib] = born->gpol_globalindex[ib]+
 210             STILL_P2*born->vsolv_globalindex[ia]/r4;
 211     }
 212
 213     /* 1-3 */
 214     for (j = 0; j < idef->il[F_GB13].nr; j += 3)
 215     {
 216         m  = idef->il[F_GB13].iatoms[j];
 217         ia = idef->il[F_GB13].iatoms[j+1];
 218         ib = idef->il[F_GB13].iatoms[j+2];
 219
 220         r  = idef->iparams[m].gb.st;
 221         r4 = r*r*r*r;
 222
 223         born->gpol_globalindex[ia] = born->gpol_globalindex[ia]+
 224             STILL_P3*born->vsolv_globalindex[ib]/r4;
 225         born->gpol_globalindex[ib] = born->gpol_globalindex[ib]+
 226             STILL_P3*born->vsolv_globalindex[ia]/r4;
 227     }
 228
 229     sfree(vsol);
 230     sfree(gp);
 231
 232     return 0;
 233 }
 234
 235 /* Initialize all GB datastructs and compute polarization energies */
 236 int init_gb(gmx_genborn_t **p_born,
 237             t_forcerec *fr, const t_inputrec *ir,
 238             const gmx_mtop_t *mtop, int gb_algorithm)
 239 {
 240     int             i, j, m, ai, aj, jj, natoms, nalloc;
 241     real            rai, sk, p, doffset;
 242
 243     t_atoms         atoms;
 244     gmx_genborn_t  *born;
 245     gmx_localtop_t *localtop;
 246
 247     natoms   = mtop->natoms;
 248
 249     atoms    = gmx_mtop_global_atoms(mtop);
 250     localtop = gmx_mtop_generate_local_top(mtop, ir);
 251
 252     snew(born, 1);
 253     *p_born = born;
 254
 255     born->nr  = natoms;
 256
 257     snew(born->drobc, natoms);
 258     snew(born->bRad,  natoms);
 259
 260     /* Allocate memory for the global data arrays */
 261     snew(born->param_globalindex, natoms+3);
 262     snew(born->gpol_globalindex,  natoms+3);
 263     snew(born->vsolv_globalindex, natoms+3);
 264     snew(born->gb_radius_globalindex, natoms+3);
 265     snew(born->use_globalindex,    natoms+3);
 266
 267     snew(fr->invsqrta, natoms);
 268     snew(fr->dvda,     natoms);
 269
 270     fr->dadx              = NULL;
 271     fr->dadx_rawptr       = NULL;
 272     fr->nalloc_dadx       = 0;
 273     born->gpol_still_work = NULL;
 274     born->gpol_hct_work   = NULL;
 275
 276     /* snew(born->asurf,natoms); */
 277     /* snew(born->dasurf,natoms); */
 278
 279     /* Initialize the gb neighbourlist */
 280     init_gb_nblist(natoms, &(fr->gblist));
 281
 282     /* Do the Vsites exclusions (if any) */
 283     for (i = 0; i < natoms; i++)
 284     {
 285         jj = atoms.atom[i].type;
 286         if (mtop->atomtypes.gb_radius[atoms.atom[i].type] > 0)
 287         {
 288             born->use_globalindex[i] = 1;
 289         }
 290         else
 291         {
 292             born->use_globalindex[i] = 0;
 293         }
 294
 295         /* If we have a Vsite, put vs_globalindex[i]=0 */
 296         if (C6 (fr->nbfp, fr->ntype, jj, jj) == 0 &&
 297             C12(fr->nbfp, fr->ntype, jj, jj) == 0 &&
 298             atoms.atom[i].q == 0)
 299         {
 300             born->use_globalindex[i] = 0;
 301         }
 302     }
 303
 304     /* Copy algorithm parameters from inputrecord to local structure */
 305     born->obc_alpha          = ir->gb_obc_alpha;
 306     born->obc_beta           = ir->gb_obc_beta;
 307     born->obc_gamma          = ir->gb_obc_gamma;
 308     born->gb_doffset         = ir->gb_dielectric_offset;
 309     born->gb_epsilon_solvent = ir->gb_epsilon_solvent;
 310     born->epsilon_r          = ir->epsilon_r;
 311
 312     doffset = born->gb_doffset;
 313
 314     /* Set the surface tension */
 315     born->sa_surface_tension = ir->sa_surface_tension;
 316
 317     /* If Still model, initialise the polarisation energies */
 318     if (gb_algorithm == egbSTILL)
 319     {
 320         init_gb_still(&(mtop->atomtypes), &(localtop->idef), &atoms,
 321                       born, natoms);
 322     }
 323
 324
 325     /* If HCT/OBC,  precalculate the sk*atype->S_hct factors */
 326     else if (gb_algorithm == egbHCT || gb_algorithm == egbOBC)
 327     {
 328
 329         snew(born->gpol_hct_work, natoms+3);
 330
 331         for (i = 0; i < natoms; i++)
 332         {
 333             if (born->use_globalindex[i] == 1)
 334             {
 335                 rai = mtop->atomtypes.gb_radius[atoms.atom[i].type]-doffset;
 336                 sk  = rai * mtop->atomtypes.S_hct[atoms.atom[i].type];
 337                 born->param_globalindex[i]     = sk;
 338                 born->gb_radius_globalindex[i] = rai;
 339             }
 340             else
 341             {
 342                 born->param_globalindex[i]     = 0;
 343                 born->gb_radius_globalindex[i] = 0;
 344             }
 345         }
 346     }
 347
 348     /* Allocate memory for work arrays for temporary use */
 349     snew(born->work, natoms+4);
 350     snew(born->count, natoms);
 351     snew(born->nblist_work, natoms);
 352
 353     /* Domain decomposition specific stuff */
 354     born->nalloc = 0;
 355
 356     return 0;
 357 }
 358
 359
 360
 361 static int
 362 calc_gb_rad_still(t_commrec *cr, t_forcerec *fr, gmx_localtop_t *top,
 363                   rvec x[], t_nblist *nl,
 364                   gmx_genborn_t *born, t_mdatoms *md)
 365 {
 366     int  i, k, n, nj0, nj1, ai, aj, type;
 367     int  shift;
 368     real shX, shY, shZ;
 369     real gpi, dr, dr2, dr4, idr4, rvdw, ratio, ccf, theta, term, rai, raj;
 370     real ix1, iy1, iz1, jx1, jy1, jz1, dx11, dy11, dz11;
 371     real rinv, idr2, idr6, vaj, dccf, cosq, sinq, prod, gpi2;
 372     real factor;
 373     real vai, prod_ai, icf4, icf6;
 374
 375     factor  = 0.5*ONE_4PI_EPS0;
 376     n       = 0;
 377
 378     for (i = 0; i < born->nr; i++)
 379     {
 380         born->gpol_still_work[i] = 0;
 381     }
 382
 383     for (i = 0; i < nl->nri; i++)
 384     {
 385         ai      = nl->iinr[i];
 386
 387         nj0     = nl->jindex[i];
 388         nj1     = nl->jindex[i+1];
 389
 390         /* Load shifts for this list */
 391         shift   = nl->shift[i];
 392         shX     = fr->shift_vec[shift][0];
 393         shY     = fr->shift_vec[shift][1];
 394         shZ     = fr->shift_vec[shift][2];
 395
 396         gpi     = 0;
 397
 398         rai     = top->atomtypes.gb_radius[md->typeA[ai]];
 399         vai     = born->vsolv[ai];
 400         prod_ai = STILL_P4*vai;
 401
 402         /* Load atom i coordinates, add shift vectors */
 403         ix1     = shX + x[ai][0];
 404         iy1     = shY + x[ai][1];
 405         iz1     = shZ + x[ai][2];
 406
 407         for (k = nj0; k < nj1 && nl->jjnr[k] >= 0; k++)
 408         {
 409             aj    = nl->jjnr[k];
 410             jx1   = x[aj][0];
 411             jy1   = x[aj][1];
 412             jz1   = x[aj][2];
 413
 414             dx11  = ix1-jx1;
 415             dy11  = iy1-jy1;
 416             dz11  = iz1-jz1;
 417
 418             dr2   = dx11*dx11+dy11*dy11+dz11*dz11;
 419             rinv  = gmx_invsqrt(dr2);
 420             idr2  = rinv*rinv;
 421             idr4  = idr2*idr2;
 422             idr6  = idr4*idr2;
 423
 424             raj = top->atomtypes.gb_radius[md->typeA[aj]];
 425
 426             rvdw  = rai + raj;
 427
 428             ratio = dr2 / (rvdw * rvdw);
 429             vaj   = born->vsolv[aj];
 430
 431             if (ratio > STILL_P5INV)
 432             {
 433                 ccf  = 1.0;
 434                 dccf = 0.0;
 435             }
 436             else
 437             {
 438                 theta = ratio*STILL_PIP5;
 439                 cosq  = cos(theta);
 440                 term  = 0.5*(1.0-cosq);
 441                 ccf   = term*term;
 442                 sinq  = 1.0 - cosq*cosq;
 443                 dccf  = 2.0*term*sinq*gmx_invsqrt(sinq)*theta;
 444             }
 445
 446             prod                       = STILL_P4*vaj;
 447             icf4                       = ccf*idr4;
 448             icf6                       = (4*ccf-dccf)*idr6;
 449             born->gpol_still_work[aj] += prod_ai*icf4;
 450             gpi                        = gpi+prod*icf4;
 451
 452             /* Save ai->aj and aj->ai chain rule terms */
 453             fr->dadx[n++]   = prod*icf6;
 454             fr->dadx[n++]   = prod_ai*icf6;
 455         }
 456         born->gpol_still_work[ai] += gpi;
 457     }
 458
 459     /* Parallel summations */
 460     if (DOMAINDECOMP(cr))
 461     {
 462         dd_atom_sum_real(cr->dd, born->gpol_still_work);
 463     }
 464
 465     /* Calculate the radii */
 466     for (i = 0; i < fr->natoms_force; i++) /* PELA born->nr */
 467     {
 468         if (born->use[i] != 0)
 469         {
 470             gpi             = born->gpol[i]+born->gpol_still_work[i];
 471             gpi2            = gpi * gpi;
 472             born->bRad[i]   = factor*gmx_invsqrt(gpi2);
 473             fr->invsqrta[i] = gmx_invsqrt(born->bRad[i]);
 474         }
 475     }
 476
 477     /* Extra communication required for DD */
 478     if (DOMAINDECOMP(cr))
 479     {
 480         dd_atom_spread_real(cr->dd, born->bRad);
 481         dd_atom_spread_real(cr->dd, fr->invsqrta);
 482     }
 483
 484     return 0;
 485
 486 }
 487
 488
 489 static int
 490 calc_gb_rad_hct(t_commrec *cr, t_forcerec *fr, gmx_localtop_t *top,
 491                 rvec x[], t_nblist *nl,
 492                 gmx_genborn_t *born, t_mdatoms *md)
 493 {
 494     int   i, k, n, ai, aj, nj0, nj1, at0, at1;
 495     int   shift;
 496     real  shX, shY, shZ;
 497     real  rai, raj, gpi, dr2, dr, sk, sk_ai, sk2, sk2_ai, lij, uij, diff2, tmp, sum_ai;
 498     real  rad, min_rad, rinv, rai_inv;
 499     real  ix1, iy1, iz1, jx1, jy1, jz1, dx11, dy11, dz11;
 500     real  lij2, uij2, lij3, uij3, t1, t2, t3;
 501     real  lij_inv, dlij, duij, sk2_rinv, prod, log_term;
 502     real  doffset, raj_inv, dadx_val;
 503     real *gb_radius;
 504
 505     doffset   = born->gb_doffset;
 506     gb_radius = born->gb_radius;
 507
 508     for (i = 0; i < born->nr; i++)
 509     {
 510         born->gpol_hct_work[i] = 0;
 511     }
 512
 513     /* Keep the compiler happy */
 514     n    = 0;
 515     prod = 0;
 516
 517     for (i = 0; i < nl->nri; i++)
 518     {
 519         ai     = nl->iinr[i];
 520
 521         nj0    = nl->jindex[i];
 522         nj1    = nl->jindex[i+1];
 523
 524         /* Load shifts for this list */
 525         shift   = nl->shift[i];
 526         shX     = fr->shift_vec[shift][0];
 527         shY     = fr->shift_vec[shift][1];
 528         shZ     = fr->shift_vec[shift][2];
 529
 530         rai     = gb_radius[ai];
 531         rai_inv = 1.0/rai;
 532
 533         sk_ai   = born->param[ai];
 534         sk2_ai  = sk_ai*sk_ai;
 535
 536         /* Load atom i coordinates, add shift vectors */
 537         ix1     = shX + x[ai][0];
 538         iy1     = shY + x[ai][1];
 539         iz1     = shZ + x[ai][2];
 540
 541         sum_ai  = 0;
 542
 543         for (k = nj0; k < nj1 && nl->jjnr[k] >= 0; k++)
 544         {
 545             aj    = nl->jjnr[k];
 546
 547             jx1   = x[aj][0];
 548             jy1   = x[aj][1];
 549             jz1   = x[aj][2];
 550
 551             dx11  = ix1 - jx1;
 552             dy11  = iy1 - jy1;
 553             dz11  = iz1 - jz1;
 554
 555             dr2   = dx11*dx11+dy11*dy11+dz11*dz11;
 556             rinv  = gmx_invsqrt(dr2);
 557             dr    = rinv*dr2;
 558
 559             sk    = born->param[aj];
 560             raj   = gb_radius[aj];
 561
 562             /* aj -> ai interaction */
 563             if (rai < dr+sk)
 564             {
 565                 lij     = 1.0/(dr-sk);
 566                 dlij    = 1.0;
 567
 568                 if (rai > dr-sk)
 569                 {
 570                     lij  = rai_inv;
 571                     dlij = 0.0;
 572                 }
 573
 574                 lij2     = lij*lij;
 575                 lij3     = lij2*lij;
 576
 577                 uij      = 1.0/(dr+sk);
 578                 uij2     = uij*uij;
 579                 uij3     = uij2*uij;
 580
 581                 diff2    = uij2-lij2;
 582
 583                 lij_inv  = gmx_invsqrt(lij2);
 584                 sk2      = sk*sk;
 585                 sk2_rinv = sk2*rinv;
 586                 prod     = 0.25*sk2_rinv;
 587
 588                 log_term = log(uij*lij_inv);
 589
 590                 tmp      = lij-uij + 0.25*dr*diff2 + (0.5*rinv)*log_term +
 591                     prod*(-diff2);
 592
 593                 if (rai < sk-dr)
 594                 {
 595                     tmp = tmp + 2.0 * (rai_inv-lij);
 596                 }
 597
 598                 t1 = 0.5*lij2 + prod*lij3 - 0.25*(lij*rinv+lij3*dr);
 599                 t2 = -0.5*uij2 - 0.25*sk2_rinv*uij3 + 0.25*(uij*rinv+uij3*dr);
 600                 t3 = 0.125*(1.0+sk2_rinv*rinv)*(-diff2)+0.25*log_term*rinv*rinv;
 601
 602                 dadx_val = (dlij*t1+t2+t3)*rinv; /* rb2 is moved to chainrule */
 603                 /* fr->dadx[n++] = (dlij*t1+duij*t2+t3)*rinv; */
 604                 /* rb2 is moved to chainrule    */
 605
 606                 sum_ai += 0.5*tmp;
 607             }
 608             else
 609             {
 610                 dadx_val = 0.0;
 611             }
 612             fr->dadx[n++] = dadx_val;
 613
 614
 615             /* ai -> aj interaction */
 616             if (raj < dr + sk_ai)
 617             {
 618                 lij     = 1.0/(dr-sk_ai);
 619                 dlij    = 1.0;
 620                 raj_inv = 1.0/raj;
 621
 622                 if (raj > dr-sk_ai)
 623                 {
 624                     lij  = raj_inv;
 625                     dlij = 0.0;
 626                 }
 627
 628                 lij2     = lij  * lij;
 629                 lij3     = lij2 * lij;
 630
 631                 uij      = 1.0/(dr+sk_ai);
 632                 uij2     = uij  * uij;
 633                 uij3     = uij2 * uij;
 634
 635                 diff2    = uij2-lij2;
 636
 637                 lij_inv  = gmx_invsqrt(lij2);
 638                 sk2      =  sk2_ai; /* sk2_ai = sk_ai * sk_ai in i loop above */
 639                 sk2_rinv = sk2*rinv;
 640                 prod     = 0.25 * sk2_rinv;
 641
 642                 /* log_term = table_log(uij*lij_inv,born->log_table,
 643                    LOG_TABLE_ACCURACY); */
 644                 log_term = log(uij*lij_inv);
 645
 646                 tmp      = lij-uij + 0.25*dr*diff2 + (0.5*rinv)*log_term +
 647                     prod*(-diff2);
 648
 649                 if (raj < sk_ai-dr)
 650                 {
 651                     tmp     = tmp + 2.0 * (raj_inv-lij);
 652                 }
 653
 654                 /* duij = 1.0 */
 655                 t1      = 0.5*lij2 + prod*lij3 - 0.25*(lij*rinv+lij3*dr);
 656                 t2      = -0.5*uij2 - 0.25*sk2_rinv*uij3 + 0.25*(uij*rinv+uij3*dr);
 657                 t3      = 0.125*(1.0+sk2_rinv*rinv)*(-diff2)+0.25*log_term*rinv*rinv;
 658
 659                 dadx_val = (dlij*t1+t2+t3)*rinv; /* rb2 is moved to chainrule    */
 660                 /* fr->dadx[n++] = (dlij*t1+duij*t2+t3)*rinv; */ /* rb2 is moved to chainrule    */
 661
 662                 born->gpol_hct_work[aj] += 0.5*tmp;
 663             }
 664             else
 665             {
 666                 dadx_val = 0.0;
 667             }
 668             fr->dadx[n++] = dadx_val;
 669         }
 670
 671         born->gpol_hct_work[ai] += sum_ai;
 672     }
 673
 674     /* Parallel summations */
 675     if (DOMAINDECOMP(cr))
 676     {
 677         dd_atom_sum_real(cr->dd, born->gpol_hct_work);
 678     }
 679
 680     for (i = 0; i < fr->natoms_force; i++) /* PELA born->nr */
 681     {
 682         if (born->use[i] != 0)
 683         {
 684             rai     = top->atomtypes.gb_radius[md->typeA[i]]-doffset;
 685             sum_ai  = 1.0/rai - born->gpol_hct_work[i];
 686             min_rad = rai + doffset;
 687             rad     = 1.0/sum_ai;
 688
 689             born->bRad[i]   = rad > min_rad ? rad : min_rad;
 690             fr->invsqrta[i] = gmx_invsqrt(born->bRad[i]);
 691         }
 692     }
 693
 694     /* Extra communication required for DD */
 695     if (DOMAINDECOMP(cr))
 696     {
 697         dd_atom_spread_real(cr->dd, born->bRad);
 698         dd_atom_spread_real(cr->dd, fr->invsqrta);
 699     }
 700
 701
 702     return 0;
 703 }
 704
 705 static int
 706 calc_gb_rad_obc(t_commrec *cr, t_forcerec *fr, gmx_localtop_t *top,
 707                 rvec x[], t_nblist *nl, gmx_genborn_t *born, t_mdatoms *md)
 708 {
 709     int   i, k, ai, aj, nj0, nj1, n, at0, at1;
 710     int   shift;
 711     real  shX, shY, shZ;
 712     real  rai, raj, gpi, dr2, dr, sk, sk2, lij, uij, diff2, tmp, sum_ai;
 713     real  rad, min_rad, sum_ai2, sum_ai3, tsum, tchain, rinv, rai_inv, lij_inv, rai_inv2;
 714     real  log_term, prod, sk2_rinv, sk_ai, sk2_ai;
 715     real  ix1, iy1, iz1, jx1, jy1, jz1, dx11, dy11, dz11;
 716     real  lij2, uij2, lij3, uij3, dlij, duij, t1, t2, t3;
 717     real  doffset, raj_inv, dadx_val;
 718     real *gb_radius;
 719
 720     /* Keep the compiler happy */
 721     n    = 0;
 722     prod = 0;
 723     raj  = 0;
 724
 725     doffset   = born->gb_doffset;
 726     gb_radius = born->gb_radius;
 727
 728     for (i = 0; i < born->nr; i++)
 729     {
 730         born->gpol_hct_work[i] = 0;
 731     }
 732
 733     for (i = 0; i < nl->nri; i++)
 734     {
 735         ai      = nl->iinr[i];
 736
 737         nj0     = nl->jindex[i];
 738         nj1     = nl->jindex[i+1];
 739
 740         /* Load shifts for this list */
 741         shift   = nl->shift[i];
 742         shX     = fr->shift_vec[shift][0];
 743         shY     = fr->shift_vec[shift][1];
 744         shZ     = fr->shift_vec[shift][2];
 745
 746         rai      = gb_radius[ai];
 747         rai_inv  = 1.0/rai;
 748
 749         sk_ai    = born->param[ai];
 750         sk2_ai   = sk_ai*sk_ai;
 751
 752         /* Load atom i coordinates, add shift vectors */
 753         ix1      = shX + x[ai][0];
 754         iy1      = shY + x[ai][1];
 755         iz1      = shZ + x[ai][2];
 756
 757         sum_ai   = 0;
 758
 759         for (k = nj0; k < nj1 && nl->jjnr[k] >= 0; k++)
 760         {
 761             aj    = nl->jjnr[k];
 762
 763             jx1   = x[aj][0];
 764             jy1   = x[aj][1];
 765             jz1   = x[aj][2];
 766
 767             dx11  = ix1 - jx1;
 768             dy11  = iy1 - jy1;
 769             dz11  = iz1 - jz1;
 770
 771             dr2   = dx11*dx11+dy11*dy11+dz11*dz11;
 772             rinv  = gmx_invsqrt(dr2);
 773             dr    = dr2*rinv;
 774
 775             /* sk is precalculated in init_gb() */
 776             sk    = born->param[aj];
 777             raj   = gb_radius[aj];
 778
 779             /* aj -> ai interaction */
 780             if (rai < dr+sk)
 781             {
 782                 lij       = 1.0/(dr-sk);
 783                 dlij      = 1.0;
 784
 785                 if (rai > dr-sk)
 786                 {
 787                     lij  = rai_inv;
 788                     dlij = 0.0;
 789                 }
 790
 791                 uij      = 1.0/(dr+sk);
 792                 lij2     = lij  * lij;
 793                 lij3     = lij2 * lij;
 794                 uij2     = uij  * uij;
 795                 uij3     = uij2 * uij;
 796
 797                 diff2    = uij2-lij2;
 798
 799                 lij_inv  = gmx_invsqrt(lij2);
 800                 sk2      = sk*sk;
 801                 sk2_rinv = sk2*rinv;
 802                 prod     = 0.25*sk2_rinv;
 803
 804                 log_term = log(uij*lij_inv);
 805
 806                 tmp      = lij-uij + 0.25*dr*diff2 + (0.5*rinv)*log_term + prod*(-diff2);
 807
 808                 if (rai < sk-dr)
 809                 {
 810                     tmp = tmp + 2.0 * (rai_inv-lij);
 811                 }
 812
 813                 /* duij    = 1.0; */
 814                 t1      = 0.5*lij2 + prod*lij3 - 0.25*(lij*rinv+lij3*dr);
 815                 t2      = -0.5*uij2 - 0.25*sk2_rinv*uij3 + 0.25*(uij*rinv+uij3*dr);
 816                 t3      = 0.125*(1.0+sk2_rinv*rinv)*(-diff2)+0.25*log_term*rinv*rinv;
 817
 818                 dadx_val = (dlij*t1+t2+t3)*rinv; /* rb2 is moved to chainrule    */
 819
 820                 sum_ai += 0.5*tmp;
 821             }
 822             else
 823             {
 824                 dadx_val = 0.0;
 825             }
 826             fr->dadx[n++] = dadx_val;
 827
 828             /* ai -> aj interaction */
 829             if (raj < dr + sk_ai)
 830             {
 831                 lij     = 1.0/(dr-sk_ai);
 832                 dlij    = 1.0;
 833                 raj_inv = 1.0/raj;
 834
 835                 if (raj > dr-sk_ai)
 836                 {
 837                     lij  = raj_inv;
 838                     dlij = 0.0;
 839                 }
 840
 841                 lij2     = lij  * lij;
 842                 lij3     = lij2 * lij;
 843
 844                 uij      = 1.0/(dr+sk_ai);
 845                 uij2     = uij  * uij;
 846                 uij3     = uij2 * uij;
 847
 848                 diff2    = uij2-lij2;
 849
 850                 lij_inv  = gmx_invsqrt(lij2);
 851                 sk2      =  sk2_ai; /* sk2_ai = sk_ai * sk_ai in i loop above */
 852                 sk2_rinv = sk2*rinv;
 853                 prod     = 0.25 * sk2_rinv;
 854
 855                 /* log_term = table_log(uij*lij_inv,born->log_table,LOG_TABLE_ACCURACY); */
 856                 log_term = log(uij*lij_inv);
 857
 858                 tmp      = lij-uij + 0.25*dr*diff2 + (0.5*rinv)*log_term + prod*(-diff2);
 859
 860                 if (raj < sk_ai-dr)
 861                 {
 862                     tmp     = tmp + 2.0 * (raj_inv-lij);
 863                 }
 864
 865                 t1      = 0.5*lij2 + prod*lij3 - 0.25*(lij*rinv+lij3*dr);
 866                 t2      = -0.5*uij2 - 0.25*sk2_rinv*uij3 + 0.25*(uij*rinv+uij3*dr);
 867                 t3      = 0.125*(1.0+sk2_rinv*rinv)*(-diff2)+0.25*log_term*rinv*rinv;
 868
 869                 dadx_val = (dlij*t1+t2+t3)*rinv; /* rb2 is moved to chainrule    */
 870
 871                 born->gpol_hct_work[aj] += 0.5*tmp;
 872
 873             }
 874             else
 875             {
 876                 dadx_val = 0.0;
 877             }
 878             fr->dadx[n++] = dadx_val;
 879
 880         }
 881         born->gpol_hct_work[ai] += sum_ai;
 882
 883     }
 884
 885     /* Parallel summations */
 886     if (DOMAINDECOMP(cr))
 887     {
 888         dd_atom_sum_real(cr->dd, born->gpol_hct_work);
 889     }
 890
 891     for (i = 0; i < fr->natoms_force; i++) /* PELA born->nr */
 892     {
 893         if (born->use[i] != 0)
 894         {
 895             rai        = top->atomtypes.gb_radius[md->typeA[i]];
 896             rai_inv2   = 1.0/rai;
 897             rai        = rai-doffset;
 898             rai_inv    = 1.0/rai;
 899             sum_ai     = rai * born->gpol_hct_work[i];
 900             sum_ai2    = sum_ai  * sum_ai;
 901             sum_ai3    = sum_ai2 * sum_ai;
 902
 903             tsum          = tanh(born->obc_alpha*sum_ai-born->obc_beta*sum_ai2+born->obc_gamma*sum_ai3);
 904             born->bRad[i] = rai_inv - tsum*rai_inv2;
 905             born->bRad[i] = 1.0 / born->bRad[i];
 906
 907             fr->invsqrta[i] = gmx_invsqrt(born->bRad[i]);
 908
 909             tchain         = rai * (born->obc_alpha-2*born->obc_beta*sum_ai+3*born->obc_gamma*sum_ai2);
 910             born->drobc[i] = (1.0-tsum*tsum)*tchain*rai_inv2;
 911         }
 912     }
 913
 914     /* Extra (local) communication required for DD */
 915     if (DOMAINDECOMP(cr))
 916     {
 917         dd_atom_spread_real(cr->dd, born->bRad);
 918         dd_atom_spread_real(cr->dd, fr->invsqrta);
 919         dd_atom_spread_real(cr->dd, born->drobc);
 920     }
 921
 922     return 0;
 923
 924 }
 925
 926
 927
 928 int calc_gb_rad(t_commrec *cr, t_forcerec *fr, t_inputrec *ir, gmx_localtop_t *top,
 929                 rvec x[], t_nblist *nl, gmx_genborn_t *born, t_mdatoms *md, t_nrnb     *nrnb)
 930 {
 931     real *p;
 932     int   cnt;
 933     int   ndadx;
 934
 935     if (fr->bAllvsAll && fr->dadx == NULL)
 936     {
 937         /* We might need up to 8 atoms of padding before and after,
 938          * and another 4 units to guarantee SSE alignment.
 939          */
 940         fr->nalloc_dadx = 2*(md->homenr+12)*(md->nr/2+1+12);
 941         snew(fr->dadx_rawptr, fr->nalloc_dadx);
 942         fr->dadx = (real *) (((size_t) fr->dadx_rawptr + 16) & (~((size_t) 15)));
 943     }
 944     else
 945     {
 946         /* In the SSE-enabled gb-loops, when writing to dadx, we
 947          * always write 2*4 elements at a time, even in the case with only
 948          * 1-3 j particles, where we only really need to write 2*(1-3)
 949          * elements. This is because we want dadx to be aligned to a 16-
 950          * byte boundary, and being able to use _mm_store/load_ps
 951          */
 952         ndadx = 2 * (nl->nrj + 3*nl->nri);
 953
 954         /* First, reallocate the dadx array, we need 3 extra for SSE */
 955         if (ndadx + 3 > fr->nalloc_dadx)
 956         {
 957             fr->nalloc_dadx = over_alloc_large(ndadx) + 3;
 958             srenew(fr->dadx_rawptr, fr->nalloc_dadx);
 959             fr->dadx = (real *) (((size_t) fr->dadx_rawptr + 16) & (~((size_t) 15)));
 960         }
 961     }
 962
 963     if (fr->bAllvsAll)
 964     {
 965         cnt = md->homenr*(md->nr/2+1);
 966
 967         if (ir->gb_algorithm == egbSTILL)
 968         {
 969             genborn_allvsall_calc_still_radii(fr, md, born, top, x[0], &fr->AllvsAll_workgb);
 970             /* 13 flops in outer loop, 47 flops in inner loop */
 971             inc_nrnb(nrnb, eNR_BORN_AVA_RADII_STILL, md->homenr*13+cnt*47);
 972         }
 973         else if (ir->gb_algorithm == egbHCT || ir->gb_algorithm == egbOBC)
 974         {
 975             genborn_allvsall_calc_hct_obc_radii(fr, md, born, ir->gb_algorithm, top, x[0], &fr->AllvsAll_workgb);
 976             /* 24 flops in outer loop, 183 in inner */
 977             inc_nrnb(nrnb, eNR_BORN_AVA_RADII_HCT_OBC, md->homenr*24+cnt*183);
 978         }
 979         else
 980         {
 981             gmx_fatal(FARGS, "Bad gb algorithm for all-vs-all interactions");
 982         }
 983         return 0;
 984     }
 985
 986     /* Switch for determining which algorithm to use for Born radii calculation */
 987 #ifdef GMX_DOUBLE
 988
 989     switch (ir->gb_algorithm)
 990     {
 991         case egbSTILL:
 992             calc_gb_rad_still(cr, fr, top, x, nl, born, md);
 993             break;
 994         case egbHCT:
 995             calc_gb_rad_hct(cr, fr, top, x, nl, born, md);
 996             break;
 997         case egbOBC:
 998             calc_gb_rad_obc(cr, fr, top, x, nl, born, md);
 999             break;
1000
1001         default:
1002             gmx_fatal(FARGS, "Unknown double precision algorithm for Born radii calculation: %d", ir->gb_algorithm);
1003     }
1004
1005 #else
1006
1007     switch (ir->gb_algorithm)
1008     {
1009         case egbSTILL:
1010             calc_gb_rad_still(cr, fr, top, x, nl, born, md);
1011             break;
1012         case egbHCT:
1013             calc_gb_rad_hct(cr, fr, top, x, nl, born, md);
1014             break;
1015         case egbOBC:
1016             calc_gb_rad_obc(cr, fr, top, x, nl, born, md);
1017             break;
1018
1019         default:
1020             gmx_fatal(FARGS, "Unknown algorithm for Born radii calculation: %d", ir->gb_algorithm);
1021     }
1022
1023 #endif /* Double or single precision */
1024
1025     if (fr->bAllvsAll == FALSE)
1026     {
1027         switch (ir->gb_algorithm)
1028         {
1029             case egbSTILL:
1030                 /* 17 flops per outer loop iteration, 47 flops per inner loop */
1031                 inc_nrnb(nrnb, eNR_BORN_RADII_STILL, nl->nri*17+nl->nrj*47);
1032                 break;
1033             case egbHCT:
1034             case egbOBC:
1035                 /* 61 (assuming 10 for tanh) flops for outer loop iteration, 183 flops per inner loop */
1036                 inc_nrnb(nrnb, eNR_BORN_RADII_HCT_OBC, nl->nri*61+nl->nrj*183);
1037                 break;
1038
1039             default:
1040                 break;
1041         }
1042     }
1043
1044     return 0;
1045 }
1046
1047
1048
1049 real gb_bonds_tab(rvec x[], rvec f[], rvec fshift[], real *charge, real *p_gbtabscale,
1050                   real *invsqrta, real *dvda, real *GBtab, t_idef *idef, real epsilon_r,
1051                   real gb_epsilon_solvent, real facel, const t_pbc *pbc, const t_graph *graph)
1052 {
1053     int      i, j, n0, m, nnn, type, ai, aj;
1054     int      ki;
1055
1056     real     isai, isaj;
1057     real     r, rsq11;
1058     real     rinv11, iq;
1059     real     isaprod, qq, gbscale, gbtabscale, Y, F, Geps, Heps2, Fp, VV, FF, rt, eps, eps2;
1060     real     vgb, fgb, vcoul, fijC, dvdatmp, fscal, dvdaj;
1061     real     vctot;
1062
1063     rvec     dx;
1064     ivec     dt;
1065
1066     t_iatom *forceatoms;
1067
1068     /* Scale the electrostatics by gb_epsilon_solvent */
1069     facel = facel * ((1.0/epsilon_r) - 1.0/gb_epsilon_solvent);
1070
1071     gbtabscale = *p_gbtabscale;
1072     vctot      = 0.0;
1073
1074     for (j = F_GB12; j <= F_GB14; j++)
1075     {
1076         forceatoms = idef->il[j].iatoms;
1077
1078         for (i = 0; i < idef->il[j].nr; )
1079         {
1080             /* To avoid reading in the interaction type, we just increment i to pass over
1081              * the types in the forceatoms array, this saves some memory accesses
1082              */
1083             i++;
1084             ai            = forceatoms[i++];
1085             aj            = forceatoms[i++];
1086
1087             ki            = pbc_rvec_sub(pbc, x[ai], x[aj], dx);
1088             rsq11         = iprod(dx, dx);
1089
1090             isai          = invsqrta[ai];
1091             iq            = (-1)*facel*charge[ai];
1092
1093             rinv11        = gmx_invsqrt(rsq11);
1094             isaj          = invsqrta[aj];
1095             isaprod       = isai*isaj;
1096             qq            = isaprod*iq*charge[aj];
1097             gbscale       = isaprod*gbtabscale;
1098             r             = rsq11*rinv11;
1099             rt            = r*gbscale;
1100             n0            = rt;
1101             eps           = rt-n0;
1102             eps2          = eps*eps;
1103             nnn           = 4*n0;
1104             Y             = GBtab[nnn];
1105             F             = GBtab[nnn+1];
1106             Geps          = eps*GBtab[nnn+2];
1107             Heps2         = eps2*GBtab[nnn+3];
1108             Fp            = F+Geps+Heps2;
1109             VV            = Y+eps*Fp;
1110             FF            = Fp+Geps+2.0*Heps2;
1111             vgb           = qq*VV;
1112             fijC          = qq*FF*gbscale;
1113             dvdatmp       = -(vgb+fijC*r)*0.5;
1114             dvda[aj]      = dvda[aj] + dvdatmp*isaj*isaj;
1115             dvda[ai]      = dvda[ai] + dvdatmp*isai*isai;
1116             vctot         = vctot + vgb;
1117             fgb           = -(fijC)*rinv11;
1118
1119             if (graph)
1120             {
1121                 ivec_sub(SHIFT_IVEC(graph, ai), SHIFT_IVEC(graph, aj), dt);
1122                 ki = IVEC2IS(dt);
1123             }
1124
1125             for (m = 0; (m < DIM); m++)             /*  15              */
1126             {
1127                 fscal               = fgb*dx[m];
1128                 f[ai][m]           += fscal;
1129                 f[aj][m]           -= fscal;
1130                 fshift[ki][m]      += fscal;
1131                 fshift[CENTRAL][m] -= fscal;
1132             }
1133         }
1134     }
1135
1136     return vctot;
1137 }
1138
1139 real calc_gb_selfcorrections(t_commrec *cr, int natoms,
1140                              real *charge, gmx_genborn_t *born, real *dvda, double facel)
1141 {
1142     int  i, ai, at0, at1;
1143     real rai, e, derb, q, q2, fi, rai_inv, vtot;
1144
1145     if (DOMAINDECOMP(cr))
1146     {
1147         at0 = 0;
1148         at1 = cr->dd->nat_home;
1149     }
1150     else
1151     {
1152         at0 = 0;
1153         at1 = natoms;
1154
1155     }
1156
1157     /* Scale the electrostatics by gb_epsilon_solvent */
1158     facel = facel * ((1.0/born->epsilon_r) - 1.0/born->gb_epsilon_solvent);
1159
1160     vtot = 0.0;
1161
1162     /* Apply self corrections */
1163     for (i = at0; i < at1; i++)
1164     {
1165         ai       = i;
1166
1167         if (born->use[ai] == 1)
1168         {
1169             rai       = born->bRad[ai];
1170             rai_inv   = 1.0/rai;
1171             q         = charge[ai];
1172             q2        = q*q;
1173             fi        = facel*q2;
1174             e         = fi*rai_inv;
1175             derb      = 0.5*e*rai_inv*rai_inv;
1176             dvda[ai] += derb*rai;
1177             vtot     -= 0.5*e;
1178         }
1179     }
1180
1181     return vtot;
1182
1183 }
1184
1185 real calc_gb_nonpolar(t_commrec *cr, t_forcerec *fr, int natoms, gmx_genborn_t *born, gmx_localtop_t *top,
1186                       real *dvda, t_mdatoms *md)
1187 {
1188     int  ai, i, at0, at1;
1189     real e, es, rai, rbi, term, probe, tmp, factor;
1190     real rbi_inv, rbi_inv2;
1191
1192     /* To keep the compiler happy */
1193     factor = 0;
1194
1195     if (DOMAINDECOMP(cr))
1196     {
1197         at0 = 0;
1198         at1 = cr->dd->nat_home;
1199     }
1200     else
1201     {
1202         at0 = 0;
1203         at1 = natoms;
1204     }
1205
1206     /* factor is the surface tension */
1207     factor = born->sa_surface_tension;
1208     /*
1209
1210        // The surface tension factor is 0.0049 for Still model, 0.0054 for HCT/OBC
1211        if(gb_algorithm==egbSTILL)
1212        {
1213           factor=0.0049*100*CAL2JOULE;
1214        }
1215        else
1216        {
1217           factor=0.0054*100*CAL2JOULE;
1218        }
1219      */
1220     /* if(gb_algorithm==egbHCT || gb_algorithm==egbOBC) */
1221
1222     es    = 0;
1223     probe = 0.14;
1224     term  = M_PI*4;
1225
1226     for (i = at0; i < at1; i++)
1227     {
1228         ai        = i;
1229
1230         if (born->use[ai] == 1)
1231         {
1232             rai       = top->atomtypes.gb_radius[md->typeA[ai]];
1233             rbi_inv   = fr->invsqrta[ai];
1234             rbi_inv2  = rbi_inv * rbi_inv;
1235             tmp       = (rai*rbi_inv2)*(rai*rbi_inv2);
1236             tmp       = tmp*tmp*tmp;
1237             e         = factor*term*(rai+probe)*(rai+probe)*tmp;
1238             dvda[ai]  = dvda[ai] - 6*e*rbi_inv2;
1239             es        = es + e;
1240         }
1241     }
1242
1243     return es;
1244 }
1245
1246
1247
1248 real calc_gb_chainrule(int natoms, t_nblist *nl, real *dadx, real *dvda, rvec x[], rvec t[], rvec fshift[],
1249                        rvec shift_vec[], int gb_algorithm, gmx_genborn_t *born)
1250 {
1251     int          i, k, n, ai, aj, nj0, nj1, n0, n1;
1252     int          shift;
1253     real         shX, shY, shZ;
1254     real         fgb, fij, rb2, rbi, fix1, fiy1, fiz1;
1255     real         ix1, iy1, iz1, jx1, jy1, jz1, dx11, dy11, dz11, rsq11;
1256     real         rinv11, tx, ty, tz, rbai, rbaj, fgb_ai;
1257     real        *rb;
1258     volatile int idx;
1259
1260     n  = 0;
1261     rb = born->work;
1262
1263     n0 = 0;
1264     n1 = natoms;
1265
1266     if (gb_algorithm == egbSTILL)
1267     {
1268         for (i = n0; i < n1; i++)
1269         {
1270             rbi   = born->bRad[i];
1271             rb[i] = (2 * rbi * rbi * dvda[i])/ONE_4PI_EPS0;
1272         }
1273     }
1274     else if (gb_algorithm == egbHCT)
1275     {
1276         for (i = n0; i < n1; i++)
1277         {
1278             rbi   = born->bRad[i];
1279             rb[i] = rbi * rbi * dvda[i];
1280         }
1281     }
1282     else if (gb_algorithm == egbOBC)
1283     {
1284         for (i = n0; i < n1; i++)
1285         {
1286             rbi   = born->bRad[i];
1287             rb[i] = rbi * rbi * born->drobc[i] * dvda[i];
1288         }
1289     }
1290
1291     for (i = 0; i < nl->nri; i++)
1292     {
1293         ai   = nl->iinr[i];
1294
1295         nj0  = nl->jindex[i];
1296         nj1  = nl->jindex[i+1];
1297
1298         /* Load shifts for this list */
1299         shift   = nl->shift[i];
1300         shX     = shift_vec[shift][0];
1301         shY     = shift_vec[shift][1];
1302         shZ     = shift_vec[shift][2];
1303
1304         /* Load atom i coordinates, add shift vectors */
1305         ix1  = shX + x[ai][0];
1306         iy1  = shY + x[ai][1];
1307         iz1  = shZ + x[ai][2];
1308
1309         fix1 = 0;
1310         fiy1 = 0;
1311         fiz1 = 0;
1312
1313         rbai = rb[ai];
1314
1315         for (k = nj0; k < nj1 && nl->jjnr[k] >= 0; k++)
1316         {
1317             aj = nl->jjnr[k];
1318
1319             jx1     = x[aj][0];
1320             jy1     = x[aj][1];
1321             jz1     = x[aj][2];
1322
1323             dx11    = ix1 - jx1;
1324             dy11    = iy1 - jy1;
1325             dz11    = iz1 - jz1;
1326
1327             rbaj    = rb[aj];
1328
1329             fgb     = rbai*dadx[n++];
1330             fgb_ai  = rbaj*dadx[n++];
1331
1332             /* Total force between ai and aj is the sum of ai->aj and aj->ai */
1333             fgb     = fgb + fgb_ai;
1334
1335             tx      = fgb * dx11;
1336             ty      = fgb * dy11;
1337             tz      = fgb * dz11;
1338
1339             fix1    = fix1 + tx;
1340             fiy1    = fiy1 + ty;
1341             fiz1    = fiz1 + tz;
1342
1343             /* Update force on atom aj */
1344             t[aj][0] = t[aj][0] - tx;
1345             t[aj][1] = t[aj][1] - ty;
1346             t[aj][2] = t[aj][2] - tz;
1347         }
1348
1349         /* Update force and shift forces on atom ai */
1350         t[ai][0] = t[ai][0] + fix1;
1351         t[ai][1] = t[ai][1] + fiy1;
1352         t[ai][2] = t[ai][2] + fiz1;
1353
1354         fshift[shift][0] = fshift[shift][0] + fix1;
1355         fshift[shift][1] = fshift[shift][1] + fiy1;
1356         fshift[shift][2] = fshift[shift][2] + fiz1;
1357
1358     }
1359
1360     return 0;
1361 }
1362
1363
1364 void
1365 calc_gb_forces(t_commrec *cr, t_mdatoms *md, gmx_genborn_t *born, gmx_localtop_t *top,
1366                rvec x[], rvec f[], t_forcerec *fr, t_idef *idef, int gb_algorithm, int sa_algorithm, t_nrnb *nrnb,
1367                const t_pbc *pbc, const t_graph *graph, gmx_enerdata_t *enerd)
1368 {
1369     real v = 0;
1370     int  cnt;
1371     int  i;
1372
1373     /* PBC or not? */
1374     const t_pbc *pbc_null;
1375
1376     if (fr->bMolPBC)
1377     {
1378         pbc_null = pbc;
1379     }
1380     else
1381     {
1382         pbc_null = NULL;
1383     }
1384
1385     if (sa_algorithm == esaAPPROX)
1386     {
1387         /* Do a simple ACE type approximation for the non-polar solvation */
1388         enerd->term[F_NPSOLVATION] += calc_gb_nonpolar(cr, fr, born->nr, born, top, fr->dvda, md);
1389     }
1390
1391     /* Calculate the bonded GB-interactions using either table or analytical formula */
1392     enerd->term[F_GBPOL]       += gb_bonds_tab(x, f, fr->fshift, md->chargeA, &(fr->gbtabscale),
1393                                                fr->invsqrta, fr->dvda, fr->gbtab.data, idef, born->epsilon_r, born->gb_epsilon_solvent, fr->epsfac, pbc_null, graph);
1394
1395     /* Calculate self corrections to the GB energies - currently only A state used! (FIXME) */
1396     enerd->term[F_GBPOL]       += calc_gb_selfcorrections(cr, born->nr, md->chargeA, born, fr->dvda, fr->epsfac);
1397
1398     /* If parallel, sum the derivative of the potential w.r.t the born radii */
1399     if (DOMAINDECOMP(cr))
1400     {
1401         dd_atom_sum_real(cr->dd, fr->dvda);
1402         dd_atom_spread_real(cr->dd, fr->dvda);
1403     }
1404
1405     if (fr->bAllvsAll)
1406     {
1407         genborn_allvsall_calc_chainrule(fr, md, born, x[0], f[0], gb_algorithm, fr->AllvsAll_workgb);
1408         cnt = md->homenr*(md->nr/2+1);
1409         /* 9 flops for outer loop, 15 for inner */
1410         inc_nrnb(nrnb, eNR_BORN_AVA_CHAINRULE, md->homenr*9+cnt*15);
1411         return;
1412     }
1413
1414     calc_gb_chainrule(fr->natoms_force, &(fr->gblist), fr->dadx, fr->dvda,
1415                       x, f, fr->fshift, fr->shift_vec, gb_algorithm, born);
1416
1417     if (!fr->bAllvsAll)
1418     {
1419         /* 9 flops for outer loop, 15 for inner */
1420         inc_nrnb(nrnb, eNR_BORN_CHAINRULE, fr->gblist.nri*9+fr->gblist.nrj*15);
1421     }
1422 }
1423
1424 static void add_j_to_gblist(gbtmpnbl_t *list, int aj)
1425 {
1426     if (list->naj >= list->aj_nalloc)
1427     {
1428         list->aj_nalloc = over_alloc_large(list->naj+1);
1429         srenew(list->aj, list->aj_nalloc);
1430     }
1431
1432     list->aj[list->naj++] = aj;
1433 }
1434
1435 static gbtmpnbl_t *find_gbtmplist(struct gbtmpnbls *lists, int shift)
1436 {
1437     int ind, i;
1438
1439     /* Search the list with the same shift, if there is one */
1440     ind = 0;
1441     while (ind < lists->nlist && shift != lists->list[ind].shift)
1442     {
1443         ind++;
1444     }
1445     if (ind == lists->nlist)
1446     {
1447         if (lists->nlist == lists->list_nalloc)
1448         {
1449             lists->list_nalloc++;
1450             srenew(lists->list, lists->list_nalloc);
1451             for (i = lists->nlist; i < lists->list_nalloc; i++)
1452             {
1453                 lists->list[i].aj        = NULL;
1454                 lists->list[i].aj_nalloc = 0;
1455             }
1456
1457         }
1458
1459         lists->list[lists->nlist].shift = shift;
1460         lists->list[lists->nlist].naj   = 0;
1461         lists->nlist++;
1462     }
1463
1464     return &lists->list[ind];
1465 }
1466
1467 static void add_bondeds_to_gblist(t_ilist *il,
1468                                   gmx_bool bMolPBC, t_pbc *pbc, t_graph *g, rvec *x,
1469                                   struct gbtmpnbls *nls)
1470 {
1471     int         ind, j, ai, aj, shift, found;
1472     rvec        dx;
1473     ivec        dt;
1474     gbtmpnbl_t *list;
1475
1476     shift = CENTRAL;
1477     for (ind = 0; ind < il->nr; ind += 3)
1478     {
1479         ai = il->iatoms[ind+1];
1480         aj = il->iatoms[ind+2];
1481
1482         shift = CENTRAL;
1483         if (g != NULL)
1484         {
1485             rvec_sub(x[ai], x[aj], dx);
1486             ivec_sub(SHIFT_IVEC(g, ai), SHIFT_IVEC(g, aj), dt);
1487             shift = IVEC2IS(dt);
1488         }
1489         else if (bMolPBC)
1490         {
1491             shift = pbc_dx_aiuc(pbc, x[ai], x[aj], dx);
1492         }
1493
1494         /* Find the list for this shift or create one */
1495         list = find_gbtmplist(&nls[ai], shift);
1496
1497         found = 0;
1498
1499         /* So that we do not add the same bond twice.
1500          * This happens with some constraints between 1-3 atoms
1501          * that are in the bond-list but should not be in the GB nb-list */
1502         for (j = 0; j < list->naj; j++)
1503         {
1504             if (list->aj[j] == aj)
1505             {
1506                 found = 1;
1507             }
1508         }
1509
1510         if (found == 0)
1511         {
1512             if (ai == aj)
1513             {
1514                 gmx_incons("ai == aj");
1515             }
1516
1517             add_j_to_gblist(list, aj);
1518         }
1519     }
1520 }
1521
1522 static int
1523 compare_int (const void * a, const void * b)
1524 {
1525     return ( *(int*)a - *(int*)b );
1526 }
1527
1528
1529
1530 int make_gb_nblist(t_commrec *cr, int gb_algorithm,
1531                    rvec x[], matrix box,
1532                    t_forcerec *fr, t_idef *idef, t_graph *graph, gmx_genborn_t *born)
1533 {
1534     int               i, l, ii, j, k, n, nj0, nj1, ai, aj, at0, at1, found, shift, s;
1535     int               apa;
1536     t_nblist         *nblist;
1537     t_pbc             pbc;
1538
1539     struct gbtmpnbls *nls;
1540     gbtmpnbl_t       *list = NULL;
1541
1542     set_pbc(&pbc, fr->ePBC, box);
1543     nls   = born->nblist_work;
1544
1545     for (i = 0; i < born->nr; i++)
1546     {
1547         nls[i].nlist = 0;
1548     }
1549
1550     if (fr->bMolPBC)
1551     {
1552         set_pbc_dd(&pbc, fr->ePBC, cr->dd, TRUE, box);
1553     }
1554
1555     switch (gb_algorithm)
1556     {
1557         case egbHCT:
1558         case egbOBC:
1559             /* Loop over 1-2, 1-3 and 1-4 interactions */
1560             for (j = F_GB12; j <= F_GB14; j++)
1561             {
1562                 add_bondeds_to_gblist(&idef->il[j], fr->bMolPBC, &pbc, graph, x, nls);
1563             }
1564             break;
1565         case egbSTILL:
1566             /* Loop over 1-4 interactions */
1567             add_bondeds_to_gblist(&idef->il[F_GB14], fr->bMolPBC, &pbc, graph, x, nls);
1568             break;
1569         default:
1570             gmx_incons("Unknown GB algorithm");
1571     }
1572
1573     /* Loop over the VDWQQ and VDW nblists to set up the nonbonded part of the GB list */
1574     for (n = 0; (n < fr->nnblists); n++)
1575     {
1576         for (i = 0; (i < eNL_NR); i++)
1577         {
1578             nblist = &(fr->nblists[n].nlist_sr[i]);
1579
1580             if (nblist->nri > 0 && (i == eNL_VDWQQ || i == eNL_QQ))
1581             {
1582                 for (j = 0; j < nblist->nri; j++)
1583                 {
1584                     ai    = nblist->iinr[j];
1585                     shift = nblist->shift[j];
1586
1587                     /* Find the list for this shift or create one */
1588                     list = find_gbtmplist(&nls[ai], shift);
1589
1590                     nj0 = nblist->jindex[j];
1591                     nj1 = nblist->jindex[j+1];
1592
1593                     /* Add all the j-atoms in the non-bonded list to the GB list */
1594                     for (k = nj0; k < nj1; k++)
1595                     {
1596                         add_j_to_gblist(list, nblist->jjnr[k]);
1597                     }
1598                 }
1599             }
1600         }
1601     }
1602
1603     /* Zero out some counters */
1604     fr->gblist.nri = 0;
1605     fr->gblist.nrj = 0;
1606
1607     fr->gblist.jindex[0] = fr->gblist.nri;
1608
1609     for (i = 0; i < fr->natoms_force; i++)
1610     {
1611         for (s = 0; s < nls[i].nlist; s++)
1612         {
1613             list = &nls[i].list[s];
1614
1615             /* Only add those atoms that actually have neighbours */
1616             if (born->use[i] != 0)
1617             {
1618                 fr->gblist.iinr[fr->gblist.nri]  = i;
1619                 fr->gblist.shift[fr->gblist.nri] = list->shift;
1620                 fr->gblist.nri++;
1621
1622                 for (k = 0; k < list->naj; k++)
1623                 {
1624                     /* Memory allocation for jjnr */
1625                     if (fr->gblist.nrj >= fr->gblist.maxnrj)
1626                     {
1627                         fr->gblist.maxnrj += over_alloc_large(fr->gblist.maxnrj);
1628
1629                         if (debug)
1630                         {
1631                             fprintf(debug, "Increasing GB neighbourlist j size to %d\n", fr->gblist.maxnrj);
1632                         }
1633
1634                         srenew(fr->gblist.jjnr, fr->gblist.maxnrj);
1635                     }
1636
1637                     /* Put in list */
1638                     if (i == list->aj[k])
1639                     {
1640                         gmx_incons("i == list->aj[k]");
1641                     }
1642                     fr->gblist.jjnr[fr->gblist.nrj++] = list->aj[k];
1643                 }
1644
1645                 fr->gblist.jindex[fr->gblist.nri] = fr->gblist.nrj;
1646             }
1647         }
1648     }
1649
1650
1651 #ifdef SORT_GB_LIST
1652     for (i = 0; i < fr->gblist.nri; i++)
1653     {
1654         nj0 = fr->gblist.jindex[i];
1655         nj1 = fr->gblist.jindex[i+1];
1656         ai  = fr->gblist.iinr[i];
1657
1658         /* Temporary fix */
1659         for (j = nj0; j < nj1; j++)
1660         {
1661             if (fr->gblist.jjnr[j] < ai)
1662             {
1663                 fr->gblist.jjnr[j] += fr->natoms_force;
1664             }
1665         }
1666         qsort(fr->gblist.jjnr+nj0, nj1-nj0, sizeof(int), compare_int);
1667         /* Fix back */
1668         for (j = nj0; j < nj1; j++)
1669         {
1670             if (fr->gblist.jjnr[j] >= fr->natoms_force)
1671             {
1672                 fr->gblist.jjnr[j] -= fr->natoms_force;
1673             }
1674         }
1675
1676     }
1677 #endif
1678
1679     return 0;
1680 }
1681
1682 void make_local_gb(const t_commrec *cr, gmx_genborn_t *born, int gb_algorithm)
1683 {
1684     int           i, at0, at1;
1685     gmx_domdec_t *dd = NULL;
1686
1687     if (DOMAINDECOMP(cr))
1688     {
1689         dd  = cr->dd;
1690         at0 = 0;
1691         at1 = dd->nat_tot;
1692     }
1693     else
1694     {
1695         /* Single node, just copy pointers and return */
1696         if (gb_algorithm == egbSTILL)
1697         {
1698             born->gpol      = born->gpol_globalindex;
1699             born->vsolv     = born->vsolv_globalindex;
1700             born->gb_radius = born->gb_radius_globalindex;
1701         }
1702         else
1703         {
1704             born->param     = born->param_globalindex;
1705             born->gb_radius = born->gb_radius_globalindex;
1706         }
1707
1708         born->use = born->use_globalindex;
1709
1710         return;
1711     }
1712
1713     /* Reallocation of local arrays if necessary */
1714     /* fr->natoms_force is equal to dd->nat_tot */
1715     if (DOMAINDECOMP(cr) && dd->nat_tot > born->nalloc)
1716     {
1717         int nalloc;
1718
1719         nalloc = dd->nat_tot;
1720
1721         /* Arrays specific to different gb algorithms */
1722         if (gb_algorithm == egbSTILL)
1723         {
1724             srenew(born->gpol,  nalloc+3);
1725             srenew(born->vsolv, nalloc+3);
1726             srenew(born->gb_radius, nalloc+3);
1727             for (i = born->nalloc; (i < nalloc+3); i++)
1728             {
1729                 born->gpol[i]      = 0;
1730                 born->vsolv[i]     = 0;
1731                 born->gb_radius[i] = 0;
1732             }
1733         }
1734         else
1735         {
1736             srenew(born->param, nalloc+3);
1737             srenew(born->gb_radius, nalloc+3);
1738             for (i = born->nalloc; (i < nalloc+3); i++)
1739             {
1740                 born->param[i]     = 0;
1741                 born->gb_radius[i] = 0;
1742             }
1743         }
1744
1745         /* All gb-algorithms use the array for vsites exclusions */
1746         srenew(born->use,    nalloc+3);
1747         for (i = born->nalloc; (i < nalloc+3); i++)
1748         {
1749             born->use[i] = 0;
1750         }
1751
1752         born->nalloc = nalloc;
1753     }
1754
1755     /* With dd, copy algorithm specific arrays */
1756     if (gb_algorithm == egbSTILL)
1757     {
1758         for (i = at0; i < at1; i++)
1759         {
1760             born->gpol[i]      = born->gpol_globalindex[dd->gatindex[i]];
1761             born->vsolv[i]     = born->vsolv_globalindex[dd->gatindex[i]];
1762             born->gb_radius[i] = born->gb_radius_globalindex[dd->gatindex[i]];
1763             born->use[i]       = born->use_globalindex[dd->gatindex[i]];
1764         }
1765     }
1766     else
1767     {
1768         for (i = at0; i < at1; i++)
1769         {
1770             born->param[i]     = born->param_globalindex[dd->gatindex[i]];
1771             born->gb_radius[i] = born->gb_radius_globalindex[dd->gatindex[i]];
1772             born->use[i]       = born->use_globalindex[dd->gatindex[i]];
1773         }
1774     }
1775 }