src/gromacs/mdlib/genborn.c

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