src/gromacs/mdlib/genborn_allvsall.c

   1 /*
   2  *                This source code is part of
   3  *
   4  *                 G   R   O   M   A   C   S
   5  *
   6  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   7  * Copyright (c) 2001-2009, The GROMACS Development Team
   8  *
   9  * Gromacs is a library for molecular simulation and trajectory analysis,
  10  * written by Erik Lindahl, David van der Spoel, Berk Hess, and others - for
  11  * a full list of developers and information, check out http://www.gromacs.org
  12  *
  13  * This program is free software; you can redistribute it and/or modify it under
  14  * the terms of the GNU Lesser General Public License as published by the Free
  15  * Software Foundation; either version 2 of the License, or (at your option) any
  16  * later version.
  17  * As a special exception, you may use this file as part of a free software
  18  * library without restriction.  Specifically, if other files instantiate
  19  * templates or use macros or inline functions from this file, or you compile
  20  * this file and link it with other files to produce an executable, this
  21  * file does not by itself cause the resulting executable to be covered by
  22  * the GNU Lesser General Public License.
  23  *
  24  * In plain-speak: do not worry about classes/macros/templates either - only
  25  * changes to the library have to be LGPL, not an application linking with it.
  26  *
  27  * To help fund GROMACS development, we humbly ask that you cite
  28  * the papers people have written on it - you can find them on the website!
  29  */
  30 #ifdef HAVE_CONFIG_H
  31 #include <config.h>
  32 #endif
  33
  34 #include <math.h>
  35 #include "types/simple.h"
  36
  37 #include "vec.h"
  38 #include "smalloc.h"
  39
  40 #include "partdec.h"
  41 #include "network.h"
  42 #include "physics.h"
  43 #include "genborn.h"
  44 #include "genborn_allvsall.h"
  45
  46
  47 typedef struct
  48 {
  49     int *      jindex_gb;
  50     int **     exclusion_mask_gb;
  51 }
  52 gmx_allvsallgb2_data_t;
  53
  54 static int
  55 calc_maxoffset(int i, int natoms)
  56 {
  57     int maxoffset;
  58
  59     if ((natoms % 2) == 1)
  60     {
  61         /* Odd number of atoms, easy */
  62         maxoffset = natoms/2;
  63     }
  64     else if ((natoms % 4) == 0)
  65     {
  66         /* Multiple of four is hard */
  67         if (i < natoms/2)
  68         {
  69             if ((i % 2) == 0)
  70             {
  71                 maxoffset = natoms/2;
  72             }
  73             else
  74             {
  75                 maxoffset = natoms/2-1;
  76             }
  77         }
  78         else
  79         {
  80             if ((i % 2) == 1)
  81             {
  82                 maxoffset = natoms/2;
  83             }
  84             else
  85             {
  86                 maxoffset = natoms/2-1;
  87             }
  88         }
  89     }
  90     else
  91     {
  92         /* natoms/2 = odd */
  93         if ((i % 2) == 0)
  94         {
  95             maxoffset = natoms/2;
  96         }
  97         else
  98         {
  99             maxoffset = natoms/2-1;
 100         }
 101     }
 102
 103     return maxoffset;
 104 }
 105
 106 static void
 107 setup_gb_exclusions_and_indices(gmx_allvsallgb2_data_t     *   aadata,
 108                                 t_ilist     *                  ilist,
 109                                 int                            natoms,
 110                                 gmx_bool                       bInclude12,
 111                                 gmx_bool                       bInclude13,
 112                                 gmx_bool                       bInclude14)
 113 {
 114     int i, j, k, tp;
 115     int a1, a2;
 116     int nj0, nj1;
 117     int max_offset;
 118     int max_excl_offset;
 119     int nj;
 120
 121     /* This routine can appear to be a bit complex, but it is mostly book-keeping.
 122      * To enable the fast all-vs-all kernel we need to be able to stream through all coordinates
 123      * whether they should interact or not.
 124      *
 125      * To avoid looping over the exclusions, we create a simple mask that is 1 if the interaction
 126      * should be present, otherwise 0. Since exclusions typically only occur when i & j are close,
 127      * we create a jindex array with three elements per i atom: the starting point, the point to
 128      * which we need to check exclusions, and the end point.
 129      * This way we only have to allocate a short exclusion mask per i atom.
 130      */
 131
 132     /* Allocate memory for jindex arrays */
 133     snew(aadata->jindex_gb, 3*natoms);
 134
 135     /* Pointer to lists with exclusion masks */
 136     snew(aadata->exclusion_mask_gb, natoms);
 137
 138     for (i = 0; i < natoms; i++)
 139     {
 140         /* Start */
 141         aadata->jindex_gb[3*i]       = i+1;
 142         max_offset                   = calc_maxoffset(i, natoms);
 143
 144         /* first check the max range of atoms to EXCLUDE */
 145         max_excl_offset = 0;
 146         if (!bInclude12)
 147         {
 148             for (j = 0; j < ilist[F_GB12].nr; j += 3)
 149             {
 150                 a1 = ilist[F_GB12].iatoms[j+1];
 151                 a2 = ilist[F_GB12].iatoms[j+2];
 152
 153                 if (a1 == i)
 154                 {
 155                     k = a2-a1;
 156                 }
 157                 else if (a2 == i)
 158                 {
 159                     k = a1+natoms-a2;
 160                 }
 161                 else
 162                 {
 163                     continue;
 164                 }
 165                 if (k > 0 && k <= max_offset)
 166                 {
 167                     max_excl_offset = (k > max_excl_offset) ? k : max_excl_offset;
 168                 }
 169             }
 170         }
 171         if (!bInclude13)
 172         {
 173             for (j = 0; j < ilist[F_GB13].nr; j += 3)
 174             {
 175                 a1 = ilist[F_GB13].iatoms[j+1];
 176                 a2 = ilist[F_GB13].iatoms[j+2];
 177
 178
 179                 if (a1 == i)
 180                 {
 181                     k = a2-a1;
 182                 }
 183                 else if (a2 == i)
 184                 {
 185                     k = a1+natoms-a2;
 186                 }
 187                 else
 188                 {
 189                     continue;
 190                 }
 191                 if (k > 0 && k <= max_offset)
 192                 {
 193                     max_excl_offset = (k > max_excl_offset) ? k : max_excl_offset;
 194                 }
 195             }
 196         }
 197         if (!bInclude14)
 198         {
 199             for (j = 0; j < ilist[F_GB14].nr; j += 3)
 200             {
 201                 a1 = ilist[F_GB14].iatoms[j+1];
 202                 a2 = ilist[F_GB14].iatoms[j+2];
 203
 204
 205                 if (a1 == i)
 206                 {
 207                     k = a2-a1;
 208                 }
 209                 else if (a2 == i)
 210                 {
 211                     k = a1+natoms-a2;
 212                 }
 213                 else
 214                 {
 215                     continue;
 216                 }
 217                 if (k > 0 && k <= max_offset)
 218                 {
 219                     max_excl_offset = (k > max_excl_offset) ? k : max_excl_offset;
 220                 }
 221             }
 222         }
 223         max_excl_offset = (max_offset < max_excl_offset) ? max_offset : max_excl_offset;
 224
 225         aadata->jindex_gb[3*i+1] = i+1+max_excl_offset;
 226
 227         snew(aadata->exclusion_mask_gb[i], max_excl_offset);
 228
 229         /* Include everything by default */
 230         for (j = 0; j < max_excl_offset; j++)
 231         {
 232             /* Use all-ones to mark interactions that should be present, compatible with SSE */
 233             aadata->exclusion_mask_gb[i][j] = 0xFFFFFFFF;
 234         }
 235         /* Go through exclusions again */
 236         if (!bInclude12)
 237         {
 238             for (j = 0; j < ilist[F_GB12].nr; j += 3)
 239             {
 240                 a1 = ilist[F_GB12].iatoms[j+1];
 241                 a2 = ilist[F_GB12].iatoms[j+2];
 242
 243                 if (a1 == i)
 244                 {
 245                     k = a2-a1;
 246                 }
 247                 else if (a2 == i)
 248                 {
 249                     k = a1+natoms-a2;
 250                 }
 251                 else
 252                 {
 253                     continue;
 254                 }
 255                 if (k > 0 && k <= max_offset)
 256                 {
 257                     aadata->exclusion_mask_gb[i][k-1] = 0;
 258                 }
 259             }
 260         }
 261         if (!bInclude13)
 262         {
 263             for (j = 0; j < ilist[F_GB13].nr; j += 3)
 264             {
 265                 a1 = ilist[F_GB13].iatoms[j+1];
 266                 a2 = ilist[F_GB13].iatoms[j+2];
 267
 268                 if (a1 == i)
 269                 {
 270                     k = a2-a1;
 271                 }
 272                 else if (a2 == i)
 273                 {
 274                     k = a1+natoms-a2;
 275                 }
 276                 else
 277                 {
 278                     continue;
 279                 }
 280                 if (k > 0 && k <= max_offset)
 281                 {
 282                     aadata->exclusion_mask_gb[i][k-1] = 0;
 283                 }
 284             }
 285         }
 286         if (!bInclude14)
 287         {
 288             for (j = 0; j < ilist[F_GB14].nr; j += 3)
 289             {
 290                 a1 = ilist[F_GB14].iatoms[j+1];
 291                 a2 = ilist[F_GB14].iatoms[j+2];
 292
 293                 if (a1 == i)
 294                 {
 295                     k = a2-a1;
 296                 }
 297                 else if (a2 == i)
 298                 {
 299                     k = a1+natoms-a2;
 300                 }
 301                 else
 302                 {
 303                     continue;
 304                 }
 305                 if (k > 0 && k <= max_offset)
 306                 {
 307                     aadata->exclusion_mask_gb[i][k-1] = 0;
 308                 }
 309             }
 310         }
 311
 312         /* End */
 313
 314         /* End */
 315         aadata->jindex_gb[3*i+2] = i+1+max_offset;
 316     }
 317 }
 318
 319
 320 static void
 321 genborn_allvsall_setup(gmx_allvsallgb2_data_t     **  p_aadata,
 322                        t_ilist     *                  ilist,
 323                        int                            natoms,
 324                        gmx_bool                       bInclude12,
 325                        gmx_bool                       bInclude13,
 326                        gmx_bool                       bInclude14)
 327 {
 328     int                     i, j, idx;
 329     gmx_allvsallgb2_data_t *aadata;
 330     real                   *p;
 331
 332     snew(aadata, 1);
 333     *p_aadata = aadata;
 334
 335     setup_gb_exclusions_and_indices(aadata, ilist, natoms, bInclude12, bInclude13, bInclude14);
 336 }
 337
 338
 339
 340 int
 341 genborn_allvsall_calc_still_radii(t_forcerec *           fr,
 342                                   t_mdatoms *            mdatoms,
 343                                   gmx_genborn_t *        born,
 344                                   gmx_localtop_t *       top,
 345                                   real *                 x,
 346                                   t_commrec *            cr,
 347                                   void *                 work)
 348 {
 349     gmx_allvsallgb2_data_t *aadata;
 350     int                     natoms;
 351     int                     ni0, ni1;
 352     int                     nj0, nj1, nj2;
 353     int                     i, j, k, n;
 354     int              *      mask;
 355
 356     real                    ix, iy, iz;
 357     real                    jx, jy, jz;
 358     real                    dx, dy, dz;
 359     real                    rsq, rinv;
 360     real                    gpi, rai, vai;
 361     real                    prod_ai;
 362     real                    irsq, idr4, idr6;
 363     real                    raj, rvdw, ratio;
 364     real                    vaj, ccf, dccf, theta, cosq;
 365     real                    term, prod, icf4, icf6, gpi2, factor, sinq;
 366
 367     natoms              = mdatoms->nr;
 368     ni0                 = mdatoms->start;
 369     ni1                 = mdatoms->start+mdatoms->homenr;
 370     factor              = 0.5*ONE_4PI_EPS0;
 371     n                   = 0;
 372
 373     aadata = *((gmx_allvsallgb2_data_t **)work);
 374
 375     if (aadata == NULL)
 376     {
 377         genborn_allvsall_setup(&aadata, top->idef.il, mdatoms->nr,
 378                                FALSE, FALSE, TRUE);
 379         *((gmx_allvsallgb2_data_t **)work) = aadata;
 380     }
 381
 382
 383     for (i = 0; i < born->nr; i++)
 384     {
 385         born->gpol_still_work[i] = 0;
 386     }
 387
 388
 389     for (i = ni0; i < ni1; i++)
 390     {
 391         /* We assume shifts are NOT used for all-vs-all interactions */
 392
 393         /* Load i atom data */
 394         ix                = x[3*i];
 395         iy                = x[3*i+1];
 396         iz                = x[3*i+2];
 397
 398         gpi               = 0.0;
 399
 400         rai     = top->atomtypes.gb_radius[mdatoms->typeA[i]];
 401         vai     = born->vsolv[i];
 402         prod_ai = STILL_P4*vai;
 403
 404         /* Load limits for loop over neighbors */
 405         nj0              = aadata->jindex_gb[3*i];
 406         nj1              = aadata->jindex_gb[3*i+1];
 407         nj2              = aadata->jindex_gb[3*i+2];
 408
 409         mask             = aadata->exclusion_mask_gb[i];
 410
 411         /* Prologue part, including exclusion mask */
 412         for (j = nj0; j < nj1; j++, mask++)
 413         {
 414             if (*mask != 0)
 415             {
 416                 k = j%natoms;
 417
 418                 /* load j atom coordinates */
 419                 jx                = x[3*k];
 420                 jy                = x[3*k+1];
 421                 jz                = x[3*k+2];
 422
 423                 /* Calculate distance */
 424                 dx                = ix - jx;
 425                 dy                = iy - jy;
 426                 dz                = iz - jz;
 427                 rsq               = dx*dx+dy*dy+dz*dz;
 428
 429                 /* Calculate 1/r and 1/r2 */
 430                 rinv              = gmx_invsqrt(rsq);
 431                 irsq              = rinv*rinv;
 432                 idr4              = irsq*irsq;
 433                 idr6              = idr4*irsq;
 434
 435                 raj = top->atomtypes.gb_radius[mdatoms->typeA[k]];
 436
 437                 rvdw  = rai + raj;
 438
 439                 ratio = rsq / (rvdw * rvdw);
 440                 vaj   = born->vsolv[k];
 441
 442
 443                 if (ratio > STILL_P5INV)
 444                 {
 445                     ccf  = 1.0;
 446                     dccf = 0.0;
 447                 }
 448                 else
 449                 {
 450                     theta = ratio*STILL_PIP5;
 451                     cosq  = cos(theta);
 452                     term  = 0.5*(1.0-cosq);
 453                     ccf   = term*term;
 454                     sinq  = 1.0 - cosq*cosq;
 455                     dccf  = 2.0*term*sinq*gmx_invsqrt(sinq)*theta;
 456                 }
 457
 458                 prod          = STILL_P4*vaj;
 459                 icf4          = ccf*idr4;
 460                 icf6          = (4*ccf-dccf)*idr6;
 461
 462                 born->gpol_still_work[k] += prod_ai*icf4;
 463                 gpi                       = gpi+prod*icf4;
 464
 465                 /* Save ai->aj and aj->ai chain rule terms */
 466                 fr->dadx[n++]   = prod*icf6;
 467                 fr->dadx[n++]   = prod_ai*icf6;
 468
 469                 /* 27 flops, plus one cos(x) - estimate at 20 flops  => 47 */
 470
 471             }
 472         }
 473
 474         /* Main part, no exclusions */
 475         for (j = nj1; j < nj2; j++)
 476         {
 477             k = j%natoms;
 478
 479             /* load j atom coordinates */
 480             jx                = x[3*k];
 481             jy                = x[3*k+1];
 482             jz                = x[3*k+2];
 483
 484             /* Calculate distance */
 485             dx                = ix - jx;
 486             dy                = iy - jy;
 487             dz                = iz - jz;
 488             rsq               = dx*dx+dy*dy+dz*dz;
 489
 490             /* Calculate 1/r and 1/r2 */
 491             rinv              = gmx_invsqrt(rsq);
 492             irsq              = rinv*rinv;
 493             idr4              = irsq*irsq;
 494             idr6              = idr4*irsq;
 495
 496             raj = top->atomtypes.gb_radius[mdatoms->typeA[k]];
 497
 498             rvdw  = rai + raj;
 499
 500             ratio = rsq / (rvdw * rvdw);
 501             vaj   = born->vsolv[k];
 502
 503             if (ratio > STILL_P5INV)
 504             {
 505                 ccf  = 1.0;
 506                 dccf = 0.0;
 507             }
 508             else
 509             {
 510                 theta = ratio*STILL_PIP5;
 511                 cosq  = cos(theta);
 512                 term  = 0.5*(1.0-cosq);
 513                 ccf   = term*term;
 514                 sinq  = 1.0 - cosq*cosq;
 515                 dccf  = 2.0*term*sinq*gmx_invsqrt(sinq)*theta;
 516             }
 517
 518             prod          = STILL_P4*vaj;
 519             icf4          = ccf*idr4;
 520             icf6          = (4*ccf-dccf)*idr6;
 521
 522             born->gpol_still_work[k] += prod_ai*icf4;
 523             gpi                       = gpi+prod*icf4;
 524
 525             /* Save ai->aj and aj->ai chain rule terms */
 526             fr->dadx[n++]   = prod*icf6;
 527             fr->dadx[n++]   = prod_ai*icf6;
 528         }
 529         born->gpol_still_work[i] += gpi;
 530     }
 531
 532     /* Parallel summations */
 533     if (PARTDECOMP(cr))
 534     {
 535         gmx_sum(natoms, born->gpol_still_work, cr);
 536     }
 537
 538     /* Calculate the radii */
 539     for (i = 0; i < natoms; i++)
 540     {
 541         if (born->use[i] != 0)
 542         {
 543             gpi             = born->gpol[i]+born->gpol_still_work[i];
 544             gpi2            = gpi * gpi;
 545             born->bRad[i]   = factor*gmx_invsqrt(gpi2);
 546             fr->invsqrta[i] = gmx_invsqrt(born->bRad[i]);
 547         }
 548     }
 549
 550     return 0;
 551 }
 552
 553
 554
 555 int
 556 genborn_allvsall_calc_hct_obc_radii(t_forcerec *           fr,
 557                                     t_mdatoms *            mdatoms,
 558                                     gmx_genborn_t *        born,
 559                                     int                    gb_algorithm,
 560                                     gmx_localtop_t *       top,
 561                                     real *                 x,
 562                                     t_commrec *            cr,
 563                                     void *                 work)
 564 {
 565     gmx_allvsallgb2_data_t *aadata;
 566     int                     natoms;
 567     int                     ni0, ni1;
 568     int                     nj0, nj1, nj2;
 569     int                     i, j, k, n;
 570     int              *      mask;
 571
 572     real                    ix, iy, iz;
 573     real                    jx, jy, jz;
 574     real                    dx, dy, dz;
 575     real                    rsq, rinv;
 576     real                    prod, raj;
 577     real                    rai, doffset, rai_inv, rai_inv2, sk_ai, sk2_ai, sum_ai;
 578     real                    dr, sk, lij, dlij, lij2, lij3, uij2, uij3, diff2, uij, log_term;
 579     real                    lij_inv, sk2, sk2_rinv, tmp, t1, t2, t3, raj_inv, sum_ai2, sum_ai3, tsum;
 580     real                    tchain;
 581     real                    dadxi, dadxj;
 582     real                    rad, min_rad;
 583
 584     natoms              = mdatoms->nr;
 585     ni0                 = mdatoms->start;
 586     ni1                 = mdatoms->start+mdatoms->homenr;
 587
 588     n       = 0;
 589     prod    = 0;
 590     raj     = 0;
 591     doffset = born->gb_doffset;
 592
 593     aadata = *((gmx_allvsallgb2_data_t **)work);
 594
 595     if (aadata == NULL)
 596     {
 597         genborn_allvsall_setup(&aadata, top->idef.il, mdatoms->nr,
 598                                TRUE, TRUE, TRUE);
 599         *((gmx_allvsallgb2_data_t **)work) = aadata;
 600     }
 601
 602     for (i = 0; i < born->nr; i++)
 603     {
 604         born->gpol_hct_work[i] = 0;
 605     }
 606
 607     for (i = ni0; i < ni1; i++)
 608     {
 609         /* We assume shifts are NOT used for all-vs-all interactions */
 610
 611         /* Load i atom data */
 612         ix                = x[3*i];
 613         iy                = x[3*i+1];
 614         iz                = x[3*i+2];
 615
 616         rai      = top->atomtypes.gb_radius[mdatoms->typeA[i]]-doffset;
 617         rai_inv  = 1.0/rai;
 618
 619         sk_ai    = born->param[i];
 620         sk2_ai   = sk_ai*sk_ai;
 621
 622         sum_ai   = 0;
 623
 624         /* Load limits for loop over neighbors */
 625         nj0              = aadata->jindex_gb[3*i];
 626         nj1              = aadata->jindex_gb[3*i+1];
 627         nj2              = aadata->jindex_gb[3*i+2];
 628
 629         mask             = aadata->exclusion_mask_gb[i];
 630
 631         /* Prologue part, including exclusion mask */
 632         for (j = nj0; j < nj1; j++, mask++)
 633         {
 634             if (*mask != 0)
 635             {
 636                 k = j%natoms;
 637
 638                 /* load j atom coordinates */
 639                 jx                = x[3*k];
 640                 jy                = x[3*k+1];
 641                 jz                = x[3*k+2];
 642
 643                 /* Calculate distance */
 644                 dx                = ix - jx;
 645                 dy                = iy - jy;
 646                 dz                = iz - jz;
 647                 rsq               = dx*dx+dy*dy+dz*dz;
 648
 649                 /* Calculate 1/r and 1/r2 */
 650                 rinv              = gmx_invsqrt(rsq);
 651                 dr                = rsq*rinv;
 652
 653                 /* sk is precalculated in init_gb() */
 654                 sk    = born->param[k];
 655                 raj   = top->atomtypes.gb_radius[mdatoms->typeA[k]]-doffset;
 656
 657                 /* aj -> ai interaction */
 658
 659
 660                 if (rai < dr+sk)
 661                 {
 662                     lij       = 1.0/(dr-sk);
 663                     dlij      = 1.0;
 664
 665                     if (rai > dr-sk)
 666                     {
 667                         lij  = rai_inv;
 668                         dlij = 0.0;
 669                     }
 670
 671                     uij      = 1.0/(dr+sk);
 672                     lij2     = lij  * lij;
 673                     lij3     = lij2 * lij;
 674                     uij2     = uij  * uij;
 675                     uij3     = uij2 * uij;
 676
 677                     diff2    = uij2-lij2;
 678
 679                     lij_inv  = gmx_invsqrt(lij2);
 680                     sk2      = sk*sk;
 681                     sk2_rinv = sk2*rinv;
 682                     prod     = 0.25*sk2_rinv;
 683
 684                     log_term = log(uij*lij_inv);
 685                     /* log_term = table_log(uij*lij_inv,born->log_table,LOG_TABLE_ACCURACY); */
 686                     tmp      = lij-uij + 0.25*dr*diff2 + (0.5*rinv)*log_term + prod*(-diff2);
 687
 688                     if (rai < sk-dr)
 689                     {
 690                         tmp = tmp + 2.0 * (rai_inv-lij);
 691                     }
 692
 693                     t1      = 0.5*lij2 + prod*lij3 - 0.25*(lij*rinv+lij3*dr);
 694                     t2      = -0.5*uij2 - prod*uij3 + 0.25*(uij*rinv+uij3*dr);
 695
 696                     t3      = 0.125*(1.0+sk2_rinv*rinv)*(-diff2)+0.25*log_term*rinv*rinv;
 697
 698                     dadxi = (dlij*t1+t2+t3)*rinv;
 699
 700                     sum_ai += 0.5*tmp;
 701                 }
 702                 else
 703                 {
 704                     dadxi = 0.0;
 705                 }
 706
 707                 /* ai -> aj interaction */
 708                 if (raj < dr + sk_ai)
 709                 {
 710                     lij     = 1.0/(dr-sk_ai);
 711                     dlij    = 1.0;
 712                     raj_inv = 1.0/raj;
 713
 714                     if (raj > dr-sk_ai)
 715                     {
 716                         lij  = raj_inv;
 717                         dlij = 0.0;
 718                     }
 719
 720                     lij2     = lij  * lij;
 721                     lij3     = lij2 * lij;
 722
 723                     uij      = 1.0/(dr+sk_ai);
 724                     uij2     = uij  * uij;
 725                     uij3     = uij2 * uij;
 726
 727                     diff2    = uij2-lij2;
 728
 729                     lij_inv  = gmx_invsqrt(lij2);
 730                     sk2      =  sk2_ai; /* sk2_ai = sk_ai * sk_ai in i loop above */
 731                     sk2_rinv = sk2*rinv;
 732                     prod     = 0.25 * sk2_rinv;
 733
 734                     /* log_term = table_log(uij*lij_inv,born->log_table,LOG_TABLE_ACCURACY); */
 735                     log_term = log(uij*lij_inv);
 736
 737                     tmp      = lij-uij + 0.25*dr*diff2 + (0.5*rinv)*log_term + prod*(-diff2);
 738
 739                     if (raj < sk_ai-dr)
 740                     {
 741                         tmp     = tmp + 2.0 * (raj_inv-lij);
 742                     }
 743
 744                     t1      = 0.5*lij2 + prod*lij3 - 0.25*(lij*rinv+lij3*dr);
 745                     t2      = -0.5*uij2 - 0.25*sk2_rinv*uij3 + 0.25*(uij*rinv+uij3*dr);
 746                     t3      = 0.125*(1.0+sk2_rinv*rinv)*(-diff2)+0.25*log_term*rinv*rinv;
 747
 748                     dadxj = (dlij*t1+t2+t3)*rinv; /* rb2 is moved to chainrule  */
 749
 750                     born->gpol_hct_work[k] += 0.5*tmp;
 751                 }
 752                 else
 753                 {
 754                     dadxj = 0.0;
 755                 }
 756                 fr->dadx[n++] = dadxi;
 757                 fr->dadx[n++] = dadxj;
 758
 759             }
 760         }
 761
 762         /* Main part, no exclusions */
 763         for (j = nj1; j < nj2; j++)
 764         {
 765             k = j%natoms;
 766
 767             /* load j atom coordinates */
 768             jx                = x[3*k];
 769             jy                = x[3*k+1];
 770             jz                = x[3*k+2];
 771
 772             /* Calculate distance */
 773             dx                = ix - jx;
 774             dy                = iy - jy;
 775             dz                = iz - jz;
 776             rsq               = dx*dx+dy*dy+dz*dz;
 777
 778             /* Calculate 1/r and 1/r2 */
 779             rinv              = gmx_invsqrt(rsq);
 780             dr                = rsq*rinv;
 781
 782             /* sk is precalculated in init_gb() */
 783             sk    = born->param[k];
 784             raj   = top->atomtypes.gb_radius[mdatoms->typeA[k]]-doffset;
 785
 786             /* aj -> ai interaction */
 787             if (rai < dr+sk)
 788             {
 789                 lij       = 1.0/(dr-sk);
 790                 dlij      = 1.0;
 791
 792                 if (rai > dr-sk)
 793                 {
 794                     lij  = rai_inv;
 795                     dlij = 0.0;
 796                 }
 797
 798                 uij      = 1.0/(dr+sk);
 799                 lij2     = lij  * lij;
 800                 lij3     = lij2 * lij;
 801                 uij2     = uij  * uij;
 802                 uij3     = uij2 * uij;
 803
 804                 diff2    = uij2-lij2;
 805
 806                 lij_inv  = gmx_invsqrt(lij2);
 807                 sk2      = sk*sk;
 808                 sk2_rinv = sk2*rinv;
 809                 prod     = 0.25*sk2_rinv;
 810
 811                 log_term = log(uij*lij_inv);
 812                 /* log_term = table_log(uij*lij_inv,born->log_table,LOG_TABLE_ACCURACY); */
 813                 tmp      = lij-uij + 0.25*dr*diff2 + (0.5*rinv)*log_term + prod*(-diff2);
 814
 815                 if (rai < sk-dr)
 816                 {
 817                     tmp = tmp + 2.0 * (rai_inv-lij);
 818                 }
 819
 820                 /* duij    = 1.0; */
 821                 t1      = 0.5*lij2 + prod*lij3 - 0.25*(lij*rinv+lij3*dr);
 822                 t2      = -0.5*uij2 - 0.25*sk2_rinv*uij3 + 0.25*(uij*rinv+uij3*dr);
 823                 t3      = 0.125*(1.0+sk2_rinv*rinv)*(-diff2)+0.25*log_term*rinv*rinv;
 824
 825                 dadxi = (dlij*t1+t2+t3)*rinv; /* rb2 is moved to chainrule      */
 826
 827                 sum_ai += 0.5*tmp;
 828             }
 829             else
 830             {
 831                 dadxi = 0.0;
 832             }
 833
 834             /* ai -> aj interaction */
 835             if (raj < dr + sk_ai)
 836             {
 837                 lij     = 1.0/(dr-sk_ai);
 838                 dlij    = 1.0;
 839                 raj_inv = 1.0/raj;
 840
 841                 if (raj > dr-sk_ai)
 842                 {
 843                     lij  = raj_inv;
 844                     dlij = 0.0;
 845                 }
 846
 847                 lij2     = lij  * lij;
 848                 lij3     = lij2 * lij;
 849
 850                 uij      = 1.0/(dr+sk_ai);
 851                 uij2     = uij  * uij;
 852                 uij3     = uij2 * uij;
 853
 854                 diff2    = uij2-lij2;
 855
 856                 lij_inv  = gmx_invsqrt(lij2);
 857                 sk2      =  sk2_ai; /* sk2_ai = sk_ai * sk_ai in i loop above */
 858                 sk2_rinv = sk2*rinv;
 859                 prod     = 0.25 * sk2_rinv;
 860
 861                 /* log_term = table_log(uij*lij_inv,born->log_table,LOG_TABLE_ACCURACY); */
 862                 log_term = log(uij*lij_inv);
 863
 864                 tmp      = lij-uij + 0.25*dr*diff2 + (0.5*rinv)*log_term + prod*(-diff2);
 865
 866                 if (raj < sk_ai-dr)
 867                 {
 868                     tmp     = tmp + 2.0 * (raj_inv-lij);
 869                 }
 870
 871                 t1      = 0.5*lij2 + prod*lij3 - 0.25*(lij*rinv+lij3*dr);
 872                 t2      = -0.5*uij2 - 0.25*sk2_rinv*uij3 + 0.25*(uij*rinv+uij3*dr);
 873                 t3      = 0.125*(1.0+sk2_rinv*rinv)*(-diff2)+0.25*log_term*rinv*rinv;
 874
 875                 dadxj = (dlij*t1+t2+t3)*rinv; /* rb2 is moved to chainrule      */
 876
 877                 born->gpol_hct_work[k] += 0.5*tmp;
 878             }
 879             else
 880             {
 881                 dadxj = 0.0;
 882             }
 883             fr->dadx[n++] = dadxi;
 884             fr->dadx[n++] = dadxj;
 885         }
 886         born->gpol_hct_work[i] += sum_ai;
 887     }
 888
 889     /* Parallel summations */
 890     if (PARTDECOMP(cr))
 891     {
 892         gmx_sum(natoms, born->gpol_hct_work, cr);
 893     }
 894
 895     if (gb_algorithm == egbHCT)
 896     {
 897         /* HCT */
 898         for (i = 0; i < natoms; i++)
 899         {
 900             if (born->use[i] != 0)
 901             {
 902                 rai     = top->atomtypes.gb_radius[mdatoms->typeA[i]]-born->gb_doffset;
 903                 sum_ai  = 1.0/rai - born->gpol_hct_work[i];
 904                 min_rad = rai + born->gb_doffset;
 905                 rad     = 1.0/sum_ai;
 906
 907                 born->bRad[i]   = rad > min_rad ? rad : min_rad;
 908                 fr->invsqrta[i] = gmx_invsqrt(born->bRad[i]);
 909             }
 910         }
 911
 912     }
 913     else
 914     {
 915         /* OBC */
 916         /* Calculate the radii */
 917         for (i = 0; i < natoms; i++)
 918         {
 919             if (born->use[i] != 0)
 920             {
 921                 rai        = top->atomtypes.gb_radius[mdatoms->typeA[i]];
 922                 rai_inv2   = 1.0/rai;
 923                 rai        = rai-doffset;
 924                 rai_inv    = 1.0/rai;
 925                 sum_ai     = rai * born->gpol_hct_work[i];
 926                 sum_ai2    = sum_ai  * sum_ai;
 927                 sum_ai3    = sum_ai2 * sum_ai;
 928
 929                 tsum          = tanh(born->obc_alpha*sum_ai-born->obc_beta*sum_ai2+born->obc_gamma*sum_ai3);
 930                 born->bRad[i] = rai_inv - tsum*rai_inv2;
 931                 born->bRad[i] = 1.0 / born->bRad[i];
 932
 933                 fr->invsqrta[i] = gmx_invsqrt(born->bRad[i]);
 934
 935                 tchain         = rai * (born->obc_alpha-2*born->obc_beta*sum_ai+3*born->obc_gamma*sum_ai2);
 936                 born->drobc[i] = (1.0-tsum*tsum)*tchain*rai_inv2;
 937             }
 938         }
 939     }
 940     return 0;
 941 }
 942
 943
 944
 945
 946
 947 int
 948 genborn_allvsall_calc_chainrule(t_forcerec *           fr,
 949                                 t_mdatoms *            mdatoms,
 950                                 gmx_genborn_t *        born,
 951                                 real *                 x,
 952                                 real *                 f,
 953                                 int                    gb_algorithm,
 954                                 void *                 work)
 955 {
 956     gmx_allvsallgb2_data_t *aadata;
 957     int                     natoms;
 958     int                     ni0, ni1;
 959     int                     nj0, nj1, nj2;
 960     int                     i, j, k, n;
 961     int                     idx;
 962     int              *      mask;
 963
 964     real                    ix, iy, iz;
 965     real                    fix, fiy, fiz;
 966     real                    jx, jy, jz;
 967     real                    dx, dy, dz;
 968     real                    tx, ty, tz;
 969     real                    rbai, rbaj, fgb, fgb_ai, rbi;
 970     real              *     rb;
 971     real              *     dadx;
 972
 973     natoms              = mdatoms->nr;
 974     ni0                 = mdatoms->start;
 975     ni1                 = mdatoms->start+mdatoms->homenr;
 976     dadx                = fr->dadx;
 977
 978     aadata = (gmx_allvsallgb2_data_t *)work;
 979
 980     n  = 0;
 981     rb = born->work;
 982
 983     /* Loop to get the proper form for the Born radius term */
 984     if (gb_algorithm == egbSTILL)
 985     {
 986         for (i = 0; i < natoms; i++)
 987         {
 988             rbi   = born->bRad[i];
 989             rb[i] = (2 * rbi * rbi * fr->dvda[i])/ONE_4PI_EPS0;
 990         }
 991     }
 992     else if (gb_algorithm == egbHCT)
 993     {
 994         for (i = 0; i < natoms; i++)
 995         {
 996             rbi   = born->bRad[i];
 997             rb[i] = rbi * rbi * fr->dvda[i];
 998         }
 999     }
1000     else if (gb_algorithm == egbOBC)
1001     {
1002         for (idx = 0; idx < natoms; idx++)
1003         {
1004             rbi     = born->bRad[idx];
1005             rb[idx] = rbi * rbi * born->drobc[idx] * fr->dvda[idx];
1006         }
1007     }
1008
1009     for (i = ni0; i < ni1; i++)
1010     {
1011         /* We assume shifts are NOT used for all-vs-all interactions */
1012
1013         /* Load i atom data */
1014         ix                = x[3*i];
1015         iy                = x[3*i+1];
1016         iz                = x[3*i+2];
1017
1018         fix               = 0;
1019         fiy               = 0;
1020         fiz               = 0;
1021
1022         rbai              = rb[i];
1023
1024         /* Load limits for loop over neighbors */
1025         nj0              = aadata->jindex_gb[3*i];
1026         nj1              = aadata->jindex_gb[3*i+1];
1027         nj2              = aadata->jindex_gb[3*i+2];
1028
1029         mask             = aadata->exclusion_mask_gb[i];
1030
1031         /* Prologue part, including exclusion mask */
1032         for (j = nj0; j < nj1; j++, mask++)
1033         {
1034             if (*mask != 0)
1035             {
1036                 k = j%natoms;
1037
1038                 /* load j atom coordinates */
1039                 jx                = x[3*k];
1040                 jy                = x[3*k+1];
1041                 jz                = x[3*k+2];
1042
1043                 /* Calculate distance */
1044                 dx                = ix - jx;
1045                 dy                = iy - jy;
1046                 dz                = iz - jz;
1047
1048                 rbaj              = rb[k];
1049
1050                 fgb     = rbai*dadx[n++];
1051                 fgb_ai  = rbaj*dadx[n++];
1052
1053                 /* Total force between ai and aj is the sum of ai->aj and aj->ai */
1054                 fgb     = fgb + fgb_ai;
1055
1056                 tx      = fgb * dx;
1057                 ty      = fgb * dy;
1058                 tz      = fgb * dz;
1059
1060                 fix     = fix + tx;
1061                 fiy     = fiy + ty;
1062                 fiz     = fiz + tz;
1063
1064                 /* Update force on atom aj */
1065                 f[3*k]   = f[3*k] - tx;
1066                 f[3*k+1] = f[3*k+1] - ty;
1067                 f[3*k+2] = f[3*k+2] - tz;
1068             }
1069         }
1070
1071         /* Main part, no exclusions */
1072         for (j = nj1; j < nj2; j++)
1073         {
1074             k = j%natoms;
1075
1076             /* load j atom coordinates */
1077             jx                = x[3*k];
1078             jy                = x[3*k+1];
1079             jz                = x[3*k+2];
1080
1081             /* Calculate distance */
1082             dx                = ix - jx;
1083             dy                = iy - jy;
1084             dz                = iz - jz;
1085
1086             rbaj              = rb[k];
1087
1088             fgb     = rbai*dadx[n++];
1089             fgb_ai  = rbaj*dadx[n++];
1090
1091             /* Total force between ai and aj is the sum of ai->aj and aj->ai */
1092             fgb     = fgb + fgb_ai;
1093
1094             tx      = fgb * dx;
1095             ty      = fgb * dy;
1096             tz      = fgb * dz;
1097
1098             fix     = fix + tx;
1099             fiy     = fiy + ty;
1100             fiz     = fiz + tz;
1101
1102             /* Update force on atom aj */
1103             f[3*k]   = f[3*k] - tx;
1104             f[3*k+1] = f[3*k+1] - ty;
1105             f[3*k+2] = f[3*k+2] - tz;
1106         }
1107         /* Update force and shift forces on atom ai */
1108         f[3*i]   = f[3*i] + fix;
1109         f[3*i+1] = f[3*i+1] + fiy;
1110         f[3*i+2] = f[3*i+2] + fiz;
1111     }
1112
1113     return 0;
1114 }