src/gromacs/mdlib/genborn.c

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