src/gromacs/mdlib/genborn.cpp

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