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