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