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 && nl->jjnr[k]>=0;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             born->gpol_still_work[aj] += prod_ai*icf4;
 574             gpi             = gpi+prod*icf4;
 575
 576             /* Save ai->aj and aj->ai chain rule terms */
 577             fr->dadx[n++]   = prod*icf6;
 578             fr->dadx[n++]   = prod_ai*icf6;
 579         }
 580         born->gpol_still_work[ai]+=gpi;
 581     }
 582
 583     /* Parallel summations */
 584     if(PARTDECOMP(cr))
 585     {
 586         gmx_sum(natoms, born->gpol_still_work, cr);
 587     }
 588     else if(DOMAINDECOMP(cr))
 589     {
 590         dd_atom_sum_real(cr->dd, born->gpol_still_work);
 591     }
 592
 593     /* Calculate the radii */
 594         for(i=0;i<fr->natoms_force;i++) /* PELA born->nr */
 595     {
 596                 if(born->use[i] != 0)
 597         {
 598             gpi  = born->gpol[i]+born->gpol_still_work[i];
 599             gpi2 = gpi * gpi;
 600             born->bRad[i]   = factor*gmx_invsqrt(gpi2);
 601             fr->invsqrta[i] = gmx_invsqrt(born->bRad[i]);
 602                 }
 603     }
 604
 605     /* Extra communication required for DD */
 606     if(DOMAINDECOMP(cr))
 607     {
 608         dd_atom_spread_real(cr->dd, born->bRad);
 609         dd_atom_spread_real(cr->dd, fr->invsqrta);
 610     }
 611
 612     return 0;
 613
 614 }
 615
 616
 617 static int
 618 calc_gb_rad_hct(t_commrec *cr,t_forcerec *fr,int natoms, gmx_localtop_t *top,
 619                 const t_atomtypes *atype, rvec x[], t_nblist *nl,
 620                 gmx_genborn_t *born,t_mdatoms *md)
 621 {
 622     int i,k,n,ai,aj,nj0,nj1,at0,at1;
 623     int shift;
 624     real shX,shY,shZ;
 625     real rai,raj,gpi,dr2,dr,sk,sk_ai,sk2,sk2_ai,lij,uij,diff2,tmp,sum_ai;
 626     real rad,min_rad,rinv,rai_inv;
 627     real ix1,iy1,iz1,jx1,jy1,jz1,dx11,dy11,dz11;
 628     real lij2, uij2, lij3, uij3, t1,t2,t3;
 629     real lij_inv,dlij,duij,sk2_rinv,prod,log_term;
 630     real doffset,raj_inv,dadx_val;
 631     real *gb_radius;
 632
 633     doffset = born->gb_doffset;
 634     gb_radius = born->gb_radius;
 635
 636     for(i=0;i<born->nr;i++)
 637     {
 638         born->gpol_hct_work[i] = 0;
 639     }
 640
 641     /* Keep the compiler happy */
 642     n    = 0;
 643     prod = 0;
 644
 645     for(i=0;i<nl->nri;i++)
 646     {
 647         ai     = nl->iinr[i];
 648
 649         nj0    = nl->jindex[i];
 650         nj1    = nl->jindex[i+1];
 651
 652         /* Load shifts for this list */
 653         shift   = nl->shift[i];
 654         shX     = fr->shift_vec[shift][0];
 655         shY     = fr->shift_vec[shift][1];
 656         shZ     = fr->shift_vec[shift][2];
 657
 658         rai     = gb_radius[ai];
 659         rai_inv = 1.0/rai;
 660
 661         sk_ai   = born->param[ai];
 662         sk2_ai  = sk_ai*sk_ai;
 663
 664         /* Load atom i coordinates, add shift vectors */
 665         ix1     = shX + x[ai][0];
 666         iy1     = shY + x[ai][1];
 667         iz1     = shZ + x[ai][2];
 668
 669         sum_ai  = 0;
 670
 671         for(k=nj0;k<nj1 && nl->jjnr[k]>=0;k++)
 672         {
 673             aj    = nl->jjnr[k];
 674
 675             jx1   = x[aj][0];
 676             jy1   = x[aj][1];
 677             jz1   = x[aj][2];
 678
 679             dx11  = ix1 - jx1;
 680             dy11  = iy1 - jy1;
 681             dz11  = iz1 - jz1;
 682
 683             dr2   = dx11*dx11+dy11*dy11+dz11*dz11;
 684             rinv  = gmx_invsqrt(dr2);
 685             dr    = rinv*dr2;
 686
 687             sk    = born->param[aj];
 688             raj   = gb_radius[aj];
 689
 690             /* aj -> ai interaction */
 691             if(rai < dr+sk)
 692             {
 693                 lij     = 1.0/(dr-sk);
 694                 dlij    = 1.0;
 695
 696                 if(rai>dr-sk)
 697                 {
 698                     lij  = rai_inv;
 699                     dlij = 0.0;
 700                 }
 701
 702                 lij2     = lij*lij;
 703                 lij3     = lij2*lij;
 704
 705                 uij      = 1.0/(dr+sk);
 706                 uij2     = uij*uij;
 707                 uij3     = uij2*uij;
 708
 709                 diff2    = uij2-lij2;
 710
 711                 lij_inv  = gmx_invsqrt(lij2);
 712                 sk2      = sk*sk;
 713                 sk2_rinv = sk2*rinv;
 714                 prod     = 0.25*sk2_rinv;
 715
 716                 log_term = log(uij*lij_inv);
 717
 718                 tmp      = lij-uij + 0.25*dr*diff2 + (0.5*rinv)*log_term +
 719                     prod*(-diff2);
 720
 721                 if(rai<sk-dr)
 722                 {
 723                     tmp = tmp + 2.0 * (rai_inv-lij);
 724                 }
 725
 726                 t1 = 0.5*lij2 + prod*lij3 - 0.25*(lij*rinv+lij3*dr);
 727                 t2 = -0.5*uij2 - 0.25*sk2_rinv*uij3 + 0.25*(uij*rinv+uij3*dr);
 728                 t3 = 0.125*(1.0+sk2_rinv*rinv)*(-diff2)+0.25*log_term*rinv*rinv;
 729
 730                 dadx_val = (dlij*t1+t2+t3)*rinv; /* rb2 is moved to chainrule */
 731                 /* fr->dadx[n++] = (dlij*t1+duij*t2+t3)*rinv; */
 732                 /* rb2 is moved to chainrule    */
 733
 734                 sum_ai += 0.5*tmp;
 735             }
 736             else
 737             {
 738                 dadx_val = 0.0;
 739             }
 740             fr->dadx[n++] = dadx_val;
 741
 742
 743             /* ai -> aj interaction */
 744             if(raj < dr + sk_ai)
 745             {
 746                 lij     = 1.0/(dr-sk_ai);
 747                 dlij    = 1.0;
 748                 raj_inv = 1.0/raj;
 749
 750                 if(raj>dr-sk_ai)
 751                 {
 752                     lij = raj_inv;
 753                     dlij = 0.0;
 754                 }
 755
 756                 lij2     = lij  * lij;
 757                 lij3     = lij2 * lij;
 758
 759                 uij      = 1.0/(dr+sk_ai);
 760                 uij2     = uij  * uij;
 761                 uij3     = uij2 * uij;
 762
 763                 diff2    = uij2-lij2;
 764
 765                 lij_inv  = gmx_invsqrt(lij2);
 766                 sk2      =  sk2_ai; /* sk2_ai = sk_ai * sk_ai in i loop above */
 767                 sk2_rinv = sk2*rinv;
 768                 prod     = 0.25 * sk2_rinv;
 769
 770                 /* log_term = table_log(uij*lij_inv,born->log_table,
 771                    LOG_TABLE_ACCURACY); */
 772                 log_term = log(uij*lij_inv);
 773
 774                 tmp      = lij-uij + 0.25*dr*diff2 + (0.5*rinv)*log_term +
 775                            prod*(-diff2);
 776
 777                 if(raj<sk_ai-dr)
 778                 {
 779                     tmp     = tmp + 2.0 * (raj_inv-lij);
 780                 }
 781
 782                 /* duij = 1.0 */
 783                 t1      = 0.5*lij2 + prod*lij3 - 0.25*(lij*rinv+lij3*dr);
 784                 t2      = -0.5*uij2 - 0.25*sk2_rinv*uij3 + 0.25*(uij*rinv+uij3*dr);
 785                 t3      = 0.125*(1.0+sk2_rinv*rinv)*(-diff2)+0.25*log_term*rinv*rinv;
 786
 787                 dadx_val = (dlij*t1+t2+t3)*rinv; /* rb2 is moved to chainrule    */
 788                 /* fr->dadx[n++] = (dlij*t1+duij*t2+t3)*rinv; */ /* rb2 is moved to chainrule    */
 789
 790                 born->gpol_hct_work[aj] += 0.5*tmp;
 791             }
 792             else
 793             {
 794                 dadx_val = 0.0;
 795             }
 796             fr->dadx[n++] = dadx_val;
 797         }
 798
 799         born->gpol_hct_work[ai] += sum_ai;
 800     }
 801
 802     /* Parallel summations */
 803     if(PARTDECOMP(cr))
 804     {
 805         gmx_sum(natoms, born->gpol_hct_work, cr);
 806     }
 807     else if(DOMAINDECOMP(cr))
 808     {
 809         dd_atom_sum_real(cr->dd, born->gpol_hct_work);
 810     }
 811
 812     for(i=0;i<fr->natoms_force;i++) /* PELA born->nr */
 813     {
 814                 if(born->use[i] != 0)
 815         {
 816             rai     = top->atomtypes.gb_radius[md->typeA[i]]-doffset;
 817             sum_ai  = 1.0/rai - born->gpol_hct_work[i];
 818             min_rad = rai + doffset;
 819             rad     = 1.0/sum_ai;
 820
 821             born->bRad[i]   = rad > min_rad ? rad : min_rad;
 822             fr->invsqrta[i] = gmx_invsqrt(born->bRad[i]);
 823         }
 824     }
 825
 826     /* Extra communication required for DD */
 827     if(DOMAINDECOMP(cr))
 828     {
 829         dd_atom_spread_real(cr->dd, born->bRad);
 830         dd_atom_spread_real(cr->dd, fr->invsqrta);
 831     }
 832
 833
 834     return 0;
 835 }
 836
 837 static int
 838 calc_gb_rad_obc(t_commrec *cr, t_forcerec *fr, int natoms, gmx_localtop_t *top,
 839                     const t_atomtypes *atype, rvec x[], t_nblist *nl, gmx_genborn_t *born,t_mdatoms *md)
 840 {
 841     int i,k,ai,aj,nj0,nj1,n,at0,at1;
 842     int shift;
 843     real shX,shY,shZ;
 844     real rai,raj,gpi,dr2,dr,sk,sk2,lij,uij,diff2,tmp,sum_ai;
 845     real rad, min_rad,sum_ai2,sum_ai3,tsum,tchain,rinv,rai_inv,lij_inv,rai_inv2;
 846     real log_term,prod,sk2_rinv,sk_ai,sk2_ai;
 847     real ix1,iy1,iz1,jx1,jy1,jz1,dx11,dy11,dz11;
 848     real lij2,uij2,lij3,uij3,dlij,duij,t1,t2,t3;
 849     real doffset,raj_inv,dadx_val;
 850     real *gb_radius;
 851
 852     /* Keep the compiler happy */
 853     n    = 0;
 854     prod = 0;
 855     raj  = 0;
 856
 857     doffset = born->gb_doffset;
 858     gb_radius = born->gb_radius;
 859
 860     for(i=0;i<born->nr;i++)
 861     {
 862         born->gpol_hct_work[i] = 0;
 863     }
 864
 865     for(i=0;i<nl->nri;i++)
 866     {
 867         ai      = nl->iinr[i];
 868
 869         nj0     = nl->jindex[i];
 870         nj1     = nl->jindex[i+1];
 871
 872         /* Load shifts for this list */
 873         shift   = nl->shift[i];
 874         shX     = fr->shift_vec[shift][0];
 875         shY     = fr->shift_vec[shift][1];
 876         shZ     = fr->shift_vec[shift][2];
 877
 878         rai      = gb_radius[ai];
 879         rai_inv  = 1.0/rai;
 880
 881         sk_ai    = born->param[ai];
 882         sk2_ai   = sk_ai*sk_ai;
 883
 884         /* Load atom i coordinates, add shift vectors */
 885         ix1      = shX + x[ai][0];
 886         iy1      = shY + x[ai][1];
 887         iz1      = shZ + x[ai][2];
 888
 889         sum_ai   = 0;
 890
 891         for(k=nj0;k<nj1 && nl->jjnr[k]>=0;k++)
 892         {
 893             aj    = nl->jjnr[k];
 894
 895             jx1   = x[aj][0];
 896             jy1   = x[aj][1];
 897             jz1   = x[aj][2];
 898
 899             dx11  = ix1 - jx1;
 900             dy11  = iy1 - jy1;
 901             dz11  = iz1 - jz1;
 902
 903             dr2   = dx11*dx11+dy11*dy11+dz11*dz11;
 904             rinv  = gmx_invsqrt(dr2);
 905             dr    = dr2*rinv;
 906
 907             /* sk is precalculated in init_gb() */
 908             sk    = born->param[aj];
 909             raj   = gb_radius[aj];
 910
 911             /* aj -> ai interaction */
 912             if(rai < dr+sk)
 913             {
 914                 lij       = 1.0/(dr-sk);
 915                 dlij      = 1.0;
 916
 917                 if(rai>dr-sk)
 918                 {
 919                     lij  = rai_inv;
 920                     dlij = 0.0;
 921                 }
 922
 923                 uij      = 1.0/(dr+sk);
 924                 lij2     = lij  * lij;
 925                 lij3     = lij2 * lij;
 926                 uij2     = uij  * uij;
 927                 uij3     = uij2 * uij;
 928
 929                 diff2    = uij2-lij2;
 930
 931                 lij_inv  = gmx_invsqrt(lij2);
 932                 sk2      = sk*sk;
 933                 sk2_rinv = sk2*rinv;
 934                 prod     = 0.25*sk2_rinv;
 935
 936                 log_term = log(uij*lij_inv);
 937
 938                 tmp      = lij-uij + 0.25*dr*diff2 + (0.5*rinv)*log_term + prod*(-diff2);
 939
 940                 if(rai < sk-dr)
 941                 {
 942                     tmp = tmp + 2.0 * (rai_inv-lij);
 943                 }
 944
 945                 /* duij    = 1.0; */
 946                 t1      = 0.5*lij2 + prod*lij3 - 0.25*(lij*rinv+lij3*dr);
 947                 t2      = -0.5*uij2 - 0.25*sk2_rinv*uij3 + 0.25*(uij*rinv+uij3*dr);
 948                 t3      = 0.125*(1.0+sk2_rinv*rinv)*(-diff2)+0.25*log_term*rinv*rinv;
 949
 950                 dadx_val = (dlij*t1+t2+t3)*rinv; /* rb2 is moved to chainrule    */
 951
 952                 sum_ai += 0.5*tmp;
 953             }
 954             else
 955             {
 956                 dadx_val = 0.0;
 957             }
 958             fr->dadx[n++] = dadx_val;
 959
 960             /* ai -> aj interaction */
 961             if(raj < dr + sk_ai)
 962             {
 963                 lij     = 1.0/(dr-sk_ai);
 964                 dlij    = 1.0;
 965                 raj_inv = 1.0/raj;
 966
 967                 if(raj>dr-sk_ai)
 968                 {
 969                     lij = raj_inv;
 970                     dlij = 0.0;
 971                 }
 972
 973                 lij2     = lij  * lij;
 974                 lij3     = lij2 * lij;
 975
 976                 uij      = 1.0/(dr+sk_ai);
 977                 uij2     = uij  * uij;
 978                 uij3     = uij2 * uij;
 979
 980                 diff2    = uij2-lij2;
 981
 982                 lij_inv  = gmx_invsqrt(lij2);
 983                 sk2      =  sk2_ai; /* sk2_ai = sk_ai * sk_ai in i loop above */
 984                 sk2_rinv = sk2*rinv;
 985                 prod     = 0.25 * sk2_rinv;
 986
 987                 /* log_term = table_log(uij*lij_inv,born->log_table,LOG_TABLE_ACCURACY); */
 988                 log_term = log(uij*lij_inv);
 989
 990                 tmp      = lij-uij + 0.25*dr*diff2 + (0.5*rinv)*log_term + prod*(-diff2);
 991
 992                 if(raj<sk_ai-dr)
 993                 {
 994                     tmp     = tmp + 2.0 * (raj_inv-lij);
 995                 }
 996
 997                 t1      = 0.5*lij2 + prod*lij3 - 0.25*(lij*rinv+lij3*dr);
 998                 t2      = -0.5*uij2 - 0.25*sk2_rinv*uij3 + 0.25*(uij*rinv+uij3*dr);
 999                 t3      = 0.125*(1.0+sk2_rinv*rinv)*(-diff2)+0.25*log_term*rinv*rinv;
1000
1001                 dadx_val = (dlij*t1+t2+t3)*rinv; /* rb2 is moved to chainrule    */
1002
1003                 born->gpol_hct_work[aj] += 0.5*tmp;
1004
1005             }
1006             else
1007             {
1008                 dadx_val = 0.0;
1009             }
1010             fr->dadx[n++] = dadx_val;
1011
1012         }
1013         born->gpol_hct_work[ai] += sum_ai;
1014
1015     }
1016
1017     /* Parallel summations */
1018     if(PARTDECOMP(cr))
1019     {
1020         gmx_sum(natoms, born->gpol_hct_work, cr);
1021     }
1022     else if(DOMAINDECOMP(cr))
1023     {
1024         dd_atom_sum_real(cr->dd, born->gpol_hct_work);
1025     }
1026
1027     for(i=0;i<fr->natoms_force;i++) /* PELA born->nr */
1028     {
1029                 if(born->use[i] != 0)
1030         {
1031             rai        = top->atomtypes.gb_radius[md->typeA[i]];
1032             rai_inv2   = 1.0/rai;
1033             rai        = rai-doffset;
1034             rai_inv    = 1.0/rai;
1035             sum_ai     = rai * born->gpol_hct_work[i];
1036             sum_ai2    = sum_ai  * sum_ai;
1037             sum_ai3    = sum_ai2 * sum_ai;
1038
1039             tsum    = tanh(born->obc_alpha*sum_ai-born->obc_beta*sum_ai2+born->obc_gamma*sum_ai3);
1040             born->bRad[i] = rai_inv - tsum*rai_inv2;
1041             born->bRad[i] = 1.0 / born->bRad[i];
1042
1043             fr->invsqrta[i] = gmx_invsqrt(born->bRad[i]);
1044
1045             tchain  = rai * (born->obc_alpha-2*born->obc_beta*sum_ai+3*born->obc_gamma*sum_ai2);
1046             born->drobc[i] = (1.0-tsum*tsum)*tchain*rai_inv2;
1047         }
1048     }
1049
1050     /* Extra (local) communication required for DD */
1051     if(DOMAINDECOMP(cr))
1052     {
1053         dd_atom_spread_real(cr->dd, born->bRad);
1054         dd_atom_spread_real(cr->dd, fr->invsqrta);
1055         dd_atom_spread_real(cr->dd, born->drobc);
1056     }
1057
1058     return 0;
1059
1060 }
1061
1062
1063
1064 int calc_gb_rad(t_commrec *cr, t_forcerec *fr, t_inputrec *ir,gmx_localtop_t *top,
1065                 const t_atomtypes *atype, rvec x[], t_nblist *nl, gmx_genborn_t *born,t_mdatoms *md,t_nrnb     *nrnb)
1066 {
1067     real *p;
1068     int   cnt;
1069     int ndadx;
1070
1071     if(fr->bAllvsAll && fr->dadx==NULL)
1072     {
1073         /* We might need up to 8 atoms of padding before and after,
1074          * and another 4 units to guarantee SSE alignment.
1075          */
1076         fr->nalloc_dadx = 2*(md->homenr+12)*(md->nr/2+1+12);
1077         snew(fr->dadx_rawptr,fr->nalloc_dadx);
1078         fr->dadx = (real *) (((size_t) fr->dadx_rawptr + 16) & (~((size_t) 15)));
1079     }
1080     else
1081     {
1082         /* In the SSE-enabled gb-loops, when writing to dadx, we
1083          * always write 2*4 elements at a time, even in the case with only
1084          * 1-3 j particles, where we only really need to write 2*(1-3)
1085          * elements. This is because we want dadx to be aligned to a 16-
1086          * byte boundary, and being able to use _mm_store/load_ps
1087          */
1088         ndadx = 2 * (nl->nrj + 3*nl->nri);
1089
1090         /* First, reallocate the dadx array, we need 3 extra for SSE */
1091         if (ndadx + 3 > fr->nalloc_dadx)
1092         {
1093             fr->nalloc_dadx = over_alloc_large(ndadx) + 3;
1094             srenew(fr->dadx_rawptr,fr->nalloc_dadx);
1095             fr->dadx = (real *) (((size_t) fr->dadx_rawptr + 16) & (~((size_t) 15)));
1096         }
1097     }
1098
1099     if(fr->bAllvsAll)
1100     {
1101         cnt = md->homenr*(md->nr/2+1);
1102
1103         if(ir->gb_algorithm==egbSTILL)
1104         {
1105 #if 0 && defined (GMX_X86_SSE2)
1106             if(fr->use_acceleration)
1107             {
1108 #  ifdef GMX_DOUBLE
1109                 genborn_allvsall_calc_still_radii_sse2_double(fr,md,born,top,x[0],cr,&fr->AllvsAll_workgb);
1110 #  else
1111                 genborn_allvsall_calc_still_radii_sse2_single(fr,md,born,top,x[0],cr,&fr->AllvsAll_workgb);
1112 #  endif
1113             }
1114             else
1115             {
1116                 genborn_allvsall_calc_still_radii(fr,md,born,top,x[0],cr,&fr->AllvsAll_workgb);
1117             }
1118 #else
1119             genborn_allvsall_calc_still_radii(fr,md,born,top,x[0],cr,&fr->AllvsAll_workgb);
1120 #endif
1121             /* 13 flops in outer loop, 47 flops in inner loop */
1122             inc_nrnb(nrnb,eNR_BORN_AVA_RADII_STILL,md->homenr*13+cnt*47);
1123         }
1124         else if(ir->gb_algorithm==egbHCT || ir->gb_algorithm==egbOBC)
1125         {
1126 #if 0 && defined (GMX_X86_SSE2)
1127             if(fr->use_acceleration)
1128             {
1129 #  ifdef GMX_DOUBLE
1130                 genborn_allvsall_calc_hct_obc_radii_sse2_double(fr,md,born,ir->gb_algorithm,top,x[0],cr,&fr->AllvsAll_workgb);
1131 #  else
1132                 genborn_allvsall_calc_hct_obc_radii_sse2_single(fr,md,born,ir->gb_algorithm,top,x[0],cr,&fr->AllvsAll_workgb);
1133 #  endif
1134             }
1135             else
1136             {
1137                 genborn_allvsall_calc_hct_obc_radii(fr,md,born,ir->gb_algorithm,top,x[0],cr,&fr->AllvsAll_workgb);
1138             }
1139 #else
1140             genborn_allvsall_calc_hct_obc_radii(fr,md,born,ir->gb_algorithm,top,x[0],cr,&fr->AllvsAll_workgb);
1141 #endif
1142             /* 24 flops in outer loop, 183 in inner */
1143             inc_nrnb(nrnb,eNR_BORN_AVA_RADII_HCT_OBC,md->homenr*24+cnt*183);
1144         }
1145         else
1146         {
1147             gmx_fatal(FARGS,"Bad gb algorithm for all-vs-all interactions");
1148         }
1149         return 0;
1150     }
1151
1152     /* Switch for determining which algorithm to use for Born radii calculation */
1153 #ifdef GMX_DOUBLE
1154
1155 #if 0 && defined (GMX_X86_SSE2)
1156     /* x86 or x86-64 with GCC inline assembly and/or SSE intrinsics */
1157     switch(ir->gb_algorithm)
1158     {
1159         case egbSTILL:
1160             if(fr->use_acceleration)
1161             {
1162                 calc_gb_rad_still_sse2_double(cr,fr,born->nr,top, atype, x[0], nl, born);
1163             }
1164             else
1165             {
1166                 calc_gb_rad_still(cr,fr,born->nr,top,atype,x,nl,born,md);
1167             }
1168             break;
1169         case egbHCT:
1170             if(fr->use_acceleration)
1171             {
1172                 calc_gb_rad_hct_obc_sse2_double(cr,fr,born->nr,top, atype, x[0], nl, born, md, ir->gb_algorithm);
1173             }
1174             else
1175             {
1176                 calc_gb_rad_hct(cr,fr,born->nr,top,atype,x,nl,born,md);
1177             }
1178             break;
1179         case egbOBC:
1180             if(fr->use_acceleration)
1181             {
1182                 calc_gb_rad_hct_obc_sse2_double(cr,fr,born->nr,top, atype, x[0], nl, born, md, ir->gb_algorithm);
1183             }
1184             else
1185             {
1186                 calc_gb_rad_obc(cr,fr,born->nr,top,atype,x,nl,born,md);
1187             }
1188             break;
1189
1190         default:
1191             gmx_fatal(FARGS, "Unknown double precision sse-enabled algorithm for Born radii calculation: %d",ir->gb_algorithm);
1192     }
1193 #else
1194     switch(ir->gb_algorithm)
1195     {
1196         case egbSTILL:
1197             calc_gb_rad_still(cr,fr,born->nr,top,atype,x,nl,born,md);
1198             break;
1199         case egbHCT:
1200             calc_gb_rad_hct(cr,fr,born->nr,top,atype,x,nl,born,md);
1201             break;
1202         case egbOBC:
1203             calc_gb_rad_obc(cr,fr,born->nr,top,atype,x,nl,born,md);
1204             break;
1205
1206         default:
1207             gmx_fatal(FARGS, "Unknown double precision algorithm for Born radii calculation: %d",ir->gb_algorithm);
1208     }
1209
1210 #endif
1211
1212 #else
1213
1214 #if 0 && defined (GMX_X86_SSE2)
1215     /* x86 or x86-64 with GCC inline assembly and/or SSE intrinsics */
1216     switch(ir->gb_algorithm)
1217     {
1218         case egbSTILL:
1219             if(fr->use_acceleration)
1220             {
1221             calc_gb_rad_still_sse2_single(cr,fr,born->nr,top, atype, x[0], nl, born);
1222             }
1223             else
1224             {
1225                 calc_gb_rad_still(cr,fr,born->nr,top,atype,x,nl,born,md);
1226             }
1227             break;
1228         case egbHCT:
1229                 if(fr->use_acceleration)
1230                 {
1231                     calc_gb_rad_hct_obc_sse2_single(cr,fr,born->nr,top, atype, x[0], nl, born, md, ir->gb_algorithm);
1232                 }
1233                 else
1234                 {
1235                     calc_gb_rad_hct(cr,fr,born->nr,top,atype,x,nl,born,md);
1236                 }
1237             break;
1238
1239         case egbOBC:
1240             if(fr->use_acceleration)
1241             {
1242                 calc_gb_rad_hct_obc_sse2_single(cr,fr,born->nr,top, atype, x[0], nl, born, md, ir->gb_algorithm);
1243             }
1244             else
1245             {
1246                 calc_gb_rad_obc(cr,fr,born->nr,top,atype,x,nl,born,md);
1247             }
1248             break;
1249
1250         default:
1251             gmx_fatal(FARGS, "Unknown sse-enabled algorithm for Born radii calculation: %d",ir->gb_algorithm);
1252     }
1253
1254 #else
1255     switch(ir->gb_algorithm)
1256     {
1257         case egbSTILL:
1258             calc_gb_rad_still(cr,fr,born->nr,top,atype,x,nl,born,md);
1259             break;
1260         case egbHCT:
1261             calc_gb_rad_hct(cr,fr,born->nr,top,atype,x,nl,born,md);
1262             break;
1263         case egbOBC:
1264             calc_gb_rad_obc(cr,fr,born->nr,top,atype,x,nl,born,md);
1265             break;
1266
1267         default:
1268             gmx_fatal(FARGS, "Unknown algorithm for Born radii calculation: %d",ir->gb_algorithm);
1269     }
1270
1271 #endif /* Single precision sse */
1272
1273 #endif /* Double or single precision */
1274
1275     if(fr->bAllvsAll==FALSE)
1276     {
1277         switch(ir->gb_algorithm)
1278         {
1279             case egbSTILL:
1280                 /* 17 flops per outer loop iteration, 47 flops per inner loop */
1281                 inc_nrnb(nrnb,eNR_BORN_RADII_STILL,nl->nri*17+nl->nrj*47);
1282                 break;
1283             case egbHCT:
1284             case egbOBC:
1285                 /* 61 (assuming 10 for tanh) flops for outer loop iteration, 183 flops per inner loop */
1286                 inc_nrnb(nrnb,eNR_BORN_RADII_HCT_OBC,nl->nri*61+nl->nrj*183);
1287                 break;
1288
1289             default:
1290                 break;
1291         }
1292     }
1293
1294     return 0;
1295 }
1296
1297
1298
1299 real gb_bonds_tab(rvec x[], rvec f[], rvec fshift[], real *charge, real *p_gbtabscale,
1300                   real *invsqrta, real *dvda, real *GBtab, t_idef *idef, real epsilon_r,
1301                   real gb_epsilon_solvent, real facel, const t_pbc *pbc, const t_graph *graph)
1302 {
1303     int i,j,n0,m,nnn,type,ai,aj;
1304         int ki;
1305
1306         real isai,isaj;
1307     real r,rsq11;
1308     real rinv11,iq;
1309     real isaprod,qq,gbscale,gbtabscale,Y,F,Geps,Heps2,Fp,VV,FF,rt,eps,eps2;
1310     real vgb,fgb,vcoul,fijC,dvdatmp,fscal,dvdaj;
1311     real vctot;
1312
1313         rvec dx;
1314         ivec dt;
1315
1316     t_iatom *forceatoms;
1317
1318     /* Scale the electrostatics by gb_epsilon_solvent */
1319     facel = facel * ((1.0/epsilon_r) - 1.0/gb_epsilon_solvent);
1320
1321     gbtabscale=*p_gbtabscale;
1322     vctot = 0.0;
1323
1324     for(j=F_GB12;j<=F_GB14;j++)
1325     {
1326         forceatoms = idef->il[j].iatoms;
1327
1328         for(i=0;i<idef->il[j].nr; )
1329         {
1330             /* To avoid reading in the interaction type, we just increment i to pass over
1331              * the types in the forceatoms array, this saves some memory accesses
1332              */
1333             i++;
1334             ai            = forceatoms[i++];
1335             aj            = forceatoms[i++];
1336
1337                         ki            = pbc_rvec_sub(pbc,x[ai],x[aj],dx);
1338                         rsq11         = iprod(dx,dx);
1339
1340                         isai          = invsqrta[ai];
1341                         iq            = (-1)*facel*charge[ai];
1342
1343             rinv11        = gmx_invsqrt(rsq11);
1344             isaj          = invsqrta[aj];
1345             isaprod       = isai*isaj;
1346             qq            = isaprod*iq*charge[aj];
1347             gbscale       = isaprod*gbtabscale;
1348             r             = rsq11*rinv11;
1349             rt            = r*gbscale;
1350             n0            = rt;
1351             eps           = rt-n0;
1352             eps2          = eps*eps;
1353             nnn           = 4*n0;
1354             Y             = GBtab[nnn];
1355             F             = GBtab[nnn+1];
1356             Geps          = eps*GBtab[nnn+2];
1357             Heps2         = eps2*GBtab[nnn+3];
1358             Fp            = F+Geps+Heps2;
1359             VV            = Y+eps*Fp;
1360             FF            = Fp+Geps+2.0*Heps2;
1361             vgb           = qq*VV;
1362             fijC          = qq*FF*gbscale;
1363             dvdatmp       = -(vgb+fijC*r)*0.5;
1364             dvda[aj]      = dvda[aj] + dvdatmp*isaj*isaj;
1365             dvda[ai]      = dvda[ai] + dvdatmp*isai*isai;
1366             vctot         = vctot + vgb;
1367             fgb           = -(fijC)*rinv11;
1368
1369                         if (graph) {
1370                                 ivec_sub(SHIFT_IVEC(graph,ai),SHIFT_IVEC(graph,aj),dt);
1371                                 ki=IVEC2IS(dt);
1372                         }
1373
1374                         for (m=0; (m<DIM); m++) {                       /*  15          */
1375                                 fscal=fgb*dx[m];
1376                                 f[ai][m]+=fscal;
1377                                 f[aj][m]-=fscal;
1378                                 fshift[ki][m]+=fscal;
1379                                 fshift[CENTRAL][m]-=fscal;
1380                         }
1381         }
1382     }
1383
1384     return vctot;
1385 }
1386
1387 real calc_gb_selfcorrections(t_commrec *cr, int natoms,
1388                  real *charge, gmx_genborn_t *born, real *dvda, t_mdatoms *md, double facel)
1389 {
1390     int i,ai,at0,at1;
1391     real rai,e,derb,q,q2,fi,rai_inv,vtot;
1392
1393     if(PARTDECOMP(cr))
1394     {
1395         pd_at_range(cr,&at0,&at1);
1396     }
1397     else if(DOMAINDECOMP(cr))
1398     {
1399         at0=0;
1400         at1=cr->dd->nat_home;
1401     }
1402     else
1403     {
1404         at0=0;
1405         at1=natoms;
1406
1407     }
1408
1409     /* Scale the electrostatics by gb_epsilon_solvent */
1410     facel = facel * ((1.0/born->epsilon_r) - 1.0/born->gb_epsilon_solvent);
1411
1412     vtot=0.0;
1413
1414     /* Apply self corrections */
1415     for(i=at0;i<at1;i++)
1416     {
1417         ai       = i;
1418
1419         if(born->use[ai]==1)
1420         {
1421             rai      = born->bRad[ai];
1422             rai_inv  = 1.0/rai;
1423             q        = charge[ai];
1424             q2       = q*q;
1425             fi       = facel*q2;
1426             e        = fi*rai_inv;
1427             derb     = 0.5*e*rai_inv*rai_inv;
1428             dvda[ai] += derb*rai;
1429             vtot     -= 0.5*e;
1430         }
1431     }
1432
1433    return vtot;
1434
1435 }
1436
1437 real calc_gb_nonpolar(t_commrec *cr, t_forcerec *fr,int natoms,gmx_genborn_t *born, gmx_localtop_t *top,
1438                       const t_atomtypes *atype, real *dvda,int gb_algorithm, t_mdatoms *md)
1439 {
1440     int ai,i,at0,at1;
1441     real e,es,rai,rbi,term,probe,tmp,factor;
1442     real rbi_inv,rbi_inv2;
1443
1444     /* To keep the compiler happy */
1445     factor=0;
1446
1447     if(PARTDECOMP(cr))
1448     {
1449         pd_at_range(cr,&at0,&at1);
1450     }
1451     else if(DOMAINDECOMP(cr))
1452     {
1453         at0 = 0;
1454         at1 = cr->dd->nat_home;
1455     }
1456     else
1457     {
1458         at0=0;
1459         at1=natoms;
1460     }
1461
1462   /* factor is the surface tension */
1463   factor = born->sa_surface_tension;
1464   /*
1465
1466     // The surface tension factor is 0.0049 for Still model, 0.0054 for HCT/OBC
1467     if(gb_algorithm==egbSTILL)
1468     {
1469         factor=0.0049*100*CAL2JOULE;
1470     }
1471     else
1472     {
1473         factor=0.0054*100*CAL2JOULE;
1474     }
1475     */
1476     /* if(gb_algorithm==egbHCT || gb_algorithm==egbOBC) */
1477
1478     es    = 0;
1479     probe = 0.14;
1480     term  = M_PI*4;
1481
1482     for(i=at0;i<at1;i++)
1483     {
1484         ai        = i;
1485
1486         if(born->use[ai]==1)
1487         {
1488             rai       = top->atomtypes.gb_radius[md->typeA[ai]];
1489             rbi_inv   = fr->invsqrta[ai];
1490             rbi_inv2  = rbi_inv * rbi_inv;
1491             tmp       = (rai*rbi_inv2)*(rai*rbi_inv2);
1492             tmp       = tmp*tmp*tmp;
1493             e         = factor*term*(rai+probe)*(rai+probe)*tmp;
1494             dvda[ai]  = dvda[ai] - 6*e*rbi_inv2;
1495             es        = es + e;
1496         }
1497     }
1498
1499     return es;
1500 }
1501
1502
1503
1504 real calc_gb_chainrule(int natoms, t_nblist *nl, real *dadx, real *dvda, rvec x[], rvec t[], rvec fshift[],
1505                        rvec shift_vec[], int gb_algorithm, gmx_genborn_t *born, t_mdatoms *md)
1506 {
1507     int i,k,n,ai,aj,nj0,nj1,n0,n1;
1508     int shift;
1509     real shX,shY,shZ;
1510     real fgb,fij,rb2,rbi,fix1,fiy1,fiz1;
1511     real ix1,iy1,iz1,jx1,jy1,jz1,dx11,dy11,dz11,rsq11;
1512     real rinv11,tx,ty,tz,rbai,rbaj,fgb_ai;
1513     real *rb;
1514     volatile int idx;
1515
1516     n  = 0;
1517     rb = born->work;
1518
1519   n0 = 0;
1520   n1 = natoms;
1521
1522     if(gb_algorithm==egbSTILL)
1523     {
1524         for(i=n0;i<n1;i++)
1525         {
1526           rbi   = born->bRad[i];
1527           rb[i] = (2 * rbi * rbi * dvda[i])/ONE_4PI_EPS0;
1528         }
1529     }
1530     else if(gb_algorithm==egbHCT)
1531     {
1532         for(i=n0;i<n1;i++)
1533         {
1534           rbi   = born->bRad[i];
1535           rb[i] = rbi * rbi * dvda[i];
1536         }
1537     }
1538     else if(gb_algorithm==egbOBC)
1539     {
1540         for(i=n0;i<n1;i++)
1541         {
1542           rbi   = born->bRad[i];
1543           rb[i] = rbi * rbi * born->drobc[i] * dvda[i];
1544         }
1545     }
1546
1547     for(i=0;i<nl->nri;i++)
1548     {
1549         ai   = nl->iinr[i];
1550
1551         nj0  = nl->jindex[i];
1552         nj1  = nl->jindex[i+1];
1553
1554         /* Load shifts for this list */
1555         shift   = nl->shift[i];
1556         shX     = shift_vec[shift][0];
1557         shY     = shift_vec[shift][1];
1558         shZ     = shift_vec[shift][2];
1559
1560         /* Load atom i coordinates, add shift vectors */
1561         ix1  = shX + x[ai][0];
1562         iy1  = shY + x[ai][1];
1563         iz1  = shZ + x[ai][2];
1564
1565         fix1 = 0;
1566         fiy1 = 0;
1567         fiz1 = 0;
1568
1569         rbai = rb[ai];
1570
1571         for(k=nj0;k<nj1 && nl->jjnr[k]>=0;k++)
1572         {
1573             aj = nl->jjnr[k];
1574
1575             jx1     = x[aj][0];
1576             jy1     = x[aj][1];
1577             jz1     = x[aj][2];
1578
1579             dx11    = ix1 - jx1;
1580             dy11    = iy1 - jy1;
1581             dz11    = iz1 - jz1;
1582
1583             rbaj    = rb[aj];
1584
1585             fgb     = rbai*dadx[n++];
1586             fgb_ai  = rbaj*dadx[n++];
1587
1588             /* Total force between ai and aj is the sum of ai->aj and aj->ai */
1589             fgb     = fgb + fgb_ai;
1590
1591             tx      = fgb * dx11;
1592             ty      = fgb * dy11;
1593             tz      = fgb * dz11;
1594
1595             fix1    = fix1 + tx;
1596             fiy1    = fiy1 + ty;
1597             fiz1    = fiz1 + tz;
1598
1599             /* Update force on atom aj */
1600             t[aj][0] = t[aj][0] - tx;
1601             t[aj][1] = t[aj][1] - ty;
1602             t[aj][2] = t[aj][2] - tz;
1603         }
1604
1605         /* Update force and shift forces on atom ai */
1606         t[ai][0] = t[ai][0] + fix1;
1607         t[ai][1] = t[ai][1] + fiy1;
1608         t[ai][2] = t[ai][2] + fiz1;
1609
1610         fshift[shift][0] = fshift[shift][0] + fix1;
1611         fshift[shift][1] = fshift[shift][1] + fiy1;
1612         fshift[shift][2] = fshift[shift][2] + fiz1;
1613
1614     }
1615
1616     return 0;
1617 }
1618
1619
1620 void
1621 calc_gb_forces(t_commrec *cr, t_mdatoms *md, gmx_genborn_t *born, gmx_localtop_t *top, const t_atomtypes *atype,
1622                rvec x[], rvec f[], t_forcerec *fr, t_idef *idef, int gb_algorithm, int sa_algorithm, t_nrnb *nrnb, gmx_bool bRad,
1623                const t_pbc *pbc, const t_graph *graph, gmx_enerdata_t *enerd)
1624 {
1625     real v=0;
1626     int  cnt;
1627     int i;
1628
1629         /* PBC or not? */
1630         const t_pbc *pbc_null;
1631
1632         if (fr->bMolPBC)
1633                 pbc_null = pbc;
1634         else
1635                 pbc_null = NULL;
1636
1637   if(sa_algorithm == esaAPPROX)
1638   {
1639     /* Do a simple ACE type approximation for the non-polar solvation */
1640     enerd->term[F_NPSOLVATION] += calc_gb_nonpolar(cr, fr,born->nr, born, top, atype, fr->dvda, gb_algorithm,md);
1641   }
1642
1643   /* Calculate the bonded GB-interactions using either table or analytical formula */
1644     enerd->term[F_GBPOL]       += gb_bonds_tab(x,f,fr->fshift, md->chargeA,&(fr->gbtabscale),
1645                                      fr->invsqrta,fr->dvda,fr->gbtab.data,idef,born->epsilon_r,born->gb_epsilon_solvent, fr->epsfac, pbc_null, graph);
1646
1647     /* Calculate self corrections to the GB energies - currently only A state used! (FIXME) */
1648     enerd->term[F_GBPOL]       += calc_gb_selfcorrections(cr,born->nr,md->chargeA, born, fr->dvda, md, fr->epsfac);
1649
1650     /* If parallel, sum the derivative of the potential w.r.t the born radii */
1651     if(PARTDECOMP(cr))
1652     {
1653         gmx_sum(md->nr,fr->dvda, cr);
1654     }
1655     else if(DOMAINDECOMP(cr))
1656     {
1657         dd_atom_sum_real(cr->dd,fr->dvda);
1658         dd_atom_spread_real(cr->dd,fr->dvda);
1659     }
1660
1661     if(fr->bAllvsAll)
1662     {
1663 #if 0 && defined (GMX_X86_SSE2)
1664         if(fr->use_acceleration)
1665         {
1666 #  ifdef GMX_DOUBLE
1667             genborn_allvsall_calc_chainrule_sse2_double(fr,md,born,x[0],f[0],gb_algorithm,fr->AllvsAll_workgb);
1668 #  else
1669             genborn_allvsall_calc_chainrule_sse2_single(fr,md,born,x[0],f[0],gb_algorithm,fr->AllvsAll_workgb);
1670 #  endif
1671         }
1672         else
1673         {
1674             genborn_allvsall_calc_chainrule(fr,md,born,x[0],f[0],gb_algorithm,fr->AllvsAll_workgb);
1675         }
1676 #else
1677         genborn_allvsall_calc_chainrule(fr,md,born,x[0],f[0],gb_algorithm,fr->AllvsAll_workgb);
1678 #endif
1679         cnt = md->homenr*(md->nr/2+1);
1680         /* 9 flops for outer loop, 15 for inner */
1681         inc_nrnb(nrnb,eNR_BORN_AVA_CHAINRULE,md->homenr*9+cnt*15);
1682         return;
1683     }
1684
1685 #if 0 && defined (GMX_X86_SSE2)
1686     if(fr->use_acceleration)
1687     {
1688 #  ifdef GMX_DOUBLE
1689         calc_gb_chainrule_sse2_double(fr->natoms_force, &(fr->gblist),fr->dadx,fr->dvda,x[0],
1690                                       f[0],fr->fshift[0],fr->shift_vec[0],gb_algorithm,born,md);
1691 #  else
1692         calc_gb_chainrule_sse2_single(fr->natoms_force, &(fr->gblist),fr->dadx,fr->dvda,x[0],
1693                                       f[0],fr->fshift[0],fr->shift_vec[0],gb_algorithm,born,md);
1694 #  endif
1695     }
1696     else
1697     {
1698         calc_gb_chainrule(fr->natoms_force, &(fr->gblist), fr->dadx, fr->dvda,
1699                           x, f, fr->fshift, fr->shift_vec, gb_algorithm, born, md);
1700     }
1701 #else
1702     calc_gb_chainrule(fr->natoms_force, &(fr->gblist), fr->dadx, fr->dvda,
1703                       x, f, fr->fshift, fr->shift_vec, gb_algorithm, born, md);
1704 #endif
1705
1706     if(!fr->bAllvsAll)
1707     {
1708         /* 9 flops for outer loop, 15 for inner */
1709         inc_nrnb(nrnb,eNR_BORN_CHAINRULE,fr->gblist.nri*9+fr->gblist.nrj*15);
1710     }
1711 }
1712
1713 static void add_j_to_gblist(gbtmpnbl_t *list,int aj)
1714 {
1715     if (list->naj >= list->aj_nalloc)
1716     {
1717         list->aj_nalloc = over_alloc_large(list->naj+1);
1718         srenew(list->aj,list->aj_nalloc);
1719     }
1720
1721     list->aj[list->naj++] = aj;
1722 }
1723
1724 static gbtmpnbl_t *find_gbtmplist(struct gbtmpnbls *lists,int shift)
1725 {
1726     int ind,i;
1727
1728     /* Search the list with the same shift, if there is one */
1729     ind = 0;
1730     while (ind < lists->nlist && shift != lists->list[ind].shift)
1731     {
1732         ind++;
1733     }
1734     if (ind == lists->nlist)
1735     {
1736         if (lists->nlist == lists->list_nalloc)
1737         {
1738             lists->list_nalloc++;
1739             srenew(lists->list,lists->list_nalloc);
1740             for(i=lists->nlist; i<lists->list_nalloc; i++)
1741             {
1742                 lists->list[i].aj        = NULL;
1743                 lists->list[i].aj_nalloc = 0;
1744             }
1745
1746         }
1747
1748         lists->list[lists->nlist].shift = shift;
1749         lists->list[lists->nlist].naj   = 0;
1750         lists->nlist++;
1751     }
1752
1753     return &lists->list[ind];
1754 }
1755
1756 static void add_bondeds_to_gblist(t_ilist *il,
1757                                   gmx_bool bMolPBC,t_pbc *pbc,t_graph *g,rvec *x,
1758                                   struct gbtmpnbls *nls)
1759 {
1760     int  ind,j,ai,aj,shift,found;
1761     rvec dx;
1762     ivec dt;
1763     gbtmpnbl_t *list;
1764
1765     shift = CENTRAL;
1766     for(ind=0; ind<il->nr; ind+=3)
1767     {
1768         ai = il->iatoms[ind+1];
1769         aj = il->iatoms[ind+2];
1770
1771         shift = CENTRAL;
1772         if (g != NULL)
1773         {
1774           rvec_sub(x[ai],x[aj],dx);
1775           ivec_sub(SHIFT_IVEC(g,ai),SHIFT_IVEC(g,aj),dt);
1776           shift = IVEC2IS(dt);
1777         }
1778         else if (bMolPBC)
1779         {
1780           shift = pbc_dx_aiuc(pbc,x[ai],x[aj],dx);
1781         }
1782
1783         /* Find the list for this shift or create one */
1784         list = find_gbtmplist(&nls[ai],shift);
1785
1786         found=0;
1787
1788         /* So that we do not add the same bond twice.
1789          * This happens with some constraints between 1-3 atoms
1790          * that are in the bond-list but should not be in the GB nb-list */
1791         for(j=0;j<list->naj;j++)
1792         {
1793             if (list->aj[j] == aj)
1794             {
1795                 found = 1;
1796             }
1797         }
1798
1799         if (found == 0)
1800         {
1801                         if(ai == aj)
1802                         {
1803                                 gmx_incons("ai == aj");
1804                         }
1805
1806             add_j_to_gblist(list,aj);
1807         }
1808     }
1809 }
1810
1811 static int
1812 compare_int (const void * a, const void * b)
1813 {
1814     return ( *(int*)a - *(int*)b );
1815 }
1816
1817
1818
1819 int make_gb_nblist(t_commrec *cr, int gb_algorithm, real gbcut,
1820                    rvec x[], matrix box,
1821                    t_forcerec *fr, t_idef *idef, t_graph *graph, gmx_genborn_t *born)
1822 {
1823     int i,l,ii,j,k,n,nj0,nj1,ai,aj,at0,at1,found,shift,s;
1824     int apa;
1825     t_nblist *nblist;
1826     t_pbc pbc;
1827
1828     struct gbtmpnbls *nls;
1829     gbtmpnbl_t *list =NULL;
1830
1831     set_pbc(&pbc,fr->ePBC,box);
1832     nls   = born->nblist_work;
1833
1834     for(i=0;i<born->nr;i++)
1835     {
1836         nls[i].nlist = 0;
1837     }
1838
1839     if (fr->bMolPBC)
1840     {
1841         set_pbc_dd(&pbc,fr->ePBC,cr->dd,TRUE,box);
1842     }
1843
1844     switch (gb_algorithm)
1845     {
1846     case egbHCT:
1847     case egbOBC:
1848         /* Loop over 1-2, 1-3 and 1-4 interactions */
1849         for(j=F_GB12;j<=F_GB14;j++)
1850         {
1851             add_bondeds_to_gblist(&idef->il[j],fr->bMolPBC,&pbc,graph,x,nls);
1852         }
1853         break;
1854     case egbSTILL:
1855         /* Loop over 1-4 interactions */
1856         add_bondeds_to_gblist(&idef->il[F_GB14],fr->bMolPBC,&pbc,graph,x,nls);
1857         break;
1858     default:
1859         gmx_incons("Unknown GB algorithm");
1860     }
1861
1862     /* Loop over the VDWQQ and VDW nblists to set up the nonbonded part of the GB list */
1863     for(n=0; (n<fr->nnblists); n++)
1864     {
1865         for(i=0; (i<eNL_NR); i++)
1866         {
1867             nblist=&(fr->nblists[n].nlist_sr[i]);
1868
1869             if (nblist->nri > 0 && (i==eNL_VDWQQ || i==eNL_QQ))
1870             {
1871                 for(j=0;j<nblist->nri;j++)
1872                 {
1873                     ai    = nblist->iinr[j];
1874                     shift = nblist->shift[j];
1875
1876                     /* Find the list for this shift or create one */
1877                     list = find_gbtmplist(&nls[ai],shift);
1878
1879                     nj0 = nblist->jindex[j];
1880                     nj1 = nblist->jindex[j+1];
1881
1882                     /* Add all the j-atoms in the non-bonded list to the GB list */
1883                     for(k=nj0;k<nj1;k++)
1884                     {
1885                         add_j_to_gblist(list,nblist->jjnr[k]);
1886                     }
1887                 }
1888             }
1889         }
1890     }
1891
1892     /* Zero out some counters */
1893         fr->gblist.nri=0;
1894     fr->gblist.nrj=0;
1895
1896         fr->gblist.jindex[0] = fr->gblist.nri;
1897
1898         for(i=0;i<fr->natoms_force;i++)
1899     {
1900         for(s=0; s<nls[i].nlist; s++)
1901         {
1902             list = &nls[i].list[s];
1903
1904             /* Only add those atoms that actually have neighbours */
1905             if (born->use[i] != 0)
1906             {
1907                 fr->gblist.iinr[fr->gblist.nri]  = i;
1908                 fr->gblist.shift[fr->gblist.nri] = list->shift;
1909                 fr->gblist.nri++;
1910
1911                 for(k=0; k<list->naj; k++)
1912                 {
1913                     /* Memory allocation for jjnr */
1914                     if(fr->gblist.nrj >= fr->gblist.maxnrj)
1915                     {
1916                         fr->gblist.maxnrj += over_alloc_large(fr->gblist.maxnrj);
1917
1918                         if (debug)
1919                         {
1920                             fprintf(debug,"Increasing GB neighbourlist j size to %d\n",fr->gblist.maxnrj);
1921                         }
1922
1923                         srenew(fr->gblist.jjnr,fr->gblist.maxnrj);
1924                     }
1925
1926                     /* Put in list */
1927                                         if(i == list->aj[k])
1928                                         {
1929                                                 gmx_incons("i == list->aj[k]");
1930                                         }
1931                     fr->gblist.jjnr[fr->gblist.nrj++] = list->aj[k];
1932                 }
1933
1934                                 fr->gblist.jindex[fr->gblist.nri] = fr->gblist.nrj;
1935             }
1936                 }
1937         }
1938
1939
1940 #ifdef SORT_GB_LIST
1941     for(i=0;i<fr->gblist.nri;i++)
1942     {
1943         nj0 = fr->gblist.jindex[i];
1944         nj1 = fr->gblist.jindex[i+1];
1945         ai  = fr->gblist.iinr[i];
1946
1947         /* Temporary fix */
1948                 for(j=nj0;j<nj1;j++)
1949                 {
1950             if(fr->gblist.jjnr[j]<ai)
1951                 fr->gblist.jjnr[j]+=fr->natoms_force;
1952         }
1953         qsort(fr->gblist.jjnr+nj0,nj1-nj0,sizeof(int),compare_int);
1954         /* Fix back */
1955         for(j=nj0;j<nj1;j++)
1956         {
1957             if(fr->gblist.jjnr[j]>=fr->natoms_force)
1958                 fr->gblist.jjnr[j]-=fr->natoms_force;
1959         }
1960
1961     }
1962 #endif
1963
1964     return 0;
1965 }
1966
1967 void make_local_gb(const t_commrec *cr, gmx_genborn_t *born, int gb_algorithm)
1968 {
1969     int i,at0,at1;
1970     gmx_domdec_t *dd=NULL;
1971
1972     if(DOMAINDECOMP(cr))
1973     {
1974         dd = cr->dd;
1975         at0 = 0;
1976         at1 = dd->nat_tot;
1977     }
1978     else
1979     {
1980         /* Single node or particle decomp (global==local), just copy pointers and return */
1981         if(gb_algorithm==egbSTILL)
1982         {
1983             born->gpol      = born->gpol_globalindex;
1984             born->vsolv     = born->vsolv_globalindex;
1985             born->gb_radius = born->gb_radius_globalindex;
1986         }
1987         else
1988         {
1989             born->param     = born->param_globalindex;
1990             born->gb_radius = born->gb_radius_globalindex;
1991         }
1992
1993         born->use = born->use_globalindex;
1994
1995         return;
1996     }
1997
1998     /* Reallocation of local arrays if necessary */
1999     /* fr->natoms_force is equal to dd->nat_tot */
2000     if (DOMAINDECOMP(cr) && dd->nat_tot > born->nalloc)
2001     {
2002         int nalloc;
2003
2004         nalloc = dd->nat_tot;
2005
2006         /* Arrays specific to different gb algorithms */
2007         if (gb_algorithm == egbSTILL)
2008         {
2009             srenew(born->gpol,  nalloc+3);
2010             srenew(born->vsolv, nalloc+3);
2011             srenew(born->gb_radius, nalloc+3);
2012             for(i=born->nalloc; (i<nalloc+3); i++)
2013             {
2014                 born->gpol[i] = 0;
2015                 born->vsolv[i] = 0;
2016                 born->gb_radius[i] = 0;
2017             }
2018         }
2019         else
2020         {
2021             srenew(born->param, nalloc+3);
2022             srenew(born->gb_radius, nalloc+3);
2023             for(i=born->nalloc; (i<nalloc+3); i++)
2024             {
2025                 born->param[i] = 0;
2026                 born->gb_radius[i] = 0;
2027             }
2028         }
2029
2030         /* All gb-algorithms use the array for vsites exclusions */
2031         srenew(born->use,    nalloc+3);
2032         for(i=born->nalloc; (i<nalloc+3); i++)
2033         {
2034             born->use[i] = 0;
2035         }
2036
2037         born->nalloc = nalloc;
2038     }
2039
2040     /* With dd, copy algorithm specific arrays */
2041     if(gb_algorithm==egbSTILL)
2042     {
2043         for(i=at0;i<at1;i++)
2044         {
2045             born->gpol[i]  = born->gpol_globalindex[dd->gatindex[i]];
2046             born->vsolv[i] = born->vsolv_globalindex[dd->gatindex[i]];
2047             born->gb_radius[i] = born->gb_radius_globalindex[dd->gatindex[i]];
2048             born->use[i]   = born->use_globalindex[dd->gatindex[i]];
2049         }
2050     }
2051     else
2052     {
2053         for(i=at0;i<at1;i++)
2054         {
2055             born->param[i]     = born->param_globalindex[dd->gatindex[i]];
2056             born->gb_radius[i] = born->gb_radius_globalindex[dd->gatindex[i]];
2057             born->use[i]       = born->use_globalindex[dd->gatindex[i]];
2058         }
2059     }
2060 }
2061