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