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