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