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                 return pbc_dx_aiuc(pbc,xi,xj,dx);
 102         }
 103         else {
 104                 rvec_sub(xi,xj,dx);
 105                 return CENTRAL;
 106         }
 107 }
 108
 109 int init_gb_nblist(int natoms, t_nblist *nl)
 110 {
 111     nl->maxnri      = natoms*4;
 112     nl->maxnrj      = 0;
 113     nl->maxlen      = 0;
 114     nl->nri         = 0;
 115     nl->nrj         = 0;
 116     nl->iinr        = NULL;
 117     nl->gid         = NULL;
 118     nl->shift       = NULL;
 119     nl->jindex      = NULL;
 120     nl->jjnr        = NULL;
 121     /*nl->nltype      = nltype;*/
 122
 123     srenew(nl->iinr,   nl->maxnri);
 124     srenew(nl->gid,    nl->maxnri);
 125     srenew(nl->shift,  nl->maxnri);
 126     srenew(nl->jindex, nl->maxnri+1);
 127
 128     nl->jindex[0] = 0;
 129
 130     return 0;
 131 }
 132
 133 void gb_pd_send(t_commrec *cr, real *send_data, int nr)
 134 {
 135 #ifdef GMX_MPI
 136     int i,cur;
 137     int *index,*sendc,*disp;
 138
 139     snew(sendc,cr->nnodes);
 140     snew(disp,cr->nnodes);
 141
 142     index = pd_index(cr);
 143     cur   = cr->nodeid;
 144
 145     /* Setup count/index arrays */
 146     for(i=0;i<cr->nnodes;i++)
 147     {
 148         sendc[i]  = index[i+1]-index[i];
 149         disp[i]   = index[i];
 150     }
 151
 152     /* Do communication */
 153     MPI_Gatherv(send_data+index[cur],sendc[cur],GMX_MPI_REAL,send_data,sendc,
 154                 disp,GMX_MPI_REAL,0,cr->mpi_comm_mygroup);
 155     MPI_Bcast(send_data,nr,GMX_MPI_REAL,0,cr->mpi_comm_mygroup);
 156
 157 #endif
 158 }
 159
 160
 161 int init_gb_plist(t_params *p_list)
 162 {
 163     p_list->nr    = 0;
 164     p_list->param = NULL;
 165
 166     return 0;
 167 }
 168
 169
 170
 171 int init_gb_still(const t_commrec *cr, t_forcerec  *fr,
 172                   const t_atomtypes *atype, t_idef *idef, t_atoms *atoms,
 173                   gmx_genborn_t *born,int natoms)
 174 {
 175
 176     int i,j,i1,i2,k,m,nbond,nang,ia,ib,ic,id,nb,idx,idx2,at;
 177     int iam,ibm;
 178     int at0,at1;
 179     real length,angle;
 180     real r,ri,rj,ri2,ri3,rj2,r2,r3,r4,rk,ratio,term,h,doffset;
 181     real p1,p2,p3,factor,cosine,rab,rbc;
 182
 183     real *vsol;
 184     real *gp;
 185
 186     snew(vsol,natoms);
 187     snew(gp,natoms);
 188     snew(born->gpol_still_work,natoms+3);
 189
 190     if(PAR(cr))
 191     {
 192         if(PARTDECOMP(cr))
 193         {
 194             pd_at_range(cr,&at0,&at1);
 195
 196             for(i=0;i<natoms;i++)
 197             {
 198                 vsol[i] = gp[i] = 0;
 199             }
 200         }
 201         else
 202         {
 203             at0 = 0;
 204             at1 = natoms;
 205         }
 206     }
 207     else
 208     {
 209         at0 = 0;
 210         at1 = natoms;
 211     }
 212
 213     doffset = born->gb_doffset;
 214
 215     for(i=0;i<natoms;i++)
 216     {
 217         born->gpol_globalindex[i]=born->vsolv_globalindex[i]=
 218             born->gb_radius_globalindex[i]=0;
 219     }
 220
 221     /* Compute atomic solvation volumes for Still method */
 222     for(i=0;i<natoms;i++)
 223     {
 224         ri=atype->gb_radius[atoms->atom[i].type];
 225         born->gb_radius_globalindex[i] = ri;
 226         r3=ri*ri*ri;
 227         born->vsolv_globalindex[i]=(4*M_PI/3)*r3;
 228     }
 229
 230     for(j=0;j<idef->il[F_GB12].nr;j+=3)
 231     {
 232         m=idef->il[F_GB12].iatoms[j];
 233         ia=idef->il[F_GB12].iatoms[j+1];
 234         ib=idef->il[F_GB12].iatoms[j+2];
 235
 236         r=1.01*idef->iparams[m].gb.st;
 237
 238         ri   = atype->gb_radius[atoms->atom[ia].type];
 239         rj   = atype->gb_radius[atoms->atom[ib].type];
 240
 241         ri2  = ri*ri;
 242         ri3  = ri2*ri;
 243         rj2  = rj*rj;
 244
 245         ratio  = (rj2-ri2-r*r)/(2*ri*r);
 246         h      = ri*(1+ratio);
 247         term   = (M_PI/3.0)*h*h*(3.0*ri-h);
 248
 249         if(PARTDECOMP(cr))
 250         {
 251             vsol[ia]+=term;
 252         }
 253         else
 254         {
 255             born->vsolv_globalindex[ia] -= term;
 256         }
 257
 258         ratio  = (ri2-rj2-r*r)/(2*rj*r);
 259         h      = rj*(1+ratio);
 260         term   = (M_PI/3.0)*h*h*(3.0*rj-h);
 261
 262         if(PARTDECOMP(cr))
 263         {
 264             vsol[ib]+=term;
 265         }
 266         else
 267         {
 268             born->vsolv_globalindex[ib] -= term;
 269         }
 270     }
 271
 272     if(PARTDECOMP(cr))
 273     {
 274         gmx_sum(natoms,vsol,cr);
 275
 276         for(i=0;i<natoms;i++)
 277         {
 278             born->vsolv_globalindex[i]=born->vsolv_globalindex[i]-vsol[i];
 279         }
 280     }
 281
 282     /* Get the self-, 1-2 and 1-3 polarization energies for analytical Still
 283        method */
 284     /* Self */
 285     for(j=0;j<natoms;j++)
 286     {
 287         if(born->use_globalindex[j]==1)
 288         {
 289             born->gpol_globalindex[j]=-0.5*ONE_4PI_EPS0/
 290                 (atype->gb_radius[atoms->atom[j].type]-doffset+STILL_P1);
 291         }
 292     }
 293
 294     /* 1-2 */
 295     for(j=0;j<idef->il[F_GB12].nr;j+=3)
 296     {
 297         m=idef->il[F_GB12].iatoms[j];
 298         ia=idef->il[F_GB12].iatoms[j+1];
 299         ib=idef->il[F_GB12].iatoms[j+2];
 300
 301         r=idef->iparams[m].gb.st;
 302
 303         r4=r*r*r*r;
 304
 305         if(PARTDECOMP(cr))
 306         {
 307             gp[ia]+=STILL_P2*born->vsolv_globalindex[ib]/r4;
 308             gp[ib]+=STILL_P2*born->vsolv_globalindex[ia]/r4;
 309         }
 310         else
 311         {
 312             born->gpol_globalindex[ia]=born->gpol_globalindex[ia]+
 313                 STILL_P2*born->vsolv_globalindex[ib]/r4;
 314             born->gpol_globalindex[ib]=born->gpol_globalindex[ib]+
 315                 STILL_P2*born->vsolv_globalindex[ia]/r4;
 316         }
 317     }
 318
 319     /* 1-3 */
 320     for(j=0;j<idef->il[F_GB13].nr;j+=3)
 321     {
 322         m=idef->il[F_GB13].iatoms[j];
 323         ia=idef->il[F_GB13].iatoms[j+1];
 324         ib=idef->il[F_GB13].iatoms[j+2];
 325
 326         r=idef->iparams[m].gb.st;
 327         r4=r*r*r*r;
 328
 329         if(PARTDECOMP(cr))
 330         {
 331             gp[ia]+=STILL_P3*born->vsolv[ib]/r4;
 332             gp[ib]+=STILL_P3*born->vsolv[ia]/r4;
 333         }
 334         else
 335         {
 336             born->gpol_globalindex[ia]=born->gpol_globalindex[ia]+
 337                 STILL_P3*born->vsolv_globalindex[ib]/r4;
 338             born->gpol_globalindex[ib]=born->gpol_globalindex[ib]+
 339                 STILL_P3*born->vsolv_globalindex[ia]/r4;
 340         }
 341     }
 342
 343     if(PARTDECOMP(cr))
 344     {
 345         gmx_sum(natoms,gp,cr);
 346
 347         for(i=0;i<natoms;i++)
 348         {
 349             born->gpol_globalindex[i]=born->gpol_globalindex[i]+gp[i];
 350         }
 351     }
 352
 353     sfree(vsol);
 354     sfree(gp);
 355
 356     return 0;
 357 }
 358
 359 /* Initialize all GB datastructs and compute polarization energies */
 360 int init_gb(gmx_genborn_t **p_born,
 361             const t_commrec *cr, t_forcerec *fr, const t_inputrec *ir,
 362             const gmx_mtop_t *mtop, real rgbradii, int gb_algorithm)
 363 {
 364     int i,j,m,ai,aj,jj,natoms,nalloc;
 365     real rai,sk,p,doffset;
 366
 367     t_atoms        atoms;
 368     gmx_genborn_t  *born;
 369     gmx_localtop_t *localtop;
 370
 371     natoms   = mtop->natoms;
 372
 373     atoms    = gmx_mtop_global_atoms(mtop);
 374     localtop = gmx_mtop_generate_local_top(mtop,ir);
 375
 376     snew(born,1);
 377     *p_born = born;
 378
 379     born->nr  = natoms;
 380
 381     snew(born->drobc, natoms);
 382     snew(born->bRad,  natoms);
 383
 384     /* Allocate memory for the global data arrays */
 385     snew(born->param_globalindex, natoms+3);
 386     snew(born->gpol_globalindex,  natoms+3);
 387     snew(born->vsolv_globalindex, natoms+3);
 388     snew(born->gb_radius_globalindex, natoms+3);
 389     snew(born->use_globalindex,    natoms+3);
 390
 391     snew(fr->invsqrta, natoms);
 392     snew(fr->dvda,     natoms);
 393
 394     fr->dadx              = NULL;
 395     fr->dadx_rawptr       = NULL;
 396     fr->nalloc_dadx       = 0;
 397     born->gpol_still_work = NULL;
 398     born->gpol_hct_work   = NULL;
 399
 400     /* snew(born->asurf,natoms); */
 401     /* snew(born->dasurf,natoms); */
 402
 403     /* Initialize the gb neighbourlist */
 404     init_gb_nblist(natoms,&(fr->gblist));
 405
 406     /* Do the Vsites exclusions (if any) */
 407     for(i=0;i<natoms;i++)
 408     {
 409         jj = atoms.atom[i].type;
 410         if (mtop->atomtypes.gb_radius[atoms.atom[i].type] > 0)
 411         {
 412             born->use_globalindex[i] = 1;
 413         }
 414         else
 415         {
 416             born->use_globalindex[i] = 0;
 417         }
 418
 419         /* If we have a Vsite, put vs_globalindex[i]=0 */
 420         if (C6 (fr->nbfp,fr->ntype,jj,jj) == 0 &&
 421             C12(fr->nbfp,fr->ntype,jj,jj) == 0 &&
 422             atoms.atom[i].q == 0)
 423         {
 424             born->use_globalindex[i]=0;
 425         }
 426     }
 427
 428     /* Copy algorithm parameters from inputrecord to local structure */
 429     born->obc_alpha  = ir->gb_obc_alpha;
 430     born->obc_beta   = ir->gb_obc_beta;
 431     born->obc_gamma  = ir->gb_obc_gamma;
 432     born->gb_doffset = ir->gb_dielectric_offset;
 433     born->gb_epsilon_solvent = ir->gb_epsilon_solvent;
 434     born->epsilon_r = ir->epsilon_r;
 435
 436     doffset = born->gb_doffset;
 437
 438     /* Set the surface tension */
 439     born->sa_surface_tension = ir->sa_surface_tension;
 440
 441     /* If Still model, initialise the polarisation energies */
 442     if(gb_algorithm==egbSTILL)
 443     {
 444         init_gb_still(cr, fr,&(mtop->atomtypes), &(localtop->idef), &atoms,
 445                       born, natoms);
 446     }
 447
 448
 449     /* If HCT/OBC,  precalculate the sk*atype->S_hct factors */
 450     else if(gb_algorithm==egbHCT || gb_algorithm==egbOBC)
 451     {
 452
 453         snew(born->gpol_hct_work, natoms+3);
 454
 455         for(i=0;i<natoms;i++)
 456         {
 457             if(born->use_globalindex[i]==1)
 458             {
 459                 rai = mtop->atomtypes.gb_radius[atoms.atom[i].type]-doffset;
 460                 sk  = rai * mtop->atomtypes.S_hct[atoms.atom[i].type];
 461                 born->param_globalindex[i] = sk;
 462                 born->gb_radius_globalindex[i] = rai;
 463             }
 464             else
 465             {
 466                 born->param_globalindex[i] = 0;
 467                 born->gb_radius_globalindex[i] = 0;
 468             }
 469         }
 470     }
 471
 472     /* Allocate memory for work arrays for temporary use */
 473     snew(born->work,natoms+4);
 474     snew(born->count,natoms);
 475     snew(born->nblist_work,natoms);
 476
 477     /* Domain decomposition specific stuff */
 478     born->nalloc = 0;
 479
 480     return 0;
 481 }
 482
 483
 484
 485 static int
 486 calc_gb_rad_still(t_commrec *cr, t_forcerec *fr,int natoms, gmx_localtop_t *top,
 487                   const t_atomtypes *atype, rvec x[], t_nblist *nl,
 488                   gmx_genborn_t *born,t_mdatoms *md)
 489 {
 490     int i,k,n,nj0,nj1,ai,aj,type;
 491     int shift;
 492     real shX,shY,shZ;
 493     real gpi,dr,dr2,dr4,idr4,rvdw,ratio,ccf,theta,term,rai,raj;
 494     real ix1,iy1,iz1,jx1,jy1,jz1,dx11,dy11,dz11;
 495     real rinv,idr2,idr6,vaj,dccf,cosq,sinq,prod,gpi2;
 496     real factor;
 497     real vai, prod_ai, icf4,icf6;
 498
 499     factor  = 0.5*ONE_4PI_EPS0;
 500     n       = 0;
 501
 502     for(i=0;i<born->nr;i++)
 503     {
 504         born->gpol_still_work[i]=0;
 505     }
 506
 507         for(i=0;i<nl->nri;i++ )
 508     {
 509         ai      = nl->iinr[i];
 510
 511         nj0     = nl->jindex[i];
 512         nj1     = nl->jindex[i+1];
 513
 514         /* Load shifts for this list */
 515         shift   = nl->shift[i];
 516         shX     = fr->shift_vec[shift][0];
 517         shY     = fr->shift_vec[shift][1];
 518         shZ     = fr->shift_vec[shift][2];
 519
 520         gpi     = 0;
 521
 522         rai     = top->atomtypes.gb_radius[md->typeA[ai]];
 523         vai     = born->vsolv[ai];
 524         prod_ai = STILL_P4*vai;
 525
 526         /* Load atom i coordinates, add shift vectors */
 527         ix1     = shX + x[ai][0];
 528         iy1     = shY + x[ai][1];
 529         iz1     = shZ + x[ai][2];
 530
 531         for(k=nj0;k<nj1;k++)
 532         {
 533             aj    = nl->jjnr[k];
 534             jx1   = x[aj][0];
 535             jy1   = x[aj][1];
 536             jz1   = x[aj][2];
 537
 538             dx11  = ix1-jx1;
 539             dy11  = iy1-jy1;
 540             dz11  = iz1-jz1;
 541
 542             dr2   = dx11*dx11+dy11*dy11+dz11*dz11;
 543             rinv  = gmx_invsqrt(dr2);
 544             idr2  = rinv*rinv;
 545             idr4  = idr2*idr2;
 546             idr6  = idr4*idr2;
 547
 548             raj = top->atomtypes.gb_radius[md->typeA[aj]];
 549
 550             rvdw  = rai + raj;
 551
 552             ratio = dr2 / (rvdw * rvdw);
 553             vaj   = born->vsolv[aj];
 554
 555             if(ratio>STILL_P5INV)
 556             {
 557                 ccf=1.0;
 558                 dccf=0.0;
 559             }
 560             else
 561             {
 562                 theta = ratio*STILL_PIP5;
 563                 cosq  = cos(theta);
 564                 term  = 0.5*(1.0-cosq);
 565                 ccf   = term*term;
 566                 sinq  = 1.0 - cosq*cosq;
 567                 dccf  = 2.0*term*sinq*gmx_invsqrt(sinq)*theta;
 568             }
 569
 570             prod          = STILL_P4*vaj;
 571             icf4          = ccf*idr4;
 572             icf6          = (4*ccf-dccf)*idr6;
 573
 574             born->gpol_still_work[aj] += prod_ai*icf4;
 575             gpi             = gpi+prod*icf4;
 576
 577             /* Save ai->aj and aj->ai chain rule terms */
 578             fr->dadx[n++]   = prod*icf6;
 579             fr->dadx[n++]   = prod_ai*icf6;
 580         }
 581         born->gpol_still_work[ai]+=gpi;
 582     }
 583
 584     /* Parallel summations */
 585     if(PARTDECOMP(cr))
 586     {
 587         gmx_sum(natoms, born->gpol_still_work, cr);
 588     }
 589     else if(DOMAINDECOMP(cr))
 590     {
 591         dd_atom_sum_real(cr->dd, born->gpol_still_work);
 592     }
 593
 594     /* Calculate the radii */
 595         for(i=0;i<fr->natoms_force;i++) /* PELA born->nr */
 596     {
 597                 if(born->use[i] != 0)
 598         {
 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;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;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             inc_nrnb(nrnb,eNR_BORN_AVA_RADII_STILL,cnt);
1124         }
1125         else if(ir->gb_algorithm==egbHCT || ir->gb_algorithm==egbOBC)
1126         {
1127 #if 0 && defined (GMX_X86_SSE2)
1128             if(fr->use_acceleration)
1129             {
1130 #  ifdef GMX_DOUBLE
1131                 genborn_allvsall_calc_hct_obc_radii_sse2_double(fr,md,born,ir->gb_algorithm,top,x[0],cr,&fr->AllvsAll_workgb);
1132 #  else
1133                 genborn_allvsall_calc_hct_obc_radii_sse2_single(fr,md,born,ir->gb_algorithm,top,x[0],cr,&fr->AllvsAll_workgb);
1134 #  endif
1135             }
1136             else
1137             {
1138                 genborn_allvsall_calc_hct_obc_radii(fr,md,born,ir->gb_algorithm,top,x[0],cr,&fr->AllvsAll_workgb);
1139             }
1140 #else
1141             genborn_allvsall_calc_hct_obc_radii(fr,md,born,ir->gb_algorithm,top,x[0],cr,&fr->AllvsAll_workgb);
1142 #endif
1143             inc_nrnb(nrnb,eNR_BORN_AVA_RADII_HCT_OBC,cnt);
1144         }
1145         else
1146         {
1147             gmx_fatal(FARGS,"Bad gb algorithm for all-vs-all interactions");
1148         }
1149         inc_nrnb(nrnb,eNR_NBKERNEL_OUTER,md->homenr);
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                 inc_nrnb(nrnb,eNR_BORN_RADII_STILL,nl->nrj);
1283                 break;
1284             case egbHCT:
1285             case egbOBC:
1286                 inc_nrnb(nrnb,eNR_BORN_RADII_HCT_OBC,nl->nrj);
1287                 break;
1288
1289             default:
1290                 break;
1291         }
1292         inc_nrnb(nrnb,eNR_NBKERNEL_OUTER,nl->nri);
1293     }
1294
1295     return 0;
1296 }
1297
1298
1299
1300 real gb_bonds_tab(rvec x[], rvec f[], rvec fshift[], real *charge, real *p_gbtabscale,
1301                   real *invsqrta, real *dvda, real *GBtab, t_idef *idef, real epsilon_r,
1302                   real gb_epsilon_solvent, real facel, const t_pbc *pbc, const t_graph *graph)
1303 {
1304     int i,j,n0,m,nnn,type,ai,aj;
1305         int ki;
1306
1307         real isai,isaj;
1308     real r,rsq11;
1309     real rinv11,iq;
1310     real isaprod,qq,gbscale,gbtabscale,Y,F,Geps,Heps2,Fp,VV,FF,rt,eps,eps2;
1311     real vgb,fgb,vcoul,fijC,dvdatmp,fscal,dvdaj;
1312     real vctot;
1313
1314         rvec dx;
1315         ivec dt;
1316
1317     t_iatom *forceatoms;
1318
1319     /* Scale the electrostatics by gb_epsilon_solvent */
1320     facel = facel * ((1.0/epsilon_r) - 1.0/gb_epsilon_solvent);
1321
1322     gbtabscale=*p_gbtabscale;
1323     vctot = 0.0;
1324
1325     for(j=F_GB12;j<=F_GB14;j++)
1326     {
1327         forceatoms = idef->il[j].iatoms;
1328
1329         for(i=0;i<idef->il[j].nr; )
1330         {
1331             /* To avoid reading in the interaction type, we just increment i to pass over
1332              * the types in the forceatoms array, this saves some memory accesses
1333              */
1334             i++;
1335             ai            = forceatoms[i++];
1336             aj            = forceatoms[i++];
1337
1338                         ki            = pbc_rvec_sub(pbc,x[ai],x[aj],dx);
1339                         rsq11         = iprod(dx,dx);
1340
1341                         isai          = invsqrta[ai];
1342                         iq            = (-1)*facel*charge[ai];
1343
1344             rinv11        = gmx_invsqrt(rsq11);
1345             isaj          = invsqrta[aj];
1346             isaprod       = isai*isaj;
1347             qq            = isaprod*iq*charge[aj];
1348             gbscale       = isaprod*gbtabscale;
1349             r             = rsq11*rinv11;
1350             rt            = r*gbscale;
1351             n0            = rt;
1352             eps           = rt-n0;
1353             eps2          = eps*eps;
1354             nnn           = 4*n0;
1355             Y             = GBtab[nnn];
1356             F             = GBtab[nnn+1];
1357             Geps          = eps*GBtab[nnn+2];
1358             Heps2         = eps2*GBtab[nnn+3];
1359             Fp            = F+Geps+Heps2;
1360             VV            = Y+eps*Fp;
1361             FF            = Fp+Geps+2.0*Heps2;
1362             vgb           = qq*VV;
1363             fijC          = qq*FF*gbscale;
1364             dvdatmp       = -(vgb+fijC*r)*0.5;
1365             dvda[aj]      = dvda[aj] + dvdatmp*isaj*isaj;
1366             dvda[ai]      = dvda[ai] + dvdatmp*isai*isai;
1367             vctot         = vctot + vgb;
1368             fgb           = -(fijC)*rinv11;
1369
1370                         if (graph) {
1371                                 ivec_sub(SHIFT_IVEC(graph,ai),SHIFT_IVEC(graph,aj),dt);
1372                                 ki=IVEC2IS(dt);
1373                         }
1374
1375                         for (m=0; (m<DIM); m++) {                       /*  15          */
1376                                 fscal=fgb*dx[m];
1377                                 f[ai][m]+=fscal;
1378                                 f[aj][m]-=fscal;
1379                                 fshift[ki][m]+=fscal;
1380                                 fshift[CENTRAL][m]-=fscal;
1381                         }
1382         }
1383     }
1384
1385     return vctot;
1386 }
1387
1388 real calc_gb_selfcorrections(t_commrec *cr, int natoms,
1389                  real *charge, gmx_genborn_t *born, real *dvda, t_mdatoms *md, double facel)
1390 {
1391     int i,ai,at0,at1;
1392     real rai,e,derb,q,q2,fi,rai_inv,vtot;
1393
1394     if(PARTDECOMP(cr))
1395     {
1396         pd_at_range(cr,&at0,&at1);
1397     }
1398     else if(DOMAINDECOMP(cr))
1399     {
1400         at0=0;
1401         at1=cr->dd->nat_home;
1402     }
1403     else
1404     {
1405         at0=0;
1406         at1=natoms;
1407
1408     }
1409
1410     /* Scale the electrostatics by gb_epsilon_solvent */
1411     facel = facel * ((1.0/born->epsilon_r) - 1.0/born->gb_epsilon_solvent);
1412
1413     vtot=0.0;
1414
1415     /* Apply self corrections */
1416     for(i=at0;i<at1;i++)
1417     {
1418         ai       = i;
1419
1420         if(born->use[ai]==1)
1421         {
1422             rai      = born->bRad[ai];
1423             rai_inv  = 1.0/rai;
1424             q        = charge[ai];
1425             q2       = q*q;
1426             fi       = facel*q2;
1427             e        = fi*rai_inv;
1428             derb     = 0.5*e*rai_inv*rai_inv;
1429             dvda[ai] += derb*rai;
1430             vtot     -= 0.5*e;
1431         }
1432     }
1433
1434    return vtot;
1435
1436 }
1437
1438 real calc_gb_nonpolar(t_commrec *cr, t_forcerec *fr,int natoms,gmx_genborn_t *born, gmx_localtop_t *top,
1439                       const t_atomtypes *atype, real *dvda,int gb_algorithm, t_mdatoms *md)
1440 {
1441     int ai,i,at0,at1;
1442     real e,es,rai,rbi,term,probe,tmp,factor;
1443     real rbi_inv,rbi_inv2;
1444
1445     /* To keep the compiler happy */
1446     factor=0;
1447
1448     if(PARTDECOMP(cr))
1449     {
1450         pd_at_range(cr,&at0,&at1);
1451     }
1452     else if(DOMAINDECOMP(cr))
1453     {
1454         at0 = 0;
1455         at1 = cr->dd->nat_home;
1456     }
1457     else
1458     {
1459         at0=0;
1460         at1=natoms;
1461     }
1462
1463   /* factor is the surface tension */
1464   factor = born->sa_surface_tension;
1465   /*
1466
1467     // The surface tension factor is 0.0049 for Still model, 0.0054 for HCT/OBC
1468     if(gb_algorithm==egbSTILL)
1469     {
1470         factor=0.0049*100*CAL2JOULE;
1471     }
1472     else
1473     {
1474         factor=0.0054*100*CAL2JOULE;
1475     }
1476     */
1477     /* if(gb_algorithm==egbHCT || gb_algorithm==egbOBC) */
1478
1479     es    = 0;
1480     probe = 0.14;
1481     term  = M_PI*4;
1482
1483     for(i=at0;i<at1;i++)
1484     {
1485         ai        = i;
1486
1487         if(born->use[ai]==1)
1488         {
1489             rai       = top->atomtypes.gb_radius[md->typeA[ai]];
1490             rbi_inv   = fr->invsqrta[ai];
1491             rbi_inv2  = rbi_inv * rbi_inv;
1492             tmp       = (rai*rbi_inv2)*(rai*rbi_inv2);
1493             tmp       = tmp*tmp*tmp;
1494             e         = factor*term*(rai+probe)*(rai+probe)*tmp;
1495             dvda[ai]  = dvda[ai] - 6*e*rbi_inv2;
1496             es        = es + e;
1497         }
1498     }
1499
1500     return es;
1501 }
1502
1503
1504
1505 real calc_gb_chainrule(int natoms, t_nblist *nl, real *dadx, real *dvda, rvec x[], rvec t[], rvec fshift[],
1506                        rvec shift_vec[], int gb_algorithm, gmx_genborn_t *born, t_mdatoms *md)
1507 {
1508     int i,k,n,ai,aj,nj0,nj1,n0,n1;
1509     int shift;
1510     real shX,shY,shZ;
1511     real fgb,fij,rb2,rbi,fix1,fiy1,fiz1;
1512     real ix1,iy1,iz1,jx1,jy1,jz1,dx11,dy11,dz11,rsq11;
1513     real rinv11,tx,ty,tz,rbai,rbaj,fgb_ai;
1514     real *rb;
1515     volatile int idx;
1516
1517     n  = 0;
1518     rb = born->work;
1519
1520   n0 = 0;
1521   n1 = natoms;
1522
1523     if(gb_algorithm==egbSTILL)
1524     {
1525         for(i=n0;i<n1;i++)
1526         {
1527           rbi   = born->bRad[i];
1528           rb[i] = (2 * rbi * rbi * dvda[i])/ONE_4PI_EPS0;
1529         }
1530     }
1531     else if(gb_algorithm==egbHCT)
1532     {
1533         for(i=n0;i<n1;i++)
1534         {
1535           rbi   = born->bRad[i];
1536           rb[i] = rbi * rbi * dvda[i];
1537         }
1538     }
1539     else if(gb_algorithm==egbOBC)
1540     {
1541         for(i=n0;i<n1;i++)
1542         {
1543           rbi   = born->bRad[i];
1544           rb[i] = rbi * rbi * born->drobc[i] * dvda[i];
1545         }
1546     }
1547
1548     for(i=0;i<nl->nri;i++)
1549     {
1550         ai   = nl->iinr[i];
1551
1552         nj0  = nl->jindex[i];
1553         nj1  = nl->jindex[i+1];
1554
1555         /* Load shifts for this list */
1556         shift   = nl->shift[i];
1557         shX     = shift_vec[shift][0];
1558         shY     = shift_vec[shift][1];
1559         shZ     = shift_vec[shift][2];
1560
1561         /* Load atom i coordinates, add shift vectors */
1562         ix1  = shX + x[ai][0];
1563         iy1  = shY + x[ai][1];
1564         iz1  = shZ + x[ai][2];
1565
1566         fix1 = 0;
1567         fiy1 = 0;
1568         fiz1 = 0;
1569
1570         rbai = rb[ai];
1571
1572         for(k=nj0;k<nj1;k++)
1573         {
1574             aj = nl->jjnr[k];
1575
1576             jx1     = x[aj][0];
1577             jy1     = x[aj][1];
1578             jz1     = x[aj][2];
1579
1580             dx11    = ix1 - jx1;
1581             dy11    = iy1 - jy1;
1582             dz11    = iz1 - jz1;
1583
1584             rbaj    = rb[aj];
1585
1586             fgb     = rbai*dadx[n++];
1587             fgb_ai  = rbaj*dadx[n++];
1588
1589             /* Total force between ai and aj is the sum of ai->aj and aj->ai */
1590             fgb     = fgb + fgb_ai;
1591
1592             tx      = fgb * dx11;
1593             ty      = fgb * dy11;
1594             tz      = fgb * dz11;
1595
1596             fix1    = fix1 + tx;
1597             fiy1    = fiy1 + ty;
1598             fiz1    = fiz1 + tz;
1599
1600             /* Update force on atom aj */
1601             t[aj][0] = t[aj][0] - tx;
1602             t[aj][1] = t[aj][1] - ty;
1603             t[aj][2] = t[aj][2] - tz;
1604         }
1605
1606         /* Update force and shift forces on atom ai */
1607         t[ai][0] = t[ai][0] + fix1;
1608         t[ai][1] = t[ai][1] + fiy1;
1609         t[ai][2] = t[ai][2] + fiz1;
1610
1611         fshift[shift][0] = fshift[shift][0] + fix1;
1612         fshift[shift][1] = fshift[shift][1] + fiy1;
1613         fshift[shift][2] = fshift[shift][2] + fiz1;
1614
1615     }
1616
1617     return 0;
1618 }
1619
1620
1621 void
1622 calc_gb_forces(t_commrec *cr, t_mdatoms *md, gmx_genborn_t *born, gmx_localtop_t *top, const t_atomtypes *atype,
1623                rvec x[], rvec f[], t_forcerec *fr, t_idef *idef, int gb_algorithm, int sa_algorithm, t_nrnb *nrnb, gmx_bool bRad,
1624                const t_pbc *pbc, const t_graph *graph, gmx_enerdata_t *enerd)
1625 {
1626     real v=0;
1627     int  cnt;
1628     int i;
1629
1630         /* PBC or not? */
1631         const t_pbc *pbc_null;
1632
1633         if (fr->bMolPBC)
1634                 pbc_null = pbc;
1635         else
1636                 pbc_null = NULL;
1637
1638   if(sa_algorithm == esaAPPROX)
1639   {
1640     /* Do a simple ACE type approximation for the non-polar solvation */
1641     enerd->term[F_NPSOLVATION] += calc_gb_nonpolar(cr, fr,born->nr, born, top, atype, fr->dvda, gb_algorithm,md);
1642   }
1643
1644   /* Calculate the bonded GB-interactions using either table or analytical formula */
1645     enerd->term[F_GBPOL]       += gb_bonds_tab(x,f,fr->fshift, md->chargeA,&(fr->gbtabscale),
1646                                      fr->invsqrta,fr->dvda,fr->gbtab.tab,idef,born->epsilon_r,born->gb_epsilon_solvent, fr->epsfac, pbc_null, graph);
1647
1648     /* Calculate self corrections to the GB energies - currently only A state used! (FIXME) */
1649     enerd->term[F_GBPOL]       += calc_gb_selfcorrections(cr,born->nr,md->chargeA, born, fr->dvda, md, fr->epsfac);
1650
1651     /* If parallel, sum the derivative of the potential w.r.t the born radii */
1652     if(PARTDECOMP(cr))
1653     {
1654         gmx_sum(md->nr,fr->dvda, cr);
1655     }
1656     else if(DOMAINDECOMP(cr))
1657     {
1658         dd_atom_sum_real(cr->dd,fr->dvda);
1659         dd_atom_spread_real(cr->dd,fr->dvda);
1660     }
1661
1662     if(fr->bAllvsAll)
1663     {
1664 #if 0 && defined (GMX_X86_SSE2)
1665         if(fr->use_acceleration)
1666         {
1667 #  ifdef GMX_DOUBLE
1668             genborn_allvsall_calc_chainrule_sse2_double(fr,md,born,x[0],f[0],gb_algorithm,fr->AllvsAll_workgb);
1669 #  else
1670             genborn_allvsall_calc_chainrule_sse2_single(fr,md,born,x[0],f[0],gb_algorithm,fr->AllvsAll_workgb);
1671 #  endif
1672         }
1673         else
1674         {
1675             genborn_allvsall_calc_chainrule(fr,md,born,x[0],f[0],gb_algorithm,fr->AllvsAll_workgb);
1676         }
1677 #else
1678         genborn_allvsall_calc_chainrule(fr,md,born,x[0],f[0],gb_algorithm,fr->AllvsAll_workgb);
1679 #endif
1680         cnt = md->homenr*(md->nr/2+1);
1681         inc_nrnb(nrnb,eNR_BORN_AVA_CHAINRULE,cnt);
1682         inc_nrnb(nrnb,eNR_NBKERNEL_OUTER,md->homenr);
1683         return;
1684     }
1685
1686 #if 0 && defined (GMX_X86_SSE2)
1687     if(fr->use_acceleration)
1688     {
1689 #  ifdef GMX_DOUBLE
1690         calc_gb_chainrule_sse2_double(fr->natoms_force, &(fr->gblist),fr->dadx,fr->dvda,x[0],
1691                                       f[0],fr->fshift[0],fr->shift_vec[0],gb_algorithm,born,md);
1692 #  else
1693         calc_gb_chainrule_sse2_single(fr->natoms_force, &(fr->gblist),fr->dadx,fr->dvda,x[0],
1694                                       f[0],fr->fshift[0],fr->shift_vec[0],gb_algorithm,born,md);
1695 #  endif
1696     }
1697     else
1698     {
1699         calc_gb_chainrule(fr->natoms_force, &(fr->gblist), fr->dadx, fr->dvda,
1700                           x, f, fr->fshift, fr->shift_vec, gb_algorithm, born, md);
1701     }
1702 #else
1703     calc_gb_chainrule(fr->natoms_force, &(fr->gblist), fr->dadx, fr->dvda,
1704                       x, f, fr->fshift, fr->shift_vec, gb_algorithm, born, md);
1705 #endif
1706
1707     if(!fr->bAllvsAll)
1708     {
1709         inc_nrnb(nrnb,eNR_BORN_CHAINRULE,fr->gblist.nrj);
1710         inc_nrnb(nrnb,eNR_NBKERNEL_OUTER,fr->gblist.nri);
1711
1712     }
1713 }
1714
1715 static void add_j_to_gblist(gbtmpnbl_t *list,int aj)
1716 {
1717     if (list->naj >= list->aj_nalloc)
1718     {
1719         list->aj_nalloc = over_alloc_large(list->naj+1);
1720         srenew(list->aj,list->aj_nalloc);
1721     }
1722
1723     list->aj[list->naj++] = aj;
1724 }
1725
1726 static gbtmpnbl_t *find_gbtmplist(struct gbtmpnbls *lists,int shift)
1727 {
1728     int ind,i;
1729
1730     /* Search the list with the same shift, if there is one */
1731     ind = 0;
1732     while (ind < lists->nlist && shift != lists->list[ind].shift)
1733     {
1734         ind++;
1735     }
1736     if (ind == lists->nlist)
1737     {
1738         if (lists->nlist == lists->list_nalloc)
1739         {
1740             lists->list_nalloc++;
1741             srenew(lists->list,lists->list_nalloc);
1742             for(i=lists->nlist; i<lists->list_nalloc; i++)
1743             {
1744                 lists->list[i].aj        = NULL;
1745                 lists->list[i].aj_nalloc = 0;
1746             }
1747
1748         }
1749
1750         lists->list[lists->nlist].shift = shift;
1751         lists->list[lists->nlist].naj   = 0;
1752         lists->nlist++;
1753     }
1754
1755     return &lists->list[ind];
1756 }
1757
1758 static void add_bondeds_to_gblist(t_ilist *il,
1759                                   gmx_bool bMolPBC,t_pbc *pbc,t_graph *g,rvec *x,
1760                                   struct gbtmpnbls *nls)
1761 {
1762     int  ind,j,ai,aj,shift,found;
1763     rvec dx;
1764     ivec dt;
1765     gbtmpnbl_t *list;
1766
1767     shift = CENTRAL;
1768     for(ind=0; ind<il->nr; ind+=3)
1769     {
1770         ai = il->iatoms[ind+1];
1771         aj = il->iatoms[ind+2];
1772
1773         shift = CENTRAL;
1774         if (g != NULL)
1775         {
1776           rvec_sub(x[ai],x[aj],dx);
1777           ivec_sub(SHIFT_IVEC(g,ai),SHIFT_IVEC(g,aj),dt);
1778           shift = IVEC2IS(dt);
1779         }
1780         else if (bMolPBC)
1781         {
1782           shift = pbc_dx_aiuc(pbc,x[ai],x[aj],dx);
1783         }
1784
1785         /* Find the list for this shift or create one */
1786         list = find_gbtmplist(&nls[ai],shift);
1787
1788         found=0;
1789
1790         /* So that we do not add the same bond twice.
1791          * This happens with some constraints between 1-3 atoms
1792          * that are in the bond-list but should not be in the GB nb-list */
1793         for(j=0;j<list->naj;j++)
1794         {
1795             if (list->aj[j] == aj)
1796             {
1797                 found = 1;
1798             }
1799         }
1800
1801         if (found == 0)
1802         {
1803                         if(ai == aj)
1804                         {
1805                                 gmx_incons("ai == aj");
1806                         }
1807
1808             add_j_to_gblist(list,aj);
1809         }
1810     }
1811 }
1812
1813 static int
1814 compare_int (const void * a, const void * b)
1815 {
1816     return ( *(int*)a - *(int*)b );
1817 }
1818
1819
1820
1821 int make_gb_nblist(t_commrec *cr, int gb_algorithm, real gbcut,
1822                    rvec x[], matrix box,
1823                    t_forcerec *fr, t_idef *idef, t_graph *graph, gmx_genborn_t *born)
1824 {
1825     int i,l,ii,j,k,n,nj0,nj1,ai,aj,at0,at1,found,shift,s;
1826     int apa;
1827     t_nblist *nblist;
1828     t_pbc pbc;
1829
1830     struct gbtmpnbls *nls;
1831     gbtmpnbl_t *list =NULL;
1832
1833     set_pbc(&pbc,fr->ePBC,box);
1834     nls   = born->nblist_work;
1835
1836     for(i=0;i<born->nr;i++)
1837     {
1838         nls[i].nlist = 0;
1839     }
1840
1841     if (fr->bMolPBC)
1842     {
1843         set_pbc_dd(&pbc,fr->ePBC,cr->dd,TRUE,box);
1844     }
1845
1846     switch (gb_algorithm)
1847     {
1848     case egbHCT:
1849     case egbOBC:
1850         /* Loop over 1-2, 1-3 and 1-4 interactions */
1851         for(j=F_GB12;j<=F_GB14;j++)
1852         {
1853             add_bondeds_to_gblist(&idef->il[j],fr->bMolPBC,&pbc,graph,x,nls);
1854         }
1855         break;
1856     case egbSTILL:
1857         /* Loop over 1-4 interactions */
1858         add_bondeds_to_gblist(&idef->il[F_GB14],fr->bMolPBC,&pbc,graph,x,nls);
1859         break;
1860     default:
1861         gmx_incons("Unknown GB algorithm");
1862     }
1863
1864     /* Loop over the VDWQQ and VDW nblists to set up the nonbonded part of the GB list */
1865     for(n=0; (n<fr->nnblists); n++)
1866     {
1867         for(i=0; (i<eNL_NR); i++)
1868         {
1869             nblist=&(fr->nblists[n].nlist_sr[i]);
1870
1871             if (nblist->nri > 0 && (i==eNL_VDWQQ || i==eNL_QQ))
1872             {
1873                 for(j=0;j<nblist->nri;j++)
1874                 {
1875                     ai    = nblist->iinr[j];
1876                     shift = nblist->shift[j];
1877
1878                     /* Find the list for this shift or create one */
1879                     list = find_gbtmplist(&nls[ai],shift);
1880
1881                     nj0 = nblist->jindex[j];
1882                     nj1 = nblist->jindex[j+1];
1883
1884                     /* Add all the j-atoms in the non-bonded list to the GB list */
1885                     for(k=nj0;k<nj1;k++)
1886                     {
1887                         add_j_to_gblist(list,nblist->jjnr[k]);
1888                     }
1889                 }
1890             }
1891         }
1892     }
1893
1894     /* Zero out some counters */
1895         fr->gblist.nri=0;
1896     fr->gblist.nrj=0;
1897
1898         fr->gblist.jindex[0] = fr->gblist.nri;
1899
1900         for(i=0;i<fr->natoms_force;i++)
1901     {
1902         for(s=0; s<nls[i].nlist; s++)
1903         {
1904             list = &nls[i].list[s];
1905
1906             /* Only add those atoms that actually have neighbours */
1907             if (born->use[i] != 0)
1908             {
1909                 fr->gblist.iinr[fr->gblist.nri]  = i;
1910                 fr->gblist.shift[fr->gblist.nri] = list->shift;
1911                 fr->gblist.nri++;
1912
1913                 for(k=0; k<list->naj; k++)
1914                 {
1915                     /* Memory allocation for jjnr */
1916                     if(fr->gblist.nrj >= fr->gblist.maxnrj)
1917                     {
1918                         fr->gblist.maxnrj += over_alloc_large(fr->gblist.maxnrj);
1919
1920                         if (debug)
1921                         {
1922                             fprintf(debug,"Increasing GB neighbourlist j size to %d\n",fr->gblist.maxnrj);
1923                         }
1924
1925                         srenew(fr->gblist.jjnr,fr->gblist.maxnrj);
1926                     }
1927
1928                     /* Put in list */
1929                                         if(i == list->aj[k])
1930                                         {
1931                                                 gmx_incons("i == list->aj[k]");
1932                                         }
1933                     fr->gblist.jjnr[fr->gblist.nrj++] = list->aj[k];
1934                 }
1935
1936                                 fr->gblist.jindex[fr->gblist.nri] = fr->gblist.nrj;
1937             }
1938                 }
1939         }
1940
1941
1942 #ifdef SORT_GB_LIST
1943     for(i=0;i<fr->gblist.nri;i++)
1944     {
1945         nj0 = fr->gblist.jindex[i];
1946         nj1 = fr->gblist.jindex[i+1];
1947         ai  = fr->gblist.iinr[i];
1948
1949         /* Temporary fix */
1950                 for(j=nj0;j<nj1;j++)
1951                 {
1952             if(fr->gblist.jjnr[j]<ai)
1953                 fr->gblist.jjnr[j]+=fr->natoms_force;
1954         }
1955         qsort(fr->gblist.jjnr+nj0,nj1-nj0,sizeof(int),compare_int);
1956         /* Fix back */
1957         for(j=nj0;j<nj1;j++)
1958         {
1959             if(fr->gblist.jjnr[j]>=fr->natoms_force)
1960                 fr->gblist.jjnr[j]-=fr->natoms_force;
1961         }
1962
1963     }
1964 #endif
1965
1966     return 0;
1967 }
1968
1969 void make_local_gb(const t_commrec *cr, gmx_genborn_t *born, int gb_algorithm)
1970 {
1971     int i,at0,at1;
1972     gmx_domdec_t *dd=NULL;
1973
1974     if(DOMAINDECOMP(cr))
1975     {
1976         dd = cr->dd;
1977         at0 = 0;
1978         at1 = dd->nat_tot;
1979     }
1980     else
1981     {
1982         /* Single node or particle decomp (global==local), just copy pointers and return */
1983         if(gb_algorithm==egbSTILL)
1984         {
1985             born->gpol      = born->gpol_globalindex;
1986             born->vsolv     = born->vsolv_globalindex;
1987             born->gb_radius = born->gb_radius_globalindex;
1988         }
1989         else
1990         {
1991             born->param     = born->param_globalindex;
1992             born->gb_radius = born->gb_radius_globalindex;
1993         }
1994
1995         born->use = born->use_globalindex;
1996
1997         return;
1998     }
1999
2000     /* Reallocation of local arrays if necessary */
2001     /* fr->natoms_force is equal to dd->nat_tot */
2002     if (DOMAINDECOMP(cr) && dd->nat_tot > born->nalloc)
2003     {
2004         int nalloc;
2005
2006         nalloc = dd->nat_tot;
2007
2008         /* Arrays specific to different gb algorithms */
2009         if (gb_algorithm == egbSTILL)
2010         {
2011             srenew(born->gpol,  nalloc+3);
2012             srenew(born->vsolv, nalloc+3);
2013             srenew(born->gb_radius, nalloc+3);
2014             for(i=born->nalloc; (i<nalloc+3); i++)
2015             {
2016                 born->gpol[i] = 0;
2017                 born->vsolv[i] = 0;
2018                 born->gb_radius[i] = 0;
2019             }
2020         }
2021         else
2022         {
2023             srenew(born->param, nalloc+3);
2024             srenew(born->gb_radius, nalloc+3);
2025             for(i=born->nalloc; (i<nalloc+3); i++)
2026             {
2027                 born->param[i] = 0;
2028                 born->gb_radius[i] = 0;
2029             }
2030         }
2031
2032         /* All gb-algorithms use the array for vsites exclusions */
2033         srenew(born->use,    nalloc+3);
2034         for(i=born->nalloc; (i<nalloc+3); i++)
2035         {
2036             born->use[i] = 0;
2037         }
2038
2039         born->nalloc = nalloc;
2040     }
2041
2042     /* With dd, copy algorithm specific arrays */
2043     if(gb_algorithm==egbSTILL)
2044     {
2045         for(i=at0;i<at1;i++)
2046         {
2047             born->gpol[i]  = born->gpol_globalindex[dd->gatindex[i]];
2048             born->vsolv[i] = born->vsolv_globalindex[dd->gatindex[i]];
2049             born->gb_radius[i] = born->gb_radius_globalindex[dd->gatindex[i]];
2050             born->use[i]   = born->use_globalindex[dd->gatindex[i]];
2051         }
2052     }
2053     else
2054     {
2055         for(i=at0;i<at1;i++)
2056         {
2057             born->param[i]     = born->param_globalindex[dd->gatindex[i]];
2058             born->gb_radius[i] = born->gb_radius_globalindex[dd->gatindex[i]];
2059             born->use[i]       = born->use_globalindex[dd->gatindex[i]];
2060         }
2061     }
2062 }
2063