src/gromacs/mdlib/genborn.c

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