src/gromacs/mdlib/genborn_allvsall_sse2_double.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 #if 0 && defined (GMX_X86_SSE2)
  48
  49 #include <gmx_sse2_double.h>
  50
  51
  52 #define SIMD_WIDTH 2
  53 #define UNROLLI    2
  54 #define UNROLLJ    2
  55
  56
  57
  58
  59
  60
  61
  62
  63
  64 typedef struct
  65 {
  66     int   *      jindex_gb;
  67     int   **     prologue_mask_gb;
  68     int   **     epilogue_mask;
  69     int   *      imask;
  70     double *     gb_radius;
  71     double *     workparam;
  72     double *     work;
  73     double *     x_align;
  74     double *     y_align;
  75     double *     z_align;
  76     double *     fx_align;
  77     double *     fy_align;
  78     double *     fz_align;
  79 }
  80 gmx_allvsallgb2_data_t;
  81
  82
  83 static int
  84 calc_maxoffset(int i, int natoms)
  85 {
  86     int maxoffset;
  87
  88     if ((natoms % 2) == 1)
  89     {
  90         /* Odd number of atoms, easy */
  91         maxoffset = natoms/2;
  92     }
  93     else if ((natoms % 4) == 0)
  94     {
  95         /* Multiple of four is hard */
  96         if (i < natoms/2)
  97         {
  98             if ((i % 2) == 0)
  99             {
 100                 maxoffset = natoms/2;
 101             }
 102             else
 103             {
 104                 maxoffset = natoms/2-1;
 105             }
 106         }
 107         else
 108         {
 109             if ((i % 2) == 1)
 110             {
 111                 maxoffset = natoms/2;
 112             }
 113             else
 114             {
 115                 maxoffset = natoms/2-1;
 116             }
 117         }
 118     }
 119     else
 120     {
 121         /* natoms/2 = odd */
 122         if ((i % 2) == 0)
 123         {
 124             maxoffset = natoms/2;
 125         }
 126         else
 127         {
 128             maxoffset = natoms/2-1;
 129         }
 130     }
 131
 132     return maxoffset;
 133 }
 134
 135 static void
 136 setup_gb_exclusions_and_indices(gmx_allvsallgb2_data_t     *   aadata,
 137                                 t_ilist     *                  ilist,
 138                                 int                            start,
 139                                 int                            end,
 140                                 int                            natoms,
 141                                 gmx_bool                       bInclude12,
 142                                 gmx_bool                       bInclude13,
 143                                 gmx_bool                       bInclude14)
 144 {
 145     int   i, j, k, tp;
 146     int   a1, a2;
 147     int   ni0, ni1, nj0, nj1, nj;
 148     int   imin, imax, iexcl;
 149     int   max_offset;
 150     int   max_excl_offset;
 151     int   firstinteraction;
 152     int   ibase;
 153     int  *pi;
 154
 155     /* This routine can appear to be a bit complex, but it is mostly book-keeping.
 156      * To enable the fast all-vs-all kernel we need to be able to stream through all coordinates
 157      * whether they should interact or not.
 158      *
 159      * To avoid looping over the exclusions, we create a simple mask that is 1 if the interaction
 160      * should be present, otherwise 0. Since exclusions typically only occur when i & j are close,
 161      * we create a jindex array with three elements per i atom: the starting point, the point to
 162      * which we need to check exclusions, and the end point.
 163      * This way we only have to allocate a short exclusion mask per i atom.
 164      */
 165
 166     ni0 = (start/UNROLLI)*UNROLLI;
 167     ni1 = ((end+UNROLLI-1)/UNROLLI)*UNROLLI;
 168
 169     /* Set the interaction mask to only enable the i atoms we want to include */
 170     snew(pi, 2*(natoms+UNROLLI+2*SIMD_WIDTH));
 171     aadata->imask = (int *) (((size_t) pi + 16) & (~((size_t) 15)));
 172     for (i = 0; i < natoms+UNROLLI; i++)
 173     {
 174         aadata->imask[2*i]   = (i >= start && i < end) ? 0xFFFFFFFF : 0;
 175         aadata->imask[2*i+1] = (i >= start && i < end) ? 0xFFFFFFFF : 0;
 176     }
 177
 178     /* Allocate memory for our modified jindex array */
 179     snew(aadata->jindex_gb, 4*(natoms+UNROLLI));
 180     for (i = 0; i < 4*(natoms+UNROLLI); i++)
 181     {
 182         aadata->jindex_gb[i] = 0;
 183     }
 184
 185     /* Create the exclusion masks for the prologue part */
 186     snew(aadata->prologue_mask_gb, natoms+UNROLLI); /* list of pointers */
 187
 188     /* First zero everything to avoid uninitialized data */
 189     for (i = 0; i < natoms+UNROLLI; i++)
 190     {
 191         aadata->prologue_mask_gb[i] = NULL;
 192     }
 193
 194     /* Calculate the largest exclusion range we need for each UNROLLI-tuplet of i atoms. */
 195     for (ibase = ni0; ibase < ni1; ibase += UNROLLI)
 196     {
 197         max_excl_offset = -1;
 198
 199         /* First find maxoffset for the next 4 atoms (or fewer if we are close to end) */
 200         imax = ((ibase+UNROLLI) < end) ? (ibase+UNROLLI) : end;
 201
 202         /* Which atom is the first we (might) interact with? */
 203         imin = natoms; /* Guaranteed to be overwritten by one of 'firstinteraction' */
 204         for (i = ibase; i < imax; i++)
 205         {
 206             /* Before exclusions, which atom is the first we (might) interact with? */
 207             firstinteraction = i+1;
 208             max_offset       = calc_maxoffset(i, natoms);
 209
 210             if (!bInclude12)
 211             {
 212                 for (j = 0; j < ilist[F_GB12].nr; j += 3)
 213                 {
 214                     a1 = ilist[F_GB12].iatoms[j+1];
 215                     a2 = ilist[F_GB12].iatoms[j+2];
 216
 217                     if (a1 == i)
 218                     {
 219                         k = a2;
 220                     }
 221                     else if (a2 == i)
 222                     {
 223                         k = a1;
 224                     }
 225                     else
 226                     {
 227                         continue;
 228                     }
 229
 230                     if (k == firstinteraction)
 231                     {
 232                         firstinteraction++;
 233                     }
 234                 }
 235             }
 236             if (!bInclude13)
 237             {
 238                 for (j = 0; j < ilist[F_GB13].nr; j += 3)
 239                 {
 240                     a1 = ilist[F_GB13].iatoms[j+1];
 241                     a2 = ilist[F_GB13].iatoms[j+2];
 242
 243                     if (a1 == i)
 244                     {
 245                         k = a2;
 246                     }
 247                     else if (a2 == i)
 248                     {
 249                         k = a1;
 250                     }
 251                     else
 252                     {
 253                         continue;
 254                     }
 255
 256                     if (k == firstinteraction)
 257                     {
 258                         firstinteraction++;
 259                     }
 260                 }
 261             }
 262             if (!bInclude14)
 263             {
 264                 for (j = 0; j < ilist[F_GB14].nr; j += 3)
 265                 {
 266                     a1 = ilist[F_GB14].iatoms[j+1];
 267                     a2 = ilist[F_GB14].iatoms[j+2];
 268                     if (a1 == i)
 269                     {
 270                         k = a2;
 271                     }
 272                     else if (a2 == i)
 273                     {
 274                         k = a1;
 275                     }
 276                     else
 277                     {
 278                         continue;
 279                     }
 280
 281                     if (k == firstinteraction)
 282                     {
 283                         firstinteraction++;
 284                     }
 285                 }
 286             }
 287             imin = (firstinteraction < imin) ? firstinteraction : imin;
 288         }
 289         /* round down to j unrolling factor */
 290         imin = (imin/UNROLLJ)*UNROLLJ;
 291
 292         for (i = ibase; i < imax; i++)
 293         {
 294             max_offset = calc_maxoffset(i, natoms);
 295
 296             if (!bInclude12)
 297             {
 298                 for (j = 0; j < ilist[F_GB12].nr; j += 3)
 299                 {
 300                     a1 = ilist[F_GB12].iatoms[j+1];
 301                     a2 = ilist[F_GB12].iatoms[j+2];
 302
 303                     if (a1 == i)
 304                     {
 305                         k = a2;
 306                     }
 307                     else if (a2 == i)
 308                     {
 309                         k = a1;
 310                     }
 311                     else
 312                     {
 313                         continue;
 314                     }
 315
 316                     if (k < imin)
 317                     {
 318                         k += natoms;
 319                     }
 320
 321                     if (k > i+max_offset)
 322                     {
 323                         continue;
 324                     }
 325
 326                     k = k - imin;
 327
 328                     if (k+natoms <= max_offset)
 329                     {
 330                         k += natoms;
 331                     }
 332                     max_excl_offset = (k > max_excl_offset) ? k : max_excl_offset;
 333                 }
 334             }
 335             if (!bInclude13)
 336             {
 337                 for (j = 0; j < ilist[F_GB13].nr; j += 3)
 338                 {
 339                     a1 = ilist[F_GB13].iatoms[j+1];
 340                     a2 = ilist[F_GB13].iatoms[j+2];
 341
 342                     if (a1 == i)
 343                     {
 344                         k = a2;
 345                     }
 346                     else if (a2 == i)
 347                     {
 348                         k = a1;
 349                     }
 350                     else
 351                     {
 352                         continue;
 353                     }
 354
 355                     if (k < imin)
 356                     {
 357                         k += natoms;
 358                     }
 359
 360                     if (k > i+max_offset)
 361                     {
 362                         continue;
 363                     }
 364
 365                     k = k - imin;
 366
 367                     if (k+natoms <= max_offset)
 368                     {
 369                         k += natoms;
 370                     }
 371                     max_excl_offset = (k > max_excl_offset) ? k : max_excl_offset;
 372                 }
 373             }
 374             if (!bInclude14)
 375             {
 376                 for (j = 0; j < ilist[F_GB14].nr; j += 3)
 377                 {
 378                     a1 = ilist[F_GB14].iatoms[j+1];
 379                     a2 = ilist[F_GB14].iatoms[j+2];
 380
 381                     if (a1 == i)
 382                     {
 383                         k = a2;
 384                     }
 385                     else if (a2 == i)
 386                     {
 387                         k = a1;
 388                     }
 389                     else
 390                     {
 391                         continue;
 392                     }
 393
 394                     if (k < imin)
 395                     {
 396                         k += natoms;
 397                     }
 398
 399                     if (k > i+max_offset)
 400                     {
 401                         continue;
 402                     }
 403
 404                     k = k - imin;
 405
 406                     if (k+natoms <= max_offset)
 407                     {
 408                         k += natoms;
 409                     }
 410                     max_excl_offset = (k > max_excl_offset) ? k : max_excl_offset;
 411                 }
 412             }
 413         }
 414
 415         /* The offset specifies the last atom to be excluded, so add one unit to get an upper loop limit */
 416         max_excl_offset++;
 417         /* round up to j unrolling factor */
 418         max_excl_offset = (max_excl_offset/UNROLLJ+1)*UNROLLJ;
 419
 420         /* Set all the prologue masks length to this value (even for i>end) */
 421         for (i = ibase; i < ibase+UNROLLI; i++)
 422         {
 423             aadata->jindex_gb[4*i]   = imin;
 424             aadata->jindex_gb[4*i+1] = imin+max_excl_offset;
 425         }
 426     }
 427
 428     /* Now the hard part, loop over it all again to calculate the actual contents of the prologue masks */
 429     for (ibase = ni0; ibase < ni1; ibase += UNROLLI)
 430     {
 431         for (i = ibase; i < ibase+UNROLLI; i++)
 432         {
 433             nj   = aadata->jindex_gb[4*i+1] - aadata->jindex_gb[4*i];
 434             imin = aadata->jindex_gb[4*i];
 435
 436             /* Allocate aligned memory */
 437             snew(pi, 2*(nj+2*SIMD_WIDTH));
 438             aadata->prologue_mask_gb[i] = (int *) (((size_t) pi + 16) & (~((size_t) 15)));
 439
 440             max_offset = calc_maxoffset(i, natoms);
 441
 442             /* Include interactions i+1 <= j < i+maxoffset */
 443             for (k = 0; k < nj; k++)
 444             {
 445                 j = imin + k;
 446
 447                 if ( (j > i) && (j <= i+max_offset) )
 448                 {
 449                     aadata->prologue_mask_gb[i][2*k]   = 0xFFFFFFFF;
 450                     aadata->prologue_mask_gb[i][2*k+1] = 0xFFFFFFFF;
 451                 }
 452                 else
 453                 {
 454                     aadata->prologue_mask_gb[i][2*k]   = 0;
 455                     aadata->prologue_mask_gb[i][2*k+1] = 0;
 456                 }
 457             }
 458
 459             /* Clear out the explicit exclusions */
 460             if (i < end)
 461             {
 462                 if (!bInclude12)
 463                 {
 464                     for (j = 0; j < ilist[F_GB12].nr; j += 3)
 465                     {
 466                         a1 = ilist[F_GB12].iatoms[j+1];
 467                         a2 = ilist[F_GB12].iatoms[j+2];
 468
 469                         if (a1 == i)
 470                         {
 471                             k = a2;
 472                         }
 473                         else if (a2 == i)
 474                         {
 475                             k = a1;
 476                         }
 477                         else
 478                         {
 479                             continue;
 480                         }
 481
 482                         if (k > i+max_offset)
 483                         {
 484                             continue;
 485                         }
 486                         k = k-i;
 487
 488                         if (k+natoms <= max_offset)
 489                         {
 490                             k += natoms;
 491                         }
 492
 493                         k = k+i-imin;
 494                         if (k >= 0)
 495                         {
 496                             aadata->prologue_mask_gb[i][2*k]   = 0;
 497                             aadata->prologue_mask_gb[i][2*k+1] = 0;
 498                         }
 499                     }
 500                 }
 501                 if (!bInclude13)
 502                 {
 503                     for (j = 0; j < ilist[F_GB13].nr; j += 3)
 504                     {
 505                         a1 = ilist[F_GB13].iatoms[j+1];
 506                         a2 = ilist[F_GB13].iatoms[j+2];
 507
 508                         if (a1 == i)
 509                         {
 510                             k = a2;
 511                         }
 512                         else if (a2 == i)
 513                         {
 514                             k = a1;
 515                         }
 516                         else
 517                         {
 518                             continue;
 519                         }
 520
 521                         if (k > i+max_offset)
 522                         {
 523                             continue;
 524                         }
 525                         k = k-i;
 526
 527                         if (k+natoms <= max_offset)
 528                         {
 529                             k += natoms;
 530                         }
 531
 532                         k = k+i-imin;
 533                         if (k >= 0)
 534                         {
 535                             aadata->prologue_mask_gb[i][2*k]   = 0;
 536                             aadata->prologue_mask_gb[i][2*k+1] = 0;
 537                         }
 538                     }
 539                 }
 540                 if (!bInclude14)
 541                 {
 542                     for (j = 0; j < ilist[F_GB14].nr; j += 3)
 543                     {
 544                         a1 = ilist[F_GB14].iatoms[j+1];
 545                         a2 = ilist[F_GB14].iatoms[j+2];
 546
 547                         if (a1 == i)
 548                         {
 549                             k = a2;
 550                         }
 551                         else if (a2 == i)
 552                         {
 553                             k = a1;
 554                         }
 555                         else
 556                         {
 557                             continue;
 558                         }
 559
 560                         if (k > i+max_offset)
 561                         {
 562                             continue;
 563                         }
 564                         k = k-i;
 565
 566                         if (k+natoms <= max_offset)
 567                         {
 568                             k += natoms;
 569                         }
 570
 571                         k = k+i-imin;
 572                         if (k >= 0)
 573                         {
 574                             aadata->prologue_mask_gb[i][2*k]   = 0;
 575                             aadata->prologue_mask_gb[i][2*k+1] = 0;
 576                         }
 577                     }
 578                 }
 579             }
 580         }
 581     }
 582
 583     /* Construct the epilogue mask - this just contains the check for maxoffset */
 584     snew(aadata->epilogue_mask, natoms+UNROLLI);
 585
 586     /* First zero everything to avoid uninitialized data */
 587     for (i = 0; i < natoms+UNROLLI; i++)
 588     {
 589         aadata->jindex_gb[4*i+2]    = aadata->jindex_gb[4*i+1];
 590         aadata->jindex_gb[4*i+3]    = aadata->jindex_gb[4*i+1];
 591         aadata->epilogue_mask[i]    = NULL;
 592     }
 593
 594     for (ibase = ni0; ibase < ni1; ibase += UNROLLI)
 595     {
 596         /* Find the lowest index for which we need to use the epilogue */
 597         imin       = ibase;
 598         max_offset = calc_maxoffset(imin, natoms);
 599
 600         imin = imin + 1 + max_offset;
 601
 602         /* Find largest index for which we need to use the epilogue */
 603         imax = ibase + UNROLLI-1;
 604         imax = (imax < end) ? imax : end;
 605
 606         max_offset = calc_maxoffset(imax, natoms);
 607         imax       = imax + 1 + max_offset + UNROLLJ - 1;
 608
 609         for (i = ibase; i < ibase+UNROLLI; i++)
 610         {
 611             /* Start of epilogue - round down to j tile limit */
 612             aadata->jindex_gb[4*i+2] = (imin/UNROLLJ)*UNROLLJ;
 613             /* Make sure we dont overlap - for small systems everything is done in the prologue */
 614             aadata->jindex_gb[4*i+2] = (aadata->jindex_gb[4*i+1] > aadata->jindex_gb[4*i+2]) ? aadata->jindex_gb[4*i+1] : aadata->jindex_gb[4*i+2];
 615             /* Round upwards to j tile limit */
 616             aadata->jindex_gb[4*i+3] = (imax/UNROLLJ)*UNROLLJ;
 617             /* Make sure we dont have a negative range for the epilogue */
 618             aadata->jindex_gb[4*i+3] = (aadata->jindex_gb[4*i+2] > aadata->jindex_gb[4*i+3]) ? aadata->jindex_gb[4*i+2] : aadata->jindex_gb[4*i+3];
 619         }
 620     }
 621
 622     /* And fill it with data... */
 623
 624     for (ibase = ni0; ibase < ni1; ibase += UNROLLI)
 625     {
 626         for (i = ibase; i < ibase+UNROLLI; i++)
 627         {
 628
 629             nj = aadata->jindex_gb[4*i+3] - aadata->jindex_gb[4*i+2];
 630
 631             /* Allocate aligned memory */
 632             snew(pi, 2*(nj+2*SIMD_WIDTH));
 633             aadata->epilogue_mask[i] = (int *) (((size_t) pi + 16) & (~((size_t) 15)));
 634
 635             max_offset = calc_maxoffset(i, natoms);
 636
 637             for (k = 0; k < nj; k++)
 638             {
 639                 j = aadata->jindex_gb[4*i+2] + k;
 640                 aadata->epilogue_mask[i][2*k]   = (j <= i+max_offset) ? 0xFFFFFFFF : 0;
 641                 aadata->epilogue_mask[i][2*k+1] = (j <= i+max_offset) ? 0xFFFFFFFF : 0;
 642             }
 643         }
 644     }
 645 }
 646
 647
 648 static void
 649 genborn_allvsall_setup(gmx_allvsallgb2_data_t     **  p_aadata,
 650                        gmx_localtop_t     *           top,
 651                        gmx_genborn_t     *            born,
 652                        t_mdatoms     *                mdatoms,
 653                        double                         radius_offset,
 654                        int                            gb_algorithm,
 655                        gmx_bool                       bInclude12,
 656                        gmx_bool                       bInclude13,
 657                        gmx_bool                       bInclude14)
 658 {
 659     int                     i, j, idx;
 660     int                     natoms;
 661     gmx_allvsallgb2_data_t *aadata;
 662     double                 *p;
 663
 664     natoms = mdatoms->nr;
 665
 666     snew(aadata, 1);
 667     *p_aadata = aadata;
 668
 669     snew(p, 2*natoms+2*SIMD_WIDTH);
 670     aadata->x_align = (double *) (((size_t) p + 16) & (~((size_t) 15)));
 671     snew(p, 2*natoms+2*SIMD_WIDTH);
 672     aadata->y_align = (double *) (((size_t) p + 16) & (~((size_t) 15)));
 673     snew(p, 2*natoms+2*SIMD_WIDTH);
 674     aadata->z_align = (double *) (((size_t) p + 16) & (~((size_t) 15)));
 675     snew(p, 2*natoms+2*SIMD_WIDTH);
 676     aadata->fx_align = (double *) (((size_t) p + 16) & (~((size_t) 15)));
 677     snew(p, 2*natoms+2*SIMD_WIDTH);
 678     aadata->fy_align = (double *) (((size_t) p + 16) & (~((size_t) 15)));
 679     snew(p, 2*natoms+2*SIMD_WIDTH);
 680     aadata->fz_align = (double *) (((size_t) p + 16) & (~((size_t) 15)));
 681
 682     snew(p, 2*natoms+UNROLLJ+SIMD_WIDTH);
 683     aadata->gb_radius = (double *) (((size_t) p + 16) & (~((size_t) 15)));
 684
 685     snew(p, 2*natoms+UNROLLJ+SIMD_WIDTH);
 686     aadata->workparam = (double *) (((size_t) p + 16) & (~((size_t) 15)));
 687
 688     snew(p, 2*natoms+UNROLLJ+SIMD_WIDTH);
 689     aadata->work = (double *) (((size_t) p + 16) & (~((size_t) 15)));
 690
 691     for (i = 0; i < mdatoms->nr; i++)
 692     {
 693         aadata->gb_radius[i] = top->atomtypes.gb_radius[mdatoms->typeA[i]] - radius_offset;
 694         if (gb_algorithm == egbSTILL)
 695         {
 696             aadata->workparam[i] = born->vsolv[i];
 697         }
 698         else if (gb_algorithm == egbOBC)
 699         {
 700             aadata->workparam[i] = born->param[i];
 701         }
 702         aadata->work[i]      = 0.0;
 703     }
 704     for (i = 0; i < mdatoms->nr; i++)
 705     {
 706         aadata->gb_radius[natoms+i] = aadata->gb_radius[i];
 707         aadata->workparam[natoms+i] = aadata->workparam[i];
 708         aadata->work[natoms+i]      = aadata->work[i];
 709     }
 710
 711     for (i = 0; i < 2*natoms+SIMD_WIDTH; i++)
 712     {
 713         aadata->x_align[i]  = 0.0;
 714         aadata->y_align[i]  = 0.0;
 715         aadata->z_align[i]  = 0.0;
 716         aadata->fx_align[i] = 0.0;
 717         aadata->fy_align[i] = 0.0;
 718         aadata->fz_align[i] = 0.0;
 719     }
 720
 721     setup_gb_exclusions_and_indices(aadata, top->idef.il, mdatoms->start, mdatoms->start+mdatoms->homenr, mdatoms->nr,
 722                                     bInclude12, bInclude13, bInclude14);
 723 }
 724
 725
 726 /*
 727  * This routine apparently hits a compiler bug visual studio has had 'forever'.
 728  * It is present both in VS2005 and VS2008, and the only way around it is to
 729  * decrease optimization. We do that with at pragma, and only for MSVC, so it
 730  * will not hurt any of the well-behaving and supported compilers out there.
 731  * MS: Fix your compiler, it sucks like a black hole!
 732  */
 733 #ifdef _MSC_VER
 734 #pragma optimize("t",off)
 735 #endif
 736
 737 int
 738 genborn_allvsall_calc_still_radii_sse2_double(t_forcerec   *           fr,
 739                                               t_mdatoms   *            mdatoms,
 740                                               gmx_genborn_t   *        born,
 741                                               gmx_localtop_t   *       top,
 742                                               double *                 x,
 743                                               t_commrec   *            cr,
 744                                               void   *                 paadata)
 745 {
 746     gmx_allvsallgb2_data_t *aadata;
 747     int                     natoms;
 748     int                     ni0, ni1;
 749     int                     nj0, nj1, nj2, nj3;
 750     int                     i, j, k, n;
 751     int              *      mask;
 752     int              *      pmask0;
 753     int              *      pmask1;
 754     int              *      emask0;
 755     int              *      emask1;
 756     double                  ix, iy, iz;
 757     double                  jx, jy, jz;
 758     double                  dx, dy, dz;
 759     double                  rsq, rinv;
 760     double                  gpi, rai, vai;
 761     double                  prod_ai;
 762     double                  irsq, idr4, idr6;
 763     double                  raj, rvdw, ratio;
 764     double                  vaj, ccf, dccf, theta, cosq;
 765     double                  term, prod, icf4, icf6, gpi2, factor, sinq;
 766     double            *     gb_radius;
 767     double            *     vsolv;
 768     double            *     work;
 769     double                  tmpsum[2];
 770     double            *     x_align;
 771     double            *     y_align;
 772     double            *     z_align;
 773     int              *      jindex;
 774     double            *     dadx;
 775
 776     __m128d                 ix_SSE0, iy_SSE0, iz_SSE0;
 777     __m128d                 ix_SSE1, iy_SSE1, iz_SSE1;
 778     __m128d                 gpi_SSE0, rai_SSE0, prod_ai_SSE0;
 779     __m128d                 gpi_SSE1, rai_SSE1, prod_ai_SSE1;
 780     __m128d                 imask_SSE0, jmask_SSE0;
 781     __m128d                 imask_SSE1, jmask_SSE1;
 782     __m128d                 jx_SSE, jy_SSE, jz_SSE;
 783     __m128d                 dx_SSE0, dy_SSE0, dz_SSE0;
 784     __m128d                 dx_SSE1, dy_SSE1, dz_SSE1;
 785     __m128d                 rsq_SSE0, rinv_SSE0, irsq_SSE0, idr4_SSE0, idr6_SSE0;
 786     __m128d                 rsq_SSE1, rinv_SSE1, irsq_SSE1, idr4_SSE1, idr6_SSE1;
 787     __m128d                 raj_SSE, vaj_SSE, prod_SSE;
 788     __m128d                 rvdw_SSE0, ratio_SSE0;
 789     __m128d                 rvdw_SSE1, ratio_SSE1;
 790     __m128d                 theta_SSE0, sinq_SSE0, cosq_SSE0, term_SSE0;
 791     __m128d                 theta_SSE1, sinq_SSE1, cosq_SSE1, term_SSE1;
 792     __m128d                 ccf_SSE0, dccf_SSE0;
 793     __m128d                 ccf_SSE1, dccf_SSE1;
 794     __m128d                 icf4_SSE0, icf6_SSE0;
 795     __m128d                 icf4_SSE1, icf6_SSE1;
 796     __m128d                 half_SSE, one_SSE, two_SSE, four_SSE;
 797     __m128d                 still_p4_SSE, still_p5inv_SSE, still_pip5_SSE;
 798
 799     natoms              = mdatoms->nr;
 800     ni0                 = (mdatoms->start/SIMD_WIDTH)*SIMD_WIDTH;
 801     ni1                 = mdatoms->start+mdatoms->homenr;
 802
 803     n = 0;
 804
 805     aadata = *((gmx_allvsallgb2_data_t **)paadata);
 806
 807
 808     if (aadata == NULL)
 809     {
 810         genborn_allvsall_setup(&aadata, top, born, mdatoms, 0.0,
 811                                egbSTILL, FALSE, FALSE, TRUE);
 812         *((gmx_allvsallgb2_data_t **)paadata) = aadata;
 813     }
 814
 815     x_align = aadata->x_align;
 816     y_align = aadata->y_align;
 817     z_align = aadata->z_align;
 818
 819     gb_radius = aadata->gb_radius;
 820     vsolv     = aadata->workparam;
 821     work      = aadata->work;
 822     jindex    = aadata->jindex_gb;
 823     dadx      = fr->dadx;
 824
 825     still_p4_SSE    = _mm_set1_pd(STILL_P4);
 826     still_p5inv_SSE = _mm_set1_pd(STILL_P5INV);
 827     still_pip5_SSE  = _mm_set1_pd(STILL_PIP5);
 828     half_SSE        = _mm_set1_pd(0.5);
 829     one_SSE         = _mm_set1_pd(1.0);
 830     two_SSE         = _mm_set1_pd(2.0);
 831     four_SSE        = _mm_set1_pd(4.0);
 832
 833     /* This will be summed, so it has to extend to natoms + buffer */
 834     for (i = 0; i < natoms+1+natoms/2; i++)
 835     {
 836         work[i] = 0;
 837     }
 838
 839     for (i = ni0; i < ni1+1+natoms/2; i++)
 840     {
 841         k           = i%natoms;
 842         x_align[i]  = x[3*k];
 843         y_align[i]  = x[3*k+1];
 844         z_align[i]  = x[3*k+2];
 845         work[i]     = 0;
 846     }
 847
 848     for (i = ni0; i < ni1; i += UNROLLI)
 849     {
 850         /* We assume shifts are NOT used for all-vs-all interactions */
 851         /* Load i atom data */
 852         ix_SSE0          = _mm_load1_pd(x_align+i);
 853         iy_SSE0          = _mm_load1_pd(y_align+i);
 854         iz_SSE0          = _mm_load1_pd(z_align+i);
 855         ix_SSE1          = _mm_load1_pd(x_align+i+1);
 856         iy_SSE1          = _mm_load1_pd(y_align+i+1);
 857         iz_SSE1          = _mm_load1_pd(z_align+i+1);
 858
 859         gpi_SSE0         = _mm_setzero_pd();
 860         gpi_SSE1         = _mm_setzero_pd();
 861
 862         rai_SSE0         = _mm_load1_pd(gb_radius+i);
 863         rai_SSE1         = _mm_load1_pd(gb_radius+i+1);
 864
 865         prod_ai_SSE0     = _mm_set1_pd(STILL_P4*vsolv[i]);
 866         prod_ai_SSE1     = _mm_set1_pd(STILL_P4*vsolv[i+1]);
 867
 868         /* Load limits for loop over neighbors */
 869         nj0              = jindex[4*i];
 870         nj1              = jindex[4*i+1];
 871         nj2              = jindex[4*i+2];
 872         nj3              = jindex[4*i+3];
 873
 874         pmask0           = aadata->prologue_mask_gb[i];
 875         pmask1           = aadata->prologue_mask_gb[i+1];
 876         emask0           = aadata->epilogue_mask[i];
 877         emask1           = aadata->epilogue_mask[i+1];
 878
 879         imask_SSE0        = _mm_load1_pd((double *)(aadata->imask+2*i));
 880         imask_SSE1        = _mm_load1_pd((double *)(aadata->imask+2*i+2));
 881
 882         /* Prologue part, including exclusion mask */
 883         for (j = nj0; j < nj1; j += UNROLLJ)
 884         {
 885             jmask_SSE0 = _mm_load_pd((double *)pmask0);
 886             jmask_SSE1 = _mm_load_pd((double *)pmask1);
 887             pmask0    += 2*UNROLLJ;
 888             pmask1    += 2*UNROLLJ;
 889
 890             /* load j atom coordinates */
 891             jx_SSE            = _mm_load_pd(x_align+j);
 892             jy_SSE            = _mm_load_pd(y_align+j);
 893             jz_SSE            = _mm_load_pd(z_align+j);
 894
 895             /* Calculate distance */
 896             dx_SSE0            = _mm_sub_pd(ix_SSE0, jx_SSE);
 897             dy_SSE0            = _mm_sub_pd(iy_SSE0, jy_SSE);
 898             dz_SSE0            = _mm_sub_pd(iz_SSE0, jz_SSE);
 899             dx_SSE1            = _mm_sub_pd(ix_SSE1, jx_SSE);
 900             dy_SSE1            = _mm_sub_pd(iy_SSE1, jy_SSE);
 901             dz_SSE1            = _mm_sub_pd(iz_SSE1, jz_SSE);
 902
 903             /* rsq = dx*dx+dy*dy+dz*dz */
 904             rsq_SSE0           = gmx_mm_calc_rsq_pd(dx_SSE0, dy_SSE0, dz_SSE0);
 905             rsq_SSE1           = gmx_mm_calc_rsq_pd(dx_SSE1, dy_SSE1, dz_SSE1);
 906
 907             /* Combine masks */
 908             jmask_SSE0         = _mm_and_pd(jmask_SSE0, imask_SSE0);
 909             jmask_SSE1         = _mm_and_pd(jmask_SSE1, imask_SSE1);
 910
 911             /* Calculate 1/r and 1/r2 */
 912             rinv_SSE0          = gmx_mm_invsqrt_pd(rsq_SSE0);
 913             rinv_SSE1          = gmx_mm_invsqrt_pd(rsq_SSE1);
 914
 915             /* Apply mask */
 916             rinv_SSE0          = _mm_and_pd(rinv_SSE0, jmask_SSE0);
 917             rinv_SSE1          = _mm_and_pd(rinv_SSE1, jmask_SSE1);
 918
 919             irsq_SSE0          = _mm_mul_pd(rinv_SSE0, rinv_SSE0);
 920             irsq_SSE1          = _mm_mul_pd(rinv_SSE1, rinv_SSE1);
 921             idr4_SSE0          = _mm_mul_pd(irsq_SSE0, irsq_SSE0);
 922             idr4_SSE1          = _mm_mul_pd(irsq_SSE1, irsq_SSE1);
 923             idr6_SSE0          = _mm_mul_pd(idr4_SSE0, irsq_SSE0);
 924             idr6_SSE1          = _mm_mul_pd(idr4_SSE1, irsq_SSE1);
 925
 926             raj_SSE            = _mm_load_pd(gb_radius+j);
 927             vaj_SSE            = _mm_load_pd(vsolv+j);
 928
 929             rvdw_SSE0          = _mm_add_pd(rai_SSE0, raj_SSE);
 930             rvdw_SSE1          = _mm_add_pd(rai_SSE1, raj_SSE);
 931
 932             ratio_SSE0         = _mm_mul_pd(rsq_SSE0, gmx_mm_inv_pd( _mm_mul_pd(rvdw_SSE0, rvdw_SSE0)));
 933             ratio_SSE1         = _mm_mul_pd(rsq_SSE1, gmx_mm_inv_pd( _mm_mul_pd(rvdw_SSE1, rvdw_SSE1)));
 934
 935             ratio_SSE0         = _mm_min_pd(ratio_SSE0, still_p5inv_SSE);
 936             ratio_SSE1         = _mm_min_pd(ratio_SSE1, still_p5inv_SSE);
 937             theta_SSE0         = _mm_mul_pd(ratio_SSE0, still_pip5_SSE);
 938             theta_SSE1         = _mm_mul_pd(ratio_SSE1, still_pip5_SSE);
 939             gmx_mm_sincos_pd(theta_SSE0, &sinq_SSE0, &cosq_SSE0);
 940             gmx_mm_sincos_pd(theta_SSE1, &sinq_SSE1, &cosq_SSE1);
 941             term_SSE0          = _mm_mul_pd(half_SSE, _mm_sub_pd(one_SSE, cosq_SSE0));
 942             term_SSE1          = _mm_mul_pd(half_SSE, _mm_sub_pd(one_SSE, cosq_SSE1));
 943             ccf_SSE0           = _mm_mul_pd(term_SSE0, term_SSE0);
 944             ccf_SSE1           = _mm_mul_pd(term_SSE1, term_SSE1);
 945             dccf_SSE0          = _mm_mul_pd(_mm_mul_pd(two_SSE, term_SSE0),
 946                                             _mm_mul_pd(sinq_SSE0, theta_SSE0));
 947             dccf_SSE1          = _mm_mul_pd(_mm_mul_pd(two_SSE, term_SSE1),
 948                                             _mm_mul_pd(sinq_SSE1, theta_SSE1));
 949
 950             prod_SSE           = _mm_mul_pd(still_p4_SSE, vaj_SSE);
 951             icf4_SSE0          = _mm_mul_pd(ccf_SSE0, idr4_SSE0);
 952             icf4_SSE1          = _mm_mul_pd(ccf_SSE1, idr4_SSE1);
 953             icf6_SSE0          = _mm_mul_pd( _mm_sub_pd( _mm_mul_pd(four_SSE, ccf_SSE0), dccf_SSE0), idr6_SSE0);
 954             icf6_SSE1          = _mm_mul_pd( _mm_sub_pd( _mm_mul_pd(four_SSE, ccf_SSE1), dccf_SSE1), idr6_SSE1);
 955
 956             _mm_store_pd(work+j, _mm_add_pd(_mm_load_pd(work+j),
 957                                             _mm_add_pd(_mm_mul_pd(prod_ai_SSE0, icf4_SSE0),
 958                                                        _mm_mul_pd(prod_ai_SSE1, icf4_SSE1))));
 959
 960
 961             gpi_SSE0           = _mm_add_pd(gpi_SSE0, _mm_mul_pd(prod_SSE, icf4_SSE0));
 962             gpi_SSE1           = _mm_add_pd(gpi_SSE1, _mm_mul_pd(prod_SSE, icf4_SSE1));
 963
 964             /* Save ai->aj and aj->ai chain rule terms */
 965             _mm_store_pd(dadx, _mm_mul_pd(prod_SSE, icf6_SSE0));
 966             dadx += 2;
 967             _mm_store_pd(dadx, _mm_mul_pd(prod_SSE, icf6_SSE1));
 968             dadx += 2;
 969
 970             _mm_store_pd(dadx, _mm_mul_pd(prod_ai_SSE0, icf6_SSE0));
 971             dadx += 2;
 972             _mm_store_pd(dadx, _mm_mul_pd(prod_ai_SSE1, icf6_SSE1));
 973             dadx += 2;
 974         }
 975
 976         /* Main part, no exclusions */
 977         for (j = nj1; j < nj2; j += UNROLLJ)
 978         {
 979
 980             /* load j atom coordinates */
 981             jx_SSE            = _mm_load_pd(x_align+j);
 982             jy_SSE            = _mm_load_pd(y_align+j);
 983             jz_SSE            = _mm_load_pd(z_align+j);
 984
 985             /* Calculate distance */
 986             dx_SSE0            = _mm_sub_pd(ix_SSE0, jx_SSE);
 987             dy_SSE0            = _mm_sub_pd(iy_SSE0, jy_SSE);
 988             dz_SSE0            = _mm_sub_pd(iz_SSE0, jz_SSE);
 989             dx_SSE1            = _mm_sub_pd(ix_SSE1, jx_SSE);
 990             dy_SSE1            = _mm_sub_pd(iy_SSE1, jy_SSE);
 991             dz_SSE1            = _mm_sub_pd(iz_SSE1, jz_SSE);
 992
 993             /* rsq = dx*dx+dy*dy+dz*dz */
 994             rsq_SSE0           = gmx_mm_calc_rsq_pd(dx_SSE0, dy_SSE0, dz_SSE0);
 995             rsq_SSE1           = gmx_mm_calc_rsq_pd(dx_SSE1, dy_SSE1, dz_SSE1);
 996
 997             /* Calculate 1/r and 1/r2 */
 998             rinv_SSE0          = gmx_mm_invsqrt_pd(rsq_SSE0);
 999             rinv_SSE1          = gmx_mm_invsqrt_pd(rsq_SSE1);
1000
1001             /* Apply mask */
1002             rinv_SSE0          = _mm_and_pd(rinv_SSE0, imask_SSE0);
1003             rinv_SSE1          = _mm_and_pd(rinv_SSE1, imask_SSE1);
1004
1005             irsq_SSE0          = _mm_mul_pd(rinv_SSE0, rinv_SSE0);
1006             irsq_SSE1          = _mm_mul_pd(rinv_SSE1, rinv_SSE1);
1007             idr4_SSE0          = _mm_mul_pd(irsq_SSE0, irsq_SSE0);
1008             idr4_SSE1          = _mm_mul_pd(irsq_SSE1, irsq_SSE1);
1009             idr6_SSE0          = _mm_mul_pd(idr4_SSE0, irsq_SSE0);
1010             idr6_SSE1          = _mm_mul_pd(idr4_SSE1, irsq_SSE1);
1011
1012             raj_SSE            = _mm_load_pd(gb_radius+j);
1013
1014             rvdw_SSE0          = _mm_add_pd(rai_SSE0, raj_SSE);
1015             rvdw_SSE1          = _mm_add_pd(rai_SSE1, raj_SSE);
1016             vaj_SSE            = _mm_load_pd(vsolv+j);
1017
1018             ratio_SSE0         = _mm_mul_pd(rsq_SSE0, gmx_mm_inv_pd( _mm_mul_pd(rvdw_SSE0, rvdw_SSE0)));
1019             ratio_SSE1         = _mm_mul_pd(rsq_SSE1, gmx_mm_inv_pd( _mm_mul_pd(rvdw_SSE1, rvdw_SSE1)));
1020
1021             ratio_SSE0         = _mm_min_pd(ratio_SSE0, still_p5inv_SSE);
1022             ratio_SSE1         = _mm_min_pd(ratio_SSE1, still_p5inv_SSE);
1023             theta_SSE0         = _mm_mul_pd(ratio_SSE0, still_pip5_SSE);
1024             theta_SSE1         = _mm_mul_pd(ratio_SSE1, still_pip5_SSE);
1025             gmx_mm_sincos_pd(theta_SSE0, &sinq_SSE0, &cosq_SSE0);
1026             gmx_mm_sincos_pd(theta_SSE1, &sinq_SSE1, &cosq_SSE1);
1027             term_SSE0          = _mm_mul_pd(half_SSE, _mm_sub_pd(one_SSE, cosq_SSE0));
1028             term_SSE1          = _mm_mul_pd(half_SSE, _mm_sub_pd(one_SSE, cosq_SSE1));
1029             ccf_SSE0           = _mm_mul_pd(term_SSE0, term_SSE0);
1030             ccf_SSE1           = _mm_mul_pd(term_SSE1, term_SSE1);
1031             dccf_SSE0          = _mm_mul_pd(_mm_mul_pd(two_SSE, term_SSE0),
1032                                             _mm_mul_pd(sinq_SSE0, theta_SSE0));
1033             dccf_SSE1          = _mm_mul_pd(_mm_mul_pd(two_SSE, term_SSE1),
1034                                             _mm_mul_pd(sinq_SSE1, theta_SSE1));
1035
1036             prod_SSE           = _mm_mul_pd(still_p4_SSE, vaj_SSE );
1037             icf4_SSE0          = _mm_mul_pd(ccf_SSE0, idr4_SSE0);
1038             icf4_SSE1          = _mm_mul_pd(ccf_SSE1, idr4_SSE1);
1039             icf6_SSE0          = _mm_mul_pd( _mm_sub_pd( _mm_mul_pd(four_SSE, ccf_SSE0), dccf_SSE0), idr6_SSE0);
1040             icf6_SSE1          = _mm_mul_pd( _mm_sub_pd( _mm_mul_pd(four_SSE, ccf_SSE1), dccf_SSE1), idr6_SSE1);
1041
1042             _mm_store_pd(work+j, _mm_add_pd(_mm_load_pd(work+j),
1043                                             _mm_add_pd(_mm_mul_pd(prod_ai_SSE0, icf4_SSE0),
1044                                                        _mm_mul_pd(prod_ai_SSE1, icf4_SSE1))));
1045
1046             gpi_SSE0           = _mm_add_pd(gpi_SSE0, _mm_mul_pd(prod_SSE, icf4_SSE0));
1047             gpi_SSE1           = _mm_add_pd(gpi_SSE1, _mm_mul_pd(prod_SSE, icf4_SSE1));
1048
1049             /* Save ai->aj and aj->ai chain rule terms */
1050             _mm_store_pd(dadx, _mm_mul_pd(prod_SSE, icf6_SSE0));
1051             dadx += 2;
1052             _mm_store_pd(dadx, _mm_mul_pd(prod_SSE, icf6_SSE1));
1053             dadx += 2;
1054
1055             _mm_store_pd(dadx, _mm_mul_pd(prod_ai_SSE0, icf6_SSE0));
1056             dadx += 2;
1057             _mm_store_pd(dadx, _mm_mul_pd(prod_ai_SSE1, icf6_SSE1));
1058             dadx += 2;
1059         }
1060         /* Epilogue part, including exclusion mask */
1061         for (j = nj2; j < nj3; j += UNROLLJ)
1062         {
1063             jmask_SSE0 = _mm_load_pd((double *)emask0);
1064             jmask_SSE1 = _mm_load_pd((double *)emask1);
1065             emask0    += 2*UNROLLJ;
1066             emask1    += 2*UNROLLJ;
1067
1068             /* load j atom coordinates */
1069             jx_SSE            = _mm_load_pd(x_align+j);
1070             jy_SSE            = _mm_load_pd(y_align+j);
1071             jz_SSE            = _mm_load_pd(z_align+j);
1072
1073             /* Calculate distance */
1074             dx_SSE0            = _mm_sub_pd(ix_SSE0, jx_SSE);
1075             dy_SSE0            = _mm_sub_pd(iy_SSE0, jy_SSE);
1076             dz_SSE0            = _mm_sub_pd(iz_SSE0, jz_SSE);
1077             dx_SSE1            = _mm_sub_pd(ix_SSE1, jx_SSE);
1078             dy_SSE1            = _mm_sub_pd(iy_SSE1, jy_SSE);
1079             dz_SSE1            = _mm_sub_pd(iz_SSE1, jz_SSE);
1080
1081             /* rsq = dx*dx+dy*dy+dz*dz */
1082             rsq_SSE0           = gmx_mm_calc_rsq_pd(dx_SSE0, dy_SSE0, dz_SSE0);
1083             rsq_SSE1           = gmx_mm_calc_rsq_pd(dx_SSE1, dy_SSE1, dz_SSE1);
1084
1085             /* Combine masks */
1086             jmask_SSE0         = _mm_and_pd(jmask_SSE0, imask_SSE0);
1087             jmask_SSE1         = _mm_and_pd(jmask_SSE1, imask_SSE1);
1088
1089             /* Calculate 1/r and 1/r2 */
1090             rinv_SSE0          = gmx_mm_invsqrt_pd(rsq_SSE0);
1091             rinv_SSE1          = gmx_mm_invsqrt_pd(rsq_SSE1);
1092
1093             /* Apply mask */
1094             rinv_SSE0          = _mm_and_pd(rinv_SSE0, jmask_SSE0);
1095             rinv_SSE1          = _mm_and_pd(rinv_SSE1, jmask_SSE1);
1096
1097             irsq_SSE0          = _mm_mul_pd(rinv_SSE0, rinv_SSE0);
1098             irsq_SSE1          = _mm_mul_pd(rinv_SSE1, rinv_SSE1);
1099             idr4_SSE0          = _mm_mul_pd(irsq_SSE0, irsq_SSE0);
1100             idr4_SSE1          = _mm_mul_pd(irsq_SSE1, irsq_SSE1);
1101             idr6_SSE0          = _mm_mul_pd(idr4_SSE0, irsq_SSE0);
1102             idr6_SSE1          = _mm_mul_pd(idr4_SSE1, irsq_SSE1);
1103
1104             raj_SSE            = _mm_load_pd(gb_radius+j);
1105             vaj_SSE            = _mm_load_pd(vsolv+j);
1106
1107             rvdw_SSE0          = _mm_add_pd(rai_SSE0, raj_SSE);
1108             rvdw_SSE1          = _mm_add_pd(rai_SSE1, raj_SSE);
1109
1110             ratio_SSE0         = _mm_mul_pd(rsq_SSE0, gmx_mm_inv_pd( _mm_mul_pd(rvdw_SSE0, rvdw_SSE0)));
1111             ratio_SSE1         = _mm_mul_pd(rsq_SSE1, gmx_mm_inv_pd( _mm_mul_pd(rvdw_SSE1, rvdw_SSE1)));
1112
1113             ratio_SSE0         = _mm_min_pd(ratio_SSE0, still_p5inv_SSE);
1114             ratio_SSE1         = _mm_min_pd(ratio_SSE1, still_p5inv_SSE);
1115             theta_SSE0         = _mm_mul_pd(ratio_SSE0, still_pip5_SSE);
1116             theta_SSE1         = _mm_mul_pd(ratio_SSE1, still_pip5_SSE);
1117             gmx_mm_sincos_pd(theta_SSE0, &sinq_SSE0, &cosq_SSE0);
1118             gmx_mm_sincos_pd(theta_SSE1, &sinq_SSE1, &cosq_SSE1);
1119             term_SSE0          = _mm_mul_pd(half_SSE, _mm_sub_pd(one_SSE, cosq_SSE0));
1120             term_SSE1          = _mm_mul_pd(half_SSE, _mm_sub_pd(one_SSE, cosq_SSE1));
1121             ccf_SSE0           = _mm_mul_pd(term_SSE0, term_SSE0);
1122             ccf_SSE1           = _mm_mul_pd(term_SSE1, term_SSE1);
1123             dccf_SSE0          = _mm_mul_pd(_mm_mul_pd(two_SSE, term_SSE0),
1124                                             _mm_mul_pd(sinq_SSE0, theta_SSE0));
1125             dccf_SSE1          = _mm_mul_pd(_mm_mul_pd(two_SSE, term_SSE1),
1126                                             _mm_mul_pd(sinq_SSE1, theta_SSE1));
1127
1128             prod_SSE           = _mm_mul_pd(still_p4_SSE, vaj_SSE);
1129             icf4_SSE0          = _mm_mul_pd(ccf_SSE0, idr4_SSE0);
1130             icf4_SSE1          = _mm_mul_pd(ccf_SSE1, idr4_SSE1);
1131             icf6_SSE0          = _mm_mul_pd( _mm_sub_pd( _mm_mul_pd(four_SSE, ccf_SSE0), dccf_SSE0), idr6_SSE0);
1132             icf6_SSE1          = _mm_mul_pd( _mm_sub_pd( _mm_mul_pd(four_SSE, ccf_SSE1), dccf_SSE1), idr6_SSE1);
1133
1134             _mm_store_pd(work+j, _mm_add_pd(_mm_load_pd(work+j),
1135                                             _mm_add_pd(_mm_mul_pd(prod_ai_SSE0, icf4_SSE0),
1136                                                        _mm_mul_pd(prod_ai_SSE1, icf4_SSE1))));
1137
1138             gpi_SSE0           = _mm_add_pd(gpi_SSE0, _mm_mul_pd(prod_SSE, icf4_SSE0));
1139             gpi_SSE1           = _mm_add_pd(gpi_SSE1, _mm_mul_pd(prod_SSE, icf4_SSE1));
1140
1141             /* Save ai->aj and aj->ai chain rule terms */
1142             _mm_store_pd(dadx, _mm_mul_pd(prod_SSE, icf6_SSE0));
1143             dadx += 2;
1144             _mm_store_pd(dadx, _mm_mul_pd(prod_SSE, icf6_SSE1));
1145             dadx += 2;
1146
1147             _mm_store_pd(dadx, _mm_mul_pd(prod_ai_SSE0, icf6_SSE0));
1148             dadx += 2;
1149             _mm_store_pd(dadx, _mm_mul_pd(prod_ai_SSE1, icf6_SSE1));
1150             dadx += 2;
1151         }
1152         GMX_MM_TRANSPOSE2_PD(gpi_SSE0, gpi_SSE1);
1153         gpi_SSE0 = _mm_add_pd(gpi_SSE0, gpi_SSE1);
1154         _mm_store_pd(work+i, _mm_add_pd(gpi_SSE0, _mm_load_pd(work+i)));
1155     }
1156
1157     /* In case we have written anything beyond natoms, move it back.
1158      * Never mind that we leave stuff above natoms; that will not
1159      * be accessed later in the routine.
1160      * In principle this should be a move rather than sum, but this
1161      * way we dont have to worry about even/odd offsets...
1162      */
1163     for (i = natoms; i < ni1+1+natoms/2; i++)
1164     {
1165         work[i-natoms] += work[i];
1166     }
1167
1168     /* Parallel summations */
1169     if (PARTDECOMP(cr))
1170     {
1171         gmx_sum(natoms, work, cr);
1172     }
1173
1174     factor  = 0.5 * ONE_4PI_EPS0;
1175     /* Calculate the radii - should we do all atoms, or just our local ones? */
1176     for (i = 0; i < natoms; i++)
1177     {
1178         if (born->use[i] != 0)
1179         {
1180             gpi             = born->gpol[i]+work[i];
1181             gpi2            = gpi * gpi;
1182             born->bRad[i]   = factor*gmx_invsqrt(gpi2);
1183             fr->invsqrta[i] = gmx_invsqrt(born->bRad[i]);
1184         }
1185     }
1186
1187     return 0;
1188 }
1189 /* Reinstate MSVC optimization */
1190 #ifdef _MSC_VER
1191 #pragma optimize("",on)
1192 #endif
1193
1194
1195 int
1196 genborn_allvsall_calc_hct_obc_radii_sse2_double(t_forcerec   *           fr,
1197                                                 t_mdatoms   *            mdatoms,
1198                                                 gmx_genborn_t   *        born,
1199                                                 int                      gb_algorithm,
1200                                                 gmx_localtop_t   *       top,
1201                                                 double *                 x,
1202                                                 t_commrec   *            cr,
1203                                                 void   *                 paadata)
1204 {
1205     gmx_allvsallgb2_data_t *aadata;
1206     int                     natoms;
1207     int                     ni0, ni1;
1208     int                     nj0, nj1, nj2, nj3;
1209     int                     i, j, k, n;
1210     int              *      mask;
1211     int              *      pmask0;
1212     int              *      pmask1;
1213     int              *      emask0;
1214     int              *      emask1;
1215     double            *     gb_radius;
1216     double            *     vsolv;
1217     double            *     work;
1218     double                  tmpsum[2];
1219     double            *     x_align;
1220     double            *     y_align;
1221     double            *     z_align;
1222     int              *      jindex;
1223     double            *     dadx;
1224     double            *     obc_param;
1225     double                  rad, min_rad;
1226     double                  rai, rai_inv, rai_inv2, sum_ai, sum_ai2, sum_ai3, tsum, tchain;
1227
1228     __m128d                 ix_SSE0, iy_SSE0, iz_SSE0;
1229     __m128d                 ix_SSE1, iy_SSE1, iz_SSE1;
1230     __m128d                 gpi_SSE0, rai_SSE0, prod_ai_SSE0;
1231     __m128d                 gpi_SSE1, rai_SSE1, prod_ai_SSE1;
1232     __m128d                 imask_SSE0, jmask_SSE0;
1233     __m128d                 imask_SSE1, jmask_SSE1;
1234     __m128d                 jx_SSE, jy_SSE, jz_SSE;
1235     __m128d                 dx_SSE0, dy_SSE0, dz_SSE0;
1236     __m128d                 dx_SSE1, dy_SSE1, dz_SSE1;
1237     __m128d                 rsq_SSE0, rinv_SSE0, irsq_SSE0, idr4_SSE0, idr6_SSE0;
1238     __m128d                 rsq_SSE1, rinv_SSE1, irsq_SSE1, idr4_SSE1, idr6_SSE1;
1239     __m128d                 raj_SSE, raj_inv_SSE, sk_aj_SSE, sk2_aj_SSE;
1240     __m128d                 ccf_SSE0, dccf_SSE0, prod_SSE0;
1241     __m128d                 ccf_SSE1, dccf_SSE1, prod_SSE1;
1242     __m128d                 icf4_SSE0, icf6_SSE0;
1243     __m128d                 icf4_SSE1, icf6_SSE1;
1244     __m128d                 oneeighth_SSE, onefourth_SSE, half_SSE, one_SSE, two_SSE, four_SSE;
1245     __m128d                 still_p4_SSE, still_p5inv_SSE, still_pip5_SSE;
1246     __m128d                 rai_inv_SSE0;
1247     __m128d                 rai_inv_SSE1;
1248     __m128d                 sk_ai_SSE0, sk2_ai_SSE0, sum_ai_SSE0;
1249     __m128d                 sk_ai_SSE1, sk2_ai_SSE1, sum_ai_SSE1;
1250     __m128d                 lij_inv_SSE0, sk2_rinv_SSE0;
1251     __m128d                 lij_inv_SSE1, sk2_rinv_SSE1;
1252     __m128d                 dr_SSE0;
1253     __m128d                 dr_SSE1;
1254     __m128d                 t1_SSE0, t2_SSE0, t3_SSE0, t4_SSE0;
1255     __m128d                 t1_SSE1, t2_SSE1, t3_SSE1, t4_SSE1;
1256     __m128d                 obc_mask1_SSE0, obc_mask2_SSE0, obc_mask3_SSE0;
1257     __m128d                 obc_mask1_SSE1, obc_mask2_SSE1, obc_mask3_SSE1;
1258     __m128d                 uij_SSE0, uij2_SSE0, uij3_SSE0;
1259     __m128d                 uij_SSE1, uij2_SSE1, uij3_SSE1;
1260     __m128d                 lij_SSE0, lij2_SSE0, lij3_SSE0;
1261     __m128d                 lij_SSE1, lij2_SSE1, lij3_SSE1;
1262     __m128d                 dlij_SSE0, diff2_SSE0, logterm_SSE0;
1263     __m128d                 dlij_SSE1, diff2_SSE1, logterm_SSE1;
1264     __m128d                 doffset_SSE, tmpSSE;
1265
1266     natoms              = mdatoms->nr;
1267     ni0                 = (mdatoms->start/SIMD_WIDTH)*SIMD_WIDTH;
1268     ni1                 = mdatoms->start+mdatoms->homenr;
1269
1270     n = 0;
1271
1272     aadata = *((gmx_allvsallgb2_data_t **)paadata);
1273
1274
1275     if (aadata == NULL)
1276     {
1277         genborn_allvsall_setup(&aadata, top, born, mdatoms, born->gb_doffset,
1278                                egbOBC, TRUE, TRUE, TRUE);
1279         *((gmx_allvsallgb2_data_t **)paadata) = aadata;
1280     }
1281
1282     x_align = aadata->x_align;
1283     y_align = aadata->y_align;
1284     z_align = aadata->z_align;
1285
1286     gb_radius = aadata->gb_radius;
1287     work      = aadata->work;
1288     jindex    = aadata->jindex_gb;
1289     dadx      = fr->dadx;
1290     obc_param = aadata->workparam;
1291
1292     oneeighth_SSE   = _mm_set1_pd(0.125);
1293     onefourth_SSE   = _mm_set1_pd(0.25);
1294     half_SSE        = _mm_set1_pd(0.5);
1295     one_SSE         = _mm_set1_pd(1.0);
1296     two_SSE         = _mm_set1_pd(2.0);
1297     four_SSE        = _mm_set1_pd(4.0);
1298     doffset_SSE     = _mm_set1_pd(born->gb_doffset);
1299
1300     for (i = 0; i < natoms; i++)
1301     {
1302         x_align[i]  = x[3*i];
1303         y_align[i]  = x[3*i+1];
1304         z_align[i]  = x[3*i+2];
1305     }
1306
1307     /* Copy again */
1308     for (i = 0; i < natoms/2+1; i++)
1309     {
1310         x_align[natoms+i]  = x_align[i];
1311         y_align[natoms+i]  = y_align[i];
1312         z_align[natoms+i]  = z_align[i];
1313     }
1314
1315     for (i = 0; i < natoms+natoms/2+1; i++)
1316     {
1317         work[i] = 0;
1318     }
1319
1320     for (i = ni0; i < ni1; i += UNROLLI)
1321     {
1322         /* We assume shifts are NOT used for all-vs-all interactions */
1323
1324         /* Load i atom data */
1325         ix_SSE0          = _mm_load1_pd(x_align+i);
1326         iy_SSE0          = _mm_load1_pd(y_align+i);
1327         iz_SSE0          = _mm_load1_pd(z_align+i);
1328         ix_SSE1          = _mm_load1_pd(x_align+i+1);
1329         iy_SSE1          = _mm_load1_pd(y_align+i+1);
1330         iz_SSE1          = _mm_load1_pd(z_align+i+1);
1331
1332         rai_SSE0         = _mm_load1_pd(gb_radius+i);
1333         rai_SSE1         = _mm_load1_pd(gb_radius+i+1);
1334         rai_inv_SSE0     = gmx_mm_inv_pd(rai_SSE0);
1335         rai_inv_SSE1     = gmx_mm_inv_pd(rai_SSE1);
1336
1337         sk_ai_SSE0       = _mm_load1_pd(obc_param+i);
1338         sk_ai_SSE1       = _mm_load1_pd(obc_param+i+1);
1339         sk2_ai_SSE0      = _mm_mul_pd(sk_ai_SSE0, sk_ai_SSE0);
1340         sk2_ai_SSE1      = _mm_mul_pd(sk_ai_SSE1, sk_ai_SSE1);
1341
1342         sum_ai_SSE0      = _mm_setzero_pd();
1343         sum_ai_SSE1      = _mm_setzero_pd();
1344
1345         /* Load limits for loop over neighbors */
1346         nj0              = jindex[4*i];
1347         nj1              = jindex[4*i+1];
1348         nj2              = jindex[4*i+2];
1349         nj3              = jindex[4*i+3];
1350
1351         pmask0           = aadata->prologue_mask_gb[i];
1352         pmask1           = aadata->prologue_mask_gb[i+1];
1353         emask0           = aadata->epilogue_mask[i];
1354         emask1           = aadata->epilogue_mask[i+1];
1355
1356         imask_SSE0        = _mm_load1_pd((double *)(aadata->imask+2*i));
1357         imask_SSE1        = _mm_load1_pd((double *)(aadata->imask+2*i+2));
1358
1359         /* Prologue part, including exclusion mask */
1360         for (j = nj0; j < nj1; j += UNROLLJ)
1361         {
1362             jmask_SSE0 = _mm_load_pd((double *)pmask0);
1363             jmask_SSE1 = _mm_load_pd((double *)pmask1);
1364             pmask0    += 2*UNROLLJ;
1365             pmask1    += 2*UNROLLJ;
1366
1367             /* load j atom coordinates */
1368             jx_SSE            = _mm_load_pd(x_align+j);
1369             jy_SSE            = _mm_load_pd(y_align+j);
1370             jz_SSE            = _mm_load_pd(z_align+j);
1371
1372             /* Calculate distance */
1373             dx_SSE0            = _mm_sub_pd(ix_SSE0, jx_SSE);
1374             dy_SSE0            = _mm_sub_pd(iy_SSE0, jy_SSE);
1375             dz_SSE0            = _mm_sub_pd(iz_SSE0, jz_SSE);
1376             dx_SSE1            = _mm_sub_pd(ix_SSE1, jx_SSE);
1377             dy_SSE1            = _mm_sub_pd(iy_SSE1, jy_SSE);
1378             dz_SSE1            = _mm_sub_pd(iz_SSE1, jz_SSE);
1379
1380             /* rsq = dx*dx+dy*dy+dz*dz */
1381             rsq_SSE0           = gmx_mm_calc_rsq_pd(dx_SSE0, dy_SSE0, dz_SSE0);
1382             rsq_SSE1           = gmx_mm_calc_rsq_pd(dx_SSE1, dy_SSE1, dz_SSE1);
1383
1384             /* Combine masks */
1385             jmask_SSE0         = _mm_and_pd(jmask_SSE0, imask_SSE0);
1386             jmask_SSE1         = _mm_and_pd(jmask_SSE1, imask_SSE1);
1387
1388             /* Calculate 1/r and 1/r2 */
1389             rinv_SSE0          = gmx_mm_invsqrt_pd(rsq_SSE0);
1390             rinv_SSE1          = gmx_mm_invsqrt_pd(rsq_SSE1);
1391
1392             /* Apply mask */
1393             rinv_SSE0          = _mm_and_pd(rinv_SSE0, jmask_SSE0);
1394             rinv_SSE1          = _mm_and_pd(rinv_SSE1, jmask_SSE1);
1395
1396             dr_SSE0            = _mm_mul_pd(rsq_SSE0, rinv_SSE0);
1397             dr_SSE1            = _mm_mul_pd(rsq_SSE1, rinv_SSE1);
1398
1399             sk_aj_SSE          = _mm_load_pd(obc_param+j);
1400             raj_SSE            = _mm_load_pd(gb_radius+j);
1401             raj_inv_SSE        = gmx_mm_inv_pd(raj_SSE);
1402
1403             /* Evaluate influence of atom aj -> ai */
1404             t1_SSE0            = _mm_add_pd(dr_SSE0, sk_aj_SSE);
1405             t1_SSE1            = _mm_add_pd(dr_SSE1, sk_aj_SSE);
1406             t2_SSE0            = _mm_sub_pd(dr_SSE0, sk_aj_SSE);
1407             t2_SSE1            = _mm_sub_pd(dr_SSE1, sk_aj_SSE);
1408             t3_SSE0            = _mm_sub_pd(sk_aj_SSE, dr_SSE0);
1409             t3_SSE1            = _mm_sub_pd(sk_aj_SSE, dr_SSE1);
1410
1411             obc_mask1_SSE0     = _mm_cmplt_pd(rai_SSE0, t1_SSE0);
1412             obc_mask1_SSE1     = _mm_cmplt_pd(rai_SSE1, t1_SSE1);
1413             obc_mask2_SSE0     = _mm_cmplt_pd(rai_SSE0, t2_SSE0);
1414             obc_mask2_SSE1     = _mm_cmplt_pd(rai_SSE1, t2_SSE1);
1415             obc_mask3_SSE0     = _mm_cmplt_pd(rai_SSE0, t3_SSE0);
1416             obc_mask3_SSE1     = _mm_cmplt_pd(rai_SSE1, t3_SSE1);
1417             obc_mask1_SSE0     = _mm_and_pd(obc_mask1_SSE0, jmask_SSE0);
1418             obc_mask1_SSE1     = _mm_and_pd(obc_mask1_SSE1, jmask_SSE1);
1419
1420             uij_SSE0           = gmx_mm_inv_pd(t1_SSE0);
1421             uij_SSE1           = gmx_mm_inv_pd(t1_SSE1);
1422             lij_SSE0           = _mm_or_pd(   _mm_and_pd(obc_mask2_SSE0, gmx_mm_inv_pd(t2_SSE0)),
1423                                               _mm_andnot_pd(obc_mask2_SSE0, rai_inv_SSE0));
1424             lij_SSE1           = _mm_or_pd(   _mm_and_pd(obc_mask2_SSE1, gmx_mm_inv_pd(t2_SSE1)),
1425                                               _mm_andnot_pd(obc_mask2_SSE1, rai_inv_SSE1));
1426             dlij_SSE0          = _mm_and_pd(one_SSE, obc_mask2_SSE0);
1427             dlij_SSE1          = _mm_and_pd(one_SSE, obc_mask2_SSE1);
1428
1429             uij2_SSE0          = _mm_mul_pd(uij_SSE0, uij_SSE0);
1430             uij2_SSE1          = _mm_mul_pd(uij_SSE1, uij_SSE1);
1431             uij3_SSE0          = _mm_mul_pd(uij2_SSE0, uij_SSE0);
1432             uij3_SSE1          = _mm_mul_pd(uij2_SSE1, uij_SSE1);
1433             lij2_SSE0          = _mm_mul_pd(lij_SSE0, lij_SSE0);
1434             lij2_SSE1          = _mm_mul_pd(lij_SSE1, lij_SSE1);
1435             lij3_SSE0          = _mm_mul_pd(lij2_SSE0, lij_SSE0);
1436             lij3_SSE1          = _mm_mul_pd(lij2_SSE1, lij_SSE1);
1437
1438             diff2_SSE0         = _mm_sub_pd(uij2_SSE0, lij2_SSE0);
1439             diff2_SSE1         = _mm_sub_pd(uij2_SSE1, lij2_SSE1);
1440             lij_inv_SSE0       = gmx_mm_invsqrt_pd(lij2_SSE0);
1441             lij_inv_SSE1       = gmx_mm_invsqrt_pd(lij2_SSE1);
1442             sk2_aj_SSE         = _mm_mul_pd(sk_aj_SSE, sk_aj_SSE);
1443             sk2_rinv_SSE0      = _mm_mul_pd(sk2_aj_SSE, rinv_SSE0);
1444             sk2_rinv_SSE1      = _mm_mul_pd(sk2_aj_SSE, rinv_SSE1);
1445             prod_SSE0          = _mm_mul_pd(onefourth_SSE, sk2_rinv_SSE0);
1446             prod_SSE1          = _mm_mul_pd(onefourth_SSE, sk2_rinv_SSE1);
1447
1448             logterm_SSE0       = gmx_mm_log_pd(_mm_mul_pd(uij_SSE0, lij_inv_SSE0));
1449             logterm_SSE1       = gmx_mm_log_pd(_mm_mul_pd(uij_SSE1, lij_inv_SSE1));
1450
1451             t1_SSE0            = _mm_sub_pd(lij_SSE0, uij_SSE0);
1452             t1_SSE1            = _mm_sub_pd(lij_SSE1, uij_SSE1);
1453             t2_SSE0            = _mm_mul_pd(diff2_SSE0,
1454                                             _mm_sub_pd(_mm_mul_pd(onefourth_SSE, dr_SSE0),
1455                                                        prod_SSE0));
1456             t2_SSE1            = _mm_mul_pd(diff2_SSE1,
1457                                             _mm_sub_pd(_mm_mul_pd(onefourth_SSE, dr_SSE1),
1458                                                        prod_SSE1));
1459
1460             t3_SSE0            = _mm_mul_pd(half_SSE, _mm_mul_pd(rinv_SSE0, logterm_SSE0));
1461             t3_SSE1            = _mm_mul_pd(half_SSE, _mm_mul_pd(rinv_SSE1, logterm_SSE1));
1462             t1_SSE0            = _mm_add_pd(t1_SSE0, _mm_add_pd(t2_SSE0, t3_SSE0));
1463             t1_SSE1            = _mm_add_pd(t1_SSE1, _mm_add_pd(t2_SSE1, t3_SSE1));
1464             t4_SSE0            = _mm_mul_pd(two_SSE, _mm_sub_pd(rai_inv_SSE0, lij_SSE0));
1465             t4_SSE1            = _mm_mul_pd(two_SSE, _mm_sub_pd(rai_inv_SSE1, lij_SSE1));
1466             t4_SSE0            = _mm_and_pd(t4_SSE0, obc_mask3_SSE0);
1467             t4_SSE1            = _mm_and_pd(t4_SSE1, obc_mask3_SSE1);
1468             t1_SSE0            = _mm_mul_pd(half_SSE, _mm_add_pd(t1_SSE0, t4_SSE0));
1469             t1_SSE1            = _mm_mul_pd(half_SSE, _mm_add_pd(t1_SSE1, t4_SSE1));
1470
1471             sum_ai_SSE0        = _mm_add_pd(sum_ai_SSE0, _mm_and_pd(t1_SSE0, obc_mask1_SSE0));
1472             sum_ai_SSE1        = _mm_add_pd(sum_ai_SSE1, _mm_and_pd(t1_SSE1, obc_mask1_SSE1));
1473
1474             t1_SSE0            = _mm_add_pd(_mm_mul_pd(half_SSE, lij2_SSE0),
1475                                             _mm_mul_pd(prod_SSE0, lij3_SSE0));
1476             t1_SSE1            = _mm_add_pd(_mm_mul_pd(half_SSE, lij2_SSE1),
1477                                             _mm_mul_pd(prod_SSE1, lij3_SSE1));
1478             t1_SSE0            = _mm_sub_pd(t1_SSE0,
1479                                             _mm_mul_pd(onefourth_SSE,
1480                                                        _mm_add_pd(_mm_mul_pd(lij_SSE0, rinv_SSE0),
1481                                                                   _mm_mul_pd(lij3_SSE0, dr_SSE0))));
1482             t1_SSE1            = _mm_sub_pd(t1_SSE1,
1483                                             _mm_mul_pd(onefourth_SSE,
1484                                                        _mm_add_pd(_mm_mul_pd(lij_SSE1, rinv_SSE1),
1485                                                                   _mm_mul_pd(lij3_SSE1, dr_SSE1))));
1486
1487             t2_SSE0            = _mm_mul_pd(onefourth_SSE,
1488                                             _mm_add_pd(_mm_mul_pd(uij_SSE0, rinv_SSE0),
1489                                                        _mm_mul_pd(uij3_SSE0, dr_SSE0)));
1490             t2_SSE1            = _mm_mul_pd(onefourth_SSE,
1491                                             _mm_add_pd(_mm_mul_pd(uij_SSE1, rinv_SSE1),
1492                                                        _mm_mul_pd(uij3_SSE1, dr_SSE1)));
1493             t2_SSE0            = _mm_sub_pd(t2_SSE0,
1494                                             _mm_add_pd(_mm_mul_pd(half_SSE, uij2_SSE0),
1495                                                        _mm_mul_pd(prod_SSE0, uij3_SSE0)));
1496             t2_SSE1            = _mm_sub_pd(t2_SSE1,
1497                                             _mm_add_pd(_mm_mul_pd(half_SSE, uij2_SSE1),
1498                                                        _mm_mul_pd(prod_SSE1, uij3_SSE1)));
1499             t3_SSE0            = _mm_mul_pd(_mm_mul_pd(onefourth_SSE, logterm_SSE0),
1500                                             _mm_mul_pd(rinv_SSE0, rinv_SSE0));
1501             t3_SSE1            = _mm_mul_pd(_mm_mul_pd(onefourth_SSE, logterm_SSE1),
1502                                             _mm_mul_pd(rinv_SSE1, rinv_SSE1));
1503             t3_SSE0            = _mm_sub_pd(t3_SSE0,
1504                                             _mm_mul_pd(_mm_mul_pd(diff2_SSE0, oneeighth_SSE),
1505                                                        _mm_add_pd(one_SSE,
1506                                                                   _mm_mul_pd(sk2_rinv_SSE0, rinv_SSE0))));
1507             t3_SSE1            = _mm_sub_pd(t3_SSE1,
1508                                             _mm_mul_pd(_mm_mul_pd(diff2_SSE1, oneeighth_SSE),
1509                                                        _mm_add_pd(one_SSE,
1510                                                                   _mm_mul_pd(sk2_rinv_SSE1, rinv_SSE1))));
1511
1512             t1_SSE0            = _mm_mul_pd(rinv_SSE0,
1513                                             _mm_add_pd(_mm_mul_pd(dlij_SSE0, t1_SSE0),
1514                                                        _mm_add_pd(t2_SSE0, t3_SSE0)));
1515             t1_SSE1            = _mm_mul_pd(rinv_SSE1,
1516                                             _mm_add_pd(_mm_mul_pd(dlij_SSE1, t1_SSE1),
1517                                                        _mm_add_pd(t2_SSE1, t3_SSE1)));
1518
1519             _mm_store_pd(dadx, _mm_and_pd(t1_SSE0, obc_mask1_SSE0));
1520             dadx += 2;
1521             _mm_store_pd(dadx, _mm_and_pd(t1_SSE1, obc_mask1_SSE1));
1522             dadx += 2;
1523
1524             /* Evaluate influence of atom ai -> aj */
1525             t1_SSE0            = _mm_add_pd(dr_SSE0, sk_ai_SSE0);
1526             t1_SSE1            = _mm_add_pd(dr_SSE1, sk_ai_SSE1);
1527             t2_SSE0            = _mm_sub_pd(dr_SSE0, sk_ai_SSE0);
1528             t2_SSE1            = _mm_sub_pd(dr_SSE1, sk_ai_SSE1);
1529             t3_SSE0            = _mm_sub_pd(sk_ai_SSE0, dr_SSE0);
1530             t3_SSE1            = _mm_sub_pd(sk_ai_SSE1, dr_SSE1);
1531
1532             obc_mask1_SSE0     = _mm_cmplt_pd(raj_SSE, t1_SSE0);
1533             obc_mask1_SSE1     = _mm_cmplt_pd(raj_SSE, t1_SSE1);
1534             obc_mask2_SSE0     = _mm_cmplt_pd(raj_SSE, t2_SSE0);
1535             obc_mask2_SSE1     = _mm_cmplt_pd(raj_SSE, t2_SSE1);
1536             obc_mask3_SSE0     = _mm_cmplt_pd(raj_SSE, t3_SSE0);
1537             obc_mask3_SSE1     = _mm_cmplt_pd(raj_SSE, t3_SSE1);
1538             obc_mask1_SSE0     = _mm_and_pd(obc_mask1_SSE0, jmask_SSE0);
1539             obc_mask1_SSE1     = _mm_and_pd(obc_mask1_SSE1, jmask_SSE1);
1540
1541             uij_SSE0           = gmx_mm_inv_pd(t1_SSE0);
1542             uij_SSE1           = gmx_mm_inv_pd(t1_SSE1);
1543             lij_SSE0           = _mm_or_pd(   _mm_and_pd(obc_mask2_SSE0, gmx_mm_inv_pd(t2_SSE0)),
1544                                               _mm_andnot_pd(obc_mask2_SSE0, raj_inv_SSE));
1545             lij_SSE1           = _mm_or_pd(   _mm_and_pd(obc_mask2_SSE1, gmx_mm_inv_pd(t2_SSE1)),
1546                                               _mm_andnot_pd(obc_mask2_SSE1, raj_inv_SSE));
1547             dlij_SSE0          = _mm_and_pd(one_SSE, obc_mask2_SSE0);
1548             dlij_SSE1          = _mm_and_pd(one_SSE, obc_mask2_SSE1);
1549
1550             uij2_SSE0          = _mm_mul_pd(uij_SSE0, uij_SSE0);
1551             uij2_SSE1          = _mm_mul_pd(uij_SSE1, uij_SSE1);
1552             uij3_SSE0          = _mm_mul_pd(uij2_SSE0, uij_SSE0);
1553             uij3_SSE1          = _mm_mul_pd(uij2_SSE1, uij_SSE1);
1554             lij2_SSE0          = _mm_mul_pd(lij_SSE0, lij_SSE0);
1555             lij2_SSE1          = _mm_mul_pd(lij_SSE1, lij_SSE1);
1556             lij3_SSE0          = _mm_mul_pd(lij2_SSE0, lij_SSE0);
1557             lij3_SSE1          = _mm_mul_pd(lij2_SSE1, lij_SSE1);
1558
1559             diff2_SSE0         = _mm_sub_pd(uij2_SSE0, lij2_SSE0);
1560             diff2_SSE1         = _mm_sub_pd(uij2_SSE1, lij2_SSE1);
1561             lij_inv_SSE0       = gmx_mm_invsqrt_pd(lij2_SSE0);
1562             lij_inv_SSE1       = gmx_mm_invsqrt_pd(lij2_SSE1);
1563             sk2_rinv_SSE0      = _mm_mul_pd(sk2_ai_SSE0, rinv_SSE0);
1564             sk2_rinv_SSE1      = _mm_mul_pd(sk2_ai_SSE1, rinv_SSE1);
1565             prod_SSE0          = _mm_mul_pd(onefourth_SSE, sk2_rinv_SSE0);
1566             prod_SSE1          = _mm_mul_pd(onefourth_SSE, sk2_rinv_SSE1);
1567
1568             logterm_SSE0       = gmx_mm_log_pd(_mm_mul_pd(uij_SSE0, lij_inv_SSE0));
1569             logterm_SSE1       = gmx_mm_log_pd(_mm_mul_pd(uij_SSE1, lij_inv_SSE1));
1570             t1_SSE0            = _mm_sub_pd(lij_SSE0, uij_SSE0);
1571             t1_SSE1            = _mm_sub_pd(lij_SSE1, uij_SSE1);
1572             t2_SSE0            = _mm_mul_pd(diff2_SSE0,
1573                                             _mm_sub_pd(_mm_mul_pd(onefourth_SSE, dr_SSE0),
1574                                                        prod_SSE0));
1575             t2_SSE1            = _mm_mul_pd(diff2_SSE1,
1576                                             _mm_sub_pd(_mm_mul_pd(onefourth_SSE, dr_SSE1),
1577                                                        prod_SSE1));
1578             t3_SSE0            = _mm_mul_pd(half_SSE, _mm_mul_pd(rinv_SSE0, logterm_SSE0));
1579             t3_SSE1            = _mm_mul_pd(half_SSE, _mm_mul_pd(rinv_SSE1, logterm_SSE1));
1580             t1_SSE0            = _mm_add_pd(t1_SSE0, _mm_add_pd(t2_SSE0, t3_SSE0));
1581             t1_SSE1            = _mm_add_pd(t1_SSE1, _mm_add_pd(t2_SSE1, t3_SSE1));
1582             t4_SSE0            = _mm_mul_pd(two_SSE, _mm_sub_pd(raj_inv_SSE, lij_SSE0));
1583             t4_SSE1            = _mm_mul_pd(two_SSE, _mm_sub_pd(raj_inv_SSE, lij_SSE1));
1584             t4_SSE0            = _mm_and_pd(t4_SSE0, obc_mask3_SSE0);
1585             t4_SSE1            = _mm_and_pd(t4_SSE1, obc_mask3_SSE1);
1586             t1_SSE0            = _mm_mul_pd(half_SSE, _mm_add_pd(t1_SSE0, t4_SSE0));
1587             t1_SSE1            = _mm_mul_pd(half_SSE, _mm_add_pd(t1_SSE1, t4_SSE1));
1588
1589             _mm_store_pd(work+j, _mm_add_pd(_mm_load_pd(work+j),
1590                                             _mm_add_pd(_mm_and_pd(t1_SSE0, obc_mask1_SSE0),
1591                                                        _mm_and_pd(t1_SSE1, obc_mask1_SSE1))));
1592
1593             t1_SSE0            = _mm_add_pd(_mm_mul_pd(half_SSE, lij2_SSE0),
1594                                             _mm_mul_pd(prod_SSE0, lij3_SSE0));
1595             t1_SSE1            = _mm_add_pd(_mm_mul_pd(half_SSE, lij2_SSE1),
1596                                             _mm_mul_pd(prod_SSE1, lij3_SSE1));
1597             t1_SSE0            = _mm_sub_pd(t1_SSE0,
1598                                             _mm_mul_pd(onefourth_SSE,
1599                                                        _mm_add_pd(_mm_mul_pd(lij_SSE0, rinv_SSE0),
1600                                                                   _mm_mul_pd(lij3_SSE0, dr_SSE0))));
1601             t1_SSE1            = _mm_sub_pd(t1_SSE1,
1602                                             _mm_mul_pd(onefourth_SSE,
1603                                                        _mm_add_pd(_mm_mul_pd(lij_SSE1, rinv_SSE1),
1604                                                                   _mm_mul_pd(lij3_SSE1, dr_SSE1))));
1605             t2_SSE0            = _mm_mul_pd(onefourth_SSE,
1606                                             _mm_add_pd(_mm_mul_pd(uij_SSE0, rinv_SSE0),
1607                                                        _mm_mul_pd(uij3_SSE0, dr_SSE0)));
1608             t2_SSE1            = _mm_mul_pd(onefourth_SSE,
1609                                             _mm_add_pd(_mm_mul_pd(uij_SSE1, rinv_SSE1),
1610                                                        _mm_mul_pd(uij3_SSE1, dr_SSE1)));
1611             t2_SSE0            = _mm_sub_pd(t2_SSE0,
1612                                             _mm_add_pd(_mm_mul_pd(half_SSE, uij2_SSE0),
1613                                                        _mm_mul_pd(prod_SSE0, uij3_SSE0)));
1614             t2_SSE1            = _mm_sub_pd(t2_SSE1,
1615                                             _mm_add_pd(_mm_mul_pd(half_SSE, uij2_SSE1),
1616                                                        _mm_mul_pd(prod_SSE1, uij3_SSE1)));
1617
1618             t3_SSE0            = _mm_mul_pd(_mm_mul_pd(onefourth_SSE, logterm_SSE0),
1619                                             _mm_mul_pd(rinv_SSE0, rinv_SSE0));
1620             t3_SSE1            = _mm_mul_pd(_mm_mul_pd(onefourth_SSE, logterm_SSE1),
1621                                             _mm_mul_pd(rinv_SSE1, rinv_SSE1));
1622
1623             t3_SSE0            = _mm_sub_pd(t3_SSE0,
1624                                             _mm_mul_pd(_mm_mul_pd(diff2_SSE0, oneeighth_SSE),
1625                                                        _mm_add_pd(one_SSE,
1626                                                                   _mm_mul_pd(sk2_rinv_SSE0, rinv_SSE0))));
1627             t3_SSE1            = _mm_sub_pd(t3_SSE1,
1628                                             _mm_mul_pd(_mm_mul_pd(diff2_SSE1, oneeighth_SSE),
1629                                                        _mm_add_pd(one_SSE,
1630                                                                   _mm_mul_pd(sk2_rinv_SSE1, rinv_SSE1))));
1631
1632
1633             t1_SSE0            = _mm_mul_pd(rinv_SSE0,
1634                                             _mm_add_pd(_mm_mul_pd(dlij_SSE0, t1_SSE0),
1635                                                        _mm_add_pd(t2_SSE0, t3_SSE0)));
1636             t1_SSE1            = _mm_mul_pd(rinv_SSE1,
1637                                             _mm_add_pd(_mm_mul_pd(dlij_SSE1, t1_SSE1),
1638                                                        _mm_add_pd(t2_SSE1, t3_SSE1)));
1639
1640             _mm_store_pd(dadx, _mm_and_pd(t1_SSE0, obc_mask1_SSE0));
1641             dadx += 2;
1642             _mm_store_pd(dadx, _mm_and_pd(t1_SSE1, obc_mask1_SSE1));
1643             dadx += 2;
1644         }
1645
1646         /* Main part, no exclusions */
1647         for (j = nj1; j < nj2; j += UNROLLJ)
1648         {
1649             /* load j atom coordinates */
1650             jx_SSE            = _mm_load_pd(x_align+j);
1651             jy_SSE            = _mm_load_pd(y_align+j);
1652             jz_SSE            = _mm_load_pd(z_align+j);
1653
1654             /* Calculate distance */
1655             dx_SSE0            = _mm_sub_pd(ix_SSE0, jx_SSE);
1656             dy_SSE0            = _mm_sub_pd(iy_SSE0, jy_SSE);
1657             dz_SSE0            = _mm_sub_pd(iz_SSE0, jz_SSE);
1658             dx_SSE1            = _mm_sub_pd(ix_SSE1, jx_SSE);
1659             dy_SSE1            = _mm_sub_pd(iy_SSE1, jy_SSE);
1660             dz_SSE1            = _mm_sub_pd(iz_SSE1, jz_SSE);
1661
1662             /* rsq = dx*dx+dy*dy+dz*dz */
1663             rsq_SSE0           = gmx_mm_calc_rsq_pd(dx_SSE0, dy_SSE0, dz_SSE0);
1664             rsq_SSE1           = gmx_mm_calc_rsq_pd(dx_SSE1, dy_SSE1, dz_SSE1);
1665
1666             /* Calculate 1/r and 1/r2 */
1667             rinv_SSE0          = gmx_mm_invsqrt_pd(rsq_SSE0);
1668             rinv_SSE1          = gmx_mm_invsqrt_pd(rsq_SSE1);
1669
1670             /* Apply mask */
1671             rinv_SSE0          = _mm_and_pd(rinv_SSE0, imask_SSE0);
1672             rinv_SSE1          = _mm_and_pd(rinv_SSE1, imask_SSE1);
1673
1674             dr_SSE0            = _mm_mul_pd(rsq_SSE0, rinv_SSE0);
1675             dr_SSE1            = _mm_mul_pd(rsq_SSE1, rinv_SSE1);
1676
1677             sk_aj_SSE          = _mm_load_pd(obc_param+j);
1678             raj_SSE            = _mm_load_pd(gb_radius+j);
1679
1680             raj_inv_SSE        = gmx_mm_inv_pd(raj_SSE);
1681
1682             /* Evaluate influence of atom aj -> ai */
1683             t1_SSE0            = _mm_add_pd(dr_SSE0, sk_aj_SSE);
1684             t1_SSE1            = _mm_add_pd(dr_SSE1, sk_aj_SSE);
1685             t2_SSE0            = _mm_sub_pd(dr_SSE0, sk_aj_SSE);
1686             t2_SSE1            = _mm_sub_pd(dr_SSE1, sk_aj_SSE);
1687             t3_SSE0            = _mm_sub_pd(sk_aj_SSE, dr_SSE0);
1688             t3_SSE1            = _mm_sub_pd(sk_aj_SSE, dr_SSE1);
1689
1690             obc_mask1_SSE0     = _mm_cmplt_pd(rai_SSE0, t1_SSE0);
1691             obc_mask1_SSE1     = _mm_cmplt_pd(rai_SSE1, t1_SSE1);
1692             obc_mask2_SSE0     = _mm_cmplt_pd(rai_SSE0, t2_SSE0);
1693             obc_mask2_SSE1     = _mm_cmplt_pd(rai_SSE1, t2_SSE1);
1694             obc_mask3_SSE0     = _mm_cmplt_pd(rai_SSE0, t3_SSE0);
1695             obc_mask3_SSE1     = _mm_cmplt_pd(rai_SSE1, t3_SSE1);
1696             obc_mask1_SSE0     = _mm_and_pd(obc_mask1_SSE0, imask_SSE0);
1697             obc_mask1_SSE1     = _mm_and_pd(obc_mask1_SSE1, imask_SSE1);
1698
1699             uij_SSE0           = gmx_mm_inv_pd(t1_SSE0);
1700             uij_SSE1           = gmx_mm_inv_pd(t1_SSE1);
1701             lij_SSE0           = _mm_or_pd(   _mm_and_pd(obc_mask2_SSE0, gmx_mm_inv_pd(t2_SSE0)),
1702                                               _mm_andnot_pd(obc_mask2_SSE0, rai_inv_SSE0));
1703             lij_SSE1           = _mm_or_pd(   _mm_and_pd(obc_mask2_SSE1, gmx_mm_inv_pd(t2_SSE1)),
1704                                               _mm_andnot_pd(obc_mask2_SSE1, rai_inv_SSE1));
1705             dlij_SSE0          = _mm_and_pd(one_SSE, obc_mask2_SSE0);
1706             dlij_SSE1          = _mm_and_pd(one_SSE, obc_mask2_SSE1);
1707
1708             uij2_SSE0          = _mm_mul_pd(uij_SSE0, uij_SSE0);
1709             uij2_SSE1          = _mm_mul_pd(uij_SSE1, uij_SSE1);
1710             uij3_SSE0          = _mm_mul_pd(uij2_SSE0, uij_SSE0);
1711             uij3_SSE1          = _mm_mul_pd(uij2_SSE1, uij_SSE1);
1712             lij2_SSE0          = _mm_mul_pd(lij_SSE0, lij_SSE0);
1713             lij2_SSE1          = _mm_mul_pd(lij_SSE1, lij_SSE1);
1714             lij3_SSE0          = _mm_mul_pd(lij2_SSE0, lij_SSE0);
1715             lij3_SSE1          = _mm_mul_pd(lij2_SSE1, lij_SSE1);
1716
1717             diff2_SSE0         = _mm_sub_pd(uij2_SSE0, lij2_SSE0);
1718             diff2_SSE1         = _mm_sub_pd(uij2_SSE1, lij2_SSE1);
1719             lij_inv_SSE0       = gmx_mm_invsqrt_pd(lij2_SSE0);
1720             lij_inv_SSE1       = gmx_mm_invsqrt_pd(lij2_SSE1);
1721             sk2_aj_SSE         = _mm_mul_pd(sk_aj_SSE, sk_aj_SSE);
1722             sk2_rinv_SSE0      = _mm_mul_pd(sk2_aj_SSE, rinv_SSE0);
1723             sk2_rinv_SSE1      = _mm_mul_pd(sk2_aj_SSE, rinv_SSE1);
1724             prod_SSE0          = _mm_mul_pd(onefourth_SSE, sk2_rinv_SSE0);
1725             prod_SSE1          = _mm_mul_pd(onefourth_SSE, sk2_rinv_SSE1);
1726
1727             logterm_SSE0       = gmx_mm_log_pd(_mm_mul_pd(uij_SSE0, lij_inv_SSE0));
1728             logterm_SSE1       = gmx_mm_log_pd(_mm_mul_pd(uij_SSE1, lij_inv_SSE1));
1729
1730             t1_SSE0            = _mm_sub_pd(lij_SSE0, uij_SSE0);
1731             t1_SSE1            = _mm_sub_pd(lij_SSE1, uij_SSE1);
1732             t2_SSE0            = _mm_mul_pd(diff2_SSE0,
1733                                             _mm_sub_pd(_mm_mul_pd(onefourth_SSE, dr_SSE0),
1734                                                        prod_SSE0));
1735             t2_SSE1            = _mm_mul_pd(diff2_SSE1,
1736                                             _mm_sub_pd(_mm_mul_pd(onefourth_SSE, dr_SSE1),
1737                                                        prod_SSE1));
1738
1739             t3_SSE0            = _mm_mul_pd(half_SSE, _mm_mul_pd(rinv_SSE0, logterm_SSE0));
1740             t3_SSE1            = _mm_mul_pd(half_SSE, _mm_mul_pd(rinv_SSE1, logterm_SSE1));
1741             t1_SSE0            = _mm_add_pd(t1_SSE0, _mm_add_pd(t2_SSE0, t3_SSE0));
1742             t1_SSE1            = _mm_add_pd(t1_SSE1, _mm_add_pd(t2_SSE1, t3_SSE1));
1743             t4_SSE0            = _mm_mul_pd(two_SSE, _mm_sub_pd(rai_inv_SSE0, lij_SSE0));
1744             t4_SSE1            = _mm_mul_pd(two_SSE, _mm_sub_pd(rai_inv_SSE1, lij_SSE1));
1745             t4_SSE0            = _mm_and_pd(t4_SSE0, obc_mask3_SSE0);
1746             t4_SSE1            = _mm_and_pd(t4_SSE1, obc_mask3_SSE1);
1747             t1_SSE0            = _mm_mul_pd(half_SSE, _mm_add_pd(t1_SSE0, t4_SSE0));
1748             t1_SSE1            = _mm_mul_pd(half_SSE, _mm_add_pd(t1_SSE1, t4_SSE1));
1749
1750             sum_ai_SSE0        = _mm_add_pd(sum_ai_SSE0, _mm_and_pd(t1_SSE0, obc_mask1_SSE0));
1751             sum_ai_SSE1        = _mm_add_pd(sum_ai_SSE1, _mm_and_pd(t1_SSE1, obc_mask1_SSE1));
1752
1753             t1_SSE0            = _mm_add_pd(_mm_mul_pd(half_SSE, lij2_SSE0),
1754                                             _mm_mul_pd(prod_SSE0, lij3_SSE0));
1755             t1_SSE1            = _mm_add_pd(_mm_mul_pd(half_SSE, lij2_SSE1),
1756                                             _mm_mul_pd(prod_SSE1, lij3_SSE1));
1757
1758             t1_SSE0            = _mm_sub_pd(t1_SSE0,
1759                                             _mm_mul_pd(onefourth_SSE,
1760                                                        _mm_add_pd(_mm_mul_pd(lij_SSE0, rinv_SSE0),
1761                                                                   _mm_mul_pd(lij3_SSE0, dr_SSE0))));
1762             t1_SSE1            = _mm_sub_pd(t1_SSE1,
1763                                             _mm_mul_pd(onefourth_SSE,
1764                                                        _mm_add_pd(_mm_mul_pd(lij_SSE1, rinv_SSE1),
1765                                                                   _mm_mul_pd(lij3_SSE1, dr_SSE1))));
1766
1767             t2_SSE0            = _mm_mul_pd(onefourth_SSE,
1768                                             _mm_add_pd(_mm_mul_pd(uij_SSE0, rinv_SSE0),
1769                                                        _mm_mul_pd(uij3_SSE0, dr_SSE0)));
1770             t2_SSE1            = _mm_mul_pd(onefourth_SSE,
1771                                             _mm_add_pd(_mm_mul_pd(uij_SSE1, rinv_SSE1),
1772                                                        _mm_mul_pd(uij3_SSE1, dr_SSE1)));
1773             t2_SSE0            = _mm_sub_pd(t2_SSE0,
1774                                             _mm_add_pd(_mm_mul_pd(half_SSE, uij2_SSE0),
1775                                                        _mm_mul_pd(prod_SSE0, uij3_SSE0)));
1776             t2_SSE1            = _mm_sub_pd(t2_SSE1,
1777                                             _mm_add_pd(_mm_mul_pd(half_SSE, uij2_SSE1),
1778                                                        _mm_mul_pd(prod_SSE1, uij3_SSE1)));
1779             t3_SSE0            = _mm_mul_pd(_mm_mul_pd(onefourth_SSE, logterm_SSE0),
1780                                             _mm_mul_pd(rinv_SSE0, rinv_SSE0));
1781             t3_SSE1            = _mm_mul_pd(_mm_mul_pd(onefourth_SSE, logterm_SSE1),
1782                                             _mm_mul_pd(rinv_SSE1, rinv_SSE1));
1783             t3_SSE0            = _mm_sub_pd(t3_SSE0,
1784                                             _mm_mul_pd(_mm_mul_pd(diff2_SSE0, oneeighth_SSE),
1785                                                        _mm_add_pd(one_SSE,
1786                                                                   _mm_mul_pd(sk2_rinv_SSE0, rinv_SSE0))));
1787             t3_SSE1            = _mm_sub_pd(t3_SSE1,
1788                                             _mm_mul_pd(_mm_mul_pd(diff2_SSE1, oneeighth_SSE),
1789                                                        _mm_add_pd(one_SSE,
1790                                                                   _mm_mul_pd(sk2_rinv_SSE1, rinv_SSE1))));
1791
1792             t1_SSE0            = _mm_mul_pd(rinv_SSE0,
1793                                             _mm_add_pd(_mm_mul_pd(dlij_SSE0, t1_SSE0),
1794                                                        _mm_add_pd(t2_SSE0, t3_SSE0)));
1795             t1_SSE1            = _mm_mul_pd(rinv_SSE1,
1796                                             _mm_add_pd(_mm_mul_pd(dlij_SSE1, t1_SSE1),
1797                                                        _mm_add_pd(t2_SSE1, t3_SSE1)));
1798
1799             _mm_store_pd(dadx, _mm_and_pd(t1_SSE0, obc_mask1_SSE0));
1800             dadx += 2;
1801             _mm_store_pd(dadx, _mm_and_pd(t1_SSE1, obc_mask1_SSE1));
1802             dadx += 2;
1803
1804             /* Evaluate influence of atom ai -> aj */
1805             t1_SSE0            = _mm_add_pd(dr_SSE0, sk_ai_SSE0);
1806             t1_SSE1            = _mm_add_pd(dr_SSE1, sk_ai_SSE1);
1807             t2_SSE0            = _mm_sub_pd(dr_SSE0, sk_ai_SSE0);
1808             t2_SSE1            = _mm_sub_pd(dr_SSE1, sk_ai_SSE1);
1809             t3_SSE0            = _mm_sub_pd(sk_ai_SSE0, dr_SSE0);
1810             t3_SSE1            = _mm_sub_pd(sk_ai_SSE1, dr_SSE1);
1811
1812             obc_mask1_SSE0     = _mm_cmplt_pd(raj_SSE, t1_SSE0);
1813             obc_mask1_SSE1     = _mm_cmplt_pd(raj_SSE, t1_SSE1);
1814             obc_mask2_SSE0     = _mm_cmplt_pd(raj_SSE, t2_SSE0);
1815             obc_mask2_SSE1     = _mm_cmplt_pd(raj_SSE, t2_SSE1);
1816             obc_mask3_SSE0     = _mm_cmplt_pd(raj_SSE, t3_SSE0);
1817             obc_mask3_SSE1     = _mm_cmplt_pd(raj_SSE, t3_SSE1);
1818             obc_mask1_SSE0     = _mm_and_pd(obc_mask1_SSE0, imask_SSE0);
1819             obc_mask1_SSE1     = _mm_and_pd(obc_mask1_SSE1, imask_SSE1);
1820
1821             uij_SSE0           = gmx_mm_inv_pd(t1_SSE0);
1822             uij_SSE1           = gmx_mm_inv_pd(t1_SSE1);
1823             lij_SSE0           = _mm_or_pd(   _mm_and_pd(obc_mask2_SSE0, gmx_mm_inv_pd(t2_SSE0)),
1824                                               _mm_andnot_pd(obc_mask2_SSE0, raj_inv_SSE));
1825             lij_SSE1           = _mm_or_pd(   _mm_and_pd(obc_mask2_SSE1, gmx_mm_inv_pd(t2_SSE1)),
1826                                               _mm_andnot_pd(obc_mask2_SSE1, raj_inv_SSE));
1827             dlij_SSE0          = _mm_and_pd(one_SSE, obc_mask2_SSE0);
1828             dlij_SSE1          = _mm_and_pd(one_SSE, obc_mask2_SSE1);
1829
1830             uij2_SSE0          = _mm_mul_pd(uij_SSE0, uij_SSE0);
1831             uij2_SSE1          = _mm_mul_pd(uij_SSE1, uij_SSE1);
1832             uij3_SSE0          = _mm_mul_pd(uij2_SSE0, uij_SSE0);
1833             uij3_SSE1          = _mm_mul_pd(uij2_SSE1, uij_SSE1);
1834             lij2_SSE0          = _mm_mul_pd(lij_SSE0, lij_SSE0);
1835             lij2_SSE1          = _mm_mul_pd(lij_SSE1, lij_SSE1);
1836             lij3_SSE0          = _mm_mul_pd(lij2_SSE0, lij_SSE0);
1837             lij3_SSE1          = _mm_mul_pd(lij2_SSE1, lij_SSE1);
1838
1839             diff2_SSE0         = _mm_sub_pd(uij2_SSE0, lij2_SSE0);
1840             diff2_SSE1         = _mm_sub_pd(uij2_SSE1, lij2_SSE1);
1841             lij_inv_SSE0       = gmx_mm_invsqrt_pd(lij2_SSE0);
1842             lij_inv_SSE1       = gmx_mm_invsqrt_pd(lij2_SSE1);
1843             sk2_rinv_SSE0      = _mm_mul_pd(sk2_ai_SSE0, rinv_SSE0);
1844             sk2_rinv_SSE1      = _mm_mul_pd(sk2_ai_SSE1, rinv_SSE1);
1845             prod_SSE0          = _mm_mul_pd(onefourth_SSE, sk2_rinv_SSE0);
1846             prod_SSE1          = _mm_mul_pd(onefourth_SSE, sk2_rinv_SSE1);
1847
1848             logterm_SSE0       = gmx_mm_log_pd(_mm_mul_pd(uij_SSE0, lij_inv_SSE0));
1849             logterm_SSE1       = gmx_mm_log_pd(_mm_mul_pd(uij_SSE1, lij_inv_SSE1));
1850             t1_SSE0            = _mm_sub_pd(lij_SSE0, uij_SSE0);
1851             t1_SSE1            = _mm_sub_pd(lij_SSE1, uij_SSE1);
1852             t2_SSE0            = _mm_mul_pd(diff2_SSE0,
1853                                             _mm_sub_pd(_mm_mul_pd(onefourth_SSE, dr_SSE0),
1854                                                        prod_SSE0));
1855             t2_SSE1            = _mm_mul_pd(diff2_SSE1,
1856                                             _mm_sub_pd(_mm_mul_pd(onefourth_SSE, dr_SSE1),
1857                                                        prod_SSE1));
1858             t3_SSE0            = _mm_mul_pd(half_SSE, _mm_mul_pd(rinv_SSE0, logterm_SSE0));
1859             t3_SSE1            = _mm_mul_pd(half_SSE, _mm_mul_pd(rinv_SSE1, logterm_SSE1));
1860             t1_SSE0            = _mm_add_pd(t1_SSE0, _mm_add_pd(t2_SSE0, t3_SSE0));
1861             t1_SSE1            = _mm_add_pd(t1_SSE1, _mm_add_pd(t2_SSE1, t3_SSE1));
1862             t4_SSE0            = _mm_mul_pd(two_SSE, _mm_sub_pd(raj_inv_SSE, lij_SSE0));
1863             t4_SSE1            = _mm_mul_pd(two_SSE, _mm_sub_pd(raj_inv_SSE, lij_SSE1));
1864             t4_SSE0            = _mm_and_pd(t4_SSE0, obc_mask3_SSE0);
1865             t4_SSE1            = _mm_and_pd(t4_SSE1, obc_mask3_SSE1);
1866             t1_SSE0            = _mm_mul_pd(half_SSE, _mm_add_pd(t1_SSE0, t4_SSE0));
1867             t1_SSE1            = _mm_mul_pd(half_SSE, _mm_add_pd(t1_SSE1, t4_SSE1));
1868
1869             _mm_store_pd(work+j, _mm_add_pd(_mm_load_pd(work+j),
1870                                             _mm_add_pd(_mm_and_pd(t1_SSE0, obc_mask1_SSE0),
1871                                                        _mm_and_pd(t1_SSE1, obc_mask1_SSE1))));
1872
1873             t1_SSE0            = _mm_add_pd(_mm_mul_pd(half_SSE, lij2_SSE0),
1874                                             _mm_mul_pd(prod_SSE0, lij3_SSE0));
1875             t1_SSE1            = _mm_add_pd(_mm_mul_pd(half_SSE, lij2_SSE1),
1876                                             _mm_mul_pd(prod_SSE1, lij3_SSE1));
1877             t1_SSE0            = _mm_sub_pd(t1_SSE0,
1878                                             _mm_mul_pd(onefourth_SSE,
1879                                                        _mm_add_pd(_mm_mul_pd(lij_SSE0, rinv_SSE0),
1880                                                                   _mm_mul_pd(lij3_SSE0, dr_SSE0))));
1881             t1_SSE1            = _mm_sub_pd(t1_SSE1,
1882                                             _mm_mul_pd(onefourth_SSE,
1883                                                        _mm_add_pd(_mm_mul_pd(lij_SSE1, rinv_SSE1),
1884                                                                   _mm_mul_pd(lij3_SSE1, dr_SSE1))));
1885             t2_SSE0            = _mm_mul_pd(onefourth_SSE,
1886                                             _mm_add_pd(_mm_mul_pd(uij_SSE0, rinv_SSE0),
1887                                                        _mm_mul_pd(uij3_SSE0, dr_SSE0)));
1888             t2_SSE1            = _mm_mul_pd(onefourth_SSE,
1889                                             _mm_add_pd(_mm_mul_pd(uij_SSE1, rinv_SSE1),
1890                                                        _mm_mul_pd(uij3_SSE1, dr_SSE1)));
1891             t2_SSE0            = _mm_sub_pd(t2_SSE0,
1892                                             _mm_add_pd(_mm_mul_pd(half_SSE, uij2_SSE0),
1893                                                        _mm_mul_pd(prod_SSE0, uij3_SSE0)));
1894             t2_SSE1            = _mm_sub_pd(t2_SSE1,
1895                                             _mm_add_pd(_mm_mul_pd(half_SSE, uij2_SSE1),
1896                                                        _mm_mul_pd(prod_SSE1, uij3_SSE1)));
1897
1898             t3_SSE0            = _mm_mul_pd(_mm_mul_pd(onefourth_SSE, logterm_SSE0),
1899                                             _mm_mul_pd(rinv_SSE0, rinv_SSE0));
1900             t3_SSE1            = _mm_mul_pd(_mm_mul_pd(onefourth_SSE, logterm_SSE1),
1901                                             _mm_mul_pd(rinv_SSE1, rinv_SSE1));
1902
1903             t3_SSE0            = _mm_sub_pd(t3_SSE0,
1904                                             _mm_mul_pd(_mm_mul_pd(diff2_SSE0, oneeighth_SSE),
1905                                                        _mm_add_pd(one_SSE,
1906                                                                   _mm_mul_pd(sk2_rinv_SSE0, rinv_SSE0))));
1907             t3_SSE1            = _mm_sub_pd(t3_SSE1,
1908                                             _mm_mul_pd(_mm_mul_pd(diff2_SSE1, oneeighth_SSE),
1909                                                        _mm_add_pd(one_SSE,
1910                                                                   _mm_mul_pd(sk2_rinv_SSE1, rinv_SSE1))));
1911
1912             t1_SSE0            = _mm_mul_pd(rinv_SSE0,
1913                                             _mm_add_pd(_mm_mul_pd(dlij_SSE0, t1_SSE0),
1914                                                        _mm_add_pd(t2_SSE0, t3_SSE0)));
1915             t1_SSE1            = _mm_mul_pd(rinv_SSE1,
1916                                             _mm_add_pd(_mm_mul_pd(dlij_SSE1, t1_SSE1),
1917                                                        _mm_add_pd(t2_SSE1, t3_SSE1)));
1918
1919             _mm_store_pd(dadx, _mm_and_pd(t1_SSE0, obc_mask1_SSE0));
1920             dadx += 2;
1921             _mm_store_pd(dadx, _mm_and_pd(t1_SSE1, obc_mask1_SSE1));
1922             dadx += 2;
1923         }
1924
1925         /* Epilogue part, including exclusion mask */
1926         for (j = nj2; j < nj3; j += UNROLLJ)
1927         {
1928             jmask_SSE0 = _mm_load_pd((double *)emask0);
1929             jmask_SSE1 = _mm_load_pd((double *)emask1);
1930             emask0    += 2*UNROLLJ;
1931             emask1    += 2*UNROLLJ;
1932
1933             /* load j atom coordinates */
1934             jx_SSE            = _mm_load_pd(x_align+j);
1935             jy_SSE            = _mm_load_pd(y_align+j);
1936             jz_SSE            = _mm_load_pd(z_align+j);
1937
1938             /* Calculate distance */
1939             dx_SSE0            = _mm_sub_pd(ix_SSE0, jx_SSE);
1940             dy_SSE0            = _mm_sub_pd(iy_SSE0, jy_SSE);
1941             dz_SSE0            = _mm_sub_pd(iz_SSE0, jz_SSE);
1942             dx_SSE1            = _mm_sub_pd(ix_SSE1, jx_SSE);
1943             dy_SSE1            = _mm_sub_pd(iy_SSE1, jy_SSE);
1944             dz_SSE1            = _mm_sub_pd(iz_SSE1, jz_SSE);
1945
1946             /* rsq = dx*dx+dy*dy+dz*dz */
1947             rsq_SSE0           = gmx_mm_calc_rsq_pd(dx_SSE0, dy_SSE0, dz_SSE0);
1948             rsq_SSE1           = gmx_mm_calc_rsq_pd(dx_SSE1, dy_SSE1, dz_SSE1);
1949
1950             /* Combine masks */
1951             jmask_SSE0         = _mm_and_pd(jmask_SSE0, imask_SSE0);
1952             jmask_SSE1         = _mm_and_pd(jmask_SSE1, imask_SSE1);
1953
1954             /* Calculate 1/r and 1/r2 */
1955             rinv_SSE0          = gmx_mm_invsqrt_pd(rsq_SSE0);
1956             rinv_SSE1          = gmx_mm_invsqrt_pd(rsq_SSE1);
1957
1958             /* Apply mask */
1959             rinv_SSE0          = _mm_and_pd(rinv_SSE0, jmask_SSE0);
1960             rinv_SSE1          = _mm_and_pd(rinv_SSE1, jmask_SSE1);
1961
1962             dr_SSE0            = _mm_mul_pd(rsq_SSE0, rinv_SSE0);
1963             dr_SSE1            = _mm_mul_pd(rsq_SSE1, rinv_SSE1);
1964
1965             sk_aj_SSE          = _mm_load_pd(obc_param+j);
1966             raj_SSE            = _mm_load_pd(gb_radius+j);
1967
1968             raj_inv_SSE        = gmx_mm_inv_pd(raj_SSE);
1969
1970             /* Evaluate influence of atom aj -> ai */
1971             t1_SSE0            = _mm_add_pd(dr_SSE0, sk_aj_SSE);
1972             t1_SSE1            = _mm_add_pd(dr_SSE1, sk_aj_SSE);
1973             t2_SSE0            = _mm_sub_pd(dr_SSE0, sk_aj_SSE);
1974             t2_SSE1            = _mm_sub_pd(dr_SSE1, sk_aj_SSE);
1975             t3_SSE0            = _mm_sub_pd(sk_aj_SSE, dr_SSE0);
1976             t3_SSE1            = _mm_sub_pd(sk_aj_SSE, dr_SSE1);
1977
1978             obc_mask1_SSE0     = _mm_cmplt_pd(rai_SSE0, t1_SSE0);
1979             obc_mask1_SSE1     = _mm_cmplt_pd(rai_SSE1, t1_SSE1);
1980             obc_mask2_SSE0     = _mm_cmplt_pd(rai_SSE0, t2_SSE0);
1981             obc_mask2_SSE1     = _mm_cmplt_pd(rai_SSE1, t2_SSE1);
1982             obc_mask3_SSE0     = _mm_cmplt_pd(rai_SSE0, t3_SSE0);
1983             obc_mask3_SSE1     = _mm_cmplt_pd(rai_SSE1, t3_SSE1);
1984             obc_mask1_SSE0     = _mm_and_pd(obc_mask1_SSE0, jmask_SSE0);
1985             obc_mask1_SSE1     = _mm_and_pd(obc_mask1_SSE1, jmask_SSE1);
1986
1987             uij_SSE0           = gmx_mm_inv_pd(t1_SSE0);
1988             uij_SSE1           = gmx_mm_inv_pd(t1_SSE1);
1989             lij_SSE0           = _mm_or_pd(   _mm_and_pd(obc_mask2_SSE0, gmx_mm_inv_pd(t2_SSE0)),
1990                                               _mm_andnot_pd(obc_mask2_SSE0, rai_inv_SSE0));
1991             lij_SSE1           = _mm_or_pd(   _mm_and_pd(obc_mask2_SSE1, gmx_mm_inv_pd(t2_SSE1)),
1992                                               _mm_andnot_pd(obc_mask2_SSE1, rai_inv_SSE1));
1993
1994             dlij_SSE0          = _mm_and_pd(one_SSE, obc_mask2_SSE0);
1995             dlij_SSE1          = _mm_and_pd(one_SSE, obc_mask2_SSE1);
1996
1997             uij2_SSE0          = _mm_mul_pd(uij_SSE0, uij_SSE0);
1998             uij2_SSE1          = _mm_mul_pd(uij_SSE1, uij_SSE1);
1999             uij3_SSE0          = _mm_mul_pd(uij2_SSE0, uij_SSE0);
2000             uij3_SSE1          = _mm_mul_pd(uij2_SSE1, uij_SSE1);
2001             lij2_SSE0          = _mm_mul_pd(lij_SSE0, lij_SSE0);
2002             lij2_SSE1          = _mm_mul_pd(lij_SSE1, lij_SSE1);
2003             lij3_SSE0          = _mm_mul_pd(lij2_SSE0, lij_SSE0);
2004             lij3_SSE1          = _mm_mul_pd(lij2_SSE1, lij_SSE1);
2005
2006             diff2_SSE0         = _mm_sub_pd(uij2_SSE0, lij2_SSE0);
2007             diff2_SSE1         = _mm_sub_pd(uij2_SSE1, lij2_SSE1);
2008             lij_inv_SSE0       = gmx_mm_invsqrt_pd(lij2_SSE0);
2009             lij_inv_SSE1       = gmx_mm_invsqrt_pd(lij2_SSE1);
2010             sk2_aj_SSE         = _mm_mul_pd(sk_aj_SSE, sk_aj_SSE);
2011             sk2_rinv_SSE0      = _mm_mul_pd(sk2_aj_SSE, rinv_SSE0);
2012             sk2_rinv_SSE1      = _mm_mul_pd(sk2_aj_SSE, rinv_SSE1);
2013             prod_SSE0          = _mm_mul_pd(onefourth_SSE, sk2_rinv_SSE0);
2014             prod_SSE1          = _mm_mul_pd(onefourth_SSE, sk2_rinv_SSE1);
2015
2016             logterm_SSE0       = gmx_mm_log_pd(_mm_mul_pd(uij_SSE0, lij_inv_SSE0));
2017             logterm_SSE1       = gmx_mm_log_pd(_mm_mul_pd(uij_SSE1, lij_inv_SSE1));
2018
2019             t1_SSE0            = _mm_sub_pd(lij_SSE0, uij_SSE0);
2020             t1_SSE1            = _mm_sub_pd(lij_SSE1, uij_SSE1);
2021             t2_SSE0            = _mm_mul_pd(diff2_SSE0,
2022                                             _mm_sub_pd(_mm_mul_pd(onefourth_SSE, dr_SSE0),
2023                                                        prod_SSE0));
2024             t2_SSE1            = _mm_mul_pd(diff2_SSE1,
2025                                             _mm_sub_pd(_mm_mul_pd(onefourth_SSE, dr_SSE1),
2026                                                        prod_SSE1));
2027
2028             t3_SSE0            = _mm_mul_pd(half_SSE, _mm_mul_pd(rinv_SSE0, logterm_SSE0));
2029             t3_SSE1            = _mm_mul_pd(half_SSE, _mm_mul_pd(rinv_SSE1, logterm_SSE1));
2030             t1_SSE0            = _mm_add_pd(t1_SSE0, _mm_add_pd(t2_SSE0, t3_SSE0));
2031             t1_SSE1            = _mm_add_pd(t1_SSE1, _mm_add_pd(t2_SSE1, t3_SSE1));
2032             t4_SSE0            = _mm_mul_pd(two_SSE, _mm_sub_pd(rai_inv_SSE0, lij_SSE0));
2033             t4_SSE1            = _mm_mul_pd(two_SSE, _mm_sub_pd(rai_inv_SSE1, lij_SSE1));
2034             t4_SSE0            = _mm_and_pd(t4_SSE0, obc_mask3_SSE0);
2035             t4_SSE1            = _mm_and_pd(t4_SSE1, obc_mask3_SSE1);
2036             t1_SSE0            = _mm_mul_pd(half_SSE, _mm_add_pd(t1_SSE0, t4_SSE0));
2037             t1_SSE1            = _mm_mul_pd(half_SSE, _mm_add_pd(t1_SSE1, t4_SSE1));
2038
2039             sum_ai_SSE0        = _mm_add_pd(sum_ai_SSE0, _mm_and_pd(t1_SSE0, obc_mask1_SSE0));
2040             sum_ai_SSE1        = _mm_add_pd(sum_ai_SSE1, _mm_and_pd(t1_SSE1, obc_mask1_SSE1));
2041
2042             t1_SSE0            = _mm_add_pd(_mm_mul_pd(half_SSE, lij2_SSE0),
2043                                             _mm_mul_pd(prod_SSE0, lij3_SSE0));
2044             t1_SSE1            = _mm_add_pd(_mm_mul_pd(half_SSE, lij2_SSE1),
2045                                             _mm_mul_pd(prod_SSE1, lij3_SSE1));
2046             t1_SSE0            = _mm_sub_pd(t1_SSE0,
2047                                             _mm_mul_pd(onefourth_SSE,
2048                                                        _mm_add_pd(_mm_mul_pd(lij_SSE0, rinv_SSE0),
2049                                                                   _mm_mul_pd(lij3_SSE0, dr_SSE0))));
2050             t1_SSE1            = _mm_sub_pd(t1_SSE1,
2051                                             _mm_mul_pd(onefourth_SSE,
2052                                                        _mm_add_pd(_mm_mul_pd(lij_SSE1, rinv_SSE1),
2053                                                                   _mm_mul_pd(lij3_SSE1, dr_SSE1))));
2054
2055             t2_SSE0            = _mm_mul_pd(onefourth_SSE,
2056                                             _mm_add_pd(_mm_mul_pd(uij_SSE0, rinv_SSE0),
2057                                                        _mm_mul_pd(uij3_SSE0, dr_SSE0)));
2058             t2_SSE1            = _mm_mul_pd(onefourth_SSE,
2059                                             _mm_add_pd(_mm_mul_pd(uij_SSE1, rinv_SSE1),
2060                                                        _mm_mul_pd(uij3_SSE1, dr_SSE1)));
2061             t2_SSE0            = _mm_sub_pd(t2_SSE0,
2062                                             _mm_add_pd(_mm_mul_pd(half_SSE, uij2_SSE0),
2063                                                        _mm_mul_pd(prod_SSE0, uij3_SSE0)));
2064             t2_SSE1            = _mm_sub_pd(t2_SSE1,
2065                                             _mm_add_pd(_mm_mul_pd(half_SSE, uij2_SSE1),
2066                                                        _mm_mul_pd(prod_SSE1, uij3_SSE1)));
2067             t3_SSE0            = _mm_mul_pd(_mm_mul_pd(onefourth_SSE, logterm_SSE0),
2068                                             _mm_mul_pd(rinv_SSE0, rinv_SSE0));
2069             t3_SSE1            = _mm_mul_pd(_mm_mul_pd(onefourth_SSE, logterm_SSE1),
2070                                             _mm_mul_pd(rinv_SSE1, rinv_SSE1));
2071             t3_SSE0            = _mm_sub_pd(t3_SSE0,
2072                                             _mm_mul_pd(_mm_mul_pd(diff2_SSE0, oneeighth_SSE),
2073                                                        _mm_add_pd(one_SSE,
2074                                                                   _mm_mul_pd(sk2_rinv_SSE0, rinv_SSE0))));
2075             t3_SSE1            = _mm_sub_pd(t3_SSE1,
2076                                             _mm_mul_pd(_mm_mul_pd(diff2_SSE1, oneeighth_SSE),
2077                                                        _mm_add_pd(one_SSE,
2078                                                                   _mm_mul_pd(sk2_rinv_SSE1, rinv_SSE1))));
2079
2080             t1_SSE0            = _mm_mul_pd(rinv_SSE0,
2081                                             _mm_add_pd(_mm_mul_pd(dlij_SSE0, t1_SSE0),
2082                                                        _mm_add_pd(t2_SSE0, t3_SSE0)));
2083             t1_SSE1            = _mm_mul_pd(rinv_SSE1,
2084                                             _mm_add_pd(_mm_mul_pd(dlij_SSE1, t1_SSE1),
2085                                                        _mm_add_pd(t2_SSE1, t3_SSE1)));
2086
2087             _mm_store_pd(dadx, _mm_and_pd(t1_SSE0, obc_mask1_SSE0));
2088             dadx += 2;
2089             _mm_store_pd(dadx, _mm_and_pd(t1_SSE1, obc_mask1_SSE1));
2090             dadx += 2;
2091
2092             /* Evaluate influence of atom ai -> aj */
2093             t1_SSE0            = _mm_add_pd(dr_SSE0, sk_ai_SSE0);
2094             t1_SSE1            = _mm_add_pd(dr_SSE1, sk_ai_SSE1);
2095             t2_SSE0            = _mm_sub_pd(dr_SSE0, sk_ai_SSE0);
2096             t2_SSE1            = _mm_sub_pd(dr_SSE1, sk_ai_SSE1);
2097             t3_SSE0            = _mm_sub_pd(sk_ai_SSE0, dr_SSE0);
2098             t3_SSE1            = _mm_sub_pd(sk_ai_SSE1, dr_SSE1);
2099
2100             obc_mask1_SSE0     = _mm_cmplt_pd(raj_SSE, t1_SSE0);
2101             obc_mask1_SSE1     = _mm_cmplt_pd(raj_SSE, t1_SSE1);
2102             obc_mask2_SSE0     = _mm_cmplt_pd(raj_SSE, t2_SSE0);
2103             obc_mask2_SSE1     = _mm_cmplt_pd(raj_SSE, t2_SSE1);
2104             obc_mask3_SSE0     = _mm_cmplt_pd(raj_SSE, t3_SSE0);
2105             obc_mask3_SSE1     = _mm_cmplt_pd(raj_SSE, t3_SSE1);
2106             obc_mask1_SSE0     = _mm_and_pd(obc_mask1_SSE0, jmask_SSE0);
2107             obc_mask1_SSE1     = _mm_and_pd(obc_mask1_SSE1, jmask_SSE1);
2108
2109             uij_SSE0           = gmx_mm_inv_pd(t1_SSE0);
2110             uij_SSE1           = gmx_mm_inv_pd(t1_SSE1);
2111             lij_SSE0           = _mm_or_pd(   _mm_and_pd(obc_mask2_SSE0, gmx_mm_inv_pd(t2_SSE0)),
2112                                               _mm_andnot_pd(obc_mask2_SSE0, raj_inv_SSE));
2113             lij_SSE1           = _mm_or_pd(   _mm_and_pd(obc_mask2_SSE1, gmx_mm_inv_pd(t2_SSE1)),
2114                                               _mm_andnot_pd(obc_mask2_SSE1, raj_inv_SSE));
2115
2116             dlij_SSE0          = _mm_and_pd(one_SSE, obc_mask2_SSE0);
2117             dlij_SSE1          = _mm_and_pd(one_SSE, obc_mask2_SSE1);
2118
2119             uij2_SSE0          = _mm_mul_pd(uij_SSE0, uij_SSE0);
2120             uij2_SSE1          = _mm_mul_pd(uij_SSE1, uij_SSE1);
2121             uij3_SSE0          = _mm_mul_pd(uij2_SSE0, uij_SSE0);
2122             uij3_SSE1          = _mm_mul_pd(uij2_SSE1, uij_SSE1);
2123             lij2_SSE0          = _mm_mul_pd(lij_SSE0, lij_SSE0);
2124             lij2_SSE1          = _mm_mul_pd(lij_SSE1, lij_SSE1);
2125             lij3_SSE0          = _mm_mul_pd(lij2_SSE0, lij_SSE0);
2126             lij3_SSE1          = _mm_mul_pd(lij2_SSE1, lij_SSE1);
2127
2128             diff2_SSE0         = _mm_sub_pd(uij2_SSE0, lij2_SSE0);
2129             diff2_SSE1         = _mm_sub_pd(uij2_SSE1, lij2_SSE1);
2130             lij_inv_SSE0       = gmx_mm_invsqrt_pd(lij2_SSE0);
2131             lij_inv_SSE1       = gmx_mm_invsqrt_pd(lij2_SSE1);
2132             sk2_rinv_SSE0      = _mm_mul_pd(sk2_ai_SSE0, rinv_SSE0);
2133             sk2_rinv_SSE1      = _mm_mul_pd(sk2_ai_SSE1, rinv_SSE1);
2134             prod_SSE0          = _mm_mul_pd(onefourth_SSE, sk2_rinv_SSE0);
2135             prod_SSE1          = _mm_mul_pd(onefourth_SSE, sk2_rinv_SSE1);
2136
2137             logterm_SSE0       = gmx_mm_log_pd(_mm_mul_pd(uij_SSE0, lij_inv_SSE0));
2138             logterm_SSE1       = gmx_mm_log_pd(_mm_mul_pd(uij_SSE1, lij_inv_SSE1));
2139             t1_SSE0            = _mm_sub_pd(lij_SSE0, uij_SSE0);
2140             t1_SSE1            = _mm_sub_pd(lij_SSE1, uij_SSE1);
2141             t2_SSE0            = _mm_mul_pd(diff2_SSE0,
2142                                             _mm_sub_pd(_mm_mul_pd(onefourth_SSE, dr_SSE0),
2143                                                        prod_SSE0));
2144             t2_SSE1            = _mm_mul_pd(diff2_SSE1,
2145                                             _mm_sub_pd(_mm_mul_pd(onefourth_SSE, dr_SSE1),
2146                                                        prod_SSE1));
2147             t3_SSE0            = _mm_mul_pd(half_SSE, _mm_mul_pd(rinv_SSE0, logterm_SSE0));
2148             t3_SSE1            = _mm_mul_pd(half_SSE, _mm_mul_pd(rinv_SSE1, logterm_SSE1));
2149             t1_SSE0            = _mm_add_pd(t1_SSE0, _mm_add_pd(t2_SSE0, t3_SSE0));
2150             t1_SSE1            = _mm_add_pd(t1_SSE1, _mm_add_pd(t2_SSE1, t3_SSE1));
2151             t4_SSE0            = _mm_mul_pd(two_SSE, _mm_sub_pd(raj_inv_SSE, lij_SSE0));
2152             t4_SSE1            = _mm_mul_pd(two_SSE, _mm_sub_pd(raj_inv_SSE, lij_SSE1));
2153             t4_SSE0            = _mm_and_pd(t4_SSE0, obc_mask3_SSE0);
2154             t4_SSE1            = _mm_and_pd(t4_SSE1, obc_mask3_SSE1);
2155             t1_SSE0            = _mm_mul_pd(half_SSE, _mm_add_pd(t1_SSE0, t4_SSE0));
2156             t1_SSE1            = _mm_mul_pd(half_SSE, _mm_add_pd(t1_SSE1, t4_SSE1));
2157
2158             _mm_store_pd(work+j, _mm_add_pd(_mm_load_pd(work+j),
2159                                             _mm_add_pd(_mm_and_pd(t1_SSE0, obc_mask1_SSE0),
2160                                                        _mm_and_pd(t1_SSE1, obc_mask1_SSE1))));
2161
2162             t1_SSE0            = _mm_add_pd(_mm_mul_pd(half_SSE, lij2_SSE0),
2163                                             _mm_mul_pd(prod_SSE0, lij3_SSE0));
2164             t1_SSE1            = _mm_add_pd(_mm_mul_pd(half_SSE, lij2_SSE1),
2165                                             _mm_mul_pd(prod_SSE1, lij3_SSE1));
2166
2167             t1_SSE0            = _mm_sub_pd(t1_SSE0,
2168                                             _mm_mul_pd(onefourth_SSE,
2169                                                        _mm_add_pd(_mm_mul_pd(lij_SSE0, rinv_SSE0),
2170                                                                   _mm_mul_pd(lij3_SSE0, dr_SSE0))));
2171             t1_SSE1            = _mm_sub_pd(t1_SSE1,
2172                                             _mm_mul_pd(onefourth_SSE,
2173                                                        _mm_add_pd(_mm_mul_pd(lij_SSE1, rinv_SSE1),
2174                                                                   _mm_mul_pd(lij3_SSE1, dr_SSE1))));
2175             t2_SSE0            = _mm_mul_pd(onefourth_SSE,
2176                                             _mm_add_pd(_mm_mul_pd(uij_SSE0, rinv_SSE0),
2177                                                        _mm_mul_pd(uij3_SSE0, dr_SSE0)));
2178             t2_SSE1            = _mm_mul_pd(onefourth_SSE,
2179                                             _mm_add_pd(_mm_mul_pd(uij_SSE1, rinv_SSE1),
2180                                                        _mm_mul_pd(uij3_SSE1, dr_SSE1)));
2181             t2_SSE0            = _mm_sub_pd(t2_SSE0,
2182                                             _mm_add_pd(_mm_mul_pd(half_SSE, uij2_SSE0),
2183                                                        _mm_mul_pd(prod_SSE0, uij3_SSE0)));
2184             t2_SSE1            = _mm_sub_pd(t2_SSE1,
2185                                             _mm_add_pd(_mm_mul_pd(half_SSE, uij2_SSE1),
2186                                                        _mm_mul_pd(prod_SSE1, uij3_SSE1)));
2187
2188             t3_SSE0            = _mm_mul_pd(_mm_mul_pd(onefourth_SSE, logterm_SSE0),
2189                                             _mm_mul_pd(rinv_SSE0, rinv_SSE0));
2190             t3_SSE1            = _mm_mul_pd(_mm_mul_pd(onefourth_SSE, logterm_SSE1),
2191                                             _mm_mul_pd(rinv_SSE1, rinv_SSE1));
2192
2193             t3_SSE0            = _mm_sub_pd(t3_SSE0,
2194                                             _mm_mul_pd(_mm_mul_pd(diff2_SSE0, oneeighth_SSE),
2195                                                        _mm_add_pd(one_SSE,
2196                                                                   _mm_mul_pd(sk2_rinv_SSE0, rinv_SSE0))));
2197             t3_SSE1            = _mm_sub_pd(t3_SSE1,
2198                                             _mm_mul_pd(_mm_mul_pd(diff2_SSE1, oneeighth_SSE),
2199                                                        _mm_add_pd(one_SSE,
2200                                                                   _mm_mul_pd(sk2_rinv_SSE1, rinv_SSE1))));
2201
2202             t1_SSE0            = _mm_mul_pd(rinv_SSE0,
2203                                             _mm_add_pd(_mm_mul_pd(dlij_SSE0, t1_SSE0),
2204                                                        _mm_add_pd(t2_SSE0, t3_SSE0)));
2205             t1_SSE1            = _mm_mul_pd(rinv_SSE1,
2206                                             _mm_add_pd(_mm_mul_pd(dlij_SSE1, t1_SSE1),
2207                                                        _mm_add_pd(t2_SSE1, t3_SSE1)));
2208
2209             _mm_store_pd(dadx, _mm_and_pd(t1_SSE0, obc_mask1_SSE0));
2210             dadx += 2;
2211             _mm_store_pd(dadx, _mm_and_pd(t1_SSE1, obc_mask1_SSE1));
2212             dadx += 2;
2213         }
2214         GMX_MM_TRANSPOSE2_PD(sum_ai_SSE0, sum_ai_SSE1);
2215         sum_ai_SSE0 = _mm_add_pd(sum_ai_SSE0, sum_ai_SSE1);
2216         _mm_store_pd(work+i, _mm_add_pd(sum_ai_SSE0, _mm_load_pd(work+i)));
2217     }
2218
2219
2220     for (i = 0; i < natoms/2+1; i++)
2221     {
2222         work[i] += work[natoms+i];
2223     }
2224
2225     /* Parallel summations */
2226
2227     if (PARTDECOMP(cr))
2228     {
2229         gmx_sum(natoms, work, cr);
2230     }
2231
2232     if (gb_algorithm == egbHCT)
2233     {
2234         /* HCT */
2235         for (i = 0; i < natoms; i++)
2236         {
2237             if (born->use[i] != 0)
2238             {
2239                 rai     = top->atomtypes.gb_radius[mdatoms->typeA[i]]-born->gb_doffset;
2240                 sum_ai  = 1.0/rai - work[i];
2241                 min_rad = rai + born->gb_doffset;
2242                 rad     = 1.0/sum_ai;
2243
2244                 born->bRad[i]   = rad > min_rad ? rad : min_rad;
2245                 fr->invsqrta[i] = gmx_invsqrt(born->bRad[i]);
2246             }
2247         }
2248
2249     }
2250     else
2251     {
2252         /* OBC */
2253
2254         /* Calculate the radii */
2255         for (i = 0; i < natoms; i++)
2256         {
2257
2258             if (born->use[i] != 0)
2259             {
2260                 rai        = top->atomtypes.gb_radius[mdatoms->typeA[i]];
2261                 rai_inv2   = 1.0/rai;
2262                 rai        = rai-born->gb_doffset;
2263                 rai_inv    = 1.0/rai;
2264                 sum_ai     = rai * work[i];
2265                 sum_ai2    = sum_ai  * sum_ai;
2266                 sum_ai3    = sum_ai2 * sum_ai;
2267
2268                 tsum          = tanh(born->obc_alpha*sum_ai-born->obc_beta*sum_ai2+born->obc_gamma*sum_ai3);
2269                 born->bRad[i] = rai_inv - tsum*rai_inv2;
2270                 born->bRad[i] = 1.0 / born->bRad[i];
2271
2272                 fr->invsqrta[i] = gmx_invsqrt(born->bRad[i]);
2273
2274                 tchain         = rai * (born->obc_alpha-2*born->obc_beta*sum_ai+3*born->obc_gamma*sum_ai2);
2275                 born->drobc[i] = (1.0-tsum*tsum)*tchain*rai_inv2;
2276             }
2277         }
2278     }
2279
2280     return 0;
2281 }
2282
2283
2284
2285
2286
2287
2288
2289
2290 int
2291 genborn_allvsall_calc_chainrule_sse2_double(t_forcerec   *           fr,
2292                                             t_mdatoms   *            mdatoms,
2293                                             gmx_genborn_t   *        born,
2294                                             double *                 x,
2295                                             double *                 f,
2296                                             int                      gb_algorithm,
2297                                             void   *                 paadata)
2298 {
2299     gmx_allvsallgb2_data_t *aadata;
2300     int                     natoms;
2301     int                     ni0, ni1;
2302     int                     nj0, nj1, nj2, nj3;
2303     int                     i, j, k, n;
2304     int                     idx;
2305     int              *      mask;
2306     int              *      pmask0;
2307     int              *      emask0;
2308     int              *      jindex;
2309
2310     double                  ix, iy, iz;
2311     double                  fix, fiy, fiz;
2312     double                  jx, jy, jz;
2313     double                  dx, dy, dz;
2314     double                  tx, ty, tz;
2315     double                  rbai, rbaj, fgb, fgb_ai, rbi;
2316     double            *     rb;
2317     double            *     dadx;
2318     double            *     x_align;
2319     double            *     y_align;
2320     double            *     z_align;
2321     double            *     fx_align;
2322     double            *     fy_align;
2323     double            *     fz_align;
2324     double                  tmpsum[2];
2325
2326     __m128d                 jmask_SSE0, jmask_SSE1;
2327     __m128d                 ix_SSE0, iy_SSE0, iz_SSE0;
2328     __m128d                 ix_SSE1, iy_SSE1, iz_SSE1;
2329     __m128d                 fix_SSE0, fiy_SSE0, fiz_SSE0;
2330     __m128d                 fix_SSE1, fiy_SSE1, fiz_SSE1;
2331     __m128d                 rbai_SSE0, rbai_SSE1;
2332     __m128d                 imask_SSE0, imask_SSE1;
2333     __m128d                 jx_SSE, jy_SSE, jz_SSE, rbaj_SSE;
2334     __m128d                 dx_SSE0, dy_SSE0, dz_SSE0;
2335     __m128d                 dx_SSE1, dy_SSE1, dz_SSE1;
2336     __m128d                 fgb_SSE0, fgb_ai_SSE0;
2337     __m128d                 fgb_SSE1, fgb_ai_SSE1;
2338     __m128d                 tx_SSE0, ty_SSE0, tz_SSE0;
2339     __m128d                 tx_SSE1, ty_SSE1, tz_SSE1;
2340     __m128d                 t1, t2, tmpSSE;
2341
2342     natoms              = mdatoms->nr;
2343     ni0                 = (mdatoms->start/SIMD_WIDTH)*SIMD_WIDTH;
2344     ni1                 = mdatoms->start+mdatoms->homenr;
2345
2346     aadata = (gmx_allvsallgb2_data_t *)paadata;
2347
2348     x_align  = aadata->x_align;
2349     y_align  = aadata->y_align;
2350     z_align  = aadata->z_align;
2351     fx_align = aadata->fx_align;
2352     fy_align = aadata->fy_align;
2353     fz_align = aadata->fz_align;
2354
2355     jindex    = aadata->jindex_gb;
2356     dadx      = fr->dadx;
2357
2358     n  = 0;
2359     rb = aadata->work;
2360
2361     /* Loop to get the proper form for the Born radius term */
2362     if (gb_algorithm == egbSTILL)
2363     {
2364         for (i = 0; i < natoms; i++)
2365         {
2366             rbi   = born->bRad[i];
2367             rb[i] = (2 * rbi * rbi * fr->dvda[i])/ONE_4PI_EPS0;
2368         }
2369     }
2370     else if (gb_algorithm == egbHCT)
2371     {
2372         for (i = 0; i < natoms; i++)
2373         {
2374             rbi   = born->bRad[i];
2375             rb[i] = rbi * rbi * fr->dvda[i];
2376         }
2377     }
2378     else if (gb_algorithm == egbOBC)
2379     {
2380         for (idx = 0; idx < natoms; idx++)
2381         {
2382             rbi     = born->bRad[idx];
2383             rb[idx] = rbi * rbi * born->drobc[idx] * fr->dvda[idx];
2384         }
2385     }
2386
2387     for (i = 0; i < 2*natoms; i++)
2388     {
2389         fx_align[i]       = 0;
2390         fy_align[i]       = 0;
2391         fz_align[i]       = 0;
2392     }
2393
2394
2395     for (i = 0; i < natoms; i++)
2396     {
2397         rb[i+natoms] = rb[i];
2398     }
2399
2400     for (i = ni0; i < ni1; i += UNROLLI)
2401     {
2402         /* We assume shifts are NOT used for all-vs-all interactions */
2403
2404         /* Load i atom data */
2405         ix_SSE0          = _mm_load1_pd(x_align+i);
2406         iy_SSE0          = _mm_load1_pd(y_align+i);
2407         iz_SSE0          = _mm_load1_pd(z_align+i);
2408         ix_SSE1          = _mm_load1_pd(x_align+i+1);
2409         iy_SSE1          = _mm_load1_pd(y_align+i+1);
2410         iz_SSE1          = _mm_load1_pd(z_align+i+1);
2411
2412         fix_SSE0         = _mm_setzero_pd();
2413         fiy_SSE0         = _mm_setzero_pd();
2414         fiz_SSE0         = _mm_setzero_pd();
2415         fix_SSE1         = _mm_setzero_pd();
2416         fiy_SSE1         = _mm_setzero_pd();
2417         fiz_SSE1         = _mm_setzero_pd();
2418
2419         rbai_SSE0        = _mm_load1_pd(rb+i);
2420         rbai_SSE1        = _mm_load1_pd(rb+i+1);
2421
2422         /* Load limits for loop over neighbors */
2423         nj0              = jindex[4*i];
2424         nj3              = jindex[4*i+3];
2425
2426         /* No masks necessary, since the stored chain rule derivatives will be zero in those cases! */
2427         for (j = nj0; j < nj3; j += UNROLLJ)
2428         {
2429             /* load j atom coordinates */
2430             jx_SSE           = _mm_load_pd(x_align+j);
2431             jy_SSE           = _mm_load_pd(y_align+j);
2432             jz_SSE           = _mm_load_pd(z_align+j);
2433
2434             /* Calculate distance */
2435             dx_SSE0          = _mm_sub_pd(ix_SSE0, jx_SSE);
2436             dy_SSE0          = _mm_sub_pd(iy_SSE0, jy_SSE);
2437             dz_SSE0          = _mm_sub_pd(iz_SSE0, jz_SSE);
2438             dx_SSE1          = _mm_sub_pd(ix_SSE1, jx_SSE);
2439             dy_SSE1          = _mm_sub_pd(iy_SSE1, jy_SSE);
2440             dz_SSE1          = _mm_sub_pd(iz_SSE1, jz_SSE);
2441
2442             rbaj_SSE         = _mm_load_pd(rb+j);
2443
2444             fgb_SSE0         = _mm_mul_pd(rbai_SSE0, _mm_load_pd(dadx));
2445             dadx            += 2;
2446             fgb_SSE1         = _mm_mul_pd(rbai_SSE1, _mm_load_pd(dadx));
2447             dadx            += 2;
2448
2449             fgb_ai_SSE0      = _mm_mul_pd(rbaj_SSE, _mm_load_pd(dadx));
2450             dadx            += 2;
2451             fgb_ai_SSE1      = _mm_mul_pd(rbaj_SSE, _mm_load_pd(dadx));
2452             dadx            += 2;
2453
2454             /* Total force between ai and aj is the sum of ai->aj and aj->ai */
2455             fgb_SSE0         = _mm_add_pd(fgb_SSE0, fgb_ai_SSE0);
2456             fgb_SSE1         = _mm_add_pd(fgb_SSE1, fgb_ai_SSE1);
2457
2458             /* Calculate temporary vectorial force */
2459             tx_SSE0            = _mm_mul_pd(fgb_SSE0, dx_SSE0);
2460             ty_SSE0            = _mm_mul_pd(fgb_SSE0, dy_SSE0);
2461             tz_SSE0            = _mm_mul_pd(fgb_SSE0, dz_SSE0);
2462             tx_SSE1            = _mm_mul_pd(fgb_SSE1, dx_SSE1);
2463             ty_SSE1            = _mm_mul_pd(fgb_SSE1, dy_SSE1);
2464             tz_SSE1            = _mm_mul_pd(fgb_SSE1, dz_SSE1);
2465
2466             /* Increment i atom force */
2467             fix_SSE0          = _mm_add_pd(fix_SSE0, tx_SSE0);
2468             fiy_SSE0          = _mm_add_pd(fiy_SSE0, ty_SSE0);
2469             fiz_SSE0          = _mm_add_pd(fiz_SSE0, tz_SSE0);
2470             fix_SSE1          = _mm_add_pd(fix_SSE1, tx_SSE1);
2471             fiy_SSE1          = _mm_add_pd(fiy_SSE1, ty_SSE1);
2472             fiz_SSE1          = _mm_add_pd(fiz_SSE1, tz_SSE1);
2473
2474             /* Decrement j atom force */
2475             _mm_store_pd(fx_align+j,
2476                          _mm_sub_pd( _mm_load_pd(fx_align+j), _mm_add_pd(tx_SSE0, tx_SSE1) ));
2477             _mm_store_pd(fy_align+j,
2478                          _mm_sub_pd( _mm_load_pd(fy_align+j), _mm_add_pd(ty_SSE0, ty_SSE1) ));
2479             _mm_store_pd(fz_align+j,
2480                          _mm_sub_pd( _mm_load_pd(fz_align+j), _mm_add_pd(tz_SSE0, tz_SSE1) ));
2481         }
2482
2483         /* Add i forces to mem */
2484         GMX_MM_TRANSPOSE2_PD(fix_SSE0, fix_SSE1);
2485         fix_SSE0 = _mm_add_pd(fix_SSE0, fix_SSE1);
2486         _mm_store_pd(fx_align+i, _mm_add_pd(fix_SSE0, _mm_load_pd(fx_align+i)));
2487
2488         GMX_MM_TRANSPOSE2_PD(fiy_SSE0, fiy_SSE1);
2489         fiy_SSE0 = _mm_add_pd(fiy_SSE0, fiy_SSE1);
2490         _mm_store_pd(fy_align+i, _mm_add_pd(fiy_SSE0, _mm_load_pd(fy_align+i)));
2491
2492         GMX_MM_TRANSPOSE2_PD(fiz_SSE0, fiz_SSE1);
2493         fiz_SSE0 = _mm_add_pd(fiz_SSE0, fiz_SSE1);
2494         _mm_store_pd(fz_align+i, _mm_add_pd(fiz_SSE0, _mm_load_pd(fz_align+i)));
2495     }
2496
2497     for (i = 0; i < natoms; i++)
2498     {
2499         f[3*i]       += fx_align[i] + fx_align[natoms+i];
2500         f[3*i+1]     += fy_align[i] + fy_align[natoms+i];
2501         f[3*i+2]     += fz_align[i] + fz_align[natoms+i];
2502     }
2503
2504     return 0;
2505 }
2506
2507 #else
2508 /* dummy variable when not using SSE */
2509 int genborn_allvsall_sse2_double_dummy;
2510
2511
2512 #endif